首页 > 代码库 > 菲波那契数列的快速幂矩阵求法

菲波那契数列的快速幂矩阵求法

时间;2014.05.15

地点:基地二楼

-----------------------------------------------------------------------

一、背景

  著名的斐波那契数列为一个这样的序列:0 1 1 2 3 5 8 13 21 34......简单的递推公式如下:

  F(0)=0,F(1)=1,当n>=1时,F(n)=F(n-1)+F(n-2)

显然,我们用直接的按公式递归的算法去计算该数列的第n项效率并不高,因为这样每次递归调用我们只是将规规模缩小了一点点,我们常用的递归都是将规模减半吧~,所以在这里这么用递归并不高效,我是终于又看到一个O(log)的算法了,还一个是二分查找,时间复杂度为O(logn)的算法并不多见,通过合理的算法设计,能将一个复杂的算法降到对数时间,真是神奇。在这里我将用到的是快速幂的思想。

-----------------------------------------------------------------------

二、快速幂算法

  前面有文章已经介绍过快速幂,通过求数幂来说明,比如求数7的5次方,我们知道5用二进制表示为101,于是我们按101这个序列从右至左,每次用7去乘当前7的幂,若是碰到1,还将当前得到的值与当前7的幂相乘。比如在这里7的5次方7(5)=7(1)*7(4)。计算过程为这样,从最右边开始,我们首先碰到1,于是将当前目标值为7(1),当前幂值平方为7(2),接下来碰到一个0,于是当前目标值不变,幂值 平方得7(4),再接下来碰到1,于是将当前目标值7(1)与当前幂值7(4)相乘,得新的目标值7(5),当前幂平方得7(8)。再接下来都是0了,就无需继续了,这可以通过将幂5的二进制形式不断右移实现。比如在这里,我们的斐波那契矩阵求法是基于等式:

|    F(n-1)     F(n)         |                           |           0       1              |  (n)

|                                  |          =               |                                    |                =B^n

|    F(n)        F(N+1)    |                           |            1       1             | 

那么我们这里即就是计算矩阵B的n次方而已,这个算法的伪代码如下:

long long CaculateFibonacci(size_t n)
{
  初始化目标矩阵为单位矩阵A,这里恰好对应F(0)=0
  初始化基矩阵为B的形式
  为了在计算矩阵乘法过程中保持原矩阵我们还需设定4个临时变量
  while(只要当前n非0)
  {
      if(当前n的最低位为1)
      {
        将目标矩阵A与当前基矩阵相乘
      }
      在每次循环中均将当前基矩阵平方处理,于是可得到A^2 A^4 A^8 A^16等
      将n右移一位
   }
}
      
三、代码实现如下:

#include<iostream>
#include<cassert>
using namespace std;
long long CaculateFibonacci(size_t n)
{
	//Precondition:
	//Postcondition:
	assert(n >= 0);
	long long ans_arr[2][2] = { 1, 0, 0, 1 };
	long long base_arr[2][2] = { 0, 1, 1, 1 };
	long long temp00 = 0, temp01 = 0, temp10 = 0,temp11 = 0; 
	while (n)
	{
		if (n & 1)  //当前最低位为1,将目标矩阵与当前基矩阵对应的幂阵相乘
		{
			temp00 = ans_arr[0][0], temp01 = ans_arr[0][1], temp10 = ans_arr[1][0], temp11 = ans_arr[1][1];
			ans_arr[0][0] = temp00 * base_arr[0][0] + temp01 * base_arr[1][0];
			ans_arr[0][1] = temp00 * base_arr[0][1] + temp01 * base_arr[1][1];
			ans_arr[1][0] = temp10 * base_arr[0][0] + temp11 * base_arr[1][0];
			ans_arr[1][1] = temp10 * base_arr[0][1] + temp11 * base_arr[1][1];
		}
	    //从基矩阵计算当前相应的幂阵
		temp00 = base_arr[0][0], temp01 = base_arr[0][1], temp10 = base_arr[1][0], temp11 = base_arr[1][1];
		base_arr[0][0] = temp00 * temp00 + temp01 * temp10;
		base_arr[0][1] = temp00 * temp01 + temp01 * temp11;
		base_arr[1][0] = temp10 * temp00 + temp11 * temp10;
		base_arr[1][1] = temp10 * temp01 + temp11 * temp11;
		n >>= 1;
	}
	return ans_arr[0][1];
}
int main()
{
	cout << "Input n: " << endl;
	int n = 0;
	cin >> n;
	long long a = CaculateFibonacci(n);
	cout <<"The result is: "<< a << endl;
}