implemented asymptotic formula for erfc (fixed both slowness for large

arguments, and call to MPFR_WARNING with return value NaN for huge arguments)


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

4 files changed, 128 insertions, 29 deletions

diff --git a/BUGS b/BUGS
index b84f1c827..4d2ccb983 100644
--- a/BUGS
+++ b/BUGS
@@ -33,12 +33,9 @@ Known bugs:
   underflow flag is not always set. Example: x = 0.11 in base 2,
   y = 1.0e38, GMP_RNDN.
 
-* On some very large values, mpfr_erfc crashes (segmentation fault,
-  probably due to lack of memory or integer overflow).
-
-* Some large input values for mpfr_eint and mpfr_erfc are currently not
-  supported. However they should be detected, so that one can have the
-  following documented behavior: NaN is returned and the erange flag is
+* Some large input values for mpfr_eint are currently not supported.
+  However they should be detected, so that one can have the following
+  documented behavior: NaN is returned and the erange flag is
   set to indicate that the input is in an unsupported domain. If MPFR
   has been built with --enable-warnings and the environment variable
   MPFR_QUIET is null, a warning is output to the standard error stream.
@@ -57,7 +54,9 @@ Potential bugs:
 
 * Possible infinite loop in some functions for particular cases: when
   the exact result is an exactly representable number or the middle of
-  consecutive two such numbers (for example, erf).
+  consecutive two such numbers. However for non-algebraic functions, it is
+  believed that no such case exists, except the well-known cases like cos(0)=1,
+  exp(0)=1, and so on, and the x^y function when y is an integer or y=1/2^k.
 
 * The mpfr_set_ld function may be quite slow if the long double type has an
   exponent of more than 15 bits.
diff --git a/algorithms.tex b/algorithms.tex
index 2ec2d26d0..f1cf2212b 100644
--- a/algorithms.tex
+++ b/algorithms.tex
@@ -1075,6 +1075,12 @@ that ${\rm erfc} \, x \leq 2^{-n}$, thus
 ${\rm erf} \, x = 1$ or ${\rm erf} \, x = 1 - 2^{-n}$ according to the
 rounding mode.
 
+More precisely, \cite[formul{\ae} 7.1.23 and 7.1.24]{AbSt73} give:
+\[ \sqrt{\pi} x e^{x^2} {\rm erfc} x \approx
+   1 + \sum_{k=1}^n (-1)^k \frac{1 \times 3 \times \cdots \times (2k-1)}
+   {(2x^2)^k}, \]
+with the error bounded in absolute value by the next term and of the same sign.
+
 \subsection{The hyperbolic cosine function}
 
 The {\tt mpfr\_cosh} ($\cosh{x}$) function implements the hyperbolic
diff --git a/erfc.c b/erfc.c
index 0a01e3420..43ce09e9d 100644
--- a/erfc.c
+++ b/erfc.c
@@ -19,13 +19,82 @@ along with the MPFR Library; see the file COPYING.LIB.  If not, write to
 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
 MA 02110-1301, USA. */
 
-#include <stdio.h>   /* for MPFR_WARNING */
-#include <stdlib.h>  /* for MPFR_WARNING */
 #define MPFR_NEED_LONGLONG_H
 #include "mpfr-impl.h"
 
 /* erfc(x) = 1 - erf(x) */
 
+/* Put in y an approximation of erfc(x) for large x, using formulae 7.1.23 and
+   7.1.24 from Abramowitz and Stegun.
+   Returns e such that the error is bounded by 2^e ulp(y). */
+static mp_exp_t
+mpfr_erfc_asympt (mpfr_ptr y, mpfr_srcptr x)
+{
+  mpfr_t t, xx, err;
+  unsigned long k;
+  mp_prec_t prec = MPFR_PREC(y);
+  mp_exp_t exp_err;
+
+  mpfr_init2 (t, prec);
+  mpfr_init2 (xx, prec);
+  mpfr_init2 (err, 31);
+  /* let u = 2^(1-p), and let us represent the error as (1+u)^err
+     with a bound for err */
+  mpfr_mul (xx, x, x, GMP_RNDD); /* err <= 1 */
+  mpfr_ui_div (xx, 1, xx, GMP_RNDU); /* upper bound for 1/(2x^2), err <= 2 */
+  mpfr_div_2exp (xx, xx, 1, GMP_RNDU); /* exact */
+  mpfr_set_ui (t, 1, GMP_RNDN); /* current term, exact */
+  mpfr_set (y, t, GMP_RNDN);    /* current sum  */
+  mpfr_set_ui (err, 0, GMP_RNDN);
+  for (k = 1; ; k++)
+    {
+      mpfr_mul_ui (t, t, 2 * k - 1, GMP_RNDU); /* err <= 4k-3 */
+      mpfr_mul (t, t, xx, GMP_RNDU);           /* err <= 4k */
+      /* for -1 < x < 1, and |nx| < 1, we have |(1+x)^n| <= 1+7/4|nx|.
+         Indeed, for x>=0: log((1+x)^n) = n*log(1+x) <= n*x. Let y=n*x < 1,
+         then exp(y) <= 1+7/4*y.
+         For x<=0, let x=-x, we can prove by induction that (1-x)^n >= 1-n*x.*/
+      mpfr_mul_2si (err, err, MPFR_GET_EXP (y) - MPFR_GET_EXP (t), GMP_RNDU);
+      mpfr_add_ui (err, err, 14 * k, GMP_RNDU); /* 2^(1-p) * t <= 2 ulp(t) */
+      mpfr_div_2si (err, err, MPFR_GET_EXP (y) - MPFR_GET_EXP (t), GMP_RNDU);
+      if (MPFR_GET_EXP (t) + (mp_exp_t) prec <= MPFR_GET_EXP (y))
+        {
+          /* the truncation error is bounded by |t| < ulp(y) */
+          mpfr_add_ui (err, err, 1, GMP_RNDU);
+          break;
+        }
+      if (k & 2 != 0)
+        mpfr_sub (y, y, t, GMP_RNDN);
+      else
+        mpfr_add (y, y, t, GMP_RNDN);
+    }
+  /* the error on y is bounded by err*ulp(y) */
+  mpfr_mul (t, x, x, GMP_RNDU); /* rel. err <= 2^(1-p) */
+  mpfr_div_2exp (err, err, 3, GMP_RNDU); /* err/8 */
+  mpfr_add (err, err, t, GMP_RNDU);      /* err/8 + xx */
+  mpfr_mul_2exp (err, err, 3, GMP_RNDU); /* err + 8*xx */
+  mpfr_exp (t, t, GMP_RNDU); /* err <= 1/2*ulp(t) + err(x*x)*t
+                                <= 1/2*ulp(t)+2*|x*x|*ulp(t)
+                                <= (2*|x*x|+1/2)*ulp(t) */
+  mpfr_mul (t, t, x, GMP_RNDN); /* err <= 1/2*ulp(t) + (4*|x*x|+1)*ulp(t)
+                                   <= (4*|x*x|+3/2)*ulp(t) */
+  mpfr_const_pi (xx, GMP_RNDZ); /* err <= ulp(Pi) */
+  mpfr_sqrt (xx, xx, GMP_RNDN); /* err <= 1/2*ulp(xx) + ulp(Pi)/2/sqrt(Pi)
+                                   <= 3/2*ulp(xx) */
+  mpfr_mul (t, t, xx, GMP_RNDN); /* err <= (8 |xx| + 13/2) * ulp(t) */
+  mpfr_div (y, y, t, GMP_RNDN); /* the relative error on input y is bounded
+                                   by (1+u)^err with u = 2^(1-p), that on
+                                   t is bounded by (1+u)^(8 |xx| + 13/2),
+                                   thus that on output y is bounded by
+                                   8 |xx| + 7 + err. */
+  mpfr_add_ui (err, err, 7, GMP_RNDU);
+  exp_err = MPFR_GET_EXP (err);
+  mpfr_clear (t);
+  mpfr_clear (xx);
+  mpfr_clear (err);
+  return exp_err;
+}
+
 int
 mpfr_erfc (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd)
 {
@@ -58,16 +127,6 @@ mpfr_erfc (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd)
       /* for x >= 38582, erfc(x) < 2^(-2^31) */
       if (mpfr_cmp_ui (x, 38582) >= 0)
         return mpfr_underflow (y, (rnd == GMP_RNDN) ? GMP_RNDZ : rnd, 1);
-      if (MPFR_GET_EXP (x) >= 12)
-        {
-          /* FIXME: Improve the algorithm to be able to compute the actual
-             value. For the time being, we regard this as a range error,
-             so that the caller can cleanly deal with the problem. */
-          MPFR_WARNING ("Too large input in mpfr_erfc, returning NaN");
-          MPFR_SET_ERANGE ();
-          MPFR_SET_NAN (y);
-          MPFR_RET_NAN;
-        }
     }
 
   if (MPFR_SIGN (x) < 0)
@@ -98,24 +157,36 @@ mpfr_erfc (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd)
                                     rnd, inex = _inexact; goto end);
 
   prec = MPFR_PREC (y) + MPFR_INT_CEIL_LOG2 (MPFR_PREC (y)) + 3;
+  if (MPFR_GET_EXP (x) > 0)
+    prec += 2 * MPFR_GET_EXP(x);
+
   mpfr_init2 (tmp, prec);
 
   MPFR_ZIV_INIT (loop, prec);            /* Initialize the ZivLoop controler */
   for (;;)                               /* Infinite loop */
     {
-      mpfr_erf (tmp, x, GMP_RNDN);
-      MPFR_ASSERTD (!MPFR_IS_SINGULAR (tmp)); /* FIXME: 0 only for x=0 ? */
-      te = MPFR_GET_EXP (tmp);
-      mpfr_ui_sub (tmp, 1, tmp, GMP_RNDN);
-      /* See error analysis of expm1 for details */
-      if (MPFR_IS_ZERO (tmp))
-        prec *=2;
+      /* use asymptotic formula only whenever x^2 >= p*log(2),
+         otherwise it will not converge */
+      if (2 * MPFR_GET_EXP (x) - 2 >= MPFR_INT_CEIL_LOG2 (prec))
+        /* we have x^2 >= p in that case */
+        err = mpfr_erfc_asympt (tmp, x);
       else
         {
-          err = prec - (MAX (te - MPFR_GET_EXP (tmp), 0) + 1);
-          if (MPFR_LIKELY (MPFR_CAN_ROUND (tmp, err, MPFR_PREC (y), rnd)))
-            break;
+          mpfr_erf (tmp, x, GMP_RNDN);
+          MPFR_ASSERTD (!MPFR_IS_SINGULAR (tmp)); /* FIXME: 0 only for x=0 ? */
+          te = MPFR_GET_EXP (tmp);
+          mpfr_ui_sub (tmp, 1, tmp, GMP_RNDN);
+          /* See error analysis in algorithms.tex for details */
+          if (MPFR_IS_ZERO (tmp))
+            {
+              prec *= 2;
+              err = prec; /* ensures MPFR_CAN_ROUND fails */
+            }
+          else
+            err = MAX (te - MPFR_GET_EXP (tmp), 0) + 1;
         }
+      if (MPFR_LIKELY (MPFR_CAN_ROUND (tmp, prec - err, MPFR_PREC (y), rnd)))
+        break;
       MPFR_ZIV_NEXT (loop, prec);        /* Increase used precision */
       mpfr_set_prec (tmp, prec);
     }
diff --git a/tests/terf.c b/tests/terf.c
index 6d5ec883b..5f6b2178c 100644
--- a/tests/terf.c
+++ b/tests/terf.c
@@ -390,11 +390,34 @@ large_arg (void)
 
   mpfr_set_si_2exp (x, -1, 173, GMP_RNDN);
   mpfr_erfc (y, x, GMP_RNDN);
+  if (mpfr_cmp_ui (y, 2) != 0)
+    {
+      printf ("mpfr_erfc failed for large x (1)\n");
+      exit (1);
+    }
+
 
   mpfr_set_prec (x, 33);
   mpfr_set_prec (y, 43);
   mpfr_set_str_binary (x, "1.11000101010111011000111100101001e6");
   mpfr_erfc (y, x, GMP_RNDD);
+  mpfr_set_prec (x, 43);
+  mpfr_set_str_binary (x, "100010011100101100001101100101011101101E-18579");
+  if (mpfr_cmp (x, y) != 0)
+    {
+      printf ("mpfr_erfc failed for large x (2)\n");
+      exit (1);
+    }
+
+  mpfr_set_prec (y, 43);
+  mpfr_set_si_2exp (x, 1, 11, GMP_RNDN);
+  mpfr_erfc (y, x, GMP_RNDN);
+  mpfr_set_str_binary (x, "0.1100000100100010101111001111010010001000110E-6051113");
+  if (mpfr_cmp (x, y) != 0)
+    {
+      printf ("mpfr_erfc failed for large x (3)\n");
+      exit (1);
+    }
 
   mpfr_clear (x);
   mpfr_clear (y);