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
# 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
# Compute a Gröbner basis for this example, which can be
# easily solved to get the critical points
show(I.groebner_basis())
We now compute with the extended critical point system for the Apéry example (in the new direction [2,1,1,1])
Requirements: None
# 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
# Compute a Gröbner basis, and print the elimination poly S(t) in t
GB = I.groebner_basis()
S = GB[-1]
S.factor()
# Define the large factor R(t) of S(t)
R = (S/(t-1)).simplify_full()
R
# This Gröbner basis is not in triangular form
[p.variables() for p in GB]
# 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()
# 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]
# 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]])
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
# 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()
# Now we have parametrized the solutions by u satisfying an explicit polynomial
print([p.variables() for p in GB])
GB[-1].factor()
# 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])])
print([SR(p).leading_coeff(v) for (p,v) in zip(GB[1:6],[w,x,y,z,t])])
# 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]])
We encode the previous example with a Kronecker representation, and observe a decrease in coefficient sizes.
Requirements: None
# 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()
# 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]
# 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])