求斐波那契 (Fibonacci) 数列第 n 项的算法

最近看算法导论的习题时联想到 SICP 的一条习题后想到的，折腾了下，得到了几种算法，不同的算法时间复杂度有 Θ(aⁿ), Θ(n), Θ(lg(n)) 几种。时间复杂度为 Θ(lg(n)) 的算法比较有意思。下面一个个说，具体算法用 Common Lisp (CL) 实现，至于为什么用 CL 我在文章最后说明。

约定

符号 ^ 表示乘方运算。F(n) 表示 Fibonacci 数列的第 n 项，n 从 0 开始。由此， Fibonacci 数列可以由如下定义给出：

F(0) = 0
F(1) = 1
F(n) = F(n - 1) + F(n - 2)

这个定义本身是递归的。要是可以用在 Wordpress 里面使用 LaTeX 输入数学公式的插件就好了，下文中还会出现 _ 来表示下标，这些习惯跟 LaTeX 中数学公式的输入一样，不过看起来就会比较丑了。

最简单的递归定义和迭代实现

递归算法

根据定义可以很方便的得到一个求 F(n) 的算法：

(defun fibonacci-rec-simple (n)
  (cond ((= n 0) 0)
  ((= n 1) 1)
  (t (+ (fibonacci-rec-simple (- n 1))
              (fibonacci-rec-simple (- n 2))))))

defun 定义一个函数，参数是 n，cond 类似于 C 中的 switch，t 表示真，函数体的最后一个 form 作为返回值返回。另外注意数学表达式都是前缀表达式的形式。就算不懂 CL 应该也能看懂吧。

这个递归的版本应该是最容易写的了，直接翻译定义就行。不过可惜的是这个算法因为太多的重复计算导致时间复杂度非常大，计算 F(n) 时需要计算 F(0), F(1) 共 F(n + 1) 次（可以用数学归纳法很容易的证明），而 F(n) 是最接近 ((1 + sqrt(5))/2)ⁿ/sqrt(5) 的整数，所以这时计算 F(n) 的时间复杂度为 Θ(aⁿ)。在我的笔记本上，当 n = 50 时，我就已经不愿意等待到结果输出了。

迭代算法

如果定义变换 T(a, b) 为：

T(a, b) = (a + b, a)

那么，从 (1, 0) 开始，进行 n 次变换后，序对中的第二个元素就是 F(n) 了。由此可以得到求 F(n) 的迭代算法：

(defun fibonacci-iter-simple (n)
  (loop repeat (+ n 1)
     for b = 0 then a
     and a = 1 then (+ a b)
     finally (return b)))

这个算法的时间复杂度为 Θ(n)，没什么好说的。稍微解释下语法吧，loop 是 CL 中构造循环的一个宏，原来我是用 do 来进行循环的，但发现用 loop 更容易让人看懂，所以这里就写了个 loop 的版本。

这一段不少内容抄袭了 SCIP 1.2.2 中的分析，尽情鄙视我吧。

对递归定义的变形提高算法性能

抱歉，上面的内容没有什么新意，非原创。下面的应该会有点新意了。

最初的递归定义每次递归都只是将问题规模从 n 变到 n – 1 和 n – 2，我的想法是把递归定义进行变换，使得每次递归能够将问题规模一下子减小许多。因此我不断将 F(n) 进行递归展开，

F(n) = 1 * F(n - 1) + 1 * F(n - 2)
     = 2 * F(n - 2) + 1 * F(n - 3)
     = 3 * F(n - 3) + 2 * F(n - 4)
     = 5 * F(n - 4) + 3 * F(n - 5)
     = 8 * F(n - 5) + 5 * F(n - 6)
     = ...
     = F(p + 1) * F(n - p) + F(p) * F(n - (p + 1))

上面的关系式可以很容易使用数学归纳法证明。很自然的想到希望能够在每一次迭代的时候将问题的规模减小到原来的一半。

当 n 为奇数时，令 p = (n – 1) / 2，则：

F(n) = F((n + 1)/2)^2 + F((n - 1)/2)^2

当 n 为偶数时，令 p = n / 2，则

F(n) = F(n/2 + 1)*F(n/2) + F(n/2)*F(n/2 - 1)
     = F(n/2)^2 + 2*F(n/2 - 1)*F(n/2)

上面两个式子合起来可以写出下面的形式：

F(2*n) = F(n)^2 + 2*F(n)*F(n - 1)
F(2*n + 1) = F(n + 1)^2 + F(n)^2

有了上面的递归定义，可以得到下面的算法：

(defun fibonacci-rec (n)
    (cond ((= n 0) 0)
    ((= n 1) 1)
    ((evenp n)
           (let ((tmp (fibonacci-rec (/ n 2))))
             (* tmp (+ tmp (* 2 (fibonacci-rec (- (/ n 2) 1)))))))
    (t (+ (expt (fibonacci-rec (/ (+ n 1) 2)) 2)
      (expt (fibonacci-rec (/ (- n 1) 2)) 2)))))

这个算法在每次递归调用时都会将问题划分为两个规模为原先一半的子问题，表示该算法运行时间的递归式 (recurence) 为：

T(n) = 2*T(n/2) + c

其中 c 为进行乘法运算和加法运算需要的时间，将其看作常量。而 n/2 代表 (n / 2) 得到的数的上或下取整得到的整数。由主方法 (master method) 和替换法进行分析都得到这个算法的时间复杂度应该为 Θ(n)，从算法导论上的例子看这个递归式的时间复杂度也确实是 Θ(n)，但是从实际的运行情况来看其效率比上面的迭代算法快不少，使用的内存也少许多，与下面的迭代版本的 Θ(lg(n)) 的算法差距不大。（性能比较是在 SBCL 中使用 time 进行的，通过求 F(50000) 进行，这个数值在我的机器上除了 fibonacci-rec-simple 以外不会进行垃圾回收，因此避免了垃圾回收导致的性能下降。）是不是我的递归式弄错了，请知道的人请告诉我。

避免重复的计算

上面的递归算法中其实还是进行了不少的重复计算，而它的复杂度也因此不能降到 Θ (lg(n))。下面分两种情况考虑上面的递归算法。

如果递归过程中问题划分为求 F(2n) 和 F(2n + 1)，F(2n) 会通过 F(n), F(n – 1) 求得，而 F(2n + 1) 通过 F(n), F(n + 1) 求得。这里的关系如下：

F(2n)     --> F(n - 1), F(n)
F(2n + 1) --> F(n),     F(n + 1)

而如果递归过程中问题划分为求 F(2n – 1) 和 F(2n) 时，有如下关系：

F(2n - 1) --> F(n - 1), F(n)
F(2n)     --> F(n - 1), F(n)

上面的关系式中，F(n) 和 F(n – 1) 都重复出现，但实际上上面的算法都对它们进行了重复计算，减少这些重复计算就可以提高算法的性能。

上面的第二中情况需要计算两个子问题，而第一种情况里面要计算三个子问题，但实际上可以利用 F(n – 1), F(n) 来得到 F(n + 1)，因此实际上也只需要计算F(n – 1), F(n) 两个子问题就可以了。由此我们得到下面的关系式：

F(2n - 1), F(2n) --> F(n - 1), F(n)
F(2n), F(2n + 1) --> F(n - 1), F(n)

注意这里的关系式与前面式子的区别，前面的式子中，每个问题划分为两个子问题，而这里我们将求解两个 Fibonacci 数列中的项看成一个问题，每个问题被划分为一个更小的子问题。利用上面的关系式就可以避免算法执行中的重复计算。由此得到下面的程序：

(defun fibonacci-rec-2 (n)
  (cond ((= n 0) 0)
  (t (multiple-value-bind (x) (fib-two n) x))))

;; Returns the n'th and n-1'th fibonacci number.
(defun fib-two (n)
  (cond ((= n 1) (values 1 0))
  ((evenp n)
   (multiple-value-bind (i i-1) (fib-two (/ n 2))
     (values (* i (+ i (* 2 i-1)))
         (+ (* i-1 i-1) (* i i)))))
  (t
   (multiple-value-bind (i i-1) (fib-two (/ (- n 1) 2))
     (values (+ (* i i) (expt (+ i-1 i) 2))
         (* i (+ i (* 2 i-1))))))))

fib-two 返回的是 F(n) 和 F(n – 1)，利用了 CL 中函数返回多个返回值的特性。这个算法的递归式如下：

T(n) = T(n/2) + c

c 代表加法乘法运算的时间，n/2 的含义与上面的相同。这个递归式的运行时间为 Θ(lg(n))。实际测试的结果，这个算法的性能是最好的，内存开销也最小（为什么内存开销小不明白，可能跟我使用的 Common Lisp 实现有关），比下面要介绍的复杂度为 Θ(lg(n)) 迭代算法更快。

对迭代的变形提高算法性能

上面的递归算法都是从改变递归关系出发的，通过修改迭代过程也可以得到性能更好的算法。这里的思路是从 SICP 的习题 1.19 中看来的。有兴趣的直接去看原书好了。

考虑之前提到的 T 变换（不是递归式里面的 T），实际上求 F(n) 就是对最初的数值做 n 次 T 变换。

F(n) = T(T(T(...T(1, 0)))) = T^n(1, 0)

改变迭代过程的思路其实使用下面的关系式求乘方是一样的：

a^n = (a^2)^{n/2} = (a^4)^{n/4} = ... = (a^n)^1

具体算法实现的时候，当 n 不是偶数的时老老实实多乘一个 a，是偶数的话就平方。对 Fibonacci 数列，我们希望能够找到一个变换序列 {T_i}（下划线表示下标），使得：

T_1(a, b) = T^2(a, b)
T_2(a, b) = T^4(a, b)
T_3(a, b) = T^8(a, b)
...
T_{k-1}(a, b) = T^{n/2}(a, b)
T_k(a, b) = T^n(a, b)

对任意的 i，下式成立：
T_{i+1}(a, b) = T_i^2(a, b)

SICP 上面用的办法是定义 T{pq}(a, b) = (bp + aq + ap, bp + aq)，然后寻找变换 T{p’q’}，使得：

T_{p'q'}(a, b) = T_{pq}^2(a, b)

通过计算可以得到

p' = p^2 + q^2
q' = q^2 + 2pq

这个式子与之前改变递归定义得到的关系式形式完全相同，不过我没有看出来它们是如何联系起来的。最后给出具体的算法实现，我把书上 Scheme 的实现转换成 CL 实现。

(defun fibonacci-iter (n)
   (fib-iter 1 0 0 1 n))

(defun fib-iter (a b p q count)
  (cond ((= count 0) b)
  ((evenp count)
   (fib-iter a b
         (+ (* p p) (* q q)) ; p
         (* q (+ q p p))     ; q
         (/ count 2)))
  (t (fib-iter (+ (* b q) (* a q) (* a p)) ; a
           (+ (* b p) (* a q))         ; b
           p q (- count 1)))))

函数 fib-iter 虽然递归调用了自己，但其实是尾递归，实际的思想是迭代。我另外写了一个没有递归调用的版本，但是用了多个 setf 赋值，不够优雅，所以还是给原来的版本。按理说这个算法也是 Θ(lg(n)) 的时间复杂度，但是实际运行中，在求 F(50000) 时仅比 fibonacci-rec 快了没多少，比 fibonacci-rec-2 慢了不少。比后者慢可能是因为这个算法在输入为奇数时只能前进一小步，而 fibonacci-rec-2 不管输入是奇数还是偶数都可以将问题规模缩小一半。

关于递归中的重复计算

为了解决递归中的重复计算问题其实有另外的解决方案，将已经计算过的值存入某个缓存中。函数调用开始计算之前先到缓存中查，如果已经计算过则可以直接利用之前的结果，否则才进行计算。如果实现了这样的缓存，那么即使使用最初的 fibonacci-rec-simple，函数的时间复杂度也会降低到 Θ(lg(n))，而 fibonacci-rec 的时间复杂度也会降到 Θ (lg(n))。这个想法也是看 SICP 时看到的，在 1.2.2 的脚注 34 中提出，练习 3.27 给出了一个实现。

为什么用 Common Lisp

用 CL 不是为了装 B，我会的语言中用 C/C++, Java 来研究算法显然太麻烦，不支持大整数，性能测试不方便。而 Ruby 的实现会导致一个程序使用不同的语法而导致巨大的性能差异。比如 inject 加 block 的性能比 while 的性能差很多，但是前者在实现累加之类的操作时显然会更优雅。而我对 Scheme 的了解仅限于 SICP 上用到的那些。这样下来就只剩下 Common Lisp 了。

而事实上，用 CL 研究算法真的是非常方便。在 SLIME 和 Emacs 的支持下，我可以使用一个功能强大的 IDE，编译、调试、代码编辑都非常方便，更重要的是与系统交互式的开发过程。我可以在写完一个函数之后非常方便的测试这个函数，简化测试过程大大的减少了我实现一个算法的时间。另外语言本身的特性也使得算法的实现变得简单，比如让函数返回多个值，CL 本身就有特性支持这个功能。大整数使得我可以使用大的输入来测试性能，time 方便的让我得到函数运行时间、内存的统计信息。使用其他语言我很难享受到这些好处。

CL 也有不爽的地方：数学式子得转换成前缀表达式，不过这个也可以通过 macro 来解决，Google 到了几个现成的，暂时还没有研究怎么用。

Random Tech Thoughts

The title above is not random

求斐波那契 (Fibonacci) 数列第 N 项的算法

约定