[PATCH 4/4] math: Improve fmodf

[Wilco Dijkstra]

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.

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);

+    }
+
+  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).

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

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

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

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

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

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

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

Cheers,
Wilco