math: Improve fmodf

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 8 (which is sizeof of uint32_t minus
MANTISSA_WIDTH plus the signal bit):

   while (ex > ey)
     {
       mx << 8;
       ex -= 8;
       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      | 17.2549  | 12.0318
  x86_64 (Ryzen 9) | normal          | 85.4096  | 49.9641
  x86_64 (Ryzen 9) | close-exponents | 19.1072  | 15.8224
  aarch64 (N1)     | subnormal       | 10.2182  | 6.81778
  aarch64 (N1)     | normal          | 60.0616  | 20.3667
  aarch64 (N1)     | close-exponents | 11.5256  | 8.39685

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, and multiplication):

  Architecture     | Input           | master   | patch
  -----------------|-----------------|----------|--------
  armhf (N1)       | subnormal       | 11.6662  | 10.8955
  armhf (N1)       | normal          | 69.2759  | 34.1524
  armhf (N1)       | close-exponents | 13.6472  | 18.2131

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>

3 files changed, 191 insertions, 93 deletions

diff --git a/math/libm-test-fmod.inc b/math/libm-test-fmod.inc
index 8841c13..43376c3 100644
--- a/math/libm-test-fmod.inc
+++ b/math/libm-test-fmod.inc
@@ -213,6 +213,10 @@ static const struct test_ff_f_data fmod_test_data[] =
     TEST_ff_f (fmod, -0x1p127L, -0x3p-148L, -0x1p-147L, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
     TEST_ff_f (fmod, -0x1p127L, 0x3p-126L, -0x1p-125L, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
     TEST_ff_f (fmod, -0x1p127L, -0x3p-126L, -0x1p-125L, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
+    TEST_ff_f (fmod, 0x1.3a3e6p-127, 0x1.8b8338p-128, 0x1.d1f31p-129, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
+    TEST_ff_f (fmod, 0x1.3a3e6p-127, -0x1.8b8338p-128, 0x1.d1f31p-129, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
+    TEST_ff_f (fmod, -0x1.3a3e6p-127, 0x1.8b8338p-128, -0x1.d1f31p-129, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
+    TEST_ff_f (fmod, -0x1.3a3e6p-127, -0x1.8b8338p-128, -0x1.d1f31p-129, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
 #if !TEST_COND_binary32
     TEST_ff_f (fmod, 0x1p1023L, 0x3p-1074L, 0x1p-1073L, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
     TEST_ff_f (fmod, 0x1p1023L, -0x3p-1074L, 0x1p-1073L, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
diff --git a/sysdeps/ieee754/flt-32/e_fmodf.c b/sysdeps/ieee754/flt-32/e_fmodf.c
index b71c4f7..4482b8c 100644
--- a/sysdeps/ieee754/flt-32/e_fmodf.c
+++ b/sysdeps/ieee754/flt-32/e_fmodf.c
@@ -1,102 +1,155 @@
-/* e_fmodf.c -- float version of e_fmod.c.
- */
-
-/*
- * ====================================================
- * 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_fmodf(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 <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 float 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 8 (which is the size of uint32_t minus MANTISSA_WIDTH plus
+   the signal bit):
+
+   mx * 2^ex == 2^8 * mx * 2^(ex - 8)
+
+   or
+
+   while (ex > ey)
+     {
+       mx << 8;
+       ex -= 8;
+       mx %= my;
+     }  */
 
 float
 __ieee754_fmodf (float x, float y)
 {
-	int32_t n,hx,hy,hz,ix,iy,sx,i;
-
-	GET_FLOAT_WORD(hx,x);
-	GET_FLOAT_WORD(hy,y);
-	sx = hx&0x80000000;		/* sign of x */
-	hx ^=sx;		/* |x| */
-	hy &= 0x7fffffff;	/* |y| */
-
-    /* purge off exception values */
-	if(hy==0||(hx>=0x7f800000)||		/* y=0,or x not finite */
-	   (hy>0x7f800000))			/* or y is NaN */
-	    return (x*y)/(x*y);
-	if(hx<hy) return x;			/* |x|<|y| return x */
-	if(hx==hy)
-	    return Zero[(uint32_t)sx>>31];	/* |x|=|y| return x*0*/
-
-    /* determine ix = ilogb(x) */
-	if(hx<0x00800000) {	/* subnormal x */
-	    for (ix = -126,i=(hx<<8); i>0; i<<=1) ix -=1;
-	} else ix = (hx>>23)-127;
-
-    /* determine iy = ilogb(y) */
-	if(hy<0x00800000) {	/* subnormal y */
-	    for (iy = -126,i=(hy<<8); i>=0; i<<=1) iy -=1;
-	} else iy = (hy>>23)-127;
-
-    /* set up {hx,lx}, {hy,ly} and align y to x */
-	if(ix >= -126)
-	    hx = 0x00800000|(0x007fffff&hx);
-	else {		/* subnormal x, shift x to normal */
-	    n = -126-ix;
-	    hx = hx<<n;
-	}
-	if(iy >= -126)
-	    hy = 0x00800000|(0x007fffff&hy);
-	else {		/* subnormal y, shift y to normal */
-	    n = -126-iy;
-	    hy = 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[(uint32_t)sx>>31];
-		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[(uint32_t)sx>>31];
-	while(hx<0x00800000) {		/* normalize x */
-	    hx = hx+hx;
-	    iy -= 1;
-	}
-	if(iy>= -126) {		/* normalize output */
-	    hx = ((hx-0x00800000)|((iy+127)<<23));
-	    SET_FLOAT_WORD(x,hx|sx);
-	} else {		/* subnormal output */
-	    n = -126 - iy;
-	    hx >>= n;
-	    SET_FLOAT_WORD(x,hx|sx);
-	    x *= one;		/* create necessary signal */
-	}
-	return x;		/* exact output */
+  uint32_t hx = asuint (x);
+  uint32_t hy = asuint (y);
+
+  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 asfloat (sx);
+    }
+
+  int ex = hx >> MANTISSA_WIDTH;
+  int ey = hy >> MANTISSA_WIDTH;
+
+  /* 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);
+    }
+
+  /* Special case, both x and y are subnormal.  */
+  if (__glibc_unlikely (ex == 0 && ey == 0))
+    return 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_mantissa (hx) | (MANTISSA_MASK + 1);
+  ex--;
+
+  uint32_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 = __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 asfloat (sx);
+
+  if (exp_diff == 0)
+    return make_float (mx, 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 23045f5..829430e 100644
--- a/sysdeps/ieee754/flt-32/math_config.h
+++ b/sysdeps/ieee754/flt-32/math_config.h
@@ -110,6 +110,47 @@ 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 uint32_t
+get_mantissa (uint32_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_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 == 0))
+    {
+      x >>= -ep;
+      ep = 0;
+    }
+  return asfloat (s + x + (ep << MANTISSA_WIDTH));
+}
+
 #define NOINLINE __attribute__ ((noinline))
 
 attribute_hidden float __math_oflowf (uint32_t);