思路

先考慮暴力\(dp\)，\(f[i][j]\)表示前\(i\)個數，數字之和模\(P\)餘\(j\)的方案數。
我們先不考慮必須有質數這個條件，先統計出全部方案。然後再減去沒有質數的方案就行了。
那麼就有\(f[i + 1][(j + k) \% p] += f[i][j](1\le k \le m)\)
然後發現這個其實並不需要\(O(m)\)的轉移，因為\((j + k) \% p\)對於每個餘數來說，最少有\(\lfloor\frac{m}{p}\rfloor\)個。然後有\(m\%p\)個數有\(\lfloor\frac{m}{p}\rfloor + 1\)個.只要對這\(m\%p\)個數單獨處理一下就行了。
然後計算沒有質數的方案數，全部用合數轉移即可。
然後發現這個是可以矩陣優化的。
轉移矩陣的第\(i\)行第\(j\)列表示從\(1 + i\)到\(m+i\)中有多少個數字模\(p\)餘\(j\)。同樣可以計算一下，然後\(O(p^2)\)的列出矩陣。
對於全部用合數轉移的情況，轉移矩陣可以\(O(p\frac{m}{lnm})\)轉移,親測吸氧可過

代碼

/*
* @Author: wxyww
* @Date:   2019-02-13 16:22:09
* @Last Modified time: 2019-02-13 20:51:05
*/
#include<cstdio>
#include<iostream>
#include<cstdlib>
#include<cstring>
#include<algorithm>
#include<queue>
#include<vector>
#include<ctime>
using namespace std;
typedef long long ll;
// #define int ll
const int mod = 20170408;
ll read() {
    ll x=0,f=1;char c=getchar();
    while(c<'0'||c>'9') {
        if(c=='-') f=-1;
        c=getchar();
    }
    while(c>='0'&&c<='9') {
        x=x*10+c-'0';
        c=getchar();
    }
    return x*f;
}
int n,m,p;
namespace BF2 {
    const int N = 110;
    struct node {
        int a[N][N],n,m;
        node() {
            memset(a,0,sizeof(a));n = 0,m = 0;
        }
        node(int nn,int mm) {
            memset(a,0,sizeof(a));n = nn;m = mm;
        }
        node(int nn) {
            n = nn;
            memset(a,0,sizeof(a));
            for(int i = 0;i <= nn;++i) a[i][i] = 1;
        }
    };
    node operator * (const node &x,const node &y) {
        int n = x.n,m = y.m;
        node c(n,m);
        int k = x.m;
        for(int i = 0;i <= k;++i) {
            for(int j = 0;j <= n;++j) {
                for(int z = 0;z <= m;++z) {
                    c.a[j][z] += 1ll * x.a[j][i] * y.a[i][z] % mod;
                    c.a[j][z] %= mod;
                }
            }
        }
        return c;
    }
    node operator ^ (node x,int y) {
        node ans(x.n);
        for(;y;y >>= 1,x = x * x)
            if(y & 1) ans = ans * x;
        return ans;
    }
    int dis[2000003],vis[20000003];
    int tot,tmp[2000003];
    int js;
    void Eur() {
        vis[1] = 1;
        for(int i = 2;i <= m;++i) {
            if(!vis[i]) dis[++js] = i;
            for(int j = 1;1ll * dis[j] * i <= m;++j) {
                vis[dis[j] * i] = 1;
                if(dis[j] % i == 0) break;
            }
        }
    }
    void Main() {
        node C(p - 1,p - 1);
        int k = m % p;
        for(int i = 0;i < p;++i) {
            for(int j = i + 1;j <= i + k;++j) {
                C.a[i][j % p]++;
            }
        }
        for(int i = 0;i < p;++i) for(int j = 0;j < p;++j) C.a[i][j] += m / p;
        node ans(0,p - 1);
        ans.a[0][0] = 1;
        ans = ans * (C ^ n);
        int anss = ans.a[0][0];
        memset(ans.a,0,sizeof(ans.a));
        memset(C.a,0,sizeof(C.a));
        ans.a[0][0] = 1;
        Eur();
        for(int i = 0;i < p;++i) {
            for(int j = i + 1;j <= i + k;++j) {
                C.a[i][j % p]++;
            }
        }
        for(int i = 0;i < p;++i) for(int j = 0;j < p;++j) C.a[i][j] += m / p;
        for(int i = 0;i < p;++i)
            for(int j = 1;j <= js;++j)
                C.a[i][(i + dis[j]) % p]--;
        ans = ans * (C ^ n);
        cout<<(anss - ans.a[0][0] + mod) % mod;
    }

}
signed main() {
    n = read(),m = read(),p = read();
    BF2::Main();
    return 0;
}