Print this page
11210 libm should be cstyle(1ONBLD) clean
*** 20,29 ****
--- 20,30 ----
*/
/*
* Copyright 2011 Nexenta Systems, Inc. All rights reserved.
*/
+
/*
* Copyright 2006 Sun Microsystems, Inc. All rights reserved.
* Use is subject to license terms.
*/
*** 97,113 ****
#pragma weak __erfcl = erfcl
#include "libm.h"
#include "longdouble.h"
! static const long double
! tiny = 1e-40L,
nearunfl = 1e-4000L,
half = 0.5L,
one = 1.0L,
onehalf = 1.5L,
! L16_3 = 16.0L/3.0L;
/*
* Coefficients for even polynomial P for erf(x)=x+x*P(x^2) on [0,0.84375]
*/
static const long double P[] = { /* 21 coeffs */
1.283791670955125738961589031215451715556e-0001L,
--- 98,114 ----
#pragma weak __erfcl = erfcl
#include "libm.h"
#include "longdouble.h"
! static const long double tiny = 1e-40L,
nearunfl = 1e-4000L,
half = 0.5L,
one = 1.0L,
onehalf = 1.5L,
! L16_3 = 16.0L / 3.0L;
!
/*
* Coefficients for even polynomial P for erf(x)=x+x*P(x^2) on [0,0.84375]
*/
static const long double P[] = { /* 21 coeffs */
1.283791670955125738961589031215451715556e-0001L,
*** 149,158 ****
--- 150,160 ----
5.390833481581033423020320734201065475098e-0004L,
-1.978853912815115495053119023517805528300e-0004L,
6.184234513953600118335017885706420552487e-0005L,
-5.331802711697810861017518515816271808286e-0006L,
};
+
static const long double Q1[] = { /* 12 bottom coeffs with leading 1.0 hidden */
9.081506296064882195280178373107623196655e-0001L,
6.821049531968204097604392183650687642520e-0001L,
4.067869178233539502315055970743271822838e-0001L,
1.702332233546316765818144723063881095577e-0001L,
*** 163,172 ****
--- 165,175 ----
3.185620255011299476196039491205159718620e-0004L,
1.273405072153008775426376193374105840517e-0005L,
4.753866999959432971956781228148402971454e-0006L,
-1.002287602111660026053981728549540200683e-0006L,
};
+
/*
* Rational erf(x) = ((float)0.95478588343) + P2(x-1.5)/Q2(x-1.5)
* on [1.25,1.75]
*/
static const long double C2 = (long double)((float)0.95478588343);
*** 182,191 ****
--- 185,195 ----
-4.289851942513144714600285769022420962418e-0005L,
8.304719841341952705874781636002085119978e-0005L,
-1.040460226177309338781902252282849903189e-0005L,
2.122913331584921470381327583672044434087e-0006L,
};
+
static const long double Q2[] = { /* 13 bottom coeffs with leading 1.0 hidden */
7.448815737306992749168727691042003832150e-0001L,
7.161813850236008294484744312430122188043e-0001L,
3.603134756584225766144922727405641236121e-0001L,
1.955811609133766478080550795194535852653e-0001L,
*** 197,206 ****
--- 201,211 ----
8.664587895570043348530991997272212150316e-0005L,
1.109201582511752087060167429397033701988e-0005L,
1.357834375781831062713347000030984364311e-0006L,
4.957746280594384997273090385060680016451e-0008L,
};
+
/*
* erfc(x) = exp(-x*x)/x * R1(1/x)/S1(1/x) on [1.75, 16/3]
*/
static const long double R1[] = { /* 14 top coeffs */
4.630195122654315016370705767621550602948e+0006L,
*** 216,225 ****
--- 221,231 ----
2.839793161868140305907004392890348777338e+0003L,
2.786687241658423601778258694498655680778e+0002L,
1.779177837102695602425897452623985786464e+0001L,
5.641895835477470769043614623819144434731e-0001L,
};
+
static const long double S1[] = { /* 15 bottom coeffs with leading 1.0 hidden */
4.630195122654331529595606896287596843110e+0006L,
1.780411093345512024324781084220509055058e+0007L,
3.250113097051800703707108623715776848283e+0007L,
3.737857099176755050912193712123489115755e+0007L,
*** 237,247 ****
};
/*
* erfc(x) = exp(-x*x)/x * R2(1/x)/S2(1/x) on [16/3, 107]
*/
! static const long double R2[] = { /* 15 top coeffs in reverse order!!*/
2.447288012254302966796326587537136931669e+0005L,
8.768592567189861896653369912716538739016e+0005L,
1.552293152581780065761497908005779524953e+0006L,
1.792075924835942935864231657504259926729e+0006L,
1.504001463155897344947500222052694835875e+0006L,
--- 243,253 ----
};
/*
* erfc(x) = exp(-x*x)/x * R2(1/x)/S2(1/x) on [16/3, 107]
*/
! static const long double R2[] = { /* 15 top coeffs in reverse order!! */
2.447288012254302966796326587537136931669e+0005L,
8.768592567189861896653369912716538739016e+0005L,
1.552293152581780065761497908005779524953e+0006L,
1.792075924835942935864231657504259926729e+0006L,
1.504001463155897344947500222052694835875e+0006L,
*** 254,263 ****
--- 260,270 ----
7.362346487427048068212968889642741734621e+0002L,
9.980359714211411423007641056580813116207e+0001L,
9.426910895135379181107191962193485174159e+0000L,
5.641895835477562869480794515623601280429e-0001L,
};
+
static const long double S2[] = { /* 16 coefficients */
2.447282203601902971246004716790604686880e+0005L,
1.153009852759385309367759460934808489833e+0006L,
2.608580649612639131548966265078663384849e+0006L,
3.766673917346623308850202792390569025740e+0006L,
*** 273,347 ****
1.773972700887629157006326333696896516769e+0002L,
1.670876451822586800422009013880457094162e+0001L,
1.000L,
};
! long double erfl(x)
! long double x;
{
! long double s,y,t;
if (!finitel(x)) {
! if (x != x) return x+x; /* NaN */
! return copysignl(one,x); /* return +-1.0 is x=Inf */
}
y = fabsl(x);
if (y <= 0.84375L) {
! if (y<=tiny) return x+P[0]*x;
! s = y*y;
! t = __poly_libmq(s,21,P);
! return x+x*t;
}
! if (y<=1.25L) {
! s = y-one;
! t = C1+__poly_libmq(s,12,P1)/(one+s*__poly_libmq(s,12,Q1));
! return (signbitl(x))? -t: t;
! } else if (y<=1.75L) {
! s = y-onehalf;
! t = C2+__poly_libmq(s,12,P2)/(one+s*__poly_libmq(s,13,Q2));
! return (signbitl(x))? -t: t;
}
! if (y<=9.0L) t = erfcl(y); else t = tiny;
! return (signbitl(x))? t-one: one-t;
}
! long double erfcl(x)
! long double x;
{
! long double s,y,t;
if (!finitel(x)) {
! if (x != x) return x+x; /* NaN */
/* return 2.0 if x= -inf; 0.0 if x= +inf */
! if (x < 0.0L) return 2.0L; else return 0.0L;
}
if (x <= 0.84375L) {
! if (x<=0.25) return one-erfl(x);
! s = x*x;
! t = half-x;
! t = t - x*__poly_libmq(s,21,P);
! return half+t;
}
! if (x<=1.25L) {
! s = x-one;
! t = one-C1;
! return t - __poly_libmq(s,12,P1)/(one+s*__poly_libmq(s,12,Q1));
! } else if (x<=1.75L) {
! s = x-onehalf;
! t = one-C2;
! return t - __poly_libmq(s,12,P2)/(one+s*__poly_libmq(s,13,Q2));
}
! if (x>=107.0L) return nearunfl*nearunfl; /* underflow */
! else if (x >= L16_3) {
! y = __poly_libmq(x,15,R2);
! t = y/__poly_libmq(x,16,S2);
} else {
! y = __poly_libmq(x,14,R1);
! t = y/__poly_libmq(x,15,S1);
}
/*
* Note that exp(-x*x+d) = exp(-x*x)*exp(d), so to compute
* exp(-x*x) with a small relative error, we need to compute
* -x*x with a small absolute error. To this end, we set y
* equal to the leading part of x but with enough trailing
--- 280,380 ----
1.773972700887629157006326333696896516769e+0002L,
1.670876451822586800422009013880457094162e+0001L,
1.000L,
};
! long double
! erfl(long double x)
{
! long double s, y, t;
if (!finitel(x)) {
! if (x != x)
! return (x + x); /* NaN */
!
! return (copysignl(one, x)); /* return +-1.0 is x=Inf */
}
y = fabsl(x);
+
if (y <= 0.84375L) {
! if (y <= tiny)
! return (x + P[0] * x);
!
! s = y * y;
! t = __poly_libmq(s, 21, P);
! return (x + x * t);
}
!
! if (y <= 1.25L) {
! s = y - one;
! t = C1 + __poly_libmq(s, 12, P1) / (one + s * __poly_libmq(s,
! 12, Q1));
! return ((signbitl(x)) ? -t : t);
! } else if (y <= 1.75L) {
! s = y - onehalf;
! t = C2 + __poly_libmq(s, 12, P2) / (one + s * __poly_libmq(s,
! 13, Q2));
! return ((signbitl(x)) ? -t : t);
}
!
! if (y <= 9.0L)
! t = erfcl(y);
! else
! t = tiny;
!
! return ((signbitl(x)) ? t - one : one - t);
}
! long double
! erfcl(long double x)
{
! long double s, y, t;
if (!finitel(x)) {
! if (x != x)
! return (x + x); /* NaN */
!
/* return 2.0 if x= -inf; 0.0 if x= +inf */
! if (x < 0.0L)
! return (2.0L);
! else
! return (0.0L);
}
if (x <= 0.84375L) {
! if (x <= 0.25)
! return (one - erfl(x));
!
! s = x * x;
! t = half - x;
! t = t - x * __poly_libmq(s, 21, P);
! return (half + t);
}
!
! if (x <= 1.25L) {
! s = x - one;
! t = one - C1;
! return (t - __poly_libmq(s, 12, P1) / (one + s * __poly_libmq(s,
! 12, Q1)));
! } else if (x <= 1.75L) {
! s = x - onehalf;
! t = one - C2;
! return (t - __poly_libmq(s, 12, P2) / (one + s * __poly_libmq(s,
! 13, Q2)));
}
!
! if (x >= 107.0L) {
! return (nearunfl * nearunfl); /* underflow */
! } else if (x >= L16_3) {
! y = __poly_libmq(x, 15, R2);
! t = y / __poly_libmq(x, 16, S2);
} else {
! y = __poly_libmq(x, 14, R1);
! t = y / __poly_libmq(x, 15, S1);
}
+
/*
* Note that exp(-x*x+d) = exp(-x*x)*exp(d), so to compute
* exp(-x*x) with a small relative error, we need to compute
* -x*x with a small absolute error. To this end, we set y
* equal to the leading part of x but with enough trailing
*** 358,366 ****
* small but not so large that the conversion to int overflows.
* When long double arithmetic is slow, however, the following
* non-portable code is preferable.
*/
y = x;
! *(2+(int*)&y) = *(3+(int*)&y) = 0;
! t *= expl(-y*y)*expl(-(x-y)*(x+y));
! return t;
}
--- 391,399 ----
* small but not so large that the conversion to int overflows.
* When long double arithmetic is slow, however, the following
* non-portable code is preferable.
*/
y = x;
! *(2 + (int *)&y) = *(3 + (int *)&y) = 0;
! t *= expl(-y * y) * expl(-(x - y) * (x + y));
! return (t);
}