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 } 278