Efficiently Evaluating Recursively Defined Functions

In this blog post, I will share a technique I have used in the past to evaluate linear recurrence relations efficiently using a little bit of linear algebra. As a motivating example, assume you have a sequence of numbers that you want to know more about: 1, 2, 3, 7, 24, 100, 448, 2051, 9444, 43549, 200889, 926770, 4275598, 19725314, 91002117, 419835526. Can you spot a pattern (hint: it's a 3rd order linear recurrence relation with constant term and a dependence on \( n \), the index into the sequence)? Once you figure it out, check out my other blog post on finding recurrence relations from integer sequences to see if you came up with the same recurrence relation. By this point, you have figured out \( f(0) = 1 \), \( f(1) = 2 \), \( f(2) = 3 \), and \( f(n) = 5 f(n-1) - 2f(n-2) + f(n-3) - 2n + 1 \) for \( n \geq 3 \). Now let's say I want top compute the last \( 8 \) digits of \( f(10^{18}) \). A naive solution would be to write a for loop, but that would many orders of magnitude too slow on modern computer architecture. In this rest of this blog post, I will outline an approach that I take to overcome this problem, and will provide a Haskell implementation of this approach. While the simple simple for loop will ultimately not be sufficiently efficient, it will serve as a good starting point, so here is a python implementation that would compute \( f(n) \) for \( n \geq 3 \) using constant memory

f0 = 3
f1 = 2
f2 = 1
for n in range(3, N+1):
  f = 5*f0 - 2*f1 + f2 - 2*n + 1
  f2 = f1
  f1 = f0
  f0 = f
print f0

Now we can create an efficient algorithm by collecting the variables in the program into a vector \( [f0, f1, f2, n, 1] \) and defining a linear transformation that can go from one vector to the next: \( [f0, f1, f2, n, 1] \rightarrow [f0', f1', f2', n+1, 1] \). From linear algebra, we know that we can think of a matrix as a linear transformation. To utilize this idea, we need to fill in the \( 5 \times 5 \) matrix below so that it satisfies the equation below. $$ \begin{bmatrix} & & & & \\ & & & & \\ & & & & \\ & & & & \\ & & & & \end{bmatrix} \begin{bmatrix} f(n-1) \\ f(n-2) \\ f(n-3) \\ n \\ 1 \end{bmatrix} = \begin{bmatrix} f(n) \\ f(n-1) \\ f(n-2) \\ n+1 \\ 1 \end{bmatrix} $$ It's fairly easy to deduce the entries in this matrix, as shown below: $$ \begin{bmatrix} 5&-2&1&-2&1\\1&0&0&0&0\\0&1&0&0&0\\0&0&0&1&1\\0&0&0&0&1 \end{bmatrix} \begin{bmatrix} f(n-1) \\ f(n-2) \\ f(n-3) \\ n \\ 1 \end{bmatrix} = \begin{bmatrix} f(n) \\ f(n-1) \\ f(n-2) \\ n+1 \\ 1 \end{bmatrix} $$ Starting with the initial vector (let's call it \(v\)) \( [f(2), f(1), f(0), 3, 1] = [3, 2, 1, 3, 1] \), we can multiply by our matrix (lets call it \(A\)) to get the next vector: $$ \begin{bmatrix} 5&-2&1&-2&1\\1&0&0&0&0\\0&1&0&0&0\\0&0&0&1&1\\0&0&0&0&1 \end{bmatrix} \begin{bmatrix} 3 \\ 2 \\ 1 \\ 3 \\ 1 \end{bmatrix} = \begin{bmatrix} 7 \\ 3 \\ 2 \\ 4 \\ 1 \end{bmatrix} $$ We can keep doing this iteratively to get the next vector. $$ \begin{bmatrix} 5&-2&1&-2&1\\1&0&0&0&0\\0&1&0&0&0\\0&0&0&1&1\\0&0&0&0&1 \end{bmatrix} \begin{bmatrix} 7 \\ 3 \\ 2 \\ 4 \\ 1 \end{bmatrix} = \begin{bmatrix} 24 \\ 7 \\ 3 \\ 5 \\ 1 \end{bmatrix} $$ $$ \begin{bmatrix} 5&-2&1&-2&1\\1&0&0&0&0\\0&1&0&0&0\\0&0&0&1&1\\0&0&0&0&1 \end{bmatrix} \begin{bmatrix} 24 \\ 7 \\ 3 \\ 5 \\ 1 \end{bmatrix} = \begin{bmatrix} 100 \\ 24 \\\ 7 \\ 6 \\ 1 \end{bmatrix} $$ Now what happens if you continue this process \( k \) times? In other words we are multiplying by \( A \) \(k\) times, or equivalently raising \( A \) to the power of \( k \). Since matrix multiplication is associative, \( A (A v) = (A A) v \), so repeating this \( k \) times is equivalent to evaluating the following expression: $$ \begin{bmatrix} 5&-2&1&-2&1\\1&0&0&0&0\\0&1&0&0&0\\0&0&0&1&1\\0&0&0&0&1 \end{bmatrix}^k \begin{bmatrix} f(2) \\ f(1) \\ f(0) \\ 3 \\ 1 \end{bmatrix} = \begin{bmatrix} f(2+k) \\ f(1+k) \\ f(k) \\ 3+k \\ 1 \end{bmatrix} $$ You might be wondering how evaluating this is any more efficient than using the for loop, as it still needs to multiply the initial vector by \( A \) \( k \) times. Luckily, there is a well known efficient exponentiation method that can compute \( x^n \) in \( O(\log{n}) \). While this method is normally introduced using real numbers or integers, it can be used for matrices in the same way. I'm not going to go into details about how it works, but I will say that it utilizes the fact that \( x^{2n+1} = x \cdot x^{2n} \) and \( x^{2n} = (x^n)^2 \). Using this efficient exponentiation method, we know \( f(n) \) is the middle term of the vector \( A^n v \). The Haskell ecosystem makes it very easy to write this sort of program. I will be using two modules: matrix to carry out the matrix operations, and data-modular to perform the operations mod \( 10^8 \) (because we wanted the last 8 digits of the answer). The resulting program is very simple:

{-# LANGUAGE DataKinds, TypeOperators #-}

import Data.Modular
import Data.Matrix

type IntMod = Int `Mod` 100000000

f :: Integral a => Integer -> a
f n = (_A^n * v) ! (3,1)               
  where _A = fromLists [[5,-2,1,-2,1], [1,0,0,0,0], [0,1,0,0,0], [0,0,0,1,1], [0,0,0,0,1]]
        v = fromLists [[3],[2],[1],[3],[1]]

We can verify that this function works as expected interactively, then use it to evaluate \( f(10^{18}) \) as desired.

Ryan> map f [1..10]
[2,3,7,24,100,448,2051,9444,43549,200889]
Ryan> f (10^18) :: IntMod
78847092

so our final answer is \( f(10^{18}) \equiv 78847092 \mod 10^8 \). One thing that I really like about Haskell is that you can have polymorphic output. For example, the output of f can be anything of type Integral, so it could be Int (64 bit int) if we want efficiency, Integer (arbitrary precision int), or IntMod (type of Int that works in \( \mathbb{Z}_{10^8} \) in our case. Then based on context the function returns the type that it needs to so that it will fit in with the rest of the program. In my case, it defaulted to Integer when I called map f [1..10], but I was able to coerce it to my own custom numeric type when I called f (10^18) :: IntMod . This language features provides really nice flexibility in these types of situations, and it's only one many reasons why I enjoy programming in Haskell so much. I hope that you found this blog post useful. This idea has helped me solve a handful of project euler problems, and maybe it can help you too some day.