1/* mpfr_erfc -- The Complementary Error Function of a floating-point number 2 3Copyright 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc. 4Contributed by the AriC and Caramel projects, INRIA. 5 6This file is part of the GNU MPFR Library. 7 8The GNU MPFR Library is free software; you can redistribute it and/or modify 9it under the terms of the GNU Lesser General Public License as published by 10the Free Software Foundation; either version 3 of the License, or (at your 11option) any later version. 12 13The GNU MPFR Library is distributed in the hope that it will be useful, but 14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16License for more details. 17 18You should have received a copy of the GNU Lesser General Public License 19along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 20http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2151 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*/ 33static mpfr_exp_t 34mpfr_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 123int 124mpfr_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} 278