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