(* This Mathematica notebook computes the asymptotic number of pi shapes of length n *)

Clear["*"] ;

*All comments begin with a * for clarity

eqn = {S == z^2 * S * T + z^2 * T, T == z^2 * S * T + 1} ;

*These equations are for (nonempty) pi shapes

solset = Simplify[Solve[eqn, S, T]]

{{S→ -(-1 + z^2 + (1 - 2 z^2 - 3 z^4)^(1/2))/(2 z^2)}, {S→ (1 - z^2 + (1 - 2 z^2 - 3 z^4)^(1/2))/(2 z^2)}}

*Here is the formula for the power series.  We choose the one that does not explode at z=0 (the one whose numerator is 0 at z=0.

sol = S/. Extract[solset, 1]

-(-1 + z^2 + (1 - 2 z^2 - 3 z^4)^(1/2))/(2 z^2)

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

z^2 + z^4 + 2 z^6 + 4 z^8 + 9 z^10 + O[z]^11

sol /.z→0.00001

0.

*The dominant singularity will occur where the square root evaluates to 0 or where the denominator is 0.  (However, a generating function is always analytic at z=0.)  So we must find roots of that polynomial.

poly = 1 - 2 z^2 - 3 z^4

1 - 2 z^2 - 3 z^4

solset = NSolve[poly == 0, z]

{{z→ -0.57735}, {z→0. - 1. }, {z→0. + 1. }, {z→0.57735}}

*We immediately get the exponential growth as 1/0.57735=1.73, as this is one of the dominant singularities.

growth = 1/z/.Extract[solset, 4]

1.73205

*However, we see a problem.  There is no unique dominant singularity.  There are 2 singularities at z=±0.57735.  This is related to the fact that the coefficients of all odd terms are 0.  Thus technically this series doesn't have an asymptotic limit.  We will only consider even terms, make the substitution w=z^2, and run through the identical analysis.  (Note, this series can be considered the sum of 2 series that are identical, except that one oscillates between positive and negative.  These correspond to the 2 dominant singularities.  If we don't make substitution, blindly applying our Flajolet-Odlyzko method or alternatively applying Bender's theorem will pick up the positive dominant singularity, and we will only be off by a constant factor of 2.)

solw = sol/.z→Sqrt[w]

-(-1 + w + (1 - 2 w - 3 w^2)^(1/2))/(2 w)

polyw = 1 - 2w - 3w^2

1 - 2 w - 3 w^2

solset = Solve[polyw == 0, w]

{{w→ -1}, {w→1/3}}

r = w/.Extract[solset, 2]

1/3

*This gives us our dominant singularity, r.  The exponential growth is 1/r=3.

growth = 1/r

3

p2 = Simplify[polyw/(1 - w/r)]

1 + w

left = -Sqrt[p2]/(2 * w)

-(1 + w)^(1/2)/(2 w)

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

K = left/.w→r

-3^(1/2)

K2 = K/Gamma[-1/2]

3/π^(1/2)/2

*So this is Gamma[-α] where α=1/2 is because the singularity is of the form (1 - w/ρ)^(1/2).

*Now we must convert back to z.  First state results in terms of S.  Since α=1/2, we know we will have an n^(-α-1) or an n^(-3/2) term.

*So we find S_2n ∼K2 (3)^n (n)^(-3/2).  This is in terms of w, that's where the 2n subscript comes from.  So, to get S_n, we need merely substitute (n/2) for n, and get
S_n=Knew (Sqrt[3])^n (n)^(-3/2).

Knew = K2 * 2^(3/2)

6/π^(1/2)

Created by Mathematica  (December 22, 2006) Valid XHTML 1.1!