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

[Wilco Dijkstra]

Hi Adhemerval,

+  uint64_t hx = asuint64 (x);
+  uint64_t hy = asuint64 (y);
+  uint64_t sx = hx & UINT64_C (0x8000000000000000);

We below change, could become sign = hx >> 63;

+  hx &= UINT64_C (0x7FFFFFFFFFFFFFFF);
+  hy &= UINT64_C (0x7FFFFFFFFFFFFFFF);
+
+  /* |y| == 0 or |x| not finite or |y| is NaN */
+  if (hy == 0
+      || (hx >= EXPONENT_MASK || hy > EXPONENT_MASK))

Surely (hy - 1) >= EXPONENT_MASK to remove the hy == 0 check.

+    return NAN;

Since NAN returns a float NaN, wouldn't C99 nan("") be better?

+  /* |y| < DBL/MAX / 2 ? */
+  if (hy < UINT64_C (0x7fdfffffffffffff))
+    x = __ieee754_fmod (x, y + y);

Can use hy <=, but hy < UINT64_C (0x7fe0000000000000) is a simpler constant
on various targets.

+  if ((hx - hy) == 0)
+    return 0.0 * x;

Is this strictly required? I think the code below will also return 0.0 with sign of
original x, or does the correction for rounding mode make it fail? Does
remainder (x, x) return a different sign of zero than remainder (x + x , x)?

If it is needed, placing if (hx == hy) before the fmod call would be better.

+  x = fabs (x);
+  y = fabs (y);
+
+  /* |y| < 2 * DBL_MIN */
+  if (hy < UINT64_C (0x20000000000000))
     {
+      if (x + x > y)
         {
+         x -= y;
+         if (x + x >= y)
+           x -= y;
         }
     }
+  else
+    {
+      double y_half = y * 0.5;
+      if (x > y_half)
+       {
+         x -= y;
+         if (x >= y_half)
+           x -= y;
+       }
+    }

Since these branches may be fairly hard to predict, it's worth checking whether
this works out better:

x = (x >= y) ? x - y : x;
x = ((x + x) >= y) ? x - y : x;

You can always start with x >= y and only use (x + x) or y * 0.5 for the 2nd part.

In either case, we can remove the if (hy < UINT64_C (0x20000000000000)) by adding
an else after the fmod to handle the case of large y.

+  hx = asuint64 (x);
+  /* Make sure x is not -0. This can occur only when x = p  and rounding
+     direction is towards negative infinity.  */
+  if (hx == UINT64_C (0x8000000000000000))
+    hx = 0;
+  return asdouble (hx ^ sx);

Isn't this equivalent to:

if (__glibc_unlikely (x == 0.0))
  x = 0.0;
return sx ? -x : x;

That avoids the float-int-float roundtrip and large constants, so that might be faster
overall.

Cheers,
Wilco