#!/usr/bin/python3
# Élimination de Gauss

import numpy as np
    
def triangulariser(M):  # Transformer une matrice n fois (n+1) en forme triangulaire supérieure
    n = M.shape[0]      # le nombre de lignes
    for i in range(n):  # une boucle sur toutes les colonnes sauf la dernière
        for k in range(i, n): # boucle sur les éléments sous de la diagonale pour trouver un pivot
            if M[k, i] != 0:        # ... cet élément va-t-il le faire?
                M[[i, k], :] = M[[k, i], :]   # si oui: échanger les lignes i et k
                pivot = M[i, i]               # ...enrégistrer sa valeur...
                break                         # ...et arrêter la boucle sur k
        else:                   # il n'y a pas d'élément non nul sous la diagonale?
            continue                # alors rien à faire pour cette colonne
        for k in range(i+1, n): # éliminer tous les éléments au-dessous de la diagonale:
            facteur = -M[k, i]/pivot           # ajouter (facteur) ...
            M[k, :] += facteur * M[i, :]       #  ... * (ligne du pivot) à la ligne k.
    

def gauss(A, b):  # Trouver la solution x du système linéaire Ax = b
    n = A.shape[0]         # le nombre de lignes
    M = np.empty((n, n+1)) # la matrice M
    M[:, :n] = np.copy(A)  # copier A dans les premières n colonnes
    M[:, n] = np.copy(b)   # copier b dans la dernière colonne
    triangulariser(M)  # après, M est triangulaire supérieure
    
    x = np.empty(n)    # on mettra la solution ici
    for i in range(n-1, -1, -1):    # boucle en arrière sur les lignes, départ = dernière ligne
        sigma = 0.                    # la somme des x que l'on connaît déjà fois les Mik correspondants
        for k in range(i+1, n):     # boucle sur les éléments non nul
            sigma += M[i, k] * x[k]
        # On pourrait remplacer ces dernières 3 lignes par une seule ligne, plus élégant mais moins clair:
        # sigma = M[i, i+1:n] @ x[i+1:]
        if M[i, i] == 0:
            if M[i, n] - sigma == 0:
                print("Attention, solution pas unique!")
                x[i] = 42
            else:
                print("Erreur, pas de solution")
                return
        else:
            x[i] = (M[i, n] - sigma)/M[i, i]    
    return x

# tester le code
if __name__ == "__main__":
    A = np.array([[5,3,-1],[-2,2,-4],[0,-1,1]], dtype=float)
    b = np.array([0, 1, 3], dtype=float)
    print("Solution x = ", gauss(A, b))
    print("A.x = ", A @ gauss(A, b))

    A = np.array([[0,2,-2,2],[-1,3,4,2],[1,-1,0,-3],[1,1,-2,0]], dtype=float)
    b = np.array([4,2,2,0], dtype=float)
    print("Solution x = ", gauss(A, b))
    print("A.x = ", A @ gauss(A, b))