Project Euler 845

题目

Prime Digit Sum

Let \(D(n)\) be the \(n\)-th positive integer that has the sum of its digits a prime.

For example, \(D(61) = 157\) and \(D(10^8) = 403539364\).

Find \(D(10^{16})\).

解决方案

令\(N=10^{16}\)。设 \(S(x)\) 为整数 \(x\) 的十进制数位和。一个数是否满足条件，只取决于它的数位和是否为质数。与其直接枚举到第 \(N\) 个，不如：先按位数统计，即有多少个 \(L\) 位正整数满足数位和为质数；找到 \(D(N)\) 属于哪一位数 \(L\)；再从高位到低位逐位构造第 \(n\) 个最小解，可见这就是典型的DP 计数 + 字典序恢复

我们定义\(F(k, s)\)表示长度为 \(k\) 的十进制串（允许前导 \(0\)）中，有多少个串 \(x\) 满足\(s + S(x)\) 是质数。这里的 \(s\) 表示已经选定的前缀数位和。

当 \(k = 0\) 时，后缀为空串，\(S(x)=0\)，因此：\(F(0, s) = 1\) 当且仅当 \(s\) 是质数；否则 \(F(0, s) = 0\)。

当 \(k \ge 1\)时，设后缀第一位为 \(d(0\le d\le 9)\)，剩余 \(k-1\) 位构成串 \(y\)，则 \(S(x) = d + S(y)\)，条件变为 \(s + d + S(y)\) 为质数。

于是对固定 \(d\)，可行后缀数就是 \(F(k-1, s+d)\)，对所有 \(d\) 求和得到：

\[ F(k,s)=\sum_{d=0}^{9} F(k-1, s+d). \]

这就是题目所需的基本递推。

对 \(L\) 位正整数（最高位不能为 \(0\)），最高位 \(a\) 只能取 \(1\le a\le 9\)，剩余 \(L-1\) 位允许前导 \(0\)。

选定最高位 \(a\) 后，前缀数位和是 \(a\)，剩余可行个数为 \(F(L-1, a)\)。因此得到\(L\)位满足条件的数的数量为\(c(L)\)：

\[ c(L)=\sum_{a=1}^{9} F(L-1,a). \]

从 \(L=1,2,3,\dots\) 依次累加，找到最小的 \(L\) 使得累计数 \(\ge n\)，即可确定 \(D(n)\) 的位数，同时把 \(n\) 减去更短位数的数量，变成在固定长度 \(L\) 内找第 \(n\) 个。

接下来我们尝试字典序逐位恢复第 \(n\) 个最小解。已经确定长度为 \(L\)，并把 \(n\) 调整为在 \(L\) 位解集合中的排名（从 \(1\) 开始）。从高位到低位构造：

• 当前还剩 \(k\) 位未填，当前前缀数位和为 \(s\)。
• 尝试当前位取 \(d\)（首位 \(d=1\le d\le 9\)，其余位 \(d=0\le d\le 9\)）：若选 \(d\)，那么剩余 \(k-1\) 位的可行数量是 \(F(k-1, s+d)\)。
• 按 \(d\) 从小到大累加，找到第一个使得累计覆盖到 \(n\) 的 \(d\)，就确定这一位；
• 若 \(n\) 超过该 \(d\) 的数量，则 \(n \mathrel{-}= F(k-1, s+d)\) 继续试更大的 \(d\)。

由于每一位最多试 \(10\) 个数字，总步数约 \(10L\)。

对于需要的质数范围，如果最多考虑到 \(\ell\) 位，则数位和最大 \(9\cdot \ell\)。只需筛到这个上界即可。本题实际会在很小位数内结束（代码会自动找到足够的 \(L\)），预计算到 \(\ell \ln \ell\) 位也几乎没有开销。

代码

N = 10 ** 16
from tools import Factorizer


def build_F(max_len: int):
    max_sum = 9 * max_len
    primes = set(Factorizer(max_sum).primes)
    F = [[0] * (max_sum + 10) for _ in range(max_len + 1)]
    for s in range(max_sum + 1):
        F[0][s] = int(s in primes)
    for k in range(1, max_len + 1):
        prev = F[k - 1]
        cur = F[k]
        for s in range(max_sum, -1, -1):
            cur[s] = sum(prev[s:s + 10])
    return F, max_sum


def D(n: int) -> int:
    max_len = int(len(str(n)) * 2)
    F, max_sum = build_F(max_len)

    nn = n
    L = 0
    for l in range(1, max_len + 1):
        c = 0
        for d in range(1, 10):
            c += F[l - 1][d]
        if nn > c:
            nn -= c
        else:
            L = l
            break

    rem = L - 1
    s = 0
    digits = []

    for d in range(1, 10):
        ways = F[rem][s + d]
        if nn > ways:
            nn -= ways
        else:
            digits.append(d)
            s += d
            break

    while rem > 0:
        rem -= 1
        for d in range(10):
            ways = F[rem][s + d]
            if nn > ways:
                nn -= ways
            else:
                digits.append(d)
                s += d
                break

    ans = 0
    for d in digits:
        ans = ans * 10 + d
    return ans


print(D(N))