SDOI2017 序列计数

最新推荐文章于 2020-04-22 15:53:09 发布

原创最新推荐文章于 2020-04-22 15:53:09 发布 · 213 阅读

0 ·

CC 4.0 BY-SA版权

本文介绍了一种使用动态规划解决特定数学问题的方法，重点在于计算包含或不包含质数的数列组合数量。通过矩阵乘法进行状态转移，算法复杂度为O(p^3*logn)，并结合了O(m)的线筛算法来预处理质数信息。

传送门

$d p [i] [j] [0 / 1]$ 表示前 $i$ 个数的和模 $p$ 余数为 $j$ ，没有质数/有质数。
$c n t [i]$ 表示 $1$ ~ $m$ 中模 $p$ 余数为 $i$ 的数的个数。
$p c n t [i]$ 表示 $1$ ~ $m$ 中模 $p$ 余数为 $i$ 的质数的个数。

下面的数组下标减法为模 $p$ 意义下的。
那么 $d p [i] [j] [0]$ 必须用非质数转移：
$dp[i][j][0]=\sum_{t=0}^{p-1} dp[i-1][t][0]*(cnt[j-t]-pcnt[j-t])$
$d p [i] [j] [1]$ 由前面的两个状态转移：
$dp[i][j][1]=\sum_{t=0}^{p-1} dp[i-1][t][0]*pcnt[j-t]+dp[i-1][t][1]*cnt[j-t]$

构造矩阵转移即可（矩阵分块）。复杂度约为 $O(p^3*logn)$ ，还要加上一个 $O (m)$ 的线筛。

#include<bits/stdc++.h>
#define re register
#define cs const

cs int M=2e7+10,P=205,mod=20170408;

int n,m,p,lp,cnt[P],pcnt[P];

inline int add(int x,int y){return x+y>=mod?x+y-mod:x+y;}
inline int mul(int x,int y){return 1ll*x*y%mod;}

inline void Add(int &x,int y){x=(x+y>=mod)?x+y-mod:x+y;}
inline int pdec(int x,int y){return x-y<0?x-y+p:x-y;}


struct matrix{
	int a[P][P];
	matrix(int t=0){
		memset(a,0,sizeof a);
		for(int re i=0;i<P;++i) a[i][i]=t;
	}
	friend inline matrix operator*(cs matrix &a,cs matrix &b){
		matrix c;
		for(int re i=0;i<lp;++i)
			for(int re k=0;k<lp;++k)
				for(int re j=0;j<lp;++j)
					Add(c.a[i][j],mul(a.a[i][k],b.a[k][j]));
		return c;
	}
	friend inline matrix operator^(matrix a,int b){
		matrix ret(1);
		for(;b;b>>=1,a=a*a) if(b&1) ret=ret*a;
		return ret;
	}
}trans,A;

int mark[M],Pri[M],tot=0;
inline void linear_sieves(){
	mark[0]=mark[1]=1;
	for(int re i=2;i<=m;++i){
		if(!mark[i]) Pri[++tot]=i,++pcnt[i%p];
		for(int re j=1;j<=tot&&i*Pri[j]<=m;++j){
			mark[i*Pri[j]]=1;
			if(i%Pri[j]==0) break;
		}
	}
}

inline void init(){
	for(int re i=0,k;i<p;++i){
		for(int re j=0;j<p;++j){
			k=pdec(i,j);
			trans.a[i][j]=cnt[k]-pcnt[k];
			trans.a[i][j+p]=pcnt[k];
			trans.a[i+p][j+p]=cnt[k];
		}
	}A.a[0][0]=1;
}

int main(){
//	freopen("2302.in","r",stdin);
	scanf("%d%d%d",&n,&m,&p);
	lp=p<<1,linear_sieves();
	for(int re i=0;i<p;++i) cnt[i]=m/p;
	for(int re i=1;i<=(m%p);++i) ++cnt[i];
	init(),printf("%d\n",(A*(trans^n)).a[0][p]);
}