1 /* mpfr_erfc -- The Complementary Error Function of a floating-point number
2
3 Copyright 2005-2023 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramba 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 https://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_sqr (xx, 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_sqr (t, 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_sqr (t, 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[%Pd]=%.*Rg rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, rnd),
136 ("y[%Pd]=%.*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_IS_POS (x))
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_STAT_STATIC_ASSERT ((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_IS_NEG (x))
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_IS_NEG (x),
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 controller */
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_IS_POS (x) &&
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 Controller */
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 }