contrib/mpfr/erf.c

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