diff options
author | bde <bde@FreeBSD.org> | 2007-05-02 15:24:49 +0000 |
---|---|---|
committer | bde <bde@FreeBSD.org> | 2007-05-02 15:24:49 +0000 |
commit | f9978195dcf7995705c06b7f08567ba6ddcaed8f (patch) | |
tree | 1b4dfe607c9ba9f340aef50d582778c8219ed5c5 /lib | |
parent | 98484906b64f28e7d3ec42156666cd1efd6b9998 (diff) | |
download | FreeBSD-src-f9978195dcf7995705c06b7f08567ba6ddcaed8f.zip FreeBSD-src-f9978195dcf7995705c06b7f08567ba6ddcaed8f.tar.gz |
Fix tgamma() on some special args:
(1) tgamma(-Inf) returned +Inf and failed to raise any exception, but
should always have raised an exception, and should behave like
tgamma(negative integer).
(2) tgamma(negative integer) returned +Inf and raised divide-by-zero,
but should return NaN and raise "invalid" on any IEEEish system.
(3) About half of the 2**52 negative intgers between -2**53 and -2**52
were misclassified as non-integers by using floor(x + 0.5) to round
to nearest, so tgamma(x) was wrong (+-0 instead of +Inf and now NaN)
on these args. The floor() expression is hard to use since rounding
of (x + 0.5) may give x or x + 1, depending on |x| and the current
rounding mode. The fixed version uses ceil(x) to classify x before
operating on x and ends up being more efficient since ceil(x) is
needed anyway.
(4) On at least the problematic args in (3), tgamma() raised a spurious
inexact.
(5) tgamma(large positive) raised divide-by-zero but should raise overflow.
(6) tgamma(+Inf) raised divide-by-zero but should not raise any exception.
(7) Raise inexact for tiny |x| in a way that has some chance of not being
optimized away.
The fix for (5) and (6), and probably for (2), also prevents -O optimizing
away the exception.
PR: 112180 (2)
Standards: Annex F in C99 (IEC 60559 binding) requires (1), (2) and (6).
Diffstat (limited to 'lib')
-rw-r--r-- | lib/msun/bsdsrc/b_tgamma.c | 29 |
1 files changed, 15 insertions, 14 deletions
diff --git a/lib/msun/bsdsrc/b_tgamma.c b/lib/msun/bsdsrc/b_tgamma.c index 784e563..5593bbd 100644 --- a/lib/msun/bsdsrc/b_tgamma.c +++ b/lib/msun/bsdsrc/b_tgamma.c @@ -49,7 +49,7 @@ __FBSDID("$FreeBSD$"); /* METHOD: * x < 0: Use reflection formula, G(x) = pi/(sin(pi*x)*x*G(x)) - * At negative integers, return +Inf and raise divide-by-zero. + * At negative integers, return NaN and raise invalid. * * x < 6.5: * Use argument reduction G(x+1) = xG(x) to reach the @@ -66,12 +66,12 @@ __FBSDID("$FreeBSD$"); * avoid premature round-off. * * Special values: - * -Inf: return +Inf (without raising any exception!); - * negative integer: return +Inf and raise divide-by-zero; + * -Inf: return NaN and raise invalid; + * negative integer: return NaN and raise invalid; * other x ~< 177.79: return +-0 and raise underflow; * +-0: return +-Inf and raise divide-by-zero; - * finite x ~> 171.63: return +Inf and raise divide-by-zero(!); - * +Inf: return +Inf and raise divide-by-zero(!); + * finite x ~> 171.63: return +Inf and raise overflow; + * +Inf: return +Inf; * NaN: return NaN. * * Accuracy: tgamma(x) is accurate to within @@ -135,7 +135,7 @@ tgamma(x) if (x >= 6) { if(x > 171.63) - return(one/zero); + return (x / zero); u = large_gam(x); return(__exp__D(u.a, u.b)); } else if (x >= 1.0 + LEFT + x0) @@ -143,12 +143,11 @@ tgamma(x) else if (x > 1.e-17) return (smaller_gam(x)); else if (x > -1.e-17) { - if (x == 0.0) - return (one/x); - one+1e-20; /* Raise inexact flag. */ + if (x != 0.0) + u.a = one - tiny; /* raise inexact */ return (one/x); } else if (!finite(x)) - return (x*x); /* x = NaN, -Inf */ + return (x - x); /* x is NaN or -Inf */ else return (neg_gam(x)); } @@ -282,11 +281,13 @@ neg_gam(x) struct Double lg, lsine; double y, z; - y = floor(x + .5); + y = ceil(x); if (y == x) /* Negative integer. */ - return (one/zero); - z = fabs(x - y); - y = .5*ceil(x); + return ((x - x) / zero); + z = y - x; + if (z > 0.5) + z = one - z; + y = 0.5 * y; if (y == ceil(y)) sgn = -1; if (z < .25) |