From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: from mail-yb1-xb32.google.com (mail-yb1-xb32.google.com [IPv6:2607:f8b0:4864:20::b32]) by sourceware.org (Postfix) with ESMTPS id 659A8385840E for ; Thu, 16 Mar 2023 00:59:07 +0000 (GMT) DMARC-Filter: OpenDMARC Filter v1.4.2 sourceware.org 659A8385840E Authentication-Results: sourceware.org; dmarc=pass (p=none dis=none) header.from=gmail.com Authentication-Results: sourceware.org; spf=pass smtp.mailfrom=gmail.com Received: by mail-yb1-xb32.google.com with SMTP id p203so127787ybb.13 for ; Wed, 15 Mar 2023 17:59:07 -0700 (PDT) DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=20210112; t=1678928347; h=content-transfer-encoding:cc:to:subject:message-id:date:from :in-reply-to:references:mime-version:from:to:cc:subject:date :message-id:reply-to; bh=PeWky+P87CBqEQCdOUYA0yNJkOnzDGR0+BuLw8trdb0=; b=lffXywxgtDbxcWVhJtjtNlEXeSzOErDTcRMUmqcKMXe8v5TGdzVfaYqpaaLvFz50m2 tcPkLx96MPb0gWkVMguiydK0YM8sxTbK0l2ABnfubzsR9hLw2MtV3vi1qmIt1y0B+nQ7 No2FRANCXxzW36UGFOD9nQNbAARsk6oudWsOGSP0cYdp16JPgaqxAirrlexUxzWkPxfv Zvw/EnFGtJ45+08binvPil7I8v/gmx4nhz/0v7QIS15e+bgBRjcvbsrA6IksnUdKuy8U DrW4OS834/Y0JSBW+ztNEQELKRX5q8akC77PR3mvH922uzkARB1WPEdQjQ/ES8K6iwbm RNUw== X-Google-DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=1e100.net; s=20210112; t=1678928347; h=content-transfer-encoding:cc:to:subject:message-id:date:from :in-reply-to:references:mime-version:x-gm-message-state:from:to:cc :subject:date:message-id:reply-to; bh=PeWky+P87CBqEQCdOUYA0yNJkOnzDGR0+BuLw8trdb0=; b=P8l7bUsRTdCAf+u8a8XNAt6ibxE6fQmpBfwisXCOgl5AFcIc02+Pe6yp+fL6vBhIV1 aDkKZkLvthBm1y4qBhT5jKKdrjXMpPJwB+pyjmJljJrqpjx0SjJ9xns6ycKOQfYNn6je 2uJe7Ymf2P5BA+5g+M9felWwSjTYb6JcCUp9MA0vHwlRINe1fF8sL6EpiOcG/EymBhiB G+uAejxsa309EgmTXv4PM9+C1SH5MdpyV23RvHqrzYpV9dnn/hhH0CLR2044VvI8cpgr WeKOouYTl748riSdK87aJ5toWFrv04A0Nj97ExaRSrKqZvwSZucrDQ9A3IEbtFUfFiOF g82Q== X-Gm-Message-State: AO0yUKVX2GH3eHp0kcpgHtJ8oIAUQVYMVE/jJ12zfAfCea1RsAdXjj5P YYlfvv2Ye9ZqRINRON8r/BnLSVyY8JtORm0/Ax0= X-Google-Smtp-Source: AK7set9KYcCpgBgvRIHAT/xHxfnJcVrz0bg5lj73dVv1KV99PMV2HzcI32XhqlIcTQajaTwA6MbtijNrpnRSOsp8Gds= X-Received: by 2002:a5b:24c:0:b0:b3b:a48c:b23e with SMTP id g12-20020a5b024c000000b00b3ba48cb23emr8751427ybp.2.1678928346593; Wed, 15 Mar 2023 17:59:06 -0700 (PDT) MIME-Version: 1.0 References: <20230315205910.4120377-1-adhemerval.zanella@linaro.org> <20230315205910.4120377-4-adhemerval.zanella@linaro.org> In-Reply-To: <20230315205910.4120377-4-adhemerval.zanella@linaro.org> From: "H.J. Lu" Date: Wed, 15 Mar 2023 17:58:30 -0700 Message-ID: Subject: Re: [PATCH v2 3/5] math: Improve fmod To: Adhemerval Zanella Cc: libc-alpha@sourceware.org, Wilco Dijkstra , kirill Content-Type: text/plain; charset="UTF-8" Content-Transfer-Encoding: quoted-printable X-Spam-Status: No, score=-3021.8 required=5.0 tests=BAYES_00,DKIM_SIGNED,DKIM_VALID,DKIM_VALID_AU,DKIM_VALID_EF,FREEMAIL_FROM,GIT_PATCH_0,KAM_ASCII_DIVIDERS,KAM_SHORT,RCVD_IN_DNSWL_NONE,SPF_HELO_NONE,SPF_PASS,TXREP autolearn=ham autolearn_force=no version=3.4.6 X-Spam-Checker-Version: SpamAssassin 3.4.6 (2021-04-09) on server2.sourceware.org List-Id: On Wed, Mar 15, 2023 at 1:59=E2=80=AFPM Adhemerval Zanella wrote: > > This uses a new algorithm similar to already proposed earlier [1]. > With x =3D mx * 2^ex and y =3D my * 2^ey (mx, my, ex, ey being integers), > the simplest implementation is: > > mx * 2^ex =3D=3D 2 * mx * 2^(ex - 1) > > while (ex > ey) > { > mx *=3D 2; > --ex; > mx %=3D 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 -=3D 11; > mx %=3D 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 I tried it with the test in https://sourceware.org/bugzilla/show_bug.cgi?id=3D30179 On Intel i7-10710U, I got time ./sse 3.13user 0.00system 0:03.13elapsed 99%CPU (0avgtext+0avgdata 512maxresident= )k 0inputs+0outputs (0major+37minor)pagefaults 0swaps time ./x87 0.24user 0.00system 0:00.24elapsed 100%CPU (0avgtext+0avgdata 512maxresiden= t)k 0inputs+0outputs (0major+37minor)pagefaults 0swaps time ./generic 0.55user 0.00system 0:00.55elapsed 99%CPU (0avgtext+0avgdata 512maxresident= )k 0inputs+0outputs (0major+37minor)pagefaults 0swaps The new generic is still slower than x87. > 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 > > [1] https://sourceware.org/pipermail/libc-alpha/2020-November/119794.html > --- > math/libm-test-fmod.inc | 6 + > sysdeps/ieee754/dbl-64/e_fmod.c | 232 ++++++++++++++++----------- > sysdeps/ieee754/dbl-64/math_config.h | 61 +++++++ > 3 files changed, 203 insertions(+), 96 deletions(-) > > diff --git a/math/libm-test-fmod.inc b/math/libm-test-fmod.inc > index 39fd02f9ef..76c1e0cd45 100644 > --- a/math/libm-test-fmod.inc > +++ b/math/libm-test-fmod.inc > @@ -217,6 +217,12 @@ static const struct test_ff_f_data fmod_test_data[] = =3D > TEST_ff_f (fmod, -0x1p1023L, 0x3p-1022L, -0x1p-1021L, NO_INEXACT_EXC= EPTION|ERRNO_UNCHANGED), > TEST_ff_f (fmod, -0x1p1023L, -0x3p-1022L, -0x1p-1021L, NO_INEXACT_EX= CEPTION|ERRNO_UNCHANGED), > #endif > +#if TEST_COND_binary64 > + TEST_ff_f (fmod, 0x0.cded0a47373e9p-1022, 0x0.3212f5b8c8c16p-1022, 0= x0.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 <=3D -16381 > TEST_ff_f (fmod, 0x1p16383L, 0x3p-16445L, 0x1p-16445L, NO_INEXACT_EX= CEPTION|ERRNO_UNCHANGED), > TEST_ff_f (fmod, 0x1p16383L, -0x3p-16445L, 0x1p-16445L, NO_INEXACT_E= XCEPTION|ERRNO_UNCHANGED), > diff --git a/sysdeps/ieee754/dbl-64/e_fmod.c b/sysdeps/ieee754/dbl-64/e_f= mod.c > index 60b8bbb9d2..d21143e0cf 100644 > --- a/sysdeps/ieee754/dbl-64/e_fmod.c > +++ b/sysdeps/ieee754/dbl-64/e_fmod.c > @@ -1,105 +1,145 @@ > -/* > - * =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D > - * 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. > - * =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D > - */ > - > -/* > - * __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 > + . */ > > -#include > -#include > -#include > #include > +#include > +#include "math_config.h" > + > +/* With x =3D mx * 2^ex and y =3D my * 2^ey (mx, my, ex, ey being intege= rs), the > + simplest implementation is: > + > + mx * 2^ex =3D=3D 2 * mx * 2^(ex - 1) > > -static const double one =3D 1.0, Zero[] =3D {0.0, -0.0,}; > + while (ex > ey) > + { > + mx *=3D 2; > + --ex; > + mx %=3D my; > + } > + > + With mx/my being mantissa of double floating pointer, on each step th= e > + argument reduction can be improved 11 (which is sizeo of uint64_t min= us > + MANTISSA_WIDTH plus the signal bit): > + > + while (ex > ey) > + { > + mx << 11; > + ex -=3D 11; > + mx %=3D 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 =3D hx&UINT64_C(0x8000000000000000); /* sign of x */ > - hx ^=3Dsx; /* |x| */ > - hy &=3D UINT64_C(0x7fffffffffffffff); /* |y| */ > - > - /* purge off exception values */ > - if(__builtin_expect(hy=3D=3D0 > - || hx >=3D UINT64_C(0x7ff0000000000000) > - || hy > UINT64_C(0x7ff0000000000000), 0)) > - /* y=3D0,or x not finite or y is NaN */ > - return (x*y)/(x*y); > - if(__builtin_expect(hx<=3Dhy, 0)) { > - if(hx - return Zero[(uint64_t)sx>>63]; /* |x|=3D|y| return x*0*/ > - } > - > - /* determine ix =3D ilogb(x) */ > - if(__builtin_expect(hx - /* subnormal x */ > - for (ix =3D -1022,i=3D(hx<<11); i>0; i<<=3D1) ix -=3D1; > - } else ix =3D (hx>>52)-1023; > - > - /* determine iy =3D ilogb(y) */ > - if(__builtin_expect(hy - for (iy =3D -1022,i=3D(hy<<11); i>0; i<<=3D1) iy -=3D1; > - } else iy =3D (hy>>52)-1023; > - > - /* set up hx, hy and align y to x */ > - if(__builtin_expect(ix >=3D -1022, 1)) > - hx =3D UINT64_C(0x0010000000000000)|(UINT64_C(0x000ffffffffff= fff)&hx); > - else { /* subnormal x, shift x to normal */ > - n =3D -1022-ix; > - hx<<=3Dn; > - } > - if(__builtin_expect(iy >=3D -1022, 1)) > - hy =3D UINT64_C(0x0010000000000000)|(UINT64_C(0x000ffffffffff= fff)&hy); > - else { /* subnormal y, shift y to normal */ > - n =3D -1022-iy; > - hy<<=3Dn; > - } > - > - /* fix point fmod */ > - n =3D ix - iy; > - while(n--) { > - hz=3Dhx-hy; > - if(hz<0){hx =3D hx+hx;} > - else { > - if(hz=3D=3D0) /* return sign(x)*0 */ > - return Zero[(uint64_t)sx>>63]; > - hx =3D hz+hz; > - } > - } > - hz=3Dhx-hy; > - if(hz>=3D0) {hx=3Dhz;} > - > - /* convert back to floating value and restore the sign */ > - if(hx=3D=3D0) /* return sign(x)*0 */ > - return Zero[(uint64_t)sx>>63]; > - while(hx - hx =3D hx+hx; > - iy -=3D 1; > - } > - if(__builtin_expect(iy>=3D -1022, 1)) { /* normalize output */ > - hx =3D ((hx-UINT64_C(0x0010000000000000))|((uint64_t)(iy+1023)<= <52)); > - INSERT_WORDS64(x,hx|sx); > - } else { /* subnormal output */ > - n =3D -1022 - iy; > - hx>>=3Dn; > - INSERT_WORDS64(x,hx|sx); > - x *=3D one; /* create necessary signal */ > - } > - return x; /* exact output */ > + uint64_t hx =3D asuint64 (x); > + uint64_t hy =3D asuint64 (y); > + > + uint64_t sx =3D hx & SIGN_MASK; > + /* Get |x| and |y|. */ > + hx ^=3D sx; > + hy &=3D ~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 =3D=3D 0 || hx >=3D EXPONENT_MASK || h= y > EXPONENT_MASK)) > + return (x * y) / (x * y); > + > + if (__glibc_unlikely (hx <=3D hy)) > + { > + if (hx < hy) > + return x; > + return asdouble (sx); > + } > + > + int ex =3D hx >> MANTISSA_WIDTH; > + int ey =3D hy >> MANTISSA_WIDTH; > + > + /* Common case where exponents are close: ey >=3D -907 and |x/y| < 2^5= 2, */ > + if (__glibc_likely (ey > MANTISSA_WIDTH && ex - ey <=3D EXPONENT_WIDTH= )) > + { > + uint64_t mx =3D (hx & MANTISSA_MASK) | (MANTISSA_MASK + 1); > + uint64_t my =3D (hy & MANTISSA_MASK) | (MANTISSA_MASK + 1); > + > + uint64_t d =3D (ex =3D=3D 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 =3D=3D 0 && ey =3D=3D 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 =3D get_mantissa (hx) | (MANTISSA_MASK + 1); > + ex--; > + uint64_t my =3D get_mantissa (hy) | (MANTISSA_MASK + 1); > + > + int lead_zeros_my =3D EXPONENT_WIDTH; > + if (__glibc_likely (ey > 0)) > + ey--; > + else > + { > + my =3D hy; > + lead_zeros_my =3D clz_uint64 (my); > + } > + > + /* Assume hy !=3D 0 */ > + int tail_zeros_my =3D ctz_uint64 (my); > + int sides_zeroes =3D lead_zeros_my + tail_zeros_my; > + int exp_diff =3D ex - ey; > + > + int right_shift =3D exp_diff < tail_zeros_my ? exp_diff : tail_zeros_m= y; > + my >>=3D right_shift; > + exp_diff -=3D right_shift; > + ey +=3D right_shift; > + > + int left_shift =3D exp_diff < EXPONENT_WIDTH ? exp_diff : EXPONENT_WID= TH; > + mx <<=3D left_shift; > + exp_diff -=3D left_shift; > + > + mx %=3D my; > + > + if (__glibc_unlikely (mx =3D=3D 0)) > + return asdouble (sx); > + > + if (exp_diff =3D=3D 0) > + return make_double (my, ey, sx); > + > + /* Assume modulo/divide operation is slow, so use multiplication with = invert > + values. */ > + uint64_t inv_hy =3D UINT64_MAX / my; > + while (exp_diff > sides_zeroes) { > + exp_diff -=3D sides_zeroes; > + uint64_t hd =3D (mx * inv_hy) >> (BIT_WIDTH - sides_zeroes); > + mx <<=3D sides_zeroes; > + mx -=3D hd * my; > + while (__glibc_unlikely (mx > my)) > + mx -=3D my; > + } > + uint64_t hd =3D (mx * inv_hy) >> (BIT_WIDTH - exp_diff); > + mx <<=3D exp_diff; > + mx -=3D hd * my; > + while (__glibc_unlikely (mx > my)) > + mx -=3D 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-6= 4/math_config.h > index 3cbaeede64..d00f629531 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) =3D=3D sizeof (unsigned long)) > + return __builtin_clzl (x); > + else > + return __builtin_clzll (x); > +} > + > +static inline int > +ctz_uint64 (uint64_t x) > +{ > + if (sizeof (uint64_t) =3D=3D 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 =3D 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 =3D clz_uint64 (x) - EXPONENT_WIDTH; > + x <<=3D lz; > + ep -=3D lz; > + > + if (__glibc_unlikely (ep < 0)) > + { > + x >>=3D -ep; > + ep =3D 0; > + } > + > + return asdouble (s + x + (ep << MANTISSA_WIDTH)); > +} > + > + > #define NOINLINE __attribute__ ((noinline)) > > /* Error handling tail calls for special cases, with a sign argument. > -- > 2.34.1 > --=20 H.J.