Fix constexpr remainder and make fmod exact for large quotients

.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

include/ccmath/internal/math/runtime/func/basic/fmod_rt.hpp

@@ -12,8 +12,8 @@
 
 #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_impl.hpp"
 
 #include <type_traits>
 
@@ -28,7 +28,12 @@ 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 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 ccm::internal::fmod(x, y);
+			}
 		}
 #endif
 	}

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
 	}

include/ccmath/math/basic/fmod.hpp

@@ -15,7 +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/nearest/trunc.hpp"
+#include "ccmath/math/basic/impl/fmod_impl.hpp"
 
 #include <limits>
 
@@ -68,7 +68,14 @@ namespace ccm
 				}
 			}
 
-			return static_cast<T>(x - (ccm::trunc<T>(x / y) * 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(static_cast<double>(x), static_cast<double>(y)));
+			}
 		}
 
 		template <typename T, typename U, typename TC = std::common_type_t<T, U>>

include/ccmath/math/basic/impl/CMakeLists.txt

@@ -1,4 +1,5 @@
 ccm_add_headers(
+        fmod_impl.hpp
         nan_double_impl.hpp
         nan_float_impl.hpp
         nan_ldouble_impl.hpp

include/ccmath/math/basic/impl/fmod_impl.hpp

@@ -0,0 +1,87 @@
+/*
+ * 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");
+			// 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);
+			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

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);
 	}
 
 	/**

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));
 		}
 	}
 

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)')

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);

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.