Project Euler 396

发表于 分类于 阅读次数： 本文字数： 5.9k 阅读时长 ≈ 5 分钟

Project Euler 396

题目

Weak Goodstein sequence

For any positive integer \(n\), the \(n\text{th}\) weak Goodstein sequence \(\{g_1, g_2, g_3, \dots\}\) is defined as:

  • \(g_1 = n\)

  • for \(k > 1\), \(g_k\) is obtained by writing \(g_{k-1}\) in base \(k\), interpreting it as a base \(k + 1\) number, and subtracting \(1\).

The sequence terminates when \(g_k\) becomes \(0\).

For example, the \(6\text{th}\) weak Goodstein sequence is \(\{6, 11, 17, 25, \dots\}\):

  • \(g_1 = 6\).

  • \(g_2 = 11\) since \(6 = 110_2, 110_3 = 12\), and \(12 - 1 = 11\).

  • \(g_3 = 17\) since \(11 = 102_3, 102_4 = 18\), and \(18 - 1 = 17\).

  • \(g_4 = 25\) since \(17 = 101_4, 101_5 = 26\), and \(26 - 1 = 25\).

and so on.

It can be shown that every weak Goodstein sequence terminates.

Let \(G(n)\) be the number of nonzero elements in the nth weak Goodstein sequence.

It can be verified that \(G(2) = 3, G(4) = 21\) and \(G(6) = 381\).

It can also be verified that \(\sum G(n) = 2517\) for \(1 \le n < 8\).

Find the last \(9\) digits of \(\sum G(n)\) for \(1 \le n < 16\).

解决方案

若一个数在 \(k\) 进制下写作\(\displaystyle{x=\overline{a_ta_{t-1}\dots a_0}_k=\sum_{i=0}^t a_i k^i, 0\le a_i<k,}\)那么“把同一串数字当作 \((k+1)\) 进制解释”就是\(\displaystyle{\overline{a_ta_{t-1}\dots a_0}_{k+1}=\sum_{i=0}^t a_i (k+1)^i.}\)因此单步操作可以看成：\(\displaystyle{x\mapsto \left(\sum_{i=0}^t a_i (k+1)^i\right)-1.}\)

困难在于：每一步的数位 \(a_i\) 是“当前数在当前底 \(k\) 下的表示”，会出现借位与进位，直接跟踪数值会迅速爆炸。

关键技巧是把“弱 Goodstein 过程的长度转化成一个只由二进制位决定的函数迭代系统。

接下来我们介绍一个计算\(G(n)\)的算法。

令\(H(n)=G(n)+3.\)，可见对二进制展开\(\displaystyle{n=\sum_{i\ge 0} b_i2^i, b_i\in\{0,1\},}\)可以从初值 \(3\) 出发，对每个为 \(1\) 的比特位做一次“变换”，最终得到 \(H(n)\)。

为此定义一族函数 \({F_i}\)：

  • 基础层：\(F_0(x)=x+1\)
  • 递归层：\(F_{m+1}(x)=F_m^{\langle x\rangle}(x)\)，其中\(F_m^{\langle 0\rangle}(x)=x,F_m^{\langle t+1\rangle}(x)=F_m(F_m^{\langle t\rangle}(x)).\)，也就是对\(F_m\)进行了\(x\)次迭代。

由定义可直接算出：

  • \(F_1(x)\)：对 \(F_0\) 做 \(x\) 次迭代：\(F_1(x)=F_0^{\langle x\rangle}(x)=x+x=2x\)
  • \(F_2(x)\)：对 \(F_1(y)=2y\) 做 \(x\) 次迭代，从 \(x\) 开始\(F_2(x)=F_1^{\langle x\rangle}(x)=x\cdot 2^x\)
  • \(F_3(x)\)：\(F_3(x)=F_2^{\langle x\rangle}(x), F_2(y)=y\cdot 2^y\)

令初值\(x=3,\)把 \(n\) 的二进制位从低到高扫描，若第 \(i\) 位为 \(1\)，则更新\(x\leftarrow F_i(x).\)扫描完成后，最终得到\(H(n)=x, G(n)=H(n)-3.\)

我们将在最后证明这两个变换出\(G(n)\)的算法和题目完全等价。

因为 \(n<16\)，二进制最多 4 位，所以只会用到 \(F_0,F_1,F_2,F_3\)。

可见，我们可以先直接使用\(F_2,F_1,F_0\)直接算出\(H(m)\)和\(G(m)\)，其中 \(0\le m<8\)：

\[\begin{array}{|c|c|c|} \hline m & H(m) & G(m) \\ \hline 0 & 3 & 0 \\ \hline 1 & 4 & 1 \\ \hline 2 & 6 & 3 \\ \hline 3 & 8 & 5 \\ \hline 4 & 24 & 21 \\ \hline 5 & 64 & 61 \\ \hline 6 & 384 & 381 \\ \hline 7 & 2048 & 2045 \\ \hline \end{array}\]

对 \(8\le n<16\)，最高位 \(2^3\) 必为 \(1\)，因此\(H(n)=F_3(H(n-8)),G(n)=F_3(H(n-8))-3.\)所以\(\displaystyle{\sum_{8\le n<16}H(n)=\sum_{0\le m< 8}F_3(H(m)).}\)

为此，对模数\(M = 10^9=2^9\cdot 5^9\).分别在模 \(2^9\) 与模 \(5^9\) 下计算，再用 CRT 合并。

注意到\(y\ge 9\Longrightarrow 2^y\equiv 0\pmod{2^9}\Longrightarrow F_2(y)=y\cdot 2^y\equiv 0\pmod{2^9}.\)

而在本题中，对任何 \(h\in\{H(m):0\le m<8\}\)，第一次应用 \(F_2\) 后数值至少变成 \(24\)，再应用一次就必然是 \(2^9\) 的倍数。因此对所有这些 \(h\)，\(F_3(h)\equiv 0\pmod{2^9},G(n)=F_3(h)-3\equiv -3\equiv 509\pmod{2^9}.\)

所以 \(8\le n<16\) 的所有 \(G(n)\) 在模 \(2^9\) 下都是同一个值 \(509\)。

令\(m_5=5^9.\)因为 \(\gcd(2,m_5)=1\)，可用欧拉定理，得\(2^y\bmod m_5=2^{y\bmod \varphi(m_5)}\bmod m_5,\varphi(5^9)=4\cdot 5^8.\)

为了能同时保留 \(y\bmod m_5\) 与 \(y\bmod \varphi(m_5)\) 的信息，取\(L=\mathrm{lcm}(m_5,\varphi(m_5))=4\cdot 5^9.\)在模 \(L\) 下迭代\(T(y)=y\cdot 2^y\pmod L\)

即可正确得到 \(y\bmod m_5\)，从而得到 \(F_3(h)\bmod m_5\)。

最后把\(G(n)\bmod 2^9\)（恒为 \(509\)），\(G(n)\bmod 5^9\)用 CRT 合并得到 \(G(n)\bmod M\)，并累加求和即可。

本文最后证明我们介绍的变换\(H(n)\)恰好就是\(G(n)=H(n)-3\)。

对于弱 Goodstein 的一步（从 \(k\) 到 \(k+1\)）是：

  1. 把当前数写成 \(k\) 进制数字串
  2. 把同一串数字当作 \((k+1)\) 进制解释
  3. 减 \(1\)

关键观察：第 2 步只是“换底解释”，数字串本身不变，第 3 步才会真正改变数字串（借位）。

因此我们可以把状态写成：当前基数 \(k\)，当前数字串（即“在基数 \(k\) 下的表示”）。把“在 \((k+1)\) 进制下减一”的借位规则写成一个纯粹的数字串操作。

从基数 \(2\) 开始的弱 Goodstein 序列中，每一步基数加一。

定义：\(B(n)\)：从初值 \(n\) 出发，第一次得到 \(0\) 时的基数。例如 \(n=2\) 时，序列是 \(2,2,1,0\)，得到 \(0\) 时的基数为 \(5\)，所以 \(B(2)=5\)。

那么有一个立即成立的关系：

引理 1： 非零项个数 \(G(n)=B(n)-2\)，因此\(H(n)=G(n)+3=B(n)+1.\)

证明： 初始基数为 \(2\)，每产生一个非零项就会做一次“基数加一”的操作；当第一次产生 \(0\) 时，基数从 \(2\) 增加了 \(G(n)\) 次，所以 \(B(n)=2+G(n)\)，即 \(G(n)=B(n)-2\)。代入 \(H(n)=G(n)+3\) 得 \(H(n)=B(n)+1\)。证毕。

所以你要证明的“\(H(n)\) 的计算规则”，等价于证明目标： \(B(n)+1\) 可以由二进制位通过 \(F_i\) 的迭代得到。

为了描述二进制位的作用，我们先研究纯幂的情况。

对任意 \(k\ge 2\)，定义：\(f_p(k)\)：从基数 \(k\) 开始，初始数是 \(k^p\)，其弱 Goodstein 过程第一次得到 \(0\) 时的基数。

也就是：\(f_p(k)=B_k(k^p),\)其中 \(B_k(\cdot)\) 表示“从基数 \(k\) 开始”的终止基数。

当 \(p=0\) 时，初始数是 \(k^0=1\)。一步操作：写成 \(k\) 进制就是 “\(1\)”，换成 \((k+1)\) 进制仍然是 \(1\)，再减 \(1\) 得 \(0\)。所以立刻有：\(f_0(k)=k+1.\)

接下来我们证明递推\(f_{p+1}(k)=f_p^{(k+1)}(k)\)。这里的 \(f_p^{(t)}\) 表示对函数 \(f_p\) 做 \(t\) 次复合迭代：\(f_p^{(1)}(k)=f_p(k),f_p^{(t+1)}(k)=f_p(f_p^{(t)}(k)).\)

引理 2： 对任意 \(p\ge 0\)，\(f_{p+1}(k)=f_p^{(k+1)}(k).\)

证明： 初始数为 \(k^{p+1}\)，它在 \(k\) 进制下是\(1\underbrace{00\cdots 0}_{p+1}.\)，即有\(p+1\)个\(0\)。

进行第一步（换底到 \(k+1\) 再减一）。换底解释：仍是 \((k+1)^{p+1}\)，随后减一得到：\((k+1)^{p+1}-1\)。而在 \((k+1)\) 进制下，有\((k+1)^{p+1}-1=\underbrace{kk\cdots k}_{p+1}\)，也就是\(p+1\)个\(k\)。

也就是说第一步之后，基数变为 \(k+1\)，数字串变为全 \(k\)。

接下来发生什么？ 在基数递增的过程中，最低位会不断减小；每当最低位从 \(0\) 再减一，就会向高位借位，并把若干低位填成“当前基数减一”，即填成“上一基数”。

因此对于形如\(\underbrace{kk\cdots k}_{p+1}\)的数字串，它的演化具有严格的“分层”结构：

  • 要让最高位从 \(k\) 下降到 \(k-1\)，必须让低 \(p\) 位经历一次“完整的归零过程”；
  • 而“低 \(p\) 位从全 \(k\) 下降到全 \(0\) 并最终触发一次借位”这一段过程，恰好与“从基数 \(k\) 开始处理 \(k^p\)”同构，只是起始基数发生了变化。

更精确地说： 每当高位减少 \(1\)，低 \(p\) 位所经历的过程会使“当前基数”从某个值 \(x\) 变为 \(f_p(x)\)。 而高位初值是 \(k\)，在第一次借位前还要多经历一次“从 \(k\) 到 \(k+1\) 的换底”效应，因此一共需要经历 \(k+1\) 次这样的“低层归零”过程。

于是总基数变化就是对映射 \(f_p\) 连续复合 \(k+1\) 次，从 \(k\) 出发：

\[ f_{p+1}(k)=f_p^{(k+1)}(k). \]

证毕。因此，这个引理是整个结构定理的核心：更高一位的 \(1\) 会迫使低一层过程重复 \(k+1\) 次。

现在定义平移变量：\(x=k+1.\)再定义\(F_p(x)=f_p(x-1)+1.\)那么由 \(f_0(k)=k+1\) 得\(F_0(x)=x+1,\)并且由引理 2：

\[ \begin{aligned} F_{p+1}(x) &=f_{p+1}(x-1)+1 \\ &=f_p^{\langle x\rangle}(x-1)+1. \end{aligned} \]

注意到 \(f_p^{\langle x\rangle}(x-1)+1\) 恰好对应 “\(F_p\) 在点 \(x\) 上迭代 \(x\) 次”（这是因为每次应用 \(f_p\) 都会让参数 \(k\) 变化，而 \(x=k+1\) 同步变化）。

因此得到：\(F_{p+1}(x)=F_p^{\langle x\rangle}(x),\)这正是我们写的递归层定义。

考虑任意 \(n\) 的二进制展开：\(\displaystyle{n=\sum_{i\ge 0} b_i2^i, b_i\in{0,1}.}\)从基数 \(2\) 开始，也就是从 \(x=3\) 开始。

当某一位 \(b_i=1\) 时，它代表在当前基数意义下存在一个“\(k^i\) 的分量”。而前面我们已经证明：处理一个 \(k^i\) 分量，终止基数的变化由 \(f_i\)（等价地由 \(F_i\)）决定。

并且由于弱 Goodstein 的减一借位机制总是从最低位开始，低位的归零必定先发生，从而更高位的“\(1\)”总是在某个新的基数上才真正开始发挥作用，于是这些作用在基数上表现为函数复合：

从 \(x=3\) 出发，按 \(i=0,1,2,\dots\) 递增扫描：

  • 若 \(b_i=0\)，此位不贡献任何“层级操作”，跳过
  • 若 \(b_i=1\)，此位触发一次 \(F_i\) 变换：\(x\leftarrow F_i(x)\)

最终得到的 \(x\) 就是 \(B(n)+1\)。结合引理 1，这就得到了严格结论：

定理： 设 \(x=3\)，按二进制位从低到高扫描，若 \(b_i=1\) 则执行 \(x\leftarrow F_i(x)\)，最后得到\(x=H(n)=G(n)+3.\)

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
from tools import CRT

MOD = 10 ** 9
M2 = 2 ** 9
M5 = 5 ** 9
L = 4 * M5
PHI5 = 4 * 5 ** 8

def pow2_mod_L(e: int) -> int:
    r5 = pow(2, e % PHI5, M5)
    r4 = 0
    r = r5 + ((r4 - r5) & 3) * M5
    return r % L


def F2_L(x: int) -> int:
    return (x % L) * pow2_mod_L(x) % L


def F3_L(x: int) -> int:
    r = x % L
    for _ in range(x):
        r = F2_L(r)
    return r


def G_mod(n: int) -> int:
    temp = 3
    if n & 1:
        temp += 1
    if n & 2:
        temp <<= 1
    if n & 4:
        temp *= 1 << temp

    if n & 8:
        xL = F3_L(temp)
        g5 = (xL % M5 - 3) % M5
        g2 = (M2 - 3) % M2
        return CRT([M2, M5], [g2, g5])

    return (temp - 3) % MOD


s = 0
for n in range(1, 16):
    s = (s + G_mod(n)) % MOD
print(s)