题目描述
Alice想要得到一个长度为nn的序列，序列中的数都是不超过mm的正整数，而且这nn个数的和是pp的倍数。
Alice还希望，这nn个数中，至少有一个数是质数。
Alice想知道，有多少个序列满足她的要求。
输入输出格式
输入格式：
一行三个数，n,m,pn,m,p。
输出格式：
一行一个数，满足Alice的要求的序列数量，答案对2017040820170408取模。
输入输出样例
说明
对20\%20%的数据，1leq n,mleq1001≤n,m≤100
对50\%50%的数据，1leq m leq 1001≤m≤100
对80\%80%的数据，1leq mleq 10^61≤m≤106
对100\%100%的数据，1leq n leq 10^9,1leq m leq 2 imes 10^7,1leq pleq 1001≤n≤109,1≤m≤2×107,1≤p≤100
时间限制：3s
空间限制：128MB
至少有一个素数的方案=所有方案-没有素数的方案
于是用容斥就变成了简单的dp，先讨论所有方案
令f[i][j]表示i个数，和%p为j的方案数
f[i][j]=∑f[i-1][(j-k+p)%p]*cnt[k]
cnt[k]是1~m中%p等于k的数量
发现显然上式可以写为矩阵
于是用矩阵快速幂就行
然后用欧拉筛把素数筛掉，再做一次
1 #include<iostream> 2 #include<cstdio> 3 #include<cstring> 4 #include<algorithm> 5 using namespace std; 6 typedef long long lol; 7 struct Matrix 8 { 9 lol a[201][201]; 10 }pre,ans1,ans2,Mat1,Mat2; 11 lol n,m,p,Mod=20170408; 12 long long cnt1[201],cnt2[201]; 13 lol tot,pri[25000001]; 14 bool vis[25000001]; 15 Matrix operator*(const Matrix a,const Matrix b) 16 { 17 lol i,j,k; 18 Matrix res; 19 memset(res.a,0,sizeof(res.a)); 20 for (k=0;k<p;k++) 21 for (i=0;i<p;i++) 22 if (a.a[i][k]) 23 { 24 for (j=0;j<p;j++) 25 { 26 res.a[i][j]+=a.a[i][k]*b.a[k][j]; 27 res.a[i][j]%=Mod; 28 } 29 } 30 return res; 31 } 32 Matrix qpow1(lol y) 33 {lol i; 34 Matrix res; 35 memset(res.a,0,sizeof(res.a)); 36 for (i=0;i<p;i++) 37 res.a[i][i]=1; 38 while (y) 39 { 40 if (y&1) res=res*Mat1; 41 Mat1=Mat1*Mat1; 42 y=y/2; 43 } 44 return res; 45 } 46 Matrix qpow2(lol y) 47 {lol i; 48 Matrix res; 49 memset(res.a,0,sizeof(res.a)); 50 for (i=0;i<p;i++) 51 res.a[i][i]=1; 52 while (y) 53 { 54 if (y&1) res=res*Mat2; 55 Mat2=Mat2*Mat2; 56 y=y/2; 57 } 58 return res; 59 } 60 int main() 61 {lol i,j; 62 cin>>n>>m>>p; 63 cnt1[1]++; 64 for (i=2;i<=m;i++) 65 { 66 cnt1[i%p]++,cnt1[i%p]%=Mod; 67 if (vis[i]==0) 68 { 69 ++tot; 70 pri[tot]=i; 71 } 72 for (j=1;j<=tot;j++) 73 { 74 if (i*pri[j]>m) break; 75 vis[i*pri[j]]=1; 76 if (i%pri[j]==0) break; 77 } 78 } 79 cnt2[1]++; 80 for(i=2;i<=m;i++) 81 if (vis[i]) cnt2[i%p]++,cnt2[i%p]%=Mod; 82 for (i=0;i<p;i++) 83 { 84 for (j=0;j<p;j++) 85 { 86 Mat1.a[i][(i+j)%p]+=cnt1[j]%Mod; 87 Mat1.a[i][(i+j)%p]%=Mod; 88 } 89 } 90 for (i=0;i<p;i++) 91 pre.a[0][i]=cnt1[i]; 92 ans1=qpow1(n-1); 93 ans1=pre*ans1; 94 for (i=0;i<p;i++) 95 { 96 for (j=0;j<p;j++) 97 { 98 Mat2.a[i][(i+j)%p]+=cnt2[j]%Mod; 99 Mat2.a[i][(i+j)%p]%=Mod; 100 } 101 } 102 for (i=0;i<p;i++) 103 pre.a[0][i]=cnt2[i]; 104 ans2=qpow2(n-1); 105 ans2=pre*ans2; 106 cout<<(ans1.a[0][0]-ans2.a[0][0]+Mod)%Mod; 107 }