注意到
p
很小,容易想到矩乘。
关于那个质数的限制,只需要
1
到
构造
1
到
扣去质数直接暴力
O(m/lnm∗P)
。
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long LL;
const int maxn=20000005,MOD=20170408;
int ans,n,m,P,p[maxn],vis[maxn];
void getP(){
for(int i=2;i<=m;i++){
if(!vis[i]) p[++p[0]]=i;
for(int j=1;j<=p[0]&&(LL)i*p[j]<=m;j++){
vis[i*p[j]]=true;
if(i%p[j]==0) break;
}
}
}
struct Matrix{
int n,m,a[105][105];
Matrix(){ memset(a,0,sizeof(a)); n=m=0; }
Matrix operator * (const Matrix &b){
Matrix c; c.n=n; c.m=b.m;
for(int i=0;i<c.n;i++)
for(int j=0;j<c.m;j++)
for(int k=0;k<m;k++) (c.a[i][j]+=(LL)a[i][k]*b.a[k][j]%MOD)%=MOD;
return c;
}
} T,Ans;
Matrix Pow(Matrix a,int b){
Matrix res; res.n=res.m=P; for(int i=0;i<P;i++) res.a[i][i]=1;
for(;b;b>>=1,a=a*a) if(b&1) res=res*a;
return res;
}
int main(){
freopen("loj2002.in","r",stdin);
freopen("loj2002.out","w",stdout);
scanf("%d%d%d",&n,&m,&P);
getP();
T=Matrix(); T.n=T.m=P;
for(int i=0;i<=P-1&&i<=m;i++)
for(int j=0;j<=P-1;j++) (T.a[j][(j+i)%P]+=(m-i)/P+(i>0))%=MOD;
Ans=Matrix(); Ans.n=1; Ans.m=P; Ans.a[0][0]=1;
Ans=Ans*Pow(T,n); ans=Ans.a[0][0];
for(int i=1;i<=p[0];i++)
for(int j=0;j<=P-1;j++) T.a[j][(j+p[i])%P]--;
Ans=Matrix(); Ans.n=1; Ans.m=P; Ans.a[0][0]=1;
Ans=Ans*Pow(T,n); ans=(ans-Ans.a[0][0])%MOD;
printf("%d\n",(ans+MOD)%MOD);
return 0;
}