[BZOJ4818][SDOI2017]序列计数（DP+容斥原理+矩乘）

最新推荐文章于 2019-10-22 18:13:47 发布

原创最新推荐文章于 2019-10-22 18:13:47 发布 · 1k 阅读

1 ·

CC 4.0 BY-SA版权

BZOJ、UOJ、LOJ 专栏收录该内容

341 篇文章

订阅专栏

这篇博客介绍了如何利用动态规划(DP)和容斥原理解决[BZOJ4818][SDOI2017]序列计数问题。首先通过容斥原理，用总方案数减去不含质数的方案数来处理「至少一个质数」的条件。接着进行预处理，计算不同模数下质数和合数的数量。然后建立朴素DP方程，并推导矩阵优化方法，以应对大规模数据。最终结果为经过矩阵运算后的f[n][0]减去g[n][0]。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

1、大意做法

首先发现， $DP$ 的统计方案数中，「至少一个质数」的条件不容易直接 $DP$ ，但是「不含任意一个质数」的条件容易 $DP$ 。所以这里使用容斥原理，用总方案数减去不含任意一个质数的方案数。

2、预处理

$cnt_0,cnt_1,cnt_2,...cnt_{p-1}$ 分别为 $1...m$ 的数中模 $p$ 为 $0,1,2,...p-1$ 的数的个数。
$pri_0,pri_1,pri_2,...pri_{p-1}$ 分别为 $1...m$ 的数中的合数模 $p$ 为 $0,1,2,...p-1$ 的数的个数。

3、建立朴素DP

设 $f[i][j]$ 表示 $i$ 个数的序列，总和模 $p$ 为 $j$ 的方案总数，
$g[i][j]$ 为 $i$ 个数的序列，总和模 $p$ 为 $j$ 且不含质数的方案总数。
可以推出 $DP$ 方程：
$f[i][j]=\sum_{k=0}^{p-1}(f[i-1][(j-k+p) \mod p]*cnt_k)$
$g[i][j]=\sum _{k=0}^{p-1}(g[i-1][(j-k+p) \mod p]*pri_k)$
边界为 $f[1][j]=cnt_j$ ， $g[1][j]=pri_j$
最后结果为 $f[n][0]-g[n][0]$ 。

4、矩阵优化

可以发现，上面的朴素 $DP$ 无法承受 $n<=10^9$ 的数据范围。
这里考虑构建 $p*p$ 的矩阵 $P$ ， $Q$ ， $p*1$ 的矩阵 $V$ ， $W$ 。
矩阵 $P$ 的内容为：

⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ c n t 0 c n t 1 c n t 2 c n t 3 . . . c n t p - 1 c n t p - 1 c n t 0 c n t 1 c n t 2 . . . c n t p - 2 c n t p - 2 c n t p - 1 c n t 0 c n t 1 . . . c n t p - 3 c n t p - 3 c n t p - 2 c n t p - 1 c n t 0 . . . c n t p - 4 c n t p - 4 c n t p - 3 c n t p - 2 c n t p - 1 . . . c n t p - 5 . . . . . . . . . . . . . . . . . . c n t 1 c n t 2 c n t 3 c n t 4 . . . c n t 0 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥

$\begin{bmatrix} cnt_0 & cnt_{p-1} & cnt_{p-2} & cnt_{p-3} & cnt_{p-4} & ... & cnt_1 \\ cnt_1 & cnt_0 & cnt_{p-1} & cnt_{p-2} & cnt_{p-3} & ... & cnt_2 \\ cnt_2 & cnt_1 & cnt_0 & cnt_{p-1} & cnt_{p-2} & ... & cnt_3 \\ cnt_3 & cnt_2 & cnt_1 & cnt_0 & cnt_{p-1} & ... & cnt_4 \\ ... & ... & ... & ... & ... & ... & ...\\ cnt_{p-1} & cnt_{p-2} & cnt_{p-3} & cnt_{p-4} & cnt_{p-5} & ... & cnt_0 \end{bmatrix}$
矩阵

Q $Q$ 的内容为：

⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ p r i 0 p r i 1 p r i 2 p r i 3 . . . p r i p - 1 p r i p - 1 p r i 0 p r i 1 p r i 2 . . . p r i p - 2 p r i p - 2 p r i p - 1 p r i 0 p r i 1 . . . p r i p - 3 p r i p - 3 p r i p - 2 p r i p - 1 p r i 0 . . . p r i p - 4 p r i p - 4 p r i p - 3 p r i p - 2 p r i p - 1 . . . p r i p - 5 . . . . . . . . . . . . . . . . . . p r i 1 p r i 2 p r i 3 p r i 4 . . . p r i 0 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥

$\begin{bmatrix} pri_0 & pri_{p-1} & pri_{p-2} & pri_{p-3} & pri_{p-4} & ... & pri_1 \\ pri_1 & pri_0 & pri_{p-1} & pri_{p-2} & pri_{p-3} & ... & pri_2 \\ pri_2 & pri_1 & pri_0 & pri_{p-1} & pri_{p-2} & ... & pri_3 \\ pri_3 & pri_2 & pri_1 & pri_0 & pri_{p-1} & ... & pri_4 \\ ... & ... & ... & ... & ... & ... & ...\\ pri_{p-1} & pri_{p-2} & pri_{p-3} & pri_{p-4} & pri_{p-5} & ... & pri_0 \end{bmatrix}$
矩阵

V $V$ 的内容为：

⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ c n t 0 c n t 1 c n t 2 . . . c n t p - 1 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥

$\begin{bmatrix} cnt_0 \\ cnt_1 \\ cnt_2 \\ ... \\ cnt_{p-1} \end{bmatrix}$
矩阵

W $W$ 的内容为：

⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ p r i 0 p r i 1 p r i 2 . . . p r i p - 1 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥

$\begin{bmatrix} pri_0 \\ pri_1 \\ pri_2 \\ ... \\ pri_{p-1} \end{bmatrix}$
然后进行矩阵乘方：

F=Pn−1∗V $F=P^{n-1}*V$

G=Qn−1∗W $G=Q^{n-1}*W$
最后结果即为

F[1][1]−G[1][1] $F[1][1]-G[1][1]$ 。
以上操作注意取模。
代码：

#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
const int N = 105, CYX = 20170408, M = 2e7 + 5;
struct cyx {
    int m, n, v[N][N];
    cyx() {}
    cyx(int _m, int _n) :
        m(_m), n(_n) {memset(v, 0, sizeof(v));}
    friend inline cyx operator * (cyx a, cyx b) {
        cyx res = cyx(a.m, b.n);
        int i, j, k;
        for (i = 1; i <= res.m; i++)
            for (j = 1; j <= res.n; j++)
                for (k = 1; k <= a.n; k++)
                    res.v[i][j] = (res.v[i][j] + 1ll * a.v[i][k]
                        * b.v[k][j] % CYX) % CYX;
        return res;
    }
    friend inline cyx operator ^ (cyx a, int b) {
        cyx c = a, res = cyx(a.m, a.n);
        int i, j; for (i = 1; i <= res.m; i++)
            res.v[i][i] = 1;
        while (b) {
            if (b & 1) res = res * c;
            c = c * c;
            b >>= 1;
        }
        return res;
    }
} P, Q, V, W;
int n0, m0, p0, cnt[N], cnt0[N];
bool prime[M];
inline int read() {
    int res = 0; bool bo = 0; char c;
    while (((c = getchar()) < '0' || c > '9') && c != '-');
    if (c == '-') bo = 1; else res = c - 48;
    while ((c = getchar()) >= '0' && c <= '9')
        res = (res << 3) + (res << 1) + (c - 48);
    return bo ? ~res + 1 : res;
}
int main() {
    int i, j; n0 = read(); m0 = read(); p0 = read();
    for (i = 0; i < p0; i++) cnt[i] = m0 / p0;
    for (i = 1; i <= m0 % p0; i++) cnt[i]++;
    P = cyx(p0, p0); P.v[1][1] = cnt[0];
    for (i = 2; i <= p0; i++) P.v[1][i] = cnt[p0 - i + 1];
    for (i = 2; i <= p0; i++) {
        for (j = 2; j <= p0; j++)
            P.v[i][j] = P.v[i - 1][j - 1];
        P.v[i][1] = P.v[i - 1][p0];
    }
    memset(prime, true, sizeof(prime));
    prime[0] = prime[1] = 0;
    for (i = 2; i <= m0; i++) {
        if (!prime[i]) continue;
        if (i * i > m0) break;
        for (j = i * i; j <= m0; j += i)
            prime[j] = 0;
    }
    for (i = 1; i <= m0; i++)
        if (!prime[i]) cnt0[i % p0]++;
    Q = cyx(p0, p0); Q.v[1][1] = cnt0[0];
    for (i = 2; i <= p0; i++) Q.v[1][i] = cnt0[p0 - i + 1];
    for (i = 2; i <= p0; i++) {
        for (j = 2; j <= p0; j++)
            Q.v[i][j] = Q.v[i - 1][j - 1];
        Q.v[i][1] = Q.v[i - 1][p0];
    }
    for (i = 1; i <= p0; i++) for (j = 1; j <= p0; j++)
        P.v[i][j] %= CYX, Q.v[i][j] %= CYX;
    V = cyx(p0, 1); W = cyx(p0, 1);
    for (i = 1; i <= p0; i++) V.v[i][1] = cnt[i - 1] % CYX;
    for (i = 1; i <= p0; i++) W.v[i][1] = cnt0[i - 1] % CYX;
    P = (P ^ n0 - 1) * V; Q = (Q ^ n0 - 1) * W;
    printf("%d\n", (P.v[1][1] - Q.v[1][1] + CYX) % CYX);
    return 0;
}