## Example 7.7 (Apéry Critical Points)¶

First, we see how to compute critical points for the Apéry example (in the [1,1,1,1] direction) using Gröbner bases.

Requirements: None

In :
# Set basic info for this example
var('w,x,y,z')
vars = [x,y,z,w]
H = 1 - w*(1 + x)*(1 + y)*(1 + z)*(x*y*z + y*z + y + z + 1)

# Form the critical point system
C = [H] + [x * diff(H,x) - v *diff(H,v) for v in vars[1:]]

# Convert this to an ideal in Sage
L = PolynomialRing(QQ,4,'x,y,z,w',order='lex')
I = L.ideal([L(k) for k in C])
I

Out:
Ideal (-x^2*y^2*z^2*w - x^2*y^2*z*w - x^2*y*z^2*w - x^2*y*z*w - 2*x*y^2*z^2*w - 3*x*y^2*z*w - x*y^2*w - 3*x*y*z^2*w - 5*x*y*z*w - 2*x*y*w - x*z^2*w - 2*x*z*w - x*w - y^2*z^2*w - 2*y^2*z*w - y^2*w - 2*y*z^2*w - 4*y*z*w - 2*y*w - z^2*w - 2*z*w - w + 1, -x^2*y*z^2*w - x^2*y*z*w + 2*x*y^2*z^2*w + 3*x*y^2*z*w + x*y^2*w - x*z^2*w - 2*x*z*w - x*w + 2*y^2*z^2*w + 4*y^2*z*w + 2*y^2*w + 2*y*z^2*w + 4*y*z*w + 2*y*w, -x^2*y^2*z*w - x^2*y*z*w + 2*x*y^2*z^2*w - x*y^2*w + 3*x*y*z^2*w - 2*x*y*w + x*z^2*w - x*w + 2*y^2*z^2*w + 2*y^2*z*w + 4*y*z^2*w + 4*y*z*w + 2*z^2*w + 2*z*w, -x^2*y^2*z^2*w - x^2*y^2*z*w - x^2*y*z^2*w - x^2*y*z*w + y^2*z^2*w + 2*y^2*z*w + y^2*w + 2*y*z^2*w + 4*y*z*w + 2*y*w + z^2*w + 2*z*w + w) of Multivariate Polynomial Ring in x, y, z, w over Rational Field
In :
# Compute a Gröbner basis for this example, which can be
# easily solved to get the critical points
show(I.groebner_basis())


## Example 7.8 (An Extended Critical Point System)¶

We now compute with the extended critical point system for the Apéry example (in the new direction [2,1,1,1])

Requirements: None

In :
# Set basic info for this example
var('w,x,y,z,λ,t')
vars = [w,x,y,z]
r = [2,1,1,1]
H = 1 - w*(1 + x)*(1 + y)*(1 + z)*(x*y*z + y*z + y + z + 1)

# Form the extended critical point system
E = [H] + [vars[k] * diff(H,vars[k]) - r[k] * λ for k in range(len(vars))] + [H.subs([v==t*v for v in vars])]

# Convert this to an ideal in Sage
M = PolynomialRing(QQ,6,'λ,w,x,y,z,t',order='lex')
I = M.ideal([M(k) for k in E])
I

Out:
Ideal (-w*x^2*y^2*z^2 - w*x^2*y^2*z - w*x^2*y*z^2 - w*x^2*y*z - 2*w*x*y^2*z^2 - 3*w*x*y^2*z - w*x*y^2 - 3*w*x*y*z^2 - 5*w*x*y*z - 2*w*x*y - w*x*z^2 - 2*w*x*z - w*x - w*y^2*z^2 - 2*w*y^2*z - w*y^2 - 2*w*y*z^2 - 4*w*y*z - 2*w*y - w*z^2 - 2*w*z - w + 1, -2*λ - w*x^2*y^2*z^2 - w*x^2*y^2*z - w*x^2*y*z^2 - w*x^2*y*z - 2*w*x*y^2*z^2 - 3*w*x*y^2*z - w*x*y^2 - 3*w*x*y*z^2 - 5*w*x*y*z - 2*w*x*y - w*x*z^2 - 2*w*x*z - w*x - w*y^2*z^2 - 2*w*y^2*z - w*y^2 - 2*w*y*z^2 - 4*w*y*z - 2*w*y - w*z^2 - 2*w*z - w, -λ - 2*w*x^2*y^2*z^2 - 2*w*x^2*y^2*z - 2*w*x^2*y*z^2 - 2*w*x^2*y*z - 2*w*x*y^2*z^2 - 3*w*x*y^2*z - w*x*y^2 - 3*w*x*y*z^2 - 5*w*x*y*z - 2*w*x*y - w*x*z^2 - 2*w*x*z - w*x, -λ - 2*w*x^2*y^2*z^2 - 2*w*x^2*y^2*z - w*x^2*y*z^2 - w*x^2*y*z - 4*w*x*y^2*z^2 - 6*w*x*y^2*z - 2*w*x*y^2 - 3*w*x*y*z^2 - 5*w*x*y*z - 2*w*x*y - 2*w*y^2*z^2 - 4*w*y^2*z - 2*w*y^2 - 2*w*y*z^2 - 4*w*y*z - 2*w*y, -λ - 2*w*x^2*y^2*z^2 - w*x^2*y^2*z - 2*w*x^2*y*z^2 - w*x^2*y*z - 4*w*x*y^2*z^2 - 3*w*x*y^2*z - 6*w*x*y*z^2 - 5*w*x*y*z - 2*w*x*z^2 - 2*w*x*z - 2*w*y^2*z^2 - 2*w*y^2*z - 4*w*y*z^2 - 4*w*y*z - 2*w*z^2 - 2*w*z, -w*x^2*y^2*z^2*t^7 - w*x^2*y^2*z*t^6 - w*x^2*y*z^2*t^6 - w*x^2*y*z*t^5 - 2*w*x*y^2*z^2*t^6 - 3*w*x*y^2*z*t^5 - w*x*y^2*t^4 - 3*w*x*y*z^2*t^5 - 5*w*x*y*z*t^4 - 2*w*x*y*t^3 - w*x*z^2*t^4 - 2*w*x*z*t^3 - w*x*t^2 - w*y^2*z^2*t^5 - 2*w*y^2*z*t^4 - w*y^2*t^3 - 2*w*y*z^2*t^4 - 4*w*y*z*t^3 - 2*w*y*t^2 - w*z^2*t^3 - 2*w*z*t^2 - w*t + 1) of Multivariate Polynomial Ring in λ, w, x, y, z, t over Rational Field
In :
# Compute a Gröbner basis, and print the elimination poly S(t) in t
GB = I.groebner_basis()
S = GB[-1]
S.factor()

Out:
(1/324) * (t - 1) * (324*t^18 + 1728*t^17 + 3105*t^16 + 9489*t^15 + 8461*t^14 - 66827*t^13 - 84029*t^12 - 419086*t^11 + 680942*t^10 - 1558672*t^9 + 3739659*t^8 + 3302612*t^7 - 5081362*t^6 - 4058824*t^5 + 309952*t^4 + 1599568*t^3 - 1073504*t^2 + 2989728*t + 20736)
In :
# Define the large factor R(t) of S(t)
R = (S/(t-1)).simplify_full()
R

Out:
t^18 + 16/3*t^17 + 115/12*t^16 + 3163/108*t^15 + 8461/324*t^14 - 66827/324*t^13 - 84029/324*t^12 - 209543/162*t^11 + 340471/162*t^10 - 389668/81*t^9 + 1246553/108*t^8 + 825653/81*t^7 - 2540681/162*t^6 - 1014706/81*t^5 + 77488/81*t^4 + 399892/81*t^3 - 268376/81*t^2 + 83048/9*t + 64
In :
# This Gröbner basis is not in triangular form
[p.variables() for p in GB]

Out:
[(λ,), (w, z), (x, z), (y, z), (z,), (z, t), (t,)]
In :
# If we add the polynomial t-1 to the ideal, we get the
# solutions corresponding to critical points
(I + M.ideal([M(t-1)])).groebner_basis()

Out:
[λ + 1/2, w + 249/4*z^2 - 2785/24*z + 1057/36, x + 3/2*z^2 - 3/4*z - 3/4, y - z, z^3 - 7/6*z^2 - 5/6*z + 1/3, t - 1]
In :
# If we add the large factor R(t) we get a more complicated
# system, but it is in triangular form
GB2 = (I + M.ideal([M(R)])).groebner_basis()
[p.variables() for p in GB2]

Out:
[(λ,), (w, t), (x, t), (y, t), (z, t), (t,)]
In :
# We end by printing the degrees and coefficient sizes for the
# polynomials defining x,y, and z in terms of t
def rat_bit_size(p):
L = [abs(k.numerator()) for k in p.coefficients()] + [abs(k.denominator()) for k in p.coefficients()]
return log(max(L),10).round()

print("Polynomial degrees: ",[k.degree() for k in GB2[2:-1]])
print("Max coefficient digit length: ", [rat_bit_size(k) for k in GB2[2:-1]])

Polynomial degrees:  [17, 17, 17]
Max coefficient digit length:  [55, 55, 55]


## Example 7.9 (Large Coefficient Sizes in a Univariate Representation)¶

We encode the extended critical point system for the Apéry example (in the direction [2,1,1,1]) by a linear form u in the variables.

Requirements: None

In :
# Set basic info for this example
var('w,x,y,z,λ,t,u')
vars = [w,x,y,z]
r = [2,1,1,1]
H = 1 - w*(1 + x)*(1 + y)*(1 + z)*(x*y*z + y*z + y + z + 1)

# Form the extended critical point system
E = [H] + [vars[k] * diff(H,vars[k]) - r[k] * λ for k in range(len(vars))] + [H.subs([v==t*v for v in vars])]

# Add the linear form u = t + x
E = E + [u-t-x]

# Convert this to an ideal in Sage
M = PolynomialRing(QQ,7,'λ,w,x,y,z,t,u',order='lex')
I = M.ideal([M(k) for k in E])
GB = I.groebner_basis()

In :
# Now we have parametrized the solutions by u satisfying an explicit polynomial
print([p.variables() for p in GB])
GB[-1].factor()

[(λ,), (w, u), (x, u), (y, u), (z, u), (t, u), (u,)]

Out:
(1/2902376448) * (12*u^3 - 19*u^2 - 10*u + 8) * (241864704*u^18 + 3345795072*u^17 + 21369754368*u^16 + 92838534912*u^15 + 337834871280*u^14 + 1051608329928*u^13 + 2617295811813*u^12 + 4911318707544*u^11 + 5774357518452*u^10 - 503520253323*u^9 - 18653245199288*u^8 - 42468890898735*u^7 - 49451567523618*u^6 - 5425655809061*u^5 + 74101038753291*u^4 + 83418755993667*u^3 + 5275556811046*u^2 - 26574394429140*u - 8505500275728)
In :
# As expected, each of w,x,y,z, and t are defined by linear polynomials in u
print([SR(p).degree(v) for (p,v) in zip(GB[1:6],[w,x,y,z,t])])

[1, 1, 1, 1, 1]
[1, 1, 1, 1, 1]

In :
# We end by printing the degrees and coefficient sizes for the
# polynomials defining x,y,z, and t in terms of u
print("Polynomial degrees: ",[k.degree() for k in GB[2:-1]])
print("Max coefficient digit length: ", [rat_bit_size(k) for k in GB[2:-1]])

Polynomial degrees:  [20, 20, 20, 20]
Max coefficient digit length:  [158, 159, 159, 158]


## Example 7.10 (A Kronecker Representation)¶

We encode the previous example with a Kronecker representation, and observe a decrease in coefficient sizes.

Requirements: None

In :
# Set basic info for this example
var('w,x,y,z,λ,t,u')
vars = [w,x,y,z,λ,t]
r = [2,1,1,1]
H = 1 - w*(1 + x)*(1 + y)*(1 + z)*(x*y*z + y*z + y + z + 1)

# Form the extended critical point system
E = [H] + [vars[k] * diff(H,vars[k]) - r[k] * λ for k in range(4)] + [H.subs([v==t*v for v in vars])]

# Add the linear form u = t + x
E = E + [u-t-x]

# Convert this to an ideal in Sage
M = PolynomialRing(QQ,7,'λ,w,x,y,z,t,u',order='lex')
I = M.ideal([M(k) for k in E])
GB = I.groebner_basis()

In :
# Compute the Q polynomials for the Kronecker representation from the Gröbner basis
P = GB[-1]
Pd = diff(P,M(u))
Q = {}
for v in vars:
[Qv] = list(filter(lambda p: SR(p).has(v), GB))
_,Q[v] = QQ[u](Pd * v.subs(solve(SR(Qv),v))).quo_rem(P)
sys = [SR(Q[v]/Pd) for v in vars]

In :
# The polynomial degrees appearing are the same as the last example,
# but the coefficient lengths are much decreased
print("Polynomial degrees: ",[Q[v].degree() for v in vars])
print("Max coefficient digit length: ", [rat_bit_size(Q[v]) for v in vars])

Polynomial degrees:  [20, 20, 20, 20, 20, 20]
Max coefficient digit length:  [20, 16, 16, 16, 16, 16]

In [ ]: