1
0
mirror of https://git.FreeBSD.org/src.git synced 2024-12-16 10:20:30 +00:00

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).
This commit is contained in:
Bruce Evans 2007-05-02 15:24:49 +00:00
parent dd936b27fc
commit e95cc9b700
Notes: svn2git 2020-12-20 02:59:44 +00:00
svn path=/head/; revision=169212

View File

@ -49,7 +49,7 @@ __FBSDID("$FreeBSD$");
/* METHOD: /* METHOD:
* x < 0: Use reflection formula, G(x) = pi/(sin(pi*x)*x*G(x)) * 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: * x < 6.5:
* Use argument reduction G(x+1) = xG(x) to reach the * Use argument reduction G(x+1) = xG(x) to reach the
@ -66,12 +66,12 @@ __FBSDID("$FreeBSD$");
* avoid premature round-off. * avoid premature round-off.
* *
* Special values: * Special values:
* -Inf: return +Inf (without raising any exception!); * -Inf: return NaN and raise invalid;
* negative integer: return +Inf and raise divide-by-zero; * negative integer: return NaN and raise invalid;
* other x ~< 177.79: return +-0 and raise underflow; * other x ~< 177.79: return +-0 and raise underflow;
* +-0: return +-Inf and raise divide-by-zero; * +-0: return +-Inf and raise divide-by-zero;
* finite x ~> 171.63: return +Inf and raise divide-by-zero(!); * finite x ~> 171.63: return +Inf and raise overflow;
* +Inf: return +Inf and raise divide-by-zero(!); * +Inf: return +Inf;
* NaN: return NaN. * NaN: return NaN.
* *
* Accuracy: tgamma(x) is accurate to within * Accuracy: tgamma(x) is accurate to within
@ -135,7 +135,7 @@ tgamma(x)
if (x >= 6) { if (x >= 6) {
if(x > 171.63) if(x > 171.63)
return(one/zero); return (x / zero);
u = large_gam(x); u = large_gam(x);
return(__exp__D(u.a, u.b)); return(__exp__D(u.a, u.b));
} else if (x >= 1.0 + LEFT + x0) } else if (x >= 1.0 + LEFT + x0)
@ -143,12 +143,11 @@ tgamma(x)
else if (x > 1.e-17) else if (x > 1.e-17)
return (smaller_gam(x)); return (smaller_gam(x));
else if (x > -1.e-17) { else if (x > -1.e-17) {
if (x == 0.0) if (x != 0.0)
return (one/x); u.a = one - tiny; /* raise inexact */
one+1e-20; /* Raise inexact flag. */
return (one/x); return (one/x);
} else if (!finite(x)) } else if (!finite(x))
return (x*x); /* x = NaN, -Inf */ return (x - x); /* x is NaN or -Inf */
else else
return (neg_gam(x)); return (neg_gam(x));
} }
@ -282,11 +281,13 @@ neg_gam(x)
struct Double lg, lsine; struct Double lg, lsine;
double y, z; double y, z;
y = floor(x + .5); y = ceil(x);
if (y == x) /* Negative integer. */ if (y == x) /* Negative integer. */
return (one/zero); return ((x - x) / zero);
z = fabs(x - y); z = y - x;
y = .5*ceil(x); if (z > 0.5)
z = one - z;
y = 0.5 * y;
if (y == ceil(y)) if (y == ceil(y))
sgn = -1; sgn = -1;
if (z < .25) if (z < .25)