\(\(\(
\)\( (*\ 
      This\ Mathematica\ notebook\ computes\ the\ asymptotic\ number\ of\ pi\ shapes\ which\ are\ compatible\ with\ an\ RNA\ sequence\ of\ length\ n, \ 
      under\ the\ assumption\ that\ any\ base\ can\ basepair\ with\ any\ other\ base\ and\ that\ there\ are\ at\ least\ 3\ unpaired\ bases\ in\ a\ hairpin\ \(\(loop\)\(.\)\)\ *) \)\(
\)\(
\)\(
\)\(Clear["*"];\)\)\)

eqn = {S == z*S + A,  
        A == z^2*A*B + z^2*B, 
B == 
          z^2*A*B + z^3}

solset = Simplify[Solve[eqn,  S,  {A,  B}]]

Here are the 2 possibilities for the generating function.  We choose the one that does not blow up at 0.

sol = S /. Part[solset,  2,  1]

The dominant singularity will occur where the square root evaluates to 0, or 
where the denominator is 0 (z = 1,  but not z = 0, actually analytic there),
whichever is closest to origin.

poly = 1 - 2 (z^5) - 4 (z^7) + (z^10)

solset = NSolve[poly == 0,  z]

r = z /. Part[solset,  9,  1]

This gives us our dominant singularity, r.

1/r

This is our exponential growth rate.

p2a = poly/(1 - z/r)

Series[p2a,  {z,  0,  10}]

This is how we have figured to divide polynomials in Mathematica.  Since r is only a numerical approximation, p2a is not actually a polynomial gotten by dividing out one root, but it is really close to one.  By cutting and pasting the above, we can get a polynomial p2 such that p2*(1-z/r) is very close to poly.  That is all we need.

p2 = 
      1 + 1.3221783432637244` z + 1.748155571395607` (z^2) + 
        2.311373437155093` (z^3) + 3.0560479018015014` (z^4) + 
        2.0406403517384906` (z^5) + 2.6980904794587004` (z^6) - 
        0.4326431998936666` (z^7) - 0.5720314692597253` (z^8) - 
        0.7563276203205369` (z^9)

Expand[p2*(1 - z/r)]

(Very close to p1, so p2 is what we want.)

sol

left = Sqrt[p2]/(2 (-1 + z)*z^2)

Remember, only the part that was next to the dominant singularity will matter.

K = 1/Gamma[-1/2]*left /. z -> r

1/r

Remember alpha=1/2 as the singularity was (1-z/r)^(1/2), and that exponent gives us alpha.  We then know we have a factor of n^(-alpha-1)

So S_n ~ 2.44251*1.32218^n * n^(-3/2)