[src/erf.c] In the computation of an error bound, replaced some

  double's (which could overflow) by mpfr_t to fix bug reported
  by Naoki Shibata:
  https://sympa.inria.fr/sympa/arc/mpfr/2018-07/msg00028.html
[tests/terf.c] Added a testcase for this bug. Increased the number
  of generic tests in order to reproduce the bug there too with the
  default seed.
(merged changesets r12946-12949 from the trunk)

git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/branches/4.0@12950 280ebfd0-de03-0410-8827-d642c229c3f4

2 files changed, 34 insertions, 10 deletions

diff --git a/src/erf.c b/src/erf.c
index 28c4c585d..c547152ae 100644
--- a/src/erf.c
+++ b/src/erf.c
@@ -186,10 +186,8 @@ mpfr_erf_0 (mpfr_ptr res, mpfr_srcptr x, double xf2, mpfr_rnd_t rnd_mode)
 {
   mpfr_prec_t n, m;
   mpfr_exp_t nuk, sigmak;
-  double tauk;
   mpfr_t y, s, t, u;
   unsigned int k;
-  int log2tauk;
   int inex;
   MPFR_GROUP_DECL (group);
   MPFR_ZIV_DECL (loop);
@@ -204,10 +202,14 @@ mpfr_erf_0 (mpfr_ptr res, mpfr_srcptr x, double xf2, mpfr_rnd_t rnd_mode)
   MPFR_ZIV_INIT (loop, m);
   for (;;)
     {
+      mpfr_t tauk;
+      mpfr_exp_t log2tauk;
+
       mpfr_mul (y, x, x, MPFR_RNDU); /* err <= 1 ulp */
       mpfr_set_ui (s, 1, MPFR_RNDN);
       mpfr_set_ui (t, 1, MPFR_RNDN);
-      tauk = 0.0;
+      mpfr_init2 (tauk, 53);
+      mpfr_set_ui (tauk, 0, MPFR_RNDU);
 
       for (k = 1; ; k++)
         {
@@ -222,12 +224,13 @@ mpfr_erf_0 (mpfr_ptr res, mpfr_srcptr x, double xf2, mpfr_rnd_t rnd_mode)
           sigmak -= MPFR_GET_EXP(s);
           nuk = MPFR_GET_EXP(u) - MPFR_GET_EXP(s);
 
-          if ((nuk < - (mpfr_exp_t) m) && ((double) k >= xf2))
+          if ((nuk < - (mpfr_exp_t) m) && (k >= xf2))
             break;
 
           /* tauk <- 1/2 + tauk * 2^sigmak + (1+8k)*2^nuk */
-          tauk = 0.5 + mul_2exp (tauk, sigmak)
-            + mul_2exp (1.0 + 8.0 * (double) k, nuk);
+          mpfr_mul_2si (tauk, tauk, sigmak, MPFR_RNDU);
+          mpfr_add_d (tauk, tauk, 0.5 + mul_2exp (1.0 + 8.0 * (double) k, nuk),
+                      MPFR_RNDU);
         }
 
       mpfr_mul (s, x, s, MPFR_RNDU);
@@ -236,8 +239,14 @@ mpfr_erf_0 (mpfr_ptr res, mpfr_srcptr x, double xf2, mpfr_rnd_t rnd_mode)
       mpfr_const_pi (t, MPFR_RNDZ);
       mpfr_sqrt (t, t, MPFR_RNDZ);
       mpfr_div (s, s, t, MPFR_RNDN);
-      tauk = 4.0 * tauk + 11.0; /* final ulp-error on s */
-      log2tauk = __gmpfr_ceil_log2 (tauk);
+      /* tauk = 4 * tauk + 11: final ulp-error on s */
+      mpfr_mul_2ui (tauk, tauk, 2, MPFR_RNDU);
+      mpfr_add_ui (tauk, tauk, 11, MPFR_RNDU);
+      /* In practice, one should not get an infinity, as it would require
+         a huge precision and a lot of time, but just in case... */
+      MPFR_ASSERTN (!MPFR_IS_INF (tauk));
+      log2tauk = MPFR_GET_EXP (tauk);
+      mpfr_clear (tauk);
 
       if (MPFR_LIKELY (MPFR_CAN_ROUND (s, m - log2tauk, n, rnd_mode)))
         break;
diff --git a/tests/terf.c b/tests/terf.c
index 5e95cee88..c01b5d546 100644
--- a/tests/terf.c
+++ b/tests/terf.c
@@ -636,6 +636,20 @@ reduced_expo_range (void)
   mpfr_set_emax (emax);
 }
 
+/* Similar to a bug reported by Naoki Shibata:
+   https://sympa.inria.fr/sympa/arc/mpfr/2018-07/msg00028.html
+*/
+static void
+bug20180723 (void)
+{
+  mpfr_t x;
+
+  mpfr_init2 (x, 256);
+  mpfr_set_ui (x, 28, MPFR_RNDN);
+  mpfr_erfc (x, x, MPFR_RNDN);
+  mpfr_clear (x);
+}
+
 int
 main (int argc, char *argv[])
 {
@@ -646,9 +660,10 @@ main (int argc, char *argv[])
   large_arg ();
   test_erfc ();
   reduced_expo_range ();
+  bug20180723 ();
 
-  test_generic_erf (MPFR_PREC_MIN, 100, 15);
-  test_generic_erfc (MPFR_PREC_MIN, 100, 15);
+  test_generic_erf (MPFR_PREC_MIN, 300, 150);
+  test_generic_erfc (MPFR_PREC_MIN, 300, 150);
 
   data_check ("data/erf",  mpfr_erf,  "mpfr_erf");
   data_check ("data/erfc", mpfr_erfc, "mpfr_erfc");