2 * ====================================================
3 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 * Developed at SunPro, a Sun Microsystems, Inc. business.
6 * Permission to use, copy, modify, and distribute this
7 * software is freely granted, provided that this notice
9 * ====================================================
13 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
15 * Permission to use, copy, modify, and distribute this software for any
16 * purpose with or without fee is hereby granted, provided that the above
17 * copyright notice and this permission notice appear in all copies.
19 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
20 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
21 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
22 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
23 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
24 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
25 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
28 /* double erf(double x)
29 * double erfc(double x)
32 * erf(x) = --------- | exp(-t*t)dt
39 * erfc(-x) = 2 - erfc(x)
42 * 1. For |x| in [0, 0.84375]
43 * erf(x) = x + x*R(x^2)
44 * erfc(x) = 1 - erf(x) if x in [-.84375,0.25]
45 * = 0.5 + ((0.5-x)-x*R) if x in [0.25,0.84375]
46 * Remark. The formula is derived by noting
47 * erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....)
49 * 2/sqrt(pi) = 1.128379167095512573896158903121545171688
50 * is close to one. The interval is chosen because the fix
51 * point of erf(x) is near 0.6174 (i.e., erf(x)=x when x is
52 * near 0.6174), and by some experiment, 0.84375 is chosen to
53 * guarantee the error is less than one ulp for erf.
55 * 2. For |x| in [0.84375,1.25], let s = |x| - 1, and
56 * c = 0.84506291151 rounded to single (24 bits)
57 * erf(x) = sign(x) * (c + P1(s)/Q1(s))
58 * erfc(x) = (1-c) - P1(s)/Q1(s) if x > 0
59 * 1+(c+P1(s)/Q1(s)) if x < 0
60 * Remark: here we use the taylor series expansion at x=1.
61 * erf(1+s) = erf(1) + s*Poly(s)
62 * = 0.845.. + P1(s)/Q1(s)
63 * Note that |P1/Q1|< 0.078 for x in [0.84375,1.25]
65 * 3. For x in [1.25,1/0.35(~2.857143)],
66 * erfc(x) = (1/x)*exp(-x*x-0.5625+R1(z)/S1(z))
68 * erf(x) = 1 - erfc(x)
70 * 4. For x in [1/0.35,107]
71 * erfc(x) = (1/x)*exp(-x*x-0.5625+R2/S2) if x > 0
72 * = 2.0 - (1/x)*exp(-x*x-0.5625+R2(z)/S2(z))
74 * = 2.0 - tiny (if x <= -6.666)
76 * erf(x) = sign(x)*(1.0 - erfc(x)) if x < 6.666, else
77 * erf(x) = sign(x)*(1.0 - tiny)
79 * To compute exp(-x*x-0.5625+R/S), let s be a single
80 * precision number and s := x; then
81 * -x*x = -s*s + (s-x)*(s+x)
82 * exp(-x*x-0.5626+R/S) =
83 * exp(-s*s-0.5625)*exp((s-x)*(s+x)+R/S);
85 * Here 4 and 5 make use of the asymptotic series
87 * erfc(x) ~ ---------- * ( 1 + Poly(1/x^2) )
90 * 5. For inf > x >= 107
91 * erf(x) = sign(x) *(1 - tiny) (raise inexact)
92 * erfc(x) = tiny*tiny (raise underflow) if x > 0
96 * erf(0) = 0, erf(inf) = 1, erf(-inf) = -1,
97 * erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2,
98 * erfc/erf(NaN) is NaN
104 #include "math_private.h"
106 static const long double
111 /* c = (float)0.84506291151 */
112 erx
= 0.845062911510467529296875L,
114 * Coefficients for approximation to erf on [0,0.84375]
117 efx
= 1.2837916709551257389615890312154517168810E-1L,
118 /* 8 * (2/sqrt(pi) - 1) */
119 efx8
= 1.0270333367641005911692712249723613735048E0L
,
122 1.122751350964552113068262337278335028553E6L
,
123 -2.808533301997696164408397079650699163276E6L
,
124 -3.314325479115357458197119660818768924100E5L
,
125 -6.848684465326256109712135497895525446398E4L
,
126 -2.657817695110739185591505062971929859314E3L
,
127 -1.655310302737837556654146291646499062882E2L
,
131 8.745588372054466262548908189000448124232E6L
,
132 3.746038264792471129367533128637019611485E6L
,
133 7.066358783162407559861156173539693900031E5L
,
134 7.448928604824620999413120955705448117056E4L
,
135 4.511583986730994111992253980546131408924E3L
,
136 1.368902937933296323345610240009071254014E2L
,
137 /* 1.000000000000000000000000000000000000000E0 */
141 * Coefficients for approximation to erf in [0.84375,1.25]
143 /* erf(x+1) = 0.845062911510467529296875 + pa(x)/qa(x)
144 -0.15625 <= x <= +.25
145 Peak relative error 8.5e-22 */
148 -1.076952146179812072156734957705102256059E0L
,
149 1.884814957770385593365179835059971587220E2L
,
150 -5.339153975012804282890066622962070115606E1L
,
151 4.435910679869176625928504532109635632618E1L
,
152 1.683219516032328828278557309642929135179E1L
,
153 -2.360236618396952560064259585299045804293E0L
,
154 1.852230047861891953244413872297940938041E0L
,
155 9.394994446747752308256773044667843200719E-2L,
159 4.559263722294508998149925774781887811255E2L
,
160 3.289248982200800575749795055149780689738E2L
,
161 2.846070965875643009598627918383314457912E2L
,
162 1.398715859064535039433275722017479994465E2L
,
163 6.060190733759793706299079050985358190726E1L
,
164 2.078695677795422351040502569964299664233E1L
,
165 4.641271134150895940966798357442234498546E0L
,
166 /* 1.000000000000000000000000000000000000000E0 */
170 * Coefficients for approximation to erfc in [1.25,1/0.35]
172 /* erfc(1/x) = x exp (-1/x^2 - 0.5625 + ra(x^2)/sa(x^2))
173 1/2.85711669921875 < 1/x < 1/1.25
174 Peak relative error 3.1e-21 */
177 1.363566591833846324191000679620738857234E-1L,
178 1.018203167219873573808450274314658434507E1L
,
179 1.862359362334248675526472871224778045594E2L
,
180 1.411622588180721285284945138667933330348E3L
,
181 5.088538459741511988784440103218342840478E3L
,
182 8.928251553922176506858267311750789273656E3L
,
183 7.264436000148052545243018622742770549982E3L
,
184 2.387492459664548651671894725748959751119E3L
,
185 2.220916652813908085449221282808458466556E2L
,
189 -1.382234625202480685182526402169222331847E1L
,
190 -3.315638835627950255832519203687435946482E2L
,
191 -2.949124863912936259747237164260785326692E3L
,
192 -1.246622099070875940506391433635999693661E4L
,
193 -2.673079795851665428695842853070996219632E4L
,
194 -2.880269786660559337358397106518918220991E4L
,
195 -1.450600228493968044773354186390390823713E4L
,
196 -2.874539731125893533960680525192064277816E3L
,
197 -1.402241261419067750237395034116942296027E2L
,
198 /* 1.000000000000000000000000000000000000000E0 */
201 * Coefficients for approximation to erfc in [1/.35,107]
203 /* erfc(1/x) = x exp (-1/x^2 - 0.5625 + rb(x^2)/sb(x^2))
204 1/6.6666259765625 < 1/x < 1/2.85711669921875
205 Peak relative error 4.2e-22 */
207 -4.869587348270494309550558460786501252369E-5L,
208 -4.030199390527997378549161722412466959403E-3L,
209 -9.434425866377037610206443566288917589122E-2L,
210 -9.319032754357658601200655161585539404155E-1L,
211 -4.273788174307459947350256581445442062291E0L
,
212 -8.842289940696150508373541814064198259278E0L
,
213 -7.069215249419887403187988144752613025255E0L
,
214 -1.401228723639514787920274427443330704764E0L
,
218 4.936254964107175160157544545879293019085E-3L,
219 1.583457624037795744377163924895349412015E-1L,
220 1.850647991850328356622940552450636420484E0L
,
221 9.927611557279019463768050710008450625415E0L
,
222 2.531667257649436709617165336779212114570E1L
,
223 2.869752886406743386458304052862814690045E1L
,
224 1.182059497870819562441683560749192539345E1L
,
225 /* 1.000000000000000000000000000000000000000E0 */
227 /* erfc(1/x) = x exp (-1/x^2 - 0.5625 + rc(x^2)/sc(x^2))
228 1/107 <= 1/x <= 1/6.6666259765625
229 Peak relative error 1.1e-21 */
231 -8.299617545269701963973537248996670806850E-5L,
232 -6.243845685115818513578933902532056244108E-3L,
233 -1.141667210620380223113693474478394397230E-1L,
234 -7.521343797212024245375240432734425789409E-1L,
235 -1.765321928311155824664963633786967602934E0L
,
236 -1.029403473103215800456761180695263439188E0L
,
240 8.413244363014929493035952542677768808601E-3L,
241 2.065114333816877479753334599639158060979E-1L,
242 1.639064941530797583766364412782135680148E0L
,
243 4.936788463787115555582319302981666347450E0L
,
244 5.005177727208955487404729933261347679090E0L
,
245 /* 1.000000000000000000000000000000000000000E0 */
251 long double R
, S
, P
, Q
, s
, y
, z
, r
;
253 u_int32_t se
, i0
, i1
;
255 GET_LDOUBLE_WORDS (se
, i0
, i1
, x
);
260 i
= ((se
& 0xffff) >> 15) << 1;
261 return (long double) (1 - i
) + one
/ x
; /* erf(+-inf)=+-1 */
264 ix
= (ix
<< 16) | (i0
>> 16);
265 if (ix
< 0x3ffed800) /* |x|<0.84375 */
267 if (ix
< 0x3fde8000) /* |x|<2**-33 */
270 return 0.125 * (8.0 * x
+ efx8
* x
); /*avoid underflow */
274 r
= pp
[0] + z
* (pp
[1]
275 + z
* (pp
[2] + z
* (pp
[3] + z
* (pp
[4] + z
* pp
[5]))));
276 s
= qq
[0] + z
* (qq
[1]
277 + z
* (qq
[2] + z
* (qq
[3] + z
* (qq
[4] + z
* (qq
[5] + z
)))));
281 if (ix
< 0x3fffa000) /* 1.25 */
282 { /* 0.84375 <= |x| < 1.25 */
284 P
= pa
[0] + s
* (pa
[1] + s
* (pa
[2]
285 + s
* (pa
[3] + s
* (pa
[4] + s
* (pa
[5] + s
* (pa
[6] + s
* pa
[7]))))));
286 Q
= qa
[0] + s
* (qa
[1] + s
* (qa
[2]
287 + s
* (qa
[3] + s
* (qa
[4] + s
* (qa
[5] + s
* (qa
[6] + s
))))));
288 if ((se
& 0x8000) == 0)
293 if (ix
>= 0x4001d555) /* 6.6666259765625 */
294 { /* inf>|x|>=6.666 */
295 if ((se
& 0x8000) == 0)
302 if (ix
< 0x4000b6db) /* 2.85711669921875 */
304 R
= ra
[0] + s
* (ra
[1] + s
* (ra
[2] + s
* (ra
[3] + s
* (ra
[4] +
305 s
* (ra
[5] + s
* (ra
[6] + s
* (ra
[7] + s
* ra
[8])))))));
306 S
= sa
[0] + s
* (sa
[1] + s
* (sa
[2] + s
* (sa
[3] + s
* (sa
[4] +
307 s
* (sa
[5] + s
* (sa
[6] + s
* (sa
[7] + s
* (sa
[8] + s
))))))));
310 { /* |x| >= 1/0.35 */
311 R
= rb
[0] + s
* (rb
[1] + s
* (rb
[2] + s
* (rb
[3] + s
* (rb
[4] +
312 s
* (rb
[5] + s
* (rb
[6] + s
* rb
[7]))))));
313 S
= sb
[0] + s
* (sb
[1] + s
* (sb
[2] + s
* (sb
[3] + s
* (sb
[4] +
314 s
* (sb
[5] + s
* (sb
[6] + s
))))));
317 GET_LDOUBLE_WORDS (i
, i0
, i1
, z
);
319 SET_LDOUBLE_WORDS (z
, i
, i0
, i1
);
321 expl (-z
* z
- 0.5625) * expl ((z
- x
) * (z
+ x
) + R
/ S
);
322 if ((se
& 0x8000) == 0)
332 long double R
, S
, P
, Q
, s
, y
, z
, r
;
333 u_int32_t se
, i0
, i1
;
335 GET_LDOUBLE_WORDS (se
, i0
, i1
, x
);
338 { /* erfc(nan)=nan */
339 /* erfc(+-inf)=0,2 */
340 return (long double) (((se
& 0xffff) >> 15) << 1) + one
/ x
;
343 ix
= (ix
<< 16) | (i0
>> 16);
344 if (ix
< 0x3ffed800) /* |x|<0.84375 */
346 if (ix
< 0x3fbe0000) /* |x|<2**-65 */
349 r
= pp
[0] + z
* (pp
[1]
350 + z
* (pp
[2] + z
* (pp
[3] + z
* (pp
[4] + z
* pp
[5]))));
351 s
= qq
[0] + z
* (qq
[1]
352 + z
* (qq
[2] + z
* (qq
[3] + z
* (qq
[4] + z
* (qq
[5] + z
)))));
354 if (ix
< 0x3ffd8000) /* x<1/4 */
356 return one
- (x
+ x
* y
);
365 if (ix
< 0x3fffa000) /* 1.25 */
366 { /* 0.84375 <= |x| < 1.25 */
368 P
= pa
[0] + s
* (pa
[1] + s
* (pa
[2]
369 + s
* (pa
[3] + s
* (pa
[4] + s
* (pa
[5] + s
* (pa
[6] + s
* pa
[7]))))));
370 Q
= qa
[0] + s
* (qa
[1] + s
* (qa
[2]
371 + s
* (qa
[3] + s
* (qa
[4] + s
* (qa
[5] + s
* (qa
[6] + s
))))));
372 if ((se
& 0x8000) == 0)
383 if (ix
< 0x4005d600) /* 107 */
387 if (ix
< 0x4000b6db) /* 2.85711669921875 */
388 { /* |x| < 1/.35 ~ 2.857143 */
389 R
= ra
[0] + s
* (ra
[1] + s
* (ra
[2] + s
* (ra
[3] + s
* (ra
[4] +
390 s
* (ra
[5] + s
* (ra
[6] + s
* (ra
[7] + s
* ra
[8])))))));
391 S
= sa
[0] + s
* (sa
[1] + s
* (sa
[2] + s
* (sa
[3] + s
* (sa
[4] +
392 s
* (sa
[5] + s
* (sa
[6] + s
* (sa
[7] + s
* (sa
[8] + s
))))))));
394 else if (ix
< 0x4001d555) /* 6.6666259765625 */
395 { /* 6.666 > |x| >= 1/.35 ~ 2.857143 */
396 R
= rb
[0] + s
* (rb
[1] + s
* (rb
[2] + s
* (rb
[3] + s
* (rb
[4] +
397 s
* (rb
[5] + s
* (rb
[6] + s
* rb
[7]))))));
398 S
= sb
[0] + s
* (sb
[1] + s
* (sb
[2] + s
* (sb
[3] + s
* (sb
[4] +
399 s
* (sb
[5] + s
* (sb
[6] + s
))))));
404 return two
- tiny
; /* x < -6.666 */
406 R
= rc
[0] + s
* (rc
[1] + s
* (rc
[2] + s
* (rc
[3] +
407 s
* (rc
[4] + s
* rc
[5]))));
408 S
= sc
[0] + s
* (sc
[1] + s
* (sc
[2] + s
* (sc
[3] +
412 GET_LDOUBLE_WORDS (hx
, i0
, i1
, z
);
415 SET_LDOUBLE_WORDS (z
, hx
, i0
, i1
);
416 r
= expl (-z
* z
- 0.5625) *
417 expl ((z
- x
) * (z
+ x
) + R
/ S
);
418 if ((se
& 0x8000) == 0)
425 if ((se
& 0x8000) == 0)