4818: [Sdoi2017]序列计数
Time Limit: 30 Sec Memory Limit: 128 MBSubmit: 754 Solved: 461
[Submit][Status][Discuss]
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
思路{
水题一道....首先容斥一下,$Ans=$所有的方案数-不含质数的方案数.
二维DP转移是很好想的 $ F [ i ] [ j ] $ 表示前$ i $个数和 $ %p $为$ j $的方案数.不含质数的同理
$ F [ i ] [ (j+x) ] + = F [ i - 1] [ j] $
这个暴力DP不对,看到n那么大,不含$Max,Min$等操作,直接矩阵快速幂了.
构造转移矩阵时可以第一行枚举m,求出来.
然后发现第二行以后都是上一行往左移动一位.
那么不用$ O( m*p) $的构造,直接$ O(m*m) $了.
在BZOJ上跑得贼慢......
}
#include<bits/stdc++.h>
#define il inline
#define RG register
#define ll long long
#define db double
#define N 110
#define mod 20170408
using namespace std;
bool vis[20000000];int P[2000000];
void pre(){
vis[1]=true;
for(int i=2;i<=20000000;++i){
if(!vis[i])P[++P[0]]=i;
for(int j=1;j<=P[0]&&P[j]*i<=20000000;++j){
vis[i*P[j]]=true;
if(!(i%P[j]))break;
}
}
}
int ss[N];
struct matrix{
int n,m;
ll ma[N][N];
matrix operator *(const matrix & a)const{
matrix c;memset(c.ma,0,sizeof(c.ma));
if(m!=a.n)return c;c.n=n,c.m=a.m;
for(RG int i=1;i<=c.n;++i)
for(RG int j=1;j<=c.m;++j){
for(RG int h=1;h<=a.n;++h){
c.ma[i][j]+=(ma[i][h]*a.ma[h][j])%mod;
if(c.ma[i][j]>=mod)c.ma[i][j]-=mod;
}
}
return c;
}
}ans,temp;
matrix qp(matrix a,ll b){
matrix Ans;Ans.n=Ans.m=a.n;
for(int i=1;i<=a.n;++i)
for(int j=1;j<=a.n;++j)
Ans.ma[i][j]=(i==j);
for(;b;b>>=1,a=a*a)
if(b&1)Ans=Ans*a;
return Ans;
}
int n,m,p;
int main(){
scanf("%d%d%d",&n,&m,&p);pre();
ans.n=p,ans.m=1;
for(int i=1;i<=m;++i){
ans.ma[(i%p)+1][1]++;
if(ans.ma[(i%p)+1][1]>=mod)ans.ma[(i%p)+1][1]-=mod;
}
for(int i=1;i<=m;++i){
int tmp=i%p;tmp=p-tmp;tmp++;
if(tmp>p)tmp=1;
temp.ma[1][tmp]++;
if(temp.ma[1][tmp]>=mod)temp.ma[1][tmp]-=mod;
}
for(int i=2;i<=p;++i){
for(int j=2;j<=p;++j)
temp.ma[i][j]=temp.ma[i-1][j-1];
temp.ma[i][1]=temp.ma[i-1][p];
}temp.n=temp.m=p;
temp=qp(temp,n-1);
ans=temp*ans;
ll Ans1=ans.ma[1][1];
memset(ans.ma,0,sizeof(ans.ma));
for(int i=1;i<=m;++i)
if(vis[i]){
ans.ma[(i%p)+1][1]++;
if(ans.ma[(i%p)+1][1]>=mod)ans.ma[(i%p)+1][1]-=mod;
}
memset(temp.ma,0,sizeof(temp.ma));
for(int i=1;i<=m;++i){
if(!vis[i])continue;
int tmp=i%p;tmp=p-tmp;tmp++;
if(tmp>p)tmp=1;
temp.ma[1][tmp]++;
if(temp.ma[1][tmp]>=mod)temp.ma[1][tmp]-=mod;
}
for(int i=2;i<=p;++i){
for(int j=2;j<=p;++j)
temp.ma[i][j]=temp.ma[i-1][j-1];
temp.ma[i][1]=temp.ma[i-1][p];
}temp.n=temp.m=p;
ans.n=p,ans.m=1;
temp=qp(temp,n-1);
ans=temp*ans;
ll Ans2=ans.ma[1][1];
cout<<(Ans1-Ans2+mod)%mod;
return 0;
}