[PATCH v2 3/4] math: Optimize dbl-64 remainder implementation

[Wilco Dijkstra]

Hi Adhemerval,

+  uint64_t hx = asuint64 (x);
+  uint64_t hy = asuint64 (y);
+  uint64_t sx = hx >> 63;
+
+  hx &= UINT64_C (0x7FFFFFFFFFFFFFFF);
+  hy &= UINT64_C (0x7FFFFFFFFFFFFFFF);

hx &= ~SIGN_MASK?

+  /* |x| not finite or |y| is NaN or 0 */
+  if ((hx >= EXPONENT_MASK || (hy - 1) >= EXPONENT_MASK))
+    return (x * y) / ( x * y);

We can move this into the else below and just check hx:

+  /* |y| < DBL_MAX / 2 ? */
+  y = fabs (y);
+  if (__glibc_likely (hy < UINT64_C (0x7fe0000000000000)))

We can add hx < EXPONENT_MASK here to skip the NaN checks above
(the y == 0.0 case is handled correctly by fmod call).

+  /* Make sure x is not -0. This can occur only when x = p  and rounding
+     direction is towards negative infinity.  */
+  if (__glibc_unlikely (x == -0.0))
+    x = 0.0;

This still seems to work out badly - I did some experiments, this moves the
zero check to the case that requires it (after the first x - y):

          x -= y;
          if (x + x >= y)
            x -= y;
          else if (x == 0.0)
            x = 0.0;

This gives 5-15% improvements to throughput and latency for the
denormal and close exponent cases.

This is probably as good as it gets without merging fmod and remainder
into a shared implementation.

Cheers,
Wilco