Integral of GeneratingFunction

Question

I know that GeneratingFunction can be used to compute the generating function $\sum_{n=0}^\infty a_n x^n$ of a sequence $(a_n)_n$ via

GeneratingFunction[a[n],n,x]

I also know that polynomial (in $n$) coefficients of the general term $a_n$ of the sequence result in an expression involving derivatives of the generating function, e.g.

GeneratingFunction[n a[n],n,x]
(* x (GeneratingFunction^(0,0,1))[a[n],n,x] *)

This does not seem to happen with rational (in $n$) coefficients, which would lead to indefinite integrals of the generating function. For example I would expect the input

GeneratingFunction[a[n]/(n + 1), n, x]

to result in

(*1/x Integrate[GeneratingFunction[a[n],n,y],{y,0,x}]*)

or

(*1/x GeneratingFunction^(0,0,-1))[a[n],n,x]*)

but this does not seem to be the case.

I thus have the following question: Is it possible to get Mathematica to compute the generating function of a linear combination of lagged values of $a_n$ with rational (in $n$) coefficients in terms of derivatives and integrals of the generating function of the sequence $(a_n)_n$?

Answer 1

In this answer, I handle coefficients of the form

$$ \frac{p(n)}{q(n)}, $$

where $p(n)$ and $q(n)$ are polynomials in the index symbol $n$, $\deg(p(n)) < \deg(q(n))$, and all of the roots of $q(n)$ are negative integers. If $\deg(p(n)) \geq \deg(q(n))$, we may have to differentiate a power series, and as I mentioned in my comment, Mathematica 9 doesn't handle that conveniently. If $q(n)$ has positive integer roots, we have to handle undefined expressions and things are more complicated.

First, some setup code

(* http://mathematica.stackexchange.com/questions/32531/how-to-unprotect-generatingfunction *)
GeneratingFunction;
protected = Unprotect[GeneratingFunction];

(* what can usually be used as a variable, according to 
    ref/message/General/ivar *)
variablePattern = Except[_String | _?NumberQ | _Plus | _Times | 
    _Sum | _Product | _^_Integer];

arg2 = Sequence[n:variablePattern, x:variablePattern];

Now, we add a downvalue for GeneratingFunction which expands $p(n)/q(n)$ into partial fractions.

GeneratingFunction[expr_, arg2, 
    opts:OptionsPattern[]] /; expr =!= Apart[expr, n] :=
    GeneratingFunction[Apart[expr, n], n, x, opts];

I will assume that the coefficients are now in the form

$$ \frac{c}{(n + k)^j}.$$

To go from here, in the next rule we add a downvalue saying

$$ \sum_{n \geq 0} \frac{1}{n+k} f(n) x^n = x^{-k} \int_0^x t^{k-1} \left( \sum_{n \geq 0} f(n) t^n \right) dt,$$

which gets used repeatedly until $j=0$.

GeneratingFunction[expr_ * (n_ + k_?Positive)^(j_?Negative), arg2, opts:OptionsPattern[]] := With[
    {t = Unique[]},
    x^(-k) * Integrate[t^(k-1) * GeneratingFunction[expr * (n + k)^(j + 1),
        n, t, opts], {t, 0, x}]
];

And we're done.

Protect[Evaluate[protected]];