Source Code for Module basic_sigbased_gb

  1  # Implementations of algorithms found in joint paper by Eder and Perry 
  2  # Copyright (C) 2010-2011, the University of Southern Mississippi 
  3  # released into the public domain 
  4  # 
  5  # Originally from http://www.math.usm.edu/perry/Research/basic_sigbased_gb.py 
  6  # slightly changed for JAS compatibility, changes are labled with JAS 
  7  # $Id: basic_sigbased_gb.py 5474 2016-03-25 16:39:42Z kredel $ 
  8   
  9  # this implementation has one significant difference from the paper: 
 10  # the paper maintains monic signatures, but 
 11  # this implementation maintains monic polynomials 
 12   
 13 -class sigbased_gb: 
 14    # the base class from which all other classes are derived 
 15     
 16 -  def basis_sig(self,F): 
 17      # incremental basis computation 
 18      # F is a container of generators of an ideal 
 19      F.sort(key=lambda x: -x.lm().degree()) # JAS 
 20      #F.sort(key=lambda x: -x.lm().totalDeg()) # JAS 
 21      #print "JAS F = " + str([ str(g) for g in F]);  
 22      G = list() 
 23      for f in F: 
 24        ff = f.reduce(G); 
 25        if (ff == 0): 
 26            continue; 
 27        G = self.incremental_basis(G,f) 
 28        Gnew = [g[1] for g in G] 
 29        R = f.parent() 
 30        print "size before reduction", len(G) 
 31        G = R.ideal(Gnew).interreduced_basis() 
 32      return G 
 33     
 34 -  def spoly_multipliers(self,f,g): 
 35      # multipliers for the s-polynomial of f and g 
 36      # returns uf,ug such that 
 37      # that is, spoly(f,g) = uf.f - ug.g 
 38      tf = f.lm(); tg = g.lm() 
 39      tfg = tf.lcm(tg) 
 40      R = self.R 
 41      return (R.monomial_quotient(tfg,tf),R.monomial_quotient(tfg,tg)) 
 42     
 43 -  def subset(self,S,criterion): 
 44      # this should be changed to use Python's filter() command 
 45      result = set() 
 46      for s in S: 
 47        if criterion(s): 
 48          result.add(s) 
 49      return result 
 50     
 51 -  def min_sig_degree(self,P): 
 52      # determines the minimal degree of a signature in P 
 53      return min([p[0].degree() for p in P]) 
 54     
 55 -  def new_pair(self,sig,p,q,G): 
 56      # creates a new critical pair from p and q, with signature sig 
 57      # it needs G for the sake of F5; see derived class below 
 58      return (sig,p,q) 
 59     
 60 -  def spoly(self,s,G): 
 61      # computes the spolynomial 
 62      # assumes that s has the form (signature, poly, poly) 
 63      f = s[1]; g = s[2] 
 64      tf = f.lm(); tg = g.lm() 
 65      tfg = tf.lcm(tg) 
 66      R = f.parent() 
 67      uf = R.monomial_quotient(tfg,tf); ug = R.monomial_quotient(tfg,tg) 
 68      return uf*f - ug*g 
 69     
 70 -  def initialize_Syz(self,F,G): 
 71      # initializes Syz; initially, this does nothing 
 72      return set() 
 73     
 74 -  def prune_P(self,P,Syz): 
 75      # prunes P using Syz; initially, this does nothing 
 76      return P 
 77       
 78 -  def prune_S(self,S,Syz,Done,G): 
 79      # prunes S using Syz, Done, and G; initially, this does nothing 
 80      # (Done is used as a shortcut) 
 81      return S 
 82     
 83 -  def update_Syz(self,Syz,sigma,r): 
 84      # updates Syz using sigma and r 
 85      # some algorithms use this; some don't 
 86      return Syz 
 87     
 88 -  def sigsafe_reduction(self,s,sigma,G,F,Syz): 
 89      # computes a complete sigma-reduction of s modulo G 
 90      # F is assumed to be a subset of G that represents the previous GB incrementally 
 91      # Syz is sent but not used (I should probably remove this) 
 92      r = s 
 93      r_sigma = sigma 
 94      R = self.R 
 95      reduced = True 
 96      while (r != 0) and reduced: 
 97        reduced = False 
 98        r = r.reduce(F) 
 99        if any(g[1] != 0 and R.monomial_divides(g[1].lm(),r.lm()) for g in G): 
100          for g in G: 
101            if g[1] != 0 and R.monomial_divides(g[1].lm(),r.lm()): 
102              u = self.R.monomial_quotient(r.lt(),g[1].lt(),coeff=True) 
103              sig_ug = u*g[0] 
104              if (sig_ug < r_sigma) or ((sig_ug.lm() == r_sigma.lm()) and (sig_ug.lc() != r_sigma.lc())): 
105                reduced = True 
106                r -= u*g[1] 
107                if (sig_ug.lm() == r_sigma.lm()): 
108                  r_sigma -= sig_ug 
109      # ensure that r is monic 
110      if r != 0: 
111        c = r.lc() 
112        r *= c**(-1) 
113        r_sigma *= c**(-1) 
114      return r_sigma, r 
115       
116 -  def sig_redundant(self,sigma,r,G): 
117      # test whether (sigma,r) is signature-redundant wrt G 
118      R = self.R 
119      return any(g[0] != 0 and R.monomial_divides(g[0].lm(),sigma.lm()) and R.monomial_quotient(sigma.lm(),g[0].lm())*g[1].lm()==r.lm() for g in G) 
120      #if any ((g[0] != 0 and R.monomial_divides(g[0].lm(),sigma.lm())) and (g[1] != 0 and R.monomial_divides(g[1].lm(),r.lm())) and not R.monomial_quotient(sigma.lm(),g[0].lm())*g[1].lm()==r.lm() for g in G): 
121      #  print "counterexample at", (sigma, r.lm()) 
122      return any ((g[0] != 0 and R.monomial_divides(g[0].lm(),sigma.lm())) and (g[1] != 0 and R.monomial_divides(g[1].lm(),r.lm())) for g in G) 
123     
124 -  def incremental_basis(self,F,g): 
125      # assuming that F is a Groebner basis of the ideal generated by F, 
126      # compute a Groebner basis of F+[g] 
127      self.R = g.parent(); R = self.R 
128      # to record a signature, we use only the leading monomial of a minimal representation 
129      # so that elements of F have "signature" 0 and g has "signature" 1 
130      G = [(R(0),F[i]) for i in xrange(len(F))] + [(R(1),g)] 
131      #print "JAS G = " + str([ str(gg[0])+","+str(gg[1]) for gg in G]);  
132      # the structure of a pair can vary, except for its first entry, 
133      # which should be the signature 
134      P = set([self.new_pair(self.spoly_multipliers(g,f)[0],g,f,G) for f in F]) 
135      Syz = self.initialize_Syz(F,G) 
136      # Done will track new polynomials computed by the algorithm 
137      Done = list() 
138      while len(P) != 0: 
139        P = self.prune_P(P,Syz) 
140        if len(P) != 0: 
141          S = list(self.subset(P,lambda x: x[0].degree() == self.min_sig_degree(P))) 
142          print "treating", len(S), "signatures of degree", self.min_sig_degree(P) 
143          P.difference_update(S) 
144          while len(S) != 0: 
145            S = self.prune_S(S,Syz,Done,G) 
146            if len(S) != 0: 
147              # sort by signature 
148              S.sort(key=lambda x:x[0]); s = S.pop(0) 
149              sigma,r = self.sigsafe_reduction(self.spoly(s,G),s[0],G,F,Syz) 
150              if (r != 0) and (not self.sig_redundant(sigma,r,G)):  
151                #print "new polynomial", str(sigma), ", ", str(r) 
152                for (tau,g) in G: 
153                  if (g != 0): 
154                    rmul,gmul = self.spoly_multipliers(r,g) 
155                    if rmul*sigma.lm() != gmul*tau.lm(): 
156                      if rmul*sigma.lm() > gmul*tau.lm(): 
157                        p = self.new_pair(rmul*sigma,r,g,G) 
158                      else: 
159                        p = self.new_pair(gmul*tau,g,r,G) 
160                      if p[0].degree() == sigma.degree(): 
161                        S.append(p) 
162                      else: 
163                        P.add(p) 
164                G.append((sigma,r)) 
165                Done.append((sigma,r)) 
166              elif r == 0: 
167                #print "zero reduction at", (sigma,r.lm())  
168                self.update_Syz(Syz,sigma,r) 
169                Done.append((sigma,r)) 
170              #else: 
171                #print "sig-redundant at", sigma 
172      return list(self.subset(G,lambda x: x[1] != 0)) 
173   
174 -class ggv(sigbased_gb): 
175    # the plugin implementation of ggv 
176     
177 -  def new_pair(self,sig,p,q,G): 
178      # creates a new critical pair from p and q, with signature sig 
179      # it needs G for the sake of F5; see derived class below 
180      i = -1; j = -1; k = 0 
181      up,uq = self.spoly_multipliers(p,q) 
182      while (i<0 or j<0) and k < len(G): 
183        if p == G[k][1]: 
184          i = k 
185        elif q == G[k][1]: 
186          j = k 
187        k += 1; 
188      if (i == -1): 
189        i=len(G) 
190      elif (j == -1): 
191        j = len(G) 
192      return (sig,i,j) 
193     
194 -  def initialize_Syz(self,F,G): 
195      # recognize trivial syzygies 
196      return set([f.lm() for f in F]) 
197       
198 -  def spoly(self,s,G): 
199      # ggv only computes part of an S-polynomial 
200      # (as if it were computing a row of the Macaulay matrix 
201      # and not subsequently triangularizing) 
202      f = G[s[1]][1]; g = G[s[2]][1] 
203      tf = f.lm(); tg = g.lm() 
204      tfg = tf.lcm(tg) 
205      uf = self.R.monomial_quotient(tfg,tf) 
206      return uf*f 
207     
208 -  def prune_P(self,P,Syz): 
209      # remove any pair whose signature is divisible by an element of Syz 
210      result = set() 
211      R = self.R 
212      for p in P: 
213        if not any(R.monomial_divides(t,p[0]) for t in Syz): 
214          result.add(p) 
215      return result 
216     
217 -  def prune_S(self,S,Syz,Done,G): 
218      # watch out for new syzygies discovered, and allow only one polynomial 
219      # per signature 
220      result = list() 
221      R = self.R 
222      for s in S: 
223        if not any(R.monomial_divides(t,s[0]) for t in Syz): 
224          if not any(s[0].lm()==sig[0].lm() and s[1]<sig[1] for sig in S): 
225            if not any(s[0].lm()==sig[0].lm() for sig in result): 
226              result.append(s) 
227      return result 
228     
229 -  def update_Syz(self,Syz,sigma,r): 
230      # add non-trivial syzygies to the basis 
231      # polynomials that reduce to zero indicate non-trivial syzygies 
232      if r == 0: 
233        Syz.add(sigma.lm()) 
234      return Syz 
235   
236 -class ggv_first_implementation(ggv): 
237     
238 -  def new_pair(self,sig,p,q,G): 
239      # creates a new critical pair from p and q, with signature sig 
240      # it needs G for the sake of F5; see derived class below 
241      i = -1; j = -1; k = 0 
242      return (sig,p,q) 
243     
244 -  def spoly(self,s,G): 
245      # ggv only computes part of an S-polynomial 
246      # (as if it were computing a row of the Macaulay matrix 
247      # and not subsequently triangularizing) 
248      # -- at least, that's how I read "(t_i,m+1)" on p. 6 of that paper 
249      f = s[1]; g = s[2] 
250      tf = f.lm(); tg = g.lm() 
251      tfg = tf.lcm(tg) 
252      uf = self.R.monomial_quotient(tfg,tf) 
253      return uf*f 
254   
255 -  def prune_S(self,S,Syz,Done,G): 
256      # watch out for new syzygies discovered, and allow only one polynomial 
257      # per signature 
258      result = list() 
259      R = self.R 
260      for s in S: 
261        if not any(R.monomial_divides(t,s[0]) for t in Syz): 
262          if not any(s[0].lm()==sig[0].lm() for sig in Done): 
263            result.append(s) 
264      return result 
265     
266 -class coeff_free_sigbased_gb(sigbased_gb): 
267    # child class of sigbased_gb that implements semi-complete reduction 
268     
269 -  def sigsafe_reduction(self,s,sigma,G,F,Syz): 
270      # see sigbased_gb.sigsafe_reduction 
271      r = s 
272      r_sigma = sigma 
273      R = self.R 
274      reduced = True 
275      while (r != 0) and reduced: 
276        reduced = False 
277        r = r.reduce(F) 
278        if any( g[1] != 0 and R.monomial_divides(g[1].lm(),r.lm()) for g in G ): 
279          for g in G: 
280            if g[1] != 0 and R.monomial_divides(g[1].lm(),r.lm()): 
281              u = self.R.monomial_quotient(r.lt(),g[1].lt(),coeff=True) 
282              sig_ug = u*g[0] 
283              if sig_ug.lm() < r_sigma.lm(): # type error: sig_ug.lm() < r_sigma  
284                reduced = True 
285                r -= u*g[1] 
286      # ensure that r is monic 
287      if r != 0: 
288        c = r.lc() 
289        r *= c**(-1) 
290      return r_sigma, r 
291       
292 -class arris_algorithm(coeff_free_sigbased_gb): 
293    # the plugin implementation of arri's algorithm 
294     
295 -  def initialize_Syz(self,F,G): 
296      # recognize trivial syzygies 
297      return set([f.lm() for f in F]) 
298       
299 -  def update_Syz(self,Syz,sigma,r): 
300      # add non-trivial syzygies to the basis 
301      # polynomials that reduce to zero indicate non-trivial syzygies 
302      if r == 0: 
303        Syz.add(sigma.lm()) 
304      return Syz 
305     
306 -  def prune_P(self,P,Syz): 
307      # remove any pair whose signature is divisible by an element of Syz 
308      result = set() 
309      for p in P: 
310        if not any(self.R.monomial_divides(t,p[0]) for t in Syz): 
311          result.add(p) 
312      return result 
313     
314 -  def prune_S(self,S,Syz,Done,G): 
315      # watch out for new syzygies discovered, and apply arri's rewritable criterion: 
316      # for any s-polynomial of a given signature, if there exists another (s-)polynomial 
317      # in S or Done of identical signature but lower lm, discard the first 
318      result = list() 
319      for s in S: 
320        if not any(self.R.monomial_divides(t,s[0]) for t in Syz): 
321          if not any(s[0]==sig[0] and s[1].lm()>sig[1].lm() for sig in S): 
322            for (sig,f) in Done: 
323              if self.R.monomial_divides(sig,s[0]): 
324                u = self.R.monomial_quotient(s[0],sig) 
325                if u*f.lm() < s[1].lm(): 
326                  break 
327            else: 
328              result.append(s) 
329      return result 
330     
331 -  def new_pair(self,sig,p,q,G): 
332      # in arri's algorithm, each pair is (sigma,s) where s is the s-polynomial 
333      # and sigma is its natural signature 
334      tp = p.lm(); tq = q.lm() 
335      tpq = tp.lcm(tq) 
336      R = p.parent() #JAS tpq.parent() 
337      up = R.monomial_quotient(tpq,tp); uq = R.monomial_quotient(tpq,tq) 
338      return (sig,up*p-uq*q) 
339     
340 -  def spoly(self,s,G): 
341      return s[1] 
342   
343 -class f5(coeff_free_sigbased_gb): 
344    # the plugin implementation of arri's algorithm 
345     
346 -  def initialize_Syz(self,F,G): 
347      # recognize trivial syzygies 
348      return set([f.lm() for f in F]) 
349       
350 -  def update_Syz(self,Syz,sigma,r): 
351      # recognize trivial syzygies 
352      # see class f5z for a more thorough update_Syz in line w/arri and ggv 
353      return Syz 
354     
355 -  def prune_P(self,P,Syz): 
356      # remove any pair whose signature is divisible by an element of Syz 
357      result = set() 
358      for p in P: 
359        if not any(self.R.monomial_divides(t,p[0]) for t in Syz): 
360          result.add(p) 
361      return result 
362     
363 -  def prune_S(self,S,Syz,Done,G): 
364      # watch out for new syzygies discovered, and apply faugere's rewritable criterion: 
365      # for any (sigma,p,q) in S, if there exists (tau,g) such that tau divides sigma 
366      # but g was generated after p, discard (sigma,p,q) 
367      result = list() 
368      for (sig,u,j,v,k) in S: 
369        if not any(self.R.monomial_divides(t,sig) for t in Syz): 
370          if G[j][0] == 0 or not any(self.R.monomial_divides(Done[i][0],G[j][0]*u) and Done[i][0] > G[j][0] for i in xrange(len(Done))): 
371            result.append((sig,u,j,v,k)) 
372      return result 
373     
374 -  def new_pair(self,sig,p,q,G): 
375      # it's easier to deal with faugere's criterion if one creates pairs 
376      # using indices rather than polynomials 
377      # note that this while look gives f5 a disadvantage 
378      i = -1; j = -1; k = 0 
379      up,uq = self.spoly_multipliers(p,q) 
380      while (i<0 or j<0) and k < len(G): 
381        if p == G[k][1]: 
382          i = k 
383        elif q == G[k][1]: 
384          j = k 
385        k += 1; 
386      if (i == -1): 
387        i=len(G) 
388      elif (j == -1): 
389        j = len(G) 
390      return (sig,up,i,uq,j) 
391       
392 -  def spoly(self,s,G): 
393      # since s has the structure (sigma,up,i,uq,j) 
394      # we have to compute the s-polynomial by looking up f and g 
395      f = G[s[2]][1]; g = G[s[4]][1] 
396      uf = s[1]; ug = s[3] 
397      return uf*f - ug*g 
398       
399 -class f5z(f5): 
400     
401 -  def update_Syz(self,Syz,sigma,r): 
402      # recognize trivial syzygies 
403      if r == 0: 
404        Syz.add(sigma.lm()) 
405      return Syz 
406       
407 -class min_size_mons(arris_algorithm): 
408    # the plugin implementation of arri's algorithm 
409     
410 -  def prune_S(self,S,Syz,Done,G): 
411      # watch out for new syzygies discovered, and apply the minimal "number of monomials" criterion: 
412      # for any s-polynomial of a given signature, if there exists another polynomial 
413      # in Done of identical signature but fewer monomials, replace this s-polynomial 
414      # by the multiple of the polynomial with fewer monomials 
415      result = list() 
416      R = G[0][1].parent() 
417      for (sigma,s) in S: 
418        if not any(R.monomial_divides(tau,sigma) for tau in Syz): 
419          if not any(tau == sigma for (tau,g) in result): 
420            for (tau,g) in Done: 
421              if tau.divides(sigma) and len(g.monomials()) < len(s.monomials()): 
422                u = R.monomial_quotient(sigma,tau) 
423                result.append((u*tau,u*g)) 
424                break 
425            else: 
426              result.append((sigma,s)) 
427      return result 
428