高斯消元

给定一组 nn 个未知数的方程组，有 nn 个方程

{a1,1x1+...+a1,nxn=b1...a1,1x1+...+an,nxn=bn\left\{\begin{matrix} a_{1,1}x_{1} + ... +a_{1,n}x_{n} = b_{1} \\ ... \\a_{1,1}x_{1} + ... +a_{n,n}x_{n} = b_{n} \end{matrix}\right.

转化到下面的表格

x1x_{1} x2x_{2} ...... xnx_{n} b1−nb_{1-n}
？ ？ ？ ？ ？
？ ？ ？ ？ ？
？ ？ ？ ？ ？
？ ？ ？ ？ ？

<注：下文中 ?≠0? \ne 0 >

唯一解：

x1x_{1} x2x_{2} ...... xnx_{n} b1−nb_{1-n}
？ 0 0 0 ？
0 ？ 0 0 ？
0 0 ？ 0 ？
0 0 0 ？ ？

无穷多组解：

x1x_{1} x2x_{2} ...... xnx_{n} b1−nb_{1-n}
？ 0 0 0 ？
0 ？ 0 0 ？
0 0 0 0 0
0 0 0 0 0

无解：

x1x_{1} x2x_{2} ...... xnx_{n} b1−nb_{1-n}
？ 0 0 0 ？
0 ？ 0 0 ？
0 0 0 0 ？
0 0 0 0 ？

枚举每一列c：

  1. 找到绝对值最大的一行
  2. 将该行变换到最上面
  3. 将该行第一个数变成1
  4. 将该行下面所有行的第c列变成0

数组从零储存， ai,0∼n−1a_{i,0\sim n - 1} 是系数，ai,na_{i,n} 是结果

第一部分，消元：
先看整体代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
 int c, r;
for (c = 0, r = 0; c < n; c ++ )
{
    int t = r;
    for (int i = r; i < n; i ++ )  // 找绝对值最大的行
        if (fabs(a[i][c]) > fabs(a[t][c]))
            t = i;

    if (fabs(a[t][c]) < eps) continue;

    for (int i = c; i <= n; i ++ ) swap(a[t][i], a[r][i]);  // 将绝对值最大的行换到最顶端
    for (int i = n; i >= c; i -- ) a[r][i] /= a[r][c];  // 将当前行的首位变成1
    for (int i = r + 1; i < n; i ++ )  // 用当前行将下面所有的列消成0
        if (fabs(a[i][c]) > eps)
            for (int j = n; j >= c; j -- )
                a[i][j] -= a[r][j] * a[i][c];

    r ++ ;
}

分析：

1
int c, r;

cc :当前列数
rr :当前梯数（即上文中 1∼x1 \sim x 列中的 ?? 数量）

1
2
3
for (int i = r; i < n; i ++ )
       if (fabs(a[i][c]) > fabs(a[t][c]))
           t = i;

找第 cc 列绝对值最大的行
<注：fabs()用于浮点数绝对值>

1
if (fabs(a[t][c]) < eps) continue;

如果第 cc 列绝对值最大的一行的第 cc 列的绝对值都为 00，则这一列无有效的阶梯，直接跳出

1
for (int i = c; i <= n; i ++ ) swap(a[t][i], a[r][i]);

将第 cc 列的绝对值最大第一行换到第 cc 行

1
for (int i = n; i >= c; i -- ) a[r][i] /= a[r][c];

将当前行的的第 cc 列的系数消成 11 ，方便下一步
<注：从后往前，避免 cc 列提前被改变>

1
2
3
4
for (int i = r + 1; i < n; i ++ )
       if (fabs(a[i][c]) > eps) // 误差
           for (int j = n; j >= c; j -- )
               a[i][j] -= a[r][j] * a[i][c];

用第 cc 行将下面所有行的第 cc 列的系数消成 00
< eps=10−6eps = 10^{-6} >
<注：因为第 cc 行的 cc 列的系数已经化为 11 ，所以可以将被消式减去第 cc 行的式子乘上被消行的第 cc 列的系数>
<注：从后往前，避免 cc 列提前被改变>

1
r ++;

有效阶梯数加一

第二部分：找解：
先看整体代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
if (r < n)
{
    for (int i = r; i < n; i ++ )
        if (fabs(a[i][n]) > eps)
            return 2; // 无解
    return 1; // 有无穷多组解
}

for (int i = n - 1; i >= 0; i -- )
    for (int j = i + 1; j < n; j ++ )
        a[i][n] -= a[i][j] * a[j][n];

return 0; // 有唯一解

分析：

1
if (r < n)

如果有效阶梯数小于 nn，说明无唯一解

1
2
3
4
 for (int i = r; i < n; i ++ )
    if (fabs(a[i][n]) > eps)
            return 2;
return 1;

一排系数为 00, 结果不为 00 ，显然无解，返回 22，否则显然无穷多组解，返回 11

1
2
3
4
5
for (int i = n - 1; i >= 0; i -- )
    for (int j = i + 1; j < n; j ++ )
        a[i][n] -= a[i][j] * a[j][n];

return 0; // 有唯一解

逐句分析：

  • for (int i = n - 1; i >= 0; i -- )从 xnx_{n} 开始，一个一个带入消元
  • for (int j = i + 1; j < n; j ++ )枚举，往回带入当前行的 xx
  • a[i][n] -= a[i][j] * a[j][n];最后直接输出结果，系数可以不动，改变结果，消剩下的结果便是 xx
    <注:第 nn 行是=时，由于数组默认初始化为 00 ，所以结果不变>
    图例：
    图例

AC CODEAC \text{ } CODE

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
#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>

using namespace std;

const int N = 110;
const double eps = 1e-8;

int n;
double a[N][N];

int gauss()  // 高斯消元，答案存于a[i][n]中，0 <= i < n
{
    int c, r;
    for (c = 0, r = 0; c < n; c ++ )
    {
        int t = r;
        for (int i = r; i < n; i ++ )  // 找绝对值最大的行
            if (fabs(a[i][c]) > fabs(a[t][c]))
                t = i;

        if (fabs(a[t][c]) < eps) continue;

        for (int i = c; i <= n; i ++ ) swap(a[t][i], a[r][i]);  // 将绝对值最大的行换到最顶端
        for (int i = n; i >= c; i -- ) a[r][i] /= a[r][c];  // 将当前行的首位变成1
        for (int i = r + 1; i < n; i ++ )  // 用当前行将下面所有的列消成0
            if (fabs(a[i][c]) > eps)
                for (int j = n; j >= c; j -- )
                    a[i][j] -= a[r][j] * a[i][c];

        r ++ ;
    }

    if (r < n)
    {
        for (int i = r; i < n; i ++ )
            if (fabs(a[i][n]) > eps)
                return 2; // 无解
        return 1; // 有无穷多组解
    }

    for (int i = n - 1; i >= 0; i -- )
        for (int j = i + 1; j < n; j ++ )
            a[i][n] -= a[i][j] * a[j][n];

    return 0; // 有唯一解
}


int main()
{
    scanf("%d", &n);
    for (int i = 0; i < n; i ++ )
        for (int j = 0; j < n + 1; j ++ )
            scanf("%lf", &a[i][j]);

    int t = gauss();
    if (t == 2) puts("No solution");
    else if (t == 1) puts("Infinite group solutions");
    else
    {
        for (int i = 0; i < n; i ++ )
        {
            if (fabs(a[i][n]) < eps) a[i][n] = 0;  // 去掉输出 -0.00 的情况
            printf("%.2lf\n", a[i][n]);
        }
    }

    return 0;
}