0. 前言
写这篇博客的缘由是感觉高斯消元有很多细节需要注意。所以不会讲原理，权当作一篇可供复习的资料吧。
1. 浮点方程组
例题。
- 需要注意输出时对于 \(0\) 的特判。因为当 \(x\) 为一个接近 \(0\) 却为负值的时候，直接输出会变成
-0.00。 - 但有些数据会专门卡你，所以在找主元时进行随机化也是不错的处理方案。
#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\n",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. 异或线性方程组
例题。
- 可以用 \(\rm bitset\) 优化。
- 题目可能询问解为 \(1\) 的最小个数。你可以从后往前枚举自由变元的值，从而依次解出 \(x\)。但是在不能确定自由变元个数的情况时，最好直接使用 \(\text{meet-in-middle}\)。
- 如何理解自由变元？若没有自由变元，高斯消元回代时只有一个未知数。当出现自由变元，你会发现在这一个方程中的 \(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,'\n');
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. 整数线性方程组
- 防止精度误差，用 \(\rm 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. 模线性方程组
- 对于模数为质数，在除以 \(a[i][i]\) 时可以用费马小定理求解逆元。
- 对于模数非质数，题目就要保证 \(a[i][i]\) 与 \(\rm mod\) 互质，然后用 \(\rm exgcd\) 求解。求 \(a\) 在 \(p\) 意义下的逆元可以这样：用 \(\rm exgcd\) 解出 \(ax+bp=1\)，移个项就能发现 \(x\) 即为所求。
void gauss(int n) {
int j,tmp;
for(int i=1;i<=n;++i) {
for(j=i; !a[j][i] && j<=n; ++j);
if(!a[j][i]) break;
if(i^j) swap(a[i],a[j]);
int Inv = inv(a[i][i]);
for(j=i+1;j<=n;++j) if(a[j][i]) {
tmp = 1ll*a[j][i]*Inv%mod;
for(int k=i;k<=n+1;++k)
a[j][k] = (a[j][k]-a[i][k]*tmp)%mod;
}
}
for(int i=n;i;--i) {
for(j=i+1;j<=n;++j)
a[i][n+1] -= a[i][j]*a[j][j]%mod;
a[i][i] = a[i][n+1]*inv(a[i][i])%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^3\log 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;
}
