mpc-1.3.1.tar.gz

/* mpc_pow -- Raise a complex number to the power of another complex number. Copyright (C) 2009, 2010, 2011, 2012, 2014, 2015, 2016, 2018, 2020, 2022 INRIA This file is part of GNU MPC. 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 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. 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/ . */ #include <stdio.h> /* for MPC_ASSERT */ #include "mpc-impl.h" /* Return non-zero iff c+i*d is an exact square (a+i*b)^2, with a, b both of the form m*2^e with m, e integers. If so, returns in a+i*b the corresponding square root, with a >= 0. The variables a, b must not overlap with c, d. We have c = a^2 - b^2 and d = 2*a*b. If one of a, b is exact, then both are (see algorithms.tex). Case 1: a <> 0 and b <> 0. Let a = m*2^e and b = n*2^f with m, e, n, f integers, m and n odd (we will treat apart the case a = 0 or b = 0). Then 2*a*b = m*n*2^(e+f+1), thus necessarily e+f >= -1. Assume e < 0, then f >= 0, then a^2 - b^2 = m^2*2^(2e) - n^2*2^(2f) cannot be an integer, since n^2*2^(2f) is an integer, and m^2*2^(2e) is not. Similarly when f < 0 (and thus e >= 0). Thus we have e, f >= 0, and a, b are both integers. Let A = 2a^2, then eliminating b between c = a^2 - b^2 and d = 2*a*b gives A^2 - 2c*A - d^2 = 0, which has solutions c +/- sqrt(c^2+d^2). We thus need c^2+d^2 to be a square, and c + sqrt(c^2+d^2) --- the solution we are interested in --- to be two times a square. Then b = d/(2a) is necessarily an integer. Case 2: a = 0. Then d is necessarily zero, thus it suffices to check whether c = -b^2, i.e., if -c is a square. Case 3: b = 0. Then d is necessarily zero, thus it suffices to check whether c = a^2, i.e., if c is a square. */ static int mpc_perfect_square_p (mpz_t a, mpz_t b, mpz_t c, mpz_t d) { if (mpz_cmp_ui (d, 0) == 0) /* case a = 0 or b = 0 */ { /* necessarily c < 0 here, since we have already considered the case where x is real non-negative and y is real */ MPC_ASSERT (mpz_cmp_ui (c, 0) < 0); mpz_neg (b, c); if (mpz_perfect_square_p (b)) /* case 2 above */ { mpz_sqrt (b, b); mpz_set_ui (a, 0); return 1; /* c + i*d = (0 + i*b)^2 */ } } else /* both a and b are non-zero */ { if (mpz_divisible_2exp_p (d, 1) == 0) return 0; /* d must be even */ mpz_mul (a, c, c); mpz_addmul (a, d, d); /* c^2 + d^2 */ if (mpz_perfect_square_p (a)) { mpz_sqrt (a, a); mpz_add (a, c, a); /* c + sqrt(c^2+d^2) */ if (mpz_divisible_2exp_p (a, 1)) { mpz_tdiv_q_2exp (a, a, 1); if (mpz_perfect_square_p (a)) { mpz_sqrt (a, a); mpz_tdiv_q_2exp (b, d, 1); /* d/2 */ mpz_divexact (b, b, a); /* d/(2a) */ return 1; } } } } 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), MPFR_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), MPFR_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]. If x^y is not exactly representable, return -1. If intermediate computations lead to numbers of more than maxprec bits, then abort and return -2 (in that case, to avoid loops, mpc_pow_exact 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) { mpfr_exp_t ec, ed, ey; 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; int inex_im; 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, MPFR_RNDN); } mpz_init (my); mpz_init (a); mpz_init (b); mpz_init (c); mpz_init (d); mpz_init (u); 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 with my odd */ if (x_imag) { mpz_set_ui (c, 0); ec = 0; } else ec = mpfr_get_z_exp (c, mpc_realref(x)); if (mpfr_zero_p (mpc_imagref(x))) { mpz_set_ui (d, 0); ed = ec; } else { ed = mpfr_get_z_exp (d, mpc_imagref(x)); if (x_imag) 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) goto end; } else if (ed < ec) { mpz_mul_2exp (c, c, (unsigned long int) (ec - ed)); if ((mpfr_prec_t) mpz_sizeinbase (c, 2) > maxprec) goto end; ec = ed; } /* now ec=ed and x = (c + I * d) * 2^ec */ /* divide by two if possible */ if (mpz_cmp_ui (c, 0) == 0) { t = mpz_scan1 (d, 0); mpz_tdiv_q_2exp (d, d, t); ec += (mpfr_exp_t) 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; } else /* neither c nor d is zero */ { unsigned long v; t = mpz_scan1 (c, 0); v = mpz_scan1 (d, 0);