Let $(f_n) = f_0,f_1,\dots$ be a complex-valued sequence with generating function
$$ F(z) = \sum_{n \geq 0} f_nz^n $$
$f_n$ is P-recursive if it satisfies a non-trivial linear recurrence with poly coeffs
$$ 0 = c_r(n)f_{n+r} + \cdots + c_1(n)f_{n+1} + c_0(n)f_n$$
$F(z)$ is D-finite if it satisfies a non-trivial linear ODE with poly coeffs
$$ 0 = p_s(z)F^{(s)}(z) + \cdots + p_1(z)F'(z) + p_0(z)F(z)$$
Question 1: How do we find asymptotics of a P-recursive sequence?
Question 2: When does a P-recursive sequence have all positive terms?
Motivation: Canham model for biomembranes
Canham Problem: For fixed genus $g$ and fixed isoperimetric ratio $\iota_0$ find the surface $S$ minimizing
$$E(S) = \int_S (\kappa_1 + \kappa_2)^2 dA$$
Marques and Neves (2014): Solution existence for genus one
Yu and Chen (2020): The genus one model has a unique solution (for a range of isoperimetric ratios) if the sequence $(d_n)$ with initial terms
$$(d_0,\dots,d_6) = \left(72, 1932, 31248, \frac{790101}{2}, \frac{17208645}{4}, \frac{338898609}{8},\frac{1551478257}{4}\right)$$
satisfying
$$ r_7(n)d_{n+7} + \cdots + r_0(n)d_{n} = 0$$
has all positive terms where the $r_j(n)$ are explicit polynomials.
M. and Mezzarobba (2021): This sequence has all positive terms.
Dong, M. and Mezzarobba (2022): Extend to general P-recursive sequences (and Sage implementation!).
Shift Algebra¶
We can model a P-recursion using the Shift algebra consisting of polynomials in $S_n$ and $n$ such that $S_n n = (n+1)S_n$.
Such polynomials act on sequences via
$$n \cdot (f_n) = (nf_n) \qquad\qquad S_n \cdot (f_n) = (f_{n+1})$$
The Shift algebra is implemented in the Sage ore_algebra package (Kauers, Jaroschek, Johansson, Mezzarobba, Verron, ...)
# Define the ring of shift operators
from ore_algebra import *
Ind.<n> = PolynomialRing(QQ); Shift.<Sn> = OreAlgebra(Ind)
# Algebraic Simplification
Sn * n^10
(n^10 + 10*n^9 + 45*n^8 + 120*n^7 + 210*n^6 + 252*n^5 + 210*n^4 + 120*n^3 + 45*n^2 + 10*n + 1)*Sn
The package can generate sequence terms efficiently
# Central binomial coefficients
rec = (n + 1)*Sn - 4*n - 2
rec.to_list([1],10)
[1, 2, 6, 20, 70, 252, 924, 3432, 12870, 48620]
[binomial(2*k,k) for k in range(10)]
[1, 2, 6, 20, 70, 252, 924, 3432, 12870, 48620]
It can also guess P-recursions
LST = [binomial(2*k,k)^2 for k in range(10)]
LST
[1, 4, 36, 400, 4900, 63504, 853776, 11778624, 165636900, 2363904400]
guess(LST,Shift)
(-n^2 - 2*n - 1)*Sn + 16*n^2 + 16*n + 4
Lattice Path Example¶
Let $s_n$ be the number of lattice paths on the steps $\{\text{North, South, East, West}\} = \{(\pm1,0),(0,\pm1)\}$ that start at the origin and stay in $\mathbb{N}^2$.
We can numerically compute terms in the counting sequence, then guess a P-recurrence it satisfies.
@CachedFunction
def NSEW(i,j,n):
if i<0 or j<0: return 0
elif n==0 and i==0 and j==0: return 1
elif n==0: return 0
else: return NSEW(i-1,j,n-1) + NSEW(i,j-1,n-1) + NSEW(i+1,j,n-1) + NSEW(i,j+1,n-1)
def s(n): return add([NSEW(i,j,n) for i in range(n+1) for j in range(n+1)])
print([s(n) for n in range(10)])
[1, 2, 6, 18, 60, 200, 700, 2450, 8820, 31752]
rec = guess([s(n) for n in range(50)], Shift)
rec
(-n^2 - 7*n - 12)*Sn^2 + (8*n + 20)*Sn + 16*n^2 + 48*n + 32
The solutions of a P-recurrence form a $\mathbb{C}$-vector space, and we can compute series expansions of an asymptotic basis.
rec = (-n^2 - 7*n - 12)*Sn^2 + (8*n + 20)*Sn + 16*n^2 + 48*n + 32
show(rec.generalized_series_solutions(n=2))
Thus, there exist $\color{red}\lambda_{\color{red}1},\color{red}\lambda_{\color{red}2} \in \mathbb{C}$ such that
$$ s_n = \color{red}\lambda_{\color{red}1} \frac{4^n}{n}\left(1 - \frac{3}{2n} + \cdots\right) + \color{red}\lambda_{\color{red}2} \frac{(-4)^n}{n^3}\left(1 - \frac{9}{2n} + \cdots\right) $$
How do we find $\lambda_1,\lambda_2$?
How do we prove $\lambda_1 \neq 0$?
Move to differential equations and use transfer theorems (Flajolet-Odlyzko, etc.).
Differential Weyl Algebra¶
We can model a D-Finite equation using polynomials in $D_z$ and $z$ such that $D_z z = zD_z + 1$.
# Define the ring of differential operators
Pols.<z> = PolynomialRing(QQ); Diff.<Dz> = OreAlgebra(Pols)
Dz * z
z*Dz + 1
# Can guess differential equations (or convert rec <-> ODE)
diff = guess([s(n) for n in range(50)],Diff)
diff
(16*z^4 - z^2)*Dz^3 + (128*z^3 + 8*z^2 - 6*z)*Dz^2 + (224*z^2 + 28*z - 6)*Dz + 64*z + 12
# We can compute a 'local' basis of solutions
diff.local_basis_expansions(0,5)
[z^(-2) - 4*z^(-1)*log(z) - 8*log(z) - 16*z*log(z) - 8*z - 48*z^2*log(z) - 28*z^2, z^(-1), 1 + 2*z + 6*z^2]
# We can represent our GF in terms of a basis at the origin
ini = [0,0,1]
Techniques from numeric analysis, the study of ODEs, and computer algebra give practical algorithms for numeric analytic continuation (Chudnovskys, van der Hoeven, Mezzarobba, ...).
We obtain singular generating function expansions with rigorously approximated coefficients.
# Potential singularities for the GF are roots of the leading ODE coefficient
diff.leading_coefficient().factor()
(16) * (z - 1/4) * (z + 1/4) * z^2
# The generating function has a singular expansion at z=1/4 which begins
ini = [0,0,1]
bas = diff.numerical_transition_matrix([0,1/4],1e-10) * vector(ini)
loc = diff.local_basis_expansions(1/4,2)
add([a*b for [a,b] in zip(bas,loc)])
([-1.2732395447+/-3.84e-11]+[+/-5.44e-15]*I)*log(z - 1/4) + ([7.639437268+/-4.45e-10]+[+/-3.27e-14]*I)*(z - 1/4)*log(z - 1/4) + ([-2.3906971441+/-3.21e-11]+[4.0000000000+/-2.92e-11]*I)*1 + ([12.890661954+/-2.40e-10]+[-24.000000000+/-1.17e-10]*I)*(z - 1/4)
# Using the coefficient extraction [z^n]log(1-4*z) = -4^n/n dominant asymptotic behaviour is
var('N')
show(-bas[0].real() * 4^N/N)
What if we want explicit error bounds?
Flajolet-Odlyzko transfer theorems in principle give explicit upper bounds -- but they can be very painful to compute.
We combine D-finite computer algebra tools with transfer theorems to make everything automatic.
# Reimport needed functions
from ore_algebra import *
from ore_algebra import (OreAlgebra, DFiniteFunctionRing, UnivariateDFiniteFunction)
from ore_algebra.analytic.singularity_analysis import (bound_coefficients, check_seq_bound, eval_bound)
Pols.<z> = PolynomialRing(QQ); Diff.<Dz> = OreAlgebra(Pols);
Lattice Path Example¶
# Enter the differential operator and initial conditions
dop = (z^2*(4*z - 1)*(4*z + 1)*Dz^3 + 2*z*(4*z+1)*(16*z-3)*Dz^2 + 2*(112*z^2 + 14*z - 3)*Dz + 4*(16*z + 3))
bound_coefficients(dop, [1, 2, 6], order=2)
1.000000000000000*4^n*(([1.27323954473516 +/- 3.09e-15] + [+/- 5.53e-38]*I)*n^(-1) + ([-1.90985931710274 +/- 4.63e-15] + [+/- 8.29e-38]*I)*n^(-2) + B([1503.374850517149 +/- 3.85e-13]*n^(-3)*log(n), n >= 7))
Thus, the number $s_n$ of quadrant walks satisfies
$$ \begin{aligned} s_n \in \frac{4^{n}}{n} \cdot \biggl( & \left[1.27 \pm 4.39 \cdot 10^{-5}\right] + \left[-1.91 \pm 5.60 \cdot 10^{-5}\right] \frac{1}{n} + \left[1503.38\pm 3.05 \cdot 10^{-3}\right] \frac{\log^2 n}{n^3} \biggr) \end{aligned} $$
when $n \geq 7$.
Genus One Biomembrane Example¶
seqini = [72, 1932, 31248, 790101/2, 17208645/4, 338898609/8, 1551478257/4]
deq = (8388593*z^2*(3*z^4 - 164*z^3 + 370*z^2 - 164*z + 3)*(z + 1)^2*(z^2 - 6*z + 1)^2*(z - 1)^3*Dz^3 + 8388593*z*(z + 1)*(z^2 - 6*z + 1)*(66*z^8 - 3943*z^7 + 18981*z^6 - 16759*z^5 - 30383*z^4 + 47123*z^3 - 17577*z^2 + 971*z - 15)*(z - 1)^2*Dz^2 + 16777186*(z - 1)*(210*z^12 - 13761*z^11 + 101088*z^10 - 178437*z^9 - 248334*z^8 + 930590*z^7 - 446064*z^6 - 694834*z^5 + 794998*z^4 - 267421*z^3 + 24144*z^2 - 649*z + 6)*Dz + 6341776308*z^12 - 427012938072*z^11 + 2435594423178*z^10 - 2400915979716*z^9 - 10724094731502*z^8 + 26272536406048*z^7 - 8496738740956*z^6 - 30570113263064*z^5 + 39394376229112*z^4 - 19173572139496*z^3 + 3825886272626*z^2 - 170758199108*z + 2701126946)
deq
(25165779*z^15 - 1702884379*z^14 + 22196217078*z^13 - 108363844374*z^12 + 197744302789*z^11 - 20375892397*z^10 - 350408306796*z^9 + 350408306796*z^8 + 20375892397*z^7 - 197744302789*z^6 + 108363844374*z^5 - 22196217078*z^4 + 1702884379*z^3 - 25165779*z^2)*Dz^3 + (553647138*z^14 - 36951752165*z^13 + 394079321954*z^12 - 1450287066584*z^11 + 1482230828728*z^10 + 2523313940179*z^9 - 6434889690300*z^8 + 3026134593192*z^7 + 3073463034898*z^6 - 3856865346575*z^5 + 1475536731514*z^4 - 205218539152*z^3 + 9026126068*z^2 - 125828895*z)*Dz^2 + (3523209060*z^13 - 234394065606*z^12 + 1926843034914*z^11 - 4689642916650*z^10 - 1172674969842*z^9 + 19779027227864*z^8 - 23096380215644*z^7 - 4173660561220*z^6 + 24995188572752*z^5 - 17824401172934*z^4 + 4891640236090*z^3 - 415956772498*z^2 + 10989056830*z - 100663116)*Dz + 6341776308*z^12 - 427012938072*z^11 + 2435594423178*z^10 - 2400915979716*z^9 - 10724094731502*z^8 + 26272536406048*z^7 - 8496738740956*z^6 - 30570113263064*z^5 + 39394376229112*z^4 - 19173572139496*z^3 + 3825886272626*z^2 - 170758199108*z + 2701126946
bound_coefficients(deq, seqini, order=5, prec=20, n0=50)
1.00000*5.828427124746190?^n*(([8.072 +/- 1.77e-4] + [+/- 4.23e-16]*I)*n^3*log(n) + ([1.371 +/- 5.79e-4] + [+/- 1.08e-4]*I)*n^3 + ([50.509 +/- 8.16e-4] + [+/- 2.65e-15]*I)*n^2*log(n) + ([29.70 +/- 2.92e-3] + [+/- 6.71e-4]*I)*n^2 + ([106.36 +/- 2.08e-3] + [+/- 5.57e-15]*I)*n*log(n) + ([118.76 +/- 5.85e-3] + [+/- 1.42e-3]*I)*n + ([72.85 +/- 4.56e-3] + [+/- 3.82e-15]*I)*log(n) + [154.7 +/- 0.0133] + [+/- 9.56e-4]*I + ([-0.28378 +/- 3.33e-6] + [+/- 1.49e-17]*I)*n^(-1)*log(n) + ([35.5 +/- 0.0224] + [+/- 3.41e-6]*I)*n^(-1) + B([3056.71 +/- 2.97e-3]*n^(-2)*log(n)^2, n >= 50))
Corollary: The sequence has all positive terms. (Check finite number by hand then use explicit bound)
Difficulty for Positivity: Cancellation¶
If there are multiple dominant singularities, it can be hard to decide if cancellation occurs (open problem even for univariate functions!). This is not a worry for $\mathbb{N}$-rational functions (basically all the ones we care about in actual enumeration problems).
Difficulty for Asymptotics: Certifying Singularities¶
# Consider the sequence satisfying the P-recurrence
ini = [1,1/4]
rec = (n+3)^2*Sn^2 - (n+2)*(3*n+11)/2*Sn + (n+4)*(n+1)/2
rec
(n^2 + 6*n + 9)*Sn^2 + (-3/2*n^2 - 17/2*n - 11)*Sn + 1/2*n^2 + 5/2*n + 2
# The GF satisfies the D-finite equation
diff = rec.to_D(Diff)
diff
(1/2*z^4 - 3/2*z^3 + z^2)*Dz^4 + (7*z^3 - 16*z^2 + 7*z)*Dz^3 + (26*z^2 - 41*z + 9)*Dz^2 + (26*z - 22)*Dz + 4
# There is an apparent singularity at z=1
diff.leading_coefficient().factor()
(1/2) * (z - 2) * (z - 1) * z^2
It is not an actual singularity but the algorithm cannot detect this (the famous "connection problem")
bound_coefficients(diff, ini, prec=1000)
1.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*([+/- 1.01e-305] + [+/- 2.15e-333]*I + B(5.080416061514054781866889243246987462043762207031250000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*(4/7)^n, n >= 11))
If the GF is known analytic at a point, this can be passed to the algorithm
bound_coefficients(diff, ini, prec=20, known_analytic=[0,1])
1.00000*(1/2)^n*(([1.0000 +/- 8.44e-6] + [+/- 5.72e-14]*I)*n^(-1) + ([-1.0000 +/- 8.44e-6] + [+/- 5.72e-14]*I)*n^(-2) + ([1.000 +/- 1.96e-5] + [+/- 1.72e-13]*I)*n^(-3) + B([3068.05 +/- 7.82e-4]*n^(-4)*log(n), n >= 9))
This issue arises in lattice path enumeration.
For instance, walks on the steps $\{\text{North}, \text{South East}, \text{South West}\}$ restricted to the quadrant satisfy the following D-finite equation.
# D-finite equation and initial coeffs of GF
ini = [1,1,2,3,8,15,39]
deq = (196608*z^16 - 4096*z^15 - 106496*z^14 + 90112*z^13 + 4096*z^12 + 16192*z^11 + 3776*z^10 - 4736*z^9 - 560*z^8 - 124*z^7 - 33*z^6 + 52*z^5 + 7*z^4 - 2*z^3)*Dz^4 + (3145728*z^15 - 462848*z^14 - 1073152*z^13 + 1252352*z^12 + 192512*z^11 + 370112*z^10 + 76800*z^9 - 58176*z^8 - 11152*z^7 - 2468*z^6 - 968*z^5 + 518*z^4 + 90*z^3 - 22*z^2)*Dz^3 + (14155776*z^14 - 3870720*z^13 - 2199552*z^12 + 4733952*z^11 + 1582080*z^10 + 2059392*z^9 + 525696*z^8 - 164256*z^7 - 47904*z^6 - 9648*z^5 - 5508*z^4 + 1236*z^3 + 306*z^2 - 60*z)*Dz^2 + (18874368*z^13 - 7544832*z^12 + 294912*z^11 + 5081088*z^10 + 3201024*z^9 + 3148416*z^8 + 1068864*z^7 - 68832*z^6 - 41136*z^5 - 5520*z^4 - 7584*z^3 + 648*z^2 + 276*z - 36)*Dz + 4718592*z^12 - 2482176*z^11 + 811008*z^10 + 964608*z^9 + 1075200*z^8 + 880128*z^7 + 387648*z^6 + 31296*z^5 + 6864*z^4 + 3024*z^3 - 1320*z^2 + 72*z + 36
deq
(196608*z^16 - 4096*z^15 - 106496*z^14 + 90112*z^13 + 4096*z^12 + 16192*z^11 + 3776*z^10 - 4736*z^9 - 560*z^8 - 124*z^7 - 33*z^6 + 52*z^5 + 7*z^4 - 2*z^3)*Dz^4 + (3145728*z^15 - 462848*z^14 - 1073152*z^13 + 1252352*z^12 + 192512*z^11 + 370112*z^10 + 76800*z^9 - 58176*z^8 - 11152*z^7 - 2468*z^6 - 968*z^5 + 518*z^4 + 90*z^3 - 22*z^2)*Dz^3 + (14155776*z^14 - 3870720*z^13 - 2199552*z^12 + 4733952*z^11 + 1582080*z^10 + 2059392*z^9 + 525696*z^8 - 164256*z^7 - 47904*z^6 - 9648*z^5 - 5508*z^4 + 1236*z^3 + 306*z^2 - 60*z)*Dz^2 + (18874368*z^13 - 7544832*z^12 + 294912*z^11 + 5081088*z^10 + 3201024*z^9 + 3148416*z^8 + 1068864*z^7 - 68832*z^6 - 41136*z^5 - 5520*z^4 - 7584*z^3 + 648*z^2 + 276*z - 36)*Dz + 4718592*z^12 - 2482176*z^11 + 811008*z^10 + 964608*z^9 + 1075200*z^8 + 880128*z^7 + 387648*z^6 + 31296*z^5 + 6864*z^4 + 3024*z^3 - 1320*z^2 + 72*z + 36
# There is an apparent singularity at z=1/3 that is not an actual singularity
bound_coefficients(deq, ini, prec=20)
1.00000*3^n*(([+/- 8.87e-16] + [+/- 8.87e-16]*I)*n^(-1/2) + ([+/- 1.83e-15] + [+/- 1.83e-15]*I)*n^(-3/2) + ([+/- 4.97e-15] + [+/- 4.97e-15]*I)*n^(-5/2) + B([2.18961e+7 +/- 68.1]*n^(-7/2), n >= 10))
Solution: A new representation for the sequence!
Analytic Combinatorics in Several Variables (ACSV)¶
Consider now a $d$-variate series
$$ F(\mathbf{z}) = \sum_{\mathbf{i}\in\mathbb{N}^d} f_{\mathbf{i}}z_1^{i_1}\cdots z_d^{i_d} = \sum_{\mathbf{i}\in\mathbb{N}^d} f_{\mathbf{i}}\mathbf{z}^{\mathbf{i}}$$
The r-diagonal sequence consists of the coefficients $(f_{n\mathbf{r}}) = f_{\mathbf{0}}, f_\mathbf{r}, f_{2\mathbf{r}},\dots$
Note $f_{n\mathbf{r}}$ is defined only if $n\mathbf{r} \in \mathbb{N}^d$.
Assume that the rational function
$$F(\mathbf{z}) = \frac{G(\mathbf{z})}{H(\mathbf{z})} = \sum_{\mathbf{i}\in\mathbb{N}^d} f_\mathbf{i} \mathbf{z}^\mathbf{i} $$
converges in a neighbourhood of the origin.
The singularities of $F(\mathbf{z})$ are given by $\mathcal{V} = \{\mathbf{z}:H(\mathbf{z})=0\}$.
Singularities closest to the origin are called minimal points.
The Cauchy integral formula has a higher-dim generalization
$$ f_{n\mathbf{r}} = \frac{1}{(2\pi i)^d} \int_\mathcal{C} F(\mathbf{z}) \frac{d\mathbf{z}}{\mathbf{z}^{n\mathbf{r}+\mathbf{1}}}.$$
Simplest case: Denominator $H$ and partial derivatives don't simultaneously vanish (smooth simple poles).
In this case, the critical points of $F$ in the direction $\mathbf{r}$ are defined by
$$ H(\mathbf{z}) = 0, \frac{z_1}{r_1} H_{z_1}(\mathbf{z}) = \cdots = \frac{z_d}{r_d} H_{z_d}(\mathbf{z}) $$
Critical points: Locally the singular set has a saddle-point
Minimal points: Cauchy integral can be easily deformed close to
Generically: A finite set of minimal critical points.
Main Theorem of Smooth ACSV¶
Suppose that the critical point system admits a finite number of solutions.
If there is exactly one minimal solution $\mathbf{w} \in \mathbb{C}_*^d$ and $H_{z_d}(\mathbf{w})$ and a certain determinant is non-zero then
$$[\mathbf{z}^{n\mathbf{r}}]\frac{G(\mathbf{z})}{H(\mathbf{z})} = \mathbf{w}^{-n\mathbf{r}} n^{(1-d)/2} \left(C_\mathbf{w} + O\left(\frac{1}{n}\right)\right)$$
If there are a finite number of minimal critical points with the same coordinate-wise moduli satisfying these conditions, then just add the contributions.
Can show uniform error as $\mathbf{r}$ varies over generic directions, which also gives central limit theorems.
Consider the function
$$F(x,y) = \frac{1}{1-x-y} = \sum_{i,j \geq 0}\binom{i+j}{i}x^jy^j.$$
Critical point equations in direction $\mathbf{r} = (1,1)$
$$ 1-x-y = 0 \qquad -x=-y $$
Unique minimal critical point
$$(x,y)=\left(\frac{1}{2},\frac{1}{2}\right)$$
Asymptotics $$ \binom{2n}{n} = [x^ny^n]F(x,y) = \frac{4^n}{\sqrt{\pi n}}\left(1+O\left(n^{-1}\right)\right)$$
Consider the function
$$F(x,y) = \frac{1}{1-x-y} = \sum_{i,j \geq 0}\binom{i+j}{i}x^jy^j.$$
Critical point equations in direction $\mathbf{r} = (r,s)$
$$ 1-x-y = 0 \qquad -sx=-ry $$
Unique minimal critical point
$$(x,y)=\left(\frac{r}{r+s},\frac{s}{r+s}\right)$$
Asymptotics $$ \binom{rn+sn}{rn} = [x^{rn}y^{sn}]F(x,y) = \left(\frac{r+s}{r}\right)^{rn}\left(\frac{r+s}{s}\right)^{sn}\frac{\sqrt{r+s}}{\sqrt{2rs\pi n}}\left(1+O\left(n^{-1}\right)\right)$$
The hardest part of asymptotics under these conditions is verifying minimality.
This is much easier when $F$ is combinatorial, meaning it has only a finite number of negative series coefficients: think $O(2^{4d})$ vs $O(2^{12d})$.
M. and Salvy (2020): Complexity and symbolic-numeric algorithms for asymptotics under generic assumptions using a Kronecker representation for algebraic varieties. Preliminary Maple implementation.
Luo, Hackl, Selover, Wong, and M. (2022): Sage package with rigorous numerics in the smooth simple pole + combinatorial case.
All assumptions needed for asymptotics verified by package, except combinatoriality.
Software for Smooth ACSV¶
The sage-acsv package for ACSV.
Not official part of sage, but tracked on PyPi -- install with sage -pip install sage-acsv
Note: distinct from the built-in mvgf Sage package, which needs serious work
# Import code from Asymptotics.sage
from sage_acsv import *
# Make variables to be used later
var('x,y,w,z')
(x, y, w, z)
# Run central binomial coefficient example
F = 1/(1-x-y)
diagonal_asy(F, as_symbolic=True)
4^n/(sqrt(pi)*sqrt(n))
In general, asymptotics depend on (arbitrarily high degree) algebraic numbers.
# Sequence Alignment Example
F = (x^2*y^2-x*y+1)/(1-(x+y+x*y-x*y^2-x^2*y+x^2*y^3+x^3*y^2))
diagonal_asy(F, as_symbolic=True)
0.9430514023983397?*4.518911369262258?^n/(sqrt(pi)*sqrt(n))
diagonal_asy(F)
[(4.518911369262258?, 1/sqrt(n), 1/sqrt(pi), 0.9430514023983397?)]
a = diagonal_asy(F)[0][-1]
a.minpoly()
x^10 - x^8 + 5099/50272*x^6 - 17/6284*x^4 + 1/50272*x^2 - 1/804352
Apéry's proof of $\zeta(3)$ irrationality depends on bounding the exponential growth of the sequence
$$ \sum_{k=0}^n \binom{n}{k}^2\binom{n+k}{k}^2 = [w^nx^ny^nz^n]\frac{1}{1-w(1+x)(1+y)(1+z)(1+y+z+yz+xyz)}.$$
F = 1/(1 - w*(1 + x)*(1 + y)*(1 + z)*(x*y*z + y*z + y + z + 1))
diagonal_asy(F, as_symbolic=True)
1.225275868941647?*33.97056274847714?^n/(pi^(3/2)*n^(3/2))
# The quantities here are algebraic of degree two, so we can represent them in terms of radicals
var('n')
asm_vals = diagonal_asy(F)
show(add([a.radical_expression()^n*b*c*d.radical_expression() for (a,b,c,d) in asm_vals]))
The sequence $s_n$ counting lattice walks in the non-negative quadrant is the main diagonal of
$$ F(x,y,z) = \frac{(1+x)(1+y)}{1-zxy(1/x+x+1/y+y)}$$
# There are two minimal critical points determining asymptotics
F = (1+x)*(1+y)/(1-z*x*y*(1/x+x+1/y+y))
diagonal_asy(F)
[(-4, 1/n, 1/pi, 0), (4, 1/n, 1/pi, 4)]
Only the first term contributes to dominant asymptotics -- the second is $O\left((-4)^nn^{-2}\right)$.
diagonal_asy(F, as_symbolic=True)
4*4^n/(pi*n)
The algorithm uses a Kronecker Representation / RUR data structure, defined by
- square-free polynomial $P \in \mathbb{Z}[u]$
- rational functions $A,B,C \in \mathbb{Z}(u)$ with denominators $P'(u)$
- list $U$ of roots of $P$
Diagonal asymptotics are determined by
$$ (2 \pi r_d n)^{(1-d)/2} \sum_{u \in U} A(u) \sqrt{B(u)} C(u)^n\left(1+O\left(\frac{1}{n}\right)\right) $$
We can work numerically with sufficient precision to rigorously verify when there are minimal critical points.
Current algorithm uses Gröbner basis computations to compute this representation in Sage (slow!).
Work in progress to upgrade to something like msolve
The ability to compute Kronecker representations may be of independent interest (although it is currently the slowest part of most computations)
from sage_acsv import kronecker
var('x, y, z, t, lambda_')
kronecker([x*y*z - x*z + x - lambda_, x*y*z - y*z + y - lambda_,
x*y*z - x*z - y*z + z - lambda_,
x*y*z - x*z - y*z + x + y + z - 1,
x*y*z*t^3 - x*z*t^2 - y*z*t^2 + x*t + y*t + z*t -1],
[x, y, z, t, lambda_], linear_form = x + t)
(u_^8 - 18*u_^7 + 146*u_^6 - 692*u_^5 + 2067*u_^4 - 3922*u_^3 + 4553*u_^2 - 2925*u_ + 790, [10*u_^7 - 153*u_^6 + 1046*u_^5 - 4081*u_^4 + 9589*u_^3 - 13270*u_^2 + 9844*u_ - 2985, 10*u_^7 - 154*u_^6 + 1061*u_^5 - 4180*u_^4 + 9954*u_^3 - 14044*u_^2 + 10714*u_ - 3380, -u_^7 + 11*u_^6 - 56*u_^5 + 157*u_^4 - 182*u_^3 - 140*u_^2 + 527*u_ - 335, 8*u_^7 - 139*u_^6 + 1030*u_^5 - 4187*u_^4 + 10021*u_^3 - 14048*u_^2 + 10631*u_ - 3335, -12*u_^7 + 181*u_^6 - 1231*u_^5 + 4801*u_^4 - 11275*u_^3 + 15548*u_^2 - 11452*u_ + 3440])
Documentation is available in Sage
Source code hosted on https://pypi.org/project/sage-acsv/
help(kronecker)
Help on function kronecker in module sage_acsv.kronecker:
kronecker(system, vs, linear_form=None)
Computes the Kronecker Representation of a system of polynomials
INPUT:
* ``system`` -- A system of polynomials in ``d`` variables
* ``vs`` -- Variables of the system
* ``linear_form`` -- (Optional) A linear combination of the
input variables that separates the critical point solutions
OUTPUT:
A polynomial ``P`` and ``d`` polynomials ``Q1, ..., Q_d`` such that
``z_i = Q_i(u)/P'(u)`` for ``u`` ranging over the roots of ``P``.
Examples::
sage: from sage_acsv import kronecker
sage: var('x,y')
(x, y)
sage: kronecker([x**3+y**3-10, y**2-2], [x,y], x+y)
(u_^6 - 6*u_^4 - 20*u_^3 + 36*u_^2 - 120*u_ + 100,
[60*u_^3 - 72*u_^2 + 360*u_ - 600, 12*u_^4 - 72*u_^2 + 240*u_])
If you have a non-combinatorial function, computations are much more expensive.
Only practical package uses Homotopy Continuation in Julia...
More details can be found in these books
