riscv
/
musl
mirror of https://git.musl-libc.org/git/musl
1
0
Fork 0
Code Issues Projects Releases Wiki Activity
Browse Source

Merge remote-tracking branch 'nsz/math'

rs-1.0
Rich Felker 14 years ago
parent
dd95916382 121e3a38a1
commit
21f01e1bf9
3 changed files with 207 additions and 297 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 155
      src/math/erf.c
  2. 154
      src/math/erff.c
  3. 195
      src/math/erfl.c

155
src/math/erf.c
View File

@ -106,13 +106,10 @@
#include "libm.h"
static const double
tiny = 1e-300,
/* c = (float)0.84506291151 */
erx = 8.45062911510467529297e-01, /* 0x3FEB0AC1, 0x60000000 */
/*
* Coefficients for approximation to erf on [0,0.84375]
*/
efx = 1.28379167095512586316e-01, /* 0x3FC06EBA, 0x8214DB69 */
efx8 = 1.02703333676410069053e+00, /* 0x3FF06EBA, 0x8214DB69 */
pp0 = 1.28379167095512558561e-01, /* 0x3FC06EBA, 0x8214DB68 */
pp1 = -3.25042107247001499370e-01, /* 0xBFD4CD7D, 0x691CB913 */
@ -177,52 +174,31 @@ sb5 = 2.55305040643316442583e+03, /* 0x40A3F219, 0xCEDF3BE6 */
sb6 = 4.74528541206955367215e+02, /* 0x407DA874, 0xE79FE763 */
sb7 = -2.24409524465858183362e+01; /* 0xC03670E2, 0x42712D62 */
double erf(double x)
static double erfc1(double x)
{
int32_t hx,ix,i;
double R,S,P,Q,s,y,z,r;
double s,P,Q;
s = fabs(x) - 1;
P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
Q = 1+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
return 1 - erx - P/Q;
}
static double erfc2(uint32_t ix, double x)
{
double s,z,R,S;
if (ix < 0x3ff40000) /* |x| < 1.25 */
return erfc1(x);
GET_HIGH_WORD(hx, x);
ix = hx & 0x7fffffff;
if (ix >= 0x7ff00000) {
/* erf(nan)=nan, erf(+-inf)=+-1 */
i = ((uint32_t)hx>>31)<<1;
return (double)(1-i) + 1.0/x;
}
if (ix < 0x3feb0000) { /* |x|<0.84375 */
if (ix < 0x3e300000) { /* |x|<2**-28 */
if (ix < 0x00800000)
/* avoid underflow */
return 0.125*(8.0*x + efx8*x);
return x + efx*x;
}
z = x*x;
r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
s = 1.0+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
y = r/s;
return x + x*y;
}
if (ix < 0x3ff40000) { /* 0.84375 <= |x| < 1.25 */
s = fabs(x)-1.0;
P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
Q = 1.0+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
if (hx >= 0)
return erx + P/Q;
return -erx - P/Q;
}
if (ix >= 0x40180000) { /* inf > |x| >= 6 */
if (hx >= 0)
return 1.0 - tiny;
return tiny - 1.0;
}
x = fabs(x);
s = 1.0/(x*x);
if (ix < 0x4006DB6E) { /* |x| < 1/0.35 */
s = 1/(x*x);
if (ix < 0x4006db6d) { /* |x| < 1/.35 ~ 2.85714 */
R = ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
ra5+s*(ra6+s*ra7))))));
S = 1.0+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
sa5+s*(sa6+s*(sa7+s*sa8)))))));
} else { /* |x| >= 1/0.35 */
} else { /* |x| > 1/.35 */
R = rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
rb5+s*rb6)))));
S = 1.0+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
@ -230,22 +206,52 @@ double erf(double x)
}
z = x;
SET_LOW_WORD(z,0);
r = exp(-z*z-0.5625)*exp((z-x)*(z+x)+R/S);
if (hx >= 0)
return 1.0 - r/x;
return r/x - 1.0;
return exp(-z*z-0.5625)*exp((z-x)*(z+x)+R/S)/x;
}
double erf(double x)
{
double r,s,z,y;
uint32_t ix;
int sign;
GET_HIGH_WORD(ix, x);
sign = ix>>31;
ix &= 0x7fffffff;
if (ix >= 0x7ff00000) {
/* erf(nan)=nan, erf(+-inf)=+-1 */
return 1-2*sign + 1/x;
}
if (ix < 0x3feb0000) { /* |x| < 0.84375 */
if (ix < 0x3e300000) { /* |x| < 2**-28 */
/* avoid underflow */
return 0.125*(8*x + efx8*x);
}
z = x*x;
r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
s = 1.0+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
y = r/s;
return x + x*y;
}
if (ix < 0x40180000) /* 0.84375 <= |x| < 6 */
y = 1 - erfc2(ix,x);
else
y = 1 - 0x1p-1022;
return sign ? -y : y;
}
double erfc(double x)
{
int32_t hx,ix;
double R,S,P,Q,s,y,z,r;
double r,s,z,y;
uint32_t ix;
int sign;
GET_HIGH_WORD(hx, x);
ix = hx & 0x7fffffff;
GET_HIGH_WORD(ix, x);
sign = ix>>31;
ix &= 0x7fffffff;
if (ix >= 0x7ff00000) {
/* erfc(nan)=nan, erfc(+-inf)=0,2 */
return (double)(((uint32_t)hx>>31)<<1) + 1.0/x;
return 2*sign + 1/x;
}
if (ix < 0x3feb0000) { /* |x| < 0.84375 */
if (ix < 0x3c700000) /* |x| < 2**-56 */
@ -254,50 +260,13 @@ double erfc(double x)
r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
s = 1.0+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
y = r/s;
if (hx < 0x3fd00000) { /* x < 1/4 */
if (sign || ix < 0x3fd00000) { /* x < 1/4 */
return 1.0 - (x+x*y);
} else {
r = x*y;
r += x - 0.5;
return 0.5 - r ;
}
}
if (ix < 0x3ff40000) { /* 0.84375 <= |x| < 1.25 */
s = fabs(x)-1.0;
P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
Q = 1.0+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
if (hx >= 0) {
z = 1.0-erx;
return z - P/Q;
} else {
z = erx+P/Q;
return 1.0 + z;
}
return 0.5 - (x - 0.5 + x*y);
}
if (ix < 0x403c0000) { /* |x| < 28 */
x = fabs(x);
s = 1.0/(x*x);
if (ix < 0x4006DB6D) { /* |x| < 1/.35 ~ 2.857143*/
R = ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
ra5+s*(ra6+s*ra7))))));
S = 1.0+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
sa5+s*(sa6+s*(sa7+s*sa8)))))));
} else { /* |x| >= 1/.35 ~ 2.857143 */
if (hx < 0 && ix >= 0x40180000) /* x < -6 */
return 2.0 - tiny;
R = rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
rb5+s*rb6)))));
S = 1.0+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
sb5+s*(sb6+s*sb7))))));
}
z = x;
SET_LOW_WORD(z, 0);
r = exp(-z*z-0.5625)*exp((z-x)*(z+x)+R/S);
if (hx > 0)
return r/x;
return 2.0 - r/x;
if (ix < 0x403c0000) { /* 0.84375 <= |x| < 28 */
return sign ? 2 - erfc2(ix,x) : erfc2(ix,x);
}
if (hx > 0)
return tiny*tiny;
return 2.0 - tiny;
return sign ? 2 - 0x1p-1022 : 0x1p-1022*0x1p-1022;
}

154
src/math/erff.c
View File

@ -16,13 +16,10 @@
#include "libm.h"
static const float
tiny = 1e-30,
/* c = (subfloat)0.84506291151 */
erx = 8.4506291151e-01, /* 0x3f58560b */
/*
* Coefficients for approximation to erf on [0,0.84375]
*/
efx = 1.2837916613e-01, /* 0x3e0375d4 */
efx8 = 1.0270333290e+00, /* 0x3f8375d4 */
pp0 = 1.2837916613e-01, /* 0x3e0375d4 */
pp1 = -3.2504209876e-01, /* 0xbea66beb */
@ -87,47 +84,26 @@ sb5 = 2.5530502930e+03, /* 0x451f90ce */
sb6 = 4.7452853394e+02, /* 0x43ed43a7 */
sb7 = -2.2440952301e+01; /* 0xc1b38712 */
float erff(float x)
static float erfc1(float x)
{
int32_t hx,ix,i;
float R,S,P,Q,s,y,z,r;
float s,P,Q;
s = fabsf(x) - 1;
P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
Q = 1+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
return 1 - erx - P/Q;
}
static float erfc2(uint32_t ix, float x)
{
float s,z,R,S;
if (ix < 0x3fa00000) /* |x| < 1.25 */
return erfc1(x);
GET_FLOAT_WORD(hx, x);
ix = hx & 0x7fffffff;
if (ix >= 0x7f800000) {
/* erf(nan)=nan, erf(+-inf)=+-1 */
i = ((uint32_t)hx>>31)<<1;
return (float)(1-i)+1.0f/x;
}
if (ix < 0x3f580000) { /* |x| < 0.84375 */
if (ix < 0x31800000) { /* |x| < 2**-28 */
if (ix < 0x04000000)
/*avoid underflow */
return 0.125f*(8.0f*x + efx8*x);
return x + efx*x;
}
z = x*x;
r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
s = 1.0f+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
y = r/s;
return x + x*y;
}
if (ix < 0x3fa00000) { /* 0.84375 <= |x| < 1.25 */
s = fabsf(x)-1.0f;
P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
Q = 1.0f+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
if (hx >= 0)
return erx + P/Q;
return -erx - P/Q;
}
if (ix >= 0x40c00000) { /* inf > |x| >= 6 */
if (hx >= 0)
return 1.0f - tiny;
return tiny - 1.0f;
}
x = fabsf(x);
s = 1.0f/(x*x);
if (ix < 0x4036DB6E) { /* |x| < 1/0.35 */
s = 1/(x*x);
if (ix < 0x4036db6d) { /* |x| < 1/0.35 */
R = ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
ra5+s*(ra6+s*ra7))))));
S = 1.0f+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
@ -139,23 +115,53 @@ float erff(float x)
sb5+s*(sb6+s*sb7))))));
}
GET_FLOAT_WORD(ix, x);
SET_FLOAT_WORD(z, ix&0xfffff000);
r = expf(-z*z - 0.5625f) * expf((z-x)*(z+x) + R/S);
if (hx >= 0)
return 1.0f - r/x;
return r/x - 1.0f;
SET_FLOAT_WORD(z, ix&0xffffe000);
return expf(-z*z - 0.5625f) * expf((z-x)*(z+x) + R/S)/x;
}
float erff(float x)
{
float r,s,z,y;
uint32_t ix;
int sign;
GET_FLOAT_WORD(ix, x);
sign = ix>>31;
ix &= 0x7fffffff;
if (ix >= 0x7f800000) {
/* erf(nan)=nan, erf(+-inf)=+-1 */
return 1-2*sign + 1/x;
}
if (ix < 0x3f580000) { /* |x| < 0.84375 */
if (ix < 0x31800000) { /* |x| < 2**-28 */
/*avoid underflow */
return 0.125f*(8*x + efx8*x);
}
z = x*x;
r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
s = 1+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
y = r/s;
return x + x*y;
}
if (ix < 0x40c00000) /* |x| < 6 */
y = 1 - erfc2(ix,x);
else
y = 1 - 0x1p-120f;
return sign ? -y : y;
}
float erfcf(float x)
{
int32_t hx,ix;
float R,S,P,Q,s,y,z,r;
float r,s,z,y;
uint32_t ix;
int sign;
GET_FLOAT_WORD(hx, x);
ix = hx & 0x7fffffff;
GET_FLOAT_WORD(ix, x);
sign = ix>>31;
ix &= 0x7fffffff;
if (ix >= 0x7f800000) {
/* erfc(nan)=nan, erfc(+-inf)=0,2 */
return (float)(((uint32_t)hx>>31)<<1) + 1.0f/x;
return 2*sign + 1/x;
}
if (ix < 0x3f580000) { /* |x| < 0.84375 */
@ -165,50 +171,12 @@ float erfcf(float x)
r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
s = 1.0f+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
y = r/s;
if (hx < 0x3e800000) { /* x<1/4 */
if (sign || ix < 0x3e800000) /* x < 1/4 */
return 1.0f - (x+x*y);
} else {
r = x*y;
r += (x-0.5f);
return 0.5f - r ;
}
}
if (ix < 0x3fa00000) { /* 0.84375 <= |x| < 1.25 */
s = fabsf(x)-1.0f;
P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
Q = 1.0f+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
if(hx >= 0) {
z = 1.0f - erx;
return z - P/Q;
} else {
z = erx + P/Q;
return 1.0f + z;
}
return 0.5f - (x - 0.5f + x*y);
}
if (ix < 0x41e00000) { /* |x| < 28 */
x = fabsf(x);
s = 1.0f/(x*x);
if (ix < 0x4036DB6D) { /* |x| < 1/.35 ~ 2.857143*/
R = ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
ra5+s*(ra6+s*ra7))))));
S = 1.0f+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
sa5+s*(sa6+s*(sa7+s*sa8)))))));
} else { /* |x| >= 1/.35 ~ 2.857143 */
if (hx < 0 && ix >= 0x40c00000) /* x < -6 */
return 2.0f-tiny;
R = rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
rb5+s*rb6)))));
S = 1.0f+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
sb5+s*(sb6+s*sb7))))));
}
GET_FLOAT_WORD(ix, x);
SET_FLOAT_WORD(z, ix&0xfffff000);
r = expf(-z*z - 0.5625f) * expf((z-x)*(z+x) + R/S);
if (hx > 0)
return r/x;
return 2.0f - r/x;
return sign ? 2 - erfc2(ix,x) : erfc2(ix,x);
}
if (hx > 0)
return tiny*tiny;
return 2.0f - tiny;
return sign ? 2 - 0x1p-120f : 0x1p-120f*0x1p-120f;
}

195
src/math/erfl.c
View File

@ -107,15 +107,11 @@ long double erfl(long double x)
}
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
static const long double
tiny = 1e-4931L,
/* c = (float)0.84506291151 */
erx = 0.845062911510467529296875L,
/*
* Coefficients for approximation to erf on [0,0.84375]
*/
/* 2/sqrt(pi) - 1 */
efx = 1.2837916709551257389615890312154517168810E-1L,
/* 8 * (2/sqrt(pi) - 1) */
efx8 = 1.0270333367641005911692712249723613735048E0L,
pp[6] = {
@ -239,84 +235,108 @@ sc[] = {
/* 1.000000000000000000000000000000000000000E0 */
};
long double erfl(long double x)
static long double erfc1(long double x)
{
long double R, S, P, Q, s, y, z, r;
int32_t ix, i;
uint32_t se, i0, i1;
long double s,P,Q;
GET_LDOUBLE_WORDS(se, i0, i1, x);
ix = se & 0x7fff;
s = fabsl(x) - 1;
P = pa[0] + s * (pa[1] + s * (pa[2] +
s * (pa[3] + s * (pa[4] + s * (pa[5] + s * (pa[6] + s * pa[7]))))));
Q = qa[0] + s * (qa[1] + s * (qa[2] +
s * (qa[3] + s * (qa[4] + s * (qa[5] + s * (qa[6] + s))))));
return 1 - erx - P / Q;
}
if (ix >= 0x7fff) { /* erf(nan)=nan */
i = ((se & 0xffff) >> 15) << 1;
return (long double)(1 - i) + 1.0 / x; /* erf(+-inf)=+-1 */
}
static long double erfc2(uint32_t ix, long double x)
{
long double s,z,R,S;
uint32_t i0,i1;
ix = (ix << 16) | (i0 >> 16);
if (ix < 0x3ffed800) { /* |x| < 0.84375 */
if (ix < 0x3fde8000) { /* |x| < 2**-33 */
if (ix < 0x00080000)
return 0.125 * (8.0 * x + efx8 * x); /* avoid underflow */
return x + efx * x;
}
z = x * x;
r = pp[0] + z * (pp[1] +
z * (pp[2] + z * (pp[3] + z * (pp[4] + z * pp[5]))));
s = qq[0] + z * (qq[1] +
z * (qq[2] + z * (qq[3] + z * (qq[4] + z * (qq[5] + z)))));
y = r / s;
return x + x * y;
}
if (ix < 0x3fffa000) { /* 0.84375 <= |x| < 1.25 */
s = fabsl(x) - 1.0;
P = pa[0] + s * (pa[1] + s * (pa[2] +
s * (pa[3] + s * (pa[4] + s * (pa[5] + s * (pa[6] + s * pa[7]))))));
Q = qa[0] + s * (qa[1] + s * (qa[2] +
s * (qa[3] + s * (qa[4] + s * (qa[5] + s * (qa[6] + s))))));
if ((se & 0x8000) == 0)
return erx + P / Q;
return -erx - P / Q;
}
if (ix >= 0x4001d555) { /* inf > |x| >= 6.6666259765625 */
if ((se & 0x8000) == 0)
return 1.0 - tiny;
return tiny - 1.0;
if (ix < 0x3fffa000) /* 0.84375 <= |x| < 1.25 */
return erfc1(x);
x = fabsl(x);
s = 1 / (x * x);
if (ix < 0x4000b6db) { /* 1.25 <= |x| < 2.857 ~ 1/.35 */
R = ra[0] + s * (ra[1] + s * (ra[2] + s * (ra[3] + s * (ra[4] +
s * (ra[5] + s * (ra[6] + s * (ra[7] + s * ra[8])))))));
S = sa[0] + s * (sa[1] + s * (sa[2] + s * (sa[3] + s * (sa[4] +
s * (sa[5] + s * (sa[6] + s * (sa[7] + s * (sa[8] + s))))))));
} else { /* 2.857 <= |x| */
R = rb[0] + s * (rb[1] + s * (rb[2] + s * (rb[3] + s * (rb[4] +
s * (rb[5] + s * (rb[6] + s * rb[7]))))));
S = sb[0] + s * (sb[1] + s * (sb[2] + s * (sb[3] + s * (sb[4] +
s * (sb[5] + s * (sb[6] + s))))));
}
x = fabsl (x);
s = 1.0 / (x * x);
if (ix < 0x4000b6db) { /* 1.25 <= |x| < 2.85711669921875 ~ 1/.35 */
R = ra[0] + s * (ra[1] + s * (ra[2] + s * (ra[3] + s * (ra[4] +
s * (ra[5] + s * (ra[6] + s * (ra[7] + s * ra[8])))))));
S = sa[0] + s * (sa[1] + s * (sa[2] + s * (sa[3] + s * (sa[4] +
s * (sa[5] + s * (sa[6] + s * (sa[7] + s * (sa[8] + s))))))));
} else { /* 2.857 <= |x| < 6.667 */
} else if (ix < 0x4001d555) { /* 6.6666259765625 > |x| >= 1/.35 ~ 2.857143 */
R = rb[0] + s * (rb[1] + s * (rb[2] + s * (rb[3] + s * (rb[4] +
s * (rb[5] + s * (rb[6] + s * rb[7]))))));
S = sb[0] + s * (sb[1] + s * (sb[2] + s * (sb[3] + s * (sb[4] +
s * (sb[5] + s * (sb[6] + s))))));
} else { /* 107 > |x| >= 6.666 */
R = rc[0] + s * (rc[1] + s * (rc[2] + s * (rc[3] +
s * (rc[4] + s * rc[5]))));
S = sc[0] + s * (sc[1] + s * (sc[2] + s * (sc[3] +
s * (sc[4] + s))));
}
z = x;
GET_LDOUBLE_WORDS(i, i0, i1, z);
GET_LDOUBLE_WORDS(ix, i0, i1, z);
i1 = 0;
SET_LDOUBLE_WORDS(z, i, i0, i1);
r = expl(-z * z - 0.5625) * expl((z - x) * (z + x) + R / S);
if ((se & 0x8000) == 0)
return 1.0 - r / x;
return r / x - 1.0;
i0 &= 0xffffff00;
SET_LDOUBLE_WORDS(z, ix, i0, i1);
return expl(-z*z - 0.5625) * expl((z - x) * (z + x) + R / S) / x;
}
long double erfcl(long double x)
long double erfl(long double x)
{
int32_t hx, ix;
long double R, S, P, Q, s, y, z, r;
uint32_t se, i0, i1;
long double r, s, z, y;
uint32_t i0, i1, ix;
int sign;
GET_LDOUBLE_WORDS(se, i0, i1, x);
ix = se & 0x7fff;
if (ix >= 0x7fff) { /* erfc(nan) = nan, erfc(+-inf) = 0,2 */
return (long double)(((se & 0xffff) >> 15) << 1) + 1.0 / x;
GET_LDOUBLE_WORDS(ix, i0, i1, x);
sign = ix >> 15;
ix &= 0x7fff;
if (ix == 0x7fff) {
/* erf(nan)=nan, erf(+-inf)=+-1 */
return 1 - 2*sign + 1/x;
}
ix = (ix << 16) | (i0 >> 16);
if (ix < 0x3ffed800) { /* |x| < 0.84375 */
if (ix < 0x3fde8000) { /* |x| < 2**-33 */
return 0.125 * (8 * x + efx8 * x); /* avoid underflow */
}
z = x * x;
r = pp[0] + z * (pp[1] +
z * (pp[2] + z * (pp[3] + z * (pp[4] + z * pp[5]))));
s = qq[0] + z * (qq[1] +
z * (qq[2] + z * (qq[3] + z * (qq[4] + z * (qq[5] + z)))));
y = r / s;
return x + x * y;
}
if (ix < 0x4001d555) /* |x| < 6.6666259765625 */
y = 1 - erfc2(ix,x);
else
y = 1 - 0x1p-16382L;
return sign ? -y : y;
}
long double erfcl(long double x)
{
long double r, s, z, y;
uint32_t i0, i1, ix;
int sign;
GET_LDOUBLE_WORDS(ix, i0, i1, x);
sign = ix>>15;
ix &= 0x7fff;
if (ix == 0x7fff)
/* erfc(nan) = nan, erfc(+-inf) = 0,2 */
return 2*sign + 1/x;
ix = (ix << 16) | (i0 >> 16);
if (ix < 0x3ffed800) { /* |x| < 0.84375 */
@ -330,57 +350,10 @@ long double erfcl(long double x)
y = r / s;
if (ix < 0x3ffd8000) /* x < 1/4 */
return 1.0 - (x + x * y);
r = x * y;
r += x - 0.5L;
return 0.5L - r;
}
if (ix < 0x3fffa000) { /* 0.84375 <= |x| < 1.25 */
s = fabsl(x) - 1.0;
P = pa[0] + s * (pa[1] + s * (pa[2] +
s * (pa[3] + s * (pa[4] + s * (pa[5] + s * (pa[6] + s * pa[7]))))));
Q = qa[0] + s * (qa[1] + s * (qa[2] +
s * (qa[3] + s * (qa[4] + s * (qa[5] + s * (qa[6] + s))))));
if ((se & 0x8000) == 0) {
z = 1.0 - erx;
return z - P / Q;
}
z = erx + P / Q;
return 1.0 + z;
return 0.5 - (x - 0.5 + x * y);
}
if (ix < 0x4005d600) { /* |x| < 107 */
x = fabsl(x);
s = 1.0 / (x * x);
if (ix < 0x4000b6db) { /* 1.25 <= |x| < 2.85711669921875 ~ 1/.35 */
R = ra[0] + s * (ra[1] + s * (ra[2] + s * (ra[3] + s * (ra[4] +
s * (ra[5] + s * (ra[6] + s * (ra[7] + s * ra[8])))))));
S = sa[0] + s * (sa[1] + s * (sa[2] + s * (sa[3] + s * (sa[4] +
s * (sa[5] + s * (sa[6] + s * (sa[7] + s * (sa[8] + s))))))));
} else if (ix < 0x4001d555) { /* 6.6666259765625 > |x| >= 1/.35 ~ 2.857143 */
R = rb[0] + s * (rb[1] + s * (rb[2] + s * (rb[3] + s * (rb[4] +
s * (rb[5] + s * (rb[6] + s * rb[7]))))));
S = sb[0] + s * (sb[1] + s * (sb[2] + s * (sb[3] + s * (sb[4] +
s * (sb[5] + s * (sb[6] + s))))));
} else { /* 107 > |x| >= 6.666 */
if (se & 0x8000)
return 2.0 - tiny;/* x < -6.666 */
R = rc[0] + s * (rc[1] + s * (rc[2] + s * (rc[3] +
s * (rc[4] + s * rc[5]))));
S = sc[0] + s * (sc[1] + s * (sc[2] + s * (sc[3] +
s * (sc[4] + s))));
}
z = x;
GET_LDOUBLE_WORDS(hx, i0, i1, z);
i1 = 0;
i0 &= 0xffffff00;
SET_LDOUBLE_WORDS(z, hx, i0, i1);
r = expl(-z * z - 0.5625) * expl((z - x) * (z + x) + R / S);
if ((se & 0x8000) == 0)
return r / x;
return 2.0 - r / x;
}
if ((se & 0x8000) == 0)
return tiny * tiny;
return 2.0 - tiny;
if (ix < 0x4005d600) /* |x| < 107 */
return sign ? 2 - erfc2(ix,x) : erfc2(ix,x);
return sign ? 2 - 0x1p-16382L : 0x1p-16382L*0x1p-16382L;
}
#endif

Write Preview
Loading…
Cancel
Save