高斯消元 复习

高斯消元法是用来解线性方程组的一种算法

板子

for(reg i = 1;i < n;i++)
	{
		int r = i;
		while(!a[r][i]) r++;
		swap(a[r],a[i]); 
		double pass = a[i][i];
		for(reg k = 1;k <= n;k++) a[i][k] /= pass;
		for(reg j = 1;j < n;j++)
		{
			if(j == i) continue;
			double div = a[j][i];
			for(reg k = 1;k <= n;k++)
				a[j][k] -= a[i][k] * div;
		}
	}

用途

这里主要说下用途

做很多关于期望的题时

\[dp[i]=f(dp[j]) \]

\[dp[j]=f(dp[i]) \]

出现了以上形式的\(dp\)方程

明显直接\(dp\)是不行的

但如果 我们能求出

\(n\)个关于\(dp[i]=f(dp[j])(i,j\in[1,n],s.t.i,j\in Z)\)

就可以高斯消元来求这些未知量

例题

P3232 [HNOI2013]游走

#include <stack>
#include <cstdio>
#include <vector>
#include <string>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
#define reg register int
#define isdigit(x) ('0' <= x&&x <= '9')
template<typename T>
inline T Read(T Type)
{
	T x = 0,f = 1;
	char a = getchar();
	while(!isdigit(a)) {if(a == '-') f = -1;a = getchar();}
	while(isdigit(a)) {x = (x << 1) + (x << 3) + (a ^ '0');a = getchar();}
	return x * f;
}
const int MAXN = 510,MAXM = MAXN * MAXN;
struct node
{
	int v,next;
}edge[MAXM << 1];
double ans,a[MAXN][MAXN],p[MAXN],edge_p[MAXM];
int cnt,n,m,head[MAXN],deg[MAXN];
inline void addedge(int u,int v)
{
	deg[u]++;
	edge[++cnt].v = v;
	edge[cnt].next = head[u];
	head[u] = cnt;
}
inline void gauss()
{
	a[1][n] = -1;
	for(reg i = 1;i < n;i++)
	{
		for(reg e = head[i];e;e = edge[e].next)
		{
			int v = edge[e].v;
			if(v == n) continue;
			a[i][v] += 1.0 / deg[v];
		}
		a[i][i] -= 1;
	}
	for(reg i = 1;i < n;i++)
	{
		int r = i;
		while(!a[r][i]) r++;
		swap(a[r],a[i]);
		double pass = a[i][i];
		for(reg k = 1;k <= n;k++) a[i][k] /= pass;
		for(reg j = 1;j < n;j++)
		{
			if(j == i) continue;
			double div = a[j][i];
			for(reg k = 1;k <= n;k++)
				a[j][k] -= a[i][k] * div;
		}
	}
	for(reg i = 1;i < n;i++) p[i] = a[i][n] / a[i][i];
	cnt = 0;
	for(reg i = 1;i <= n;i++)
		for(reg e = head[i];e;e = edge[e].next)
		{
			int v = edge[e].v;
			if(v < i) continue;
			edge_p[++cnt] = (i == n?0:p[i] / deg[i]) + (v == n?0:p[v] / deg[v]);
		}
}
int main()
{
	n = Read(1),m = Read(1);
	for(reg i = 1;i <= m;i++)
	{
		int u = Read(1),v = Read(1);
		addedge(u,v),addedge(v,u);
	}
	gauss();
	sort(edge_p + 1,edge_p + 1 + cnt);
	for(reg i = 1;i <= m;i++) ans += edge_p[i] * (m - i + 1);
	printf("%.3f",ans);
	return 0;
}

P3211 [HNOI2011]XOR和路径

#include <stack>
#include <cstdio>
#include <vector>
#include <string>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
#define reg register int
#define isdigit(x) ('0' <= x&&x <= '9')
template<typename T>
inline T Read(T Type)
{
	T x = 0,f = 1;
	char a = getchar();
	while(!isdigit(a)) {if(a == '-') f = -1;a = getchar();}
	while(isdigit(a)) {x = (x << 1) + (x << 3) + (a ^ '0');a = getchar();}
	return x * f;
}
const int MAXN = 110,MAXM = 10010;
struct node
{
	int v,w,next;
}edge[MAXM << 1];
double ans,a[MAXN][MAXN];
int n,m,maxt,cnt,head[MAXN],deg[MAXN];
inline void addedge(int u,int v,int w)
{
	edge[++cnt].v = v;
	edge[cnt].w = w;
	edge[cnt].next = head[u];
	head[u] = cnt;
}
inline int getlen(int x) {int res = 0;while(x) x >>= 1,res++;return res;}
inline int getx(int x,int t) {return (x >> (t - 1) & 1);}
inline double fabs(double x) {if(x < 0) return -x;return x;}
inline void gauss(int t)
{
	for(reg i = 1;i <= n;i++)
		for(reg j = 1;j <= n + 1;j++)
			a[i][j] = 0;
	for(reg i = 1;i < n;i++)
	{
		for(reg j = head[i];j;j = edge[j].next)
		{
			int v = edge[j].v,w = getx(edge[j].w,t);
			if(w) a[i][v] -= 1.0,a[i][n + 1] -= 1.0;
			else a[i][v] += 1;
		}
		a[i][i] -= deg[i];
	}
	a[n][n] = 1;
	for(reg i = 1;i <= n;i++)
	{
		int r = i;
		while(!a[r][i]) r++;
		swap(a[r],a[i]);
		double pass = a[i][i];
		for(reg j = 1;j <= n + 1;j++) a[i][j] /= pass;
		for(reg j = 1;j <= n;j++)
		{
			if(j == i) continue;
			double pass = a[j][i];
			for(reg k = 1;k <= n + 1;k++)
				a[j][k] -= a[i][k] * pass;
		}
	}
	ans += a[1][n + 1] / a[1][1] * (1 << (t - 1));
}
int main()
{
	n = Read(1),m = Read(1);
	for(reg i = 1;i <= m;i++)
	{
		int u = Read(1),v = Read(1),w = Read(1);
		addedge(u,v,w),deg[u]++;
		maxt = max(maxt,getlen(w));
		if(u != v) addedge(v,u,w),deg[v]++;
	}
	for(reg i = 1;i <= maxt;i++) gauss(i);
	printf("%.3f",ans);
	return 0;
}

P4035 [JSOI2008]球形空间产生器

#include <stack>
#include <cstdio>
#include <vector>
#include <string>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
#define reg register int
#define isdigit(x) ('0' <= x&&x <= '9')
template<typename T>
inline T Read(T Type)
{
	T x = 0,f = 1;
	char a = getchar();
	while(!isdigit(a)) {if(a == '-') f = -1;a = getchar();}
	while(isdigit(a)) {x = (x << 1) + (x << 3) + (a ^ '0');a = getchar();}
	return x * f;
}
const int MAXN = 15;
double pow[MAXN],a[MAXN][MAXN];
inline double fabs(double x) {if(x < 0) return -x;return x;}
struct Matrix
{
	int n,m;
	double c[MAXN][MAXN];
}T;
int main()
{
	int n = Read(1);
	for(reg i = 1;i <= n + 1;i++)
	{
		for(reg j = 1;j <= n;j++)
		{
			scanf("%lf",&a[i][j]);
			pow[i] += a[i][j] * a[i][j];
			if(i != 1) T.c[i - 1][j] = 2 * (a[i - 1][j] - a[i][j]);
		}
		if(i != 1) T.c[i - 1][n + 1] = pow[i - 1] - pow[i];
	}
	for(reg i = 1;i <= n;i++)
	{
		int r = i;
		for(reg j = i + 1;j <= n;j++)
			if(fabs(T.c[r][i]) > fabs(T.c[j][i])) r = j;
		swap(T.c[r],T.c[i]);
		double pass = T.c[i][i];
		for(reg j = 1;j <= n + 1;j++) T.c[i][j] /= pass;
		for(reg j = 1;j <= n;j++)
		{
			if(i == j) continue;
			double pass = T.c[j][i];
			for(reg k = 1;k <= n + 1;k++) T.c[j][k] -= T.c[i][k] * pass;
		}
	}
	for(reg i = 1;i <= n;i++)
		printf("%.3f ",T.c[i][n + 1]);
	return 0;
}