P7722 [Ynoi2007] tmpq 题解

前言

这篇题解是 zak 的做法，我实现得很丑所以常数巨大，但我尽量把我的实现讲清楚（）。

做法

转化 & 暴力

首先，令 \(b'_i=b_{a_i},c'_i=c_{a_i}\)，就可以将原题条件转换为 \(b'_i=a_i=c'_i\)，那么对 \(a_i\) 进行单点修改，其实就相当于对 \(b'_i\) 和 \(c'_i\) 都进行单点修改。

对于每一种 \(b'_i=a_i=c'_i=w\) 的 \(w\)，他们的答案数是互相独立互不干扰的，所以可以对每种 \(w\) 单独考虑对答案的贡献。

设 \(dp_{i,1/2/3}\) 表示考虑前 \(i\) 个数中，\(b'_p=w/b'_p=a_q=w/b'_p=a_q=c'_r=w\) 的方案数（\(p<q<r\le i\)），其中 \(dp_{i,0}=1\)，则有

\[dp_{i,1}=dp_{i-1,1}+dp_{i-1,0}\times [b'_i=w] \]

\[dp_{i,2}=dp_{i-1,2}+dp_{i-1,1}\times [a_i=w] \]

\[dp_{i,3}=dp_{i-1,3}+dp_{i-1,2}\times [c'_i=w] \]

共有 \(n\) 种 \(w\)，每次修改可能导致 \(6\) 种 \(w\) 的答案改变，暴力做一次 DP 是 \(O(n)\) 的，所以总复杂度是 \(O(n^2+6nm)\)

由于 \(dp_i\) 只依赖于 \(dp_{i-1}\)，所以可以空间压缩成 \(dp_{1/2/3}\)，注意转移顺序。

根号分治：小于等于根号的部分

对 \(w\) 进行根号分治，统计出 \(w\) 在 \(a_i,b'_i,c'_i\) 以及修改后，\(w\) 的总出现次数 \(cnt_w\)（显然 \(\sum cnt_w=3nm\)）。按照 \(cnt_w\) 的大小根号分治。

对于出现次数 \(<B\) 的 \(w\)，能够有效更改 \(dp\) 的位置不超过 \(B\) 个，因此每次涉及到修改 \(w\) 的操作时，暴力重做这一次 DP。

使用动态数组 \(vec_w\) 存下这 \(B\) 个位置，当进行修改操作时，相当于往 \(vec_w\) 里添加或删除位置，暴力找到这个位置进行操作，时间复杂度 \(O(B)\)。

接着考虑维护前缀和，令 \(f_i\) 表示 \(c'_i=w\) 且符合题意的方案数，则 \(f_i=dp_{i,2}\times[c'_i=w]\)，要维护的就是 \(f_i\) 的前缀和。

对序列进行分块，考虑到只需要进行 \(m\) 次查询，而暴力重新做 DP 本身已经占据了 \(O(B)\) 的时间，所以需要 \(O(1)\) 修改，\(O(\sqrt n)\) 查询（这里的根号只是象征义）。令 \(s_i\) 表示块 \(i\) 的 \(f\) 的和就可以做到。

根号分治：大于根号的部分

这个 DP 只有四个状态，所以考虑序列动态 DP。

\[\begin{bmatrix} dp_0 & dp_1 & dp_2 & dp_3 \end{bmatrix} \times \begin{bmatrix} 1 & b'_i=w & 0 & 0 \\ 0 & 1 & a_i=w & 0 \\ 0 & 0 & 1 & c'_i=w \\ 0 & 0 & 0 & 1 \end{bmatrix} = \begin{bmatrix} dp'_0 & dp'_1 & dp'_2 & dp'_3 \end{bmatrix} \]

由于出现次数 \(>B\) 的 \(w\) 小于 \(\frac{n+m}{B}\) 个，所以单独考虑每一个 \(w\)，\(O(n)\) 维护矩阵前缀乘。

由于枚举每种 \(w\) 已经消耗了 \(\sqrt n\) 的时间，查询数有 \(m\) 个，但是总修改数一定也是 \(m\) 的级别的（和上面 \(O(6mn)\) 一样的道理），所以需要做到 \(O(1)\) 查询，\(O(\sqrt n)\) 修改。设 \(SB_i\) 表示前 \(i\) 个块的前缀乘，\(S_i\) 表示块的左端点到 \(i\) 的前缀乘，那么查询 \(r\) 就是 \(SB_{bel_r-1}\times S_r\)，修改按照定义修改即可。

取矩阵中的哪个数呢？初始的 \(\begin{bmatrix}dp_0 & dp_1 & dp_2 & dp_3\end{bmatrix}\) 是 \(\begin{bmatrix}1&0&0&0\end{bmatrix}\)，最终希望得到 \(dp_3\)，带进去算可以得到是要矩阵 \((0,3)\) 这个位置的数。

总时间复杂度 \(O((n+m)B+4^3\times \frac{n+m}{B}\times(n+m))\)，空间复杂度线性，\(B\) 可以取 \(\sqrt{n+m}\)，但我取的 \(B=1000\) 才能通过（）。

我错的一些细节

可能会出现 \(cnt_{c'_i}\le B\)，但是修改后的 \(cnt_{c'_i}>B\)，这样他就不会重做 \(c'_i\) 的 DP，导致 \(f_i\) 不能被消除贡献。

代码细节有点多，都打注释了。

#include <bits/stdc++.h>
#define ll long long
#define pii pair<int, int>
#define fi first
#define se second
#define IOS ios::sync_with_stdio(0), cin.tie(0), cout.tie(0)
using namespace std;
const int N = 2e5 + 10, B = 210, M = 5e4 + 10;
int n, m, a[N], aa[N], b[N], c[N], cnt[N];
ll ans[M];
struct QUES {
    int opt, x, y;
} q[M];

int block, nb, L[B], R[B], bel[N];
// 预处理分块
// 我根号分治和序列分块共用了 block
inline void init() {
    block = 1000;
    for (int i = 1; i <= n; i++)
        bel[i] = (i + block - 1) / block;
    nb = bel[n];
    for (int i = 1; i <= nb; i++)
        L[i] = R[i - 1] + 1, R[i] = i * block;
    R[nb] = n;
}

namespace Small {
    ll dp[4], f[N], s[N];
    vector<pii> vec[N];
    // 转移顺序：不能同一位的 1 到 2
    bool cmp(pii x, pii y) {
        return x.fi != y.fi ? x.fi < y.fi : x.se > y.se;
    }
    // 暴力寻找一个数并添加/删除
    inline void add(int id, pii x) {
        auto it = lower_bound(vec[id].begin(), vec[id].end(), x, cmp);
        vec[id].insert(it, x);
    }
    inline void del(int id, pii x) {
        // 提前把贡献消除
        if (x.se == 3) {
            s[bel[x.fi]] -= f[x.fi];
            f[x.fi] = 0;
        }
        auto it = lower_bound(vec[id].begin(), vec[id].end(), x, cmp);
        vec[id].erase(it);
    }
    // 重做 w=id 的 DP
    inline void work(int id) {
        dp[0] = 1;
        dp[1] = dp[2] = dp[3] = 0;
        for (auto [x, y] : vec[id]) {
            dp[y] += dp[y - 1];
            s[bel[x]] -= f[x];
            if (y == 3)
                f[x] = dp[2];
            s[bel[x]] += f[x];
        }
    }
    // 求 f 的前缀和
    inline ll query(int x) {
        ll ret = 0;
        for (int i = 1; i < bel[x]; i++)
            ret += s[i];
        for (int i = L[bel[x]]; i <= x; i++)
            ret += f[i];
        return ret;
    }
    inline void solve() {
        for (int i = 1; i <= n; i++) {
            if (cnt[b[a[i]]] <= block)
                add(b[a[i]], {i, 1});
            if (cnt[a[i]] <= block)
                add(a[i], {i, 2});
            if (cnt[c[a[i]]] <= block)
                add(c[a[i]], {i, 3});
        }
        for (int i = 1; i <= n; i++)
            work(i);
        for (int i = 1; i <= m; i++) {
            auto [opt, x, y] = q[i];
            if (opt == 1) {
                if (cnt[b[a[x]]] <= block)
                    del(b[a[x]], {x, 1}), work(b[a[x]]);
                if (cnt[a[x]] <= block)
                    del(a[x], {x, 2}), work(a[x]);
                if (cnt[c[a[x]]] <= block)
                    del(c[a[x]], {x, 3}), work(c[a[x]]);
                a[x] = y;
                if (cnt[b[a[x]]] <= block)
                    add(b[a[x]], {x, 1}), work(b[a[x]]);
                if (cnt[a[x]] <= block)
                    add(a[x], {x, 2}), work(a[x]);
                if (cnt[c[a[x]]] <= block)
                    add(c[a[x]], {x, 3}), work(c[a[x]]);
            } else
                ans[i] += query(x);
        }
    }
}

namespace Big {
    // cnt[w] > B 的数集
    vector<int> big;
    struct matrix {
        ll d[4][4];
        // 因为是前缀乘，所以初始化为单位矩阵
        matrix() {
            memset(d, 0, sizeof d);
            for (int i = 0; i < 4; i++)
                d[i][i] = 1;
        }
        matrix(bool I) {
            memset(d, 0, sizeof d);
            for (int i = 0; I and i < 4; i++)
                d[i][i] = 1;
        }
    } mat, S[N], SB[B];
    // 实现太垃圾不展开过不了 QAQ
    inline matrix operator*(const matrix &x, const matrix &y) {
        matrix ret;
        ret.d[0][0] = x.d[0][0] * y.d[0][0] + x.d[0][1] * y.d[1][0] + x.d[0][2] * y.d[2][0] + x.d[0][3] * y.d[3][0];
        ret.d[0][1] = x.d[0][0] * y.d[0][1] + x.d[0][1] * y.d[1][1] + x.d[0][2] * y.d[2][1] + x.d[0][3] * y.d[3][1];
        ret.d[0][2] = x.d[0][0] * y.d[0][2] + x.d[0][1] * y.d[1][2] + x.d[0][2] * y.d[2][2] + x.d[0][3] * y.d[3][2];
        ret.d[0][3] = x.d[0][0] * y.d[0][3] + x.d[0][1] * y.d[1][3] + x.d[0][2] * y.d[2][3] + x.d[0][3] * y.d[3][3];

        ret.d[1][0] = x.d[1][0] * y.d[0][0] + x.d[1][1] * y.d[1][0] + x.d[1][2] * y.d[2][0] + x.d[1][3] * y.d[3][0];
        ret.d[1][1] = x.d[1][0] * y.d[0][1] + x.d[1][1] * y.d[1][1] + x.d[1][2] * y.d[2][1] + x.d[1][3] * y.d[3][1];
        ret.d[1][2] = x.d[1][0] * y.d[0][2] + x.d[1][1] * y.d[1][2] + x.d[1][2] * y.d[2][2] + x.d[1][3] * y.d[3][2];
        ret.d[1][3] = x.d[1][0] * y.d[0][3] + x.d[1][1] * y.d[1][3] + x.d[1][2] * y.d[2][3] + x.d[1][3] * y.d[3][3];

        ret.d[2][0] = x.d[2][0] * y.d[0][0] + x.d[2][1] * y.d[1][0] + x.d[2][2] * y.d[2][0] + x.d[2][3] * y.d[3][0];
        ret.d[2][1] = x.d[2][0] * y.d[0][1] + x.d[2][1] * y.d[1][1] + x.d[2][2] * y.d[2][1] + x.d[2][3] * y.d[3][1];
        ret.d[2][2] = x.d[2][0] * y.d[0][2] + x.d[2][1] * y.d[1][2] + x.d[2][2] * y.d[2][2] + x.d[2][3] * y.d[3][2];
        ret.d[2][3] = x.d[2][0] * y.d[0][3] + x.d[2][1] * y.d[1][3] + x.d[2][2] * y.d[2][3] + x.d[2][3] * y.d[3][3];

        ret.d[3][0] = x.d[3][0] * y.d[0][0] + x.d[3][1] * y.d[1][0] + x.d[3][2] * y.d[2][0] + x.d[3][3] * y.d[3][0];
        ret.d[3][1] = x.d[3][0] * y.d[0][1] + x.d[3][1] * y.d[1][1] + x.d[3][2] * y.d[2][1] + x.d[3][3] * y.d[3][1];
        ret.d[3][2] = x.d[3][0] * y.d[0][2] + x.d[3][1] * y.d[1][2] + x.d[3][2] * y.d[2][2] + x.d[3][3] * y.d[3][2];
        ret.d[3][3] = x.d[3][0] * y.d[0][3] + x.d[3][1] * y.d[1][3] + x.d[3][2] * y.d[2][3] + x.d[3][3] * y.d[3][3];
        return ret;
    }
    // 先 O(n) 做一次 w=x 的
    inline void work(int x) {
        for (int i = 1; i <= n; i++) {
            int bi = bel[i];
            mat.d[0][1] = b[a[i]] == x;
            mat.d[1][2] = a[i] == x;
            mat.d[2][3] = c[a[i]] == x;
            if (i == L[bi])
                S[i] = mat;
            else
                S[i] = S[i - 1] * mat;
        }
        for (int i = 1; i <= nb; i++)
            SB[i] = SB[i - 1] * S[R[i]];
    }
    // O(sqrt(n)) 更新 w=x 位于 id 处的修改
    inline void update(int x, int id) {
        int bid = bel[id];
        for (int i = id; i <= R[bid]; i++) {
            mat.d[0][1] = b[a[i]] == x;
            mat.d[1][2] = a[i] == x;
            mat.d[2][3] = c[a[i]] == x;
            if (i == L[bid])
                S[i] = mat;
            else
                S[i] = S[i - 1] * mat;
        }
        for (int i = bid; i <= nb; i++)
            SB[i] = SB[i - 1] * S[R[i]];
    }
    // 前缀答案
    inline ll query(int x) {
        return (SB[bel[x] - 1] * S[x]).d[0][3];
    }
    inline void solve() {
        for (int i = 1; i <= n; i++)
            if (cnt[i] > block)
                big.push_back(i);
        int cc = 0;
        for (int x : big) {
            // 先清零
            for (int i = 1; i <= nb; i++)
                S[i] = SB[i] = matrix();
            for (int i = 1; i <= n; i++)
                a[i] = aa[i];
            work(x);
            for (int i = 1; i <= m; i++) {
                auto [opt, id, v] = q[i];
                if (q[i].opt == 1) {
                    // 有关联的才更新，不然复杂度不对
                    if (a[id] == x or b[a[id]] == x or c[a[id]] == x or v == x or b[v] == x or c[v] == x) {
                        a[id] = v;
                        update(x, id);
                        cc++;
                    } else
                        a[id] = v;
                } else
                    ans[i] += query(id);
            }
        }
    }
}

signed main() {
    IOS;
    cin >> n >> m;
    init();
    for (int i = 1; i <= n; i++) {
        cin >> a[i];
        aa[i] = a[i];
    }
    for (int i = 1; i <= n; i++)
        cin >> b[i];
    for (int i = 1; i <= n; i++)
        cin >> c[i];
    for (int i = 1; i <= n; i++) {
        cnt[a[i]]++, cnt[b[a[i]]]++, cnt[c[a[i]]]++;
    }
    for (int i = 1; i <= m; i++) {
        cin >> q[i].opt >> q[i].x;
        if (q[i].opt == 1) {
            cin >> q[i].y;
            cnt[q[i].y]++;
            cnt[b[q[i].y]]++;
            cnt[c[q[i].y]]++;
        }
    }
    Small::solve();
    Big::solve();
    for (int i = 1; i <= m; i++)
        if (q[i].opt == 2)
            cout << ans[i] << "\n";
    return 0;
}