diff --git a/include/boost/decimal/detail/cmath/impl/sqrt128_impl.hpp b/include/boost/decimal/detail/cmath/impl/sqrt128_impl.hpp index 0dad4f249..ba6af1d41 100644 --- a/include/boost/decimal/detail/cmath/impl/sqrt128_impl.hpp +++ b/include/boost/decimal/detail/cmath/impl/sqrt128_impl.hpp @@ -13,14 +13,14 @@ // Algorithm (inspired by SoftFloat f128_sqrt): // 1. Caller passes gx in [1, 10); get sig_gx = gx * 10^33 as u256 // 2. Use approx_recip_sqrt64 to get initial r ≈ 10^16/sqrt(gx) (~48 bits) -// 3. Compute sig_z = sig_gx * r / scale as initial sqrt approximation +// 3. Compute sig_z = sig_gx * r / scale as initial sqrt approximation, ×√10 if exp was odd // 4. Remainder-based refinement using u256 arithmetic: // rem = sig_gx * scale - sig_z² // q = rem / (2 * sig_z) (next correction) // sig_z = sig_z + q // 5. Repeat refinement to reach 34 decimal digits -// 6. Final rounding check (ensure sig_z² ≤ sig_gx * scale) -// 7. Rescale by 10^(exp/2) and ×√10 if exp was odd +// 6. Final rounding check (ensure sig_z² ≤ sig_gx * scale) in the current rounding mode +// 7. Rescale by 10^(exp/2) // // Key: ALL arithmetic uses u256/i256_sub, no floating-point rounding errors // ============================================================================ @@ -30,7 +30,7 @@ #include #include #include -#include +#include #ifndef BOOST_DECIMAL_BUILD_MODULE #include @@ -83,8 +83,16 @@ constexpr auto sqrt128_impl(T x, int exp10val) noexcept -> T // r_scaled is 64-bit; use mul128_by_64 (SoftFloat-style) instead of full umul256 u256 sig_z = mul128_by_64(gx_sig, r_scaled) / scale16; - // Precompute target = sig_gx * 10^33 (avoids recomputing in each Newton iteration) - const u256 target = umul256(gx_sig, scale33_128); + // If exp is odd, take √(10 × gx) instead so the result is rounded once (no ×√10 later) + // sig_z ≈ √gx × 10^33 × √10, then Newton corrects it + const bool odd = (exp10val & 1) != 0; + if (odd) + { + sig_z = mul128_by_64(static_cast(sig_z), 31622776601683793ULL) / scale16; + } + + // Precompute target = sig_gx * 10^33, ×10 if exp was odd (avoids recomputing in each Newton iteration) + const u256 target = umul256(gx_sig, odd ? scale33_128 * 10U : scale33_128); // ---------- Newton corrections using u256 ---------- // Newton: sig_z_new = sig_z + (sig_gx * 10^33 - sig_z²) / (2 * sig_z) @@ -159,9 +167,9 @@ constexpr auto sqrt128_impl(T x, int exp10val) noexcept -> T } } - // ---------- Final rounding (round-to-nearest) ---------- - // Find sig_z such that sig_z is the closest integer to sqrt(target) - // First ensure sig_z² ≤ target, then check if sig_z+1 is closer + // ---------- Final rounding (in the current rounding mode) ---------- + // Find sig_z such that sig_z is the correctly rounded integer of √target + // First ensure sig_z² ≤ target, then check if the mode picks sig_z+1 { u256 sig_z_sq = umul256(static_cast(sig_z), static_cast(sig_z)); @@ -175,14 +183,14 @@ constexpr auto sqrt128_impl(T x, int exp10val) noexcept -> T sig_z_sq = umul256(static_cast(sig_z), static_cast(sig_z)); } - // Step 2: Round-to-nearest check + // Step 2: Round-to-nearest check (sqrt_steps_up applies the rounding mode) // If (sig_z + 0.5)² < target, then sig_z+1 is closer // Equivalent: sig_z² + sig_z + 0.25 < target // Since we work with integers: if target - sig_z² > sig_z, round up u256 rem; i256_sub(target, sig_z_sq, rem); // rem = target - sig_z², guaranteed non-negative - if (rem > sig_z) + if (sqrt_steps_up(rem != u256{0}, rem > sig_z)) { u256 one{1}; sig_z = sig_z + one; @@ -203,16 +211,12 @@ constexpr auto sqrt128_impl(T x, int exp10val) noexcept -> T // = sig_z_hi * 10^-16 + sig_z_lo * 10^-33 T z = T{sig_z_hi, -16} + T{sig_z_lo, -33}; - // ---------- Rescale: sqrt(x) = z × 10^(e/2), ×√10 when e odd ---------- + // ---------- Rescale: sqrt(x) = z × 10^(e/2) ---------- const int half_exp = (exp10val >= 0) ? (exp10val / 2) : ((exp10val - 1) / 2); if (half_exp != 0) { z *= T{1, half_exp}; } - if ((exp10val & 1) != 0) - { - z *= numbers::sqrt10_v; - } return z; } diff --git a/include/boost/decimal/detail/cmath/impl/sqrt32_impl.hpp b/include/boost/decimal/detail/cmath/impl/sqrt32_impl.hpp index cf6dd0a34..ee47d7117 100644 --- a/include/boost/decimal/detail/cmath/impl/sqrt32_impl.hpp +++ b/include/boost/decimal/detail/cmath/impl/sqrt32_impl.hpp @@ -13,10 +13,10 @@ // Algorithm: // 1. Caller passes gx in [1, 10); get sig_gx = gx * 10^6 as integer // 2. Call approx_recip_sqrt32 → r_scaled ≈ 10^7 / sqrt(gx) (integer, ~24 bits) -// 3. Compute sig_z = sig_gx * r_scaled / 10^7 ≈ sqrt(gx) * 10^6 +// 3. Compute sig_z = sig_gx * r_scaled / 10^7 ≈ sqrt(gx) * 10^6, ×√10 if exp was odd // 4. Newton correction using exact integer remainder -// 5. Final rounding check (integer) -// 6. Rescale by 10^(exp/2) and ×√10 if exp was odd +// 5. Final rounding check (integer) in the current rounding mode +// 6. Rescale by 10^(exp/2) // // Key improvement: ALL arithmetic is integer, no floating-point until final result // ============================================================================ @@ -24,7 +24,7 @@ #include #include #include -#include +#include #ifndef BOOST_DECIMAL_BUILD_MODULE #include @@ -64,8 +64,16 @@ constexpr auto sqrt32_impl(T x, int exp10val) noexcept -> T std::uint64_t product = static_cast(sig_gx) * r_scaled; std::uint32_t sig_z = static_cast(product / scale7); - // Precompute target = sig_gx * 10^6 (avoids recomputing in Newton and rounding) - const std::uint64_t target = static_cast(sig_gx) * scale6; + // If exp is odd, take √(10 × gx) instead so the result is rounded once (no ×√10 later) + // sig_z ≈ √gx × 10^6 × √10, then Newton corrects it + const bool odd = (exp10val & 1) != 0; + if (odd) + { + sig_z = static_cast(static_cast(sig_z) * 31622777ULL / scale7); + } + + // Precompute target = sig_gx * 10^6, ×10 if exp was odd (avoids recomputing in Newton and rounding) + const std::uint64_t target = static_cast(sig_gx) * (odd ? scale7 : scale6); // ---------- Newton correction with exact integer remainder ---------- // rem = target - sig_z² (exact integer) @@ -84,25 +92,27 @@ constexpr auto sqrt32_impl(T x, int exp10val) noexcept -> T } // ---------- Final rounding check ---------- - // Ensure z² ≤ gx (z is a lower bound) - if (rem < 0) + // Ensure z² ≤ gx (z is a lower bound); Newton can land a few steps above it + while (rem < 0) { --sig_z; + rem += 2 * static_cast(sig_z) + 1; + } + // Round up if the mode asks for it: target - sig_z² > sig_z means (sig_z + 0.5)² < target + if (sqrt_steps_up(rem != 0, rem > static_cast(sig_z))) + { + ++sig_z; } // Convert back to decimal type T z{sig_z, -6}; // sig_z * 10^-6 - // ---------- Rescale: sqrt(x) = z × 10^(e/2), ×√10 when e odd ---------- + // ---------- Rescale: sqrt(x) = z × 10^(e/2) ---------- const int half_exp = (exp10val >= 0) ? (exp10val / 2) : ((exp10val - 1) / 2); if (half_exp != 0) { z *= T{1, half_exp}; } - if ((exp10val & 1) != 0) - { - z *= numbers::sqrt10_v; - } return z; } diff --git a/include/boost/decimal/detail/cmath/impl/sqrt64_impl.hpp b/include/boost/decimal/detail/cmath/impl/sqrt64_impl.hpp index 72953cc36..ab83a4f39 100644 --- a/include/boost/decimal/detail/cmath/impl/sqrt64_impl.hpp +++ b/include/boost/decimal/detail/cmath/impl/sqrt64_impl.hpp @@ -13,10 +13,10 @@ // Algorithm: // 1. Caller passes gx in [1, 10); get sig_gx = gx * 10^15 as integer // 2. Call approx_recip_sqrt64 → r_scaled ≈ 10^16 / sqrt(gx) (integer, ~48 bits) -// 3. Compute sig_z = sig_gx * r_scaled / 10^16 ≈ sqrt(gx) * 10^15 +// 3. Compute sig_z = sig_gx * r_scaled / 10^16 ≈ sqrt(gx) * 10^15, ×√10 if exp was odd // 4. Newton corrections using exact integer remainder (needs 128-bit) -// 5. Final rounding check (integer) -// 6. Rescale by 10^(exp/2) and ×√10 if exp was odd +// 5. Final rounding check (integer) in the current rounding mode +// 6. Rescale by 10^(exp/2) // // Key improvement: ALL arithmetic is integer, no floating-point until final result // Requires: 128-bit integer support (__uint128_t or equivalent) for full precision @@ -27,7 +27,7 @@ #include #include #include -#include +#include #ifndef BOOST_DECIMAL_BUILD_MODULE #include @@ -68,8 +68,16 @@ constexpr auto sqrt64_impl(T x, int exp10val) noexcept -> T int128::uint128_t product = static_cast(sig_gx) * r_scaled; std::uint64_t sig_z = static_cast(product / scale16); - // Precompute target = sig_gx * 10^15 (avoids recomputing in each Newton iteration and rounding) - const int128::uint128_t target = static_cast(sig_gx) * scale15; + // If exp is odd, take √(10 × gx) instead so the result is rounded once (no ×√10 later) + // sig_z ≈ √gx × 10^15 × √10, then Newton corrects it + const bool odd = (exp10val & 1) != 0; + if (odd) + { + sig_z = static_cast(static_cast(sig_z) * 31622776601683793ULL / scale16); + } + + // Precompute target = sig_gx * 10^15, ×10 if exp was odd (avoids recomputing in each Newton iteration and rounding) + const int128::uint128_t target = static_cast(sig_gx) * (odd ? scale16 : scale15); // ---------- Newton correction with exact integer remainder ---------- // rem = target - sig_z² (exact 128-bit integer) @@ -96,30 +104,33 @@ constexpr auto sqrt64_impl(T x, int exp10val) noexcept -> T } } - // Final rounding check + // Final rounding check (in the current rounding mode) { int128::uint128_t z_squared = static_cast(sig_z) * sig_z; int128::int128_t rem = static_cast(target) - static_cast(z_squared); - if (rem < 0) + // Ensure z² ≤ gx (z is a lower bound); Newton can land a few steps above it + while (rem < 0) { --sig_z; + rem += static_cast(2 * sig_z + 1); + } + // Round up if the mode asks for it: target - sig_z² > sig_z means (sig_z + 0.5)² < target + if (sqrt_steps_up(rem != 0, rem > static_cast(sig_z))) + { + ++sig_z; } } // Convert back to decimal type T z{sig_z, -15}; // sig_z * 10^-15 - // ---------- Rescale: sqrt(x) = z × 10^(e/2), ×√10 when e odd ---------- + // ---------- Rescale: sqrt(x) = z × 10^(e/2) ---------- const int half_exp = (exp10val >= 0) ? (exp10val / 2) : ((exp10val - 1) / 2); if (half_exp != 0) { z *= T{1, half_exp}; } - if ((exp10val & 1) != 0) - { - z *= numbers::sqrt10_v; - } return z; } diff --git a/include/boost/decimal/detail/cmath/sqrt.hpp b/include/boost/decimal/detail/cmath/sqrt.hpp index 5865105a0..96f814989 100644 --- a/include/boost/decimal/detail/cmath/sqrt.hpp +++ b/include/boost/decimal/detail/cmath/sqrt.hpp @@ -25,7 +25,6 @@ #include #include #include -#include // Implementation files (like SoftFloat's separate .c files) #include @@ -92,6 +91,7 @@ constexpr auto sqrt_impl(T x) noexcept auto sig = frexp10(x, &exp10val); // ---------- Fast path: pure powers of 10 ---------- + // Only even powers return here; an odd power goes through the kernel so √10 is rounded in the current mode const auto zeros_removal = remove_trailing_zeros(sig); const bool is_pure = (zeros_removal.trimmed_number == 1U); @@ -99,23 +99,10 @@ constexpr auto sqrt_impl(T x) noexcept { const int p10 = exp10val + static_cast(zeros_removal.number_of_removed_zeros); - if (p10 == 0) + if ((p10 % 2) == 0) { - return T{1}; + return T{1, p10 / 2}; } - - const int p10_mod2 = (p10 % 2); - T result = T{1, p10 / 2}; - - if (p10_mod2 == 1) - { - result *= numbers::sqrt10_v; - } - else if (p10_mod2 == -1) - { - result /= numbers::sqrt10_v; - } - return result; } // ---------- Dispatch to precision-specific implementation (C++14 compatible) ---------- diff --git a/include/boost/decimal/detail/fenv_rounding.hpp b/include/boost/decimal/detail/fenv_rounding.hpp index 8f36166f2..e98473564 100644 --- a/include/boost/decimal/detail/fenv_rounding.hpp +++ b/include/boost/decimal/detail/fenv_rounding.hpp @@ -660,6 +660,28 @@ BOOST_DECIMAL_CUDA_CONSTEXPR auto overflow_is_finite(const bool is_negative) noe round == (is_negative ? rounding_mode::fe_dec_upward : rounding_mode::fe_dec_downward); } +// IEEE 754-2019 4.3: the floor of an inexact root steps up in the upward mode, past the half +// point in the nearest modes, and never in the modes toward zero. A root is never a tie. +BOOST_DECIMAL_CUDA_CONSTEXPR auto sqrt_steps_up(const bool inexact, const bool past_half) noexcept -> bool +{ + auto round {_boost_decimal_global_rounding_mode}; + #ifndef BOOST_DECIMAL_NO_CONSTEVAL_DETECTION + if (!BOOST_DECIMAL_IS_CONSTANT_EVALUATED(inexact)) + { + round = _boost_decimal_global_runtime_rounding_mode; + } + #endif + if (round == rounding_mode::fe_dec_upward) + { + return inexact; + } + if (round == rounding_mode::fe_dec_toward_zero || round == rounding_mode::fe_dec_downward) + { + return false; + } + return past_half; +} + #if defined(__clang__) # pragma clang diagnostic push # pragma clang diagnostic ignored "-Wsign-conversion" diff --git a/test/Jamfile b/test/Jamfile index 589b2c10e..fc3d39abb 100644 --- a/test/Jamfile +++ b/test/Jamfile @@ -219,6 +219,8 @@ run test_sin_cos.cpp ; run test_sinh.cpp ; run test_snprintf.cpp ; run test_sqrt.cpp ; +run test_sqrt_rounding.cpp ; +run test_sqrt_rounding_upward.cpp ; run test_string_construction.cpp ; run test_string_locale_conversion.cpp ; run test_strtod.cpp ; diff --git a/test/test_sqrt.cpp b/test/test_sqrt.cpp index 60aa7e5ad..877d85139 100644 --- a/test/test_sqrt.cpp +++ b/test/test_sqrt.cpp @@ -163,15 +163,10 @@ namespace local { result_val_p10_is_ok = (val_p10 == decimal_type { 1, np / 2 }); } - else if(np_mod2 == -1) + else { - decimal_type val_p10_ctrl = decimal_type { 1, np / 2 } / boost::decimal::numbers::sqrt10_v; - - result_val_p10_is_ok = (val_p10 == val_p10_ctrl); - } - else if(np_mod2 == 1) - { - decimal_type val_p10_ctrl = decimal_type { 1, np / 2 } * boost::decimal::numbers::sqrt10_v; + // sqrt(10^np) is sqrt(10) * 10^((np - 1) / 2), and the scale of the nearest sqrt(10) is exact + decimal_type val_p10_ctrl = decimal_type { 1, (np - 1) / 2 } * boost::decimal::numbers::sqrt10_v; result_val_p10_is_ok = (val_p10 == val_p10_ctrl); } diff --git a/test/test_sqrt_rounding.cpp b/test/test_sqrt_rounding.cpp new file mode 100644 index 000000000..51398475a --- /dev/null +++ b/test/test_sqrt_rounding.cpp @@ -0,0 +1,167 @@ +// Copyright 2026 Matt Borland +// Distributed under the Boost Software License, Version 1.0. +// https://www.boost.org/LICENSE_1_0.txt +// +// IEEE 754-2019 4.3 and 5.4.1: sqrt is correctly rounded in the current mode. The 32 and 64 +// bit kernels truncated the root, an odd exponent multiplied by a rounded sqrt(10), and no +// kernel read the mode. + +#include +#include +#include + +using namespace boost::decimal; +using namespace boost::decimal::literals; + +template struct fast_of; +template <> struct fast_of { using type = decimal_fast32_t; }; +template <> struct fast_of { using type = decimal_fast64_t; }; +template <> struct fast_of { using type = decimal_fast128_t; }; + +// The value of the mode from the floor and the ceiling of the root. A root is never a tie. +template +constexpr auto expected(const rounding_mode mode, const T floor, const T ceil, const bool near_up) -> T +{ + return mode == rounding_mode::fe_dec_upward ? ceil : + mode == rounding_mode::fe_dec_toward_zero || mode == rounding_mode::fe_dec_downward ? floor : + near_up ? ceil : floor; +} + +template +void check_mode(const T x, const rounding_mode mode, const T value) +{ + using F = typename fast_of::type; + fesetround(mode); + BOOST_TEST_EQ(sqrt(x), value); + BOOST_TEST_EQ(sqrt(static_cast(x)), static_cast(value)); +} + +// floor and ceil come from an integer square root at the precision of the type +template +void check(const T x, const T floor, const T ceil, const bool near_up) +{ + #ifdef BOOST_DECIMAL_NO_CONSTEVAL_DETECTION + // fesetround has no effect here, thus only the compile-time mode is tested + check_mode(x, _boost_decimal_global_rounding_mode, expected(_boost_decimal_global_rounding_mode, floor, ceil, near_up)); + #else + check_mode(x, rounding_mode::fe_dec_to_nearest, expected(rounding_mode::fe_dec_to_nearest, floor, ceil, near_up)); + check_mode(x, rounding_mode::fe_dec_to_nearest_from_zero, expected(rounding_mode::fe_dec_to_nearest_from_zero, floor, ceil, near_up)); + check_mode(x, rounding_mode::fe_dec_toward_zero, floor); + check_mode(x, rounding_mode::fe_dec_downward, floor); + check_mode(x, rounding_mode::fe_dec_upward, ceil); + #endif +} + +template +void check(const int n, const T floor, const T ceil, const bool near_up) +{ + check(T {n, 0}, floor, ceil, near_up); +} + +int main() +{ + constexpr auto mode {_boost_decimal_global_rounding_mode}; + static_assert(sqrt(decimal32_t {11, 0}) == expected(mode, decimal32_t {3316624, -6}, decimal32_t {3316625, -6}, true), ""); + static_assert(sqrt(decimal64_t {26, 0}) == expected(mode, decimal64_t {UINT64_C(5099019513592784), -15}, decimal64_t {UINT64_C(5099019513592785), -15}, true), ""); + + // A failure prints every digit of the two values + std::cerr.precision(std::numeric_limits::digits10); + + check(2, 1.414213_df, 1.414214_df, true); + check(3, 1.732050_df, 1.732051_df, true); + check(4, 2_df, 2_df, false); + check(5, 2.236067_df, 2.236068_df, true); + check(6, 2.449489_df, 2.449490_df, true); + check(7, 2.645751_df, 2.645752_df, false); + check(8, 2.828427_df, 2.828428_df, false); + check(9, 3_df, 3_df, false); + check(10, 3.162277_df, 3.162278_df, true); + check(11, 3.316624_df, 3.316625_df, true); + check(12, 3.464101_df, 3.464102_df, true); + check(13, 3.605551_df, 3.605552_df, false); + check(15, 3.872983_df, 3.872984_df, false); + check(17, 4.123105_df, 4.123106_df, true); + check(19, 4.358898_df, 4.358899_df, true); + check(26, 5.099019_df, 5.099020_df, true); + check(29, 5.385164_df, 5.385165_df, true); + check(48, 6.928203_df, 6.928204_df, false); + check(50, 7.071067_df, 7.071068_df, true); + check(61, 7.810249_df, 7.810250_df, true); + check(73, 8.544003_df, 8.544004_df, true); + check(83, 9.110433_df, 9.110434_df, true); + check(97, 9.848857_df, 9.848858_df, true); + check(99, 9.949874_df, 9.949875_df, false); + check(100, 10_df, 10_df, false); + check(101, 10.04987_df, 10.04988_df, true); + check(500, 22.36067_df, 22.36068_df, true); + check(999, 31.60696_df, 31.60697_df, false); + check(1000, 31.62277_df, 31.62278_df, true); + // The Newton step lands two above the floor here, and one step down was not enough + check(6.028760e-17_df, 7.764508e-9_df, 7.764509e-9_df, true); + check(4.472279e-67_df, 6.687509e-34_df, 6.687510e-34_df, true); + check(5.428215e+70_df, 2.329852e+35_df, 2.329853e+35_df, true); + + check(2, 1.414213562373095_dd, 1.414213562373096_dd, false); + check(3, 1.732050807568877_dd, 1.732050807568878_dd, false); + check(4, 2_dd, 2_dd, false); + check(5, 2.236067977499789_dd, 2.236067977499790_dd, true); + check(6, 2.449489742783178_dd, 2.449489742783179_dd, false); + check(7, 2.645751311064590_dd, 2.645751311064591_dd, true); + check(8, 2.828427124746190_dd, 2.828427124746191_dd, false); + check(9, 3_dd, 3_dd, false); + check(10, 3.162277660168379_dd, 3.162277660168380_dd, false); + check(11, 3.316624790355399_dd, 3.316624790355400_dd, true); + check(12, 3.464101615137754_dd, 3.464101615137755_dd, true); + check(13, 3.605551275463989_dd, 3.605551275463990_dd, false); + check(15, 3.872983346207416_dd, 3.872983346207417_dd, true); + check(17, 4.123105625617660_dd, 4.123105625617661_dd, true); + check(19, 4.358898943540673_dd, 4.358898943540674_dd, true); + check(26, 5.099019513592784_dd, 5.099019513592785_dd, true); + check(29, 5.385164807134504_dd, 5.385164807134505_dd, false); + check(48, 6.928203230275509_dd, 6.928203230275510_dd, false); + check(50, 7.071067811865475_dd, 7.071067811865476_dd, false); + check(61, 7.810249675906654_dd, 7.810249675906655_dd, false); + check(73, 8.544003745317531_dd, 8.544003745317532_dd, false); + check(83, 9.110433579144298_dd, 9.110433579144299_dd, true); + check(97, 9.848857801796104_dd, 9.848857801796105_dd, true); + check(99, 9.949874371066199_dd, 9.949874371066200_dd, true); + check(100, 10_dd, 10_dd, false); + check(101, 10.04987562112089_dd, 10.04987562112090_dd, false); + check(500, 22.36067977499789_dd, 22.36067977499790_dd, true); + check(999, 31.60696125855821_dd, 31.60696125855822_dd, true); + check(1000, 31.62277660168379_dd, 31.62277660168380_dd, false); + + check(2, 1.414213562373095048801688724209698_dl, 1.414213562373095048801688724209699_dl, false); + check(3, 1.732050807568877293527446341505872_dl, 1.732050807568877293527446341505873_dl, false); + check(4, 2_dl, 2_dl, false); + check(5, 2.236067977499789696409173668731276_dl, 2.236067977499789696409173668731277_dl, false); + check(6, 2.449489742783178098197284074705891_dl, 2.449489742783178098197284074705892_dl, false); + check(7, 2.645751311064590590501615753639260_dl, 2.645751311064590590501615753639261_dl, false); + check(8, 2.828427124746190097603377448419396_dl, 2.828427124746190097603377448419397_dl, false); + check(9, 3_dl, 3_dl, false); + check(10, 3.162277660168379331998893544432718_dl, 3.162277660168379331998893544432719_dl, true); + check(11, 3.316624790355399849114932736670686_dl, 3.316624790355399849114932736670687_dl, true); + check(12, 3.464101615137754587054892683011744_dl, 3.464101615137754587054892683011745_dl, true); + check(13, 3.605551275463989293119221267470495_dl, 3.605551275463989293119221267470496_dl, true); + check(15, 3.872983346207416885179265399782399_dl, 3.872983346207416885179265399782400_dl, true); + check(17, 4.123105625617660549821409855974077_dl, 4.123105625617660549821409855974078_dl, false); + check(19, 4.358898943540673552236981983859615_dl, 4.358898943540673552236981983859616_dl, true); + check(26, 5.099019513592784830028224109022781_dl, 5.099019513592784830028224109022782_dl, true); + check(29, 5.385164807134504031250710491540329_dl, 5.385164807134504031250710491540330_dl, true); + check(48, 6.928203230275509174109785366023489_dl, 6.928203230275509174109785366023490_dl, false); + check(50, 7.071067811865475244008443621048490_dl, 7.071067811865475244008443621048491_dl, false); + check(61, 7.810249675906654394129722735759101_dl, 7.810249675906654394129722735759102_dl, false); + check(73, 8.544003745317531167871648326239706_dl, 8.544003745317531167871648326239707_dl, false); + check(83, 9.110433579144298881945626104688669_dl, 9.110433579144298881945626104688670_dl, false); + check(97, 9.848857801796104721746211414917624_dl, 9.848857801796104721746211414917625_dl, false); + check(99, 9.949874371066199547344798210012060_dl, 9.949874371066199547344798210012061_dl, false); + check(100, 10_dl, 10_dl, false); + check(101, 10.04987562112089027021926491275957_dl, 10.04987562112089027021926491275958_dl, true); + check(500, 22.36067977499789696409173668731276_dl, 22.36067977499789696409173668731277_dl, false); + check(999, 31.60696125855821654520421398569900_dl, 31.60696125855821654520421398569901_dl, false); + check(1000, 31.62277660168379331998893544432718_dl, 31.62277660168379331998893544432719_dl, true); + + fesetround(rounding_mode::fe_dec_to_nearest); + + return boost::report_errors(); +} diff --git a/test/test_sqrt_rounding_upward.cpp b/test/test_sqrt_rounding_upward.cpp new file mode 100644 index 000000000..0dc7f2e58 --- /dev/null +++ b/test/test_sqrt_rounding_upward.cpp @@ -0,0 +1,11 @@ +// Copyright 2026 Matt Borland +// Distributed under the Boost Software License, Version 1.0. +// https://www.boost.org/LICENSE_1_0.txt +// +// The same checks under the compile-time upward mode, and without constant evaluation +// detection, where only that mode exists + +#define BOOST_DECIMAL_FE_DEC_UPWARD +#define BOOST_DECIMAL_NO_CONSTEVAL_DETECTION + +#include "test_sqrt_rounding.cpp"