From e1d8c942fb732f0b21ccc78825f7af4df8d17f5e Mon Sep 17 00:00:00 2001
From: Eylon Krause <eylon1909@gmail.com>
Date: Tue, 11 Aug 2026 19:18:42 +0300
Subject: [PATCH] contrib/math: preserve negative zero in Expm1

Expm1(-0.0) returned +0.0 instead of -0.0. On the |x| < ln2/2 polynomial
path, ExpPoly forms 0.5 * (+0) + (-0) = +0, dropping the sign bit, and
Expm1 -- unlike the odd functions in this file -- never reapplies the
sign of x. C11 Annex F (F.10.3.3) requires expm1(+/-0) = +/-0.

Reapply x's sign on the polynomial path via Or(y, And(x, kNegZero)),
using the sign-bit mask already defined in this function. Since expm1(x)
shares the sign of x, this only affects the x == -0.0 case; every other
input (including +0, small nonzero x, and the large-|x| branch) is
unchanged.

Add a bit-exact signed-zero regression test for Expm1; the existing ULP
comparator treats +0 == -0 and cannot catch this.
---
 hwy/contrib/math/math-inl.h   |  5 ++++-
 hwy/contrib/math/math_test.cc | 20 ++++++++++++++++++++
 2 files changed, 24 insertions(+), 1 deletion(-)

diff --git a/hwy/contrib/math/math-inl.h b/hwy/contrib/math/math-inl.h
index b3d9b86de3..24878de6aa 100644
--- a/hwy/contrib/math/math-inl.h
+++ b/hwy/contrib/math/math-inl.h
@@ -2660,7 +2660,10 @@ HWY_INLINE V Expm1(const D d, V x) {
 
   // Reduce, approximate, and then reconstruct.
   const V y = impl.ExpPoly(d, impl.ExpReduce(d, x, q));
-  const V z = IfThenElse(Lt(Abs(x), kLn2Over2), y,
+  // Reapply the sign of x on the polynomial path so a negative-zero input is
+  // preserved: expm1(-0) = -0 (C11 F.10.3.3), which the polynomial drops when
+  // it forms 0.5 * (+0) + (-0) = +0.
+  const V z = IfThenElse(Lt(Abs(x), kLn2Over2), Or(y, And(x, kNegZero)),
                          Sub(impl.LoadExpShortRange(d, Add(y, kOne), q), kOne));
   return IfThenElse(Lt(x, kLowerBound), kNegOne, z);
 }
diff --git a/hwy/contrib/math/math_test.cc b/hwy/contrib/math/math_test.cc
index 45d3999428..a25c7ba6a9 100644
--- a/hwy/contrib/math/math_test.cc
+++ b/hwy/contrib/math/math_test.cc
@@ -191,6 +191,25 @@ HWY_NOINLINE void TestAllPow() {
   ForFloat3264Types(ForPartialVectors<TestPow>());
 }
 
+// expm1(+/-0) must preserve the sign of zero (C11 F.10.3.3). The ULP comparator
+// used by the other math tests treats +0 == -0, so check the bit pattern.
+struct TestExpm1SignedZero {
+  template <class T, class D>
+  HWY_NOINLINE void operator()(T, D d) {
+    using TU = MakeUnsigned<T>;
+    const T pos0 = ConvertScalarTo<T>(0.0);
+    const T neg0 = ConvertScalarTo<T>(-0.0);
+    const T got_pos = GetLane(CallExpm1(d, Set(d, pos0)));
+    const T got_neg = GetLane(CallExpm1(d, Set(d, neg0)));
+    HWY_ASSERT_EQ(BitCastScalar<TU>(pos0), BitCastScalar<TU>(got_pos));
+    HWY_ASSERT_EQ(BitCastScalar<TU>(neg0), BitCastScalar<TU>(got_neg));
+  }
+};
+
+HWY_NOINLINE void TestAllExpm1SignedZero() {
+  ForFloat3264Types(ForPartialVectors<TestExpm1SignedZero>());
+}
+
 }  // namespace
 // NOLINTNEXTLINE(google-readability-namespace-comments)
 }  // namespace HWY_NAMESPACE
@@ -213,6 +232,7 @@ HWY_EXPORT_AND_TEST_P(HwyMathTest, TestAllCbrt);
 HWY_EXPORT_AND_TEST_P(HwyMathTest, TestAllTgamma);
 HWY_EXPORT_AND_TEST_P(HwyMathTest, TestAllLogGamma);
 HWY_EXPORT_AND_TEST_P(HwyMathTest, TestAllPow);
+HWY_EXPORT_AND_TEST_P(HwyMathTest, TestAllExpm1SignedZero);
 HWY_AFTER_TEST();
 }  // namespace
 }  // namespace hwy