Merged changesets r12026-12045 from the trunk (bug fixes and tests).

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

7 files changed, 154 insertions, 8 deletions

diff --git a/src/hypot.c b/src/hypot.c
index 07ceccfb7..246f55d0b 100644
--- a/src/hypot.c
+++ b/src/hypot.c
@@ -98,7 +98,13 @@ mpfr_hypot (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mpfr_rnd_t rnd_mode)
           /* If z > abs(x), then it was already rounded up; otherwise
              z = abs(x), and we need to add one ulp due to y. */
           if (mpfr_abs (z, x, rnd_mode) == 0)
-            mpfr_nexttoinf (z);
+            {
+              mpfr_nexttoinf (z);
+              /* since mpfr_nexttoinf does not set the overflow flag,
+                 we have to check manually for overflow */
+              if (MPFR_UNLIKELY (MPFR_IS_INF (z)))
+                MPFR_SET_OVERFLOW ();
+            }
           MPFR_RET (1);
         }
       else /* MPFR_RNDZ, MPFR_RNDD, MPFR_RNDN */
diff --git a/src/lngamma.c b/src/lngamma.c
index 147ef6bfc..5511fd1dc 100644
--- a/src/lngamma.c
+++ b/src/lngamma.c
@@ -350,7 +350,10 @@ GAMMA_FUNC (mpfr_ptr y, mpfr_srcptr z0, mpfr_rnd_t rnd)
              ulp(u)/2 + (2-z0)*max(1,log(2-z0))*2^(1-w)
              = (1/2 + (2-z0)*max(1,log(2-z0))*2^(1-E(u))) ulp(u) */
           d = (double) MPFR_GET_EXP(s) * 0.694; /* upper bound for log(2-z0) */
-          err_u = MPFR_GET_EXP(s) + __gmpfr_ceil_log2 (d) + 1 - MPFR_GET_EXP(u);
+          if (MPFR_IS_ZERO(u)) /* in that case the error on u is zero */
+            err_u = 0;
+          else
+            err_u = MPFR_GET_EXP(s) + __gmpfr_ceil_log2 (d) + 1 - MPFR_GET_EXP(u);
           err_u = (err_u >= 0) ? err_u + 1 : 0;
           /* now the error on u is bounded by 2^err_u ulps */
 
@@ -397,11 +400,15 @@ GAMMA_FUNC (mpfr_ptr y, mpfr_srcptr z0, mpfr_rnd_t rnd)
             }
           else
             {
-              err_s += 1 - MPFR_GET_EXP(v);
+              /* if v = 0 here, it was 1 before the call to mpfr_log,
+                 thus the error on v was zero */
+              if (!MPFR_IS_ZERO(v))
+                err_s += 1 - MPFR_GET_EXP(v);
               err_s = (err_s >= 0) ? err_s + 1 : 0;
               /* the error on v is bounded by 2^err_s ulps */
               err_u += MPFR_GET_EXP(u); /* absolute error on u */
-              err_s += MPFR_GET_EXP(v); /* absolute error on v */
+              if (!MPFR_IS_ZERO(v)) /* same as above */
+                err_s += MPFR_GET_EXP(v); /* absolute error on v */
               mpfr_sub (s, v, u, MPFR_RNDN);
               /* the total error on s is bounded by ulp(s)/2 + 2^(err_u-w)
                  + 2^(err_s-w) <= ulp(s)/2 + 2^(max(err_u,err_s)+1-w) */
diff --git a/src/sin.c b/src/sin.c
index 142fdf834..ea920c477 100644
--- a/src/sin.c
+++ b/src/sin.c
@@ -148,9 +148,8 @@ mpfr_sin (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
       sign = MPFR_SIGN(xx);
       /* now that the argument is reduced, precision m is enough */
       mpfr_set_prec (c, m);
-      mpfr_cos (c, xx, MPFR_RNDZ);    /* can't be exact */
-      mpfr_nexttoinf (c);           /* now c = cos(x) rounded away */
-      mpfr_mul (c, c, c, MPFR_RNDU); /* away */
+      mpfr_cos (c, xx, MPFR_RNDA);    /* c = cos(x) rounded away */
+      mpfr_mul (c, c, c, MPFR_RNDU);  /* away */
       mpfr_ui_sub (c, 1, c, MPFR_RNDZ);
       mpfr_sqrt (c, c, MPFR_RNDZ);
       if (MPFR_IS_NEG_SIGN(sign))
diff --git a/src/subnormal.c b/src/subnormal.c
index 11a6af3ca..0a135aa9f 100644
--- a/src/subnormal.c
+++ b/src/subnormal.c
@@ -144,7 +144,7 @@ mpfr_subnormalize (mpfr_ptr y, int old_inexact, mpfr_rnd_t rnd)
                 {
                   if (SAME_SIGN (inexact, MPFR_INT_SIGN (y)))
                     mpfr_nexttozero (dest);
-                  else
+                  else  /* subnormal range, thus no overflow */
                     mpfr_nexttoinf (dest);
                   inexact = -inexact;
                 }
diff --git a/tests/thypot.c b/tests/thypot.c
index cb9cefe64..531dd2744 100644
--- a/tests/thypot.c
+++ b/tests/thypot.c
@@ -310,11 +310,40 @@ alltst (void)
     }
 }
 
+/* Test failing with GMP_CHECK_RANDOMIZE=1513841234 (overflow flag not set).
+   The bug was in fact in mpfr_nexttoinf which didn't set the overflow flag. */
+static void
+bug20171221 (void)
+{
+  mpfr_t x, u, y;
+  int inex;
+  mpfr_exp_t emax;
+
+  mpfr_init2 (x, 12);
+  mpfr_init2 (u, 12);
+  mpfr_init2 (y, 11);
+  mpfr_set_str_binary (x, "0.111111111110E0");
+  mpfr_set_str_binary (u, "0.111011110100E-177");
+  emax = mpfr_get_emax ();
+  mpfr_set_emax (0);
+  mpfr_clear_flags ();
+  inex = mpfr_hypot (y, x, u, MPFR_RNDU);
+  MPFR_ASSERTN(mpfr_inf_p (y) && mpfr_sgn (y) > 0);
+  MPFR_ASSERTN(inex > 0);
+  MPFR_ASSERTN(mpfr_inexflag_p ());
+  MPFR_ASSERTN(mpfr_overflow_p ());
+  mpfr_set_emax (emax);
+  mpfr_clear (x);
+  mpfr_clear (u);
+  mpfr_clear (y);
+}
+
 int
 main (int argc, char *argv[])
 {
   tests_start_mpfr ();
 
+  bug20171221 ();
   special ();
 
   test_large ();
diff --git a/tests/tj1.c b/tests/tj1.c
index 534b4486c..c3aab3e78 100644
--- a/tests/tj1.c
+++ b/tests/tj1.c
@@ -27,6 +27,51 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
 #define REDUCE_EMAX 262143 /* otherwise arg. reduction is too expensive */
 #include "tgeneric.c"
 
+/* test for small input, where j1(x) = x/2 - x^3/16 + ... */
+static void
+test_small (void)
+{
+  mpfr_t x, y;
+  int inex, sign;
+  mpfr_exp_t e, emin;
+
+  mpfr_init2 (x, 10);
+  mpfr_init2 (y, 9);
+  emin = mpfr_get_emin ();
+  for (e = -4; e >= -30; e--)
+    {
+      if (e == -30)
+        {
+          e = mpfr_get_emin_min () - 1;
+          mpfr_set_emin (e + 1);
+        }
+      for (sign = -1; sign <= 1; sign += 2)
+        {
+          mpfr_set_si_2exp (x, sign, e, MPFR_RNDN);
+          mpfr_nexttoinf (x);
+          inex = mpfr_j1 (y, x, MPFR_RNDN);
+          if (e >= -29)
+            {
+              /* since |x| is just above 2^e, |j1(x)| is just above 2^(e-1),
+                 thus y should be 2^(e-1) and the inexact flag should be
+                 of opposite sign of x */
+              MPFR_ASSERTN(mpfr_cmp_si_2exp (y, sign, e - 1) == 0);
+              MPFR_ASSERTN(VSIGN (inex) * sign < 0);
+            }
+          else
+            {
+              /* here |y| should be 0.5*2^emin and the inexact flag should
+                 have the sign of x */
+              MPFR_ASSERTN(mpfr_cmp_si_2exp (y, sign, e) == 0);
+              MPFR_ASSERTN(VSIGN (inex) * sign > 0);
+            }
+        }
+    }
+  mpfr_set_emin (emin);
+  mpfr_clear (x);
+  mpfr_clear (y);
+}
+
 int
 main (int argc, char *argv[])
 {
@@ -34,6 +79,8 @@ main (int argc, char *argv[])
 
   tests_start_mpfr ();
 
+  test_small ();
+
   mpfr_init (x);
   mpfr_init (y);
 
diff --git a/tests/tlngamma.c b/tests/tlngamma.c
index 377097807..703da3af7 100644
--- a/tests/tlngamma.c
+++ b/tests/tlngamma.c
@@ -252,11 +252,69 @@ special (void)
   mpfr_clear (y);
 }
 
+/* test failing with GMP_CHECK_RANDOMIZE=1513869588 */
+static void
+bug20171220 (void)
+{
+  mpfr_t x, y, z;
+  int inex;
+
+  mpfr_init2 (x, 15);
+  mpfr_init2 (y, 15);
+  mpfr_init2 (z, 15);
+
+  mpfr_set_str (x, "1.01111e+00", 10, MPFR_RNDN); /* x = 8283/8192 */
+  inex = mpfr_lngamma (y, x, MPFR_RNDN);
+  mpfr_set_str (z, "-0.0063109971733698154140545190234", 10, MPFR_RNDN);
+  MPFR_ASSERTN(mpfr_equal_p (y, z));
+  MPFR_ASSERTN(inex > 0);
+
+  mpfr_set_prec (x, 43);
+  mpfr_set_prec (y, 1);
+  mpfr_set_prec (z, 1);
+  mpfr_set_str (x, "9.8888652212918e-01", 10, MPFR_RNDN);
+  /* lngamma(x) = 0.00651701007222520, should be rounded up to 0.0078125 */
+  inex = mpfr_lngamma (y, x, MPFR_RNDU);
+  mpfr_set_ui_2exp (z, 1, -7, MPFR_RNDN);
+  MPFR_ASSERTN(mpfr_equal_p (y, z));
+  MPFR_ASSERTN(inex > 0);
+
+  mpfr_clear (x);
+  mpfr_clear (y);
+  mpfr_clear (z);
+}
+
+/* taway failing with GMP_CHECK_RANDOMIZE=1513888822 */
+static void
+bug20171220a (void)
+{
+  mpfr_t x, y, z;
+  int inex;
+
+  mpfr_init2 (x, 198);
+  mpfr_init2 (y, 1);
+  mpfr_init2 (z, 1);
+
+  mpfr_set_str (x, "9.999962351340362288585900348170984233205352566408878552154832e-01", 10, MPFR_RNDN);
+  inex = mpfr_lngamma (y, x, MPFR_RNDA);
+  /* lngamma(x) ~ 2.1731512683e0-6 ~ 2^-18.81, should be rounded to 2^-18 */
+  mpfr_set_si_2exp (z, 1, -18, MPFR_RNDN);
+  MPFR_ASSERTN(mpfr_equal_p (y, z));
+  MPFR_ASSERTN(inex > 0);
+
+  mpfr_clear (x);
+  mpfr_clear (y);
+  mpfr_clear (z);
+}
+
 int
 main (void)
 {
   tests_start_mpfr ();
 
+  bug20171220 ();
+  bug20171220a ();
+
   special ();
   test_generic (MPFR_PREC_MIN, 100, 2);