[src/mul.c] with Mulders', convert RNDF to RNDZ

[tests/tmul_d.c] added a non-regression test


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

2 files changed, 26 insertions, 14 deletions

diff --git a/src/mul.c b/src/mul.c
index c22a49222..a9d78141b 100644
--- a/src/mul.c
+++ b/src/mul.c
@@ -981,22 +981,14 @@ mpfr_mul (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
         /* now the approximation is in tmp[tn-n]...tmp[tn-1] */
         MPFR_ASSERTD (MPFR_LIMB_MSB (tmp[tn-1]) != 0);
 
+        /* for RNDF, we simply use RNDZ, since anyway here we multiply numbers
+           with large precisions, thus the overhead of RNDZ is small */
+        if (rnd_mode == MPFR_RNDF)
+          rnd_mode = MPFR_RNDZ;
+        
         /* if the most significant bit b1 is zero, we have only p-1 correct
            bits */
-        if (rnd_mode == MPFR_RNDF && p + b1 - 1 >= aq
-            && ((aq % GMP_NUMB_BITS) != 0))
-          {
-            /* Warning: while we know the error is bounded by
-               ufp(tmp)/2^(p+b1-1), and we want aq correct bits,
-               simply truncating tmp might give a result which is
-               1 ulp too small, in case tmp ends with 111...111.
-               But if we add one to tmp[tn-n], and then truncate, we
-               cannot exceed what would be obtained by RNDU,
-               since tmp has more bits than a */
-            mpn_add_1 (tmp + tn - n, tmp + tn - n, n, MPFR_LIMB_ONE);
-            /* we can round */
-          }
-        else if (MPFR_UNLIKELY (!mpfr_round_p (tmp, tn, p + b1 - 1,
+        if (MPFR_UNLIKELY (!mpfr_round_p (tmp, tn, p + b1 - 1,
                                                aq + (rnd_mode == MPFR_RNDN))))
           {
             tmp -= k - tn; /* tmp may have changed, FIX IT!!!!! */
diff --git a/tests/tmul_d.c b/tests/tmul_d.c
index d29976146..42787b014 100644
--- a/tests/tmul_d.c
+++ b/tests/tmul_d.c
@@ -70,6 +70,24 @@ check_nans (void)
   mpfr_clear (y);
 }
 
+static void
+bug20171218 (void)
+{
+  mpfr_t a, b;
+  mpfr_init2 (a, 22);
+  mpfr_init2 (b, 1312);
+  mpfr_set_str (b, "1.ffffffffffffffffffc3e0007ffe000700fff8001fff800000000007e0fffe0007fffffff00001f800000003ffffffe1fc00000003fffe7c03ffffffffffffffe000000fffffffff83f0000007ffc03ffffffffc0007ff0000ffcfffffe00000f80000000003c007ffffffff3ff807ffffffffff000000000000001fc000fffe000600000ff0003ffe00fffffffc00001ffc0fffffffffff00000807fe03ffffffc01ffe", 16, MPFR_RNDN);
+  mpfr_mul_d (a, b, 0.5, MPFR_RNDF);
+  /* a should be 1 or nextbelow(1) */
+  if (mpfr_cmp_ui (a, 1) != 0)
+    {
+      mpfr_nextabove (a);
+      MPFR_ASSERTN(mpfr_cmp_ui (a, 1) == 0);
+    }
+  mpfr_clear (a);
+  mpfr_clear (b);
+}
+
 #define TEST_FUNCTION mpfr_mul_d
 #define DOUBLE_ARG2
 #define RAND_FUNCTION(x) mpfr_random2(x, MPFR_LIMB_SIZE (x), 1, RANDS)
@@ -84,6 +102,8 @@ main (void)
 
   tests_start_mpfr ();
 
+  bug20171218 ();
+
   /* check with enough precision */
   mpfr_init2 (x, IEEE_DBL_MANT_DIG);
   mpfr_init2 (y, IEEE_DBL_MANT_DIG);