ore_algebra.examples.cbt

Asymptotic enumeration of Compacted Binary Trees

after Genitrini et al.

This example shows how to compute rigorous enclosures of the constant factor in front of the asymptotic expansion of a P-recursive sequence using numerical_transition_matrix().

We consider the example of compacted binary trees of right height at most k, as studied by Genitrini, Gittenberger, Kauers and Wallner (Asymptotic Enumeration of Compacted Binary Trees, arXiv:1703.10031 [math.CO], 2017). Thanks to Michael Wallner for the data and his explanations.

Our goal is to determine a numerical enclosure of the constant κ_{k} appearing in Theorem 3.4 of the paper. Let us take for example k=10. Our starting point is the operator M_{10} given in Theorem 8.6 of the paper, which annihilates the exponential generating series of compacted binary trees of right height at most 10:

sage: from ore_algebra.examples import cbt
sage: dop = cbt.dop[10]
sage: dop
(z^6 - 21*z^5 + 70*z^4 - 84*z^3 + 45*z^2 - 11*z + 1)*Dz^11 + ...

The singularities are:

sage: dop.leading_coefficient().roots(QQbar, multiplicities=False)
[0.2651878342412026?,
 0.3188643842942825?,
 0.4462147547781043?,
 0.7747192223207199?,
 1.988156536964752?,
 17.20685726740094?]

and the dominant one is:

sage: s = _[0]
sage: s
0.2651878342412026?

The origin is an ordinary point, so any solution is characterized by the coefficients in its series expansions of:

sage: dop.local_basis_monomials(0)
[1, z, z^2, z^3, z^4, z^5, z^6, z^7, z^8, z^9, z^10]

while at s, the relevant coefficients are those of:

sage: dop.local_basis_monomials(s)
[1,
 z - 0.2651878342412026?,
 (z - 0.2651878342412026?)^2,
 (z - 0.2651878342412026?)^3,
 (z - 0.2651878342412026?)^4,
 (z - 0.2651878342412026?)^4.260514654474679?,
 (z - 0.2651878342412026?)^5,
 (z - 0.2651878342412026?)^6,
 (z - 0.2651878342412026?)^7,
 (z - 0.2651878342412026?)^8,
 (z - 0.2651878342412026?)^9]

This is the matrix that sends the coefficients of a solution in the basis (1 + O(z^11), z + O(z^11), …) to its coefficients in the basis:

(1 + O((z - 0.265)^10),
  ...,
 (z - 0.265...)^4.26... + + O((z - 0.265)^10),
 ...)

sage: mat = dop.numerical_transition_matrix([0,s], 1e-30) # long time (3.6 s)

We are interested in the fifth row:

sage: mat.row(5) # long time
(..., [4.923936...e-6 +/- ...] + [-5.260478...e-6 +/- ...]*I, ...)

To get the constant corresponding to the particular solution we are interested in, we need to multiply this (row) vector by the (column of) initial values at the origin corresponding to the first 11 coefficients of the generating series of compacted trees (of right height at most 10, but all compacted trees of size at most 10 have right height at most 10):

sage: QQ[['z']](cbt.egf)
1 + z + 3/2*z^2 + 5/2*z^3 + 37/8*z^4 + ...

This is an exponential generating series, but the numerical_transition_matrix() method uses the Taylor coefficients as initial values, not the derivatives, so that we can simply set:

sage: ini = list(cbt.egf)[:dop.order()]

The decomposition in the local basis at the dominant singularity s is hence:

sage: coef = mat*vector(ini) # long time

and we already saw that we were interested in the coefficient of (z - s)^4.260514654474679? + O((z - s)^10), i.e.

sage: coef[5] # long time
[-15159.961154304779924257964877...] + [16196.11585885375838162771522...]*I

Sage’s AsymptoticRing apparently doesn’t digest the algebraic exponent yet, so let’s do the singularity analysis by hand. Slightly annoyingly, the local basis used by the analytic continuation code is defined in terms of (z-s)^α, while extraction of coefficients and application of transfer theorems are easier with expressions of the form (1-z/s)^α. Since:

[z^n] (1-z/s)^α ~ (1/s)^n·n^(-α-1)/Γ(-α),

the constant κ_{10} we are looking for is given by

sage: alpha = QQbar(dop.local_basis_monomials(s)[5].op[1])
sage: C = ComplexBallField(100)
sage: C(-s)^alpha*coef[5]/C(-alpha).gamma() # long time
[645.8284296998659315345812...] + [+/- ...]*I

We glossed over a subtle point, though: in general, we may need to pay attention to the choice of branch of the complex logarithm in the above transformation. What happens in our case is the following. The element of the local basis at s involving (z-s)^α has a branch cut to the left of s. Our analytic continuation path [0,s] arrives at s “above” the branch cut, where log(z-s) = log(s-z) + iπ, and hence:

(z-s)^α = exp(α·i·π)·(s-z)^α = (-1)^α·(s-z)^α.

The rhs, with its branch cut now to the right of s, provides an analytic continuation to a Δ-domain “to the left” of s of the branch of the local solution to which the output of numerical_transition_matrix() refer. This is exactly what we need in the context of singularity analysis, and hence everything should be fine with the above formula. If however we had used, say, an analytic continuation path of the form [0, -i, s], then the correction needed may have been something else than (-1)^α.

With a quick and dirty implementation of some cases of singularity analysis, we can automate the process for other values of k:

sage: def asy(dop, prec=100):
....:     real_sing = dop.leading_coefficient().roots(AA,
....:                                                  multiplicities=False)
....:     s = min(real_sing, key=abs)
....:     ini = list(cbt.egf[:dop.order()])
....:     mat = dop.numerical_transition_matrix([0, s], 10^(-prec))
....:     cmb = mat*vector(ini)
....:     loc = dop.local_basis_monomials(s)
....:     C = cmb[0].parent()
....:     res = []
....:     for pos, mon in enumerate(loc):
....:         if mon.nops() == 0 or mon.operator() is operator.add:
....:             pass
....:         elif mon.operator() is operator.pow:
....:             if mon.op[1] in ZZ:
....:                 pass
....:             else:
....:                 expo = mon.op[1].pyobject()
....:                 res.append(C(-s)^expo*cmb[pos]/C(-expo).gamma())
....:         elif mon.operator() is operator.mul:
....:             assert str(mon.op[1].operator()) == 'log'
....:             assert mon.op[0].operator is operator.pow
....:             expo = mon.op[0,1].pyobject()
....:             res.append(C(-s)^expo*cmb[pos]/C(-expo).gamma())
....:
....:     nonan = [pos for pos, mon in enumerate(loc)
....:              if not (mon.nops() == 0
....:                      or mon.operator() is not operator.pow
....:                      or mon.op[1] in ZZ)]
....:     if len(nonan) != 1:
....:         return CBF(NaN)
....:     pos = nonan[0]
....:     expo = loc[pos].op[1].pyobject()
....:     return (C(-s)^expo*cmb[pos]/C(-expo).gamma()).real()

We obtain:

sage: for k, dop in cbt.dop.items(): # long time (10 s)
....:     print("{}\t{}".format(k, asy(dop, 50)))
2  [0.5613226189564568270393235883810334361992061196...]
3  [0.6049645385653542644762421366594344081594004532...]
4  [0.8724609101215661991266866210277371543438236597...]
5  [1.6248570260792824202355889447707451896274397227...]
6  [3.7818686091669627122667627632166635874791894574...]
7  [10.708084931657092542368755716629129442143758729...]
8  [36.087875288239535234327556576384625828336477172...]
9  [142.21543933025087303695985127830779667104241363...]
10 [645.82842969986593153458120613640308394397361391...]