[Sdoi2017]序列计数 矩阵优化dp

题目

https://www.lydsy.com/JudgeOnline/problem.php?id=4818

思路

先考虑没有质数限制
dp是在同余系下的，所以(f[i][j])表示前i个点，和为j的方案数
转移就是(f[i][j]=f[i-1][k]+g[(j-k)\%p])
g[i]是x%p==i出现的个数
有质数的话
用tot-无质数
无质数就在g[i]上删去质数出现的个数，再跑一边
但是！！
n很大，应该是带个log的
矩阵优化吧

代码

#include <bits/stdc++.h>
using namespace std;
const int mod = 20170408;
const int N = 1e5 + 7;
int read() {
    int x = 0, f = 1;
    char s = getchar();
    for (; s > '9' || s < '0'; s = getchar())
        if (s == '-')
            f = -1;
    for (; s >= '0' && s <= '9'; s = getchar()) x = x * 10 + s - '0';
    return x * f;
}
int n, m, p;
int pri[2000008], cnt;
bool vis[20000008];
void get_pri(int n) {
    vis[1] = 1;
    for (int i = 2; i <= n; ++i) {
        if (!vis[i])
            pri[++cnt] = i;
        for (int j = 1; i * pri[j] <= n && j <= cnt; ++j) {
            vis[i * pri[j]] = 1;
            if (i % pri[j] == 0)
                break;
        }
    }
}
struct node {
    int ch[201][201];
} a, b, c;
node mul(node a, node b) {
    node c = {};
    for (int i = 0; i < p; ++i)
        for (int k = 0; k < p; ++k)
            for (int j = 0; j < p; ++j) c.ch[i][j] = (c.ch[i][j] + 1LL * a.ch[i][k] * b.ch[k][j] % mod) % mod;
    return c;
}
node q_pow(node a, int nb) {
    node ans = {};
    for (int i = 0; i < p; ++i) ans.ch[i][i] = 1;
    while (nb) {
        if (nb & 1)
            ans = mul(ans, a);
        a = mul(a, a);
        nb >>= 1;
    }
    return ans;
}
int g[200];
int main() {
    // freopen("count.in","r",stdin);
    // freopen("count.out","w",stdout);
    n = read(), m = read(), p = read();
    get_pri(m);
    // one
    for (int i = 1; i <= m; ++i) g[i % p]++;
    for (int i = 0; i < p; ++i)
        for (int j = 0; j < p; ++j) a.ch[i][j] = g[(i - j + p) % p];
    for (int i = 0; i < p; ++i) b.ch[i][0] = g[i];
    c = mul(q_pow(a, n - 1), b);
    int ans = c.ch[0][0];
    memset(g, 0, sizeof(g));
    memset(a.ch, 0, sizeof(a.ch));
    memset(b.ch, 0, sizeof(b.ch));
    memset(c.ch, 0, sizeof(c.ch));

    // two
    for (int i = 1; i <= m; ++i)
        if (vis[i])
            g[i % p]++;
    for (int i = 0; i < p; ++i)
        for (int j = 0; j < p; ++j) a.ch[i][j] = g[(i - j + p) % p];
    for (int i = 0; i < p; ++i) b.ch[i][0] = g[i];
    c = mul(q_pow(a, n - 1), b);
    ans -= c.ch[0][0];
    ans = (ans % mod + mod) % mod;
    printf("%d",ans);
    return 0;
}
/*
3 5 3
*/