--- execution: timeout: 300 jupytext: notebook_metadata_filter: all text_representation: extension: .md format_name: myst format_version: 0.13 jupytext_version: 1.13.0 kernelspec: display_name: Python 3 (ipykernel) language: python name: python3 --- ```{code-cell} ipython3 :tags: [hide-cell] import mmf_setup; mmf_setup.nbinit() import logging logging.getLogger("matplotlib").setLevel(logging.CRITICAL) %matplotlib inline import numpy as np, matplotlib.pyplot as plt ``` (sec:RecursionAndInvariants)= # Recursion and Invariants Consider the famous Fibonacci numbers defined by \begin{align*} F_0 &= 0,\\ F_1 &= 1,\\ F_{n} &= F_{n-1} + F_{n-2}, \qquad \text{for } n>1. \end{align*} We code this directly in terms of recursion: ```{code-cell} def fib0(n): """Return the n'th Fibonnacci number.""" if n == 0: return 0 if n == 1: return 1 return fib0(n-1) + fib0(n-2) [fib0(n) for n in range(15)] ``` If you try timing these, however, you will notice a disturbing trend: :::{margin} Challenge! If you play a little, you should find that $t \propto 1.6^{n}$. Where does $1.6$ come from? Can you provide some uncertainties for this? ::: ```{code-cell} :tags: [hide-input] from timeit import timeit ns = np.arange(12) number = 20 ts = np.array( [timeit("fib(n)", globals=dict(fib=fib0, n=n), number=number)*1000/number for n in ns]) fig, axs = plt.subplots(1, 2, figsize=(8,3), constrained_layout=True) for ax in axs: ax.plot(ns, ts, '+-') ax.set(xlabel="$n$", ylabel="$t$ [ms]") axs[1].set_yscale('log'); m, b = np.polyfit(ns[4:], np.log10(ts)[4:], deg=1) axs[1].plot(ns, 10**(m*ns + b), label=fr'$t={1000*10**b:.2g}μs×{10**m:.2g}^{{n}}$'); axs[1].legend(); ``` ```{code-cell} :tags: [margin] t_h = 0.13e-6 * (1.6**55) / 60**2 print(f"{t_h:.2g}h on my computer") ``` If this trend continues -- and it seems pretty likely -- then `fib(55)` will take about 6 hours to compute. :::{admonition} Solution :class: dropdown How many times must `fib0` recurse? Call this $S(n)$: we have the base-cases $S(0) = S(1) = 1$. Then, by recursion, we must have \begin{gather*} S(N) = S(N-1) + S(N-2), \qquad \text{for } N > 1. \end{gather*} This means that that $S(N)$ is just the Fibonacci series shifted by one: \begin{gather*} S(N) = F_{N+1}. \end{gather*} Thus, the time grows in the same way as the Fibonacci series. Can we find an explicit formula? We give two solutions here. 1. The Fibonacci numbers satisfy a **2-term homogeneous and linear recurrence relation**: \begin{gather*} F_{n} - F_{n-1} - F_{n-1} = 0. \end{gather*} There is an exact analogy with **2nd-order homogeneous and linear ordinary differential equations** (ODEs). The following strategy generally solves these. As with ODEs, we guess a solution of the form $F_{n} = e^{an}$. Plugging this in, we have \begin{gather*} e^{an} = e^{a(n-1)} + e^{a(n-2)} = e^{an}(e^{-a} + e^{-2a}) = e^{an}\bigl(e^{-a} + (e^{-a})^2\bigr). \end{gather*} Letting $x = e^{a}$ so that $F_{n} = e^{an} = x^n$, this becomes a quadratic equation: \begin{gather*} 1 = \frac{1}{x} + \frac{1}{x^2}, \qquad x^2 - x - 1 = 0, \\ x_{\pm} = \frac{1\pm \sqrt{5}}{2}. \end{gather*} Since this is a two-term recurrence, there are two independent solutions, just as with a 2nd order ODE. The general solution can be expressed as \begin{gather*} F_{n} = c_+ x_+^n + c_-x_-^n \end{gather*} where $c+{\pm}$ ate chosen to match the initial conditions $F_{0} = c_{+} + c_{-} = 0$, and $F_{1} = c_+x_+ + c_-x_- = 1$. Solving this gives the complete solution: \begin{gather*} F_{n} = \frac{x_+^n - x_-^n}{\sqrt{5}} = \frac{\varphi^{n} - (-\varphi)^{-n}}{\sqrt{5}}. \end{gather*} Note that $x_+ = \varphi$ is the famous [Golden ratio][]: \begin{gather*} x_+ = \frac{1+\sqrt{5}}{2} = \varphi = 1.618\dots,\qquad x_- = \frac{-1}{\varphi} = 1-\varphi = -.618\dots. \end{gather*} Asymptotically, $F_{n} \sim \varphi^{n}$, explaining the asymptotic time-dependence of $1.6^{n}$ seen in the cost of `fib0`. In fact: \begin{gather*} \DeclarePairedDelimiter{\round}⌊⌉ \DeclarePairedDelimiter{\round}\lfloor\rceil F_{n} = \round[\Big]{\varphi^{n}/\sqrt{5}} \end{gather*} where $\round{x}$ means "round $x$ to the nearest integer. 2. Alternatively (see {cite}`Graham:1994`), one can form the generating fuction \begin{alignat*}{6} f(z) &= &&F_0 + &&F_1 z + &&F_2 z^2 + && F_3 z^3 + \cdots\\ zf(z) &= && &&F_0z + &&F_1 z^2 + && F_2 z^3 + \cdots\\ z^2f(z) &=&& && &&F_0 z^2 + && F_1 z^3 + \cdots. \end{alignat*} Subtracting the second two from the first, and using the recurrence $F_{n} - F_{n-1} - F_{n-2} = 0$, we have \begin{gather*} (1-z-z^2)f(z) = F_0 + (F_1 - F_0) z = z, \\ f(z) = \frac{z}{1-z-z^2} = \sum_{n=0}^{\infty} F_{n}z^n. \end{gather*} Now we simply need to find the [Maclaurin series][] of this function, which can be done using partial fractions and then geometric series: \begin{gather*} z^2 + z - 1 = (z + \varphi)(z + 1-\varphi), \qquad \varphi = \frac{1 + \sqrt{5}}{2}\\ f(z) = \frac{z}{1-z-z^2} = \frac{-z}{(z + \varphi)(z + 1 - \varphi)} = \frac{z/\sqrt{5}}{z+1-\varphi} - \frac{z/\sqrt{5}}{z + \varphi}\\ = \frac{1}{\sqrt{5}} \left(\sum_{n=1}^{\infty} \frac{z^n}{(\varphi-1)^n} - \sum_{n=1}^{\infty} \frac{z^n}{(-\varphi)^n}\right)\\ = \sum_{n=1}^{\infty} \underbrace{\frac{1}{\sqrt{5}} \left(\frac{1}{(\varphi-1)^n} - \frac{1}{(-\varphi)^n}\right) }_{F_{n}}z^{n}. \end{gather*} Using the fact that $\varphi-1 = \phi^{-1}$, we have the same result as above: \begin{gather*} F_{n} = \frac{\varphi^{n} - (-\varphi)^{-n}}{\sqrt{5}}. \end{gather*} ::: ```{code-cell} :tags: [margin, hide-input] # Check xp = phi = (1+np.sqrt(5))/2 xm = -1/phi x_ = np.array([xp, xm]) assert np.allclose(x_**2 - x_, 1) print(xp, xm) ``` :::{margin} There is actually a really simple fix using [memoization][] with {py:func}`functools.cache`: ::: ```{code-cell} :tags: [margin] from functools import cache @cache def fib0a(n): """Return the n'th Fibonnacci number.""" if n == 0: return 0 if n == 1: return 1 return fib0a(n-1) + fib0a(n-2) %timeit fib0a(100) ``` :::{margin} The {py:func}`functools.cache` **decorator** adds a cache that stores previous values so that they do not need to be recomputed, leading to linear runtime at the cost of linear space complexity. ::: ## From Recursion to Iteration One way to improve this is to accumulate the results as we go along. This turns our recursive program into an iterative program. ```{code-cell} def fib1(n): F0 = 0 F1 = 1 while n > 0: F0, F1 = F1, F0 + F1 n -= 1 return F1 ``` Is this correct? Perhaps we should write a test: ```{code-cell} def test_fib(fib): res = list(map(fib, range(15))) exact = [0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377] if res == exact: return True else: print(f"{res} != {exact}") return False test_fib(fib0) ``` Our test seems to work with `fib0`. What about `fib1`? ```{code-cell} test_fib(fib1) ``` Looking at the output, we can probably guess how to fix this: ```{code-cell} def fib1a(n): F0 = 0 F1 = 1 while n > 0: F0, F1 = F1, F0 + F1 n -= 1 return F0 test_fib(fib1a) ``` But how can we code in a way that we know our code is correct? A general approach for this strategy involves **invariants**. ## Invariants Here is how we might structure our code, adding **invariants** so that we can ascertain the correctness of our program: ```{code-cell} def fib2(n): # Invariant: F0 = F_{m}, F1 = F_{m+1} m = 0 F0 = 0 F1 = 1 # Invariant satisfied: F0 = F_{0} = 0, F1 = F_{1} = 1 # Check that n is a positive integer... assert n >= 0 and n % 1 == 0 while m < n: # At the start of the loop we assume the invariant is satisfied: # F0 = F_{m}, F1 = F_{m+1} F0, F1 = F1, F0 + F1 # F0 = F_{m+1}, F1 = F_{m} + F_{m+1} m += 1 # F0 = F_{m}, F1 = F_{m-1} + F_{m} # At the end of the loop, we check the invarient: # F0 = F_{m}, F1 = F_{m-1} + F_{m} = F_{m+1} # Invariant satisfied! # The loop will terminate when m == n (as long as n is an integer) # thus, F0 contains F_{m} return F0 # Now let's check to be sure test_fib(fib2) ``` [Golden ratio]: [Maclaurin series]: [Memoization]: