{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Differential operators"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "from ore_algebra import OreAlgebra"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "Pols.<z> = PolynomialRing(QQ)\n",
    "DiffOps.<Dz> = OreAlgebra(Pols)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DiffOps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Dz(sin(z))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "Dz*z"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Toy Examples 1: Numerical Evaluation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(Dz - 1).numerical_solution(ini=[1], path=[0, 1]) # variants..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "⤷ Note that the output is an interval!<br/><br/>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "dop = Dz*z*Dz\n",
    "dop(log(z))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dop.numerical_solution(ini=[0, 1], path=[1, 2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#dop.numerical_solution(ini=[0, 1], path=[1, -1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "dop.numerical_solution(ini=[0, 1], path=[1, i, -1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dop.numerical_solution(ini=[0, 1], path=[1, -i, -1])"
   ]
  },
  {
   "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": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from ore_algebra.examples import polya\n",
    "dim = 15\n",
    "polya.dop[dim]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "polya.dop[dim].leading_coefficient()(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "polya.dop[dim].local_basis_monomials(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "ini = [0]*(dim - 1) + [1]\n",
    "polya.dop[dim].numerical_solution(ini, [0,1/(2*dim)], 1e-100)"
   ]
  },
  {
   "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": [
    "### Local monodromy matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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\n",
    "dop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dop.leading_coefficient().factor()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dop.numerical_transition_matrix([1/2, i/2, -1/2, -i/2, 1/2], 1e-5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "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": [
    "### 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": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "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": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "dop = cbt.dop[5]; dop"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Singular points:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "sing = dop.leading_coefficient().roots(QQbar, multiplicities=False)\n",
    "sing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "1/(4*(cos(pi/8)^2)).n()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Local behaviors at the dominant singularity:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "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": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "dop.local_basis_monomials(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "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": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "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": [
    "# Bonus Examples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Local monodromy matrices (ctd)\n",
    "\n",
    "Local monodromy = formal monodromy + singular connection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dop.numerical_transition_matrix([1/2, I/2, -1/2, -I/2, 1/2], 1e-5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dop.local_basis_monomials(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "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": [
    "### 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": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from ore_algebra.examples import iint\n",
    "i = 69; iint.word[i]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "dop = iint.diffop(iint.word[i])\n",
    "dop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "iint.ini[i]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "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": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "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": [
    "### Asymptotics of Apéry Numbers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mat = dop.numerical_transition_matrix([0, sing[1]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b0 = 0; b1 = 6\n",
    "d1 = b0*mat[1,2] + b1*mat[1,3]\n",
    "d1/c1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "zeta(3.)"
   ]
  },
  {
   "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": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "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": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "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": null,
   "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": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "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": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "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": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "den_ratios"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "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": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "SageMath 8.9.beta1",
   "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
}