# -*- Mode: Python -*-
# # vi:si:et:sw=4:sts=4:ts=4
#

"""
F4

AUTHOR: Martin Albrecht <malb@informatik.uni-bremen.de> 
"""

from sage.rings.all import *
from sage.misc.misc import exists
from sage.rings.ideal import is_Ideal
from sage.matrix.matrix_modn_sparse import Matrix_modn_sparse
from sage.matrix.matrix_space import MatrixSpace

def CoeffMatrix_modint(F):
    """
    See MQ.matrix()
    """

    if is_Ideal(F):
        F = F.gens()
    if isinstance(F,set):
        F = list(F)
        
    R = F[0].parent()
    k = R.base_ring()

    m = set([mon for f in F for mon in f.monomials()])
    m = tuple(sorted(m,reverse=True))

    #construct dictionary for fast lookups
    v = dict( map(lambda x,y:(x,y), m , range(0,len(m)) ) )

    MS = MatrixSpace(R.base_ring(), len(F), len(v), sparse=True)
    A  = CoeffMatrix_modint_impl(MS,{},False,False)
        
    for x in range( 0 , len(F) ):
        poly = F[x]
        for y in poly.monomials():
            A[ x , v[y] ] = poly.monomial_coefficient(y)
    print "% 4d x % 4d, % 4d, % 4d"%(A.nrows(), A.ncols(), A.rank(), A.nrows()-A.rank())
    return ( A , m )

class CoeffMatrix_modint_impl(Matrix_modn_sparse):
    """
    Sparse Matrix with some extra functions suitable
    for a coefficient matrix.
    """
    
    def rows_dict(self):
        """
        Returns a dictionary of whichs keys are row numbers and whichs
        values are ordered lists of column numbers. Those row/column
        number pairs represent non-zero entries.
        """
        try:
            return self.__md
        except AttributeError:
            d = self._dict()
            md = {}
            for point in d.iterkeys():
                if not md.has_key(point[0]):
                    md[point[0]] = [point[1]]
                else:
                    md[point[0]].append(point[1])
            for values in md.itervalues():
                values.sort()
            self.__md = md
            return md

    def polynomial(self,v,i):

        R = v[0].parent()
        poly_dict = dict()
        rd = self.rows_dict()
        f = R(0)
        try:
            for col in rd[self.nrows()-i-1]:
                f += self[self.nrows()-i-1, col] * v[col]
        except KeyError:
            pass

        return f

    def __mul__(self,right):
        """
        The product of a coefficent matrix with a monomial vector
        is a list of poylnomials.
        """
        if not isinstance(right, tuple):
            raise NotImplementedError

        R = right[0].parent()

        rd = self.rows_dict()
        rows = sorted(rd.iterkeys())
        nrows = self.nrows()
        return [self.polynomial(right,nrows-i-1) for i in rows]


class F4_orig:
    """
    Original F4 as described by Faugere
    """
    def __init__(self):
        pass

    def __call__(self, F, sel=None):
        if is_Ideal(F):
            F = F.gens()

        self.ring = F[0].parent()
        self.rr_bases = []
        G = list(F)
        F0p = F
        i = 0
        P = set([self.pair(f,g) for f in G for g in G if f<g ] )

        if sel==None:
            sel = self.normal_strategy

        while P != set():
            i += 1
            Pd, d = sel(P)
            print "% 2d"%(d,), 
            P = P.difference(Pd)
            Ld = set(self.left(Pd)).union(set(self.right(Pd)))
            Fdp = self.reduction(Ld,G)
            for h in Fdp:
                P = P.union(set([self.pair(h,g) for g in G ]))
                G.append(h)
            sys.stdout.flush()
        return G

    def reduction(self,L,G):
        F = self.symbolic_preprocessing(L,G)
        Ft = self.row_echelon(F)
        LMF = LM(F)
        Ftp = set([f for f in Ft if f.lm() not in LMF])
        return list(Ftp)

    def symbolic_preprocessing(self,L,G):
        """
        """ 
        G = G
        F = set([t*f for (t,f) in L ])
        Done = LM(F)
        M = set([m  for f in F for m in f.monomials()])
        R = self.ring
        while M != Done:
            m = M.difference(Done).pop()
            Done.add(m)
            t,g = self.ring.monomial_reduce(m,G)
            if t!=R(0): F.add(t*g)
            M = set([m  for f in F for m in f.monomials()])
        return F

    def pair(self,f,g):
        lcm = self.ring.monomial_lcm(f.lm(),g.lm())
        # it seems better speed-wise to calculate those on the fly
        #tf = LCMdLM(lcm,f.lm())
        #tg = LCMdLM(lcm,g.lm())
        return (lcm,f,g)
    
    def left(self,p):
        if isinstance(p,(list,set,tuple)):
            s = set()
            for f in p:
                s.add((self.ring.monomial_quotient(f[0],f[1].lm()),f[1]))
            return s
        else:
            return (self.ring.monomial_quotient(p[0],p[1].lm()),p[1])

    def right(self,p):
        if isinstance(p,(list,set,tuple)):
            s = set()
            for f in p:
                s.add((self.ring.monomial_quotient(f[0],f[2].lm()),f[2]))
            return s
        else:
            return (self.ring.monomial_quotient(p[0],p[2].lm()),p[2])

    def row_echelon(self, F):
        """
        """
        A,v=CoeffMatrix_modint(F)
        if self.protocol:
            A.visualize_structure(maxsize=None)

        A.echelonize()

        if self.protocol:
            A.visualize_structure(maxsize=None)

        F = A*v
        return F

    # Selection Strategies

    def normal_strategy(self,P):
        """

        The normal selection strategy

        INPUT:
            P -- a list of critical pairs

        OUTPUT:
            a sublist of P
        
        """ 
        d = min(set([ lcm.total_degree() for (lcm,fi,fj) in P ]))
        return set([ (lcm,fi,fj) for (lcm,fi,fj) in P if lcm.total_degree()==d]), d

    # Update Strategies

    def update_pairsGF2(self,G,B,h):
        """
        Following Becker, 'Groebner Bases', Springer 1993 as suggested
        by Faugere in his F4 paper. Also works in the quotient ring.

        INPUT:
            G -- an intermediate Groebner basis
            B -- a list of critical pairs
            h -- a polynomial

        OUTPUT:
            an intermediate Groebner basis, a list of critical pairs

        """

        R = self.ring

        hlm = h.lm()

        G_new = list()


        # if G is a set then C only contains unique elements
        C = list([self.pair(h,g) for g in G]) # 1.86
        D = list() # only adding elements of C, thus unique


        # Criterion F & M 
        while C!=list():
            (lcmhg1,h,g1) = C.pop()

            # will be removed in next loop
            if R.monomial_pairwise_prime(hlm,g1.lm()):
                D.append((lcmhg1,h,g1))
                continue

            found = 0
            for c in C:
                if R.monomial_divides( c[0], lcmhg1 ):
                    found=1; break
            if found: continue

            found = 0
            for d in D:
                if R.monomial_divides( d[0], lcmhg1 ):
                    found=1; break
            if found: continue

            D.append((lcmhg1,h,g1))


        E = list() #only adding elements of D, thus unique


        # Buchberger criterion 1
        for (lcmhg,h,g) in D:
            # if LM(h) and LM(g) are not disjoint
            if not R.monomial_pairwise_prime(hlm,g.lm()):
                E.append((lcmhg,h,g))


        B_new = set()

        # Criterion B_k
        for (lcmg1g2,g1,g2) in B:
            if not self.ring.monomial_divides( hlm, lcmg1g2 ) or \
                   self.ring.monomial_lcm(g1.lm(), hlm) == lcmg1g2 or \
                   self.ring.monomial_lcm( hlm,g2.lm()) == lcmg1g2 :
                B_new.add((lcmg1g2,g1,g2))

        B_new = B_new.union(E)

        r=[]
        for gen in hlm.variables():
            # consider F -- always true
            # consider M
            for f in G:
                if R.monomial_divides( hlm, f.lm() ):
                    r.append( (gen,f) )
                    break
        r.append((R(1),h))
        G_new = self.reduction(r,[],[],True)[0]

        
        for g in G: 
            if not R.monomial_divides(hlm, g.lm()):
                G_new.append(g)
        G_new.append(h)

        return G_new,B_new

    def update_pairs(self,G,B,h):
        """
        Following Becker, 'Groebner Bases', Springer 1993 as suggested
        by Faugere in his F4 paper.

        INPUT:
            G -- an intermediate Groebner basis
            B -- a list of critical pairs
            h -- a polynomial

        OUTPUT:
            an intermediate Groebner basis, a list of critical pairs

        """

        R = self.ring

        # if G is a set then C only contains unique elements
        C = [self.pair(h,g) for g in G]
        D = list() # only adding elements of C, thus unique

        # Criterion M
        
        while C!=list():
            (lcmhg1,h,g1) = C.pop()
            
            lcm_divides = lambda lcmhg2: R.monomial_divides(  lcmhg2[0], lcmhg1 )

            # if LM(h) and LM(g_1) are disjoint
            if R.monomial_pairwise_prime(h.lm(),g.lm()) or \
               (\
                   not exists(C, lcm_divides )[0] \
                   and \
                   not exists(D, lcm_divides )[0]\
                ):
                D.append((lcmhg1,h,g1))

        E = list() #only adding elements of D, thus unique

        # Criterion F
        
        while D != list():
            (lcmhg,h,g) = D.pop()
            # if LM(h) and LM(g) are not disjoint
            if not R.monomial_pairwise_prime(h.lm(),g.lm()):
                E.append((lcmhg,h,g))

        B_new = set()

        # Criterion B_k

        while B != set():
            lcmg1g2,g1,g2 = B.pop()
            if not self.ring.monomial_divides( h.lm(), lcmg1g2 ) or \
                   self.ring.monomial_lcm(g1.lm(), h.lm()) == lcmg1g2 or \
                   self.ring.monomial_lcm( h.lm(),g2.lm()) == lcmg1g2 :
                B_new.add((lcmg1g2,g1,g2))

        B_new = B_new.union(E)

        G_new = list()

        while G != list():
            g = G.pop()
            if not R.monomial_divides( h.lm(), g.lm() ):
                G_new.append(g)

        G_new.append(h)

        return G_new,B_new

    def update_simple(self,G, P, h):
        """
        Adding all critical pairs

        INPUT:
            G -- an intermediate Groebner basis
            B -- a list of critical pairs
            h -- a polynomial

        OUTPUT:
            an intermediate Groebner basis, a list of critical pairs
        
        """
        return G+[h],P.union([self.pair(g,h) for g in G])


class F4(F4_orig):
    """
    The improved F4 as described in Faugere's paper.
    """

    def __call__(self, F, Sel=None, Update=None, protocol=False):
        """
        INPUT:
            F      -- a finite subset of R[x]
            Sel    -- selection strategy
            Update -- update pairs to select critical pairs to compute

        OUTPUT:
            G -- a Groebner Basis for F
        """
        self.protocol = protocol
        
        if is_Ideal(F):
            F = F.gens()

        # pretty looking code
        Left = self.left
        Right = self.right
        Reduction = self.reduction
        first = self.first
        if Sel==None: Sel = self.normal_strategy
        if Update==None: Update = self.update_pairs
        
        self.ring = F[0].parent()
        self.ring._singular_().set_ring()
        self.term_order = self.ring.term_order()

        # We maintain a list of dictionaries which contain f.lm() => f
        # maps for the sets $F_j^~$ to allow O(1) lookups for this code:
        #"$F_j^~$ is the row echelon form of F_j w.r.t. < there exists a
        # (unique) $p \in F_j^~ such that LM(p) = LM(u*f)"
        self.Ftd = [[]]
        
        F = list(F) #
        Fd = dict()

        G = list()
        P = set()
        i = 0

        while F != list():
            f = first(F)
            F.remove(f)
            G,P = Update(G,P,f)
        sys.stdout.flush()

        while P != set(): 
            i += 1

            Pd, d = Sel(P)
            print "%2d"%(d),
            if self.protocol:
                print "P",sorted(P)
                print "P%d"%d,sorted(Pd)
                print "G",sorted(G)
            P = P.difference(Pd)
            Ld = Left(Pd).union( Right(Pd) )
            if self.protocol:
                print "L%d"%d,sorted(Ld)
            Fdp,Fd[i] = Reduction(Ld,G,Fd)
            for h in Fdp:
                G,P = Update(G,P,h)
        return G

    def reduction(self,L,G,Fset,no_update=False):
        """
        INPUT:
            L -- a finite subset of M x R[x]
            G -- a finite subset of R[x]
            F -- (F_k)k=1,\dots,(d-1), where F_k is finite subset of R[x]

        OUTPUT:
            F~+,F 
        """
        F = self.symbolic_preprocessing(L,G,Fset)
        #print max([len(f.dict()) for f in F])
        if self.protocol:
            print " F", sorted(F)

        Ft = self.row_echelon(F)
        if self.protocol:
            print " Ft", sorted(Ft)

        LMF = LM(F)
        Ftp = list(set([f for f in Ft if f.lm() not in LMF]))

        if self.protocol:
            print " Ftp", sorted(Ftp)
        if not no_update:
            # maintain the f.lm()=>f dictionary
            self.Ftd.append( dict([(f.lm(),f) for f in Ft]) )
        return Ftp,F

    def symbolic_preprocessing(self,L,G,Fset):
        """
        INPUT:
            L -- a finite subset of M x R[x]
            G -- a finite subset of R[x
            F -- (F_k)k=1,\dots,(d-1), where F_k is finite subset of R[x]

        OUTPUT:
            a finite subset of R[x]
        """
        Simplify = self.simplify
        R = self.ring
        Mul = lambda (m,f): m*f

        F = set([Mul(Simplify(m,f,Fset)) for (m,f) in L])
        if self.protocol:
            print "  F",sorted(F)

        Done = LM(F)

        if self.protocol:
            print "  Done",sorted(Done)
        
        M = set([m for f in F for m in f.monomials()])

        if self.protocol:
            print "  T(F)",sorted(M)

        MdivDone = M.difference(Done)
        zero = R(0)
        G = tuple(G)

        while MdivDone != set():#M != Done
            #m = M.difference(Done).pop()
            m = MdivDone.pop()
            Done.add(m)
            t,g = self.ring.monomial_reduce(m,G)
            if t!=zero:
                tg = Mul(Simplify(t,g,Fset))
                F.add(tg)
                # M = set([m for f in F for m in f.monomials()])
                for tgm in tg.monomials():
                    M.add(tgm)
                    if tgm not in Done:
                        MdivDone.add(tgm)
        return F

    def simplify(self,t,f,F):
        """
        INPUT:
            t -- \in M a monomial
            f -- \in R[x] a polynomial
            F -- (F_k)k=1,\dots,(d-1), where F_k is finite subset of R[x]

        OUTPUT:
            a non evaluated product, i.e. an element of M x R[x]
        """
        for u in sorted(self.ring.monomial_all_divisors(t),reverse=True):
            uf = u*f
            for j in F:
                if uf in F[j]:
                    # F~_j is the row echelon form of F_j w.r.t. <
                    # there exists a (unique) p \in F~_j such that LM(p) = LM(u*f)
                    p = self.Ftd[j][uf.lm()]
                    if u!=t:
                        return self.simplify(self.ring.monomial_quotient(t,u),p,F) #t/u
                    else:
                        return (self.ring(1),p)
        return (t,f)
    

    def first(self,G):
        """
        Returns the largest element of G.

        INPUT:
            G -- a finite subset of G

        OUTPUT:
            a polynomial \in G
        
        """
        mg = G[0]
        mm = mg.lm()
        for g in G:
            if g.lm() > mm:
                mm = g.lm()
                mg = g
        return mg

def LM(F):
    """
    """
    if isinstance(F,(list,set,tuple)):
        return set([f.lm() for f in F])
    else:
        return F.lm()