斐波那契数列在诸多方面，从头状花序到通用编码（如禁忌编码）表现出种种有趣的性质。研究高效计算斐氏数列的过程充满乐趣。
斐波那契数列的经典公式表示为如下递推形式：
由此引出其直观（朴素）的递归实现：
1
2
3
4
5
6
7
| uint64_t fibo_rec(uint64_t n)
{
if (n <= 2)
return 1;
else
return fibo_rec(n - 1) + fibo_rec(n - 2);
}
|
看上去很完美，可是计算的过程中产生了次递归调用！运行时间呈指数增长（因为，其中为黄金比例）。
通过显式地优化尾递可以消除递归调用：
1
2
3
4
5
6
7
8
9
10
11
12
13
| uint64_t fibo_iter(uint64_t n)
{
uint64_t n_1 = 1, n_2 = 0;
for (uint64_t i = 0; i < n; i++)
{
uint64_t t = n_1 + n_2;
n_2 = n_1;
n_1 = t;
}
return n_1;
}
|
这次计算的时间复杂度下降到，大大优于原始递归版本。迭代版本借助临时变量保存和，并使用移位寄存器将后继的斐波那契数列分别保存到n_2和n_1。有些人可能会忌讳临时变量，那么可以将上述代码重写为：
1
2
3
4
5
6
7
8
9
10
11
12
| uint64_t fibo_iter_var(uint64_t n)
{
uint64_t n_1 = 1, n_2 = 0;
for (uint64_t i = 0; i < n; i++)
{
n_1 = n_1 + n_2;
n_2 = n_1 - n_2;
}
return n_1;
}
|
n_2 = n_1 - n_2部分化简得n_1（，即，结果为）。此处我们通过额外操作抵消掉临时变量，好坏别当别论。
⁂
是否有进一步优化的空间？根据以往的经验，线性时间似乎是最好的结果。无独有偶，有些聪明的家伙发现矩阵
具有“平移”并增加一个矢量分量的性质。比如
看上去有几分眼熟吧。事实上我们将a代换为，b代换为，得
有点眉目了！但这是基于已知和的前提下，一旦成立意味着存在（递归）分解形式：
反复如此迭代到初始条件和，那么有
矩阵运算似乎代价不菲，至少要进行线性次数的矩阵乘法。幸运的是不需要那么多次，我们有办法在时间复杂度内计算出（此处矩阵维数为常数，故可视矩阵乘积为常量）。如何计算呢？且观察式子和。整数情况时计算过程如下：
1
2
3
4
5
6
7
8
9
10
11
12
13
| int expo(int x, int n)
{
int t = 1;
int y = x;
while (n)
{
if (n & 1)
t *= y;
y *= y;
n >>= 1;
}
return t;
}
|
矩阵计算方式同理可得。以上代码中t为乘积当前值，y为乘积“平方”。给定两个矩阵
和
第一个初始化为单位阵，第二个初始化为斐波那契矩阵。根据矩阵的对称性，和各自恒等于和，这两个变量可以分别消掉。
和
其中虚线位置表示对称部分。将上述式子代入幂函数，将有
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
| uint64_t fibo_expo(uint64_t n)
{
uint64_t a_t = 1, b_t = 0, d_t = 1, // "1"
a_y = 0, b_y = 1, d_y = 1; // "y"
while (n)
{
if (n & 1)
{
// t*=y;
uint64_t a = a_t * a_y + b_t * b_y;
uint64_t b = a_t * b_y + b_t * d_y;
uint64_t d = b_t * b_y + d_t * d_y;
a_t = a;
b_t = b;
d_t = d;
}
//y*=y;
uint64_t a = a_y * a_y + b_y * b_y;
uint64_t b = b_y * (a_y + d_y);
uint64_t d = b_y * b_y + d_y * d_y;
a_y = a;
b_y = b;
d_y = d;
n >>= 1;
}
return b_t;
}
|
提取出公共子式，改变赋值顺序，消除多余的临时变量，得
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
| uint64_t fibo_expo_var(uint64_t n)
{
uint64_t a_t = 1, b_t = 0, d_t = 1, // "1"
a_y = 0, b_y = 1, d_y = 1; // "y"
while (n)
{
if (n & 1)
{
// t*=y;
uint64_t bx = b_t * b_y;
b_t = a_t * b_y + b_t * d_y;
a_t = a_t * a_y + bx;
d_t = bx + d_t * d_y;
}
//y*=y;
uint64_t b2 = b_y * b_y;
b_y *= (a_y + d_y);
a_y = a_y * a_y + b2;
d_y = b2 + d_y * d_y;
n >>= 1;
}
return b_t;
}
|
太棒了，我们找着了一个的算法来计算斐波那契数列。不过仍有稍加改进的余地。再怎么改进呢？且注意，
这样一来消去变量d，得到仅包含变量a和b的方程。修改乘积项和平方项使其只包含变量a和b：
以及
其中虚线表示无需关心的部分。利用上述结果，可得：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
| uint64_t fibo_expo_var_2(uint64_t n)
{
uint64_t a_t = 1, b_t = 0, // "1"
a_y = 0, b_y = 1; // "y"
while (n)
{
if (n & 1)
{
// t*=y;
uint64_t bx = b_t * b_y;
b_t = a_t * b_y + b_t * a_y + bx;
a_t = a_t * a_y + bx;
}
//y*=y;
uint64_t b2 = b_y * b_y;
b_y = 2 * a_y * b_y + b2;
a_y = a_y * a_y + b2;
n >>= 1;
}
return b_t;
}
|
这几个算法就运行速度比较结果如何？原始递归算法效率低得惊人，生成需要一分多钟，而其他算法只花费几微秒。进行10000000（1000万）次迭代测试（使用相同的随机输入，随机数n位于1到70之间），结果如下
| iter |
3.29529s |
| iter_var |
3.30153s |
| expo |
2.28118s |
| expo_var |
2.2531s |
| expo_var_2 |
2.22531s |
最后三个版本的确快很多，但是快速乘幂算法版本之间的效率差距并不显著。假如改进程度相比迭代版本没有那么大的话，我们应该意识到刚只是计算相对较小的斐波那契数列。64位无符号整数最大仅允许存储（因为，，所以两边取对数解出n，推出），不过通过大整数（抽象为类）来计算，远比迭代方法和递归算法快得多得多。