题解：P3702 [SDOI2017] 序列计数

· · 题解

>题面链接<

更佳观感： \mathtt{My \ Blog} ——题解：P3702 [SDOI2017] 序列计数

首先看到“至少有一个数是质数”就觉得不好算，一眼容斥成 n 个数的和是 p 的倍数的方案数减去 n 个非质数的和是 p 的倍数的方案数。

那我们就考虑 n 个数的和是 p 的倍数的方案数怎么求（另一个再加个筛素数就可以了）。

状态定义： f_{i,j} 表示 i 个数的和模 p 为 j 的方案数。

转移： 定义 \text{cnt}_x = \displaystyle\sum_{i \in [1..m]}\big[i \bmod p = x\big]

那么 f_{i,j} = \displaystyle\sum_{k \in [0..p)} f_{i - 1, k} \cdot \text{cnt}_{j - k}

我们发现这个东西好像可以矩阵快速幂！

因为 f_{1, i} = \displaystyle \text{cnt}_i，所以构造

\mathbf{A} = \begin{bmatrix} \text{cnt}_0 \\ \text{cnt}_1 \\ \vdots \\ \text{cnt}_{p-1} \\ \end{bmatrix}_{n \times 1}

随后假设已经求出来了 \begin{bmatrix} f_{i-1,0} \\ f_{i-1,1} \\ \vdots \\ f_{i-1,p-1} \\ \end{bmatrix}_{n \times 1}

那么相应的就应构造

\mathbf{B} = \begin{bmatrix} \text{cnt}_{0 - 0} & \text{cnt}_{0 - 1} & \cdots & \text{cnt}_{0 - (p-1)} \\ \text{cnt}_{1 - 0} & \text{cnt}_{1 - 1} & \cdots & \text{cnt}_{1 - (p-1)} \\ \cdots & \cdots & \cdots & \cdots \\ \text{cnt}_{(p - 1) - 0} & \text{cnt}_{(p - 1) - 1} & \cdots & \text{cnt}_{(p - 1) - (p - 1)}\end{bmatrix}_{n \times n}

同理将 \text{cnt} 换成 \text{cnt}_x = \displaystyle\sum_{i \in [1..m]}\big[i \bmod p = x \land \text{is\_not\_prime}(i) \big] 可以求出 n 个非质数的和是 p 的倍数的方案数。

然后写✍代码~~~

代码

\textbf{\tiny{如直接复制代码CE后果自负}}

#include <bits/stdc++.h>‍‌
using namespace std;‭‌
typedef long long LL;‭‌
const int N = 100 + 10, M = 2e7 + 10, mod = 20170408;‭‌
int n, m, p, primes[M], cnt, cnt1[N], cnt2[N];‌
bool st[M];‭‌
struct Matrix{‭‌
    int m, n, mat[N][N];‌
    Matrix(‭) {‭}‌
    Matrix(int _m, int _n) {‭
        m = _m, n = _n;‍‌
        memset(mat, 0, sizeof mat);‍‌
    }‭‍‌
    inline Matrix operator * (Matrix b) {‭
        if(n != b.m) throw "Invalid Matrix Operation.";‍‌
        Matrix res(Matrix(m, b.n));‍‍‌‌
        for(int i = 0; i < m; i ++)‍‍‌‌
            for(int j = 0; j < b.n; j ++)‍‍‌‌
                for(int k = 0; k < n; k ++)‍‍‌‌
                    res.mat[i][j] = (res.mat[i][j] + (LL)mat[i][k] * b.mat[k][j] % mod) % mod;‍‌
        return res;‍‍‌‌
    }‭
    inline Matrix operator ^ (int k) {‭
        if(m != n) throw "Invalid Matrix Operation.";
        Matrix tmp(*this), res(Matrix(m, n));‭
        for(int i = 0; i < m; i ++) res.mat[i][i] = 1;‭
        while(k) {‭
            if(k & 1) res = res * tmp;
            tmp = tmp * tmp;
            k >>= 1;
        }‭‍‌
        return res;‍‌
    }‭
}A, B;‭‍‌
int get_ans(int *cnt) {‭
    A = Matrix(p, p);‭
    for(int i = 0; i < p; i ++)‭
        for(int j = 0; j < p; j ++)‭
            A.mat[i][j] = cnt[(i - j + p) % p];‭
    B = Matrix(p, 1);‭
    for(int i = 0; i < p; i ++) B.mat[i][0] = cnt[i];‭
    return ((A ^ (n - 1)) * B).mat[0][0];‭
}
int main(‍‌)‍‍‌‌
{‍‌
    scanf("%d%d%d", &n, &m, &p);‭
    for(int i = 1; i <= m; i ++) cnt1[i % p] ++;
    st[1] = true;‭ // 自然1不是质数
    for(int i = 1; i <= m; i ++) {‭
        if(!st[i]) primes[++ cnt] = i;‭
        for(int j = 1; j <= cnt && primes[j] <= m / i; j ++) {‭
            st[i * primes[j]] = true;
            if(!(i % primes[j])) break;
        }‭
    }‭
    for(int i = 1; i <= m; i ++)
        if(st[i])‭
            cnt2[i % p] ++;‭
    printf("%d\n", (get_ans(cnt1) - get_ans(cnt2) + mod) % mod);‭
    return 0;‭
}‭