「Bzoj 4818」「SDOI2017」序列计数 (DP + 矩阵快速幂优化)

BZOJ 4818
题意：求$n$长序列值域$[1,m]$，求序列和是$p$倍数并且序列至少有一个质数的序列个数。

转化题意，则
$\text{ans}=$满足和为$p$的倍数的方案数$−$满足和为$p$的倍数且不含质数的方案数

那么我们先考虑满足和为$p$的倍数的方案数怎么算，显然我们可以设$dp(i,j)$为前$i$个和模$p=j$的方案，然后可以这样转移

$$
dp(i,j)=\sum_{k=1}^m dp(i-1,(j-k) \mod p)
$$

考虑优化，我们想办法让它与$m$复杂度无关，我们发现$i$位填$u$和填$u+pv$效果是一样的，所以我们直接记$cnt_i$表示$1∼m$中$\mod p$为$i$的个数，那么转移为

$$
dp(i,j)=\sum_{k=0}^{p-1} dp(i-1,(j-k) \mod p) \times cnt_k
$$

那么现在复杂度为$O(np^2)$

发现$n$很大，考虑矩阵优化到$O(\log n \times p^3)$

矩阵构造就循环一下构造就行，然后就是常规的优化方法

注意vis数组要开成bool

#include<cstdio> 
#include<cstring>
#include<algorithm>
#include<iostream>
#include<cmath>
#include<map>
#include<queue>
#include<set>
#define ms(i, j) memset(i, j, sizeof i)
#define LL long long
#define db double
#define fir first
#define sec second
#define mp make_pair
using namespace std;

namespace flyinthesky {

    const LL MO = 20170408;

    LL n, m, p;
    LL tot, pri[1400000 + 2];
    bool vis[20000000 + 2];
    LL gg1[100 + 2], gg2[100 + 2];

    struct matrix {
        LL a[100 + 2][100 + 2];
        matrix() {for (int i = 0; i <= 100; ++i) for (int j = 0; j <= 100; ++j) a[i][j] = 0;}
        void init() {for (int i = 0; i <= 100; ++i) a[i][i] = 1;}
        void clean() {for (int i = 0; i <= 100; ++i) for (int j = 0; j <= 100; ++j) a[i][j] = 0;}
    }f, ma1, ma2, orz1, orz2;

    void mul(matrix &a, matrix &b, matrix &ret) {
        ret.clean();
        for (LL i = 0; i <= p; ++i) {
            for (LL j = 0; j <= p; ++j) {
                for (LL k = 0; k <= p; ++k) {
                    (ret.a[i][j] += a.a[i][k] * b.a[k][j] % MO) %= MO;
                }
            }
        }
    }
    matrix ksm(matrix &a, LL b) {
        matrix ans, bs = a, tmp;
        ans.init();
        int fl = 0;
        while (b) {
            if (b & 1) {
                if (!fl) fl = 1, ans = bs;
                else mul(ans, bs, tmp), ans = tmp;
            }
            mul(bs, bs, tmp), bs = tmp;
            b >>= 1;
        }
        return ans;
    }

    void clean() {
    }
    int solve() {

        clean();
        cin >> n >> m >> p;

        vis[1] = 1;
        for (int i = 2; i <= m; ++i) {
            if (!vis[i]) pri[++tot] = i;
            for (int j = 1; j <= tot && pri[j] * i <= m; ++j) {
                vis[pri[j] * i] = 1;
                if (i % pri[j] == 0) {
                    break ;
                }
            }
        }

        gg1[0] = m / p;
        for (int i = 1; i < p; ++i) gg1[i] = m / p + (m % p >= i);
        for (int i = 1; i <= m; ++i) if (vis[i]) gg2[i % p]++;

        f.a[0][0] = 1;
        for (int i = 0; i < p; ++i) {
            for (int j = 0; j < p; ++j) {
                ma1.a[j][i] = gg1[((p + i - j) % p + p) % p];
                ma2.a[j][i] = gg2[((p + i - j) % p + p) % p];
            }
        }

        ma1 = ksm(ma1, n);
        ma2 = ksm(ma2, n);

        mul(f, ma1, orz1), mul(f, ma2, orz2);

        cout << ((orz1.a[0][0] - orz2.a[0][0]) % MO + MO) % MO;

        return 0;
    }  
}
int main() {
    flyinthesky::solve();
    return 0;
}