Project Euler 467

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

Project Euler 467

题目

Superinteger

An integer s is called a superinteger of another integer \(n\) if the digits of \(n\) form a subsequence of the digits of \(s\).

For example, \(2718281828\) is a superinteger of \(18828\), while \(314159\) is not a superinteger of \(151\).

Let \(p(n)\) be the \(n\text{th}\) prime number, and let \(c(n)\) be the \(n\text{th}\) composite number. For example, \(p(1) = 2, p(10) = 29, c(1) = 4\) and \(c(10) = 18\).

\(\{p(i) : i \ge 1\} = \{2, 3, 5, 7, 11, 13, 17, 19, 23, 29, \dots\}\)

\(\{c(i) : i \ge 1\} = \{4, 6, 8, 9, 10, 12, 14, 15, 16, 18, \dots\}\)

Let \(P^D\) the sequence of the digital roots of \(\{p(i)\}\) (\(C^D\) is defined similarly for \(\{c(i)\}\)):

\(P^D = \{2, 3, 5, 7, 2, 4, 8, 1, 5, 2, \dots\}\)

\(C^D = \{4, 6, 8, 9, 1, 3, 5, 6, 7, 9, \dots\}\)

Let \(P_n\) be the integer formed by concatenating the first \(n\) elements of \(P^D\) (\(C_n\) is defined similarly for \(C^D\)).

\(P_{10} = 2357248152\)

\(C_{10} = 4689135679\)

Let \(f(n)\) be the smallest positive integer that is a common superinteger of \(P_n\) and \(C_n\).

For example, \(f(10) = 2357246891352679\), and \(f(100) \bmod 1000000007 = 771661825\).

Find \(f(10000) \bmod 1000000007\).

解决方案

可见，数字根有闭式：

\[ \mathrm{dr}(x)=1+((x-1)\bmod 9), (x\ge 1) \]

那么\(P^D = \{\mathrm{dr}(p(i))\}_{i\ge 1},C^D = \{\mathrm{dr}(c(i))\}_{i\ge 1}.\)

由于这里每一位都在 \(1\sim 9\)，不会出现前导 \(0\)。因此：数值最小与“长度最短优先，长度相同则字典序最小”完全一致。

于是问题严格等价为：给定两个长度为 \(n\) 的字符串 \(A=P_n\)、\(B=C_n\)，求它们的 最短公共超序列（Shortest Common Supersequence, SCS）中，字典序最小的那一个。

设 \(A\) 与 \(B\) 长度都为 \(n\)，下标从 \(0\) 开始。定义：\(f(i,j)\)表示字符串\(A[i..n-1]\)与\(B[j..n-1]\)的最短公共超序列长度。那么不难写下状态转移方程：

\[ f(i,j)= \left \{\begin{aligned} &n-j & & \text{if}\quad i=n \\ &n-i & & \text{else if}\quad j=n \\ &f(i+1,j+1)+1 & & \text{else if}\quad A[i] = B[j] \\ &\min(f(i+1,j),f(i,j+1))+1 & & \text{else} \end{aligned}\right. \]

对于第三行的转移，当 \(A[i]=B[j]\)时，这一位必须只取一次；对于第四行的转移，若 \(A[i]\ne B[j]\)，则这一位只能取 \(A[i]\) 或 \(B[j]\)。这样迭代可得到最短长度。

若只求长度，上述 DP 已够。但题目要求“最小的 superinteger”，即最短长度中的字典序最小者。

当 \(A[i]\ne B[j]\) 且两种选择都会得到同样的最短长度\(f(i+1,j) = f(i,j+1)\)时，此时两条路的最短长度都为：$ 1+f(i+1,j)=1+f(i,j+1)$

那么字典序比较只取决于第一位，直接选更小的数字即可：若 \(A[i]<B[j]\) 取 \(A[i]\)否则取 \(B[j]\)。后缀 DP 的好处在于：当我们决定当前第一位以后，剩下部分已经由状态 \((i+1,j)\) 或 \((i,j+1)\) 定义为“同样最短且字典序最小”，因此这条“首位最小”规则是全局正确的。

最终答案字符串长度最多为 \(2n\)，不能直接拼接成整数。

定义\(v(i,j)\) 为该最优 SCS（最短且字典序最小）的数值对 \(M=10^9+7\) 取模的结果

若在状态 \((i,j)\) 选择把数字 \(d\) 放在最前面，后面接最优子问题的字符串，设子问题长度为 \(L\)、模值为 \(R\)，则新模值为：\((d || \text{suffix}) \bmod M=(d\cdot 10^L + R)\bmod M\)。

因此只需要在 DP 过程中同步维护长度与模值，而不需要真正构造字符串。

因此，转移时的“长度 + 模值”一起维护。

当 \(A[i]=B[j]\)时，有\(f(i,j) = 1 + f(i+1,j+1),v = (A[i]\cdot 10^{f(i+1,j+1)} + v(i+1,j+1))\bmod M\)

当 \(A[i]\ne B[j]\)：

  • 取 \(A[i]\) 的候选：\(f_A = 1 + f(i+1,j),v_A = (A[i]\cdot 10^{f(i+1,j)} + v(i+1,j))\bmod M\)
  • 取 \(B[j]\) 的候选：\(f_B = 1 + f(i,j+1),v_B = (B[j]\cdot 10^{f(i,j+1)} + v(i,j+1))\bmod M\)

选择规则则是：若 \(f_A<f_B\) 选 \(A\)；若 \(f_A>f_B\) 选 \(B\)；若 \(f_A=f_B\)，选 \(A[i],B[j]\)更小的那一个（即字典序更小的首位）。

在实现时计算顺序使用后缀 DP，它可以用滚动数组完成。只需要保存“下一行” \(i+1\) 的一维数组，以及当前行 \(i\) 的一维数组即可。

代码

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
103
104
105
106
107
108
109
110
111
112
113
114
#include <bits/stdc++.h>
#include "tools.h"
using namespace std;

typedef long long ll;
const int N = 10000;
ll MOD = 1000000007LL;
pair<vector<int>, vector<int>> build_sequences(int n)
{
    int L = n * log(n) * 2 + 40;
    Factorizer fact(L);
    vector<int> A(n), B(n);
    int pa = 0, pb = 0;
    for (int x = 2; x <= L; x++)
    {
        int d = (x - 1) % 9 + 1;
        if (fact.spf[x] == x)
        {
            if (pa < n)
                A[pa++] = d;
        }
        else
        {
            if (x >= 4 && pb < n)
                B[pb++] = d;
        }
        if (pa == n && pb == n)
            break;
    }
    return {A, B};
}

ll solve(int n)
{
    auto [A, B] = build_sequences(n);

    vector<ll> pow10(2 * n + 1, 1);
    for (int i = 1; i <= 2 * n; i++)
        pow10[i] = (pow10[i - 1] * 10) % MOD;

    vector<int> len_next(n + 1, 0), len_cur(n + 1, 0);
    vector<ll> val_next(n + 1, 0), val_cur(n + 1, 0);

    len_next[n] = 0;
    val_next[n] = 0;
    for (int j = n - 1; j >= 0; j--)
    {
        len_next[j] = len_next[j + 1] + 1;
        val_next[j] = (B[j] * pow10[len_next[j + 1]] + val_next[j + 1]) % MOD;
    }

    for (int i = n - 1; i >= 0; i--)
    {
        int ai = A[i];

        len_cur[n] = len_next[n] + 1;
        val_cur[n] = (ai * pow10[len_next[n]] + val_next[n]) % MOD;

        for (int j = n - 1; j >= 0; j--)
        {
            int bj = B[j];
            if (ai == bj)
            {
                int l = len_next[j + 1];
                len_cur[j] = l + 1;
                val_cur[j] = (ai * pow10[l] + val_next[j + 1]) % MOD;
            }
            else
            {
                int lA = len_next[j] + 1;
                ll vA = (ai * pow10[len_next[j]] + val_next[j]) % MOD;

                int lB = len_cur[j + 1] + 1;
                ll vB = (bj * pow10[len_cur[j + 1]] + val_cur[j + 1]) % MOD;

                if (lA < lB)
                {
                    len_cur[j] = lA;
                    val_cur[j] = vA;
                }
                else if (lA > lB)
                {
                    len_cur[j] = lB;
                    val_cur[j] = vB;
                }
                else
                {
                    if (ai < bj)
                    {
                        len_cur[j] = lA;
                        val_cur[j] = vA;
                    }
                    else
                    {
                        len_cur[j] = lB;
                        val_cur[j] = vB;
                    }
                }
            }
        }

        swap(len_next, len_cur);
        swap(val_next, val_cur);
    }

    return val_next[0];
}

int main()
{
    cout << solve(N) << "\n";
    return 0;
}