From: Wilco Dijkstra <Wilco.Dijkstra@arm.com>
To: Adhemerval Zanella <adhemerval.zanella@linaro.org>,
"libc-alpha@sourceware.org" <libc-alpha@sourceware.org>,
"H . J . Lu" <hjl.tools@gmail.com>
Cc: kirill <kirill.okhotnikov@gmail.com>
Subject: Re: [PATCH 4/4] math: Improve fmodf
Date: Tue, 14 Mar 2023 16:42:37 +0000 [thread overview]
Message-ID: <PAWPR08MB89827C35016A4F09DD92BF3983BE9@PAWPR08MB8982.eurprd08.prod.outlook.com> (raw)
In-Reply-To: <20230310175900.2388957-5-adhemerval.zanella@linaro.org>
Hi Adhemerval,
This looks good overall. I guess we still needs error handling and the wrapper
disabled since it spends 1/3 of the time in the useless wrapper but that can be
a separate patch.
A few comments below, there are a few things that look wrong, and I think most
of the helper functions just add extra complexity and branches for no gain:
+ 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 sx ? -0.0 : 0.0;
This should be return asfloat (sx);
+ }
+
+ int ex = get_unbiased_exponent (hx);
+ int ey = get_unbiased_exponent (hy);
Should be hx >> MANTISSA_WIDTH since we cleared the sign bits (now we don't
need to add get_unbiased_exponent).
+ /* Common case where exponents are close: ey >= -103 and |x/y| < 2^8, */
+ if (__glibc_likely (ey > MANTISSA_WIDTH && ex - ey <= EXPONENT_WIDTH))
+ {
+ uint32_t mx = get_explicit_mantissa (hx);
+ uint32_t my = get_explicit_mantissa (hy);
Note this is equivalent to:
mx = (hx & MANTISSA_MASK) | (MANTISSA_MASK + 1);
So we don't need get_explicit_mantissa (or we could change it to do the above).
If we do this before the if statement, we don't need to repeat it below.
+ uint32_t d = (ex == ey) ? (mx - my) : (mx << (ex - ey)) % my;
+ if (d == 0)
+ return 0.0;
Looks like a bug, should be asfloat (sx)?
+ return make_float (d, ey - 1, sx);
+ }
+
+ /* Special case, both x and y are subnormal. */
+ if (__glibc_unlikely (ex == 0 && ey == 0))
+ return asfloat (hx % hy);
Similarly, shouldn't this be 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_explicit_mantissa (hx);
+ ex--;
+
+ uint32_t my = get_explicit_mantissa (hy);
If we set mx/my above then this isn't needed.
+ int lead_zeros_my = EXPONENT_WIDTH;
+ if (__glibc_likely (ey > 0))
+ ey--;
+ else
+ {
+ my = get_mantissa (hy);
This is really my = hy; since we know hy is positive denormal.
+ 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 sx ? -0.0 : 0.0;
Should be 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)
diff --git a/sysdeps/ieee754/flt-32/math_config.h b/sysdeps/ieee754/flt-32/math_config.h
index 23045f59d6..cdab3a36ef 100644
--- a/sysdeps/ieee754/flt-32/math_config.h
+++ b/sysdeps/ieee754/flt-32/math_config.h
@@ -110,6 +110,95 @@ 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 bool
+is_quiet_nan (uint32_t x)
+{
+ return (x & EXP_MANT_MASK) == (EXPONENT_MASK | QUIET_NAN_MASK);
+}
+
+static inline bool
+is_inf_or_nan (uint32_t x)
+{
+ return (x & EXPONENT_MASK) == EXPONENT_MASK;
+}
+
+static inline uint16_t
+get_unbiased_exponent (uint32_t x)
+{
+ return (x & EXPONENT_MASK) >> MANTISSA_WIDTH;
+}
+
+/* Return mantissa with the implicit bit set iff X is a normal number. */
+static inline uint32_t
+get_explicit_mantissa (uint32_t x)
+{
+ uint32_t p1 = (get_unbiased_exponent (x) > 0 && !is_inf_or_nan (x)
+ ? (MANTISSA_MASK + 1) : 0);
+ uint32_t p2 = (x & MANTISSA_MASK);
+ return p1 | p2;
+}
I don't think we need this (and anything called by it).
+static inline uint32_t
+set_mantissa (uint32_t x, uint32_t m)
+{
+ m &= MANTISSA_MASK;
+ x &= ~(MANTISSA_MASK);
+ return x |= m;
+}
+
+static inline uint32_t
+get_mantissa (uint32_t x)
+{
+ return x & MANTISSA_MASK;
+}
+
+static inline uint32_t
+set_unbiased_exponent (uint32_t x, uint32_t e)
+{
+ e = (e << MANTISSA_WIDTH) & EXPONENT_MASK;
+ x &= ~(EXPONENT_MASK);
+ return x |= e;
+}
+
+/* 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;
+
+ uint32_t r = 0;
+ if (__glibc_likely (ep >= 0))
+ {
+ r = set_mantissa (r, x);
+ r = set_unbiased_exponent (r, ep + 1);
+ }
+ else
+ r = set_mantissa (r, x >> -ep);
+
+ return asfloat (r | s);
+}
This is overly complex, make_float is as trivial as:
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));
}
And now you don't need set_unbiased_exponent and set_mantissa either.
Cheers,
Wilco
next prev parent reply other threads:[~2023-03-14 16:43 UTC|newest]
Thread overview: 10+ messages / expand[flat|nested] mbox.gz Atom feed top
2023-03-10 17:58 [PATCH 0/4] Improve fmod and fmodf Adhemerval Zanella
2023-03-10 17:58 ` [PATCH 1/4] benchtests: Add fmod benchmark Adhemerval Zanella
2023-03-10 17:58 ` [PATCH 2/4] benchtests: Add fmodf benchmark Adhemerval Zanella
2023-03-10 17:58 ` [PATCH 3/4] math: Improve fmod Adhemerval Zanella
2023-03-10 17:59 ` [PATCH 4/4] math: Improve fmodf Adhemerval Zanella
2023-03-10 23:17 ` H.J. Lu
2023-03-13 15:19 ` Matt Turner
2023-03-13 16:38 ` Adhemerval Zanella Netto
2023-03-14 16:42 ` Wilco Dijkstra [this message]
2023-03-15 17:50 ` Adhemerval Zanella Netto
Reply instructions:
You may reply publicly to this message via plain-text email
using any one of the following methods:
* Save the following mbox file, import it into your mail client,
and reply-to-all from there: mbox
Avoid top-posting and favor interleaved quoting:
https://en.wikipedia.org/wiki/Posting_style#Interleaved_style
* Reply using the --to, --cc, and --in-reply-to
switches of git-send-email(1):
git send-email \
--in-reply-to=PAWPR08MB89827C35016A4F09DD92BF3983BE9@PAWPR08MB8982.eurprd08.prod.outlook.com \
--to=wilco.dijkstra@arm.com \
--cc=adhemerval.zanella@linaro.org \
--cc=hjl.tools@gmail.com \
--cc=kirill.okhotnikov@gmail.com \
--cc=libc-alpha@sourceware.org \
/path/to/YOUR_REPLY
https://kernel.org/pub/software/scm/git/docs/git-send-email.html
* If your mail client supports setting the In-Reply-To header
via mailto: links, try the mailto: link
Be sure your reply has a Subject: header at the top and a blank line
before the message body.
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).