From 43bd35ec940648a1424c6f48f0dafd54faafec5a Mon Sep 17 00:00:00 2001 From: zimmerma Date: Tue, 14 Oct 2014 20:15:19 +0000 Subject: fixed bug in mpc_pow (cf http://lists.gforge.inria.fr/pipermail/mpc-discuss/2014-October/001315.html) git-svn-id: svn://scm.gforge.inria.fr/svn/mpc/trunk@1455 211d60ee-9f03-0410-a15a-8952a2c7a4e4 --- src/pow.c | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/src/pow.c b/src/pow.c index 5a2b819..9047a60 100644 --- a/src/pow.c +++ b/src/pow.c @@ -645,7 +645,7 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) pi = mpfr_get_prec (mpc_imagref(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 */ + Ziv's strategy; probably wrong now since q is not computed */ if (p < 64) p = 64; mpc_init2 (u, p); @@ -658,18 +658,17 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) { int ret_exp; mpfr_exp_t dr, di; - mpfr_prec_t q=0; - /* to avoid warning message, real initialization below */ + mpfr_prec_t q; 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)); - } + /* Compute q such that |Re (y log x)|, |Im (y log x)| < 2^q. + We recompute it at each loop since we might get different + bounds if the precision is not enough. */ + 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 (); -- cgit v1.2.1