思路：
矩阵快速幂；
代码：
#include <bits/stdc++.h> using namespace std; #define mod 20170408 #define ll long long struct MatrixType { int n,m; ll ai[105][105]; void mem(int n_,int m_) { n=n_,m=m_; for(int i=0;i<=n;i++) for(int v=0;v<=m;v++) ai[i][v]=0; } MatrixType operator*(const MatrixType pos)const { MatrixType res;res.mem(this->n,pos.m); for(int i=0;i<=res.n;i++) for(int v=0;v<=res.m;v++) for(int x=0;x<=this->m;x++) res.ai[i][v]=(res.ai[i][v]+(this->ai[i][x]*pos.ai[x][v])%mod)%mod; return res; } void debug() { puts(""); printf("%d %d ",n,m); for(int i=0;i<=2;i++) { for(int v=0;v<=2;v++) printf("%d ",this->ai[i][v]); putchar(' '); } puts(""); } }; int n,m,p,pi[2000005],cnt,ai[1001],bi[1001]; bool if_p[20000005]; MatrixType poww(MatrixType pos,int mi) { MatrixType res=pos;mi--; while(mi) { if(mi&1) res=res*pos; mi=mi>>1,pos=pos*pos; } return res; } void ouler(int lit) { for(int i=2;i<=lit;i++) { if(!if_p[i]) pi[++cnt]=i; for(int j=1;pi[j]*i<=lit&&j<=cnt;j++) { if_p[i*pi[j]]=true; if(i%pi[j]==0) break; } } } int main() { //freopen("ans3.txt","w",stdout); scanf("%d%d%d",&n,&m,&p),ouler(m); MatrixType sta1;sta1.mem(0,p-1); sta1.ai[0][0]=1; MatrixType pos1,pos2,ans1,ans2; pos1.mem(p-1,p-1),pos2.mem(p-1,p-1); for(int i=1;i<=m;i++) ai[i%p]++; for(int i=0;i<p;i++) bi[i]=ai[i]; for(int i=1;i<=cnt;i++) bi[pi[i]%p]--; for(int i=0;i<p;i++) { for(int v=0;v<p;v++) { if(i>v) pos1.ai[i][v]=ai[i-v],pos2.ai[i][v]=bi[i-v]; if(i==v) pos1.ai[i][v]=ai[0],pos2.ai[i][v]=bi[0]; if(i<v) pos1.ai[i][v]=ai[p-v+i],pos2.ai[i][v]=bi[p-v+i]; } } pos1=poww(pos1,n),pos2=poww(pos2,n); ans1=sta1*pos1,ans2=sta1*pos2; cout<<(ans1.ai[0][0]-ans2.ai[0][0]+mod)%mod; return 0; }