"""
F5/Matrix (first attempt)

AUTHOR:
    -- 20081015 Martin Albrecht

EXAMPLE:
    sage: P.<a,b,c,d,e,h> = PolynomialRing(GF(127))
    sage: I = P * [29*b^2 - 12*d^2 + 21*a*e + 43*e^2 - 43*h^2, \
                   -20*a*b - 15*c^2 + 30*d^2 - 2*d*e - 24*h^2, \
                   24*b^2 + 13*a*d - 49*a*e - 17*a*h - 21*h^2, \
                   29*c*d + 50*d^2 - 26*c*e + 30*a*h,  \
                   30*b*d - 37*c*d - 63*c*e - 3*e*h - 29*h^2]
    sage: f5 = F5Matrix()
    sage: f5(I,5)
"""

from itertools import ifilter

from sage.rings.polynomial.multi_polynomial_ideal import MPolynomialIdeal

class F5Matrix:
    """
    A toy implementation of F5/Matrix, the matrix version of F5.
    """
    def __init__(self):
        self.M = [[]] # macaulays matrices
        self.L = [[]] # list of labeled polynomials

    def __call__(self, F, dmax):
        if isinstance(F, MPolynomialIdeal):
            F = F.gens()
        else:
            F = list(F)

        if not all(f.is_homogeneous() for f in F):
            F = Ideal(F).homogenize()
            F = F.gens()

        return self.basis(sorted(F,reverse=False), dmax)

    def basis(self, F, D):
        r"""
        INPUT: 
        
        - ``F`` - `f_1, ..., f_m` homogeneous polynomials of degrees `d_1 <= d_2 <= ... d_m`
        - ``D`` - maximal degree (``>0``)

        OUTPUT: 
        
            A ``D``-Groebner basis for the ideal spanned by `f_1, ... f_m`
        """
        gauss_elimination = self.gauss_elimination
        add_signature = self.add_signature
        self.__init__()

        M = self.M
        L = self.L

        R = F[0].parent()
        variables = sorted(R.gens())

        for d in range(1, D+1):
            if get_verbose() >= 0:
                print "d %2d"%d,
            L.append([])
            M.append([])

            for i,f in enumerate(F):
                di = f.degree()

                if di == d:
                    add_signature( MSignature(R(1), i, len(M[d])), d)
                    M[d].append(f)
                    continue
                    
                for (mult, m, r) in ifilter(lambda elem: elem[1] == i, L[d-1]):
                    assert(m == i)

                    for x in variables:
                        # Avoid to consider xy and yx as multiplier.
                        if mult != 1 and x < max(mult.variables()):
                           continue
                        
                        # For all j < m, if we have a row labeled (t, f_j)
                        # in the matrix `M_{D-dm}(f_1,...f_m-1)` that has
                        # leading t' then you can remove the row (t', f_m)
                        # in `M_{D}(f_1,...,f_m)`

                        found = False
                        for (t, j , rbar) in L[d-di]:
                            if j>=i:
                                continue
                            if M[d-di][rbar].lm() == x*mult:
                                found = True
                                break
                        if not found:
                            add_signature( MSignature(x*mult, i, len(M[d])), d )
                            M[d].append( x*M[d-1][r] )
                        

            Mtilde = gauss_elimination(M[d])
            Lbar = []
            Mbar = []

            nz = 0
            for i,f in enumerate(Mtilde):
                mult, idx, row = L[d][i]
                if f == 0:
                    verbose("reduction to zero %s %d %d"%(mult, idx, row),level=0)
                    pass
                else:
                    Lbar.append( (mult, idx, nz) )
                    Mbar.append( f )
                    nz += 1

            M[d] = Mbar
            L[d] = Lbar

        return Ideal([f for d in range(len(self.M)) for f in self.M[d]]).interreduced_basis()
        
    def add_signature(self, s, d):
        self.L[d].append(s)

    def gauss_elimination(self, F):
        check = True
        if len(F) == 0:
            print
            return F
        A,v = mq.MPolynomialSystem(F).coefficient_matrix()
        if check:
            E = A.echelon_form()
        print "%4d x %4d, %4d, %4d"%(A.nrows(), A.ncols(), A.rank(), A.nrows()-A.rank())
        nrows, ncols = A.nrows(), A.ncols()

        if False:
            start_row = 0
            for c in range(ncols):
                for r in range(start_row, nrows):
                    if A[r, c]:
                        a_inverse = ~A[r,c]
                        A.rescale_row(r, a_inverse, c)
                        A.swap_rows(r, start_row)
                        for i in range(0, nrows):
                            if i == start_row:
                                continue
                            if A[i,c]:
                                minus_b = -A[i, c]
                                A.add_multiple_of_row(i, start_row, minus_b, c)
                        start_row = start_row + 1
                        break
            return (A*v).list()
        for c in xrange(ncols):
            for r in xrange(nrows):
                if A[r,c] != 0:
                    if any(A[r,i] for i in xrange(c)):
                       continue

                    a_inverse = ~A[r,c]
                    A.rescale_row(r, a_inverse, c)
                    for i in xrange(r+1,nrows):
                        if A[i,c] != 0:
                            #if any(A[i,_] for _ in xrange(c)):
                            #    continue
                            minus_b = -A[i,c]
                            A.add_multiple_of_row(i, r, minus_b, c)
                    break
        if check:
            Al = set([f.lm() for f in (A*v).list()])
            El = set([f.lm() for f in (E*v).list()])
            assert(Al == El)
        return (A*v).list()

from UserList import UserList
class MSignature(UserList):
    def __init__(self, multiplier, polynomial, row):
        UserList.__init__(self, [multiplier, polynomial, row])

f5 = F5Matrix()