序列计数

/*
    先求出全部的再删去没有质数的，矩阵乘法随便转移一下。
    f[i][j]表示选i个，总和%p=j的方案数。 
*/
#include<iostream>
#include<cstdio>
#include<cstring>
#define N 110
#define M 20000010
#define lon long long
#define mod 20170408
using namespace std;
int n,m,p,mark[M],prime[M],num,s1[N],s2[N];
struct node{
    node(){memset(a,0,sizeof(a));}
    lon a[N][N];
};node jc;
node operator *(node x,node y){
    node z;
    for(int i=0;i<p;i++)
        for(int j=0;j<p;j++)
            for(int k=0;k<p;k++)
                z.a[i][j]=(z.a[i][j]+x.a[i][k]*y.a[k][j])%mod;
                //注意这里只mod一次就好，mod多了会超时！ 
    return z;
}
node operator ^(node x,int y){
    node z;
    for(int i=0;i<p;i++) z.a[i][i]=1;
    while(y){
        if(y&1) z=z*x;
        x=x*x;
        y>>=1;
    }
    return z;
}
void init(){
    for(int i=2;i<=m;i++) mark[i]=1;
    for(int i=2;i<=m;i++){
        if(mark[i]) prime[++num]=i;
        for(int j=1;j<=num&&i*prime[j]<=m;j++){
            mark[i*prime[j]]=0;
            if(i%prime[j]==0) break;
        }
    }
    int mm=m/p;
    for(int i=1;i<=mm+1;i++)
       for(int j=1;j<=p&&(i-1)*p+j<=m;j++){
           if(!mark[(i-1)*p+j]) s2[j]++;//!pri 
           s1[j]++;//all
       }
    s1[0]=s1[p];s2[0]=s2[p];
}
int main(){
    scanf("%d%d%d",&n,&m,&p);
    init();
    for(int i=0;i<p;i++)
        for(int j=0;j<p;j++)
            jc.a[i][j]=s1[(i-j+p)%p];
    jc=jc^n;
    lon ans=jc.a[0][0];
    for(int i=0;i<p;i++)
        for(int j=0;j<p;j++)
            jc.a[i][j]=s2[(i-j+p)%p];
    jc=jc^n;
    cout<<((ans-jc.a[0][0])%mod+mod)%mod;
    return 0;
}