20250902

T1

相似子矩阵

显然只需要对 1 的总数的每个因数求答案。他都数据随机了，还放 T1，那感性理解一下显然直接 \(\mathcal{O}(n^3d(s))\) 暴力就是能过的。

代码

#include <iostream>
#include <string.h>
using namespace std;
int n;
int a[205][205], b[205][205];
int p[205], cnt[40005];
int d[105], ans[40005], dcnt;
int main() {
    freopen("maware.in", "r", stdin);
    freopen("maware.out", "w", stdout);
    int tc, s;
    cin >> tc;
    while (tc--) {
        memset(ans, 0, sizeof ans);
        cin >> n; s = dcnt = 0;
        for (int i = 1; i <= n; i++) for (int j = 1; j <= n; j++) cin >> a[i][j], s += a[i][j], b[i][j] = b[i - 1][j] + a[i][j];
        for (int i = 1; i < s; i++) if (s % i == 0) d[++dcnt] = i;
        for (int i = 0; i < n; i++) {
            for (int j = i + 2; j <= n + 1; j++) {
                int t = 0; ++cnt[0];
                for (int k = 1; k <= n; k++) t += b[i][k] + b[n][k] - b[j - 1][k], p[k] = p[k - 1] + b[j - 1][k] - b[i][k], ++cnt[p[k]];
                for (int k = 1, c = 0; k <= n; k++) {
                    if (b[j - 1][k] - b[i][k]) c = 0;
                    else c++, ans[s] += c;
                }
                for (int k = 1; k <= dcnt; k++) {
                    if (d[k] < t) continue;
                    int x = t; ans[d[k]] += cnt[d[k] - t];
                    for (int r = n; r; r--) {
                        x += b[j - 1][r] - b[i][r];
                        if (x > d[k]) break;
                        ans[d[k]] += cnt[d[k] - x];
                    }
                }
                for (int k = 0; k <= n; k++) --cnt[p[k]];
            }
        }
        for (int i = 0; i <= s; i++) if (s % (i + 1) == 0 && ans[s / (i + 1)]) cout << i << " " << ans[s / (i + 1)] << "\n";
    }
    return 0;
}

T2

因子操作

观察到最短路不增，每次更新直接枚举入度尝试更新之前点的答案，如果更新了就继续枚举那个点的入度，否则不管。显然这样做的复杂度是 \(\sum dist_xind_x\) 的，其中 \(dist_x\) 为初始最短路，\(ind_x\) 为入度。搜索发现这个东西的和似乎好像不是很会超过 \(10^8\)，反正过了，跑得飞快。

代码

#include <iostream>
#include <vector>
#include <queue>
using namespace std;
int n, w, q;
queue<int> Q;
vector<int> G[50005];
bool vis[50005];
int dist[50005];
int v[50005];
void dfs(int x) {
    dist[x] = 2147483647;
    vis[x] = 1;
    for (int l = 2, r; l <= min(x, w); l = r + 1) {
        r = (x / (x / l));
        G[x / l].emplace_back(x);
        if (!vis[x / l]) dfs(x / l);
        dist[x] = min(dist[x], dist[x / l] + v[x]);
    }
    if (w > x) dist[x] = min(dist[x], v[x]);
}
int main() {
    freopen("yaya.in", "r", stdin);
    freopen("yaya.out", "w", stdout);
    ios::sync_with_stdio(false);
    cin.tie(0);
    cout.tie(0);
    cin >> n >> w >> q;
    for (int i = 1; i <= n; i++) {
        for (int j = i; j <= n; j += i) 
            ++v[j];
    }
    for (int i = 1; i <= n; i++) dfs(i);
    while (q--) {
        int op, x;
        cin >> op >> x;
        if (op == 1) {
            --dist[x]; --v[x];
            Q.push(x);
            while (Q.size()) {
                x = Q.front(); Q.pop();
                for (auto y : G[x]) if (dist[y] > dist[x] + v[y]) dist[y] = dist[x] + v[y], Q.push(y);
            }
        } else cout << dist[x] << "\n";
    }
    return 0;
}

T3

数点

考虑经典 trick，我们把三次方转化为在矩形中任选三个可重点的方案数。于是只需要对每三个点求包含它们的矩形个数的和即可。我们考虑分治，每次统计左边两个点右边一个点的三元组的贡献。贡献分六种，按照最左边点在右边点上面还是下面分两类各 3 种，并由对称性只讨论其中一类。中间点有三种情况：

1. 在中间。从左往右扫，BIT 维护每个点左上点的 \(\sum i(n - a_i + 1)\)，右边直接查即可。

2. 在上面。从右往左扫左边，BIT 维护每个点右上点的 \(\sum (n - a_* + 1)\)，再乘上自己的 \(i\)，右边直接查即可。

3. 在下面。从下往上扫右边点，动态维护这类贡献。当一个左边点从在上面变成在下面，删去它和它右下点的贡献，并加上它和它左上点的贡献即可。

另一类对称，再做一遍即可。对于左边一个右边两个的三元点组，翻转序列重做。注意我们希望整个分治和原来呈镜像对称，于是区间长度为奇数时需要调整中点位置。总复杂度 \(\mathcal{O}(n \log^2n)\)。

代码

#include <iostream>
#include <algorithm>
#include <string.h>
#include <vector>
#define lowbit(x) ((x) & (-(x)))
#define int long long
using namespace std;
const int P = 998244353;
inline void Madd(int &x, int y) { (x += y) >= P ? (x -= P) : 0; }
int n, K;
int a[100005], _a[100005];
namespace _1 {
    int work() {
        int ret = 0;
        for (int i = 1; i <= n; i++) ret += a[i] * (n - a[i] + 1) % P * i % P * (n - i + 1) % P;
        return ret % P;
    }
}
struct BIT {
    int bit[100005];
    vector<pair<int, int> > vec;
    void add(int x, int y, int z = 0, int w = 1) {
        z ? (x = n - x + 1) : 0;
        if (w) vec.emplace_back(x, y);
        for (; x <= 100000; x += lowbit(x)) Madd(bit[x], y);
    }
    int query(int x, int z = 0) {
        int ret = 0;
        for ((z ? x = n - x + 1 : 0); x; x -= lowbit(x)) Madd(ret, bit[x]);
        return ret;
    }
    void Re() { while (vec.size()) add(vec.back().first, P - vec.back().second, 0, 0), vec.pop_back(); }
} bit1, bit2, bit3, bit4;
namespace _2 {
    int work() {
        int ret = 0;
        for (int i = 1; i <= n; i++) {
            ret += (n - i + 1) * (n - a[i] + 1) % P * bit1.query(a[i]) % P;
            ret += (n - i + 1) * a[i] % P * bit2.query(a[i], 1) % P;
            bit1.add(a[i], a[i] * i % P);
            bit2.add(a[i], i * (n - a[i] + 1) % P, 1);
        }
        return 2 * ret % P;
    }
}
namespace _3 {
    int ret;
    // 0 LU, 1 RD, 2 RU, 3 LD
    int val[4][100005], asdf;
    int ar[100005], al[100005];
    void Solve(int l, int r) {
        if (r - l + 1 <= 2) return;
        int mid = (l + r) >> 1;
        if (asdf && ((r - l + 1) & 1)) --mid;
        Solve(l, mid), Solve(mid + 1, r);
        for (int i = l; i <= mid; i++) {
            bit3.add(a[i], val[0][i] = bit1.query(a[i], 1), 1);
            bit1.add(a[i], i * (n - a[i] + 1) % P, 1);
        }
        for (int i = mid; i >= l; i--) {
            bit4.add(a[i], val[2][i] = i * bit2.query(a[i], 1) % P, 1);
            bit2.add(a[i], n - a[i] + 1, 1);
        }
        for (int i = mid + 1; i <= r; i++) {
            Madd(ret, (n - i + 1) * a[i] % P * bit3.query(a[i], 1) % P);
            Madd(ret, (n - i + 1) * a[i] % P * bit4.query(a[i], 1) % P);
        }
        bit1.Re(), bit2.Re(), bit3.Re(), bit4.Re();

        for (int i = l; i <= mid; i++) {
            bit3.add(a[i], val[3][i] = bit1.query(a[i]));
            bit1.add(a[i], i * a[i] % P);
        }
        for (int i = mid; i >= l; i--) {
            bit4.add(a[i], val[1][i] = i * bit2.query(a[i]) % P);
            bit2.add(a[i], a[i]);
        }
        for (int i = mid + 1; i <= r; i++) {
            Madd(ret, (n - i + 1) * (n - a[i] + 1) % P * bit3.query(a[i]) % P);
            Madd(ret, (n - i + 1) * (n - a[i] + 1) % P * bit4.query(a[i]) % P);
        }
        bit1.Re(), bit2.Re(), bit3.Re(), bit4.Re();

        for (int i = l; i <= mid; i++) al[i] = a[i];
        for (int i = mid + 1; i <= r; i++) ar[i] = a[i];
        sort(al + l, al + mid + 1), sort(ar + mid + 1, ar + r + 1);
        int cur = 0;
        for (int i = mid + 1, j = l; i <= r; i++) {
            while (j <= mid && al[j] < ar[i]) {
                Madd(cur, P - val[1][_a[al[j]]] * (n - al[j] + 1) % P);
                Madd(cur, val[0][_a[al[j]]] * al[j] % P);

                Madd(cur, P - val[3][_a[al[j]]] * (n - al[j] + 1) % P);
                Madd(cur, val[2][_a[al[j]]] * al[j] % P);
                ++j;
            }
            Madd(ret, cur * (n - _a[ar[i]] + 1) % P);
        }
        cur = 0;
    }
    int work() {
        Solve(1, n);
        reverse(a + 1, a + n + 1);
        for (int i = 1; i <= n; i++) _a[a[i]] = i;
        asdf = 1;
        Solve(1, n);
        return ret;
    }
}
signed main() {
    freopen("points.in", "r", stdin);
    freopen("points.out", "w", stdout);
    cin >> n >> K;
    for (int i = 1; i <= n; i++) cin >> a[i], _a[a[i]] = i;
    if (K == 1) cout << _1::work() << "\n";
    else if (K == 2) cout << (_2::work() + _1::work()) % P << "\n";
    else cout << (6 * _3::work() + 3 * _2::work() + _1::work()) % P << "\n";
    return 0;
}

T4

有向无环图

首先观察到负数不会被乘，而且 \(|a_i| \le i\)，每个点只会连向更大点。于是就发现一个负数变成正数的最早时间就是其某个后继变成正数的最早时间 \(+1\)。从而可以将每个点的贡献独立开，视为每个点 \(i\) 从 \(d_i\) 开始会将其权值随 DAG 传播。注意到传播的过程可以用矩阵刻画。将每一天第 \(i\) 个点（\(d_i \ge t\)）对每个点的贡献用一个列向量表示，则每过一天相当于给这个贡献左乘一个矩阵 \(T\)。容易发现矩阵 \(T\) 就是原图的邻接矩阵，将每一列的每个元素乘以列号，得到的东西。那么每次询问的答案可以写成 \(\sum\limits_{p = l}^r\sum\limits_{d_i \le k}(T^{k - d_i}\vec{e_i}a_i)_p + \sum\limits_{p = l}^{r} [d_i \le k]a_i\)。后一半每次询问暴力即可。我们考虑前一半。（\(\vec{e_i}\) 表示单位阵的第 \(i\) 列。）

接下来引入特征值与特征向量的概念。对于一个矩阵 \(A\)，其的一对特征值与特征向量是一对满足 \(T\vec{x} = \lambda\vec{x}\) 的 \(\lambda\) 和 \(\vec{x}\)。它们还满足：\(T^d\vec{x} = \lambda^d\vec{x}, T^{-1}\vec{x} = \lambda^{-1}\vec{x}\)。

那么回过来看本题。由于本题的矩阵为上三角矩阵，于是其 \(n\) 个特征值即为其对角线上的所有 \(n\) 个元素。接下来要对每个特征值求出对应的特征向量，只需要将 \((T - \lambda_iE)\vec{x_i} = \vec{0}\) 视为线性方程组解方程即可。由于是上三角矩阵，因此可以直接执行回代，单次复杂度为 \(\mathcal{O}(n^2)\)，做 \(n\) 次的总复杂度是三次方。注意这里 \((T - \lambda_iE)\) 对角线上必有一个 \(0\)。这个东西是一个自由变量，为了不让高斯消元解出来一个平凡解，我们任意给它钦定一个非零值就可以了。其他都正常。

求出 \(n\) 个特征向量之后，我们考虑用 \(n\) 个特征向量来表示每个 \(e_i\)。设 \(e_i = \sum\limits_{j = 1}^nc_{i, j}x_j\)，发现这玩意就是个矩阵乘法，即 \(E = CX\)。那么 \(C = X^{-1}\)。由于特征向量都是线性无关的（我不会证），因此 \(X\) 满秩，高消求逆即可。

接下来再考虑原来的式子。

\(\begin{equation}\begin{split} \sum\limits_{p = l}^r\sum\limits_{d_i \le k}(T^{k - d_i}\vec{e_i}a_i)_p &= \sum\limits_{p = l}^r\sum\limits_{d_i \le k}(T^{k - d_i}\sum\limits_{j = 1}^nc_{i, j}x_ja_i)_p \\ &= \sum\limits_{p = l}^r\sum\limits_{d_i \le k}a_i(\sum\limits_{j = 1}^nc_{i, j}\lambda_j^{k - d_i}x_j)_p \\ &= \sum_{d_i \le k}a_i(\sum_{j = 1}^nc_{i, j}\lambda^{k - d_i})\sum\limits_{p = l}^rx_{j, p} \\ &= \sum\limits_{j = 1}^n\lambda_j^k\sum\limits_{d_i \le k}c_{i, j}\lambda_j^{-d_i}a_i\sum\limits_{p = l}^rx_{j, p} \end{split}\end{equation}\)

于是对每个 \(j\) 开 BIT 维护 \(c_{i, j}\lambda^{-d_i}a_i\)，最后那个区间和用前缀和搞掉，单次询问即为 \(\mathcal{O}(n \log n)\)。对于修改，若其没有影响任何 \(a_i\) 是否为正，则只需要在每个 BIT 上分别更新一下 \(d_i\) 位置。否则先重新 dfs 求出每个 \(d_i\)，然后发现对每个 \(d_i\) 的改变都枚举每个 BIT 并修改很浪费，于是直接重构每棵 BIT 就好了。由于这样的修改总共只有 \(\mathcal{O}(n)\) 个，因此重构的总复杂度为 \(\mathcal{O}(n^3)\)，需要预处理每个 \(\lambda_j\) 的 \(-d_i\) 次方（注意到 \(d_i \le n\)）。于是做完了，总复杂度 \(\mathcal{O}(n^3 + qn\log n)\)。

代码

#include <iostream>
#include <string.h>
#include <vector>
#define int long long
#define lowbit(x) ((x) & (-(x)))
using namespace std;
const int P = 1000000007;
inline void Madd(int &x, int y) { (x += y) >= P ? (x -= P) : 0; }
inline int Msum(int x, int y) { return Madd(x, y), x; }
int qpow(int x, int y = P - 2) {
    int ret = 1;
    while (y) {
        if (y & 1) 
            ret = ret * x % P;
        y >>= 1;
        x = x * x % P;
    }
    return ret;
}
int A[305][605];
void Gauss(int n, int m, bool f = 1) {
    if (f) for (int i = 1; i <= n; i++) {
        if (!(A[i][i] % P)) for (int j = i + 1; j <= n; j++) if (A[j][i] % P) { swap(A[i], A[j]); break; }
        int coe = qpow(A[i][i]);
        for (int j = i + 1; j <= n; j++) for (int k = m; k >= i; k--) Madd(A[j][k], P - A[i][k] * coe % P * A[j][i] % P);
    }
    for (int i = n; i; i--) {
        int t = qpow(A[i][i]);
        for (int j = m; j >= i; j--) A[i][j] = A[i][j] * t % P;
        for (int j = 1; j < i; j++) {
            for (int k = m; k >= i; k--) Madd(A[j][k], P - A[j][i] * A[i][k] % P);
        }
    }
}
int n, q;
int a[305];
int T[305][305];
vector<int> G[305];
int X[305][305], C[305][305];
int pre[305][305];
int d[305];
bool vis[305], N[305];
int dfs(int x) {
    vis[x] = 1; d[x] = (N[x] ? P : 0);
    for (int v : G[x]) {
        if (!vis[v]) d[x] = min(d[x], dfs(v) + 1);
        else d[x] = min(d[x], d[v] + 1);
    }
    return d[x];
}
struct BIT {
    int bit[305];
    void build(int *a) {
        for (int i = 1; i <= n + 1; i++) bit[i] = Msum(bit[i - 1], a[i]);
        for (int i = n + 1; i; i--) bit[i] = Msum(bit[i], P - bit[i - lowbit(i)]);
    }
    void add(int x, int y) { for (++x; x <= n + 1; x += lowbit(x)) Madd(bit[x], y); }
    int query(int x) {
        int ret = 0;
        for (++x; x; x -= lowbit(x)) Madd(ret, bit[x]);
        return ret;
    }
} bit[305];
int ipw[305][305];
int tmp[305];
void Re() {
    memset(vis, 0, sizeof vis);
    for (int i = 1; i <= n; i++) dfs(i);
    for (int j = 1; j <= n; j++) {
        memset(tmp, 0, sizeof tmp);
        for (int i = 1; i <= n; i++) d[i] != P ? Madd(tmp[d[i] + 1], C[i][j] * ipw[j][d[i]] % P * (P + a[i]) % P) : void();
        bit[j].build(tmp);
    }
}
int asdf[305][305];
signed main() {
    freopen("dag.in", "r", stdin);
    freopen("dag.out", "w", stdout);
    ios::sync_with_stdio(false);
    cin.tie(0);
    cout.tie(0);
    cin >> n;
    for (int i = 1; i <= n; i++) {
        ipw[i][0] = 1, ipw[i][1] = qpow(i);
        for (int j = 2; j <= n + 1; j++) ipw[i][j] = ipw[i][j - 1] * ipw[i][1] % P;
    }
    for (int i = 1; i <= n; i++) cin >> a[i], N[i] = (a[i] <= 0);
    for (int i = 1; i <= n; i++) {
        int k, x;
        cin >> k;
        while (k--) cin >> x, T[i][x] = x, G[i].emplace_back(x);
    }
    for (int i = 1; i <= n; i++) T[i][i] = i;
    for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= n; j++) A[j][n + 1] = 0;
        for (int x = 1; x <= n; x++) {
            for (int y = 1; y <= n; y++) 
                A[x][y] = Msum(T[x][y], P - (x == y) * T[i][i]);
        }
        A[i][i] = A[i][n + 1] = 1;
        Gauss(n, n + 1, 0);
        for (int j = 1; j <= n; j++) X[i][j] = A[j][n + 1];
    }
    memset(A, 0, sizeof A);
    for (int i = 1; i <= n; i++) A[i][i + n] = 1;
    for (int x = 1; x <= n; x++) for (int y = 1; y <= n; y++) A[x][y] = X[x][y];
    Gauss(n, n * 2);
    for (int i = 1; i <= n; i++) {
        for (int j = n + 1; j <= n * 2; j++) 
            C[i][j - n] = A[i][j];
    }
    for (int i = 1; i <= n; i++) for (int j = 1; j <= n; j++) pre[i][j] = Msum(pre[i][j - 1], X[i][j]);
    Re();
    cin >> q;
    while (q--) {
        int op, x, y, k;
        cin >> op;
        if (op == 1) {
            cin >> k >> x >> y;
            int ans = 0;
            for (int j = 1; j <= n; j++) Madd(ans, qpow(j, k) * bit[j].query(min(n, k)) % P * Msum(pre[j][y], P - pre[j][x - 1]) % P);
            for (int i = x; i <= y; i++) if (d[i] > k) Madd(ans, P + a[i]);
            cout << ans << "\n";
        } else {
            cin >> x >> y; a[x] += y;
            if (N[x] && a[x] > 0) N[x] = 0, Re();
            else for (int j = 1; j <= n; j++) d[x] != P ? bit[j].add(d[x], ipw[j][d[x]] * C[x][j] % P * y % P) : void();
        }
    }
    return 0;
}

---

第四次因为 memset 没加 #include <string.h> CE，已经在 snippets 里加上了。以后再也不会了。

T4，式子中若有矩阵乘向量，可考虑使用特征值与特征向量，尝试用特征向量的线性组合表示被乘向量来消去矩阵乘法在式子中的存在。

DAG 中注意是否保证 \(1\) 能走到所有点！！！