[docs]defgroebner(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``. """ifmethodisNone:method=query('groebner')_groebner_methods={'buchberger':_buchberger,'f5b':_f5b,}try:_groebner=_groebner_methods[method]exceptKeyError:raiseValueError("'%s' is not a valid Groebner bases algorithm (valid are 'buchberger' and 'f5b')"%method)domain,orig=ring.domain,Noneifnotdomain.has_Fieldornotdomain.has_assoc_Field:try:orig,ring=ring,ring.clone(domain=domain.get_field())exceptDomainError:raiseDomainError("can't compute a Groebner basis over %s"%domain)else:seq=[s.set_ring(ring)forsinseq]G=_groebner(seq,ring)iforigisnotNone:G=[g.clear_denoms()[1].set_ring(orig)forginG]returnG

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.orderdomain=ring.domainmonomial_mul=ring.monomial_mulmonomial_div=ring.monomial_divmonomial_lcm=ring.monomial_lcmdefselect(P):# normal selection strategy# select the pair with minimum LCM(LM(f), LM(g))pr=min(P,key=lambdapair:order(monomial_lcm(f[pair[0]].LM,f[pair[1]].LM)))returnprdefnormal(g,J):h=g.rem([f[j]forjinJ])ifnoth:returnNoneelse:h=h.monic()ifnothinI:I[h]=len(f)f.append(h)returnh.LM,I[h]defupdate(G,B,ih):# update G using the set of critical pairs B and h# [BW] page 230h=f[ih]mh=h.LM# filter new pairs (h, g), g in GC=G.copy()D=set()whileC:# select a pair (h, g) by popping an element from Cig=C.pop()g=f[ig]mg=g.LMLCMhg=monomial_lcm(mh,mg)deflcm_divides(ip):# LCM(LM(h), LM(p)) divides LCM(LM(h), LM(g))m=monomial_lcm(mh,f[ip].LM)returnmonomial_div(LCMhg,m)# HT(h) and HT(g) disjoint: mh*mg == LCMhgifmonomial_mul(mh,mg)==LCMhgor(notany(lcm_divides(ipx)foripxinC)andnotany(lcm_divides(pr[1])forprinD)):D.add((ih,ig))E=set()whileD:# select h, g from D (h the same as above)ih,ig=D.pop()mg=f[ig].LMLCMhg=monomial_lcm(mh,mg)ifnotmonomial_mul(mh,mg)==LCMhg:E.add((ih,ig))# filter old pairsB_new=set()whileB:# select g1, g2 from B (-> CP)ig1,ig2=B.pop()mg1=f[ig1].LMmg2=f[ig2].LMLCM12=monomial_lcm(mg1,mg2)# if HT(h) does not divide lcm(HT(g1), HT(g2))ifnotmonomial_div(LCM12,mh)or \
monomial_lcm(mg1,mh)==LCM12or \
monomial_lcm(mg2,mh)==LCM12:B_new.add((ig1,ig2))B_new|=E# filter polynomialsG_new=set()whileG:ig=G.pop()mg=f[ig].LMifnotmonomial_div(mg,mh):G_new.add(ig)G_new.add(ih)returnG_new,B_new# end of update ################################ifnotf:return[]# replace f with a reduced list of initial polynomials; see [BW] page 203f1=f[:]whileTrue:f=f1[:]f1=[]foriinrange(len(f)):p=f[i]r=p.rem(f[:i])ifr:f1.append(r.monic())iff==f1:breakI={}# ip = I[p]; p = f[ip]F=set()# set of indices of polynomialsG=set()# set of indices of intermediate would-be Groebner basisCP=set()# set of pairs of indices of critical pairsfori,hinenumerate(f):I[h]=iF.add(i)###################################### algorithm GROEBNERNEWS2 in [BW] page 232whileF:# select p with minimum monomial according to the monomial orderingh=min([f[x]forxinF],key=lambdaf:order(f.LM))ih=I[h]F.remove(ih)G,CP=update(G,CP,ih)# count the number of critical pairs which reduce to zeroreductions_to_zero=0whileCP:ig1,ig2=select(CP)CP.remove((ig1,ig2))h=spoly(f[ig1],f[ig2],ring)# ordering divisors is on average more efficient [Cox] page 111G1=sorted(G,key=lambdag:order(f[g].LM))ht=normal(h,G1)ifht:G,CP=update(G,CP,ht[1])else:reductions_to_zero+=1####################################### now G is a Groebner basis; reduce itGr=set()foriginG:ht=normal(f[ig],G-set([ig]))ifht:Gr.add(ht[1])Gr=[f[ig]foriginGr]# order according to the monomial orderingGr=sorted(Gr,key=lambdaf:order(f.LM),reverse=True)returnGr

defSign(f):returnf[0]defPolyn(f):returnf[1]defNum(f):returnf[2]defsig(monomial,index):return(monomial,index)deflbp(signature,polynomial,number):return(signature,polynomial,number)# signature functionsdefsig_cmp(u,v,order):""" Compare two signatures by extending the term order to K[X]^n. u < v iff - the index of v is greater than the index of u or - the index of v is equal to the index of u and u[0] < v[0] w.r.t. order u > v otherwise """ifu[1]>v[1]:return-1ifu[1]==v[1]:#if u[0] == v[0]:# return 0iforder(u[0])<order(v[0]):return-1return1defsig_key(s,order):""" Key for comparing two signatures. s = (m, k), t = (n, l) s < t iff [k > l] or [k == l and m < n] s > t otherwise """return(-s[1],order(s[0]))defsig_mult(s,m):""" Multiply a signature by a monomial. The product of a signature (m, i) and a monomial n is defined as (m * t, i). """returnsig(monomial_mul(s[0],m),s[1])# labeled polynomial functionsdeflbp_sub(f,g):""" Subtract labeled polynomial g from f. The signature and number of the difference of f and g are signature and number of the maximum of f and g, w.r.t. lbp_cmp. """ifsig_cmp(Sign(f),Sign(g),Polyn(f).ring.order)<0:max_poly=gelse:max_poly=fret=Polyn(f)-Polyn(g)returnlbp(Sign(max_poly),ret,Num(max_poly))deflbp_mul_term(f,cx):""" Multiply a labeled polynomial with a term. The product of a labeled polynomial (s, p, k) by a monomial is defined as (m * s, m * p, k). """returnlbp(sig_mult(Sign(f),cx[0]),Polyn(f).mul_term(cx),Num(f))deflbp_cmp(f,g):""" Compare two labeled polynomials. f < g iff - Sign(f) < Sign(g) or - Sign(f) == Sign(g) and Num(f) > Num(g) f > g otherwise """ifsig_cmp(Sign(f),Sign(g),Polyn(f).ring.order)==-1:return-1ifSign(f)==Sign(g):ifNum(f)>Num(g):return-1#if Num(f) == Num(g):# return 0return1deflbp_key(f):""" Key for comparing two labeled polynomials. """return(sig_key(Sign(f),Polyn(f).ring.order),-Num(f))# algorithm and helper functionsdefcritical_pair(f,g,ring):""" Compute the critical pair corresponding to two labeled polynomials. A critical pair is a tuple (um, f, vm, g), where um and vm are terms such that um * f - vm * g is the S-polynomial of f and g (so, wlog assume um * f > vm * g). For performance sake, a critical pair is represented as a tuple (Sign(um * f), um, f, Sign(vm * g), vm, g), since um * f creates a new, relatively expensive object in memory, whereas Sign(um * f) and um are lightweight and f (in the tuple) is a reference to an already existing object in memory. """domain=ring.domainltf=Polyn(f).LTltg=Polyn(g).LTlt=(monomial_lcm(ltf[0],ltg[0]),domain.one)um=term_div(lt,ltf,domain)vm=term_div(lt,ltg,domain)# The full information is not needed (now), so only the product# with the leading term is considered:fr=lbp_mul_term(lbp(Sign(f),Polyn(f).leading_term(),Num(f)),um)gr=lbp_mul_term(lbp(Sign(g),Polyn(g).leading_term(),Num(g)),vm)# return in proper order, such that the S-polynomial is just# u_first * f_first - u_second * f_second:iflbp_cmp(fr,gr)==-1:return(Sign(gr),vm,g,Sign(fr),um,f)else:return(Sign(fr),um,f,Sign(gr),vm,g)defcp_cmp(c,d):""" Compare two critical pairs c and d. c < d iff - lbp(c[0], _, Num(c[2]) < lbp(d[0], _, Num(d[2])) (this corresponds to um_c * f_c and um_d * f_d) or - lbp(c[0], _, Num(c[2]) >< lbp(d[0], _, Num(d[2])) and lbp(c[3], _, Num(c[5])) < lbp(d[3], _, Num(d[5])) (this corresponds to vm_c * g_c and vm_d * g_d) c > d otherwise """zero=Polyn(c[2]).ring.zeroc0=lbp(c[0],zero,Num(c[2]))d0=lbp(d[0],zero,Num(d[2]))r=lbp_cmp(c0,d0)ifr==-1:return-1ifr==0:c1=lbp(c[3],zero,Num(c[5]))d1=lbp(d[3],zero,Num(d[5]))r=lbp_cmp(c1,d1)ifr==-1:return-1#if r == 0:# return 0return1defcp_key(c,ring):""" Key for comparing critical pairs. """return(lbp_key(lbp(c[0],ring.zero,Num(c[2]))),lbp_key(lbp(c[3],ring.zero,Num(c[5]))))defs_poly(cp):""" Compute the S-polynomial of a critical pair. The S-polynomial of a critical pair cp is cp[1] * cp[2] - cp[4] * cp[5]. """returnlbp_sub(lbp_mul_term(cp[2],cp[1]),lbp_mul_term(cp[5],cp[4]))defis_rewritable_or_comparable(sign,num,B):""" Check if a labeled polynomial is redundant by checking if its signature and number imply rewritability or comparability. (sign, num) is comparable if there exists a labeled polynomial h in B, such that sign[1] (the index) is less than Sign(h)[1] and sign[0] is divisible by the leading monomial of h. (sign, num) is rewritable if there exists a labeled polynomial h in B, such thatsign[1] is equal to Sign(h)[1], num < Num(h) and sign[0] is divisible by Sign(h)[0]. """forhinB:# comparableifsign[1]<Sign(h)[1]:ifmonomial_divides(Polyn(h).LM,sign[0]):returnTrue# rewritableifsign[1]==Sign(h)[1]:ifnum<Num(h):ifmonomial_divides(Sign(h)[0],sign[0]):returnTruereturnFalsedeff5_reduce(f,B):""" F5-reduce a labeled polynomial f by B. Continously searches for non-zero labeled polynomial h in B, such that the leading term lt_h of h divides the leading term lt_f of f and Sign(lt_h * h) < Sign(f). If such a labeled polynomial h is found, f gets replaced by f - lt_f / lt_h * h. If no such h can be found or f is 0, f is no further F5-reducible and f gets returned. A polynomial that is reducible in the usual sense need not be F5-reducible, e.g.: >>> from sympy.polys.groebnertools import lbp, sig, f5_reduce, Polyn >>> from sympy.polys import ring, QQ, lex >>> R, x,y,z = ring("x,y,z", QQ, lex) >>> f = lbp(sig((1, 1, 1), 4), x, 3) >>> g = lbp(sig((0, 0, 0), 2), x, 2) >>> Polyn(f).rem([Polyn(g)]) 0 >>> f5_reduce(f, [g]) (((1, 1, 1), 4), x, 3) """order=Polyn(f).ring.orderdomain=Polyn(f).ring.domainifnotPolyn(f):returnfwhileTrue:g=fforhinB:ifPolyn(h):ifmonomial_divides(Polyn(h).LM,Polyn(f).LM):t=term_div(Polyn(f).LT,Polyn(h).LT,domain)ifsig_cmp(sig_mult(Sign(h),t[0]),Sign(f),order)<0:# The following check need not be done and is in general slower than without.#if not is_rewritable_or_comparable(Sign(gp), Num(gp), B):hp=lbp_mul_term(h,t)f=lbp_sub(f,hp)breakifg==fornotPolyn(f):returnfdef_f5b(F,ring):""" Computes a reduced Groebner basis for the ideal generated by F. f5b is an implementation of the F5B algorithm by Yao Sun and Dingkang Wang. Similarly to Buchberger's algorithm, the algorithm proceeds by computing critical pairs, computing the S-polynomial, reducing it and adjoining the reduced S-polynomial if it is not 0. Unlike Buchberger's algorithm, each polynomial contains additional information, namely a signature and a number. The signature specifies the path of computation (i.e. from which polynomial in the original basis was it derived and how), the number says when the polynomial was added to the basis. With this information it is (often) possible to decide if an S-polynomial will reduce to 0 and can be discarded. Optimizations include: Reducing the generators before computing a Groebner basis, removing redundant critical pairs when a new polynomial enters the basis and sorting the critical pairs and the current basis. Once a Groebner basis has been found, it gets reduced. ** References ** Yao Sun, Dingkang Wang: "A New Proof for the Correctness of F5 (F5-Like) Algorithm", http://arxiv.org/abs/1004.0084 (specifically v4) Thomas Becker, Volker Weispfenning, Groebner bases: A computational approach to commutative algebra, 1993, p. 203, 216 """order=ring.orderdomain=ring.domain# reduce polynomials (like in Mario Pernici's implementation) (Becker, Weispfenning, p. 203)B=FwhileTrue:F=BB=[]foriinxrange(len(F)):p=F[i]r=p.rem(F[:i])ifr:B.append(r)ifF==B:break# basisB=[lbp(sig(ring.zero_monom,i+1),F[i],i+1)foriinxrange(len(F))]B.sort(key=lambdaf:order(Polyn(f).LM),reverse=True)# critical pairsCP=[critical_pair(B[i],B[j],ring)foriinxrange(len(B))forjinxrange(i+1,len(B))]CP.sort(key=lambdacp:cp_key(cp,ring),reverse=True)k=len(B)reductions_to_zero=0whilelen(CP):cp=CP.pop()# discard redundant critical pairs:ifis_rewritable_or_comparable(cp[0],Num(cp[2]),B):continueifis_rewritable_or_comparable(cp[3],Num(cp[5]),B):continues=s_poly(cp)p=f5_reduce(s,B)p=lbp(Sign(p),Polyn(p).monic(),k+1)ifPolyn(p):# remove old critical pairs, that become redundant when adding p:indices=[]fori,cpinenumerate(CP):ifis_rewritable_or_comparable(cp[0],Num(cp[2]),[p]):indices.append(i)elifis_rewritable_or_comparable(cp[3],Num(cp[5]),[p]):indices.append(i)foriinreversed(indices):delCP[i]# only add new critical pairs that are not made redundant by p:forginB:ifPolyn(g):cp=critical_pair(p,g,ring)ifis_rewritable_or_comparable(cp[0],Num(cp[2]),[p]):continueelifis_rewritable_or_comparable(cp[3],Num(cp[5]),[p]):continueCP.append(cp)# sort (other sorting methods/selection strategies were not as successful)CP.sort(key=lambdacp:cp_key(cp,ring),reverse=True)# insert p into B:m=Polyn(p).LMiforder(m)<=order(Polyn(B[-1]).LM):B.append(p)else:fori,qinenumerate(B):iforder(m)>order(Polyn(q).LM):B.insert(i,p)breakk+=1#print(len(B), len(CP), "%d critical pairs removed" % len(indices))else:reductions_to_zero+=1# reduce Groebner basis:H=[Polyn(g).monic()forginB]H=red_groebner(H,ring)returnsorted(H,key=lambdaf:order(f.LM),reverse=True)

[docs]defis_groebner(G,ring):""" Check if G is a Groebner basis. """foriinxrange(len(G)):forjinxrange(i+1,len(G)):s=spoly(G[i],G[j])s=s.rem(G)ifs:returnFalsereturnTrue

[docs]defis_minimal(G,ring):""" Checks if G is a minimal Groebner basis. """order=ring.orderdomain=ring.domainG.sort(key=lambdag:order(g.LM))fori,ginenumerate(G):ifg.LC!=domain.one:returnFalseforhinG[:i]+G[i+1:]:ifmonomial_divides(h.LM,g.LM):returnFalsereturnTrue

[docs]defis_reduced(G,ring):""" Checks if G is a reduced Groebner basis. """order=ring.orderdomain=ring.domainG.sort(key=lambdag:order(g.LM))fori,ginenumerate(G):ifg.LC!=domain.one:returnFalseforterming:forhinG[:i]+G[i+1:]:ifmonomial_divides(h.LM,term[0]):returnFalsereturnTrue

defgroebner_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]_ """iff.ring!=g.ring:raiseValueError("Values should be equal")ring=f.ringdomain=ring.domainifnotfornotg:returnring.zeroiflen(f)<=1andlen(g)<=1:monom=monomial_lcm(f.LM,g.LM)coeff=domain.lcm(f.LC,g.LC)returnring.term_new(monom,coeff)fc,f=f.primitive()gc,g=g.primitive()lcm=domain.lcm(fc,gc)f_terms=[((1,)+monom,coeff)formonom,coeffinf.terms()]g_terms=[((0,)+monom,coeff)formonom,coeffing.terms()] \
+[((1,)+monom,-coeff)formonom,coeffing.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)defis_independent(h,j):returnall(notmonom[j]formonominh.monoms())H=[hforhinbasisifis_independent(h,0)]h_terms=[(monom[1:],coeff*lcm)formonom,coeffinH[0].terms()]h=ring.from_terms(h_terms)returnhdefgroebner_gcd(f,g):"""Computes GCD of two polynomials using Groebner bases. """iff.ring!=g.ring:raiseValueError("Values should be equal")domain=f.ring.domainifnotdomain.has_Field:fc,f=f.primitive()gc,g=g.primitive()gcd=domain.gcd(fc,gc)H=(f*g).quo([groebner_lcm(f,g)])iflen(H)!=1:raiseValueError("Length should be 1")h=H[0]ifnotdomain.has_Field:returngcd*helse:returnh.monic()