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)modp]。
如果不考虑包含质数的情况，因为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 9212897479*/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;  }