题意

Alice 想要得到一个长度为 \(n\) 的序列,序列中的数都是不超过 \(m\) 的正整数,而且这 \(n\) 个数的和是 \(p\) 的倍数。Alice 还希望,这 \(n\) 个数中,至少有一个数是质数。Alice 想知道,有多少个序列满足她的要求,答案对 \(20170408\) 取模。

其中,\(1 \leq n \leq 10^9\)\(1 \leq m \leq 2 \times 10^7\)\(1 \leq p \leq 100\)

题解

由于 \(p\) 的取值范围较小,考虑 DP。

\(a(i, j)\) 代表序列前 \(i\) 个数的和模 \(p\)\(j\) 的序列个数。

\(b(i, j)\) 代表序列前 \(i\) 个数的和模 \(p\)\(j\) 且这 \(i\) 个数不为质数的序列个数。

采取矩阵快速幂加速,答案为 \(a(n, 0) - b(n, 0)\)

代码

 1 2 3 4 5 6 7 8 91011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
int n, m, p;
const int MOD = 20170408;
const int MAX_P = 105;
int ansa[MAX_P], ansb[MAX_P];
const int SIV_N = 2e7 + 5;
int not_prime[SIV_N];
vector<int> prime;
void linear_sieve() {
    memset(not_prime, 0, sizeof(not_prime));
    prime.clear();
    for (int i = 2; i <= m; i++) {
        if (!not_prime[i]) {
            prime.push_back(i);
            ansb[i % p]--;
        }
        for (int j = 0; j < prime.size() && i * prime[j] <= m; j++) {
            not_prime[i * prime[j]] = 1;
            if (i % prime[j] == 0)
                break;
        }
    }
}
const int MAT_N = 105;
const int MAT_MOD = 20170408;
struct Matrix {
    int v[MAT_N][MAT_N], n;
    Matrix(int _n = 0, int op = 0): n(_n) {
        memset(v, 0, sizeof(v));
        if (op)
            for (int i = 0; i < n; i++)
                v[i][i] = 1;
    }
    Matrix operator*(const Matrix &mat) const {
        Matrix ret(n);
        for (int i = 0; i < n; i++)
            for (int j = 0; j < n; j++)
                for (int k = 0; k < n; k++)
                    ret.v[i][j] = (ret.v[i][j] + (ll)v[i][k] * mat.v[k][j] % MAT_MOD) % MAT_MOD;
        return ret;
    }
    friend Matrix operator^(Matrix mat, ll n) {
        Matrix ret(mat.n, 1);
        while (n) {
            if (n & 1)
                ret = ret * mat;
            n >>= 1;
            mat = mat * mat;
        }
        return ret;
    }
};
int get(int x, int y) {
    if (x <= y)
        return y - x;
    return p - x + y;
}
int main() {
    scanf("%d%d%d", &n, &m, &p);
    for (int t = m % p, i = 0; i < p; i++) {
        ansa[i] = m / p;
        if (i != 0 && t >= i)
            ansa[i]++;
        ansb[i] = ansa[i];
    }
    linear_sieve();
    if (n == 1) {
        printf("%d\n", ansa[0] - ansb[0]);
        return 0;
    }
    int ans;
    Matrix fcta(p);
    for (int i = 0; i < p; i++)
        fcta.v[i][0] = ansa[i];
    Matrix fctb(p);
    for (int i = 0; i < p; i++)
        fctb.v[i][0] = ansb[i];
    Matrix a(p);
    for (int i = 0; i < p; i++)
        for (int j = 0; j < p; j++)
            a.v[i][j] = ansa[get(j, i)];
    a = a ^ (n - 1);
    a = a * fcta;
    Matrix b(p);
    for (int i = 0; i < p; i++)
        for (int j = 0; j < p; j++)
            b.v[i][j] = ansb[get(j, i)];
    b = b ^ (n - 1);
    b = b * fctb;
    ans = ((a.v[0][0] - b.v[0][0]) % MOD + MOD) % MOD;
    printf("%d\n", ans);
    return 0;
}
comments powered by Disqus