dns response handling: don't treat too many addresses as an error

[musl] / src / math / fmaf.c

diff --git a/src/math/fmaf.c b/src/math/fmaf.c
index 745ee39..7c65acf 100644 (file)
--- 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.
@@ -39,21 +40,35 @@ 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 */
        {
                /*
-               TODO: underflow is not raised correctly, example in
-               downward rouding: fmaf(0x1p-120f, 0x1p-120f, 0x1p-149f)
+               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;
        }
@@ -62,14 +77,16 @@ float fmaf(float x, float y, float 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);
-       z = adjusted_result;
+       double err;
+       int neg = u.i >> 63;
+       if (neg == (z > xy))
+               err = xy - result + z;
+       else
+               err = z - result + xy;
+       if (neg == (err < 0))
+               u.i++;
+       else
+               u.i--;
+       z = u.f;
        return z;
 }