tetration
(modulo/tetration.hpp)

• View this file on GitHub
• Last update: 2020-12-05 07:59:51+09:00
• Include: #include "modulo/tetration.hpp"

Tetration

$a \uparrow \uparrow b \mod m$を$\mathrm{O}(\sqrt{M})$で計算するライブラリ。$a \uparrow \uparrow b$とは、

\[a \uparrow \uparrow b =\begin{cases} 1 & \mathrm{if}\ \ b = 0\newline a^{a \uparrow \uparrow (b-1)} & \mathrm{otherwise}\end{cases}\]

として定められる二項演算である。

概要

(この内容はkoprickyさんのライブラリを参考に書きました。)

$a\uparrow \uparrow b$と$a\uparrow \uparrow (b-1)$の関係式を元に再帰的に計算する。

$a$と$m$が互いに素の時は、オイラーの定理$a ^ {\phi(m)}\equiv 1 \mod m$より、

\[a \uparrow \uparrow b \equiv a^{a \uparrow \uparrow (b-1)} \equiv a^{(a \uparrow \uparrow (b-1) \mod \phi(m))} \mod m\]

という関係式を得ることができる。

次に$a$と$m$が互いに素でない場合を考える。$m=kl$を$l$と$a$が互いに素で$k$が$a$の素因数の積で表されるように取る。

この時、$a^i \mod k$は$i$が十分に大きいときに$0$に等しくなり、$a^i \mod l$はオイラーの定理より$\phi(l)$周期を持つ。よって中国剰余定理と合わせて$a^i \mod m$は$i$が十分に大きい時に$\phi(m)$周期になることがわかる。

ここで具体的に$m=kl$を具体的に計算する必要はないことに留意する。というのも、トーシェント関数の乗法性より$a^i \mod m$は$\phi(m)$周期を持つためである。また、$\lfloor\log_2 m\rfloor \leq \phi(m)$ より$a^{\phi(m)} \equiv 0 \mod k$になることから、$i$が$\phi(m)$未満の時は愚直に計算して、それより$i$が大きいときは周期性を利用して計算すればよい。(実装の際は$a\uparrow \uparrow b$が$m$以上かどうかをフラグとして持っておけばよい。)

以上より、$a\uparrow \uparrow b \mod m$は$a\uparrow \uparrow (b-1) \mod \phi(m)$の値から導出できることがわかった。計算量はアルゴリズム全体でのトーシェント関数の計算量に依存するので、計算量を求めるために数列$m,\phi(m),\phi(\phi(m))\ldots$の性質を考察する。

ここで、$m$が偶数の時$\phi(m)$は$\frac{m}{2}$以下になる。また、$m$が奇数の時、$\phi(m)$は偶数なので$\phi(\phi(m))$は$\frac{\phi(m)}{2}$以下になる。以上より、上記の数列は少なくとも$2$個おきに半減して、$\mathrm{O}(\log_2 m)$回で$1$に収束する。よって全体の計算量は$T(m) = \mathrm{O}(\sqrt{m})+T(\frac{m}{2})$を解いて$\mathrm{O}(\sqrt{m})$となる。

Verified with

• verify/verify-yosupo-math/yosupo-tetration-mod.test.cpp

Code

#pragma once



template <typename T = uint32_t, typename U = uint64_t, T MAX = 1000000000>
T tetration(uint64_t a, uint64_t b, uint64_t m) {
  auto prime_table = [](T m) -> vector<int> {
    T thres = sqrt(m) + 1;
    vector<bool> flg(thres + 1);
    for (T i = 2; i * i <= thres; ++i) {
      if (!flg[i])
        for (T j = i * i; j <= thres; j += i) flg[j] = 1;
    }
    vector<int> ret;
    for (T i = 2; i <= thres; i++)
      if (!flg[i]) ret.push_back(i);
    return ret;
  };

  static const vector<int> ps = prime_table(MAX);

  auto totient = [&](T n) -> T {
    T ret = n;
    for (auto& p : ps) {
      if (p * p > n) break;
      if (n % p == 0) {
        ret = ret / p * (p - 1);
        while (n % p == 0) n /= p;
      }
    }
    if (n != 1) ret = ret / n * (n - 1);
    return ret;
  };

  auto mpow = [](U a, U p, T m) -> pair<T, int> {
    U ret = 1 % m;
    int flg = true;
    for (; p; p >>= 1) {
      if (p & 1)
        if ((ret *= a) >= m) flg = 0, ret %= m;
      if (p == 1) break;
      if ((a *= a) >= m) flg = 0, a %= m;
    }
    return {ret, flg};
  };

  auto calc = [&](auto rec, T a, U b, T m) -> pair<T, int> {
    if (a == 0) return {!(b & 1), true};
    if (a == 1) return {1, true};
    if (m == 1) return {0, false};
    if (b == 0) return {1, true};
    if (b == 1) return {a % m, a < m};
    T phi_m = totient(m), pre, ret;
    int flg1, flg2;
    tie(pre, flg1) = rec(rec, a, b - 1, phi_m);
    tie(ret, flg2) = mpow(a % m, U(pre) + (flg1 ? 0 : phi_m), m);
    return {ret, flg1 && flg2};
  };

  return calc(calc, a, b, m).first % m;
}

/**
 * @brief tetration
 * @docs docs/modulo/tetration.md
 */

#line 2 "modulo/tetration.hpp"



template <typename T = uint32_t, typename U = uint64_t, T MAX = 1000000000>
T tetration(uint64_t a, uint64_t b, uint64_t m) {
  auto prime_table = [](T m) -> vector<int> {
    T thres = sqrt(m) + 1;
    vector<bool> flg(thres + 1);
    for (T i = 2; i * i <= thres; ++i) {
      if (!flg[i])
        for (T j = i * i; j <= thres; j += i) flg[j] = 1;
    }
    vector<int> ret;
    for (T i = 2; i <= thres; i++)
      if (!flg[i]) ret.push_back(i);
    return ret;
  };

  static const vector<int> ps = prime_table(MAX);

  auto totient = [&](T n) -> T {
    T ret = n;
    for (auto& p : ps) {
      if (p * p > n) break;
      if (n % p == 0) {
        ret = ret / p * (p - 1);
        while (n % p == 0) n /= p;
      }
    }
    if (n != 1) ret = ret / n * (n - 1);
    return ret;
  };

  auto mpow = [](U a, U p, T m) -> pair<T, int> {
    U ret = 1 % m;
    int flg = true;
    for (; p; p >>= 1) {
      if (p & 1)
        if ((ret *= a) >= m) flg = 0, ret %= m;
      if (p == 1) break;
      if ((a *= a) >= m) flg = 0, a %= m;
    }
    return {ret, flg};
  };

  auto calc = [&](auto rec, T a, U b, T m) -> pair<T, int> {
    if (a == 0) return {!(b & 1), true};
    if (a == 1) return {1, true};
    if (m == 1) return {0, false};
    if (b == 0) return {1, true};
    if (b == 1) return {a % m, a < m};
    T phi_m = totient(m), pre, ret;
    int flg1, flg2;
    tie(pre, flg1) = rec(rec, a, b - 1, phi_m);
    tie(ret, flg2) = mpow(a % m, U(pre) + (flg1 ? 0 : phi_m), m);
    return {ret, flg1 && flg2};
  };

  return calc(calc, a, b, m).first % m;
}

/**
 * @brief tetration
 * @docs docs/modulo/tetration.md
 */

Back to top page