[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\n",(ans-S[0][0]+Mod)%Mod);
return 0;
}