B3646
一扶苏一
·
·
题解
数列前缀和 3
这是一道洛谷夏令营的作业题,出这道题的目的是强调非交换群的区间和必须注意左右乘的区别。
考虑维护矩阵的前缀积:s_i 表示前 i 个矩阵的乘积。 注意对于一个 [l, r] 的询问,\prod\limits_{i = l}^r a_i \neq s_r \times s_{l-1}^{-1} 。例如:CD \neq ABCD \times (AB)^{-1}。
但是注意到 (AB)^{-1} = B^{-1} A^{-1}。我们只要将 (AB)^{-1} 左乘上 ABCD,就有 (AB)^{-1} \times ABCD = (B^{-1}(A^{-1}A)B)CD = CD。
所以我们有 \prod\limits_{i = l}^r a_i = s_{l - 1}^{-1} \times s_r。维护前缀积,做一个矩阵求逆即可。时间复杂度 O((n+q) k^3)。
矩阵求逆可以用高斯消元法,但因为 k 只有 2 和 3,所以也可以手搓一下伴随矩阵,下面介绍一下伴随矩阵法求逆矩阵。
对 n 阶矩阵 A=(a_{ij})_{n \times n},划去第 i 行第 j 列后剩余的 (n-1) 阶行列式的值称为元素 a_{i,j} 的余子式,记为 M_{i, j}。记 A_{i, j}=(-1)^{i+j}M_{i,j} 为 a_{i, j} 的代数余子式。
对于二阶矩阵 $A = \begin{pmatrix} a & b \\ c & d\end{pmatrix}$,直接套用上述结论得到 $A^{-1} = \frac{1}{ad - bc}\begin{pmatrix} d & -b \\ -c & a\end{pmatrix}$;对于三阶矩阵,其行列式只有六项,所有的代数余子式均为二阶行列式,都可以轻松算出。
```cpp
#include <array>
#include <iostream>
#include <algorithm>
const int p = 1145141;
typedef long long int ll;
int n, k, q;
std::array<int, p> inv;
struct Matrix {
ll A[3][3];
inline Matrix operator*(const Matrix& o) const {
Matrix ret;
for (int i = 0; i < k; ++i) {
for (int j = 0; j < k; ++j) {
ret.A[i][j] = 0;
for(int h = 0; h < k; ++h) {
ret.A[i][j] += A[i][h] * o.A[h][j];
}
ret.A[i][j] %= p;
}
}
return ret;
}
inline int Det2(ll a, ll b, ll c, ll d) {
int ret = (a * d - b * c) % p;
if (ret < 0) ret += p;
return ret;
}
inline int Det() {
if (k == 2) { return Det2(A[0][0], A[0][1], A[1][0], A[1][1]); }
else {
ll ret = 0;
(ret += A[0][0] * A[1][1] * A[2][2]) %= p;
(ret += A[0][1] * A[1][2] * A[2][0]) %= p;
(ret += A[0][2] * A[1][0] * A[2][1]) %= p;
(ret -= A[0][2] * A[1][1] * A[2][0]) %= p;
(ret -= A[0][0] * A[1][2] * A[2][1]) %= p;
(ret -= A[0][1] * A[1][0] * A[2][2]) %= p;
(ret += p) %= p;
return ret;
}
}
inline Matrix operator~() {
ll d = inv[Det()];
Matrix ret;
if (k == 2) {
ret.A[0][0] = A[1][1] * d % p;
ret.A[1][0] = (p - A[1][0]) * d % p;
ret.A[1][1] = A[0][0] * d % p;
ret.A[0][1] = (p - A[0][1]) * d % p;
} else {
ret.A[0][0] = Det2(A[1][1], A[1][2], A[2][1], A[2][2]) * d % p;
ret.A[1][0] = Det2(A[1][0], A[1][2], A[2][0], A[2][2]) * d % p;
ret.A[2][0] = Det2(A[1][0], A[1][1], A[2][0], A[2][1]) * d % p;
ret.A[0][1] = Det2(A[0][1], A[0][2], A[2][1], A[2][2]) * d % p;
ret.A[1][1] = Det2(A[0][0], A[0][2], A[2][0], A[2][2]) * d % p;
ret.A[2][1] = Det2(A[0][0], A[0][1], A[2][0], A[2][1]) * d % p;
ret.A[0][2] = Det2(A[0][1], A[0][2], A[1][1], A[1][2]) * d % p;
ret.A[1][2] = Det2(A[0][0], A[0][2], A[1][0], A[1][2]) * d % p;
ret.A[2][2] = Det2(A[0][0], A[0][1], A[1][0], A[1][1]) * d % p;
for (int i = 0; i < 3; ++i) {
for (int j = 0; j < 3; ++j) if ((i + j) & 1) {
ret.A[i][j] = p - ret.A[i][j];
}
}
}
return ret;
}
};
Matrix a[p];
void get_inv(const int x) {
inv[1] = 1;
for (int i = 2; i < x; ++i) inv[i] = 1ll * (p - p / i) * inv[p % i] % p;
}
int main() {
std::ios::sync_with_stdio(false);
std::cin.tie(nullptr);
get_inv(p);
std::cin >> n >> k >> q;
for (int i = 0; i < k; ++i) a[0].A[i][i] = 1;
for (int i = 1; i <= n; ++i) {
for (int j = 0; j < k; ++j) {
for (int h = 0; h < k; ++h) {
std::cin >> a[i].A[j][h];
}
}
a[i] = a[i - 1] * a[i];
}
int ans = 0;
for (int l, r; q; --q) {
std::cin >> l >> r;
Matrix mul = (~a[l - 1]) * a[r];
for (int i = 0; i < k; ++i) {
for (int j = 0; j < k; ++j) {
ans ^= mul.A[i][j];
}
}
}
std::cout << ans << std::endl;
}
```