import numpy as np

def echelonize(a, eps = 1e-8):
    m = a.shape[0]
    n = a.shape[1]
    i = 0
    j = 0
    while i < m and j < n:
        # print("Iteration:", i, j)    
        # print(a)
        # Find maximal element in j-th column, starting from i-th row
        indMax = i
        vMax = abs(a[i, j])
        for k in range(i+1, m):
            v = abs(a[k, j])
            if v > vMax:
                vMax = v
                indMax = k
        if vMax <= eps:
            for k in range(i, m):
                a[k, j] = 0.
            j += 1
            continue
        # print("max. element in column", i, "is in row", indMax)
        if indMax != i:
            # Swap rows i, indMax
            a[[i, indMax]] = a[[indMax, i]]
            a[indMax] = -a[indMax]  # To save a sign of determinant            
            # print("swapping rows", i, indMax)
            # print(a)
        # Annihilate j-th column starting from i+1-th row
        r = (-1./a[i, j])
        for k in range(i+1, m):
            a[k] += a[i]*(a[k, j]*r)
            a[k, j] = 0.
        # print("annihilating column", j)
        # print(a)
        i += 1
        j += 1
    return i                
            
def reversePass(a, r, eps = 1e-8):
    '''a is the matrix in the echelonized form,
    r is the rank of matrix a.
    Reduce the matrix a to diagonal form with units at the beginnings of rows
    and seroes in the columns above the first nonzero elements of rows'''
    m = a.shape[0]      # Number of rows
    n = a.shape[1]      # Number of columns
    for i in range(r-1, -1, -1):
        # Find the first nonzero element in i-th row
        j0 = n
        for j in range(n):
            if abs(a[i, j]) > eps:
                j0 = j
                break
        if j0 >= n:        
            assert "Incorrect echelon form of matrix"                
        # multiplyMatrixRow(a, i, 1./a[i][j0])
        a[i] *= (1./a[i, j0])
        for k in range(i):
            # addMatrixRows(a, k, i, -a[k][j0])
            a[k] -= a[i]*a[k, j0]
            
def inverseMatrix(a, eps = 1e-8):
    '''Return the inverse matrix to the matrix a of real numbers.
    The matrix a is represented as 2-dimensional numpy-array'''

    # To do...
    pass
            
def solveLinearSystem(a, b, eps = 1e-8):
    '''Solve a linear system a*x = b,
    where the matrix a is represented by a 2-dimensional numpy-array
    of real numbers, and b is a linear numpy-array or any list/tuple'''

    # To do...
    pass

def determinant(a, eps = 1e-8):
    m = a.shape[0]
    n = a.shape[1]
    if m != n:
        raise ValueError("Determinant of non-square matrix")
    b = a.copy()
    r = echelonize(b, eps)
    if r < m:
        return 0.
    d = 1.
    for i in range(m):
        d *= b[i, i]
    return d