X-Git-Url: http://nsz.repo.hu/git/?a=blobdiff_plain;f=src%2Fmath%2Ffmaf.c;h=80f5cd8a332fff28bb3022d6303546d216c78985;hb=2fab90a71acd3698954c08b9062db67188443dd7;hp=0dccf108471875e048b210e31408ef9525a23744;hpb=b69f695acedd4ce2798ef9ea28d834ceccc789bd;p=musl

diff --git a/src/math/fmaf.c b/src/math/fmaf.c
index 0dccf108..80f5cd8a 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.
@@ -37,28 +38,56 @@
  */
 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 */
-		result - xy == z ||                 /* exact */
+	if ((u.i & 0x1fffffff) != 0x10000000 || /* not a halfway case */
+		e == 0x7ff ||                   /* NaN */
+		(result - xy == z && result - z == xy) || /* exact */
 		fegetround() != FE_TONEAREST)       /* not round-to-nearest */
-		return (result);
+	{
+		/*
+		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;
+	}
 
 	/*
 	 * If result is inexact, and exactly halfway between two float values,
 	 * we need to adjust the low-order bit in the direction of the error.
 	 */
+#ifdef FE_TOWARDZERO
 	fesetround(FE_TOWARDZERO);
+#endif
 	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);
-	return (adjusted_result);
+	if (result == adjusted_result) {
+		u.f = adjusted_result;
+		u.i++;
+		adjusted_result = u.f;
+	}
+	z = adjusted_result;
+	return z;
 }