[PATCH v2 4/5] math: Improve fmodf

[Adhemerval Zanella Netto]

On 16/03/23 15:11, Wilco Dijkstra wrote:
> Hi Adhemerval,
> 
> 
> +  /* Common case where exponents are close: ey >= -103 and |x/y| < 2^8,  */
> +  if (__glibc_likely (ey > MANTISSA_WIDTH && ex - ey <= EXPONENT_WIDTH))
> +    {
> +      uint64_t mx = (hx & MANTISSA_MASK) | (MANTISSA_MASK + 1);
> +      uint64_t my = (hy & MANTISSA_MASK) | (MANTISSA_MASK + 1);
> +
> +      uint32_t d = (ex == ey) ? (mx - my) : (mx << (ex - ey)) % my;
> +      return make_float (d, ey - 1, sx);
> +    }
> 
> So we do need to handle zero case here, I think fmodf (5, 1) will now
> return 0.5... We could add support in make_float (see below).

Ok, although at least for fmodf (5, 1) current algorithm does return 0.

> 
> +  /* Special case, both x and y are subnormal.  */
> +  if (__glibc_unlikely (ex == 0 && ey == 0))
> +    return asfloat (sx | hx % hy);
> 
> This will handle zero case correctly.
> 
> +  /* 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_mantissa (hx) | (MANTISSA_MASK + 1);
> +  ex--;
> +
> +  uint32_t my = get_mantissa (hy) | (MANTISSA_MASK + 1);
> 
> Odd uses of get_mantissa here but not in the fast path above... They are computing
> the same value so can be shared.
> 
> +  mx %= my;
> +
> +  if (__glibc_unlikely (mx == 0))
> +    return asfloat (sx);
> +
> +  if (exp_diff == 0)
> +    return make_float (my, ey, sx);
> 
> That should be make_float (mx, ey, sx) - no testcase caught this?!?
> 
> (note both issues are in double version too)

Unfortunately no, I will add proper tests for this.

> 
> --- a/sysdeps/ieee754/flt-32/math_config.h
> +++ b/sysdeps/ieee754/flt-32/math_config.h
> 
> +   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;
> +
> +  if (__glibc_unlikely (ep < 0))
> 
> We could do (ep < 0 || x == 0) to handle the zero case correctly (this works
> eventhough clz (0) is not well defined).

Ack.

> 
> +    {
> +      x >>= -ep;
> +      ep = 0;
> +    }
> +  return asfloat (s + x + (ep << MANTISSA_WIDTH));
> +}
> +
>  
> Cheers,
> Wilco