From 18b60cc288c853f36d77211049f79af75c9ad9f5 Mon Sep 17 00:00:00 2001
From: ianpike <ianpike98@gmail.com>
Date: Fri, 19 Jun 2026 15:32:55 -0400
Subject: [PATCH 1/7] Fix constexpr remainder and make fmod exact for large
 quotients

---
 .../math/runtime/func/basic/fmod_rt.hpp       |  11 +-
 .../math/runtime/func/basic/remainder_rt.hpp  |  20 +--
 include/ccmath/math/basic/fmod.hpp            |  13 +-
 include/ccmath/math/basic/impl/CMakeLists.txt |   2 +
 .../math/basic/impl/fmod_double_impl.hpp      | 121 ++++++++++++++++++
 .../math/basic/impl/fmod_float_impl.hpp       | 121 ++++++++++++++++++
 include/ccmath/math/basic/remainder.hpp       |  30 ++---
 tests/src/math/basic/fmod_test.cpp            |  15 +++
 tests/src/math/basic/remainder_test.cpp       | 120 +++++++++++++++--
 tests/src/math/basic/remquo_test.cpp          |   5 +
 10 files changed, 410 insertions(+), 48 deletions(-)
 create mode 100644 include/ccmath/math/basic/impl/fmod_double_impl.hpp
 create mode 100644 include/ccmath/math/basic/impl/fmod_float_impl.hpp

diff --git a/include/ccmath/internal/math/runtime/func/basic/fmod_rt.hpp b/include/ccmath/internal/math/runtime/func/basic/fmod_rt.hpp
index 2b68f332..1e3d6ff5 100644
--- a/include/ccmath/internal/math/runtime/func/basic/fmod_rt.hpp
+++ b/include/ccmath/internal/math/runtime/func/basic/fmod_rt.hpp
@@ -12,8 +12,9 @@
 
 #include "ccmath/internal/math/generic/builtins/basic/fmod.hpp"
 #include "ccmath/internal/math/runtime/func/detail/system_math.hpp"
-#include "ccmath/internal/math/runtime/func/detail/trunc_scalar.hpp"
 #include "ccmath/internal/math/runtime/func/rt_dispatch.hpp"
+#include "ccmath/math/basic/impl/fmod_double_impl.hpp"
+#include "ccmath/math/basic/impl/fmod_float_impl.hpp"
 
 #include <type_traits>
 
@@ -28,7 +29,13 @@ namespace ccm::rt
 		if constexpr (ccm::builtin::has_runtime_fmod<T>) { return ccm::builtin::fmod_rt(x, y); }
 		else
 		{
-			return static_cast<T>(x - (detail::trunc_scalar<T>(x / y) * y));
+			// No builtin and no system math, so reuse the exact fdlibm bit-reduction. It is exact for
+			// every magnitude, unlike the trunc formula the other no-builtin basic fallbacks share.
+			if constexpr (std::is_same_v<T, float>) { return ccm::internal::fmod_float(x, y); }
+			else
+			{
+				return static_cast<T>(ccm::internal::fmod_double(static_cast<double>(x), static_cast<double>(y)));
+			}
 		}
 #endif
 	}
diff --git a/include/ccmath/internal/math/runtime/func/basic/remainder_rt.hpp b/include/ccmath/internal/math/runtime/func/basic/remainder_rt.hpp
index 786e52e8..c9747052 100644
--- a/include/ccmath/internal/math/runtime/func/basic/remainder_rt.hpp
+++ b/include/ccmath/internal/math/runtime/func/basic/remainder_rt.hpp
@@ -12,12 +12,9 @@
 
 #include "ccmath/internal/math/generic/builtins/basic/remainder.hpp"
 #include "ccmath/internal/math/runtime/func/detail/system_math.hpp"
-#include "ccmath/internal/math/runtime/func/detail/trunc_scalar.hpp"
 #include "ccmath/internal/math/runtime/func/rt_dispatch.hpp"
-#include "ccmath/internal/predef/unlikely.hpp"
-#include "ccmath/internal/support/fp/fp_bits.hpp"
+#include "ccmath/math/basic/remquo.hpp"
 
-#include <limits>
 #include <type_traits>
 
 namespace ccm::rt
@@ -31,16 +28,11 @@ namespace ccm::rt
 		if constexpr (ccm::builtin::has_runtime_remainder<T>) { return ccm::builtin::remainder_rt(x, y); }
 		else
 		{
-			using FPBits_t = typename ccm::support::fp::FPBits<T>;
-			const FPBits_t x_bits(x);
-			const FPBits_t y_bits(y);
-			const bool x_is_nan = x_bits.is_nan();
-			const bool y_is_nan = y_bits.is_nan();
-			if (CCM_UNLIKELY((x_bits.is_inf() && !y_is_nan) || (y_bits.is_zero() && !x_is_nan) || (x_is_nan || y_is_nan)))
-			{
-				return -std::numeric_limits<T>::quiet_NaN();
-			}
-			return static_cast<T>(x - (detail::trunc_scalar<T>(x / y) * y));
+			// No builtin and no system math, so reuse the exact remquo reduction. It rounds the
+			// quotient to nearest, ties-to-even, and covers the special cases, unlike the trunc
+			// formula that the other no-builtin basic fallbacks share.
+			int quotient = 0;
+			return ccm::remquo<T>(x, y, &quotient);
 		}
 #endif
 	}
diff --git a/include/ccmath/math/basic/fmod.hpp b/include/ccmath/math/basic/fmod.hpp
index ba9f8ba7..85e68503 100644
--- a/include/ccmath/math/basic/fmod.hpp
+++ b/include/ccmath/math/basic/fmod.hpp
@@ -15,7 +15,8 @@
 #include "ccmath/internal/predef/unlikely.hpp"
 #include "ccmath/internal/support/fp/fp_bits.hpp"
 #include "ccmath/internal/support/is_constant_evaluated.hpp"
-#include "ccmath/math/nearest/trunc.hpp"
+#include "ccmath/math/basic/impl/fmod_double_impl.hpp"
+#include "ccmath/math/basic/impl/fmod_float_impl.hpp"
 
 #include <limits>
 
@@ -68,7 +69,15 @@ namespace ccm
 				}
 			}
 
-			return static_cast<T>(x - (ccm::trunc<T>(x / y) * y));
+			// Exact, magnitude-independent reduction via the fdlibm integer bit-reduction. The old
+			// x - trunc(x / y) * y formula was only exact while x / y stayed representable, so it lost
+			// low bits once abs(x / y) reached 2^53. long double delegates to the double kernel, matching
+			// the remquol convention.
+			if constexpr (std::is_same_v<T, float>) { return internal::fmod_float(x, y); }
+			else
+			{
+				return static_cast<T>(internal::fmod_double(static_cast<double>(x), static_cast<double>(y)));
+			}
 		}
 
 		template <typename T, typename U, typename TC = std::common_type_t<T, U>>
diff --git a/include/ccmath/math/basic/impl/CMakeLists.txt b/include/ccmath/math/basic/impl/CMakeLists.txt
index 44fe9fa9..77099abe 100644
--- a/include/ccmath/math/basic/impl/CMakeLists.txt
+++ b/include/ccmath/math/basic/impl/CMakeLists.txt
@@ -1,4 +1,6 @@
 ccm_add_headers(
+        fmod_double_impl.hpp
+        fmod_float_impl.hpp
         nan_double_impl.hpp
         nan_float_impl.hpp
         nan_ldouble_impl.hpp
diff --git a/include/ccmath/math/basic/impl/fmod_double_impl.hpp b/include/ccmath/math/basic/impl/fmod_double_impl.hpp
new file mode 100644
index 00000000..9e6b98ed
--- /dev/null
+++ b/include/ccmath/math/basic/impl/fmod_double_impl.hpp
@@ -0,0 +1,121 @@
+/*
+ * Copyright (c) Ian Pike
+ * Copyright (c) CCMath contributors
+ *
+ * CCMath is provided under the Apache-2.0 License WITH LLVM-exception.
+ * See LICENSE for more information.
+ *
+ * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+ */
+
+#pragma once
+
+#include "ccmath/internal/predef/unlikely.hpp"
+#include "ccmath/internal/support/bits.hpp"
+
+#include <cstdint>
+
+namespace ccm::internal
+{
+	namespace impl
+	{
+		// Exact floating-point remainder for double using the fdlibm integer bit-reduction.
+		// The quotient is reduced one binary digit at a time on the raw significand integers, so
+		// the result is exact for every magnitude and is independent of the active rounding mode,
+		// unlike the x - trunc(x / y) * y formula which loses bits once x / y is not representable.
+		// NOLINTNEXTLINE(readability-function-cognitive-complexity)
+		constexpr double fmod_double_impl(double x, double y) noexcept
+		{
+			std::int64_t hx = support::double_to_int64(x);
+			std::int64_t hy = support::double_to_int64(y);
+
+			// Sign of the result follows x.
+			const std::int64_t sx = hx & static_cast<std::int64_t>(0x8000000000000000ULL);
+
+			hx &= 0x7fffffffffffffffLL; // |x|
+			hy &= 0x7fffffffffffffffLL; // |y|
+
+			// Purge exceptional inputs: y is zero, x is inf or NaN, or y is NaN. The (x * y) / (x * y)
+			// form yields the IEEE-mandated NaN for these cases.
+			if (CCM_UNLIKELY(hy == 0 || hx >= 0x7ff0000000000000LL || hy > 0x7ff0000000000000LL)) { return (x * y) / (x * y); }
+
+			// |x| < |y| leaves x unchanged (this also covers x == 0). |x| == |y| gives signed zero.
+			if (hx < hy) { return x; }
+			if (hx == hy) { return support::int64_to_double(sx); }
+
+			// Unbiased exponent of x, handling subnormals by normalizing the leading significand bit.
+			int ix = 0;
+			if (hx < 0x0010000000000000LL)
+			{
+				ix = -1022;
+				for (std::int64_t i = hx << 11; i > 0; i <<= 1) { ix -= 1; }
+			}
+			else
+			{
+				ix = static_cast<int>(hx >> 52) - 1023;
+			}
+
+			// Unbiased exponent of y.
+			int iy = 0;
+			if (hy < 0x0010000000000000LL)
+			{
+				iy = -1022;
+				for (std::int64_t i = hy << 11; i > 0; i <<= 1) { iy -= 1; }
+			}
+			else
+			{
+				iy = static_cast<int>(hy >> 52) - 1023;
+			}
+
+			// Promote both significands to integers with the implicit bit made explicit (subnormals are
+			// shifted up so their leading set bit sits in the same position a normal significand would).
+			if (ix >= -1022) { hx = 0x0010000000000000LL | (0x000fffffffffffffLL & hx); }
+			else
+			{
+				hx <<= (-1022 - ix);
+			}
+			if (iy >= -1022) { hy = 0x0010000000000000LL | (0x000fffffffffffffLL & hy); }
+			else
+			{
+				hy <<= (-1022 - iy);
+			}
+
+			// Fixed-point remainder: shift-and-subtract for each binary digit of the quotient.
+			int n = ix - iy;
+			while (n--)
+			{
+				const std::int64_t hz = hx - hy;
+				if (hz < 0) { hx = hx + hx; }
+				else
+				{
+					if (hz == 0) { return support::int64_to_double(sx); }
+					hx = hz + hz;
+				}
+			}
+			const std::int64_t hz = hx - hy;
+			if (hz >= 0) { hx = hz; }
+
+			// An exact zero remainder keeps the sign of x.
+			if (hx == 0) { return support::int64_to_double(sx); }
+
+			// Renormalize the remainder back into a floating-point significand.
+			while (hx < 0x0010000000000000LL)
+			{
+				hx = hx + hx;
+				iy -= 1;
+			}
+			if (iy >= -1022) // normal result
+			{
+				hx = (hx - 0x0010000000000000LL) | (static_cast<std::int64_t>(iy + 1023) << 52);
+			}
+			else // subnormal result
+			{
+				hx >>= (-1022 - iy);
+			}
+			return support::int64_to_double(hx | sx);
+		}
+	} // namespace impl
+
+	constexpr double fmod_double(double x, double y) noexcept
+	{ return impl::fmod_double_impl(x, y); }
+} // namespace ccm::internal
diff --git a/include/ccmath/math/basic/impl/fmod_float_impl.hpp b/include/ccmath/math/basic/impl/fmod_float_impl.hpp
new file mode 100644
index 00000000..314e2630
--- /dev/null
+++ b/include/ccmath/math/basic/impl/fmod_float_impl.hpp
@@ -0,0 +1,121 @@
+/*
+ * Copyright (c) Ian Pike
+ * Copyright (c) CCMath contributors
+ *
+ * CCMath is provided under the Apache-2.0 License WITH LLVM-exception.
+ * See LICENSE for more information.
+ *
+ * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+ */
+
+#pragma once
+
+#include "ccmath/internal/predef/unlikely.hpp"
+#include "ccmath/internal/support/bits.hpp"
+
+#include <cstdint>
+
+namespace ccm::internal
+{
+	namespace impl
+	{
+		// Exact floating-point remainder for float using the fdlibm integer bit-reduction.
+		// The quotient is reduced one binary digit at a time on the raw significand integers, so
+		// the result is exact for every magnitude and is independent of the active rounding mode,
+		// unlike the x - trunc(x / y) * y formula which loses bits once x / y is not representable.
+		// NOLINTNEXTLINE(readability-function-cognitive-complexity)
+		constexpr float fmod_float_impl(float x, float y) noexcept
+		{
+			std::int32_t hx = support::float_to_int32(x);
+			std::int32_t hy = support::float_to_int32(y);
+
+			// Sign of the result follows x.
+			const std::int32_t sx = hx & static_cast<std::int32_t>(0x80000000U);
+
+			hx &= 0x7fffffff; // |x|
+			hy &= 0x7fffffff; // |y|
+
+			// Purge exceptional inputs: y is zero, x is inf or NaN, or y is NaN. The (x * y) / (x * y)
+			// form yields the IEEE-mandated NaN for these cases.
+			if (CCM_UNLIKELY(hy == 0 || hx >= 0x7f800000 || hy > 0x7f800000)) { return (x * y) / (x * y); }
+
+			// |x| < |y| leaves x unchanged (this also covers x == 0). |x| == |y| gives signed zero.
+			if (hx < hy) { return x; }
+			if (hx == hy) { return support::int32_to_float(sx); }
+
+			// Unbiased exponent of x, handling subnormals by normalizing the leading significand bit.
+			int ix = 0;
+			if (hx < 0x00800000)
+			{
+				ix = -126;
+				for (std::int32_t i = hx << 8; i > 0; i <<= 1) { ix -= 1; }
+			}
+			else
+			{
+				ix = (hx >> 23) - 127;
+			}
+
+			// Unbiased exponent of y.
+			int iy = 0;
+			if (hy < 0x00800000)
+			{
+				iy = -126;
+				for (std::int32_t i = hy << 8; i > 0; i <<= 1) { iy -= 1; }
+			}
+			else
+			{
+				iy = (hy >> 23) - 127;
+			}
+
+			// Promote both significands to integers with the implicit bit made explicit (subnormals are
+			// shifted up so their leading set bit sits in the same position a normal significand would).
+			if (ix >= -126) { hx = 0x00800000 | (0x007fffff & hx); }
+			else
+			{
+				hx <<= (-126 - ix);
+			}
+			if (iy >= -126) { hy = 0x00800000 | (0x007fffff & hy); }
+			else
+			{
+				hy <<= (-126 - iy);
+			}
+
+			// Fixed-point remainder: shift-and-subtract for each binary digit of the quotient.
+			int n = ix - iy;
+			while (n--)
+			{
+				const std::int32_t hz = hx - hy;
+				if (hz < 0) { hx = hx + hx; }
+				else
+				{
+					if (hz == 0) { return support::int32_to_float(sx); }
+					hx = hz + hz;
+				}
+			}
+			const std::int32_t hz = hx - hy;
+			if (hz >= 0) { hx = hz; }
+
+			// An exact zero remainder keeps the sign of x.
+			if (hx == 0) { return support::int32_to_float(sx); }
+
+			// Renormalize the remainder back into a floating-point significand.
+			while (hx < 0x00800000)
+			{
+				hx = hx + hx;
+				iy -= 1;
+			}
+			if (iy >= -126) // normal result
+			{
+				hx = (hx - 0x00800000) | ((iy + 127) << 23);
+			}
+			else // subnormal result
+			{
+				hx >>= (-126 - iy);
+			}
+			return support::int32_to_float(hx | sx);
+		}
+	} // namespace impl
+
+	constexpr float fmod_float(float x, float y) noexcept
+	{ return impl::fmod_float_impl(x, y); }
+} // namespace ccm::internal
diff --git a/include/ccmath/math/basic/remainder.hpp b/include/ccmath/math/basic/remainder.hpp
index 44a35170..8147d740 100644
--- a/include/ccmath/math/basic/remainder.hpp
+++ b/include/ccmath/math/basic/remainder.hpp
@@ -10,12 +10,9 @@
 
 #pragma once
 
-#include "ccmath/internal/math/generic/builtins/basic/remainder.hpp"
 #include "ccmath/internal/math/runtime/func/basic/remainder_rt.hpp"
-#include "ccmath/internal/predef/unlikely.hpp"
-#include "ccmath/internal/support/fp/fp_bits.hpp"
 #include "ccmath/internal/support/is_constant_evaluated.hpp"
-#include "ccmath/math/nearest/trunc.hpp"
+#include "ccmath/math/basic/remquo.hpp"
 
 namespace ccm
 {
@@ -31,23 +28,14 @@ namespace ccm
 	{
 		if (!ccm::support::is_constant_evaluated()) { return ccm::rt::remainder_rt(x, y); }
 
-		using FPBits_t = typename ccm::support::fp::FPBits<T>;
-		const FPBits_t x_bits(x);
-		const FPBits_t y_bits(y);
-
-		const bool x_is_nan = x_bits.is_nan();
-		const bool y_is_nan = y_bits.is_nan();
-
-		// If x is ±∞ and y is not NaN, NaN is returned.
-		// If y is ±0 and x is not NaN, NaN is returned.
-		// If either argument is NaN, NaN is returned.
-		if (CCM_UNLIKELY((x_bits.is_inf() && !y_is_nan) || (y_bits.is_zero() && !x_is_nan) || (x_is_nan || y_is_nan)))
-		{
-			// All major compilers return -NaN.
-			return -std::numeric_limits<T>::quiet_NaN();
-		}
-
-		return static_cast<T>(x - (ccm::trunc<T>(x / y) * y));
+		// remainder is the remainder component of remquo, so reuse the exact iterative reduction
+		// already implemented there. The quotient rounds to nearest, ties-to-even, which keeps the
+		// result in the closed interval from -abs(y) / 2 to +abs(y) / 2. The result is an exact
+		// subtraction, so it does not depend on the active rounding mode, and the remquo kernel also
+		// covers the special cases (remainder of a finite value by an infinity is that finite value,
+		// and NaN is returned for an infinite dividend, a zero divisor, or a NaN argument).
+		int quotient = 0;
+		return ccm::remquo<T>(x, y, &quotient);
 	}
 
 	/**
diff --git a/tests/src/math/basic/fmod_test.cpp b/tests/src/math/basic/fmod_test.cpp
index 6297e937..78e03aa2 100644
--- a/tests/src/math/basic/fmod_test.cpp
+++ b/tests/src/math/basic/fmod_test.cpp
@@ -17,6 +17,21 @@
 #include <cmath>
 #include <limits>
 
+TEST(CcmathBasicTests, FmodLargeQuotientCompileTime)
+{
+	// The exact fdlibm bit-reduction must reduce these at compile time even though abs(x / y) is far
+	// above 2^53. The old constexpr x - trunc(x / y) * y formula collapsed to 0 here.
+	static_assert(ccm::fmod(1e30, 3.0) == 1.0, "fmod(1e30, 3) must be 1");
+	static_assert(ccm::fmod(-1e30, 3.0) == -1.0, "fmod(-1e30, 3) must be -1");
+	static_assert(ccm::fmod(1e300, 7.0) == 1.0, "fmod(1e300, 7) must reduce exactly to 1");
+	static_assert(ccm::fmod(1e30F, 3.0F) == 0.0F, "fmodf(1e30, 3) must reduce exactly to 0");
+
+	// Small / normal cases stay exact.
+	static_assert(ccm::fmod(10.0, 3.0) == 1.0, "fmod(10, 3) must be 1");
+	static_assert(ccm::fmod(-10.0, 3.0) == -1.0, "fmod(-10, 3) must be -1");
+	static_assert(ccm::fmod(7.5, 2.0) == 1.5, "fmod(7.5, 2) must be 1.5");
+}
+
 TEST(CcmathBasicTests, Fmod)
 {
 
diff --git a/tests/src/math/basic/remainder_test.cpp b/tests/src/math/basic/remainder_test.cpp
index 4cb7d91b..69da67f4 100644
--- a/tests/src/math/basic/remainder_test.cpp
+++ b/tests/src/math/basic/remainder_test.cpp
@@ -17,20 +17,122 @@
 #include <cmath>
 #include <limits>
 
+namespace
+{
+	constexpr double make_remainder(double x, double y)
+	{ return ccm::remainder(x, y); }
+
+	// Pulls the result of remainder(finite, inf) into a constant context.
+	constexpr double remainder_by_inf(double x)
+	{ return ccm::remainder(x, std::numeric_limits<double>::infinity()); }
+} // namespace
+
+TEST(CcmathBasicTests, RemainderCompileTime)
+{
+	// remainder rounds the quotient to nearest, ties-to-even, so these are the values the standard
+	// requires. The old trunc/fmod formula returned 2.0 and 1.5 here.
+	static_assert(make_remainder(5.0, 3.0) == -1.0, "remainder(5, 3) must be -1");
+	static_assert(make_remainder(7.5, 2.0) == -0.5, "remainder(7.5, 2) must be -0.5");
+
+	// A half-way quotient rounds to even.
+	static_assert(make_remainder(2.0, 1.0) == 0.0, "remainder(2, 1) must round the tie to 0");
+	static_assert(make_remainder(3.0, 2.0) == -1.0, "remainder(3, 2) must be -1");
+
+	// remainder of a finite value by an infinity is that finite value.
+	static_assert(remainder_by_inf(3.0) == 3.0, "remainder(3, inf) must be 3");
+
+	// abs(x / y) well above 2^53 must still reduce exactly at compile time. This is the residual that
+	// slipped through before fmod's generic kernel became the exact fdlibm reduction (remainder routes
+	// through remquo, which reduces with ccm::fmod).
+	static_assert(make_remainder(1e30, 3.0) == 1.0, "remainder(1e30, 3) must be 1");
+	static_assert(make_remainder(-1e30, 3.0) == -1.0, "remainder(-1e30, 3) must be -1");
+
+	static_assert(ccm::remainder(1.0, 1.0) == 0.0, "remainder(1, 1) must be 0");
+	static_assert(ccm::remainderf(5.0F, 3.0F) == -1.0F, "remainderf(5, 3) must be -1");
+}
+
+TEST(CcmathBasicTests, RemainderConfirmedCases)
+{
+	// The three confirmed-wrong cases from the H1 bug report, checked at runtime.
+	EXPECT_EQ(ccm::remainder(5.0, 3.0), -1.0);
+	EXPECT_EQ(ccm::remainder(7.5, 2.0), -0.5);
+
+	const double inf = std::numeric_limits<double>::infinity();
+	EXPECT_EQ(ccm::remainder(3.0, inf), 3.0);
+	EXPECT_EQ(ccm::remainder(-3.0, inf), -3.0);
+	EXPECT_EQ(ccm::remainder(3.0, -inf), 3.0);
+}
+
+TEST(CcmathBasicTests, RemainderSignOfZero)
+{
+	// remainder(+/-0, y) is +/-0 for finite nonzero y, sign of x preserved.
+	const double pos = ccm::remainder(0.0, 3.0);
+	const double neg = ccm::remainder(-0.0, 3.0);
+	EXPECT_EQ(pos, 0.0);
+	EXPECT_EQ(neg, 0.0);
+	EXPECT_FALSE(std::signbit(pos));
+	EXPECT_TRUE(std::signbit(neg));
+
+	const float posf = ccm::remainderf(0.0F, 3.0F);
+	const float negf = ccm::remainderf(-0.0F, 3.0F);
+	EXPECT_EQ(posf, 0.0F);
+	EXPECT_EQ(negf, 0.0F);
+	EXPECT_FALSE(std::signbit(posf));
+	EXPECT_TRUE(std::signbit(negf));
+}
+
+TEST(CcmathBasicTests, RemainderSpecialPropagation)
+{
+	const double inf = std::numeric_limits<double>::infinity();
+	const double nan = std::numeric_limits<double>::quiet_NaN();
+
+	// NaN is returned for an infinite dividend, a zero divisor, or a NaN argument.
+	EXPECT_TRUE(std::isnan(ccm::remainder(inf, 3.0)));
+	EXPECT_TRUE(std::isnan(ccm::remainder(-inf, 3.0)));
+	EXPECT_TRUE(std::isnan(ccm::remainder(3.0, 0.0)));
+	EXPECT_TRUE(std::isnan(ccm::remainder(3.0, -0.0)));
+	EXPECT_TRUE(std::isnan(ccm::remainder(nan, 3.0)));
+	EXPECT_TRUE(std::isnan(ccm::remainder(3.0, nan)));
+	EXPECT_TRUE(std::isnan(ccm::remainder(nan, nan)));
+}
+
+TEST(CcmathBasicTests, RemainderLargeQuotient)
+{
+	// |x / y| well above 2^53 must still reduce exactly. The trunc formula collapsed here.
+	ccm::test::ExpectBinaryMatchesStd(1e30, 3.0, ccm::remainder<double>, static_cast<double (*)(double, double)>(std::remainder));
+	ccm::test::ExpectBinaryMatchesStd(-1e30, 3.0, ccm::remainder<double>, static_cast<double (*)(double, double)>(std::remainder));
+	ccm::test::ExpectBinaryMatchesStd(1e300, 7.0, ccm::remainder<double>, static_cast<double (*)(double, double)>(std::remainder));
+	ccm::test::ExpectBinaryMatchesStd(1e30F, 3.0F, ccm::remainder<float>, static_cast<float (*)(float, float)>(std::remainder));
+}
+
 TEST(CcmathBasicTests, Remainder)
 {
 	static_assert(ccm::remainder(1.0, 1.0) == 0.0, "remainder has failed testing that it is static_assert-able!");
 
-	ccm::test::ExpectBinaryMatchesStd(1.0, 1.0, ccm::remainder<double>, static_cast<double (*)(double, double)>(std::remainder));
-	ccm::test::ExpectBinaryMatchesStd(1.0, 0.0, ccm::remainder<double>, static_cast<double (*)(double, double)>(std::remainder));
-	ccm::test::ExpectBinaryMatchesStd(0.0, 1.0, ccm::remainder<double>, static_cast<double (*)(double, double)>(std::remainder));
-	ccm::test::ExpectBinaryMatchesStd(0.0, 0.0, ccm::remainder<double>, static_cast<double (*)(double, double)>(std::remainder));
-	ccm::test::ExpectBinaryMatchesStd(-1.0, 1.0, ccm::remainder<double>, static_cast<double (*)(double, double)>(std::remainder));
-	ccm::test::ExpectBinaryMatchesStd(1.0, -1.0, ccm::remainder<double>, static_cast<double (*)(double, double)>(std::remainder));
-	ccm::test::ExpectBinaryMatchesStd(-1.0, -1.0, ccm::remainder<double>, static_cast<double (*)(double, double)>(std::remainder));
-	ccm::test::ExpectBinaryMatchesStd(-1.0, 0.0, ccm::remainder<double>, static_cast<double (*)(double, double)>(std::remainder));
+	using double_fn = double (*)(double, double);
+	using float_fn	= float (*)(float, float);
 
+	// A spread of normal values, both signs, for double.
+	const double doubles[][2] = {
+		{ 1.0, 1.0 },  { 1.0, 0.0 }, { 0.0, 1.0 },	 { 0.0, 0.0 },	  { -1.0, 1.0 },	{ 1.0, -1.0 }, { -1.0, -1.0 },
+		{ -1.0, 0.0 }, { 5.0, 3.0 }, { 7.5, 2.0 },	 { -5.0, 3.0 },	  { 5.0, -3.0 },	{ -7.5, 2.0 }, { 9.0, 4.0 },
+		{ 2.0, 1.0 },  { 3.0, 2.0 }, { 10.5, 3.25 }, { -10.5, 3.25 }, { 123.456, 7.0 }, { 0.1, 0.03 }, { -0.1, 0.03 },
+	};
+	for (const auto & c : doubles) { ccm::test::ExpectBinaryMatchesStd(c[0], c[1], ccm::remainder<double>, static_cast<double_fn>(std::remainder)); }
+
+	// A spread of normal values, both signs, for float.
+	const float floats[][2] = {
+		{ 1.0F, 1.0F }, { 5.0F, 3.0F }, { 7.5F, 2.0F },	  { -5.0F, 3.0F },	  { 5.0F, -3.0F },
+		{ 9.0F, 4.0F }, { 2.0F, 1.0F }, { 10.5F, 3.25F }, { 123.456F, 7.0F }, { 0.1F, 0.03F },
+	};
+	for (const auto & c : floats) { ccm::test::ExpectBinaryMatchesStd(c[0], c[1], ccm::remainder<float>, static_cast<float_fn>(std::remainder)); }
+
+	// Subnormal operands.
 	constexpr double subnormal_dividend = -5.0166534782602e-198;
 	constexpr double subnormal_divisor	= 5.27085811e-315;
-	ccm::test::ExpectBinaryMatchesStd(subnormal_dividend, subnormal_divisor, ccm::remainder<double>, static_cast<double (*)(double, double)>(std::remainder));
+	ccm::test::ExpectBinaryMatchesStd(subnormal_dividend, subnormal_divisor, ccm::remainder<double>, static_cast<double_fn>(std::remainder));
+
+	// Long double goes through the same public entry point.
+	ccm::test::ExpectBinaryMatchesStd(5.0L, 3.0L, ccm::remainderl, static_cast<long double (*)(long double, long double)>(std::remainder));
+	ccm::test::ExpectBinaryMatchesStd(7.5L, 2.0L, ccm::remainderl, static_cast<long double (*)(long double, long double)>(std::remainder));
 }
diff --git a/tests/src/math/basic/remquo_test.cpp b/tests/src/math/basic/remquo_test.cpp
index 60a3c9ca..cddf8659 100644
--- a/tests/src/math/basic/remquo_test.cpp
+++ b/tests/src/math/basic/remquo_test.cpp
@@ -55,6 +55,11 @@ TEST(CcmathBasicTests, Remquo)
 	static_assert(sa_quotient == -4, "sa_quotient == -4");
 	static_assert(sa_remainder == 1, "sa_quotient == 1");
 
+	// abs(x / y) well above 2^53 must reduce exactly at compile time. remquo coarse-reduces with
+	// ccm::fmod, so this locks in that the exact fdlibm fmod reduction propagates through remquo.
+	static_assert(get_ccm_rem(1e30, 3.0) == 1.0, "remquo(1e30, 3) remainder must be 1");
+	static_assert(get_ccm_rem(-1e30, 3.0) == -1.0, "remquo(-1e30, 3) remainder must be -1");
+
 	// Test with positive values
 	ccm::test::ExpectRemquoMatchesStd(7.0, 2.0);
 

From 327035189775c9ccee7c3a5f75ceb7518f0cb85f Mon Sep 17 00:00:00 2001
From: ianpike <ianpike98@gmail.com>
Date: Fri, 19 Jun 2026 16:03:30 -0400
Subject: [PATCH 2/7] Set a read-only default token permission for the lint
 workflow

---
 .github/workflows/lint.yml | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/.github/workflows/lint.yml b/.github/workflows/lint.yml
index ef24bd6e..899f67fc 100644
--- a/.github/workflows/lint.yml
+++ b/.github/workflows/lint.yml
@@ -24,6 +24,9 @@ on:
     branches:
       - '**'
 
+permissions:
+  contents: read
+
 concurrency:
   group: lint-${{ github.event.pull_request.number || github.ref }}
   cancel-in-progress: true

From 36002ba31c7a2ffa4c0661cd85f49517780905e4 Mon Sep 17 00:00:00 2001
From: ianpike <ianpike98@gmail.com>
Date: Fri, 19 Jun 2026 16:03:43 -0400
Subject: [PATCH 3/7] Normalize spacing in the freestanding meson option

---
 meson_options.txt | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/meson_options.txt b/meson_options.txt
index 518b27a0..82eda77f 100644
--- a/meson_options.txt
+++ b/meson_options.txt
@@ -4,6 +4,6 @@
 option('enable_runtime_simd', type: 'boolean', value: true, description: 'Enable SIMD optimization for runtime evaluation (does not affect compile-time)')
 option('disable_errno', type: 'boolean', value: false, description: 'Disable the use of errno in ccmath during runtime (may lead to faster evaluations but is non-standard)')
 option('disable_fenv', type: 'boolean', value: false, description: 'Completely disable the host floating-point environment: drop every <cfenv> and <fenv.h> include, assume round-to-nearest at runtime, and signal no fp-exceptions (auto-enabled when no host fenv header exists)')
-option('freestanding', type : 'boolean', value : false, description : 'Build for a freestanding C++ environment: restrict to freestanding-conformant headers, dropping the test-only <bitset> path in the runtime SIMD layer and the MSVC <math.h> system-math path (auto-enabled when __STDC_HOSTED__ is 0)')
+option('freestanding', type: 'boolean', value: false, description: 'Build for a freestanding C++ environment: restrict to freestanding-conformant headers, dropping the test-only <bitset> path in the runtime SIMD layer and the MSVC <math.h> system-math path (auto-enabled when __STDC_HOSTED__ is 0)')
 option('disable_reduced_precision_powl', type: 'boolean', value: false, description: 'Return quiet NaN from powl on non-binary80 long double instead of the default reduced-precision double fallback')
 option('deterministic', type: 'boolean', value: false, description: 'Produce bit-identical cross-hardware math: route transcendentals through the generic kernels (no libm), force the correctly-rounded FMA path, disable runtime SIMD, and evaluate long double in double precision (GCC/Clang)')

From e5a117f8255343149ca534d68de89764c47a9020 Mon Sep 17 00:00:00 2001
From: ianpike <ianpike98@gmail.com>
Date: Fri, 19 Jun 2026 19:57:59 -0400
Subject: [PATCH 4/7] Fix long double narrowing in the remquo dispatch

---
 include/ccmath/math/basic/remquo.hpp | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/include/ccmath/math/basic/remquo.hpp b/include/ccmath/math/basic/remquo.hpp
index d0eb7c14..e7636696 100644
--- a/include/ccmath/math/basic/remquo.hpp
+++ b/include/ccmath/math/basic/remquo.hpp
@@ -49,7 +49,9 @@ namespace ccm
 		if constexpr (std::is_same_v<T, float>) { return internal::remquo_float(x, y, quo); }
 		else
 		{
-			return internal::remquo_double(x, y, quo);
+			// double runs natively. long double delegates to the double kernel, so the narrowing
+			// is made explicit to stay clean under -Wconversion.
+			return static_cast<T>(internal::remquo_double(static_cast<double>(x), static_cast<double>(y), quo));
 		}
 	}
 

From 707f775aef52fa0cb325510f9d64f1e475098140 Mon Sep 17 00:00:00 2001
From: ianpike <ianpike98@gmail.com>
Date: Sat, 20 Jun 2026 01:43:16 -0400
Subject: [PATCH 5/7] Rework fmod

---
 .../math/runtime/func/basic/fmod_rt.hpp       |  10 +-
 include/ccmath/math/basic/fmod.hpp            |  14 +-
 include/ccmath/math/basic/impl/CMakeLists.txt |   3 +-
 .../math/basic/impl/fmod_double_impl.hpp      | 121 ------------------
 .../math/basic/impl/fmod_float_impl.hpp       | 121 ------------------
 include/ccmath/math/basic/impl/fmod_impl.hpp  |  83 ++++++++++++
 tests/src/math/basic/fmod_test.cpp            |   3 +-
 7 files changed, 95 insertions(+), 260 deletions(-)
 delete mode 100644 include/ccmath/math/basic/impl/fmod_double_impl.hpp
 delete mode 100644 include/ccmath/math/basic/impl/fmod_float_impl.hpp
 create mode 100644 include/ccmath/math/basic/impl/fmod_impl.hpp

diff --git a/include/ccmath/internal/math/runtime/func/basic/fmod_rt.hpp b/include/ccmath/internal/math/runtime/func/basic/fmod_rt.hpp
index 1e3d6ff5..ebbb7f3a 100644
--- a/include/ccmath/internal/math/runtime/func/basic/fmod_rt.hpp
+++ b/include/ccmath/internal/math/runtime/func/basic/fmod_rt.hpp
@@ -13,8 +13,7 @@
 #include "ccmath/internal/math/generic/builtins/basic/fmod.hpp"
 #include "ccmath/internal/math/runtime/func/detail/system_math.hpp"
 #include "ccmath/internal/math/runtime/func/rt_dispatch.hpp"
-#include "ccmath/math/basic/impl/fmod_double_impl.hpp"
-#include "ccmath/math/basic/impl/fmod_float_impl.hpp"
+#include "ccmath/math/basic/impl/fmod_impl.hpp"
 
 #include <type_traits>
 
@@ -29,12 +28,11 @@ namespace ccm::rt
 		if constexpr (ccm::builtin::has_runtime_fmod<T>) { return ccm::builtin::fmod_rt(x, y); }
 		else
 		{
-			// No builtin and no system math, so reuse the exact fdlibm bit-reduction. It is exact for
-			// every magnitude, unlike the trunc formula the other no-builtin basic fallbacks share.
-			if constexpr (std::is_same_v<T, float>) { return ccm::internal::fmod_float(x, y); }
+			// No builtin and no system math, so use the exact FPBits reduction. It is exact for every magnitude.
+			if constexpr (std::is_same_v<T, long double>) { return static_cast<T>(ccm::internal::fmod(static_cast<double>(x), static_cast<double>(y))); }
 			else
 			{
-				return static_cast<T>(ccm::internal::fmod_double(static_cast<double>(x), static_cast<double>(y)));
+				return ccm::internal::fmod(x, y);
 			}
 		}
 #endif
diff --git a/include/ccmath/math/basic/fmod.hpp b/include/ccmath/math/basic/fmod.hpp
index 85e68503..3dc9615f 100644
--- a/include/ccmath/math/basic/fmod.hpp
+++ b/include/ccmath/math/basic/fmod.hpp
@@ -15,8 +15,7 @@
 #include "ccmath/internal/predef/unlikely.hpp"
 #include "ccmath/internal/support/fp/fp_bits.hpp"
 #include "ccmath/internal/support/is_constant_evaluated.hpp"
-#include "ccmath/math/basic/impl/fmod_double_impl.hpp"
-#include "ccmath/math/basic/impl/fmod_float_impl.hpp"
+#include "ccmath/math/basic/impl/fmod_impl.hpp"
 
 #include <limits>
 
@@ -69,14 +68,13 @@ namespace ccm
 				}
 			}
 
-			// Exact, magnitude-independent reduction via the fdlibm integer bit-reduction. The old
-			// x - trunc(x / y) * y formula was only exact while x / y stayed representable, so it lost
-			// low bits once abs(x / y) reached 2^53. long double delegates to the double kernel, matching
-			// the remquol convention.
-			if constexpr (std::is_same_v<T, float>) { return internal::fmod_float(x, y); }
+			// Exact, magnitude-independent reduction over FPBits that gives the same result in every
+			// rounding mode. long double reduces through the double kernel, matching the fmodl and
+			// remquol convention.
+			if constexpr (std::is_same_v<T, float>) { return internal::fmod(x, y); }
 			else
 			{
-				return static_cast<T>(internal::fmod_double(static_cast<double>(x), static_cast<double>(y)));
+				return static_cast<T>(internal::fmod(static_cast<double>(x), static_cast<double>(y)));
 			}
 		}
 
diff --git a/include/ccmath/math/basic/impl/CMakeLists.txt b/include/ccmath/math/basic/impl/CMakeLists.txt
index 77099abe..80969919 100644
--- a/include/ccmath/math/basic/impl/CMakeLists.txt
+++ b/include/ccmath/math/basic/impl/CMakeLists.txt
@@ -1,6 +1,5 @@
 ccm_add_headers(
-        fmod_double_impl.hpp
-        fmod_float_impl.hpp
+        fmod_impl.hpp
         nan_double_impl.hpp
         nan_float_impl.hpp
         nan_ldouble_impl.hpp
diff --git a/include/ccmath/math/basic/impl/fmod_double_impl.hpp b/include/ccmath/math/basic/impl/fmod_double_impl.hpp
deleted file mode 100644
index 9e6b98ed..00000000
--- a/include/ccmath/math/basic/impl/fmod_double_impl.hpp
+++ /dev/null
@@ -1,121 +0,0 @@
-/*
- * Copyright (c) Ian Pike
- * Copyright (c) CCMath contributors
- *
- * CCMath is provided under the Apache-2.0 License WITH LLVM-exception.
- * See LICENSE for more information.
- *
- * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
- */
-
-#pragma once
-
-#include "ccmath/internal/predef/unlikely.hpp"
-#include "ccmath/internal/support/bits.hpp"
-
-#include <cstdint>
-
-namespace ccm::internal
-{
-	namespace impl
-	{
-		// Exact floating-point remainder for double using the fdlibm integer bit-reduction.
-		// The quotient is reduced one binary digit at a time on the raw significand integers, so
-		// the result is exact for every magnitude and is independent of the active rounding mode,
-		// unlike the x - trunc(x / y) * y formula which loses bits once x / y is not representable.
-		// NOLINTNEXTLINE(readability-function-cognitive-complexity)
-		constexpr double fmod_double_impl(double x, double y) noexcept
-		{
-			std::int64_t hx = support::double_to_int64(x);
-			std::int64_t hy = support::double_to_int64(y);
-
-			// Sign of the result follows x.
-			const std::int64_t sx = hx & static_cast<std::int64_t>(0x8000000000000000ULL);
-
-			hx &= 0x7fffffffffffffffLL; // |x|
-			hy &= 0x7fffffffffffffffLL; // |y|
-
-			// Purge exceptional inputs: y is zero, x is inf or NaN, or y is NaN. The (x * y) / (x * y)
-			// form yields the IEEE-mandated NaN for these cases.
-			if (CCM_UNLIKELY(hy == 0 || hx >= 0x7ff0000000000000LL || hy > 0x7ff0000000000000LL)) { return (x * y) / (x * y); }
-
-			// |x| < |y| leaves x unchanged (this also covers x == 0). |x| == |y| gives signed zero.
-			if (hx < hy) { return x; }
-			if (hx == hy) { return support::int64_to_double(sx); }
-
-			// Unbiased exponent of x, handling subnormals by normalizing the leading significand bit.
-			int ix = 0;
-			if (hx < 0x0010000000000000LL)
-			{
-				ix = -1022;
-				for (std::int64_t i = hx << 11; i > 0; i <<= 1) { ix -= 1; }
-			}
-			else
-			{
-				ix = static_cast<int>(hx >> 52) - 1023;
-			}
-
-			// Unbiased exponent of y.
-			int iy = 0;
-			if (hy < 0x0010000000000000LL)
-			{
-				iy = -1022;
-				for (std::int64_t i = hy << 11; i > 0; i <<= 1) { iy -= 1; }
-			}
-			else
-			{
-				iy = static_cast<int>(hy >> 52) - 1023;
-			}
-
-			// Promote both significands to integers with the implicit bit made explicit (subnormals are
-			// shifted up so their leading set bit sits in the same position a normal significand would).
-			if (ix >= -1022) { hx = 0x0010000000000000LL | (0x000fffffffffffffLL & hx); }
-			else
-			{
-				hx <<= (-1022 - ix);
-			}
-			if (iy >= -1022) { hy = 0x0010000000000000LL | (0x000fffffffffffffLL & hy); }
-			else
-			{
-				hy <<= (-1022 - iy);
-			}
-
-			// Fixed-point remainder: shift-and-subtract for each binary digit of the quotient.
-			int n = ix - iy;
-			while (n--)
-			{
-				const std::int64_t hz = hx - hy;
-				if (hz < 0) { hx = hx + hx; }
-				else
-				{
-					if (hz == 0) { return support::int64_to_double(sx); }
-					hx = hz + hz;
-				}
-			}
-			const std::int64_t hz = hx - hy;
-			if (hz >= 0) { hx = hz; }
-
-			// An exact zero remainder keeps the sign of x.
-			if (hx == 0) { return support::int64_to_double(sx); }
-
-			// Renormalize the remainder back into a floating-point significand.
-			while (hx < 0x0010000000000000LL)
-			{
-				hx = hx + hx;
-				iy -= 1;
-			}
-			if (iy >= -1022) // normal result
-			{
-				hx = (hx - 0x0010000000000000LL) | (static_cast<std::int64_t>(iy + 1023) << 52);
-			}
-			else // subnormal result
-			{
-				hx >>= (-1022 - iy);
-			}
-			return support::int64_to_double(hx | sx);
-		}
-	} // namespace impl
-
-	constexpr double fmod_double(double x, double y) noexcept
-	{ return impl::fmod_double_impl(x, y); }
-} // namespace ccm::internal
diff --git a/include/ccmath/math/basic/impl/fmod_float_impl.hpp b/include/ccmath/math/basic/impl/fmod_float_impl.hpp
deleted file mode 100644
index 314e2630..00000000
--- a/include/ccmath/math/basic/impl/fmod_float_impl.hpp
+++ /dev/null
@@ -1,121 +0,0 @@
-/*
- * Copyright (c) Ian Pike
- * Copyright (c) CCMath contributors
- *
- * CCMath is provided under the Apache-2.0 License WITH LLVM-exception.
- * See LICENSE for more information.
- *
- * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
- */
-
-#pragma once
-
-#include "ccmath/internal/predef/unlikely.hpp"
-#include "ccmath/internal/support/bits.hpp"
-
-#include <cstdint>
-
-namespace ccm::internal
-{
-	namespace impl
-	{
-		// Exact floating-point remainder for float using the fdlibm integer bit-reduction.
-		// The quotient is reduced one binary digit at a time on the raw significand integers, so
-		// the result is exact for every magnitude and is independent of the active rounding mode,
-		// unlike the x - trunc(x / y) * y formula which loses bits once x / y is not representable.
-		// NOLINTNEXTLINE(readability-function-cognitive-complexity)
-		constexpr float fmod_float_impl(float x, float y) noexcept
-		{
-			std::int32_t hx = support::float_to_int32(x);
-			std::int32_t hy = support::float_to_int32(y);
-
-			// Sign of the result follows x.
-			const std::int32_t sx = hx & static_cast<std::int32_t>(0x80000000U);
-
-			hx &= 0x7fffffff; // |x|
-			hy &= 0x7fffffff; // |y|
-
-			// Purge exceptional inputs: y is zero, x is inf or NaN, or y is NaN. The (x * y) / (x * y)
-			// form yields the IEEE-mandated NaN for these cases.
-			if (CCM_UNLIKELY(hy == 0 || hx >= 0x7f800000 || hy > 0x7f800000)) { return (x * y) / (x * y); }
-
-			// |x| < |y| leaves x unchanged (this also covers x == 0). |x| == |y| gives signed zero.
-			if (hx < hy) { return x; }
-			if (hx == hy) { return support::int32_to_float(sx); }
-
-			// Unbiased exponent of x, handling subnormals by normalizing the leading significand bit.
-			int ix = 0;
-			if (hx < 0x00800000)
-			{
-				ix = -126;
-				for (std::int32_t i = hx << 8; i > 0; i <<= 1) { ix -= 1; }
-			}
-			else
-			{
-				ix = (hx >> 23) - 127;
-			}
-
-			// Unbiased exponent of y.
-			int iy = 0;
-			if (hy < 0x00800000)
-			{
-				iy = -126;
-				for (std::int32_t i = hy << 8; i > 0; i <<= 1) { iy -= 1; }
-			}
-			else
-			{
-				iy = (hy >> 23) - 127;
-			}
-
-			// Promote both significands to integers with the implicit bit made explicit (subnormals are
-			// shifted up so their leading set bit sits in the same position a normal significand would).
-			if (ix >= -126) { hx = 0x00800000 | (0x007fffff & hx); }
-			else
-			{
-				hx <<= (-126 - ix);
-			}
-			if (iy >= -126) { hy = 0x00800000 | (0x007fffff & hy); }
-			else
-			{
-				hy <<= (-126 - iy);
-			}
-
-			// Fixed-point remainder: shift-and-subtract for each binary digit of the quotient.
-			int n = ix - iy;
-			while (n--)
-			{
-				const std::int32_t hz = hx - hy;
-				if (hz < 0) { hx = hx + hx; }
-				else
-				{
-					if (hz == 0) { return support::int32_to_float(sx); }
-					hx = hz + hz;
-				}
-			}
-			const std::int32_t hz = hx - hy;
-			if (hz >= 0) { hx = hz; }
-
-			// An exact zero remainder keeps the sign of x.
-			if (hx == 0) { return support::int32_to_float(sx); }
-
-			// Renormalize the remainder back into a floating-point significand.
-			while (hx < 0x00800000)
-			{
-				hx = hx + hx;
-				iy -= 1;
-			}
-			if (iy >= -126) // normal result
-			{
-				hx = (hx - 0x00800000) | ((iy + 127) << 23);
-			}
-			else // subnormal result
-			{
-				hx >>= (-126 - iy);
-			}
-			return support::int32_to_float(hx | sx);
-		}
-	} // namespace impl
-
-	constexpr float fmod_float(float x, float y) noexcept
-	{ return impl::fmod_float_impl(x, y); }
-} // namespace ccm::internal
diff --git a/include/ccmath/math/basic/impl/fmod_impl.hpp b/include/ccmath/math/basic/impl/fmod_impl.hpp
new file mode 100644
index 00000000..661347dc
--- /dev/null
+++ b/include/ccmath/math/basic/impl/fmod_impl.hpp
@@ -0,0 +1,83 @@
+/*
+ * Copyright (c) Ian Pike
+ * Copyright (c) CCMath contributors
+ *
+ * CCMath is provided under the Apache-2.0 License WITH LLVM-exception.
+ * See LICENSE for more information.
+ *
+ * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+ */
+
+#pragma once
+
+#include "ccmath/internal/predef/unlikely.hpp"
+#include "ccmath/internal/support/fp/fp_bits.hpp"
+#include "ccmath/math/fmanip/copysign.hpp"
+#include "ccmath/math/fmanip/scalbn.hpp"
+
+#include <cstdint>
+#include <limits>
+#include <type_traits>
+
+namespace ccm::internal
+{
+	namespace impl
+	{
+		// Exact floating-point remainder for any IEEE type, read straight off FPBits. Each magnitude is
+		// an integer significand scaled by a power of two, abs(v) == get_explicit_mantissa(v) *
+		// 2^(get_explicit_exponent(v) - mant_bits), an identity FPBits gives for normals and subnormals
+		// alike. The remainder of abs(x) by abs(y) is then
+		// ((sig_x * 2^gap) mod sig_y) * 2^(exp_y - mant_bits), with gap the difference of the two binade
+		// exponents. The modular product is built a chunk of bits at a time with 64-bit arithmetic, so
+		// the loop runs in gap / chunk steps instead of one step per bit, and the result stays exact for
+		// every magnitude and independent of the active rounding mode.
+		template <typename T>
+		constexpr T fmod_impl(T x, T y) noexcept
+		{
+			static_assert(std::is_floating_point_v<T>, "fmod_impl requires a floating-point type");
+			using FPBits_t = support::fp::FPBits<T>;
+
+			const FPBits_t x_bits(x);
+			const FPBits_t y_bits(y);
+
+			// NaN for an infinite dividend, a zero divisor, or a NaN argument.
+			if (CCM_UNLIKELY(x_bits.is_inf() || y_bits.is_zero() || x_bits.is_nan() || y_bits.is_nan())) { return FPBits_t::quiet_nan().get_val(); }
+
+			// Compare magnitudes through the exponent and significand bits. For finite values that integer
+			// order matches the floating-point order, so this also covers fmod(x, +/-inf) = x.
+			const auto x_mag = x_bits.abs().uintval();
+			const auto y_mag = y_bits.abs().uintval();
+
+			// Any abs(x) < abs(y) returns x unchanged.
+			if (x_mag < y_mag) { return x; }
+			// abs(x) == abs(y) leaves a zero that keeps the sign of x.
+			if (x_mag == y_mag) { return ccm::copysign(static_cast<T>(0), x); }
+
+			constexpr int mant_bits = std::numeric_limits<T>::digits - 1;
+			// Largest shift that keeps a reduced significand (below 2^(mant_bits + 1)) inside 64 bits.
+			constexpr int chunk = 63 - mant_bits;
+
+			const int exp_x			  = x_bits.get_explicit_exponent();
+			const int exp_y			  = y_bits.get_explicit_exponent();
+			const std::uint64_t sig_y = static_cast<std::uint64_t>(y_bits.get_explicit_mantissa());
+			std::uint64_t r			  = static_cast<std::uint64_t>(x_bits.get_explicit_mantissa()) % sig_y;
+
+			for (int gap = exp_x - exp_y; gap > 0;)
+			{
+				const int step = gap < chunk ? gap : chunk;
+				r			   = (r << step) % sig_y;
+				gap -= step;
+			}
+
+			if (r == 0) { return ccm::copysign(static_cast<T>(0), x); }
+			return ccm::copysign(ccm::scalbn(static_cast<T>(r), exp_y - mant_bits), x);
+		}
+	} // namespace impl
+
+	// Thin wrapper so the public dispatch reads the same as the other basic functions. float and
+	// double evaluate natively. long double is reduced through the double kernel by the caller,
+	// matching the existing fmodl and remquol convention.
+	template <typename T>
+	constexpr T fmod(T x, T y) noexcept
+	{ return impl::fmod_impl(x, y); }
+} // namespace ccm::internal
diff --git a/tests/src/math/basic/fmod_test.cpp b/tests/src/math/basic/fmod_test.cpp
index 78e03aa2..be251dd5 100644
--- a/tests/src/math/basic/fmod_test.cpp
+++ b/tests/src/math/basic/fmod_test.cpp
@@ -19,8 +19,7 @@
 
 TEST(CcmathBasicTests, FmodLargeQuotientCompileTime)
 {
-	// The exact fdlibm bit-reduction must reduce these at compile time even though abs(x / y) is far
-	// above 2^53. The old constexpr x - trunc(x / y) * y formula collapsed to 0 here.
+	// Large quotients must still reduce exactly at compile time, even when abs(x / y) is far above 2^53.
 	static_assert(ccm::fmod(1e30, 3.0) == 1.0, "fmod(1e30, 3) must be 1");
 	static_assert(ccm::fmod(-1e30, 3.0) == -1.0, "fmod(-1e30, 3) must be -1");
 	static_assert(ccm::fmod(1e300, 7.0) == 1.0, "fmod(1e300, 7) must reduce exactly to 1");

From d615122181c23f3547cc4ed437a6f5467df096e8 Mon Sep 17 00:00:00 2001
From: ianpike <ianpike98@gmail.com>
Date: Sun, 21 Jun 2026 13:16:35 -0400
Subject: [PATCH 6/7] Add fmod conformance tests and a binary64 kernel guard

---
 include/ccmath/math/basic/impl/fmod_impl.hpp  |   4 +
 tests/conformance/fenv_exception_test.cpp     |   3 +
 .../conformance/rounding_conformance_test.cpp |  21 +++
 tests/src/math/basic/fmod_test.cpp            | 166 +++++++++++++++---
 tests/src/math/basic/remquo_test.cpp          |  56 ++----
 5 files changed, 189 insertions(+), 61 deletions(-)

diff --git a/include/ccmath/math/basic/impl/fmod_impl.hpp b/include/ccmath/math/basic/impl/fmod_impl.hpp
index 661347dc..76424a19 100644
--- a/include/ccmath/math/basic/impl/fmod_impl.hpp
+++ b/include/ccmath/math/basic/impl/fmod_impl.hpp
@@ -35,6 +35,10 @@ namespace ccm::internal
 		constexpr T fmod_impl(T x, T y) noexcept
 		{
 			static_assert(std::is_floating_point_v<T>, "fmod_impl requires a floating-point type");
+			// The reduction lifts the significand into a uint64_t and shifts it by chunk = 63 - mant_bits
+			// bits, so it only supports types up to binary64. Wider types (80-bit or 128-bit long double)
+			// would give a non-positive chunk and must reduce through the double kernel at the call site.
+			static_assert(std::numeric_limits<T>::digits <= 53, "fmod_impl only supports types up to binary64");
 			using FPBits_t = support::fp::FPBits<T>;
 
 			const FPBits_t x_bits(x);
diff --git a/tests/conformance/fenv_exception_test.cpp b/tests/conformance/fenv_exception_test.cpp
index 5df4738a..9dc3b28a 100644
--- a/tests/conformance/fenv_exception_test.cpp
+++ b/tests/conformance/fenv_exception_test.cpp
@@ -42,6 +42,9 @@ TEST(CcmathFenvExceptionTests, DomainErrorsRaiseInvalidLikeStd)
 	ccm::test::ExpectFenvFlagsMatchStd([] { consume(ccm::log1p(runtime_value(-2.0))); }, [] { consume(std::log1p(runtime_value(-2.0))); }, FE_INVALID);
 	ccm::test::ExpectFenvFlagsMatchStd(
 		[] { consume(ccm::fmod(runtime_value(1.0), runtime_value(0.0))); }, [] { consume(std::fmod(runtime_value(1.0), runtime_value(0.0))); }, FE_INVALID);
+	ccm::test::ExpectFenvFlagsMatchStd([] { consume(ccm::fmod(runtime_value(std::numeric_limits<double>::infinity()), runtime_value(1.0))); },
+									   [] { consume(std::fmod(runtime_value(std::numeric_limits<double>::infinity()), runtime_value(1.0))); },
+									   FE_INVALID);
 	ccm::test::ExpectFenvFlagsMatchStd([] { consume(ccm::remainder(runtime_value(1.0), runtime_value(0.0))); },
 									   [] { consume(std::remainder(runtime_value(1.0), runtime_value(0.0))); },
 									   FE_INVALID);
diff --git a/tests/conformance/rounding_conformance_test.cpp b/tests/conformance/rounding_conformance_test.cpp
index b23b1e19..3721431d 100644
--- a/tests/conformance/rounding_conformance_test.cpp
+++ b/tests/conformance/rounding_conformance_test.cpp
@@ -168,6 +168,27 @@ TEST(CcmathRoundingConformanceTests, CeilTruncRoundFloorIndependentOfMode)
 	ccm::test::ExpectFpUnaryOverMatchesStdAllModes(ccm::test::samples::kNearbyIntProbeDouble, ccm::floor<double>, static_cast<double (*)(double)>(std::floor));
 }
 
+TEST(CcmathRoundingConformanceTests, FmodRemainderIndependentOfRoundingMode)
+{
+	// fmod and remainder are exact operations, so the result must not vary with the active rounding
+	// mode. Both operands go through runtime_value so the kernels read the dynamic mode at run time.
+	ccm::test::ForEachRoundingModeOrSkip(
+		[&](int mode)
+		{
+			ccm::test::ScopedRoundingMode scope(mode);
+			ASSERT_TRUE(scope.active());
+
+			ccm::test::ExpectFpEq(ccm::fmod(runtime_value(10.0), runtime_value(3.0)), 1.0);
+			ccm::test::ExpectFpEq(ccm::fmod(runtime_value(7.5), runtime_value(2.0)), 1.5);
+			ccm::test::ExpectFpEq(ccm::fmod(runtime_value(-10.0), runtime_value(3.0)), -1.0);
+			ccm::test::ExpectFpEq(ccm::fmod(runtime_value(1e30), runtime_value(3.0)), 1.0);
+			ccm::test::ExpectFpEq(ccm::internal::fmod(runtime_value(1e300), runtime_value(7.0)), 1.0);
+
+			ccm::test::ExpectFpEq(ccm::remainder(runtime_value(5.0), runtime_value(3.0)), -1.0);
+			ccm::test::ExpectFpEq(ccm::remainder(runtime_value(7.5), runtime_value(2.0)), -0.5);
+		});
+}
+
 // The following pin the FE_TONEAREST guard in the expo runtime headers: outside round to nearest
 // the runtime dispatch must fall back to the generic kernel, so it matches the generic wrapper bit
 // for bit over normal inputs.
diff --git a/tests/src/math/basic/fmod_test.cpp b/tests/src/math/basic/fmod_test.cpp
index be251dd5..890fa86e 100644
--- a/tests/src/math/basic/fmod_test.cpp
+++ b/tests/src/math/basic/fmod_test.cpp
@@ -9,13 +9,49 @@
  */
 
 #include "utils/std_compare.hpp"
+#include "utils/test_runtime.hpp"
 
 #include <gtest/gtest.h>
 
 #include <ccmath/ccmath.hpp>
 
 #include <cmath>
+#include <cstdint>
+#include <cstring>
 #include <limits>
+#include <random>
+
+using ccm::test::runtime_value;
+
+namespace
+{
+	// fmod is an exact operation, so every result must match the host libm bit for bit, not within a
+	// ULP band. ccm::fmod is the public entry (builtin or kernel); ccm::internal::fmod is the kernel
+	// the constexpr and no-builtin paths actually run.
+	void expect_public_exact(double x, double y)
+	{
+		SCOPED_TRACE(testing::Message() << "fmod(" << x << ", " << y << ")");
+		ccm::test::ExpectFpEq(ccm::fmod(x, y), std::fmod(x, y));
+	}
+	void expect_public_exact_f(float x, float y)
+	{
+		SCOPED_TRACE(testing::Message() << "fmodf(" << x << ", " << y << ")");
+		ccm::test::ExpectFpEq(ccm::fmod(x, y), std::fmod(x, y));
+	}
+	// Exercises the chunked-modulo kernel directly at runtime (the public path uses __builtin_fmod on
+	// GCC/Clang, so without this the kernel only has compile-time coverage). runtime_value defeats
+	// constant folding so the kernel really executes.
+	void expect_kernel_exact(double x, double y)
+	{
+		SCOPED_TRACE(testing::Message() << "kernel fmod(" << x << ", " << y << ")");
+		ccm::test::ExpectFpEq(ccm::internal::fmod(runtime_value(x), runtime_value(y)), std::fmod(x, y));
+	}
+	void expect_kernel_exact_f(float x, float y)
+	{
+		SCOPED_TRACE(testing::Message() << "kernel fmodf(" << x << ", " << y << ")");
+		ccm::test::ExpectFpEq(ccm::internal::fmod(runtime_value(x), runtime_value(y)), std::fmod(x, y));
+	}
+} // namespace
 
 TEST(CcmathBasicTests, FmodLargeQuotientCompileTime)
 {
@@ -29,35 +65,123 @@ TEST(CcmathBasicTests, FmodLargeQuotientCompileTime)
 	static_assert(ccm::fmod(10.0, 3.0) == 1.0, "fmod(10, 3) must be 1");
 	static_assert(ccm::fmod(-10.0, 3.0) == -1.0, "fmod(-10, 3) must be -1");
 	static_assert(ccm::fmod(7.5, 2.0) == 1.5, "fmod(7.5, 2) must be 1.5");
+	static_assert(ccm::fmod(1, 2) == 1, "fmod must be usable in a static_assert");
+}
+
+TEST(CcmathBasicTests, FmodBitExactVsStd)
+{
+	for (double y : { 3.0, -3.0, 2.5, 0.5, 1024.0, 1e-300 })
+	{
+		for (double x : { 10.0, -10.0, 7.5, -7.5, 0.0, -0.0, 1.5, 1e30, -1e30, 1e300 })
+		{
+			expect_public_exact(x, y);
+			expect_kernel_exact(x, y);
+		}
+	}
+	for (float y : { 3.0F, -3.0F, 2.5F, 0.5F })
+	{
+		for (float x : { 10.0F, -10.0F, 7.5F, 0.0F, -0.0F, 1e30F })
+		{
+			expect_public_exact_f(x, y);
+			expect_kernel_exact_f(x, y);
+		}
+	}
+}
+
+TEST(CcmathBasicTests, FmodKernelRandomBitExactVsStd)
+{
+	// Stress the runtime kernel against the host libm on random finite pairs, bit for bit.
+	std::mt19937_64 rng(0xF110D);
+	std::uniform_int_distribution<std::uint64_t> bd(0, ~std::uint64_t(0));
+	int checked = 0;
+	while (checked < 6000)
+	{
+		double x;
+		double y;
+		const std::uint64_t bx = bd(rng);
+		const std::uint64_t by = bd(rng);
+		std::memcpy(&x, &bx, sizeof x);
+		std::memcpy(&y, &by, sizeof y);
+		if (!std::isfinite(x) || !std::isfinite(y) || y == 0.0) { continue; }
+		ccm::test::ExpectFpEq(ccm::internal::fmod(x, y), std::fmod(x, y));
+		++checked;
+	}
+
+	std::mt19937 rngf(0xF110F);
+	std::uniform_int_distribution<std::uint32_t> bf(0, ~std::uint32_t(0));
+	checked = 0;
+	while (checked < 6000)
+	{
+		float x;
+		float y;
+		const std::uint32_t bx = bf(rngf);
+		const std::uint32_t by = bf(rngf);
+		std::memcpy(&x, &bx, sizeof x);
+		std::memcpy(&y, &by, sizeof y);
+		if (!std::isfinite(x) || !std::isfinite(y) || y == 0.0F) { continue; }
+		ccm::test::ExpectFpEq(ccm::internal::fmod(x, y), std::fmod(x, y));
+		++checked;
+	}
 }
 
-TEST(CcmathBasicTests, Fmod)
+TEST(CcmathBasicTests, FmodSignedZeroResult)
 {
+	// A zero remainder (exact multiple, or equal magnitudes) keeps the sign of x.
+	ccm::test::ExpectSignedZero(ccm::fmod(6.0, 3.0), false);
+	ccm::test::ExpectSignedZero(ccm::fmod(-6.0, 3.0), true);
+	ccm::test::ExpectSignedZero(ccm::fmod(3.0, 3.0), false);
+	ccm::test::ExpectSignedZero(ccm::fmod(-3.0, 3.0), true);
+	ccm::test::ExpectSignedZero(ccm::internal::fmod(runtime_value(6.0), runtime_value(3.0)), false);
+	ccm::test::ExpectSignedZero(ccm::internal::fmod(runtime_value(-6.0), runtime_value(3.0)), true);
+	ccm::test::ExpectSignedZero(ccm::fmod(6.0F, 3.0F), false);
+	ccm::test::ExpectSignedZero(ccm::fmod(-6.0F, 3.0F), true);
 
-	// Test that fmod works with static_assert
-	static_assert(ccm::fmod(1, 2) == 1, "fmod has failed testing that it is static_assert-able!");
+	// fmod(+/-0, y) returns +/-0 for non-zero y.
+	ccm::test::ExpectSignedZero(ccm::fmod(0.0, 3.0), false);
+	ccm::test::ExpectSignedZero(ccm::fmod(-0.0, 3.0), true);
+}
 
-	ccm::test::ExpectBinaryMatchesStd(10.0F, 3.0F, ccm::fmod<float>, static_cast<float (*)(float, float)>(std::fmod));
-	ccm::test::ExpectBinaryMatchesStd(10.0F, -3.0F, ccm::fmod<float>, static_cast<float (*)(float, float)>(std::fmod));
-	ccm::test::ExpectBinaryMatchesStd(-10.0F, 3.0F, ccm::fmod<float>, static_cast<float (*)(float, float)>(std::fmod));
-	ccm::test::ExpectBinaryMatchesStd(-10.0F, -3.0F, ccm::fmod<float>, static_cast<float (*)(float, float)>(std::fmod));
-	ccm::test::ExpectBinaryMatchesStd(0.0F, 3.0F, ccm::fmod<float>, static_cast<float (*)(float, float)>(std::fmod));
+TEST(CcmathBasicTests, FmodSpecialCases)
+{
+	constexpr double inf = std::numeric_limits<double>::infinity();
+	constexpr double nan = std::numeric_limits<double>::quiet_NaN();
 
-	// This is a tough test as it forces rounding precision issues.
-	// EXPECT_FLOAT_EQ(ccm::fmod(30.508474576271183309f, 6.1016949152542370172f), std::fmod(30.508474576271183309f, 6.1016949152542370172f));
+	// fmod(x, +/-inf) == x for finite x.
+	expect_public_exact(10.0, inf);
+	expect_public_exact(10.0, -inf);
+	expect_public_exact(-7.5, inf);
 
-	// Test fmod with integer numbers
-	ccm::test::ExpectSameAsStd(ccm::fmod(10, 3), std::fmod(10, 3));
+	// Domain errors yield NaN: fmod(x, 0), fmod(+/-inf, y), and NaN in either argument.
+	EXPECT_TRUE(std::isnan(ccm::fmod(10.0, 0.0)));
+	EXPECT_TRUE(std::isnan(ccm::fmod(inf, 3.0)));
+	EXPECT_TRUE(std::isnan(ccm::fmod(-inf, 3.0)));
+	EXPECT_TRUE(std::isnan(ccm::fmod(nan, 3.0)));
+	EXPECT_TRUE(std::isnan(ccm::fmod(3.0, nan)));
+	EXPECT_TRUE(std::isnan(ccm::internal::fmod(runtime_value(inf), runtime_value(3.0))));
+	EXPECT_TRUE(std::isnan(ccm::internal::fmod(runtime_value(10.0), runtime_value(0.0))));
+}
 
-	ccm::test::ExpectBinaryMatchesStd(1.5, 2.781342323457793e-309, ccm::fmod<double>, static_cast<double (*)(double, double)>(std::fmod));
+TEST(CcmathBasicTests, FmodSubnormal)
+{
+	constexpr double dmin = std::numeric_limits<double>::denorm_min();
+	expect_public_exact(1.5, 2.781342323457793e-309); // normal dividend, subnormal divisor
+	expect_kernel_exact(1.5, 2.781342323457793e-309);
+	expect_kernel_exact(5.0 * dmin, 2.0 * dmin); // subnormal dividend and divisor
+	expect_kernel_exact(7.0 * dmin, 3.0 * dmin);
+	expect_kernel_exact(1e300, dmin); // largest gap (normal over min subnormal)
+	expect_kernel_exact(std::numeric_limits<double>::max(), dmin);
+}
 
-	/// Test Edge Cases
+TEST(CcmathBasicTests, FmodTypesAndOverloads)
+{
+	// Integer overload promotes to double.
+	ccm::test::ExpectFpEq(ccm::fmod(10, 3), std::fmod(10, 3));
+	static_assert(std::is_same_v<decltype(ccm::fmod(10, 3)), double>, "integer fmod must return double");
 
-	ccm::test::ExpectBinaryMatchesStd(0.0F, 1.0F, ccm::fmod<float>, static_cast<float (*)(float, float)>(std::fmod));
-	ccm::test::ExpectBinaryMatchesStd(-0.0F, 1.0F, ccm::fmod<float>, static_cast<float (*)(float, float)>(std::fmod));
-	ccm::test::ExpectBinaryMatchesStd(10.0F, 0.0F, ccm::fmod<float>, static_cast<float (*)(float, float)>(std::fmod));
-	ccm::test::ExpectBinaryMatchesStd(10.0F, std::numeric_limits<float>::infinity(), ccm::fmod<float>, static_cast<float (*)(float, float)>(std::fmod));
-	ccm::test::ExpectBinaryMatchesStd(10.0F, -std::numeric_limits<float>::infinity(), ccm::fmod<float>, static_cast<float (*)(float, float)>(std::fmod));
-	ccm::test::ExpectBinaryMatchesStd(std::numeric_limits<float>::quiet_NaN(), 10.0F, ccm::fmod<float>, static_cast<float (*)(float, float)>(std::fmod));
-	ccm::test::ExpectBinaryMatchesStd(10.0F, std::numeric_limits<float>::quiet_NaN(), ccm::fmod<float>, static_cast<float (*)(float, float)>(std::fmod));
+	// fmodf / fmodl wrappers. long double delegates to the double kernel by design, so use values that
+	// are exact through double.
+	ccm::test::ExpectFpEq(ccm::fmodf(10.0F, 3.0F), std::fmod(10.0F, 3.0F));
+	EXPECT_EQ(ccm::fmodl(10.0L, 3.0L), 1.0L);
+	EXPECT_EQ(ccm::fmodl(7.5L, 2.0L), 1.5L);
+	EXPECT_EQ(ccm::fmodl(-10.0L, 3.0L), -1.0L);
 }
diff --git a/tests/src/math/basic/remquo_test.cpp b/tests/src/math/basic/remquo_test.cpp
index cddf8659..14ad3faf 100644
--- a/tests/src/math/basic/remquo_test.cpp
+++ b/tests/src/math/basic/remquo_test.cpp
@@ -56,7 +56,7 @@ TEST(CcmathBasicTests, Remquo)
 	static_assert(sa_remainder == 1, "sa_quotient == 1");
 
 	// abs(x / y) well above 2^53 must reduce exactly at compile time. remquo coarse-reduces with
-	// ccm::fmod, so this locks in that the exact fdlibm fmod reduction propagates through remquo.
+	// ccm::fmod, so this locks in that the exact fmod reduction propagates through remquo.
 	static_assert(get_ccm_rem(1e30, 3.0) == 1.0, "remquo(1e30, 3) remainder must be 1");
 	static_assert(get_ccm_rem(-1e30, 3.0) == -1.0, "remquo(-1e30, 3) remainder must be -1");
 
@@ -69,43 +69,19 @@ TEST(CcmathBasicTests, Remquo)
 	// Test with zero values
 	ccm::test::ExpectRemquoMatchesStd(0.0, 2.0);
 
-	/* TODO(IanP): Infinity and NaN remquo cases fail on CI but pass locally. Investigate.
-	// Test with infinity
-	bool isCcmLeftInfinityNegative = (std::signbit(get_ccm_rem(std::numeric_limits<double>::infinity(), 2.0)) == true); // NOLINT
-	bool isStdLeftInfinityNegative = (std::signbit(get_std_rem(std::numeric_limits<double>::infinity(), 2.0)) == true); // NOLINT
-	bool didCcmLeftInfinityReturnNan = std::isnan(get_ccm_rem(std::numeric_limits<double>::infinity(), 2.0));
-	bool didStdLeftInfinityReturnNan = std::isnan(get_std_rem(std::numeric_limits<double>::infinity(), 2.0));
-	EXPECT_EQ(isCcmLeftInfinityNegative, isStdLeftInfinityNegative);
-	EXPECT_EQ(didCcmLeftInfinityReturnNan, didStdLeftInfinityReturnNan);
-	EXPECT_EQ(get_ccm_quo(std::numeric_limits<double>::infinity(), 2.0), get_std_quo(std::numeric_limits<double>::infinity(), 2.0));
-
-	// Test with negative infinity
-	bool isCcmLeftNegativeInfinityNegative = (std::signbit(get_ccm_rem(-std::numeric_limits<double>::infinity(), 2.0)) == true); // NOLINT
-	bool isStdLeftNegativeInfinityNegative = (std::signbit(get_std_rem(-std::numeric_limits<double>::infinity(), 2.0)) == true); // NOLINT
-	bool didCcmLeftNegativeInfinityReturnNan = std::isnan(get_ccm_rem(-std::numeric_limits<double>::infinity(), 2.0));
-	bool didStdLeftNegativeInfinityReturnNan = std::isnan(get_std_rem(-std::numeric_limits<double>::infinity(), 2.0));
-	EXPECT_EQ(isCcmLeftNegativeInfinityNegative, isStdLeftNegativeInfinityNegative);
-	EXPECT_EQ(didCcmLeftNegativeInfinityReturnNan, didStdLeftNegativeInfinityReturnNan);
-	EXPECT_EQ(get_ccm_quo(-std::numeric_limits<double>::infinity(), 2.0), get_std_quo(-std::numeric_limits<double>::infinity(), 2.0));
-
-	// Test with NaN
-	bool isCcmLeftNanNegative = (std::signbit(get_ccm_rem(std::numeric_limits<double>::quiet_NaN(), 2.0)) == true &&
-	std::isnan(get_ccm_rem(std::numeric_limits<double>::quiet_NaN(), 2.0)) == true); // NOLINT bool isStdLeftNanNegative =
-	(std::signbit(get_std_rem(std::numeric_limits<double>::quiet_NaN(), 2.0)) == true && std::isnan(get_std_rem(std::numeric_limits<double>::quiet_NaN(), 2.0))
-	== true); // NOLINT bool didCcmLeftNanReturnNan = std::isnan(get_ccm_rem(std::numeric_limits<double>::quiet_NaN(), 2.0)); bool didStdLeftNanReturnNan =
-	std::isnan(get_std_rem(std::numeric_limits<double>::quiet_NaN(), 2.0)); EXPECT_EQ(isCcmLeftNanNegative, isStdLeftNanNegative);
-	EXPECT_EQ(didCcmLeftNanReturnNan, didStdLeftNanReturnNan);
-	EXPECT_EQ(get_ccm_quo(std::numeric_limits<double>::quiet_NaN(), 2.0), get_std_quo(std::numeric_limits<double>::quiet_NaN(), 2.0));
-
-	// Test with negative NaN
-	bool isCcmLeftNegativeNanNegative = (std::signbit(get_ccm_rem(-std::numeric_limits<double>::quiet_NaN(), 2.0)) == true &&
-	std::isnan(get_ccm_rem(-std::numeric_limits<double>::quiet_NaN(), 2.0)) == true); // NOLINT bool isStdLeftNegativeNanNegative =
-	(std::signbit(get_std_rem(-std::numeric_limits<double>::quiet_NaN(), 2.0)) == true &&
-	std::isnan(get_std_rem(-std::numeric_limits<double>::quiet_NaN(), 2.0)) == true); // NOLINT bool didCcmLeftNegativeNanReturnNan =
-	std::isnan(get_ccm_rem(-std::numeric_limits<double>::quiet_NaN(), 2.0)); bool didStdLeftNegativeNanReturnNan =
-	std::isnan(get_std_rem(-std::numeric_limits<double>::quiet_NaN(), 2.0)); EXPECT_EQ(isCcmLeftNegativeNanNegative, isStdLeftNegativeNanNegative);
-	EXPECT_EQ(didCcmLeftNegativeNanReturnNan, didStdLeftNegativeNanReturnNan);
-	EXPECT_EQ(get_ccm_quo(-std::numeric_limits<double>::quiet_NaN(), 2.0), get_std_quo(-std::numeric_limits<double>::quiet_NaN(), 2.0));
-
-	 */
+	// Domain cases: the standard returns a NaN remainder when x is infinite or y is zero, and
+	// propagates a NaN argument. The quotient (and the NaN sign) are unspecified in these cases, so only
+	// the NaN result is asserted. The earlier version compared quo and signbit here, which is what made
+	// it pass locally but fail on other platforms.
+	constexpr double inf = std::numeric_limits<double>::infinity();
+	constexpr double nan = std::numeric_limits<double>::quiet_NaN();
+	EXPECT_TRUE(std::isnan(get_ccm_rem(inf, 2.0)));
+	EXPECT_TRUE(std::isnan(get_ccm_rem(-inf, 2.0)));
+	EXPECT_TRUE(std::isnan(get_ccm_rem(2.0, 0.0)));
+	EXPECT_TRUE(std::isnan(get_ccm_rem(nan, 2.0)));
+	EXPECT_TRUE(std::isnan(get_ccm_rem(2.0, nan)));
+
+	// std agrees that these are NaN, confirming the shared contract.
+	EXPECT_TRUE(std::isnan(get_std_rem(inf, 2.0)));
+	EXPECT_TRUE(std::isnan(get_std_rem(2.0, 0.0)));
 }

From 5f6c1f6037c98f9d764bf54d387ae88dc047075c Mon Sep 17 00:00:00 2001
From: ianpike <ianpike98@gmail.com>
Date: Sun, 21 Jun 2026 13:48:58 -0400
Subject: [PATCH 7/7] Avoid function-style casts in the fmod test

---
 tests/src/math/basic/fmod_test.cpp | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/tests/src/math/basic/fmod_test.cpp b/tests/src/math/basic/fmod_test.cpp
index 890fa86e..fe79df19 100644
--- a/tests/src/math/basic/fmod_test.cpp
+++ b/tests/src/math/basic/fmod_test.cpp
@@ -92,7 +92,7 @@ TEST(CcmathBasicTests, FmodKernelRandomBitExactVsStd)
 {
 	// Stress the runtime kernel against the host libm on random finite pairs, bit for bit.
 	std::mt19937_64 rng(0xF110D);
-	std::uniform_int_distribution<std::uint64_t> bd(0, ~std::uint64_t(0));
+	std::uniform_int_distribution<std::uint64_t> bd(0, std::numeric_limits<std::uint64_t>::max());
 	int checked = 0;
 	while (checked < 6000)
 	{
@@ -108,7 +108,7 @@ TEST(CcmathBasicTests, FmodKernelRandomBitExactVsStd)
 	}
 
 	std::mt19937 rngf(0xF110F);
-	std::uniform_int_distribution<std::uint32_t> bf(0, ~std::uint32_t(0));
+	std::uniform_int_distribution<std::uint32_t> bf(0, std::numeric_limits<std::uint32_t>::max());
 	checked = 0;
 	while (checked < 6000)
 	{