Source code for sympy.polys.groebnertools

"""Groebner bases algorithms. """

from __future__ import print_function, division

from sympy.polys.monomials import monomial_mul, monomial_div, monomial_lcm, monomial_divides, term_div
from sympy.polys.orderings import lex
from sympy.polys.polyerrors import DomainError
from sympy.polys.polyconfig import query
from sympy.core.symbol import Dummy
from sympy.core.compatibility import xrange

[docs]def groebner(seq, ring, method=None): """ Computes Groebner basis for a set of polynomials in K[X]. Wrapper around the (default) improved Buchberger and the other algorithms for computing Groebner bases. The choice of algorithm can be changed via method argument or :func:setup from :mod:sympy.polys.polyconfig, where method can be either buchberger or f5b. """ if method is None: method = query('groebner') _groebner_methods = { 'buchberger': _buchberger, 'f5b': _f5b, } try: _groebner = _groebner_methods[method] except KeyError: raise ValueError("'%s' is not a valid Groebner bases algorithm (valid are 'buchberger' and 'f5b')" % method) domain, orig = ring.domain, None if not domain.has_Field or not domain.has_assoc_Field: try: orig, ring = ring, ring.clone(domain=domain.get_field()) except DomainError: raise DomainError("can't compute a Groebner basis over %s" % domain) else: seq = [ s.set_ring(ring) for s in seq ] G = _groebner(seq, ring) if orig is not None: G = [ g.clear_denoms()[1].set_ring(orig) for g in G ] return G
def _buchberger(f, ring): """ Computes Groebner basis for a set of polynomials in K[X]. Given a set of multivariate polynomials F, finds another set G, such that Ideal F = Ideal G and G is a reduced Groebner basis. The resulting basis is unique and has monic generators if the ground domains is a field. Otherwise the result is non-unique but Groebner bases over e.g. integers can be computed (if the input polynomials are monic). Groebner bases can be used to choose specific generators for a polynomial ideal. Because these bases are unique you can check for ideal equality by comparing the Groebner bases. To see if one polynomial lies in an ideal, divide by the elements in the base and see if the remainder vanishes. They can also be used to solve systems of polynomial equations as, by choosing lexicographic ordering, you can eliminate one variable at a time, provided that the ideal is zero-dimensional (finite number of solutions). References ========== 1. [Bose03]_ 2. [Giovini91]_ 3. [Ajwa95]_ 4. [Cox97]_ Algorithm used: an improved version of Buchberger's algorithm as presented in T. Becker, V. Weispfenning, Groebner Bases: A Computational Approach to Commutative Algebra, Springer, 1993, page 232. """ order = ring.order domain = ring.domain monomial_mul = ring.monomial_mul monomial_div = ring.monomial_div monomial_lcm = ring.monomial_lcm def select(P): # normal selection strategy # select the pair with minimum LCM(LM(f), LM(g)) pr = min(P, key=lambda pair: order(monomial_lcm(f[pair[0]].LM, f[pair[1]].LM))) return pr def normal(g, J): h = g.rem([ f[j] for j in J ]) if not h: return None else: h = h.monic() if not h in I: I[h] = len(f) f.append(h) return h.LM, I[h] def update(G, B, ih): # update G using the set of critical pairs B and h # [BW] page 230 h = f[ih] mh = h.LM # filter new pairs (h, g), g in G C = G.copy() D = set() while C: # select a pair (h, g) by popping an element from C ig = C.pop() g = f[ig] mg = g.LM LCMhg = monomial_lcm(mh, mg) def lcm_divides(ip): # LCM(LM(h), LM(p)) divides LCM(LM(h), LM(g)) m = monomial_lcm(mh, f[ip].LM) return monomial_div(LCMhg, m) # HT(h) and HT(g) disjoint: mh*mg == LCMhg if monomial_mul(mh, mg) == LCMhg or ( not any(lcm_divides(ipx) for ipx in C) and not any(lcm_divides(pr[1]) for pr in D)): D.add((ih, ig)) E = set() while D: # select h, g from D (h the same as above) ih, ig = D.pop() mg = f[ig].LM LCMhg = monomial_lcm(mh, mg) if not monomial_mul(mh, mg) == LCMhg: E.add((ih, ig)) # filter old pairs B_new = set() while B: # select g1, g2 from B (-> CP) ig1, ig2 = B.pop() mg1 = f[ig1].LM mg2 = f[ig2].LM LCM12 = monomial_lcm(mg1, mg2) # if HT(h) does not divide lcm(HT(g1), HT(g2)) if not monomial_div(LCM12, mh) or \ monomial_lcm(mg1, mh) == LCM12 or \ monomial_lcm(mg2, mh) == LCM12: B_new.add((ig1, ig2)) B_new |= E # filter polynomials G_new = set() while G: ig = G.pop() mg = f[ig].LM if not monomial_div(mg, mh): G_new.add(ig) G_new.add(ih) return G_new, B_new # end of update ################################ if not f: return [] # replace f with a reduced list of initial polynomials; see [BW] page 203 f1 = f[:] while True: f = f1[:] f1 = [] for i in range(len(f)): p = f[i] r = p.rem(f[:i]) if r: f1.append(r.monic()) if f == f1: break I = {} # ip = I[p]; p = f[ip] F = set() # set of indices of polynomials G = set() # set of indices of intermediate would-be Groebner basis CP = set() # set of pairs of indices of critical pairs for i, h in enumerate(f): I[h] = i F.add(i) ##################################### # algorithm GROEBNERNEWS2 in [BW] page 232 while F: # select p with minimum monomial according to the monomial ordering h = min([f[x] for x in F], key=lambda f: order(f.LM)) ih = I[h] F.remove(ih) G, CP = update(G, CP, ih) # count the number of critical pairs which reduce to zero reductions_to_zero = 0 while CP: ig1, ig2 = select(CP) CP.remove((ig1, ig2)) h = spoly(f[ig1], f[ig2], ring) # ordering divisors is on average more efficient [Cox] page 111 G1 = sorted(G, key=lambda g: order(f[g].LM)) ht = normal(h, G1) if ht: G, CP = update(G, CP, ht[1]) else: reductions_to_zero += 1 ###################################### # now G is a Groebner basis; reduce it Gr = set() for ig in G: ht = normal(f[ig], G - set([ig])) if ht: Gr.add(ht[1]) Gr = [f[ig] for ig in Gr] # order according to the monomial ordering Gr = sorted(Gr, key=lambda f: order(f.LM), reverse=True) return Gr
[docs]def spoly(p1, p2, ring): """ Compute LCM(LM(p1), LM(p2))/LM(p1)*p1 - LCM(LM(p1), LM(p2))/LM(p2)*p2 This is the S-poly provided p1 and p2 are monic """ LM1 = p1.LM LM2 = p2.LM LCM12 = ring.monomial_lcm(LM1, LM2) m1 = ring.monomial_div(LCM12, LM1) m2 = ring.monomial_div(LCM12, LM2) s1 = p1.mul_monom(m1) s2 = p2.mul_monom(m2) s = s1 - s2 return s # F5B # convenience functions
[docs]def red_groebner(G, ring): """ Compute reduced Groebner basis, from BeckerWeispfenning93, p. 216 Selects a subset of generators, that already generate the ideal and computes a reduced Groebner basis for them. """ def reduction(P): """ The actual reduction algorithm. """ Q = [] for i, p in enumerate(P): h = p.rem(P[:i] + P[i + 1:]) if h: Q.append(h) return [p.monic() for p in Q] F = G H = [] while F: f0 = F.pop() if not any(monomial_divides(f.LM, f0.LM) for f in F + H): H.append(f0) # Becker, Weispfenning, p. 217: H is Groebner basis of the ideal generated by G. return reduction(H)
[docs]def is_groebner(G, ring): """ Check if G is a Groebner basis. """ for i in xrange(len(G)): for j in xrange(i + 1, len(G)): s = spoly(G[i], G[j]) s = s.rem(G) if s: return False return True
[docs]def is_minimal(G, ring): """ Checks if G is a minimal Groebner basis. """ order = ring.order domain = ring.domain G.sort(key=lambda g: order(g.LM)) for i, g in enumerate(G): if g.LC != domain.one: return False for h in G[:i] + G[i + 1:]: if monomial_divides(h.LM, g.LM): return False return True
[docs]def is_reduced(G, ring): """ Checks if G is a reduced Groebner basis. """ order = ring.order domain = ring.domain G.sort(key=lambda g: order(g.LM)) for i, g in enumerate(G): if g.LC != domain.one: return False for term in g: for h in G[:i] + G[i + 1:]: if monomial_divides(h.LM, term[0]): return False return True
def groebner_lcm(f, g): """ Computes LCM of two polynomials using Groebner bases. The LCM is computed as the unique generater of the intersection of the two ideals generated by f and g. The approach is to compute a Groebner basis with respect to lexicographic ordering of t*f and (1 - t)*g, where t is an unrelated variable and then filtering out the solution that doesn't contain t. References ========== 1. [Cox97]_ """ if f.ring != g.ring: raise ValueError("Values should be equal") ring = f.ring domain = ring.domain if not f or not g: return ring.zero if len(f) <= 1 and len(g) <= 1: monom = monomial_lcm(f.LM, g.LM) coeff = domain.lcm(f.LC, g.LC) return ring.term_new(monom, coeff) fc, f = f.primitive() gc, g = g.primitive() lcm = domain.lcm(fc, gc) f_terms = [ ((1,) + monom, coeff) for monom, coeff in f.terms() ] g_terms = [ ((0,) + monom, coeff) for monom, coeff in g.terms() ] \ + [ ((1,) + monom,-coeff) for monom, coeff in g.terms() ] t = Dummy("t") t_ring = ring.clone(symbols=(t,) + ring.symbols, order=lex) F = t_ring.from_terms(f_terms) G = t_ring.from_terms(g_terms) basis = groebner([F, G], t_ring) def is_independent(h, j): return all(not monom[j] for monom in h.monoms()) H = [ h for h in basis if is_independent(h, 0) ] h_terms = [ (monom[1:], coeff*lcm) for monom, coeff in H[0].terms() ] h = ring.from_terms(h_terms) return h def groebner_gcd(f, g): """Computes GCD of two polynomials using Groebner bases. """ if f.ring != g.ring: raise ValueError("Values should be equal") domain = f.ring.domain if not domain.has_Field: fc, f = f.primitive() gc, g = g.primitive() gcd = domain.gcd(fc, gc) H = (f*g).quo([groebner_lcm(f, g)]) if len(H) != 1: raise ValueError("Length should be 1") h = H[0] if not domain.has_Field: return gcd*h else: return h.monic()