contrib/mpfr/src/erf.c

   1 /* mpfr_erf -- error function of a floating-point number
   2
   3 Copyright 2001, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
   4 Contributed by the AriC and Caramel projects, INRIA.
   5
   6 This file is part of the GNU MPFR Library.
   7
   8 The GNU MPFR Library is free software; you can redistribute it and/or modify
   9 it under the terms of the GNU Lesser General Public License as published by
  10 the Free Software Foundation; either version 3 of the License, or (at your
  11 option) any later version.
  12
  13 The GNU MPFR Library is distributed in the hope that it will be useful, but
  14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  15 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  16 License for more details.
  17
  18 You should have received a copy of the GNU Lesser General Public License
  19 along with the GNU MPFR Library; see the file COPYING.LESSER.  If not, see
  20 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
  21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
  22
  23 #define MPFR_NEED_LONGLONG_H
  24 #include "mpfr-impl.h"
  25
  26 #define EXP1 2.71828182845904523536 /* exp(1) */
  27
  28 static int mpfr_erf_0 (mpfr_ptr, mpfr_srcptr, double, mpfr_rnd_t);
  29
  30 int
  31 mpfr_erf (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
  32 {
  33   mpfr_t xf;
  34   int inex, large;
  35   MPFR_SAVE_EXPO_DECL (expo);
  36
  37   MPFR_LOG_FUNC
  38     (("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, rnd_mode),
  39      ("y[%Pu]=%.*Rg inexact=%d", mpfr_get_prec (y), mpfr_log_prec, y, inex));
  40
  41   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
  42     {
  43       if (MPFR_IS_NAN (x))
  44         {
  45           MPFR_SET_NAN (y);
  46           MPFR_RET_NAN;
  47         }
  48       else if (MPFR_IS_INF (x)) /* erf(+inf) = +1, erf(-inf) = -1 */
  49         return mpfr_set_si (y, MPFR_INT_SIGN (x), MPFR_RNDN);
  50       else /* erf(+0) = +0, erf(-0) = -0 */
  51         {
  52           MPFR_ASSERTD (MPFR_IS_ZERO (x));
  53           return mpfr_set (y, x, MPFR_RNDN); /* should keep the sign of x */
  54         }
  55     }
  56
  57   /* now x is neither NaN, Inf nor 0 */
  58
  59   /* first try expansion at x=0 when x is small, or asymptotic expansion
  60      where x is large */
  61
  62   MPFR_SAVE_EXPO_MARK (expo);
  63
  64   /* around x=0, we have erf(x) = 2x/sqrt(Pi) (1 - x^2/3 + ...),
  65      with 1 - x^2/3 <= sqrt(Pi)*erf(x)/2/x <= 1 for x >= 0. This means that
  66      if x^2/3 < 2^(-PREC(y)-1) we can decide of the correct rounding,
  67      unless we have a worst-case for 2x/sqrt(Pi). */
  68   if (MPFR_EXP(x) < - (mpfr_exp_t) (MPFR_PREC(y) / 2))
  69     {
  70       /* we use 2x/sqrt(Pi) (1 - x^2/3) <= erf(x) <= 2x/sqrt(Pi) for x > 0
  71          and 2x/sqrt(Pi) <= erf(x) <= 2x/sqrt(Pi) (1 - x^2/3) for x < 0.
  72          In both cases |2x/sqrt(Pi) (1 - x^2/3)| <= |erf(x)| <= |2x/sqrt(Pi)|.
  73          We will compute l and h such that l <= |2x/sqrt(Pi) (1 - x^2/3)|
  74          and |2x/sqrt(Pi)| <= h. If l and h round to the same value to
  75          precision PREC(y) and rounding rnd_mode, then we are done. */
  76       mpfr_t l, h; /* lower and upper bounds for erf(x) */
  77       int ok, inex2;
  78
  79       mpfr_init2 (l, MPFR_PREC(y) + 17);
  80       mpfr_init2 (h, MPFR_PREC(y) + 17);
  81       /* first compute l */
  82       mpfr_mul (l, x, x, MPFR_RNDU);
  83       mpfr_div_ui (l, l, 3, MPFR_RNDU); /* upper bound on x^2/3 */
  84       mpfr_ui_sub (l, 1, l, MPFR_RNDZ); /* lower bound on 1 - x^2/3 */
  85       mpfr_const_pi (h, MPFR_RNDU); /* upper bound of Pi */
  86       mpfr_sqrt (h, h, MPFR_RNDU); /* upper bound on sqrt(Pi) */
  87       mpfr_div (l, l, h, MPFR_RNDZ); /* lower bound on 1/sqrt(Pi) (1 - x^2/3) */
  88       mpfr_mul_2ui (l, l, 1, MPFR_RNDZ); /* 2/sqrt(Pi) (1 - x^2/3) */
  89       mpfr_mul (l, l, x, MPFR_RNDZ); /* |l| is a lower bound on
  90                                        |2x/sqrt(Pi) (1 - x^2/3)| */
  91       /* now compute h */
  92       mpfr_const_pi (h, MPFR_RNDD); /* lower bound on Pi */
  93       mpfr_sqrt (h, h, MPFR_RNDD); /* lower bound on sqrt(Pi) */
  94       mpfr_div_2ui (h, h, 1, MPFR_RNDD); /* lower bound on sqrt(Pi)/2 */
  95       /* since sqrt(Pi)/2 < 1, the following should not underflow */
  96       mpfr_div (h, x, h, MPFR_IS_POS(x) ? MPFR_RNDU : MPFR_RNDD);
  97       /* round l and h to precision PREC(y) */
  98       inex = mpfr_prec_round (l, MPFR_PREC(y), rnd_mode);
  99       inex2 = mpfr_prec_round (h, MPFR_PREC(y), rnd_mode);
 100       /* Caution: we also need inex=inex2 (inex might be 0). */
 101       ok = SAME_SIGN (inex, inex2) && mpfr_cmp (l, h) == 0;
 102       if (ok)
 103         mpfr_set (y, h, rnd_mode);
 104       mpfr_clear (l);
 105       mpfr_clear (h);
 106       if (ok)
 107         goto end;
 108       /* this test can still fail for small precision, for example
 109          for x=-0.100E-2 with a target precision of 3 bits, since
 110          the error term x^2/3 is not that small. */
 111     }
 112
 113   mpfr_init2 (xf, 53);
 114   mpfr_const_log2 (xf, MPFR_RNDU);
 115   mpfr_div (xf, x, xf, MPFR_RNDZ); /* round to zero ensures we get a lower
 116                                      bound of |x/log(2)| */
 117   mpfr_mul (xf, xf, x, MPFR_RNDZ);
 118   large = mpfr_cmp_ui (xf, MPFR_PREC (y) + 1) > 0;
 119   mpfr_clear (xf);
 120
 121   /* when x goes to infinity, we have erf(x) = 1 - 1/sqrt(Pi)/exp(x^2)/x + ...
 122      and |erf(x) - 1| <= exp(-x^2) is true for any x >= 0, thus if
 123      exp(-x^2) < 2^(-PREC(y)-1) the result is 1 or 1-epsilon.
 124      This rewrites as x^2/log(2) > p+1. */
 125   if (MPFR_UNLIKELY (large))
 126     /* |erf x| = 1 or 1- */
 127     {
 128       mpfr_rnd_t rnd2 = MPFR_IS_POS (x) ? rnd_mode : MPFR_INVERT_RND(rnd_mode);
 129       if (rnd2 == MPFR_RNDN || rnd2 == MPFR_RNDU || rnd2 == MPFR_RNDA)
 130         {
 131           inex = MPFR_INT_SIGN (x);
 132           mpfr_set_si (y, inex, rnd2);
 133         }
 134       else /* round to zero */
 135         {
 136           inex = -MPFR_INT_SIGN (x);
 137           mpfr_setmax (y, 0); /* warning: setmax keeps the old sign of y */
 138           MPFR_SET_SAME_SIGN (y, x);
 139         }
 140     }
 141   else  /* use Taylor */
 142     {
 143       double xf2;
 144
 145       /* FIXME: get rid of doubles/mpfr_get_d here */
 146       xf2 = mpfr_get_d (x, MPFR_RNDN);
 147       xf2 = xf2 * xf2; /* xf2 ~ x^2 */
 148       inex = mpfr_erf_0 (y, x, xf2, rnd_mode);
 149     }
 150
 151  end:
 152   MPFR_SAVE_EXPO_FREE (expo);
 153   return mpfr_check_range (y, inex, rnd_mode);
 154 }
 155
 156 /* return x*2^e */
 157 static double
 158 mul_2exp (double x, mpfr_exp_t e)
 159 {
 160   if (e > 0)
 161     {
 162       while (e--)
 163         x *= 2.0;
 164     }
 165   else
 166     {
 167       while (e++)
 168         x /= 2.0;
 169     }
 170
 171   return x;
 172 }
 173
 174 /* evaluates erf(x) using the expansion at x=0:
 175
 176    erf(x) = 2/sqrt(Pi) * sum((-1)^k*x^(2k+1)/k!/(2k+1), k=0..infinity)
 177
 178    Assumes x is neither NaN nor infinite nor zero.
 179    Assumes also that e*x^2 <= n (target precision).
 180  */
 181 static int
 182 mpfr_erf_0 (mpfr_ptr res, mpfr_srcptr x, double xf2, mpfr_rnd_t rnd_mode)
 183 {
 184   mpfr_prec_t n, m;
 185   mpfr_exp_t nuk, sigmak;
 186   double tauk;
 187   mpfr_t y, s, t, u;
 188   unsigned int k;
 189   int log2tauk;
 190   int inex;
 191   MPFR_ZIV_DECL (loop);
 192
 193   n = MPFR_PREC (res); /* target precision */
 194
 195   /* initial working precision */
 196   m = n + (mpfr_prec_t) (xf2 / LOG2) + 8 + MPFR_INT_CEIL_LOG2 (n);
 197
 198   mpfr_init2 (y, m);
 199   mpfr_init2 (s, m);
 200   mpfr_init2 (t, m);
 201   mpfr_init2 (u, m);
 202
 203   MPFR_ZIV_INIT (loop, m);
 204   for (;;)
 205     {
 206       mpfr_mul (y, x, x, MPFR_RNDU); /* err <= 1 ulp */
 207       mpfr_set_ui (s, 1, MPFR_RNDN);
 208       mpfr_set_ui (t, 1, MPFR_RNDN);
 209       tauk = 0.0;
 210
 211       for (k = 1; ; k++)
 212         {
 213           mpfr_mul (t, y, t, MPFR_RNDU);
 214           mpfr_div_ui (t, t, k, MPFR_RNDU);
 215           mpfr_div_ui (u, t, 2 * k + 1, MPFR_RNDU);
 216           sigmak = MPFR_GET_EXP (s);
 217           if (k % 2)
 218             mpfr_sub (s, s, u, MPFR_RNDN);
 219           else
 220             mpfr_add (s, s, u, MPFR_RNDN);
 221           sigmak -= MPFR_GET_EXP(s);
 222           nuk = MPFR_GET_EXP(u) - MPFR_GET_EXP(s);
 223
 224           if ((nuk < - (mpfr_exp_t) m) && ((double) k >= xf2))
 225             break;
 226
 227           /* tauk <- 1/2 + tauk * 2^sigmak + (1+8k)*2^nuk */
 228           tauk = 0.5 + mul_2exp (tauk, sigmak)
 229             + mul_2exp (1.0 + 8.0 * (double) k, nuk);
 230         }
 231
 232       mpfr_mul (s, x, s, MPFR_RNDU);
 233       MPFR_SET_EXP (s, MPFR_GET_EXP (s) + 1);
 234
 235       mpfr_const_pi (t, MPFR_RNDZ);
 236       mpfr_sqrt (t, t, MPFR_RNDZ);
 237       mpfr_div (s, s, t, MPFR_RNDN);
 238       tauk = 4.0 * tauk + 11.0; /* final ulp-error on s */
 239       log2tauk = __gmpfr_ceil_log2 (tauk);
 240
 241       if (MPFR_LIKELY (MPFR_CAN_ROUND (s, m - log2tauk, n, rnd_mode)))
 242         break;
 243
 244       /* Actualisation of the precision */
 245       MPFR_ZIV_NEXT (loop, m);
 246       mpfr_set_prec (y, m);
 247       mpfr_set_prec (s, m);
 248       mpfr_set_prec (t, m);
 249       mpfr_set_prec (u, m);
 250
 251     }
 252   MPFR_ZIV_FREE (loop);
 253
 254   inex = mpfr_set (res, s, rnd_mode);
 255
 256   mpfr_clear (y);
 257   mpfr_clear (t);
 258   mpfr_clear (u);
 259   mpfr_clear (s);
 260
 261   return inex;
 262 }