高斯消元：主元选择、自由变量与方程组求解

一、引入：我们为什么敢对方程做加减？

有两种商品，买两个甲和一个乙花 8 元，买一个甲和两个乙花 7 元。甲、乙各多少钱？

你大概会先把其中一个式子乘二，再相减，让一个未知数“消失”。

这就是消元。到了程序里，我们仍然做同一件事，只不过未知数可能有几十个，方程也不一定和未知数一样多。我们要把“挑一个未知数消掉”整理成固定流程，还要判断：答案只有一个，还是根本没有答案，或者有无穷多个答案。

设甲单价为 $x$ ，乙单价为 $y$ ，问题变成：

\begin{cases}2x+y=8,\\x+2y=7.\end{cases}

第二行乘二再减第一行，得到 $3y=6$ ，于是 $y=2$ ，再代回得到 $x=3$ 。我们没有猜价格，而是通过新等式把原来互相牵扯的两个量拆开了。

在消元过程中，我们允许三种基本操作：

交换两行。
一行乘上非零数。
把一行的某个倍数加到另一行。

物理意义：为什么这样做不会改变解集？交换只是改顺序；非零倍乘可以除回来；加上去的那一行也可以再减回来。每一步都能撤销，所以旧答案不会消失，新答案也不会冒出来。

这里有两个绝不能省的动作：第一，等号右边也要一起变化，不能只改左边系数。第二，一行乘零会把信息抹掉，这不属于允许的非零倍乘。别把一行式子直接清空，当成“消掉一个方程”。

为了方便程序处理，我们把未知数的名字暂时省略，只存系数和右侧常数，这就得到了增广矩阵：

\left[\begin{array}{cc|c}2&1&8\\1&2&7\end{array}\right]

最后一列并不是多出来的未知数，它是等号右边的常数。如果程序处理

n

个未知数，那每一行就要存

n+1

个数。

二、核心思想：先找一个可靠的“主角”（主元）

准备处理 $x$ 这一列时，就先找一行，借它表达 $x$ 与其余量的关系。这一行叫主元行，被选中的 $x$ 系数叫主元（Pivot）。

在上例中选第一行，除以 2，得到 $x+0.5y=4$ 。再让第二行减去这一行，变成 $1.5y=3$ 。第二列也选到主元以后，把第二行除以 1.5，得到 $y=2$ ；再从第一行消去半个 $y$ ，得到 $x=3$ 。

最后两行就像一张答案表：每个未知数各占一行，自己的系数为 1，其他未知数的系数为 0。这种把主元列“上下都消干净”的写法，常称为高斯—约旦消元。

高斯消元的四个增广矩阵状态：从2x+y=8、x+2y=7出发，依次归一和消元，得到x=3、y=2。

1. 为什么还要交换行？

假如第一行是 $0x+y=2$ ，第二行是 $x+y=5$ 。如果照着第一行就除以主元，程序会发生“除以零”的惨剧。我们只需要把两行交换，先使用第二行处理 $x$ ，问题就顺理成章地解决了。

所以，程序不能死板地只盯着当前位置。处理某一列时，要在尚未使用的行中寻找可用主元。找到以后，换到当前待处理行，再进行消元。

2. 降维打击：不只避开零，还要尽量避开很小的数

考虑 $0.00000001x+y=1$ 和 $x+y=2$ 。如果执意用第一行的很小系数作主元，除完后就会出现极大的中间数。后面再做两个大数的减法，原本很小的舍入误差会被无限放大。

因此，我们在当前列里，要选择剩余各行中绝对值最大的系数作为主元，再把它交换到前面。这就叫部分选主元。“部分”是指只在当前列的剩余行中选，不用连未知数的列也一起乱换。它通常能让计算更加稳定。

三、灵魂数组：`where`（打破行与列的绑定）

在最简单的二元题里，处理第一列用第一行，处理第二列用第二行，好像只需要一个下标。但一般输入有 $m$ 个方程、 $n$ 个未知数，还可能遇到整列都没有主元的情况。

所以，我们必须把两个量分开：

col：表示正在尝试消去哪个未知数。
row：表示已经找到了多少个主元（当前正在安排第几行）。

物理推导：每次都向右换列考察未知数（col++），但只有真的找到主元时，才让 row 向下走一格（row++）。

比如只有 $y+z=3$ 、 $2y+2z=6$ 两条方程，却有 $x,y,z$ 三个未知数。第一列的 $x$ 系数全为零，跳过这一列，但不能顺便浪费掉第一行。下一列用第一行处理 $y$ ，第二行被消成全零，最后只找到一个主元。

我们需要一个灵魂数组 where[col]，它的物理意义是：记录是哪一个行号负责了这个未知数（列）。如果没有主元，就记为 -1。这样无论跳过几列，最终读答案时，都不会错把第几行直接当成第几个未知数。

注：找到的主元个数，也就是最后成功推进的 row，在数学上叫系数矩阵的秩。它代表有多少条真正独立的约束。

四、终局判定：三种结果，按什么顺序判断？

先别急着数主元，我们应该先找矛盾。

两条式子 $x+y=3$ 和 $2x+2y=7$ ，第二条减去第一条的两倍后，变成 $0=1$ 。这根本不可能成立，这就意味着整个方程组无解。

如果得到的是 $0=0$ ，那就完全不同了。删去它不会影响其他条件，但它也不能当作一个新约束。

1. 判定逻辑拆解

第一步（判矛盾）：有没有一行，所有未知数系数都接近 0，但等号右端常数明显不为 0？如果有，果断判定无解 (NONE)。即使其他行能解出一些变量，也不能忽略这一行的绝对冲突。
第二步（判自由变量）：如果没有矛盾，再去数主元。如果主元个数（row）小于未知数个数（ $n$ ），说明有些变量没有被约束，得到无穷多解 (INF)。
第三步（唯一解）：主元个数等于未知数个数，每个未知数都被独立约束住，得到唯一解 (UNIQUE)。

几何上的反直觉：假设有两个二元方程 $x+y=1$ 和 $x+y=2$ 。在这个增广矩阵里，你只能找到 1 个非零主元。如果我们跳过判矛盾，直接算 $n - \text{row}$ ，就会得出 $2 - 1 = 1$ ，误以为“存在 1 个自由变量，有无穷多解”。但在几何上，这是两条平行的直线，它们根本没有任何交点。既然连一个立足点都没有，又怎么能顺着直线的方向“自由”滑动呢？

所以，只有在方程组一致（不存在 $0 = c$ 且 $c \neq 0$ 的矛盾）的情况下，计算自由变量才是有意义的。只要扫描到任何一行出现了绝对冲突，整个系统就被判定失效，直接无解。

2. 自由变量是不是就不能往下算了？

当然不是。比如消完后剩下两条约束：第一个未知数加第二个再加第三个等于 4，第二个未知数与第三个相等。你可以把第三个未知数设成自由实数 $t$ ，其他变量也跟着用 $t$ 表达出来。不过在本课的基础模板中，只需要报告无穷多解，不需要输出参数表达式。遇到自由变量时直接返回分类即可。

五、浮点警报：EPS 是尺子的刻度，不是数学中的零

计算机里的小数通常不能精确保存。例如一个理论上应该抵消为零的结果，算出来可能留下一个极小的尾巴（比如 $10^{-16}$ ）。如果直接用 ==0，程序就会把舍入误差当作新的约束。

我们设一个容差常数 EPS = 1e-10：绝对值小于它时，按接近零处理。

要把浮点数当成尺子上的刻度来看待：一条方程两边同时乘上一万，数学上的解没有变，但这把尺子的尺度就变了。固定的绝对容差对不同尺度的意义并不相同。若实际题目系数跨度极大，只靠定死的 EPS 也会出问题。课堂测试中，我们采用绝对值不超过 $10^4$ 、尺度适中的数据。

六、核心模板代码实现

以下代码适用于 $m$ 个方程、 $n$ 个未知数的情况。

输入约定：第一行 m n（ $1\le m,n\le 100$ ）。之后每行给出 $n$ 个系数和一个右端常数。
输出约定：无解输出 NONE；无穷多解输出 INF；唯一解先输出 UNIQUE，再按未知数顺序输出一行数值，保留十位小数。
精度控制：使用 long double 防止中途溢出或掉精度，并在最终输出时处理负零。

实战输入示例

输出结果

UNIQUE
3.0000000000 2.0000000000

#include<bits/stdc++.h>
using namespace std;

// 高斯消元强烈建议使用 long double 和 EPS
using ld = long double;
const ld EPS = 1e-10L;
const int N = 105;

ld a[N][N];
int where[N];
int m, n;

void solve(){
    cin >> m >> n;
    for(int i = 0; i < m; i++){
        for(int j = 0; j <= n; j++){
            cin >> a[i][j];
        }
    }
    
    // 初始化灵魂数组
    memset(where, -1, sizeof(where));
    
    int row = 0; // row 代表已经找到的主元个数，也就是正在安排的行
    for(int col = 0; col < n && row < m; col++){
        // 1. 部分选主元：在当前列的剩余行中，找绝对值最大的系数
        int p = row;
        for(int i = row + 1; i < m; i++){
            if(fabsl(a[i][col]) > fabsl(a[p][col])) p = i;
        }
        
        // 如果这一列全是 0，说明没有主元，跳过，row 不动，去下一列
        if(fabsl(a[p][col]) < EPS) continue;
        
        // 2. 将主元行交换上来 (只换从 col 到 n 的部分即可)
        for(int j = col; j <= n; j++) swap(a[p][j], a[row][j]);
        
        // 3. 归一化：让主元系数变成 1
        ld pivot = a[row][col];
        for(int j = col; j <= n; j++) a[row][j] /= pivot;
        a[row][col] = 1; 
        
        // 4. 消元：用当前行去消掉其他所有行的当前列
        for(int i = 0; i < m; i++){
            if(i == row) continue;
            ld t = a[i][col];
            for(int j = col + 1; j <= n; j++){
                a[i][j] -= t * a[row][j];
            }
            a[i][col] = 0; // 显式置 0，防止微小精度残留
        }
        
        // 记录：第 col 个未知数，是由第 row 行负责出来的
        where[col] = row;
        row++; // 主元找到，行指针往下走
    }
    
    // 第一步判定：找矛盾 (无解)
    for(int i = 0; i < m; i++){
        bool zero = true;
        for(int j = 0; j < n; j++){
            if(fabsl(a[i][j]) >= EPS) zero = false;
        }
        // 左边系数全是 0，但右边常数不为 0，出现 0 = 1 的矛盾
        if(zero && fabsl(a[i][n]) >= EPS){
            cout << "NONE\n";
            return;
        }
    }
    
    // 第二步判定：找自由变量 (无穷多解)
    if(row < n){
        cout << "INF\n";
        return;
    }
    
    // 第三步判定：唯一解输出
    cout << "UNIQUE\n" << fixed << setprecision(10);
    for(int j = 0; j < n; j++){
        ld x = a[where[j]][n];
        // 避免输出 -0.0000000000
        if(fabsl(x) < 5e-11L) x = 0; 
        if(j > 0) cout << " ";
        cout << x;
    }
    cout << '\n';
}

signed main(){
    // 基础优化
    ios::sync_with_stdio(0), cin.tie(0);
    solve();
    return 0;
}

复杂度分析：存储空间为 $O(mn)$ 。处理每个主元时要遍历行与剩余列，主元最多有 $\min(m,n)$ 个，因此总时间复杂度为 $O(mn\min(m,n))$ 。如果 $m$ 和 $n$ 相等，也就是方阵时，这就是我们熟悉的 $O(N^3)$ 。

七、模质数域的高斯消元与逆元（选学）

先修提醒：阅读本节需要掌握基础的取模运算，以及如何通过费马小定理或扩展欧几里得求得“乘法逆元”。

问题场景：如果我们不是求实数解，而是要求所有运算都在模某个质数 $p$ （例如 $10^9+7$ 或 $998244353$ ）的规则下进行，比如 $(2x + 3y) \equiv 5 \pmod 7$ 。

逻辑切换：把除法变成乘法：在实数域里，把主元系数变成 1 的操作是“除以主元”。但在模质数意义下没有直接的实数除法，我们要把它替换为：乘以主元在模 $p$ 下的乘法逆元（记作 $a^{-1}$ ）。

这种“模域消元”有非常省心的特点：彻底告别精度误差，并且选主元只要不为 0 即可，不需要找绝对值最大的。但从实数域模板改造过来，需要做全套的整数化修改，不能只改归一化：

类型与零判断：撤掉 long double 和 EPS，数组用 long long，判零真刀真枪写 == 0。
输入预处理：读入系数时，可能遇到负数或超大数，需统一转成正余数：val = (val % mod + mod) % mod。
归一化（除变乘）：求出主元的逆元，当前行同乘逆元并取模。
消元减法：普通减法会产生负数，必须改成安全的模减法形式。
收尾与输出：系数全为 0、右端却不为 0，仍然判无解；无矛盾时，若有 n-row 个自由变量，解数是 $mod^{n-row}$ ，不是无穷多。输出标签按题意改，唯一解直接输出整数，不再沿用实数版的小数格式。

关键代码变动清单：模数在代码中记为 mod，不要和原模板的主元行号 p 混用。快速幂 qpow(x,mod-2) 复用《数学基础》·组合数 $O(1)$ 查询中的实现，并把 mod 设为当前质数。以下片段替换原主元循环中的选主元、归一化和消元部分：

// 1. 找一个非零主元，换到当前行
int pos = row;
while(pos < m && a[pos][col] == 0) pos++;
if(pos == m) continue;
for(int j = col; j <= n; j++) swap(a[pos][j], a[row][j]);

// 2. 归一化：除法变乘法
long long inv_pivot = qpow(a[row][col], mod - 2);
for(int j = col; j <= n; j++){
    a[row][j] = (a[row][j] * inv_pivot) % mod;
}

// 3. 消元减法：注意用 + mod 抵消负数
for(int i = 0; i < m; i++){
    if(i == row) continue;
    long long t = a[i][col];
    for(int j = col + 1; j <= n; j++){
        // 模域下的减法安全写法
        a[i][j] = (a[i][j] - t * a[row][j] % mod + mod) % mod;
    }
    a[i][col] = 0;
}

八、异或方程组与 GF(2) 的极简加速（选学）

问题场景：经典算法游戏“关灯问题”：每个开关不仅控制自己，还会触发相邻的灯。开关只有“按”或“不按”两种动作，灯也只有“亮（1）”或“灭（0）”两种状态。问怎么按开关能把所有灯熄灭？

在这个系统里，所有系数和状态都是 0 和 1。

加减法变成了异或运算（ $1 \oplus 1 = 0, 1 \oplus 0 = 1$ ），即不进位的加法。
乘法变成了与运算（ $1 \times 1 = 1, 1 \times 0 = 0$ ）。这种只由 0 和 1 构成且运算规则契合位运算的体系，叫作 GF(2) 域（二阶伽罗瓦域）。

位运算加速操作：bitset 压位 既然矩阵每一行全是 0 和 1，在 C++ 里，我们不用普通的 int 二维数组去存它，而是直接用 std::bitset。它可以把 64 个位压入一块连续内存。

对 bitset 的下标和异或操作还不熟，可以先回看《bitset 专题》中的基本操作。

当我们用主元行去消掉其他行时：原本我们要写一整个 for 循环挨个对齐相减；现在只需要一行代码：a[i] ^= a[row];！底层会借助位运算，以 64 位（或 32 位）为一批次瞬间处理完毕，将高斯消元的核心时间复杂度从 $O(N^3)$ 强行除以 64，实现常数级的提速。

解的物理意义修正：如果在 GF(2) 域中发现了自由变量（即找到的主元个数 $\text{row} < n$ ），它意味着有多组解。但与实数域的“无穷多解”不同，二元域的状态只有 0 和 1。 $n - \text{row}$ 个自由变量，每个变量有 2 种选择，因此解的总数是精确的 $2^{n - \text{row}}$ 个。为了标签不引起误解，遇到多解时我们可以输出 MULTIPLE。

完整模板：GF(2) 异或消元

💡【实战输入示例】 第一行是方程数 $m$ 和未知数 $n$ （约束 $1 \le m, n \le 100$ ）。接下来的 $m$ 行，前 $n$ 个数是开关对某盏灯的影响系数（0 或 1），最后第 $n+1$ 个数是目标常数。

输出结果

UNIQUE
1 0 0

仅按下 1 号开关，可满足三盏灯的所有触发条件。

#include<bits/stdc++.h>
using namespace std;

const int N = 105;
// 每一行是一个 bitset，下标 0~n-1 存未知数系数，下标 n 存等式右侧常数
bitset<N> a[N];
int where[N];
int m, n;

void solve(){
    cin >> m >> n;
    for(int i = 0; i < m; i++){
        for(int j = 0; j <= n; j++){
            int val;
            cin >> val;
            a[i][j] = val;
        }
    }
    
    memset(where, -1, sizeof(where));
    int row = 0;
    
    for(int col = 0; col < n && row < m; col++){
        // 1. 找主元：系数只能是 0 或 1，遇到 1 就可以抓来当主元
        int p = row;
        while(p < m && a[p][col] == 0) p++;
        
        // 若当前列全是 0，说明没有主元，跳过处理下一列
        if(p == m) continue;
        
        // 2. 交换整行
        swap(a[row], a[p]);
        
        // 3. 极简消元：不需要乘除归一，直接用当前主元行去异或其余目标行
        for(int i = 0; i < m; i++){
            if(i != row && a[i][col] == 1){
                a[i] ^= a[row]; // 常数暴降的批量处理
            }
        }
        
        where[col] = row;
        row++;
    }
    
    // 判定：是否有 0 = 1 的矛盾？
    for(int i = 0; i < m; i++){
        // 如果左边未知数的系数全被异或成了 0，唯独常数项（第 n 位）是 1，即出现矛盾
        // bitset 的 count() 函数统计整行里 1 的个数
        if(a[i].count() == 1 && a[i][n] == 1){
            cout << "NONE\n";
            return;
        }
    }
    
    // 判定：是否有自由变量
    if(row < n){
        // 在 GF(2) 域下，有 2^(n-row) 个解，而不是实数域的无穷多解
        cout << "MULTIPLE\n";
        return;
    }
    
    // 输出唯一解
    cout << "UNIQUE\n";
    for(int j = 0; j < n; j++){
        cout << a[where[j]][n] << (j == n - 1 ? "" : " ");
    }
    cout << '\n';
}

signed main(){
    ios::sync_with_stdio(0), cin.tie(0);
    solve();
    return 0;
}

九、渐进式实战练习

要彻底掌握这段代码，不要只跑通过的样例，请手动把三种结果都“捏”出来遇到一次。

这里用的是本讲自己的输入输出约定。拿去做在线模板题时，先对照题面调整方程数的读法，以及无解、多解和小数的输出格式。

第一题：检查无解 输入第一行 2 2，两行分别为 1 1 3、2 2 7，输出应为 NONE。

训练指引：单步调试，请指出哪一步出现了矛盾行，感受 0 = 1 是如何在程序中被 zero && fabsl(a[i][n]) >= EPS 逮到的。

第二题：检查自由变量 输入第一行 2 3，两行分别为 0 1 1 3、0 2 2 6，输出应为 INF。

训练指引：未知数一开始就有整列为零。这个样例能检查出你的行号 row 是否错误地跟着列号 col 一起无脑增加了。

第三题：检查换行操作 输入第一行 2 2，两行分别为 0 1 2、1 1 5，输出应为 UNIQUE 及答案 3.0000000000 2.0000000000。

训练指引：如果跑出了除零错误，说明主元搜索和 swap 操作还没写好。

进阶挑战：自己造数据 先选一个整数向量当答案，再随意写几行小整数系数，右端用这组答案代入算出来。这样至少知道方程组有解；随后故意添加重复行、全零行，或把同一条方程的右边乱改掉，去检查程序的冗余处理、自由度探测和冲突捕捉。

高斯消元