From 8f438115f2c12b40fa7b1884f87db72857af67f6 Mon Sep 17 00:00:00 2001
From: Szabolcs Nagy <nsz@port70.net>
Date: Mon, 7 Oct 2013 18:46:43 +0000
Subject: math: fix rare underflow issue in fma

the issue is described in commits 1e5eb73545ca6cfe8b918798835aaf6e07af5beb
and ffd8ac2dd50f99c3c83d7d9d845df9874ec3e7d5
---
 src/math/fmaf.c | 36 +++++++++++++++++++++++++++---------
 1 file changed, 27 insertions(+), 9 deletions(-)

(limited to 'src/math/fmaf.c')

diff --git a/src/math/fmaf.c b/src/math/fmaf.c
index 745ee393..aa57feb6 100644
--- a/src/math/fmaf.c
+++ b/src/math/fmaf.c
@@ -26,7 +26,8 @@
  */
 
 #include <fenv.h>
-#include "libm.h"
+#include <math.h>
+#include <stdint.h>
 
 /*
  * Fused multiply-add: Compute x * y + z with a single rounding error.
@@ -39,21 +40,35 @@ float fmaf(float x, float y, float z)
 {
 	#pragma STDC FENV_ACCESS ON
 	double xy, result;
-	uint32_t hr, lr;
+	union {double f; uint64_t i;} u;
+	int e;
 
 	xy = (double)x * y;
 	result = xy + z;
-	EXTRACT_WORDS(hr, lr, result);
+	u.f = result;
+	e = u.i>>52 & 0x7ff;
 	/* Common case: The double precision result is fine. */
-	if ((lr & 0x1fffffff) != 0x10000000 ||  /* not a halfway case */
-		(hr & 0x7ff00000) == 0x7ff00000 ||  /* NaN */
+	if ((u.i & 0x1fffffff) != 0x10000000 || /* not a halfway case */
+		e == 0x7ff ||                   /* NaN */
 		result - xy == z ||                 /* exact */
 		fegetround() != FE_TONEAREST)       /* not round-to-nearest */
 	{
 		/*
-		TODO: underflow is not raised correctly, example in
-		downward rouding: fmaf(0x1p-120f, 0x1p-120f, 0x1p-149f)
+		underflow may not be raised correctly, example:
+		fmaf(0x1p-120f, 0x1p-120f, 0x1p-149f)
 		*/
+#if defined(FE_INEXACT) && defined(FE_UNDERFLOW)
+		if (e < 0x3ff-126 && e >= 0x3ff-149 && fetestexcept(FE_INEXACT)) {
+			feclearexcept(FE_INEXACT);
+			/* TODO: gcc and clang bug workaround */
+			volatile float vz = z;
+			result = xy + vz;
+			if (fetestexcept(FE_INEXACT))
+				feraiseexcept(FE_UNDERFLOW);
+			else
+				feraiseexcept(FE_INEXACT);
+		}
+#endif
 		z = result;
 		return z;
 	}
@@ -68,8 +83,11 @@ float fmaf(float x, float y, float z)
 	volatile double vxy = xy;  /* XXX work around gcc CSE bug */
 	double adjusted_result = vxy + z;
 	fesetround(FE_TONEAREST);
-	if (result == adjusted_result)
-		SET_LOW_WORD(adjusted_result, lr + 1);
+	if (result == adjusted_result) {
+		u.f = adjusted_result;
+		u.i++;
+		adjusted_result = u.f;
+	}
 	z = adjusted_result;
 	return z;
 }
-- 
cgit v1.2.1