BZOJ 4818: [Sdoi2017]序列计数 （动态规划+矩阵乘法）

来源：http://blog.csdn.net/qq_33229466/article/details/70055284
题目：http://www.lydsy.com/JudgeOnline/problem.php?id=4818
分析

一眼容斥，用所有方案减去不含质数的方案。
设f[i,j]表示序列前i个数模p的余数为j时的方案数。f[i,j]=∑f[i−1,(j−k)modp]
构建矩阵然后快速幂即可。
求不含质数的方案同理。
在一开始建矩阵的时候我的复杂度是mp，发现会超时。但是注意到矩阵的每一列都是循环同构的，于是就可以先O(m)处理好一列，然后再构建其他列即可。
代码：

#include<iostream>#include<cstdio>#include<cstdlib>#include<cstring>#include<algorithm>using namespace std;typedef long long LL;const int M=20000005;const int P=105;const int MOD=20170408;int n,m,p,prime[M/10],tot,f[P];struct arr{int a[P][P];}ans,a;bool not_prime[M];void get_prime(int n){    not_prime[1]=1;    for (int i=2;i<=n;i++)    {        if (!not_prime[i]) prime[++tot]=i;        for (int j=1;j<=tot&&i*prime[j]<=n;j++)        {            not_prime[i*prime[j]]=1;            if (i%prime[j]==0) break;        }    }}void mul(arr &c,arr a,arr b){    memset(c.a,0,sizeof(c.a));    for (int i=0;i<p;i++)        for (int j=0;j<p;j++)            for (int k=0;k<p;k++)                c.a[i][j]=(c.a[i][j]+(LL)a.a[i][k]*b.a[k][j]%MOD)%MOD;}arr ksm(arr x,int y){    memset(ans.a,0,sizeof(ans.a));    for (int i=0;i<p;i++) ans.a[i][i]=1;    while (y)    {        if (y&1) mul(ans,ans,x);        mul(x,x,x);y>>=1;    }    return ans;}int solve1(){    for (int i=1;i<=m;i++) f[i%p]++;    for (int j=1;j<=m;j++) a.a[(-j%p+p)%p][0]++;    for (int i=1;i<p;i++)        for (int j=0;j<p;j++)            a.a[j][i]=a.a[(j-1+p)%p][i-1];    a=ksm(a,n-1);    int ans=0;    for (int i=0;i<p;i++) ans=(ans+(LL)f[i]*a.a[i][0]%MOD)%MOD;    return ans;}int solve2(){    memset(f,0,sizeof(f));    for (int i=1;i<=m;i++) if (not_prime[i]) f[i%p]++;    memset(a.a,0,sizeof(a.a));    for (int j=1;j<=m;j++)        if (not_prime[j]) a.a[(-j%p+p)%p][0]++;    for (int i=1;i<p;i++)        for (int j=0;j<p;j++)            a.a[j][i]=a.a[(j-1+p)%p][i-1];    a=ksm(a,n-1);    int ans=0;    for (int i=0;i<p;i++) ans=(ans+(LL)f[i]*a.a[i][0]%MOD)%MOD;    return ans;}int main(){    scanf("%d%d%d",&n,&m,&p);    get_prime(m);    int x=solve1(),y=solve2();    printf("%d",(x-y+MOD)%MOD);    return 0;}