[PATCH 5/5] math: New generic fmaf implementation

[Wilco Dijkstra]

Hi Adhemerval,


+  uint64_t u = asuint64 (result);
+  /* Common case: The double precision result is fine. */
+  if ((u & 0x1fffffff) != 0x10000000 ||          /* not a halfway case */
+       (result - xy == z && result - z == xy) || /* exact */
+        __fegetround () != FE_TONEAREST)         /* not round-to-nearest */
+    {
+      /* Underflow may not be raised correctly, example:
+        fmaf(0x1p-120f, 0x1p-120f, 0x1p-149f)  */
+      int e = u >> MANTISSA_WIDTH & 0x7ff;
+      if (__glibc_unlikely (e <= EXPONENT_BIAS - 126
+                           && e >= EXPONENT_BIAS - 149))
+       return fmaf_fallback (xy, z);
+      return result;
+    }

So that is a bit complex... We don't need the fmaf_fallback code at all - and we
should avoid that kind of code at any cost. It is at least 10 times slower than the
alternative...

We should return immediately if the result is inexact, and handle the IEEE underflow
case by reusing the same fallback code:

  if (__glibc_likely ((u & 0xfffffff) != 0))
     return result;

  if ((u & 0x10000000) == 0
      && ((u >> MANTISSA_WIDTH) & 0x7ff) > EXPONENT_BIAS - 126)
    return result;

  if (result - xy == z && result - z == xy)
    return result;

Now do the standard correction for exact results. No need for checking the rounding
mode as it doesn't improve the common case.

This passes the test suite and results in a 20% speedup.

Cheers,
Wilco