【LOJ#2002】【SDOI2017】—序列计数（倍增dp）

传送门

$f[i][j][0/1]$

表示$i$个数，模$p$为$j$,有无质数的方案

倍增$dp$$O(m^2log)$就完了啊
为什么要$m^3log$的矩乘啊

#include<bits/stdc++.h>
using namespace std;
const int RLEN=(1<<20)|5;
inline char gc(){
	static char ibuf[RLEN],*ib,*ob;
	(ib==ob)&&(ob=(ib=ibuf)+fread(ibuf,1,RLEN,stdin));
	return (ib==ob)?EOF:*ib++;
}
#define gc getchar
inline int read(){
	char ch=gc();
	int res=0,f=1;
	while(!isdigit(ch))f^=ch=='-',ch=gc();
	while(isdigit(ch))res=(res+(res<<2)<<1)+(ch^48),ch=gc();
	return f?res:-res;
}
#define pb push_back
#define re register
#define fi first
#define se second
#define pii pair<int,int>
#define cs const
#define bg begin
#define ll long long
#define poly vector<int>
template<class tp>inline void chemx(tp &a,tp b){a<b?a=b:0;}
template<class tp>inline void chemn(tp &a,tp b){a>b?a=b:0;}
cs int mod=20170408;
inline int add(int a,int b){a+=b-mod;return a+(a>>31&mod);}
inline void Add(int &a,int b){a+=b-mod;a+=a>>31&mod;}
inline int dec(int a,int b){a-=b;return a+(a>>31&mod);}
inline void Dec(int &a,int b){a-=b;a+=a>>31&mod;}
inline int mul(int a,int b){static ll r;r=1ll*a*b;return r>=mod?r%mod:r;}
inline void Mul(int &a,int b){static ll r;r=1ll*a*b;a=r>=mod?r%mod:r;}
inline int ksm(int a,int b,int res=1){for(;b;b>>=1,Mul(a,a))(b&1)&&(Mul(res,a),1);return res;}
inline int Inv(int x){return ksm(x,mod-2);}
cs int N=20000005;
char xxx;
bitset<N> vis;
int pr[2000005],tot;
int n,m,p;
struct node{
	int a0,a1;
	node(int _a0=0,int _a1=0):a0(_a0),a1(_a1){}
	inline void operator +=(cs node &x){
		Add(a0,x.a0),Add(a1,x.a1);
	}
	friend inline node operator *(cs node &a,cs node &b){
		return node(mul(a.a0,b.a0),(1ll*a.a0*b.a1+1ll*b.a0*a.a1+1ll*a.a1*b.a1)%mod);
	}
}f[105],res[105],g[105];
inline void Mult(node *a,node *b,node *tmp){
	for(int i=0;i<p;i++)g[i]=node();
	for(int i=0;i<p;i++)
	for(int j=0;j<p;j++)
	g[(i+j)>=p?(i+j-p):(i+j)]+=a[i]*b[j];
	for(int i=0;i<p;i++)tmp[i]=g[i];
}
char yyy;
inline void init(int l){
	cs int len=l;
	vis[1]=1;
	for(int i=2;i<=len;i++){
		if(!vis[i])pr[++tot]=i;
		for(int j=1;j<=tot&&i*pr[j]<=len;j++){
			vis[i*pr[j]]=1;
			if(i-i/pr[j]*pr[j]==0)break;
		}
	}
}
int main(){
	n=read(),m=read(),p=read();
	init(m);
	for(int i=1;i<=m;i++)if(vis[i])f[i-i/p*p].a0++;else f[i-i/p*p].a1++;
	res[0].a0=1;
	for(;n;n>>=1,Mult(f,f,f))if(n&1)Mult(res,f,res);
	cout<<res[0].a1<<'
';
}