Project Euler 543

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

Project Euler 543

题目

Prime-Sum Numbers

Define function \(P(n,k) = 1\) if \(n\) can be written as the sum of \(k\) prime numbers (with repetitions allowed), and \(P(n,k) = 0\) otherwise.

For example, \(P(10,2) = 1\) because \(10\) can be written as either \(3 + 7\) or \(5 + 5\), but \(P(11,2) = 0\) because no two primes can sum to \(11\).

Let \(S(n)\) be the sum of all \(P(i,k)\) over \(1 \le i,k \le n\).

For example, \(S(10) = 20, S(100) = 2402,\) and \(S(1000) = 248838\).

Let \(F(k)\) be the \(k\text{th}\) Fibonacci number (with \(F(0) = 0\) and \(F(1) = 1\)).

Find the sum of all \(S(F(k))\) over \(3 \le k \le 44\)

解决方案

令\(\displaystyle{T(n,k)=\sum_{i=1}^{n} P(i,k),}\)即：在 \(1\le i\le n\) 中，有多少个整数可以表示成 \(k\) 个素数之和。则\(\displaystyle{S(n)=\sum_{k=1}^{n} T(n,k).}\)所以问题转化为：求出所有 \(T(n,k)\) 的闭式表达，再把它们加起来。

本题只会用到如下三个结论：

  1. 偶数的二素数分解（强 Goldbach 猜想）：任意偶数 \(m>2\) 都能写成两个素数之和，因此对所有偶数 \(m\ge 4\)，\(P(m,2)=1.\)
  2. 奇数的三素数分解（弱 Goldbach 猜想）：任意奇数 \(m\ge 7\) 都能写成三个素数之和，因此对所有奇数 \(m\ge 7\)，\(P(m,3)=1.\)
  3. 允许重复 + 可以补 \(2\) 因为 \(2\) 是素数，所以当已知某个数能用 \(r\) 个素数表示时，就可以在表达式中额外加入若干个 \(2\)，把“素数个数”抬高。

这些事实将使得 \(k\ge 3\) 的部分完全变成一个简单的线性计数。

接下来我们分别对\(k=1,2,\ge 3\)三种不同的情况分别讨论\(T(n,k)\)。

\(T(n,1)\)

\(P(i,1)=1\) 当且仅当 \(i\) 本身是素数，因此

\[T(n,1)=\pi(n),\]

其中 \(\pi(n)\) 表示不超过 \(n\) 的素数个数。

\(T(n,2)\)

对于偶数部分，根据强 Goldbach，所有偶数 \(i\ge 4\) 都满足 \(P(i,2)=1\)。在 \(1\le i\le n\) 中，偶数 \(\ge 4\) 的个数是\(\left\lfloor\dfrac n2\right\rfloor-1.\)因此偶数贡献为\(\max\left(\left\lfloor\dfrac n2\right\rfloor-1,0\right).\)

对于偶数部分，若 \(i\) 是奇数且 \(P(i,2)=1\)，那么两个素数之和为奇数，必定有且仅有一个是 \(2\)，于是\(i=2+p,\)且\(p\)为素数。

因此奇数 \(i\le n\) 可写成两素数之和，当且仅当 \(i-2\) 是素数。这样的 \(i\) 的个数正好等于 \(n-2\) 以内的素数个数去掉素数 \(2\)：\(\max(\pi(n-2)-1,0).\)综上

\[ T(n,2)=\max\left(\left\lfloor\frac n2\right\rfloor-1,0\right)+\max(\pi(n-2)-1,0). \]

\(T(n,k),k\ge 3\)

关键观察：当 \(k\ge 3\) 时，几乎所有足够大的整数都能表示为 \(k\) 个素数之和，而且判定条件只剩下“是否至少为 \(2k\)”。解析来我们证明这个性质。

首先证明充分性：若 \(i\ge 2k\)，则 \(P(i,k)=1\)

  • 若 \(i\) 为偶数且 \(i\ge 2k\)，令\(i' = i-2(k-2).\)则 \(i'\) 仍为偶数，且\(i'\ge 4.\)由强 Goldbach，\(i'\) 是两个素数之和，因此 \(i\) 是 \(k\) 个素数之和（再补 \(k-2\) 个 \(2\)）。
  • 若 \(i\) 为奇数且 \(i\ge 2k\)，由于奇数至少为 \(2k+1\)，令\(i' = i-2(k-3).\)则 \(i'\) 为奇数且\(i'\ge 7.\)由弱 Goldbach，\(i'\) 是三个素数之和，因此 \(i\) 是 \(k\) 个素数之和（再补 \(k-3\) 个 \(2\)）。

因此当 \(k\ge 3\) 时，只要 \(i\ge 2k\) 就必定可表示。

接下来证明必要性：若 \(i<2k\)，则 \(P(i,k)=0.\)因为每个素数至少为 \(2\)，\(k\) 个素数之和至少为 \(2k\)，所以 \(i<2k\) 不可能。于是对 \(k\ge 3\) 有严格等价：\(P(i,k)=1\iff i\ge 2k.\)

因此有

\[ T(n,k)=\sum_{i=1}^{n} P(i,k)=\#\{i\le n:\ i\ge 2k\}=\max(n-2k+1,0). \]

由定义可分得\(\displaystyle{S(n)=T(n,1)+T(n,2)+\sum_{k=3}^{n}T(n,k).}\)其中

\[\begin{aligned} T(n,1)&=\pi(n),\\ T(n,2)&=\max\left(\left\lfloor\frac n2\right\rfloor-1,0\right)+\max(\pi(n-2)-1,0),\\ \sum_{k=3}^{n}T(n,k)&=\sum_{k=3}^{\lfloor n/2\rfloor} (n-2k+1). \end{aligned} \]

设\(m=\left\lfloor\dfrac n2\right\rfloor.\)那么可以写：

\[ \sum_{k=3}^{n}T(n,k)= \begin{cases} 0,& m<3,\\ (m-2)\bigl(m-2+(n\bmod 2)\bigr),& m\ge 3. \end{cases} \]

因此 \(S(n)\) 完全由 \(\pi(n)\)、\(\pi(n-2)\) 与一个二次多项式组成。为了高效计算\(\pi\)，题目642提供了min-25筛法的实现。

代码

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
from math import isqrt
from functools import lru_cache
from tools import Factorizer, int_sqrt

N = 45
SIEVE_MAX = 1_000_000
factorizer = Factorizer(SIEVE_MAX)
primes = factorizer.primes
spf = factorizer.spf
pi_small = [0] * (SIEVE_MAX + 1)
c = 0
for i in range(2, SIEVE_MAX + 1):
    if spf[i] == i:
        c += 1
    pi_small[i] = c


def iroot3(n: int) -> int:
    x = int(round(n ** (1 / 3)))
    while (x + 1) ** 3 <= n:
        x += 1
    while x ** 3 > n:
        x -= 1
    return x


def iroot4(n: int) -> int:
    return int_sqrt(int_sqrt(n))


@lru_cache(maxsize=None)
def phi(x: int, a: int) -> int:
    if a == 0:
        return x
    if a == 1:
        return x - x // 2
    return phi(x, a - 1) - phi(x // primes[a - 1], a - 1)


@lru_cache(maxsize=None)
def prime_pi(x: int) -> int:
    if x < SIEVE_MAX:
        return pi_small[x]
    a = prime_pi(iroot4(x))
    b = prime_pi(isqrt(x))
    c = prime_pi(iroot3(x))
    res = phi(x, a) + ((b + a - 2) * (b - a + 1)) // 2
    for i in range(a, b):
        p = primes[i]
        w = x // p
        res -= prime_pi(w)
        if i < c:
            lim = prime_pi(isqrt(w))
            for j in range(i, lim):
                res -= prime_pi(w // primes[j]) - j
    return res


def S(n: int) -> int:
    if n == 2:
        return 1
    if n == 3:
        return 2
    if n == 5:
        return 5
    t1 = prime_pi(n)
    t2 = 0
    if n >= 4:
        t2 += n // 2 - 1
    if n > 2:
        t2 += prime_pi(n - 2) - 1
    m = n // 2
    t3 = 0
    if m >= 3:
        t3 = (m - 2) * (m - 2 + (n & 1))
    return t1 + t2 + t3


def solve(N) -> int:
    F = [0] * N
    F[0] = 0
    F[1] = 1
    for i in range(2, N):
        F[i] = F[i - 1] + F[i - 2]
    ans = 0
    for k in range(3, N):
        ans += S(F[k])
    return ans


if __name__ == "__main__":
    print(solve(N))