[hihocdoer#1195 : 高斯消元·一] 高斯消元模板

题目链接：[hihocdoer#1195 : 高斯消元·一]

题解思路：

小Ho：<吧唧><吧唧><吧唧>

小Hi：小Ho，你还吃呢。想好了么?

小Ho：肿抢着呢(正想着呢)...<吞咽>...我记得这个问题上课有提到过，应该是一元一次方程组吧。

我们把每一件商品的价格看作是x[1]..x[n]，第i个组合中第j件商品数量记为a[i][j]，其价格记作y[i]，则可以列出方程式：

a[1][1] * x[1] + a[1][2] * x[2] + ... + a[1][n] * x[n] = y[1]a[2][1] * x[1] + a[2][2] * x[2] + ... + a[2][n] * x[n] = y[2]                                  ...a[m][1] * x[1] + a[m][2] * x[2] + ... + a[m][n] * x[n] = y[m]

我们可以对方程组进行3种操作而不改变方程组的解集：

1. 交换两行。

2. 把第i行乘以一个非0系数k。即对于j = 1..n， 令a[i][j] = k*a[i][j], y[i]=k*y[i]

3. 把第p行乘以一个非0系数k之后加在第i行上。即对于j=1..n， 令a[i][j] = a[i][j]+k*a[p][j]，y[i]=y[i]+k*p[i]

以上三个操作叫做初等行变换。

我们可以使用它们，对这个方程组中的a[i][j]进行加减乘除变换，举个例子：

a[1][1] * x[1] + a[1][2] * x[2] + ... + a[1][n] * x[n] = y[1]    式子(1)a[2][1] * x[1] + a[2][2] * x[2] + ... + a[2][n] * x[n] = y[2]    式子(2)

我们可以通过 式子(1) - 式子(2) * (a[1][1] / a[2][1])，将第1行第1列的a[1][1]变换为0。

对整个方程组进行多次变换之后，可以使得a[i][j]满足：

a[i][j] = 1 (i == j)a[i][j] = 0 (i != j)

则整个方程组变成了：

x[1] = y'[1]x[2] = y'[2]...x[n] = y'[n]0 = y'[n + 1]0 = y'[n + 2]...0 = y'[m]

这样的话，y'[1] .. y'[n]就是我们要求的x[1]..x[n]

小Hi：挺不错的嘛，继续？

小Ho：好，关于如何变换，我们可以利用一个叫高斯消元的算法。高斯消元分成了2个步骤：

首先我们要计算出上三角矩阵，也就是将方程组变为：

a[1][1] * x[1] + a[1][2] * x[2] + ... + a[1][n] * x[n] = y'[1]      0 * x[1] + a[2][2] * x[2] + ... + a[2][n] * x[n] = y'[2]      0 * x[1] +       0 * x[2] + ... + a[3][n] * x[n] = y'[3]                                   ...      0 * x[1] +       0 * x[2] + ... + a[n][n] * x[n] = y'[n]      0 * x[1] +       0 * x[2] + ... +       0 * x[n] = y'[n + 1]                               ...      0 * x[1] +       0 * x[2] + ... +       0 * x[n] = y'[m]

也就是通过变换，将所有a[i][j](i>j)变换为0。同时要保证对角线上的元素a[i][i]不为0。

方法也很见简单，从第1行开始，我们利用当前行第i列不为0，就可以通过变换将i+1..M行第一列全部变换为0，接着对于第2行，我们用同样的方法将第3..M行第2列也变换为0...不断重复直到第n行为止。

假如计算到第i行时，第i列已经为0，则我们需要在第i+1..M行中找到一行第i列不为0的行k，并交换第i行和第k行，来保证a[i][i] != 0。但这时候还有可能出现一个情况，就是第i..M行中的i列均为0，此时可以判定，该方程组有多解。

当得到上三角矩阵后，就可以从第n行开始逆推，一步一步将a[i][j](i<j)也变换为0.

因为第n行为a[n][n] * x[n] = y'[n]，则x[n] = y'[n] / a[n][n]。

第n-1行为a[n-1][n-1] * x[n - 1] + a[n][n] * x[n] = y'[n - 1]。我们将得到的x[n]代入，即可计算出x[n-1]。

同样的依次类推就可以得到所有的x[1]..x[n]。

而对于多解和无解的判定：

当在求出的上三角矩阵中出现了 a[i][1] = a[i][2] = ... = a[i][n] = 0, 但是y'[i] != 0时，产生了矛盾，即出现了无解的情况。

而多解的证明如下：

假设n=3,m=3，而我们计算出了上三角矩阵为：

a * x[1] + b * x[2] + c * x[3] = d                      e * x[3] = f                             0 = 0

当我们在第一个式子中消去x[3]后，有a * x[1] + b * x[2] = g，显然x[1]和x[2]有无穷多种可能的取值。

小Hi：既然小Ho你都已经把整个算法讲了，那么我就只能给出伪代码了：

// 处理出上三角矩阵For i = 1..N    Flag ← False    For j = i..M                // 从第i行开始，找到第i列不等于0的行j        If a[j][i] != 0 Then            Swap(j, i)          // 交换第i行和第j行            Flag ← True            Break        End If    End For    // 若无法找到，则存在多个解    If (not Flag) Then        manySolutionsFlag ← True // 出现了秩小于N的情况        continue;    End If    // 消除第i+1行到第M行的第i列    For j = i+1 .. M        a[j][] ← a[j][] - a[i][] * (a[j][i] / a[i][i])        b[j] ← b[j] - b[i] * (a[j][i] / a[i][i])    End ForEnd For // 检查是否无解，即存在 0 = x 的情况For i = 1..M    If (第i行系数均为0 and (b[i] != 0)) Then        return "No solutions"    End IfEnd For// 判定多解If (manySolutionsFlag) Thenreturn "Many solutions"End If// 此时存在唯一解// 由于每一行都比前一行少一个系数，所以在M行中只有前N行有系数// 解析来从第N行开始处理每一行的解For i = N .. 1    // 利用已经计算出的结果，将第i行中第i+1列至第N列的系数消除    For j = i + 1 .. N        b[i] ← b[i] - a[i][j] * value[j]        a[i][j] ← 0    End For    value[i] ← b[i] / a[i][i]End For

那最后能够拜托你实现一下这个算法么？

小Ho：没问题，等我吃完这包薯片就去！

#include <set>#include <stack>#include <queue>#include <cmath>#include <cstdio>#include <string>#include <cstring>#include <iostream>#include <algorithm>using namespace std;//#pragma comment(linker, "/STACK:1024000000,1024000000")#define FIN             freopen("input.txt","r",stdin)#define FOUT            freopen("output.txt","w",stdout)#define fst             first#define snd             second//typedef __int64 LL;typedef long long LL;typedef pair<int, int> PII;const double eps = 1e-6;const int MAXN = 1000 + 5;const int MAXE = 1000 + 5;int T, N, M, K;struct Gauss {    // 等式的个数和变元的个数    int equ, var;    // 系数矩阵    double a[MAXN][MAXN];    // 常数列    double b[MAXN];    // 答案    double x[MAXN];    void read(int _equ, int _var) {        equ = _equ, var = _var;        for (int i = 1; i <= equ; i++) {            for (int j = 1; j <= var; j++) {                scanf("%lf", &a[i][j]);            }            scanf("%lf", &b[i]);        }    }    /**     * 交换r1, r2两行     */    void chg_row(int r1, int r2) {        for (int i = 1; i <= var; i++) {            swap(a[r1][i], a[r2][i]);        }        swap(b[r1], b[r2]);    }    /**     * 求解多元方程组     * @return [-1——无解, 0——无穷解, 1——唯一解]     */    int run() {        bool flag, manysolutionsflag = false;        // 处理出上三角矩阵        for (int i = 1; i <= var; i++) {            flag = false;            // 从第i行开始，找到第i列不等于0的行j            for (int j = i; j <= equ; j++) {                if (fabs(a[j][i]) > eps) {                    chg_row(i, j);                    flag = true;                    break;                }            }            // 若无法找到,即出现了秩小于N的情况,则存在多个解            if (!flag) {                manysolutionsflag = true;                continue;            }            // 消除第i+1行到第M行的第i列            for (int j = i + 1; j <= equ; j++) {                double t = a[j][i] / a[i][i];                for (int k = 1; k <= var; k++) {                    a[j][k] -= a[i][k] * t;                }                b[j] -= b[i] * t;            }        }        // 检查是否无解，即存在 0 = x 的情况        for (int i = 1; i <= equ; i++) {            flag = false;            // 判断第 i 行系数是否全为0, 并且b[i] != 0 的情况            for (int j = 1; j <= var; j++) {                if (fabs(a[i][j]) > eps) {                    flag = true;                    break;                }            }            if (!flag && fabs(b[i]) > eps) {                return -1;            }        }        // 判定多解        if (manysolutionsflag) return 0;        // 此时存在唯一解        // 由于每一行都比前一行少一个系数，所以在M行中只有前N行有系数        // 解析来从第N行开始处理每一行的解        for (int i = var; i >= 1; i--) {            // 利用已经计算出的结果，将第i行中第i+1列至第N列的系数消除            for (int j = i + 1; j <= var; j++) {                b[i] -= a[i][j] * x[j];                a[i][j] = 0;            }            x[i] = b[i] / a[i][i];        }        return 1;    }} gauss;int main() {#ifndef ONLINE_JUDGE    FIN;#endif // ONLINE_JUDGE    while (~scanf ("%d %d", &N, &M) ) {        gauss.read(M, N);        int ret = gauss.run();        switch (ret) {        case -1:            printf("No solutions\n");            break;        case 0:            printf("Many solutions\n");            break;        case 1:            for (int i = 1; i <= gauss.var; i++) {                printf("%d\n", (int)(gauss.x[i] + 0.5));            }            break;        }    }    return 0;}