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