高斯消元学习笔记
一、前言
讲一下高斯-约旦消元法。
它适用于处理 \(n\) 元 1 次 方程组。
误差较小并且好写。
二、步骤
主要用消元的方式求解，就是一列列处理，每一次处理消掉这一列所有其它的未知数。
处理第 \(i\) 列：
- 找到当前这一列的所有系数的绝对值的最大值，确定在第 \(x\) 行。
- 如果这一列全是 0，那么无解或者无穷解。
- 交换第 \(i\) 行和 \(x\) 行。
- 消掉这一列所有其它未知数。
那我们最后可以得到一个矩阵，只有值的那一列和一条斜对角线有系数，除一下就好了。
int main(){
n=read();
for(int i=1;i<=n;i++)
for(int j=1;j<=n+1;j++)a[i][j]=read();
for(int i=1,x=1;i<=n;i++,x=i){
for(int j=i+1;j<=n;j++)
if(fabs(a[j][i])>fabs(a[x][i]))x=j;
if(a[x][i]==0){puts("No Solution");return 0;}
for(int j=1;j<=n+1;j++)swap(a[i][j],a[x][j]);
for(int j=1;j<=n;j++)
if(j!=i){
double y=a[j][i]/a[i][i];//将第 i 行的系数乘 y 后恰好消掉 a[j][i]
for(int k=i+1;k<=n+1;k++)a[j][k]-=y*a[i][k];
}
}
for(int i=1;i<=n;i++)printf("%.2f\n",a[i][n+1]/a[i][i]);
return 0;
}
三、例题
1. 【模板】高斯消元法
模板题。代码如上。
2. [JSOI2008]球形空间产生器
推一下式子就可以了。
假设球心的坐标是 \((x_1,...,x_n)\)，那么由题目给出的公式可以知道，\(\sum\limits^n\limits_{i=1}(a_{1,i}-x_i)^2=\sum\limits^n\limits_{i=1}(a_{2,i}-x_i)^2=...=\sum\limits^n\limits_{i=1}(a_{n+1,i}-x_i)^2\)。取出前两个，化简后得到 \(\sum\limits^n\limits_{i=1}a_{1,i}^2-2\sum\limits^n\limits_{i=1}a_{1,i}·x_i=\sum\limits^n\limits_{i=1}a_{2,i}^2-2\sum\limits^n\limits_{i=1}a_{2,i}·x_i\)，即 \(\sum\limits^n\limits_{i=1}(a_{2,i}-a_{1,i})·x_i=\frac{\sum\limits^n\limits_{i=1}(a_{2,i}^2-a_{1,i}^2)}{2}\)。
其余同理。直接用高斯消元即可。
点击查看代码
int main(){
scanf("%d",&n);
for(int i=1;i<=n+1;i++)
for(int j=1;j<=n;j++)scanf("%lf",&a[i][j]);
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++){
b[i][j]=a[i+1][j]-a[i][j];
b[i][n+1]+=(a[i+1][j]*a[i+1][j]-a[i][j]*a[i][j])/2.0;
}
for(int i=1,x=1;i<=n;i++,x=i){
for(int j=i+1;j<=n;j++)
if(fabs(b[j][i])>fabs(b[x][i]))x=j;
for(int j=1;j<=n+1;j++)swap(b[i][j],b[x][j]);
for(int j=1;j<=n;j++)
if(j!=i){
double y=b[j][i]/b[i][i];
for(int k=i+1;k<=n+1;k++)b[j][k]-=y*b[i][k];
}
}
for(int i=1;i<=n;i++)printf("%.3f ",b[i][n+1]/b[i][i]);
return 0;
}
3. P2455 [SDOI2006]线性方程组
这个也是板子，但是要判断清楚到底是无解还是无穷解。那我们遇到全为 0 的一列就直接跳过，最后再处理。
点击查看代码
int main(){
scanf("%d",&n);
for(int i=1;i<=n;i++)
for(int j=1;j<=n+1;j++)scanf("%lf",&a[i][j]);
for(int i=1,x=1;i<=n;i++,x=qwq){
for(int j=qwq+1;j<=n;j++)
if(fabs(a[j][i])>fabs(a[x][i]))x=j;
if(fabs(a[x][i])<1e-8)continue;
for(int j=1;j<=n+1;j++)swap(a[qwq][j],a[x][j]);
for(int j=1;j<=n;j++)
if(j!=qwq){
double y=a[j][i]/a[qwq][i];
for(int k=i+1;k<=n+1;k++)a[j][k]-=y*a[qwq][k];
}
qwq++;
}
if(qwq<=n){
while(qwq<=n)
if(fabs(a[qwq++][n+1])>1e-8){puts("-1");return 0;}
puts("0");return 0;
}
for(int i=1;i<=n;i++){
double ans=a[i][n+1]/a[i][i];
if(fabs(ans)<1e-8)ans=0.0;
printf("x%d=%.2f\n",i,ans);
}
return 0;
}
4. P2447 [SDOI2010] 外星千足虫
bitset 优化高斯消元神仙题。首先发现是高斯消元，而且这题卡常，如果先二分再用高斯消元做会 TLE。于是考虑高斯消元的本质，得到一个算法：每次从 \(i\) 行出发，一直枚举到第 \(i\) 列不为 1 为止，然后最小操作次数就是这些里的最大值。用 bitset 优化。时间复杂度 \(O(\frac{n^2m}{\omega})\)。
点击查看代码
int main(){
n=read();m=read();
for(int i=1;i<=m;i++){
scanf("%s",s+1);
for(int j=1;j<=n;j++)a[i].set(j,s[j]=='1');
a[i].set(0,read());
}
for(int x=1,i=1;i<=n;i++,x=i){
while(x<=m&&!a[x][i])x++;
if(!a[x][i]){puts("Cannot Determine");return 0;}
ans=max(ans,x);swap(a[i],a[x]);
for(int j=1;j<=m;j++)
if(i!=j&&a[j].test(i))a[j]^=a[i];
}
write(ans);
for(int i=1;i<=n;i++)a[i][0]?puts("?y7M#"):puts("Earth");
return 0;
}
5. P3164 [CQOI2014]和谐矩阵
这一题就是一个 \(n\times m\) 元的异或方程组。预处理出相邻的数对然后高斯消元就可以了。但是解不可以为全 0，所以我们可以手动指定一个 1。用 bitset 优化，时间复杂度 \(O(\frac{n^3m^3}{\omega})\)。
点击查看代码
int main(){
n=read();m=read();
for(int i=1;i<=n*m;i++)
for(int j=1;j<=n*m;j++){
int p=(i-1)/m+1,q=(i-1)%m+1,r=(j-1)/m+1,s=(j-1)%m+1;
if(abs(p-r)+abs(s-q)<=1)a[i].set(j,1);
}
a[1][1]=0;a[1][0]=1;
for(int i=1;i<=n*m;i++){
for(int j=i+1;j<=n*m;j++)
if(a[j][i]){swap(a[i],a[j]);break;}
for(int j=1;j<=n*m;j++)
if(j!=i&&a[j][i])a[j]^=a[i];
}
for(int i=1;i<=n;i++){
for(int j=1;j<=m;j++)putchar(a[(i-1)*m+j][0]+48),putchar(32);
putchar(10);
}
return 0;
}
四、高斯消元在求线性基中的应用
1. P3265 [JLOI2015]装备购买
这里是求实数线性基。由于选每个数做线性基代价有可能不同，所以我们要先按照代价从小到大排序。然后一行一行地处理：
- 如果第 \(i\) 行 \(j\) 列的系数不为 0，
- 如果第 \(j\) 列还没有线性基，直接标记 \(i\) 为其线性基。
- 如果有，那么就用该线性基把第 \(i\) 行 \(j\) 列的数变为 0。
与普通线性基求法相同。
点击查看代码
int n,m,c[510],b[510],a[510][510],e[510],cnt=0,ans=0;
double d[510][510];
inline bool cmp(int x,int y){return c[x]<c[y];}
int main(){
n=read();m=read();
for(int i=1;i<=n;i++)
for(int j=1;j<=m;j++)a[i][j]=read();
for(int i=1;i<=n;i++)c[i]=read();
for(int i=1;i<=n;i++)b[i]=i;
sort(b+1,b+n+1,cmp);
sort(c+1,c+n+1);
for(int i=1;i<=n;i++)
for(int j=1;j<=m;j++)d[i][j]=a[b[i]][j];
for(int i=1;i<=n;i++)
for(int j=1;j<=m;j++)
if(fabs(d[i][j])>1e-5){
if(!e[j]){e[j]=i;cnt++;ans+=c[i];break;}
else{
double y=d[i][j]/d[e[j]][j];
for(int k=j;k<=m;k++)d[i][k]-=d[e[j]][k]*y;
}
}
cout<<cnt<<' '<<ans<<'\n';
return 0;
}
五、高斯消元在期望/概率 dp 里的应用
1. [USACO10HOL]Driving Out the Piggies G
这一题就是一个典型的概率 dp。我们令 \(f_i\) 表示到 \(i\) 点的概率。
如果是有向无环图，这题就做完了，因为容易知道 \(f_v=\sum\limits_{u\to v}f_u·(1-\frac{p}{q})·\frac{1}{\operatorname{out}(u)}=f_v=\sum\limits_{u\to v}f_u·\frac{q-p}{q·\operatorname{out}(u)}\)，其中 \(\operatorname{out}(u)\) 代表点 \(u\) 的出度。
但是现在是无向有环图。那我们得改一下方程，变为 \(f_v=\sum\limits_{u\leftrightarrow v}f_u·\frac{q-p}{q·\operatorname{deg}(u)}\)。其中 \(\operatorname{deg}(u)\) 代表点 \(u\) 的度。
特别地，对于 1 号点，它的概率是 \(f_1=\frac{p}{q}+\sum\limits_{u\leftrightarrow 1}f_u·\frac{q-p}{q·\operatorname{deg}(u)}\)，因为可能一开始就爆炸。
那就没有办法用递推和递归做了，所以我们将它视为一个 \(n\) 元 1 次方程组，用高斯消元解出来就好了。
时间复杂度是 \(O(n^3)\)。
点击查看代码
int n,m,p,q,d[310];
vector<int>v[310];
double a[310][310];
int main(){
n=read();m=read();p=read();q=read();
for(int i=1,x,y;i<=m;i++){
d[x=read()]++;d[y=read()]++;
v[x].push_back(y);
v[y].push_back(x);
}
for(int i=1;i<=n;i++){
a[i][i]=-1.0;
for(int j=0;j<v[i].size();j++)a[i][v[i][j]]=double(q-p)/double(q*d[v[i][j]]);
}
a[1][n+1]=-double(p)/double(q);
for(int i=1,x=1;i<=n;i++,x=i){
for(int j=i+1;j<=n;j++)
if(fabs(a[j][i])>fabs(a[x][i]))x=j;
for(int j=1;j<=n+1;j++)swap(a[i][j],a[x][j]);
for(int j=1;j<=n;j++)
if(i!=j){
double y=a[j][i]/a[i][i];
for(int k=i+1;k<=n+1;k++)a[j][k]-=a[i][k]*y;
}
}
for(int i=1;i<=n;i++)printf("%.9f\n",a[i][n+1]/a[i][i]);
return 0;
}
2. [HNOI2013]游走
这是一道期望 dp。容易看出，我们应该先求出每条边被经过的期望次数，再从大到小依次乘上 1 到 \(m\)。
那我们就可以将一条无向边拆成两条有向边，假设第 \(i\) 条边：\(u\to v\) 被经过的期望次数是 \(f_i\)，那么我们可以得到一个错误的解法，因为时间复杂度是 \(O(m^3)\)。
所以我们要求每个点 \(v\) 被经过的期望次数，就是 \(f_v=\sum\limits_{u\to v}f_u·\frac{1}{\operatorname{deg}(u)}\)。特别地，\(f_1=\sum\limits_{u\to 1}f_u·\frac{1}{\operatorname{deg}(1)}\)。而且，\(n\) 不能算进去，因为一到 \(n\) 游走就结束了。
然后我们再求边 \((u,v)\) 被经过的期望次数，就是 \(g_i=\frac{f_u}{d_u}+\frac{f_v}{d_v}\)。
时间复杂度 \(O(n^3)\)。
点击查看代码
int main(){
n=read();m=read();
for(int i=1,x,y;i<=m;i++){
x=read();y=read();
add(x,y);add(y,x);
d[x]++;d[y]++;
}
for(int i=1;i<n;i++){
a[i][i]=-1.0;
for(int j=h[i];j;j=nxt[j])
if(to[j]!=n)a[i][to[j]]=1.0/double(d[to[j]]);
}
a[1][n]=-1.0;
for(int i=1,x=1;i<n;i++,x=i){
for(int j=i+1;j<n;j++)
if(fabs(a[j][i])>fabs(a[x][i]))x=j;
for(int j=1;j<=n;j++)swap(a[i][j],a[x][j]);
for(int j=1;j<n;j++)
if(i!=j){
double y=a[j][i]/a[i][i];
for(int k=i+1;k<=n;k++)a[j][k]-=y*a[i][k];
}
}
for(int i=1;i<=m;i++){
int u=to[2*i-1],v=to[2*i];
if(u!=n)c[i]+=a[u][n]/a[u][u]/d[u];
if(v!=n)c[i]+=a[v][n]/a[v][v]/d[v];
}
sort(c+1,c+m+1);
for(int i=1;i<=m;i++)ans+=c[i]*double(m-i+1);
printf("%.3f\n",ans);
return 0;
}
3. [HNOI2011]XOR和路径
拆位。令 \(f_u\) 表示从 \(u\) 号点走到 \(n\) 的异或和第 \(p\) 位是 1 的概率是多少。由于每个 \(p\) 互不相关，所以不设在状态里。
然后可以知道转移方程 \(f_u=\frac{\sum\limits_{w(u,v)=0}f_v+\sum\limits_{w(u,v)=1}(1-f_v)}{\operatorname{deg}(u)}\)。改一下形式，变为 \(\operatorname{deg}(u)·f_u+\sum\limits_{w(u,v)=1}f_v-\sum\limits_{w(u,v)=0}f_v=\sum\limits_{w(u,v)=1}1\)。
注意到了点 \(n\) 后就不回去了。所以 \(f_n=0\)。我因为这个调了好久。
点击查看代码
int n,m,d[110],cnt=0;
vector<int>v[110][110];
double f[110][110],ans=0;
int main(){
n=read();m=read();
for(int i=1,x,y,z;i<=m;i++){
x=read();y=read();z=read();
v[x][y].push_back(z);d[x]++;
if(x^y)v[y][x].push_back(z),d[y]++;
}
for(int p=0;p<=30;p++){
for(int i=1;i<=n;i++)
for(int j=1;j<=n+1;j++)f[i][j]=0;
for(int i=1;i<n;i++){
f[i][i]=d[i];
for(int j=1;j<=n;j++)
for(int k=0;k<v[i][j].size();k++)
if(!((v[i][j][k]>>p)&1))f[i][j]-=1;
else f[i][n+1]+=1,f[i][j]+=1;
}
f[n][n]=-1;
for(int i=1,x=i;i<=n;i++,x=i){
for(int j=i+1;j<=n;j++)
if(fabs(f[j][i])>fabs(f[x][i]))x=j;
for(int j=1;j<=n+1;j++)swap(f[i][j],f[x][j]);
for(int j=1;j<=n;j++)
if(j!=i){
double y=f[j][i]/f[i][i];
for(int k=i+1;k<=n+1;k++)f[j][k]-=y*f[i][k];
}
}
ans+=f[1][n+1]/f[1][1]*double(1<<p);
}
printf("%.3f\n",ans);
return 0;
}
