diff options
Diffstat (limited to 'src')
95 files changed, 655 insertions, 806 deletions
diff --git a/src/internal/libm.h b/src/internal/libm.h index 2092763e..67c42b98 100644 --- a/src/internal/libm.h +++ b/src/internal/libm.h @@ -135,6 +135,7 @@ float __tandf(double,int); float __expo2f(float); float complex __ldexp_cexpf(float complex,int); +int __rem_pio2l(long double, long double *); long double __sinl(long double, long double, int); long double __cosl(long double, long double); long double __tanl(long double, long double, int); diff --git a/src/math/__cos.c b/src/math/__cos.c index ba439857..8699c1d5 100644 --- a/src/math/__cos.c +++ b/src/math/__cos.c @@ -51,7 +51,6 @@ #include "libm.h" static const double -one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ C1 = 4.16666666666666019037e-02, /* 0x3FA55555, 0x5555554C */ C2 = -1.38888888888741095749e-03, /* 0xBF56C16C, 0x16C15177 */ C3 = 2.48015872894767294178e-05, /* 0x3EFA01A0, 0x19CB1590 */ @@ -67,6 +66,6 @@ double __cos(double x, double y) w = z*z; r = z*(C1+z*(C2+z*C3)) + w*w*(C4+z*(C5+z*C6)); hz = 0.5*z; - w = one-hz; - return w + (((one-w)-hz) + (z*r-x*y)); + w = 1.0-hz; + return w + (((1.0-w)-hz) + (z*r-x*y)); } diff --git a/src/math/__cosdf.c b/src/math/__cosdf.c index a3b399e6..a65f7f21 100644 --- a/src/math/__cosdf.c +++ b/src/math/__cosdf.c @@ -18,7 +18,6 @@ /* |cos(x) - c(x)| < 2**-34.1 (~[-5.37e-11, 5.295e-11]). */ static const double -one = 1.0, C0 = -0x1ffffffd0c5e81.0p-54, /* -0.499999997251031003120 */ C1 = 0x155553e1053a42.0p-57, /* 0.0416666233237390631894 */ C2 = -0x16c087e80f1e27.0p-62, /* -0.00138867637746099294692 */ @@ -32,5 +31,5 @@ float __cosdf(double x) z = x*x; w = z*z; r = C2+z*C3; - return ((one+z*C0) + w*C1) + (w*z)*r; + return ((1.0+z*C0) + w*C1) + (w*z)*r; } diff --git a/src/math/__cosl.c b/src/math/__cosl.c index 80036ddb..9d325768 100644 --- a/src/math/__cosl.c +++ b/src/math/__cosl.c @@ -41,8 +41,6 @@ * almost for free from the complications needed to search for the best * higher coefficients. */ -static const double one = 1.0; - static const long double C1 = 0.0416666666666666666136L; /* 0xaaaaaaaaaaaaaa9b.0p-68 */ @@ -61,7 +59,7 @@ long double __cosl(long double x, long double y) z = x*x; r = z*(C1+z*(C2+z*(C3+z*(C4+z*(C5+z*(C6+z*C7)))))); hz = 0.5*z; - w = one-hz; - return w + (((one-w)-hz) + (z*r-x*y)); + w = 1.0-hz; + return w + (((1.0-w)-hz) + (z*r-x*y)); } #endif diff --git a/src/math/__rem_pio2.c b/src/math/__rem_pio2.c index a7d779e0..0ef57fb5 100644 --- a/src/math/__rem_pio2.c +++ b/src/math/__rem_pio2.c @@ -29,7 +29,6 @@ * pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3) */ static const double -zero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */ two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */ invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */ pio2_1 = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */ @@ -163,7 +162,7 @@ medium: } tx[2] = z; 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); if (hx < 0) { y[0] = -ty[0]; diff --git a/src/math/__rem_pio2_large.c b/src/math/__rem_pio2_large.c index 35835f83..27b619cc 100644 --- a/src/math/__rem_pio2_large.c +++ b/src/math/__rem_pio2_large.c @@ -271,8 +271,6 @@ static const double PIo2[] = { }; static const double -zero = 0.0, -one = 1.0, two24 = 1.67772160000000000000e+07, /* 0x41700000, 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] */ j = jv-jx; m = jx+jk; 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] */ for (i=0; i<=jk; i++) { @@ -346,14 +344,14 @@ recompute: } } if (ih == 2) { - z = one - z; + z = 1.0 - z; if (carry != 0) - z -= scalbn(one,q0); + z -= scalbn(1.0,q0); } } /* check if recomputation is needed */ - if (z == zero) { + if (z == 0.0) { j = 0; for (i=jz-1; i>=jk; i--) j |= iq[i]; if (j == 0) { /* need recomputation */ @@ -391,7 +389,7 @@ recompute: } /* convert integer "bit" chunk to floating-point value */ - fw = scalbn(one,q0); + fw = scalbn(1.0,q0); for (i=jz; i>=0; i--) { q[i] = fw*(double)iq[i]; fw *= twon24; diff --git a/src/math/__rem_pio2l.h b/src/math/__rem_pio2l.c index 11123c3f..f0cb99ac 100644 --- a/src/math/__rem_pio2l.h +++ b/src/math/__rem_pio2l.c @@ -32,7 +32,6 @@ * pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3) */ static const double -zero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */ two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */ pio2_1 = 1.57079632679597125389e+00, /* 0x3FF921FB, 0x54444000 */ pio2_2 = -1.07463465549783099519e-12, /* -0x12e7b967674000.0p-92 */ @@ -44,7 +43,7 @@ pio2_1t = -1.07463465549719416346e-12L, /* -0x973dcb3b399d747f.0p-103 */ pio2_2t = 6.36831716351095013979e-25L, /* 0xc51701b839a25205.0p-144 */ pio2_3t = -2.75299651904407171810e-37L; /* -0xbb5bf6c7ddd660ce.0p-185 */ -static inline int __rem_pio2l(long double x, long double *y) +int __rem_pio2l(long double x, long double *y) { union IEEEl2bits u,u1; long double z,w,t,r,fn; @@ -114,7 +113,7 @@ static inline int __rem_pio2l(long double x, long double *y) } tx[2] = z; nx = 3; - while (tx[nx-1] == zero) + while (tx[nx-1] == 0.0) nx--; /* skip zero term */ n = __rem_pio2_large(tx,ty,e0,nx,2); r = (long double)ty[0] + ty[1]; diff --git a/src/math/__sin.c b/src/math/__sin.c index 80f3273c..9aead04b 100644 --- a/src/math/__sin.c +++ b/src/math/__sin.c @@ -42,7 +42,6 @@ #include "libm.h" static const double -half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */ S1 = -1.66666666666666324348e-01, /* 0xBFC55555, 0x55555549 */ S2 = 8.33333333332248946124e-03, /* 0x3F811111, 0x1110F8A6 */ S3 = -1.98412698298579493134e-04, /* 0xBF2A01A0, 0x19C161D5 */ @@ -61,5 +60,5 @@ double __sin(double x, double y, int iy) if (iy == 0) return x + v*(S1 + z*r); else - return x - ((z*(half*y - v*r) - y) - v*S1); + return x - ((z*(0.5*y - v*r) - y) - v*S1); } diff --git a/src/math/__sinl.c b/src/math/__sinl.c index 67c4bdc5..068adffb 100644 --- a/src/math/__sinl.c +++ b/src/math/__sinl.c @@ -24,8 +24,6 @@ * See __cosl.c for more details about the polynomial. */ -static const double half = 0.5; - static const long double 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))))); if (iy == 0) 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 diff --git a/src/math/__tan.c b/src/math/__tan.c index f1be2ec8..01e3fe48 100644 --- a/src/math/__tan.c +++ b/src/math/__tan.c @@ -59,13 +59,9 @@ static const double T[] = { 7.14072491382608190305e-05, /* 3F12B80F, 32F0A7E9 */ -1.85586374855275456654e-05, /* BEF375CB, DB605373 */ 2.59073051863633712884e-05, /* 3EFB2A70, 74BF7AD4 */ -/* one */ 1.00000000000000000000e+00, /* 3FF00000, 00000000 */ -/* pio4 */ 7.85398163397448278999e-01, /* 3FE921FB, 54442D18 */ -/* pio4lo */ 3.06161699786838301793e-17 /* 3C81A626, 33145C07 */ -}; -#define one T[13] -#define pio4 T[14] -#define pio4lo T[15] +}, +pio4 = 7.85398163397448278999e-01, /* 3FE921FB, 54442D18 */ +pio4lo = 3.06161699786838301793e-17; /* 3C81A626, 33145C07 */ double __tan(double x, double y, int iy) { diff --git a/src/math/__tanl.c b/src/math/__tanl.c index e39e9df4..50ba2140 100644 --- a/src/math/__tanl.c +++ b/src/math/__tanl.c @@ -51,8 +51,7 @@ long double __tanl(long double x, long double y, int iy) { int i; iy = iy == 1 ? -1 : 1; /* XXX recover original interface */ - // FIXME: this is wrong, use copysign, signbit or union bithack - osign = x >= 0 ? 1.0 : -1.0; /* XXX slow, probably wrong for -0 */ + osign = copysignl(1.0, x); if (fabsl(x) >= 0.67434) { if (x < 0) { x = -x; diff --git a/src/math/acos.c b/src/math/acos.c index 456a2219..54d266ee 100644 --- a/src/math/acos.c +++ b/src/math/acos.c @@ -36,8 +36,7 @@ #include "libm.h" 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 */ static const volatile double pio2_lo = 6.12323399573676603587e-17; /* 0x3C91A626, 0x33145C07 */ @@ -75,25 +74,25 @@ double acos(double x) return pio2_hi + pio2_lo; z = x*x; 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; return pio2_hi - (x - (pio2_lo-x*r)); } 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))))); - q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4))); + q = 1.0+z*(qS1+z*(qS2+z*(qS3+z*qS4))); s = sqrt(z); r = p/q; w = r*s-pio2_lo; return pi - 2.0*(s+w); } else { /* x > 0.5 */ - z = (one-x)*0.5; + z = (1.0-x)*0.5; s = sqrt(z); df = s; SET_LOW_WORD(df,0); c = (z-df*df)/(s+df); 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; w = r*s+c; return 2.0*(df+w); diff --git a/src/math/acosf.c b/src/math/acosf.c index 945347cb..d875f3ef 100644 --- a/src/math/acosf.c +++ b/src/math/acosf.c @@ -16,8 +16,7 @@ #include "libm.h" static const float -one = 1.0000000000e+00, /* 0x3F800000 */ -pi = 3.1415925026e+00, /* 0x40490fda */ +pi = 3.1415925026e+00, /* 0x40490fda */ pio2_hi = 1.5707962513e+00; /* 0x3fc90fda */ static const volatile float pio2_lo = 7.5497894159e-08; /* 0x33a22168 */ @@ -46,13 +45,13 @@ float acosf(float x) return pio2_hi + pio2_lo; z = x*x; p = z*(pS0+z*(pS1+z*pS2)); - q = one+z*qS1; + q = 1.0f+z*qS1; r = p/q; return pio2_hi - (x - (pio2_lo-x*r)); } 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)); - q = one+z*qS1; + q = 1.0f+z*qS1; s = sqrtf(z); r = p/q; w = r*s-pio2_lo; @@ -60,14 +59,14 @@ float acosf(float x) } else { /* x > 0.5 */ int32_t idf; - z = (one-x)*0.5f; + z = (1.0f-x)*0.5f; s = sqrtf(z); df = s; GET_FLOAT_WORD(idf,df); SET_FLOAT_WORD(df,idf&0xfffff000); c = (z-df*df)/(s+df); p = z*(pS0+z*(pS1+z*pS2)); - q = one+z*qS1; + q = 1.0f+z*qS1; r = p/q; w = r*s+c; return 2.0f*(df+w); diff --git a/src/math/acosh.c b/src/math/acosh.c index a7c87e3c..15f51c6e 100644 --- a/src/math/acosh.c +++ b/src/math/acosh.c @@ -27,7 +27,6 @@ #include "libm.h" static const double -one = 1.0, ln2 = 6.93147180559945286227e-01; /* 0x3FE62E42, 0xFEFA39EF */ double acosh(double x) @@ -47,9 +46,9 @@ double acosh(double x) return 0.0; /* acosh(1) = 0 */ } else if (hx > 0x40000000) { /* 2**28 > x > 2 */ 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 */ - t = x-one; + t = x-1.0; return log1p(t + sqrt(2.0*t+t*t)); } } diff --git a/src/math/acoshf.c b/src/math/acoshf.c index 57ce5cb8..0f7aae2a 100644 --- a/src/math/acoshf.c +++ b/src/math/acoshf.c @@ -16,7 +16,6 @@ #include "libm.h" static const float -one = 1.0, ln2 = 6.9314718246e-01; /* 0x3f317218 */ float acoshf(float x) @@ -32,12 +31,12 @@ float acoshf(float x) return x + x; return logf(x) + ln2; /* acosh(huge)=log(2x) */ } else if (hx == 0x3f800000) { - return 0.0; /* acosh(1) = 0 */ + return 0.0f; /* acosh(1) = 0 */ } else if (hx > 0x40000000) { /* 2**28 > x > 2 */ 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 */ - t = x-one; + t = x-1.0f; return log1pf(t + sqrtf(2.0f*t+t*t)); } } diff --git a/src/math/acoshl.c b/src/math/acoshl.c index d8310a73..a4024516 100644 --- a/src/math/acoshl.c +++ b/src/math/acoshl.c @@ -32,7 +32,6 @@ long double acoshl(long double x) } #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 static const long double -one = 1.0, ln2 = 6.931471805599453094287e-01L; /* 0x3FFE, 0xB17217F7, 0xD1CF79AC */ long double acoshl(long double x) @@ -51,10 +50,10 @@ long double acoshl(long double x) return 0.0; /* acosh(1) = 0 */ } else if (se > 0x4000) { /* x > 2 */ 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 */ - t = x - one; + t = x - 1.0; return log1pl(t + sqrtl(2.0*t + t*t)); } #endif diff --git a/src/math/acosl.c b/src/math/acosl.c index 170520fe..cc565336 100644 --- a/src/math/acosl.c +++ b/src/math/acosl.c @@ -25,7 +25,6 @@ long double acosl(long double x) #include "__invtrigl.h" static const long double -one = 1.00000000000000000000e+00, pi = 3.14159265358979323846264338327950280e+00L; long double acosl(long double x) @@ -55,7 +54,7 @@ long double acosl(long double x) r = p / q; return pio2_hi - (x - (pio2_lo - x * r)); } else if (expsign < 0) { /* x < -0.5 */ - z = (one + x) * 0.5; + z = (1.0 + x) * 0.5; p = P(z); q = Q(z); s = sqrtl(z); @@ -63,7 +62,7 @@ long double acosl(long double x) w = r * s - pio2_lo; return pi - 2.0 * (s + w); } else { /* x > 0.5 */ - z = (one - x) * 0.5; + z = (1.0 - x) * 0.5; s = sqrtl(z); u.e = s; u.bits.manl = 0; diff --git a/src/math/asin.c b/src/math/asin.c index 04bd0c14..e3d8b250 100644 --- a/src/math/asin.c +++ b/src/math/asin.c @@ -42,7 +42,6 @@ #include "libm.h" static const double -one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ huge = 1.000e+300, pio2_hi = 1.57079632679489655800e+00, /* 0x3FF921FB, 0x54442D18 */ 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 */ } else if (ix < 0x3fe00000) { /* |x|<0.5 */ 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*/ } t = x*x; 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; return x + x*w; } /* 1 > |x| >= 0.5 */ - w = one - fabs(x); + w = 1.0 - fabs(x); t = w*0.5; 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); if (ix >= 0x3FEF3333) { /* if |x| > 0.975 */ w = p/q; diff --git a/src/math/asinf.c b/src/math/asinf.c index 7865a45b..c1c42c7b 100644 --- a/src/math/asinf.c +++ b/src/math/asinf.c @@ -16,7 +16,6 @@ #include "libm.h" static const float -one = 1.0000000000e+00, /* 0x3F800000 */ huge = 1.000e+30, /* coefficients for R(x^2) */ pS0 = 1.6666586697e-01, @@ -41,20 +40,20 @@ float asinf(float x) return (x-x)/(x-x); /* asin(|x|>1) is NaN */ } else if (ix < 0x3f000000) { /* |x|<0.5 */ if (ix < 0x39800000) { /* |x| < 2**-12 */ - if (huge+x > one) + if (huge+x > 1.0f) return x; /* return x with inexact if x!=0 */ } t = x*x; p = t*(pS0+t*(pS1+t*pS2)); - q = one+t*qS1; + q = 1.0f+t*qS1; w = p/q; return x + x*w; } /* 1 > |x| >= 0.5 */ - w = one - fabsf(x); + w = 1.0f - fabsf(x); t = w*0.5f; p = t*(pS0+t*(pS1+t*pS2)); - q = one+t*qS1; + q = 1.0f+t*qS1; s = sqrt(t); w = p/q; t = pio2-2.0*(s+s*w); diff --git a/src/math/asinh.c b/src/math/asinh.c index 92aa9446..11bbd71a 100644 --- a/src/math/asinh.c +++ b/src/math/asinh.c @@ -23,7 +23,6 @@ #include "libm.h" static const double -one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ ln2 = 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */ huge= 1.00000000000000000000e+300; @@ -38,17 +37,17 @@ double asinh(double x) return x+x; if (ix < 0x3e300000) { /* |x| < 2**-28 */ /* return x inexact except 0 */ - if (huge+x > one) + if (huge+x > 1.0) return x; } if (ix > 0x41b00000) { /* |x| > 2**28 */ w = log(fabs(x)) + ln2; } else if (ix > 0x40000000) { /* 2**28 > |x| > 2.0 */ 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 */ 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) return w; diff --git a/src/math/asinhf.c b/src/math/asinhf.c index cb6e5b89..efe3af94 100644 --- a/src/math/asinhf.c +++ b/src/math/asinhf.c @@ -16,7 +16,6 @@ #include "libm.h" static const float -one = 1.0000000000e+00, /* 0x3F800000 */ ln2 = 6.9314718246e-01, /* 0x3f317218 */ huge= 1.0000000000e+30; @@ -31,17 +30,17 @@ float asinhf(float x) return x+x; if (ix < 0x31800000) { /* |x| < 2**-28 */ /* return x inexact except 0 */ - if (huge+x > one) + if (huge+x > 1.0f) return x; } if (ix > 0x4d800000) { /* |x| > 2**28 */ w = logf(fabsf(x)) + ln2; } else if (ix > 0x40000000) { /* 2**28 > |x| > 2.0 */ 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 */ 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) return w; diff --git a/src/math/asinhl.c b/src/math/asinhl.c index b2edf904..dc5dd71f 100644 --- a/src/math/asinhl.c +++ b/src/math/asinhl.c @@ -29,7 +29,6 @@ long double asinhl(long double x) } #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 static const long double -one = 1.000000000000000000000e+00L, /* 0x3FFF, 0x00000000, 0x00000000 */ ln2 = 6.931471805599453094287e-01L, /* 0x3FFE, 0xB17217F7, 0xD1CF79AC */ huge = 1.000000000000000000e+4900L; @@ -44,17 +43,17 @@ long double asinhl(long double x) return x + x; /* x is inf or NaN */ if (ix < 0x3fde) { /* |x| < 2**-34 */ /* return x, raise inexact if x != 0 */ - if (huge+x > one) + if (huge+x > 1.0) return x; } if (ix > 0x4020) { /* |x| > 2**34 */ w = logl(fabsl(x)) + ln2; } else if (ix > 0x4000) { /* 2**34 > |x| > 2.0 */ 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 */ 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) return -w; diff --git a/src/math/asinl.c b/src/math/asinl.c index 370997bc..ddd807e2 100644 --- a/src/math/asinl.c +++ b/src/math/asinl.c @@ -23,9 +23,7 @@ long double asinl(long double x) } #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 #include "__invtrigl.h" -static const long double -one = 1.00000000000000000000e+00, -huge = 1.000e+300; +static const long double huge = 1.000e+300; long double asinl(long double x) { @@ -45,7 +43,7 @@ long double asinl(long double x) } else if (expt < BIAS-1) { /* |x|<0.5 */ if (expt < ASIN_LINEAR) { /* if |x| is small, asinl(x)=x */ /* return x with inexact if x!=0 */ - if (huge+x > one) + if (huge+x > 1.0) return x; } t = x*x; @@ -55,7 +53,7 @@ long double asinl(long double x) return x + x*w; } /* 1 > |x| >= 0.5 */ - w = one - fabsl(x); + w = 1.0 - fabsl(x); t = w*0.5; p = P(t); q = Q(t); diff --git a/src/math/atan.c b/src/math/atan.c index d31782c2..f34551c7 100644 --- a/src/math/atan.c +++ b/src/math/atan.c @@ -60,9 +60,7 @@ static const double aT[] = { 1.62858201153657823623e-02, /* 0x3F90AD3A, 0xE322DA11 */ }; -static const double -one = 1.0, -huge = 1.0e300; +static const double huge = 1.0e300; double atan(double x) { @@ -86,7 +84,7 @@ double atan(double x) if (ix < 0x3fdc0000) { /* |x| < 0.4375 */ if (ix < 0x3e400000) { /* |x| < 2^-27 */ /* raise inexact */ - if (huge+x > one) + if (huge+x > 1.0) return x; } id = -1; @@ -95,15 +93,15 @@ double atan(double x) if (ix < 0x3ff30000) { /* |x| < 1.1875 */ if (ix < 0x3fe60000) { /* 7/16 <= |x| < 11/16 */ 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 */ id = 1; - x = (x-one)/(x+one); + x = (x-1.0)/(x+1.0); } } else { if (ix < 0x40038000) { /* |x| < 2.4375 */ 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 */ id = 3; x = -1.0/x; diff --git a/src/math/atan2.c b/src/math/atan2.c index d928f0ed..143c3834 100644 --- a/src/math/atan2.c +++ b/src/math/atan2.c @@ -42,7 +42,6 @@ static const volatile double tiny = 1.0e-300; static const double -zero = 0.0, pi_o_4 = 7.8539816339744827900E-01, /* 0x3FE921FB, 0x54442D18 */ pi_o_2 = 1.5707963267948965580E+00, /* 0x3FF921FB, 0x54442D18 */ pi = 3.1415926535897931160E+00; /* 0x400921FB, 0x54442D18 */ @@ -89,8 +88,8 @@ double atan2(double y, double x) } } else { switch(m) { - case 0: return zero; /* atan(+...,+INF) */ - case 1: return -zero; /* atan(-...,+INF) */ + case 0: return 0.0; /* atan(+...,+INF) */ + case 1: return -0.0; /* atan(-...,+INF) */ case 2: return pi+tiny; /* atan(+...,-INF) */ case 3: return -pi-tiny; /* atan(-...,-INF) */ } diff --git a/src/math/atan2f.c b/src/math/atan2f.c index 19b134dc..96839cba 100644 --- a/src/math/atan2f.c +++ b/src/math/atan2f.c @@ -18,7 +18,6 @@ static const volatile float tiny = 1.0e-30; static const float -zero = 0.0, pi_o_4 = 7.8539818525e-01, /* 0x3f490fdb */ pi_o_2 = 1.5707963705e+00, /* 0x3fc90fdb */ pi = 3.1415927410e+00; /* 0x40490fdb */ @@ -63,8 +62,8 @@ float atan2f(float y, float x) } } else { switch (m) { - case 0: return zero; /* atan(+...,+INF) */ - case 1: return -zero; /* atan(-...,+INF) */ + case 0: return 0.0f; /* atan(+...,+INF) */ + case 1: return -0.0f; /* atan(-...,+INF) */ case 2: return pi+tiny; /* atan(+...,-INF) */ case 3: return -pi-tiny; /* atan(-...,-INF) */ } diff --git a/src/math/atan2l.c b/src/math/atan2l.c index 48abc058..0fc901c8 100644 --- a/src/math/atan2l.c +++ b/src/math/atan2l.c @@ -27,7 +27,6 @@ long double atan2l(long double y, long double x) static const volatile long double tiny = 1.0e-300; static const long double -zero = 0.0, pi = 3.14159265358979323846264338327950280e+00L; long double atan2l(long double y, long double x) @@ -75,8 +74,8 @@ long double atan2l(long double y, long double x) } } else { switch(m) { - case 0: return zero; /* atan(+...,+INF) */ - case 1: return -zero; /* atan(-...,+INF) */ + case 0: return 0.0; /* atan(+...,+INF) */ + case 1: return -0.0; /* atan(-...,+INF) */ case 2: return pi+tiny; /* atan(+...,-INF) */ case 3: return -pi-tiny; /* atan(-...,-INF) */ } diff --git a/src/math/atanf.c b/src/math/atanf.c index 73cebb5e..4e31b902 100644 --- a/src/math/atanf.c +++ b/src/math/atanf.c @@ -38,9 +38,7 @@ static const float aT[] = { 6.1687607318e-02, }; -static const float -one = 1.0, -huge = 1.0e30; +static const float huge = 1.0e30; float atanf(float x) { @@ -60,7 +58,7 @@ float atanf(float x) if (ix < 0x3ee00000) { /* |x| < 0.4375 */ if (ix < 0x39800000) { /* |x| < 2**-12 */ /* raise inexact */ - if(huge+x>one) + if(huge+x>1.0f) return x; } id = -1; @@ -69,15 +67,15 @@ float atanf(float x) if (ix < 0x3f980000) { /* |x| < 1.1875 */ if (ix < 0x3f300000) { /* 7/16 <= |x| < 11/16 */ 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 */ id = 1; - x = (x - one)/(x + one); + x = (x - 1.0f)/(x + 1.0f); } } else { if (ix < 0x401c0000) { /* |x| < 2.4375 */ 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 */ id = 3; x = -1.0f/x; diff --git a/src/math/atanh.c b/src/math/atanh.c index 29290463..dbe241d1 100644 --- a/src/math/atanh.c +++ b/src/math/atanh.c @@ -30,8 +30,7 @@ #include "libm.h" -static const double one = 1.0, huge = 1e300; -static const double zero = 0.0; +static const double huge = 1e300; double atanh(double x) { @@ -44,15 +43,15 @@ double atanh(double x) if ((ix | ((lx|-lx)>>31)) > 0x3ff00000) /* |x| > 1 */ return (x-x)/(x-x); if (ix == 0x3ff00000) - return x/zero; - if (ix < 0x3e300000 && (huge+x) > zero) /* x < 2**-28 */ + return x/0.0; + if (ix < 0x3e300000 && (huge+x) > 0.0) /* x < 2**-28 */ return x; SET_HIGH_WORD(x, ix); if (ix < 0x3fe00000) { /* x < 0.5 */ t = x+x; - t = 0.5*log1p(t + t*x/(one-x)); + t = 0.5*log1p(t + t*x/(1.0-x)); } else - t = 0.5*log1p((x+x)/(one-x)); + t = 0.5*log1p((x+x)/(1.0-x)); if (hx >= 0) return t; return -t; diff --git a/src/math/atanhf.c b/src/math/atanhf.c index 9c65dc7f..2be780bb 100644 --- a/src/math/atanhf.c +++ b/src/math/atanhf.c @@ -15,8 +15,7 @@ #include "libm.h" -static const float one = 1.0, huge = 1e30; -static const float zero = 0.0; +static const float huge = 1e30; float atanhf(float x) { @@ -28,15 +27,15 @@ float atanhf(float x) if (ix > 0x3f800000) /* |x| > 1 */ return (x-x)/(x-x); if (ix == 0x3f800000) - return x/zero; - if (ix < 0x31800000 && huge+x > zero) /* x < 2**-28 */ + return x/0.0f; + if (ix < 0x31800000 && huge+x > 0.0f) /* x < 2**-28 */ return x; SET_FLOAT_WORD(x, ix); if (ix < 0x3f000000) { /* x < 0.5 */ t = x+x; - t = 0.5f*log1pf(t + t*x/(one-x)); + t = 0.5f*log1pf(t + t*x/(1.0f-x)); } else - t = 0.5f*log1pf((x+x)/(one-x)); + t = 0.5f*log1pf((x+x)/(1.0f-x)); if (hx >= 0) return t; return -t; diff --git a/src/math/atanhl.c b/src/math/atanhl.c index af0f856d..931bae32 100644 --- a/src/math/atanhl.c +++ b/src/math/atanhl.c @@ -34,7 +34,7 @@ long double atanhl(long double x) return atanh(x); } #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) { @@ -48,15 +48,15 @@ long double atanhl(long double x) /* |x| > 1 */ return (x-x)/(x-x); if (ix == 0x3fff) - return x/zero; - if (ix < 0x3fe3 && huge+x > zero) /* x < 2**-28 */ + return x/0.0; + if (ix < 0x3fe3 && huge+x > 0.0) /* x < 2**-28 */ return x; SET_LDOUBLE_EXP(x, ix); if (ix < 0x3ffe) { /* x < 0.5 */ t = x + x; - t = 0.5*log1pl(t + t*x/(one-x)); + t = 0.5*log1pl(t + t*x/(1.0 - x)); } else - t = 0.5*log1pl((x + x)/(one - x)); + t = 0.5*log1pl((x + x)/(1.0 - x)); if (se <= 0x7fff) return t; return -t; diff --git a/src/math/atanl.c b/src/math/atanl.c index 4e99955e..36072c17 100644 --- a/src/math/atanl.c +++ b/src/math/atanl.c @@ -23,9 +23,7 @@ long double atanl(long double x) } #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 #include "__invtrigl.h" -static const long double -one = 1.0, -huge = 1.0e300; +static const long double huge = 1.0e300; 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 (expt < ATAN_LINEAR) { /* if |x| is small, atanl(x)~=x */ /* raise inexact */ - if (huge+x > one) + if (huge+x > 1.0) return x; } id = -1; @@ -62,15 +60,15 @@ long double atanl(long double x) if (expman < (BIAS << 8) + 0x30) { /* |x| < 1.1875 */ if (expman < ((BIAS - 1) << 8) + 0x60) { /* 7/16 <= |x| < 11/16 */ 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 */ id = 1; - x = (x-one)/(x+one); + x = (x-1.0)/(x+1.0); } } else { if (expman < ((BIAS + 1) << 8) + 0x38) { /* |x| < 2.4375 */ 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 */ id = 3; x = -1.0/x; diff --git a/src/math/cbrtl.c b/src/math/cbrtl.c index d138b9f2..5297d68f 100644 --- a/src/math/cbrtl.c +++ b/src/math/cbrtl.c @@ -118,11 +118,7 @@ long double cbrtl(long double x) * Round it away from zero to 32 bits (32 so that t*t is exact, and * away from zero for technical reasons). */ - volatile double vd2 = 0x1.0p32; - volatile double vd1 = 0x1.0p-31; - #define vd ((long double)vd2 + vd1) - - t = dt + vd - 0x1.0p32; + t = dt + (0x1.0p32L + 0x1.0p-32L) - 0x1.0p32; #elif LDBL_MANT_DIG == 113 /* * Round dt away from zero to 47 bits. Since we don't trust the 47, diff --git a/src/math/cosh.c b/src/math/cosh.c index 5f38b276..a1f7dbc7 100644 --- a/src/math/cosh.c +++ b/src/math/cosh.c @@ -32,7 +32,7 @@ #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) { @@ -49,21 +49,21 @@ double cosh(double x) /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */ if (ix < 0x3fd62e43) { t = expm1(fabs(x)); - w = one+t; + w = 1.0+t; if (ix < 0x3c800000) 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; */ if (ix < 0x40360000) { 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) - return half*exp(fabs(x)); + return 0.5*exp(fabs(x)); /* |x| in [log(maxdouble), overflowthresold] */ if (ix <= 0x408633CE) diff --git a/src/math/coshf.c b/src/math/coshf.c index 9e87afcd..97318f10 100644 --- a/src/math/coshf.c +++ b/src/math/coshf.c @@ -15,7 +15,7 @@ #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) { @@ -32,21 +32,21 @@ float coshf(float x) /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */ if (ix < 0x3eb17218) { t = expm1f(fabsf(x)); - w = one+t; + w = 1.0f+t; if (ix<0x39800000) - return one; /* cosh(tiny) = 1 */ - return one + (t*t)/(w+w); + return 1.0f; /* cosh(tiny) = 1 */ + return 1.0f + (t*t)/(w+w); } /* |x| in [0.5*ln2,9], return (exp(|x|)+1/exp(|x|))/2; */ if (ix < 0x41100000) { 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) - return half*expf(fabsf(x)); + return 0.5f*expf(fabsf(x)); /* |x| in [log(maxfloat), overflowthresold] */ if (ix <= 0x42b2d4fc) diff --git a/src/math/coshl.c b/src/math/coshl.c index bcc9128a..36b52c02 100644 --- a/src/math/coshl.c +++ b/src/math/coshl.c @@ -38,7 +38,7 @@ long double coshl(long double x) return cosh(x); } #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) { @@ -56,27 +56,27 @@ long double coshl(long double x) /* |x| in [0,0.5*ln2], return 1+expm1l(|x|)^2/(2*expl(|x|)) */ if (ex < 0x3ffd || (ex == 0x3ffd && mx < 0xb17217f7u)) { t = expm1l(fabsl(x)); - w = one + t; + w = 1.0 + t; 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; */ if (ex < 0x4003 || (ex == 0x4003 && mx < 0xb0000000u)) { 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)) - return half*expl(fabsl(x)); + return 0.5*expl(fabsl(x)); /* |x| in [log(maxdouble), log(2*maxdouble)) */ if (ex == 0x400c && (mx < 0xb174ddc0u || (mx == 0xb174ddc0u && lx < 0x31aec0ebu))) { - w = expl(half*fabsl(x)); - t = half*w; + w = expl(0.5*fabsl(x)); + t = 0.5*w; return t*w; } diff --git a/src/math/cosl.c b/src/math/cosl.c index 2c650cdc..25d9005a 100644 --- a/src/math/cosl.c +++ b/src/math/cosl.c @@ -36,8 +36,6 @@ long double cosl(long double x) { return cos(x); } #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 -#include "__rem_pio2l.h" - long double cosl(long double x) { union IEEEl2bits z; diff --git a/src/math/erf.c b/src/math/erf.c index 18ee01cf..dd8c3a77 100644 --- a/src/math/erf.c +++ b/src/math/erf.c @@ -107,9 +107,6 @@ static const double 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 */ erx = 8.45062911510467529297e-01, /* 0x3FEB0AC1, 0x60000000 */ /* @@ -190,7 +187,7 @@ double erf(double x) if (ix >= 0x7ff00000) { /* erf(nan)=nan, erf(+-inf)=+-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 < 0x3e300000) { /* |x|<2**-28 */ @@ -201,42 +198,42 @@ double erf(double x) } z = x*x; r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4))); - s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5)))); + s = 1.0+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5)))); y = r/s; return x + x*y; } if (ix < 0x3ff40000) { /* 0.84375 <= |x| < 1.25 */ - s = fabs(x)-one; + s = fabs(x)-1.0; 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) return erx + P/Q; return -erx - P/Q; } if (ix >= 0x40180000) { /* inf > |x| >= 6 */ if (hx >= 0) - return one-tiny; - return tiny-one; + return 1.0 - tiny; + return tiny - 1.0; } x = fabs(x); - s = one/(x*x); + s = 1.0/(x*x); if (ix < 0x4006DB6E) { /* |x| < 1/0.35 */ R = ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*( ra5+s*(ra6+s*ra7)))))); - S = one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*( + S = 1.0+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*( sa5+s*(sa6+s*(sa7+s*sa8))))))); } else { /* |x| >= 1/0.35 */ R = rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*( rb5+s*rb6))))); - S = one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*( + S = 1.0+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*( sb5+s*(sb6+s*sb7)))))); } z = x; SET_LOW_WORD(z,0); r = exp(-z*z-0.5625)*exp((z-x)*(z+x)+R/S); if (hx >= 0) - return one-r/x; - return r/x-one; + return 1.0 - r/x; + return r/x - 1.0; } double erfc(double x) @@ -248,49 +245,49 @@ double erfc(double x) ix = hx & 0x7fffffff; if (ix >= 0x7ff00000) { /* 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 < 0x3c700000) /* |x| < 2**-56 */ - return one - x; + return 1.0 - x; z = x*x; r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4))); - s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5)))); + s = 1.0+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5)))); y = r/s; if (hx < 0x3fd00000) { /* x < 1/4 */ - return one - (x+x*y); + return 1.0 - (x+x*y); } else { r = x*y; - r += x-half; - return half - r ; + r += x - 0.5; + return 0.5 - r ; } } 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))))); - 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) { - z = one-erx; + z = 1.0-erx; return z - P/Q; } else { z = erx+P/Q; - return one+z; + return 1.0 + z; } } if (ix < 0x403c0000) { /* |x| < 28 */ x = fabs(x); - s = one/(x*x); + s = 1.0/(x*x); if (ix < 0x4006DB6D) { /* |x| < 1/.35 ~ 2.857143*/ R = ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*( ra5+s*(ra6+s*ra7)))))); - S = 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))))))); } else { /* |x| >= 1/.35 ~ 2.857143 */ 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*( 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)))))); } z = x; @@ -298,9 +295,9 @@ double erfc(double x) r = exp(-z*z-0.5625)*exp((z-x)*(z+x)+R/S); if (hx > 0) return r/x; - return two-r/x; + return 2.0 - r/x; } if (hx > 0) return tiny*tiny; - return two-tiny; + return 2.0 - tiny; } diff --git a/src/math/erff.c b/src/math/erff.c index eef4851e..a0a304e5 100644 --- a/src/math/erff.c +++ b/src/math/erff.c @@ -17,9 +17,6 @@ static const float tiny = 1e-30, -half = 5.0000000000e-01, /* 0x3F000000 */ -one = 1.0000000000e+00, /* 0x3F800000 */ -two = 2.0000000000e+00, /* 0x40000000 */ /* c = (subfloat)0.84506291151 */ erx = 8.4506291151e-01, /* 0x3f58560b */ /* @@ -100,7 +97,7 @@ float erff(float x) if (ix >= 0x7f800000) { /* erf(nan)=nan, erf(+-inf)=+-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 < 0x31800000) { /* |x| < 2**-28 */ @@ -111,42 +108,42 @@ float erff(float x) } z = x*x; r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4))); - s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5)))); + s = 1.0f+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5)))); y = r/s; return x + x*y; } if (ix < 0x3fa00000) { /* 0.84375 <= |x| < 1.25 */ - s = fabsf(x)-one; + s = fabsf(x)-1.0f; 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) return erx + P/Q; return -erx - P/Q; } if (ix >= 0x40c00000) { /* inf > |x| >= 6 */ if (hx >= 0) - return one - tiny; - return tiny - one; + return 1.0f - tiny; + return tiny - 1.0f; } x = fabsf(x); - s = one/(x*x); + s = 1.0f/(x*x); if (ix < 0x4036DB6E) { /* |x| < 1/0.35 */ R = ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*( ra5+s*(ra6+s*ra7)))))); - S = one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*( + S = 1.0f+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*( sa5+s*(sa6+s*(sa7+s*sa8))))))); } else { /* |x| >= 1/0.35 */ R = rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*( rb5+s*rb6))))); - S = one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*( + S = 1.0f+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*( sb5+s*(sb6+s*sb7)))))); } GET_FLOAT_WORD(ix, x); SET_FLOAT_WORD(z, ix&0xfffff000); r = expf(-z*z - 0.5625f) * expf((z-x)*(z+x) + R/S); if (hx >= 0) - return one - r/x; - return r/x - one; + return 1.0f - r/x; + return r/x - 1.0f; } float erfcf(float x) @@ -158,50 +155,50 @@ float erfcf(float x) ix = hx & 0x7fffffff; if (ix >= 0x7f800000) { /* 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 < 0x23800000) /* |x| < 2**-56 */ - return one - x; + return 1.0f - x; z = x*x; r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4))); - s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5)))); + s = 1.0f+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5)))); y = r/s; if (hx < 0x3e800000) { /* x<1/4 */ - return one - (x+x*y); + return 1.0f - (x+x*y); } else { r = x*y; - r += (x-half); - return half - r ; + r += (x-0.5f); + return 0.5f - r ; } } 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))))); - 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) { - z = one - erx; + z = 1.0f - erx; return z - P/Q; } else { z = erx + P/Q; - return one + z; + return 1.0f + z; } } if (ix < 0x41e00000) { /* |x| < 28 */ x = fabsf(x); - s = one/(x*x); + s = 1.0f/(x*x); if (ix < 0x4036DB6D) { /* |x| < 1/.35 ~ 2.857143*/ R = ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*( ra5+s*(ra6+s*ra7)))))); - S = 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))))))); } else { /* |x| >= 1/.35 ~ 2.857143 */ 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*( 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)))))); } 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); if (hx > 0) return r/x; - return two - r/x; + return 2.0f - r/x; } if (hx > 0) return tiny*tiny; - return two - tiny; + return 2.0f - tiny; } diff --git a/src/math/erfl.c b/src/math/erfl.c index a80c2ce1..f7449859 100644 --- a/src/math/erfl.c +++ b/src/math/erfl.c @@ -108,9 +108,6 @@ long double erfl(long double x) #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 static const long double tiny = 1e-4931L, -half = 0.5L, -one = 1.0L, -two = 2.0L, /* c = (float)0.84506291151 */ erx = 0.845062911510467529296875L, @@ -248,12 +245,12 @@ long double erfl(long double x) int32_t ix, i; uint32_t se, i0, i1; - GET_LDOUBLE_WORDS (se, i0, i1, x); + GET_LDOUBLE_WORDS(se, i0, i1, x); ix = se & 0x7fff; if (ix >= 0x7fff) { /* erf(nan)=nan */ 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); @@ -272,7 +269,7 @@ long double erfl(long double x) return x + x * y; } 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] + s * (pa[3] + s * (pa[4] + s * (pa[5] + s * (pa[6] + s * pa[7])))))); 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 ((se & 0x8000) == 0) - return one - tiny; - return tiny - one; + return 1.0 - tiny; + return tiny - 1.0; } x = fabsl (x); - s = one / (x * x); + s = 1.0 / (x * x); if (ix < 0x4000b6db) { /* 1.25 <= |x| < 2.85711669921875 ~ 1/.35 */ R = ra[0] + s * (ra[1] + s * (ra[2] + s * (ra[3] + s * (ra[4] + s * (ra[5] + s * (ra[6] + s * (ra[7] + s * ra[8]))))))); @@ -305,8 +302,8 @@ long double erfl(long double x) SET_LDOUBLE_WORDS(z, i, i0, i1); r = expl(-z * z - 0.5625) * expl((z - x) * (z + x) + R / S); if ((se & 0x8000) == 0) - return one - r / x; - return r / x - one; + return 1.0 - r / x; + return r / x - 1.0; } long double erfcl(long double x) @@ -315,16 +312,16 @@ long double erfcl(long double x) long double R, S, P, Q, s, y, z, r; uint32_t se, i0, i1; - GET_LDOUBLE_WORDS (se, i0, i1, x); + GET_LDOUBLE_WORDS(se, i0, i1, x); ix = se & 0x7fff; if (ix >= 0x7fff) { /* erfc(nan) = nan, erfc(+-inf) = 0,2 */ - return (long double)(((se & 0xffff) >> 15) << 1) + one / x; + return (long double)(((se & 0xffff) >> 15) << 1) + 1.0 / x; } ix = (ix << 16) | (i0 >> 16); if (ix < 0x3ffed800) { /* |x| < 0.84375 */ if (ix < 0x3fbe0000) /* |x| < 2**-65 */ - return one - x; + return 1.0 - x; z = x * x; r = pp[0] + z * (pp[1] + 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))))); y = r / s; if (ix < 0x3ffd8000) /* x < 1/4 */ - return one - (x + x * y); + return 1.0 - (x + x * y); r = x * y; - r += x - half; - return half - r; + r += x - 0.5L; + return 0.5L - r; } 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] + s * (pa[3] + s * (pa[4] + s * (pa[5] + s * (pa[6] + s * pa[7])))))); Q = qa[0] + s * (qa[1] + s * (qa[2] + s * (qa[3] + s * (qa[4] + s * (qa[5] + s * (qa[6] + s)))))); if ((se & 0x8000) == 0) { - z = one - erx; + z = 1.0 - erx; return z - P / Q; } z = erx + P / Q; - return one + z; + return 1.0 + z; } if (ix < 0x4005d600) { /* |x| < 107 */ - x = fabsl (x); - s = one / (x * x); + x = fabsl(x); + s = 1.0 / (x * x); if (ix < 0x4000b6db) { /* 1.25 <= |x| < 2.85711669921875 ~ 1/.35 */ R = ra[0] + s * (ra[1] + s * (ra[2] + s * (ra[3] + s * (ra[4] + s * (ra[5] + s * (ra[6] + s * (ra[7] + s * ra[8]))))))); @@ -365,26 +362,25 @@ long double erfcl(long double x) s * (sb[5] + s * (sb[6] + s)))))); } else { /* 107 > |x| >= 6.666 */ 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] + s * (rc[4] + s * rc[5])))); S = sc[0] + s * (sc[1] + s * (sc[2] + s * (sc[3] + s * (sc[4] + s)))); } z = x; - GET_LDOUBLE_WORDS (hx, i0, i1, z); + GET_LDOUBLE_WORDS(hx, i0, i1, z); i1 = 0; i0 &= 0xffffff00; - SET_LDOUBLE_WORDS (z, hx, i0, i1); - r = expl (-z * z - 0.5625) * - expl ((z - x) * (z + x) + R / S); + SET_LDOUBLE_WORDS(z, hx, i0, i1); + r = expl(-z * z - 0.5625) * expl((z - x) * (z + x) + R / S); if ((se & 0x8000) == 0) return r / x; - return two - r / x; + return 2.0 - r / x; } if ((se & 0x8000) == 0) return tiny * tiny; - return two - tiny; + return 2.0 - tiny; } #endif diff --git a/src/math/exp.c b/src/math/exp.c index a538b8cd..29bf9609 100644 --- a/src/math/exp.c +++ b/src/math/exp.c @@ -74,7 +74,6 @@ #include "libm.h" static const double -one = 1.0, halF[2] = {0.5,-0.5,}, huge = 1.0e+300, o_threshold = 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */ @@ -134,8 +133,8 @@ double exp(double x) STRICT_ASSIGN(double, x, hi - lo); } else if(hx < 0x3e300000) { /* |x| < 2**-28 */ /* raise inexact */ - if (huge+x > one) - return one+x; + if (huge+x > 1.0) + return 1.0+x; } else k = 0; @@ -147,8 +146,8 @@ double exp(double x) INSERT_WORDS(twopk, 0x3ff00000+((k+1000)<<20), 0); c = x - t*(P1+t*(P2+t*(P3+t*(P4+t*P5)))); if (k == 0) - return one - ((x*c)/(c-2.0) - x); - y = one-((lo-(x*c)/(2.0-c))-hi); + return 1.0 - ((x*c)/(c-2.0) - x); + y = 1.0-((lo-(x*c)/(2.0-c))-hi); if (k < -1021) return y*twopk*twom1000; if (k == 1024) diff --git a/src/math/expf.c b/src/math/expf.c index f706ac5d..3c3e2ab9 100644 --- a/src/math/expf.c +++ b/src/math/expf.c @@ -16,7 +16,6 @@ #include "libm.h" static const float -one = 1.0, halF[2] = {0.5,-0.5,}, huge = 1.0e+30, o_threshold = 8.8721679688e+01, /* 0x42b17180 */ @@ -72,8 +71,8 @@ float expf(float x) STRICT_ASSIGN(float, x, hi - lo); } else if(hx < 0x39000000) { /* |x|<2**-14 */ /* raise inexact */ - if (huge+x > one) - return one + x; + if (huge+x > 1.0f) + return 1.0f + x; } else k = 0; @@ -85,8 +84,8 @@ float expf(float x) SET_FLOAT_WORD(twopk, 0x3f800000+((k+100)<<23)); c = x - t*(P1+t*P2); if (k == 0) - return one - ((x*c)/(c - 2.0f) - x); - y = one - ((lo - (x*c)/(2.0f - c)) - hi); + return 1.0f - ((x*c)/(c - 2.0f) - x); + y = 1.0f - ((lo - (x*c)/(2.0f - c)) - hi); if (k < -125) return y*twopk*twom100; if (k == 128) diff --git a/src/math/expl.c b/src/math/expl.c index 9507fd2e..b289ffec 100644 --- a/src/math/expl.c +++ b/src/math/expl.c @@ -102,13 +102,13 @@ long double expl(long double x) if (x > MAXLOGL) return INFINITY; if (x < MINLOGL) - return 0.0L; + return 0.0; /* Express e**x = e**g 2**n * = e**g e**(n loge(2)) * = e**(g + n loge(2)) */ - px = floorl(LOG2EL * x + 0.5L); /* floor() truncates toward -infinity. */ + px = floorl(LOG2EL * x + 0.5); /* floor() truncates toward -infinity. */ n = px; x -= px * C1; x -= px * C2; @@ -120,8 +120,8 @@ long double expl(long double x) xx = x * x; px = x * __polevll(xx, P, 2); x = px/(__polevll(xx, Q, 3) - px); - x = 1.0L + ldexpl(x, 1); - x = ldexpl(x, n); + x = 1.0 + 2.0 * x; + x = scalbnl(x, n); return x; } #endif diff --git a/src/math/expm1.c b/src/math/expm1.c index ffa82264..f8f32c46 100644 --- a/src/math/expm1.c +++ b/src/math/expm1.c @@ -107,7 +107,6 @@ #include "libm.h" static const double -one = 1.0, huge = 1.0e+300, tiny = 1.0e-300, 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 */ /* raise inexact */ 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 */ hfx = 0.5*x; 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; e = hxs*((r1-t)/(6.0 - x*t)); if (k == 0) /* c is 0 */ @@ -195,17 +194,17 @@ double expm1(double x) if (k == 1) { if (x < -0.25) 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 */ - y = one - (e-x); + y = 1.0 - (e-x); if (k == 1024) y = y*2.0*0x1p1023; else y = y*twopk; - return y - one; + return y - 1.0; } - t = one; + t = 1.0; if (k < 20) { SET_HIGH_WORD(t, 0x3ff00000 - (0x200000>>k)); /* t=1-2^-k */ y = t-(e-x); @@ -213,7 +212,7 @@ double expm1(double x) } else { SET_HIGH_WORD(t, ((0x3ff-k)<<20)); /* 2^-k */ y = x-(e+t); - y += one; + y += 1.0; y = y*twopk; } return y; diff --git a/src/math/expm1f.c b/src/math/expm1f.c index a8b79e88..d9568466 100644 --- a/src/math/expm1f.c +++ b/src/math/expm1f.c @@ -16,7 +16,6 @@ #include "libm.h" static const float -one = 1.0, huge = 1.0e+30, tiny = 1.0e-30, o_threshold = 8.8721679688e+01, /* 0x42b17180 */ @@ -54,7 +53,7 @@ float expm1f(float x) if (xsb != 0) { /* x < -27*ln2 */ /* raise inexact */ 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 */ hfx = 0.5f*x; hxs = x*hfx; - r1 = one+hxs*(Q1+hxs*Q2); + r1 = 1.0f+hxs*(Q1+hxs*Q2); t = 3.0f - r1*hfx; e = hxs*((r1-t)/(6.0f - x*t)); if (k == 0) /* c is 0 */ @@ -100,17 +99,17 @@ float expm1f(float x) if (k == 1) { if (x < -0.25f) 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 */ - y = one - (e - x); + y = 1.0f - (e - x); if (k == 128) y = y*2.0f*0x1p127f; else y = y*twopk; - return y - one; + return y - 1.0f; } - t = one; + t = 1.0f; if (k < 23) { SET_FLOAT_WORD(t, 0x3f800000 - (0x1000000>>k)); /* t=1-2^-k */ y = t - (e - x); @@ -118,7 +117,7 @@ float expm1f(float x) } else { SET_FLOAT_WORD(t, (0x7f-k)<<23); /* 2^-k */ y = x - (e + t); - y += one; + y += 1.0f; y = y*twopk; } return y; diff --git a/src/math/expm1l.c b/src/math/expm1l.c index 2f94dfa2..ca0d141f 100644 --- a/src/math/expm1l.c +++ b/src/math/expm1l.c @@ -97,11 +97,11 @@ long double expm1l(long double x) return x; /* Minimum value.*/ if (x < minarg) - return -1.0L; + return -1.0; xx = C1 + C2; /* Express x = ln 2 (k + remainder), remainder not exceeding 1/2. */ - px = floorl (0.5 + x / xx); + px = floorl(0.5 + x / xx); k = px; /* remainder times ln 2 */ x -= px * C1; @@ -116,7 +116,7 @@ long double expm1l(long double x) /* exp(x) = exp(k ln 2) exp(remainder ln 2) = 2^k exp(remainder ln 2). We have qx = exp(remainder ln 2) - 1, so exp(x) - 1 = 2^k (qx + 1) - 1 = 2^k qx + 2^k - 1. */ - px = ldexpl(1.0L, k); + px = scalbnl(1.0, k); x = px * qx + (px - 1.0); return x; } diff --git a/src/math/fma.c b/src/math/fma.c index 87d450c7..5fb95406 100644 --- a/src/math/fma.c +++ b/src/math/fma.c @@ -247,7 +247,7 @@ static inline double add_and_denormalize(double a, double b, int scale) INSERT_WORD64(sum.hi, hibits); } } - return (ldexp(sum.hi, scale)); + return scalbn(sum.hi, scale); } /* @@ -364,7 +364,7 @@ double fma(double x, double y, double z) } } if (spread <= DBL_MANT_DIG * 2) - zs = ldexp(zs, -spread); + zs = scalbn(zs, -spread); else zs = copysign(DBL_MIN, zs); @@ -390,7 +390,7 @@ double fma(double x, double y, double z) */ fesetround(oround); volatile double vzs = zs; /* XXX gcc CSE bug workaround */ - return (xy.hi + vzs + ldexp(xy.lo, spread)); + return xy.hi + vzs + scalbn(xy.lo, spread); } if (oround != FE_TONEAREST) { @@ -400,13 +400,13 @@ double fma(double x, double y, double z) */ fesetround(oround); adj = r.lo + xy.lo; - return (ldexp(r.hi + adj, spread)); + return scalbn(r.hi + adj, spread); } adj = add_adjusted(r.lo, xy.lo); if (spread + ilogb(r.hi) > -1023) - return (ldexp(r.hi + adj, spread)); + return scalbn(r.hi + adj, spread); else - return (add_and_denormalize(r.hi, adj, spread)); + return add_and_denormalize(r.hi, adj, spread); } #endif diff --git a/src/math/fmal.c b/src/math/fmal.c index cbaf46eb..be64f145 100644 --- a/src/math/fmal.c +++ b/src/math/fmal.c @@ -115,7 +115,7 @@ static inline long double add_and_denormalize(long double a, long double b, int if (bits_lost != 1 ^ (int)(u.bits.manl & 1)) sum.hi = nextafterl(sum.hi, INFINITY * sum.lo); } - return (ldexp(sum.hi, scale)); + return scalbnl(sum.hi, scale); } /* @@ -228,7 +228,7 @@ long double fmal(long double x, long double y, long double z) } } if (spread <= LDBL_MANT_DIG * 2) - zs = ldexpl(zs, -spread); + zs = scalbnl(zs, -spread); else zs = copysignl(LDBL_MIN, zs); @@ -254,7 +254,7 @@ long double fmal(long double x, long double y, long double z) */ fesetround(oround); volatile long double vzs = zs; /* XXX gcc CSE bug workaround */ - return (xy.hi + vzs + ldexpl(xy.lo, spread)); + return xy.hi + vzs + scalbnl(xy.lo, spread); } if (oround != FE_TONEAREST) { @@ -264,13 +264,13 @@ long double fmal(long double x, long double y, long double z) */ fesetround(oround); adj = r.lo + xy.lo; - return (ldexpl(r.hi + adj, spread)); + return scalbnl(r.hi + adj, spread); } adj = add_adjusted(r.lo, xy.lo); if (spread + ilogbl(r.hi) > -16383) - return (ldexpl(r.hi + adj, spread)); + return scalbnl(r.hi + adj, spread); else - return (add_and_denormalize(r.hi, adj, spread)); + return add_and_denormalize(r.hi, adj, spread); } #endif diff --git a/src/math/fmod.c b/src/math/fmod.c index 6856844e..84a1b4ac 100644 --- a/src/math/fmod.c +++ b/src/math/fmod.c @@ -17,7 +17,7 @@ #include "libm.h" -static const double one = 1.0, Zero[] = {0.0, -0.0,}; +static const double Zero[] = {0.0, -0.0,}; double fmod(double x, double y) { @@ -140,7 +140,6 @@ double fmod(double x, double y) lx = hx>>(n-32); hx = sx; } INSERT_WORDS(x, hx|sx, lx); - x *= one; /* create necessary signal */ } return x; /* exact output */ } diff --git a/src/math/fmodf.c b/src/math/fmodf.c index 4b50a3d3..837e9371 100644 --- a/src/math/fmodf.c +++ b/src/math/fmodf.c @@ -20,7 +20,7 @@ #include "libm.h" -static const float one = 1.0, Zero[] = {0.0, -0.0,}; +static const float Zero[] = {0.0, -0.0,}; float fmodf(float x, float y) { @@ -99,7 +99,6 @@ float fmodf(float x, float y) n = -126 - iy; hx >>= n; SET_FLOAT_WORD(x, hx|sx); - x *= one; /* create necessary signal */ } return x; /* exact output */ } diff --git a/src/math/fmodl.c b/src/math/fmodl.c index 2e3eec1f..b930c49d 100644 --- a/src/math/fmodl.c +++ b/src/math/fmodl.c @@ -48,7 +48,7 @@ typedef uint32_t manh_t; #define MANL_SHIFT (LDBL_MANL_SIZE - 1) -static const long double one = 1.0, Zero[] = {0.0, -0.0,}; +static const long double Zero[] = {0.0, -0.0,}; /* * fmodl(x,y) @@ -153,7 +153,6 @@ long double fmodl(long double x, long double y) } else { ux.bits.exp = iy + BIAS; } - x = ux.e * one; /* create necessary signal */ - return x; /* exact output */ + return ux.e; /* exact output */ } #endif diff --git a/src/math/j0.c b/src/math/j0.c index b5490641..986968e8 100644 --- a/src/math/j0.c +++ b/src/math/j0.c @@ -60,7 +60,6 @@ static double pzero(double), qzero(double); static const double huge = 1e300, -one = 1.0, invsqrtpi = 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */ tpi = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */ /* R0/S0 on [0, 2.00] */ @@ -73,8 +72,6 @@ S02 = 1.16926784663337450260e-04, /* 0x3F1EA6D2, 0xDD57DBF4 */ S03 = 5.13546550207318111446e-07, /* 0x3EA13B54, 0xCE84D5A9 */ S04 = 1.16614003333790000205e-09; /* 0x3E1408BC, 0xF4745D8F */ -static const double zero = 0.0; - double j0(double x) { double z, s,c,ss,cc,r,u,v; @@ -83,7 +80,7 @@ double j0(double x) GET_HIGH_WORD(hx, x); ix = hx & 0x7fffffff; if (ix >= 0x7ff00000) - return one/(x*x); + return 1.0/(x*x); x = fabs(x); if (ix >= 0x40000000) { /* |x| >= 2.0 */ s = sin(x); @@ -92,7 +89,7 @@ double j0(double x) cc = s+c; if (ix < 0x7fe00000) { /* make sure x+x does not overflow */ z = -cos(x+x); - if ((s*c) < zero) + if (s*c < 0.0) cc = z/ss; else ss = z/cc; @@ -112,20 +109,20 @@ double j0(double x) } if (ix < 0x3f200000) { /* |x| < 2**-13 */ /* raise inexact if x != 0 */ - if (huge+x > one) { + if (huge+x > 1.0) { if (ix < 0x3e400000) /* |x| < 2**-27 */ - return one; - return one - 0.25*x*x; + return 1.0; + return 1.0 - 0.25*x*x; } } z = x*x; 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 */ - return one + z*(-0.25+(r/s)); + return 1.0 + z*(-0.25+(r/s)); } else { 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; /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0 */ if (ix >= 0x7ff00000) - return one/(x+x*x); + return 1.0/(x+x*x); if ((ix|lx) == 0) - return -one/zero; + return -1.0/0.0; if (hx < 0) - return zero/zero; + return 0.0/0.0; if (ix >= 0x40000000) { /* |x| >= 2.0 */ /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0)) * where x0 = x-pi/4 @@ -178,7 +175,7 @@ double y0(double x) */ if (ix < 0x7fe00000) { /* make sure x+x does not overflow */ z = -cos(x+x); - if (s*c < zero) + if (s*c < 0.0) cc = z/ss; else ss = z/cc; @@ -197,7 +194,7 @@ double y0(double x) } z = x*x; 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)); } @@ -286,10 +283,10 @@ static double pzero(double x) else if (ix >= 0x40122E8B){p = pR5; q = pS5;} else if (ix >= 0x4006DB6D){p = pR3; q = pS3;} 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])))); - s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); - return one + r/s; + s = 1.0+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); + 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 >= 0x4006DB6D){p = qR3; q = qS3;} 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])))); - 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; } diff --git a/src/math/j0f.c b/src/math/j0f.c index 52cb0ab4..2ee2824b 100644 --- a/src/math/j0f.c +++ b/src/math/j0f.c @@ -19,7 +19,6 @@ static float pzerof(float), qzerof(float); static const float huge = 1e30, -one = 1.0, invsqrtpi = 5.6418961287e-01, /* 0x3f106ebb */ tpi = 6.3661974669e-01, /* 0x3f22f983 */ /* R0/S0 on [0, 2.00] */ @@ -32,8 +31,6 @@ S02 = 1.1692678527e-04, /* 0x38f53697 */ S03 = 5.1354652442e-07, /* 0x3509daa6 */ S04 = 1.1661400734e-09; /* 0x30a045e8 */ -static const float zero = 0.0; - float j0f(float x) { float z, s,c,ss,cc,r,u,v; @@ -42,7 +39,7 @@ float j0f(float x) GET_FLOAT_WORD(hx, x); ix = hx & 0x7fffffff; if (ix >= 0x7f800000) - return one/(x*x); + return 1.0f/(x*x); x = fabsf(x); if (ix >= 0x40000000) { /* |x| >= 2.0 */ s = sinf(x); @@ -51,7 +48,7 @@ float j0f(float x) cc = s+c; if (ix < 0x7f000000) { /* make sure x+x does not overflow */ z = -cosf(x+x); - if (s*c < zero) + if (s*c < 0.0f) cc = z/ss; else ss = z/cc; @@ -71,20 +68,20 @@ float j0f(float x) } if (ix < 0x39000000) { /* |x| < 2**-13 */ /* raise inexact if x != 0 */ - if (huge+x > one) { + if (huge+x > 1.0f) { if (ix < 0x32000000) /* |x| < 2**-27 */ - return one; - return one - 0.25f*x*x; + return 1.0f; + return 1.0f - 0.25f*x*x; } } z = x*x; 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 */ - return one + z*(-0.25f + (r/s)); + return 1.0f + z*(-0.25f + (r/s)); } else { 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; /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0 */ if (ix >= 0x7f800000) - return one/(x+x*x); + return 1.0f/(x+x*x); if (ix == 0) - return -one/zero; + return -1.0f/0.0f; if (hx < 0) - return zero/zero; + return 0.0f/0.0f; if (ix >= 0x40000000) { /* |x| >= 2.0 */ /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0)) * where x0 = x-pi/4 @@ -137,7 +134,7 @@ float y0f(float x) */ if (ix < 0x7f000000) { /* make sure x+x not overflow */ z = -cosf(x+x); - if (s*c < zero) + if (s*c < 0.0f) cc = z/ss; else ss = z/cc; @@ -156,7 +153,7 @@ float y0f(float x) } z = x*x; 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)); } @@ -244,10 +241,10 @@ static float pzerof(float x) else if (ix >= 0x40f71c58){p = pR5; q = pS5;} else if (ix >= 0x4036db68){p = pR3; q = pS3;} 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])))); - s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); - return one + r/s; + s = 1.0f+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); + 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 >= 0x4036db68){p = qR3; q = qS3;} 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])))); - 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; } diff --git a/src/math/j1.c b/src/math/j1.c index 29ccff0c..4a2cba8d 100644 --- a/src/math/j1.c +++ b/src/math/j1.c @@ -60,7 +60,6 @@ static double pone(double), qone(double); static const double huge = 1e300, -one = 1.0, invsqrtpi = 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */ tpi = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */ /* R0/S0 on [0,2] */ @@ -74,8 +73,6 @@ s03 = 1.17718464042623683263e-06, /* 0x3EB3BFF8, 0x333F8498 */ s04 = 5.04636257076217042715e-09, /* 0x3E35AC88, 0xC97DFF2C */ s05 = 1.23542274426137913908e-11; /* 0x3DAB2ACF, 0xCFB97ED8 */ -static const double zero = 0.0; - double j1(double x) { double z,s,c,ss,cc,r,u,v,y; @@ -84,7 +81,7 @@ double j1(double x) GET_HIGH_WORD(hx, x); ix = hx & 0x7fffffff; if (ix >= 0x7ff00000) - return one/x; + return 1.0/x; y = fabs(x); if (ix >= 0x40000000) { /* |x| >= 2.0 */ s = sin(y); @@ -93,7 +90,7 @@ double j1(double x) cc = s-c; if (ix < 0x7fe00000) { /* make sure y+y not overflow */ z = cos(y+y); - if (s*c > zero) + if (s*c > 0.0) cc = z/ss; else ss = z/cc; @@ -116,12 +113,12 @@ double j1(double x) } if (ix < 0x3e400000) { /* |x| < 2**-27 */ /* raise inexact if x!=0 */ - if (huge+x > one) + if (huge+x > 1.0) return 0.5*x; } z = x*x; 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; return x*0.5 + r/s; } @@ -151,11 +148,11 @@ double y1(double x) ix = 0x7fffffff & hx; /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */ if (ix >= 0x7ff00000) - return one/(x+x*x); + return 1.0/(x+x*x); if ((ix|lx) == 0) - return -one/zero; + return -1.0/0.0; if (hx < 0) - return zero/zero; + return 0.0/0.0; if (ix >= 0x40000000) { /* |x| >= 2.0 */ s = sin(x); c = cos(x); @@ -163,7 +160,7 @@ double y1(double x) cc = s-c; if (ix < 0x7fe00000) { /* make sure x+x not overflow */ z = cos(x+x); - if (s*c > zero) + if (s*c > 0.0) cc = z/ss; else ss = z/cc; @@ -192,8 +189,8 @@ double y1(double x) return -tpi/x; z = x*x; 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])))); - return x*(u/v) + tpi*(j1(x)*log(x)-one/x); + 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)-1.0/x); } /* 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 >= 0x4006DB6D){p = pr3; q = ps3;} 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])))); - s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); - return one+ r/s; + s = 1.0+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); + return 1.0+ r/s; } /* 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 >= 0x4006DB6D){p = qr3; q = qs3;} 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])))); - 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; } diff --git a/src/math/j1f.c b/src/math/j1f.c index 2d617b67..de459e0d 100644 --- a/src/math/j1f.c +++ b/src/math/j1f.c @@ -19,7 +19,6 @@ static float ponef(float), qonef(float); static const float huge = 1e30, -one = 1.0, invsqrtpi = 5.6418961287e-01, /* 0x3f106ebb */ tpi = 6.3661974669e-01, /* 0x3f22f983 */ /* R0/S0 on [0,2] */ @@ -33,8 +32,6 @@ s03 = 1.1771846857e-06, /* 0x359dffc2 */ s04 = 5.0463624390e-09, /* 0x31ad6446 */ s05 = 1.2354227016e-11; /* 0x2d59567e */ -static const float zero = 0.0; - float j1f(float x) { float z,s,c,ss,cc,r,u,v,y; @@ -43,7 +40,7 @@ float j1f(float x) GET_FLOAT_WORD(hx, x); ix = hx & 0x7fffffff; if (ix >= 0x7f800000) - return one/x; + return 1.0f/x; y = fabsf(x); if (ix >= 0x40000000) { /* |x| >= 2.0 */ s = sinf(y); @@ -52,7 +49,7 @@ float j1f(float x) cc = s-c; if (ix < 0x7f000000) { /* make sure y+y not overflow */ z = cosf(y+y); - if (s*c > zero) + if (s*c > 0.0f) cc = z/ss; else ss = z/cc; @@ -74,12 +71,12 @@ float j1f(float x) } if (ix < 0x32000000) { /* |x| < 2**-27 */ /* raise inexact if x!=0 */ - if (huge+x > one) + if (huge+x > 1.0f) return 0.5f*x; } z = x*x; 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; return 0.5f*x + r/s; } @@ -108,11 +105,11 @@ float y1f(float x) ix = 0x7fffffff & hx; /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */ if (ix >= 0x7f800000) - return one/(x+x*x); + return 1.0f/(x+x*x); if (ix == 0) - return -one/zero; + return -1.0f/0.0f; if (hx < 0) - return zero/zero; + return 0.0f/0.0f; if (ix >= 0x40000000) { /* |x| >= 2.0 */ s = sinf(x); c = cosf(x); @@ -120,7 +117,7 @@ float y1f(float x) cc = s-c; if (ix < 0x7f000000) { /* make sure x+x not overflow */ z = cosf(x+x); - if (s*c > zero) + if (s*c > 0.0f) cc = z/ss; else ss = z/cc; @@ -149,8 +146,8 @@ float y1f(float x) return -tpi/x; z = x*x; 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])))); - return x*(u/v) + tpi*(j1f(x)*logf(x)-one/x); + 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)-1.0f/x); } /* 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 >= 0x4036db68){p = pr3; q = ps3;} 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])))); - s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); - return one + r/s; + s = 1.0f+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); + return 1.0f + r/s; } /* 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 >= 0x4036db68){p = qr3; q = qs3;} 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])))); - 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; } diff --git a/src/math/jn.c b/src/math/jn.c index 082a17bc..d95af15d 100644 --- a/src/math/jn.c +++ b/src/math/jn.c @@ -37,12 +37,7 @@ #include "libm.h" -static const double -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; +static const double invsqrtpi = 5.64189583547756279280e-01; /* 0x3FE20DD7, 0x50429B6D */ 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) */ x = fabs(x); if ((ix|lx) == 0 || ix >= 0x7ff00000) /* if x is 0 or inf */ - b = zero; + b = 0.0; else if ((double)n <= x) { /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */ if (ix >= 0x52D00000) { /* x > 2**302 */ @@ -108,11 +103,11 @@ double jn(int n, double x) * J(n,x) = 1/n!*(x/2)^n - ... */ if (n > 33) /* underflow */ - b = zero; + b = 0.0; else { temp = x*0.5; b = temp; - for (a=one,i=2; i<=n; i++) { + for (a=1.0,i=2; i<=n; i++) { a *= (double)i; /* a = n! */ b *= temp; /* b = (x/2)^n */ } @@ -165,10 +160,10 @@ double jn(int n, double x) q1 = tmp; } m = n+n; - for (t=zero, i = 2*(n+k); i>=m; i -= 2) - t = one/(i/x-t); + for (t=0.0, i = 2*(n+k); i>=m; i -= 2) + t = 1.0/(i/x-t); a = t; - b = one; + b = 1.0; /* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n) * Hence, if n*(log(2n/x)) > ... * single 8.8722839355e+01 @@ -178,7 +173,7 @@ double jn(int n, double x) * likely underflow to zero */ tmp = n; - v = two/x; + v = 2.0/x; tmp = tmp*log(fabs(v*tmp)); if (tmp < 7.09782712893383973096e+02) { for (i=n-1,di=(double)(i+i); i>0; i--) { @@ -186,7 +181,7 @@ double jn(int n, double x) b *= di; b = b/x - a; a = temp; - di -= two; + di -= 2.0; } } else { for (i=n-1,di=(double)(i+i); i>0; i--) { @@ -194,12 +189,12 @@ double jn(int n, double x) b *= di; b = b/x - a; a = temp; - di -= two; + di -= 2.0; /* scale b to avoid spurious overflow */ if (b > 1e100) { a /= 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) return x+x; if ((ix|lx) == 0) - return -one/zero; + return -1.0/0.0; if (hx < 0) - return zero/zero; + return 0.0/0.0; sign = 1; if (n < 0) { n = -n; @@ -242,7 +237,7 @@ double yn(int n, double x) if (n == 1) return sign*y1(x); if (ix == 0x7ff00000) - return zero; + return 0.0; if (ix >= 0x52D00000) { /* x > 2**302 */ /* (x >> n**2) * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi) diff --git a/src/math/jnf.c b/src/math/jnf.c index b0b36e6b..fd291103 100644 --- a/src/math/jnf.c +++ b/src/math/jnf.c @@ -16,12 +16,6 @@ #define _GNU_SOURCE #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) { 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) */ x = fabsf(x); if (ix == 0 || ix >= 0x7f800000) /* if x is 0 or inf */ - b = zero; + b = 0.0f; else if((float)n <= x) { /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */ a = j0f(x); @@ -63,11 +57,11 @@ float jnf(int n, float x) * J(n,x) = 1/n!*(x/2)^n - ... */ if (n > 33) /* underflow */ - b = zero; + b = 0.0f; else { temp = 0.5f * x; b = temp; - for (a=one,i=2; i<=n; i++) { + for (a=1.0f,i=2; i<=n; i++) { a *= (float)i; /* a = n! */ b *= temp; /* b = (x/2)^n */ } @@ -121,10 +115,10 @@ float jnf(int n, float x) q1 = tmp; } m = n+n; - for (t=zero, i = 2*(n+k); i>=m; i -= 2) - t = one/(i/x-t); + for (t=0.0f, i = 2*(n+k); i>=m; i -= 2) + t = 1.0f/(i/x-t); a = t; - b = one; + b = 1.0f; /* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n) * Hence, if n*(log(2n/x)) > ... * single 8.8722839355e+01 @@ -134,7 +128,7 @@ float jnf(int n, float x) * likely underflow to zero */ tmp = n; - v = two/x; + v = 2.0f/x; tmp = tmp*logf(fabsf(v*tmp)); if (tmp < 88.721679688f) { for (i=n-1,di=(float)(i+i); i>0; i--) { @@ -142,7 +136,7 @@ float jnf(int n, float x) b *= di; b = b/x - a; a = temp; - di -= two; + di -= 2.0f; } } else { for (i=n-1,di=(float)(i+i); i>0; i--){ @@ -150,12 +144,12 @@ float jnf(int n, float x) b *= di; b = b/x - a; a = temp; - di -= two; + di -= 2.0f; /* scale b to avoid spurious overflow */ if (b > 1e10f) { a /= b; t /= b; - b = one; + b = 1.0f; } } } @@ -183,9 +177,9 @@ float ynf(int n, float x) if (ix > 0x7f800000) return x+x; if (ix == 0) - return -one/zero; + return -1.0f/0.0f; if (hx < 0) - return zero/zero; + return 0.0f/0.0f; sign = 1; if (n < 0) { n = -n; @@ -196,7 +190,7 @@ float ynf(int n, float x) if (n == 1) return sign*y1f(x); if (ix == 0x7f800000) - return zero; + return 0.0f; a = y0f(x); b = y1f(x); diff --git a/src/math/lgamma_r.c b/src/math/lgamma_r.c index a8ef1956..e3ed1733 100644 --- a/src/math/lgamma_r.c +++ b/src/math/lgamma_r.c @@ -82,8 +82,6 @@ static const double two52= 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */ -half= 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */ -one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ pi = 3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */ a0 = 7.72156649015328655494e-02, /* 0x3FB3C467, 0xE37DB0C8 */ a1 = 3.22467033424113591611e-01, /* 0x3FD4A34C, 0xC4A60FAD */ @@ -148,8 +146,6 @@ w4 = -5.95187557450339963135e-04, /* 0xBF4380CB, 0x8C0FE741 */ w5 = 8.36339918996282139126e-04, /* 0x3F4B67BA, 0x4CDAD5D1 */ w6 = -1.63092934096575273989e-03; /* 0xBF5AB89D, 0x0B9E43E4 */ -static const double zero = 0.00000000000000000000e+00; - static double sin_pi(double x) { double y,z; @@ -159,7 +155,7 @@ static double sin_pi(double x) ix &= 0x7fffffff; if (ix < 0x3fd00000) - return __sin(pi*x, zero, 0); + return __sin(pi*x, 0.0, 0); y = -x; /* negative x is assumed */ @@ -174,7 +170,7 @@ static double sin_pi(double x) n = (int)(y*4.0); } else { if (ix >= 0x43400000) { - y = zero; /* y must be even */ + y = 0.0; /* y must be even */ n = 0; } else { if (ix < 0x43300000) @@ -186,14 +182,14 @@ static double sin_pi(double x) } } switch (n) { - case 0: y = __sin(pi*y, zero, 0); break; + case 0: y = __sin(pi*y, 0.0, 0); break; 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 4: y = __sin(pi*(one-y), zero, 0); break; + case 4: y = __sin(pi*(1.0-y), 0.0, 0); break; case 5: - case 6: y = -__cos(pi*(y-1.5), zero); break; - default: y = __sin(pi*(y-2.0), zero, 0); break; + case 6: y = -__cos(pi*(y-1.5), 0.0); break; + default: y = __sin(pi*(y-2.0), 0.0, 0); break; } return -y; } @@ -213,7 +209,7 @@ double __lgamma_r(double x, int *signgamp) if (ix >= 0x7ff00000) return x*x; if ((ix|lx) == 0) - return one/zero; + return 1.0/0.0; if (ix < 0x3b900000) { /* |x|<2**-70, return -log(|x|) */ if(hx < 0) { *signgamp = -1; @@ -223,12 +219,12 @@ double __lgamma_r(double x, int *signgamp) } if (hx < 0) { if (ix >= 0x43300000) /* |x|>=2**52, must be -integer */ - return one/zero; + return 1.0/0.0; t = sin_pi(x); - if (t == zero) /* -integer */ - return one/zero; + if (t == 0.0) /* -integer */ + return 1.0/0.0; nadj = log(pi/fabs(t*x)); - if (t < zero) + if (t < 0.0) *signgamp = -1; x = -x; } @@ -241,17 +237,17 @@ double __lgamma_r(double x, int *signgamp) if (ix <= 0x3feccccc) { /* lgamma(x) = lgamma(x+1)-log(x) */ r = -log(x); if (ix >= 0x3FE76944) { - y = one - x; + y = 1.0 - x; i = 0; } else if (ix >= 0x3FCDA661) { - y = x - (tc-one); + y = x - (tc-1.0); i = 1; } else { y = x; i = 2; } } else { - r = zero; + r = 0.0; if (ix >= 0x3FFBB4C3) { /* [1.7316,2] */ y = 2.0 - x; i = 0; @@ -259,7 +255,7 @@ double __lgamma_r(double x, int *signgamp) y = x - tc; i = 1; } else { - y = x - one; + y = x - 1.0; i = 2; } } @@ -282,16 +278,16 @@ double __lgamma_r(double x, int *signgamp) break; case 2: 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; } } else if (ix < 0x40200000) { /* x < 8.0 */ i = (int)x; y = x - (double)i; 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))))); - r = half*y+p/q; - z = one; /* lgamma(1+s) = log(s) + lgamma(s) */ + q = 1.0+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6))))); + r = 0.5*y+p/q; + z = 1.0; /* lgamma(1+s) = log(s) + lgamma(s) */ switch (i) { case 7: z *= y + 6.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 */ t = log(x); - z = one/x; + z = 1.0/x; y = z*z; 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 */ - r = x*(log(x)-one); + r = x*(log(x)-1.0); if (hx < 0) r = nadj - r; return r; diff --git a/src/math/lgammaf_r.c b/src/math/lgammaf_r.c index f1adcf69..976986aa 100644 --- a/src/math/lgammaf_r.c +++ b/src/math/lgammaf_r.c @@ -17,8 +17,6 @@ static const float two23= 8.3886080000e+06, /* 0x4b000000 */ -half= 5.0000000000e-01, /* 0x3f000000 */ -one = 1.0000000000e+00, /* 0x3f800000 */ pi = 3.1415927410e+00, /* 0x40490fdb */ a0 = 7.7215664089e-02, /* 0x3d9e233f */ a1 = 3.2246702909e-01, /* 0x3ea51a66 */ @@ -83,8 +81,6 @@ w4 = -5.9518753551e-04, /* 0xba1c065c */ w5 = 8.3633989561e-04, /* 0x3a5b3dd2 */ w6 = -1.6309292987e-03; /* 0xbad5c4e8 */ -static const float zero = 0.0000000000e+00; - static float sin_pif(float x) { float y,z; @@ -109,7 +105,7 @@ static float sin_pif(float x) n = (int)(y*4.0f); } else { if (ix >= 0x4b800000) { - y = zero; /* y must be even */ + y = 0.0f; /* y must be even */ n = 0; } else { if (ix < 0x4b000000) @@ -125,7 +121,7 @@ static float sin_pif(float x) case 1: case 2: y = __cosdf(pi*(0.5f - y)); break; case 3: - case 4: y = __sindf(pi*(one - y)); break; + case 4: y = __sindf(pi*(1.0f - y)); break; case 5: case 6: y = -__cosdf(pi*(y - 1.5f)); break; default: y = __sindf(pi*(y - 2.0f)); break; @@ -148,7 +144,7 @@ float __lgammaf_r(float x, int *signgamp) if (ix >= 0x7f800000) return x*x; if (ix == 0) - return one/zero; + return 1.0f/0.0f; if (ix < 0x35000000) { /* |x| < 2**-21, return -log(|x|) */ if (hx < 0) { *signgamp = -1; @@ -158,12 +154,12 @@ float __lgammaf_r(float x, int *signgamp) } if (hx < 0) { if (ix >= 0x4b000000) /* |x| >= 2**23, must be -integer */ - return one/zero; + return 1.0f/0.0f; t = sin_pif(x); - if (t == zero) /* -integer */ - return one/zero; + if (t == 0.0f) /* -integer */ + return 1.0f/0.0f; nadj = logf(pi/fabsf(t*x)); - if (t < zero) + if (t < 0.0f) *signgamp = -1; x = -x; } @@ -176,17 +172,17 @@ float __lgammaf_r(float x, int *signgamp) if (ix <= 0x3f666666) { /* lgamma(x) = lgamma(x+1)-log(x) */ r = -logf(x); if (ix >= 0x3f3b4a20) { - y = one - x; + y = 1.0f - x; i = 0; } else if (ix >= 0x3e6d3308) { - y = x - (tc-one); + y = x - (tc-1.0f); i = 1; } else { y = x; i = 2; } } else { - r = zero; + r = 0.0f; if (ix >= 0x3fdda618) { /* [1.7316,2] */ y = 2.0f - x; i = 0; @@ -194,7 +190,7 @@ float __lgammaf_r(float x, int *signgamp) y = x - tc; i = 1; } else { - y = x - one; + y = x - 1.0f; i = 2; } } @@ -217,16 +213,16 @@ float __lgammaf_r(float x, int *signgamp) break; case 2: 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; } } else if (ix < 0x41000000) { /* x < 8.0 */ i = (int)x; y = x - (float)i; 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))))); - r = half*y+p/q; - z = one; /* lgamma(1+s) = log(s) + lgamma(s) */ + q = 1.0f+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6))))); + r = 0.5f*y+p/q; + z = 1.0f; /* lgamma(1+s) = log(s) + lgamma(s) */ switch (i) { case 7: z *= y + 6.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 */ t = logf(x); - z = one/x; + z = 1.0f/x; y = z*z; 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 */ - r = x*(logf(x)-one); + r = x*(logf(x)-1.0f); if (hx < 0) r = nadj - r; return r; diff --git a/src/math/lgammal.c b/src/math/lgammal.c index ec7c9a04..8fae1be8 100644 --- a/src/math/lgammal.c +++ b/src/math/lgammal.c @@ -95,8 +95,6 @@ long double __lgammal_r(long double x, int *sg) } #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 static const long double -half = 0.5L, -one = 1.0L, pi = 3.14159265358979323846264L, two63 = 9.223372036854775808e18L, @@ -200,8 +198,6 @@ w5 = 8.412723297322498080632E-4L, w6 = -1.880801938119376907179E-3L, w7 = 4.885026142432270781165E-3L; -static const long double zero = 0.0L; - static long double sin_pi(long double x) { long double y, z; @@ -226,7 +222,7 @@ static long double sin_pi(long double x) n = (int) (y*4.0); } else { if (ix >= 0x403f8000) { /* 2^64 */ - y = zero; /* y must be even */ + y = 0.0; /* y must be even */ n = 0; } else { if (ix < 0x403e8000) /* 2^63 */ @@ -244,11 +240,11 @@ static long double sin_pi(long double x) break; case 1: case 2: - y = cosl(pi * (half - y)); + y = cosl(pi * (0.5 - y)); break; case 3: case 4: - y = sinl(pi * (one - y)); + y = sinl(pi * (1.0 - y)); break; case 5: case 6: @@ -273,7 +269,7 @@ long double __lgammal_r(long double x, int *sg) { if ((ix | i0 | i1) == 0) { if (se & 0x8000) *sg = -1; - return one / fabsl(x); + return 1.0 / fabsl(x); } ix = (ix << 16) | (i0 >> 16); @@ -291,10 +287,10 @@ long double __lgammal_r(long double x, int *sg) { } if (se & 0x8000) { t = sin_pi (x); - if (t == zero) - return one / fabsl(t); /* -integer */ + if (t == 0.0) + return 1.0 / fabsl(t); /* -integer */ nadj = logl(pi / fabsl(t * x)); - if (t < zero) + if (t < 0.0) *sg = -1; x = -x; } @@ -306,19 +302,19 @@ long double __lgammal_r(long double x, int *sg) { else if (ix < 0x40008000) { /* x < 2.0 */ if (ix <= 0x3ffee666) { /* 8.99993896484375e-1 */ /* lgamma(x) = lgamma(x+1) - log(x) */ - r = -logl (x); + r = -logl(x); if (ix >= 0x3ffebb4a) { /* 7.31597900390625e-1 */ - y = x - one; + y = x - 1.0; i = 0; } else if (ix >= 0x3ffced33) { /* 2.31639862060546875e-1 */ - y = x - (tc - one); + y = x - (tc - 1.0); i = 1; } else { /* x < 0.23 */ y = x; i = 2; } } else { - r = zero; + r = 0.0; if (ix >= 0x3fffdda6) { /* 1.73162841796875 */ /* [1.7316,2] */ y = x - 2.0; @@ -329,7 +325,7 @@ long double __lgammal_r(long double x, int *sg) { i = 1; } else { /* [0.9, 1.23] */ - y = x - one; + y = x - 1.0; i = 2; } } @@ -337,7 +333,7 @@ long double __lgammal_r(long double x, int *sg) { case 0: p1 = a0 + y * (a1 + y * (a2 + y * (a3 + y * (a4 + y * a5)))); 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; case 1: 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: 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))))); - r += (-half * y + p1 / p2); + r += (-0.5 * y + p1 / p2); } } else if (ix < 0x40028000) { /* 8.0 */ /* x < 8.0 */ i = (int)x; - t = zero; + t = 0.0; y = x - (double)i; 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)))))); - r = half * y + p / q; - z = one;/* lgamma(1+s) = log(s) + lgamma(s) */ + r = 0.5 * y + p / q; + z = 1.0;/* lgamma(1+s) = log(s) + lgamma(s) */ switch (i) { case 7: z *= (y + 6.0); /* FALLTHRU */ @@ -370,18 +366,18 @@ long double __lgammal_r(long double x, int *sg) { z *= (y + 3.0); /* FALLTHRU */ case 3: z *= (y + 2.0); /* FALLTHRU */ - r += logl (z); + r += logl(z); break; } } else if (ix < 0x40418000) { /* 2^66 */ /* 8.0 <= x < 2**66 */ - t = logl (x); - z = one / x; + t = logl(x); + z = 1.0 / x; y = z * z; 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 */ - r = x * (logl (x) - one); + r = x * (logl(x) - 1.0); if (se & 0x8000) r = nadj - r; return r; diff --git a/src/math/log.c b/src/math/log.c index 1bb006a3..98051205 100644 --- a/src/math/log.c +++ b/src/math/log.c @@ -74,8 +74,6 @@ Lg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */ Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */ Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */ -static const double zero = 0.0; - double log(double x) { double hfsq,f,s,z,R,w,t1,t2,dk; @@ -87,9 +85,9 @@ double log(double x) k = 0; if (hx < 0x00100000) { /* x < 2**-1022 */ if (((hx&0x7fffffff)|lx) == 0) - return -two54/zero; /* log(+-0)=-inf */ + return -two54/0.0; /* log(+-0)=-inf */ if (hx < 0) - return (x-x)/zero; /* log(-#) = NaN */ + return (x-x)/0.0; /* log(-#) = NaN */ /* subnormal number, scale up x */ k -= 54; x *= two54; @@ -104,9 +102,9 @@ double log(double x) k += i>>20; f = x - 1.0; if ((0x000fffff&(2+hx)) < 3) { /* -2**-20 <= f < 2**-20 */ - if (f == zero) { + if (f == 0.0) { if (k == 0) { - return zero; + return 0.0; } dk = (double)k; return dk*ln2_hi + dk*ln2_lo; diff --git a/src/math/log10.c b/src/math/log10.c index 5422599a..ed65d9be 100644 --- a/src/math/log10.c +++ b/src/math/log10.c @@ -27,8 +27,6 @@ ivln10lo = 2.50829467116452752298e-11, /* 0x3dbb9438, 0xca9aadd5 */ log10_2hi = 3.01029995663611771306e-01, /* 0x3FD34413, 0x509F6000 */ log10_2lo = 3.69423907715893078616e-13; /* 0x3D59FEF3, 0x11F12B36 */ -static const double zero = 0.0; - double log10(double x) { double f,hfsq,hi,lo,r,val_hi,val_lo,w,y,y2; @@ -40,9 +38,9 @@ double log10(double x) k = 0; if (hx < 0x00100000) { /* x < 2**-1022 */ if (((hx&0x7fffffff)|lx) == 0) - return -two54/zero; /* log(+-0)=-inf */ + return -two54/0.0; /* log(+-0)=-inf */ if (hx<0) - return (x-x)/zero; /* log(-#) = NaN */ + return (x-x)/0.0; /* log(-#) = NaN */ /* subnormal number, scale up x */ k -= 54; x *= two54; @@ -51,7 +49,7 @@ double log10(double x) if (hx >= 0x7ff00000) return x+x; if (hx == 0x3ff00000 && lx == 0) - return zero; /* log(1) = +0 */ + return 0.0; /* log(1) = +0 */ k += (hx>>20) - 1023; hx &= 0x000fffff; i = (hx+0x95f64)&0x100000; diff --git a/src/math/log10f.c b/src/math/log10f.c index 7e4ac9a8..e10749b5 100644 --- a/src/math/log10f.c +++ b/src/math/log10f.c @@ -23,8 +23,6 @@ ivln10lo = -3.1689971365e-05, /* 0xb804ead9 */ log10_2hi = 3.0102920532e-01, /* 0x3e9a2080 */ log10_2lo = 7.9034151668e-07; /* 0x355427db */ -static const float zero = 0.0; - float log10f(float x) { float f,hfsq,hi,lo,r,y; @@ -35,9 +33,9 @@ float log10f(float x) k = 0; if (hx < 0x00800000) { /* x < 2**-126 */ if ((hx&0x7fffffff) == 0) - return -two25/zero; /* log(+-0)=-inf */ + return -two25/0.0f; /* log(+-0)=-inf */ if (hx < 0) - return (x-x)/zero; /* log(-#) = NaN */ + return (x-x)/0.0f; /* log(-#) = NaN */ /* subnormal number, scale up x */ k -= 25; x *= two25; @@ -46,7 +44,7 @@ float log10f(float x) if (hx >= 0x7f800000) return x+x; if (hx == 0x3f800000) - return zero; /* log(1) = +0 */ + return 0.0f; /* log(1) = +0 */ k += (hx>>23) - 127; hx &= 0x007fffff; i = (hx+(0x4afb0d))&0x800000; diff --git a/src/math/log10l.c b/src/math/log10l.c index b954cc77..f0eeeafb 100644 --- a/src/math/log10l.c +++ b/src/math/log10l.c @@ -123,9 +123,9 @@ long double log10l(long double x) if (isnan(x)) return x; - if(x <= 0.0L) { - if(x == 0.0L) - return -1.0L / (x - x); + if(x <= 0.0) { + if(x == 0.0) + return -1.0 / (x - x); return (x - x) / (x - x); } if (x == INFINITY) @@ -142,12 +142,12 @@ long double log10l(long double x) if (e > 2 || e < -2) { if (x < SQRTH) { /* 2(2x-1)/(2x+1) */ e -= 1; - z = x - 0.5L; - y = 0.5L * z + 0.5L; + z = x - 0.5; + y = 0.5 * z + 0.5; } else { /* 2 (x-1)/(x+1) */ - z = x - 0.5L; - z -= 0.5L; - y = 0.5L * x + 0.5L; + z = x - 0.5; + z -= 0.5; + y = 0.5 * x + 0.5; } x = z / y; z = x*x; @@ -158,13 +158,13 @@ long double log10l(long double x) /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */ if (x < SQRTH) { e -= 1; - x = ldexpl(x, 1) - 1.0L; /* 2x - 1 */ + x = 2.0*x - 1.0; } else { - x = x - 1.0L; + x = x - 1.0; } z = x*x; y = x * (z * __polevll(x, P, 6) / __p1evll(x, Q, 7)); - y = y - ldexpl(z, -1); /* -0.5x^2 + ... */ + y = y - 0.5*z; done: /* Multiply log of fraction by log10(e) diff --git a/src/math/log1p.c b/src/math/log1p.c index f7154d0c..6c67249c 100644 --- a/src/math/log1p.c +++ b/src/math/log1p.c @@ -88,8 +88,6 @@ Lp5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */ Lp6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */ Lp7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */ -static const double zero = 0.0; - double log1p(double x) { double hfsq,f,c,s,z,R,u; @@ -102,12 +100,12 @@ double log1p(double x) if (hx < 0x3FDA827A) { /* 1+x < sqrt(2)+ */ if (ax >= 0x3ff00000) { /* 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 */ } if (ax < 0x3e200000) { /* |x| < 2**-29 */ /* 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 - x*x*0.5; } @@ -151,9 +149,9 @@ double log1p(double x) } hfsq = 0.5*f*f; if (hu == 0) { /* |f| < 2**-20 */ - if (f == zero) { + if (f == 0.0) { if(k == 0) - return zero; + return 0.0; c += k*ln2_lo; return k*ln2_hi + c; } diff --git a/src/math/log1pf.c b/src/math/log1pf.c index 75eeb371..39832d28 100644 --- a/src/math/log1pf.c +++ b/src/math/log1pf.c @@ -27,8 +27,6 @@ Lp5 = 1.8183572590e-01, /* 3E3A3325 */ Lp6 = 1.5313838422e-01, /* 3E1CD04F */ Lp7 = 1.4798198640e-01; /* 3E178897 */ -static const float zero = 0.0; - float log1pf(float x) { float hfsq,f,c,s,z,R,u; @@ -41,12 +39,12 @@ float log1pf(float x) if (hx < 0x3ed413d0) { /* 1+x < sqrt(2)+ */ if (ax >= 0x3f800000) { /* x <= -1.0 */ 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 */ } if (ax < 0x38000000) { /* |x| < 2**-15 */ /* 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 - x*x*0.5f; } @@ -91,9 +89,9 @@ float log1pf(float x) } hfsq = 0.5f * f * f; if (hu == 0) { /* |f| < 2**-20 */ - if (f == zero) { + if (f == 0.0f) { if (k == 0) - return zero; + return 0.0f; c += k*ln2_lo; return k*ln2_hi+c; } diff --git a/src/math/log1pl.c b/src/math/log1pl.c index 17eb4cef..1400d365 100644 --- a/src/math/log1pl.c +++ b/src/math/log1pl.c @@ -118,11 +118,11 @@ long double log1pl(long double xm1) if (xm1 == 0.0) return xm1; - x = xm1 + 1.0L; + x = xm1 + 1.0; /* Test for domain errors. */ - if (x <= 0.0L) { - if (x == 0.0L) + if (x <= 0.0) { + if (x == 0.0) return -INFINITY; return NAN; } @@ -136,12 +136,12 @@ long double log1pl(long double xm1) if (e > 2 || e < -2) { if (x < SQRTH) { /* 2(2x-1)/(2x+1) */ e -= 1; - z = x - 0.5L; - y = 0.5L * z + 0.5L; + z = x - 0.5; + y = 0.5 * z + 0.5; } else { /* 2 (x-1)/(x+1) */ - z = x - 0.5L; - z -= 0.5L; - y = 0.5L * x + 0.5L; + z = x - 0.5; + z -= 0.5; + y = 0.5 * x + 0.5; } x = z / y; z = x*x; @@ -156,12 +156,12 @@ long double log1pl(long double xm1) if (x < SQRTH) { e -= 1; if (e != 0) - x = 2.0 * x - 1.0L; + x = 2.0 * x - 1.0; else x = xm1; } else { if (e != 0) - x = x - 1.0L; + x = x - 1.0; else x = xm1; } diff --git a/src/math/log2.c b/src/math/log2.c index a5b8abdd..1974215d 100644 --- a/src/math/log2.c +++ b/src/math/log2.c @@ -27,8 +27,6 @@ two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */ ivln2hi = 1.44269504072144627571e+00, /* 0x3ff71547, 0x65200000 */ ivln2lo = 1.67517131648865118353e-10; /* 0x3de705fc, 0x2eefa200 */ -static const double zero = 0.0; - double log2(double x) { double f,hfsq,hi,lo,r,val_hi,val_lo,w,y; @@ -40,9 +38,9 @@ double log2(double x) k = 0; if (hx < 0x00100000) { /* x < 2**-1022 */ if (((hx&0x7fffffff)|lx) == 0) - return -two54/zero; /* log(+-0)=-inf */ + return -two54/0.0; /* log(+-0)=-inf */ if (hx < 0) - return (x-x)/zero; /* log(-#) = NaN */ + return (x-x)/0.0; /* log(-#) = NaN */ /* subnormal number, scale up x */ k -= 54; x *= two54; @@ -51,7 +49,7 @@ double log2(double x) if (hx >= 0x7ff00000) return x+x; if (hx == 0x3ff00000 && lx == 0) - return zero; /* log(1) = +0 */ + return 0.0; /* log(1) = +0 */ k += (hx>>20) - 1023; hx &= 0x000fffff; i = (hx+0x95f64) & 0x100000; diff --git a/src/math/log2f.c b/src/math/log2f.c index c9b9282c..e0d6a9e4 100644 --- a/src/math/log2f.c +++ b/src/math/log2f.c @@ -21,8 +21,6 @@ two25 = 3.3554432000e+07, /* 0x4c000000 */ ivln2hi = 1.4428710938e+00, /* 0x3fb8b000 */ ivln2lo = -1.7605285393e-04; /* 0xb9389ad4 */ -static const float zero = 0.0; - float log2f(float x) { float f,hfsq,hi,lo,r,y; @@ -33,9 +31,9 @@ float log2f(float x) k = 0; if (hx < 0x00800000) { /* x < 2**-126 */ if ((hx&0x7fffffff) == 0) - return -two25/zero; /* log(+-0)=-inf */ + return -two25/0.0f; /* log(+-0)=-inf */ if (hx < 0) - return (x-x)/zero; /* log(-#) = NaN */ + return (x-x)/0.0f; /* log(-#) = NaN */ /* subnormal number, scale up x */ k -= 25; x *= two25; @@ -44,7 +42,7 @@ float log2f(float x) if (hx >= 0x7f800000) return x+x; if (hx == 0x3f800000) - return zero; /* log(1) = +0 */ + return 0.0f; /* log(1) = +0 */ k += (hx>>23) - 127; hx &= 0x007fffff; i = (hx+(0x4afb0d))&0x800000; diff --git a/src/math/log2l.c b/src/math/log2l.c index 4339c033..8ebce9c4 100644 --- a/src/math/log2l.c +++ b/src/math/log2l.c @@ -121,8 +121,8 @@ long double log2l(long double x) return x; if (x == INFINITY) return x; - if (x <= 0.0L) { - if (x == 0.0L) + if (x <= 0.0) { + if (x == 0.0) return -INFINITY; return NAN; } @@ -139,12 +139,12 @@ long double log2l(long double x) if (e > 2 || e < -2) { if (x < SQRTH) { /* 2(2x-1)/(2x+1) */ e -= 1; - z = x - 0.5L; - y = 0.5L * z + 0.5L; + z = x - 0.5; + y = 0.5 * z + 0.5; } else { /* 2 (x-1)/(x+1) */ - z = x - 0.5L; - z -= 0.5L; - y = 0.5L * x + 0.5L; + z = x - 0.5; + z -= 0.5; + y = 0.5 * x + 0.5; } x = z / y; z = x*x; @@ -155,13 +155,13 @@ long double log2l(long double x) /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */ if (x < SQRTH) { e -= 1; - x = ldexpl(x, 1) - 1.0L; /* 2x - 1 */ + x = 2.0*x - 1.0; } else { - x = x - 1.0L; + x = x - 1.0; } z = x*x; y = x * (z * __polevll(x, P, 6) / __p1evll(x, Q, 7)); - y = y - ldexpl(z, -1); /* -0.5x^2 + ... */ + y = y - 0.5*z; done: /* Multiply log of fraction by log2(e) diff --git a/src/math/logf.c b/src/math/logf.c index a4ed697b..c7f7dbe6 100644 --- a/src/math/logf.c +++ b/src/math/logf.c @@ -25,8 +25,6 @@ Lg2 = 0xccce13.0p-25, /* 0.40000972152 */ Lg3 = 0x91e9ee.0p-25, /* 0.28498786688 */ Lg4 = 0xf89e26.0p-26; /* 0.24279078841 */ -static const float zero = 0.0; - float logf(float x) { float hfsq,f,s,z,R,w,t1,t2,dk; @@ -37,9 +35,9 @@ float logf(float x) k = 0; if (ix < 0x00800000) { /* x < 2**-126 */ if ((ix & 0x7fffffff) == 0) - return -two25/zero; /* log(+-0)=-inf */ + return -two25/0.0f; /* log(+-0)=-inf */ if (ix < 0) - return (x-x)/zero; /* log(-#) = NaN */ + return (x-x)/0.0f; /* log(-#) = NaN */ /* subnormal number, scale up x */ k -= 25; x *= two25; @@ -54,9 +52,9 @@ float logf(float x) k += i>>23; f = x - 1.0f; if ((0x007fffff & (0x8000 + ix)) < 0xc000) { /* -2**-9 <= f < 2**-9 */ - if (f == zero) { + if (f == 0.0f) { if (k == 0) - return zero; + return 0.0f; dk = (float)k; return dk*ln2_hi + dk*ln2_lo; } diff --git a/src/math/logl.c b/src/math/logl.c index ee7ca64a..ffd81038 100644 --- a/src/math/logl.c +++ b/src/math/logl.c @@ -119,8 +119,8 @@ long double logl(long double x) return x; if (x == INFINITY) return x; - if (x <= 0.0L) { - if (x == 0.0L) + if (x <= 0.0) { + if (x == 0.0) return -INFINITY; return NAN; } @@ -137,12 +137,12 @@ long double logl(long double x) if (e > 2 || e < -2) { if (x < SQRTH) { /* 2(2x-1)/(2x+1) */ e -= 1; - z = x - 0.5L; - y = 0.5L * z + 0.5L; + z = x - 0.5; + y = 0.5 * z + 0.5; } else { /* 2 (x-1)/(x+1) */ - z = x - 0.5L; - z -= 0.5L; - y = 0.5L * x + 0.5L; + z = x - 0.5; + z -= 0.5; + y = 0.5 * x + 0.5; } x = z / y; z = x*x; @@ -156,14 +156,14 @@ long double logl(long double x) /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */ if (x < SQRTH) { e -= 1; - x = ldexpl(x, 1) - 1.0L; /* 2x - 1 */ + x = 2.0*x - 1.0; } else { - x = x - 1.0L; + x = x - 1.0; } z = x*x; y = x * (z * __polevll(x, P, 6) / __p1evll(x, Q, 6)); y = y + e * C2; - z = y - ldexpl(z, -1); /* y - 0.5 * z */ + z = y - 0.5*z; /* Note, the sum of above terms does not exceed x/4, * so it contributes at most about 1/4 lsb to the error. */ diff --git a/src/math/modf.c b/src/math/modf.c index ff85b2a3..cca3b652 100644 --- a/src/math/modf.c +++ b/src/math/modf.c @@ -21,8 +21,6 @@ #include "libm.h" -static const double one = 1.0; - double modf(double x, double *iptr) { int32_t i0,i1,j0; @@ -51,7 +49,7 @@ double modf(double x, double *iptr) *iptr = x; return 0.0 / x; } - *iptr = x*one; + *iptr = x; GET_HIGH_WORD(high, x); INSERT_WORDS(x, high & 0x80000000, 0); /* return +-0 */ return x; diff --git a/src/math/modff.c b/src/math/modff.c index d535314c..bf6e4ced 100644 --- a/src/math/modff.c +++ b/src/math/modff.c @@ -1,51 +1,33 @@ -/* origin: FreeBSD /usr/src/lib/msun/src/s_modff.c */ -/* - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - #include "libm.h" -static const float one = 1.0; - float modff(float x, float *iptr) { - int32_t i0,j0; - uint32_t i; + uint32_t u, mask; + int e; - GET_FLOAT_WORD(i0, x); - j0 = ((i0>>23) & 0xff) - 0x7f; /* exponent of x */ - if (j0 < 23) { /* integer part in x */ - if (j0 < 0) { /* |x| < 1 */ - SET_FLOAT_WORD(*iptr, i0 & 0x80000000); /* *iptr = +-0 */ - return x; - } - i = 0x007fffff >> j0; - if ((i0&i) == 0) { /* x is integral */ - uint32_t ix; - *iptr = x; - GET_FLOAT_WORD(ix, x); - SET_FLOAT_WORD(x, ix & 0x80000000); /* return +-0 */ - return x; - } - SET_FLOAT_WORD(*iptr, i0&~i); - return x - *iptr; - } else { /* no fraction part */ - uint32_t ix; - *iptr = x*one; - if (x != x) /* NaN */ + GET_FLOAT_WORD(u, x); + e = (int)(u>>23 & 0xff) - 0x7f; + + /* no fractional part */ + if (e >= 23) { + *iptr = x; + if (e == 0x80 && u<<9 != 0) /* nan */ return x; - GET_FLOAT_WORD(ix, x); - SET_FLOAT_WORD(x, ix & 0x80000000); /* return +-0 */ + SET_FLOAT_WORD(x, u & 0x80000000); + return x; + } + /* no integral part */ + if (e < 0) { + SET_FLOAT_WORD(*iptr, u & 0x80000000); + return x; + } + + mask = 0x007fffff>>e; + if ((u & mask) == 0) { + *iptr = x; + SET_FLOAT_WORD(x, u & 0x80000000); return x; } + SET_FLOAT_WORD(*iptr, u & ~mask); + return x - *iptr; } diff --git a/src/math/modfl.c b/src/math/modfl.c index 2ca67b11..6520a1c2 100644 --- a/src/math/modfl.c +++ b/src/math/modfl.c @@ -54,7 +54,7 @@ long double modfl(long double x, long double *iptr) /* The number of fraction bits in manh, not counting the integer bit */ #define HIBITS (LDBL_MANT_DIG - LDBL_MANL_SIZE) -static const long double zero[] = { 0.0L, -0.0L }; +static const long double zero[] = { 0.0, -0.0 }; long double modfl(long double x, long double *iptr) { @@ -81,7 +81,7 @@ long double modfl(long double x, long double *iptr) return x - ux.e; } else if (e >= LDBL_MANT_DIG - 1) { /* x has no fraction part. */ *iptr = x; - if (x != x) /* Handle NaNs. */ + if (e == LDBL_MAX_EXP && (ux.bits.manh|ux.bits.manl)) /* nan */ return x; return zero[ux.bits.sign]; } else { /* Fraction part is in manl. */ diff --git a/src/math/pow.c b/src/math/pow.c index 5bcedd5b..052825a6 100644 --- a/src/math/pow.c +++ b/src/math/pow.c @@ -59,9 +59,6 @@ static const double bp[] = {1.0, 1.5,}, dp_h[] = { 0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */ dp_l[] = { 0.0, 1.35003920212974897128e-08,}, /* 0x3E4CFDEB, 0x43CFD006 */ -zero = 0.0, -one = 1.0, -two = 2.0, two53 = 9007199254740992.0, /* 0x43400000, 0x00000000 */ huge = 1.0e300, tiny = 1.0e-300, @@ -101,15 +98,15 @@ double pow(double x, double y) ix = hx & 0x7fffffff; iy = hy & 0x7fffffff; - /* y == zero: x**0 = 1 */ + /* y == 0.0: x**0 = 1 */ if ((iy|ly) == 0) - return one; + return 1.0; /* x == 1: 1**y = 1, even if y is NaN */ 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) || iy > 0x7ff00000 || (iy == 0x7ff00000 && ly != 0)) return (x+0.0) + (y+0.0); @@ -141,15 +138,15 @@ double pow(double x, double y) if (ly == 0) { if (iy == 0x7ff00000) { /* y is +-inf */ if (((ix-0x3ff00000)|lx) == 0) /* (-1)**+-inf is 1 */ - return one; + return 1.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 */ - return hy < 0 ? -y : zero; + return hy < 0 ? -y : 0.0; } if (iy == 0x3ff00000) { /* y is +-1 */ if (hy < 0) - return one/x; + return 1.0/x; return x; } 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 */ z = ax; if (hy < 0) /* z = (1/|x|) */ - z = one/z; + z = 1.0/z; if (hx < 0) { if (((ix-0x3ff00000)|yisint) == 0) { z = (z-z)/(z-z); /* (-1)**non-int is NaN */ @@ -187,9 +184,9 @@ double pow(double x, double y) if ((n|yisint) == 0) 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) - s = -one;/* (-ve)**(odd int) */ + s = -1.0;/* (-ve)**(odd int) */ /* |y| is huge */ 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; /* now |1-x| is tiny <= 2**-20, suffice to compute 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)); u = ivln2_h*t; /* ivln2_h has 21 sig. bits */ 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) */ 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; s_h = ss; SET_LOW_WORD(s_h, 0); /* t_h=ax+bp[k] High */ - t_h = zero; + t_h = 0.0; SET_HIGH_WORD(t_h, ((ix>>1)|0x20000000) + 0x00080000 + (k<<18)); t_l = ax - (t_h-bp[k]); 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] */ n = j + (0x00100000>>(k+1)); k = ((n&0x7fffffff)>>20) - 0x3ff; /* new k for n */ - t = zero; + t = 0.0; SET_HIGH_WORD(t, n & ~(0x000fffff>>k)); n = ((n&0x000fffff)|0x00100000)>>(20-k); if (j < 0) @@ -314,8 +311,8 @@ double pow(double x, double y) w = v - (z-u); t = z*z; t1 = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5)))); - r = (z*t1)/(t1-two) - (w + z*w); - z = one - (r-z); + r = (z*t1)/(t1-2.0) - (w + z*w); + z = 1.0 - (r-z); GET_HIGH_WORD(j, z); j += n<<20; if ((j>>20) <= 0) /* subnormal output */ diff --git a/src/math/powf.c b/src/math/powf.c index 01aced0e..c101d95c 100644 --- a/src/math/powf.c +++ b/src/math/powf.c @@ -19,9 +19,6 @@ static const float bp[] = {1.0, 1.5,}, dp_h[] = { 0.0, 5.84960938e-01,}, /* 0x3f15c000 */ dp_l[] = { 0.0, 1.56322085e-06,}, /* 0x35d1cfdc */ -zero = 0.0, -one = 1.0, -two = 2.0, two24 = 16777216.0, /* 0x4b800000 */ huge = 1.0e30, tiny = 1.0e-30, @@ -60,15 +57,15 @@ float powf(float x, float y) ix = hx & 0x7fffffff; iy = hy & 0x7fffffff; - /* y == zero: x**0 = 1 */ + /* y == 0: x**0 = 1 */ if (iy == 0) - return one; + return 1.0f; /* x == 1: 1**y = 1, even if y is NaN */ 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) return (x+0.0f) + (y+0.0f); @@ -92,15 +89,15 @@ float powf(float x, float y) /* special value of y */ if (iy == 0x7f800000) { /* y is +-inf */ if (ix == 0x3f800000) /* (-1)**+-inf is 1 */ - return one; + return 1.0f; 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 */ - return hy < 0 ? -y : zero; + return hy < 0 ? -y : 0.0f; } if (iy == 0x3f800000) { /* y is +-1 */ if (hy < 0) - return one/x; + return 1.0f/x; return x; } 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 */ z = ax; if (hy < 0) /* z = (1/|x|) */ - z = one/z; + z = 1.0f/z; if (hx < 0) { if (((ix-0x3f800000)|yisint) == 0) { z = (z-z)/(z-z); /* (-1)**non-int is NaN */ @@ -131,9 +128,9 @@ float powf(float x, float y) if ((n|yisint) == 0) 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) */ - sn = -one; + sn = -1.0f; /* |y| is huge */ 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) */ 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_h = s; GET_FLOAT_WORD(is, s_h); @@ -257,8 +254,8 @@ float powf(float x, float y) w = v - (z - u); t = z*z; t1 = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5)))); - r = (z*t1)/(t1-two) - (w+z*w); - z = one - (r - z); + r = (z*t1)/(t1-2.0f) - (w+z*w); + z = 1.0f - (r - z); GET_FLOAT_WORD(j, z); j += n<<23; if ((j>>23) <= 0) /* subnormal output */ diff --git a/src/math/powl.c b/src/math/powl.c index a6ee275f..a1d2f076 100644 --- a/src/math/powl.c +++ b/src/math/powl.c @@ -203,44 +203,44 @@ long double powl(long double x, long double y) volatile long double z=0; long double w=0, W=0, Wa=0, Wb=0, ya=0, yb=0, u=0; - if (y == 0.0L) - return 1.0L; + if (y == 0.0) + return 1.0; if (isnan(x)) return x; if (isnan(y)) return y; - if (y == 1.0L) + if (y == 1.0) return x; // FIXME: this is wrong, see pow special cases in c99 F.9.4.4 - if (!isfinite(y) && (x == -1.0L || x == 1.0L) ) + if (!isfinite(y) && (x == -1.0 || x == 1.0) ) return y - y; /* +-1**inf is NaN */ - if (x == 1.0L) - return 1.0L; + if (x == 1.0) + return 1.0; if (y >= LDBL_MAX) { - if (x > 1.0L) + if (x > 1.0) return INFINITY; - if (x > 0.0L && x < 1.0L) - return 0.0L; - if (x < -1.0L) + if (x > 0.0 && x < 1.0) + return 0.0; + if (x < -1.0) return INFINITY; - if (x > -1.0L && x < 0.0L) - return 0.0L; + if (x > -1.0 && x < 0.0) + return 0.0; } if (y <= -LDBL_MAX) { - if (x > 1.0L) - return 0.0L; - if (x > 0.0L && x < 1.0L) + if (x > 1.0) + return 0.0; + if (x > 0.0 && x < 1.0) return INFINITY; - if (x < -1.0L) - return 0.0L; - if (x > -1.0L && x < 0.0L) + if (x < -1.0) + return 0.0; + if (x > -1.0 && x < 0.0) return INFINITY; } if (x >= LDBL_MAX) { - if (y > 0.0L) + if (y > 0.0) return INFINITY; - return 0.0L; + return 0.0; } w = floorl(y); @@ -253,29 +253,29 @@ long double powl(long double x, long double y) yoddint = 0; if (iyflg) { ya = fabsl(y); - ya = floorl(0.5L * ya); - yb = 0.5L * fabsl(w); + ya = floorl(0.5 * ya); + yb = 0.5 * fabsl(w); if( ya != yb ) yoddint = 1; } if (x <= -LDBL_MAX) { - if (y > 0.0L) { + if (y > 0.0) { if (yoddint) return -INFINITY; return INFINITY; } - if (y < 0.0L) { + if (y < 0.0) { if (yoddint) - return -0.0L; + return -0.0; return 0.0; } } nflg = 0; /* flag = 1 if x<0 raised to integer power */ - if (x <= 0.0L) { - if (x == 0.0L) { + if (x <= 0.0) { + if (x == 0.0) { if (y < 0.0) { if (signbit(x) && yoddint) return -INFINITY; @@ -283,12 +283,12 @@ long double powl(long double x, long double y) } if (y > 0.0) { if (signbit(x) && yoddint) - return -0.0L; + return -0.0; return 0.0; } - if (y == 0.0L) - return 1.0L; /* 0**0 */ - return 0.0L; /* 0**y */ + if (y == 0.0) + return 1.0; /* 0**0 */ + return 0.0; /* 0**y */ } if (iyflg == 0) return (x - x) / (x - x); /* (x<0)**(non-int) is NaN */ @@ -343,7 +343,7 @@ long double powl(long double x, long double y) */ z = x*x; w = x * (z * __polevll(x, P, 3) / __p1evll(x, Q, 3)); - w = w - ldexpl(z, -1); /* w - 0.5 * z */ + w = w - 0.5*z; /* Convert to base 2 logarithm: * multiply by log2(e) = 1 + LOG2EA @@ -355,7 +355,8 @@ long double powl(long double x, long double y) /* Compute exponent term of the base 2 logarithm. */ w = -i; - w = ldexpl(w, -LNXT); /* divide by NXT */ + // TODO: use w * 0x1p-5; + w = scalbnl(w, -LNXT); /* divide by NXT */ w += e; /* Now base 2 log of x is w + z. */ @@ -380,7 +381,7 @@ long double powl(long double x, long double y) H = Fb + Gb; Ha = reducl(H); - w = ldexpl( Ga+Ha, LNXT ); + w = scalbnl( Ga+Ha, LNXT ); /* Test the power of 2 for overflow */ if (w > MEXP) @@ -391,9 +392,9 @@ long double powl(long double x, long double y) e = w; Hb = H - Ha; - if (Hb > 0.0L) { + if (Hb > 0.0) { e += 1; - Hb -= 1.0L/NXT; /*0.0625L;*/ + Hb -= 1.0/NXT; /*0.0625L;*/ } /* Now the product y * log2(x) = Hb + e/NXT. @@ -415,16 +416,16 @@ long double powl(long double x, long double y) w = douba(e); z = w * z; /* 2**-e * ( 1 + (2**Hb-1) ) */ z = z + w; - z = ldexpl(z, i); /* multiply by integer power of 2 */ + z = scalbnl(z, i); /* multiply by integer power of 2 */ if (nflg) { /* For negative x, * find out if the integer exponent * is odd or even. */ - w = ldexpl(y, -1); + w = 0.5*y; w = floorl(w); - w = ldexpl(w, 1); + w = 2.0*w; if (w != y) z = -z; /* odd exponent */ } @@ -438,9 +439,9 @@ static long double reducl(long double x) { long double t; - t = ldexpl(x, LNXT); + t = scalbnl(x, LNXT); t = floorl(t); - t = ldexpl(t, -LNXT); + t = scalbnl(t, -LNXT); return t; } @@ -483,18 +484,18 @@ static long double powil(long double x, int nn) long double s; int n, e, sign, asign, lx; - if (x == 0.0L) { + if (x == 0.0) { if (nn == 0) - return 1.0L; + return 1.0; else if (nn < 0) return LDBL_MAX; - return 0.0L; + return 0.0; } if (nn == 0) - return 1.0L; + return 1.0; - if (x < 0.0L) { + if (x < 0.0) { asign = -1; x = -x; } else @@ -516,7 +517,7 @@ static long double powil(long double x, int nn) e = (lx - 1)*n; if ((e == 0) || (e > 64) || (e < -64)) { s = (s - 7.0710678118654752e-1L) / (s + 7.0710678118654752e-1L); - s = (2.9142135623730950L * s - 0.5L + lx) * nn * LOGE2L; + s = (2.9142135623730950L * s - 0.5 + lx) * nn * LOGE2L; } else { s = LOGE2L * e; } @@ -530,8 +531,8 @@ static long double powil(long double x, int nn) * since roundoff error in 1.0/x will be amplified. * The precise demarcation should be the gradual underflow threshold. */ - if (s < -MAXLOGL+2.0L) { - x = 1.0L/x; + if (s < -MAXLOGL+2.0) { + x = 1.0/x; sign = -sign; } @@ -539,7 +540,7 @@ static long double powil(long double x, int nn) if (n & 1) y = x; else { - y = 1.0L; + y = 1.0; asign = 0; } @@ -555,7 +556,7 @@ static long double powil(long double x, int nn) if (asign) y = -y; /* odd power of negative number */ if (sign < 0) - y = 1.0L/y; + y = 1.0/y; return y; } diff --git a/src/math/remainder.c b/src/math/remainder.c index db176c88..2dede3a8 100644 --- a/src/math/remainder.c +++ b/src/math/remainder.c @@ -20,8 +20,6 @@ #include "libm.h" -static const double zero = 0.0; - double remainder(double x, double p) { int32_t hx,hp; @@ -35,17 +33,15 @@ double remainder(double x, double p) hx &= 0x7fffffff; /* purge off exception values */ - if ((hp|lp) == 0) /* p = 0 */ - return (x*p)/(x*p); - if (hx >= 0x7ff00000 || /* x not finite */ + if ((hp|lp) == 0 || /* p = 0 */ + hx >= 0x7ff00000 || /* x not finite */ (hp >= 0x7ff00000 && (hp-0x7ff00000 | lp) != 0)) /* p is NaN */ - // FIXME: why long double? - return ((long double)x*p)/((long double)x*p); + return (x*p)/(x*p); if (hp <= 0x7fdfffff) x = fmod(x, p+p); /* now x < 2p */ if (((hx-hp)|(lx-lp)) == 0) - return zero*x; + return 0.0*x; x = fabs(x); p = fabs(p); if (hp < 0x00200000) { diff --git a/src/math/remainderf.c b/src/math/remainderf.c index 61c3c660..77671039 100644 --- a/src/math/remainderf.c +++ b/src/math/remainderf.c @@ -15,8 +15,6 @@ #include "libm.h" -static const float zero = 0.0; - float remainderf(float x, float p) { int32_t hx,hp; @@ -30,16 +28,13 @@ float remainderf(float x, float p) hx &= 0x7fffffff; /* purge off exception values */ - if (hp == 0) /* p = 0 */ + if (hp == 0 || hx >= 0x7f800000 || hp > 0x7f800000) /* p = 0, x not finite, p is NaN */ return (x*p)/(x*p); - if (hx >= 0x7f800000 || hp > 0x7f800000) /* x not finite, p is NaN */ - // FIXME: why long double? - return ((long double)x*p)/((long double)x*p); if (hp <= 0x7effffff) x = fmodf(x, p + p); /* now x < 2p */ if (hx - hp == 0) - return zero*x; + return 0.0f*x; x = fabsf(x); p = fabsf(p); if (hp < 0x01000000) { diff --git a/src/math/remquol.c b/src/math/remquol.c index dd18f35c..721231b4 100644 --- a/src/math/remquol.c +++ b/src/math/remquol.c @@ -48,7 +48,7 @@ typedef uint32_t manh_t; #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 diff --git a/src/math/rintl.c b/src/math/rintl.c index 1cc35df5..6a311a9d 100644 --- a/src/math/rintl.c +++ b/src/math/rintl.c @@ -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 * 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 x; diff --git a/src/math/scalbnl.c b/src/math/scalbnl.c index a5d0adba..c605b8da 100644 --- a/src/math/scalbnl.c +++ b/src/math/scalbnl.c @@ -29,7 +29,7 @@ long double scalbnl(long double x, int n) return x * 0x1p-16382L; } } - scale.e = 1.0L; + scale.e = 1.0; scale.bits.exp = 0x3fff + n; return x * scale.e; } diff --git a/src/math/sincosl.c b/src/math/sincosl.c index 378dc979..e14129a2 100644 --- a/src/math/sincosl.c +++ b/src/math/sincosl.c @@ -9,8 +9,6 @@ void sincosl(long double x, long double *sin, long double *cos) *cos = c; } #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 -#include "__rem_pio2l.h" - void sincosl(long double x, long double *sin, long double *cos) { union IEEEl2bits u; diff --git a/src/math/sinh.c b/src/math/sinh.c index 935879c5..0c67ba39 100644 --- a/src/math/sinh.c +++ b/src/math/sinh.c @@ -29,7 +29,7 @@ #include "libm.h" -static const double one = 1.0, huge = 1.0e307; +static const double huge = 1.0e307; double sinh(double x) { @@ -50,12 +50,12 @@ double sinh(double x) if (ix < 0x40360000) { /* |x|<22 */ if (ix < 0x3e300000) /* |x|<2**-28 */ /* raise inexact, return x */ - if (huge+x > one) + if (huge+x > 1.0) return x; t = expm1(fabs(x)); if (ix < 0x3ff00000) - return h*(2.0*t - t*t/(t+one)); - return h*(t + t/(t+one)); + return h*(2.0*t - t*t/(t+1.0)); + return h*(t + t/(t+1.0)); } /* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */ diff --git a/src/math/sinhf.c b/src/math/sinhf.c index fd11b849..b8d88224 100644 --- a/src/math/sinhf.c +++ b/src/math/sinhf.c @@ -15,7 +15,7 @@ #include "libm.h" -static const float one = 1.0, huge = 1.0e37; +static const float huge = 1.0e37; float sinhf(float x) { @@ -36,12 +36,12 @@ float sinhf(float x) if (ix < 0x41100000) { /* |x|<9 */ if (ix < 0x39800000) /* |x|<2**-12 */ /* raise inexact, return x */ - if (huge+x > one) + if (huge+x > 1.0f) return x; t = expm1f(fabsf(x)); if (ix < 0x3f800000) - return h*(2.0f*t - t*t/(t+one)); - return h*(t + t/(t+one)); + return h*(2.0f*t - t*t/(t+1.0f)); + return h*(t + t/(t+1.0f)); } /* |x| in [9, logf(maxfloat)] return 0.5*exp(|x|) */ diff --git a/src/math/sinhl.c b/src/math/sinhl.c index 2252dec9..8a54677e 100644 --- a/src/math/sinhl.c +++ b/src/math/sinhl.c @@ -35,7 +35,7 @@ long double sinhl(long double x) return sinh(x); } #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) { @@ -55,12 +55,12 @@ long double sinhl(long double x) /* |x| in [0,25], return sign(x)*0.5*(E+E/(E+1))) */ if (ix < 0x4003 || (ix == 0x4003 && i0 <= 0xc8000000)) { /* |x| < 25 */ if (ix < 0x3fdf) /* |x|<2**-32 */ - if (huge + x > one) + if (huge + x > 1.0) return x;/* sinh(tiny) = tiny with inexact */ t = expm1l(fabsl(x)); if (ix < 0x3fff) - return h*(2.0*t - t*t/(t + one)); - return h*(t + t/(t + one)); + return h*(2.0*t - t*t/(t + 1.0)); + return h*(t + t/(t + 1.0)); } /* |x| in [25, log(maxdouble)] return 0.5*exp(|x|) */ diff --git a/src/math/sinl.c b/src/math/sinl.c index 0b1aeb75..7e0b44f4 100644 --- a/src/math/sinl.c +++ b/src/math/sinl.c @@ -34,8 +34,6 @@ long double sinl(long double x) return sin(x); } #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 -#include "__rem_pio2l.h" - long double sinl(long double x) { union IEEEl2bits z; diff --git a/src/math/sqrt.c b/src/math/sqrt.c index 2ebd022b..b2775673 100644 --- a/src/math/sqrt.c +++ b/src/math/sqrt.c @@ -78,7 +78,7 @@ #include "libm.h" -static const double one = 1.0, tiny = 1.0e-300; +static const double tiny = 1.0e-300; double sqrt(double x) { @@ -161,13 +161,13 @@ double sqrt(double x) /* use floating add to find out rounding direction */ if ((ix0|ix1) != 0) { - z = one - tiny; /* raise inexact flag */ - if (z >= one) { - z = one + tiny; + z = 1.0 - tiny; /* raise inexact flag */ + if (z >= 1.0) { + z = 1.0 + tiny; if (q1 == (uint32_t)0xffffffff) { q1 = 0; q++; - } else if (z > one) { + } else if (z > 1.0) { if (q1 == (uint32_t)0xfffffffe) q++; q1 += 2; diff --git a/src/math/sqrtf.c b/src/math/sqrtf.c index 35c24e50..28cb4ad3 100644 --- a/src/math/sqrtf.c +++ b/src/math/sqrtf.c @@ -15,7 +15,7 @@ #include "libm.h" -static const float one = 1.0, tiny = 1.0e-30; +static const float tiny = 1.0e-30; float sqrtf(float x) { @@ -68,10 +68,10 @@ float sqrtf(float x) /* use floating add to find out rounding direction */ if (ix != 0) { - z = one - tiny; /* raise inexact flag */ - if (z >= one) { - z = one + tiny; - if (z > one) + z = 1.0f - tiny; /* raise inexact flag */ + if (z >= 1.0f) { + z = 1.0f + tiny; + if (z > 1.0f) q += 2; else q += q & 1; diff --git a/src/math/tanh.c b/src/math/tanh.c index 957c43e9..21138643 100644 --- a/src/math/tanh.c +++ b/src/math/tanh.c @@ -35,7 +35,7 @@ #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) { @@ -48,26 +48,26 @@ double tanh(double x) /* x is INF or NaN */ if (ix >= 0x7ff00000) { if (jx >= 0) - return one/x + one; /* tanh(+-inf)=+-1 */ + return 1.0f/x + 1.0f; /* tanh(+-inf)=+-1 */ else - return one/x - one; /* tanh(NaN) = NaN */ + return 1.0f/x - 1.0f; /* tanh(NaN) = NaN */ } if (ix < 0x40360000) { /* |x| < 22 */ if (ix < 0x3e300000) { /* |x| < 2**-28 */ /* tanh(tiny) = tiny with inexact */ - if (huge+x > one) + if (huge+x > 1.0f) return x; } if (ix >= 0x3ff00000) { /* |x| >= 1 */ - t = expm1(two*fabs(x)); - z = one - two/(t+two); + t = expm1(2.0f*fabs(x)); + z = 1.0f - 2.0f/(t+2.0f); } else { - t = expm1(-two*fabs(x)); - z= -t/(t+two); + t = expm1(-2.0f*fabs(x)); + z= -t/(t+2.0f); } } else { /* |x| >= 22, return +-1 */ - z = one - tiny; /* raise inexact */ + z = 1.0f - tiny; /* raise inexact */ } return jx >= 0 ? z : -z; } diff --git a/src/math/tanhf.c b/src/math/tanhf.c index 97d0eb53..7cb459d0 100644 --- a/src/math/tanhf.c +++ b/src/math/tanhf.c @@ -15,7 +15,9 @@ #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) { @@ -28,26 +30,26 @@ float tanhf(float x) /* x is INF or NaN */ if(ix >= 0x7f800000) { if (jx >= 0) - return one/x + one; /* tanh(+-inf)=+-1 */ + return 1.0f/x + 1.0f; /* tanh(+-inf)=+-1 */ else - return one/x - one; /* tanh(NaN) = NaN */ + return 1.0f/x - 1.0f; /* tanh(NaN) = NaN */ } if (ix < 0x41100000) { /* |x| < 9 */ if (ix < 0x39800000) { /* |x| < 2**-12 */ /* tanh(tiny) = tiny with inexact */ - if (huge+x > one) + if (huge+x > 1.0f) return x; } if (ix >= 0x3f800000) { /* |x|>=1 */ - t = expm1f(two*fabsf(x)); - z = one - two/(t+two); + t = expm1f(2.0f*fabsf(x)); + z = 1.0f - 2.0f/(t+2.0f); } else { - t = expm1f(-two*fabsf(x)); - z = -t/(t+two); + t = expm1f(-2.0f*fabsf(x)); + z = -t/(t+2.0f); } } else { /* |x| >= 9, return +-1 */ - z = one - tiny; /* raise inexact */ + z = 1.0f - tiny; /* raise inexact */ } return jx >= 0 ? z : -z; } diff --git a/src/math/tanhl.c b/src/math/tanhl.c index e62be59b..92efb20d 100644 --- a/src/math/tanhl.c +++ b/src/math/tanhl.c @@ -41,7 +41,7 @@ long double tanhl(long double x) return tanh(x); } #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) { @@ -57,8 +57,8 @@ long double tanhl(long double x) if (ix == 0x7fff) { /* for NaN it's not important which branch: tanhl(NaN) = NaN */ if (se & 0x8000) - return one/x-one; /* tanhl(-inf)= -1; */ - return one/x+one; /* tanhl(+inf)= +1 */ + return 1.0/x-1.0; /* tanhl(-inf)= -1; */ + return 1.0/x+1.0; /* tanhl(+inf)= +1 */ } /* |x| < 23 */ @@ -66,17 +66,17 @@ long double tanhl(long double x) if ((ix|jj0|jj1) == 0) /* x == +- 0 */ return x; 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 */ - t = expm1l(two*fabsl(x)); - z = one - two/(t+two); + t = expm1l(2.0*fabsl(x)); + z = 1.0 - 2.0/(t+2.0); } else { - t = expm1l(-two*fabsl(x)); - z = -t/(t+two); + t = expm1l(-2.0*fabsl(x)); + z = -t/(t+2.0); } /* |x| > 23, return +-1 */ } else { - z = one - tiny; /* raise inexact flag */ + z = 1.0 - tiny; /* raise inexact flag */ } return se & 0x8000 ? -z : z; } diff --git a/src/math/tanl.c b/src/math/tanl.c index 462ead91..0194eaf7 100644 --- a/src/math/tanl.c +++ b/src/math/tanl.c @@ -38,8 +38,6 @@ long double tanl(long double x) return tan(x); } #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 -#include "__rem_pio2l.h" - long double tanl(long double x) { union IEEEl2bits z; diff --git a/src/math/tgammal.c b/src/math/tgammal.c index 1b8fddea..1b460da5 100644 --- a/src/math/tgammal.c +++ b/src/math/tgammal.c @@ -183,18 +183,18 @@ static long double stirf(long double x) { long double y, w, v; - w = 1.0L/x; + w = 1.0/x; /* For large x, use rational coefficients from the analytical expansion. */ - if (x > 1024.0L) + if (x > 1024.0) w = (((((6.97281375836585777429E-5L * w + 7.84039221720066627474E-4L) * w - 2.29472093621399176955E-4L) * w - 2.68132716049382716049E-3L) * w + 3.47222222222222222222E-3L) * w + 8.33333333333333333333E-2L) * w - + 1.0L; + + 1.0; else - w = 1.0L + w * __polevll(w, STIR, 8); + w = 1.0 + w * __polevll(w, STIR, 8); y = expl(x); if (x > MAXSTIR) { /* Avoid overflow in pow() */ v = powl(x, 0.5L * x - 0.25L); @@ -219,10 +219,10 @@ long double tgammal(long double x) if (x == -INFINITY) return x - x; q = fabsl(x); - if (q > 13.0L) { + if (q > 13.0) { if (q > MAXGAML) goto goverf; - if (x < 0.0L) { + if (x < 0.0) { p = floorl(q); if (p == q) return (x - x) / (x - x); @@ -231,7 +231,7 @@ long double tgammal(long double x) signgam = -1; z = q - p; if (z > 0.5L) { - p += 1.0L; + p += 1.0; z = q - p; } z = q * sinl(PIL * z); @@ -247,25 +247,25 @@ goverf: return signgam * z; } - z = 1.0L; - while (x >= 3.0L) { - x -= 1.0L; + z = 1.0; + while (x >= 3.0) { + x -= 1.0; z *= x; } while (x < -0.03125L) { z /= x; - x += 1.0L; + x += 1.0; } if (x <= 0.03125L) goto small; - while (x < 2.0L) { + while (x < 2.0) { z /= x; - x += 1.0L; + x += 1.0; } - if (x == 2.0L) + if (x == 2.0) return z; - x -= 2.0L; + x -= 2.0; p = __polevll(x, P, 7); q = __polevll(x, Q, 8); z = z * p / q; @@ -274,9 +274,9 @@ goverf: return z; small: - if (x == 0.0L) + if (x == 0.0) return (x - x) / (x - x); - if (x < 0.0L) { + if (x < 0.0) { x = -x; q = z / (x * __polevll(x, SN, 8)); signgam = -1; |