斐波拉契

斐波纳契数列的定义

斐波纳契数列以意大利数学家雷纳尔多·波纳契（也称为斐波纳契）的名字命名，可在多个学科中找到它的应用，比如数学、建筑学、几何学、计算机科学和自然中。它是从 0, 1 对（或 1, 1）开始的整数数列。数列中的后一个数是前两个数的和。用数学中的项来表示就是，斐波纳契数列 (Fn) 定义由与前两项（F0、F1 或 F1、F2）的递归关系来定义，如下所示：

1 2 3 4 5 6 7

F0 = 0 F1 = 1 Fn = F(n-2) + F(n-1)  if n >1 or F1 = 1 F2 = 1 Fn = F(n-2) + F(n-1)  if n >2

斐波纳契数的另一种增强形式是，负斐波纳契数列，其中的索引扩展到负数。在斐波纳契数的三种定义中，本文的分析基于该数列的更加现代的形式 F0, F1, …, Fn。表 1 中给出了前几个斐波纳契数。

表 1：前几个斐波纳契数的列表

F0	F1	F2	F3	F4	F5	F6	F7	F8	F9	F10	F11
0	1	1	2	3	5	8	13	21	34	55	89

递归算法

基于斐波纳契数列的递归定义，计算第 n 个斐波纳契数的递归算法是一种直观的实现。

清单 1：递归算法的伪代码

1 2 3 4 5 6

Fib(n): if n in range [0,1]     return n else     return Fib(n-2) + Fib(n-1) end if

递归算法与该数列的递归定义具有自然的对应关系。递归步骤的运行时复杂度可用以下等式来表达：

1 2 3

T(n) = T(n-1) + T(n-2)+ Theta(1)      >= 2 T(n-2)      = Theta(2<sup>n/2</sup>)

因此，此实现的复杂度是指数时间级的。

清单 2：计算第 n 个斐波纳契数的递归算法

1	long long fib_r(int n)<br>{<br> return n <= 1 ? n : fib_r(n-1) + fib_r(n-2) ; <br>}

图 1. 斐波纳契递归算法的递归树

从图 1 中显示的算法的递归树，可以看到一些表达式在过程中计算了多次。例如，在 Fib(n) 和 Fib(n-1) 执行的调用中，都计算了表达式 Fib(n-2)。反复计算是递归算法具有指数时间的主要原因。如果可以避免所有多余的计算，时间复杂度将获得改善。对于斐波纳契数列，这个问题可以使用动态编程解决。通过利用记忆 保留计算的值供以后使用，动态编程技术消除了在递归算法中发生的反复计算，加快了执行速度。

动态编程算法

为了消除多余的计算，动态编程利用了默记法 (Memoization) 来保留计算的值供以后使用。由于每一项只计算一次，所以动态编程方法的时间复杂度是线性的。本文将讨论计算斐波纳契数列的两种著名的动态编程算法：默记法版本和自底向上版本。

动态编程方法的默记法版本

从该技术的名称可以看出，该方法通过将计算的值存储在某种形式的数据结构中来记住它们。只有在新计算的结果不在算法的记忆中时，才需要计算它。尽管默记法方法仍依赖于递归步骤，但默记法查找可以确保每个表达式在该过程中只计算一次。记忆查找操作被认为具有不变的时间复杂度。

清单 3：默记法动态编程的伪代码

1 2 3 4 5 6 7 8 9 10 11 12

Mem = { }                                 //Data structure storing computed values Fib(n):     if n in Mem                       //Return existing value from Memory         return Mem[n]             //Early exit when handling known values     end if if n if n in range [0,1]                 //Base step of the recursion     fib ? n else                                    //Recursive steps     fib ?  Fib(n-2) + Fib(n-1)     end if     Mem[n] ? fib                   //Memorize new value     return fib                      //Return the new computed value

因为斐波纳契数的每个值只计算一次，所以此方法的运行时复杂度为线性时间。如图 2 所示，动态编程方法的默记法版本的递归树不会反复计算已知值。

图 2. 动态编程算法的递归树

尽管默记法版本可以降低运行时复杂度，但它仍依赖于递归函数调用。设置和管理函数调用的成本已在复杂度计算中忽略，但它们实际上并不是微不足道的。下面即将讨论的技术称为自底向上遍历，可以帮助避免支持函数调用的额外成本。

动态编程的反向遍历（自底向上）版本

反向遍历方法不从递归树顶部的根开始遍历，而是自底向上遍历递归树。在计算斐波纳契数的特殊情况下，该算法将从索引 0 开始，然后向上遍历到索引 n。它仍使用一种数据结构来记住计算的值，但它将递归函数调用替换为查找操作。根据定义，第 n 个斐波纳契数计算为 Fn = F(n-2) + F(n-1)。因为前面的两项 F(n-2) 和 F(n-1) 可从算法的记忆中进行检索，所以不需要递归函数调用。

清单 4：反向遍历动态编程的伪代码

1 2 3 4 5 6 7 8 9 10

Fib(n):     Mem = {}                                         //the memory     for k in range [0,1,…,n]                         //traverse bottom-up         if k in range [0,1]                         //compute each term             Mem[k] ? k                             //starting with indices 0 and 1         else:             Mem[k] ? Mem[k-1] + Mem[k-2]           //Retrieve known values         end if     end for     return Mem[n]

由于记忆查找可预防反复计算已知的项，所以斐波纳契数列的每一项都只计算了一次。因此该算法具有线性时间复杂度。此外，因为递归函数调用被替换成了记忆查找操作，所以消除了设置和管理递归函数调用的额外成本。

清单 5：计算第 n 个斐波纳契数的动态编程方法

1 2 3 4 5 6 7 8 9 10 11 12

long long fib_d(int n) {    int k;    long long memory[n];    for ( k = 0 ; k <= n ; ++k ) {       if ( k <= 1 )          memory[k] = k;       else          memory[k] = memory[k-1] + memory[k-2];    }    return memory[n]; }

迭代算法

类似于反向遍历方法，迭代算法迭代从索引 0 一直到索引 n 的整个数列。在此过程中，它使用了前两个数的和更新了斐波纳契数列的每一项。但是，迭代算法没有记住计算的值。图 3 提供了用于计算斐波纳契数的迭代算法的可视化表示。

图 3. 迭代算法的可视化

点击查看大图

因为为该程序设定的需求是计算第 n 个斐波纳契数，所以在实际实现中，可以跳过使用计算的值 res 更新斐波纳契数列的实际值的步骤。在需要时，很容易实现该算法。

清单 6：计算第 n 个斐波纳契数的迭代方法

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

long long fib_i(int n) {    int k;    long long e0 = 0, e1 = 1, res = 0;    for ( k = 0 ; k <= n ; ++k ) {       //iterate from 0 to n       if ( k <= 1 )          res = k;                      //base case for F0 and F1       else {          res = e0 + e1;               //res is the sum of 2 previous terms             e0 = e1;                     //update the 2 previous terms          e1 = res;                    // -------- ditto --------       }    }    return res;                      //return Fn }

因为迭代算法仅迭代整个斐波纳契数列一次，所以它的复杂度为线性时间级。类似于反向遍历算法，它不会调用递归函数。

优化的矩阵乘方算法

这很可能是具有最佳性能的算法：它对斐波纳契数的计算具有对数时间复杂度。该算法取决于以下数学运算。

图 4. 矩阵乘方算法的数学基础

第 n 个斐波纳契数将是将矩阵 A 扩展到 (n-1) 乘方后的结果矩阵的元素 (0, 0)。因为关注的所有矩阵都是二维的，所以矩阵运算可在固定的时间内完成。清单 7 中显示了一段示例代码。

清单 7：固定时间内的二维矩阵乘法

1 2 3 4 5 6 7 8 9 10 11 12

void multiply(long long M[2][2], long long N[2][2]) {   long long a =  M[0][0]*N[0][0] + M[0][1]*N[1][0];   long long b =  M[0][0]*N[0][1] + M[0][1]*N[1][1];   long long c =  M[1][0]*N[0][0] + M[1][1]*N[1][0];   long long d =  M[1][0]*N[0][1] + M[1][1]*N[1][1];     M[0][0] = a;   M[0][1] = b;   M[1][0] = c;   M[1][1] = d; }

为了将矩阵 A 扩展到 n 次方，一种直观的 算法将具有线性的运行时间，因为它必须调用乘法运算 (n-1) 次。

An = A A … A

这可以使用以下技术进一步优化：

An = A(n/2) A(n/2)

因此，无需将矩阵 A 扩展到 n 次方，该算法将它扩展到 (n/2) 次方并将 An/2 与自身固定次数的乘法操作。

清单 8：具有对数时间的矩阵乘方算法

1 2 3 4 5 6 7 8 9 10 11

void power(long long M[2][2], int n) {   if( n == 0 || n == 1) return;               //no operation needed for base case   long long A[2][2] = {{1,1},{1,0}};     power(M, n/2);                             //recursively call power on n/2   multiply(M, M);                           //a constant time multiplication     if (n%2 != 0)                             //To handle odd n      multiply(M, A); }

通过在 n/2 上调用乘方函数，输入的大小会每次减半。相应地，复杂度与对数时间有关。清单 9 显示了一种使用优化的矩阵乘方方法的示例实现。

清单 9：计算第 n 个斐波纳契数的优化的矩阵乘方方法

1 2 3 4 5 6 7 8

long long fib_m(int n) {   long long A[2][2] = {{1,1},{1,0}};              //matrix A   if (n == 0)       return 0;                                  //base case F0 = 0   power(A, n-1);                                 //raise A to power (n-1)   return A[0][0];                               //Fn is element (0,0) }

内联汇编实现

内联汇编实现将基于迭代算法。从图 3 中显示的可视化形式可以看到，迭代算法计算第 n 斐波纳契数的方法是将它之前的两项 (n-2)和 (n-1)相加。(n-2) 项在此计算后可丢弃，因为后续计算不会再使用它。此观察发现可用来实现内联汇编方法。

算法

内联汇编实现将两个最新的斐波纳契数保留在两个单独的寄存器中。在每次迭代时，该算法都会计算最新的项的和，然后将计算的值存储在保存较小的项的寄存器中。这将实际覆盖后一项。因为未来的计算不需要较小的项，所以删除此值是可接受的。因此较大的项是这次迭代的结果。图 4 显示了内联汇编方法的可视化表示。

图 5. 内联汇编方法的可视化表示

点击查看大图

图 4 中已经非常直观地显示，该实现在斐波纳契数列中的每一项上迭代一次。每次迭代的结果将存储在图 4 中红色的寄存器中。该实现的伪代码如清单 10 所示。

清单 10：内联汇编实现的伪代码

1 2 3 4 5 6 7 8 9

Fib(n):     load F0 to register 0                              //smaller term in register 0     load F1 to register 1                              //greater term in register 1     for k in range [2,3,…,n]         swap the values in two registers              //register 1 holds smaller term         sum up the values in two register             //sum up 2 terms         store the sum to register 1                   //register 1 holds the sum     end for     return register 1                                 //return the sum

通过三异或运算 ( triple exclusive-or operation) 交换两个值

可通过各种方式交换存储在两个寄存器中的值，一种是使用临时寄存器或变量的简单技术。但是，还有一种不需要经历临时存储而交换值的更快方式：使用三次连续的异或 (XOR) 运算的异或交换算法。表 2 中给出了三重 XOR 运算的事实表。

表 2：三异或运算的事实表

初始值		R0 = R0 XOR R1		R1 = R0 XOR R1		R0 = R0 XOR R1
R0	R1	R0	R1	R0	R1	R0	R1
0	0	0	0	0	0	0	0
0	1	1	1	1	0	1	0
1	0	1	0	1	1	0	1
1	1	0	1	0	1	1	1

该事实表显示，两个寄存器 R0 和 R1 中的值在执行三次连续 XOR 运算后已交换。

汇编指令

有许多汇编指令可用来实现该算法。有关可用指令的完整列表，请参阅图书 z/Architecture 操作原理（IBM 出版物编号 SA22-7832-10）。在 IBM z Systems™ 上的长长的汇编指令列表中，本文选择使用以下指令来实现斐波纳契数列计算。

XGR R0, R1：

• 存储在第一个和第二个寄存器中的值的异或结果放在第一个操作数位置 (R0)。
• 存储在第 2 个操作数 (R1) 中的值保持不变。
• 这些操作数和结果的长度为 64 位。

AGR R0, R1：

• 将第二个操作数与第一个操作数相加，将两数的和放在第一个操作数所在的位置上。
• 存储在第 2 个操作数 (R1) 中的值保持不变。
• 这些操作数和结果被视为 64 位有符号二进制整数。

BRCT R0，分支地址：

• 从第一个操作数中减去值 1，将结果放回到第一个操作数位置。
• 结果为 0 时，继续执行正常的指令数列。
• 结果不为 0 时，当前程序状态字 (PSW) 中的指令地址替换为分支地址。

内联汇编代码

计算第 n 个斐波纳契数的实际的内联汇编代码如清单 11 所示。

清单 11：计算第 n 个斐波纳契数的内联汇编方法

1 2 3 4 5 6 7 8 9 10 11 12 13 14

long long fib_asm(int upto) {     if ( upto < 2 ) return upto;                 //F0 = 0, F1 = 1     int mycount = upto - 1;     long long f0 = 0L, f1 = 1L;     asm ( "0:                 "                  //back to this when mycount > 0           "XGR   %0, %1       "                  //1st XOR           "XGR   %1, %0       "                  //2nd XOR           "XGR   %0, %1       "                  //3rd XOR           "AGR   %1, %0       "                  //Add 2 terms, store to f1           "BRCT  %2, 0b       "                  //back to 0 if mycount > 0          :"+r"(f0), "+r"(f1), "+r"(mycount)         );     return f1;                                     //return f1  }

在循环的开头，在值 f0 和 f1 中分别加载 0 和 1。变量 mycount 用于保存所需的斐波纳契数的索引；它也是要执行的迭代次数。在每次迭代的开头，f0 保存更小的项，f1 保存斐波纳契数列的更大的项。在每次迭代期间，该算法通过三次 XOR 运算来交换 f0 和 f1 的值。作为交换的结果，f0 保存更大的项，f1 保存更小的项。AGR 运算对 f0 和 f1 的值求和，然后将该和存储回 f1 中。在任何迭代结束时，f1 保存更大的项，f0 保存更小的项。当 mycount 大于 0 时，BRCT 运算从 mycount 减去 1 并将分支重置回标签 0。mycount 变为 0 时，汇编指令结束。程序会返回存储在 f1 中的值，这是第 n 个斐波纳契数。

因为内联汇编仅迭代斐波纳契数列的每一项一次，所以复杂度为线性时间级。

性能预期

在上面给出的 5 种实现中，具有对数时间复杂度的矩阵乘方算法在性能上排名第一。动态编程方法、迭代算法和内联汇编实现位居中间，具有线性时间复杂度。递归算法由于其指数时间复杂度而位于最底部。

表 3：运行时复杂度

方法	递归算法	迭代算法	动态编程	内联汇编	矩阵乘方
复杂度	指数	线性	线性	线性	对数

基于运行时分析，预计在采样大小足够大时，实际性能与复杂度呈正比。

测试程序

每个实现都可以封装在一个函数中。测试程序使用了以下两个变量：

volatile clock_t start、end；

它们记录每个函数的开始时间和返回时间。使用了 5 个不同的变量来保存每个实现的运行时间。

volatile float elapsedASM、elapsedI、elapsedR、elapsedM、elapsedD；

清单 12 给出了运行时间的计算方法。

清单 12：运行时间的计算

1 2 3 4

start = clock(); while ( myCount-- ) { resASM = fib_asm(limit); } end = clock(); elapsedASM =  ((float) (end - start))*1000 / CLOCKS_PER_SEC ;

除了递归算法之外，所有实现都将对 92 项执行 10,000,000 次计算。大量的循环次数用于延长执行时间，以便正确记录运行时间。另外，由于大于斐波纳契数列中第 100 个数的项超出了整数数据类型长度，所以大于第 90 个数的项实际上是具有 64 位整数的斐波纳契计算的上限。递归函数将使用小得多的参数来进行测试，它只计算斐波纳契数列的第 40 项两次。由于递归算法具有指数时间复杂度，所以应该选择较小的参数来确保程序在其他方法的时间范围内完成。

请参阅完整的 测试程序 了解更多细节。

实际性能

图 6 中给出了在 IBM 多伦多软件实验室的一台运行 Linux 的 z System 机器上执行测试的实际性能结果。

图 6. 相对性能

点击查看大图

从图 6 可以明显地看到，在计算最大的斐波纳契数时，内联汇编实现优于其他所有方法。它大约比动态编程和迭代方法快 4 倍。令人惊奇的是，它的性能也优于经过优化的矩阵乘方算法。这一性能优势源于内联汇编实现较小的资源占用。

从理论上讲，如果输入的大小足够大，优化的矩阵乘方算法的性能可能超过内联汇编实现。但实际上，当应用程序的输入受到变量大小限制时，从内联汇编实现生成的紧凑代码的实际速度可能超过某种更快的算法的性能。这是使用内联汇编进一步提高性能，超越最佳算法可以提供的性能的一个示例。

结束语

在计算斐波纳契数列的特殊情况下，内联汇编方法的性能优于其他所有算法，即使已证明理论上速度将会更慢。这个示例证明，理论性能可能不切实际。适当地使用内联汇编来调优最注重性能的代码节，这是超越算法所提供的默认速度的一条途径。

一定要注意的是，使用内联汇编所获得的性能只有通过仔细计划和大量测试才能实现。因为 IBM XL 编译器基于多维分析来执行高度复杂的优化，所以在大多数情况下，使用编译器所提供的适当的优化水平可能是最佳选择。必须根据特定系统上的经验测试数据来决定是否使用内联汇编。此外，内联汇编仅应用于最注重性能的代码节。成功实现的前提是：应用程序的固有知识和对系统上可用的汇编程序指令集的全面理解。copy的，好好复习一下。