テトレーション（tetration）
テトレーションは余りなじみのない言葉ですがハイパー演算子とクヌースの矢印（高校数学の美しい物語）に分かりやすい説明があります。要はハイパー演算子の$n=4$の場合で、階乗を繰り返す操作になります。
クヌースの矢印を用いると以下のように表されます。
\begin{align}
& 2 \uparrow \uparrow 3 = 2^{2^2} = 2^4 = 16\\
&漸化式にすると \\
& a \uparrow \uparrow n = a ^{a \uparrow \uparrow (n-1)} \ \ \ \cdots (1) \\
& 表されます
\end{align}
実際の数はとてつもなく大きくなってしまうので、この剰余計算を行います。
テトレーションでaが素数の場合
【例題１】$2027 \uparrow \uparrow 2027 $ の下８桁を求めよ
指数部の剰余計算に以下のオイラーの定理（整数論）を使いたいので、まず$a$が素数の場合を考えます。
\begin{align}
& n を自然数，a をn と互いに素な正の整数としたとき\\
&\large a ^{\phi(m)}= 1 \pmod m\\
\end{align}
式（１）にこの定理を適用してプログラムにします。$\phi(m)$はsympyのtotient関数を使いました。
import sympy as sp
def tetration(a,n,M):
if M <= 2 or n == 1:
return pow(a,n,M)
return pow(a,tetration(a,n-1,int(sp.totient(M))),M)
M = 10**8
print(tetration(2027, 2027, M))
# 24536083
このプログラムでは$M$のtotientを次々にとって2になるまで続けるのですが24回再帰呼び出しが行われます。
totientの代わりにカーマイケルのラムダ関数を使う
要はtotient(m)は$a$を何乗したら$\pmod m$で１になるかを求めるのですが、必ずしも最小値ではないようです。そこで巷で話題のカーマイケル数・カーマイケルの定理についてを参考にしてtotientの代わりに「カーマイケルのラムダ関数」を使ってみます。
定義はやや分かりづらいですが、幸いSympyにはカーマイケルのラムダ関数を計算してくれるreduced_totiont関数がありましたのでそれを使うと10回で収束しました。
それほど大きな改善ではありませんが、覚えておくと良いかもしれません。
from sympy.functions.combinatorial.numbers import reduced_totient
M = 10**8
L = [M]
while True:
M = reduced_totient(M)
L.append(M)
if M <=2: break
print(len(L), L)
# 10 [100000000, 5000000, 250000, 12500, 2500, 500, 100, 20, 4, 2]
テトレーションでaが2の場合
ここで、なぜいきなり$a=2$の場合になるかと言うと、アッカーマン関数(Wikikedia)によれば
\begin{align}
& アッカーマン関数はクヌースの矢印表記を用いると、\\
&\large A(m,n) = 2 \uparrow ^{m-2}(n+3)-3\\
& と表せるので例えば \\
& A(4,4) = 2 \uparrow \uparrow7-3 \\
& とテトレーションの形にできるので \\
& 2 \uparrow \uparrow n \pmod m \\
& \\
\end{align}
が計算できればアッカーマン関数の剰余計算もできることになります
【例題１】$2 \uparrow \uparrow 7 $ の下８桁を求めよ
以下の方針で解いて行きます。
- $M=m \times 2^k$というように2のべき乗とそれ以外に分ける
- $m$は2と互いに素なので前と同様に$\phi(m)$を次の法として漸化式を続けその答えを$y$とする
- $2^q \equiv 0 \pmod {2^k} ただしq \ge k$
ここから以下の連立合同式を中国剰余定理を使って解く
\begin{align}
& x = y \pmod m \\
& x = 0 \pmod {2^k}\\
& これを満たすxがM=m \times 2^kを法として存在する
\end{align}
これをプログラムにすると以下になります。中国剰余定理を解くのはsympyのcrtを使っています。また指数部の剰余を求めるのは、今回はあまり差がなかったのでtotientを用いています。
import sympy as sp
def sep2a(M): # M = m * 2^k
fct = sp.factorint(M)
k = 0 if 2 not in fct else fct[2]
return M // 2**k, k
def tet2(n, M):
if n == 0: return 1 % M
m, k = sep2(M) # M = m * 2^k
y = pow(2,tet2(n-1,int(sp.totient(m))), m)
return sp.ntheory.modular.crt([m, 2**k],[y, 0])[0]
M = 10**8
n = 7
print(f"tetration(2,{n},{M}) = {tet2(n, M)}")
# tetration(2,7,100000000) = 8948736
（開発環境：Google Colab）
この考え方は以下のProject Eulerの問題を解くのに役に立ちます。
- Problem 188: Hyperexponentiation
- Problem 282: The Ackermann Function