
# Define the rings of differential and shift operators
from ore_algebra import *
Ind.<n> = PolynomialRing(QQ); Shift.<Sn> = OreAlgebra(Ind)
Pols.<z> = PolynomialRing(QQ); Diff.<Dz> = OreAlgebra(Pols)


# Define differential equation for the generating function
diff = z^2*(4*z-1)*(4*z+1)*Dz^3 + 2*z*(4*z+1)*(16*z-3)*Dz^2 + 2*(112*z^2 + 14*z - 3)*Dz + 64*z + 12
diff

(16*z^4 - z^2)*Dz^3 + (128*z^3 + 8*z^2 - 6*z)*Dz^2 + (224*z^2 + 28*z - 6)*Dz + 64*z + 12


# Convert to recurrence for the sequence a_n
rec = diff.to_S(Shift).primitive_part()
rec

(-n^2 - 7*n - 12)*Sn^2 + (8*n + 20)*Sn + 16*n^2 + 48*n + 32


# There is a basis of solutions to this recurrence, whose asymptotic expansions as n->infinity begin
rec.generalized_series_solutions()

[4^n*n^(-1)*(1 - 3/2*n^(-1) + 19/8*n^(-2) - 63/16*n^(-3) + 871/128*n^(-4) + O(n^(-5))),
 (-4)^n*n^(-3)*(1 - 9/2*n^(-1) + 107/8*n^(-2) - 525/16*n^(-3) + 9263/128*n^(-4) + O(n^(-5)))]


# We can use the recurrence to generate terms of the counting sequence, from the initial values a_0 = 1 and a_1 = 2
rec.to_list([1,2],20)

[1,
 2,
 6,
 18,
 60,
 200,
 700,
 2450,
 8820,
 31752,
 116424,
 426888,
 1585584,
 5889312,
 22084920,
 82818450,
 312869700,
 1181952200,
 4491418360,
 17067389768]


# There is a basis of solutions to the D-finite equation defined by diff whose series expansions near the origin begin
diff.generalized_series_solutions(5)

[1 + 2*z + 6*z^2 + 18*z^3 + 60*z^4 + O(z^5),
 z^(-1)*(1 + 2*z + 4*z^2 + 12*z^3 + 36*z^4 + O(z^5)),
 z^(-2)*((z + 2*z^2 + 4*z^3 + 12*z^4 + O(z^5))*log(z) - 1/4 - 3/2*z - 2*z^2 - 2*z^3 - 5*z^4 + O(z^5))]

Example 2.23 (D-Finite and P-Recursive Bases of Solutions)¶