Tetration
(math/tetration.hpp)

• View this file on GitHub
• Last update: 2026-06-24 14:35:02+09:00
• Include: #include "math/tetration.hpp"

Overview

This header computes power towers quickly under a modulus.

For a repeated base, tetration means:

tetration(base, 0) = 1
tetration(base, 1) = base
tetration(base, 2) = base^base
tetration(base, 3) = base^(base^base)

The implementation works for non-prime and non-coprime moduli. It recursively uses Euler’s totient function and lifts the exponent when the true exponent is large enough. Totients are computed with the existing 64-bit Pollard-Rho factorization library.

The ordinary modular exponent convention 0^0 = 1 is used. Therefore, for example, power_tower_mod({2, 0}, mod) is 1 mod mod.

Functions

All functions are in namespace m1une::math. Bases must be non-negative integers. Signed input types are accepted, but negative values are invalid and asserted against in debug builds.

Function signature	Description	Complexity
uint64_t tetration_mod(base, height, mod)	Returns base^^height mod mod.	Totient-chain factorization and logarithmic modular powers
uint64_t tetration_bounded(base, height, limit)	Returns min(base^^height, limit).	O(log limit) effective recursion for base >= 2
uint64_t power_tower_mod(vector<T> bases, mod)	Returns bases[0]^(bases[1]^(...)) mod mod. Empty tower is 1.	Totient-chain factorization and logarithmic modular powers
uint64_t power_tower_bounded(vector<T> bases, limit)	Returns min(tower, limit). Empty tower is 1.	Bounded by the tower length and log limit

mod must be positive. If mod == 1, modular functions return 0.

The bounded functions are useful for comparisons without big integers. For example, tetration_bounded(3, 4, 1000000) returns 1000000, meaning the true value is at least that limit.

Example

#include "math/tetration.hpp"
#include <iostream>
#include <vector>

int main() {
    std::cout << m1une::math::tetration_mod(2ULL, 4, 1000) << '\n';  // 536
    std::cout << m1une::math::tetration_mod(3ULL, 3, 100) << '\n';   // 87

    std::vector<unsigned long long> bases;
    bases.push_back(2);
    bases.push_back(3);
    bases.push_back(4);
    std::cout << m1une::math::power_tower_mod(bases, 1000000007) << '\n';

    if (m1une::math::tetration_bounded(2ULL, 5, 1000000) == 1000000) {
        std::cout << "large\n";
    }
}

Notes

• height = 0 is defined as an empty tower and equals 1.
• An empty vector passed to power_tower_mod or power_tower_bounded also represents an empty tower and equals 1.
• The return type is always uint64_t.
• The modulus may be composite; coprimality between the base and modulus is not required.

Depends on

• 64-bit Prime Factorization (math/prime_factorization.hpp)

Required by

• Math All (math/all.hpp)

Verified with

• verify/math/math_algorithms.test.cpp
• verify/math/tetration.test.cpp

Code

#ifndef M1UNE_MATH_TETRATION_HPP
#define M1UNE_MATH_TETRATION_HPP 1

#include <cassert>
#include <concepts>
#include <cstdint>
#include <type_traits>
#include <vector>

#include "prime_factorization.hpp"

namespace m1une {
namespace math {

namespace tetration_detail {

template <std::integral T>
requires(!std::same_as<std::remove_cv_t<T>, bool>)
uint64_t to_uint64(T value) {
    if constexpr (std::signed_integral<T>) {
        assert(value >= 0);
    }
    return static_cast<uint64_t>(value);
}

inline uint64_t multiply_mod(uint64_t first, uint64_t second, uint64_t mod) {
    return static_cast<uint64_t>(
        static_cast<__uint128_t>(first) * second % mod
    );
}

inline uint64_t pow_mod(uint64_t base, __uint128_t exponent, uint64_t mod) {
    assert(mod >= 1);
    if (mod == 1) return 0;
    base %= mod;
    uint64_t result = 1 % mod;
    while (exponent > 0) {
        if ((exponent & 1) != 0) result = multiply_mod(result, base, mod);
        base = multiply_mod(base, base, mod);
        exponent >>= 1;
    }
    return result;
}

inline uint64_t pow_bounded(uint64_t base, uint64_t exponent, uint64_t limit) {
    if (limit == 0) return 0;
    __uint128_t result = 1;
    for (uint64_t i = 0; i < exponent; i++) {
        result *= base;
        if (result >= limit) return limit;
    }
    return static_cast<uint64_t>(result);
}

inline uint64_t exponent_threshold(uint64_t base, uint64_t limit) {
    assert(base >= 2);
    if (limit <= 1) return 0;

    uint64_t exponent = 0;
    uint64_t value = 1;
    while (value < limit) {
        exponent++;
        if (value > limit / base) return exponent;
        value *= base;
    }
    return exponent;
}

inline uint64_t tetration_bounded_unsigned(uint64_t base, uint64_t height, uint64_t limit) {
    if (limit == 0) return 0;
    if (height == 0) return limit < 1 ? limit : 1;
    if (height == 1) return base < limit ? base : limit;

    if (base == 0) {
        const uint64_t value = (height & 1) == 0 ? 1 : 0;
        return value < limit ? value : limit;
    }
    if (base == 1) return limit < 1 ? limit : 1;

    const uint64_t threshold = exponent_threshold(base, limit);
    const uint64_t exponent = tetration_bounded_unsigned(base, height - 1, threshold);
    if (exponent >= threshold) return limit;
    return pow_bounded(base, exponent, limit);
}

inline uint64_t tetration_mod_unsigned(uint64_t base, uint64_t height, uint64_t mod) {
    assert(mod >= 1);
    if (mod == 1) return 0;
    if (height == 0) return 1 % mod;
    if (height == 1) return base % mod;
    if (base == 0) return (height & 1) == 0 ? 1 % mod : 0;
    if (base == 1) return 1 % mod;

    const uint64_t phi = euler_phi(mod);
    uint64_t reduced_exponent = tetration_mod_unsigned(base, height - 1, phi);
    __uint128_t exponent = reduced_exponent;
    if (tetration_bounded_unsigned(base, height - 1, phi) >= phi) {
        exponent += phi;
    }
    return pow_mod(base, exponent, mod);
}

inline uint64_t power_tower_bounded_unsigned(
    const std::vector<uint64_t>& bases,
    int index,
    uint64_t limit
) {
    if (limit == 0) return 0;
    if (index == int(bases.size())) return limit < 1 ? limit : 1;

    const uint64_t base = bases[index];
    if (index + 1 == int(bases.size())) return base < limit ? base : limit;

    if (base == 0) {
        const uint64_t exponent = power_tower_bounded_unsigned(bases, index + 1, 1);
        const uint64_t value = exponent == 0 ? 1 : 0;
        return value < limit ? value : limit;
    }
    if (base == 1) return limit < 1 ? limit : 1;

    const uint64_t threshold = exponent_threshold(base, limit);
    const uint64_t exponent = power_tower_bounded_unsigned(bases, index + 1, threshold);
    if (exponent >= threshold) return limit;
    return pow_bounded(base, exponent, limit);
}

inline uint64_t power_tower_mod_unsigned(
    const std::vector<uint64_t>& bases,
    int index,
    uint64_t mod
) {
    assert(mod >= 1);
    if (mod == 1) return 0;
    if (index == int(bases.size())) return 1 % mod;
    if (index + 1 == int(bases.size())) return bases[index] % mod;

    const uint64_t phi = euler_phi(mod);
    uint64_t reduced_exponent = power_tower_mod_unsigned(bases, index + 1, phi);
    __uint128_t exponent = reduced_exponent;
    if (power_tower_bounded_unsigned(bases, index + 1, phi) >= phi) {
        exponent += phi;
    }
    return pow_mod(bases[index], exponent, mod);
}

template <std::integral T>
requires(!std::same_as<std::remove_cv_t<T>, bool>)
std::vector<uint64_t> normalize_bases(const std::vector<T>& bases) {
    std::vector<uint64_t> result;
    result.reserve(bases.size());
    for (T base : bases) result.push_back(to_uint64(base));
    return result;
}

}  // namespace tetration_detail

template <std::integral T>
requires(!std::same_as<std::remove_cv_t<T>, bool>)
uint64_t tetration_mod(T base, uint64_t height, uint64_t mod) {
    assert(mod >= 1);
    return tetration_detail::tetration_mod_unsigned(
        tetration_detail::to_uint64(base),
        height,
        mod
    );
}

template <std::integral T>
requires(!std::same_as<std::remove_cv_t<T>, bool>)
uint64_t tetration_bounded(T base, uint64_t height, uint64_t limit) {
    return tetration_detail::tetration_bounded_unsigned(
        tetration_detail::to_uint64(base),
        height,
        limit
    );
}

template <std::integral T>
requires(!std::same_as<std::remove_cv_t<T>, bool>)
uint64_t power_tower_mod(const std::vector<T>& bases, uint64_t mod) {
    assert(mod >= 1);
    std::vector<uint64_t> normalized = tetration_detail::normalize_bases(bases);
    return tetration_detail::power_tower_mod_unsigned(normalized, 0, mod);
}

template <std::integral T>
requires(!std::same_as<std::remove_cv_t<T>, bool>)
uint64_t power_tower_bounded(const std::vector<T>& bases, uint64_t limit) {
    std::vector<uint64_t> normalized = tetration_detail::normalize_bases(bases);
    return tetration_detail::power_tower_bounded_unsigned(normalized, 0, limit);
}

}  // namespace math
}  // namespace m1une

#endif  // M1UNE_MATH_TETRATION_HPP

#line 1 "math/tetration.hpp"



#include <cassert>
#include <concepts>
#include <cstdint>
#include <type_traits>
#include <vector>

#line 1 "math/prime_factorization.hpp"



#include <algorithm>
#line 7 "math/prime_factorization.hpp"
#include <numeric>
#include <utility>
#line 10 "math/prime_factorization.hpp"

namespace m1une {
namespace math {

namespace internal {

inline uint64_t multiply_mod(uint64_t a, uint64_t b, uint64_t mod) {
    return static_cast<uint64_t>(static_cast<unsigned __int128>(a) * b % mod);
}

inline uint64_t power_mod(uint64_t base, uint64_t exponent, uint64_t mod) {
    uint64_t result = 1;
    while (exponent > 0) {
        if (exponent & 1) result = multiply_mod(result, base, mod);
        base = multiply_mod(base, base, mod);
        exponent >>= 1;
    }
    return result;
}

inline uint64_t pollard_random() {
    static uint64_t state = 0x123456789abcdef0ULL;
    state += 0x9e3779b97f4a7c15ULL;
    uint64_t value = state;
    value = (value ^ (value >> 30)) * 0xbf58476d1ce4e5b9ULL;
    value = (value ^ (value >> 27)) * 0x94d049bb133111ebULL;
    return value ^ (value >> 31);
}

}  // namespace internal

inline bool is_prime(uint64_t value) {
    if (value < 2) return false;
    for (uint64_t prime : {2ULL, 3ULL, 5ULL, 7ULL, 11ULL, 13ULL, 17ULL, 19ULL, 23ULL, 29ULL, 31ULL, 37ULL}) {
        if (value % prime == 0) return value == prime;
    }

    uint64_t odd_part = value - 1;
    int power_of_two = 0;
    while ((odd_part & 1) == 0) {
        odd_part >>= 1;
        power_of_two++;
    }

    for (uint64_t base : {2ULL, 325ULL, 9375ULL, 28178ULL, 450775ULL, 9780504ULL, 1795265022ULL}) {
        if (base % value == 0) continue;
        uint64_t x = internal::power_mod(base % value, odd_part, value);
        if (x == 1 || x == value - 1) continue;

        bool composite = true;
        for (int i = 1; i < power_of_two; i++) {
            x = internal::multiply_mod(x, x, value);
            if (x == value - 1) {
                composite = false;
                break;
            }
        }
        if (composite) return false;
    }
    return true;
}

namespace internal {

inline uint64_t pollard_rho(uint64_t value) {
    for (uint64_t prime : {2ULL, 3ULL, 5ULL, 7ULL, 11ULL, 13ULL, 17ULL, 19ULL, 23ULL, 29ULL, 31ULL, 37ULL}) {
        if (value % prime == 0) return prime;
    }

    while (true) {
        const uint64_t constant = pollard_random() % (value - 1) + 1;
        uint64_t y = pollard_random() % (value - 1) + 1;
        uint64_t x = 0;
        uint64_t saved_y = 0;
        uint64_t gcd = 1;
        uint64_t segment_length = 1;

        auto advance = [&](uint64_t current) {
            return static_cast<uint64_t>(
                (static_cast<unsigned __int128>(multiply_mod(current, current, value)) + constant) % value);
        };

        while (gcd == 1) {
            x = y;
            for (uint64_t i = 0; i < segment_length; i++) y = advance(y);

            for (uint64_t offset = 0; offset < segment_length && gcd == 1; offset += 128) {
                saved_y = y;
                uint64_t product = 1;
                const uint64_t block = std::min<uint64_t>(128, segment_length - offset);
                for (uint64_t i = 0; i < block; i++) {
                    y = advance(y);
                    const uint64_t difference = x > y ? x - y : y - x;
                    product = multiply_mod(product, difference, value);
                }
                gcd = std::gcd(product, value);
            }
            segment_length <<= 1;
        }

        if (gcd == value) {
            do {
                saved_y = advance(saved_y);
                const uint64_t difference = x > saved_y ? x - saved_y : saved_y - x;
                gcd = std::gcd(difference, value);
            } while (gcd == 1);
        }
        if (gcd != value) return gcd;
    }
}

inline void factor_recursively(uint64_t value, std::vector<uint64_t>& factors) {
    if (value == 1) return;
    if (is_prime(value)) {
        factors.push_back(value);
        return;
    }
    const uint64_t divisor = pollard_rho(value);
    factor_recursively(divisor, factors);
    factor_recursively(value / divisor, factors);
}

}  // namespace internal

inline std::vector<uint64_t> prime_factors(uint64_t value) {
    assert(value >= 1);
    std::vector<uint64_t> result;
    internal::factor_recursively(value, result);
    std::sort(result.begin(), result.end());
    return result;
}

inline std::vector<std::pair<uint64_t, int>> prime_factorize(uint64_t value) {
    std::vector<uint64_t> factors = prime_factors(value);
    std::vector<std::pair<uint64_t, int>> result;
    for (uint64_t prime : factors) {
        if (result.empty() || result.back().first != prime) {
            result.emplace_back(prime, 1);
        } else {
            result.back().second++;
        }
    }
    return result;
}

inline std::vector<uint64_t> divisors(uint64_t value) {
    std::vector<uint64_t> result = {1};
    for (const auto& factor : prime_factorize(value)) {
        const int current_size = int(result.size());
        uint64_t power = 1;
        for (int exponent = 1; exponent <= factor.second; exponent++) {
            power *= factor.first;
            for (int i = 0; i < current_size; i++) {
                result.push_back(result[i] * power);
            }
        }
    }
    std::sort(result.begin(), result.end());
    return result;
}

inline uint64_t euler_phi(uint64_t value) {
    assert(value >= 1);
    uint64_t result = value;
    for (const auto& factor : prime_factorize(value)) {
        result = result / factor.first * (factor.first - 1);
    }
    return result;
}

inline int mobius(uint64_t value) {
    assert(value >= 1);
    int result = 1;
    for (const auto& factor : prime_factorize(value)) {
        if (factor.second >= 2) return 0;
        result = -result;
    }
    return result;
}

}  // namespace math
}  // namespace m1une


#line 11 "math/tetration.hpp"

namespace m1une {
namespace math {

namespace tetration_detail {

template <std::integral T>
requires(!std::same_as<std::remove_cv_t<T>, bool>)
uint64_t to_uint64(T value) {
    if constexpr (std::signed_integral<T>) {
        assert(value >= 0);
    }
    return static_cast<uint64_t>(value);
}

inline uint64_t multiply_mod(uint64_t first, uint64_t second, uint64_t mod) {
    return static_cast<uint64_t>(
        static_cast<__uint128_t>(first) * second % mod
    );
}

inline uint64_t pow_mod(uint64_t base, __uint128_t exponent, uint64_t mod) {
    assert(mod >= 1);
    if (mod == 1) return 0;
    base %= mod;
    uint64_t result = 1 % mod;
    while (exponent > 0) {
        if ((exponent & 1) != 0) result = multiply_mod(result, base, mod);
        base = multiply_mod(base, base, mod);
        exponent >>= 1;
    }
    return result;
}

inline uint64_t pow_bounded(uint64_t base, uint64_t exponent, uint64_t limit) {
    if (limit == 0) return 0;
    __uint128_t result = 1;
    for (uint64_t i = 0; i < exponent; i++) {
        result *= base;
        if (result >= limit) return limit;
    }
    return static_cast<uint64_t>(result);
}

inline uint64_t exponent_threshold(uint64_t base, uint64_t limit) {
    assert(base >= 2);
    if (limit <= 1) return 0;

    uint64_t exponent = 0;
    uint64_t value = 1;
    while (value < limit) {
        exponent++;
        if (value > limit / base) return exponent;
        value *= base;
    }
    return exponent;
}

inline uint64_t tetration_bounded_unsigned(uint64_t base, uint64_t height, uint64_t limit) {
    if (limit == 0) return 0;
    if (height == 0) return limit < 1 ? limit : 1;
    if (height == 1) return base < limit ? base : limit;

    if (base == 0) {
        const uint64_t value = (height & 1) == 0 ? 1 : 0;
        return value < limit ? value : limit;
    }
    if (base == 1) return limit < 1 ? limit : 1;

    const uint64_t threshold = exponent_threshold(base, limit);
    const uint64_t exponent = tetration_bounded_unsigned(base, height - 1, threshold);
    if (exponent >= threshold) return limit;
    return pow_bounded(base, exponent, limit);
}

inline uint64_t tetration_mod_unsigned(uint64_t base, uint64_t height, uint64_t mod) {
    assert(mod >= 1);
    if (mod == 1) return 0;
    if (height == 0) return 1 % mod;
    if (height == 1) return base % mod;
    if (base == 0) return (height & 1) == 0 ? 1 % mod : 0;
    if (base == 1) return 1 % mod;

    const uint64_t phi = euler_phi(mod);
    uint64_t reduced_exponent = tetration_mod_unsigned(base, height - 1, phi);
    __uint128_t exponent = reduced_exponent;
    if (tetration_bounded_unsigned(base, height - 1, phi) >= phi) {
        exponent += phi;
    }
    return pow_mod(base, exponent, mod);
}

inline uint64_t power_tower_bounded_unsigned(
    const std::vector<uint64_t>& bases,
    int index,
    uint64_t limit
) {
    if (limit == 0) return 0;
    if (index == int(bases.size())) return limit < 1 ? limit : 1;

    const uint64_t base = bases[index];
    if (index + 1 == int(bases.size())) return base < limit ? base : limit;

    if (base == 0) {
        const uint64_t exponent = power_tower_bounded_unsigned(bases, index + 1, 1);
        const uint64_t value = exponent == 0 ? 1 : 0;
        return value < limit ? value : limit;
    }
    if (base == 1) return limit < 1 ? limit : 1;

    const uint64_t threshold = exponent_threshold(base, limit);
    const uint64_t exponent = power_tower_bounded_unsigned(bases, index + 1, threshold);
    if (exponent >= threshold) return limit;
    return pow_bounded(base, exponent, limit);
}

inline uint64_t power_tower_mod_unsigned(
    const std::vector<uint64_t>& bases,
    int index,
    uint64_t mod
) {
    assert(mod >= 1);
    if (mod == 1) return 0;
    if (index == int(bases.size())) return 1 % mod;
    if (index + 1 == int(bases.size())) return bases[index] % mod;

    const uint64_t phi = euler_phi(mod);
    uint64_t reduced_exponent = power_tower_mod_unsigned(bases, index + 1, phi);
    __uint128_t exponent = reduced_exponent;
    if (power_tower_bounded_unsigned(bases, index + 1, phi) >= phi) {
        exponent += phi;
    }
    return pow_mod(bases[index], exponent, mod);
}

template <std::integral T>
requires(!std::same_as<std::remove_cv_t<T>, bool>)
std::vector<uint64_t> normalize_bases(const std::vector<T>& bases) {
    std::vector<uint64_t> result;
    result.reserve(bases.size());
    for (T base : bases) result.push_back(to_uint64(base));
    return result;
}

}  // namespace tetration_detail

template <std::integral T>
requires(!std::same_as<std::remove_cv_t<T>, bool>)
uint64_t tetration_mod(T base, uint64_t height, uint64_t mod) {
    assert(mod >= 1);
    return tetration_detail::tetration_mod_unsigned(
        tetration_detail::to_uint64(base),
        height,
        mod
    );
}

template <std::integral T>
requires(!std::same_as<std::remove_cv_t<T>, bool>)
uint64_t tetration_bounded(T base, uint64_t height, uint64_t limit) {
    return tetration_detail::tetration_bounded_unsigned(
        tetration_detail::to_uint64(base),
        height,
        limit
    );
}

template <std::integral T>
requires(!std::same_as<std::remove_cv_t<T>, bool>)
uint64_t power_tower_mod(const std::vector<T>& bases, uint64_t mod) {
    assert(mod >= 1);
    std::vector<uint64_t> normalized = tetration_detail::normalize_bases(bases);
    return tetration_detail::power_tower_mod_unsigned(normalized, 0, mod);
}

template <std::integral T>
requires(!std::same_as<std::remove_cv_t<T>, bool>)
uint64_t power_tower_bounded(const std::vector<T>& bases, uint64_t limit) {
    std::vector<uint64_t> normalized = tetration_detail::normalize_bases(bases);
    return tetration_detail::power_tower_bounded_unsigned(normalized, 0, limit);
}

}  // namespace math
}  // namespace m1une

Back to top page