BZOJ4818 序列计数

4818: [Sdoi2017]序列计数

Time Limit: 30 Sec  Memory Limit: 128 MB

Description

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

，这n个数中，至少有一个数是质数。Alice想知道，有多少个序列满足她的要求。

Input

一行三个数，n,m,p。

1<=n<=10^9,1<=m<=2×10^7,1<=p<=100

Output

一行一个数，满足Alice的要求的序列数量，答案对20170408取模。

Sample Input

3 5 3

Sample Output

33

---

月考后日常颓水题，好久没有1A了，可能SDOI比较良心

SBDP，D[i]表示和为i的方案数，用所有的减去没有素数的，写出方程发现可以矩阵转移，然后乱敲

---

#include<bits/stdc++.h>
using namespace std;
template <class _T> inline void read(_T &_x) {
    int _t; bool flag = false;
    while ((_t = getchar()) != '-' && (_t < '0' || _t > '9')) ;
    if (_t == '-') _t = getchar(), flag = true; _x = _t - '0';
    while ((_t = getchar()) >= '0' && _t <= '9') _x = _x * 10 + _t - '0';
    if (flag) _x = -_x;
}
using namespace std;
typedef long long LL;
const int mod = 20170408;
const int maxv = 200;
struct Mat {
    int n, m, a[maxv][maxv];
    Mat() {}
    Mat(int x, int y):n(x), m(y) {
        for (register int i = 0, j; i < n; ++i)
            for (j = 0; j < m; ++j)
                a[i][j] = 0;
    }
    inline void init() {for (int i = 0; i < n; ++i) a[i][i] = 1; }
    inline Mat operator * (Mat B) {
        Mat C(n, B.m);
        for (register int i = 0, j, k; i < n; ++i)
            for (j = 0; j < B.m; ++j)
                for (k = 0; k < m; ++k) {
                    C.a[i][j] += (int)((LL)a[i][k] * B.a[k][j] % mod);
                    if (C.a[i][j] >= mod) C.a[i][j] -= mod;
                }
        return C;
    }
    inline Mat operator ^ (int t) {
        Mat res(n, m), tmp = *this; res.init();
        while (t) {
            if (t & 1) res = res * tmp;
            tmp = tmp * tmp, t >>= 1;
        }
        return res;
    }
};
const int maxp = 210;
const int maxm = 20000010;
int n, m, p;
int cnt_p[maxp], cnt_n[maxp];
bool vis[maxm];
int prime[maxm / 10], pcnt;
inline void Init() {
    cnt_p[1 % p] = 1;
    for (register int i = 2, j; i <= m; ++i) {
        if (!vis[i]) {
            prime[++pcnt] = i;
        } else {
            ++cnt_p[i % p];
        }
        for (j = 1; j <= pcnt && i * prime[j] <= m; ++j) {
            vis[i * prime[j]] = true;
            if (i % prime[j] == 0) break;
        }
    }
    int tmpa = m / p, tmpb = m % p;
    for (register int i = 0; i < p; ++i) {
        cnt_n[i] = tmpa;
        cnt_n[i] += (i && i <= tmpb);
    }
}
inline int getres(Mat &a) {
    Mat x(p, 1);
    x.a[0][0] = 1;
    return ((a ^ n) * x).a[0][0];
}
int main() {
    //freopen();
    //freopen();
    read(n), read(m), read(p);
    Init();
    Mat a(p, p), b(p, p);
    for (int i = 0, j, to; i < p; ++i) {
        for (j = 0; j < p; ++j) {
            to = i + j;
            if (to >= p) to -= p;
            a.a[i][to] = cnt_n[j];
            b.a[i][to] = cnt_p[j];
        }
    }
    int ans = getres(a) - getres(b);
    if (ans < 0) ans += mod;
    cout << ans << endl;
    return 0;
}