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

[Adhemerval Zanella Netto]

On 29/09/25 15:26, Wilco Dijkstra wrote:
> 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?

Ack.

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

Yeah, we can simplify the check for the likely case.  But we still need
the fully check on the else case to avoid FE_INVALID for qNaN inputs
(triggered by the FP operations).

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

We need on both, the lower case is required for the corner case

  remainder_downward (max_value, max_value)
  remainder_downward (max_value, -max_value)
  remainder_downward (-max_value, max_value)
  remainder_downward (-max_value, -max_value)

And it is still an improvement.

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

Thanks, I will send an updated version.