[frames] | no frames]

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 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):
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.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; g = s
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 != 0 and R.monomial_divides(g.lm(),r.lm()) for g in G):
100          for g in G:
101            if g != 0 and R.monomial_divides(g.lm(),r.lm()):
102              u = self.R.monomial_quotient(r.lt(),g.lt(),coeff=True)
103              sig_ug = u*g
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
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 and R.monomial_divides(g.lm(),sigma.lm()) and R.monomial_quotient(sigma.lm(),g.lm())*g.lm()==r.lm() for g in G)
120      #if any ((g != 0 and R.monomial_divides(g.lm(),sigma.lm())) and (g != 0 and R.monomial_divides(g.lm(),r.lm())) and not R.monomial_quotient(sigma.lm(),g.lm())*g.lm()==r.lm() for g in G):
121      #  print "counterexample at", (sigma, r.lm())
122      return any ((g != 0 and R.monomial_divides(g.lm(),sigma.lm())) and (g != 0 and R.monomial_divides(g.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)+","+str(gg) 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),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.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); s = S.pop(0)
149              sigma,r = self.sigsafe_reduction(self.spoly(s,G),s,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.degree() == sigma.degree():
161                        S.append(p)
162                      else:
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 != 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]:
184          i = k
185        elif q == G[k]:
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]; g = G[s]
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) for t in Syz):
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) for t in Syz):
224          if not any(s.lm()==sig.lm() and s<sig for sig in S):
225            if not any(s.lm()==sig.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:
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; g = s
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) for t in Syz):
262          if not any(s.lm()==sig.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 != 0 and R.monomial_divides(g.lm(),r.lm()) for g in G ):
279          for g in G:
280            if g != 0 and R.monomial_divides(g.lm(),r.lm()):
281              u = self.R.monomial_quotient(r.lt(),g.lt(),coeff=True)
282              sig_ug = u*g
283              if sig_ug.lm() < r_sigma.lm(): # type error: sig_ug.lm() < r_sigma
284                reduced = True
285                r -= u*g
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:
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) for t in Syz):
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) for t in Syz):
321          if not any(s==sig and s.lm()>sig.lm() for sig in S):
322            for (sig,f) in Done:
323              if self.R.monomial_divides(sig,s):
324                u = self.R.monomial_quotient(s,sig)
325                if u*f.lm() < s.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
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) for t in Syz):
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 or not any(self.R.monomial_divides(Done[i],G[j]*u) and Done[i] > G[j] 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]:
382          i = k
383        elif q == G[k]:
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]; g = G[s]
396      uf = s; ug = s
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:
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.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
<!--
expandto(location.href);
// -->

```

 Generated by Epydoc 3.0.1 on Sun Apr 21 14:11:38 2019 http://epydoc.sourceforge.net