Examples

// Copyright 2025 Matt Borland
// Distributed under the Boost Software License, Version 1.0.
// https://www.boost.org/LICENSE_1_0.txt

#include <boost/int128/int128.hpp>
#include <boost/int128/literals.hpp>
#include <boost/int128/iostream.hpp>
#include <iostream>
#include <limits>
#include <sstream>

int main()
{
    using boost::int128::uint128;
    using boost::int128::int128;

    std::cout << "=== uint128 Construction ===" << std::endl;

    // 1) From a builtin integer type
    constexpr uint128 from_builtin {42U};
    std::cout << "From builtin (42U): " << from_builtin << std::endl;

    // 2) From high and low 64-bit values (high, low)
    constexpr uint128 from_parts {UINT64_C(0x1), UINT64_C(0x0)};  // 2^64
    std::cout << "From parts (1, 0) = 2^64: " << from_parts << std::endl;

    constexpr uint128 max_value {UINT64_C(0xFFFFFFFFFFFFFFFF), UINT64_C(0xFFFFFFFFFFFFFFFF)};
    std::cout << "From parts (max, max): " << max_value << std::endl;
    std::cout << "  Equals numeric_limits max? " << std::boolalpha
              << (max_value == std::numeric_limits<uint128>::max()) << std::endl;

    // 3) From user-defined literals.
    // The library provides only string-form UDLs
    // For small values like this a string is still parsed rather than direct construction
    // Using the constructors for values that fit in (unsigned) long long should be preferred for performance
    using namespace boost::int128::literals;
    const auto small_literal {12345_U128};
    std::cout << "From literal 12345_U128: " << small_literal << std::endl;

    // 4) From macro (like UINT64_C but for 128-bit), good for values that exceed unsigned long long
    const auto from_macro {BOOST_INT128_UINT128_C(340282366920938463463374607431768211455)};
    std::cout << "From BOOST_INT128_UINT128_C(max): " << from_macro << std::endl;

    // 5) From input stream
    std::stringstream ss;
    ss.str("12345678901234567890123456789");
    uint128 from_stream;
    ss >> from_stream;
    std::cout << "From stringstream: " << from_stream << std::endl;

    std::cout << "\n=== int128 Construction ===" << std::endl;

    // Signed from builtin
    constexpr int128 signed_builtin {-42};
    std::cout << "From builtin (-42): " << signed_builtin << std::endl;

    // Signed from parts. Both words are stored unsigned, but the constructor takes
    // the high word signed; read it back with signed_high().
    constexpr int128 min_value {INT64_MIN, 0};
    std::cout << "From parts (INT64_MIN, 0): " << min_value << std::endl;
    std::cout << "  Equals numeric_limits min? "
              << (min_value == std::numeric_limits<int128>::min()) << std::endl;

    // Signed literals. Values that fit in unsigned long long can be written
    // directly; the leading minus is parsed as a unary operator on the
    // literal result (lowercase and uppercase suffixes both work):
    const auto negative_literal {-12345_i128};
    std::cout << "From literal -12345_i128: " << negative_literal << std::endl;

    const auto positive_literal {12345_I128};
    std::cout << "From literal 12345_I128: " << positive_literal << std::endl;

    // For magnitudes beyond unsigned long long you can use the macro or a string literal
    const auto large_signed {BOOST_INT128_INT128_C(-99999999999999999999)};
    std::cout << "From BOOST_INT128_INT128_C(-99999999999999999999): " << large_signed << std::endl;

    const auto large_signed_string {"-99999999999999999999"_i128};
    std::cout << "From string literal: " << large_signed_string << std::endl;

    // Signed macro
    const auto from_signed_macro {BOOST_INT128_INT128_C(-170141183460469231731687303715884105728)};
    std::cout << "From BOOST_INT128_INT128_C(min): " << from_signed_macro << std::endl;

    std::cout << "\n=== Default and Copy Construction ===" << std::endl;

    // Default construction (zero-initialized)
    constexpr uint128 default_constructed {};
    std::cout << "Default constructed: " << default_constructed << std::endl;

    // Copy construction
    const uint128 copied {from_macro};
    std::cout << "Copy constructed: " << copied << std::endl;

    std::cout << "\n=== Floating-Point Construction ===" << std::endl;

    // Floating-point construction truncates toward zero, matching the behavior of
    // a static_cast from a floating-point type to a built-in integer.
    constexpr uint128 from_double {12345.9};
    std::cout << "uint128 from 12345.9 (truncated): " << from_double << std::endl;

    constexpr int128 from_negative_double {-12345.9};
    std::cout << "int128 from -12345.9 (truncated toward zero): " << from_negative_double << std::endl;

    // Values that exceed the 64-bit range are routed through the full 128-bit decomposition.
    const double two_to_the_100 {1.2676506002282294e30};  // 2^100
    const uint128 large_from_double {two_to_the_100};
    std::cout << "uint128 from 2^100: " << large_from_double << std::endl;

    std::cout << "\n=== Floating-Point Edge Cases ===" << std::endl;

    // NaN yields zero for both signed and unsigned (mirrors libgcc's __fix(uns)Xfti).
    const double nan_value {std::numeric_limits<double>::quiet_NaN()};
    const uint128 unsigned_from_nan {nan_value};
    const int128 signed_from_nan {nan_value};
    std::cout << "uint128 from NaN: " << unsigned_from_nan << std::endl;
    std::cout << "int128 from NaN: " << signed_from_nan << std::endl;

    // Negative values are clamped to zero when constructing uint128.
    const uint128 unsigned_from_negative {-1.0};
    std::cout << "uint128 from -1.0 (clamped to zero): " << unsigned_from_negative << std::endl;

    // Positive overflow saturates: anything >= 2^128 (including +infinity) becomes UINT128_MAX.
    const double infinity {std::numeric_limits<double>::infinity()};
    const uint128 saturated_unsigned {infinity};
    std::cout << "uint128 from +infinity (saturates to UINT128_MAX): " << saturated_unsigned << std::endl;

    // For int128, values >= 2^127 saturate to INT128_MAX and values <= -2^127 saturate to INT128_MIN.
    const double huge {1e40};  // Well beyond 2^127 (~ 1.7e38)
    const int128 saturated_positive {huge};
    const int128 saturated_negative {-huge};
    std::cout << "int128 from 1e40 (saturates to INT128_MAX): " << saturated_positive << std::endl;
    std::cout << "int128 from -1e40 (saturates to INT128_MIN): " << saturated_negative << std::endl;

    return 0;
}

=== uint128 Construction ===
From builtin (42U): 42
From parts (1, 0) = 2^64: 18446744073709551616
From parts (max, max): 340282366920938463463374607431768211455
  Equals numeric_limits max? true
From literal 12345_U128: 12345
From BOOST_INT128_UINT128_C(max): 340282366920938463463374607431768211455
From stringstream: 12345678901234567890123456789

=== int128 Construction ===
From builtin (-42): -42
From parts (INT64_MIN, 0): -170141183460469231731687303715884105728
  Equals numeric_limits min? true
From literal -12345_i128: -12345
From literal 12345_I128: 12345
From BOOST_INT128_INT128_C(-99999999999999999999): -99999999999999999999
From string literal: -99999999999999999999
From BOOST_INT128_INT128_C(min): -170141183460469231731687303715884105728

=== Default and Copy Construction ===
Default constructed: 0
Copy constructed: 340282366920938463463374607431768211455

=== Floating-Point Construction ===
uint128 from 12345.9 (truncated): 12345
int128 from -12345.9 (truncated toward zero): -12345
uint128 from 2^100: 1267650600228229401496703205376

=== Floating-Point Edge Cases ===
uint128 from NaN: 0
int128 from NaN: 0
uint128 from -1.0 (clamped to zero): 0
uint128 from +infinity (saturates to UINT128_MAX): 340282366920938463463374607431768211455
int128 from 1e40 (saturates to INT128_MAX): 170141183460469231731687303715884105727
int128 from -1e40 (saturates to INT128_MIN): -170141183460469231731687303715884105728

// Copyright 2025 Matt Borland
// Distributed under the Boost Software License, Version 1.0.
// https://www.boost.org/LICENSE_1_0.txt

#include <boost/int128/int128.hpp>
#include <boost/int128/iostream.hpp>
#include <iostream>

int main()
{
    // The types of this library support all arithmetic operations one would expect.
    // They can be between values of the same type, or the same signedness by default.
    // See `mixed_type_arithmetic.cpp` for operations with different signedness.

    using boost::int128::int128;
    using boost::int128::uint128;

    // Basic arithmetic with signed 128-bit integers
    int128 x {1000000000000LL};  // 1 trillion
    int128 y {999999999999LL};   // Just under 1 trillion

    std::cout << "=== Signed 128-bit Arithmetic ===" << std::endl;
    std::cout << "x = " << x << std::endl;
    std::cout << "y = " << y << std::endl;

    // Addition and subtraction
    std::cout << "\nAddition and Subtraction:" << std::endl;
    std::cout << "x + y = " << (x + y) << std::endl;
    std::cout << "x - y = " << (x - y) << std::endl;

    // Multiplication - results that exceed 64-bit range
    std::cout << "\nMultiplication (exceeds 64-bit range):" << std::endl;
    std::cout << "x * y = " << (x * y) << std::endl;

    // Division and modulo
    std::cout << "\nDivision and Modulo:" << std::endl;
    std::cout << "x / 7 = " << (x / 7) << std::endl;
    std::cout << "x % 7 = " << (x % 7) << std::endl;

    // Comparisons
    std::cout << "\nComparisons:" << std::endl;
    std::cout << "x > y:  " << std::boolalpha << (x > y) << std::endl;
    std::cout << "x == y: " << (x == y) << std::endl;
    std::cout << "x != y: " << (x != y) << std::endl;

    // Negative values and absolute value
    int128 negative {-42};
    std::cout << "\nNegative values:" << std::endl;
    std::cout << "negative = " << negative << std::endl;
    std::cout << "abs(negative) = " << boost::int128::abs(negative) << std::endl;

    // Compound assignment operators
    std::cout << "\nCompound assignment operators:" << std::endl;
    int128 z {100};
    std::cout << "z = " << z << std::endl;

    z += 50;
    std::cout << "z += 50: " << z << std::endl;

    z -= 25;
    std::cout << "z -= 25: " << z << std::endl;

    z *= 2;
    std::cout << "z *= 2:  " << z << std::endl;

    z /= 5;
    std::cout << "z /= 5:  " << z << std::endl;

    z %= 7;
    std::cout << "z %= 7:  " << z << std::endl;

    // Unsigned 128-bit arithmetic - useful for very large positive values
    std::cout << "\n=== Unsigned 128-bit Arithmetic ===" << std::endl;
    uint128 large {UINT64_C(0x1), UINT64_C(0x0)};  // 2^64
    std::cout << "large (2^64) = " << large << std::endl;
    std::cout << "large * 2 = " << (large * 2U) << std::endl;
    std::cout << "large + large = " << (large + large) << std::endl;

    // Increment and decrement
    std::cout << "\nIncrement and Decrement:" << std::endl;
    int128 counter {10};
    std::cout << "counter = " << counter << std::endl;
    std::cout << "++counter = " << ++counter << std::endl;
    std::cout << "counter++ = " << counter++ << std::endl;
    std::cout << "counter = " << counter << std::endl;
    std::cout << "--counter = " << --counter << std::endl;

    return 0;
}

=== Signed 128-bit Arithmetic ===
x = 1000000000000
y = 999999999999

Addition and Subtraction:
x + y = 1999999999999
x - y = 1

Multiplication (exceeds 64-bit range):
x * y = 999999999999000000000000

Division and Modulo:
x / 7 = 142857142857
x % 7 = 1

Comparisons:
x > y:  true
x == y: false
x != y: true

Negative values:
negative = -42
abs(negative) = 42

Compound assignment operators:
z = 100
z += 50: 150
z -= 25: 125
z *= 2:  250
z /= 5:  50
z %= 7:  1

=== Unsigned 128-bit Arithmetic ===
large (2^64) = 18446744073709551616
large * 2 = 36893488147419103232
large + large = 36893488147419103232

Increment and Decrement:
counter = 10
++counter = 11
counter++ = 11
counter = 12
--counter = 11

// Copyright 2025 Matt Borland
// Distributed under the Boost Software License, Version 1.0.
// https://www.boost.org/LICENSE_1_0.txt

#include <boost/int128/int128.hpp>
#include <boost/int128/iostream.hpp>
#include <iostream>
#include <sstream>

int main()
{
    using boost::int128::int128;
    using boost::int128::uint128;

    std::cout << "=== Basic Streaming ===" << std::endl;

    // Both types allow streaming as one would expect from a regular builtin-type
    constexpr int128 signed_value {-42};
    std::cout << "Signed value: " << signed_value << std::endl;

    // We can also use <iomanip> to change the output format
    constexpr uint128 unsigned_value {0x1, UINT64_MAX};
    std::cout << "Unsigned value (dec): " << unsigned_value << '\n'
              << "Unsigned value (hex): " << std::hex << unsigned_value << '\n'
              << "Unsigned value (oct): " << std::oct << unsigned_value << std::endl;

    // Hex also can be manipulated to be uppercase
    std::cout << "Upper unsigned value: " << std::hex << std::uppercase << unsigned_value << std::endl;

    // And returned to default formatting
    std::cout << "Lower unsigned value: " << std::dec << std::nouppercase << unsigned_value << std::endl;

    // Large values that exceed 64-bit range
    std::cout << "\n=== Large Values (Beyond 64-bit) ===" << std::endl;

    // 2^64 = 18446744073709551616 (first value that doesn't fit in uint64_t)
    constexpr uint128 two_to_64 {1, 0};
    std::cout << "2^64 = " << two_to_64 << std::endl;

    // 2^100 = a very large number
    constexpr uint128 two_to_100 {uint128{1} << 100};
    std::cout << "2^100 = " << two_to_100 << std::endl;

    // Maximum uint128 value
    constexpr auto uint_max {std::numeric_limits<uint128>::max()};
    std::cout << "uint128 max = " << uint_max << std::endl;

    // Minimum and maximum int128 values
    constexpr auto int_min {std::numeric_limits<int128>::min()};
    constexpr auto int_max {std::numeric_limits<int128>::max()};
    std::cout << "int128 min = " << int_min << std::endl;
    std::cout << "int128 max = " << int_max << std::endl;

    // String conversion using stringstream
    std::cout << "\n=== String Conversion with std::stringstream ===" << std::endl;

    // Convert uint128 to string
    std::ostringstream oss;
    oss << two_to_100;
    auto str {oss.str()};
    std::cout << "uint128 to string: \"" << str << "\"" << std::endl;

    // Convert string to uint128
    std::istringstream iss {"123456789012345678901234567890"};
    uint128 parsed_value {};
    iss >> parsed_value;
    std::cout << "String to uint128: " << parsed_value << std::endl;

    // Round-trip: value -> string -> value
    std::cout << "\n=== Round-trip Conversion ===" << std::endl;
    constexpr uint128 original {0xDEADBEEF, 0xCAFEBABE12345678};
    std::ostringstream oss2;
    oss2 << original;
    auto original_str {oss2.str()};

    std::istringstream iss2 {original_str};
    uint128 round_tripped {};
    iss2 >> round_tripped;

    std::cout << "Original:     " << original << std::endl;
    std::cout << "As string:    \"" << original_str << "\"" << std::endl;
    std::cout << "Round-tripped: " << round_tripped << std::endl;
    std::cout << "Match: " << std::boolalpha << (original == round_tripped) << std::endl;

    return 0;
}

=== Basic Streaming ===
Signed value: -42
Unsigned value (dec): 36893488147419103231
Unsigned value (hex): 1ffffffffffffffff
Unsigned value (oct): 3777777777777777777777
Upper unsigned value: 1FFFFFFFFFFFFFFFF
Lower unsigned value: 36893488147419103231

=== Large Values (Beyond 64-bit) ===
2^64 = 18446744073709551616
2^100 = 1267650600228229401496703205376
uint128 max = 340282366920938463463374607431768211455
int128 min = -170141183460469231731687303715884105728
int128 max = 170141183460469231731687303715884105727

=== String Conversion with std::stringstream ===
uint128 to string: "1267650600228229401496703205376"
String to uint128: 123456789012345678901234567890

=== Round-trip Conversion ===
Original:     68915718020162848918556923512
As string:    "68915718020162848918556923512"
Round-tripped: 68915718020162848918556923512
Match: true

// Copyright 2025 Matt Borland
// Distributed under the Boost Software License, Version 1.0.
// https://www.boost.org/LICENSE_1_0.txt
//
// This example demonstrates the rollover behavior for both signed and unsigned int128

#include <boost/int128.hpp>
#include <iostream>
#include <iomanip>
#include <limits>

int main()
{
    using boost::int128::uint128;
    using boost::int128::int128;

    constexpr uint128 max_unsigned_value {std::numeric_limits<uint128>::max()};
    constexpr uint128 min_unsigned_value {std::numeric_limits<uint128>::min()};

    std::cout << "=== uint128 behavior ===" << std::endl;

    std::cout << "Max of uint128: " << max_unsigned_value << '\n'
              << "Max + 1U: " << max_unsigned_value + 1U << "\n\n";

    std::cout << "Min of uint128: " << min_unsigned_value << '\n'
              << "Min - 1U: " << min_unsigned_value - 1U << "\n\n";

    constexpr int128 max_signed_value {std::numeric_limits<int128>::max()};
    constexpr int128 min_signed_value {std::numeric_limits<int128>::min()};

    std::cout << "=== int128 behavior ===" << std::endl;

    std::cout << "Max of int128: " << max_signed_value << '\n'
              << "Max + 1: " << max_signed_value + 1 << "\n";

    std::cout << "\nMin of int128: " << min_signed_value << '\n'
              << "Min - 1: " << min_signed_value - 1 << '\n' << std::endl;
}

=== uint128 behavior ===
Max of uint128: 340282366920938463463374607431768211455
Max + 1U: 0

Min of uint128: 0
Min - 1U: 340282366920938463463374607431768211455

=== int128 behavior ===
Max of int128: 170141183460469231731687303715884105727
Max + 1: -170141183460469231731687303715884105728

Min of int128: -170141183460469231731687303715884105728
Min - 1: 170141183460469231731687303715884105727

// Copyright 2025 Matt Borland
// Distributed under the Boost Software License, Version 1.0.
// https://www.boost.org/LICENSE_1_0.txt

// Individual headers

#include <boost/int128/int128.hpp>
#include <boost/int128/bit.hpp>

// Or you can do a single header

// #include <boost/int128.hpp>

int main()
{
    // The functions from bit are only available for uint128

    constexpr boost::int128::uint128 x {1U};

    // All the functions are constexpr

    // Does the value have only a single bit set?
    static_assert(boost::int128::has_single_bit(x), "Should have one bit");

    // How many zeros from the left
    static_assert(boost::int128::countl_zero(x) == 127U, "Should be 127");

    // The bit width of the value
    // 1 + 1 is 10 in binary which is 2 bits wide
    static_assert(boost::int128::bit_width(x + x) == 2U, "2 bits wide");

    // The smallest power of two not greater than the input value
    static_assert(boost::int128::bit_floor(3U * x) == 2U, "2 < 3");

    // The smallest power of two not Smaller than the input value
    static_assert(boost::int128::bit_ceil(5U * x) == 8U, "8 > 5");

    // How many zeros from the right?
    static_assert(boost::int128::countr_zero(2U * x) == 1, "1 zero to the right of 10");

    // How many 1-bits in the value
    static_assert(boost::int128::popcount(7U * x) == 3, "111");

    // Swap the bytes
    // Create a value with distinct byte pattern
    constexpr boost::int128::uint128 original{
        0x0123456789ABCDEFULL,
        0xFEDCBA9876543210ULL
    };

    // Expected result after byteswap
    constexpr boost::int128::uint128 expected{
        0x1032547698BADCFEULL,
        0xEFCDAB8967452301ULL
    };

    static_assert(boost::int128::byteswap(original) == expected, "Mismatched byteswap");
    static_assert(boost::int128::byteswap(expected) == original, "Mismatched byteswap");

    return 0;
}

// Copyright 2025 Matt Borland
// Distributed under the Boost Software License, Version 1.0.
// https://www.boost.org/LICENSE_1_0.txt

// Individual headers

#include <boost/int128/int128.hpp>
#include <boost/int128/numeric.hpp>
#include <boost/int128/iostream.hpp>

// Or you can do a single header

// #include <boost/int128.hpp>

#include <limits>
#include <type_traits>
#include <iostream>

int main()
{
    using boost::int128::uint128;
    using boost::int128::int128;

    // std::numeric_limits is overloaded for both types
    constexpr auto uint_max {std::numeric_limits<uint128>::max()};
    static_assert(std::is_same<decltype(uint_max), const uint128>::value, "Types should match");

    constexpr auto int_max {std::numeric_limits<int128>::max()};
    constexpr auto int_min {std::numeric_limits<int128>::min()};

    std::cout << "=== Saturating Arithmetic ===" << std::endl;
    std::cout << "uint128 max = " << uint_max << std::endl;
    std::cout << "int128 max  = " << int_max << std::endl;
    std::cout << "int128 min  = " << int_min << std::endl;

    // Saturating arithmetic returns max on overflow, or min on underflow rather than rolling over
    std::cout << "\n=== Saturating Addition and Subtraction ===" << std::endl;
    std::cout << "saturating_add(uint_max, uint_max) = " << boost::int128::saturating_add(uint_max, uint_max)
              << " (saturates to uint_max)" << std::endl;
    std::cout << "saturating_sub(0, uint_max) = " << boost::int128::saturating_sub(uint128{0}, uint_max)
              << " (saturates to 0, not underflow)" << std::endl;

    // This is especially useful for signed types since rollover is undefined
    std::cout << "\n=== Saturating Multiplication ===" << std::endl;
    std::cout << "saturating_mul(int_max, 2) = " << boost::int128::saturating_mul(int_max, 2)
              << " (saturates to int_max)" << std::endl;
    std::cout << "saturating_mul(-(int_max - 2), 5) = " << boost::int128::saturating_mul(-(int_max - 2), 5)
              << " (saturates to int_min)" << std::endl;

    // The only case in the library where saturating_div overflows is x = int_min and y = -1
    std::cout << "\n=== Saturating Division ===" << std::endl;
    std::cout << "saturating_div(int_min, -1) = " << boost::int128::saturating_div(int_min, -1)
              << " (saturates to int_max; normally this overflows)" << std::endl;

    // saturating_cast allows types to be safely converted without rollover behavior
    std::cout << "\n=== Saturating Casts ===" << std::endl;
    std::cout << "saturating_cast<int128>(uint_max) = " << boost::int128::saturating_cast<int128>(uint_max)
              << " (saturates to int_max)" << std::endl;

    // You can also cast to builtin types
    std::cout << "saturating_cast<int64_t>(int_max) = " << boost::int128::saturating_cast<std::int64_t>(int_max)
              << " (saturates to INT64_MAX)" << std::endl;

    // Even of different signedness as this is treated like a static cast
    std::cout << "saturating_cast<int32_t>(uint_max) = " << boost::int128::saturating_cast<std::int32_t>(uint_max)
              << " (saturates to INT32_MAX)" << std::endl;

    return 0;
}

=== Saturating Arithmetic ===
uint128 max = 340282366920938463463374607431768211455
int128 max  = 170141183460469231731687303715884105727
int128 min  = -170141183460469231731687303715884105728

=== Saturating Addition and Subtraction ===
saturating_add(uint_max, uint_max) = 340282366920938463463374607431768211455 (saturates to uint_max)
saturating_sub(0, uint_max) = 0 (saturates to 0, not underflow)

=== Saturating Multiplication ===
saturating_mul(int_max, 2) = 170141183460469231731687303715884105727 (saturates to int_max)
saturating_mul(-(int_max - 2), 5) = -170141183460469231731687303715884105728 (saturates to int_min)

=== Saturating Division ===
saturating_div(int_min, -1) = 170141183460469231731687303715884105727 (saturates to int_max; normally this overflows)

=== Saturating Casts ===
saturating_cast<int128>(uint_max) = 170141183460469231731687303715884105727 (saturates to int_max)
saturating_cast<int64_t>(int_max) = 9223372036854775807 (saturates to INT64_MAX)
saturating_cast<int32_t>(uint_max) = 2147483647 (saturates to INT32_MAX)

// Copyright 2025 Matt Borland
// Distributed under the Boost Software License, Version 1.0.
// https://www.boost.org/LICENSE_1_0.txt

#include <boost/int128/int128.hpp>
#include <boost/int128/numeric.hpp>
#include <boost/int128/iostream.hpp>
#include <iostream>

int main()
{
    using boost::int128::uint128;
    using boost::int128::int128;

    std::cout << "=== Greatest Common Divisor (gcd) ===" << std::endl;

    // Basic gcd
    constexpr uint128 a {48};
    constexpr uint128 b {18};
    std::cout << "gcd(" << a << ", " << b << ") = " << boost::int128::gcd(a, b) << std::endl;

    // gcd with larger values
    constexpr uint128 large_a {123456789012345678ULL};
    constexpr uint128 large_b {987654321098765432ULL};
    std::cout << "gcd(" << large_a << ", " << large_b << ") = "
              << boost::int128::gcd(large_a, large_b) << std::endl;

    // gcd with 128-bit values
    constexpr uint128 huge_a {uint128{1} << 100};
    constexpr uint128 huge_b {uint128{1} << 80};
    std::cout << "gcd(2^100, 2^80) = " << boost::int128::gcd(huge_a, huge_b) << " (= 2^80)" << std::endl;

    // Signed gcd (always returns positive)
    constexpr int128 neg_a {-48};
    constexpr int128 neg_b {18};
    std::cout << "gcd(" << neg_a << ", " << neg_b << ") = "
              << boost::int128::gcd(neg_a, neg_b) << " (always positive)" << std::endl;

    std::cout << "\n=== Least Common Multiple (lcm) ===" << std::endl;

    // Basic lcm
    constexpr uint128 x {12};
    constexpr uint128 y {18};
    std::cout << "lcm(" << x << ", " << y << ") = " << boost::int128::lcm(x, y) << std::endl;

    // lcm with coprime numbers
    constexpr uint128 p {7};
    constexpr uint128 q {11};
    std::cout << "lcm(" << p << ", " << q << ") = " << boost::int128::lcm(p, q)
              << " (coprime: lcm = p * q)" << std::endl;

    // Relationship: gcd(a,b) * lcm(a,b) = a * b
    std::cout << "\nVerifying gcd * lcm = a * b:" << std::endl;
    auto g {boost::int128::gcd(x, y)};
    auto l {boost::int128::lcm(x, y)};
    std::cout << "gcd(" << x << ", " << y << ") * lcm(" << x << ", " << y << ") = "
              << (g * l) << std::endl;
    std::cout << x << " * " << y << " = " << (x * y) << std::endl;

    std::cout << "\n=== Midpoint ===" << std::endl;

    // Unsigned midpoint
    constexpr uint128 low {10};
    constexpr uint128 high {20};
    std::cout << "midpoint(" << low << ", " << high << ") = "
              << boost::int128::midpoint(low, high) << std::endl;

    // Midpoint with odd sum (rounds toward first argument)
    constexpr uint128 odd_low {10};
    constexpr uint128 odd_high {21};
    std::cout << "midpoint(" << odd_low << ", " << odd_high << ") = "
              << boost::int128::midpoint(odd_low, odd_high) << " (rounds toward first arg)" << std::endl;
    std::cout << "midpoint(" << odd_high << ", " << odd_low << ") = "
              << boost::int128::midpoint(odd_high, odd_low) << " (rounds toward first arg)" << std::endl;

    // Midpoint avoids overflow (unlike (a+b)/2)
    std::cout << "\n--- Overflow-safe midpoint ---" << std::endl;
    constexpr auto uint_max {std::numeric_limits<uint128>::max()};
    constexpr auto uint_max_minus_10 {uint_max - 10U};
    std::cout << "midpoint(uint128_max, uint128_max - 10) = "
              << boost::int128::midpoint(uint_max, uint_max_minus_10) << std::endl;
    std::cout << "(This would overflow if computed as (a + b) / 2)" << std::endl;

    // Signed midpoint
    std::cout << "\n--- Signed midpoint ---" << std::endl;
    constexpr int128 neg {-100};
    constexpr int128 pos {100};
    std::cout << "midpoint(" << neg << ", " << pos << ") = "
              << boost::int128::midpoint(neg, pos) << std::endl;

    constexpr int128 neg2 {-100};
    constexpr int128 neg3 {-50};
    std::cout << "midpoint(" << neg2 << ", " << neg3 << ") = "
              << boost::int128::midpoint(neg2, neg3) << std::endl;

    return 0;
}

=== Greatest Common Divisor (gcd) ===
gcd(48, 18) = 6
gcd(123456789012345678, 987654321098765432) = 2
gcd(2^100, 2^80) = 1208925819614629174706176 (= 2^80)
gcd(-48, 18) = 6 (always positive)

=== Least Common Multiple (lcm) ===
lcm(12, 18) = 36
lcm(7, 11) = 77 (coprime: lcm = p * q)

Verifying gcd * lcm = a * b:
gcd(12, 18) * lcm(12, 18) = 216
12 * 18 = 216

=== Midpoint ===
midpoint(10, 20) = 15
midpoint(10, 21) = 15 (rounds toward first arg)
midpoint(21, 10) = 16 (rounds toward first arg)

--- Overflow-safe midpoint ---
midpoint(uint128_max, uint128_max - 10) = 340282366920938463463374607431768211450
(This would overflow if computed as (a + b) / 2)

--- Signed midpoint ---
midpoint(-100, 100) = 0
midpoint(-100, -50) = -75

// Copyright 2026 Matt Borland
// Distributed under the Boost Software License, Version 1.0.
// https://www.boost.org/LICENSE_1_0.txt

#include <boost/int128/int128.hpp>
#include <boost/int128/numeric.hpp>
#include <boost/int128/iostream.hpp>
#include <iostream>
#include <limits>

// tag::exclude[]
// Only present with MSVC 14.1
#ifdef _MSC_VER
#  pragma warning(push)
#  pragma warning(disable : 4307) // integral constant overflow
#  pragma warning(disable : 4146) // unary minus operator applied to unsigned type, result still unsigned
#endif
// end::exclude[]

int main()
{
    using boost::int128::uint128;
    using boost::int128::int128;

    std::cout << "=== Every rounding mode on -12 / 5 ===" << std::endl;

    // The exact quotient is -2.4, so the nearest integer is -2 and no mode ties
    constexpr int128 x {-12};
    constexpr int128 y {5};

    std::cout << "div_to_zero        = " << boost::int128::div_to_zero(x, y) << std::endl;
    std::cout << "div_away_zero      = " << boost::int128::div_away_zero(x, y) << std::endl;
    std::cout << "div_to_pos_inf     = " << boost::int128::div_to_pos_inf(x, y) << std::endl;
    std::cout << "div_to_neg_inf     = " << boost::int128::div_to_neg_inf(x, y) << std::endl;
    std::cout << "div_euclid         = " << boost::int128::div_euclid(x, y) << std::endl;
    std::cout << "div_ties_to_zero   = " << boost::int128::div_ties_to_zero(x, y) << std::endl;

    std::cout << "\n=== Tie breaking on -7 / 2 ===" << std::endl;

    // The exact quotient is -3.5, so every tie-breaking rule picks a different side
    constexpr int128 tie_x {-7};
    constexpr int128 tie_y {2};

    std::cout << "div_ties_to_zero     = " << boost::int128::div_ties_to_zero(tie_x, tie_y) << std::endl;
    std::cout << "div_ties_away_zero   = " << boost::int128::div_ties_away_zero(tie_x, tie_y) << std::endl;
    std::cout << "div_ties_to_pos_inf  = " << boost::int128::div_ties_to_pos_inf(tie_x, tie_y) << std::endl;
    std::cout << "div_ties_to_neg_inf  = " << boost::int128::div_ties_to_neg_inf(tie_x, tie_y) << std::endl;
    std::cout << "div_ties_to_odd      = " << boost::int128::div_ties_to_odd(tie_x, tie_y) << std::endl;
    std::cout << "div_ties_to_even     = " << boost::int128::div_ties_to_even(tie_x, tie_y) << std::endl;

    std::cout << "\n=== Quotient and remainder from one division ===" << std::endl;

    // Each div_rem_ function performs a single division and returns both halves
    constexpr auto floored {boost::int128::div_rem_to_neg_inf(x, y)};
    std::cout << "div_rem_to_neg_inf(-12, 5): quotient = " << floored.quotient
              << ", remainder = " << floored.remainder << std::endl;

    constexpr auto truncated {boost::int128::div_rem_to_zero(x, y)};
    std::cout << "div_rem_to_zero(-12, 5):    quotient = " << truncated.quotient
              << ", remainder = " << truncated.remainder << std::endl;

    // The remainder always satisfies x == quotient * y + remainder
    std::cout << "quotient * y + remainder = " << (floored.quotient * y + floored.remainder) << std::endl;

    std::cout << "\n=== Euclidean remainder is never negative ===" << std::endl;

    // operator% takes its sign from the dividend, which makes it a poor fit for
    // wrapping an offset into a range. rem_euclid always lands in [0, abs(y)).
    constexpr int128 modulus {7};
    for (int128 offset {-9}; offset <= -6; ++offset)
    {
        std::cout << offset << " % 7 = " << (offset % modulus)
                  << ", rem_euclid(" << offset << ", 7) = " << boost::int128::rem_euclid(offset, modulus) << std::endl;
    }

    std::cout << "\n=== Ceiling division without overflow ===" << std::endl;

    // Counting fixed size blocks needed to cover a length is the classic use for
    // rounding towards positive infinity. The usual (length + block - 1) / block
    // overflows here, while div_to_pos_inf does not.
    constexpr auto length {(std::numeric_limits<uint128>::max)()};
    constexpr uint128 block {1000};

    std::cout << "length = " << length << std::endl;
    std::cout << "div_to_pos_inf(length, 1000) = " << boost::int128::div_to_pos_inf(length, block) << std::endl;
    std::cout << "div_to_zero(length, 1000)    = " << boost::int128::div_to_zero(length, block) << std::endl;
    std::cout << "(length + 999) / 1000        = " << ((length + 999U) / block) << " (wrapped)" << std::endl;

    std::cout << "\n=== Unbiased rounding of a scaled value ===" << std::endl;

    // Rounding half to even keeps a long running sum from drifting upwards, which is
    // what operator/ plus a manual half-adjustment would do
    constexpr int128 scale {1000};
    const int128 samples[] {int128{1500}, int128{2500}, int128{3500}, int128{-1500}};

    for (const auto sample : samples)
    {
        std::cout << sample << " / 1000: ties_to_even = " << boost::int128::div_ties_to_even(sample, scale)
                  << ", ties_away_zero = " << boost::int128::div_ties_away_zero(sample, scale) << std::endl;
    }

    return 0;
}

=== Every rounding mode on -12 / 5 ===
div_to_zero        = -2
div_away_zero      = -3
div_to_pos_inf     = -2
div_to_neg_inf     = -3
div_euclid         = -3
div_ties_to_zero   = -2

=== Tie breaking on -7 / 2 ===
div_ties_to_zero     = -3
div_ties_away_zero   = -4
div_ties_to_pos_inf  = -3
div_ties_to_neg_inf  = -4
div_ties_to_odd      = -3
div_ties_to_even     = -4

=== Quotient and remainder from one division ===
div_rem_to_neg_inf(-12, 5): quotient = -3, remainder = 3
div_rem_to_zero(-12, 5):    quotient = -2, remainder = -2
quotient * y + remainder = -12

=== Euclidean remainder is never negative ===
-9 % 7 = -2, rem_euclid(-9, 7) = 5
-8 % 7 = -1, rem_euclid(-8, 7) = 6
-7 % 7 = 0, rem_euclid(-7, 7) = 0
-6 % 7 = -6, rem_euclid(-6, 7) = 1

=== Ceiling division without overflow ===
length = 340282366920938463463374607431768211455
div_to_pos_inf(length, 1000) = 340282366920938463463374607431768212
div_to_zero(length, 1000)    = 340282366920938463463374607431768211
(length + 999) / 1000        = 0 (wrapped)

=== Unbiased rounding of a scaled value ===
1500 / 1000: ties_to_even = 2, ties_away_zero = 2
2500 / 1000: ties_to_even = 2, ties_away_zero = 3
3500 / 1000: ties_to_even = 4, ties_away_zero = 4
-1500 / 1000: ties_to_even = -2, ties_away_zero = -2

// Copyright 2026 Matt Borland
// Distributed under the Boost Software License, Version 1.0.
// https://www.boost.org/LICENSE_1_0.txt

// Individual headers

#include <boost/int128/utilities.hpp>
#include <boost/int128/iostream.hpp>

// Or you can do a single header

// #include <boost/int128.hpp>

#include <cstdint>
#include <limits>
#include <iostream>

int main()
{
    using boost::int128::uint128;
    using boost::int128::int128;
    using boost::int128::ckd_add;
    using boost::int128::ckd_sub;
    using boost::int128::ckd_mul;

    std::cout << std::boolalpha;

    // ckd_add, ckd_sub, and ckd_mul implement the C23 stdckdint.h contract: the
    // operation is evaluated as if both operands had infinite range, the result
    // is written to *result wrapped to that type's width, and the function
    // returns true when the exact result did not fit.
    constexpr auto u_max {std::numeric_limits<uint128>::max()};
    constexpr auto i_max {std::numeric_limits<int128>::max()};
    constexpr auto i_min {std::numeric_limits<int128>::min()};

    // A result that fits returns false and holds the exact value.
    std::cout << "=== Results That Fit ===" << std::endl;
    int128 r {};
    bool overflow {ckd_add(&r, int128{20}, int128{22})};
    std::cout << "ckd_add(20, 22): overflow=" << overflow << ", result=" << r << std::endl;

    // Addition that exceeds the type wraps modulo 2^128 and reports overflow.
    std::cout << "\n=== Addition Overflow ===" << std::endl;
    uint128 u {};
    overflow = ckd_add(&u, u_max, uint128{1});
    std::cout << "ckd_add(UINT128_MAX, 1): overflow=" << overflow << ", wrapped=" << u << std::endl;

    // Subtracting below zero in an unsigned type wraps to the top of the range.
    std::cout << "\n=== Subtraction Underflow ===" << std::endl;
    overflow = ckd_sub(&u, uint128{0}, uint128{1});
    std::cout << "ckd_sub(0, 1): overflow=" << overflow << ", wrapped=" << u << std::endl;

    // Multiplication detects overflow that operator* would silently roll over,
    // including INT128_MIN * -1, whose true result is not representable.
    std::cout << "\n=== Multiplication Overflow ===" << std::endl;
    overflow = ckd_mul(&r, i_max, int128{2});
    std::cout << "ckd_mul(INT128_MAX, 2): overflow=" << overflow << ", wrapped=" << r << std::endl;
    overflow = ckd_mul(&r, i_min, int128{-1});
    std::cout << "ckd_mul(INT128_MIN, -1): overflow=" << overflow << ", wrapped=" << r << std::endl;

    // The result type and the two operand types are independent: they may differ
    // in width and signedness, and the exact mathematical value is always used.
    std::cout << "\n=== Mixed Types ===" << std::endl;
    std::int64_t narrow {};
    overflow = ckd_add(&narrow, uint128{5}, int128{-3});
    std::cout << "ckd_add<int64_t>(uint128{5}, int128{-3}): overflow=" << overflow
              << ", result=" << narrow << std::endl;

    // Narrow targets make the wrap-around easy to see (400 modulo 256 is 144).
    std::uint8_t byte {};
    overflow = ckd_mul(&byte, std::uint8_t{20}, std::uint8_t{20});
    std::cout << "ckd_mul<uint8_t>(20, 20): overflow=" << overflow
              << ", wrapped=" << static_cast<int>(byte) << std::endl;

    return 0;
}

=== Results That Fit ===
ckd_add(20, 22): overflow=false, result=42

=== Addition Overflow ===
ckd_add(UINT128_MAX, 1): overflow=true, wrapped=0

=== Subtraction Underflow ===
ckd_sub(0, 1): overflow=true, wrapped=340282366920938463463374607431768211455

=== Multiplication Overflow ===
ckd_mul(INT128_MAX, 2): overflow=true, wrapped=-2
ckd_mul(INT128_MIN, -1): overflow=true, wrapped=-170141183460469231731687303715884105728

=== Mixed Types ===
ckd_add<int64_t>(uint128{5}, int128{-3}): overflow=false, result=2
ckd_mul<uint8_t>(20, 20): overflow=true, wrapped=144

// Copyright 2025 Matt Borland
// Distributed under the Boost Software License, Version 1.0.
// https://www.boost.org/LICENSE_1_0.txt

#include <boost/int128.hpp>
#include <iostream>

int main()
{
    // Mixed-sign comparisons and arithmetic between int128, uint128, and
    // built-in integer types of opposite signedness follow the C++ usual
    // arithmetic conversions, identical to the built-in __int128 /
    // unsigned __int128 types.

    std::cout << "=== Mixed Type Arithmetic with uint128 ===" << std::endl;

    constexpr boost::int128::uint128 unsigned_value {3};
    std::cout << "unsigned_value = " << unsigned_value << std::endl;

    constexpr auto greater_unsigned_value {unsigned_value + 5};

    std::cout << "unsigned_value + 1 = " << (unsigned_value + 1) << std::endl;
    std::cout << "unsigned_value - 1 = " << (unsigned_value - 1) << std::endl;
    std::cout << "unsigned_value * 2 = " << (unsigned_value * 2) << std::endl;
    std::cout << "unsigned_value / 3 = " << (unsigned_value / 3) << std::endl;
    std::cout << "unsigned_value % 3 = " << (unsigned_value % 3) << std::endl;
    std::cout << "unsigned_value + 5 = " << (unsigned_value + 5)
              << " (same as greater_unsigned_value: " << greater_unsigned_value << ")" << std::endl;

    std::cout << "\n=== Mixed Type Arithmetic with int128 ===" << std::endl;

    constexpr boost::int128::int128 signed_value {-3};
    std::cout << "signed_value = " << signed_value << std::endl;

    std::cout << "signed_value + 1U = " << (signed_value + 1U) << std::endl;
    std::cout << "signed_value - 4U = " << (signed_value - 4U) << std::endl;
    std::cout << "signed_value * 2 = " << (signed_value * 2) << std::endl;
    std::cout << "signed_value / 4U = " << (signed_value / 4U) << std::endl;

    return 0;
}

=== Mixed Type Arithmetic with uint128 ===
unsigned_value = 3
unsigned_value + 1 = 4
unsigned_value - 1 = 2
unsigned_value * 2 = 6
unsigned_value / 3 = 1
unsigned_value % 3 = 0
unsigned_value + 5 = 8 (same as greater_unsigned_value: 8)

=== Mixed Type Arithmetic with int128 ===
signed_value = -3
signed_value + 1U = -2
signed_value - 4U = -7
signed_value * 2 = -6
signed_value / 4U = 0

// Copyright 2025 Matt Borland
// Distributed under the Boost Software License, Version 1.0.
// https://www.boost.org/LICENSE_1_0.txt

#include <boost/int128.hpp>
#include <boost/int128/random.hpp> // Not included in the convenience header, but needed for boost.random interop

#include <boost/math/statistics/univariate_statistics.hpp>
#include <boost/random/uniform_int_distribution.hpp>
#include <iostream>
#include <limits>
#include <array>
#include <random>

int main()
{
    std::cout << "=== uint128 ===" << '\n';

    // Setup our rng and distribution
    std::mt19937_64 rng {42};
    boost::random::uniform_int_distribution<boost::int128::uint128> dist {0, (std::numeric_limits<boost::int128::uint128>::max)()};

    // Create a dataset for ourselves of random uint128s using our dist and rng from above
    std::array<boost::int128::uint128, 10000> data_set;
    for (auto& value : data_set)
    {
        value = dist(rng);
    }

    // Perform some rudimentary statistical analysis on our dataset
    std::cout << "    Mean: " << boost::math::statistics::mean(data_set) << '\n';
    std::cout << "Variance: " << boost::math::statistics::variance(data_set) << '\n';
    std::cout << "  Median: " << boost::math::statistics::median(data_set) << '\n';

    std::cout << "=== int128 ===" << '\n';

    // We can also generate random signed integers using int128
    boost::random::uniform_int_distribution<boost::int128::int128> signed_dist {std::numeric_limits<boost::int128::int128>::min(), std::numeric_limits<boost::int128::int128>::max()};

    std::cout << "Random int128: " << signed_dist(rng) << std::endl;

    return 0;
}

=== uint128 ===
    Mean: 22125900135088040520646253247977468
Variance: 15183108029620265677746188314852225
  Median: 169775281866460752209725324063124732284
=== int128 ===
Random int128: 45422201008201503618595888886744218664

// Copyright 2025 Matt Borland
// Distributed under the Boost Software License, Version 1.0.
// https://www.boost.org/LICENSE_1_0.txt

#include <boost/int128/int128.hpp>
#include <boost/int128/iostream.hpp>
#include <boost/int128/charconv.hpp>
#include <boost/charconv.hpp>
#include <iostream>
#include <limits>
#include <cstring>

int main()
{
    using boost::int128::uint128;
    using boost::int128::int128;

    char buffer[64];

    // === to_chars: Convert integers to character strings ===
    std::cout << "=== to_chars ===" << std::endl;

    // Unsigned 128-bit to decimal string
    constexpr uint128 max_u128 {std::numeric_limits<uint128>::max()};
    auto result {boost::charconv::to_chars(buffer, buffer + sizeof(buffer), max_u128)};
    *result.ptr = '\0';
    std::cout << "uint128 max (decimal): " << buffer << std::endl;

    // Signed 128-bit to decimal string
    constexpr int128 min_i128 {std::numeric_limits<int128>::min()};
    result = boost::charconv::to_chars(buffer, buffer + sizeof(buffer), min_i128);
    *result.ptr = '\0';
    std::cout << "int128 min (decimal):  " << buffer << std::endl;

    // Hexadecimal output (base 16)
    uint128 hex_value {UINT64_C(0xDEADBEEF), UINT64_C(0xCAFEBABE12345678)};
    result = boost::charconv::to_chars(buffer, buffer + sizeof(buffer), hex_value, 16);
    *result.ptr = '\0';
    std::cout << "uint128 (hex): 0x" << buffer << std::endl;

    // Octal output (base 8)
    result = boost::charconv::to_chars(buffer, buffer + sizeof(buffer), int128{511}, 8);
    *result.ptr = '\0';
    std::cout << "int128 511 (octal): 0" << buffer << std::endl;

    // === from_chars: Parse character strings to integers ===
    std::cout << "\n=== from_chars ===" << std::endl;

    // Parse decimal string to uint128
    const char* decimal_str {"340282366920938463463374607431768211455"};
    uint128 parsed_unsigned;
    boost::charconv::from_chars(decimal_str, decimal_str + std::strlen(decimal_str), parsed_unsigned);
    std::cout << "Parsed \"" << decimal_str << "\"" << std::endl;
    std::cout << "  Result: " << parsed_unsigned << std::endl;
    std::cout << "  Equals max? " << std::boolalpha << (parsed_unsigned == max_u128) << std::endl;

    // Parse negative decimal string to int128
    const char* negative_str {"-170141183460469231731687303715884105728"};
    int128 parsed_signed;
    boost::charconv::from_chars(negative_str, negative_str + std::strlen(negative_str), parsed_signed);
    std::cout << "Parsed \"" << negative_str << "\"" << std::endl;
    std::cout << "  Result: " << parsed_signed << std::endl;
    std::cout << "  Equals min? " << (parsed_signed == min_i128) << std::endl;

    // Parse hexadecimal string (base 16)
    const char* hex_str {"DEADBEEFCAFEBABE12345678"};
    uint128 parsed_hex;
    boost::charconv::from_chars(hex_str, hex_str + std::strlen(hex_str), parsed_hex, 16);
    std::cout << "Parsed hex \"" << hex_str << "\"" << std::endl;
    std::cout << "  Result: " << parsed_hex << std::endl;

    return 0;
}

=== to_chars ===
uint128 max (decimal): 340282366920938463463374607431768211455
int128 min (decimal):  -170141183460469231731687303715884105728
uint128 (hex): 0xdeadbeefcafebabe12345678
int128 511 (octal): 0777

=== from_chars ===
Parsed "340282366920938463463374607431768211455"
  Result: 340282366920938463463374607431768211455
  Equals max? true
Parsed "-170141183460469231731687303715884105728"
  Result: -170141183460469231731687303715884105728
  Equals min? true
Parsed hex "DEADBEEFCAFEBABE12345678"
  Result: 68915718020162848918556923512

// Copyright 2026 Matt Borland
// Distributed under the Boost Software License, Version 1.0.
// https://www.boost.org/LICENSE_1_0.txt

// This example demonstrates Boost.ContainerHash integration with int128 types.
// Including <boost/int128/hash.hpp> injects hash_value overloads for the library
// types, so boost::hash, boost::hash_combine, boost::hash_range, and the
// boost::unordered containers (which default to boost::hash) all work with no
// extra configuration.

#include <boost/int128/int128.hpp>
#include <boost/int128/hash.hpp>
#include <boost/container_hash/hash.hpp>
#include <boost/unordered/unordered_map.hpp>
#include <boost/unordered/unordered_flat_map.hpp>
#include <cstddef>
#include <functional>
#include <iostream>
#include <string>

using boost::int128::uint128;
using boost::int128::int128;

// A user-defined composite key that holds 128-bit fields. Providing a hash_value
// overload in the type's own namespace lets Boost.ContainerHash find it via ADL,
// and boost::hash_combine reuses the int128 hashes supplied by hash.hpp.
struct point
{
    int128 x;
    int128 y;
};

bool operator==(const point& lhs, const point& rhs)
{
    return lhs.x == rhs.x && lhs.y == rhs.y;
}

std::size_t hash_value(const point& p)
{
    std::size_t seed {0};
    boost::hash_combine(seed, p.x);
    boost::hash_combine(seed, p.y);
    return seed;
}

int main()
{
    std::cout << "=== boost::hash on 128-bit types ===" << std::endl;

    // boost::hash<T> dispatches to the hash_value overloads from hash.hpp, which
    // delegate to std::hash, so the two functors always agree.
    const uint128 big {UINT64_C(0xDEADBEEF), UINT64_C(0xCAFEBABE12345678)};
    const int128 neg {-123456789012345678LL};

    std::cout << "boost::hash matches std::hash (uint128): "
              << std::boolalpha << (boost::hash<uint128>{}(big) == std::hash<uint128>{}(big)) << std::endl;
    std::cout << "boost::hash matches std::hash (int128):  "
              << (boost::hash<int128>{}(neg) == std::hash<int128>{}(neg)) << std::endl;

    std::cout << "\n=== boost::unordered_map<uint128, ...> ===" << std::endl;

    // boost::unordered_map defaults to boost::hash<Key>, so uint128 keys need
    // no explicit hasher.
    boost::unordered_map<uint128, std::string> labels {};
    labels[uint128{1, 0}] = "two to the sixty-fourth";                   // 2^64
    labels[uint128{UINT64_C(0x8000000000000000), 0}] = "two to the one hundred twenty-seventh"; // 2^127
    labels[uint128{42}] = "forty-two";

    std::cout << "Entries: " << labels.size() << std::endl;
    std::cout << "Label at 2^64: " << labels[uint128{1, 0}] << std::endl;
    std::cout << "Contains 42: " << (labels.find(uint128{42}) != labels.end()) << std::endl;

    std::cout << "\n=== hash_combine for a composite key ===" << std::endl;

    // The point hasher combines two int128 fields; boost::hash<point> finds it
    // via ADL, letting point be used as a key directly.
    boost::unordered_map<point, long> populations {};
    populations[point{10, 20}] = 5000000;
    populations[point{-30, 40}] = 250000;

    std::cout << "Cities stored: " << populations.size() << std::endl;
    std::cout << "Population at (10, 20): " << populations[point{10, 20}] << std::endl;
    std::cout << "Same coordinate hashes equal: "
              << (boost::hash<point>{}(point{10, 20}) == boost::hash<point>{}(point{10, 20})) << std::endl;

    std::cout << "\n=== boost::unordered_flat_map<int128, ...> ===" << std::endl;

    // The modern flat container also defaults to boost::hash.
    boost::unordered_flat_map<int128, int> counts {};
    counts[int128{-1}] = 1;
    counts[int128{0}] = 2;
    counts[int128{1}] = 3;

    std::cout << "Flat map size: " << counts.size() << std::endl;
    std::cout << "counts[-1] = " << counts[int128{-1}] << std::endl;

    return 0;
}

=== boost::hash on 128-bit types ===
boost::hash matches std::hash (uint128): true
boost::hash matches std::hash (int128):  true

=== boost::unordered_map<uint128, ...> ===
Entries: 3
Label at 2^64: two to the sixty-fourth
Contains 42: true

=== hash_combine for a composite key ===
Cities stored: 2
Population at (10, 20): 5000000
Same coordinate hashes equal: true

=== boost::unordered_flat_map<int128, ...> ===
Flat map size: 3
counts[-1] = 1

// Copyright 2025 Matt Borland
// Distributed under the Boost Software License, Version 1.0.
// https://www.boost.org/LICENSE_1_0.txt

#include <boost/int128/int128.hpp>
#include <boost/int128/string.hpp>
#include <boost/int128/literals.hpp>
#include <iostream>
#include <string>
#include <cstdint>

int main()
{
    using boost::int128::uint128;
    using boost::int128::int128;
    using boost::int128::to_string;
    using namespace boost::int128::literals;

    std::cout << "=== to_string with uint128 ===" << std::endl;

    // Compare against std::to_string for values that fit in 64 bits
    constexpr uint128 u_small {UINT64_C(1234567890)};
    const auto u_small_str {to_string(u_small)};
    const auto u_small_std {std::to_string(std::uint64_t{1234567890})};
    std::cout << "uint128 to_string(1234567890): " << u_small_str << std::endl;
    std::cout << "std::to_string(uint64_t 1234567890): " << u_small_std << std::endl;
    std::cout << "Match: " << std::boolalpha << (u_small_str == u_small_std) << std::endl;

    constexpr uint128 u_max64 {UINT64_MAX};
    const auto u_max64_str {to_string(u_max64)};
    const auto u_max64_std {std::to_string(UINT64_MAX)};
    std::cout << "\nuint128 to_string(UINT64_MAX): " << u_max64_str << std::endl;
    std::cout << "std::to_string(UINT64_MAX):      " << u_max64_std << std::endl;
    std::cout << "Match: " << (u_max64_str == u_max64_std) << std::endl;

    // Values beyond 64-bit range
    const auto large_unsigned {340282366920938463463374607431768211455_U128};
    std::cout << "\nuint128 max: " << to_string(large_unsigned) << std::endl;

    std::cout << "\n=== to_string with int128 ===" << std::endl;

    // Compare against std::to_string for values that fit in 64 bits
    constexpr int128 s_negative {-42};
    const auto s_neg_str {to_string(s_negative)};
    const auto s_neg_std {std::to_string(std::int64_t{-42})};
    std::cout << "int128 to_string(-42): " << s_neg_str << std::endl;
    std::cout << "std::to_string(int64_t -42): " << s_neg_std << std::endl;
    std::cout << "Match: " << (s_neg_str == s_neg_std) << std::endl;

    constexpr int128 s_large {INT64_MAX};
    const auto s_large_str {to_string(s_large)};
    const auto s_large_std {std::to_string(INT64_MAX)};
    std::cout << "\nint128 to_string(INT64_MAX): " << s_large_str << std::endl;
    std::cout << "std::to_string(INT64_MAX):     " << s_large_std << std::endl;
    std::cout << "Match: " << (s_large_str == s_large_std) << std::endl;

    // Values beyond 64-bit range.
    // INT128_MIN via a negated numeric literal fails: the sign is not part of the
    // literal, so the positive magnitude 2^127 is read first and rejected as out of
    // range. A string literal keeps the sign with the digits, so it parses cleanly.
    const auto large_negative {"-170141183460469231731687303715884105728"_i128};
    std::cout << "\nint128 min with string literal: " << to_string(large_negative) << std::endl;

    const auto large_negative_c {BOOST_INT128_INT128_C(-170141183460469231731687303715884105728)};
    std::cout << "\nint128 min with INT128_C macro: " << to_string(large_negative_c) << std::endl;

    const auto large_positive {std::numeric_limits<int128>::max()};
    std::cout << "int128 max: " << to_string(large_positive) << std::endl;

    return 0;
}

=== to_string with uint128 ===
uint128 to_string(1234567890): 1234567890
std::to_string(uint64_t 1234567890): 1234567890
Match: true

uint128 to_string(UINT64_MAX): 18446744073709551615
std::to_string(UINT64_MAX):      18446744073709551615
Match: true

uint128 max: 340282366920938463463374607431768211455

=== to_string with int128 ===
int128 to_string(-42): -42
std::to_string(int64_t -42): -42
Match: true

int128 to_string(INT64_MAX): 9223372036854775807
std::to_string(INT64_MAX):     9223372036854775807
Match: true

int128 min with string literal: 0

int128 min with INT128_C macro: -170141183460469231731687303715884105728
int128 max: 170141183460469231731687303715884105727

// Copyright 2025 Matt Borland
// Distributed under the Boost Software License, Version 1.0.
// https://www.boost.org/LICENSE_1_0.txt

// This example demonstrates {fmt} library integration with int128 types.
// Requires {fmt} to be installed: https://github.com/fmtlib/fmt
//
// For C++20 std::format support, use <boost/int128/format.hpp> instead,
// which provides the same formatting capabilities with std::format.

#include <boost/int128/int128.hpp>
#include <boost/int128/fmt_format.hpp>

#include <fmt/format.h>
#include <iostream>

int main()
{
    using boost::int128::int128;
    using boost::int128::uint128;

    std::cout << "=== Basic Formatting ===" << std::endl;

    constexpr uint128 unsigned_value {0xDEADBEEF, 0xCAFEBABE12345678};
    constexpr int128 signed_value {-123456789012345678};

    // Default decimal formatting
    std::cout << fmt::format("Default (decimal): {}", unsigned_value) << std::endl;
    std::cout << fmt::format("Signed value: {}", signed_value) << std::endl;

    std::cout << "\n=== Base Specifiers ===" << std::endl;

    // Different bases: binary, octal, decimal, hex
    constexpr uint128 value {255};
    std::cout << fmt::format("Binary:      {:b}", value) << std::endl;
    std::cout << fmt::format("Octal:       {:o}", value) << std::endl;
    std::cout << fmt::format("Decimal:     {:d}", value) << std::endl;
    std::cout << fmt::format("Hexadecimal: {:x}", value) << std::endl;
    std::cout << fmt::format("Hex (upper): {:X}", value) << std::endl;

    std::cout << "\n=== Alternate Form (Prefixes) ===" << std::endl;

    // Using # for alternate form adds base prefixes
    std::cout << fmt::format("Binary with prefix:  {:#b}", value) << std::endl;
    std::cout << fmt::format("Octal with prefix:   {:#o}", value) << std::endl;
    std::cout << fmt::format("Hex with prefix:     {:#x}", value) << std::endl;
    std::cout << fmt::format("Hex upper prefix:    {:#X}", value) << std::endl;

    std::cout << "\n=== Sign Options ===" << std::endl;

    constexpr int128 positive {42};
    constexpr int128 negative {-42};

    // Sign specifiers: + (always show), - (default), space (space for positive)
    std::cout << fmt::format("Plus sign:  {:+} and {:+}", positive, negative) << std::endl;
    std::cout << fmt::format("Minus only: {} and {}", positive, negative) << std::endl;
    std::cout << fmt::format("Space sign: {: } and {: }", positive, negative) << std::endl;

    std::cout << "\n=== Zero Padding ===" << std::endl;

    // Padding with zeros (no alignment specifier)
    std::cout << fmt::format("8-digit padding:  {:08}", value) << std::endl;
    std::cout << fmt::format("16-digit padding: {:016}", value) << std::endl;

    std::cout << "\n=== Alignment ===" << std::endl;

    // Left, right, and center alignment with default fill (space)
    std::cout << fmt::format("Left align:   '{:<10}'", positive) << std::endl;
    std::cout << fmt::format("Right align:  '{:>10}'", positive) << std::endl;
    std::cout << fmt::format("Center align: '{:^10}'", positive) << std::endl;

    std::cout << "\n=== Alignment with Fill Characters ===" << std::endl;

    // Custom fill characters
    std::cout << fmt::format("Left with *:   '{:*<10}'", positive) << std::endl;
    std::cout << fmt::format("Right with 0:  '{:0>10}'", positive) << std::endl;
    std::cout << fmt::format("Center with -: '{:-^10}'", positive) << std::endl;

    std::cout << "\n=== Alignment with Sign ===" << std::endl;

    // Alignment combined with sign specifiers
    std::cout << fmt::format("Right align +:  '{:>+10}'", positive) << std::endl;
    std::cout << fmt::format("Left align +:   '{:<+10}'", positive) << std::endl;
    std::cout << fmt::format("Center align +: '{:^+11}'", positive) << std::endl;
    std::cout << fmt::format("Right align -:  '{:*>10}'", negative) << std::endl;

    std::cout << "\n=== Alignment with Hex and Prefix ===" << std::endl;

    // Alignment with base specifiers and prefixes
    std::cout << fmt::format("Right align hex:  '{:>10x}'", value) << std::endl;
    std::cout << fmt::format("Left align hex:   '{:<10x}'", value) << std::endl;
    std::cout << fmt::format("Center with prefix: '{:*^#12x}'", value) << std::endl;

    std::cout << "\n=== Large Values ===" << std::endl;

    // Demonstrate with values beyond 64-bit range
    constexpr auto uint_max {std::numeric_limits<uint128>::max()};
    constexpr auto int_min {std::numeric_limits<int128>::min()};

    std::cout << fmt::format("uint128 max: {}", uint_max) << std::endl;
    std::cout << fmt::format("uint128 max (hex): {:#x}", uint_max) << std::endl;
    std::cout << fmt::format("int128 min: {}", int_min) << std::endl;

    std::cout << "\n=== Combined Format Specifiers ===" << std::endl;

    // Combining multiple specifiers
    std::cout << fmt::format("Hex with prefix, uppercase, padded: {:#016X}", unsigned_value) << std::endl;
    std::cout << fmt::format("Signed with plus, padded: {:+020}", signed_value) << std::endl;

    return 0;
}

=== Basic Formatting ===
Default (decimal): 68915718020162848918556923512
Signed value: -123456789012345678

=== Base Specifiers ===
Binary:      11111111
Octal:       377
Decimal:     255
Hexadecimal: ff
Hex (upper): FF

=== Alternate Form (Prefixes) ===
Binary with prefix:  0b11111111
Octal with prefix:   0377
Hex with prefix:     0xff
Hex upper prefix:    0XFF

=== Sign Options ===
Plus sign:  +42 and -42
Minus only: 42 and -42
Space sign:  42 and -42

=== Zero Padding ===
8-digit padding:  00000255
16-digit padding: 0000000000000255

=== Alignment ===
Left align:   '42        '
Right align:  '        42'
Center align: '    42    '

=== Alignment with Fill Characters ===
Left with *:   '42********'
Right with 0:  '0000000042'
Center with -: '----42----'

=== Alignment with Sign ===
Right align +:  '       +42'
Left align +:   '+42       '
Center align +: '    +42    '
Right align -:  '*******-42'

=== Alignment with Hex and Prefix ===
Right align hex:  '        ff'
Left align hex:   'ff        '
Center with prefix: '****0xff****'

=== Large Values ===
uint128 max: 340282366920938463463374607431768211455
uint128 max (hex): 0xffffffffffffffffffffffffffffffff
int128 min: -170141183460469231731687303715884105728

=== Combined Format Specifiers ===
Hex with prefix, uppercase, padded: 0XDEADBEEFCAFEBABE12345678
Signed with plus, padded: -00123456789012345678

// Copyright 2025 Matt Borland
// Distributed under the Boost Software License, Version 1.0.
// https://www.boost.org/LICENSE_1_0.txt

#include <boost/int128/int128.hpp>
#include <boost/int128/cstdlib.hpp>
#include <boost/int128/iostream.hpp>
#include <iostream>

int main()
{
    using boost::int128::uint128;
    using boost::int128::int128;

    std::cout << "=== div() Function ===" << std::endl;
    std::cout << "Returns both quotient and remainder in a single operation" << std::endl;

    // Unsigned division
    std::cout << "\n--- Unsigned Division ---" << std::endl;

    constexpr uint128 dividend {1000000000000000000ULL};
    constexpr uint128 divisor {7};

    auto uresult {boost::int128::div(dividend, divisor)};
    std::cout << dividend << " / " << divisor << " = " << uresult.quot
              << " remainder " << uresult.rem << std::endl;

    // Verify: quot * divisor + rem == dividend
    std::cout << "Verification: " << uresult.quot << " * " << divisor
              << " + " << uresult.rem << " = " << (uresult.quot * divisor + uresult.rem) << std::endl;

    // Large value division
    std::cout << "\n--- Large Value Division ---" << std::endl;

    constexpr uint128 large_dividend {uint128{1} << 100};  // 2^100
    constexpr uint128 large_divisor {uint128{1} << 50};    // 2^50

    auto large_result {boost::int128::div(large_dividend, large_divisor)};
    std::cout << "2^100 / 2^50 = " << large_result.quot
              << " remainder " << large_result.rem << std::endl;

    // Signed division
    std::cout << "\n--- Signed Division ---" << std::endl;

    constexpr int128 signed_dividend {-100};
    constexpr int128 signed_divisor {7};

    auto sresult {boost::int128::div(signed_dividend, signed_divisor)};
    std::cout << signed_dividend << " / " << signed_divisor << " = " << sresult.quot
              << " remainder " << sresult.rem << std::endl;

    // Different sign combinations
    std::cout << "\n--- Sign Combinations ---" << std::endl;

    constexpr int128 pos {17};
    constexpr int128 neg {-17};
    constexpr int128 div_pos {5};
    constexpr int128 div_neg {-5};

    auto pp {boost::int128::div(pos, div_pos)};
    auto pn {boost::int128::div(pos, div_neg)};
    auto np {boost::int128::div(neg, div_pos)};
    auto nn {boost::int128::div(neg, div_neg)};

    std::cout << " 17 /  5 = " << pp.quot << " remainder " << pp.rem << std::endl;
    std::cout << " 17 / -5 = " << pn.quot << " remainder " << pn.rem << std::endl;
    std::cout << "-17 /  5 = " << np.quot << " remainder " << np.rem << std::endl;
    std::cout << "-17 / -5 = " << nn.quot << " remainder " << nn.rem << std::endl;

    // Edge case: dividend smaller than divisor
    std::cout << "\n--- Edge Cases ---" << std::endl;

    auto small_div {boost::int128::div(uint128{3}, uint128{10})};
    std::cout << "3 / 10 = " << small_div.quot << " remainder " << small_div.rem << std::endl;

    return 0;
}

=== div() Function ===
Returns both quotient and remainder in a single operation

--- Unsigned Division ---
1000000000000000000 / 7 = 142857142857142857 remainder 1
Verification: 142857142857142857 * 7 + 1 = 1000000000000000000

--- Large Value Division ---
2^100 / 2^50 = 1125899906842624 remainder 0

--- Signed Division ---
-100 / 7 = -14 remainder -2

--- Sign Combinations ---
 17 /  5 = 3 remainder 2
 17 / -5 = -3 remainder 2
-17 /  5 = -3 remainder -2
-17 / -5 = 3 remainder -2

--- Edge Cases ---
3 / 10 = 0 remainder 3

//  Copyright Matt Borland 2026.
//  Use, modification and distribution are subject to the
//  Boost Software License, Version 1.0. (See accompanying file
//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

#include <iostream>
#include <vector>
#include <random>
#include <limits>
#include <boost/int128.hpp>
#include <boost/int128/numeric.hpp>
#include <boost/int128/random.hpp>
#include <boost/random/uniform_int_distribution.hpp>

#include <cuda_runtime.h>

using test_type = boost::int128::uint128;

// Calculates the GCD of 2 values on device
__global__ void cuda_gcd(const test_type* in1, const test_type* in2, test_type* out, int numElements)
{
    int i = blockDim.x * blockIdx.x + threadIdx.x;

    if (i < numElements)
    {
        out[i] = boost::int128::gcd(in1[i], in2[i]);
    }
}

// Allocate managed space so that the arrays can be used on both host and device
void allocate(test_type** in, int numElements)
{
    cudaError_t err = cudaSuccess;
    err = cudaMallocManaged(in, numElements * sizeof(test_type));
    if (err != cudaSuccess)
    {
        throw std::runtime_error(cudaGetErrorString(err));
    }

    cudaDeviceSynchronize();
}

void cleanup(test_type** in1, test_type** in2, test_type** out)
{
    if (*in1 != nullptr)
    {
        cudaFree(*in1);
        *in1 = nullptr;
    }

    if (*in2 != nullptr)
    {
        cudaFree(*in2);
        *in2 = nullptr;
    }

    if (*out != nullptr)
    {
        cudaFree(*out);
        *out = nullptr;
    }

    cudaDeviceReset();
}

int main()
{
    std::mt19937_64 rng {42};

    const int numElements = 50000;
    std::cout << "[Vector operation on " << numElements << " elements]" << std::endl;

    // Allocate managed space for our inputs and GPU outputs
    // We then fill them with random numbers

    test_type* in1 = nullptr;
    test_type* in2 = nullptr;
    test_type* out = nullptr;

    allocate(&in1, numElements);
    allocate(&in2, numElements);
    allocate(&out, numElements);

    boost::random::uniform_int_distribution<test_type> dist {(std::numeric_limits<test_type>::min)(), (std::numeric_limits<test_type>::max)()};
    for (std::size_t i = 0; i < numElements; ++i)
    {
        in1[i] = dist(rng);
        in2[i] = dist(rng);
    }

    const int threadsPerBlock = 256;
    const int blocksPerGrid = (numElements + threadsPerBlock - 1) / threadsPerBlock;
    std::cout << "CUDA kernel launch with " << blocksPerGrid << " blocks of " << threadsPerBlock << " threads" << std::endl;

    // Launch the CUDA kernel and check for errors

    cuda_gcd<<<blocksPerGrid, threadsPerBlock>>>(in1, in2, out, numElements);
    cudaDeviceSynchronize();

    cudaError_t err = cudaGetLastError();
    if (err != cudaSuccess)
    {
        std::cerr << "Failed to launch kernel (error code " << cudaGetErrorString(err) << ")!" << std::endl;
        cleanup(&in1, &in2, &out);
        return EXIT_FAILURE;
    }

    // We now will perform the same operation using the same inputs on CPU,
    // to compare the results for equality

    std::vector<test_type> results;
    results.reserve(numElements);

    for (int i = 0; i < numElements; ++i)
    {
        results.emplace_back(boost::int128::gcd(in1[i], in2[i]));
    }

    // We can now compare that our operation on GPU and the same operation on CPU have identical results

    for (int i = 0; i < numElements; ++i)
    {
        if (out[i] != results[i])
        {
            std::cerr << "Result verification failed at element: " << i << "!" << std::endl;
            cleanup(&in1, &in2, &out);
            return EXIT_FAILURE;
        }
    }

    cleanup(&in1, &in2, &out);

    std::cout << "All CPU and GPU computed elements match!" << std::endl;

    return 0;
}

[Vector operation on 50000 elements]
CUDA kernel launch with 196 blocks of 256 threads
All CPU and GPU computed elements match!

//  Copyright Matt Borland 2026.
//  Use, modification and distribution are subject to the
//  Boost Software License, Version 1.0. (See accompanying file
//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

#include <sycl/sycl.hpp>
#include <boost/int128.hpp>
#include <boost/int128/numeric.hpp>
#include <iostream>
#include <vector>
#include <random>
#include <cstdint>
#include <cstdlib>

using test_type = boost::int128::uint128;

// Calculates the GCD of two values on the SYCL device and verifies against the host
int main()
{
    std::mt19937_64 rng {42};

    const int numElements {50000};
    std::cout << "[Vector operation on " << numElements << " elements]" << std::endl;

    sycl::queue q;
    std::cout << "SYCL device: " << q.get_device().get_info<sycl::info::device::name>() << std::endl;

    // Allocate shared (USM) memory so the arrays are usable on both host and device
    test_type* in1 {sycl::malloc_shared<test_type>(numElements, q)};
    test_type* in2 {sycl::malloc_shared<test_type>(numElements, q)};
    test_type* out {sycl::malloc_shared<test_type>(numElements, q)};

    for (int i {0}; i < numElements; ++i)
    {
        in1[i] = test_type{rng(), rng()};
        in2[i] = test_type{rng(), rng()};
    }

    // Launch the SYCL kernel: each work item computes one gcd
    q.submit([&](sycl::handler& h)
    {
        h.parallel_for(sycl::range<1>(numElements), [=](sycl::id<1> idx)
        {
            const int i {static_cast<int>(idx[0])};
            out[i] = boost::int128::gcd(in1[i], in2[i]);
        });
    }).wait();

    // Perform the same operation on the host and compare
    std::vector<test_type> results;
    results.reserve(numElements);
    for (int i {0}; i < numElements; ++i)
    {
        results.emplace_back(boost::int128::gcd(in1[i], in2[i]));
    }

    int ret {EXIT_SUCCESS};
    for (int i {0}; i < numElements; ++i)
    {
        if (out[i] != results[i])
        {
            std::cerr << "Result verification failed at element: " << i << "!" << std::endl;
            ret = EXIT_FAILURE;
            break;
        }
    }

    if (ret == EXIT_SUCCESS)
    {
        std::cout << "All CPU and GPU computed elements match!" << std::endl;
    }

    sycl::free(in1, q);
    sycl::free(in2, q);
    sycl::free(out, q);

    return ret;
}

[Vector operation on 50000 elements]
SYCL device: Intel(R) ...
All CPU and GPU computed elements match!

//  Copyright Matt Borland 2026.
//  Use, modification and distribution are subject to the
//  Boost Software License, Version 1.0. (See accompanying file
//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

// Deterministic Miller-Rabin primality testing on the GPU.
//
// Every step of the test reduces a product of two residues modulo n. For a
// 64-bit n both operands can approach 2^64, so the product needs a full 128
// bits before it can be reduced. uint128 supplies that intermediate in
// device code, where unsigned __int128 is not portably available.

#include <boost/int128.hpp>
#include <cstdint>
#include <iostream>
#include <limits>
#include <stdexcept>
#include <vector>

#include <cuda_runtime.h>

using boost::int128::uint128;

// Number of bases that make Miller-Rabin deterministic for every 64-bit n
constexpr int num_bases {12};

// Computes (a * b) % n for 64-bit operands. The 128-bit product is formed
// first so that no high-order bits are discarded before the reduction.
__host__ __device__ std::uint64_t mulmod(const std::uint64_t a, const std::uint64_t b, const std::uint64_t n) noexcept
{
    return static_cast<std::uint64_t>((uint128{a} * uint128{b}) % uint128{n});
}

// Deterministic Miller-Rabin. Exact for the whole 64-bit range.
__host__ __device__ bool is_prime(const std::uint64_t n) noexcept
{
    const std::uint64_t bases[num_bases] {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};

    for (int i {0}; i < num_bases; ++i)
    {
        if (n == bases[i])
        {
            return true;
        }

        if (n % bases[i] == 0U)
        {
            return false;
        }
    }

    if (n < 41U)
    {
        return false;
    }

    // Write n - 1 as d * 2^s with d odd
    std::uint64_t d {n - 1U};
    int s {0};
    while ((d & 1U) == 0U)
    {
        d >>= 1U;
        ++s;
    }

    for (int i {0}; i < num_bases; ++i)
    {
        // powm keeps every intermediate in a uint128, so a^d mod n is exact
        auto x {static_cast<std::uint64_t>(boost::int128::powm(uint128{bases[i]}, uint128{d}, uint128{n}))};

        if (x == 1U || x == n - 1U)
        {
            continue;
        }

        bool composite {true};
        for (int r {1}; r < s; ++r)
        {
            x = mulmod(x, x, n);

            if (x == n - 1U)
            {
                composite = false;
                break;
            }
        }

        if (composite)
        {
            return false;
        }
    }

    return true;
}

// One candidate per thread. Writes the candidate back when it is prime, 0 otherwise.
__global__ void cuda_find_primes(const std::uint64_t first, std::uint64_t* out, const int numCandidates)
{
    const auto i {static_cast<int>(blockDim.x * blockIdx.x + threadIdx.x)};

    if (i < numCandidates)
    {
        const std::uint64_t candidate {first + 2ULL * static_cast<std::uint64_t>(i)};
        out[i] = is_prime(candidate) ? candidate : 0ULL;
    }
}

// Allocate managed space so that the array can be used on both host and device
void allocate(std::uint64_t** out, int numElements)
{
    const cudaError_t err {cudaMallocManaged(out, numElements * sizeof(std::uint64_t))};
    if (err != cudaSuccess)
    {
        throw std::runtime_error(cudaGetErrorString(err));
    }

    cudaDeviceSynchronize();
}

void cleanup(std::uint64_t** out)
{
    if (*out != nullptr)
    {
        cudaFree(*out);
        *out = nullptr;
    }

    cudaDeviceReset();
}

int main()
{
    // A 64-bit multiply is not enough, even though both operands are below n
    const std::uint64_t a {11400714819323198485ULL};
    const std::uint64_t b {14029467366897019727ULL};
    const std::uint64_t n {18446744073709551557ULL};

    std::cout << "Full 128-bit product a * b: " << uint128{a} * uint128{b} << std::endl;
    std::cout << "  (a * b) % n, 64-bit math: " << (a * b) % n << "  <- the product wrapped" << std::endl;
    std::cout << "  (a * b) % n via uint128: " << mulmod(a, b, n) << std::endl;
    std::cout << std::endl;

    // Test the largest odd values that fit in 64 bits
    const int numCandidates {1 << 16};
    const std::uint64_t last {(std::numeric_limits<std::uint64_t>::max)()};
    const std::uint64_t first {last - 2ULL * static_cast<std::uint64_t>(numCandidates - 1)};

    std::cout << "[Miller-Rabin on " << numCandidates << " odd candidates from " << first << " to " << last << "]" << std::endl;

    std::uint64_t* out {nullptr};
    allocate(&out, numCandidates);

    const int threadsPerBlock {256};
    const int blocksPerGrid {(numCandidates + threadsPerBlock - 1) / threadsPerBlock};
    std::cout << "CUDA kernel launch with " << blocksPerGrid << " blocks of " << threadsPerBlock << " threads" << std::endl;

    cuda_find_primes<<<blocksPerGrid, threadsPerBlock>>>(first, out, numCandidates);
    cudaDeviceSynchronize();

    const cudaError_t err {cudaGetLastError()};
    if (err != cudaSuccess)
    {
        std::cerr << "Failed to launch kernel (error code " << cudaGetErrorString(err) << ")!" << std::endl;
        cleanup(&out);
        return EXIT_FAILURE;
    }

    // Run the identical test on the host for comparison
    std::vector<std::uint64_t> results;
    results.reserve(numCandidates);

    for (int i {0}; i < numCandidates; ++i)
    {
        const std::uint64_t candidate {first + 2ULL * static_cast<std::uint64_t>(i)};
        results.emplace_back(is_prime(candidate) ? candidate : 0ULL);
    }

    for (int i {0}; i < numCandidates; ++i)
    {
        if (out[i] != results[i])
        {
            std::cerr << "Result verification failed at element: " << i << "!" << std::endl;
            cleanup(&out);
            return EXIT_FAILURE;
        }
    }

    int found {0};
    std::uint64_t largest {0};
    for (int i {0}; i < numCandidates; ++i)
    {
        if (out[i] != 0U)
        {
            ++found;
            largest = out[i];
        }
    }

    std::cout << "Primes found on device: " << found << std::endl;
    std::cout << "Largest prime below 2^64: " << largest << std::endl;

    cleanup(&out);

    std::cout << "All CPU and GPU computed elements match!" << std::endl;

    return 0;
}

Full 128-bit product a * b: 159945956516994065446882290339781513595
  (a * b) % n, 64-bit math: 17693923505768731003  <- the product wrapped
  (a * b) % n via uint128: 12755662373092413040

[Miller-Rabin on 65536 odd candidates from 18446744073709420545 to 18446744073709551615]
CUDA kernel launch with 256 blocks of 256 threads
Primes found on device: 2879
Largest prime below 2^64: 18446744073709551557
All CPU and GPU computed elements match!

//  Copyright Matt Borland 2026.
//  Use, modification and distribution are subject to the
//  Boost Software License, Version 1.0. (See accompanying file
//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

// Exact, order-independent reduction on the GPU.
//
// A dot product of 64-bit integers overflows a 64-bit accumulator, and a
// floating-point accumulator loses the low-order bits. That second problem is
// worse on a GPU than on a CPU: a parallel reduction visits the elements in
// whatever order the hardware schedules, and floating-point addition is not
// associative, so the answer changes from run to run and from device to
// device. A uint128 accumulator is wide enough to be exact here, and
// integer addition is associative, so the device result is bit-for-bit equal
// to the host result no matter how the work is partitioned.

#include <sycl/sycl.hpp>
#include <boost/int128.hpp>
#include <cstdint>
#include <cstdlib>
#include <iomanip>
#include <iostream>
#include <random>

using boost::int128::uint128;

int main()
{
    const int numElements {1 << 20};
    const int numPartials {1024};

    std::cout << "[Exact dot product of " << numElements << " pairs of 40-bit integers]" << std::endl;

    sycl::queue q;
    std::cout << "SYCL device: " << q.get_device().get_info<sycl::info::device::name>() << std::endl;

    // Shared (USM) allocations are readable from both the host and the device
    std::uint64_t* lhs {sycl::malloc_shared<std::uint64_t>(numElements, q)};
    std::uint64_t* rhs {sycl::malloc_shared<std::uint64_t>(numElements, q)};
    uint128* partials {sycl::malloc_shared<uint128>(numPartials, q)};

    std::mt19937_64 rng {42};
    for (int i {0}; i < numElements; ++i)
    {
        lhs[i] = rng() >> 24U;
        rhs[i] = rng() >> 24U;
    }

    // Each work item accumulates a strided slice into its own 128-bit partial sum
    q.submit([&](sycl::handler& h)
    {
        h.parallel_for(sycl::range<1>(numPartials), [=](sycl::id<1> idx)
        {
            const int p {static_cast<int>(idx[0])};

            uint128 acc {0};
            for (int i {p}; i < numElements; i += numPartials)
            {
                acc += uint128{lhs[i]} * uint128{rhs[i]};
            }

            partials[p] = acc;
        });
    }).wait();

    uint128 device_total {0};
    for (int i {0}; i < numPartials; ++i)
    {
        device_total += partials[i];
    }

    // The same sum on the host, walking the data in both directions
    uint128 host_ascending {0};
    for (int i {0}; i < numElements; ++i)
    {
        host_ascending += uint128{lhs[i]} * uint128{rhs[i]};
    }

    uint128 host_descending {0};
    for (int i {numElements - 1}; i >= 0; --i)
    {
        host_descending += uint128{lhs[i]} * uint128{rhs[i]};
    }

    // What the same reduction produces with a 64-bit accumulator
    std::uint64_t narrow {0};
    for (int i {0}; i < numElements; ++i)
    {
        narrow += lhs[i] * rhs[i];
    }

    std::cout << std::endl;
    std::cout << "uint128 accumulator, exact and " << boost::int128::bit_width(host_ascending) << " bits wide" << std::endl;
    std::cout << "  device, " << numPartials << " strided partials : " << device_total << std::endl;
    std::cout << "  host, ascending order         : " << host_ascending << std::endl;
    std::cout << "  host, descending order        : " << host_descending << std::endl;
    std::cout << "  all three agree, so the partitioning of the reduction does not matter" << std::endl;

    std::cout << std::endl;
    std::cout << "uint64_t accumulator, silently wrong" << std::endl;
    std::cout << "  result                        : " << narrow << std::endl;
    std::cout << "  times it wrapped past 2^64    : " << (host_ascending >> 64U) << std::endl;

    // Widening to double is not a fix either: a single product already needs
    // more significand bits than double has.
    const uint128 exact_product {uint128{lhs[0]} * uint128{rhs[0]}};
    const double rounded_product {static_cast<double>(lhs[0]) * static_cast<double>(rhs[0])};
    const uint128 rounded_as_integer {rounded_product};
    const uint128 rounding_error {rounded_as_integer > exact_product ?
                                    rounded_as_integer - exact_product :
                                    exact_product - rounded_as_integer};

    std::cout << std::endl;
    std::cout << "double, too few significand bits for even one product" << std::endl;
    std::cout << "  lhs[0] * rhs[0] as uint128  : " << exact_product
              << " (" << boost::int128::bit_width(exact_product) << " bits)" << std::endl;
    std::cout << "  lhs[0] * rhs[0] as double     : " << std::setprecision(17) << rounded_product << std::endl;
    std::cout << "  rounding error                : " << rounding_error << std::endl;

    const bool exact {device_total == host_ascending && device_total == host_descending};

    std::cout << std::endl;
    if (exact)
    {
        std::cout << "All CPU and GPU computed elements match!" << std::endl;
    }
    else
    {
        std::cerr << "Result verification failed!" << std::endl;
    }

    sycl::free(lhs, q);
    sycl::free(rhs, q);
    sycl::free(partials, q);

    return exact ? EXIT_SUCCESS : EXIT_FAILURE;
}

[Exact dot product of 1048576 pairs of 40-bit integers]
SYCL device: Intel(R) ...

uint128 accumulator, exact and 99 bits wide
  device, 1024 strided partials : 317117374481213326859021202374
  host, ascending order         : 317117374481213326859021202374
  host, descending order        : 317117374481213326859021202374
  all three agree, so the partitioning of the reduction does not matter

uint64_t accumulator, silently wrong
  result                        : 5008986077640668102
  times it wrapped past 2^64    : 17190967317

double, too few significand bits for even one product
  lhs[0] * rhs[0] as uint128  : 583389027108597783863934 (79 bits)
  lhs[0] * rhs[0] as double     : 5.8338902710859781e+23
  rounding error                : 23584130

All CPU and GPU computed elements match!