summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorzimmerma <zimmerma@211d60ee-9f03-0410-a15a-8952a2c7a4e4>2009-06-22 12:09:01 +0000
committerzimmerma <zimmerma@211d60ee-9f03-0410-a15a-8952a2c7a4e4>2009-06-22 12:09:01 +0000
commit49981fe688dfd8575766a3838ffbae2b47646a57 (patch)
treea411cc5f92a1396dd3d189c10cac4908add494b9
parent27bb54a848f3e28e7d1e731bab9caee88df4976e (diff)
downloadmpc-49981fe688dfd8575766a3838ffbae2b47646a57.tar.gz
[pow.c] save/restore the exponent range in mpc_pow_exact
detect zero, Inf, NaN after mpfr_can_round git-svn-id: svn://scm.gforge.inria.fr/svn/mpc/trunk@613 211d60ee-9f03-0410-a15a-8952a2c7a4e4
-rw-r--r--src/pow.c16
-rw-r--r--tests/pow.dat2
2 files changed, 17 insertions, 1 deletions
diff --git a/src/pow.c b/src/pow.c
index e077b98..6891ad3 100644
--- a/src/pow.c
+++ b/src/pow.c
@@ -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