add complex log, exp, pow

This commit is contained in:
Fredrik Johansson 2012-11-09 17:45:25 +01:00
parent 1dd795be42
commit 45ad953f69
6 changed files with 295 additions and 1 deletions

View file

@ -123,6 +123,15 @@ Precision and comparisons
Returns nonzero iff zero is contained in *x*.
Complex parts
-------------------------------------------------------------------------------
.. function:: void fmpcb_arg(fmprb_t r, const fmpcb_t z, long prec)
Sets *r* to a real interval containing the complex argument of *z*. We
define the complex argument have a discontinuity on `(-\infty,0]`, with
the special value `\operatorname{arg}(0) = 0`, and
`\operatorname{arg}(x+0i) = \pi` for `x < 0`.
Arithmetic
-------------------------------------------------------------------------------
@ -175,9 +184,27 @@ Arithmetic
Sets *z* to the quotient of *x* and *y*.
Elementary functions
-------------------------------------------------------------------------------
.. function:: void fmpcb_log(fmpcb_t y, const fmpcb_t z, long prec)
Sets *y* to the principal branch of the natural logarithm of *z*,
computed as
`\log(a+bi) = \frac{1}{2} \log(a^2 + b^2) + i \operatorname{arg}(a+bi)`.
.. function:: void fmpcb_exp(fmpcb_t y, const fmpcb_t z, long prec)
Sets *y* to the exponential function of *z*, computed as
`\exp(a+bi) = \exp(a) \left( \cos(b) + \sin(b) i \right)`.
.. function:: void fmpcb_pow_fmpz(fmpcb_t y, const fmpcb_t b, const fmpz_t e, long prec)
.. function:: void fmpcb_pow_ui(fmpcb_t y, const fmpcb_t b, ulong e, long prec)
Sets *y* to *b* raised to the power *e*.
Sets *y* to *b* raised to the power *e*, computed using binary exponentiation.
.. function:: void fmpcb_pow(fmpcb_t r, const fmpcb_t x, const fmpcb_t y, long prec)
Sets *r* to *x* raised to the power *y*, computed as `x^y = \exp(y \log x)`.

View file

@ -198,6 +198,8 @@ fmpcb_get_abs_ubound_fmpr(fmpr_t u, const fmpcb_t z, long prec)
fmpr_clear(v);
}
void fmpcb_arg(fmprb_t r, const fmpcb_t z, long prec);
static __inline__ void
fmpcb_add(fmpcb_t z, const fmpcb_t x, const fmpcb_t y, long prec)
{
@ -314,6 +316,12 @@ void fmpcb_pow_fmpz(fmpcb_t y, const fmpcb_t b, const fmpz_t e, long prec);
void fmpcb_pow_ui(fmpcb_t y, const fmpcb_t b, ulong e, long prec);
void fmpcb_log(fmpcb_t r, const fmpcb_t z, long prec);
void fmpcb_exp(fmpcb_t r, const fmpcb_t z, long prec);
void fmpcb_pow(fmpcb_t r, const fmpcb_t x, const fmpcb_t y, long prec);
static __inline__ void
_fmpcb_vec_zero(fmpcb_struct * A, long n)
{

137
fmpcb/arg.c Normal file
View file

@ -0,0 +1,137 @@
/*=============================================================================
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2012 Fredrik Johansson
******************************************************************************/
#include "fmpcb.h"
void
fmpcb_arg(fmprb_t r, const fmpcb_t z, long prec)
{
#define a fmpcb_realref(z)
#define b fmpcb_imagref(z)
#define am fmprb_midref(a)
#define ar fmprb_radref(a)
#define bm fmprb_midref(b)
#define br fmprb_radref(b)
/* a real number */
if (fmprb_is_zero(b))
{
/* exactly zero */
if (fmprb_is_zero(a))
{
/* define arg(0) = 0 by convention */
fmprb_zero(r);
}
/* interval contains only nonnegative numbers */
else if (fmpr_sgn(am) > 0 && fmpr_cmpabs(am, ar) >= 0)
{
fmprb_zero(r);
}
/* interval contains only positive numbers */
else if (fmpr_sgn(am) < 0 && fmpr_cmpabs(am, ar) > 0)
{
fmprb_const_pi(r, prec);
}
else
{
/* both positive and negative -- argument will be in [0, pi] */
fmprb_t t;
fmprb_init(t);
fmprb_const_pi(t, prec);
fmprb_mul_2exp_si(t, t, -1);
fmprb_zero(r);
fmprb_add_error(r, t);
fmprb_clear(t);
}
}
/* an imaginary number */
else if (fmprb_is_zero(a))
{
/* interval contains only positive numbers */
if (fmpr_sgn(bm) > 0 && fmpr_cmpabs(bm, br) > 0)
{
fmprb_const_pi(r, prec);
fmprb_mul_2exp_si(r, r, -1);
}
/* interval contains only negative numbers */
else if (fmpr_sgn(bm) < 0 && fmpr_cmpabs(bm, br) > 0)
{
fmprb_const_pi(r, prec);
fmprb_neg(r, r);
fmprb_mul_2exp_si(r, r, -1);
}
else
{
/* both positive and negative -- argument will be in 0 +/- pi/2 */
fmprb_t t;
fmprb_const_pi(t, FMPRB_RAD_PREC);
fmprb_mul_2exp_si(t, t, -1);
fmprb_zero(r);
fmprb_add_error(r, t);
fmprb_clear(t);
}
}
/* strictly in the right half-plane -- atan(b/a) */
else if (fmpr_sgn(am) > 0 && fmpr_cmpabs(am, ar) > 0)
{
fmprb_div(r, b, a, prec);
fmprb_atan(r, r, prec);
}
/* strictly in the upper half-plane -- pi/2 - atan(a/b) */
else if (fmpr_sgn(bm) > 0 && fmpr_cmpabs(bm, br) > 0)
{
fmprb_t t;
fmprb_init(t);
fmprb_div(r, a, b, prec);
fmprb_atan(r, r, prec);
fmprb_const_pi(t, prec);
fmprb_mul_2exp_si(t, t, -1);
fmprb_sub(r, t, r, prec);
fmprb_clear(t);
}
/* strictly in the lower half-plane -- -pi/2 - atan(a/b) */
else if (fmpr_sgn(bm) < 0 && fmpr_cmpabs(bm, br) > 0)
{
fmprb_t t;
fmprb_init(t);
fmprb_div(r, a, b, prec);
fmprb_atan(r, r, prec);
fmprb_const_pi(t, prec);
fmprb_mul_2exp_si(t, t, -1);
fmprb_add(r, t, r, prec);
fmprb_neg(r, r);
fmprb_clear(t);
}
/* overlaps the nonpositive half-axis -- [-pi, pi] */
else
{
fmprb_t t;
fmprb_const_pi(t, FMPRB_RAD_PREC);
fmprb_zero(r);
fmprb_add_error(r, t);
fmprb_clear(t);
}
}

50
fmpcb/exp.c Normal file
View file

@ -0,0 +1,50 @@
/*=============================================================================
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2012 Fredrik Johansson
******************************************************************************/
#include "fmpcb.h"
void
fmpcb_exp(fmpcb_t r, const fmpcb_t z, long prec)
{
#define a fmpcb_realref(z)
#define b fmpcb_imagref(z)
fmprb_t t, u, v;
fmprb_init(t);
fmprb_init(u);
fmprb_init(v);
fmprb_exp(t, a, prec);
fmprb_sin_cos(u, v, b, prec);
fmprb_mul(fmpcb_realref(r), t, v, prec);
fmprb_mul(fmpcb_imagref(r), t, u, prec);
fmprb_clear(t);
fmprb_clear(u);
fmprb_clear(v);
}

61
fmpcb/log.c Normal file
View file

@ -0,0 +1,61 @@
/*=============================================================================
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2012 Fredrik Johansson
******************************************************************************/
#include "fmpcb.h"
void
fmpcb_log(fmpcb_t r, const fmpcb_t z, long prec)
{
#define a fmpcb_realref(z)
#define b fmpcb_imagref(z)
fmprb_t t, u;
fmprb_init(t);
fmprb_init(u);
fmprb_mul(t, a, a, prec);
fmprb_mul(u, b, b, prec);
fmprb_add(t, t, u, prec);
if (fmprb_contains_zero(t) || fmpr_sgn(fmprb_midref(t)) < 0)
{
fmpr_zero(fmprb_midref(t));
fmpr_pos_inf(fmprb_radref(t));
}
else
{
fmprb_log(t, t, prec);
}
fmpcb_arg(u, z, prec);
fmprb_mul_2exp_si(fmpcb_realref(r), t, -1);
fmprb_set(fmpcb_imagref(r), u);
fmprb_clear(t);
fmprb_clear(u);
}

View file

@ -78,3 +78,14 @@ fmpcb_pow_ui(fmpcb_t y, const fmpcb_t b, ulong e, long prec)
fmpz_clear(f);
}
void
fmpcb_pow(fmpcb_t r, const fmpcb_t x, const fmpcb_t y, long prec)
{
fmpcb_t t;
fmpcb_init(t);
fmpcb_log(t, x, prec);
fmpcb_mul(t, t, y, prec);
fmpcb_exp(r, t, prec);
fmpcb_clear(t);
}