高斯消元

---

简单的讲，高斯消元就是模拟小学生解多元一次方程组的过程。只不过这种方法更有规律可循，更适合计算机去解决。

对于方程组

$$
\begin{Bmatrix}
k_{11}x + k_{12}y + k_{13}z = d_1 \\
k_{21}x + k_{22}y + k_{23}z = d_2 \\
k_{31}x + k_{32}y + k_{33}z = d_3 \\
\end{Bmatrix}$$

我们把它保存为一个增广矩阵, 这个矩阵的意思就是每个未知数在每个方程中的系数，添上一列为未知数与系数积的和。

$$
\begin{Bmatrix}
k_{11} & k_{12} & k_{13} & | & d_1 \\
k_{21} & k_{22} & k_{23} & | & d_2 \\
k_{31} & k_{32} & k_{33} & | & d_3 \\
\end{Bmatrix}
$$

然后通过以下几种简单的操作，将这个矩阵变换

1. 将矩阵的某一行同时乘一个数

2. 将矩阵的某一行同时减去另一行

3. 交换两行的顺序（即交换方程组中两个方程的顺序）

通过这些操作，我们的目标是将矩阵变换为如下的形式 (就简称三角矩阵吧)

$$
\begin{Bmatrix}
a_{11} & a_{12} & a_{13} & | & d_1 \\
0 & a_{22} & a_{23} & | & d_2 \\
0 & 0 & a_{33} & | & d_3 \\
\end{Bmatrix}
$$

这样，可知 $z=\frac{d_3}{a_{33}}$

将 $z$代入 2 式又知 $y=\frac{d_2-a_{23}z}{a_{22}}$

将 $y$,$z$带入 1 式又知 $x=\frac{d_1-a_{13}z-a_{12}y}{a_{11}}$

这个过程也可以称为 “回代”

然而，如何获得一个三角矩阵呢？

这其实就是我们常用的加减消元法，对于 A 式和 B 式，通过给 A 式乘一个数，使两式主元系数相同，再由 B 式减去 A 式即可。这样原矩阵可以变成：(用 2,3 式减去 1 式)

$$
\begin{Bmatrix}
x & x & x & | & x \\
0 & x & x & | & x \\
0 & x & x & | & x \\
\end{Bmatrix}
$$

进而变成 (用 3 式减去 2 式)

$$
\begin{Bmatrix}
x & x & x & | & x \\
0 & x & x & | & x \\
0 & 0 & x & | & x \\
\end{Bmatrix}
$$

由此获得了一个三角矩阵。这个过程也称为 “消元”。

---

总结一下，高斯消元的总体过程就是：

消元->回代 _(废话)

但是这也会遇到几个问题：1. 由于不断的乘除操作，未知数系数的误差会逐渐增大，以至结果十分离谱。2. 如果方程无解或有无数解，我们又该如何知道呢?

对于问题 1，我们可以采用分子分母表示法计算，或用实数计算时采用 “列主消元”，就是每次用主元（要消的元）最大的方程去减其余的方程。请自行百度。

对于问题 2，我们分情况讨论

1. 可以变换成为严格的三角矩阵：有唯一解

如果无法构造出严格的三角矩阵，那么至少可以构造出如下的阶梯矩阵。

$$
\begin{Bmatrix}
x & x & x & x & … & x & | & x_1 \\
0 & x & x & x & … & x & | & x_2 \\
0 & 0 & 0 & x & … & x & | & x_3 \\
… & … & … & … & … & x & | & … \\
0 & 0 & 0 & 0 & … & x & | & x_{n-2} \\
0 & 0 & 0 & 0 & … & 0 & | & x_{n-1} \\
0 & 0 & 0 & 0 & … & 0 & | & x_{n} \\
\end{Bmatrix}
$$

2. 这个矩阵中如果系数全部为 0 的几行中（即 $x_{n-2},x_{n-1},x_{n}$所在这几行），存在 $x_{i}$不为 0，那么方程组无解。

3. 如果 $x_{i}$都为 0，那么方程组有无数组解。

在这里给出一个没有解个数判断的，不完备的，模块化的，分数表示的，不知道对不对的，贼 jb 长的：高斯消元模板。

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cstring>

#define MX 1001

using namespace std;

typedef long long ll;
typedef struct frc
{
    ll up,dn;
}frac;

int n;
frac mat[MX][MX];
frac ans[MX];

ll gcd(ll a,ll b)
{
    return (b==0?a:gcd(b,a%b));
}

frac mul(frac a,frac b)
{
    frac nxt;
    nxt.up=a.up*b.up;
    nxt.dn=a.dn*b.dn;
    ll gd=gcd(abs(nxt.up),abs(nxt.dn));
    nxt.up/=gd,nxt.dn/=gd;
    return nxt;
}

frac chu(frac up,frac down)
{
    frac nxt;
    nxt.up=up.up*down.dn;
    nxt.dn=up.dn*down.up;
    ll gd=gcd(abs(nxt.up),abs(nxt.dn));
    nxt.up/=gd,nxt.dn/=gd;
    return nxt;
}

frac mns(frac a,frac b)
{
    frac nxt;
    nxt.up=a.up*b.dn-a.dn*b.up;
    nxt.dn=a.dn*b.dn;
    ll gd=gcd(abs(nxt.dn),abs(nxt.up));
    nxt.up/=gd,nxt.dn/=gd;
    return nxt;
}

void outf(frac a)
{
    printf("%.3f ",((double)a.up)/((double)a.dn));
}

void onelization(frac (*eqt),int yuan,int len)
{
    eqt[0]=chu(eqt[0],eqt[yuan]);
    for(int i=len;i>=yuan;i--)eqt[i]=chu(eqt[i],eqt[yuan]);
}

void elimination(frac (*use),frac (*tar),int len)
{
    for(int i=0;i<=len;i++)tar[i]=mns(tar[i],use[i]);
}

void trilization(frac (*tar)[MX],int num)
{
    int mxpos;
    frac tme;
    for(int i=1;i<=num;i++)
    {
        mxpos=i;
        for(int j=i;j<=num;j++)if(tar[mxpos][i].up*tar[j][i].dn<tar[j][i].up*tar[mxpos][i].dn)mxpos=j;
        swap(tar[mxpos],tar[i]);
        onelization(tar[i],i,num);
        for(int j=i+1;j<=num;j++)
        {
            tme=tar[j][i];
            for(int k=0;k<=num;k++)
                tar[j][k]=mns(tar[j][k],mul(tar[i][k],tme));
        }
    }
}

void anslization(frac (*tar)[MX],int num)
{
    frac now;
    for(int i=num;i>=1;i--)
    {
        now=tar[i][0];
        for(int j=num;j>i;j--)now=mns(now,mul(tar[i][j],ans[j]));
        ans[i]=now;
    }
}

void input()
{
    scanf("%d",&n);
    for(int i=1;i<=n;i++)
    {
        for(int j=1;j<=n;j++)scanf("%I64d",&mat[i][j].up),mat[i][j].dn=1;
        scanf("%I64d",&mat[i][0].up),mat[i][0].dn=1;
    }
}

int main()
{
    input();
    trilization(mat,n);
    anslization(mat,n);
    for(int i=1;i<=n;i++)outf(ans[i]);
    cout<<endl;
    return 0;
}