bzoj4818 [Sdoi2017]序列计数

传送门

矩阵优化dp、容斥原理。
先写出dp转移方程：f[i&1][(j+k)%p][1]=f[(i&1)^1][j][0]+f[(i&1)^1][j][1] (这是我当时在考场上写的20分暴力)
上面的方程对有无素数进行了分类，其实可以进行容斥，用所有方案数减去没有素数的方案数。又因为题目要求序列之和为p的倍数，所以上面的方程对序列之和模p的余数进行了存储。而余数让我们想到了什么呢？矩阵！矩阵可以对不同的余数进行记录，并且可以很好地处理转移。题目中p的最大值为100，所以我们可以开100*100的矩阵，矩阵中存储余数为某一定值的时候的方案数。（这种方法是我自己yy的，所以我还是举个例子比较好）

以p=4为例：

上面左边的矩阵是转移矩阵，中间的是初始矩阵，右边的是结果矩阵，转移矩阵中写的数字是序列之和对p取模后的余数，但事实上转移矩阵中存储的是序列之和为该余数时的方案数，初始矩阵和结果矩阵中的f[i]表示序列之和模p后余数为i的方案数
我们可以看出，通过将初始矩阵不断地用转移矩阵进行矩阵乘法，就可以得到将序列增长后可行的方案数。所以我们对转移矩阵进行n-1次乘法，然后乘上初始矩阵就可以得出答案。这样求出所有方案数和没有素数的方案数就好了。

CODE：

#include<cstdio>#include<cstring>#define mod 20170408#define N 20000005int prime[1500000];bool b[N];int n,m,p,tot,ans;struct Matrix{    int a[105][105];    Matrix(){memset(a,0,sizeof(a));}    inline Matrix operator *(const Matrix &x)const    {        Matrix ans;        for(int i=0;i<p;i++)          for(int j=0;j<p;j++)            if(a[i][j]) for(int k=0;k<p;k++)              ans.a[i][k]=(1ll*ans.a[i][k]+1ll*a[i][j]*x.a[j][k])%mod;        return ans;    }}m1,m2,tmp1,tmp2;inline void euler(){    for(int i=2;i<=m;i++)    {        if(!b[i]) prime[++tot]=i;        for(int j=1;j<=tot&&i*prime[j]<=m;j++)        {            b[i*prime[j]]=1;            if(i%prime[j]==0) break;        }    }}inline void Matrix_init(){    int num=m/p,rest=m%p;    for(int i=0;i<p;i++)      m1.a[0][i]=num;    for(int i=p-1;rest;i--,rest--)      m1.a[0][i]++;    for(int i=1;i<p;i++)    {        m1.a[i][0]=m1.a[i-1][p-1];        for(int j=1;j<p;j++)          m1.a[i][j]=m1.a[i-1][j-1];    }    m2=m1;    for(int i=1;i<=tot;i++)    {        int rest=prime[i]%p;        int pos=p-rest;        if(!rest) pos=0;        for(int j=0;j<p;j++)        {            m2.a[j][pos]--;            pos++;            if(pos==p) pos=0;        }    }    tmp1.a[0][0]=m1.a[0][0];    tmp2.a[0][0]=m2.a[0][0];    for(int i=1,pos=p-1;i<p;i++,pos--)      tmp1.a[i][0]=m1.a[0][pos],      tmp2.a[i][0]=m2.a[0][pos];}inline Matrix ksm(Matrix a,int b){    Matrix ans=a;    for(b--;b;b>>=1,a=a*a)      if(b&1) ans=ans*a;    return ans;}int main(){    scanf("%d%d%d",&n,&m,&p);    euler(),Matrix_init();    m1=ksm(m1,n-1),m2=ksm(m2,n-1);    m1=m1*tmp1,m2=m2*tmp2;    ans=m1.a[0][0]-m2.a[0][0];    if(ans<0) ans+=mod;    printf("%d",ans);    return 0;}