9

How would the Fibonacci's closed form code look like in haskell?

enter image description here

Don Stewart
  • 137,316
  • 36
  • 365
  • 468
edgarmtze
  • 24,683
  • 80
  • 235
  • 386
  • Choose a language already; you really expect one answer to give code for all of those? – ildjarn May 17 '11 at 22:07
  • 4
    Just write it yourself. The expression's given right there... do you really need other people to work for you? – erjiang May 17 '11 at 22:08
  • 10
    The formula involves exponentiation, subtraction, and division. Can you please be more specific about which part of that you are having trouble with? – Rob Kennedy May 17 '11 at 22:08
  • @ildjan: Not really... Wwhatever language you want, I added `haskell` because I think it could be easier, but a `c code` would be cool ... – edgarmtze May 17 '11 at 22:09
  • Are you expecting the resulting code to give the exact n-th Fibonacci number, or an approximation? If it's the latter, then in my view there's not a lot to discuss. The former, on the other hand, is more interesting. – NPE May 17 '11 at 22:10
  • 1
    Seeing as you now narrowed this down to a specific language, I can live with that. – Jeff Mercado May 17 '11 at 22:16
  • 2
    I've just verified, and a direct implementation using doubles gives correct (exact) answers to n=70. I'm voting to re-open hoping to see a discussion around whether Binet's formula can be used to get exact answers for n>70. – NPE May 17 '11 at 22:19
  • Amusingly, I recently wrote http://stackoverflow.com/questions/4327846/calculating-fibonacci/4328194#4328194 one (horrible) way of re-jiggering Binet's formula into giving exact answers. Theoretically. In practice you'll get bored past 20 ;-) – ephemient May 18 '11 at 01:00

2 Answers2

50

Here's a straightforward translation of the formula to Haskell: 

fib n = round $ (phi^n - (1 - phi)^n) / sqrt 5
  where phi = (1 + sqrt 5) / 2

This gives correct values only up to n = 75, because it uses Double precision floating-point arithmetic.

However, we can avoid floating-point arithmetic by working with numbers of the form a + b * sqrt 5! Let's make a data type for them:

data Ext = Ext !Integer !Integer
    deriving (Eq, Show)

instance Num Ext where
    fromInteger a = Ext a 0
    negate (Ext a b) = Ext (-a) (-b)
    (Ext a b) + (Ext c d) = Ext (a+c) (b+d)
    (Ext a b) * (Ext c d) = Ext (a*c + 5*b*d) (a*d + b*c) -- easy to work out on paper
    -- remaining instance methods are not needed

We get exponentiation for free since it is implemented in terms of the Num methods. Now, we have to rearrange the formula slightly to use this.

fib n = divide $ twoPhi^n - (2-twoPhi)^n
  where twoPhi = Ext 1 1
        divide (Ext 0 b) = b `div` 2^n -- effectively divides by 2^n * sqrt 5

This gives an exact answer.

Daniel Fischer points out that we can use the formula phi^n = fib(n-1) + fib(n)*phi and work with numbers of the form a + b * phi (i.e. ℤ[φ]). This avoids the clumsy division step, and uses only one exponentiation. This gives a much nicer implementation:

data ZPhi = ZPhi !Integer !Integer
  deriving (Eq, Show)

instance Num ZPhi where
  fromInteger n = ZPhi n 0
  negate (ZPhi a b) = ZPhi (-a) (-b)
  (ZPhi a b) + (ZPhi c d) = ZPhi (a+c) (b+d)
  (ZPhi a b) * (ZPhi c d) = ZPhi (a*c+b*d) (a*d+b*c+b*d)

fib n = let ZPhi _ x = phi^n in x
  where phi = ZPhi 0 1
hammar
  • 138,522
  • 17
  • 304
  • 385
  • Is this the same trick used [here](http://stackoverflow.com/questions/4327846/calculating-fibonacci/4328194#4328194) ? – Brian May 18 '11 at 13:59
  • @Brian: It's quite similar. The first one is based on [the matrix formula](http://en.wikipedia.org/wiki/Fibonacci_number#Matrix_form), but uses the same trick to get exponentiation for free. The second one does indeed use the same approach as I did, but it also uses [the binomial expansion](http://en.wikipedia.org/wiki/Binomial_expansion) to calculate only the required terms of (a + b√5)^n. – hammar May 18 '11 at 16:00
  • 1
    As a side note: this comes pre-defined in [Math.Algebra.Field.Extension](http://hackage.haskell.org/packages/archive/HaskellForMaths/0.1.9/doc/html/Math-Algebra-Field-Extension.html) (most of which I don't understand, but I used it successfully for this problem) – sleepyMonad May 18 '11 at 18:46
  • @Brian: Interestingly, it turns out that his second approach is much slower. It starts to choke at `n = 20`, while mine is virtually instant up to `n = 10^6` – hammar May 18 '11 at 19:07
  • @sleepyMonad: Nice! Although that's ℚ[√5], it might result in smaller intermediate numbers. I should try it and compare :) – hammar May 18 '11 at 20:35
  • @hammar: I have to say, this is **awesome**. I have [some questions](http://stackoverflow.com/questions/6082640) that pertains to your answer so maybe you might be able to answer them. – Jeff Mercado May 21 '11 at 15:29
  • 2
    What you need for this is `Z[(1+sqrt 5)/2]`, the ring of algebraic integers in `Q[sqrt 5]`. Nice is, with `phi = (1+sqrt 5)/2`, we have `phi^n = fib(n-1) + fib(n)*phi`. – Daniel Fischer Nov 20 '11 at 01:24
  • @DanielFischer: Ah, yes. That avoids the silly division which was annoying me. Also, that gives a nice way of calculating subsequent Fibonacci numbers since you get both `fib(n-1)` and `fib(n)` from it. – hammar Nov 20 '11 at 01:55
  • @hammar can you please elaborate more on the `ZPhi` multiplication case? I do not understand why is it any different than with the `Ext` case. – haskelline Jan 29 '14 at 22:35
  • @haskelline: The core of the difference is that while `(√5)² = 5`, an integer, instead `φ² = 1 + φ`, which has both an integral and a `φ` component. More lengthily: we have `(a₁ + b₁φ)(a₂ + b₂φ) = a₁a₂ + (a₁b₂ + a₂b₁)φ + b₁b₂φ²` by distributivity. Since `φ² = 1 + φ`, the last term is `b₁b₂(1+φ)`, and so we can combine coefficients to get `(a₁a₂ + b₁b₂) + (a₁b₂ + a₂b₁ + b₁b₂)φ`. In the `Ext` case, we instead have `(a₁ + b₁√5)(a₂ + b₂√5) = a₁a₂ + (a₁b₂ + a₂b₁)√5 + b₁b₂(√5)²`, and `(√5)²` is just `5`, so combining coefficients gets us `(a₁a₂ + 5b₁b₂) + (a₁b₂ + a₂b₁)√5`. Does this help? – Antal Spector-Zabusky Jun 29 '14 at 05:04
16

Trivially, Binet's formula, from the Haskell wiki page is given in Haskell as:

fib n = round $ phi ^ n / sq5
  where
    sq5 = sqrt 5 
    phi = (1 + sq5) / 2

Which includes sharing of the result of the square root. For example:

*Main> fib 1000
4346655768693891486263750038675
5014010958388901725051132915256
4761122929200525397202952340604
5745805780073202508613097599871
6977051839168242483814062805283
3118210513272735180508820756626
59534523370463746326528

For arbitrary integers, you'll need to be a bit more careful about the conversion to floating point values. Note that Binet's value differs from the recursive formula by quite a bit at this point:

*Main> let fibs = 0 : 1 : zipWith (+) fibs (tail fibs) 
*Main> fibs !!   1000
4346655768693745643568852767504
0625802564660517371780402481729
0895365554179490518904038798400
7925516929592259308032263477520
9689623239873322471161642996440
9065331879382989696499285160037
04476137795166849228875

You may need more precision :-)

Don Stewart
  • 137,316
  • 36
  • 365
  • 468