高斯消元

高斯消元
参考资料

回到顶部

高斯消元

这里介绍的是高斯-约旦消元法。

相对于传统的高斯消元，约旦消元法的精度更好、代码更简单，没有回带的过程。

约旦消元法大致思路如下：

1.选择一个尚未被选过的未知数作为主元，选择一个包含这个主元的方程。

2.将这个方程主元的系数化为1。

3.通过加减消元，消掉其它方程中的这个未知数。

4.重复以上步骤，直到把每一行都变成只有一项有系数。

注意，高斯消元可能会有精度问题，因此我们要取主元系数最大的式子去更新其他式子

例题一P3389 【模板】高斯消元法

在这里，无解和无穷解我们都输出No Solution

代码

(不够严谨的，无法具体判断无解还是无数解)

#include<bits/stdc++.h>
using namespace std;
const int maxn=105;
int n;
double a[maxn][maxn],ANS[maxn];
void gauss(){
	for(int i=1;i<=n;i++){
		int maxx=i;
		for(int j=i+1;j<=n;j++)
			if(fabs(a[maxx][i])<fabs(a[j][i]))maxx=j;
		if(i!=maxx)swap(a[i],a[maxx]);
		if(!a[i][i])continue;
		for(int j=1;j<=n;j++){
			if(j==i)continue;
			double tmp=a[j][i]/a[i][i];
			for(int k=i;k<=n+1;k++)
				a[j][k]-=a[i][k]*tmp;
		}
	}
	return ;
}
int main(){
	cin>>n;
	for(int i=1;i<=n;i++)
		for(int j=1;j<=n+1;j++)
			cin>>a[i][j];
	gauss();
	for(int i=1;i<=n;i++){
		if(a[i][i]==0){
			cout<<"No Solution\n";
			return 0;
		}
		else
		ANS[i]=a[i][n+1]/a[i][i];
	}
	for(int i=1;i<=n;i++)printf("%.2f\n",ANS[i]);
	return 0;
}

例题二P2455 [SDOI2006]线性方程组

这一次，我们要具体判断无解与无穷解

题意：在高斯消元的基础上，判断是正常的还是“无解”还是“有无数解”。

简述一下高斯消元的思想：让第i个式子对应第i元，并用第i个式子消去其它式子上的第i元，最终使方程组呈现出类似一条对角线的样子，即第i 个式子只在第 i元和等式右边有值，其余都是0。

一般的高斯消元每次会选择该项系数绝对值最大的来对应第 i元。这样，如果发现有一元系数为0，就能说明方程没有正常的解了。

那么我们来研究，没有正常的解，究竟是“无解”还是“有无数解”。

一个 $n a i v e$ 的想法是，在完成消元后枚举每一项，如果发现系数是0则讨论等式右边：如果为0，则有无数解；如果不为0，则无解。

其中无解优先级更高，即如果一个方程有一项无解且有一项有无数解，那么判定无解。

交上去喜提90。

存在漏洞：

2
0 2 3
0 0 0

答案显然是0，但是我们的程序会说这是-1。

我们发现，我们可以根据第一个式子来求出第二元，但是程序会用它来算第一元，并且以后算第二元的时候不会再考虑该式，

因为我们认为它已经对应第一元了！

也就是说，我们的答案受到消元顺序影响。

注意到，我们的症结在于把还有用的式子放到了用不上它的地方，并且以后也不来看它是否可用。

那么，我们就来看它是否可用。

原本的高斯消元中，我们认为只有 i之后的式子是可用的，因为我们不用管无解还是无数解，系数为0直接判掉。

但这里，我们有可能会在系数为0之后继续做下去。这就是受到消元顺序影响的原因。

那么，可用的就不仅仅是之后的式子，还有之前系数为0的式子。

于是解决了。

代码

#include<bits/stdc++.h>
using namespace std;
const int maxn=105;
const double eps=1e-7;
int n;
double a[maxn][maxn],ANS[maxn];
void gauss(){
	for(int i=1;i<=n;i++){
		int maxx=i;
		for(int j=1;j<=n;j++){//从i-n改为了1-n 
			if(fabs(a[j][j])>=eps&&j<i)continue;//增加了这句话 
			if(fabs(a[maxx][i])<fabs(a[j][i]))maxx=j;
		}
		if(i!=maxx)swap(a[i],a[maxx]);
		if(!a[i][i])continue;
		for(int j=1;j<=n;j++){
			if(j==i)continue;
			double tmp=a[j][i]/a[i][i];//这句话不能放进里面去，因为a[j][i]会修改 
			for(int k=i;k<=n+1;k++)
				a[j][k]-=a[i][k]*tmp;
		}
	}
	return ;
}
int main(){
	cin>>n;
	for(int i=1;i<=n;i++)
		for(int j=1;j<=n+1;j++)
			cin>>a[i][j];
	gauss();
	bool flag=false;
	for(int i=1;i<=n;i++)if(fabs(a[i][i])<=eps&&fabs(a[i][n+1])>eps)flag=true;
	if(flag){cout<<-1<<endl;return 0;}//先判断无解 
	flag=false;
	for(int i=1;i<=n;i++)if(fabs(a[i][i])<=eps&&fabs(a[i][n+1])<=eps)flag=true;
	if(flag){cout<<0<<endl;return 0;}//再判断无数解,并注意误差 
	for(int i=1;i<=n;i++)
		ANS[i]=a[i][n+1]/a[i][i];
	for(int i=1;i<=n;i++)printf("x%d=%.2f\n",i,ANS[i]);
	return 0;
}

例题三P2447 [SDOI2010] 外星千足虫

高斯消元解异或方程组裸题

乍一看，这好像并不好做，最暴力的方法是O(2^n)枚举

但是有了高斯消元，我们就可以不这么暴力了，而是可以借助高斯消元的思想来解了

举个例子：

现在我们有一个方程组长这样：

x 1 \oplus x 3 = 0 x 1 \oplus x 2 = 0 x 1 \oplus x 2 \oplus x 3 = 1

（当然，给出的例子是很简单的了）

但是我们可以借助这种思想：我们首先观察第一个方程：我们希望只有这一个方程中含有 $x_{1}$ （这也是高斯消元的思想），这样在其他变量已知的情况下我们就可以求出 $x_{1}$ 了

那么根据异或的性质，我们用第一个方程逐个与含有 $x_{1}$ 的方程异或，就可以消去所有的 $x_{1}$ ！

于是做出来的结果就长这样：

x 1 \oplus x 3 = 0 x 2 \oplus x 3 = 0 x 2 = 1

接下来就是第二个方程：我们要把 $x_{2}$ 消掉，于是我们的方法同上

结果长这样：

x 1 \oplus x 3 = 0 x 2 \oplus x 3 = 0 x 3 = 1

最后一个方程

x 1 = 1 x 2 = 1 x 3 = 1

这就是高斯消元解异或方程组的过程及原理

但是这里会出现一个问题：“自由元”问题

按道理，给出n个方程组应该是可以解出每个答案的，但是由于异或运算特殊的性质，有时n个方程是难以解出每个值的

为什么？

举个例子：

x 1 \oplus x 2 = 1 x 2 \oplus x 3 = 1 x 1 \oplus x 3 = 0

不难发现，这组方程本身解就不唯一（显然，(1,0,1)与(0,1,0)都是该方程组的解）

（如果去探索本质的原因：上面的异或方程在逻辑上等同于以下两个方程： $x_{1} == x_{3}$ 且 $x_{1}! = x_{2}$ ，那么考虑取值的方式，不难发现解不是唯一的）

那么如果用正常的方式来解，最后会变成这种形式：

x 1 \oplus x 2 = 1 x 2 \oplus x 3 = 1 x 2 \oplus x 3 = 1

这也能看出来，因为丢失了一个方程啊！

这时我们就称 $x_{3}$ 是一个自由元，也就是他的值既可以为 $0$ 又可以为 $1$ ，此时该异或方程组有很多组解！

同时我们可以发现，自由元的取值会直接影响不是自由元的取值（就如上面的例子）

对于自由元问题，目前仅有的办法就是 $d f s$ ，枚举每个自由元的状态然后回代求解

但是对于这道题，我们没有必要这么做

由于对于解不唯一的情况，他只要求输出一句话，所以我们只需在解不唯一的情况下（也就是解了全部方程仍不能得到唯一解时）特判即可

接下来我们就可以研究如何输出最小方程数了

我们直接像上面一样一个一个元解，不过我们现在不需要找到最大的元，因为这题不会出现精度问题

因此我们找到第一个式子，我们没有使用过，并且该元的系数不为0，记录下该位置

如果我们找不到，我们就返回解不唯一

否则我们用该式子更新所有式子

另外，如果我们使用常规的写法的时间复杂度是 $O (n^{2} m)$ 的，会超时

因此我们可以使用 $b i t s e t$ 进行优化，会给复杂度除以32，所以最终的时间复杂度是 $O (\frac{n^{2} m}{32})$

代码

#include<bits/stdc++.h>
using namespace std;
const int maxn=1e3+5;
int n,m,maxx=0; 
bitset<maxn>a[2*maxn];
string x;
char y;
bool Gauss(){
	for(int i=1;i<=n;i++){
		for(int j=1;j<=m;j++){
			if(a[j][j]&&j<i||!a[j][i])continue;
			if(i!=j)
			swap(a[i],a[j]);
			maxx=max(maxx,j);
			break;
		}
		if(!a[i][i])return false;
		for(int j=1;j<=m;j++){
			if(i==j)continue;
			if(a[j][i])
			for(int k=i;k<=n+1;k++)
				a[j][k]=(a[j][k]^a[i][k]);
		}
	}
	return true;
}
int main(){
	cin>>n>>m;
	for(int i=1;i<=m;i++){
		cin>>x>>y;
		for(int j=1;j<=n;j++)a[i][j]=x[j-1]-'0';
		a[i][n+1]=y-'0';
	} 
	if(!Gauss()){
		cout<<"Cannot Determine\n";
		return 0;
	}
	cout<<maxx<<endl;
	for(int i=1;i<=n;i++){
		if(a[i][n+1])
			printf("?y7M#\n");
		else printf("Earth\n");
	}
	return 0;
}

例题四P2962 [USACO09NOV]Lights G

对于这一道题来说，因为题目要求我们求最少的操作次数，所以我们不仅要解异或方程，对于自由元我们不能直接设为0，我们要dfs一下求最小值

代码

#include<bits/stdc++.h>
using namespace std;
const int maxn=40;
int minn=1e9;
int n,m,x,y,d[maxn];
bitset<maxn>a[maxn];
void Gauss(){
	for(int i=1;i<=n;i++){
		for(int j=1;j<=n;j++){
			if(a[j][j]&&j<i||!a[j][i])continue;
			if(i!=j)swap(a[i],a[j]);
			break;
		}
		if(!a[i][i])continue;
		for(int j=1;j<=n;j++){
			if(i==j||!a[j][i])continue;
			a[j]=(a[j]^a[i]);
		}
	}
	return ;
}
void dfs(int now,int num){
	if(num>minn)return ;
	if(now==0){
		minn=min(minn,num);
		return ;
	}
	if(a[now][now]==0){
		d[now]=0;
		dfs(now-1,num);
		d[now]=1;
		dfs(now-1,num+1);
	}
	else{
		d[now]=a[now][n+1];
		for(int i=now+1;i<=n;i++)d[now]=(d[now]^(a[now][i]&d[i]));
		dfs(now-1,num+d[now]);
	}
	return ;
}
int main(){
	cin>>n>>m;
	for(int i=1;i<=m;i++){
		cin>>x>>y;
		a[x][y]=1;
		a[y][x]=1;
	}
	for(int i=1;i<=n;i++)a[i][i]=1,a[i][n+1]=1;
	Gauss();
	dfs(n,0);
	cout<<minn;
	return 0;
}

例题五P6126 [JSOI2012]始祖鸟

我们考虑要满足第 $i$ 只始祖鸟需要满足什么条件，我们可以把每一只始祖鸟分为奇数个朋友和偶数个朋友 2 类，并且定义$ x_i=1 $表示第 i 只始祖鸟去了上游，$ x_i=0$ 表示第 $i$ 只始祖鸟去了下游，并定义 $i$ 的朋友分别为 $a_{i, 1}, a_{i, 2}, \dots, a_{i, k}$

当第 $i$ 只始祖鸟的朋友数为偶数时，显然我们必须保证$x_{a_i,1} \oplus x_{a_i,2} \oplus \cdots \oplus x_{a_i,k}=1 $且$ x_i=0 $或者$ x_{a_i,1} \oplus x_{a_i,2} \oplus \cdots \oplus x_{a_i,k}=0 $且$ x_{a_i,1} \oplus x_{a_i,2} \oplus \cdots \oplus x_{a_i,k}=1 $且$ x_i=1$

然后就是异或高斯消元的模板题了。

这题对于自由元，我们可以选择直接令它为0，这样就不会对前面的产生影响了

代码

#include<bits/stdc++.h>
using namespace std;
const int maxn=2*1e3+5;
bitset<maxn>d[maxn];
int n,m,x,ANS[maxn],cnt=0;
void Gauss(){
	for(int i=1;i<=n;i++){
		for(int j=1;j<=n;j++){
			if(d[j][j]&&j<i||!d[j][i])continue;
			if(i!=j)swap(d[i],d[j]);
			break;
		}
		if(!d[i][i])continue;
		for(int j=1;j<=n;j++){
			if(i==j||!d[j][i])continue;
			d[j]=(d[j]^d[i]);
		}
	}
	return ;
}
int main(){
	cin>>n;
	for(int i=1;i<=n;i++){
		cin>>m;
		for(int j=1;j<=m;j++){
			cin>>x;
			d[i][x]=1;
		}
		if(m&1)d[i][i]=1,d[i][n+1]=1;
		else d[i][n+1]=0;
	}
	Gauss();
	bool flag=true;
	for(int i=1;i<=n;i++)
		if(!d[i][i]&&d[i][n+1])flag=false;
	if(!flag){
		printf("Impossible\n");
		return 0;
	}
	for(int i=1;i<=n;i++){
		if(d[i][i]&&d[i][n+1])ANS[++cnt]=i;
	}
	cout<<cnt<<endl;
	for(int i=1;i<=cnt;i++)
		cout<<ANS[i]<<" ";
	return 0;
}