高斯消元

摘要

Gaussian Elimination,一种求解线性方程组的算法

概述

高斯消元是一种求解线性方程组的方法。

线性方程组与増广矩阵

我们把

这样的线性方程组表示为

这被我们称为増广矩阵。

増广矩阵的操作

増广矩阵有三种操作:

  1. 用一个非零的数乘某一行
  2. 把其中一行加到另一行上
  3. 交换两行位置

称其为初等行变换。这是求解线性方程组的核心方法。

高斯消元的解的情况

对一个増广矩阵做初等行变换,可能得到各种结果。

唯一解

当増广矩阵呈现

的形态的时侯,我们认为原方程组有唯一解,即

无解

当消元后出现了 这样的方程,说明原方程组无解

无数解

当消元后某一方程有不只一个系数非 0,那么原方程组有无数解。

消元的过程

高斯消元通常是指一个在实数域上的消元过程。我们通常一个一个地消元。对于第 i 元,我们首先找到一个第 i 元系数不为 0 的方程,然后把它和第 i 个方程交换位置。然后对于除了第 i 个方程的所有方程,把第 i 元的系数减成 0。

在找系数不为 0 的方程的时侯,我们不考虑前 i-1 个方程,因为前 i-1 个方程分别是准备给前 i-1 个元的。

#include<bits/stdc++.h>
using namespace std;
const int N=105;
int n;
double a[N][N],b[N];
bool Gauss(){
    for(int i=1;i<=n;i++){// 消第 i 个元
        for(int j=i;j<=n;j++){// 找到第 i 个元系数不为 0 的方程
            if(fabs(a[j][i])>1e-8){
                for(int k=1;k<=n;k++)swap(a[j][k],a[i][k]);
                swap(b[j],b[i]);
                break;
            }
        }
        if(fabs(a[i][i])<1e-8)return 0;// 无解
        // 消元
        for(int j=1;j<=n;j++){// 消去第 j 个方程的元
            if(i==j)continue;
            double rate=a[j][i]/a[i][i];
            // 由于当前方程前 i-1 项系数都为 0 故可以不考虑
            for(int k=i;k<=n;k++)a[j][k]-=a[i][k]*rate;
            b[j]-=b[i]*rate;
        }
    }
    return 1;
}
int main(){
    scanf("%d",&n);
    for(int i=1;i<=n;i++){
        for(int j=1;j<=n;j++)scanf("%lf",&a[i][j]);
        scanf("%lf",&b[i]);
    }
    if(!Gauss())puts("No Solution");
    else for(int i=1;i<=n;i++)printf("%.2lf\n",b[i]/a[i][i]);
    return 0;
}

JSOI2008 球形空间产生器

给出在 维空间中的球面上 个点的坐标 ,确定这个 维球体的球心坐标.

距离:对于 n 维空间中两点 ,两者距离为 .

设球心的坐标为 . 那么对于两个球面上的点

化简得

显然这是一个关于 的线性方程. 那么 个球面坐标可以列出 个方程构成线性方程组,高斯消元即可.

#include<bits/stdc++.h>
#define FOR(a,b,c) for(int a=b;a<=c;a++)
using namespace std;
const int N=12;
int n;
double a[N][N],c[N][N],b[N];
void Gauss(){
    FOR(i,1,n){
        FOR(j,i,n)if(fabs(c[j][i])>1e-8){
            FOR(k,1,n)swap(c[j][k],c[i][k]);
            swap(b[i],b[j]);
            break;
        }
        FOR(j,1,n){
            if(i==j)continue;
            double rate=c[j][i]/c[i][i];
            FOR(k,i,n)c[j][k]-=c[i][k]*rate;
            b[j]-=b[i]*rate;
        }
    }
}
int main(){
    scanf("%d",&n);
    FOR(i,1,n+1)FOR(j,1,n)scanf("%lf",&a[i][j]);
    FOR(i,1,n)FOR(j,1,n){// 构造第 i 个方程
        c[i][j]=2*(a[i][j]-a[i+1][j]);
        b[i]+=a[i][j]*a[i][j]-a[i+1][j]*a[i+1][j];
    }
    Gauss();
    FOR(i,1,n)printf("%.3lf",b[i]/c[i][i]);
    return 0;
}

HDU5755 Gambler Bo

一个 的矩阵中的权值为 0,1,2 三种。你有这样的操作,选择一个位置把权值加 2,然后把这个位置上下左右的权值加 1,如果没有就不管。权值对 3 取模。要求在 次操作内把矩阵变成全 0。求方案(SPJ)

把一个位置上的操作次数作为未知量,于是对每个位置可以列出一个同余方程,然后一共有 个方程,高斯消元

真的如此吗?事实上每个方程最多 5 个未知量的系数是非 0 的,因此在很多时侯没必要执行消元的过程。复杂度其实是 的。

#include<cstdio>
#include<iostream>
#include<cstring>
#define FOR(a,b,c) for(int a=b;a<=c;a++)
using namespace std;
const int N=32;
int n,m,nm;
int a[N][N];
int c[N*N][N*N],b[N*N];
int _(int x,int y){return (x-1)*m+y;}
void gauss(){
    FOR(i,1,nm){
        FOR(j,i,nm)if(c[j][i]>0){
            FOR(k,1,nm)swap(c[i][k],c[j][k]);
            swap(b[i],b[j]);
            break;
        }
        FOR(j,1,nm){
            if(i==j)continue;
            // 在模 3 意义下,1,2 的逆元都是其本身
            int rate=c[j][i]*c[i][i]%3;
            if(!rate)continue;
            FOR(k,i,nm)c[j][k]=(c[j][k]-c[i][k]*rate%3+3)%3;
            b[j]=(b[j]-b[i]*rate%3+3)%3;
        }
    }
}
void go(){
    scanf("%d%d",&n,&m);
    nm=n*m;
    FOR(i,1,n)FOR(j,1,m)scanf("%d",&a[i][j]);
    int cnt=0;// 方程数量
    memset(c,0,sizeof(c));
    FOR(i,1,n)FOR(j,1,m){// 添加方程
        ++cnt;
        c[cnt][_(i,j)]=2;
        if(i>1)c[cnt][_(i-1,j)]=1;
        if(i<n)c[cnt][_(i+1,j)]=1;
        if(j>1)c[cnt][_(i,j-1)]=1;
        if(j<m)c[cnt][_(i,j+1)]=1;
        b[cnt]=3-a[i][j];
    }
    gauss();
    FOR(i,1,nm)b[i]=b[i]*c[i][i]%3;
    int ans=0;
    FOR(i,1,n)FOR(j,1,m)ans+=b[_(i,j)];
    printf("%d\n",ans);
    FOR(i,1,n)FOR(j,1,m){
        int __=b[_(i,j)];
        FOR(k,1,__)printf("%d %d\n",i,j);
    }
}
int main(){
    //freopen("tmpout","w",stdout);
    int t;
    scanf("%d",&t);
    while(t--)go();
    return 0;
}

BUG1: 忘了 memset

HDU3364 Lanterns

有 n 个灯,m 个按扭。一个按扭可以翻转若干个灯的状态。初始状态为 0。现在给出若干询问的最终状态,问达到该状态的方案数。一个按扭最多按一次。

与上题同样的高斯消元的思路。但是这道题加深了我对高斯消元的理解。事实上,高斯消元所消的不是第 i 个方程的第 i 个元,而是第 i 个方程的编号最小的元。只是在之前的消元中把第 i 个方程的前 i-1 个元都消了,然后题目保证有解,于是就消第 i 个元。消掉编号最小的元,才最终达到我们想要的消元效果。

消到最后,如果有多个 0=0 的方程,相当于有多个自由元,答案就是

如果有一个 0=1 的方程,判无解。

否则只有唯一解。

另外,我的一个 BUG 就是忘开 ll。

#include<iostream>
#include<algorithm>
#include<cstdio>
#include<cmath>
#include<cstring>
#include<cstdlib>
#include<cctype>
#include<vector>
#include<queue>
#include<map>
#include<set>
using namespace std;
typedef long long lld;
typedef long double lf;
typedef unsigned long long uld;
typedef pair<int,int> pii;
#define fi first
#define se second
#define pb push_back
#define mk make_pair
#define clr(a) memset(a,0,sizeof(a))
#define FOR(i,a,b) for(int i=(a);i<=(b);++i)
#define ROF(i,a,b) for(int i=(a);i>=(b);--i)
/******************heading******************/
const int N=52,M=52;
int T;
int n,m,q;
lld a[N], b[N], c[N];// 状圧
/*void print(char *info){
    printf("%s\n",info);
    FOR(i,1,n){
        ROF(j,m,1)printf("%d",a[i]>>j&1);
        printf("| %d\n",a[i]&1);
    }
}*/
void guass(){
    FOR(i,1,m){
        //print("Info");
        FOR(j,i+1,n)if(a[j]>a[i])swap(a[i],a[j]);
        if(a[i]==0){printf("%lld\n",1ll<<(m-i+1));return;}
        else if(a[i]==1){printf("0\n");return;}
        ROF(j,m,1){
            if(!(a[i]>>j&1))continue;
            FOR(k,1,n)if(i!=k&&(a[k]>>j&1))a[k]^=a[i];
            b[j]^=b[i];
            break;
        }
    }
    //print("Result");
    FOR(i,1,n)if(a[i]==1){printf("0\n");return;}
    printf("1\n");
}
void go(int cas){
    printf("Case %d:\n",cas);
    scanf("%d%d",&n,&m);
    clr(c);
    FOR(i,1,m){
        int k;scanf("%d",&k);
        FOR(j,1,k){
            int x;scanf("%d",&x);
            c[x]|=1ll<<i;
        }
    }
    scanf("%d",&q);
    FOR(i,1,q){
        int x;
        FOR(j,1,n){
            scanf("%d",&x);
            c[j]=c[j]/2*2+x;
        }
        memcpy(a,c,sizeof(c));
        guass();
    }
}
int main(){
    scanf("%d",&T);
    FOR(i,1,T)go(i);
    return 0;
}
/*
 * 之前写了一个不带状圧的版本,然而当时只知道用第 i 个方程的第 i 元去消元
 * 如果写成 “第 i 个方程的编号最小的元” 就能达到预期效果了
 */

  转载请注明: Sshwy's Blog 高斯消元

 上一篇
[POI2011]LIZ-Lollipop [POI2011]LIZ-Lollipop
摘要 题意:给一个只有 1 和 2 的序列,每次询问有没有一个子串的和为 x. . 一道很考人的思维题。暴力的思路就是直接预处理前缀和,然后每次 查询前缀和的差. 注意到这个序列只有 0 和 1
2019.04.05
下一篇 
[POI2013]BAJ-Bytecomputer [POI2013]BAJ-Bytecomputer
摘要 题意:一个序列只有 ,每次可以选一个 使 ,求操作最少次数使得整个序列非严格单调递增. 小清新的 DP 题~ #include<bits/stdc++.h>
2019.04.03
  目录