题意

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 9101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596
#include <cstdio>
#include <cstring>
#include <vector>
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