【模板】高斯消元

我们从一个高斯消元开始：

struct row{
	double a[105];
	row operator -(const row& other)const{
		row b;
		for(int i=1;i<=m;i++) b.a[i]=a[i]-other.a[i];
		return b;
	}
	row operator *(const double& div)const{
		row b;
		for(int i=1;i<=m;i++) b.a[i]=a[i]*div;
		return b;
	}
	row operator /(const double& div)const{
		row b;
		for(int i=1;i<=m;i++) b.a[i]=a[i]/div;
		return b;
	}
	bool iszero(){
		bool ret=1;
		for(int i=1;i<m;i++) if(a[i]) ret=0;
		return ret;
	}
}r[105];
int Gauss_Elimination() 
{
	for(int i=1;i<m;i++){
		if(!r[i].a[i]){
			int f=0;
			for(int j=i+1;j<=n;j++)
				if(r[j].a[i])
					f=j;
			if(f) swap(r[i],r[f]);
		}
		if(r[i].iszero()){
			for(int j=i;j<=n;j++)
				if(r[j].a[m]!=0) return NO_SOLUTION;
			return UNLIMITED_SOLUTIONS;
		}
		int fir=0;
		for(int j=1;j<=m;j++)
			if(fabs(r[i].a[j])>1e-9) {fir=j;break;}
		for(int j=i+1;j<=n;j++) r[j]=r[j]-r[i]*(r[j].a[fir]/r[i].a[fir]);
		r[i]=r[i]/(r[i].a[fir]);
	}
	for(int i=m-1;i>=1;i--) 
		for(int j=m;j>i;j--)
			r[i]=r[i]-r[j]*(r[i].a[j]);
	return SOLVED;
}

值得注意的是，我发现很多人的高斯消元板子比我短很多，仔细观看后发现主要有两原因：
1.我是完全按照传统线代教材：先化行阶梯矩阵，最后最简化；但是实际可以优化成直接最简；
2.实际没有必要判断该行是不是头非零，如果是零去找人更换，直接从这行开始找，非零更换即可（等价）
3.实际可以按列遍历，这列没有行就看下一列，此时行还是这一行。如果找到行，才在下一列换下一行。（这个倒没太大所谓）

于是我们最后的模板为

struct row{
	double a[105];
	row operator -(const row& other)const{
		row b;
		for(int i=1;i<=m;i++) b.a[i]=a[i]-other.a[i];
		return b;
	}
	row operator *(const double& div)const{
		row b;
		for(int i=1;i<=m;i++) b.a[i]=a[i]*div;
		return b;
	}
	row operator /(const double& div)const{
		row b;
		for(int i=1;i<=m;i++) b.a[i]=a[i]/div;
		return b;
	}
	bool iszero(){
		bool ret=1;
		for(int i=1;i<m;i++) if(a[i]) ret=0;
		return ret;
	} 
}r[105];
int Gauss_Elimination() 
{
	for(int i=1;i<m;i++)
	{
		for(int j=i;j<=n;j++)
			if(r[j].a[i])
			{
				swap(r[i],r[j]);
				break;
			}
		if(r[i].iszero())
		{
			for(int j=i;j<=n;j++)
				if(r[j].a[m]!=0) return NO_SOLUTION;
			return UNLIMITED_SOLUTIONS;
		}
		int fir=0;
		for(int j=1;j<=m;j++)
			if(fabs(r[i].a[j])>1e-6)
			{
				fir=j;
				break;
			}
		for(int j=1;j<=n;j++)
		{
			if(j==i) continue;
			r[j]=r[j]-r[i]*(r[j].a[fir]/r[i].a[fir]);
		}
		r[i]=r[i]/(r[i].a[fir]);
	}
	return SOLVED;
}