[BZOJ4818][矩阵快速幂优化DP]SDOI2017：序列计数

BZOJ4818

状态比较好想：$f[i][j]$表示第$i$个数，它们的和$\%p$ 为 $j$的方案数
转移方程$f[i][(j+k)\%p]=\sum f[i-1][j]\ast s[k]$，s[k]表示mod p余k的数的个数
用矩阵快速幂优化转移即可
关于质数的限定：我们先求总方案数，然后求一个质数都不选的方案数，相减就完了

Code：

#include<bits/stdc++.h>
#define mod 20170408
#define ll long long
using namespace std;
inline int read(){
	int res=0,f=1;char ch=getchar();
	while(!isdigit(ch)) {if(ch=='-') f=-f;ch=getchar();}
	while(isdigit(ch)) {res=(res<<1)+(res<<3)+(ch^48);ch=getchar();}
	return res*f;
}
const int N=105,M=2e7+5;
int n,m,p;
int pt[M],pri[M],tot=0;
ll s[N];
struct Ma{
	ll a[N][N];
	Ma(){memset(a,0,sizeof(a));}
	Ma operator * (const Ma b)const{
		Ma res;
		for(int i=0;i<p;i++)
			for(int j=0;j<p;j++)
				for(int k=0;k<p;k++)
					res.a[i][j]=(res.a[i][j]+a[i][k]*b.a[k][j]%mod)%mod;
		return res;
	}
};
inline void init(){
	pt[1]=1;
	for(int i=2;i<=M-5;i++){
		if(!pt[i]) pri[++tot]=i;
		for(int j=1;j<=tot && i*pri[j]<=M-5;j++){
			pt[i*pri[j]]=1;
			if(i%pri[j]==0) break;
		}
	}
}
inline Ma ksm(Ma a,int b){
	Ma ans;
	for(int i=0;i<p;i++) ans.a[i][i]=1;
	for(;b;b>>=1){
		if(b&1) ans=ans*a;
		a=a*a;
	}
	return ans;
}
inline int solve(){
	Ma ret;
	for(int i=0;i<p;i++)
		for(int j=0;j<p;j++) ret.a[i][(i+j)%p]=s[j];
	ret=ksm(ret,n);
	return ret.a[0][0];
}
int main(){
	init();
	n=read();m=read();p=read();
	for(int i=1;i<=m;i++) ++s[i%p];ll ans=solve();
	for(int i=1;i<=m;i++) if(!pt[i]) --s[i%p];ans-=solve();
	cout<<(ans+mod)%mod;
	return 0;
}