Template -「高斯消元」

#include <bits/stdc++.h>
using namespace std;

typedef double Doub;
typedef long long LL;
typedef pair <int, int> PII;
// template <typename Tem> inline Tem Abs (Tem x) { return x < 0 ? -x : x; }
template <typename Tem> inline Tem Max (Tem x, Tem y) { return x > y ? x : y; }
template <typename Tem> inline Tem Min (Tem x, Tem y) { return x < y ? x : y; }

template <typename Tem>
inline void Read (Tem &x) {
    bool Flag = false; x = 0;
    char s = getchar ();
    while (s < '0' || s > '9') {
        if (s == '-')
           Flag = true;
        s = getchar ();
    }
    while ('0' <= s && s <= '9')
        x = (x << 3) + (x << 1) + (s ^ 48), s = getchar ();
    if (Flag)
        x = -x;
}

template <typename Tem>
inline void Write (Tem x) {
    if (x < 0)
        putchar ('-'), x = -x;
    if (x > 9)
       Write (x / 10);
    putchar (x % 10 + '0');
}

template <typename Tem>
inline void Print (Tem x, char s) { Write (x), putchar (s); }

const int Maxn = 4e2 + 5;

struct Eliminate {
	bool Free[Maxn];
	int n, m, Rank, Opt;
	double a[Maxn][Maxn], x[Maxn], Eps;	
	Eliminate () { Eps = 1e-16, memset (a, 0, sizeof a); }
	Eliminate (int N, int M) { n = N, m = M; }

    double Abs (double x) { return x < Eps ? -x : x; }
	
	void Calc () {
		int r = 1, c = 1, Pos;
		for (; r <= n && c <= m; r++, c++) {
			Pos = r;
			for (int i = r + 1; i <= n; i++)
				if (Abs (a[i][c]) > Abs (a[Pos][c]))
					Pos = i;
			if (Abs (a[Pos][c]) < Eps) {
				r--;
				continue;
			}
			if (Pos != r)
				for (int i = c; i <= m; i++)
					swap (a[r][i], a[Pos][i]);
			double x;
			for (int i = 1; i <= n; i++)
				if (i != r && Abs (a[i][c]) > Eps) {
					x = a[i][c] / a[r][c];
					for (int j = m; j >= c; j--)
						a[i][j] -= x * a[r][j];
				}
		}
		Rank = r;
	}
	
	void Check () {
        for (int i = 1; i <= m; i++) 
            Free[i] = 1;
        for (int i = Rank - 1, Pos, Cnt; i >= 0; i--) {
            Pos = 0, Cnt = 0;
            for (int j = 0; j < m; j++) 
                if ((Abs (a[i][j]) > Eps) && Free[j]) 
                    ++Cnt, Pos = j;
            if (Cnt == 1) 
                Free[Pos] = false, x[Pos] = a[i][m] / a[i][Pos];
        }
	}
} q;

int main () {
	int n, m; Read (n), Read (m);
	q.n = n, q.m = m + 1;
	for (int i = 1; i <= n; i++) 
		for (int j = 1; j <= m + 1; j++)	
			scanf ("%lf", &q.a[i][j]);
	q.Calc (), q.Check ();
	if (q.Opt == -1)
		printf ("No Solution\n");
	else
		for (int i = 1; i <= n; i++)
            if (q.Free [i])
                printf ("X[%d] not determined\n", i);
            else
			    printf ("X[%d] = %.4f\n", i, q.x[i]);
	return 0;
}
posted @ 2021-10-18 15:42  STrAduts  阅读(36)  评论(0编辑  收藏  举报