Add modular exponentiation function

doc/modules/ROOT/nav.adoc

@@ -27,6 +27,7 @@
 *** xref:api_reference.adoc#api_ios[`<ios>`]
 *** xref:api_reference.adoc#api_numeric[`<numeric>`]
 *** xref:api_reference.adoc#api_string[`<string>`]
+*** xref:api_reference.adoc#api_utilities[Utilities]
 ** xref:api_reference.adoc#api_macros[Macros]
 *** xref:api_reference.adoc#api_macro_literals[Literals]
 *** xref:api_reference.adoc#api_macro_configuration[Configuration]
@@ -55,6 +56,8 @@
 * xref:stream.adoc[]
 * xref:numeric.adoc[]
 * xref:string.adoc[]
+* xref:utilities.adoc[]
+** xref:utilities.adoc#powm[Modular Exponentiation]
 * Benchmarks
 ** xref:u128_benchmarks.adoc[]
 *** xref:u128_benchmarks.adoc#u128_linux[Linux]

doc/modules/ROOT/pages/api_reference.adoc

@@ -257,6 +257,17 @@ Listed by analogous STL header.
 | `std::string` conversion of base-10 values
 |===
 
+[#api_utilities]
+=== xref:utilities.adoc[Utilities]
+
+[cols="1,2", options="header"]
+|===
+| Function | Description
+
+| xref:utilities.adoc#powm[`powm`]
+| Modular exponentiation `(base ^ exp) mod m`
+|===
+
 [#api_macros]
 == Macros
 
@@ -377,4 +388,7 @@ Listed by analogous STL header.
 | `<boost/int128/random.hpp>`
 | Required for usage of Boost.Random
 
+| xref:utilities.adoc[`<boost/int128/utilities.hpp>`]
+| Modular exponentiation and other library-specific utilities
+
 |===

doc/modules/ROOT/pages/utilities.adoc

@@ -0,0 +1,65 @@
+////
+Copyright 2026 Matt Borland
+Distributed under the Boost Software License, Version 1.0.
+https://www.boost.org/LICENSE_1_0.txt
+////
+
+[#utilities]
+= Utilities
+:idprefix: utilities_
+
+The `<boost/int128/utilities.hpp>` header collects helpers that operate on the library types directly and would not fit naturally into the analogous STL-style headers.
+The functions are tuned specifically for `uint128_t` and `int128_t` rather than being template generalizations, which allows the library to dispatch to a fast path based on the shape of the modulus.
+
+[source, c++]
+----
+#include <boost/int128/utilities.hpp>
+----
+
+[#powm]
+== Modular Exponentiation
+
+Computes `(base ^ exp) mod m`.
+The naive expression `pow(base, exp) % m` is unusable for 128-bit inputs because `base ^ exp` overflows almost immediately; `powm` performs the reduction inside the exponentiation loop and selects an algorithm based on the modulus:
+
+* If `has_single_bit(m)` is `true`, modular reduction collapses to a bitmask and no division is performed.
+* If the modulus fits in 64 bits (`m.high == 0`), the loop runs on 64-bit lanes. Each squaring is a single 64x64 -> 128 multiply followed by a 128-by-64 reduction.
+* Otherwise the modulus uses the full 128 bits, and `powm` uses a shift-and-add inner multiply so that no intermediate value ever exceeds 128 bits. This avoids forming the 256-bit product that a naive square-and-multiply implementation would require.
+
+[source, c++]
+----
+namespace boost {
+namespace int128 {
+
+BOOST_INT128_HOST_DEVICE constexpr uint128_t powm(uint128_t base, uint128_t exp, uint128_t m) noexcept;
+
+BOOST_INT128_HOST_DEVICE constexpr int128_t powm(int128_t base, int128_t exp, int128_t m) noexcept;
+
+} // namespace int128
+} // namespace boost
+----
+
+The signed overload returns the non-negative residue in the range `[0, m)`, matching the convention used by `pow(a, b, m)` in Python and most arbitrary-precision libraries.
+Negative bases are reduced before exponentiation; `(std::numeric_limits<int128_t>::min)()` is handled correctly even though its magnitude is not representable in `int128_t`.
+
+=== Special Cases
+
+[cols="1,1", options="header"]
+|===
+| Input | Result
+
+| `m == 0`
+| `0` (consistent with the library's convention for division by zero)
+
+| `m == 1`
+| `0`
+
+| `exp == 0`
+| `1` (including `powm(0, 0, m)`, which follows the conventional definition `0^0 == 1`)
+
+| `base == 0` and `exp > 0`
+| `0`
+
+| Signed overload with `m <= 0` or `exp < 0`
+| `0` (modular exponentiation requires a positive modulus; a negative exponent would require a modular inverse, which this interface does not provide)
+|===

include/boost/int128.hpp

@@ -15,5 +15,6 @@
 #include <boost/int128/climits.hpp>
 #include <boost/int128/cstdlib.hpp>
 #include <boost/int128/string.hpp>
+#include <boost/int128/utilities.hpp>
 
 #endif // BOOST_INT128_HPP

include/boost/int128/utilities.hpp

@@ -0,0 +1,177 @@
+// Copyright 2026 Matt Borland
+// Distributed under the Boost Software License, Version 1.0.
+// https://www.boost.org/LICENSE_1_0.txt
+
+#ifndef BOOST_INT128_UTILITIES_HPP
+#define BOOST_INT128_UTILITIES_HPP
+
+#include <boost/int128/int128.hpp>
+#include <boost/int128/bit.hpp>
+#include <boost/int128/detail/config.hpp>
+
+#ifndef BOOST_INT128_BUILD_MODULE
+
+#include <cstdint>
+
+#endif
+
+namespace boost {
+namespace int128 {
+
+namespace detail {
+
+// Modular addition for 128-bit operands assuming 0 <= a, b < m
+BOOST_INT128_HOST_DEVICE constexpr uint128_t addmod(const uint128_t a, const uint128_t b, const uint128_t m) noexcept
+{
+    const uint128_t s {a + b};
+
+    if (s < a || s >= m)
+    {
+        return s - m;
+    }
+
+    return s;
+}
+
+// Modular multiplication via shift-and-add for the full 128-bit modulus case
+BOOST_INT128_HOST_DEVICE constexpr uint128_t mulmod_shift(uint128_t a, uint128_t b, const uint128_t m) noexcept
+{
+    uint128_t result {0};
+
+    while (b != 0U)
+    {
+        if (static_cast<bool>(b.low & 1U))
+        {
+            result = addmod(result, a, m);
+        }
+
+        a = addmod(a, a, m);
+        b >>= 1;
+    }
+
+    return result;
+}
+
+// Modular multiplication when the modulus fits in 64 bits
+BOOST_INT128_HOST_DEVICE constexpr std::uint64_t mulmod_word(const std::uint64_t a, const std::uint64_t b, const std::uint64_t m) noexcept
+{
+    return ((uint128_t{a} * uint128_t{b}) % uint128_t{m}).low;
+}
+
+} // namespace detail
+
+// Computes (base ^ exp) mod m using fast modular exponentiation with
+// optimizations specific to the boost::int128 library types
+BOOST_INT128_EXPORT BOOST_INT128_HOST_DEVICE constexpr uint128_t powm(uint128_t base, uint128_t exp, const uint128_t m) noexcept
+{
+    if (BOOST_INT128_UNLIKELY(m == 0U))
+    {
+        return uint128_t{0};
+    }
+
+    if (m == 1U)
+    {
+        return uint128_t{0};
+    }
+
+    if (exp == 0U)
+    {
+        return uint128_t{1};
+    }
+
+    base %= m;
+
+    if (base == 0U)
+    {
+        return uint128_t{0};
+    }
+
+    // Power-of-two modulus: reduction is just a bitmask.
+    if (has_single_bit(m))
+    {
+        const uint128_t mask {m - 1U};
+        uint128_t result {1};
+
+        while (exp != 0U)
+        {
+            if (static_cast<bool>(exp.low & 1U))
+            {
+                result = (result * base) & mask;
+            }
+
+            base = (base * base) & mask;
+            exp >>= 1;
+        }
+
+        return result;
+    }
+
+    // Modulus fits in 64 bits: stay in 64-bit lanes.
+    if (m.high == 0U)
+    {
+        const auto mm {m.low};
+        std::uint64_t result {1};
+        auto b {base.low};
+
+        while (exp != 0U)
+        {
+            if (static_cast<bool>(exp.low & 1U))
+            {
+                result = detail::mulmod_word(result, b, mm);
+            }
+
+            b = detail::mulmod_word(b, b, mm);
+            exp >>= 1;
+        }
+
+        return uint128_t{result};
+    }
+
+    // General 128-bit modulus: shift-and-add for each squaring keeps every
+    // intermediate strictly below m without needing a 256-bit product.
+    uint128_t result {1};
+
+    while (exp != 0U)
+    {
+        if (static_cast<bool>(exp.low & 1U))
+        {
+            result = detail::mulmod_shift(result, base, m);
+        }
+
+        base = detail::mulmod_shift(base, base, m);
+        exp >>= 1;
+    }
+
+    return result;
+}
+
+// Signed overload. Returns the non-negative residue in [0, m)
+BOOST_INT128_EXPORT BOOST_INT128_HOST_DEVICE constexpr int128_t powm(const int128_t base, const int128_t exp, const int128_t m) noexcept
+{
+    if (BOOST_INT128_UNLIKELY(m <= 0 || exp < 0))
+    {
+        return int128_t{0};
+    }
+
+    const uint128_t um {static_cast<uint128_t>(m)};
+
+    uint128_t ub {};
+
+    if (base.high < 0)
+    {
+        const uint128_t magnitude {static_cast<uint128_t>(abs(base))};
+        const uint128_t r {magnitude % um};
+        ub = r == 0U ? uint128_t{0} : static_cast<uint128_t>(um - r);
+    }
+    else
+    {
+        ub = static_cast<uint128_t>(base) % um;
+    }
+
+    return static_cast<int128_t>(powm(ub, static_cast<uint128_t>(exp), um));
+}
+
+} // namespace int128
+} // namespace boost
+
+#endif // BOOST_INT128_UTILITIES_HPP

test/Jamfile

@@ -78,6 +78,7 @@ run test_x64_msvc_div.cpp ;
 
 run test_gcd_lcm.cpp ;
 run test_midpoint.cpp ;
+run test_powm.cpp ;
 
 run test_format.cpp ;
 run test_fmt_format.cpp ;
@@ -128,3 +129,4 @@ compile compile_tests/literals_compile.cpp ;
 compile compile_tests/numeric_compile.cpp ;
 compile compile_tests/string_compile.cpp ;
 compile compile_tests/random_compile.cpp ;
+compile compile_tests/utilities_compile.cpp ;

test/compile_tests/utilities_compile.cpp

@@ -0,0 +1,10 @@
+// Copyright 2026 Matt Borland
+// Distributed under the Boost Software License, Version 1.0.
+// https://www.boost.org/LICENSE_1_0.txt
+
+#include <boost/int128/utilities.hpp>
+
+int main()
+{
+    return 0;
+}