Fibonacci 数列及其计算方法

斐波那契数列（Fibonacci sequence），又称黄金分割数列，这个数列最早是由印度数学家提出来的。不过更多的人学习到这个数列是因为意大利数学家列昂纳多·斐波那契（Leonardoda Fibonacci）和他的《Liber Abaci》一书。在这本书中，列昂纳多·斐波那契以兔子繁殖为例子引出了这个序列，因此这个序列又称为“兔子数列”。
这个序列的前几项是这样的：0,1,1,2,3,5,8,13,21,34,⋯
在数学上，斐波纳契数列以如下被以递归的方法定义：

Fibonacci 数列的通项公式

Fibonacci 数列除了递归形式之外，当然还可以写出通项公式。下面就来算算 F(n)。

是典型的线性差分方程，可以用经典的待定系数法来解，当然也可以用 Z 变换解。考虑到不是每个人都学过 Z 变换，这里就给个最基本的待定系数法。假设 F(n)=C×Wn，C 和 W 是两个待定系数，那么有：

化简一下，得到：

很显然，W 有两个解：

那么 F(n)=C1Wn1+C2Wn2 也满足 F(n)=F(n−1)+F(n−2)。
C1,C2 可以通过起始条件来确定。

这是一个简单的二元一次方程，计算后可以得到：

所以：

当 n 较大时，C1Wn1=−0.4472×(−0.618)n≈0
所以 n 较大时，

有些没有学过差分方程理论的同学可能会问为什么假设 F(n)=C×Wn。简单的说，这是猜的。当然也不是无缘无故的猜测。我们知道 F(n) 的差分方程是这样的：

我们再构造两个辅助的差分方程：

那么当起始条件相同时，明显有 F2(n)<F(n)<F1(n)。
F1(n),F2(n) 都是等比数列，通项公式如下：

所以 2√nC2<F(n)<2nC1。
所以我们猜测 F(n)=Wn，其中 2√<W<2。

程序实现

斐波纳契数列的定义是递归形式的。自然，用递归函数最容易实现。

uint64_t Fibonacci(unsigned char n)
{
if( n == 0 ) return 0;
if( n == 1 ) return 1;
return Fibonacci(n - 1) + Fibonacci(n - 2);
}

这个代码虽然简单，但是计算复杂度却太高。当 n 稍微大一点时，计算时间就会非常长。

我们可以简单的算一下，当计算 F(n) 时，F(n) 内部会调用 2+F(n−1)+F(n−2) 次自己。

我们用 B(n) 表示 F(n) 中多少次调用自己。那么就有递归表达式：

那么这个 B(n) 具体是多少呢？ 显而易见，当 n≥2 时，B(n)≥F(n)。

我们又知道 F(n)<B(n)<2F(n)，所以递归算法的时间复杂度是 O(1.618n)。1.618100=7.90409×1020。也就是说当 n=100 时，递归函数要反复调用自己 1020 次左右，这个计算时间是不能忍受的。

如果将递归算法改为递推，时间复杂度会降低很多。

uint64_t Fibonacci(unsigned char n)
{
if( n == 0 ) return 0;
if( n == 1 ) return 1;
uint64_t fn, fnn = 1, fnnn = 0;
for(int i = 2; i <= n; i ++)
{
fn = fnn + fnnn;
fnn = fn;
fnnn = fnn;
}
return fn;
}

这个程序很简答，for 循环中计算了n−2 次， 所以时间复杂度为 O(n)。

那么原来的递归算法就没有改进的余地了吗？递归算法之所以算的慢，是因为计算F(n) 时，对于小于 n 的数 m，重复计算了很多次F(m)。下面这幅图给出了计算 F(6) 时的递归调用关系。

可以看到 F(4),F(3) 都重复计算了好几遍。其实，只要将计算过的F(m) 保存下来，下次用时直接读取就好了，这样就省去了反复计算 F(m) 的问题。下面是改进后的代码，时间复杂度降低为 O(n)：

uint64_t Fibonacci(unsigned char n)
{
static uint64_t fib[256] = {0, 1};
if(n == 0) return 0;
if( fib[n] != 0 ) return fib[n];
fib[n] = Fibonacci(n - 1) + Fibonacci(n - 2);
return fib[n];
}

这个代码还用到了 static 型局部数据变量的一个特性：没有指定值的元素自动初始化为 0。如果没有这个特性，我们的程序中就要写上 255 个 0 了。

除了上面给出的两种 O(n) 算法之外，还有没有更快的计算 F(n) 的方法呢？ 答案是肯定的。
我们将递推表达式变变型就能得到：

因此，计算 F(n) 就变成了计算一个矩阵的 n−1 次方的问题了。而这个问题是有快速算法的。
当 n 为偶数时，Mn=Mn/2×Mn/2 所以只要计算出了 Mn/2 了之后再自乘一下就行了。
当 n 为 奇数时，Mn=M(n−1)/2×M(n−1)/2+M 所以只要计算出了 M(n−1)/2 了之后自乘一下在加一下自己就可以了。

也就是说当 n=2m 时，只需要进行 m 次乘法计算就可以了。当 n<m 时，最多也就是进行 m 次乘法计算 和 m 次加法计算。所以算法复杂度是 O(log(n))。当 n 较大时，log(n)≪n，所以当 n 较大时这样计算会快的多。这里就不写出这种算法的具体实现了，感兴趣的同学可以自己写写。

这里还有一个问题一直没有解决，就是我们的程序可以正确计算到 n 为多少。 我们知道一个 m 位的无符号整数可以表示的最大的整数是 2m−1。当 F(n)>2m−1 自然计算结果就是错误的了。 估计这个最大的 n 还要用到我们上面给出的通项公式的近似结果：

简单的计算可知：

所以：

所以 32 位无符号整数最大可以表示到 F(47)，64 位无符号整数最大可以表示到 F(93)。