高斯消元

目录

• 0. 前言
• 1. 浮点方程组
• 2. 异或线性方程组
• 3. 整数线性方程组
• 4. 模线性方程组
• 5. 矩阵行列式

0. 前言

写这篇博客的缘由是感觉高斯消元有很多细节需要注意。所以不会讲原理，权当作一篇可供复习的资料吧。

1. 浮点方程组

例题。

1. 需要注意输出时对于 (0) 的特判。因为当 (x) 为一个接近 (0) 却为负值的时候，直接输出会变成 -0.00。
2. 但有些数据会专门卡你，所以在找主元时进行随机化也是不错的处理方案。

#include <cstdio>

#define print(x,y) write(x),putchar(y) 

template <class T> inline T read(const T sample) {
    T x=0; int f=1; char s;
    while((s=getchar())>'9'||s<'0') if(s=='-') f=-1;
    while(s>='0'&&s<='9') x=(x<<1)+(x<<3)+(s^48),s=getchar();
    return x*f;
}
template <class T> inline void write(const T x) {
    if(x<0) return (void) (putchar('-'),write(-x));
    if(x>9) write(x/10);
    putchar(x%10^48);
}

#include <cmath>
#include <iostream>
using namespace std;

const int maxn=55;
const double eps=1e-9;

double a[maxn][maxn];

void Gauss(int equ,int var) {
	int pos,row,col;
	bool che=0;
	for(row=col=1;row<=equ and col<var;++row,++col) {
		pos=row;
		for(int i=row+1;i<=equ;++i)
			if(fabs(a[i][col])-fabs(a[pos][col])>eps)
				pos=i;
		if(pos^row)
			for(int i=col;i<=var;++i)
				swap(a[pos][i],a[row][i]);
		if(fabs(a[row][col])<eps) {
			che=1,--row;
			continue;
		}
		for(int i=row+1;i<=equ;++i)
			if(fabs(a[i][col])>eps) {
				double tmp=a[i][col]/a[row][col];
				for(int j=col;j<=var;++j)
					a[i][j]-=tmp*a[row][j];
				a[i][col]=0;
				// 直接赋值为 0，减小精度误差的影响 
			}
	}
	for(int i=row;i<=equ;++i)
		if(a[i][var]>eps) return (void)(puts("-1"));
	if(che) return (void)(puts("0"));
	/*
		equ-row 就是此时的自由变元个数 
	*/
	for(int i=equ;i>=1;--i) {
		for(int j=i+1;j<=equ;++j)
			a[i][var]-=a[i][j]*a[j][var];
		a[i][var]/=a[i][i];
	}
	for(int i=1;i<=equ;++i) {
		if(fabs(a[i][var])<eps)
			a[i][var]=0;
		printf("x%d=%.2f
",i,a[i][var]);
	}
}

int main() {
	int n=read(9);
	for(int i=1;i<=n;++i)
		for(int j=1;j<=n+1;++j)
			a[i][j]=read(9);
	Gauss(n,n+1);
	return 0;
}

2. 异或线性方程组

例题。

1. 可以用 ( m bitset) 优化。
2. 题目可能询问解为 (1) 的最小个数。你可以从后往前枚举自由变元的值，从而依次解出 (x)。但是在不能确定自由变元个数的情况时，最好直接使用 ( ext{meet-in-middle})。
3. 如何理解自由变元？若没有自由变元，高斯消元回代时只有一个未知数。当出现自由变元，你会发现在这一个方程中的 (2) 个未知数都可以被视为自由变元，我们可以任取一个来赋任意值使得未知数个数又变成 (1)。

const int maxn=1005;

int n,m,x[maxn];
bitset <maxn> a[maxn<<1];

void Gauss() {
	int tmp,row=1,col=1,ans=0;
	for(;row<=m && col<n;++row,++col) {
		tmp=row;
		while(tmp<m && !a[tmp][col]) ++tmp;
		if(!a[tmp][col]) return (void)puts("Cannot Determine");
		ans=Max(ans,tmp);
		if(tmp^row) swap(a[tmp],a[row]);
		rep(i,row+1,m) {
			if(!a[i][col]) continue;
			a[i]^=a[row];
		} 
	}
	fep(i,n-1,1) {
		x[i]=a[i][n];
		rep(j,i+1,n-1) {
			if(a[i][j]) x[i]^=x[j];
		}
	}
	print(ans,'
');
	rep(i,1,n-1) puts((x[i]&1)?"?y7M#":"Earth");
}

int main() {
	n=read(9)+1,m=read(9);
	rep(i,1,m) {
		int x;
		rep(j,1,n-1) scanf("%1d",&x),a[i][j]=x;
		a[i][n]=read(9);
	}
	Gauss();
	return 0;
}

3. 整数线性方程组

1. 防止精度误差，用 ( m lcm) 来消元。

for(int i=row+1;i<equ;i++)
	if(a[i][col]!=0){
		int lcm=LCM(fabs(a[i][col]),fabs(a[row][col]));
		int ta=lcm/abs(a[i][col]),tb=lcm/abs(a[row][col]);
		if(a[i][col]*a[row][col]<0) tb=-tb; // 注意正负
		for(int j=col;j<=var;j++)
			a[i][j]=a[i][j]*ta-a[row][j]*tb;
	}

4. 模线性方程组

1. 对于模数为质数，在除以 (a[i][i]) 时可以用费马小定理求解逆元。
2. 对于模数非质数，题目就要保证 (a[i][i]) 与 ( m mod) 互质，然后用 ( m exgcd) 求解。求 (a) 在 (p) 意义下的逆元可以这样：用 ( m exgcd) 解出 (ax+bp=1)，移个项就能发现 (x) 即为所求。

void gauss() {
	int tmp;
	for(int i = 0; i <= n; ++ i) {
		/*
			用 lcm 来消元，所以主元只需要不为 0 即可
		*/
		if(! a[i][i]) {
			tmp = 0;
			for(int j = i + 1; j <= n; ++ j) if(a[j][i]) {tmp = j; break;}
			if(! tmp) continue;
			for(int j = i; j <= n + 1; ++ j) swap(a[tmp][j], a[i][j]);
		}
		for(int j = i + 1; j <= n; ++ j) {
			tmp = a[j][i];
			if(! a[j][i]) continue;
			for(int k = i; k <= n + 1; ++ k) a[j][k] = (0ll + 1ll * a[j][k] * a[i][i] % mod - 1ll * tmp * a[i][k] % mod + mod) % mod;
		}
	}
	for(int i = n; ~i; -- i) {
		for(int j = i + 1; j <= n; ++ j) a[i][n + 1] = (0ll + a[i][n + 1] - 1ll * ans[j] * a[i][j] % mod + mod) % mod;
		ans[i] = 1ll * a[i][n + 1] * qkpow(a[i][i], mod - 2) % mod;
	}
}

5. 矩阵行列式

戳这。

int Gauss() {
	int j,tmp,r=1; bool f=0;
	for(int i=1;i<=n;++i) {
		for(j=i;j<=n;++j)
			if(a[j][i]) break;
		if(j>n) return 0;
		if(i^j) swap(a[i],a[j]),f^=1;
		r=1ll*r*a[i][i]%mod;
		tmp=inv(a[i][i]);
		for(j=i;j<=n;++j)
			a[i][j]=1ll*a[i][j]*tmp%mod;
		for(j=i+1;j<=n;++j) {
			a[j][i]=mod-a[j][i];
			for(int k=i+1;k<=n;++k)
				a[j][k]=inc(a[j][k],1ll*a[j][i]*a[i][k]%mod);
			a[j][i]=0;
		}
	}
	return f?mod-r:r;
}

如果模数不为质数呢？可以模拟辗转相除法：将 (a_{i,i}) 和 (a_{j,i}) 对消，直到其中之一为 (0) 即可。所以实际上这是 (mathcal O(n^3log v))，可能稍微会慢一点。

最后 if(!a[i][i]) return 0; 其实等价于判断是否第 (i) 列全为 (0)。

int Det() {
	int pos,ret=1,tmp; bool f=0;
	for(int i=1;i<=n;++i) {
		for(int j=i+1;j<=n;++j)
			while(a[j][i]) {
				int t=a[i][i]/a[j][i];
				for(int k=i;k<=n;++k) {
					a[i][k]=(a[i][k]-1ll*a[j][k]*t%mod+mod)%mod;
					swap(a[i][k],a[j][k]);
				}
				f^=1;
			}
		if(!a[i][i]) return 0;
		ret=1ll*ret*a[i][i]%mod;
	}
	return f?mod-ret:ret;
}