蒟蒻林荫小复习——高斯消元
在无数次逃避之后，林荫终于鼓起了勇气，迈出了向数论进击的第一步！
高斯消元：众所周知就是高斯解方程用的消元方法。
至于啥是消元法：将方程组中的一方程的未知数用含有另一未知数的代数式表示，并将其代入到另一方程中，这就消去了一未知数，得到一解；或将方程组中的一方程倍乘某个常数加到另外一方程中去，也可达到消去一未知数的目的。消元法主要用于二元一次方程组的求解。
听着是不是很强？
然而这就是我们在初一学的解二元一次方程的扩展：尝试将一个未知量用其他的未知量表示，直到最后得到一元一次方程，解开后再反向带入回去即可。
前置芝士：无！
实际上这个是需要懂得一些矩阵的知识的，但是我们可以用普通的方式来解释它（NOIP范围内）
题目：给定N个N元一次多项式，求出解
先找个样例：
3
1 3 4 5
1 4 7 3
9 3 2 2
每一行前3个数分别是3个元（x,y,z）的系数，最后一个数是等式的值
下面我们可以列出一个方阵（不是矩阵）：
1 3 4 5
1 4 7 3
9 3 2 2
观察如下方阵，我们可知：
- 上下交换任意两整个横行对方阵没有影响（相当于把方程式给出的顺序变一下）
- 整个横行中所有数同时乘k（实数）对方阵没有影响（相当于把方程式左右同时乘k）
那么由于我们的目的是将前面n*n的方阵消成只有从左上到右下的对角线为1，其他位置为0，
那么这个时候第i行的第n+1列就是第i个未知数的值
现在我们来确定一下消元的步骤:
- 现在是确定第i个元,找到一个之前没用过的方程记为第i个等式（你总不能通过一个等式同时将两个元用其他元表示吧，如果这样你一定会得到一个恒等式。这也就是为啥要解n元一次方程一定要有n个不同的n元一次方程），并将其用之前提到的性质1转换到第i横行的位置（方便下面继续寻找）
- 将这个等式同时除以第i个元的系数，使得可以用其他元表示第i个元（假设2x+3y+4z=10，当前元是y，那么等式会变成2/3x+y+4/3z=10/3，实际上这种情况是不可能存在的，因为在消去y之前，x一定已经被消去，x的系数会变为0，这里举出只是为了说明怎么除）
- 然后开始用这个等式对其他等式进行操作（实际上相当于将第i+1个等式的第i+1个元用其他元表示）
- 当对每一个元进行过操作之后，矩阵的第n+1列就是每个元的值。
好了，现在小朋友们可能有一个问题，这样算来算去元消在哪里了呢？
我们现在开始考虑这个问题：
当我们对第一个元进行过操作之后，除了第一个式子以外，其他式子中的第一个元都已经由其他元的组合代替。
第二个元的操作也是一样，只有第二个式子中第二个元的系数为1，其他式子（包括第1个式子）的第二个元已经全部有3，4，5，6......等元替代
因为在操作之前第二个式子中的第一个元的系数已经为0(被第1个式子消去了）.这样的话，第i个式子在算法最后只有第i个元的系数为1，其余都是0，这样就达到了我们的目的。
1 3 4 5 0 1 3 -2 0 -24 -34 -43 1 0 -5 11 0 1 3 -2 0 0 38 -91 1 0 0 -0.973684 0 1 0 5.18421 0 0 1 -2.39474
把样例运行之后，输出对每个元操作结束后的方阵就得到了上面的东西
可以看到，算法每次将一个纵列上除了当前第i个元所在位置以外的常数变成0
这样的话，我们最后就可以求出解。
但是，既然是个多元方程，就有可能出现无解和有无穷多个解的情况！
那么怎么判断无解呢？
当算法进行到某个元i时，无法找到一个没有用过且该元系数不为0的方程
换句话说，如果给了n个n元一次方程，其中有一个系数为0，就相当于给了一个n-1元一次方程，可以直接判定无解。
先来份代码！
#include<iostream> #include<cstdio> using namespace std; int n; double wws[110][110]; int main() { scanf("%d",&n); for(int i=1;i<=n;i++) { for(int j=1;j<=n+1;j++) { cin>>wws[i][j]; } } for(int i=1;i<=n;i++)//枚举每个元 { int pl=i; while(pl<=n&&wws[pl][i]==0) { pl++; } if(pl>n) { cout<<"No Solution"<<endl; return 0; } for(int j=1;j<=n+1;j++)//把第一个合法的行翻到这个元的位置 { swap(wws[pl][j],wws[i][j]); } double k=wws[i][i]; for(int j=1;j<=n+1;j++) { wws[i][j]=wws[i][j]/k;//对自己行的处理，确定当前目标元的常数=1 } for(int m=1;m<=n;m++) { if(m==i) continue; double ki=wws[m][i]; for(int k=1;k<=n+1;k++) { wws[m][k]-=ki*wws[i][k];//因为ws[i]中存有当前元常数为1时每个元对应的常数 //因此消去时枚举每个方程，ki为该方程中当前元常数，对其他元依次消去 } } } for(int i=1;i<=n;i++) { printf("%.2lf\n",wws[i][n+1]); } return 0; }
无解判完了那怎样判断无穷解？
这个先等我研究研究，马上回来填坑！
再填个小坑
wws[m][k]-=ki*wws[i][k];
代码中这一行可能有点难度哈
假定我们现在有个方程：x+2y+3z=10
那么x用y和z表示就是x=10-2y-3z
那么我们现在尝试将上述x带入2x+3y+z=15中
可得20-4y-6z+3y+z=15
转化一下（3-4）y+（1-6）z=15-20
然后和上面代码结合一下：
for(int m=1;m<=n;m++)
{//枚举每个横行
if(m==i)//是第i行就不管(前面已经处理过）
continue;
double ki=wws[m][i];//ki是在第m行中元i的系数（相当于上面2x+3y+z=15中的2，此时元i为x）
for(int k=1;k<=n+1;k++)
{//wws[i][k]可以认为是在用其他元表示元i的表达式中元k的个数的相反数，也是式子i中元k的系数（假设当前元k为上面的y，那么wws[i][k]==2）
wws[m][k]-=ki*wws[i][k];//因为前面存的是相反数，所以要用减法。
//因此消去时枚举每个方程，ki为该方程中当前元常数，对其他元依次消去
}
}
判断误解实际上可能还有一种情况,就是式子足够n个,但是有些式子本质是一样的
那么我们消元消完后,可以发现,如果有某一行系数全部是0了,但是常数并不是0,这个时候就代表无解
至于多解,在确定有解后,找一下看有多少个式子系数和常数全是0
这样就代表这个元可以随便取值,怎么取都可以.
那自然就多解了.
UPD: 之前的消元方法不能很好的判断是否有多解或无解
下面给出新的代码
#include<iostream> #include<cstdio> using namespace std; double A[55][55],EPS=1e-8; int n,t; double ABS(double x) { return max(x,-x); } void Work() { //freopen("gaess.in","r",stdin); //freopen("gaess.out","w",stdout); scanf("%d",&n); for(int i=1;i<=n;i++) { for(int j=1;j<=n+1;j++) scanf("%lf",&A[i][j]); } int row,i; for(i=1,row=1;row<=n&&i<=n;i++,row++) { //row代表行,i代表决策元 int t=row; for(int j=row+1;j<=n;j++) { if(ABS(A[j][i])>ABS(A[t][i])) t=j; } if(t!=row) for(int j=1;j<=n+1;j++) { swap(A[t][j],A[row][j]); } if(ABS(A[row][i])<EPS) { row--;continue; } //代表这个元没能找到合适求解的式子,先放一放,式子下面还能用 for(int j=row+1;j<=n;j++) { if(ABS(A[j][i])>EPS) { double Kel=A[j][i]/A[row][i]; for(int k=1;k<=n+1;k++) { A[j][k]-=A[row][k]*Kel; } } } } row--; for(int j=row;j<=n;j++) { if(ABS(A[j][n])<EPS&&ABS(A[j][n+1]>EPS)) { printf("%d\n",-1); return; } } if(row<n) { printf("%d\n",0); return; } for(int j=row;j>=1;j--) { for(int k=j+1;k<=n;k++) { A[j][n+1]-=A[j][k]*A[k][n+1]; } A[j][n+1]/=A[j][j]; } for(int i=1;i<=n;i++) { if(A[i][n+1]) { cout<<"x"<<i<<"="; printf("%.2lf\n",A[i][n+1]); } else { cout<<"x"<<i<<"="<<0<<endl; } } return; } int main() { Work(); return 0; }
多解输出0,无解输出-1
完结撒花！
