Demo: SageMath's Asymptotic Expansion Module

Demo - ACSV.ipynb

Demo - Asymptotic Expansions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "feb0e9ed",
   "metadata": {},
   "source": [
    "# Asymptotic Expansions in SageMath\n",
    "*Benjamin Hackl, Thomas Hagelmayer, Clemens Heuberger, Daniel Krenn, ...*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b8c5bc09",
   "metadata": {},
   "source": [
    "- Overview: https://doc.sagemath.org/html/en/reference/asymptotic/index.html\n",
    "- Documentation of the main structure, `AsymptoticRing`, including examples: https://doc.sagemath.org/html/en/reference/asymptotic/sage/rings/asymptotic/asymptotic_ring.html"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c15a5d71",
   "metadata": {},
   "source": [
    "An asymptotic expansion is an expression like\n",
    "$$ 5z^3 + 4z^2 + O(z) $$\n",
    "as $z\\to\\infty$, or\n",
    "$$ 3 x^42 y^2 + 7 x^3 y^3 + O(x^2) + O(y) $$\n",
    "as $x, y\\to\\infty$. It is a truncated series (after a finite number of terms), which approximates a function.\n",
    "\n",
    "This module allows to create and to carry out computations with such expressions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3a0e14cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "%display typeset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e46a3cd2",
   "metadata": {},
   "outputs": [],
   "source": [
    "assume(SR.an_element() > 0)  # required for some computations"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ce733f5",
   "metadata": {},
   "source": [
    "## Creating an \"Asymptotic Ring\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ff5bd31a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# first we create a suitable coefficient ring,\n",
    "# the symbolic ring without variables\n",
    "SCR = SR.subring(no_variables=True)\n",
    "SCR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e0a674e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "A = AsymptoticRing(\n",
    "    growth_group='QQ^n * n^QQ * log(n)^QQ',  # \"structure\" of the expansions\n",
    "    coefficient_ring=SCR,\n",
    "    default_prec=7,  # default precision\n",
    ")\n",
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b03c1b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "n = A.gen()  # variable n in the asymptotic ring\n",
    "n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9aa83006",
   "metadata": {},
   "source": [
    "Computations work intuitively:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f31e9b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "expansion = (\n",
    "    42 * (1/3)^n * n^(5/2) * log(n)^2 \n",
    "    + O((1/3)^n * n * log(n)^(1/2))\n",
    ")\n",
    "expansion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c59d6fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# with typesetting enabled, a string-based\n",
    "# representation can be printed:\n",
    "print(expansion)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6aead5cb",
   "metadata": {},
   "source": [
    "## Basic Expansion Arithmetic\n",
    "\n",
    "Addition, subtraction, multiplication and absorption (via $O$-terms) is easy!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1afee8b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "(n^2 + n + 1) + (2*n + O(n^(-1)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cde74baf",
   "metadata": {},
   "outputs": [],
   "source": [
    "n^2 + O(n) - O(n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23ca50d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "O(n^2 + n + 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb9760a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "(1 + 1/n) * (n^3 - n^2 + n + O(n^0))  # O(n^0) = O(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e358934",
   "metadata": {},
   "source": [
    "Furthermore, division, applying the exponential function and the logarithm are possible as well. If necessary, the expressions are expanded as series automatically---the precision of this expansion depends on the default_prec-parameter specified during creation of the asymptotic ring.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0ee364a",
   "metadata": {},
   "outputs": [],
   "source": [
    "(n - 1) / (n^3 + n^2 + O(n))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9be4b3a",
   "metadata": {},
   "outputs": [],
   "source": [
    "exp(1 + 1/n^2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3c2e2132",
   "metadata": {},
   "outputs": [],
   "source": [
    "log(n+1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a570d2da",
   "metadata": {},
   "source": [
    "Advanced computations are possible as well, e.g. computing the power of some asymptotic expansion to another asymptotic expansion:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4fec221b",
   "metadata": {},
   "outputs": [],
   "source": [
    "(1 + 1/n)^n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d41a91c8",
   "metadata": {},
   "source": [
    "### Experimental: Bounded Error Terms / B-Terms"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a0437645",
   "metadata": {},
   "source": [
    "An initial implementation of B-Terms (bounded error terms; \"O-Terms with explicit bound\") has been the goal of a Google Summer of Code project in 2021."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf4b1bb3",
   "metadata": {},
   "outputs": [],
   "source": [
    "A.<n> = AsymptoticRing(\"n^QQ\", SCR, default_prec=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca3e52a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "A.B(42*n^2, valid_from=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a23ea90e",
   "metadata": {},
   "source": [
    "A B-Term $B_{n\\geq 5}(42 n^2)$ is a term representing all growth functions that are bounded by $42 n^2$ given that $n\\geq 5$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c85f5d75",
   "metadata": {},
   "outputs": [],
   "source": [
    "A.B(42*n^2, valid_from=5) - A.B(5*n^2, valid_from=42)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc05e727",
   "metadata": {},
   "source": [
    "O-Terms *cannot* be absorbed by higher-order B-Terms!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "627f4c0a",
   "metadata": {},
   "outputs": [],
   "source": [
    "A.B(42*n^2, valid_from=5) + O(n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "01c46c66",
   "metadata": {},
   "outputs": [],
   "source": [
    "A.B(42*n^2, valid_from=5) + n^2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c27e7522",
   "metadata": {},
   "outputs": [],
   "source": [
    "A.B(42*n^2, valid_from=5) + A.B(5*n, valid_from=100)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b2a94f1",
   "metadata": {},
   "source": [
    "## Asymptotic Expansion Generators\n",
    "\n",
    "For the sake of convenience, asymptotic expansions of certain well-known objects (like binomial coefficients $\\binom{kn}{n}$) can be accessed easily via so-called generators. Generators are called via `asymptotic_expansion`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5b77d134",
   "metadata": {},
   "outputs": [],
   "source": [
    "A = AsymptoticRing(\n",
    "    \"QQ^n * n^QQ * log(n)^QQ\",\n",
    "    coefficient_ring=SR,\n",
    "    default_prec=5\n",
    ")\n",
    "n = A.gen()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4fe00e66",
   "metadata": {},
   "source": [
    "#### Harmonic Numbers\n",
    "... are defined as $H_n = \\sum_{k=1}^n \\frac{1}{k}$. We have implemented their expansion directly via a generator:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e3f2701a",
   "metadata": {},
   "outputs": [],
   "source": [
    "asymptotic_expansions.HarmonicNumber(n, precision=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fff050b3",
   "metadata": {},
   "source": [
    "If you did not define a suitable asymptotic ring before, the string name of the variable name to be used can be passed:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9072edf1",
   "metadata": {},
   "outputs": [],
   "source": [
    "asymptotic_expansions.HarmonicNumber('n', precision=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2732cf9",
   "metadata": {},
   "source": [
    "#### Catalan Numbers\n",
    "... are given by $C_n = \\frac{1}{n+1} \\binom{2n}{n}$, and can be obtained by dividing the asymptotic expansion of the central binomial coefficient (available via `asymptotic_expansions.Binomial_kn_over_n`) by $n+1$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79ac23f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "catalan_n = (\n",
    "    asymptotic_expansions.Binomial_kn_over_n(n, k=2, precision=5)\n",
    "    / (n + 1)\n",
    ")\n",
    "catalan_n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a640fb2",
   "metadata": {},
   "source": [
    "## Singularity Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7bebda52",
   "metadata": {},
   "source": [
    "For explicity given generating functions, our module can extract the coefficient growth by means of singularity analysis directly: this is implemented as a method of the underlying asymptotic ring.\n",
    "\n",
    "Consider the well-known\n",
    "$$ \\mathrm{Catalan}(z) = \\frac{1 - \\sqrt{1 - 4z}}{2z}. $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fdced1b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "def catalan(z):\n",
    "    return (1 - sqrt(1 - 4*z)) / (2*z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f9ee4a54",
   "metadata": {},
   "outputs": [],
   "source": [
    "C_n = A.coefficients_of_generating_function(\n",
    "    function=catalan,\n",
    "    singularities=(1/4,),\n",
    "    precision=5\n",
    ")\n",
    "C_n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05c3bfdd",
   "metadata": {},
   "outputs": [],
   "source": [
    "C_n - catalan_n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae00bf53",
   "metadata": {},
   "source": [
    "Harmonic numbers have the generating function\n",
    "$$ \\sum_{n\\geq 0} H_n z^n = -\\frac{\\log(1 - z)}{1 - z}. $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "84a0b126",
   "metadata": {},
   "outputs": [],
   "source": [
    "def harmonic_gf(z):\n",
    "    return - log(1 - z) / (1 - z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "abcd1600",
   "metadata": {},
   "outputs": [],
   "source": [
    "H_n = A.coefficients_of_generating_function(\n",
    "    function=harmonic_gf,\n",
    "    singularities=(1,),\n",
    "    precision=10\n",
    ")\n",
    "H_n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb8015d2",
   "metadata": {},
   "outputs": [],
   "source": [
    "asymptotic_expansions.HarmonicNumber('n', precision=10) - H_n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08ca9c54",
   "metadata": {},
   "source": [
    "#### Inverse Functions\n",
    "There is also a generator that extracts the coefficient growth of a function $y(z)$ satisfying the implicit equation $y(z) = z \\Phi(y(z))$ (\"singular inversion\" -- Flajolet–Sedgewick, Analytic Combinatorics, Chapter VI.7); the *usual conditions* with respect to $\\Phi$ need to hold."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8c1301d",
   "metadata": {},
   "source": [
    "*Example:* the generating function $y(z)$ of unary-binary-ternary trees satisfies\n",
    "$$ y = z (1 + 3y + 3y^2 + y^3). $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f10de465",
   "metadata": {},
   "outputs": [],
   "source": [
    "def unary_binary_ternary_trees(y):\n",
    "    return 1 + 3*y + 3*y^2 + y^3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31657b01",
   "metadata": {},
   "outputs": [],
   "source": [
    "tree_asy = asymptotic_expansions.InverseFunctionAnalysis(\n",
    "    var=n,\n",
    "    phi=unary_binary_ternary_trees,\n",
    "    precision=7,\n",
    ")\n",
    "tree_asy"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d2358438",
   "metadata": {},
   "source": [
    "## Bootstrapping\n",
    "\n",
    "**Bootstrapping** in order to find the dominant singularity is easily possible as well. As an example, consider the longest run of words over a two letter alphabet (see *Flajolet--Sedgewick*: *Analytic Combinatorics*, Example V.4). The generating function counting runs where one of the two letters has less than $n$ consecutive repetitions is\n",
    "$$ \\frac{1 - z^n}{1 - 2z + z^{n+1}}. $$\n",
    "Then, the dominant singularity satisfies the fixed-point equation $z = \\frac{1 + z^{n+1}}{2} =: f(z)$. Applying $f$ repeatedly to the approximation $z = \\frac{1}{2} + O\\Big(\\Big(\\frac{3}{5}\\Big)^n\\Big)$ increases the precision of the expansion of the dominant singularity:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41a78e5f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# change default_prec to get less/more precise expansion\n",
    "A.<n> = AsymptoticRing(\"QQ^n * n^QQ\", QQ, default_prec=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "543e06ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(z):\n",
    "    return (1 + z^(n+1))/2\n",
    "\n",
    "z_0 = 1/2 + O((3/5)^n)\n",
    "z_0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0de856f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "z_1 = f(z_0)\n",
    "z_1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c23e176b",
   "metadata": {},
   "outputs": [],
   "source": [
    "z_2 = f(z_1)\n",
    "z_2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ee25795",
   "metadata": {},
   "outputs": [],
   "source": [
    "z_3 = f(z_2)\n",
    "z_3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "501ffd17",
   "metadata": {},
   "outputs": [],
   "source": [
    "z_4 = f(z_3)\n",
    "z_4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3465155f",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 9.7",
   "language": "sage",
   "name": "sagemath"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}

Dockerfile

FROM sagemath/sagemath:9.7

ARG NB_UID=1000
ARG NB_USER=sage

USER root
RUN apt update && apt install -y python3 python3-pip

USER ${NB_UID}
ENV PATH="${PATH}:${HOME}/.local/bin"
RUN pip3 install notebook
RUN sage -pip install sage_acsv
RUN ln -s $(sage -sh -c 'ls -d $SAGE_VENV/share/jupyter/kernels/sagemath') $HOME/.local/share/jupyter/kernels/sagemath-dev

# partially superfluous -- create separate directory to hold notebooks
WORKDIR ${HOME}/notebooks
COPY --chown=${NB_UID}:${NB_UID} . .
USER root
RUN chown -R ${NB_UID}:${NB_UID} .
USER ${NB_UID}

ENTRYPOINT []