import numpy as np

def echelonize(a, eps=1e-8):
    """Reduce a matrix to row echelon form.
    Return the rank of matrix"""
    m, n = a.shape  # number of rows, columns
    i = 0; j = 0
    while i < m and j < n:
        # Find the maximal element in j-th column,
        # starting from i-th row
        k = i; maxVal = abs(a[i, j])
        for l in range(i+1, m):
            if abs(a[l, j]) > maxVal:
                maxVal = abs(a[l, j])
                k = l
        if maxVal <= eps:
            # Zero column
            for l in range(i, m):
                a[l, j] = 0.
            j += 1
            continue
        if k != i:
            a[[i, k]] = a[[k, i]]   # Swap rows i, k. To save a determinant,
            a[k] *= (-1.)           # change a sign of one row.
        r = a[i, j]
        for l in range(i+1, m):
            a[l] -= a[i]*(a[l, j]/r)
            a[l, j] = 0.
        i += 1
        j += 1
    return i    # Return the rank of matrix = number of nonzero rows

def determinant(a):
    if a.shape[0] != a.shape[1]:
        raise ValueError("determinant of non-square matrix")
    b = a.copy()
    r = echelonize(b)   # r is the rank of the matrix
    if r < b.shape[0]:
        # print("Singular matrix")
        return 0.   # Determinant of singular matrix = 0
    d = 1.
    for i in range(b.shape[0]):
        d *= b[i, i]
    return d

def reversePass(a, r, eps = 1e-8):
    '''a is an echelonized matrix, r is its rank.
    Annihilate elements in columns above the
    first nonzero elements of rows, and normalize
    rows so that first nonzero elements = 1.
    Example:
    [2 3 0 1]       [1 0 -7.5 0]
    [0 1 5 2]  -->  [0 1   5  0]
    [0 0 0 1]       [0 0   0  1]
    '''
    m, n = a.shape
    assert r <= m
    for i in range(r-1, -1, -1):
        j = 0
        while j < n and abs(a[i, j]) <= eps:
            j += 1
        if j >= n:
            continue
        assert abs(a[i, j]) > eps
        a[i] *= 1./a[i, j]
        for k in range(i):
            a[k] -= a[i]*a[k, j]
    
def invMatrix(a, eps = 1e-8):
    """Return an inverse matrix of a matrix a"""
    pass        # To do...

def solveLinearSystem(a, b, eps = 1e-8):
    """Solve a linear system
        a x = b
    where a is non-singular square matrix
    and b is array of free terms. Return an array x"""
    pass        # To do...