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