xref: /illumos-gate/usr/src/lib/libm/common/LD/erfl.c (revision ddc0e0b5)
125c28e83SPiotr Jasiukajtis /*
225c28e83SPiotr Jasiukajtis  * CDDL HEADER START
325c28e83SPiotr Jasiukajtis  *
425c28e83SPiotr Jasiukajtis  * The contents of this file are subject to the terms of the
525c28e83SPiotr Jasiukajtis  * Common Development and Distribution License (the "License").
625c28e83SPiotr Jasiukajtis  * You may not use this file except in compliance with the License.
725c28e83SPiotr Jasiukajtis  *
825c28e83SPiotr Jasiukajtis  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
925c28e83SPiotr Jasiukajtis  * or http://www.opensolaris.org/os/licensing.
1025c28e83SPiotr Jasiukajtis  * See the License for the specific language governing permissions
1125c28e83SPiotr Jasiukajtis  * and limitations under the License.
1225c28e83SPiotr Jasiukajtis  *
1325c28e83SPiotr Jasiukajtis  * When distributing Covered Code, include this CDDL HEADER in each
1425c28e83SPiotr Jasiukajtis  * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
1525c28e83SPiotr Jasiukajtis  * If applicable, add the following below this CDDL HEADER, with the
1625c28e83SPiotr Jasiukajtis  * fields enclosed by brackets "[]" replaced with your own identifying
1725c28e83SPiotr Jasiukajtis  * information: Portions Copyright [yyyy] [name of copyright owner]
1825c28e83SPiotr Jasiukajtis  *
1925c28e83SPiotr Jasiukajtis  * CDDL HEADER END
2025c28e83SPiotr Jasiukajtis  */
2125c28e83SPiotr Jasiukajtis 
2225c28e83SPiotr Jasiukajtis /*
2325c28e83SPiotr Jasiukajtis  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
2425c28e83SPiotr Jasiukajtis  */
2525c28e83SPiotr Jasiukajtis /*
2625c28e83SPiotr Jasiukajtis  * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
2725c28e83SPiotr Jasiukajtis  * Use is subject to license terms.
2825c28e83SPiotr Jasiukajtis  */
2925c28e83SPiotr Jasiukajtis 
3025c28e83SPiotr Jasiukajtis /* long double function erf,erfc (long double x)
3125c28e83SPiotr Jasiukajtis  * K.C. Ng, September, 1989.
3225c28e83SPiotr Jasiukajtis  *			     x
3325c28e83SPiotr Jasiukajtis  *		      2      |\
3425c28e83SPiotr Jasiukajtis  *     erf(x)  =  ---------  | exp(-t*t)dt
35*ddc0e0b5SRichard Lowe  *	 	   sqrt(pi) \|
3625c28e83SPiotr Jasiukajtis  *			     0
3725c28e83SPiotr Jasiukajtis  *
3825c28e83SPiotr Jasiukajtis  *     erfc(x) =  1-erf(x)
3925c28e83SPiotr Jasiukajtis  *
4025c28e83SPiotr Jasiukajtis  * method:
4125c28e83SPiotr Jasiukajtis  * 	Since erf(-x) = -erf(x), we assume x>=0.
4225c28e83SPiotr Jasiukajtis  *	For x near 0, we have the expansion
43*ddc0e0b5SRichard Lowe  *
4425c28e83SPiotr Jasiukajtis  *     	    erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....).
4525c28e83SPiotr Jasiukajtis  *
4625c28e83SPiotr Jasiukajtis  * 	Since 2/sqrt(pi) = 1.128379167095512573896158903121545171688,
4725c28e83SPiotr Jasiukajtis  *	we use x + x*P(x^2) to approximate erf(x). This formula will
4825c28e83SPiotr Jasiukajtis  *	guarantee the error less than one ulp where x is not too far
4925c28e83SPiotr Jasiukajtis  *	away from 0. We note that erf(x)=x at x = 0.6174...... After
50*ddc0e0b5SRichard Lowe  *	some experiment, we choose the following approximation on
5125c28e83SPiotr Jasiukajtis  *	interval [0,0.84375].
5225c28e83SPiotr Jasiukajtis  *
5325c28e83SPiotr Jasiukajtis  *	For x in [0,0.84375]
5425c28e83SPiotr Jasiukajtis  *		   2		    2        4		     40
5525c28e83SPiotr Jasiukajtis  *	   P = 	P(x ) = (p0 + p1 * x + p2 * x + ... + p20 * x  )
5625c28e83SPiotr Jasiukajtis  *
5725c28e83SPiotr Jasiukajtis  *	   erf(x)  = x + x*P
5825c28e83SPiotr Jasiukajtis  *	   erfc(x) = 1 - erf(x) 	  if x<=0.25
5925c28e83SPiotr Jasiukajtis  *		   = 0.5 + ((0.5-x)-x*P)  if x in [0.25,0.84375]
6025c28e83SPiotr Jasiukajtis  *	precision: |P(x^2)-(erf(x)-x)/x| <= 2**-122.50
6125c28e83SPiotr Jasiukajtis  *
62*ddc0e0b5SRichard Lowe  *	For x in [0.84375,1.25], let s = x - 1, and
6325c28e83SPiotr Jasiukajtis  *	c = 0.84506291151 rounded to single (24 bits)
6425c28e83SPiotr Jasiukajtis  *	   erf(x)  = c  + P1(s)/Q1(s)
65*ddc0e0b5SRichard Lowe  *	   erfc(x) = (1-c)  - P1(s)/Q1(s)
6625c28e83SPiotr Jasiukajtis  *	precision: |P1/Q1 - (erf(x)-c)| <= 2**-118.41
6725c28e83SPiotr Jasiukajtis  *
68*ddc0e0b5SRichard Lowe  *
69*ddc0e0b5SRichard Lowe  *	For x in [1.25,1.75], let s = x - 1.5, and
7025c28e83SPiotr Jasiukajtis  *	c = 0.95478588343 rounded to single (24 bits)
7125c28e83SPiotr Jasiukajtis  *	   erf(x)  = c  + P2(s)/Q2(s)
72*ddc0e0b5SRichard Lowe  *	   erfc(x) = (1-c)  - P2(s)/Q2(s)
7325c28e83SPiotr Jasiukajtis  *	precision: |P1/Q1 - (erf(x)-c)| <= 2**-123.83
7425c28e83SPiotr Jasiukajtis  *
7525c28e83SPiotr Jasiukajtis  *
7625c28e83SPiotr Jasiukajtis  *	For x in [1.75,16/3]
7725c28e83SPiotr Jasiukajtis  *	   erfc(x) = exp(-x*x)*(1/x)*R1(1/x)/S1(1/x)
7825c28e83SPiotr Jasiukajtis  *	   erf(x)  = 1 - erfc(x)
7925c28e83SPiotr Jasiukajtis  *	precision: absolute error of R1/S1 is bounded by 2**-124.03
8025c28e83SPiotr Jasiukajtis  *
8125c28e83SPiotr Jasiukajtis  *	For x in [16/3,107]
8225c28e83SPiotr Jasiukajtis  *	   erfc(x) = exp(-x*x)*(1/x)*R2(1/x)/S2(1/x)
8325c28e83SPiotr Jasiukajtis  *	   erf(x)  = 1 - erfc(x) (if x>=9 simple return erf(x)=1 with inexact)
8425c28e83SPiotr Jasiukajtis  *	precision: absolute error of R2/S2 is bounded by 2**-120.07
8525c28e83SPiotr Jasiukajtis  *
8625c28e83SPiotr Jasiukajtis  *	Else if inf > x >= 107
87*ddc0e0b5SRichard Lowe  *	   erf(x)  = 1 with inexact
8825c28e83SPiotr Jasiukajtis  *	   erfc(x) = 0 with underflow
89*ddc0e0b5SRichard Lowe  *
9025c28e83SPiotr Jasiukajtis  *	Special case:
9125c28e83SPiotr Jasiukajtis  *	   erf(inf)  = 1
9225c28e83SPiotr Jasiukajtis  *	   erfc(inf) = 0
9325c28e83SPiotr Jasiukajtis  */
9425c28e83SPiotr Jasiukajtis 
95*ddc0e0b5SRichard Lowe #pragma weak __erfl = erfl
96*ddc0e0b5SRichard Lowe #pragma weak __erfcl = erfcl
9725c28e83SPiotr Jasiukajtis 
9825c28e83SPiotr Jasiukajtis #include "libm.h"
9925c28e83SPiotr Jasiukajtis #include "longdouble.h"
10025c28e83SPiotr Jasiukajtis 
10125c28e83SPiotr Jasiukajtis static long double
10225c28e83SPiotr Jasiukajtis tiny	    = 1e-40L,
10325c28e83SPiotr Jasiukajtis nearunfl    = 1e-4000L,
10425c28e83SPiotr Jasiukajtis half	    = 0.5L,
10525c28e83SPiotr Jasiukajtis one	    = 1.0L,
10625c28e83SPiotr Jasiukajtis onehalf	    = 1.5L,
10725c28e83SPiotr Jasiukajtis L16_3	    = 16.0L/3.0L;
10825c28e83SPiotr Jasiukajtis /*
10925c28e83SPiotr Jasiukajtis  * Coefficients for even polynomial P for erf(x)=x+x*P(x^2) on [0,0.84375]
11025c28e83SPiotr Jasiukajtis  */
11125c28e83SPiotr Jasiukajtis static long double P[] = { 	/* 21 coeffs */
11225c28e83SPiotr Jasiukajtis    1.283791670955125738961589031215451715556e-0001L,
11325c28e83SPiotr Jasiukajtis   -3.761263890318375246320529677071815594603e-0001L,
11425c28e83SPiotr Jasiukajtis    1.128379167095512573896158903121205899135e-0001L,
11525c28e83SPiotr Jasiukajtis   -2.686617064513125175943235483344625046092e-0002L,
11625c28e83SPiotr Jasiukajtis    5.223977625442187842111846652980454568389e-0003L,
11725c28e83SPiotr Jasiukajtis   -8.548327023450852832546626271083862724358e-0004L,
11825c28e83SPiotr Jasiukajtis    1.205533298178966425102164715902231976672e-0004L,
11925c28e83SPiotr Jasiukajtis   -1.492565035840625097674944905027897838996e-0005L,
12025c28e83SPiotr Jasiukajtis    1.646211436588924733604648849172936692024e-0006L,
12125c28e83SPiotr Jasiukajtis   -1.636584469123491976815834704799733514987e-0007L,
12225c28e83SPiotr Jasiukajtis    1.480719281587897445302529007144770739305e-0008L,
12325c28e83SPiotr Jasiukajtis   -1.229055530170782843046467986464722047175e-0009L,
12425c28e83SPiotr Jasiukajtis    9.422759064320307357553954945760654341633e-0011L,
12525c28e83SPiotr Jasiukajtis   -6.711366846653439036162105104991433380926e-0012L,
12625c28e83SPiotr Jasiukajtis    4.463224090341893165100275380693843116240e-0013L,
12725c28e83SPiotr Jasiukajtis   -2.783513452582658245422635662559779162312e-0014L,
12825c28e83SPiotr Jasiukajtis    1.634227412586960195251346878863754661546e-0015L,
12925c28e83SPiotr Jasiukajtis   -9.060782672889577722765711455623117802795e-0017L,
13025c28e83SPiotr Jasiukajtis    4.741341801266246873412159213893613602354e-0018L,
13125c28e83SPiotr Jasiukajtis   -2.272417596497826188374846636534317381203e-0019L,
13225c28e83SPiotr Jasiukajtis    8.069088733716068462496835658928566920933e-0021L,
13325c28e83SPiotr Jasiukajtis };
13425c28e83SPiotr Jasiukajtis 
135*ddc0e0b5SRichard Lowe /*
13625c28e83SPiotr Jasiukajtis  * Rational erf(x) = ((float)0.84506291151) + P1(x-1)/Q1(x-1) on [0.84375,1.25]
13725c28e83SPiotr Jasiukajtis  */
13825c28e83SPiotr Jasiukajtis static long double C1   = (long double)((float)0.84506291151);
13925c28e83SPiotr Jasiukajtis static long double P1[] = { 	/*  12 top coeffs */
14025c28e83SPiotr Jasiukajtis   -2.362118560752659955654364917390741930316e-0003L,
14125c28e83SPiotr Jasiukajtis    4.129623379624420034078926610650759979146e-0001L,
14225c28e83SPiotr Jasiukajtis   -3.973857505403547283109417923182669976904e-0002L,
14325c28e83SPiotr Jasiukajtis    4.357503184084022439763567513078036755183e-0002L,
14425c28e83SPiotr Jasiukajtis    8.015593623388421371247676683754171456950e-0002L,
14525c28e83SPiotr Jasiukajtis   -1.034459310403352486685467221776778474602e-0002L,
14625c28e83SPiotr Jasiukajtis    5.671850295381046679675355719017720821383e-0003L,
14725c28e83SPiotr Jasiukajtis    1.219262563232763998351452194968781174318e-0003L,
14825c28e83SPiotr Jasiukajtis    5.390833481581033423020320734201065475098e-0004L,
14925c28e83SPiotr Jasiukajtis   -1.978853912815115495053119023517805528300e-0004L,
15025c28e83SPiotr Jasiukajtis    6.184234513953600118335017885706420552487e-0005L,
15125c28e83SPiotr Jasiukajtis   -5.331802711697810861017518515816271808286e-0006L,
15225c28e83SPiotr Jasiukajtis };
15325c28e83SPiotr Jasiukajtis static long double Q1[] = { 	/*  12 bottom coeffs with leading 1.0 hidden */
15425c28e83SPiotr Jasiukajtis    9.081506296064882195280178373107623196655e-0001L,
15525c28e83SPiotr Jasiukajtis    6.821049531968204097604392183650687642520e-0001L,
15625c28e83SPiotr Jasiukajtis    4.067869178233539502315055970743271822838e-0001L,
15725c28e83SPiotr Jasiukajtis    1.702332233546316765818144723063881095577e-0001L,
15825c28e83SPiotr Jasiukajtis    7.498098377690553934266423088708614219356e-0002L,
15925c28e83SPiotr Jasiukajtis    2.050154396918178697056927234366372760310e-0002L,
16025c28e83SPiotr Jasiukajtis    7.012988534031999899054782333851905939379e-0003L,
16125c28e83SPiotr Jasiukajtis    1.149904787014400354649843451234570731076e-0003L,
16225c28e83SPiotr Jasiukajtis    3.185620255011299476196039491205159718620e-0004L,
16325c28e83SPiotr Jasiukajtis    1.273405072153008775426376193374105840517e-0005L,
16425c28e83SPiotr Jasiukajtis    4.753866999959432971956781228148402971454e-0006L,
16525c28e83SPiotr Jasiukajtis   -1.002287602111660026053981728549540200683e-0006L,
16625c28e83SPiotr Jasiukajtis };
167*ddc0e0b5SRichard Lowe /*
16825c28e83SPiotr Jasiukajtis  * Rational erf(x) = ((float)0.95478588343) + P2(x-1.5)/Q2(x-1.5)
16925c28e83SPiotr Jasiukajtis  * on [1.25,1.75]
17025c28e83SPiotr Jasiukajtis  */
17125c28e83SPiotr Jasiukajtis static long double C2   = (long double)((float)0.95478588343);
17225c28e83SPiotr Jasiukajtis static long double P2[] = { 	/*  12 top coeffs */
17325c28e83SPiotr Jasiukajtis    1.131926304864446730135126164594785863512e-0002L,
17425c28e83SPiotr Jasiukajtis    1.273617996967754151544330055186210322832e-0001L,
17525c28e83SPiotr Jasiukajtis   -8.169980734667512519897816907190281143423e-0002L,
17625c28e83SPiotr Jasiukajtis    9.512267486090321197833634271787944271746e-0002L,
17725c28e83SPiotr Jasiukajtis   -2.394251569804872160005274999735914368170e-0002L,
17825c28e83SPiotr Jasiukajtis    1.108768660227528667525252333184520222905e-0002L,
17925c28e83SPiotr Jasiukajtis    3.527435492933902414662043314373277494221e-0004L,
18025c28e83SPiotr Jasiukajtis    4.946116273341953463584319006669474625971e-0004L,
18125c28e83SPiotr Jasiukajtis   -4.289851942513144714600285769022420962418e-0005L,
18225c28e83SPiotr Jasiukajtis    8.304719841341952705874781636002085119978e-0005L,
18325c28e83SPiotr Jasiukajtis   -1.040460226177309338781902252282849903189e-0005L,
18425c28e83SPiotr Jasiukajtis    2.122913331584921470381327583672044434087e-0006L,
18525c28e83SPiotr Jasiukajtis };
18625c28e83SPiotr Jasiukajtis static long double Q2[] = { 	/*  13 bottom coeffs with leading 1.0 hidden */
18725c28e83SPiotr Jasiukajtis    7.448815737306992749168727691042003832150e-0001L,
18825c28e83SPiotr Jasiukajtis    7.161813850236008294484744312430122188043e-0001L,
18925c28e83SPiotr Jasiukajtis    3.603134756584225766144922727405641236121e-0001L,
19025c28e83SPiotr Jasiukajtis    1.955811609133766478080550795194535852653e-0001L,
19125c28e83SPiotr Jasiukajtis    7.253059963716225972479693813787810711233e-0002L,
19225c28e83SPiotr Jasiukajtis    2.752391253757421424212770221541238324978e-0002L,
19325c28e83SPiotr Jasiukajtis    7.677654852085240257439050673446546828005e-0003L,
19425c28e83SPiotr Jasiukajtis    2.141102244555509687346497060326630061069e-0003L,
19525c28e83SPiotr Jasiukajtis    4.342123013830957093949563339130674364271e-0004L,
19625c28e83SPiotr Jasiukajtis    8.664587895570043348530991997272212150316e-0005L,
19725c28e83SPiotr Jasiukajtis    1.109201582511752087060167429397033701988e-0005L,
19825c28e83SPiotr Jasiukajtis    1.357834375781831062713347000030984364311e-0006L,
19925c28e83SPiotr Jasiukajtis    4.957746280594384997273090385060680016451e-0008L,
20025c28e83SPiotr Jasiukajtis };
20125c28e83SPiotr Jasiukajtis /*
20225c28e83SPiotr Jasiukajtis  * erfc(x) = exp(-x*x)/x * R1(1/x)/S1(1/x) on [1.75, 16/3]
20325c28e83SPiotr Jasiukajtis  */
20425c28e83SPiotr Jasiukajtis static long double R1[] = { 	/*  14 top coeffs */
20525c28e83SPiotr Jasiukajtis    4.630195122654315016370705767621550602948e+0006L,
20625c28e83SPiotr Jasiukajtis    1.257949521746494830700654204488675713628e+0007L,
20725c28e83SPiotr Jasiukajtis    1.704153822720260272814743497376181625707e+0007L,
20825c28e83SPiotr Jasiukajtis    1.502600568706061872381577539537315739943e+0007L,
20925c28e83SPiotr Jasiukajtis    9.543710793431995284827024445387333922861e+0006L,
21025c28e83SPiotr Jasiukajtis    4.589344808584091011652238164935949522427e+0006L,
21125c28e83SPiotr Jasiukajtis    1.714660662941745791190907071920671844289e+0006L,
21225c28e83SPiotr Jasiukajtis    5.034802147768798894307672256192466283867e+0005L,
21325c28e83SPiotr Jasiukajtis    1.162286400443554670553152110447126850725e+0005L,
21425c28e83SPiotr Jasiukajtis    2.086643834548901681362757308058660399137e+0004L,
21525c28e83SPiotr Jasiukajtis    2.839793161868140305907004392890348777338e+0003L,
21625c28e83SPiotr Jasiukajtis    2.786687241658423601778258694498655680778e+0002L,
21725c28e83SPiotr Jasiukajtis    1.779177837102695602425897452623985786464e+0001L,
21825c28e83SPiotr Jasiukajtis    5.641895835477470769043614623819144434731e-0001L,
21925c28e83SPiotr Jasiukajtis };
22025c28e83SPiotr Jasiukajtis static long double S1[] = { 	/* 15 bottom coeffs with leading 1.0 hidden */
22125c28e83SPiotr Jasiukajtis    4.630195122654331529595606896287596843110e+0006L,
22225c28e83SPiotr Jasiukajtis    1.780411093345512024324781084220509055058e+0007L,
22325c28e83SPiotr Jasiukajtis    3.250113097051800703707108623715776848283e+0007L,
22425c28e83SPiotr Jasiukajtis    3.737857099176755050912193712123489115755e+0007L,
22525c28e83SPiotr Jasiukajtis    3.029787497516578821459174055870781168593e+0007L,
22625c28e83SPiotr Jasiukajtis    1.833850619965384765005769632103205777227e+0007L,
22725c28e83SPiotr Jasiukajtis    8.562719999736915722210391222639186586498e+0006L,
22825c28e83SPiotr Jasiukajtis    3.139684562074658971315545539760008136973e+0006L,
22925c28e83SPiotr Jasiukajtis    9.106421313731384880027703627454366930945e+0005L,
23025c28e83SPiotr Jasiukajtis    2.085108342384266508613267136003194920001e+0005L,
23125c28e83SPiotr Jasiukajtis    3.723126272693120340730491416449539290600e+0004L,
23225c28e83SPiotr Jasiukajtis    5.049169878567344046145695360784436929802e+0003L,
23325c28e83SPiotr Jasiukajtis    4.944274532748010767670150730035392093899e+0002L,
23425c28e83SPiotr Jasiukajtis    3.153510608818213929982940249162268971412e+0001L,
23525c28e83SPiotr Jasiukajtis    1.0e00L,
23625c28e83SPiotr Jasiukajtis };
23725c28e83SPiotr Jasiukajtis 
23825c28e83SPiotr Jasiukajtis /*
23925c28e83SPiotr Jasiukajtis  * erfc(x) = exp(-x*x)/x * R2(1/x)/S2(1/x) on [16/3, 107]
24025c28e83SPiotr Jasiukajtis  */
24125c28e83SPiotr Jasiukajtis static long double R2[] = { 	/*  15 top coeffs in reverse order!!*/
24225c28e83SPiotr Jasiukajtis    2.447288012254302966796326587537136931669e+0005L,
24325c28e83SPiotr Jasiukajtis    8.768592567189861896653369912716538739016e+0005L,
24425c28e83SPiotr Jasiukajtis    1.552293152581780065761497908005779524953e+0006L,
24525c28e83SPiotr Jasiukajtis    1.792075924835942935864231657504259926729e+0006L,
24625c28e83SPiotr Jasiukajtis    1.504001463155897344947500222052694835875e+0006L,
24725c28e83SPiotr Jasiukajtis    9.699485556326891411801230186016013019935e+0005L,
24825c28e83SPiotr Jasiukajtis    4.961449933661807969863435013364796037700e+0005L,
24925c28e83SPiotr Jasiukajtis    2.048726544693474028061176764716228273791e+0005L,
25025c28e83SPiotr Jasiukajtis    6.891532964330949722479061090551896886635e+0004L,
25125c28e83SPiotr Jasiukajtis    1.888014709010307507771964047905823237985e+0004L,
25225c28e83SPiotr Jasiukajtis    4.189692064988957745054734809642495644502e+0003L,
25325c28e83SPiotr Jasiukajtis    7.362346487427048068212968889642741734621e+0002L,
25425c28e83SPiotr Jasiukajtis    9.980359714211411423007641056580813116207e+0001L,
25525c28e83SPiotr Jasiukajtis    9.426910895135379181107191962193485174159e+0000L,
25625c28e83SPiotr Jasiukajtis    5.641895835477562869480794515623601280429e-0001L,
25725c28e83SPiotr Jasiukajtis };
25825c28e83SPiotr Jasiukajtis static long double S2[] = { 	/* 16 coefficients */
25925c28e83SPiotr Jasiukajtis    2.447282203601902971246004716790604686880e+0005L,
26025c28e83SPiotr Jasiukajtis    1.153009852759385309367759460934808489833e+0006L,
26125c28e83SPiotr Jasiukajtis    2.608580649612639131548966265078663384849e+0006L,
26225c28e83SPiotr Jasiukajtis    3.766673917346623308850202792390569025740e+0006L,
26325c28e83SPiotr Jasiukajtis    3.890566255138383910789924920541335370691e+0006L,
26425c28e83SPiotr Jasiukajtis    3.052882073900746207613166259994150527732e+0006L,
26525c28e83SPiotr Jasiukajtis    1.885574519970380988460241047248519418407e+0006L,
26625c28e83SPiotr Jasiukajtis    9.369722034759943185851450846811445012922e+0005L,
26725c28e83SPiotr Jasiukajtis    3.792278350536686111444869752624492443659e+0005L,
26825c28e83SPiotr Jasiukajtis    1.257750606950115799965366001773094058720e+0005L,
26925c28e83SPiotr Jasiukajtis    3.410830600242369370645608634643620355058e+0004L,
27025c28e83SPiotr Jasiukajtis    7.513984469742343134851326863175067271240e+0003L,
27125c28e83SPiotr Jasiukajtis    1.313296320593190002554779998138695507840e+0003L,
27225c28e83SPiotr Jasiukajtis    1.773972700887629157006326333696896516769e+0002L,
27325c28e83SPiotr Jasiukajtis    1.670876451822586800422009013880457094162e+0001L,
27425c28e83SPiotr Jasiukajtis    1.000L,
27525c28e83SPiotr Jasiukajtis };
27625c28e83SPiotr Jasiukajtis 
erfl(x)277*ddc0e0b5SRichard Lowe long double erfl(x)
27825c28e83SPiotr Jasiukajtis long double x;
27925c28e83SPiotr Jasiukajtis {
28025c28e83SPiotr Jasiukajtis 	long double erfcl(long double),s,y,t;
28125c28e83SPiotr Jasiukajtis 
28225c28e83SPiotr Jasiukajtis 	if (!finitel(x)) {
28325c28e83SPiotr Jasiukajtis 	    if (x != x) return x+x; 	/* NaN */
28425c28e83SPiotr Jasiukajtis 	    return copysignl(one,x);	/* return +-1.0 is x=Inf */
28525c28e83SPiotr Jasiukajtis 	}
28625c28e83SPiotr Jasiukajtis 
28725c28e83SPiotr Jasiukajtis 	y = fabsl(x);
28825c28e83SPiotr Jasiukajtis 	if (y <= 0.84375L) {
28925c28e83SPiotr Jasiukajtis 	    if (y<=tiny) return x+P[0]*x;
29025c28e83SPiotr Jasiukajtis 	    s = y*y;
29125c28e83SPiotr Jasiukajtis 	    t = __poly_libmq(s,21,P);
29225c28e83SPiotr Jasiukajtis 	    return  x+x*t;
29325c28e83SPiotr Jasiukajtis 	}
29425c28e83SPiotr Jasiukajtis 	if (y<=1.25L) {
29525c28e83SPiotr Jasiukajtis 	    s = y-one;
29625c28e83SPiotr Jasiukajtis 	    t = C1+__poly_libmq(s,12,P1)/(one+s*__poly_libmq(s,12,Q1));
29725c28e83SPiotr Jasiukajtis 	    return (signbitl(x))? -t: t;
29825c28e83SPiotr Jasiukajtis 	} else if (y<=1.75L) {
29925c28e83SPiotr Jasiukajtis 	    s = y-onehalf;
30025c28e83SPiotr Jasiukajtis 	    t = C2+__poly_libmq(s,12,P2)/(one+s*__poly_libmq(s,13,Q2));
30125c28e83SPiotr Jasiukajtis 	    return (signbitl(x))? -t: t;
302*ddc0e0b5SRichard Lowe 	}
30325c28e83SPiotr Jasiukajtis 	if (y<=9.0L) t = erfcl(y); else t = tiny;
30425c28e83SPiotr Jasiukajtis 	return (signbitl(x))? t-one: one-t;
30525c28e83SPiotr Jasiukajtis }
30625c28e83SPiotr Jasiukajtis 
erfcl(x)30725c28e83SPiotr Jasiukajtis long double erfcl(x)
30825c28e83SPiotr Jasiukajtis long double x;
30925c28e83SPiotr Jasiukajtis {
31025c28e83SPiotr Jasiukajtis 	long double erfl(long double),s,y,t;
31125c28e83SPiotr Jasiukajtis 
31225c28e83SPiotr Jasiukajtis 	if (!finitel(x)) {
31325c28e83SPiotr Jasiukajtis 	    if (x != x) return x+x; 	/* NaN */
31425c28e83SPiotr Jasiukajtis 	    				/* return 2.0 if x= -inf
31525c28e83SPiotr Jasiukajtis 						  0.0 if x= +inf */
31625c28e83SPiotr Jasiukajtis 	    if (x<0.0L) return 2.0L; else return 0.0L;
31725c28e83SPiotr Jasiukajtis 	}
31825c28e83SPiotr Jasiukajtis 
31925c28e83SPiotr Jasiukajtis 	if (x <= 0.84375L) {
32025c28e83SPiotr Jasiukajtis 	    if (x<=0.25) return one-erfl(x);
32125c28e83SPiotr Jasiukajtis 	    s = x*x;
32225c28e83SPiotr Jasiukajtis 	    t = half-x;
32325c28e83SPiotr Jasiukajtis 	    t = t - x*__poly_libmq(s,21,P);
32425c28e83SPiotr Jasiukajtis 	    return  half+t;
32525c28e83SPiotr Jasiukajtis 	}
32625c28e83SPiotr Jasiukajtis 	if (x<=1.25L) {
32725c28e83SPiotr Jasiukajtis 	    s = x-one;
32825c28e83SPiotr Jasiukajtis 	    t = one-C1;
32925c28e83SPiotr Jasiukajtis 	    return t - __poly_libmq(s,12,P1)/(one+s*__poly_libmq(s,12,Q1));
33025c28e83SPiotr Jasiukajtis 	} else if (x<=1.75L) {
33125c28e83SPiotr Jasiukajtis 	    s = x-onehalf;
33225c28e83SPiotr Jasiukajtis 	    t = one-C2;
33325c28e83SPiotr Jasiukajtis 	    return t - __poly_libmq(s,12,P2)/(one+s*__poly_libmq(s,13,Q2));
334*ddc0e0b5SRichard Lowe 	}
33525c28e83SPiotr Jasiukajtis 	if (x>=107.0L) return nearunfl*nearunfl;		/* underflow */
33625c28e83SPiotr Jasiukajtis 	else if (x >= L16_3) {
33725c28e83SPiotr Jasiukajtis 	    y = __poly_libmq(x,15,R2);
33825c28e83SPiotr Jasiukajtis 	    t = y/__poly_libmq(x,16,S2);
33925c28e83SPiotr Jasiukajtis 	} else {
34025c28e83SPiotr Jasiukajtis 	    y = __poly_libmq(x,14,R1);
34125c28e83SPiotr Jasiukajtis 	    t = y/__poly_libmq(x,15,S1);
34225c28e83SPiotr Jasiukajtis 	}
34325c28e83SPiotr Jasiukajtis 	/* see comment in ../Q/erfl.c */
34425c28e83SPiotr Jasiukajtis 	y = x;
34525c28e83SPiotr Jasiukajtis 	*(int*)&y = 0;
34625c28e83SPiotr Jasiukajtis 	t *= expl(-y*y)*expl(-(x-y)*(x+y));
34725c28e83SPiotr Jasiukajtis 	return  t;
34825c28e83SPiotr Jasiukajtis }
349