o.g.f for linear homogeneous recurrences with constant coefficients for Maple

[Max Alekseyev]

On Nov 28, 2007 1:44 PM, Richard Mathar <mathar at strw.leidenuniv.nl> wrote:
>
> A Maple implementation of calculating the ordinary generating function
> for sequences with linear homogeneous recurrences with constant coefficients
> (and some very simple inhomogeneous cases) is in
> http://www.strw.leidenuniv.nl/~mathar/progs/GenFLinRec.mp .
> A short introduction on this is in
> http://www.strw.leidenuniv.nl/~mathar/public/mathar20071126.pdf .
> No new mathematics is involved; it is merely useful if one wishes to compile
> ordinary generating function on some industrial scale.

Richard,

I did not look thoroughly but it looks like your routines can be
optimized, at least in the homogeneous case. It is known that the
o.g.f. of linear recurrent homogeneous sequence can be easily derived
from the generator and characteristic polynomial of this sequence.
This is my PARI/GP code for two procedures: one computes the generator
from the given characteristic polynomial f and vector u of initial
terms (the length of u must be equal to the degree of f), and the
other computes the o.g.f.:

\\ Given a characteristic polynomial f and an initial vector u,
compute the generator
{ lrgen(f,u) = local(m=#u);
  if( m!=poldegree(f), error("length(u) != deg(f)") );
  sum(k=0,m-1,sum(i=0,k,u[k+1-i]*polcoeff(f,m-i))*Pol([1,0])^(m-1-k))
}

\\ o.g.f. is x*lrgen(f,u)/f with follow-up substitution x -> x^(-1)
{ lrogf(f,u) = subst( x*lrgen(f,u)/f, x, 1/x) }

For example, for Fibonacci sequence we have:

? lrogf(x^2-x-1,[0,1])
%1 = x/(-x^2 - x + 1)

And this is a wrapper that offers an interface similar to your
GenFLinRec() function (but only first 3 arguments are supported):

{ GenFLinRec(initL, recurL, off) = x^off * lrogf(
Pol(concat([[-1],vector(#initL-#recurL,i,0),recurL])), initL) }

Homogeneous examples from your code work well:

> # Example 3:
> # GenFLinRec( [1,2,2], [0,0,1], 0, 0, true,[]) returns the generating
> # function specified by a(0)=1, a(1)=a(2)=2, a(n)=a(n-3) (so it is a periodic
> # sequence of the a(n) in the taylor series a(n)*x^n with period 3).

? GenFLinRec( [1,2,2], [0,0,1], 0 )
%2 = (-2*x^2 - 2*x - 1)/(x^3 - 1)

? %2 + O(x^10)
%3 = 1 + 2*x + 2*x^2 + x^3 + 2*x^4 + 2*x^5 + x^6 + 2*x^7 + 2*x^8 + x^9 + O(x^10)

> # Example 4:
> # GenFLinRec( [0,1], [3,k], 0, 0, true,[]) returns the generating function
> # for a(0)=0, a(1)=1, a(n)=3a(n-1)+k*a(n-2), see A083857 in the OEIS.

? GenFLinRec( [0,1], [3,k], 0 )
%4 = -x/(k*x^2 + 3*x - 1)

? subst(%4,k,2) + O(x^10)
%5 = x + 3*x^2 + 11*x^3 + 39*x^4 + 139*x^5 + 495*x^6 + 1763*x^7 +
6279*x^8 + 22363*x^9 + O(x^10)
(this corresponds to the last row in the example in A08385)

Regards,
Max