Coverage for /usr/lib/python3/dist-packages/sympy/polys/groebnertools.py: 9%
374 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1"""Groebner bases algorithms. """
4from sympy.core.symbol import Dummy
5from sympy.polys.monomials import monomial_mul, monomial_lcm, monomial_divides, term_div
6from sympy.polys.orderings import lex
7from sympy.polys.polyerrors import DomainError
8from sympy.polys.polyconfig import query
10def groebner(seq, ring, method=None):
11 """
12 Computes Groebner basis for a set of polynomials in `K[X]`.
14 Wrapper around the (default) improved Buchberger and the other algorithms
15 for computing Groebner bases. The choice of algorithm can be changed via
16 ``method`` argument or :func:`sympy.polys.polyconfig.setup`, where
17 ``method`` can be either ``buchberger`` or ``f5b``.
19 """
20 if method is None:
21 method = query('groebner')
23 _groebner_methods = {
24 'buchberger': _buchberger,
25 'f5b': _f5b,
26 }
28 try:
29 _groebner = _groebner_methods[method]
30 except KeyError:
31 raise ValueError("'%s' is not a valid Groebner bases algorithm (valid are 'buchberger' and 'f5b')" % method)
33 domain, orig = ring.domain, None
35 if not domain.is_Field or not domain.has_assoc_Field:
36 try:
37 orig, ring = ring, ring.clone(domain=domain.get_field())
38 except DomainError:
39 raise DomainError("Cannot compute a Groebner basis over %s" % domain)
40 else:
41 seq = [ s.set_ring(ring) for s in seq ]
43 G = _groebner(seq, ring)
45 if orig is not None:
46 G = [ g.clear_denoms()[1].set_ring(orig) for g in G ]
48 return G
50def _buchberger(f, ring):
51 """
52 Computes Groebner basis for a set of polynomials in `K[X]`.
54 Given a set of multivariate polynomials `F`, finds another
55 set `G`, such that Ideal `F = Ideal G` and `G` is a reduced
56 Groebner basis.
58 The resulting basis is unique and has monic generators if the
59 ground domains is a field. Otherwise the result is non-unique
60 but Groebner bases over e.g. integers can be computed (if the
61 input polynomials are monic).
63 Groebner bases can be used to choose specific generators for a
64 polynomial ideal. Because these bases are unique you can check
65 for ideal equality by comparing the Groebner bases. To see if
66 one polynomial lies in an ideal, divide by the elements in the
67 base and see if the remainder vanishes.
69 They can also be used to solve systems of polynomial equations
70 as, by choosing lexicographic ordering, you can eliminate one
71 variable at a time, provided that the ideal is zero-dimensional
72 (finite number of solutions).
74 Notes
75 =====
77 Algorithm used: an improved version of Buchberger's algorithm
78 as presented in T. Becker, V. Weispfenning, Groebner Bases: A
79 Computational Approach to Commutative Algebra, Springer, 1993,
80 page 232.
82 References
83 ==========
85 .. [1] [Bose03]_
86 .. [2] [Giovini91]_
87 .. [3] [Ajwa95]_
88 .. [4] [Cox97]_
90 """
91 order = ring.order
93 monomial_mul = ring.monomial_mul
94 monomial_div = ring.monomial_div
95 monomial_lcm = ring.monomial_lcm
97 def select(P):
98 # normal selection strategy
99 # select the pair with minimum LCM(LM(f), LM(g))
100 pr = min(P, key=lambda pair: order(monomial_lcm(f[pair[0]].LM, f[pair[1]].LM)))
101 return pr
103 def normal(g, J):
104 h = g.rem([ f[j] for j in J ])
106 if not h:
107 return None
108 else:
109 h = h.monic()
111 if h not in I:
112 I[h] = len(f)
113 f.append(h)
115 return h.LM, I[h]
117 def update(G, B, ih):
118 # update G using the set of critical pairs B and h
119 # [BW] page 230
120 h = f[ih]
121 mh = h.LM
123 # filter new pairs (h, g), g in G
124 C = G.copy()
125 D = set()
127 while C:
128 # select a pair (h, g) by popping an element from C
129 ig = C.pop()
130 g = f[ig]
131 mg = g.LM
132 LCMhg = monomial_lcm(mh, mg)
134 def lcm_divides(ip):
135 # LCM(LM(h), LM(p)) divides LCM(LM(h), LM(g))
136 m = monomial_lcm(mh, f[ip].LM)
137 return monomial_div(LCMhg, m)
139 # HT(h) and HT(g) disjoint: mh*mg == LCMhg
140 if monomial_mul(mh, mg) == LCMhg or (
141 not any(lcm_divides(ipx) for ipx in C) and
142 not any(lcm_divides(pr[1]) for pr in D)):
143 D.add((ih, ig))
145 E = set()
147 while D:
148 # select h, g from D (h the same as above)
149 ih, ig = D.pop()
150 mg = f[ig].LM
151 LCMhg = monomial_lcm(mh, mg)
153 if not monomial_mul(mh, mg) == LCMhg:
154 E.add((ih, ig))
156 # filter old pairs
157 B_new = set()
159 while B:
160 # select g1, g2 from B (-> CP)
161 ig1, ig2 = B.pop()
162 mg1 = f[ig1].LM
163 mg2 = f[ig2].LM
164 LCM12 = monomial_lcm(mg1, mg2)
166 # if HT(h) does not divide lcm(HT(g1), HT(g2))
167 if not monomial_div(LCM12, mh) or \
168 monomial_lcm(mg1, mh) == LCM12 or \
169 monomial_lcm(mg2, mh) == LCM12:
170 B_new.add((ig1, ig2))
172 B_new |= E
174 # filter polynomials
175 G_new = set()
177 while G:
178 ig = G.pop()
179 mg = f[ig].LM
181 if not monomial_div(mg, mh):
182 G_new.add(ig)
184 G_new.add(ih)
186 return G_new, B_new
187 # end of update ################################
189 if not f:
190 return []
192 # replace f with a reduced list of initial polynomials; see [BW] page 203
193 f1 = f[:]
195 while True:
196 f = f1[:]
197 f1 = []
199 for i in range(len(f)):
200 p = f[i]
201 r = p.rem(f[:i])
203 if r:
204 f1.append(r.monic())
206 if f == f1:
207 break
209 I = {} # ip = I[p]; p = f[ip]
210 F = set() # set of indices of polynomials
211 G = set() # set of indices of intermediate would-be Groebner basis
212 CP = set() # set of pairs of indices of critical pairs
214 for i, h in enumerate(f):
215 I[h] = i
216 F.add(i)
218 #####################################
219 # algorithm GROEBNERNEWS2 in [BW] page 232
221 while F:
222 # select p with minimum monomial according to the monomial ordering
223 h = min([f[x] for x in F], key=lambda f: order(f.LM))
224 ih = I[h]
225 F.remove(ih)
226 G, CP = update(G, CP, ih)
228 # count the number of critical pairs which reduce to zero
229 reductions_to_zero = 0
231 while CP:
232 ig1, ig2 = select(CP)
233 CP.remove((ig1, ig2))
235 h = spoly(f[ig1], f[ig2], ring)
236 # ordering divisors is on average more efficient [Cox] page 111
237 G1 = sorted(G, key=lambda g: order(f[g].LM))
238 ht = normal(h, G1)
240 if ht:
241 G, CP = update(G, CP, ht[1])
242 else:
243 reductions_to_zero += 1
245 ######################################
246 # now G is a Groebner basis; reduce it
247 Gr = set()
249 for ig in G:
250 ht = normal(f[ig], G - {ig})
252 if ht:
253 Gr.add(ht[1])
255 Gr = [f[ig] for ig in Gr]
257 # order according to the monomial ordering
258 Gr = sorted(Gr, key=lambda f: order(f.LM), reverse=True)
260 return Gr
262def spoly(p1, p2, ring):
263 """
264 Compute LCM(LM(p1), LM(p2))/LM(p1)*p1 - LCM(LM(p1), LM(p2))/LM(p2)*p2
265 This is the S-poly provided p1 and p2 are monic
266 """
267 LM1 = p1.LM
268 LM2 = p2.LM
269 LCM12 = ring.monomial_lcm(LM1, LM2)
270 m1 = ring.monomial_div(LCM12, LM1)
271 m2 = ring.monomial_div(LCM12, LM2)
272 s1 = p1.mul_monom(m1)
273 s2 = p2.mul_monom(m2)
274 s = s1 - s2
275 return s
277# F5B
279# convenience functions
282def Sign(f):
283 return f[0]
286def Polyn(f):
287 return f[1]
290def Num(f):
291 return f[2]
294def sig(monomial, index):
295 return (monomial, index)
298def lbp(signature, polynomial, number):
299 return (signature, polynomial, number)
301# signature functions
304def sig_cmp(u, v, order):
305 """
306 Compare two signatures by extending the term order to K[X]^n.
308 u < v iff
309 - the index of v is greater than the index of u
310 or
311 - the index of v is equal to the index of u and u[0] < v[0] w.r.t. order
313 u > v otherwise
314 """
315 if u[1] > v[1]:
316 return -1
317 if u[1] == v[1]:
318 #if u[0] == v[0]:
319 # return 0
320 if order(u[0]) < order(v[0]):
321 return -1
322 return 1
325def sig_key(s, order):
326 """
327 Key for comparing two signatures.
329 s = (m, k), t = (n, l)
331 s < t iff [k > l] or [k == l and m < n]
332 s > t otherwise
333 """
334 return (-s[1], order(s[0]))
337def sig_mult(s, m):
338 """
339 Multiply a signature by a monomial.
341 The product of a signature (m, i) and a monomial n is defined as
342 (m * t, i).
343 """
344 return sig(monomial_mul(s[0], m), s[1])
346# labeled polynomial functions
349def lbp_sub(f, g):
350 """
351 Subtract labeled polynomial g from f.
353 The signature and number of the difference of f and g are signature
354 and number of the maximum of f and g, w.r.t. lbp_cmp.
355 """
356 if sig_cmp(Sign(f), Sign(g), Polyn(f).ring.order) < 0:
357 max_poly = g
358 else:
359 max_poly = f
361 ret = Polyn(f) - Polyn(g)
363 return lbp(Sign(max_poly), ret, Num(max_poly))
366def lbp_mul_term(f, cx):
367 """
368 Multiply a labeled polynomial with a term.
370 The product of a labeled polynomial (s, p, k) by a monomial is
371 defined as (m * s, m * p, k).
372 """
373 return lbp(sig_mult(Sign(f), cx[0]), Polyn(f).mul_term(cx), Num(f))
376def lbp_cmp(f, g):
377 """
378 Compare two labeled polynomials.
380 f < g iff
381 - Sign(f) < Sign(g)
382 or
383 - Sign(f) == Sign(g) and Num(f) > Num(g)
385 f > g otherwise
386 """
387 if sig_cmp(Sign(f), Sign(g), Polyn(f).ring.order) == -1:
388 return -1
389 if Sign(f) == Sign(g):
390 if Num(f) > Num(g):
391 return -1
392 #if Num(f) == Num(g):
393 # return 0
394 return 1
397def lbp_key(f):
398 """
399 Key for comparing two labeled polynomials.
400 """
401 return (sig_key(Sign(f), Polyn(f).ring.order), -Num(f))
403# algorithm and helper functions
406def critical_pair(f, g, ring):
407 """
408 Compute the critical pair corresponding to two labeled polynomials.
410 A critical pair is a tuple (um, f, vm, g), where um and vm are
411 terms such that um * f - vm * g is the S-polynomial of f and g (so,
412 wlog assume um * f > vm * g).
413 For performance sake, a critical pair is represented as a tuple
414 (Sign(um * f), um, f, Sign(vm * g), vm, g), since um * f creates
415 a new, relatively expensive object in memory, whereas Sign(um *
416 f) and um are lightweight and f (in the tuple) is a reference to
417 an already existing object in memory.
418 """
419 domain = ring.domain
421 ltf = Polyn(f).LT
422 ltg = Polyn(g).LT
423 lt = (monomial_lcm(ltf[0], ltg[0]), domain.one)
425 um = term_div(lt, ltf, domain)
426 vm = term_div(lt, ltg, domain)
428 # The full information is not needed (now), so only the product
429 # with the leading term is considered:
430 fr = lbp_mul_term(lbp(Sign(f), Polyn(f).leading_term(), Num(f)), um)
431 gr = lbp_mul_term(lbp(Sign(g), Polyn(g).leading_term(), Num(g)), vm)
433 # return in proper order, such that the S-polynomial is just
434 # u_first * f_first - u_second * f_second:
435 if lbp_cmp(fr, gr) == -1:
436 return (Sign(gr), vm, g, Sign(fr), um, f)
437 else:
438 return (Sign(fr), um, f, Sign(gr), vm, g)
441def cp_cmp(c, d):
442 """
443 Compare two critical pairs c and d.
445 c < d iff
446 - lbp(c[0], _, Num(c[2]) < lbp(d[0], _, Num(d[2])) (this
447 corresponds to um_c * f_c and um_d * f_d)
448 or
449 - lbp(c[0], _, Num(c[2]) >< lbp(d[0], _, Num(d[2])) and
450 lbp(c[3], _, Num(c[5])) < lbp(d[3], _, Num(d[5])) (this
451 corresponds to vm_c * g_c and vm_d * g_d)
453 c > d otherwise
454 """
455 zero = Polyn(c[2]).ring.zero
457 c0 = lbp(c[0], zero, Num(c[2]))
458 d0 = lbp(d[0], zero, Num(d[2]))
460 r = lbp_cmp(c0, d0)
462 if r == -1:
463 return -1
464 if r == 0:
465 c1 = lbp(c[3], zero, Num(c[5]))
466 d1 = lbp(d[3], zero, Num(d[5]))
468 r = lbp_cmp(c1, d1)
470 if r == -1:
471 return -1
472 #if r == 0:
473 # return 0
474 return 1
477def cp_key(c, ring):
478 """
479 Key for comparing critical pairs.
480 """
481 return (lbp_key(lbp(c[0], ring.zero, Num(c[2]))), lbp_key(lbp(c[3], ring.zero, Num(c[5]))))
484def s_poly(cp):
485 """
486 Compute the S-polynomial of a critical pair.
488 The S-polynomial of a critical pair cp is cp[1] * cp[2] - cp[4] * cp[5].
489 """
490 return lbp_sub(lbp_mul_term(cp[2], cp[1]), lbp_mul_term(cp[5], cp[4]))
493def is_rewritable_or_comparable(sign, num, B):
494 """
495 Check if a labeled polynomial is redundant by checking if its
496 signature and number imply rewritability or comparability.
498 (sign, num) is comparable if there exists a labeled polynomial
499 h in B, such that sign[1] (the index) is less than Sign(h)[1]
500 and sign[0] is divisible by the leading monomial of h.
502 (sign, num) is rewritable if there exists a labeled polynomial
503 h in B, such thatsign[1] is equal to Sign(h)[1], num < Num(h)
504 and sign[0] is divisible by Sign(h)[0].
505 """
506 for h in B:
507 # comparable
508 if sign[1] < Sign(h)[1]:
509 if monomial_divides(Polyn(h).LM, sign[0]):
510 return True
512 # rewritable
513 if sign[1] == Sign(h)[1]:
514 if num < Num(h):
515 if monomial_divides(Sign(h)[0], sign[0]):
516 return True
517 return False
520def f5_reduce(f, B):
521 """
522 F5-reduce a labeled polynomial f by B.
524 Continuously searches for non-zero labeled polynomial h in B, such
525 that the leading term lt_h of h divides the leading term lt_f of
526 f and Sign(lt_h * h) < Sign(f). If such a labeled polynomial h is
527 found, f gets replaced by f - lt_f / lt_h * h. If no such h can be
528 found or f is 0, f is no further F5-reducible and f gets returned.
530 A polynomial that is reducible in the usual sense need not be
531 F5-reducible, e.g.:
533 >>> from sympy.polys.groebnertools import lbp, sig, f5_reduce, Polyn
534 >>> from sympy.polys import ring, QQ, lex
536 >>> R, x,y,z = ring("x,y,z", QQ, lex)
538 >>> f = lbp(sig((1, 1, 1), 4), x, 3)
539 >>> g = lbp(sig((0, 0, 0), 2), x, 2)
541 >>> Polyn(f).rem([Polyn(g)])
542 0
543 >>> f5_reduce(f, [g])
544 (((1, 1, 1), 4), x, 3)
546 """
547 order = Polyn(f).ring.order
548 domain = Polyn(f).ring.domain
550 if not Polyn(f):
551 return f
553 while True:
554 g = f
556 for h in B:
557 if Polyn(h):
558 if monomial_divides(Polyn(h).LM, Polyn(f).LM):
559 t = term_div(Polyn(f).LT, Polyn(h).LT, domain)
560 if sig_cmp(sig_mult(Sign(h), t[0]), Sign(f), order) < 0:
561 # The following check need not be done and is in general slower than without.
562 #if not is_rewritable_or_comparable(Sign(gp), Num(gp), B):
563 hp = lbp_mul_term(h, t)
564 f = lbp_sub(f, hp)
565 break
567 if g == f or not Polyn(f):
568 return f
571def _f5b(F, ring):
572 """
573 Computes a reduced Groebner basis for the ideal generated by F.
575 f5b is an implementation of the F5B algorithm by Yao Sun and
576 Dingkang Wang. Similarly to Buchberger's algorithm, the algorithm
577 proceeds by computing critical pairs, computing the S-polynomial,
578 reducing it and adjoining the reduced S-polynomial if it is not 0.
580 Unlike Buchberger's algorithm, each polynomial contains additional
581 information, namely a signature and a number. The signature
582 specifies the path of computation (i.e. from which polynomial in
583 the original basis was it derived and how), the number says when
584 the polynomial was added to the basis. With this information it
585 is (often) possible to decide if an S-polynomial will reduce to
586 0 and can be discarded.
588 Optimizations include: Reducing the generators before computing
589 a Groebner basis, removing redundant critical pairs when a new
590 polynomial enters the basis and sorting the critical pairs and
591 the current basis.
593 Once a Groebner basis has been found, it gets reduced.
595 References
596 ==========
598 .. [1] Yao Sun, Dingkang Wang: "A New Proof for the Correctness of F5
599 (F5-Like) Algorithm", https://arxiv.org/abs/1004.0084 (specifically
600 v4)
602 .. [2] Thomas Becker, Volker Weispfenning, Groebner bases: A computational
603 approach to commutative algebra, 1993, p. 203, 216
604 """
605 order = ring.order
607 # reduce polynomials (like in Mario Pernici's implementation) (Becker, Weispfenning, p. 203)
608 B = F
609 while True:
610 F = B
611 B = []
613 for i in range(len(F)):
614 p = F[i]
615 r = p.rem(F[:i])
617 if r:
618 B.append(r)
620 if F == B:
621 break
623 # basis
624 B = [lbp(sig(ring.zero_monom, i + 1), F[i], i + 1) for i in range(len(F))]
625 B.sort(key=lambda f: order(Polyn(f).LM), reverse=True)
627 # critical pairs
628 CP = [critical_pair(B[i], B[j], ring) for i in range(len(B)) for j in range(i + 1, len(B))]
629 CP.sort(key=lambda cp: cp_key(cp, ring), reverse=True)
631 k = len(B)
633 reductions_to_zero = 0
635 while len(CP):
636 cp = CP.pop()
638 # discard redundant critical pairs:
639 if is_rewritable_or_comparable(cp[0], Num(cp[2]), B):
640 continue
641 if is_rewritable_or_comparable(cp[3], Num(cp[5]), B):
642 continue
644 s = s_poly(cp)
646 p = f5_reduce(s, B)
648 p = lbp(Sign(p), Polyn(p).monic(), k + 1)
650 if Polyn(p):
651 # remove old critical pairs, that become redundant when adding p:
652 indices = []
653 for i, cp in enumerate(CP):
654 if is_rewritable_or_comparable(cp[0], Num(cp[2]), [p]):
655 indices.append(i)
656 elif is_rewritable_or_comparable(cp[3], Num(cp[5]), [p]):
657 indices.append(i)
659 for i in reversed(indices):
660 del CP[i]
662 # only add new critical pairs that are not made redundant by p:
663 for g in B:
664 if Polyn(g):
665 cp = critical_pair(p, g, ring)
666 if is_rewritable_or_comparable(cp[0], Num(cp[2]), [p]):
667 continue
668 elif is_rewritable_or_comparable(cp[3], Num(cp[5]), [p]):
669 continue
671 CP.append(cp)
673 # sort (other sorting methods/selection strategies were not as successful)
674 CP.sort(key=lambda cp: cp_key(cp, ring), reverse=True)
676 # insert p into B:
677 m = Polyn(p).LM
678 if order(m) <= order(Polyn(B[-1]).LM):
679 B.append(p)
680 else:
681 for i, q in enumerate(B):
682 if order(m) > order(Polyn(q).LM):
683 B.insert(i, p)
684 break
686 k += 1
688 #print(len(B), len(CP), "%d critical pairs removed" % len(indices))
689 else:
690 reductions_to_zero += 1
692 # reduce Groebner basis:
693 H = [Polyn(g).monic() for g in B]
694 H = red_groebner(H, ring)
696 return sorted(H, key=lambda f: order(f.LM), reverse=True)
699def red_groebner(G, ring):
700 """
701 Compute reduced Groebner basis, from BeckerWeispfenning93, p. 216
703 Selects a subset of generators, that already generate the ideal
704 and computes a reduced Groebner basis for them.
705 """
706 def reduction(P):
707 """
708 The actual reduction algorithm.
709 """
710 Q = []
711 for i, p in enumerate(P):
712 h = p.rem(P[:i] + P[i + 1:])
713 if h:
714 Q.append(h)
716 return [p.monic() for p in Q]
718 F = G
719 H = []
721 while F:
722 f0 = F.pop()
724 if not any(monomial_divides(f.LM, f0.LM) for f in F + H):
725 H.append(f0)
727 # Becker, Weispfenning, p. 217: H is Groebner basis of the ideal generated by G.
728 return reduction(H)
731def is_groebner(G, ring):
732 """
733 Check if G is a Groebner basis.
734 """
735 for i in range(len(G)):
736 for j in range(i + 1, len(G)):
737 s = spoly(G[i], G[j], ring)
738 s = s.rem(G)
739 if s:
740 return False
742 return True
745def is_minimal(G, ring):
746 """
747 Checks if G is a minimal Groebner basis.
748 """
749 order = ring.order
750 domain = ring.domain
752 G.sort(key=lambda g: order(g.LM))
754 for i, g in enumerate(G):
755 if g.LC != domain.one:
756 return False
758 for h in G[:i] + G[i + 1:]:
759 if monomial_divides(h.LM, g.LM):
760 return False
762 return True
765def is_reduced(G, ring):
766 """
767 Checks if G is a reduced Groebner basis.
768 """
769 order = ring.order
770 domain = ring.domain
772 G.sort(key=lambda g: order(g.LM))
774 for i, g in enumerate(G):
775 if g.LC != domain.one:
776 return False
778 for term in g.terms():
779 for h in G[:i] + G[i + 1:]:
780 if monomial_divides(h.LM, term[0]):
781 return False
783 return True
785def groebner_lcm(f, g):
786 """
787 Computes LCM of two polynomials using Groebner bases.
789 The LCM is computed as the unique generator of the intersection
790 of the two ideals generated by `f` and `g`. The approach is to
791 compute a Groebner basis with respect to lexicographic ordering
792 of `t*f` and `(1 - t)*g`, where `t` is an unrelated variable and
793 then filtering out the solution that does not contain `t`.
795 References
796 ==========
798 .. [1] [Cox97]_
800 """
801 if f.ring != g.ring:
802 raise ValueError("Values should be equal")
804 ring = f.ring
805 domain = ring.domain
807 if not f or not g:
808 return ring.zero
810 if len(f) <= 1 and len(g) <= 1:
811 monom = monomial_lcm(f.LM, g.LM)
812 coeff = domain.lcm(f.LC, g.LC)
813 return ring.term_new(monom, coeff)
815 fc, f = f.primitive()
816 gc, g = g.primitive()
818 lcm = domain.lcm(fc, gc)
820 f_terms = [ ((1,) + monom, coeff) for monom, coeff in f.terms() ]
821 g_terms = [ ((0,) + monom, coeff) for monom, coeff in g.terms() ] \
822 + [ ((1,) + monom,-coeff) for monom, coeff in g.terms() ]
824 t = Dummy("t")
825 t_ring = ring.clone(symbols=(t,) + ring.symbols, order=lex)
827 F = t_ring.from_terms(f_terms)
828 G = t_ring.from_terms(g_terms)
830 basis = groebner([F, G], t_ring)
832 def is_independent(h, j):
833 return not any(monom[j] for monom in h.monoms())
835 H = [ h for h in basis if is_independent(h, 0) ]
837 h_terms = [ (monom[1:], coeff*lcm) for monom, coeff in H[0].terms() ]
838 h = ring.from_terms(h_terms)
840 return h
842def groebner_gcd(f, g):
843 """Computes GCD of two polynomials using Groebner bases. """
844 if f.ring != g.ring:
845 raise ValueError("Values should be equal")
846 domain = f.ring.domain
848 if not domain.is_Field:
849 fc, f = f.primitive()
850 gc, g = g.primitive()
851 gcd = domain.gcd(fc, gc)
853 H = (f*g).quo([groebner_lcm(f, g)])
855 if len(H) != 1:
856 raise ValueError("Length should be 1")
857 h = H[0]
859 if not domain.is_Field:
860 return gcd*h
861 else:
862 return h.monic()