OEIS/Generating Functions and Recurrences: Difference between revisions

From tehowiki
Jump to navigation Jump to search
imported>Gfis
No edit summary
combstruct reference
 
(24 intermediate revisions by 2 users not shown)
Line 136: Line 136:
       = 2 - (numbers of generators in the name);
       = 2 - (numbers of generators in the name);
  ''trms'' = number of terms desired (with a default number of 20)
  ''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'').
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
[https://oeis.org/A162760 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
[https://oeis.org/A025225 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);
 
[https://oreis.org/A007750 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);
[https://oreis.org/A007788 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
[https://oeis.org/A162760 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===
* [https://oeis.org/A158454 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
* [https://oeis.org/A159764 A159764] Riordan array (1/(1+4x+x^2), x/(1+4x+x^2)).
G.f.: 1/(1+4*x+x^2-y*x).
* [https://oeis.org/A129818 A129818] Riordan array (1/(1+x), x/(1+x)^2)
O.g.f.: (1+x)/(1+(2-y)*x+x^2).
* [https://oeis.org/A109956 A109956] Inverse of Riordan array (1/(1-x), x/(1-x)^3), [https://oeis.org/A109955 A109955].
* MMA for exponential Riordan arrays in [https://oeis.org/A256893 A256893]
* '''more''': [https://oeis.org/A260780 A260780] with some diagonals
* [https://oeis.org/A207543 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===
* [https://oeis.org/A046816 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]); # [https://oeis.org/A299958 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]); # [https://oeis.org/A001263 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 ([https://oeis.org/A005325 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}]
 
===Cf. [[OEIS/Square Root Recurrences]]===
===Further reading===
* [https://algo.inria.fr/flajolet/Publications/book.pdf Analytic combinatorics] by Philippe Flajolet and Robert Sedgewick. Cambridge: Cambridge University Press, 2009, pp. xiv+810.

Latest revision as of 12:07, 20 October 2023

Linear 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)
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}]

Cf. OEIS/Square Root Recurrences

Further reading

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