BZOJ4818 [Sdoi2017]序列计数

看数据范围一眼矩乘

f[i][j]表示前i个数，和模p余j

用任意取的方案减去不取质数的方案即得到至少有一个质数的方案

对任取和不取质数分别随便做个转移矩阵快速幂即可

#include<iostream>
#include<cstring>
#include<ctime>
#include<cmath>
#include<algorithm>
#include<iomanip>
#include<cstdlib>
#include<cstdio>
#include<map>
#include<bitset>
#include<set>
#include<stack>
#include<vector>
#include<queue>
using namespace std;
#define MAXN 100
#define MAXM 20000010
#define ll long long
#define eps 1e-8
#define MOD 20170408
#define INF 1000000000
int P;
struct mat{
	ll x[MAXN];
	mat(){
		memset(x,0,sizeof(x));
	}
	mat operator *(const mat &y){
		mat z;
		int i,j;
		for(i=0;i<P;i++){
			for(j=0;j<P;j++){
				z.x[i]+=x[j]*y.x[(i-j+P)%P];
				if(z.x[i]>1e18){
					z.x[i]%=MOD;
				}
			}
			z.x[i]%=MOD;
		}
		return z;
	}
};
int n,m;
int p[MAXM],tot;
bool np[MAXM];
mat Xa,Xn,Ra,Rn;
int ca[MAXN],cn[MAXN];
int ans;
void su(){
	int i,j;
	np[1]=1;
	for(i=2;i<MAXM;i++){
		if(!np[i]){
			p[++tot]=i;
		}
		for(j=1;(ll)i*p[j]<MAXM;j++){
			np[i*p[j]]=1;
			if(!(i%p[j])){
				break;
			}
		}
	}
}
int main(){
	int i,j;
	su();
	scanf("%d%d%d",&n,&m,&P);
	for(i=1;i<=m;i++){
		Xa.x[i%P]++;
		if(np[i]){
			Xn.x[i%P]++;
		}
	}
	Ra.x[0]=Rn.x[0]=1;
	while(n){
		if(n&1){
			Ra=Ra*Xa;
			Rn=Rn*Xn;
		}
		Xa=Xa*Xa;
		Xn=Xn*Xn;
		n>>=1;
	}
	printf("%d\n",(Ra.x[0]-Rn.x[0]+MOD)%MOD);
	return 0;
}

/*

*/