""" The method QuadraticForm::basis_of_short_vectors currently doesn't return 
    a basis all the time, the following can be used to correct this defect.

    An example of this problem is the E8 lattice:
        Q = QuadraticForm( matrix( [[2,0,0,0,0,0,0,1],
                                    [0,2,1,1,1,1,1,1],
                                    [0,1,2,1,1,1,1,1],
                                    [0,1,1,2,1,1,1,1],
                                    [0,1,1,1,2,1,1,1],
                                    [0,1,1,1,1,2,1,1],
                                    [0,1,1,1,1,1,2,0],
                                    [1,1,1,1,1,1,0,2]] ))
        B = Q.basis_of_short_vectors()
        matrix(B).det()
    the result of this will be -2, which indicates this is not a basis.
    We can use the following code to correct this by:
        B = fix_basis(Q, matrix(B).transpose()).columns()
        # B now is a basis
        matrix(B).det()
    The determinant will be -1.
    We remark that the code may not return an optimally short solution.
"""

def find_missed_1(M,L):
    """
        Look for vector v in L (shortest) such that v is not in the image of M
        and there is a substitution of v for a column of M such that the image
        of M' would is strictly larger than that of M
    """
    M2 = M.inverse()
    for Li in L:
        for v in Li:
            v2 = M2*vector(v)
            div = lcm([ x.denominator() for x in v2])
            if div != 1:
                for x in (M2*vector(v)):
                    if x.denominator() == div and x.numerator().abs() == 1:
                        return v
    return 0

def find_missed_2(M,L):
    """
        As above, except only require the image of M' to have smaller index
        than that of M
    """
    M2 = M.inverse()
    for Li in L:
        for v in Li:
            v2 = M2*vector(v)
            div = lcm([ x.denominator() for x in v2])
            if div != 1:
                for x in (M2*vector(v)):
                    if x.denominator() == div and x.numerator().abs() < div:
                        return v
    return 0

def fix_basis(Q,M):
    """
        Iteratively attempt to convert the matrix M into a change of basis
        matrix of Q by substituting small vectors
    """
    new = M
    depth = 1
    while new.det().abs() != 1:
        # Outer loop simply controls generating too many vectors to test
        depth = depth + 1
        L = Q.short_vector_list_up_to_length(depth)
        while new.det().abs() != 1:
            # Inner loop iterates over possible substitutions
            v = find_missed_1(new,L)
            C = new.columns();
            if v == 0:
                # We didn't find a clean substitute, look for a non-clean one
                v = find_missed_2(new,L)
                if v == 0:
                    break
                div = lcm([ x.denominator() for x in (new.inverse()*v)])
                v2 = new.inverse()*v*div
                for i in range(Q.dim()):
                    if v2[i].abs() < div:
                        C[i] = v
                        new = (matrix(C)).transpose()
                        break
            else:
                div = lcm([ x.denominator() for x in (new.inverse()*v)])
                v2 = new.inverse()*v*div
                for i in range(Q.dim()):
                    if v2[i].abs() == 1:
                        C[i] = v
                        new = (matrix(C)).transpose()
                        break
    return new