#include "Gauss.h"

#include "Matrix.h"

#include <algorithm>
#include <execution>
#include <ranges>

namespace Gauss {

static void SwapLines(Matrix& mat, std::size_t line1, std::size_t line2) {
	std::swap_ranges(
		std::execution::par_unseq, mat.GetLineIterator(line1), mat.GetLineIterator(line1 + 1), mat.GetLineIterator(line2));
}

static void DivideLine(Matrix& mat, std::size_t line, Matrix::Element number) {
	std::transform(std::execution::par_unseq, mat.GetLineIterator(line), mat.GetLineIterator(line + 1), mat.GetLineIterator(line),
		[number](Matrix::Element e) { return e /= 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);

	auto range = std::views::iota(static_cast<std::size_t>(0), mat.GetColumnCount());

	std::for_each(std::execution::par_unseq, range.begin(), range.end(), [&mat, pivot, anul, line, pivot_line](std::size_t j) {
		mat.at(line, j) = mat.at(line, j) * pivot - mat.at(pivot_line, j) * anul;
	});
}

static void GaussJordanReduced(Matrix& a_Matrix, bool a_Normalise) {
	int indice_ligne_pivot = -1;

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

		int indice_ligne_pivot_trouve = FirstNotNullElementIndexOnColumn(a_Matrix, 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(a_Matrix, indice_ligne_pivot_trouve, indice_ligne_pivot);
		}

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

		if (a_Normalise) {
			DivideLine(a_Matrix, indice_ligne_pivot, pivot);
		}

		auto range = std::views::iota(static_cast<std::size_t>(0), a_Matrix.GetRawCount());

		// On simplifie les autres lignes
		std::for_each(std::execution::par_unseq, range.begin(), range.end(), [&a_Matrix, j, indice_ligne_pivot](std::size_t i) {
			if (i != static_cast<std::size_t>(indice_ligne_pivot)) {
				SimplifyLine(a_Matrix, i, indice_ligne_pivot, j);
			}
		});
	}
}

static void GaussJordanTriangular(Matrix& a_Matrix, bool a_Normalise) {
	int indice_ligne_pivot = -1;

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

		int indice_ligne_pivot_trouve = FirstNotNullElementIndexOnColumn(a_Matrix, 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(a_Matrix, indice_ligne_pivot_trouve, indice_ligne_pivot);
		}

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

		if (a_Normalise) {
			DivideLine(a_Matrix, indice_ligne_pivot, pivot);
		}

		auto range = std::views::iota(static_cast<std::size_t>(indice_ligne_pivot + 1), a_Matrix.GetRawCount());

		// On simplifie les autres lignes après la ligne du pivot
		std::for_each(std::execution::par_unseq, range.begin(), range.end(),
			[&a_Matrix, indice_ligne_pivot, j](std::size_t i) { 
				SimplifyLine(a_Matrix, i, indice_ligne_pivot, j); 
			});
	}
}

void GaussJordan(Matrix& a_Matrix, bool a_Reduite, bool a_Normalise) {
	if (a_Reduite)
		GaussJordanReduced(a_Matrix, a_Normalise);
	else
		GaussJordanTriangular(a_Matrix, a_Normalise);
}

}  // namespace Gauss