X-Git-Url: http://nsz.repo.hu/git/?p=musl;a=blobdiff_plain;f=src%2Fmath%2Ffma.c;h=89def795943573d1568c5328ee95b2c9efc27b36;hp=7076d4b9308fbf9b7a0c60801b347c1c5e6c4e40;hb=ffd8ac2dd50f99c3c83d7d9d845df9874ec3e7d5;hpb=e5fb6820a42a1f675ba09c15273953e1ace65777

diff --git a/src/math/fma.c b/src/math/fma.c
index 7076d4b9..89def795 100644
--- a/src/math/fma.c
+++ b/src/math/fma.c
@@ -89,6 +89,7 @@ static int getexp(long double x)
 
 double fma(double x, double y, double z)
 {
+	#pragma STDC FENV_ACCESS ON
 	long double hi, lo1, lo2, xy;
 	int round, ez, exy;
 
@@ -118,9 +119,17 @@ double fma(double x, double y, double z)
 	} else if (ez > exy - 12) {
 		add(&hi, &lo2, xy, z);
 		if (hi == 0) {
+			/*
+			xy + z is 0, but it should be calculated with the
+			original rounding mode so the sign is correct, if the
+			compiler does not support FENV_ACCESS ON it does not
+			know about the changed rounding mode and eliminates
+			the xy + z below without the volatile memory access
+			*/
+			volatile double z_;
 			fesetround(round);
-			/* make sure that the sign of 0 is correct */
-			return (xy + z) + lo1;
+			z_ = z;
+			return (xy + z_) + lo1;
 		}
 	} else {
 		/*
@@ -134,10 +143,36 @@ double fma(double x, double y, double z)
 		hi = xy;
 		lo2 = z;
 	}
+	/*
+	the result is stored before return for correct precision and exceptions
+
+	one corner case is when the underflow flag should be raised because
+	the precise result is an inexact subnormal double, but the calculated
+	long double result is an exact subnormal double
+	(so rounding to double does not raise exceptions)
+
+	in nearest rounding mode dadd takes care of this: the last bit of the
+	result is adjusted so rounding sees an inexact value when it should
+
+	in non-nearest rounding mode fenv is used for the workaround
+	*/
 	fesetround(round);
 	if (round == FE_TONEAREST)
-		return dadd(hi, dadd(lo1, lo2));
-	return hi + (lo1 + lo2);
+		z = dadd(hi, dadd(lo1, lo2));
+	else {
+#if defined(FE_INEXACT) && defined(FE_UNDERFLOW)
+		int e = fetestexcept(FE_INEXACT);
+		feclearexcept(FE_INEXACT);
+#endif
+		z = hi + (lo1 + lo2);
+#if defined(FE_INEXACT) && defined(FE_UNDERFLOW)
+		if (getexp(z) < 0x3fff-1022 && fetestexcept(FE_INEXACT))
+			feraiseexcept(FE_UNDERFLOW);
+		else if (e)
+			feraiseexcept(FE_INEXACT);
+#endif
+	}
+	return z;
 }
 #else
 /* origin: FreeBSD /usr/src/lib/msun/src/s_fma.c */
@@ -306,6 +341,7 @@ static inline struct dd dd_mul(double a, double b)
  */
 double fma(double x, double y, double z)
 {
+	#pragma STDC FENV_ACCESS ON
 	double xs, ys, zs, adj;
 	struct dd xy, r;
 	int oround;