Combinatorics with SageMath: Generating functions¶

Nadia Lafrenière¶

Dartmouth College (🇺🇸) $\xrightarrow{🚌}$ Concordia University (🇨🇦)¶

EAUMP-ICTP School on Enumerative Combinatorics¶

Arusha, Tanzania, July 26, 2023¶

Acknowledgements:

  • Duncan Levear, from Boston College, made a part of this worksheet (in 2019), and gave me permission to use it
  • Steven Melczer's worksheet on the Fibonacci numbers

A motivating example¶

Professor Balázs gave you this exercise:

a. In how many ways can I pay 10 shillings if I am only allowed 1 and 2-shilling coins? How about 11 shillings? 20 shillings?

b. In how many ways can I pay $n$ shillings, if I am only allowed 1 and 2-shilling coins? Can you find a general formula for the number of ways? Can you prove your formula?

I would like you to say how you would do it with combinatorial structures, as we saw them yesterday.

We want to partition the integer 10 in any way with only parts 1 and 2!

In [1]:
Partitions(10, max_part=2).cardinality()
Out[1]:
6

We could also check how many that is for all the numbers between 10 and 20.

In [2]:
for n in range(10, 21):
    print(n, Partitions(n, max_part=2).cardinality())
10 6
11 6
12 7
13 7
14 8
15 8
16 9
17 9
18 10
19 10
20 11

It seems that the sequence is [1,2,2,3,3,4,4,,..].

If you did that exercise the way Professor Balász suggested you, you would get the generating function $$f(x) = \frac{-x^2 + x + 1}{x^3 - x^2 - x + 1}.$$

Can we expand it as to get the answer for 10, 11 and 20 shillings all at once?

In [3]:
f(x) = (-x^2 + x + 1)/(x^3 - x^2 - x + 1)
f.series(x,21)
Out[3]:
x |--> 1 + 2*x + 2*x^2 + 3*x^3 + 3*x^4 + 4*x^5 + 4*x^6 + 5*x^7 + 5*x^8 + 6*x^9 + 6*x^10 + 7*x^11 + 7*x^12 + 8*x^13 + 8*x^14 + 9*x^15 + 9*x^16 + 10*x^17 + 10*x^18 + 11*x^19 + 11*x^20 + Order(x^21)

Generating functions might seem like a complicated technique for this, in comparison with the partition approach. The power of generating function comes from the fact that one could get much more information about them. For example, we could simultaneously keep track of the number of coins needed for each way, and the amount. Then, the coefficients of the bivariate generating functions are the number of ways to pay a fixed amount using a fixed number of coins.

Bivariate generating functions¶

[Lecture at the board: bivariate generating functions]¶

Motivating example: the chicken nuggets problem¶

McDonald's, a famous chain of fast food in America, used to sell chicken nuggets in packs of 6, 9 or 20. What is the largest number of nuggets that one cannot order?

In [1]:
S = (1/((1-x^6)*(1-x^9)*(1-x^20))).series(x,100); S
Out[1]:
1 + 1*x^6 + 1*x^9 + 1*x^12 + 1*x^15 + 2*x^18 + 1*x^20 + 1*x^21 + 2*x^24 + 1*x^26 + 2*x^27 + 1*x^29 + 2*x^30 + 1*x^32 + 2*x^33 + 1*x^35 + 3*x^36 + 2*x^38 + 2*x^39 + 1*x^40 + 1*x^41 + 3*x^42 + 2*x^44 + 3*x^45 + 1*x^46 + 2*x^47 + 3*x^48 + 1*x^49 + 2*x^50 + 3*x^51 + 1*x^52 + 2*x^53 + 4*x^54 + 1*x^55 + 3*x^56 + 3*x^57 + 2*x^58 + 2*x^59 + 5*x^60 + 1*x^61 + 3*x^62 + 4*x^63 + 2*x^64 + 3*x^65 + 5*x^66 + 2*x^67 + 3*x^68 + 5*x^69 + 2*x^70 + 3*x^71 + 6*x^72 + 2*x^73 + 4*x^74 + 5*x^75 + 3*x^76 + 3*x^77 + 7*x^78 + 2*x^79 + 5*x^80 + 6*x^81 + 3*x^82 + 4*x^83 + 7*x^84 + 3*x^85 + 5*x^86 + 7*x^87 + 3*x^88 + 5*x^89 + 8*x^90 + 3*x^91 + 6*x^92 + 7*x^93 + 4*x^94 + 5*x^95 + 9*x^96 + 3*x^97 + 7*x^98 + 8*x^99 + Order(x^100)

I would like to know what powers do not appear in this expansion, without having to inspect it manually.

In [2]:
nonuggets = [t[1] for t in S.coefficients() if t[0]==0]
print("You cannot make packs of " + str(nonuggets))
You cannot make packs of [1, 2, 3, 4, 5, 7, 8, 10, 11, 13, 14, 16, 17, 19, 22, 23, 25, 28, 31, 34, 37, 43]

Two more combinatorial questions¶

How many ways are there to make 60 chicken nuggets?

How many of those ways use $7$ packs?

In [5]:
# How many ways are there to make 60 chicken nuggets?
# In class

There are 5 ways to get 60 chicken nuggets.

To get the number of ways that use 7 packs, I need a variable that keep tracks of the number of packs. Let $y$ be this variable.

Then, the generating function needs to be something like $$\sum_{n\geq 0}\sum_{k\geq 0}\sum_{P\text{ partition of $n$ nuggets into $k$ packs of 6, 9 or 20}} y^k x^n.$$

This is the same as $$ \frac{1}{(1-yx^6) (1-yx^9)(1-yx^{20})}.$$

In [6]:
# In class

So there is a single way to make 7 packs of nuggets, for a total of 60.

We can also add variables that keep track of all the types of packs (either six ($s$), nine ($n$) or twenty ($t$)).

In [ ]:
# In class

Multivariate generating functions: Takeaway¶

Actually, there is no limit on the number of variables one can define, but it might be hard to keep track of all the information we are interested in.

Getting the generating function¶

I think you have seen two types of recursions to play with generating functions: a simpler one is linear recursion, and another one uses the symbolic method.

Catalan generating function using the symbolic method¶

In Hans's lecture on Monday, he did a long computation to get the Catalan numbers from the symbolic method.

Let $F$ be the generating function for Dyck paths, and let $x$ track the number of pairs of steps (one up, one down). Using the first return, we saw that the generating function $$C = \sum_{P \text{ Dyck Paths}} x^{\# \text{up-steps}(P)}$$ satisifes the equation $C = 1 + xC^2$.

Pass this generating function to Sage.

In [19]:
# In class

Now, let's solve this generating function. As we saw in Hans's lecture, there are two answers, but only one has nice coefficients.

In [21]:
# In class
solve(sys, C)
Out[21]:
[C == -1/2*(sqrt(-4*x + 1) - 1)/x, C == 1/2*(sqrt(-4*x + 1) + 1)/x]
In [22]:
sol = _

There are two solutions, only one of which is right. Expand them to figure out how to do that.

In [19]:
# In class

The latter has a lot of negative numbers, so we should avoid it. Then, the Dyck Paths would be counted by the generating function 'C1' here,

In [32]:
C1 = sol[0].rhs()
C1.series(x, 10)
Out[32]:
1 + 1*x + 2*x^2 + 5*x^3 + 14*x^4 + 42*x^5 + 132*x^6 + 429*x^7 + 1430*x^8 + 4862*x^9 + Order(x^10)
In [22]:
oeis([1, 1, 2, 5, 14, 42, 132, 429, 1430])  # Requires internet.
Out[22]:
0: A000108: Catalan numbers: C(n) = binomial(2n,n)/(n+1) = (2n)!/(n!(n+1)!).
1: A120588: G.f. is 1 + x*c(x), where c(x) is the g.f. of the Catalan numbers (A000108).
2: A211216: Expansion of (1-8*x+21*x^2-20*x^3+5*x^4)/(1-9*x+28*x^2-35*x^3+15*x^4-x^5).

We can inspect the OEIS entry A000108 to see that what we got is indeed the generating function for the Catalan numbers.

In [23]:
oeis('A000108').formulas()[0:4]
Out[23]:
0: a(n) = binomial(2*n, n)/(n+1) = (2*n)!/(n!*(n+1)!) = A000984(n)/(n+1).
1: Recurrence: a(n) = 2*(2*n-1)*a(n-1)/(n+1) with a(0) = 1.
2: Recurrence: a(n) = Sum_{k=0..n-1} a(k)a(n-1-k).
3: G.f.: A(x) = (1 - sqrt(1 - 4*x)) / (2*x), and satisfies A(x) = 1 + x*A(x)^2.

Generating functions: What Sage can and cannot do.¶

What Sage easily does for us:

  • Taking the recurrence
  • Solving it, giving the generating function
  • Expanding it; giving any (finite) number of terms

What Sage struggles with; it's something you can do, but it requires effort and knowledge:

  • Giving a closed form for the coefficients
  • Expanding infinite sums

Some other computer algebra systems do it better than SageMath. Some optional packages can be added to SageMath to do it, but it is still not that simple.

Linear recurrence¶

The simpler case is a linear recurrence: something of the form $A_{n} = a\cdot a_{n-1}+b\cdot a_{n-2}+c\cdot a_{n-3}+\ldots$.

A linear recurrence is also called a C-finite sequence. Hence, we can explore the space of CFiniteSequences.

Fibonacci numbers¶

Recall that the Fibonacci numbers are defined using $F_0 = 1$, $F_1 =1$ and $F_n = F_{n-1}+F_{n-2}$ for $n\geq 2$.

In [24]:
# Define the ring of C-finite sequences over the rationals; a C-finite sequence is a recurrence
C.<x> = CFiniteSequences(QQ)

# Give the recurrence to SageMath. This code can be modified for other recurrences.
fibo = C.from_recurrence([1,1],[1,1]) # Recurrence F_{n+2} = F_{n+1} + F_n with v_0 = v_1 = 1
fibo
Out[24]:
C-finite sequence, generated by -1/(x^2 + x - 1)

We got a generating function for the sequence. Can we get other informations, like the coefficients or the partial fraction decomposition?

In [27]:
fibo.series(10).coefficients()
Out[27]:
[1, 1, 2, 3, 5, 8, 13, 21, 34, 55]
In [28]:
fibo[1000]
Out[28]:
70330367711422815821835254877183549770181269836358732742604905087154537118196933579742249494562611733487750449241765991088186363265450223647106012053374121273867339111198139373125598767690091902245245323403501

We can even get easily the partial fraction decomposition...

In [27]:
# In class
Out[27]:
(0, [-1/(x^2 + x - 1)])

Here, however, the denominator cannot we factored over $\mathbb{Q}$, so it is already decomposed.

In [33]:
factor(x^2+x-1)
Out[33]:
x^2 + x - 1
In [34]:
# In comparison...
factor(x^2-1)
Out[34]:
(x - 1) * (x + 1)

Decomposition over another field

It is possible to decompose it over other fields than $\mathbb{Q}$, but we must let SageMath know that we want to do this.

In [34]:
R.<x> = PolynomialRing(RR)  # CC stands for complex, RR for real, QQ for rational, ZZ for integers, QQbar for algebraic, etc. 
factor(x^2+x-1)
Out[34]:
(x - 0.618033988749895) * (x + 1.61803398874989)

There is even a way to 'guess' a generating function from a few coefficients. However, the sequence must be a linear recurrence.

In [35]:
C.guess([1,2,4,8,16,32])
Out[35]:
C-finite sequence, generated by -1/2/(x - 1/2)

For this to work, we must feed the 'guess' function enough terms.

In [37]:
# I secretly generated a sequence from a recurrence. Here are the first terms.
mystery = _.coefficients(); mystery
Out[37]:
[7, 9, 23, 76, 167, 434, 1148, 2851, 7317, 18759]
In [38]:
C.guess([7, 9, 23, 76, 167, 434])
Out[38]:
0
In [39]:
C.guess([7, 9, 23, 76, 167, 434, 1148, 2851, 7317])
Out[39]:
C-finite sequence, generated by (-2/5*x - 7/5)/(x^3 + 2/5*x^2 + 1/5*x - 1/5)

Now that we have the sequence, we can ask for the recurrence. This is exactly what I inputed SageMath (in secret)!

In [40]:
_.recurrence_repr()
Out[40]:
'homogeneous linear recurrence with constant coefficients of degree 3: a(n+3) = a(n+2) + 2*a(n+1) + 5*a(n), starting a(0...) = [7, 9, 23]'

Takeaway¶

  • SageMath knows a lot of tools to make our life with generating functions simpler.
  • SageMath has a lot of things here and there for generating functions, but would benefit from a uniform framework and a greater ease of use!

Open problem: Make a uniform framework for generating functions in SageMath!