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 <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 |
|
|
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 |
|
|
GET_LDOUBLE_WORDS (se, i0, i1, x); |
256 |
|
|
ix = se & 0x7fff; |
257 |
|
|
|
258 |
|
|
if (ix >= 0x7fff) |
259 |
|
|
{ /* erf(nan)=nan */ |
260 |
|
|
i = ((se & 0xffff) >> 15) << 1; |
261 |
|
|
return (long double) (1 - i) + one / x; /* erf(+-inf)=+-1 */ |
262 |
|
|
} |
263 |
|
|
|
264 |
|
|
ix = (ix << 16) | (i0 >> 16); |
265 |
|
|
if (ix < 0x3ffed800) /* |x|<0.84375 */ |
266 |
|
|
{ |
267 |
|
|
if (ix < 0x3fde8000) /* |x|<2**-33 */ |
268 |
|
|
{ |
269 |
|
|
if (ix < 0x00080000) |
270 |
|
|
return 0.125 * (8.0 * x + efx8 * x); /*avoid underflow */ |
271 |
|
|
return x + efx * x; |
272 |
|
|
} |
273 |
|
|
z = x * x; |
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))))); |
278 |
|
|
y = r / s; |
279 |
|
|
return x + x * y; |
280 |
|
|
} |
281 |
|
|
if (ix < 0x3fffa000) /* 1.25 */ |
282 |
|
|
{ /* 0.84375 <= |x| < 1.25 */ |
283 |
|
|
s = fabsl (x) - one; |
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) |
289 |
|
|
return erx + P / Q; |
290 |
|
|
else |
291 |
|
|
return -erx - P / Q; |
292 |
|
|
} |
293 |
|
|
if (ix >= 0x4001d555) /* 6.6666259765625 */ |
294 |
|
|
{ /* inf>|x|>=6.666 */ |
295 |
|
|
if ((se & 0x8000) == 0) |
296 |
|
|
return one - tiny; |
297 |
|
|
else |
298 |
|
|
return tiny - one; |
299 |
|
|
} |
300 |
|
|
x = fabsl (x); |
301 |
|
|
s = one / (x * x); |
302 |
|
|
if (ix < 0x4000b6db) /* 2.85711669921875 */ |
303 |
|
|
{ |
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)))))))); |
308 |
|
|
} |
309 |
|
|
else |
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)))))); |
315 |
|
|
} |
316 |
|
|
z = x; |
317 |
|
|
GET_LDOUBLE_WORDS (i, i0, i1, z); |
318 |
|
|
i1 = 0; |
319 |
|
|
SET_LDOUBLE_WORDS (z, i, i0, i1); |
320 |
|
|
r = |
321 |
|
|
expl (-z * z - 0.5625) * expl ((z - x) * (z + x) + R / S); |
322 |
|
|
if ((se & 0x8000) == 0) |
323 |
|
|
return one - r / x; |
324 |
|
|
else |
325 |
|
|
return r / x - one; |
326 |
|
|
} |
327 |
|
|
DEF_STD(erfl); |
328 |
|
|
|
329 |
|
|
long double |
330 |
|
|
erfcl(long double x) |
331 |
|
|
{ |
332 |
|
|
int32_t hx, ix; |
333 |
|
|
long double R, S, P, Q, s, y, z, r; |
334 |
|
|
u_int32_t se, i0, i1; |
335 |
|
|
|
336 |
|
|
GET_LDOUBLE_WORDS (se, i0, i1, x); |
337 |
|
|
ix = se & 0x7fff; |
338 |
|
|
if (ix >= 0x7fff) |
339 |
|
|
{ /* erfc(nan)=nan */ |
340 |
|
|
/* erfc(+-inf)=0,2 */ |
341 |
|
|
return (long double) (((se & 0xffff) >> 15) << 1) + one / x; |
342 |
|
|
} |
343 |
|
|
|
344 |
|
|
ix = (ix << 16) | (i0 >> 16); |
345 |
|
|
if (ix < 0x3ffed800) /* |x|<0.84375 */ |
346 |
|
|
{ |
347 |
|
|
if (ix < 0x3fbe0000) /* |x|<2**-65 */ |
348 |
|
|
return one - x; |
349 |
|
|
z = x * x; |
350 |
|
|
r = pp[0] + z * (pp[1] |
351 |
|
|
+ z * (pp[2] + z * (pp[3] + z * (pp[4] + z * pp[5])))); |
352 |
|
|
s = qq[0] + z * (qq[1] |
353 |
|
|
+ z * (qq[2] + z * (qq[3] + z * (qq[4] + z * (qq[5] + z))))); |
354 |
|
|
y = r / s; |
355 |
|
|
if (ix < 0x3ffd8000) /* x<1/4 */ |
356 |
|
|
{ |
357 |
|
|
return one - (x + x * y); |
358 |
|
|
} |
359 |
|
|
else |
360 |
|
|
{ |
361 |
|
|
r = x * y; |
362 |
|
|
r += (x - half); |
363 |
|
|
return half - r; |
364 |
|
|
} |
365 |
|
|
} |
366 |
|
|
if (ix < 0x3fffa000) /* 1.25 */ |
367 |
|
|
{ /* 0.84375 <= |x| < 1.25 */ |
368 |
|
|
s = fabsl (x) - one; |
369 |
|
|
P = pa[0] + s * (pa[1] + s * (pa[2] |
370 |
|
|
+ s * (pa[3] + s * (pa[4] + s * (pa[5] + s * (pa[6] + s * pa[7])))))); |
371 |
|
|
Q = qa[0] + s * (qa[1] + s * (qa[2] |
372 |
|
|
+ s * (qa[3] + s * (qa[4] + s * (qa[5] + s * (qa[6] + s)))))); |
373 |
|
|
if ((se & 0x8000) == 0) |
374 |
|
|
{ |
375 |
|
|
z = one - erx; |
376 |
|
|
return z - P / Q; |
377 |
|
|
} |
378 |
|
|
else |
379 |
|
|
{ |
380 |
|
|
z = erx + P / Q; |
381 |
|
|
return one + z; |
382 |
|
|
} |
383 |
|
|
} |
384 |
|
|
if (ix < 0x4005d600) /* 107 */ |
385 |
|
|
{ /* |x|<107 */ |
386 |
|
|
x = fabsl (x); |
387 |
|
|
s = one / (x * x); |
388 |
|
|
if (ix < 0x4000b6db) /* 2.85711669921875 */ |
389 |
|
|
{ /* |x| < 1/.35 ~ 2.857143 */ |
390 |
|
|
R = ra[0] + s * (ra[1] + s * (ra[2] + s * (ra[3] + s * (ra[4] + |
391 |
|
|
s * (ra[5] + s * (ra[6] + s * (ra[7] + s * ra[8]))))))); |
392 |
|
|
S = sa[0] + s * (sa[1] + s * (sa[2] + s * (sa[3] + s * (sa[4] + |
393 |
|
|
s * (sa[5] + s * (sa[6] + s * (sa[7] + s * (sa[8] + s)))))))); |
394 |
|
|
} |
395 |
|
|
else if (ix < 0x4001d555) /* 6.6666259765625 */ |
396 |
|
|
{ /* 6.666 > |x| >= 1/.35 ~ 2.857143 */ |
397 |
|
|
R = rb[0] + s * (rb[1] + s * (rb[2] + s * (rb[3] + s * (rb[4] + |
398 |
|
|
s * (rb[5] + s * (rb[6] + s * rb[7])))))); |
399 |
|
|
S = sb[0] + s * (sb[1] + s * (sb[2] + s * (sb[3] + s * (sb[4] + |
400 |
|
|
s * (sb[5] + s * (sb[6] + s)))))); |
401 |
|
|
} |
402 |
|
|
else |
403 |
|
|
{ /* |x| >= 6.666 */ |
404 |
|
|
if (se & 0x8000) |
405 |
|
|
return two - tiny; /* x < -6.666 */ |
406 |
|
|
|
407 |
|
|
R = rc[0] + s * (rc[1] + s * (rc[2] + s * (rc[3] + |
408 |
|
|
s * (rc[4] + s * rc[5])))); |
409 |
|
|
S = sc[0] + s * (sc[1] + s * (sc[2] + s * (sc[3] + |
410 |
|
|
s * (sc[4] + s)))); |
411 |
|
|
} |
412 |
|
|
z = x; |
413 |
|
|
GET_LDOUBLE_WORDS (hx, i0, i1, z); |
414 |
|
|
i1 = 0; |
415 |
|
|
i0 &= 0xffffff00; |
416 |
|
|
SET_LDOUBLE_WORDS (z, hx, i0, i1); |
417 |
|
|
r = expl (-z * z - 0.5625) * |
418 |
|
|
expl ((z - x) * (z + x) + R / S); |
419 |
|
|
if ((se & 0x8000) == 0) |
420 |
|
|
return r / x; |
421 |
|
|
else |
422 |
|
|
return two - r / x; |
423 |
|
|
} |
424 |
|
|
else |
425 |
|
|
{ |
426 |
|
|
if ((se & 0x8000) == 0) |
427 |
|
|
return tiny * tiny; |
428 |
|
|
else |
429 |
|
|
return two - tiny; |
430 |
|
|
} |
431 |
|
|
} |
432 |
|
|
DEF_STD(erfcl); |