We run our algorithm and compute the Kronecker representation "by hand" using Sage's capability to compute Gröbner bases. Note that Example 7.1 uses the Maple package of Melczer and Salvy to automatically compute this.
Requirements: None
# We will need our procedure to get the Hessian matrix appearing in asymptotics
# This code copied from Chapter 5
def getHes(H,vars):
dd = len(vars)
V = zero_vector(SR,dd)
U = matrix(SR,dd)
M = matrix(SR,dd-1)
for j in range(dd):
V[j] = 1
for i in range(dd):
U[i,j] = vars[i]*vars[j]*diff(H,vars[i],vars[j])/vars[-1]/diff(H,vars[-1])
for i in range(dd-1):
for j in range(dd-1):
M[i,j] = V[i]*V[j] + U[i,j] - V[j]*U[i,-1] - V[i]*U[j,-1] + V[i]*V[j]*U[-1,-1]
if i == j: M[i,j] = M[i,j] + V[i]
return M
# Set basic info for the Apéry example
var('w,x,y,z,λ,t,u')
vars = [w,x,y,z,λ,t]
H = 1 - w*(1 + x)*(1 + y)*(1 + z)*(x*y*z + y*z + y + z + 1)
# Form the extended critical point system
C = [H] + [v * diff(H,v) - λ for v in vars[0:-2]]
S = C + [H.subs([v==t*v for v in vars])]
S
# Define the multivariate polynomial ring and compute a Gröbner basis
R = PolynomialRing(QQbar,7,'w,x,y,z,λ,t,u',order='lex')
S = S + [u - (w+t)]
GB = R.ideal([R(k) for k in S]).groebner_basis()
# The Kronecker is representation parametrized by the roots of the following polynomial
P = QQ[u](GB[-1])
P
# The polynomial P is square-free: its gcd with its derivative is zero
# (Note: we need to cast between polynomial and symbolic rings to compute desired quantities)
QQ[u](P).gcd(QQ[u](SR(P).diff(u)))
# Compute the Q polynomials for the Kronecker representation from the Gröbner basis
Pd = diff(P,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]
# For instance, here is the first variable in the Kronecker representation
show("In the Kronecker rep. the variable \t",vars[0], " is represented by\t",Q[vars[0]]," divided by\t",Pd)
# All polynomials in the Kronecker representation have degree close to degree(P) = 14
[Q[v].degree() for v in vars]
# They have coefficients that are at most 18 digits long
# As expected, this is an *order of magnitude* smaller than the constants
# appearing in the Gröbner basis, which have constants with ~100 digits
[Q[v].norm(infinity).log(10).round() for v in vars]
max(map(lambda p: abs(QQ(p).numerator()),GB[0].coefficients())).log(10).round()
# Note that P factors
P.factor()
# The critical points, where t=1, are given by the values of u which are roots of
Pt = gcd(P,Pd-Q[t])
Pt
# There is a critical point with positive coordinates (and one with negative coordinates)
for uval in Pt.roots(QQbar, multiplicities=False):
print("There is a critical point where the coordinates {} equal".format(vars[:-2]))
print([p.subs(u == uval) for p in sys[:-2]])
# We can compute the algebraic values of t at solutions of the system
tvals = [sys[-1].subs(u == uval) for uval in P.roots(QQbar, multiplicities=False)]
print("At the solutions of the system, t takes the values {}".format(tvals))
# And we can verify none lie in (0,1)
# This means the point with positive coordinates is *minimal*
num_in_zo = list(filter(lambda c: c.is_real() and c>0 and c<1, tvals))
print("Here are the values lying in (0,1): {}".format(num_in_zo))
# We can simplify the system to the points where t = 1
# First, compute reciprocal of the new derivative modulo the new P
newP = Pt
newPd = diff(newP,u)
_,invPd,_ = Pd.xgcd(newP)
# Then update the Q polynomials using polynomial division
newQ = {}
for v in vars[0:-1]:
_,newQ[v] = QQ[u](Q[v] * newPd * invPd).quo_rem(newP)
# Print the new representation
print("The Kronecker representation for the critical points (where t=1 in the extended system) is given by")
show([v == newQ[v]/newPd for v in vars[0:-2]])
print("where u satisfies {}=0.".format(newP))
# Compute the Hessian determinant Htilde for final asymptotic formula
hes = vars[-3]*diff(H,vars[-3])*getHes(H,vars[0:-2])
hesdet = hes.determinant().factor()
# Add the Hessian, the product of the variables, and the negation of the numerator
# to the Kronecker rep. We use a straightforward but inefficient approach by
# computing another Gröbner basis!
var('T,h,g')
PolySys = [newP, T - w*x*y*z, h - hesdet, g+1] + S
R2 = PolynomialRing(QQ,10,'g,T,h,w,x,y,z,λ,t,u',order='lex')
GB2 = R2.ideal([R2(str(k)) for k in PolySys]).groebner_basis()
for v in [T,h,g]:
[Qv] = list(filter(lambda p: SR(p).has(v), GB2))
_,newQ[v] = QQ[u](newPd * v.subs(solve(SR(Qv),v))).quo_rem(newP)
# Print the new elements of the representation
print("The Kronecker representation for the critical points (where t=1 in the extended system) is given by")
show([v == newQ[v]/newPd for v in [T,h,g]])
print("where u satisfies {}=0.".format(newP))
# Get final asymptotic constants
d = 4
A = newQ[g]/newQ[λ]
B = newQ[λ]^(d-1)/newQ[h]*newPd^(2-d)
C = newPd/newQ[T]
# Find root of P that gives critical point with positive coordinates
[posu] = filter(lambda uval: sys[0].subs(u == uval)>0, newP.roots(QQbar, multiplicities=False))
# Print result of algorithm
var('n')
ASM = ((n*2*pi)^(-3/2) * A * sqrt(B) * SR(C)^n).simplify()
print("Asymptotics are given by evaluating")
show(ASM)
print("At the root {} of the polynomial {}.".format(posu,newP))
print("This gives dominant asymptotic behaviour")
show(ASM.subs(u=posu))
# Because u is the root of a quadratic polynomial, we can express asymptotics in radicals
maxima_calculus('algebraic: true;') # This helps simplify roots
show(ASM.subs(u=posu.radical_expression()).simplify_full())