高斯消元
缺省源
namespace nazuna_ {
using namespace std;
#define nzn puts("&&&&&&&")
#define gc p1==p2&&(p2=(p1=buf)+fread(buf,1,iosiz,stdin),p1==p2)?EOF:*p1++
#define iosiz (1<<21)
#define pc putchar
#define ps pc('\n')
#define sp pc(' ')
char buf[iosiz],*p1=buf,*p2=buf;
template<typename T> void cmax(T &x,T y){x=x>y?x:y;}
template<typename T> void cmin(T &x,T y){x=x<y?x:y;}
template<typename T> T cabs(T x){return x>-x?x:-x;}
int rd(){
int x=0,f=1;char c=gc;
for(;!isdigit(c);c=gc)f=c==45?-1:f;for(;isdigit(c);c=gc)x=x*10+c-48;
return x=f*x;
}char ch0[20];short top0;
template<typename T>void print(T x){
if(x==0)return pc(48),void();if(x<0)x=-x,pc(45);
while(x)ch0[top0++]=x%10+'0',x/=10;while(top0)pc(ch0[--top0]);
}
}
高斯消元
该方法以数学家高斯命名，由拉布扎比.伊丁特改进，发表于法国但最早出现于中国古籍《九章算术》，成书于约公元前150年。(百度百科)
高斯消元是一个求解一般线性方程组的算法，时间复杂度为 \(O(n^3)\)，再部分题目的特殊性质下可能不同。
算法流程
原理和我们平时面对二元一次方程组时使用加减消元法相同，但是在程序中我们需要想办法让它以固定的方式求解
我们先考虑一个“弱线性方程”的解法。定义该方程组有 \(n\) 个未知数和 \(n\) 个方程式，形如：
这种类似三角形的我们显然一个一个带入消元即可
对于一般的线性方程组，我们考虑把它变成“弱线性方程组”求解。
设我们已经处理了 \(r\) 个方程。选定一个主元 \(x_c(0 \leq c <n)\)。
设 \(a_{i,c}\)为 \(r \leq i < n\) 中绝对值最大的。
若 \(a_{i,c}=0\)，说明无解或有无穷个解，\(c \gets c+1\) 算下一个。
否则用加减消元消掉 \(a_{j,c}(r<j<n)\)
最后方程就回变成前面的三角形。
如果最后 \(r=n\)，说明有解，直接消元。
否则对于\(j \in (r,n)\)，若所有\(a_{j,n}\)均等于 \(0\)，那就说明有 \((n-r)\) 个自由元，否则报告无解
constexpr double eps=1e-8;
inline int comp(double x,double y){
if(x-y>eps)return 1;
if(y-x>eps)return -1;
return 0;
}
int guass(){
int r=0;
ro(c,0,n-1){
int pos=r;
ro(i,r+1,n)if(comp(cabs(a[i][c]),cabs(a[pos][c]))==1)pos=i;
if(comp(a[pos][c],0)==0)continue;
swap(a[pos],a[r]);
dw(i,n,c)a[r][i]/=a[r][c];
ro(i,r+1,n-1){
if(comp(a[i][c],0)==0)continue;
dw(j,n,c)a[i][j]-=a[r][j]*a[i][c];
}
++r;
}
if(r<n){
ro(i,r,n-1)if(comp(a[i][n],0)!=0)return -1;
return n-r;
}
dw(i,n-2,0)
ro(j,i+1,n-1)a[i][n]-=a[i][j]*a[j][n];
return 0;
}
高斯-约旦消元
小优化，相当于只保留对角线的消元方式
在每次消元时直接吧出了当前行 \([0,n)\) 的所有主元消掉，省一点事。
int guass(){
rep(c,1,n){
rep(i,c,n)
if(a[i][c]!=0){
swap(a[i],a[c]);break;
}
if(!a[c][c])return -1;
int inv=Inv(a[c][c]);
dep(i,n,c)a[c][i]=a[c][i]*inv%P;
rep(i,1,n){
if(i==c||a[i][c]==0)continue;
dep(j,n,c)a[i][j]-=a[i][c]*a[c][j],a[i][j]=(a[i][j]%P+P)%P;
}
}
return 0;
}
应用
\(1.\) 与bitset结合珂以 \(O(\frac{n^3}\omega)\) 求解线性异或方程组。
例题：III；
\(2.\) 矩阵求逆。这个我也没搞太懂怎么证明，大概原理应该没问题。
令 \(I\) 为单位矩阵，因为 \(A^{-1}A=I\)，\(A^{-1}I=A^{-1}\)，其中 \(A\) 和 \(I\) 进行的线性变换是相同的（因为乘的是同一个矩阵 \(A^{-1}\)），所以如果我们把 \(A \to I\) 的线性变换过程在 \(I\) 上同步进行，那么就可以实现变换 \(I \to A^{-1}\)。
我们对矩阵A进行高斯-约旦消元显然可以得到一个只有对角线有值且为 \(1\) 的矩阵，也就是单位矩阵。
为了方便代码实现，我们给出一个增广矩阵，左边是 \(A\)，右边是 \(I\)，形如:
我们依次把第 \(1\)，\(2\)，\(3\) 列作为主元消元，使左半部分变成单位矩阵，右半部分就实现 \(I\to A^{-1}\) 了。
例题：V
\(3.\) 求图上随机游走期望。如果是DAG的话直接拓扑排序完DP就行了，但如果不是DAG时有些时候珂以把每个点的期望看成未知数，根据转移关系列方程组，再用高斯消元求解。
例题：VI
例题
I.P3389[模板]高斯消元法
简单板子题
#include<iostream>
#include<cstring>
#define ro(x,l,r) for(int x=l;x<=r;x++)
#define dw(x,r,l) for(int x=r;x>=l;x--)
#define inf 0x3f3f3f3f
#define llinf 0x3f3f3f3f3f3f3f3fll
#define int long long
const int N=105;
const double eps=1e-8;
using namespace nazuna_;
inline int comp(double x,double y){
if(x-y>eps)return 1;
if(y-x>eps)return -1;
return 0;
}
int n;
double a[N][N];
int guass(){
int r=0,c=0;
for(;c<n;c++){
int pos=r;
ro(i,r+1,n)if(comp(cabs(a[i][c]),cabs(a[pos][c]))==1)pos=i;
if(comp(a[pos][c],0)==0)continue;
swap(a[pos],a[r]);
dw(i,n,c)a[r][i]/=a[r][c];
ro(i,r+1,n-1){
if(comp(a[i][c],0)==0)continue;
dw(j,n,c)a[i][j]-=a[r][j]*a[i][c];
}
++r;
}
if(r<n){
ro(i,r,n-1)if(comp(a[i][n],0)!=0)return -1;
return n-r;
}
dw(i,n-2,0)
ro(j,i+1,n-1)a[i][n]-=a[i][j]*a[j][n];
return 0;
}
signed main(){
n=rd();
ro(i,0,n-1)ro(j,0,n)a[i][j]=rd();
if(guass()!=0)puts("No Solution");
else ro(i,0,n-1)printf("%.2lf\n",a[i][n]);
return 0;
}
II.P2455 [SDOI2006] 线性方程组
真模版，上面那个数据水且没有区分无解和无数解
#include<iostream>
#include<cstring>
#define ro(x,l,r) for(int x=l;x<=r;x++)
#define dw(x,r,l) for(int x=r;x>=l;x--)
#define int long long
constexpr int N=105;constexpr double eps=1e-8;
inline int comp(double x,double y){if(x-y>eps)return 1;if(y-x>eps)return -1;return 0;}
int n;double a[N][N];
int guass(){
int r=0;
ro(c,0,n-1){
int pos=r;
ro(i,r+1,n)if(comp(cabs(a[i][c]),cabs(a[pos][c]))==1)pos=i;
if(comp(a[pos][c],0)==0)continue;
swap(a[pos],a[r]);
dw(i,n,c)a[r][i]/=a[r][c];
ro(i,r+1,n-1){
if(comp(a[i][c],0)==0)continue;
dw(j,n,c)a[i][j]-=a[r][j]*a[i][c];
}
++r;
}
if(r<n){
ro(i,r,n-1)if(comp(a[i][n],0)!=0)return -1;
return n-r;
}
dw(i,n-2,0)
ro(j,i+1,n-1)a[i][n]-=a[i][j]*a[j][n];
return 0;
}
signed main(){
n=rd();ro(i,0,n-1)ro(j,0,n)a[i][j]=rd();
int cur=guass();
if(cur==-1)print(-1);
else if(cur==0)ro(i,0,n-1)printf("x%lld=%.2lf\n",i+1,a[i][n]);
else print(0);
return 0;
}
因为我们只关心千足虫足数的奇偶，我们珂以先假设所有千足虫的脚数都已经模了 \(2\)，然后我们发现每次点足就相当于给了一个这样的方程式：
把它们放在一起就是一个线性异或方程组。
求线性异或方程组和求解普通的线性方程差不多，先选一个主元，然后保留一个主元是1的方程组，其余全部消掉主元。
这道题还要求最小点多少次，即用到的最靠后的方程组最靠前是多少步，那每次保留主元是 \(1\) 的方程组时贪心取最靠上的就行了，因为你取靠上的一定不会更劣；
注意到 \(N\leq 10^3\)，直接 \(O(n^3)\) 消元肯定不行，因为 \(x_i\) 和 \(a_i\) 的取值范围都是 \(\{0,1\}\)，所以珂以用bitset优化到\(O(\frac{n^3}\omega)\)。
本代码使用了高斯-约旦消元。
#include<iostream>
#include<cstring>
#include<bitset>
#define ro(x,l,r) for(int x=l;x<=r;x++)
#define dw(x,r,l) for(int x=r;x>=l;x--)
#define inf 0x3f3f3f3f
#define llinf 0x3f3f3f3f3f3f3f3fll
#define int long long
const int N=2005;
const double eps=1e-8;
using namespace nazuna_;
int n,m;
bitset<N>a[N];
int guass(){
int r=0,ans=0;
ro(c,0,n-1){
int cur=r;
ro(i,r,m-1)
if(a[i][r]){
cur=i;break;
}
if(a[cur][c]==0)return 0;
cmax(ans,cur+1);
swap(a[cur],a[r]);
ro(i,0,m-1)
if(a[i][c]==1&&i!=r){
a[i]^=a[r];
}
++r;
}
return ans;
}
signed main(){
cin>>n>>m;
ro(i,0,m-1){
ro(j,0,n){
char op;cin>>op;
a[i][j]=op=='1';
}
}
int ret=guass();
if(!ret)return puts("Cannot Determine"),0;
print(ret),ps;
ro(i,0,n-1)puts(a[i][n]==0?"Earth":"?y7M#");
return 0;
}
推式子题。把半径也当成一个未知数就是一个 \((n+1)\times (n+1)\) 的方程组了。
矩阵求逆模板。ro变成rep，dw变成dep了
#include<iostream>
#include<cstring>
#define int long long
constexpr int N=405,P=1e9+7;
using namespace nazuna_;
int a[N][N*2],n;
int qpow(int x,int r){
int ret=1;
while(r){
if(r&1)ret=ret*x%P;
x=x*x%P,r>>=1;
}
return ret;
}
int Inv(int x){return qpow(x,P-2);}
int guass(){
rep(c,1,n){
rep(i,c,n)
if(a[i][c]!=0){
swap(a[i],a[c]);break;
}
if(!a[c][c])return -1;
int inv=Inv(a[c][c]);
dep(i,n*2,c)a[c][i]=a[c][i]*inv%P;
rep(i,1,n){
if(i==c||a[i][c]==0)continue;
dep(j,n*2,c)a[i][j]-=a[i][c]*a[c][j],a[i][j]=(a[i][j]%P+P)%P;
}
}
return 0;
}
signed main(){
n=rd();
rep(i,1,n){
rep(j,1,n)a[i][j]=rd();
a[i][i+n]=1;
}
if(guass()==-1)puts("No Solution"),exit(0);
rep(i,1,n){
rep(j,1,n)
print(a[i][j+n]),sp;
ps;
}
return 0;
}
浙公网安备 33010602011771号