# -*- coding: utf-8 -*- """ Jean-Charles Faugère's F5 Algorithm. These implementations are heavily inspired by John Perry's pseudocode and Singular implementation of these algorithms. See http://www.math.usm.edu/perry/Research/ for details. The docstrings are almost verbatim copies from Just Gash's explanations for each F5 function in his thesis: "On Efficient Computation of Gröbner Bases". Note that Justin begins at f_m while we begin at f_0, e.g. the first GB we calculate is while Justin calculates first. AUTHOR: -- 20081013 Martin Albrecht (initial version based on John Perry's pseudocode) -- 20081013 John Perry (loop from 0 to m-1 instead of m-1 to 0) -- 20090112 Martin Albrecht (F5SansRewriting) -- 20090124 Martin Albrecht and John Perry (F4F5) -- 20090126 John Perry (correction to compute_spols) -- 20090414 John Perry (new_reduction, alt_reduction, and associated functions) EXAMPLE: sage: execfile('f5.py') sage: R. = PolynomialRing(GF(29)) sage: I = R* [3*x^4*y + 18*x*y^4 + 4*x^3*y*z + 20*x*y^3*z + 3*x^2*z^3, \ 3*x^3*y^2 + 7*x^2*y^3 + 24*y^2*z^3, \ 12*x*y^4 + 17*x^4*z + 27*y^4*z + 11*x^3*z^2] sage: J = I.homogenize() sage: f5 = F5() # original F5 sage: gb = f5(J) Increment 1 Processing 1 pairs of degree 7 Processing 1 pairs of degree 9 Processing 1 pairs of degree 11 Ended with 4 polynomials Increment 2 Processing 1 pairs of degree 6 Processing 2 pairs of degree 7 Processing 5 pairs of degree 8 Processing 7 pairs of degree 9 Processing 11 pairs of degree 10 Processing 10 pairs of degree 11 verbose 0 (...: f5.py, top_reduction) Polynomial 29 reduced to zero. verbose 0 (...: f5.py, top_reduction) Polynomial 27 reduced to zero. verbose 0 (...: f5.py, top_reduction) Polynomial 25 reduced to zero. Processing 3 pairs of degree 12 Processing 4 pairs of degree 13 Processing 1 pairs of degree 14 Ended with 27 polynomials sage: f5.zero_reductions, len(gb) (3, 18) sage: f5 = F5R() # F5 with interreduced B sage: gb = f5(J) Increment 1 Processing 1 pairs of degree 7 Processing 1 pairs of degree 9 Processing 1 pairs of degree 11 Ended with 4 polynomials Increment 2 Processing 1 pairs of degree 6 Processing 2 pairs of degree 7 Processing 5 pairs of degree 8 Processing 7 pairs of degree 9 Processing 11 pairs of degree 10 Processing 10 pairs of degree 11 verbose 0 (...: f5.py, top_reduction) Polynomial 29 reduced to zero. verbose 0 (...: f5.py, top_reduction) Polynomial 27 reduced to zero. verbose 0 (...: f5.py, top_reduction) Polynomial 25 reduced to zero. Processing 3 pairs of degree 12 Processing 4 pairs of degree 13 Processing 1 pairs of degree 14 Ended with 27 polynomials sage: f5.zero_reductions, len(gb) (3, 18) sage: f5 = F5C() # F5 with interreduced B and Gprev sage: gb = f5(J) Increment 1 Processing 1 pairs of degree 7 Processing 1 pairs of degree 9 Processing 1 pairs of degree 11 Ended with 4 polynomials Increment 2 Processing 1 pairs of degree 6 Processing 2 pairs of degree 7 Processing 5 pairs of degree 8 Processing 7 pairs of degree 9 Processing 11 pairs of degree 10 Processing 10 pairs of degree 11 verbose 0 (...: f5.py, top_reduction) Polynomial 29 reduced to zero. verbose 0 (...: f5.py, top_reduction) Polynomial 27 reduced to zero. verbose 0 (...: f5.py, top_reduction) Polynomial 25 reduced to zero. Processing 3 pairs of degree 12 Processing 4 pairs of degree 13 Processing 1 pairs of degree 14 Ended with 27 polynomials sage: f5.zero_reductions, len(gb) (3, 18) sage: Ideal(gb).basis_is_groebner() True sage: f5 = F4F5() # F5-style F5 sage: gb = f5(J) Increment 1 Processing 1 pairs of degree 7 5 x 13, 5, 0 Processing 1 pairs of degree 9 14 x 29, 14, 0 Processing 1 pairs of degree 11 Ended with 4 polynomials Increment 2 Processing 1 pairs of degree 6 6 x 18, 6, 0 Processing 2 pairs of degree 7 11 x 23, 11, 0 Processing 5 pairs of degree 8 18 x 27, 18, 0 Processing 7 pairs of degree 9 19 x 23, 19, 0 Processing 11 pairs of degree 10 15 x 15, 15, 0 Processing 10 pairs of degree 11 14 x 11, 11, 3 Processing 3 pairs of degree 12 Processing 4 pairs of degree 13 Processing 1 pairs of degree 14 Ended with 27 polynomials sage: f5.zero_reductions, len(gb) (3, 18) sage: Ideal(gb).basis_is_groebner() True For additional diagnostics there are a number of commented commands. To count the number of reductions, one can uncomment commands to "interreduce" and comment out commands with "reduced_basis"; also uncomment commands with "normal_form" and comment out commands with "reduce". """ from sage.misc.cachefunc import cached_method divides = lambda x,y: x.parent().monomial_divides(x,y) LCM = lambda f,g: f.parent().monomial_lcm(f,g) LM = lambda f: f.lm() LT = lambda f: f.lt() def compare_by_degree(f,g): if f.total_degree() > g.total_degree(): return 1 elif f.total_degree() < g.total_degree(): return -1 else: return cmp(f.lm(),g.lm()) class F5: """ Jean-Charles Faugère's F5 Algorithm. """ def __init__(self, F=None): if F is not None: self.Rules = [[]] self.L = [0] self.zero_reductions = 0 self.reductions = 0 def poly(self, i): return self.L[i][1] def sig(self, i): return self.L[i][0] def __call__(self, F): if isinstance(F, sage.rings.polynomial.multi_polynomial_ideal.MPolynomialIdeal): F = F.interreduced_basis() else: F = Ideal(list(F)).interreduced_basis() if not all(f.is_homogeneous() for f in F): F = Ideal(F).homogenize() F = F.gens() return self.basis(F) def basis(self, F): """ F5's main routine. Computes a Gröbner basis for F. INPUT: F -- a list of polynomials OUTPUT: G -- a list of polynomials; a Gröbner basis for """ poly = self.poly incremental_basis = self.incremental_basis self.__init__(F) Rules = self.Rules L = self.L m = len(F) F = sorted(F, cmp=compare_by_degree) f0 = F[0] L[0] = (Signature(1, 0), f0*f0.lc()**(-1)) Rules.append(list()) Gprev = set([0]) B = [f0] for i in xrange(1,m): print "Increment", i f = F[i] L.append( (Signature(1,i), f*f.lc()**(-1)) ) Gcurr = incremental_basis(i, B, Gprev) if any(poly(lambd) == 1 for lambd in Gcurr): return set(1) Gprev = Gcurr B = [poly(l) for l in Gprev] #return B return Ideal([poly(l) for l in Gprev]).interreduced_basis() #return self.interreduce(B) def interreduce(self, RF): F = list(RF) for each in xrange(len(F)): F[each] = F[each]*F[each].lc()**(-1) i = 0 while i < len(F): reduceme = F.pop(0) reduced = False for j in xrange(len(F)): quo, rem = self.divide(reduceme,F[j]) reduceme = rem if (quo != 0) and (rem != 0): reduceme = rem*rem.lc()**(-1) j = -1 reduced = True self.reductions += 1 if (reduceme != 0): F.append(reduceme) if reduced: i = -1 i = i + 1 return F def minimal_basis(self, RF): F = list(RF) for each in xrange(len(F)): F[each] = F[each]*F[each].lc()**(-1) i = 0 R = F[0].parent() while i < len(F): reduceme = F.pop(0) top_reducible = False for j in xrange(len(F)): if R.monomial_divides(F[j].lm(), reduceme.lm()): top_reducible = True break if not top_reducible: F.append(reduceme) else: i = -1 i = i + 1 return F def normal_form(self,f,B): remainder = f quotient = [0 for each in B] i = 0 while i < len(B): quo, rem = self.divide(remainder, B[i]) remainder = rem if quo != 0: i = -1 i = i + 1 return remainder def divide(self,dividend,divisor): remainder = dividend quotient = 0 mons = remainder.monomials() coeffs = remainder.coefficients() t = divisor.lm() c = divisor.lc() i = 0 while (remainder != 0) and (i < len(mons)): if t.divides(mons[i]): self.reductions += 1 quotient = quotient + (mons[i]/t*coeffs[i]/c).numerator() remainder = remainder - (mons[i]/t*coeffs[i]/c*divisor).numerator() mons = remainder.monomials() coeffs = remainder.coefficients() else: i = i + 1 return quotient, remainder def incremental_basis(self, i, B, Gprev): """ adapted from Justin Gash (p.49): 'This is the portion of the algorithm that is called (m-1) times and at the end of each call a new Gröbner basis is produced. After the first call to this algorithm, the Gröbner basis for is returned. In general, after the k-th call to this algorithm, the Gröbner basis for . This is why F5 is called an iterative algorithm. The process used by this algorithm is similar to many Gröbner basis algorithms: it moves degree-by-degree; it generates a new set of S- polynomials S to consider; it reduces this new set S of S-polynomials by G_curr and G_prev.' """ L = self.L critical_pair = self.critical_pair compute_spols = self.compute_spols reduction = self.reduction Rules = self.Rules curr_idx = len(L) - 1 Gcurr = Gprev.union([curr_idx]) Rules.append( list() ) P = reduce(lambda x,y: x.union(y), [critical_pair(curr_idx, j, i, Gprev) for j in Gprev], set()) print len(P), "critical pairs" while len(P) != 0: d = min(t.degree() for (t,k,u,l,v) in P) Pd = [(t,k,u,l,v) for (t,k,u,l,v) in P if t.degree() == d] print "Processing", len(Pd), "pairs of degree", d, "of", len(P), "total" #print Pd P = P.difference(Pd) S = compute_spols(Pd) print len(S), "polynomials generated" R = reduction(S, B, Gprev, Gcurr) # These alternate reduction algorithms do not seem to improve things #R = self.new_reduction(S, B, Gprev, Gcurr) #R = self.alt_reduction(S, B, Gprev, Gcurr) print len(R), "polynomials left" for k in R: P = reduce(lambda x,y: x.union(y), [critical_pair(j, k, i, Gprev) for j in Gcurr], P) #print k, ":", self.sig(k), "lm:", self.poly(k).lm() Gcurr.add(k) print "Ended with", len(Gcurr), "polynomials" #print "Leading monomials:" #print [ self.poly(each).lm() for each in Gcurr ] return Gcurr def critical_pair(self, k, l, i, Gprev): """ 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_0*F_{e_0} is not top-reducible by (2) S(l) = m_1*F_{e_1} is not top-reducible by (3) S(m_0*k) > S(m_1*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 #print "crit_pair(%s,%s,%s,%s)"%(k, l, i, Gprev) #print self.L tk = poly(k).lt() tl = poly(l).lt() t = LCM(tk, tl) u0 = t//tk u1 = t//tl m0, e0 = sig(k) m1, e1 = sig(l) if e0 == e1 and u0*m0 == u1*m1: return set() # Stegers and Gash leave out the == i check, Faugere and Perry # have it. It is unclear for now, whether the check is # necessary. if e0 == i and is_top_reducible(u0*m0, Gprev): #if e0>0 and is_top_reducible(u0*m0, e0-1): return set() if e1 == i and is_top_reducible(u1*m1, Gprev): #if e1>0 and is_top_reducible(u1*m1, e1-1): return set() # This check was introduced by Stegers, it isn't strictly # necessary if is_rewritable(u0, k) or is_rewritable(u1, l): return set() if u0 * sig(k) < u1 * sig(l): u0, u1 = u1, u0 k, l = l, k # selecting pairs by signature, rather than by lcm, seems to improve things return set([(t,k,u0,l,u1)]) #return set([(u0*sig(k)[0],k,u0,l,u1)]) def compute_spols(self, P): """ adapted from Justin Gash (p.51): 'Though at first glance this subroutine may look complicated, compute_spols essentially does one thing: form the new S-polynomials output from critical_pairs as admissible signed polynomials. We note that, because critical_pairs ensured that S(u*k) < S(v*l), we know that the signature of all new polynomials will always be of the form u_L*S(r_{i_L}) in compute_spols.' """ poly = self.poly sig = self.sig spol = self.spol is_rewritable = self.is_rewritable add_rule = self.add_rule L = self.L S = list() P = sorted(P, key=lambda x: x[0]) for (t,k,u,l,v) in P: if not is_rewritable(u,k) and not is_rewritable(v,l): s = spol(poly(k), poly(l)) verbose("generated" + str(len(L)) + "from S-poly of" + str(k) + "and" + str(l), level=1) #print "generated", len(L), "from S-poly of", k, sig(k), "and", l, sig(l) L.append( (u * sig(k), s.lc()**-1 * s) ) add_rule(u * sig(k), len(L)-1) if s != 0: S.append(len(L)-1) S = sorted(S, key=lambda x: sig(x)) return S def spol(self, f, g): return LCM(LM(f),LM(g)) // LT(f) * f - LCM(LM(f),LM(g)) // LT(g) * g def reduction(self, S, B, Gprev, Gcurr): """ adapted from Justin Gash (p.54ff): 'Let's begin our discussion by focusing our attention to the outer layer of the reduction subroutine(s): reduction. The signed polynomial with smallest signature, denoted h, is grabbed and removed from the todo list of polynomial to be reduced. It's normal form, with respect to the previous Gröbner basis, and other information is sent to the sub-subroutine top_reduction. If top_reduction determines that the signed polynomial can be reduced, then nothing will be added to completed and the reduced (still signed) version of h will be placed back into todo. If no top reduction is possible, h is made monic by K-multiplication and the resulting signed polynomial is placed in completed. This description of reduction seems very similar to other reduction routines from other algorithms. The difference lies in the phrase, "If top_reduction determines that the signed polynomial can be reduced ..."' """ L = self.L sig = self.sig poly = self.poly top_reduction = self.top_reduction to_do = S completed = set() while len(to_do): k, to_do = to_do[0], to_do[1:] #print "Reducing", str(k), L[k] h = poly(k).reduce(B) #h = self.normal_form(poly(k),B) L[k] = (sig(k), h) newly_completed, redo = top_reduction(k, Gprev, Gcurr.union(completed)) completed = completed.union( newly_completed ) #if k in newly_completed: # print "completed", k, "lm", poly(k).lt() for j in redo: # insert j in to_do, sorted by increasing signature to_do.append(j) to_do.sort(key=lambda x: sig(x)) return completed def top_reduction(self, k, Gprev, Gcurr): """ adapted from Justin Gash (p.55ff): 'We will go through top_reduction step-by-step. If the signed polynomial being examined has polynomial part 0, then there is no data left in that particular signed polynomial - an empty ordered pair is returned. Otherwise top_reduction calls upon another sub-subroutine find_reductor. Essentially, if find_reductor comes back negative, the current signed polynomial is made monic and returned to reduction to be placed in completed. If a top-reduction is deemed possible, then there are two possible cases: either the reduction will increase the signature of polynomial or it won't. In the latter case, the signature of r_{k_0} is maintained, the polynomial portion is top-reduced and the signed polynomial is returned to reduction to be added back into todo; this case corresponds to top-reduction in previous algorithms. In the former case, however, the signature will change. This is marked by adding a new polynomial r_N (our notation here describes N after N was incremented) with appropriate signature based upon the reductor, not S(r_{k_0}). A new rule is added (as I mentioned previously, this will be explained later) and then both r_{k_0} and r_N are sent back to reduction to be added back into todo. This is done because r_N has a different signature than r_{k_0} and r_{k_0} might still be reducible by another signed polynomial. """ from sage.misc.misc import verbose find_reductor = self.find_reductor add_rule = self.add_rule poly = self.poly sig = self.sig L = self.L #print "in top-reduction for", k, sig(k) if poly(k) == 0: verbose("Polynomial " + str(k) + " reduced to zero.",level=0) self.zero_reductions += 1 return set(),set() p = poly(k) J = find_reductor(k, Gprev, Gcurr) if J == set(): L[k] = ( sig(k), p * p.lc()**(-1) ) return set([k]),set() j = J.pop() q = poly(j) u = p.lt()//q.lt() #print "reducing by", j, sig(j), u*sig(j) p = p - u*q self.reductions += 1 if p != 0: p = p * p.lc()**(-1) if u * sig(j) < sig(k): L[k] = (sig(k), p) return set(), set([k]) else: verbose("generated" + str(len(L)) + "from reduction of" + str(k) + "by" + str(j), level=1) #print "Generated", len(L), "from reduction of", k, sig(k), "by", j, sig(j) L.append((u.lm() * sig(j), p)) add_rule(u.lm() * sig(j), len(L)-1) return set(), set([k, len(L)-1]) def find_reductor(self, k, Gprev, Gcurr): """ adapted from Justin Gash (p.56ff): 'For a previously added signed polynomial in G_curr to become a reductor of r_{k_0}, it must meet four requirements: (1) u = HT(r_{k_0})/HT(r_{k_j}) \in T (2) NF(u_{t_j}, G_curr) = u_{t_j} (3) not is_rewriteable(u, r_{k_j}) (4) u_{t_j} F_{k_j} = S(r_{k_0}) We will go through each requirement one-by-one. Requirement (1) is simply the normal top-reduction requirement. The only thing of note here is that, in testing for the top-reducibility, u is assigned a particular value to be used in subsequent tests. 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 re- ductor wound 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.) Requirement (3) will be discussed when we discuss is_rewriteable. That discussion is approaching rapidly. Requirement (4) is a check that makes sure we don't reduce by something that has the same signature as r_{k_0} . 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.) """ is_rewritable = self.is_rewritable is_top_reducible = self.is_top_reducible poly = self.poly sig = self.sig t = poly(k).lt() for j in Gcurr: tprime = poly(j).lt() if divides(tprime,t): u = t // tprime mj, ej = sig(j) if u * sig(j) != sig(k) and not is_rewritable(u, j) \ and not is_top_reducible(u*mj, Gprev): #and not is_top_reducible(u*mj, -1): return set([j]) return set() def new_reduction(self, S, B, Gprev, Gcurr): """ For each k in S, this reduces poly(k) with respect to Gcurr, with the following exceptions: (1) new_reduction calls find_reducers instead of find_reductor. (2) find_reducers tries to identify all polynomials whose leading monomial might divide a monomials generated during a reduction pass of each polynomial. (3) Those that are signature-safe are returned to new_reduction; those that are not signature-safe are ignored. (4) poly(k) is immediately reduced as far as possible by all the reducers in B and reducers. This terminates for the same reason that reduction and top_reduction terminate. However, it seems to take much longer, even though the reductions are all passed off to Singular. The reason appears to be that the reductions that are *not* type-safe cause F5 to dawdle longer at a given degree, generating S-polynomials that correspond to these reductions. This process tends to generate more polynomials than the usual F5 routine. """ from sage.misc.misc import verbose L = self.L sig = self.sig poly = self.poly top_reduction = self.top_reduction to_do = S completed = set() Gnew = Gcurr.difference(Gprev) #print "updates done with", Gnew while len(to_do): k, to_do = to_do[0], to_do[1:] h = poly(k) old_t = 0 while old_t != h.lm(): old_t = h.lm() # METHOD 1 reducers = self.find_reducers(k, Gcurr.union(completed)) # METHOD 2 #too_large = self.update_reducers(k) #reducers = self.reducers #reducer_list = self.reducer_list # METHOD 3 #reducers = self.symbolic_preprocessing(k, Gcurr) #print "reducing", k #h = h.reduce(B.union(reducers)) h = h.reduce(reducers) L[k] = (sig(k), h) #print "finished" if h != 0: completed = completed.union([k]) #self.reducer_list.append((1, k)) #self.reducers.append(poly(k)) else: verbose("Polynomial " + str(k) + " reduced to zero.",level=0) #for each in too_large: # self.reducer_list.append(each[0]) # self.reducers.append(each[1]) return completed def find_reducers(self, k, Gcurr): """ Computes a list of reducers p for poly(k) such that p.lm() is a monomial that *might* appear in the reduction of poly(k) and reduction by poly(j) is signature-safe (smaller than sig(k), not rewritable, not top-reducible by Gprev). """ poly = self.poly sig = self.sig result = [] # print k R = poly(k).parent() monomial_list = set(poly(k).monomials()) while len(monomial_list) != 0: #mon = monomial_list.pop() mon = max(monomial_list, key=lambda x: x.lm()) monomial_list = monomial_list.difference(set([mon])) for j in Gcurr: if R.monomial_divides(poly(j).lm(), mon): u = R.monomial_quotient(mon, poly(j).lm()) mj, ej = sig(j) if u * sig(j) < sig(k) and not self.is_rewritable(u, j) and not self.is_top_reducible(u*mj, -1): #result.append((u*poly(j))) result.append(self.polynomial_multiple(u, j)) monomial_list = monomial_list.union(result[-1].monomials()[1:]) break return result def alt_reduction(self, S, B, Gprev, Gcurr): """ For each k in S, this reduces poly(k) with respect to Gcurr, with the following exceptions: (1) alt_reduction calls alt_find_reducers instead of find_reductor. (2) alt_find_reducers tries to identify all polynomials whose leading monomial might divide a monomials generated during a reduction pass of each polynomial. (3) These are divided into two lists: safe_reducers, polynomials that can be used safely in the reduction of poly(k), and unsafe_reducers, a set of integers j such that poly(j) *might* be usable. (4) poly(k) is immediately reduced as far as possible by all the reducers in B and safe_reducers. (5) Then we test whether poly(k).lm() is top-reducible by any j indexed by unsafe_reducers: for each such j we compute the reduction and add a new polynomial with the larger signature to the basis; this new polynomial is also added to the list of polys to be reduced. This terminates for the same reason that reduction and top_reduction terminate. It terminates more quickly because safe top-reductions of the same signature index are all passed off to Singular, rather than treated one-by-one. It does not generate more polynomials than the usual F5 routine because it does exactly the same thing, only in a way optimized for Singular. """ from sage.misc.misc import verbose L = self.L sig = self.sig poly = self.poly top_reduction = self.top_reduction to_do = S completed = set() while len(to_do): k, to_do = to_do[0], to_do[1:] h = poly(k).reduce(B) L[k] = (sig(k), h) safe_reducers, unsafe_reducers = self.alt_find_reducers(k, Gcurr.union(completed)) h = h.reduce(safe_reducers) L[k] = (sig(k), h) # check for top-reduction by elements of unsafe_reducers if h != 0: completed = completed.union([k]) R = h.parent() # when we loop through unsafe_reducers, we do not stop at the first one that might # divide h.lm(), but generate a new S-poly with each of them for j in unsafe_reducers: if R.monomial_divides(poly(j).lm(), h.lm()): u = R.monomial_quotient(h.lm(), poly(j).lm()) # j might be in unsafe_reducers for a different monomial than u, # so we have to check this reduction again against the criteria if not self.is_rewritable(u, j) and not self.is_top_reducible(u*(sig(j)[0]), -1): c = h.lc() * poly(j).lc()**(-1) L.append( (u * sig(j), h - c * u * poly(j)) ) self.add_rule(u*sig(j), len(L) - 1) verbose("generated" + str(len(L) - 1) + "from top-reduction of" + str(k) + "by" + str(j), level=1) #print "generated", len(L)-1, sig(len(L)-1), "from top-reduction of", k, sig(k), "by", j, sig(j) # we have to reduce the new polynomial, too to_do.append(len(L) - 1) to_do.sort(key=lambda x: sig(x)) #print "finished" else: verbose("Polynomial " + str(k) + " reduced to zero.",level=0) return completed def alt_find_reducers(self, k, Gcurr): """ Computes two lists of reducers of poly(k). Elements of smaller_reducers are polynomials p such that p.lm() is a monomial that *might* appear in the reduction of poly(k) and reduction by poly(j) is signature-safe (smaller than sig(k), not rewritable, not top-reducible by Gprev). Elements of larger_reducers are integers indexing the polynomials p such that p.lm() is a monomial that *might* appear in the reduction of poly(k) and although reduction by poly(j) is not rewritable nor top-reducible by Gprev, it is not signature-safe because u*sig(j) > sig(k). """ poly = self.poly sig = self.sig safe_reducers = [] unsafe_reducers = [] # print k R = poly(k).parent() monomial_list = set(poly(k).monomials()) while len(monomial_list) != 0: #mon = monomial_list.pop() mon = max(monomial_list, key=lambda x: x.lm()) monomial_list = monomial_list.difference(set([mon])) for j in Gcurr: if R.monomial_divides(poly(j).lm(), mon): u = R.monomial_quotient(mon, poly(j).lm()) mj, ej = sig(j) mk, ek = sig(k) # polynomials of smaller signature index are always okay if ej < ek: #safe_reducers.append(u*poly(j)) safe_reducers.append(self.polynomial_multiple(u, j)) monomial_list = monomial_list.union(safe_reducers[-1].monomials()[1:]) break # for polynomials of the current signature index, we have to test against # the criteria if not self.is_rewritable(u, j) and not self.is_top_reducible(u*mj, -1): # if the signature of the reductor is too large, we cannot reduce, so we # put it aside for now; otherwise we add it to list of safe reducers if u*sig(j) < sig(k): #safe_reducers.append(u*poly(j)) safe_reducers.append(self.polynomial_multiple(u, j)) monomial_list = monomial_list.union(safe_reducers[-1].monomials()[1:]) break else: unsafe_reducers.append(j) return safe_reducers, unsafe_reducers #@cached_method # do NOT cache this method if you are using F5C: # the interreduction of the basis wreaks havoc with the cache def polynomial_multiple(self, t, j): return t * self.poly(j) def update_reducers(self, k): """ This computes all multiples of degree d of the polynomials currently in the basis and usable for reduction. This technique is of no real use: it slows F5 down dreadfully, and generates far too many intermediate polynomials. """ reducers = self.reducers reducer_list = self.reducer_list sig = self.sig poly = self.poly too_large = [] # reducer_list contains tuples (t, j) corresponding to t*poly(j) which # appears in reducers # 1. increase degree of reducers if necessary print "checking degree" while poly(k).degree() > reducers[0].degree(): R = poly(k).parent() for each in range(len(reducer_list)): (t, j) = reducer_list.pop(0) p = reducers.pop(0) u, ej = sig(j) for x in R.gens(): if not self.is_rewritable(x*t, j): if not self.is_top_reducible(x*t*u, -1): if not (x*t, j) in reducer_list: # this should never happen reducer_list.append((x*u, j)) reducers.append(x*p) #print "added", x, " * poly(", each, ") with signature", x*sig(each) #else: #print each, x*sig(each), "in list" #else: #print each, x*sig(each), "top-reducible" #else: #print each, x*sig(each), "rewritable by", self.find_rewriting(x*u, each), sig(self.find_rewriting(x*u, each)) print "done with degree" # 2. prune rewritables and signatures that are too large i = 0 print "Pruning reducers @ signature", sig(k) while i < len(reducer_list): t, j = reducer_list[i] if self.is_rewritable(t, j): reducer_list.pop(i) reducers.pop(i) #print "pruned", t*sig(j), "because it's now rewritable by", self.find_rewriting(t, j), sig(self.find_rewriting(t, j)) elif t*sig(j) > sig(k): too_large.append((reducer_list.pop(i),reducers.pop(i))) #print "pruned", t*sig(j), "because it's too large" else: #print "kept", t*sig(j) i = i + 1 print "pruned", len(too_large) return too_large def symbolic_preprocessing(self, k, Gcurr): """ Add polynomials to the set S such that all possible reductors for all elements in S are available. This takes too long and generates far too many polynomials, which strongly implies that I goofed up somewhere. It is meant to be the same as find_reducers, which performs reasonably well, but is superseded by alt_reduction and alt_find_reducers, which seem to accomplish what we want. INPUT: k -- the index of an S-polynomials Gcurr -- the Gröbner basis computed so far indexed in L """ poly = self.poly L = self.L find_reductor = self.find_reductor F = [] Done = set() # the set of all monomials M = set([poly(k).monomials()]) while M != Done: m = M.difference(Done).pop() Done.add(m) t, g = self.identify_reducer(m, self.sig(k), Gcurr) if t!=0: F.append(t*g) M = M.union(set([m for f in F for m in f.monomials()])) return F def identify_reducer(self, m, max_sig, Gcurr): r""" Find a reductor $g_i$ for $m$ in $G_{prev}$ and $G_{curr}$ subject to the following contraint. * the leading monomial of $g_i$ divides $m$ * $g_i \in G_{prev}$ is preferred over $g_i \in G_{curr}$ * if $g_i \in G_{curr}$ then * $g_i$ is not rewritable * $g_i$ is not top reducible by $G_prev$ INPUT: m -- a monomial max_sig -- the maximum signature permissible Ccurr -- the Gröbner basis computed so far indexed in L """ 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 Gcurr: if R.monomial_divides(poly(k).lm(),m): t = R.monomial_quotient(m,poly(k).lm()) if t*sig(k) >= max_sig: continue if is_rewritable(t, k): continue if is_top_reducible(t * sig(k)[0], -1): continue return t, poly(k) return 0, -1 def is_top_reducible(self, t, l): """ Note, that this function test traditional top reduction and not top_reduction as implemented in the function with the same name of this class. """ R = t.parent() poly = self.poly for g in 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): """ adapted from Justin Gash (p.57): 'find_rewriting gives us information to be used as an additional criterion for eliminating critical pairs. Proof of this fact is given in section 3.4.3. In short, we could remove all discussion of rules and find_rewriting and F5 would work fine. (But it would work much more slowly.) So we will treat these final four subroutines as a separate module that works in conjunction with F5, but is not an official part of the F5 criteria per se.' """ 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 class F5R(F5): def basis(self, F): """ F5's main routine. Computes a Gröbner basis for F. INPUT: F -- a list of polynomials OUTPUT: G -- a list of polynomials; a Gröbner basis for """ poly = self.poly incremental_basis = self.incremental_basis self.__init__(F) Rules = self.Rules L = self.L m = len(F) F = sorted(F, cmp=compare_by_degree) f0 = F[0] L[0] = (Signature(1, 0), f0*f0.lc()**(-1)) Rules.append(list()) Gprev = set([0]) B = [f0] for i in xrange(1, m): print "Increment", i f = F[i] L.append( (Signature(1,i), f*f.lc()**(-1)) ) Gcurr = incremental_basis(i, B, Gprev) if any(poly(lambd) == 1 for lambd in Gcurr): return set(1) Gprev = Gcurr B = Ideal([poly(l) for l in Gprev]).interreduced_basis() #B = self.interreduce([poly(l) for l in Gprev]) return B class F5C(F5): def basis(self, F): """ F5's main routine. Computes a Gröbner basis for F. INPUT: F -- a list of polynomials OUTPUT: G -- a list of polynomials; a Gröbner basis for """ incremental_basis = self.incremental_basis poly = self.poly self.__init__(F) Rules = self.Rules L = self.L m = len(F) F = sorted(F, cmp=compare_by_degree) f0 = F[0] L[0] = (Signature(1, 0), f0*f0.lc()**(-1)) Rules.append(list()) Gprev = set([0]) B = set([f0]) for i in xrange(1, m): print "Increment", i f = F[i] L.append( (Signature(1,len(L)), f*f.lc()**(-1)) ) Gcurr = incremental_basis(len(L)-1, B, Gprev) if any(poly(lambd) == 1 for lambd in Gcurr): return set(1) B = Ideal([poly(l) for l in Gcurr]).interreduced_basis() #B = self.interreduce([poly(l) for l in Gcurr]) #B = self.minimal_basis([poly(l) for l in Gcurr]) if i != m-1: Gprev = self.setup_reduced_basis(B) return B def setup_reduced_basis(self, B): """ Update the global L and Rules to match the reduced basis B. OUTPUT: Gcurr -- index set for B """ add_rule = self.add_rule L = self.L Rules = self.Rules # we don't want to replace L but modify it L[:] = [(Signature(1,i), f) for i,f in enumerate(B)] Rules[:] = [[] for _ in xrange(len(B))] Gcurr = set() for i,f in enumerate(B): Gcurr.add( i ) # Eder proved that the following was unnecessary #t = B[i].lt() #for j in xrange(i+1, len(B)): # fjlt = B[j].lt() # u = LCM(t, fjlt)//fjlt # add_rule( Signature(u, j), -1 ) return Gcurr class F4F5(F5C): #class F4F5(F5): """ F4-Style F5 Till Steger's calls this F4.5. We don't know how Jean-Charles Faugère calls it. """ def reduction(self, S, B, Gprev, Gcurr): """ INPUT: S -- a list of components of S-polynomials B -- ignored Gprev -- the previous Gröbner basis indexed in L Ccurr -- the Gröbner basis computed so far indexed in L """ L = self.L add_rule = self.add_rule poly = self.poly S = self.symbolic_preprocessing(S, Gprev, Gcurr) St = self.gauss_elimination(S) Ret = [] for k, (s,p,idx) in enumerate(St): if (s,p,idx) == S[k] and idx == -1: continue # ignore unchanged new polynomials if idx >= 0: L[idx] = L[idx][0], p # update p if p != 0: Ret.append(idx) 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, Gprev, Gcurr): """ 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 Gprev -- the previous Gröbner basis indexed in L Ccurr -- the Gröbner basis computed so far indexed in L """ poly = self.poly 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() # the set of all monomials M = set([m for f in F for m in f[1].monomials()]) while M != Done: m = M.difference(Done).pop() Done.add(m) t, g = find_reductor(m, Gprev, Gcurr) 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 find_reductor(self, m, Gprev, Gcurr): r""" Find a reductor $g_i$ for $m$ in $G_{prev}$ and $G_{curr}$ subject to the following contraint. is a * the leading monomial of $g_i$ divides $m$ * $g_i \in G_{prev}$ is preferred over $g_i \in G_{curr}$ * if $g_i \in G_{curr}$ then * $g_i$ is not rewritable * $g_i$ is not top reducible by $G_prev$ INPUT: m -- a monomial Gprev -- the previous Gröbner basis indexed in L Ccurr -- the Gröbner basis computed so far indexed in L """ 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 Gprev: if R.monomial_divides(poly(k).lm(),m): return R.monomial_quotient(m,poly(k).lm()), L[k] for k in Gcurr: if R.monomial_divides(poly(k).lm(),m): t = R.monomial_quotient(m,poly(k).lm()) if is_rewritable(t, k): continue if is_top_reducible(t * sig(k)[0], Gprev): continue return t, L[k] return 0, -1 def gauss_elimination(self, F1): """ Perform permuted 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(0,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 F = (A*v).list() return [(F1[i][0],F[i],F1[i][2]) for i in xrange(len(F))] class F5SansRewriting(F5): """ A variant of F5 which does not use the rewriting rule. This is motivated by the following observation by Justin Gash (p.57): 'Rewritten gives us information to be used as an additional criterion for eliminating critical pairs. Proof of this fact is given in section 3.4.3. In short, we could remove all discussion of rules and Rewritten and F5 would work fine. (But it would work much more slowly.) So we will treat these final four subroutines as a separate module that works in conjunction with F5, but is not an official part of the F5 criteria per se.' """ def critical_pair(self, k, l, i, Gprev): """ 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_0*F_{e_0} is not top-reducible by (2) S(l) = m_1*F_{e_1} is not top-reducible by (3) S(m_0*k) > S(m_1*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 tk = poly(k).lt() tl = poly(l).lt() t = LCM(tk, tl) u0 = t//tk u1 = t//tl m0, e0 = sig(k) m1, e1 = sig(l) if e0 == e1 and u0*m0 == u1*m1: return set() # Stegers and Gash leave out the == i check, Faugere and Perry # have it. It is unclear for now, whether the check is # necessary. if e0 == i and is_top_reducible(u0*m0, Gprev): return set() if e1 == i and is_top_reducible(u1*m1, Gprev): return set() if u0 * sig(k) < u1 * sig(l): u0, u1 = u1, u0 k, l = l, k return set([(t,k,u0,l,u1)]) def compute_spols(self, P): """ adapted from Justin Gash (p.51): 'Though at first glance this subroutine may look complicated, compute_spols essentially does one thing: form the new S-polynomials output from critical_pairs as admissible signed polynomials. We note that, because critical_pairs ensured that S(u*k) < S(v*l), we know that the signature of all new polynomials will always be of the form u_L*S(r_{i_L}) in compute_spols.' """ poly = self.poly sig = self.sig spol = self.spol L = self.L S = list() P = sorted(P, key=lambda x: x[0]) for (t,k,u,l,v) in P: s = spol(poly(k), poly(l)) if s != 0: L.append( (u * sig(k), s) ) S.append(len(L)-1) S = sorted(S, key=lambda x: sig(x)) return S def top_reduction(self, k, Gprev, Gcurr): """ adapted from Justin Gash (p.55ff): 'We will go through top_reduction step-by-step. If the signed polynomial being examined has polynomial part 0, then there is no data left in that particular signed polynomial - an empty ordered pair is returned. Otherwise top_reduction calls upon another sub-subroutine find_reductor. Essentially, if find_reductor comes back negative, the current signed polynomial is made monic and returned to reduction to be placed in completed. If a top-reduction is deemed possible, then there are two possible cases: either the reduction will increase the signature of polynomial or it won't. In the latter case, the signature of r_{k_0} is maintained, the polynomial portion is top-reduced and the signed polynomial is returned to reduction to be added back into todo; this case corresponds to top-reduction in previous algorithms. In the former case, however, the signature will change. This is marked by adding a new polynomial r_N (our notation here describes N after N was incremented) with appropriate signature based upon the reductor, not S(r_{k_0}). A new rule is added (as I mentioned previously, this will be explained later) and then both r_{k_0} and r_N are sent back to reduction to be added back into todo. This is done because r_N has a different signature than r_{k_0} and r_{k_0} might still be reducible by another signed polynomial. """ from sage.misc.misc import verbose find_reductor = self.find_reductor poly = self.poly sig = self.sig L = self.L if poly(k) == 0: verbose("reduction to zero.",level=0) self.zero_reductions += 1 return set(),set() p = poly(k) J = find_reductor(k, Gprev, Gcurr) if J == set(): L[k] = ( sig(k), p * p.lc()**(-1) ) return set([k]),set() j = J.pop() q = poly(j) u = p.lt()//q.lt() p = p - u*q if p != 0: p = p * p.lc()**(-1) if u * sig(j) < sig(k): L[k] = (sig(k), p) return set(), set([k]) else: L.append((u * sig(j), p)) return set(), set([k, len(L)-1]) def find_reductor(self, k, Gprev, Gcurr): """ adapted from Justin Gash (p.56ff): 'For a previously added signed polynomial in G_curr to become a reductor of r_{k_0}, it must meet four requirements: (1) u = HT(r_{k_0})/HT(r_{k_j}) \in T (2) NF(u_{t_j}, G_curr) = u_{t_j} ... (4) u_{t_j} F_{k_j} = S(r_{k_0}) We will go through each requirement one-by-one. Requirement (1) is simply the normal top-reduction requirement. The only thing of note here is that, in testing for the top-reducibility, u is assigned a particular value to be used in subsequent tests. 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 re- ductor wound 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.) ... Requirement (4) is a check that makes sure we don't reduce by something that has the same signature as r_{k_0} . 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.) """ is_top_reducible = self.is_top_reducible poly = self.poly sig = self.sig t = poly(k).lt() for j in Gcurr: tprime = poly(j).lt() if divides(tprime,t): u = t // tprime mj, ej = sig(j) if u * sig(j) != sig(k) and not is_top_reducible(u*mj, Gprev): return set([j]) return set() 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 = F5() f5r = F5R() f5c = F5C() ff = F4F5()