luogu3702-[SDOI2017]序列计数
Description
Alice想要得到一个长度为nn的序列,序列中的数都是不超过mm的正整数,而且这nn个数的和是pp的倍数。
Alice还希望,这nn个数中,至少有一个数是质数。
Alice想知道,有多少个序列满足她的要求。
输入输出格式
输入格式:
一行三个数,n,m,p。
输出格式:
一行一个数,满足Alice的要求的序列数量,答案对20170408取模。
输入输出样例
输入样例#1:
3 5 3
输出样例#1:
33
说明
对20%的数据,1≤n,m≤100
对50%的数据,1≤m≤100
对80%的数据,1≤m≤10^6
对100%的数据,1≤n≤10^9 ,1≤m≤2×10^7,1≤p≤100
时间限制:3s
空间限制:128MB
Solution
洛谷P3702 [SDOI2017]序列计数 - 蒟蒻的博客 - CSDN博客
母函数.
设 \(f_i\) 为多项式f的 \(x^i\) 项系数.
设 \(f_i\) 表示\(a\)个数, 和 mod p为 \(i\) 方案数.
多项式快速幂即可.
Code
#include<cstdio>
#include<iostream>
#include<cmath>
#include<cstring>
#include<algorithm>
#include<set>
#include<map>
using namespace std;
#define rep(i,l,r) for(register int i=(l);i<=(r);++i)
#define repdo(i,l,r) for(register int i=(l);i>=(r);--i)
#define il inline
typedef double db;
typedef long long ll;
//---------------------------------------
const int psz=150,msz=2e7+50,nmod=20170408;
int n,m,p;
int nopr[msz],pr[msz],pp=0;
void init(int bnd){
nopr[1]=1;//a
rep(i,2,bnd){
if(nopr[i]==0)pr[++pp]=i;//b
rep(j,1,pp){
if(i*pr[j]>bnd)break;
nopr[i*pr[j]]=1;
if(i%pr[j]==0)break;
}
}
}
void add(ll &a,ll b){a=(a+b)%nmod;}
struct tpo{
ll val[psz];
void cl(){memset(val,0,sizeof(val));}
tpo(){cl();}
ll& operator[](int p){return val[p];}
const ll& operator[](int p)const{return val[p];}
}f1,f2;
typedef const tpo& ctpo;
tpo operator*(ctpo a,ctpo b){
tpo res;
rep(i,0,p-1)rep(j,0,p-1)add(res[(i+j)%p],a[i]*b[j]%nmod);
return res;
}
tpo operator^(tpo a,ll p){
tpo res;res[0]=1;
for(;p;a=a*a,p>>=1){
if(p&1)res=res*a;
}
return res;
}
int sol(){
init(2e7+5);
rep(i,1,m){
++f1[i%p];
if(nopr[i])++f2[i%p];
}
return ((((f1^n)[0]-(f2^n)[0])%nmod)+nmod)%nmod;
}
int main(){
ios::sync_with_stdio(0),cin.tie(0);
cin>>n>>m>>p;
cout<<sol()<<'\n';
return 0;
}