{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Differential operators"
   ]
  },
  {
   "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": [
    "Pols.<z> = PolynomialRing(QQ)\n",
    "DiffOps.<Dz> = OreAlgebra(Pols)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Univariate Ore algebra in Dz over Univariate Polynomial Ring in z over Rational Field"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "DiffOps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "cos(z)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Dz(sin(z))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "z*Dz + 1"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Dz*z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "RecOps.<Sz> = OreAlgebra(Pols)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(z + 1)*Sz"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Sz*z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "sin(z + 1)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Sz(sin(z))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Toy Examples 1: Numerical Evaluation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[2.7182818285 +/- 4.37e-11]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(Dz - 1).numerical_solution(ini=[1], path=[0, 1], eps=1e-10) # variants..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "⤷ Note that the output is an interval!<br/><br/>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[2.7182818285 +/- 4.37e-11]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "my_interval = _\n",
    "my_interval"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[13.731185225 +/- 2.46e-10]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(Dz - 1).numerical_solution(ini=[my_interval/3], path=[0, my_interval])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.17675339323955288 +/- 4.56e-18]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dop = Dz^2 - z  # Airy Ai\n",
    "ini = [1/3*3^(1/3)/gamma(2/3), -1/2*3^(1/6)*gamma(2/3)/pi]\n",
    "dop.numerical_solution(ini, [0, -100])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Toy Examples 2: Paths"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dop = Dz*z*Dz\n",
    "dop(log(z))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.69314718055994531 +/- 1.23e-18]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dop.numerical_solution(ini=[0, 1], path=[1, 2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "#dop.numerical_solution(ini=[0, 1], path=[1, -1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[+/- 6.62e-28] + [3.1415926535897932 +/- 3.88e-17]*I"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dop.numerical_solution(ini=[0, 1], path=[1, i, -1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[+/- 6.62e-28] + [-3.1415926535897932 +/- 3.88e-17]*I"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dop.numerical_solution(ini=[0, 1], path=[1, -i, -1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Local monodromy matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-z^2 + z)*Dz^2 + (-11/6*z + 1)*Dz - 1/6"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a, b, c = 1/2, 1/3, 1   # ₂F₁(a,b;c;z)\n",
    "dop = z*(1-z)*Dz^2 + (c - (a + b + 1)*z)*Dz - a*b\n",
    "dop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-1) * (z - 1) * z"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dop.leading_coefficient().factor()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[1.00000 +/- 8.31e-8] + [-0.646840 +/- 3.98e-7]*I           [+/- 1.37e-7] + [2.19209 +/- 1.97e-6]*I]\n",
       "[        [+/- 1.10e-7] + [-0.190869 +/- 4.45e-7]*I  [1.00000 +/- 2.38e-7] + [0.646840 +/- 5.60e-7]*I]"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dop.numerical_transition_matrix([1/2, i/2, -1/2, -i/2, 1/2], 1e-5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[ [1.38287 +/- 2.35e-6] + [-0.663154 +/- 3.06e-7]*I   [1.03532 +/- 3.33e-6] + [-1.79323 +/- 2.53e-6]*I]\n",
       "[[-0.326494 +/- 5.88e-7] + [0.565505 +/- 1.43e-7]*I   [0.117128 +/- 4.56e-7] + [1.52918 +/- 5.37e-7]*I]"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dop.numerical_transition_matrix([1/2, 1-i/2, 3/2, 1+i/2,1/2], 1e-5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Local monodromy by singular connection\n",
    "\n",
    "Local monodromy = formal monodromy + singular connection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "a, b, c = 1/2, 1/3, 1   # ₂F₁(a,b;c;z)\n",
    "dop = z*(1-z)*Dz^2 + (c - (a + b + 1)*z)*Dz - a*b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[1.00000 +/- 8.31e-8] + [-0.646840 +/- 3.98e-7]*I           [+/- 1.37e-7] + [2.19209 +/- 1.97e-6]*I]\n",
       "[        [+/- 1.10e-7] + [-0.190869 +/- 4.45e-7]*I  [1.00000 +/- 2.38e-7] + [0.646840 +/- 5.60e-7]*I]"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dop.numerical_transition_matrix([1/2, I/2, -1/2, -I/2, 1/2], 1e-5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[log(z) + 1/6*z*log(z) + 1/2*z + 1/12*z^2*log(z) + 41/144*z^2 + 35/648*z^3*log(z) + 167/864*z^3,\n",
       " 1 + 1/6*z + 1/12*z^2 + 35/648*z^3]"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dop.local_basis_expansions(0, 4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[1.000000 +/- 2.58e-8] + [-0.6468403 +/- 4.39e-8]*I           [+/- 7.44e-9] + [2.1920882 +/- 7.55e-8]*I]\n",
       "[        [+/- 1.81e-8] + [-0.19086933 +/- 8.92e-9]*I  [1.000000 +/- 3.44e-8] + [0.6468403 +/- 4.17e-8]*I]"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mat = dop.numerical_transition_matrix([0, 1/2], 1e-7)\n",
    "mon = matrix([[1, 0], [CBF(2*pi*I), 1]])\n",
    "mat*mon*~mat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Pólya walks in dimension 15\n",
    "\n",
    "(Thanks to B. Salvy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "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": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from ore_algebra.examples import polya\n",
    "dim = 15\n",
    "polya.dop[dim]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "polya.dop[dim].leading_coefficient()(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "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": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "polya.dop[dim].local_basis_monomials(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.03586962312535651422940427078451975465049 +/- 2.68e-42]"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ini = [0]*(dim - 1) + [1]\n",
    "v = polya.dop[dim].numerical_solution(ini, [0,1/(2*dim)], 1e-60)\n",
    "(v - 1)/v"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(The evaluation point is singular, but the solution has a finite limit.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Face-centered cubic lattices\n",
    "[Guttmann 2009, Broadhurst 2009, Koutschan 2013]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(27122036833024*z^43 + 8208413201024064*z^42 + 1028987679702510976*z^41 + 75518451137118783792*z^40 + 3743195619381989907184*z^39 + 135369638077546936261428*z^38 + 3745615314367420203992832*z^37 + 81811619367860049045984675*z^36 + 1440466637248203913774334250*z^35 + 20724331113040275023719172850*z^34 + 245446627541652046097792768214*z^33 + 2395828801191215780780578117794*z^32 + 19147407470673111231862249418166*z^31 + 122863963621496746370188659696702*z^30 + 602621255648485924378700672331054*z^29 + 1935192664301617476137337671088360*z^28 + 694152712036783264243644290673234*z^27 - 39030042885818935455901289133872622*z^26 - 297645962803933196564873733670191774*z^25 - 1329742929728007215704002549281591538*z^24 - 3903989614825648819224432657208727646*z^23 - 6038015534019664017777438417359311914*z^22 + 7565280951156009750992823479550694170*z^21 + 83328126336960183101771239549883786325*z^20 + 297859836972471180382017327162905955900*z^19 + 681226694393685252017130073908325840500*z^18 + 1055504316932586226613390044310017920000*z^17 + 974982625144110654834660688990434600000*z^16 + 80921481727424794623135930623472000000*z^15 - 1245692778975371208980497936649580000000*z^14 - 1877612972166046542841525891548000000000*z^13 - 1186201691981014544058180007080000000000*z^12 - 4424139263701398029187830400000000000*z^11 + 526588498855023559134177312000000000000*z^10 + 410999234738834010247469568000000000000*z^9 + 174333012213810958051184640000000000000*z^8 + 29950530354743793153638400000000000000*z^7 + 2276991208061142220800000000000000000*z^6)*Dz^8 + (1600200173148416*z^42 + 478782978712278912*z^41 + 58815380786135567104*z^40 + 4218590040421804170816*z^39 + 204216444469816446653424*z^38 + 7214118624119937529541160*z^37 + 195106070712453547506798168*z^36 + 4168870319524368197533959000*z^35 + 71874412795312940511795668940*z^34 + 1013528039913249207367842378270*z^33 + 11775924181048893848357395670676*z^32 + 112862055818213392356279768225402*z^31 + 886340494836475569866741139358344*z^30 + 5592675973186567437279733685351646*z^29 + 26983635759333711243427828354079724*z^28 + 85059388463264142313662526542420618*z^27 + 24816956833181644480403313150735864*z^26 - 1739731529923503295984796806526752758*z^25 - 13215685421423157401833903137021991092*z^24 - 60101514732517779329542749898893453858*z^23 - 187406121933740017212741167478185137320*z^22 - 367088786736715063908412462166156515566*z^21 - 136331238303988349001415414181532146340*z^20 + 2052937632229799753666758504303681446150*z^19 + 8942220864711302092023950168348534856300*z^18 + 22112779083456047399791690319673356808000*z^17 + 36662299830964853548300895468723502480000*z^16 + 38663936209054739955701649076784708400000*z^15 + 15575841209632684184725074680551176000000*z^14 - 23399775927110778754739301057544560000000*z^13 - 45957581844555068108338692961807200000000*z^12 - 32525005285459811112066289505232000000000*z^11 - 3661218292337523929046664304640000000000*z^10 + 10970395301506611292814537164800000000000*z^9 + 9713112405197935942595533824000000000000*z^8 + 4260978719420783892503777280000000000000*z^7 + 756467079122535026688000000000000000000*z^6 + 59201771409589697740800000000000000000*z^5)*Dz^7 + (35882454730090752*z^41 + 10612604051614486656*z^40 + 1276532600942212775168*z^39 + 89393980129433032096320*z^38 + 4221606838983473228197008*z^37 + 145494567985766484898923048*z^36 + 3840828004490920060950969480*z^35 + 80160062388267727172211985080*z^34 + 1350855094398006902682870922050*z^33 + 18631082892630536824222949409585*z^32 + 211815796834464054711973645322142*z^31 + 1986708322085667572665525016037411*z^30 + 15263082383031406770429022758762048*z^29 + 94068732852089205756130773605094705*z^28 + 441055376229095921513357130918811338*z^27 + 1319636945498761264973744224282378779*z^26 - 137626809673226795399591264079041112*z^25 - 31072001737970299221405533198706303141*z^24 - 226886176666918560987240200768631693150*z^23 - 1033954017266382248984767586852072344191*z^22 - 3356732946224373601649087937349109785896*z^21 - 7573126212785007618891225542456994124245*z^20 - 9076459539413303184641722134776573895810*z^19 + 10278671248090335377408918358815408788425*z^18 + 85149274357043292385925033653294291853550*z^17 + 240689360358498296007939096187740586134000*z^16 + 429409878921957648790555775268242743350000*z^15 + 495779225046771906420255540348281344800000*z^14 + 287121363379312616871562346484465378000000*z^13 - 119682652007548350954457856750250720000000*z^12 - 395683465592680867401293480616198000000000*z^11 - 327383462755042385949747691240824000000000*z^10 - 86642575450501391066787202019520000000000*z^9 + 59704683972170679548931977222400000000000*z^8 + 72511610277412390990839363072000000000000*z^7 + 33882896755872071956886261760000000000000*z^6 + 6311156771304917325766656000000000000000*z^5 + 512323021813756999680000000000000000000*z^4)*Dz^6 + (390720062616543744*z^40 + 114216661424360307456*z^39 + 13440822351615963069696*z^38 + 917965180366474611870720*z^37 + 42237673932263775988570560*z^36 + 1418218839310481932976078400*z^35 + 36487206836408001197689910640*z^34 + 742504152062765602237759452720*z^33 + 12205694666919011642462560650930*z^32 + 164251022394443778986763736539405*z^31 + 1821767115612836434053737404755054*z^30 + 16657938302824007267243724365434191*z^29 + 124533460849620200009711445328730256*z^28 + 743593796442908540070532245488378205*z^27 + 3336004607088107531634361889061221370*z^26 + 9020222161854736084202629824390547047*z^25 - 9198824722943404205447421299404277112*z^24 - 279449868802402514175677041789492570017*z^23 - 1907863427661939885576723126598906643790*z^22 - 8574083646475710050757565542672979674555*z^21 - 28405587296847231070183606856583770811720*z^20 - 69574258175312955514440713973653616428745*z^19 - 114991436892487711669937849824912517430330*z^18 - 70378017201579863364495432167182725333675*z^17 + 261641966501089147843656083216157842879550*z^16 + 1049410824795136384837467209810025539400000*z^15 + 2089663780964272997600159898811800513390000*z^14 + 2597860026796363803313313450929745040000000*z^13 + 1759136834585156085432113720072647266000000*z^12 - 134675614991223713108740928290811280000000*z^11 - 1578996098791370746284707453439169200000000*z^10 - 1556811681322720025894531955998040000000000*z^9 - 678908917613499441761342434795200000000000*z^8 + 3149221592935046990103881875200000000000*z^7 + 192407344459752425261121833472000000000000*z^6 + 104126840210444782060628951040000000000000*z^5 + 20982276643045393599012864000000000000000*z^4 + 1787438098327996643328000000000000000000*z^3)*Dz^5 + (2192816677949990400*z^39 + 633490213477308768000*z^38 + 72864986011484455353600*z^37 + 4847486869795537260532800*z^36 + 217014017048761645614816000*z^35 + 7088016995800124707996090560*z^34 + 177409657131610270482016190640*z^33 + 3513089912736584156549238676620*z^32 + 56201587732740449959670675451690*z^31 + 735847759326730529024504240988015*z^30 + 7934411063073314432988482485900500*z^29 + 70405171912630286359945571896774110*z^28 + 508882813920850610699235633677324220*z^27 + 2913386772290646501282812655546011475*z^26 + 12237989774964463385062890926963215950*z^25 + 27303163874555616052155898475524386210*z^24 - 84400724272601405065271773264397209530*z^23 - 1309329548085768562973129072537724164955*z^22 - 8229269199062442444264260234977847805360*z^21 - 35948740918844475140318574840001115213670*z^20 - 119246681320333134593914488403328142970080*z^19 - 304273308297438630162837099041285050546455*z^18 - 577145030895901790697311126386896036767490*z^17 - 707778167790136602728038144967670916837350*z^16 - 153341406553907334245470038125935935813900*z^15 + 1611202920628825942940219406515876419542000*z^14 + 4188993616205017046899739124211544543460000*z^13 + 5699626392082018037453259396194906388000000*z^12 + 4145140187203309836183311398252469964000000*z^11 + 95068892397133773199630362250506960000000*z^10 - 3224476947136810362409243540972029600000000*z^9 - 3636835528138928302767664987399536000000000*z^8 - 2103918711236593285373196111532800000000000*z^7 - 551590602425414384164117858252800000000000*z^6 + 116215925694410902420898178048000000000000*z^5 + 112838527684482371017981562880000000000000*z^4 + 26640440911690813973554790400000000000000*z^3 + 2427272627793177607372800000000000000000*z^2)*Dz^4 + (6219625486549063680*z^38 + 1775531336308022522880*z^37 + 199409996635132589752320*z^36 + 12904862497592448920163840*z^35 + 561222248755128125708191680*z^34 + 17798695421072697669468739680*z^33 + 432530168604725658189210596640*z^32 + 8315189920333341531658617695280*z^31 + 129103723904595771409928232487740*z^30 + 1639190738531986170699647097803790*z^29 + 17111040709840823035760757618682440*z^28 + 146518901443861803658771329866897880*z^27 + 1015534278806669843159745327151252620*z^26 + 5496338276053076075068754467310102760*z^25 + 20890574209714927539267068315744951640*z^24 + 30164609970591947189827076234922007050*z^23 - 289127416281529376142095631015519267120*z^22 - 3053669279873063346793150591937974700130*z^21 - 17566486109105161467894504789161406270600*z^20 - 73673650638461574538679743097050051115220*z^19 - 240407438967557398913317296975336574702980*z^18 - 619168293687639511251067273975020197114100*z^17 - 1241225303290460623795990859905959226579320*z^16 - 1835553795134837646262350779261931882894750*z^15 - 1670314602837141110845640706031012073555700*z^14 + 30668748618622961862750876075275710932000*z^13 + 2931594155313390328935716187001614568260000*z^12 + 4919456458899666498684069708388548285600000*z^11 + 3777365646243762653104795884206143332000000*z^10 + 68195639154415674514017863641593600000000*z^9 - 3118810937523253726235666782666096800000000*z^8 - 3771833787399616704258908808294288000000000*z^7 - 2488444831930996824908954989144320000000000*z^6 - 893959103422093803323305262745600000000000*z^5 - 100083705719332806676962561024000000000000*z^4 + 25642030845677374764180172800000000000000*z^3 + 10280761031833373040014131200000000000000*z^2 + 1033754008459758568243200000000000000000*z)*Dz^3 + (8133356405487237120*z^37 + 2294131782043664317440*z^36 + 251295328534762193633280*z^35 + 15795453015240816970091520*z^34 + 666093618246765502077439680*z^33 + 20469375712416843040909376160*z^32 + 481817338267338639783328749120*z^31 + 8967973126212020032517991216960*z^30 + 134696914854304536722281866954300*z^29 + 1651833654984876079820125885678650*z^28 + 16607490026343429532811575311949230*z^27 + 136263869454304799146859253346813455*z^26 + 895865319327471447638111289873238710*z^25 + 4489711074264384906529925990254793265*z^24 + 14491852283494577826654003932547711690*z^23 - 168509066471194546983174648133542750*z^22 - 386265894549826881229123104470731096440*z^21 - 3163259131060568584546113343781987561220*z^20 - 16636182069413821170544684047556220568150*z^19 - 66246740089393676080981537130378090658525*z^18 - 209080245872850631566312137449619561543730*z^17 - 529097465740104776391772834675033946593335*z^16 - 1065038620757313929391639361032363453750930*z^15 - 1653651644685620142167009422124022555221700*z^14 - 1829383474513975929874027770563298831967800*z^13 - 1088709896836905806660284149277621568328000*z^12 + 437384067337328886944483963336952904080000*z^11 + 1678380365365432006625451473236269012000000*z^10 + 1564385355592027935922683162898655112000000*z^9 + 303607398715325954207032303663107840000000*z^8 - 929554263384554771136483801584745600000000*z^7 - 1334658535726482371536908049179648000000000*z^6 - 941977534006524837182879263564800000000000*z^5 - 350297653852787677589275016601600000000000*z^4 - 64286241473892148234640584704000000000000*z^3 - 3745850575154857509858508800000000000000*z^2 + 668706777983396247404544000000000000000*z + 72863718657956551065600000000000000000)*Dz^2 + (3964156903514787840*z^36 + 1104718489963413534720*z^35 + 117871088739930352834560*z^34 + 7183287516644479615795200*z^33 + 293105835218942903781855360*z^32 + 8706572378734984776799502400*z^31 + 197949776138115866133849254880*z^30 + 3555494624146318548046453851120*z^29 + 51457898672013865098111291247320*z^28 + 606522834979531840684521625237020*z^27 + 5835366836846027182876920856348950*z^26 + 45455202501826358974606219981974015*z^25 + 279153404467502062948531557838260750*z^24 + 1252275399372837134061507042628908795*z^23 + 2943182802923552038161307584706940070*z^22 - 10483513115206289398510413216920199750*z^21 - 169948182933507479161257565568616530700*z^20 - 1154969594776277160649077785983820553870*z^19 - 5548694490781020038019823355124585193590*z^18 - 20745229517577451272377158241970915439245*z^17 - 62232963928794638659423069651761724690290*z^16 - 150810045901978932864163493046405461262105*z^15 - 292528626523005661629390236883046859976150*z^14 - 441395096063183148839008172248580337780300*z^13 - 484123578764537043031861206473715269343000*z^12 - 312976584649334763451810663858004196420000*z^11 + 31415133499909950234831915395869293600000*z^10 + 330738460087674555468491482629558468000000*z^9 + 391096978918364972225128472061480072000000*z^8 + 232460170425948027345434850305279520000000*z^7 + 18060134934884299834847099345568000000000*z^6 - 100213400891192102370293326036992000000000*z^5 - 83859064515985136495903099458560000000000*z^4 - 29487051768804049058487328972800000000000*z^3 - 5516482719743042228472840192000000000000*z^2 - 379553789663505987184558080000000000000*z + 8379327645665003372544000000000000000)*Dz + 410085196915322880*z^35 + 112905266474211563520*z^34 + 11716692637614961489920*z^33 + 690817401287078917363200*z^32 + 27204862643846611522761600*z^31 + 778811406918247228618497600*z^30 + 17044384124115240781429792800*z^29 + 294245234066850000428339092800*z^28 + 4083424587805117060272476125800*z^27 + 45973029549119796235386142827300*z^26 + 419695598890898253203455876749930*z^25 + 3064297176174091671717985958725620*z^24 + 17169584489259696388755804636033570*z^23 + 64581771961684810077279475394020500*z^22 + 51714221934272099420476126216766700*z^21 - 1473967391504437899277380487903179960*z^20 - 14237554341321335335392023192872385940*z^19 - 83216340134393115016834220980384454340*z^18 - 364019154328107562568847906822488063550*z^17 - 1261571478513401088177035093275526304300*z^16 - 3528341032098896995323439017117956856150*z^15 - 7964369518593778029521056070442794466900*z^14 - 14280500726162786254712841163875001728600*z^13 - 19534653115686342543580831960941978918000*z^12 - 18398783334222380084238012428704731960000*z^11 - 7553741785990309357234054786177488000000*z^10 + 8887432309419522403983976171775697600000*z^9 + 21137039158366320685856256980012112000000*z^8 + 22682693553934804690446647295508800000000*z^7 + 14938183428146261190546354671616000000000*z^6 + 4690246528584816329940199400448000000000*z^5 - 872829008634785738926162452480000000000*z^4 - 1104327940779745890150773145600000000000*z^3 - 353898708207580856772919296000000000000*z^2 - 52079342910774444874137600000000000000*z - 2428790621931885035520000000000000000"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from ore_algebra.examples import fcc\n",
    "fcc.dop6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1.02774910062749883985936367927396850209243990900114872425172166 +/- 4.37e-63] + [+/- 5.14e-73]*I"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fcc.dop6.numerical_solution([0, 0, 0, 0, 0, 1, 0, 0], [0, 1], 1e-60, assume_analytic=True)"
   ]
  },
  {
   "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": 32,
   "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": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from ore_algebra.examples import iint\n",
    "i = 69; iint.word[i]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "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": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dop = iint.diffop(iint.word[i])\n",
    "dop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "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": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "iint.ini[i]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.25000000000000000?]"
      ]
     },
     "execution_count": 35,
     "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": 36,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[-0.97080469562493124056011987537954344846933233520382808407829772093922068543247104239693854364491332570 +/- 9.81e-102] + [+/- 2.77e-102]*I"
      ]
     },
     "execution_count": 36,
     "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",
    "dop.numerical_solution(iint.ini[i], path, 1e-100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Volumes of semi-algebraic sets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "· PF eq order = 3\n",
      "· interval: [-3, -1]\n",
      "·· slice #1: ρ = -2971/1024\n",
      "··· PF eq order = 2\n",
      "··· interval: [-0.7629341015437308?, 0.7629341015437308?]\n",
      "···· slice #1: ρ = -333/512\n",
      "····· slice length = [0.458464611594969357299489644474 +/- 7.93e-31]\n",
      "···· slice #2: ρ = -3/16\n",
      "····· slice length = [0.84045215889193560180330863340 +/- 3.84e-30]\n",
      "···· integrating PF equation over [-0.7629341015437308?, 0.7629341015437308?]...\n",
      "···· ...piece volume = [1.042502360545567477912266674 +/- 6.60e-28] \n",
      "··· slice volume = [1.042502360545567477912266674 +/- 6.60e-28]\n",
      "·· slice #2: ρ = -1609/1024\n",
      "··· PF eq order = 2\n",
      "··· interval: [-2.555592041400167?, 2.555592041400167?]\n",
      "···· slice #1: ρ = -813/512\n",
      "····· slice length = [1.94451736178767524728375194202 +/- 3.84e-30]\n",
      "···· slice #2: ρ = 1841/1024\n",
      "····· slice length = [1.84355214687635799781072006275 +/- 4.16e-30]\n",
      "···· integrating PF equation over [-2.555592041400167?, 2.555592041400167?]...\n",
      "···· ...piece volume = [8.9321241233688167691817030 +/- 8.10e-26] \n",
      "··· slice volume = [8.9321241233688167691817030 +/- 8.10e-26]\n",
      "·· slice #3: ρ = -1447/1024\n",
      "··· PF eq order = 2\n",
      "··· interval: [-2.646353743028272?, 2.646353743028272?]\n",
      "···· slice #1: ρ = 1079/1024\n",
      "····· slice length = [1.94287384559422284208929621305 +/- 4.71e-30]\n",
      "···· slice #2: ρ = 687/512\n",
      "····· slice length = [1.99736127598458294060695738287 +/- 4.40e-30]\n",
      "···· integrating PF equation over [-2.646353743028272?, 2.646353743028272?]...\n",
      "···· ...piece volume = [9.05299012102810075746420 +/- 2.37e-24] \n",
      "··· slice volume = [9.05299012102810075746420 +/- 2.37e-24]\n",
      "·· integrating PF equation over [-3, -1]...\n",
      "·· ...piece volume = [13.03207391697359297216659 +/- 2.58e-24] \n",
      "· interval: [-1, 0]\n",
      "·· slice #1: ρ = -1021/1024\n",
      "··· PF eq order = 2\n",
      "··· interval: [-2.829461219372189?, -0.07649046954459323?]\n",
      "···· slice #1: ρ = -5761/2048\n",
      "····· slice length = [0.35112251808149490976859102440 +/- 5.00e-30]\n",
      "···· slice #2: ρ = -1375/512\n",
      "····· slice length = [1.00469443765590664671527499614 +/- 4.94e-30]\n",
      "···· integrating PF equation over [-2.829461219372189?, -0.07649046954459323?]...\n",
      "···· ...piece volume = [3.981064601636449324301241352 +/- 2.41e-28] \n",
      "··· interval: [-0.07649046954459323?, 0.07649046954459323?]\n",
      "···· slice #1: ρ = -29/1024\n",
      "····· slice length = 0\n",
      "···· slice #2: ρ = 45/2048\n",
      "····· slice length = 0\n",
      "···· integrating PF equation over [-0.07649046954459323?, 0.07649046954459323?]...\n",
      "···· ...piece volume = 0 \n",
      "··· interval: [0.07649046954459323?, 2.829461219372189?]\n",
      "···· slice #1: ρ = 1867/1024\n",
      "····· slice length = [1.99389627375147032942435597344 +/- 5.69e-30]\n",
      "···· slice #2: ρ = 3895/2048\n",
      "····· slice length = [1.97816248534009683807163846683 +/- 7.26e-30]\n",
      "···· integrating PF equation over [0.07649046954459323?, 2.829461219372189?]...\n",
      "···· ...piece volume = [3.981064601636449324301241 +/- 6.24e-25] \n",
      "··· slice volume = [7.962129203272898648602483 +/- 5.68e-25]\n",
      "·· slice #2: ρ = -233/512\n",
      "··· PF eq order = 2\n",
      "··· interval: [-2.965283106239012?, -0.8904515147645516?]\n",
      "···· slice #1: ρ = -1489/512\n",
      "····· slice length = [0.66221351268681499597164867676 +/- 1.75e-30]\n",
      "···· slice #2: ρ = -83/32\n",
      "····· slice length = [1.54769896172188521586787386448 +/- 4.71e-30]\n",
      "···· integrating PF equation over [-2.965283106239012?, -0.8904515147645516?]...\n",
      "···· ...piece volume = [3.249392395759612579727661331 +/- 5.53e-28] \n",
      "··· interval: [-0.8904515147645516?, 0.8904515147645516?]\n",
      "···· slice #1: ρ = -581/1024\n",
      "····· slice length = 0\n",
      "···· slice #2: ρ = 207/512\n",
      "····· slice length = 0\n",
      "···· integrating PF equation over [-0.8904515147645516?, 0.8904515147645516?]...\n",
      "···· ...piece volume = 0 \n",
      "··· interval: [0.8904515147645516?, 2.965283106239012?]\n",
      "···· slice #1: ρ = 17/16\n",
      "····· slice length = [1.07223199863158322943862248025 +/- 7.23e-30]\n",
      "···· slice #2: ρ = 1193/512\n",
      "····· slice length = [1.85477515115246255326540162547 +/- 2.45e-30]\n",
      "···· integrating PF equation over [0.8904515147645516?, 2.965283106239012?]...\n",
      "···· ...piece volume = [3.24939239575961257972766133 +/- 6.63e-27] \n",
      "··· slice volume = [6.49878479151922515945532266 +/- 7.54e-27]\n",
      "·· slice #3: ρ = -267/1024\n",
      "··· PF eq order = 2\n",
      "··· interval: [-2.988647438500854?, -0.9654084688139601?]\n",
      "···· slice #1: ρ = -1293/512\n",
      "····· slice length = [1.68484755897555280649424888121 +/- 3.69e-30]\n",
      "···· slice #2: ρ = -643/512\n",
      "····· slice length = [1.39340849025983296334216722480 +/- 7.01e-30]\n",
      "···· integrating PF equation over [-2.988647438500854?, -0.9654084688139601?]...\n",
      "···· ...piece volume = [3.1753548383096139183976395 +/- 5.26e-26] \n",
      "··· interval: [-0.9654084688139601?, 0.9654084688139601?]\n",
      "···· slice #1: ρ = -885/1024\n",
      "····· slice length = 0\n",
      "···· slice #2: ρ = -53/64\n",
      "····· slice length = 0\n",
      "···· integrating PF equation over [-0.9654084688139601?, 0.9654084688139601?]...\n",
      "···· ...piece volume = 0 \n",
      "··· interval: [0.9654084688139601?, 2.988647438500854?]\n",
      "···· slice #1: ρ = 599/256\n",
      "····· slice length = [1.87024322618968416215255633412 +/- 5.09e-30]\n",
      "···· slice #2: ρ = 715/256\n",
      "····· slice length = [1.18624194459034992382051850380 +/- 5.05e-30]\n",
      "···· integrating PF equation over [0.9654084688139601?, 2.988647438500854?]...\n",
      "···· ...piece volume = [3.175354838309613918397639546 +/- 8.66e-28] \n",
      "··· slice volume = [6.3507096766192278367952791 +/- 1.47e-26]\n",
      "·· integrating PF equation over [-1, 0]...\n",
      "·· ...piece volume = [6.707134885205124265502393 +/- 3.38e-25] \n",
      "· interval: [0, 1]\n",
      "·· slice #1: ρ = 203/1024\n",
      "··· PF eq order = 2\n",
      "··· interval: [-2.993442839790868?, -0.980153067176355?]\n",
      "···· slice #1: ρ = -757/512\n",
      "····· slice length = [1.72241534088442947993465158425 +/- 5.70e-30]\n",
      "···· slice #2: ρ = -129/128\n",
      "····· slice length = [0.462663836618875501846523585533 +/- 7.65e-31]\n",
      "···· integrating PF equation over [-2.993442839790868?, -0.980153067176355?]...\n",
      "···· ...piece volume = [3.160930248136342774702468365 +/- 7.77e-28] \n",
      "··· interval: [-0.980153067176355?, 0.980153067176355?]\n",
      "···· slice #1: ρ = -123/512\n",
      "····· slice length = 0\n",
      "···· slice #2: ρ = 11/1024\n",
      "····· slice length = 0\n",
      "···· integrating PF equation over [-0.980153067176355?, 0.980153067176355?]...\n",
      "···· ...piece volume = 0 \n",
      "··· interval: [0.980153067176355?, 2.993442839790868?]\n",
      "···· slice #1: ρ = 1263/1024\n",
      "····· slice length = [1.32112383947474290084045758787 +/- 4.32e-30]\n",
      "···· slice #2: ρ = 167/64\n",
      "····· slice length = [1.57409139853110946142179387741 +/- 7.08e-30]\n",
      "···· integrating PF equation over [0.980153067176355?, 2.993442839790868?]...\n",
      "···· ...piece volume = [3.16093024813634277470246836 +/- 9.07e-27] \n",
      "··· slice volume = [6.32186049627268554940493673 +/- 5.44e-27]\n",
      "·· slice #2: ρ = 697/1024\n",
      "··· PF eq order = 2\n",
      "··· interval: [-2.921762556064575?, -0.7325956825022900?]\n",
      "···· slice #1: ρ = -403/256\n",
      "····· slice length = [1.91709719706520579780672597983 +/- 7.20e-30]\n",
      "···· slice #2: ρ = -619/512\n",
      "····· slice length = [1.58082418408033357321882491781 +/- 4.99e-30]\n",
      "···· integrating PF equation over [-2.921762556064575?, -0.7325956825022900?]...\n",
      "···· ...piece volume = [3.4081996526641283573756781 +/- 8.27e-26] \n",
      "··· interval: [-0.7325956825022900?, 0.7325956825022900?]\n",
      "···· slice #1: ρ = -461/1024\n",
      "····· slice length = 0\n",
      "···· slice #2: ρ = 331/512\n",
      "····· slice length = 0\n",
      "···· integrating PF equation over [-0.7325956825022900?, 0.7325956825022900?]...\n",
      "···· ...piece volume = 0 \n",
      "··· interval: [0.7325956825022900?, 2.921762556064575?]\n",
      "···· slice #1: ρ = 193/128\n",
      "····· slice length = [1.87671102837248445873724227967 +/- 4.30e-30]\n",
      "···· slice #2: ρ = 2087/1024\n",
      "····· slice length = [1.97775167783927169187063719727 +/- 5.57e-30]\n",
      "···· integrating PF equation over [0.7325956825022900?, 2.921762556064575?]...\n",
      "···· ...piece volume = [3.4081996526641283573756781 +/- 7.74e-26] \n",
      "··· slice volume = [6.8163993053282567147513561 +/- 7.98e-26]\n",
      "·· slice #3: ρ = 799/1024\n",
      "··· PF eq order = 2\n",
      "··· interval: [-2.896752209411764?, -0.6254385363342537?]\n",
      "···· slice #1: ρ = -715/512\n",
      "····· slice length = [1.83275605410443441052480387291 +/- 5.91e-30]\n",
      "···· slice #2: ρ = -479/512\n",
      "····· slice length = [1.24712375723136725932529481765 +/- 6.41e-30]\n",
      "···· integrating PF equation over [-2.896752209411764?, -0.6254385363342537?]...\n",
      "···· ...piece volume = [3.51682354126377608771849909 +/- 7.00e-27] \n",
      "··· interval: [-0.6254385363342537?, 0.6254385363342537?]\n",
      "···· slice #1: ρ = -263/1024\n",
      "····· slice length = 0\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "···· slice #2: ρ = 631/1024\n",
      "····· slice length = 0\n",
      "···· integrating PF equation over [-0.6254385363342537?, 0.6254385363342537?]...\n",
      "···· ...piece volume = 0 \n",
      "··· interval: [0.6254385363342537?, 2.896752209411764?]\n",
      "···· slice #1: ρ = 859/1024\n",
      "····· slice length = [1.03941172948742980836593367454 +/- 4.30e-30]\n",
      "···· slice #2: ρ = 1149/512\n",
      "····· slice length = [1.85330459474314893076143365725 +/- 6.81e-30]\n",
      "···· integrating PF equation over [0.6254385363342537?, 2.896752209411764?]...\n",
      "···· ...piece volume = [3.51682354126377608771849909 +/- 7.88e-27] \n",
      "··· slice volume = [7.0336470825275521754369982 +/- 2.35e-26]\n",
      "·· integrating PF equation over [0, 1]...\n",
      "·· ...piece volume = [6.7071348852051242655023929 +/- 6.65e-26] \n",
      "· interval: [1, 3]\n",
      "·· slice #1: ρ = 1371/1024\n",
      "··· PF eq order = 2\n",
      "··· interval: [-2.684666581576908?, 2.684666581576908?]\n",
      "···· slice #1: ρ = -3001/2048\n",
      "····· slice length = [1.99977148445789850522530524014 +/- 4.53e-30]\n",
      "···· slice #2: ρ = 65/2048\n",
      "····· slice length = [1.50120032896689418013671677887 +/- 3.18e-30]\n",
      "···· integrating PF equation over [-2.684666581576908?, 2.684666581576908?]...\n",
      "···· ...piece volume = [9.0372366052727530219173661 +/- 8.24e-26] \n",
      "··· slice volume = [9.0372366052727530219173661 +/- 8.24e-26]\n",
      "·· slice #2: ρ = 549/256\n",
      "··· PF eq order = 2\n",
      "··· interval: [-2.097852644437507?, 2.097852644437507?]\n",
      "···· slice #1: ρ = -3371/2048\n",
      "····· slice length = [1.42160828066475677698884481060 +/- 5.46e-30]\n",
      "···· slice #2: ρ = -171/128\n",
      "····· slice length = [1.70021762587569529176377085343 +/- 2.41e-30]\n",
      "···· integrating PF equation over [-2.097852644437507?, 2.097852644437507?]...\n",
      "···· ...piece volume = [6.96835655296237511793174390 +/- 9.05e-27] \n",
      "··· slice volume = [6.96835655296237511793174390 +/- 9.05e-27]\n",
      "·· slice #3: ρ = 1411/512\n",
      "··· PF eq order = 2\n",
      "··· interval: [-1.185427815273714?, 1.185427815273714?]\n",
      "···· slice #1: ρ = -119/256\n",
      "····· slice length = [1.21377403589794682406668928262 +/- 2.43e-30]\n",
      "···· slice #2: ρ = 169/512\n",
      "····· slice length = [1.26255694914163068706242855126 +/- 6.90e-30]\n",
      "···· integrating PF equation over [-1.185427815273714?, 1.185427815273714?]...\n",
      "···· ...piece volume = [2.467431098167482394844393 +/- 3.53e-25] \n",
      "··· slice volume = [2.467431098167482394844393 +/- 3.53e-25]\n",
      "·· integrating PF equation over [1, 3]...\n",
      "·· ...piece volume = [13.032073916973592972166589 +/- 1.45e-25] \n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[39.47841760435743447533796 +/- 6.02e-24]"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import volume\n",
    "with open(\"torus.data\") as f:\n",
    "    vol = volume.doit(f, 100)\n",
    "vol"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Bonus Examples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Asymptotics of Apéry Numbers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0, 0.02943725152285942?, 33.97056274847714?]"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dop = (z^2*(z^2-34*z+1)*Dz^4 + 5*z*(2*z^2-51*z+1)*Dz^3 + (25*z^2-418*z+4)*Dz^2 + (15*z-117)*Dz + 1)\n",
    "sing = dop.leading_coefficient().roots(AA, multiplicities=False); sing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1/2*log(z)^2, log(z), 1, z]"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dop.local_basis_monomials(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "⤷ Solutions $y = y_0 + y_1 z + \\dots$ analytic at 0 are characterized by $y_0, y_1$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1,\n",
       " sqrt(z - 0.02943725152285942?),\n",
       " z - 0.02943725152285942?,\n",
       " (z - 0.02943725152285942?)^2]"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dop.local_basis_monomials(sing[1]) # sing[1] = α⁻¹ ≈ 0.029"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "⤷ A hyperplane of analytic solutions. Does $a(z)$ lie on it?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "mat = dop.numerical_transition_matrix([0, sing[1]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[4.5463762475228446 +/- 7.72e-17]*I"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a0 = 1; a1 = 5\n",
    "c1 = a0*mat[1,2] + a1*mat[1,3]\n",
    "c1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "⤷ Thus $a(z) \\sim 4.54\\dots \\sqrt{z - \\alpha^{-1}}$.\n",
    "What about $b(z)$?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1.2020569031595943 +/- 5.17e-17]"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b0 = 0; b1 = 6\n",
    "d1 = b0*mat[1,2] + b1*mat[1,3]\n",
    "d1/c1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.20205690315959"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "zeta(3.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Compacted binary trees of bounded 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": 45,
   "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": 45,
     "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": 46,
   "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": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dop = cbt.dop[5]; dop"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Singular points:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.2928932188134525?, 0.50000000000000000?, 1.707106781186548?]"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sing = dop.leading_coefficient().roots(QQbar, multiplicities=False)\n",
    "sing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.292893218813452"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "1/(4*(cos(pi/8)^2)).n()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Local behaviors at the dominant singularity:"
   ]
  },
  {
   "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 = sing[0]\n",
    "dop.local_basis_monomials(s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Singularity analysis implies\n",
    "\n",
    "$$\\#\\{\\text{cbt of rh ≤ 5}\\} \\sim \\kappa_5 \\, n! α^n n^β, \\qquad α = 4 \\cos²(\\pi/8), \\quad β=-(α+5)/8$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1, z, z^2, z^3, z^4, z^5]"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dop.local_basis_monomials(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[31.42732179257817 +/- 3.99e-15] + [27.45408401646515 +/- 4.94e-15]*I"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ini = list(cbt.egf[:6])\n",
    "# 2 = index of the singular element of the local basis\n",
    "c = (dop.numerical_transition_matrix([0,s])*vector(ini))[2]\n",
    "c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1.6248570260793 +/- 3.26e-14]"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expo = dop.local_basis_monomials(s)[2].op[1]\n",
    "(CBF(-s)^expo*c/CBF(-expo).gamma()).real()"
   ]
  },
  {
   "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": "code",
   "execution_count": 53,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(494, 9960)"
      ]
     },
     "execution_count": 53,
     "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": 54,
   "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": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Pol.<x> = QQ[]\n",
    "Dop.<Dx> = OreAlgebra(Pol)\n",
    "t0 = 2^(-sz)\n",
    "L = Dx * (x*(x-1)*(x-t0)) * Dx + x\n",
    "L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "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": 56,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "9.556619453472961320858516937971887816879016410846290172126354839510701168391942496813468949482302488162252696886520244471966592218665291723177173168199392132406276973789705340073690947932553910124796231579916084856421141195367398954924504111462358438957864101423930407260331976176725065395140834421424427590908128799409817584336622391673407243198108939261965265090478623967959224212373223620909718920076884504325493512957093231290793326009149024368397362080608773817261344914045941037215573589929466261952514844246399651185700946374898667516902390725364103154432169875542859847614933659911701076729537808003318825249648022560871556004599073286958174424349446361789734983296900487540550803091366045082974239981262926261619181189568107287215173880557039860815668766673099266893087110588547841898487650663644043980236519485686759761011681600245421501121353018152348931578768245258254780107760973263563775281863663012073352722917223246017584235308522688246510732993890785886249250725937793549943881255036341563534916875647477037886097657249792008428528870297161215131300381890686100703684816317908759813352398646880833758301441006787007561756610067089326048633636601655581680294114920232193088316072162202643820964032906843841817185410145450874756553779319769681938557795157701492901383950831572956075682891800333143265314437926822965689897336891281115810393634163501579557476432032208140913208912647082900790120127010816341829130398447859634549497723253118914899666245483559245438168221086949951270207323864052877486421440537213113570849672571099832804454703454119688641196701003291562809385222820524421206480989940681005905092328271286077626602437533505092641703146035760742921680842263781843854873991144874044147403548564686612183559041672386026367636029777757668704070453267111227724964978924932012765810717992589937434808382663799797773043364171093373260596575966020347101598724092360916823714857848399257003615082945323894832711059841119689648486033169179927891587569070063188570258161818486202606231937307760766826338422121070834991233113786550649092555451055378950142673684858167685590606256915969226571204633960531410689676714870414271577763534727803754179156714701688996696004086261971352341976630738528832473896535081911494009668569113201638788353744280283060834315636625364415819665437542451448080606771694562471371975196008988383248568008014743171957640899885078174330701068772917303994305295373221848857890630191210892023261909597212215222580922194025813235619871426751892331277173873023952664284270669470966003915290659729911588342236100212184270452170452531046709945743782625626814519771418724775958330968324506098946055581698099590764388244028515969168881038955796413001853465395680656243512452699320906509054087659901425563057011885470050958906345477915933201856109902467292184386952372661834542279554976803521112299906835150188799637765061595534103641294316535688856568007252802956390958368379318710244523774844628336072323710306061629273180267967815264454788464436467987614675704690957518274321357948708194043185e-299"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tr = mat.trace().real()\n",
    "Pol.<la> = 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": 57,
   "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": 57,
     "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": 58,
   "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": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "den_ratios"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAHWCAYAAAD6oMSKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi40LCBodHRwOi8vbWF0cGxvdGxpYi5vcmcv7US4rQAAIABJREFUeJzt3Xd4VFX+x/FPCiSUEHoVBNZGjBUCiiLsT8S1y6ogKghrYxdUxC66uhaiyCqWWMCGIogisKyyCDZAQUUQQbCi0kMnCS11fn8cb4YoJWVmzi3v1/PM870ZkplPMiT55txz7okLhUIhAQAAwPPibQcAAABAZNDYAQAA+ASNHQAAgE/Q2AEAAPgEjR0AAIBP0NgBAAD4BI0dAACAT9DYAQAA+ASNHQAAgE/Q2AEAAPgEjR0AAIBP0NgB8L1QKKTc3FyxNTYAv6OxA+B7eXl5Sk1NVV5enu0oABBVNHYAAAA+QWMHAADgEzR2AHwrKytLaWlpysjIsB0FAGIiLsRsYgA+l5ubq9TUVOXk5KhOnTq24wBA1DBiBwAA4BM0dgAAAD5BYwcAAOATNHZAQIRC0hdfSLfcIrVqJcXFhW8nnyxNn27eBwDgXSyeAHxu6lSpZ8+Kfcy0adJ550Unjw0sngAQFIzYAT71xRdmNG7vpu7qq6WZM6WCAjM659xWrZIGDw6/3/nnm499663Y5wYAVB4jdoDPFBdLHTpIixebt7t2lf77XyklpXwfn5MjdekiLV0avu+bb6Sjj4581lhhxA5AUDBiB/jIhg1SYmK4qfv5Z+njj8vf1ElSaqq0ZImUlyc51/VNT5cGDIh4XABAhNHYAT7x5ZdS06bm+JprzCnWNm0q/3i1a5vTuc7I3SuvSDVqSCUlVY4KAIgSGjvAB5YtC4+uTZggjR4ducdOTzendxMTpT17pIQEc7oWAOA+NHaAx61YYZovyVyy5NJLI/8c8fFSYaF04YXm7bp1peXLI/88AICqobEDPGzdOumww8zxxInSWWdF9/mmTJFGjTLHRx9t5u8BANyDxg7wqOJiqUULczx6tNSrV2ye98YbpdmzzfGf/2zm4QEA3IHGDvComjVN/cc/zGKJWDrtNGnuXHPcqZOZ4wcAsI/GDvCgW281FxmOi5OysuxkOPVU6d13zXF6urRmjZ0cAIAwGjvAY1aulEaONMcFBXaznH22NH68OW7Z0n4eAAg6GjvAY1q3NvXTT80lSGzr00e6+WZznJRkNwsABB2NHeAhgwaZ2rmzubnFyJHS4Yeb41NOsZsFAIKMvWIBj/jlF6ltW3NcUmLm17lJKGSudydJDz4oDRtmN8/e2CsWQFAwYgd4hNPUzZ/vvqZOMpl27jTHd98d3q8WABA7NHaABzgLFFq3lk46yWqUA6pZM7y37AknSEVFdvNkZWUpLS1NGc5+awDgc5yKBVxu71OcO3ZItWrZzVMed94pPfywlJoqbd9uOw2nYgEEByN2gMsNHGjqtdd6o6mTpMxMU3NypMces5sFAIKEETvAxXbulGrXNsduXDBxILt3h3fH+PlnqU0be1kYsQMQFIzYAS7mXNLk9de91dRJUo0a0pw55rhtW3NKGQAQXTR2gEvl5UlLlpjjyy6zm6WyunQxFzCWpB497GYBgCCgsQNcqmtXU99+226OqnJW9L7/vrRwod0sAOB3zLEDXGjHDiklxRz74Tt0xQrpsMPMsY3Phzl2AIKCETvAhf78Z1Pfestujkj505+kv/7VHPfvbzUKAPgaI3aAy+zaFb6siZ++O/e+Ht/KlVKrVrF7bkbsAAQFI3aAy/TqZeq4cXZzRFpcnPTpp+b40EPtZgEAv6KxA1wkFJLefdccX3653SzR0LmzdMQR5viRR+xmAQA/orEDXOShh0y96Sa7OaLpm29MveMOKTfXbhYA8Bvm2AEu4lyEuLg4PB/NjyZOlC691BzH4icQc+wABIWPf3UA3jJjhqmdO/u7qZOk3r3Dx5Mm2csBAH7DiB3gEs5o3bZtUt26drPEQk5O+POM9j64jNgBCAqfjwsA3pCTEz4OQlMnSamp4bmEF15oNwsA+AWNHeACV1xh6rRpdnPE2mOPmTptmrRhg90sAOAHnIoFXMA5DRnE78aPPw7vtBGtz59TsQCCghE7wLK33za1Z0+7OWzp1i18/P771mIAgC8wYgdY5ozW5eZKKSl2s9iyaZPUuLE5jsZPJEbsAAQFI3aARVu3ho+D2tRJUqNG0vnnm+N77rGbBQC8jMYOsKhvX1P/9z+7OdxgyhRTH3xQKiyMzGNmZWUpLS1NGRkZkXlAAHA5TsUCFgV50cS+jBplLoFyyinSJ59E7nE5FQsgKBixAyxZtMhUBpPChgwx9dNPpfXr7WYBAC+isQMsOfdcU998024Ot5kzx9Tmze3mAAAvorEDLHFGpFq3thrDdbp0CR9/9ZW9HADgRTR2gAWvvWbqtdfazeFWv/xi6okn2s0BAF7D4gnAAmfRxM6dUs2adrO41SGHSGvXmtWyVd1LlsUTAIKCxg6IsaIiqVo1c8x33/7l5kqpqea4ql8nGjsAQcGpWCDGHnnE1JEj7eZwuzp1pLPOMsePPmo3CwB4BSN2QIw5p2GLiqSEBLtZ3G7v0c2SkvDXrqIYsQMQFIzYATGUnx8+pqk7uMRE6fbbzfGVV9rNAgBeQGMHxND995v69NN2c3jJww+b+tprZRtjAMAf0dgBMTR8uKl//7vdHF7z/POmnnmm3RwA4HbMsQNipKBASkoyx3zXVZwzvy4nxyysqAjm2AEICkbsgBhxVsE+/rjdHF41ZYqp7dvbzQEAbsaIHRAjrIatOudrmJ0tNWlS/o9jxA5AUDBiB8TA3n8+0dRV3pw5prZtazcHALgVjR0QA87esDffbDeH13XpYuquXdJPP9nNAgBuRGMHxIBzDbZ//ctuDrd45pln1KZNGyUnJ6t9+/aaO3duuT7ujTfekBQn6UIdfnhUIwKAJ9HYATFUq5btBPZNnDhRQ4YM0bBhw/TVV1+pS5cuOuuss7Rq1aoDftzKlSt1yy23qIszbCdpyZJopwUAb6GxA6Ls/fdN7d3bbg63eOyxx3TVVVfp6quvVrt27TRq1Ci1bNlSzz777H4/pri4WJdffrn+9a9/qW3bture3dx/3HExCg0AHkFjB0RZ//6mjhplNYYrFBQUaOHCherRo0eZ+3v06KF58+bt9+Puv/9+NWrUSFdddZUkM/LpXBPws8+iFhcAPCfRdgDA79auNbVpU7s53GDz5s0qLi5Wk99dq6RJkybKzs7e58d8+umnevHFF7V48eIy9//8s9SihXTyyX+84HN+fr7y99p/LDc3NzKfAAC4HCN2QBQ5Kzc7drSbw23inAvS/SYUCv3hPknKy8vTFVdcoTFjxqhhw4Zl/q158/C17D78sOzHZWZmKjU1tfTWsmXLiOYHALdixA6IoltvNfXJJ+3mcIuGDRsqISHhD6NzGzdu/MMoniStWLFCv/76q84777zS+0pKSiRJiYmJmj//e3Xs+CedfnrZUbs777xTQ4cOLX07NzeX5g5AIDBiB0TR1KmmdupkN4dbVK9eXe3bt9esWbPK3D9r1ix17tz5D+9/1FFHaenSpVq8eHHp7fzzz9ef//xnLV68WMcd11Jt2pj3nT49/HFJSUmqU6dOmRsABAEjdkCUFBfbTuBOQ4cOVd++fdWhQwedfPLJGj16tFatWqWBAwdKkvr166cWLVooMzNTycnJSk9PL/PxdevWlaTS+xctkurVk845549z7QAgaGjsgCh56SVT77vPagzX6d27t7Zs2aL7779f69evV3p6uqZPn65DDz1UkrRq1SrFx5f/ZELdulJamrR8ufSf/0gXXBCt5ADgfnGhEH/jAtFQo4a0Z4+0c6dUs6btNP6Wmyulpprjff1Ey83NVWpqqnJycjgtC8DXmGMHRMmePabS1EVfnTrSCSeY47fespsFAGxyzanYJ56QVq2S2rSR+vUzP6gBr1qxwlQucxI7c+dKtWtLvXox1w5AcFkdsQuFpGHDpLg4acgQ6bHHpOuvN6dUTjqJyefwrttuM/WJJ+zmCJJatczPDUkaP95uFgCwxdocu1BI2nt+9NVXm62XFiyQbropfP+uXWauEuAlzrV2GTmKrV27TIMnlf3aM8cOQFBYG7GrV8/Udu2kkhJpzBjplFPMyF0oJPXsaf69Zk2psNBWSqDiGGm2p2ZNqWtXc/zKK1ajAIAVVhq7s86ScnLM8fLl4dGNvU2eLF13nTl2NvsGvMBpKP75T6sxAmvGDFMHDLCbAwBsiHljN3Fi+AfvbzsD7ddzz0mJiWYE7+67o58NiIQbbjDVmWeH2EpOlv7v/8zxCy/YzQIAsRbTOXZ7z6vLzZVSUg7+MYWFUvXq5jgnh9WycD/m19mXn28aPMm8DsyxAxAUMR2xu/xyU6+/vnxNnSRVqyZNmmSOW7SITi4gUn75xdQOHezmCLqkJKlHD3M8erTdLAAQSzEbsdu6VWrQwBxX5hmdUZBFi8IXIgXcpndv6c03pU8+MYuBYM/eo3Y5OYzYAQiGmI3YnXaaqe++W7mPd0ZCTjwxMnmAaHjzTVNp6uxLSpL+8hdz/OKLdrMAQKzEZMSuqMicUpWqNu/IGbX7/nvpiCOqnguIpL3nkDK/zh0KCpxV9bmSGLED4H8xGbFzLlvy5JNVe5zly0098siqPQ4QDR98YOpVV9nNgbDq1aVzzrGdAgBiJyYjdpFcJeg81qZNUsOGVX88IFI6d5bmz5dWrpRatbKdBg4zaseIHYBgiPqI3fPPm+qsiK2qWbNMveSSyDweECnz55tKU+cu1auH59qxQhaA30V9xO6BB8wV+PPzw9ejqyquEwa3YX6du23ZkquGDVMl5SgUYsQOgH9FfcTunnvML7pINXWSdNllpo4bF7nHBKpi+nRTBw+2mwNlZWVlKS0tTZ07Z5Te99RTFgMBQJTFdOeJSNm1S6pVyxx7Lz38qH17c43FNWu4kLYbOTtPSDmS6vBzA4BvxXyv2EioWTN8nJNjLwfgWLTIVJo6dzvvPFMff9xuDgCIFk82dpI0ebKpQ4bYzQEw+uMdr7xi6tChVmMAQNR4trHr2dNU5wc1YMvUqabedJPdHDi4xESpVy9z/O9/280CANHgyTl2jjp1pLw8acMGqXFj22kQVMccI33zjZSdLTVpYjsN9sWZY5eTk6OaNetEZCccAHAjz47YSeHTsQMH2s2BYPvmG1Np6rwhMVHq08ccjxhhNwsARJqnR+wkrmkHu0pKpIQEc8z/Qffae8SuTp06Ki42DZ7E6wbAXzw9YidJzZubunGj3RwIprfeMvX22+3mQMUkJEh9+5rjhx6ymwUAIsnzI3bTp5tNvq++WhozxnYaBM2RR0o//MDexW73+xE7SYzaAfAlz4/YnX22qS+8YDcHgumHH0ylqfOehATpqqvM8T//aTcLAESK50fspPA8u+Li8H6dQLQVFYnVlR6xrxE7qewcyZKS8M8SAPAqX7RBd9xh6muv2c2BYJkwwdS777abA5UXHx/e3/fWW+1mAYBI8MWI3Y4dUkqK2Wps507baRAUbdpIv/4qbd0q1atnOw0OZH8jdpIZbXVG+hm1A+B1vhixq13b1F277OZAsPz6q6k0dd4WFyfddps5dkbvAMCrfDFiJ0ldukiffCItXy61a2c7DfyO69d5y4FG7KSyo3bM1QXgZb758fXUU6YOGWI3B4LB2R+WeVn+EBcn3XuvOb76artZAKAqfDNiJ7ELBWInI0P68ktp3TqpWTPbaXAwBxuxczg/Q4qKwiOyAOAlvhmx21thoe0E8LsvvzSVps5fHn7Y1Msvt5sDACrLV43dgw+a+sorVmMA8Chna7iJE/kDEYA3+aqxu/FGU++8024O+NsXX5h62WV2cyA6nnzS1J497eYAgMrw1Rw7iXl2iL6ePc3iiSVLpGOOsZ0G5VHeOXYO5+fIrl1SjRpRDgcAEeSrETtJatrU1NxcuzngX86KWJo698vKylJaWpoyMjIq9HHOdI7TT498JgCIJt+N2L38svS3v0n33Re+fAEQSYwKe09FR+yk8Oucm2t2tgEAL/DdiF3fvqbef7/dHPCnlStNPe00uzkQfZMnm1rBwT4AsMp3jV1ioqklJXZzwJ9GjTLV2YIK/uUsnvj+e2nTJrtZAKC8fNfYSeFri+3YYTcH/Mdp7M46y24OxMYHH5h61FF2cwBAefmysbvnHlOfe85uDvgXe4kGw//9n6lbt4ZPwwOAm/lu8YQk7dljLlGQmipt3247DfwiP19KTpbq1TO/6OEdlVk84Vi4UOrQwRz776clAL/x5bhDcrKpOTl2c8BfJkww9ZZb7OZAbLVvHz7++mt7OQCgPHzZ2ElS7dqm5ufbzQH/GDnS1IED7eZA7K1YYerxx9vNAQAH49vG7q67TB0/3m4O+MeyZabWr283B2KvbVupbl1z7CyoAAA38uUcO0nats38Am7XTlq+3HYa+AEXJvauqsyxc2zeLDVqZI75PwDArXw7Ylevnqnffms3B/zhq69M7dXLbg7Y07ChdNxx5pgzAQDcyreNHRBJ//63qTffbDcH7Jo3z9TLL7ebAwD2x9eNXZ8+pi5YYDcHvO/1103t2NFuDthVs6Z0/vnmeMQIu1kAYF98O8dOMqfPTjzRnD6bONF2GngZ8+u8LRJz7BzFxWW3LnT+bwCAG/h6xO6EE0x98027OeBtztZ0bdrYzQF3SEgIX8vwqqvsZgGA3/N1YwdEwquvmjpkiN0ccI9HHzX15ZfNTjcA4Ba+b+yczbtzc+3mgHeNGmXqgAF2c8BdXnjB1FNPtZsDAPbm+8bu+utNHTvWbg54148/mpqSYjcH3MU5DbtwoZSdbTcLADh839j17Wvqk0/azQHAf+bMMbVZM7s5AMDh+8bOGWX56Se7OeBNzjZiPXvazQF36tIlfLxwob0cAODwfWMHVMUTT5jKwgnsz6+/mtqhg9UYACApII3dueea+v33dnPAe8aMMXXvkRlgb4ceKrVta465tBIA2wLR2N1wg6lPP203B7yLi9B6U1ZWltLS0pSRkRHV5/n6a1N7947q0wDAQfl65wlHKCTFx4ePgfIoKpKqVZNq15by8mynQVVEcueJ/enTR3rjDXPxYuc6dwAQa4EYsWO0BZXxzjumDhpkNwe8Yfx4U0eOlPLz7WYBEFyBaOz2xogdyss5dT9woN0c8Ia4OGn0aHPcsaPdLACCKzCN3TXXmDpvnt0c8I4PPjC1dWurMeAhzs+ZJUuklSvtZgEQTIFp7K67ztTnn7ebA4C/ffmlqfxBAMCGwDR27dub+tprdnPAG5zFEocfbjcHvKd9eykpyRy/957dLACCJzCNHVAR48aZOniw3RzwpnXrTP3LX+zmABA8gWrsEhJMZQEFDiYry9Qrr7SbA95Uv354G7rbbrObBUCwBKqxc1Y3fvSR3RxwP2eP2NRUuzngXW+/beqjj0q7dtnNAiA4AtXYOQsonnvObg4A/hcXJ73yijlu185qFAABEoidJ/bmXKw4WJ81KmLtWumQQ6TTTpNmz7adBpEQi50n9sf5mbN0qZSeHtOnBhBAgRqxA8rDuSQOO04gEr7/3tRjjrGbA0AwBLaxY8QO++MsnHAmvwNVccQR4abu8cftZgHgf4Fr7AYMMPWLL+zmgHtt3WpqtWp2c8A/nIsWDx0qFRTYzQLA3wLX2F11lakvvWQ3B4DgqF5deuopc3zccXazAPC3wC2eCIWk+HgpMVEqLLSdBm7z9dfS8cdLvXtLb7xhOw0ixebiib05Cym+/lo69lhrMQD4WOBG7JwfrEVFdnPAnZz5dew4EV3PPPOM2rRpo+TkZLVv315z587d7/uOGTNGXbp0Ub169VSvXj11795dX3h0LoWzkIJROwDRErjGDjiQMWNMPeUUuzn8bOLEiRoyZIiGDRumr776Sl26dNFZZ52lVatW7fP9P/74Y/Xp00cfffSR5s+fr1atWqlHjx5au3ZtjJNX3RFHSJ06meP77rMaBYBPBe5UrCSdd570zjvSjz9Khx1mOw3chOscRl+nTp104okn6tlnny29r127drrwwguVmZl50I8vLi5WvXr19PTTT6tfv37lek63nIqVpOJiMxVEkrZvZ3cTAJEVyBG7v/3N1JdftpsD7lJSYqrzSxeRV1BQoIULF6pHjx5l7u/Ro4fmzZtXrsfYtWuXCgsLVb9+/WhEjLqEBOnNN81x3bp2swDwn0A2dueeayorY7G3Dz809e9/t5vDzzZv3qzi4mI1adKkzP1NmjRRdnZ2uR7jjjvuUIsWLdS9e/f9vk9+fr5yc3PL3NzkkkvCx5Mn28sBwH8C2dg51ycr5+8RBISz44SzpzCiJ8455/2bUCj0h/v2ZcSIEZowYYImT56s5OTk/b5fZmamUlNTS28tW7ascuZI27bN1IsuMqdnASASAtnYAfsyaZKpRx9tN4efNWzYUAkJCX8Yndu4ceMfRvF+b+TIkRo+fLhmzpypYw9yrZA777xTOTk5pbfVq1dXOXuk1a0r3X23Oe7SxW4WAP4R2MbOWZm2ZYvdHECQVK9eXe3bt9esWbPK3D9r1ix17tx5vx/36KOP6oEHHtCMGTPUoUOHgz5PUlKS6tSpU+bmRg88YOr8+dI339jNAsAfAtvYXXGFqW+9ZTcH3ME5FVarlt0cQTB06FC98MILeumll/Ttt9/qpptu0qpVqzRw4EBJUr9+/XTnnXeWvv+IESN0991366WXXlLr1q2VnZ2t7Oxs7dixw9anEFHffWeqs58sAFRFYBu73r1NHTfObg64w3vvmfpbb4Eo6t27t0aNGqX7779fxx9/vObMmaPp06fr0EMPlSStWrVK69evL33/Z555RgUFBbr44ovVrFmz0tvIkSNtfQoRdeSR4QVdzop9AKisQF7HzsE1y+A4/3zpv/81OwMccYTtNIg0N13Hbl+crQ4l6ZdfpNatrcYB4GE0dqKxA/8X/M7tjZ0kffFFeO4v/w8BVFZgT8VK4QvR8kMUgG0dO0onnWSOb7zRbhYA3hXoxq5vX1O//NJuDtjlLJxgayfY9umnpj75pDklCwAVFejGzlkZywKKYHMWTlx9td0cQHy8ufSJJLVtazcLAG8KdGPXtaupNHbB9sILptLYwQ1OOknq1s0c838SQEUFevGExKR58H8gCLyweGJve6+SXbZMSkuzmweAdwR6xA4A3CguTlq82BwffTR/dAAov8A3ds7V3vPy7OaAHc4vzGrV7OYAfu+446R+/cyxM20EAA4m8I2ds4BiyhS7OWDH7NmmMpcJbjR2rKlz50ozZ9rNAsAbAt/Y9eljKgsogslZOHHNNXZzAPuzdq2pZ54p7d5tNwsA9wv84gmJyfNBxmsfDF5bPPF7Tz4Zvmgx/1cBHEjgR+wAwO1uuEGqVcsc33+/3SwA3I3GDoHFyAe8ZNs2U++9V/r2W7tZALgXjZ2kSy4xddkyuzkQWwsWmHrllXZzAOVRrVr4/2xamlRUZDcPAHeisVN4Zezrr9vNgdhi4QS8pkMH6eabzXH9+nazAHAnFk9IKiiQkpKkli2lVatsp0GsOAsnSkrCx/Anry+e+D3n/+tDD0l33WU3CwB3YcROUvXqpq5ebTcH7KCp86+srCylpaUpIyPDdpSI2rPH1GHDpOXL7WYB4C6M2P2Gy14ED695cPhtxE6SvvxScvrVggJ2TwFgMGL3mzZtTM3Pt5sDsbF0qam9etnNAVRWhw7SHXeYY+esAwDQ2P3m0ktNfe89uzkQGy++aCpbicHLMjOlhg3Ncd++drMAcAcau984jd3EiXZzIDbGjDH19NPt5gCqasMGU8eNkyZNspsFgH3MsftNKCTFx0sJCVwfKgiYXxcsfpxjt7fsbKlZM3O8cqXUqpXdPADsYcTuN84v+uJiuzkAoKKaNpWmTzfHhx4qFRbazQPAHho7BM5PP5l67rl2cwCRdNZZ0uDB5pjFFEBw0djt5YQTTM3NtZsD0fXSS6ay4wT85qmnzOidJPXoYTcLADto7PbiLKCYNs1uDkSXs3Di7LPt5gCiYd06U2fNMqtmAQQLjd1enGuavfGG3RyIrs2bTU1MtJsDiIa4OCkvzxzfdZf07rt28wCILVbF/g6rJf2P1zh4/L4qdl9++kk6/HBzvGSJdMwxdvMAiA1G7BAo69eb2rWr3RxAtB12mPTRR+b42GOljRvt5gEQGzR2CJRx40wdMMBuDiAWunUL77LSpIm0Z4/VOABigMbud/78Z1P569afnBWxF19sNwcQK3/7m3Tzzea4Rg2mIAB+R2P3O87K2LfespsD0fHdd6bWqmU3BxBLI0eG/2iN56c+4Gt8i/9Oz56m0tgB8JMPP5SaNzfHDRrYzQIgemjsfqdRI1Nnz7abA5G3a5epRx5pNwdgy9q1pm7dKp14ot0sAKKDxg6B8fbbprJwAkFWUmLqV1+xOwXgRzR2CAxn4US/fnZzADbFxUlFReZ41iyzxywA/6Cx24fu3U11rnkGf/j4Y1ObNbMaA7AuIUEqLDTHM2ZI551nNw+AyKGx2wdnazHn1B0A+E1iYri5e+cd6YIL7OYBEBk0dvvAylj/KS42NSC7SeE3WVlZSktLU0ZGhu0orpSYKBUUmONp06Qzz7SbB0DVsVfsfrCfqL/MmmUmig8ZIj3+uO00iLUg7hVbEYWFUvXq5vjEE6WFC+3mAVB5jNghEF5+2dT+/a3GAFypWrXwgopFi6R69ezmAVB5jNjtByN2/sLrGWyM2JVPKFR2Z4riYnaqALyGb9n9cK7vtG6d3RwAECtxcaa5a9vWvJ2QEL6wNwBvoLHbj0suMZWVsQCCZsUK6dxzzXGtWlz6CfASGrv9uPBCU998024OVN1335l60UV2cwBe8t//SsOGmePmzaU5c+zmAVA+NHb70bChqZ98YjcHqu6VV0xl4QRQMQ8+KE2caI67dpVGjLCbB8DBsXjiAJhw7w9Nm0obNpjrdVWrZjsNbGDxRNV8953Urp057to1vIsLAPehsTsAGjt/4HUEjV3V7dyjXLBBAAAgAElEQVQp1a4dfjs/P3ztOwDuwanYAzjjDFPXrrWbAwBsq1VLKimR/vQn83ZSkvTrr1YjAdgHGrsDuPhiU6dMsZsDlbdtm6nsKAVUXVyc9NNP0tCh5u02bcJzWAG4A43dATgrYydNspsDlTdhgqksnAAi59//lj74wBwPGGC2IWOqA+AOzLE7COZneVvHjtKCBdKWLVL9+rbTwBbm2EXH9u1ltx9buVJq1cpeHgCM2MHnFiwwlaYOiLy6dc0fvb17m7cPPVTKzLSbCQg6RuwOghE7b+P1g8SIXSzMni116xZ+OzdXSkmxFgcILEbsDqJLF1M3bbKbAxVXVGRqo0Z2cwBB0LWruQRKrVrm7Tp1pBdesJsJCCIau4NwVsb+5z92c6DiZs0ylYUTQGxUry7t2CGNHWvevuYaM2qem2s3FxAkNHYH8de/mvr223ZzoOKcXy5XXmk3BxA0/fqZBs+ZCpGaKt17r91MQFAwx64cmKflTbxucDDHzp433wwvrpCkZcuktDR7eQC/Y8QOABA1vXqZ+a6dO5u3jz5aOuEEs3czgMijsQMARFVCgvTpp2a0TpIWLzZbknF6Fog8GrtycLajysmxmwPl9/33pl50kd0cAMLS0szUiKeeMm/ff7+ZMjF5st1cgJ/Q2JWDszL2nXfs5kD5OQsnWBEbbFlZWUpLS1MGmwW7yuDBUnGx1LOnefuii0yD98kndnMBfsDiiXJYsUI67DCzd+yUKbbToDyaN5fWrzfzeKpVs50GtrF4wr3y8sycuxUrwvfNnSudeqq9TICX0diVEyssvYXXC3ujsXO/9eul9HRp69bwfbNnS6edZi8T4EWcigUAWNesmbRli5SdLTVsaO7r2tX8keZMrQBwcDR28B1nkcvxx9vNAaDimjQxWzhu3Ci1aWPu69/fNHhDhpi5eQD2j8aunI46ytQ9e+zmwMFNmmQqO04A3tWokfTzz9Lu3eFFFk88ISUmmjm0zsp3AGXR2JWTc9kMZ/9RuJdz2uayy+zmAFB1ycnmciglJVJmprlv/Xrzx3ZcnPTAA+bfABg0duXkNHbsGet+c+ea2rix3RwAIicuTrrjDrMgavlyqWlTc/8//2kugJyaKs2bZzcj4AY0duXkzNeisQMAu9q1M6N2xcWmsZOk3FzplFNMA9i+PadqEVw0duXkXD5jxw67OXBgzimZWrXs5gAQffHx0r/+ZUbx1q6VunUz9y9aFD5V2769tGSJ1ZhATNHYwVecK9ezcAIIlubNpY8+Mk3e0qXmmniSafKOO840eXXrSq++ypw8+BuNXQU0b24qy+3dy1k4QWMHBFd6umnuQiHpu+/CFznOyTE/GxISTKN3ySVld7wA/IDGrgKcBRTsZ+heTmPH1qAAJOnII80OFqGQmYd3113hf5s0yWwXGRdnbgMHSitX2ssKRAKNXQWwMtb9nNFUZ04kADhSUqSHHjJNXigkffqp2d3C8fzzUuvW4Uave3fp3Xc5dQtvYa/YCiguDl8cc+1a22mwL+wRi31hr1gcTChkrlM6fLgZ4dufiy6Srr3WNH3xDI3AhfhvWQEJCaauW2c3B/ZtzRpTu3e3mwOA98TFST16SB9/HB7R+/lnaehQqXr18Pu9/bZ05pnheXpxceZyWA8+KH37rbX4QCkaO/jGuHGmsnACQCS0aSP9+99Sfn642fvmG+nGG6XatcPv9/XX0j33SGlp4WYvLk5q0MDsgPPKK5zlQexwKraCUlLMtexKSpjH5Tbt2pkVcDt2cB07lMWpWETTjh3S//5nRvMmTSr/lROOOMJcVPmUU6TOnc1CD07voqpo7Cqof3+z8nLRIumEE2ynwd6YX4f9obGDDaGQ2QFj5kxze+89qaioYo8RFycde6y5paebUcF27cwiD2d6ELA3GrsK+u9/pfPPl+6+22w+DfegscP+0NjBjYqLzandTz6RPv/c3H74ITKP3bSpaf4OOcTcWrQwt0MOMf/WrJk5s8GZJ/+hsaugPXukGjXMX0zLl9tOA0d+vpScLB16qPTrr7bTwG1o7OBlRUXSjz9Ky5aZ3zvLl5uRwJ9+iv02l4mJZgePOnWk1FRTU1LCtXbt8K1mTXOrVcvcatQI35KTw7VaNfMzfNcuM80pMdGMRjrVucXH//E4Pp7m9PcSK/NBoVBIeXl5kc7iKd9+ay52CXf4z39MvfRSXhdI+fn5ys/PL33b+XmVy38OeJQz4tajR/k/prBQys6W1q83V3NYv94s4tiwwdy/caO5bd1a/scsKpI2bzY3VNz06WZOZUWkpKQorgLda6VG7Jy/fgEAABA9FT3TUKnGrqIjdhkZGVqwYEFFn6ZKovmcTk+bkxO+Lzc3Vy1bttTq1atjeqrHb1/byj7nvl6TquD19PZz/n7Ebv369erYsaOWL1+uFi1aRO15f8+PX1s3PCffn/56Tl7PA6voiF2lTsXGxcVV6IufkJAQ83kt0XzOnj2lKVOkTZukP/2p7L/VqVMnpp+r3762VX3OSMfi9fTPc0rmBySvpz+eU+L700/PKfF6RkpMrpgzaNCgWDxNzJ7T2TN28uSoPUW5+e1r66bntCEoX1teT57Ti4LyteX19PZzsiq2EnJyzKqgTp2kzz4z97Hqzp6vvzZb+vTpI40fH5nH5PX0lzVr1pSe6jnkkENsx0EV8f3pL7yekcU1rivBmc/1+efh+5KSknTvvfcqKSnJTqgAe/VVUyO5lRivp784ryOvpz/w/ekvvJ6RxYhdJXExXPdo2FDassUs7U+s1KxR+B0jAgCCghE7eN6WLabS1AEAgo7GrpJOP93U7Gy7OQAAABw0dpXkrIydOtVujqDbts3U9u3t5gAAwA1o7CqpZ09T337bbo6ge/NNUyO5cAIAAK+isaukpk1Nff/9P/7bddddp7i4OI0aNSq2oQJo7FhT+/SJ3GMWFhbq9ttv1zHHHKNatWqpefPm6tevn9atWxe5JwFQYZmZmcrIyFBKSooaN26sCy+8UN9//73tWIiAzMxMxcXFaciQIbajeB6NXYRNnTpVn3/+uZo3b247SiDMn29qw4aRe8xdu3Zp0aJFuueee7Ro0SJNnjxZP/zwg84///zIPQmACps9e7YGDRqkzz77TLNmzVJRUZF69OihnTt32o6GKliwYIFGjx6tY4891nYUX2AdYQStXbtWgwcP1nvvvadzzjnHdhxUUmpqqmbNmlXmvqeeekodO3bUqlWr1KpVK0vJgGCbMWNGmbdffvllNW7cWAsXLtRpp51mKRWqYseOHbr88ss1ZswYPfjgg7bj+AIjdlXQsaOpOTlSSUmJ+vbtq1tvvVVHH3203WABUVJiakpK9J8rJydHcXFxqlu3bvSfDEC55OTkSJLq169vOQkqa9CgQTrnnHPUvXt321F8gxG7KrjoIumLL6R33pFWrXpEiYmJuuGGG2zHCoyPPza1f//oPs+ePXt0xx136LLLLuPitoBLhEIhDR06VKeeeqrS09Ntx0ElvPHGG1q0aJEWLFhgO4qv0NhVwuuvv67rrruudNeJ559/Vz/88IQWLVqkOGdLCkSds3Ciqo2d83o6/ve//6lLly6SzEKKSy+9VCUlJXrmmWeq9kQAImbw4MFasmSJPvnkE9tRUAmrV6/WjTfeqJkzZyo5Odl2HF9hS7FKyMvL04YNGyRJhx8uSW8pLm6Y4uPDZ7aLi4sVHx+vli1b6tdff7WS0++cHrqkJHxcGXu/npLUokUL1ahRQ4WFherVq5d+/vlnffjhh2rQoEEVE8MWthTzl+uvv15Tp07VnDlz1KZNG9txUAlTp05Vz549lZCQUHpfcXGx4uLiFB8fr/z8/DL/hvKjsasi01Bs0dKl68vcf+aZZ6pv374aMGCAjjzySCvZ/C6a+/U6Td2PP/6ojz76SI0aNYr8kyBmaOz8IRQK6frrr9eUKVP08ccf63DzlzU8KC8vTytXrixz34ABA3TUUUfp9ttv5/R6FXAqNiIaKD297GhOtWrV1LRpU5o6DyoqKtLFF1+sRYsW6Z133lFxcbGyf9s7rn79+qpevbrlhEAwDRo0SOPHj9d//vMfpaSklH5fpqamqkaNGpbToSJSUlL+0LzVqlVLDRo0oKmrIlbFVtFRR5m6Z4/dHEHjnN0+++zIP/aaNWs0bdo0rVmzRscff7yaNWtWeps3b17knxBAuTz77LPKyclRt27dynxfTpw40XY0wDUYsauiiy6SHnpImjVLOu+88P3Mq4uuV181NRpbibVu3VrMUADch+9Lf/vYudQBqoQRuyq6+GJTJ02ymyNonBWxbAaBA8nKylJaWpoyMjJsRwGAmGDxRBWFQlJ8vFSrlrRjh+00wRHNhRPwHxZPAAgKRuyqyGkw2KoQAADYRmMHz3Ga6Hbt7OYAAMBtaOwioGVLUwsL7eYIirffNjXaW4kBAOA1NHYR4CygmD3bbo6geOUVU/v2tRoDAADXobGLAFbGxtZHH5narJndHAAAuA2rYiOgpERKSJAaNpQ2bbKdxv9YEYuKYlUsgKBgxC4C4n/7Km7ebDdHEJSUmFq7tt0cAAC4EY0dPGXOHFMHDLCbAwAAN6Kxi5BGjUx1RpQQHS+/bCorYgEA+CMauwhxFlB89pndHH7n7BF7wgl2cwAA4EY0dhFy0UWmvvWW3RxB4SygAAAAYTR2EdK1q6lc8gQAANhCYxchiYmmrlljN4ef/fijqRdcYDcHAABuRWMHzxg71lRWxAIAsG9coDiC6taVcnLMyljmgEVeixbSunVSQYFUrZrtNPASLlAMICgYsYugSy4xdcECuzn8at06U2nqAADYNxq7COrVy1RWxgIAABto7CKoWzdTaewib/t2U0880W4OAADcjMYugpxThCtX2s3hRxMnmsrCCQAA9o/GDp7gbCXWp4/dHAAAuBmNXYSlpprKWuPI+vxzUxs0sJsDAAA3o7GLsN69Tf3iC7s5AEhZWVlKS0tTRkaG7SgAEBM0dhHmrIx98027OfykqMjURo3s5oD3DBo0SMuXL9cCrkEEICBo7CLM2TPWmeyPqps501QWTgAAcGA0dhHm7Bm7dq3dHH7iLJzo399qDAAAXI8txaLA2U6Mr2xk8PVEVbGlGICgYMQuCho3NrW42G4OAAAQLDR2UeAsoPjkE7s5AABAsNDYRQErYyPnq69MveIKuzkAAPACGrsoOOUUU1kZW3UvvmjqVVfZzQEAgBeweCJKmPAfGTVqSHv2SCUl4a8pUFEsngAQFIzYwdX27DGVpg4AgIOjsYuSQw4x1dk1AQAAINpo7KLE2TP2ww/t5vCyn3829eyz7eYAAMAraOyixGnsWEBRec6OEyycAACgfFg8ESWhkBQfL9WqJe3YYTuNNzVvLq1fL+XnS9Wr204DL2PxBICgYMQuSpzJ/jt32s3hZevXm0pTBwBA+dDYAQAA+ASNXRQdfbSpjNpV3KZNpp50kt0cAAB4CY1dFPXpY+q0aXZzeNHYsaZee63dHAAAeAmNXRQ5jd348XZzeNGYMaY6++4CAICDY1VslLG1WOXwdUMksSoWQFAwYgcAAOATNHZwna1bTT3xRLs54H1ZWVlKS0tTRkaG7SgAEBM0dlF2xhmmrlljN4eXvPaaqSycQFUNGjRIy5cv14IFC2xHAYCYoLGLMmcBBVuLlZ+zcML52gEAgPJh8USU5eRIdetK7dtLX35pO403sHACkcbiCQBBwYhdlKWmmrpwod0cAADA/2js4Co5OaYec4zdHAAAeBGNXQxUq2YqpxYPzrmY8zXX2M0BAIAX0djFgLMI4Kuv7ObwgtGjTb3iCrs5AADwIhq7GLjsMlMnTLCbwwsWLza1Xj27OQAA8CJWxcZAUZE5Hdu8ubR2re007saKWEQDq2IBBAUjdjGQmGjqunV2c7jd9u2msnACAIDKobGDa7z6qql//7vdHAAAeBWNXYy0a2fqjh12c7jZc8+ZysIJAAAqh8YuRq680tQpU+zmcLNvvzU1JcVuDgAAvIrGLkacUaixY+3mAAAA/kVjFyMtWpj6wQd2c7jVhg2mnnyy3RwAAHgZjR1c4aWXTGXhBAAAlcd17GKoWjVzTbuSkvD12mC0aiWtXi3t3i0lJ9tOA7/hOnYAgoIRuxhyFlB8/rndHG60erWpNHX+FgqFdN9996l58+aqUaOGunXrpmXLlh3wYzIzM5WRkaGUlBQ1btxYF154ob7//vsYJQYAb6GxiyGnsWMBBYJqxIgReuyxx/T0009rwYIFatq0qc444wzl5eXt92Nmz56tQYMG6bPPPtOsWbNUVFSkHj16aOfOnTFMDgDewKnYGCopkRISpJo1JX4nha1cKbVuLZ1xhjRzpu00iJZQKKTmzZtryJAhuv322yVJ+fn5atKkiR555BFdd9115XqcTZs2qXHjxpo9e7ZOO+20cn0Mp2IBBAUjdjEU/9tXe9cuuzncZswYU1k44W+//PKLsrOz1aNHj9L7kpKS1LVrV82bN6/cj5OTkyNJql+//n7fJz8/X7m5uWVuABAENHawztlx4rzz7OZAdGVnZ0uSmjRpUub+Jk2alP7bwYRCIQ0dOlSnnnqq0tPT9/t+mZmZSk1NLb21bNmy8sEBwENo7GLMOXPkXLcN0pYtpiYm2s2ByHr99ddVu3bt0lthYaEkKe53S8JDodAf7tufwYMHa8mSJZowYcIB3+/OO+9UTk5O6W21szoHAHyOX6UxduWV0pw50vjx0k032U4DRM/555+vTp06lb6dn58vyYzcNWvWrPT+jRs3/mEUb1+uv/56TZs2TXPmzNEhhxxywPdNSkpSUlJSJZMDgHcxYhdjF19sKitjja+/NrV3b7s5EHkpKSk67LDDSm9paWlq2rSpZs2aVfo+BQUFmj17tjp37rzfxwmFQho8eLAmT56sDz/8UG3atIlFfADwJBq7GHMW5DkNTdA9+aSpN95oNweiLy4uTkOGDNHw4cM1ZcoUffPNN+rfv79q1qypyy67rPT9Tj/9dD399NOlbw8aNEjjxo3T+PHjlZKSouzsbGVnZ2v37t02Pg0AcDVOxcIqZyuxk06ymwOxcdttt2n37t36xz/+oW3btqlTp06aOXOmUlJSSt9nxYoV2rx5c+nbzz77rCSpW7duZR7r5ZdfVv/+/WMRGwA8g+vYWdCsmZSdLRUUmG3GgsyZM8//QkQT17EDEBScirXA2YHivffs5rDtt0WSOsDlyAAAQAXQ2FnA1mLG5Mmm3nCD3RwAAPgFp2It4RSk1LmzNH++tH691LSp7TTwM07FAggKRuxgzfz5ptLUAQAQGTR2lgV5xA4AAEQWjZ0lV19tqjNqFTTr15t6gOvSAgCACqKxs+Taa019/nm7OWx57jlTuTAxAACRw+IJS0IhKT4+fBw09epJ27dzLT/EBosnAAQFI3aWOKtig2r7dlNp6gAAiBwaOxcI2ohd0D5fAABihcbOImcBxbx5dnPEmvP5Op8/AACIDBo7i5wFFKNH280Ra088Yer119vNAQCA37B4wqKgLqBg1w3EGosnAAQFI3YWBX0BBQAAiCwaO5cIyujVjh2mtmljNwcAAH5EY2fZNdeYGpQFFC+9ZOott9jNAQCAH9HYWXbddaYGZQeKkSNN7d/fagwERFZWltLS0pSRkWE7CgDEBIsnXCBIiwmC9LnCPVg8ASAoGLFDzNDMAQAQXTR2LuBc8qSkxG6OaJs509TBg+3mAADAr2jsXOCGG0x97z27OaJt+HBTWTgBAEB0MMfOBVaskA47TDrzTGnGDNtpoof5dbCFOXYAgoLGziWC0PQE4XOEO9HYAQgKTsUiJhYtMvWSS+zmAADAz2jsXCI93dRt2+zmiJbMTFPvustuDgAA/IzGziVuvNHUF1+0myNaJk0y9fjj7eYAAMDPmGPnErt3SzVrSoccIq1ebTtN5DG/DjYxxw5AUDBi5xI1api6Zo3dHNHw66+mnnaa1RgAAPgejR2i7tFHTR02zG4OAAD8jsbORXr1MnXxYrs5Iu2ZZ0w94wy7OQAA8DsaOxdxFlA8+aTdHNHizLMDAADRweIJFwmFwvvG+uVV2bpVatBAOuoo6dtvbadBULF4AkBQMGLnIn4c0XrqKVO5fh0AANHHiJ3LtG0r/fKLlJMj+WFgISlJKiiQCgulxETbaRBUjNgBCApG7FzGWTnqLDjwuoICU2nqAACIPkbsXKagwIxy1a4t5eXZTlM1OTlS3bpSy5bSqlW20yDIGLEDEBSM2LlM9eqm7thhN0ckPP64qQ8+aDcHAABBwYidC9WuLe3cKeXnhxs9L3IWgxQVSQkJdrMg2BixAxAUjNi5kDPP7vXX7eaIFJo62JKVlaW0tDRlZGTYjgIAMcGInQs5c9MOO0z68UfbaSonO1tq1kw65hhpyRLbaRB0jNgBCApG7FwoNdXUn36ym6MqHn7YVObXAQAQO4zYuZQzP62kxJsXLvZ6fvgLI3YAgoIRO5caMsTU6dPt5qgqmjoAAGKHxs6lbr7Z1OHD7eaojK+/NvWcc+zmAAAgaDgV62LOaJfXXqG//EV67z1p6VIpPd12GoBTsQCCg8bOxbza2Hk1N/yLxg5AUHAq1sUGDDD1o4/s5qiIkhLbCQAACC4aOxe7915T77rLbo6KGD/eVC9lBgDALzgV63JeO63ZpIm0caO0fXv4enyAbZyKBRAUjNh5RHGx7QTls3GjqTR1AADEHo2dy91zj6kTJtjNUR6rVpl67LF2cwAAEFQ0di7nXM/uzjvt5iiPW281NSvLbg4AAIKKOXYe4JV5dl7JieBhjh2AoGDEzgMaNzZ1+3a7OQ6kqMh2AgAAQGPnAU8+aeoDD9jNcSBjxpj64IN2cwAAEGScivWAUEiKjw8fu5FzGnb3bik52W4W4Pc4FQsgKBix8wCnafICmjoAAOyhsfOI/v1NnTnTaox9evttUwcOtJsDAICg41SsR6xbJ7VoIR15pPTdd7bTlFW9ulRYKOXkSJzlghtxKhZAUDBi5xHNm5v6/fd2c+xLYaGp/L4EAMAuGjsPad3a1HXrrMYoY8YMU6+4wm4OYF+ysrKUlpamjIwM21EAICY4FeshH34onX661KuXNHGi7TRG/frStm3S5s1Sgwa20wD7xqlYAEFBY+cxbtvdwW15gH2hsQMQFJyK9ai8PNsJpDlzTO3Z024OAABg0Nh5zBtvmDp0qN0cknT55aZmZdnNAQAADE7FepBbTn+6JQdwMJyKBRAUjNh5UK1apm7aZC/De++ZesEF9jIAAICyaOw8aMoUU6+7zl6Gv/zF1JdftpcBAACUxalYj7J5GrS4WEpMtPf8QEVxKhZAUDBi51GtWpm6cmXsn3v4cFMzM2P/3AAAYP8YsfOoL7+UMjKk006TZs+O7XM7o4XFxVI8fxrAAxixAxAU/Fr2qA4dTHWuJRcrP/5oamoqTR0AAG7Dr2YPcxYwTJsWu+fs2tXUTz+N3XMCAIDy4VSsh+XmmpEzKTaLGFg0Aa/iVCyAoGDEzsP2/v20cWP0n+/220199NHoPxcAAKg4Ruw8bu5cs4CiQwdpwYLoPpezaKKkJHwMeAEjdgCCghE7j+vSxdQvv4zu6dFXXzX1jDNo6gAAcCsaOx+45RZTR4yI3nNceaWpsVyoAQAAKobGzgceftjUO+6IzuNPnWpq+/ZScnJ0ngPBEAqFdN9996l58+aqUaOGunXrpmXLlpX74zMzMxUXF6chQ4ZEMSUAeBeNnQ8kJIQvQzJ2bOQfv2dPUz/6KPKPjWAZMWKEHnvsMT399NNasGCBmjZtqjPOOEN5eXkH/dgFCxZo9OjROvbYY2OQFAC8icbOJ2bMMLV//8g+rjNa17q1lJIS2cdGsIRCIY0aNUrDhg3TX//6V6Wnp2vs2LHatWuXxo8ff8CP3bFjhy6//HKNGTNG9erVi1FiAPAeGjufSE4Oj9o9/3xkHjMUCo/WLV4cmcdEcP3yyy/Kzs5Wjx49Su9LSkpS165dNW/evAN+7KBBg3TOOeeoe/fu0Y4JAJ6WaDsAImfGDKlGDWngQOnqq80p2qpwRv+uvTZ8IWSgsrKzsyVJTZo0KXN/kyZNtHLlyv1+3BtvvKFFixZpQQWu55Ofn6/8/PzSt3NzcyuYFgC8iRE7H0lOlv75T3N81llVe6ytW8OXOHnuuao9FoLp9ddfV+3atUtvhYWFkqS4310vJxQK/eE+x+rVq3XjjTdq3LhxSq7Ayp3MzEylpqaW3lq2bFn5TwQAPIQLFPuQ8zty/nzppJOq9hjvvCOdc05kciFY8vLytGHDhtK38/PzlZ6erkWLFumEE04ovf+CCy5Q3bp1NXYfK3+mTp2qnj17KmGv4efi4mLFxcUpPj5e+fn5Zf5t7+f6/Yhdy5YtuUAxAN/jVKwPffONlJ4unXyylJNTduux8nDm1bVtS1OHyktJSVHKXituQqGQmjZtqlmzZpU2dgUFBZo9e7YeeeSRfT7G6aefrqVLl5a5b8CAATrqqKN0++2377Opk8zcvaSkpAh9JgDgHTR2PnT00dIzz0j/+IeZG1eRMdkJE8IrYX/6KTr5EEzO9eeGDx+uww8/XIcffriGDx+umjVr6rLLLit9v9NPP109e/bU4MGDlZKSovT09DKPU6tWLTVo0OAP9wMAaOx86+9/N03a3LlSRkb59pFdskRyfr9u2sTWYYi82267Tbt379Y//vEPbdu2TZ06ddLMmTPLjOytWLFCmzdvtpgSALyLOXY+5zRnRx8tLV26/2btf/+Tzj7bHM+cafaEBfwiNzdXqampzLED4HusivW53btNXbZMio83TdveduyQWrUKN3WzZ9PUAQDgVZyK9bnkZDPHrnt36YMPpDPP3P/7rlxpmjwAAOBNjNgFxPvvS1u2SHtdZaLU5Mmm+aOpAwDA2xixC5D69aVFi2ynAITBbwcAAAChSURBVAAA0cKIHQAAgE/Q2AEAAPgEjR0AAIBP0NgBAAD4BI0dAACAT9DYAQAA+ASNHQAAgE/Q2AEAAPhEXCgUCtkOAQDRFAqFlJeXp5SUFMXFxdmOAwBRQ2MHAADgE5yKBQAA8AkaOwAAAJ+gsQMAAPAJGjsAAACfoLEDAADwCRo7AAAAn6CxAwAA8AkaOwAAAJ+gsQMAAPAJGjsAAACf+H9RZZG6vyS+1wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "Graphics object consisting of 1 graphics primitive"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from ore_algebra.analytic.function import DFiniteFunction\n",
    "P.<x> = QQ[]\n",
    "A.<Dx> = 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": {
  "kernelspec": {
   "display_name": "SageMath 8.9.beta8",
   "language": "sage",
   "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,
   "scroll": true,
   "slideNumber": false,
   "transition": "none"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}