diff options
-rw-r--r-- | src/pow.c | 16 | ||||
-rw-r--r-- | tests/pow.dat | 2 |
2 files changed, 17 insertions, 1 deletions
@@ -109,7 +109,7 @@ mpc_perfect_square_p (mpz_t a, mpz_t b, mpz_t c, mpz_t d) static int mpc_pow_exact (mpc_ptr z, mpc_srcptr x, mpfr_srcptr y, mpc_rnd_t rnd) { - mp_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; @@ -253,10 +253,18 @@ mpc_pow_exact (mpc_ptr z, mpc_srcptr x, mpfr_srcptr y, mpc_rnd_t rnd) ed *= 2; } + /* 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); @@ -480,6 +488,12 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) 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_RE(u))); + /* idem for Im(u) */ + MPC_ASSERT (z_real || mpfr_number_p (MPC_IM(u))); if (loop == 0) /* first iteration of Ziv's algorithm */ { diff --git a/tests/pow.dat b/tests/pow.dat index e4a8bda..b945121 100644 --- a/tests/pow.dat +++ b/tests/pow.dat @@ -179,3 +179,5 @@ - 0 5 -25 2 +0 2 +0 53 0xCCCCCCCCCCCCDp-54 2 -2 2 -0 N N - - 53 -0x85649E3220691p-63 53 -0x14A25D455A9D0Dp-60 3 5 2 3 2 -3 2 +0 N N + 0 53 0xABCC77118461Dp-74 2 +0 3 5 3 5 2 -8 2 0 N N ++ 0 53 -0x127DB86014739Dp-93 2 +0 2 -1 2 -0 2 1 4 -9 N N ++ + 24 0xC1F98Dp-21 24 0x12FF89p-2 24 -7 24 +0 24 0xCFFFF3p-21 24 +0 N N |