Gauss 消元

Gauss 消元・$\mathbb{R}$

std::optional<vec2d<f64>> gauss(vec2d<f64> v, f64 eps = 1E-6) {
	for (int i = 0; i < v.x; i++) {
		int mi = i;
		for (int j = i + 1; j < v.x; j++)
			if (std::abs(v[j][i]) > std::abs(v[mi][i]))
				mi = j;
		for (int j = 0; j < v.y; j++)
			std::swap(v[i][j], v[mi][j]);
		if (std::abs(v[i][i]) < eps)
			return std::nullopt;
		auto ivi = 1 / v[i][i];
		for (int k = v.y - 1; k >= i; k--)
			v[i][k] *= ivi;
		for (int j = 0; j < v.x; j++) {
			if (j != i) {
				auto tmp = v[j][i];
				for (int k = 0; k < v.y; k++)
					v[j][k] -= v[i][k] * tmp;
			}
		}
	}
	return v;
}

Gauss 消元・${\mathbb{Z}}_p$

template <class Z>
std::optional<vec2d<Z>> gauss(vec2d<Z> v) {
	for (int i = 0; i < v.x; i++) {
		int mi = i;
		for (int j = i + 1; j < v.x; j++)
			if (v[j][i].v != 0)
				mi = j;
		for (int j = 0; j < v.y; j++)
			std::swap(v[i][j], v[mi][j]);
		if (v[i][i].v == 0)
			return std::nullopt;
		auto ivi = v[i][i].inv();
		for (int k = v.y - 1; k >= i; k--)
			v[i][k] *= ivi;
		for (int j = 0; j < v.x; j++) {
			if (j != i) {
				auto tmp = v[j][i];
				for (int k = 0; k < v.y; k++)
					v[j][k] -= v[i][k] * tmp;
			}
		}
	}
	return v;
}