This is an extended version of my entry in the Lockdown Mathoff at the Aperiodical
Binet’s formula1 is a lovely way to generate the $n$th Fibonacci number, $F_n$. If $\phi = \frac{1}{2}\left(\sqrt{5} + 1\right)$, then
$$F_n = \frac{ \phi^n - (-\phi)^{-n} }{\sqrt{5}}$$
Haskell and computation
The main reason I’m writing about this is this post. Haskell is something I’ve been meaning to play with for a while - it has a reputation of having a fearsome learning curve, and of being a very mathematical language. I’m not sure what that means, precisely, but it interests me. Implementing Binet’s formula struck me as a nice First Proper Problem to solve in Haskell.
I didn’t want to do it the same way as Gabriel did - I wouldn’t know a monoid from a monorail - but I did take heed of his warning not to do it with floating point arithmetic. Instead, I figured it would be nice to work with elements of $\bb{Q}(\sqrt{5})$, which is a highfalutin’ way of writing ’numbers of the form $a + b\sqrt{5}$, where $a$ and $b$ are integers.
If I rewrite Binet’s formula as $\frac{ (1 + \sqrt{5})^n - (1 - \sqrt{5})^n }{2^n \sqrt{5} }$, I get everything in terms of those nice elements. The question now is, how to work out that denominator efficiently?
Multiplying 2-tuples
I figured that I could use 2-tuples for this - I can represent $a + b\sqrt{5}$ with the tuple (a,b) and work out how to multiply such tuples together. Since $(a + b\sqrt{5})(c + d\sqrt{5}) = (ac + 5bd) + (ad + bc)\sqrt{5}$, I can define a function to multiply two tuples together:
-- Tuple multiplication: (a + b√5)(c + d√5) = (ac + 5bd) + (ad + bc)√5 fmul :: (Integer, Integer) -> (Integer, Integer) -> (Integer, Integer) fmul (a,b) (c,d) = (a*c + 5*b*d, a*d + b*c)
The first line is a comment explaining what I’m doing, like I did above.
The second is a type description, which I’m (generally) a bit hazy on, details-wise; here, it says fmul takes a 2-tuple of Integers, a second 2-tuple of Integers, and returns a third 2-tuple of Integers. (Haskell Integers go as big as you need them to, as opposed to Ints, which are limited to - I think - 64 bits.)
The third line tells the compiler how to compute the function: if you fmul two 2-tuples together, you get a third 2-tuple defined exactly as I set out before.
In principle, I imagine I could figure out the $n$th Fibonacci number just using what I’ve done so far (is it a fold? I’ll check another time), but I wanted to do it cleverly and to explore recursion, and there’s a natural way to work out $x^n$ in general: if $n$ is even, square $x^{n/2}$, and if it’s odd, work out $(x)\left(x^{(n-1)/2}\right)$ - and as long as I know what $x^1$ looks like, it all comes out with relatively few calculations (somewhere in the region of $\log_2(n)$, I think.)
Working out powers
So, it’ll be helpful if I tell Haskell how to square one of my 2-tuples:
-- Tuple squaring (for efficiency) fsq :: (Integer, Integer) -> (Integer, Integer) fsq (a,b) = fmul (a,b) (a,b)
This time, fsq is a function that takes a single 2-tuple and returns a single 2-tuple, and it’s defined in the way you’d expect: you use fmul to multiply the tuple by itself.
Now to implement the recursive powers:
-- Powers of the tuples - divide and conquer fpow :: Integer -> (Integer, Integer) fpow 0 = (1,0) fpow 1 = (1,1) fpow n \| mod n 2 == 0 = fsq $ fpow $ div n 2 \| otherwise = fmul (1, 1) (fsq $ fpow $ div (n-1) 2)
fpow is going to work out the $n$th power of (1,1), so it takes an Integer ( $n$ ) and returns a 2-tuple.
This time, the definition is a bit more complicated: I need to explain what to return in several different cases. The first two are straightforward base cases: (1,1)$^0 = $(1,0) (which I don’t think I need unless someone explicitly calls fpow 0, which is conceivable), and (1,1)$^1 = $(1,1), which I will need.
The grunt work is in the line that starts fpow n. It has a guard next to it, the | symbol, which tells Haskell to do different things under different conditions. Here, mod 2 n == 0 asks “is the number even?” - in which case, I want to square fpow applied to $\frac{n}{2}$ - even though I know $n$ is even, I have to use careful integer division because ‘normal’ division in Haskell returns a Float, and fpow needs an Integer argument2 . The **&dollars;**s? They’re there to separate the arguments, by which I mean “the code didn’t work until I put them in.”
And otherwise means… well, otherwise. Otherwise, I need to multiply (1,1) by the appropriate square - again, the dollars are there to make it work.
Getting the Fibonacci numbers
So far so good. However, I need to work out (1,1)$^n - $(1,-1)$^n$. Luckily, there’s a clever hack to bring to play: if $(a+b\sqrt{5})^n = A + B\sqrt{5}$, then $(a-b\sqrt{5})^n = A - B\sqrt{5}$, and their difference is just $2B\sqrt{5}$ - or double the second element of the 2-tuple. I still have to divide the whole thing by $2^n \sqrt{5}$, though; the $n$th Fibonacci number turns out to be the second element of the 2-tuple divided by $2^{n-1}$. Or, using Haskell’s handy snd to retrieve the second element of a 2-tuple:
-- fib n is the nth Fibonacci number fib :: Integer -> Integer fib n = div (snd $ fpow n) (2 ^ (n-1))
And it’s quick! On my machine, working out the hundred-millionth Fibonacci number (which begins 4737103… and ends …0546875, and is about 20 million digits long) takes a handful of seconds to work out (and considerably longer to print).
One last ‘rather neatly’
Rather neatly, the first element of the 2-tuple, divided by $2^{n-1}$, gives the $n$th Lucas number. That, however, is a story for another day.