Branch data Line data Source code
1 : : /*
2 : : * ====================================================
3 : : * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
4 : : *
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
8 : : * is preserved.
9 : : * ====================================================
10 : : */
11 : :
12 : : /*
13 : : * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
14 : : *
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.
18 : : *
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.
26 : : */
27 : :
28 : : /* double erf(double x)
29 : : * double erfc(double x)
30 : : * x
31 : : * 2 |\
32 : : * erf(x) = --------- | exp(-t*t)dt
33 : : * sqrt(pi) \|
34 : : * 0
35 : : *
36 : : * erfc(x) = 1-erf(x)
37 : : * Note that
38 : : * erf(-x) = -erf(x)
39 : : * erfc(-x) = 2 - erfc(x)
40 : : *
41 : : * Method:
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 + ....)
48 : : * and that
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.
54 : : *
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]
64 : : *
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))
67 : : * z=1/x^2
68 : : * erf(x) = 1 - erfc(x)
69 : : *
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))
73 : : * if -6.666<x<0
74 : : * = 2.0 - tiny (if x <= -6.666)
75 : : * z=1/x^2
76 : : * erf(x) = sign(x)*(1.0 - erfc(x)) if x < 6.666, else
77 : : * erf(x) = sign(x)*(1.0 - tiny)
78 : : * Note1:
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);
84 : : * Note2:
85 : : * Here 4 and 5 make use of the asymptotic series
86 : : * exp(-x*x)
87 : : * erfc(x) ~ ---------- * ( 1 + Poly(1/x^2) )
88 : : * x*sqrt(pi)
89 : : *
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
93 : : * = 2 - tiny if x<0
94 : : *
95 : : * 7. Special case:
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
99 : : */
100 : :
101 : :
102 : : #include <openlibm_math.h>
103 : :
104 : : #include "math_private.h"
105 : :
106 : : static const long double
107 : : tiny = 1e-4931L,
108 : : half = 0.5L,
109 : : one = 1.0L,
110 : : two = 2.0L,
111 : : /* c = (float)0.84506291151 */
112 : : erx = 0.845062911510467529296875L,
113 : : /*
114 : : * Coefficients for approximation to erf on [0,0.84375]
115 : : */
116 : : /* 2/sqrt(pi) - 1 */
117 : : efx = 1.2837916709551257389615890312154517168810E-1L,
118 : : /* 8 * (2/sqrt(pi) - 1) */
119 : : efx8 = 1.0270333367641005911692712249723613735048E0L,
120 : :
121 : : pp[6] = {
122 : : 1.122751350964552113068262337278335028553E6L,
123 : : -2.808533301997696164408397079650699163276E6L,
124 : : -3.314325479115357458197119660818768924100E5L,
125 : : -6.848684465326256109712135497895525446398E4L,
126 : : -2.657817695110739185591505062971929859314E3L,
127 : : -1.655310302737837556654146291646499062882E2L,
128 : : },
129 : :
130 : : qq[6] = {
131 : : 8.745588372054466262548908189000448124232E6L,
132 : : 3.746038264792471129367533128637019611485E6L,
133 : : 7.066358783162407559861156173539693900031E5L,
134 : : 7.448928604824620999413120955705448117056E4L,
135 : : 4.511583986730994111992253980546131408924E3L,
136 : : 1.368902937933296323345610240009071254014E2L,
137 : : /* 1.000000000000000000000000000000000000000E0 */
138 : : },
139 : :
140 : : /*
141 : : * Coefficients for approximation to erf in [0.84375,1.25]
142 : : */
143 : : /* erf(x+1) = 0.845062911510467529296875 + pa(x)/qa(x)
144 : : -0.15625 <= x <= +.25
145 : : Peak relative error 8.5e-22 */
146 : :
147 : : pa[8] = {
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,
156 : : },
157 : :
158 : : qa[7] = {
159 : : 4.559263722294508998149925774781887811255E2L,
160 : : 3.289248982200800575749795055149780689738E2L,
161 : : 2.846070965875643009598627918383314457912E2L,
162 : : 1.398715859064535039433275722017479994465E2L,
163 : : 6.060190733759793706299079050985358190726E1L,
164 : : 2.078695677795422351040502569964299664233E1L,
165 : : 4.641271134150895940966798357442234498546E0L,
166 : : /* 1.000000000000000000000000000000000000000E0 */
167 : : },
168 : :
169 : : /*
170 : : * Coefficients for approximation to erfc in [1.25,1/0.35]
171 : : */
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 */
175 : :
176 : : ra[] = {
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,
186 : : },
187 : :
188 : : sa[] = {
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 */
199 : : },
200 : : /*
201 : : * Coefficients for approximation to erfc in [1/.35,107]
202 : : */
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 */
206 : : rb[] = {
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,
215 : : },
216 : :
217 : : sb[] = {
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 */
226 : : },
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 */
230 : : rc[] = {
231 : : -8.299617545269701963973537248996670806850E-5L,
232 : : -6.243845685115818513578933902532056244108E-3L,
233 : : -1.141667210620380223113693474478394397230E-1L,
234 : : -7.521343797212024245375240432734425789409E-1L,
235 : : -1.765321928311155824664963633786967602934E0L,
236 : : -1.029403473103215800456761180695263439188E0L,
237 : : },
238 : :
239 : : sc[] = {
240 : : 8.413244363014929493035952542677768808601E-3L,
241 : : 2.065114333816877479753334599639158060979E-1L,
242 : : 1.639064941530797583766364412782135680148E0L,
243 : : 4.936788463787115555582319302981666347450E0L,
244 : : 5.005177727208955487404729933261347679090E0L,
245 : : /* 1.000000000000000000000000000000000000000E0 */
246 : : };
247 : :
248 : : long double
249 : 0 : erfl(long double x)
250 : : {
251 : : long double R, S, P, Q, s, y, z, r;
252 : : int32_t ix, i;
253 : : u_int32_t se, i0, i1;
254 : :
255 : 0 : GET_LDOUBLE_WORDS (se, i0, i1, x);
256 : 0 : ix = se & 0x7fff;
257 : :
258 [ # # ]: 0 : if (ix >= 0x7fff)
259 : : { /* erf(nan)=nan */
260 : 0 : i = ((se & 0xffff) >> 15) << 1;
261 : 0 : return (long double) (1 - i) + one / x; /* erf(+-inf)=+-1 */
262 : : }
263 : :
264 : 0 : ix = (ix << 16) | (i0 >> 16);
265 [ # # ]: 0 : if (ix < 0x3ffed800) /* |x|<0.84375 */
266 : : {
267 [ # # ]: 0 : if (ix < 0x3fde8000) /* |x|<2**-33 */
268 : : {
269 [ # # ]: 0 : if (ix < 0x00080000)
270 : 0 : return 0.125 * (8.0 * x + efx8 * x); /*avoid underflow */
271 : 0 : return x + efx * x;
272 : : }
273 : 0 : z = x * x;
274 : 0 : r = pp[0] + z * (pp[1]
275 : 0 : + z * (pp[2] + z * (pp[3] + z * (pp[4] + z * pp[5]))));
276 : 0 : s = qq[0] + z * (qq[1]
277 : 0 : + z * (qq[2] + z * (qq[3] + z * (qq[4] + z * (qq[5] + z)))));
278 : 0 : y = r / s;
279 : 0 : return x + x * y;
280 : : }
281 [ # # ]: 0 : if (ix < 0x3fffa000) /* 1.25 */
282 : : { /* 0.84375 <= |x| < 1.25 */
283 : 0 : s = fabsl (x) - one;
284 : 0 : P = pa[0] + s * (pa[1] + s * (pa[2]
285 : 0 : + s * (pa[3] + s * (pa[4] + s * (pa[5] + s * (pa[6] + s * pa[7]))))));
286 : 0 : Q = qa[0] + s * (qa[1] + s * (qa[2]
287 : 0 : + s * (qa[3] + s * (qa[4] + s * (qa[5] + s * (qa[6] + s))))));
288 [ # # ]: 0 : if ((se & 0x8000) == 0)
289 : 0 : return erx + P / Q;
290 : : else
291 : 0 : return -erx - P / Q;
292 : : }
293 [ # # ]: 0 : if (ix >= 0x4001d555) /* 6.6666259765625 */
294 : : { /* inf>|x|>=6.666 */
295 [ # # ]: 0 : if ((se & 0x8000) == 0)
296 : 0 : return one - tiny;
297 : : else
298 : 0 : return tiny - one;
299 : : }
300 : 0 : x = fabsl (x);
301 : 0 : s = one / (x * x);
302 [ # # ]: 0 : if (ix < 0x4000b6db) /* 2.85711669921875 */
303 : : {
304 : 0 : R = ra[0] + s * (ra[1] + s * (ra[2] + s * (ra[3] + s * (ra[4] +
305 : 0 : s * (ra[5] + s * (ra[6] + s * (ra[7] + s * ra[8])))))));
306 : 0 : S = sa[0] + s * (sa[1] + s * (sa[2] + s * (sa[3] + s * (sa[4] +
307 : 0 : s * (sa[5] + s * (sa[6] + s * (sa[7] + s * (sa[8] + s))))))));
308 : : }
309 : : else
310 : : { /* |x| >= 1/0.35 */
311 : 0 : R = rb[0] + s * (rb[1] + s * (rb[2] + s * (rb[3] + s * (rb[4] +
312 : 0 : s * (rb[5] + s * (rb[6] + s * rb[7]))))));
313 : 0 : S = sb[0] + s * (sb[1] + s * (sb[2] + s * (sb[3] + s * (sb[4] +
314 : 0 : s * (sb[5] + s * (sb[6] + s))))));
315 : : }
316 : 0 : z = x;
317 : 0 : GET_LDOUBLE_WORDS (i, i0, i1, z);
318 : 0 : i1 = 0;
319 : 0 : SET_LDOUBLE_WORDS (z, i, i0, i1);
320 : 0 : r =
321 : 0 : expl (-z * z - 0.5625) * expl ((z - x) * (z + x) + R / S);
322 [ # # ]: 0 : if ((se & 0x8000) == 0)
323 : 0 : return one - r / x;
324 : : else
325 : 0 : return r / x - one;
326 : : }
327 : :
328 : : long double
329 : 0 : erfcl(long double x)
330 : : {
331 : : int32_t hx, ix;
332 : : long double R, S, P, Q, s, y, z, r;
333 : : u_int32_t se, i0, i1;
334 : :
335 : 0 : GET_LDOUBLE_WORDS (se, i0, i1, x);
336 : 0 : ix = se & 0x7fff;
337 [ # # ]: 0 : if (ix >= 0x7fff)
338 : : { /* erfc(nan)=nan */
339 : : /* erfc(+-inf)=0,2 */
340 : 0 : return (long double) (((se & 0xffff) >> 15) << 1) + one / x;
341 : : }
342 : :
343 : 0 : ix = (ix << 16) | (i0 >> 16);
344 [ # # ]: 0 : if (ix < 0x3ffed800) /* |x|<0.84375 */
345 : : {
346 [ # # ]: 0 : if (ix < 0x3fbe0000) /* |x|<2**-65 */
347 : 0 : return one - x;
348 : 0 : z = x * x;
349 : 0 : r = pp[0] + z * (pp[1]
350 : 0 : + z * (pp[2] + z * (pp[3] + z * (pp[4] + z * pp[5]))));
351 : 0 : s = qq[0] + z * (qq[1]
352 : 0 : + z * (qq[2] + z * (qq[3] + z * (qq[4] + z * (qq[5] + z)))));
353 : 0 : y = r / s;
354 [ # # ]: 0 : if (ix < 0x3ffd8000) /* x<1/4 */
355 : : {
356 : 0 : return one - (x + x * y);
357 : : }
358 : : else
359 : : {
360 : 0 : r = x * y;
361 : 0 : r += (x - half);
362 : 0 : return half - r;
363 : : }
364 : : }
365 [ # # ]: 0 : if (ix < 0x3fffa000) /* 1.25 */
366 : : { /* 0.84375 <= |x| < 1.25 */
367 : 0 : s = fabsl (x) - one;
368 : 0 : P = pa[0] + s * (pa[1] + s * (pa[2]
369 : 0 : + s * (pa[3] + s * (pa[4] + s * (pa[5] + s * (pa[6] + s * pa[7]))))));
370 : 0 : Q = qa[0] + s * (qa[1] + s * (qa[2]
371 : 0 : + s * (qa[3] + s * (qa[4] + s * (qa[5] + s * (qa[6] + s))))));
372 [ # # ]: 0 : if ((se & 0x8000) == 0)
373 : : {
374 : 0 : z = one - erx;
375 : 0 : return z - P / Q;
376 : : }
377 : : else
378 : : {
379 : 0 : z = erx + P / Q;
380 : 0 : return one + z;
381 : : }
382 : : }
383 [ # # ]: 0 : if (ix < 0x4005d600) /* 107 */
384 : : { /* |x|<107 */
385 : 0 : x = fabsl (x);
386 : 0 : s = one / (x * x);
387 [ # # ]: 0 : if (ix < 0x4000b6db) /* 2.85711669921875 */
388 : : { /* |x| < 1/.35 ~ 2.857143 */
389 : 0 : R = ra[0] + s * (ra[1] + s * (ra[2] + s * (ra[3] + s * (ra[4] +
390 : 0 : s * (ra[5] + s * (ra[6] + s * (ra[7] + s * ra[8])))))));
391 : 0 : S = sa[0] + s * (sa[1] + s * (sa[2] + s * (sa[3] + s * (sa[4] +
392 : 0 : s * (sa[5] + s * (sa[6] + s * (sa[7] + s * (sa[8] + s))))))));
393 : : }
394 [ # # ]: 0 : else if (ix < 0x4001d555) /* 6.6666259765625 */
395 : : { /* 6.666 > |x| >= 1/.35 ~ 2.857143 */
396 : 0 : R = rb[0] + s * (rb[1] + s * (rb[2] + s * (rb[3] + s * (rb[4] +
397 : 0 : s * (rb[5] + s * (rb[6] + s * rb[7]))))));
398 : 0 : S = sb[0] + s * (sb[1] + s * (sb[2] + s * (sb[3] + s * (sb[4] +
399 : 0 : s * (sb[5] + s * (sb[6] + s))))));
400 : : }
401 : : else
402 : : { /* |x| >= 6.666 */
403 [ # # ]: 0 : if (se & 0x8000)
404 : 0 : return two - tiny; /* x < -6.666 */
405 : :
406 : 0 : R = rc[0] + s * (rc[1] + s * (rc[2] + s * (rc[3] +
407 : 0 : s * (rc[4] + s * rc[5]))));
408 : 0 : S = sc[0] + s * (sc[1] + s * (sc[2] + s * (sc[3] +
409 : 0 : s * (sc[4] + s))));
410 : : }
411 : 0 : z = x;
412 : 0 : GET_LDOUBLE_WORDS (hx, i0, i1, z);
413 : 0 : i1 = 0;
414 : 0 : i0 &= 0xffffff00;
415 : 0 : SET_LDOUBLE_WORDS (z, hx, i0, i1);
416 : 0 : r = expl (-z * z - 0.5625) *
417 : 0 : expl ((z - x) * (z + x) + R / S);
418 [ # # ]: 0 : if ((se & 0x8000) == 0)
419 : 0 : return r / x;
420 : : else
421 : 0 : return two - r / x;
422 : : }
423 : : else
424 : : {
425 [ # # ]: 0 : if ((se & 0x8000) == 0)
426 : 0 : return tiny * tiny;
427 : : else
428 : 0 : return two - tiny;
429 : : }
430 : : }
|
