Haskell中的斐波那契闭式表达式

Question

在 Haskell 中，斐波那契数列的闭合形式代码会是什么样子呢？

- edgarmtze

---

选择一种语言吧，你真的期望一个答案为所有这些语言提供代码吗？ - ildjarn

你可以自己写，表达方式已经给出了...你真的需要别人为你工作吗？ - erjiang

公式涉及指数、减法和除法。您能否更具体地说明您在哪一部分有困难？ - Rob Kennedy

既然你现在把它缩小到了特定的语言，我可以接受这个。 - Jeff Mercado

我刚刚验证了一下，使用双精度浮点数的直接实现可以得到n=70的正确（精确）答案。我投票支持重新开放讨论，希望能够讨论一下是否可以使用Binet公式来获得n>70的精确答案。 - NPE

Answer 1

这里是将公式直接翻译成 Haskell 的简单方法：

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

这段代码仅在n = 75时提供正确的值，因为它使用Double精度浮点算术。

然而，我们可以通过处理形式为a + b * sqrt 5的数字来避免浮点算术！让我们为它们创建一个数据类型：

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

由于指数运算是以Num方法的形式实现的，因此我们可以免费使用它。现在，我们需要稍微重新排列一下公式才能使用它。

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

这给出了一个精确的答案。

---

Daniel Fischer指出，我们可以使用公式 phi^n = fib(n-1) + fib(n)*phi 并使用形如 a + b * phi (即 ℤ[φ]) 的数字进行计算。这避免了笨重的除法步骤，并且只使用了一个指数运算。这样做可以得到一个更好的实现：

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

---

顺便提一下：这在Math.Algebra.Field.Extension中预定义（其中大部分我不理解，但我成功地用它解决了这个问题）。 - sleepyMonad

@sleepyMonad：不错！虽然那是ℚ[√5]，但可能会导致更小的中间数字。我应该试一下并比较一下 :) - hammar

你需要的是Z[(1+sqrt 5)/2]，即在Q[sqrt 5]中的代数整数环。很好的是，对于phi = (1+sqrt 5)/2，我们有phi^n = fib(n-1) + fib(n)*phi。 - Daniel Fischer

@hammar，您能否详细解释一下“ZPhi”乘法案例？我不明白为什么它与“Ext”案例有所不同。 - haskelline

@haskelline：区别的核心在于，虽然 (√5)² = 5 是一个整数，但是 φ² = 1 + φ 则同时具有整数和 φ 组成部分。更详细地说，我们通过分配律得到 (a₁ + b₁φ)(a₂ + b₂φ) = a₁a₂ + (a₁b₂ + a₂b₁)φ + b₁b₂φ²。由于 φ² = 1 + φ，最后一项为 b₁b₂(1+φ)，因此我们可以合并系数得到 (a₁a₂ + b₁b₂) + (a₁b₂ + a₂b₁ + b₁b₂)φ。在 Ext 的情况下，我们有 (a₁ + b₁√5)(a₂ + b₂√5) = a₁a₂ + (a₁b₂ + a₂b₁)√5 + b₁b₂(√5)²，而 (√5)² 就是 5，因此合并系数得到 (a₁a₂ + 5b₁b₂) + (a₁b₂ + a₂b₁)√5。这有帮助吗？ - Antal Spector-Zabusky

Answer 2

从Haskell wiki页面得知，Binet公式可以以Haskell语言表示为：

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

这包括分享平方根的结果。例如：

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

对于任意整数，转换为浮点数需要更加小心。请注意，此时Binet的值与递归公式相差很大：相当大。

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

你可能需要更高的精度 :-)

- Don Stewart