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 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..ebbb7f3a 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,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 @@ -28,7 +28,12 @@ namespace ccm::rt if constexpr (ccm::builtin::has_runtime_fmod) { return ccm::builtin::fmod_rt(x, y); } else { - return static_cast(x - (detail::trunc_scalar(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) { return static_cast(ccm::internal::fmod(static_cast(x), static_cast(y))); } + else + { + return ccm::internal::fmod(x, 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 #include namespace ccm::rt @@ -31,16 +28,11 @@ namespace ccm::rt if constexpr (ccm::builtin::has_runtime_remainder) { return ccm::builtin::remainder_rt(x, y); } else { - using FPBits_t = typename ccm::support::fp::FPBits; - 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::quiet_NaN(); - } - return static_cast(x - (detail::trunc_scalar(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(x, y, "ient); } #endif } diff --git a/include/ccmath/math/basic/fmod.hpp b/include/ccmath/math/basic/fmod.hpp index ba9f8ba7..3dc9615f 100644 --- a/include/ccmath/math/basic/fmod.hpp +++ b/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 @@ -68,7 +68,14 @@ namespace ccm } } - return static_cast(x - (ccm::trunc(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) { return internal::fmod(x, y); } + else + { + return static_cast(internal::fmod(static_cast(x), static_cast(y))); + } } template > diff --git a/include/ccmath/math/basic/impl/CMakeLists.txt b/include/ccmath/math/basic/impl/CMakeLists.txt index 44fe9fa9..80969919 100644 --- a/include/ccmath/math/basic/impl/CMakeLists.txt +++ b/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 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..76424a19 --- /dev/null +++ b/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 +#include +#include + +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 + constexpr T fmod_impl(T x, T y) noexcept + { + static_assert(std::is_floating_point_v, "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::digits <= 53, "fmod_impl only supports types up to binary64"); + using FPBits_t = support::fp::FPBits; + + 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(0), x); } + + constexpr int mant_bits = std::numeric_limits::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(y_bits.get_explicit_mantissa()); + std::uint64_t r = static_cast(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(0), x); } + return ccm::copysign(ccm::scalbn(static_cast(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 + constexpr T fmod(T x, T y) noexcept + { return impl::fmod_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; - 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::quiet_NaN(); - } - - return static_cast(x - (ccm::trunc(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(x, y, "ient); } /** 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) { 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(internal::remquo_double(static_cast(x), static_cast(y), quo)); } } 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 and 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 path in the runtime SIMD layer and the MSVC 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 path in the runtime SIMD layer and the MSVC 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)') 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::infinity()), runtime_value(1.0))); }, + [] { consume(std::fmod(runtime_value(std::numeric_limits::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, static_cast(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 6297e937..fe79df19 100644 --- a/tests/src/math/basic/fmod_test.cpp +++ b/tests/src/math/basic/fmod_test.cpp @@ -9,41 +9,179 @@ */ #include "utils/std_compare.hpp" +#include "utils/test_runtime.hpp" #include #include #include +#include +#include #include +#include -TEST(CcmathBasicTests, Fmod) +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) +{ + // 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"); + 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"); + 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 bd(0, std::numeric_limits::max()); + 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; + } - // Test that fmod works with static_assert - static_assert(ccm::fmod(1, 2) == 1, "fmod has failed testing that it is static_assert-able!"); + std::mt19937 rngf(0xF110F); + std::uniform_int_distribution bf(0, std::numeric_limits::max()); + 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, 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); + + // 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, static_cast(std::fmod)); - ccm::test::ExpectBinaryMatchesStd(10.0F, -3.0F, ccm::fmod, static_cast(std::fmod)); - ccm::test::ExpectBinaryMatchesStd(-10.0F, 3.0F, ccm::fmod, static_cast(std::fmod)); - ccm::test::ExpectBinaryMatchesStd(-10.0F, -3.0F, ccm::fmod, static_cast(std::fmod)); - ccm::test::ExpectBinaryMatchesStd(0.0F, 3.0F, ccm::fmod, static_cast(std::fmod)); +TEST(CcmathBasicTests, FmodSpecialCases) +{ + constexpr double inf = std::numeric_limits::infinity(); + constexpr double nan = std::numeric_limits::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, static_cast(std::fmod)); +TEST(CcmathBasicTests, FmodSubnormal) +{ + constexpr double dmin = std::numeric_limits::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::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, "integer fmod must return double"); - ccm::test::ExpectBinaryMatchesStd(0.0F, 1.0F, ccm::fmod, static_cast(std::fmod)); - ccm::test::ExpectBinaryMatchesStd(-0.0F, 1.0F, ccm::fmod, static_cast(std::fmod)); - ccm::test::ExpectBinaryMatchesStd(10.0F, 0.0F, ccm::fmod, static_cast(std::fmod)); - ccm::test::ExpectBinaryMatchesStd(10.0F, std::numeric_limits::infinity(), ccm::fmod, static_cast(std::fmod)); - ccm::test::ExpectBinaryMatchesStd(10.0F, -std::numeric_limits::infinity(), ccm::fmod, static_cast(std::fmod)); - ccm::test::ExpectBinaryMatchesStd(std::numeric_limits::quiet_NaN(), 10.0F, ccm::fmod, static_cast(std::fmod)); - ccm::test::ExpectBinaryMatchesStd(10.0F, std::numeric_limits::quiet_NaN(), ccm::fmod, static_cast(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/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 #include +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::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::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::infinity(); + const double nan = std::numeric_limits::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, static_cast(std::remainder)); + ccm::test::ExpectBinaryMatchesStd(-1e30, 3.0, ccm::remainder, static_cast(std::remainder)); + ccm::test::ExpectBinaryMatchesStd(1e300, 7.0, ccm::remainder, static_cast(std::remainder)); + ccm::test::ExpectBinaryMatchesStd(1e30F, 3.0F, ccm::remainder, static_cast(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, static_cast(std::remainder)); - ccm::test::ExpectBinaryMatchesStd(1.0, 0.0, ccm::remainder, static_cast(std::remainder)); - ccm::test::ExpectBinaryMatchesStd(0.0, 1.0, ccm::remainder, static_cast(std::remainder)); - ccm::test::ExpectBinaryMatchesStd(0.0, 0.0, ccm::remainder, static_cast(std::remainder)); - ccm::test::ExpectBinaryMatchesStd(-1.0, 1.0, ccm::remainder, static_cast(std::remainder)); - ccm::test::ExpectBinaryMatchesStd(1.0, -1.0, ccm::remainder, static_cast(std::remainder)); - ccm::test::ExpectBinaryMatchesStd(-1.0, -1.0, ccm::remainder, static_cast(std::remainder)); - ccm::test::ExpectBinaryMatchesStd(-1.0, 0.0, ccm::remainder, static_cast(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, static_cast(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, static_cast(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, static_cast(std::remainder)); + ccm::test::ExpectBinaryMatchesStd(subnormal_dividend, subnormal_divisor, ccm::remainder, static_cast(std::remainder)); + + // Long double goes through the same public entry point. + ccm::test::ExpectBinaryMatchesStd(5.0L, 3.0L, ccm::remainderl, static_cast(std::remainder)); + ccm::test::ExpectBinaryMatchesStd(7.5L, 2.0L, ccm::remainderl, static_cast(std::remainder)); } diff --git a/tests/src/math/basic/remquo_test.cpp b/tests/src/math/basic/remquo_test.cpp index 60a3c9ca..14ad3faf 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 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); @@ -64,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::infinity(), 2.0)) == true); // NOLINT - bool isStdLeftInfinityNegative = (std::signbit(get_std_rem(std::numeric_limits::infinity(), 2.0)) == true); // NOLINT - bool didCcmLeftInfinityReturnNan = std::isnan(get_ccm_rem(std::numeric_limits::infinity(), 2.0)); - bool didStdLeftInfinityReturnNan = std::isnan(get_std_rem(std::numeric_limits::infinity(), 2.0)); - EXPECT_EQ(isCcmLeftInfinityNegative, isStdLeftInfinityNegative); - EXPECT_EQ(didCcmLeftInfinityReturnNan, didStdLeftInfinityReturnNan); - EXPECT_EQ(get_ccm_quo(std::numeric_limits::infinity(), 2.0), get_std_quo(std::numeric_limits::infinity(), 2.0)); - - // Test with negative infinity - bool isCcmLeftNegativeInfinityNegative = (std::signbit(get_ccm_rem(-std::numeric_limits::infinity(), 2.0)) == true); // NOLINT - bool isStdLeftNegativeInfinityNegative = (std::signbit(get_std_rem(-std::numeric_limits::infinity(), 2.0)) == true); // NOLINT - bool didCcmLeftNegativeInfinityReturnNan = std::isnan(get_ccm_rem(-std::numeric_limits::infinity(), 2.0)); - bool didStdLeftNegativeInfinityReturnNan = std::isnan(get_std_rem(-std::numeric_limits::infinity(), 2.0)); - EXPECT_EQ(isCcmLeftNegativeInfinityNegative, isStdLeftNegativeInfinityNegative); - EXPECT_EQ(didCcmLeftNegativeInfinityReturnNan, didStdLeftNegativeInfinityReturnNan); - EXPECT_EQ(get_ccm_quo(-std::numeric_limits::infinity(), 2.0), get_std_quo(-std::numeric_limits::infinity(), 2.0)); - - // Test with NaN - bool isCcmLeftNanNegative = (std::signbit(get_ccm_rem(std::numeric_limits::quiet_NaN(), 2.0)) == true && - std::isnan(get_ccm_rem(std::numeric_limits::quiet_NaN(), 2.0)) == true); // NOLINT bool isStdLeftNanNegative = - (std::signbit(get_std_rem(std::numeric_limits::quiet_NaN(), 2.0)) == true && std::isnan(get_std_rem(std::numeric_limits::quiet_NaN(), 2.0)) - == true); // NOLINT bool didCcmLeftNanReturnNan = std::isnan(get_ccm_rem(std::numeric_limits::quiet_NaN(), 2.0)); bool didStdLeftNanReturnNan = - std::isnan(get_std_rem(std::numeric_limits::quiet_NaN(), 2.0)); EXPECT_EQ(isCcmLeftNanNegative, isStdLeftNanNegative); - EXPECT_EQ(didCcmLeftNanReturnNan, didStdLeftNanReturnNan); - EXPECT_EQ(get_ccm_quo(std::numeric_limits::quiet_NaN(), 2.0), get_std_quo(std::numeric_limits::quiet_NaN(), 2.0)); - - // Test with negative NaN - bool isCcmLeftNegativeNanNegative = (std::signbit(get_ccm_rem(-std::numeric_limits::quiet_NaN(), 2.0)) == true && - std::isnan(get_ccm_rem(-std::numeric_limits::quiet_NaN(), 2.0)) == true); // NOLINT bool isStdLeftNegativeNanNegative = - (std::signbit(get_std_rem(-std::numeric_limits::quiet_NaN(), 2.0)) == true && - std::isnan(get_std_rem(-std::numeric_limits::quiet_NaN(), 2.0)) == true); // NOLINT bool didCcmLeftNegativeNanReturnNan = - std::isnan(get_ccm_rem(-std::numeric_limits::quiet_NaN(), 2.0)); bool didStdLeftNegativeNanReturnNan = - std::isnan(get_std_rem(-std::numeric_limits::quiet_NaN(), 2.0)); EXPECT_EQ(isCcmLeftNegativeNanNegative, isStdLeftNegativeNanNegative); - EXPECT_EQ(didCcmLeftNegativeNanReturnNan, didStdLeftNegativeNanReturnNan); - EXPECT_EQ(get_ccm_quo(-std::numeric_limits::quiet_NaN(), 2.0), get_std_quo(-std::numeric_limits::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::infinity(); + constexpr double nan = std::numeric_limits::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))); }