X-Git-Url: http://nsz.repo.hu/git/?p=musl;a=blobdiff_plain;f=src%2Fmath%2Ffma.c;h=850a4a6c7f0dcc164521082be8fef0fac986c11f;hp=5fb9540601c520dac6ce9a519e8c9abe91aed7d7;hb=6a4cfbdbe718a115a22629ad0cb2ae21391a0454;hpb=2786c7d21611b9fa3b2fe356542cf213e7dd0ba4

diff --git a/src/math/fma.c b/src/math/fma.c
index 5fb95406..850a4a6c 100644
--- a/src/math/fma.c
+++ b/src/math/fma.c
@@ -41,7 +41,7 @@ static void mul(long double *hi, long double *lo, long double x, long double y)
 
 /*
 assume (long double)(hi+lo) == hi
-return an adjusted hi so that rounding it to double is correct
+return an adjusted hi so that rounding it to double (or less) precision is correct
 */
 static long double adjust(long double hi, long double lo)
 {
@@ -55,17 +55,25 @@ static long double adjust(long double hi, long double lo)
 	ulo.x = lo;
 	if (uhi.bits.s == ulo.bits.s)
 		uhi.bits.m++;
-	else
+	else {
 		uhi.bits.m--;
+		/* handle underflow and take care of ld80 implicit msb */
+		if (uhi.bits.m == (uint64_t)-1/2) {
+			uhi.bits.m *= 2;
+			uhi.bits.e--;
+		}
+	}
 	return uhi.x;
 }
 
+/* adjusted add so the result is correct when rounded to double (or less) precision */
 static long double dadd(long double x, long double y)
 {
 	add(&x, &y, x, y);
 	return adjust(x, y);
 }
 
+/* adjusted mul so the result is correct when rounded to double (or less) precision */
 static long double dmul(long double x, long double y)
 {
 	mul(&x, &y, x, y);
@@ -81,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;
 
@@ -110,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 {
 		/*
@@ -126,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 */
@@ -298,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;
@@ -397,6 +441,8 @@ double fma(double x, double y, double z)
 		/*
 		 * There is no need to worry about double rounding in directed
 		 * rounding modes.
+		 * TODO: underflow is not raised properly, example in downward rounding:
+		 * fma(0x1.000000001p-1000, 0x1.000000001p-30, -0x1p-1066)
 		 */
 		fesetround(oround);
 		adj = r.lo + xy.lo;