bzoj4818: [Sdoi2017]序列计数（矩阵快速幂优化dp）

78 阅读 0 评论 52 点赞

我是靠谱客的博主高挑钢笔，这篇文章主要介绍bzoj4818: [Sdoi2017]序列计数（矩阵快速幂优化dp），现在分享给大家，希望可以做个参考。

传送门
题意简述：问有多少长度为n的序列，序列中的数都是不超过m的正整数，而且这n个数的和是p的倍数，且其中至少有一个数是质数，答案对 $20170408$ 取模 $(n \leq 1 e 9, m \leq 2 e 7, p \leq 100)$ 。

思路：
首先因为只需要是 $p$ 的倍数，因此可以看成全局和对 $p$ 取模为 $0$ 方案数。
设状态 $f_{0/1,i,j}$ 表示不限制选出的数/选出的数不能是质数，允许选 $i$ 个数，全局和对 $p$ 取模等于 $j$ 方案数，答案就是 $f_{0,n,0}-f_{1,n,0}$ 。
用矩阵快速幂来优化这个 $d p$ 即可。
代码：

复制代码

#include<bits/stdc++.h>
#define ri register int
using namespace std;
typedef long long ll;
const int mod=20170408,M=2e7+5,N=105;
inline int add(const int&a,const int&b){return a+b>=mod?a+b-mod:a+b;}
inline int dec(const int&a,const int&b){return a>=b?a-b:a-b+mod;}
inline int mul(const int&a,const int&b){return (ll)a*b%mod;}
int n,m,p,s1[105],s2[N],pri[M],tot=0;
bool vis[M];
inline void init(){
	vis[1]=1;
	for(ri i=2;i<=m;++i){
		if(!vis[i])pri[++tot]=i;
		for(ri j=1;j<=tot&&i*pri[j]<=m;++j){
			vis[i*pri[j]]=1;
			if(i%pri[j]==0)break;
		}
	}
}
struct Mat{
	int a[N][N];
	Mat(int v=0){memset(a,0,sizeof(a));for(ri i=0;i<p;++i)a[i][i]=v;}
	friend inline Mat operator*(const Mat&a,const Mat&b){
		Mat ret;
		for(ri i=0;i<p;++i)for(ri k=0;k<p;++k)if(a.a[i][k])
			for(ri j=0;j<p;++j)ret.a[i][j]=add(ret.a[i][j],mul(a.a[i][k],b.a[k][j]));
		return ret;
	}
	friend inline Mat operator^(Mat&a,int p){
		Mat ret(1);
		for(;p;p>>=1,a=a*a)if(p&1)ret=ret*a;
		return ret;
	}
};
int main(){
	cin>>n>>m>>p;
	init();
	Mat f0,f1;
	for(ri i=1;i<=m;++i){
		++s1[i%p];
		if(vis[i])++s2[i%p];
	}
	f0.a[0][0]=f1.a[0][0]=1;
	for(ri i=0;i<p;++i)for(ri j=0;j<p;++j)f0.a[i][(i+j)%p]=s1[j]%mod,f1.a[i][(i+j)%p]=s2[j]%mod;
	f0=f0^n,f1=f1^n;
	cout<<dec(f0.a[0][0],f1.a[0][0]);
	return 0;
}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
#include<bits/stdc++.h>
#define ri register int
using namespace std;
typedef long long ll;
const int mod=20170408,M=2e7+5,N=105;
inline int add(const int&a,const int&b){return a+b>=mod?a+b-mod:a+b;}
inline int dec(const int&a,const int&b){return a>=b?a-b:a-b+mod;}
inline int mul(const int&a,const int&b){return (ll)a*b%mod;}
int n,m,p,s1[105],s2[N],pri[M],tot=0;
bool vis[M];
inline void init(){
	vis[1]=1;
	for(ri i=2;i<=m;++i){
		if(!vis[i])pri[++tot]=i;
		for(ri j=1;j<=tot&&i*pri[j]<=m;++j){
			vis[i*pri[j]]=1;
			if(i%pri[j]==0)break;
		}
	}
}
struct Mat{
	int a[N][N];
	Mat(int v=0){memset(a,0,sizeof(a));for(ri i=0;i<p;++i)a[i][i]=v;}
	friend inline Mat operator*(const Mat&a,const Mat&b){
		Mat ret;
		for(ri i=0;i<p;++i)for(ri k=0;k<p;++k)if(a.a[i][k])
			for(ri j=0;j<p;++j)ret.a[i][j]=add(ret.a[i][j],mul(a.a[i][k],b.a[k][j]));
		return ret;
	}
	friend inline Mat operator^(Mat&a,int p){
		Mat ret(1);
		for(;p;p>>=1,a=a*a)if(p&1)ret=ret*a;
		return ret;
	}
};
int main(){
	cin>>n>>m>>p;
	init();
	Mat f0,f1;
	for(ri i=1;i<=m;++i){
		++s1[i%p];
		if(vis[i])++s2[i%p];
	}
	f0.a[0][0]=f1.a[0][0]=1;
	for(ri i=0;i<p;++i)for(ri j=0;j<p;++j)f0.a[i][(i+j)%p]=s1[j]%mod,f1.a[i][(i+j)%p]=s2[j]%mod;
	f0=f0^n,f1=f1^n;
	cout<<dec(f0.a[0][0],f1.a[0][0]);
	return 0;
}