# -*- coding: utf-8 -*-
"""
Jean-Charles Faugère's F5 Algorithm in F4-Style.

This variant of F5 proceed degree-by-degree in the outer loop instead
of by index of the generators. Futhermore, this variant uses linear
algebra to perform the top reductions.

AUTHORS:
- Martin Albrecht and John Perry (2009-01): use linear algebra to
  perform reductions in F5 proper
- Martin Albrecht (2009-04): proceed degree by degree, performance
  improvements, clean-ups, documentation

EXAMPLE::

    sage: execfile('f5_2.py')
    sage: P = PolynomialRing(GF(32003),3,'x')
    sage: I = sage.rings.ideal.Katsura(P).homogenize()
    sage: gb = f5_2(I); f5_2.print_stats()
     3    1    4 x    7,    4,    0
                   |L|:    4
               L is GB: True
    reductions to zero:    0
           max. degree:    3

    sage: I = sage.rings.ideal.Cyclic(P).homogenize()
    sage: gb = f5_2(I); f5_2.print_stats()
     0    0
                   |L|:    3
               L is GB: True
    reductions to zero:    0
           max. degree:    3

    sage: P = PolynomialRing(GF(32003),4,'x')
    sage: I = sage.rings.ideal.Katsura(P).homogenize()
    sage: gb = f5_2(I); f5_2.print_stats()
     3    2    8 x   15,    8,    0
     4    2   12 x   19,   12,    0
                   |L|:    9
               L is GB: True
    reductions to zero:    0
           max. degree:    4

    sage: I = sage.rings.ideal.Cyclic(P).homogenize()
    sage: gb = f5_2(I); f5_2.print_stats()
     4    1    3 x    5,    2,    1
     5    2    5 x   11,    5,    0
     6    1    2 x    7,    2,    0
                   |L|:    8
               L is GB: True
    reductions to zero:    1
           max. degree:    6

    sage: P = PolynomialRing(GF(32003),5,'x')
    sage: I = sage.rings.ideal.Katsura(P).homogenize()
    sage: gb = f5_2(I); f5_2.print_stats()
     3    4   15 x   29,   15,    0
     4    3   27 x   42,   27,    0
     5    2   27 x   42,   27,    0
                   |L|:   16
               L is GB: True
    reductions to zero:    0
           max. degree:    5

    sage: I = sage.rings.ideal.Cyclic(P).homogenize()
    sage: gb = f5_2(I); f5_2.print_stats()
     4    1    6 x   22,    6,    0
     5    2   13 x   31,   13,    0
     6    6   19 x   38,   19,    0
     7    6   21 x   39,   21,    0
     8    7   27 x   45,   27,    0
     9    8   15 x   33,   15,    0
    10    4
    11    6    7 x   21,    7,    0
    12    3    9 x   23,    9,    0
    13    3    7 x   21,    7,    0
                   |L|:   39
               L is GB: True
    reductions to zero:    0
           max. degree:   13

    sage: P = PolynomialRing(GF(32003),6,'x')
    sage: I = sage.rings.ideal.Katsura(P).homogenize()
    sage: gb = f5_2(I); f5_2.print_stats()
     3    5   23 x   48,   23,    0
     4    7   62 x   92,   62,    0
     5    4   94 x  125,   94,    0
     6    2   56 x   87,   56,    0
                   |L|:   31
               L is GB: True
    reductions to zero:    0
           max. degree:    6

    sage: I = sage.rings.ideal.Cyclic(P).homogenize()
    sage: gb = f5_2(I); f5_2.print_stats()
     4    1    7 x   37,    7,    0
     5    2   15 x   57,   15,    0
     6    6   40 x   98,   40,    0
     7   18   55 x  118,   55,    0
     8   28  100 x  165,   99,    1
     9   48  153 x  211,  149,    4
    10   72  150 x  207,  143,    7
    11   92  166 x  225,  162,    4
    12   55   78 x  132,   78,    0
    13   30   56 x  104,   56,    0
    14   23   54 x   80,   54,    0
    15   34   34 x   60,   34,    0
    16   26   41 x   67,   41,    0
    17   15   35 x   61,   35,    0
    18    2
    19    2
                   |L|:  257
               L is GB: True
    reductions to zero:   16
           max. degree:   17

    sage: P = PolynomialRing(GF(7),3,'x')
    sage: I = sage.rings.ideal.Katsura(P).homogenize()
    sage: gb = f5_2(I); f5_2.print_stats()
     3    1    4 x    7,    4,    0
                   |L|:    4
               L is GB: True
    reductions to zero:    0
           max. degree:    3

    sage: I = sage.rings.ideal.Cyclic(P).homogenize()
    sage: gb = f5_2(I); f5_2.print_stats()
     0    0
                   |L|:    3
               L is GB: True
    reductions to zero:    0
           max. degree:    3

    sage: P = PolynomialRing(GF(7),4,'x')
    sage: I = sage.rings.ideal.Katsura(P).homogenize()
    sage: gb = f5_2(I); f5_2.print_stats()
     3    2   10 x   17,   10,    0
     4    2   16 x   22,   16,    0
                   |L|:    8
               L is GB: True
    reductions to zero:    0
           max. degree:    4

    sage: I = sage.rings.ideal.Cyclic(P).homogenize()
    sage: gb = f5_2(I); f5_2.print_stats()
     4    1    3 x    5,    2,    1
     5    2    5 x   11,    5,    0
     6    1    2 x    7,    2,    0
                   |L|:    8
               L is GB: True
    reductions to zero:    1
           max. degree:    6

    sage: P = PolynomialRing(GF(7),5,'x')
    sage: I = sage.rings.ideal.Katsura(P).homogenize()
    sage: gb = f5_2(I); f5_2.print_stats()
     3    4   15 x   29,   15,    0
     4    3   27 x   42,   27,    0
     5    2   26 x   41,   26,    0
                   |L|:   16
               L is GB: True
    reductions to zero:    0
           max. degree:    5

    sage: I = sage.rings.ideal.Cyclic(P).homogenize()
    sage: gb = f5_2(I); f5_2.print_stats()
     4    1    6 x   22,    6,    0
     5    2   13 x   31,   13,    0
     6    6   19 x   38,   19,    0
     7    6   21 x   39,   21,    0
     8    7   27 x   45,   27,    0
     9   10   14 x   32,   14,    0
    10    3
    11    6    7 x   21,    7,    0
    12    3    8 x   21,    8,    0
    13    3    7 x   21,    7,    0
                   |L|:   39
               L is GB: True
    reductions to zero:    0
           max. degree:   13

    sage: P = PolynomialRing(GF(7),6,'x')
    sage: I = sage.rings.ideal.Katsura(P).homogenize()
    sage: gb = f5_2(I); f5_2.print_stats()
     3    5   23 x   48,   23,    0
     4    7   62 x   92,   62,    0
     5    7   72 x  102,   72,    0
     6    5   92 x  121,   92,    0
     7    1   85 x  114,   85,    0
     8    1
                   |L|:   41
               L is GB: True
    reductions to zero:    0
           max. degree:    7

    sage: I = sage.rings.ideal.Cyclic(P).homogenize()
    sage: gb = f5_2(I); f5_2.print_stats()
     4    1    7 x   37,    7,    0
     5    2   15 x   57,   15,    0
     6    6   40 x   98,   40,    0
     7   18   53 x  116,   53,    0
     8   28  100 x  165,   99,    1
     9   48  153 x  211,  149,    4
    10   72  144 x  201,  137,    7
    11   86  163 x  222,  159,    4
    12   55   71 x  125,   71,    0
    13   29   56 x  104,   56,    0
    14   22   54 x   80,   54,    0
    15   30   34 x   60,   34,    0
    16   21   39 x   65,   39,    0
    17   16   39 x   65,   39,    0
    18    4
    19    1
                   |L|:  248
               L is GB: True
    reductions to zero:   16
           max. degree:   17

    sage: sr = mq.SR(1,1,1,4)
    sage: F,s = sr.polynomial_system()
    sage: I = Ideal(F.gens()).homogenize()
    sage: gb = f5_2(I); f5_2.print_stats()
     3   22   46 x   35,   31,   15
     4   33   64 x   35,   33,   31
     5   14   22 x   17,   15,    7
                   |L|:   88
               L is GB: True
    reductions to zero:   53
           max. degree:    5
"""

from itertools import chain

def compare_by_degree(f,g):
    """
    Compare ``f`` and ``g`` with respect to their degree first and
    only if those match w.r.t. the monomial ordering.

    INPUT:
    
    - ``f`` - a polynomial
    - ``g`` - a polynomial
    """
    if f.total_degree() > g.total_degree():
        return 1
    elif f.total_degree() < g.total_degree():
        return -1
    else:
        return cmp(f, g)

class F5_2:
    """
    Jean-Charles Faugère's F5 Algorithm in F4-Style.

    This variant of F5 proceed degree-by-degree in the outer loop
    instead of by index of the generators. Futhermore, this variant
    uses linear algebra to perform the top reductions.
    """
    def __init__(self):
        self.L = [] # labels
        self.Rules = [] # rewriting rules
        self.verbose = 0
        self.zero_reductions = 0 # we count reductions to zero
        
    def __call__(self, F, D=None, proof=False):
        """
        Compute a Gröbner basis for the input system ``F``. 

        If ``D`` is not ``None`` a D-Gröbner Basis is computed
        instead. If ``proof`` is ``True`` we invoke the Buchberger
        criterion whether the given basis is a Gröbner basis to force
        termination. However, we expect F5 to terminate on all inputs
        and thus ``proof`` defaults to ``False`` for now.

        INPUT:

        - ``F`` - a list of polynomials
        - ``D`` - a maximal degree (default: ``None``)
        - ``proof`` - whether to provably force termination (default: ``False``)
        """
        self.__init__()

        L, Rules = self.L, self.Rules
        critical_pair = self.critical_pair
        select = self.select
        compute_spols = self.compute_spols
        sig, poly = self.sig, self.poly
        is_top_reducible = self.is_top_reducible

        try:
            F = F.gens()
        except AttributeError:
            pass

        F = Ideal(F).interreduced_basis()
        F = sorted(F, cmp=compare_by_degree)
        m = len(F)

        assert(all(f == f.homogenize() for f in F))

        R = F[0].parent()

        # First, we encode that each f_i is generated by itself
        P,G = set(), []
        for i,f in enumerate(F):
            L.append( (Signature(R(1),i), f*f.lc()**(-1)) )
            Rules.append([])
            self.add_rule(Signature(R(1),i),i)
            P = P.union(set([critical_pair(i, g, G) for g in chain(*G)]))
            G.append([i])

        d = 0
        while P: # as long as we have S-polynomials
            # TODO: this is a technicality, get rid of this
            P = set([p for p in P if p != tuple()])

            d = self.minimal_d(P) # get the minimal d
            Pd = select(P,d) # and select all pairs with minimal degree
            print "%2d %4d"%(d, len(Pd)),

            P = P.difference(Pd)

            # we allow a degree bound D
            if D is not None and d > D:
                print
                continue
            
            # compute the S-polynomials
            S = compute_spols(Pd)

            if len(S):
                self.degree_reached = d

            # and perform top-reductions
            Stilde = self.reduction(S, G)
            print

            # we add each new polynomial, and repeat
            for h in Stilde:
                Pnew = set([critical_pair(h, g, G) for g in chain(*G)])
                P = P.union(Pnew)
                G[sig(h)[1]].append(h)
                
            P = set([p for p in P if p != tuple()])
            sys.stdout.flush()
            if proof and self.terminate(P, G):
                break

        return [poly(f) for f in chain(*G)]

    def critical_pair(self, k, l, G):
        """
        Return the critical pair for the polynomials ``k`` and ``l``
        indexed in ``L`` iff the F5 criteria pass.

        Otherwise, return the empty tuple.

        INPUT:
        
        - ``k`` - a polynomial index for ``L``
        - ``l`` - a polynomial index for ``L``
        - ``G`` - the intermediate Gröbner basis indexed in ``L``

        adapted from Justin Gash (p.51): 
       
        'It is the subroutine critical_pair that is responsible for
         imposing the F5 Criterion from Theorem 3.3.1. Note that in
         condition (3) of Theorem 3.3.1, it is required that all pairs
         (r_i, r_j) be normalized. The reader will recall from
         Definition 3.2.2 that a pair is normalized if: 
         
         (1) S(k) = m_k*F_{e_k} is not top-reducible by <f_0, ..., f_{e_k}-1> 

         (2) S(l) = m_l*F_{e_l} is not top-reducible by <f_0, ..., f_{e_l}-1>
         
         (3) S(m_k*k) > S(m_l*l)

         If these three conditions are not met in ``critical_pair()``
         (note that the third condition will always be met because
         ``cirtical_pair()`` forces it to be met), the nominated
         critical pair is dropped and () is returned.

         Once we have collected the nominated critical pairs that pass
         the F5 criterion test of ``critical_pair(()``, we send them
         to ``compute_spols()``.'
        """
        poly = self.poly
        sig = self.sig
        is_top_reducible = self.is_top_reducible
        is_rewritable = self.is_rewritable
        LCM = lambda f,g: f.parent().monomial_lcm(f,g)

        tk = poly(k).lt()
        tl = poly(l).lt()
        t = LCM(tk, tl)
        uk = t//tk
        ul = t//tl
        mk, ek = sig(k)
        ml, el = sig(l)

        # they are are same
        if ek == el and uk*mk == ul*ml:
            return tuple()
        
        if is_top_reducible(uk*mk, G[:ek]):
            return tuple()

        if is_top_reducible(ul*ml, G[:el]):
            return tuple()

        # this check is in compute_spols() again, but we can filter
        # some stuff out here already
        if is_rewritable(uk, k) or is_rewritable(ul, l):
            return tuple()

        # preserve order
        if uk * sig(k) < ul * sig(l):
            uk, ul = ul, uk
            k, l = l, k
        return (t,uk,k,ul,l)

    def minimal_d(self, P):
        if len(P) == 0:
            return 0
        d = iter(P).next()[0].total_degree()
        for (t,_,_,_,_) in P:
            if t.total_degree() < d:
                d = t.total_degree()
        return d

    def select(self, P, d):
        return set([p for p in P if p[0].total_degree() == d])

    def compute_spols(self, P):
        poly = self.poly
        sig = self.sig
        spol = self.spol
        is_rewritable = self.is_rewritable
        is_top_reducible = self.is_top_reducible
        add_rule = self.add_rule

        L = self.L

        S = list()
        P = sorted(P, key=lambda x: x[0])
        for (t,u,k,v,l) in P:
            if not is_rewritable(u,k) and not is_rewritable(v,l):
                s = spol(poly(k), poly(l))
                if s != 0:
                    s = s.lc()**-1 * s
                    L.append( (u * sig(k), s) )
                add_rule(u * sig(k), len(L)-1)
                if s != 0:
                    S.append(len(L)-1)
        S = sorted(S, key=lambda f: sig(f))
        return S

    def reduction(self, S, G):
        """
        INPUT:

        - ``S`` - a list of components of S-polynomials
        - ``G`` - the intermediate Gröbner basis
        """
        L = self.L
        add_rule = self.add_rule
        poly = self.poly

        F = self.symbolic_preprocessing(S, G)
        Ft = self.gauss_elimination(F)

        Ret = []

        for k, (s, p, i) in enumerate(Ft):
            if i < 0 and p.lm() == F[k][1].lm():
                continue # ignore unchanged new polynomials
            elif i >= 0:
                assert(L[i][0] == s)
                L[i] = s, p # update p
                if p != 0:
                    Ret.append(i)
            else:
                L.append( (s,p) ) # we have a new polynomial
                add_rule( s, len(L)-1 )
                if p != 0:
                    Ret.append(len(L)-1)
        return Ret
        
    def symbolic_preprocessing(self, S, G):
        """
        Add polynomials to the set ``S`` such that all possible
        reductors for all elements in ``S`` are available.

        INPUT:

        - ``S`` - a list of components of S-polynomials
        - ``G`` - the intermediate Gröbner basis indexed in ``L``
        """
        poly = self.poly
        sig = self.sig
        L = self.L
        find_reductor = self.find_reductor

        # We add a new marker for each polynomial which encodes
        # whether the polynomial was added by this routine or is an
        # original input polynomial.
        F = [L[k]+(k,) for k in S]
        Done = set() # we already added the S-polynomials

        # the set of all monomials
        M = set([m for (sig_f, poly_f, i_f) in F for m in poly_f.monomials()])

        while M != Done:
            M = sorted(M)
            for i,m in enumerate(M):
                if m not in Done:
                    break
            M = set(M)
            Done.add(m)

            # we need to find the polynomial with the minimal
            # signature which has the monomial m, alternatively we
            # could just use the signature of the polynomial it comes
            # from.
            ms = self.minimal_signature(m, F)
            
            t, g = find_reductor(m, ms, G)
            if t!=0:
                F.append( (t*g[0], t*g[1], -1) )
                M = M.union((t*g[1]).monomials())
        return sorted(F, key=lambda f: f[0]) # sort by signature

    def minimal_signature(self, m, F):
        ms = (1, 10**20)
        for (sig_f, poly_f, _) in F:
            if m in poly_f.monomials() and sig_f < ms:
                ms = sig_f
        return ms

    def find_reductor(self, m, sig_m, G):
        """
        Find a reductor `g_i` for `m` in `G` subject to the F5
        constaints.

        INPUT:
        - ``m`` - a monomial
        - ``sig_m`` - the signature of the smalles f which contains ``m``
        - ``G`` - the intermediate Gröbner basis
        """
        is_rewritable = self.is_rewritable
        is_top_reducible = self.is_top_reducible
        sig = self.sig
        poly = self.poly 

        L = self.L
        R = m.parent()
        for k in chain(*G):
            # Requirement (1) is the normal top-reduction criterion.
            if not R.monomial_divides(poly(k).lm(),m):
                continue
            t =  R.monomial_quotient(m, poly(k).lm())

            # Requirement (2) is making sure that the signature of the
            # reductor is normalized.  Recall that we only want
            # signatures of our polynomials to be normalized - we are
            # discarding non-normalized S-polynomials. If we ignored
            # this condition and our reductor would up having larger
            # signature than S(r_{k_0}), then top_reduction would
            # create a new signed polynomial with our reductor's
            # non-normalized signature. (We might add that, if the
            # reductor had smaller signature than S(r_{k_0}), it would
            # be fine to reduce by it; however, F5 doesn't miss
            # anything by forgoing this opportunity because, by Lemma
            # 3.2.1 (The Normalization Lemma), there will be another
            # normalized reductor with the same head term and smaller
            # signature.
            if is_top_reducible(t * sig(k)[0], G[:sig(k)[1]]):
                continue

            # Requirement (3)
            if is_rewritable(t, k):
                continue

            #  Requirement (4) is a check that makes sure we don't
            #  reduce by something that has the same signature as
            #  m. Recall that we want all signed polynomials used
            #  during the run of F5 to be admissible. If we reduced by
            #  a polynomial that has the same signature, we would be
            #  left with a new polynomial for which we would have no
            #  idea what the signature is. The act of reduction would
            #  have certainly lowered the signature, thus causing
            #  admissibility to be lost. (We will comment on this
            #  requirement later in subsection 3.5. With a little
            #  care, we can loosen this requirement.)
            if sig_m == t*sig(k):
                continue
            return t, L[k]
        return 0, -1
        
    def gauss_elimination(self, F1):
        """
        Perform permuted F5-style Gaussian elimination on ``F1``.

        INPUT:

        - ``F1`` - a list of tuples ``(sig, poly, idx)``
        """
        F = [f[1] for f in F1]

        if len(F) == 0:
            return F
        A,v = mq.MPolynomialSystem(F).coefficient_matrix()
        self.zero_reductions += A.nrows()-A.rank()
        print "%4d x %4d, %4d, %4d"%(A.nrows(), A.ncols(), A.rank(), A.nrows()-A.rank()),
        nrows, ncols = A.nrows(), A.ncols()
        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:
                            minus_b = -A[i,c]
                            A.add_multiple_of_row(i, r, minus_b, c)
                    break

        F = (A*v).list()
        return [(F1[i][0],F[i],F1[i][2]) for i in xrange(len(F))]

    def poly(self, i):
        return self.L[i][1]

    def sig(self, i):
        return self.L[i][0]

    def spol(self, f, g):
        LM = lambda f: f.lm()
        LT = lambda f: f.lt()
        LCM = lambda f,g: f.parent().monomial_lcm(f,g)
        return LCM(LM(f),LM(g)) // LT(f) * f - LCM(LM(f),LM(g)) // LT(g) * g

    def is_top_reducible(self, t, l):
        R = t.parent()
        poly = self.poly
        for g in chain(*l):
            if R.monomial_divides(poly(g).lm(),t):
                return True
        return False

    def add_rule(self, s, k):
        self.Rules[s[1]].append( (s[0],k) )

    def is_rewritable(self, u, k):
        j = self.find_rewriting(u, k)
        return j != k

    def find_rewriting(self, u, k):
        """
        INPUT:
        
        - ``k`` - an index in L
        - ``u`` - a monomial
        """

        divides = lambda x,y: x.parent().monomial_divides(x,y)
        Rules = self.Rules
        mk, v = self.sig(k)
        for ctr in reversed(xrange(len(Rules[v]))):
            mj, j = Rules[v][ctr]
            if divides(mj, u * mk):
                return j
        return k

    def terminate(self, P, G):
        I = Ideal([self.poly(f) for f in chain(*G)]).interreduced_basis()
        for (t,u,k,v,l) in P:
            if not self.is_rewritable(u,k) and not self.is_rewritable(v,l):
                s = self.spol(self.poly(k), self.poly(l))
                if s.reduce(I) != 0:
                    return False
        return True

    def print_stats(self):
        print "               |L|: %4d"%len(self.L)
        print "           L is GB: %s"%Ideal([f[1] for f in self.L]).basis_is_groebner()
        print "reductions to zero: %4d"%self.zero_reductions
        print "       max. degree: %4d"%self.degree_reached

from UserList import UserList

class Signature(UserList):
    def __init__(self, multiplier, index):
        """
        Create a new signature from the mulitplier and the index.
        """
        UserList.__init__(self, (multiplier, index))
         
    def __lt__(self, other):
        """
        """
        if self[1] < other[1]:
            return True
        elif self[1] > other[1]:
            return False
        else:
            return (self[0] < other[0])

    def __eq__(self, other):
        return self[0] == other[0] and self[1] == other[1]
    
    def __neq__(self, other):
        return self[0] != other[0] or self[1] != other[1]
  
    def __rmul__(self, other):
        if isinstance(self, Signature):
            return Signature(other * self[0], self[1])
        else:
            raise TypeError

    def __hash__(self):
        return hash(self[0]) + hash(self[1])

f5_2  = F5_2()