diff options
author | Paul Zimmermann <paul.zimmermann@inria.fr> | 2014-10-14 20:15:19 +0000 |
---|---|---|
committer | Andreas Enge <andreas.enge@inria.fr> | 2014-11-20 15:34:24 +0100 |
commit | 25c5d6c72ae569acf64835ccc05ce82af659dd7b (patch) | |
tree | c810346d4105a423e8a3620d22772056a86e29c6 | |
parent | 0131b49656e0d3d052c835b6787653e9126f9d37 (diff) | |
download | mpc-git-25c5d6c72ae569acf64835ccc05ce82af659dd7b.tar.gz |
fixed bug in mpc_pow
(cf http://lists.gforge.inria.fr/pipermail/mpc-discuss/2014-October/001315.html)
git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/mpc/trunk@1455 211d60ee-9f03-0410-a15a-8952a2c7a4e4
-rw-r--r-- | src/pow.c | 17 |
1 files changed, 8 insertions, 9 deletions
@@ -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 initialisation 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 (); |