Fix the double rounding problem with subnormals, and

remove the XXX comments, which no longer apply.

2 files changed, 36 insertions, 16 deletions

diff --git a/lib/msun/src/s_fma.c b/lib/msun/src/s_fma.c
index 5403070..c6d27bd 100644
--- a/lib/msun/src/s_fma.c
+++ b/lib/msun/src/s_fma.c
@@ -45,15 +45,8 @@ __FBSDID("$FreeBSD$");
  * to be stored in FP registers in order to avoid incorrect results.
  * This is the default on FreeBSD, but not on many other systems.
  *
- * Tests on an Itanium 2 indicate that the hardware's FMA instruction
- * is almost twice as fast as this implementation.  The hardware
- * instruction should be used on platforms that support it.
- *
- * XXX May incur an absolute error of 0x1p-1074 for subnormal results
- *     due to double rounding induced by the final scaling operation.
- *
- * XXX On machines supporting quad precision, we should use that, but
- *     see the caveat in s_fmaf.c.
+ * Hardware instructions should be used on architectures that support it,
+ * since this implementation will likely be several times slower.
  */
 #if LDBL_MANT_DIG != 113
 double
@@ -174,8 +167,23 @@ fma(double x, double y, double z)
 	s = r - c;
 	rr = (c - (r - s)) + (zs - s) + cc;
 
-	fesetround(oround);
-	return (ldexp(r + rr, ex + ey));
+	spread = ex + ey;
+	if (spread + ilogb(r) > -1023) {
+		fesetround(oround);
+		r = r + rr;
+	} else {
+		/*
+		 * The result is subnormal, so we round before scaling to
+		 * avoid double rounding.
+		 */
+		p = ldexp(copysign(0x1p-1022, r), -spread);
+		c = r + p;
+		s = c - r;
+		cc = (r - (c - s)) + (p - s) + rr;
+		fesetround(oround);
+		r = (c + cc) - p;
+	}
+	return (ldexp(r, spread));
 }
 #else	/* LDBL_MANT_DIG == 113 */
 /*
diff --git a/lib/msun/src/s_fmal.c b/lib/msun/src/s_fmal.c
index c1b6413..1604778 100644
--- a/lib/msun/src/s_fmal.c
+++ b/lib/msun/src/s_fmal.c
@@ -39,9 +39,6 @@ __FBSDID("$FreeBSD$");
  *
  *	Dekker, T.  A Floating-Point Technique for Extending the
  *	Available Precision.  Numer. Math. 18, 224-242 (1971).
- *
- * XXX May incur a small error for subnormal results due to double
- *     rounding induced by the final scaling operation.
  */
 long double
 fmal(long double x, long double y, long double z)
@@ -165,6 +162,21 @@ fmal(long double x, long double y, long double z)
 	s = r - c;
 	rr = (c - (r - s)) + (zs - s) + cc;
 
-	fesetround(oround);
-	return (ldexpl(r + rr, ex + ey));
+	spread = ex + ey;
+	if (spread + ilogbl(r) > -16383) {
+		fesetround(oround);
+		r = r + rr;
+	} else {
+		/*
+		 * The result is subnormal, so we round before scaling to
+		 * avoid double rounding.
+		 */
+		p = ldexpl(copysignl(0x1p-16382L, r), -spread);
+		c = r + p;
+		s = c - r;
+		cc = (r - (c - s)) + (p - s) + rr;
+		fesetround(oround);
+		r = (c + cc) - p;
+	}
+	return (ldexpl(r, spread));
 }