code cleanup of named constants

zero, one, two, half are replaced by const literals
The policy was to use the f suffix for float consts (1.0f),
but don't use suffix for long double consts (these consts
can be exactly represented as double).
This commit is contained in:
nsz
2012-03-19 23:41:19 +01:00
parent b03255af77
commit 0cbb654791
73 changed files with 513 additions and 623 deletions

View File

@ -51,7 +51,6 @@
#include "libm.h" #include "libm.h"
static const double static const double
one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
C1 = 4.16666666666666019037e-02, /* 0x3FA55555, 0x5555554C */ C1 = 4.16666666666666019037e-02, /* 0x3FA55555, 0x5555554C */
C2 = -1.38888888888741095749e-03, /* 0xBF56C16C, 0x16C15177 */ C2 = -1.38888888888741095749e-03, /* 0xBF56C16C, 0x16C15177 */
C3 = 2.48015872894767294178e-05, /* 0x3EFA01A0, 0x19CB1590 */ C3 = 2.48015872894767294178e-05, /* 0x3EFA01A0, 0x19CB1590 */
@ -67,6 +66,6 @@ double __cos(double x, double y)
w = z*z; w = z*z;
r = z*(C1+z*(C2+z*C3)) + w*w*(C4+z*(C5+z*C6)); r = z*(C1+z*(C2+z*C3)) + w*w*(C4+z*(C5+z*C6));
hz = 0.5*z; hz = 0.5*z;
w = one-hz; w = 1.0-hz;
return w + (((one-w)-hz) + (z*r-x*y)); return w + (((1.0-w)-hz) + (z*r-x*y));
} }

View File

@ -18,7 +18,6 @@
/* |cos(x) - c(x)| < 2**-34.1 (~[-5.37e-11, 5.295e-11]). */ /* |cos(x) - c(x)| < 2**-34.1 (~[-5.37e-11, 5.295e-11]). */
static const double static const double
one = 1.0,
C0 = -0x1ffffffd0c5e81.0p-54, /* -0.499999997251031003120 */ C0 = -0x1ffffffd0c5e81.0p-54, /* -0.499999997251031003120 */
C1 = 0x155553e1053a42.0p-57, /* 0.0416666233237390631894 */ C1 = 0x155553e1053a42.0p-57, /* 0.0416666233237390631894 */
C2 = -0x16c087e80f1e27.0p-62, /* -0.00138867637746099294692 */ C2 = -0x16c087e80f1e27.0p-62, /* -0.00138867637746099294692 */
@ -32,5 +31,5 @@ float __cosdf(double x)
z = x*x; z = x*x;
w = z*z; w = z*z;
r = C2+z*C3; r = C2+z*C3;
return ((one+z*C0) + w*C1) + (w*z)*r; return ((1.0+z*C0) + w*C1) + (w*z)*r;
} }

View File

@ -41,8 +41,6 @@
* almost for free from the complications needed to search for the best * almost for free from the complications needed to search for the best
* higher coefficients. * higher coefficients.
*/ */
static const double one = 1.0;
static const long double static const long double
C1 = 0.0416666666666666666136L; /* 0xaaaaaaaaaaaaaa9b.0p-68 */ C1 = 0.0416666666666666666136L; /* 0xaaaaaaaaaaaaaa9b.0p-68 */
@ -61,7 +59,7 @@ long double __cosl(long double x, long double y)
z = x*x; z = x*x;
r = z*(C1+z*(C2+z*(C3+z*(C4+z*(C5+z*(C6+z*C7)))))); r = z*(C1+z*(C2+z*(C3+z*(C4+z*(C5+z*(C6+z*C7))))));
hz = 0.5*z; hz = 0.5*z;
w = one-hz; w = 1.0-hz;
return w + (((one-w)-hz) + (z*r-x*y)); return w + (((1.0-w)-hz) + (z*r-x*y));
} }
#endif #endif

View File

@ -29,7 +29,6 @@
* pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3) * pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3)
*/ */
static const double static const double
zero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */ two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */ invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
pio2_1 = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */ pio2_1 = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */
@ -163,7 +162,7 @@ medium:
} }
tx[2] = z; tx[2] = z;
nx = 3; nx = 3;
while (tx[nx-1] == zero) nx--; /* skip zero term */ while (tx[nx-1] == 0.0) nx--; /* skip zero term */
n = __rem_pio2_large(tx,ty,e0,nx,1); n = __rem_pio2_large(tx,ty,e0,nx,1);
if (hx < 0) { if (hx < 0) {
y[0] = -ty[0]; y[0] = -ty[0];

View File

@ -271,8 +271,6 @@ static const double PIo2[] = {
}; };
static const double static const double
zero = 0.0,
one = 1.0,
two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */ two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
twon24 = 5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */ twon24 = 5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */
@ -293,7 +291,7 @@ int __rem_pio2_large(double *x, double *y, int e0, int nx, int prec)
/* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */ /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */
j = jv-jx; m = jx+jk; j = jv-jx; m = jx+jk;
for (i=0; i<=m; i++,j++) for (i=0; i<=m; i++,j++)
f[i] = j<0 ? zero : (double)ipio2[j]; f[i] = j<0 ? 0.0 : (double)ipio2[j];
/* compute q[0],q[1],...q[jk] */ /* compute q[0],q[1],...q[jk] */
for (i=0; i<=jk; i++) { for (i=0; i<=jk; i++) {
@ -346,14 +344,14 @@ recompute:
} }
} }
if (ih == 2) { if (ih == 2) {
z = one - z; z = 1.0 - z;
if (carry != 0) if (carry != 0)
z -= scalbn(one,q0); z -= scalbn(1.0,q0);
} }
} }
/* check if recomputation is needed */ /* check if recomputation is needed */
if (z == zero) { if (z == 0.0) {
j = 0; j = 0;
for (i=jz-1; i>=jk; i--) j |= iq[i]; for (i=jz-1; i>=jk; i--) j |= iq[i];
if (j == 0) { /* need recomputation */ if (j == 0) { /* need recomputation */
@ -391,7 +389,7 @@ recompute:
} }
/* convert integer "bit" chunk to floating-point value */ /* convert integer "bit" chunk to floating-point value */
fw = scalbn(one,q0); fw = scalbn(1.0,q0);
for (i=jz; i>=0; i--) { for (i=jz; i>=0; i--) {
q[i] = fw*(double)iq[i]; q[i] = fw*(double)iq[i];
fw *= twon24; fw *= twon24;

View File

@ -32,7 +32,6 @@
* pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3) * pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3)
*/ */
static const double static const double
zero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */ two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
pio2_1 = 1.57079632679597125389e+00, /* 0x3FF921FB, 0x54444000 */ pio2_1 = 1.57079632679597125389e+00, /* 0x3FF921FB, 0x54444000 */
pio2_2 = -1.07463465549783099519e-12, /* -0x12e7b967674000.0p-92 */ pio2_2 = -1.07463465549783099519e-12, /* -0x12e7b967674000.0p-92 */
@ -114,7 +113,7 @@ int __rem_pio2l(long double x, long double *y)
} }
tx[2] = z; tx[2] = z;
nx = 3; nx = 3;
while (tx[nx-1] == zero) while (tx[nx-1] == 0.0)
nx--; /* skip zero term */ nx--; /* skip zero term */
n = __rem_pio2_large(tx,ty,e0,nx,2); n = __rem_pio2_large(tx,ty,e0,nx,2);
r = (long double)ty[0] + ty[1]; r = (long double)ty[0] + ty[1];

View File

@ -42,7 +42,6 @@
#include "libm.h" #include "libm.h"
static const double static const double
half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
S1 = -1.66666666666666324348e-01, /* 0xBFC55555, 0x55555549 */ S1 = -1.66666666666666324348e-01, /* 0xBFC55555, 0x55555549 */
S2 = 8.33333333332248946124e-03, /* 0x3F811111, 0x1110F8A6 */ S2 = 8.33333333332248946124e-03, /* 0x3F811111, 0x1110F8A6 */
S3 = -1.98412698298579493134e-04, /* 0xBF2A01A0, 0x19C161D5 */ S3 = -1.98412698298579493134e-04, /* 0xBF2A01A0, 0x19C161D5 */
@ -61,5 +60,5 @@ double __sin(double x, double y, int iy)
if (iy == 0) if (iy == 0)
return x + v*(S1 + z*r); return x + v*(S1 + z*r);
else else
return x - ((z*(half*y - v*r) - y) - v*S1); return x - ((z*(0.5*y - v*r) - y) - v*S1);
} }

View File

@ -24,8 +24,6 @@
* See __cosl.c for more details about the polynomial. * See __cosl.c for more details about the polynomial.
*/ */
static const double half = 0.5;
static const long double static const long double
S1 = -0.166666666666666666671L; /* -0xaaaaaaaaaaaaaaab.0p-66 */ S1 = -0.166666666666666666671L; /* -0xaaaaaaaaaaaaaaab.0p-66 */
@ -47,6 +45,6 @@ long double __sinl(long double x, long double y, int iy)
r = S2+z*(S3+z*(S4+z*(S5+z*(S6+z*(S7+z*S8))))); r = S2+z*(S3+z*(S4+z*(S5+z*(S6+z*(S7+z*S8)))));
if (iy == 0) if (iy == 0)
return x+v*(S1+z*r); return x+v*(S1+z*r);
return x-((z*(half*y-v*r)-y)-v*S1); return x-((z*(0.5*y-v*r)-y)-v*S1);
} }
#endif #endif

View File

@ -59,13 +59,9 @@ static const double T[] = {
7.14072491382608190305e-05, /* 3F12B80F, 32F0A7E9 */ 7.14072491382608190305e-05, /* 3F12B80F, 32F0A7E9 */
-1.85586374855275456654e-05, /* BEF375CB, DB605373 */ -1.85586374855275456654e-05, /* BEF375CB, DB605373 */
2.59073051863633712884e-05, /* 3EFB2A70, 74BF7AD4 */ 2.59073051863633712884e-05, /* 3EFB2A70, 74BF7AD4 */
/* one */ 1.00000000000000000000e+00, /* 3FF00000, 00000000 */ },
/* pio4 */ 7.85398163397448278999e-01, /* 3FE921FB, 54442D18 */ pio4 = 7.85398163397448278999e-01, /* 3FE921FB, 54442D18 */
/* pio4lo */ 3.06161699786838301793e-17 /* 3C81A626, 33145C07 */ pio4lo = 3.06161699786838301793e-17; /* 3C81A626, 33145C07 */
};
#define one T[13]
#define pio4 T[14]
#define pio4lo T[15]
double __tan(double x, double y, int iy) double __tan(double x, double y, int iy)
{ {

View File

@ -36,7 +36,6 @@
#include "libm.h" #include "libm.h"
static const double static const double
one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
pi = 3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */ pi = 3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */
pio2_hi = 1.57079632679489655800e+00; /* 0x3FF921FB, 0x54442D18 */ pio2_hi = 1.57079632679489655800e+00; /* 0x3FF921FB, 0x54442D18 */
static const volatile double static const volatile double
@ -75,25 +74,25 @@ double acos(double x)
return pio2_hi + pio2_lo; return pio2_hi + pio2_lo;
z = x*x; z = x*x;
p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5))))); p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4))); q = 1.0+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
r = p/q; r = p/q;
return pio2_hi - (x - (pio2_lo-x*r)); return pio2_hi - (x - (pio2_lo-x*r));
} else if (hx < 0) { /* x < -0.5 */ } else if (hx < 0) { /* x < -0.5 */
z = (one+x)*0.5; z = (1.0+x)*0.5;
p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5))))); p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4))); q = 1.0+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
s = sqrt(z); s = sqrt(z);
r = p/q; r = p/q;
w = r*s-pio2_lo; w = r*s-pio2_lo;
return pi - 2.0*(s+w); return pi - 2.0*(s+w);
} else { /* x > 0.5 */ } else { /* x > 0.5 */
z = (one-x)*0.5; z = (1.0-x)*0.5;
s = sqrt(z); s = sqrt(z);
df = s; df = s;
SET_LOW_WORD(df,0); SET_LOW_WORD(df,0);
c = (z-df*df)/(s+df); c = (z-df*df)/(s+df);
p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5))))); p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4))); q = 1.0+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
r = p/q; r = p/q;
w = r*s+c; w = r*s+c;
return 2.0*(df+w); return 2.0*(df+w);

View File

@ -16,7 +16,6 @@
#include "libm.h" #include "libm.h"
static const float static const float
one = 1.0000000000e+00, /* 0x3F800000 */
pi = 3.1415925026e+00, /* 0x40490fda */ pi = 3.1415925026e+00, /* 0x40490fda */
pio2_hi = 1.5707962513e+00; /* 0x3fc90fda */ pio2_hi = 1.5707962513e+00; /* 0x3fc90fda */
static const volatile float static const volatile float
@ -46,13 +45,13 @@ float acosf(float x)
return pio2_hi + pio2_lo; return pio2_hi + pio2_lo;
z = x*x; z = x*x;
p = z*(pS0+z*(pS1+z*pS2)); p = z*(pS0+z*(pS1+z*pS2));
q = one+z*qS1; q = 1.0f+z*qS1;
r = p/q; r = p/q;
return pio2_hi - (x - (pio2_lo-x*r)); return pio2_hi - (x - (pio2_lo-x*r));
} else if (hx < 0) { /* x < -0.5 */ } else if (hx < 0) { /* x < -0.5 */
z = (one+x)*0.5f; z = (1.0f+x)*0.5f;
p = z*(pS0+z*(pS1+z*pS2)); p = z*(pS0+z*(pS1+z*pS2));
q = one+z*qS1; q = 1.0f+z*qS1;
s = sqrtf(z); s = sqrtf(z);
r = p/q; r = p/q;
w = r*s-pio2_lo; w = r*s-pio2_lo;
@ -60,14 +59,14 @@ float acosf(float x)
} else { /* x > 0.5 */ } else { /* x > 0.5 */
int32_t idf; int32_t idf;
z = (one-x)*0.5f; z = (1.0f-x)*0.5f;
s = sqrtf(z); s = sqrtf(z);
df = s; df = s;
GET_FLOAT_WORD(idf,df); GET_FLOAT_WORD(idf,df);
SET_FLOAT_WORD(df,idf&0xfffff000); SET_FLOAT_WORD(df,idf&0xfffff000);
c = (z-df*df)/(s+df); c = (z-df*df)/(s+df);
p = z*(pS0+z*(pS1+z*pS2)); p = z*(pS0+z*(pS1+z*pS2));
q = one+z*qS1; q = 1.0f+z*qS1;
r = p/q; r = p/q;
w = r*s+c; w = r*s+c;
return 2.0f*(df+w); return 2.0f*(df+w);

View File

@ -27,7 +27,6 @@
#include "libm.h" #include "libm.h"
static const double static const double
one = 1.0,
ln2 = 6.93147180559945286227e-01; /* 0x3FE62E42, 0xFEFA39EF */ ln2 = 6.93147180559945286227e-01; /* 0x3FE62E42, 0xFEFA39EF */
double acosh(double x) double acosh(double x)
@ -47,9 +46,9 @@ double acosh(double x)
return 0.0; /* acosh(1) = 0 */ return 0.0; /* acosh(1) = 0 */
} else if (hx > 0x40000000) { /* 2**28 > x > 2 */ } else if (hx > 0x40000000) { /* 2**28 > x > 2 */
t = x*x; t = x*x;
return log(2.0*x - one/(x+sqrt(t-one))); return log(2.0*x - 1.0/(x+sqrt(t-1.0)));
} else { /* 1 < x < 2 */ } else { /* 1 < x < 2 */
t = x-one; t = x-1.0;
return log1p(t + sqrt(2.0*t+t*t)); return log1p(t + sqrt(2.0*t+t*t));
} }
} }

View File

@ -16,7 +16,6 @@
#include "libm.h" #include "libm.h"
static const float static const float
one = 1.0,
ln2 = 6.9314718246e-01; /* 0x3f317218 */ ln2 = 6.9314718246e-01; /* 0x3f317218 */
float acoshf(float x) float acoshf(float x)
@ -32,12 +31,12 @@ float acoshf(float x)
return x + x; return x + x;
return logf(x) + ln2; /* acosh(huge)=log(2x) */ return logf(x) + ln2; /* acosh(huge)=log(2x) */
} else if (hx == 0x3f800000) { } else if (hx == 0x3f800000) {
return 0.0; /* acosh(1) = 0 */ return 0.0f; /* acosh(1) = 0 */
} else if (hx > 0x40000000) { /* 2**28 > x > 2 */ } else if (hx > 0x40000000) { /* 2**28 > x > 2 */
t = x*x; t = x*x;
return logf(2.0f*x - one/(x+sqrtf(t-one))); return logf(2.0f*x - 1.0f/(x+sqrtf(t-1.0f)));
} else { /* 1 < x < 2 */ } else { /* 1 < x < 2 */
t = x-one; t = x-1.0f;
return log1pf(t + sqrtf(2.0f*t+t*t)); return log1pf(t + sqrtf(2.0f*t+t*t));
} }
} }

View File

@ -32,7 +32,6 @@ long double acoshl(long double x)
} }
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
static const long double static const long double
one = 1.0,
ln2 = 6.931471805599453094287e-01L; /* 0x3FFE, 0xB17217F7, 0xD1CF79AC */ ln2 = 6.931471805599453094287e-01L; /* 0x3FFE, 0xB17217F7, 0xD1CF79AC */
long double acoshl(long double x) long double acoshl(long double x)
@ -51,10 +50,10 @@ long double acoshl(long double x)
return 0.0; /* acosh(1) = 0 */ return 0.0; /* acosh(1) = 0 */
} else if (se > 0x4000) { /* x > 2 */ } else if (se > 0x4000) { /* x > 2 */
t = x*x; t = x*x;
return logl(2.0*x - one/(x + sqrtl(t - one))); return logl(2.0*x - 1.0/(x + sqrtl(t - 1.0)));
} }
/* 1 < x <= 2 */ /* 1 < x <= 2 */
t = x - one; t = x - 1.0;
return log1pl(t + sqrtl(2.0*t + t*t)); return log1pl(t + sqrtl(2.0*t + t*t));
} }
#endif #endif

View File

@ -25,7 +25,6 @@ long double acosl(long double x)
#include "__invtrigl.h" #include "__invtrigl.h"
static const long double static const long double
one = 1.00000000000000000000e+00,
pi = 3.14159265358979323846264338327950280e+00L; pi = 3.14159265358979323846264338327950280e+00L;
long double acosl(long double x) long double acosl(long double x)
@ -55,7 +54,7 @@ long double acosl(long double x)
r = p / q; r = p / q;
return pio2_hi - (x - (pio2_lo - x * r)); return pio2_hi - (x - (pio2_lo - x * r));
} else if (expsign < 0) { /* x < -0.5 */ } else if (expsign < 0) { /* x < -0.5 */
z = (one + x) * 0.5; z = (1.0 + x) * 0.5;
p = P(z); p = P(z);
q = Q(z); q = Q(z);
s = sqrtl(z); s = sqrtl(z);
@ -63,7 +62,7 @@ long double acosl(long double x)
w = r * s - pio2_lo; w = r * s - pio2_lo;
return pi - 2.0 * (s + w); return pi - 2.0 * (s + w);
} else { /* x > 0.5 */ } else { /* x > 0.5 */
z = (one - x) * 0.5; z = (1.0 - x) * 0.5;
s = sqrtl(z); s = sqrtl(z);
u.e = s; u.e = s;
u.bits.manl = 0; u.bits.manl = 0;

View File

@ -42,7 +42,6 @@
#include "libm.h" #include "libm.h"
static const double static const double
one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
huge = 1.000e+300, huge = 1.000e+300,
pio2_hi = 1.57079632679489655800e+00, /* 0x3FF921FB, 0x54442D18 */ pio2_hi = 1.57079632679489655800e+00, /* 0x3FF921FB, 0x54442D18 */
pio2_lo = 6.12323399573676603587e-17, /* 0x3C91A626, 0x33145C07 */ pio2_lo = 6.12323399573676603587e-17, /* 0x3C91A626, 0x33145C07 */
@ -76,20 +75,20 @@ double asin(double x)
return (x-x)/(x-x); /* asin(|x|>1) is NaN */ return (x-x)/(x-x); /* asin(|x|>1) is NaN */
} else if (ix < 0x3fe00000) { /* |x|<0.5 */ } else if (ix < 0x3fe00000) { /* |x|<0.5 */
if (ix < 0x3e500000) { /* if |x| < 2**-26 */ if (ix < 0x3e500000) { /* if |x| < 2**-26 */
if (huge+x > one) if (huge+x > 1.0)
return x; /* return x with inexact if x!=0*/ return x; /* return x with inexact if x!=0*/
} }
t = x*x; t = x*x;
p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5))))); p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5)))));
q = one+t*(qS1+t*(qS2+t*(qS3+t*qS4))); q = 1.0+t*(qS1+t*(qS2+t*(qS3+t*qS4)));
w = p/q; w = p/q;
return x + x*w; return x + x*w;
} }
/* 1 > |x| >= 0.5 */ /* 1 > |x| >= 0.5 */
w = one - fabs(x); w = 1.0 - fabs(x);
t = w*0.5; t = w*0.5;
p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5))))); p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5)))));
q = one+t*(qS1+t*(qS2+t*(qS3+t*qS4))); q = 1.0+t*(qS1+t*(qS2+t*(qS3+t*qS4)));
s = sqrt(t); s = sqrt(t);
if (ix >= 0x3FEF3333) { /* if |x| > 0.975 */ if (ix >= 0x3FEF3333) { /* if |x| > 0.975 */
w = p/q; w = p/q;

View File

@ -16,7 +16,6 @@
#include "libm.h" #include "libm.h"
static const float static const float
one = 1.0000000000e+00, /* 0x3F800000 */
huge = 1.000e+30, huge = 1.000e+30,
/* coefficients for R(x^2) */ /* coefficients for R(x^2) */
pS0 = 1.6666586697e-01, pS0 = 1.6666586697e-01,
@ -41,20 +40,20 @@ float asinf(float x)
return (x-x)/(x-x); /* asin(|x|>1) is NaN */ return (x-x)/(x-x); /* asin(|x|>1) is NaN */
} else if (ix < 0x3f000000) { /* |x|<0.5 */ } else if (ix < 0x3f000000) { /* |x|<0.5 */
if (ix < 0x39800000) { /* |x| < 2**-12 */ if (ix < 0x39800000) { /* |x| < 2**-12 */
if (huge+x > one) if (huge+x > 1.0f)
return x; /* return x with inexact if x!=0 */ return x; /* return x with inexact if x!=0 */
} }
t = x*x; t = x*x;
p = t*(pS0+t*(pS1+t*pS2)); p = t*(pS0+t*(pS1+t*pS2));
q = one+t*qS1; q = 1.0f+t*qS1;
w = p/q; w = p/q;
return x + x*w; return x + x*w;
} }
/* 1 > |x| >= 0.5 */ /* 1 > |x| >= 0.5 */
w = one - fabsf(x); w = 1.0f - fabsf(x);
t = w*0.5f; t = w*0.5f;
p = t*(pS0+t*(pS1+t*pS2)); p = t*(pS0+t*(pS1+t*pS2));
q = one+t*qS1; q = 1.0f+t*qS1;
s = sqrt(t); s = sqrt(t);
w = p/q; w = p/q;
t = pio2-2.0*(s+s*w); t = pio2-2.0*(s+s*w);

View File

@ -23,7 +23,6 @@
#include "libm.h" #include "libm.h"
static const double static const double
one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
ln2 = 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */ ln2 = 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */
huge= 1.00000000000000000000e+300; huge= 1.00000000000000000000e+300;
@ -38,17 +37,17 @@ double asinh(double x)
return x+x; return x+x;
if (ix < 0x3e300000) { /* |x| < 2**-28 */ if (ix < 0x3e300000) { /* |x| < 2**-28 */
/* return x inexact except 0 */ /* return x inexact except 0 */
if (huge+x > one) if (huge+x > 1.0)
return x; return x;
} }
if (ix > 0x41b00000) { /* |x| > 2**28 */ if (ix > 0x41b00000) { /* |x| > 2**28 */
w = log(fabs(x)) + ln2; w = log(fabs(x)) + ln2;
} else if (ix > 0x40000000) { /* 2**28 > |x| > 2.0 */ } else if (ix > 0x40000000) { /* 2**28 > |x| > 2.0 */
t = fabs(x); t = fabs(x);
w = log(2.0*t + one/(sqrt(x*x+one)+t)); w = log(2.0*t + 1.0/(sqrt(x*x+1.0)+t));
} else { /* 2.0 > |x| > 2**-28 */ } else { /* 2.0 > |x| > 2**-28 */
t = x*x; t = x*x;
w =log1p(fabs(x) + t/(one+sqrt(one+t))); w =log1p(fabs(x) + t/(1.0+sqrt(1.0+t)));
} }
if (hx > 0) if (hx > 0)
return w; return w;

View File

@ -16,7 +16,6 @@
#include "libm.h" #include "libm.h"
static const float static const float
one = 1.0000000000e+00, /* 0x3F800000 */
ln2 = 6.9314718246e-01, /* 0x3f317218 */ ln2 = 6.9314718246e-01, /* 0x3f317218 */
huge= 1.0000000000e+30; huge= 1.0000000000e+30;
@ -31,17 +30,17 @@ float asinhf(float x)
return x+x; return x+x;
if (ix < 0x31800000) { /* |x| < 2**-28 */ if (ix < 0x31800000) { /* |x| < 2**-28 */
/* return x inexact except 0 */ /* return x inexact except 0 */
if (huge+x > one) if (huge+x > 1.0f)
return x; return x;
} }
if (ix > 0x4d800000) { /* |x| > 2**28 */ if (ix > 0x4d800000) { /* |x| > 2**28 */
w = logf(fabsf(x)) + ln2; w = logf(fabsf(x)) + ln2;
} else if (ix > 0x40000000) { /* 2**28 > |x| > 2.0 */ } else if (ix > 0x40000000) { /* 2**28 > |x| > 2.0 */
t = fabsf(x); t = fabsf(x);
w = logf(2.0f*t + one/(sqrtf(x*x+one)+t)); w = logf(2.0f*t + 1.0f/(sqrtf(x*x+1.0f)+t));
} else { /* 2.0 > |x| > 2**-28 */ } else { /* 2.0 > |x| > 2**-28 */
t = x*x; t = x*x;
w =log1pf(fabsf(x) + t/(one+sqrtf(one+t))); w =log1pf(fabsf(x) + t/(1.0f+sqrtf(1.0f+t)));
} }
if (hx > 0) if (hx > 0)
return w; return w;

View File

@ -29,7 +29,6 @@ long double asinhl(long double x)
} }
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
static const long double static const long double
one = 1.000000000000000000000e+00L, /* 0x3FFF, 0x00000000, 0x00000000 */
ln2 = 6.931471805599453094287e-01L, /* 0x3FFE, 0xB17217F7, 0xD1CF79AC */ ln2 = 6.931471805599453094287e-01L, /* 0x3FFE, 0xB17217F7, 0xD1CF79AC */
huge = 1.000000000000000000e+4900L; huge = 1.000000000000000000e+4900L;
@ -44,17 +43,17 @@ long double asinhl(long double x)
return x + x; /* x is inf or NaN */ return x + x; /* x is inf or NaN */
if (ix < 0x3fde) { /* |x| < 2**-34 */ if (ix < 0x3fde) { /* |x| < 2**-34 */
/* return x, raise inexact if x != 0 */ /* return x, raise inexact if x != 0 */
if (huge+x > one) if (huge+x > 1.0)
return x; return x;
} }
if (ix > 0x4020) { /* |x| > 2**34 */ if (ix > 0x4020) { /* |x| > 2**34 */
w = logl(fabsl(x)) + ln2; w = logl(fabsl(x)) + ln2;
} else if (ix > 0x4000) { /* 2**34 > |x| > 2.0 */ } else if (ix > 0x4000) { /* 2**34 > |x| > 2.0 */
t = fabsl(x); t = fabsl(x);
w = logl(2.0*t + one/(sqrtl(x*x + one) + t)); w = logl(2.0*t + 1.0/(sqrtl(x*x + 1.0) + t));
} else { /* 2.0 > |x| > 2**-28 */ } else { /* 2.0 > |x| > 2**-28 */
t = x*x; t = x*x;
w =log1pl(fabsl(x) + t/(one + sqrtl(one + t))); w =log1pl(fabsl(x) + t/(1.0 + sqrtl(1.0 + t)));
} }
if (hx & 0x8000) if (hx & 0x8000)
return -w; return -w;

View File

@ -23,9 +23,7 @@ long double asinl(long double x)
} }
#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
#include "__invtrigl.h" #include "__invtrigl.h"
static const long double static const long double huge = 1.000e+300;
one = 1.00000000000000000000e+00,
huge = 1.000e+300;
long double asinl(long double x) long double asinl(long double x)
{ {
@ -45,7 +43,7 @@ long double asinl(long double x)
} else if (expt < BIAS-1) { /* |x|<0.5 */ } else if (expt < BIAS-1) { /* |x|<0.5 */
if (expt < ASIN_LINEAR) { /* if |x| is small, asinl(x)=x */ if (expt < ASIN_LINEAR) { /* if |x| is small, asinl(x)=x */
/* return x with inexact if x!=0 */ /* return x with inexact if x!=0 */
if (huge+x > one) if (huge+x > 1.0)
return x; return x;
} }
t = x*x; t = x*x;
@ -55,7 +53,7 @@ long double asinl(long double x)
return x + x*w; return x + x*w;
} }
/* 1 > |x| >= 0.5 */ /* 1 > |x| >= 0.5 */
w = one - fabsl(x); w = 1.0 - fabsl(x);
t = w*0.5; t = w*0.5;
p = P(t); p = P(t);
q = Q(t); q = Q(t);

View File

@ -60,9 +60,7 @@ static const double aT[] = {
1.62858201153657823623e-02, /* 0x3F90AD3A, 0xE322DA11 */ 1.62858201153657823623e-02, /* 0x3F90AD3A, 0xE322DA11 */
}; };
static const double static const double huge = 1.0e300;
one = 1.0,
huge = 1.0e300;
double atan(double x) double atan(double x)
{ {
@ -86,7 +84,7 @@ double atan(double x)
if (ix < 0x3fdc0000) { /* |x| < 0.4375 */ if (ix < 0x3fdc0000) { /* |x| < 0.4375 */
if (ix < 0x3e400000) { /* |x| < 2^-27 */ if (ix < 0x3e400000) { /* |x| < 2^-27 */
/* raise inexact */ /* raise inexact */
if (huge+x > one) if (huge+x > 1.0)
return x; return x;
} }
id = -1; id = -1;
@ -95,15 +93,15 @@ double atan(double x)
if (ix < 0x3ff30000) { /* |x| < 1.1875 */ if (ix < 0x3ff30000) { /* |x| < 1.1875 */
if (ix < 0x3fe60000) { /* 7/16 <= |x| < 11/16 */ if (ix < 0x3fe60000) { /* 7/16 <= |x| < 11/16 */
id = 0; id = 0;
x = (2.0*x-one)/(2.0+x); x = (2.0*x-1.0)/(2.0+x);
} else { /* 11/16 <= |x| < 19/16 */ } else { /* 11/16 <= |x| < 19/16 */
id = 1; id = 1;
x = (x-one)/(x+one); x = (x-1.0)/(x+1.0);
} }
} else { } else {
if (ix < 0x40038000) { /* |x| < 2.4375 */ if (ix < 0x40038000) { /* |x| < 2.4375 */
id = 2; id = 2;
x = (x-1.5)/(one+1.5*x); x = (x-1.5)/(1.0+1.5*x);
} else { /* 2.4375 <= |x| < 2^66 */ } else { /* 2.4375 <= |x| < 2^66 */
id = 3; id = 3;
x = -1.0/x; x = -1.0/x;

View File

@ -42,7 +42,6 @@
static const volatile double static const volatile double
tiny = 1.0e-300; tiny = 1.0e-300;
static const double static const double
zero = 0.0,
pi_o_4 = 7.8539816339744827900E-01, /* 0x3FE921FB, 0x54442D18 */ pi_o_4 = 7.8539816339744827900E-01, /* 0x3FE921FB, 0x54442D18 */
pi_o_2 = 1.5707963267948965580E+00, /* 0x3FF921FB, 0x54442D18 */ pi_o_2 = 1.5707963267948965580E+00, /* 0x3FF921FB, 0x54442D18 */
pi = 3.1415926535897931160E+00; /* 0x400921FB, 0x54442D18 */ pi = 3.1415926535897931160E+00; /* 0x400921FB, 0x54442D18 */
@ -89,8 +88,8 @@ double atan2(double y, double x)
} }
} else { } else {
switch(m) { switch(m) {
case 0: return zero; /* atan(+...,+INF) */ case 0: return 0.0; /* atan(+...,+INF) */
case 1: return -zero; /* atan(-...,+INF) */ case 1: return -0.0; /* atan(-...,+INF) */
case 2: return pi+tiny; /* atan(+...,-INF) */ case 2: return pi+tiny; /* atan(+...,-INF) */
case 3: return -pi-tiny; /* atan(-...,-INF) */ case 3: return -pi-tiny; /* atan(-...,-INF) */
} }

View File

@ -18,7 +18,6 @@
static const volatile float static const volatile float
tiny = 1.0e-30; tiny = 1.0e-30;
static const float static const float
zero = 0.0,
pi_o_4 = 7.8539818525e-01, /* 0x3f490fdb */ pi_o_4 = 7.8539818525e-01, /* 0x3f490fdb */
pi_o_2 = 1.5707963705e+00, /* 0x3fc90fdb */ pi_o_2 = 1.5707963705e+00, /* 0x3fc90fdb */
pi = 3.1415927410e+00; /* 0x40490fdb */ pi = 3.1415927410e+00; /* 0x40490fdb */
@ -63,8 +62,8 @@ float atan2f(float y, float x)
} }
} else { } else {
switch (m) { switch (m) {
case 0: return zero; /* atan(+...,+INF) */ case 0: return 0.0f; /* atan(+...,+INF) */
case 1: return -zero; /* atan(-...,+INF) */ case 1: return -0.0f; /* atan(-...,+INF) */
case 2: return pi+tiny; /* atan(+...,-INF) */ case 2: return pi+tiny; /* atan(+...,-INF) */
case 3: return -pi-tiny; /* atan(-...,-INF) */ case 3: return -pi-tiny; /* atan(-...,-INF) */
} }

View File

@ -27,7 +27,6 @@ long double atan2l(long double y, long double x)
static const volatile long double static const volatile long double
tiny = 1.0e-300; tiny = 1.0e-300;
static const long double static const long double
zero = 0.0,
pi = 3.14159265358979323846264338327950280e+00L; pi = 3.14159265358979323846264338327950280e+00L;
long double atan2l(long double y, long double x) long double atan2l(long double y, long double x)
@ -75,8 +74,8 @@ long double atan2l(long double y, long double x)
} }
} else { } else {
switch(m) { switch(m) {
case 0: return zero; /* atan(+...,+INF) */ case 0: return 0.0; /* atan(+...,+INF) */
case 1: return -zero; /* atan(-...,+INF) */ case 1: return -0.0; /* atan(-...,+INF) */
case 2: return pi+tiny; /* atan(+...,-INF) */ case 2: return pi+tiny; /* atan(+...,-INF) */
case 3: return -pi-tiny; /* atan(-...,-INF) */ case 3: return -pi-tiny; /* atan(-...,-INF) */
} }

View File

@ -38,9 +38,7 @@ static const float aT[] = {
6.1687607318e-02, 6.1687607318e-02,
}; };
static const float static const float huge = 1.0e30;
one = 1.0,
huge = 1.0e30;
float atanf(float x) float atanf(float x)
{ {
@ -60,7 +58,7 @@ float atanf(float x)
if (ix < 0x3ee00000) { /* |x| < 0.4375 */ if (ix < 0x3ee00000) { /* |x| < 0.4375 */
if (ix < 0x39800000) { /* |x| < 2**-12 */ if (ix < 0x39800000) { /* |x| < 2**-12 */
/* raise inexact */ /* raise inexact */
if(huge+x>one) if(huge+x>1.0f)
return x; return x;
} }
id = -1; id = -1;
@ -69,15 +67,15 @@ float atanf(float x)
if (ix < 0x3f980000) { /* |x| < 1.1875 */ if (ix < 0x3f980000) { /* |x| < 1.1875 */
if (ix < 0x3f300000) { /* 7/16 <= |x| < 11/16 */ if (ix < 0x3f300000) { /* 7/16 <= |x| < 11/16 */
id = 0; id = 0;
x = (2.0f*x - one)/(2.0f + x); x = (2.0f*x - 1.0f)/(2.0f + x);
} else { /* 11/16 <= |x| < 19/16 */ } else { /* 11/16 <= |x| < 19/16 */
id = 1; id = 1;
x = (x - one)/(x + one); x = (x - 1.0f)/(x + 1.0f);
} }
} else { } else {
if (ix < 0x401c0000) { /* |x| < 2.4375 */ if (ix < 0x401c0000) { /* |x| < 2.4375 */
id = 2; id = 2;
x = (x - 1.5f)/(one + 1.5f*x); x = (x - 1.5f)/(1.0f + 1.5f*x);
} else { /* 2.4375 <= |x| < 2**26 */ } else { /* 2.4375 <= |x| < 2**26 */
id = 3; id = 3;
x = -1.0f/x; x = -1.0f/x;

View File

@ -30,8 +30,7 @@
#include "libm.h" #include "libm.h"
static const double one = 1.0, huge = 1e300; static const double huge = 1e300;
static const double zero = 0.0;
double atanh(double x) double atanh(double x)
{ {
@ -44,15 +43,15 @@ double atanh(double x)
if ((ix | ((lx|-lx)>>31)) > 0x3ff00000) /* |x| > 1 */ if ((ix | ((lx|-lx)>>31)) > 0x3ff00000) /* |x| > 1 */
return (x-x)/(x-x); return (x-x)/(x-x);
if (ix == 0x3ff00000) if (ix == 0x3ff00000)
return x/zero; return x/0.0;
if (ix < 0x3e300000 && (huge+x) > zero) /* x < 2**-28 */ if (ix < 0x3e300000 && (huge+x) > 0.0) /* x < 2**-28 */
return x; return x;
SET_HIGH_WORD(x, ix); SET_HIGH_WORD(x, ix);
if (ix < 0x3fe00000) { /* x < 0.5 */ if (ix < 0x3fe00000) { /* x < 0.5 */
t = x+x; t = x+x;
t = 0.5*log1p(t + t*x/(one-x)); t = 0.5*log1p(t + t*x/(1.0-x));
} else } else
t = 0.5*log1p((x+x)/(one-x)); t = 0.5*log1p((x+x)/(1.0-x));
if (hx >= 0) if (hx >= 0)
return t; return t;
return -t; return -t;

View File

@ -15,8 +15,7 @@
#include "libm.h" #include "libm.h"
static const float one = 1.0, huge = 1e30; static const float huge = 1e30;
static const float zero = 0.0;
float atanhf(float x) float atanhf(float x)
{ {
@ -28,15 +27,15 @@ float atanhf(float x)
if (ix > 0x3f800000) /* |x| > 1 */ if (ix > 0x3f800000) /* |x| > 1 */
return (x-x)/(x-x); return (x-x)/(x-x);
if (ix == 0x3f800000) if (ix == 0x3f800000)
return x/zero; return x/0.0f;
if (ix < 0x31800000 && huge+x > zero) /* x < 2**-28 */ if (ix < 0x31800000 && huge+x > 0.0f) /* x < 2**-28 */
return x; return x;
SET_FLOAT_WORD(x, ix); SET_FLOAT_WORD(x, ix);
if (ix < 0x3f000000) { /* x < 0.5 */ if (ix < 0x3f000000) { /* x < 0.5 */
t = x+x; t = x+x;
t = 0.5f*log1pf(t + t*x/(one-x)); t = 0.5f*log1pf(t + t*x/(1.0f-x));
} else } else
t = 0.5f*log1pf((x+x)/(one-x)); t = 0.5f*log1pf((x+x)/(1.0f-x));
if (hx >= 0) if (hx >= 0)
return t; return t;
return -t; return -t;

View File

@ -34,7 +34,7 @@ long double atanhl(long double x)
return atanh(x); return atanh(x);
} }
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
static const long double zero = 0.0, one = 1.0, huge = 1e4900L; static const long double huge = 1e4900L;
long double atanhl(long double x) long double atanhl(long double x)
{ {
@ -48,15 +48,15 @@ long double atanhl(long double x)
/* |x| > 1 */ /* |x| > 1 */
return (x-x)/(x-x); return (x-x)/(x-x);
if (ix == 0x3fff) if (ix == 0x3fff)
return x/zero; return x/0.0;
if (ix < 0x3fe3 && huge+x > zero) /* x < 2**-28 */ if (ix < 0x3fe3 && huge+x > 0.0) /* x < 2**-28 */
return x; return x;
SET_LDOUBLE_EXP(x, ix); SET_LDOUBLE_EXP(x, ix);
if (ix < 0x3ffe) { /* x < 0.5 */ if (ix < 0x3ffe) { /* x < 0.5 */
t = x + x; t = x + x;
t = 0.5*log1pl(t + t*x/(one-x)); t = 0.5*log1pl(t + t*x/(1.0 - x));
} else } else
t = 0.5*log1pl((x + x)/(one - x)); t = 0.5*log1pl((x + x)/(1.0 - x));
if (se <= 0x7fff) if (se <= 0x7fff)
return t; return t;
return -t; return -t;

View File

@ -23,9 +23,7 @@ long double atanl(long double x)
} }
#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
#include "__invtrigl.h" #include "__invtrigl.h"
static const long double static const long double huge = 1.0e300;
one = 1.0,
huge = 1.0e300;
long double atanl(long double x) long double atanl(long double x)
{ {
@ -53,7 +51,7 @@ long double atanl(long double x)
if (expman < ((BIAS - 2) << 8) + 0xc0) { /* |x| < 0.4375 */ if (expman < ((BIAS - 2) << 8) + 0xc0) { /* |x| < 0.4375 */
if (expt < ATAN_LINEAR) { /* if |x| is small, atanl(x)~=x */ if (expt < ATAN_LINEAR) { /* if |x| is small, atanl(x)~=x */
/* raise inexact */ /* raise inexact */
if (huge+x > one) if (huge+x > 1.0)
return x; return x;
} }
id = -1; id = -1;
@ -62,15 +60,15 @@ long double atanl(long double x)
if (expman < (BIAS << 8) + 0x30) { /* |x| < 1.1875 */ if (expman < (BIAS << 8) + 0x30) { /* |x| < 1.1875 */
if (expman < ((BIAS - 1) << 8) + 0x60) { /* 7/16 <= |x| < 11/16 */ if (expman < ((BIAS - 1) << 8) + 0x60) { /* 7/16 <= |x| < 11/16 */
id = 0; id = 0;
x = (2.0*x-one)/(2.0+x); x = (2.0*x-1.0)/(2.0+x);
} else { /* 11/16 <= |x| < 19/16 */ } else { /* 11/16 <= |x| < 19/16 */
id = 1; id = 1;
x = (x-one)/(x+one); x = (x-1.0)/(x+1.0);
} }
} else { } else {
if (expman < ((BIAS + 1) << 8) + 0x38) { /* |x| < 2.4375 */ if (expman < ((BIAS + 1) << 8) + 0x38) { /* |x| < 2.4375 */
id = 2; id = 2;
x = (x-1.5)/(one+1.5*x); x = (x-1.5)/(1.0+1.5*x);
} else { /* 2.4375 <= |x| < 2^ATAN_CONST */ } else { /* 2.4375 <= |x| < 2^ATAN_CONST */
id = 3; id = 3;
x = -1.0/x; x = -1.0/x;

View File

@ -32,7 +32,7 @@
#include "libm.h" #include "libm.h"
static const double one = 1.0, half = 0.5, huge = 1.0e300; static const double huge = 1.0e300;
double cosh(double x) double cosh(double x)
{ {
@ -49,21 +49,21 @@ double cosh(double x)
/* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */ /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
if (ix < 0x3fd62e43) { if (ix < 0x3fd62e43) {
t = expm1(fabs(x)); t = expm1(fabs(x));
w = one+t; w = 1.0+t;
if (ix < 0x3c800000) if (ix < 0x3c800000)
return w; /* cosh(tiny) = 1 */ return w; /* cosh(tiny) = 1 */
return one + (t*t)/(w+w); return 1.0 + (t*t)/(w+w);
} }
/* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|))/2; */ /* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|))/2; */
if (ix < 0x40360000) { if (ix < 0x40360000) {
t = exp(fabs(x)); t = exp(fabs(x));
return half*t + half/t; return 0.5*t + 0.5/t;
} }
/* |x| in [22, log(maxdouble)] return half*exp(|x|) */ /* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */
if (ix < 0x40862E42) if (ix < 0x40862E42)
return half*exp(fabs(x)); return 0.5*exp(fabs(x));
/* |x| in [log(maxdouble), overflowthresold] */ /* |x| in [log(maxdouble), overflowthresold] */
if (ix <= 0x408633CE) if (ix <= 0x408633CE)

View File

@ -15,7 +15,7 @@
#include "libm.h" #include "libm.h"
static const float one = 1.0, half = 0.5, huge = 1.0e30; static const float huge = 1.0e30;
float coshf(float x) float coshf(float x)
{ {
@ -32,21 +32,21 @@ float coshf(float x)
/* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */ /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
if (ix < 0x3eb17218) { if (ix < 0x3eb17218) {
t = expm1f(fabsf(x)); t = expm1f(fabsf(x));
w = one+t; w = 1.0f+t;
if (ix<0x39800000) if (ix<0x39800000)
return one; /* cosh(tiny) = 1 */ return 1.0f; /* cosh(tiny) = 1 */
return one + (t*t)/(w+w); return 1.0f + (t*t)/(w+w);
} }
/* |x| in [0.5*ln2,9], return (exp(|x|)+1/exp(|x|))/2; */ /* |x| in [0.5*ln2,9], return (exp(|x|)+1/exp(|x|))/2; */
if (ix < 0x41100000) { if (ix < 0x41100000) {
t = expf(fabsf(x)); t = expf(fabsf(x));
return half*t + half/t; return 0.5f*t + 0.5f/t;
} }
/* |x| in [9, log(maxfloat)] return half*exp(|x|) */ /* |x| in [9, log(maxfloat)] return 0.5f*exp(|x|) */
if (ix < 0x42b17217) if (ix < 0x42b17217)
return half*expf(fabsf(x)); return 0.5f*expf(fabsf(x));
/* |x| in [log(maxfloat), overflowthresold] */ /* |x| in [log(maxfloat), overflowthresold] */
if (ix <= 0x42b2d4fc) if (ix <= 0x42b2d4fc)

View File

@ -38,7 +38,7 @@ long double coshl(long double x)
return cosh(x); return cosh(x);
} }
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
static const long double one = 1.0, half = 0.5, huge = 1.0e4900L; static const long double huge = 1.0e4900L;
long double coshl(long double x) long double coshl(long double x)
{ {
@ -56,27 +56,27 @@ long double coshl(long double x)
/* |x| in [0,0.5*ln2], return 1+expm1l(|x|)^2/(2*expl(|x|)) */ /* |x| in [0,0.5*ln2], return 1+expm1l(|x|)^2/(2*expl(|x|)) */
if (ex < 0x3ffd || (ex == 0x3ffd && mx < 0xb17217f7u)) { if (ex < 0x3ffd || (ex == 0x3ffd && mx < 0xb17217f7u)) {
t = expm1l(fabsl(x)); t = expm1l(fabsl(x));
w = one + t; w = 1.0 + t;
if (ex < 0x3fbc) return w; /* cosh(tiny) = 1 */ if (ex < 0x3fbc) return w; /* cosh(tiny) = 1 */
return one+(t*t)/(w+w); return 1.0+(t*t)/(w+w);
} }
/* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */ /* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */
if (ex < 0x4003 || (ex == 0x4003 && mx < 0xb0000000u)) { if (ex < 0x4003 || (ex == 0x4003 && mx < 0xb0000000u)) {
t = expl(fabsl(x)); t = expl(fabsl(x));
return half*t + half/t; return 0.5*t + 0.5/t;
} }
/* |x| in [22, ln(maxdouble)] return half*exp(|x|) */ /* |x| in [22, ln(maxdouble)] return 0.5*exp(|x|) */
if (ex < 0x400c || (ex == 0x400c && mx < 0xb1700000u)) if (ex < 0x400c || (ex == 0x400c && mx < 0xb1700000u))
return half*expl(fabsl(x)); return 0.5*expl(fabsl(x));
/* |x| in [log(maxdouble), log(2*maxdouble)) */ /* |x| in [log(maxdouble), log(2*maxdouble)) */
if (ex == 0x400c && (mx < 0xb174ddc0u || if (ex == 0x400c && (mx < 0xb174ddc0u ||
(mx == 0xb174ddc0u && lx < 0x31aec0ebu))) (mx == 0xb174ddc0u && lx < 0x31aec0ebu)))
{ {
w = expl(half*fabsl(x)); w = expl(0.5*fabsl(x));
t = half*w; t = 0.5*w;
return t*w; return t*w;
} }

View File

@ -107,9 +107,6 @@
static const double static const double
tiny = 1e-300, tiny = 1e-300,
half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
two = 2.00000000000000000000e+00, /* 0x40000000, 0x00000000 */
/* c = (float)0.84506291151 */ /* c = (float)0.84506291151 */
erx = 8.45062911510467529297e-01, /* 0x3FEB0AC1, 0x60000000 */ erx = 8.45062911510467529297e-01, /* 0x3FEB0AC1, 0x60000000 */
/* /*
@ -190,7 +187,7 @@ double erf(double x)
if (ix >= 0x7ff00000) { if (ix >= 0x7ff00000) {
/* erf(nan)=nan, erf(+-inf)=+-1 */ /* erf(nan)=nan, erf(+-inf)=+-1 */
i = ((uint32_t)hx>>31)<<1; i = ((uint32_t)hx>>31)<<1;
return (double)(1-i) + one/x; return (double)(1-i) + 1.0/x;
} }
if (ix < 0x3feb0000) { /* |x|<0.84375 */ if (ix < 0x3feb0000) { /* |x|<0.84375 */
if (ix < 0x3e300000) { /* |x|<2**-28 */ if (ix < 0x3e300000) { /* |x|<2**-28 */
@ -201,42 +198,42 @@ double erf(double x)
} }
z = x*x; z = x*x;
r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4))); r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5)))); s = 1.0+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
y = r/s; y = r/s;
return x + x*y; return x + x*y;
} }
if (ix < 0x3ff40000) { /* 0.84375 <= |x| < 1.25 */ if (ix < 0x3ff40000) { /* 0.84375 <= |x| < 1.25 */
s = fabs(x)-one; s = fabs(x)-1.0;
P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6))))); 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))))); Q = 1.0+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
if (hx >= 0) if (hx >= 0)
return erx + P/Q; return erx + P/Q;
return -erx - P/Q; return -erx - P/Q;
} }
if (ix >= 0x40180000) { /* inf > |x| >= 6 */ if (ix >= 0x40180000) { /* inf > |x| >= 6 */
if (hx >= 0) if (hx >= 0)
return one-tiny; return 1.0 - tiny;
return tiny-one; return tiny - 1.0;
} }
x = fabs(x); x = fabs(x);
s = one/(x*x); s = 1.0/(x*x);
if (ix < 0x4006DB6E) { /* |x| < 1/0.35 */ if (ix < 0x4006DB6E) { /* |x| < 1/0.35 */
R = ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*( R = ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
ra5+s*(ra6+s*ra7)))))); ra5+s*(ra6+s*ra7))))));
S = one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*( S = 1.0+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
sa5+s*(sa6+s*(sa7+s*sa8))))))); sa5+s*(sa6+s*(sa7+s*sa8)))))));
} else { /* |x| >= 1/0.35 */ } else { /* |x| >= 1/0.35 */
R = rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*( R = rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
rb5+s*rb6))))); rb5+s*rb6)))));
S = one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*( S = 1.0+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
sb5+s*(sb6+s*sb7)))))); sb5+s*(sb6+s*sb7))))));
} }
z = x; z = x;
SET_LOW_WORD(z,0); SET_LOW_WORD(z,0);
r = exp(-z*z-0.5625)*exp((z-x)*(z+x)+R/S); r = exp(-z*z-0.5625)*exp((z-x)*(z+x)+R/S);
if (hx >= 0) if (hx >= 0)
return one-r/x; return 1.0 - r/x;
return r/x-one; return r/x - 1.0;
} }
double erfc(double x) double erfc(double x)
@ -248,49 +245,49 @@ double erfc(double x)
ix = hx & 0x7fffffff; ix = hx & 0x7fffffff;
if (ix >= 0x7ff00000) { if (ix >= 0x7ff00000) {
/* erfc(nan)=nan, erfc(+-inf)=0,2 */ /* erfc(nan)=nan, erfc(+-inf)=0,2 */
return (double)(((uint32_t)hx>>31)<<1) + one/x; return (double)(((uint32_t)hx>>31)<<1) + 1.0/x;
} }
if (ix < 0x3feb0000) { /* |x| < 0.84375 */ if (ix < 0x3feb0000) { /* |x| < 0.84375 */
if (ix < 0x3c700000) /* |x| < 2**-56 */ if (ix < 0x3c700000) /* |x| < 2**-56 */
return one - x; return 1.0 - x;
z = x*x; z = x*x;
r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4))); r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5)))); s = 1.0+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
y = r/s; y = r/s;
if (hx < 0x3fd00000) { /* x < 1/4 */ if (hx < 0x3fd00000) { /* x < 1/4 */
return one - (x+x*y); return 1.0 - (x+x*y);
} else { } else {
r = x*y; r = x*y;
r += x-half; r += x - 0.5;
return half - r ; return 0.5 - r ;
} }
} }
if (ix < 0x3ff40000) { /* 0.84375 <= |x| < 1.25 */ if (ix < 0x3ff40000) { /* 0.84375 <= |x| < 1.25 */
s = fabs(x)-one; s = fabs(x)-1.0;
P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6))))); 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))))); Q = 1.0+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
if (hx >= 0) { if (hx >= 0) {
z = one-erx; z = 1.0-erx;
return z - P/Q; return z - P/Q;
} else { } else {
z = erx+P/Q; z = erx+P/Q;
return one+z; return 1.0 + z;
} }
} }
if (ix < 0x403c0000) { /* |x| < 28 */ if (ix < 0x403c0000) { /* |x| < 28 */
x = fabs(x); x = fabs(x);
s = one/(x*x); s = 1.0/(x*x);
if (ix < 0x4006DB6D) { /* |x| < 1/.35 ~ 2.857143*/ if (ix < 0x4006DB6D) { /* |x| < 1/.35 ~ 2.857143*/
R = ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*( R = ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
ra5+s*(ra6+s*ra7)))))); ra5+s*(ra6+s*ra7))))));
S = one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*( S = 1.0+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
sa5+s*(sa6+s*(sa7+s*sa8))))))); sa5+s*(sa6+s*(sa7+s*sa8)))))));
} else { /* |x| >= 1/.35 ~ 2.857143 */ } else { /* |x| >= 1/.35 ~ 2.857143 */
if (hx < 0 && ix >= 0x40180000) /* x < -6 */ if (hx < 0 && ix >= 0x40180000) /* x < -6 */
return two-tiny; return 2.0 - tiny;
R = rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*( R = rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
rb5+s*rb6))))); rb5+s*rb6)))));
S = one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*( S = 1.0+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
sb5+s*(sb6+s*sb7)))))); sb5+s*(sb6+s*sb7))))));
} }
z = x; z = x;
@ -298,9 +295,9 @@ double erfc(double x)
r = exp(-z*z-0.5625)*exp((z-x)*(z+x)+R/S); r = exp(-z*z-0.5625)*exp((z-x)*(z+x)+R/S);
if (hx > 0) if (hx > 0)
return r/x; return r/x;
return two-r/x; return 2.0 - r/x;
} }
if (hx > 0) if (hx > 0)
return tiny*tiny; return tiny*tiny;
return two-tiny; return 2.0 - tiny;
} }

View File

@ -17,9 +17,6 @@
static const float static const float
tiny = 1e-30, tiny = 1e-30,
half = 5.0000000000e-01, /* 0x3F000000 */
one = 1.0000000000e+00, /* 0x3F800000 */
two = 2.0000000000e+00, /* 0x40000000 */
/* c = (subfloat)0.84506291151 */ /* c = (subfloat)0.84506291151 */
erx = 8.4506291151e-01, /* 0x3f58560b */ erx = 8.4506291151e-01, /* 0x3f58560b */
/* /*
@ -100,7 +97,7 @@ float erff(float x)
if (ix >= 0x7f800000) { if (ix >= 0x7f800000) {
/* erf(nan)=nan, erf(+-inf)=+-1 */ /* erf(nan)=nan, erf(+-inf)=+-1 */
i = ((uint32_t)hx>>31)<<1; i = ((uint32_t)hx>>31)<<1;
return (float)(1-i)+one/x; return (float)(1-i)+1.0f/x;
} }
if (ix < 0x3f580000) { /* |x| < 0.84375 */ if (ix < 0x3f580000) { /* |x| < 0.84375 */
if (ix < 0x31800000) { /* |x| < 2**-28 */ if (ix < 0x31800000) { /* |x| < 2**-28 */
@ -111,42 +108,42 @@ float erff(float x)
} }
z = x*x; z = x*x;
r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4))); r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5)))); s = 1.0f+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
y = r/s; y = r/s;
return x + x*y; return x + x*y;
} }
if (ix < 0x3fa00000) { /* 0.84375 <= |x| < 1.25 */ if (ix < 0x3fa00000) { /* 0.84375 <= |x| < 1.25 */
s = fabsf(x)-one; s = fabsf(x)-1.0f;
P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6))))); 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))))); Q = 1.0f+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
if (hx >= 0) if (hx >= 0)
return erx + P/Q; return erx + P/Q;
return -erx - P/Q; return -erx - P/Q;
} }
if (ix >= 0x40c00000) { /* inf > |x| >= 6 */ if (ix >= 0x40c00000) { /* inf > |x| >= 6 */
if (hx >= 0) if (hx >= 0)
return one - tiny; return 1.0f - tiny;
return tiny - one; return tiny - 1.0f;
} }
x = fabsf(x); x = fabsf(x);
s = one/(x*x); s = 1.0f/(x*x);
if (ix < 0x4036DB6E) { /* |x| < 1/0.35 */ if (ix < 0x4036DB6E) { /* |x| < 1/0.35 */
R = ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*( R = ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
ra5+s*(ra6+s*ra7)))))); ra5+s*(ra6+s*ra7))))));
S = one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*( S = 1.0f+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
sa5+s*(sa6+s*(sa7+s*sa8))))))); sa5+s*(sa6+s*(sa7+s*sa8)))))));
} else { /* |x| >= 1/0.35 */ } else { /* |x| >= 1/0.35 */
R = rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*( R = rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
rb5+s*rb6))))); rb5+s*rb6)))));
S = one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*( S = 1.0f+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
sb5+s*(sb6+s*sb7)))))); sb5+s*(sb6+s*sb7))))));
} }
GET_FLOAT_WORD(ix, x); GET_FLOAT_WORD(ix, x);
SET_FLOAT_WORD(z, ix&0xfffff000); SET_FLOAT_WORD(z, ix&0xfffff000);
r = expf(-z*z - 0.5625f) * expf((z-x)*(z+x) + R/S); r = expf(-z*z - 0.5625f) * expf((z-x)*(z+x) + R/S);
if (hx >= 0) if (hx >= 0)
return one - r/x; return 1.0f - r/x;
return r/x - one; return r/x - 1.0f;
} }
float erfcf(float x) float erfcf(float x)
@ -158,50 +155,50 @@ float erfcf(float x)
ix = hx & 0x7fffffff; ix = hx & 0x7fffffff;
if (ix >= 0x7f800000) { if (ix >= 0x7f800000) {
/* erfc(nan)=nan, erfc(+-inf)=0,2 */ /* erfc(nan)=nan, erfc(+-inf)=0,2 */
return (float)(((uint32_t)hx>>31)<<1) + one/x; return (float)(((uint32_t)hx>>31)<<1) + 1.0f/x;
} }
if (ix < 0x3f580000) { /* |x| < 0.84375 */ if (ix < 0x3f580000) { /* |x| < 0.84375 */
if (ix < 0x23800000) /* |x| < 2**-56 */ if (ix < 0x23800000) /* |x| < 2**-56 */
return one - x; return 1.0f - x;
z = x*x; z = x*x;
r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4))); r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5)))); s = 1.0f+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
y = r/s; y = r/s;
if (hx < 0x3e800000) { /* x<1/4 */ if (hx < 0x3e800000) { /* x<1/4 */
return one - (x+x*y); return 1.0f - (x+x*y);
} else { } else {
r = x*y; r = x*y;
r += (x-half); r += (x-0.5f);
return half - r ; return 0.5f - r ;
} }
} }
if (ix < 0x3fa00000) { /* 0.84375 <= |x| < 1.25 */ if (ix < 0x3fa00000) { /* 0.84375 <= |x| < 1.25 */
s = fabsf(x)-one; s = fabsf(x)-1.0f;
P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6))))); 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))))); Q = 1.0f+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
if(hx >= 0) { if(hx >= 0) {
z = one - erx; z = 1.0f - erx;
return z - P/Q; return z - P/Q;
} else { } else {
z = erx + P/Q; z = erx + P/Q;
return one + z; return 1.0f + z;
} }
} }
if (ix < 0x41e00000) { /* |x| < 28 */ if (ix < 0x41e00000) { /* |x| < 28 */
x = fabsf(x); x = fabsf(x);
s = one/(x*x); s = 1.0f/(x*x);
if (ix < 0x4036DB6D) { /* |x| < 1/.35 ~ 2.857143*/ if (ix < 0x4036DB6D) { /* |x| < 1/.35 ~ 2.857143*/
R = ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*( R = ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
ra5+s*(ra6+s*ra7)))))); ra5+s*(ra6+s*ra7))))));
S = one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*( S = 1.0f+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
sa5+s*(sa6+s*(sa7+s*sa8))))))); sa5+s*(sa6+s*(sa7+s*sa8)))))));
} else { /* |x| >= 1/.35 ~ 2.857143 */ } else { /* |x| >= 1/.35 ~ 2.857143 */
if (hx < 0 && ix >= 0x40c00000) /* x < -6 */ if (hx < 0 && ix >= 0x40c00000) /* x < -6 */
return two-tiny; return 2.0f-tiny;
R = rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*( R = rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
rb5+s*rb6))))); rb5+s*rb6)))));
S = one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*( S = 1.0f+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
sb5+s*(sb6+s*sb7)))))); sb5+s*(sb6+s*sb7))))));
} }
GET_FLOAT_WORD(ix, x); GET_FLOAT_WORD(ix, x);
@ -209,9 +206,9 @@ float erfcf(float x)
r = expf(-z*z - 0.5625f) * expf((z-x)*(z+x) + R/S); r = expf(-z*z - 0.5625f) * expf((z-x)*(z+x) + R/S);
if (hx > 0) if (hx > 0)
return r/x; return r/x;
return two - r/x; return 2.0f - r/x;
} }
if (hx > 0) if (hx > 0)
return tiny*tiny; return tiny*tiny;
return two - tiny; return 2.0f - tiny;
} }

View File

@ -108,9 +108,6 @@ long double erfl(long double x)
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
static const long double static const long double
tiny = 1e-4931L, tiny = 1e-4931L,
half = 0.5L,
one = 1.0L,
two = 2.0L,
/* c = (float)0.84506291151 */ /* c = (float)0.84506291151 */
erx = 0.845062911510467529296875L, erx = 0.845062911510467529296875L,
@ -253,7 +250,7 @@ long double erfl(long double x)
if (ix >= 0x7fff) { /* erf(nan)=nan */ if (ix >= 0x7fff) { /* erf(nan)=nan */
i = ((se & 0xffff) >> 15) << 1; i = ((se & 0xffff) >> 15) << 1;
return (long double)(1 - i) + one / x; /* erf(+-inf)=+-1 */ return (long double)(1 - i) + 1.0 / x; /* erf(+-inf)=+-1 */
} }
ix = (ix << 16) | (i0 >> 16); ix = (ix << 16) | (i0 >> 16);
@ -272,7 +269,7 @@ long double erfl(long double x)
return x + x * y; return x + x * y;
} }
if (ix < 0x3fffa000) { /* 0.84375 <= |x| < 1.25 */ if (ix < 0x3fffa000) { /* 0.84375 <= |x| < 1.25 */
s = fabsl (x) - one; s = fabsl(x) - 1.0;
P = pa[0] + s * (pa[1] + s * (pa[2] + P = pa[0] + s * (pa[1] + s * (pa[2] +
s * (pa[3] + s * (pa[4] + s * (pa[5] + s * (pa[6] + s * pa[7])))))); s * (pa[3] + s * (pa[4] + s * (pa[5] + s * (pa[6] + s * pa[7]))))));
Q = qa[0] + s * (qa[1] + s * (qa[2] + Q = qa[0] + s * (qa[1] + s * (qa[2] +
@ -283,11 +280,11 @@ long double erfl(long double x)
} }
if (ix >= 0x4001d555) { /* inf > |x| >= 6.6666259765625 */ if (ix >= 0x4001d555) { /* inf > |x| >= 6.6666259765625 */
if ((se & 0x8000) == 0) if ((se & 0x8000) == 0)
return one - tiny; return 1.0 - tiny;
return tiny - one; return tiny - 1.0;
} }
x = fabsl (x); x = fabsl (x);
s = one / (x * x); s = 1.0 / (x * x);
if (ix < 0x4000b6db) { /* 1.25 <= |x| < 2.85711669921875 ~ 1/.35 */ 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] + 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 * (ra[5] + s * (ra[6] + s * (ra[7] + s * ra[8])))))));
@ -305,8 +302,8 @@ long double erfl(long double x)
SET_LDOUBLE_WORDS(z, i, i0, i1); SET_LDOUBLE_WORDS(z, i, i0, i1);
r = expl(-z * z - 0.5625) * expl((z - x) * (z + x) + R / S); r = expl(-z * z - 0.5625) * expl((z - x) * (z + x) + R / S);
if ((se & 0x8000) == 0) if ((se & 0x8000) == 0)
return one - r / x; return 1.0 - r / x;
return r / x - one; return r / x - 1.0;
} }
long double erfcl(long double x) long double erfcl(long double x)
@ -318,13 +315,13 @@ long double erfcl(long double x)
GET_LDOUBLE_WORDS(se, i0, i1, x); GET_LDOUBLE_WORDS(se, i0, i1, x);
ix = se & 0x7fff; ix = se & 0x7fff;
if (ix >= 0x7fff) { /* erfc(nan) = nan, erfc(+-inf) = 0,2 */ if (ix >= 0x7fff) { /* erfc(nan) = nan, erfc(+-inf) = 0,2 */
return (long double)(((se & 0xffff) >> 15) << 1) + one / x; return (long double)(((se & 0xffff) >> 15) << 1) + 1.0 / x;
} }
ix = (ix << 16) | (i0 >> 16); ix = (ix << 16) | (i0 >> 16);
if (ix < 0x3ffed800) { /* |x| < 0.84375 */ if (ix < 0x3ffed800) { /* |x| < 0.84375 */
if (ix < 0x3fbe0000) /* |x| < 2**-65 */ if (ix < 0x3fbe0000) /* |x| < 2**-65 */
return one - x; return 1.0 - x;
z = x * x; z = x * x;
r = pp[0] + z * (pp[1] + r = pp[0] + z * (pp[1] +
z * (pp[2] + z * (pp[3] + z * (pp[4] + z * pp[5])))); z * (pp[2] + z * (pp[3] + z * (pp[4] + z * pp[5]))));
@ -332,27 +329,27 @@ long double erfcl(long double x)
z * (qq[2] + z * (qq[3] + z * (qq[4] + z * (qq[5] + z))))); z * (qq[2] + z * (qq[3] + z * (qq[4] + z * (qq[5] + z)))));
y = r / s; y = r / s;
if (ix < 0x3ffd8000) /* x < 1/4 */ if (ix < 0x3ffd8000) /* x < 1/4 */
return one - (x + x * y); return 1.0 - (x + x * y);
r = x * y; r = x * y;
r += x - half; r += x - 0.5L;
return half - r; return 0.5L - r;
} }
if (ix < 0x3fffa000) { /* 0.84375 <= |x| < 1.25 */ if (ix < 0x3fffa000) { /* 0.84375 <= |x| < 1.25 */
s = fabsl (x) - one; s = fabsl(x) - 1.0;
P = pa[0] + s * (pa[1] + s * (pa[2] + P = pa[0] + s * (pa[1] + s * (pa[2] +
s * (pa[3] + s * (pa[4] + s * (pa[5] + s * (pa[6] + s * pa[7])))))); s * (pa[3] + s * (pa[4] + s * (pa[5] + s * (pa[6] + s * pa[7]))))));
Q = qa[0] + s * (qa[1] + s * (qa[2] + Q = qa[0] + s * (qa[1] + s * (qa[2] +
s * (qa[3] + s * (qa[4] + s * (qa[5] + s * (qa[6] + s)))))); s * (qa[3] + s * (qa[4] + s * (qa[5] + s * (qa[6] + s))))));
if ((se & 0x8000) == 0) { if ((se & 0x8000) == 0) {
z = one - erx; z = 1.0 - erx;
return z - P / Q; return z - P / Q;
} }
z = erx + P / Q; z = erx + P / Q;
return one + z; return 1.0 + z;
} }
if (ix < 0x4005d600) { /* |x| < 107 */ if (ix < 0x4005d600) { /* |x| < 107 */
x = fabsl(x); x = fabsl(x);
s = one / (x * x); s = 1.0 / (x * x);
if (ix < 0x4000b6db) { /* 1.25 <= |x| < 2.85711669921875 ~ 1/.35 */ 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] + 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 * (ra[5] + s * (ra[6] + s * (ra[7] + s * ra[8])))))));
@ -365,7 +362,7 @@ long double erfcl(long double x)
s * (sb[5] + s * (sb[6] + s)))))); s * (sb[5] + s * (sb[6] + s))))));
} else { /* 107 > |x| >= 6.666 */ } else { /* 107 > |x| >= 6.666 */
if (se & 0x8000) if (se & 0x8000)
return two - tiny;/* x < -6.666 */ return 2.0 - tiny;/* x < -6.666 */
R = rc[0] + s * (rc[1] + s * (rc[2] + s * (rc[3] + R = rc[0] + s * (rc[1] + s * (rc[2] + s * (rc[3] +
s * (rc[4] + s * rc[5])))); s * (rc[4] + s * rc[5]))));
S = sc[0] + s * (sc[1] + s * (sc[2] + s * (sc[3] + S = sc[0] + s * (sc[1] + s * (sc[2] + s * (sc[3] +
@ -376,15 +373,14 @@ long double erfcl(long double x)
i1 = 0; i1 = 0;
i0 &= 0xffffff00; i0 &= 0xffffff00;
SET_LDOUBLE_WORDS(z, hx, i0, i1); SET_LDOUBLE_WORDS(z, hx, i0, i1);
r = expl (-z * z - 0.5625) * r = expl(-z * z - 0.5625) * expl((z - x) * (z + x) + R / S);
expl ((z - x) * (z + x) + R / S);
if ((se & 0x8000) == 0) if ((se & 0x8000) == 0)
return r / x; return r / x;
return two - r / x; return 2.0 - r / x;
} }
if ((se & 0x8000) == 0) if ((se & 0x8000) == 0)
return tiny * tiny; return tiny * tiny;
return two - tiny; return 2.0 - tiny;
} }
#endif #endif

View File

@ -74,7 +74,6 @@
#include "libm.h" #include "libm.h"
static const double static const double
one = 1.0,
halF[2] = {0.5,-0.5,}, halF[2] = {0.5,-0.5,},
huge = 1.0e+300, huge = 1.0e+300,
o_threshold = 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */ o_threshold = 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */
@ -134,8 +133,8 @@ double exp(double x)
STRICT_ASSIGN(double, x, hi - lo); STRICT_ASSIGN(double, x, hi - lo);
} else if(hx < 0x3e300000) { /* |x| < 2**-28 */ } else if(hx < 0x3e300000) { /* |x| < 2**-28 */
/* raise inexact */ /* raise inexact */
if (huge+x > one) if (huge+x > 1.0)
return one+x; return 1.0+x;
} else } else
k = 0; k = 0;
@ -147,8 +146,8 @@ double exp(double x)
INSERT_WORDS(twopk, 0x3ff00000+((k+1000)<<20), 0); INSERT_WORDS(twopk, 0x3ff00000+((k+1000)<<20), 0);
c = x - t*(P1+t*(P2+t*(P3+t*(P4+t*P5)))); c = x - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
if (k == 0) if (k == 0)
return one - ((x*c)/(c-2.0) - x); return 1.0 - ((x*c)/(c-2.0) - x);
y = one-((lo-(x*c)/(2.0-c))-hi); y = 1.0-((lo-(x*c)/(2.0-c))-hi);
if (k < -1021) if (k < -1021)
return y*twopk*twom1000; return y*twopk*twom1000;
if (k == 1024) if (k == 1024)

View File

@ -16,7 +16,6 @@
#include "libm.h" #include "libm.h"
static const float static const float
one = 1.0,
halF[2] = {0.5,-0.5,}, halF[2] = {0.5,-0.5,},
huge = 1.0e+30, huge = 1.0e+30,
o_threshold = 8.8721679688e+01, /* 0x42b17180 */ o_threshold = 8.8721679688e+01, /* 0x42b17180 */
@ -72,8 +71,8 @@ float expf(float x)
STRICT_ASSIGN(float, x, hi - lo); STRICT_ASSIGN(float, x, hi - lo);
} else if(hx < 0x39000000) { /* |x|<2**-14 */ } else if(hx < 0x39000000) { /* |x|<2**-14 */
/* raise inexact */ /* raise inexact */
if (huge+x > one) if (huge+x > 1.0f)
return one + x; return 1.0f + x;
} else } else
k = 0; k = 0;
@ -85,8 +84,8 @@ float expf(float x)
SET_FLOAT_WORD(twopk, 0x3f800000+((k+100)<<23)); SET_FLOAT_WORD(twopk, 0x3f800000+((k+100)<<23));
c = x - t*(P1+t*P2); c = x - t*(P1+t*P2);
if (k == 0) if (k == 0)
return one - ((x*c)/(c - 2.0f) - x); return 1.0f - ((x*c)/(c - 2.0f) - x);
y = one - ((lo - (x*c)/(2.0f - c)) - hi); y = 1.0f - ((lo - (x*c)/(2.0f - c)) - hi);
if (k < -125) if (k < -125)
return y*twopk*twom100; return y*twopk*twom100;
if (k == 128) if (k == 128)

View File

@ -107,7 +107,6 @@
#include "libm.h" #include "libm.h"
static const double static const double
one = 1.0,
huge = 1.0e+300, huge = 1.0e+300,
tiny = 1.0e-300, tiny = 1.0e-300,
o_threshold = 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */ o_threshold = 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */
@ -148,7 +147,7 @@ double expm1(double x)
if (xsb != 0) { /* x < -56*ln2, return -1.0 with inexact */ if (xsb != 0) { /* x < -56*ln2, return -1.0 with inexact */
/* raise inexact */ /* raise inexact */
if(x+tiny<0.0) if(x+tiny<0.0)
return tiny-one; /* return -1 */ return tiny-1.0; /* return -1 */
} }
} }
@ -182,7 +181,7 @@ double expm1(double x)
/* x is now in primary range */ /* x is now in primary range */
hfx = 0.5*x; hfx = 0.5*x;
hxs = x*hfx; hxs = x*hfx;
r1 = one+hxs*(Q1+hxs*(Q2+hxs*(Q3+hxs*(Q4+hxs*Q5)))); r1 = 1.0+hxs*(Q1+hxs*(Q2+hxs*(Q3+hxs*(Q4+hxs*Q5))));
t = 3.0-r1*hfx; t = 3.0-r1*hfx;
e = hxs*((r1-t)/(6.0 - x*t)); e = hxs*((r1-t)/(6.0 - x*t));
if (k == 0) /* c is 0 */ if (k == 0) /* c is 0 */
@ -195,17 +194,17 @@ double expm1(double x)
if (k == 1) { if (k == 1) {
if (x < -0.25) if (x < -0.25)
return -2.0*(e-(x+0.5)); return -2.0*(e-(x+0.5));
return one+2.0*(x-e); return 1.0+2.0*(x-e);
} }
if (k <= -2 || k > 56) { /* suffice to return exp(x)-1 */ if (k <= -2 || k > 56) { /* suffice to return exp(x)-1 */
y = one - (e-x); y = 1.0 - (e-x);
if (k == 1024) if (k == 1024)
y = y*2.0*0x1p1023; y = y*2.0*0x1p1023;
else else
y = y*twopk; y = y*twopk;
return y - one; return y - 1.0;
} }
t = one; t = 1.0;
if (k < 20) { if (k < 20) {
SET_HIGH_WORD(t, 0x3ff00000 - (0x200000>>k)); /* t=1-2^-k */ SET_HIGH_WORD(t, 0x3ff00000 - (0x200000>>k)); /* t=1-2^-k */
y = t-(e-x); y = t-(e-x);
@ -213,7 +212,7 @@ double expm1(double x)
} else { } else {
SET_HIGH_WORD(t, ((0x3ff-k)<<20)); /* 2^-k */ SET_HIGH_WORD(t, ((0x3ff-k)<<20)); /* 2^-k */
y = x-(e+t); y = x-(e+t);
y += one; y += 1.0;
y = y*twopk; y = y*twopk;
} }
return y; return y;

View File

@ -16,7 +16,6 @@
#include "libm.h" #include "libm.h"
static const float static const float
one = 1.0,
huge = 1.0e+30, huge = 1.0e+30,
tiny = 1.0e-30, tiny = 1.0e-30,
o_threshold = 8.8721679688e+01, /* 0x42b17180 */ o_threshold = 8.8721679688e+01, /* 0x42b17180 */
@ -54,7 +53,7 @@ float expm1f(float x)
if (xsb != 0) { /* x < -27*ln2 */ if (xsb != 0) { /* x < -27*ln2 */
/* raise inexact */ /* raise inexact */
if (x+tiny < 0.0f) if (x+tiny < 0.0f)
return tiny-one; /* return -1 */ return tiny-1.0f; /* return -1 */
} }
} }
@ -87,7 +86,7 @@ float expm1f(float x)
/* x is now in primary range */ /* x is now in primary range */
hfx = 0.5f*x; hfx = 0.5f*x;
hxs = x*hfx; hxs = x*hfx;
r1 = one+hxs*(Q1+hxs*Q2); r1 = 1.0f+hxs*(Q1+hxs*Q2);
t = 3.0f - r1*hfx; t = 3.0f - r1*hfx;
e = hxs*((r1-t)/(6.0f - x*t)); e = hxs*((r1-t)/(6.0f - x*t));
if (k == 0) /* c is 0 */ if (k == 0) /* c is 0 */
@ -100,17 +99,17 @@ float expm1f(float x)
if (k == 1) { if (k == 1) {
if (x < -0.25f) if (x < -0.25f)
return -2.0f*(e-(x+0.5f)); return -2.0f*(e-(x+0.5f));
return one + 2.0f*(x-e); return 1.0f + 2.0f*(x-e);
} }
if (k <= -2 || k > 56) { /* suffice to return exp(x)-1 */ if (k <= -2 || k > 56) { /* suffice to return exp(x)-1 */
y = one - (e - x); y = 1.0f - (e - x);
if (k == 128) if (k == 128)
y = y*2.0f*0x1p127f; y = y*2.0f*0x1p127f;
else else
y = y*twopk; y = y*twopk;
return y - one; return y - 1.0f;
} }
t = one; t = 1.0f;
if (k < 23) { if (k < 23) {
SET_FLOAT_WORD(t, 0x3f800000 - (0x1000000>>k)); /* t=1-2^-k */ SET_FLOAT_WORD(t, 0x3f800000 - (0x1000000>>k)); /* t=1-2^-k */
y = t - (e - x); y = t - (e - x);
@ -118,7 +117,7 @@ float expm1f(float x)
} else { } else {
SET_FLOAT_WORD(t, (0x7f-k)<<23); /* 2^-k */ SET_FLOAT_WORD(t, (0x7f-k)<<23); /* 2^-k */
y = x - (e + t); y = x - (e + t);
y += one; y += 1.0f;
y = y*twopk; y = y*twopk;
} }
return y; return y;

View File

@ -60,7 +60,6 @@ static double pzero(double), qzero(double);
static const double static const double
huge = 1e300, huge = 1e300,
one = 1.0,
invsqrtpi = 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */ invsqrtpi = 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
tpi = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */ tpi = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
/* R0/S0 on [0, 2.00] */ /* R0/S0 on [0, 2.00] */
@ -73,8 +72,6 @@ S02 = 1.16926784663337450260e-04, /* 0x3F1EA6D2, 0xDD57DBF4 */
S03 = 5.13546550207318111446e-07, /* 0x3EA13B54, 0xCE84D5A9 */ S03 = 5.13546550207318111446e-07, /* 0x3EA13B54, 0xCE84D5A9 */
S04 = 1.16614003333790000205e-09; /* 0x3E1408BC, 0xF4745D8F */ S04 = 1.16614003333790000205e-09; /* 0x3E1408BC, 0xF4745D8F */
static const double zero = 0.0;
double j0(double x) double j0(double x)
{ {
double z, s,c,ss,cc,r,u,v; double z, s,c,ss,cc,r,u,v;
@ -83,7 +80,7 @@ double j0(double x)
GET_HIGH_WORD(hx, x); GET_HIGH_WORD(hx, x);
ix = hx & 0x7fffffff; ix = hx & 0x7fffffff;
if (ix >= 0x7ff00000) if (ix >= 0x7ff00000)
return one/(x*x); return 1.0/(x*x);
x = fabs(x); x = fabs(x);
if (ix >= 0x40000000) { /* |x| >= 2.0 */ if (ix >= 0x40000000) { /* |x| >= 2.0 */
s = sin(x); s = sin(x);
@ -92,7 +89,7 @@ double j0(double x)
cc = s+c; cc = s+c;
if (ix < 0x7fe00000) { /* make sure x+x does not overflow */ if (ix < 0x7fe00000) { /* make sure x+x does not overflow */
z = -cos(x+x); z = -cos(x+x);
if ((s*c) < zero) if (s*c < 0.0)
cc = z/ss; cc = z/ss;
else else
ss = z/cc; ss = z/cc;
@ -112,20 +109,20 @@ double j0(double x)
} }
if (ix < 0x3f200000) { /* |x| < 2**-13 */ if (ix < 0x3f200000) { /* |x| < 2**-13 */
/* raise inexact if x != 0 */ /* raise inexact if x != 0 */
if (huge+x > one) { if (huge+x > 1.0) {
if (ix < 0x3e400000) /* |x| < 2**-27 */ if (ix < 0x3e400000) /* |x| < 2**-27 */
return one; return 1.0;
return one - 0.25*x*x; return 1.0 - 0.25*x*x;
} }
} }
z = x*x; z = x*x;
r = z*(R02+z*(R03+z*(R04+z*R05))); r = z*(R02+z*(R03+z*(R04+z*R05)));
s = one+z*(S01+z*(S02+z*(S03+z*S04))); s = 1.0+z*(S01+z*(S02+z*(S03+z*S04)));
if (ix < 0x3FF00000) { /* |x| < 1.00 */ if (ix < 0x3FF00000) { /* |x| < 1.00 */
return one + z*(-0.25+(r/s)); return 1.0 + z*(-0.25+(r/s));
} else { } else {
u = 0.5*x; u = 0.5*x;
return (one+u)*(one-u) + z*(r/s); return (1.0+u)*(1.0-u) + z*(r/s);
} }
} }
@ -151,11 +148,11 @@ double y0(double x)
ix = 0x7fffffff & hx; ix = 0x7fffffff & hx;
/* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0 */ /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0 */
if (ix >= 0x7ff00000) if (ix >= 0x7ff00000)
return one/(x+x*x); return 1.0/(x+x*x);
if ((ix|lx) == 0) if ((ix|lx) == 0)
return -one/zero; return -1.0/0.0;
if (hx < 0) if (hx < 0)
return zero/zero; return 0.0/0.0;
if (ix >= 0x40000000) { /* |x| >= 2.0 */ if (ix >= 0x40000000) { /* |x| >= 2.0 */
/* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0)) /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0))
* where x0 = x-pi/4 * where x0 = x-pi/4
@ -178,7 +175,7 @@ double y0(double x)
*/ */
if (ix < 0x7fe00000) { /* make sure x+x does not overflow */ if (ix < 0x7fe00000) { /* make sure x+x does not overflow */
z = -cos(x+x); z = -cos(x+x);
if (s*c < zero) if (s*c < 0.0)
cc = z/ss; cc = z/ss;
else else
ss = z/cc; ss = z/cc;
@ -197,7 +194,7 @@ double y0(double x)
} }
z = x*x; z = x*x;
u = u00+z*(u01+z*(u02+z*(u03+z*(u04+z*(u05+z*u06))))); u = u00+z*(u01+z*(u02+z*(u03+z*(u04+z*(u05+z*u06)))));
v = one+z*(v01+z*(v02+z*(v03+z*v04))); v = 1.0+z*(v01+z*(v02+z*(v03+z*v04)));
return u/v + tpi*(j0(x)*log(x)); return u/v + tpi*(j0(x)*log(x));
} }
@ -286,10 +283,10 @@ static double pzero(double x)
else if (ix >= 0x40122E8B){p = pR5; q = pS5;} else if (ix >= 0x40122E8B){p = pR5; q = pS5;}
else if (ix >= 0x4006DB6D){p = pR3; q = pS3;} else if (ix >= 0x4006DB6D){p = pR3; q = pS3;}
else if (ix >= 0x40000000){p = pR2; q = pS2;} else if (ix >= 0x40000000){p = pR2; q = pS2;}
z = one/(x*x); z = 1.0/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); s = 1.0+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
return one + r/s; return 1.0 + r/s;
} }
@ -382,8 +379,8 @@ static double qzero(double x)
else if (ix >= 0x40122E8B){p = qR5; q = qS5;} else if (ix >= 0x40122E8B){p = qR5; q = qS5;}
else if (ix >= 0x4006DB6D){p = qR3; q = qS3;} else if (ix >= 0x4006DB6D){p = qR3; q = qS3;}
else if (ix >= 0x40000000){p = qR2; q = qS2;} else if (ix >= 0x40000000){p = qR2; q = qS2;}
z = one/(x*x); z = 1.0/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5]))))); s = 1.0+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
return (-.125 + r/s)/x; return (-.125 + r/s)/x;
} }

View File

@ -19,7 +19,6 @@ static float pzerof(float), qzerof(float);
static const float static const float
huge = 1e30, huge = 1e30,
one = 1.0,
invsqrtpi = 5.6418961287e-01, /* 0x3f106ebb */ invsqrtpi = 5.6418961287e-01, /* 0x3f106ebb */
tpi = 6.3661974669e-01, /* 0x3f22f983 */ tpi = 6.3661974669e-01, /* 0x3f22f983 */
/* R0/S0 on [0, 2.00] */ /* R0/S0 on [0, 2.00] */
@ -32,8 +31,6 @@ S02 = 1.1692678527e-04, /* 0x38f53697 */
S03 = 5.1354652442e-07, /* 0x3509daa6 */ S03 = 5.1354652442e-07, /* 0x3509daa6 */
S04 = 1.1661400734e-09; /* 0x30a045e8 */ S04 = 1.1661400734e-09; /* 0x30a045e8 */
static const float zero = 0.0;
float j0f(float x) float j0f(float x)
{ {
float z, s,c,ss,cc,r,u,v; float z, s,c,ss,cc,r,u,v;
@ -42,7 +39,7 @@ float j0f(float x)
GET_FLOAT_WORD(hx, x); GET_FLOAT_WORD(hx, x);
ix = hx & 0x7fffffff; ix = hx & 0x7fffffff;
if (ix >= 0x7f800000) if (ix >= 0x7f800000)
return one/(x*x); return 1.0f/(x*x);
x = fabsf(x); x = fabsf(x);
if (ix >= 0x40000000) { /* |x| >= 2.0 */ if (ix >= 0x40000000) { /* |x| >= 2.0 */
s = sinf(x); s = sinf(x);
@ -51,7 +48,7 @@ float j0f(float x)
cc = s+c; cc = s+c;
if (ix < 0x7f000000) { /* make sure x+x does not overflow */ if (ix < 0x7f000000) { /* make sure x+x does not overflow */
z = -cosf(x+x); z = -cosf(x+x);
if (s*c < zero) if (s*c < 0.0f)
cc = z/ss; cc = z/ss;
else else
ss = z/cc; ss = z/cc;
@ -71,20 +68,20 @@ float j0f(float x)
} }
if (ix < 0x39000000) { /* |x| < 2**-13 */ if (ix < 0x39000000) { /* |x| < 2**-13 */
/* raise inexact if x != 0 */ /* raise inexact if x != 0 */
if (huge+x > one) { if (huge+x > 1.0f) {
if (ix < 0x32000000) /* |x| < 2**-27 */ if (ix < 0x32000000) /* |x| < 2**-27 */
return one; return 1.0f;
return one - 0.25f*x*x; return 1.0f - 0.25f*x*x;
} }
} }
z = x*x; z = x*x;
r = z*(R02+z*(R03+z*(R04+z*R05))); r = z*(R02+z*(R03+z*(R04+z*R05)));
s = one+z*(S01+z*(S02+z*(S03+z*S04))); s = 1.0f+z*(S01+z*(S02+z*(S03+z*S04)));
if(ix < 0x3F800000) { /* |x| < 1.00 */ if(ix < 0x3F800000) { /* |x| < 1.00 */
return one + z*(-0.25f + (r/s)); return 1.0f + z*(-0.25f + (r/s));
} else { } else {
u = 0.5f*x; u = 0.5f*x;
return (one+u)*(one-u) + z*(r/s); return (1.0f+u)*(1.0f-u) + z*(r/s);
} }
} }
@ -110,11 +107,11 @@ float y0f(float x)
ix = 0x7fffffff & hx; ix = 0x7fffffff & hx;
/* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0 */ /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0 */
if (ix >= 0x7f800000) if (ix >= 0x7f800000)
return one/(x+x*x); return 1.0f/(x+x*x);
if (ix == 0) if (ix == 0)
return -one/zero; return -1.0f/0.0f;
if (hx < 0) if (hx < 0)
return zero/zero; return 0.0f/0.0f;
if (ix >= 0x40000000) { /* |x| >= 2.0 */ if (ix >= 0x40000000) { /* |x| >= 2.0 */
/* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0)) /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0))
* where x0 = x-pi/4 * where x0 = x-pi/4
@ -137,7 +134,7 @@ float y0f(float x)
*/ */
if (ix < 0x7f000000) { /* make sure x+x not overflow */ if (ix < 0x7f000000) { /* make sure x+x not overflow */
z = -cosf(x+x); z = -cosf(x+x);
if (s*c < zero) if (s*c < 0.0f)
cc = z/ss; cc = z/ss;
else else
ss = z/cc; ss = z/cc;
@ -156,7 +153,7 @@ float y0f(float x)
} }
z = x*x; z = x*x;
u = u00+z*(u01+z*(u02+z*(u03+z*(u04+z*(u05+z*u06))))); u = u00+z*(u01+z*(u02+z*(u03+z*(u04+z*(u05+z*u06)))));
v = one+z*(v01+z*(v02+z*(v03+z*v04))); v = 1.0f+z*(v01+z*(v02+z*(v03+z*v04)));
return u/v + tpi*(j0f(x)*logf(x)); return u/v + tpi*(j0f(x)*logf(x));
} }
@ -244,10 +241,10 @@ static float pzerof(float x)
else if (ix >= 0x40f71c58){p = pR5; q = pS5;} else if (ix >= 0x40f71c58){p = pR5; q = pS5;}
else if (ix >= 0x4036db68){p = pR3; q = pS3;} else if (ix >= 0x4036db68){p = pR3; q = pS3;}
else if (ix >= 0x40000000){p = pR2; q = pS2;} else if (ix >= 0x40000000){p = pR2; q = pS2;}
z = one/(x*x); z = 1.0f/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); s = 1.0f+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
return one + r/s; return 1.0f + r/s;
} }
@ -340,8 +337,8 @@ static float qzerof(float x)
else if (ix >= 0x40f71c58){p = qR5; q = qS5;} else if (ix >= 0x40f71c58){p = qR5; q = qS5;}
else if (ix >= 0x4036db68){p = qR3; q = qS3;} else if (ix >= 0x4036db68){p = qR3; q = qS3;}
else if (ix >= 0x40000000){p = qR2; q = qS2;} else if (ix >= 0x40000000){p = qR2; q = qS2;}
z = one/(x*x); z = 1.0f/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5]))))); s = 1.0f+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
return (-.125f + r/s)/x; return (-.125f + r/s)/x;
} }

View File

@ -60,7 +60,6 @@ static double pone(double), qone(double);
static const double static const double
huge = 1e300, huge = 1e300,
one = 1.0,
invsqrtpi = 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */ invsqrtpi = 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
tpi = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */ tpi = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
/* R0/S0 on [0,2] */ /* R0/S0 on [0,2] */
@ -74,8 +73,6 @@ s03 = 1.17718464042623683263e-06, /* 0x3EB3BFF8, 0x333F8498 */
s04 = 5.04636257076217042715e-09, /* 0x3E35AC88, 0xC97DFF2C */ s04 = 5.04636257076217042715e-09, /* 0x3E35AC88, 0xC97DFF2C */
s05 = 1.23542274426137913908e-11; /* 0x3DAB2ACF, 0xCFB97ED8 */ s05 = 1.23542274426137913908e-11; /* 0x3DAB2ACF, 0xCFB97ED8 */
static const double zero = 0.0;
double j1(double x) double j1(double x)
{ {
double z,s,c,ss,cc,r,u,v,y; double z,s,c,ss,cc,r,u,v,y;
@ -84,7 +81,7 @@ double j1(double x)
GET_HIGH_WORD(hx, x); GET_HIGH_WORD(hx, x);
ix = hx & 0x7fffffff; ix = hx & 0x7fffffff;
if (ix >= 0x7ff00000) if (ix >= 0x7ff00000)
return one/x; return 1.0/x;
y = fabs(x); y = fabs(x);
if (ix >= 0x40000000) { /* |x| >= 2.0 */ if (ix >= 0x40000000) { /* |x| >= 2.0 */
s = sin(y); s = sin(y);
@ -93,7 +90,7 @@ double j1(double x)
cc = s-c; cc = s-c;
if (ix < 0x7fe00000) { /* make sure y+y not overflow */ if (ix < 0x7fe00000) { /* make sure y+y not overflow */
z = cos(y+y); z = cos(y+y);
if (s*c > zero) if (s*c > 0.0)
cc = z/ss; cc = z/ss;
else else
ss = z/cc; ss = z/cc;
@ -116,12 +113,12 @@ double j1(double x)
} }
if (ix < 0x3e400000) { /* |x| < 2**-27 */ if (ix < 0x3e400000) { /* |x| < 2**-27 */
/* raise inexact if x!=0 */ /* raise inexact if x!=0 */
if (huge+x > one) if (huge+x > 1.0)
return 0.5*x; return 0.5*x;
} }
z = x*x; z = x*x;
r = z*(r00+z*(r01+z*(r02+z*r03))); r = z*(r00+z*(r01+z*(r02+z*r03)));
s = one+z*(s01+z*(s02+z*(s03+z*(s04+z*s05)))); s = 1.0+z*(s01+z*(s02+z*(s03+z*(s04+z*s05))));
r *= x; r *= x;
return x*0.5 + r/s; return x*0.5 + r/s;
} }
@ -151,11 +148,11 @@ double y1(double x)
ix = 0x7fffffff & hx; ix = 0x7fffffff & hx;
/* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */ /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
if (ix >= 0x7ff00000) if (ix >= 0x7ff00000)
return one/(x+x*x); return 1.0/(x+x*x);
if ((ix|lx) == 0) if ((ix|lx) == 0)
return -one/zero; return -1.0/0.0;
if (hx < 0) if (hx < 0)
return zero/zero; return 0.0/0.0;
if (ix >= 0x40000000) { /* |x| >= 2.0 */ if (ix >= 0x40000000) { /* |x| >= 2.0 */
s = sin(x); s = sin(x);
c = cos(x); c = cos(x);
@ -163,7 +160,7 @@ double y1(double x)
cc = s-c; cc = s-c;
if (ix < 0x7fe00000) { /* make sure x+x not overflow */ if (ix < 0x7fe00000) { /* make sure x+x not overflow */
z = cos(x+x); z = cos(x+x);
if (s*c > zero) if (s*c > 0.0)
cc = z/ss; cc = z/ss;
else else
ss = z/cc; ss = z/cc;
@ -192,8 +189,8 @@ double y1(double x)
return -tpi/x; return -tpi/x;
z = x*x; z = x*x;
u = U0[0]+z*(U0[1]+z*(U0[2]+z*(U0[3]+z*U0[4]))); u = U0[0]+z*(U0[1]+z*(U0[2]+z*(U0[3]+z*U0[4])));
v = one+z*(V0[0]+z*(V0[1]+z*(V0[2]+z*(V0[3]+z*V0[4])))); v = 1.0+z*(V0[0]+z*(V0[1]+z*(V0[2]+z*(V0[3]+z*V0[4]))));
return x*(u/v) + tpi*(j1(x)*log(x)-one/x); return x*(u/v) + tpi*(j1(x)*log(x)-1.0/x);
} }
/* For x >= 8, the asymptotic expansions of pone is /* For x >= 8, the asymptotic expansions of pone is
@ -282,10 +279,10 @@ static double pone(double x)
else if (ix >= 0x40122E8B){p = pr5; q = ps5;} else if (ix >= 0x40122E8B){p = pr5; q = ps5;}
else if (ix >= 0x4006DB6D){p = pr3; q = ps3;} else if (ix >= 0x4006DB6D){p = pr3; q = ps3;}
else if (ix >= 0x40000000){p = pr2; q = ps2;} else if (ix >= 0x40000000){p = pr2; q = ps2;}
z = one/(x*x); z = 1.0/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); s = 1.0+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
return one+ r/s; return 1.0+ r/s;
} }
/* For x >= 8, the asymptotic expansions of qone is /* For x >= 8, the asymptotic expansions of qone is
@ -378,8 +375,8 @@ static double qone(double x)
else if (ix >= 0x40122E8B){p = qr5; q = qs5;} else if (ix >= 0x40122E8B){p = qr5; q = qs5;}
else if (ix >= 0x4006DB6D){p = qr3; q = qs3;} else if (ix >= 0x4006DB6D){p = qr3; q = qs3;}
else if (ix >= 0x40000000){p = qr2; q = qs2;} else if (ix >= 0x40000000){p = qr2; q = qs2;}
z = one/(x*x); z = 1.0/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5]))))); s = 1.0+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
return (.375 + r/s)/x; return (.375 + r/s)/x;
} }

View File

@ -19,7 +19,6 @@ static float ponef(float), qonef(float);
static const float static const float
huge = 1e30, huge = 1e30,
one = 1.0,
invsqrtpi = 5.6418961287e-01, /* 0x3f106ebb */ invsqrtpi = 5.6418961287e-01, /* 0x3f106ebb */
tpi = 6.3661974669e-01, /* 0x3f22f983 */ tpi = 6.3661974669e-01, /* 0x3f22f983 */
/* R0/S0 on [0,2] */ /* R0/S0 on [0,2] */
@ -33,8 +32,6 @@ s03 = 1.1771846857e-06, /* 0x359dffc2 */
s04 = 5.0463624390e-09, /* 0x31ad6446 */ s04 = 5.0463624390e-09, /* 0x31ad6446 */
s05 = 1.2354227016e-11; /* 0x2d59567e */ s05 = 1.2354227016e-11; /* 0x2d59567e */
static const float zero = 0.0;
float j1f(float x) float j1f(float x)
{ {
float z,s,c,ss,cc,r,u,v,y; float z,s,c,ss,cc,r,u,v,y;
@ -43,7 +40,7 @@ float j1f(float x)
GET_FLOAT_WORD(hx, x); GET_FLOAT_WORD(hx, x);
ix = hx & 0x7fffffff; ix = hx & 0x7fffffff;
if (ix >= 0x7f800000) if (ix >= 0x7f800000)
return one/x; return 1.0f/x;
y = fabsf(x); y = fabsf(x);
if (ix >= 0x40000000) { /* |x| >= 2.0 */ if (ix >= 0x40000000) { /* |x| >= 2.0 */
s = sinf(y); s = sinf(y);
@ -52,7 +49,7 @@ float j1f(float x)
cc = s-c; cc = s-c;
if (ix < 0x7f000000) { /* make sure y+y not overflow */ if (ix < 0x7f000000) { /* make sure y+y not overflow */
z = cosf(y+y); z = cosf(y+y);
if (s*c > zero) if (s*c > 0.0f)
cc = z/ss; cc = z/ss;
else else
ss = z/cc; ss = z/cc;
@ -74,12 +71,12 @@ float j1f(float x)
} }
if (ix < 0x32000000) { /* |x| < 2**-27 */ if (ix < 0x32000000) { /* |x| < 2**-27 */
/* raise inexact if x!=0 */ /* raise inexact if x!=0 */
if (huge+x > one) if (huge+x > 1.0f)
return 0.5f*x; return 0.5f*x;
} }
z = x*x; z = x*x;
r = z*(r00+z*(r01+z*(r02+z*r03))); r = z*(r00+z*(r01+z*(r02+z*r03)));
s = one+z*(s01+z*(s02+z*(s03+z*(s04+z*s05)))); s = 1.0f+z*(s01+z*(s02+z*(s03+z*(s04+z*s05))));
r *= x; r *= x;
return 0.5f*x + r/s; return 0.5f*x + r/s;
} }
@ -108,11 +105,11 @@ float y1f(float x)
ix = 0x7fffffff & hx; ix = 0x7fffffff & hx;
/* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */ /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
if (ix >= 0x7f800000) if (ix >= 0x7f800000)
return one/(x+x*x); return 1.0f/(x+x*x);
if (ix == 0) if (ix == 0)
return -one/zero; return -1.0f/0.0f;
if (hx < 0) if (hx < 0)
return zero/zero; return 0.0f/0.0f;
if (ix >= 0x40000000) { /* |x| >= 2.0 */ if (ix >= 0x40000000) { /* |x| >= 2.0 */
s = sinf(x); s = sinf(x);
c = cosf(x); c = cosf(x);
@ -120,7 +117,7 @@ float y1f(float x)
cc = s-c; cc = s-c;
if (ix < 0x7f000000) { /* make sure x+x not overflow */ if (ix < 0x7f000000) { /* make sure x+x not overflow */
z = cosf(x+x); z = cosf(x+x);
if (s*c > zero) if (s*c > 0.0f)
cc = z/ss; cc = z/ss;
else else
ss = z/cc; ss = z/cc;
@ -149,8 +146,8 @@ float y1f(float x)
return -tpi/x; return -tpi/x;
z = x*x; z = x*x;
u = U0[0]+z*(U0[1]+z*(U0[2]+z*(U0[3]+z*U0[4]))); u = U0[0]+z*(U0[1]+z*(U0[2]+z*(U0[3]+z*U0[4])));
v = one+z*(V0[0]+z*(V0[1]+z*(V0[2]+z*(V0[3]+z*V0[4])))); v = 1.0f+z*(V0[0]+z*(V0[1]+z*(V0[2]+z*(V0[3]+z*V0[4]))));
return x*(u/v) + tpi*(j1f(x)*logf(x)-one/x); return x*(u/v) + tpi*(j1f(x)*logf(x)-1.0f/x);
} }
/* For x >= 8, the asymptotic expansions of pone is /* For x >= 8, the asymptotic expansions of pone is
@ -239,10 +236,10 @@ static float ponef(float x)
else if (ix >= 0x40f71c58){p = pr5; q = ps5;} else if (ix >= 0x40f71c58){p = pr5; q = ps5;}
else if (ix >= 0x4036db68){p = pr3; q = ps3;} else if (ix >= 0x4036db68){p = pr3; q = ps3;}
else if (ix >= 0x40000000){p = pr2; q = ps2;} else if (ix >= 0x40000000){p = pr2; q = ps2;}
z = one/(x*x); z = 1.0f/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); s = 1.0f+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
return one + r/s; return 1.0f + r/s;
} }
/* For x >= 8, the asymptotic expansions of qone is /* For x >= 8, the asymptotic expansions of qone is
@ -335,8 +332,8 @@ static float qonef(float x)
else if (ix >= 0x40f71c58){p = qr5; q = qs5;} else if (ix >= 0x40f71c58){p = qr5; q = qs5;}
else if (ix >= 0x4036db68){p = qr3; q = qs3;} else if (ix >= 0x4036db68){p = qr3; q = qs3;}
else if (ix >= 0x40000000){p = qr2; q = qs2;} else if (ix >= 0x40000000){p = qr2; q = qs2;}
z = one/(x*x); z = 1.0f/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5]))))); s = 1.0f+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
return (.375f + r/s)/x; return (.375f + r/s)/x;
} }

View File

@ -37,12 +37,7 @@
#include "libm.h" #include "libm.h"
static const double static const double invsqrtpi = 5.64189583547756279280e-01; /* 0x3FE20DD7, 0x50429B6D */
invsqrtpi = 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
two = 2.00000000000000000000e+00, /* 0x40000000, 0x00000000 */
one = 1.00000000000000000000e+00; /* 0x3FF00000, 0x00000000 */
static const double zero = 0.00000000000000000000e+00;
double jn(int n, double x) double jn(int n, double x)
{ {
@ -69,7 +64,7 @@ double jn(int n, double x)
sgn = (n&1)&(hx>>31); /* even n -- 0, odd n -- sign(x) */ sgn = (n&1)&(hx>>31); /* even n -- 0, odd n -- sign(x) */
x = fabs(x); x = fabs(x);
if ((ix|lx) == 0 || ix >= 0x7ff00000) /* if x is 0 or inf */ if ((ix|lx) == 0 || ix >= 0x7ff00000) /* if x is 0 or inf */
b = zero; b = 0.0;
else if ((double)n <= x) { else if ((double)n <= x) {
/* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */ /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
if (ix >= 0x52D00000) { /* x > 2**302 */ if (ix >= 0x52D00000) { /* x > 2**302 */
@ -108,11 +103,11 @@ double jn(int n, double x)
* J(n,x) = 1/n!*(x/2)^n - ... * J(n,x) = 1/n!*(x/2)^n - ...
*/ */
if (n > 33) /* underflow */ if (n > 33) /* underflow */
b = zero; b = 0.0;
else { else {
temp = x*0.5; temp = x*0.5;
b = temp; b = temp;
for (a=one,i=2; i<=n; i++) { for (a=1.0,i=2; i<=n; i++) {
a *= (double)i; /* a = n! */ a *= (double)i; /* a = n! */
b *= temp; /* b = (x/2)^n */ b *= temp; /* b = (x/2)^n */
} }
@ -165,10 +160,10 @@ double jn(int n, double x)
q1 = tmp; q1 = tmp;
} }
m = n+n; m = n+n;
for (t=zero, i = 2*(n+k); i>=m; i -= 2) for (t=0.0, i = 2*(n+k); i>=m; i -= 2)
t = one/(i/x-t); t = 1.0/(i/x-t);
a = t; a = t;
b = one; b = 1.0;
/* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n) /* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
* Hence, if n*(log(2n/x)) > ... * Hence, if n*(log(2n/x)) > ...
* single 8.8722839355e+01 * single 8.8722839355e+01
@ -178,7 +173,7 @@ double jn(int n, double x)
* likely underflow to zero * likely underflow to zero
*/ */
tmp = n; tmp = n;
v = two/x; v = 2.0/x;
tmp = tmp*log(fabs(v*tmp)); tmp = tmp*log(fabs(v*tmp));
if (tmp < 7.09782712893383973096e+02) { if (tmp < 7.09782712893383973096e+02) {
for (i=n-1,di=(double)(i+i); i>0; i--) { for (i=n-1,di=(double)(i+i); i>0; i--) {
@ -186,7 +181,7 @@ double jn(int n, double x)
b *= di; b *= di;
b = b/x - a; b = b/x - a;
a = temp; a = temp;
di -= two; di -= 2.0;
} }
} else { } else {
for (i=n-1,di=(double)(i+i); i>0; i--) { for (i=n-1,di=(double)(i+i); i>0; i--) {
@ -194,12 +189,12 @@ double jn(int n, double x)
b *= di; b *= di;
b = b/x - a; b = b/x - a;
a = temp; a = temp;
di -= two; di -= 2.0;
/* scale b to avoid spurious overflow */ /* scale b to avoid spurious overflow */
if (b > 1e100) { if (b > 1e100) {
a /= b; a /= b;
t /= b; t /= b;
b = one; b = 1.0;
} }
} }
} }
@ -229,9 +224,9 @@ double yn(int n, double x)
if ((ix|((uint32_t)(lx|-lx))>>31) > 0x7ff00000) if ((ix|((uint32_t)(lx|-lx))>>31) > 0x7ff00000)
return x+x; return x+x;
if ((ix|lx) == 0) if ((ix|lx) == 0)
return -one/zero; return -1.0/0.0;
if (hx < 0) if (hx < 0)
return zero/zero; return 0.0/0.0;
sign = 1; sign = 1;
if (n < 0) { if (n < 0) {
n = -n; n = -n;
@ -242,7 +237,7 @@ double yn(int n, double x)
if (n == 1) if (n == 1)
return sign*y1(x); return sign*y1(x);
if (ix == 0x7ff00000) if (ix == 0x7ff00000)
return zero; return 0.0;
if (ix >= 0x52D00000) { /* x > 2**302 */ if (ix >= 0x52D00000) { /* x > 2**302 */
/* (x >> n**2) /* (x >> n**2)
* Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi) * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)

View File

@ -16,12 +16,6 @@
#define _GNU_SOURCE #define _GNU_SOURCE
#include "libm.h" #include "libm.h"
static const float
two = 2.0000000000e+00, /* 0x40000000 */
one = 1.0000000000e+00; /* 0x3F800000 */
static const float zero = 0.0000000000e+00;
float jnf(int n, float x) float jnf(int n, float x)
{ {
int32_t i,hx,ix, sgn; int32_t i,hx,ix, sgn;
@ -47,7 +41,7 @@ float jnf(int n, float x)
sgn = (n&1)&(hx>>31); /* even n -- 0, odd n -- sign(x) */ sgn = (n&1)&(hx>>31); /* even n -- 0, odd n -- sign(x) */
x = fabsf(x); x = fabsf(x);
if (ix == 0 || ix >= 0x7f800000) /* if x is 0 or inf */ if (ix == 0 || ix >= 0x7f800000) /* if x is 0 or inf */
b = zero; b = 0.0f;
else if((float)n <= x) { else if((float)n <= x) {
/* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */ /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
a = j0f(x); a = j0f(x);
@ -63,11 +57,11 @@ float jnf(int n, float x)
* J(n,x) = 1/n!*(x/2)^n - ... * J(n,x) = 1/n!*(x/2)^n - ...
*/ */
if (n > 33) /* underflow */ if (n > 33) /* underflow */
b = zero; b = 0.0f;
else { else {
temp = 0.5f * x; temp = 0.5f * x;
b = temp; b = temp;
for (a=one,i=2; i<=n; i++) { for (a=1.0f,i=2; i<=n; i++) {
a *= (float)i; /* a = n! */ a *= (float)i; /* a = n! */
b *= temp; /* b = (x/2)^n */ b *= temp; /* b = (x/2)^n */
} }
@ -121,10 +115,10 @@ float jnf(int n, float x)
q1 = tmp; q1 = tmp;
} }
m = n+n; m = n+n;
for (t=zero, i = 2*(n+k); i>=m; i -= 2) for (t=0.0f, i = 2*(n+k); i>=m; i -= 2)
t = one/(i/x-t); t = 1.0f/(i/x-t);
a = t; a = t;
b = one; b = 1.0f;
/* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n) /* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
* Hence, if n*(log(2n/x)) > ... * Hence, if n*(log(2n/x)) > ...
* single 8.8722839355e+01 * single 8.8722839355e+01
@ -134,7 +128,7 @@ float jnf(int n, float x)
* likely underflow to zero * likely underflow to zero
*/ */
tmp = n; tmp = n;
v = two/x; v = 2.0f/x;
tmp = tmp*logf(fabsf(v*tmp)); tmp = tmp*logf(fabsf(v*tmp));
if (tmp < 88.721679688f) { if (tmp < 88.721679688f) {
for (i=n-1,di=(float)(i+i); i>0; i--) { for (i=n-1,di=(float)(i+i); i>0; i--) {
@ -142,7 +136,7 @@ float jnf(int n, float x)
b *= di; b *= di;
b = b/x - a; b = b/x - a;
a = temp; a = temp;
di -= two; di -= 2.0f;
} }
} else { } else {
for (i=n-1,di=(float)(i+i); i>0; i--){ for (i=n-1,di=(float)(i+i); i>0; i--){
@ -150,12 +144,12 @@ float jnf(int n, float x)
b *= di; b *= di;
b = b/x - a; b = b/x - a;
a = temp; a = temp;
di -= two; di -= 2.0f;
/* scale b to avoid spurious overflow */ /* scale b to avoid spurious overflow */
if (b > 1e10f) { if (b > 1e10f) {
a /= b; a /= b;
t /= b; t /= b;
b = one; b = 1.0f;
} }
} }
} }
@ -183,9 +177,9 @@ float ynf(int n, float x)
if (ix > 0x7f800000) if (ix > 0x7f800000)
return x+x; return x+x;
if (ix == 0) if (ix == 0)
return -one/zero; return -1.0f/0.0f;
if (hx < 0) if (hx < 0)
return zero/zero; return 0.0f/0.0f;
sign = 1; sign = 1;
if (n < 0) { if (n < 0) {
n = -n; n = -n;
@ -196,7 +190,7 @@ float ynf(int n, float x)
if (n == 1) if (n == 1)
return sign*y1f(x); return sign*y1f(x);
if (ix == 0x7f800000) if (ix == 0x7f800000)
return zero; return 0.0f;
a = y0f(x); a = y0f(x);
b = y1f(x); b = y1f(x);

View File

@ -82,8 +82,6 @@
static const double static const double
two52= 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */ two52= 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
half= 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
pi = 3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */ pi = 3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */
a0 = 7.72156649015328655494e-02, /* 0x3FB3C467, 0xE37DB0C8 */ a0 = 7.72156649015328655494e-02, /* 0x3FB3C467, 0xE37DB0C8 */
a1 = 3.22467033424113591611e-01, /* 0x3FD4A34C, 0xC4A60FAD */ a1 = 3.22467033424113591611e-01, /* 0x3FD4A34C, 0xC4A60FAD */
@ -148,8 +146,6 @@ w4 = -5.95187557450339963135e-04, /* 0xBF4380CB, 0x8C0FE741 */
w5 = 8.36339918996282139126e-04, /* 0x3F4B67BA, 0x4CDAD5D1 */ w5 = 8.36339918996282139126e-04, /* 0x3F4B67BA, 0x4CDAD5D1 */
w6 = -1.63092934096575273989e-03; /* 0xBF5AB89D, 0x0B9E43E4 */ w6 = -1.63092934096575273989e-03; /* 0xBF5AB89D, 0x0B9E43E4 */
static const double zero = 0.00000000000000000000e+00;
static double sin_pi(double x) static double sin_pi(double x)
{ {
double y,z; double y,z;
@ -159,7 +155,7 @@ static double sin_pi(double x)
ix &= 0x7fffffff; ix &= 0x7fffffff;
if (ix < 0x3fd00000) if (ix < 0x3fd00000)
return __sin(pi*x, zero, 0); return __sin(pi*x, 0.0, 0);
y = -x; /* negative x is assumed */ y = -x; /* negative x is assumed */
@ -174,7 +170,7 @@ static double sin_pi(double x)
n = (int)(y*4.0); n = (int)(y*4.0);
} else { } else {
if (ix >= 0x43400000) { if (ix >= 0x43400000) {
y = zero; /* y must be even */ y = 0.0; /* y must be even */
n = 0; n = 0;
} else { } else {
if (ix < 0x43300000) if (ix < 0x43300000)
@ -186,14 +182,14 @@ static double sin_pi(double x)
} }
} }
switch (n) { switch (n) {
case 0: y = __sin(pi*y, zero, 0); break; case 0: y = __sin(pi*y, 0.0, 0); break;
case 1: case 1:
case 2: y = __cos(pi*(0.5-y), zero); break; case 2: y = __cos(pi*(0.5-y), 0.0); break;
case 3: case 3:
case 4: y = __sin(pi*(one-y), zero, 0); break; case 4: y = __sin(pi*(1.0-y), 0.0, 0); break;
case 5: case 5:
case 6: y = -__cos(pi*(y-1.5), zero); break; case 6: y = -__cos(pi*(y-1.5), 0.0); break;
default: y = __sin(pi*(y-2.0), zero, 0); break; default: y = __sin(pi*(y-2.0), 0.0, 0); break;
} }
return -y; return -y;
} }
@ -213,7 +209,7 @@ double __lgamma_r(double x, int *signgamp)
if (ix >= 0x7ff00000) if (ix >= 0x7ff00000)
return x*x; return x*x;
if ((ix|lx) == 0) if ((ix|lx) == 0)
return one/zero; return 1.0/0.0;
if (ix < 0x3b900000) { /* |x|<2**-70, return -log(|x|) */ if (ix < 0x3b900000) { /* |x|<2**-70, return -log(|x|) */
if(hx < 0) { if(hx < 0) {
*signgamp = -1; *signgamp = -1;
@ -223,12 +219,12 @@ double __lgamma_r(double x, int *signgamp)
} }
if (hx < 0) { if (hx < 0) {
if (ix >= 0x43300000) /* |x|>=2**52, must be -integer */ if (ix >= 0x43300000) /* |x|>=2**52, must be -integer */
return one/zero; return 1.0/0.0;
t = sin_pi(x); t = sin_pi(x);
if (t == zero) /* -integer */ if (t == 0.0) /* -integer */
return one/zero; return 1.0/0.0;
nadj = log(pi/fabs(t*x)); nadj = log(pi/fabs(t*x));
if (t < zero) if (t < 0.0)
*signgamp = -1; *signgamp = -1;
x = -x; x = -x;
} }
@ -241,17 +237,17 @@ double __lgamma_r(double x, int *signgamp)
if (ix <= 0x3feccccc) { /* lgamma(x) = lgamma(x+1)-log(x) */ if (ix <= 0x3feccccc) { /* lgamma(x) = lgamma(x+1)-log(x) */
r = -log(x); r = -log(x);
if (ix >= 0x3FE76944) { if (ix >= 0x3FE76944) {
y = one - x; y = 1.0 - x;
i = 0; i = 0;
} else if (ix >= 0x3FCDA661) { } else if (ix >= 0x3FCDA661) {
y = x - (tc-one); y = x - (tc-1.0);
i = 1; i = 1;
} else { } else {
y = x; y = x;
i = 2; i = 2;
} }
} else { } else {
r = zero; r = 0.0;
if (ix >= 0x3FFBB4C3) { /* [1.7316,2] */ if (ix >= 0x3FFBB4C3) { /* [1.7316,2] */
y = 2.0 - x; y = 2.0 - x;
i = 0; i = 0;
@ -259,7 +255,7 @@ double __lgamma_r(double x, int *signgamp)
y = x - tc; y = x - tc;
i = 1; i = 1;
} else { } else {
y = x - one; y = x - 1.0;
i = 2; i = 2;
} }
} }
@ -282,16 +278,16 @@ double __lgamma_r(double x, int *signgamp)
break; break;
case 2: case 2:
p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5))))); p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5)))));
p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*v5)))); p2 = 1.0+y*(v1+y*(v2+y*(v3+y*(v4+y*v5))));
r += -0.5*y + p1/p2; r += -0.5*y + p1/p2;
} }
} else if (ix < 0x40200000) { /* x < 8.0 */ } else if (ix < 0x40200000) { /* x < 8.0 */
i = (int)x; i = (int)x;
y = x - (double)i; y = x - (double)i;
p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6)))))); p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6))))));
q = one+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6))))); q = 1.0+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))));
r = half*y+p/q; r = 0.5*y+p/q;
z = one; /* lgamma(1+s) = log(s) + lgamma(s) */ z = 1.0; /* lgamma(1+s) = log(s) + lgamma(s) */
switch (i) { switch (i) {
case 7: z *= y + 6.0; /* FALLTHRU */ case 7: z *= y + 6.0; /* FALLTHRU */
case 6: z *= y + 5.0; /* FALLTHRU */ case 6: z *= y + 5.0; /* FALLTHRU */
@ -303,12 +299,12 @@ double __lgamma_r(double x, int *signgamp)
} }
} else if (ix < 0x43900000) { /* 8.0 <= x < 2**58 */ } else if (ix < 0x43900000) { /* 8.0 <= x < 2**58 */
t = log(x); t = log(x);
z = one/x; z = 1.0/x;
y = z*z; y = z*z;
w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6))))); w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6)))));
r = (x-half)*(t-one)+w; r = (x-0.5)*(t-1.0)+w;
} else /* 2**58 <= x <= inf */ } else /* 2**58 <= x <= inf */
r = x*(log(x)-one); r = x*(log(x)-1.0);
if (hx < 0) if (hx < 0)
r = nadj - r; r = nadj - r;
return r; return r;

View File

@ -17,8 +17,6 @@
static const float static const float
two23= 8.3886080000e+06, /* 0x4b000000 */ two23= 8.3886080000e+06, /* 0x4b000000 */
half= 5.0000000000e-01, /* 0x3f000000 */
one = 1.0000000000e+00, /* 0x3f800000 */
pi = 3.1415927410e+00, /* 0x40490fdb */ pi = 3.1415927410e+00, /* 0x40490fdb */
a0 = 7.7215664089e-02, /* 0x3d9e233f */ a0 = 7.7215664089e-02, /* 0x3d9e233f */
a1 = 3.2246702909e-01, /* 0x3ea51a66 */ a1 = 3.2246702909e-01, /* 0x3ea51a66 */
@ -83,8 +81,6 @@ w4 = -5.9518753551e-04, /* 0xba1c065c */
w5 = 8.3633989561e-04, /* 0x3a5b3dd2 */ w5 = 8.3633989561e-04, /* 0x3a5b3dd2 */
w6 = -1.6309292987e-03; /* 0xbad5c4e8 */ w6 = -1.6309292987e-03; /* 0xbad5c4e8 */
static const float zero = 0.0000000000e+00;
static float sin_pif(float x) static float sin_pif(float x)
{ {
float y,z; float y,z;
@ -109,7 +105,7 @@ static float sin_pif(float x)
n = (int)(y*4.0f); n = (int)(y*4.0f);
} else { } else {
if (ix >= 0x4b800000) { if (ix >= 0x4b800000) {
y = zero; /* y must be even */ y = 0.0f; /* y must be even */
n = 0; n = 0;
} else { } else {
if (ix < 0x4b000000) if (ix < 0x4b000000)
@ -125,7 +121,7 @@ static float sin_pif(float x)
case 1: case 1:
case 2: y = __cosdf(pi*(0.5f - y)); break; case 2: y = __cosdf(pi*(0.5f - y)); break;
case 3: case 3:
case 4: y = __sindf(pi*(one - y)); break; case 4: y = __sindf(pi*(1.0f - y)); break;
case 5: case 5:
case 6: y = -__cosdf(pi*(y - 1.5f)); break; case 6: y = -__cosdf(pi*(y - 1.5f)); break;
default: y = __sindf(pi*(y - 2.0f)); break; default: y = __sindf(pi*(y - 2.0f)); break;
@ -148,7 +144,7 @@ float __lgammaf_r(float x, int *signgamp)
if (ix >= 0x7f800000) if (ix >= 0x7f800000)
return x*x; return x*x;
if (ix == 0) if (ix == 0)
return one/zero; return 1.0f/0.0f;
if (ix < 0x35000000) { /* |x| < 2**-21, return -log(|x|) */ if (ix < 0x35000000) { /* |x| < 2**-21, return -log(|x|) */
if (hx < 0) { if (hx < 0) {
*signgamp = -1; *signgamp = -1;
@ -158,12 +154,12 @@ float __lgammaf_r(float x, int *signgamp)
} }
if (hx < 0) { if (hx < 0) {
if (ix >= 0x4b000000) /* |x| >= 2**23, must be -integer */ if (ix >= 0x4b000000) /* |x| >= 2**23, must be -integer */
return one/zero; return 1.0f/0.0f;
t = sin_pif(x); t = sin_pif(x);
if (t == zero) /* -integer */ if (t == 0.0f) /* -integer */
return one/zero; return 1.0f/0.0f;
nadj = logf(pi/fabsf(t*x)); nadj = logf(pi/fabsf(t*x));
if (t < zero) if (t < 0.0f)
*signgamp = -1; *signgamp = -1;
x = -x; x = -x;
} }
@ -176,17 +172,17 @@ float __lgammaf_r(float x, int *signgamp)
if (ix <= 0x3f666666) { /* lgamma(x) = lgamma(x+1)-log(x) */ if (ix <= 0x3f666666) { /* lgamma(x) = lgamma(x+1)-log(x) */
r = -logf(x); r = -logf(x);
if (ix >= 0x3f3b4a20) { if (ix >= 0x3f3b4a20) {
y = one - x; y = 1.0f - x;
i = 0; i = 0;
} else if (ix >= 0x3e6d3308) { } else if (ix >= 0x3e6d3308) {
y = x - (tc-one); y = x - (tc-1.0f);
i = 1; i = 1;
} else { } else {
y = x; y = x;
i = 2; i = 2;
} }
} else { } else {
r = zero; r = 0.0f;
if (ix >= 0x3fdda618) { /* [1.7316,2] */ if (ix >= 0x3fdda618) { /* [1.7316,2] */
y = 2.0f - x; y = 2.0f - x;
i = 0; i = 0;
@ -194,7 +190,7 @@ float __lgammaf_r(float x, int *signgamp)
y = x - tc; y = x - tc;
i = 1; i = 1;
} else { } else {
y = x - one; y = x - 1.0f;
i = 2; i = 2;
} }
} }
@ -217,16 +213,16 @@ float __lgammaf_r(float x, int *signgamp)
break; break;
case 2: case 2:
p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5))))); p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5)))));
p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*v5)))); p2 = 1.0f+y*(v1+y*(v2+y*(v3+y*(v4+y*v5))));
r += -0.5f*y + p1/p2; r += -0.5f*y + p1/p2;
} }
} else if (ix < 0x41000000) { /* x < 8.0 */ } else if (ix < 0x41000000) { /* x < 8.0 */
i = (int)x; i = (int)x;
y = x - (float)i; y = x - (float)i;
p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6)))))); p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6))))));
q = one+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6))))); q = 1.0f+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))));
r = half*y+p/q; r = 0.5f*y+p/q;
z = one; /* lgamma(1+s) = log(s) + lgamma(s) */ z = 1.0f; /* lgamma(1+s) = log(s) + lgamma(s) */
switch (i) { switch (i) {
case 7: z *= y + 6.0f; /* FALLTHRU */ case 7: z *= y + 6.0f; /* FALLTHRU */
case 6: z *= y + 5.0f; /* FALLTHRU */ case 6: z *= y + 5.0f; /* FALLTHRU */
@ -238,12 +234,12 @@ float __lgammaf_r(float x, int *signgamp)
} }
} else if (ix < 0x5c800000) { /* 8.0 <= x < 2**58 */ } else if (ix < 0x5c800000) { /* 8.0 <= x < 2**58 */
t = logf(x); t = logf(x);
z = one/x; z = 1.0f/x;
y = z*z; y = z*z;
w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6))))); w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6)))));
r = (x-half)*(t-one)+w; r = (x-0.5f)*(t-1.0f)+w;
} else /* 2**58 <= x <= inf */ } else /* 2**58 <= x <= inf */
r = x*(logf(x)-one); r = x*(logf(x)-1.0f);
if (hx < 0) if (hx < 0)
r = nadj - r; r = nadj - r;
return r; return r;

View File

@ -95,8 +95,6 @@ long double __lgammal_r(long double x, int *sg)
} }
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
static const long double static const long double
half = 0.5L,
one = 1.0L,
pi = 3.14159265358979323846264L, pi = 3.14159265358979323846264L,
two63 = 9.223372036854775808e18L, two63 = 9.223372036854775808e18L,
@ -200,8 +198,6 @@ w5 = 8.412723297322498080632E-4L,
w6 = -1.880801938119376907179E-3L, w6 = -1.880801938119376907179E-3L,
w7 = 4.885026142432270781165E-3L; w7 = 4.885026142432270781165E-3L;
static const long double zero = 0.0L;
static long double sin_pi(long double x) static long double sin_pi(long double x)
{ {
long double y, z; long double y, z;
@ -226,7 +222,7 @@ static long double sin_pi(long double x)
n = (int) (y*4.0); n = (int) (y*4.0);
} else { } else {
if (ix >= 0x403f8000) { /* 2^64 */ if (ix >= 0x403f8000) { /* 2^64 */
y = zero; /* y must be even */ y = 0.0; /* y must be even */
n = 0; n = 0;
} else { } else {
if (ix < 0x403e8000) /* 2^63 */ if (ix < 0x403e8000) /* 2^63 */
@ -244,11 +240,11 @@ static long double sin_pi(long double x)
break; break;
case 1: case 1:
case 2: case 2:
y = cosl(pi * (half - y)); y = cosl(pi * (0.5 - y));
break; break;
case 3: case 3:
case 4: case 4:
y = sinl(pi * (one - y)); y = sinl(pi * (1.0 - y));
break; break;
case 5: case 5:
case 6: case 6:
@ -273,7 +269,7 @@ long double __lgammal_r(long double x, int *sg) {
if ((ix | i0 | i1) == 0) { if ((ix | i0 | i1) == 0) {
if (se & 0x8000) if (se & 0x8000)
*sg = -1; *sg = -1;
return one / fabsl(x); return 1.0 / fabsl(x);
} }
ix = (ix << 16) | (i0 >> 16); ix = (ix << 16) | (i0 >> 16);
@ -291,10 +287,10 @@ long double __lgammal_r(long double x, int *sg) {
} }
if (se & 0x8000) { if (se & 0x8000) {
t = sin_pi (x); t = sin_pi (x);
if (t == zero) if (t == 0.0)
return one / fabsl(t); /* -integer */ return 1.0 / fabsl(t); /* -integer */
nadj = logl(pi / fabsl(t * x)); nadj = logl(pi / fabsl(t * x));
if (t < zero) if (t < 0.0)
*sg = -1; *sg = -1;
x = -x; x = -x;
} }
@ -308,17 +304,17 @@ long double __lgammal_r(long double x, int *sg) {
/* lgamma(x) = lgamma(x+1) - log(x) */ /* lgamma(x) = lgamma(x+1) - log(x) */
r = -logl(x); r = -logl(x);
if (ix >= 0x3ffebb4a) { /* 7.31597900390625e-1 */ if (ix >= 0x3ffebb4a) { /* 7.31597900390625e-1 */
y = x - one; y = x - 1.0;
i = 0; i = 0;
} else if (ix >= 0x3ffced33) { /* 2.31639862060546875e-1 */ } else if (ix >= 0x3ffced33) { /* 2.31639862060546875e-1 */
y = x - (tc - one); y = x - (tc - 1.0);
i = 1; i = 1;
} else { /* x < 0.23 */ } else { /* x < 0.23 */
y = x; y = x;
i = 2; i = 2;
} }
} else { } else {
r = zero; r = 0.0;
if (ix >= 0x3fffdda6) { /* 1.73162841796875 */ if (ix >= 0x3fffdda6) { /* 1.73162841796875 */
/* [1.7316,2] */ /* [1.7316,2] */
y = x - 2.0; y = x - 2.0;
@ -329,7 +325,7 @@ long double __lgammal_r(long double x, int *sg) {
i = 1; i = 1;
} else { } else {
/* [0.9, 1.23] */ /* [0.9, 1.23] */
y = x - one; y = x - 1.0;
i = 2; i = 2;
} }
} }
@ -337,7 +333,7 @@ long double __lgammal_r(long double x, int *sg) {
case 0: case 0:
p1 = a0 + y * (a1 + y * (a2 + y * (a3 + y * (a4 + y * a5)))); p1 = a0 + y * (a1 + y * (a2 + y * (a3 + y * (a4 + y * a5))));
p2 = b0 + y * (b1 + y * (b2 + y * (b3 + y * (b4 + y)))); p2 = b0 + y * (b1 + y * (b2 + y * (b3 + y * (b4 + y))));
r += half * y + y * p1/p2; r += 0.5 * y + y * p1/p2;
break; break;
case 1: case 1:
p1 = g0 + y * (g1 + y * (g2 + y * (g3 + y * (g4 + y * (g5 + y * g6))))); p1 = g0 + y * (g1 + y * (g2 + y * (g3 + y * (g4 + y * (g5 + y * g6)))));
@ -348,17 +344,17 @@ long double __lgammal_r(long double x, int *sg) {
case 2: case 2:
p1 = y * (u0 + y * (u1 + y * (u2 + y * (u3 + y * (u4 + y * (u5 + y * u6)))))); p1 = y * (u0 + y * (u1 + y * (u2 + y * (u3 + y * (u4 + y * (u5 + y * u6))))));
p2 = v0 + y * (v1 + y * (v2 + y * (v3 + y * (v4 + y * (v5 + y))))); p2 = v0 + y * (v1 + y * (v2 + y * (v3 + y * (v4 + y * (v5 + y)))));
r += (-half * y + p1 / p2); r += (-0.5 * y + p1 / p2);
} }
} else if (ix < 0x40028000) { /* 8.0 */ } else if (ix < 0x40028000) { /* 8.0 */
/* x < 8.0 */ /* x < 8.0 */
i = (int)x; i = (int)x;
t = zero; t = 0.0;
y = x - (double)i; y = x - (double)i;
p = y * (s0 + y * (s1 + y * (s2 + y * (s3 + y * (s4 + y * (s5 + y * s6)))))); p = y * (s0 + y * (s1 + y * (s2 + y * (s3 + y * (s4 + y * (s5 + y * s6))))));
q = r0 + y * (r1 + y * (r2 + y * (r3 + y * (r4 + y * (r5 + y * (r6 + y)))))); q = r0 + y * (r1 + y * (r2 + y * (r3 + y * (r4 + y * (r5 + y * (r6 + y))))));
r = half * y + p / q; r = 0.5 * y + p / q;
z = one;/* lgamma(1+s) = log(s) + lgamma(s) */ z = 1.0;/* lgamma(1+s) = log(s) + lgamma(s) */
switch (i) { switch (i) {
case 7: case 7:
z *= (y + 6.0); /* FALLTHRU */ z *= (y + 6.0); /* FALLTHRU */
@ -376,12 +372,12 @@ long double __lgammal_r(long double x, int *sg) {
} else if (ix < 0x40418000) { /* 2^66 */ } else if (ix < 0x40418000) { /* 2^66 */
/* 8.0 <= x < 2**66 */ /* 8.0 <= x < 2**66 */
t = logl(x); t = logl(x);
z = one / x; z = 1.0 / x;
y = z * z; y = z * z;
w = w0 + z * (w1 + y * (w2 + y * (w3 + y * (w4 + y * (w5 + y * (w6 + y * w7)))))); w = w0 + z * (w1 + y * (w2 + y * (w3 + y * (w4 + y * (w5 + y * (w6 + y * w7))))));
r = (x - half) * (t - one) + w; r = (x - 0.5) * (t - 1.0) + w;
} else /* 2**66 <= x <= inf */ } else /* 2**66 <= x <= inf */
r = x * (logl (x) - one); r = x * (logl(x) - 1.0);
if (se & 0x8000) if (se & 0x8000)
r = nadj - r; r = nadj - r;
return r; return r;

View File

@ -74,8 +74,6 @@ Lg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */
Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */ Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */
Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */ Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
static const double zero = 0.0;
double log(double x) double log(double x)
{ {
double hfsq,f,s,z,R,w,t1,t2,dk; double hfsq,f,s,z,R,w,t1,t2,dk;
@ -87,9 +85,9 @@ double log(double x)
k = 0; k = 0;
if (hx < 0x00100000) { /* x < 2**-1022 */ if (hx < 0x00100000) { /* x < 2**-1022 */
if (((hx&0x7fffffff)|lx) == 0) if (((hx&0x7fffffff)|lx) == 0)
return -two54/zero; /* log(+-0)=-inf */ return -two54/0.0; /* log(+-0)=-inf */
if (hx < 0) if (hx < 0)
return (x-x)/zero; /* log(-#) = NaN */ return (x-x)/0.0; /* log(-#) = NaN */
/* subnormal number, scale up x */ /* subnormal number, scale up x */
k -= 54; k -= 54;
x *= two54; x *= two54;
@ -104,9 +102,9 @@ double log(double x)
k += i>>20; k += i>>20;
f = x - 1.0; f = x - 1.0;
if ((0x000fffff&(2+hx)) < 3) { /* -2**-20 <= f < 2**-20 */ if ((0x000fffff&(2+hx)) < 3) { /* -2**-20 <= f < 2**-20 */
if (f == zero) { if (f == 0.0) {
if (k == 0) { if (k == 0) {
return zero; return 0.0;
} }
dk = (double)k; dk = (double)k;
return dk*ln2_hi + dk*ln2_lo; return dk*ln2_hi + dk*ln2_lo;

View File

@ -27,8 +27,6 @@ ivln10lo = 2.50829467116452752298e-11, /* 0x3dbb9438, 0xca9aadd5 */
log10_2hi = 3.01029995663611771306e-01, /* 0x3FD34413, 0x509F6000 */ log10_2hi = 3.01029995663611771306e-01, /* 0x3FD34413, 0x509F6000 */
log10_2lo = 3.69423907715893078616e-13; /* 0x3D59FEF3, 0x11F12B36 */ log10_2lo = 3.69423907715893078616e-13; /* 0x3D59FEF3, 0x11F12B36 */
static const double zero = 0.0;
double log10(double x) double log10(double x)
{ {
double f,hfsq,hi,lo,r,val_hi,val_lo,w,y,y2; double f,hfsq,hi,lo,r,val_hi,val_lo,w,y,y2;
@ -40,9 +38,9 @@ double log10(double x)
k = 0; k = 0;
if (hx < 0x00100000) { /* x < 2**-1022 */ if (hx < 0x00100000) { /* x < 2**-1022 */
if (((hx&0x7fffffff)|lx) == 0) if (((hx&0x7fffffff)|lx) == 0)
return -two54/zero; /* log(+-0)=-inf */ return -two54/0.0; /* log(+-0)=-inf */
if (hx<0) if (hx<0)
return (x-x)/zero; /* log(-#) = NaN */ return (x-x)/0.0; /* log(-#) = NaN */
/* subnormal number, scale up x */ /* subnormal number, scale up x */
k -= 54; k -= 54;
x *= two54; x *= two54;
@ -51,7 +49,7 @@ double log10(double x)
if (hx >= 0x7ff00000) if (hx >= 0x7ff00000)
return x+x; return x+x;
if (hx == 0x3ff00000 && lx == 0) if (hx == 0x3ff00000 && lx == 0)
return zero; /* log(1) = +0 */ return 0.0; /* log(1) = +0 */
k += (hx>>20) - 1023; k += (hx>>20) - 1023;
hx &= 0x000fffff; hx &= 0x000fffff;
i = (hx+0x95f64)&0x100000; i = (hx+0x95f64)&0x100000;

View File

@ -23,8 +23,6 @@ ivln10lo = -3.1689971365e-05, /* 0xb804ead9 */
log10_2hi = 3.0102920532e-01, /* 0x3e9a2080 */ log10_2hi = 3.0102920532e-01, /* 0x3e9a2080 */
log10_2lo = 7.9034151668e-07; /* 0x355427db */ log10_2lo = 7.9034151668e-07; /* 0x355427db */
static const float zero = 0.0;
float log10f(float x) float log10f(float x)
{ {
float f,hfsq,hi,lo,r,y; float f,hfsq,hi,lo,r,y;
@ -35,9 +33,9 @@ float log10f(float x)
k = 0; k = 0;
if (hx < 0x00800000) { /* x < 2**-126 */ if (hx < 0x00800000) { /* x < 2**-126 */
if ((hx&0x7fffffff) == 0) if ((hx&0x7fffffff) == 0)
return -two25/zero; /* log(+-0)=-inf */ return -two25/0.0f; /* log(+-0)=-inf */
if (hx < 0) if (hx < 0)
return (x-x)/zero; /* log(-#) = NaN */ return (x-x)/0.0f; /* log(-#) = NaN */
/* subnormal number, scale up x */ /* subnormal number, scale up x */
k -= 25; k -= 25;
x *= two25; x *= two25;
@ -46,7 +44,7 @@ float log10f(float x)
if (hx >= 0x7f800000) if (hx >= 0x7f800000)
return x+x; return x+x;
if (hx == 0x3f800000) if (hx == 0x3f800000)
return zero; /* log(1) = +0 */ return 0.0f; /* log(1) = +0 */
k += (hx>>23) - 127; k += (hx>>23) - 127;
hx &= 0x007fffff; hx &= 0x007fffff;
i = (hx+(0x4afb0d))&0x800000; i = (hx+(0x4afb0d))&0x800000;

View File

@ -88,8 +88,6 @@ Lp5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */
Lp6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */ Lp6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */
Lp7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */ Lp7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
static const double zero = 0.0;
double log1p(double x) double log1p(double x)
{ {
double hfsq,f,c,s,z,R,u; double hfsq,f,c,s,z,R,u;
@ -102,12 +100,12 @@ double log1p(double x)
if (hx < 0x3FDA827A) { /* 1+x < sqrt(2)+ */ if (hx < 0x3FDA827A) { /* 1+x < sqrt(2)+ */
if (ax >= 0x3ff00000) { /* x <= -1.0 */ if (ax >= 0x3ff00000) { /* x <= -1.0 */
if (x == -1.0) if (x == -1.0)
return -two54/zero; /* log1p(-1)=+inf */ return -two54/0.0; /* log1p(-1)=+inf */
return (x-x)/(x-x); /* log1p(x<-1)=NaN */ return (x-x)/(x-x); /* log1p(x<-1)=NaN */
} }
if (ax < 0x3e200000) { /* |x| < 2**-29 */ if (ax < 0x3e200000) { /* |x| < 2**-29 */
/* raise inexact */ /* raise inexact */
if (two54 + x > zero && ax < 0x3c900000) /* |x| < 2**-54 */ if (two54 + x > 0.0 && ax < 0x3c900000) /* |x| < 2**-54 */
return x; return x;
return x - x*x*0.5; return x - x*x*0.5;
} }
@ -151,9 +149,9 @@ double log1p(double x)
} }
hfsq = 0.5*f*f; hfsq = 0.5*f*f;
if (hu == 0) { /* |f| < 2**-20 */ if (hu == 0) { /* |f| < 2**-20 */
if (f == zero) { if (f == 0.0) {
if(k == 0) if(k == 0)
return zero; return 0.0;
c += k*ln2_lo; c += k*ln2_lo;
return k*ln2_hi + c; return k*ln2_hi + c;
} }

View File

@ -27,8 +27,6 @@ Lp5 = 1.8183572590e-01, /* 3E3A3325 */
Lp6 = 1.5313838422e-01, /* 3E1CD04F */ Lp6 = 1.5313838422e-01, /* 3E1CD04F */
Lp7 = 1.4798198640e-01; /* 3E178897 */ Lp7 = 1.4798198640e-01; /* 3E178897 */
static const float zero = 0.0;
float log1pf(float x) float log1pf(float x)
{ {
float hfsq,f,c,s,z,R,u; float hfsq,f,c,s,z,R,u;
@ -41,12 +39,12 @@ float log1pf(float x)
if (hx < 0x3ed413d0) { /* 1+x < sqrt(2)+ */ if (hx < 0x3ed413d0) { /* 1+x < sqrt(2)+ */
if (ax >= 0x3f800000) { /* x <= -1.0 */ if (ax >= 0x3f800000) { /* x <= -1.0 */
if (x == -1.0f) if (x == -1.0f)
return -two25/zero; /* log1p(-1)=+inf */ return -two25/0.0f; /* log1p(-1)=+inf */
return (x-x)/(x-x); /* log1p(x<-1)=NaN */ return (x-x)/(x-x); /* log1p(x<-1)=NaN */
} }
if (ax < 0x38000000) { /* |x| < 2**-15 */ if (ax < 0x38000000) { /* |x| < 2**-15 */
/* raise inexact */ /* raise inexact */
if (two25 + x > zero && ax < 0x33800000) /* |x| < 2**-24 */ if (two25 + x > 0.0f && ax < 0x33800000) /* |x| < 2**-24 */
return x; return x;
return x - x*x*0.5f; return x - x*x*0.5f;
} }
@ -91,9 +89,9 @@ float log1pf(float x)
} }
hfsq = 0.5f * f * f; hfsq = 0.5f * f * f;
if (hu == 0) { /* |f| < 2**-20 */ if (hu == 0) { /* |f| < 2**-20 */
if (f == zero) { if (f == 0.0f) {
if (k == 0) if (k == 0)
return zero; return 0.0f;
c += k*ln2_lo; c += k*ln2_lo;
return k*ln2_hi+c; return k*ln2_hi+c;
} }

View File

@ -118,11 +118,11 @@ long double log1pl(long double xm1)
if (xm1 == 0.0) if (xm1 == 0.0)
return xm1; return xm1;
x = xm1 + 1.0L; x = xm1 + 1.0;
/* Test for domain errors. */ /* Test for domain errors. */
if (x <= 0.0L) { if (x <= 0.0) {
if (x == 0.0L) if (x == 0.0)
return -INFINITY; return -INFINITY;
return NAN; return NAN;
} }
@ -136,12 +136,12 @@ long double log1pl(long double xm1)
if (e > 2 || e < -2) { if (e > 2 || e < -2) {
if (x < SQRTH) { /* 2(2x-1)/(2x+1) */ if (x < SQRTH) { /* 2(2x-1)/(2x+1) */
e -= 1; e -= 1;
z = x - 0.5L; z = x - 0.5;
y = 0.5L * z + 0.5L; y = 0.5 * z + 0.5;
} else { /* 2 (x-1)/(x+1) */ } else { /* 2 (x-1)/(x+1) */
z = x - 0.5L; z = x - 0.5;
z -= 0.5L; z -= 0.5;
y = 0.5L * x + 0.5L; y = 0.5 * x + 0.5;
} }
x = z / y; x = z / y;
z = x*x; z = x*x;
@ -156,12 +156,12 @@ long double log1pl(long double xm1)
if (x < SQRTH) { if (x < SQRTH) {
e -= 1; e -= 1;
if (e != 0) if (e != 0)
x = 2.0 * x - 1.0L; x = 2.0 * x - 1.0;
else else
x = xm1; x = xm1;
} else { } else {
if (e != 0) if (e != 0)
x = x - 1.0L; x = x - 1.0;
else else
x = xm1; x = xm1;
} }

View File

@ -27,8 +27,6 @@ two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
ivln2hi = 1.44269504072144627571e+00, /* 0x3ff71547, 0x65200000 */ ivln2hi = 1.44269504072144627571e+00, /* 0x3ff71547, 0x65200000 */
ivln2lo = 1.67517131648865118353e-10; /* 0x3de705fc, 0x2eefa200 */ ivln2lo = 1.67517131648865118353e-10; /* 0x3de705fc, 0x2eefa200 */
static const double zero = 0.0;
double log2(double x) double log2(double x)
{ {
double f,hfsq,hi,lo,r,val_hi,val_lo,w,y; double f,hfsq,hi,lo,r,val_hi,val_lo,w,y;
@ -40,9 +38,9 @@ double log2(double x)
k = 0; k = 0;
if (hx < 0x00100000) { /* x < 2**-1022 */ if (hx < 0x00100000) { /* x < 2**-1022 */
if (((hx&0x7fffffff)|lx) == 0) if (((hx&0x7fffffff)|lx) == 0)
return -two54/zero; /* log(+-0)=-inf */ return -two54/0.0; /* log(+-0)=-inf */
if (hx < 0) if (hx < 0)
return (x-x)/zero; /* log(-#) = NaN */ return (x-x)/0.0; /* log(-#) = NaN */
/* subnormal number, scale up x */ /* subnormal number, scale up x */
k -= 54; k -= 54;
x *= two54; x *= two54;
@ -51,7 +49,7 @@ double log2(double x)
if (hx >= 0x7ff00000) if (hx >= 0x7ff00000)
return x+x; return x+x;
if (hx == 0x3ff00000 && lx == 0) if (hx == 0x3ff00000 && lx == 0)
return zero; /* log(1) = +0 */ return 0.0; /* log(1) = +0 */
k += (hx>>20) - 1023; k += (hx>>20) - 1023;
hx &= 0x000fffff; hx &= 0x000fffff;
i = (hx+0x95f64) & 0x100000; i = (hx+0x95f64) & 0x100000;

View File

@ -21,8 +21,6 @@ two25 = 3.3554432000e+07, /* 0x4c000000 */
ivln2hi = 1.4428710938e+00, /* 0x3fb8b000 */ ivln2hi = 1.4428710938e+00, /* 0x3fb8b000 */
ivln2lo = -1.7605285393e-04; /* 0xb9389ad4 */ ivln2lo = -1.7605285393e-04; /* 0xb9389ad4 */
static const float zero = 0.0;
float log2f(float x) float log2f(float x)
{ {
float f,hfsq,hi,lo,r,y; float f,hfsq,hi,lo,r,y;
@ -33,9 +31,9 @@ float log2f(float x)
k = 0; k = 0;
if (hx < 0x00800000) { /* x < 2**-126 */ if (hx < 0x00800000) { /* x < 2**-126 */
if ((hx&0x7fffffff) == 0) if ((hx&0x7fffffff) == 0)
return -two25/zero; /* log(+-0)=-inf */ return -two25/0.0f; /* log(+-0)=-inf */
if (hx < 0) if (hx < 0)
return (x-x)/zero; /* log(-#) = NaN */ return (x-x)/0.0f; /* log(-#) = NaN */
/* subnormal number, scale up x */ /* subnormal number, scale up x */
k -= 25; k -= 25;
x *= two25; x *= two25;
@ -44,7 +42,7 @@ float log2f(float x)
if (hx >= 0x7f800000) if (hx >= 0x7f800000)
return x+x; return x+x;
if (hx == 0x3f800000) if (hx == 0x3f800000)
return zero; /* log(1) = +0 */ return 0.0f; /* log(1) = +0 */
k += (hx>>23) - 127; k += (hx>>23) - 127;
hx &= 0x007fffff; hx &= 0x007fffff;
i = (hx+(0x4afb0d))&0x800000; i = (hx+(0x4afb0d))&0x800000;

View File

@ -25,8 +25,6 @@ Lg2 = 0xccce13.0p-25, /* 0.40000972152 */
Lg3 = 0x91e9ee.0p-25, /* 0.28498786688 */ Lg3 = 0x91e9ee.0p-25, /* 0.28498786688 */
Lg4 = 0xf89e26.0p-26; /* 0.24279078841 */ Lg4 = 0xf89e26.0p-26; /* 0.24279078841 */
static const float zero = 0.0;
float logf(float x) float logf(float x)
{ {
float hfsq,f,s,z,R,w,t1,t2,dk; float hfsq,f,s,z,R,w,t1,t2,dk;
@ -37,9 +35,9 @@ float logf(float x)
k = 0; k = 0;
if (ix < 0x00800000) { /* x < 2**-126 */ if (ix < 0x00800000) { /* x < 2**-126 */
if ((ix & 0x7fffffff) == 0) if ((ix & 0x7fffffff) == 0)
return -two25/zero; /* log(+-0)=-inf */ return -two25/0.0f; /* log(+-0)=-inf */
if (ix < 0) if (ix < 0)
return (x-x)/zero; /* log(-#) = NaN */ return (x-x)/0.0f; /* log(-#) = NaN */
/* subnormal number, scale up x */ /* subnormal number, scale up x */
k -= 25; k -= 25;
x *= two25; x *= two25;
@ -54,9 +52,9 @@ float logf(float x)
k += i>>23; k += i>>23;
f = x - 1.0f; f = x - 1.0f;
if ((0x007fffff & (0x8000 + ix)) < 0xc000) { /* -2**-9 <= f < 2**-9 */ if ((0x007fffff & (0x8000 + ix)) < 0xc000) { /* -2**-9 <= f < 2**-9 */
if (f == zero) { if (f == 0.0f) {
if (k == 0) if (k == 0)
return zero; return 0.0f;
dk = (float)k; dk = (float)k;
return dk*ln2_hi + dk*ln2_lo; return dk*ln2_hi + dk*ln2_lo;
} }

View File

@ -21,8 +21,6 @@
#include "libm.h" #include "libm.h"
static const double one = 1.0;
double modf(double x, double *iptr) double modf(double x, double *iptr)
{ {
int32_t i0,i1,j0; int32_t i0,i1,j0;
@ -51,7 +49,7 @@ double modf(double x, double *iptr)
*iptr = x; *iptr = x;
return 0.0 / x; return 0.0 / x;
} }
*iptr = x*one; *iptr = x;
GET_HIGH_WORD(high, x); GET_HIGH_WORD(high, x);
INSERT_WORDS(x, high & 0x80000000, 0); /* return +-0 */ INSERT_WORDS(x, high & 0x80000000, 0); /* return +-0 */
return x; return x;

View File

@ -59,9 +59,6 @@ static const double
bp[] = {1.0, 1.5,}, bp[] = {1.0, 1.5,},
dp_h[] = { 0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */ dp_h[] = { 0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */
dp_l[] = { 0.0, 1.35003920212974897128e-08,}, /* 0x3E4CFDEB, 0x43CFD006 */ dp_l[] = { 0.0, 1.35003920212974897128e-08,}, /* 0x3E4CFDEB, 0x43CFD006 */
zero = 0.0,
one = 1.0,
two = 2.0,
two53 = 9007199254740992.0, /* 0x43400000, 0x00000000 */ two53 = 9007199254740992.0, /* 0x43400000, 0x00000000 */
huge = 1.0e300, huge = 1.0e300,
tiny = 1.0e-300, tiny = 1.0e-300,
@ -101,15 +98,15 @@ double pow(double x, double y)
ix = hx & 0x7fffffff; ix = hx & 0x7fffffff;
iy = hy & 0x7fffffff; iy = hy & 0x7fffffff;
/* y == zero: x**0 = 1 */ /* y == 0.0: x**0 = 1 */
if ((iy|ly) == 0) if ((iy|ly) == 0)
return one; return 1.0;
/* x == 1: 1**y = 1, even if y is NaN */ /* x == 1: 1**y = 1, even if y is NaN */
if (hx == 0x3ff00000 && lx == 0) if (hx == 0x3ff00000 && lx == 0)
return one; return 1.0;
/* y != zero: result is NaN if either arg is NaN */ /* y != 0.0: result is NaN if either arg is NaN */
if (ix > 0x7ff00000 || (ix == 0x7ff00000 && lx != 0) || if (ix > 0x7ff00000 || (ix == 0x7ff00000 && lx != 0) ||
iy > 0x7ff00000 || (iy == 0x7ff00000 && ly != 0)) iy > 0x7ff00000 || (iy == 0x7ff00000 && ly != 0))
return (x+0.0) + (y+0.0); return (x+0.0) + (y+0.0);
@ -141,15 +138,15 @@ double pow(double x, double y)
if (ly == 0) { if (ly == 0) {
if (iy == 0x7ff00000) { /* y is +-inf */ if (iy == 0x7ff00000) { /* y is +-inf */
if (((ix-0x3ff00000)|lx) == 0) /* (-1)**+-inf is 1 */ if (((ix-0x3ff00000)|lx) == 0) /* (-1)**+-inf is 1 */
return one; return 1.0;
else if (ix >= 0x3ff00000) /* (|x|>1)**+-inf = inf,0 */ else if (ix >= 0x3ff00000) /* (|x|>1)**+-inf = inf,0 */
return hy >= 0 ? y : zero; return hy >= 0 ? y : 0.0;
else /* (|x|<1)**+-inf = 0,inf */ else /* (|x|<1)**+-inf = 0,inf */
return hy < 0 ? -y : zero; return hy < 0 ? -y : 0.0;
} }
if (iy == 0x3ff00000) { /* y is +-1 */ if (iy == 0x3ff00000) { /* y is +-1 */
if (hy < 0) if (hy < 0)
return one/x; return 1.0/x;
return x; return x;
} }
if (hy == 0x40000000) /* y is 2 */ if (hy == 0x40000000) /* y is 2 */
@ -166,7 +163,7 @@ double pow(double x, double y)
if (ix == 0x7ff00000 || ix == 0 || ix == 0x3ff00000) { /* x is +-0,+-inf,+-1 */ if (ix == 0x7ff00000 || ix == 0 || ix == 0x3ff00000) { /* x is +-0,+-inf,+-1 */
z = ax; z = ax;
if (hy < 0) /* z = (1/|x|) */ if (hy < 0) /* z = (1/|x|) */
z = one/z; z = 1.0/z;
if (hx < 0) { if (hx < 0) {
if (((ix-0x3ff00000)|yisint) == 0) { if (((ix-0x3ff00000)|yisint) == 0) {
z = (z-z)/(z-z); /* (-1)**non-int is NaN */ z = (z-z)/(z-z); /* (-1)**non-int is NaN */
@ -187,9 +184,9 @@ double pow(double x, double y)
if ((n|yisint) == 0) if ((n|yisint) == 0)
return (x-x)/(x-x); return (x-x)/(x-x);
s = one; /* s (sign of result -ve**odd) = -1 else = 1 */ s = 1.0; /* s (sign of result -ve**odd) = -1 else = 1 */
if ((n|(yisint-1)) == 0) if ((n|(yisint-1)) == 0)
s = -one;/* (-ve)**(odd int) */ s = -1.0;/* (-ve)**(odd int) */
/* |y| is huge */ /* |y| is huge */
if (iy > 0x41e00000) { /* if |y| > 2**31 */ if (iy > 0x41e00000) { /* if |y| > 2**31 */
@ -206,7 +203,7 @@ double pow(double x, double y)
return hy > 0 ? s*huge*huge : s*tiny*tiny; return hy > 0 ? s*huge*huge : s*tiny*tiny;
/* now |1-x| is tiny <= 2**-20, suffice to compute /* now |1-x| is tiny <= 2**-20, suffice to compute
log(x) by x-x^2/2+x^3/3-x^4/4 */ log(x) by x-x^2/2+x^3/3-x^4/4 */
t = ax - one; /* t has 20 trailing zeros */ t = ax - 1.0; /* t has 20 trailing zeros */
w = (t*t)*(0.5 - t*(0.3333333333333333333333-t*0.25)); w = (t*t)*(0.5 - t*(0.3333333333333333333333-t*0.25));
u = ivln2_h*t; /* ivln2_h has 21 sig. bits */ u = ivln2_h*t; /* ivln2_h has 21 sig. bits */
v = t*ivln2_l - w*ivln2; v = t*ivln2_l - w*ivln2;
@ -239,12 +236,12 @@ double pow(double x, double y)
/* compute ss = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */ /* compute ss = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
u = ax - bp[k]; /* bp[0]=1.0, bp[1]=1.5 */ u = ax - bp[k]; /* bp[0]=1.0, bp[1]=1.5 */
v = one/(ax+bp[k]); v = 1.0/(ax+bp[k]);
ss = u*v; ss = u*v;
s_h = ss; s_h = ss;
SET_LOW_WORD(s_h, 0); SET_LOW_WORD(s_h, 0);
/* t_h=ax+bp[k] High */ /* t_h=ax+bp[k] High */
t_h = zero; t_h = 0.0;
SET_HIGH_WORD(t_h, ((ix>>1)|0x20000000) + 0x00080000 + (k<<18)); SET_HIGH_WORD(t_h, ((ix>>1)|0x20000000) + 0x00080000 + (k<<18));
t_l = ax - (t_h-bp[k]); t_l = ax - (t_h-bp[k]);
s_l = v*((u-s_h*t_h)-s_h*t_l); s_l = v*((u-s_h*t_h)-s_h*t_l);
@ -299,7 +296,7 @@ double pow(double x, double y)
if (i > 0x3fe00000) { /* if |z| > 0.5, set n = [z+0.5] */ if (i > 0x3fe00000) { /* if |z| > 0.5, set n = [z+0.5] */
n = j + (0x00100000>>(k+1)); n = j + (0x00100000>>(k+1));
k = ((n&0x7fffffff)>>20) - 0x3ff; /* new k for n */ k = ((n&0x7fffffff)>>20) - 0x3ff; /* new k for n */
t = zero; t = 0.0;
SET_HIGH_WORD(t, n & ~(0x000fffff>>k)); SET_HIGH_WORD(t, n & ~(0x000fffff>>k));
n = ((n&0x000fffff)|0x00100000)>>(20-k); n = ((n&0x000fffff)|0x00100000)>>(20-k);
if (j < 0) if (j < 0)
@ -314,8 +311,8 @@ double pow(double x, double y)
w = v - (z-u); w = v - (z-u);
t = z*z; t = z*z;
t1 = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5)))); t1 = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
r = (z*t1)/(t1-two) - (w + z*w); r = (z*t1)/(t1-2.0) - (w + z*w);
z = one - (r-z); z = 1.0 - (r-z);
GET_HIGH_WORD(j, z); GET_HIGH_WORD(j, z);
j += n<<20; j += n<<20;
if ((j>>20) <= 0) /* subnormal output */ if ((j>>20) <= 0) /* subnormal output */

View File

@ -19,9 +19,6 @@ static const float
bp[] = {1.0, 1.5,}, bp[] = {1.0, 1.5,},
dp_h[] = { 0.0, 5.84960938e-01,}, /* 0x3f15c000 */ dp_h[] = { 0.0, 5.84960938e-01,}, /* 0x3f15c000 */
dp_l[] = { 0.0, 1.56322085e-06,}, /* 0x35d1cfdc */ dp_l[] = { 0.0, 1.56322085e-06,}, /* 0x35d1cfdc */
zero = 0.0,
one = 1.0,
two = 2.0,
two24 = 16777216.0, /* 0x4b800000 */ two24 = 16777216.0, /* 0x4b800000 */
huge = 1.0e30, huge = 1.0e30,
tiny = 1.0e-30, tiny = 1.0e-30,
@ -60,15 +57,15 @@ float powf(float x, float y)
ix = hx & 0x7fffffff; ix = hx & 0x7fffffff;
iy = hy & 0x7fffffff; iy = hy & 0x7fffffff;
/* y == zero: x**0 = 1 */ /* y == 0: x**0 = 1 */
if (iy == 0) if (iy == 0)
return one; return 1.0f;
/* x == 1: 1**y = 1, even if y is NaN */ /* x == 1: 1**y = 1, even if y is NaN */
if (hx == 0x3f800000) if (hx == 0x3f800000)
return one; return 1.0f;
/* y != zero: result is NaN if either arg is NaN */ /* y != 0: result is NaN if either arg is NaN */
if (ix > 0x7f800000 || iy > 0x7f800000) if (ix > 0x7f800000 || iy > 0x7f800000)
return (x+0.0f) + (y+0.0f); return (x+0.0f) + (y+0.0f);
@ -92,15 +89,15 @@ float powf(float x, float y)
/* special value of y */ /* special value of y */
if (iy == 0x7f800000) { /* y is +-inf */ if (iy == 0x7f800000) { /* y is +-inf */
if (ix == 0x3f800000) /* (-1)**+-inf is 1 */ if (ix == 0x3f800000) /* (-1)**+-inf is 1 */
return one; return 1.0f;
else if (ix > 0x3f800000) /* (|x|>1)**+-inf = inf,0 */ else if (ix > 0x3f800000) /* (|x|>1)**+-inf = inf,0 */
return hy >= 0 ? y : zero; return hy >= 0 ? y : 0.0f;
else /* (|x|<1)**+-inf = 0,inf */ else /* (|x|<1)**+-inf = 0,inf */
return hy < 0 ? -y : zero; return hy < 0 ? -y : 0.0f;
} }
if (iy == 0x3f800000) { /* y is +-1 */ if (iy == 0x3f800000) { /* y is +-1 */
if (hy < 0) if (hy < 0)
return one/x; return 1.0f/x;
return x; return x;
} }
if (hy == 0x40000000) /* y is 2 */ if (hy == 0x40000000) /* y is 2 */
@ -115,7 +112,7 @@ float powf(float x, float y)
if (ix == 0x7f800000 || ix == 0 || ix == 0x3f800000) { /* x is +-0,+-inf,+-1 */ if (ix == 0x7f800000 || ix == 0 || ix == 0x3f800000) { /* x is +-0,+-inf,+-1 */
z = ax; z = ax;
if (hy < 0) /* z = (1/|x|) */ if (hy < 0) /* z = (1/|x|) */
z = one/z; z = 1.0f/z;
if (hx < 0) { if (hx < 0) {
if (((ix-0x3f800000)|yisint) == 0) { if (((ix-0x3f800000)|yisint) == 0) {
z = (z-z)/(z-z); /* (-1)**non-int is NaN */ z = (z-z)/(z-z); /* (-1)**non-int is NaN */
@ -131,9 +128,9 @@ float powf(float x, float y)
if ((n|yisint) == 0) if ((n|yisint) == 0)
return (x-x)/(x-x); return (x-x)/(x-x);
sn = one; /* s (sign of result -ve**odd) = -1 else = 1 */ sn = 1.0f; /* s (sign of result -ve**odd) = -1 else = 1 */
if ((n|(yisint-1)) == 0) /* (-ve)**(odd int) */ if ((n|(yisint-1)) == 0) /* (-ve)**(odd int) */
sn = -one; sn = -1.0f;
/* |y| is huge */ /* |y| is huge */
if (iy > 0x4d000000) { /* if |y| > 2**27 */ if (iy > 0x4d000000) { /* if |y| > 2**27 */
@ -178,7 +175,7 @@ float powf(float x, float y)
/* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */ /* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
u = ax - bp[k]; /* bp[0]=1.0, bp[1]=1.5 */ u = ax - bp[k]; /* bp[0]=1.0, bp[1]=1.5 */
v = one/(ax+bp[k]); v = 1.0f/(ax+bp[k]);
s = u*v; s = u*v;
s_h = s; s_h = s;
GET_FLOAT_WORD(is, s_h); GET_FLOAT_WORD(is, s_h);
@ -257,8 +254,8 @@ float powf(float x, float y)
w = v - (z - u); w = v - (z - u);
t = z*z; t = z*z;
t1 = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5)))); t1 = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
r = (z*t1)/(t1-two) - (w+z*w); r = (z*t1)/(t1-2.0f) - (w+z*w);
z = one - (r - z); z = 1.0f - (r - z);
GET_FLOAT_WORD(j, z); GET_FLOAT_WORD(j, z);
j += n<<23; j += n<<23;
if ((j>>23) <= 0) /* subnormal output */ if ((j>>23) <= 0) /* subnormal output */

View File

@ -48,7 +48,7 @@ typedef uint32_t manh_t;
#define MANL_SHIFT (LDBL_MANL_SIZE - 1) #define MANL_SHIFT (LDBL_MANL_SIZE - 1)
static const long double Zero[] = {0.0L, -0.0L}; static const long double Zero[] = {0.0, -0.0};
/* /*
* Return the IEEE remainder and set *quo to the last n bits of the * Return the IEEE remainder and set *quo to the last n bits of the

View File

@ -79,7 +79,7 @@ long double rintl(long double x)
* If the result is +-0, then it must have the same sign as x, but * If the result is +-0, then it must have the same sign as x, but
* the above calculation doesn't always give this. Fix up the sign. * the above calculation doesn't always give this. Fix up the sign.
*/ */
if (ex < BIAS && x == 0.0L) if (ex < BIAS && x == 0.0)
return zero[sign]; return zero[sign];
return x; return x;

View File

@ -29,7 +29,7 @@ long double scalbnl(long double x, int n)
return x * 0x1p-16382L; return x * 0x1p-16382L;
} }
} }
scale.e = 1.0L; scale.e = 1.0;
scale.bits.exp = 0x3fff + n; scale.bits.exp = 0x3fff + n;
return x * scale.e; return x * scale.e;
} }

View File

@ -29,7 +29,7 @@
#include "libm.h" #include "libm.h"
static const double one = 1.0, huge = 1.0e307; static const double huge = 1.0e307;
double sinh(double x) double sinh(double x)
{ {
@ -50,12 +50,12 @@ double sinh(double x)
if (ix < 0x40360000) { /* |x|<22 */ if (ix < 0x40360000) { /* |x|<22 */
if (ix < 0x3e300000) /* |x|<2**-28 */ if (ix < 0x3e300000) /* |x|<2**-28 */
/* raise inexact, return x */ /* raise inexact, return x */
if (huge+x > one) if (huge+x > 1.0)
return x; return x;
t = expm1(fabs(x)); t = expm1(fabs(x));
if (ix < 0x3ff00000) if (ix < 0x3ff00000)
return h*(2.0*t - t*t/(t+one)); return h*(2.0*t - t*t/(t+1.0));
return h*(t + t/(t+one)); return h*(t + t/(t+1.0));
} }
/* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */ /* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */

View File

@ -15,7 +15,7 @@
#include "libm.h" #include "libm.h"
static const float one = 1.0, huge = 1.0e37; static const float huge = 1.0e37;
float sinhf(float x) float sinhf(float x)
{ {
@ -36,12 +36,12 @@ float sinhf(float x)
if (ix < 0x41100000) { /* |x|<9 */ if (ix < 0x41100000) { /* |x|<9 */
if (ix < 0x39800000) /* |x|<2**-12 */ if (ix < 0x39800000) /* |x|<2**-12 */
/* raise inexact, return x */ /* raise inexact, return x */
if (huge+x > one) if (huge+x > 1.0f)
return x; return x;
t = expm1f(fabsf(x)); t = expm1f(fabsf(x));
if (ix < 0x3f800000) if (ix < 0x3f800000)
return h*(2.0f*t - t*t/(t+one)); return h*(2.0f*t - t*t/(t+1.0f));
return h*(t + t/(t+one)); return h*(t + t/(t+1.0f));
} }
/* |x| in [9, logf(maxfloat)] return 0.5*exp(|x|) */ /* |x| in [9, logf(maxfloat)] return 0.5*exp(|x|) */

View File

@ -35,7 +35,7 @@ long double sinhl(long double x)
return sinh(x); return sinh(x);
} }
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
static const long double one = 1.0, huge = 1.0e4931L; static const long double huge = 1.0e4931L;
long double sinhl(long double x) long double sinhl(long double x)
{ {
@ -55,12 +55,12 @@ long double sinhl(long double x)
/* |x| in [0,25], return sign(x)*0.5*(E+E/(E+1))) */ /* |x| in [0,25], return sign(x)*0.5*(E+E/(E+1))) */
if (ix < 0x4003 || (ix == 0x4003 && i0 <= 0xc8000000)) { /* |x| < 25 */ if (ix < 0x4003 || (ix == 0x4003 && i0 <= 0xc8000000)) { /* |x| < 25 */
if (ix < 0x3fdf) /* |x|<2**-32 */ if (ix < 0x3fdf) /* |x|<2**-32 */
if (huge + x > one) if (huge + x > 1.0)
return x;/* sinh(tiny) = tiny with inexact */ return x;/* sinh(tiny) = tiny with inexact */
t = expm1l(fabsl(x)); t = expm1l(fabsl(x));
if (ix < 0x3fff) if (ix < 0x3fff)
return h*(2.0*t - t*t/(t + one)); return h*(2.0*t - t*t/(t + 1.0));
return h*(t + t/(t + one)); return h*(t + t/(t + 1.0));
} }
/* |x| in [25, log(maxdouble)] return 0.5*exp(|x|) */ /* |x| in [25, log(maxdouble)] return 0.5*exp(|x|) */

View File

@ -78,7 +78,7 @@
#include "libm.h" #include "libm.h"
static const double one = 1.0, tiny = 1.0e-300; static const double tiny = 1.0e-300;
double sqrt(double x) double sqrt(double x)
{ {
@ -161,13 +161,13 @@ double sqrt(double x)
/* use floating add to find out rounding direction */ /* use floating add to find out rounding direction */
if ((ix0|ix1) != 0) { if ((ix0|ix1) != 0) {
z = one - tiny; /* raise inexact flag */ z = 1.0 - tiny; /* raise inexact flag */
if (z >= one) { if (z >= 1.0) {
z = one + tiny; z = 1.0 + tiny;
if (q1 == (uint32_t)0xffffffff) { if (q1 == (uint32_t)0xffffffff) {
q1 = 0; q1 = 0;
q++; q++;
} else if (z > one) { } else if (z > 1.0) {
if (q1 == (uint32_t)0xfffffffe) if (q1 == (uint32_t)0xfffffffe)
q++; q++;
q1 += 2; q1 += 2;

View File

@ -15,7 +15,7 @@
#include "libm.h" #include "libm.h"
static const float one = 1.0, tiny = 1.0e-30; static const float tiny = 1.0e-30;
float sqrtf(float x) float sqrtf(float x)
{ {
@ -68,10 +68,10 @@ float sqrtf(float x)
/* use floating add to find out rounding direction */ /* use floating add to find out rounding direction */
if (ix != 0) { if (ix != 0) {
z = one - tiny; /* raise inexact flag */ z = 1.0f - tiny; /* raise inexact flag */
if (z >= one) { if (z >= 1.0f) {
z = one + tiny; z = 1.0f + tiny;
if (z > one) if (z > 1.0f)
q += 2; q += 2;
else else
q += q & 1; q += q & 1;

View File

@ -35,7 +35,7 @@
#include "libm.h" #include "libm.h"
static const double one = 1.0, two = 2.0, tiny = 1.0e-300, huge = 1.0e300; static const double tiny = 1.0e-300, huge = 1.0e300;
double tanh(double x) double tanh(double x)
{ {
@ -48,26 +48,26 @@ double tanh(double x)
/* x is INF or NaN */ /* x is INF or NaN */
if (ix >= 0x7ff00000) { if (ix >= 0x7ff00000) {
if (jx >= 0) if (jx >= 0)
return one/x + one; /* tanh(+-inf)=+-1 */ return 1.0f/x + 1.0f; /* tanh(+-inf)=+-1 */
else else
return one/x - one; /* tanh(NaN) = NaN */ return 1.0f/x - 1.0f; /* tanh(NaN) = NaN */
} }
if (ix < 0x40360000) { /* |x| < 22 */ if (ix < 0x40360000) { /* |x| < 22 */
if (ix < 0x3e300000) { /* |x| < 2**-28 */ if (ix < 0x3e300000) { /* |x| < 2**-28 */
/* tanh(tiny) = tiny with inexact */ /* tanh(tiny) = tiny with inexact */
if (huge+x > one) if (huge+x > 1.0f)
return x; return x;
} }
if (ix >= 0x3ff00000) { /* |x| >= 1 */ if (ix >= 0x3ff00000) { /* |x| >= 1 */
t = expm1(two*fabs(x)); t = expm1(2.0f*fabs(x));
z = one - two/(t+two); z = 1.0f - 2.0f/(t+2.0f);
} else { } else {
t = expm1(-two*fabs(x)); t = expm1(-2.0f*fabs(x));
z= -t/(t+two); z= -t/(t+2.0f);
} }
} else { /* |x| >= 22, return +-1 */ } else { /* |x| >= 22, return +-1 */
z = one - tiny; /* raise inexact */ z = 1.0f - tiny; /* raise inexact */
} }
return jx >= 0 ? z : -z; return jx >= 0 ? z : -z;
} }

View File

@ -15,7 +15,9 @@
#include "libm.h" #include "libm.h"
static const float one = 1.0, two = 2.0, tiny = 1.0e-30, huge = 1.0e30; static const float
tiny = 1.0e-30,
huge = 1.0e30;
float tanhf(float x) float tanhf(float x)
{ {
@ -28,26 +30,26 @@ float tanhf(float x)
/* x is INF or NaN */ /* x is INF or NaN */
if(ix >= 0x7f800000) { if(ix >= 0x7f800000) {
if (jx >= 0) if (jx >= 0)
return one/x + one; /* tanh(+-inf)=+-1 */ return 1.0f/x + 1.0f; /* tanh(+-inf)=+-1 */
else else
return one/x - one; /* tanh(NaN) = NaN */ return 1.0f/x - 1.0f; /* tanh(NaN) = NaN */
} }
if (ix < 0x41100000) { /* |x| < 9 */ if (ix < 0x41100000) { /* |x| < 9 */
if (ix < 0x39800000) { /* |x| < 2**-12 */ if (ix < 0x39800000) { /* |x| < 2**-12 */
/* tanh(tiny) = tiny with inexact */ /* tanh(tiny) = tiny with inexact */
if (huge+x > one) if (huge+x > 1.0f)
return x; return x;
} }
if (ix >= 0x3f800000) { /* |x|>=1 */ if (ix >= 0x3f800000) { /* |x|>=1 */
t = expm1f(two*fabsf(x)); t = expm1f(2.0f*fabsf(x));
z = one - two/(t+two); z = 1.0f - 2.0f/(t+2.0f);
} else { } else {
t = expm1f(-two*fabsf(x)); t = expm1f(-2.0f*fabsf(x));
z = -t/(t+two); z = -t/(t+2.0f);
} }
} else { /* |x| >= 9, return +-1 */ } else { /* |x| >= 9, return +-1 */
z = one - tiny; /* raise inexact */ z = 1.0f - tiny; /* raise inexact */
} }
return jx >= 0 ? z : -z; return jx >= 0 ? z : -z;
} }

View File

@ -41,7 +41,7 @@ long double tanhl(long double x)
return tanh(x); return tanh(x);
} }
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
static const long double one=1.0, two=2.0, tiny = 1.0e-4900L; static const long double tiny = 1.0e-4900L;
long double tanhl(long double x) long double tanhl(long double x)
{ {
@ -57,8 +57,8 @@ long double tanhl(long double x)
if (ix == 0x7fff) { if (ix == 0x7fff) {
/* for NaN it's not important which branch: tanhl(NaN) = NaN */ /* for NaN it's not important which branch: tanhl(NaN) = NaN */
if (se & 0x8000) if (se & 0x8000)
return one/x-one; /* tanhl(-inf)= -1; */ return 1.0/x-1.0; /* tanhl(-inf)= -1; */
return one/x+one; /* tanhl(+inf)= +1 */ return 1.0/x+1.0; /* tanhl(+inf)= +1 */
} }
/* |x| < 23 */ /* |x| < 23 */
@ -66,17 +66,17 @@ long double tanhl(long double x)
if ((ix|jj0|jj1) == 0) /* x == +- 0 */ if ((ix|jj0|jj1) == 0) /* x == +- 0 */
return x; return x;
if (ix < 0x3fc8) /* |x| < 2**-55 */ if (ix < 0x3fc8) /* |x| < 2**-55 */
return x*(one+tiny); /* tanh(small) = small */ return x*(1.0+tiny); /* tanh(small) = small */
if (ix >= 0x3fff) { /* |x| >= 1 */ if (ix >= 0x3fff) { /* |x| >= 1 */
t = expm1l(two*fabsl(x)); t = expm1l(2.0*fabsl(x));
z = one - two/(t+two); z = 1.0 - 2.0/(t+2.0);
} else { } else {
t = expm1l(-two*fabsl(x)); t = expm1l(-2.0*fabsl(x));
z = -t/(t+two); z = -t/(t+2.0);
} }
/* |x| > 23, return +-1 */ /* |x| > 23, return +-1 */
} else { } else {
z = one - tiny; /* raise inexact flag */ z = 1.0 - tiny; /* raise inexact flag */
} }
return se & 0x8000 ? -z : z; return se & 0x8000 ? -z : z;
} }

View File

@ -183,18 +183,18 @@ static long double stirf(long double x)
{ {
long double y, w, v; long double y, w, v;
w = 1.0L/x; w = 1.0/x;
/* For large x, use rational coefficients from the analytical expansion. */ /* For large x, use rational coefficients from the analytical expansion. */
if (x > 1024.0L) if (x > 1024.0)
w = (((((6.97281375836585777429E-5L * w w = (((((6.97281375836585777429E-5L * w
+ 7.84039221720066627474E-4L) * w + 7.84039221720066627474E-4L) * w
- 2.29472093621399176955E-4L) * w - 2.29472093621399176955E-4L) * w
- 2.68132716049382716049E-3L) * w - 2.68132716049382716049E-3L) * w
+ 3.47222222222222222222E-3L) * w + 3.47222222222222222222E-3L) * w
+ 8.33333333333333333333E-2L) * w + 8.33333333333333333333E-2L) * w
+ 1.0L; + 1.0;
else else
w = 1.0L + w * __polevll(w, STIR, 8); w = 1.0 + w * __polevll(w, STIR, 8);
y = expl(x); y = expl(x);
if (x > MAXSTIR) { /* Avoid overflow in pow() */ if (x > MAXSTIR) { /* Avoid overflow in pow() */
v = powl(x, 0.5L * x - 0.25L); v = powl(x, 0.5L * x - 0.25L);
@ -219,10 +219,10 @@ long double tgammal(long double x)
if (x == -INFINITY) if (x == -INFINITY)
return x - x; return x - x;
q = fabsl(x); q = fabsl(x);
if (q > 13.0L) { if (q > 13.0) {
if (q > MAXGAML) if (q > MAXGAML)
goto goverf; goto goverf;
if (x < 0.0L) { if (x < 0.0) {
p = floorl(q); p = floorl(q);
if (p == q) if (p == q)
return (x - x) / (x - x); return (x - x) / (x - x);
@ -231,7 +231,7 @@ long double tgammal(long double x)
signgam = -1; signgam = -1;
z = q - p; z = q - p;
if (z > 0.5L) { if (z > 0.5L) {
p += 1.0L; p += 1.0;
z = q - p; z = q - p;
} }
z = q * sinl(PIL * z); z = q * sinl(PIL * z);
@ -247,25 +247,25 @@ goverf:
return signgam * z; return signgam * z;
} }
z = 1.0L; z = 1.0;
while (x >= 3.0L) { while (x >= 3.0) {
x -= 1.0L; x -= 1.0;
z *= x; z *= x;
} }
while (x < -0.03125L) { while (x < -0.03125L) {
z /= x; z /= x;
x += 1.0L; x += 1.0;
} }
if (x <= 0.03125L) if (x <= 0.03125L)
goto small; goto small;
while (x < 2.0L) { while (x < 2.0) {
z /= x; z /= x;
x += 1.0L; x += 1.0;
} }
if (x == 2.0L) if (x == 2.0)
return z; return z;
x -= 2.0L; x -= 2.0;
p = __polevll(x, P, 7); p = __polevll(x, P, 7);
q = __polevll(x, Q, 8); q = __polevll(x, Q, 8);
z = z * p / q; z = z * p / q;
@ -274,9 +274,9 @@ goverf:
return z; return z;
small: small:
if (x == 0.0L) if (x == 0.0)
return (x - x) / (x - x); return (x - x) / (x - x);
if (x < 0.0L) { if (x < 0.0) {
x = -x; x = -x;
q = z / (x * __polevll(x, SN, 8)); q = z / (x * __polevll(x, SN, 8));
signgam = -1; signgam = -1;