[SDOI2017]序列计数

首先将问题转换成一个容斥的形式，把质数的限制去掉。

显然有一个(Theta(NMP))的dp转移对吧。

我们不难发现每一次(dp)的转移都是一样的，所以直接矩阵快速幂优化即可。

/*
  mail: mleautomaton@foxmail.com
  author: MLEAutoMaton
  This Code is made by MLEAutoMaton
*/
#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<math.h>
#include<algorithm>
#include<queue>
#include<set>
#include<map>
#include<iostream>
using namespace std;
#define ll long long
#define re register
#define file(a) freopen(a".in","r",stdin);freopen(a".out","w",stdout)
inline int gi(){
	int f=1,sum=0;char ch=getchar();
	while(ch>'9' || ch<'0'){if(ch=='-')f=-1;ch=getchar();}
	while(ch>='0' && ch<='9'){sum=(sum<<3)+(sum<<1)+ch-'0';ch=getchar();}
	return f*sum;
}
const int N=20000010,P=110,Mod=20170408;
int n,m,p;
struct matrix{
	int a[P][P];
	void clear(){memset(a,0,sizeof(a));}
	int*operator[](int x){return a[x];}
	matrix operator*(matrix b){
		matrix c;
		for(int i=0;i<p;i++)
			for(int j=0;j<p;j++){
				c[i][j]=0;
				for(int k=0;k<p;k++)
					c[i][j]=(c[i][j]+1ll*a[i][k]*b[k][j]%Mod)%Mod;
			}
		return c;
	}
}S,T;
int ha[N];
bool is_prime[N];
int prime[5000010],s;
void init(){
	is_prime[1]=1;
	for(int i=2;i<=m;i++){
		if(!is_prime[i])prime[++s]=i;
		for(int j=1;j<=s && i*prime[j]<=m;j++){
			is_prime[i*prime[j]]=1;
			if(i%prime[j]==0)break;
		}
	}
}
int main(){
#ifndef ONLINE_JUDGE
	freopen("in.in","r",stdin);
#endif
	n=gi();m=gi();p=gi();init();
	for(int i=1;i<=m;i++)ha[i%p]++;
	for(int i=0;i<p;i++)S[0][i]=ha[i]%Mod;
	for(int i=0;i<p;i++)
		for(int j=0;j<p;j++)
			T[i][(i+j)%p]=(T[i][(i+j)%p]+ha[j])%Mod;
	int b=n-1;
	while(b){
		if(b&1)S=S*T;
		T=T*T;b>>=1;
	}
	int ans=S[0][0];
	memset(ha,0,sizeof(ha));T.clear();S.clear();
	for(int i=1;i<=m;i++)if(is_prime[i])ha[i%p]++;
	for(int i=0;i<=min(m,p-1);i++)S[0][i]=ha[i]%Mod;
	for(int i=0;i<p;i++)
		for(int j=0;j<p;j++)
			T[i][(i+j)%p]=(T[i][(i+j)%p]+ha[j])%Mod;
	b=n-1;
	while(b){
		if(b&1)S=S*T;
		T=T*T;b>>=1;
	}
	printf("%d
",(ans-S[0][0]+Mod)%Mod);
	return 0;
}