Clean up code.

Fix bug with Exponent range.
Add ZivLoop controller.


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

1 files changed, 50 insertions, 60 deletions

diff --git a/erf.c b/erf.c
index 3a7bd042d..dea027ecb 100644
--- a/erf.c
+++ b/erf.c
@@ -27,30 +27,27 @@ MA 02111-1307, USA. */
 
 #define EXP1 2.71828182845904523536 /* exp(1) */
 
-static int mpfr_erf_0 (mpfr_ptr, mpfr_srcptr, mp_rnd_t);
+static int mpfr_erf_0 (mpfr_ptr, mpfr_srcptr, double, mp_rnd_t);
 
 int
 mpfr_erf (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode)
 {
   double xf;
-  int sign_x;
-  mp_rnd_t rnd2;
-  double n = (double) MPFR_PREC(y);
   int inex;
+  MPFR_SAVE_EXPO_DECL (expo);
 
-  sign_x = MPFR_SIGN (x);
-  if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(x) ))
+  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
     {
-      if (MPFR_IS_NAN(x))
+      if (MPFR_IS_NAN (x))
 	{
-	  MPFR_SET_NAN(y);
+	  MPFR_SET_NAN (y);
 	  MPFR_RET_NAN;
 	}
-      else if (MPFR_IS_INF(x)) /* erf(+inf) = +1, erf(-inf) = -1 */
-	return mpfr_set_si (y, MPFR_FROM_SIGN_TO_INT(sign_x), GMP_RNDN);
+      else if (MPFR_IS_INF (x)) /* erf(+inf) = +1, erf(-inf) = -1 */
+	return mpfr_set_si (y, MPFR_INT_SIGN (x), GMP_RNDN);
       else /* erf(+0) = +0, erf(-0) = -0 */
         {
-          MPFR_ASSERTD(MPFR_IS_ZERO(x));
+          MPFR_ASSERTD (MPFR_IS_ZERO (x));
           return mpfr_set (y, x, GMP_RNDN); /* should keep the sign of x */
         }
     }
@@ -60,39 +57,31 @@ mpfr_erf (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode)
   xf = mpfr_get_d (x, GMP_RNDN);
   xf = xf * xf; /* xf ~ x^2 */
 
-  rnd2 = MPFR_IS_POS_SIGN(sign_x) ? rnd_mode : MPFR_INVERT_RND(rnd_mode);
-
   /* use expansion at x=0 when e*x^2 <= n (target precision)
      otherwise use asymptotic expansion */
+  MPFR_SAVE_EXPO_MARK (expo);
 
-  if (xf > n * LOG2) /* |erf x| = 1 or 1- */
+  if (MPFR_UNLIKELY (xf > LOG2 * ((double) MPFR_PREC (y)))) 
+    /* |erf x| = 1 or 1- */
     {
+      mp_rnd_t rnd2 = MPFR_IS_POS (x) ? rnd_mode : MPFR_INVERT_RND(rnd_mode);
       if (rnd2 == GMP_RNDN || rnd2 == GMP_RNDU)
         {
-          if (MPFR_IS_POS_SIGN(sign_x))
-            {
-              mpfr_set_ui (y, 1, rnd2);
-              inex = 1;
-            }
-          else
-            {
-              mpfr_set_si (y, -1, rnd2);
-              inex = -1;
-            }
+	  inex = MPFR_INT_SIGN (x);
+	  mpfr_set_si (y, inex, rnd2);
         }
       else /* round to zero */
         {
+          inex = -MPFR_INT_SIGN (x);
           mpfr_setmax (y, 0); /* warning: setmax keeps the old sign of y */
-          MPFR_SET_SAME_SIGN(y, x);
-          inex = MPFR_IS_POS_SIGN(sign_x) ? -1 : 1;
+          MPFR_SET_SAME_SIGN (y, x);
         }
     }
   else  /* use Taylor */
-    {
-      inex = mpfr_erf_0 (y, x, rnd_mode);
-    }
+    inex = mpfr_erf_0 (y, x, xf, rnd_mode);
 
-  return inex;
+  MPFR_SAVE_EXPO_FREE (expo);
+  return mpfr_check_range (y, inex, rnd_mode);
 }
 
 /* return x*2^e */
@@ -121,41 +110,33 @@ mul_2exp (double x, mp_exp_t e)
    Assumes also that e*x^2 <= n (target precision).
  */
 static int
-mpfr_erf_0 (mpfr_ptr res, mpfr_srcptr x, mp_rnd_t rnd_mode)
+mpfr_erf_0 (mpfr_ptr res, mpfr_srcptr x, double xf2, mp_rnd_t rnd_mode)
 {
   mp_prec_t n, m;
   mp_exp_t nuk, sigmak;
-  double xf, tauk;
+  double tauk;
   mpfr_t y, s, t, u;
   unsigned int k;
-  long log2tauk;
-  int ok;
+  int log2tauk;
   int inex;
+  MPFR_ZIV_DECL (loop);
 
-  n = MPFR_PREC(res); /* target precision */
-  xf = mpfr_get_d (x, GMP_RNDN);
+  n = MPFR_PREC (res); /* target precision */
 
   /* initial working precision */
-  m = n + (mp_prec_t) (xf * xf / LOG2) + 8;
+  m = n + (mp_prec_t) (xf2 / LOG2) + 8 + MPFR_INT_CEIL_LOG2 (n);
 
-  mpfr_init2 (y, 2);
-  mpfr_init2 (s, 2);
-  mpfr_init2 (t, 2);
-  mpfr_init2 (u, 2);
+  mpfr_init2 (y, m);
+  mpfr_init2 (s, m);
+  mpfr_init2 (t, m);
+  mpfr_init2 (u, m);
 
-  do
+  MPFR_ZIV_INIT (loop, m);
+  for (;;)
     {
-
-      m += MPFR_INT_CEIL_LOG2 (n);
-
-      mpfr_set_prec (y, m);
-      mpfr_set_prec (s, m);
-      mpfr_set_prec (t, m);
-      mpfr_set_prec (u, m);
-
       mpfr_mul (y, x, x, GMP_RNDU); /* err <= 1 ulp */
-      mpfr_set_ui (s, 1, GMP_RNDN); /* exact */
-      mpfr_set_ui (t, 1, GMP_RNDN); /* exact */
+      mpfr_set_ui (s, 1, GMP_RNDN);
+      mpfr_set_ui (t, 1, GMP_RNDN);
       tauk = 0.0;
 
       for (k = 1; ; k++)
@@ -163,15 +144,15 @@ mpfr_erf_0 (mpfr_ptr res, mpfr_srcptr x, mp_rnd_t rnd_mode)
           mpfr_mul (t, y, t, GMP_RNDU);
           mpfr_div_ui (t, t, k, GMP_RNDU);
           mpfr_div_ui (u, t, 2 * k + 1, GMP_RNDU);
-          sigmak = MPFR_EXP(s);
+          sigmak = MPFR_GET_EXP (s);
           if (k % 2)
             mpfr_sub (s, s, u, GMP_RNDN);
           else
             mpfr_add (s, s, u, GMP_RNDN);
-          sigmak -= MPFR_EXP(s);
-          nuk = MPFR_EXP(u) - MPFR_EXP(s);
+          sigmak -= MPFR_GET_EXP(s);
+          nuk = MPFR_GET_EXP(u) - MPFR_GET_EXP(s);
 
-          if ((nuk < - (mp_exp_t) m) && ((double) k >= xf * xf))
+          if ((nuk < - (mp_exp_t) m) && ((double) k >= xf2))
             break;
 
           /* tauk <- 1/2 + tauk * 2^sigmak + (1+8k)*2^nuk */
@@ -180,7 +161,7 @@ mpfr_erf_0 (mpfr_ptr res, mpfr_srcptr x, mp_rnd_t rnd_mode)
         }
 
       mpfr_mul (s, x, s, GMP_RNDU);
-      MPFR_EXP(s) ++;
+      MPFR_SET_EXP (s, MPFR_GET_EXP (s) + 1);
 
       mpfr_const_pi (t, GMP_RNDZ);
       mpfr_sqrt (t, t, GMP_RNDZ);
@@ -188,10 +169,19 @@ mpfr_erf_0 (mpfr_ptr res, mpfr_srcptr x, mp_rnd_t rnd_mode)
       tauk = 4.0 * tauk + 11.0; /* final ulp-error on s */
       log2tauk = __gmpfr_ceil_log2 (tauk);
 
-      ok = mpfr_can_round (s, m - log2tauk, GMP_RNDN, GMP_RNDZ,
-                           n + (rnd_mode == GMP_RNDN));
+      if (MPFR_LIKELY (mpfr_can_round (s, m - log2tauk, GMP_RNDN, GMP_RNDZ,
+				       n + (rnd_mode == GMP_RNDN))))
+	break;
+
+      /* Actualisation of the precision */
+      MPFR_ZIV_NEXT (loop, m);
+      mpfr_set_prec (y, m);
+      mpfr_set_prec (s, m);
+      mpfr_set_prec (t, m);
+      mpfr_set_prec (u, m);
+
     }
-  while (ok == 0);
+  MPFR_ZIV_FREE (loop);
 
   inex = mpfr_set (res, s, rnd_mode);