14a238c70SJohn Marino /* mpfr_erfc -- The Complementary Error Function of a floating-point number
24a238c70SJohn Marino
3*ab6d115fSJohn Marino Copyright 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
4*ab6d115fSJohn Marino Contributed by the AriC and Caramel projects, INRIA.
54a238c70SJohn Marino
64a238c70SJohn Marino This file is part of the GNU MPFR Library.
74a238c70SJohn Marino
84a238c70SJohn Marino The GNU MPFR Library is free software; you can redistribute it and/or modify
94a238c70SJohn Marino it under the terms of the GNU Lesser General Public License as published by
104a238c70SJohn Marino the Free Software Foundation; either version 3 of the License, or (at your
114a238c70SJohn Marino option) any later version.
124a238c70SJohn Marino
134a238c70SJohn Marino The GNU MPFR Library is distributed in the hope that it will be useful, but
144a238c70SJohn Marino WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
154a238c70SJohn Marino or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
164a238c70SJohn Marino License for more details.
174a238c70SJohn Marino
184a238c70SJohn Marino You should have received a copy of the GNU Lesser General Public License
194a238c70SJohn Marino along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
204a238c70SJohn Marino http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
214a238c70SJohn Marino 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
224a238c70SJohn Marino
234a238c70SJohn Marino #define MPFR_NEED_LONGLONG_H
244a238c70SJohn Marino #include "mpfr-impl.h"
254a238c70SJohn Marino
264a238c70SJohn Marino /* erfc(x) = 1 - erf(x) */
274a238c70SJohn Marino
284a238c70SJohn Marino /* Put in y an approximation of erfc(x) for large x, using formulae 7.1.23 and
294a238c70SJohn Marino 7.1.24 from Abramowitz and Stegun.
304a238c70SJohn Marino Returns e such that the error is bounded by 2^e ulp(y),
314a238c70SJohn Marino or returns 0 in case of underflow.
324a238c70SJohn Marino */
334a238c70SJohn Marino static mpfr_exp_t
mpfr_erfc_asympt(mpfr_ptr y,mpfr_srcptr x)344a238c70SJohn Marino mpfr_erfc_asympt (mpfr_ptr y, mpfr_srcptr x)
354a238c70SJohn Marino {
364a238c70SJohn Marino mpfr_t t, xx, err;
374a238c70SJohn Marino unsigned long k;
384a238c70SJohn Marino mpfr_prec_t prec = MPFR_PREC(y);
394a238c70SJohn Marino mpfr_exp_t exp_err;
404a238c70SJohn Marino
414a238c70SJohn Marino mpfr_init2 (t, prec);
424a238c70SJohn Marino mpfr_init2 (xx, prec);
434a238c70SJohn Marino mpfr_init2 (err, 31);
444a238c70SJohn Marino /* let u = 2^(1-p), and let us represent the error as (1+u)^err
454a238c70SJohn Marino with a bound for err */
464a238c70SJohn Marino mpfr_mul (xx, x, x, MPFR_RNDD); /* err <= 1 */
474a238c70SJohn Marino mpfr_ui_div (xx, 1, xx, MPFR_RNDU); /* upper bound for 1/(2x^2), err <= 2 */
484a238c70SJohn Marino mpfr_div_2ui (xx, xx, 1, MPFR_RNDU); /* exact */
494a238c70SJohn Marino mpfr_set_ui (t, 1, MPFR_RNDN); /* current term, exact */
504a238c70SJohn Marino mpfr_set (y, t, MPFR_RNDN); /* current sum */
514a238c70SJohn Marino mpfr_set_ui (err, 0, MPFR_RNDN);
524a238c70SJohn Marino for (k = 1; ; k++)
534a238c70SJohn Marino {
544a238c70SJohn Marino mpfr_mul_ui (t, t, 2 * k - 1, MPFR_RNDU); /* err <= 4k-3 */
554a238c70SJohn Marino mpfr_mul (t, t, xx, MPFR_RNDU); /* err <= 4k */
564a238c70SJohn Marino /* for -1 < x < 1, and |nx| < 1, we have |(1+x)^n| <= 1+7/4|nx|.
574a238c70SJohn Marino Indeed, for x>=0: log((1+x)^n) = n*log(1+x) <= n*x. Let y=n*x < 1,
584a238c70SJohn Marino then exp(y) <= 1+7/4*y.
594a238c70SJohn Marino For x<=0, let x=-x, we can prove by induction that (1-x)^n >= 1-n*x.*/
604a238c70SJohn Marino mpfr_mul_2si (err, err, MPFR_GET_EXP (y) - MPFR_GET_EXP (t), MPFR_RNDU);
614a238c70SJohn Marino mpfr_add_ui (err, err, 14 * k, MPFR_RNDU); /* 2^(1-p) * t <= 2 ulp(t) */
624a238c70SJohn Marino mpfr_div_2si (err, err, MPFR_GET_EXP (y) - MPFR_GET_EXP (t), MPFR_RNDU);
634a238c70SJohn Marino if (MPFR_GET_EXP (t) + (mpfr_exp_t) prec <= MPFR_GET_EXP (y))
644a238c70SJohn Marino {
654a238c70SJohn Marino /* the truncation error is bounded by |t| < ulp(y) */
664a238c70SJohn Marino mpfr_add_ui (err, err, 1, MPFR_RNDU);
674a238c70SJohn Marino break;
684a238c70SJohn Marino }
694a238c70SJohn Marino if (k & 1)
704a238c70SJohn Marino mpfr_sub (y, y, t, MPFR_RNDN);
714a238c70SJohn Marino else
724a238c70SJohn Marino mpfr_add (y, y, t, MPFR_RNDN);
734a238c70SJohn Marino }
744a238c70SJohn Marino /* the error on y is bounded by err*ulp(y) */
754a238c70SJohn Marino mpfr_mul (t, x, x, MPFR_RNDU); /* rel. err <= 2^(1-p) */
764a238c70SJohn Marino mpfr_div_2ui (err, err, 3, MPFR_RNDU); /* err/8 */
774a238c70SJohn Marino mpfr_add (err, err, t, MPFR_RNDU); /* err/8 + xx */
784a238c70SJohn Marino mpfr_mul_2ui (err, err, 3, MPFR_RNDU); /* err + 8*xx */
794a238c70SJohn Marino mpfr_exp (t, t, MPFR_RNDU); /* err <= 1/2*ulp(t) + err(x*x)*t
804a238c70SJohn Marino <= 1/2*ulp(t)+2*|x*x|*ulp(t)
814a238c70SJohn Marino <= (2*|x*x|+1/2)*ulp(t) */
824a238c70SJohn Marino mpfr_mul (t, t, x, MPFR_RNDN); /* err <= 1/2*ulp(t) + (4*|x*x|+1)*ulp(t)
834a238c70SJohn Marino <= (4*|x*x|+3/2)*ulp(t) */
844a238c70SJohn Marino mpfr_const_pi (xx, MPFR_RNDZ); /* err <= ulp(Pi) */
854a238c70SJohn Marino mpfr_sqrt (xx, xx, MPFR_RNDN); /* err <= 1/2*ulp(xx) + ulp(Pi)/2/sqrt(Pi)
864a238c70SJohn Marino <= 3/2*ulp(xx) */
874a238c70SJohn Marino mpfr_mul (t, t, xx, MPFR_RNDN); /* err <= (8 |xx| + 13/2) * ulp(t) */
884a238c70SJohn Marino mpfr_div (y, y, t, MPFR_RNDN); /* the relative error on input y is bounded
894a238c70SJohn Marino by (1+u)^err with u = 2^(1-p), that on
904a238c70SJohn Marino t is bounded by (1+u)^(8 |xx| + 13/2),
914a238c70SJohn Marino thus that on output y is bounded by
924a238c70SJohn Marino 8 |xx| + 7 + err. */
934a238c70SJohn Marino
944a238c70SJohn Marino if (MPFR_IS_ZERO(y))
954a238c70SJohn Marino {
964a238c70SJohn Marino /* If y is zero, most probably we have underflow. We check it directly
974a238c70SJohn Marino using the fact that erfc(x) <= exp(-x^2)/sqrt(Pi)/x for x >= 0.
984a238c70SJohn Marino We compute an upper approximation of exp(-x^2)/sqrt(Pi)/x.
994a238c70SJohn Marino */
1004a238c70SJohn Marino mpfr_mul (t, x, x, MPFR_RNDD); /* t <= x^2 */
1014a238c70SJohn Marino mpfr_neg (t, t, MPFR_RNDU); /* -x^2 <= t */
1024a238c70SJohn Marino mpfr_exp (t, t, MPFR_RNDU); /* exp(-x^2) <= t */
1034a238c70SJohn Marino mpfr_const_pi (xx, MPFR_RNDD); /* xx <= sqrt(Pi), cached */
1044a238c70SJohn Marino mpfr_mul (xx, xx, x, MPFR_RNDD); /* xx <= sqrt(Pi)*x */
1054a238c70SJohn Marino mpfr_div (y, t, xx, MPFR_RNDN); /* if y is zero, this means that the upper
1064a238c70SJohn Marino approximation of exp(-x^2)/sqrt(Pi)/x
1074a238c70SJohn Marino is nearer from 0 than from 2^(-emin-1),
1084a238c70SJohn Marino thus we have underflow. */
1094a238c70SJohn Marino exp_err = 0;
1104a238c70SJohn Marino }
1114a238c70SJohn Marino else
1124a238c70SJohn Marino {
1134a238c70SJohn Marino mpfr_add_ui (err, err, 7, MPFR_RNDU);
1144a238c70SJohn Marino exp_err = MPFR_GET_EXP (err);
1154a238c70SJohn Marino }
1164a238c70SJohn Marino
1174a238c70SJohn Marino mpfr_clear (t);
1184a238c70SJohn Marino mpfr_clear (xx);
1194a238c70SJohn Marino mpfr_clear (err);
1204a238c70SJohn Marino return exp_err;
1214a238c70SJohn Marino }
1224a238c70SJohn Marino
1234a238c70SJohn Marino int
mpfr_erfc(mpfr_ptr y,mpfr_srcptr x,mpfr_rnd_t rnd)1244a238c70SJohn Marino mpfr_erfc (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd)
1254a238c70SJohn Marino {
1264a238c70SJohn Marino int inex;
1274a238c70SJohn Marino mpfr_t tmp;
1284a238c70SJohn Marino mpfr_exp_t te, err;
1294a238c70SJohn Marino mpfr_prec_t prec;
1304a238c70SJohn Marino mpfr_exp_t emin = mpfr_get_emin ();
1314a238c70SJohn Marino MPFR_SAVE_EXPO_DECL (expo);
1324a238c70SJohn Marino MPFR_ZIV_DECL (loop);
1334a238c70SJohn Marino
1344a238c70SJohn Marino MPFR_LOG_FUNC
1354a238c70SJohn Marino (("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, rnd),
1364a238c70SJohn Marino ("y[%Pu]=%.*Rg inexact=%d", mpfr_get_prec (y), mpfr_log_prec, y, inex));
1374a238c70SJohn Marino
1384a238c70SJohn Marino if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
1394a238c70SJohn Marino {
1404a238c70SJohn Marino if (MPFR_IS_NAN (x))
1414a238c70SJohn Marino {
1424a238c70SJohn Marino MPFR_SET_NAN (y);
1434a238c70SJohn Marino MPFR_RET_NAN;
1444a238c70SJohn Marino }
1454a238c70SJohn Marino /* erfc(+inf) = 0+, erfc(-inf) = 2 erfc (0) = 1 */
1464a238c70SJohn Marino else if (MPFR_IS_INF (x))
1474a238c70SJohn Marino return mpfr_set_ui (y, MPFR_IS_POS (x) ? 0 : 2, rnd);
1484a238c70SJohn Marino else
1494a238c70SJohn Marino return mpfr_set_ui (y, 1, rnd);
1504a238c70SJohn Marino }
1514a238c70SJohn Marino
1524a238c70SJohn Marino if (MPFR_SIGN (x) > 0)
1534a238c70SJohn Marino {
1544a238c70SJohn Marino /* by default, emin = 1-2^30, thus the smallest representable
1554a238c70SJohn Marino number is 1/2*2^emin = 2^(-2^30):
1564a238c70SJohn Marino for x >= 27282, erfc(x) < 2^(-2^30-1), and
1574a238c70SJohn Marino for x >= 1787897414, erfc(x) < 2^(-2^62-1).
1584a238c70SJohn Marino */
1594a238c70SJohn Marino if ((emin >= -1073741823 && mpfr_cmp_ui (x, 27282) >= 0) ||
1604a238c70SJohn Marino mpfr_cmp_ui (x, 1787897414) >= 0)
1614a238c70SJohn Marino {
1624a238c70SJohn Marino /* May be incorrect if MPFR_EMAX_MAX >= 2^62. */
1634a238c70SJohn Marino MPFR_ASSERTN ((MPFR_EMAX_MAX >> 31) >> 31 == 0);
1644a238c70SJohn Marino return mpfr_underflow (y, (rnd == MPFR_RNDN) ? MPFR_RNDZ : rnd, 1);
1654a238c70SJohn Marino }
1664a238c70SJohn Marino }
1674a238c70SJohn Marino
1684a238c70SJohn Marino /* Init stuff */
1694a238c70SJohn Marino MPFR_SAVE_EXPO_MARK (expo);
1704a238c70SJohn Marino
1714a238c70SJohn Marino if (MPFR_SIGN (x) < 0)
1724a238c70SJohn Marino {
1734a238c70SJohn Marino mpfr_exp_t e = MPFR_EXP(x);
1744a238c70SJohn Marino /* For x < 0 going to -infinity, erfc(x) tends to 2 by below.
1754a238c70SJohn Marino More precisely, we have 2 + 1/sqrt(Pi)/x/exp(x^2) < erfc(x) < 2.
1764a238c70SJohn Marino Thus log2 |2 - erfc(x)| <= -log2|x| - x^2 / log(2).
1774a238c70SJohn Marino If |2 - erfc(x)| < 2^(-PREC(y)) then the result is either 2 or
1784a238c70SJohn Marino nextbelow(2).
1794a238c70SJohn Marino For x <= -27282, -log2|x| - x^2 / log(2) <= -2^30.
1804a238c70SJohn Marino */
1814a238c70SJohn Marino if ((MPFR_PREC(y) <= 7 && e >= 2) || /* x <= -2 */
1824a238c70SJohn Marino (MPFR_PREC(y) <= 25 && e >= 3) || /* x <= -4 */
1834a238c70SJohn Marino (MPFR_PREC(y) <= 120 && mpfr_cmp_si (x, -9) <= 0) ||
1844a238c70SJohn Marino mpfr_cmp_si (x, -27282) <= 0)
1854a238c70SJohn Marino {
1864a238c70SJohn Marino near_two:
1874a238c70SJohn Marino mpfr_set_ui (y, 2, MPFR_RNDN);
1884a238c70SJohn Marino mpfr_set_inexflag ();
1894a238c70SJohn Marino if (rnd == MPFR_RNDZ || rnd == MPFR_RNDD)
1904a238c70SJohn Marino {
1914a238c70SJohn Marino mpfr_nextbelow (y);
1924a238c70SJohn Marino inex = -1;
1934a238c70SJohn Marino }
1944a238c70SJohn Marino else
1954a238c70SJohn Marino inex = 1;
1964a238c70SJohn Marino goto end;
1974a238c70SJohn Marino }
1984a238c70SJohn Marino else if (e >= 3) /* more accurate test */
1994a238c70SJohn Marino {
2004a238c70SJohn Marino mpfr_t t, u;
2014a238c70SJohn Marino int near_2;
2024a238c70SJohn Marino mpfr_init2 (t, 32);
2034a238c70SJohn Marino mpfr_init2 (u, 32);
2044a238c70SJohn Marino /* the following is 1/log(2) rounded to zero on 32 bits */
2054a238c70SJohn Marino mpfr_set_str_binary (t, "1.0111000101010100011101100101001");
2064a238c70SJohn Marino mpfr_sqr (u, x, MPFR_RNDZ);
2074a238c70SJohn Marino mpfr_mul (t, t, u, MPFR_RNDZ); /* t <= x^2/log(2) */
2084a238c70SJohn Marino mpfr_neg (u, x, MPFR_RNDZ); /* 0 <= u <= |x| */
2094a238c70SJohn Marino mpfr_log2 (u, u, MPFR_RNDZ); /* u <= log2(|x|) */
2104a238c70SJohn Marino mpfr_add (t, t, u, MPFR_RNDZ); /* t <= log2|x| + x^2 / log(2) */
2114a238c70SJohn Marino /* Taking into account that mpfr_exp_t >= mpfr_prec_t */
2124a238c70SJohn Marino mpfr_set_exp_t (u, MPFR_PREC (y), MPFR_RNDU);
2134a238c70SJohn Marino near_2 = mpfr_cmp (t, u) >= 0; /* 1 if PREC(y) <= u <= t <= ... */
2144a238c70SJohn Marino mpfr_clear (t);
2154a238c70SJohn Marino mpfr_clear (u);
2164a238c70SJohn Marino if (near_2)
2174a238c70SJohn Marino goto near_two;
2184a238c70SJohn Marino }
2194a238c70SJohn Marino }
2204a238c70SJohn Marino
2214a238c70SJohn Marino /* erfc(x) ~ 1, with error < 2^(EXP(x)+1) */
2224a238c70SJohn Marino MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, __gmpfr_one, - MPFR_GET_EXP (x) - 1,
2234a238c70SJohn Marino 0, MPFR_SIGN(x) < 0,
2244a238c70SJohn Marino rnd, inex = _inexact; goto end);
2254a238c70SJohn Marino
2264a238c70SJohn Marino prec = MPFR_PREC (y) + MPFR_INT_CEIL_LOG2 (MPFR_PREC (y)) + 3;
2274a238c70SJohn Marino if (MPFR_GET_EXP (x) > 0)
2284a238c70SJohn Marino prec += 2 * MPFR_GET_EXP(x);
2294a238c70SJohn Marino
2304a238c70SJohn Marino mpfr_init2 (tmp, prec);
2314a238c70SJohn Marino
2324a238c70SJohn Marino MPFR_ZIV_INIT (loop, prec); /* Initialize the ZivLoop controler */
2334a238c70SJohn Marino for (;;) /* Infinite loop */
2344a238c70SJohn Marino {
2354a238c70SJohn Marino /* use asymptotic formula only whenever x^2 >= p*log(2),
2364a238c70SJohn Marino otherwise it will not converge */
2374a238c70SJohn Marino if (MPFR_SIGN (x) > 0 &&
2384a238c70SJohn Marino 2 * MPFR_GET_EXP (x) - 2 >= MPFR_INT_CEIL_LOG2 (prec))
2394a238c70SJohn Marino /* we have x^2 >= p in that case */
2404a238c70SJohn Marino {
2414a238c70SJohn Marino err = mpfr_erfc_asympt (tmp, x);
2424a238c70SJohn Marino if (err == 0) /* underflow case */
2434a238c70SJohn Marino {
2444a238c70SJohn Marino mpfr_clear (tmp);
2454a238c70SJohn Marino MPFR_SAVE_EXPO_FREE (expo);
2464a238c70SJohn Marino return mpfr_underflow (y, (rnd == MPFR_RNDN) ? MPFR_RNDZ : rnd, 1);
2474a238c70SJohn Marino }
2484a238c70SJohn Marino }
2494a238c70SJohn Marino else
2504a238c70SJohn Marino {
2514a238c70SJohn Marino mpfr_erf (tmp, x, MPFR_RNDN);
2524a238c70SJohn Marino MPFR_ASSERTD (!MPFR_IS_SINGULAR (tmp)); /* FIXME: 0 only for x=0 ? */
2534a238c70SJohn Marino te = MPFR_GET_EXP (tmp);
2544a238c70SJohn Marino mpfr_ui_sub (tmp, 1, tmp, MPFR_RNDN);
2554a238c70SJohn Marino /* See error analysis in algorithms.tex for details */
2564a238c70SJohn Marino if (MPFR_IS_ZERO (tmp))
2574a238c70SJohn Marino {
2584a238c70SJohn Marino prec *= 2;
2594a238c70SJohn Marino err = prec; /* ensures MPFR_CAN_ROUND fails */
2604a238c70SJohn Marino }
2614a238c70SJohn Marino else
2624a238c70SJohn Marino err = MAX (te - MPFR_GET_EXP (tmp), 0) + 1;
2634a238c70SJohn Marino }
2644a238c70SJohn Marino if (MPFR_LIKELY (MPFR_CAN_ROUND (tmp, prec - err, MPFR_PREC (y), rnd)))
2654a238c70SJohn Marino break;
2664a238c70SJohn Marino MPFR_ZIV_NEXT (loop, prec); /* Increase used precision */
2674a238c70SJohn Marino mpfr_set_prec (tmp, prec);
2684a238c70SJohn Marino }
2694a238c70SJohn Marino MPFR_ZIV_FREE (loop); /* Free the ZivLoop Controler */
2704a238c70SJohn Marino
2714a238c70SJohn Marino inex = mpfr_set (y, tmp, rnd); /* Set y to the computed value */
2724a238c70SJohn Marino mpfr_clear (tmp);
2734a238c70SJohn Marino
2744a238c70SJohn Marino end:
2754a238c70SJohn Marino MPFR_SAVE_EXPO_FREE (expo);
2764a238c70SJohn Marino return mpfr_check_range (y, inex, rnd);
2774a238c70SJohn Marino }
278