高斯消元

基本高斯消元法

对于行列式：

∣a11a12a13...a1na21a22a23...a2na31a32a33...a3n...............an1an2an3...ann∣\left | \begin{array}{cccc} a_{1 1} & a_{1 2} & a_{1 3} & ... & a_{1 n} \\ a_{2 1} & a_{2 2} & a_{2 3} & ... & a_{2 n} \\ a_{3 1} & a_{3 2} & a_{3 3} & ... & a_{3 n} \\ ... & ... & ... & ... & ... \\ a_{n 1} & a_{n 2} & a_{n 3} & ... & a_{n n} \end{array} \right |

如果把行列式变为：

∣a11a12a13...a1n0a22a23...a2n00a33...a3n...............000...ann∣\left | \begin{array}{cccc} a_{1 1} & a_{1 2} & a_{1 3} & ... & a_{1 n} \\ 0 & a_{2 2} & a_{2 3} & ... & a_{2 n} \\ 0 & 0 & a_{3 3} & ... & a_{3 n} \\ ... & ... & ... & ... & ...\\ 0 & 0 & 0 & ... & a_{n n} \end{array} \right |

由于求行列式要从每行每列都取一个数并乘起来，所以行列式的值就显然等于从左上到右下对角线的乘积。我们可以通过高斯消元来达到这个效果。

  1. 对于上面的行列式，根据引理6，将 aija_{i j} 加上 −aija1j×a1j-\frac{a_{i j}}{a_{1 j}} \times a_{1 j} 即可把 ai1a_{i 1} 消为0。对 n−1n - 1 个行均执行一遍，即可把第一列全部消为0。

  2. 由于第一列已经处理完了，所以观察 a22a_{2 2} 到 anna_{n n} 的行列式，将 ai2a_{i 2} 消为0可将 aija_{i j} 加上 −aija2j×a2j-\frac{a_{i j}}{a_{2 j}} \times a_{2 j} ，对 n−2n - 2 个行均执行一遍，即可把第二列全部消为0。

  3. 如此将行列式的左下部分全变为0。

  • 核心思路：用第 ii 行把第 i+1i + 1 行到第 nn 行的第 ii 列消为0。

主元高斯消元法

根据语言特性，除以 10510^5 的精度比除以 10−510^{-5} 的精度要大些，而本质就在于 10510^5 绝对值更大，所以我们除的时候可以找一个当前列绝对值最大的一行与当前行交换，来提高精度。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
double z[110][110];
int n;
double gauss() {
    double x = 1.0;
    for(int i = 1; i <= n; i++) {
        for(int j = i + 1; j <= n; j++)
            if(fabs(z[j][i]) > fabs(z[i][i])) {
                x = -x;
                for(int l = 1; l <= n; l++)
                    swap(z[j][l], z[i][l]);
            }
        if(fabs(z[i][i] <= 1e-8)) return 0.0;
        for(int j = i + 1; j <= n; j++) {
            double k = z[j][i] / z[i][i];
            for(int l = 1; l <= n; l++) {
                z[j][l] = z[j][l] - z[i][l] * k;
            }
        }
    }
    for(int i = 1; i <= n; i++) x *= z[i][i];
    return x;
}

辗转相消法

对于所有数都为整数的行列式，虽然主元高斯消元法降低了误差， 但是并没有消除精度误差，计算过程中仍然在使用实数，所以可以使用类似辗转相除法的方法解决问题：用第 i+1i + 1 行减去第 ii 行，然后交换两行位置，直到第 ii 行等于0。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
int n, z[110][110];
int gauss() {
    int ans = 1;
    for(int i = 1; i <= n; i++) {
        for(int j = i + 1; j <= n; j++) {
            while(z[j][i]) {
                int k = z[i][i] / z[j][i];
                for(int l = 1; l <= n; l++) {
                    z[i][l] = z[i][l] - z[j][l] * k; //比取模快
                }
                for(int l = 1; l <= n; l++) swap(z[i][l], z[j][l]);
                ans = -ans;
            }
        }
    }
    for(int i = 1; i <= n; i++) ans *= z[i][i];
    return ans;
}