APM_2F005 - Algorithms For Discrete Mathematics Bachelor 2
Lecturer: Lucas Gerin

Symbolic computing 2: Generating functions with SymPy

Table of contents

Basics of generating functions

Let us first explain how we will manipulate generating functions with SymPy. We consider the example of

$$ f(x)=\frac{1}{1-2x}=1+2x+4x^2+8x^3+16x^4+ \dots $$
We first introduce a symbolic variable $x$ and an expression $f$ as follows:

One can extract $n$-th coefficient as follows:

Exercise 1. Fibonacci generating function

Let $f_n$ be the Fibonacci sequence defined by $f_0=1$, $f_1=1$ and for all $n\geq 2$, $$ f_n=f_{n-1}+f_{n-2}. $$

  1. Write a recursive function `Fibonacci(n)` which returns the $n$-th Fibonacci number.
  2. Let $F$ be the generating function associated to the sequence $(f_n)$. Show that $$ F(x)=\frac{1}{1-x-x^2}. $$
  3. Write another function `FibonacciGF(n)` which also returns the $n$-th Fibonacci number by extracting the $n$-th coefficient in $F(x)$.
    Let $n\geq 2$ and $x>0$: \begin{align*} \sum_{n\geq 2} f_n x^n &= \sum_{n\geq 2} f_{n-1} x^n + \sum_{n\geq 2} f_{n-2}x^n \\ F(x)-1-x &= \sum_{p\geq 1} f_{p} x^{p+1} + \sum_{q\geq 0} f_{q}x^{q+2} \qquad \text{(we put }n-2=p\text{ and }n-1=q)\\ F(x)-1-x &= x\sum_{p\geq 1} f_{p} x^{p} + x^2\sum_{q\geq 0} f_{q}x^{q}\\ F(x)-1-x &= x(F(x)-1) + x^2F(x). \end{align*} We solve this last equation with the following script:

Exercise 2. Recurrence of order two and asymptotics

Let $j_n$ be defined by

\begin{align} j_0&=0, \notag\\ j_1&=1, \notag\\ j_2&=2,\notag\\ j_n&=2j_{n-2}+5 \quad (\text{ for every }n\geq 3). \qquad \text{(#)} \end{align}
(Theory) Find the expression for the generating function $J(x)$ of the $j_n$'s.
  1. We multiply eq. $(\#)$ by $x^n$ and sum the resulting expression for all $n\geq 3$: \begin{align*} \sum_{n\geq 3} j_n x^n &= 2 \sum_{n\geq 3} j_{n-2} x^n + 5\sum_{n\geq 3} x^n \\ J(x)-j_1x-j_2x^2 &= 2 \sum_{p\geq 1} j_{p} x^{p+2} + 5(x^3+x^4+x^5+\dots) \qquad \text{(we put }n-2=p)\\ J(x)-x-2x^2 &= 2 x^2J(x) + 5\frac{x^3}{1-x}. \end{align*} We solve this last equation with the following script:
We find $$ J(x)= \frac{x \left(3 x^{2} + x + 1\right)}{2 x^{3} - 2 x^{2} - x + 1} $$
  1. Write a function which extracts the $n$-th coefficient in $J(x)$.
    2. Compare your results with a recursive function which computes the $j_n$'s.
  1. What is the radius of convergence of $J(x)$? (You can ask help to SymPy.)
  2. What does it imply for the asymptotic behaviour of $j_n$? (Apply the "exponential growth formula", that we saw in class.)
  1. $J$ is a rational fraction, it explodes exactly when the denominator is zero. The radius of convergence $\rho$ of $J(x)$ is therefore the smallest positive solution of $2 \rho^{3} - 2 \rho^{2} - \rho + 1=0$, i.e. $$ \rho= \sqrt{2}/2=1/\sqrt{2} $$ according to the previous script.
  2. The exponential growth formula says that for every $\varepsilon$, for $n$ large enough $$ C_1(\sqrt{2}-\varepsilon)^n =C_1(1/\rho-\varepsilon)^n \leq J_n \leq C_2(1/\rho+\varepsilon)^n =C_2 (\sqrt{2}+\varepsilon)^n. $$ (The exponential growth formula only ensures that the left inequality holds for infinitely many $n$'s.)
With a plot, find an approximation of $r$ such that $j_n$ grows like $\text{const}\times r^n$. Compare with the previous question.
We plot $n\mapsto (j_n)^{1/n}$ which seems to converge to $\sqrt{2}$. This is consistent with the exponential growth formula.

Exercise 3. A pair of generating functions

Let $a_n,b_n$ be defined by $a_1=b_1=1$ and, for every $n\geq 1$, \begin{equation} \begin{cases} a_{n+1}&= a_n +2b_n,\\ b_{n+1}&= a_n +b_n. \end{cases} \tag{&} \end{equation}

  1. Find a $2\times 2$ system whose solutions are $A(x),B(x)$ , where $A,B$ are the generating functions of sequences $(a_n)_{n\geq 1},(b_n)_{n\geq 1}$. (Coefficients of this system should depend on $x$.)
  2. Solve this system with `solve` and write a script which uses function $A$ to return $a_1,\dots,a_{20}$.
  1. We multiply both terms of both equations in eq.(&) by $x^{n+1}$ and take the sum for $n\geq 1$. We obtain $$ \begin{cases} \sum_{n\geq 1}a_{n+1}x^{n+1}&= \sum_{n\geq 1}a_{n}x^{n+1} +2\sum_{n\geq 1}b_{n}x^{n+1},\\ \sum_{n\geq 1}b_{n+1}x^{n+1}&= \sum_{n\geq 1}a_{n}x^{n+1} +\sum_{n\geq 1}b_{n}x^{n+1}.\\ \end{cases} $$ We find $$ \begin{cases} A(x)-x&= xA(x) +2xB(x),\\ B(x)-x&= xA(x)+xB(x).\\ \end{cases} $$ In the above script we solve this sytem of equations. We find $$ \begin{cases} A(x)&=- \frac{x \left(x + 1\right)}{x^{2} + 2 x - 1}\\ B(x)&=- \frac{x}{x^{2} + 2 x - 1} \end{cases}. $$

Automatic decomposition of fractions

In class we saw that for GFs it is useful to decompose fractions like this: $$ \frac{1-x+x^2}{(1-2x)(1-x)^2} = \frac{3}{1-2x} - \frac{1}{1-x} - \frac{1}{(1-x)^2}. $$ Here are examples on how to do that with SymPy.

Exercise 3 (continued)

The goal of the exercise is to find coefficients $\alpha,\beta,a,b,c$ such that $$ A(x)=\frac{a}{x-\alpha}+\frac{b}{x-\beta}+c, $$ where $$ A(x)=\frac{-x \left(x + 1\right)}{x^2+2x-1} $$ was defined in the previous exercise.
  1. (Theory) Compute $\lim_{x\to +\infty} A(x)$ and deduce $c$.
  2. (Theory + SymPy) Use`SymPy` to find coefficients $\alpha,\beta$.
  3. (Theory + SymPy) Use`SymPy` again to find coefficients $a,b$.
Question 1. Taking the limit $(x\to +\infty)$ in the equation $$ \frac{a}{x-\alpha}+\frac{b}{x-\beta}+c=\frac{-x \left(x + 1\right)}{x^2+2x-1} $$ yields $0+0+c=-1$.
Question 2.
The short code above shows that $$ A(x)= \frac{-x \left(x + 1\right)}{x^2+2x-1}= \frac{-x \left(x + 1\right)}{(x-(\sqrt{2}-1))(x-(-\sqrt{2}-1))} $$ so we must have that $\alpha=\sqrt{2}-1$ and $\beta=-\sqrt{2}-1$ (or the contrary) to ensure that $A$ has the proper definition domain.
In order to find $a,b$ we solve with the previous code the system $$ \begin{cases} A(0)&=\frac{a}{0-\alpha}+\frac{b}{0-\beta}+(-1)\\ A(1)&=\frac{a}{1-\alpha}+\frac{b}{1-\beta}+(-1)\\ \end{cases} $$ with $\alpha=\sqrt{2}-1$ and $\beta=-\sqrt{2}-1$. We find `{b: 1/2 + sqrt(2)/2, a: -sqrt(2)/2 + 1/2}` and finally $$ A(x)=-1 + \frac{- \frac{\sqrt{2}}{2} + \frac{1}{2}}{x - \sqrt{2} + 1} + \frac{\frac{1}{2} + \frac{\sqrt{2}}{2}}{x + 1 + \sqrt{2}}. $$
(Bonus: Theory) Deduce a proof of the formula $$ a_n=\frac{1}{2} \left(1 + \sqrt{2}\right)^{n} + \frac{1}{2} \left(- \sqrt{2} + 1\right)^{n} $$
(Hint: Use the formula $$ \frac{1}{x-\rho}= -\frac{1/\rho}{1-x/\rho} = -1/\rho \sum_{n\geq 0}x^n(1/\rho)^n. \tag{E} $$
We have that $$ A(x)=-1 + \frac{- \frac{\sqrt{2}}{2} + \frac{1}{2}}{x - \alpha} + \frac{\frac{1}{2} + \frac{\sqrt{2}}{2}}{x -\beta} $$ with $\alpha=\sqrt{2}-1$ and $\beta=-\sqrt{2}-1$. Using equation (E) with $\rho= \alpha $ and $\rho=\beta$ we obtain \begin{align*} A(x)&=-1 + \frac{1-\sqrt{2}}{2} (-1/\alpha \sum_{n\geq 0}x^n(1/\alpha)^n)+ \frac{1+\sqrt{2}}{2}(-1/\rho \sum_{n\geq 0}x^n(1/\beta)^n)\\ &= -1 + \frac{1}{2} \sum_{n\geq 0}x^n(1/\alpha)^n+ \frac{1}{2} \sum_{n\geq 0}x^n(1/\beta)^n\\ &= -1 + \frac{1}{2} \sum_{n\geq 0}x^n(1+\sqrt{2})^n+ \frac{1}{2} \sum_{n\geq 0}x^n(1-\sqrt{2})^n \end{align*} since $1/\alpha=1+\sqrt{2}$ and $1/\beta = 1-\sqrt{2}$. By identification of coefficients we obtain $$ a_n=\frac{1}{2} \left(1 + \sqrt{2}\right)^{n} + \frac{1}{2} \left(- \sqrt{2} + 1\right)^{n}. $$

Exercise 4. A quadratic GF

Let $(g_n)_{n\geq 0}$ be the sequence of non-negative integers such that the generating function $G(x)=\sum_{n\geq 0}g_n x^n$ satisfies the equation $$ G(x)=x+x G(x)^2+1. $$
  1. Compute $g_0,g_1,\dots,g_{20}$.
  2. **(More difficult)** Find the radius of convergence of $G$.
  3. What does it imply for the growth of the sequence $(g_n)$?
  1. According to the above script, the only $G$ solution which has non-negative coefficients is $$ G(x)=\frac{1-\sqrt{1-4x^2-4x}}{2x}. $$ According to SymPy $$ g_{19}=83015133184 ,g_{20}= 37209012224 ,g_{21}=1673660915712 $$
  2. We want to find when $G$ explodes. We might be worried by the fact that the denominator goes to. zero, but actually $G$ is well-defined when $x\to 0$. Indeed, using $\sqrt{1+h}=1+\tfrac{1}{2}h+o(h)$ we obtain $$ G(x)=\frac{1-\sqrt{1-4x^2-4x}}{2x}=\frac{1}{2x}(1-\tfrac{1}{2}4x+o(x))\to 1. $$ The series will explode when what is inside the square root is less than zero, i.e. when $1-4\rho^2-4\rho<0$. According to the code below we have $\rho=\frac{1}{2}(\sqrt{2}-1)$.
  1. Observe that $1/\rho=2+2\sqrt{2}$. The exponential growth formula tells us that for every $\varepsilon> 0$, for large enough $n$ then $$ (2+2\sqrt{2}-\varepsilon)^n \leq g_n \leq (2+2\sqrt{2}+\varepsilon)^n. $$ (The exponential growth formula only ensures that the left inequality holds for infinitely many $n$'s.)

Exercise 5. A continued fraction (Taken from 2022 Test)

We set $$ u_1=\frac{1}{3+1}, \quad u_2=\frac{1}{3+\frac{1}{3+1}},\quad u_3=\frac{1}{3+\frac{1}{3+\frac{1}{3+1}}}, \quad u_4=\frac{1}{3+\frac{1}{3+\frac{1}{3+\frac{1}{3+1}}}}, \dots $$ Using SymPy, write $u_{25}$ as a rational fraction $a/b$.
The sequence $(u_n)$ can be defined by $u_1=1/4$ and $$ u_{n+1}=\frac{1}{3+u_n}. $$ Using the script below allows one finds: $$ u_{25}=\frac{3387464586001}{11188035508324}. $$