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

4 files changed, 20 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;
                 }