Fibonacci's Closed-form expression in Haskell

How would the Fibonacci's closed form 开发者_如何学Ccode look like in haskell?

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

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