We can use recursion to compute the factorial \(n!\) of the integer \(n\), like so.

def f(n):
    if n == 1:
        return 1
    return n * f(n - 1)

print(f(4))

24

Note that, when using recursion, always think about the stopping condition. If you miss a corner case, your code will keep running.

Note further that in the code, we first check the boundary condition, i.e., n == 1, and only when we fell through all these conditions, we call the recursion.

Now we ask how often this function is called to compute for instance \(5!\).

def f(n):
    global no_times_called
    no_times_called += 1
    if n == 1:
        return 1
    return n * f(n - 1)


no_times_called = 0
f(5)
print(no_times_called)

5

Obviously, to compute \(5!\), we have to call the function f also \(5\) times. Of course, in this example it’s trivial, but this does not hold always.

The number of Rabbit pairs increases per generation, supposedly, according to Fibonacci’s rule:

Let’s implement this as a recursion, and see how often the function is called.

def Fibb(n):
    global no_times_called
    no_times_called += 1
    if n == 0 or n == 1:
        return 1
    return Fibb(n - 1) + Fibb(n - 2)


for n in range(20):
    no_times_called = 0
    Fibb(n)
    print(f"{n=:5d}, {no_times_called=:8d}")

n=    0, no_times_called=       1
n=    1, no_times_called=       1
n=    2, no_times_called=       3
n=    3, no_times_called=       5
n=    4, no_times_called=       9
n=    5, no_times_called=      15
n=    6, no_times_called=      25
n=    7, no_times_called=      41
n=    8, no_times_called=      67
n=    9, no_times_called=     109
n=   10, no_times_called=     177
n=   11, no_times_called=     287
n=   12, no_times_called=     465
n=   13, no_times_called=     753
n=   14, no_times_called=    1219
n=   15, no_times_called=    1973
n=   16, no_times_called=    3193
n=   17, no_times_called=    5167
n=   18, no_times_called=    8361
n=   19, no_times_called=   13529

Compared to the recursive algorithm to compute \(n!\), this algorithm clearly explodes.

What if we would store intermediate results rather then recompute them? The concept of storing intermediate results is known as memoization, or cashing, or hashing. Python provides the function cache in the functools library to exactly do this. Let’s see how to use it, and what it does.

from functools import cache


@cache
def Fibb(n):
    global no_times_called
    no_times_called += 1
    if n == 0 or n == 1:
        return 1
    return Fibb(n - 1) + Fibb(n - 2)


no_times_called = 0
Fibb(19)
print(f"{no_times_called=:8d}")

no_times_called=      20

Rather than doing 13529 computations, we now just need 20 calls! Amazing, isn’t it? If you don’t believe how useful memoization is, try your hand at F(300) with and without memoization. ³³Without memoization, you’ll be dead before knowing the answer.

There are some remarks in order.

• There is a closed form solution to compute the \(n\)th Fibonacci number. Of course we could have used that, but here we wanted to see how memoization serves to dramatically reduce numerical work.
• Here we don’t explain how memoization is implemented. It’s not hard, though, and there are excellent explanations on the web.

We like to compute the number of possibilities \(\P{n, N}\) for \(x = (x_1, \ldots x_n)\) such that \(x_i \in \{0,1,\ldots, N\}\) and \(\sum_i x_i \leq N\). It is easy to see that \(\P{n, N}\) satisfies the recursion: \[\P{n, N} = \sum_{i=0}^N \P{n-1, N-i},\] with boundary conditions \(\P{1, N} = N+1\) for all \(N\). Note that by using the summation above, this condition can be replaced by \(\P{0, N} = 1\) for all \(N\).

Computing \(\P{n,N}\) is easy when we use memoization. In fact, we can code the compuation in nearly one line!

from functools import lru_cache

@lru_cache(maxsize=128)
def P(n, N):
   return (n==0) or sum( P(n-1, N-i) for i in range(N+1) )

n=5
N=80
print(P(n = 5, N = 80))

32801517

We throw multiple times with a coin that lands heads with probability \(p\). What is the probability \(\P{n,k}\) to see at least \(k\) heads in row when you throw a coin \(n\) times?

A bit of thinking shows that \(\P{n,k}\) must satisfy the recursion \[\P{n,k} = p^k + \sum_{i=1}^k p^{i-1} q\, \P{n-i,k},\] because it is possible to throw \(k\) times heads from the first throw, but otherwise you throw \(i\), \(i < k\), times a heads, then a tails, after which you have to start all over again.

Reasoning similarly, the expected number times \(\E{n,k}\) to see at least \(k\) heads in row when you throw a coin \(n\) times must satisfy the recursion \[\E{n,k} = p^k(1+\E{n-k,k}) + \sum_{i=1}^k p^{i-1} q\, \E{n-i,k}.\]

p = 0.5
q = 1. - p

@lru_cache(maxsize=128)
def P(n, k):
    """
    probability to see at least k heads in row when a coin is thrown n times
    """
    if n < k:
        return 0
    else:
        return sum(P(n-i,k) * p**(i-1) * q for i in range(1,k+1)) + p**k

@lru_cache(maxsize=128)
def E(n, k):
    """
    expected number of times to see at least k heads in row when a coin is thrown n times
    """
    if n < k:
        return 0
    else:
        tot = sum(E(n-i,k) * p**(i-1) *q for i in range(1,k+1))
        tot += p**k * (1 + E(n-k,k))
        return tot



k = 2

for n in range(k,10):
    print(n, P(n,k), E(n,k))

The next exercise challenges you generalize what we discussed above.

from math import comb
from functools import cache

N = 45

@cache
def T(n, m):
    if n >= N:
        return 0
    res = comb(N, m)
    for k in range(1, min(m, N - n) + 1):
        P = comb(n, m - k)
        P *= comb(N - n, k)
        res += P * T(n + k, m)
    return res / (comb(N, m) - comb(n, m))

print(T(0, 6))

31.497085595869386