patch from Patrick to solve efficiency problem when one operand is sparse

(e.g. from ui_pow_ui)


git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@3739 280ebfd0-de03-0410-8827-d642c229c3f4

1 files changed, 49 insertions, 17 deletions

diff --git a/mul.c b/mul.c
index 86c026309..236cc3726 100644
--- a/mul.c
+++ b/mul.c
@@ -175,19 +175,15 @@ mpfr_mul (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode)
   if (mpfr_cmp (ta, a) || inexact1*inexact2 < 0
       || (inexact1*inexact2 == 0 && (inexact1|inexact2) != 0))
     {
-      fprintf (stderr, "mpfr_mul return different values for %s\n"
-               "Prec_a = %lu, Prec_b = %lu, Prec_c = %lu\nB = ",
-               mpfr_print_rnd_mode (rnd_mode),
-               MPFR_PREC (a), MPFR_PREC (b), MPFR_PREC (c));
-      mpfr_out_str (stderr, 16, 0, tb, GMP_RNDN);
-      fprintf (stderr, "\nC = ");
-      mpfr_out_str (stderr, 16, 0, tc, GMP_RNDN);
-      fprintf (stderr, "\nOldMul: ");
-      mpfr_out_str (stderr, 16, 0, ta, GMP_RNDN);
-      fprintf (stderr, "\nNewMul: ");
-      mpfr_out_str (stderr, 16, 0, a, GMP_RNDN);
-      fprintf (stderr, "\nNewInexact = %d | OldInexact = %d\n",
-               inexact1, inexact2);
+      printf ("mpfr_mul return different values for %s\n"
+              "Prec_a= %lu, Prec_b= %lu, Prec_c= %lu\nB= ",
+              mpfr_print_rnd_mode (rnd_mode),
+              MPFR_PREC (a), MPFR_PREC (b), MPFR_PREC (c));
+      mpfr_out_str (stdout, 16, 0, tb, GMP_RNDN);
+      printf("\nC="); mpfr_out_str (stdout, 16, 0, tc, GMP_RNDN);
+      printf("\nOldMul: "); mpfr_out_str (stdout, 16, 0, ta, GMP_RNDN);
+      printf("\nNewMul: "); mpfr_out_str (stdout, 16, 0, a, GMP_RNDN);
+      printf("\nNewInexact = %d | OldInexact = %d\n", inexact1, inexact2);
       MPFR_ASSERTN(0);
     }
 
@@ -341,7 +337,6 @@ mpfr_mul (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode)
           /* Sum those two partial products */
           add_ssaaaa (tmp[2], tmp[1], tmp[2], tmp[1], t1, t2);
           tmp[3] += (tmp[2] < t1);
-
           b1 = tmp[3];
         }
       b1 >>= (BITS_PER_MP_LIMB - 1);
@@ -357,15 +352,52 @@ mpfr_mul (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode)
       mp_size_t n;
       mp_prec_t p;
 
+      /* First check if we can reduce the precision of b or c:
+         exact values are a nightmare for the short product trick */
+      bp = MPFR_MANT (b);
+      cp = MPFR_MANT (c);
+      MPFR_ASSERTN (MPFR_MUL_THRESHOLD >= 1);
+      if (MPFR_UNLIKELY ((bp[0] == 0 && bp[1] == 0)
+                         || (cp[0] == 0 && cp[1] == 0)))
+        {
+          mpfr_t b_tmp, c_tmp;
+          /* Check for b */
+          while (*bp == 0)
+            {
+              bp ++;
+              bn --;
+              MPFR_ASSERTD (bn > 0);
+            } /* This must end since the MSL is != 0 */
+          /* Check for c too */
+          while (*cp == 0)
+            {
+              cp ++;
+              cn --;
+              MPFR_ASSERTD (cn > 0);
+            } /* This must end since the MSL is != 0 */
+          /* It is not the faster way, but it is safer */
+          MPFR_SET_SAME_SIGN (b_tmp, b);
+          MPFR_SET_EXP (b_tmp, MPFR_GET_EXP (b));
+          MPFR_PREC (b_tmp) = bn * BITS_PER_MP_LIMB;
+          MPFR_MANT (b_tmp) = bp;
+          
+          MPFR_SET_SAME_SIGN (c_tmp, c);
+          MPFR_SET_EXP (c_tmp, MPFR_GET_EXP (c));
+          MPFR_PREC (c_tmp) = cn * BITS_PER_MP_LIMB;
+          MPFR_MANT (c_tmp) = cp;
+          /* Recall again mpfr_mul with the fixed arguments */
+          return mpfr_mul (a, b_tmp, c_tmp, rnd_mode);
+        }
+
       /* Compute estimated precision of mulhigh.
          We could use `+ (n < cn) + (n < bn)' instead of `+ 2',
          but does it worth it? */
       n = MPFR_LIMB_SIZE (a) + 1;
       n = MIN (n, cn);
       MPFR_ASSERTD (n >= 1 && 2*n <= k && n <= cn && n <= bn);
-      p = n*BITS_PER_MP_LIMB - MPFR_INT_CEIL_LOG2 (n + 2);
-      bp = MPFR_MANT (b) + bn - n;
-      cp = MPFR_MANT (c) + cn - n;
+      p = n * BITS_PER_MP_LIMB - MPFR_INT_CEIL_LOG2 (n + 2);
+      bp += bn - n;
+      cp += cn - n;
       /* Check if MulHigh can produce a roundable result.
          We may lost 1 bit due to RNDN, 1 due to final shift. */
       if (MPFR_UNLIKELY (MPFR_PREC (a) > p - 5))