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
与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…