summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorPaul Zimmermann <paul.zimmermann@inria.fr>2014-10-14 20:15:19 +0000
committerAndreas Enge <andreas.enge@inria.fr>2014-11-20 15:34:24 +0100
commit25c5d6c72ae569acf64835ccc05ce82af659dd7b (patch)
treec810346d4105a423e8a3620d22772056a86e29c6
parent0131b49656e0d3d052c835b6787653e9126f9d37 (diff)
downloadmpc-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.c17
1 files changed, 8 insertions, 9 deletions
diff --git a/src/pow.c b/src/pow.c
index 892f467..2525644 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 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 ();