diff --git a/src/sage/rings/polynomial/toy_d_basis.py b/src/sage/rings/polynomial/toy_d_basis.py index adca6cd8cca..e6745c87d59 100644 --- a/src/sage/rings/polynomial/toy_d_basis.py +++ b/src/sage/rings/polynomial/toy_d_basis.py @@ -6,12 +6,12 @@ .. NOTE:: - The notion of 'term' and 'monomial' in [BW93]_ is swapped from the - notion of those words in Sage (or the other way around, however you - prefer it). In Sage a term is a monomial multiplied by a - coefficient, while in [BW93]_ a monomial is a term multiplied by a - coefficient. Also, what is called LM (the leading monomial) in - Sage is called HT (the head term) in [BW93]_. + The notion of 'term' and 'monomial' in [BW93]_ is swapped from the + notion of those words in Sage (or the other way around, however you + prefer it). In Sage a term is a monomial multiplied by a + coefficient, while in [BW93]_ a monomial is a term multiplied by a + coefficient. Also, what is called LM (the leading monomial) in + Sage is called HT (the head term) in [BW93]_. EXAMPLES:: @@ -76,7 +76,7 @@ We first form a certain ideal `I` in `\ZZ[x, y, z]`, and note that the Groebner basis of `I` over `\QQ` contains 1, so there are no solutions over `\QQ` or an algebraic closure of it (this is not surprising as -there are 4 equations in 3 unknowns).:: +there are 4 equations in 3 unknowns). :: sage: P. = PolynomialRing(IntegerRing(), 3, order='degneglex') sage: I = ideal( x^2 - 3*y, y^3 - x*y, z^3 - x, x^4 - y*z + 1 ) @@ -87,10 +87,7 @@ note that there is a certain integer in the ideal which is not 1:: sage: gb = d_basis(I); gb - [z - 107196348594952664476180297953816049406949517534824683390654620424703403993052759002989622, - y + 84382748470495086324437828161121754084154498572003307352857967748090984550697850484197972764799434672877850291328840, - x + 105754645239745824529618668609551113725317621921665293762587811716173, - 282687803443] + [z ..., y ..., x ..., 282687803443] Now for each prime `p` dividing this integer 282687803443, the Groebner basis of I modulo `p` will be non-trivial and will thus give a solution @@ -126,7 +123,7 @@ def spol(g1, g2): """ - Return S-Polynomial of ``g_1`` and ``g_2``. + Return the S-Polynomial of ``g_1`` and ``g_2``. Let `a_i t_i` be `LT(g_i)`, `b_i = a/a_i` with `a = LCM(a_i,a_j)`, and `s_i = t/t_i` with `t = LCM(t_i,t_j)`. Then the S-Polynomial @@ -134,8 +131,8 @@ def spol(g1, g2): INPUT: - - ``g1`` - polynomial - - ``g2`` - polynomial + - ``g1`` -- polynomial + - ``g2`` -- polynomial EXAMPLES:: @@ -146,20 +143,20 @@ def spol(g1, g2): sage: spol(f,g) x*z - 2*y """ - a1,a2 = g1.lc(),g2.lc() + a1, a2 = g1.lc(), g2.lc() a = a1.lcm(a2) - b1,b2 = a//a1, a//a2 + b1, b2 = a // a1, a // a2 - t1,t2 = g1.lm(), g2.lm() - t = t1.parent().monomial_lcm(t1,t2) - s1,s2 = t//t1, t//t2 + t1, t2 = g1.lm(), g2.lm() + t = t1.parent().monomial_lcm(t1, t2) + s1, s2 = t // t1, t // t2 - return b1*s1*g1 - b2*s2*g2 + return b1 * s1 * g1 - b2 * s2 * g2 def gpol(g1, g2): """ - Return G-Polynomial of ``g_1`` and ``g_2``. + Return the G-Polynomial of ``g_1`` and ``g_2``. Let `a_i t_i` be `LT(g_i)`, `a = a_i*c_i + a_j*c_j` with `a = GCD(a_i,a_j)`, and `s_i = t/t_i` with `t = LCM(t_i,t_j)`. Then the @@ -167,8 +164,8 @@ def gpol(g1, g2): INPUT: - - ``g1`` - polynomial - - ``g2`` - polynomial + - ``g1`` -- polynomial + - ``g2`` -- polynomial EXAMPLES:: @@ -179,18 +176,22 @@ def gpol(g1, g2): sage: gpol(f,g) x^2*y - y """ - a1,a2 = g1.lc(),g2.lc() - a, c1, c2 = xgcd(a1,a2) + a1, a2 = g1.lc(), g2.lc() + a, c1, c2 = xgcd(a1, a2) - t1,t2 = g1.lm(), g2.lm() - t = t1.parent().monomial_lcm(t1,t2) - s1,s2 = t//t1, t//t2 + t1, t2 = g1.lm(), g2.lm() + t = t1.parent().monomial_lcm(t1, t2) + s1, s2 = t // t1, t // t2 - return c1*s1*g1 + c2*s2*g2 + return c1 * s1 * g1 + c2 * s2 * g2 -LM = lambda f: f.lm() -LC = lambda f: f.lc() +def LM(f): + return f.lm() + + +def LC(f): + return f.lc() def d_basis(F, strat=True): @@ -199,8 +200,8 @@ def d_basis(F, strat=True): INPUT: - - ``F`` - an ideal - - ``strat`` - use update strategy (default: ``True``) + - ``F`` -- an ideal + - ``strat`` -- use update strategy (default: ``True``) EXAMPLES:: @@ -214,7 +215,6 @@ def d_basis(F, strat=True): [x - 2020, y - 11313, 22627] """ R = F.ring() - K = R.base_ring() G = set(inter_reduction(F.gens())) B = set((f1, f2) for f1 in G for f2 in G if f1 != f2) @@ -223,46 +223,48 @@ def d_basis(F, strat=True): LCM = R.monomial_lcm divides = R.monomial_divides - divides_ZZ = lambda x, y: ZZ(x).divides(ZZ(y)) - - while B!=set(): - while C!=set(): - f1,f2 = select(C) - C.remove( (f1,f2) ) - lcm_lmf1_lmf2 = LCM(LM(f1),LM(f2) ) - if not any( divides(LM(g), lcm_lmf1_lmf2) and \ - divides_ZZ( LC(g), LC(f1) ) and \ - divides_ZZ( LC(g), LC(f2) ) \ - for g in G): - h = gpol(f1,f2) + + def divides_ZZ(x, y): + return ZZ(x).divides(ZZ(y)) + + while B: + while C: + f1, f2 = select(C) + C.remove((f1, f2)) + lcm_lmf1_lmf2 = LCM(LM(f1), LM(f2)) + if not any(divides(LM(g), lcm_lmf1_lmf2) and + divides_ZZ(LC(g), LC(f1)) and + divides_ZZ(LC(g), LC(f2)) + for g in G): + h = gpol(f1, f2) h0 = h.reduce(G) if h0.lc() < 0: h0 *= -1 if not strat: - D = D.union( [(g,h0) for g in G] ) + D = D.union([(g, h0) for g in G]) G.add(h0) else: - G, D = update(G,D,h0) + G, D = update(G, D, h0) G = inter_reduction(G) - f1,f2 = select(B) - B.remove((f1,f2)) - h = spol(f1,f2) - h0 = h.reduce( G ) + f1, f2 = select(B) + B.remove((f1, f2)) + h = spol(f1, f2) + h0 = h.reduce(G) if h0 != 0: if h0.lc() < 0: - h0 *= -1 + h0 *= -1 if not strat: - D = D.union( [(g,h0) for g in G] ) - G.add( h0 ) + D = D.union([(g, h0) for g in G]) + G.add(h0) else: - G, D = update(G,D,h0) + G, D = update(G, D, h0) B = B.union(D) C = D D = set() - return Sequence(sorted(inter_reduction(G),reverse=True)) + return Sequence(sorted(inter_reduction(G), reverse=True)) def select(P): @@ -271,10 +273,11 @@ def select(P): INPUT: - - ``P`` - a list of critical pairs + - ``P`` -- a list of critical pairs OUTPUT: - an element of P + + an element of P EXAMPLES:: @@ -289,12 +292,12 @@ def select(P): (-2*y - 1, 3*x^2 + 7) """ min_d = 2**20 - min_pair = 0,0 - for fi,fj in sorted(P): - d = fi.parent().monomial_lcm(fi.lm(),fj.lm()).total_degree() + min_pair = 0, 0 + for fi, fj in sorted(P): + d = fi.parent().monomial_lcm(fi.lm(), fj.lm()).total_degree() if d < min_d: min_d = d - min_pair = fi,fj + min_pair = fi, fj return min_pair @@ -308,12 +311,13 @@ def update(G, B, h): INPUT: - - ``G`` - an intermediate Groebner basis - - ``B`` - a list of critical pairs - - ``h`` - a polynomial + - ``G`` -- an intermediate Groebner basis + - ``B`` -- a list of critical pairs + - ``h`` -- a polynomial OUTPUT: - ``G,B`` where ``G`` and ``B`` are updated + + ``G,B`` where ``G`` and ``B`` are updated EXAMPLES:: @@ -331,49 +335,51 @@ def update(G, B, h): R = h.parent() LCM = R.monomial_lcm - lt_divides = lambda x,y: (R.monomial_divides( LM(h), LM(g) ) and LC(h).divides(LC(g)) ) - lt_pairwise_prime = lambda x,y : R.monomial_pairwise_prime(LM(x),LM(y)) and gcd(LC(x),LC(y)) == 1 - lcm_divides = lambda f,g1,h: R.monomial_divides( LCM(LM(h),LM(f[1])), LCM(LM(h),LM(g1)) ) and \ - lcm(LC(h),LC(f[1])).divides(lcm(LC(h),LC(g1))) + def lt_divides(x, y): + return R.monomial_divides(LM(h), LM(g)) and LC(h).divides(LC(g)) - C = set([(h,g) for g in G]) + def lt_pairwise_prime(x, y): + return (R.monomial_pairwise_prime(LM(x), LM(y)) + and gcd(LC(x), LC(y)) == 1) - D = set() - while C != set(): - (h,g1) = C.pop() + def lcm_divides(f, g1, h): + return (R.monomial_divides(LCM(LM(h), LM(f[1])), LCM(LM(h), LM(g1))) + and lcm(LC(h), LC(f[1])).divides(lcm(LC(h), LC(g1)))) + C = set((h, g) for g in G) + + D = set() + while C: + (h, g1) = C.pop() - if lt_pairwise_prime(h,g) or \ - (\ - not any( lcm_divides(f,g1,h) for f in C ) \ - and - not any( lcm_divides(f,g1,h) for f in D ) \ - ): - D.add( (h,g1) ) + if (lt_pairwise_prime(h, g1) or + (not any(lcm_divides(f, g1, h) for f in C) and + not any(lcm_divides(f, g1, h) for f in D))): + D.add((h, g1)) E = set() - while D != set(): - (h,g) = D.pop() - if not lt_pairwise_prime(h,g): - E.add( (h,g) ) + while D: + (h, g) = D.pop() + if not lt_pairwise_prime(h, g): + E.add((h, g)) B_new = set() - while B != set(): - g1,g2 = B.pop() + while B: + g1, g2 = B.pop() - lcm_lg1_lg2 = lcm(LC(g1),LC(g2)) * LCM(LM(g1),LM(g2)) - if not lt_divides(lcm_lg1_lg2, h) or \ - lcm(LC(g1),LC( h)) * R.monomial_lcm(LM(g1),LM( h)) == lcm_lg1_lg2 or \ - lcm(LC( h),LC(g2)) * R.monomial_lcm(LM( h),LM(g2)) == lcm_lg1_lg2 : - B_new.add( (g1,g2) ) + lcm_12 = lcm(LC(g1), LC(g2)) * LCM(LM(g1), LM(g2)) + if (not lt_divides(lcm_12, h) or + lcm(LC(g1), LC(h)) * R.monomial_lcm(LM(g1), LM(h)) == lcm_12 or + lcm(LC(h), LC(g2)) * R.monomial_lcm(LM(h), LM(g2)) == lcm_12): + B_new.add((g1, g2)) - B_new = B_new.union( E ) + B_new = B_new.union(E) G_new = set() - while G != set(): + while G: g = G.pop() - if not lt_divides(g,h): + if not lt_divides(g, h): G_new.add(g) G_new.add(h)