矩阵 - 题目列表

矩阵

GitHub跳转原题关系图 Canvas 关系图返回列表

在 GF(2) 中提取 Krylov 状态序列的最小递推，用多项式快速幂求 A^k b。

OJ: shumeng

题目 ID: CSP201512E

难度:省选/NOI-

标签:线性代数GF(2)矩阵快速幂

日期: 2026-07-31 16:21

形式化题目

在 GF(2) 上（加法为异或、乘法为与）给定一个可逆的 $m \times m$ 矩阵 $A$ 和初始向量 $b$ 。回答 $n$ 组询问，每组给定非负整数 $k$ ，输出 $A^k b$ （其中 $A^0 b = b$ ）。

思路

直接对每个询问做矩阵快速幂需要 $O(m^3 \log k)$ ，而 $m$ 最大 $1000$ ，不可接受。本题的关键是：只针对给定的一个向量 $b$ 建立 Krylov 序列，而不是对整张矩阵做幂。

Krylov 序列与最小递推

考虑向量序列：

b,\; Ab,\; A^2b,\; \ldots

这些向量都在 $m$ 维空间中，最多 $m$ 个线性无关，因此必存在最小的正整数 $d$ 使得 $A^d b$ 可被 $A^0b, \ldots, A^{d-1}b$ 线性表出：

A^d b = \sum_{i=0}^{d-1} c_i A^i b

这个线性关系对应一个递推多项式 $P(x) = x^d - \sum c_i x^i$ ，满足 $P(A)b = 0$ 。

高斯消元找递推

逐项生成状态 $state[i] = A^i b$ ，对每个新向量做“按最高位”的位集高斯消元：

若线性无关，加入基，记录它对应的 $A$ 的幂次表达式；
若消到零，说明找到了首个线性关系，得到 $d$ 和系数 $c_i$ ，停止。

多项式快速幂

对每个询问，求 $x^k \bmod P(x)$ 的多项式系数（在 GF(2) 上做多项式乘法并对 $P(x)$ 取模），答案就是：

A^k b = \sum_i [x^i]\big(x^k \bmod P(x)\big) \cdot A^i b

位集让每个多项式乘法只需 $O(d^2/64)$ 的时间。

代码

cpp

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113/**
 * Author by Rainboy blog: https://rainboylv.com github: https://github.com/rainboylvx
 * rbook: -> https://rbook.roj.ac.cn  https://rbook2.roj.ac.cn
 * rainboy的学习导航网站: https://idx.roj.ac.cn
 * create_at: 2026-07-31 16:21
 * update_at: 2026-08-17 23:01
 */
#include <bits/stdc++.h>
using namespace std;

const int MAXM = 1005; // 向量/矩阵最大维数
const int MAXP = 2010; // 递推多项式最高次数（2*m 量级）

int m;
bitset<MAXM> matrix_row[MAXM];    // 矩阵 A 的每一行，用于把向量左乘 A
bitset<MAXM> linear_basis[MAXM];  // 高斯消元后的向量基（按最高位归位）
bitset<MAXM> linear_expression[MAXM]; // 基向量对应的多项式系数（列）
bitset<MAXM> state[MAXM];         // state[i] = A^i * b
bitset<MAXP> relation;            // 最小递推多项式 P(x)，满足 x^d = sum relation[i]*x^i

// 计算 A * vector_value（GF(2) 下，行向量点积取异或）。
bitset<MAXM> multiply_vector(const bitset<MAXM> &vector_value) {
    bitset<MAXM> result;
    for (int i = 0; i < m; i++) {
        result[i] = (matrix_row[i] & vector_value).count() % 2;
    }
    return result;
}

// GF(2) 多项式乘法并对 P(x) 取模（次数为 degree）。
bitset<MAXP> multiply_polynomial(const bitset<MAXP> &left, const bitset<MAXP> &right,
                                 int degree) {
    bitset<MAXP> result;
    // 先做普通多项式乘法，再对次数 >= degree 的项用 x^degree = P(x) 归约。
    for (int i = 0; i < degree; i++) {
        if (left[i]) result ^= right << i;
    }
    for (int i = 2 * degree - 2; i >= degree; i--) {
        if (result[i]) result ^= relation << (i - degree);
    }
    return result;
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    cin >> m;
    for (int i = 0; i < m; i++) {
        string line;
        cin >> line;
        for (int j = 0; j < m; j++) matrix_row[i][j] = line[j] == '1';
    }
    string initial;
    cin >> initial;
    bitset<MAXM> current; // 当前状态向量 b
    for (int i = 0; i < m; i++) current[i] = initial[i] == '1';

    // 逐项生成 Krylov 序列 b, Ab, A^2b, ...，用高斯消元找第一个线性相关。
    int degree = 0;
    for (int step = 0; step <= m; step++) {
        bitset<MAXM> value = current;
        bitset<MAXM> expression;
        expression[step] = 1; // 初始表示 A^step * b = 1 * (A^step * b)
        // 按最高位从高到低消元。
        for (int bit = m - 1; bit >= 0; bit--) {
            if (!value[bit]) continue;
            if (linear_basis[bit].none()) {
                // 当前向量线性无关，加入基并记录它对应的 A 的幂次。
                linear_basis[bit] = value;
                linear_expression[bit] = expression;
                state[step] = current;
                degree = step + 1;
                break;
            }
            value ^= linear_basis[bit];
            expression ^= linear_expression[bit];
        }
        if (value.none()) {
            // 找到首个线性关系：A^step * b = sum expression[i] * A^i * b。
            degree = step;
            for (int i = 0; i < degree; i++) relation[i] = expression[i];
            relation[degree] = 1; // P(x) = x^degree - sum expression[i]*x^i
            break;
        }
        current = multiply_vector(current);
    }

    int q;
    cin >> q;
    while (q--) {
        long long k;
        cin >> k;
        // 用多项式快速幂计算 x^k mod P(x)，结果系数决定 state 的线性组合。
        bitset<MAXP> result, power;
        result[0] = 1;
        if (degree == 1) power[0] = relation[0]; // 一次递推的特例：x = relation[0]
        else power[1] = 1;                        // 否则 x^1
        while (k > 0) {
            if (k & 1) result = multiply_polynomial(result, power, degree);
            power = multiply_polynomial(power, power, degree);
            k >>= 1;
        }
        // 答案 = sum result[i] * state[i]。
        bitset<MAXM> answer;
        for (int i = 0; i < degree; i++) {
            if (result[i]) answer ^= state[i];
        }
        for (int i = 0; i < m; i++) cout << answer[i];
        cout << '\n';
    }
    return 0;
}

复杂度

设 $d$ 为递推阶数， $d \leqslant m$ 。

预处理：生成 $m$ 个状态并消元，位集操作下约为 $O(m^3/64)$ 。
每次询问：多项式快速幂， $O(d^2 \log k / 64)$ 。

总结

不需要显式计算 $m \times m$ 矩阵的幂。只针对给定初始向量建立 Krylov 序列，维数最多为 $m$ ，用位集高斯消元提取最小递推多项式后，任意高次幂都能通过多项式取模快速求出。这是“线性递推 + 快速幂”思想在向量空间中的推广。

时间限制： 1.0 秒

空间限制： 256 MB

问题描述

创造一个世界只需要定义一个初状态和状态转移规则。

宏观世界的物体运动规律始终跟物体当前的状态有关，也就是说只要知道物体足够多的状态信息，例如位置、速度等，我们就能知道物体之后任意时刻的状态。

现在小 M 创造了一个简化的世界。

这个世界中，时间是离散的，物理规律是线性的：世界的初始状态可以用一个 $m$ 维向量 $b^{(0)}$ 表示，状态的转移方式用 $m \times m$ 的矩阵 $A$ 表示。

若已知这个世界当前的状态是 $b$ ，那么下一时刻就等于 $b$ 左乘状态转移矩阵 $A$ ，即 $Ab$ 。

这个世界中，物体的状态也是离散的，也就是说可以用整数表示。再进一步，整数都可以用二进制编码拆分为有限位 $0$ 和 $1$ 。因此，这里的矩阵 $A$ 和向量 $b$ 的每个元素都是 $0$ 或 $1$ ，矩阵乘法中的加法运算视为异或运算（ $\mathrm{xor}$ ），乘法运算视为与运算（ $\mathrm{and}$ ）。

具体地，设矩阵 $A$ 第 $i$ 行第 $j$ 列的元素为 $a_{i, j}$ ，向量 $b$ 的第 $i$ 个元素为 $b_i$ 。那么乘法 $Ab$ 所得的第 $k$ 个元素为

(a_{k,1}\text{ and }b_1)\text{ xor }(a_{k,2}\text{ and }b_2)\text{ xor } \cdots\text{ xor }(a_{k,m}\text{ and }b_m)

矩阵和矩阵的乘法也有类似的表达。

小 M 发现，这样的矩阵运算也有乘法结合律，例如有 $A(Ab) = (AA)b = A^2b$ 。

为了保证自己创造的世界维度不轻易下降，小 M 保证了矩阵 $A$ 可逆，也就是说存在一个矩阵 $A^{-1}$ ，使得对任意向量 $d$ ，都有 $A^{-1}Ad=d$ 。

小 M 想了解自己创造的世界是否合理，他希望知道这个世界在不同时刻的状态。

具体地，小 M 有 $n$ 组询问，每组询问会给出一个非负整数 $k$ ，小 M 希望你帮他求出 $A^kb$ 。

输入格式

从标准输入读入数据。

输入第一行包含一个整数 $m$ ，表示矩阵和向量的规模。

接下来 $m$ 行，每行包含一个长度为 $m$ 的 $01$ 串，表示矩阵 $A$ 。

接下来一行，包含一个长度为 $m$ 的 $01$ 串，表示初始向量 $b^{(0)}$ 。（ $b^{(0)}$ 是列向量，这里表示它的转置）

注意： $01$ 串两个相邻的数字之间均没有空格。

接下来一行，包含一个正整数 $n$ ，表示询问的个数。

最后 $n$ 行，每行包含一个非负整数 $k$ ，表示询问 $A^kb(0)$ 。

注意： $k$ 可能为 $0$ ，此时是求 $A^0b^{(0)} = b^{(0)}$ 。

输出格式

输出到标准输出。

输出 $n$ 行，每行包含一个 $01$ 串，表示对应询问中 $A^kb^{(0)}$ 的结果。

注意： $01$ 串两个相邻的数字之间不要输出空格。

input1

output1

评测用例规模与约定

测试点编号	$m=$	$n=$	$k\le $
1	$10$	$100$	$10^3$
2	$10$	$100$	$10^4$
3	$30$	$100$	$10^5$
4	$180$	$100$	$10^5$
5	$10$	$100$	$10^9$
6	$30$	$100$	$10^9$
7	$180$	$100$	$10^9$
8	$600$	$100$	$10^9$
9	$800$	$100$	$10^9$
10	$1000$	$100$	$10^9$

矩阵

目录

形式化题目

思路

Krylov 序列与最小递推

高斯消元找递推

多项式快速幂

代码

复杂度

总结