From eed74046737673f8f4ca2aea9d21241e79efa2ed Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Fri, 15 May 2026 15:11:55 -0400 Subject: [PATCH 1/5] Implement isqrt --- include/boost/int128/utilities.hpp | 35 ++++++++++++++++++++++++++++++ 1 file changed, 35 insertions(+) diff --git a/include/boost/int128/utilities.hpp b/include/boost/int128/utilities.hpp index 2be954c8..0464b548 100644 --- a/include/boost/int128/utilities.hpp +++ b/include/boost/int128/utilities.hpp @@ -218,6 +218,41 @@ BOOST_INT128_EXPORT BOOST_INT128_HOST_DEVICE constexpr int128_t ipow(int128_t ba return result; } +// Integer square root: returns floor(sqrt(n)). +BOOST_INT128_EXPORT BOOST_INT128_HOST_DEVICE constexpr uint128_t isqrt(const uint128_t n) noexcept +{ + if (n < 2U) + { + return n; + } + + // 2^ceil(bit_width(n)/2) is the smallest power of two whose square exceeds n. + uint128_t x {uint128_t{1} << ((bit_width(n) + 1) / 2)}; + + while (true) + { + const uint128_t y {(x + n / x) >> 1}; + + if (y >= x) + { + return x; + } + + x = y; + } +} + +// Signed overload. Negative inputs are documented to return 0. +BOOST_INT128_EXPORT BOOST_INT128_HOST_DEVICE constexpr int128_t isqrt(const int128_t n) noexcept +{ + if (BOOST_INT128_UNLIKELY(n < 0)) + { + return int128_t{0}; + } + + return static_cast(isqrt(static_cast(n))); +} + } // namespace int128 } // namespace boost From 717d89e252cee2662e50899f3a48504d44777e19 Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Fri, 15 May 2026 15:12:01 -0400 Subject: [PATCH 2/5] Test isqrt --- test/Jamfile | 1 + test/test_isqrt.cpp | 258 ++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 259 insertions(+) create mode 100644 test/test_isqrt.cpp diff --git a/test/Jamfile b/test/Jamfile index e8b163c1..25e598d7 100644 --- a/test/Jamfile +++ b/test/Jamfile @@ -81,6 +81,7 @@ run test_gcd_lcm.cpp ; run test_midpoint.cpp ; run test_powm.cpp ; run test_ipow.cpp ; +run test_isqrt.cpp ; run test_format.cpp ; run test_fmt_format.cpp ; diff --git a/test/test_isqrt.cpp b/test/test_isqrt.cpp new file mode 100644 index 00000000..369f996c --- /dev/null +++ b/test/test_isqrt.cpp @@ -0,0 +1,258 @@ +// Copyright 2026 Matt Borland +// Distributed under the Boost Software License, Version 1.0. +// https://www.boost.org/LICENSE_1_0.txt + +#ifndef BOOST_INT128_BUILD_MODULE + +#include + +#else + +import boost.int128; + +#endif + +#include +#include +#include + +using namespace boost::int128; + +namespace { + +// Naive bit-by-bit reference. Independent of the Newton implementation under +// test, so it serves as a ground truth for cross-checking. +constexpr uint128_t isqrt_ref(uint128_t n) noexcept +{ + if (n < 2U) + { + return n; + } + + uint128_t res {0}; + uint128_t bit {uint128_t{1} << 126}; + + while (bit > n) + { + bit >>= 2; + } + + while (bit != 0U) + { + if (n >= res + bit) + { + n -= res + bit; + res = (res >> 1) + bit; + } + else + { + res >>= 1; + } + + bit >>= 2; + } + + return res; +} + +// Verify the defining invariant: r == isqrt(n) iff r*r <= n < (r+1)^2. +// Computed with overflow-safe comparisons so callers can pass values near the +// uint128 upper bound. +void check_invariant(const uint128_t n) +{ + const uint128_t r {isqrt(n)}; + + BOOST_TEST(r * r <= n); + + const uint128_t r_plus_1 {r + 1U}; + + // (r+1)^2 may overflow uint128 when r is close to 2^64; in that case the + // invariant n < (r+1)^2 is trivially satisfied since n fits in 128 bits. + if (r_plus_1 != 0U && r_plus_1 <= ((std::numeric_limits::max)() / r_plus_1)) + { + BOOST_TEST(n < r_plus_1 * r_plus_1); + } +} + +} // namespace + +void test_uint128_isqrt_small() +{ + BOOST_TEST_EQ(isqrt(uint128_t{0}), uint128_t{0}); + BOOST_TEST_EQ(isqrt(uint128_t{1}), uint128_t{1}); + BOOST_TEST_EQ(isqrt(uint128_t{2}), uint128_t{1}); + BOOST_TEST_EQ(isqrt(uint128_t{3}), uint128_t{1}); + BOOST_TEST_EQ(isqrt(uint128_t{4}), uint128_t{2}); + BOOST_TEST_EQ(isqrt(uint128_t{5}), uint128_t{2}); + BOOST_TEST_EQ(isqrt(uint128_t{8}), uint128_t{2}); + BOOST_TEST_EQ(isqrt(uint128_t{9}), uint128_t{3}); + BOOST_TEST_EQ(isqrt(uint128_t{10}), uint128_t{3}); + BOOST_TEST_EQ(isqrt(uint128_t{15}), uint128_t{3}); + BOOST_TEST_EQ(isqrt(uint128_t{16}), uint128_t{4}); + BOOST_TEST_EQ(isqrt(uint128_t{99}), uint128_t{9}); + BOOST_TEST_EQ(isqrt(uint128_t{100}), uint128_t{10}); + BOOST_TEST_EQ(isqrt(uint128_t{101}), uint128_t{10}); + + // Exhaustive cross-check against the reference for every small n. + for (std::uint64_t i {0}; i < 200; ++i) + { + BOOST_TEST_EQ(isqrt(uint128_t{i}), isqrt_ref(uint128_t{i})); + } +} + +void test_uint128_isqrt_perfect_squares() +{ + // Squares spanning the full 64-bit range, including the largest k whose + // square still fits in 128 bits (k = 2^64 - 1). + for (std::uint64_t k {0}; k < 10000; ++k) + { + const uint128_t kk {k}; + BOOST_TEST_EQ(isqrt(kk * kk), kk); + + if (k > 0) + { + BOOST_TEST_EQ(isqrt(kk * kk - 1U), kk - 1U); + BOOST_TEST_EQ(isqrt(kk * kk + 2U * kk), kk); // (k+1)^2 - 1 + } + } + + // Powers of two as bases: k = 2^i for i in [0, 63] - largest exact square + // is (2^63)^2 = 2^126. + for (int i {0}; i < 64; ++i) + { + const uint128_t k {uint128_t{1} << i}; + BOOST_TEST_EQ(isqrt(k * k), k); + } + + // Largest representable perfect square: (2^64 - 1)^2 = 2^128 - 2^65 + 1. + const uint128_t k_max {(std::numeric_limits::max)()}; + BOOST_TEST_EQ(isqrt(k_max * k_max), k_max); +} + +void test_uint128_isqrt_bit_boundaries() +{ + // 2^(2k) has integer square root 2^k. + for (int k {0}; k < 64; ++k) + { + const uint128_t n {uint128_t{1} << (2 * k)}; + BOOST_TEST_EQ(isqrt(n), uint128_t{1} << k); + } + + // 2^(2k+1) has integer square root floor(2^(k+0.5)) = floor(sqrt(2) * 2^k). + // Check the invariant holds rather than hard-coding the value. + for (int k {0}; k < 63; ++k) + { + check_invariant(uint128_t{1} << (2 * k + 1)); + } + + // Just below and just above bit boundaries. + for (int k {2}; k < 128; ++k) + { + const uint128_t boundary {uint128_t{1} << k}; + check_invariant(boundary - 1U); + check_invariant(boundary); + check_invariant(boundary + 1U); + } +} + +void test_uint128_isqrt_extreme() +{ + // (uint128 max). (2^64)^2 = 2^128 wraps, so isqrt(2^128 - 1) = 2^64 - 1. + const uint128_t u128_max {(std::numeric_limits::max)()}; + const uint128_t u64_max {(std::numeric_limits::max)()}; + BOOST_TEST_EQ(isqrt(u128_max), u64_max); + + // 2^128 - 2^65 + 1 = (2^64 - 1)^2 - exact square at the very top. + BOOST_TEST_EQ(isqrt(u64_max * u64_max), u64_max); + + // One above the largest representable square: still has isqrt = 2^64 - 1. + BOOST_TEST_EQ(isqrt(u64_max * u64_max + 1U), u64_max); + + // A handful of large hand-picked values, cross-checked against the bit-by- + // bit reference. + const uint128_t big_a {UINT64_C(0x0123456789ABCDEF), UINT64_C(0xFEDCBA9876543210)}; + const uint128_t big_b {UINT64_C(0xDEADBEEFCAFEBABE), UINT64_C(0x0123456789ABCDEF)}; + const uint128_t big_c {UINT64_C(0x8000000000000000), 0U}; + + BOOST_TEST_EQ(isqrt(big_a), isqrt_ref(big_a)); + BOOST_TEST_EQ(isqrt(big_b), isqrt_ref(big_b)); + BOOST_TEST_EQ(isqrt(big_c), isqrt_ref(big_c)); + + check_invariant(big_a); + check_invariant(big_b); + check_invariant(big_c); + check_invariant(u128_max); +} + +void test_int128_isqrt() +{ + BOOST_TEST_EQ(isqrt(int128_t{0}), int128_t{0}); + BOOST_TEST_EQ(isqrt(int128_t{1}), int128_t{1}); + BOOST_TEST_EQ(isqrt(int128_t{2}), int128_t{1}); + BOOST_TEST_EQ(isqrt(int128_t{100}), int128_t{10}); + BOOST_TEST_EQ(isqrt(int128_t{144}), int128_t{12}); + BOOST_TEST_EQ(isqrt(int128_t{INT64_C(1000000000000000000)}), int128_t{INT64_C(1000000000)}); + + // int128 max = 2^127 - 1. floor(sqrt) = floor(2^63.5) = 6074001000.7e9 + // Use the unsigned implementation as the source of truth. + const int128_t i128_max {(std::numeric_limits::max)()}; + BOOST_TEST_EQ(isqrt(i128_max), static_cast(isqrt(static_cast(i128_max)))); + + // Negative inputs are documented to return 0. + BOOST_TEST_EQ(isqrt(int128_t{-1}), int128_t{0}); + BOOST_TEST_EQ(isqrt(int128_t{-100}), int128_t{0}); + BOOST_TEST_EQ(isqrt((std::numeric_limits::min)()), int128_t{0}); +} + +void test_isqrt_against_ipow() +{ + // isqrt(ipow(k, 2)) == k for any k whose square fits. + for (std::uint64_t k {0}; k < 1000; ++k) + { + BOOST_TEST_EQ(isqrt(ipow(uint128_t{k}, 2U)), uint128_t{k}); + } + + // isqrt is monotonically non-decreasing. + uint128_t prev {0}; + + for (std::uint64_t i {0}; i < 1000; ++i) + { + const uint128_t curr {isqrt(uint128_t{i})}; + BOOST_TEST(curr >= prev); + prev = curr; + } +} + +void test_constexpr_isqrt() +{ + constexpr uint128_t r1 {isqrt(uint128_t{0})}; + static_assert(r1 == uint128_t{0}, "isqrt(0) constexpr"); + + constexpr uint128_t r2 {isqrt(uint128_t{1})}; + static_assert(r2 == uint128_t{1}, "isqrt(1) constexpr"); + + constexpr uint128_t r3 {isqrt(uint128_t{100})}; + static_assert(r3 == uint128_t{10}, "isqrt(100) constexpr"); + + constexpr uint128_t r4 {isqrt(uint128_t{UINT64_C(1000000000000000000)})}; + static_assert(r4 == uint128_t{UINT64_C(1000000000)}, "isqrt(10^18) constexpr"); + + constexpr int128_t r5 {isqrt(int128_t{-5})}; + static_assert(r5 == int128_t{0}, "isqrt negative constexpr"); + + constexpr int128_t r6 {isqrt(int128_t{12321})}; + static_assert(r6 == int128_t{111}, "isqrt(12321) constexpr"); +} + +int main() +{ + test_uint128_isqrt_small(); + test_uint128_isqrt_perfect_squares(); + test_uint128_isqrt_bit_boundaries(); + test_uint128_isqrt_extreme(); + test_int128_isqrt(); + test_isqrt_against_ipow(); + test_constexpr_isqrt(); + + return boost::report_errors(); +} From 477622b5d5c1b340886da9afada47c9fd716b7c0 Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Fri, 15 May 2026 15:26:32 -0400 Subject: [PATCH 3/5] Exclude gcc-7 x32 --- test/test_isqrt.cpp | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/test/test_isqrt.cpp b/test/test_isqrt.cpp index 369f996c..2c009329 100644 --- a/test/test_isqrt.cpp +++ b/test/test_isqrt.cpp @@ -2,6 +2,14 @@ // Distributed under the Boost Software License, Version 1.0. // https://www.boost.org/LICENSE_1_0.txt +#if defined(__GNUC__) && __GNUC__ == 7 && defined(__i386__) + +// 32-bit GCC-7 fails with: "error: constexpr loop iteration count exceeds limit of 262144" + +int main() { return 0; } + +#else + #ifndef BOOST_INT128_BUILD_MODULE #include @@ -256,3 +264,5 @@ int main() return boost::report_errors(); } + +#endif From 6c85c620a4061f9e84790922a061b64a8ee53c9a Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Tue, 19 May 2026 13:59:19 -0400 Subject: [PATCH 4/5] Fix copy-paste error --- .drone.jsonnet | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/.drone.jsonnet b/.drone.jsonnet index 4aaae745..d7b0c56f 100644 --- a/.drone.jsonnet +++ b/.drone.jsonnet @@ -1,8 +1,9 @@ # Copyright 2022, 2023 Peter Dimov +# Copyright 2025 - 2026 Matt Borland # Distributed under the Boost Software License, Version 1.0. # https://www.boost.org/LICENSE_1_0.txt -local library = "decimal"; +local library = "int128"; local triggers = { From c8a4cd608a81a5af043f011f9a1fd6bf4cd429c4 Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Tue, 19 May 2026 14:01:51 -0400 Subject: [PATCH 5/5] Ignore MSVC warnings from 14.2 --- test/test_isqrt.cpp | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/test/test_isqrt.cpp b/test/test_isqrt.cpp index 2c009329..5d485e31 100644 --- a/test/test_isqrt.cpp +++ b/test/test_isqrt.cpp @@ -231,6 +231,12 @@ void test_isqrt_against_ipow() } } +#ifdef _MSC_VER +# pragma warning(push) +# pragma warning(disable : 4307) // integral constant overflow +# pragma warning(disable : 4308) // negative integral constant converted to unsigned type +#endif + void test_constexpr_isqrt() { constexpr uint128_t r1 {isqrt(uint128_t{0})}; @@ -252,6 +258,10 @@ void test_constexpr_isqrt() static_assert(r6 == int128_t{111}, "isqrt(12321) constexpr"); } +#ifdef _MSC_VER +# pragma warning(pop) +#endif + int main() { test_uint128_isqrt_small();