# -*- coding: utf-8 -*-
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


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

if __name__ == '__main__':
    n = 5
    M = numpy.random.random([n, n])
    x = numpy.random.random([n])
    v = numpy.dot(M, x)
    x_solved = solve(M, v)
    print(abs(x - x_solved).max())