题意分析
首先，“至少一个数是质数”就很不好判断，考虑容斥掉。答案即 不考虑是否为质数的答案 减去 没有质数的答案。
观察到 $1\leq p\leq100$，考虑复杂度与 $p$ 正相关的算法。
设 $f_{i,j}$ 表示第 $1\sim i$ 个数，和模 $p$ 为 $j$ 时的不考虑是否为质数答案。没有质数的答案为 $g_{i,j}$。答案即：
\[(f_{n,0}-g_{n,0})\bmod20170408\]设 $\textit{cntF}{i}$ 表示 $1\sim m$ 中模 $p$ 余 $i$ 的数的个数，$\textit{cntG}{i}$ 表示 $1\sim m$ 中模 $p$ 余 $i$ 的合数的个数。
则有：
\[f_{i,j}=\sum_{k=0}^{p-1}f_{i-1,k}\cdot \textit{cntF}_{j-k}\\ g_{i,j}=\sum_{k=0}^{p-1}g_{i-1,k}\cdot \textit{cntG}_{j-k}\\\]即将两个序列对位相加，满足和模 $p$ 余 $0$ 即可。
但是这样复杂度会炸掉（$\mathcal O\left(np^2\right)$）。容易发现，序列里的数不定，因此一半的答案和后一半的答案一样，因此可以分治。
即枚举前一半的序列和后一半的序列，只要相加的和模 $p$ 余 $0$ 即可：
\[f_{i,j}=\sum_{k=0}^{p-1}f_{\frac i2,k}f_{\frac i2,(j-k)\bmod p}\\ g_{i,j}=\sum_{k=0}^{p-1}g_{\frac i2,k}g_{\frac i2,(j-k)\bmod p}\\\]注意 $i$ 为奇数时需要特判。
同时，滚动数组可以将空间复杂度优化至线性。
AC 代码
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
//#include<bits/stdc++.h>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<iomanip>
#include<cstdio>
#include<string>
#include<vector>
#include<cmath>
#include<ctime>
#include<deque>
#include<queue>
#include<stack>
#include<list>
using namespace std;
constexpr const int N=1e9,M=2e7,P=20170408;
int vis[M+1],prime[M+1],primeSize;
void pre(){
for(int i=2;i<=M;i++){
if(!vis[i]){
vis[i]=i;
prime[++primeSize]=i;
}
for(int j=1;j<=primeSize&&1ll*i*prime[j]<=M;j++){
vis[i*prime[j]]=prime[j];
if(i%prime[j]==0){
break;
}
}
}
}
int p,f[2][M+1],g[2][M+1],cntF[M+1],cntG[M+1];
void solve(int x,bool mode){
if(x==1){
for(int i=0;i<p;i++){
f[mode][i]=cntF[i];
g[mode][i]=cntG[i];
}
return;
}
solve(x>>1,!mode);
for(int i=0;i<p;i++){
f[mode][i]=g[mode][i]=0;
}
for(int i=0;i<p;i++){
for(int j=0;j<p;j++){
int pl=i-j;
if(i-j<0){
pl+=p;
}
f[mode][i]=(f[mode][i]+1ll*f[!mode][j]*f[!mode][pl]%P)%P;
g[mode][i]=(g[mode][i]+1ll*g[!mode][j]*g[!mode][pl]%P)%P;
}
}
if(x&1){
for(int i=0;i<p;i++){
f[!mode][i]=f[mode][i];
g[!mode][i]=g[mode][i];
}
for(int i=0;i<p;i++){
f[mode][i]=g[mode][i]=0;
}
for(int i=0;i<p;i++){
for(int j=0;j<p;j++){
int pl=i-j;
if(pl<0){
pl+=p;
}
f[mode][i]=(f[mode][i]+1ll*f[!mode][j]*cntF[pl]%P)%P;
g[mode][i]=(g[mode][i]+1ll*g[!mode][j]*cntG[pl]%P)%P;
}
}
}
}
int query(int n,int m,int p){
for(int i=1;i<=m;i++){
cntF[i%p]++;
cntG[i%p]+=(vis[i]!=i);
}
solve(n,0);
int ans=(f[0][0]-g[0][0])%P;
if(ans<0){
ans+=P;
}
return ans;
}
int main(){
/*freopen("test.in","r",stdin);
freopen("test.out","w",stdout);*/
ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);
pre();
int n,m;
cin>>n>>m>>p;
cout<<query(n,m,p)<<'\n';
cout.flush();
/*fclose(stdin);
fclose(stdout);*/
return 0;
}
