contrib/mpfr/erfc.c

   1 /* mpfr_erfc -- The Complementary Error Function of a floating-point number
   2
   3 Copyright 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
   4 Contributed by the Arenaire and Cacao 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 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 /* erfc(x) = 1 - erf(x) */
  27
  28 /* Put in y an approximation of erfc(x) for large x, using formulae 7.1.23 and
  29    7.1.24 from Abramowitz and Stegun.
  30    Returns e such that the error is bounded by 2^e ulp(y),
  31    or returns 0 in case of underflow.
  32 */
  33 static mp_exp_t
  34 mpfr_erfc_asympt (mpfr_ptr y, mpfr_srcptr x)
  35 {
  36   mpfr_t t, xx, err;
  37   unsigned long k;
  38   mp_prec_t prec = MPFR_PREC(y);
  39   mp_exp_t exp_err;
  40
  41   mpfr_init2 (t, prec);
  42   mpfr_init2 (xx, prec);
  43   mpfr_init2 (err, 31);
  44   /* let u = 2^(1-p), and let us represent the error as (1+u)^err
  45      with a bound for err */
  46   mpfr_mul (xx, x, x, GMP_RNDD); /* err <= 1 */
  47   mpfr_ui_div (xx, 1, xx, GMP_RNDU); /* upper bound for 1/(2x^2), err <= 2 */
  48   mpfr_div_2ui (xx, xx, 1, GMP_RNDU); /* exact */
  49   mpfr_set_ui (t, 1, GMP_RNDN); /* current term, exact */
  50   mpfr_set (y, t, GMP_RNDN);    /* current sum  */
  51   mpfr_set_ui (err, 0, GMP_RNDN);
  52   for (k = 1; ; k++)
  53     {
  54       mpfr_mul_ui (t, t, 2 * k - 1, GMP_RNDU); /* err <= 4k-3 */
  55       mpfr_mul (t, t, xx, GMP_RNDU);           /* err <= 4k */
  56       /* for -1 < x < 1, and |nx| < 1, we have |(1+x)^n| <= 1+7/4|nx|.
  57          Indeed, for x>=0: log((1+x)^n) = n*log(1+x) <= n*x. Let y=n*x < 1,
  58          then exp(y) <= 1+7/4*y.
  59          For x<=0, let x=-x, we can prove by induction that (1-x)^n >= 1-n*x.*/
  60       mpfr_mul_2si (err, err, MPFR_GET_EXP (y) - MPFR_GET_EXP (t), GMP_RNDU);
  61       mpfr_add_ui (err, err, 14 * k, GMP_RNDU); /* 2^(1-p) * t <= 2 ulp(t) */
  62       mpfr_div_2si (err, err, MPFR_GET_EXP (y) - MPFR_GET_EXP (t), GMP_RNDU);
  63       if (MPFR_GET_EXP (t) + (mp_exp_t) prec <= MPFR_GET_EXP (y))
  64         {
  65           /* the truncation error is bounded by |t| < ulp(y) */
  66           mpfr_add_ui (err, err, 1, GMP_RNDU);
  67           break;
  68         }
  69       if (k & 1)
  70         mpfr_sub (y, y, t, GMP_RNDN);
  71       else
  72         mpfr_add (y, y, t, GMP_RNDN);
  73     }
  74   /* the error on y is bounded by err*ulp(y) */
  75   mpfr_mul (t, x, x, GMP_RNDU); /* rel. err <= 2^(1-p) */
  76   mpfr_div_2ui (err, err, 3, GMP_RNDU);  /* err/8 */
  77   mpfr_add (err, err, t, GMP_RNDU);      /* err/8 + xx */
  78   mpfr_mul_2ui (err, err, 3, GMP_RNDU);  /* err + 8*xx */
  79   mpfr_exp (t, t, GMP_RNDU); /* err <= 1/2*ulp(t) + err(x*x)*t
  80                                 <= 1/2*ulp(t)+2*|x*x|*ulp(t)
  81                                 <= (2*|x*x|+1/2)*ulp(t) */
  82   mpfr_mul (t, t, x, GMP_RNDN); /* err <= 1/2*ulp(t) + (4*|x*x|+1)*ulp(t)
  83                                    <= (4*|x*x|+3/2)*ulp(t) */
  84   mpfr_const_pi (xx, GMP_RNDZ); /* err <= ulp(Pi) */
  85   mpfr_sqrt (xx, xx, GMP_RNDN); /* err <= 1/2*ulp(xx) + ulp(Pi)/2/sqrt(Pi)
  86                                    <= 3/2*ulp(xx) */
  87   mpfr_mul (t, t, xx, GMP_RNDN); /* err <= (8 |xx| + 13/2) * ulp(t) */
  88   mpfr_div (y, y, t, GMP_RNDN); /* the relative error on input y is bounded
  89                                    by (1+u)^err with u = 2^(1-p), that on
  90                                    t is bounded by (1+u)^(8 |xx| + 13/2),
  91                                    thus that on output y is bounded by
  92                                    8 |xx| + 7 + err. */
  93
  94   if (MPFR_IS_ZERO(y))
  95     {
  96       /* If y is zero, most probably we have underflow. We check it directly
  97          using the fact that erfc(x) <= exp(-x^2)/sqrt(Pi)/x for x >= 0.
  98          We compute an upper approximation of exp(-x^2)/sqrt(Pi)/x.
  99       */
 100       mpfr_mul (t, x, x, GMP_RNDD); /* t <= x^2 */
 101       mpfr_neg (t, t, GMP_RNDU);    /* -x^2 <= t */
 102       mpfr_exp (t, t, GMP_RNDU);    /* exp(-x^2) <= t */
 103       mpfr_const_pi (xx, GMP_RNDD); /* xx <= sqrt(Pi), cached */
 104       mpfr_mul (xx, xx, x, GMP_RNDD); /* xx <= sqrt(Pi)*x */
 105       mpfr_div (y, t, xx, GMP_RNDN); /* if y is zero, this means that the upper
 106                                         approximation of exp(-x^2)/sqrt(Pi)/x
 107                                         is nearer from 0 than from 2^(-emin-1),
 108                                         thus we have underflow. */
 109       exp_err = 0;
 110     }
 111   else
 112     {
 113       mpfr_add_ui (err, err, 7, GMP_RNDU);
 114       exp_err = MPFR_GET_EXP (err);
 115     }
 116
 117   mpfr_clear (t);
 118   mpfr_clear (xx);
 119   mpfr_clear (err);
 120   return exp_err;
 121 }
 122
 123 int
 124 mpfr_erfc (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd)
 125 {
 126   int inex;
 127   mpfr_t tmp;
 128   mp_exp_t te, err;
 129   mp_prec_t prec;
 130   MPFR_SAVE_EXPO_DECL (expo);
 131   MPFR_ZIV_DECL (loop);
 132
 133   MPFR_LOG_FUNC (("x[%#R]=%R rnd=%d", x, x, rnd),
 134                  ("y[%#R]=%R inexact=%d", y, y, inex));
 135
 136   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
 137     {
 138       if (MPFR_IS_NAN (x))
 139         {
 140           MPFR_SET_NAN (y);
 141           MPFR_RET_NAN;
 142         }
 143       /* erfc(+inf) = 0+, erfc(-inf) = 2 erfc (0) = 1 */
 144       else if (MPFR_IS_INF (x))
 145         return mpfr_set_ui (y, MPFR_IS_POS (x) ? 0 : 2, rnd);
 146       else
 147         return mpfr_set_ui (y, 1, rnd);
 148     }
 149
 150   if (MPFR_SIGN (x) > 0)
 151     {
 152       /* for x >= 27282, erfc(x) < 2^(-2^30-1) */
 153       if (mpfr_cmp_ui (x, 27282) >= 0)
 154         return mpfr_underflow (y, (rnd == GMP_RNDN) ? GMP_RNDZ : rnd, 1);
 155     }
 156
 157   if (MPFR_SIGN (x) < 0)
 158     {
 159       mp_exp_t e = MPFR_EXP(x);
 160       /* For x < 0 going to -infinity, erfc(x) tends to 2 by below.
 161          More precisely, we have 2 + 1/sqrt(Pi)/x/exp(x^2) < erfc(x) < 2.
 162          Thus log2 |2 - erfc(x)| <= -log2|x| - x^2 / log(2).
 163          If |2 - erfc(x)| < 2^(-PREC(y)) then the result is either 2 or
 164          nextbelow(2).
 165          For x <= -27282, -log2|x| - x^2 / log(2) <= -2^30.
 166       */
 167       if ((MPFR_PREC(y) <= 7 && e >= 2) ||  /* x <= -2 */
 168           (MPFR_PREC(y) <= 25 && e >= 3) || /* x <= -4 */
 169           (MPFR_PREC(y) <= 120 && mpfr_cmp_si (x, -9) <= 0) ||
 170           mpfr_cmp_si (x, -27282) <= 0)
 171         {
 172         near_two:
 173           mpfr_set_ui (y, 2, GMP_RNDN);
 174           mpfr_set_inexflag ();
 175           if (rnd == GMP_RNDZ || rnd == GMP_RNDD)
 176             {
 177               mpfr_nextbelow (y);
 178               return -1;
 179             }
 180           else
 181             return 1;
 182         }
 183       else if (e >= 3) /* more accurate test */
 184         {
 185           mpfr_t t, u;
 186           int near_2;
 187           mpfr_init2 (t, 32);
 188           mpfr_init2 (u, 32);
 189           /* the following is 1/log(2) rounded to zero on 32 bits */
 190           mpfr_set_str_binary (t, "1.0111000101010100011101100101001");
 191           mpfr_sqr (u, x, GMP_RNDZ);
 192           mpfr_mul (t, t, u, GMP_RNDZ); /* t <= x^2/log(2) */
 193           mpfr_neg (u, x, GMP_RNDZ); /* 0 <= u <= |x| */
 194           mpfr_log2 (u, u, GMP_RNDZ); /* u <= log2(|x|) */
 195           mpfr_add (t, t, u, GMP_RNDZ); /* t <= log2|x| + x^2 / log(2) */
 196           near_2 = mpfr_cmp_ui (t, MPFR_PREC(y)) >= 0;
 197           mpfr_clear (t);
 198           mpfr_clear (u);
 199           if (near_2)
 200             goto near_two;
 201         }
 202     }
 203
 204   /* Init stuff */
 205   MPFR_SAVE_EXPO_MARK (expo);
 206
 207   /* erfc(x) ~ 1, with error < 2^(EXP(x)+1) */
 208   MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, __gmpfr_one, - MPFR_GET_EXP (x) - 1,
 209                                     0, MPFR_SIGN(x) < 0,
 210                                     rnd, inex = _inexact; goto end);
 211
 212   prec = MPFR_PREC (y) + MPFR_INT_CEIL_LOG2 (MPFR_PREC (y)) + 3;
 213   if (MPFR_GET_EXP (x) > 0)
 214     prec += 2 * MPFR_GET_EXP(x);
 215
 216   mpfr_init2 (tmp, prec);
 217
 218   MPFR_ZIV_INIT (loop, prec);            /* Initialize the ZivLoop controler */
 219   for (;;)                               /* Infinite loop */
 220     {
 221       /* use asymptotic formula only whenever x^2 >= p*log(2),
 222          otherwise it will not converge */
 223       if (MPFR_SIGN (x) > 0 &&
 224           2 * MPFR_GET_EXP (x) - 2 >= MPFR_INT_CEIL_LOG2 (prec))
 225         /* we have x^2 >= p in that case */
 226         {
 227           err = mpfr_erfc_asympt (tmp, x);
 228           if (err == 0) /* underflow case */
 229             {
 230               mpfr_clear (tmp);
 231               MPFR_SAVE_EXPO_FREE (expo);
 232               return mpfr_underflow (y, (rnd == GMP_RNDN) ? GMP_RNDZ : rnd, 1);
 233             }
 234         }
 235       else
 236         {
 237           mpfr_erf (tmp, x, GMP_RNDN);
 238           MPFR_ASSERTD (!MPFR_IS_SINGULAR (tmp)); /* FIXME: 0 only for x=0 ? */
 239           te = MPFR_GET_EXP (tmp);
 240           mpfr_ui_sub (tmp, 1, tmp, GMP_RNDN);
 241           /* See error analysis in algorithms.tex for details */
 242           if (MPFR_IS_ZERO (tmp))
 243             {
 244               prec *= 2;
 245               err = prec; /* ensures MPFR_CAN_ROUND fails */
 246             }
 247           else
 248             err = MAX (te - MPFR_GET_EXP (tmp), 0) + 1;
 249         }
 250       if (MPFR_LIKELY (MPFR_CAN_ROUND (tmp, prec - err, MPFR_PREC (y), rnd)))
 251         break;
 252       MPFR_ZIV_NEXT (loop, prec);        /* Increase used precision */
 253       mpfr_set_prec (tmp, prec);
 254     }
 255   MPFR_ZIV_FREE (loop);                  /* Free the ZivLoop Controler */
 256
 257   inex = mpfr_set (y, tmp, rnd);    /* Set y to the computed value */
 258   mpfr_clear (tmp);
 259
 260  end:
 261   MPFR_SAVE_EXPO_FREE (expo);
 262   return mpfr_check_range (y, inex, rnd);
 263 }