Example 7.6 (Distribution of Leaves in Planar Trees)

We use Sturm's theorem to verify minimality of a critical point in a parametrized direction.
Requirements: None

# Recall from Example 3.11 in Chapter 3 that the number of binary trees on n nodes with
# k leaves is the coefficient of u^k*y^n*z^n in the power series expansion of the following
var('u,y,z,s,λ,t')
F = (u*y*z - y*z - 2*y + 1)*y/(u*y*z - u*z - y*z - y + 1)
G,H = F.numerator_denominator()
r = [s,1,1]
show(F)

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{{\left(u y z - y z - 2 \, y + 1\right)} y}{u y z - u z - y z - y + 1}$$

# We split the variables into real/imaginary components and define the extended critical point system
var('a1,a2,a3,b1,b2,b3,lambdaR,lambdaI')
assume(a1,'real'); assume(a2,'real'); assume(a3,'real')
assume(b1,'real'); assume(b2,'real'); assume(b3,'real')
Hc = H.subs(u=a1+I*b1, y=a2+I*b2, z=a3+I*b3)
Hr = Hc.real()
Hi = Hc.imag()
CPab1 = [diff(Hr,v)*v + diff(Hr,w)*w - R*lambdaR for (v,w,R) in zip([a1,a2,a3],[b1,b2,b3],r)]
CPab2 = [diff(Hi,v)*v + diff(Hi,w)*w - R*lambdaI for (v,w,R) in zip([a1,a2,a3],[b1,b2,b3],r)]

# Compute a Gröbner basis for the critical point equations in the a and b variables
R = PolynomialRing(QQ,9,'a1,a2,a3,b1,b2,b3,lambdaR,lambdaI,s',order='lex')
J = R.ideal([R(k) for k in [Hr,Hi] + CPab1 + CPab2])
GB = R.ideal([R(k) for k in [Hr,Hi] + CPab1 + CPab2]).groebner_basis()

# Note there are a finite number of (complex) solutions! 
[cp] = solve([SR(k) for k in GB],[a1,a2,a3,b1,b2,b3,lambdaR,lambdaI])
print("The ideal generated by {} has (for fixed s) a single solution:".format(GB))
show(cp)

The ideal generated by [a1*a3 - s, a1*s^2 - 2*a1*s + a1 - s^2, a2 - s, a3*s - s^2 + 2*s - 1, b1, b2, b3, lambdaR - s + 1, lambdaI] has (for fixed s) a single solution:

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left[a_{1} = \frac{s^{2}}{s^{2} - 2 \, s + 1}, a_{2} = s, a_{3} = \frac{s^{2} - 2 \, s + 1}{s}, b_{1} = 0, b_{2} = 0, b_{3} = 0, \mathit{lambdaR} = s - 1, \mathit{lambdaI} = 0\right]$$

# To determine minimality, we add the extra equations using x and y variables
# First, define xj + I*yj as a root of H
# Then encode that the coordinate-wise moduli of xj + I*yj is a multiple of crit pt aj + I*bj
# Finally, add the extra equations derived from critical point method
var('x1,x2,x3,y1,y2,y3,ν')
Hrxy = Hr.subs(a1==x1,a2==x2,a3==x3,b1==y1,b2==y2,b3==y3)
Hixy = Hi.subs(a1==x1,a2==x2,a3==x3,b1==y1,b2==y2,b3==y3)
mod = [X^2+Y^2 - t*(A^2+B^2) for (X,Y,A,B) in zip([x1,x2,x3],[y1,y2,y3],[a1,a2,a3],[b1,b2,b3])]
nuxy = [(Y - ν*X)*diff(Hrxy,X) - (X + ν*Y)*diff(Hrxy,Y) for (X,Y) in zip([x1,x2,x3],[y1,y2,y3])]

# Define the ideal with all these equations
S = PolynomialRing(QQ,17,'a1,a2,a3,b1,b2,b3,x1,x2,x3,y1,y2,y3,lambdaR,lambdaI,ν,s,t',order='lex')
Jxy = S.ideal([S(k) for k in list(GB) + [Hrxy,Hixy] + mod + nuxy])

# Compute the corresponding Gröbner basis. This can take a minute or so.
# (Note: the same Gröbner computation is faster in Maple)
GB2 = Jxy.groebner_basis()

# The values of t at the solutions of the system satisfy the following
tPol = GB2[-1].factor()
show(tPol)

$$\newcommand{\Bold}[1]{\mathbf{#1}}(t - 1) \cdot (s^{4} t^{2} - 3 s^{2} t + 2 s t - t + 1) \cdot (s^{4} t^{3} - s^{4} t^{2} + 4 s^{3} t^{2} - 11 s^{2} t^{2} - s^{2} t + 6 s t^{2} + 6 s t - t^{2} - 2 t - 1) \cdot (s^{4} t^{3} - s^{4} t^{2} + 4 s^{3} t^{2} - 3 s^{2} t^{2} - s^{2} t + 2 s t^{2} + 2 s t - t^{2} - 2 t - 1) \cdot (s^{4} t^{3} - s^{4} t^{2} + 8 s^{3} t^{2} - 11 s^{2} t^{2} + 3 s^{2} t + 6 s t^{2} - 6 s t - t^{2} + 2 t - 1)$$

# To prove minimality for any s in (0,1) we must show that this polynomial
# has no root with t in (0,1). We do this using Sturm sequences, examining
# each factor -- other than t-1, which trivially has no root in (0,1) -- 
# seperately. Details of Sturm sequences are given in the textbook.

# Function to compute a Sturm sequence
def Sturm(poly):
    ST = []
    q = [poly, diff(poly,t)]
    while q[-1] != 0:
        _,r = q[-2].quo_rem(q[-1])
        q += [-r]
    return q[:-1]

# Function to count sign alternations of sequence at a point
def alternations(L,pt):
    alt = 0
    L2 = [p.subs(s=pt) for p in L if p.subs(s=pt) != 0]
    for k in range(len(L2)-1):
        if sign(L2[k]*L2[k+1])<0: alt += 1
    return alt

First Factor

# First factor 
p1 = QQ[s].fraction_field()[t](tPol[1][0])
q = Sturm(p1)
C0 = [k.subs(t=0) for k in q]
C1 = [k.subs(t=1) for k in q]
print("This factor is",p1)
show("Sturm sequence at t=0 is", C0)
show("Sturm sequence at t=1 is", C1)

This factor is s^4*t^2 + (-3*s^2 + 2*s - 1)*t + 1

$$\newcommand{\Bold}[1]{\mathbf{#1}}\verb|Sturm|\phantom{\verb!x!}\verb|sequence|\phantom{\verb!x!}\verb|at|\phantom{\verb!x!}\verb|t=0|\phantom{\verb!x!}\verb|is| \left[1, -3 s^{2} + 2 s - 1, \frac{\frac{5}{4} s^{4} - 3 s^{3} + \frac{5}{2} s^{2} - s + \frac{1}{4}}{s^{4}}\right]$$

$$\newcommand{\Bold}[1]{\mathbf{#1}}\verb|Sturm|\phantom{\verb!x!}\verb|sequence|\phantom{\verb!x!}\verb|at|\phantom{\verb!x!}\verb|t=1|\phantom{\verb!x!}\verb|is| \left[s^{4} - 3 s^{2} + 2 s, 2 s^{4} - 3 s^{2} + 2 s - 1, \frac{\frac{5}{4} s^{4} - 3 s^{3} + \frac{5}{2} s^{2} - s + \frac{1}{4}}{s^{4}}\right]$$

# No roots of a numerator or denominator for either Sturm sequence lie in (0,1)
rts = flatten([QQ[s](SR(k).numerator()).roots(AA, multiplicities=False) for k in C0])
rts += flatten([QQ[s](SR(k).denominator()).roots(AA, multiplicities=False) for k in C0])
rts += flatten([QQ[s](SR(k).numerator()).roots(AA, multiplicities=False) for k in C1])
rts += flatten([QQ[s](SR(k).denominator()).roots(AA, multiplicities=False) for k in C1])
rts01 = list(filter(lambda x: x>0 and x<1, rts))
rts01

[]

# Thus, the signs are the same for all s in (0,1), including both endpoints 
# This means the number of roots of this factor with t in (0,1] for s in (0,1) is
alternations(C1,1/2) - alternations(C0,1/2)

0

Second Factor

p2 = QQ[s].fraction_field()[t](tPol[2][0])
q = Sturm(p2)
C0 = [k.subs(t=0) for k in q]
C1 = [k.subs(t=1) for k in q]
print("This factor is",p2)

This factor is s^4*t^3 + (-s^4 + 4*s^3 - 11*s^2 + 6*s - 1)*t^2 + (-s^2 + 6*s - 2)*t - 1

# The sign conditions of the Sturm sequence can only change at roots of the numerator
# or denominator of its entries -- i.e., for s in (0,1) it can only change at these values
rts = flatten([QQ[s](SR(k).numerator()).roots(AA, multiplicities=False) for k in C0])
rts += flatten([QQ[s](SR(k).denominator()).roots(AA, multiplicities=False) for k in C0])
rts += flatten([QQ[s](SR(k).numerator()).roots(AA, multiplicities=False) for k in C1])
rts += flatten([QQ[s](SR(k).denominator()).roots(AA, multiplicities=False) for k in C1])
rts01 = list(filter(lambda x: x>0 and x<1, rts))
rts01 = sorted(set(rts01))
rts01

[0.3542486889354094?, 0.3851648071345041?, 0.3941724074329110?]

# For s in (0,1) between these points the sign conditions are constant,
# so we can test the sign conditions by picking values of s between each point
# This shows there are no roots for t in (0,1) and s in (0,1) when s not in rts01
pts = [0.1, 0.36, 0.39, 0.5]
[alternations(C1,pt) - alternations(C0,pt) for pt in pts]

[0, 0, 0, 0]

# Manually show there are no roots for t in (0,1) when s in rts01 by solving the poly
[AA[t](SR(p2).subs(s=rt)).roots(multiplicities=False) for rt in rts01]

[[7.140984766613703?], [5.388399446915442?], [5.018578328500075?]]

Third Factor

p3 = QQ[s].fraction_field()[t](tPol[3][0])
q = Sturm(p3)
C0 = [k.subs(t=0) for k in q]
C1 = [k.subs(t=1) for k in q]
print("This factor is",p3)

This factor is s^4*t^3 + (-s^4 + 4*s^3 - 3*s^2 + 2*s - 1)*t^2 + (-s^2 + 2*s - 2)*t - 1

# Find roots in (0,1)
rts = flatten([QQ[s](SR(k).numerator()).roots(AA, multiplicities=False) for k in C0])
rts += flatten([QQ[s](SR(k).denominator()).roots(AA, multiplicities=False) for k in C0])
rts += flatten([QQ[s](SR(k).numerator()).roots(AA, multiplicities=False) for k in C1])
rts += flatten([QQ[s](SR(k).denominator()).roots(AA, multiplicities=False) for k in C1])
rts01 = list(filter(lambda x: x>0 and x<1, rts))
rts01 = sorted(set(rts01))
rts01

[0.7254470427491280?]

# For s in (0,1) on each side of this point the sign conditions are constant,
# so we can test the sign conditions by picking two values of s in (0,1)
# This shows there are no roots for t in (0,1) and s in (0,1) when s not the element of rts01
pts = [0.1,0.9]
[alternations(C1,pt) - alternations(C0,pt) for pt in pts]

[0, 0]

# Manually show there are no roots for t in (0,1) when s in rts01 by solving the poly
[AA[t](SR(p3).subs(s=rt)).roots(multiplicities=False) for rt in rts01]

[[2.148303747649492?]]

Fourth (Final) Factor

# Final factor 
p4 = QQ[s].fraction_field()[t](tPol[4][0])
q = Sturm(p4)
C0 = [k.subs(t=0) for k in q]
C1 = [k.subs(t=1) for k in q]
print("This factor is",p4)

This factor is s^4*t^3 + (-s^4 + 8*s^3 - 11*s^2 + 6*s - 1)*t^2 + (3*s^2 - 6*s + 2)*t - 1

# Find roots in (0,1)
rts = flatten([QQ[s](SR(k).numerator()).roots(AA, multiplicities=False) for k in C0])
rts += flatten([QQ[s](SR(k).denominator()).roots(AA, multiplicities=False) for k in C0])
rts += flatten([QQ[s](SR(k).numerator()).roots(AA, multiplicities=False) for k in C1])
rts += flatten([QQ[s](SR(k).denominator()).roots(AA, multiplicities=False) for k in C1])
rts01 = list(filter(lambda x: x>0 and x<1, rts))
rts01 = sorted(set(rts01))
rts01

[0.2370545200314687?,
 0.2403702750125248?,
 0.2548844570965663?,
 0.2616705580491203?,
 0.3720842206893895?,
 0.4226497308103743?]

# For s in (0,1) on each side of this point the sign conditions are constant,
# so we can test the sign conditions by picking two values of s in (0,1)
# This shows there are no roots for t in (0,1) and s in (0,1) when s not the element of rts01
pts = [0.1, 0.24, 0.25, 0.26, 0.3, 0.4, 0.9]
[alternations(C1,pt) - alternations(C0,pt) for pt in pts]

[0, 0, 0, 0, 0, 0, 0]

# Manually show there are no roots for t in (0,1) when s in rts01 by solving the poly
[AA[t](SR(p4).subs(s=rt)).roots(multiplicities=False) for rt in rts01]

[[1.6630008475?, 13.799307214583437?],
 [1.6745190182?],
 [1.725018869517811?],
 [1.748599118643231?],
 [2.084201751623?],
 [2.170435189413383?]]

Computing Asymptotics

Combining everything above, we have shown that the single critical point is minimal for all $s \in (0,1)$. We end by computing asymptotics after importing our code for smooth ACSV from Chapter 5. If you just want to compute asymptotics assuming minimality, you only need to run the code in this section.

# Import code for smooth asymptotics
# Set a parameter to help simplify some algebraic numbers
maxima_calculus('algebraic: true;')

# Procedure to get Hessian appearing in asymptotics
# Input: H, member of the symbolic ring
#        r, direction vector (which can contain symbolic entries)
#        vars, vector of variables
#        CP, a dictionary mapping elements of vars
# Output: The Hessian H defined in Lemma 5.5 of the textbook at the point w defined by CP
def getHes(H,r,vars,CP):
    dd = len(vars)
    V = zero_vector(SR,dd)
    U = matrix(SR,dd)
    M = matrix(SR,dd-1)

    for j in range(dd):
        V[j] = r[j]/r[-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(CP)

# Procedure to apply differential operator to f and set all variables to zero
# Input: dop, element of a DifferentialWeylAlgebra over a polynomial ring
#        f, an element of the base polynomial ring of dop
# Output: dop(f) evaluated when all variables are zero
def eval_op(dop, f):
    if len(f.parent().gens()) == 1:
        return add([prod([factorial(k) for k in E[0][1]])*E[1]*f[E[0][1][0]] for E in dop])
    else:
        return add([prod([factorial(k) for k in E[0][1]])*E[1]*f[(v for v in E[0][1])] for E in dop])

# Procedure to get critical points of rational function with denominator H, in direction r
# Input: H, member of the symbolic ring
#        r, direction vector (which can contain symbolic entries)
#        vars, vector of variables
# Output: Solutions (if found by solve) of the smooth critical point equations of H in the direction r
def critpt(H,r,vars):
    d = len(vars)
    criteqs = [H] + [r[j]*vars[0]*diff(H,vars[0]) - r[0]*vars[j]*diff(H,vars[j]) for j in range(1,d)]
    return solve(criteqs,vars,solution_dict=true)

# Procedure to compute asymptotic contribution of a strictly minimal contributing point
# Input: G, member of the symbolic ring
#        H, member of the symbolic ring
#        r, direction vector (which can contain symbolic entries)
#        vars, vector of variables
#        CP, a dictionary mapping elements of vars to coordinates of a strictly minimal contributing point
#        M, positive integer describing the number of terms in the asymptotic expansion to compute
#        g, parametrization of variable vars[-1] near CP, in terms of the remaining variables
# Output: ex, pw, se such that ex*pw*(se+O(n^(M-1)) gives an asymptotic expansion of the r-diagonal of 
#         G/H in the variables vars, to order M.
# NOTE: Unlike the textbook, M here refers to the number of terms in the expansion
#       (not the order of the expansion, so M should be at least 1)
def smoothContrib(G,H,r,vars,CP,M,g):
    # Preliminary definitions
    dd = len(vars)
    field = SR
    tvars = list(var('t%d'%i) for i in range(dd-1))
    dvars = list(var('dt%d'%i) for i in range(dd-1))

    # Define differential Weyl algebra and set variable names
    W = DifferentialWeylAlgebra(PolynomialRing(field,tvars))
    WR = W.base_ring()
    T = PolynomialRing(field,tvars).gens()
    D = list(W.differentials())

    # Compute Hessian matrix and differential operator Epsilon
    HES = getHes(H,r,vars,CP)
    HESinv = HES.inverse()
    v = matrix(W,[D[k] for k in range(dd-1)])
    Epsilon = -(v * HESinv.change_ring(W) * v.transpose())[0,0]

    # Define quantities for calculating asymptotics
    tsubs = [v == v.subs(CP)*exp(I*t) for [v,t] in zip(vars,tvars)]
    tsubs += [vars[-1]==g.subs(tsubs)]
    P = (-G/g/diff(H,vars[-1])).subs(tsubs)
    psi = log(g.subs(tsubs)/g.subs(CP)) + I * add([r[k]*tvars[k] for k in range(dd-1)])/r[-1]
    v = matrix(SR,[tvars[k] for k in range(dd-1)])
    psiTilde = psi - (v * HES * v.transpose())[0,0]/2

    # Recursive function to convert symbolic expression to polynomial in t variables
    def to_poly(p,k):
        if k == 0:
            return add([a*T[k]^int(b) for [a,b] in p.coefficients(tvars[k])])
        return add([to_poly(a,k-1)*T[k]^int(b) for [a,b] in p.coefficients(tvars[k])])

    # Compute Taylor expansions to sufficient orders
    N = 2*M
    PsiSeries = to_poly(taylor(psiTilde,*((v,0) for v in tvars), N),dd-2)
    PSeries = to_poly(taylor(P,*((v,0) for v in tvars), N),dd-2)

    # Precompute products used for asymptotics
    EE = [Epsilon^k for k in range(3*M-2)]
    PP = [PSeries] + [0 for k in range(2*M-2)]
    for k in range(1,2*M-1):
        PP[k] = PP[k-1]*PsiSeries
    
    # Function to compute constants appearing in asymptotic expansion
    def Clj(l,j):
        return (-1)^j*SR(eval_op(EE[l+j],PP[l]))/(2^(l+j)*factorial(l)*factorial(l+j))
    
    # Compute different parts of asymptotic expansion
    var('n')
    ex = (prod([1/v^k for (v,k) in zip(vars,r)]).subs(CP).canonicalize_radical())^n
    pw = (r[-1]*n)^((1-dd)/2)
    se = sqrt((2*pi)^(1-dd)/HES.det()) * add([add([Clj(l,j) for l in range(2*j+1)])/(r[-1]*n)^j for j in range(M)])
    
    return ex, pw, se.canonicalize_radical()

# Procedure to aid in printing an asymptotic expansion
# Procedure to get critical points of rational function with denominator H, in direction r
# Input: ex,pw,se as returned by smoothContrib(G,H,r,vars,CP,M,g)
# Output: None (function pretty prints the asymptotic expression defined by ex,pw,se, and M)
def disp_asm(ex,pw,se,M):
    show(ex*pw,LatexExpr("\\Bigg("), se, LatexExpr("+ O\\Bigg("), n^(-M), LatexExpr("\\Bigg)\\Bigg)"))

# Restate basic parameters to make this section independent of above
var('u,y,z,s')
F = (u*y*z - y*z - 2*y + 1)*y/(u*y*z - u*z - y*z - y + 1)
show(F)
G,H = F.numerator_denominator()
r = [s,1,1]
vars = [u,y,z]

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{{\left(u y z - y z - 2 \, y + 1\right)} y}{u y z - u z - y z - y + 1}$$

# Get critical point (shown minimal above)
[CP] = critpt(H,r,vars)

# Get and print asymptotics
M = 2
g = vars[-1].subs(solve(H,vars[-1]))
ex,pw,se = smoothContrib(G,H,r,vars,CP,M,g)
print("The asymptotic expansion for the [s,1,1] diagonal begins")
disp_asm(ex,pw,se,M)

The asymptotic expansion for the [s,1,1] diagonal begins

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{\left(\frac{{\left(s - 1\right)}^{2 \, s}}{{\left(s^{2} - 2 \, s + 1\right)} s^{2 \, s}}\right)^{n}}{n} \Bigg( \frac{1}{2 \, \pi n} + O\Bigg( \frac{1}{n^{2}} \Bigg)\Bigg)$$