题意

长为 $n$ 的序列，每个位置可以填 $[1, m]$ 的正整数，问有多少中填法，满足所有数之和是 $p$ 的倍数并且至少填了一个质数。

$n \leq 10^9, m \leq 2 \times 10^7, p \leq 100$

题目分析

假设序列第 $i$ 个数为 $a_i$，相当于要满足 $\sum a_i \equiv 0 \pmod p \implies \sum (a_i \bmod p) \equiv 0 \pmod p$。
要至少存在填一个质数，可以由简单的容斥转化为：任意填的方案数 - 不填质数的方案数。分成两次分别计算即可。

统计方案数，自然想到 DP。首先不难设计出这样的转移方程：

$$f(i, j) \times c (j) \to f(i+1, (j+k) \bmod p)$$

$f(i, j)$ 表示填到第 $i$ 个位置，当前模 $p$ 余 $j$ 的方案数，$c (j)$ 表示 $[1, m]$ 中 $x \bmod p = j$ 的 $x$ 个数。注意上面的容斥，求不填质数的方案数时 $c (j)$ 不考虑所有质数。

由于 $n$ 很大，考虑使用矩阵快速幂加速转移。
具体地：

• 构造转移矩阵 $trans$，满足 $trans_{i, j} = \sum_{(k+j) \bmod p = i} c (j)$。
• 维护一个答案向量 $F_i$，有：

$$F_i= \begin{bmatrix} f_{i, 0}\\ f_{i, 1}\\ \dots\\ f_{i, p-1} \end{bmatrix} $$

• 那么转移方程可改写为：

$$F_{i+1}= \begin{bmatrix} f_{i+1, 0}\\ f_{i+1, 1}\\ \dots\\ f_{i+1, p-1} \end{bmatrix} = trans \times \left( F_i= \begin{bmatrix} f_{i, 0}\\ f_{i, 1}\\ \dots\\ f_{i, p-1} \end{bmatrix} \right) $$

初始向量 $F_0=\begin{bmatrix}1\\0\\\dots\\0\end{bmatrix}$，将其左乘 $trans$ 共 $n$ 即得到最终答案向量 $F_n$，题目所需为 $f_{n, 0}$。

代码实现

#include <cstdio>
#include <iostream>
#include <cstring>

using std::cin;
using std::cout;

const int M = 2e7;
const int P = 100;
const int MOD = 20170408;

int n, m, p;

bool Vis[M + 5];

struct Mat {
  int n, m;
  int C[P + 5][P + 5];
}trans, zero, st;

int Ct[P + 5];

void Init();
int Solve (bool);
Mat Qpow (int);
Mat Mul (Mat, Mat);

int main() {
  std::ios::sync_with_stdio (0), cin.tie (0);

  cin >> n >> m >> p;
  Init();

  cout << ((Solve (true) - Solve (false)) % MOD + MOD) % MOD << '\n';
  return 0;

}void Init() {
  Vis[1] = true;
  for (int i = 2; i <= m; ++i)
    if (!Vis[i]) for (int j = i + i; j <= m; j += i) Vis[j] = true;

}int Solve (bool flag) {
  memset (Ct, 0, sizeof (Ct));
  memset (trans.C, 0, sizeof (trans.C));
  for (int i = 1; i <= m; ++i)
    if ((Vis[i]) || ((!Vis[i]) && flag)) ++Ct[i % p];
  
  for (int i = 0; i < p; ++i)
    for (int j = 0; j < p; ++j) trans.C[(i + j) % p][i] += Ct[j];
  trans.n = trans.m = zero.n = zero.m = p;
  st.n = p, st.m = st.C[0][0] = 1;

  Mat ans = Qpow (n);
  return ans.C[0][0];

}Mat Qpow (int k) {
  Mat res = st, base = trans;
  while (k) {
    if (k & 1) res = Mul (base, res);
    base = Mul (base, base), k >>= 1;
  }return res;

}Mat Mul (Mat u, Mat v) {
  Mat w = zero;
  w.n = u.n, w.m = v.m;
  for (int i = 0; i < w.n; ++i)
    for (int k = 0; k < u.m; ++k)
      if (u.C[i][k])
        for (int j = 0; j < v.m; ++j)
          w.C[i][j] = (w.C[i][j] + 1ll * u.C[i][k] * v.C[k][j] % MOD) % MOD;
  return w;
}