Project Euler 341

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

Project Euler 341

题目

Golomb’s self-describing sequence

The Golomb’s self-describing sequence \((G(n))\) is the only nondecreasing sequence of natural numbers such that \(n\) appears exactly \(G(n)\) times in the sequence. The values of \(G(n)\) for the first few \(n\) are

\[\begin{matrix} n & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10 & 11 & 12 & 13 & 14 & 15 & \ldots \\ G(n) & 1 & 2 & 2 & 3 & 3 & 4 & 4 & 4 & 5 & 5 & 5 & 6 & 6 & 6 & 6 & \ldots \end{matrix}\]

You are given that \(G(10^3) = 86\), \(G(10^6) = 6137\).

You are also given that \(\sum G(n^3) = 153506976\) for \(1 \le n \lt 10^3\).

Find \(\sum G(n^3)\) for \(1 \le n \lt 10^6\).

解决方案

定义前缀和\(\displaystyle{F(m)=\sum_{i=1}^{m} G(i), F(0)=0.}\)

由于“整数 \(m\) 出现 \(G(m)\) 次”且序列非递减，所有值不超过 \(m\) 的项必然构成序列的最前一段，因此：

  • 值为 \(m\) 的项，正好占据下标区间\((F(m-1)+1,F(m-1)+2,\dots,F(m)).\)
  • 等价地，对任意下标 \(n\)，\(G(n)=m \iff F(m-1)<n\le F(m).\)

这给出一个非常重要的“反演”描述：\(G(n)=\min\{m:F(m)\ge n\}.\)

接着定义一个加权前缀和\(\displaystyle{S(m)=\sum_{i=1}^{m} i\cdot G(i), S(0)=0.}\)

注意到序列前 \(F(m)\) 项的取值恰好是 \(1,2,\dots,m\)，其中 \(i\) 出现 \(G(i)\) 次，因此前 \(F(m)\) 项的数值和为\(\displaystyle{\sum_{k=1}^{F(m)} G(k)=\sum_{i=1}^{m} i\cdot G(i)=S(m).}\)

而左侧正是 \(F\) 在点 \(F(m)\) 的取值，于是得到关键恒等式：\(F(F(m))=S(m).\)

固定一个 \(m\)。在下标区间\((F(m-1),F(m)]\)内，序列恒等于 \(m\)，即 \(G(t)=m\)。因此 \(F(t)\) 在这段上每次增加 \(m\)，并且起点为\(F(F(m-1))=S(m-1).\)

于是对任意 \(q=1,2,\dots,G(m),F(F(m-1)+q)=S(m-1)+m\cdot q.\)

也就是说，\(F\) 在这一整段上是严格等差的。现在要计算 \(G(N)\)，根据 \(G(N)=\min\{t:F(t)\ge N\}\)，如果我们找到唯一的 \(m\) 使得\(S(m-1)<N\le S(m),\)那么 \(N\) 必落在这段等差范围内。令 \(s=S(m-1)\)，则最小的 \(q\) 满足\(s+m\cdot q\ge N\)为\(q=\left\lceil\dfrac{N-s}{m}\right\rceil.\)

对应的下标 \(t\)（也就是 \(G(N)\)）为\(G(N)=F(m-1)+q.\)

这就是可以直接用于扫描的计算公式：若\(S(m-1)<N\le S(m)\)，那么\(G(N)=F(m-1)+\left\lceil\dfrac{N-S(m-1)}{m}\right\rceil.\)

Golomb 序列有经典递推（可由“每个整数出现次数由自身给出”与上面的反演结构推导得到）： \[G(1)=1,G(n)=1+G\bigl(n-G(G(n-1))\bigr), (n>1).\]

该递推只依赖更小下标的值，因此可以线性生成 \(G(1\dots M)\)。

最终，我们仅需要维护：

  • \(G(m)\)：当前递推得到；
  • \(F(m) = F(m-1) + G(m)\)；
  • \(S(m) = S(m-1) + mG(m)\)。

每当某个 \(n\) 落入 \((S(m-1),S(m)]\)，用\(q=\left\lceil\dfrac{N-S(m-1)}{m}\right\rceil\)并累加\(G(N)=F(m-1)+q.\)

代码

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
from array import array

def sum_golomb_cubes(limit_n: int = 10**6) -> int:
    """
    计算 sum_{1<=n<limit_n} G(n^3)。

    关键：
      F(m) = sum_{i=1..m} G(i)
      S(m) = sum_{i=1..m} i*G(i) = F(F(m))

    若 S(m-1) < N <= S(m)，则
      G(N) = F(m-1) + ceil((N - S(m-1)) / m)
    """

    # 1-indexed: g[n] = G(n)
    g = array('I', [0, 1])  # g[1] = 1

    # F(m-1) 与 S(m-1)
    F_prev = 0              # F(m-1) = sum_{i<=m-1} G(i)
    S_prev = 0              # S(m-1) = sum_{i<=m-1} i*G(i)

    ans = 0

    # 立方点 N = n^3, n = 1..limit_n-1
    n = 1
    target = 1  # 1^3

    m = 1
    while n < limit_n:
        # 生成 G(m)
        if m == 1:
            Gm = 1
        else:
            # Golomb 递推：G(m) = 1 + G(m - G(G(m-1)))
            Gm = 1 + g[m - g[g[m - 1]]]
            g.append(Gm)

        # 当前段右端：S(m) = S(m-1) + m*G(m)
        S_curr = S_prev + m * Gm

        # 处理所有落在 (S_prev, S_curr] 的立方点
        while n < limit_n and target <= S_curr:
            # q = ceil((target - S_prev) / m)
            q = (target - S_prev + m - 1) // m
            # G(target) = F(m-1) + q
            ans += F_prev + q

            n += 1
            target = n * n * n  # 下一个立方点

        # 推进到下一个 m：更新 F(m), S(m)
        F_prev += Gm
        S_prev = S_curr
        m += 1

    return ans


if __name__ == "__main__":
    # 题目给出的校验：
    # G(10^3) = 86, G(10^6) = 6137
    # sum_{1<=n<10^3} G(n^3) = 153506976
    #
    # 直接跑主答案：
    print(sum_golomb_cubes(10**6))  # 56098610614277014