This is an extended version of my entry in the Lockdown Mathoff at the Aperiodical

Binet’s formula ((As ordained by Stigler’s law (which is attributed to Merton), Binet’s formula was known at least a century before Binet wrote about it)) 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 argument ((Ridiculously pedantic, you say? I refer the reader to the comment about Haskell being a very mathematical language.)) . 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.

* The code is available here.