高斯消元

高斯消元是用来解 \(n\) 元 \(1\) 次方程组的办法。

P3389 【模板】高斯消元法

对于人而言，解方程组可能有一些技巧，但是对于计算机，我们需要通式通法。
对于一次方程组，一个比较通用的思路是加减消元。对于每个式子，我们将其化为只剩下一个未知数有系数，其余未知数的系数都为 \(0\) 的式子。

考虑如何去做呢？我们以样例为例：

\[\begin{bmatrix} 1 & 3 & 4& 5\\ 1 &3 & 7 &3 \\ 9& 3 & 2 &2 \end{bmatrix} \]

我们将所有方程装在一个 \(n\times (n+1)\) 的矩阵中，最后一列的数表示等号右边的常数项。

高斯-约旦消元法

我们这里只介绍比较常用的高斯-约旦消元法。这种消元的目标是将与上面类似的矩阵化为以下形式：

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

发现化为以下形式后，每一行都只有一个未知数带有系数。这样对于每一行，我们就可以直接通过常数项除以这个系数得到这一行对应的这个未知数。

考虑如何得到这个形式。我们直接暴力一点，以当前行的系数为基底，去消去其余所有行的对应的系数。

注意到，这样的操作方式可以保证已经操作的行的前缀系数全部为 0，所以直接操作并不会改变这个性质，因为当前行对应的系数也为 0。

唯一一个不是很直观的问题是，我们要找到当前对应的未知数系数最大的那一行，将那一行换过来。可以证明这样误差最小，还可以避免当前行对应系数恰好是 0 的问题。
可以自己列一下式子，比较直观这里就不赘述了。

code

注意到无解或者说有无限解的特征是对于当前这个未知数，所有式子关于其的系数都为 0。这个还是比较好理解的，因为这样所有方程都跟这个数没有关系了。可能导致最后剩下一个方程是 \(0\times x_i=b\)，而这个 \(b\) 为 \(0\) 或其他实数时分别对应无限解和无解的情况。
时间复杂度比较显然，\(O(n^3)\)。所以说这个相对暴力。

#include<bits/stdc++.h>
using namespace std;
#define ld double
const int N=105;
int n;double a[N][N];
signed main(){
	ios::sync_with_stdio(false),cin.tie(0),cout.tie(0);
	cin>>n;for(int i=1;i<=n;i++)for(int j=1;j<=n+1;j++) cin>>a[i][j];
	for(int i=1;i<=n;i++){
		int m=i;
		for(int j=i+1;j<=n;j++) if(fabs(a[m][i])<fabs(a[j][i]))m=j;
		for(int j=1;j<=n+1;j++) swap(a[i][j],a[m][j]);
		if(a[i][i]==0) {cout<<"No Solution";return 0;}
		for(int j=1;j<=n;j++){
			if(i!=j){
				double d=a[j][i]/a[i][i];
				for(int k=i+1;k<=n+1;k++) a[j][k]-=a[i][k]*d;
			}
		}
	}
	for(int i=1;i<=n;i++){
		cout<<fixed<<setprecision(2)<<a[i][n+1]/a[i][i]<<'\n';
	}
	return 0;
}

P3232 [HNOI2013] 游走

先转化为求每个点的期望经过次数 \(f_u\)。显然有

\[f_u=\frac{\sum_{v\ne n}f_v}{tot_u}+[u=1],u\le n-1 \]

其中 \(v\) 为与 \(u\) 连边的点，\(tot_u\) 表示 \(u\) 的出边条数。这个东西明显不能直接算，有后效性。
由于 \(n\) 比较小，因此考虑高斯消元，将 \(f_u\) 挪到右边，常数项挪到左边，变成类似模板题的标准方程后直接消。
然后考虑边的期望经过次数 \(g_{u,v}\)，其中 \(u,v\) 是这条边的两个端点。
显然有

\[g_{u,v}=\frac{f_u}{tot_u}\times [u\ne n]+\frac{f_v}{tot_v}\times [v\ne n] \]

然后贪心直接排序计算贡献即可。

code

#include<bits/stdc++.h>
using namespace std;
#define ld long double
const int N=1e3+7,M=2e5+7;const ld eps=-1e12;
int n,m,tot[N],mp[N][N],idcnt=0;
ld a[N][N],f[N],g[M];
vector <int> q[N];
bool eq(ld x,ld y){return fabs(x-y)<eps;}
void gaosi(){
	for(int i=1;i<n;i++){
		int m=i;
		for(int j=i+1;j<=n;j++) if(fabs(a[j][i])-fabs(a[m][i])>0.0) m=j;
		for(int j=1;j<=n+1;j++) swap(a[i][j],a[m][j]);
		for(int j=1;j<=n;j++){
			if(j==i) continue;ld d=a[j][i]/a[i][i];
			for(int k=1;k<=n+1;k++) a[j][k]-=a[i][k]*d;
		}
	}
}
signed main(){
	ios::sync_with_stdio(false),cin.tie(0),cout.tie(0);
	cin>>n>>m;for(int i=1,u,v;i<=m;i++) cin>>u>>v,q[u].push_back(v),tot[u]++,q[v].push_back(u),tot[v]++;
	for(int i=1;i<n;i++) {a[i][i]=-1;for(int j=0;j<tot[i];j++) a[i][q[i][j]]=1.0/(ld)tot[q[i][j]];}
	a[1][n+1]=-1;gaosi();
	for(int i=1;i<n;i++) f[i]=a[i][n+1]/a[i][i];
	for(int i=1;i<n;i++){for(int j=0;j<tot[i];j++){int u=i,v=q[u][j],t=mp[u][v];if(!mp[u][v])t=++idcnt,mp[u][v]=mp[v][u]=t;g[t]+=f[u]/(ld)tot[u];}}
	sort(g+1,g+m+1);ld ans=0.0;
	for(int i=1;i<=m;i++) ans+=1.0*(ld)(m-i+1)*g[i];
	cout<<fixed<<setprecision(3)<<ans<<'\n';return 0;
}