#include "Gauss.h"

#include "Matrix.h"

namespace Gauss {

static void SwapLines(Matrix& mat, std::size_t line1, std::size_t line2) {
	for (std::size_t k = 0; k < mat.GetColumnCount(); k++) {
		std::swap(mat.at(line1, k), mat.at(line2, k));
	}
}

static void DivideLine(Matrix& mat, std::size_t line, Matrix::Element number) {
	for (std::size_t j = 0; j < mat.GetColumnCount(); j++) {
		mat.at(line, j) /= number;
	}
}

static int FirstNotNullElementIndexOnColumn(Matrix& mat, std::size_t column, std::size_t startLine = 0) {
	for (std::size_t i = startLine; i < mat.GetRawCount(); i++) {
		if (!IsEqualZero(mat.at(i, column))) {
			return i;
		}
	}
	return -1;
}

static void SimplifyLine(Matrix& mat, std::size_t line, std::size_t pivot_line, std::size_t pivot_column) {
	const Matrix::Element pivot = mat.at(pivot_line, pivot_column);
	const Matrix::Element anul = mat.at(line, pivot_column);

	for (std::size_t j = 0; j < mat.GetColumnCount(); j++) {
		mat.at(line, j) = mat.at(line, j) * pivot - mat.at(pivot_line, j) * anul;
	}
}

void GaussJordan(Matrix& mat, bool reduite, bool normalise) {
	int indice_ligne_pivot = -1;

	for (std::size_t j = 0; j < mat.GetColumnCount(); j++) {

		int indice_ligne_pivot_trouve = FirstNotNullElementIndexOnColumn(mat, j, indice_ligne_pivot + 1);

		if (indice_ligne_pivot_trouve < 0)	// colonne de 0
			continue;						// on regarde la prochaine colonne

		indice_ligne_pivot++;

		if (indice_ligne_pivot_trouve != indice_ligne_pivot) {
			SwapLines(mat, indice_ligne_pivot_trouve, indice_ligne_pivot);
		}

		Matrix::Element pivot = mat.at(indice_ligne_pivot, j);

		if (normalise) {
			DivideLine(mat, indice_ligne_pivot, pivot);
		}

		// On simplifie les autres lignes
		for (std::size_t i = (reduite ? 0 : j); i < mat.GetRawCount(); i++) {
			// Pour les lignes autre que la ligne pivot
			if (i != static_cast<std::size_t>(indice_ligne_pivot)) {
				SimplifyLine(mat, i, indice_ligne_pivot, j);
			}
		}
	}
}

}  // namespace Gauss