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;
}

原文地址:https://www.cnblogs.com/Mychael/p/8781495.html