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/* Long double expansions are
13 Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
14 and are incorporated herein by permission of the author. The author
15 reserves the right to distribute this material elsewhere under different
16 copying permissions. These modifications are distributed here under
17 the following terms:
18
19 This library is free software; you can redistribute it and/or
20 modify it under the terms of the GNU Lesser General Public
21 License as published by the Free Software Foundation; either
22 version 2.1 of the License, or (at your option) any later version.
23
24 This library is distributed in the hope that it will be useful,
25 but WITHOUT ANY WARRANTY; without even the implied warranty of
26 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
27 Lesser General Public License for more details.
28
29 You should have received a copy of the GNU Lesser General Public
30 License along with this library; if not, see
31 <http://www.gnu.org/licenses/>. */
32
33/* double erf(double x)
34 * double erfc(double x)
35 * x
36 * 2 |\
37 * erf(x) = --------- | exp(-t*t)dt
38 * sqrt(pi) \|
39 * 0
40 *
41 * erfc(x) = 1-erf(x)
42 * Note that
43 * erf(-x) = -erf(x)
44 * erfc(-x) = 2 - erfc(x)
45 *
46 * Method:
47 * 1. For |x| in [0, 0.84375]
48 * erf(x) = x + x*R(x^2)
49 * erfc(x) = 1 - erf(x) if x in [-.84375,0.25]
50 * = 0.5 + ((0.5-x)-x*R) if x in [0.25,0.84375]
51 * Remark. The formula is derived by noting
52 * erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....)
53 * and that
54 * 2/sqrt(pi) = 1.128379167095512573896158903121545171688
55 * is close to one. The interval is chosen because the fix
56 * point of erf(x) is near 0.6174 (i.e., erf(x)=x when x is
57 * near 0.6174), and by some experiment, 0.84375 is chosen to
58 * guarantee the error is less than one ulp for erf.
59 *
60 * 2. For |x| in [0.84375,1.25], let s = |x| - 1, and
61 * c = 0.84506291151 rounded to single (24 bits)
62 * erf(x) = sign(x) * (c + P1(s)/Q1(s))
63 * erfc(x) = (1-c) - P1(s)/Q1(s) if x > 0
64 * 1+(c+P1(s)/Q1(s)) if x < 0
65 * Remark: here we use the taylor series expansion at x=1.
66 * erf(1+s) = erf(1) + s*Poly(s)
67 * = 0.845.. + P1(s)/Q1(s)
68 * Note that |P1/Q1|< 0.078 for x in [0.84375,1.25]
69 *
70 * 3. For x in [1.25,1/0.35(~2.857143)],
71 * erfc(x) = (1/x)*exp(-x*x-0.5625+R1(z)/S1(z))
72 * z=1/x^2
73 * erf(x) = 1 - erfc(x)
74 *
75 * 4. For x in [1/0.35,107]
76 * erfc(x) = (1/x)*exp(-x*x-0.5625+R2/S2) if x > 0
77 * = 2.0 - (1/x)*exp(-x*x-0.5625+R2(z)/S2(z))
78 * if -6.666<x<0
79 * = 2.0 - tiny (if x <= -6.666)
80 * z=1/x^2
81 * erf(x) = sign(x)*(1.0 - erfc(x)) if x < 6.666, else
82 * erf(x) = sign(x)*(1.0 - tiny)
83 * Note1:
84 * To compute exp(-x*x-0.5625+R/S), let s be a single
85 * precision number and s := x; then
86 * -x*x = -s*s + (s-x)*(s+x)
87 * exp(-x*x-0.5626+R/S) =
88 * exp(-s*s-0.5625)*exp((s-x)*(s+x)+R/S);
89 * Note2:
90 * Here 4 and 5 make use of the asymptotic series
91 * exp(-x*x)
92 * erfc(x) ~ ---------- * ( 1 + Poly(1/x^2) )
93 * x*sqrt(pi)
94 *
95 * 5. For inf > x >= 107
96 * erf(x) = sign(x) *(1 - tiny) (raise inexact)
97 * erfc(x) = tiny*tiny (raise underflow) if x > 0
98 * = 2 - tiny if x<0
99 *
100 * 7. Special case:
101 * erf(0) = 0, erf(inf) = 1, erf(-inf) = -1,
102 * erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2,
103 * erfc/erf(NaN) is NaN
104 */
105
106
107#include <errno.h>
108#include <float.h>
109#include <math.h>
110#include <math_private.h>
111#include <libm-alias-ldouble.h>
112
113static const long double
114tiny = 1e-4931L,
115 half = 0.5L,
116 one = 1.0L,
117 two = 2.0L,
118 /* c = (float)0.84506291151 */
119 erx = 0.845062911510467529296875L,
120/*
121 * Coefficients for approximation to erf on [0,0.84375]
122 */
123 /* 2/sqrt(pi) - 1 */
124 efx = 1.2837916709551257389615890312154517168810E-1L,
125
126 pp[6] = {
127 1.122751350964552113068262337278335028553E6L,
128 -2.808533301997696164408397079650699163276E6L,
129 -3.314325479115357458197119660818768924100E5L,
130 -6.848684465326256109712135497895525446398E4L,
131 -2.657817695110739185591505062971929859314E3L,
132 -1.655310302737837556654146291646499062882E2L,
133 },
134
135 qq[6] = {
136 8.745588372054466262548908189000448124232E6L,
137 3.746038264792471129367533128637019611485E6L,
138 7.066358783162407559861156173539693900031E5L,
139 7.448928604824620999413120955705448117056E4L,
140 4.511583986730994111992253980546131408924E3L,
141 1.368902937933296323345610240009071254014E2L,
142 /* 1.000000000000000000000000000000000000000E0 */
143 },
144
145/*
146 * Coefficients for approximation to erf in [0.84375,1.25]
147 */
148/* erf(x+1) = 0.845062911510467529296875 + pa(x)/qa(x)
149 -0.15625 <= x <= +.25
150 Peak relative error 8.5e-22 */
151
152 pa[8] = {
153 -1.076952146179812072156734957705102256059E0L,
154 1.884814957770385593365179835059971587220E2L,
155 -5.339153975012804282890066622962070115606E1L,
156 4.435910679869176625928504532109635632618E1L,
157 1.683219516032328828278557309642929135179E1L,
158 -2.360236618396952560064259585299045804293E0L,
159 1.852230047861891953244413872297940938041E0L,
160 9.394994446747752308256773044667843200719E-2L,
161 },
162
163 qa[7] = {
164 4.559263722294508998149925774781887811255E2L,
165 3.289248982200800575749795055149780689738E2L,
166 2.846070965875643009598627918383314457912E2L,
167 1.398715859064535039433275722017479994465E2L,
168 6.060190733759793706299079050985358190726E1L,
169 2.078695677795422351040502569964299664233E1L,
170 4.641271134150895940966798357442234498546E0L,
171 /* 1.000000000000000000000000000000000000000E0 */
172 },
173
174/*
175 * Coefficients for approximation to erfc in [1.25,1/0.35]
176 */
177/* erfc(1/x) = x exp (-1/x^2 - 0.5625 + ra(x^2)/sa(x^2))
178 1/2.85711669921875 < 1/x < 1/1.25
179 Peak relative error 3.1e-21 */
180
181 ra[] = {
182 1.363566591833846324191000679620738857234E-1L,
183 1.018203167219873573808450274314658434507E1L,
184 1.862359362334248675526472871224778045594E2L,
185 1.411622588180721285284945138667933330348E3L,
186 5.088538459741511988784440103218342840478E3L,
187 8.928251553922176506858267311750789273656E3L,
188 7.264436000148052545243018622742770549982E3L,
189 2.387492459664548651671894725748959751119E3L,
190 2.220916652813908085449221282808458466556E2L,
191 },
192
193 sa[] = {
194 -1.382234625202480685182526402169222331847E1L,
195 -3.315638835627950255832519203687435946482E2L,
196 -2.949124863912936259747237164260785326692E3L,
197 -1.246622099070875940506391433635999693661E4L,
198 -2.673079795851665428695842853070996219632E4L,
199 -2.880269786660559337358397106518918220991E4L,
200 -1.450600228493968044773354186390390823713E4L,
201 -2.874539731125893533960680525192064277816E3L,
202 -1.402241261419067750237395034116942296027E2L,
203 /* 1.000000000000000000000000000000000000000E0 */
204 },
205/*
206 * Coefficients for approximation to erfc in [1/.35,107]
207 */
208/* erfc(1/x) = x exp (-1/x^2 - 0.5625 + rb(x^2)/sb(x^2))
209 1/6.6666259765625 < 1/x < 1/2.85711669921875
210 Peak relative error 4.2e-22 */
211 rb[] = {
212 -4.869587348270494309550558460786501252369E-5L,
213 -4.030199390527997378549161722412466959403E-3L,
214 -9.434425866377037610206443566288917589122E-2L,
215 -9.319032754357658601200655161585539404155E-1L,
216 -4.273788174307459947350256581445442062291E0L,
217 -8.842289940696150508373541814064198259278E0L,
218 -7.069215249419887403187988144752613025255E0L,
219 -1.401228723639514787920274427443330704764E0L,
220 },
221
222 sb[] = {
223 4.936254964107175160157544545879293019085E-3L,
224 1.583457624037795744377163924895349412015E-1L,
225 1.850647991850328356622940552450636420484E0L,
226 9.927611557279019463768050710008450625415E0L,
227 2.531667257649436709617165336779212114570E1L,
228 2.869752886406743386458304052862814690045E1L,
229 1.182059497870819562441683560749192539345E1L,
230 /* 1.000000000000000000000000000000000000000E0 */
231 },
232/* erfc(1/x) = x exp (-1/x^2 - 0.5625 + rc(x^2)/sc(x^2))
233 1/107 <= 1/x <= 1/6.6666259765625
234 Peak relative error 1.1e-21 */
235 rc[] = {
236 -8.299617545269701963973537248996670806850E-5L,
237 -6.243845685115818513578933902532056244108E-3L,
238 -1.141667210620380223113693474478394397230E-1L,
239 -7.521343797212024245375240432734425789409E-1L,
240 -1.765321928311155824664963633786967602934E0L,
241 -1.029403473103215800456761180695263439188E0L,
242 },
243
244 sc[] = {
245 8.413244363014929493035952542677768808601E-3L,
246 2.065114333816877479753334599639158060979E-1L,
247 1.639064941530797583766364412782135680148E0L,
248 4.936788463787115555582319302981666347450E0L,
249 5.005177727208955487404729933261347679090E0L,
250 /* 1.000000000000000000000000000000000000000E0 */
251 };
252
253long double
254__erfl (long double x)
255{
256 long double R, S, P, Q, s, y, z, r;
257 int32_t ix, i;
258 uint32_t se, i0, i1;
259
260 GET_LDOUBLE_WORDS (se, i0, i1, x);
261 ix = se & 0x7fff;
262
263 if (ix >= 0x7fff)
264 { /* erf(nan)=nan */
265 i = ((se & 0xffff) >> 15) << 1;
266 return (long double) (1 - i) + one / x; /* erf(+-inf)=+-1 */
267 }
268
269 ix = (ix << 16) | (i0 >> 16);
270 if (ix < 0x3ffed800) /* |x|<0.84375 */
271 {
272 if (ix < 0x3fde8000) /* |x|<2**-33 */
273 {
274 if (ix < 0x00080000)
275 {
276 /* Avoid spurious underflow. */
277 long double ret = 0.0625 * (16.0 * x + (16.0 * efx) * x);
278 math_check_force_underflow (ret);
279 return ret;
280 }
281 return x + efx * x;
282 }
283 z = x * x;
284 r = pp[0] + z * (pp[1]
285 + z * (pp[2] + z * (pp[3] + z * (pp[4] + z * pp[5]))));
286 s = qq[0] + z * (qq[1]
287 + z * (qq[2] + z * (qq[3] + z * (qq[4] + z * (qq[5] + z)))));
288 y = r / s;
289 return x + x * y;
290 }
291 if (ix < 0x3fffa000) /* 1.25 */
292 { /* 0.84375 <= |x| < 1.25 */
293 s = fabsl (x) - one;
294 P = pa[0] + s * (pa[1] + s * (pa[2]
295 + s * (pa[3] + s * (pa[4] + s * (pa[5] + s * (pa[6] + s * pa[7]))))));
296 Q = qa[0] + s * (qa[1] + s * (qa[2]
297 + s * (qa[3] + s * (qa[4] + s * (qa[5] + s * (qa[6] + s))))));
298 if ((se & 0x8000) == 0)
299 return erx + P / Q;
300 else
301 return -erx - P / Q;
302 }
303 if (ix >= 0x4001d555) /* 6.6666259765625 */
304 { /* inf>|x|>=6.666 */
305 if ((se & 0x8000) == 0)
306 return one - tiny;
307 else
308 return tiny - one;
309 }
310 x = fabsl (x);
311 s = one / (x * x);
312 if (ix < 0x4000b6db) /* 2.85711669921875 */
313 {
314 R = ra[0] + s * (ra[1] + s * (ra[2] + s * (ra[3] + s * (ra[4] +
315 s * (ra[5] + s * (ra[6] + s * (ra[7] + s * ra[8])))))));
316 S = sa[0] + s * (sa[1] + s * (sa[2] + s * (sa[3] + s * (sa[4] +
317 s * (sa[5] + s * (sa[6] + s * (sa[7] + s * (sa[8] + s))))))));
318 }
319 else
320 { /* |x| >= 1/0.35 */
321 R = rb[0] + s * (rb[1] + s * (rb[2] + s * (rb[3] + s * (rb[4] +
322 s * (rb[5] + s * (rb[6] + s * rb[7]))))));
323 S = sb[0] + s * (sb[1] + s * (sb[2] + s * (sb[3] + s * (sb[4] +
324 s * (sb[5] + s * (sb[6] + s))))));
325 }
326 z = x;
327 GET_LDOUBLE_WORDS (i, i0, i1, z);
328 i1 = 0;
329 SET_LDOUBLE_WORDS (z, i, i0, i1);
330 r =
331 __ieee754_expl (-z * z - 0.5625) * __ieee754_expl ((z - x) * (z + x) +
332 R / S);
333 if ((se & 0x8000) == 0)
334 return one - r / x;
335 else
336 return r / x - one;
337}
338
339libm_alias_ldouble (__erf, erf)
340long double
341__erfcl (long double x)
342{
343 int32_t hx, ix;
344 long double R, S, P, Q, s, y, z, r;
345 uint32_t se, i0, i1;
346
347 GET_LDOUBLE_WORDS (se, i0, i1, x);
348 ix = se & 0x7fff;
349 if (ix >= 0x7fff)
350 { /* erfc(nan)=nan */
351 /* erfc(+-inf)=0,2 */
352 return (long double) (((se & 0xffff) >> 15) << 1) + one / x;
353 }
354
355 ix = (ix << 16) | (i0 >> 16);
356 if (ix < 0x3ffed800) /* |x|<0.84375 */
357 {
358 if (ix < 0x3fbe0000) /* |x|<2**-65 */
359 return one - x;
360 z = x * x;
361 r = pp[0] + z * (pp[1]
362 + z * (pp[2] + z * (pp[3] + z * (pp[4] + z * pp[5]))));
363 s = qq[0] + z * (qq[1]
364 + z * (qq[2] + z * (qq[3] + z * (qq[4] + z * (qq[5] + z)))));
365 y = r / s;
366 if (ix < 0x3ffd8000) /* x<1/4 */
367 {
368 return one - (x + x * y);
369 }
370 else
371 {
372 r = x * y;
373 r += (x - half);
374 return half - r;
375 }
376 }
377 if (ix < 0x3fffa000) /* 1.25 */
378 { /* 0.84375 <= |x| < 1.25 */
379 s = fabsl (x) - one;
380 P = pa[0] + s * (pa[1] + s * (pa[2]
381 + s * (pa[3] + s * (pa[4] + s * (pa[5] + s * (pa[6] + s * pa[7]))))));
382 Q = qa[0] + s * (qa[1] + s * (qa[2]
383 + s * (qa[3] + s * (qa[4] + s * (qa[5] + s * (qa[6] + s))))));
384 if ((se & 0x8000) == 0)
385 {
386 z = one - erx;
387 return z - P / Q;
388 }
389 else
390 {
391 z = erx + P / Q;
392 return one + z;
393 }
394 }
395 if (ix < 0x4005d600) /* 107 */
396 { /* |x|<107 */
397 x = fabsl (x);
398 s = one / (x * x);
399 if (ix < 0x4000b6db) /* 2.85711669921875 */
400 { /* |x| < 1/.35 ~ 2.857143 */
401 R = ra[0] + s * (ra[1] + s * (ra[2] + s * (ra[3] + s * (ra[4] +
402 s * (ra[5] + s * (ra[6] + s * (ra[7] + s * ra[8])))))));
403 S = sa[0] + s * (sa[1] + s * (sa[2] + s * (sa[3] + s * (sa[4] +
404 s * (sa[5] + s * (sa[6] + s * (sa[7] + s * (sa[8] + s))))))));
405 }
406 else if (ix < 0x4001d555) /* 6.6666259765625 */
407 { /* 6.666 > |x| >= 1/.35 ~ 2.857143 */
408 R = rb[0] + s * (rb[1] + s * (rb[2] + s * (rb[3] + s * (rb[4] +
409 s * (rb[5] + s * (rb[6] + s * rb[7]))))));
410 S = sb[0] + s * (sb[1] + s * (sb[2] + s * (sb[3] + s * (sb[4] +
411 s * (sb[5] + s * (sb[6] + s))))));
412 }
413 else
414 { /* |x| >= 6.666 */
415 if (se & 0x8000)
416 return two - tiny; /* x < -6.666 */
417
418 R = rc[0] + s * (rc[1] + s * (rc[2] + s * (rc[3] +
419 s * (rc[4] + s * rc[5]))));
420 S = sc[0] + s * (sc[1] + s * (sc[2] + s * (sc[3] +
421 s * (sc[4] + s))));
422 }
423 z = x;
424 GET_LDOUBLE_WORDS (hx, i0, i1, z);
425 i1 = 0;
426 i0 &= 0xffffff00;
427 SET_LDOUBLE_WORDS (z, hx, i0, i1);
428 r = __ieee754_expl (-z * z - 0.5625) *
429 __ieee754_expl ((z - x) * (z + x) + R / S);
430 if ((se & 0x8000) == 0)
431 {
432 long double ret = r / x;
433 if (ret == 0)
434 __set_errno (ERANGE);
435 return ret;
436 }
437 else
438 return two - r / x;
439 }
440 else
441 {
442 if ((se & 0x8000) == 0)
443 {
444 __set_errno (ERANGE);
445 return tiny * tiny;
446 }
447 else
448 return two - tiny;
449 }
450}
451
452libm_alias_ldouble (__erfc, erfc)
453