Integer Roots and Powers
(math/integer_arithmetic.hpp)

• View this file on GitHub
• Last update: 2026-07-15 01:09:58+09:00
• Include: #include "math/integer_arithmetic.hpp"

Overview

This header provides exact integer roots and powers without converting through floating point.

isqrt(value);
ceil_sqrt(value);
floor_kth_root(value, degree);
ipow(base, exponent);
checked_ipow(base, exponent);

The descriptive aliases floor_sqrt, integer_pow, and checked_integer_pow are also available.

All functions accept standard integral types except bool. Root inputs must be non-negative. The degree passed to floor_kth_root must be positive, and power exponents must be unsigned integers.

Square Root

Function	Result	Complexity
isqrt(value)	$\lfloor\sqrt{\mathrm{value}}\rfloor$	$O(\log \mathrm{value})$
floor_sqrt(value)	Alias of isqrt	$O(\log \mathrm{value})$
ceil_sqrt(value)	$\lceil\sqrt{\mathrm{value}}\rceil$	$O(\log \mathrm{value})$

The implementation compares with division rather than multiplying candidate roots, so it remains correct near the maximum value of the integer type.

K-th Root

Function	Exact signature	Result	Complexity
floor_kth_root	template <std::integral T, std::integral Degree> constexpr T floor_kth_root(T value, Degree degree)	$\lfloor\mathrm{value}^{1/\mathrm{degree}}\rfloor$	$O(\log \mathrm{value})$ integer operations

value must be non-negative and degree must be positive. Both signed and unsigned degree types are accepted. In particular, degree 1 returns value. For positive value, any degree at least the number of value bits returns 1.

Candidate powers are compared to value using a division bound. No intermediate multiplication can overflow, including for std::numeric_limits<std::uint64_t>::max().

Integer Power

Function	Result	Complexity
ipow(base, exponent)	Exact base raised to exponent	$O(\log \mathrm{exponent})$
integer_pow(base, exponent)	Alias of ipow	$O(\log \mathrm{exponent})$
checked_ipow(base, exponent)	The power, or std::nullopt on overflow	$O(\log \mathrm{exponent})$
checked_integer_pow(base, exponent)	Alias of checked_ipow	$O(\log \mathrm{exponent})$

ipow requires the result to fit in the base type and asserts this condition in debug builds. Use checked_ipow when overflow is possible.

As usual, every nonzero integer to exponent zero is one. This library also defines zero to exponent zero as one, which is convenient for binary exponentiation and combinatorial formulas.

Example

#include "math/integer_arithmetic.hpp"

#include <iostream>

int main() {
    std::cout << m1une::math::isqrt(20LL) << "\n";      // 4
    std::cout << m1une::math::ceil_sqrt(20LL) << "\n"; // 5
    std::cout << m1une::math::floor_kth_root(1000ULL, 3) << "\n"; // 10
    std::cout << m1une::math::ipow(3LL, 10U) << "\n";  // 59049

    auto large = m1une::math::checked_ipow(10LL, 19U);
    if (!large) std::cout << "overflow\n";
}

Required by

• Math All (math/all.hpp)
• Modular Kth Root (math/modular_kth_root.hpp)
• Squarefree Count (math/squarefree_count.hpp)
• Sum of Two Squares (math/two_square_sum.hpp)

Verified with

• verify/math/integer_arithmetic.test.cpp
• verify/math/math_algorithms.test.cpp
• verify/math/modular_kth_root.test.cpp
• verify/math/squarefree_count.test.cpp
• verify/math/two_square_sum.test.cpp
• verify/math/two_square_sum.test.cpp

Code

#ifndef M1UNE_MATH_INTEGER_ARITHMETIC_HPP
#define M1UNE_MATH_INTEGER_ARITHMETIC_HPP 1

#include <cassert>
#include <concepts>
#include <limits>
#include <optional>
#include <type_traits>

namespace m1une {
namespace math {

namespace integer_arithmetic_detail {

template <std::integral T>
requires(!std::same_as<std::remove_cv_t<T>, bool>)
constexpr std::optional<T> checked_multiply(T first, T second) {
    constexpr T minimum = std::numeric_limits<T>::min();
    constexpr T maximum = std::numeric_limits<T>::max();

    if constexpr (std::unsigned_integral<T>) {
        if (second != 0 && maximum / second < first) return std::nullopt;
    } else {
        if (0 < first) {
            if (0 < second) {
                if (maximum / second < first) return std::nullopt;
            } else if (second < minimum / first) {
                return std::nullopt;
            }
        } else if (first < 0) {
            if (0 < second) {
                if (first < minimum / second) return std::nullopt;
            } else if (second < maximum / first) {
                return std::nullopt;
            }
        }
    }
    return T(first * second);
}

template <std::unsigned_integral T>
constexpr bool kth_power_leq(T base, unsigned exponent, T limit) {
    assert(exponent > 0);
    if (base <= 1) return base <= limit;

    const T multiplication_limit = limit / base;
    T product = 1;
    for (unsigned i = 0; i < exponent; i++) {
        if (product > multiplication_limit) return false;
        product *= base;
    }
    return true;
}

}  // namespace integer_arithmetic_detail

// Returns floor(sqrt(value)) exactly, without floating-point arithmetic.
template <std::integral T>
requires(!std::same_as<std::remove_cv_t<T>, bool>)
constexpr T isqrt(T value) {
    if constexpr (std::signed_integral<T>) assert(0 <= value);
    if (value <= 1) return value;

    T low = 1;
    T high = value / 2 + 1;
    while (low < high) {
        T middle = low + (high - low + 1) / 2;
        if (middle <= value / middle) {
            low = middle;
        } else {
            high = middle - 1;
        }
    }
    return low;
}

template <std::integral T>
requires(!std::same_as<std::remove_cv_t<T>, bool>)
constexpr T floor_sqrt(T value) {
    return isqrt(value);
}

// Returns ceil(sqrt(value)) exactly, without floating-point arithmetic.
template <std::integral T>
requires(!std::same_as<std::remove_cv_t<T>, bool>)
constexpr T ceil_sqrt(T value) {
    T result = isqrt(value);
    if (result == 0) return 0;
    if (result != 0 && value / result == result && value % result == 0) {
        return result;
    }
    return result + 1;
}

// Returns floor(value^(1 / degree)) exactly, without floating-point arithmetic.
template <std::integral T, std::integral Degree>
requires(
    !std::same_as<std::remove_cv_t<T>, bool>
    && !std::same_as<std::remove_cv_t<Degree>, bool>
)
constexpr T floor_kth_root(T value, Degree degree) {
    if constexpr (std::signed_integral<T>) {
        assert(0 <= value);
        if (value < 0) return T();
    }
    assert(0 < degree);
    if (degree <= 0) return T();
    if (value <= 1 || degree == 1) return value;
    if (degree == 2) return isqrt(value);

    using U = std::make_unsigned_t<T>;
    using UDegree = std::make_unsigned_t<Degree>;
    constexpr int digits = std::numeric_limits<U>::digits;
    const UDegree unsigned_degree = static_cast<UDegree>(degree);
    if (unsigned_degree >= static_cast<UDegree>(digits)) return T(1);
    const unsigned exponent = static_cast<unsigned>(unsigned_degree);
    const U unsigned_value = static_cast<U>(value);

    int bit_width = 0;
    for (U remaining = unsigned_value; remaining != 0; remaining >>= 1) {
        bit_width++;
    }
    const int root_bits =
        (bit_width + static_cast<int>(exponent) - 1) /
        static_cast<int>(exponent);

    U low = 1;
    U high = U(1) << root_bits;
    while (high - low > 1) {
        const U middle = low + (high - low) / 2;
        if (
            integer_arithmetic_detail::kth_power_leq(
                middle, exponent, unsigned_value
            )
        ) {
            low = middle;
        } else {
            high = middle;
        }
    }
    return static_cast<T>(low);
}

// Returns base^exponent, or nullopt when the result does not fit in T.
template <std::integral T, std::unsigned_integral Exponent>
requires(
    !std::same_as<std::remove_cv_t<T>, bool>
    && !std::same_as<std::remove_cv_t<Exponent>, bool>
)
constexpr std::optional<T> checked_ipow(T base, Exponent exponent) {
    T result = 1;
    while (exponent != 0) {
        if (exponent & 1) {
            auto product =
                integer_arithmetic_detail::checked_multiply(result, base);
            if (!product.has_value()) return std::nullopt;
            result = *product;
        }
        exponent >>= 1;
        if (exponent != 0) {
            auto square =
                integer_arithmetic_detail::checked_multiply(base, base);
            if (!square.has_value()) return std::nullopt;
            base = *square;
        }
    }
    return result;
}

template <std::integral T, std::unsigned_integral Exponent>
requires(
    !std::same_as<std::remove_cv_t<T>, bool>
    && !std::same_as<std::remove_cv_t<Exponent>, bool>
)
constexpr std::optional<T> checked_integer_pow(T base, Exponent exponent) {
    return checked_ipow(base, exponent);
}

// Returns base^exponent. The result must be representable by T.
template <std::integral T, std::unsigned_integral Exponent>
requires(
    !std::same_as<std::remove_cv_t<T>, bool>
    && !std::same_as<std::remove_cv_t<Exponent>, bool>
)
constexpr T ipow(T base, Exponent exponent) {
    std::optional<T> result = checked_ipow(base, exponent);
    assert(result.has_value());
    return result.value_or(T());
}

template <std::integral T, std::unsigned_integral Exponent>
requires(
    !std::same_as<std::remove_cv_t<T>, bool>
    && !std::same_as<std::remove_cv_t<Exponent>, bool>
)
constexpr T integer_pow(T base, Exponent exponent) {
    return ipow(base, exponent);
}

}  // namespace math
}  // namespace m1une

#endif  // M1UNE_MATH_INTEGER_ARITHMETIC_HPP

#line 1 "math/integer_arithmetic.hpp"



#include <cassert>
#include <concepts>
#include <limits>
#include <optional>
#include <type_traits>

namespace m1une {
namespace math {

namespace integer_arithmetic_detail {

template <std::integral T>
requires(!std::same_as<std::remove_cv_t<T>, bool>)
constexpr std::optional<T> checked_multiply(T first, T second) {
    constexpr T minimum = std::numeric_limits<T>::min();
    constexpr T maximum = std::numeric_limits<T>::max();

    if constexpr (std::unsigned_integral<T>) {
        if (second != 0 && maximum / second < first) return std::nullopt;
    } else {
        if (0 < first) {
            if (0 < second) {
                if (maximum / second < first) return std::nullopt;
            } else if (second < minimum / first) {
                return std::nullopt;
            }
        } else if (first < 0) {
            if (0 < second) {
                if (first < minimum / second) return std::nullopt;
            } else if (second < maximum / first) {
                return std::nullopt;
            }
        }
    }
    return T(first * second);
}

template <std::unsigned_integral T>
constexpr bool kth_power_leq(T base, unsigned exponent, T limit) {
    assert(exponent > 0);
    if (base <= 1) return base <= limit;

    const T multiplication_limit = limit / base;
    T product = 1;
    for (unsigned i = 0; i < exponent; i++) {
        if (product > multiplication_limit) return false;
        product *= base;
    }
    return true;
}

}  // namespace integer_arithmetic_detail

// Returns floor(sqrt(value)) exactly, without floating-point arithmetic.
template <std::integral T>
requires(!std::same_as<std::remove_cv_t<T>, bool>)
constexpr T isqrt(T value) {
    if constexpr (std::signed_integral<T>) assert(0 <= value);
    if (value <= 1) return value;

    T low = 1;
    T high = value / 2 + 1;
    while (low < high) {
        T middle = low + (high - low + 1) / 2;
        if (middle <= value / middle) {
            low = middle;
        } else {
            high = middle - 1;
        }
    }
    return low;
}

template <std::integral T>
requires(!std::same_as<std::remove_cv_t<T>, bool>)
constexpr T floor_sqrt(T value) {
    return isqrt(value);
}

// Returns ceil(sqrt(value)) exactly, without floating-point arithmetic.
template <std::integral T>
requires(!std::same_as<std::remove_cv_t<T>, bool>)
constexpr T ceil_sqrt(T value) {
    T result = isqrt(value);
    if (result == 0) return 0;
    if (result != 0 && value / result == result && value % result == 0) {
        return result;
    }
    return result + 1;
}

// Returns floor(value^(1 / degree)) exactly, without floating-point arithmetic.
template <std::integral T, std::integral Degree>
requires(
    !std::same_as<std::remove_cv_t<T>, bool>
    && !std::same_as<std::remove_cv_t<Degree>, bool>
)
constexpr T floor_kth_root(T value, Degree degree) {
    if constexpr (std::signed_integral<T>) {
        assert(0 <= value);
        if (value < 0) return T();
    }
    assert(0 < degree);
    if (degree <= 0) return T();
    if (value <= 1 || degree == 1) return value;
    if (degree == 2) return isqrt(value);

    using U = std::make_unsigned_t<T>;
    using UDegree = std::make_unsigned_t<Degree>;
    constexpr int digits = std::numeric_limits<U>::digits;
    const UDegree unsigned_degree = static_cast<UDegree>(degree);
    if (unsigned_degree >= static_cast<UDegree>(digits)) return T(1);
    const unsigned exponent = static_cast<unsigned>(unsigned_degree);
    const U unsigned_value = static_cast<U>(value);

    int bit_width = 0;
    for (U remaining = unsigned_value; remaining != 0; remaining >>= 1) {
        bit_width++;
    }
    const int root_bits =
        (bit_width + static_cast<int>(exponent) - 1) /
        static_cast<int>(exponent);

    U low = 1;
    U high = U(1) << root_bits;
    while (high - low > 1) {
        const U middle = low + (high - low) / 2;
        if (
            integer_arithmetic_detail::kth_power_leq(
                middle, exponent, unsigned_value
            )
        ) {
            low = middle;
        } else {
            high = middle;
        }
    }
    return static_cast<T>(low);
}

// Returns base^exponent, or nullopt when the result does not fit in T.
template <std::integral T, std::unsigned_integral Exponent>
requires(
    !std::same_as<std::remove_cv_t<T>, bool>
    && !std::same_as<std::remove_cv_t<Exponent>, bool>
)
constexpr std::optional<T> checked_ipow(T base, Exponent exponent) {
    T result = 1;
    while (exponent != 0) {
        if (exponent & 1) {
            auto product =
                integer_arithmetic_detail::checked_multiply(result, base);
            if (!product.has_value()) return std::nullopt;
            result = *product;
        }
        exponent >>= 1;
        if (exponent != 0) {
            auto square =
                integer_arithmetic_detail::checked_multiply(base, base);
            if (!square.has_value()) return std::nullopt;
            base = *square;
        }
    }
    return result;
}

template <std::integral T, std::unsigned_integral Exponent>
requires(
    !std::same_as<std::remove_cv_t<T>, bool>
    && !std::same_as<std::remove_cv_t<Exponent>, bool>
)
constexpr std::optional<T> checked_integer_pow(T base, Exponent exponent) {
    return checked_ipow(base, exponent);
}

// Returns base^exponent. The result must be representable by T.
template <std::integral T, std::unsigned_integral Exponent>
requires(
    !std::same_as<std::remove_cv_t<T>, bool>
    && !std::same_as<std::remove_cv_t<Exponent>, bool>
)
constexpr T ipow(T base, Exponent exponent) {
    std::optional<T> result = checked_ipow(base, exponent);
    assert(result.has_value());
    return result.value_or(T());
}

template <std::integral T, std::unsigned_integral Exponent>
requires(
    !std::same_as<std::remove_cv_t<T>, bool>
    && !std::same_as<std::remove_cv_t<Exponent>, bool>
)
constexpr T integer_pow(T base, Exponent exponent) {
    return ipow(base, exponent);
}

}  // namespace math
}  // namespace m1une

Back to top page