高斯消元

更新日志

2025/02/20：开工。

---

思路

考虑一组线性方程组：

\[\begin{cases} a_{1,1}x_1+a_{1,2}x_2+\dots+a_{1,n}x_n=b_1\\ a_{2,1}x_1+a_{2,2}x_2+\dots+a_{2,n}x_n=b_2\\ \vdots\\ a_{n,1}x_1+a_{n,2}x_2+\dots+a_{n,n}x_n=b_n \end{cases} \]

不难想到使用消元法求得所有 \(x\) 的值。

实现

高斯消元

我们考虑使用增广矩阵表示上面的方程组。

\[\begin{bmatrix} a_{1,1}&a_{1,2}&\dots&a_{1,n}&b_1\\ a_{2,1}&a_{2,2}&\dots&a_{2,n}&b_2\\ \vdots\\ a_{n,1}&a_{n,2}&\dots&a_{n,n}&b_n\\ \end{bmatrix} \]

我们考虑将整个矩阵消成下面这样倒三角的形式（\(1\) 表示此位不一定为 \(0\)）（不包括最右边的那一列）：

\[\begin{bmatrix} 1&1&1\\ 0&1&1\\ 0&0&1 \end{bmatrix} \]

对于每一个 \(x_i\)，我们都选择一个 \(a_{j,i}\ne0\) 的行 \(j\) 令 \(a_{j,i}\) 作为主元，然后对于它下面每一行，都通过整行加减把 \(a_{k,i}\) 消成 \(0\)。

然后我们就可以倒着依次得到 \(x\) 的值，具体地，得到当前 \(x\) 后通过代入法得到上一个 \(x\) 的值。

无解和无数解很好判断，如果代入完后当前 \(x\) 系数为 \(0\)：

• \(b_i=0\)，那么无数解
• 否则，无解。

码量较大，不推荐使用，故不附模板。

高斯-约旦消元

考虑去除代入操作。我们把矩阵消成下面这样（同高斯消元部分）：

\[\begin{bmatrix} 1&0&0\\ 0&1&0\\ 0&0&1 \end{bmatrix} \]

答案 \(x_i=\frac{b_i}{a_{i,i}}\)。

我们只需在选完主元以后，把所有非主元所在行这一列都消成 \(0\) 即可。

考虑如何判断无解和无数解。

为了方便，我们把所有全为 \(0\) 的行放在最后。具体的，如果当前列全是 \(0\)，就跳过当前列，下一列仍然以这一行作为主元所在行。否则，下一列就以下一行作为主元所在行。这样，前面的几行都是非全 \(0\) 的。

我们只需要判断最后几行即可，原理和高斯消元是一样的。

细节

首先是除法要判断除数不为 \(0\)。

然后是如果每次选绝对值最大的数作为主元可以减小误差。

模板

LG2455

下面是主体部分（mtx 储存增广矩阵）：

    int r=1;
    rep(i,1,n){
        int chs=0;
        rep(j,r,n)if(!chs||fabs(mtx[j][i])>fabs(mtx[chs][i]))chs=j;
        if(fabs(mtx[chs][i])<eps)continue;
        swap(mtx[chs],mtx[r]);
        rep(j,1,n)if(j!=r&&fabs(mtx[r][i])>eps)mtx[j]=mtx[j]-mtx[r]*(mtx[j][i]/mtx[r][i]);
        r++;
    }
    if(r==n+1)rep(i,1,n)cout<<'x'<<i<<'='<<mtx[i][n+1]/mtx[i][i]<<"\n";
    else{
        rep(i,r,n)if(fabs(mtx[i][n+1])>eps)cout<<-1,exit(0);//无解
        cout<<0;//无数解
    }

以及一些定义（这里对 poly 运算进行了简化，因为长度必然相等）：

const ld eps=1e-12;
#define poly vec<ld>
poly operator*(poly a,ld b){
    repl(i,0,a.size())a[i]*=b;
    return a;
}
poly operator+(poly a,poly b){
    repl(i,0,a.size())a[i]+=b[i];
    return a;
}
poly operator-(poly a,poly b){
    repl(i,0,a.size())a[i]-=b[i];
    return a;
}

以及完整代码：

#include<bits/stdc++.h>
using namespace std;

typedef long long ll;
typedef unsigned long long ull;
typedef __int128 i128;
typedef double db;
typedef long double ld;
typedef pair<int,int> pii;
typedef pair<ll,ll> pll;
typedef pair<int,ll> pil;
typedef pair<ll,int> pli;
template <typename Type>
using vec=vector<Type>;
template <typename Type>
using grheap=priority_queue<Type>;
template <typename Type>
using lrheap=priority_queue<Type,vector<Type>,greater<Type> >;
#define fir first
#define sec second
#define pub push_back
#define pob pop_back
#define puf push_front
#define pof pop_front
#define chmax(a,b) a=max(a,b)
#define chmin(a,b) a=min(a,b)
#define rep(i,x,y) for(int i=(x);i<=(y);i++)
#define per(i,x,y) for(int i=(x);i>=(y);i--)
#define repl(i,x,y) for(int i=(x);i<(y);i++)
#define file(f) freopen(#f".in","r",stdin);freopen(#f".out","w",stdout);

const int inf=0x3f3f3f3f;
const ll INF=0x3f3f3f3f3f3f3f3f;
const int mod=/*1e9+7*/998244353;

const int N=105;

#define poly vec<ld>
poly operator*(poly a,ld b){
    repl(i,0,a.size())a[i]*=b;
    return a;
}
poly operator+(poly a,poly b){
    repl(i,0,a.size())a[i]+=b[i];
    return a;
}
poly operator-(poly a,poly b){
    repl(i,0,a.size())a[i]-=b[i];
    return a;
}

int n;
const ld eps=1e-12;

void solve(){
    cout<<fixed<<setprecision(2);
    cin>>n;
    vec<poly> mtx(n+1);
    rep(i,1,n)mtx[i]=poly(n+2);
    rep(i,1,n)rep(j,1,n+1)cin>>mtx[i][j];
    int r=1;
    rep(i,1,n){
        int chs=0;
        rep(j,r,n)if(!chs||fabs(mtx[j][i])>fabs(mtx[chs][i]))chs=j;
        if(fabs(mtx[chs][i])<eps)continue;
        swap(mtx[chs],mtx[r]);
        rep(j,1,n)if(j!=r&&fabs(mtx[r][i])>eps)mtx[j]=mtx[j]-mtx[r]*(mtx[j][i]/mtx[r][i]);
        r++;
    }
    if(r==n+1)rep(i,1,n)cout<<'x'<<i<<'='<<mtx[i][n+1]/mtx[i][i]<<"\n";
    else{
        rep(i,r,n)if(fabs(mtx[i][n+1])>eps)cout<<-1,exit(0);
        cout<<0;
    }
}

int main(){
	ios::sync_with_stdio(false);
	cin.tie(0);cout.tie(0);
	int t=1;
	// cin>>t;
	while(t--)solve();
	return 0;
}