Round sqrt correctly in the current rounding mode

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 <boost/decimal/detail/remove_trailing_zeros.hpp>
 #include <boost/decimal/detail/u256.hpp>
 #include <boost/decimal/detail/i256.hpp>
-#include <boost/decimal/numbers.hpp>
+#include <boost/decimal/detail/fenv_rounding.hpp>
 
 #ifndef BOOST_DECIMAL_BUILD_MODULE
 #include <limits>
@@ -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<int128::uint128_t>(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<int128::uint128_t>(sig_z), static_cast<int128::uint128_t>(sig_z));
 
@@ -175,14 +183,14 @@ constexpr auto sqrt128_impl(T x, int exp10val) noexcept -> T
             sig_z_sq = umul256(static_cast<int128::uint128_t>(sig_z), static_cast<int128::uint128_t>(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<T>;
-    }
 
     return z;
 }

include/boost/decimal/detail/cmath/impl/sqrt32_impl.hpp

@@ -13,18 +13,18 @@
 // 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
 // ============================================================================
 
 #include <boost/decimal/detail/cmath/impl/approx_recip_sqrt_impl.hpp>
 #include <boost/decimal/detail/cmath/frexp10.hpp>
 #include <boost/decimal/detail/remove_trailing_zeros.hpp>
-#include <boost/decimal/numbers.hpp>
+#include <boost/decimal/detail/fenv_rounding.hpp>
 
 #ifndef BOOST_DECIMAL_BUILD_MODULE
 #include <limits>
@@ -64,8 +64,16 @@ constexpr auto sqrt32_impl(T x, int exp10val) noexcept -> T
     std::uint64_t product = static_cast<std::uint64_t>(sig_gx) * r_scaled;
     std::uint32_t sig_z = static_cast<std::uint32_t>(product / scale7);
 
-    // Precompute target = sig_gx * 10^6 (avoids recomputing in Newton and rounding)
-    const std::uint64_t target = static_cast<std::uint64_t>(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<std::uint32_t>(static_cast<std::uint64_t>(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<std::uint64_t>(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<std::int64_t>(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<std::int64_t>(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<T>;
-    }
 
     return z;
 }

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 <boost/decimal/detail/cmath/impl/approx_recip_sqrt_impl.hpp>
 #include <boost/decimal/detail/cmath/frexp10.hpp>
 #include <boost/decimal/detail/remove_trailing_zeros.hpp>
-#include <boost/decimal/numbers.hpp>
+#include <boost/decimal/detail/fenv_rounding.hpp>
 
 #ifndef BOOST_DECIMAL_BUILD_MODULE
 #include <limits>
@@ -68,8 +68,16 @@ constexpr auto sqrt64_impl(T x, int exp10val) noexcept -> T
     int128::uint128_t product = static_cast<int128::uint128_t>(sig_gx) * r_scaled;
     std::uint64_t sig_z = static_cast<std::uint64_t>(product / scale16);
 
-    // Precompute target = sig_gx * 10^15 (avoids recomputing in each Newton iteration and rounding)
-    const int128::uint128_t target = static_cast<int128::uint128_t>(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<std::uint64_t>(static_cast<int128::uint128_t>(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<int128::uint128_t>(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<int128::uint128_t>(sig_z) * sig_z;
         int128::int128_t rem = static_cast<int128::int128_t>(target) - static_cast<int128::int128_t>(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<int128::int128_t>(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<int128::int128_t>(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<T>;
-    }
 
     return z;
 }

include/boost/decimal/detail/cmath/sqrt.hpp

@@ -25,7 +25,6 @@
 #include <boost/decimal/detail/config.hpp>
 #include <boost/decimal/detail/cmath/frexp10.hpp>
 #include <boost/decimal/detail/remove_trailing_zeros.hpp>
-#include <boost/decimal/numbers.hpp>
 
 // Implementation files (like SoftFloat's separate .c files)
 #include <boost/decimal/detail/cmath/impl/sqrt_lookup.hpp>
@@ -92,30 +91,18 @@ 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);
 
     if (is_pure)
     {
         const int p10 = exp10val + static_cast<int>(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<T>;
-        }
-        else if (p10_mod2 == -1)
-        {
-            result /= numbers::sqrt10_v<T>;
-        }
-        return result;
     }
 
     // ---------- Dispatch to precision-specific implementation (C++14 compatible) ----------

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"

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 ;

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<decimal_type>;
-
-          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<decimal_type>;
+          // 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<decimal_type>;
 
           result_val_p10_is_ok = (val_p10 == val_p10_ctrl);
         }