【BZOJ4818】[Sdoi2017]序列计数
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
题解：至少包含1个质数的数量=总数-不包含质数的数量 （这种补集法不是一次两次见到了吧？）
于是我们考虑用DP求解，先快筛出1..m内的质数，1..m内除以P模为j的数的个数，1..m内除以P模为j的合数的个数
然后设f[i][j]表示i个数，总和除以P模j的方案数，g[i][j]表示i个合数，总和除以P模j的方案数，容易得出
f[i+1][(j+k)%P]+=f[i][j]+1..m内除以P模为j的数的个数
g[i+1][(j+k)%P]+=g[i][j]+1..m内除以P模为j的合数的个数
发现时间复杂度O(np)，用矩乘快速幂优化一下就好啦
#include <cstdio>
#include <cstring>
#include <iostream>
#define mod 20170408
using namespace std;
typedef long long ll;
int np[20000010],cnt[110],sum[110],pri[10000010];
int n,m,p,tot;
typedef struct matrix
{
ll v[110][110];
}M;
M x,ans,emp;
ll ans1;
M mmul(M a,M b)
{
M c=emp;
int i,j,k;
for(i=0;i<p;i++)
for(j=0;j<p;j++)
for(k=0;k<p;k++)
c.v[i][j]=(c.v[i][j]+a.v[i][k]*b.v[k][j])%mod;
return c;
}
void pm(int y)
{
while(y)
{
if(y&1) ans=mmul(ans,x);
x=mmul(x,x),y>>=1;
}
}
int main()
{
scanf("%d%d%d",&n,&m,&p);
int i,j;
np[1]=cnt[1]=sum[1]=1;
for(i=2;i<=m;i++)
{
sum[i%p]=(sum[i%p]+1)%mod;
if(!np[i]) pri[++tot]=i;
else cnt[i%p]=(cnt[i%p]+1)%mod;
for(j=1;j<=tot&&i*pri[j]<=m;j++)
{
np[i*pri[j]]=1;
if(i%pri[j]==0) break;
}
}
for(i=0;i<p;i++)
for(j=0;j<p;j++)
x.v[i][(i+j)%p]=(x.v[i][(i+j)%p]+sum[j])%mod;
ans.v[0][0]=1;
pm(n);
ans1=ans.v[0][0];
memset(ans.v,0,sizeof(ans.v)),memset(x.v,0,sizeof(x.v));
ans.v[0][0]=1;
for(i=0;i<p;i++)
for(j=0;j<p;j++)
x.v[i][(i+j)%p]=(x.v[i][(i+j)%p]+cnt[j])%mod;
pm(n);
printf("%lld",(ans1-ans.v[0][0]+mod)%mod);
return 0;
}