洛谷 P3389 【模板】高斯消元法 × 洛谷 P2455 [SDOI2006]线性方程组

一、传送门

https://www.luogu.com.cn/problem/P3389
https://www.luogu.com.cn/problem/P2455

二、代码

日……两道题合起来搞了我十几个小时，最后把我气得从别人的 AC 代码开始再一句一句改成适配自己模板的同义的语句。
有些细节还是不太懂，暂时先把模板放上来。以后有时间再出详解。
————————————————————————————————————————————————————
高斯消元函数 gauss_elim() 的返回值：
0，唯一解；1，格式错（矩阵不是n行×(n+1)列）；2，无数解；3，无解。
————————————————————————————————————————————————————
P3389：

#include<cstdio>
#include<cmath>
#include<algorithm>
#include<cstdint>
#include<vector>
#pragma warning(disable:4996)
using namespace std;
template<class _Ty = double> class matrix {
private:
	vector<vector<_Ty>> m; size_t row, col, _size;
public:
	matrix() { row = col = _size = 0; }
	matrix(const size_t& r, const size_t& c) { init(r, c); }
	virtual ~matrix() { del(); }
	void clear() { fill(m[0], m[0] + _size, 0); }
	void del() { m.resize(0); }
	void init(const size_t& r, const size_t& c) {
		row = r; col = c; _size = row * col; m.resize(r); for (size_t i = 0; i != r; ++i)m[i].resize(c);
	}
	void resize(const size_t& r, const size_t& c) const { del(); init(r, c); }
	void resize_directly(const size_t& r, const size_t& c) { row = r; col = c; _size = row * col; }
	size_t _row() const { return row; }
	size_t _col() const { return col; }
	size_t size() const { return _size; }
	void row_mul(const size_t& row, const _Ty& val) {
		for (size_t i = 0; i != col; ++i)m[row][i] *= val;
	}
	void row_mul_add_to(const size_t& rs, const _Ty& val, const size_t& rd) {
		for (size_t i = 0; i != col; ++i)m[rd][i] += m[rs][i] * val;
	}
	void interchange(const size_t& r1, const size_t& r2) { swap(m[r1], m[r2]); }
	vector<_Ty>& operator[](const size_t& r) { return m.at(r); }
	const vector<_Ty>& operator[](const size_t& r) const { return m.at(r); }
};
template<class _Ty = double> inline unsigned gauss_elim(matrix<_Ty>& a, _Ty* const x) {
	_Ty y, b; size_t n = a._row(), p, c; bool inf = false;
	for (size_t i = 0; i < n; ++i) {
		p = i; b = abs(a[p][i]);
		for (size_t j = i + 1; j < n; ++j)if (abs(a[j][i]) > b) { p = j; }
		if (p != i)a.interchange(p, i);
		if (a[i][i] == 0)continue;
		b = 1 / a[i][i];
		for (size_t j = 0; j < i; ++j)a.row_mul_add_to(i, -a[j][i] * b, j);
		for (size_t j = i + 1; j < n; ++j)a.row_mul_add_to(i, -a[j][i] * b, j);
	}
	for (size_t i = 0; i < n; ++i) {
		c = 0;
		for (size_t j = 0; j <= n; ++j) { if (a[i][j] == 0)++c; }
		if (c == n && a[i][n]) { return 3; }
		if (c == a._col()) { inf = true; }
	}
	if (inf)return 2;
	for (size_t i = n - 1; i != SIZE_MAX; --i) {
		y = a[i][n]; for (size_t j = n - 1; j > i; --j) { y -= a[i][j] * x[j]; }
		x[i] = y / a[i][i];
	}
	return 0;
}
matrix<double> m; size_t n; double x[100];
int main() {
	scanf("%llu", &n); m.init(n, n + 1);
	for (size_t i = 0; i < n; ++i)
		for (size_t j = 0; j <= n; ++j)scanf("%lf", &m[i][j]);
	if (!gauss_elim(m, x))for (size_t i = 0; i < n; ++i)printf("%.2lf\n", x[i]);
	else puts("No Solution");
	return 0;
}

P2455：

#include<cstdio>
#include<cmath>
#include<algorithm>
#include<cstdint>
#include<vector>
#pragma warning(disable:4996)
using namespace std;
template<class _Ty = double> class matrix {
private:
	vector<vector<_Ty>> m; size_t row, col, _size;
public:
	matrix() { row = col = _size = 0; }
	matrix(const size_t& r, const size_t& c) { init(r, c); }
	virtual ~matrix() { del(); }
	void clear() { fill(m[0], m[0] + _size, 0); }
	void del() { m.resize(0); }
	void init(const size_t& r, const size_t& c) {
		row = r; col = c; _size = row * col; m.resize(r); for (size_t i = 0; i != r; ++i)m[i].resize(c);
	}
	void resize(const size_t& r, const size_t& c) const { del(); init(r, c); }
	void resize_directly(const size_t& r, const size_t& c) { row = r; col = c; _size = row * col; }
	size_t _row() const { return row; }
	size_t _col() const { return col; }
	size_t size() const { return _size; }
	void row_mul(const size_t& row, const _Ty& val) { for (size_t i = 0; i != col; ++i)m[row][i] *= val; }
	void row_mul_add_to(const size_t& rs, const _Ty& val, const size_t& rd) {
		for (size_t i = 0; i != col; ++i)m[rd][i] += m[rs][i] * val;
	}
	void interchange(const size_t& r1, const size_t& r2) { swap(m[r1], m[r2]); }
	vector<_Ty>& operator[](const size_t& r) { return m.at(r); }
	const vector<_Ty>& operator[](const size_t& r) const { return m.at(r); }
};
template<class _Ty = double> inline unsigned gauss_elim(matrix<_Ty>& a, _Ty* const x) {
	_Ty y, b; size_t n = a._row(), p, c; bool inf = false;
    for (size_t i = 0; i < n; ++i) {
		p = i; b = abs(a[p][i]);
		for (size_t j = i + 1; j < n; ++j)if (abs(a[j][i]) > b) { p = j; }
		if (p != i)a.interchange(p, i);
		if (a[i][i] == 0)continue;
		b = 1 / a[i][i];
        for (size_t j = 0; j < i; ++j)a.row_mul_add_to(i, -a[j][i] * b, j);
        for (size_t j = i + 1; j < n; ++j)a.row_mul_add_to(i, -a[j][i] * b, j);
    }
    for (size_t i = 0; i < n; ++i) {
        c = 0;
        for (size_t j = 0; j <= n; ++j) { if (a[i][j] == 0)++c; }
		if (c == n && a[i][n]) { return 3; }
		if (c == a._col()) { inf = true; }
    }
    if (inf)return 2;
	for (size_t i = n - 1; i != SIZE_MAX; --i) {
		y = a[i][n]; for (size_t j = n - 1; j > i; --j) { y -= a[i][j] * x[j]; }
		x[i] = y / a[i][i];
	}
	return 0;
}
matrix<double> m; size_t n; double x[50];
int main() {
	scanf("%llu", &n); m.init(n, n + 1);
	for (size_t i = 0; i < n; ++i)
		for (size_t j = 0; j <= n; ++j)scanf("%lf", &m[i][j]);
	switch (gauss_elim(m, x)) {
	case 3:puts("-1"); return 0;
	case 2:puts("0"); return 0;
	default:for (size_t i = 0; i < n; ++i)printf("x%llu=%.2lf\n", i + 1, x[i]); return 0;
	}
}