From 6a8cdf6fedaae32990f8fa3c347c703c2a99fa16 Mon Sep 17 00:00:00 2001 From: das Date: Thu, 3 Apr 2008 06:14:51 +0000 Subject: Fix some corner cases: - fma(x, y, z) returns z, not NaN, if z is infinite, x and y are finite, x*y overflows, and x*y and z have opposite signs. - fma(x, y, z) doesn't generate an overflow, underflow, or inexact exception if z is NaN or infinite, as per IEEE 754R. - If the rounding mode is set to FE_DOWNWARD, fma(1.0, 0.0, -0.0) is -0.0, not +0.0. --- lib/msun/src/s_fma.c | 15 ++++++++++----- lib/msun/src/s_fmal.c | 15 ++++++++++----- 2 files changed, 20 insertions(+), 10 deletions(-) (limited to 'lib/msun/src') diff --git a/lib/msun/src/s_fma.c b/lib/msun/src/s_fma.c index c6d27bd..ad1fc4a 100644 --- a/lib/msun/src/s_fma.c +++ b/lib/msun/src/s_fma.c @@ -60,14 +60,19 @@ fma(double x, double y, double z) int ex, ey, ez; int spread; - if (z == 0.0) - return (x * y); + /* + * Handle special cases. The order of operations and the particular + * return values here are crucial in handling special cases involving + * infinities, NaNs, overflows, and signed zeroes correctly. + */ if (x == 0.0 || y == 0.0) return (x * y + z); - - /* Results of frexp() are undefined for these cases. */ - if (!isfinite(x) || !isfinite(y) || !isfinite(z)) + if (z == 0.0) + return (x * y); + if (!isfinite(x) || !isfinite(y)) return (x * y + z); + if (!isfinite(z)) + return (z); xs = frexp(x, &ex); ys = frexp(y, &ey); diff --git a/lib/msun/src/s_fmal.c b/lib/msun/src/s_fmal.c index 1604778..4d5d114 100644 --- a/lib/msun/src/s_fmal.c +++ b/lib/msun/src/s_fmal.c @@ -55,14 +55,19 @@ fmal(long double x, long double y, long double z) int ex, ey, ez; int spread; - if (z == 0.0) - return (x * y); + /* + * Handle special cases. The order of operations and the particular + * return values here are crucial in handling special cases involving + * infinities, NaNs, overflows, and signed zeroes correctly. + */ if (x == 0.0 || y == 0.0) return (x * y + z); - - /* Results of frexp() are undefined for these cases. */ - if (!isfinite(x) || !isfinite(y) || !isfinite(z)) + if (z == 0.0) + return (x * y); + if (!isfinite(x) || !isfinite(y)) return (x * y + z); + if (!isfinite(z)) + return (z); xs = frexpl(x, &ex); ys = frexpl(y, &ey); -- cgit v1.1