Alice 想要得到一个长度为 n n 的序列，序列中的数都是不超过 m m 的正整数，而且这 n n 个数的和是 p p 的倍数。

Alice 还希望，这 n n 个数中，至少有一个数是质数。

Alice 想知道，有多少个序列满足她的要求。

题解

考虑一个比较显然的 dp，设 f(i,j) f(i, j) 表示序列长度为 i i ，序列中元素模 p p 的值为 j j 时的方案数。这样考虑不到需要有质数的要求，那么考虑容斥，设 g(i,j) g(i, j) 表示的意义与 f(i,j) f(i, j) 相同，但是不允许选出的序列中的元素有质数。

这样最后的答案就是 f(n,0)−g(n,0) f(n, 0) - g(n, 0) 。

刚才 dp 的转移方程是平凡的：

可以观察到这是一个一维转移，可以使用矩阵乘法来优化转移。

考虑一个转移矩阵 A A 和初始矩阵 F F 。

考虑构造转移矩阵，若 F F 矩阵的形式如下：

我们以转移矩阵的第一行为例，若转移矩阵第一行形式如下：

我们应使得

∑j=0p−1aj f(i,j)=f(i+1,0) \sum\limits_{j = 0}^{p - 1} a_j\ f(i, j) = f(i + 1, 0)

也就是说，转移矩阵第一行的每个位置对应的应该是对于每个 j j 能够使得 j+k,k∈[1,m]=0 j + k, k \in [1, m] = 0 的方案数。余下几行同理。

在这里注意一个细节，我们并不需要枚举 [1,m] [1, m] 中的数来构造转移矩阵，我们可以预处理处 [1,m] [1, m] 中的数模 p p 后的结果，并统计每个结果的出现次数，这样枚举的复杂度就从 O(m) O(m) 减小到 O(p) O(p) 了。原来构造转移矩阵需要 O(mp) O(mp) 这样优化以后只需要 O(p2) O(p^2) 。

这是对于 f(i,j) f(i, j) 的处理。同理，如果要处理 g(i,j) g(i, j) 的话，只需要在构造转移矩阵的时候排除掉加上质数转移的情况，处理的方式是线筛出 [1,m] [1, m] 中的质数，然后在上一步统计 [1,m] [1, m] 中的数模 p p 所得结果出现次数时分质数和非质数来讨论就可以了。

代码

#include <cstdio>
#include <cstring>
#include <climits>
#include <algorithm>
const int MAX_N = 1e9;
const int MAX_M = 2e7;
const int MAX_P = 100;
const int MAX_PRIME = 1270607;
const int MOD = 20170408;
bool isNotPrime[MAX_M + 1];
int primes[MAX_PRIME + 1], cnt;
int n, m, p;
inline int get()
{
    int res = 1, Q = 1;
    char c;
    while ((c = getchar()) < 48 || c > 57)
    if (c == '-') Q = -1;
    if (Q) res = c-48;
    while ((c=getchar()) >= 48 && c <= 57)
    res = res * 10 + c - 48;
    return res * Q;
}
inline void sieve() {
    isNotPrime[0] = isNotPrime[1] = true;
    for (int i = 1; i <= MAX_M; i++) {
        if (!isNotPrime[i]) primes[cnt++] = i;
        for (int j = 0; j < cnt && (long long)i * primes[j] <= MAX_M; j++) {
            isNotPrime[i * primes[j]] = true;
            if (i % primes[j] == 0) break;
        }
    }
}
struct Matrix {
    long long a[MAX_P][MAX_P];
    int n, m;
    Matrix(int n, int m, bool unit = false) : n(n), m(m) {
        memset(a, 0, sizeof(a));
        if (unit) for (int i = 0; i < n; i++) a[i][i] = 1;
    }
};
Matrix operator*(const Matrix &a, const Matrix &b) {
    Matrix res(a.n, b.m);
    for (int i = 0; i < a.n; i++) {
        for (int j = 0; j < b.m; j++) {
            for (int k = 0; k < a.m; k++) {
                res.a[i][j] += a.a[i][k] * b.a[k][j] % MOD;
            }
        }
    }
    return res;
}
Matrix pow(Matrix a, int n) {
    Matrix res(a.n, a.n, true);
    for (; n; n >>= 1, a = a * a) if (n & 1) res = res * a;
    return res;
}
int primeCnt[MAX_P], nonPrimeCnt[MAX_P];
inline int calc(bool flag) {
    Matrix A(p, p), F(p, 1);
    for (int j = 0; j < p; j++) {
        for (int k = 0; k < p; k++) {
            A.a[(j + k) % p][j] = (primeCnt[k] * flag + nonPrimeCnt[k]) % MOD;
        }
    }
    A = pow(A, n);
    F.a[0][0] = 1;
    A = A * F;
    return A.a[0][0];
}
int main() {
    n = get(), m = get(), p = get();
    sieve();
    for (int i = 1; i <= m; i++) {
        if (isNotPrime[i]) nonPrimeCnt[i % p]++;
        else primeCnt[i % p]++;
    }
    printf("%d\n", (calc(true) - calc(false) + MOD) % MOD);
    return 0;
}