Written by Colin+ in arithmetic, haskell.

*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}}$$

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?

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

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 **$**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.

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

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.