NZMATH  1.2.0
About: NZMATH is a Python based number theory oriented calculation system.
  Fossies Dox: NZMATH-1.2.0.tar.gz  ("inofficial" and yet experimental doxygen-generated source code documentation)  

groebner.py
Go to the documentation of this file.
1 """
2 Groebner basis
3 """
4 
5 from __future__ import division
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)[0], g) for g in gbasis])
85  lbs = sorted([order.leading_term(g)[0] 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)[0]
111  for j, g in enumerate(groebner):
112  if j <= i:
113  continue
114  lb_g = order.leading_term(g)[0]
115  lcm_f_g = lb_f.lcm(lb_g)
116  if lcm_f_g == lb_f + lb_g: # disjoint
117  treat.add((i, j))
118  treat.add((j, i))
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)
129  treat.add((i, j))
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])[0]
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)[0]
140  hindex = len(groebner)
141  for i, f in enumerate(groebner):
142  lb_f = order.leading_term(f)[0]
143  lcm_f_h = lb_f.lcm(lb_h)
144  if lcm_f_h == lb_f + lb_h: # disjoint
145  treat.add((i, hindex))
146  else:
147  k = order.bisect(lcms, lcm_f_h)
148  pairs.insert(k, (i, hindex))
149  lcms.insert(k, lcm_f_h)
150  groebner.append(h)
151 
152  # finish
153  return groebner
nzmath.poly.groebner.step_reduce
def step_reduce(reducee, reducer, order)
Definition: groebner.py:22
nzmath.bigrange.range
def range(start, stop=None, step=None)
Definition: bigrange.py:19
nzmath.ring
Definition: ring.py:1
nzmath.poly.groebner.normal_strategy
def normal_strategy(generating, order)
Definition: groebner.py:99
nzmath.poly.groebner.s_polynomial
def s_polynomial(f, g, order)
Definition: groebner.py:10
nzmath.poly.groebner.reduce_closure
def reduce_closure(f, reducers, order)
Definition: groebner.py:35
nzmath.poly.groebner.reduce_groebner
def reduce_groebner(gbasis, order)
Definition: groebner.py:75
nzmath.compatibility
Definition: compatibility.py:1
nzmath.poly.groebner.buchberger
def buchberger(generating, order)
Definition: groebner.py:48