先记录一段流水帐，要不然会不记得自己是怎么找到这些东西的，没兴趣看流水帐的跳过。

我想看看用 Common Lisp 怎么写插入法排序（我对语言本身还不熟悉，所以想看别人怎么写的），于是 Google “common lisp insertion sort”，来到了这里，LiteratePrograms，一个记录了很多程序的 Wiki，很不错的地方。看了一会儿就想也许 Fibonacci 数列的算法也会有，于是 Google “common lisp fibonacci”，这次同样也有 LiteratePrograms 上面的文章出现在结果中，不过位置比较靠后，而且上面的算法也很普通。第三个结果来自 CLiki，文章在这里，这里的算法除了那个 generator 以外我在上一篇文章里面都有给出，包括那两个 Θ(lg(n)) 的算法。而最后在 LiteratePrograms 上逛的时候又看到了一个用矩阵乘法的算法。

利用 generator 来求 Fibonacci 数列

利用 generator 的算法比较特别，有点像 continuation。看过 Wikipedia 的解释以后才知道什么是 generator，而 continuation 可以用来实现 generator。代码如下：

(let ((a 0)
      (b 1))
  (defun fib ()
    (prog1 a (psetf a b b (+ a b)))))

这里的实现使用了 closure，let 返回了一个函数 fib，fib 中使用到的 a, b 在 let 创建的环境中，每一次调用 fib 都会返回 Fibonacci 数列中的下一个值，并且更新环境中的 a, b。prog1 对第一个 form 求值，然后对后面的 form 进行求值，最后返回第一个 form 的值。而 psetf 在赋值之前对所有的需要求值的 form 进行求值，然后在赋值，这样就避免了要使用一个变量来记录 a 的值的需要。

这个算法只是用函数的形式表达迭代而已，除此以外没有什么特别的。但是注意使用 generator 我们相当于得到了一个长度无限的 Fibonacci 数列，我们方便地可以控制迭代的过程，在需要的时候才进行迭代。

关于 closure 再多说两句。这个实现里的 fib 函数实际上带有自己的状态，其状态就是用 let 定义时候的 a,b 变量，注意这两个变量在函数之外是无法访问到的。使用 C++ 当然也可以实现相同的效果，但是你得定义一个类，里面封装这两个状态变量，然后重载括号符，使这个类成为一个 functor，实在是麻烦，有兴趣的可以去 Wikipedia 的 generator 条目里面看看 C++ 的实现。而使用 closure 的话，创建带状态的函数一切都是那么自然，不知道 Java 7 里面最终是否会加入 closure。

Dijkstra 的笔记，In honour of Fibonacci

在 CLiki 上发现我所使用的 Θ(lg(n)) 的算法原来 Dijkstra 也有给出过，在他的一篇名为 In honour of Fibonacci 的笔记中记录着。去看了看他的笔记，原来我用的算法和 Dijkstra 记录下来的一样，不过他在计数的时候从 1 开始而不是从 0 开始。（Dijkstra 是不是也是数学家？）Dijkstra 还把这个关系式泛化了，不过有几处推导过程没有看懂，也不知道泛化以后有什么用。

用矩阵乘法求 Fibonacci 数列

另外在 LiteratePrograms 上看到求 Fibonacci 数列的 Ruby 版本中，有一个使用矩阵的算法。方法是对矩阵 ((1 1) (1 0)) 做 n – 1 次乘方，即可得到 Fibonacci 数列第 n 项。代码如下：

require 'matrix'

FIB_MATRIX = Matrix[[1,1],[1,0]]
def fib(n)
    (FIB_MATRIX**(n-1))[0,0]
end

这个算法的优化就可以交给矩阵乘法去做了，矩阵乘法如果是 Θ(lg(n)) 的时间复杂度那么这个算法也是 Θ(lg(n)) 的时间复杂度了。

后来又到 Wikipedia 上搜了一下，在 Fibonacci number 条目下有关于 Fibonacci 数列的非常之多的数学关系和等式的证明，还列举了非常多的参考资料，包括 Dijkstra 的笔记。矩阵形式在上面也有，而且通过矩阵乘法的性质很容易的推导出了复杂度为 Θ(lg(n)) 算法中利用的等式。众人编辑的力量果然强。

为什么不直接用通式求 Fibonacci 数列的第 n 项

Fibonacci 数列的通式我当然知道，但是直接使用它来求数列的第 n 项在精度上会出现问题。通式中有无理数，因此计算过程中会产生误差。根据 CLiki 上的说明，在 CLISP 2.33.2，32 bit 的机器上，对单精度浮点数，当 n < 32 时，使用 truncate 可以得到正确结果；对双精度浮点数，n < 76 时，使用 rounding 可以得到正确结果。超出这个范围都将在整数级别上出现错误。

除了计算上的问题，我觉得这个问题本身也比较有意思。就像 MIT 的教授在上算法导论时说的，”Speed is fun!“，我也这么觉得。而且，当发现自己做着 Dijkstra 曾经也做过的事情的时候还是会小有成就感的。