p2455才是高斯消元的真模板题好吧

高斯消元

题面：传送门

对不起我真的是太懒了

先来明确高斯消元是干什么的：

高斯消元是一种求解线性方程组的方法。

所谓线性方程组，就是由$M$个$N$元一次方程组共同构成的。

线性方程组的所有系数可以写成一个$M$行$N$列的系数矩阵，再加上每个方程等号右侧的常数，可以写成一个$M$行$N+1$列的“增广矩阵”。

$\begin{cases}x_1+2x_2-x_3=-6\2x_1+x_2-3x_3=-9\{-x_1-x_2+2x_3=7}\end{cases}\Rightarrow\left[\begin{array}{ccc|c}1 &2&-1&-6\2&1&-3&-9\{-1}&-1&2&7 \end{array}\right]
$

相信各位小学或者初中一定学过二元一次方程组吧

如果没学过那我觉得宁可以退役了

有个消元的方法相信你记忆犹新：加减消元法；

其实高斯消元法干的事情和加减消元法是差不多的；

步骤大概就是：

我们将一行加上另一行的若干倍，这样就可以消去其中一行的某一个系数；

我们可以来看一下下面的过程：

$\begin{bmatrix} 1&2&-1&-6\2&1&-3&-9\{-1}&-1&2&7\end{bmatrix}$

我们将第二行每一项分别加上第一行的$-2$倍，于是我们发现，第二行的第一项系数被消成$0$了。如下

$\begin{bmatrix} 1&2&-1&-6\0&-3&-1&3\{-1}&-1&2&7\end{bmatrix}$

接下来我们都这么操作。

对于第$n$行，我们让该行的第$n$个元素作为“主元”，并且将其他行的同列的系数都消掉，最后我们可以得到下面的形式：

$\left[
\begin{array}{ccc|c}
1&0&0&1\0&1&0&-2\0&0&1&3
\end{array}
\right]$

是的，成了这样以后，发现主对角线上的数和常数一般都不是$0$，其他位置都是$0$了（后面还有一些特殊情况与判定技巧），我们就可以愉快的写出答案：

$x_1=1$，$x_2=-2$，$x_3=3$；

是的，这就是高斯消元算法。

但是在实际应用中，我们也许会碰到很多奇奇怪怪的情况。

就例如你只掌握上面说到的东西是远远不够通过这道题的。

题目中还有一些更高端的东西：什么时候有无穷多的实数解？什么时候无解？这都是我们需要考虑的；

首先，在高斯消元的过程中，可能出现$0=d$这样的方程，其中$d$是一个非零常数，这表明这些方程中出现了矛盾，方程组无解。

其次，我们也可能出现这样的情况：

我们高斯消元结束以后，发现某一行系数全都是$0$，并且常数也是$0$.

啊就是出现了$0==0$的情况？？

是的，此时方程组有无数解，如下图：

$\left[
\begin{array}{ccc|c}
1&2&0&4\0&0&1&1\0&0&0&0
\end{array}
\right]$

那么我们总结一下：

1.当高斯消元结束以后，如果存在系数全是$0$，但是常数不是$0$的行，那么方程组无解

2.在高斯消元结束以后，如果存在系数全是$0$，并且常数也是$0$的行，那么方程组有无数解

学会了这些就能通过这些题了吗？

当然不是！

还有一个东西你没有注意到：精度损失

精度损失是怎么来的?如何降低精度损失呢?

我们在加减消元的时候，有的时候需要加的也许不是一个整数，可能是某一行的几分之几倍。

于是，考虑到这，我们发现精度损失的原因了。

所以，我们就要让产生几分之几倍（确切的说：当式子一的系数是$a_1$，式子二的系数是$a_2$时，我们加的这个数就是$a_2/a_1$倍（注意我们主元系数是在分母上））的这个过程有较小的精度损失

我们可以思考一下，我们在消元的时候一定要刻意按照题目中给的顺序一行一行来嘛？

也许可以不这样。

在什么情况下，我们能产生较小的精度损失呢？

我们举个例子来看：

考虑这样的一个方程组：

$1000x_1+ax_2+bx_3=c$;（这几个都是x不是乘号）

$0.5x_1+dx_2+ex_3=f$

我们设精度误差$eps=1e-8$；

如果我们选择$1000$作为主元，

那么理论值$ima_1=0.5/1000$

实际值$rea_1=0.5+eps/1000+eps$

误差值Δ在$rea_1-ima_1=$一个超级小的数

如果我们选择$0.5$作为主元，那么理论值$ima_2=1000/0.5$

实际值$rea_2=1000+eps/0.5+eps$；

误差值Δ在$rea_2-ima_2=40$左右

发现差距了吗？？

知道该怎么做了吗？？

是的，我们每次选择较大的作为主元就好了！

代码实现嘛，长这样的;

double maxx;//最大值
int maxline;//最大值所在行
for(re int i=1;i<=n;++i){//枚举每一列
	maxx=0,maxline=i;
	for(re int j=1;j<=n;++j)//枚举每一行的同一列位置
		if(!vis[j]&&abs(a[j][i])>maxx){//如果没来过并且很大（因为我们加的数正负都可以，所以我们取绝对值
			maxx=abs(a[j][i]),maxline=j;//记录最大值以及所在的行数
		}
           
	//........
       
	if(maxline!=i)swap(a[i],a[maxline]);//如果发现我们在消元过程中，这个位置不是需要的最大的，我们就把这一行和最大的那一行换一下位置
	vis[i]=1;//来过了
}

掌握了这些，这个题就可以$A$了

CODE：

#include<bits/stdc++.h>
using namespace std;
 
#define re register
 
const int N=105;
const double eps=1e-8;//精度误差
 
double a[N][N];//矩阵
bool vis[N];
int n;
 
inline int read(){
    int x=0,f=1;char ch=getchar();
    while(!isdigit(ch)){if(ch=='-')f=-1;ch=getchar();}
    while(isdigit(ch)){x=x*10+ch-48;ch=getchar();}
    return x*f;
}
 
int main(){
    n=read();
    
    for(re int i=1;i<=n;++i){
        for(re int j=1;j<=n+1;++j){//读入矩阵，把常数也读进去所以是n+1
 
            a[i][j]=read();
        }
    }
    
    double maxx;
    int maxline;
    for(re int i=1;i<=n;++i){
        maxx=0,maxline=i;
        for(re int j=1;j<=n;++j)
            if(!vis[j]&&abs(a[j][i])>maxx)maxx=abs(a[j][i]),maxline=j;//找最大的
        
        if(maxx<=eps)continue;//如果太小了，我们可以认为是0，并且把他跳过去，因为已经是0了不需要再操作
        if(maxline!=i)swap(a[i],a[maxline]);//将系数最大的换到i行
        vis[i]=1;
        
        for(re int j=1;j<=n;++j)//对其他行关于i列进行消元
            if(j!=i){
                double t=a[j][i]/a[i][i];
                for(re int k=i;k<=n+1;++k){
                    a[j][k]-=t*a[i][k];
                }
            }
    } 
    
    bool flag=0;//因为在消元完成后，只有主对角线上系数不是零，所以我们只需要检查一下，如果主对角线上出现了0的那一行的常数是不是0就可以了
    for(re int i=1;i<=n;++i)
        if(abs(a[i][i])<=eps){
            flag=1;//如果发现主对角线上有0存在，那么一定出了特殊情况，要么是无解，要么是无穷解
            abs(a[i][n+1])>eps?puts("-1"):puts("0");//如果常数不是零，那么可以认定就是无解，如果主对角线上有0，并且常数是0，那么就一定有无穷多解
                return 0;
        }
        
    
    for(re int i=1;i<=n;i++){//不然我们就要输出每一项的解是多少
        printf("x%d=%.2lf\n",i,a[i][n+1]/a[i][i]);//这里解释一下，为什么要用常数除以系数：因为我们在消元以后，系数并不一定变成了1，也可能是其他的常数，那么就成了kx=b的形式，所以我们拿b/k就是x了
    }
    
    return 0;
}

附上经典高斯：

inline void Gauss(){
	int maxline=0;
	double maxx=0;
	for(re int i=1;i<=n;i++){
		maxline=i,maxx=0;
		for(re int j=i+1;j<=n;j++){
			if(fabs(a[j][i])>maxx)maxx=fabs(a[j][i]),maxline=j;
		}
		if(fabs(a[maxline][i])<eps){
			puts("No Solution");
			exit(0);
		}
		if(maxline!=i)swap(a[i],a[maxline]);
		for(re int j=n+1;j>=i;j--){
			for(re int k=i+1;k<=n;k++){
				a[k][j]-=a[k][i]/a[i][i]*a[i][j];
			}
		}
	}

	for(re int i=n;i;i--){
		for(re int j=i+1;j<=n;j++){
			a[i][n+1]-=a[j][n+1]*a[i][j];
		}
		a[i][n+1]/=a[i][i];
	}
}