@@ -216,6 +216,11 @@ static const struct test_ff_f_data fmod_test_data[] =
TEST_ff_f (fmod, -0x1p1023L, -0x3p-1073L, -0x1p-1073L, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
TEST_ff_f (fmod, -0x1p1023L, 0x3p-1022L, -0x1p-1021L, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
TEST_ff_f (fmod, -0x1p1023L, -0x3p-1022L, -0x1p-1021L, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
+#else
+ 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),
#endif
#if TEST_COND_binary64
TEST_ff_f (fmod, 0x0.cded0a47373e9p-1022, 0x0.3212f5b8c8c16p-1022, 0x0.05a1336414391p-1022, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
@@ -1,102 +1,144 @@
-/* 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)
-static const float one = 1.0, Zero[] = {0.0, -0.0,};
+ 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;
+ } */
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 (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)
@@ -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 >>= -ep;
+ ep = 0;
+ }
+ return asfloat (s + x + (ep << MANTISSA_WIDTH));
+}
+
#define NOINLINE __attribute__ ((noinline))
attribute_hidden float __math_oflowf (uint32_t);
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 --- math/libm-test-fmod.inc | 5 + sysdeps/ieee754/flt-32/e_fmodf.c | 228 ++++++++++++++++----------- sysdeps/ieee754/flt-32/math_config.h | 41 +++++ 3 files changed, 181 insertions(+), 93 deletions(-)