diff options
| -rw-r--r-- | src/math/__signbit.c | 2 | ||||
| -rw-r--r-- | src/math/__signbitf.c | 2 | ||||
| -rw-r--r-- | src/math/cbrtl.c | 33 | ||||
| -rw-r--r-- | src/math/nearbyint.c | 23 | ||||
| -rw-r--r-- | src/math/nearbyintf.c | 14 | ||||
| -rw-r--r-- | src/math/nearbyintl.c | 11 | ||||
| -rw-r--r-- | src/math/pow.c | 65 | ||||
| -rw-r--r-- | src/math/powf.c | 35 | ||||
| -rw-r--r-- | src/math/powl.c | 139 | ||||
| -rw-r--r-- | src/math/remquo.c | 1 |
10 files changed, 127 insertions, 198 deletions
diff --git a/src/math/__signbit.c b/src/math/__signbit.c index ffe717ce..e700b6b7 100644 --- a/src/math/__signbit.c +++ b/src/math/__signbit.c @@ -1,6 +1,6 @@ #include "libm.h" -// FIXME: macro +// FIXME: macro in math.h int __signbit(double x) { union { diff --git a/src/math/__signbitf.c b/src/math/__signbitf.c index ff3e81ff..40ad3cfd 100644 --- a/src/math/__signbitf.c +++ b/src/math/__signbitf.c @@ -1,6 +1,6 @@ #include "libm.h" -// FIXME +// FIXME: macro in math.h int __signbitf(float x) { union { diff --git a/src/math/cbrtl.c b/src/math/cbrtl.c index 5297d68f..0af65ef5 100644 --- a/src/math/cbrtl.c +++ b/src/math/cbrtl.c @@ -23,9 +23,9 @@ long double cbrtl(long double x) return cbrt(x); } #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 -#define BIAS (LDBL_MAX_EXP - 1) -static const unsigned -B1 = 709958130; /* B1 = (127-127.0/3-0.03306235651)*2**23 */ + +#define BIAS (LDBL_MAX_EXP - 1) +static const unsigned B1 = 709958130; /* B1 = (127-127.0/3-0.03306235651)*2**23 */ long double cbrtl(long double x) { @@ -48,25 +48,10 @@ long double cbrtl(long double x) if (k == BIAS + LDBL_MAX_EXP) return x + x; -// FIXME: extended precision is default on linux.. -#undef __i386__ -#ifdef __i386__ - fp_prec_t oprec; - - oprec = fpgetprec(); - if (oprec != FP_PE) - fpsetprec(FP_PE); -#endif - if (k == 0) { /* If x = +-0, then cbrt(x) = +-0. */ - if ((u.bits.manh | u.bits.manl) == 0) { -#ifdef __i386__ - if (oprec != FP_PE) - fpsetprec(oprec); -#endif - return (x); - } + if ((u.bits.manh | u.bits.manl) == 0) + return x; /* Adjust subnormal numbers. */ u.e *= 0x1.0p514; k = u.bits.exp; @@ -118,7 +103,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). */ - t = dt + (0x1.0p32L + 0x1.0p-32L) - 0x1.0p32; + t = dt + (0x1.0p32L + 0x1.0p-31L) - 0x1.0p32; #elif LDBL_MANT_DIG == 113 /* * Round dt away from zero to 47 bits. Since we don't trust the 47, @@ -129,8 +114,6 @@ long double cbrtl(long double x) * dt is much more accurate than needed. */ t = dt + 0x2.0p-46 + 0x1.0p60L - 0x1.0p60; -#else -#error "Unsupported long double format" #endif /* @@ -144,10 +127,6 @@ long double cbrtl(long double x) t = t+t*r; /* error <= 0.5 + 0.5/3 + epsilon */ t *= v.e; -#ifdef __i386__ - if (oprec != FP_PE) - fpsetprec(oprec); -#endif return t; } #endif diff --git a/src/math/nearbyint.c b/src/math/nearbyint.c index 714c55ca..7a4c58cf 100644 --- a/src/math/nearbyint.c +++ b/src/math/nearbyint.c @@ -1,20 +1,19 @@ #include <fenv.h> #include <math.h> -/* -rint may raise inexact (and it should not alter the fenv otherwise) -nearbyint must not raise inexact +/* nearbyint is the same as rint, but it must not raise the inexact exception */ -(according to ieee754r section 7.9 both functions should raise invalid -when the input is signaling nan, but c99 does not define snan so saving -and restoring the entire fenv should be fine) -*/ +double nearbyint(double x) +{ +#ifdef FE_INEXACT + int e; -double nearbyint(double x) { - fenv_t e; - - fegetenv(&e); + e = fetestexcept(FE_INEXACT); +#endif x = rint(x); - fesetenv(&e); +#ifdef FE_INEXACT + if (!e) + feclearexcept(FE_INEXACT); +#endif return x; } diff --git a/src/math/nearbyintf.c b/src/math/nearbyintf.c index 07df8f54..39c3d73b 100644 --- a/src/math/nearbyintf.c +++ b/src/math/nearbyintf.c @@ -1,11 +1,17 @@ #include <fenv.h> #include <math.h> -float nearbyintf(float x) { - fenv_t e; +float nearbyintf(float x) +{ +#ifdef FE_INEXACT + int e; - fegetenv(&e); + e = fetestexcept(FE_INEXACT); +#endif x = rintf(x); - fesetenv(&e); +#ifdef FE_INEXACT + if (!e) + feclearexcept(FE_INEXACT); +#endif return x; } diff --git a/src/math/nearbyintl.c b/src/math/nearbyintl.c index 2906f383..0ff4b1f9 100644 --- a/src/math/nearbyintl.c +++ b/src/math/nearbyintl.c @@ -10,11 +10,16 @@ long double nearbyintl(long double x) #include <fenv.h> long double nearbyintl(long double x) { - fenv_t e; +#ifdef FE_INEXACT + int e; - fegetenv(&e); + e = fetestexcept(FE_INEXACT); +#endif x = rintl(x); - fesetenv(&e); +#ifdef FE_INEXACT + if (!e) + feclearexcept(FE_INEXACT); +#endif return x; } #endif diff --git a/src/math/pow.c b/src/math/pow.c index 052825a6..f257814e 100644 --- a/src/math/pow.c +++ b/src/math/pow.c @@ -21,24 +21,28 @@ * * Special cases: * 1. (anything) ** 0 is 1 - * 2. (anything) ** 1 is itself - * 3. (anything except 1) ** NAN is NAN, 1 ** NAN is 1 + * 2. 1 ** (anything) is 1 + * 3. (anything except 1) ** NAN is NAN * 4. NAN ** (anything except 0) is NAN * 5. +-(|x| > 1) ** +INF is +INF * 6. +-(|x| > 1) ** -INF is +0 * 7. +-(|x| < 1) ** +INF is +0 * 8. +-(|x| < 1) ** -INF is +INF - * 9. +-1 ** +-INF is 1 + * 9. -1 ** +-INF is 1 * 10. +0 ** (+anything except 0, NAN) is +0 * 11. -0 ** (+anything except 0, NAN, odd integer) is +0 - * 12. +0 ** (-anything except 0, NAN) is +INF - * 13. -0 ** (-anything except 0, NAN, odd integer) is +INF - * 14. -0 ** (odd integer) = -( +0 ** (odd integer) ) - * 15. +INF ** (+anything except 0,NAN) is +INF - * 16. +INF ** (-anything except 0,NAN) is +0 - * 17. -INF ** (anything) = -0 ** (-anything) - * 18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer) - * 19. (-anything except 0 and inf) ** (non-integer) is NAN + * 12. +0 ** (-anything except 0, NAN) is +INF, raise divbyzero + * 13. -0 ** (-anything except 0, NAN, odd integer) is +INF, raise divbyzero + * 14. -0 ** (+odd integer) is -0 + * 15. -0 ** (-odd integer) is -INF, raise divbyzero + * 16. +INF ** (+anything except 0,NAN) is +INF + * 17. +INF ** (-anything except 0,NAN) is +0 + * 18. -INF ** (+odd integer) is -INF + * 19. -INF ** (anything) = -0 ** (-anything), (anything except odd integer) + * 20. (anything) ** 1 is (anything) + * 21. (anything) ** -1 is 1/(anything) + * 22. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer) + * 23. (-anything except 0 and inf) ** (non-integer) is NAN * * Accuracy: * pow(x,y) returns x**y nearly rounded. In particular @@ -98,18 +102,16 @@ double pow(double x, double y) ix = hx & 0x7fffffff; iy = hy & 0x7fffffff; - /* y == 0.0: x**0 = 1 */ + /* x**0 = 1, even if x is NaN */ if ((iy|ly) == 0) return 1.0; - - /* x == 1: 1**y = 1, even if y is NaN */ + /* 1**y = 1, even if y is NaN */ if (hx == 0x3ff00000 && lx == 0) return 1.0; - - /* y != 0.0: result is NaN if either arg is NaN */ + /* 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); + return x + y; /* determine if y is an odd int when x < 0 * yisint = 0 ... y is not an integer @@ -142,13 +144,10 @@ double pow(double x, double y) else if (ix >= 0x3ff00000) /* (|x|>1)**+-inf = inf,0 */ return hy >= 0 ? y : 0.0; else /* (|x|<1)**+-inf = 0,inf */ - return hy < 0 ? -y : 0.0; - } - if (iy == 0x3ff00000) { /* y is +-1 */ - if (hy < 0) - return 1.0/x; - return x; + return hy >= 0 ? 0.0 : -y; } + if (iy == 0x3ff00000) /* y is +-1 */ + return hy >= 0 ? x : 1.0/x; if (hy == 0x40000000) /* y is 2 */ return x*x; if (hy == 0x3fe00000) { /* y is 0.5 */ @@ -174,19 +173,13 @@ double pow(double x, double y) } } - /* CYGNUS LOCAL + fdlibm-5.3 fix: This used to be - n = (hx>>31)+1; - but ANSI C says a right shift of a signed negative quantity is - implementation defined. */ - n = ((uint32_t)hx>>31) - 1; - - /* (x<0)**(non-int) is NaN */ - if ((n|yisint) == 0) - return (x-x)/(x-x); - - s = 1.0; /* s (sign of result -ve**odd) = -1 else = 1 */ - if ((n|(yisint-1)) == 0) - s = -1.0;/* (-ve)**(odd int) */ + s = 1.0; /* sign of result */ + if (hx < 0) { + if (yisint == 0) /* (x<0)**(non-int) is NaN */ + return (x-x)/(x-x); + if (yisint == 1) /* (x<0)**(odd int) */ + s = -1.0; + } /* |y| is huge */ if (iy > 0x41e00000) { /* if |y| > 2**31 */ diff --git a/src/math/powf.c b/src/math/powf.c index c101d95c..427c8965 100644 --- a/src/math/powf.c +++ b/src/math/powf.c @@ -57,17 +57,15 @@ float powf(float x, float y) ix = hx & 0x7fffffff; iy = hy & 0x7fffffff; - /* y == 0: x**0 = 1 */ + /* x**0 = 1, even if x is NaN */ if (iy == 0) return 1.0f; - - /* x == 1: 1**y = 1, even if y is NaN */ + /* 1**y = 1, even if y is NaN */ if (hx == 0x3f800000) return 1.0f; - - /* y != 0: result is NaN if either arg is NaN */ + /* NaN if either arg is NaN */ if (ix > 0x7f800000 || iy > 0x7f800000) - return (x+0.0f) + (y+0.0f); + return x + y; /* determine if y is an odd int when x < 0 * yisint = 0 ... y is not an integer @@ -93,13 +91,10 @@ float powf(float x, float y) else if (ix > 0x3f800000) /* (|x|>1)**+-inf = inf,0 */ return hy >= 0 ? y : 0.0f; else /* (|x|<1)**+-inf = 0,inf */ - return hy < 0 ? -y : 0.0f; - } - if (iy == 0x3f800000) { /* y is +-1 */ - if (hy < 0) - return 1.0f/x; - return x; + return hy >= 0 ? 0.0f: -y; } + if (iy == 0x3f800000) /* y is +-1 */ + return hy >= 0 ? x : 1.0f/x; if (hy == 0x40000000) /* y is 2 */ return x*x; if (hy == 0x3f000000) { /* y is 0.5 */ @@ -122,15 +117,13 @@ float powf(float x, float y) return z; } - n = ((uint32_t)hx>>31) - 1; - - /* (x<0)**(non-int) is NaN */ - if ((n|yisint) == 0) - return (x-x)/(x-x); - - sn = 1.0f; /* s (sign of result -ve**odd) = -1 else = 1 */ - if ((n|(yisint-1)) == 0) /* (-ve)**(odd int) */ - sn = -1.0f; + sn = 1.0f; /* sign of result */ + if (hx < 0) { + if (yisint == 0) /* (x<0)**(non-int) is NaN */ + return (x-x)/(x-x); + if (yisint == 1) /* (x<0)**(odd int) */ + sn = -1.0f; + } /* |y| is huge */ if (iy > 0x4d000000) { /* if |y| > 2**27 */ diff --git a/src/math/powl.c b/src/math/powl.c index a1d2f076..2ee0b10b 100644 --- a/src/math/powl.c +++ b/src/math/powl.c @@ -78,8 +78,6 @@ long double powl(long double x, long double y) /* Table size */ #define NXT 32 -/* log2(Table size) */ -#define LNXT 5 /* log(1+x) = x - .5x^2 + x^3 * P(z)/Q(z) * on the domain 2^(-1/32) - 1 <= x <= 2^(1/32) - 1 @@ -203,38 +201,35 @@ 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.0) - return 1.0; - if (isnan(x)) + /* make sure no invalid exception is raised by nan comparision */ + if (isnan(x)) { + if (!isnan(y) && y == 0.0) + return 1.0; return x; - if (isnan(y)) + } + if (isnan(y)) { + if (x == 1.0) + return 1.0; return y; + } + if (x == 1.0) + return 1.0; /* 1**y = 1, even if y is nan */ + if (x == -1.0 && !isfinite(y)) + return 1.0; /* -1**inf = 1 */ + if (y == 0.0) + return 1.0; /* x**0 = 1, even if x is nan */ 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.0 || x == 1.0) ) - return y - y; /* +-1**inf is NaN */ - if (x == 1.0) - return 1.0; if (y >= LDBL_MAX) { - if (x > 1.0) + if (x > 1.0 || x < -1.0) return INFINITY; - if (x > 0.0 && x < 1.0) - return 0.0; - if (x < -1.0) - return INFINITY; - if (x > -1.0 && x < 0.0) + if (x != 0.0) return 0.0; } if (y <= -LDBL_MAX) { - if (x > 1.0) + if (x > 1.0 || x < -1.0) return 0.0; - if (x > 0.0 && x < 1.0) - return INFINITY; - if (x < -1.0) - return 0.0; - if (x > -1.0 && x < 0.0) + if (x != 0.0) return INFINITY; } if (x >= LDBL_MAX) { @@ -244,6 +239,7 @@ long double powl(long double x, long double y) } w = floorl(y); + /* Set iyflg to 1 if y is an integer. */ iyflg = 0; if (w == y) @@ -271,43 +267,33 @@ long double powl(long double x, long double y) return 0.0; } } - - - nflg = 0; /* flag = 1 if x<0 raised to integer power */ + nflg = 0; /* (x<0)**(odd int) */ if (x <= 0.0) { if (x == 0.0) { if (y < 0.0) { if (signbit(x) && yoddint) - return -INFINITY; - return INFINITY; + /* (-0.0)**(-odd int) = -inf, divbyzero */ + return -1.0/0.0; + /* (+-0.0)**(negative) = inf, divbyzero */ + return 1.0/0.0; } - if (y > 0.0) { - if (signbit(x) && yoddint) - return -0.0; - return 0.0; - } - if (y == 0.0) - return 1.0; /* 0**0 */ - return 0.0; /* 0**y */ + if (signbit(x) && yoddint) + return -0.0; + return 0.0; } if (iyflg == 0) return (x - x) / (x - x); /* (x<0)**(non-int) is NaN */ - nflg = 1; + /* (x<0)**(integer) */ + if (yoddint) + nflg = 1; /* negate result */ + x = -x; } - - /* Integer power of an integer. */ - if (iyflg) { - i = w; - w = floorl(x); - if (w == x && fabsl(y) < 32768.0) { - w = powil(x, (int)y); - return w; - } + /* (+integer)**(integer) */ + if (iyflg && floorl(x) == x && fabsl(y) < 32768.0) { + w = powil(x, (int)y); + return nflg ? -w : w; } - if (nflg) - x = fabsl(x); - /* separate significand from exponent */ x = frexpl(x, &i); e = i; @@ -354,9 +340,7 @@ long double powl(long double x, long double y) z += x; /* Compute exponent term of the base 2 logarithm. */ - w = -i; - // TODO: use w * 0x1p-5; - w = scalbnl(w, -LNXT); /* divide by NXT */ + w = -i / NXT; w += e; /* Now base 2 log of x is w + z. */ @@ -381,7 +365,7 @@ long double powl(long double x, long double y) H = Fb + Gb; Ha = reducl(H); - w = scalbnl( Ga+Ha, LNXT ); + w = (Ga + Ha) * NXT; /* Test the power of 2 for overflow */ if (w > MEXP) @@ -418,18 +402,8 @@ long double powl(long double x, long double y) z = z + w; 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 = 0.5*y; - w = floorl(w); - w = 2.0*w; - if (w != y) - z = -z; /* odd exponent */ - } - + if (nflg) + z = -z; return z; } @@ -439,15 +413,14 @@ static long double reducl(long double x) { long double t; - t = scalbnl(x, LNXT); + t = x * NXT; t = floorl(t); - t = scalbnl(t, -LNXT); + t = t / NXT; return t; } -/* powil.c - * - * Real raised to integer power, long double precision +/* + * Positive real raised to integer power, long double precision * * * SYNOPSIS: @@ -460,7 +433,7 @@ static long double reducl(long double x) * * DESCRIPTION: * - * Returns argument x raised to the nth power. + * Returns argument x>0 raised to the nth power. * The routine efficiently decomposes n as a sum of powers of * two. The desired power is a product of two-to-the-kth * powers of x. Thus to compute the 32767 power of x requires @@ -482,25 +455,11 @@ static long double powil(long double x, int nn) { long double ww, y; long double s; - int n, e, sign, asign, lx; - - if (x == 0.0) { - if (nn == 0) - return 1.0; - else if (nn < 0) - return LDBL_MAX; - return 0.0; - } + int n, e, sign, lx; if (nn == 0) return 1.0; - if (x < 0.0) { - asign = -1; - x = -x; - } else - asign = 0; - if (nn < 0) { sign = -1; n = -nn; @@ -539,10 +498,8 @@ static long double powil(long double x, int nn) /* First bit of the power */ if (n & 1) y = x; - else { + else y = 1.0; - asign = 0; - } ww = x; n >>= 1; @@ -553,8 +510,6 @@ static long double powil(long double x, int nn) n >>= 1; } - if (asign) - y = -y; /* odd power of negative number */ if (sign < 0) y = 1.0/y; return y; diff --git a/src/math/remquo.c b/src/math/remquo.c index 79c9a55e..e92984ed 100644 --- a/src/math/remquo.c +++ b/src/math/remquo.c @@ -35,7 +35,6 @@ double remquo(double x, double y, int *quo) hy &= 0x7fffffff; /* |y| */ /* purge off exception values */ - // FIXME: signed shift if ((hy|ly) == 0 || hx >= 0x7ff00000 || /* y = 0, or x not finite */ (hy|((ly|-ly)>>31)) > 0x7ff00000) /* or y is NaN */ return (x*y)/(x*y); |