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)