LOj #2002. 「SDOI2017」序列计数 （容斥+dp+矩阵快速幂）

题目链接：
LOj 2002

题意：
要求得到一个长度为 $n$ 的序列，序列中的数都是不超过 $m$ 的正整数，而且这 $n$个数的和是 $p$ 的倍数。这$n$个数中，至少有一个数是质数。问你有多少个序列满足要求。

题解：
根据容斥原理，因为题目要求至少有一个数是素数，用所有方案减去不含质数的方案就是答案。
设 $dp[i][j]$表示序列前 $i$ 个数模 $p$ 的余数为 $j$ 时的方案数。$dp[i][j]=∑dp[i−1] [(j−k) \mod p]$。
如果不考虑包含质数的情况，因为$p$很小，我们不妨构建$p*p$的矩阵，然后对于$(i,j)$填有多少个数$x$使得$(i+x)\%p=j$，最后矩阵快速幂一下就行了。对于不含质数的情况，我们也同样构建这样一个矩阵，然后矩阵快速幂一下就可以了。

AC代码：

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int mod = 20170408;
const int sz = 110;
int n,m,p;
struct matrix
{
    int ma[sz][sz];
    matrix friend operator *(const matrix a,const matrix b)
    {
        matrix ret;
        memset(ret.ma,0,sizeof(ret.ma));
        for(int i=0;i<p;i++){
            for(int j=0;j<p;j++)
            {
                for(int k=0;k<p;k++)
                {
                    ret.ma[i][j] = (ret.ma[i][j] + 1LL * a.ma[i][k] * b.ma[k][j]%mod)%mod;
                }
            }
        }
        return ret;
    }
}A,ans;


matrix multipow(matrix x, int k)
{
    memset(ans.ma,0,sizeof(ans.ma));
    for(int i=0;i<p;i++){
        ans.ma[i][i] = 1;
    }
    while(k)
    {
        if(k&1)ans=ans*x;
        k>>=1;
        x=x*x;
    } 
    return ans;
}

bool isprime[20000010];
int prime[2000010];

void sieve()
{
    isprime[0] = true; isprime[1] = true;
    for(int i=2,cnt=0;i<=m;i++){
        if(isprime[i]==false) prime[++cnt] = i;
        for(int j=1;j<=cnt&&i*prime[j]<=m;j++)
        {
            isprime[i*prime[j]] = true;
            if(i%prime[j]==0)break;
        }
    }
}

int dp[110];

int solve1()
{
    for (int i=1;i<=m;i++)  
    {
        dp[i%p]++;
    }
    for (int j=1;j<=m;j++)
    {   
        A.ma[(-j%p+p)%p][0]++;
    }
    for (int i=1;i<p;i++)
    {
        for (int j=0;j<p;j++)
        {
            A.ma[j][i] = A.ma[(j-1+p)%p][i-1];
        }         
    }
//  for(int i=1;i<=m;i++)
//  {
    //  for(int j=1;j<=m;j++){
    //      printf("%d ",A.ma[i][j]);
    //  }
    //  cout<<endl;
//  }
    A=multipow(A,n-1);
    int ans=0;
    for (int i=0;i<p;i++) 
    {
        ans = (ans + 1LL * dp[i] * A.ma[i][0] % mod ) % mod;
    }
    //cout<<"ans1="<<ans<<endl;
    return ans;
}

int solve2()
{
    memset(dp,0,sizeof(dp));
    for (int i=1;i<=m;i++)
    {
        if (isprime[i]) 
        {
            dp[i%p]++;
        }
    }

    memset(A.ma,0,sizeof(A.ma));

    for (int j=1;j<=m;j++)
    {
        if (isprime[j]) 
        {
            A.ma[(-j%p+p)%p][0]++;
        }
    }

    for (int i=1;i<p;i++)
    {
        for (int j=0;j<p;j++)
        {
            A.ma[j][i] = A.ma[(j-1+p)%p][i-1];
        }           
    }
    //for(int i=1;i<=m;i++)
//  {
    //  for(int j=1;j<=m;j++){
    //      printf("%d ",A.ma[i][j]);
    //  }
    //  cout<<endl;
//  }       
    A=multipow(A,n-1);
    int ans=0;
    for (int i=0;i<p;i++) 
    {
        ans = (ans + 1LL * dp[i] * A.ma[i][0] % mod ) % mod;
    }
   // cout<<"ans2="<<ans<<endl;
    return ans;
}
/*
59 74 92
12897479
*/
int main()
{
    //freopen("in.txt","r",stdin);
    scanf("%d%d%d",&n,&m,&p);
    sieve();
    int x=solve1();
    int y=solve2();
    printf("%d\n",(x-y+mod)%mod);
    return 0; 
 }