# How to find the asymptotics of a linear two-dimensional recurrence relation

Let $$d$$ be a positive number. There is a two dimensional recurrence relation as follow: $$R(n,m) = R(n-1,m-1) + R(n,m-d)$$ where $$R(0,m) = 1$$ and $$R(n,0) = R(n,1) = \cdots = R(n, d-1) = 1$$ for all $$n,m>0$$.

How to analyze the asymptotics of $$R(n, kn)$$ for fixed $$k$$? It is easy to see that $$R(n, kn) = O\left( c_{k,d}^{n} \cdot (n+k+d)^{O(1)} \right)$$ Is there a way (or an algorithm) to find $$c_{k,d}$$ given $$k$$ and $$d$$?

PS: I have calcuated the bivariate generating function of $$R(\cdot, \cdot)$$: \begin{align} f(x,y) &= \frac{1 - xy - y^{d} + xy^{d}}{(1 - x)(1 - y)(1 - xy - y^{d})} \\ &= \frac{1}{(1 - x)(1 - y)} + \frac{xy^{d}}{(1 - x)(1 - y)(1 - xy - y^{d})} \\ \end{align}

• For $k=1$ the gf seems to be $\frac1{1-t-t^d}$, right? Feb 16 at 8:51
• @მამუკაჯიბლაძე You mean the generating function of $T(n) = R(n, n)$? I am not sure. Feb 16 at 8:58
• @მამუკაჯიბლაძე Yes, that is right, there was a question when $d=1$, math.stackexchange.com/questions/2065067/… Feb 16 at 9:05
• If you add $y-y$ to the numerator of $f$ you can split it into $$f(x,y) = \frac{1}{(1 - x)(1 - xy - y^d)} + \frac{\sum_{i=1}^{d-1} y^i}{1 - xy - y^d}$$ The second term is easy to extract coefficients from, but the first one is messy. Feb 16 at 12:03
• More generally,$$R(n,m)=\sum_{j=0}^n\binom{\left[\frac{m-j+1}d\right]+j-1}j$$ Feb 16 at 14:31

Your generating function $$f(x,y)$$ is convergent on a polydisk $$|x|<\epsilon$$, $$|y|< \delta$$, for some $$\epsilon, \delta < 1$$. We can reduce $$\epsilon$$ so that $$\epsilon \ll \delta^k$$ and the domain of $$f(x,y)$$ to the product of the polydisk $$|x|<\epsilon$$ and the annulus $$\delta' < |y| < \delta$$, with $$\delta' \sim \delta$$. Then $$g(z,y) = f(z/y^k, y)$$ is well-defined and analytic on the product domain of $$|z|\lesssim \epsilon/\delta^k$$ and $$\delta' < |y| < \delta$$. The generating function $$g(z,y)$$ has a unique expansion that is a power series in $$z$$ and a Laurent series in $$y$$. Since $$x^n y^m = (z/y^k)^n y^m = z^n y^{m-kn}$$, the coefficients $$R(n,kn)$$ are the coefficients of $$z^n y^0$$ in this double expansion. Hence, using the Cauchy integral formula, the generating function for $$R(n,kn)$$ is now given by the contour integral $$\sum_{n=0}^\infty R(n,kn) z^n = \frac{1}{2\pi i} \oint_{|y|\sim \delta} g(z,y) \frac{dy}{y} .$$ Of course, the integration contour is now allowed to deform, as long as it doesn't hit any singularities of the integrand.

The above integrand has a convenient partial fraction expansion expansion with respect to $$z$$, $$\frac{g(z,y)}{y} = \frac{y^k}{(1-y-y^d)(y^k-z)} - \frac{y^{d+k-1}(1-y^d)}{(1-y)(1-y-y^d) (y^{k-1} (1-y^d)-z)} .$$ Expanding in $$z$$ under the integral gives the formula $$R(n,kn) = \frac{1}{2\pi i} \oint \frac{dy}{y^{kn} (1-y-y^d)} - \frac{1}{2\pi i} \oint \frac{y^d dy}{(1-y)(1-y-y^d) y^{n(k-1)} (1-y^d)^n} ,$$ whose leading asymptotics can be estimated by residues or steepest descent.

Using residues on the first term, all the contributions come from the roots of $$1-y-y^d=0$$. The root with the smallest magnitude $$y_*$$ will give the leading contribution. Some experimentation shows that $$y_*$$ is the unique positive real root. The leading asymptotic term then looks like $$C_* y_*^{-kn} .$$ In the special case $$d=1$$, $$y_* = 1/2$$. In the special case $$k=1$$, the second term has no poles inside its contour and hence evaluates to zero. So in the sequel we can assume that $$k>1$$.

Applying steepest descent to the second integral, we find the stationary phase points at the roots $$y_\star$$ of $$y^d = (k-1)/(k+d-1) + O(1/n)$$. Let $$y_\star$$ be that root (actually we just need its $$n\to \infty$$ limit) which minimizes the magnitude of $$y_\star^{k-1} (1-y_\star^d)$$. It is easy too see that $$y_\star = \sqrt[d]{(k-1)/(d+k-1)}$$. Then the leading asymptotic term is $$C_\star n^{-1/2} y_\star^{-n(k-1)} (1-y_\star^d)^{-n} = C_\star n^{-1/2} w_\star^{-kn} ,$$ where $$w_\star = \left(\frac{k-1}{d+k-1}\right)^{\frac{1-1/k}{d}} \left(\frac{d}{d+k-1}\right)^{1/k} = \frac{s^{s/k}}{(s+1)^{(s+1)/k}}$$ with $$s=(k-1)/d$$, as pointed out in the comments.

The conclusion is that $$R(n,kn) = O\left((\min(y_*,w_\star)^{-k})^n\right) .$$

Experimentation suggests that $$y_* < w_\star$$ in all cases, except for $$(k,d)=(2,1)$$, when $$y_* = w_\star = 1/2$$. Otherwise the $$y_*$$ contribution always seems to dominate over the $$w_\star$$ contribution. At least that's how it seemed for not too large values of $$d$$.

Doing some rudimentary asymptotic calculations for large $$d$$, it seems that $$y_* \sim 1- \frac{\log d}{d} + o(\log(d)/d)$$. A better approximation is $$y_* \sim 1 - \frac{W(d)}{d} + \frac{W(d)^3}{2(W(d)+1) d^2} + O(W(d)^3/d^3) ,$$ where $$W(d)$$ is the Lambert W function. On the other hand, $$w_\star$$ has a minimum as a function of $$k$$ (treating it as a continuous variable) around $$k_\star \sim W(d) + 1 + O(W(d)^2/d)$$, the minimum reaches roughly $$w_\star \sim 1 - \frac{W(d)}{d-1} \lesssim y_* .$$ One should push the asymptotics of $$k_\star$$ one more order to get a better estimate of the difference, but experiments do show that the minimum of $$w_\star$$ does dip below $$y_*$$, for instance this happens near $$(k,d) = (5, 205)$$ or $$(6,700)$$. So it looks like the $$w_\star$$ asymptotic contribution will dominate in small ranges of $$k \sim k_\star(d)$$ if it happens to fall near an integer. To get the size of that window, one could go to a quadratic Taylor approximation of $$w_\star$$.

With a bit of extra work, one could also extract the coefficients $$C_*$$ and $$C_\star$$ from the integral formula.

• Thanks, I mean the asymptotics is $C^{n} \cdot n^{O(1)}$. I am considering the exact expression of $C$ which depends on $k$ and $d$ only. Feb 16 at 15:07
• @Blanco I see, my bad. Sorry, I missed the exponent $n$ on your $c_{k,d}^n$. Then this approach suggests that the largest of $y_*^{-k}$ and $y_\star^{-k} (1-y_\star^d)^{-1}$ is the base of the leading exponential asymptotic. Feb 16 at 16:20
• Thanks a lot, very well explained. The root of $1-y-y^d$ makes sense, since it is the characteristic equation of the linear 1-dimensional recurrence relation of $R(n) = R(n-1) + R(n-d)$. As for the second root, how to caluate it, does it mean $y_{\star} = \sqrt[d]{(k-1)/(k+d-1)}$? Feb 17 at 6:22
• @Blanco Yes, that's right. Feb 17 at 7:39
• According to Iosif Pinelis's and my results, the second base should be $y_{\star}^{-(k-1)} \cdot (1 - y_{\star}^d)^{-1}$ which is equal to $(s+1)^{s+1}/s^s$ where $s = (k-1)/d$. Do I miss something? Feb 17 at 8:22

This is to complement Blanco's answer by showing that $$\begin{equation*} R(n,kn)=\exp\{(C_{k,d}+o(1))\,n\} \tag{0}\label{0} \end{equation*}$$ (as $$n\to\infty$$), where $$k\ge1$$ and $$d\ge1$$ (are fixed), $$\begin{equation*} C_{k,d}:=\frac k{1+y_{k,d}\,d}\,\big(\ln(1+y_{k,d})+y_{k,d}\,\ln(1+1/y_{k,d})\big), \end{equation*}$$ $$\begin{equation*} y_{k,d}:=\max\Big(x_d,\frac{k-1}d\Big), \end{equation*}$$ and $$x_d$$ is the unique positive root of the equation $$\begin{equation*} x_d(1+x_d)^{d-1}=1. \tag{0.5}\label{0.5} \end{equation*}$$

In view of Blanco's answer, it is enough to show that $$\begin{equation*} B(n):=B(n,kn)=\exp\{(C_{k,d}+o(1))\,n\}, \tag{1}\label{1} \end{equation*}$$ where $$\begin{equation*} B(n):=\max\{c_{a,b}\colon (a,b)\in E_{n,k,d}\}, \end{equation*}$$ $$\begin{equation*} E_{n,k,d}:=\{(a,b)\colon 0\le a\le n-1,b\ge0,bd+a\le kn-d,\ a,b \text{ are integers}\}, \end{equation*}$$ $$\begin{equation*} c_{a,b}:=\binom{a+b}b=\binom{a+b}a. \end{equation*}$$ Note that $$c_{a,b}$$ is increasing in $$a$$ and in $$b$$.

Note also that $$c_{a,b}=(a+b)^{O(1)}=n^{O(1)}$$ for $$(a,b)\in E_{n,k,d}$$ if $$a=O(1)$$ or $$b=O(1)$$.

So, it remains to consider the case when $$a\to\infty$$ and $$b\to\infty$$. Then, by Stirling's formula, \begin{equation*} \begin{aligned} \ln c_{a,b}\sim a\ln(1+b/a)+b\ln(1+a/b). \end{aligned} \tag{2}\label{2} \end{equation*}

Also, \begin{equation*} \begin{aligned} \frac{c_{a+d,b-1}}{c_{a,b}}&=b\frac{(a+d+1)\cdots(a+d+b-1)}{(a+1)\cdots(a+b)} \\ &= b\frac{(a+b+1)\cdots(a+b+d-1)}{(a+1)\cdots(a+d)} \\ &\sim b\frac{(a+b)^{d-1}}{a^d} =\frac ba\Big(1+\frac ba\Big)^{d-1}. \end{aligned} \end{equation*} So, for a fixed value of $$bd+a$$, the maximum of $$c_{a,b}$$ occurs when $$b/a\to x_d$$ (recall \eqref{0.5}).

Let now $$(a,b)\in E_{n,k,d}$$ be a maximizer of $$c_{a,b}$$ (such that $$a\to\infty$$ and $$b\to\infty$$). Then, since $$c_{a,b}$$ is increasing in $$a$$ and in $$b$$, we have $$\begin{equation*} bd+a\sim kn. \end{equation*}$$

The conditions $$b/a\to x_d$$ and $$bd+a\sim kn$$ imply $$\begin{equation*} a\sim k\frac1{1+x_d\,d}\,n, \tag{3}\label{3} \end{equation*}$$ and the latter condition is compatible with condition $$a\le n-1$$ only if $$\begin{equation*} k\frac1{1+x_d\,d}\le1. \tag{4}\label{4} \end{equation*}$$

If this is the case, then $$\begin{equation*} b\sim k\frac{x_d}{1+x_d\,d}\,n, \end{equation*}$$ so that, by \eqref{2}, \begin{equation*} \begin{aligned} \frac{\ln c_{a,b}}n\to \frac k{1+x_d\,d}\,(\ln(1+x_d)+x_d\,\ln(1+1/x_d)). \end{aligned} \tag{5}\label{5} \end{equation*} Also, \eqref{4} implies $$x_d\ge(k-1)/d$$, so that $$y_{k,d}=x_d$$. So, we have proved \eqref{1}, and thus \eqref{0} -- in the case when condition (3) is compatible with condition $$a\le n-1$$.

Otherwise, we have $$a=n-1\sim n$$ and still $$bd+a\sim kn$$, whence $$b\sim(k-1)n/d$$. So, by \eqref{2}, here \begin{equation*} \begin{aligned} \frac{\ln c_{a,b}}n\to \ln\Big(1+\frac{k-1}d\Big)+\frac{k-1}d\,\ln\Big(1+\frac d{k-1}\Big). \end{aligned} \tag{6}\label{6} \end{equation*} Also, in this "incompatibility" case, we have $$k\frac1{1+x_d\,d}\ge1$$ -- cf. \eqref{4}. So, here $$x_d\le(k-1)/d$$ and hence $$y_{k,d}=\frac{k-1}d$$. So, we have proved \eqref{1}, and thus \eqref{0} -- in the case when condition (3) is incompatible with condition $$a\le n-1$$.

Thus, in either case, \eqref{0} is proved. $$\quad\Box$$

The book Analytic Combinatorics in Several Variables by Pemantle and Wilson covers problems like this extensively.

This is a part of my approach, and I have not finished.

Firstly, $$\frac{1}{(1-x)(1-y)} = \sum_{i,j \geq 0} x^{i}y^{j}$$ and $$\frac{1}{1 - xy - y^{d}} = \sum_{i \geq 0} \frac{x^{i}y^{i}}{(1 - y^{d})^{i + 1}} = \sum_{i, j \geq 0} \binom{i + j}{j} x^{i}y^{jd + i}$$

Thus \begin{align} \frac{1}{(1 - x)(1 - y)(1 - xy - y^{d})} = &\left( \sum_{i,j \geq 0} x^{i}y^{j} \right) \cdot \left( \sum_{i, j \geq 0} \binom{i + j}{j} x^{i}y^{jd + i} \right) \\ = &\sum_{i,j \geq 0} \sum_{a \leq i, bd + a \leq j} \binom{a + b}{b} x^{i}y^{j} \\ \end{align} We have $$f(x, y) = \sum_{i,j \geq 0} \left( 1 + \sum_{a \leq i - 1, bd + a \leq j - d} \binom{a + b}{b} \right) x^{i}y^{j}$$ and $$R(n, m) = 1 + \sum_{a + 1 \leq n, (b + 1)d + a \leq m} \binom{a + b}{b}$$

Let $$B(n,m) = \max_{a + 1 \leq n, (b + 1)d + a \leq m} \binom{a + b}{b}$$ then we have $$B(n,m) \leq R(n,m) \leq nmB(n,m)$$ Thus $$R(n, m) = O(B(n,m))\cdot (nm)^{O(1)}$$