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!
Partitions(10, max_part=2).cardinality()
6
We could also check how many that is for all the numbers between 10 and 20.
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?
f(x) = (-x^2 + x + 1)/(x^3 - x^2 - x + 1)
f.series(x,21)
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.
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?
S = (1/((1-x^6)*(1-x^9)*(1-x^20))).series(x,100); S
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.
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]
How many ways are there to make 60 chicken nuggets?
How many of those ways use $7$ packs?
# 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 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 class
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.
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.
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 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 class
solve(sys, C)
[C == -1/2*(sqrt(-4*x + 1) - 1)/x, C == 1/2*(sqrt(-4*x + 1) + 1)/x]
sol = _
There are two solutions, only one of which is right. Expand them to figure out how to do that.
# 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,
C1 = sol[0].rhs()
C1.series(x, 10)
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)
oeis([1, 1, 2, 5, 14, 42, 132, 429, 1430]) # Requires internet.
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.
oeis('A000108').formulas()[0:4]
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.
What Sage easily does for us:
What Sage struggles with; it's something you can do, but it requires effort and knowledge:
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.
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
.
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$.
# 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
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?
fibo.series(10).coefficients()
[1, 1, 2, 3, 5, 8, 13, 21, 34, 55]
fibo[1000]
70330367711422815821835254877183549770181269836358732742604905087154537118196933579742249494562611733487750449241765991088186363265450223647106012053374121273867339111198139373125598767690091902245245323403501
We can even get easily the partial fraction decomposition...
# In class
(0, [-1/(x^2 + x - 1)])
Here, however, the denominator cannot we factored over $\mathbb{Q}$, so it is already decomposed.
factor(x^2+x-1)
x^2 + x - 1
# In comparison...
factor(x^2-1)
(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.
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)
(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.
C.guess([1,2,4,8,16,32])
C-finite sequence, generated by -1/2/(x - 1/2)
For this to work, we must feed the 'guess' function enough terms.
# I secretly generated a sequence from a recurrence. Here are the first terms.
mystery = _.coefficients(); mystery
[7, 9, 23, 76, 167, 434, 1148, 2851, 7317, 18759]
C.guess([7, 9, 23, 76, 167, 434])
0
C.guess([7, 9, 23, 76, 167, 434, 1148, 2851, 7317])
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)!
_.recurrence_repr()
'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]'
Open problem: Make a uniform framework for generating functions in SageMath!