1 /* mpfr_erf -- error function of a floating-point number
3 Copyright 2001, 2003-2015 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 #define EXP1 2.71828182845904523536 /* exp(1) */
28 static int mpfr_erf_0 (mpfr_ptr
, mpfr_srcptr
, double, mpfr_rnd_t
);
31 mpfr_erf (mpfr_ptr y
, mpfr_srcptr x
, mpfr_rnd_t rnd_mode
)
35 MPFR_SAVE_EXPO_DECL (expo
);
38 (("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (x
), mpfr_log_prec
, x
, rnd_mode
),
39 ("y[%Pu]=%.*Rg inexact=%d", mpfr_get_prec (y
), mpfr_log_prec
, y
, inex
));
41 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x
)))
48 else if (MPFR_IS_INF (x
)) /* erf(+inf) = +1, erf(-inf) = -1 */
49 return mpfr_set_si (y
, MPFR_INT_SIGN (x
), MPFR_RNDN
);
50 else /* erf(+0) = +0, erf(-0) = -0 */
52 MPFR_ASSERTD (MPFR_IS_ZERO (x
));
53 return mpfr_set (y
, x
, MPFR_RNDN
); /* should keep the sign of x */
57 /* now x is neither NaN, Inf nor 0 */
59 /* first try expansion at x=0 when x is small, or asymptotic expansion
62 MPFR_SAVE_EXPO_MARK (expo
);
64 /* around x=0, we have erf(x) = 2x/sqrt(Pi) (1 - x^2/3 + ...),
65 with 1 - x^2/3 <= sqrt(Pi)*erf(x)/2/x <= 1 for x >= 0. This means that
66 if x^2/3 < 2^(-PREC(y)-1) we can decide of the correct rounding,
67 unless we have a worst-case for 2x/sqrt(Pi). */
68 if (MPFR_EXP(x
) < - (mpfr_exp_t
) (MPFR_PREC(y
) / 2))
70 /* we use 2x/sqrt(Pi) (1 - x^2/3) <= erf(x) <= 2x/sqrt(Pi) for x > 0
71 and 2x/sqrt(Pi) <= erf(x) <= 2x/sqrt(Pi) (1 - x^2/3) for x < 0.
72 In both cases |2x/sqrt(Pi) (1 - x^2/3)| <= |erf(x)| <= |2x/sqrt(Pi)|.
73 We will compute l and h such that l <= |2x/sqrt(Pi) (1 - x^2/3)|
74 and |2x/sqrt(Pi)| <= h. If l and h round to the same value to
75 precision PREC(y) and rounding rnd_mode, then we are done. */
76 mpfr_t l
, h
; /* lower and upper bounds for erf(x) */
79 mpfr_init2 (l
, MPFR_PREC(y
) + 17);
80 mpfr_init2 (h
, MPFR_PREC(y
) + 17);
82 mpfr_mul (l
, x
, x
, MPFR_RNDU
);
83 mpfr_div_ui (l
, l
, 3, MPFR_RNDU
); /* upper bound on x^2/3 */
84 mpfr_ui_sub (l
, 1, l
, MPFR_RNDZ
); /* lower bound on 1 - x^2/3 */
85 mpfr_const_pi (h
, MPFR_RNDU
); /* upper bound of Pi */
86 mpfr_sqrt (h
, h
, MPFR_RNDU
); /* upper bound on sqrt(Pi) */
87 mpfr_div (l
, l
, h
, MPFR_RNDZ
); /* lower bound on 1/sqrt(Pi) (1 - x^2/3) */
88 mpfr_mul_2ui (l
, l
, 1, MPFR_RNDZ
); /* 2/sqrt(Pi) (1 - x^2/3) */
89 mpfr_mul (l
, l
, x
, MPFR_RNDZ
); /* |l| is a lower bound on
90 |2x/sqrt(Pi) (1 - x^2/3)| */
92 mpfr_const_pi (h
, MPFR_RNDD
); /* lower bound on Pi */
93 mpfr_sqrt (h
, h
, MPFR_RNDD
); /* lower bound on sqrt(Pi) */
94 mpfr_div_2ui (h
, h
, 1, MPFR_RNDD
); /* lower bound on sqrt(Pi)/2 */
95 /* since sqrt(Pi)/2 < 1, the following should not underflow */
96 mpfr_div (h
, x
, h
, MPFR_IS_POS(x
) ? MPFR_RNDU
: MPFR_RNDD
);
97 /* round l and h to precision PREC(y) */
98 inex
= mpfr_prec_round (l
, MPFR_PREC(y
), rnd_mode
);
99 inex2
= mpfr_prec_round (h
, MPFR_PREC(y
), rnd_mode
);
100 /* Caution: we also need inex=inex2 (inex might be 0). */
101 ok
= SAME_SIGN (inex
, inex2
) && mpfr_cmp (l
, h
) == 0;
103 mpfr_set (y
, h
, rnd_mode
);
108 /* this test can still fail for small precision, for example
109 for x=-0.100E-2 with a target precision of 3 bits, since
110 the error term x^2/3 is not that small. */
114 mpfr_const_log2 (xf
, MPFR_RNDU
);
115 mpfr_div (xf
, x
, xf
, MPFR_RNDZ
); /* round to zero ensures we get a lower
116 bound of |x/log(2)| */
117 mpfr_mul (xf
, xf
, x
, MPFR_RNDZ
);
118 large
= mpfr_cmp_ui (xf
, MPFR_PREC (y
) + 1) > 0;
121 /* when x goes to infinity, we have erf(x) = 1 - 1/sqrt(Pi)/exp(x^2)/x + ...
122 and |erf(x) - 1| <= exp(-x^2) is true for any x >= 0, thus if
123 exp(-x^2) < 2^(-PREC(y)-1) the result is 1 or 1-epsilon.
124 This rewrites as x^2/log(2) > p+1. */
125 if (MPFR_UNLIKELY (large
))
126 /* |erf x| = 1 or 1- */
128 mpfr_rnd_t rnd2
= MPFR_IS_POS (x
) ? rnd_mode
: MPFR_INVERT_RND(rnd_mode
);
129 if (rnd2
== MPFR_RNDN
|| rnd2
== MPFR_RNDU
|| rnd2
== MPFR_RNDA
)
131 inex
= MPFR_INT_SIGN (x
);
132 mpfr_set_si (y
, inex
, rnd2
);
134 else /* round to zero */
136 inex
= -MPFR_INT_SIGN (x
);
137 mpfr_setmax (y
, 0); /* warning: setmax keeps the old sign of y */
138 MPFR_SET_SAME_SIGN (y
, x
);
141 else /* use Taylor */
145 /* FIXME: get rid of doubles/mpfr_get_d here */
146 xf2
= mpfr_get_d (x
, MPFR_RNDN
);
147 xf2
= xf2
* xf2
; /* xf2 ~ x^2 */
148 inex
= mpfr_erf_0 (y
, x
, xf2
, rnd_mode
);
152 MPFR_SAVE_EXPO_FREE (expo
);
153 return mpfr_check_range (y
, inex
, rnd_mode
);
158 mul_2exp (double x
, mpfr_exp_t e
)
174 /* evaluates erf(x) using the expansion at x=0:
176 erf(x) = 2/sqrt(Pi) * sum((-1)^k*x^(2k+1)/k!/(2k+1), k=0..infinity)
178 Assumes x is neither NaN nor infinite nor zero.
179 Assumes also that e*x^2 <= n (target precision).
182 mpfr_erf_0 (mpfr_ptr res
, mpfr_srcptr x
, double xf2
, mpfr_rnd_t rnd_mode
)
185 mpfr_exp_t nuk
, sigmak
;
191 MPFR_ZIV_DECL (loop
);
193 n
= MPFR_PREC (res
); /* target precision */
195 /* initial working precision */
196 m
= n
+ (mpfr_prec_t
) (xf2
/ LOG2
) + 8 + MPFR_INT_CEIL_LOG2 (n
);
203 MPFR_ZIV_INIT (loop
, m
);
206 mpfr_mul (y
, x
, x
, MPFR_RNDU
); /* err <= 1 ulp */
207 mpfr_set_ui (s
, 1, MPFR_RNDN
);
208 mpfr_set_ui (t
, 1, MPFR_RNDN
);
213 mpfr_mul (t
, y
, t
, MPFR_RNDU
);
214 mpfr_div_ui (t
, t
, k
, MPFR_RNDU
);
215 mpfr_div_ui (u
, t
, 2 * k
+ 1, MPFR_RNDU
);
216 sigmak
= MPFR_GET_EXP (s
);
218 mpfr_sub (s
, s
, u
, MPFR_RNDN
);
220 mpfr_add (s
, s
, u
, MPFR_RNDN
);
221 sigmak
-= MPFR_GET_EXP(s
);
222 nuk
= MPFR_GET_EXP(u
) - MPFR_GET_EXP(s
);
224 if ((nuk
< - (mpfr_exp_t
) m
) && ((double) k
>= xf2
))
227 /* tauk <- 1/2 + tauk * 2^sigmak + (1+8k)*2^nuk */
228 tauk
= 0.5 + mul_2exp (tauk
, sigmak
)
229 + mul_2exp (1.0 + 8.0 * (double) k
, nuk
);
232 mpfr_mul (s
, x
, s
, MPFR_RNDU
);
233 MPFR_SET_EXP (s
, MPFR_GET_EXP (s
) + 1);
235 mpfr_const_pi (t
, MPFR_RNDZ
);
236 mpfr_sqrt (t
, t
, MPFR_RNDZ
);
237 mpfr_div (s
, s
, t
, MPFR_RNDN
);
238 tauk
= 4.0 * tauk
+ 11.0; /* final ulp-error on s */
239 log2tauk
= __gmpfr_ceil_log2 (tauk
);
241 if (MPFR_LIKELY (MPFR_CAN_ROUND (s
, m
- log2tauk
, n
, rnd_mode
)))
244 /* Actualisation of the precision */
245 MPFR_ZIV_NEXT (loop
, m
);
246 mpfr_set_prec (y
, m
);
247 mpfr_set_prec (s
, m
);
248 mpfr_set_prec (t
, m
);
249 mpfr_set_prec (u
, m
);
252 MPFR_ZIV_FREE (loop
);
254 inex
= mpfr_set (res
, s
, rnd_mode
);