fixed tiny bugs for corner cases

improved test coverage


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

2 files changed, 112 insertions, 8 deletions

diff --git a/exp.c b/exp.c
index 721ff8e76..1b1050e60 100644
--- a/exp.c
+++ b/exp.c
@@ -50,10 +50,11 @@ mpfr_exp (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode)
 	  MPFR_SET_POS(y);
 	  MPFR_RET(0);
 	}
-      else if (MPFR_IS_ZERO(x))
-	return mpfr_set_ui (y, 1, GMP_RNDN);
       else
-	MPFR_ASSERTN(0);
+        {
+          MPFR_ASSERTD(MPFR_IS_ZERO(x));
+          return mpfr_set_ui (y, 1, GMP_RNDN);
+        }
     }
   MPFR_CLEAR_FLAGS(y);
 
@@ -69,14 +70,19 @@ mpfr_exp (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode)
   /* result is 0 when exp(x) < 1/2*2^(__gmpfr_emin), i.e.
      x < (__gmpfr_emin-1) * LOG2 */
   if (d < ((double) __gmpfr_emin - 1.0) * LOG2)
-    return mpfr_set_underflow (y, rnd_mode, 1);
+    {
+      /* warning: mpfr_set_underflow rounds away for RNDN */
+      if (rnd_mode == GMP_RNDN && d < ((double) __gmpfr_emin - 2.0) * LOG2)
+        rnd_mode = GMP_RNDZ;
+      return mpfr_set_underflow (y, rnd_mode, 1);
+    }
 
   /* if x < 2^(-precy), then exp(x) i.e. gives 1 +/- 1 ulp(1) */
   if (expx < -precy)
     {
       int signx = MPFR_SIGN(x);
 
-      if (MPFR_IS_NEG_SIGN(signx) && (rnd_mode == GMP_RNDD || rnd_mode == GMP_RNDZ))
+      if (MPFR_IS_NEG_SIGN(signx) && (rnd_mode == GMP_RNDD))
         {
           MPFR_SET_POS(y);
           mpfr_setmax (y, 0);  /* y = 1 - epsilon */
diff --git a/tests/texp.c b/tests/texp.c
index d574a8c55..a632aadba 100644
--- a/tests/texp.c
+++ b/tests/texp.c
@@ -151,6 +151,17 @@ check_worst_cases (void)
       exit (1);
     }
 
+  mpfr_set_prec (x, 13001);
+  mpfr_set_prec (y, 13001);
+  mpfr_random (x);
+  mpfr_exp (y, x, GMP_RNDN);
+  mpfr_exp_2 (x, x, GMP_RNDN);
+  if (mpfr_cmp (x, y))
+    {
+      printf ("mpfr_exp_2 and mpfr_exp3 differ for prec=13001\n");
+      exit (1);
+    }
+
   mpfr_clear (x);
   mpfr_clear (y);
   return 0;
@@ -200,15 +211,102 @@ static void
 check_special ()
 {
   mpfr_t x, y, z;
+  mp_exp_t emin, emax;
+
+  mpfr_init (x);
+  mpfr_init (y);
+  mpfr_init (z);
+
+  /* check exp(NaN) = NaN */
+  mpfr_set_nan (x);
+  mpfr_exp (y, x, GMP_RNDN);
+  if (!mpfr_nan_p (y))
+    {
+      printf ("Error for exp(NaN)\n");
+      exit (1);
+    }
+
+  /* check exp(+inf) = +inf */
+  mpfr_set_inf (x, 1);
+  mpfr_exp (y, x, GMP_RNDN);
+  if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0)
+    {
+      printf ("Error for exp(+inf)\n");
+      exit (1);
+    }
+
+  /* check exp(-inf) = +0 */
+  mpfr_set_inf (x, -1);
+  mpfr_exp (y, x, GMP_RNDN);
+  if (mpfr_cmp_ui (y, 0) || mpfr_sgn (y) < 0)
+    {
+      printf ("Error for exp(-inf)\n");
+      exit (1);
+    }
+
+  /* check overflow */
+  emax = mpfr_get_emax ();
+  mpfr_set_emax (10);
+  mpfr_set_ui (x, 7, GMP_RNDN);
+  mpfr_exp (y, x, GMP_RNDN);
+  if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0)
+    {
+      printf ("Error for exp(7) for emax=10\n");
+      exit (1);
+    }
+  mpfr_set_emax (emax);
+
+  /* check underflow */
+  emin = mpfr_get_emin ();
+  mpfr_set_emin (-10);
+  mpfr_set_si (x, -9, GMP_RNDN);
+  mpfr_exp (y, x, GMP_RNDN);
+  if (mpfr_cmp_ui (y, 0) || mpfr_sgn (y) < 0)
+    {
+      printf ("Error for exp(-9) for emin=-10\n");
+      printf ("Expected +0\n");
+      printf ("Got      "); mpfr_print_binary (y); puts ("");
+      exit (1);
+    }
+  mpfr_set_emin (emin);
+
+  /* check case EXP(x) < -precy */
+  mpfr_set_prec (y, 2);
+  mpfr_set_str_binary (x, "-0.1E-3");
+  mpfr_exp (y, x, GMP_RNDD);
+  if (mpfr_cmp_ui_2exp (y, 3, -2))
+    {
+      printf ("Error for exp(-1/16), prec=2, RNDD\n");
+      exit (1);
+    }
+  mpfr_exp (y, x, GMP_RNDZ);
+  if (mpfr_cmp_ui (y, 1))
+    {
+      printf ("Error for exp(-1/16), prec=2, RNDZ\n");
+      exit (1);
+    }
+  mpfr_set_str_binary (x, "0.1E-3");
+  mpfr_exp (y, x, GMP_RNDN);
+  if (mpfr_cmp_ui (y, 1))
+    {
+      printf ("Error for exp(1/16), prec=2, RNDN\n");
+      exit (1);
+    }
+  mpfr_exp (y, x, GMP_RNDU);
+  if (mpfr_cmp_ui_2exp (y, 3, -1))
+    {
+      printf ("Error for exp(1/16), prec=2, RNDU\n");
+      exit (1);
+    }
 
   /* bug reported by Franky Backeljauw, 28 Mar 2003 */
-  mpfr_init2 (x, 53);
-  mpfr_init2 (y, 53);
+  mpfr_set_prec (x, 53);
+  mpfr_set_prec (y, 53);
   mpfr_set_str_binary (x, "1.1101011000111101011110000111010010101001101001110111e28");
   mpfr_exp (y, x, GMP_RNDN);
 
   mpfr_set_prec (x, 153);
-  mpfr_init2 (z, 153);
+  mpfr_set_prec (z, 153);
   mpfr_set_str_binary (x, "1.1101011000111101011110000111010010101001101001110111e28");
   mpfr_exp (z, x, GMP_RNDN);
   mpfr_prec_round (z, 53, GMP_RNDN);