tetration_mod.hpp¶
计算固定底数的幂塔模正整数，支持完整 uint64_t 范围的底数、高度和模数，并处理不互质与 \(0^0\)。
\[
\displaystyle a\uparrow\uparrow h \bmod m
\]
Complexity: Time: \(O(\log \mathrm{modulus})\) factorizations and \(O(\log^{2} \mathrm{modulus})\) modular multiplications. Space: \(O(\log \mathrm{modulus})\).
Source: Library Checker
AC 记录：tetration_mod。
Implementation¶
当前头文件，省略 include guard；依赖见 #include。
// SPDX-License-Identifier: Apache-2.0
// Adapted from Library Checker, with changes for libra.
// Sources and license: third_party/README.md.
/// @complexity Time: O(log modulus) factorizations and O(log^2 modulus)
/// modular multiplications. Space: O(log modulus).
#include "noya/math/factorize.hpp"
#include <cassert>
#include <cstdint>
namespace noya {
namespace tetration_mod_detail {
inline std::uint64_t totient(std::uint64_t val) {
std::uint64_t res = val;
std::uint64_t p = 2;
// Bounded trial division avoids factoring the small values in the phi chain.
for (; p <= 1024 && p * p <= val; p += 1 + (p & 1)) {
if (val % p == 0) {
res = res / p * (p - 1);
do { val /= p; } while (val % p == 0);
}
}
if (p * p > val) {
if (val > 1) res = res / val * (val - 1);
} else {
for (auto [q, e] : factorize(val)) res = res / q * (q - 1);
}
return res;
}
inline std::uint64_t power(std::uint64_t bas, unsigned __int128 exp,
std::uint64_t mod, bool &lar) {
if (bas == 1 || exp == 0) {
return 1;
}
std::uint64_t res = 1;
lar = lar || bas >= mod;
bas %= mod;
while (true) {
if (exp & 1) {
unsigned __int128 prd = static_cast<unsigned __int128>(res) * bas;
if (prd >= mod) {
lar = true;
}
res = std::uint64_t(prd % mod);
}
exp >>= 1;
if (exp == 0) {
break;
}
unsigned __int128 squ = static_cast<unsigned __int128>(bas) * bas;
if (squ >= mod) {
lar = true;
}
bas = std::uint64_t(squ % mod);
}
return res;
}
// Residue plus its large-value marker can require 65 bits.
inline unsigned __int128 tower(std::uint64_t bas, std::uint64_t hei,
std::uint64_t mod, bool &lar) {
if (bas == 0) {
return (~hei) & 1ULL;
}
if (mod == 1) {
lar = true;
return 1;
}
if (bas == 1 || hei == 0) {
return 1;
}
if (hei == 1) {
lar = lar || bas >= mod;
return bas % mod + static_cast<unsigned __int128>(lar) * mod;
}
auto exp = hei == 2 ? static_cast<unsigned __int128>(bas)
: tower(bas, hei - 1, totient(mod), lar);
std::uint64_t res = power(bas, exp, mod, lar);
return res + static_cast<unsigned __int128>(lar) * mod;
}
} // namespace tetration_mod_detail
/// @brief Compute a power tower of `hei` copies of `bas` modulo
/// `mod`, with `0^0 = 1`. Recursing through Euler-totient moduli keeps
/// the exponent small; an extra "large" bit distinguishes a reduced exponent
/// from a genuinely small one when the base is not coprime to the modulus.
inline std::uint64_t tetration_mod(std::uint64_t bas, std::uint64_t hei,
std::uint64_t mod) {
assert(mod >= 1);
bool lar = false;
return tetration_mod_detail::tower(bas, hei, mod, lar) % mod;
}
} // namespace noya
// NOYA_TETRATION_MOD_HPP
#ifndef NOYA_TETRATION_MOD_HPP
#define NOYA_TETRATION_MOD_HPP 1
// SPDX-License-Identifier: Apache-2.0
// Adapted from Library Checker, with changes for libra.
// Sources and license: third_party/README.md.
/// @complexity Time: O(log modulus) factorizations and O(log^2 modulus)
/// modular multiplications. Space: O(log modulus).
#include "noya/math/factorize.hpp"
#include <cassert>
#include <cstdint>
namespace noya {
namespace tetration_mod_detail {
inline std::uint64_t totient(std::uint64_t val) {
std::uint64_t res = val;
std::uint64_t p = 2;
// Bounded trial division avoids factoring the small values in the phi chain.
for (; p <= 1024 && p * p <= val; p += 1 + (p & 1)) {
if (val % p == 0) {
res = res / p * (p - 1);
do { val /= p; } while (val % p == 0);
}
}
if (p * p > val) {
if (val > 1) res = res / val * (val - 1);
} else {
for (auto [q, e] : factorize(val)) res = res / q * (q - 1);
}
return res;
}
inline std::uint64_t power(std::uint64_t bas, unsigned __int128 exp,
std::uint64_t mod, bool &lar) {
if (bas == 1 || exp == 0) {
return 1;
}
std::uint64_t res = 1;
lar = lar || bas >= mod;
bas %= mod;
while (true) {
if (exp & 1) {
unsigned __int128 prd = static_cast<unsigned __int128>(res) * bas;
if (prd >= mod) {
lar = true;
}
res = std::uint64_t(prd % mod);
}
exp >>= 1;
if (exp == 0) {
break;
}
unsigned __int128 squ = static_cast<unsigned __int128>(bas) * bas;
if (squ >= mod) {
lar = true;
}
bas = std::uint64_t(squ % mod);
}
return res;
}
// Residue plus its large-value marker can require 65 bits.
inline unsigned __int128 tower(std::uint64_t bas, std::uint64_t hei,
std::uint64_t mod, bool &lar) {
if (bas == 0) {
return (~hei) & 1ULL;
}
if (mod == 1) {
lar = true;
return 1;
}
if (bas == 1 || hei == 0) {
return 1;
}
if (hei == 1) {
lar = lar || bas >= mod;
return bas % mod + static_cast<unsigned __int128>(lar) * mod;
}
auto exp = hei == 2 ? static_cast<unsigned __int128>(bas)
: tower(bas, hei - 1, totient(mod), lar);
std::uint64_t res = power(bas, exp, mod, lar);
return res + static_cast<unsigned __int128>(lar) * mod;
}
} // namespace tetration_mod_detail
/// @brief Compute a power tower of `hei` copies of `bas` modulo
/// `mod`, with `0^0 = 1`. Recursing through Euler-totient moduli keeps
/// the exponent small; an extra "large" bit distinguishes a reduced exponent
/// from a genuinely small one when the base is not coprime to the modulus.
inline std::uint64_t tetration_mod(std::uint64_t bas, std::uint64_t hei,
std::uint64_t mod) {
assert(mod >= 1);
bool lar = false;
return tetration_mod_detail::tower(bas, hei, mod, lar) % mod;
}
} // namespace noya
#endif // NOYA_TETRATION_MOD_HPP
// SPDX-License-Identifier: Apache-2.0
// Adapted from Library Checker, with changes for libra.
// Sources and license: third_party/README.md.
/// @complexity Time: O(log modulus) factorizations and O(log^2 modulus)
/// modular multiplications. Space: O(log modulus).
/// @complexity Time: O(log n) modular multiplications for primality testing;
/// Pollard-rho factorization is expected about O(n^(1/4)).
/// Space: O(log n) recursion and factors.
#include <algorithm>
#include <array>
#include <cassert>
#include <cstdint>
#include <numeric>
#include <utility>
#include <vector>
namespace noya {
namespace factorize_internal {
using u64 = std::uint64_t;
using u128 = unsigned __int128;
inline u64 multiply_mod(u64 a, u64 b, u64 mod) {
return u64(u128(a) * b % mod);
}
inline u64 power_mod(u64 a, u64 exp, u64 mod) {
u64 res = 1;
while (exp > 0) {
if (exp & 1) {
res = multiply_mod(res, a, mod);
}
a = multiply_mod(a, a, mod);
exp >>= 1;
}
return res;
}
inline bool miller_rabin(u64 n) {
if (n < 2) {
return false;
}
for (u64 p :
std::array<u64, 12>{2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37}) {
if (n % p == 0) {
return n == p;
}
}
int shf = __builtin_ctzll(n - 1);
u64 odd = (n - 1) >> shf;
for (u64 bas :
std::array<u64, 7>{2, 325, 9375, 28178, 450775, 9780504, 1795265022}) {
if (bas % n == 0) {
continue;
}
u64 val = power_mod(bas % n, odd, n);
if (val == 1 || val == n - 1) {
continue;
}
bool cmp = true;
for (int i = 1; i < shf; i++) {
val = multiply_mod(val, val, n);
if (val == n - 1) {
cmp = false;
break;
}
}
if (cmp) {
return false;
}
}
return true;
}
inline u64 splitmix64(u64 &st) {
u64 z = (st += 0x9e3779b97f4a7c15ULL);
z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9ULL;
z = (z ^ (z >> 27)) * 0x94d049bb133111ebULL;
return z ^ (z >> 31);
}
inline u64 pollard_rho(u64 n) {
if (n % 2 == 0) {
return 2;
}
if (n % 3 == 0) {
return 3;
}
static u64 st = 0x123456789abcdef0ULL;
while (true) {
u64 y = splitmix64(st) % (n - 1) + 1;
u64 c = splitmix64(st) % (n - 1) + 1;
constexpr u64 blk = 128;
u64 g = 1;
u64 r = 1;
u64 q = 1;
u64 x = 0;
u64 sy = 0;
auto nxt = [&](u64 val) {
val = multiply_mod(val, val, n);
// Add c modulo n without overflowing 64 bits.
return val >= n - c ? val - (n - c) : val + c;
};
while (g == 1) {
x = y;
for (u64 i = 0; i < r; i++) {
y = nxt(y);
}
for (u64 off = 0; off < r && g == 1; off += blk) {
sy = y;
for (u64 i = 0; i < std::min(blk, r - off); i++) {
y = nxt(y);
u64 dif = x > y ? x - y : y - x;
q = multiply_mod(q, dif, n);
}
g = std::gcd(q, n);
}
r <<= 1;
}
if (g == n) {
do {
sy = nxt(sy);
u64 dif = x > sy ? x - sy : sy - x;
g = std::gcd(dif, n);
} while (g == 1);
}
if (g != n) {
return g;
}
}
}
inline void collect_factors(u64 n, std::vector<u64> &res) {
if (n == 1) {
return;
}
if (miller_rabin(n)) {
res.push_back(n);
return;
}
u64 fct = pollard_rho(n);
collect_factors(fct, res);
collect_factors(n / fct, res);
}
} // namespace factorize_internal
/// @brief Deterministic Miller-Rabin primality test for unsigned 64-bit
/// integers.
inline bool is_prime(std::uint64_t n) {
return factorize_internal::miller_rabin(n);
}
/// @brief Return the prime factors of n with multiplicity in increasing order.
inline std::vector<std::uint64_t> prime_factors(std::uint64_t n) {
assert(n >= 1);
std::vector<std::uint64_t> res;
factorize_internal::collect_factors(n, res);
std::sort(res.begin(), res.end());
return res;
}
/// @brief Return the prime factorization of n as (prime, exponent) pairs.
inline std::vector<std::pair<std::uint64_t, int>> factorize(std::uint64_t n) {
std::vector<std::pair<std::uint64_t, int>> res;
for (std::uint64_t p : prime_factors(n)) {
if (res.empty() || res.back().first != p) {
res.emplace_back(p, 1);
} else {
res.back().second++;
}
}
return res;
}
} // namespace noya
#include <cassert>
#include <cstdint>
namespace noya {
namespace tetration_mod_detail {
inline std::uint64_t totient(std::uint64_t val) {
std::uint64_t res = val;
std::uint64_t p = 2;
// Bounded trial division avoids factoring the small values in the phi chain.
for (; p <= 1024 && p * p <= val; p += 1 + (p & 1)) {
if (val % p == 0) {
res = res / p * (p - 1);
do { val /= p; } while (val % p == 0);
}
}
if (p * p > val) {
if (val > 1) res = res / val * (val - 1);
} else {
for (auto [q, e] : factorize(val)) res = res / q * (q - 1);
}
return res;
}
inline std::uint64_t power(std::uint64_t bas, unsigned __int128 exp,
std::uint64_t mod, bool &lar) {
if (bas == 1 || exp == 0) {
return 1;
}
std::uint64_t res = 1;
lar = lar || bas >= mod;
bas %= mod;
while (true) {
if (exp & 1) {
unsigned __int128 prd = static_cast<unsigned __int128>(res) * bas;
if (prd >= mod) {
lar = true;
}
res = std::uint64_t(prd % mod);
}
exp >>= 1;
if (exp == 0) {
break;
}
unsigned __int128 squ = static_cast<unsigned __int128>(bas) * bas;
if (squ >= mod) {
lar = true;
}
bas = std::uint64_t(squ % mod);
}
return res;
}
// Residue plus its large-value marker can require 65 bits.
inline unsigned __int128 tower(std::uint64_t bas, std::uint64_t hei,
std::uint64_t mod, bool &lar) {
if (bas == 0) {
return (~hei) & 1ULL;
}
if (mod == 1) {
lar = true;
return 1;
}
if (bas == 1 || hei == 0) {
return 1;
}
if (hei == 1) {
lar = lar || bas >= mod;
return bas % mod + static_cast<unsigned __int128>(lar) * mod;
}
auto exp = hei == 2 ? static_cast<unsigned __int128>(bas)
: tower(bas, hei - 1, totient(mod), lar);
std::uint64_t res = power(bas, exp, mod, lar);
return res + static_cast<unsigned __int128>(lar) * mod;
}
} // namespace tetration_mod_detail
/// @brief Compute a power tower of `hei` copies of `bas` modulo
/// `mod`, with `0^0 = 1`. Recursing through Euler-totient moduli keeps
/// the exponent small; an extra "large" bit distinguishes a reduced exponent
/// from a genuinely small one when the base is not coprime to the modulus.
inline std::uint64_t tetration_mod(std::uint64_t bas, std::uint64_t hei,
std::uint64_t mod) {
assert(mod >= 1);
bool lar = false;
return tetration_mod_detail::tower(bas, hei, mod, lar) % mod;
}
} // namespace noya
// NOYA_TETRATION_MOD_HPP