[PATCH v2 4/5] math: Improve fmodf

[Wilco Dijkstra]

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

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

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

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