Binet’s formula and Haskell

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

Binet’s formula ¹ is a lovely way to generate the $n$ th Fibonacci number, $F_{n}$ . If $ϕ = \frac{1}{2} (\sqrt{5} + 1)$ , then

F_n = f r a c p h i^{n} - (- p h i)^{- n} s q r t 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 $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}) = (a c + 5 b d) + (a d + b c) \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) (x^{(n - 1) / 2})$ - 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 ² . 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 $2 B \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.

Footnotes:

1. 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

2. Ridiculously pedantic, you say? I refer the reader to the comment about Haskell being a very mathematical language.