SDOI2017 序列计数

传送门

$dp[i][j][0/1]$表示前$i$个数的和模$p$余数为$j$，没有质数/有质数。
$cnt[i]$表示$1$~$m$中模$p$余数为$i$的数的个数。
$pcnt[i]$表示$1$~$m$中模$p$余数为$i$的质数的个数。

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

构造矩阵转移即可（矩阵分块）。复杂度约为$O(p^3\ast 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]);
}