## Attachment 'f5.py'

```   1 """
2 Faugere's F5
3
4 These implementations are heavily inspired by John Perry's pseudocode
5 and Singular implementation of these algorithms.
6
7 See http://www.math.usm.edu/perry/Research/ for details.
8
9 This implementation runs faster than the Singular script
10 implementation but uses much more RAM for some reason.
11
12 AUTHOR:
13     -- 20081013 Martin Albrecht (initial version based on John Perry's pseudocode)
14     -- 20081013 John Perry (loop from 0 to m-1 instead of m-1 to 0)
15
16 EXAMPLE:
17     sage: R.<x,y,z> = PolynomialRing(GF(29))
18     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, \
19                    3*x^3*y^2 + 7*x^2*y^3 + 24*y^2*z^3,
20                    12*x*y^4 + 17*x^4*z + 27*y^4*z + 11*x^3*z^2]
21     sage: J = I.homogenize()
22
23     sage: f5 = F5() # original F5
24     sage: gb = f5(J)
25     sage: f5.zero_reductions, len(gb)
26     d 7
27     d 9
28     d 11
29     d 6
30     d 7
31     d 8
32     d 9
33     d 10
34     d 11
35     verbose 0 (179: 283.py, top_reduction) reduction to zero.
36     verbose 0 (179: 283.py, top_reduction) reduction to zero.
37     verbose 0 (179: 283.py, top_reduction) reduction to zero.
38     d 12
39     d 13
40     d 14
41     d 16
42     (3, 27)
43
44     sage: f5 = F5R() # F5 with interreduced B
45     sage: gb = f5(J)
46     sage: f5.zero_reductions, len(gb)
47     d 7
48     d 9
49     d 11
50     d 6
51     d 7
52     d 8
53     d 9
54     d 10
55     d 11
56     verbose 0 (179: 283.py, top_reduction) reduction to zero.
57     verbose 0 (179: 283.py, top_reduction) reduction to zero.
58     verbose 0 (179: 283.py, top_reduction) reduction to zero.
59     d 12
60     d 13
61     d 14
62     d 16
63     (3, 18)
64
65     sage: f5 = F5C() # F5 with interreduced B and Gprev
66     sage: gb = f5(J)
67     sage: f5.zero_reductions, len(gb)
68     d 7
69     d 9
70     d 11
71     d 6
72     d 7
73     d 8
74     d 9
75     d 10
76     d 11
77     verbose 0 (179: 283.py, top_reduction) reduction to zero.
78     verbose 0 (179: 283.py, top_reduction) reduction to zero.
79     verbose 0 (179: 283.py, top_reduction) reduction to zero.
80     d 12
81     d 13
82     d 14
83     d 16
84     (3, 18)
85 """
86
87 divides = lambda x,y: x.parent().monomial_divides(x,y)
88 LCM = lambda f,g: f.parent().monomial_lcm(f,g)
89 LM = lambda f: f.lm()
90 LT = lambda f: f.lt()
91
92 def compare_by_degree(f,g):
93     if f.total_degree() > g.total_degree():
94         return 1
95     elif f.total_degree() < g.total_degree():
96         return -1
97     else:
98         return cmp(f.lm(),g.lm())
99
100 class F5:
101     def __init__(self, F=None):
102         if F is not None:
103             self.Rules = [[] for _ in range(len(F))]
104             self.L = [0]
105             self.zero_reductions = 0
106
107     def poly(self, i):
108         return self.L[i][1]
109
110     def sig(self, i):
111         return self.L[i][0]
112
113     def __call__(self, F):
114         if isinstance(F, sage.rings.polynomial.multi_polynomial_ideal.MPolynomialIdeal):
115             F = F.reduced_basis()
116         else:
117             F = Ideal(list(F)).reduced_basis()
118         if not all(f.is_homogeneous() for f in F):
119             F = Ideal(F).homogenize()
120             F = F.gens()
121         return self.basis(F)
122
123     def basis(self, F):
124         poly = self.poly
125         incremental_basis = self.incremental_basis
126
127         self.__init__(F)
128
129         Rules = self.Rules
130         L = self.L
131
132         m = len(F)
133         F = sorted(F, cmp=compare_by_degree)
134
135         f0 = F[0]
136         L[0] = (Signature(1, 0), f0*f0.lc()**(-1))
137         Rules[0] = list()
138
139         Gprev = set([0])
140         B = [f0]
141
142         for i in xrange(1,m):
143             Gcurr = incremental_basis(F[i], i, B, Gprev)
144             if any(poly(lambd) == 1 for lambd in Gcurr):
145                 return set(1)
146             Gprev = Gcurr
147             B = [poly(l) for l in Gprev]
148         return B
149
150     def incremental_basis(self, f, i, B, Gprev):
151         L = self.L
152         critical_pair = self.critical_pair
153         compute_spols = self.compute_spols
154         reduction = self.reduction
155         Rules = self.Rules
156
157         L.append( (Signature(1,i), f*f.lc()**(-1)) )
158         curr_idx = len(L) - 1
159         Gcurr = Gprev.union([curr_idx])
160         Rules[i] = list()
161
162         P = reduce(lambda x,y: x.union(y), [critical_pair(curr_idx, j, i, Gprev) for j in Gprev], set())
163         while len(P) != 0:
164             d = min(t.degree() for (t,k,u,l,v) in P)
165             print "d", d
166             Pd = [(t,k,u,l,v) for (t,k,u,l,v) in P if t.degree() == d]
167             P = P.difference(Pd)
168             S = compute_spols(Pd)
169             R = reduction(S, B, Gprev, Gcurr)
170             for k in R:
171                 P = reduce(lambda x,y: x.union(y), [critical_pair(j, k, i, Gprev) for j in Gcurr], P)
173         return Gcurr
174
175     def critical_pair(self, k, l, i, Gprev):
176         poly = self.poly
177         sig = self.sig
178         is_top_reducible = self.is_top_reducible
179         is_rewritable = self.is_rewritable
180
181         #print "crit_pair(%s,%s,%s,%s)"%(k, l, i, Gprev)
182         #print self.L
183         tk = poly(k).lt()
184         tl = poly(l).lt()
185         t = LCM(tk, tl)
186         u0 = t//tk
187         u1 = t//tl
188         m0, e0 = sig(k)
189         m1, e1 = sig(l)
190         if e0 == e1 and u0*m0 == u1*m1:
191             return set()
192         #print "test1", e0, i, u0, m0
193         #print "test2", e1, i, u1, m1
194         if e0 == i and is_top_reducible(u0*m0, Gprev):
195             #print "test1 done"
196             return set()
197         if e1 == i and is_top_reducible(u1*m1, Gprev):
198             #print "test2 done"
199             return set()
200         if is_rewritable(u0, k) or is_rewritable(u1, l):
201             #print "test3", is_rewritable(u0, k), is_rewritable(u1, l)
202             return set()
203         if u0 * sig(k) < u1 * sig(l):
204             u0, u1 = u1, u0
205             k, l = l, k
206         return set([(t,k,u0,l,u1)])
207
208     def compute_spols(self, P):
209         poly = self.poly
210         sig = self.sig
211         spol = self.spol
212         is_rewritable = self.is_rewritable
214
215         L = self.L
216
217         S = list()
218         P = sorted(P, key=lambda x: x[0])
219         for (t,k,u,l,v) in P:
220             if not is_rewritable(u,k) and not is_rewritable(v,l):
221                 s = spol(poly(k), poly(l))
222                 if s != 0:
223                     L.append( (u * sig(k), s) )
224                     add_rule(u * sig(k), len(L)-1)
225                     S.append(len(L)-1)
226         S = sorted(S, key=lambda x: sig(x))
227         return S
228
229     def spol(self, f, g):
230         return LCM(LM(f),LM(g)) // LT(f) * f - LCM(LM(f),LM(g)) // LT(g) * g
231
232     def reduction(self, S, B, Gprev, Gcurr):
233         L = self.L
234         sig = self.sig
235         poly = self.poly
236         top_reduction = self.top_reduction
237
238         to_do = S
239         completed = set()
240         while len(to_do):
241             k, to_do = to_do[0], to_do[1:]
242             h = poly(k).reduce(B)
243             L[k] = (sig(k), h)
244             newly_completed, redo = top_reduction(k, Gprev, Gcurr.union(completed))
245             completed = completed.union( newly_completed )
246             for j in redo:
247                 # insert j in to_do, sorted by increasing signature
248                 to_do.append(j)
249                 to_do.sort(key=lambda x: sig(x))
250         return completed
251
252     def top_reduction(self, k, Gprev, Gcurr):
253         find_reductor = self.find_reductor
255         poly = self.poly
256         sig = self.sig
257         L = self.L
258
259         if poly(k) == 0:
260             verbose("reduction to zero.",level=0)
261             self.zero_reductions += 1
262             return set(),set()
263         p = poly(k)
264         J = find_reductor(k, Gprev, Gcurr)
265         if J == set():
266             L[k] = ( sig(k), p * p.lc()**(-1) )
267             return set([k]),set()
268         j = J.pop()
269         q = poly(j)
270         u = p.lt()//q.lt()
271         p = p - u*q
272         if p != 0:
273             p = p * p.lc()**(-1)
274         if u * sig(j) < sig(k):
275             L[k] = (sig(k), p)
276             return set(), set([k])
277         else:
278             L.append((u * sig(j), p))
279             add_rule(u * sig(j), len(L)-1)
280             return set(), set([k, len(L)-1])
281
282     def find_reductor(self, k, Gprev, Gcurr):
283         is_rewritable = self.is_rewritable
284         is_top_reducible = self.is_top_reducible
285
286         poly = self.poly
287         sig = self.sig
288         t = poly(k).lt()
289         for j in Gcurr:
290             tprime = poly(j).lt()
291             if divides(tprime,t):
292                 u = t // tprime
293                 mj, ej = sig(j)
294                 if u * sig(j) != sig(k) and not is_rewritable(u, j) \
295                         and not is_top_reducible(u*mj, Gprev):
296                     return set([j])
297         return set()
298
299     def add_rule(self, s, k):
300         self.Rules[s[1]].append( (s[0],k) )
301
302     def is_rewritable(self, u, k):
303         j = self.find_rewriting(u, k)
304         return j != k
305
306     def find_rewriting(self, u, k):
307         Rules = self.Rules
308         mk, v = self.sig(k)
309         for ctr in reversed(xrange(len(Rules[v]))):
310             mj, j = Rules[v][ctr]
311             if divides(mj, u * mk):
312                 return j
313         return k
314
315     def is_top_reducible(self, t, l):
316         poly = self.poly
317         for g in l:
318             if divides(poly(g).lt(), t):
319                 return True
320         return False
321
322 class F5R(F5):
323     def basis(self, F):
324         poly = self.poly
325         incremental_basis = self.incremental_basis
326
327         self.__init__(F)
328
329         Rules = self.Rules
330         L = self.L
331
332         m = len(F)
333         F = sorted(F, cmp=compare_by_degree)
334
335         f0 = F[0]
336         L[0] = (Signature(1, 0), f0*f0.lc()**(-1))
337         Rules[0] = list()
338
339         Gprev = set([0])
340         B = [f0]
341
342         for i in xrange(1,m):
343             Gcurr = incremental_basis(F[i], i, B, Gprev)
344             if any(poly(lambd) == 1 for lambd in Gcurr):
345                 return set(1)
346             Gprev = Gcurr
347             B = Ideal([poly(l) for l in Gprev]).reduced_basis()
348
349         return B
350
351 class F5C(F5):
352     def basis(self, F):
353         poly = self.poly
354
355         self.__init__(F)
356
357         Rules = self.Rules
358         L = self.L
359
360         m = len(F)
361         F = sorted(F, cmp=compare_by_degree)
362
363         f0 = F[0]
364         L[0] = (Signature(1, 0), f0*f0.lc()**(-1))
365         Rules[0] = list()
366
367         Gprev = set([0])
368         B = set([f0])
369
370         for i in xrange(1,m):
371             Gcurr = self.incremental_basis(F[i], B, Gprev)
372             if any(poly(lambd) == 1 for lambd in Gcurr):
373                 return set(1)
374             B = Ideal([poly(l) for l in Gcurr]).reduced_basis()
375             if i != m-1:
376                 Gprev = self.setup_reduced_basis(B)
377         return B
378
379     def incremental_basis(self, f, B, Gprev):
380         L = self.L
381         critical_pair = self.critical_pair
382         compute_spols = self.compute_spols
383         reduction = self.reduction
384         Rules = self.Rules
385
386         i = len(L)
387         L.append( (Signature(1,i), f*f.lc()**(-1)) )
388         Rules.append( list() )
389         Gcurr = Gprev.union([i])
390         P = reduce(lambda x,y: x.union(y), [critical_pair(i, j, i, Gprev) for j in Gprev], set())
391         while len(P) != 0:
392             d = min(t.degree() for (t,k,u,l,v) in P)
393             print "d", d
394             Pd = [(t,k,u,l,v) for (t,k,u,l,v) in P if t.degree() == d]
395             P = P.difference(Pd)
396             S = compute_spols(Pd)
397             R = reduction(S, B, Gprev, Gcurr)
398             for k in R:
399                 P = reduce(lambda x,y: x.union(y), [critical_pair(j, k, i, Gprev) for j in Gcurr], P)
401         return Gcurr
402
403     def setup_reduced_basis(self, B):
405         self.L = range(len(B))
406         self.Rules = [[] for _ in range(len(B))]
407
408         L = self.L
409         Rules = self.Rules
410         Gcurr = set()
411
412         for i in range(len(B)):
413             L[i] = (Signature(1,i), B[i])
414             Gcurr.add( i )
415             Rules[i] = []
416             t = B[i].lt()
417             for j in range(i+1, len(B)):
418                 u = LCM(t, B[j].lt())//B[j].lt()
419                 add_rule( Signature(u, j), -1 )
420         return Gcurr
421
422 from UserList import UserList
423
424 class Signature(UserList):
425      def __init__(self, monomial, index):
426          UserList.__init__(self,[monomial, index])
427
428      def __lt__(self, other):
429          if self[1] < other[1]:
430              return True
431          elif self[1] > other[1]:
432              return False
433          else:
434              return (self[0] < other[0])
435
436      def __eq__(self, other):
437          return self[0] == other[0] and self[1] == other[1]
438
439      def __neq__(self, other):
440          return self[0] != other[0] or self[1] != other[1]
441
442      def __rmul__(self, other):
443          if isinstance(self, Signature):
444              return Signature(other * self[0], self[1])
445          else:
446              raise TypeError
447
448 f5 = F5C()
```

## 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

You are not allowed to attach a file to this page.