[BZOJ]4818: [Sdoi2017]序列计数

Time Limit: 30 Sec Memory Limit: 128 MB

Description

　　Alice想要得到一个长度为n的序列，序列中的数都是不超过m的正整数，而且这n个数的和是p的倍数。Alice还希望，这n个数中，至少有一个数是质数。Alice想知道，有多少个序列满足她的要求。

Input

　　一行三个数，n,m,p。
　　1<=n<=10^9,1<=m<=2×10^7,1<=p<=100

Output

　　一行一个数，满足Alice的要求的序列数量，答案对20170408取模。

Sample Input

　　3 5 3

Sample Output

　　33

Solution

　　f1[i][j]表示任意选i个1~m内的数，和模p为j的方案数，f2[i][j]表示选i个1~m内的不是质数的数，和模p为j的方案数，f1
[0]-f2
[0]就是答案，矩阵乘法加速DP就可以了，复杂度$O(m+p^3logn)$。

Code

#include<cstdio>
#include<cstring>
#define MN 20000000
#define MP 1270610
#define MK 100
#define MOD 20170408
struct mat
{
int z[MK][MK];
mat(){memset(z,0,sizeof(z));}
mat operator*(const mat&b)
{
mat c;int i,j,k;
for(i=0;i<MK;++i)for(k=0;k<MK;++k)if(z[i][k])for(j=0;j<MK;++j)
c.z[i][j]=(c.z[i][j]+1LL*z[i][k]*b.z[k][j])%MOD;
return c;
}
}a,x;
bool u[MN+5];
int p[MP+5],pn,s[MK];
int main()
{
int n,m,k,i,j,ans;
for(u[1]=1,i=2;i<=MN;++i)
{
if(!u[i])p[++pn]=i;
for(j=1;i*p[j]<=MN;++j){u[i*p[j]]=1;if(i%p[j]==0)break;}
}
scanf("%d%d%d",&n,&m,&k);
for(i=1;i<=m;++i)++s[i%k];
a.z[0][0]=1;
for(i=0;i<k;++i)for(j=0;j<k;++j)x.z[i][(i+j)%k]=s[j]%MOD;
for(i=n;i;i>>=1,x=x*x)if(i&1)a=a*x;
ans=a.z[0][0];
memset(s,0,sizeof(s));memset(a.z,0,sizeof(a.z));
for(i=1;i<=m;++i)if(u[i])++s[i%k];
a.z[0][0]=1;
for(i=0;i<k;++i)for(j=0;j<k;++j)x.z[i][(i+j)%k]=s[j]%MOD;
for(i=n;i;i>>=1,x=x*x)if(i&1)a=a*x;
printf("%d",(ans-a.z[0][0]+MOD)%MOD);
}