OEIS/Generating Functions and Recurrences

Ordinary generating function o.g.f.

An o.g.f. consists of a fraction of two univariate polynomials (the numerator and the denominator polynomial). This rational function is expanded into a power series by polynomial division, beginning with the lowest power of the variable (we always use x), and the correspondant coefficients of the powers x are the terms of the integer sequences defined by the o.g.f.

Example for the polynomial division:

 1 / (1-x-x^2) = 1 + x + 2x^2 + 3x^3 + 5x^5 ... (Fibonnacci, A000045)
-1 + x + x^2    <---- denominator * (-1)
------------    <---- add
 0 + x +  x^2  
    -x +  x^2 + x^3
    --------------
     0 + 2x^2 + x^3
        -2x^2 + 2x^3 + 2x^4
       -------------------
          0   + 3x^2 + 2x^4
               -3x^3 + 3x^4 + 3x^5
               -------------------
                 0   + 5x^4 + 3x^5 ...

Corresponding Mathematica functions:

 CoefficientList[Series[numerator_polynomial / denominator_polynomial, {x, 0, number of terms}], x]
 CoefficientList[Series[-(x/(-1 + x + x^2)), {x, 0, 20}], x]
 FindGeneratingFunction[{terms}]

Linear Recurrence with constant coefficients

There is a one-to-one relation between an o.g.f. and a linear recurrence. The denominator's coefficients (without the constant 1) are the signature of the linear recurrence. The number of elements in the signature is the order of the linear recurrence.

In our example the linear recurrence is the following:

 a(0) = 1, a(1) = 1, a(n) = a(n - 1) + a(n - 2) for n > 1. 

The number of initial terms must be greater than or equal to the order.

Corresponding Mathematica functions:

 LinearRecurrence[{signature}, {initial terms}, number of terms]
 LinearRecurrence[{1, 1}, {1, 1}, 8] for the Fibonnacci sequence
 FindLinearRecurrence[{terms}]

Linear recurrences can be implemented very easily by additions and multiplications only. An array of length order may be used as a rolling buffer where all indexes are taken mod order.

Initial terms

OEIS sequences often start out with some initial terms before the recurrence really starts. For example:

A049221 Number of horizontally convex n-ominoes in which the top row has exactly 1 square, 
        which is not above the rightmost square in the second row.
Offset 1
1, 0, 1, 4, 14, 46, 148, 474, 1518, 4864, 15590, 49974, 160196
a(n) = 5a(n-1) - 7a(n-2) + 4a(n-3) for n >= 6
Linear recurrence with constant coefficients, signature (5,-7,4)

The condition n>=6 and the order 3 show that the recurrence only starts at a(3) = 1, and that a(1) = 1 and a(2) = 0 were prefixed to the sequence. The second parameter of the Mathematica call

LinearRecurrence[{5, -7, 4}, {1, 0, 1, 4, 14}, 40]

is convenient because such initial terms may simply be written before the terms belonging to the recurrence.

How do we see such initial terms in the o.g.f.? The numerator and the denominator polynomials should have the same degree. If the numerator has a higher degree, we can transform the fraction into mixed form:

G.f.: x*(1-x)^2*(1-3x+x^2)/(1-5x+7x^2-4x^3) 
  = (x-5x^2+6x^3-5x^4+x^5)/(1-5x+7x^2-4x^3)
  = 0  + 1x + 0x^2 + x^3(1-x+x^2)/(1-5x+7x^2-4x^3)
    a(0) a(1) a(2)   a(3 ...)

Now we can directly recognize the additional initial terms.

Offset

The power series of generated by an o.g.f. usually has the general form

c0*x^0 + c1*x^1 + c2*x^2 + c3*x^3 ...

The integer sequence defined by this o.g.f. then simply is

c0, c1, c2, c3, ...

Offset 0

All OEIS sequences are assigned an offset that relates them to the context for which they are defined. The most natural offset for a sequence defined by an o.g.f. is 0, that is the sequence starts with a(0) = c0.

Positive Offsets

But many OEIS sequences have offset 1 or higher? The convention in the OEIS is that for an offset i > 0 the expansion should still start with c0 + c1*x + c2*x^2 + ..., but that all constants cj with j < i do not belong to the sequence and are ignored.

As a rule of thumb, the offset can be shifted from 0 to i by multiplying the numerator polynomial with x^i.

Negative Offsets

There are even OEIS sequences with negative offsets defined by g.f.s, for example:

A066280 a(n) = 1^n + 2^(n+1) + 3^(n+2).
5, 12, 32, 90, 260, 762, 2252, 6690, 19940, ...
Offset -1
Linear recurrence with constant coefficients, signature (6, -11, 6).
a(n) = 6*a(n-1) - 11*a(n-2) + 6*a(n-3).
G.f.: (-18*x+15*x^2+5)/((1-x) * (3*x-1) * (2*x-1) * x)
 = (5-18*x+15*x^2)/(x-6x^2+11x^3-6x^4)

Here the denominator polynomial has no constant, therefore the division first yields a factor 5/x = 5*x^(-1).

Computation of the o.g.f. from the linear recurrence

Example (from youtube.com):

a(n) = 2*a(n-1) + 4*a(n-2); a(0) = 1, a(1) = 3; [1, 3, 10, 32, ...]
a(n) - 2*a(n-1) - 4*a(n-2) = 0

 A       = 1 + 3x + 10x^2 + 32x^3 + ...
-(2x*A   =     2x +  6x^2 + 20x^3 + ...)
-(4x^2*A =           4x^2 + 12x^3 + ...)
----------------------------------------
(1 - 2x - 4x^2)*A = 1+x

A = (1 + x) / (1 - 2x - 4x^2) 

Example (from youtube.com):

a(n) = a(n-1) + 2*a(n-2) + 3; a(0) = 2, a(1) = 2; [2, 2, 9, 16, 37, 72, 149, 296, ...]
a(n) - a(n-1) - 2*a(n-2) = 3

 A       = 2 + 2x + 9x^2 + 16x^3 + 37x^4 + ...
-(x*A    =     2x + 2x^2 +  9x^3 + 16x^4 + ...)
-(2x^2*A =          4x^2 +  4x^3 + 18x^4 + ...)
----------------------------------------
(1 - x - 2x^2)*A = 2 + 3x^2 +3x^3 + 3x^4 + ... = 2 + 3x^2/(1-x) = (2 - 2x + 3x^2)/(1-x)

A = (2 - 2x + 3x^2) / ((1 - x - 2x^2)*(1 - x))

Example:

A173652: 0, 0, 0,  1,   4,   8,  33,  68, 245, 549, 1815 ...
                                            x^3  x^4  x^5  x^6  x^7  x^8  x^9  x^10
LinearRecurrence[{0,    9,     3,    -17,    -8,     6,     7,    -1,    1}, {0, 0, 0, 1, 4, 8, 33, 68, 245}, 31]
denominator: (1 + 0x^1 - 9x^2 - 3x^3 + 17x^4 + 8x^5 - 6x^6 - 7x^7 + x^8 - x^9)
numerator: x^3

Coxeter Group Sequences

There is a group of 2302 similiar sequences for Coxeter groups. The main entry is:

A169452 
Number of reduced words of length n in Coxeter group on 7 generators S_i 
with relations (S_i)^2 = (S_i S_j)^33 = I.
1, 7, 42, 252, 1512, 9072, 54432, 326592, 1959552
signature(5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,-15).
G.f.:(t^33+2*t^32+2*t^31+2*t^30
+2*t^29+2*t^28+2*t^27+2*t^26+2*t^25+2*t^24+2*t^23+2*t^22+2*t^21+2*t^20
+2*t^19+2*t^18+2*t^17+2*t^16+2*t^15+2*t^14+2*t^13+2*t^12+2*t^11+2*t^10
+2*t^9+2*t^8+2*t^7+2*t^6+2*t^5+2*t^4+2*t^3+2*t^2+2*t+1) / 
(15*t^33-5*t^32-5*t^31-5*t^30
-5*t^29-5*t^28-5*t^27-5*t^26-5*t^25-5*t^24-5*t^23-5*t^22-5*t^21-5*t^20
-5*t^19-5*t^18-5*t^17-5*t^16-5*t^15-5*t^14-5*t^13-5*t^12-5*t^11-5*t^10
-5*t^9-5*t^8-5*t^7-5*t^6-5*t^5-5*t^4-5*t^3-5*t^2-5*t+1)

• Harvey P. Dale gave this Mathematica program for the g.f. on Aug 16 2014:

coxG[{pwr_, c1_, c2_, trms_:20}] := Module[{
num=Total[  2t^Range[pwr-1]]+   t^pwr+1, 
den=Total[c2*t^Range[pwr-1]]+c1*t^pwr+1}, 
CoefficientList[Series[num/den, {t, 0, trms}], t]]; 
coxG[{33, 15, -5, 30}] (* for A169452 *)

• with the parameters:

pwr  = largest exponent in the o.g.f. and of "(S_i S_j)" in the name,
c1   = first coefficient in the denominator of the g.f., 
     = triangular(-c2) = binomial(-c2 + 1, 2) 
c2   = negative of second coefficient in the denominator of the o.g.f., 
     = 2 - (numbers of generators in the name);
trms = number of terms desired (with a default number of 20)

So all these sequences have a o.g.f which looks complicated, but which can be easily derived by from only two parameters (pwr, c2).

Generating functions in your CAS

• PARI

A162760 (PARI) my(x='x+O('x^20)); Vec((1+x)*(1-x^3)/(1-10*x+54*x^3-45*x^4)) \\ _G. C. Greubel_, Apr 26 2019
A025225: a(n)=polcoeff((1-sqrt(1-8*x+x*O(x^n)))/2, n)
for(n=0,32,print(n" ",polcoeff(x*(1-2*x^2*(1+x-x^3-x^4-x^5)+O(x^n))/(1+x),n)))

• Mathematica

CoefficientList[Series[polx1/polx2, {x, 0, 29}], x]

• Maple

seq(A(n,9), n=0..47); gfun[guessgf]([%],x)[1];
seq(coeff(series( (1+x*exp(x))/(1-x), x, n+1)*n!, x, n), n = 0..30);
with(StringTools):
DeleteSpace(convert([seq(coeff(expand((6*n-5)*(6*n-4)*(6*n-3)*(6*n-2)*(6*n-1)*n),n,j),j=0..6)],string));

ogf := (x^4 + 2*x^3 + 7*x^2 - 6*x + 1)^(-1/2):
series(ogf, x, 27): seq(coeff(%, x, n), n=0..25); # A299443

S:= series((1 + sqrt(1-2*sin(2*x)))/(2*cos(x)), x, 51):
seq(coeff(S, x, n)*n!, n=0..16); # A318599, exponential g.f.

# A128096
g:=4*(1-z^2-sqrt((1+z+z^2)*(1-3*z+z^2)))/(1-z+z^2+sqrt((1+z+z^2)*(1-3*z+z^2)))^2: 
gser:=series(g, z=0, 38): seq(coeff(gser, z, n), n=1..35);

read "C:\\Program Files\\Maple 2019\\FPS.mpl":
interface(prettyprint=0):
with(gfun):
timelimit(4,diffeqtorec(HolonomicDE(1/(1+x-sqrt(x)),F(x)),F(x),a(k)));
(2*k^2+5*k+3)*a(k)+(-6*k-12)*a(k+1)+(-6*k-12)*a(k+2)+(-2*k^2-11*k-15)*a(k+3)

a := n -> 3^n*hypergeom([-n/3, (1-n)/3, (2-n)/3], [1, 1], 1):
seq(simplify(a(n)), n=0..24);

A007750 m:=30; s:=series(x*(1+6*x+x^2)/((1-x)*(1-16*x^2+x^4)), x, m+1): seq(coeff(s, x, j), j=0..m); 
A007788 m:=20; s:=series( (1-sin(3*x))^(-1/3), x, m+1): seq(j!*coeff(s, x, j), j=0..m); # G. C. Greubel

• Sage

A162760 (Sage) ((1+x)*(1-x^3)/(1-10*x+54*x^3-45*x^4)).series(x, 20).coefficients(x, sparse=False)

Triangles, Riordan arrays, g.f. with two variables

• A158454 Riordan array (1/(1-x^2), x/(1+x)^2).

1, 0, 1, 1, -2, 1, 0, 4, -4, 1, 1, -6, 11, -6, 1, 0, 9, -24, 22
O.g.f. column k with leading zeros: (1/(1-x^2))*(x/(1+x)^2)^k, k >= 0.
O.g.f. row polynomials (rising powers in y): ((1+x)/(1-x))/(1+(2-y)*x+x^2)
With[{m = 15}, CoefficientList[CoefficientList[Series[(1+x)/((1-x)*(1 + x)^2 - y*x*(1-x)), {x, 0, m}, {y, 0, m}], x], y]] // Flatten

• A159764 Riordan array (1/(1+4x+x^2), x/(1+4x+x^2)).

G.f.: 1/(1+4*x+x^2-y*x).

• A129818 Riordan array (1/(1+x), x/(1+x)^2)

O.g.f.: (1+x)/(1+(2-y)*x+x^2).

• A109956 Inverse of Riordan array (1/(1-x), x/(1-x)^3), A109955.
• MMA for exponential Riordan arrays in A256893
• more: A260780 with some diagonals
• A207543 Triangle read by rows, expansion of (1+y*x)/(1-2*y*x+y*(y-1)*x^2).

s := (1+y*x)/(1-2*y*x+y*(y-1)*x^2): t := series(s, x, 12):
seq(print(seq(coeff(coeff(t, x, n), y, m), m=0..n)), n=0..11);

Higher dimensions

• A046816 Pascal's tetrahedron: entries in 3-dimensional version of Pascal triangle, or trinomial coefficients.

G.f.: 1/(1-x-x*y-x*y*z) 
Mathematica: Flatten[CoefficientList[CoefficientList[CoefficientList[Series[1/(1-x-x*y-x*y*z), {x, 0, 6}], x], y], z]]

Holonomic Recurrences

Coefficients of a[n-k] are polynomials in n.

• Maple

f:= gfun:-rectoproc({(25000*(2*n+7))*(4*n-1)*(4*n+9)*(n+1)*a(n)+(n+5)*(n+4)*(n+3)*(n+2)*a(n+5), 
a(0)=1, a(1)=-5, a(2)=-25, a(3)=-125, a(4)=0}, a(n), remember): map(f, [$0..40]); # A299958
n:=8; f:= gfun:-rectoproc({k*(k+1)*a(k) - (k-n)*(k-n-1)*a(k-1) = 0, a(0)=1}, a(k), remember): map(f, [$0..n-1]); # A001263

restart; # A094914
gf := sqrt((x^2 + x - 1354363)^2);
diffeqtorec(HolonomicDE(gf,F(x)),F(x),a(k));
(-2+k)*a(k)+k*a(k+1)+(-1354363*k-2708726)*a(k+2)

• Mathematica (A005325)

RecurrenceTable[{3(-1+n)*n*a[-2+n]+n*(1+2n)*a[-1+n]-(-5+n)*(7+n)*a[n]==0, a[5]==1,a[6]==6}, a,{n,5,20}]

Further reading

• Analytic combinatorics by Philippe Flajolet and Robert Sedgewick. Cambridge: Cambridge University Press, 2009, pp. xiv+810.