[<prev] [next>] [thread-next>] [day] [month] [year] [list]
Message-ID: <20260828044649.1165602-1-splitterblue@gmail.com>
Date: Thu, 27 Aug 2026 23:46:49 -0500
From: SplitterBlue <splitterblue@...il.com>
To: musl@...ts.openwall.com
Cc: SplitterBlue <splitterblue@...il.com>
Subject: [PATCH] math: avoid four-way quadrant branch in atan2/atan2f
the quadrant fixup selected one of four expressions from the signs of
both arguments. the two sign tests are independent: sign(x) decides
whether the result is the reference angle or its reflection about pi,
and sign(y) only negates. folding sign(y) into the operands rather than
the result removes the dependence on sign(y) and leaves one branch on
sign(x).
the sign has to be applied by xor on the operands, not by negating the
rounded result. -(pi-a) and (a-pi) are the same real number but not the
same computation: nearest and toward-zero are sign-symmetric, the
directed modes are not, so negating after rounding rounds the wrong way
in quadrant 3 under FE_UPWARD and FE_DOWNWARD. (-pi)-(-a) is one
rounding of the same exact value as (a-pi), so results and exception
flags are unchanged in all four rounding modes.
---
tested by building musl twice from f21a9653, with and without this
patch, and comparing the two libcs: identical result bits and identical
fetestexcept() flags for atan2 and atan2f in all four rounding modes,
over 1.2M random bit patterns for both operands, a special-value cross
product (+-0, +-inf, nan, subnormals, DBL_MIN, DBL_MAX), and values
around the |y/x| < 0x1p-64 cutoff. libc-test math gives the same 21
pre-existing failures on both builds and no new ones.
the obvious form of this change, reflect and then copysign(z,y), is not
equivalent: it differs by 1 ulp in quadrant 3 under the directed modes,
for the reason above.
timings against the real atan2, best of 5, three repeats, on a core
ultra 9 285H, as percent change in ns/op:
atan2 atan2f
random signs -12.3% -14.0%
gaussian x and y -10.9% -13.3%
x>0, random sign(y) -22.7% -21.2%
y>0, random sign(x) +2.5% +1.8%
90% first quadrant -4.8% -1.9%
sorted by quadrant +1.0% +2.7%
the gain comes from removing a mispredicted branch, so it depends on
how much the caller's quadrants vary. with x>0 and random sign(y),
branch-misses per call fall from 1.53 to 1.00 and 21 cycles go away.
with y>0 and random sign(x) both versions mispredict on sign(x), the
misses are unchanged at 1.53, and the extra xor is a small loss.
also measured on an i5-10400F, and on a cortex-a53 with gcc 14.2 where
the gain is 2-7%: atan2 costs about 180ns there and the mispredict
penalty is small. gcc 13.3 and clang 18.1, at -O2, -O3 and -Os, with
and without -ffreestanding -fno-builtin -frounding-math. code size is
the same or smaller everywhere measured, by 13-32 bytes on x86_64.
the patch uses only asuint64/asdouble, which atan2.c already uses via
EXTRACT_WORDS, so it adds no representation assumption that the file
did not already make. only x86_64 and aarch64 were measured. on 32-bit
targets the 64-bit xor in atan2 needs an extra instruction pair, which
was not measured; atan2f is unaffected, it stays 32-bit.
src/math/atan2.c | 16 ++++++++++------
src/math/atan2f.c | 16 ++++++++++------
2 files changed, 20 insertions(+), 12 deletions(-)
diff --git a/src/math/atan2.c b/src/math/atan2.c
index 5a1903c..b21153c 100644
--- a/src/math/atan2.c
+++ b/src/math/atan2.c
@@ -47,6 +47,7 @@ double atan2(double y, double x)
{
double z;
uint32_t m,lx,ly,ix,iy;
+ uint64_t sign;
if (isnan(x) || isnan(y))
return x+y;
@@ -97,11 +98,14 @@ double atan2(double y, double x)
z = 0;
else
z = atan(fabs(y/x));
- switch (m) {
- case 0: return z; /* atan(+,+) */
- case 1: return -z; /* atan(-,+) */
- case 2: return pi - (z-pi_lo); /* atan(+,-) */
- default: /* case 3 */
- return (z-pi_lo) - pi; /* atan(-,-) */
+ /* Fold sign(y) into both operands rather than selecting one of four
+ expressions. (asdouble(-pi) - asdouble(-a)) has the same exact value as
+ (a - pi), so the final subtraction is a single rounding of the same
+ quantity as before under every rounding mode. */
+ sign = (uint64_t)(m&1) << 63;
+ if (m & 2) {
+ z -= pi_lo;
+ return asdouble(asuint64(pi) ^ sign) - asdouble(asuint64(z) ^ sign);
}
+ return asdouble(asuint64(z) ^ sign);
}
diff --git a/src/math/atan2f.c b/src/math/atan2f.c
index c634d00..824d12c 100644
--- a/src/math/atan2f.c
+++ b/src/math/atan2f.c
@@ -23,6 +23,7 @@ float atan2f(float y, float x)
{
float z;
uint32_t m,ix,iy;
+ uint32_t sign;
if (isnan(x) || isnan(y))
return x+y;
@@ -73,11 +74,14 @@ float atan2f(float y, float x)
z = 0.0;
else
z = atanf(fabsf(y/x));
- switch (m) {
- case 0: return z; /* atan(+,+) */
- case 1: return -z; /* atan(-,+) */
- case 2: return pi - (z-pi_lo); /* atan(+,-) */
- default: /* case 3 */
- return (z-pi_lo) - pi; /* atan(-,-) */
+ /* Fold sign(y) into both operands rather than selecting one of four
+ expressions. (asfloat(-pi) - asfloat(-a)) has the same exact value as
+ (a - pi), so the final subtraction is a single rounding of the same
+ quantity as before under every rounding mode. */
+ sign = (uint32_t)(m&1) << 31;
+ if (m & 2) {
+ z -= pi_lo;
+ return asfloat(asuint(pi) ^ sign) - asfloat(asuint(z) ^ sign);
}
+ return asfloat(asuint(z) ^ sign);
}
--
2.43.0
Powered by blists - more mailing lists
Confused about mailing lists and their use?
Read about mailing lists on Wikipedia
and check out these
guidelines on proper formatting of your messages.
