Automatic Asymptotics with Sage

20 Jan 2022 - Tags: sage

Recurrence relations show up all the time in combinatorics and computer science, and even simple objects give recurrences that are difficult or impossible to solve. Thankfully, we can often find a generating function for our objects, and through the power of complex analysis and algebraic geometry, we can use the singularities of the generating function in order to get good asymptotic estimates. Software like Maple can automatically compute asymptotic expansions for a lot of generating functions using the asympt function… Can sage?

The answer is “yes”, which is why I’m writing this post, but the longer answer is “with some work”. Which is the real reason I’m writing this post. I had to slog through a lot of kind of crummy documentation to get this working, and I want to make sure that the next people looking into this have an easier time.

There are two modules (that I can find) which provide features for computing asymptotics:

• multivariate generating functions
• asymptotic rings

These both have pros and cons:

The former is slightly easier to work with, and allows generating functions with multiple variables. The downside is that it only works for functions with polynomial denominators¹.

The latter has nicer documentation, works for more general functions, and gives more detailed asymptotic expansions. The downside is that it only works for functions of a single variable, and it’s kind of annoying to use because it doesn’t accept symbolic inputs natively. You also have to manually provide a list of singularities².

---

Let’s start with the multivariate case, since it’s more complicated. I’ll give a code dump, then explain what’s going on.

The tl;dr is that it works for functions of the form $h / p$ where $h$ is a holomorphic function and $p$ is a polynomial. The number of variables is irrelevant.

This is pretty obviously cribbed from the examples here, but the basic idea is this:

We need the denominator of our generating function to be a factored polynomial, so we build a polynomial ring with the variables in the denominator and factor over that ring.

Then, we put this into FFPD to get an object that the module knows how to deal with.

We choose a multi-index $\alpha$ which controls the “direction” in which we take our asymptotics³. For instance, if our generating function is $f(x,y) = \sum_{m,n} F_{m,n} x^m y^n$, then

• $\alpha = (1,1)$ will give us the asymptotics of $F_{n \alpha} = F_{n,n}$ as $n \to \infty$
• $\alpha = (2,1)$ will give us the asymptotics of $F_{n \alpha} = F_{2n,n}$ as $n \to \infty$
• etc.

If we view the $F_{m,n}$ as being located at the lattice points $(m,n)$ in the plane, then $\alpha$ picks out the direction in which we move.

Next we look at the critical points of the denominator. That is, the points where its gradient is either undefined or vanishing. The former is the singular case, and the latter is the smooth case.

In the one variable case, we know the asymptotics are controlled by the singularity closest to the origin. See chapter $5$ of Generatingfunctionology, for instance.

From skimming these articles, we now have a whole algebraic variety of singularities, and for reasons the asymptotics are now governed by the critical points on this variety. In particular, we’re on the hunt for the critical point which is closest to the origin (which we’ll call the Dominant Singularity).

So now, what does the code do? Well we look for the location of the minimal smooth point and the minimal singular point. Then take the smaller one and run the asymptotic function provided by the module.

It seems to work quite well, and is surprisingly fast too. For a simple example, we can try

which outputs⁴

\[\frac{1}{5} \, \sqrt{5} \left(\frac{2}{\sqrt{5} - 1}\right)^{r}\]

Of course, \(\frac{z}{1-z-z^2}\) is the generating function for the fibonacci numbers, and we’ve successfully recovered the asymptotic formula (though it requires some massaging to make it look like its usual presentation).

This really shines when it comes to multivariate series, though. For instance, we know that

\[\binom{n}{k} = [x^n y^k] (1 - x(1+y))^{-1}\]

So if we’re interested in the asymptotics of $\binom{3n}{n}$ we can compute

which quickly gives

\[\frac{1}{41472} \, \left(\frac{27}{4}\right)^{r} {\left(\frac{10368 \, \sqrt{6} \sqrt{2}}{\sqrt{\pi} \sqrt{r}} - \frac{1008 \, \sqrt{6} \sqrt{2}}{\sqrt{\pi} r^{\frac{3}{2}}} + \frac{49 \, \sqrt{6} \sqrt{2}}{\sqrt{\pi} r^{\frac{5}{2}}}\right)}\]

Now $\binom{15}{5} = 3003$ and evaluating the above at $r=5$ gives $3002.931$.

---

But what if we’re interested in generating functions whose numerator isn’t holomorphic, or whose denominator isn’t a polynomial? It’s clear that we should be interested in such things, since the famous catalan numbers already have the generating function

\[\frac{1 - \sqrt{1 - 4z}}{2z}\]

For cases like this (which must be of a single variable), we work with the second asymptotics package here.

This is nice because it returns an asymptotic expansion, which is quite intuitive to work with. The major downside, though, might be surprising:

It only works with (callable) python functions.

This is annoying⁵, but it really isn’t that much of a hassle to rewrite your function using def. We also need to provide the dominant singularities by hand, which is also a bit of work⁶.

If and when I get around to automatically finding singularities, I’ll update this code, but for now there’s much less boilerplate here than in the previous case:

The documentation for this function is also quite good, so I’ll stick to one example. For the catalan numbers, we’ll have:

This looks like it has singularities at $0$ (which makes the denominator vanish) and $1/4$ (which makes the $\sqrt{\cdot}$ vanish), but actually it only has the singularity at $1/4$. The singularity at $0$ is removable. So we call

which outputs

\[\frac{1}{\sqrt{\pi}} 4^{n} n^{-\frac{3}{2}} - \frac{9}{8 \, \sqrt{\pi}} 4^{n} n^{-\frac{5}{2}} + \frac{145}{128 \, \sqrt{\pi}} 4^{n} n^{-\frac{7}{2}} - \frac{1155}{1024 \, \sqrt{\pi}} 4^{n} n^{-\frac{9}{2}} + \frac{36939}{32768 \, \sqrt{\pi}} 4^{n} n^{-\frac{11}{2}} + O\!\left(4^{n} n^{-6}\right)\]

A quick massage turns this into the (maybe not so) familiar

\[\frac{4^n}{\sqrt{\pi}} \left ( n^{-3/2} - \frac{9}{8} n^{-5/2} + \frac{145}{128} n^{-7/2} - \frac{1155}{1024} n^{-9/2} + \frac{36939}{32768} n^{-11/2} \pm O \left ( n^{-13/2} \right ) \right )\]

which you can find as formula (20) here.

---

Lastly, we can combine these different asymptotic approximations into one function. In my init.sage file, which you can find here, I’ve defined the function

This is definitely too brittle to give to other people, but it works fine for my purposes. One day I’d like to make it a bit more robust (this is not how try/except statements are supposed to be used), and ideally make it a bit more automatic (in particular following the discussion of singlevar_asy above). Another option might be to try and get it working with the singularlity_analysis function, which would require a certain amount of parsing… Anyways, as it stands it’s a plenty servicable substitute for maple’s asympt function, and I’m glad I took the time to code it up.

If anyone has any other ideas for making this more usable, or encounters any problems, definitely let me know!

If not, take care and stay warm!

I’ll see you all in the next one ^_^.

---