VUA 11828 Back to Kernighan-Ritchie（状态机，高斯消元法解线性方程组）

题目链接

初次做线性方程组的题。

可用高斯-约当消元法将原增广矩阵化为对角阵。需要注意无穷解、不定解等特殊情况。。

用一组数据解释一些特殊情况：

5

1 2

2 3

3 2

4 5

5 4

0 0

5

1 2 3 4 5

0

化简后的对角阵为：

1.000000  0.000000  0.000000  0.000000  0.000000  1.000000
0.000000  1.000000  -1.000000  0.000000  0.000000  1.000000
0.000000  0.000000  0.000000  0.000000  0.000000  1.000000
0.000000  0.000000  0.000000  1.000000  -1.000000  0.000000
0.000000  0.000000  0.000000  0.000000  0.000000  0.000000

明显答案为： f(1) = 1 , f(2) = f(3 ) = inf  , f(4)  = f(5 ) = 0;

判断方法为：

1. 无穷解：" m[ i ][ i ] = 0  &&  m[ i ][ n ] != 0 "  || "m[ i ][ j ] != 0  &&  f( j ) = inf  "  

2. 不定解："m[ i ][ i ] = m[ i ][ n ] = 0"  （ 注意这题应输出0.000 ，看上面数据中的5）

3.其他情况有唯一解： f( i ) = m[ i ][ n ] / m[ i ][ i ];

参考代码：

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;
#define MEM(a) memset(a,0,sizeof(a))

const int maxn = 110;
const double eps = 1e-8;
typedef double Mat[maxn][maxn];

int head[maxn] , next[maxn*maxn] , pre[maxn*maxn] ,en;
int inf[maxn] , d[maxn] , N;
Mat m;

void init(){
    MEM(head) , MEM(m) ,MEM(d) ,MEM(inf);
    en = 1;
}
void add_fa(int u,int v){
    next[en] = head[u] , pre[en] = v , head[u] = en++;
}
void Debug(){
    printf("\nMat:\n");
    for(int i=0;i<N;i++){
        for(int j=0;j<=N;j++){
            printf("%f  ",m[i][j]);
        }
        printf("\n");
    }
}

void gause_jordan(Mat &m ,int n){
    int i , j , k ,r;
    double K;
    for(i=0;i<n;i++){
        r = i;
        for(j=i+1 ;j<n;j++){
            if(abs(m[j][i]) > abs(m[r][i])) r = j;
        }
        if(abs(m[r][i]) < eps) continue;
        if(r!=i) for(int j=0;j<=n;j++) swap(m[i][j] , m[r][j]);
        for(k=0;k<n;k++) {
            if(k==i) continue;
            K = m[k][i] / m[i][i];
            for(j=n;j>=i;j--)
                m[k][j] -= K*m[i][j];
        }
    }
}

int main()
{
    int cas = 1;
    while(scanf("%d",&N)!=EOF && N){
        init();
        int u , v;
        do{
            scanf("%d%d",&u,&v);
            add_fa(v-1,u-1);
            d[u-1]++;
        }while(u+v);

        m[0][N] = 1;
        for(int i=0;i<N;i++){
            m[i][i] = 1.0;
            for(int p=head[i] ; p ; p=next[p]){
                int f = pre[p];
                m[i][f] -= 1.0 / d[f];
            }
        }

        gause_jordan(m , N);
        for(int i=N-1 ;i>=0;i--){
            if(abs(m[i][i] < eps && abs(m[i][N]) > eps)) {
                inf[i] = 1; continue;
            }
            for(int j=i+1;j<N;j++) {
                if(abs(m[i][j]) > eps && inf[j]) {
                    inf[i] = 1; break;
                }
            }
        }
        //Debug();
        int Q,q;
        printf("Case #%d:\n",cas++);
        scanf("%d",&Q);
        while(Q--){
            scanf("%d",&q); q--;
            if(inf[q]) printf("infinity\n");
            else printf("%.3f\n",abs(m[q][q]) < eps ? 0.0 : m[q][N] / m[q][q]);
        }
    }
    return 0;
}