Project Euler 368

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

Project Euler 368

题目

A Kempner-like series

The harmonic series \(1 + \dfrac{1}{2} + \dfrac{1}{3} + \dfrac{1}{4} + \dots\) is well known to be divergent.

If we however omit from this series every term where the denominator has a \(9\) in it, the series remarkably enough converges to approximately \(22.9206766193\).

This modified harmonic series is called the Kempner series.

Let us now consider another modified harmonic series by omitting from the harmonic series every term where the denominator has \(3\) or more equal consecutive digits.

One can verify that out of the first \(1200\) terms of the harmonic series, only \(20\) terms will be omitted.

These \(20\) omitted terms are: \[\dfrac{1}{111}, \dfrac{1}{222}, \dfrac{1}{333}, \dfrac{1}{444}, \dfrac{1}{555}, \dfrac{1}{666}, \dfrac{1}{777}, \dfrac{1}{888}, \dfrac{1}{999}, \dfrac{1}{1000}, \dfrac{1}{1110},\] \[\dfrac{1}{1111}, \dfrac{1}{1112}, \dfrac{1}{1113}, \dfrac{1}{1114}, \dfrac{1}{1115}, \dfrac{1}{1116}, \dfrac{1}{1117}, \dfrac{1}{1118}, \dfrac{1}{1119}\]

This series converges as well.

Find the value the series converges to.

Give your answer rounded to \(10\) digits behind the decimal point.

解决方案

限制条件只与“连续相同数字的长度”有关，因此用末尾两位是否相同来刻画状态最方便。

对 \(n\ge 2\)、\(d\in\{0,\dots,9\}\)，定义两类 合法的 \(n\) 位正整数集合（首位非零，且不含三连相同）：

  • \(S_1(n,d)\)：末位为 \(d\)，倒数第二位 不是 \(d\)（即末尾同数字连段长度为 \(1\)）。
  • \(S_2(n,d)\)：末两位都是 \(d\)（即末尾同数字连段长度为 \(2\)）。

再定义幂和（\(j\ge 1\)）：\(\displaystyle{f_1(n,d,j)=\sum_{x\in S_1(n,d)}\frac1{x^j},f_2(n,d,j)=\sum_{x\in S_2(n,d)}\frac1{x^j}.}\)

那么目标就是把所有位数的 \(j=1\) 汇总：

\[ S=\sum_{n\ge 1}\sum_{d=0}^9\bigl(f_1(n,d,1)+f_2(n,d,1)\bigr), \]

其中 \(n=1\) 时只有一位数 \(1,\dots,9\)。此时 \(S_2\) 无定义，可约定 \(S_2(1,d)=\varnothing\)（等价于 \(f_2(1,d,j)=0\)），因此一位数部分直接暴力计入即可。

把一个 \((n-1)\) 位数 \(y\) 追加一位 \(d\) 得到\(x=10y+d.\)我们可以基于动态规划的思想来计算。

  • 因此，若要得到 \(S_2(n,d)\)：末尾要变成 “\(dd\)”，因此 \((n-1)\) 位数必须以 单个 \(d\) 结尾，即来自 \(S_1(n-1,d)\)。有\(S_2(n,d)=\{10y+d:y\in S_1(n-1,d)\}.\)
  • 若要得到 \(S_1(n,d)\)：末尾连段长度为 \(1\)，说明 \((n-1)\) 位数的末位必须是 \(d'\ne d\)，它可以处于 \(S_1\) 或 \(S_2\)，有\(\displaystyle{S_1(n,d)=\bigcup_{d'\ne d}\{10y+d:y\in S_1(n-1,d')\cup S_2(n-1,d')\}.}\)

到这里为止，集合递推是完全精确的；难点是如何把它变成对 \(f_1,f_2\) 的快速递推，而不枚举所有整数。

接下来使用论文中的按位追加、用幂和递推思想求解。

对任意 \(j\ge 1\)，有\(\displaystyle{\frac{1}{(10y+d)^j}=\frac1{10^j y^j}\left(1+\frac{d}{10y}\right)^{-j}.}\)利用负二项级数：\(\displaystyle{(1+u)^{-j}=\sum_{t\ge 0}(-1)^t\binom{j+t-1}{t}u^t,}\)代入 \(u=\dfrac{d}{10y}\) 得到

\[ \frac1{(10y+d)^j} =\sum_{t\ge 0}(-1)^t\binom{j+t-1}{t}\frac{d^t}{10^{j+t}}\cdot \frac1{y^{j+t}}. \]

这正是论文中用于“按位追加、用幂和递推”的核心展开方式；并且由于这里是绝对收敛的幂级数，允许交换求和次序从而把“对 \(y\) 求和”移入每一项系数中。

我们针对项\(\dfrac{1}{y^{j+t}}\)，定义系数\(c(d;j,t)=(-1)^t\dbinom{j+t-1}{t}\dfrac{d^t}{10^{j+t}}.\)

由 \(S_2(n,d)=\{10y+d:y\in S_1(n-1,d)\}\)，得到

\[ f_2(n,d,j)=\sum_{y\in S_1(n-1,d)}\frac1{(10y+d)^j}=\sum_{t\ge 0}c(d;j,t)f_1(n-1,d,j+t). \]

由 \(S_1(n,d)\) 来自所有末位 \(d'\ne d\) 的 \(S_1\cup S_2\)：

\[ f_1(n,d,j)=\sum_{d'\ne d}\sum_{y\in S_1(n-1,d')\cup S_2(n-1,d')}\frac1{(10y+d)^j} =\sum_{t\ge 0}c(d;j,t)\sum_{d'\ne d}\bigl(f_1(n-1,d',j+t)+f_2(n-1,d',j+t)\bigr). \]

为实现上再做一步合并，记\(\displaystyle{g(n-1,d,k)=f_1(n-1,d,k)+f_2(n-1,d,k),G(n-1,k)=\sum_{d=0}^9 g(n-1,d,k),}\)则\(\displaystyle{\sum_{d'\ne d}g(n-1,d',k)=G(n-1,k)-g(n-1,d,k),}\)从而每个 \(d\) 的“排除自身”只需一次减法即可完成。

把所有合法正整数按十进制位数分块。令 \(T_n\) 表示合法的 \(n\) 位正整数个数，\(P_n\) 表示这些整数的集合（首位非零且不含三连）。对每个 \(n\)，所有 \(n\) 位数都满足 \(x\ge 10^{n-1}\)，因此有粗界\(\displaystyle{\sum_{x\in P_n}\frac1x\le\frac{T_n}{10^{n-1}}.}\)

关键是控制 \(T_n\) 的增长率：禁止三连是局部约束，可以用有限状态自动机计数，因此 \(T_n\) 满足固定维度的线性递推，从而存在常数 \(\lambda<10\) 使得\(T_n=O(\lambda^n).\)

把这两点合在一起就得到\(\displaystyle{\sum_{x\in P_n}\frac1x=O\left(\left(\frac{\lambda}{10}\right)^n\right),}\)也就是说按位数分块的贡献呈几何级数衰减。于是对任意 \(K\)，尾项满足

\[ \sum_{n>K}\sum_{x\in P_n}\frac1x=O\left(\left(\frac{\lambda}{10}\right)^K\right), \]

因此，选取足够大的 \(K\) 就能把尾项压到远小于最终需要的 \(10^{-10}\) 精度之外。

为了给出一个具体的 \(\lambda\)，可以用“末尾连续相同数字的长度”作状态。对 \(d\in{0,\dots,9}\) 定义：

  • \(U_n(d)\)：长度为 \(n\) 的合法串，末位为 \(d\) 且末尾连段长度为 \(1\) 的数量；
  • \(V_n(d)\)：长度为 \(n\) 的合法串，末位为 \(d\) 且末尾连段长度为 \(2\) 的数量。

（这里允许前导零；只影响常数项，不影响指数增长率。）

追加一位时有递推\(\displaystyle{V_{n+1}(d)=U_n(d),U_{n+1}(d)=\sum_{d'\ne d}(U_n(d')+V_n(d')).}\)把它们对 \(d\) 求和，令\(\displaystyle{U_n=\sum_{d=0}^9 U_n(d), V_n=\sum_{d=0}^9 V_n(d),}\)则得到 \(2\times2\) 的汇总递推

\[ \begin{pmatrix} U_{n+1}\\V_{n+1} \end{pmatrix}= \begin{pmatrix} 9 & 9\\ 1 & 0 \end{pmatrix} \begin{pmatrix} U_{n}\\V_{n} \end{pmatrix}. \]

因此增长率由该矩阵的最大特征值给出：\(\lambda=\dfrac{9+\sqrt{117}}2\approx 9.908<10,\)从而位数尾项确实按 \(\left(\dfrac{\lambda}{10}\right)^n\) 级别快速衰减。

递推 \(f_1,f_2\) 时要用到展开\(\displaystyle{\frac1{(10y+d)^j}=\sum_{t\ge 0}(-1)^t\binom{j+t-1}{t}\frac{d^t}{10^{j+t}}\cdot \frac1{y^{j+t}}.}\)这里的高阶项（大 \(t\)）之所以可以截断，关键在于参与递推的 \(y\) 足够大，从而使展开参数的绝对值很小。

当我们从 \((n-1)\) 位数 \(y\) 生成 \(n\) 位数 \(10y+d\) 时，总有 \(y\ge 10^{n-2}\)，因此\(\left|\dfrac{d}{10y}\right|\le \dfrac9{10^{n-1}}.\)

对 \(n\ge 4\)，右侧已经不超过 \(0.009\)。把展开写成\(\displaystyle{\frac1{(10y+d)^j}\frac1{10^j y^j}\left(1+\frac{d}{10y}\right)^{-j},}\)

可以看出这就是在 \(u=\dfrac{d}{10y}\) 处对 \((1+u)^{-j}\) 做幂级数展开；而当 \(|u|\) 极小时，级数的第 \(t\) 项会按 \(|u|^t\) 的速度快速衰减。因此只保留到某个最大阶 \(t\le J-j\)（等价于只维护到最高幂次 \(j+t\le J\)）时，被舍去的尾项之和将非常小，并且随着 \(n\) 增大衰减得更快。配合位数方向的截断，这保证了在固定的 \((K,J)\) 下即可达到所需的数值精度。

由此，在正式计算时，

  • 直接把 \(1\) 到 \(99\) 的调和和加入（两位以内不可能出现三连）。
  • 直接枚举 \(100\) 到 \(999\)，过滤掉 \(111,222,\dots,999\)，把它们的 \(\dfrac1n\) 加入总和；同时用这些三位数精确初始化 \(f_1(3,d,j), f_2(3,d,j)\)（对 \(1\le j\le J\)）。
  • 对 \(n=4,5,\dots,K\) 使用上面的递推计算新的 \(f_1,f_2\)，每步把 \(\sum_d\bigl(f_1(n,d,1)+f_2(n,d,1)\bigr)\) 加到总和。

代码

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
98
99
100
101
102
from tools import comb

def has_triple_repeat(s: str) -> bool:
    # contains some ddd
    for i in range(len(s) - 2):
        if s[i] == s[i + 1] == s[i + 2]:
            return True
    return False


def solve(J: int = 20, K: int = 4000) -> float:
    """
    Compute the convergent value of:
        sum_{n>=1, n has no 3 equal consecutive decimal digits} 1/n

    Parameters
    ----------
    J : int
        Max power tracked in f1/f2 (Taylor / negative-binomial truncation).
    K : int
        Max number of digits iterated (tail is negligible).

    Returns
    -------
    float
        Approximate limit value.
    """
    # ---- initial exact contribution for 1..99 ----
    total = 0.0
    for n in range(1, 100):
        total += 1.0 / n

    # ---- add all valid 3-digit numbers, and also build exact f1,f2 for 3 digits ----
    # f1[d][j] = sum_{x in S1(3,d)} x^{-j}, f2 similarly
    f1 = [[0.0] * (J + 1) for _ in range(10)]
    f2 = [[0.0] * (J + 1) for _ in range(10)]

    for n in range(100, 1000):
        if has_triple_repeat(str(n)):
            continue
        total += 1.0 / n
        d = n % 10
        d2 = (n // 10) % 10
        target = f2 if (d == d2) else f1  # S2 if last two digits equal, else S1
        inv = 1.0 / n
        v = inv
        for j in range(1, J + 1):
            target[d][j] += v
            v *= inv

    # ---- precompute coefficients c(d; j, t) = (-1)^t * C(j+t-1, t) * d^t / 10^(j+t) ----
    coeff = [[[0.0] * (J + 1) for _ in range(J + 1)] for __ in range(10)]
    for d in range(10):
        dpow = [1.0] * (J + 1)
        for t in range(1, J + 1):
            dpow[t] = dpow[t - 1] * d
        for j in range(1, J + 1):
            for t in range(0, J - j + 1):
                c = comb(j + t - 1, t)
                coeff[d][j][t] = ((-1.0) ** t) * c * dpow[t] / (10.0 ** (j + t))

    # ---- iterate digits from 4 to K ----
    for nd in range(4, K + 1):
        # g[d][j] = f1[d][j] + f2[d][j]
        g = [[f1[d][j] + f2[d][j] for j in range(J + 1)] for d in range(10)]
        total_g = [0.0] * (J + 1)
        for j in range(1, J + 1):
            s = 0.0
            for d in range(10):
                s += g[d][j]
            total_g[j] = s

        new1 = [[0.0] * (J + 1) for _ in range(10)]
        new2 = [[0.0] * (J + 1) for _ in range(10)]

        for d in range(10):
            cd = coeff[d]
            gd = g[d]
            f1d_prev = f1[d]
            for j in range(1, J + 1):
                s1 = 0.0
                s2 = 0.0
                max_t = J - j
                for t in range(0, max_t + 1):
                    co = cd[j][t]
                    idx = j + t
                    # f2: comes only from S1(n-1, d)
                    s2 += co * f1d_prev[idx]
                    # f1: comes from all last digits != d
                    s1 += co * (total_g[idx] - gd[idx])
                new1[d][j] = s1
                new2[d][j] = s2

        f1, f2 = new1, new2
        total += sum(f1[d][1] + f2[d][1] for d in range(10))

    return total


if __name__ == "__main__":
    ans = solve(J=20, K=4000)
    print(f"{ans:.10f}")