computing $n^{th}$ Schröder number

Question

We can compute $n^{th}$ Schröder number in $O(n^2)$ using the recurrence relation:

$S_{n+3} = 3 S_{n+2} + \sum_{k=0}^n S_{k+1} * S_{n-k+1}$.

My question is how to compute them faster?

Answer 1

The Schroder numbers $S_n$ can be expressed in terms of Legendre polynomials evaluated at $3$ by the following formula \begin{eqnarray*} S_n= \frac{1}{2} (-P_{n-1}(3)+6P_n(3)-P_{n+1}(3)). \end{eqnarray*} Legendre polynomials satisfy the recurrence relation \begin{eqnarray*} (n+1)P_{n+1}(x)=(2n+1)xP_n(x)-nP_{n-1}(x). \end{eqnarray*} Substitute $x=3$ in this formaula and generate the sequence $P_n(3)$ (this gives the sequence $1,3,13,63,\cdots$ which are the central Delannoy numbers.) And now use the first equation to compute the Schroder numbers $1,2,6,22,90,\cdots$.

Answer 2

Delannoy numbers can be defined as $$ D_n = [x^n y^n]\frac{1}{1-x-y-xy}=[x^n y^n]\frac{1}{2-(1+x)(1+y)}=\sum_{k\geq 0}\frac{1}{2^{k+1}}\binom{k}{n}^2 = P_n(3)\tag{1} $$ and due to the orthogonality relations in $L^2(0,2\pi)$ the RHS of $(1)$ can be represented as $$ D_n = P_n(3) = \frac{1}{2\pi}\int_{0}^{2\pi}\frac{d\theta}{(3-2\sqrt{2}\cos\theta)^{n+1}} \tag{2} $$ giving both the generating function and the asymptotic behaviour of Delannoy numbers: $$ \sum_{n\geq 0}D_n x^n = \frac{1}{\sqrt{1-6x+x^2}},\qquad D_n\sim\frac{C}{\sqrt{n}}(1+\sqrt{2})^{2n}. \tag{3} $$ Schroeder numbers are related to Delannoy numbers via $S_n = \frac{1}{2}\left(-D_{n-1}+6\,D_n-D_{n+1}\right)$, hence they have a similar generating function, a similar asymptotic behaviour and they can be computed by tackling the integral $\frac{1}{2\pi}\int_{0}^{2\pi}\frac{d\theta}{(3-2\sqrt{2}\cos\theta)^{n}}$ through Gaussian quadrature and rounding to the closest integer.

The existence of an algorithm by repeat squaring (similar to the one allowing a fast computation of Fibonacci or Lucas numbers) does not seem to be investigated in the literature (or, at least, I have not been able to find it) but it is for sure an interesting question, which deserves further thoughts.

Here it is one: for any $R\in(-1,1)$ (in particular for $R=\frac{2\sqrt{2}}{3}$) we have the following identity: $$ \frac{1}{1-R\cos\theta} = \frac{1}{\sqrt{1-R^2}}\sum_{h\geq 0}\left(\frac{R}{1+\sqrt{1-R^2}}\right)^h \cos(h\theta)\tag{4} $$ i.e. we have a Fourier (cosine) series whose coefficients have an exponential decay. A similar identity holds for $\frac{1}{(1-R\cos\theta)^2}$, and given the Fourier cosine series of $\frac{1}{(1-R\cos\theta)^n}$ and $\frac{1}{(1-R\cos\theta)^{n+1}}$, to compute the Fourier cosine series of $\frac{1}{(1-R\cos\theta)^{2n}}$ and $\frac{1}{(1-R\cos\theta)^{2n+1}}$ is simple by convolution. In particular the RHS of $(2)$ can be computed through a repeat-squaring approach: the hardest part of this approach is to understand where the exponential decay allows to truncate the involved Fourier cosine series.

This is also essentially equivalent to understanding how squaring works in the binomial base.
A concrete example is given by the computation of $D_5$:

$$ D_5=\sum_{k\geq 0}\frac{\binom{k}{5}^2}{2^{k+1}}=\sum_{k\geq 0}\frac{a_{10}\binom{k}{10}+a_9\binom{k}{9}+\ldots+a_1\binom{k}{1}+a_0}{2^{k+1}}=a_{10}+a_9+\ldots+a_1+a_0 $$ where we also know that $D_5=\sum_{k=0}^{5}\binom{5}{k}^2 2^{5-k}$.