# Alfonso Cevallos. May 19 2015.
# Please report any mistakes.
# We work with Python numbers, and built-in elementary operations.
# We implement the simplified HNF algorithm, where the numbers can
# potentially grow exponentially large.
# The goal here is simply to understand the higher level idea of the algorithms.

import math
import os
import random
import time


def div(a,b): #computes (q,r) such that a = qb + r, with r < b
    r = a%b
    q = (a-r)/b
    return q,r


def exgcd(a,b): #computes (g,x,y) such that gcd(a,b) = g = xa + yb
    if a<b:
        g,x,y = exgcd(b,a)
        return g,y,x
    if b<0:
        g,x,y = exgcd(a,-b)
        return g,x,-y
    if b==0: return a,1,0

    q,r = div(a,b)
    g,x,y = exgcd(b,r)
    return g,y,x-y*q


def matrixMult(a1,a2,b1,b2):
    # Computes the matrix multiplication C = A*B,
    # where A = [a1, a2], B = [b1, b2], C = c1, c2.
    # Returns c1, c2.
    # The column vectors a1, a2, c1, c2 are of the same length,
    # and the column vectors b1, b2 are of length 2.

    length = len(a1)
    c1, c2 = [], []
    for i in range(length):
        e1 = a1[i]*b1[0] + a2[i]*b1[1]
        c1.append(e1)
        e2 = a1[i]*b2[0] + a2[i]*b2[1]
        c2.append(e2)
    return c1, c2


def HNF(A):
    # Takes as input a full row-rank matrix A (returns -1 if not full row rank).
    # A is m x n. Remark: A[j][i] corresponds to entry in i-th row, j-th column.
    # In other words, we see matrices as a list of column vectors, indexed by j.
    # Algorithm returns H and U, such that AU = H = [triangle|0] is HNF of A.

    n = len(A)
    m = len(A[0])
    if n < m: return -1

    H, U = [], []    # we start by defining H = S and U = I_n
    for j in range(n):
        colU = [0]*n
        colU[j] = 1
        U.append(colU)
        colH = A[j][:]
        H.append(colH)

    for i in range(m):  #main loop over every row

        for j in range(i+1,n):  # make zeroes in position after diagonal, using gcd.
            if H[j][i]!=0:
                g,x,y = exgcd(H[i][i], H[j][i])
                b1, b2 = [x, y], [-H[j][i]/g, H[i][i]/g]
                H[i], H[j] = matrixMult(H[i], H[j], b1, b2)
                U[i], U[j] = matrixMult(U[i], U[j], b1, b2)

        for j in range(i):      # reduce the positions before diagonal, using div with remainder
            q,r = div(H[j][i], H[i][i])
            b1, b2 = [1, -q], [0, 1]
            H[j], H[i] = matrixMult(H[j], H[i], b1, b2)
            U[j], U[i] = matrixMult(U[j], U[i], b1, b2)

    return H,U


def Tests():
    # The following two matrices are given as examples on the lecture notes 9 of last year
    # Link to this file: http://disopt.epfl.ch/files/content/sites/disopt/files/shared/Cal14/Lecture-09.pdf.pdf
    # Remark: All matrices are printed transposed. So columns are printed as rows.

    A1 = [[2,2],[3,4],[4,6]]
    H1, U1 = HNF(A1)

    A2 = [[4,-8], [3,9], [2,1]]
    H2, U2 = HNF(A2)

    print "For matrix A:"
    for a in A1:
        print a
    print "\nIts HNF is:"
    for h in H1:
        print h
    print "\nAnd the corresponding unimodular matrix U is:"
    for u in U1:
        print u

    print "\n\nAnd for matrix A:"
    for a in A2:
        print a
    print "\nIts HNF is:"
    for h in H2:
        print h
    print "\nAnd the corresponding unimodular matrix U is:"
    for u in U2:
        print u


Tests()