## Attachment 'f5.pyx'

   1 """
2 Faugere's F_5 algorithm
3
4 This implementation is along the lines of John Perry's pseudocode
5 and Singular implementation. It was inspired by Martin Albrecht's
6 implementation and discussions with him and Ludovic Perret.
7
8 See http://www.math.usm.edu/perry/Research/ for details.
9
10 AUTHOR:
11     -- 20081022 Simon King
12 EXAMPLE:
13     sage: R.<x,y,z> = PolynomialRing(GF(29))
14     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, \
15                    3*x^3*y^2 + 7*x^2*y^3 + 24*y^2*z^3,
16                    12*x*y^4 + 17*x^4*z + 27*y^4*z + 11*x^3*z^2]
17     sage: J = I.homogenize()
18
19     sage: f5 = F5() # original F5
20     sage: gb = f5(J)
21     Generator 2/3
22     1 critical pairs in degree 7
23     1 critical pairs in degree 9
24     1 critical pairs in degree 11
25     Generator 3/3
26     1 critical pairs in degree 6
27     2 critical pairs in degree 7
28     5 critical pairs in degree 8
29     7 critical pairs in degree 9
30     11 critical pairs in degree 10
31     10 critical pairs in degree 11
32     WARNING: top-reduction to zero
33     WARNING: reduction to zero
34     WARNING: reduction to zero
35     3 critical pairs in degree 12
36     5 critical pairs in degree 13
37     2 critical pairs in degree 14
38     1 critical pairs in degree 16
39     Note that we encountered 3 S-polynomials reducing to zero
40     sage: f5.zeroes(), len(gb)
41     (3, 18)
42 """
43
44 import sage
45 import sage.all
46 from sage.all import copy
47 import sys
48 from sage.structure.sequence import Sequence
49 from sage.rings.ideal import Ideal
50 #from sage.rings.polynomial.multi_polynomial_libsingular cimport MPolynomial_libsingular
51
52 class CounterClass:
53     def __init__(self):
54         self.Gen = 0
55         self.Del = 0
56     def __call__(self,i):
57         if i>0:
58             self.Gen += 1
59         else:
60             self.Del += 1
61 Counter = CounterClass()
62 CounterB = CounterClass()
63
64 def dehomogenize(I,R):
65     "dehomogenize(I,R): dehomogenize an ideal I that was obtained by homogenizing an ideal of ring R"
66     h = I.parent()('h')
67     return ([R(X(h=1)) for X in I]*R).reduced_basis()
68
69 cdef class DecoratedPolynomial:
70     """
71     A decorated polynomial is a polynomial p in an ideal I = <f_1,...,f_m> of a polynomial
72     ring R, together with a 'signature' (u,i), where u is a monomial in R and 1 <= i <= m
73     is an integer. The signature has the following meaning:
74     * One can express p as an ideal combination $\sum_{k=1}^i t_k\cdot f_k$, and u is the leading
75       term of $t_i$ of f_1,...,f_i.
76     The signature allows for avoiding some useless critical pairs in the Buchberger algorithm.
77     """
78     cdef int nu
79     cdef object mu  # the multiplier
80     cdef tuple Emu  # exponents of the multiplier
81     #cdef MPolynomial_libsingular p
82     cdef object p
83
84     def __init__(self,mu, int nu,  p):
85         self.mu = mu
86         self.Emu = mu.degrees()
87         self.nu = nu
88         self.p = p
89         #self.rule = -1
90         #self.hist = []
91     #def __init__(self,mu, int nu, object p):
92     #    pass
93     def sig(DecoratedPolynomial self):
94         return (self.mu,self.nu)
95     def poly(DecoratedPolynomial self):
96         return self.p
97     #def info(self):
98     #    return self.hist
99
100     def __richcmp__(DecoratedPolynomial self, DecoratedPolynomial other, op):
101         # < 0, <= 1, == 2, != 3, > 4, >= 5
102         # Idea for applications: When considering a pair of decorated monomials and
103         # reducing them, we always reduce the "bigger" one,
104         if op == 0:
105             if self.nu<other.nu:
106                 return True
107             if (self.nu==other.nu) and (self.mu<other.mu):
108                 return True
109             return False
110         if op == 1:
111             if self.nu<other.nu:
112                 return True
113             if (self.nu==other.nu) and (self.mu<=other.mu):
114                 return True
115             return False
116         if op == 2:
117             return (self.nu==other.nu) and (self.Emu==other.Emu) and (self.p==other.p)
118         if op == 3:
119             return not ((self.nu==other.nu) and (self.Emu==other.Emu) and (self.p==other.p))
120         if op == 4:
121             if self.nu>other.nu:
122                 return True
123             if (self.nu==other.nu) and (self.mu>other.mu):
124                 return True
125             return False
126         if op == 5:
127             if self.nu>other.nu:
128                 return True
129             if (self.nu==other.nu) and (self.mu>=other.mu):
130                 return True
131             return False
132
133
134 cdef class F5:
135     cdef dict Rules       # rewriting rules (indexed by leading monomial, yielding a monomial ideal)
136     cdef object Rgb       # reduced GB of the first few input generators
137     cdef object RgbMon    # their leading ideal
138     cdef list Inputgen    # the ideal generators we got, interreduced and increasingly sorted
139     cdef list Basis       # list of decorated polynomials that yields self.Rgb
140     cdef object Ring      # the polynomial ring we started with
141     cdef object HRing     # homogenized ring
142     cdef object One       # "One" in the ring we work with (may be homogenisation of self.Ring
143     cdef int nvars        # number of variables of HRing
144     cdef int homog        # 1, if the input was homogeneous; 0 otherwise (then we homogenize/dehomogenize)
145     cdef int Zeroes       # Number of S-polynomials reducing to zero, to test the efficiency of
146                           # our implementation of the F5 criterion.
147                           # We are not keeping track of the corresponding critical pairs.
148     def __init__(F5 self):
149         "F = F5() set up the machinery to compute Groebner bases with the F5 algorithm"
150         self.Rules={}
151         self.Rgb = None
152         self.RgbMon = None
153         self.Inputgen = []
154         self.Basis = []
155         self.Ring = None
156         self.HRing = None
157         self.One = None
158         self.nvars = 0
159         self.homog = 0
160         self.Zeroes = 0
161
162     #def __init__(self):
163     #    pass
164
165     #def __dealloc__(self):
166     #    self.Rules={}
167     #    self.Rgb = None
168     #    self.RgbMon = None
169     #    self.Inputgen = []
170     #    self.Basis = []
171     #    self.Ring = None
172     #    self.HRing = None
173     #    self.One = None
174     #    self.nvars = 0
175     #    self.homog = 0
176     #    self.Zeroes = 0
177     #    print "Deleted"
178
179     def __call__(F5 self, I):
180         """
181         F(I): Apply the F5 algorithm to an ideal I and return its reduced Groebner basis
182         """
183         if not isinstance(I, sage.rings.polynomial.multi_polynomial_ideal.MPolynomialIdeal):
184             I = Ideal(I)
185         self.Ring = I.ring()
186         self.Zeroes = 0 # S-polynomials reducing to zero.
187                          # We do not keep track of the corresponding critical pair.
188         self.Rgb            = Sequence([I.ring()(0)],I.ring()) # (reduced) GB of the first few input generators
189         self.RgbMon = self.Rgb
190         # TODO: is homogenization really needed?
191         # Certainly it is needed in the F5 matrix version. But here??
192         if not I.is_homogeneous():
193             print "homogenizing the input"
194             J = I.homogenize()
195             self.homog = 0
196         else:
197             J = I
198             self.homog = 1
199         if J.gens()==(0,):
200             self.Inputgen = [J.ring()(0)]     # reduced_basis would yield a segfault
201         else:
202             self.Inputgen       = list(J.reduced_basis())   # the ideal generators we got, sorted increasingly
203         self.HRing = J.ring()
204         self.Inputgen.sort()
205         self.Rules          = {} # rewriting rules
206         self.Basis          = [] # List of decorated polynomials that will eventually form a
207                                  # a Groebner basis of J
208         self.One            = J.ring()(1) # self.One.parent() =self.HRing != self.Ring, if there was a homogenisation
209         self.nvars          = len(J.ring().gens())
210         return self.basis()
211
212     ## Some methods that reveal the C-attributes of self.
213     def rules(self):
214         return self.Rules
215     def zeroes(self):
216         return self.Zeroes
217
218     #############
219     ## Main procedure
220     ## Return the previously computed basis, or do the computation when first called
221     def basis(slf, maxgen=None):
222         """
223         F.basis(): Compute a reduced Groebner basis using the F5 algorithm.
224         """
225         cdef F5 self = slf
226         if self.Basis:
227             if self.homog:
228                 return self.Rgb # or [X.poly for X in self.Basis]*self.HRing, currently without reduction
229             else:
230                 return dehomogenize(self.Rgb,self.Ring)
231         if not self.homog:
232             h = self.HRing('h')
233         if self.Inputgen == 0:
234             return Sequence([self.HRing(0)], self.HRing)
235         cdef DecoratedPolynomial newDP = DecoratedPolynomial(self.One,0,self.Inputgen/self.Inputgen.lc())
236
237         self.Rules=[]
238         self.Basis.append(newDP)
239         cdef int i
240         cdef int UpTo
241         if maxgen is None:
242             UpTo = len(self.Inputgen)
243         else:
244             UpTo = maxgen
245         cdef int laststop = 0
246         cdef DecoratedPolynomial X,Y
247         for i in range(1,UpTo):
248             print "Generator %d/%d"%(i+1,UpTo)
249             # by induction, we have a reduced GB for the generators 1,...,i-1.
250             self.Rgb = ([Y.p for Y in self.Basis]*self.HRing).reduced_basis()
251             self.RgbMon = ([Y.p.lm() for Y in self.Basis]*self.HRing).reduced_basis()
252             self.insert(i)
253             for X in self.Basis[laststop:]:
254                 # it it contains a constant polynomial, the ideal is trivial
255                 if self.homog:
256                     if X.p.degree()==0:
257                         self.Rgb = Sequence([self.One], self.HRing)
258                         return Sequence([self.Ring(1)], self.Ring)
259                 else:
260                     if X.p(h=1).degree()==0:
261                         self.Rgb = Sequence([self.One], self.HRing)
262                         return Sequence([self.Ring(1)], self.Ring)
263             laststop = len(self.Basis)
264         # Now, there only remains to toss out the reduced GB
265         if self.Zeroes:
266             print "Note that we encountered %d S-polynomials reducing to zero"%(self.Zeroes)
267         self.Rgb = ([Y.p for Y in self.Basis]*self.HRing).reduced_basis()
268         if self.homog:
269             return self.Rgb
270         else:
271             return dehomogenize(self.Rgb,self.Ring)
272
273     cdef insert(F5 self,int i):
274         """F.insert(i): Compute a Groebner basis for the ideal spanned by generators #0,...,i,
275            when a Groebner basis for generators #0,...,i-1 is already known.
276         """
277         # We may reduce the i-th generator f_i by everything that is in the "old" ideal (f_0,...,f_{i-1})
278         f = self.Inputgen[i].reduce(self.Rgb)
279         self.Rules[i]=[]
280         if f.lc()==0:
281             return
282         cdef DecoratedPolynomial newDP = DecoratedPolynomial(self.One, i, f/f.lc())
283
284         cdef int newindices = len(self.Basis)  # self.Basis[:newindices] was our starting point
285         self.Basis.append(newDP)
286         cdef int j,d
287         cdef list P
288         cdef list Pd,tmpL,X
289         cdef DecoratedPolynomial Y,Z
290         # generate new S-polynomials
291         # The argument i is used for the F5 criterion (it appears in the signature of
292         # polynomials we want to reduce with)
293         P=self.crit_pairs(newDP,newindices,i,newindices)  # only crit pairs involving newDP and old GB elements
294         if not self.homog:
295             h = self.HRing('h')
296         while P:
297             #print i
298             #print "len P =",len(P)
299             d = min([X.degree() for X in P])
300             tmpL,Pd = [],[]
301             for X in P:
302                 if X.degree() == d:
303                     Pd.append(X)
304                 else:
305                     tmpL.append(X)
306             print "%d critical pairs in degree %d"%(len(Pd),d)
307             P = tmpL
308             S = self.SPolys(Pd)
309             R = self.reduction(S, i, newindices)
310             for Y in R:
311                 self.Basis.append(Y)
312                 if self.homog:
313                     if Y.p.degree()==0:
314                         return
315                 else:
316                     if Y.p(h=1).degree()==0:
317                         return
318                 P.extend(self.crit_pairs(Y,len(self.Basis),i,newindices))   # crit pairs for all of self.Basis (including the new ones)
319
320
321     cdef list crit_pair(self, DecoratedPolynomial X, DecoratedPolynomial Y, int i, int newindices):
322         """
323         F.crit_pair(X,Y,i,newindices) returns (t,u_X,X,u_Y,Y) corresponding to a critical pair
324             X,Y when necessary for the computation of a Groebner basis of (f_1,...,f_i). Here,
325             the F5 criterion for label i is used. By switching X and Y, if necessary, we will
326             have X<Y in the output.
327             Assumption: F.Basis[:newindices] (a private attribute) corresponds to a GB of
328                         (f_1,...,f_{i-1}).
329         """
330         R = self.HRing
331         t_X = X.p.lm()
332         t_Y = Y.p.lm()
333         t = t_X.lcm(t_Y)
334         u_X = R.monomial_quotient(t,t_X)
335         u_Y = R.monomial_quotient(t,t_Y)
336         umuX = X.mu*u_X
337         if X.nu==i and self.topreducible(umuX):
338             #X.hist.append(('topred',umuX,i))
339             return [] # since the row with signature (umuX,nu) is kicked out by F5, we don't need the S-poly
340         umuY = Y.mu*u_Y
341         if (X.nu==Y.nu) and (umuX==umuY):  # found in Perry's f5_library.lib
342             # I guess this is part of "rewrite criterion".
343             # But I confess I don't understand why this is correct
344             return []
345         if Y.nu==i and self.topreducible(umuY):
346             #Y.hist.append(('topred',umuY,i))
347             return []
348         # adopt minor optimization (check is_rewritable)
349         if self.is_rewritable(u_X,X) or self.is_rewritable(u_Y,Y):
350             return []
351         # Eventually we try to reduce from "top" to "bottom".
352         # We return the pair X,Y, if u_X*X is on top of u_Y*Y.
353         # Hence, later we will replace u_Y*Y by the S-poly of X and Y
354         if (X.nu>Y.nu) or ((X.nu==Y.nu) and (umuX>umuY)):
355             return [[t,Y,u_Y,X,u_X]]
356         return [[t,X,u_X,Y,u_Y]]
357
358     cdef list crit_pairs(self, DecoratedPolynomial X, int ub, int i, int newindices):
359         "compare with crit_pair, but here Y runs on self.Basis[:ub]"
360         cdef list CP = []
361         cdef DecoratedPolynomial Y
362         for Y in self.Basis[:ub]:
363             CP.extend(self.crit_pair(X,Y,i,newindices))
364         return CP
365
366     cdef list SPolys(self,list P):
367         """
368         F.SPolys(P), P a list of critical pairs, given by a 5-tuple of
369           1) a monomial t (the least common multiple of the leading monomials of X.poly() and Y.poly()),
370           2) the cofactor u_X = t/lm(X.poly()) for X,
371           3) a decorated polynomial X,
372           4) the cofactor u_Y = t/lm(X.poly()) for Y,
373           5) a decorated polynomial Y.
374         Moreover, (u_X*X)<(u_Y*Y).
375         The S-polynomials are computed and tested for rewritability. The non-rewritable S-polynomials
376         are added to the rewrite rules and are returned in a list.
377         """
378         # P is a list of tuples (lcm, u_X,X, u_Y,Y), with u_X*lm(X) = u_Y*lm(Y) = lcm
379         P.sort() # which means increasing by lcm
380         cdef list S = []
381         cdef DecoratedPolynomial X,Y,Z
382         for (t,X,u_X,Y,u_Y) in P:
383             if (not self.is_rewritable(u_X,X)) and (not self.is_rewritable(u_Y,Y)):
384                 s = (u_X*X.p*Y.p.lc() - u_Y*Y.p*X.p.lc()) # it will be reduced soon, in a different method...
385                 if s!=0:
386                     # We know that u_X*X is on top of u_Y*Y.
387                     # Hence, thinking abouz the F5 matrix version, we would
388                     # create rows labeled (X.mu*u_X,X.nu) and (Y.mu*u_Y,Y.nu),
389                     # and a reduction from top to bottom would yield
390                     # the S-polynomial at (Y.mu*u_Y,Y.nu).
391                     # TODO: Is this thinking correct?
392                     Z = DecoratedPolynomial(Y.mu*u_Y,Y.nu, s/s.lc())
393                     #Z.hist.append(('S-poly',X,X.sig(),Y,Y.sig()))
395                     S.append(Z)
396         S.sort() # by increasing signature, hence, from top to bottom in the F5 matrix version
397         return S
398
399     cdef list reduction(self,list S, int i, int newindices):
400         """
401         F.reduction(S,i,newindices): S a list of decorated polynomials, and F.Basis[:newindices] (a private
402         attribute) provides a Groebner basis for the ideal generatd by the generators #0,...,i-1.
403         It is assumed that the list is sorted by increasing signature.
404         We will reduce each element of S by the GB and the elements of S of lower signature.
405         Returns the list of completely reduced decorated polynomials.
406         """
407         cdef DecoratedPolynomial X,Y
408         cdef list truncGB = list(self.Rgb)
409         # We use self.Rgb, since it is a reduced GB for the ideal of the first i-1 indices.
410         # Later on, we add to it new reduced S-polynomials.
411         cdef list completed  = []
412         cdef list newly_completed = []
413         cdef list redo = []
414         #sys.stdout.write('r')
415         #sys.stdout.flush()
416         while S:
417             #sys.stdout.write('.')
418             #sys.stdout.flush()
419             X = S.pop(0) # this is the one with the smallest signature
420             # By construction, the signature of all our input polynomials has index i:
421             #sys.stdout.write('r')
422             #sys.stdout.flush()
423             X.p = X.p.reduce(truncGB)
424             #sys.stdout.write('.')
425             #sys.stdout.flush()
426             #X.hist.append('reduction')
427             if X.p.lc()!=0:
428                 X.p = X.p/X.p.lc()
429             newly_completed, redo = self.top_reduction(X,newindices,completed)
430             completed.extend(newly_completed)
431             truncGB.extend([Y.p for Y in newly_completed]) # hence basis is "old basis" plus "completed"
432             S.extend(redo)
433             S.sort() # by increasing signature, since reductions must be done from top to bottom.
434             #print "S",len(S)
435             #print "rule",len(self.Rules[i])
436             #print "Basis", len(self.Basis)
437         #sys.stdout.write('\n')
438         #sys.stdout.flush()
439         return completed
440
441     cdef tuple top_reduction(self, DecoratedPolynomial X, int newindices,list completed):
442         if X.p.lc()==0:
443             print "WARNING: reduction to zero"
444             self.Zeroes = self.Zeroes+1
445             return ([],[])
446         #X.hist.append(('test',X.p))
447         #sys.stdout.write('t')
448         #sys.stdout.flush()
449         J = self.find_reductor(X,newindices,completed)
450         #sys.stdout.write('.')
451         #sys.stdout.flush()
452         if not J:
453             return ([X],[])
454         cdef DecoratedPolynomial Y = J
455
456         u = J
457         p = X.p-u*Y.p
458         #sys.stdout.write('\n')
459         #sys.stdout.flush()
460         if p.lc()!=0:
461             p = p/p.lc()
462         else:
463             print "WARNING: top-reduction to zero"
464             self.Zeroes = self.Zeroes+1
465             return ([],[])
466         newmono = Y.mu*u
467         # TODO: is it newmono<X.mu or newmono<=X.mu?
468         # At least it seems there is an infinite loop if it is "<" only
469         if (Y.nu<X.nu) or (Y.nu==X.nu and newmono<=X.mu):
470             # i.e., (signature of Y)*u is not bigger, hence, it is on top of sig(X)
471             # hence, X ought to be replaced
472             X.p = p
473             #X.hist.append(('top-red by',Y,'with',u,X.p.lm(),Y.p.lm(),X.p.lm()))
474             #X.hist.append('redo1')
475             return ([],[X])
476         # Otherwise, (signature of Y)*u is below sig(X). Hence,
477         # the corresponding row (that is now being created) must be reduced
478         # (which is done right on the spot).
479         cdef DecoratedPolynomial Z = DecoratedPolynomial(newmono,Y.nu,p)
480
481         #Z.hist.append(('top-red obtained from',X,Y))
482         #X.hist.append(('top-red 2',Y,Z))
484         return ([],[X,Z])
485
486     cdef list find_reductor(self, DecoratedPolynomial X, int newindices, list completed):
487         """
488         F.find_reductor(X,i,C): returns [], or [[u,Y]] if there is a decorated polynomial Y
489         in self.Basis+C (private attribute) and X is *safely* top-reducible by Y with multiplier u.
490         That's to say, u*Y has the same leading monomial as X, Y is neither detected by the F5
491         criterion nor the rewritten criterion, and u*Y has smaller signature than X.
492         """
493         cdef DecoratedPolynomial Y
494         t = X.p.lm()
495         R = self.HRing
496         for Y in self.Basis:
497             if Y<X:
498                 tt = Y.p.lm()
499                 if R.monomial_divides(tt,t):
500                     u = R.monomial_quotient(t,tt)
501                     umu = Y.mu*u
502                     # Y is a reductor of X if u*Y has different signature from X,
503                     # (u,Y) is not rewritable, and umu is not topreducible.
504                     if ((Y.nu!=X.nu) or (umu!=X.mu)) and (not self.topreducible(umu)) and (not self.is_rewritable(u,Y)):
505                         return [u,Y]
506         for Y in completed:
507             if Y<X:
508                 tt = Y.p.lm()
509                 if R.monomial_divides(tt,t):
510                     u = R.monomial_quotient(t,tt)
511                     umu = Y.mu*u
512                     if ((Y.nu!=X.nu) or (umu!=X.mu)) and (not self.topreducible(umu)) and (not self.is_rewritable(u,Y)):
513                         return [u,Y]
514         return []
515
516     def topreducible(self, u):
517         """
518         F.topreducible(u): Tests whether the F5 criterion would apply to decorated polynomials
519         with signature (u,i), where u is a monomial and self.Rgb (private attribute) provides
520         a Groebner basis for the ideal spanned by the generators #0,...,i-1.
521         """
522         return u.reduce(self.RgbMon)==0
523
525         "F.add_rule(X) adds the rewriting rule corresponding to the decorated polynomial X"
527         self.Rules[X.nu].append(X)
528
529     def is_rewritable(self, u, DecoratedPolynomial X):
530         """
531         F.is_rewritable(u,X) tests if the is a rewriting (corresponding to Gaussian
532         elimination "from top") for u*X, where u is a monomial and X a decorated polynomial.
533         """
534         cdef DecoratedPolynomial Y
535         cdef list R = self.Rules[X.nu]
536         cdef int l = len(R)
537         HR = self.HRing
538         umu = X.mu*u
539         while (1):
540             l-=1
541             if l<0:
542                 return False
543             Y = R[l]
544             if Y is X:
545                 return False
546             if HR.monomial_divides(Y.mu,umu):
547                 return True
548         #return (X.mu*u).reduce([Y.mu for Y in self.Rules[X.nu][X.rule+1:]])==0


## Attached Files

To refer to attachments on a page, use attachment:filename, as shown below in the list of files. Do NOT use the URL of the [get] link, since this is subject to change and can break easily.
• [get | view] (2008-10-14 14:37:00, 30.0 KB) [[attachment:[email protected]]]
• [get | view] (2008-10-14 15:39:57, 12.4 KB) [[attachment:f5.py]]
• [get | view] (2008-10-29 13:41:03, 22.7 KB) [[attachment:f5.pyx]]
• [get | view] (2008-10-23 05:36:34, 4.8 KB) [[attachment:gema.sage]]
• [get | view] (2008-10-15 10:01:59, 18.7 KB) [[attachment:jgd_solve.patch]]
• [get | view] (2008-10-15 09:30:20, 31.2 KB) [[attachment:m4ri_trsm_UL_LR.patch]]
• [get | view] (2008-10-16 12:57:59, 10.8 KB) [[attachment:sage-view.el]]
All files | Selected Files: delete move to page copy to page