Some computational experiments with the generating function for rooted 3-ary trees.
Requirements: ore_algebra package
Note: See the Maple worksheet corresponding to this example for a better treatment of algebraic functions
# Define the minimal polynomial P(z,y) of the generating function y(z)
var('z y')
p = z*y^3-y+1
p
y^3*z - y + 1
# There are three solutions -- the generating function is the one defined at z=0
parent(p)
Symbolic Ring
As noted in Example 2.8 on binary trees, there is no good built-in method to compute series solutions to bivariate algebraic equations. The Maple gfun package of Salvy and Zimmerman has better support for solving algebraic equations. In Sage the best workaround is to use the fact that any algebraic function is D-finite, and use the ore_algebra package which can compute with D-finite functions.
# Define the ring of linear differential operators with polynomial coefficients
from ore_algebra import *
Pols.<z> = PolynomialRing(QQ); Diff.<Dz> = OreAlgebra(Pols)
# The GF here satisfies (27*z^2 - 4*z)*y''(z) + (54*z - 6)*y'(z) + 6*y(z) = 0
deq = (27*z^2 - 4*z)*Dz^2 + (54*z - 6)*Dz + 6
show(deq)
# A basis of series spanning the solutions to the differential equation at the origin can be computed
# Here there are two basis elements (note the original equation has three solutions!)
deq.generalized_series_solutions(5)
[1 + z + 3*z^2 + 12*z^3 + 55*z^4 + O(z^5), z^(-1/2)*(1 - 3/8*z - 105/128*z^2 - 3003/1024*z^3 - 415701/32768*z^4 + O(z^5))]
# The previous two functions span the C-vector space of solutions
# The generating function is the first basis element
# The other two solutions are linear combinations of the basis elements
# (This can be determined after using the Newton polynom method to compute some initial terms of the solutions)
# First, convert from ore_algebra data structure to PuiseuxSeriesRing
order = 20
[bas1, bas2] = deq.generalized_series_solutions(order)
P.<Z> = PuiseuxSeriesRing(QQ, 'Z')
bas1 = add(bas1[k]*Z^k for k in range(order))
bas2 = Z^(-1/2)*add(bas2[k]*Z^k for k in range(order))
# Define (truncations of) the solutions of the original algebraic equation
sol1, sol2, sol3 = bas1, -bas1/2+bas2, -bas1/2-bas2
print("The generating function expansion begins {}".format(sol1 + O(Z^4)))
print("The first extraneous root expansion begins {}".format(sol2 + O(Z^4)))
print("The second extraneous root expansion begins {}".format(sol3 + O(Z^4)))
The generating function expansion begins 1 + Z + 3*Z^2 + 12*Z^3 + O(Z^4) The first extraneous root expansion begins Z^(-1/2) - 1/2 - 3/8*Z^(1/2) - 1/2*Z - 105/128*Z^(3/2) - 3/2*Z^2 - 3003/1024*Z^(5/2) - 6*Z^3 - 415701/32768*Z^(7/2) + O(Z^4) The second extraneous root expansion begins -Z^(-1/2) - 1/2 + 3/8*Z^(1/2) - 1/2*Z + 105/128*Z^(3/2) - 3/2*Z^2 + 3003/1024*Z^(5/2) - 6*Z^3 + 415701/32768*Z^(7/2) + O(Z^4)
# Verify these truncations satisfy the original polynomial equation
print(QQ['y,z'](p).subs(z=Z, y=sol1) + O(Z^order))
print(QQ['y,z'](p).subs(z=Z, y=sol2) + O(Z^(order-1/2)))
print(QQ['y,z'](p).subs(z=Z, y=sol3) + O(Z^(order-1/2)))
O(Z^20) O(Z^(39/2)) O(Z^(39/2))
# The coefficients of an algebraic function satisfy a P-recurrence
# This computation shows (4*n^2 + 10*n + 6)*c_{n+1} = (27*n^2 + 27*n + 6)*c_n
Ind.<n> = PolynomialRing(QQ); Shift.<Sn> = OreAlgebra(Ind)
rec = deq.to_S(Shift)
rec
(-4*n^2 - 10*n - 6)*Sn + 27*n^2 + 27*n + 6
# The recurrence can be used to compute many terms of the sequence efficiently
# We must give the initial term in the sequence (to specify which solution we care about)
rec.to_list([1],100)
[1, 1, 3, 12, 55, 273, 1428, 7752, 43263, 246675, 1430715, 8414640, 50067108, 300830572, 1822766520, 11124755664, 68328754959, 422030545335, 2619631042665, 16332922290300, 102240109897695, 642312451217745, 4048514844039120, 25594403741131680, 162250238001816900, 1031147983159782228, 6568517413771094628, 41932353590942745504, 268225186597703313816, 1718929965542850284040, 11034966795189838872624, 70956023048640039202464, 456949965738717944767791, 2946924270225408943665279, 19030649059639789214206725, 123052100237542105872786180, 796607831560617902288322405, 5162879946168545215371343587, 33496962712940417760973884708, 217550867863011281855594752680, 1414282077098335379544565517191, 9202600068524372703278082352971, 59932899605936040714626166584475, 390645234961546222075767026462400, 2548271840422037344975860237738000, 16635641296703864308483436321233200, 108679967966404768511803431114031200, 710496811856518520237299474629019200, 4647985909007237458743106679064711300, 30426054945480277365983787382745806500, 199293672373583488061784498821082334140, 1306164582615216509992597891759417970640, 8565425435995939036682228482499260153620, 56200126088515763642375073109736372959980, 368938646289542831658405169038432458575200, 2423210784425567315944099994775276777059520, 15923573648667567272266673527613832841413720, 104687476184489521220187993237607647032909880, 688567758918141068396974784875829607672693720, 4530954345822264258644101276351421147198184000, 29827534579873756985313709683286979461162921200, 196437153964454478730548625427197056659262085200, 1294204010506457507610556492169378285551215751008, 8530003035876456718702048001257342146013974679872, 56241607952753745334610705298987817310837053258095, 370956832348855971723751149936236380931126414813199, 2447605128273073612473937013215645786143652277643069, 16155005655931729647837925891539419948344504253763540, 106663615595824898981200406146704399055680776314377585, 704473645075333982479839404375246601833187235495177575, 4654236823512905605830331092249349585849895972864188380, 30758394797997531797498626576159446322402868116541053080, 203332434576992304310763190915100018001976431327881300605, 1344540697525747546508940943448381528040466460776733688345, 8893284334350876887982683022464217743514531016184218020225, 58839382658211258583908739382534119830998391168269411727984, 389393448030386221345676747652575518004370415855851052692084, 2577631282919200956633028190053768646831611621675032460656732, 17067177441141637299333044076622552421990396597411263436215128, 113033874001160162959279988594397749859516049115408248169631120, 748787983804424579516665098357127267955239740580277791815018327, 4961464757037224531534985442670344980877697358568843600799306847, 32882080174687947061547830350646698775706069633619231757847290833, 217973959654500204011489852387188606183318411529865773814157742652, 1445252109923910376344144995487752011623333719727290501567878101375, 9584601763460302229579715758517432380901823683912811373762124102225, 63576097515282316603895885789526036214977764196665380350439340193600, 421794522825528197054812704341545150612473166601324799290501001836160, 2798938115352065287582716121471351131365858630376626215322733246583440, 18576719705401457735393791293940952349191454540384945883946343595976400, 123317944475304041534070828910017600069789165812654835972878434854616080, 818772502536805216137581421420934726746919383028500564694002449108580480, 5437236316053020537834908352297364717424827453983329983815489236809477600, 36113491439880878000389421474642993568226346638040392335147107484021962400, 239902862756655477969490010268904992734553271992656034188565891324117528000, 1593949260304043599506518415982289783085823828321335142998353636056928131200, 10592174804426530661811370192734642970794251252936514832220762884943715581700, 70398903279856267626202294675106298761480269429707320229058060062421316645700, 467967101083449406896355810879265085737895415777067193322072779449434230312300, 3111229424897731247920647891425112478044238442641503526791791921845941601254000]