31  *     erf(x)  =  ---------  | exp(-t*t)dt
35  *     erfc(x) =  1-erf(x)
37  *              erf(-x) = -erf(x)
38  *              erfc(-x) = 2 - erfc(x)
43  *          erfc(x) = 1 - erf(x)           if x in [-.84375,0.25]
44  *                  = 0.5 + ((0.5-x)-x*R)  if x in [0.25,0.84375]
46  *          erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....)
54  *      2. For |x| in [0.84375,1.25], let s = |x| - 1, and
57  *      erfc(x) = (1-c)  - P1(s)/Q1(s) if x > 0
65  *      erfc(x) = (1/x)*exp(-x*x-0.5625+R1(z)/S1(z))
67  *      erf(x)  = 1 - erfc(x)
70  *      erfc(x) = (1/x)*exp(-x*x-0.5625+R2/S2) if x > 0
71  *                      = 2.0 - (1/x)*exp(-x*x-0.5625+R2(z)/S2(z))
72  *                             if -6.666<x<0
73  *                      = 2.0 - tiny            (if x <= -6.666)
75  *      erf(x)  = sign(x)*(1.0 - erfc(x)) if x < 6.666, else
76  *      erf(x)  = sign(x)*(1.0 - tiny)
78  *         To compute exp(-x*x-0.5625+R/S), let s be a single
80  *              -x*x = -s*s + (s-x)*(s+x)
81  *              exp(-x*x-0.5626+R/S) =
82  *                      exp(-s*s-0.5625)*exp((s-x)*(s+x)+R/S);
85  *                        exp(-x*x)
86  *              erfc(x) ~ ---------- * ( 1 + Poly(1/x^2) )
90  *      erf(x)  = sign(x) *(1 - tiny)  (raise inexact)
92  *                      = 2 - tiny if x<0
95  *      erf(0)  = 0, erf(inf)  = 1, erf(-inf) = -1,
96  *      erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2,
119 /* 8 * (2/sqrt(pi) - 1) */
123 	-2.808533301997696164408397079650699163276E6L,
124 	-3.314325479115357458197119660818768924100E5L,
125 	-6.848684465326256109712135497895525446398E4L,
126 	-2.657817695110739185591505062971929859314E3L,
127 	-1.655310302737837556654146291646499062882E2L,
129 qq[6] = {  variable
143    -0.15625 <= x <= +.25
144    Peak relative error 8.5e-22  */
146 	-1.076952146179812072156734957705102256059E0L,
148 	-5.339153975012804282890066622962070115606E1L,
151 	-2.360236618396952560064259585299045804293E0L,
153 	 9.394994446747752308256773044667843200719E-2L,
169 /* erfc(1/x) = x exp (-1/x^2 - 0.5625 + ra(x^2)/sa(x^2))
171    Peak relative error 3.1e-21  */
173 	1.363566591833846324191000679620738857234E-1L,
184 	-1.382234625202480685182526402169222331847E1L,
185 	-3.315638835627950255832519203687435946482E2L,
186 	-2.949124863912936259747237164260785326692E3L,
187 	-1.246622099070875940506391433635999693661E4L,
188 	-2.673079795851665428695842853070996219632E4L,
189 	-2.880269786660559337358397106518918220991E4L,
190 	-1.450600228493968044773354186390390823713E4L,
191 	-2.874539731125893533960680525192064277816E3L,
192 	-1.402241261419067750237395034116942296027E2L,
199 /* erfc(1/x) = x exp (-1/x^2 - 0.5625 + rb(x^2)/sb(x^2))
201    Peak relative error 4.2e-22  */
203 	-4.869587348270494309550558460786501252369E-5L,
204 	-4.030199390527997378549161722412466959403E-3L,
205 	-9.434425866377037610206443566288917589122E-2L,
206 	-9.319032754357658601200655161585539404155E-1L,
207 	-4.273788174307459947350256581445442062291E0L,
208 	-8.842289940696150508373541814064198259278E0L,
209 	-7.069215249419887403187988144752613025255E0L,
210 	-1.401228723639514787920274427443330704764E0L,
213 	4.936254964107175160157544545879293019085E-3L,
214 	1.583457624037795744377163924895349412015E-1L,
222 /* erfc(1/x) = x exp (-1/x^2 - 0.5625 + rc(x^2)/sc(x^2))
224    Peak relative error 1.1e-21  */
226 	-8.299617545269701963973537248996670806850E-5L,
227 	-6.243845685115818513578933902532056244108E-3L,
228 	-1.141667210620380223113693474478394397230E-1L,
229 	-7.521343797212024245375240432734425789409E-1L,
230 	-1.765321928311155824664963633786967602934E0L,
231 	-1.029403473103215800456761180695263439188E0L,
234 	8.413244363014929493035952542677768808601E-3L,
235 	2.065114333816877479753334599639158060979E-1L,
246 	s = fabsl(x) - 1;  in erfc1()
251 	return 1 - erx - P / Q;  in erfc1()
281 	u.i.m &= -1ULL << 40;  in erfc2()
283 	return expl(-z*z - 0.5625) * expl((z - x) * (z + x) + R / S) / x;  in erfc2()
294 		/* erf(nan)=nan, erf(+-inf)=+-1 */  in erfl()
295 		return 1 - 2*sign + 1/x;  in erfl()
297 		if (ix < 0x3fde8000) {  /* |x| < 2**-33 */  in erfl()
303 		s = qq[0] + z * (qq[1] +  in erfl()
304 		     z * (qq[2] + z * (qq[3] + z * (qq[4] + z * (qq[5] + z)))));  in erfl()
309 		y = 1 - erfc2(ix,x);  in erfl()
311 		y = 1 - 0x1p-16382L;  in erfl()
312 	return sign ? -y : y;  in erfl()
323 		/* erfc(nan) = nan, erfc(+-inf) = 0,2 */  in erfcl()
326 		if (ix < 0x3fbe0000)  /* |x| < 2**-65 */  in erfcl()
327 			return 1.0 - x;  in erfcl()
331 		s = qq[0] + z * (qq[1] +  in erfcl()
332 		     z * (qq[2] + z * (qq[3] + z * (qq[4] + z * (qq[5] + z)))));  in erfcl()
335 			return 1.0 - (x + x * y);  in erfcl()
336 		return 0.5 - (x - 0.5 + x * y);  in erfcl()
339 		return sign ? 2 - erfc2(ix,x) : erfc2(ix,x);  in erfcl()
340 	y = 0x1p-16382L;  in erfcl()
341 	return sign ? 2 - y : y*y;  in erfcl()