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

[Adhemerval Zanella Netto]

On 11/09/25 18:02, Wilco Dijkstra wrote:
> 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.

Ack.

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

Right, I think it should ok and not create a nan call for all targets.

> 
> +  /* |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.

Ack.

> 
> +  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)?

I think you correct, this check seems unnecessary and it does not trigger any
regressions.  I will remove from the float version as well.

> 
> 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.

I did not see much difference neither on x86_64 (Ryzen 9) or aarch64 (N1) with
recent gcc.
> 
> 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.

I am not sure if we can remove it, since it to handle small 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.

I think you mean x == -0.0 and yeah I think it should better indeed.

> 
> Cheers,
> Wilco