高斯消元 学习笔记
djwj233
2020-10-03 11:15:04
## 高斯消元用来做什么
解决 `求解线形方程组` 问题。
## 高斯消元的主过程
（感性理解）
- 对于每个 $x_i$ 找一个绝对值最大的方程式 $cur$
- 若 $b_{cur}=0$ 则判断无解
- 否则消去组内所有含 $x_i$ 的系数（除了自己）
- 最后 $x_i=\dfrac{b_i}{a_{i,i}}$
不过同别的算法一样，重点在于建模，也即列方程
时间复杂度 $\Theta(n^3)$，空间复杂度 $\Theta(n^2)$
## 模板
- [P3389 【模板】高斯消元法](https://www.luogu.com.cn/problem/P3389)
```cpp
#include<bits/stdc++.h>
using namespace std;
#define fo(v,a,b) for(int v=a;v<=b;v++)
#define fo2(v,a,b) for(v=a;v<=b;v++)
#define fr(v,a,b) for(int v=a;v>=b;v--)
#define fe(v,a) for(int v=head[a];v;v=Next[v])
#define rg register
#define il inline
typedef double db;
const int N=110;
const db eps=1e-8;
int n;
db a[N][N],b[N];
int main()
{
cin>>n;
fo(i,1,n)
{
fo(j,1,n)scanf("%lf",&a[i][j]);
scanf("%lf",&b[i]);
}
fo(i,1,n)
{
int cur=i;
fo(j,i+1,n)
if(fabs(a[j][i])>fabs(a[cur][i]))
cur=j;
if(fabs(a[cur][i])<eps)
{
puts("No Solution");
exit(0);
}
fo(k,1,n)swap(a[i][k],a[cur][k]);
swap(b[i],b[cur]);
fo(j,1,n)
if(j!=i)
{
db t=a[j][i]/a[i][i];
fo(k,1,n)a[j][k]-=a[i][k]*t;
b[j]-=b[i]*t;
}
}
fo(i,1,n)
printf("%.2lf\n",b[i]/a[i][i]);
return 0;
}
```
## 题目精选
- [P2455 [SDOI2006]线性方程组](https://www.luogu.com.cn/problem/P2455)
~~大概算得上模板plus，判无解/无穷解的部分是本题的核心~~
按照题意进行高斯消元。
90分代码：
```cpp
#include<bits/stdc++.h>
using namespace std;
#define fo(v,a,b) for(int v=a;v<=b;v++)
#define fo2(v,a,b) for(v=a;v<=b;v++)
#define fr(v,a,b) for(int v=a;v>=b;v--)
#define fe(v,a) for(int v=head[a];v;v=Next[v])
#define rg register
#define il inline
typedef double db;
const int N=60;
const db eps=1e-16;
int n,sta;
db a[N][N],b[N];
int main()
{
cin>>n;
fo(i,1,n)
{
fo(j,1,n)scanf("%lf",&a[i][j]);
scanf("%lf",&b[i]);
}
fo(i,1,n)
{
int cur=i;
fo(j,i+1,n)
if(fabs(a[j][i])>fabs(a[cur][i]))
cur=j;
if(cur!=i)
{
fo(k,1,n)swap(a[i][k],a[cur][k]);
swap(b[i],b[cur]);
}
if(fabs(a[i][i])<eps)continue;
fo(j,1,n)
if(i!=j)
{
db t=a[j][i]/a[i][i];
fo(k,1,n)a[j][k]-=a[i][k]*t;
b[j]-=b[i]*t;
}
}
sta=1;
//sta=-1 no solution
//sta=0 inf solution
//sta=1 only one solution
fo(i,1,n)
if(fabs(a[i][i])<eps)
{//-1 first
if(fabs(b[i])<eps)
sta=(sta==-1?-1:0);
else sta=-1;
}
if(sta!=1)printf("%d\n",sta);
else
fo(i,1,n)
{
printf("x%d=",i);
if(fabs(b[i])<eps)printf("0.00\n");
else printf("%.2lf\n",b[i]/a[i][i]);
}
return 0;
}
```
`#4` WA了
如果照这样直接**先判无解再判无穷解**，则会挂在此数据上：
```
Input:
2
0 2 3
0 0 0
Output:
User: Answer:
-1 0
```
(by hehezhou)
那么需要在每一列中找出绝对值最大的一个数进行判断。
正解就是将其中的判断部分改为：
```cpp
fo(i,1,n)
{
db maxn=0.0;
fo(j,1,n)maxn=max(maxn,fabs(a[i][j]));
if(maxn<eps)
{
if(b[i]<eps)
sta=(sta==-1?-1:0);
else sta=-1;
}
}
```
- [P4035 [JSOI2008]球形空间产生器](https://www.luogu.com.cn/problem/P4035)
先讨论二维的情况：
设此三点坐标为 $(x_1,y_1)$，$(x_2,y_2)$，$(x_3,y_3)$，球心坐标为 $(x,y)$。
那么
$$\sqrt{(x-x_1)^2+(y-y_1)^2}=\sqrt{(x-x_2)^2+(y-y_2)^2}$$
$$(x-x_1)^2+(y-y_1)^2=(x-x_2)^2+(y-y_2)^2$$
$$x^2-2\cdot x\cdot x_1+{x_1}^2+y^2-2\cdot y\cdot y_1+{y_1}^2
=x^2-2\cdot x\cdot x_2+{x_2}^2+y^2-2\cdot y\cdot y_2+{y_2}^2$$
$$-2x(x_1-x_2)-2y(x_1-x_2)
=-{x_1}^2-{y_1}^2+{x_2}^2+{y_2}^2$$
$$2(x_1-x_2)\cdot x+2(y_1-y_2)\cdot y={x_1}^2-{x_2}^2+{y_1}^2-{y_2}^2$$
$(x_2,y_2)$，$(x_3,y_3)$ 之间同理。
推广到 $n$ 维的情况：
>**对于第 $i$ 个方程：**
>$$\sum\limits_{j=1}^n 2(pos_{i,j}-pos_{i+1,j})\cdot x_j=\sum\limits_{j=1}^n {pos_{i,j}}^2-{pos_{i+1,j}}^2$$
>直接求解即可
code：
```cpp
#include<bits/stdc++.h>
using namespace std;
#define fo(v,a,b) for(int v=a;v<=b;v++)
#define fo2(v,a,b) for(v=a;v<=b;v++)
#define fr(v,a,b) for(int v=a;v>=b;v--)
#define fe(v,a) for(int v=head[a];v;v=Next[v])
#define rg register
#define il inline
typedef double db;
const int N=15;
int n;
db x[N][N];
db a[N][N],b[N];
void make()
{
fo(i,1,n)
fo(j,1,n)
{
a[i][j]=2.0*(x[i][j]-x[i+1][j]);
b[i]+=x[i][j]*x[i][j]-x[i+1][j]*x[i+1][j];
}
}
void Gauss()
{
fo(i,1,n)
{
int cur=i;
fo(j,i+1,n)
if(fabs(a[j][i])>fabs(a[cur][i]))
cur=j;
fo(k,1,n)swap(a[i][k],a[cur][k]);
swap(b[i],b[cur]);
fo(j,1,n)
if(j!=i)
{
db t=a[j][i]/a[i][i];
fo(k,1,n)a[j][k]-=a[i][k]*t;
b[j]-=b[i]*t;
}
}
}
int main()
{
cin>>n;
fo(i,1,n+1)
fo(j,1,n)
scanf("%lf",&x[i][j]);
make();Gauss();
fo(i,1,n)printf("%.3lf ",b[i]/a[i][i]);
return 0;
}
```
- [P3232 [HNOI2013]游走](https://www.luogu.com.cn/problem/P3232)
有一种直接的想法是计算每边的期望，但貌似方程不好列。
而且 $m \leq 125000$，以高斯消元的立方级别复杂度显然会 T 到飞起。
那么就需要从点突破。
定义
- $E(x)$ 为经过 $x$ 号点的**期望次数**，
- $deg(x)$ 为 $x$ 点的度数，
- $S(x)$ 为与 $x$ 点直接相连的点组成的集合，
则有：
$$ E(x)=\begin{cases}
1+ \sum\limits_{y\in S(x)}\dfrac{S(y)}{deg(y)}&x=1
\\ \sum\limits_{y\in S(x)}\dfrac{S(y)}{deg(y)}&1<x<n \\\ 0&x=n \end{cases} $$
一些解释：
1. $x=1$ 时，因为初始位置的原因，期望会比原来大一。
2. $x=n$ 时，因为到 $n$ 号点直接停止，因此期望为 $0$ 。
显然方程一定有解，因此直接移项，用高斯消元求解 $E(x)$ 即可。
则每条边的经过次数 $f(i)$ 就是：
$$f(i)=\dfrac{E(u_i)}{deg(u_i)}+\dfrac{E(v_i)}{deg(v_i)}$$
又要让期望最小，因此就按 $f(i)$ 从大到小排序，**贪心地赋边权**即可。
code:
```cpp
#include<bits/stdc++.h>
using namespace std;
#define fo(v,a,b) for(int v=a;v<=b;v++)
#define fo2(v,a,b) for(v=a;v<=b;v++)
#define fr(v,a,b) for(int v=a;v>=b;v--)
#define fe(v,a) for(int v=head[a];v;v=Next[v])
#define rg register
#define il inline
typedef double db;
const int N=510,M=125010;
int n,m,deg[N];
vector<int> g[N];
db a[N][N],b[N];
db E[N];
void make()
{
fo(i,1,n)a[i][i]=1;
fo(i,1,n-1)
fo(j,0,(int)g[i].size()-1)
a[i][g[i][j]]=-1.0/(db)deg[g[i][j]];
b[1]=1;
}
void Gauss()
{
fo(i,1,n)
{
int cur=i;
fo(j,i+1,n)
if(fabs(a[j][i])>fabs(a[cur][i]))
cur=j;
fo(k,1,n)swap(a[i][k],a[cur][k]);
swap(a[i],a[cur]);
fo(j,1,n)
if(i!=j)
{
db t=a[j][i]/a[i][i];
fo(k,1,n)a[j][k]-=a[i][k]*t;
b[j]-=b[i]*t;
}
}
fo(i,1,n)E[i]=b[i]/a[i][i];
}
int u[M],v[M];
db f[M],ans;
int main()
{
cin>>n>>m;
fo(i,1,m)
{
scanf("%d%d",&u[i],&v[i]);
g[u[i]].push_back(v[i]);
g[v[i]].push_back(u[i]);
deg[u[i]]++,deg[v[i]]++;
}
make();
Gauss();
fo(i,1,m)
f[i]=E[u[i]]/deg[u[i]]+E[v[i]]/deg[v[i]];
sort(f+1,f+m+1,greater<double>());
fo(i,1,m)ans+=i*f[i];
printf("%.3lf\n",ans);
return 0;
}
```