#include "Solver.h"

#include "Gauss.h"

Vect Solver::Image(Matrix&& a_Matrix) const {
	a_Matrix.Transpose();
	Gauss::GaussJordan(a_Matrix, false, false);
	a_Matrix.Transpose();
	return {std::move(a_Matrix)};
}

// https://en.wikipedia.org/wiki/Kernel_(linear_algebra)#Computation_by_Gaussian_elimination
Vect Solver::Kernel(Matrix&& a_Matrix) const {
	std::size_t matrixRawCount = a_Matrix.GetRawCount();
	std::size_t matrixColumnCount = a_Matrix.GetColumnCount();

	a_Matrix.Transpose();
	a_Matrix.Augment(Matrix::Identity(a_Matrix.GetRawCount()));
	Gauss::GaussJordan(a_Matrix, false, false);
	a_Matrix.Transpose();

	// nombre de colonnes non nulles
	std::size_t origine_colonne = Vect(a_Matrix.SubMatrix(0, 0, matrixRawCount, matrixColumnCount)).GetCardinal();

	return {a_Matrix.SubMatrix(
		matrixRawCount, origine_colonne, a_Matrix.GetRawCount() - matrixRawCount, a_Matrix.GetColumnCount() - origine_colonne)};
}

VectAffine Solver::RectangularSystem(Matrix&& a_MatrixA, const Matrix& a_VectorB) const {
	Matrix mat = a_MatrixA;
	mat.Augment(a_VectorB);
	Gauss::GaussJordan(mat, true, true);

	Solver solver;

	Vect noyau = solver.Kernel(std::move(a_MatrixA));
	Matrix origin = mat.SubMatrix(0, mat.GetColumnCount() - 1, mat.GetRawCount(), 1);

	// on rajoute des 0 si il faut

	Matrix fullOrigin {mat.GetColumnCount() - 1, 1};
	for (std::size_t i = 0; i < mat.GetRawCount(); i++) {
		fullOrigin.at(i, 0) = origin.at(i, 0);
	}

	for (std::size_t i = mat.GetRawCount(); i < mat.GetColumnCount() - 1; i++) {
		fullOrigin.at(i, 0) = 0;
	}

	return {noyau, fullOrigin};
}

std::size_t Solver::Rank(Matrix&& a_Matrix) const {
	return Image(std::move(a_Matrix)).GetCardinal();
}