{ "cells": [ { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Recent improvements to ore_algebra\n", "\n", "Marc Mezzarobba \n", "CNRS, Sorbonne Université\n", "\n", "FastRelax Final Workshop, Lyon, 2019-05-23" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### ore_algebra\n", "\n", "Sage implementation of Ore polynomials and D-finite functions\n", "\n", "**Source code:** \n", "\n", " $ sage -pip install git+https://github.com/mkauers/ore_algebra/\n", " \n", "**Documentation:** \n", "\n", "### History & Contributors\n", "* **Kauers, Jaroschek, Johansson 2013:** Initial implementation\n", "* **Mezzarobba ~2015–:** Numerical solutions of ODEs\n", "* **Kauers 2015–, Mezzarobba 2018–:** Multivariate operators\n", "* **Schwaiger 2017, Hofstadler 2017–2018:** “D-finite function” objects" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Pre-existing features\n", "\n", "* Basic arithmetic (diff, shift, qdiff, qshift, custom)\n", "* Gcrd, lclm, D-finite closure properties\n", "* Conversions (diff. eq. ↔ rec., etc.)\n", "* Polynomial & rational solutions\n", "* First-order right factors\n", "* Formal solutions at singularities\n", "* Desingularization\n", "* Guessing\n", "* Numerical solutions & connection matrices\n", "\n", "### Recent improvements\n", "* Multivariate closure properties\n", "* Creative telescoping\n", "* D-finite functions (alpha quality)\n", "* Numerics: faster, better support for reg. sing. points, large operators, high precision" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Differential operators: Basic arithmetic" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "scrolled": true, "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "from ore_algebra import OreAlgebra" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "scrolled": true }, "outputs": [], "source": [ "Pol. = PolynomialRing(QQ)\n", "Dop. = OreAlgebra(Pol)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Univariate Ore algebra in Dx over Univariate Polynomial Ring in x over Rational Field" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Dop" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "x*Dx + 1" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Dx*x" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "cos(x)" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Dx(sin(x))" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "x*cos(x)" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "(x*Dx)(sin(x))" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Recurrence operators" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/plain": [ "Univariate Ore algebra in Sn over Univariate Polynomial Ring in n over Finite Field of size 17" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Poln. = PolynomialRing(GF(17))\n", "Rop. = OreAlgebra(Poln)\n", "Rop" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(n + 1)*Sn" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Sn*n" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Sn^10 + (7*n + 6)*Sn^9 + (11*n^2 + 3*n + 2)*Sn^8 + (16*n^3 + 15*n^2 + 4*n + 2)*Sn^7 + (6*n^4 + 4*n^3 + 2*n^2 + 8*n + 13)*Sn^6 + (3*n^5 + 12*n^4 + 13*n^3 + 10*n^2 + 3*n + 9)*Sn^5 + (6*n^6 + 4*n^5 + 16*n^4 + 10*n^3 + 11*n^2 + 10*n + 3)*Sn^4 + (16*n^7 + 15*n^6 + 7*n^5 + 4*n^4 + 5*n^3 + 16*n^2 + 10*n + 3)*Sn^3 + (11*n^8 + 3*n^7 + 8*n^6 + 6*n^5 + 16*n^4 + 12*n^3 + 3*n^2 + 1)*Sn^2 + (7*n^9 + 6*n^8 + 16*n^7 + 11*n^6 + 3*n^5 + 11*n^4 + 16*n^3 + 6*n^2 + 7*n + 16)*Sn + n^10" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "(Sn - n)**10" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### ⚐ D-Finite functions as objects \n", "C. Hofstadler, experimental" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "scrolled": true, "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "from ore_algebra.dfinite_function import DFiniteFunctionRing, UnivariateDFiniteFunction" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "scrolled": true, "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/plain": [ "Ring of D-finite functions over Univariate Polynomial Ring in x over Rational Field" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Fun = DFiniteFunctionRing(Dop); Fun" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "scrolled": true, "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "my_ai = Fun(airy_ai(x))" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "-1/136080*3^(2/3)/gamma(1/3)" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "my_ai[10]" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "scrolled": true, "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/plain": [ "[0.13529241631288141552414742351546630617494414298833 +/- 1.49e-51]" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "my_ai(1)" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "scrolled": true, "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/plain": [ "Univariate D-finite function defined by the annihilating operator x*Dx^2 - Dx - 4*x^5 and the coefficient sequence defined by (n^10 + 22*n^9 + 186*n^8 + 708*n^7 + 777*n^6 - 2562*n^5 - 7876*n^4 - 3928*n^3 + 6912*n^2 + 5760*n)*Sn^6 - 4*n^8 - 48*n^7 - 168*n^6 + 924*n^4 + 1008*n^3 - 752*n^2 - 960*n and {0: 1/3*3^(1/3)/gamma(2/3), 1: 0, 2: -1/3*3^(2/3)/gamma(1/3), 3: 0, 4: 0, 5: 0, 6: 1/18*3^(1/3)/gamma(2/3), 7: 0, 8: -1/36*3^(2/3)/gamma(1/3)}" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "my_ai(x^2)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "scrolled": true, "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "data": { "text/plain": [ "IVP(pt=0, dop=1/3*Dx^3 - x^2*Dx^2 - 4*x*Dx - 2, ini=[1, 1/gamma(4/3), 1/gamma(5/3)])" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from ore_algebra.examples.stdfun import mittag_leffler_e\n", "ivp = mittag_leffler_e(1/3, 1)\n", "ivp" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "scrolled": true, "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "my_e = UnivariateDFiniteFunction(Fun, ivp.dop, ivp.ini)" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/plain": [ "Univariate D-finite function defined by the annihilating operator (-81*x^9 - 414*x^6 - 52*x^3 - 24)*Dx^6 + (486*x^11 + 3213*x^8 + 2796*x^5 + 300*x^2)*Dx^5 + (-729*x^13 - 3483*x^10 + 5148*x^7 - 3732*x^4 + 2064*x)*Dx^4 + (-6075*x^12 - 47574*x^9 - 33384*x^6 - 7464*x^3 + 1968)*Dx^3 + (1458*x^14 + 4293*x^11 - 53604*x^8 - 2244*x^5 - 41088*x^2)*Dx^2 + (7047*x^13 + 60075*x^10 + 71808*x^7 + 136452*x^4 - 31872*x)*Dx - 729*x^15 - 2430*x^12 + 31752*x^9 - 22028*x^6 + 46720*x^3 - 2160 and the coefficient sequence defined by (-24*n^33 - 7848*n^32 - 1217904*n^31 - 119265120*n^30 - 8265240096*n^29 - 430836655392*n^28 - 17522874893376*n^27 - 569149922442240*n^26 - 14980851481012560*n^25 - 322206414583126320*n^24 - 5678126299596310560*n^23 - 81721683804628296000*n^22 - 949652588246790733920*n^21 - 8681713167577264503840*n^20 - 58748680805012448505920*n^19 - 242135956214206048051200*n^18 + 108488260436315928944040*n^17 + 10476807873210919884061080*n^16 + 84066023678942093501383440*n^15 + 328407072430892969322928800*n^14 + 161957855017977822122938944*n^13 - 5556526268492872609705821312*n^12 - 30808184921580886798491192576*n^11 - 61848743568136626694433518080*n^10 + 96516578682202812595845129216*n^9 + 853030845899725268959565893632*n^8 + 1755693662806275221362303696896*n^7 - 191121075488552116917636956160*n^6 - 6918611018576047183218416025600*n^5 - 10305536133084558124200099840000*n^4 - 651883440675078087990312960000*n^3 + 9710692999529424209156505600000*n^2 + 5748846330021297285758976000000*n)*Sn^21 + (-52*n^33 - 15768*n^32 - 2266192*n^31 - 205228464*n^30 - 13131424864*n^29 - 630793023888*n^28 - 23590821575184*n^27 - 702719853964272*n^26 - 16907977517369640*n^25 - 331025204263890240*n^24 - 5279792730717476400*n^23 - 68200672898298404880*n^22 - 701622354982097190720*n^21 - 5530245187115752498800*n^20 - 30124828013909392092720*n^19 - 68576800943396360742480*n^18 + 571745721721581523264380*n^17 + 7492473681331898408673240*n^16 + 41649036729364135964121280*n^15 + 99768141174052190820326400*n^14 - 273482617801028411195658368*n^13 - 3270682878463161988366755072*n^12 - 11711448970574997636645705728*n^11 - 10059256312341487086991773696*n^10 + 76801979904398494328865844224*n^9 + 322846351399709844933782550528*n^8 + 422861006998310521457169874944*n^7 - 488850815277845481558737092608*n^6 - 2352232797247927527470849064960*n^5 - 2540273963213858698740695040000*n^4 + 482772293714613829168005120000*n^3 + 2719501179843522132317306880000*n^2 + 1381740258268276716050841600000*n)*Sn^18 + (-414*n^33 - 117678*n^32 - 15852576*n^31 - 1345594824*n^30 - 80700092856*n^29 - 3633815655240*n^28 - 127399183240896*n^27 - 3557731583055600*n^26 - 80245669721651460*n^25 - 1472279089838630580*n^24 - 21987171811970146080*n^23 - 265389681264853335840*n^22 - 2539134234244059547320*n^21 - 18385725342062668599720*n^20 - 88138188290507684371200*n^19 - 108534146024035898174160*n^18 + 2443632611990285662853490*n^17 + 24833637843205336319220210*n^16 + 121245809437122843644887680*n^15 + 227257412847909330770915880*n^14 - 1106099032265374878229720416*n^13 - 9933584593498812982942176672*n^12 - 31972810234198476727554597504*n^11 - 18761909942962878443753678976*n^10 + 235289950832015737087257086976*n^9 + 902805274335383262934308679680*n^8 + 1108702158744622782158950680576*n^7 - 1476049778974732201091769876480*n^6 - 6591499425554580260012924928000*n^5 - 7038297709246683060184350720000*n^4 + 1380775208407863970913648640000*n^3 + 7639985744558539462921420800000*n^2 + 3899687418093882578042880000000*n)*Sn^15 + (-81*n^33 - 18900*n^32 - 1999656*n^31 - 124129380*n^30 - 4754325036*n^29 - 94831186068*n^28 + 798409333576*n^27 + 136202820180876*n^26 + 5626859890107690*n^25 + 150244980230360340*n^24 + 2958076803935621160*n^23 + 44707003575753198420*n^22 + 522317644089199847700*n^21 + 4622690166053310278820*n^20 + 28843536943661955489720*n^19 + 94346507950084965030660*n^18 - 297866764199743660754145*n^17 - 6149478723439194914314320*n^16 - 39608446972567185637607040*n^15 - 117688111266949538549650560*n^14 + 137113268102502964083036576*n^13 + 2913079673575808379081154560*n^12 + 12078165530462456499372270336*n^11 + 15374326616664284229238287360*n^10 - 63952753134498320456998269184*n^9 - 329179698517056051282299154432*n^8 - 512742818335596794644189188096*n^7 + 357125936983452978364669722624*n^6 + 2474569756741412395763355156480*n^5 + 3018711856234964620121210880000*n^4 - 254978118045128898082897920000*n^3 - 3064821762415666751880560640000*n^2 - 1655071469079834867322060800000*n)*Sn^12 + (486*n^32 + 118503*n^31 + 13524948*n^30 + 958340016*n^29 + 47116511700*n^28 + 1699711744468*n^27 + 46308339416268*n^26 + 962659458170280*n^25 + 15123563382620160*n^24 + 171523512133221690*n^23 + 1186100656517839140*n^22 - 37805196339922440*n^21 - 119715775715274835860*n^20 - 1652851915028270167500*n^19 - 12217664497859617376700*n^18 - 42469132523979420144360*n^17 + 126657823154351913282330*n^16 + 2472822893446539301290135*n^15 + 13578493106461078407258360*n^14 + 24139238444251948125814680*n^13 - 138749005788261364047779136*n^12 - 1053511283509579824159117888*n^11 - 2542721346034360611629054208*n^10 + 2052017446773332083506041984*n^9 + 28003740422874185244995320320*n^8 + 61582361453291997219785870592*n^7 - 480533425570793441674487808*n^6 - 226204572707464956969431900160*n^5 - 341180558039329149831045120000*n^4 - 22763074540137562621992960000*n^3 + 316325128579097767572602880000*n^2 + 186360211691558827956633600000*n)*Sn^9 + (-729*n^31 - 176661*n^30 - 20099340*n^29 - 1425944763*n^28 - 70643005308*n^27 - 2592496126863*n^26 - 72915039882360*n^25 - 1602498019959135*n^24 - 27763983229969350*n^23 - 378691305544396845*n^22 - 4005949233962800800*n^21 - 31524964604749905705*n^20 - 163123050618207140940*n^19 - 252729347108177982645*n^18 + 4264442956387338816600*n^17 + 44947121910163022424555*n^16 + 212775455886957615862935*n^15 + 314133074661371762672550*n^14 - 2470847682946988384781300*n^13 - 17927143197524940083201640*n^12 - 46978390051656573808514736*n^11 + 13491956267542245956894496*n^10 + 462244206112205932904523840*n^9 + 1276871492416154316862946688*n^8 + 775173181989109336044208128*n^7 - 3421420839768170941727404032*n^6 - 8336006423977666738296176640*n^5 - 5406695301365082643599360000*n^4 + 4263545613333384992194560000*n^3 + 7555321040085118538219520000*n^2 + 2884275787535563515494400000*n)*Sn^6 + (1458*n^29 + 329265*n^28 + 34682661*n^27 + 2260935909*n^26 + 102022226865*n^25 + 3374719185825*n^24 + 84459402937845*n^23 + 1624738560888105*n^22 + 24090362505771705*n^21 + 271805463177566175*n^20 + 2238594275860798635*n^19 + 11830705538422983615*n^18 + 16440329540618410995*n^17 - 340948354047289111125*n^16 - 3395596953711720654585*n^15 - 14616943639216277213565*n^14 - 10632331877126123028735*n^13 + 230072441021389691864100*n^12 + 1274975782424092590278580*n^11 + 2318852894373189748894320*n^10 - 4894627815672088825635408*n^9 - 34165273162434107755314240*n^8 - 59650619585806830886883136*n^7 + 28654842091496143828271616*n^6 + 258699952696179776376453120*n^5 + 320170609559343027732480000*n^4 - 27066930120236510300160000*n^3 - 317194158035945496821760000*n^2 - 168348741731156341555200000*n)*Sn^3 - 729*n^27 - 157464*n^26 - 15779205*n^25 - 972340200*n^24 - 41152593105*n^23 - 1264380020280*n^22 - 29021590766025*n^21 - 503177814196200*n^20 - 6550692787196415*n^19 - 62011070347323240*n^18 - 386150232022854975*n^17 - 950333662778830200*n^16 + 8397278950113060165*n^15 + 105798884663200766040*n^14 + 499294024107256828125*n^13 + 406201588063165369800*n^12 - 8081520464744734885020*n^11 - 43218065849227900693920*n^10 - 60308274941946097981200*n^9 + 247545163236333885436800*n^8 + 1148925647418746051335104*n^7 + 1122767140648429530468864*n^6 - 2853333449859192239646720*n^5 - 7415089937678104381440000*n^4 - 2697145321907074590720000*n^3 + 6087484712019018792960000*n^2 + 4469435621181141811200000*n and {0: 1/3*3^(1/3)/gamma(2/3), 1: -1/3*3^(2/3)/gamma(1/3) + 1/3*3^(1/3)/(gamma(4/3)*gamma(2/3)), 2: -1/3*3^(2/3)/(gamma(4/3)*gamma(1/3)) + 1/3*3^(1/3)/(gamma(5/3)*gamma(2/3)), 3: 7/18*3^(1/3)/gamma(2/3) - 1/3*3^(2/3)/(gamma(5/3)*gamma(1/3)), 4: -13/36*3^(2/3)/gamma(1/3) + 11/36*3^(1/3)/(gamma(4/3)*gamma(2/3)), 5: -5/18*3^(2/3)/(gamma(4/3)*gamma(1/3)) + 23/90*3^(1/3)/(gamma(5/3)*gamma(2/3)), 6: 121/540*3^(1/3)/gamma(2/3) - 41/180*3^(2/3)/(gamma(5/3)*gamma(1/3)), 7: -295/1512*3^(2/3)/gamma(1/3) + 1139/7560*3^(1/3)/(gamma(4/3)*gamma(2/3)), 8: -389/3024*3^(2/3)/(gamma(4/3)*gamma(1/3)) + 119/1080*3^(1/3)/(gamma(5/3)*gamma(2/3)), 9: 3313/38880*3^(1/3)/gamma(2/3) - 349/3780*3^(2/3)/(gamma(5/3)*gamma(1/3)), 10: -1363/19440*3^(2/3)/gamma(1/3) + 1999/38880*3^(1/3)/(gamma(4/3)*gamma(2/3)), 11: -2263/54432*3^(2/3)/(gamma(4/3)*gamma(1/3)) + 72901/2138400*3^(1/3)/(gamma(5/3)*gamma(2/3)), 12: 24737/1026432*3^(1/3)/gamma(2/3) - 81157/2993760*3^(2/3)/(gamma(5/3)*gamma(1/3)), 13: -400297/21228480*3^(2/3)/gamma(1/3) + 1563307/116756640*3^(1/3)/(gamma(4/3)*gamma(2/3)), 14: -766349/74299680*3^(2/3)/(gamma(4/3)*gamma(1/3)) + 1477331/179625600*3^(1/3)/(gamma(5/3)*gamma(2/3)), 15: 5835271/1077753600*3^(1/3)/gamma(2/3) - 3643349/583783200*3^(2/3)/(gamma(5/3)*gamma(1/3)), 16: -20629681/5094835200*3^(2/3)/gamma(1/3) + 78854057/28021593600*3^(1/3)/(gamma(4/3)*gamma(2/3)), 17: -10572301/5094835200*3^(2/3)/(gamma(4/3)*gamma(1/3)) + 7436603/4580452800*3^(1/3)/(gamma(5/3)*gamma(2/3)), 18: 332256637/329792601600*3^(1/3)/gamma(2/3) - 1125195523/952734182400*3^(2/3)/(gamma(5/3)*gamma(1/3)), 19: -1260215623/1742433638400*3^(2/3)/gamma(1/3) + 80668819493/162917545190400*3^(1/3)/(gamma(4/3)*gamma(2/3)), 20: -8560868333/24394070937600*3^(2/3)/(gamma(4/3)*gamma(1/3)) + 127635493/471132288000*3^(1/3)/(gamma(5/3)*gamma(2/3)), 21: 22139176921/138512892672000*3^(1/3)/gamma(2/3) - 2162693257679/11404228163328000*3^(2/3)/(gamma(5/3)*gamma(1/3)), 22: -17797586453/161000868188160*3^(2/3)/gamma(1/3) + 29105422369/388780505568000*3^(1/3)/(gamma(4/3)*gamma(2/3)), 23: -3295558349/64400347275264*3^(2/3)/(gamma(4/3)*gamma(1/3)) + 17135676017/438047023075200*3^(1/3)/(gamma(5/3)*gamma(2/3)), 24: 241225777519/10922730964992000*3^(1/3)/gamma(2/3) - 83210444770261/3147566973078528000*3^(2/3)/(gamma(5/3)*gamma(1/3)), 25: -7129346024401/483002604564480000*3^(2/3)/gamma(1/3) + 187269180969109/18885401838471168000*3^(1/3)/(gamma(4/3)*gamma(2/3)), 26: -41167459351903/6279033859338240000*3^(2/3)/(gamma(4/3)*gamma(1/3)) + 271991502318703/54668268479784960000*3^(1/3)/(gamma(5/3)*gamma(2/3))}" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "my_ai*my_e" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "scrolled": true, "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "data": { "text/plain": [ "[7.28355749378782728852876426824018405440330223126660 +/- 8.00e-51]" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "my_e(1)" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[7.28355749378782728852876426824018405440330223126660 +/- 8.00e-51]" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ivp.dop.numerical_solution(ivp.ini, [0,1], 1e-50)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Kreweras walks\n", "Adapted from an exercise by B. Salvy & A. Bostan" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "scrolled": true, "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "@cached_function\n", "def w(n, i, j):\n", " if i < 0 or j < 0:\n", " return 0\n", " elif n == 0:\n", " return 1 if i == j == 0 else 0\n", " else:\n", " return w(n-1,i-1,j-1)+w(n-1,i,j+1)+w(n-1,i+1,j)" ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "[1,\n", " 0,\n", " x,\n", " 2,\n", " 2*x^2,\n", " 8*x,\n", " 5*x^3 + 16,\n", " 30*x^2,\n", " 14*x^4 + 96*x,\n", " 112*x^3 + 192,\n", " 42*x^5 + 480*x^2,\n", " 420*x^4 + 1408*x,\n", " 132*x^6 + 2240*x^3 + 2816,\n", " 1584*x^5 + 8320*x^2,\n", " 429*x^7 + 10080*x^4 + 23296*x,\n", " 6006*x^6 + 44800*x^3 + 46592,\n", " 1430*x^8 + 44352*x^5 + 153600*x^2,\n", " 22880*x^7 + 228480*x^4 + 417792*x,\n", " 4862*x^9 + 192192*x^6 + 913920*x^3 + 835584,\n", " 87516*x^8 + 1123584*x^5 + 2976768*x^2,\n", " 16796*x^10 + 823680*x^7 + 5107200*x^4 + 7938048*x,\n", " 335920*x^9 + 5381376*x^6 + 19066880*x^3 + 15876096,\n", " 58786*x^11 + 3500640*x^8 + 27320832*x^5 + 59924480*x^2,\n", " 1293292*x^10 + 25259520*x^7 + 114250752*x^4 + 157515776*x,\n", " 208012*x^12 + 14780480*x^9 + 141453312*x^6 + 406224896*x^3 + 315031552,\n", " 4992288*x^11 + 116688000*x^8 + 652861440*x^5 + 1243545600*x^2,\n", " 742900*x^13 + 62078016*x^10 + 713856000*x^7 + 2571878400*x^4 + 3233218560*x,\n", " 19315400*x^12 + 532097280*x^9 + 3597834240*x^6 + 8817868800*x^3 + 6466437120,\n", " 2674440*x^14 + 259598976*x^11 + 3528645120*x^8 + 15498362880*x^5 + 26453606400*x^2,\n", " 74884320*x^13 + 2400349952*x^10 + 19262251008*x^7 + 58370457600*x^4 + 68191518720*x,\n", " 9694845*x^15 + 1081662400*x^12 + 17145356800*x^9 + 89890504704*x^6 + 194568192000*x^3 + 136383037440,\n", " 290845350*x^14 + 10730091008*x^11 + 100706411520*x^8 + 367464480768*x^5 + 574439424000*x^2,\n", " 35357670*x^16 + 4493059200*x^13 + 82108522496*x^10 + 505506693120*x^7 + 1336234475520*x^4 + 1470564925440*x,\n", " 1131445440*x^15 + 47593145600*x^12 + 516110712832*x^9 + 2224229449728*x^6 + 4355134586880*x^3 + 2941129850880,\n", " 129644790*x^17 + 18614102400*x^14 + 388360068096*x^11 + 2770699616256*x^8 + 8725823225856*x^5 + 12692106510336*x^2,\n", " 4407922860*x^16 + 209676096000*x^13 + 2600323934208*x^10 + 13038586429440*x^7 + 30848869990400*x^4 + 32307180208128*x,\n", " 477638700*x^18 + 76938289920*x^15 + 1817192832000*x^12 + 14858993909760*x^9 + 54762063003648*x^6 + 98716383969280*x^3 + 64614360416256,\n", " 17194993200*x^17 + 918295718400*x^14 + 12909337878528*x^11 + 74406691307520*x^8 + 207815008321536*x^5 + 284610873262080*x^2,\n", " 1767263190*x^19 + 317370445920*x^16 + 8422988313600*x^13 + 78204394684416*x^10 + 332641678786560*x^7 + 717906392383488*x^4 + 721014212263936*x,\n", " 67156001220*x^18 + 4000791075840*x^15 + 63266001111040*x^12 + 414962094243840*x^9 + 1345350789758976*x^6 + 2262492872966144*x^3 + 1442028424527872,\n", " 6564120420*x^20 + 1306819483200*x^17 + 38717332992000*x^14 + 404902407110656*x^11 + 1965609920102400*x^8 + 4967449069879296*x^5 + 6464265351331840*x^2,\n", " 262564816800*x^19 + 17349584376960*x^16 + 306534470860800*x^13 + 2268459448532992*x^10 + 8427713121484800*x^7 + 16831852220252160*x^4 + 16309838732591104*x,\n", " 24466267020*x^21 + 5372480097600*x^18 + 176650313656320*x^15 + 2066269396172800*x^12 + 11342297242664960*x^9 + 33036635436220416*x^6 + 52365762463006720*x^3 + 32619677465182208,\n", " 1027583214840*x^20 + 74924317036800*x^17 + 1470186481397760*x^14 + 12185085810573312*x^11 + 51338819098378240*x^8 + 119209118077550592*x^5 + 148465568301711360*x^2,\n", " 91482563640*x^22 + 22055444611200*x^19 + 800619844907520*x^16 + 10409596236103680*x^13 + 64104147090407424*x^10 + 212603109677989888*x^7 + 397363726925168640*x^4 + 373284857444302848*x,\n", " 4025232800160*x^21 + 322348805856000*x^18 + 6987227737374720*x^15 + 64440357652070400*x^12 + 305257843287654400*x^9 + 811757327861415936*x^6 + 1222657621308211200*x^3 + 746569714888605696,\n", " 343059613650*x^23 + 90427322905920*x^20 + 3606821773632000*x^17 + 51840721922457600*x^14 + 355710774239428608*x^11 + 1330281548643041280*x^8 + 2872372083201933312*x^5 + 3443403096745574400*x^2,\n", " 15780742227900*x^22 + 1382141195635200*x^19 + 32935435395793920*x^16 + 336070886945587200*x^13 + 1776835461659754496*x^10 + 5349581414863994880*x^7 + 9440663490244116480*x^4 + 8631463762508906496*x,\n", " 1289904147324*x^24 + 370321417614720*x^21 + 16160420133580800*x^18 + 255499135191613440*x^15 + 1941742902352281600*x^12 + 8122676396158877696*x^9 + 19971770615492247552*x^6 + 28771545875029688320*x^3 + 17262927525017812992,\n", " 61915399071552*x^23 + 5907918429853440*x^20 + 154091789814251520*x^17 + 1730800593233510400*x^14 + 10148842902961258496*x^11 + 34278471681397751808*x^8 + 69487580603014447104*x^5 + 80560328450083127296*x^2,\n", " 4861946401452*x^25 + 1514951253878400*x^22 + 72047785729920000*x^19 + 1247409727067750400*x^16 + 10444486338478080000*x^13 + 48537944318510366720*x^10 + 134425379142736281600*x^7 + 225609027931865088000*x^4 + 201400821125207818240*x,\n", " 243097320072600*x^24 + 25181856397800960*x^21 + 716044147408128000*x^18 + 8812868721102028800*x^15 + 56993369896386560000*x^12 + 214323390497318502400*x^9 + 492203695938019000320*x^6 + 681840617749636710400*x^3 + 402801642250415636480,\n", " 18367353072152*x^26 + 6191539907155200*x^23 + 319751013795333120*x^20 + 6037993891657728000*x^17 + 55435787116609536000*x^14 + 284510902522761707520*x^11 + 879853918883728588800*x^8 + 1687555528930350858240*x^5 + 1899413149445416550400*x^2,\n", " 955102359751904*x^25 + 107056555274073600*x^22 + 3306693410956615680*x^19 + 44410061194968268800*x^16 + 315198115253059584000*x^13 + 1311224159452727869440*x^10 + 3376091055264261734400*x^7 + 5420632911109611847680*x^4 + 4737359855087391866880*x,\n", " 69533550916004*x^27 + 25282121287550400*x^24 + 1413146529617771520*x^21 + 28997157603773399040*x^18 + 290684036912519577600*x^15 + 1639030199315909836800*x^12 + 5619532111940262297600*x^9 + 12153927798951342243840*x^6 + 16261898733328835543040*x^3 + 9474719710174783733760,\n", " 3754811749464216*x^26 + 454046259858048000*x^23 + 15183109225195591680*x^20 + 221677112183286988800*x^17 + 1719099143031029760000*x^14 + 7867344956716367216640*x^11 + 22523630732068257792000*x^8 + 41136371011835312209920*x^5 + 45096021697466518732800*x^2,\n", " 263747951750360*x^28 + 103151054853205632*x^25 + 6221399815927296000*x^22 + 138252701887959859200*x^19 + 1507404362846351523840*x^16 + 9294991228526395392000*x^13 + 35118003574907842068480*x^10 + 84794845108962852864000*x^7 + 130888453219475993395200*x^4 + 112238987335916668846080*x,\n", " 14769885298020160*x^27 + 1921441217853830400*x^24 + 69347869948202926080*x^21 + 1096906052341835366400*x^18 + 9257594470813956833280*x^15 + 46381067342343831552000*x^12 + 146646608334779999846400*x^9 + 300739050653121584824320*x^6 + 390098919399222568550400*x^3 + 224477974671833337692160,\n", " 1002242216651368*x^29 + 420538915939992192*x^26 + 27292307910740121600*x^23 + 654772911603962511360*x^20 + 7737634585438352179200*x^17 + 51961981868439628677120*x^14 + 215208152468475378401280*x^11 + 575560372562068871577600*x^8 + 1006319131031599149219840*x^5 + 1077416063102614713139200*x^2,\n", " 58130048565779344*x^28 + 8114549648452176384*x^25 + 315197121998930595840*x^22 + 5384195370680667340800*x^19 + 49275222153616934830080*x^16 + 269094652183894252584960*x^13 + 934248100013648967106560*x^10 + 2130701928229776528506880*x^7 + 3175017579190606941388800*x^4 + 2676528325181232340008960*x,\n", " 3814986502092304*x^30 + 1713306694570338560*x^27 + 119331612477237888000*x^24 + 3081927415100654714880*x^21 + 39346043093435645952000*x^18 + 286692201621043984465920*x^15 + 1295640917922453808742400*x^12 + 3813257551076118233088000*x^9 + 7457456748804217849774080*x^6 + 9407459493898094641152000*x^3 + 5353056650362464680017920,\n", " 228899190125538240*x^29 + 34203831829786031616*x^26 + 1426134536054500147200*x^23 + 26232219393647433154560*x^20 + 259471203102656692224000*x^17 + 1538553457379796162969600*x^14 + 5836364011810684233842688*x^11 + 14691076459935360771686400*x^8 + 24700716471152431973007360*x^5 + 25888196802757162696704000*x^2,\n", " 14544636039226909*x^31 + 6975605827893521280*x^28 + 520156385939764933632*x^25 + 14423176229884519219200*x^22 + 198341171025139128729600*x^19 + 1562757988972572306309120*x^16 + 7675082764400362468147200*x^13 + 24722982956862649984352256*x^10 + 53579220030352492226150400*x^7 + 77345677838962160723558400*x^4 + 64202728070837763487825920*x,\n", " 901767434432068358*x^30 + 143917762343908439040*x^27 + 6425461238079449180160*x^24 + 126923950822983769128960*x^21 + 1352788499812487390822400*x^18 + 8679118748674120577187840*x^15 + 35817052900535024851353600*x^12 + 98891931827450599937409024*x^9 + 185321066928513326052802560*x^6 + 227966208367467421079961600*x^3 + 128405456141675526975651840,\n", " 55534064877048198*x^32 + 28383499575566741760*x^29 + 2260817212093398024192*x^26 + 67139513344830162862080*x^23 + 991777848291222009937920*x^20 + 8423850550183705266094080*x^17 + 44795451606059977172582400*x^14 + 157185695014919423347654656*x^11 + 374748373240865431341760512*x^8 + 608233245303838608583557120*x^5 + 625278742950767783533608960*x^2,\n", " 3554180152131084672*x^31 + 604552505084105177600*x^28 + 28836083686323340836864*x^25 + 610173692708334914764800*x^22 + 6988136600425954541568000*x^19 + 48355090171184386072903680*x^16 + 216253904305117131177984000*x^13 + 651523315569086305469988864*x^10 + 1348575460797578022821560320*x^7 + 1891634494964091366408192000*x^4 + 1548309268259044035416555520*x,\n", " 212336130412243110*x^33 + 115426231607304749824*x^30 + 9800114292942336563200*x^27 + 310977373087800734515200*x^24 + 4922067787847234979102720*x^21 + 44939093830431523051929600*x^18 + 257893814246316725722152960*x^15 + 981724073512119039950848000*x^12 + 2559555882592839057203527680*x^9 + 4615124910285044789211561984*x^6 + 5548794518561334674797363200*x^3 + 3096618536518088070833111040,\n", " 14014184607208045260*x^32 + 2535592628750628930560*x^29 + 128933867315874086129664*x^26 + 2915752862537513592422400*x^23 + 35789919263571367522467840*x^20 + 266322099063884112042393600*x^17 + 1286269396116865058812723200*x^14 + 4209632827219966443309236224*x^11 + 9556731871228990411725864960*x^8 + 15022511853721360366502412288*x^5 + 15174254397698343804547891200*x^2,\n", " 812944042149730764*x^34 + 469151780081303176704*x^31 + 42374480202510510602240*x^28 + 1433679732166448957718528*x^25 + 24256582324514208822067200*x^22 + 237435561943692974783201280*x^19 + 1466039745323095397719080960*x^16 + 6032159926617022344776908800*x^13 + 17113072580220298367365808128*x^10 + 33979491097703077019469742080*x^7 + 46433218456956932041916547072*x^4 + 37521792692490450134882058240*x,\n", " 55280194866181691952*x^33 + 10619213307872036983808*x^30 + 574508742114037238480896*x^27 + 13854888167995094969548800*x^24 + 181834528240062069095792640*x^21 + 1451179588383130489234391040*x^18 + 7545491276627959529519185920*x^15 + 26720234341607254534641418240*x^12 + 66151372999170901252002283520*x^9 + 115172590878530429476518494208*x^6 + 135614479302858341201788010496*x^3 + 75043585384980900269764116480,\n", " 3116285494907301262*x^35 + 1905929106580294155360*x^32 + 182789737266649816934400*x^29 + 6580736500578972004417536*x^26 + 118756184297100814024704000*x^23 + 1243240727966936007306117120*x^20 + 8236424690823173047006003200*x^17 + 36510441661103029981544448000*x^14 + 112224984234750469045493956608*x^11 + 243715584733787530928429465600*x^8 + 372096062838329079847213596672*x^5 + 369857670825977294186694574080*x^2,\n", " 218139984643511088340*x^34 + 44413035181030034061312*x^31 + 2551620806996826936053760*x^28 + 65488083989589123019702272*x^25 + 916921271665085717741568000*x^22 + 7828832965468067097005260800*x^19 + 43698305942081625792247234560*x^16 + 166856914993730628973127270400*x^13 + 448325896098926426212433657856*x^10 + 857154502345457074720420659200*x^7 + 1143671881451141327668924907520*x^4 + 913387639257196100426271817728*x,\n", " 11959798385860453492*x^36 + 7739227281265436873280*x^33 + 786745194635389174800384*x^30 + 30082266356173117561896960*x^27 + 577836035202256967820902400*x^24 + 6455125752522203452900638720*x^21 + 45768561951967161490184601600*x^18 + 217924019243627848106791403520*x^15 + 723046631639499392216884838400*x^12 + 1707908175614957814142604410880*x^9 + 2880039127880735771060613414912*x^6 + 3327045473312411135036872458240*x^3 + 1826775278514392200852543635456,\n", " 861105483781952651424*x^35 + 185510433040481964455040*x^32 + 11298176893452474051231744*x^29 + 308011153184505011503890432*x^26 + 4591241422559429513025945600*x^23 + 41842316188653437180958474240*x^20 + 250061914365425572549283020800*x^17 + 1026351832566763413664243384320*x^14 + 2980653408546830435679958204416*x^11 + 6216596517335331212696349573120*x^8 + 9241444234518404891754935353344*x^5 + 9051216877707056379789752401920*x^2,\n", " 45950804324621742364*x^37 + 31412157788665596720960*x^34 + 3379143887998933014073344*x^31 + 136982415894627171265216512*x^28 + 2795346503429186991195684864*x^25 + 33252310217855953324128337920*x^22 + 251734259996776776535847731200*x^19 + 1283991380700838245987747102720*x^16 + 4583191804048133174811017871360*x^13 + 11721023548585120650451623084032*x^10 + 21648383166250094575977876160512*x^7 + 28258961708857932313630794055680*x^4 + 22326334965010739070148055924736*x,\n", " 3400359520022008934936*x^36 + 773922728126543687328000*x^33 + 49882600251412820683939840*x^30 + 1441920167311864960686489600*x^27 + 22837798230630612673167360000*x^24 + 221682068119039688827522252800*x^21 + 1415074242585431584077250560000*x^18 + 6225412754913155132061804134400*x^15 + 19471076291707755317824585728000*x^12 + 44063998302951581392675274752000*x^9 + 72161277220833648586592920535040*x^6 + 81910033938718644387335634944000*x^3 + 44652669930021478140296111849472,\n", " 176733862787006701400*x^38 + 127443611599728992410752*x^35 + 14485061210010235580736000*x^32 + 621488134279897438029414400*x^29 + 13449182651472485906039439360*x^26 + 170025077439470357207580672000*x^23 + 1371335584178245516933044633600*x^20 + 7474214532265908907404165120000*x^17 + 28616816695971761494155067392000*x^14 + 78922762569055434888248987484160*x^11 + 158630393890625693013630989107200*x^8 + 230108688200700305982422180167680*x^5 + 222327234976522034765625294848000*x^2,\n", " 13431773571812509306400*x^37 + 3224981532969667930018560*x^34 + 219638097362862895205744640*x^31 + 6720498807975840092589260800*x^28 + 112894187288251223412540702720*x^25 + 1164852658202328830230659072000*x^22 + 7924415758478416870832603136000*x^19 + 37271416467565999084922103398400*x^16 + 125147782548712206453140619264000*x^13 + 305938580759679420367720331345920*x^10 + 547430378916669058243118707507200*x^7 + 700330790176044409511719678771200*x^4 + 547814306982150293662500726505472*x,\n", " 680425371729975800390*x^39 + 516854647043345358110272*x^36 + 61975732068373618481226240*x^33 + 2809973118642341167235399680*x^30 + 64375304371136994571118182400*x^27 + 863308491027803473154723020800*x^24 + 7403285783241467676577077657600*x^21 + 43018256974597120155948417024000*x^18 + 176192150573948359310540852428800*x^15 + 522221611129194145446438633472000*x^12 + 1136343299964523561365818373570560*x^9 + 1811496890233341247277229177569280*x^6 + 2023177838286350516367190183116800*x^3 + 1095628613964300587325001453010944,\n", " 53073178994938112430420*x^38 + 13424060421838120533932544*x^35 + 964601394043283721674966016*x^32 + 31192675836919665416149401600*x^29 + 554798077671253371394727608320*x^26 + 6073591461924732412249925222400*x^23 + 43942937797486922631275033395200*x^20 + 220439500604962540042373509939200*x^17 + 792363181898022033332033814528000*x^14 + 2084488788844109683718837155921920*x^11 + 4049824993858527579604495556935680*x^8 + 5743475090104372917788673111490560*x^5 + 5479921687701315112903132267413504*x^2,\n", " 2622127042276492108820*x^40 + 2095356677202751451798400*x^37 + 264699782965822095035289600*x^34 + 12663484967952852961476476928*x^31 + 306639864159549253243502592000*x^28 + 4354641515306441556607673303040*x^25 + 39629107411140098009006604288000*x^22 + 244978050438951485052404367360000*x^19 + 1070706145795532337348671333990400*x^16 + 3400179171363159683263899893760000*x^13 + 7975435366012245746402507379179520*x^10 + 13860363615344693320571535595929600*x^7 + 17404469970013251266026282156032000*x^4 + 13489038000495544893300017889017856*x,\n", " 209770163382119368705600*x^39 + 55820301880681298675909376*x^36 + 4225989578480429273780797440*x^33 + 144208665553422284744977022976*x^30 + 2711341956779172344468865024000*x^27 + 31437251581089283964814753792000*x^24 + 241432715920176597101024850739200*x^21 + 1288961434617252429044958363648000*x^18 + 4947005935654224275236641885388800*x^15 + 13958630282438234489188641669120000*x^12 + 29297517671065392537805129148006400*x^9 + 45558412579132991958052525697925120*x^6 + 50124873513638163646155692609372160*x^3 + 26978076000991089786600035778035712,\n", " 10113918591637898134020*x^41 + 8491708639190097988867200*x^38 + 1128640624327199956789619712*x^35 + 56893173131184286641048944640*x^32 + 1453907037956634510133784739840*x^29 + 21828767608396682075105698775040*x^26 + 210437112624434390621617127424000*x^23 + 1381219723636359136903537518182400*x^20 + 6427388775321164139426886975488000*x^17 + 21809381006647655406957238419456000*x^14 + 54941168791676890949446493609656320*x^11 + 103452478666250046377512848379084800*x^8 + 143684224288034820790781042585763840*x^5 + 135502405322582288537959344856104960*x^2,\n", " 829341324514307646989640*x^40 + 231886138943771160665689600*x^37 + 18471555006594455630838564864*x^34 + 664212166313958089312641744896*x^31 + 13180995255599319268896496680960*x^28 + 161600276154099622360473749422080*x^25 + 1314973643208004097206766272512000*x^22 + 7456340621906849324422348878643200*x^19 + 30484186762951807061281806798028800*x^16 + 91987298474499366907021092716544000*x^13 + 207707192988824022512814280188887040*x^10 + 351368009383069211225925991425638400*x^7 + 433665113305705095477630055804305408*x^4 + 333235544941461628108166685127606272*x,\n", " 39044429911904443959240*x^42 + 34402306794667576467718400*x^39 + 4804680798914938448993088512*x^36 + 254853918351853938558816141312*x^33 + 6863525718577566922897298030592*x^30 + 108777897688314382176998456819712*x^27 + 1109021503018330741689525731328000*x^24 + 7714512040153624036946362132070400*x^21 + 38142050104369652313391246186905600*x^18 + 137948643129317268317719691368857600*x^15 + 372037518274641883935063086098022400*x^12 + 755298883595723718228415564323225600*x^9 + 1147802163984692756671358238657085440*x^6 + 1245397248467665915217809391027748864*x^3 + 666471089882923256216333370255212544,\n", " 3279732112599973292576160*x^41 + 962393645774877772071616000*x^38 + 80560675313313488514623840256*x^35 + 3048465846383583358880508641280*x^32 + 63759528533234228245494025420800*x^29 + 825273637007211593706484160004096*x^26 + 7103308841938444938921923051520000*x^23 + 42698927105966570251005446219366400*x^20 + 185555919426663173416497954422784000*x^17 + 597230967707570516485883213905920000*x^14 + 1445631499581465606147673705980887040*x^11 + 2644410278950703178236786758385664000*x^8 + 3602332945428881882476262779785314304*x^5 + 3360595749833384215667104705947893760*x^2,\n", " 150853479205085351660700*x^43 + 139329342518403684694259520*x^40 + 20422743080469484150195072000*x^37 + 1138439590108326762577548165120*x^34 + 32266838497106236167842922233856*x^31 + 539038183599614119472074913218560*x^28 + 5802867460214229948200310005563392*x^25 + 42706215469161836137348218224640000*x^22 + 223909007994702746438199291150336000*x^19 + 861215092704068442968000029733683200*x^16 + 2479538224551430696031046170836992000*x^13 + 5405404737565480092552171248450273280*x^10 + 8918403293716096993269163185143808000*x^7 + 10832189975765169296956594372781015040*x^4 + 8257463842447744072782028706043396096*x,\n", " 12973399211637340242820200*x^42 + 3990667588181438870255334400*x^39 + 350617653205500103890548996096*x^36 + 13944117215488325315546117898240*x^33 + 306961775755539749363923672891392*x^30 + 4188507225673556602979184971612160*x^27 + 38073112295523228166925110896230400*x^24 + 242191693060668901827805895353958400*x^21 + 1116505551177205550112649859039232000*x^18 + 3823959837461126866384229797095014400*x^15 + 9891914398475019813795813824397312000*x^12 + 19472886632223468656399126236653158400*x^9 + 28966973897989883034138242025347088384*x^6 + 31025531782438509591283085117101178880*x^3 + 16514927684895488145564057412086792192,\n", " 583300119592996693088040*x^44 + 564113923367195406323099520*x^41 + 86683361788852520270356377600*x^38 + 5071948517602850817923558080512*x^35 + 151095956693201853419201217822720*x^32 + 2656980616375819469904126217814016*x^29 + 30157252024849607541450131795607552*x^26 + 234432691452317953027830420445593600*x^23 + 1301076304581732937726120042482892800*x^20 + 5310945324518599373508820951105536000*x^17 + 16282667694995766011700591394081996800*x^14 + 37984951290144076084975925085685678080*x^11 + 67642658827723627964333806927321497600*x^8 + 90500368154666616816715927747829956608*x^5 + 83578983577181291143864637458313379840*x^2,\n", " 51330410524183708991747520*x^43 + 16533748645517237250385463040*x^40 + 1522925348415060901581014384640*x^37 + 63577946206570946872562911150080*x^34 + 1471183879187141465086786558492672*x^31 + 21133026104363605706695068623044608*x^28 + 202565692846159628014268809796911104*x^25 + 1361372253625233616589982611779092480*x^22 + 6645382559987043412202277404934144000*x^19 + 24166797822546288277049161200068198400*x^16 + 66627927579522904599602419957392998400*x^13 + 140594583414181665358001628199948124160*x^10 + 226642696872184579344026684857660735488*x^7 + 271196389419708043659519110086089768960*x^4 + 205200814851562342394591799552824573952*x]" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Pol. = ZZ[]\n", "terms = [Pol({(i,j): w(k,i,j) for i in range(k+1) for j in range(k+1)})(y=0)\n", " for k in range(90)]\n", "terms" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "### Guessing differential operators" ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "scrolled": true, "slideshow": { "slide_type": "-" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 564 ms, sys: 4.14 ms, total: 568 ms\n", "Wall time: 566 ms\n" ] }, { "data": { "text/plain": [ "((31104*x^6 - 3888*x^3 - 972)*t^12 + (7776*x^4 + 3888*x)*t^11 + (-6912*x^8 - 8208*x^5 - 4320*x^2)*t^10 + (3168*x^6 - 936*x^3 + 36)*t^9 + (4572*x^4 - 144*x)*t^8 + (256*x^8 - 2504*x^5 + 160*x^2)*t^7 + (272*x^6 + 40*x^3)*t^6 - 180*x^4*t^5 + 104*x^5*t^4 - 16*x^6*t^3)*Dt^4 + ((497664*x^6 - 62208*x^3 - 15552)*t^11 + (77760*x^4 + 59292*x)*t^10 + (-124416*x^8 - 93312*x^5 - 57348*x^2)*t^9 + (80928*x^6 - 28692*x^3 + 414)*t^8 + (-12096*x^7 + 78228*x^4 - 1548*x)*t^7 + (3456*x^8 - 40248*x^5 + 1404*x^2)*t^6 + (3712*x^6 + 944*x^3)*t^5 + (448*x^7 - 2146*x^4)*t^4 + 1100*x^5*t^3 - 168*x^6*t^2)*Dt^3 + ((2239488*x^6 - 279936*x^3 - 69984)*t^10 + (139968*x^4 + 253692*x)*t^9 + (-635904*x^8 - 262656*x^5 - 205632*x^2)*t^8 + (525888*x^6 - 186696*x^3 + 1170)*t^7 + (-117504*x^7 + 362700*x^4 - 4005*x)*t^6 + (12672*x^8 - 171528*x^5 + 2400*x^2)*t^5 + (13008*x^6 + 4497*x^3)*t^4 + (3360*x^7 - 6486*x^4)*t^3 + 2856*x^5*t^2 - 432*x^6*t)*Dt^2 + ((2985984*x^6 - 373248*x^3 - 93312)*t^9 + (-93312*x^4 + 320760*x)*t^8 + (-967680*x^8 - 158976*x^5 - 204984*x^2)*t^7 + (953280*x^6 - 320904*x^3 + 720)*t^6 + (-253440*x^7 + 478080*x^4 - 2160*x)*t^5 + (13056*x^8 - 198744*x^5 + 135*x^2)*t^4 + (10656*x^6 + 4473*x^3)*t^3 + (5376*x^7 - 4506*x^4)*t^2 + 1608*x^5*t - 240*x^6)*Dt + (746496*x^6 - 93312*x^3 - 23328)*t^8 + (-93312*x^4 + 75816*x)*t^7 + (-276480*x^8 + 3456*x^5 - 34128*x^2)*t^6 + (311040*x^6 - 97200*x^3)*t^5 + (-92160*x^7 + 112680*x^4)*t^4 + (2304*x^8 - 35136*x^5)*t^3 - 288*x^6*t^2 + 1344*x^7*t" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from ore_algebra.guessing import guess\n", "%time dop = guess(terms, OreAlgebra(ZZ['x']['t'], 'Dt'))\n", "dop" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "### Guessing algebraic equations" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "scrolled": true, "slideshow": { "slide_type": "-" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 950 ms, sys: 189 ms, total: 1.14 s\n", "Wall time: 1.13 s\n" ] }, { "data": { "text/plain": [ "16*x^6*t^10*C^6 + (96*x^4*t^9 - 48*x^5*t^8)*C^5 + ((48*x^5 + 192*x^2)*t^8 - 192*x^3*t^7 + 56*x^4*t^6)*C^4 + ((192*x^3 + 128)*t^7 + (-96*x^4 - 192*x)*t^6 + 128*x^2*t^5 - 32*x^3*t^4)*C^3 + ((48*x^4 + 192*x)*t^6 - 264*x^2*t^5 + (64*x^3 + 32)*t^4 - 32*x*t^3 + 9*x^2*t^2)*C^2 + (96*x^2*t^5 + (-48*x^3 - 144)*t^4 + 104*x*t^3 - 16*x^2*t^2 + 2*t - x)*C + (16*x^3 + 108)*t^4 - 72*x*t^3 + 8*x^2*t^2 - 2*t + x" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%time algeq = guess(terms, OreAlgebra(ZZ['x']['t'], 'C'))\n", "algeq" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### The Voigt profile\n", "\n", "$$V(x) = \\frac{1}{σ \\sqrt{2π}} \\frac{λ}{π} \\int_{-∞}^{+∞} \\frac{\\exp\\bigl(-(u-x)²/(2σ²)\\bigr)}{u² + λ²} \\mathrm du,\n", "\\qquad σ = 1, λ = 1/2$$" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [], "source": [ "P. = PolynomialRing(QQ)\n", "A. = OreAlgebra(P)\n", "f = exp(-(u-x)^2/2)/(u^2 + 1/4)" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "op1 = Dx + x - u\n", "op1(f)" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "op2 = (4*u^2 + 1)*Du + (4*u^3 - 4*u^2*x + 9*u - x)\n", "op2(f).simplify_full()" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "### Creative telescoping\n", "\n", "$$V(x) = \\frac{1}{(2π)^{3/2}} \\int_{-∞}^{+∞} f(x,u) \\, \\mathrm du, \\qquad f(x,u) = \\frac{\\exp\\bigl(-(u-x)²/2\\bigr)}{u² + 1/4}$$" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Left Ideal (Dx + x - u, (4*u^2 + 1)*Du - 4*x*u^2 + 4*u^3 - x + 9*u) of Multivariate Ore algebra in Dx, Du over Fraction Field of Multivariate Polynomial Ring in x, u over Rational Field" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ideal = A.ideal([op1, op2]); ideal" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Idea:**\n", "Find a pair $(P,Q)$ such that\n", "$P - D_u Q \\in I$\n", "and $P$ does not depend on $u$.\n", "\n", "**Then**\n", "$P\\left(\\int f \\, \\mathrm du \\right) = \\int P(f) \\, \\mathrm du = \\int (D_u Q(f)) \\, \\mathrm du = 0$." ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(-4*Dx^3 - 8*x*Dx^2 + (-4*x^2 - 13)*Dx - 8*x, 4*u^2 + 1)" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "[tel], [cert] = ideal.ct(Du)\n", "tel, cert" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ "(tel-Du*cert).reduce(ideal)" ] }, { "cell_type": "code", "execution_count": 31, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "data": { "text/plain": [ "[0.08242408278858694 +/- 2.82e-18]" ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ini = [(1/2*exp(1/8)*erfc(1/4*2^(1/2))*2^(1/2)/pi^(1/2)), 0,\n", " 1/2*(1/2/pi-5/8*exp(1/8)*erfc(1/4*2^(1/2))*2^(1/2)/pi^(1/2))]\n", "tel.numerical_solution(ini, [0,2])" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Another walk model\n", "After A. Bostan, F. Chyzak, M. van Hoeij, M. Kauers and L. Pech, *[Hypergeometric expressions for generating functions of walks with small steps in the quarter plane](https://arxiv.org/pdf/1606.02982)*, 2017." ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(-u^3*v^2 - u^2*v^2 + u^3 - u*v^2 + 2*u^2 - v^2 + 2*u + 1)/(u^4*v^3*t + u^3*v^3*t + 2*u^2*v^3*t - u^3*v^2 - u^4*t + u^3*v*t + u*v^3*t + u^3*v - u^2*v^2 - 2*u^3*t + u^2*v*t + v^3*t + u^2*v - u*v^2 - 3*u^2*t + u*v*t + u*v - 2*u*t - t)" ] }, "execution_count": 32, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from ore_algebra.examples import ssw\n", "Rat = Frac(PolynomialRing(ZZ,'u,v,t'))\n", "q = Rat(ssw.rat[10](x=1,y=1))\n", "q" ] }, { "cell_type": "code", "execution_count": 33, "metadata": { "scrolled": true, "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "A. = OreAlgebra(Rat)" ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Left Ideal ((-u^7*v^4*t - u^7*v^3*t - 2*u^6*v^4*t - 2*u^6*v^3*t - 4*u^5*v^4*t + u^6*v^3 + u^7*v*t - 4*u^5*v^3*t - 5*u^4*v^4*t + 2*u^5*v^3 + u^7*t + 3*u^6*v*t - 5*u^4*v^3*t - 5*u^3*v^4*t - u^6*v + 3*u^4*v^3 + 4*u^6*t + 6*u^5*v*t - 5*u^3*v^3*t - 4*u^2*v^4*t - 3*u^5*v + 3*u^3*v^3 + 9*u^5*t + 8*u^4*v*t - 4*u^2*v^3*t - 2*u*v^4*t - 5*u^4*v + 2*u^2*v^3 + 13*u^4*t + 8*u^3*v*t - 2*u*v^3*t - v^4*t - 5*u^3*v + u*v^3 + 13*u^3*t + 6*u^2*v*t - v^3*t - 3*u^2*v + 9*u^2*t + 3*u*v*t - u*v + 4*u*t + v*t + t)*Du - u^6*v^4*t - u^6*v^3*t - 2*u^5*v^4*t - 2*u^5*v^3*t - 2*u^4*v^4*t + u^6*v*t + 2*u^5*v^2*t - 2*u^4*v^3*t - 4*u^3*v^4*t + u^6*t + 4*u^5*v*t + 4*u^4*v^2*t - 4*u^3*v^3*t - u^2*v^4*t + 4*u^5*t + 6*u^4*v*t + 2*u^3*v^2*t - u^2*v^3*t - 2*u*v^4*t - u^4*v + 3*u^2*v^3 + 7*u^4*t + 6*u^3*v*t - 2*u^2*v^2*t - 2*u*v^3*t - 2*u^3*v + 2*u*v^3 + 8*u^3*t + 2*u^2*v*t - 4*u*v^2*t - 3*u^2*v + v^3 + 5*u^2*t - 2*v^2*t - 2*u*v + 2*u*t - v*t - v, (-u^5*v^5*t - u^4*v^5*t + u^5*v^3*t - 2*u^3*v^5*t + u^4*v^4 + u^5*v^2*t + u^4*v^3*t - 2*u^2*v^5*t - u^4*v^3 + u^3*v^4 + 2*u^4*v^2*t + 2*u^3*v^3*t - u*v^5*t - u^4*v^2 - u^3*v^3 + u^2*v^4 - u^5*t + u^4*v*t + 3*u^3*v^2*t + 2*u^2*v^3*t - v^5*t + u^4*v - 2*u^3*v^2 - u^2*v^3 + u*v^4 - 3*u^4*t + 2*u^3*v*t + 3*u^2*v^2*t + u*v^3*t + 2*u^3*v - 2*u^2*v^2 - u*v^3 - 5*u^3*t + 2*u^2*v*t + 2*u*v^2*t + v^3*t + 2*u^2*v - u*v^2 - 5*u^2*t + u*v*t + v^2*t + u*v - 3*u*t - t)*Dv - u^5*v^4*t - u^4*v^4*t + 3*u^5*v^2*t - 2*u^3*v^4*t - 2*u^5*v*t + 7*u^4*v^2*t - 2*u^2*v^4*t + u^4*v^2 - 4*u^4*v*t + 10*u^3*v^2*t - u*v^4*t - 2*u^4*v + u^3*v^2 + u^4*t - 6*u^3*v*t + 10*u^2*v^2*t - v^4*t + u^4 - 4*u^3*v + u^2*v^2 + 2*u^3*t - 6*u^2*v*t + 7*u*v^2*t + 2*u^3 - 4*u^2*v + u*v^2 + 2*u^2*t - 4*u*v*t + 3*v^2*t + 2*u^2 - 2*u*v + u*t - 2*v*t + u, (-u^5*v^4*t - u^5*v^3*t - u^4*v^4*t - u^4*v^3*t - 2*u^3*v^4*t + u^4*v^3 + u^5*v*t - 2*u^3*v^3*t - 2*u^2*v^4*t + u^3*v^3 + u^5*t + 2*u^4*v*t - 2*u^2*v^3*t - u*v^4*t - u^4*v + u^2*v^3 + 3*u^4*t + 3*u^3*v*t - u*v^3*t - v^4*t - 2*u^3*v + u*v^3 + 5*u^3*t + 3*u^2*v*t - v^3*t - 2*u^2*v + 5*u^2*t + 2*u*v*t - u*v + 3*u*t + v*t + t)*Dt - u^5*v^4 - u^5*v^3 - u^4*v^4 - u^4*v^3 - 2*u^3*v^4 + u^5*v - 2*u^3*v^3 - 2*u^2*v^4 + u^5 + 2*u^4*v - 2*u^2*v^3 - u*v^4 + 3*u^4 + 3*u^3*v - u*v^3 - v^4 + 5*u^3 + 3*u^2*v - v^3 + 5*u^2 + 2*u*v + 3*u + v + 1) of Multivariate Ore algebra in Du, Dv, Dt over Fraction Field of Multivariate Polynomial Ring in u, v, t over Integer Ring" ] }, "execution_count": 34, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ideal = A.ideal([q*D - D(q) for D in Du,Dv,Dt])\n", "ideal" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "### ⚐ Fast creative telescoping" ] }, { "cell_type": "code", "execution_count": 35, "metadata": { "scrolled": true, "slideshow": { "slide_type": "-" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 1.87 s, sys: 52.6 ms, total: 1.93 s\n", "Wall time: 1.96 s\n" ] }, { "data": { "text/plain": [ "[(-9*u^9*t^4 - 24*u^8*t^4 - 3*u^8*t^3 - 37*u^7*t^4 - u^7*t^3 - 37*u^6*t^4 + 5*u^7*t^2 - u^6*t^3 - 15*u^5*t^4 + 6*u^6*t^2 - 3*u^5*t^3 + 15*u^4*t^4 - u^6*t + u^5*t^2 + 3*u^4*t^3 + 37*u^3*t^4 - u^5*t - u^4*t^2 + u^3*t^3 + 37*u^2*t^4 + u^4*t - 6*u^3*t^2 + u^2*t^3 + 24*u*t^4 + u^3*t - 5*u^2*t^2 + 3*u*t^3 + 9*t^4)*Dt^2 + (-2*u^8 - 4*u^7 - 6*u^6 - 6*u^5 - 4*u^4 - 2*u^3)*Du + (-36*u^9*t^3 - 96*u^8*t^3 - 12*u^8*t^2 - 148*u^7*t^3 - 6*u^7*t^2 - 148*u^6*t^3 + 16*u^7*t - 6*u^6*t^2 - 60*u^5*t^3 + 20*u^6*t - 12*u^5*t^2 + 60*u^4*t^3 - 2*u^6 + 4*u^5*t + 12*u^4*t^2 + 148*u^3*t^3 - 2*u^5 - 4*u^4*t + 6*u^3*t^2 + 148*u^2*t^3 + 2*u^4 - 20*u^3*t + 6*u^2*t^2 + 96*u*t^3 + 2*u^3 - 16*u^2*t + 12*u*t^2 + 36*t^3)*Dt - 18*u^9*t^2 - 48*u^8*t^2 - 6*u^8*t - 74*u^7*t^2 - 4*u^7*t - 74*u^6*t^2 + 4*u^7 - 4*u^6*t - 30*u^5*t^2 + 4*u^6 - 6*u^5*t + 30*u^4*t^2 + 6*u^4*t + 74*u^3*t^2 - 6*u^4 + 4*u^3*t + 74*u^2*t^2 - 8*u^3 + 4*u^2*t + 48*u*t^2 - 6*u^2 + 6*u*t + 18*t^2,\n", " (-9*u^10*t^3 - 33*u^9*t^3 - 3*u^9*t^2 - 79*u^8*t^3 - 4*u^8*t^2 - 131*u^7*t^3 + 5*u^8*t - 8*u^7*t^2 - 168*u^6*t^3 + 11*u^7*t - 9*u^6*t^2 - 168*u^5*t^3 - u^7 + 17*u^6*t - 9*u^5*t^2 - 131*u^4*t^3 - 2*u^6 + 17*u^5*t - 8*u^4*t^2 - 79*u^3*t^3 - 2*u^5 + 11*u^4*t - 4*u^3*t^2 - 33*u^2*t^3 - u^4 + 5*u^3*t - 3*u^2*t^2 - 9*u*t^3)*Du*Dt + (-9*u^10*t^2 - 33*u^9*t^2 - 6*u^9*t - 79*u^8*t^2 - 15*u^8*t - 131*u^7*t^2 + u^8 - 30*u^7*t - 168*u^6*t^2 + 2*u^7 - 39*u^6*t - 168*u^5*t^2 + 3*u^6 - 39*u^5*t - 131*u^4*t^2 + 3*u^5 - 30*u^4*t - 79*u^3*t^2 + 2*u^4 - 15*u^3*t - 33*u^2*t^2 + u^3 - 6*u^2*t - 9*u*t^2)*Du + (-18*u^9*t^3 - 60*u^8*t^3 - 3*u^8*t^2 - 106*u^7*t^3 - u^7*t^2 - 137*u^6*t^3 + 8*u^7*t + 4*u^6*t^2 - 109*u^5*t^3 + 17*u^6*t - 7*u^5*t^2 - 59*u^4*t^3 - u^6 + 10*u^5*t - 2*u^4*t^2 + 6*u^3*t^3 - 2*u^5 + 7*u^4*t - 12*u^3*t^2 + 27*u^2*t^3 - 6*u^3*t - 3*u^2*t^2 + 27*u*t^3 - 3*u^2*t + 9*t^3)*Dt - 18*u^9*t^2 - 60*u^8*t^2 - 6*u^8*t - 106*u^7*t^2 - 12*u^7*t - 137*u^6*t^2 + 4*u^7 - 12*u^6*t - 109*u^5*t^2 + 8*u^6 - 24*u^5*t - 59*u^4*t^2 + 4*u^5 - 15*u^4*t + 6*u^3*t^2 - u^4 - 18*u^3*t + 27*u^2*t^2 - 6*u^3 - 3*u^2*t + 27*u*t^2 - 3*u^2 + 9*t^2,\n", " (9*u^12*t^3 + 33*u^11*t^3 + 3*u^11*t^2 + 70*u^10*t^3 + 4*u^10*t^2 + 98*u^9*t^3 - 5*u^10*t + 5*u^9*t^2 + 89*u^8*t^3 - 11*u^9*t + 5*u^8*t^2 + 37*u^7*t^3 + u^9 - 12*u^8*t + u^7*t^2 - 37*u^6*t^3 + 2*u^8 - 6*u^7*t - u^6*t^2 - 89*u^5*t^3 + u^7 + 6*u^6*t - 5*u^5*t^2 - 98*u^4*t^3 - u^6 + 12*u^5*t - 5*u^4*t^2 - 70*u^3*t^3 - 2*u^5 + 11*u^4*t - 4*u^3*t^2 - 33*u^2*t^3 - u^4 + 5*u^3*t - 3*u^2*t^2 - 9*u*t^3)*Du^2 + (36*u^11*t^3 + 126*u^10*t^3 + 18*u^10*t^2 + 172*u^9*t^3 + 34*u^9*t^2 + 98*u^8*t^3 - 14*u^9*t + 8*u^8*t^2 - 158*u^7*t^3 - 30*u^8*t - 8*u^7*t^2 - 442*u^6*t^3 + 2*u^8 + 6*u^7*t - 64*u^6*t^2 - 590*u^5*t^3 + 4*u^7 + 36*u^6*t - 68*u^5*t^2 - 514*u^4*t^3 - 4*u^6 + 60*u^5*t - 28*u^4*t^2 - 294*u^3*t^3 - 8*u^5 + 54*u^4*t - 12*u^3*t^2 - 108*u^2*t^3 - 4*u^4 + 14*u^3*t + 18*u^2*t^2 - 6*u*t^3 - 2*u^3 + 6*u^2*t + 6*u*t^2)*Du + (-6*u^9*t^5 - 6*u^9*t^4 - 12*u^8*t^5 - 6*u^9*t^3 - 12*u^8*t^4 + 6*u^7*t^5 - 12*u^8*t^3 + 6*u^7*t^4 + 30*u^6*t^5 + 6*u^7*t^3 + 30*u^6*t^4 + 18*u^5*t^5 + 30*u^6*t^3 + 18*u^5*t^4 - 18*u^4*t^5 + 18*u^5*t^3 - 18*u^4*t^4 - 30*u^3*t^5 - 18*u^4*t^3 - 30*u^3*t^4 - 6*u^2*t^5 - 30*u^3*t^3 - 6*u^2*t^4 + 12*u*t^5 - 6*u^2*t^3 + 12*u*t^4 + 6*t^5 + 12*u*t^3 + 6*t^4 + 6*t^3)*Dt + 18*u^10*t^3 - 6*u^9*t^4 + 54*u^9*t^3 - 12*u^8*t^4 + 6*u^9*t^2 + 38*u^8*t^3 + 6*u^7*t^4 + 14*u^8*t^2 - 68*u^7*t^3 + 30*u^6*t^4 - 4*u^8*t - 2*u^7*t^2 - 148*u^6*t^3 + 18*u^5*t^4 - 8*u^7*t - 16*u^6*t^2 - 272*u^5*t^3 - 18*u^4*t^4 + 8*u^6*t + 10*u^5*t^2 - 244*u^4*t^3 - 30*u^3*t^4 + 46*u^5*t - 76*u^4*t^2 - 188*u^3*t^3 - 6*u^2*t^4 + 2*u^4*t - 2*u^3*t^2 - 30*u^2*t^3 + 12*u*t^4 - 6*u^4 + 22*u^3*t - 6*u*t^3 + 6*t^4 + 12*u*t^2 + 6*t^3 + 6*t^2]" ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%time tel1, _ = ideal.ct(Dv)\n", "tel1" ] }, { "cell_type": "code", "execution_count": 36, "metadata": { "scrolled": true, "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 12.9 s, sys: 90.3 ms, total: 13 s\n", "Wall time: 13.1 s\n" ] }, { "data": { "text/plain": [ "[(1892352000*t^24 + 9122611200*t^23 - 18640527360*t^22 - 58034765824*t^21 - 38532706304*t^20 + 41965621248*t^19 + 94071750144*t^18 + 56013030400*t^17 - 5294450432*t^16 - 21996570880*t^15 - 10952369952*t^14 - 2543280896*t^13 - 537441232*t^12 - 276140992*t^11 - 112877022*t^10 - 19917936*t^9 + 1121240*t^8 + 1056920*t^7 + 134035*t^6 - 5391*t^5 - 2254*t^4 - 103*t^3 + 3*t^2)*Dt^5 + (47308800000*t^23 + 227888332800*t^22 - 542727905280*t^21 - 1484019662848*t^20 - 767620100096*t^19 + 1085290090496*t^18 + 1896743070208*t^17 + 918495748096*t^16 - 215512785664*t^15 - 427218085376*t^14 - 200103936864*t^13 - 53308965120*t^12 - 16198105488*t^11 - 7684582384*t^10 - 2788409498*t^9 - 526917856*t^8 - 11674372*t^7 + 14725960*t^6 + 2406665*t^5 + 42072*t^4 - 17460*t^3 - 836*t^2 + 27*t)*Dt^4 + (378470400000*t^22 + 1821691084800*t^21 - 4964213882880*t^20 - 12009478135808*t^19 - 4381901037568*t^18 + 8936076752896*t^17 + 11847464069120*t^16 + 4115996062720*t^15 - 2219866571776*t^14 - 2632683934976*t^13 - 1121810512000*t^12 - 331856938400*t^11 - 129562335520*t^10 - 61796701736*t^9 - 21296426392*t^8 - 4208335952*t^7 - 300047680*t^6 + 45411880*t^5 + 10529090*t^4 + 535888*t^3 - 18106*t^2 - 1020*t + 48)*Dt^3 + (1135411200000*t^21 + 5460826521600*t^20 - 16785832181760*t^19 - 36057594863616*t^18 - 7287121182720*t^17 + 27418579642368*t^16 + 26356460402688*t^15 + 4411389139968*t^14 - 7504283639808*t^13 - 5923884943104*t^12 - 2207397164160*t^11 - 718441676064*t^10 - 351342324480*t^9 - 173136422904*t^8 - 57960776808*t^7 - 11637201048*t^6 - 1161271776*t^5 - 2986416*t^4 + 10037526*t^3 + 711192*t^2 + 15354*t + 792)*Dt^2 + (1135411200000*t^20 + 5456579788800*t^19 - 18705041326080*t^18 - 35704801148928*t^17 - 845976502272*t^16 + 27816503623680*t^15 + 17768620793856*t^14 - 1927605172224*t^13 - 7625766414336*t^12 - 4160593317888*t^11 - 1282076116992*t^10 - 484808936256*t^9 - 302844173952*t^8 - 153142047744*t^7 - 49288460160*t^6 - 9650682192*t^5 - 1077476208*t^4 - 60418992*t^3 - 1875024*t^2 - 158760*t - 9216)*Dt + 227082240000*t^19 + 1090466611200*t^18 - 4130053816320*t^17 - 6994180227072*t^16 + 1235817283584*t^15 + 5584717234176*t^14 + 1907260735488*t^13 - 1376741382144*t^12 - 1399425761280*t^11 - 495281498112*t^10 - 109934770176*t^9 - 64845759168*t^8 - 53271954240*t^7 - 26356354176*t^6 - 7768879584*t^5 - 1370419584*t^4 - 140485008*t^3 - 8567520*t^2 - 451296*t - 22464]" ] }, "execution_count": 36, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%time tel2, _ = tel1[0].parent().ideal(tel1).ct(Du)\n", "tel2 # HolonomicFunctions.m: ≥ 1 min" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Pólya walks in dimension 15\n", "\n", "(Thanks to B. Salvy)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Task:** Compute a certain solution at $z=1$" ] }, { "cell_type": "code", "execution_count": 37, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "(269276305858560000*z^30 - 80950584950784000*z^28 + 3629735201193984*z^26 - 57080630763520*z^24 + 409981265408*z^22 - 1499829760*z^20 + 2871232*z^18 - 2720*z^16 + z^14)*Dz^15 + (60587168818176000000*z^29 - 16999622839664640000*z^27 + 707798364232826880*z^25 - 10274513537433600*z^23 + 67646908792320*z^21 - 224974464000*z^19 + 387616320*z^17 - 326400*z^15 + 105*z^13)*Dz^14 + (5937542544181248000000*z^28 - 1551046657578762240000*z^26 + 59789931630571929600*z^24 - 798278106947641344*z^22 + 4796169345443840*z^20 - 14416495588352*z^18 + 22181540160*z^16 - 16423904*z^14 + 4550*z^12)*Dz^13 + (334481563322210304000000*z^27 - 81131900756901396480000*z^25 + 2886264219850496409600*z^23 - 35303727569401454592*z^21 + 192591205900620800*z^19 - 519787439069184*z^17 + 707833651200*z^15 - 454991264*z^13 + 106470*z^11)*Dz^12 + (12041336279599570944000000*z^26 - 2704250458682470563840000*z^24 + 88474218684064461619200*z^22 - 987113006841956179968*z^20 + 4862067116732345856*z^18 - 11694546001056256*z^16 + 13948983786816*z^14 - 7666274304*z^12 + 1479478*z^10)*Dz^11 + (291400337966309616844800000*z^25 - 60403894924097334804480000*z^23 + 1810516178566181073715200*z^21 - 18336647935059261603840*z^19 + 81033820655292578304*z^17 - 172214449187123200*z^15 + 177735512066112*z^13 - 81994323840*z^11 + 12662650*z^9)*Dz^10 + (4856672299438493614080000000*z^24 - 926086890189285634867200000*z^22 + 25324208695139969934950400*z^20 - 231569354145921664204800*z^18 + 911597715337047246336*z^16 - 1694812596013786624*z^14 + 1491624282894720*z^12 - 564676102848*z^10 + 67128490*z^8)*Dz^9 + (56198636607788283248640000000*z^23 - 9821600336936930947891200000*z^21 + 243895079188460801654784000*z^19 - 2001474041629081060638720*z^17 + 6961148053301631762432*z^15 - 11190413197947013632*z^13 + 8252818589658240*z^11 - 2491958589120*z^9 + 216627840*z^7)*Dz^8 + (449589092862306265989120000000*z^22 - 71725429608240620136038400000*z^20 + 1609161818725523357171712000*z^18 - 11770087561727055499100160*z^16 + 35825937728946311725056*z^14 - 49112927471281826304*z^12 + 29704887544668864*z^10 - 6898246766112*z^8 + 408741333*z^6)*Dz^7 + (2447762838917000781496320000000*z^21 - 354907365563571111670579200000*z^19 + 7152475405487676025208832000*z^17 - 46268511513163153122263040*z^15 + 121891943669163953209344*z^13 - 140140992502944986112*z^11 + 67667087104034880*z^9 - 11517430544256*z^7 + 420693273*z^5)*Dz^6 + (8811946220101202813386752000000*z^20 - 1155571219953730314672537600000*z^18 + 20785112953613635862082355200*z^16 - 117852071549595952781721600*z^14 + 265216641208694799482880*z^12 - 250559836896574725120*z^10 + 93303962552021952*z^8 - 10897207743264*z^6 + 210766920*z^4)*Dz^5 + (20027150500230006394060800000000*z^19 - 2362573936582134301458432000000*z^17 + 37651520030507075820847104000*z^15 - 185180175371863397892096000*z^13 + 350420651229762451537920*z^11 - 265169341550020423680*z^9 + 72752913986864640*z^7 - 5304232959840*z^5 + 42355950*z^3)*Dz^4 + (26702867333640008525414400000000*z^18 - 2816819552897219199762432000000*z^16 + 39443088818742150876364800000*z^14 - 166217350344029675008819200*z^12 + 259441279222310968688640*z^10 - 152179348773380567040*z^8 + 28891552538695680*z^6 - 1138125841504*z^4 + 2375101*z^2)*Dz^3 + (18486600461750775132979200000000*z^17 - 1732080379790685408067584000000*z^15 + 21106053778472440493506560000*z^13 - 75098351566004361206169600*z^11 + 94387689349096767160320*z^9 - 41091696296267489280*z^7 + 4930837635294720*z^5 - 82548913344*z^3 + 16383*z)*Dz^2 + (5281885846214507180851200000000*z^16 - 436217453243899431616512000000*z^14 + 4573788149507776189562880000*z^12 - 13497729262079414265446400*z^10 + 13245706827488316948480*z^8 - 4031063063164477440*z^6 + 265858264373760*z^4 - 1204325184*z^2 + 1)*Dz + 352125723080967145390080000000*z^15 - 25412342171473281024000000000*z^13 + 226231973017884794290176000*z^11 - 541574695562332078080000*z^9 + 398335457415580876800*z^7 - 77667722566041600*z^5 + 2222757642240*z^3 - 983040*z" ] }, "execution_count": 37, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from ore_algebra.examples import polya\n", "dim = 15\n", "polya.dop[dim]" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "### ⚐ Numerics: Logarithms" ] }, { "cell_type": "code", "execution_count": 38, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(2*z - 1) * (2*z + 1) * (6*z - 1) * (6*z + 1) * (10*z - 1) * (10*z + 1) * (14*z - 1) * (14*z + 1) * (18*z - 1) * (18*z + 1) * (22*z - 1) * (22*z + 1) * (26*z - 1) * (26*z + 1) * (30*z - 1) * (30*z + 1) * z^14" ] }, "execution_count": 38, "metadata": {}, "output_type": "execute_result" } ], "source": [ "polya.dop[dim].leading_coefficient().factor()" ] }, { "cell_type": "code", "execution_count": 39, "metadata": { "scrolled": true, "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/plain": [ "[1/87178291200*log(z)^14,\n", " 1/6227020800*log(z)^13,\n", " 1/479001600*log(z)^12,\n", " 1/39916800*log(z)^11,\n", " 1/3628800*log(z)^10,\n", " 1/362880*log(z)^9,\n", " 1/40320*log(z)^8,\n", " 1/5040*log(z)^7,\n", " 1/720*log(z)^6,\n", " 1/120*log(z)^5,\n", " 1/24*log(z)^4,\n", " 1/6*log(z)^3,\n", " 1/2*log(z)^2,\n", " log(z),\n", " 1]" ] }, "execution_count": 39, "metadata": {}, "output_type": "execute_result" } ], "source": [ "polya.dop[dim].local_basis_monomials(0)" ] }, { "cell_type": "code", "execution_count": 40, "metadata": { "scrolled": false, "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 4.3 s, sys: 43.5 ms, total: 4.34 s\n", "Wall time: 4.31 s\n" ] }, { "data": { "text/plain": [ "[1.03720412092152168870320272201726012686492785463956609556680543456793165591178031741211 +/- 7.28e-87]" ] }, "execution_count": 40, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ini = [0]*(dim - 1) + [1]\n", "# 90 s in Nov. 2017\n", "%time polya.dop[dim].numerical_solution(ini, [0,1/(2*dim)], 1e-100)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Iterated integrals\n", "After J. Ablinger, J. Blümlein, C. G. Raab, and C. Schneider,\n", "*[Iterated Binomial Sums and their Associated Iterated Integrals](http://arxiv.org/pdf/1407.1822)*,\n", "Journal of Mathematical Physics 55(11), 2014.\n", "\n", "$$\\int_{0}^1 dx_1 \\, w_1(x_1) \\int_{x_1}^1 dx_2 \\, w_2(x_2) \\, \\cdots \\! \\int_{x_{n-1}}^1 dx_n \\, w_n(x_n)$$" ] }, { "cell_type": "code", "execution_count": 41, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "[x/(x - 1), 2/(sqrt(4*x - 1)*x), 2/(sqrt(4*x - 1)*x), -1/(x - 1), -1/(x - 1)]" ] }, "execution_count": 41, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from ore_algebra.examples import iint\n", "i = 69; iint.word[i]" ] }, { "cell_type": "code", "execution_count": 42, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "(x^9 - 13/4*x^8 + 15/4*x^7 - 7/4*x^6 + 1/4*x^5)*Dx^6 + (27/2*x^8 - 35*x^7 + 30*x^6 - 9*x^5 + 1/2*x^4)*Dx^5 + (101/2*x^7 - 397/4*x^6 + 57*x^5 - 17/2*x^4 + 1/4*x^3)*Dx^4 + (111/2*x^6 - 303/4*x^5 + 45/2*x^4 + 3/4*x^3 - 3/4*x^2)*Dx^3 + (12*x^5 - 37/4*x^4 + 1/4*x^3 - 3/2*x^2 + 3/2*x)*Dx^2 + (-1/4*x^2 + 3/2*x - 3/2)*Dx" ] }, "execution_count": 42, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dop = iint.diffop(iint.word[i])\n", "dop" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "### ⚐ Numerics: Advanced interface" ] }, { "cell_type": "code", "execution_count": 43, "metadata": { "scrolled": true, "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/plain": [ "[0, 0, 0, -1/3, 1/3*I*pi + 2/3, -2/3*I*pi + 1/6*pi^2 - 11/12]" ] }, "execution_count": 43, "metadata": {}, "output_type": "execute_result" } ], "source": [ "iint.ini[i]" ] }, { "cell_type": "code", "execution_count": 44, "metadata": { "scrolled": true, "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/plain": [ "[0.25000000000000000?]" ] }, "execution_count": 44, "metadata": {}, "output_type": "execute_result" } ], "source": [ "roots = dop.leading_coefficient().roots(AA, multiplicities=False)\n", "sing = list(reversed(sorted([s for s in roots if 0 < s < 1])))\n", "sing" ] }, { "cell_type": "code", "execution_count": 45, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 745 ms, sys: 41.5 ms, total: 786 ms\n", "Wall time: 732 ms\n" ] }, { "data": { "text/plain": [ "[-0.97080469562493124056011987537954344846933233520382808407829772093922068543247104239693854364491332570 +/- 9.81e-102] + [+/- 2.77e-102]*I" ] }, "execution_count": 45, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from ore_algebra.analytic.path import Point\n", "path = [1] + [Point(s, dop, outgoing_branch=(0,-1)) for s in sing] + [0]\n", "%time dop.numerical_solution(iint.ini[i], path, 1e-100)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Bonus examples" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Transition matrices\n", "" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Compacted binary trees of fixed right height\n", "After A. Genitrini, B. Gittenberger, M. Kauers, M. Wallner, *[Asymptotic Enumeration of Compacted Binary Trees](https://arxiv.org/abs/1703.10031)*, 2017" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Exponential generating function of compacted binary trees:" ] }, { "cell_type": "code", "execution_count": 46, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "1 + z + 3/2*z^2 + 5/2*z^3 + 37/8*z^4 + 373/40*z^5 + 4829/240*z^6 + 76981/1680*z^7 + 293057/2688*z^8 + 32536277/120960*z^9 + 827662693/1209600*z^10" ] }, "execution_count": 46, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from ore_algebra.examples import cbt\n", "QQ[['z']](cbt.egf)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Annihilator of EGF of CBT of right height ≤ 5:" ] }, { "cell_type": "code", "execution_count": 47, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "(-4*z^3 + 10*z^2 - 6*z + 1)*Dz^6 + (9*z^3 - 58*z^2 + 75*z - 21)*Dz^5 + (-6*z^3 + 78*z^2 - 184*z + 95)*Dz^4 + (z^3 - 33*z^2 + 141*z - 110)*Dz^3 + (3*z^2 - 32*z + 40)*Dz^2 + (z - 3)*Dz" ] }, "execution_count": 47, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dop = cbt.dop[5]; dop" ] }, { "cell_type": "code", "execution_count": 48, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "[0.2928932188134525?, 0.50000000000000000?, 1.707106781186548?]" ] }, "execution_count": 48, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dop.leading_coefficient().roots(QQbar, multiplicities=False)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "### ⚐ Numerics: Algebraic points and exponents" ] }, { "cell_type": "code", "execution_count": 49, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/plain": [ "[1,\n", " z - 0.2928932188134525?,\n", " (z - 0.2928932188134525?)^1.771446609406727?,\n", " (z - 0.2928932188134525?)^2,\n", " (z - 0.2928932188134525?)^3,\n", " (z - 0.2928932188134525?)^4]" ] }, "execution_count": 49, "metadata": {}, "output_type": "execute_result" } ], "source": [ "s = _[0]\n", "dop.local_basis_monomials(s)" ] }, { "cell_type": "code", "execution_count": 50, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "0.292893218813452" ] }, "execution_count": 50, "metadata": {}, "output_type": "execute_result" } ], "source": [ "1/(4*(cos(pi/8)^2)).n()" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "-" } }, "source": [ "Singularity analysis implies\n", "\n", "$$\\#\\{\\text{cbt of rh ≤ 5}\\} \\sim \\kappa \\, n! α^n n^β, \\qquad α = 4 \\cos²(\\pi/8), \\quad β=-(α+5)/8$$\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Singularity analysis implies\n", "\n", "$$\\#\\{\\text{cbt of rh ≤ 5}\\} \\sim \\kappa_k n! α^n n^β, \\qquad α = 4 \\cos²(\\pi/8), \\quad β=-(α+5)/8$$\n" ] }, { "cell_type": "code", "execution_count": 51, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "[1, z, z^2, z^3, z^4, z^5]" ] }, "execution_count": 51, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dop.local_basis_monomials(0)" ] }, { "cell_type": "code", "execution_count": 52, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 671 ms, sys: 16.5 ms, total: 688 ms\n", "Wall time: 662 ms\n" ] }, { "data": { "text/plain": [ "[31.42732179257817 +/- 3.99e-15] + [27.45408401646515 +/- 4.94e-15]*I" ] }, "execution_count": 52, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ini = list(cbt.egf[:6])\n", "%time c = (dop.numerical_transition_matrix([0,s])*vector(ini))[2]\n", "c" ] }, { "cell_type": "code", "execution_count": 53, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "[1.6248570260793 +/- 3.26e-14] + [+/- 9.96e-15]*I" ] }, "execution_count": 53, "metadata": {}, "output_type": "execute_result" } ], "source": [ "expo = dop.local_basis_monomials(s)[2].op[1]\n", "CBF(-s)^expo*c/CBF(-expo).gamma()" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### A conjecture of M. Kontsevich\n", "(via D. van Straten and A. Bostan)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* $L = D_x\\,x\\,(x-1)\\,(x-t)\\,D_x + x$ for $t > 0$ (small)\n", "* $M(t)$ = transition matrix along a simple loop around $\\{0, t\\}$\n", "* $λ(t), \\barλ(t)$ = eigenvalues of $M(t)$\n", "\n", "Then\n", "\n", "$$t \\mapsto \\left(\\frac{\\log(λ(t))}{2πi}\\right)^2$$\n", "\n", "is analytic at $0$, and its Taylor coefficients are rationals with small denominators.\n", "\n", "**Task:** Compute that series (heuristically!)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "### ⚐ Numerics: Decent high precision evaluation at singularities" ] }, { "cell_type": "code", "execution_count": 54, "metadata": { "scrolled": true, "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/plain": [ "(494, 9960)" ] }, "execution_count": 54, "metadata": {}, "output_type": "execute_result" } ], "source": [ "terms = 20\n", "sz = ceil((terms + 1)^2 * sqrt(ZZ(terms + 1).nbits())/2)\n", "prec = terms*(sz + 4)\n", "hprec = prec + 100\n", "C = ComplexField(hprec)\n", "sz, prec" ] }, { "cell_type": "code", "execution_count": 55, "metadata": { "scrolled": true, "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/plain": [ "(x^3 - 51146728248377216718956089012931236753385031969422887335676427626502090568823039920051095192592252455482604439493126109519019633529459266458258243585/51146728248377216718956089012931236753385031969422887335676427626502090568823039920051095192592252455482604439493126109519019633529459266458258243584*x^2 + 1/51146728248377216718956089012931236753385031969422887335676427626502090568823039920051095192592252455482604439493126109519019633529459266458258243584*x)*Dx^2 + (3*x^2 - 51146728248377216718956089012931236753385031969422887335676427626502090568823039920051095192592252455482604439493126109519019633529459266458258243585/25573364124188608359478044506465618376692515984711443667838213813251045284411519960025547596296126227741302219746563054759509816764729633229129121792*x + 1/51146728248377216718956089012931236753385031969422887335676427626502090568823039920051095192592252455482604439493126109519019633529459266458258243584)*Dx + x" ] }, "execution_count": 55, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Pol. = QQ[]\n", "Dop. = OreAlgebra(Pol)\n", "t0 = 2^(-sz)\n", "L = Dx * (x*(x-1)*(x-t0)) * Dx + x\n", "L" ] }, { "cell_type": "code", "execution_count": 56, "metadata": { "scrolled": true, "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "m1 = L.numerical_transition_matrix([0, t0/2], 2^(-prec)).change_ring(C)\n", "m2 = L.numerical_transition_matrix([t0, t0/2], 2^(-prec)).change_ring(C)\n", "delta = matrix(C, [[1, 0], [2*pi*I, 1]])\n", "mat = m1*delta*~m1*m2*delta*~m2" ] }, { "cell_type": "code", "execution_count": 57, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "9.556619453472961320858516937971887816879016410846290172126354839510701168391942496813468949482302488162252696886520244471966592218665291723177173168199392132406276973789705340073690947932553910124796231579916084856421141195367398954924504111462358438957864101423930407260331976176725065395140834421424427590908128799409817584336622391673407243198108939261965265090478623967959224212373223620909718920076884504325493512957093231290793326009149024368397362080608773817261344914045941037215573589929466261952514844246399651185700946374898667516902390725364103154432169875542859847614933659911701076729537808003318825249648022560871556004599073286958174424349446361789734983296900487540550803091366045082974239981262926261619181189568107287215173880557039860815668766673099266893087110588547841898487650663644043980236519485686759761011681600245421501121353018152348931578768245258254780107760973263563775281863663012073352722917223246017584235308522688246510732993890785886249250725937793549943881255036341563534916875647477037886097657249792008428528870297161215131300381890686100703684816317908759813352398646880833758301441006787007561756610067089326048633636601655581680294114920232193088316072162202643820964032906843841817185410145450874756553779319769681938557795157701492901383950831572956075682891800333143265314437926822965689897336891281115810393634163501579557476432032208140913208912647082900790120127010816341829130398447859634549497723253118914899666245483559245438168221086949951270207323864052877486421440537213113570849672571099832804454703454119688641196701003291562809385222820524421206480989940681005905092328271286077626602437533505092641703146035760742921680842263781843854873991144874044147403548564686612183559041672386026367636029777757668704070453267111227724964978924932012765810717992589937434808382663799797773043364171093373260596575966020347101598724092360916823714857848399257003615082945323894832711059841119689648486033169179927891587569070063188570258161818486202606231937307760766826338422121070834991233113786550649092555451055378950142673684858167685590606256915969226571204633960531410689676714870414271577763534727803754179156714701688996696004086261971352341976630738528832473896535081911494009668569113201638788353744280283060834315636625364415819665437542451448080606771694562471371975196008988383248568008014743171957640899885078174330701068772917303994305295373221848857890630191210892023261909597212215222580922194025813235619871426751892331277173873023952664284270669470966003915290659729911588342236100212184270452170452531046709945743782625626814519771418724775958330968324506098946055581698099590764388244028515969168881038955796413001853465395680656243512452699320906509054087659901425563057011885470050958906345477915933201856109902467292184386952372661834542279554976803521112299906835150188799637765061595534103641294316535688856568007252802956390958368379318710244523774844628336072323710306061629273180267967815264454788464436467987614675704690957518274321357948708194043185e-299" ] }, "execution_count": 57, "metadata": {}, "output_type": "execute_result" } ], "source": [ "tr = mat.trace().real()\n", "Pol. = C[]\n", "char = (la+1/la-tr).numerator()\n", "rt = char.roots(multiplicities=False)[0]\n", "a = (rt.log()/C(2*I*pi))^2\n", "a" ] }, { "cell_type": "code", "execution_count": 58, "metadata": { "scrolled": true, "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "data": { "text/plain": [ "[0,\n", " 0,\n", " 1/4,\n", " 1/24,\n", " 101/576,\n", " 239/17280,\n", " 19153/115200,\n", " -1516283/72576000,\n", " 23167560743/121927680000,\n", " -5350452180523/76814438400000,\n", " 8122785754979827/32262064128000000,\n", " -10922037427834714189/74525368135680000000,\n", " 257615133666208067057417/688614401573683200000000,\n", " -5719273111411836892974100997/20679090479257706496000000000,\n", " 749577901737131766662423170141529/1241986174184217852149760000000000,\n", " -56589264915861458106843233145366159317/111890534432256186300171878400000000000,\n", " 166060654964554378941594941573970938041843381/161283491952031357180519752400896000000000000,\n", " -227464761949125949651795174312351837947239351830471/247010506429294584462681416394544250880000000000000,\n", " 1379581218345710272912991297578804281398062616806207006247/756608001823315069884260939301472713100492800000000000000,\n", " -37154931571287997286276581805653700634917295203190508379056930613/22016589207616772850617000970999305581601157021696000000000000000]" ] }, "execution_count": 58, "metadata": {}, "output_type": "execute_result" } ], "source": [ "coeffs = []; den_ratios = []; cur = a\n", "for k in range(terms):\n", " rat = QQ(pari.bestappr(cur, 2^(sz/2+10)))\n", " cur = (cur - rat)/t0\n", " if k >= 2:\n", " den_ratios.append(rat.denom()/coeffs[-1].denom())\n", " coeffs.append(rat)\n", "coeffs" ] }, { "cell_type": "code", "execution_count": 59, "metadata": { "scrolled": true, "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "data": { "text/plain": [ "[4,\n", " 6,\n", " 24,\n", " 30,\n", " 20/3,\n", " 630,\n", " 1680,\n", " 630,\n", " 420,\n", " 2310,\n", " 9240,\n", " 30030,\n", " 60060,\n", " 90090,\n", " 1441440,\n", " 1531530,\n", " 3063060,\n", " 29099070]" ] }, "execution_count": 59, "metadata": {}, "output_type": "execute_result" } ], "source": [ "den_ratios" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Desingularization" ] }, { "cell_type": "code", "execution_count": 60, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "(t^14 - 280163/1250*t^13 - 77823473621/3750000*t^12 + 15467808452531/58593750*t^11 - 25196309686447913/19531250000*t^10 - 668120393063773277/915527343750*t^9 + 1960569648056681994201/61035156250000*t^8 - 4671428124563366603059/30517578125000*t^7 + 41197718959797815058559/95367431640625*t^6 - 164817439510375233136933/190734863281250*t^5 + 5927091704553171594142627/4768371582031250*t^4 - 8659852913093432184666047/7152557373046875*t^3 + 5049362927700459914731448/7152557373046875*t^2 - 685309305030323264047472/11920928955078125*t - 18187783435081723962509824/35762786865234375)*Dt^5 + (29/2*t^13 - 8613043/2500*t^12 - 174431849123/500000*t^11 + 1255822655957817/312500000*t^10 - 20280225824735266/1220703125*t^9 - 49190254079706055749/976562500000*t^8 + 78414082510831851698897/122070312500000*t^7 - 18568806332381461314869/7629394531250*t^6 + 4003810338808418269811661/762939453125000*t^5 - 27771814401571351451778/3814697265625*t^4 + 11867297379816696715186737/1907348632812500*t^3 - 307053942520796638296658/476837158203125*t^2 - 1871856137277813387262933/476837158203125*t + 6163568693739391139347276/2384185791015625)*Dt^4 + (475/8*t^12 - 30050227/2000*t^11 - 16988836335053/10000000*t^10 + 21958673179999589/1250000000*t^9 - 9131560294338258463/156250000000*t^8 - 401564615212972941641/976562500000*t^7 + 104294200346934559073119/30517578125000*t^6 - 637131873161718618633403/61035156250000*t^5 + 4868593843047491215960457/305175781250000*t^4 - 426971833579701247214697/38146972656250*t^3 - 775947208574144990906417/953674316406250*t^2 + 9713861736043758399614103/953674316406250*t - 3186441544992332414987113/476837158203125)*Dt^3 + (1185/16*t^11 - 80462511/4000*t^10 - 52523948963397/20000000*t^9 + 120582837906110277/5000000000*t^8 - 18344488178046465327/312500000000*t^7 - 6879971713312897388329/7812500000000*t^6 + 2687398835047244023617471/488281250000000*t^5 - 130498326095692372331277/9765625000000*t^4 + 1786945003773501791092069/152587890625000*t^3 + 2499771123849829666054977/610351562500000*t^2 - 100956317858140812799725681/7629394531250000*t + 8521916368516547564271043/953674316406250)*Dt^2 + (5145/256*t^10 - 380739891/64000*t^9 - 305932960216191/320000000*t^8 + 153851731551379239/20000000000*t^7 - 6655819554769778547/625000000000*t^6 - 6447339574295055155727/15625000000000*t^5 + 157321670664515867630799/78125000000000*t^4 - 145611308765873254822287/39062500000000*t^3 - 29971841950167023004357/122070312500000*t^2 + 3523442028617667203992833/610351562500000*t - 14063874730153571703882243/3814697265625000)*Dt + 105/512*t^9 - 8901123/128000*t^8 - 11004229733343/640000000*t^7 + 4759356296189817/40000000000*t^6 + 11832518636433327/625000000000*t^5 - 57291975299178874233/6250000000000*t^4 + 214707950483595841419/6250000000000*t^3 - 3411921799561771135773/78125000000000*t^2 - 3177202651293055513191/19531250000000*t + 111251064741292878982587/610351562500000" ] }, "execution_count": 60, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from ore_algebra.examples import periods\n", "dop = periods.allODEs[2][0].numerator()\n", "dop" ] }, { "cell_type": "code", "execution_count": 61, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "[(-79.81779428070589?, 1),\n", " (-5.251936447658327?, 1),\n", " (-0.5142673948169302?, 1),\n", " (2.501571821316740?, 1),\n", " (2.754410642927720?, 1),\n", " (292.1291388587215?, 1),\n", " (0.1227946785361402? - 1.519364185445195?*I, 1),\n", " (0.1227946785361402? + 1.519364185445195?*I, 1),\n", " (0.2812140893416301? - 1.755844000088624?*I, 1),\n", " (0.2812140893416301? + 1.755844000088624?*I, 1),\n", " (1.749115902402969? - 1.266768315130795?*I, 1),\n", " (1.749115902402969? + 1.266768315130795?*I, 1),\n", " (4.011513729826854? - 4.25537902618774?*I, 1),\n", " (4.011513729826854? + 4.25537902618774?*I, 1)]" ] }, "execution_count": 61, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dop.leading_coefficient().roots(QQbar)" ] }, { "cell_type": "code", "execution_count": 62, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[(2.501571821316740?, 1),\n", " (2.754410642927720?, 1),\n", " (0.1227946785361402? - 1.519364185445195?*I, 1),\n", " (0.1227946785361402? + 1.519364185445195?*I, 1),\n", " (0.2812140893416301? - 1.755844000088624?*I, 1),\n", " (0.2812140893416301? + 1.755844000088624?*I, 1)]" ] }, "execution_count": 62, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dop1 = dop.desingularize().numerator()\n", "dop1.leading_coefficient().roots(QQbar)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Plots" ] }, { "cell_type": "code", "execution_count": 63, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAHWCAYAAAD6oMSKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xd4VFX+x/HPJEASIESKBCkr7IpiFpFiQFSUXRFXUaQpNlAsawFX7KKufUFcdN2VoIKCroii9HVdhZ8KqCgiRRAVbEiNNE1CS53fH8ebIdSUmTm3vF/PM8/3ZkhmPqR+59x7zgmFw+GwAAAA4HkJtgMAAAAgOmjsAAAAfILGDgAAwCdo7AAAAHyCxg4AAMAnaOwAAAB8gsYOAADAJ2jsAAAAfILGDgAAwCdo7AAAAHyCxg4AAMAnaOwA+F44HFZubq7YGhuA39HYAfC9vLw8paWlKS8vz3YUAIgpGjsAAACfoLEDAADwCRo7AL6VlZWljIwMZWZm2o4CAHERCnM1MQCfy83NVVpamnJyclSnTh3bcQAgZhixAwAA8AkaOwAAAJ+gsQMAAPAJGjsgIHbulKZNk668UqpXTwqFzC01VerVS8rKknbtsp0SAFAVTJ4AfO6tt6QePSr2MW+8IfXrF5s8NjB5AkBQMGIH+NTatWZEbu+m7qGHpKVLpZISKRyO3H78UfrXvyLvd+GF5mNHjox/bgBA5TFiB/hMcbF08snSZ5+Zt5OTpY0bpbp1y/fxO3dKvXtLc+ZE7lu7VmrWLPpZ44UROwBBwYgd4CM7d0rVqkWauoULpd27y9/USVKtWtLs2VJRkXTiiea+3/xGGjIk+nkBANFFYwf4xKZNUu3a5rhPH3OKtWPHyj9eYqK0bJn0wQfm7awsc3p2z56qZwUAxAaNHeADK1ZIjRub4/vuk6ZOjd5jn3aauSavZUvzdkqKtH599B4fABA9NHaAx332mdSmjTl+8UXpkUei/xyhkLR6tfTAA+btZs2kL76I/vMAAKqGxg7wsM2bJWd/+3ffla64IrbP9+CD0ssvm+MTTpDmz4/t8wEAKobGDvCocFhKTzfHkydLf/xjfJ738sul//s/c3zGGeY6PACAO9DYAR516qmmnnuudNFF8X3uM8+U3n/fHLdrJ33/fXyfHwBwYDR2gAe9/LL08cfm+L//tZOha1dpyhRz/LvfmdPCAAC7aOwAj9m4URo40Bxv3243S9++0pgx5jg9nb1mAcA2GjvAQ8JhqUkTc/zWWxVbeDhWbrhBuv9+c1yrlt0sABB0NHaAh1x8sal9+kjnnGM3y94eekg66SRz3KmT3SwAEGTsFQt4xIYNUtOm5titP7WhkKl//7t0++12s+yNvWIBBAUjdoBHOE2dm9eO27HD1DvukD75xG4WAAgiGjvAA/79b1Pr1ZO6dLGb5VBq1TLbm0lS587295XNyspSRkaGMp1VnAHA5zgVC7hcSYmUmGiOd+zwxgSFv/9duvNOc+yG3zCcigUQFIzYAS7nTJi47z5vNHWSORXrZL3nHrtZACBIGLEDXGznTql2bXPstZ/UwkKpRg1z/PXX0nHH2cvCiB2AoGDEDnAxZwmR116zm6MyqleXFi82x61a2c0CAEFBYwe4VE6OGemSpP797WaprPbtpZ49zXG/fnazAEAQ0NgBLnXKKab+5z92c1TVjBmmTp0qLV9uNwsA+B3X2AEulJsrpaWZYz/8hH71lZSRYY5LSiILGccL19gBCApG7AAXOv10U6dPt5sjWo4/XhowwByffbbdLADgZ4zYAS7j5Zmwh+OM1C1dKrVtG7/nZcQOQFAwYge4zEUXmfrqq3ZzxMKqVaa2a2c3BwD4FY0d4CLhsPTWW+bYWZjYT449Vjr/fHM8ZIjdLADgRzR2gIsMH27q0KF2c8TSzJmmZmVJ27fbzQIAfsM1doCLONegFRdLCT5+2fXmm5GRu3j8BuIaOwBB4eM/HYC3zJ5taseO/m7qJOm88yLHXl+nDwDchBE7wCWc0bpt26R69exmiYft26X69c1xrNe2Y8QOQFD4fFwA8IaffoocB6Gpk8z/05lA0aeP3SwA4Bc0doALOPupzptnN0e8Pf20qTNmlG1uAQCVQ2MHuMCnn5rq7DgRJHPnmtqokdUYAOALNHaAZW+8Yepll9nNYcsZZ0SO333XXg4A8AMmTwCWOZMG8vIiW4kFzZYtUsOG5jgWv5GYPAEgKBixAyzKz48cB7Wpk6Qjj5R69DDHf/2r3SwA4GU0doBFd9xh6rhxdnO4waxZpj76aNmGtyqysrKUkZGhzMzM6DwgALgcp2IBi5zTsLFex80rnn1WuuEGqVMn6ZNPove4nIoFEBSM2AGWfP+9qU2a0NQ5rr/e1IULWf4EACqDxg6wpFcvU996y24Ot3GWPznqKKsxAMCTOBULWOKM0vETuD/nc7N4sdS+fdUfj1OxAIKCETvAgqlTTQ3q2nWH88MPpnboYDcHAHgNjR1gQb9+pj77rN0cbtW8ufTb35rjyZOtRgEAT+FULBBnhYVSjRrmmJ++g9uxQ0pNNcdV/TxxKhZAUDBiB8TZyJGmPvGE3RxuV7u21L+/Ob7/frtZAMArGLED4syZGFBUJCUm2s3idiUlkc9RcbGUUMmXoozYAQgKRuyAOCooiBzT1B1eQoJ0333m+Oqr7WYBAC+gsQPiaPhwU//5T7s5vOSRR0x98cWyjTEAYH80dkAcPfSQqUOG2M3hNWPGmHrOOXZzAIDbcY0dECf5+VJysjnmp67inGsTc3Mjs2XLi2vsAAQFI3ZAnDinFJ95xm4Or3IWdT7pJLs5AMDNGLED4sQZcSopiRyjYpzPW3a2lJ5e/o9jxA5AUDBiB8RBSUnkmKau8ubNM9XZlQIAUBaNHRAHL75o6j33WI3heaefbuquXdI339jNAgBuRGMHxIGzBtu999rN4RZjxoxRixYtlJycrA4dOuiDDz4o18e99tprkkKSeunYY2MaEQA8icYOiKOaNW0nsG/y5MkaOnSo7r33Xi1dulRdunTROeeco7Vr1x7y43788Ufdfvvt6tKlS+l9S5bEOi0AeAuTJ4AY++gj6bTTpN69pWnTbKexr1OnTmrfvr2e2Wt68PHHH69evXppxIgRB/yY4uJinXHGGRo0aJA++OADbdz4i+bMmSGpfEvHMHkCQFAwYgfE2DXXmPrUU3ZzuEFBQYEWL16s7t27l7m/e/fuWrBgwUE/7uGHH9aRRx6pq389p12zplS7tvm3Dz+MWVwA8JxqtgMAfvf116b+5jd2c7jB1q1bVVxcrPR91ipJT09Xdnb2AT/mo48+0gsvvKBly5aVuf/bb6VGjaQuXfYftcvPz1d+fn7p27m5udH5DwCAyzFiB8TQunWmtmplN4fbhPZZ8yUcDu93nyTl5eXp8ssv17hx49SgQYMy/5aeLjVtao7nzCn7cSNGjFBaWlrprVmzZlHNDwBuxYgdEEM332zquHF2c7hFgwYNlJiYuN/o3ObNm/cbxZOk7777TmvWrNH5559fel/Jr4sCVqtWTQsXrtJJJ/1O3buXHbUbNmyYbr311tK3c3Nzae4ABAIjdkAMTZ9u6mmn2c3hFjVq1FCHDh00Z58htjlz5uiUU07Z7/1btWqlFStWaNmyZaW3nj176g9/+IOWLVumE05oVjoaOmtW5OOSkpJUp06dMjcACAJG7IAY2bXLdgJ3uvXWWzVgwACddNJJ6ty5s8aOHau1a9fq+uuvlyQNHDhQTZo00YgRI5ScnKzWrVuX+fgjjjhCkkrvX7hQSkuTLrigfDNkAcDPaOyAGBk50tTnnrObw2369++vbdu26eGHH9amTZvUunVrvfXWWzr66KMlSWvXrlVCQvlPJtSpI7VrJy1dKk2ZIvXrF6vkAOB+rGMHxIgzF6CoSEpMtJvF73bskFJTzfGBfqOxjh2AoOAaOyDGaOpir3ZtqVMnczxpkt0sAGCTKxq7nBzp8celPn2kv/1NKiy0nQiomo8/NrVPH7s5guS990y97DK7OQDAJqunYouLpbPOkt5/f/9/O/54aeXKyOkswEs6dZI+/dQsovu739lOExxdu0rz5knjx0uDBkXu51QsgKCw1tjt2SOlpETeHjNG+tOfpJkzpVtuidy/ZYu0z7qkgOs5L0i4gjW+9v69svfnnsYOQFBYOxXr/PK9+GLzC/iGG6QWLaShQ83bF11k/v3II83IHuAVzjInycl2cwRRcrLUrZs5HjvWbhYAsMFKY9eli6nHHy+9+uqB32fy5MiplGosygIP+cc/ylbE15tvmnrddXZzAIANcW/sJkyQPvzQHK9ceej3HT9eql/fHD/wQGxzAdFy332mXn213RxBlZRkLuuQpGeftZsFAOItrtfYhcOSs+5oXp5ZouBwiosjI3bl/RjAJq6vs6+gwDR4kvk6cI0dgKCI64hd//6m3nFH+Ru0xERp4kRzzB7ecLvvvjP1ANueIo5q1JDOPdccjxljNwsAxFPcRux++UWqW9ccV+YZnVGQTz+VMjOjlwuIposvNteHfvihdOqpttME296jdjk5jNgBCIa4jdidcIKp//lP5T5+3TpTO3aMTh4gFiZPNpWmzr4aNaSePc0x19oBCIq4jNgdbG2pinJG7Vatko49tuq5gGja+xpSrq9zh8JC0+BJuZIYsQPgf3EZsRswwNQJE6r2OF9+aepxx1XtcYBYmD3b1GuvtZsDEdWrS717204BAPETlxG7aM4SdB5r40bpqKOq/nhAtDjbiK1dy0QfNykqkqpXZ8QOQDDEfMTuqadMveaa6DzeRx+Z6sx4A9zi009Npalzl2rVpAsuMMdZWXazAECsxXzEbuRI6e67zbUu0dpBwhm1KymJHAM2cX2du23fnqv69dMk5SgcZsQOgH/FfMTurrvMH7pobgt2ww2mPvdc9B4TqApnG6u//MVuDpSVlZWljIwMde4cWSNp1CiLgQAgxuK680S07LuqPGBbu3bSsmVc++lWzs4TUo6kOvzeAOBbcd8rNhrM8gVGTo69HIBj2TJTaerc7cILTX3sMbs5ACBWPNnYSdK0aaZy6gu2lZTYToDyci7fGDbMbg4AiBXPNnbO2lT//rfdHMD06abefrvdHDi8xMTIupqPPmo3CwDEgievsXOkp0ubN0vr10tNmthOg6Bq3VpauVL66SepYUPbaXAgzjV2OTk5qlWrTulkLu/+9gOAA/PsiJ0UGSkZNMhuDgTbypWm0tR5Q2KidOWV5vihh6xGAYCo8/SInRTdXS2AiiouFqM/HrD3iF2dOnVUUmIaPImvGwB/8fSInSQ1bWrqTz/ZzYFgeuMNU++5x24OVExCgnT11eb4r3+1mwUAosnzI3azZ0tnny1dfrn08su20yBojjtOWr1a2rZNqlfPdhoczL4jdpLKjNqxiw0Av/D8iF337qZOnGg3B4Jp9WpTaeq8JyFBuvFGc3zXXXazAEC0eH7EToq80t6zJ7IjBRBrRUVS9erm2Ps/Rf52oBE7qewev4zaAfADz4/YSdLo0aY+/rjdHAiWSZNMfeABuzlQeaFQZP1BFjsH4Ae+GLFjZiJsaNFCWrNG+vln6YgjbKfBoRxsxE4qO2pXXBw5BgAv8sWvMOcCaCCe1qwxlabO20KhyMzYa6+1mwUAqsoXjZ0k9elj6qef2s2BYCgutp0A0fTww6aOH8/XFoC3+aaxe+wxU2+7zW4OBMOUKabed5/dHIie4cNNdfaSBQAv8sU1dg52oUC8nHiitHy5tGWL1KCB7TQ4nENdY7c353dIYWHkul0A8BLfjNjtraDAdgL43fLlptLU+ctTT5nar5/dHABQWb5q7EaNMvXZZ+3mgL8xIuxfN99s6syZUn6+3SwAUBm+auwGDzaV654QS/Pnm+rsNQp/GTvW1B497OYAgMrw1TV2EtfZIfbOO0/673+lr782e8XC/cp7jZ3D+T2yY4dUq1aMwwFAFPlqxE6SjjrK1Lw8uzngX//9r6k0de6XlZWljIwMZWZmVujjnF1FunSJQSgAiCHfjdi98IJ0zTXS3/4m3XOP7TTwI0aFvaeiI3ZS5Ou8fbtUt24MwwFAFPluxG7gQFOdleSBaNq40dSTT7abA7HnjMyecILdHABQEb5r7KpXN7WkxG4O+NO//mXqHXfYzYHYO/dcUzdsiDT0AOB2vjsVK0ktW0rffsspFERfYqJ50VBUxB7FXlKZU7GS9PHH0imnmBeMrI8JwAt8N2InSffea+ozz9jNAf9xRoJp6oKhc2dTCwulb76xmwUAysOXI3YFBVJSkpSaKuXm2k4DvygslGrUkGrXZta111R2xE6Svvgicp2d/35bAvAbX47Y1ahhKn98EU2TJ5vK9XXB0rp15HjxYns5AKA8fDliJ0nJyWZLIDbzRrS0bSt9/rm0ZQt7xHpNVUbsJOnHH6Xmzc2xP39jAvALX47YSdLdd5v6+ut2c8A/Pv/cVJq64Dn66MjXffZsu1kA4FB8O2K3dat05JHm2pjly22ngR+wMLF3VXXETpK2bYs0d3wPAHAr347YOb+AV6ywmwP+sGyZqRdeaDcH7KlfX3J2JpswwW4WADgY3zZ2QDSNGmXq7bfbzQG75s839aqr7OYAgIPxdWM3YICpH39sNwe875VXTO3Y0W4O2JWcLF12mTm+7z67WQDgQHx7jZ0UWX+qZ09p5kzbaeBlXF/nbdG4xs5RUhJZoLq4WErw9ctjAF7j619JzvpTs2bZzQFvcxa5btnSbg64Q0KC9NBD5viSS+xmAYB9+bqxA6LhxRdNHTrUagy4yP33m/r669LOnXazAMDefN/YtWplak6O3RzwrqeeMnXgQLs54C7OdZft29vNAQB7831jd/PNprI8ASrrhx9MrV3bbg64y6WXmrp6tbR2rd0sAODw9eQJyZwmqV3bbAfk/IEGKoKJE94XzckTe/vss8jadnx/AHAD34/Y1apl6po1VmPAo5xtxFiYGAdy0kmR4w8+sJcDABy+b+yAqvjnP01l4gQOZtMmU08/3W4OAJAC0tj16WPqF1/YzQHvca7N7NzZbg5UTlZWljIyMpTpnC+NgUaNpA4dzPG4cTF7GgAoF99fYyeZbYDOOEO69lpp7FjbaeAlXF/nD7G6xs6Rn292pZDMAsbO9w0AxFsgRuy6dDGVV9OoiIICU+vVs5sD7peUJN1wgzm+8Ua7WQAEWyAaO149ozKmTTP1L3+xmwPekJVl6rPPSrt3280CILgC0djtjVNqKK/Ro0299lq7OeANoVBkl5K2ba1GARBggWnsrrvOVJYkQHl99JGpjRvbzQHvuOIKU1evlr77zm4WAMEUmMbu+utNfeYZuzkA+Juz9uExx9jNASCYAtPYOadGXnvNbg54w88/m3rCCXZzwHvatJGOOMIcz5xpNwuA4AlMYwdUhLN+3eDBdnPAm5ydbnr1shoDQAAFqrGrWdPU4mK7OeB+zgzHyy6zmwPelJYmDRhgjnlxACCeAtXYOetMvf223Rxwv++/N7V2bbs54F0vvWTqmDFSbq7dLACCI1CN3Z//bOrzz9vNAcD/QiFpyhRzzMxqAPESqMbu2GNNnTHDbg64mzNad9ZZdnPA+/r2NXXnTunDD+1mARAMgWrsgPJ49llThwyxmwP+sH69qc7WhgAQS4Ft7NiBAgfj7DjRo4fdHPCHJk2kc84xx8OG2c0CwP8C19hddZWpn35qNwfcy9nnMzHRbg74x5tvmvrYY+a0LADESmAbuxdesJsD7sRILmIhIUF6+WVz/Lvf2c0CwN8C19idcoqpNHY4kIULTR040G4O+M/ll5v600/SokV2swDwr8A1dqGQqSUldnPAnZyFiVlUNrbGjBmjFi1aKDk5WR06dNAHH3xw0PcdN26cunTporp166pu3brq1q2bPvXotRQ//GBqx452cwDwr8A1dsChTJxoKn94Y2fy5MkaOnSo7r33Xi1dulRdunTROeeco7Vr1x7w/efOnatLLrlE77//vj7++GP95je/Uffu3bVhw4Y4J6+65s2ls882x7fcYjUKAJ8KhcPBu6roggukWbOk1aulli1tp4GbOCO6wfupiJ9OnTqpffv2euaZZ0rvO/7449WrVy+NGDHisB9fXFysunXravTo0RpYznPmubm5SktLU05OjurUqVPp7NFQUhKZmLNtm1SvntU4AHwmkCN2zgSKF1+0GgMu4+whnJJiN4efFRQUaPHixerevXuZ+7t3764FCxaU6zF27dqlwsJC1fNoR5SQYF5YSlL9+nazAPCfQDZ2555r6vjxdnPAXWbPNvX66+3m8LOtW7equLhY6enpZe5PT09XdnZ2uR7j7rvvVpMmTdStW7eDvk9+fr5yc3PL3Nzk/PMjx2+8YS8HAP8JZGNXvbqp5fw7goB47jlTr7vObo4gCDnnvH8VDof3u+9AHn/8cb366quaNm2akpOTD/p+I0aMUFpaWumtWbNmVc4cbb/8YupFF0lFRXazAPCPQDZ2wIHMnGnqccfZzeFnDRo0UGJi4n6jc5s3b95vFG9fo0aN0vDhwzV79my1adPmkO87bNgw5eTklN7WrVtX5ezRlpYmPfSQOc7MtJsFgH8EtrFz1rPbssVuDiBIatSooQ4dOmjOnDll7p8zZ45OcX4oD+Dvf/+7HnnkEb399ts66aSTDvs8SUlJqlOnTpmbG91/v6nLlrEbDoDoCGxj50ygcJa3QLA5p8Lq1rWbIwhuvfVWPf/88xo/fry++uor3XLLLVq7dq2u//XixoEDB2rYXpuqPv7447rvvvs0fvx4NW/eXNnZ2crOztaOHTts/Rei6vvvTe3UidnYAKousI3dhRea6mzzg2Bz9vJk4kTs9e/fX0899ZQefvhhtW3bVvPnz9dbb72lo48+WpK0du1abdq0qfT9x4wZo4KCAvXr109HHXVU6W3UqFG2/gtR1aKFNGCAOe7Xz24WAN4XyHXsHKxZBsef/iS9847ZGaB5c9tpEG1uWsfuYJzfR198If3+93azAPAuGjvR2IHvBb/zQmO3YoXkzAnh+xBAZQX2VKwUWfHdWZgWAGw54QSpb19zzClZAJUV6MbO2Y1o7lyrMWBZYaGpDRvazQFMmWLq1KnS8uV2swDwJho7Sf/+t90csMuZOPHnP9vNAUjSypWmnngip2QBVFygG7u2bU2lsQu2sWNNvfpquzkAScrIkC6+2Bz37Gk3CwDvCfTkCYmL5sH3QBB4YfLEvpzvy0WLpHKsyQwAkgI+YgcAbrVqlamZmVJJid0sALwj8I3dmWeaumGD3Ryww9lxwpkhDbjFscdKQ4ea44wMu1kAeEfgG7vLLzd18mS7OWDH//5n6rXX2s0BHMg//mHqqlXStGl2swDwhsA3dr17m/rqq3ZzwA4mTsDttm41tW9f6Zdf7GYB4H6BnzwhcfF8kPG1DwYvTp7Y28svR5Zn4nsVwKEEfsQOANxuwACpVStzzOgygEOhsZNUrZqpvBIOFmcrudRUuzmA8nAWLh4/XvrkE7tZALgXjZ2kSy4x9fPP7eZAfM2ebeo119jNAZRHQoL0zTfmuHNnKT/fbh4A7kRjp0hjxwSKYHnhBVOZEQuvOOYY6W9/M8fJyXazAHAnJk/IbAJfo4bUtKm0bp3tNIgXJk4Eh9cnT+wrMdEsWnzffdIjj9hOA8BNGLGTVL26qevX280BILqysrKUkZGhzMxM21GiavduUx99VPrsM7tZALgLI3a/YvQmWMJhc82Scwx/89uInSR98YV0wgnmeOdOqWZNu3kAuAMjdr9q0sTUwkK7ORAfzqzCq66ymwOorNatpSefNMe1atnNAsA9aOx+5UygePdduzkQH87ECdYEg5fdcovZU1aS+vSxmwWAO9DY/YqZscHiNHadO9vNAVTV11+bOn26NGGC3SwA7OMau18511xVry4VFNhOg1jjmspg8eM1dnvbvl2qX98cL18eufYOQPAwYvcr5w8919gB8Jp69SLXjbZpI+Xl2c0DwB4aOwTOV1+Z2ru33RxANHXqJP3jH+a4Th1Go4GgorHbS/v2pubk2M2B2Bo/3lQmTsBvhg6V/vhHc+xMqgAQLDR2e3EmUMyYYTcHYsuZOHH22XZzALHgzOz/9lvpttvsZgEQfzR2e+nf39RJk+zmQGz9/LOp1arZzQHESn6+qU8+Kb38st0sAOKLWbH7YLak//E1Dh6/z4o9kJ9+kho1MscffyydfLLdPADigxE7BMqGDaZ27Wo1BhBz6enSkiXmuHNn9sIGgoLGDoHy0kumMnECQdCunfT66+a4WTNpxw67eQDEHo3dPrp1M3XjRrs5EBvPP29q3752cwDxcuGF0kMPmePUVKmoyG4eALFFY7ePiy4ydepUuzkQGz/8YGpKit0cQDzdf7902WXmuHp1ri8F/IzGbh/OorXO6QsA8IOJE6WTTjLHCfzmB3yLH+99NGhg6ocf2s2B6MvNNZV9NBFUixZJjRub4yOOsJsFQGzQ2CEwnFHYQYPs5gBscmaG5+RIrVrZzQIg+mjsEBgTJpjqXGsEBFVJiamrVpk9ZgH4B43dAZxzjqnr1tnNgehasMDUhg3t5gBsC4Wk4mJz/OmnLF4M+AmN3QE4W4tNmWI3BwDESkJCpLlbuFDKzLSbB0B00NgdwAUXmMrMWP9w/oDVrWs3B+IrKytLGRkZyqRrOaC9m7vPPpOOP95uHgBVx16xB8F+ov7y9tvmFPutt0pPPGE7DeItiHvFVkQ4HFkCpV49ads2u3kAVB4jdgiEF1809corbaYA3CkUikyo2L7dvM2LWsCbGLE7CEbs/IWvZ7AxYld+zZtLP/5ojvPzpRo1rMYBUEGM2B1Ez56mrlljNQYAxNWaNdIZZ5jjpCROywJeQ2N3EM6esW+8YTcHAMTb3LnS4MHmuEEDafVqq3EAVACN3UGcf76pzIz1vi++MNVp1gEc3ujRUlaWOT7uOOm//7WbB0D50NgdhHMZzmef2c2BqnMmTrCVGFAxN94ozZljjs87T3rwQatxAJQDkycOgQvu/aFBA3OdUFGRlJhoOw1sYPJE1axZI7VoYY47dZI++cRqHABlKvxDAAAgAElEQVSHwIgdfM+5+JumDqic5s2lXbvM8cKF5kXvzp1WIwE4CBq7Q3D2jF2/3m4OALAtJcWcvejY0bxdu7a0bJndTAD2R2N3CP36mTptmt0cqLwtW0w95RS7OQC/WLhQGjXKHLdrJz3wgN08AMqisTsEZ8/YKVPs5kDlTZpkKhMngOi57bbIaN3DD0s1a5prWAHYx+SJw2AChbe1bSt9/rmUkxOZ6YzgYfJEbBQUmEWMHUuXmp85APYwYgdf+/xzU/lbDkRfjRrmRe8dd5i327WTrrvObiYg6BixOwxG7LyNrx8kRuziYeVKqXXryNsbNkiNG9vLAwQVI3aH8Yc/mLp5s90cqLiCAlObNbObAwiC3/9eKi6W2rQxbzdpIt19t91MQBDR2B2GMzN2xgy7OVBxb71l6hVX2M0BBEVCgrn84f/+z7w9cqQZNWfJKCB+aOwOo3dvU5kZ6z0vvWTqwIF2cwBBc+aZZpbsiSeat5s149o7IF64xq4cuE7Lm/i6wcE1dvbMmyd17Rp5+6OPWFcSiCVG7AAAMXPGGVJJidSnj3n71FPNrhU//2w3F+BXNHYAgJgKhaSpUyPX2u3cKdWrJ/XsaZo+ANFDY1cOnTubyitM71ixwtRLLrGbA0BEkybm0oj//Me8/Z//SImJ0hNP2M0F+AmNXTk4M2OdX0ZwP2fixJVXWo0By7KyspSRkaHMzEzbUbCX884zDd7995u3b7/djOq98ordXIAfMHmiHH78UWreXDr/fGnWLNtpUB4NGkjbtpmZeYmJttPANiZPuFdRkWn03nknct/06VKvXvYyAV5GY1dOzLD0Fr5e2BuNnfvt2mUmWnz2WeS+F16QrrrKXibAizgVCwCwrmZNadEiKScnsv7d1VebF2kPPcSLNKC8aOzgO9u3m9qxo90cACquTh1p2TJp925z+YskPfig2dXi7LOlHTusxgNcj8aunJz9D3futJsDhzd5sqlsJQZ4V3Kyuaa5uFgaOtTcN3u2lJpqRvHmzrUaD3AtGrty6tvX1LfftpsDh+fMiL34Yrs5AFRdQoL0j3+YU7FTp0bu/8MfTIM3aJC0Z4+9fIDb0NiVk9PY7f2LBe60cKGp9erZzQEguvr0MQ3e5s1mBwtJevFFKSXFNHkTJ3ItHsCs2HIKh80rxxo1pPx822lwKMyIxb6YFetfEybsP3M2JUWaOVM66yw7mQCbGLErJ6dZKCiwmwOHVlxsKqN1QDAMGmRexO3cKQ0ebO7bvVvq3t383q5ZU3r1VV7oITho7OAr771nKjtOAMFSs6Y0erRp4LZsiVw+s3u3dOml5oxLKGRm2O7ebTUqEFM0dhXQooWphYV2c+DgnIkTzIgFgqtBA2nKFNPk5eREZtVKZk28mjVNk9e6tfTmm4zmwV9o7CrAeQX4/vt2c+DgnL0mneVpAARbnTqRWbVFRdIzz0jVq5t/W7nSrJXnjOb94Q/Su+/S6MHbaOwqgJmxAOBdiYnS9deba6XDYWnTJmnIkMi/z50rdesWafQaN5Yee8zsOw14BbNiK6CkxPxiqF9f2rrVdhocCDNicSDMikV5fP+9aeTGjTv4+7RvL11zjdS/P5O04E6M2FVAwq+fLV69udN335l63nl2cwDwpt/+Vho71rwwDIfNJIsJE0wz51iyRLrxRvMCPxQyt6Qk6fLLzRIru3bZyw9INHbwkZdfNpWJEwCiITnZzLBfvDjS7G3fbmbfduoUeb+CAnN9b69eUq1akYYvFJIyM6W//lX64APWQEV8cCq2gho0MCN2xcWRETy4Q4sW0po1ZnuhpCTbaeAmnIpFLIXD0vLl0rRp0vTp0ooV5fu4hATplFOkzp3NLTNTatIkckkJUBk0dhV0/fXSc89JCxaYH0S4B9fX4WBo7GDLjh3SvHnSnDmmLltW8ceoVUvKyCh7a9VKat5cqlYt6pHhcTR2FTRnjlnR/LbbpFGjbKfB3mjscDA0dnCrvDxp0SLp449NXbRI2rgxOo9dp47UtGnk1qSJuTVuLDVqZG7p6WarTPgHjV0FFRaaH4LmzaUffrCdBo5du8yr2mOPlVatsp0GbkNjB68rKJC+/Vb68kuz/t7XX5u3v/nGLMJsQ1KSlJZmbqmpUu3akVvNmuZ3snOrWdPcUlIit+Rkc0tKimwHmZhoRiETEw99S0jYvx7q5lzzGASVGsQNh8PKy8uLdhZPWbNGys21nQKOyZNNvfhivi6Q8vPzlb/XlerO76tcvjngYc7IW/fu5f+YcNj8Tty4MXLLzpY2bzZ1y5bIraINYn6+eZzNmyv2cUH2xhsV+/pJUmpqqkIV6EorNWLnvPoFAABA7FT0TEOlGruKjthlZmZq0aJFFX2aKonlc9avb7am2fvVTW5urpo1a6Z169bF9VSP3z63lX1O53VGtE5J8PX09nPuO2K3adMmdezYUV9++aWaNGkSs+fdlx8/t254Tn4+/fWcfD0PraIjdpU6FRsKhSr0yU9MTIz7dS2xfM5+/aTXXpPWrzezk/ZWp06duP5f/fa5repzRjsWX0//PKdkfkHy9fTHc0r8fPrpOSW+ntESl5XYBg8eHI+nidtzumnPWL99bivznM6Ysx/WFXTb59ZPz2lDUD63fD15Ti/y6+eWWbGVsHOnmfVz4omRNYmYdWfPp5+aVeAHDZLGj4/OY/L19Jf169eXnupp2rSp7TioIn4+/YWvZ3T5YIwj/mrVMvXzzyP3JSUl6YEHHlASWx7E3UsvmRrNrcT4evqL83Xk6+kP/Hz6C1/P6GLErpJYDNc9atc2o6hs84aDYUQAQFDwZxCet3OnqTR1AICg409hJZ17rqnr1tnNAQAA4KCxq6Q+fUydPt1ujqBzVjw/7TS7OQAAcAMau0q64AJT3bDkSZBNmmRqNCdOAADgVTR2ldSgganz5+//b9ddd51CoZCeeuqp+IYKIGdG7IUXRu8xCwsLddddd+mEE05QrVq11LhxYw0cOFAbN26M3pMAqLARI0YoMzNTqampatiwoXr16qVVq1bZjoUoGDFihEKhkIYOHWo7iufR2EXZjBkztHDhQjVu3Nh2lEBw1hGM5tbFu3bt0pIlS/TXv/5VS5Ys0bRp07R69Wr17Nkzek8CoMLmzZunwYMH65NPPtGcOXNUVFSk7t27a6czgwqetGjRIo0dO1Zt2rSxHcUXKrWlGA5sw4YNGjJkiN555x316NHDdhxUUlpamubMmVPmvqefflodO3bU2rVr9Zvf/MZSMiDY3n777TJvT5gwQQ0bNtTixYt1+umnW0qFqtixY4cuu+wyjRs3To8++qjtOL7AiF0VOBfsb9smlZSUaMCAAbrjjjv0+9//3m6wgCgsNDU9PfbPlZOTo1AopCOOOCL2TwagXHJyciRJ9erVs5wElTV48GD16NFD3bp1sx3FNxixq4K+faUPP5RmzpR++mmkqlWrpr/85S+2YwXGO++YGuuJE3v27NHdd9+tSy+9lMVtAZcIh8O69dZbddppp6l169a246ASXnvtNS1ZskSLFi2yHcVXaOwq4ZVXXtF1111XuuvE2LH/1Zo1/9SSJUsUcrakQMxFaysx5+vp+N///qcuXbpIMhMpLr74YpWUlGjMmDFVeyIAUTNkyBAtX75cH374oe0oqIR169bp5ptv1uzZs5WcnGw7jq+wpVgl5OXl6aeffpIktWwpSW8oFLpXCXttfVBcXKyEhAQ1a9ZMa9assZLT76K1rdveX09JatKkiVJSUlRYWKiLLrpI33//vd577z3Vr1+/ak8Ea9hSzF9uuukmzZgxQ/Pnz1eLFi1sx0ElzJgxQ71791ZiYmLpfcXFxQqFQkpISFB+fn6Zf0P50dhVkWkutmnFik1l7j/77LM1YMAADRo0SMcdd5yVbH4Xy/16nabum2++0fvvv68jjzwy+k+CuKGx84dwOKybbrpJ06dP19y5c9XSvLKGB+Xl5enHH38sc9+gQYPUqlUr3XXXXZxerwJOxUZFfbVuXXY0p3r16mrUqBFNnQcVFRWpX79+WrJkid58800VFxcrOztbkrlIu0aNGpYTAsE0ePBgTZo0STNnzlRqamrpz2VaWppSUlIsp0NFpKam7te81apVS/Xr16epqyJmxVZRu3am7thhN0fQfP21qb17R/+x169fr1mzZmn9+vVq27atjjrqqNLbggULov+EAMrlmWeeUU5Ojrp27Vrm53Ly5Mm2owGuwYhdFfXtKy1dKr31lnTRRZH7ua4utv79b1NjMSO2efPm4goFwH34ufS3uXPn2o7gC4zYVVHfvqayZ2x8vfiiqeeeazUGXC4rK0sZGRnKzMy0HQUA4oLJE1EQCkmJiVJRke0kwRHLiRPwHyZPAAgKRuyipLjYdgIAABB0NHbwnF93EdKJJ9rNAQCA29DYRYGzlFJ+vt0cQfHGG6ZeeaXVGAAAuA6NXRQ4EyjefddujqBwJk5ceqnVGAAAuA6NXRQ4jZ0zkoTY+ugjUxs2tJsDAAC3YVZsFITDUkKCVLu2lJdnO43/MSMWFcWsWABBwYhdFDiNBrtPxJ4z+5i/zQAA7I/GDp7y/vumxmLHCQAAvI7GLkqaNze1sNBqDN+bMMHUQYPs5gAAwI1o7KLkwgtNZWZsbE2aZGq7dnZzAADgRjR2UeI0dsyMBQAAtjArNkqYGRsfzIhFZTArFkBQMGIXJcyMjb2vvjK1Tx+7OQAAcCsaO3iGs+PEVVdZjQEAgGvR2EWRMzO2oMBqDN9yZsSefbbdHAAAuBWNXRRddJGpzIyNjS1bTK1WzW4OAADcisYuipgZCwAAbKKxi6IOHUylsYu+rVtNPflkuzkAAHAzGrsoYmZs7DgLE7PjBAAAB0djB09wJk707283BwAAbkZjF2XMjI2NZctMTUuzmwMAADejsYsyZsYC7pGVlaWMjAxlZmbajgIAccGWYlG2eLF00knSlVdGTh+iagoKpKQkqVkzae1a22ngRWwpBiAoGLGLsvbtTWVmbPS8+aapTJwAAODQaOyizJkZu3On3Rx+4ox8XnGF3RwAALgdp2JjwGnu+MxGB59PVBWnYgEEBSN2MdCihanMjAUAAPFEYxcDztZizIytOmeUrnp1uzkAAPACGrsYcJY8ef11uzn84MMPTb36ars5AADwAhq7GHBmxtLYVd0LL5hKYwcAwOExeSJGuOA/OpzPY0lJ5BioKCZPAAgKRuzgCTR1AAAcHo1djPz2t6bm59vNAQAAgoPGLkYuucTUd96xm8PLVqwwtX9/uzkAAPAKGrsYufhiU1991W4OLxs/3lQmTgAAUD5MnoghJlBUTZ06Ul6eVFwsJfASBFXA5AkAQcGfS7hWXp6pNHUAAJQPfzIBAAB8gsYuhjIzTf3lF7s5vGjtWlO7d7ebAwAAL6GxiyFnZuy0aXZzeNGECaZedZXdHAAAeAmTJ2JowwapaVPprLOk2bNtp/GWJk2kjRulPXukpCTbaeB1TJ4AEBQ0djHGzNjK4fOGaKKxAxAUnIoFAADwCRo7uM6GDaZ27Wo1BnwgKytLGRkZynRmMgGAz9HYxdg555jqzPLE4Tk7Tvz5z3ZzwPsGDx6sL7/8UosWLbIdBQDigsYuxpyZsZMn283hJWPHmtq7t90cAAB4DY1djF1wgansGVt+69ebmpxsNwcAAF5DYxdjzgS8pUvt5gAAAP5HYwdX2bjR1NNPt5sDAAAvorGLg1q1TGVNtsNzdpxg4gQAABVHYxcHAwaYumCB3Rxe4Eyc6NvXbg4AALyIxi4OLr/c1IkT7ebwAmdZGCZOAABQcWwpFgfhsJSQINWuLeXl2U7jbmwlhlhgSzEAQcGIXRw4zcqOHXZzuN2mTaaedprdHAAAeBWNHVzjxRdNZeIEAACVQ2MXJ85Wlb/8YjeHmzkTJ/r1s5sDAACvorGLE2cCxdSpdnO42Zo1pqakWI0BAIBn0djFSf/+pjIzFgAAxAqNXZykp5s6d67VGK7lTJxgxwkAACqPxg6uMG6cqTfeaDcHAABexjp2cZSQYNZnKymJLIECIz1d2rxZys+XatSwnQZ+wzp2AIKCEbs4ciZQLFtmN4cbbd5sKk2dv4XDYT344INq3LixUlJS1LVrV61cufKQHzNixAhlZmYqNTVVDRs2VK9evbRq1ao4JQYAb6Gxi6MrrzTVWa8NCJrHH39cTz75pEaPHq1FixapUaNGOuuss5R3iC1Z5s2bp8GDB+uTTz7RnDlzVFRUpO7du2vnzp1xTA4A3sCp2DgqKZESE6XUVCk313Ya9/jqKykjQ+rVS5o+3XYaxEo4HFbjxo01dOhQ3XXXXZKk/Px8paena+TIkbruuuvK9ThbtmxRw4YNNW/ePJ1eztk2nIoFEBSM2MVRwq+fbfaLLevZZ01l4oS//fDDD8rOzlb37t1L70tKStIZZ5yhBQsWlPtxcnJyJEn16tU76Pvk5+crNze3zA0AgoDGDtaNGWPqmWfazYHYys7OliSlO2v//Co9Pb303w4nHA7r1ltv1WmnnabWrVsf9P1GjBihtLS00luzZs0qHxwAPITGLs66djXVWbcNUlGRqQl8N/rKK6+8otq1a5feCgsLJUmhfaaEh8Ph/e47mCFDhmj58uV69dVXD/l+w4YNU05OTult3bp1lftPAIDHVLMdIGgGDTKLFE+cKN1xh+009nGFp3/17NlTnTp1Kn07Pz9fkhm5O+qoo0rv37x5836jeAdy0003adasWZo/f76aNm16yPdNSkpSUlJSJZMDgHcxRhJnffuaOmGC3Rxu8cknpl51ld0ciL7U1FQdc8wxpbeMjAw1atRIc+bMKX2fgoICzZs3T6eccspBHyccDmvIkCGaNm2a3nvvPbVo0SIe8QHAk2js4qxWLVO/+spuDrfIyjKViRP+FwqFNHToUA0fPlzTp0/XF198oSuvvFI1a9bUpZdeWvp+Z555pkaPHl369uDBgzVx4kRNmjRJqampys7OVnZ2tnbv3m3jvwEArsapWFj1yiumduhgNwfi484779Tu3bt144036ueff1anTp00e/Zspaamlr7Pd999p61bt5a+/cwzz0iSujoXqP5qwoQJutJZHBIAIIl17Kw49ljpm2+knTulmjVtp7HLuWae70LEEuvYAQgKTsVaMGiQqTNm2M1h26+TJHWI5cgAAEAF0NhZMHCgqUGfQDFtmqk33WQ3BwAAfsGpWEs4BSmddpr00UfSxo3SXqtfAFHHqVgAQcGIHaz56CNTaeoAAIgOGjtLUlJMLSmxmwMAAPgHjZ0lzgSK+fPt5rDF2VLt5JPt5gAAwE9o7Cy5+mpTX3jBbg5bnIWJb7vNbg4AAPyEyRMWBXkCRUqKtGePVFQkJSbaTgO/Y/IEgKBgxA5W7NljKk0dAADRQ2PnAkEbsSsuNjWB7z4AAKKKP60WXXutqc6yH0Hx5pumcn0dAADRRWNn0Z//bOpzz9nNEW9PPmnqX/5iNwcAAH7D5AmLwuHI6cggfRWCPGkEdjB5AkBQMGJnkdPgAAAARAONnUsEZfRq7VpTu3a1GgMAAF+isbPMuc5u3jy7OeLlX/8y9dZb7eYAAMCPuMbOshUrpDZtpL59pSlTbKeJPef0c3Exy50g9rKyspSVlaXi4mKtXr2aa+wA+B6NnQsEaTJBkP6vcA8mTwAICsZMEDdFRabWrGk3BwAAfkVj5wKNG5u6e7fdHLE2daqpXF8HAEBs0Ni5wE03mTp5st0csfbEE6YOHmw3BwAAfsU1di6wdat05JFS+/bS4sW208QO19fBFq6xAxAUjNi5QIMGpi5ZYjcHAADwNho7xMWnn5rav7/dHAAA+BmNnUuceqqpGzfazRErI0aYes89dnMAAOBnNHYu4UygGDvWbo5YmTHD1DZt7OYAAMDPmDzhEoWFUo0aUr160rZtttNEHxMnYBOTJwAEBSN2LlG9uqnbt9vNEQvffmvqmWfazQEAgN/R2CHmRo409d577eYAAMDvaOxc5PLLTV240G6OaHv+eVO7drUaAwAA36Oxc5HbbjN11Ci7OWLFuc4OAADEBo2di7Rta+qUKXZzRNPmzaaeeKLdHAAABAGNHWLqySdN5fo6AABij+VOXKZDB7O1WHa2lJ5uO03VOadfi4ulBF5GwBKWOwEQFPypdZm77jL16aft5og2mjoAAGKPETuXKSoya9pVq2YWLfaytWulo4+WMjMje8UCNjBiByAoGEdxmWrVTC0qspsjGh59tGwFAACxxYidC6WmSjt2SAUFkR0pvMi5vq6khKVOYBcjdgCCghE7F7rzTlNfe81ujmihqYMtWVlZysjIUGZmpu0oABAXjNi50LZtUoMG0vHHS19+aTtN5axcKbVuLZ19tvT227bTIOgYsQMQFIzYuVD9+qZ+9ZXdHFXxwAOmPvKI3RwAAAQJI3Yu5fX135z8fHfBDRixAxAUHmwZguH2202dMcNuDgAA4B00di7lNHYPPWQ3R2UsWGBq//52cwAAEDScinUxr57O7NTJLEj8zTfSMcfYTgNwKhZAcDBi5wFea+ycXSZo6gAAiC8aOxe75RZTZ82ym6Midu0y1csLKwMA4FU0di52112mDhtmN0dFPPGEqf/8p90cAAAEEdfYuZzXrrNz8hYWRva9BWzjGjsAQcGIncs5pzTz8+3mqCiaOgAA4o/GzuUee8zUZ5+1m6M8Fi829bzz7OYAACCoOBXrcnv2SCkpUnKytHu37TSHdvrp0gcfSF9/LR13nO00QASnYgEEBY2dB3jlOjuv5ETw0NgBCApOxXqAM/qVnW03x6F47RpAAAD8iMbOA556ytR77rGb41BGjjT16aft5gAAIMg4FesRbj/N6eQrKpISE+1mAfbFqVgAQcGInce4sbHbOxNNHQAA9tDYecQdd5j6yit2cxzIyy+betttdnMAABB0nIr1iNxcKS3N3H75xXaashISzKhdXp5Uu7btNMD+OBULICgYsfMI529RTo7dHAfivDSgqQMAwC4aOw/JzDR11Sq7Ofb2/vum9u1rNwdwIFlZWcrIyFCm88MDAD7HqVgPWb5cOvFEqUsXaf5822mMunXNqeFNm6RGjWynAQ6MU7EAgoLGzmPctuyJ2/IAB0JjByAoOBXrMUcfberGjXZzSNK0aaZecYXdHAAAwKCx85iXXjJ1yBC7OaTIdXX//KfdHAAAwOBUrAe54fRnOGyWObGdAygPTsUCCApG7DyoZk1Tt2yxl2HcOFNZlBgAAPdgxM6D3n1X6tZN6tFDevNNOxmcUcPdu6XkZDsZgPJixA5AUNDYeZTN07G7dkm1atl7fqCiaOwABAWnYj2qbVtTly+P/3PfcIOpzh6xAADAHRix86hvv5VatpRatJC+/z6+z+2GyRtARTBiByAoGLHzqGOOMfWHH+L7vAsWmNqxY3yfFwAAHB6NnYddfrmpr70Wv+fs0sXUGTPi95wAAKB8OBXrYTt3SrVrm+N4fBX37JFSUuL3fEC0cCoWQFAwYudhzsxUKT5bjF16qanjx8f+uQAAQMUxYudxCxdKJ58s/f730hdfxO552GkCXsaIHYCgYMTO4zp1MnXlSqmkJHbP8/DDpt5yS+yeAwAAVA2NnQ/cf7+pQ4fG7jkefNDUUaNi9xwAAKBqaOx8wGm6nn46NqdJX33V1D/+MXI6FqiMcDisBx98UI0bN1ZKSoq6du2qlStXlvvjR4wYoVAopKGxfBUDAB7Gn2kfCIWkwYPN8QMPRP/xnUkTs2ZF/7ERLI8//riefPJJjR49WosWLVKjRo101llnKS8v77Afu2jRIo0dO1Zt2rSJQ1IA8CYaO5/4179MfeSR2DzuOeeUnYULVFQ4HNZTTz2le++9V3369FHr1q310ksvadeuXZo0adIhP3bHjh267LLLNG7cONWtWzdOiQHAe2jsfCIhQbrqKnPsjN5VVXGxdPPN5njmzOg8JoLrhx9+UHZ2trp37156X1JSks444wwtcLY0OYjBgwerR48e6tatW7meKz8/X7m5uWVuABAENHY+8vzzpo4ZI0Xj71jnzqY+8IBUvXrVHw/Blp2dLUlKT08vc396enrpvx3Ia6+9piVLlmjEiBHlfq4RI0YoLS2t9NasWbPKhQYAj6Gx85FQSJo61Rw3bVq1x1q6VFq0yBw7kzOAinjllVdUu3bt0lthYaEkKRQKlXm/cDi8332OdevW6eabb9bEiROVnJxc7uceNmyYcnJySm/r1q2r/H8EADyEBYp9yPkb+fbb0tlnV+0xvvlGOuaY6ORCsOTl5emnn34qfTs/P1+tW7fWkiVL1K5du9L7L7jgAh1xxBF66aWX9nuMGTNmqHfv3kpMTCy9r7i4WKFQSAkJCcrPzy/zbwfDAsUAgqKa7QCIvvXrzYjdn/4kbdggNW5csY8/9VRTzzuPpg6Vl5qaqtTU1NK3w+GwGjVqpDlz5pQ2dgUFBZo3b55Gjhx5wMc488wztWLFijL3DRo0SK1atdJdd91VrqYOAIKExs6HmjSRpk+Xevc2x0VFUnn//t19t+Rcx87yJogmZ/254cOHq2XLlmrZsqWGDx+umjVr6lJnTR2ZZq53794aMmSIUlNT1bp16zKPU6tWLdWvX3+/+wEANHa+1auXdP310rPPStWqlW/h4tdfl5yBk927I6djgWi58847tXv3bt144436+eef1alTJ82ePbvMyN53332nrVu3WkwJAN7FNXY+l5oq7dhhjgsLTZN3IMOHS/fea47XrzcjfYBfcI0dgKBgVqzP5eZKzkoP1aub5q2kJPLvH3xgRuacpm71apo6AAC8isbO50Ihae1ayVkCbPhwc71dKGRup59u7j/hBGnXLqllS3tZAQBA1dDYBcTdd5udJIYNk5KSIvf36CFt3CgtXy6lpNjLBwAAqo5r7AD4HtfYAQgKRuwAAAB8gsYOAKnpOD8AAACPSURBVADAJ2jsAAAAfILGDgAAwCdo7AAAAHyCxg4AAMAnaOwAAAB8gsYOAADAJ1igGIDvhcNh5eXlKTU1VaFQyHYcAIgZGjsAAACf4FQsAACAT9DYAQAA+ASNHQAAgE/Q2AEAAPgEjR0AAIBP0NgBAAD4BI0dAACAT9DYAQAA+ASNHQAAgE/Q2AEAAPjE/wNDFcxPHaAvYwAAAABJRU5ErkJggg==\n", "text/plain": [ "Graphics object consisting of 1 graphics primitive" ] }, "execution_count": 63, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from ore_algebra.analytic.function import DFiniteFunction\n", "P. = QQ[]\n", "A. = OreAlgebra(P)\n", "f = DFiniteFunction(Dx^2 - x,\n", " [1/(gamma(2/3)*3^(2/3)), -1/(gamma(1/3)*3^(1/3))],\n", " name='my_Ai')\n", "f.plot((-5,5))" ] } ], "metadata": { "celltoolbar": "Slideshow", "kernelspec": { "display_name": "SageMath 8.8.beta5", "language": "", "name": "sagemath" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.15" }, "rise": { "controls": false, "progress": false, "slideNumber": false, "transition": "none" } }, "nbformat": 4, "nbformat_minor": 2 }