{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "### Example 2.9 (Rooted 3-Ary Trees)\n", "Some computational experiments with the generating function for rooted 3-ary trees. \n", "*Requirements: [ore_algebra package](https://github.com/mkauers/ore_algebra/)* \n", "**Note: See the Maple worksheet corresponding to this example for a better treatment of algebraic functions**" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "y^3*z - y + 1" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Define the minimal polynomial P(z,y) of the generating function y(z)\n", "var('z y')\n", "p = z*y^3-y+1\n", "p" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Symbolic Ring" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# There are three solutions -- the generating function is the one defined at z=0\n", "parent(p)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# Define the ring of linear differential operators with polynomial coefficients\n", "from ore_algebra import *\n", "Pols. = PolynomialRing(QQ); Diff. = OreAlgebra(Pols)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(27 z^{2} - 4 z\\right) \\mathit{Dz}^{2} + \\left(54 z - 6\\right) \\mathit{Dz} + 6\n", "\\end{math}" ], "text/plain": [ "(27*z^2 - 4*z)*Dz^2 + (54*z - 6)*Dz + 6" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# The GF here satisfies (27*z^2 - 4*z)*y''(z) + (54*z - 6)*y'(z) + 6*y(z) = 0\n", "deq = (27*z^2 - 4*z)*Dz^2 + (54*z - 6)*Dz + 6\n", "show(deq)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[1 + z + 3*z^2 + 12*z^3 + 55*z^4 + O(z^5),\n", " z^(-1/2)*(1 - 3/8*z - 105/128*z^2 - 3003/1024*z^3 - 415701/32768*z^4 + O(z^5))]" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# A basis of series spanning the solutions to the differential equation at the origin can be computed\n", "# Here there are two basis elements (note the original equation has three solutions!)\n", "deq.generalized_series_solutions(5)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The generating function expansion begins 1 + Z + 3*Z^2 + 12*Z^3 + O(Z^4)\n", "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)\n", "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)\n" ] } ], "source": [ "# The previous two functions span the C-vector space of solutions\n", "# The generating function is the first basis element\n", "# The other two solutions are linear combinations of the basis elements\n", "# (This can be determined after using the Newton polynom method to compute some initial terms of the solutions)\n", "\n", "# First, convert from ore_algebra data structure to PuiseuxSeriesRing\n", "order = 20\n", "[bas1, bas2] = deq.generalized_series_solutions(order)\n", "P. = PuiseuxSeriesRing(QQ, 'Z')\n", "bas1 = add(bas1[k]*Z^k for k in range(order))\n", "bas2 = Z^(-1/2)*add(bas2[k]*Z^k for k in range(order))\n", "\n", "# Define (truncations of) the solutions of the original algebraic equation\n", "sol1, sol2, sol3 = bas1, -bas1/2+bas2, -bas1/2-bas2\n", "print(\"The generating function expansion begins {}\".format(sol1 + O(Z^4)))\n", "print(\"The first extraneous root expansion begins {}\".format(sol2 + O(Z^4)))\n", "print(\"The second extraneous root expansion begins {}\".format(sol3 + O(Z^4)))" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "O(Z^20)\n", "O(Z^(39/2))\n", "O(Z^(39/2))\n" ] } ], "source": [ "# Verify these truncations satisfy the original polynomial equation\n", "print(QQ['y,z'](p).subs(z=Z, y=sol1) + O(Z^order))\n", "print(QQ['y,z'](p).subs(z=Z, y=sol2) + O(Z^(order-1/2)))\n", "print(QQ['y,z'](p).subs(z=Z, y=sol3) + O(Z^(order-1/2)))" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(-4*n^2 - 10*n - 6)*Sn + 27*n^2 + 27*n + 6" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# The coefficients of an algebraic function satisfy a P-recurrence\n", "# This computation shows (4*n^2 + 10*n + 6)*c_{n+1} = (27*n^2 + 27*n + 6)*c_n\n", "Ind. = PolynomialRing(QQ); Shift. = OreAlgebra(Ind)\n", "rec = deq.to_S(Shift)\n", "rec" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[1,\n", " 1,\n", " 3,\n", " 12,\n", " 55,\n", " 273,\n", " 1428,\n", " 7752,\n", " 43263,\n", " 246675,\n", " 1430715,\n", " 8414640,\n", " 50067108,\n", " 300830572,\n", " 1822766520,\n", " 11124755664,\n", " 68328754959,\n", " 422030545335,\n", " 2619631042665,\n", " 16332922290300,\n", " 102240109897695,\n", " 642312451217745,\n", " 4048514844039120,\n", " 25594403741131680,\n", " 162250238001816900,\n", " 1031147983159782228,\n", " 6568517413771094628,\n", " 41932353590942745504,\n", " 268225186597703313816,\n", " 1718929965542850284040,\n", " 11034966795189838872624,\n", " 70956023048640039202464,\n", " 456949965738717944767791,\n", " 2946924270225408943665279,\n", " 19030649059639789214206725,\n", " 123052100237542105872786180,\n", " 796607831560617902288322405,\n", " 5162879946168545215371343587,\n", " 33496962712940417760973884708,\n", " 217550867863011281855594752680,\n", " 1414282077098335379544565517191,\n", " 9202600068524372703278082352971,\n", " 59932899605936040714626166584475,\n", " 390645234961546222075767026462400,\n", " 2548271840422037344975860237738000,\n", " 16635641296703864308483436321233200,\n", " 108679967966404768511803431114031200,\n", " 710496811856518520237299474629019200,\n", " 4647985909007237458743106679064711300,\n", " 30426054945480277365983787382745806500,\n", " 199293672373583488061784498821082334140,\n", " 1306164582615216509992597891759417970640,\n", " 8565425435995939036682228482499260153620,\n", " 56200126088515763642375073109736372959980,\n", " 368938646289542831658405169038432458575200,\n", " 2423210784425567315944099994775276777059520,\n", " 15923573648667567272266673527613832841413720,\n", " 104687476184489521220187993237607647032909880,\n", " 688567758918141068396974784875829607672693720,\n", " 4530954345822264258644101276351421147198184000,\n", " 29827534579873756985313709683286979461162921200,\n", " 196437153964454478730548625427197056659262085200,\n", " 1294204010506457507610556492169378285551215751008,\n", " 8530003035876456718702048001257342146013974679872,\n", " 56241607952753745334610705298987817310837053258095,\n", " 370956832348855971723751149936236380931126414813199,\n", " 2447605128273073612473937013215645786143652277643069,\n", " 16155005655931729647837925891539419948344504253763540,\n", " 106663615595824898981200406146704399055680776314377585,\n", " 704473645075333982479839404375246601833187235495177575,\n", " 4654236823512905605830331092249349585849895972864188380,\n", " 30758394797997531797498626576159446322402868116541053080,\n", " 203332434576992304310763190915100018001976431327881300605,\n", " 1344540697525747546508940943448381528040466460776733688345,\n", " 8893284334350876887982683022464217743514531016184218020225,\n", " 58839382658211258583908739382534119830998391168269411727984,\n", " 389393448030386221345676747652575518004370415855851052692084,\n", " 2577631282919200956633028190053768646831611621675032460656732,\n", " 17067177441141637299333044076622552421990396597411263436215128,\n", " 113033874001160162959279988594397749859516049115408248169631120,\n", " 748787983804424579516665098357127267955239740580277791815018327,\n", " 4961464757037224531534985442670344980877697358568843600799306847,\n", " 32882080174687947061547830350646698775706069633619231757847290833,\n", " 217973959654500204011489852387188606183318411529865773814157742652,\n", " 1445252109923910376344144995487752011623333719727290501567878101375,\n", " 9584601763460302229579715758517432380901823683912811373762124102225,\n", " 63576097515282316603895885789526036214977764196665380350439340193600,\n", " 421794522825528197054812704341545150612473166601324799290501001836160,\n", " 2798938115352065287582716121471351131365858630376626215322733246583440,\n", " 18576719705401457735393791293940952349191454540384945883946343595976400,\n", " 123317944475304041534070828910017600069789165812654835972878434854616080,\n", " 818772502536805216137581421420934726746919383028500564694002449108580480,\n", " 5437236316053020537834908352297364717424827453983329983815489236809477600,\n", " 36113491439880878000389421474642993568226346638040392335147107484021962400,\n", " 239902862756655477969490010268904992734553271992656034188565891324117528000,\n", " 1593949260304043599506518415982289783085823828321335142998353636056928131200,\n", " 10592174804426530661811370192734642970794251252936514832220762884943715581700,\n", " 70398903279856267626202294675106298761480269429707320229058060062421316645700,\n", " 467967101083449406896355810879265085737895415777067193322072779449434230312300,\n", " 3111229424897731247920647891425112478044238442641503526791791921845941601254000]" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# The recurrence can be used to compute many terms of the sequence efficiently\n", "# We must give the initial term in the sequence (to specify which solution we care about)\n", "rec.to_list([1],100)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "SageMath 9.2", "language": "sage", "name": "sagemath" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.5" } }, "nbformat": 4, "nbformat_minor": 4 }