1 /* mpfr_erf -- error function of a floating-point number
2
3 Copyright 2001, 2003-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 static int mpfr_erf_0 (mpfr_ptr, mpfr_srcptr, double, mpfr_rnd_t);
27
28 int
29 mpfr_erf (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
30 {
31 mpfr_t xf;
32 mp_limb_t xf_limb[(53 - 1) / GMP_NUMB_BITS + 1];
33 int inex, large;
34 MPFR_SAVE_EXPO_DECL (expo);
35
36 MPFR_LOG_FUNC
37 (("x[%Pd]=%.*Rg rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, rnd_mode),
38 ("y[%Pd]=%.*Rg inexact=%d", mpfr_get_prec (y), mpfr_log_prec, y, inex));
39
40 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
41 {
42 if (MPFR_IS_NAN (x))
43 {
44 MPFR_SET_NAN (y);
45 MPFR_RET_NAN;
46 }
47 else if (MPFR_IS_INF (x)) /* erf(+inf) = +1, erf(-inf) = -1 */
48 return mpfr_set_si (y, MPFR_INT_SIGN (x), MPFR_RNDN);
49 else /* erf(+0) = +0, erf(-0) = -0 */
50 {
51 MPFR_ASSERTD (MPFR_IS_ZERO (x));
52 return mpfr_set (y, x, MPFR_RNDN); /* should keep the sign of x */
53 }
54 }
55
56 /* now x is neither NaN, Inf nor 0 */
57
58 /* first try expansion at x=0 when x is small, or asymptotic expansion
59 where x is large */
60
61 MPFR_SAVE_EXPO_MARK (expo);
62
63 /* around x=0, we have erf(x) = 2x/sqrt(Pi) (1 - x^2/3 + ...),
64 with 1 - x^2/3 <= sqrt(Pi)*erf(x)/2/x <= 1 for x >= 0. This means that
65 if x^2/3 < 2^(-PREC(y)-1) we can decide of the correct rounding,
66 unless we have a worst-case for 2x/sqrt(Pi). */
67 if (MPFR_EXP(x) < - (mpfr_exp_t) (MPFR_PREC(y) / 2))
68 {
69 /* we use 2x/sqrt(Pi) (1 - x^2/3) <= erf(x) <= 2x/sqrt(Pi) for x > 0
70 and 2x/sqrt(Pi) <= erf(x) <= 2x/sqrt(Pi) (1 - x^2/3) for x < 0.
71 In both cases |2x/sqrt(Pi) (1 - x^2/3)| <= |erf(x)| <= |2x/sqrt(Pi)|.
72 We will compute l and h such that l <= |2x/sqrt(Pi) (1 - x^2/3)|
73 and |2x/sqrt(Pi)| <= h. If l and h round to the same value to
74 precision PREC(y) and rounding rnd_mode, then we are done. */
75 mpfr_t l, h; /* lower and upper bounds for erf(x) */
76 int ok, inex2;
77
78 mpfr_init2 (l, MPFR_PREC(y) + 17);
79 mpfr_init2 (h, MPFR_PREC(y) + 17);
80 /* first compute l */
81 mpfr_sqr (l, x, MPFR_RNDU);
82 mpfr_div_ui (l, l, 3, MPFR_RNDU); /* upper bound on x^2/3 */
83 mpfr_ui_sub (l, 1, l, MPFR_RNDZ); /* lower bound on 1 - x^2/3 */
84 mpfr_const_pi (h, MPFR_RNDU); /* upper bound of Pi */
85 mpfr_sqrt (h, h, MPFR_RNDU); /* upper bound on sqrt(Pi) */
86 mpfr_div (l, l, h, MPFR_RNDZ); /* lower bound on 1/sqrt(Pi) (1 - x^2/3) */
87 mpfr_mul_2ui (l, l, 1, MPFR_RNDZ); /* 2/sqrt(Pi) (1 - x^2/3) */
88 mpfr_mul (l, l, x, MPFR_RNDZ); /* |l| is a lower bound on
89 |2x/sqrt(Pi) (1 - x^2/3)| */
90 /* now compute h */
91 mpfr_const_pi (h, MPFR_RNDD); /* lower bound on Pi */
92 mpfr_sqrt (h, h, MPFR_RNDD); /* lower bound on sqrt(Pi) */
93 mpfr_div_2ui (h, h, 1, MPFR_RNDD); /* lower bound on sqrt(Pi)/2 */
94 /* since sqrt(Pi)/2 < 1, the following should not underflow */
95 mpfr_div (h, x, h, MPFR_IS_POS(x) ? MPFR_RNDU : MPFR_RNDD);
96 /* round l and h to precision PREC(y) */
97 inex = mpfr_prec_round (l, MPFR_PREC(y), rnd_mode);
98 inex2 = mpfr_prec_round (h, MPFR_PREC(y), rnd_mode);
99 /* Caution: we also need inex=inex2 (inex might be 0). */
100 ok = SAME_SIGN (inex, inex2) && mpfr_equal_p (l, h);
101 if (ok)
102 mpfr_set (y, h, rnd_mode);
103 mpfr_clear (l);
104 mpfr_clear (h);
105 if (ok)
106 goto end;
107 /* this test can still fail for small precision, for example
108 for x=-0.100E-2 with a target precision of 3 bits, since
109 the error term x^2/3 is not that small. */
110 }
111
112 MPFR_TMP_INIT1(xf_limb, xf, 53);
113 mpfr_div (xf, x, __gmpfr_const_log2_RNDU, MPFR_RNDZ); /* round to zero
114 ensures we get a lower bound of |x/log(2)| */
115 mpfr_mul (xf, xf, x, MPFR_RNDZ);
116 large = mpfr_cmp_ui (xf, MPFR_PREC (y) + 1) > 0;
117
118 /* when x goes to infinity, we have erf(x) = 1 - 1/sqrt(Pi)/exp(x^2)/x + ...
119 and |erf(x) - 1| <= exp(-x^2) is true for any x >= 0, thus if
120 exp(-x^2) < 2^(-PREC(y)-1) the result is 1 or 1-epsilon.
121 This rewrites as x^2/log(2) > p+1. */
122 if (MPFR_UNLIKELY (large))
123 /* |erf x| = 1 or 1- */
124 {
125 mpfr_rnd_t rnd2 = MPFR_IS_POS (x) ? rnd_mode : MPFR_INVERT_RND(rnd_mode);
126 if (rnd2 == MPFR_RNDN || rnd2 == MPFR_RNDU || rnd2 == MPFR_RNDA)
127 {
128 inex = MPFR_INT_SIGN (x);
129 mpfr_set_si (y, inex, rnd2);
130 }
131 else /* round to zero */
132 {
133 inex = -MPFR_INT_SIGN (x);
134 mpfr_setmax (y, 0); /* warning: setmax keeps the old sign of y */
135 MPFR_SET_SAME_SIGN (y, x);
136 }
137 }
138 else /* use Taylor */
139 {
140 double xf2;
141
142 /* FIXME: get rid of doubles/mpfr_get_d here */
143 xf2 = mpfr_get_d (x, MPFR_RNDN);
144 xf2 = xf2 * xf2; /* xf2 ~ x^2 */
145 inex = mpfr_erf_0 (y, x, xf2, rnd_mode);
146 }
147
148 end:
149 MPFR_SAVE_EXPO_FREE (expo);
150 return mpfr_check_range (y, inex, rnd_mode);
151 }
152
153 /* return x*2^e */
154 static double
155 mul_2exp (double x, mpfr_exp_t e)
156 {
157 /* Most of the times, the argument is negative */
158 if (MPFR_UNLIKELY (e > 0))
159 {
160 while (e--)
161 x *= 2.0;
162 }
163 else
164 {
165 while (e <= -16)
166 {
167 x *= (1.0 / 65536.0);
168 e += 16;
169 }
170 while (e++)
171 x *= 0.5;
172 }
173
174 return x;
175 }
176
177 /* evaluates erf(x) using the expansion at x=0:
178
179 erf(x) = 2/sqrt(Pi) * sum((-1)^k*x^(2k+1)/k!/(2k+1), k=0..infinity)
180
181 Assumes x is neither NaN nor infinite nor zero.
182 Assumes also that e*x^2 <= n (target precision).
183 */
184 static int
185 mpfr_erf_0 (mpfr_ptr res, mpfr_srcptr x, double xf2, mpfr_rnd_t rnd_mode)
186 {
187 mpfr_prec_t n, m;
188 mpfr_exp_t nuk, sigmak;
189 mpfr_t y, s, t, u;
190 unsigned int k;
191 int inex;
192 MPFR_GROUP_DECL (group);
193 MPFR_ZIV_DECL (loop);
194
195 n = MPFR_PREC (res); /* target precision */
196
197 /* initial working precision */
198 m = n + (mpfr_prec_t) (xf2 / LOG2) + 8 + MPFR_INT_CEIL_LOG2 (n);
199
200 MPFR_GROUP_INIT_4(group, m, y, s, t, u);
201
202 MPFR_ZIV_INIT (loop, m);
203 for (;;)
204 {
205 mpfr_t tauk;
206 mpfr_exp_t log2tauk;
207
208 mpfr_sqr (y, x, MPFR_RNDU); /* err <= 1 ulp */
209 mpfr_set_ui (s, 1, MPFR_RNDN);
210 mpfr_set_ui (t, 1, MPFR_RNDN);
211 mpfr_init2 (tauk, 53);
212 mpfr_set_ui (tauk, 0, MPFR_RNDU);
213
214 for (k = 1; ; k++)
215 {
216 mpfr_mul (t, y, t, MPFR_RNDU);
217 mpfr_div_ui (t, t, k, MPFR_RNDU);
218 mpfr_div_ui (u, t, 2 * k + 1, MPFR_RNDU);
219 sigmak = MPFR_GET_EXP (s);
220 if (k % 2)
221 mpfr_sub (s, s, u, MPFR_RNDN);
222 else
223 mpfr_add (s, s, u, MPFR_RNDN);
224 sigmak -= MPFR_GET_EXP(s);
225 nuk = MPFR_GET_EXP(u) - MPFR_GET_EXP(s);
226
227 if ((nuk < - (mpfr_exp_t) m) && (k >= xf2))
228 break;
229
230 /* tauk <- 1/2 + tauk * 2^sigmak + (1+8k)*2^nuk */
231 mpfr_mul_2si (tauk, tauk, sigmak, MPFR_RNDU);
232 mpfr_add_d (tauk, tauk, 0.5 + mul_2exp (1.0 + 8.0 * (double) k, nuk),
233 MPFR_RNDU);
234 }
235
236 mpfr_mul (s, x, s, MPFR_RNDU);
237 MPFR_SET_EXP (s, MPFR_GET_EXP (s) + 1);
238
239 mpfr_const_pi (t, MPFR_RNDZ);
240 mpfr_sqrt (t, t, MPFR_RNDZ);
241 mpfr_div (s, s, t, MPFR_RNDN);
242 /* tauk = 4 * tauk + 11: final ulp-error on s */
243 mpfr_mul_2ui (tauk, tauk, 2, MPFR_RNDU);
244 mpfr_add_ui (tauk, tauk, 11, MPFR_RNDU);
245 /* In practice, one should not get an infinity, as it would require
246 a huge precision and a lot of time, but just in case... */
247 MPFR_ASSERTN (!MPFR_IS_INF (tauk));
248 log2tauk = MPFR_GET_EXP (tauk);
249 mpfr_clear (tauk);
250
251 if (MPFR_LIKELY (MPFR_CAN_ROUND (s, m - log2tauk, n, rnd_mode)))
252 break;
253
254 /* Actualisation of the precision */
255 MPFR_ZIV_NEXT (loop, m);
256 MPFR_GROUP_REPREC_4 (group, m, y, s, t, u);
257 }
258 MPFR_ZIV_FREE (loop);
259
260 inex = mpfr_set (res, s, rnd_mode);
261
262 MPFR_GROUP_CLEAR (group);
263
264 return inex;
265 }