[PATCH 4/4] math: Improve fmodf

[Adhemerval Zanella Netto]

On 14/03/23 13:42, Wilco Dijkstra wrote:
> Hi Adhemerval,
> 
> This looks good overall. I guess we still needs error handling and the wrapper
> disabled since it spends 1/3 of the time in the useless wrapper but that can be
> a separate patch.

The wrapper should be simple to add.

> 
> A few comments below, there are a few things that look wrong, and I think most
> of the helper functions just add extra complexity and branches for no gain:
> 
> 
> +  uint32_t sx = hx & SIGN_MASK;
> +  /* Get |x| and |y|.  */
> +  hx ^= sx;
> +  hy &= ~SIGN_MASK;
> +
> +  /* Special cases:
> +     - If x or y is a Nan, NaN is returned.
> +     - If x is an inifinity, a NaN is returned.
> +     - If y is zero, Nan is returned.
> +     - If x is +0/-0, and y is not zero, +0/-0 is returned.  */
> +  if (__glibc_unlikely (hy == 0        || hx >= EXPONENT_MASK || hy > EXPONENT_MASK))
> +    return (x * y) / (x * y);
> +
> +  if (__glibc_unlikely (hx <= hy))
> +    {
> +      if (hx < hy)
> +       return x;
> +      return sx ? -0.0 : 0.0;
> 
> This should be return asfloat (sx);

Ack.

> 
> +    }
> +
> +  int ex = get_unbiased_exponent (hx);
> +  int ey = get_unbiased_exponent (hy);
> 
> Should be hx >> MANTISSA_WIDTH since we cleared the sign bits (now we don't
> need to add get_unbiased_exponent).

Ack.

> 
> +  /* Common case where exponents are close: ey >= -103 and |x/y| < 2^8,  */
> +  if (__glibc_likely (ey > MANTISSA_WIDTH && ex - ey <= EXPONENT_WIDTH))
> +    {
> +      uint32_t mx = get_explicit_mantissa (hx);
> +      uint32_t my = get_explicit_mantissa (hy);
> 
> Note this is equivalent to:
> 
> mx = (hx & MANTISSA_MASK) | (MANTISSA_MASK + 1);

Ack.

> 
> So we don't need get_explicit_mantissa (or we could change it to do the above).
> If we do this before the if statement, we don't need to repeat it below.
> 
> +      uint32_t d = (ex == ey) ? (mx - my) : (mx << (ex - ey)) % my;
> +      if (d == 0)
> +       return 0.0;
> 
> Looks like a bug, should be asfloat (sx)?

In fact I think this check is not really required since we already if
inputs are equal.

> 
> +      return make_float (d, ey - 1, sx);
> +    }
> +
> +  /* Special case, both x and y are subnormal.  */
> +  if (__glibc_unlikely (ex == 0 && ey == 0))
> +    return asfloat (hx % hy);
> 
> Similarly, shouldn't this be asfloat (sx | (hx % hy))?

Yes, I will add tests to check for it.

> 
> +  /* Convert |x| and |y| to 'mx + 2^ex' and 'my + 2^ey'.  Assume that hx is
> +     not subnormal by conditions above.  */
> +  uint32_t mx = get_explicit_mantissa (hx);
> +  ex--;
> +
> +  uint32_t my = get_explicit_mantissa (hy);
> 
> If we set mx/my above then this isn't needed.

Right, I think we can mask out the mantissa and set the implicit bit then.

> 
> +  int lead_zeros_my = EXPONENT_WIDTH;
> +  if (__glibc_likely (ey > 0))
> +    ey--;
> +  else
> +    {
> +      my = get_mantissa (hy);
> 
> This is really my = hy; since we know hy is positive denormal.

Indeed.

> 
> +      lead_zeros_my = __builtin_clz (my);
> +    }
> +
> +  int tail_zeros_my = __builtin_ctz (my);
> +  int sides_zeroes = lead_zeros_my + tail_zeros_my;
> +  int exp_diff = ex - ey;
> +
> +  int right_shift = exp_diff < tail_zeros_my ? exp_diff : tail_zeros_my;
> +  my >>= right_shift;
> +  exp_diff -= right_shift;
> +  ey += right_shift;
> +
> +  int left_shift = exp_diff < EXPONENT_WIDTH ? exp_diff : EXPONENT_WIDTH;
> +  mx <<= left_shift;
> +  exp_diff -= left_shift;
> +
> +  mx %= my;
> +
> +  if (__glibc_unlikely (mx == 0))
> +    return sx ? -0.0 : 0.0;
> 
> Should be asfloat (sx);

Ack.

> 
> +  if (exp_diff == 0)
> +    return make_float (my, ey, sx);
> +
> +  /* Assume modulo/divide operation is slow, so use multiplication with invert
> +     values.  */
> +  uint32_t inv_hy = UINT32_MAX / my;
> +  while (exp_diff > sides_zeroes) {
> +    exp_diff -= sides_zeroes;
> +    uint32_t hd = (mx * inv_hy) >> (BIT_WIDTH - sides_zeroes);
> +    mx <<= sides_zeroes;
> +    mx -= hd * my;
> +    while (__glibc_unlikely (mx > my))
> +      mx -= my;
> +  }
> +  uint32_t hd = (mx * inv_hy) >> (BIT_WIDTH - exp_diff);
> +  mx <<= exp_diff;
> +  mx -= hd * my;
> +  while (__glibc_unlikely (mx > my))
> +    mx -= my;
> +
> +  return make_float (mx, ey, sx);
>  }
>  libm_alias_finite (__ieee754_fmodf, __fmodf)
> diff --git a/sysdeps/ieee754/flt-32/math_config.h b/sysdeps/ieee754/flt-32/math_config.h
> index 23045f59d6..cdab3a36ef 100644
> --- a/sysdeps/ieee754/flt-32/math_config.h
> +++ b/sysdeps/ieee754/flt-32/math_config.h
> @@ -110,6 +110,95 @@ issignalingf_inline (float x)
>    return 2 * (ix ^ 0x00400000) > 2 * 0x7fc00000UL;
>  }
>  
> +#define BIT_WIDTH       32
> +#define MANTISSA_WIDTH  23
> +#define EXPONENT_WIDTH  8
> +#define MANTISSA_MASK   0x007fffff
> +#define EXPONENT_MASK   0x7f800000
> +#define EXP_MANT_MASK   0x7fffffff
> +#define QUIET_NAN_MASK  0x00400000
> +#define SIGN_MASK       0x80000000
> +
> +static inline bool
> +is_nan (uint32_t x)
> +{
> +  return (x & EXP_MANT_MASK) > EXPONENT_MASK;
> +}
> +
> +static inline bool
> +is_quiet_nan (uint32_t x)
> +{
> +   return (x & EXP_MANT_MASK) == (EXPONENT_MASK | QUIET_NAN_MASK);
> +}
> +
> +static inline bool
> +is_inf_or_nan (uint32_t x)
> +{
> +  return (x & EXPONENT_MASK) == EXPONENT_MASK;
> +}
> +
> +static inline uint16_t
> +get_unbiased_exponent (uint32_t x)
> +{
> +  return (x & EXPONENT_MASK) >> MANTISSA_WIDTH;
> +}
> +
> +/* Return mantissa with the implicit bit set iff X is a normal number.  */
> +static inline uint32_t
> +get_explicit_mantissa (uint32_t x)
> +{
> +  uint32_t p1 = (get_unbiased_exponent (x) > 0 && !is_inf_or_nan (x)
> +    ? (MANTISSA_MASK + 1) : 0);
> +  uint32_t p2 = (x & MANTISSA_MASK);
> +  return p1 | p2;
> +}
> 
> I don't think we need this (and anything called by it).
> 
> +static inline uint32_t
> +set_mantissa (uint32_t x, uint32_t m)
> +{
> +  m &= MANTISSA_MASK;
> +  x &= ~(MANTISSA_MASK);
> +  return x |= m;
> +}
> +
> +static inline uint32_t
> +get_mantissa (uint32_t x)
> +{
> +  return x & MANTISSA_MASK;
> +}
> +
> +static inline uint32_t
> +set_unbiased_exponent (uint32_t x, uint32_t e)
> +{
> +  e = (e << MANTISSA_WIDTH) & EXPONENT_MASK;
> +  x &= ~(EXPONENT_MASK);
> +  return x |= e;
> +}
> +
> +/* Convert integer number X, unbiased exponent EP, and sign S to double:
> +
> +   result = X * 2^(EP+1 - exponent_bias)
> +
> +   NB: zero is not supported.  */
> +static inline double
> +make_float (uint32_t x, int ep, uint32_t s)
> +{
> +  int lz = __builtin_clz (x) - EXPONENT_WIDTH;
> +  x <<= lz;
> +  ep -= lz;
> +
> +  uint32_t r = 0;
> +  if (__glibc_likely (ep >= 0))
> +    {
> +      r = set_mantissa (r, x);
> +      r = set_unbiased_exponent (r, ep + 1);
> +    }
> +  else
> +    r = set_mantissa (r, x >> -ep);
> +
> +  return asfloat (r | s);
> +}
> 
> This is overly complex, make_float is as trivial as:
> 
> static inline double
> make_float (uint32_t x, int ep, uint32_t s)
> {
>   int lz = __builtin_clz (x) - EXPONENT_WIDTH;
>   x <<= lz;
>   ep -= lz;
> 
>   if (__glibc_unlikely (ep < 0))
>   {
>     x >> -ep;
>     ep = 0;
>   }
>   return asfloat (s + x + (ep << MANTISSA_WIDTH));
> }
> 
> And now you don't need set_unbiased_exponent and set_mantissa either.

Indeed it is simpler, thanks.

> 
> Cheers,
> Wilco