Power tower
(math/power_tower.hpp)

View this file on GitHub
Last update: 2022-12-02 10:32:05+00:00
Include: #include "math/power_tower.hpp"

Solves ${ a_0 }^{ { a_1 }^{ { \cdots }^{ a_{n-1} } } }\pmod m$ for positive integers $a_i$ and $m$.
Algorithm implementation is based on maspy’s code.

power_tower

template <typename T, typename Promote = unsigned long long>
T power_tower(const std::vector<T> &a, T m);

It returns ${ a_0 }^{ { a_1 }^{ { \cdots }^{ a_{ n-1 } } } }\pmod m$.

Providing first $\min(2\lfloor\lg{m}\rfloor, n)$ elements of $a$ is sufficient to calculate ${ a_0 }^{ { a_1 }^{ {\cdots }^{ a_{n-1} } } }\pmod m$. In other words ${a_0}^{ {a_1}^{ {\cdots}^{a_{n-1} } } } \equiv {a_0}^{ {a_1}^{ {\cdots}^{a_{k-1} } } } \pmod m$ where $k = \min(2\lfloor\lg{m}\rfloor, n)$. Notice that ${a_0}^{ {a_1}^{ {\cdots}^{a_{n-1} } } } \equiv {a_0}^{ {a_1}^{ {\cdots}^{a_{x-1} } } } \pmod m$ is true for any integer $x$ such that $\varphi^x(m) = 1$.

T is an integral type such as int, unsigned int, long long, and unsigned long long. Promote is an integral type that can express values between $0$ and $(2m-1)^2$, inclusive.

There are various methods to calculate Euler’s phi function. Miller-Rabin and Pollard’s rho and precalculated array are possible alternatives. It affects the complexity of power_tower.

Constraints

$0 < a[i]$
$0 < m$

Complexity

$\mathcal{O}\left(\sqrt{m}\right)$

Verified with

Code

#ifndef POWER_TOWER_HPP
#define POWER_TOWER_HPP

#include <cassert>
#include <limits>
#include <type_traits>
#include <vector>

template <typename T> T euler_phi(T n) {
    T phi = n;
    for (T i = 2; i * i <= n; i++) {
        if (n % i == 0) {
            phi -= phi / i;
            while (n % i == 0) {
                n /= i;
            }
        }
    }
    if (n != 1) {
        phi -= phi / n;
    }
    return phi;
}

template <typename T, typename Promote = unsigned long long>
T power_tower(const std::vector<T> &a, T m) {
    static_assert(std::is_integral_v<T>);
    assert(m > 0);
    std::vector<unsigned long long> mod_chain{
        static_cast<unsigned long long>(m)};
    while (mod_chain.back() > 1) {
        const auto phi = euler_phi(mod_chain.back());
        mod_chain.push_back(phi);
    }
    const auto f = [](Promote x, Promote n, Promote mod) {
        Promote r = 1;
        if (x > mod) {
            x = x % mod + mod;
        }
        while (n) {
            if (n & 1) {
                r *= x;
                if (r > mod) {
                    r = r % mod + mod;
                }
            }
            x *= x;
            if (x > mod) {
                x = x % mod + mod;
            }
            n >>= 1;
        }
        return r;
    };
    Promote r = 1;
    const auto k = static_cast<int>(std::min(a.size(), mod_chain.size()));
    for (auto i = k - 1; i >= 0; --i) {
        assert(a[i] > 0);
        r = f(static_cast<Promote>(a[i]), r, mod_chain[i]);
    }
    return static_cast<T>(r % static_cast<Promote>(m));
}

#endif // POWER_TOWER_HPP

#line 1 "math/power_tower.hpp"



#include <cassert>
#include <limits>
#include <type_traits>
#include <vector>

template <typename T> T euler_phi(T n) {
    T phi = n;
    for (T i = 2; i * i <= n; i++) {
        if (n % i == 0) {
            phi -= phi / i;
            while (n % i == 0) {
                n /= i;
            }
        }
    }
    if (n != 1) {
        phi -= phi / n;
    }
    return phi;
}

template <typename T, typename Promote = unsigned long long>
T power_tower(const std::vector<T> &a, T m) {
    static_assert(std::is_integral_v<T>);
    assert(m > 0);
    std::vector<unsigned long long> mod_chain{
        static_cast<unsigned long long>(m)};
    while (mod_chain.back() > 1) {
        const auto phi = euler_phi(mod_chain.back());
        mod_chain.push_back(phi);
    }
    const auto f = [](Promote x, Promote n, Promote mod) {
        Promote r = 1;
        if (x > mod) {
            x = x % mod + mod;
        }
        while (n) {
            if (n & 1) {
                r *= x;
                if (r > mod) {
                    r = r % mod + mod;
                }
            }
            x *= x;
            if (x > mod) {
                x = x % mod + mod;
            }
            n >>= 1;
        }
        return r;
    };
    Promote r = 1;
    const auto k = static_cast<int>(std::min(a.size(), mod_chain.size()));
    for (auto i = k - 1; i >= 0; --i) {
        assert(a[i] > 0);
        r = f(static_cast<Promote>(a[i]), r, mod_chain[i]);
    }
    return static_cast<T>(r % static_cast<Promote>(m));
}