Compute frexp with integer arithmetic

include/boost/decimal/detail/countl.hpp

@@ -102,7 +102,9 @@ BOOST_DECIMAL_CUDA_CONSTEXPR auto bit_scan_reverse(std::uint64_t bb) noexcept ->
 template <typename T>
 BOOST_DECIMAL_CUDA_CONSTEXPR int countl_impl(T x) noexcept
 {
-    return x ? bit_scan_reverse(static_cast<std::uint64_t>(x)) ^ 63 : std::numeric_limits<T>::digits;
+    // bit_scan_reverse ^ 63 counts the leading zeros in 64 bits, which is 64 - digits more than T has
+    return x ? (bit_scan_reverse(static_cast<std::uint64_t>(x)) ^ 63) - (64 - std::numeric_limits<T>::digits)
+             : std::numeric_limits<T>::digits;
 }
 
 #endif

include/boost/decimal/detail/u256.hpp

@@ -939,6 +939,18 @@ BOOST_DECIMAL_CUDA_CONSTEXPR u256 operator*(const UnsignedInteger lhs, const u25
 }
 BOOST_DECIMAL_CUDA_CONSTEXPR u256 mul128_by_64(const int128::uint128_t& a, const std::uint64_t b) noexcept;
 
+// Returns the high 128 bits of a uint128 * uint128 -> u256 product
+BOOST_DECIMAL_CUDA_CONSTEXPR int128::uint128_t umul256_hi(const int128::uint128_t& a, const int128::uint128_t& b) noexcept
+{
+    const int128::uint128_t ll {int128::uint128_t{a.low} * b.low};
+    const int128::uint128_t lh {int128::uint128_t{a.low} * b.high};
+    const int128::uint128_t hl {int128::uint128_t{a.high} * b.low};
+    const int128::uint128_t hh {int128::uint128_t{a.high} * b.high};
+    // All sums are of two uint128_t, because GCC adds a uint64_t to a uint128_t with a compare and a branch.
+    const int128::uint128_t mid {int128::uint128_t{ll.high} + int128::uint128_t{lh.low} + int128::uint128_t{hl.low}};
+    return hh + int128::uint128_t{lh.high} + int128::uint128_t{hl.high} + int128::uint128_t{mid.high};
+}
+
 BOOST_DECIMAL_CUDA_CONSTEXPR u256 umul256(const int128::uint128_t& a, const int128::uint128_t& b) noexcept
 {
     if (BOOST_DECIMAL_UNLIKELY(b.high == 0U))
@@ -983,44 +995,21 @@ BOOST_DECIMAL_CUDA_CONSTEXPR u256 umul256(const int128::uint128_t& a, const int1
 // Returns the high 256 bits of a u256 * u256 -> u512 product
 BOOST_DECIMAL_CUDA_CONSTEXPR u256 umul512_hi(const u256& a, const u256& b) noexcept
 {
-    // Decompose each operand into two uint128 halves.
-    const int128::uint128_t a_lo {a.bytes[1], a.bytes[0]};
-    const int128::uint128_t a_hi {a.bytes[3], a.bytes[2]};
-    const int128::uint128_t b_lo {b.bytes[1], b.bytes[0]};
-    const int128::uint128_t b_hi {b.bytes[3], b.bytes[2]};
-
-    // Four uint128 * uint128 -> u256 partial products.
-    const u256 p_ll {umul256(a_lo, b_lo)};
-    const u256 p_lh {umul256(a_lo, b_hi)};
-    const u256 p_hl {umul256(a_hi, b_lo)};
-    const u256 p_hh {umul256(a_hi, b_hi)};
-
-    const int128::uint128_t p_ll_hi {p_ll.bytes[3], p_ll.bytes[2]};
-    const int128::uint128_t p_lh_lo {p_lh.bytes[1], p_lh.bytes[0]};
-    const int128::uint128_t p_lh_hi {p_lh.bytes[3], p_lh.bytes[2]};
-    const int128::uint128_t p_hl_lo {p_hl.bytes[1], p_hl.bytes[0]};
-    const int128::uint128_t p_hl_hi {p_hl.bytes[3], p_hl.bytes[2]};
-    const int128::uint128_t p_hh_lo {p_hh.bytes[1], p_hh.bytes[0]};
-    const int128::uint128_t p_hh_hi {p_hh.bytes[3], p_hh.bytes[2]};
-
-    int128::uint128_t w1 {p_ll_hi};
-    w1 += p_lh_lo;
-    std::uint64_t carry_w1 {(w1 < p_lh_lo) ? UINT64_C(1) : UINT64_C(0)};
-    w1 += p_hl_lo;
-    carry_w1 += (w1 < p_hl_lo) ? UINT64_C(1) : UINT64_C(0);
-
-    int128::uint128_t w2 {p_lh_hi};
-    w2 += p_hl_hi;
-    std::uint64_t carry_w2 {(w2 < p_hl_hi) ? UINT64_C(1) : UINT64_C(0)};
-    w2 += p_hh_lo;
-    carry_w2 += (w2 < p_hh_lo) ? UINT64_C(1) : UINT64_C(0);
-    const int128::uint128_t w2_before_carry {w2};
-    w2 += int128::uint128_t{0, carry_w1};
-    carry_w2 += (w2 < w2_before_carry) ? UINT64_C(1) : UINT64_C(0);
-
-    const int128::uint128_t w3 {p_hh_hi + int128::uint128_t{0, carry_w2}};
-
-    return u256{w3, w2};
+    // Schoolbook product of 64-bit limbs. It is faster than four umul256 products,
+    // because umul256 adds branches and carry compares.
+    std::uint64_t r[8] {};
+    for (int i {0}; i < 4; ++i)
+    {
+        std::uint64_t carry {0U};
+        for (int j {0}; j < 4; ++j)
+        {
+            const int128::uint128_t t {int128::uint128_t{a.bytes[i]} * b.bytes[j] + r[i + j] + carry};
+            r[i + j] = t.low;
+            carry = t.high;
+        }
+        r[i + 4] = carry;
+    }
+    return u256{r[7], r[6], r[5], r[4]};
 }
 
 // 128x64 -> 256 multiplication (SoftFloat-style lightweight primitive)

test/Jamfile

@@ -145,6 +145,9 @@ run github_issue_1473.cpp ;
 run github_issue_1476.cpp ;
 run github_issue_1482.cpp ;
 run github_issue_1485.cpp ;
+# The test sets the evaluation method before the include, and the pch would hide it
+run github_issue_1487_eval_method.cpp
+  : : : <pch>off ;
 
 run link_1.cpp link_2.cpp link_3.cpp ;
 run quick.cpp ;

test/github_issue_1487_eval_method.cpp

@@ -0,0 +1,37 @@
+// Copyright 2026 Shen-Ta Hsieh
+// Distributed under the Boost Software License, Version 1.0.
+// https://www.boost.org/LICENSE_1_0.txt
+//
+// https://github.com/boostorg/decimal/issues/1487
+//
+// frexp with the widest evaluation type. Each value is just below a power of two, so the fraction
+// in decimal128_t is just below 1, and it rounds to 1 in the type of the value.
+
+#define BOOST_DECIMAL_DEC_EVAL_METHOD 2
+
+#include <boost/decimal.hpp>
+#include <boost/core/lightweight_test.hpp>
+
+using namespace boost::decimal;
+using namespace boost::decimal::literals;
+
+template <typename T>
+void test(const T v, const T frac, const int expon)
+{
+    int e {};
+    BOOST_TEST_EQ(frexp(v, &e), frac);
+    BOOST_TEST_EQ(e, expon);
+}
+
+int main()
+{
+    // 9223372e12 = 2^63 * 0.999999996, 9536743e-13 = 2^-20 * 0.99999998,
+    // and 4722366482869645e6 = 2^72 * 0.99999999999999995
+    test(9223372e12_DF, 0.5_DF, 64);
+    test(9536743e-13_DF, 0.5_DF, -19);
+    test(static_cast<decimal_fast32_t>(9223372e12_DF), static_cast<decimal_fast32_t>(0.5_DF), 64);
+    test(4722366482869645e6_DD, 0.5_DD, 73);
+    test(static_cast<decimal_fast64_t>(4722366482869645e6_DD), static_cast<decimal_fast64_t>(0.5_DD), 73);
+
+    return boost::report_errors();
+}

test/test_frexp_ldexp.cpp

@@ -260,6 +260,23 @@ namespace local
       BOOST_TEST(result_is_ok);
     }
 
+    // frexp keeps the sign of a zero.
+    frexp_dec = frexp(-zero, &n_dec);
+    BOOST_TEST_EQ(frexp_dec, zero);
+    BOOST_TEST(signbit(frexp_dec));
+    BOOST_TEST_EQ(n_dec, 0);
+
+    // The narrow pass of decimal64 cannot decide these fractions, and the wide pass rounds them.
+    {
+      using namespace boost::decimal::literals;
+
+      auto n_dd = int { };
+      BOOST_TEST_EQ(frexp(262e200_DD, &n_dd), 0.6685196997600160_DD);
+      BOOST_TEST_EQ(n_dd, 673);
+      BOOST_TEST_EQ(frexp(110610e-22_DD, &n_dd), 0.7970290476535209_DD);
+      BOOST_TEST_EQ(n_dd, -56);
+    }
+
     for(auto index = static_cast<unsigned>(UINT8_C(0)); index < static_cast<unsigned>(UINT8_C(4)); ++index)
     {
       static_cast<void>(index);