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

diff --git a/src/math/fmaf.c b/src/math/fmaf.c
index 0dccf108..745ee393 100644
--- a/src/math/fmaf.c
+++ b/src/math/fmaf.c
@@ -37,6 +37,7 @@
  */
 float fmaf(float x, float y, float z)
 {
+	#pragma STDC FENV_ACCESS ON
 	double xy, result;
 	uint32_t hr, lr;
 
@@ -48,17 +49,27 @@ float fmaf(float x, float y, float z)
 		(hr & 0x7ff00000) == 0x7ff00000 ||  /* NaN */
 		result - xy == z ||                 /* exact */
 		fegetround() != FE_TONEAREST)       /* not round-to-nearest */
-		return (result);
+	{
+		/*
+		TODO: underflow is not raised correctly, example in
+		downward rouding: fmaf(0x1p-120f, 0x1p-120f, 0x1p-149f)
+		*/
+		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);
+	z = adjusted_result;
+	return z;
 }