User:Peter Luschny/GeneratingFunctions

G F U N
A Sage Toolbox for the OEIS
Part 1: Basic operations

Introduction

A generating function is a formal power series $\scriptstyle a(x)=\sum _{n=0}^{\infty }a_{n}x^{n}$ which lets us calculate with the sequence $\scriptstyle (a_{0},a_{1},a_{2},\ldots )$ as if it were an analytic function. Compared with the ring of sequences the ring of power series has two great advantages: It is an integral domain (which the ring of sequences is not) and its multiplication, the Cauchy product, has a natural combinatorial interpretation.

Gfun is a classic Maple package for the manipulation of generating functions by Bruno Salvy and Paul Zimmermann written in 1994 swmath. We borrow the name gfun and intend to provide in a series of posts on this blog a similar toolbox, implemented in Sage and geared for the use with the OEIS. In this post we start with the most basic operations and examples.

Our implementation is object oriented and is hopefully easier and more intuitive. We illustrate this with an example, the computation of A048172.

The Maple version:

with(gfun): 
f := series((ln(1+x)-x^2/(1+x)), x, 21):
egf := seriestoseries(f, 'revogf'): 
seriestolist(egf, 'Laplace');

The same with our Sage toolbox:

gf = log(1+x)-x^2/(1+x)
gfun(gf, 'rev').as_egf()

The constructor gfun

The constructor gfun (g, op) has two arguments, g a symbolic expression in one or two variables and op the name of an operation according to the table below. It returns an object representing a formal power series which can be accessed through the methods listed below.

The operations op:

• None    ... g unmodified (the default case with no parameter).
• 'rev'   ... compositional inverse of g.
• 'inv'   ... multiplicative inverse of g.
• 'exp'   ... exponential of g.
• 'log'   ... logarithm of g.
• 'int'   ... formal integral of g.
• 'der'   ... formal derivative of g.
• 'logder' .. logarithmic derivative of g.

The methods

The methods return the sequence of coefficients of the represented power series up to the given precision prec (which is by default 12). The type of the power series is indicated by the method name. The types are 'ogf' (ordinary power series), 'egf' (exponential power series), 'lgf' (logarithmic power series) and 'bgf' (binary power series).

• as_ogf    (prec=12, search=false)
• as_egf    (prec=12, search=false)
• as_lgf    (prec=12, search=false)
• as_bgf    (prec=12, search=false)
• as_series (prec=12, typ='ogf')

Setting the parameter 'search=true' (default is false) will look up related information in the OEIS provided you are connected with the Internet.

In the bivariate case (i. e. the given symbolic expression depends on two variables x and t) 'coeffs' returns a triangular table giving the coefficients of those polynomials which are the coefficients of the power series when expanded in t. 'values' returns a rectangular array displaying the values of these polynomials evaluated at the nonnegative integers. The values method evaluates the n-th polynomial in the n-th column of the rectangular array. The 'colgenerator' returns a rational generating function of the n-th column of the values array.

• coeffs (typ='ogf', rows = 8)
• values (typ='ogf', rows = 8, cols = 8)
• colgenerator (n, typ='ogf')
• associates (typ='ogf', prec=8)

The method 'associates' is nonstandard. It lists the numerators of the partial fraction decomposition of the generating functions of the columns of the value array as an triangular array. Many examples show that this leads to interesting and meaningful integer triangles. Apart from a few random examples they seem to be almost completely absent in the OEIS. I think they are worth to be included in the OEIS in order to promote their further study.

Example use

To load the toolbox we have to execute the next command first. The path must of course be adapted to the computing environment.

attach('~/.sage/gfun.sage')

The univariate case.

x = SR.var('x')
cato = (1-2*x-sqrt(1-4*x))/(2*x)

Note: The type of input gfun expects is a symbolic expression. ('SR' in 'x = SR.var('x')' stands for Symbolic Ring.)

gf = gfun((1-2*x-sqrt(1-4*x))/(2*x))
print gf.as_ogf()
print gf.as_ogf(6)

[0, 1, 2, 5, 14, 42, 132, 429, 1430, 4862, 16796, 58786] [0, 1, 2, 5, 14, 42]

If you want more (or fewer) terms of the sequence call as_ogf with this value. You can also call as_ogf without a value. Then the default length (= 12) is used.

You don't have to re-execute gfun to get a different length of the sequence. The returned object behaves as if it can provide unlimited access to the sequence. This is in contrast to the way some similar (Maple, Mathematica) libraries work which can be found in the OEIS. These programs work on lists; there you have to assign a fixed length at the start which you can't change afterwards. More generally the objects in these libraries behave like 1-based lists whereas the objects here behave like 0-based sequences.

You also can execute different 'as_xgf' methods on the same object. In combination with the operations this gives a multitude of different representations of sequences.

R = gfun(cato, 'rev')
print R.as_ogf()     # A181983
print R.as_egf(9)    # A001563
print R.as_lgf()     # A000290
print R.as_bgf(10)   # A036289

[0, 1, -2, 3, -4, 5, -6, 7, -8, 9, -10, 11] [0, 1, -4, 18, -96, 600, -4320, 35280, -322560] [0, 1, 4, 9, 16, 25, 36, 49, 64, 81, 100, 121] [0, 2, -8, 24, -64, 160, -384, 896, -2048, 4608]

S = gfun(1 - cato, 'inv') 
print S.as_ogf(9) # A088218 
print S.as_egf(9) # A000407 
print S.as_lgf(9) # A002457 
print S.as_bgf(9) # A069723

[1, 1, 3, 10, 35, 126, 462, 1716, 6435] [1, 1, 6, 60, 840, 15120, 332640, 8648640, 259459200] [0, 1, -6, 30, -140, 630, -2772, 12012, -51480] [1, 2, 12, 80, 560, 4032, 29568, 219648, 1647360]

The bivariate case.

The Motzkin polynomials, A097610.

x, t = SR.var('x, t')
motzkin = gfun((1-t*x-sqrt((x^2-4)*t^2-2*t*x+1))/(2*t^2))
motzkin.as_ogf(5)

[1, x, x^2 + 1, x^3 + 3*x, x^4 + 6*x^2 + 2]

motzkin.coeffs()

[[1], [0, 1], [1, 0, 1], [0, 3, 0, 1], [2, 0, 6, 0, 1], [0, 10, 0, 10, 0, 1], [5, 0, 30, 0, 15, 0, 1], [0, 35, 0, 70, 0, 21, 0, 1]]

Motzkin polynomials evaluated at the integers, A247495.

motzkin.values()

[n\k][0][1] [2] [3] [4] [5] [6] [7] [0] 1, 0, 1, 0, 2, 0, 5, 0,... A126120 [1] 1, 1, 2, 4, 9, 21, 51, 127,... A001006 [2] 1, 2, 5, 14, 42, 132, 429, 1430,... A000108 [3] 1, 3, 10, 36, 137, 543, 2219, 9285,... A002212 [4] 1, 4, 17, 76, 354, 1704, 8421, 42508,... A005572 [5] 1, 5, 26, 140, 777, 4425, 25755, 152675,... A182401 [6] 1, 6, 37, 234, 1514, 9996, 67181, 458562,... A025230 A000012,A001477,A002522,A079908, ...

Using the rational generating functions of the columns gives the array with rows and columns exchanged.

for n in range(8):
    gf = motzkin.colgenerators(n)
    print gfun(gf).as_ogf(7)

[1, 1, 1, 1, 1, 1, 1] [0, 1, 2, 3, 4, 5, 6] [1, 2, 5, 10, 17, 26, 37] [0, 4, 14, 36, 76, 140, 234] [2, 9, 42, 137, 354, 777, 1514] [0, 21, 132, 543, 1704, 4425, 9996] [5, 51, 429, 2219, 8421, 25755, 67181] [0, 127, 1430, 9285, 42508, 152675, 458562]

The column generating functions are given in the form of partial fractions.

for n in range(6):
    print motzkin.colgenerators(n)

-1/(x-1) 1/(x-1) + 1/(x-1)^2 -2/(x-1) - 3/(x-1)^2 - 2/(x-1)^3 4/(x-1) + 10/(x-1)^2 + 12/(x-1)^3 + 6/(x-1)^4 -9/(x-1) - 33/(x-1)^2 - 62/(x-1)^3 - 60/(x-1)^4 - 24/(x-1)^5 21/(x-1)+111/(x-1)^2+300/(x-1)^3+450/(x-1)^4+360/(x-1)^5+120/(x-1)^6

Collecting the numerators of these partial fractions gives (apart from the signs) the 'associated' number triangle A247497.

motzkin.associates()

[[ -1], [ 1, 1], [ -2, -3, -2], [ 4, 10, 12, 6], [ -9, -33, -62, -60, -24], [ 21, 111, 300, 450, 360, 120], [-51, -378, -1412, -3000, -3720, -2520, -720], [127, 1303, 6552, 18816, 32760, 34440, 20160, 5040]]

Looking at the absolute values of these numbers the first column lists the Motzkin numbers and the second column A058987, Catalan(n) - Motzkin(n-1). The last column shows (again apart from signs) the factorial numbers and the next to last column is A001710, the order of the alternating group.

Summary: We started with the Motzkin polynomials A097610, looked at the Motzkin values A247495 and arrived via the partial fractions of the column generating functions at the associated Motzkin numbers A247497 which include the Motzkin numbers A001006 as a special case. It is this pattern which we will follow in the examples below.

To give a second example we look at Pascal's triangle and its exponential generating function.

x, t = SR.var('x, t')
binomial = gfun(exp(t+x*t))
binomial.as_egf(5)

[1, x+1, x^2+2*x+1, x^3+3*x^2+3*x+1, x^4+4*x^3+6*x^2+4*x+1]

binomial.coeffs('egf')

[[1], [1, 1], [1, 2, 1], [1, 3, 3, 1], [1, 4, 6, 4, 1], [1, 5, 10, 10, 5, 1], [1, 6, 15, 20, 15, 6, 1], [1, 7, 21, 35, 35, 21, 7, 1]]

binomial.values('egf')

[[1, 1, 1, 1, 1, 1, 1, 1], [1, 2, 4, 8, 16, 32, 64, 128], [1, 3, 9, 27, 81, 243, 729, 2187], [1, 4, 16, 64, 256, 1024, 4096, 16384], [1, 5, 25, 125, 625, 3125, 15625, 78125], [1, 6, 36, 216, 1296, 7776, 46656, 279936], [1, 7, 49, 343, 2401, 16807, 117649, 823543], [1, 8, 64, 512, 4096, 32768, 262144, 2097152]]

binomial.associates('egf')

[[-1], [ 0, 1], [ 0, -1, -2], [ 0, 1, 6, 6], [ 0, -1, -14, -36, -24], [ 0, 1, 30, 150, 240, 120], [ 0, -1, -62, -540, -1560, -1800, -720], [ 0, 1, 126, 1806, 8400, 16800, 15120, 5040]]

The associated binomial numbers are in their signless form easily identified as k!*Stirling2(n, k) = A000142(k) *A048993(n, k) = A131689(n, k).

It appears appropriate now to consider how the relationships between the last three tables are represented in mathematical notation.

$\sum _{j=0}^{n}{\binom {n}{j}}k^{j}=(k+1)^{n}=[x^{k}]\sum _{j\geq 0}{\frac {(-1)^{n+1}j!}{(x-1)^{j+1}}}\left\{{n \atop j}\right\}$

More generally let us denote the coefficients by C(n, k), the values by V(n, k) and the associates by A(n, k). Then we can write:

$\sum _{j=0}^{n}C(n,j)x^{j}{\biggr |}_{x=k}=V(n,k)=[x^{k}]\sum _{j\geq 0}{\frac {A(n,j)}{(x-1)^{j+1}}}$

 

---

Next let us look at the Stirling numbers of the first kind:

x, t = SR.var('x, t')
stirling1 = gfun(exp(x*log(1+t)))
stirling1.as_egf(5)

[1, x, x^2 - x, x^3 - 3*x^2 + 2*x, x^4 - 6*x^3 + 11*x^2 - 6*x]

stirling1.coeffs('egf')

[[1], [0, 1], [0, -1, 1], [0, 2, -3, 1], [0, -6, 11, -6, 1], [0, 24, -50, 35, -10, 1], [0, -120, 274, -225, 85, -15, 1], [0, 720, -1764, 1624, -735, 175, -21, 1]]

stirling1.values('egf')

[[1, 0, 0, 0, 0, 0, 0, 0], [1, 1, 0, 0, 0, 0, 0, 0], [1, 2, 2, 0, 0, 0, 0, 0], [1, 3, 6, 6, 0, 0, 0, 0], [1, 4, 12, 24, 24, 0, 0, 0], [1, 5, 20, 60, 120, 120, 0, 0], [1, 6, 30, 120, 360, 720, 720, 0], [1, 7, 42, 210, 840, 2520, 5040, 5040]]

stirling1.associates('egf')

[[ -1], [ 1, 1], [ -2, -4, -2], [ 6, 18, 18, 6], [ -24, -96, -144, -96, -24], [ 120, 600, 1200, 1200, 600, 120], [-720, -4320, -10800, -14400, -10800, -4320, -720], [5040, 35280, 105840, 176400, 176400, 105840, 35280, 5040]]

The coefficients are the Stirling numbers of the first kind A048994 (which we prefer over Knuth's choice A132393 and strongly prefer over A008275 or the variant A130534). The nonzero values are the permutation coefficients A008279 and the associate numbers are apart from signs the coefficients in the expansion of x^n in terms of Laguerre polynomials A021012.

So we have C(n, k) = Stirling1(n, k), V(n, k) = n!/(n-k)! if k ≤ n else 0 and A(n, k) = (-1)^(n+1)*n!*binomial(n,k). If we plug this into our formula template given above we arrive at:

$\sum _{j=0}^{n}\left[{n \atop j}\right]k^{j}={\frac {n!}{(n-k)!}}[k\leq n]=[x^{k}]\sum _{j\geq 0}{\frac {(-1)^{n+1}j!}{(x-1)^{j+1}}}{\binom {n}{j}}$

Comparing this identity with the identity in the Pascal case we observe that the binomial coefficient has moved from the left hand side to the right hand side taking the place of the Stirling2 numbers which moved from the right hand side to the left hand side thereby transforming into the Stirling1 numbers.

The comparison now arouses our curiosity: what sum on the left hand side will we get if we plug the Stirling1 numbers in the right hand side in place of the binomial coefficients respectively the Stirling2 numbers?

$[x^{k}]\sum _{j\geq 0}{\frac {(-1)^{n+1}j!}{(x-1)^{j+1}}}\left[{n \atop j}\right]=\sum _{j=0}^{n}\quad ?\quad k^{j}$

The problem is that now we have to reverse the direction of the computation. What we know is A(n, k) = (-1)^(n+1)*A225479(n, k) and (by looking at the coefficients of the expansions) the values V(n, k):

[[1, 1, 3, 14, 88, 694, 6578, 72792], [1, 2, 8, 46, 342, 3108, 33324, 411360], [1, 3, 15, 102, 870, 8892, 105708, 1431168], [1, 4, 24, 188, 1804, 20416, 265640, 3901320], [1, 5, 35, 310, 3300, 40890, 576870, 9116760], [1, 6, 48, 474, 5538, 74484, 1131108, 19115832], [1, 7, 63, 686, 8722, 126448, 2054864, 36948240], [1, 8, 80, 952, 13080, 203232, 3517008, 66998448]]

Since we know that the columns are polynomial sequences it is easy to find the polynomials by well known algorithms, for example by the Newton polynomial interpolation.

1, 1 + 1*k, 3 + 4*k + 1*k^2, 14 + 22*k + 9*k^2 + 1*k^3, 88 + 155*k + 82*k^2 + 16*k^3 + 1*k^4, 694 + 1333*k + 835*k^2 + 220*k^3 + 25*k^4 + 1*k^5, 6578 + 13541*k + 9598*k^2 + 3085*k^3 + 485*k^4 + 36*k^5 + 1*k^6.

Extracting the coefficients gives us C(n, k). So it's time to use a well known online tool for identifying sequences and search for instance for [694 1333 835 220 25 1]. But the oracle answered: "Sorry, but the terms do not match anything in the table." This was much to my surprise; I thought somewhat naively that a classic sequence would show up.

Meanwhile I have entered the sequence into the database. It is now A253165. And using the methods of 'gfun' you can easily explore this new number triangle starting from the bivariate exponential generating function $\scriptstyle (1+\log(1-t))^{-x-1}$.

x, t = SR.var('x, t')
A253165 = gfun(1/(1+log(1-t))^(x+1))

The 'explore' methods

• explore_values(g, typ='ogf', rows=6, cols=6)

This methods is a convenience function. Given a bivariate generating function g it searches the OEIS for information about sequences in the rows and columns of the value table (you have to be online). It helps to quickly get an overview of the situation and to build a cross-reference table like the one given above.

For instance the call

gfun.explore_values(motzkin)        

tries to find information about the first 6 rows and the first 6 columns of the array of Motzkin values.

• explore_genfun(g)

Given an univariate generating function g this method searches for all combinations of the above operations with the above methods which make mathematical sense and give an integer sequence. If such sequences exist it will search the OEIS for information about these sequences.

The information provided by explore_genfun() is often very useful while studying a sequence, sometimes even more useful and inspiring than the information returned by the SuperSeeker. This is because it does not only hint to information that can be found in the OEIS but also indicates what possible valuable information is missing. Note that you have to be connected to the internet to use this method and it may take some time to finish.

After the method has finished it displays a plot which is a kind of a fingerprint of the generating function. A green path indicates that there are related sequences in the OEIS (which are listed in the output), a red path indicates that in this case a meaningful integer sequence exists which is apparently not in the database. In the other cases for one reason or another there is no related integer sequence.

This plot is a snapshot of the current situation for the Riordan numbers. It was generated by the commands:

x = SR.var('x')
riordan = (1+x-sqrt(1-2*x-3*x^2))/(2*x*(1+x))
gfun.explore_genfun(riordan)

Conventions

All sequences here are by definition functions from N = {0,1,2,3,...} to some codomain. Triangles and rectangles are functions from N X N to some codomain. In other words all sequences have 'offset' 0 and all two-dimensional arrays are (0, 0) based. Lists in the OEIS might deviate from this: often they have (1, 1) or (1,0) or (0,1) as base-point of enumeration. This means that references to the sequences in the OEIS are always modulo the different conventions.

The software

An extended version of this post with more details and many examples can be found in this Notebook on Gfun. The file gfun.sage can be downloaded from github. Even better: clone it and extend it!