Project Euler 282

题目

The Ackermann function

For non-negative integers \(m, n\), the Ackermann function \(A(m, n)\) is defined as follows:

\[A(m,n)= \left \{\begin{aligned} &n+1 & & \text{if } m=0 \\ &A(m-1,1) & & \text{if } m>0\text{ and } n=0\\ &A(m-1,A(m,n-1)) & & \text{if } m>0\text{ and } n>0\\ \end{aligned}\right.\]

For example \(A(1, 0) = 2, A(2, 2) = 7\) and \(A(3, 4) = 125\).

Find \(\displaystyle\sum_{n=0}^6 A(n,n)\) and give your answer mod \(14^8\).

欧拉降幂公式

欧拉降幂公式，用于高效计算\(a^b\% m\)的值（尤其\(b\)很大时）：

\[ a^b\equiv \left \{\begin{aligned} &a^{b\% \varphi(m)} & & \text{if}\quad \gcd(a,m)=1 \\ &a^b & & \text{else if}\quad b<\varphi(m) &\pmod m\\ &a^{b\%\varphi(m)+\varphi(m)} & & \text{else} \end{aligned}\right. \]

解决方案

令\(M=14^8\)。

该页面介绍了阿克曼函数的一部分内容：

\(\begin{aligned} &A(0,n)=n+1;\\ &A(1,n)=n+2\\ &A(2,n)=2n+3;\\ &A(3,n)=2^{n+3}-3;\\ &A(m,n)=2\uparrow^{m-2} (n+3)-3; \end{aligned}\)

那么可以知道

• \(A(5,5)=2\uparrow^3 8-3=2\uparrow^2 (2\uparrow^3 7)-3;\)
• \(A(6,6)=2\uparrow^{4}9-3=2\uparrow^3 (2\uparrow^4 8)-3=2\uparrow^2 (2\uparrow^3 ((2\uparrow^4 8)-1))-3.\)

可以知道，无论是\(2\uparrow^3 7\)还是\(2\uparrow^3 ((2\uparrow^4 8)-1)\)都是非常大的数。我们使用欧拉降幂公式来计算\(2\uparrow^{2} n\% M\)的值。

通过进一步打表发现，当\(n\)超过某个很小的正整数\(n_0\)后，数值\(2\uparrow^{2} n\% M\)的大小是恒定的。（打表程序和表内容如下所示）

for n in range(30):
    print(n, dfs(2, n, mod)[0])

0 1
1 2
2 4
3 16
4 65536
5 804023040
6 1381536256
7 915627008
8 106739712
9 407921152
10 829575168
11 829575168
12 829575168
13 829575168
14 829575168
15 829575168
16 829575168
17 829575168
18 829575168
19 829575168
20 829575168
21 829575168
22 829575168
23 829575168
24 829575168
25 829575168
26 829575168
27 829575168
28 829575168
29 829575168

由此我们可以判断出\(A(5,5)\% M\)的值和\(A(6,6)\% M\)的值是一样的，均为\(2\uparrow^{10}\%M-3\)。由此可以直接计算出相应的答案。

代码

from tools import phi
from math import log

mod = 14 ** 8


def dfs(a: int, f: int, mod: int):
    if f <= 1:
        return pow(a, f, mod), False
    elif mod == 1:
        return 0, True
    else:
        ph = phi(mod)
        val, flag = dfs(a, f - 1, ph)
        if flag or val * log(a) >= log(mod):
            return pow(a, val + ph, mod), True
        else:
            return pow(a, val, mod), False


def cal(n):
    if n == 0:
        return n + 1
    elif n == 1:
        return n + 2
    elif n == 2:
        return 2 * n + 3
    elif n == 3:
        return pow(2, n + 3, mod) - 3
    elif n == 4:
        return dfs(2, n + 3, mod)[0] - 3
    else:
        return dfs(2, 10, mod)[0] - 3


ans = sum(cal(x) for x in range(1, 7)) % mod
print(ans)