pow_ui.c: added code for a second binary trial in certain cases

pow_ui.dat: example needing a second pass


git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/mpc/trunk@756 211d60ee-9f03-0410-a15a-8952a2c7a4e4

3 files changed, 69 insertions, 47 deletions

diff --git a/src/pow_ui.c b/src/pow_ui.c
index 4337294..c1a3cb1 100644
--- a/src/pow_ui.c
+++ b/src/pow_ui.c
@@ -27,7 +27,6 @@ mpc_pow_ui_naive (mpc_ptr z, mpc_srcptr x, unsigned long y, mpc_rnd_t rnd)
 {
    int inex;
    mpc_t t;
-
    mpc_init3 (t, sizeof (unsigned long) * CHAR_BIT, MPFR_PREC_MIN);
    mpc_set_ui (t, y, MPC_RNDNN);   /* exact */
    inex = mpc_pow (z, x, t, rnd);
@@ -42,10 +41,10 @@ mpc_pow_ui (mpc_ptr z, mpc_srcptr x, unsigned long y, mpc_rnd_t rnd)
 {
    int inex;
    mpc_t t, x3;
-   mp_prec_t p, er, ei;
+   mp_prec_t p;
    unsigned long l, l0, u;
-   mp_exp_t diff;
    int has3; /* non-zero if y has '11' in its binary representation */
+   int loop, done;
 
    if (!mpc_fin_p (x) || mpfr_zero_p (MPC_RE (x)) || mpfr_zero_p (MPC_IM(x))
        || y == 0)
@@ -76,55 +75,76 @@ mpc_pow_ui (mpc_ptr z, mpc_srcptr x, unsigned long y, mpc_rnd_t rnd)
    if (has3)
       mpc_init2 (x3, p);
 
-   mpc_sqr (t, x, MPC_RNDNN);
-   if (has3) {
-      mpc_mul (x3, t, x, MPC_RNDNN);
-      if ((y >> l) & 1) /* y starts with 11... */
-      mpc_set (t, x3, MPC_RNDNN);
-   }
-   while (l-- > 0) {
-      mpc_sqr (t, t, MPC_RNDNN);
-      if ((y >> l) & 1) {
-         if ((l > 0) && ((y >> (l-1)) & 1)) /* implies has3 <> 0 */ {
-            l --;
-            mpc_sqr (t, t, MPC_RNDNN);
-            mpc_mul (t, t, x3, MPC_RNDNN);
+   loop = 0;
+   done = 0;
+   while (!done) {
+      loop++;
+
+      mpc_sqr (t, x, MPC_RNDNN);
+      if (has3) {
+         mpc_mul (x3, t, x, MPC_RNDNN);
+         if ((y >> l) & 1) /* y starts with 11... */
+            mpc_set (t, x3, MPC_RNDNN);
+      }
+      while (l-- > 0) {
+         mpc_sqr (t, t, MPC_RNDNN);
+         if ((y >> l) & 1) {
+            if ((l > 0) && ((y >> (l-1)) & 1)) /* implies has3 <> 0 */ {
+               l--;
+               mpc_sqr (t, t, MPC_RNDNN);
+               mpc_mul (t, t, x3, MPC_RNDNN);
+            }
+            else
+               mpc_mul (t, t, x, MPC_RNDNN);
          }
-         else
-            mpc_mul (t, t, x, MPC_RNDNN);
       }
-   }
-   if (has3)
-      mpc_clear (x3);
 
-   /* the absolute error on the real and imaginary parts is bounded
-      by |x|^y (|1+2^{-p}|^{y-1}-1) [see algorithms.tex].
-      For em <= 1, (1+e)^m - 1 <= 2em since
-      (1+e)^m - 1 = exp(m*log(1+e))-1 <= exp(em)-1 <= 2em for em <= 1.
-      We apply this for e=2^{-p} and m=y-1, thus the absolute error is
-      bounded by |x|^y 2^{1-p} (y-1) < 2^{l0+1-p} */
-   if (mpfr_zero_p (MPC_RE(t)) || mpfr_zero_p (MPC_IM(t)))
-      inex = mpc_pow_ui_naive (z, x, y, rnd);
-         /* since mpfr_get_exp() is not defined for zero */
-   else {
-      diff = mpfr_get_exp (MPC_RE(t)) - mpfr_get_exp (MPC_IM(t));
-      /* the factor on the real part is 2+2^(-diff+2) <= 4 for diff >= 1
-         and <= 2^(-diff+3) for diff <= 0 */
-      er = (diff >= 1) ? l0 + 3 : l0 + (-diff) + 4;
-      /* the factor on the imaginary part is 2+2^(diff+2) <= 4 for diff <= -1
-         and <= 2^(diff+3) for diff >= 0 */
-      ei = (diff <= -1) ? l0 + 3 : l0 + diff + 4;
-      if (mpfr_can_round (MPC_RE(t), p - er, GMP_RNDZ, GMP_RNDZ,
-                           MPFR_PREC(MPC_RE(z)) + (MPC_RND_RE(rnd) == GMP_RNDN))
-            && mpfr_can_round (MPC_IM(t), p - ei, GMP_RNDZ, GMP_RNDZ,
-                           MPFR_PREC(MPC_IM(z)) + (MPC_RND_IM(rnd) == GMP_RNDN)))
-         inex = mpc_set (z, t, rnd);
-      else
+      /* the absolute error on the real and imaginary parts is bounded
+         by |x|^y (|1+2^{-p}|^{y-1}-1) [see algorithms.tex].
+         For em <= 1, (1+e)^m - 1 <= 2em since
+         (1+e)^m - 1 = exp(m*log(1+e))-1 <= exp(em)-1 <= 2em for em <= 1.
+         We apply this for e=2^{-p} and m=y-1, thus the absolute error is
+         bounded by |x|^y 2^{1-p} (y-1) < 2^{l0+1-p} */
+      if (mpfr_zero_p (MPC_RE(t)) || mpfr_zero_p (MPC_IM(t))) {
          inex = mpc_pow_ui_naive (z, x, y, rnd);
+            /* since mpfr_get_exp() is not defined for zero */
+         done = 1;
+      }
+      else {
+         mp_exp_t diff;
+         mp_prec_t er, ei;
+         diff = mpfr_get_exp (MPC_RE(t)) - mpfr_get_exp (MPC_IM(t));
+         /* the factor on the real part is 2+2^(-diff+2) <= 4 for diff >= 1
+            and <= 2^(-diff+3) for diff <= 0 */
+         er = (diff >= 1) ? l0 + 3 : l0 + (-diff) + 4;
+         /* the factor on the imaginary part is 2+2^(diff+2) <= 4 for diff <= -1
+            and <= 2^(diff+3) for diff >= 0 */
+         ei = (diff <= -1) ? l0 + 3 : l0 + diff + 4;
+         if (mpfr_can_round (MPC_RE(t), p - er, GMP_RNDZ, GMP_RNDZ,
+                              MPFR_PREC(MPC_RE(z)) + (MPC_RND_RE(rnd) == GMP_RNDN))
+               && mpfr_can_round (MPC_IM(t), p - ei, GMP_RNDZ, GMP_RNDZ,
+                              MPFR_PREC(MPC_IM(z)) + (MPC_RND_IM(rnd) == GMP_RNDN))) {
+            inex = mpc_set (z, t, rnd);
+            done = 1;
+         }
+         else if (loop == 1 && SAFE_ABS(mp_prec_t, diff) < MPC_MAX_PREC(x)) {
+            /* common case, make a second trial at higher precision */
+            p += MPC_MAX_PREC(x);
+            mpc_set_prec (t, p);
+            if (has3)
+               mpc_set_prec (x3, p);
+            l = l0 - 2;
+         }
+         else {
+            inex = mpc_pow_ui_naive (z, x, y, rnd);
+            done = 1;
+         }
+      }
    }
 
    mpc_clear (t);
+   if (has3)
+      mpc_clear (x3);
 
    return inex;
 }
-
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 29f81a5..6907381 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -4,10 +4,10 @@ AM_CPPFLAGS = -I$(top_srcdir)/src
 LDADD = libmpc-tests.la $(top_builddir)/src/libmpc.la
 LOADLIBES=$(DEFS) -I$(top_srcdir)/src $(CPPFLAGS) $(CFLAGS) $(top_builddir)/src/.libs/libmpc.a $(LIBS)
 
-check_PROGRAMS = tpow_ui tabs tacos tacosh tadd tadd_fr tadd_ui targ tasin tasinh     \
+check_PROGRAMS = tabs tacos tacosh tadd tadd_fr tadd_ui targ tasin tasinh     \
 tatan tatanh tconj tcos tcosh tdiv tdiv_2exp tdiv_fr tdiv_ui texp tfr_div     \
 tfr_sub timag tio_str tlog tmul tmul_2exp tmul_fr tmul_i tmul_si \
-tmul_ui tneg tnorm tpow tpow_ld tpow_d tpow_fr tpow_si tpow_z tprec   \
+tmul_ui tneg tnorm tpow tpow_ld tpow_d tpow_fr tpow_si tpow_ui tpow_z tprec   \
 tproj treal treimref tset tsin tsinh tsqr tsqrt tstrtoc tsub tsub_fr tsub_ui  \
 ttan ttanh tui_div tui_ui_sub tget_version
 
diff --git a/tests/pow_ui.dat b/tests/pow_ui.dat
index be8c369..744bd32 100644
--- a/tests/pow_ui.dat
+++ b/tests/pow_ui.dat
@@ -91,3 +91,5 @@
 ? ? 53 +inf 53 +inf    53 1e100000000  53 1e100000000   1000000000 N N
 # underflow
 ? ? 53    0 53    0    53 1e-100000000 53 1e-100000000  1000000000 N N
+# cannot round after one loop
+? ? 420  -0x1.c3fb41a71665f9a144927e70cbc2dc899e9e30880c0b5aa924ad8a538b4cd06e503f38bdbb7cfcfded29f7504fe0c91ecd4230984@-187 420 -0xc.82a09ac98133eb05b2643c98eb1c8e1a1609e75f682b14098176abd6c8b4b3c6c72dadaf8929f9bd87f8c78d03361bacb9fb13140@-292  420 0x1.cf13ce58adc4e639fd1c3063ffc9291433647999951bc04ba6797ec4de0335336ad0a28df18573d3b6322ebab662c08eadaed4a8e@-8 420 0x3.cf71d602ca6f754ebd6af522154f3ee1c46da0a52deb1f60016fca4b1e0b4b447b752169e837bb1866aa3734850cd158a7e3ca33c@-9  24 N N