math: Improve fmod

This uses a new algorithm similar to already proposed earlier [1].
With x = mx * 2^ex and y = my * 2^ey (mx, my, ex, ey being integers),
the simplest implementation is:

   mx * 2^ex == 2 * mx * 2^(ex - 1)

   while (ex > ey)
     {
       mx *= 2;
       --ex;
       mx %= my;
     }

With mx/my being mantissa of double floating pointer, on each step the
argument reduction can be improved 11 (which is sizeo of uint64_t minus
MANTISSA_WIDTH plus the signal bit):

   while (ex > ey)
     {
       mx << 11;
       ex -= 11;
       mx %= my;
     }  */

The implementation uses builtin clz and ctz, along with shifts to
convert hx/hy back to doubles.  Different than the original patch,
this path assume modulo/divide operation is slow, so use multiplication
with invert values.

I see the following performance improvements using fmod benchtests
(result only show the 'mean' result):

  Architecture     | Input           | master   | patch
  -----------------|-----------------|----------|--------
  x86_64 (Ryzen 9) | subnormals      | 19.1584  | 12.5049
  x86_64 (Ryzen 9) | normal          | 1016.51  | 296.939
  x86_64 (Ryzen 9) | close-exponents | 18.4428  | 16.0244
  aarch64 (N1)     | subnormal       | 11.153   | 6.81778
  aarch64 (N1)     | normal          | 528.649  | 155.62
  aarch64 (N1)     | close-exponents | 11.4517  | 8.21306

I also see similar improvements on arm-linux-gnueabihf when running on
the N1 aarch64 chips, where it a lot of soft-fp implementation (for
modulo, clz, ctz, and multiplication):

  Architecture     | Input           | master   | patch
  -----------------|-----------------|----------|--------
  armhf (N1)       | subnormal       | 15.908   | 15.1083
  armhf (N1)       | normal          | 837.525  | 244.833
  armhf (N1)       | close-exponents | 16.2111  | 21.8182

Instead of using the math_private.h definitions, I used the
math_config.h instead which is used on newer math implementations.

Co-authored-by: kirill <kirill.okhotnikov@gmail.com>

[1] https://sourceware.org/pipermail/libc-alpha/2020-November/119794.html
Reviewed-by: Wilco Dijkstra  <Wilco.Dijkstra@arm.com>

2 files changed, 208 insertions, 96 deletions

diff --git a/sysdeps/ieee754/dbl-64/e_fmod.c b/sysdeps/ieee754/dbl-64/e_fmod.c
index 60b8bbb9d2..e661ca1ff8 100644
--- a/sysdeps/ieee754/dbl-64/e_fmod.c
+++ b/sysdeps/ieee754/dbl-64/e_fmod.c
@@ -1,105 +1,156 @@
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
-
-/*
- * __ieee754_fmod(x,y)
- * Return x mod y in exact arithmetic
- * Method: shift and subtract
- */
+/* Floating-point remainder function.
+   Copyright (C) 2023 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, see
+   <https://www.gnu.org/licenses/>.  */
 
-#include <math.h>
-#include <math_private.h>
-#include <stdint.h>
 #include <libm-alias-finite.h>
+#include <math.h>
+#include "math_config.h"
+
+/* With x = mx * 2^ex and y = my * 2^ey (mx, my, ex, ey being integers), the
+   simplest implementation is:
+
+   mx * 2^ex == 2 * mx * 2^(ex - 1)
+
+   or
 
-static const double one = 1.0, Zero[] = {0.0, -0.0,};
+   while (ex > ey)
+     {
+       mx *= 2;
+       --ex;
+       mx %= my;
+     }
+
+   With the mathematical equivalence of:
+
+   r == x % y == (x % (N * y)) % y
+
+   And with mx/my being mantissa of double floating point number (which uses
+   less bits than the storage type), on each step the argument reduction can
+   be improved by 11 (which is the size of uint64_t minus MANTISSA_WIDTH plus
+   the signal bit):
+
+   mx * 2^ex == 2^11 * mx * 2^(ex - 11)
+
+   or
+
+   while (ex > ey)
+     {
+       mx << 11;
+       ex -= 11;
+       mx %= my;
+     }  */
 
 double
 __ieee754_fmod (double x, double y)
 {
-	int32_t n,ix,iy;
-	int64_t hx,hy,hz,sx,i;
-
-	EXTRACT_WORDS64(hx,x);
-	EXTRACT_WORDS64(hy,y);
-	sx = hx&UINT64_C(0x8000000000000000);	/* sign of x */
-	hx ^=sx;				/* |x| */
-	hy &= UINT64_C(0x7fffffffffffffff);	/* |y| */
-
-    /* purge off exception values */
-	if(__builtin_expect(hy==0
-			    || hx >= UINT64_C(0x7ff0000000000000)
-			    || hy > UINT64_C(0x7ff0000000000000), 0))
-	  /* y=0,or x not finite or y is NaN */
-	    return (x*y)/(x*y);
-	if(__builtin_expect(hx<=hy, 0)) {
-	    if(hx<hy) return x;	/* |x|<|y| return x */
-	    return Zero[(uint64_t)sx>>63];	/* |x|=|y| return x*0*/
-	}
-
-    /* determine ix = ilogb(x) */
-	if(__builtin_expect(hx<UINT64_C(0x0010000000000000), 0)) {
-	  /* subnormal x */
-	  for (ix = -1022,i=(hx<<11); i>0; i<<=1) ix -=1;
-	} else ix = (hx>>52)-1023;
-
-    /* determine iy = ilogb(y) */
-	if(__builtin_expect(hy<UINT64_C(0x0010000000000000), 0)) {	/* subnormal y */
-	  for (iy = -1022,i=(hy<<11); i>0; i<<=1) iy -=1;
-	} else iy = (hy>>52)-1023;
-
-    /* set up hx, hy and align y to x */
-	if(__builtin_expect(ix >= -1022, 1))
-	    hx = UINT64_C(0x0010000000000000)|(UINT64_C(0x000fffffffffffff)&hx);
-	else {		/* subnormal x, shift x to normal */
-	    n = -1022-ix;
-	    hx<<=n;
-	}
-	if(__builtin_expect(iy >= -1022, 1))
-	    hy = UINT64_C(0x0010000000000000)|(UINT64_C(0x000fffffffffffff)&hy);
-	else {		/* subnormal y, shift y to normal */
-	    n = -1022-iy;
-	    hy<<=n;
-	}
-
-    /* fix point fmod */
-	n = ix - iy;
-	while(n--) {
-	    hz=hx-hy;
-	    if(hz<0){hx = hx+hx;}
-	    else {
-		if(hz==0)		/* return sign(x)*0 */
-		    return Zero[(uint64_t)sx>>63];
-		hx = hz+hz;
-	    }
-	}
-	hz=hx-hy;
-	if(hz>=0) {hx=hz;}
-
-    /* convert back to floating value and restore the sign */
-	if(hx==0)			/* return sign(x)*0 */
-	    return Zero[(uint64_t)sx>>63];
-	while(hx<UINT64_C(0x0010000000000000)) {	/* normalize x */
-	    hx = hx+hx;
-	    iy -= 1;
-	}
-	if(__builtin_expect(iy>= -1022, 1)) {	/* normalize output */
-	  hx = ((hx-UINT64_C(0x0010000000000000))|((uint64_t)(iy+1023)<<52));
-	    INSERT_WORDS64(x,hx|sx);
-	} else {		/* subnormal output */
-	    n = -1022 - iy;
-	    hx>>=n;
-	    INSERT_WORDS64(x,hx|sx);
-	    x *= one;		/* create necessary signal */
-	}
-	return x;		/* exact output */
+  uint64_t hx = asuint64 (x);
+  uint64_t hy = asuint64 (y);
+
+  uint64_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 asdouble (sx);
+    }
+
+  int ex = hx >> MANTISSA_WIDTH;
+  int ey = hy >> MANTISSA_WIDTH;
+
+  /* Common case where exponents are close: ey >= -907 and |x/y| < 2^52,  */
+  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);
+
+      uint64_t d = (ex == ey) ? (mx - my) : (mx << (ex - ey)) % my;
+      return make_double (d, ey - 1, sx);
+    }
+
+  /* Special case, both x and y are subnormal.  */
+  if (__glibc_unlikely (ex == 0 && ey == 0))
+    return asdouble (sx | hx % hy);
+
+  /* Convert |x| and |y| to 'mx + 2^ex' and 'my + 2^ey'.  Assume that hx is
+     not subnormal by conditions above.  */
+  uint64_t mx = get_mantissa (hx) | (MANTISSA_MASK + 1);
+  ex--;
+  uint64_t my = get_mantissa (hy) | (MANTISSA_MASK + 1);
+
+  int lead_zeros_my = EXPONENT_WIDTH;
+  if (__glibc_likely (ey > 0))
+    ey--;
+  else
+    {
+      my = hy;
+      lead_zeros_my = clz_uint64 (my);
+    }
+
+  /* Assume hy != 0  */
+  int tail_zeros_my = ctz_uint64 (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 asdouble (sx);
+
+  if (exp_diff == 0)
+    return make_double (mx, ey, sx);
+
+  /* Assume modulo/divide operation is slow, so use multiplication with invert
+     values.  */
+  uint64_t inv_hy = UINT64_MAX / my;
+  while (exp_diff > sides_zeroes) {
+    exp_diff -= sides_zeroes;
+    uint64_t hd = (mx * inv_hy) >> (BIT_WIDTH - sides_zeroes);
+    mx <<= sides_zeroes;
+    mx -= hd * my;
+    while (__glibc_unlikely (mx > my))
+      mx -= my;
+  }
+  uint64_t hd = (mx * inv_hy) >> (BIT_WIDTH - exp_diff);
+  mx <<= exp_diff;
+  mx -= hd * my;
+  while (__glibc_unlikely (mx > my))
+    mx -= my;
+
+  return make_double (mx, ey, sx);
 }
 libm_alias_finite (__ieee754_fmod, __fmod)
diff --git a/sysdeps/ieee754/dbl-64/math_config.h b/sysdeps/ieee754/dbl-64/math_config.h
index 3cbaeede64..2049cea3f7 100644
--- a/sysdeps/ieee754/dbl-64/math_config.h
+++ b/sysdeps/ieee754/dbl-64/math_config.h
@@ -43,6 +43,24 @@
 # define TOINT_INTRINSICS 0
 #endif
 
+static inline int
+clz_uint64 (uint64_t x)
+{
+  if (sizeof (uint64_t) == sizeof (unsigned long))
+    return __builtin_clzl (x);
+  else
+    return __builtin_clzll (x);
+}
+
+static inline int
+ctz_uint64 (uint64_t x)
+{
+  if (sizeof (uint64_t) == sizeof (unsigned long))
+    return __builtin_ctzl (x);
+  else
+    return __builtin_ctzll (x);
+}
+
 #if TOINT_INTRINSICS
 /* Round x to nearest int in all rounding modes, ties have to be rounded
    consistently with converttoint so the results match.  If the result
@@ -88,6 +106,49 @@ issignaling_inline (double x)
   return 2 * (ix ^ 0x0008000000000000) > 2 * 0x7ff8000000000000ULL;
 }
 
+#define BIT_WIDTH       64
+#define MANTISSA_WIDTH  52
+#define EXPONENT_WIDTH  11
+#define MANTISSA_MASK   UINT64_C(0x000fffffffffffff)
+#define EXPONENT_MASK   UINT64_C(0x7ff0000000000000)
+#define EXP_MANT_MASK   UINT64_C(0x7fffffffffffffff)
+#define QUIET_NAN_MASK  UINT64_C(0x0008000000000000)
+#define SIGN_MASK	UINT64_C(0x8000000000000000)
+
+static inline bool
+is_nan (uint64_t x)
+{
+  return (x & EXP_MANT_MASK) > EXPONENT_MASK;
+}
+
+static inline uint64_t
+get_mantissa (uint64_t x)
+{
+  return x & MANTISSA_MASK;
+}
+
+/* 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_double (uint64_t x, int64_t ep, uint64_t s)
+{
+  int lz = clz_uint64 (x) - EXPONENT_WIDTH;
+  x <<= lz;
+  ep -= lz;
+
+  if (__glibc_unlikely (ep < 0 || x == 0))
+    {
+      x >>= -ep;
+      ep = 0;
+    }
+
+  return asdouble (s + x + (ep << MANTISSA_WIDTH));
+}
+
+
 #define NOINLINE __attribute__ ((noinline))
 
 /* Error handling tail calls for special cases, with a sign argument.