高斯消元模板

算法目的：
            高斯消元，一般用于求解线性方程组AX = B（或 模线性方程组AX mod P = B），以四个未知数，四个方程为例，AX=B表示成4x4的矩        阵和4x1的矩阵相乘的形式：  

其中A和B(b0 b1 b2 b3)已知，要求列向量X(x0 x1 x2 x3)的值。

算法核心思想：
            对于n个方程，m个未知数的方程组，消元的具体步骤如下：

1、枚举第i (0 <= i < n) 行，初始化列为col = 0，每次从[i, n)行中找到第col列中元素绝对值最大的行和第i行进行交换(找到最大的行是为了在消元的时候把浮点数的误差降到最小)；

a) 如果第col列的元素全为0，放弃这一列的处理，col+1，i不变，转1)；

b) 否则，对于所有的行j (i < j < n)，如果a[j][col]不为0，则需要进行消元，以期第i行以下的第col列的所有元素都消为0（这一步就是线性代数中所说的初等行变换，具体的步骤就是将第j行的所有元素减去第i行的所有元素乘上一个系数，这个系数即a[j][col] / a[i][col]）。

2、重复步骤1) 直到n个方程枚举完毕或者列col == m。

3、判断解的情况：

a) 如果出现某一行，系数矩阵全为0，增广矩阵不全为0，则无解（即出现[0 0 0 0 0 b]，其中b不等于0的情况）；

b) 如果是严格上三角，则表明有唯一解；

c) 如果增广矩阵有k (k > 0)行全为0，那么表明有k个变量可以任意取值，这几个变量即自由变量；对于这种情况，一般解的范围是给定的，令解的取值有T个，自由变量有V个，那么解的个数就是 TV。

以上转自：http://www.cppblog.com/menjitianya/archive/2014/06/08/207226.html

我理解的自由变元个数就是矩阵的变量数减去矩阵的秩。嗯。就是这样。看了一天的线性代数，还是请教了陈蜜子弄懂了原理，断断续续的敲完了高斯消元。明天开始应用。

#include<limits>
#include<queue>
#include<vector>
#include<list>
#include<map>
#include<set>
#include<deque>
#include<stack>
#include<bitset>
#include<algorithm>
#include<functional>
#include<numeric>
#include<utility>
#include<sstream>
#include<iostream>
#include<iomanip>
#include<cstdio>
#include<cmath>
#include<cstdlib>
#include<cstring>
#include<ctime>

#define LL __int64
#define eps 1e-8
#define pi acos(-1)
#define INF 0x7fffffff
#define delta 0.98 //模拟退火递增变量
using namespace std;
int n,m;
int p[1010][1010];
bool freex[1010];
double x[1010];
void Debug(){
	int i,j;
	for (i=0;i<n;i++){
		for (j=0;j<=m;j++)
			printf("%5d",p[i][j]);
		cout<<endl;
	}
	return;
}
int Gauss(int n,int m){
	int n1=0,m1=0;
	for (;n1<n && m1<m;n1++,m1++){
		int maxn=n1;
		for (int i=n1+1;i<n;i++)
			if(abs(p[i][m1])>abs(p[maxn][m1])) maxn=i;
		if (maxn!=n1)
			for (int i=m1;i<=m;i++)
				swap(p[n1][i],p[maxn][i]);
		if (p[n1][m1]==0){
			n1--;
			continue;
		}
		for (int i=n1+1;i<n;i++){
			if (p[i][m1]){
				LL ta=p[n1][m1];
				LL tb=p[i][m1];
				for (int j=m1;j<m+1;j++)
					p[i][j]=p[i][j]*ta-p[n1][j]*tb;		
			}
		} 
	}
	memset(freex,0,sizeof(freex));
	int freex_t;
	//Debug();
	// 无解 
	for (int i=n1;i<n;i++)
		if (p[i][m1]!=0) return -1;
	// 无穷解 
	if(n1<m){
		for (int i=n1-1;i>=0;i--){
			int freex_num=0;
			for (int j=0;j<m;j++)
				if (p[i][j]!=0 && !freex[j]){
					freex_num++;
					freex_t=j;
				}
			if (freex_num>1) continue; //无法求出变元 
			int t=p[i][m];
			for (int j=0;j<m;j++)
				if (!p[i][j]!=0 && j!=freex_t) t-=p[i][j]*x[j];
			if (t!=0)
				x[freex_t]=t*1.0/p[i][freex_t];
			freex[freex_t]=1;   // 这个值是确定了 
		}
		return m-n1; //自由变元个数 
	}
	// 有唯一解
	for (int i=m-1;i>=0;i--){
		int t=p[i][m];
		for (int j=i+1;j<m;j++)
			if (p[i][j]!=0) t-=p[i][j]*x[j];
		if (t!=0) 
			x[i]=t*1.0/p[i][i];
	}
	return 0;
}

int main(){
	while (~scanf("%d%d",&n,&m)){
		for (int i=0;i<n;i++)
			for (int j=0;j<=m;j++)
				scanf("%d",&p[i][j]);
		memset(x,0,sizeof(x));
		int flag=Gauss(n,m);
		if (flag==-1) cout<<"无解！"<<endl;
		else if (flag==0){
			for (int i=0;i<m;i++)
				printf("x%d:%.0f\n",i+1,x[i]);
		}
		else{
			printf("有无穷解，自由变元个数为:%d\n",flag);
			for (int i=0;i<m;i++)
				if(!freex[i]) printf("x%d是不确定的\n",i+1);
				else
				printf("x%d:%.1f\n",i+1,x[i]);
		} 
	}
	return 0;
}

有些高斯消元解方程要求解出最后解，也就是单位元。以下就是杭电4818的代码。

#include<limits>
#include<queue>
#include<vector>
#include<list>
#include<map>
#include<set>
#include<deque>
#include<stack>
#include<bitset>
#include<algorithm>
#include<functional>
#include<numeric>
#include<utility>
#include<sstream>
#include<iostream>
#include<iomanip>
#include<cstdio>
#include<cmath>
#include<cstdlib>
#include<cstring>
#include<ctime>

#define LL __int64
#define eps 1e-8
#define pi acos(-1)
#define INF 0x7fffffff
#define delta 0.98 //模拟退火递增变量
using namespace std;

double a[210][210];
double x[210];
int du[210];
int g[210][210];
int add[210];
int equ,var;
vector<int>v[210];

int dcmp(double x){
	if (fabs(x)<=eps) return 0;
	if (x>0) return 1;
	else return -1;
} 
int Gauss(){
	int r,col,i,j;
	for (r=0,col=0;r<equ && col<var;r++,col++){
		int maxr=r;
		for (i=r+1;i<equ;i++)
			if (dcmp(fabs(a[i][col])-fabs(a[maxr][col]))>0) maxr=i;
		if (dcmp(a[maxr][col])==0) return 0;//无解，有自由变元 
		if (r!=maxr){
			for (i=col;i<var;i++)
				swap(a[maxr][i],a[r][i]);
			swap(x[r],x[maxr]);
		}
		x[r]/=a[r][col];
		for (i=col+1;i<var;i++)
			a[r][i]=a[r][i]/a[r][col];
		a[r][col] = 1;
		for (i=0;i<equ;i++)
		if (i!=r){
			x[i]=x[i]-x[r]*a[i][col];  
			for (j=col+1;j<var;j++)
				a[i][j]-=a[r][j]*a[i][col];
			a[i][col]=0;
		}	
	}
	return 1;
}

int main(){
	int T,i,j,n,m;
	scanf("%d",&T);
	while (T--){
		scanf("%d%d",&n,&m);  // n表示有n个人，m代表m个关系 
		memset(a,0,sizeof(a));
		memset(g,0,sizeof(g));
		memset(x,0,sizeof(x));
		for (i=0;i<n;i++){
			v[i].clear();
			du[i]=0;
		}
		for (i=1;i<=m;i++){
			int u,v;
			scanf("%d%d",&u,&v);
			g[u][v]=1;
		}
		for (i=0;i<n;i++)
			for (j=0;j<n;j++)
				if (i!=j && g[i][j]){
					du[i]++;
					v[j].push_back(i); 
				}
		equ=var=n;
		for (i=0;i<n;i++){
			a[i][i]=-1;
			int sz=v[i].size();
			for (j=0;j<sz;j++){
				int k=v[i][j];
				if (i==k) continue;
				a[i][k]=1.0/du[k];
			}
		}
		for (i=0;i<n;i++)
			a[n-1][i]=1;
		x[n-1]=1;
		/*for (i=0;i<n;i++){
			for (j=0;j<n;j++)
				printf("%5.2f",a[i][j]);
			cout<<endl;
		}*/
		for (int k=0;k<n-1;k++)
			if (g[n-1][k]==0){
				for (i=0;i<n-1;i++)
					if (g[n-1][i]) a[i][var]=1.0/(du[n-1]+1);
					else a[i][var]=0;
				a[k][var]=1.0/(du[n-1]+1);
				a[n-1][var]=1;
				add[var]=k;
				var++;
			}
		if (!Gauss()){
			printf("INF\n");
			continue;
		}
		double now=x[n-1];
		int ans=-1;
		for (i=n;i<var;i++){
			if (dcmp(x[n-1]/a[n-1][i]-now)>0){
				ans=add[i];
				now=x[n-1]/a[n-1][i];
			}
		}
		printf("%d %d\n",1,ans);
	}
	return 0;
}