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

最新推荐文章于 2019-03-20 17:04:53 发布

原创最新推荐文章于 2019-03-20 17:04:53 发布 · 535 阅读

0 ·

CC 4.0 BY-SA版权

文章标签：

#高斯消元

ACM.数学专栏收录该内容

17 篇文章

订阅专栏

本文详细介绍了高斯-约当消元法的基本原理及其与高斯消元法的区别，通过具体代码实现展示了如何求解线性方程组，并特别关注于简化阶梯形矩阵的构建过程。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

高斯消元和高斯-约当消元法（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][n] = 1;

        gauss_jorgan(A, n);
/*
        F(i, n) {
            if (fabs(A[i][i]) < eps && fabs(A[i][n]) > 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][n]) > 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][n]) < eps)
                    cout << 0.0 << endl;
                else
                    cout << A[tmp][n]/A[tmp][tmp] << endl;
            }

        }
    }
    return 0;
}