From 1e2fea632b0257745ece5379378a1485f119602a Mon Sep 17 00:00:00 2001 From: nsz Date: Tue, 20 Mar 2012 15:17:15 +0100 Subject: fix a cbrtl.c regression and remove x87 precision setting --- src/math/cbrtl.c | 33 ++++++--------------------------- 1 file changed, 6 insertions(+), 27 deletions(-) (limited to 'src') 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 -- cgit v1.2.1 From 615bbd365f322cd8d0b1a4383ee0a6ce8e35b33b Mon Sep 17 00:00:00 2001 From: nsz Date: Tue, 20 Mar 2012 19:59:50 +0100 Subject: clean up powl.c fix special cases, use multiplication instead of scalbnl --- src/math/powl.c | 139 +++++++++++++++++++------------------------------------- 1 file changed, 47 insertions(+), 92 deletions(-) (limited to 'src') 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; -- cgit v1.2.1 From f1347a3a453ce38bf55cd92f5922311235810fd1 Mon Sep 17 00:00:00 2001 From: nsz Date: Tue, 20 Mar 2012 20:04:53 +0100 Subject: clean up pow.c and powf.c fix comments about special cases --- src/math/pow.c | 65 +++++++++++++++++++++++++-------------------------------- src/math/powf.c | 35 +++++++++++++------------------ 2 files changed, 43 insertions(+), 57 deletions(-) (limited to 'src') 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 */ -- cgit v1.2.1 From 8c6fc860a97f79146bf5c092d5cfb90fa6d9355a Mon Sep 17 00:00:00 2001 From: nsz Date: Tue, 20 Mar 2012 20:08:35 +0100 Subject: remove a fixme comment --- src/math/__signbit.c | 2 +- src/math/__signbitf.c | 2 +- src/math/remquo.c | 1 - 3 files changed, 2 insertions(+), 3 deletions(-) (limited to 'src') 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/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); -- cgit v1.2.1 From 91c28f61f43ba029166772e8ac25808ea3c3dc98 Mon Sep 17 00:00:00 2001 From: nsz Date: Tue, 20 Mar 2012 22:49:19 +0100 Subject: nearbyint optimization (only clear inexact when necessary) old code saved/restored the fenv (the new code is only as slow as that when inexact is not set before the call, but some other flag is set and the rounding is inexact, which is rare) before: bench_nearbyint_exact 5000000 N 261 ns/op bench_nearbyint_inexact_set 5000000 N 262 ns/op bench_nearbyint_inexact_unset 5000000 N 261 ns/op after: bench_nearbyint_exact 10000000 N 94.99 ns/op bench_nearbyint_inexact_set 25000000 N 65.81 ns/op bench_nearbyint_inexact_unset 10000000 N 94.97 ns/op --- src/math/nearbyint.c | 23 +++++++++++------------ src/math/nearbyintf.c | 14 ++++++++++---- src/math/nearbyintl.c | 11 ++++++++--- 3 files changed, 29 insertions(+), 19 deletions(-) (limited to 'src') 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 #include -/* -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 #include -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 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 -- cgit v1.2.1