Project Euler 759

发表于 2022-08-21 更新于 2023-07-03 分类于 Project Euler 阅读次数： 111 本文字数： 2.4k 阅读时长 ≈ 2 分钟

Project Euler 759

题目

A squared recurrence relation

The function $f$ is defined for all positive integers as follows:

$\begin{aligned} f (1) & = 1 \\ f (2 n) & = 2 f (n) \\ f (2 n + 1) & = 2 n + 1 + 2 f (n) + \frac{1}{n} f (n) \end{aligned}$

It can be proven that $f (n)$ is integer for all values of $n$ .

The function $S (n)$ is defined as $S (n) = \sum_{i = 1}^{n} f (i)^{2}$ .

For example, $S (10) = 1530$ and $S (10^{2}) = 4798445$ .

Find $S (10^{16})$ . Give your answer modulo $1 000 000 007$ .

解决方案

暴力枚举出 $f$ 的前几项，在 OEIS 中查询得到结果为 A245788。

令 $b (n)$ 表示 $n$ 的二进制表示中有多少个 $1$ ，那么按照题意， $f (n) = n \cdot b (n)$ 。

将 $S$ 进行改写：

$S (n) = \sum_{i = 1}^{⌊ \log_{2} (n + 1) ⌋} (i^{2} \cdot \sum_{\begin{matrix} m \leq n, b (m) = i \end{matrix}} m^{2})$

也就是说，问题转化成枚举 $i$ ，求 $n$ 以内的所有数中，满足 $b (m) = i$ 的所有 $m^{2}$ 之和。

将 $n$ 写成一个长度为 $l$ 的二进制数字符串 $d_{1} d_{2} \dots d_{l}$ ，也就是有 $l = ⌊ \log_{2} (n + 1) ⌋$ 。

令中间状态 $c_{1} (i, j) (1 \leq i \leq l, 0 \leq j \leq 9 i)$ 分别表示如下含义：有多少个 $i$ 位有前导 0 的数，其二进制有 $j$ 个 $1$ ，并等于由字符串 $d_{1} d_{2} \dots d_{i}$ 表示的数。令状态 $s_{1} (i, j)$ 表示 $c_{1} (i, j)$ 中的所有数之和。令状态 $t_{1} (i, j)$ 表示 $c_{1} (i, j)$ 中的所有数的平方和。

不难发现，如果 $j = \sum_{k = 1}^{i} d_{k}$ ，那么 $c_{1} (i, j) = 1$ ， $s_{1} (i, j)$ 的值则为 $d_{1} d_{2} \dots d_{i}$ 这一个字符串的二进制数，并且 $t_{1} (i, j) = s_{1} (i, j)^{2}$ ；否则 $c_{1} (i, j) = s_{1} (i, j) = t_{1} (i, j) = 0$ 。

状态 $c_{1}, s_{1}, t_{1}$ 的意义在于规定了答案严格在界限上的情况。而接下来定义的状态 $c_{0}, s_{0}, t_{0}$ 则是严格不在界限上的情况。

令状态 $c_{0} (i, j)$ 表示有多少个 $i$ 位有前导 0 的数，其数位和为 $j$ ，并严格小于由字符串 $d_{1} d_{2} \dots d_{i}$ 表示的数。令状态 $s_{0} (i, j)$ 表示 $c_{0} (i, j)$ 中的所有数之和， $t_{0} (i, j)$ 表示 $c_{0} (i, j)$ 。中的所有数的平方和。

那么可以先写出 $c_{0}$ 的状态转移方程：

$c_{0} (i, j) = {\begin{aligned} 1 & if j = 0 \\ 0 & else if i = 1 \\ c_{0} (i - 1, j - 1) + c_{0} (i - 1, j) + [d_{i} = 1] \cdot c_{1} (i - 1, j) & else \end{aligned}$

其中， $[]$ 表示示性函数，如果 $[]$ 内的表达式为真，那么值为 $1$ ，否则为 $0$ 。方程的前两项，如果之前的状态已经是严格小于 $d_{1} d_{2} \dots d_{i - 1}$ 的数，那么接下来无论在后面拼接一个 $0$ 还是一个 $1$ ，都是严格小于的；对于第三项而言，如果之前取的值都是严格等于的，那么接下来如果 $d_{i} = 1$ ，那么可以拼接一个 $0$ ，从而达到严格小于的情况。

根据上面的思想，不难写出 $s_{0} (i, j)$ 和 $t_{0} (i, j)$ ：

$s_{0} (i, j) = {\begin{aligned} 0 & if i = 1 \lor j = 0 \\ 2 \cdot s_{0} (i - 1, j - 1) + c_{0} (i - 1, j - 1) + 2 \cdot s_{0} (i - 1, j) + \\ [d_{i} = 1] \cdot 2 \cdot s_{1} (i - 1, j) & else \end{aligned}$

$t_{0} (i, j) = {\begin{aligned} 0 & if i = 1 \lor j = 0 \\ 4 \cdot t_{0} (i - 1, j - 1) + 4 \cdot s_{0} (i - 1, j - 1) + c_{0} (i - 1, j - 1) + 4 \cdot t_{0} (i - 1, j) + \\ [d_{i} = 1] \cdot 2 \cdot s_{1} (i - 1, j) & else \end{aligned}$

需要注意的是，对于 $t_{0}$ 的状态转移方程最后一行的前三项，系数是 $4, 4, 1$ ，来自于平方式 $(2 n + 1)^{2} = 4 n^{2} + 4 n + 1$ 。

计算完成后，得到最终答案为：

$S (n) = \sum_{i = 1}^{l} i^{2} \cdot (t_{0} (l, i) + t_{1} (l, i))$

代码

N = 10 ** 16
mod = 1000000007
d = list(map(int, list(bin(N)[2:])))
l = len(d)
c = [[[0, 0] for _ in range(l + 1)] for _ in range(l)]
s = [[[0, 0] for _ in range(l + 1)] for _ in range(l)]
t = [[[0, 0] for _ in range(l + 1)] for _ in range(l)]
c[0][1][1] = s[0][1][1] = t[0][1][1] = 1
c[0][0][0] = 1
cnt = 1
for i in range(1, l):
    cnt += d[i]
    c[i][cnt][1] = 1
    s[i][cnt][1] = s[i - 1][cnt - d[i]][1] * 2 + d[i]
    # 因为d[i]=0,1,d[i] ** 2 = d[i]
    t[i][cnt][1] = t[i - 1][cnt - d[i]][1] * 4 + s[i - 1][cnt - d[i]][1] * d[i] * 4 + d[i]
    c[i][0][0] = 1
    s[i][0][0] = t[i][0][0] = 0
    for j in range(1, l + 1):
        c[i][j][0] = c[i - 1][j - 1][0] + c[i - 1][j][0]
        s[i][j][0] = s[i - 1][j - 1][0] * 2 + c[i - 1][j - 1][0] + s[i - 1][j][0] * 2
        t[i][j][0] = t[i - 1][j - 1][0] * 4 + s[i - 1][j - 1][0] * 4 + c[i - 1][j - 1][0] + t[i - 1][j][0] * 4
        if d[i] == 1:
            c[i][j][0] += c[i - 1][j][1]
            s[i][j][0] += s[i - 1][j][1] * 2
            t[i][j][0] += t[i - 1][j][1] * 4

ans = 0
for i in range(l + 1):
    ans += sum(t[-1][i]) * i * i
ans %= mod
print(ans)