]>
Commit | Line | Data |
---|---|---|
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 | /* Modifications and expansions for 128-bit long double 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 | <https://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. erf(x) = x + x*R(x^2) for |x| in [0, 7/8] | |
48 | * Remark. The formula is derived by noting | |
49 | * erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....) | |
50 | * and that | |
51 | * 2/sqrt(pi) = 1.128379167095512573896158903121545171688 | |
52 | * is close to one. | |
53 | * | |
54 | * 1a. erf(x) = 1 - erfc(x), for |x| > 1.0 | |
55 | * erfc(x) = 1 - erf(x) if |x| < 1/4 | |
56 | * | |
57 | * 2. For |x| in [7/8, 1], let s = |x| - 1, and | |
58 | * c = 0.84506291151 rounded to single (24 bits) | |
59 | * erf(s + c) = sign(x) * (c + P1(s)/Q1(s)) | |
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/4, 5/4], | |
66 | * erfc(s + const) = erfc(const) + s P1(s)/Q1(s) | |
67 | * for const = 1/4, 3/8, ..., 9/8 | |
68 | * and 0 <= s <= 1/8 . | |
69 | * | |
70 | * 4. For x in [5/4, 107], | |
71 | * erfc(x) = (1/x)*exp(-x*x-0.5625 + R(z)) | |
72 | * z=1/x^2 | |
73 | * The interval is partitioned into several segments | |
74 | * of width 1/8 in 1/x. | |
75 | * | |
76 | * Note1: | |
77 | * To compute exp(-x*x-0.5625+R/S), let s be a single | |
78 | * precision number and s := x; then | |
79 | * -x*x = -s*s + (s-x)*(s+x) | |
80 | * exp(-x*x-0.5626+R/S) = | |
81 | * exp(-s*s-0.5625)*exp((s-x)*(s+x)+R/S); | |
82 | * Note2: | |
83 | * Here 4 and 5 make use of the asymptotic series | |
84 | * exp(-x*x) | |
85 | * erfc(x) ~ ---------- * ( 1 + Poly(1/x^2) ) | |
86 | * x*sqrt(pi) | |
87 | * | |
88 | * 5. For inf > x >= 107 | |
89 | * erf(x) = sign(x) *(1 - tiny) (raise inexact) | |
90 | * erfc(x) = tiny*tiny (raise underflow) if x > 0 | |
91 | * = 2 - tiny if x<0 | |
92 | * | |
93 | * 7. Special case: | |
94 | * erf(0) = 0, erf(inf) = 1, erf(-inf) = -1, | |
95 | * erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2, | |
96 | * erfc/erf(NaN) is NaN | |
97 | */ | |
98 | ||
99 | #include <errno.h> | |
100 | #include <float.h> | |
101 | #include <math.h> | |
102 | #include <math_private.h> | |
103 | #include <math-underflow.h> | |
104 | #include <libm-alias-ldouble.h> | |
105 | ||
106 | /* Evaluate P[n] x^n + P[n-1] x^(n-1) + ... + P[0] */ | |
107 | ||
108 | static _Float128 | |
109 | neval (_Float128 x, const _Float128 *p, int n) | |
110 | { | |
111 | _Float128 y; | |
112 | ||
113 | p += n; | |
114 | y = *p--; | |
115 | do | |
116 | { | |
117 | y = y * x + *p--; | |
118 | } | |
119 | while (--n > 0); | |
120 | return y; | |
121 | } | |
122 | ||
123 | ||
124 | /* Evaluate x^n+1 + P[n] x^(n) + P[n-1] x^(n-1) + ... + P[0] */ | |
125 | ||
126 | static _Float128 | |
127 | deval (_Float128 x, const _Float128 *p, int n) | |
128 | { | |
129 | _Float128 y; | |
130 | ||
131 | p += n; | |
132 | y = x + *p--; | |
133 | do | |
134 | { | |
135 | y = y * x + *p--; | |
136 | } | |
137 | while (--n > 0); | |
138 | return y; | |
139 | } | |
140 | ||
141 | ||
142 | ||
143 | static const _Float128 | |
144 | tiny = L(1e-4931), | |
145 | one = 1, | |
146 | two = 2, | |
147 | /* 2/sqrt(pi) - 1 */ | |
148 | efx = L(1.2837916709551257389615890312154517168810E-1); | |
149 | ||
150 | ||
151 | /* erf(x) = x + x R(x^2) | |
152 | 0 <= x <= 7/8 | |
153 | Peak relative error 1.8e-35 */ | |
154 | #define NTN1 8 | |
155 | static const _Float128 TN1[NTN1 + 1] = | |
156 | { | |
157 | L(-3.858252324254637124543172907442106422373E10), | |
158 | L(9.580319248590464682316366876952214879858E10), | |
159 | L(1.302170519734879977595901236693040544854E10), | |
160 | L(2.922956950426397417800321486727032845006E9), | |
161 | L(1.764317520783319397868923218385468729799E8), | |
162 | L(1.573436014601118630105796794840834145120E7), | |
163 | L(4.028077380105721388745632295157816229289E5), | |
164 | L(1.644056806467289066852135096352853491530E4), | |
165 | L(3.390868480059991640235675479463287886081E1) | |
166 | }; | |
167 | #define NTD1 8 | |
168 | static const _Float128 TD1[NTD1 + 1] = | |
169 | { | |
170 | L(-3.005357030696532927149885530689529032152E11), | |
171 | L(-1.342602283126282827411658673839982164042E11), | |
172 | L(-2.777153893355340961288511024443668743399E10), | |
173 | L(-3.483826391033531996955620074072768276974E9), | |
174 | L(-2.906321047071299585682722511260895227921E8), | |
175 | L(-1.653347985722154162439387878512427542691E7), | |
176 | L(-6.245520581562848778466500301865173123136E5), | |
177 | L(-1.402124304177498828590239373389110545142E4), | |
178 | L(-1.209368072473510674493129989468348633579E2) | |
179 | /* 1.0E0 */ | |
180 | }; | |
181 | ||
182 | ||
183 | /* erf(z+1) = erf_const + P(z)/Q(z) | |
184 | -.125 <= z <= 0 | |
185 | Peak relative error 7.3e-36 */ | |
186 | static const _Float128 erf_const = L(0.845062911510467529296875); | |
187 | #define NTN2 8 | |
188 | static const _Float128 TN2[NTN2 + 1] = | |
189 | { | |
190 | L(-4.088889697077485301010486931817357000235E1), | |
191 | L(7.157046430681808553842307502826960051036E3), | |
192 | L(-2.191561912574409865550015485451373731780E3), | |
193 | L(2.180174916555316874988981177654057337219E3), | |
194 | L(2.848578658049670668231333682379720943455E2), | |
195 | L(1.630362490952512836762810462174798925274E2), | |
196 | L(6.317712353961866974143739396865293596895E0), | |
197 | L(2.450441034183492434655586496522857578066E1), | |
198 | L(5.127662277706787664956025545897050896203E-1) | |
199 | }; | |
200 | #define NTD2 8 | |
201 | static const _Float128 TD2[NTD2 + 1] = | |
202 | { | |
203 | L(1.731026445926834008273768924015161048885E4), | |
204 | L(1.209682239007990370796112604286048173750E4), | |
205 | L(1.160950290217993641320602282462976163857E4), | |
206 | L(5.394294645127126577825507169061355698157E3), | |
207 | L(2.791239340533632669442158497532521776093E3), | |
208 | L(8.989365571337319032943005387378993827684E2), | |
209 | L(2.974016493766349409725385710897298069677E2), | |
210 | L(6.148192754590376378740261072533527271947E1), | |
211 | L(1.178502892490738445655468927408440847480E1) | |
212 | /* 1.0E0 */ | |
213 | }; | |
214 | ||
215 | ||
216 | /* erfc(x + 0.25) = erfc(0.25) + x R(x) | |
217 | 0 <= x < 0.125 | |
218 | Peak relative error 1.4e-35 */ | |
219 | #define NRNr13 8 | |
220 | static const _Float128 RNr13[NRNr13 + 1] = | |
221 | { | |
222 | L(-2.353707097641280550282633036456457014829E3), | |
223 | L(3.871159656228743599994116143079870279866E2), | |
224 | L(-3.888105134258266192210485617504098426679E2), | |
225 | L(-2.129998539120061668038806696199343094971E1), | |
226 | L(-8.125462263594034672468446317145384108734E1), | |
227 | L(8.151549093983505810118308635926270319660E0), | |
228 | L(-5.033362032729207310462422357772568553670E0), | |
229 | L(-4.253956621135136090295893547735851168471E-2), | |
230 | L(-8.098602878463854789780108161581050357814E-2) | |
231 | }; | |
232 | #define NRDr13 7 | |
233 | static const _Float128 RDr13[NRDr13 + 1] = | |
234 | { | |
235 | L(2.220448796306693503549505450626652881752E3), | |
236 | L(1.899133258779578688791041599040951431383E2), | |
237 | L(1.061906712284961110196427571557149268454E3), | |
238 | L(7.497086072306967965180978101974566760042E1), | |
239 | L(2.146796115662672795876463568170441327274E2), | |
240 | L(1.120156008362573736664338015952284925592E1), | |
241 | L(2.211014952075052616409845051695042741074E1), | |
242 | L(6.469655675326150785692908453094054988938E-1) | |
243 | /* 1.0E0 */ | |
244 | }; | |
245 | /* erfc(0.25) = C13a + C13b to extra precision. */ | |
246 | static const _Float128 C13a = L(0.723663330078125); | |
247 | static const _Float128 C13b = L(1.0279753638067014931732235184287934646022E-5); | |
248 | ||
249 | ||
250 | /* erfc(x + 0.375) = erfc(0.375) + x R(x) | |
251 | 0 <= x < 0.125 | |
252 | Peak relative error 1.2e-35 */ | |
253 | #define NRNr14 8 | |
254 | static const _Float128 RNr14[NRNr14 + 1] = | |
255 | { | |
256 | L(-2.446164016404426277577283038988918202456E3), | |
257 | L(6.718753324496563913392217011618096698140E2), | |
258 | L(-4.581631138049836157425391886957389240794E2), | |
259 | L(-2.382844088987092233033215402335026078208E1), | |
260 | L(-7.119237852400600507927038680970936336458E1), | |
261 | L(1.313609646108420136332418282286454287146E1), | |
262 | L(-6.188608702082264389155862490056401365834E0), | |
263 | L(-2.787116601106678287277373011101132659279E-2), | |
264 | L(-2.230395570574153963203348263549700967918E-2) | |
265 | }; | |
266 | #define NRDr14 7 | |
267 | static const _Float128 RDr14[NRDr14 + 1] = | |
268 | { | |
269 | L(2.495187439241869732696223349840963702875E3), | |
270 | L(2.503549449872925580011284635695738412162E2), | |
271 | L(1.159033560988895481698051531263861842461E3), | |
272 | L(9.493751466542304491261487998684383688622E1), | |
273 | L(2.276214929562354328261422263078480321204E2), | |
274 | L(1.367697521219069280358984081407807931847E1), | |
275 | L(2.276988395995528495055594829206582732682E1), | |
276 | L(7.647745753648996559837591812375456641163E-1) | |
277 | /* 1.0E0 */ | |
278 | }; | |
279 | /* erfc(0.375) = C14a + C14b to extra precision. */ | |
280 | static const _Float128 C14a = L(0.5958709716796875); | |
281 | static const _Float128 C14b = L(1.2118885490201676174914080878232469565953E-5); | |
282 | ||
283 | /* erfc(x + 0.5) = erfc(0.5) + x R(x) | |
284 | 0 <= x < 0.125 | |
285 | Peak relative error 4.7e-36 */ | |
286 | #define NRNr15 8 | |
287 | static const _Float128 RNr15[NRNr15 + 1] = | |
288 | { | |
289 | L(-2.624212418011181487924855581955853461925E3), | |
290 | L(8.473828904647825181073831556439301342756E2), | |
291 | L(-5.286207458628380765099405359607331669027E2), | |
292 | L(-3.895781234155315729088407259045269652318E1), | |
293 | L(-6.200857908065163618041240848728398496256E1), | |
294 | L(1.469324610346924001393137895116129204737E1), | |
295 | L(-6.961356525370658572800674953305625578903E0), | |
296 | L(5.145724386641163809595512876629030548495E-3), | |
297 | L(1.990253655948179713415957791776180406812E-2) | |
298 | }; | |
299 | #define NRDr15 7 | |
300 | static const _Float128 RDr15[NRDr15 + 1] = | |
301 | { | |
302 | L(2.986190760847974943034021764693341524962E3), | |
303 | L(5.288262758961073066335410218650047725985E2), | |
304 | L(1.363649178071006978355113026427856008978E3), | |
305 | L(1.921707975649915894241864988942255320833E2), | |
306 | L(2.588651100651029023069013885900085533226E2), | |
307 | L(2.628752920321455606558942309396855629459E1), | |
308 | L(2.455649035885114308978333741080991380610E1), | |
309 | L(1.378826653595128464383127836412100939126E0) | |
310 | /* 1.0E0 */ | |
311 | }; | |
312 | /* erfc(0.5) = C15a + C15b to extra precision. */ | |
313 | static const _Float128 C15a = L(0.4794921875); | |
314 | static const _Float128 C15b = L(7.9346869534623172533461080354712635484242E-6); | |
315 | ||
316 | /* erfc(x + 0.625) = erfc(0.625) + x R(x) | |
317 | 0 <= x < 0.125 | |
318 | Peak relative error 5.1e-36 */ | |
319 | #define NRNr16 8 | |
320 | static const _Float128 RNr16[NRNr16 + 1] = | |
321 | { | |
322 | L(-2.347887943200680563784690094002722906820E3), | |
323 | L(8.008590660692105004780722726421020136482E2), | |
324 | L(-5.257363310384119728760181252132311447963E2), | |
325 | L(-4.471737717857801230450290232600243795637E1), | |
326 | L(-4.849540386452573306708795324759300320304E1), | |
327 | L(1.140885264677134679275986782978655952843E1), | |
328 | L(-6.731591085460269447926746876983786152300E0), | |
329 | L(1.370831653033047440345050025876085121231E-1), | |
330 | L(2.022958279982138755020825717073966576670E-2), | |
331 | }; | |
332 | #define NRDr16 7 | |
333 | static const _Float128 RDr16[NRDr16 + 1] = | |
334 | { | |
335 | L(3.075166170024837215399323264868308087281E3), | |
336 | L(8.730468942160798031608053127270430036627E2), | |
337 | L(1.458472799166340479742581949088453244767E3), | |
338 | L(3.230423687568019709453130785873540386217E2), | |
339 | L(2.804009872719893612081109617983169474655E2), | |
340 | L(4.465334221323222943418085830026979293091E1), | |
341 | L(2.612723259683205928103787842214809134746E1), | |
342 | L(2.341526751185244109722204018543276124997E0), | |
343 | /* 1.0E0 */ | |
344 | }; | |
345 | /* erfc(0.625) = C16a + C16b to extra precision. */ | |
346 | static const _Float128 C16a = L(0.3767547607421875); | |
347 | static const _Float128 C16b = L(4.3570693945275513594941232097252997287766E-6); | |
348 | ||
349 | /* erfc(x + 0.75) = erfc(0.75) + x R(x) | |
350 | 0 <= x < 0.125 | |
351 | Peak relative error 1.7e-35 */ | |
352 | #define NRNr17 8 | |
353 | static const _Float128 RNr17[NRNr17 + 1] = | |
354 | { | |
355 | L(-1.767068734220277728233364375724380366826E3), | |
356 | L(6.693746645665242832426891888805363898707E2), | |
357 | L(-4.746224241837275958126060307406616817753E2), | |
358 | L(-2.274160637728782675145666064841883803196E1), | |
359 | L(-3.541232266140939050094370552538987982637E1), | |
360 | L(6.988950514747052676394491563585179503865E0), | |
361 | L(-5.807687216836540830881352383529281215100E0), | |
362 | L(3.631915988567346438830283503729569443642E-1), | |
363 | L(-1.488945487149634820537348176770282391202E-2) | |
364 | }; | |
365 | #define NRDr17 7 | |
366 | static const _Float128 RDr17[NRDr17 + 1] = | |
367 | { | |
368 | L(2.748457523498150741964464942246913394647E3), | |
369 | L(1.020213390713477686776037331757871252652E3), | |
370 | L(1.388857635935432621972601695296561952738E3), | |
371 | L(3.903363681143817750895999579637315491087E2), | |
372 | L(2.784568344378139499217928969529219886578E2), | |
373 | L(5.555800830216764702779238020065345401144E1), | |
374 | L(2.646215470959050279430447295801291168941E1), | |
375 | L(2.984905282103517497081766758550112011265E0), | |
376 | /* 1.0E0 */ | |
377 | }; | |
378 | /* erfc(0.75) = C17a + C17b to extra precision. */ | |
379 | static const _Float128 C17a = L(0.2888336181640625); | |
380 | static const _Float128 C17b = L(1.0748182422368401062165408589222625794046E-5); | |
381 | ||
382 | ||
383 | /* erfc(x + 0.875) = erfc(0.875) + x R(x) | |
384 | 0 <= x < 0.125 | |
385 | Peak relative error 2.2e-35 */ | |
386 | #define NRNr18 8 | |
387 | static const _Float128 RNr18[NRNr18 + 1] = | |
388 | { | |
389 | L(-1.342044899087593397419622771847219619588E3), | |
390 | L(6.127221294229172997509252330961641850598E2), | |
391 | L(-4.519821356522291185621206350470820610727E2), | |
392 | L(1.223275177825128732497510264197915160235E1), | |
393 | L(-2.730789571382971355625020710543532867692E1), | |
394 | L(4.045181204921538886880171727755445395862E0), | |
395 | L(-4.925146477876592723401384464691452700539E0), | |
396 | L(5.933878036611279244654299924101068088582E-1), | |
397 | L(-5.557645435858916025452563379795159124753E-2) | |
398 | }; | |
399 | #define NRDr18 7 | |
400 | static const _Float128 RDr18[NRDr18 + 1] = | |
401 | { | |
402 | L(2.557518000661700588758505116291983092951E3), | |
403 | L(1.070171433382888994954602511991940418588E3), | |
404 | L(1.344842834423493081054489613250688918709E3), | |
405 | L(4.161144478449381901208660598266288188426E2), | |
406 | L(2.763670252219855198052378138756906980422E2), | |
407 | L(5.998153487868943708236273854747564557632E1), | |
408 | L(2.657695108438628847733050476209037025318E1), | |
409 | L(3.252140524394421868923289114410336976512E0), | |
410 | /* 1.0E0 */ | |
411 | }; | |
412 | /* erfc(0.875) = C18a + C18b to extra precision. */ | |
413 | static const _Float128 C18a = L(0.215911865234375); | |
414 | static const _Float128 C18b = L(1.3073705765341685464282101150637224028267E-5); | |
415 | ||
416 | /* erfc(x + 1.0) = erfc(1.0) + x R(x) | |
417 | 0 <= x < 0.125 | |
418 | Peak relative error 1.6e-35 */ | |
419 | #define NRNr19 8 | |
420 | static const _Float128 RNr19[NRNr19 + 1] = | |
421 | { | |
422 | L(-1.139180936454157193495882956565663294826E3), | |
423 | L(6.134903129086899737514712477207945973616E2), | |
424 | L(-4.628909024715329562325555164720732868263E2), | |
425 | L(4.165702387210732352564932347500364010833E1), | |
426 | L(-2.286979913515229747204101330405771801610E1), | |
427 | L(1.870695256449872743066783202326943667722E0), | |
428 | L(-4.177486601273105752879868187237000032364E0), | |
429 | L(7.533980372789646140112424811291782526263E-1), | |
430 | L(-8.629945436917752003058064731308767664446E-2) | |
431 | }; | |
432 | #define NRDr19 7 | |
433 | static const _Float128 RDr19[NRDr19 + 1] = | |
434 | { | |
435 | L(2.744303447981132701432716278363418643778E3), | |
436 | L(1.266396359526187065222528050591302171471E3), | |
437 | L(1.466739461422073351497972255511919814273E3), | |
438 | L(4.868710570759693955597496520298058147162E2), | |
439 | L(2.993694301559756046478189634131722579643E2), | |
440 | L(6.868976819510254139741559102693828237440E1), | |
441 | L(2.801505816247677193480190483913753613630E1), | |
442 | L(3.604439909194350263552750347742663954481E0), | |
443 | /* 1.0E0 */ | |
444 | }; | |
445 | /* erfc(1.0) = C19a + C19b to extra precision. */ | |
446 | static const _Float128 C19a = L(0.15728759765625); | |
447 | static const _Float128 C19b = L(1.1609394035130658779364917390740703933002E-5); | |
448 | ||
449 | /* erfc(x + 1.125) = erfc(1.125) + x R(x) | |
450 | 0 <= x < 0.125 | |
451 | Peak relative error 3.6e-36 */ | |
452 | #define NRNr20 8 | |
453 | static const _Float128 RNr20[NRNr20 + 1] = | |
454 | { | |
455 | L(-9.652706916457973956366721379612508047640E2), | |
456 | L(5.577066396050932776683469951773643880634E2), | |
457 | L(-4.406335508848496713572223098693575485978E2), | |
458 | L(5.202893466490242733570232680736966655434E1), | |
459 | L(-1.931311847665757913322495948705563937159E1), | |
460 | L(-9.364318268748287664267341457164918090611E-2), | |
461 | L(-3.306390351286352764891355375882586201069E0), | |
462 | L(7.573806045289044647727613003096916516475E-1), | |
463 | L(-9.611744011489092894027478899545635991213E-2) | |
464 | }; | |
465 | #define NRDr20 7 | |
466 | static const _Float128 RDr20[NRDr20 + 1] = | |
467 | { | |
468 | L(3.032829629520142564106649167182428189014E3), | |
469 | L(1.659648470721967719961167083684972196891E3), | |
470 | L(1.703545128657284619402511356932569292535E3), | |
471 | L(6.393465677731598872500200253155257708763E2), | |
472 | L(3.489131397281030947405287112726059221934E2), | |
473 | L(8.848641738570783406484348434387611713070E1), | |
474 | L(3.132269062552392974833215844236160958502E1), | |
475 | L(4.430131663290563523933419966185230513168E0) | |
476 | /* 1.0E0 */ | |
477 | }; | |
478 | /* erfc(1.125) = C20a + C20b to extra precision. */ | |
479 | static const _Float128 C20a = L(0.111602783203125); | |
480 | static const _Float128 C20b = L(8.9850951672359304215530728365232161564636E-6); | |
481 | ||
482 | /* erfc(1/x) = 1/x exp (-1/x^2 - 0.5625 + R(1/x^2)) | |
483 | 7/8 <= 1/x < 1 | |
484 | Peak relative error 1.4e-35 */ | |
485 | #define NRNr8 9 | |
486 | static const _Float128 RNr8[NRNr8 + 1] = | |
487 | { | |
488 | L(3.587451489255356250759834295199296936784E1), | |
489 | L(5.406249749087340431871378009874875889602E2), | |
490 | L(2.931301290625250886238822286506381194157E3), | |
491 | L(7.359254185241795584113047248898753470923E3), | |
492 | L(9.201031849810636104112101947312492532314E3), | |
493 | L(5.749697096193191467751650366613289284777E3), | |
494 | L(1.710415234419860825710780802678697889231E3), | |
495 | L(2.150753982543378580859546706243022719599E2), | |
496 | L(8.740953582272147335100537849981160931197E0), | |
497 | L(4.876422978828717219629814794707963640913E-2) | |
498 | }; | |
499 | #define NRDr8 8 | |
500 | static const _Float128 RDr8[NRDr8 + 1] = | |
501 | { | |
502 | L(6.358593134096908350929496535931630140282E1), | |
503 | L(9.900253816552450073757174323424051765523E2), | |
504 | L(5.642928777856801020545245437089490805186E3), | |
505 | L(1.524195375199570868195152698617273739609E4), | |
506 | L(2.113829644500006749947332935305800887345E4), | |
507 | L(1.526438562626465706267943737310282977138E4), | |
508 | L(5.561370922149241457131421914140039411782E3), | |
509 | L(9.394035530179705051609070428036834496942E2), | |
510 | L(6.147019596150394577984175188032707343615E1) | |
511 | /* 1.0E0 */ | |
512 | }; | |
513 | ||
514 | /* erfc(1/x) = 1/x exp (-1/x^2 - 0.5625 + R(1/x^2)) | |
515 | 0.75 <= 1/x <= 0.875 | |
516 | Peak relative error 2.0e-36 */ | |
517 | #define NRNr7 9 | |
518 | static const _Float128 RNr7[NRNr7 + 1] = | |
519 | { | |
520 | L(1.686222193385987690785945787708644476545E1), | |
521 | L(1.178224543567604215602418571310612066594E3), | |
522 | L(1.764550584290149466653899886088166091093E4), | |
523 | L(1.073758321890334822002849369898232811561E5), | |
524 | L(3.132840749205943137619839114451290324371E5), | |
525 | L(4.607864939974100224615527007793867585915E5), | |
526 | L(3.389781820105852303125270837910972384510E5), | |
527 | L(1.174042187110565202875011358512564753399E5), | |
528 | L(1.660013606011167144046604892622504338313E4), | |
529 | L(6.700393957480661937695573729183733234400E2) | |
530 | }; | |
531 | #define NRDr7 9 | |
532 | static const _Float128 RDr7[NRDr7 + 1] = | |
533 | { | |
534 | L(-1.709305024718358874701575813642933561169E3), | |
535 | L(-3.280033887481333199580464617020514788369E4), | |
536 | L(-2.345284228022521885093072363418750835214E5), | |
537 | L(-8.086758123097763971926711729242327554917E5), | |
538 | L(-1.456900414510108718402423999575992450138E6), | |
539 | L(-1.391654264881255068392389037292702041855E6), | |
540 | L(-6.842360801869939983674527468509852583855E5), | |
541 | L(-1.597430214446573566179675395199807533371E5), | |
542 | L(-1.488876130609876681421645314851760773480E4), | |
543 | L(-3.511762950935060301403599443436465645703E2) | |
544 | /* 1.0E0 */ | |
545 | }; | |
546 | ||
547 | /* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2)) | |
548 | 5/8 <= 1/x < 3/4 | |
549 | Peak relative error 1.9e-35 */ | |
550 | #define NRNr6 9 | |
551 | static const _Float128 RNr6[NRNr6 + 1] = | |
552 | { | |
553 | L(1.642076876176834390623842732352935761108E0), | |
554 | L(1.207150003611117689000664385596211076662E2), | |
555 | L(2.119260779316389904742873816462800103939E3), | |
556 | L(1.562942227734663441801452930916044224174E4), | |
557 | L(5.656779189549710079988084081145693580479E4), | |
558 | L(1.052166241021481691922831746350942786299E5), | |
559 | L(9.949798524786000595621602790068349165758E4), | |
560 | L(4.491790734080265043407035220188849562856E4), | |
561 | L(8.377074098301530326270432059434791287601E3), | |
562 | L(4.506934806567986810091824791963991057083E2) | |
563 | }; | |
564 | #define NRDr6 9 | |
565 | static const _Float128 RDr6[NRDr6 + 1] = | |
566 | { | |
567 | L(-1.664557643928263091879301304019826629067E2), | |
568 | L(-3.800035902507656624590531122291160668452E3), | |
569 | L(-3.277028191591734928360050685359277076056E4), | |
570 | L(-1.381359471502885446400589109566587443987E5), | |
571 | L(-3.082204287382581873532528989283748656546E5), | |
572 | L(-3.691071488256738343008271448234631037095E5), | |
573 | L(-2.300482443038349815750714219117566715043E5), | |
574 | L(-6.873955300927636236692803579555752171530E4), | |
575 | L(-8.262158817978334142081581542749986845399E3), | |
576 | L(-2.517122254384430859629423488157361983661E2) | |
577 | /* 1.00 */ | |
578 | }; | |
579 | ||
580 | /* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2)) | |
581 | 1/2 <= 1/x < 5/8 | |
582 | Peak relative error 4.6e-36 */ | |
583 | #define NRNr5 10 | |
584 | static const _Float128 RNr5[NRNr5 + 1] = | |
585 | { | |
586 | L(-3.332258927455285458355550878136506961608E-3), | |
587 | L(-2.697100758900280402659586595884478660721E-1), | |
588 | L(-6.083328551139621521416618424949137195536E0), | |
589 | L(-6.119863528983308012970821226810162441263E1), | |
590 | L(-3.176535282475593173248810678636522589861E2), | |
591 | L(-8.933395175080560925809992467187963260693E2), | |
592 | L(-1.360019508488475978060917477620199499560E3), | |
593 | L(-1.075075579828188621541398761300910213280E3), | |
594 | L(-4.017346561586014822824459436695197089916E2), | |
595 | L(-5.857581368145266249509589726077645791341E1), | |
596 | L(-2.077715925587834606379119585995758954399E0) | |
597 | }; | |
598 | #define NRDr5 9 | |
599 | static const _Float128 RDr5[NRDr5 + 1] = | |
600 | { | |
601 | L(3.377879570417399341550710467744693125385E-1), | |
602 | L(1.021963322742390735430008860602594456187E1), | |
603 | L(1.200847646592942095192766255154827011939E2), | |
604 | L(7.118915528142927104078182863387116942836E2), | |
605 | L(2.318159380062066469386544552429625026238E3), | |
606 | L(4.238729853534009221025582008928765281620E3), | |
607 | L(4.279114907284825886266493994833515580782E3), | |
608 | L(2.257277186663261531053293222591851737504E3), | |
609 | L(5.570475501285054293371908382916063822957E2), | |
610 | L(5.142189243856288981145786492585432443560E1) | |
611 | /* 1.0E0 */ | |
612 | }; | |
613 | ||
614 | /* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2)) | |
615 | 3/8 <= 1/x < 1/2 | |
616 | Peak relative error 2.0e-36 */ | |
617 | #define NRNr4 10 | |
618 | static const _Float128 RNr4[NRNr4 + 1] = | |
619 | { | |
620 | L(3.258530712024527835089319075288494524465E-3), | |
621 | L(2.987056016877277929720231688689431056567E-1), | |
622 | L(8.738729089340199750734409156830371528862E0), | |
623 | L(1.207211160148647782396337792426311125923E2), | |
624 | L(8.997558632489032902250523945248208224445E2), | |
625 | L(3.798025197699757225978410230530640879762E3), | |
626 | L(9.113203668683080975637043118209210146846E3), | |
627 | L(1.203285891339933238608683715194034900149E4), | |
628 | L(8.100647057919140328536743641735339740855E3), | |
629 | L(2.383888249907144945837976899822927411769E3), | |
630 | L(2.127493573166454249221983582495245662319E2) | |
631 | }; | |
632 | #define NRDr4 10 | |
633 | static const _Float128 RDr4[NRDr4 + 1] = | |
634 | { | |
635 | L(-3.303141981514540274165450687270180479586E-1), | |
636 | L(-1.353768629363605300707949368917687066724E1), | |
637 | L(-2.206127630303621521950193783894598987033E2), | |
638 | L(-1.861800338758066696514480386180875607204E3), | |
639 | L(-8.889048775872605708249140016201753255599E3), | |
640 | L(-2.465888106627948210478692168261494857089E4), | |
641 | L(-3.934642211710774494879042116768390014289E4), | |
642 | L(-3.455077258242252974937480623730228841003E4), | |
643 | L(-1.524083977439690284820586063729912653196E4), | |
644 | L(-2.810541887397984804237552337349093953857E3), | |
645 | L(-1.343929553541159933824901621702567066156E2) | |
646 | /* 1.0E0 */ | |
647 | }; | |
648 | ||
649 | /* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2)) | |
650 | 1/4 <= 1/x < 3/8 | |
651 | Peak relative error 8.4e-37 */ | |
652 | #define NRNr3 11 | |
653 | static const _Float128 RNr3[NRNr3 + 1] = | |
654 | { | |
655 | L(-1.952401126551202208698629992497306292987E-6), | |
656 | L(-2.130881743066372952515162564941682716125E-4), | |
657 | L(-8.376493958090190943737529486107282224387E-3), | |
658 | L(-1.650592646560987700661598877522831234791E-1), | |
659 | L(-1.839290818933317338111364667708678163199E0), | |
660 | L(-1.216278715570882422410442318517814388470E1), | |
661 | L(-4.818759344462360427612133632533779091386E1), | |
662 | L(-1.120994661297476876804405329172164436784E2), | |
663 | L(-1.452850765662319264191141091859300126931E2), | |
664 | L(-9.485207851128957108648038238656777241333E1), | |
665 | L(-2.563663855025796641216191848818620020073E1), | |
666 | L(-1.787995944187565676837847610706317833247E0) | |
667 | }; | |
668 | #define NRDr3 10 | |
669 | static const _Float128 RDr3[NRDr3 + 1] = | |
670 | { | |
671 | L(1.979130686770349481460559711878399476903E-4), | |
672 | L(1.156941716128488266238105813374635099057E-2), | |
673 | L(2.752657634309886336431266395637285974292E-1), | |
674 | L(3.482245457248318787349778336603569327521E0), | |
675 | L(2.569347069372696358578399521203959253162E1), | |
676 | L(1.142279000180457419740314694631879921561E2), | |
677 | L(3.056503977190564294341422623108332700840E2), | |
678 | L(4.780844020923794821656358157128719184422E2), | |
679 | L(4.105972727212554277496256802312730410518E2), | |
680 | L(1.724072188063746970865027817017067646246E2), | |
681 | L(2.815939183464818198705278118326590370435E1) | |
682 | /* 1.0E0 */ | |
683 | }; | |
684 | ||
685 | /* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2)) | |
686 | 1/8 <= 1/x < 1/4 | |
687 | Peak relative error 1.5e-36 */ | |
688 | #define NRNr2 11 | |
689 | static const _Float128 RNr2[NRNr2 + 1] = | |
690 | { | |
691 | L(-2.638914383420287212401687401284326363787E-8), | |
692 | L(-3.479198370260633977258201271399116766619E-6), | |
693 | L(-1.783985295335697686382487087502222519983E-4), | |
694 | L(-4.777876933122576014266349277217559356276E-3), | |
695 | L(-7.450634738987325004070761301045014986520E-2), | |
696 | L(-7.068318854874733315971973707247467326619E-1), | |
697 | L(-4.113919921935944795764071670806867038732E0), | |
698 | L(-1.440447573226906222417767283691888875082E1), | |
699 | L(-2.883484031530718428417168042141288943905E1), | |
700 | L(-2.990886974328476387277797361464279931446E1), | |
701 | L(-1.325283914915104866248279787536128997331E1), | |
702 | L(-1.572436106228070195510230310658206154374E0) | |
703 | }; | |
704 | #define NRDr2 10 | |
705 | static const _Float128 RDr2[NRDr2 + 1] = | |
706 | { | |
707 | L(2.675042728136731923554119302571867799673E-6), | |
708 | L(2.170997868451812708585443282998329996268E-4), | |
709 | L(7.249969752687540289422684951196241427445E-3), | |
710 | L(1.302040375859768674620410563307838448508E-1), | |
711 | L(1.380202483082910888897654537144485285549E0), | |
712 | L(8.926594113174165352623847870299170069350E0), | |
713 | L(3.521089584782616472372909095331572607185E1), | |
714 | L(8.233547427533181375185259050330809105570E1), | |
715 | L(1.072971579885803033079469639073292840135E2), | |
716 | L(6.943803113337964469736022094105143158033E1), | |
717 | L(1.775695341031607738233608307835017282662E1) | |
718 | /* 1.0E0 */ | |
719 | }; | |
720 | ||
721 | /* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2)) | |
722 | 1/128 <= 1/x < 1/8 | |
723 | Peak relative error 2.2e-36 */ | |
724 | #define NRNr1 9 | |
725 | static const _Float128 RNr1[NRNr1 + 1] = | |
726 | { | |
727 | L(-4.250780883202361946697751475473042685782E-8), | |
728 | L(-5.375777053288612282487696975623206383019E-6), | |
729 | L(-2.573645949220896816208565944117382460452E-4), | |
730 | L(-6.199032928113542080263152610799113086319E-3), | |
731 | L(-8.262721198693404060380104048479916247786E-2), | |
732 | L(-6.242615227257324746371284637695778043982E-1), | |
733 | L(-2.609874739199595400225113299437099626386E0), | |
734 | L(-5.581967563336676737146358534602770006970E0), | |
735 | L(-5.124398923356022609707490956634280573882E0), | |
736 | L(-1.290865243944292370661544030414667556649E0) | |
737 | }; | |
738 | #define NRDr1 8 | |
739 | static const _Float128 RDr1[NRDr1 + 1] = | |
740 | { | |
741 | L(4.308976661749509034845251315983612976224E-6), | |
742 | L(3.265390126432780184125233455960049294580E-4), | |
743 | L(9.811328839187040701901866531796570418691E-3), | |
744 | L(1.511222515036021033410078631914783519649E-1), | |
745 | L(1.289264341917429958858379585970225092274E0), | |
746 | L(6.147640356182230769548007536914983522270E0), | |
747 | L(1.573966871337739784518246317003956180750E1), | |
748 | L(1.955534123435095067199574045529218238263E1), | |
749 | L(9.472613121363135472247929109615785855865E0) | |
750 | /* 1.0E0 */ | |
751 | }; | |
752 | ||
753 | ||
754 | _Float128 | |
755 | __erfl (_Float128 x) | |
756 | { | |
757 | _Float128 a, y, z; | |
758 | int32_t i, ix, sign; | |
759 | ieee854_long_double_shape_type u; | |
760 | ||
761 | u.value = x; | |
762 | sign = u.parts32.w0; | |
763 | ix = sign & 0x7fffffff; | |
764 | ||
765 | if (ix >= 0x7fff0000) | |
766 | { /* erf(nan)=nan */ | |
767 | i = ((sign & 0xffff0000) >> 31) << 1; | |
768 | return (_Float128) (1 - i) + one / x; /* erf(+-inf)=+-1 */ | |
769 | } | |
770 | ||
771 | if (ix >= 0x3fff0000) /* |x| >= 1.0 */ | |
772 | { | |
773 | if (ix >= 0x40030000 && sign > 0) | |
774 | return one; /* x >= 16, avoid spurious underflow from erfc. */ | |
775 | y = __erfcl (x); | |
776 | return (one - y); | |
777 | /* return (one - __erfcl (x)); */ | |
778 | } | |
779 | u.parts32.w0 = ix; | |
780 | a = u.value; | |
781 | z = x * x; | |
782 | if (ix < 0x3ffec000) /* a < 0.875 */ | |
783 | { | |
784 | if (ix < 0x3fc60000) /* |x|<2**-57 */ | |
785 | { | |
786 | if (ix < 0x00080000) | |
787 | { | |
788 | /* Avoid spurious underflow. */ | |
789 | _Float128 ret = 0.0625 * (16.0 * x + (16.0 * efx) * x); | |
790 | math_check_force_underflow (ret); | |
791 | return ret; | |
792 | } | |
793 | return x + efx * x; | |
794 | } | |
795 | y = a + a * neval (z, TN1, NTN1) / deval (z, TD1, NTD1); | |
796 | } | |
797 | else | |
798 | { | |
799 | a = a - one; | |
800 | y = erf_const + neval (a, TN2, NTN2) / deval (a, TD2, NTD2); | |
801 | } | |
802 | ||
803 | if (sign & 0x80000000) /* x < 0 */ | |
804 | y = -y; | |
805 | return( y ); | |
806 | } | |
807 | ||
808 | libm_alias_ldouble (__erf, erf) | |
809 | _Float128 | |
810 | __erfcl (_Float128 x) | |
811 | { | |
812 | _Float128 y, z, p, r; | |
813 | int32_t i, ix, sign; | |
814 | ieee854_long_double_shape_type u; | |
815 | ||
816 | u.value = x; | |
817 | sign = u.parts32.w0; | |
818 | ix = sign & 0x7fffffff; | |
819 | u.parts32.w0 = ix; | |
820 | ||
821 | if (ix >= 0x7fff0000) | |
822 | { /* erfc(nan)=nan */ | |
823 | /* erfc(+-inf)=0,2 */ | |
824 | return (_Float128) (((uint32_t) sign >> 31) << 1) + one / x; | |
825 | } | |
826 | ||
827 | if (ix < 0x3ffd0000) /* |x| <1/4 */ | |
828 | { | |
829 | if (ix < 0x3f8d0000) /* |x|<2**-114 */ | |
830 | return one - x; | |
831 | return one - __erfl (x); | |
832 | } | |
833 | if (ix < 0x3fff4000) /* 1.25 */ | |
834 | { | |
835 | x = u.value; | |
836 | i = 8.0 * x; | |
837 | switch (i) | |
838 | { | |
839 | case 2: | |
840 | z = x - L(0.25); | |
841 | y = C13b + z * neval (z, RNr13, NRNr13) / deval (z, RDr13, NRDr13); | |
842 | y += C13a; | |
843 | break; | |
844 | case 3: | |
845 | z = x - L(0.375); | |
846 | y = C14b + z * neval (z, RNr14, NRNr14) / deval (z, RDr14, NRDr14); | |
847 | y += C14a; | |
848 | break; | |
849 | case 4: | |
850 | z = x - L(0.5); | |
851 | y = C15b + z * neval (z, RNr15, NRNr15) / deval (z, RDr15, NRDr15); | |
852 | y += C15a; | |
853 | break; | |
854 | case 5: | |
855 | z = x - L(0.625); | |
856 | y = C16b + z * neval (z, RNr16, NRNr16) / deval (z, RDr16, NRDr16); | |
857 | y += C16a; | |
858 | break; | |
859 | case 6: | |
860 | z = x - L(0.75); | |
861 | y = C17b + z * neval (z, RNr17, NRNr17) / deval (z, RDr17, NRDr17); | |
862 | y += C17a; | |
863 | break; | |
864 | case 7: | |
865 | z = x - L(0.875); | |
866 | y = C18b + z * neval (z, RNr18, NRNr18) / deval (z, RDr18, NRDr18); | |
867 | y += C18a; | |
868 | break; | |
869 | case 8: | |
870 | z = x - 1; | |
871 | y = C19b + z * neval (z, RNr19, NRNr19) / deval (z, RDr19, NRDr19); | |
872 | y += C19a; | |
873 | break; | |
874 | default: /* i == 9. */ | |
875 | z = x - L(1.125); | |
876 | y = C20b + z * neval (z, RNr20, NRNr20) / deval (z, RDr20, NRDr20); | |
877 | y += C20a; | |
878 | break; | |
879 | } | |
880 | if (sign & 0x80000000) | |
881 | y = 2 - y; | |
882 | return y; | |
883 | } | |
884 | /* 1.25 < |x| < 107 */ | |
885 | if (ix < 0x4005ac00) | |
886 | { | |
887 | /* x < -9 */ | |
888 | if ((ix >= 0x40022000) && (sign & 0x80000000)) | |
889 | return two - tiny; | |
890 | ||
891 | x = fabsl (x); | |
892 | z = one / (x * x); | |
893 | i = 8.0 / x; | |
894 | switch (i) | |
895 | { | |
896 | default: | |
897 | case 0: | |
898 | p = neval (z, RNr1, NRNr1) / deval (z, RDr1, NRDr1); | |
899 | break; | |
900 | case 1: | |
901 | p = neval (z, RNr2, NRNr2) / deval (z, RDr2, NRDr2); | |
902 | break; | |
903 | case 2: | |
904 | p = neval (z, RNr3, NRNr3) / deval (z, RDr3, NRDr3); | |
905 | break; | |
906 | case 3: | |
907 | p = neval (z, RNr4, NRNr4) / deval (z, RDr4, NRDr4); | |
908 | break; | |
909 | case 4: | |
910 | p = neval (z, RNr5, NRNr5) / deval (z, RDr5, NRDr5); | |
911 | break; | |
912 | case 5: | |
913 | p = neval (z, RNr6, NRNr6) / deval (z, RDr6, NRDr6); | |
914 | break; | |
915 | case 6: | |
916 | p = neval (z, RNr7, NRNr7) / deval (z, RDr7, NRDr7); | |
917 | break; | |
918 | case 7: | |
919 | p = neval (z, RNr8, NRNr8) / deval (z, RDr8, NRDr8); | |
920 | break; | |
921 | } | |
922 | u.value = x; | |
923 | u.parts32.w3 = 0; | |
924 | u.parts32.w2 &= 0xfe000000; | |
925 | z = u.value; | |
926 | r = __ieee754_expl (-z * z - 0.5625) * | |
927 | __ieee754_expl ((z - x) * (z + x) + p); | |
928 | if ((sign & 0x80000000) == 0) | |
929 | { | |
930 | _Float128 ret = r / x; | |
931 | if (ret == 0) | |
932 | __set_errno (ERANGE); | |
933 | return ret; | |
934 | } | |
935 | else | |
936 | return two - r / x; | |
937 | } | |
938 | else | |
939 | { | |
940 | if ((sign & 0x80000000) == 0) | |
941 | { | |
942 | __set_errno (ERANGE); | |
943 | return tiny * tiny; | |
944 | } | |
945 | else | |
946 | return two - tiny; | |
947 | } | |
948 | } | |
949 | ||
950 | libm_alias_ldouble (__erfc, erfc) |