Skip to content
This repository was archived by the owner on Jan 30, 2023. It is now read-only.

Commit

Permalink
py3: cleanup of toy_d_basis (plus pep corrections)
Browse files Browse the repository at this point in the history
  • Loading branch information
Frédéric Chapoton committed Nov 9, 2018
1 parent 0082a23 commit 05b18cd
Showing 1 changed file with 102 additions and 96 deletions.
198 changes: 102 additions & 96 deletions src/sage/rings/polynomial/toy_d_basis.py
Original file line number Diff line number Diff line change
Expand Up @@ -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::
Expand Down Expand Up @@ -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.<x, y, z> = PolynomialRing(IntegerRing(), 3, order='degneglex')
sage: I = ideal( x^2 - 3*y, y^3 - x*y, z^3 - x, x^4 - y*z + 1 )
Expand All @@ -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
Expand Down Expand Up @@ -126,16 +123,16 @@

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
is defined as: `b_1s_1g_1 - b_2s_2g_2`.
INPUT:
- ``g1`` - polynomial
- ``g2`` - polynomial
- ``g1`` -- polynomial
- ``g2`` -- polynomial
EXAMPLES::
Expand All @@ -146,29 +143,29 @@ 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
G-Polynomial is defined as: `c_1s_1g_1 - c_2s_2g_2`.
INPUT:
- ``g1`` - polynomial
- ``g2`` - polynomial
- ``g1`` -- polynomial
- ``g2`` -- polynomial
EXAMPLES::
Expand All @@ -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):
Expand All @@ -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::
Expand All @@ -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)
Expand All @@ -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):
Expand All @@ -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::
Expand All @@ -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


Expand All @@ -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::
Expand All @@ -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)
Expand Down

0 comments on commit 05b18cd

Please sign in to comment.