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()<<'
';
	return 0;
}