BZOJ4818 [SDOI2017]序列计数 【生成函数 + 快速幂】

题目

Alice想要得到一个长度为n的序列，序列中的数都是不超过m的正整数，而且这n个数的和是p的倍数。Alice还希望
，这n个数中，至少有一个数是质数。Alice想知道，有多少个序列满足她的要求。

输入格式

一行三个数，n,m,p。
1<=n<=10^{9,1<=m<=2×10}7,1<=p<=100

输出格式

一行一个数，满足Alice的要求的序列数量，答案对20170408取模。

输入样例

3 5 3

输出样例

33

题解

由题目中“至少一个是质数的条件”容易想到翻转一下，用总方案减去全是合数的方案

\(n\)个数凑出\(p\)的倍数，容易想到生成函数
我们构造一个生成函数\(G(x)\)，第\(i\)次项的系数表示凑出模\(p\)意义下得\(i\)的数有多少种方案
显然我们枚举第一个数算出第一个\(G(x)\)的系数，只需要求出\(G(x)^n\)，对应\(0\)项系数就是总方案
全是合数的方案，可以再构造一个生成函数\(F(x)\)，只需要第一次枚举合数即可

甚至不用上fft，直接乘即可
\(O(m + p^2logn)\)

#include<iostream>
#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
#include<bitset>
#define LL long long int
#define Redge(u) for (int k = h[u],to; k; k = ed[k].nxt)
#define REP(i,n) for (int i = 1; i <= (n); i++)
#define BUG(s,n) for (int i = 1; i <= (n); i++) cout<<s[i]<<' '; puts("");
using namespace std;
const int maxn = 20000005,maxm = 100005,INF = 1000000000,P = 20170408;
bitset<maxn> isn;
int pr[maxn >> 1],pi,n,m,p;
void init(){
	isn[1] = true;
	for (int i = 2; i <= m; i++){
		if (!isn[i]) pr[++pi] = i;
		for (int j = 1; j <= pi && i * pr[j] <= m; j++){
			isn[i * pr[j]] = true;
			if (i % pr[j] == 0) break;
		}
	}
}
struct node{
	LL a[105];
	node(){memset(a,0,sizeof(a));}
}F,G;
inline node operator *(const node& a,const node& b){
	node c;
	for (int i = 0; i < p; i++)
		for (int j = 0; j < p; j++){
			int t = (i + j) % p;
			c.a[t] = (c.a[t] + a.a[i] * b.a[j] % P) % P;
		}
	return c;
}
inline node operator ^(node a,LL b){
	node ans;
	ans.a[0] = 1;
	for (; b; b >>= 1,a = a * a)
		if (b & 1) ans = ans * a;
	return ans;
}
LL qpow(LL a,LL b){
	LL ans = 1;
	for (; b; b >>= 1,a = a * a % P)
		if (b & 1) ans = ans * a % P;
	return ans;
}
int main(){
	scanf("%d%d%d",&n,&m,&p);
	init();
	for (int i = 1; i <= m; i++){
		F.a[i % p]++;
		if (isn[i]) G.a[i % p]++;
	}
	node A = F^n,B = G^n;
	LL ans;
	ans = ((A.a[0] - B.a[0]) % P + P) % P;
	cout << ans << endl;
	return 0;
}