Difference between revisions of "OEIS/Generating Functions and Recurrences"
m (→Holonomic Recurrences) |
(→Generating functions in your CAS) |
||
Line 147: | Line 147: | ||
seq(A(n,9), n=0..47); gfun[guessgf]([%],x)[1]; | 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); | 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)); | ||
read "C:\\Program Files\\Maple 2019\\FPS.mpl": | read "C:\\Program Files\\Maple 2019\\FPS.mpl": |
Revision as of 20:56, 6 May 2021
Contents
- 1 Linear Recurrences
- 2 Ordinary generating function o.g.f.
- 3 Linear Recurrence with constant coefficients
- 4 Initial terms
- 5 Offset
- 6 Computation of the o.g.f. from the linear recurrence
- 7 Coxeter Group Sequences
- 8 Generating functions in your CAS
- 9 Triangles, Riordan arrays, g.f. with two variables
- 10 Higher dimensions
- 11 Holonomic Recurrences
- 12 Cf. OEIS/Square Root Recurrences
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));
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);
- 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).
- more: A260780 with some diagonals
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}]