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