import numpy

def solve_triangular(A, b):
    n = len(b)
    x = numpy.zeros([n])
    for i in range(n - 1, -1, -1):
        x[i] = (b[i] - numpy.dot(A[i, :], x)) / A[i, i]
    return x

def pivot_index(A,i):
    n = A.shape[0] # nombre de lignes
    j = i
    for k in range(i+1,n):
        if abs(A[k,i]) > abs(A[j,i]):
            j = k
    return j

def swap_lines(A, i, j):
    tmp = A[i, :].copy()
    A[i, :] = A[j, :]
    A[j, :] = tmp

def transvection_lines(A, i, k, factor):
    A[k, :] = A[k, :] + A[i, :] * factor

def gauss(A0, b0):
    A = A0.copy() # pour ne pas détruire A
    b = b0.copy()
    n = len(b)  # on pourrait aussi verifier que a a les bonnes dimensions
    for i in range(n - 1):
        # choix du pivot
        ipiv = pivot_index(A, i)
        # echanges
        if ipiv != i:
            swap_lines(A, i, ipiv)
            tmp = b[ipiv]
            b[ipiv] = b[i]
            b[i] = tmp
        # pivotage
        for k in range(i + 1, n):
            factor = -A[k, i] / A[i, i]
            transvection_lines(A, i, k, factor)
            b[k] = b[k] + b[i] * factor
    return [A, b] # une liste [matrice,vec]

def solve(A0, b0):
    Ab = gauss(A0, b0)
    #Ab[0] contient la mat triang sup
    #Ab[1] contient le vec apres gauss
    return solve_triangular(Ab[0], Ab[1])

n = 5
M = numpy.random.random([n,n])
x = numpy.random.random([n])
v = numpy.dot(M,x)
x_solved = solve(M,v)
abs(x - x_solved).max()