Print this page
11210 libm should be cstyle(1ONBLD) clean
*** 20,38 ****
*/
/*
* Copyright 2011 Nexenta Systems, Inc. All rights reserved.
*/
/*
* Copyright 2006 Sun Microsystems, Inc. All rights reserved.
* Use is subject to license terms.
*/
#pragma weak __erf = erf
#pragma weak __erfc = erfc
! /* INDENT OFF */
/*
* double erf(double x)
* double erfc(double x)
* x
* 2 |\
--- 20,39 ----
*/
/*
* Copyright 2011 Nexenta Systems, Inc. All rights reserved.
*/
+
/*
* Copyright 2006 Sun Microsystems, Inc. All rights reserved.
* Use is subject to license terms.
*/
#pragma weak __erf = erf
#pragma weak __erfc = erfc
!
/*
* double erf(double x)
* double erfc(double x)
* x
* 2 |\
*** 123,143 ****
* 7. Special case:
* erf(0) = 0, erf(inf) = 1, erf(-inf) = -1,
* erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2,
* erfc/erf(NaN) is NaN
*/
- /* INDENT ON */
#include "libm_macros.h"
#include <math.h>
static const double xxx[] = {
! /* tiny */ 1e-300,
/* half */ 5.00000000000000000000e-01, /* 3FE00000, 00000000 */
/* one */ 1.00000000000000000000e+00, /* 3FF00000, 00000000 */
/* two */ 2.00000000000000000000e+00, /* 40000000, 00000000 */
/* erx */ 8.45062911510467529297e-01, /* 3FEB0AC1, 60000000 */
/*
* Coefficients for approximation to erf on [0,0.84375]
*/
/* efx */ 1.28379167095512586316e-01, /* 3FC06EBA, 8214DB69 */
/* efx8 */ 1.02703333676410069053e+00, /* 3FF06EBA, 8214DB69 */
--- 124,145 ----
* 7. Special case:
* erf(0) = 0, erf(inf) = 1, erf(-inf) = -1,
* erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2,
* erfc/erf(NaN) is NaN
*/
#include "libm_macros.h"
#include <math.h>
static const double xxx[] = {
! /* tiny */
! 1e-300,
/* half */ 5.00000000000000000000e-01, /* 3FE00000, 00000000 */
/* one */ 1.00000000000000000000e+00, /* 3FF00000, 00000000 */
/* two */ 2.00000000000000000000e+00, /* 40000000, 00000000 */
/* erx */ 8.45062911510467529297e-01, /* 3FEB0AC1, 60000000 */
+
/*
* Coefficients for approximation to erf on [0,0.84375]
*/
/* efx */ 1.28379167095512586316e-01, /* 3FC06EBA, 8214DB69 */
/* efx8 */ 1.02703333676410069053e+00, /* 3FF06EBA, 8214DB69 */
*** 149,158 ****
--- 151,161 ----
/* qq1 */ 3.97917223959155352819e-01, /* 3FD97779, CDDADC09 */
/* qq2 */ 6.50222499887672944485e-02, /* 3FB0A54C, 5536CEBA */
/* qq3 */ 5.08130628187576562776e-03, /* 3F74D022, C4D36B0F */
/* qq4 */ 1.32494738004321644526e-04, /* 3F215DC9, 221C1A10 */
/* qq5 */ -3.96022827877536812320e-06, /* BED09C43, 42A26120 */
+
/*
* Coefficients for approximation to erf in [0.84375,1.25]
*/
/* pa0 */ -2.36211856075265944077e-03, /* BF6359B8, BEF77538 */
/* pa1 */ 4.14856118683748331666e-01, /* 3FDA8D00, AD92B34D */
*** 165,174 ****
--- 168,178 ----
/* qa2 */ 5.40397917702171048937e-01, /* 3FE14AF0, 92EB6F33 */
/* qa3 */ 7.18286544141962662868e-02, /* 3FB2635C, D99FE9A7 */
/* qa4 */ 1.26171219808761642112e-01, /* 3FC02660, E763351F */
/* qa5 */ 1.36370839120290507362e-02, /* 3F8BEDC2, 6B51DD1C */
/* qa6 */ 1.19844998467991074170e-02, /* 3F888B54, 5735151D */
+
/*
* Coefficients for approximation to erfc in [1.25,1/0.35]
*/
/* ra0 */ -9.86494403484714822705e-03, /* BF843412, 600D6435 */
/* ra1 */ -6.93858572707181764372e-01, /* BFE63416, E4BA7360 */
*** 184,193 ****
--- 188,198 ----
/* sa4 */ 6.45387271733267880336e+02, /* 40842B19, 21EC2868 */
/* sa5 */ 4.29008140027567833386e+02, /* 407AD021, 57700314 */
/* sa6 */ 1.08635005541779435134e+02, /* 405B28A3, EE48AE2C */
/* sa7 */ 6.57024977031928170135e+00, /* 401A47EF, 8E484A93 */
/* sa8 */ -6.04244152148580987438e-02, /* BFAEEFF2, EE749A62 */
+
/*
* Coefficients for approximation to erfc in [1/.35,28]
*/
/* rb0 */ -9.86494292470009928597e-03, /* BF843412, 39E86F4A */
/* rb1 */ -7.99283237680523006574e-01, /* BFE993BA, 70C285DE */
*** 208,217 ****
--- 213,223 ----
#define tiny xxx[0]
#define half xxx[1]
#define one xxx[2]
#define two xxx[3]
#define erx xxx[4]
+
/*
* Coefficients for approximation to erf on [0,0.84375]
*/
#define efx xxx[5]
#define efx8 xxx[6]
*** 223,232 ****
--- 229,239 ----
#define qq1 xxx[12]
#define qq2 xxx[13]
#define qq3 xxx[14]
#define qq4 xxx[15]
#define qq5 xxx[16]
+
/*
* Coefficients for approximation to erf in [0.84375,1.25]
*/
#define pa0 xxx[17]
#define pa1 xxx[18]
*** 239,248 ****
--- 246,256 ----
#define qa2 xxx[25]
#define qa3 xxx[26]
#define qa4 xxx[27]
#define qa5 xxx[28]
#define qa6 xxx[29]
+
/*
* Coefficients for approximation to erfc in [1.25,1/0.35]
*/
#define ra0 xxx[30]
#define ra1 xxx[31]
*** 258,267 ****
--- 266,276 ----
#define sa4 xxx[41]
#define sa5 xxx[42]
#define sa6 xxx[43]
#define sa7 xxx[44]
#define sa8 xxx[45]
+
/*
* Coefficients for approximation to erfc in [1/.35,28]
*/
#define rb0 xxx[46]
#define rb1 xxx[47]
*** 277,402 ****
#define sb5 xxx[57]
#define sb6 xxx[58]
#define sb7 xxx[59]
double
! erf(double x) {
int hx, ix, i;
double R, S, P, Q, s, y, z, r;
! hx = ((int *) &x)[HIWORD];
ix = hx & 0x7fffffff;
if (ix >= 0x7ff00000) { /* erf(nan)=nan */
#if defined(FPADD_TRAPS_INCOMPLETE_ON_NAN)
if (ix >= 0x7ff80000) /* assumes sparc-like QNaN */
return (x);
#endif
! i = ((unsigned) hx >> 31) << 1;
! return ((double) (1 - i) + one / x); /* erf(+-inf)=+-1 */
}
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 = one +
! 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) - one;
! P = pa0 + s * (pa1 + s * (pa2 + s * (pa3 + s * (pa4 +
! s * (pa5 + s * pa6)))));
! Q = one + s * (qa1 + s * (qa2 + s * (qa3 + s * (qa4 +
! s * (qa5 + s * qa6)))));
if (hx >= 0)
return (erx + P / Q);
else
return (-erx - P / Q);
}
if (ix >= 0x40180000) { /* inf > |x| >= 6 */
if (hx >= 0)
return (one - tiny);
else
return (tiny - one);
}
x = fabs(x);
s = one / (x * x);
if (ix < 0x4006DB6E) { /* |x| < 1/0.35 */
! R = ra0 + s * (ra1 + s * (ra2 + s * (ra3 + s * (ra4 +
! s * (ra5 + s * (ra6 + s * ra7))))));
! S = one + s * (sa1 + s * (sa2 + s * (sa3 + s * (sa4 +
! s * (sa5 + s * (sa6 + s * (sa7 + s * sa8)))))));
} else { /* |x| >= 1/0.35 */
! R = rb0 + s * (rb1 + s * (rb2 + s * (rb3 + s * (rb4 +
! s * (rb5 + s * rb6)))));
! S = one + s * (sb1 + s * (sb2 + s * (sb3 + s * (sb4 +
! s * (sb5 + s * (sb6 + s * sb7))))));
}
z = x;
! ((int *) &z)[LOWORD] = 0;
r = exp(-z * z - 0.5625) * exp((z - x) * (z + x) + R / S);
if (hx >= 0)
return (one - r / x);
else
return (r / x - one);
}
double
! erfc(double x) {
int hx, ix;
double R, S, P, Q, s, y, z, r;
! hx = ((int *) &x)[HIWORD];
ix = hx & 0x7fffffff;
if (ix >= 0x7ff00000) { /* erfc(nan)=nan */
#if defined(FPADD_TRAPS_INCOMPLETE_ON_NAN)
if (ix >= 0x7ff80000) /* assumes sparc-like QNaN */
return (x);
#endif
/* erfc(+-inf)=0,2 */
! return ((double) (((unsigned) hx >> 31) << 1) + one / x);
}
if (ix < 0x3feb0000) { /* |x| < 0.84375 */
if (ix < 0x3c700000) /* |x| < 2**-56 */
return (one - x);
z = x * x;
r = pp0 + z * (pp1 + z * (pp2 + z * (pp3 + z * pp4)));
! s = one +
! z * (qq1 + z * (qq2 + z * (qq3 + z * (qq4 + z * qq5))));
y = r / s;
if (hx < 0x3fd00000) { /* x < 1/4 */
return (one - (x + x * y));
} else {
r = x * y;
r += (x - half);
return (half - r);
}
}
if (ix < 0x3ff40000) { /* 0.84375 <= |x| < 1.25 */
s = fabs(x) - one;
! P = pa0 + s * (pa1 + s * (pa2 + s * (pa3 + s * (pa4 +
! s * (pa5 + s * pa6)))));
! Q = one + s * (qa1 + s * (qa2 + s * (qa3 + s * (qa4 +
! s * (qa5 + s * qa6)))));
if (hx >= 0) {
z = one - erx;
return (z - P / Q);
} else {
z = erx + P / Q;
return (one + z);
}
}
if (ix < 0x403c0000) { /* |x|<28 */
x = fabs(x);
s = one / (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 = one + s * (sa1 + s * (sa2 + s * (sa3 + s * (sa4 +
s * (sa5 + s * (sa6 + s * (sa7 + s * sa8)))))));
--- 286,430 ----
#define sb5 xxx[57]
#define sb6 xxx[58]
#define sb7 xxx[59]
double
! erf(double x)
! {
int hx, ix, i;
double R, S, P, Q, s, y, z, r;
! hx = ((int *)&x)[HIWORD];
ix = hx & 0x7fffffff;
+
if (ix >= 0x7ff00000) { /* erf(nan)=nan */
#if defined(FPADD_TRAPS_INCOMPLETE_ON_NAN)
if (ix >= 0x7ff80000) /* assumes sparc-like QNaN */
return (x);
#endif
! i = ((unsigned)hx >> 31) << 1;
! return ((double)(1 - i) + one / x); /* erf(+-inf)=+-1 */
}
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 = one + 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) - one;
! P = pa0 + s * (pa1 + s * (pa2 + s * (pa3 + s * (pa4 + s * (pa5 +
! s * pa6)))));
! Q = one + s * (qa1 + s * (qa2 + s * (qa3 + s * (qa4 + s * (qa5 +
! s * qa6)))));
!
if (hx >= 0)
return (erx + P / Q);
else
return (-erx - P / Q);
}
+
if (ix >= 0x40180000) { /* inf > |x| >= 6 */
if (hx >= 0)
return (one - tiny);
else
return (tiny - one);
}
+
x = fabs(x);
s = one / (x * x);
+
if (ix < 0x4006DB6E) { /* |x| < 1/0.35 */
! R = ra0 + s * (ra1 + s * (ra2 + s * (ra3 + s * (ra4 + s * (ra5 +
! s * (ra6 + s * ra7))))));
! S = one + s * (sa1 + s * (sa2 + s * (sa3 + s * (sa4 + s * (sa5 +
! s * (sa6 + s * (sa7 + s * sa8)))))));
} else { /* |x| >= 1/0.35 */
! R = rb0 + s * (rb1 + s * (rb2 + s * (rb3 + s * (rb4 + s * (rb5 +
! s * rb6)))));
! S = one + s * (sb1 + s * (sb2 + s * (sb3 + s * (sb4 + s * (sb5 +
! s * (sb6 + s * sb7))))));
}
+
z = x;
! ((int *)&z)[LOWORD] = 0;
r = exp(-z * z - 0.5625) * exp((z - x) * (z + x) + R / S);
+
if (hx >= 0)
return (one - r / x);
else
return (r / x - one);
}
double
! erfc(double x)
! {
int hx, ix;
double R, S, P, Q, s, y, z, r;
! hx = ((int *)&x)[HIWORD];
ix = hx & 0x7fffffff;
+
if (ix >= 0x7ff00000) { /* erfc(nan)=nan */
#if defined(FPADD_TRAPS_INCOMPLETE_ON_NAN)
if (ix >= 0x7ff80000) /* assumes sparc-like QNaN */
return (x);
#endif
/* erfc(+-inf)=0,2 */
! return ((double)(((unsigned)hx >> 31) << 1) + one / x);
}
if (ix < 0x3feb0000) { /* |x| < 0.84375 */
if (ix < 0x3c700000) /* |x| < 2**-56 */
return (one - x);
+
z = x * x;
r = pp0 + z * (pp1 + z * (pp2 + z * (pp3 + z * pp4)));
! s = one + z * (qq1 + z * (qq2 + z * (qq3 + z * (qq4 + z *
! qq5))));
y = r / s;
+
if (hx < 0x3fd00000) { /* x < 1/4 */
return (one - (x + x * y));
} else {
r = x * y;
r += (x - half);
return (half - r);
}
}
+
if (ix < 0x3ff40000) { /* 0.84375 <= |x| < 1.25 */
s = fabs(x) - one;
! P = pa0 + s * (pa1 + s * (pa2 + s * (pa3 + s * (pa4 + s * (pa5 +
! s * pa6)))));
! Q = one + s * (qa1 + s * (qa2 + s * (qa3 + s * (qa4 + s * (qa5 +
! s * qa6)))));
!
if (hx >= 0) {
z = one - erx;
return (z - P / Q);
} else {
z = erx + P / Q;
return (one + z);
}
}
+
if (ix < 0x403c0000) { /* |x|<28 */
x = fabs(x);
s = one / (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 = one + s * (sa1 + s * (sa2 + s * (sa3 + s * (sa4 +
s * (sa5 + s * (sa6 + s * (sa7 + s * sa8)))))));
*** 408,420 ****
R = rb0 + s * (rb1 + s * (rb2 + s * (rb3 + s * (rb4 +
s * (rb5 + s * rb6)))));
S = one + s * (sb1 + s * (sb2 + s * (sb3 + s * (sb4 +
s * (sb5 + s * (sb6 + s * sb7))))));
}
z = x;
! ((int *) &z)[LOWORD] = 0;
r = exp(-z * z - 0.5625) * exp((z - x) * (z + x) + R / S);
if (hx > 0)
return (r / x);
else
return (two - r / x);
} else {
--- 436,450 ----
R = rb0 + s * (rb1 + s * (rb2 + s * (rb3 + s * (rb4 +
s * (rb5 + s * rb6)))));
S = one + s * (sb1 + s * (sb2 + s * (sb3 + s * (sb4 +
s * (sb5 + s * (sb6 + s * sb7))))));
}
+
z = x;
! ((int *)&z)[LOWORD] = 0;
r = exp(-z * z - 0.5625) * exp((z - x) * (z + x) + R / S);
+
if (hx > 0)
return (r / x);
else
return (two - r / x);
} else {