matrix_integer_dense_hnf.py

"""
Modular algorithm to compute Hermite normal forms of integer matrices.

AUTHORS:
    -- Clement Pernet and William Stein (2008-02-07): initial version
"""

from copy import copy

from sage.misc.misc import verbose, cputime
from sage.matrix.constructor import random_matrix, matrix, matrix, identity_matrix

from sage.rings.all import ZZ, previous_prime, next_prime, CRT_list, RR
import math

def max_det_prime(n):
    """
    Return the largest prime so that it is reasonably efficiency to
    compute modulo that prime with n x n matrices in LinBox.

    INPUT:
        n -- a positive integer
        
    OUTPUT:
        a prime number

    EXAMPLES:
        sage: from sage.matrix.matrix_integer_dense_hnf import max_det_prime
        sage: max_det_prime(10000)
        524309
        sage: max_det_prime(1000)
        2097169
        sage: max_det_prime(10)
        16777259
    """
    k = int(26 - math.ceil(math.log(n)*0.7213475205))
    return next_prime(2**k)

def det_from_modp_and_divisor(A, d, p, z_mod, moduli, z_so_far=ZZ(1), N_so_far=ZZ(1)):
    """
    This is used for internal purposes for computing determinants
    quickly (with the hybrid p-adic / multimodular algorithm).
    
    INPUT:
        A -- a square matrix
        d -- a divisor of the determinant of A
        p -- a prime
        z_mod -- values of det/d (mod ...)
        moduli -- the moduli so far
        z_so_far -- for a modulus p in the list moduli,
            (z_so_far mod p) is the determinant of A modulo p.
        N_so_far -- N_so_far is the product over the primes in the list moduli.

    OUTPUT:
        A triple (det bound, new z_so_far, new N_so_far).
        
    EXAMPLES:
        sage: a = matrix(ZZ, 3, [6, 1, 2, -56, -2, -1, -11, 2, -3])
        sage: factor(a.det())
        -1 * 13 * 29
        sage: d = 13
        sage: import sage.matrix.matrix_integer_dense_hnf as matrix_integer_dense_hnf
        sage: matrix_integer_dense_hnf.det_from_modp_and_divisor(a, d, 97, [], [])
        (-377, -29, 97)
        sage: a.det()
        -377
    """
    tm = verbose("Multimodular stage of det calculation -- using p = %s"%p, level=2)
    z = A._linbox_modn_det(p) / d
    z = z.lift()
    z_mod.append(z)
    moduli.append(p)
    z = CRT_list([z_so_far, z], [N_so_far, p])
    N = N_so_far*p

    if z > N//2:
        z = z - N
    verbose("Finished multimodular det for p = %s"%p, tm, level=2)
    return (d * z, z, N)

def det_given_divisor(A, d, proof=True, stabilize=2):
    """
    Given a divisor d of the determinant of A, compute the
    determinant of A.

    INPUT:
        A -- a square integer matrix
        d -- a nonzero integer that is assumed to divide the determinant of A
        proof -- bool (default True) compute det modulo enough primes
                 so that the determinant is computed provably correctly
                 (via the Hadamard bound).  It would be VERY hard for
                 det to fail even with proof=False. 
        stabilize -- int (default: 2) if proof = False, then compute det
                 mod p until stabilize successive modulo det computations
                 stabilize.

    OUTPUT:
        integer -- determinant

    EXAMPLES:
        sage: import sage.matrix.matrix_integer_dense_hnf as matrix_integer_dense_hnf
        sage: a = matrix(ZZ,3,[-1, -1, -1, -20, 4, 1, -1, 1, 2])
        sage: matrix_integer_dense_hnf.det_given_divisor(a, 3)
        -30
        sage: matrix_integer_dense_hnf.det_given_divisor(a, 3, proof=False)
        -30
        sage: matrix_integer_dense_hnf.det_given_divisor(a, 3, proof=False, stabilize=1)
        -30
        sage: a.det()
        -30

    Here we illustrate proof=False giving a wrong answer:
        sage: p = matrix_integer_dense_hnf.max_det_prime(2)
        sage: q = previous_prime(p)
        sage: a = matrix(ZZ, 2, [p, 0, 0, q])
        sage: p * q
        1125899772623531
        sage: matrix_integer_dense_hnf.det_given_divisor(a, 1, proof=False, stabilize=2)
        0

    This still works, because we don't work modulo primes that divide
    the determinant bound, which is found using a p-adic algorithm.
        sage: a.det(proof=False, stabilize=2)
        1125899772623531


    3 primes is enough:
        sage: matrix_integer_dense_hnf.det_given_divisor(a, 1, proof=False, stabilize=3)
        1125899772623531
        sage: matrix_integer_dense_hnf.det_given_divisor(a, 1, proof=False, stabilize=5)
        1125899772623531
        sage: matrix_integer_dense_hnf.det_given_divisor(a, 1, proof=True)
        1125899772623531

    TESTS:
        sage: m = diagonal_matrix(ZZ, 68, [2]*66 + [1,1])
        sage: m.det()
        73786976294838206464
    """
    p = max_det_prime(A.nrows())
    z_mod = []
    moduli = []
    assert d != 0
    z_so_far = 1
    N_so_far = 1
    if proof:
        N = 1
        B = (2 * 10**A.hadamard_bound()) // d + 1
        dd = d
        # bad verbose statement, since computing the log overflows!
        est = int(RR(B).log() / RR(p).log()) + 1
        cnt = 1
        verbose("Multimodular det -- need to use about %s primes."%est, level=1)
        while N < B:
            if d % p != 0:
                tm = cputime()
                dd, z_so_far, N_so_far = det_from_modp_and_divisor(A, d, p, z_mod, moduli, z_so_far, N_so_far)
                N *= p
                verbose("computed det mod p=%s which is %s (of about %s)"%(p, cnt, est), tm)
            p = previous_prime(p)
            cnt += 1            
        return dd
    else:
        val = []
        while True:
            if d % p != 0:
                tm = cputime()
                dd, z_so_far, N_so_far = det_from_modp_and_divisor(A, d, p, z_mod, moduli, z_so_far, N_so_far)
                verbose("computed det mod %s"%p, tm)
                val.append(dd)
                if len(val) >= stabilize and len(set(val[-stabilize:])) == 1:
                    return val[-1]
            p = previous_prime(p)

def det_padic(A, proof=True, stabilize=2):
    """
    Return the determinant of A, computed using a p-adic/multimodular
    algorithm.

    INPUTS:
        A -- a square matrix
        proof -- boolean
        stabilize (default: 2) -- if proof False, number of successive primes so that
                     CRT det must stabilize.

    EXAMPLES:
        sage: import sage.matrix.matrix_integer_dense_hnf as h
        sage: a = matrix(ZZ, 3, [1..9])
        sage: h.det_padic(a)
        0
        sage: a = matrix(ZZ, 3, [1,2,5,-7,8,10,192,5,18])
        sage: h.det_padic(a)
        -3669
        sage: a.determinant(algorithm='ntl')
        -3669
    """
    if not A.is_square():
        raise ValueError("A must be a square matrix")
    r = A.rank()
    if r < A.nrows():
        return ZZ(0)
    v = random_matrix(ZZ, A.nrows(), 1)
    d = A._solve_right_nonsingular_square(v, check_rank=False).denominator()
    return det_given_divisor(A, d, proof=proof, stabilize=stabilize)
    
def double_det (A, b, c, proof):
    """
    Compute the determinants of the stacked integer matrices
    A.stack(b) and A.stack(c).

    INPUT:
        A -- an (n-1) x n matrix
        b -- an 1 x n matrix
        c -- an 1 x n matrix
        proof -- whether or not to compute the det modulo enough times
                 to provably compute the determinant.

    OUTPUT:
        a pair of two integers.

    EXAMPLES:
        sage: from sage.matrix.matrix_integer_dense_hnf import double_det
        sage: A = matrix(ZZ, 2, 3, [1,2,3, 4,-2,5])
        sage: b = matrix(ZZ, 1, 3, [1,-2,5])
        sage: c = matrix(ZZ, 1, 3, [8,2,10])
        sage: A.stack(b).det()
        -48
        sage: A.stack(c).det()
        42
        sage: double_det(A, b, c, False)
        (-48, 42)
    """
    # We use the "two for the price of one" algorithm, which I made up. (William Stein)

    # This is a clever trick!  First we transpose everything.  Then
    # we use that if [A|b]*v = c then [A|c]*w = b with w easy to write down!
    # In fact w is got from v by dividing all entries by -v[n], where n is the
    # number of rows of v, and *also* dividing the last entry of w by v[n] again.
    # See this as an algebra exercise where you have to think of matrix vector
    # multiply as "linear combination of columns".
    A = A.transpose()
    b = b.transpose()
    c = c.transpose()
    t = verbose('starting double det')
    B = A.augment(b)
    v = B.solve_right(-c)

    db = det_given_divisor(B, v.denominator(), proof=proof)

    n = v.nrows()
    vn = v[n-1,0]
    w = (-1/vn)*v
    w[n-1] = w[n-1]/vn
    dc = det_given_divisor(A.augment(c), w.denominator(), proof=proof)

    verbose('finished double det', t)

    return (db, dc)

def add_column_fallback(B, a, proof):
    """
    Simplistic version of add_column, in case the powerful clever one
    fails (e.g., B is singular).  

    INPUT:
        B -- a square matrix (may be singular)
        a -- an n x 1 matrix, where B has n rows
        proof -- bool; whether to prove result correct

    OUTPUT:
        x   -- a vector such that H' = H_B.augment(x) is the HNF of A = B.augment(a).

    EXAMPLES:
        sage: B = matrix(ZZ,3, [-1, -1, 1, -3, 8, -2, -1, -1, -1])
        sage: a = matrix(ZZ,3,1, [1,2,3])
        sage: import sage.matrix.matrix_integer_dense_hnf as matrix_integer_dense_hnf
        sage: matrix_integer_dense_hnf.add_column_fallback(B, a, True)
        [-3]
        [-7]
        [-2]
        sage: matrix_integer_dense_hnf.add_column_fallback(B, a, False)
        [-3]
        [-7]
        [-2]
        sage: B.augment(a).hermite_form()
        [ 1  1  1 -3]
        [ 0 11  1 -7]
        [ 0  0  2 -2]
    """
    tt = verbose('add column fallback...')
    W = B.augment(matrix(ZZ,B.nrows(),a.list()))
    H, _ = hnf(W, proof)
    C = H.matrix_from_columns([H.ncols()-1])
    verbose('finished add column fallback', tt)
    return C

def solve_system_with_difficult_last_row(B, a):
    """
    Solve B*x = a when the last row of $B$ contains huge entries using
    a clever trick that reduces the problem to solve C*x = a where $C$
    is $B$ but with the last row replaced by something small, along
    with one easy null space computation.  The latter are both solved
    $p$-adically.
    
    INPUT:
        B -- a square n x n nonsingular matrix with painful big bottom row.
        a -- an n x 1 column matrix
    OUTPUT:
        the unique solution to B*x = a.

    EXAMPLES:
        sage: from sage.matrix.matrix_integer_dense_hnf import solve_system_with_difficult_last_row
        sage: B = matrix(ZZ, 3, [1,2,4, 3,-4,7, 939082,2930982,132902384098234])
        sage: a = matrix(ZZ,3,1, [1,2,5])
        sage: z = solve_system_with_difficult_last_row(B, a)
        sage: z
        [ 106321906985474/132902379815497]
        [132902385037291/1329023798154970]
        [        -5221794/664511899077485]
        sage: B*z
        [1]
        [2]
        [5]
    """
    # Here's how:
    # 1. We make a copy of B but with the last *nasty* row of B replaced
    #    by a random very nice row.
    C = copy(B)
    while True:
        C[C.nrows()-1] = random_matrix(ZZ,1,C.ncols()).row(0)
        # 2. Then we find the unique solution to C * x = a
        try:
            x = C.solve_right(a)
        except ValueError:
            verbose("Try difficult solve again with different random vector")
        else:
            break
        

    # 3. We next delete the last row of B and find a basis vector k
    #    for the 1-dimensional kernel.
    D = B.matrix_from_rows(range(C.nrows()-1))
    N = D._rational_kernel_iml()
    if N.ncols() != 1:
        verbose("Try difficult solve again with different random vector")
        return solve_system_with_difficult_last_row(B, a)
    
    k = N.matrix_from_columns([0])

    # 4. The sought for solution z to B*z = a is some linear combination
    #
    #       z = x + alpha*k
    #
    #  of x and k, where k is the above fixed basis for the kernel of D.
    #  Setting w to be the last row of B, this column vector z satisfies
    #
    #       w * z = a'
    #
    # where a' is the last entry of a.  Thus
    #
    #       w * (x + alpha*k) = a'
    #
    # so    w * x + alpha*w*k = a'
    # so    alpha*w*k  = a' - w*x.

    w = B[-1]  # last row of B
    a_prime = a[-1]
    lhs = w*k
    rhs = a_prime - w * x
    
    if lhs[0] == 0:
        verbose("Try difficult solve again with different random vector")
        return solve_system_with_difficult_last_row(B, a)
    
    alpha = rhs[0] / lhs[0]
    z = x + alpha*k
    return z
    
def add_column(B, H_B, a, proof):
    """
    The add column procedure.

    INPUT:
        B   -- a square matrix (may be singular)
        H_B -- the Hermite normal form of B
        a -- an n x 1 matrix, where B has n rows
        proof -- bool; whether to prove result correct, in case we use fallback method.
        
    OUTPUT:
        x   -- a vector such that H' = H_B.augment(x) is the HNF of A = B.augment(a).

    EXAMPLES:
        sage: B = matrix(ZZ, 3, 3, [1,2,5, 0,-5,3, 1,1,2])
        sage: H_B = B.echelon_form()
        sage: a = matrix(ZZ, 3, 1, [1,8,-2])
        sage: import sage.matrix.matrix_integer_dense_hnf as hnf
        sage: x = hnf.add_column(B, H_B, a, True); x
        [18]
        [ 3]
        [23]
        sage: H_B.augment(x)
        [ 1  0 17 18]
        [ 0  1  3  3]
        [ 0  0 18 23]
        sage: B.augment(a).echelon_form()
        [ 1  0 17 18]
        [ 0  1  3  3]
        [ 0  0 18 23]
    """
    t0 = verbose('starting add_column')

    if B.rank() < B.nrows():
        return add_column_fallback(B, a, proof)
    else:
        z = solve_system_with_difficult_last_row(B, a)
        
    zd, d = z._clear_denom()
    x = H_B * zd
    if d != 1:
        for i in range(x.nrows()):
            x[i,0] = x[i,0]/d

    return x

def add_row(A, b, pivots, include_zero_rows):
    """
    The add row procedure.

    INPUT:
        A -- a matrix in Hermite normal form with n column
        b -- an n x 1 row matrix
        pivots -- sorted list of integers; the pivot positions of A.

    OUTPUT:
        H -- the Hermite normal form of A.stack(b).
        new_pivots -- the pivot columns of H.

    EXAMPLES:
        sage: import sage.matrix.matrix_integer_dense_hnf as hnf
        sage: A = matrix(ZZ, 2, 3, [-21, -7, 5, 1,20,-7])
        sage: b = matrix(ZZ, 1,3, [-1,1,-1])
        sage: hnf.add_row(A, b, A.pivots(), True)
        (
        [ 1  6 29]           
        [ 0  7 28]           
        [ 0  0 46], [0, 1, 2]
        )
        sage: A.stack(b).echelon_form()
        [ 1  6 29]
        [ 0  7 28]
        [ 0  0 46]
    """
    t = verbose('add hnf row')
    v = b.row(0)
    H, pivs = A._add_row_and_maintain_echelon_form(b.row(0), pivots)
    if include_zero_rows and H.nrows() != A.nrows()+1:
        H = H.matrix_from_rows(range(A.nrows()+1))
    verbose('finished add hnf row', t)
    return H, pivs

def pivots_of_hnf_matrix(H):
    """
    Return the pivot columns of a matrix H assumed to be in HNF.

    INPUT:
        H -- a matrix that must be HNF

    OUTPUT:
        list -- list of pivots

    EXAMPLES:
        sage: H = matrix(ZZ, 3, 5, [1, 0, 0, 45, -36, 0, 1, 0, 131, -107, 0, 0, 0, 178, -145]); H
        [   1    0    0   45  -36]
        [   0    1    0  131 -107]
        [   0    0    0  178 -145]
        sage: import sage.matrix.matrix_integer_dense_hnf as matrix_integer_dense_hnf             
        sage: matrix_integer_dense_hnf.pivots_of_hnf_matrix(H)
        [0, 1, 3]
    """
    pivots = []
    r = -1
    for j in xrange(H.ncols()):
        # Find first nonzero position (counting from bottom) in the j-th column
        for i in reversed(xrange(H.nrows())):
            if H[i,j]:
                if i > r:
                    pivots.append(j)
                    r = i
                elif i <= r:
                    break
    return pivots

def hnf_square(A, proof):
    """
    INPUT:
        a nonsingular n x n matrix A over the integers.
    OUTPUT:
        the Hermite normal form of A.

    EXAMPLES:
        sage: import sage.matrix.matrix_integer_dense_hnf as hnf
        sage: A = matrix(ZZ, 3, [-21, -7, 5, 1,20,-7, -1,1,-1])
        sage: hnf.hnf_square(A, False)
        [ 1  6 29]
        [ 0  7 28]
        [ 0  0 46]
        sage: A.echelon_form()
        [ 1  6 29]
        [ 0  7 28]
        [ 0  0 46]
    """
    n = A.nrows()
    m = A.ncols()
    if n != m:
        raise ValueError("A must be square.")

    # Small cases -- don't use this algorithm
    if n <= 3:
        return A.echelon_form(algorithm="pari")

    if A.rank() < A.nrows():
        raise ValueError("matrix must have full rank")
                           
    

    t = verbose("starting slicings")
    B = A.matrix_from_rows(range(m-2)).matrix_from_columns(range(n-1))
    c = A.matrix_from_rows([m-2]).matrix_from_columns (range(n-1))
    d = A.matrix_from_rows([m-1]).matrix_from_columns (range(n-1))
    b = A.matrix_from_columns([n-1]).matrix_from_rows(range(m-2))
    verbose("done slicing", t)

    try:
        (d1,d2) = double_det (B,c,d, proof=proof)
    except (ValueError, ZeroDivisionError), msg:
        d1 = B.stack(c).det(proof=proof)
        d2 = B.stack(d).det(proof=proof)
        
    (g,k,l) = d1._xgcd (d2, minimal=True)
    
    W = B.stack (k*c + l*d)
    verbose("submatrix det: g=%s"%g)
    CUTOFF = 2**30
    if g == 0:
        # Big trouble -- matrix isn't invertible
        # Since we have no good conditioning code at present,
        # in this case we just fall back to using pari.
        H = W.echelon_form(algorithm='pari')
    elif 2*g > CUTOFF:
        # Unlikely that g will be large on even slightly random input
        # if it is, we fallback to the traditional algorithm.
        # A nasty example is A = n*random_matrix(ZZ,m), where
        # this algorithm gets killed.  This is not random input though.
        f = W.gcd()
        g = g / (f**W.nrows())
        if 2*g <= CUTOFF:
            verbose("Found common factor of %s -- dividing out; get new g = %s"%(f,g))
            W0 = (W/f).change_ring(ZZ)
            H = W0._hnf_mod(2*g)
            H *= f
        else:
            verbose("Falling back to PARI HNF since input matrix is ill conditioned for p-adic hnf algorithm.")
            # We need more clever preconditioning?
            # It is important to *not* just do the submatrix, since
            # the whole rest of the algorithm will likely be very slow in
            # weird cases where the det is large.
            # E.g., matrix all of whose rows but 1 are multiplied by some
            # fixed scalar n.
            raise NotImplementedError("fallback to PARI!")
            #H = W.hermite_form(algorithm='pari')
    else:
        H = W._hnf_mod(2*g)
        
    x = add_column(W, H, b.stack(matrix(1,1,[k*A[m-2,m-1] + l*A[m-1,m-1]])), proof)
    Hprime = H.augment(x)
    pivots = pivots_of_hnf_matrix(Hprime)
    
    Hprime, pivots = add_row(Hprime, A.matrix_from_rows([m-2]), pivots, include_zero_rows=False)
    Hprime, pivots = add_row(Hprime, A.matrix_from_rows([m-1]), pivots, include_zero_rows=False)
    H = Hprime.matrix_from_rows(range(m))
    return H

def interleave_matrices(A, B, cols1, cols2):
    """
    INPUT:
        A, B -- matrices with the same number of rows
        cols1, cols2 -- disjoint lists of integers
    OUTPUT:
        construct a new matrix C by sticking the columns
        of A at the positions specified by cols1 and the
        columns of B at the positions specified by cols2.

    EXAMPLES:
        sage: A = matrix(ZZ, 2, [1,2,3,4]); B = matrix(ZZ, 2, [-1,5,2,3])
        sage: A
        [1 2]
        [3 4]
        sage: B
        [-1  5]
        [ 2  3]
        sage: import sage.matrix.matrix_integer_dense_hnf as hnf
        sage: hnf.interleave_matrices(A, B, [1,3], [0,2])
        [-1  1  5  2]
        [ 2  3  3  4]
    """
    D = A.augment(B)
    w = cols1 + cols2
    v = [w.index(i) for i in range(len(cols1) + len(cols2))]
    return D.matrix_from_columns(v)

def probable_pivot_rows(A):
    """
    Return rows of A that are very likely to be pivots.

    This really finds the pivots of A modulo a random prime. 
    
    INPUT:
        A -- a matrix
    OUTPUT:
        a tuple of integers

    EXAMPLES:
        sage: import sage.matrix.matrix_integer_dense_hnf as matrix_integer_dense_hnf
        sage: a = matrix(ZZ,3,[0, -1, -1, 0, -20, 1, 0, 1, 2])
        sage: a
        [  0  -1  -1]
        [  0 -20   1]
        [  0   1   2]
        sage: matrix_integer_dense_hnf.probable_pivot_rows(a)
        (0, 1)
    """
    return probable_pivot_columns(A.transpose())

def probable_pivot_columns(A):
    """
    INPUT:
        A -- a matrix
    OUTPUT:
        a tuple of integers

    EXAMPLES:
        sage: import sage.matrix.matrix_integer_dense_hnf as matrix_integer_dense_hnf
        sage: a = matrix(ZZ,3,[0, -1, -1, 0, -20, 1, 0, 1, 2])
        sage: a
        [  0  -1  -1]
        [  0 -20   1]
        [  0   1   2]
        sage: matrix_integer_dense_hnf.probable_pivot_columns(a)
        (1, 2)
    """
    p = ZZ.random_element(10007, 46000).next_prime()
    pivots = A._reduce(p).pivots()
    return pivots

def ones(H, pivots):
    """
    Find all 1 pivot columns of the matrix H in Hermite form, along
    with the corresponding rows, and also the non 1 pivot columns and
    non-pivot rows.  Here a 1 pivot column is a pivot column so that
    the leading bottom entry is 1.

    INPUT:
        H -- matrix in Hermite form
        pivots -- list of integers (all pivot positions of H).

    OUTPUT:
        4-tuple of integer lists: onecol, onerow, non_oneol, non_onerow

    EXAMPLES:
        sage: H = matrix(ZZ, 3, 5, [1, 0, 0, 45, -36, 0, 1, 0, 131, -107, 0, 0, 0, 178, -145]); H
        [   1    0    0   45  -36]
        [   0    1    0  131 -107]
        [   0    0    0  178 -145]
        sage: import sage.matrix.matrix_integer_dense_hnf as matrix_integer_dense_hnf 
        sage: matrix_integer_dense_hnf.ones(H, [0,1,3])
        ([0, 1], [0, 1], [2], [2])
    """
    # Find the "onecol" pivot columns of H, i.e., the columns
    # that contain exactly one "1" entry and all other entries 0.
    onecol = []
    onerow = []
    i = 0
    for c in pivots:
        if H[i,c] == 1:
            onecol.append(c)
            onerow.append(i)
        i += 1
    onecol_set = set(onecol)
    non_onerow = [i for i in range(len(pivots)) if i not in onerow]
    non_onecol = [i for i in range(H.ncols()) if i not in onecol_set][:len(non_onerow)]
    return onecol, onerow, non_onecol, non_onerow

def extract_ones_data(H, pivots):
    """
    Compute ones data and corresponding submatrices of H.  This is
    used to optimized the add_row function.
    
    INPUT:
        H -- a matrix in HNF
        pivots -- list of all pivot column positions of H
        
    OUTPUT:
        C, D, E, onecol, onerow, non_onecol, non_onerow
        where onecol, onerow, non_onecol, non_onerow are as for
        the ones function, and C, D, E are matrices:
            C -- submatrix of all non-onecol columns and onecol rows
            D -- all non-onecol columns and other rows
            E -- inverse of D
        If D isn't invertible or there are 0 or more than 2 non onecols,
        then C, D, and E are set to None.

    EXAMPLES:
        sage: H = matrix(ZZ, 3, 4, [1, 0, 0, 7, 0, 1, 5, 2, 0, 0, 6, 6])
        sage: import sage.matrix.matrix_integer_dense_hnf as matrix_integer_dense_hnf 
        sage: matrix_integer_dense_hnf.extract_ones_data(H, [0,1,2])
        (
        [0]                                      
        [5], [6], [1/6], [0, 1], [0, 1], [2], [2]
        )
    
    Here we get None's since the (2,2) position submatrix is not invertible.
        sage: H = matrix(ZZ, 3, 5, [1, 0, 0, 45, -36, 0, 1, 0, 131, -107, 0, 0, 0, 178, -145]); H
        [   1    0    0   45  -36]
        [   0    1    0  131 -107]
        [   0    0    0  178 -145]
        sage: import sage.matrix.matrix_integer_dense_hnf as matrix_integer_dense_hnf         
        sage: matrix_integer_dense_hnf.extract_ones_data(H, [0,1,3])
        (None, None, None, [0, 1], [0, 1], [2], [2])
    """
    onecol, onerow, non_onecol, non_onerow = ones(H, pivots)
    verbose('extract_ones -- got submatrix of size %s'%len(non_onecol))
    if len(non_onecol) in [1,2]:
        # Extract submatrix of all non-onecol columns and onecol rows
        C = H.matrix_from_rows_and_columns(onerow, non_onecol)
        # Extract submatrix of all non-onecol columns and other rows 
        D = H.matrix_from_rows_and_columns(non_onerow, non_onecol).transpose()
        tt = verbose("extract ones -- INVERT %s x %s"%(len(non_onerow), len(non_onecol)), level=1)
        try:
            E = D**(-1)
        except ZeroDivisionError:
            C = D = E = None
        verbose("done inverting", tt, level=1)
        return C, D, E, onecol, onerow, non_onecol, non_onerow 
    else:
        return None, None, None, onecol, onerow, non_onecol, non_onerow

def is_in_hnf_form(H, pivots):
    """
    Return True precisely if the matrix H is in Hermite normal form
    with given pivot columns.

    INPUT:
        H -- matrix
        pivots -- sorted list of integers
        
    OUTPUT:
        bool -- True or False

    EXAMPLES:
        sage: a = matrix(ZZ,3,5,[-2, -6, -3, -17, -1, 2, -1, -1, -2, -1, -2, -2, -6, 9, 2])
        sage: import sage.matrix.matrix_integer_dense_hnf as matrix_integer_dense_hnf         
        sage: matrix_integer_dense_hnf.is_in_hnf_form(a,range(3))
        False
        sage: e = a.hermite_form(); p = a.pivots()
        sage: matrix_integer_dense_hnf.is_in_hnf_form(e, p)
        True
    """
    tt = verbose('testing if matrix is in HNF')
    r = 0
    pivots_set = set(pivots)
    for j in xrange(H.ncols()):
        if j in pivots_set:
            for i in xrange(r+1, H.nrows()):
                if H[i,j]:
                    verbose('not HNF because nonzeros below pivot position',tt)
                    return False
            for i in xrange(r):
                if H[i,j] < 0 or H[i,j] >= H[r,j]:
                    verbose('not HNF because negative or too big above pivot position',tt)                    
                    return False
            r += 1
        else:
            for i in xrange(r,H.nrows()):
                if H[i,j]:
                    verbose('not HNF nonzero in wrong place in nonpivot column',tt)                    
                    return False
    verbose('done verifying in HNF -- yes', tt)
    return True

def probable_hnf(A, include_zero_rows, proof):
    """
    Return the HNF of A or raise an exception if something involving
    the randomized nature of the algorithm goes wrong along the way.
    Calling this function again a few times should result it in it
    working, at least if proof=True.
        
    INPUT:
        A -- a matrix
        include_zero_rows -- bool
        proof -- bool
        
    OUTPUT:
        the Hermite normal form of A.
        cols -- pivot columns

    EXAMPLES:
        sage: a = matrix(ZZ,4,3,[-1, -1, -1, -20, 4, 1, -1, 1, 2,1,2,3])
        sage: import sage.matrix.matrix_integer_dense_hnf as matrix_integer_dense_hnf         
        sage: matrix_integer_dense_hnf.probable_hnf(a, True, True)
        (
        [1 0 0]           
        [0 1 0]           
        [0 0 1]           
        [0 0 0], [0, 1, 2]
        )
        sage: matrix_integer_dense_hnf.probable_hnf(a, False, True)
        (
        [1 0 0]           
        [0 1 0]           
        [0 0 1], [0, 1, 2]
        )
        sage: matrix_integer_dense_hnf.probable_hnf(a, False, False)
        (
        [1 0 0]           
        [0 1 0]           
        [0 0 1], [0, 1, 2]
        )
    """
    # Find left-most full rank submatrix by working modulo a prime
    rows = list(probable_pivot_rows(A))
    B    = A.matrix_from_rows(rows)
    cols = list(probable_pivot_columns(B))
    C   = B.matrix_from_columns(cols)
    # Now C is a submatrix of A that has full rank and is square.

    # We compute the HNF of C, which is a square nonsingular matrix.
    try:
        H = hnf_square(C, proof=proof)
    except NotImplementedError:
        # raise
        # this signals that we must fallback to pari
        verbose("generic random modular HNF algorithm failed -- we fall back to PARI")
        H = A.hermite_form(algorithm='pari', include_zero_rows=include_zero_rows, proof=proof)
        return H, H.pivots()

    # The transformation matrix to HNF is the unique
    # matrix U such that U * C = H, i.e., U = H*C^(-1).

    if len(cols) < B.ncols():
        # We compute the HNF of B by multiplying the matrix D
        # got from the columns not in C by U:
        # We want to compute X = U*D.  But U = H*C^(-1),
        # so X = U*D = H*C^(-1)*D.
        # So C*H^(-1)*X = D

        # find y s.t C*y = D
        #   H^(-1)*X = y ===> X = H*y
        #   
        cols_set = set(cols)
        cols2 = [i for i in range(B.ncols()) if not i in cols_set]
        D = B.matrix_from_columns(cols2)
        Y = C.solve_right(D)
        H2 = H*Y
        H2 = H2.change_ring(ZZ)
        
        # The HNF of B is got by assembling together
        # the matrices H and H2.
        H = interleave_matrices(H, H2, cols, cols2)

    pivots = pivots_of_hnf_matrix(H)

    # Now H is the HNF of the matrix B.
    # Finally we add all remaining rows of A to H using
    # the add_row function.

    C, D, E, onecol, onerow, non_onecol, non_onerow = extract_ones_data(H, cols)
    if not proof and len(non_onecol) == 0:
        # Identity matrix -- done
        verbose("hnf -- got identity matrix -- early abort (0)")
        if include_zero_rows: H = pad_zeros(H, A.nrows())
        return H, pivots
        
    rows_set = set(rows)
    for i in range(A.nrows()):
        if not i in rows_set:
            v = A.matrix_from_rows([i])
            if v == 0: continue
            if E is None:
                H, pivots = add_row(H, v, pivots, include_zero_rows=False)
                C, D, E, onecol, onerow, non_onecol, non_onerow = extract_ones_data(H, pivots)
                if not proof and len(non_onecol) == 0:
                    # Identity matrix -- done
                    verbose("hnf -- got identity matrix -- early abort (1)")
                    if include_zero_rows: H = pad_zeros(H, A.nrows())                    
                    return H, pivots
            else:
                z = A.matrix_from_rows_and_columns([i], non_onecol)                
                w = A.matrix_from_rows_and_columns([i], onecol)
                tt = verbose("checking denom (%s x %s)"%(D.nrows(), D.ncols()))
                Y = (z - w*C).transpose()
                k = E*Y
                verbose("done checking denom",tt)
                if k.denominator() != 1:
                    H, pivots = add_row(H, v, pivots, include_zero_rows=False)
                    D = H.matrix_from_rows_and_columns(non_onerow, non_onecol).transpose()
                nn = ones(H, pivots)
                if not proof and len(nn[2]) == 0:
                    verbose("hnf -- got identity matrix -- early abort (2)")
                    if include_zero_rows: H = pad_zeros(H, A.nrows())
                    return H, pivots

    if include_zero_rows: H = pad_zeros(H, A.nrows())
    return H, pivots

def pad_zeros(A, nrows):
    """
    Add zeros to the bottom of A so that the
    resulting matrix has nrows.

    INPUT:
        A -- a matrix
        nrows -- an integer that is at least as big as the number of rows of A.

    OUTPUT:
        a matrix with nrows rows.

    EXAMPLES:
        sage: import sage.matrix.matrix_integer_dense_hnf as matrix_integer_dense_hnf 
        sage: a = matrix(ZZ, 2, 4, [1, 0, 0, 7, 0, 1, 5, 2])
        sage: matrix_integer_dense_hnf.pad_zeros(a, 4)
        [1 0 0 7]
        [0 1 5 2]
        [0 0 0 0]
        [0 0 0 0]
        sage: matrix_integer_dense_hnf.pad_zeros(a, 2)
        [1 0 0 7]
        [0 1 5 2]
    """
    nz = nrows - A.nrows()
    if nz == 0:
        return A
    if nz < 0:
        return A.matrix_from_rows(range(nrows))
    return A.stack(matrix(ZZ, nz, A.ncols()))


def hnf(A, include_zero_rows=True, proof=True):
    """
    Return the Hermite Normal Form of a general integer matrix A,
    along with the pivot columns. 
    
    INPUT:
        A -- an n x m matrix A over the integers.
        include_zero_rows -- bool (default: True) whether or not to
                             include zero rows in the output matrix
        proof -- whether or not to prove the result correct.  
        
    OUTPUT:
        matrix -- the Hermite normal form of A
        pivots -- the pivot column positions of A

    EXAMPLES:
        sage: import sage.matrix.matrix_integer_dense_hnf as matrix_integer_dense_hnf 
        sage: a = matrix(ZZ,3,5,[-2, -6, -3, -17, -1, 2, -1, -1, -2, -1, -2, -2, -6, 9, 2])
        sage: matrix_integer_dense_hnf.hnf(a)
        (
        [   2    0   26  -75  -10]           
        [   0    1   27  -73   -9]           
        [   0    0   37 -106  -13], [0, 1, 2]
        )
        sage: matrix_integer_dense_hnf.hnf(a.transpose())
        (
        [1 0 0]           
        [0 1 0]           
        [0 0 1]           
        [0 0 0]           
        [0 0 0], [0, 1, 2]
        )
        sage: matrix_integer_dense_hnf.hnf(a.transpose(), include_zero_rows=False)
        (
        [1 0 0]           
        [0 1 0]           
        [0 0 1], [0, 1, 2]
        )
    """
    if A.nrows() <= 1:
        np = A.nonzero_positions()
        if len(np) == 0:
            pivots = []
            if not include_zero_rows:
                A = A.new_matrix(0)  # 0 rows
        else:
            i,j = np[0]
            if A[i,j] < 0:
                A = -A
            pivots = [j]
        return A, pivots
    
    if proof == False:
        H, pivots = probable_hnf(A, include_zero_rows = include_zero_rows, proof=False)
        if not include_zero_rows and len(pivots) > H.nrows():
            return H.matrix_from_rows(range(len(pivots))), pivots
    
    while True:
        try:
            H, pivots = probable_hnf(A, include_zero_rows = include_zero_rows, proof=True)
        except (AssertionError, ZeroDivisionError, TypeError):
            raise
            #verbose("Assertion occurred when computing HNF; guessed pivot columns likely wrong.")
            #continue
        else:
            if is_in_hnf_form(H, pivots):
                if not include_zero_rows and len(pivots) > H.nrows():
                    H = H.matrix_from_rows(range(len(pivots)))
                return H, pivots
            verbose("After attempt the return matrix is not in HNF form since pivots must have been wrong.  We try again.")
    
def hnf_with_transformation(A, proof=True):
    """
    Compute the HNF H of A along with a transformation matrix U
    such that U*A = H.  Also return the pivots of H.
    
    INPUT:
        A -- an n x m matrix A over the integers.
        proof -- whether or not to prove the result correct.  
        
    OUTPUT:
        matrix -- the Hermite normal form H of A
        U -- a unimodular matrix such that U * A = H
        pivots -- the pivot column positions of A

    EXAMPLES:
        sage: import sage.matrix.matrix_integer_dense_hnf as matrix_integer_dense_hnf     
        sage: A = matrix(ZZ, 2, [1, -5, -10, 1, 3, 197]); A
        [  1  -5 -10]
        [  1   3 197]
        sage: H, U, pivots = matrix_integer_dense_hnf.hnf_with_transformation(A)
        sage: H
        [  1   3 197]
        [  0   8 207]
        sage: U
        [ 0  1]
        [-1  1]
        sage: U*A
        [  1   3 197]
        [  0   8 207]
    """
    # All we do is augment the input matrix with the identity matrix of the appropriate rank on the right.
    C = A.augment(identity_matrix(ZZ, A.nrows()))
    H, pivots = hnf(C, include_zero_rows=True, proof=proof)
    U = H.matrix_from_columns(range(A.ncols(), H.ncols()))
    H2 = H.matrix_from_columns(range(A.ncols()))
    return H2, U, pivots

def hnf_with_transformation_tests(n=10, m=5, trials=10):
    """
    Use this to randomly test that hnf with transformation matrix
    is working.
    
    EXAMPLES:
        sage: from sage.matrix.matrix_integer_dense_hnf import hnf_with_transformation_tests
        sage: hnf_with_transformation_tests(n=15,m=10, trials=10)
        0 1 2 3 4 5 6 7 8 9
    """
    import sys
    for i in range(trials):
        print i,
        sys.stdout.flush()
        a = random_matrix(ZZ, n, m)
        w = hnf_with_transformation(a)
        assert w[0] == w[1]*a
        w = hnf_with_transformation(a, proof=False)
        assert w[0] == w[1]*a
    

#################################################################################################
# Code for testing and benchmarking
#################################################################################################   
def benchmark_hnf(nrange, bits=4):
    """
    Run benchmark program.
    
    EXAMPLES:
        sage: import sage.matrix.matrix_integer_dense_hnf as hnf
        sage: hnf.benchmark_hnf([50,100],32)   
        ('sage', 50, 32, ...),
        ('sage', 100, 32, ...),
    """
    b = 2**bits
    for n in nrange:
        a = random_matrix(ZZ, n, x=-b,y=b)
        t = cputime()
        h,_ = hnf(a, proof=False)
        tm = cputime(t)
        print '%s,'%(('sage', n, bits, tm),)
   
def benchmark_magma_hnf(nrange, bits=4):
    """
    EXAMPLES:
        sage: import sage.matrix.matrix_integer_dense_hnf as hnf
        sage: hnf.benchmark_magma_hnf([50,100],32)     # optional - magma
        ('magma', 50, 32, ...),
        ('magma', 100, 32, ...),
    """
    from sage.interfaces.all import magma
    b = 2**bits
    for n in nrange:
        a = magma('MatrixAlgebra(IntegerRing(),%s)![Random(%s,%s) : i in [1..%s]]'%(n,-b,b,n**2))
        t = magma.cputime()
        h = a.EchelonForm()
        tm = magma.cputime(t)
        print '%s,'%(('magma', n, bits, tm),)


global sanity
def sanity_checks(times=50, n=8, m=5, proof=True, stabilize=2, check_using_magma = True):
    """
    Run random sanity checks on the modular p-adic HNF with tall and wide matrices
    both dense and sparse.

    INPUT:
        times -- number of times to randomly try matrices with each shape
        n -- number of rows
        m -- number of columns
        proof -- test with proof true
        stabilize -- parameter to pass to hnf algorithm when proof is False
        check_using_magma -- if True use Magma instead of PARI to
                             check correctness of computed HNF's.
                             Since PARI's HNF is buggy and slow (as of
                             2008-02-16 non-pivot entries sometimes
                             aren't normalized to be nonnegative) the
                             default is Magma.
    
    EXAMPLES:
        sage: import sage.matrix.matrix_integer_dense_hnf as matrix_integer_dense_hnf     
        sage: matrix_integer_dense_hnf.sanity_checks(times=5, check_using_magma=False)
        small 8 x 5
        0 1 2 3 4  (done)
        big 8 x 5
        0 1 2 3 4  (done)
        small 5 x 8
        0 1 2 3 4  (done)
        big 5 x 8
        0 1 2 3 4  (done)
        sparse 8 x 5
        0 1 2 3 4  (done)
        sparse 5 x 8
        0 1 2 3 4  (done)
        ill conditioned -- 1000*A -- 8 x 5
        0 1 2 3 4  (done)
        ill conditioned -- 1000*A but one row -- 8 x 5
        0 1 2 3 4  (done)
    """
    import sys
    def __do_check(v):
        """
        This is used internally by the sanity check code.
        """
        for i,a in enumerate(v):
            global sanity
            sanity = a
            print i,
            sys.stdout.flush()
            if check_using_magma:
                if magma(hnf(a)[0]) != magma(a).EchelonForm():
                    print "bug computing hnf of a matrix"
                    print 'a = matrix(ZZ, %s, %s, %s)'%(a.nrows(), a.ncols(), a.list())
                    return
            else:
                if hnf(a)[0] != a.echelon_form('pari'):
                    print "bug computing hnf of a matrix"
                    print 'a = matrix(ZZ, %s, %s, %s)'%(a.nrows(), a.ncols(), a.list())
                    return
        print " (done)"

    print "small %s x %s"%(n,m)
    __do_check([random_matrix(ZZ, n, m, x=-1,y=1) for _ in range(times)])
    print "big %s x %s"%(n,m)
    __do_check([random_matrix(ZZ, n, m, x=-2^32,y=2^32) for _ in range(times)])

    print "small %s x %s"%(m,n)
    __do_check([random_matrix(ZZ, m, n, x=-1,y=1) for _ in range(times)])
    print "big %s x %s"%(m,n)
    __do_check([random_matrix(ZZ, m, n, x=-2^32,y=2^32) for _ in range(times)])    

    print "sparse %s x %s"%(n,m)    
    __do_check([random_matrix(ZZ, n, m, density=0.1) for _ in range(times)])
    print "sparse %s x %s"%(m,n)    
    __do_check([random_matrix(ZZ, m, n, density=0.1) for _ in range(times)])

    print "ill conditioned -- 1000*A -- %s x %s"%(n,m)
    __do_check([1000*random_matrix(ZZ, n, m, x=-1,y=1) for _ in range(times)])

    print "ill conditioned -- 1000*A but one row -- %s x %s"%(n,m)
    v = []
    for _ in range(times):
        a = 1000*random_matrix(ZZ, n, m, x=-1,y=1)
        a[a.nrows()-1] = a[a.nrows()-1]/1000
        v.append(a)
    __do_check(v)