Rigorous Asymptotic Expansions & ACSV in SageMath
Benjamin Hackl • June 29, 2023

Rigorous Asymptotic Expansions and Analytic Combinatorics in Several Variables in SageMath

A "hands-on" software session • AofA2023 @ Taipei
This work is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License.

Interactive Demo

behackl.dev/asy
Notebook sources: behackl.dev/asy/src

Asymptotic Expansions in SageMath

Joint work with Thomas Hagelmayer, Clemens Heuberger, Daniel Krenn, ...
  • Want: framework for computations with generalized power series with rigorous control over error terms 
sage: A.<n> = AsymptoticRing('n^QQ', SR, default_prec=4)
sage: (1 + 1/n)^n
e - 1/2*e*n^(-1) + 11/24*e*n^(-2) - 7/16*e*n^(-3) + O(n^(-4))
  • Initial implementation as a Google Summer of Code project in 2015
  • Gradual improvement and refactors over the years
  • $B$-terms (explicitly bounded $O$-terms) in GSoC'21

Basic Design Ideas

  • GrowthGroup $\ni$ GrowthElement, holds only the growth
  • TermMonoid $\ni$ AsymptoticTerm, summands of the expansion, $42\log(n)^3$, $O(n^{42})$
  • AsymptoticRing $\ni$ AsymptoticExpansion, sum of (partially ordered) terms
Development principles of SageMath:
  • every (new) class / function needs documentation ...
  • ... and tests!

An $\varepsilon$ more about $B$-terms

  • Concept appears / is used occassionally (e.g.: in de Bruijn's [Asymptotic Methods in Analysis] as a didactic vehicle, Study of solutions of Diophantine equations, ...)
  • Formally: $f(n) = B_{n\geq n_0}(c\cdot g(n)) \iff \forall n\geq n_0: |f(n)| \leq c |g(n)|$
  • $B$-Term arithmetic, pop quiz: simplify!
    1. $11n + B_{n\geq 1}(3n)$
    2. $n^2 + B_{n\geq 10}(7 n) + B_{n\geq 42}(19/n) + O(1/n)$
    3. $B_{n\geq 7}(5n^3) + B_{n\geq 10}(3n)$

Something that works out of the box ...

The average size of $r$-fold cut plane trees ORTs: \[ \mathbb{E} X_{n,r} = \frac{1}{C_{n-1}} [z^n] \frac{u^{r+1}}{(1+u) (1 - u^{r+1})}, \qquad z = u/(1+u)^2 \] 

>>> catalan_shifted_asy = asymptotic_expansions.Binomial_kn_over_n(n, k=2, precision=3)
>>> catalan_shifted_asy = catalan_shifted_asy.subs(n=n-1) / n
>>> catalan_shifted_asy
1/4/sqrt(pi)*4^n*n^(-3/2) + 3/32/sqrt(pi)*4^n*n^(-5/2) + 25/512/sqrt(pi)*4^n*n^(-7/2) + O(4^n*n^(-9/2))
>>> z_asy = 1/4 * (1 - Z^-1)
>>> u_asy = upsilon(z_asy) + O(Z^-2); u_asy
1 - 2*Z^(-1/2) + 2*Z^(-1) - 2*Z^(-3/2) + 2*Z^(-2) - 2*Z^(-5/2) + O(Z^(-2))
>>> singular_expansion = fast_callable(u^(r+1) / (1 + u) / (1 - u^(r+1)), vars=(u, r))(u_asy, r)
>>> singular_expansion
(1/4/(r + 1))*Z^(1/2) + 1/4*((r + 1)^2 + r + 1)/(r + 1)^2 - 1/2 + (1/2*r - 1/2*((r + 1)^2 + r + 1)/(r + 1) + 1/12*(3*((r + 1)^2 + r + 1)^2/(r + 1)^2 - (2*(r + 1)^3 + 3*(r + 1)^2 + 4*r + 4)/(r + 1))/(r + 1) + 1/2)*Z^(-1/2) + O(Z^(-1))
>>> expectation = singular_expansion._singularity_analysis_(n, zeta=1/4, precision=3) / catalan_shifted_asy
>>> expectation.map_coefficients(lambda t: t.simplify_rational())
(1/(r + 1))*n - 1/6*(r^2 - r)/(r + 1) + O(n^(-1/2))
Theorem.[H.-Heuberger-Kropf-Prodinger, 2018]

After $r$ "cuts", the average number of remaining nodes is \[ \mathbb{E} X_{n,r} = \frac{n}{r+1} - \frac{r(r-1)}{6(r+1)} + O(n^{-1/2}). \]

... and something that does not (yet!)

$$ F(n) = \sum_{k=1}^n k \sigma(k) (k^2-3n+2)(2k^2-n) \binom{2n}{n-k} < 0 $$ for all $n \geq 5$?

Maybe at next year's AofA? 😄

Rigorous ACSV in SageMath: sage_acsv

Joint work with Andrew Luo, Stephen Melczer, Jesse Selover, Elaine Wong
  • Input: multivariate, rational, combinatorial generating function \[ F(\mathbf{z}) = \sum_{i_1, i_2, ..., i_d \geq 0} f_{i_1, ..., i_d} z_1^{i_1} \cdots z_d^{i_d}, \] and direction vector $\mathbf{r} = (r_1, ..., r_d) \in \mathbb{N}^d$
  • Output: asymptotic expansion for diagonal $f_{n \mathbf{r}}$ for $n\to\infty$

Installation & Demo

How does it work?

  • $F(\mathbf{z}) = G(\mathbf{z})/H(\mathbf{z})$, with $G$ and $H$ (multivariate) polynomials
  • Multivariate CIF: $f_{n\mathbf{r}} = \frac{1}{(2\pi i)^d} \int_{\mathcal{C}} \frac{G(\mathbf{z})}{H(\mathbf{z})} \frac{d\mathbf{z}}{\mathbf{z}^{n \mathbf{r} + 1}}$
  • Singularity set is algebraic variety, has infinitely many "minimal" points...
  • Under "mild" conditions on $H$ ($H(0)\neq 0$; $H$ and all partial derivatives do not vanish simultaneously; finitely many smooth critical points, at least one of which is minimal; phase Hessian non-singular at these points) the special ("critical") minimal points determine asymptotics cf. [Pemantle-Wilson], [Baryshnikov-Pemantle] ...

Algorithmic Approach

  • Find critical points $\mathbf{w} \in \mathbb{C}_{*}^d$ (via $H(\mathbf{w}) = 0$ and $\mathbf{0}\neq \operatorname{diag}(\mathbf{w}) \nabla H(\mathbf{w}) \| \mathbf{r}$)
  • Check which critical points are minimal

Use all-in-one strategyvia [Melczer-Salvy]:

  • Consider system $H(t \mathbf{w}) = 0$, $\operatorname{diag}(w) \nabla H(\mathbf{w}) = \lambda \mathbf{r}$, encode solutions with Kronecker representationlinear form $u - \kappa \mathbf{z} = 0$, polynomials $P$, $Q_1$, ..., $Q_d$ s.t. $Q_j(u)/P'(u) = z_j$ determine components of solution vector for $u$ iterating through roots of $P$ -- bounded coefficients / degrees!
  • Check positive ($\sim$ Pringsheim!) solutions for $t = 1$
  • Use solutions for $t\in (0,1)$ to eliminate non-minimal $\mathbf{w}$

Thank you!