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.

2 files changed, 20 insertions, 10 deletions

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);