傳送門

不考慮質數的條件，可以考慮到一個很明顯的$DP:$設$f_{i,j}$表示選$i$個數，和$mod\ p=j$的方案數，顯然是可以矩陣優化$DP$的。

而且轉移矩陣是迴圈矩陣，所以可以只用第一行的數字代替整個矩陣。當然了這道題$p \leq 100$矩陣比較小也可以直接做。

然後考慮至少要一個質數的條件，發現就是所有數參與$DP$的答案減去所有合數參與$DP$的答案，兩次算出來相減即可。

 1 #include<bits/stdc++.h>
 2 #define ll long long
 3 //This code is written by Itst
 

 4 using namespace std;
 5 
 6 inline int read(){
 7     int a = 0;
 8     char c = getchar();
 9     bool f = 0;
10     while(!isdigit(c)){
11         if(c == '-')
12             f = 1;
13         c = getchar();
14     }
15     while(isdigit(c)){
16         a = (a << 3) + (a << 1) + (c ^ '
 
0');
17         c = getchar();
18     }
19     return f ? -a : a;
20 }
21 
22 const int MOD = 20170408;
23 int N , M , P , ans;
24 bool nprime[(int)2e7 + 10];
25 struct matrix{
26     ll a[110];
27     matrix(){memset(a , 0 , sizeof(a));}
28     inline ll& operator [](int x){return a[x];}
29     matrix operator
 
 *(matrix b){
30         matrix c;
31         for(int i = 0 ; i < P ; ++i)
32             for(int j = 0 ; j < P ; ++j)
33                 c[i] += a[j] * b[i - j < 0 ? i - j + P : i - j];
34         for(int j = 0 ; j < P ; ++j)
35             c[j] %= MOD;
36         return c;
37     }
38 }S , T , G;
39 
40 int main(){
41 #ifndef ONLINE_JUDGE
42     freopen("in" , "r" , stdin);
43     //freopen("out" , "w" , stdout);
44 #endif
45     N = read();
46     M = read();
47     P = read();
48     for(int i = 0 ; i < P && i <= M ; ++i)
49         G[i % P] = (M - i) / P + (bool)i;
50     S[0] = 1;
51     T = G;
52     int K = N;
53     while(K){
54         if(K & 1)
55             S = S * T;
56         T = T * T;
57         K >>= 1;
58     }
59     ans = S[0];
60     for(int i = 2 ; i <= M ; ++i)
61         if(!nprime[i]){
62             --G[i % P];
63             for(int j = i ; j <= M / i ; ++j)
64                 nprime[i * j] = 1;
65         }
66     T = G;
67     S = matrix();
68     S[0] = 1;
69     K = N;
70     while(K){
71         if(K & 1)
72             S = S * T;
73         T = T * T;
74         K >>= 1;
75     }
76     cout << (ans - S[0] + MOD) % MOD;
77     return 0;
78 }