5 * Common Development and Distribution License (the "License").
6 * You may not use this file except in compliance with the License.
7 *
8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9 * or http://www.opensolaris.org/os/licensing.
10 * See the License for the specific language governing permissions
11 * and limitations under the License.
12 *
13 * When distributing Covered Code, include this CDDL HEADER in each
14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15 * If applicable, add the following below this CDDL HEADER, with the
16 * fields enclosed by brackets "[]" replaced with your own identifying
17 * information: Portions Copyright [yyyy] [name of copyright owner]
18 *
19 * CDDL HEADER END
20 */
21
22 /*
23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
24 */
25 /*
26 * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
27 * Use is subject to license terms.
28 */
29
30 /* long double function erf,erfc (long double x)
31 * K.C. Ng, September, 1989.
32 * x
33 * 2 |\
34 * erf(x) = --------- | exp(-t*t)dt
35 * sqrt(pi) \|
36 * 0
37 *
38 * erfc(x) = 1-erf(x)
39 *
40 * method:
41 * Since erf(-x) = -erf(x), we assume x>=0.
42 * For x near 0, we have the expansion
43 *
44 * erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....).
45 *
46 * Since 2/sqrt(pi) = 1.128379167095512573896158903121545171688,
47 * we use x + x*P(x^2) to approximate erf(x). This formula will
48 * guarantee the error less than one ulp where x is not too far
49 * away from 0. We note that erf(x)=x at x = 0.6174...... After
50 * some experiment, we choose the following approximation on
74 *
75 *
76 * For x in [1.75,16/3]
77 * erfc(x) = exp(-x*x)*(1/x)*R1(1/x)/S1(1/x)
78 * erf(x) = 1 - erfc(x)
79 * precision: absolute error of R1/S1 is bounded by 2**-124.03
80 *
81 * For x in [16/3,107]
82 * erfc(x) = exp(-x*x)*(1/x)*R2(1/x)/S2(1/x)
83 * erf(x) = 1 - erfc(x) (if x>=9 simple return erf(x)=1 with inexact)
84 * precision: absolute error of R2/S2 is bounded by 2**-120.07
85 *
86 * Else if inf > x >= 107
87 * erf(x) = 1 with inexact
88 * erfc(x) = 0 with underflow
89 *
90 * Special case:
91 * erf(inf) = 1
92 * erfc(inf) = 0
93 */
94
95 #pragma weak __erfl = erfl
96 #pragma weak __erfcl = erfcl
97
98 #include "libm.h"
99 #include "longdouble.h"
100
101 static long double
102 tiny = 1e-40L,
103 nearunfl = 1e-4000L,
104 half = 0.5L,
105 one = 1.0L,
106 onehalf = 1.5L,
107 L16_3 = 16.0L/3.0L;
108 /*
109 * Coefficients for even polynomial P for erf(x)=x+x*P(x^2) on [0,0.84375]
110 */
111 static long double P[] = { /* 21 coeffs */
112 1.283791670955125738961589031215451715556e-0001L,
113 -3.761263890318375246320529677071815594603e-0001L,
114 1.128379167095512573896158903121205899135e-0001L,
115 -2.686617064513125175943235483344625046092e-0002L,
116 5.223977625442187842111846652980454568389e-0003L,
117 -8.548327023450852832546626271083862724358e-0004L,
118 1.205533298178966425102164715902231976672e-0004L,
119 -1.492565035840625097674944905027897838996e-0005L,
120 1.646211436588924733604648849172936692024e-0006L,
121 -1.636584469123491976815834704799733514987e-0007L,
122 1.480719281587897445302529007144770739305e-0008L,
123 -1.229055530170782843046467986464722047175e-0009L,
124 9.422759064320307357553954945760654341633e-0011L,
125 -6.711366846653439036162105104991433380926e-0012L,
126 4.463224090341893165100275380693843116240e-0013L,
127 -2.783513452582658245422635662559779162312e-0014L,
133 };
134
135 /*
136 * Rational erf(x) = ((float)0.84506291151) + P1(x-1)/Q1(x-1) on [0.84375,1.25]
137 */
138 static long double C1 = (long double)((float)0.84506291151);
139 static long double P1[] = { /* 12 top coeffs */
140 -2.362118560752659955654364917390741930316e-0003L,
141 4.129623379624420034078926610650759979146e-0001L,
142 -3.973857505403547283109417923182669976904e-0002L,
143 4.357503184084022439763567513078036755183e-0002L,
144 8.015593623388421371247676683754171456950e-0002L,
145 -1.034459310403352486685467221776778474602e-0002L,
146 5.671850295381046679675355719017720821383e-0003L,
147 1.219262563232763998351452194968781174318e-0003L,
148 5.390833481581033423020320734201065475098e-0004L,
149 -1.978853912815115495053119023517805528300e-0004L,
150 6.184234513953600118335017885706420552487e-0005L,
151 -5.331802711697810861017518515816271808286e-0006L,
152 };
153 static long double Q1[] = { /* 12 bottom coeffs with leading 1.0 hidden */
154 9.081506296064882195280178373107623196655e-0001L,
155 6.821049531968204097604392183650687642520e-0001L,
156 4.067869178233539502315055970743271822838e-0001L,
157 1.702332233546316765818144723063881095577e-0001L,
158 7.498098377690553934266423088708614219356e-0002L,
159 2.050154396918178697056927234366372760310e-0002L,
160 7.012988534031999899054782333851905939379e-0003L,
161 1.149904787014400354649843451234570731076e-0003L,
162 3.185620255011299476196039491205159718620e-0004L,
163 1.273405072153008775426376193374105840517e-0005L,
164 4.753866999959432971956781228148402971454e-0006L,
165 -1.002287602111660026053981728549540200683e-0006L,
166 };
167 /*
168 * Rational erf(x) = ((float)0.95478588343) + P2(x-1.5)/Q2(x-1.5)
169 * on [1.25,1.75]
170 */
171 static long double C2 = (long double)((float)0.95478588343);
172 static long double P2[] = { /* 12 top coeffs */
173 1.131926304864446730135126164594785863512e-0002L,
174 1.273617996967754151544330055186210322832e-0001L,
175 -8.169980734667512519897816907190281143423e-0002L,
176 9.512267486090321197833634271787944271746e-0002L,
177 -2.394251569804872160005274999735914368170e-0002L,
178 1.108768660227528667525252333184520222905e-0002L,
179 3.527435492933902414662043314373277494221e-0004L,
180 4.946116273341953463584319006669474625971e-0004L,
181 -4.289851942513144714600285769022420962418e-0005L,
182 8.304719841341952705874781636002085119978e-0005L,
183 -1.040460226177309338781902252282849903189e-0005L,
184 2.122913331584921470381327583672044434087e-0006L,
185 };
186 static long double Q2[] = { /* 13 bottom coeffs with leading 1.0 hidden */
187 7.448815737306992749168727691042003832150e-0001L,
188 7.161813850236008294484744312430122188043e-0001L,
189 3.603134756584225766144922727405641236121e-0001L,
190 1.955811609133766478080550795194535852653e-0001L,
191 7.253059963716225972479693813787810711233e-0002L,
192 2.752391253757421424212770221541238324978e-0002L,
193 7.677654852085240257439050673446546828005e-0003L,
194 2.141102244555509687346497060326630061069e-0003L,
195 4.342123013830957093949563339130674364271e-0004L,
196 8.664587895570043348530991997272212150316e-0005L,
197 1.109201582511752087060167429397033701988e-0005L,
198 1.357834375781831062713347000030984364311e-0006L,
199 4.957746280594384997273090385060680016451e-0008L,
200 };
201 /*
202 * erfc(x) = exp(-x*x)/x * R1(1/x)/S1(1/x) on [1.75, 16/3]
203 */
204 static long double R1[] = { /* 14 top coeffs */
205 4.630195122654315016370705767621550602948e+0006L,
206 1.257949521746494830700654204488675713628e+0007L,
207 1.704153822720260272814743497376181625707e+0007L,
208 1.502600568706061872381577539537315739943e+0007L,
209 9.543710793431995284827024445387333922861e+0006L,
210 4.589344808584091011652238164935949522427e+0006L,
211 1.714660662941745791190907071920671844289e+0006L,
212 5.034802147768798894307672256192466283867e+0005L,
213 1.162286400443554670553152110447126850725e+0005L,
214 2.086643834548901681362757308058660399137e+0004L,
215 2.839793161868140305907004392890348777338e+0003L,
216 2.786687241658423601778258694498655680778e+0002L,
217 1.779177837102695602425897452623985786464e+0001L,
218 5.641895835477470769043614623819144434731e-0001L,
219 };
220 static long double S1[] = { /* 15 bottom coeffs with leading 1.0 hidden */
221 4.630195122654331529595606896287596843110e+0006L,
222 1.780411093345512024324781084220509055058e+0007L,
223 3.250113097051800703707108623715776848283e+0007L,
224 3.737857099176755050912193712123489115755e+0007L,
225 3.029787497516578821459174055870781168593e+0007L,
226 1.833850619965384765005769632103205777227e+0007L,
227 8.562719999736915722210391222639186586498e+0006L,
228 3.139684562074658971315545539760008136973e+0006L,
229 9.106421313731384880027703627454366930945e+0005L,
230 2.085108342384266508613267136003194920001e+0005L,
231 3.723126272693120340730491416449539290600e+0004L,
232 5.049169878567344046145695360784436929802e+0003L,
233 4.944274532748010767670150730035392093899e+0002L,
234 3.153510608818213929982940249162268971412e+0001L,
235 1.0e00L,
236 };
237
238 /*
239 * erfc(x) = exp(-x*x)/x * R2(1/x)/S2(1/x) on [16/3, 107]
240 */
241 static long double R2[] = { /* 15 top coeffs in reverse order!!*/
242 2.447288012254302966796326587537136931669e+0005L,
243 8.768592567189861896653369912716538739016e+0005L,
244 1.552293152581780065761497908005779524953e+0006L,
245 1.792075924835942935864231657504259926729e+0006L,
246 1.504001463155897344947500222052694835875e+0006L,
247 9.699485556326891411801230186016013019935e+0005L,
248 4.961449933661807969863435013364796037700e+0005L,
249 2.048726544693474028061176764716228273791e+0005L,
250 6.891532964330949722479061090551896886635e+0004L,
251 1.888014709010307507771964047905823237985e+0004L,
252 4.189692064988957745054734809642495644502e+0003L,
253 7.362346487427048068212968889642741734621e+0002L,
254 9.980359714211411423007641056580813116207e+0001L,
255 9.426910895135379181107191962193485174159e+0000L,
256 5.641895835477562869480794515623601280429e-0001L,
257 };
258 static long double S2[] = { /* 16 coefficients */
259 2.447282203601902971246004716790604686880e+0005L,
260 1.153009852759385309367759460934808489833e+0006L,
261 2.608580649612639131548966265078663384849e+0006L,
262 3.766673917346623308850202792390569025740e+0006L,
263 3.890566255138383910789924920541335370691e+0006L,
264 3.052882073900746207613166259994150527732e+0006L,
265 1.885574519970380988460241047248519418407e+0006L,
266 9.369722034759943185851450846811445012922e+0005L,
267 3.792278350536686111444869752624492443659e+0005L,
268 1.257750606950115799965366001773094058720e+0005L,
269 3.410830600242369370645608634643620355058e+0004L,
270 7.513984469742343134851326863175067271240e+0003L,
271 1.313296320593190002554779998138695507840e+0003L,
272 1.773972700887629157006326333696896516769e+0002L,
273 1.670876451822586800422009013880457094162e+0001L,
274 1.000L,
275 };
276
277 long double erfl(x)
278 long double x;
279 {
280 long double erfcl(long double),s,y,t;
281
282 if (!finitel(x)) {
283 if (x != x) return x+x; /* NaN */
284 return copysignl(one,x); /* return +-1.0 is x=Inf */
285 }
286
287 y = fabsl(x);
288 if (y <= 0.84375L) {
289 if (y<=tiny) return x+P[0]*x;
290 s = y*y;
291 t = __poly_libmq(s,21,P);
292 return x+x*t;
293 }
294 if (y<=1.25L) {
295 s = y-one;
296 t = C1+__poly_libmq(s,12,P1)/(one+s*__poly_libmq(s,12,Q1));
297 return (signbitl(x))? -t: t;
298 } else if (y<=1.75L) {
299 s = y-onehalf;
300 t = C2+__poly_libmq(s,12,P2)/(one+s*__poly_libmq(s,13,Q2));
301 return (signbitl(x))? -t: t;
302 }
303 if (y<=9.0L) t = erfcl(y); else t = tiny;
304 return (signbitl(x))? t-one: one-t;
305 }
306
307 long double erfcl(x)
308 long double x;
309 {
310 long double erfl(long double),s,y,t;
311
312 if (!finitel(x)) {
313 if (x != x) return x+x; /* NaN */
314 /* return 2.0 if x= -inf
315 0.0 if x= +inf */
316 if (x<0.0L) return 2.0L; else return 0.0L;
317 }
318
319 if (x <= 0.84375L) {
320 if (x<=0.25) return one-erfl(x);
321 s = x*x;
322 t = half-x;
323 t = t - x*__poly_libmq(s,21,P);
324 return half+t;
325 }
326 if (x<=1.25L) {
327 s = x-one;
328 t = one-C1;
329 return t - __poly_libmq(s,12,P1)/(one+s*__poly_libmq(s,12,Q1));
330 } else if (x<=1.75L) {
331 s = x-onehalf;
332 t = one-C2;
333 return t - __poly_libmq(s,12,P2)/(one+s*__poly_libmq(s,13,Q2));
334 }
335 if (x>=107.0L) return nearunfl*nearunfl; /* underflow */
336 else if (x >= L16_3) {
337 y = __poly_libmq(x,15,R2);
338 t = y/__poly_libmq(x,16,S2);
339 } else {
340 y = __poly_libmq(x,14,R1);
341 t = y/__poly_libmq(x,15,S1);
342 }
343 /* see comment in ../Q/erfl.c */
344 y = x;
345 *(int*)&y = 0;
346 t *= expl(-y*y)*expl(-(x-y)*(x+y));
347 return t;
348 }
|
5 * Common Development and Distribution License (the "License").
6 * You may not use this file except in compliance with the License.
7 *
8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9 * or http://www.opensolaris.org/os/licensing.
10 * See the License for the specific language governing permissions
11 * and limitations under the License.
12 *
13 * When distributing Covered Code, include this CDDL HEADER in each
14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15 * If applicable, add the following below this CDDL HEADER, with the
16 * fields enclosed by brackets "[]" replaced with your own identifying
17 * information: Portions Copyright [yyyy] [name of copyright owner]
18 *
19 * CDDL HEADER END
20 */
21
22 /*
23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
24 */
25
26 /*
27 * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
28 * Use is subject to license terms.
29 */
30
31 /* BEGIN CSTYLED */
32 /*
33 * long double function erf,erfc (long double x)
34 * K.C. Ng, September, 1989.
35 * x
36 * 2 |\
37 * erf(x) = --------- | exp(-t*t)dt
38 * sqrt(pi) \|
39 * 0
40 *
41 * erfc(x) = 1-erf(x)
42 *
43 * method:
44 * Since erf(-x) = -erf(x), we assume x>=0.
45 * For x near 0, we have the expansion
46 *
47 * erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....).
48 *
49 * Since 2/sqrt(pi) = 1.128379167095512573896158903121545171688,
50 * we use x + x*P(x^2) to approximate erf(x). This formula will
51 * guarantee the error less than one ulp where x is not too far
52 * away from 0. We note that erf(x)=x at x = 0.6174...... After
53 * some experiment, we choose the following approximation on
77 *
78 *
79 * For x in [1.75,16/3]
80 * erfc(x) = exp(-x*x)*(1/x)*R1(1/x)/S1(1/x)
81 * erf(x) = 1 - erfc(x)
82 * precision: absolute error of R1/S1 is bounded by 2**-124.03
83 *
84 * For x in [16/3,107]
85 * erfc(x) = exp(-x*x)*(1/x)*R2(1/x)/S2(1/x)
86 * erf(x) = 1 - erfc(x) (if x>=9 simple return erf(x)=1 with inexact)
87 * precision: absolute error of R2/S2 is bounded by 2**-120.07
88 *
89 * Else if inf > x >= 107
90 * erf(x) = 1 with inexact
91 * erfc(x) = 0 with underflow
92 *
93 * Special case:
94 * erf(inf) = 1
95 * erfc(inf) = 0
96 */
97 /* END CSTYLED */
98
99 #pragma weak __erfl = erfl
100 #pragma weak __erfcl = erfcl
101
102 #include "libm.h"
103 #include "longdouble.h"
104
105 static long double tiny = 1e-40L,
106 nearunfl = 1e-4000L,
107 half = 0.5L,
108 one = 1.0L,
109 onehalf = 1.5L,
110 L16_3 = 16.0L / 3.0L;
111
112 /*
113 * Coefficients for even polynomial P for erf(x)=x+x*P(x^2) on [0,0.84375]
114 */
115 static long double P[] = { /* 21 coeffs */
116 1.283791670955125738961589031215451715556e-0001L,
117 -3.761263890318375246320529677071815594603e-0001L,
118 1.128379167095512573896158903121205899135e-0001L,
119 -2.686617064513125175943235483344625046092e-0002L,
120 5.223977625442187842111846652980454568389e-0003L,
121 -8.548327023450852832546626271083862724358e-0004L,
122 1.205533298178966425102164715902231976672e-0004L,
123 -1.492565035840625097674944905027897838996e-0005L,
124 1.646211436588924733604648849172936692024e-0006L,
125 -1.636584469123491976815834704799733514987e-0007L,
126 1.480719281587897445302529007144770739305e-0008L,
127 -1.229055530170782843046467986464722047175e-0009L,
128 9.422759064320307357553954945760654341633e-0011L,
129 -6.711366846653439036162105104991433380926e-0012L,
130 4.463224090341893165100275380693843116240e-0013L,
131 -2.783513452582658245422635662559779162312e-0014L,
137 };
138
139 /*
140 * Rational erf(x) = ((float)0.84506291151) + P1(x-1)/Q1(x-1) on [0.84375,1.25]
141 */
142 static long double C1 = (long double)((float)0.84506291151);
143 static long double P1[] = { /* 12 top coeffs */
144 -2.362118560752659955654364917390741930316e-0003L,
145 4.129623379624420034078926610650759979146e-0001L,
146 -3.973857505403547283109417923182669976904e-0002L,
147 4.357503184084022439763567513078036755183e-0002L,
148 8.015593623388421371247676683754171456950e-0002L,
149 -1.034459310403352486685467221776778474602e-0002L,
150 5.671850295381046679675355719017720821383e-0003L,
151 1.219262563232763998351452194968781174318e-0003L,
152 5.390833481581033423020320734201065475098e-0004L,
153 -1.978853912815115495053119023517805528300e-0004L,
154 6.184234513953600118335017885706420552487e-0005L,
155 -5.331802711697810861017518515816271808286e-0006L,
156 };
157
158 static long double Q1[] = { /* 12 bottom coeffs with leading 1.0 hidden */
159 9.081506296064882195280178373107623196655e-0001L,
160 6.821049531968204097604392183650687642520e-0001L,
161 4.067869178233539502315055970743271822838e-0001L,
162 1.702332233546316765818144723063881095577e-0001L,
163 7.498098377690553934266423088708614219356e-0002L,
164 2.050154396918178697056927234366372760310e-0002L,
165 7.012988534031999899054782333851905939379e-0003L,
166 1.149904787014400354649843451234570731076e-0003L,
167 3.185620255011299476196039491205159718620e-0004L,
168 1.273405072153008775426376193374105840517e-0005L,
169 4.753866999959432971956781228148402971454e-0006L,
170 -1.002287602111660026053981728549540200683e-0006L,
171 };
172
173 /*
174 * Rational erf(x) = ((float)0.95478588343) + P2(x-1.5)/Q2(x-1.5)
175 * on [1.25,1.75]
176 */
177 static long double C2 = (long double)((float)0.95478588343);
178 static long double P2[] = { /* 12 top coeffs */
179 1.131926304864446730135126164594785863512e-0002L,
180 1.273617996967754151544330055186210322832e-0001L,
181 -8.169980734667512519897816907190281143423e-0002L,
182 9.512267486090321197833634271787944271746e-0002L,
183 -2.394251569804872160005274999735914368170e-0002L,
184 1.108768660227528667525252333184520222905e-0002L,
185 3.527435492933902414662043314373277494221e-0004L,
186 4.946116273341953463584319006669474625971e-0004L,
187 -4.289851942513144714600285769022420962418e-0005L,
188 8.304719841341952705874781636002085119978e-0005L,
189 -1.040460226177309338781902252282849903189e-0005L,
190 2.122913331584921470381327583672044434087e-0006L,
191 };
192
193 static long double Q2[] = { /* 13 bottom coeffs with leading 1.0 hidden */
194 7.448815737306992749168727691042003832150e-0001L,
195 7.161813850236008294484744312430122188043e-0001L,
196 3.603134756584225766144922727405641236121e-0001L,
197 1.955811609133766478080550795194535852653e-0001L,
198 7.253059963716225972479693813787810711233e-0002L,
199 2.752391253757421424212770221541238324978e-0002L,
200 7.677654852085240257439050673446546828005e-0003L,
201 2.141102244555509687346497060326630061069e-0003L,
202 4.342123013830957093949563339130674364271e-0004L,
203 8.664587895570043348530991997272212150316e-0005L,
204 1.109201582511752087060167429397033701988e-0005L,
205 1.357834375781831062713347000030984364311e-0006L,
206 4.957746280594384997273090385060680016451e-0008L,
207 };
208
209 /*
210 * erfc(x) = exp(-x*x)/x * R1(1/x)/S1(1/x) on [1.75, 16/3]
211 */
212 static long double R1[] = { /* 14 top coeffs */
213 4.630195122654315016370705767621550602948e+0006L,
214 1.257949521746494830700654204488675713628e+0007L,
215 1.704153822720260272814743497376181625707e+0007L,
216 1.502600568706061872381577539537315739943e+0007L,
217 9.543710793431995284827024445387333922861e+0006L,
218 4.589344808584091011652238164935949522427e+0006L,
219 1.714660662941745791190907071920671844289e+0006L,
220 5.034802147768798894307672256192466283867e+0005L,
221 1.162286400443554670553152110447126850725e+0005L,
222 2.086643834548901681362757308058660399137e+0004L,
223 2.839793161868140305907004392890348777338e+0003L,
224 2.786687241658423601778258694498655680778e+0002L,
225 1.779177837102695602425897452623985786464e+0001L,
226 5.641895835477470769043614623819144434731e-0001L,
227 };
228
229 static long double S1[] = { /* 15 bottom coeffs with leading 1.0 hidden */
230 4.630195122654331529595606896287596843110e+0006L,
231 1.780411093345512024324781084220509055058e+0007L,
232 3.250113097051800703707108623715776848283e+0007L,
233 3.737857099176755050912193712123489115755e+0007L,
234 3.029787497516578821459174055870781168593e+0007L,
235 1.833850619965384765005769632103205777227e+0007L,
236 8.562719999736915722210391222639186586498e+0006L,
237 3.139684562074658971315545539760008136973e+0006L,
238 9.106421313731384880027703627454366930945e+0005L,
239 2.085108342384266508613267136003194920001e+0005L,
240 3.723126272693120340730491416449539290600e+0004L,
241 5.049169878567344046145695360784436929802e+0003L,
242 4.944274532748010767670150730035392093899e+0002L,
243 3.153510608818213929982940249162268971412e+0001L,
244 1.0e00L,
245 };
246
247 /*
248 * erfc(x) = exp(-x*x)/x * R2(1/x)/S2(1/x) on [16/3, 107]
249 */
250 static long double R2[] = { /* 15 top coeffs in reverse order!! */
251 2.447288012254302966796326587537136931669e+0005L,
252 8.768592567189861896653369912716538739016e+0005L,
253 1.552293152581780065761497908005779524953e+0006L,
254 1.792075924835942935864231657504259926729e+0006L,
255 1.504001463155897344947500222052694835875e+0006L,
256 9.699485556326891411801230186016013019935e+0005L,
257 4.961449933661807969863435013364796037700e+0005L,
258 2.048726544693474028061176764716228273791e+0005L,
259 6.891532964330949722479061090551896886635e+0004L,
260 1.888014709010307507771964047905823237985e+0004L,
261 4.189692064988957745054734809642495644502e+0003L,
262 7.362346487427048068212968889642741734621e+0002L,
263 9.980359714211411423007641056580813116207e+0001L,
264 9.426910895135379181107191962193485174159e+0000L,
265 5.641895835477562869480794515623601280429e-0001L,
266 };
267
268 static long double S2[] = { /* 16 coefficients */
269 2.447282203601902971246004716790604686880e+0005L,
270 1.153009852759385309367759460934808489833e+0006L,
271 2.608580649612639131548966265078663384849e+0006L,
272 3.766673917346623308850202792390569025740e+0006L,
273 3.890566255138383910789924920541335370691e+0006L,
274 3.052882073900746207613166259994150527732e+0006L,
275 1.885574519970380988460241047248519418407e+0006L,
276 9.369722034759943185851450846811445012922e+0005L,
277 3.792278350536686111444869752624492443659e+0005L,
278 1.257750606950115799965366001773094058720e+0005L,
279 3.410830600242369370645608634643620355058e+0004L,
280 7.513984469742343134851326863175067271240e+0003L,
281 1.313296320593190002554779998138695507840e+0003L,
282 1.773972700887629157006326333696896516769e+0002L,
283 1.670876451822586800422009013880457094162e+0001L,
284 1.000L,
285 };
286
287 long double
288 erfl(long double x)
289 {
290 long double erfcl(long double), s, y, t;
291
292 if (!finitel(x)) {
293 if (x != x)
294 return (x + x); /* NaN */
295
296 return (copysignl(one, x)); /* return +-1.0 is x=Inf */
297 }
298
299 y = fabsl(x);
300
301 if (y <= 0.84375L) {
302 if (y <= tiny)
303 return (x + P[0] * x);
304
305 s = y * y;
306 t = __poly_libmq(s, 21, P);
307 return (x + x * t);
308 }
309
310 if (y <= 1.25L) {
311 s = y - one;
312 t = C1 + __poly_libmq(s, 12, P1) / (one + s * __poly_libmq(s,
313 12, Q1));
314 return ((signbitl(x)) ? -t : t);
315 } else if (y <= 1.75L) {
316 s = y - onehalf;
317 t = C2 + __poly_libmq(s, 12, P2) / (one + s * __poly_libmq(s,
318 13, Q2));
319 return ((signbitl(x)) ? -t : t);
320 }
321
322 if (y <= 9.0L)
323 t = erfcl(y);
324 else
325 t = tiny;
326
327 return ((signbitl(x)) ? t - one : one - t);
328 }
329
330 long double
331 erfcl(long double x)
332 {
333 long double s, y, t;
334
335 if (!finitel(x)) {
336 if (x != x)
337 return (x + x); /* NaN */
338
339 /*
340 * return 2.0 if x = -inf
341 * 0.0 if x = +inf
342 */
343 if (x < 0.0L)
344 return (2.0L);
345 else
346 return (0.0L);
347 }
348
349 if (x <= 0.84375L) {
350 if (x <= 0.25)
351 return (one - erfl(x));
352
353 s = x * x;
354 t = half - x;
355 t = t - x * __poly_libmq(s, 21, P);
356 return (half + t);
357 }
358
359 if (x <= 1.25L) {
360 s = x - one;
361 t = one - C1;
362 return (t - __poly_libmq(s, 12, P1) / (one + s * __poly_libmq(s,
363 12, Q1)));
364 } else if (x <= 1.75L) {
365 s = x - onehalf;
366 t = one - C2;
367 return (t - __poly_libmq(s, 12, P2) / (one + s * __poly_libmq(s,
368 13, Q2)));
369 }
370
371 if (x >= 107.0L) {
372 return (nearunfl * nearunfl); /* underflow */
373 } else if (x >= L16_3) {
374 y = __poly_libmq(x, 15, R2);
375 t = y / __poly_libmq(x, 16, S2);
376 } else {
377 y = __poly_libmq(x, 14, R1);
378 t = y / __poly_libmq(x, 15, S1);
379 }
380
381 /* see comment in ../Q/erfl.c */
382 y = x;
383 *(int *)&y = 0;
384 t *= expl(-y * y) * expl(-(x - y) * (x + y));
385 return (t);
386 }
|