【Bzoj4818】【Sdoi2017】序列计数

首先，想到容斥原理
答案是 所有情况 - 没有一个素数的
两部分分别处理
dp方程非常好列
可以用矩阵快速幂优化

P.S. 矩阵这一块需要提高啊！

#include<cstdio>
#include<algorithm>
#define ll long long
using namespace std;
const ll mod=20170408;
const ll maxm=20000007;
ll n,m,p,v[110],g[110][110],now[110][110],ret[110][110];
bool is_prime[maxm];
void get_prime(){
	fill(is_prime,is_prime+m+1,true);
	is_prime[0]=is_prime[1]=false;
	//is_prime[1]=true;
	for (ll i=2; i<=m; i++) if (is_prime[i]) {
		for (ll j=2*i; j<=m; j+=i) is_prime[j]=false;
	}
}
void mult(){
//ret and g
	for (ll i=0; i<p; i++)
		for (ll j=0; j<p; j++) {
			now[i][j]=0;
			for (ll k=0; k<p; k++)
				now[i][j]+=(ret[i][k]*g[k][j])%mod,now[i][j]%=mod;
		}
	for (ll i=0; i<p; i++) 
		for (ll j=0; j<p; j++) ret[i][j]=now[i][j];
}
void calc(){
	for (ll i=0; i<p; i++)
		for (ll j=0; j<p; j++) {
			now[i][j]=0;
			for (ll k=0; k<p; k++)
				now[i][j]+=(g[i][k]*g[k][j])%mod,now[i][j]%=mod;
		}
	for (ll i=0; i<p; i++) 
		for (ll j=0; j<p; j++) g[i][j]=now[i][j];
}
void qpow(ll x){
	for (ll i=0; i<p; i++) 
		for (ll j=0; j<p; j++) ret[i][j]=0;
	for (ll i=0; i<p; i++) ret[i][i]=1;
	while (x){
		if (x & 1) mult();
		calc();
		x >>= 1;
	}
}
int main(){
	scanf("%lld%lld%lld",&n,&m,&p); 
	get_prime();
	for (ll i=1; i<=m; i++) if (!is_prime[i]) ++v[i % p];
	for (ll i=0; i<p; i++) 
		for (ll j=0; j<p; j++) g[i][j]=v[(i-j+p) % p]; 
	qpow(n); 
	ll ans2=ret[0][0];
	fill(v,v+p+1,0);
	for (ll i=1; i<=m; i++) ++v[i % p];
	for (ll i=0; i<p; i++) 
		for (ll j=0; j<p; j++) 
		g[i][j]=v[(i-j+p) % p];
	qpow(n);
	ll ans1=ret[0][0];
	printf("%lld\n",(ans1-ans2+mod)%mod);
	return 0;
}