Help with deriving information from a generating function

Regarding how the formulas like the one on the OEIS can be found:

A sequence will always grow with exponential order $R^{-n}$, where $R$ is the absolute value of the singularity of $A$ closest to the origin. In the $c_0 = c_1 = c_2 = c_3 = 1$ case, Wolfram Alpha gives this for $A(x)$. Taking the root of the expression under the square root (as there are no other singularities), i.e. solving $$16x^{12}+8x^{11}+11x^{10}-4x^9=0$$ yields a value which has a reciprocal with the same minimal polynomial as found for $d$ on OEIS. The $3/2$ exponent on $n$ in the asymptotic formula results from the fact the singularity comes from a square root.

A similar process can probably be done for the general case, although the results look fairly horrific very quickly. If the dominant singularity is a pole, the subexponential factor is a polynomial of degree one less than the order, instead of being $n^{-3/2}$.

Specifics are much better explained by this reference from a similar question.

An asymptotic follows from Theorem VI.6 ("Singular Inversion") in Flajolet and Sedgewick's Analytic Combinatorics (freely legally available online at the link). The theorem concerns generating functions which solve equations of the form $Y(z) = z \phi(Y(z))$, where $\phi$ is a polynomial or more generally a power series; these count trees of various sorts. Setting $Y(z) = z A(z)$ lets us put your problem in this framework, with $\phi(u) = c_0 + c_1 u + c_2 u^2 + c_3 u^3$. We need $\phi(u)$ to satisfy three sets of conditions (equations (35) and (36) on pages 402 and 403, together with a third condition detailed on page 404):

• $\phi(0) \neq 0, [u^n] \phi(u) \ge 0, \phi(u) \neq \phi_0 + \phi_1 u$, and
• within the open disc of convergence of $\phi$, there is a (necessarily unique) positive real $\tau$ such that $\phi(\tau) - \tau \phi'(\tau) = 0$ (the characteristic equation).
• the exponents of the nonzero terms of $\phi(u)$ are not all divisible by some $p \ge 2$ ($\phi(u)$ is aperiodic).

Then the radius of convergence of $Y(z)$ is $\rho = \frac{\tau}{\phi(\tau)} = \frac{1}{\phi'(\tau)}$. For example, when $c_i = 1$ the characteristic equation is

$$(\tau^3 + \tau^2 + \tau + 1) - \tau (3 \tau^2 + 2 \tau + 1) = - 2 \tau^3 - \tau^2 + 1 = 0$$

which has a unique positive real root $\tau \approx 0.657$ (an exact value can be written down using the cubic formula) and which gives $\rho \approx 0.277$.

Theorem: With the above hypotheses, $y_n = [z^n] Y(z)$ has an symptotic expansion beginning

$$y_n \sim \sqrt{ \frac{\phi(\tau)}{2 \phi''(\tau)} } \frac{\rho^{-n}}{\sqrt{\pi n^3}}.$$

Again with $c_i = 1$ we have $\rho^{-1} \approx 3.61 \dots$ so this is the $d$ which appears in the result you cited from the OEIS, and the constant factor is

$$\sqrt{ \frac{\tau^3 + \tau^2 + \tau + 1}{(12 \tau + 4)\pi} } \approx 0.252$$

but since $Y(z) = z A(z)$ our indexing is off by $1$ and $y_n = a_{n-1}$. So to get $a_n = y_{n+1}$ we multiply by $\rho^{-1}$ and get a constant factor $\approx 0.910$ which matches up with the result you cited from the OEIS again. Using the cubic formula and a bunch of algebraic manipulations it would be possible to match the exact results using radicals but it doesn't sound fun.

All of these computations simplify substantially if $c_3 = 0$, since then all of the cubic equations become quadratic. In this case we can write down an explicit closed form for $Y(z)$ using the quadratic formula and analyze that closed form using the binomial theorem to get asymptotics. This case is an elaboration on that one; the asymptotics come from a square-root singularity although it's harder to tell that that's what's happening.

---

The first formula you give can be obtained from Lagrange inversion, which gives that if $\frac{Y(z)}{\phi(Y(z))} = z$ then $Y$ (which is the compositional inverse of $\frac{w}{\phi(w)}$) has coefficients given by

$$[z^n] Y(z) = \frac{1}{n} [w^{n-1}] \phi(w)^n.$$

When $c_i = 1$ this gives

$$\begin{eqnarray*} a_n &=& \frac{1}{n+1} [w^n] (1 + w + w^2 + w^3)^{n+1} \\ &=& \frac{1}{n+1} [w^n] (1 + w)^{n+1} (1 + w^2)^{n+1} \\ &=& \frac{1}{n+1} [w^n] \left( \sum_{i=0}^{n+1} {n+1 \choose i} w^i \right) \left( \sum_{j=0}^{n+1} {n+1 \choose j} w^{2j} \right) \\ &=& \frac{1}{n+1} \left( \sum_{i+2j = n} {n+1 \choose i} {n+1 \choose j} \right) \end{eqnarray*}$$

which is equivalent to the desired result. A nicer special case is $c_0 = c_2 = 1, c_1 = c_3 = 0$, or $\phi(w) = 1 + w^2$, which corresponds to binary trees and gives $a_{2n+1} = 0$ and

$$a_{2n} = \frac{1}{2n+1} [w^{2n}] (1 + w^2)^{2n+1} = \frac{1}{2n+1} {2n+1 \choose n} = \frac{1}{n+1} {2n \choose n}$$

which is the Catalan numbers. These have a similar asymptotic expansion because their generating function also has a square-root singularity, but the base of the exponent is $4$. It would be good to understand this case as a warmup exercise if you haven't already done so.