[PATCH] fix underflow exception in fma and fmal

[Szabolcs Nagy <nsz@port70.net>]

another corner case in the freebsd fma code where signaling underflow
may be missed for an inexact subnormal result.
(fmaf and x86 fma are not affected)
---
 src/math/fma.c  | 7 +++++++
 src/math/fmal.c | 8 ++++++++
 2 files changed, 15 insertions(+)

diff --git a/src/math/fma.c b/src/math/fma.c
index 741ccd75..c69918d1 100644
--- a/src/math/fma.c
+++ b/src/math/fma.c
@@ -279,6 +279,13 @@ static inline double add_and_denormalize(double a, double b, int scale)
 			uhi.i += 1 - (((uhi.i ^ ulo.i) >> 62) & 2);
 			sum.hi = uhi.f;
 		}
+#ifdef FE_UNDERFLOW
+		/*
+		 * Raise underflow manually because scalbn won't do it if all
+		 * lost bits are 0: fma(-0x1p-1000, 0x1.000001p-74, 0x1p-1022)
+		 */
+		feraiseexcept(FE_UNDERFLOW);
+#endif
 	}
 	return scalbn(sum.hi, scale);
 }
diff --git a/src/math/fmal.c b/src/math/fmal.c
index 4506aac6..feb19ec1 100644
--- a/src/math/fmal.c
+++ b/src/math/fmal.c
@@ -121,6 +121,14 @@ static inline long double add_and_denormalize(long double a, long double b, int
 		bits_lost = -u.i.se - scale + 1;
 		if ((bits_lost != 1) ^ LASTBIT(u))
 			sum.hi = nextafterl(sum.hi, INFINITY * sum.lo);
+#ifdef FE_UNDERFLOW
+		/*
+		 * Raise underflow manually because scalbnl won't do it if all
+		 * lost bits are 0, e.g.:
+		 * fmal(-0x1p-10000L, 0x1.0000000000001p-6445L, 0x1p-16382L)
+		 */
+		feraiseexcept(FE_UNDERFLOW);
+#endif
 	}
 	return scalbnl(sum.hi, scale);
 }
-- 
2.11.0