echelon.py

from fractions import Fraction
import numpy as np


def print_matrix(matrix):
    def tostr(i):
        return str(i) if i.denominator == 1 else '{}\\{}'.format(i.numerator, i.denominator)
    print '\n'.join(['\t'.join([tostr(col) for col in row]) for row in matrix])


def to_latex(matrix, line):
    if not isinstance(matrix, np.ndarray):
        matrix = np.array(matrix)
    def tostr(i):
        return str(i) if i.denominator == 1 else '\\tfrac{{{}}}{{{}}}'.format(i.numerator, i.denominator)
    out = '\\Rightarrow \\begin{bmatrix}' + '\\\\'.join([' & '.join([tostr(col) for col in row]) for row in matrix]) + '\\end{bmatrix}'
    if line % 4 == 1:
        out = '&' + out
    if line % 4 == 0:
        out += '\\\\'
    print out
    return line + 1


def rref(matrix):
    if not isinstance(matrix, np.ndarray):
        matrix = np.array(matrix)
    matrix = np.vectorize(Fraction)(matrix)
    h = 0
    k = 0
    pivots = []
    line = 1
    # Row echelon form
    while h < matrix.shape[0] and k < matrix.shape[1]:
        # Find max abs pivot for numerical stability, not really needed since using fractions
        i_max = np.argmax(np.abs(matrix[h:, k])) + h
        line = to_latex(matrix, line)
        if matrix[i_max, k] == 0:
            k += 1
        else:
            matrix[[h, i_max]] = matrix[[i_max, h]]  # Swap rows
            for i in xrange(h, matrix.shape[0]):
                if matrix[i, k] != 0:
                    matrix[i] /= matrix[i, k]
                    if i > h:
                        matrix[i] -= matrix[h]
            pivots.append(k)
            h += 1
            k += 1
    # Reduced row echelon form
    for i in xrange(len(pivots)):
        line = to_latex(matrix, line)
        for j in xrange(i + 1, len(pivots)):
            matrix[i] -= matrix[j] * matrix[i, pivots[j]]
    return matrix


if __name__ == '__main__':
    A = [[17, 1, 5, 0], [6, 2, 4, 0], [-13, 3, 3, 0], [1, 4, 2, 0], [2, 5, 1, 0]]
    rref(A)