mirror of
https://sourceware.org/git/glibc.git
synced 2024-11-21 12:30:06 +00:00
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>
This commit is contained in:
Adhemerval Zanella Netto
Adhemerval Zanella
parent
5c11701c51
commit
34b9f8bc17
@ -152,6 +152,16 @@ static const struct test_ff_f_data fmod_test_data[] =
|
||||
TEST_ff_f (fmod, 6.5, -2.25L, 2.0L, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
|
||||
TEST_ff_f (fmod, -6.5, -2.25L, -2.0L, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
|
||||
|
||||
TEST_ff_f (fmod, 0x1.e848p+18, 0x1.6p+3, 0x1.8p+2, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
|
||||
TEST_ff_f (fmod, 0x1.e848p+18, -0x1.6p+3, 0x1.8p+2, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
|
||||
TEST_ff_f (fmod, -0x1.e848p+18, 0x1.6p+3, -0x1.8p+2, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
|
||||
TEST_ff_f (fmod, -0x1.e848p+18, -0x1.6p+3, -0x1.8p+2, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
|
||||
|
||||
TEST_ff_f (fmod, 0x1.4p+2, 0x1p+0, 0x0p+0, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
|
||||
TEST_ff_f (fmod, 0x1.4p+2, -0x1p+0, 0x0p+0, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
|
||||
TEST_ff_f (fmod, -0x1.4p+2, 0x1p+0, -0x0p+0, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
|
||||
TEST_ff_f (fmod, -0x1.4p+2, -0x1p+0, -0x0p+0, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
|
||||
|
||||
TEST_ff_f (fmod, max_value, max_value, 0, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
|
||||
TEST_ff_f (fmod, max_value, -max_value, 0, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
|
||||
TEST_ff_f (fmod, max_value, min_value, 0, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
|
||||
@ -216,6 +226,10 @@ 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),
|
||||
TEST_ff_f (fmod, 0x0.cded0a47373e9p-1022, 0x0.3212f5b8c8c16p-1022, 0x0.05a1336414391p-1022, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
|
||||
TEST_ff_f (fmod, 0x0.cded0a47373e9p-1022, -0x0.3212f5b8c8c16p-1022, 0x0.05a1336414391p-1022, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
|
||||
TEST_ff_f (fmod, -0x0.cded0a47373e9p-1022, 0x0.3212f5b8c8c16p-1022, -0x0.05a1336414391p-1022, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
|
||||
TEST_ff_f (fmod, -0x0.cded0a47373e9p-1022, -0x0.3212f5b8c8c16p-1022, -0x0.05a1336414391p-1022, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
|
||||
#endif
|
||||
#if MIN_EXP <= -16381
|
||||
TEST_ff_f (fmod, 0x1p16383L, 0x3p-16445L, 0x1p-16445L, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
|
||||
|
||||
@ -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.
|
||||
* ====================================================
|
||||
*/
|
||||
/* Floating-point remainder function.
|
||||
Copyright (C) 2023 Free Software Foundation, Inc.
|
||||
This file is part of the GNU C Library.
|
||||
|
||||
/*
|
||||
* __ieee754_fmod(x,y)
|
||||
* Return x mod y in exact arithmetic
|
||||
* Method: shift and subtract
|
||||
*/
|
||||
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"
|
||||
|
||||
static const double one = 1.0, Zero[] = {0.0, -0.0,};
|
||||
/* 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
|
||||
|
||||
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;
|
||||
uint64_t hx = asuint64 (x);
|
||||
uint64_t hy = asuint64 (y);
|
||||
|
||||
EXTRACT_WORDS64(hx,x);
|
||||
EXTRACT_WORDS64(hy,y);
|
||||
sx = hx&UINT64_C(0x8000000000000000); /* sign of x */
|
||||
hx ^=sx; /* |x| */
|
||||
hy &= UINT64_C(0x7fffffffffffffff); /* |y| */
|
||||
uint64_t sx = hx & SIGN_MASK;
|
||||
/* Get |x| and |y|. */
|
||||
hx ^= sx;
|
||||
hy &= ~SIGN_MASK;
|
||||
|
||||
/* 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*/
|
||||
}
|
||||
/* 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);
|
||||
|
||||
/* 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;
|
||||
if (__glibc_unlikely (hx <= hy))
|
||||
{
|
||||
if (hx < hy)
|
||||
return x;
|
||||
return asdouble (sx);
|
||||
}
|
||||
|
||||
/* 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;
|
||||
int ex = hx >> MANTISSA_WIDTH;
|
||||
int ey = hy >> MANTISSA_WIDTH;
|
||||
|
||||
/* 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;
|
||||
}
|
||||
/* 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);
|
||||
|
||||
/* 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;}
|
||||
uint64_t d = (ex == ey) ? (mx - my) : (mx << (ex - ey)) % my;
|
||||
return make_double (d, ey - 1, sx);
|
||||
}
|
||||
|
||||
/* 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 */
|
||||
/* 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)
|
||||
|
||||
@ -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.
|
||||
|
||||
Loading…
Reference in New Issue
Block a user