这是一道数学分析背景的数位 DP 经典题……

为什么 BZOJ 上没人做呢……

题意

BZOJ 3452：给定 $n$ 个不含前导零的数字串 $s_1, s_2, \cdots, s_n$，对于所有满足十进制表示中不含任何一个给定数字串的正整数 $x$，求它们的倒数和，四舍五入精确到小数点后第四位。保证 $n, |s_i| \leq 3~(i = 1, 2, \cdots, n)$。

ProjectEuler 368：对于所有满足十进制表示中不含连续三个相同数位的正整数 $x$，求它们的倒数和，四舍五入精确到小数点后第十位。

注意：Tyvj 上的数据是错的

分析

对于任意一个正整数 $x$，我们可以利用字符串匹配检查它是否对答案产生贡献，而且可以发现这个正整数的任意前缀数字都是可以产生贡献的，因此可以考虑数位 DP 套字符串匹配，按照数位划分状态的同时将字符串匹配的信息存储在状态中。

由于题目涉及到多个字符串间的匹配，因此对于 BZOJ 3452 我们采用 AC 自动机进行字符串匹配，而对于 ProjectEuler 368 我们则只需要记录最后一个数位以及它是否与倒数第二个数位相同。这里假设读者已经具备相关知识储备。

在数位 DP 时，我们需要考虑在一些不含前导 $0$ 的数字 $x~(x \in S)$ 后添加一个数位 $d$ 并维护倒数和，注意到有 $$\frac{1}{10 x + d} = \sum_{k \geq 1}{\frac{(-d)^{k - 1}}{(10 x)^k}},$$ 于是有 $$\sum_{x \in S}{\frac{1}{(10 x + d)^n}} = \sum_{k \geq n}{{k - 1 \choose n - 1} \frac{(-d)^{k - n}}{10^k} \sum_{x \in S}{\frac{1}{x^k}}},$$ 那么当 $x$ 和 $k$ 足够大时， $\sum\limits_{x \in S}{\dfrac{1}{x^k}}$ 是极小的，因此可以为 $x$ 确定长度下界、为 $k$ 确定上界，对答案进行近似计算。

具体来说，注意到 $x = 1$ 时 $k$ 无法取到足够好的上界，因此我们可以设定阈值 $L$，对于 $x < 10^L$ 的情况采用枚举的方法，对于 $x \geq 10^L$ 的情况采用数位 DP 近似。设定 $L$ 后，可以直观地设定 $k$ 的上限 $K$，只需让 $\sum\limits_{x \in S}{\dfrac{1}{x^K}}$ 的最大值也远远小于答案允许的误差即可。

现在继续考虑如何用数位 DP 计算 $x \geq 10^L$ 的贡献，我们定义 $f(i, j, k)$ 表示长度为 $i$、匹配状态为 $j$ 的所有能够产生贡献的数字的 $k$ 次方倒数和，若可行的匹配状态集合为 $P$，则对答案的贡献为 $$\sum_{i > L} \sum_{j \in P}{f(i, j, 1)},$$ 可以再次设定长度的上界进行近似计算，不过这里为了保证精度更好可以采用更好的方法。

设 $M = |P|$，数位 DP 的初始状态 $f(L, j, k)$ 可以写成一个长度为 $M K$ 的行向量 $A$，而 $f(i + 1, j, k)$ 与 $f(i, j, k)$ 之间的转移也可以写成一个大小为 $M K \times M K$ 的矩阵 $B$，则 $x \geq 10^L$ 的全部 DP 信息为 $$\sum_{k \geq 1}{A B^k} = A B (I - B)^{-1},$$ 其中 $I$ 是大小为 $M K \times M K$ 的单位矩阵，所求即结果矩阵的 $M$ 个特定元素之和。

于是本题就在 $\mathcal{O}(10^L + M^3 K^3)$ 的时间复杂度中解决了。

代码

贴一个 ProjectEuler 368 的代码仅供交流学习。

#include <bits/stdc++.h>
typedef double DB;
const int maxd = 10, maxs = 101, maxl = 6;
int m, e, sz;
DB lpw[maxd + 1][maxs], c[maxs][maxs];
DB f[maxs], mat[maxs][maxs << 1 | 1], ans;
void dfs(int dep, int val, int msk) {
    if(dep == maxl) {
        DB prd = 1;
        for(int i = 0; i < e; ++i)
            f[i * m + msk] += 1 / (prd *= val);
        return;
    }
    ans += 1.0 / val;
    for(int i = 0; i < maxd; ++i) {
        if(msk & 1 && (msk >> 1) == i)
            continue;
        dfs(dep + 1, val * maxd + i, i << 1 | ((msk >> 1) == i));
    }
}
int main() {
    m = maxd << 1;
    e = (maxs - 1) / m;
    sz = e * m;
    lpw[0][0] = 0;
    for(int i = 1; i <= maxd; ++i) {
        lpw[i][0] = 0;
        lpw[i][1] = log(i);
        for(int j = 2; j < e; ++j)
            lpw[i][j] = lpw[i][j - 1] + lpw[i][1];
    }
    lpw[maxd][e] = lpw[maxd][e - 1] + lpw[maxd][1];
    for(int i = 0; i < e; ++i) {
        c[i][0] = c[i][i] = 1;
        for(int j = 1; j < i; ++j)
            c[i][j] = c[i - 1][j - 1] + c[i - 1][j];
    }
    for(int i = 0; i < sz; ++i)
        mat[i][i] = mat[i][sz + i] = 1;
    for(int i = 0; i < m; ++i)
        for(int d = 0; d < maxd; ++d) {
            if(i & 1 && (i >> 1) == d)
                continue;
            int ii = d << 1 | ((i >> 1) == d);
            for(int j = 0; j < e; ++j)
                for(int jj = 0; jj <= j; ++jj) {
                    if(!d && jj < j)
                        continue;
                    mat[j * m + i][jj * m + ii] -= ((j - jj) & 1 ? -1 : 1) * exp(log(c[j][jj]) + lpw[d][j - jj] - lpw[maxd][j + 1]);
                }
        }
    for(int i = 0, j, k; i < sz; ++i) {
        for(j = i, k = j + 1; k < sz; ++k)
            if(fabs(mat[k][i]) > fabs(mat[j][i]))
                j = k;
//        assert(fabs(mat[j][i]) > 1e-20);
        DB prd = mat[j][i];
        for(k = i; k < sz << 1; ++k) {
            std::swap(mat[j][k], mat[i][k]);
            mat[i][k] /= prd;
        }
        for(j = 0; j < sz; ++j) {
            if(i == j)
                continue;
            DB prd = mat[j][i];
            for(int k = i; k < sz << 1; ++k)
                mat[j][k] -= prd * mat[i][k];
        }
    }
    for(int i = 1; i < maxd; ++i)
        dfs(1, i, i << 1);
    for(int i = 0; i < sz; ++i)
        for(int j = 0; j < m; ++j)
            ans += f[i] * mat[i][sz + j];
    printf("%.10f\n", (double)ans);
    return 0;
}

Written with StackEdit.