题解 P3702 【[SDOI2017]序列计数】

· · 题解

我觉得这篇题解下很多矩阵解法没有写清楚矩阵怎么构造...所以我来详细地写一下

同步更新于我的博客：「LUOGU P3702」[SDOI2017]序列计数

首先这个形式显然非常不好计算，容斥一波：\text{ans}=\text{满足和为}p\text{的倍数的方案数}-\text{满足和为}p\text{的倍数且不含质数的方案数}

我们很容易得到20\rm pts的优秀做法：令f_{i,j}表示i个数\mod p等于j的方案数，则转移为f_{i,j} = \sum_{k} f_{i-1,k} \times \text{cnt}_{(j - k + p)\bmod \; p}，其中\text{cnt}_i表示1\sim m中\bmod P为i的个数。令g_{i,j}表示i个合数\mod p等于j的方案数，则转移为g_{i,j} = \sum_{k} g_{i-1,k} \times \text{compo}_{(j - k + p)\bmod \; p}，其中\text{compo}_i表示1\sim m中的合数\bmod P为i的个数。 考虑用矩阵乘法优化。对于f，矩阵如下所示：

\mathbf P = \begin{bmatrix} \text{cnt}_0 & \text{cnt}_{p-1} & \text{cnt}_{p-2} & \text{cnt}_{p-3} & \text{cnt}_{p-4} & ... & \text{cnt}_1 \\ \text{cnt}_1 & \text{cnt}_0 & \text{cnt}_{p-1} & \text{cnt}_{p-2} & \text{cnt}_{p-3} & ... & \text{cnt}_2 \\ \text{cnt}_2 & \text{cnt}_1 & \text{cnt}_0 & \text{cnt}_{p-1} & \text{cnt}_{p-2} & ... & \text{cnt}_3 \\ \text{cnt}_3 & \text{cnt}_2 & \text{cnt}_1 & \text{cnt}_0 & \text{cnt}_{p-1} & ... & \text{cnt}_4 \\ ... & ... & ... & ... & ... & ... & ...\\ \text{cnt}_{p-1} & \text{cnt}_{p-2} & \text{cnt}_{p-3} & \text{cnt}_{p-4} & \text{cnt}_{p-5} & ... & \text{cnt}_0 \end{bmatrix}

向量长这样：

\mathbf V = \begin{bmatrix} \text{cnt}_0 \\ \text{cnt}_1 \\ \text{cnt}_2 \\ ... \\ \text{cnt}_{p-1} \end{bmatrix}

最后答案矩阵为\mathbf F = \mathbf P ^ {n-1}\times \mathbf V。 对于g，矩阵如下所示：

\mathbf Q = \begin{bmatrix} \text{compo}_0 & \text{compo}_{p-1} & \text{compo}_{p-2} & \text{compo}_{p-3} & \text{compo}_{p-4} & ... & \text{compo}_1 \\ \text{compo}_1 & \text{compo}_0 & \text{compo}_{p-1} & \text{compo}_{p-2} & \text{compo}_{p-3} & ... & \text{compo}_2 \\ \text{compo}_2 & \text{compo}_1 & \text{compo}_0 & \text{compo}_{p-1} & \text{compo}_{p-2} & ... & \text{compo}_3 \\ \text{compo}_3 & \text{compo}_2 & \text{compo}_1 & \text{compo}_0 & \text{compo}_{p-1} & ... & \text{compo}_4 \\ ... & ... & ... & ... & ... & ... & ...\\ \text{compo}_{p-1} & \text{compo}_{p-2} & \text{compo}_{p-3} & \text{compo}_{p-4} & \text{compo}_{p-5} & ... & \text{compo}_0 \end{bmatrix}

向量长这样：

\mathbf W = \begin{bmatrix} \text{compo}_0 \\ \text{compo}_1 \\ \text{compo}_2 \\ ... \\ \text{compo}_{p-1} \end{bmatrix}

最后答案矩阵为\mathbf G = \mathbf Q ^ {n-1}\times \mathbf W。 最终答案为\mathbf F_{1,1} - \mathbf G_{1,1}

#include <bits/stdc++.h>
using namespace std;
const int MAXN = 1e2 + 10, mo = 20170408, MAXM = 2e7 + 10;
struct mat
{
    int m, n, ma[MAXN][MAXN];
    mat() {}
    mat(int _m, int _n) :
        m(_m), n(_n)
    {
        memset(ma, 0, sizeof(ma));
    }
    friend inline mat operator * (mat a, mat b)
    {
        mat res = mat(a.m, b.n);
        for(int i = 1; i <= res.m; i++)
            for(int j = 1; j <= res.n; j++)
                for(int k = 1; k <= a.n; k++)
                    res.ma[i][j] = (res.ma[i][j] + 1ll * a.ma[i][k]
                                    * b.ma[k][j] % mo) % mo;
        return res;
    }
    friend inline mat operator ^ (mat a, int b)
    {
        mat c = a, res = mat(a.m, a.n);
        for(int i = 1; i <= res.m; i++)
            res.ma[i][i] = 1;
        while(b)
        {
            if(b & 1) res = res * c;
            c = c * c;
            b >>= 1;
        }
        return res;
    }
} P, Q, V, W;
int n, m, p, cnt[MAXN], compo[MAXN];
bool prime[MAXM];
int pr[MAXM];
template <class T>
inline void _read(T &x)
{
    x = 0;
    char t = getchar();
    while(!isdigit(t) && t != '-') t = getchar();
    if(t == '-')
    {
        _read(x);
        x *= -1;
        return ;
    }
    while(isdigit(t))
    {
        x = x * 10 + t - '0';
        t = getchar();
    }
}
int ptot = 0;
int main()
{
    _read(n), _read(m), _read(p);
    for(int i = 0; i < p; i++) cnt[i] = m / p;
    for(int i = 1; i <= m % p; i++) cnt[i]++;
    P = mat(p, p);
    P.ma[1][1] = cnt[0];
    for(int i = 2; i <= p; i++) P.ma[1][i] = cnt[p - i + 1];
    for(int i = 2; i <= p; i++)
    {
        for(int j = 2; j <= p; j++)
            P.ma[i][j] = P.ma[i - 1][j - 1];
        P.ma[i][1] = P.ma[i - 1][p];
    }
    memset(prime, 1, sizeof prime);
    prime[1] = false;
    for(int i = 2; i <= m; i++)
    {
        if(prime[i] == true)
            pr[++ptot] = i;
        for(int j = 1; i * pr[j] <= m && j <= ptot; ++j)
        {
            prime[i * pr[j]] = false;
            if(i % pr[j] == 0)
                break;
        }
    }
    for(int i = 1; i <= m; i++)
        if(!prime[i]) compo[i % p]++;
    Q = mat(p, p);
    Q.ma[1][1] = compo[0];
    for(int i = 2; i <= p; i++) Q.ma[1][i] = compo[p - i + 1];
    for(int i = 2; i <= p; i++)
    {
        for(int j = 2; j <= p; j++)
            Q.ma[i][j] = Q.ma[i - 1][j - 1];
        Q.ma[i][1] = Q.ma[i - 1][p];
    }
    for(int i = 1; i <= p; i++)
        for(int j = 1; j <= p; j++)
            P.ma[i][j] %= mo, Q.ma[i][j] %= mo;
    V = mat(p, 1);
    W = mat(p, 1);
    for(int i = 1; i <= p; i++) V.ma[i][1] = cnt[i - 1] % mo;
    for(int i = 1; i <= p; i++) W.ma[i][1] = compo[i - 1] % mo;
    P = (P ^ n - 1) * V;
    Q = (Q ^ n - 1) * W;
    printf("%d\n", (P.ma[1][1] - Q.ma[1][1] + mo) % mo);
    return 0;
}