/* mpc_pow -- Raise a complex number to the power of another complex number.
-Copyright (C) 2009, 2010, 2011, 2012 INRIA
+Copyright (C) 2009 Paul Zimmermann
-This file is part of GNU MPC.
+This file is part of the MPC Library.
-GNU MPC is free software; you can redistribute it and/or modify it under
-the terms of the GNU Lesser General Public License as published by the
-Free Software Foundation; either version 3 of the License, or (at your
+The MPC Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published by
+the Free Software Foundation; either version 2.1 of the License, or (at your
option) any later version.
-GNU MPC 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 Lesser General Public License for
-more details.
+The MPC Library 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 Lesser General Public
+License for more details.
You should have received a copy of the GNU Lesser General Public License
-along with this program. If not, see http://www.gnu.org/licenses/ .
-*/
+along with the MPC Library; see the file COPYING.LIB. If not, write to
+the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
+MA 02111-1307, USA. */
#include <stdio.h> /* for MPC_ASSERT */
#include "mpc-impl.h"
return 0; /* not a square */
}
-/* fix the sign of Re(z) or Im(z) in case it is zero,
- and Re(x) is zero.
- sign_eps is 0 if Re(x) = +0, 1 if Re(x) = -0
- sign_a is the sign bit of Im(x).
- Assume y is an integer (does nothing otherwise).
-*/
-static void
-fix_sign (mpc_ptr z, int sign_eps, int sign_a, mpfr_srcptr y)
-{
- int ymod4 = -1;
- mpfr_exp_t ey;
- mpz_t my;
- unsigned long int t;
-
- mpz_init (my);
-
- ey = mpfr_get_z_exp (my, y);
- /* normalize so that my is odd */
- t = mpz_scan1 (my, 0);
- ey += (mpfr_exp_t) t;
- mpz_tdiv_q_2exp (my, my, t);
- /* y = my*2^ey */
-
- /* compute y mod 4 (in case y is an integer) */
- if (ey >= 2)
- ymod4 = 0;
- else if (ey == 1)
- ymod4 = mpz_tstbit (my, 0) * 2; /* correct if my < 0 */
- else if (ey == 0)
- {
- ymod4 = mpz_tstbit (my, 1) * 2 + mpz_tstbit (my, 0);
- if (mpz_cmp_ui (my , 0) < 0)
- ymod4 = 4 - ymod4;
- }
- else /* y is not an integer */
- goto end;
-
- if (mpfr_zero_p (mpc_realref(z)))
- {
- /* we assume y is always integer in that case (FIXME: prove it):
- (eps+I*a)^y = +0 + I*a^y for y = 1 mod 4 and sign_eps = 0
- (eps+I*a)^y = -0 - I*a^y for y = 3 mod 4 and sign_eps = 0 */
- MPC_ASSERT (ymod4 == 1 || ymod4 == 3);
- if ((ymod4 == 3 && sign_eps == 0) ||
- (ymod4 == 1 && sign_eps == 1))
- mpfr_neg (mpc_realref(z), mpc_realref(z), GMP_RNDZ);
- }
- else if (mpfr_zero_p (mpc_imagref(z)))
- {
- /* we assume y is always integer in that case (FIXME: prove it):
- (eps+I*a)^y = a^y - 0*I for y = 0 mod 4 and sign_a = sign_eps
- (eps+I*a)^y = -a^y +0*I for y = 2 mod 4 and sign_a = sign_eps */
- MPC_ASSERT (ymod4 == 0 || ymod4 == 2);
- if ((ymod4 == 0 && sign_a == sign_eps) ||
- (ymod4 == 2 && sign_a != sign_eps))
- mpfr_neg (mpc_imagref(z), mpc_imagref(z), GMP_RNDZ);
- }
-
- end:
- mpz_clear (my);
-}
-
/* If x^y is exactly representable (with maybe a larger precision than z),
round it in z and return the (mpc) inexact flag in [0, 10].
should be called again with a larger value of maxprec).
Assume one of Re(x) or Im(x) is non-zero, and y is non-zero (y is real).
-
- Warning: z and x might be the same variable, same for Re(z) or Im(z) and y.
-
- In case -1 or -2 is returned, z is not modified.
*/
static int
mpc_pow_exact (mpc_ptr z, mpc_srcptr x, mpfr_srcptr y, mpc_rnd_t rnd,
- mpfr_prec_t maxprec)
+ mp_prec_t maxprec)
{
- mpfr_exp_t ec, ed, ey;
+ mp_exp_t ec, ed, ey, emin, emax;
mpz_t my, a, b, c, d, u;
unsigned long int t;
int ret = -2;
- int sign_rex = mpfr_signbit (mpc_realref(x));
- int sign_imx = mpfr_signbit (mpc_imagref(x));
- int x_imag = mpfr_zero_p (mpc_realref(x));
- int z_is_y = 0;
- mpfr_t copy_of_y;
-
- if (mpc_realref (z) == y || mpc_imagref (z) == y)
- {
- z_is_y = 1;
- mpfr_init2 (copy_of_y, mpfr_get_prec (y));
- mpfr_set (copy_of_y, y, GMP_RNDN);
- }
mpz_init (my);
mpz_init (a);
ey = mpfr_get_z_exp (my, y);
/* normalize so that my is odd */
t = mpz_scan1 (my, 0);
- ey += (mpfr_exp_t) t;
+ ey += t;
mpz_tdiv_q_2exp (my, my, t);
- /* y = my*2^ey with my odd */
- if (x_imag)
+ if (mpfr_zero_p (MPC_RE(x)))
{
mpz_set_ui (c, 0);
ec = 0;
}
else
- ec = mpfr_get_z_exp (c, mpc_realref(x));
- if (mpfr_zero_p (mpc_imagref(x)))
+ ec = mpfr_get_z_exp (c, MPC_RE(x));
+ if (mpfr_zero_p (MPC_IM(x)))
{
mpz_set_ui (d, 0);
ed = ec;
}
else
{
- ed = mpfr_get_z_exp (d, mpc_imagref(x));
- if (x_imag)
+ ed = mpfr_get_z_exp (d, MPC_IM(x));
+ if (mpfr_zero_p (MPC_RE(x)))
ec = ed;
}
/* x = c*2^ec + I * d*2^ed */
/* equalize the exponents of x */
if (ec < ed)
{
- mpz_mul_2exp (d, d, (unsigned long int) (ed - ec));
- if ((mpfr_prec_t) mpz_sizeinbase (d, 2) > maxprec)
+ mpz_mul_2exp (d, d, ed - ec);
+ if (mpz_sizeinbase (d, 2) > maxprec)
goto end;
+ ed = ec;
}
else if (ed < ec)
{
- mpz_mul_2exp (c, c, (unsigned long int) (ec - ed));
- if ((mpfr_prec_t) mpz_sizeinbase (c, 2) > maxprec)
+ mpz_mul_2exp (c, c, ec - ed);
+ if (mpz_sizeinbase (c, 2) > maxprec)
goto end;
ec = ed;
}
{
t = mpz_scan1 (d, 0);
mpz_tdiv_q_2exp (d, d, t);
- ec += (mpfr_exp_t) t;
+ ec += t;
}
else if (mpz_cmp_ui (d, 0) == 0)
{
t = mpz_scan1 (c, 0);
mpz_tdiv_q_2exp (c, c, t);
- ec += (mpfr_exp_t) t;
+ ec += t;
}
else /* neither c nor d is zero */
{
t = v;
mpz_tdiv_q_2exp (c, c, t);
mpz_tdiv_q_2exp (d, d, t);
- ec += (mpfr_exp_t) t;
+ ec += t;
}
/* now either one of c, d is odd */
}
/* replace (c,d) by (c/(c^2+d^2), -d/(c^2+d^2)) */
mpz_neg (d, d);
- ec = -ec - (mpfr_exp_t) t;
+ ec = -ec - t;
mpz_neg (my, my);
}
/* invariant: (a + I*b) * 2^ed = ((c + I*d) * 2^ec)^trunc(my/2^t) */
while (t-- > 0)
{
- unsigned long int v, w;
+ unsigned long v, w;
/* square a + I*b */
mpz_mul (u, a, b);
mpz_mul (a, a, a);
{
w = mpz_scan1 (b, 0);
mpz_tdiv_q_2exp (b, b, w);
- ed += (mpfr_exp_t) w;
+ ed += w;
}
else if (mpz_cmp_ui (b, 0) == 0)
{
w = mpz_scan1 (a, 0);
mpz_tdiv_q_2exp (a, a, w);
- ed += (mpfr_exp_t) w;
+ ed += w;
}
else
{
w = v;
mpz_tdiv_q_2exp (a, a, w);
mpz_tdiv_q_2exp (b, b, w);
- ed += (mpfr_exp_t) w;
+ ed += w;
}
- if ( (mpfr_prec_t) mpz_sizeinbase (a, 2) > maxprec
- || (mpfr_prec_t) mpz_sizeinbase (b, 2) > maxprec)
+ if (mpz_sizeinbase (a, 2) > maxprec || mpz_sizeinbase (b, 2) > maxprec)
goto end;
}
/* now a+I*b = (c+I*d)^my */
sa = (sa <= sb) ? sa : sb;
mpz_tdiv_q_2exp (a, a, sa);
mpz_tdiv_q_2exp (b, b, sa);
- ed += (mpfr_exp_t) sa;
+ ed += sa;
- if ( (mpfr_prec_t) mpz_sizeinbase (a, 2) > maxprec
- || (mpfr_prec_t) mpz_sizeinbase (b, 2) > maxprec)
+ if (mpz_sizeinbase (a, 2) > maxprec || mpz_sizeinbase (b, 2) > maxprec)
goto end;
}
- ret = mpfr_set_z (mpc_realref(z), a, MPC_RND_RE(rnd));
- ret = MPC_INEX(ret, mpfr_set_z (mpc_imagref(z), b, MPC_RND_IM(rnd)));
- mpfr_mul_2si (mpc_realref(z), mpc_realref(z), ed, MPC_RND_RE(rnd));
- mpfr_mul_2si (mpc_imagref(z), mpc_imagref(z), ed, MPC_RND_IM(rnd));
+ /* save emin, emax */
+ emin = mpfr_get_emin ();
+ emax = mpfr_get_emax ();
+ mpfr_set_emin (mpfr_get_emin_min ());
+ mpfr_set_emax (mpfr_get_emax_max ());
+ ret = mpfr_set_z (MPC_RE(z), a, MPC_RND_RE(rnd));
+ ret = MPC_INEX(ret, mpfr_set_z (MPC_IM(z), b, MPC_RND_IM(rnd)));
+ mpfr_mul_2si (MPC_RE(z), MPC_RE(z), ed, MPC_RND_RE(rnd));
+ mpfr_mul_2si (MPC_IM(z), MPC_IM(z), ed, MPC_RND_IM(rnd));
+ /* restore emin, emax */
+ mpfr_set_emin (emin);
+ mpfr_set_emax (emax);
end:
mpz_clear (my);
mpz_clear (d);
mpz_clear (u);
- if (ret >= 0 && x_imag)
- fix_sign (z, sign_rex, sign_imx, (z_is_y) ? copy_of_y : y);
-
- if (z_is_y)
- mpfr_clear (copy_of_y);
-
return ret;
}
*/
#define MPFR_LIMB_HIGHBIT ((mp_limb_t) 1 << (BITS_PER_MP_LIMB - 1))
static int
-is_odd (mpfr_srcptr y, mpfr_exp_t k)
+is_odd (mpfr_srcptr y, mp_exp_t k)
{
- mpfr_exp_t expo;
- mpfr_prec_t prec;
+ mp_exp_t expo;
+ mp_prec_t prec;
mp_size_t yn;
mp_limb_t *yp;
int
mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd)
{
- int ret = -2, loop, x_real, x_imag, y_real, z_real = 0, z_imag = 0;
+ int ret = -2, loop, x_real, y_real, z_real = 0, z_imag = 0;
mpc_t t, u;
- mpfr_prec_t p, pr, pi, maxprec;
- int saved_underflow, saved_overflow;
-
- /* save the underflow or overflow flags from MPFR */
- saved_underflow = mpfr_underflow_p ();
- saved_overflow = mpfr_overflow_p ();
+ mp_prec_t p, q, pr, pi, maxprec;
+ long Q;
- x_real = mpfr_zero_p (mpc_imagref(x));
- y_real = mpfr_zero_p (mpc_imagref(y));
+ x_real = mpfr_zero_p (MPC_IM(x));
+ y_real = mpfr_zero_p (MPC_IM(y));
- if (y_real && mpfr_zero_p (mpc_realref(y))) /* case y zero */
+ if (y_real && mpfr_zero_p (MPC_RE(y))) /* case y zero */
{
- if (x_real && mpfr_zero_p (mpc_realref(x)))
+ if (x_real && mpfr_zero_p (MPC_RE(x))) /* 0^0 = NaN +i*NaN */
{
- /* we define 0^0 to be (1, +0) since the real part is
- coherent with MPFR where 0^0 gives 1, and the sign of the
- imaginary part cannot be determined */
- mpc_set_ui_ui (z, 1, 0, MPC_RNDNN);
+ mpfr_set_nan (MPC_RE(z));
+ mpfr_set_nan (MPC_IM(z));
return 0;
}
else /* x^0 = 1 +/- i*0 even for x=NaN see algorithms.tex for the
if (cx1 == 0 && inex != 0)
cx1 = -inex;
- sign_zi = (cx1 < 0 && mpfr_signbit (mpc_imagref (y)) == 0)
+ sign_zi = (cx1 < 0 && mpfr_signbit (MPC_IM (y)) == 0)
|| (cx1 == 0
- && mpfr_signbit (mpc_imagref (x)) != mpfr_signbit (mpc_realref (y)))
- || (cx1 > 0 && mpfr_signbit (mpc_imagref (y)));
+ && mpfr_signbit (MPC_IM (x)) != mpfr_signbit (MPC_RE (y)))
+ || (cx1 > 0 && mpfr_signbit (MPC_IM (y)));
/* warning: mpc_set_ui_ui does not set Im(z) to -0 if Im(rnd)=RNDD */
ret = mpc_set_ui_ui (z, 1, 0, rnd);
if (x_real) /* case x real */
{
- if (mpfr_zero_p (mpc_realref(x))) /* x is zero */
+ if (mpfr_zero_p (MPC_RE(x))) /* x is zero */
{
/* special values: exp(y*log(x)) */
mpc_init2 (u, 2);
}
/* Special case 1^y = 1 */
- if (mpfr_cmp_ui (mpc_realref(x), 1) == 0)
+ if (mpfr_cmp_ui (MPC_RE(x), 1) == 0)
{
int s1, s2;
- s1 = mpfr_signbit (mpc_realref (y));
- s2 = mpfr_signbit (mpc_imagref (x));
+ s1 = mpfr_signbit (MPC_RE (y));
+ s2 = mpfr_signbit (MPC_IM (x));
ret = mpc_set_ui (z, +1, rnd);
/* the sign of the zero imaginary part is known in some cases (see
/* x^y is real when:
(a) x is real and y is integer
(b) x is real non-negative and y is real */
- if (y_real && (mpfr_integer_p (mpc_realref(y)) ||
- mpfr_cmp_ui (mpc_realref(x), 0) >= 0))
+ if (y_real && (mpfr_integer_p (MPC_RE(y)) ||
+ mpfr_cmp_ui (MPC_RE(x), 0) >= 0))
{
int s1, s2;
- s1 = mpfr_signbit (mpc_realref (y));
- s2 = mpfr_signbit (mpc_imagref (x));
+ s1 = mpfr_signbit (MPC_RE (y));
+ s2 = mpfr_signbit (MPC_IM (x));
- ret = mpfr_pow (mpc_realref(z), mpc_realref(x), mpc_realref(y), MPC_RND_RE(rnd));
- ret = MPC_INEX(ret, mpfr_set_ui (mpc_imagref(z), 0, MPC_RND_IM(rnd)));
+ ret = mpfr_pow (MPC_RE(z), MPC_RE(x), MPC_RE(y), MPC_RND_RE(rnd));
+ ret = MPC_INEX(ret, mpfr_set_ui (MPC_IM(z), 0, MPC_RND_IM(rnd)));
/* the sign of the zero imaginary part is known in some cases
(see algorithm.tex). In such cases we have (x +s*0i)^(y+/-0i)
because mpfr_set_ui(z_i, 0, rnd) always sets z_i to +0.
*/
if (MPC_RND_IM(rnd) == GMP_RNDD || s1 != s2)
- mpfr_neg (mpc_imagref(z), mpc_imagref(z), MPC_RND_IM(rnd));
+ mpfr_neg (MPC_IM(z), MPC_IM(z), MPC_RND_IM(rnd));
goto end;
}
/* (-1)^(n+I*t) is real for n integer and t real */
- if (mpfr_cmp_si (mpc_realref(x), -1) == 0 && mpfr_integer_p (mpc_realref(y)))
+ if (mpfr_cmp_si (MPC_RE(x), -1) == 0 && mpfr_integer_p (MPC_RE(y)))
z_real = 1;
/* for x real, x^y is imaginary when:
(a) x is negative and y is half-an-integer
(b) x = -1 and Re(y) is half-an-integer
*/
- if ((mpfr_cmp_ui (mpc_realref(x), 0) < 0) && is_odd (mpc_realref(y), 1)
- && (y_real || mpfr_cmp_si (mpc_realref(x), -1) == 0))
+ if (mpfr_cmp_ui (MPC_RE(x), 0) < 0 && is_odd (MPC_RE(y), 1) &&
+ (y_real || mpfr_cmp_si (MPC_RE(x), -1) == 0))
z_imag = 1;
}
else /* x non real */
I^(n+t*I) and (-I)^(n+t*I) are imaginary for n odd and t real
(s*I)^n is real for n even and imaginary for n odd */
if ((mpc_cmp_si_si (x, 0, 1) == 0 || mpc_cmp_si_si (x, 0, -1) == 0 ||
- (mpfr_cmp_ui (mpc_realref(x), 0) == 0 && y_real)) &&
- mpfr_integer_p (mpc_realref(y)))
+ (mpfr_cmp_ui (MPC_RE(x), 0) == 0 && y_real)) &&
+ mpfr_integer_p (MPC_RE(y)))
{ /* x is I or -I, and Re(y) is an integer */
- if (is_odd (mpc_realref(y), 0))
+ if (is_odd (MPC_RE(y), 0))
z_imag = 1; /* Re(y) odd: z is imaginary */
else
z_real = 1; /* Re(y) even: z is real */
}
else /* (t+/-t*I)^(2n) is imaginary for n odd and real for n even */
- if (mpfr_cmpabs (mpc_realref(x), mpc_imagref(x)) == 0 && y_real &&
- mpfr_integer_p (mpc_realref(y)) && is_odd (mpc_realref(y), 0) == 0)
+ if (mpfr_cmpabs (MPC_RE(x), MPC_IM(x)) == 0 && y_real &&
+ mpfr_integer_p (MPC_RE(y)) && is_odd (MPC_RE(y), 0) == 0)
{
- if (is_odd (mpc_realref(y), -1)) /* y/2 is odd */
+ if (is_odd (MPC_RE(y), -1)) /* y/2 is odd */
z_imag = 1;
else
z_real = 1;
}
- pr = mpfr_get_prec (mpc_realref(z));
- pi = mpfr_get_prec (mpc_imagref(z));
+ /* first bound |Re(y log(x))|, |Im(y log(x)| < 2^q */
+ mpc_init2 (t, 64);
+ mpc_log (t, x, MPC_RNDNN);
+ mpc_mul (t, t, y, MPC_RNDNN);
+
+ /* the default maximum exponent for MPFR is emax=2^30-1, thus if
+ t > log(2^emax) = emax*log(2), then exp(t) will overflow */
+ if (mpfr_cmp_ui_2exp (MPC_RE(t), 372130558, 1) > 0)
+ goto overflow;
+
+ /* the default minimum exponent for MPFR is emin=-2^30+1, thus the
+ smallest representable value is 2^(emin-1), and if
+ t < log(2^(emin-1)) = (emin-1)*log(2), then exp(t) will underflow */
+ if (mpfr_cmp_si_2exp (MPC_RE(t), -372130558, 1) < 0)
+ goto underflow;
+
+ q = mpfr_get_exp (MPC_RE(t)) > 0 ? mpfr_get_exp (MPC_RE(t)) : 0;
+ if (mpfr_get_exp (MPC_IM(t)) > (mp_exp_t) q)
+ q = mpfr_get_exp (MPC_IM(t));
+
+ pr = mpfr_get_prec (MPC_RE(z));
+ pi = mpfr_get_prec (MPC_IM(z));
p = (pr > pi) ? pr : pi;
p += 12; /* experimentally, seems to give less than 10% of failures in
- Ziv's strategy; probably wrong now since q is not computed */
- if (p < 64)
- p = 64;
+ Ziv's strategy */
mpc_init2 (u, p);
- mpc_init2 (t, p);
pr += MPC_RND_RE(rnd) == GMP_RNDN;
pi += MPC_RND_IM(rnd) == GMP_RNDN;
- maxprec = MPC_MAX_PREC (z);
- x_imag = mpfr_zero_p (mpc_realref(x));
+ maxprec = MPFR_PREC(MPC_RE(z));
+ if (MPFR_PREC(MPC_IM(z)) > maxprec)
+ maxprec = MPFR_PREC(MPC_IM(z));
for (loop = 0;; loop++)
{
- int ret_exp;
- mpfr_exp_t dr, di;
- mpfr_prec_t q=0;
- /* to avoid warning message, real initialisation below */
-
- mpc_log (t, x, MPC_RNDNN);
- mpc_mul (t, t, y, MPC_RNDNN);
-
- if (loop == 0) {
- /* compute q such that |Re (y log x)|, |Im (y log x)| < 2^q */
- q = mpfr_get_exp (mpc_realref(t)) > 0 ? mpfr_get_exp (mpc_realref(t)) : 0;
- if (mpfr_get_exp (mpc_imagref(t)) > (mpfr_exp_t) q)
- q = mpfr_get_exp (mpc_imagref(t));
- }
-
- mpfr_clear_overflow ();
- mpfr_clear_underflow ();
- ret_exp = mpc_exp (u, t, MPC_RNDNN);
- if (mpfr_underflow_p () || mpfr_overflow_p ()) {
- /* under- and overflow flags are set by mpc_exp */
- mpc_set (z, u, MPC_RNDNN);
- ret = ret_exp;
- goto exact;
- }
+ mp_exp_t dr, di;
+ if (p + q > 64) /* otherwise we reuse the initial approximation
+ t of y*log(x), avoiding two computations */
+ {
+ mpc_set_prec (t, p + q);
+ mpc_log (t, x, MPC_RNDNN);
+ mpc_mul (t, t, y, MPC_RNDNN);
+ }
+ mpc_exp (u, t, MPC_RNDNN);
/* Since the error bound is global, we have to take into account the
exponent difference between the real and imaginary parts. We assume
either the real or the imaginary part of u is not zero.
*/
- dr = mpfr_zero_p (mpc_realref(u)) ? mpfr_get_exp (mpc_imagref(u))
- : mpfr_get_exp (mpc_realref(u));
- di = mpfr_zero_p (mpc_imagref(u)) ? dr : mpfr_get_exp (mpc_imagref(u));
+ dr = mpfr_zero_p (MPC_RE(u)) ? mpfr_get_exp (MPC_IM(u))
+ : mpfr_get_exp (MPC_RE(u));
+ di = mpfr_zero_p (MPC_IM(u)) ? dr : mpfr_get_exp (MPC_IM(u));
if (dr > di)
{
di = dr - di;
(see algorithms.tex) plus one due to the exponent difference: if
z = a + I*b, where the relative error on z is at most 2^(-p), and
EXP(a) = EXP(b) + k, the relative error on b is at most 2^(k-p) */
- if ((z_imag || (p > q + 3 + dr && mpfr_can_round (mpc_realref(u), p - q - 3 - dr, GMP_RNDN, GMP_RNDZ, pr))) &&
- (z_real || (p > q + 3 + di && mpfr_can_round (mpc_imagref(u), p - q - 3 - di, GMP_RNDN, GMP_RNDZ, pi))))
+ if ((z_imag || mpfr_can_round (MPC_RE(u), p - 3 - dr, GMP_RNDN, GMP_RNDZ, pr)) &&
+ (z_real || mpfr_can_round (MPC_IM(u), p - 3 - di, GMP_RNDN, GMP_RNDZ, pi)))
break;
/* if Re(u) is not known to be zero, assume it is a normal number, i.e.,
neither zero, Inf or NaN, otherwise we might enter an infinite loop */
- MPC_ASSERT (z_imag || mpfr_number_p (mpc_realref(u)));
+ MPC_ASSERT (z_imag || mpfr_number_p (MPC_RE(u)));
/* idem for Im(u) */
- MPC_ASSERT (z_real || mpfr_number_p (mpc_imagref(u)));
+ MPC_ASSERT (z_real || mpfr_number_p (MPC_IM(u)));
if (ret == -2) /* we did not yet call mpc_pow_exact, or it aborted
because intermediate computations had > maxprec bits */
if (y_real)
{
maxprec *= 2;
- ret = mpc_pow_exact (z, x, mpc_realref(y), rnd, maxprec);
+ ret = mpc_pow_exact (z, x, MPC_RE(y), rnd, maxprec);
if (ret != -1 && ret != -2)
goto exact;
}
}
else
p += p / 2;
- mpc_set_prec (t, p);
+ mpc_set_prec (t, p + q);
mpc_set_prec (u, p);
}
*/
mpfr_t n;
int inex, cx1;
- int sign_zi, sign_rex, sign_imx;
+ int sign_zi;
/* cx1 < 0 if |x| < 1
cx1 = 0 if |x| = 1
cx1 > 0 if |x| > 1
*/
-
- sign_rex = mpfr_signbit (mpc_realref (x));
- sign_imx = mpfr_signbit (mpc_imagref (x));
mpfr_init (n);
inex = mpc_norm (n, x, GMP_RNDN);
cx1 = mpfr_cmp_ui (n, 1);
if (cx1 == 0 && inex != 0)
cx1 = -inex;
- sign_zi = (cx1 < 0 && mpfr_signbit (mpc_imagref (y)) == 0)
- || (cx1 == 0 && sign_imx != mpfr_signbit (mpc_realref (y)))
- || (cx1 > 0 && mpfr_signbit (mpc_imagref (y)));
+ sign_zi = (cx1 < 0 && mpfr_signbit (MPC_IM (y)) == 0)
+ || (cx1 == 0
+ && mpfr_signbit (MPC_IM (x)) != mpfr_signbit (MPC_RE (y)))
+ || (cx1 > 0 && mpfr_signbit (MPC_IM (y)));
- /* copy RE(y) to n since if z==y we will destroy Re(y) below */
- mpfr_set_prec (n, mpfr_get_prec (mpc_realref (y)));
- mpfr_set (n, mpc_realref (y), GMP_RNDN);
- ret = mpfr_set (mpc_realref(z), mpc_realref(u), MPC_RND_RE(rnd));
- if (y_real && (x_real || x_imag))
- {
- /* FIXME: with y_real we assume Im(y) is really 0, which is the case
- for example when y comes from pow_fr, but in case Im(y) is +0 or
- -0, we might get different results */
- mpfr_set_ui (mpc_imagref (z), 0, MPC_RND_IM (rnd));
- fix_sign (z, sign_rex, sign_imx, n);
- ret = MPC_INEX(ret, 0); /* imaginary part is exact */
- }
- else
- {
- ret = MPC_INEX (ret, mpfr_set_ui (mpc_imagref (z), 0, MPC_RND_IM (rnd)));
- /* warning: mpfr_set_ui does not set Im(z) to -0 if Im(rnd) = RNDD */
- if (MPC_RND_IM (rnd) == GMP_RNDD || sign_zi)
- mpc_conj (z, z, MPC_RNDNN);
- }
+ ret = mpfr_set (MPC_RE(z), MPC_RE(u), MPC_RND_RE(rnd));
+ /* warning: mpfr_set_ui does not set Im(z) to -0 if Im(rnd) = RNDD */
+ ret = MPC_INEX (ret, mpfr_set_ui (MPC_IM (z), 0, MPC_RND_IM (rnd)));
+
+ if (MPC_RND_IM (rnd) == GMP_RNDD || sign_zi)
+ mpc_conj (z, z, MPC_RNDNN);
mpfr_clear (n);
}
else if (z_imag)
{
- ret = mpfr_set (mpc_imagref(z), mpc_imagref(u), MPC_RND_IM(rnd));
- /* if z is imaginary and y real, then x cannot be real */
- if (y_real && x_imag)
- {
- int sign_rex = mpfr_signbit (mpc_realref (x));
-
- /* If z overlaps with y we set Re(z) before checking Re(y) below,
- but in that case y=0, which was dealt with above. */
- mpfr_set_ui (mpc_realref (z), 0, MPC_RND_RE (rnd));
- /* Note: fix_sign only does something when y is an integer,
- then necessarily y = 1 or 3 (mod 4), and in that case the
- sign of Im(x) is irrelevant. */
- fix_sign (z, sign_rex, 0, mpc_realref (y));
- ret = MPC_INEX(0, ret);
- }
- else
- ret = MPC_INEX(mpfr_set_ui (mpc_realref(z), 0, MPC_RND_RE(rnd)), ret);
+ ret = mpfr_set (MPC_IM(z), MPC_IM(u), MPC_RND_IM(rnd));
+ ret = MPC_INEX(mpfr_set_ui (MPC_RE(z), 0, MPC_RND_RE(rnd)), ret);
}
else
ret = mpc_set (z, u, rnd);
mpc_clear (t);
mpc_clear (u);
- /* restore underflow and overflow flags from MPFR */
- if (saved_underflow)
- mpfr_set_underflow ();
- if (saved_overflow)
- mpfr_set_overflow ();
-
end:
return ret;
+
+ underflow:
+ /* If we have an underflow, we know that |z| is too small to be
+ represented, but depending on arg(z), we should return +/-0 +/- I*0.
+ We assume t is the approximation of y*log(x), thus we want
+ exp(t) = exp(Re(t))+exp(I*Im(t)).
+ FIXME: this part of code is not 100% rigorous, since we don't consider
+ rounding errors.
+ */
+ mpc_init2 (u, 64);
+ mpfr_const_pi (MPC_RE(u), GMP_RNDN);
+ mpfr_div_2exp (MPC_RE(u), MPC_RE(u), 1, GMP_RNDN); /* Pi/2 */
+ mpfr_remquo (MPC_RE(u), &Q, MPC_IM(t), MPC_RE(u), GMP_RNDN);
+ if (mpfr_sgn (MPC_RE(u)) < 0)
+ Q--; /* corresponds to positive remainder */
+ mpfr_set_ui (MPC_RE(z), 0, GMP_RNDN);
+ mpfr_set_ui (MPC_IM(z), 0, GMP_RNDN);
+ switch (Q & 3)
+ {
+ case 0: /* first quadrant: round to (+0 +0) */
+ ret = MPC_INEX(-1, -1);
+ break;
+ case 1: /* second quadrant: round to (-0 +0) */
+ mpfr_neg (MPC_RE(z), MPC_RE(z), GMP_RNDN);
+ ret = MPC_INEX(1, -1);
+ break;
+ case 2: /* third quadrant: round to (-0 -0) */
+ mpfr_neg (MPC_RE(z), MPC_RE(z), GMP_RNDN);
+ mpfr_neg (MPC_IM(z), MPC_IM(z), GMP_RNDN);
+ ret = MPC_INEX(1, 1);
+ break;
+ case 3: /* fourth quadrant: round to (+0 -0) */
+ mpfr_neg (MPC_IM(z), MPC_IM(z), GMP_RNDN);
+ ret = MPC_INEX(-1, 1);
+ break;
+ }
+ goto clear_t_and_u;
+
+ overflow:
+ /* If we have an overflow, we know that |z| is too large to be
+ represented, but depending on arg(z), we should return +/-Inf +/- I*Inf.
+ We assume t is the approximation of y*log(x), thus we want
+ exp(t) = exp(Re(t))+exp(I*Im(t)).
+ FIXME: this part of code is not 100% rigorous, since we don't consider
+ rounding errors.
+ */
+ mpc_init2 (u, 64);
+ mpfr_const_pi (MPC_RE(u), GMP_RNDN);
+ mpfr_div_2exp (MPC_RE(u), MPC_RE(u), 1, GMP_RNDN); /* Pi/2 */
+ /* the quotient is rounded to the nearest integer in mpfr_remquo */
+ mpfr_remquo (MPC_RE(u), &Q, MPC_IM(t), MPC_RE(u), GMP_RNDN);
+ if (mpfr_sgn (MPC_RE(u)) < 0)
+ Q--; /* corresponds to positive remainder */
+ switch (Q & 3)
+ {
+ case 0: /* first quadrant */
+ mpfr_set_inf (MPC_RE(z), 1);
+ mpfr_set_inf (MPC_IM(z), 1);
+ ret = MPC_INEX(1, 1);
+ break;
+ case 1: /* second quadrant */
+ mpfr_set_inf (MPC_RE(z), -1);
+ mpfr_set_inf (MPC_IM(z), 1);
+ ret = MPC_INEX(-1, 1);
+ break;
+ case 2: /* third quadrant */
+ mpfr_set_inf (MPC_RE(z), -1);
+ mpfr_set_inf (MPC_IM(z), -1);
+ ret = MPC_INEX(-1, -1);
+ break;
+ case 3: /* fourth quadrant */
+ mpfr_set_inf (MPC_RE(z), 1);
+ mpfr_set_inf (MPC_IM(z), -1);
+ ret = MPC_INEX(1, -1);
+ break;
+ }
+
+ clear_t_and_u:
+ mpc_clear (t);
+ mpc_clear (u);
+ return ret;
}
projects / platform / upstream / mpc.git / blobdiff