1 /* @(#)s_erf.c 5.1 93/09/24 */
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
6 * Developed at SunPro, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
10 * ====================================================
13 //__FBSDID("$FreeBSD$");
16 * See s_erf.c for complete comments.
18 * Converted to long double by Steven G. Kargl.
24 #include "math_private.h"
26 #if defined(__x86_64__)
27 #include "x86_64/ieeefp.h"
28 #elif defined(__i386__)
29 #include "i386/ieeefp.h"
32 /* XXX Prevent compilers from erroneously constant folding: */
33 static const volatile long double tiny
__attribute__ ((__section__(".rodata,\"a\" " SECTIONCOMMENT
))) = 0x1p
-10000L;
40 * In the domain [0, 2**-34], only the first term in the power series
41 * expansion of erf(x) is used. The magnitude of the first neglected
42 * terms is less than 2**-102.
44 static const union IEEEl2bits
45 efxu
= LD80C(0x8375d410a6db446c, -3, 1.28379167095512573902e-1L),
46 efx8u
= LD80C(0x8375d410a6db446c, 0, 1.02703333676410059122e+0L),
48 * Domain [0, 0.84375], range ~[-1.423e-22, 1.423e-22]:
49 * |(erf(x) - x)/x - pp(x)/qq(x)| < 2**-72.573
51 pp0u
= LD80C(0x8375d410a6db446c, -3, 1.28379167095512573902e-1L),
52 pp1u
= LD80C(0xa46c7d09ec3d0cec, -2, -3.21140201054840180596e-1L),
53 pp2u
= LD80C(0x9b31e66325576f86, -5, -3.78893851760347812082e-2L),
54 pp3u
= LD80C(0x804ac72c9a0b97dd, -7, -7.83032847030604679616e-3L),
55 pp4u
= LD80C(0x9f42bcbc3d5a601d, -12, -3.03765663857082048459e-4L),
56 pp5u
= LD80C(0x9ec4ad6193470693, -16, -1.89266527398167917502e-5L),
57 qq1u
= LD80C(0xdb4b8eb713188d6b, -2, 4.28310832832310510579e-1L),
58 qq2u
= LD80C(0xa5750835b2459bd1, -4, 8.07896272074540216658e-2L),
59 qq3u
= LD80C(0x8b85d6bd6a90b51c, -7, 8.51579638189385354266e-3L),
60 qq4u
= LD80C(0x87332f82cff4ff96, -11, 5.15746855583604912827e-4L),
61 qq5u
= LD80C(0x83466cb6bf9dca00, -16, 1.56492109706256700009e-5L),
62 qq6u
= LD80C(0xf5bf98c2f996bf63, -24, 1.14435527803073879724e-7L);
64 #define efx8 (efx8u.e)
77 static const union IEEEl2bits
78 erxu
= LD80C(0xd7bb3d0000000000, -1, 8.42700779438018798828e-1L),
80 * Domain [0.84375, 1.25], range ~[-8.132e-22, 8.113e-22]:
81 * |(erf(x) - erx) - pa(x)/qa(x)| < 2**-71.762
83 pa0u
= LD80C(0xe8211158da02c692, -27, 1.35116960705131296711e-8L),
84 pa1u
= LD80C(0xd488f89f36988618, -2, 4.15107507167065612570e-1L),
85 pa2u
= LD80C(0xece74f8c63fa3942, -4, -1.15675565215949226989e-1L),
86 pa3u
= LD80C(0xc8d31e020727c006, -4, 9.80589241379624665791e-2L),
87 pa4u
= LD80C(0x985d5d5fafb0551f, -5, 3.71984145558422368847e-2L),
88 pa5u
= LD80C(0xa5b6c4854d2f5452, -8, -5.05718799340957673661e-3L),
89 pa6u
= LD80C(0x85c8d58fe3993a47, -8, 4.08277919612202243721e-3L),
90 pa7u
= LD80C(0xddbfbc23677b35cf, -13, 2.11476292145347530794e-4L),
91 qa1u
= LD80C(0xb8a977896f5eff3f, -1, 7.21335860303380361298e-1L),
92 qa2u
= LD80C(0x9fcd662c3d4eac86, -1, 6.24227891731886593333e-1L),
93 qa3u
= LD80C(0x9d0b618eac67ba07, -2, 3.06727455774491855801e-1L),
94 qa4u
= LD80C(0x881a4293f6d6c92d, -3, 1.32912674218195890535e-1L),
95 qa5u
= LD80C(0xbab144f07dea45bf, -5, 4.55792134233613027584e-2L),
96 qa6u
= LD80C(0xa6c34ba438bdc900, -7, 1.01783980070527682680e-2L),
97 qa7u
= LD80C(0x8fa866dc20717a91, -9, 2.19204436518951438183e-3L);
114 static const union IEEEl2bits
116 * Domain [1.25,2.85715], range ~[-2.334e-22,2.334e-22]:
117 * |log(x*erfc(x)) + x**2 + 0.5625 - ra(x)/sa(x)| < 2**-71.860
119 ra0u
= LD80C(0xa1a091e0fb4f335a, -7, -9.86494298915814308249e-3L),
120 ra1u
= LD80C(0xc2b0d045ae37df6b, -1, -7.60510460864878271275e-1L),
121 ra2u
= LD80C(0xf2cec3ee7da636c5, 3, -1.51754798236892278250e+1L),
122 ra3u
= LD80C(0x813cc205395adc7d, 7, -1.29237335516455333420e+2L),
123 ra4u
= LD80C(0x8737c8b7b4062c2f, 9, -5.40871625829510494776e+2L),
124 ra5u
= LD80C(0x8ffe5383c08d4943, 10, -1.15194769466026108551e+3L),
125 ra6u
= LD80C(0x983573e64d5015a9, 10, -1.21767039790249025544e+3L),
126 ra7u
= LD80C(0x92a794e763a6d4db, 9, -5.86618463370624636688e+2L),
127 ra8u
= LD80C(0xd5ad1fae77c3d9a3, 6, -1.06838132335777049840e+2L),
128 ra9u
= LD80C(0x934c1a247807bb9c, 2, -4.60303980944467334806e+0L),
129 sa1u
= LD80C(0xd342f90012bb1189, 4, 2.64077014928547064865e+1L),
130 sa2u
= LD80C(0x839be13d9d5da883, 8, 2.63217811300123973067e+2L),
131 sa3u
= LD80C(0x9f8cba6d1ae1b24b, 10, 1.27639775710344617587e+3L),
132 sa4u
= LD80C(0xcaa83f403713e33e, 11, 3.24251544209971162003e+3L),
133 sa5u
= LD80C(0x8796aff2f3c47968, 12, 4.33883591261332837874e+3L),
134 sa6u
= LD80C(0xb6ef97f9c753157b, 11, 2.92697460344182158454e+3L),
135 sa7u
= LD80C(0xe02aee5f83773d1c, 9, 8.96670799139389559818e+2L),
136 sa8u
= LD80C(0xc82b83855b88e07e, 6, 1.00084987800048510018e+2L),
137 sa9u
= LD80C(0x92f030aefadf28ad, 1, 2.29591004455459083843e+0L);
158 * Domain [2.85715,7], range ~[-8.323e-22,8.390e-22]:
159 * |log(x*erfc(x)) + x**2 + 0.5625 - rb(x)/sb(x)| < 2**-70.326
161 static const union IEEEl2bits
162 rb0u
= LD80C(0xa1a091cf43abcd26, -7, -9.86494292470284646962e-3L),
163 rb1u
= LD80C(0xd19d2df1cbb8da0a, -1, -8.18804618389296662837e-1L),
164 rb2u
= LD80C(0x9a4dd1383e5daf5b, 4, -1.92879967111618594779e+1L),
165 rb3u
= LD80C(0xbff0ae9fc0751de6, 7, -1.91940164551245394969e+2L),
166 rb4u
= LD80C(0xdde08465310b472b, 9, -8.87508080766577324539e+2L),
167 rb5u
= LD80C(0xe796e1d38c8c70a9, 10, -1.85271506669474503781e+3L),
168 rb6u
= LD80C(0xbaf655a76e0ab3b5, 10, -1.49569795581333675349e+3L),
169 rb7u
= LD80C(0x95d21e3e75503c21, 8, -2.99641547972948019157e+2L),
170 sb1u
= LD80C(0x814487ed823c8cbd, 5, 3.23169247732868256569e+1L),
171 sb2u
= LD80C(0xbe4bfbb1301304be, 8, 3.80593618534539961773e+2L),
172 sb3u
= LD80C(0x809c4ade46b927c7, 11, 2.05776827838541292848e+3L),
173 sb4u
= LD80C(0xa55284359f3395a8, 12, 5.29031455540062116327e+3L),
174 sb5u
= LD80C(0xbcfa72da9b820874, 12, 6.04730608102312640462e+3L),
175 sb6u
= LD80C(0x9d09a35988934631, 11, 2.51260238030767176221e+3L),
176 sb7u
= LD80C(0xd675bbe542c159fa, 7, 2.14459898308561015684e+2L);
193 * Domain [7,108], range ~[-4.422e-22,4.422e-22]:
194 * |log(x*erfc(x)) + x**2 + 0.5625 - rc(x)/sc(x)| < 2**-70.938
196 static const union IEEEl2bits
197 /* err = -4.422092275318925082e-22 -70.937689 */
198 rc0u
= LD80C(0xa1a091cf437a17ad, -7, -9.86494292470008707260e-3L),
199 rc1u
= LD80C(0xbe79c5a978122b00, -1, -7.44045595049165939261e-1L),
200 rc2u
= LD80C(0xdb26f9bbe31a2794, 3, -1.36970155085888424425e+1L),
201 rc3u
= LD80C(0xb5f69a38f5747ac8, 6, -9.09816453742625888546e+1L),
202 rc4u
= LD80C(0xd79676d970d0a21a, 7, -2.15587750997584074147e+2L),
203 rc5u
= LD80C(0xfe528153c45ec97c, 6, -1.27161142938347796666e+2L),
204 sc1u
= LD80C(0xc5e8cd46d5604a96, 4, 2.47386727842204312937e+1L),
205 sc2u
= LD80C(0xc5f0f5a5484520eb, 7, 1.97941248254913378865e+2L),
206 sc3u
= LD80C(0x964e3c7b34db9170, 9, 6.01222441484087787522e+2L),
207 sc4u
= LD80C(0x99be1b89faa0596a, 9, 6.14970430845978077827e+2L),
208 sc5u
= LD80C(0xf80dfcbf37ffc5ea, 6, 1.24027318931184605891e+2L);
224 long double ax
,R
,S
,P
,Q
,s
,y
,z
,r
;
225 __unused
uint64_t lx
;
229 EXTRACT_LDBL80_WORDS(hx
, lx
, x
);
231 if((hx
& 0x7fff) == 0x7fff) { /* erfl(nan)=nan */
233 return (1-i
)+one
/x
; /* erfl(+-inf)=+-1 */
242 RETURNI((8*x
+efx8
*x
)/8); /* avoid spurious underflow */
246 r
= pp0
+z
*(pp1
+z
*(pp2
+z
*(pp3
+z
*(pp4
+z
*pp5
))));
247 s
= one
+z
*(qq1
+z
*(qq2
+z
*(qq3
+z
*(qq4
+z
*(qq5
+z
*qq6
)))));
253 P
= pa0
+s
*(pa1
+s
*(pa2
+s
*(pa3
+s
*(pa4
+s
*(pa5
+s
*(pa6
+s
*pa7
))))));
254 Q
= one
+s
*(qa1
+s
*(qa2
+s
*(qa3
+s
*(qa4
+s
*(qa5
+s
*(qa6
+s
*qa7
))))));
255 if(x
>=0) RETURNI(erx
+ P
/Q
); else RETURNI(-erx
- P
/Q
);
257 if(ax
>= 7) { /* inf>|x|>= 7 */
258 if(x
>=0) RETURNI(one
-tiny
); else RETURNI(tiny
-one
);
261 if(ax
< 2.85715) { /* |x| < 2.85715 */
262 R
=ra0
+s
*(ra1
+s
*(ra2
+s
*(ra3
+s
*(ra4
+s
*(ra5
+s
*(ra6
+s
*(ra7
+
263 s
*(ra8
+s
*ra9
))))))));
264 S
=one
+s
*(sa1
+s
*(sa2
+s
*(sa3
+s
*(sa4
+s
*(sa5
+s
*(sa6
+s
*(sa7
+
265 s
*(sa8
+s
*sa9
))))))));
266 } else { /* |x| >= 2.85715 */
267 R
=rb0
+s
*(rb1
+s
*(rb2
+s
*(rb3
+s
*(rb4
+s
*(rb5
+s
*(rb6
+s
*rb7
))))));
268 S
=one
+s
*(sb1
+s
*(sb2
+s
*(sb3
+s
*(sb4
+s
*(sb5
+s
*(sb6
+s
*sb7
))))));
271 r
=expl(-z
*z
-0.5625)*expl((z
-ax
)*(z
+ax
)+R
/S
);
272 if(x
>=0) RETURNI(one
-r
/ax
); else RETURNI(r
/ax
-one
);
278 long double ax
,R
,S
,P
,Q
,s
,y
,z
,r
;
279 __unused
uint64_t lx
;
282 EXTRACT_LDBL80_WORDS(hx
, lx
, x
);
284 if((hx
& 0x7fff) == 0x7fff) { /* erfcl(nan)=nan */
285 /* erfcl(+-inf)=0,2 */
286 return ((hx
>>15)<<1)+one
/x
;
296 r
= pp0
+z
*(pp1
+z
*(pp2
+z
*(pp3
+z
*(pp4
+z
*pp5
))));
297 s
= one
+z
*(qq1
+z
*(qq2
+z
*(qq3
+z
*(qq4
+z
*(qq5
+z
*qq6
)))));
299 if(ax
< 0.25L) { /* x<1/4 */
300 RETURNI(one
-(x
+x
*y
));
309 P
= pa0
+s
*(pa1
+s
*(pa2
+s
*(pa3
+s
*(pa4
+s
*(pa5
+s
*(pa6
+s
*pa7
))))));
310 Q
= one
+s
*(qa1
+s
*(qa2
+s
*(qa3
+s
*(qa4
+s
*(qa5
+s
*(qa6
+s
*qa7
))))));
312 z
= one
-erx
; RETURNI(z
- P
/Q
);
314 z
= (erx
+P
/Q
); RETURNI(one
+z
);
318 if(ax
< 108) { /* |x| < 108 */
320 if(ax
< 2.85715) { /* |x| < 2.85715 */
321 R
=ra0
+s
*(ra1
+s
*(ra2
+s
*(ra3
+s
*(ra4
+s
*(ra5
+s
*(ra6
+s
*(ra7
+
322 s
*(ra8
+s
*ra9
))))))));
323 S
=one
+s
*(sa1
+s
*(sa2
+s
*(sa3
+s
*(sa4
+s
*(sa5
+s
*(sa6
+s
*(sa7
+
324 s
*(sa8
+s
*sa9
))))))));
325 } else if(ax
< 7) { /* | |x| < 7 */
326 R
=rb0
+s
*(rb1
+s
*(rb2
+s
*(rb3
+s
*(rb4
+s
*(rb5
+s
*(rb6
+s
*rb7
))))));
327 S
=one
+s
*(sb1
+s
*(sb2
+s
*(sb3
+s
*(sb4
+s
*(sb5
+s
*(sb6
+s
*sb7
))))));
329 if(x
< -7) RETURNI(two
-tiny
);/* x < -7 */
330 R
=rc0
+s
*(rc1
+s
*(rc2
+s
*(rc3
+s
*(rc4
+s
*rc5
))));
331 S
=one
+s
*(sc1
+s
*(sc2
+s
*(sc3
+s
*(sc4
+s
*sc5
))));
334 r
= expl(-z
*z
-0.5625)*expl((z
-ax
)*(z
+ax
)+R
/S
);
335 if(x
>0) RETURNI(r
/ax
); else RETURNI(two
-r
/ax
);
337 if(x
>0) RETURNI(tiny
*tiny
); else RETURNI(two
-tiny
);