contrib/mpfr/src/erfc.c

   1 /* mpfr_erfc -- The Complementary Error Function of a floating-point number
   2
   3 Copyright 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc.
   4 Contributed by the Arenaire 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 /* 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 mpfr_exp_t
  34 mpfr_erfc_asympt (mpfr_ptr y, mpfr_srcptr x)
  35 {
  36   mpfr_t t, xx, err;
  37   unsigned long k;
  38   mpfr_prec_t prec = MPFR_PREC(y);
  39   mpfr_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, MPFR_RNDD); /* err <= 1 */
  47   mpfr_ui_div (xx, 1, xx, MPFR_RNDU); /* upper bound for 1/(2x^2), err <= 2 */
  48   mpfr_div_2ui (xx, xx, 1, MPFR_RNDU); /* exact */
  49   mpfr_set_ui (t, 1, MPFR_RNDN); /* current term, exact */
  50   mpfr_set (y, t, MPFR_RNDN);    /* current sum  */
  51   mpfr_set_ui (err, 0, MPFR_RNDN);
  52   for (k = 1; ; k++)
  53     {
  54       mpfr_mul_ui (t, t, 2 * k - 1, MPFR_RNDU); /* err <= 4k-3 */
  55       mpfr_mul (t, t, xx, MPFR_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), MPFR_RNDU);
  61       mpfr_add_ui (err, err, 14 * k, MPFR_RNDU); /* 2^(1-p) * t <= 2 ulp(t) */
  62       mpfr_div_2si (err, err, MPFR_GET_EXP (y) - MPFR_GET_EXP (t), MPFR_RNDU);
  63       if (MPFR_GET_EXP (t) + (mpfr_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, MPFR_RNDU);
  67           break;
  68         }
  69       if (k & 1)
  70         mpfr_sub (y, y, t, MPFR_RNDN);
  71       else
  72         mpfr_add (y, y, t, MPFR_RNDN);
  73     }
  74   /* the error on y is bounded by err*ulp(y) */
  75   mpfr_mul (t, x, x, MPFR_RNDU); /* rel. err <= 2^(1-p) */
  76   mpfr_div_2ui (err, err, 3, MPFR_RNDU);  /* err/8 */
  77   mpfr_add (err, err, t, MPFR_RNDU);      /* err/8 + xx */
  78   mpfr_mul_2ui (err, err, 3, MPFR_RNDU);  /* err + 8*xx */
  79   mpfr_exp (t, t, MPFR_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, MPFR_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, MPFR_RNDZ); /* err <= ulp(Pi) */
  85   mpfr_sqrt (xx, xx, MPFR_RNDN); /* err <= 1/2*ulp(xx) + ulp(Pi)/2/sqrt(Pi)
  86                                    <= 3/2*ulp(xx) */
  87   mpfr_mul (t, t, xx, MPFR_RNDN); /* err <= (8 |xx| + 13/2) * ulp(t) */
  88   mpfr_div (y, y, t, MPFR_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, MPFR_RNDD); /* t <= x^2 */
 101       mpfr_neg (t, t, MPFR_RNDU);    /* -x^2 <= t */
 102       mpfr_exp (t, t, MPFR_RNDU);    /* exp(-x^2) <= t */
 103       mpfr_const_pi (xx, MPFR_RNDD); /* xx <= sqrt(Pi), cached */
 104       mpfr_mul (xx, xx, x, MPFR_RNDD); /* xx <= sqrt(Pi)*x */
 105       mpfr_div (y, t, xx, MPFR_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, MPFR_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, mpfr_rnd_t rnd)
 125 {
 126   int inex;
 127   mpfr_t tmp;
 128   mpfr_exp_t te, err;
 129   mpfr_prec_t prec;
 130   mpfr_exp_t emin = mpfr_get_emin ();
 131   MPFR_SAVE_EXPO_DECL (expo);
 132   MPFR_ZIV_DECL (loop);
 133
 134   MPFR_LOG_FUNC
 135     (("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, rnd),
 136      ("y[%Pu]=%.*Rg inexact=%d", mpfr_get_prec (y), mpfr_log_prec, y, inex));
 137
 138   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
 139     {
 140       if (MPFR_IS_NAN (x))
 141         {
 142           MPFR_SET_NAN (y);
 143           MPFR_RET_NAN;
 144         }
 145       /* erfc(+inf) = 0+, erfc(-inf) = 2 erfc (0) = 1 */
 146       else if (MPFR_IS_INF (x))
 147         return mpfr_set_ui (y, MPFR_IS_POS (x) ? 0 : 2, rnd);
 148       else
 149         return mpfr_set_ui (y, 1, rnd);
 150     }
 151
 152   if (MPFR_SIGN (x) > 0)
 153     {
 154       /* by default, emin = 1-2^30, thus the smallest representable
 155          number is 1/2*2^emin = 2^(-2^30):
 156          for x >= 27282, erfc(x) < 2^(-2^30-1), and
 157          for x >= 1787897414, erfc(x) < 2^(-2^62-1).
 158       */
 159       if ((emin >= -1073741823 && mpfr_cmp_ui (x, 27282) >= 0) ||
 160           mpfr_cmp_ui (x, 1787897414) >= 0)
 161         {
 162           /* May be incorrect if MPFR_EMAX_MAX >= 2^62. */
 163           MPFR_ASSERTN ((MPFR_EMAX_MAX >> 31) >> 31 == 0);
 164           return mpfr_underflow (y, (rnd == MPFR_RNDN) ? MPFR_RNDZ : rnd, 1);
 165         }
 166     }
 167
 168   /* Init stuff */
 169   MPFR_SAVE_EXPO_MARK (expo);
 170
 171   if (MPFR_SIGN (x) < 0)
 172     {
 173       mpfr_exp_t e = MPFR_EXP(x);
 174       /* For x < 0 going to -infinity, erfc(x) tends to 2 by below.
 175          More precisely, we have 2 + 1/sqrt(Pi)/x/exp(x^2) < erfc(x) < 2.
 176          Thus log2 |2 - erfc(x)| <= -log2|x| - x^2 / log(2).
 177          If |2 - erfc(x)| < 2^(-PREC(y)) then the result is either 2 or
 178          nextbelow(2).
 179          For x <= -27282, -log2|x| - x^2 / log(2) <= -2^30.
 180       */
 181       if ((MPFR_PREC(y) <= 7 && e >= 2) ||  /* x <= -2 */
 182           (MPFR_PREC(y) <= 25 && e >= 3) || /* x <= -4 */
 183           (MPFR_PREC(y) <= 120 && mpfr_cmp_si (x, -9) <= 0) ||
 184           mpfr_cmp_si (x, -27282) <= 0)
 185         {
 186         near_two:
 187           mpfr_set_ui (y, 2, MPFR_RNDN);
 188           mpfr_set_inexflag ();
 189           if (rnd == MPFR_RNDZ || rnd == MPFR_RNDD)
 190             {
 191               mpfr_nextbelow (y);
 192               inex = -1;
 193             }
 194           else
 195             inex = 1;
 196           goto end;
 197         }
 198       else if (e >= 3) /* more accurate test */
 199         {
 200           mpfr_t t, u;
 201           int near_2;
 202           mpfr_init2 (t, 32);
 203           mpfr_init2 (u, 32);
 204           /* the following is 1/log(2) rounded to zero on 32 bits */
 205           mpfr_set_str_binary (t, "1.0111000101010100011101100101001");
 206           mpfr_sqr (u, x, MPFR_RNDZ);
 207           mpfr_mul (t, t, u, MPFR_RNDZ); /* t <= x^2/log(2) */
 208           mpfr_neg (u, x, MPFR_RNDZ); /* 0 <= u <= |x| */
 209           mpfr_log2 (u, u, MPFR_RNDZ); /* u <= log2(|x|) */
 210           mpfr_add (t, t, u, MPFR_RNDZ); /* t <= log2|x| + x^2 / log(2) */
 211           /* Taking into account that mpfr_exp_t >= mpfr_prec_t */
 212           mpfr_set_exp_t (u, MPFR_PREC (y), MPFR_RNDU);
 213           near_2 = mpfr_cmp (t, u) >= 0;  /* 1 if PREC(y) <= u <= t <= ... */
 214           mpfr_clear (t);
 215           mpfr_clear (u);
 216           if (near_2)
 217             goto near_two;
 218         }
 219     }
 220
 221   /* erfc(x) ~ 1, with error < 2^(EXP(x)+1) */
 222   MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, __gmpfr_one, - MPFR_GET_EXP (x) - 1,
 223                                     0, MPFR_SIGN(x) < 0,
 224                                     rnd, inex = _inexact; goto end);
 225
 226   prec = MPFR_PREC (y) + MPFR_INT_CEIL_LOG2 (MPFR_PREC (y)) + 3;
 227   if (MPFR_GET_EXP (x) > 0)
 228     prec += 2 * MPFR_GET_EXP(x);
 229
 230   mpfr_init2 (tmp, prec);
 231
 232   MPFR_ZIV_INIT (loop, prec);            /* Initialize the ZivLoop controler */
 233   for (;;)                               /* Infinite loop */
 234     {
 235       /* use asymptotic formula only whenever x^2 >= p*log(2),
 236          otherwise it will not converge */
 237       if (MPFR_SIGN (x) > 0 &&
 238           2 * MPFR_GET_EXP (x) - 2 >= MPFR_INT_CEIL_LOG2 (prec))
 239         /* we have x^2 >= p in that case */
 240         {
 241           err = mpfr_erfc_asympt (tmp, x);
 242           if (err == 0) /* underflow case */
 243             {
 244               mpfr_clear (tmp);
 245               MPFR_SAVE_EXPO_FREE (expo);
 246               return mpfr_underflow (y, (rnd == MPFR_RNDN) ? MPFR_RNDZ : rnd, 1);
 247             }
 248         }
 249       else
 250         {
 251           mpfr_erf (tmp, x, MPFR_RNDN);
 252           MPFR_ASSERTD (!MPFR_IS_SINGULAR (tmp)); /* FIXME: 0 only for x=0 ? */
 253           te = MPFR_GET_EXP (tmp);
 254           mpfr_ui_sub (tmp, 1, tmp, MPFR_RNDN);
 255           /* See error analysis in algorithms.tex for details */
 256           if (MPFR_IS_ZERO (tmp))
 257             {
 258               prec *= 2;
 259               err = prec; /* ensures MPFR_CAN_ROUND fails */
 260             }
 261           else
 262             err = MAX (te - MPFR_GET_EXP (tmp), 0) + 1;
 263         }
 264       if (MPFR_LIKELY (MPFR_CAN_ROUND (tmp, prec - err, MPFR_PREC (y), rnd)))
 265         break;
 266       MPFR_ZIV_NEXT (loop, prec);        /* Increase used precision */
 267       mpfr_set_prec (tmp, prec);
 268     }
 269   MPFR_ZIV_FREE (loop);                  /* Free the ZivLoop Controler */
 270
 271   inex = mpfr_set (y, tmp, rnd);    /* Set y to the computed value */
 272   mpfr_clear (tmp);
 273
 274  end:
 275   MPFR_SAVE_EXPO_FREE (expo);
 276   return mpfr_check_range (y, inex, rnd);
 277 }