math: fix two fma issues (only affects non-nearest rounding mode, x86)

[musl] / src / math / fma.c

diff --git a/src/math/fma.c b/src/math/fma.c
index 87d450c..89def79 100644 (file)
--- 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 */
@@ -247,7 +290,7 @@ static inline double add_and_denormalize(double a, double b, int scale)
                        INSERT_WORD64(sum.hi, hibits);
                }
        }
-       return (ldexp(sum.hi, scale));
+       return scalbn(sum.hi, scale);
 }
 
 /*
@@ -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;
@@ -364,7 +408,7 @@ double fma(double x, double y, double z)
                }
        }
        if (spread <= DBL_MANT_DIG * 2)
-               zs = ldexp(zs, -spread);
+               zs = scalbn(zs, -spread);
        else
                zs = copysign(DBL_MIN, zs);
 
@@ -390,7 +434,7 @@ double fma(double x, double y, double z)
                 */
                fesetround(oround);
                volatile double vzs = zs; /* XXX gcc CSE bug workaround */
-               return (xy.hi + vzs + ldexp(xy.lo, spread));
+               return xy.hi + vzs + scalbn(xy.lo, spread);
        }
 
        if (oround != FE_TONEAREST) {
@@ -400,13 +444,13 @@ double fma(double x, double y, double z)
                 */
                fesetround(oround);
                adj = r.lo + xy.lo;
-               return (ldexp(r.hi + adj, spread));
+               return scalbn(r.hi + adj, spread);
        }
 
        adj = add_adjusted(r.lo, xy.lo);
        if (spread + ilogb(r.hi) > -1023)
-               return (ldexp(r.hi + adj, spread));
+               return scalbn(r.hi + adj, spread);
        else
-               return (add_and_denormalize(r.hi, adj, spread));
+               return add_and_denormalize(r.hi, adj, spread);
 }
 #endif