1 /* mpfr_erfc -- The Complementary Error Function of a floating-point number
3 Copyright 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramel projects, INRIA.
6 This file is part of the GNU MPFR Library.
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.
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.
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. */
23 #define MPFR_NEED_LONGLONG_H
24 #include "mpfr-impl.h"
26 /* erfc(x) = 1 - erf(x) */
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.
34 mpfr_erfc_asympt (mpfr_ptr y
, mpfr_srcptr x
)
38 mpfr_prec_t prec
= MPFR_PREC(y
);
42 mpfr_init2 (xx
, prec
);
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
);
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
))
65 /* the truncation error is bounded by |t| < ulp(y) */
66 mpfr_add_ui (err
, err
, 1, MPFR_RNDU
);
70 mpfr_sub (y
, y
, t
, MPFR_RNDN
);
72 mpfr_add (y
, y
, t
, MPFR_RNDN
);
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)
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
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.
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. */
113 mpfr_add_ui (err
, err
, 7, MPFR_RNDU
);
114 exp_err
= MPFR_GET_EXP (err
);
124 mpfr_erfc (mpfr_ptr y
, mpfr_srcptr x
, mpfr_rnd_t rnd
)
130 mpfr_exp_t emin
= mpfr_get_emin ();
131 MPFR_SAVE_EXPO_DECL (expo
);
132 MPFR_ZIV_DECL (loop
);
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
));
138 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x
)))
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
);
149 return mpfr_set_ui (y
, 1, rnd
);
152 if (MPFR_SIGN (x
) > 0)
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).
159 if ((emin
>= -1073741823 && mpfr_cmp_ui (x
, 27282) >= 0) ||
160 mpfr_cmp_ui (x
, 1787897414) >= 0)
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);
169 MPFR_SAVE_EXPO_MARK (expo
);
171 if (MPFR_SIGN (x
) < 0)
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
179 For x <= -27282, -log2|x| - x^2 / log(2) <= -2^30.
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)
187 mpfr_set_ui (y
, 2, MPFR_RNDN
);
188 mpfr_set_inexflag ();
189 if (rnd
== MPFR_RNDZ
|| rnd
== MPFR_RNDD
)
198 else if (e
>= 3) /* more accurate test */
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 <= ... */
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,
224 rnd
, inex
= _inexact
; goto end
);
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
);
230 mpfr_init2 (tmp
, prec
);
232 MPFR_ZIV_INIT (loop
, prec
); /* Initialize the ZivLoop controler */
233 for (;;) /* Infinite loop */
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 */
241 err
= mpfr_erfc_asympt (tmp
, x
);
242 if (err
== 0) /* underflow case */
245 MPFR_SAVE_EXPO_FREE (expo
);
246 return mpfr_underflow (y
, (rnd
== MPFR_RNDN
) ? MPFR_RNDZ
: rnd
, 1);
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
))
259 err
= prec
; /* ensures MPFR_CAN_ROUND fails */
262 err
= MAX (te
- MPFR_GET_EXP (tmp
), 0) + 1;
264 if (MPFR_LIKELY (MPFR_CAN_ROUND (tmp
, prec
- err
, MPFR_PREC (y
), rnd
)))
266 MPFR_ZIV_NEXT (loop
, prec
); /* Increase used precision */
267 mpfr_set_prec (tmp
, prec
);
269 MPFR_ZIV_FREE (loop
); /* Free the ZivLoop Controler */
271 inex
= mpfr_set (y
, tmp
, rnd
); /* Set y to the computed value */
275 MPFR_SAVE_EXPO_FREE (expo
);
276 return mpfr_check_range (y
, inex
, rnd
);