|
Packit Service |
2e9770 |
/* mpc_pow -- Raise a complex number to the power of another complex number.
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
Copyright (C) 2009-2015 INRIA
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
This file is part of GNU MPC.
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
GNU MPC is free software; you can redistribute it and/or modify it under
|
|
Packit Service |
2e9770 |
the terms of the GNU Lesser General Public License as published by the
|
|
Packit Service |
2e9770 |
Free Software Foundation; either version 3 of the License, or (at your
|
|
Packit Service |
2e9770 |
option) any later version.
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
GNU MPC is distributed in the hope that it will be useful, but WITHOUT ANY
|
|
Packit Service |
2e9770 |
WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
|
|
Packit Service |
2e9770 |
FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
|
|
Packit Service |
2e9770 |
more details.
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
You should have received a copy of the GNU Lesser General Public License
|
|
Packit Service |
2e9770 |
along with this program. If not, see http://www.gnu.org/licenses/ .
|
|
Packit Service |
2e9770 |
*/
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
#include <stdio.h> /* for MPC_ASSERT */
|
|
Packit Service |
2e9770 |
#include "mpc-impl.h"
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* Return non-zero iff c+i*d is an exact square (a+i*b)^2,
|
|
Packit Service |
2e9770 |
with a, b both of the form m*2^e with m, e integers.
|
|
Packit Service |
2e9770 |
If so, returns in a+i*b the corresponding square root, with a >= 0.
|
|
Packit Service |
2e9770 |
The variables a, b must not overlap with c, d.
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
We have c = a^2 - b^2 and d = 2*a*b.
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
If one of a, b is exact, then both are (see algorithms.tex).
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
Case 1: a <> 0 and b <> 0.
|
|
Packit Service |
2e9770 |
Let a = m*2^e and b = n*2^f with m, e, n, f integers, m and n odd
|
|
Packit Service |
2e9770 |
(we will treat apart the case a = 0 or b = 0).
|
|
Packit Service |
2e9770 |
Then 2*a*b = m*n*2^(e+f+1), thus necessarily e+f >= -1.
|
|
Packit Service |
2e9770 |
Assume e < 0, then f >= 0, then a^2 - b^2 = m^2*2^(2e) - n^2*2^(2f) cannot
|
|
Packit Service |
2e9770 |
be an integer, since n^2*2^(2f) is an integer, and m^2*2^(2e) is not.
|
|
Packit Service |
2e9770 |
Similarly when f < 0 (and thus e >= 0).
|
|
Packit Service |
2e9770 |
Thus we have e, f >= 0, and a, b are both integers.
|
|
Packit Service |
2e9770 |
Let A = 2a^2, then eliminating b between c = a^2 - b^2 and d = 2*a*b
|
|
Packit Service |
2e9770 |
gives A^2 - 2c*A - d^2 = 0, which has solutions c +/- sqrt(c^2+d^2).
|
|
Packit Service |
2e9770 |
We thus need c^2+d^2 to be a square, and c + sqrt(c^2+d^2) --- the solution
|
|
Packit Service |
2e9770 |
we are interested in --- to be two times a square. Then b = d/(2a) is
|
|
Packit Service |
2e9770 |
necessarily an integer.
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
Case 2: a = 0. Then d is necessarily zero, thus it suffices to check
|
|
Packit Service |
2e9770 |
whether c = -b^2, i.e., if -c is a square.
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
Case 3: b = 0. Then d is necessarily zero, thus it suffices to check
|
|
Packit Service |
2e9770 |
whether c = a^2, i.e., if c is a square.
|
|
Packit Service |
2e9770 |
*/
|
|
Packit Service |
2e9770 |
static int
|
|
Packit Service |
2e9770 |
mpc_perfect_square_p (mpz_t a, mpz_t b, mpz_t c, mpz_t d)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
if (mpz_cmp_ui (d, 0) == 0) /* case a = 0 or b = 0 */
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
/* necessarily c < 0 here, since we have already considered the case
|
|
Packit Service |
2e9770 |
where x is real non-negative and y is real */
|
|
Packit Service |
2e9770 |
MPC_ASSERT (mpz_cmp_ui (c, 0) < 0);
|
|
Packit Service |
2e9770 |
mpz_neg (b, c);
|
|
Packit Service |
2e9770 |
if (mpz_perfect_square_p (b)) /* case 2 above */
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpz_sqrt (b, b);
|
|
Packit Service |
2e9770 |
mpz_set_ui (a, 0);
|
|
Packit Service |
2e9770 |
return 1; /* c + i*d = (0 + i*b)^2 */
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
else /* both a and b are non-zero */
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
if (mpz_divisible_2exp_p (d, 1) == 0)
|
|
Packit Service |
2e9770 |
return 0; /* d must be even */
|
|
Packit Service |
2e9770 |
mpz_mul (a, c, c);
|
|
Packit Service |
2e9770 |
mpz_addmul (a, d, d); /* c^2 + d^2 */
|
|
Packit Service |
2e9770 |
if (mpz_perfect_square_p (a))
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpz_sqrt (a, a);
|
|
Packit Service |
2e9770 |
mpz_add (a, c, a); /* c + sqrt(c^2+d^2) */
|
|
Packit Service |
2e9770 |
if (mpz_divisible_2exp_p (a, 1))
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpz_tdiv_q_2exp (a, a, 1);
|
|
Packit Service |
2e9770 |
if (mpz_perfect_square_p (a))
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpz_sqrt (a, a);
|
|
Packit Service |
2e9770 |
mpz_tdiv_q_2exp (b, d, 1); /* d/2 */
|
|
Packit Service |
2e9770 |
mpz_divexact (b, b, a); /* d/(2a) */
|
|
Packit Service |
2e9770 |
return 1;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
return 0; /* not a square */
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* fix the sign of Re(z) or Im(z) in case it is zero,
|
|
Packit Service |
2e9770 |
and Re(x) is zero.
|
|
Packit Service |
2e9770 |
sign_eps is 0 if Re(x) = +0, 1 if Re(x) = -0
|
|
Packit Service |
2e9770 |
sign_a is the sign bit of Im(x).
|
|
Packit Service |
2e9770 |
Assume y is an integer (does nothing otherwise).
|
|
Packit Service |
2e9770 |
*/
|
|
Packit Service |
2e9770 |
static void
|
|
Packit Service |
2e9770 |
fix_sign (mpc_ptr z, int sign_eps, int sign_a, mpfr_srcptr y)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
int ymod4 = -1;
|
|
Packit Service |
2e9770 |
mpfr_exp_t ey;
|
|
Packit Service |
2e9770 |
mpz_t my;
|
|
Packit Service |
2e9770 |
unsigned long int t;
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
mpz_init (my);
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
ey = mpfr_get_z_exp (my, y);
|
|
Packit Service |
2e9770 |
/* normalize so that my is odd */
|
|
Packit Service |
2e9770 |
t = mpz_scan1 (my, 0);
|
|
Packit Service |
2e9770 |
ey += (mpfr_exp_t) t;
|
|
Packit Service |
2e9770 |
mpz_tdiv_q_2exp (my, my, t);
|
|
Packit Service |
2e9770 |
/* y = my*2^ey */
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* compute y mod 4 (in case y is an integer) */
|
|
Packit Service |
2e9770 |
if (ey >= 2)
|
|
Packit Service |
2e9770 |
ymod4 = 0;
|
|
Packit Service |
2e9770 |
else if (ey == 1)
|
|
Packit Service |
2e9770 |
ymod4 = mpz_tstbit (my, 0) * 2; /* correct if my < 0 */
|
|
Packit Service |
2e9770 |
else if (ey == 0)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
ymod4 = mpz_tstbit (my, 1) * 2 + mpz_tstbit (my, 0);
|
|
Packit Service |
2e9770 |
if (mpz_cmp_ui (my , 0) < 0)
|
|
Packit Service |
2e9770 |
ymod4 = 4 - ymod4;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
else /* y is not an integer */
|
|
Packit Service |
2e9770 |
goto end;
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
if (mpfr_zero_p (mpc_realref(z)))
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
/* we assume y is always integer in that case (FIXME: prove it):
|
|
Packit Service |
2e9770 |
(eps+I*a)^y = +0 + I*a^y for y = 1 mod 4 and sign_eps = 0
|
|
Packit Service |
2e9770 |
(eps+I*a)^y = -0 - I*a^y for y = 3 mod 4 and sign_eps = 0 */
|
|
Packit Service |
2e9770 |
MPC_ASSERT (ymod4 == 1 || ymod4 == 3);
|
|
Packit Service |
2e9770 |
if ((ymod4 == 3 && sign_eps == 0) ||
|
|
Packit Service |
2e9770 |
(ymod4 == 1 && sign_eps == 1))
|
|
Packit Service |
2e9770 |
mpfr_neg (mpc_realref(z), mpc_realref(z), MPFR_RNDZ);
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
else if (mpfr_zero_p (mpc_imagref(z)))
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
/* we assume y is always integer in that case (FIXME: prove it):
|
|
Packit Service |
2e9770 |
(eps+I*a)^y = a^y - 0*I for y = 0 mod 4 and sign_a = sign_eps
|
|
Packit Service |
2e9770 |
(eps+I*a)^y = -a^y +0*I for y = 2 mod 4 and sign_a = sign_eps */
|
|
Packit Service |
2e9770 |
MPC_ASSERT (ymod4 == 0 || ymod4 == 2);
|
|
Packit Service |
2e9770 |
if ((ymod4 == 0 && sign_a == sign_eps) ||
|
|
Packit Service |
2e9770 |
(ymod4 == 2 && sign_a != sign_eps))
|
|
Packit Service |
2e9770 |
mpfr_neg (mpc_imagref(z), mpc_imagref(z), MPFR_RNDZ);
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
end:
|
|
Packit Service |
2e9770 |
mpz_clear (my);
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* If x^y is exactly representable (with maybe a larger precision than z),
|
|
Packit Service |
2e9770 |
round it in z and return the (mpc) inexact flag in [0, 10].
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
If x^y is not exactly representable, return -1.
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
If intermediate computations lead to numbers of more than maxprec bits,
|
|
Packit Service |
2e9770 |
then abort and return -2 (in that case, to avoid loops, mpc_pow_exact
|
|
Packit Service |
2e9770 |
should be called again with a larger value of maxprec).
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
Assume one of Re(x) or Im(x) is non-zero, and y is non-zero (y is real).
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
Warning: z and x might be the same variable, same for Re(z) or Im(z) and y.
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
In case -1 or -2 is returned, z is not modified.
|
|
Packit Service |
2e9770 |
*/
|
|
Packit Service |
2e9770 |
static int
|
|
Packit Service |
2e9770 |
mpc_pow_exact (mpc_ptr z, mpc_srcptr x, mpfr_srcptr y, mpc_rnd_t rnd,
|
|
Packit Service |
2e9770 |
mpfr_prec_t maxprec)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpfr_exp_t ec, ed, ey;
|
|
Packit Service |
2e9770 |
mpz_t my, a, b, c, d, u;
|
|
Packit Service |
2e9770 |
unsigned long int t;
|
|
Packit Service |
2e9770 |
int ret = -2;
|
|
Packit Service |
2e9770 |
int sign_rex = mpfr_signbit (mpc_realref(x));
|
|
Packit Service |
2e9770 |
int sign_imx = mpfr_signbit (mpc_imagref(x));
|
|
Packit Service |
2e9770 |
int x_imag = mpfr_zero_p (mpc_realref(x));
|
|
Packit Service |
2e9770 |
int z_is_y = 0;
|
|
Packit Service |
2e9770 |
mpfr_t copy_of_y;
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
if (mpc_realref (z) == y || mpc_imagref (z) == y)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
z_is_y = 1;
|
|
Packit Service |
2e9770 |
mpfr_init2 (copy_of_y, mpfr_get_prec (y));
|
|
Packit Service |
2e9770 |
mpfr_set (copy_of_y, y, MPFR_RNDN);
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
mpz_init (my);
|
|
Packit Service |
2e9770 |
mpz_init (a);
|
|
Packit Service |
2e9770 |
mpz_init (b);
|
|
Packit Service |
2e9770 |
mpz_init (c);
|
|
Packit Service |
2e9770 |
mpz_init (d);
|
|
Packit Service |
2e9770 |
mpz_init (u);
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
ey = mpfr_get_z_exp (my, y);
|
|
Packit Service |
2e9770 |
/* normalize so that my is odd */
|
|
Packit Service |
2e9770 |
t = mpz_scan1 (my, 0);
|
|
Packit Service |
2e9770 |
ey += (mpfr_exp_t) t;
|
|
Packit Service |
2e9770 |
mpz_tdiv_q_2exp (my, my, t);
|
|
Packit Service |
2e9770 |
/* y = my*2^ey with my odd */
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
if (x_imag)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpz_set_ui (c, 0);
|
|
Packit Service |
2e9770 |
ec = 0;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
else
|
|
Packit Service |
2e9770 |
ec = mpfr_get_z_exp (c, mpc_realref(x));
|
|
Packit Service |
2e9770 |
if (mpfr_zero_p (mpc_imagref(x)))
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpz_set_ui (d, 0);
|
|
Packit Service |
2e9770 |
ed = ec;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
else
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
ed = mpfr_get_z_exp (d, mpc_imagref(x));
|
|
Packit Service |
2e9770 |
if (x_imag)
|
|
Packit Service |
2e9770 |
ec = ed;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
/* x = c*2^ec + I * d*2^ed */
|
|
Packit Service |
2e9770 |
/* equalize the exponents of x */
|
|
Packit Service |
2e9770 |
if (ec < ed)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpz_mul_2exp (d, d, (unsigned long int) (ed - ec));
|
|
Packit Service |
2e9770 |
if ((mpfr_prec_t) mpz_sizeinbase (d, 2) > maxprec)
|
|
Packit Service |
2e9770 |
goto end;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
else if (ed < ec)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpz_mul_2exp (c, c, (unsigned long int) (ec - ed));
|
|
Packit Service |
2e9770 |
if ((mpfr_prec_t) mpz_sizeinbase (c, 2) > maxprec)
|
|
Packit Service |
2e9770 |
goto end;
|
|
Packit Service |
2e9770 |
ec = ed;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
/* now ec=ed and x = (c + I * d) * 2^ec */
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* divide by two if possible */
|
|
Packit Service |
2e9770 |
if (mpz_cmp_ui (c, 0) == 0)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
t = mpz_scan1 (d, 0);
|
|
Packit Service |
2e9770 |
mpz_tdiv_q_2exp (d, d, t);
|
|
Packit Service |
2e9770 |
ec += (mpfr_exp_t) t;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
else if (mpz_cmp_ui (d, 0) == 0)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
t = mpz_scan1 (c, 0);
|
|
Packit Service |
2e9770 |
mpz_tdiv_q_2exp (c, c, t);
|
|
Packit Service |
2e9770 |
ec += (mpfr_exp_t) t;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
else /* neither c nor d is zero */
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
unsigned long v;
|
|
Packit Service |
2e9770 |
t = mpz_scan1 (c, 0);
|
|
Packit Service |
2e9770 |
v = mpz_scan1 (d, 0);
|
|
Packit Service |
2e9770 |
if (v < t)
|
|
Packit Service |
2e9770 |
t = v;
|
|
Packit Service |
2e9770 |
mpz_tdiv_q_2exp (c, c, t);
|
|
Packit Service |
2e9770 |
mpz_tdiv_q_2exp (d, d, t);
|
|
Packit Service |
2e9770 |
ec += (mpfr_exp_t) t;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* now either one of c, d is odd */
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
while (ey < 0)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
/* check if x is a square */
|
|
Packit Service |
2e9770 |
if (ec & 1)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpz_mul_2exp (c, c, 1);
|
|
Packit Service |
2e9770 |
mpz_mul_2exp (d, d, 1);
|
|
Packit Service |
2e9770 |
ec --;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
/* now ec is even */
|
|
Packit Service |
2e9770 |
if (mpc_perfect_square_p (a, b, c, d) == 0)
|
|
Packit Service |
2e9770 |
break;
|
|
Packit Service |
2e9770 |
mpz_swap (a, c);
|
|
Packit Service |
2e9770 |
mpz_swap (b, d);
|
|
Packit Service |
2e9770 |
ec /= 2;
|
|
Packit Service |
2e9770 |
ey ++;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
if (ey < 0)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
ret = -1; /* not representable */
|
|
Packit Service |
2e9770 |
goto end;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* Now ey >= 0, it thus suffices to check that x^my is representable.
|
|
Packit Service |
2e9770 |
If my > 0, this is always true. If my < 0, we first try to invert
|
|
Packit Service |
2e9770 |
(c+I*d)*2^ec.
|
|
Packit Service |
2e9770 |
*/
|
|
Packit Service |
2e9770 |
if (mpz_cmp_ui (my, 0) < 0)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
/* If my < 0, 1 / (c + I*d) = (c - I*d)/(c^2 + d^2), thus a sufficient
|
|
Packit Service |
2e9770 |
condition is that c^2 + d^2 is a power of two, assuming |c| <> |d|.
|
|
Packit Service |
2e9770 |
Assume a prime p <> 2 divides c^2 + d^2,
|
|
Packit Service |
2e9770 |
then if p does not divide c or d, 1 / (c + I*d) cannot be exact.
|
|
Packit Service |
2e9770 |
If p divides both c and d, then we can write c = p*c', d = p*d',
|
|
Packit Service |
2e9770 |
and 1 / (c + I*d) = 1/p * 1/(c' + I*d'). This shows that if 1/(c+I*d)
|
|
Packit Service |
2e9770 |
is exact, then 1/(c' + I*d') is exact too, and we are back to the
|
|
Packit Service |
2e9770 |
previous case. In conclusion, a necessary and sufficient condition
|
|
Packit Service |
2e9770 |
is that c^2 + d^2 is a power of two.
|
|
Packit Service |
2e9770 |
*/
|
|
Packit Service |
2e9770 |
/* FIXME: we could first compute c^2+d^2 mod a limb for example */
|
|
Packit Service |
2e9770 |
mpz_mul (a, c, c);
|
|
Packit Service |
2e9770 |
mpz_addmul (a, d, d);
|
|
Packit Service |
2e9770 |
t = mpz_scan1 (a, 0);
|
|
Packit Service |
2e9770 |
if (mpz_sizeinbase (a, 2) != 1 + t) /* a is not a power of two */
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
ret = -1; /* not representable */
|
|
Packit Service |
2e9770 |
goto end;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
/* replace (c,d) by (c/(c^2+d^2), -d/(c^2+d^2)) */
|
|
Packit Service |
2e9770 |
mpz_neg (d, d);
|
|
Packit Service |
2e9770 |
ec = -ec - (mpfr_exp_t) t;
|
|
Packit Service |
2e9770 |
mpz_neg (my, my);
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* now ey >= 0 and my >= 0, and we want to compute
|
|
Packit Service |
2e9770 |
[(c + I * d) * 2^ec] ^ (my * 2^ey).
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
We first compute [(c + I * d) * 2^ec]^my, then square ey times. */
|
|
Packit Service |
2e9770 |
t = mpz_sizeinbase (my, 2) - 1;
|
|
Packit Service |
2e9770 |
mpz_set (a, c);
|
|
Packit Service |
2e9770 |
mpz_set (b, d);
|
|
Packit Service |
2e9770 |
ed = ec;
|
|
Packit Service |
2e9770 |
/* invariant: (a + I*b) * 2^ed = ((c + I*d) * 2^ec)^trunc(my/2^t) */
|
|
Packit Service |
2e9770 |
while (t-- > 0)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
unsigned long int v, w;
|
|
Packit Service |
2e9770 |
/* square a + I*b */
|
|
Packit Service |
2e9770 |
mpz_mul (u, a, b);
|
|
Packit Service |
2e9770 |
mpz_mul (a, a, a);
|
|
Packit Service |
2e9770 |
mpz_submul (a, b, b);
|
|
Packit Service |
2e9770 |
mpz_mul_2exp (b, u, 1);
|
|
Packit Service |
2e9770 |
ed *= 2;
|
|
Packit Service |
2e9770 |
if (mpz_tstbit (my, t)) /* multiply by c + I*d */
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpz_mul (u, a, c);
|
|
Packit Service |
2e9770 |
mpz_submul (u, b, d); /* ac-bd */
|
|
Packit Service |
2e9770 |
mpz_mul (b, b, c);
|
|
Packit Service |
2e9770 |
mpz_addmul (b, a, d); /* bc+ad */
|
|
Packit Service |
2e9770 |
mpz_swap (a, u);
|
|
Packit Service |
2e9770 |
ed += ec;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
/* remove powers of two in (a,b) */
|
|
Packit Service |
2e9770 |
if (mpz_cmp_ui (a, 0) == 0)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
w = mpz_scan1 (b, 0);
|
|
Packit Service |
2e9770 |
mpz_tdiv_q_2exp (b, b, w);
|
|
Packit Service |
2e9770 |
ed += (mpfr_exp_t) w;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
else if (mpz_cmp_ui (b, 0) == 0)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
w = mpz_scan1 (a, 0);
|
|
Packit Service |
2e9770 |
mpz_tdiv_q_2exp (a, a, w);
|
|
Packit Service |
2e9770 |
ed += (mpfr_exp_t) w;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
else
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
w = mpz_scan1 (a, 0);
|
|
Packit Service |
2e9770 |
v = mpz_scan1 (b, 0);
|
|
Packit Service |
2e9770 |
if (v < w)
|
|
Packit Service |
2e9770 |
w = v;
|
|
Packit Service |
2e9770 |
mpz_tdiv_q_2exp (a, a, w);
|
|
Packit Service |
2e9770 |
mpz_tdiv_q_2exp (b, b, w);
|
|
Packit Service |
2e9770 |
ed += (mpfr_exp_t) w;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
if ( (mpfr_prec_t) mpz_sizeinbase (a, 2) > maxprec
|
|
Packit Service |
2e9770 |
|| (mpfr_prec_t) mpz_sizeinbase (b, 2) > maxprec)
|
|
Packit Service |
2e9770 |
goto end;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
/* now a+I*b = (c+I*d)^my */
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
while (ey-- > 0)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
unsigned long sa, sb;
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* square a + I*b */
|
|
Packit Service |
2e9770 |
mpz_mul (u, a, b);
|
|
Packit Service |
2e9770 |
mpz_mul (a, a, a);
|
|
Packit Service |
2e9770 |
mpz_submul (a, b, b);
|
|
Packit Service |
2e9770 |
mpz_mul_2exp (b, u, 1);
|
|
Packit Service |
2e9770 |
ed *= 2;
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* divide by largest 2^n possible, to avoid many loops for e.g.,
|
|
Packit Service |
2e9770 |
(2+2*I)^16777216 */
|
|
Packit Service |
2e9770 |
sa = mpz_scan1 (a, 0);
|
|
Packit Service |
2e9770 |
sb = mpz_scan1 (b, 0);
|
|
Packit Service |
2e9770 |
sa = (sa <= sb) ? sa : sb;
|
|
Packit Service |
2e9770 |
mpz_tdiv_q_2exp (a, a, sa);
|
|
Packit Service |
2e9770 |
mpz_tdiv_q_2exp (b, b, sa);
|
|
Packit Service |
2e9770 |
ed += (mpfr_exp_t) sa;
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
if ( (mpfr_prec_t) mpz_sizeinbase (a, 2) > maxprec
|
|
Packit Service |
2e9770 |
|| (mpfr_prec_t) mpz_sizeinbase (b, 2) > maxprec)
|
|
Packit Service |
2e9770 |
goto end;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
ret = mpfr_set_z (mpc_realref(z), a, MPC_RND_RE(rnd));
|
|
Packit Service |
2e9770 |
ret = MPC_INEX(ret, mpfr_set_z (mpc_imagref(z), b, MPC_RND_IM(rnd)));
|
|
Packit Service |
2e9770 |
mpfr_mul_2si (mpc_realref(z), mpc_realref(z), ed, MPC_RND_RE(rnd));
|
|
Packit Service |
2e9770 |
mpfr_mul_2si (mpc_imagref(z), mpc_imagref(z), ed, MPC_RND_IM(rnd));
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
end:
|
|
Packit Service |
2e9770 |
mpz_clear (my);
|
|
Packit Service |
2e9770 |
mpz_clear (a);
|
|
Packit Service |
2e9770 |
mpz_clear (b);
|
|
Packit Service |
2e9770 |
mpz_clear (c);
|
|
Packit Service |
2e9770 |
mpz_clear (d);
|
|
Packit Service |
2e9770 |
mpz_clear (u);
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
if (ret >= 0 && x_imag)
|
|
Packit Service |
2e9770 |
fix_sign (z, sign_rex, sign_imx, (z_is_y) ? copy_of_y : y);
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
if (z_is_y)
|
|
Packit Service |
2e9770 |
mpfr_clear (copy_of_y);
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
return ret;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* Return 1 if y*2^k is an odd integer, 0 otherwise.
|
|
Packit Service |
2e9770 |
Adapted from MPFR, file pow.c.
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
Examples: with k=0, check if y is an odd integer,
|
|
Packit Service |
2e9770 |
with k=1, check if y is half-an-integer,
|
|
Packit Service |
2e9770 |
with k=-1, check if y/2 is an odd integer.
|
|
Packit Service |
2e9770 |
*/
|
|
Packit Service |
2e9770 |
#define MPFR_LIMB_HIGHBIT ((mp_limb_t) 1 << (BITS_PER_MP_LIMB - 1))
|
|
Packit Service |
2e9770 |
static int
|
|
Packit Service |
2e9770 |
is_odd (mpfr_srcptr y, mpfr_exp_t k)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpfr_exp_t expo;
|
|
Packit Service |
2e9770 |
mpfr_prec_t prec;
|
|
Packit Service |
2e9770 |
mp_size_t yn;
|
|
Packit Service |
2e9770 |
mp_limb_t *yp;
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
expo = mpfr_get_exp (y) + k;
|
|
Packit Service |
2e9770 |
if (expo <= 0)
|
|
Packit Service |
2e9770 |
return 0; /* |y| < 1 and not 0 */
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
prec = mpfr_get_prec (y);
|
|
Packit Service |
2e9770 |
if ((mpfr_prec_t) expo > prec)
|
|
Packit Service |
2e9770 |
return 0; /* y is a multiple of 2^(expo-prec), thus not odd */
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* 0 < expo <= prec:
|
|
Packit Service |
2e9770 |
y = 1xxxxxxxxxt.zzzzzzzzzzzzzzzzzz[000]
|
|
Packit Service |
2e9770 |
expo bits (prec-expo) bits
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
We have to check that:
|
|
Packit Service |
2e9770 |
(a) the bit 't' is set
|
|
Packit Service |
2e9770 |
(b) all the 'z' bits are zero
|
|
Packit Service |
2e9770 |
*/
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
prec = ((prec - 1) / BITS_PER_MP_LIMB + 1) * BITS_PER_MP_LIMB - expo;
|
|
Packit Service |
2e9770 |
/* number of z+0 bits */
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
yn = prec / BITS_PER_MP_LIMB;
|
|
Packit Service |
2e9770 |
/* yn is the index of limb containing the 't' bit */
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
yp = y->_mpfr_d;
|
|
Packit Service |
2e9770 |
/* if expo is a multiple of BITS_PER_MP_LIMB, t is bit 0 */
|
|
Packit Service |
2e9770 |
if (expo % BITS_PER_MP_LIMB == 0 ? (yp[yn] & 1) == 0
|
|
Packit Service |
2e9770 |
: yp[yn] << ((expo % BITS_PER_MP_LIMB) - 1) != MPFR_LIMB_HIGHBIT)
|
|
Packit Service |
2e9770 |
return 0;
|
|
Packit Service |
2e9770 |
while (--yn >= 0)
|
|
Packit Service |
2e9770 |
if (yp[yn] != 0)
|
|
Packit Service |
2e9770 |
return 0;
|
|
Packit Service |
2e9770 |
return 1;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* Put in z the value of x^y, rounded according to 'rnd'.
|
|
Packit Service |
2e9770 |
Return the inexact flag in [0, 10]. */
|
|
Packit Service |
2e9770 |
int
|
|
Packit Service |
2e9770 |
mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
int ret = -2, loop, x_real, x_imag, y_real, z_real = 0, z_imag = 0;
|
|
Packit Service |
2e9770 |
mpc_t t, u;
|
|
Packit Service |
2e9770 |
mpfr_prec_t p, pr, pi, maxprec;
|
|
Packit Service |
2e9770 |
int saved_underflow, saved_overflow;
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* save the underflow or overflow flags from MPFR */
|
|
Packit Service |
2e9770 |
saved_underflow = mpfr_underflow_p ();
|
|
Packit Service |
2e9770 |
saved_overflow = mpfr_overflow_p ();
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
x_real = mpfr_zero_p (mpc_imagref(x));
|
|
Packit Service |
2e9770 |
y_real = mpfr_zero_p (mpc_imagref(y));
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
if (y_real && mpfr_zero_p (mpc_realref(y))) /* case y zero */
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
if (x_real && mpfr_zero_p (mpc_realref(x)))
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
/* we define 0^0 to be (1, +0) since the real part is
|
|
Packit Service |
2e9770 |
coherent with MPFR where 0^0 gives 1, and the sign of the
|
|
Packit Service |
2e9770 |
imaginary part cannot be determined */
|
|
Packit Service |
2e9770 |
mpc_set_ui_ui (z, 1, 0, MPC_RNDNN);
|
|
Packit Service |
2e9770 |
return 0;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
else /* x^0 = 1 +/- i*0 even for x=NaN see algorithms.tex for the
|
|
Packit Service |
2e9770 |
sign of zero */
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpfr_t n;
|
|
Packit Service |
2e9770 |
int inex, cx1;
|
|
Packit Service |
2e9770 |
int sign_zi;
|
|
Packit Service |
2e9770 |
/* cx1 < 0 if |x| < 1
|
|
Packit Service |
2e9770 |
cx1 = 0 if |x| = 1
|
|
Packit Service |
2e9770 |
cx1 > 0 if |x| > 1
|
|
Packit Service |
2e9770 |
*/
|
|
Packit Service |
2e9770 |
mpfr_init (n);
|
|
Packit Service |
2e9770 |
inex = mpc_norm (n, x, MPFR_RNDN);
|
|
Packit Service |
2e9770 |
cx1 = mpfr_cmp_ui (n, 1);
|
|
Packit Service |
2e9770 |
if (cx1 == 0 && inex != 0)
|
|
Packit Service |
2e9770 |
cx1 = -inex;
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
sign_zi = (cx1 < 0 && mpfr_signbit (mpc_imagref (y)) == 0)
|
|
Packit Service |
2e9770 |
|| (cx1 == 0
|
|
Packit Service |
2e9770 |
&& mpfr_signbit (mpc_imagref (x)) != mpfr_signbit (mpc_realref (y)))
|
|
Packit Service |
2e9770 |
|| (cx1 > 0 && mpfr_signbit (mpc_imagref (y)));
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* warning: mpc_set_ui_ui does not set Im(z) to -0 if Im(rnd)=RNDD */
|
|
Packit Service |
2e9770 |
ret = mpc_set_ui_ui (z, 1, 0, rnd);
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
if (MPC_RND_IM (rnd) == MPFR_RNDD || sign_zi)
|
|
Packit Service |
2e9770 |
mpc_conj (z, z, MPC_RNDNN);
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
mpfr_clear (n);
|
|
Packit Service |
2e9770 |
return ret;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
if (!mpc_fin_p (x) || !mpc_fin_p (y))
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
/* special values: exp(y*log(x)) */
|
|
Packit Service |
2e9770 |
mpc_init2 (u, 2);
|
|
Packit Service |
2e9770 |
mpc_log (u, x, MPC_RNDNN);
|
|
Packit Service |
2e9770 |
mpc_mul (u, u, y, MPC_RNDNN);
|
|
Packit Service |
2e9770 |
ret = mpc_exp (z, u, rnd);
|
|
Packit Service |
2e9770 |
mpc_clear (u);
|
|
Packit Service |
2e9770 |
goto end;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
if (x_real) /* case x real */
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
if (mpfr_zero_p (mpc_realref(x))) /* x is zero */
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
/* special values: exp(y*log(x)) */
|
|
Packit Service |
2e9770 |
mpc_init2 (u, 2);
|
|
Packit Service |
2e9770 |
mpc_log (u, x, MPC_RNDNN);
|
|
Packit Service |
2e9770 |
mpc_mul (u, u, y, MPC_RNDNN);
|
|
Packit Service |
2e9770 |
ret = mpc_exp (z, u, rnd);
|
|
Packit Service |
2e9770 |
mpc_clear (u);
|
|
Packit Service |
2e9770 |
goto end;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* Special case 1^y = 1 */
|
|
Packit Service |
2e9770 |
if (mpfr_cmp_ui (mpc_realref(x), 1) == 0)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
int s1, s2;
|
|
Packit Service |
2e9770 |
s1 = mpfr_signbit (mpc_realref (y));
|
|
Packit Service |
2e9770 |
s2 = mpfr_signbit (mpc_imagref (x));
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
ret = mpc_set_ui (z, +1, rnd);
|
|
Packit Service |
2e9770 |
/* the sign of the zero imaginary part is known in some cases (see
|
|
Packit Service |
2e9770 |
algorithm.tex). In such cases we have
|
|
Packit Service |
2e9770 |
(x +s*0i)^(y+/-0i) = x^y + s*sign(y)*0i
|
|
Packit Service |
2e9770 |
where s = +/-1. We extend here this rule to fix the sign of the
|
|
Packit Service |
2e9770 |
zero part.
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
Note that the sign must also be set explicitly when rnd=RNDD
|
|
Packit Service |
2e9770 |
because mpfr_set_ui(z_i, 0, rnd) always sets z_i to +0.
|
|
Packit Service |
2e9770 |
*/
|
|
Packit Service |
2e9770 |
if (MPC_RND_IM (rnd) == MPFR_RNDD || s1 != s2)
|
|
Packit Service |
2e9770 |
mpc_conj (z, z, MPC_RNDNN);
|
|
Packit Service |
2e9770 |
goto end;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* x^y is real when:
|
|
Packit Service |
2e9770 |
(a) x is real and y is integer
|
|
Packit Service |
2e9770 |
(b) x is real non-negative and y is real */
|
|
Packit Service |
2e9770 |
if (y_real && (mpfr_integer_p (mpc_realref(y)) ||
|
|
Packit Service |
2e9770 |
mpfr_cmp_ui (mpc_realref(x), 0) >= 0))
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
int s1, s2;
|
|
Packit Service |
2e9770 |
s1 = mpfr_signbit (mpc_realref (y));
|
|
Packit Service |
2e9770 |
s2 = mpfr_signbit (mpc_imagref (x));
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
ret = mpfr_pow (mpc_realref(z), mpc_realref(x), mpc_realref(y), MPC_RND_RE(rnd));
|
|
Packit Service |
2e9770 |
ret = MPC_INEX(ret, mpfr_set_ui (mpc_imagref(z), 0, MPC_RND_IM(rnd)));
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* the sign of the zero imaginary part is known in some cases
|
|
Packit Service |
2e9770 |
(see algorithm.tex). In such cases we have (x +s*0i)^(y+/-0i)
|
|
Packit Service |
2e9770 |
= x^y + s*sign(y)*0i where s = +/-1.
|
|
Packit Service |
2e9770 |
We extend here this rule to fix the sign of the zero part.
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
Note that the sign must also be set explicitly when rnd=RNDD
|
|
Packit Service |
2e9770 |
because mpfr_set_ui(z_i, 0, rnd) always sets z_i to +0.
|
|
Packit Service |
2e9770 |
*/
|
|
Packit Service |
2e9770 |
if (MPC_RND_IM(rnd) == MPFR_RNDD || s1 != s2)
|
|
Packit Service |
2e9770 |
mpfr_neg (mpc_imagref(z), mpc_imagref(z), MPC_RND_IM(rnd));
|
|
Packit Service |
2e9770 |
goto end;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* (-1)^(n+I*t) is real for n integer and t real */
|
|
Packit Service |
2e9770 |
if (mpfr_cmp_si (mpc_realref(x), -1) == 0 && mpfr_integer_p (mpc_realref(y)))
|
|
Packit Service |
2e9770 |
z_real = 1;
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* for x real, x^y is imaginary when:
|
|
Packit Service |
2e9770 |
(a) x is negative and y is half-an-integer
|
|
Packit Service |
2e9770 |
(b) x = -1 and Re(y) is half-an-integer
|
|
Packit Service |
2e9770 |
*/
|
|
Packit Service |
2e9770 |
if ((mpfr_cmp_ui (mpc_realref(x), 0) < 0) && is_odd (mpc_realref(y), 1)
|
|
Packit Service |
2e9770 |
&& (y_real || mpfr_cmp_si (mpc_realref(x), -1) == 0))
|
|
Packit Service |
2e9770 |
z_imag = 1;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
else /* x non real */
|
|
Packit Service |
2e9770 |
/* I^(t*I) and (-I)^(t*I) are real for t real,
|
|
Packit Service |
2e9770 |
I^(n+t*I) and (-I)^(n+t*I) are real for n even and t real, and
|
|
Packit Service |
2e9770 |
I^(n+t*I) and (-I)^(n+t*I) are imaginary for n odd and t real
|
|
Packit Service |
2e9770 |
(s*I)^n is real for n even and imaginary for n odd */
|
|
Packit Service |
2e9770 |
if ((mpc_cmp_si_si (x, 0, 1) == 0 || mpc_cmp_si_si (x, 0, -1) == 0 ||
|
|
Packit Service |
2e9770 |
(mpfr_cmp_ui (mpc_realref(x), 0) == 0 && y_real)) &&
|
|
Packit Service |
2e9770 |
mpfr_integer_p (mpc_realref(y)))
|
|
Packit Service |
2e9770 |
{ /* x is I or -I, and Re(y) is an integer */
|
|
Packit Service |
2e9770 |
if (is_odd (mpc_realref(y), 0))
|
|
Packit Service |
2e9770 |
z_imag = 1; /* Re(y) odd: z is imaginary */
|
|
Packit Service |
2e9770 |
else
|
|
Packit Service |
2e9770 |
z_real = 1; /* Re(y) even: z is real */
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
else /* (t+/-t*I)^(2n) is imaginary for n odd and real for n even */
|
|
Packit Service |
2e9770 |
if (mpfr_cmpabs (mpc_realref(x), mpc_imagref(x)) == 0 && y_real &&
|
|
Packit Service |
2e9770 |
mpfr_integer_p (mpc_realref(y)) && is_odd (mpc_realref(y), 0) == 0)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
if (is_odd (mpc_realref(y), -1)) /* y/2 is odd */
|
|
Packit Service |
2e9770 |
z_imag = 1;
|
|
Packit Service |
2e9770 |
else
|
|
Packit Service |
2e9770 |
z_real = 1;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
pr = mpfr_get_prec (mpc_realref(z));
|
|
Packit Service |
2e9770 |
pi = mpfr_get_prec (mpc_imagref(z));
|
|
Packit Service |
2e9770 |
p = (pr > pi) ? pr : pi;
|
|
Packit Service |
2e9770 |
p += 12; /* experimentally, seems to give less than 10% of failures in
|
|
Packit Service |
2e9770 |
Ziv's strategy; probably wrong now since q is not computed */
|
|
Packit Service |
2e9770 |
if (p < 64)
|
|
Packit Service |
2e9770 |
p = 64;
|
|
Packit Service |
2e9770 |
mpc_init2 (u, p);
|
|
Packit Service |
2e9770 |
mpc_init2 (t, p);
|
|
Packit Service |
2e9770 |
pr += MPC_RND_RE(rnd) == MPFR_RNDN;
|
|
Packit Service |
2e9770 |
pi += MPC_RND_IM(rnd) == MPFR_RNDN;
|
|
Packit Service |
2e9770 |
maxprec = MPC_MAX_PREC (z);
|
|
Packit Service |
2e9770 |
x_imag = mpfr_zero_p (mpc_realref(x));
|
|
Packit Service |
2e9770 |
for (loop = 0;; loop++)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
int ret_exp;
|
|
Packit Service |
2e9770 |
mpfr_exp_t dr, di;
|
|
Packit Service |
2e9770 |
mpfr_prec_t q;
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
mpc_log (t, x, MPC_RNDNN);
|
|
Packit Service |
2e9770 |
mpc_mul (t, t, y, MPC_RNDNN);
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* Compute q such that |Re (y log x)|, |Im (y log x)| < 2^q.
|
|
Packit Service |
2e9770 |
We recompute it at each loop since we might get different
|
|
Packit Service |
2e9770 |
bounds if the precision is not enough. */
|
|
Packit Service |
2e9770 |
q = mpfr_get_exp (mpc_realref(t)) > 0 ? mpfr_get_exp (mpc_realref(t)) : 0;
|
|
Packit Service |
2e9770 |
if (mpfr_get_exp (mpc_imagref(t)) > (mpfr_exp_t) q)
|
|
Packit Service |
2e9770 |
q = mpfr_get_exp (mpc_imagref(t));
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* if q >= p, we get an error of order 1 on the imaginary part of t,
|
|
Packit Service |
2e9770 |
which is not enough to get the correct sign of exp(t) */
|
|
Packit Service |
2e9770 |
if (q >= p)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
p = p + 64;
|
|
Packit Service |
2e9770 |
goto try_again;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
mpfr_clear_overflow ();
|
|
Packit Service |
2e9770 |
mpfr_clear_underflow ();
|
|
Packit Service |
2e9770 |
ret_exp = mpc_exp (u, t, MPC_RNDNN);
|
|
Packit Service |
2e9770 |
if (mpfr_underflow_p () || mpfr_overflow_p ()) {
|
|
Packit Service |
2e9770 |
int inex_re, inex_im;
|
|
Packit Service |
2e9770 |
/* under- and overflow flags are set by mpc_exp */
|
|
Packit Service |
2e9770 |
mpc_set (z, u, MPC_RNDNN);
|
|
Packit Service |
2e9770 |
inex_re = MPC_INEX_RE(ret_exp);
|
|
Packit Service |
2e9770 |
inex_im = MPC_INEX_IM(ret_exp);
|
|
Packit Service |
2e9770 |
if (mpfr_inf_p (mpc_realref (z)))
|
|
Packit Service |
2e9770 |
inex_re = mpc_fix_inf (mpc_realref (z), MPC_RND_RE(rnd));
|
|
Packit Service |
2e9770 |
if (mpfr_inf_p (mpc_imagref (z)))
|
|
Packit Service |
2e9770 |
inex_im = mpc_fix_inf (mpc_imagref (z), MPC_RND_IM(rnd));
|
|
Packit Service |
2e9770 |
ret = MPC_INEX(inex_re,inex_im);
|
|
Packit Service |
2e9770 |
goto exact;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* Since the error bound is global, we have to take into account the
|
|
Packit Service |
2e9770 |
exponent difference between the real and imaginary parts. We assume
|
|
Packit Service |
2e9770 |
either the real or the imaginary part of u is not zero.
|
|
Packit Service |
2e9770 |
*/
|
|
Packit Service |
2e9770 |
dr = mpfr_zero_p (mpc_realref(u)) ? mpfr_get_exp (mpc_imagref(u))
|
|
Packit Service |
2e9770 |
: mpfr_get_exp (mpc_realref(u));
|
|
Packit Service |
2e9770 |
di = mpfr_zero_p (mpc_imagref(u)) ? dr : mpfr_get_exp (mpc_imagref(u));
|
|
Packit Service |
2e9770 |
if (dr > di)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
di = dr - di;
|
|
Packit Service |
2e9770 |
dr = 0;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
else
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
dr = di - dr;
|
|
Packit Service |
2e9770 |
di = 0;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
/* the term -3 takes into account the factor 4 in the complex error
|
|
Packit Service |
2e9770 |
(see algorithms.tex) plus one due to the exponent difference: if
|
|
Packit Service |
2e9770 |
z = a + I*b, where the relative error on z is at most 2^(-p), and
|
|
Packit Service |
2e9770 |
EXP(a) = EXP(b) + k, the relative error on b is at most 2^(k-p) */
|
|
Packit Service |
2e9770 |
if ((z_imag || (p > q + 3 + dr && mpfr_can_round (mpc_realref(u), p - q - 3 - dr, MPFR_RNDN, MPFR_RNDZ, pr))) &&
|
|
Packit Service |
2e9770 |
(z_real || (p > q + 3 + di && mpfr_can_round (mpc_imagref(u), p - q - 3 - di, MPFR_RNDN, MPFR_RNDZ, pi))))
|
|
Packit Service |
2e9770 |
break;
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* if Re(u) is not known to be zero, assume it is a normal number, i.e.,
|
|
Packit Service |
2e9770 |
neither zero, Inf or NaN, otherwise we might enter an infinite loop */
|
|
Packit Service |
2e9770 |
MPC_ASSERT (z_imag || mpfr_number_p (mpc_realref(u)));
|
|
Packit Service |
2e9770 |
/* idem for Im(u) */
|
|
Packit Service |
2e9770 |
MPC_ASSERT (z_real || mpfr_number_p (mpc_imagref(u)));
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
if (ret == -2) /* we did not yet call mpc_pow_exact, or it aborted
|
|
Packit Service |
2e9770 |
because intermediate computations had > maxprec bits */
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
/* check exact cases (see algorithms.tex) */
|
|
Packit Service |
2e9770 |
if (y_real)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
maxprec *= 2;
|
|
Packit Service |
2e9770 |
ret = mpc_pow_exact (z, x, mpc_realref(y), rnd, maxprec);
|
|
Packit Service |
2e9770 |
if (ret != -1 && ret != -2)
|
|
Packit Service |
2e9770 |
goto exact;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
p += dr + di + 64;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
else
|
|
Packit Service |
2e9770 |
p += p / 2;
|
|
Packit Service |
2e9770 |
try_again:
|
|
Packit Service |
2e9770 |
mpc_set_prec (t, p);
|
|
Packit Service |
2e9770 |
mpc_set_prec (u, p);
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
if (z_real)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
/* When the result is real (see algorithm.tex for details),
|
|
Packit Service |
2e9770 |
Im(x^y) =
|
|
Packit Service |
2e9770 |
+ sign(imag(y))*0i, if |x| > 1
|
|
Packit Service |
2e9770 |
+ sign(imag(x))*sign(real(y))*0i, if |x| = 1
|
|
Packit Service |
2e9770 |
- sign(imag(y))*0i, if |x| < 1
|
|
Packit Service |
2e9770 |
*/
|
|
Packit Service |
2e9770 |
mpfr_t n;
|
|
Packit Service |
2e9770 |
int inex, cx1;
|
|
Packit Service |
2e9770 |
int sign_zi, sign_rex, sign_imx;
|
|
Packit Service |
2e9770 |
/* cx1 < 0 if |x| < 1
|
|
Packit Service |
2e9770 |
cx1 = 0 if |x| = 1
|
|
Packit Service |
2e9770 |
cx1 > 0 if |x| > 1
|
|
Packit Service |
2e9770 |
*/
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
sign_rex = mpfr_signbit (mpc_realref (x));
|
|
Packit Service |
2e9770 |
sign_imx = mpfr_signbit (mpc_imagref (x));
|
|
Packit Service |
2e9770 |
mpfr_init (n);
|
|
Packit Service |
2e9770 |
inex = mpc_norm (n, x, MPFR_RNDN);
|
|
Packit Service |
2e9770 |
cx1 = mpfr_cmp_ui (n, 1);
|
|
Packit Service |
2e9770 |
if (cx1 == 0 && inex != 0)
|
|
Packit Service |
2e9770 |
cx1 = -inex;
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
sign_zi = (cx1 < 0 && mpfr_signbit (mpc_imagref (y)) == 0)
|
|
Packit Service |
2e9770 |
|| (cx1 == 0 && sign_imx != mpfr_signbit (mpc_realref (y)))
|
|
Packit Service |
2e9770 |
|| (cx1 > 0 && mpfr_signbit (mpc_imagref (y)));
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* copy RE(y) to n since if z==y we will destroy Re(y) below */
|
|
Packit Service |
2e9770 |
mpfr_set_prec (n, mpfr_get_prec (mpc_realref (y)));
|
|
Packit Service |
2e9770 |
mpfr_set (n, mpc_realref (y), MPFR_RNDN);
|
|
Packit Service |
2e9770 |
ret = mpfr_set (mpc_realref(z), mpc_realref(u), MPC_RND_RE(rnd));
|
|
Packit Service |
2e9770 |
if (y_real && (x_real || x_imag))
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
/* FIXME: with y_real we assume Im(y) is really 0, which is the case
|
|
Packit Service |
2e9770 |
for example when y comes from pow_fr, but in case Im(y) is +0 or
|
|
Packit Service |
2e9770 |
-0, we might get different results */
|
|
Packit Service |
2e9770 |
mpfr_set_ui (mpc_imagref (z), 0, MPC_RND_IM (rnd));
|
|
Packit Service |
2e9770 |
fix_sign (z, sign_rex, sign_imx, n);
|
|
Packit Service |
2e9770 |
ret = MPC_INEX(ret, 0); /* imaginary part is exact */
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
else
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
ret = MPC_INEX (ret, mpfr_set_ui (mpc_imagref (z), 0, MPC_RND_IM (rnd)));
|
|
Packit Service |
2e9770 |
/* warning: mpfr_set_ui does not set Im(z) to -0 if Im(rnd) = RNDD */
|
|
Packit Service |
2e9770 |
if (MPC_RND_IM (rnd) == MPFR_RNDD || sign_zi)
|
|
Packit Service |
2e9770 |
mpc_conj (z, z, MPC_RNDNN);
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
mpfr_clear (n);
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
else if (z_imag)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
ret = mpfr_set (mpc_imagref(z), mpc_imagref(u), MPC_RND_IM(rnd));
|
|
Packit Service |
2e9770 |
/* if z is imaginary and y real, then x cannot be real */
|
|
Packit Service |
2e9770 |
if (y_real && x_imag)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
int sign_rex = mpfr_signbit (mpc_realref (x));
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* If z overlaps with y we set Re(z) before checking Re(y) below,
|
|
Packit Service |
2e9770 |
but in that case y=0, which was dealt with above. */
|
|
Packit Service |
2e9770 |
mpfr_set_ui (mpc_realref (z), 0, MPC_RND_RE (rnd));
|
|
Packit Service |
2e9770 |
/* Note: fix_sign only does something when y is an integer,
|
|
Packit Service |
2e9770 |
then necessarily y = 1 or 3 (mod 4), and in that case the
|
|
Packit Service |
2e9770 |
sign of Im(x) is irrelevant. */
|
|
Packit Service |
2e9770 |
fix_sign (z, sign_rex, 0, mpc_realref (y));
|
|
Packit Service |
2e9770 |
ret = MPC_INEX(0, ret);
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
else
|
|
Packit Service |
2e9770 |
ret = MPC_INEX(mpfr_set_ui (mpc_realref(z), 0, MPC_RND_RE(rnd)), ret);
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
else
|
|
Packit Service |
2e9770 |
ret = mpc_set (z, u, rnd);
|
|
Packit Service |
2e9770 |
exact:
|
|
Packit Service |
2e9770 |
mpc_clear (t);
|
|
Packit Service |
2e9770 |
mpc_clear (u);
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* restore underflow and overflow flags from MPFR */
|
|
Packit Service |
2e9770 |
if (saved_underflow)
|
|
Packit Service |
2e9770 |
mpfr_set_underflow ();
|
|
Packit Service |
2e9770 |
if (saved_overflow)
|
|
Packit Service |
2e9770 |
mpfr_set_overflow ();
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
end:
|
|
Packit Service |
2e9770 |
return ret;
|
|
Packit Service |
2e9770 |
}
|