## "Fossies" - the Fresh Open Source Software Archive ### Member "NZMATH-1.2.0/nzmath/poly/groebner.py" (19 Nov 2012, 4802 Bytes) of package /linux/misc/old/NZMATH-1.2.0.tar.gz:

As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) Python source code syntax highlighting (style: standard) with prefixed line numbers. Alternatively you can here view or download the uninterpreted source code file. For more information about "groebner.py" see the Fossies "Dox" file reference documentation.

```    1 """
2 Groebner basis
3 """
4
5 from __future__ import division
6 import nzmath.compatibility
7 import nzmath.ring as ring
8
9
10 def s_polynomial(f, g, order):
11     """
12     Return S-polynomial of f and g with respect to the order.
13
14     S(f, g) = (lc(g)*T/lb(f))*f - (lc(f)*T/lb(g))*g,
15     where T = lcm(lb(f), lb(g)).
16     """
17     f_lb, f_lc = order.leading_term(f) # term = (base, coeff)
18     g_lb, g_lc = order.leading_term(g)
19     t = f_lb.lcm(g_lb)
20     return f.term_mul((t - f_lb, g_lc)) - g.term_mul((t - g_lb, f_lc))
21
22 def step_reduce(reducee, reducer, order):
23     """
24     Return the reduced polynomial of reducee by reducer.  That is, if
25     one of reducee's bases is divisible by the leading base of reducer
26     with respect to the order, the returned polynomial is the result
27     of canceling out the term.
28     """
29     lb, lc = order.leading_term(reducer) # term = (base, coeff)
30     for b, c in reducee.iterterms():
31         if b.lcm(lb) == b:
32             return reducee + reducer.term_mul((b - lb, -c / lc)), True
33     return reducee, False
34
35 def reduce_closure(f, reducers, order):
36     """
37     Return normalized form of f with respect to reducer, a set of
38     polynomials, and order.
39     """
40     while True:
41         for reducer in reducers:
42             f, reduced = step_reduce(f, reducer, order)
43             if reduced:
44                 break
45         else:
46             return f
47
48 def buchberger(generating, order):
49     """
50     Return a Groebner basis of the ideal generated by given generating
51     set of polynomials with respect to the order.
52
53     Be careful, this implementation is very naive.
54     """
55     pairs = []
56     for f in generating:
57         for g in generating:
58             if f is g:
59                 continue
60             pairs.append((f, g))
61     groebner = list(generating)
62
63     # main loop
64     while pairs:
65         f, g = pairs.pop()
66         h = reduce_closure(s_polynomial(f, g, order), groebner, order)
67         if h:
68             for f in groebner:
69                 pairs.append((f, h))
70             groebner.append(h)
71
72     # finish
73     return groebner
74
75 def reduce_groebner(gbasis, order):
76     """
77     Return the reduced Groebner basis constructed from a Groebner
78     basis.
79
80     1) lb(f) divides lb(g) => g is not in reduced Groebner basis
81     2) monic
82     """
83     reduced_basis = []
84     lb_rel = dict([(order.leading_term(g), g) for g in gbasis])
85     lbs = sorted([order.leading_term(g) for g in gbasis])[::-1]
86     for i, lbi in enumerate(lbs):
87         for lbj in lbs[(len(lbs) - 1):i:-1]:
88             if lbi == lbj.lcm(lbi):
89                 # divisor found
90                 break
91         else:
92             g = lb_rel[lbi]
93             if g[lbi] != ring.getRing(g[lbi]).one:
94                 # make it monic
95                 g = g.scalar_mul(ring.inverse(g[lbi]))
96             reduced_basis.append(g)
97     return reduced_basis
98
99 def normal_strategy(generating, order):
100     """
101     Return a Groebner basis of the ideal generated by given generating
102     set of polynomials with respect to the order.
103
104     This function uses the 'normal strategy'.
105     """
106     groebner = list(generating) # omit reduction
107     pairs, lcms = [], []
108     treat = set()
109     for i, f in enumerate(groebner):
110         lb_f = order.leading_term(f)
111         for j, g in enumerate(groebner):
112             if j <= i:
113                 continue
114             lb_g = order.leading_term(g)
115             lcm_f_g = lb_f.lcm(lb_g)
116             if lcm_f_g == lb_f + lb_g: # disjoint
119             else:
120                 # keep lcms sorted, and so pairs in paralell
121                 k = order.bisect(lcms, lcm_f_g)
122                 pairs.insert(k, (i, j))
123                 lcms.insert(k, lcm_f_g)
124
125     # main loop
126     while pairs:
127         i, j = pairs.pop(0)
128         lcm_f_g = lcms.pop(0)
130         for p in range(len(groebner)):
131             if (i, p) in treat and (p, j) in treat:
132                 pivot = order.leading_term(groebner[p])
133                 if pivot.lcm(lcm_f_g) == lcm_f_g:
134                     break
135         else:
136             f, g = groebner[i], groebner[j]
137             h = reduce_closure(s_polynomial(f, g, order), groebner, order)
138             if h:
139                 lb_h = order.leading_term(h)
140                 hindex = len(groebner)
141                 for i, f in enumerate(groebner):
142                     lb_f = order.leading_term(f)
143                     lcm_f_h = lb_f.lcm(lb_h)
144                     if lcm_f_h == lb_f + lb_h: # disjoint