uva 10828 Back to Kernighan-Ritchie(高斯消元)

高斯消元 和 高斯-约当消元法（gauss-jordan matrix elimination）的区别在于：

前者使用阶梯型矩阵，采用代入法求解。

后者采用了简化阶梯形矩阵，直接得到解，并且可以方便地处理无穷解，无解的情况。

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <vector>
#include <queue>
#include <stack>
#include <cassert>
#include <algorithm>
#include <cmath>
#include <limits>
#include <set>
#include <map>

using namespace std;

#define MIN(a, b) a < b ? a : b
#define MAX(a, b) a > b ? a : b
#define F(i, n) for (int i=0;i<(n);++i)
#define REP(i, s, t) for(int i=s;i<=t;++i)
#define IREP(i, s, t) for(int i=s;i>=t;--i)
#define REPOK(i, s, t, o) for(int i=s;i<=t && o;++i)
#define MEM0(addr, size) memset(addr, 0, size)
#define LBIT(x) x&-x

#define PI 3.1415926535897932384626433832795
#define HALF_PI 1.5707963267948966192313216916398
#define eps 1e-8

#define MAXN 100 + 10
#define MAXM 100
#define MOD 20071027

typedef long long LL;

const double maxdouble = numeric_limits<double>::max();
const int INF = 0x7FFFFFFF;

typedef double Matrix[MAXN + 1][MAXN + 1];

void gauss_jorgan(Matrix &A, int n) {
int i, j, k, r;
F(i, n) {
r = i;
REP(j, i+1, n-1)
if (fabs(A[j][i]) > fabs(A[r][i]))
r = j;
if (fabs(A[r][i]) < eps) // 此行对角元素为0，跳过
continue;
if (r != i)
REP(j, 0, n)
swap(A[i][j], A[r][j]);
// 用A[i][i]同一列的其他元素
// 最终得到化简阶梯式
REP(k, 0, n-1)
if (k != i)
IREP(j, n, i)
A[k][j] -= A[k][i]/A[i][i]*A[i][j];
}
}

Matrix A;
int n;
int d[MAXN + 1];// 记录出度
vector<int> _prev[MAXN + 1];// 前继集合
int inf[MAXN + 1];

void print(int n) {
REP(i, 0, n-1) {
REP(j, 0, n)
//printf("%.3lf ", A[i][j]);
cout << A[i][j] << " ";
printf("\n");
}
}
int main()
{
freopen("input.in", "r", stdin);
int n, u, v, q, tmp;
int ncase = 0;
std::cout.setf( std::ios::fixed, std:: ios::floatfield );
std::cout.precision(3);
while(scanf("%d", &n) && n) {
REP(i, 0, n-1)
_prev[i].clear();
memset(inf, 0, sizeof(inf));
memset(d, 0, sizeof(d));
while(scanf("%d%d", &u, &v) == 2 && u && v) {
v--;
u--;
d[u]++;
_prev[v].push_back(u);
}

// 构造方程
F(i, n) {
REP(j, 0, n)
A[i][j] = 0;
int _size = (double)_prev[i].size();
F(j, _size)
A[i][_prev[i][j]] -= 1.0/d[_prev[i][j]];
A[i][i] = 1;
}
A[0]
= 1;

gauss_jorgan(A, n);
/*
F(i, n) {
if (fabs(A[i][i]) < eps && fabs(A[i]
) > eps)
inf[i] = 1;
}
F(i, n)
REP(j, i+1, n-1)
if (fabs(A[i][j]) > eps && inf[j]) {
inf[i] = 1;
break;
}
*/
IREP(i, n-1, 0) {
if (fabs(A[i][i]) < eps && fabs(A[i]
) > eps)
inf[i] = 1;
REP(j, i+1, n-1)
if (fabs(A[i][j]) > eps && inf[j])
inf[i] = 1;
}

scanf("%d", &q);
printf("Case #%d:\n", ++ncase);
F(i, q) {
scanf("%d", &tmp);
--tmp;
if (inf[tmp])
printf("infinity\n");
else {
if (fabs(A[tmp]
) < eps)
cout << 0.0 << endl;
else
cout << A[tmp]
/A[tmp][tmp] << endl;
}

}
}
return 0;
}