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