Project Euler 308

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

<– more –>

Project Euler 308

题目

An amazing Prime-generating Automaton

A program written in the programming language Fractran consists of a list of fractions.

The internal state of the Fractran Virtual Machine is a positive integer, which is initially set to a seed value. Each iteration of a Fractran program multiplies the state integer by the first fraction in the list which will leave it an integer.

For example, one of the Fractran programs that John Horton Conway wrote for prime-generation consists of the following \(14\) fractions:

\(\dfrac{17}{91}, \dfrac{78}{85}, \dfrac{19}{51}, \dfrac{23}{38}, \dfrac{29}{33}, \dfrac{77}{29}, \dfrac{95}{23}, \dfrac{77}{19}, \dfrac{1}{17}, \dfrac{11}{13}, \dfrac{13}{11}, \dfrac{15}{2}, \dfrac{1}{7}, \dfrac{55}{1}\)

Starting with the seed integer \(2\), successive iterations of the program produce the sequence:

\(15, 825, 725, 1925, 2275, 425, \dots, 68, \mathbf{4}, 30, \dots, 136, \mathbf{8}, 60, \dots, 544, \mathbf{32}, 240, \dots\)

The powers of \(2\) that appear in this sequence are \(2^2, 2^3, 2^5, \dots\)

It can be shown that all the powers of \(2\) in this sequence have prime exponents and that all the primes appear as exponents of powers of \(2\), in proper order!

If someone uses the above Fractran program to solve Project Euler Problem 7 (find the \(10001^\text{st}\) prime), how many iterations would be needed until the program produces \(2^{10001\text{st prime}}\) ?

解决方案

FRACTRAN 的关键技巧在于：把当前整数写成素因子分解\(\displaystyle{N = \prod_{i=1}^t p_i^{e_i}.}\)乘以一个分数等价于对若干 \(e_i\) 做加减；而“能否整除”只取决于分母所需的素因子指数是否足够。

因此我们把指数统一存进长度为 10 的数组e[10]={e[0],e[1],...,e[9]}

因此，每一步按顺序找第一条满足条件的规则，执行其更新。

  1. \(\dfrac{17}{91}=\dfrac{17}{7\cdot 13}\)

    • 条件：e[3] >= 1 and e[5] >= 1
    • 更新：e[3] -= 1, e[5] -= 1, e[6] += 1

  2. \(\dfrac{78}{85}=\dfrac{2\cdot 3\cdot 13}{5\cdot 17}\)

    • 条件：e[2] >= 1 and e[6] >= 1
    • 更新：e[2] -= 1, e[6] -= 1, e[0] += 1, e[1] += 1, e[5] += 1

  3. \(\dfrac{19}{51}=\dfrac{19}{3\cdot 17}\)

    • 条件：e[1] >= 1 and e[6] >= 1
    • 更新：e[1] -= 1, e[6] -= 1, e[7] += 1

  4. \(\dfrac{23}{38}=\dfrac{23}{2\cdot 19}\)

    • 条件：e[0] >= 1 and e[7] >= 1
    • 更新：e[0] -= 1, e[7] -= 1, e[8] += 1

  5. \(\dfrac{29}{33}=\dfrac{29}{3\cdot 11}\)

    • 条件：e[1] >= 1 and e[4] >= 1
    • 更新：e[1] -= 1, e[4] -= 1, e[9] += 1

  6. \(\dfrac{77}{29}=\dfrac{7\cdot 11}{29}\)

    • 条件：e[9] >= 1
    • 更新：e[9] -= 1, e[3] += 1, e[4] += 1

  7. \(\dfrac{95}{23}=\dfrac{5\cdot 19}{23}\)

    • 条件：e[8] >= 1
    • 更新：e[8] -= 1, e[2] += 1, e[7] += 1

  8. \(\dfrac{77}{19}=\dfrac{7\cdot 11}{19}\)

    • 条件：e[7] >= 1
    • 更新：e[7] -= 1, e[3] += 1, e[4] += 1

  9. \(\dfrac{1}{17}\)

    • 条件：e[6] >= 1
    • 更新：e[6] -= 1

  10. \(\dfrac{11}{13}\)

  • 条件：e[5] >= 1
  • 更新：e[5] -= 1, e[4] += 1
  1. \(\dfrac{13}{11}\)
  • 条件：e[4] >= 1
  • 更新：e[4] -= 1, e[5] += 1
  1. \(\dfrac{15}{2}=\dfrac{3\cdot 5}{2}\)
  • 条件：e[0] >= 1
  • 更新：e[0] -= 1, e[1] += 1, e[2] += 1
  1. \(\dfrac{1}{7}\)
  • 条件：e[3] >= 1
  • 更新：e[3] -= 1
  1. \(\dfrac{55}{1}=5\cdot 11\)
  • 条件：True（永远可用）
  • 更新：e[2] += 1, e[4] += 1

初始种子 \(N=2\) 对应[1,0,0,0,0,0,0,0,0,0].因此“当前是否为纯 2 的幂” 等价于 e[1]=e[2]=...e[9]==0.

这是一份模拟这个过程的 python 代码：

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
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
from typing import List, Tuple

# p[0..9] = (2,3,5,7,11,13,17,19,23,29)
# e[0..9] = exponent vector for these primes

RULES = [
    # 1) 17/91 = 17/(7*13)
    ({3: 1, 5: 1}, {3: -1, 5: -1, 6: 1}),

    # 2) 78/85 = (2*3*13)/(5*17)
    ({2: 1, 6: 1}, {2: -1, 6: -1, 0: 1, 1: 1, 5: 1}),

    # 3) 19/51 = 19/(3*17)
    ({1: 1, 6: 1}, {1: -1, 6: -1, 7: 1}),

    # 4) 23/38 = 23/(2*19)
    ({0: 1, 7: 1}, {0: -1, 7: -1, 8: 1}),

    # 5) 29/33 = 29/(3*11)
    ({1: 1, 4: 1}, {1: -1, 4: -1, 9: 1}),

    # 6) 77/29 = (7*11)/29
    ({9: 1}, {9: -1, 3: 1, 4: 1}),

    # 7) 95/23 = (5*19)/23
    ({8: 1}, {8: -1, 2: 1, 7: 1}),

    # 8) 77/19 = (7*11)/19
    ({7: 1}, {7: -1, 3: 1, 4: 1}),

    # 9) 1/17
    ({6: 1}, {6: -1}),

    # 10) 11/13
    ({5: 1}, {5: -1, 4: 1}),

    # 11) 13/11
    ({4: 1}, {4: -1, 5: 1}),

    # 12) 15/2 = (3*5)/2
    ({0: 1}, {0: -1, 1: 1, 2: 1}),

    # 13) 1/7
    ({3: 1}, {3: -1}),

    # 14) 55/1 = 5*11  (always applicable)
    ({}, {2: 1, 4: 1}),
]


def step_fractran(e: List[int]) -> int:
    """
    Execute exactly ONE FRACTRAN iteration on exponent vector e[10],
    applying the first applicable rule among RULES.

    Returns:
        idx (int): which rule (1..14) was applied.
    """
    for idx, (need, delta) in enumerate(RULES, start=1):
        ok = True
        for i, req in need.items():
            if e[i] < req:
                ok = False
                break
        if not ok:
            continue

        for i, dv in delta.items():
            e[i] += dv
        return idx

    raise RuntimeError("No rule applicable (should be impossible because last rule is always applicable).")


def run_fractran(steps: int, e0: List[int] | None = None) -> Tuple[List[int], List[int]]:
    if e0 is None:
        e = [1, 0, 0, 0, 0, 0, 0, 0, 0, 0]  # seed N=2
    else:
        if len(e0) != 10:
            raise ValueError("e0 must have length 10.")
        e = e0[:]

    applied: List[int] = []
    for _ in range(steps):
        applied.append(step_fractran(e))
    return e, applied


def is_power_of_two_state(e: List[int]) -> bool:
    return all(v == 0 for v in e[1:])


if __name__ == "__main__":
    e, applied = run_fractran(50)
    print("final e =", e)
    print("applied rules =", applied)
    print("is power of 2?", is_power_of_two_state(e))

关键观察：我们可以发现，

  • 含 e[4], ..., e[9] 的分母（或分子）只在少数规则里出现；
  • 而 e[0], ..., e[3] 几乎贯穿所有“搬运/计数/循环”规则，是典型的“数据寄存器”。

并且模拟运行代码查看结果后发现，e[4], ..., e[9]不会超过 2。因此，e[4], ..., e[9] 很像“状态位/标签”；而 e[0], ..., e[3] 很像“数值/计数器”。考虑从这个角度进行分析。

在可达轨道上，e[8]和e[9] 是“瞬时态”，必成对消去。因此，压缩成 FSM 的关键，就是两条“不会被别的规则插队”的引理：一旦 e[9]>0，下一步一定用规则 6 把 e[9] 消掉；一旦 e[8]>0，下一步一定用规则 7 把 e[8] 消掉。

首先证明如果e[9]>0，下一步一定用规则 6 把 e[9] 消掉。要证明“下一步一定选 6”，只需证明规则 1–5 在这一时刻都不可用。在本题的可达轨道里，e[9] 的产生只能来自规则 5。当我们第一次看到 e[9] 从 \(0\) 变为 \(1\) 时，这一步一定是刚执行完规则 5。执行规则 5 的更新是而在可达轨道中，规则 5 被触发的场景是稳定态 \(S11\)（后面会定义），即执行规则 5 之前满足e[4]=1,e[5]=e[6]=e[7]=e[8]=0。因此执行规则 5 之后立刻得到 e[9]=1,e[4]=0,e[5]=e[6]=e[7]=e[8]=0。因此，根据这些条件，我们检查规则 1-5，可以发现它们并不可用。因此，规则 5 执行完（一旦 e[9]>0）立会跟着规则 6。

首先证明如果e[8]>0，下一步一定用规则 7 把 e[8] 消掉。可见，在本题可达轨道里，e[8] 的产生只能来自规则 4：而在可达轨道中，规则 4 被触发的场景是稳定态 \(S19\)（后面会定义），即执行规则 4 之前满足e[7]=1,e[4]=e[5]=e[6]=0.因此执行完规则 4 后立刻得到e[8]=1,e[7]=0,e[4]=e[5]=e[6]=0。因此，根据这些条件，我们检查规则 1-6，可以发现它们并不可用。因此，规则 4 执行完（一旦 e[9]>0）立会跟着规则 7。

综上所述，规则 5+6 可以一个 2 步的宏，并吸收掉规则 6：

  • 条件：e[1] >= 1 and e[4] >= 1
  • 更新：e[1] -= 1, e[3] += 1

规则 4+7 可以一个 2 步的宏，并吸收掉规则 7：

  • 条件：e[0] >= 1 and e[7] >= 1
  • 更新：e[0] -= 1, e[2] += 1

通过合并宏后，我们可以仅针对各类状态，建立一个状态机。这个状态机只包含 5 类状态，并忽略e[8], e[9]的存在。

状态 \(S0\)（无标记：e[4]=e[5]=e[6]=e[7]=0）

那么可见，目前规则 1–11 都不可用；可用的只有 12、13、14。

一开始，我们连续执行规则 12 共 e[0] 次，每次更新：e[0] -= 1, e[1] += 1, e[2] += 1。完成后，e[0] = 0，此时规则 12 不可用。一共进行了e[0]次步骤。

接下来，连续执行规则 13 共 e[3] 次：每次更新 e[3] -= 1。执行完后 e[3] = 0。一共进行了e[3]次步骤。

最后，执行一次规则 14，更新e[2] += 1, e[4] += 1。一共进行了 1 次步骤。

那么现在的状态标记变成了：e[4]=1,e[5]=e[6]=e[7]=0，进入了状态\(S11\)。

那么从\(S0\)状态转移可以总结成：e'[1] = e[1] + e[0]; e'[2] = e[2] + e[0] + 1; e'[0] = e'[3] = 0, e'[4] = 1，并进入了状态\(S11\)，一共花费了e[0] + e[3] + 1个步骤。

状态 \(S11\)（e[4]=1）

可见，目前规则 1–4 不可用；会被选中的规则只有 5+6 和 11。

一开始，我们连续执行规则 5+6 共 e[1] 次，每次更新：e[1] -= 1, e[3] += 1。完成后，e[1] = 0，此时规则 5+6 不可用。一共进行了e[1] * 2次步骤。

接下来，执行 1 次规则 11 共更新 e[4] -= 1, e[5] += 1。一共进行了 1 次步骤。

那么现在的状态标记变成了：e[5]=1,e[4]=e[6]=e[7]=0，进入了状态\(S13\)。

那么从\(S11\)状态转移可以总结成：e'[3] = e[3] + e[1]; e'[1] = e'[4] = 0, e'[5] = 1，并进入了状态\(S13\)，一共花费了e[1] * 2 + 1个步骤。

状态 \(S13\)（e[5]=1）

当 e[2]>0 and e[3]>0 时，会发生确定的 2 步规则：

首先执行规则 1 的内容：e[3] -= 1, e[5] -= 1, e[6] += 1

因为这时e[6] = 1，所以规则 2 紧接着执行：e[2] -= 1, e[6] -= 1,e[0] += 1, e[1] += 1, e[5] += 1。

因此，规则 1+2 发生的变化如下：e[2] -= 1, e[3] -= 1, e[0] += 1, e[1] += 1。并且e[5] = 1, e[6] = 0依旧保持不变。

令m = min(e[2], e[3])，当进行m次规则 1+2 后，结果就变成：e[2] -= m, e[3] -= m, e[0] += m, e[1] += m

也就是说，如果e[3] == 0，那么规则 1 不可用，下一条可用的是规则 10，更新：e[5] -= 1, e[4] += 1，那么贡献了 1 步数，并进入状态\(S11\)。

否则，如果e[3] > 0，那么此时必有 e[2]=0。规则 2 不可用，但规则 1 仍可用。执行一次规则 1：e[3] -= 1, e[5] -= 1, e[6] = 1，那么贡献了 1 步数，并进入状态\(S17\)。

那么从\(S13\)状态转移可以总结成：

  • 如果 e[3] <= e[2]，那么有e'[3] = 0, e'[2] = e[2] - e[3], e'[0] = e[0] + e[3], e'[1] = e[1] + e[3], e'[5] = 0, e'[4] = 1，总共花费e[3] * 2 + 1的步数。并进入状态\(S11\)。

  • 如果 e[3] > e[2]，那么有e'[3] = e[3] - e[2] - 1, e'[2] = 0, e'[0] = e[0] + e[2], e'[1] = e[1] + e[2], e'[5] = 0, e'[6] = 1，总共花费e[2] * 2 + 1的步数。并进入状态\(S17\)。

状态 \(S17\)（e[6]=1）

在可达轨道上进入 \(S17\) 时满足 e[2]=0（因为它来自 \(S13\) 的分支 2，且该分支恰是 e[2]=0 才发生），因此规则 2 这时不会插队。此时只有 3 和 9 规则可用，即只有两种情况：

  • 若 e[1]>0，规则 3 可用，那么有e[1] -= 1, e[6] = 0, e[7] = 1，花费 1 步数，进入\(S19\)。
  • 若 e[1]==0，规则 9 可用，那么有e[6] = 0，花费 1 步数，进入\(S0\)。

状态 \(S19\)（e[7]=1）

当 e[0]>0 时，可以使用规则 4+7 进行处理。每执行一次宏，步数加 2，净效应：e[0] -= 1, e[2] += 1。我们连续执行规则 4+7 共 e[0] 次。完成后，e[2] += e[0], e[0] = 0。一共进行了e[0] * 2次步骤。

此时规则 8 成为最先可用，更新后，有e[7]=0,e[3]+=1,e[4]+=1,步数加 1，进入 \(S11\)。

那么从\(S19\)状态转移可以总结成：e'[2] = e[2] + e[0]; e'[0] = 0; e'[7] = 0, e'[3] = e[3] + 1, e'[4] = 1，并进入了状态\(S11\)，一共花费了e[0] * 2 + 1个步骤。

那么这时代码可以加速简化成这样子：

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
from sympy import prime

def batch_step(e):
    a, b, c, d, s11, s13, s17, s19 = e  # 2,3,5,7,11,13,17,19

    if s11:  # S11
        return 2 * b + 1, [a, 0, c, d + b, 0, 1, 0, 0]

    if s13:  # S13
        m = min(c, d)
        a += m
        b += m
        c -= m
        d -= m
        if d == 0:
            return 2 * m + 1, [a, b, c, 0, 1, 0, 0, 0]      # -> S11 (rule10)
        else:
            return 2 * m + 1, [a, b, 0, d - 1, 0, 0, 1, 0]  # -> S17 (extra rule1)

    if s17:  # S17
        if b:
            return 1, [a, b - 1, c, d, 0, 0, 0, 1]          # -> S19 (rule3)
        else:
            return 1, [a, b, c, d, 0, 0, 0, 0]              # -> S0  (rule9)

    if s19:  # S19
        return 2 * a + 1, [0, b, c + a, d + 1, 1, 0, 0, 0]  # (4+7)*a then rule8 -> S11

    # S0
    return a + d + 1, [0, b + a, c + a + 1, 0, 1, 0, 0, 0]  # 12^a,13^d,14 -> S11


def run_fractran(prime_index=10001):
    target_exp = int(prime(prime_index))
    e = [1, 0, 0, 0, 0, 0, 0, 0]  # seed = 2
    cost = 0
    while True:
        if e[0] == target_exp and all(x == 0 for x in e[1:]):
            return cost
        dc, e = batch_step(e)
        cost += dc


if __name__ == "__main__":
    print(run_fractran(3))

为了后续方便，我们接下来用a,b,c,d分别替代e[0],e[1],e[2],e[3]。

首先，我们注意到 \(S13\) 与 \(S11\) 之间形成了一个确定的状态循环；并且 \(S13\) 的入口唯一，只能由 \(S11\) 转入。因此，一旦进入 \(S13\)，本轮可执行的二步宏次数就被固定为 b = min(c, d)。当 d 被减到 0 时，b 相当于对 d 的一次“备份”：系统回到 \(S11\) 后会利用这份备份把除数恢复出来，并继续重复该循环，不断消耗 c，直到 c = 0 为止。若最终 c == 0，则进入 \(S17\) 后，b 恰好变为余数意义下的 c % b。因此，假设起点是“刚进入 \(S11\)”，那么 \(S11\rightarrow S13\rightarrow S11\rightarrow S13\rightarrow \dots\rightarrow S17\) 这一整段过程可以被进一步简化为如下代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
# 起点：刚进入 S11（标记11在），此时 a,b,c,d 是当前寄存器值
# S11 -> S13
cost += 2 * b        # (5+6) 宏 b 次
d += b
b = 0
cost += 1            # 规则 11：11 -> 13，进入 S13


# 现在确实处于 “进入 S13 的时刻”（b==0 且 d>0）
while c >= d:
    c -= d
    a += d
    cost += 4 * d + 2      # 一次完整的 “减掉一个 d”：S13(2d步) ->S11(+1)->S11(2d步)->S13(+1)

# 现在 c 就是余数 r (0 <= r < d)
r = c
a += r
b = r                      # 进入 S17 时 b 记录余数 r
c = 0
d = d - r - 1              # 先做 r 次(1+2)，再做 1 次规则1：d := d - r - 1
cost += 2 * r + 1          # = 2r（r 次(1+2)） + 1（规则1进 S17）

进入到 \(S17\) 之后，若 b > 0，则会转入 \(S11\)，从而开启下一轮同构的循环；若 b = 0，则保持 a 的值不变，此时同时有 b = 0, c = 0，系统进入 \(S0\)，并且 d -= 1。

因此大致代码进一步变成：

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
while True:
    # 起点：刚进入 S11（标记11在），此时 a,b,c,d 是当前寄存器值

    # S11 -> S13
    cost += 2 * b          # (5+6) 宏 b 次
    d += b
    b = 0
    cost += 1              # 规则11：11 -> 13，进入 S13

    # S13：做 q 次“完整减法”（每次净效应 c-=d, a+=d，且 d 被恢复）
    while c >= d:
        c -= d
        a += d
        cost += 4 * d + 2  # = 2d(1+2) +1(10) +2d(5+6) +1(11)

    # 现在 c 就是余数 r (0 <= r < d)
    r = c
    a += r                 # 使 a 变成 n = qd+r
    b = r                  # 把余数记到 b
    c = 0
    d = d - r - 1          # r 次(1+2) 抽干 c，再 1 次规则1 进 S17
    cost += 2 * r + 1      # = 2r + 1

    # 此时：刚进入 S17（还没执行 S17 的那一步）
    if b > 0:              # r>0：执行规则3，进 S19
        cost += 1          # 规则3
        b -= 1

        # S19：执行 (4+7) 共 a 次，把 a 的值搬到 c；再执行规则8 回 S11
        cost += 2 * a + 1  # 2a(4+7) + 1(8)
        c += a
        a = 0
        d += 1             # 规则8：d += 1，同时回到 S11（标记11恢复）

        # 循环继续：回到 while True 顶部（起点仍是“刚进入 S11”）
    else:                  # r=0：执行规则9，进 S0，然后 S0 批处理回到 S11
        cost += 1          # 规则9：离开 S17 进入 S0
        break

接下来纳入 \(S0\) 状态，因此外层可以再嵌套一个循环，得到完整的代码：

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
49
50
51
a = 1
b = c = d = 0
cost = 0
N = 997
while True:
    if a == N:
        break
    cost += a + d + 1
    b += a
    c += a + 1
    a = 0
    d = 0
    while True:
        # 起点：刚进入 S11（标记11在），此时 a,b,c,d 是当前寄存器值

        # S11 -> S13
        cost += 2 * b          # (5+6) 宏 b 次
        d += b
        b = 0
        cost += 1              # 规则11：11 -> 13，进入 S13

        # S13：做 q 次“完整减法”（每次净效应 c-=d, a+=d，且 d 被恢复）
        while c >= d:
            c -= d
            a += d
            cost += 4 * d + 2  # = 2d(1+2) +1(10) +2d(5+6) +1(11)

        # 现在 c 就是余数 r (0 <= r < d)
        r = c
        a += r                 # 使 a 变成 n = qd+r
        b = r                  # 把余数记到 b
        c = 0
        d = d - r - 1          # r 次(1+2) 抽干 c，再 1 次规则1 进 S17
        cost += 2 * r + 1      # = 2r + 1

        # 此时：刚进入 S17（还没执行 S17 的那一步）
        if b > 0:              # r>0：执行规则3，进 S19
            cost += 1          # 规则3
            b -= 1

            # S19：执行 (4+7) 共 a 次，把 a 的值搬到 c；再执行规则8 回 S11
            cost += 2 * a + 1  # 2a(4+7) + 1(8)
            c += a
            a = 0
            d += 1             # 规则8：d += 1，同时回到 S11（标记11恢复）

        # 循环继续：回到 while True 顶部（起点仍是“刚进入 S11”）
        else:                  # r=0：执行规则9，进 S0，然后 S0 批处理回到 S11
            cost += 1          # 规则9：离开 S17 进入 S0
            break
print(cost)

现在开始精简这份代码。可见，最内层的 while 循环本质上是求个商值，c 就变成了余数，因此可以去掉一个循环，变成：

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
49
a = 1
b = c = d = 0
cost = 0
N = 997
while True:
    if a == N:
        break
    cost += a + d + 1
    b += a
    c += a + 1
    a = 0
    d = 0
    while True:
        # 起点：刚进入 S11（标记11在），此时 a,b,c,d 是当前寄存器值

        # S11 -> S13
        cost += 2 * b          # (5+6) 宏 b 次
        d += b
        b = 0
        cost += 1              # 规则11：11 -> 13，进入 S13

        q, r = divmod(c, d)
        cost += (4 * d + 2) * q
        a += d * q
        c = r
        # 现在 c 就是余数 r (0 <= r < d)
        a += r                 # 使 a 变成 n = qd+r
        b = r                  # 把余数记到 b
        c = 0
        d = d - r - 1          # r 次(1+2) 抽干 c，再 1 次规则1 进 S17
        cost += 2 * r + 1      # = 2r + 1

        # 此时：刚进入 S17（还没执行 S17 的那一步）
        if b > 0:              # r>0：执行规则3，进 S19
            cost += 1          # 规则3
            b -= 1

            # S19：执行 (4+7) 共 a 次，把 a 的值搬到 c；再执行规则8 回 S11
            cost += 2 * a + 1  # 2a(4+7) + 1(8)
            c += a
            a = 0
            d += 1             # 规则8：d += 1，同时回到 S11（标记11恢复）

        # 循环继续：回到 while True 顶部（起点仍是“刚进入 S11”）
        else:                  # r=0：执行规则9，进 S0，然后 S0 批处理回到 S11
            cost += 1          # 规则9：离开 S17 进入 S0
            break
print(cost)

随后 a 就变成原来的 c 的值，并且丢弃掉一下用不掉的初始化为 0 的逻辑，合并一些 cost 计算逻辑，因此进一步简化逻辑：

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
a = 1
b = c = d = 0
cost = 0
N = 997
while True:
    if a == N:
        break
    cost += a + d + 1
    b += a
    c += a + 1
    d = 0
    while True:
        # 起点：刚进入 S11（标记11在），此时 a,b,c,d 是当前寄存器值
        d += b
        q, r = divmod(c, d)
        cost += (4 * d + 2) * q + 2 * r + 2 + 2 * b + 1
        a = c
        b = r                  # 把余数记到 b
        c = 0
        d = d - r - 1          # r 次(1+2) 抽干 c，再 1 次规则1 进 S17
        # 此时：刚进入 S17（还没执行 S17 的那一步）
        if b > 0:              # r>0：执行规则3，进 S19
            b -= 1
            # S19：执行 (4+7) 共 a 次，把 a 的值搬到 c；再执行规则8 回 S11
            cost += 2 * a + 1  # 2a(4+7) + 1(8)
            c += a
            d += 1             # 规则8：d += 1，同时回到 S11（标记11恢复）

        # 循环继续：回到 while True 顶部（起点仍是“刚进入 S11”）
        else:                  # r=0：执行规则9，进 S0，然后 S0 批处理回到 S11
            break
print(cost)

由于进入\(S0\)一开始b, c必定为 0，因此改为使用=赋值号。而在内层while循环中，a总是作为c的备份来恢复，所以备份逻辑可以简化（甚至去掉，直接用原来的值a操作），到这里，c已经被彻底去掉：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
a = 1
b = c = d = 0
cost = 0
N = 997
while True:
    if a == N:
        break
    cost += a + d + 1
    b = a
    a += 1
    d = 0
    while True:
        d += b
        q, r = divmod(a, d)
        cost += (4 * d + 2) * q + 2 * r + 2 + 2 * b + 1
        b = r
        d = d - r - 1
        if b > 0:
            b -= 1
            cost += 2 * a + 1
            d += 1
        else:
            break
print(cost)

实际上，由于d每次在实际内层循环中都会减去 1，因此我们可以进一步简化代码以理清计算值的含义：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
a = 1
b = c = d0 = 0
cost = 0

N = 997
while True:
    if a == N:
        break
    a += 1
    cost += a + d0
    b = a - 1
    for d in range(b, 0, -1):
        q, r = divmod(a, d)
        cost += (4 * d + 2) * q + 2 * r + 2 + 2 * b + 1
        if r == 0:
            d0 = d - 1
            b = 0
            break
        b = r - 1
        cost += 2 * a + 1
print(cost)

也就是说，现在的程序可以理解成：对 \(n = 2,3,4,\dots\) 依次做一次“从大到小试除”的过程，并把 FRACTRAN 的步数累计到 cost。

令\(B(n)\)表示\(n\)的最大真因子。那么对于第一部分的 cost，它的总步数为：

\[S_1(N)=1+\sum_{n=3}^{N-1} n-1+B(n-1)=\dfrac{N(N-1)}{2}+\sum_{n=2}^{N-2} B(n)\]

对于第三部分，它的总步数为：

\[S_3(N)=\sum_{n=2}^N (2n+1)(n-1-B(n))\]

下面把上一条里的所有公式都改成用$...$或$$...$$包裹（内容不变）。

对固定的\(n\)，内层会从\(d=n-1\)一路试到\(d=B(n)\)才第一次命中并停止，因此第二条语句对该\(n\)的总贡献是

\[ C_2(n)=\sum_{d=B(n)}^{n-1}\left((4d+2)q_d+2r_d+2b_d+3\right), \]

其中\(q_d=\left\lfloor \dfrac{n}{d}\right\rfloor,\ r_d=n\bmod d\)，而\(b_d\)是执行到该轮试除\(d\)时代码里的\(b\)。

由代码更新规则可得：

  • 第一次试除\(d=n-1\)时\(b_{n-1}=n-1\)；
  • 若在\(d\)处失败（\(r_d>0\)），则下一轮试除\(d-1\)时用到\(b_{d-1}=r_d-1\)。

因此等价地，对\(B=B(n)\)，有\(b_{n-1}=n-1, b_d=r_{d+1}-1 (B\le d\le n-2).\)

因此\(\displaystyle{\sum_{d=B}^{n-1} b_d=(n-1)+\sum_{d=B}^{n-2}(r_{d+1}-1)= B+\sum_{d=B}^{n-1} r_d}\)，并且\(r_B=0\)（因为\(B\mid n\)）。

把\(C_2(n)\)拆项：

\[ C_2(n)= 4\sum_{d=B}^{n-1} d q_d +2\sum_{d=B}^{n-1} q_d +2\sum_{d=B}^{n-1} r_d +2\sum_{d=B}^{n-1} b_d +3(n-B). \]

用恒等式\(r_d=n-dq_d\)得到

\[ \sum_{d=B}^{n-1} r_d = n(n-B)-\sum_{d=B}^{n-1} d q_d. \]

再代入\(\displaystyle{\sum_{d=B}^{n-1} b_d= B+\sum_{d=B}^{n-1} r_d}\)，得

\[ 2\sum_{d=B}^{n-1} b_d =2B+2\sum_{d=B}^{n-1} r_d =2B+2n(n-B)-2\sum_{d=B}^{n-1} d q_d. \]

代回\(C_2(n)\)后，\(\sum d q_d\)项完全抵消，最终得到

\[ C_2(n)= 2\sum_{d=B(n)}^{n-1}\left\lfloor\frac{n}{d}\right\rfloor +4n\left(n-B(n)\right) +3n-B(n). \]

如果外层让\(n\)从\(2\)跑到\(N\)，那么第二条语句的总贡献就是

\[ \begin{aligned}S_2(n)&=\sum_{n=2}^{N} \left( 2\sum_{d=B(n)}^{n-1}\left\lfloor\frac{n}{d}\right\rfloor +4n(n-B(n)) +3n-B(n) \right) \\&=2\sum_{n=2}^{N}\sum_{d=B(n)}^{n-1}\left\lfloor\frac{n}{d}\right\rfloor +\left(\frac{N(N+1)(8N+13)}{6}-7\right) -\sum_{n=2}^{N}(4n+1)B(n). \end{aligned} \]

那么最终的结果为：

\[ \begin{aligned}S(n)&=S_1(n)+S_2(n)+S_3(n)\\ &= \frac{4N^3+9N^2+N-12}{2} + 2\sum_{n=2}^{N}\sum_{d=B(n)}^{n-1}\Bigl\lfloor\frac{n}{d}\Bigr\rfloor + \sum_{n=2}^{N-2}B(n) - \sum_{n=2}^{N}(6n+2)B(n)\\ &= \frac{4N^3+9N^2+N-12}{2} + 2\sum_{n=2}^{N}\sum_{d=B(n)}^{n-1}\Bigl\lfloor\frac{n}{d}\Bigr\rfloor -\sum_{n=2}^{N-2}(6n+1)B(n)- (6N-4)B(N-1)-(6N+2)B(N). \end{aligned} \]

注意，第二块可以使用数论分块即可高效求解。

代码

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
49
50
51
from sympy import prime

from tools import *


def sum_floor_range(n: int, L: int, R: int) -> int:
    """sum_{d=L..R} floor(n/d) via divisor blocking."""
    if L > R:
        return 0
    s = 0
    i = L
    while i <= R:
        q = n // i
        j = n // q
        if j > R:
            j = R
        s += q * (j - i + 1)
        i = j + 1
    return s


def compute_S(N: int) -> int:
    if N < 2:
        return 0

    spf = Factorizer(N).spf

    # B(n) = largest proper divisor = n / spf(n)
    B = [0 for _ in range(N+1)]
    for n in range(2, N + 1):
        B[n] = n // spf[n]

    # polynomial part: (4N^3 + 9N^2 + N - 12) / 2
    poly = (4 * N * N * N + 9 * N * N + N - 12) // 2

    # F = sum_{n=2..N} sum_{d=B(n)..n-1} floor(n/d)
    F = 0
    for n in range(2, N + 1):
        F += sum_floor_range(n, B[n], n - 1)

    # B-term: -sum_{n=2..N-1} (6n+1)B(n) -(6N+2)B(N)
    Bsum = 0
    for n in range(2, N):
        Bsum += (6 * n + 1) * B[n]

    return poly + 2 * F - Bsum - (6 * N + 2) * B[N]


if __name__ == "__main__":
    print(compute_S(prime(10001)))  # 1274952