Sage companion to Asymptotics of coefficients of algebraic series via embedding into rational series by Torin Greenwood, Stephen Melczer, Tiadora Ruza, and Mark C. Wilson.

The examples listed below are detailed in the supplement to that paper.

Sage Code for Embedding from Proposition 3 of our Paper

var('Y,x,y,z,u,v,w,p,q,r,s,t,a,b,c,p1,p2,p3,q1,q2,q3,s1,s2')
# Y will always be the additional variable.
# P will be the minimal polynomial
# var will be the variable to which Y is attached.
def Furstenberg(P, f, varToAttach, verbose=false, numTerms = 20):
    vars = [var for var in f.variables()]
    varIndex = -1
    for i in range(len(f.variables())):
        if varToAttach == f.variables()[i]:
            varIndex = i
    if varIndex == -1:
        print("The variable to attach isn't a variable of f.")
        return -1
    KK = PolynomialRing(QQ, vars + [Y])
    if not diff(P,Y).subs({a:0 for a in P.variables()}) == 0:
        if P.subs({Y:f}).canonicalize_radical() == 0:
            print("The minimal polynomial is correct.")
            PartialSeries = f.taylor([v for v in f.variables()], 0, 20)
            if PartialSeries.subs({varToAttach:0}) == 0:
                if verbose == true:
                    print("f appears to have no constant terms with respect to ", varToAttach)
                F = SR((diff(P,Y)/P).subs({varToAttach:varToAttach*Y})*Y^2)
                F = simplify(F).factor()
                print("The Furstenberg embedding is:")
                show(F)
#                 for v in vars:
#                     assume(v>0)
                TaylorSeriesFList = list(KK(F.taylor([v for v in F.variables()], 0, 2*numTerms + 2)))
                TaylorSeriesGoalList = list(KK(f.taylor([v for v in f.variables()], 0, numTerms)))
                if verbose==true:
                    print("The Taylor series for f has terms:")
                    print(TaylorSeriesGoalList)
                    print("The Taylor series for F has terms:")
                    print(TaylorSeriesFList)
                count = 0
                for term in TaylorSeriesGoalList:
                    GoalList = list(term[1].exponents()[0])
                    GoalList[-1] = term[1].exponents()[0][varIndex]
                    GoalCount = 0
                    for searchterm in TaylorSeriesFList:
                        if list(searchterm[1].exponents()[0]) == GoalList:
                            if searchterm[0] == term[0]:
                                GoalCount = 1
                            else:
                                GoalCount = -1
                    if GoalCount < 1:
                        count = count + 1
                if count == 0:
                    print("Out of terms that were checked, all terms agree.")
#                     print(latex(F))
                    return F
                else:
                    printf("Some terms disagreed: this embedding is not accurate.")
            else:
                print("subbing in ", varToAttach, "= 0 yields:")
                print(PartialSeries.subs({varToAttach:0}))
                print("f has a non-zero constant term with respect to ", varToAttach, ", so Furstenberg cannot be applied without additional manipulations.")
        else:
            print("The minimal polynomial was incorrect.")
    else:
        print("The derivative of the minimal polynomial is zero.  Check for extraneous factors in the minimal polynomial, try  a substitution, or use Safonov's algorithm.")

Imported Code for Multivariate Asymptotics

# This cell imports the code for smooth asymptotics from Chapter 5
# of An Invitation to Analytic Combinatorics: From One to Several Variables
# by Stephen Melczer (available at melczer.ca/textbook/)

# 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 = [r[j]*vars[0]*diff(H,vars[0]) - r[0]*vars[j]*diff(H,vars[j]) for j in range(1,d)]+[H]
    CritSols = solve(criteqs,vars)
#     The following code factors the solutions before displaying them.
    CPList = []
    for solution in CritSols:
        CPDict = [(vars[0], factor(solution[0].rhs()))]
        for i in range(1,len(vars)):
            CPDict = CPDict + [(vars[i], factor(solution[i].rhs()))]
        CPList = CPList + [dict(CPDict)]
    return CPList

# This searches for critical points that are not smooth.
# The final argument params should include all parameters in the generating function and the direction vector.
# For most of the applications below, the only relevant parameter is p from the direction r = [p, (1-p), p].
def nonsmoothpt(H,r,vars,params=[p]):
    d = len(vars)
    criteqs = [r[j]*vars[0]*diff(H,vars[0]) - r[0]*vars[j]*diff(H,vars[j]) for j in range(1,d)]+[H]+[diff(H,vars[0])]
    return solve(criteqs,vars+params,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)"))

# Test if all coordinates in a list of numbers are real and positive
def pos_coords(L):
    return all(map(lambda c: (c>0) and c.is_real(), L))

# Find critical points with real positive coordinates
def walk_pos_critpt(L,vars):
    return filter(lambda l: pos_coords([v.subs(l) for v in vars]), L)

Example 1: Saunders; analysis of a random tree algorithm

f = -(2*z^2*u^2*p^2+2*z^2*u^2*p*q+2*z^2*p*q-2*z*u^2*p+2*z^2*q^2-2*z*u*q+sqrt(-4*p*q*u^2*z^2+8*p*q*u*z^2-4*p*q*z^2+u^2-2*u+1)+u-1)/(2*(p^2*u^2*z^2+p*q*u^2*z^2+p*q*z^2-p*u^2*z+q^2*z^2-p*u*z-q*u*z-q*z+u))
P = Y^2*z^2*u^2*p^2 + Y^2*z^2*u^2*p*q + 2*Y*z^2*u^2*p^2 + 2*Y*z^2*u^2*p*q - Y^2*z*u^2*p + z^2*u^2*p^2 + Y^2*z^2*p*q + z^2*u^2*p*q + Y^2*z^2*q^2 - Y^2*z*u*p - 2*Y*z*u^2*p - Y^2*z*u*q + 2*Y*z^2*p*q + 2*Y*z^2*q^2 - z*u^2*p - Y^2*z*q - 2*Y*z*u*q + z^2*p*q + z^2*q^2 + Y^2*u + z*u*p - z*u*q + Y*u + z*q - Y
Furstenberg(P, f, z)

The minimal polynomial is correct.
The Furstenberg embedding is:

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{{\left(2 \, Y^{3} p^{2} u^{2} z^{2} + 2 \, Y^{3} p q u^{2} z^{2} + 2 \, Y^{2} p^{2} u^{2} z^{2} + 2 \, Y^{2} p q u^{2} z^{2} + 2 \, Y^{3} p q z^{2} + 2 \, Y^{3} q^{2} z^{2} - 2 \, Y^{2} p u^{2} z + 2 \, Y^{2} p q z^{2} + 2 \, Y^{2} q^{2} z^{2} - 2 \, Y^{2} p u z - 2 \, Y^{2} q u z - 2 \, Y p u^{2} z - 2 \, Y^{2} q z - 2 \, Y q u z + 2 \, Y u + u - 1\right)} Y}{Y^{3} p^{2} u^{2} z^{2} + Y^{3} p q u^{2} z^{2} + 2 \, Y^{2} p^{2} u^{2} z^{2} + 2 \, Y^{2} p q u^{2} z^{2} + Y^{3} p q z^{2} + Y^{3} q^{2} z^{2} + Y p^{2} u^{2} z^{2} + Y p q u^{2} z^{2} - Y^{2} p u^{2} z + 2 \, Y^{2} p q z^{2} + 2 \, Y^{2} q^{2} z^{2} - Y^{2} p u z - Y^{2} q u z - 2 \, Y p u^{2} z + Y p q z^{2} + Y q^{2} z^{2} - Y^{2} q z - 2 \, Y q u z - p u^{2} z + p u z - q u z + Y u + q z + u - 1}$$

Out of terms that were checked, all terms agree.

(2*Y^3*p^2*u^2*z^2 + 2*Y^3*p*q*u^2*z^2 + 2*Y^2*p^2*u^2*z^2 + 2*Y^2*p*q*u^2*z^2 + 2*Y^3*p*q*z^2 + 2*Y^3*q^2*z^2 - 2*Y^2*p*u^2*z + 2*Y^2*p*q*z^2 + 2*Y^2*q^2*z^2 - 2*Y^2*p*u*z - 2*Y^2*q*u*z - 2*Y*p*u^2*z - 2*Y^2*q*z - 2*Y*q*u*z + 2*Y*u + u - 1)*Y/(Y^3*p^2*u^2*z^2 + Y^3*p*q*u^2*z^2 + 2*Y^2*p^2*u^2*z^2 + 2*Y^2*p*q*u^2*z^2 + Y^3*p*q*z^2 + Y^3*q^2*z^2 + Y*p^2*u^2*z^2 + Y*p*q*u^2*z^2 - Y^2*p*u^2*z + 2*Y^2*p*q*z^2 + 2*Y^2*q^2*z^2 - Y^2*p*u*z - Y^2*q*u*z - 2*Y*p*u^2*z + Y*p*q*z^2 + Y*q^2*z^2 - Y^2*q*z - 2*Y*q*u*z - p*u^2*z + p*u*z - q*u*z + Y*u + q*z + u - 1)

# No obvious additive substitutions
assume(u > 0, z > 0)
print(limit(f,u=0))
print(limit(f,z=0))

-1/2*(2*(p*q + q^2)*z^2 + sqrt(-4*p*q*z^2 + 1) - 1)/((p*q + q^2)*z^2 - q*z)
-1/2*(u + sqrt(u^2 - 2*u + 1) - 1)/u

Example 2: Bilateral Schroeder paths

var('x, y, z, u, v, t, Y')
g = 1/(sqrt(1-2*z-4*t*z+z^2)) 
f = (g-1)
P = ((Y+1)^2)*(1-2*z-4*t*z+z^2) - 1
Furstenberg(P, f, z)

The minimal polynomial is correct.
The Furstenberg embedding is:

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{2 \, {\left(Y^{2} z^{2} - 4 \, Y t z - 2 \, Y z + 1\right)} {\left(Y + 1\right)} Y}{Y^{3} z^{2} - 4 \, Y^{2} t z + 2 \, Y^{2} z^{2} - 2 \, Y^{2} z - 8 \, Y t z + Y z^{2} - 4 \, Y z - 4 \, t z + Y - 2 \, z + 2}$$

Out of terms that were checked, all terms agree.

2*(Y^2*z^2 - 4*Y*t*z - 2*Y*z + 1)*(Y + 1)*Y/(Y^3*z^2 - 4*Y^2*t*z + 2*Y^2*z^2 - 2*Y^2*z - 8*Y*t*z + Y*z^2 - 4*Y*z - 4*t*z + Y - 2*z + 2)

# No obvious additive substitution
print(limit(f, t=0)) # A polynomial additive substitution in t is ruled out by Proposition 5 because this limit is not a polynomial
print(limit(f, z=0))

-(z^2 - 2*z - sqrt(z^2 - 2*z + 1) + 1)/(z^2 - 2*z + 1)
0

Example 3: Orthogonal polynomials

g = 1/(sqrt(1-2*x*z+z^2))
f = (g-1)
var('y')
P = ((Y+1)^2)*(1-2*x*z+z^2) - 1
Furstenberg(P, f, z)

The minimal polynomial is correct.
The Furstenberg embedding is:

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{2 \, {\left(Y^{2} z^{2} - 2 \, Y x z + 1\right)} {\left(Y + 1\right)} Y}{Y^{3} z^{2} - 2 \, Y^{2} x z + 2 \, Y^{2} z^{2} - 4 \, Y x z + Y z^{2} - 2 \, x z + Y + 2}$$

Out of terms that were checked, all terms agree.

2*(Y^2*z^2 - 2*Y*x*z + 1)*(Y + 1)*Y/(Y^3*z^2 - 2*Y^2*x*z + 2*Y^2*z^2 - 4*Y*x*z + Y*z^2 - 2*x*z + Y + 2)

# No obvious additive substitution
print(limit(f, x=0))
print(limit(f, z=0))

-(z^2 - sqrt(z^2 + 1) + 1)/(z^2 + 1)
0

Example 4: Callan; parametrized univariate family

g = (1-a*x - sqrt(1 - 2*a*x + (a^2 - 4*b)*x^2))/(2*b*x^2)
f = (g-1)
P = b*x^2*Y^2+2*b*x^2*Y+a*x*Y+b*x^2+a*x-Y
F = Furstenberg(P, f, x)

The minimal polynomial is correct.
The Furstenberg embedding is:

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{{\left(2 \, Y^{3} b x^{2} + 2 \, Y^{2} b x^{2} + Y a x - 1\right)} Y}{Y^{3} b x^{2} + 2 \, Y^{2} b x^{2} + Y b x^{2} + Y a x + a x - 1}$$

Out of terms that were checked, all terms agree.

# To improve runtime, replace parameters with squared versions to eliminate all square roots
var('x,Y,A,B')
a = A^2
b = B^2

g = (1-a*x - sqrt(1 - 2*a*x + (a^2 - 4*b)*x^2))/(2*b*x^2)
f = (g-1)
P = b*x^2*Y^2+2*b*x^2*Y+a*x*Y+b*x^2+a*x-Y
F = Furstenberg(P, f, x)

G = F.numerator()
H = F.denominator()
# Now, we compute critical points:
CritSols = critpt(H, [1,1], [x, Y])
show(CritSols)
CP1 = CritSols[0]
CP2 = CritSols[1]

The minimal polynomial is correct.
The Furstenberg embedding is:

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{{\left(2 \, B^{2} Y^{3} x^{2} + 2 \, B^{2} Y^{2} x^{2} + A^{2} Y x - 1\right)} Y}{B^{2} Y^{3} x^{2} + 2 \, B^{2} Y^{2} x^{2} + B^{2} Y x^{2} + A^{2} Y x + A^{2} x - 1}$$

Out of terms that were checked, all terms agree.

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left[\left\{x : \frac{B}{{\left(A^{2} + 2 \, B\right)} {\left(A^{2} + B\right)}}, Y : \frac{A^{2} + B}{B}\right\}, \left\{x : -\frac{B}{{\left(A^{2} - B\right)} {\left(A^{2} - 2 \, B\right)}}, Y : -\frac{A^{2} - B}{B}\right\}\right]$$

# Only CP1 has positive coordinates, and it is clearly minimal.
# Next, we check for smoothness:
nonsmoothpt(H, [1, 1], [x, Y], [p, A, B])

[]

# Now, we compute asymptotics:
r = [1,1]
vars = [Y,x]

gOptions = solve(H, x)
g1 = gOptions[0].rhs()
g2 = gOptions[1].rhs()
# Check: which solution g1 or g2 actually gives the correct critical point?  The following should be zero for the correct solution (g2 in this case):
# print((g1.subs(y=(sqrt(b)+a)/sqrt(b)) + (3*a*b - (a^2 + 2*b)*sqrt(b))/(a^4 - 5*a^2*b + 4*b^2)).subs(a=1, b=4))
# print((g2.subs(y=(sqrt(b)+a)/sqrt(b)) + (3*a*b - (a^2 + 2*b)*sqrt(b))/(a^4 - 5*a^2*b + 4*b^2)).subs(a=1, b=4))

M=2
ex,pw,se = smoothContrib(G,H,r,vars,CP1,M,g2)
print("The asymptotics:")
disp_asm(ex,pw,se.factor(),M)

The asymptotics:

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

Example 5: Butler et. al; edge flipping in the complete graph

f = (1 - x*(1 + y))/sqrt(1 - 2*x*(1 + y) - x^2*(1 - y)^2) - 1
P = (Y+1)^2*(1-2*x*(1+y) - x^2*(1-y)^2) - (1 - x*(1+y))^2
Furstenberg(P, f, x)

The minimal polynomial is correct.
The Furstenberg embedding is:

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{2 \, {\left(Y^{2} x^{2} y^{2} - 2 \, Y^{2} x^{2} y + Y^{2} x^{2} + 2 \, Y x y + 2 \, Y x - 1\right)} {\left(Y + 1\right)} Y}{Y^{3} x^{2} y^{2} - 2 \, Y^{3} x^{2} y + 2 \, Y^{2} x^{2} y^{2} + Y^{3} x^{2} - 4 \, Y^{2} x^{2} y + 2 \, Y x^{2} y^{2} + 2 \, Y^{2} x^{2} + 2 \, Y^{2} x y + 2 \, Y^{2} x + 2 \, Y x^{2} + 4 \, Y x y + 4 \, Y x - Y - 2}$$

Out of terms that were checked, all terms agree.

2*(Y^2*x^2*y^2 - 2*Y^2*x^2*y + Y^2*x^2 + 2*Y*x*y + 2*Y*x - 1)*(Y + 1)*Y/(Y^3*x^2*y^2 - 2*Y^3*x^2*y + 2*Y^2*x^2*y^2 + Y^3*x^2 - 4*Y^2*x^2*y + 2*Y*x^2*y^2 + 2*Y^2*x^2 + 2*Y^2*x*y + 2*Y^2*x + 2*Y*x^2 + 4*Y*x*y + 4*Y*x - Y - 2)

print(limit(f,x=0))
print(limit(f,y=0))
# No additive polynomial substitution in y can work, but x is not ruled out

0
-(x^2 - sqrt(-x^2 - 2*x + 1)*(x - 1) + 2*x - 1)/(x^2 + 2*x - 1)

Example 6: Bostan et al.; Ising model integrals

f = (1 - 2*w + sqrt((1 - 2*w)^2 - 4*w^2*t^2))/(2*sqrt(1 - t^2)*sqrt((1 - 2*w)^2 - 4*w^2*t^2)) - 1
P = 256*Y^4*t^8*w^4+1024*Y^3*t^8*w^4-1024*Y^4*t^6*w^4+1536*Y^2*t^8*w^4+512*Y^4*t^6*w^3-4096*Y^3*t^6*w^4+1024*Y*t^8*w^4-128*Y^4*t^6*w^2+1536*Y^4*t^4*w^4+2048*Y^3*t^6*w^3-6016*Y^2*t^6*w^4+256*t^8*w^4-1536*Y^4*t^4*w^3-512*Y^3*t^6*w^2+6144*Y^3*t^4*w^4+3072*Y^2*t^6*w^3-3840*Y*t^6*w^4+640*Y^4*t^4*w^2-1024*Y^4*t^2*w^4-6144*Y^3*t^4*w^3-768*Y^2*t^6*w^2+8704*Y^2*t^4*w^4+2048*Y*t^6*w^3-896*t^6*w^4-128*Y^4*t^4*w+1536*Y^4*t^2*w^3+2560*Y^3*t^4*w^2-4096*Y^3*t^2*w^4-8832*Y^2*t^4*w^3-512*Y*t^6*w^2+5120*Y*t^4*w^4+512*t^6*w^3+16*Y^4*t^4-896*Y^4*t^2*w^2+256*Y^4*w^4-512*Y^3*t^4*w+6144*Y^3*t^2*w^3+3744*Y^2*t^4*w^2-5504*Y^2*t^2*w^4-5376*Y*t^4*w^3-128*t^6*w^2+1040*t^4*w^4+256*Y^4*t^2*w-512*Y^4*w^3+64*Y^3*t^4-3584*Y^3*t^2*w^2+1024*Y^3*w^4-768*Y^2*t^4*w+8320*Y^2*t^2*w^3+2368*Y*t^4*w^2-2816*Y*t^2*w^4-1152*t^4*w^3-32*Y^4*t^2+384*Y^4*w^2+1024*Y^3*t^2*w-2048*Y^3*w^3+96*Y^2*t^4-4896*Y^2*t^2*w^2+1280*Y^2*w^4-512*Y*t^4*w+4352*Y*t^2*w^3+544*t^4*w^2-384*t^2*w^4-128*Y^4*w-128*Y^3*t^2+1536*Y^3*w^2+1408*Y^2*t^2*w-2560*Y^2*w^3+64*Y*t^4-2624*Y*t^2*w^2+512*Y*w^4-128*t^4*w+640*t^2*w^3+16*Y^4-512*Y^3*w-176*Y^2*t^2+1920*Y^2*w^2+768*Y*t^2*w-1024*Y*w^3+16*t^4-416*t^2*w^2+64*Y^3-640*Y^2*w-96*Y*t^2+768*Y*w^2+128*t^2*w+80*Y^2-256*Y*w-16*t^2+32*Y
Furstenberg(P, f, t)

The minimal polynomial is correct.
The Furstenberg embedding is:

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{2 \, {\left(8 \, Y^{6} t^{4} w^{2} + 16 \, Y^{5} t^{4} w^{2} + 8 \, Y^{4} t^{4} w^{2} - 16 \, Y^{4} t^{2} w^{2} + 8 \, Y^{4} t^{2} w - 32 \, Y^{3} t^{2} w^{2} - 2 \, Y^{4} t^{2} + 16 \, Y^{3} t^{2} w - 14 \, Y^{2} t^{2} w^{2} - 4 \, Y^{3} t^{2} + 8 \, Y^{2} t^{2} w - 2 \, Y^{2} t^{2} + 8 \, Y^{2} w^{2} - 8 \, Y^{2} w + 16 \, Y w^{2} + 2 \, Y^{2} - 16 \, Y w + 4 \, w^{2} + 4 \, Y - 4 \, w + 1\right)} {\left(2 \, Y t w + 2 \, w - 1\right)} {\left(2 \, Y t w - 2 \, w + 1\right)} {\left(Y t + 1\right)} {\left(Y t - 1\right)} {\left(Y + 1\right)} Y}{16 \, Y^{11} t^{8} w^{4} + 64 \, Y^{10} t^{8} w^{4} + 96 \, Y^{9} t^{8} w^{4} + 64 \, Y^{8} t^{8} w^{4} - 64 \, Y^{9} t^{6} w^{4} + 16 \, Y^{7} t^{8} w^{4} + 32 \, Y^{9} t^{6} w^{3} - 256 \, Y^{8} t^{6} w^{4} - 8 \, Y^{9} t^{6} w^{2} + 128 \, Y^{8} t^{6} w^{3} - 376 \, Y^{7} t^{6} w^{4} - 32 \, Y^{8} t^{6} w^{2} + 192 \, Y^{7} t^{6} w^{3} - 240 \, Y^{6} t^{6} w^{4} - 48 \, Y^{7} t^{6} w^{2} + 128 \, Y^{6} t^{6} w^{3} + 96 \, Y^{7} t^{4} w^{4} - 56 \, Y^{5} t^{6} w^{4} - 32 \, Y^{6} t^{6} w^{2} - 96 \, Y^{7} t^{4} w^{3} + 32 \, Y^{5} t^{6} w^{3} + 384 \, Y^{6} t^{4} w^{4} + 40 \, Y^{7} t^{4} w^{2} - 8 \, Y^{5} t^{6} w^{2} - 384 \, Y^{6} t^{4} w^{3} + 544 \, Y^{5} t^{4} w^{4} - 8 \, Y^{7} t^{4} w + 160 \, Y^{6} t^{4} w^{2} - 552 \, Y^{5} t^{4} w^{3} + 320 \, Y^{4} t^{4} w^{4} + Y^{7} t^{4} - 32 \, Y^{6} t^{4} w + 234 \, Y^{5} t^{4} w^{2} - 336 \, Y^{4} t^{4} w^{3} - 64 \, Y^{5} t^{2} w^{4} + 65 \, Y^{3} t^{4} w^{4} + 4 \, Y^{6} t^{4} - 48 \, Y^{5} t^{4} w + 148 \, Y^{4} t^{4} w^{2} + 96 \, Y^{5} t^{2} w^{3} - 72 \, Y^{3} t^{4} w^{3} - 256 \, Y^{4} t^{2} w^{4} + 6 \, Y^{5} t^{4} - 32 \, Y^{4} t^{4} w - 56 \, Y^{5} t^{2} w^{2} + 34 \, Y^{3} t^{4} w^{2} + 384 \, Y^{4} t^{2} w^{3} - 344 \, Y^{3} t^{2} w^{4} + 4 \, Y^{4} t^{4} + 16 \, Y^{5} t^{2} w - 8 \, Y^{3} t^{4} w - 224 \, Y^{4} t^{2} w^{2} + 520 \, Y^{3} t^{2} w^{3} - 176 \, Y^{2} t^{2} w^{4} - 2 \, Y^{5} t^{2} + Y^{3} t^{4} + 64 \, Y^{4} t^{2} w - 306 \, Y^{3} t^{2} w^{2} + 272 \, Y^{2} t^{2} w^{3} + 16 \, Y^{3} w^{4} - 24 \, Y t^{2} w^{4} - 8 \, Y^{4} t^{2} + 88 \, Y^{3} t^{2} w - 164 \, Y^{2} t^{2} w^{2} - 32 \, Y^{3} w^{3} + 40 \, Y t^{2} w^{3} + 64 \, Y^{2} w^{4} - 11 \, Y^{3} t^{2} + 48 \, Y^{2} t^{2} w + 24 \, Y^{3} w^{2} - 26 \, Y t^{2} w^{2} - 128 \, Y^{2} w^{3} + 80 \, Y w^{4} - 6 \, Y^{2} t^{2} - 8 \, Y^{3} w + 8 \, Y t^{2} w + 96 \, Y^{2} w^{2} - 160 \, Y w^{3} + 32 \, w^{4} + Y^{3} - Y t^{2} - 32 \, Y^{2} w + 120 \, Y w^{2} - 64 \, w^{3} + 4 \, Y^{2} - 40 \, Y w + 48 \, w^{2} + 5 \, Y - 16 \, w + 2}$$

Out of terms that were checked, all terms agree.

2*(8*Y^6*t^4*w^2 + 16*Y^5*t^4*w^2 + 8*Y^4*t^4*w^2 - 16*Y^4*t^2*w^2 + 8*Y^4*t^2*w - 32*Y^3*t^2*w^2 - 2*Y^4*t^2 + 16*Y^3*t^2*w - 14*Y^2*t^2*w^2 - 4*Y^3*t^2 + 8*Y^2*t^2*w - 2*Y^2*t^2 + 8*Y^2*w^2 - 8*Y^2*w + 16*Y*w^2 + 2*Y^2 - 16*Y*w + 4*w^2 + 4*Y - 4*w + 1)*(2*Y*t*w + 2*w - 1)*(2*Y*t*w - 2*w + 1)*(Y*t + 1)*(Y*t - 1)*(Y + 1)*Y/(16*Y^11*t^8*w^4 + 64*Y^10*t^8*w^4 + 96*Y^9*t^8*w^4 + 64*Y^8*t^8*w^4 - 64*Y^9*t^6*w^4 + 16*Y^7*t^8*w^4 + 32*Y^9*t^6*w^3 - 256*Y^8*t^6*w^4 - 8*Y^9*t^6*w^2 + 128*Y^8*t^6*w^3 - 376*Y^7*t^6*w^4 - 32*Y^8*t^6*w^2 + 192*Y^7*t^6*w^3 - 240*Y^6*t^6*w^4 - 48*Y^7*t^6*w^2 + 128*Y^6*t^6*w^3 + 96*Y^7*t^4*w^4 - 56*Y^5*t^6*w^4 - 32*Y^6*t^6*w^2 - 96*Y^7*t^4*w^3 + 32*Y^5*t^6*w^3 + 384*Y^6*t^4*w^4 + 40*Y^7*t^4*w^2 - 8*Y^5*t^6*w^2 - 384*Y^6*t^4*w^3 + 544*Y^5*t^4*w^4 - 8*Y^7*t^4*w + 160*Y^6*t^4*w^2 - 552*Y^5*t^4*w^3 + 320*Y^4*t^4*w^4 + Y^7*t^4 - 32*Y^6*t^4*w + 234*Y^5*t^4*w^2 - 336*Y^4*t^4*w^3 - 64*Y^5*t^2*w^4 + 65*Y^3*t^4*w^4 + 4*Y^6*t^4 - 48*Y^5*t^4*w + 148*Y^4*t^4*w^2 + 96*Y^5*t^2*w^3 - 72*Y^3*t^4*w^3 - 256*Y^4*t^2*w^4 + 6*Y^5*t^4 - 32*Y^4*t^4*w - 56*Y^5*t^2*w^2 + 34*Y^3*t^4*w^2 + 384*Y^4*t^2*w^3 - 344*Y^3*t^2*w^4 + 4*Y^4*t^4 + 16*Y^5*t^2*w - 8*Y^3*t^4*w - 224*Y^4*t^2*w^2 + 520*Y^3*t^2*w^3 - 176*Y^2*t^2*w^4 - 2*Y^5*t^2 + Y^3*t^4 + 64*Y^4*t^2*w - 306*Y^3*t^2*w^2 + 272*Y^2*t^2*w^3 + 16*Y^3*w^4 - 24*Y*t^2*w^4 - 8*Y^4*t^2 + 88*Y^3*t^2*w - 164*Y^2*t^2*w^2 - 32*Y^3*w^3 + 40*Y*t^2*w^3 + 64*Y^2*w^4 - 11*Y^3*t^2 + 48*Y^2*t^2*w + 24*Y^3*w^2 - 26*Y*t^2*w^2 - 128*Y^2*w^3 + 80*Y*w^4 - 6*Y^2*t^2 - 8*Y^3*w + 8*Y*t^2*w + 96*Y^2*w^2 - 160*Y*w^3 + 32*w^4 + Y^3 - Y*t^2 - 32*Y^2*w + 120*Y*w^2 - 64*w^3 + 4*Y^2 - 40*Y*w + 48*w^2 + 5*Y - 16*w + 2)

print(limit(f,w=0))
print(limit(f,t=0))
# No additive polynomial substitution can work (Proposition 5 in paper)

-(t^2 + sqrt(-t^2 + 1) - 1)/(t^2 - 1)
-1/2*(2*w + sqrt(4*w^2 - 4*w + 1) - 1)/(2*w - 1)

Example 7: Cossali; generalization of the Catalan numbers

L = ((1 - z) - sqrt((1-z)^2 - 4*x))/(2*x)
f = L - 1
# multiplicative substitution
fb = f.subs(z=x*z)
P = x*Y^2 + x*Y*z + 2*x*Y + x*z + x - Y
F = Furstenberg(P, fb, x)

The minimal polynomial is correct.
The Furstenberg embedding is:

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

Out of terms that were checked, all terms agree.

# Find critical points:
G = F.numerator()
H = F.denominator()
CritSols = critpt(H, [1, (1-p), 1], [x, z, Y])
show(CritSols)
CP = CritSols[0]

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left[\left\{x : \frac{p^{3}}{{\left(p + 1\right)}^{2}}, z : -\frac{{\left(p + 1\right)} {\left(p - 1\right)}}{p^{2}}, Y : \frac{1}{p}\right\}\right]$$

# Check for non-smooth points:
nonsmoothpt(H, [1, (1-p), 1], [x, z, Y])

[]

# Compute asymptotics:
r = [1-p,1,1]
vars = [z, Y, x]
gOptions = solve(H, x)
g = gOptions[0].rhs()

M=2
ex,pw,se = smoothContrib(G,H,r,vars,CP,M,g)
print("The asymptotics:")
disp_asm(ex,pw,se.factor(),M)

The asymptotics:

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

Example 8: Eu; refinement of noncrossing partitions

P = u*x^2*(Y)^2 + 2*u*x^2*(Y) + u*x^2 + x*(Y)^2 + x*(Y) - (Y)
f = (1 + x - sqrt(1 - 2*x + x^2 - 4*x^2*u))/(2*x*(1 + x*u)) - 1
# Two options for lifting here.  Both are combinatorial!
F = Furstenberg(P, f, x)
Furstenberg(P, f, u)

The minimal polynomial is correct.
The Furstenberg embedding is:

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{{\left(2 \, Y^{3} u x^{2} + 2 \, Y^{2} u x^{2} + 2 \, Y^{2} x + Y x - 1\right)} Y}{Y^{3} u x^{2} + 2 \, Y^{2} u x^{2} + Y u x^{2} + Y^{2} x + Y x - 1}$$

Out of terms that were checked, all terms agree.
The minimal polynomial is correct.
The Furstenberg embedding is:

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

Out of terms that were checked, all terms agree.

(2*Y^2*u*x^2 + 2*Y*u*x^2 + 2*Y*x + x - 1)*Y/(Y^2*u*x^2 + 2*Y*u*x^2 + u*x^2 + Y*x + x - 1)

# For the first embedding, we search for critical points:
G = F.numerator()
H = F.denominator()
CritSols = critpt(H, [p, (1-p), p], [x, u, Y])
show(CritSols)
CP = CritSols[0]

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left[\left\{x : -\frac{{\left(3 \, p - 2\right)} {\left(2 \, p - 1\right)}}{{\left(p - 1\right)} p}, u : \frac{{\left(p - 1\right)}^{2}}{{\left(3 \, p - 2\right)}^{2}}, Y : -\frac{p - 1}{2 \, p - 1}\right\}\right]$$

# Check for non-smooth points:
nonsmoothpt(H, [p, (1-p), p], [x, u, Y])
# Only when p = 0, which is not relevant combinatorially

[{x: 1/2/(4*r1^2 - r1), u: r1, Y: -4*r1, p: 0}]

# Now, we compute asymptotics:
r = [p, p, 1-p]
vars = [Y, x, u]
gOptions = solve(H, u)
g1 = gOptions[0].rhs()

M=2
ex,pw,se = smoothContrib(G,H,r,vars,CP,M,g1)
print("The asymptotics:")
disp_asm(ex,pw,se.factor(),M)

The asymptotics:

$$\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{\left(\frac{{\left(9 \, {\left(p - 1\right)}^{2 \, p} p^{2} - 12 \, {\left(p - 1\right)}^{2 \, p} p + 4 \, {\left(p - 1\right)}^{2 \, p}\right)} p^{p}}{{\left(p^{2} {\left(-3 \, p + 2\right)}^{p} - 2 \, p {\left(-3 \, p + 2\right)}^{p} + {\left(-3 \, p + 2\right)}^{p}\right)} \left(-1\right)^{p} {\left(3 \, p - 2\right)}^{2 \, p}}\right)^{n}}{n {\left(p - 1\right)}} \Bigg( \frac{{\left(p - 1\right)} \sqrt{p}}{2 \, \pi n \sqrt{3 \, p - 2} {\left(2 \, p - 1\right)}^{2}} + O\Bigg( \frac{1}{n^{2}} \Bigg)\Bigg)$$

Example 9: Doslic et al.; RNA secondary structures

var('x, y, Y')
#l is a parameter in the RNA corresponding to the minimum distance between base pairs.  l = 3 is common, although other values are considered.
l = 3
omega = x - 1 - x^2*(1 - pow(x, l))/(1-x)
Omega = (1 - x)*(1-y) - y*omega
S = 1/(2*x^2*y)*(Omega - sqrt(Omega^2 - 4*x^2*y))
f = S - 1
P = -Y*x^4*y+Y^2*x^2*y-Y*x^3*y-x^4*y+Y*x^2*y-x^3*y+Y*x-Y+x
Furstenberg(P, f, x)

The minimal polynomial is correct.
The Furstenberg embedding is:

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{{\left(Y^{4} x^{4} y + Y^{3} x^{3} y - 2 \, Y^{3} x^{2} y - Y^{2} x^{2} y - Y x + 1\right)} Y}{Y^{4} x^{4} y + Y^{3} x^{4} y + Y^{3} x^{3} y - Y^{3} x^{2} y + Y^{2} x^{3} y - Y^{2} x^{2} y - Y x - x + 1}$$

Out of terms that were checked, all terms agree.

(Y^4*x^4*y + Y^3*x^3*y - 2*Y^3*x^2*y - Y^2*x^2*y - Y*x + 1)*Y/(Y^4*x^4*y + Y^3*x^4*y + Y^3*x^3*y - Y^3*x^2*y + Y^2*x^3*y - Y^2*x^2*y - Y*x - x + 1)

# To make combinatorial, we use an additive substitution:
f1 = f - x - x^2
P1 = P.subs(Y = Y + x + x^2)
F = Furstenberg(P1, f1, x)

The minimal polynomial is correct.
The Furstenberg embedding is:

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{{\left(Y^{4} x^{4} y + Y^{3} x^{3} y + 2 \, Y^{3} x^{2} y + Y^{2} x^{2} y + Y x - 1\right)} Y}{Y^{4} x^{4} y + Y^{3} x^{3} y + Y^{3} x^{2} y + Y^{2} x^{3} + Y^{2} x^{2} y + Y x - 1}$$

Out of terms that were checked, all terms agree.

Example 10: Drmota; dissections

f = -(1/2)*(2*x*y^2+x*y+sqrt(x^2*y^2-4*x*y^2-2*x*y+1)-1)/(x*y*(y+1))
P = -Y + x*y^2*(1 + Y)^2 + x*y*(1 + Y)*Y
F = Furstenberg(P, f, x)

The minimal polynomial is correct.
The Furstenberg embedding is:

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

Out of terms that were checked, all terms agree.

# Find critical points:
G = F.numerator()
H = F.denominator()
CritSols = critpt(H, [p, (1-p), p], [x, y, Y])
show(CritSols)
CP = CritSols[0]

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left[\left\{x : -\frac{{\left(3 \, p - 1\right)}^{2} p}{{\left(2 \, p - 1\right)}^{3}}, y : -\frac{{\left(2 \, p - 1\right)}^{2}}{{\left(3 \, p - 1\right)} {\left(p - 1\right)}}, Y : -\frac{2 \, p - 1}{p}\right\}\right]$$

# Check for non-smooth points:
nonsmoothpt(H, [p, (1-p), p], [x, z, Y])

[]

# Compute asymptotics:
r = [(1-p), p, p]
vars = [y, Y, x]
gOptions = solve(H, x)
g1 = gOptions[0].rhs()

M=2
ex,pw,se = smoothContrib(G,H,r,vars,CP,M,g1)
print("The asymptotics:")
disp_asm(ex,pw,se.factor(),M)

The asymptotics:

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

Example 11: Roitner; ascents in Schroeder paths

A = 1 - x*v
B = sqrt(1 - 2*x*(v+2) + x^2*(v-2)^2)
C = 2*x*(1 + x*(v-1))
E = (A - B)/C
f = E - 1
P = Y^2*v*x^2-Y^2*x^2+2*Y*v*x^2+Y^2*x+Y*v*x-2*Y*x^2+v*x^2+2*Y*x+v*x-x^2-Y+x
Furstenberg(P, f, x)

The minimal polynomial is correct.
The Furstenberg embedding is:

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{{\left(2 \, Y^{3} v x^{2} - 2 \, Y^{3} x^{2} + 2 \, Y^{2} v x^{2} - 2 \, Y^{2} x^{2} + 2 \, Y^{2} x + Y v x + 2 \, Y x - 1\right)} Y}{Y^{3} v x^{2} - Y^{3} x^{2} + 2 \, Y^{2} v x^{2} - 2 \, Y^{2} x^{2} + Y v x^{2} + Y^{2} x + Y v x - Y x^{2} + 2 \, Y x + v x + x - 1}$$

Out of terms that were checked, all terms agree.

(2*Y^3*v*x^2 - 2*Y^3*x^2 + 2*Y^2*v*x^2 - 2*Y^2*x^2 + 2*Y^2*x + Y*v*x + 2*Y*x - 1)*Y/(Y^3*v*x^2 - Y^3*x^2 + 2*Y^2*v*x^2 - 2*Y^2*x^2 + Y*v*x^2 + Y^2*x + Y*v*x - Y*x^2 + 2*Y*x + v*x + x - 1)

# No obvious additive substitution
limit(E, v=0) # can't do additive substitution using variable v
limit(E, x=0) # not ruled out with variable x

1

Example 12: Asinowski and Banderier; patterns in lattice paths

f = sqrt((1 + (1 - a*b)*t^2 + (1-a)*(1-b)*t^4)/(1 + (-3-a*b)*t^2 + (1-a)*(1-b)*t^4)) - 1
P = (Y+1)^2*(1 + (-3-a*b)*t^2 + (1-a)*(1-b)*t^4) - (1 + (1 - a*b)*t^2 + (1-a)*(1-b)*t^4)
Furstenberg(P, f, t)

The minimal polynomial is correct.
The Furstenberg embedding is:

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{2 \, {\left(A^{2} B^{2} Y^{4} t^{4} - A^{2} Y^{4} t^{4} - B^{2} Y^{4} t^{4} - A^{2} B^{2} Y^{2} t^{2} + Y^{4} t^{4} - 3 \, Y^{2} t^{2} + 1\right)} {\left(Y + 1\right)} Y}{A^{2} B^{2} Y^{5} t^{4} + 2 \, A^{2} B^{2} Y^{4} t^{4} - A^{2} Y^{5} t^{4} - B^{2} Y^{5} t^{4} - 2 \, A^{2} Y^{4} t^{4} - 2 \, B^{2} Y^{4} t^{4} - A^{2} B^{2} Y^{3} t^{2} + Y^{5} t^{4} - 2 \, A^{2} B^{2} Y^{2} t^{2} + 2 \, Y^{4} t^{4} - 3 \, Y^{3} t^{2} - 6 \, Y^{2} t^{2} - 4 \, Y t^{2} + Y + 2}$$

Out of terms that were checked, all terms agree.

2*(A^2*B^2*Y^4*t^4 - A^2*Y^4*t^4 - B^2*Y^4*t^4 - A^2*B^2*Y^2*t^2 + Y^4*t^4 - 3*Y^2*t^2 + 1)*(Y + 1)*Y/(A^2*B^2*Y^5*t^4 + 2*A^2*B^2*Y^4*t^4 - A^2*Y^5*t^4 - B^2*Y^5*t^4 - 2*A^2*Y^4*t^4 - 2*B^2*Y^4*t^4 - A^2*B^2*Y^3*t^2 + Y^5*t^4 - 2*A^2*B^2*Y^2*t^2 + 2*Y^4*t^4 - 3*Y^3*t^2 - 6*Y^2*t^2 - 4*Y*t^2 + Y + 2)

# No obvious additive substitution
print(limit(f, t=0)) # additive substitution with t not ruled out
print(limit(f, a=0)) # can't do additive substitution with a or b
print(limit(f, b=0))

0
sqrt(A^2*B^2*t^4/(A^2*B^2*t^4 - A^2*B^2*t^2 - A^2*t^4 - B^2*t^4 + t^4 - 3*t^2 + 1) - A^2*B^2*t^2/(A^2*B^2*t^4 - A^2*B^2*t^2 - A^2*t^4 - B^2*t^4 + t^4 - 3*t^2 + 1) - A^2*t^4/(A^2*B^2*t^4 - A^2*B^2*t^2 - A^2*t^4 - B^2*t^4 + t^4 - 3*t^2 + 1) - B^2*t^4/(A^2*B^2*t^4 - A^2*B^2*t^2 - A^2*t^4 - B^2*t^4 + t^4 - 3*t^2 + 1) + t^4/(A^2*B^2*t^4 - A^2*B^2*t^2 - A^2*t^4 - B^2*t^4 + t^4 - 3*t^2 + 1) + t^2/(A^2*B^2*t^4 - A^2*B^2*t^2 - A^2*t^4 - B^2*t^4 + t^4 - 3*t^2 + 1) + 1/(A^2*B^2*t^4 - A^2*B^2*t^2 - A^2*t^4 - B^2*t^4 + t^4 - 3*t^2 + 1)) - 1
sqrt(A^2*B^2*t^4/(A^2*B^2*t^4 - A^2*B^2*t^2 - A^2*t^4 - B^2*t^4 + t^4 - 3*t^2 + 1) - A^2*B^2*t^2/(A^2*B^2*t^4 - A^2*B^2*t^2 - A^2*t^4 - B^2*t^4 + t^4 - 3*t^2 + 1) - A^2*t^4/(A^2*B^2*t^4 - A^2*B^2*t^2 - A^2*t^4 - B^2*t^4 + t^4 - 3*t^2 + 1) - B^2*t^4/(A^2*B^2*t^4 - A^2*B^2*t^2 - A^2*t^4 - B^2*t^4 + t^4 - 3*t^2 + 1) + t^4/(A^2*B^2*t^4 - A^2*B^2*t^2 - A^2*t^4 - B^2*t^4 + t^4 - 3*t^2 + 1) + t^2/(A^2*B^2*t^4 - A^2*B^2*t^2 - A^2*t^4 - B^2*t^4 + t^4 - 3*t^2 + 1) + 1/(A^2*B^2*t^4 - A^2*B^2*t^2 - A^2*t^4 - B^2*t^4 + t^4 - 3*t^2 + 1)) - 1

Example 13: c.f. Flajolet and Sedgewick; patterns in trees

# Here, there is a parameter, m.  We set it to 3 for now:
m = 3
f = 1/(2*z)*(1 - sqrt(1 - 4*z - 4*(u-1)*z^(m+1))) - 1
P = u*z^3+Y^2*z-z^3+2*Y*z-Y+z
Furstenberg(P, f, z)

The minimal polynomial is correct.
The Furstenberg embedding is:

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

Out of terms that were checked, all terms agree.

(2*Y^2*z + 2*Y*z - 1)*Y/(Y^2*u*z^3 - Y^2*z^3 + Y^2*z + 2*Y*z + z - 1)

# Here, there is a parameter, m.  We'll arbitrarily set it to 3 for now:
m = 3
f1 = f - z
P1 = P.subs(Y = Y + z)
Furstenberg(P1, f1, z)

The minimal polynomial is correct.
The Furstenberg embedding is:

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

Out of terms that were checked, all terms agree.

(2*Y^2*z^2 + 2*Y^2*z + 2*Y*z - 1)*Y/(Y^2*u*z^3 + 2*Y^2*z^2 + Y^2*z + 2*Y*z^2 + 2*Y*z - 1)

Example 14: Narayana numbers

f = (1 - x - x*t - sqrt((1-x-x*t)^2 - 4*x^2*t))/(2*x)
P = Y^2*x + Y*x*t + Y*x + x*t - Y
F = Furstenberg(P, f, x)

The minimal polynomial is correct.
The Furstenberg embedding is:

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{{\left(2 \, Y^{2} x + Y t x + Y x - 1\right)} Y}{Y^{2} x + Y t x + Y x + t x - 1}$$

Out of terms that were checked, all terms agree.

# Find critical points:
G = F.numerator()
H = F.denominator()
CritSols = critpt(H, [p, (1-p), p], [x, t, Y])
show(CritSols)
CP = CritSols[0]

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left[\left\{x : -\frac{{\left(2 \, p - 1\right)}^{3}}{{\left(p - 1\right)} p^{2}}, t : \frac{{\left(p - 1\right)}^{2}}{{\left(2 \, p - 1\right)}^{2}}, Y : -\frac{p - 1}{2 \, p - 1}\right\}\right]$$

# This is positive if p > 1/2.
# Search for non-smooth points:
nonsmoothpt(H, [p, (1-p), p], [x, t, Y])

[]

# Compute asymptotics:
r = [(1-p), p, p]
vars = [t, Y, x]
gOptions = solve(H, x)
g1 = gOptions[0].rhs()

M=2
ex,pw,se = smoothContrib(G,H,r,vars,CP,M,g1)
print("The asymptotics:")
disp_asm(ex,pw,se.factor(),M)

The asymptotics:

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

Example 15: Bona and Vince; assembly trees

f = 1 - sqrt((1 - x)^2 + (1 - y)^2 - 1)
P = (Y - 1)^2 - ((1 - x)^2 + (1 - y)^2 - 1)
# There are non-constant terms with respect to both x and y.  So, we attempt a multiplicative substitution:
f1 = f.subs({y:x*y})
P1 = P.subs({y:x*y})
Furstenberg(P1, f1, x)
# Embedding is non-combinatorial

The minimal polynomial is correct.
The Furstenberg embedding is:

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

Out of terms that were checked, all terms agree.

-2*(Y - 1)*Y/(Y*x^2*y^2 + Y*x^2 - 2*x*y - Y - 2*x + 2)

# Instead, we subtract off the first few terms from the series for f:
f2 = f - x - y
P2 = P.subs(Y = Y + x + y)
F = Furstenberg(P2, f2, x)

The minimal polynomial is correct.
The Furstenberg embedding is:

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

Out of terms that were checked, all terms agree.

# Now, we search for critical points:
G = F.numerator()
H = F.denominator()
CritSols = critpt(H, [p, (1-p), p], [x, y, Y])
show(CritSols)
CP1 = CritSols[0]
CP2 = CritSols[1]

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left[\left\{x : -\frac{2 \, p - \sqrt{-4 \, p^{2} + 4 \, p + 1} - 1}{4 \, {\left(p - 1\right)}}, y : -p + \frac{1}{2} \, \sqrt{-4 \, p^{2} + 4 \, p + 1} + \frac{3}{2}, Y : -\sqrt{-4 \, p^{2} + 4 \, p + 1} - 1\right\}, \left\{x : -\frac{2 \, p + \sqrt{-4 \, p^{2} + 4 \, p + 1} - 1}{4 \, {\left(p - 1\right)}}, y : -p - \frac{1}{2} \, \sqrt{-4 \, p^{2} + 4 \, p + 1} + \frac{3}{2}, Y : \sqrt{-4 \, p^{2} + 4 \, p + 1} - 1\right\}\right]$$

# Check for non-smooth points:
nonsmoothpt(H, [p, (1-p), p], [x, z, Y])

[]

# Only CP2 has positive coordinates.  We use it to compute asymptotics.  Note: this block takes about 20 minutes to run on Cocalc.
var('p, x, Y, y')
r = [p,p,1-p]
vars = [x, Y, y]
gOptions = solve(H, y)
g1 = gOptions[0].rhs()

M=2
ex,pw,se = smoothContrib(G,H,r,vars,CP2,M,g1)
print("The asymptotics:")
disp_asm(ex,pw,se.factor(),M)

The asymptotics:

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

Example 16: new and old leaves

f = -(1/2)*(z*y+sqrt(y^2*z^2-4*x*z^2+2*y*z^2-2*y*z+z^2-2*z+1)-z-1)/z - 1
P = ((Y+1) - 1)*(1 - z*((Y+1) - 1 + y)) - z*((Y+1) - 1 + x)
F = Furstenberg(P, f, x)

The minimal polynomial is correct.
The Furstenberg embedding is:

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

Out of terms that were checked, all terms agree.

# To find critical points:
G = F.numerator()
H = F.denominator()
CritSols = critpt(H, [p, q, 1-p-q, p], [x, y, z, Y])
show(CritSols)
CP = CritSols[0]

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left[\left\{x : -\frac{p}{3 \, p + 2 \, q - 1}, y : -\frac{q}{3 \, p + 2 \, q - 1}, z : \frac{3 \, p + 2 \, q - 1}{p + q - 1}, Y : -\frac{p}{3 \, p + 2 \, q - 1}\right\}\right]$$

# To check for non-smooth critical points:
nonsmoothpt(H, [p, q, 1-p-q, p], [x, y, z, Y], [p,q])

[]

# Now, we compute asymptotics:
r = [p, q, 1-p-q, p]
vars = [x, y, z, Y]
gOptions = solve(H, Y)
g1 = gOptions[0].rhs()

M=2
ex,pw,se = smoothContrib(G,H,r,vars,CP,M,g1)
print("The asymptotics:")
disp_asm(ex,pw,se.factor(),M)

The asymptotics:

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

Example 17: Schroeder trees by leaves and vertices

P = (Y - x*y)*(1 - Y) - y*Y^2
f = -(1/2)*(-x*y+sqrt(x^2*y^2-4*x*y^2-2*x*y+1)-1)/(y+1)
Furstenberg(P, f, x)

The minimal polynomial is correct.
The Furstenberg embedding is:

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

Out of terms that were checked, all terms agree.

(Y*x*y - 2*Y*y - 2*Y + 1)*Y/(Y*x*y - Y*y - x*y - Y + 1)

# Additive substitution to make embedding combinatorial:
f2 = f - x*y
P2 = P.subs(Y = Y + x*y)
F = Furstenberg(P2, f2, x)
Furstenberg(P2, f2, y)
# Both embeddings are combinatorial

The minimal polynomial is correct.
The Furstenberg embedding is:

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{{\left(2 \, Y x y^{2} + Y x y + 2 \, Y y + 2 \, Y - 1\right)} Y}{Y x^{2} y^{3} + 2 \, Y x y^{2} + Y x y + Y y + Y - 1}$$

Out of terms that were checked, all terms agree.
The minimal polynomial is correct.
The Furstenberg embedding is:

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{{\left(2 \, Y^{2} x y^{2} + 2 \, Y^{2} y + Y x y + 2 \, Y - 1\right)} Y}{Y^{2} x^{2} y^{3} + 2 \, Y^{2} x y^{2} + Y^{2} y + Y x y + Y - 1}$$

Out of terms that were checked, all terms agree.

(2*Y^2*x*y^2 + 2*Y^2*y + Y*x*y + 2*Y - 1)*Y/(Y^2*x^2*y^3 + 2*Y^2*x*y^2 + Y^2*y + Y*x*y + Y - 1)

# We look for critical points in the first embedding:
G = F.numerator()
H = F.denominator()
CritSols = critpt(H, [p, (1-p), p], [x, y, Y])
show(CritSols)
CP = CritSols[0]

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left[\left\{x : \frac{{\left(3 \, p - 1\right)} {\left(p - 1\right)} p}{{\left(2 \, p - 1\right)}^{3}}, y : -\frac{{\left(2 \, p - 1\right)}^{2}}{{\left(3 \, p - 1\right)} {\left(p - 1\right)}}, Y : \frac{{\left(3 \, p - 1\right)} {\left(2 \, p - 1\right)}}{{\left(p - 1\right)} p}\right\}\right]$$

# Check for non-smooth points:
nonsmoothpt(H, [p, (1-p), p], [x, y, Y])
# Only non-smooth if p = 0, which is an irrelevant case combinatorially.

[{x: -1/2*(2*r2 + 1)/r2^2, y: r2, Y: -4*r2, p: 0},
 {x: r3, y: 0, Y: 1, p: 0},
 {x: 0, y: -1/2, Y: 2, p: 0}]

# Now, we compute asymptotics:
r = [p,1-p,p]
vars = [x, y, Y]
gOptions = solve(H, Y)
g1 = gOptions[0].rhs()

M=2
ex,pw,se = smoothContrib(G,H,r,vars,CP,M,g1)
print("The asymptotics:")
disp_asm(ex,pw,se.factor(),M)

The asymptotics:

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

Example 18: Elizalde; degree of symmetry of lattice paths

f = (1 - x - y - sqrt((1 - x - y)^2 - 4*x*y))/(2*x*y) - 1
P = ((((Y+1)*2*x*y - (1 - x - y))^2 - ((1 - x - y)^2 - 4*x*y))/(4*x*y)).simplify()
# Multiplicative substitution
f = f.subs({y:x*y})
P = P.subs({y:x*y})
F = Furstenberg(P, f, x)

The minimal polynomial is correct.
The Furstenberg embedding is:

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{{\left(2 \, Y^{3} x^{2} y + 2 \, Y^{2} x^{2} y + Y x y + Y x - 1\right)} Y}{Y^{3} x^{2} y + 2 \, Y^{2} x^{2} y + Y x^{2} y + Y x y + Y x + x y + x - 1}$$

Out of terms that were checked, all terms agree.

# We search for critical points:
G = F.numerator()
H = F.denominator()
CritSols = critpt(H, [1, (1-p), 1], [x, y, Y])
show(CritSols)
CP = CritSols[0]

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left[\left\{x : -\frac{{\left(p - 1\right)} p^{3}}{p^{2} - p + 1}, y : \frac{{\left(p - 1\right)}^{2}}{p^{2}}, Y : -\frac{p^{2} - p + 1}{{\left(p - 1\right)} p}\right\}\right]$$

# Check for non-smooth points:
nonsmoothpt(H, [1, (1-p), 1], [x, y, Y])

[]

# Compute asymptotics:
r = [1, 1, 1-p]
vars = [x, Y, y]
gOptions = solve(H, y)
g1 = gOptions[0].rhs()

M=2
ex,pw,se = smoothContrib(G,H,r,vars,CP,M,g1)
print("The asymptotics:")
disp_asm(ex,pw,se.factor(),M)

The asymptotics:

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

Example 19: Bousquet-Melou and Rechnitzer; bar graphs

P = Y - x*y - (x + y + x*y)*Y - x*Y^2
f = -(1/2)*(x*y+sqrt(x^2*y^2-2*x^2*y+2*x*y^2+x^2+y^2-2*x-2*y+1)+x+y-1)/x
F = Furstenberg(P, f, y)

The minimal polynomial is correct.
The Furstenberg embedding is:

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

Out of terms that were checked, all terms agree.

# Following Melczer Chapter 5 Example 5.10, page 26
# See https://melczer.ca/files/TextbookCode/Chapter5/Example5-SmoothASM.html for a more complete analysis of these critical points
G = F.numerator()
H = F.denominator()
CritSols = critpt(H, [p, 1, 1], [x, y, Y])
show(CritSols)
CP1 = CritSols[0]
CP2 = CritSols[1]

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left[\left\{x : \frac{\sqrt{p^{2} + 4} - 2}{p}, y : -\frac{1}{2} \, p + \frac{1}{2} \, \sqrt{p^{2} + 4}, Y : -\frac{1}{2} \, p + \frac{1}{2} \, \sqrt{p^{2} + 4}\right\}, \left\{x : -\frac{\sqrt{p^{2} + 4} + 2}{p}, y : -\frac{1}{2} \, p - \frac{1}{2} \, \sqrt{p^{2} + 4}, Y : -\frac{1}{2} \, p - \frac{1}{2} \, \sqrt{p^{2} + 4}\right\}\right]$$

# Check for non-smooth points:
nonsmoothpt(H, [p, 1, 1], [x, y, Y])

[{x: r4, y: -1, Y: -1, p: 0}]

# Compute asymptotics (takes a long time to run):
r = [p, 1, 1]
vars = [x, Y, y]
gOptions = solve(H, y)
g1 = gOptions[0].rhs()

M=2
ex,pw,se = smoothContrib(G,H,r,vars,CP1,M,g1)
print("The asymptotics:")
disp_asm(ex,pw,se.factor(),M)

The asymptotics:

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

Example 20: Goulden and Jackson; planar maps and the quadratic method

f = -1/6*sqrt(1/3)*sqrt(2*(12*y - 1)*sqrt(-12*y + 1)/(x*y^2) + 81*(x^2*y - x + 1)^2/(x^4*y^2) - 2*(27*x^4*y^2 - x^3*(90*y + 1) + 27*x^2*(4*y + 1) - 54*x + 27)/(x^4*y^2)) + 1/2*(x^2*y - x + 1)/(x^2*y)
P = 27*x^4*y^3*Y^4+(-54*x^4*y^3+54*x^3*y^2-54*x^2*y^2)*Y^3+(27*x^4*y^3-90*x^3*y^2-x^3*y+108*x^2*y^2+27*x^2*y-54*x*y+27*y)*Y^2+(36*x^3*y^2+x^3*y-54*x^2*y^2-36*x^2*y-x^2+90*x*y+x-54*y)*Y+16*x^2*y^2+8*x^2*y+x^2-36*x*y-x+27*y
Furstenberg(P, f, x)

The derivative of the minimal polynomial is zero.  Check for extraneous factors in the minimal polynomial, try  a substitution, or use Safonov's algorithm.

f1 = (f-1)/x - y
P1 = (P.subs(Y=(Y+y)*x + 1)/x^2).full_simplify()
Furstenberg(P1, f1, y)

The minimal polynomial is correct.
The Furstenberg embedding is:

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{{\left(54 \, Y^{4} x^{3} y^{4} + 108 \, Y^{4} x^{3} y^{3} + 54 \, Y^{4} x^{3} y^{2} + 54 \, Y^{3} x^{2} y^{3} + 54 \, Y^{3} x^{2} y^{2} + 54 \, Y^{2} x y^{2} + 54 \, Y^{2} x y - 54 \, Y^{2} y^{2} - 54 \, Y^{2} y + 18 \, Y y - 1\right)} {\left(2 \, Y^{2} x^{3} y^{2} + 2 \, Y^{2} x^{3} y + Y x^{2} y + x - 1\right)} Y}{27 \, Y^{6} x^{6} y^{7} + 108 \, Y^{6} x^{6} y^{6} + 162 \, Y^{6} x^{6} y^{5} + 108 \, Y^{6} x^{6} y^{4} + 54 \, Y^{5} x^{5} y^{6} + 27 \, Y^{6} x^{6} y^{3} + 162 \, Y^{5} x^{5} y^{5} + 162 \, Y^{5} x^{5} y^{4} + 54 \, Y^{5} x^{5} y^{3} + 81 \, Y^{4} x^{4} y^{5} + 216 \, Y^{4} x^{4} y^{4} - 54 \, Y^{4} x^{3} y^{5} + 189 \, Y^{4} x^{4} y^{3} - 162 \, Y^{4} x^{3} y^{4} + 54 \, Y^{4} x^{4} y^{2} - 162 \, Y^{4} x^{3} y^{3} + 72 \, Y^{3} x^{3} y^{4} - 54 \, Y^{4} x^{3} y^{2} + 144 \, Y^{3} x^{3} y^{3} - 54 \, Y^{3} x^{2} y^{4} + 72 \, Y^{3} x^{3} y^{2} - 108 \, Y^{3} x^{2} y^{3} - Y^{2} x^{3} y^{3} - 54 \, Y^{3} x^{2} y^{2} - 2 \, Y^{2} x^{3} y^{2} + 45 \, Y^{2} x^{2} y^{3} - Y^{2} x^{3} y + 72 \, Y^{2} x^{2} y^{2} - 54 \, Y^{2} x y^{3} + 27 \, Y^{2} x^{2} y - 108 \, Y^{2} x y^{2} - Y x^{2} y^{2} + 27 \, Y^{2} y^{3} - 54 \, Y^{2} x y - Y x^{2} y + 54 \, Y^{2} y^{2} + 18 \, Y x y^{2} + 27 \, Y^{2} y + 18 \, Y x y - 2 \, Y y^{2} - 18 \, Y y - x y - x + 1}$$

Out of terms that were checked, all terms agree.

(54*Y^4*x^3*y^4 + 108*Y^4*x^3*y^3 + 54*Y^4*x^3*y^2 + 54*Y^3*x^2*y^3 + 54*Y^3*x^2*y^2 + 54*Y^2*x*y^2 + 54*Y^2*x*y - 54*Y^2*y^2 - 54*Y^2*y + 18*Y*y - 1)*(2*Y^2*x^3*y^2 + 2*Y^2*x^3*y + Y*x^2*y + x - 1)*Y/(27*Y^6*x^6*y^7 + 108*Y^6*x^6*y^6 + 162*Y^6*x^6*y^5 + 108*Y^6*x^6*y^4 + 54*Y^5*x^5*y^6 + 27*Y^6*x^6*y^3 + 162*Y^5*x^5*y^5 + 162*Y^5*x^5*y^4 + 54*Y^5*x^5*y^3 + 81*Y^4*x^4*y^5 + 216*Y^4*x^4*y^4 - 54*Y^4*x^3*y^5 + 189*Y^4*x^4*y^3 - 162*Y^4*x^3*y^4 + 54*Y^4*x^4*y^2 - 162*Y^4*x^3*y^3 + 72*Y^3*x^3*y^4 - 54*Y^4*x^3*y^2 + 144*Y^3*x^3*y^3 - 54*Y^3*x^2*y^4 + 72*Y^3*x^3*y^2 - 108*Y^3*x^2*y^3 - Y^2*x^3*y^3 - 54*Y^3*x^2*y^2 - 2*Y^2*x^3*y^2 + 45*Y^2*x^2*y^3 - Y^2*x^3*y + 72*Y^2*x^2*y^2 - 54*Y^2*x*y^3 + 27*Y^2*x^2*y - 108*Y^2*x*y^2 - Y*x^2*y^2 + 27*Y^2*y^3 - 54*Y^2*x*y - Y*x^2*y + 54*Y^2*y^2 + 18*Y*x*y^2 + 27*Y^2*y + 18*Y*x*y - 2*Y*y^2 - 18*Y*y - x*y - x + 1)