矩阵快速幂求斐波那契通项(矩阵乘法优化线性递推式)

矩阵乘法求Fibo数列第N项(hihoCoder第41周)

题目描述：

骨牌，一种古老的玩具。今天我们要研究的是骨牌的覆盖问题：

我们有一个2xN的长条形棋盘，然后用1x2的骨牌去覆盖整个棋盘。对于这个棋盘，一共有多少种不同的覆盖方法呢？

举个例子，当N=1、2、3时

解法：

对于骨牌覆盖，我们可以考虑从已经覆盖了一些骨牌的情况推广到覆盖更多骨牌的情况

最右边的一种情况是不可能发生的，否则会始终多一个格子没有办法放置骨牌。或者说灰色部分的格子数为奇数，不可能通过1x2个骨牌放置出来。

那么通过对上面的观察，我们可以发现：

在任何一个放置方案最后，一定满足前面两种情况。而灰色的部分又正好对应了长度为N-1和N-2时的放置方案。由此，我们可以得到递推公式:

f[n] = f[n-1] + f[n-2];

这个公式是不是看上去很眼熟？没错，这正是我们的费波拉契数列。[1]

因此，问题被转化成为求斐波拉契数列第N项的问题。

对于数值较小的fibo数列，我们可以直接用递归求解:

int fibo(int x)
{
	if (x < 3) return (1);
	return (fibo(x - 1) + fibo(x - 2));
}

但是，对于数值较大的斐波那契数列，这种方法是无法完美处理的，因为fibo数列第100项的大小超过了long long的最大值 2^63-1。所以必须使用一种乘法来求解斐波那契数列的通项

因为乘法具有这样一个性质：

(a * b) mod N = ((a mod N) * (b * mod N)) mod N

因此，使用乘法处理斐波拉契数列的通项可以保证数值不越界。

实际上，对于线性递推式，我们可以使用矩阵乘法来求解。对于斐波那契数列，我们希望找到一个 2 * 2 的矩阵 M ，使得(a, b) X M = (b, a + b)

对于这样一个矩阵，我们可以设

M = [ x, y ; z, w]

得到方程

(a , b) X M = (a * x + b * z , a * y + b * w)

转化为方程组

a * x + b * z  = b;

a * y + b * w = a + b;

得到: 

x=0;y=1;z=1;w=1;

 

因此

进一步可以知道：

那么接下来的问题是，能不能快速的计算出M^n？我们先来分析一下幂运算。由于乘法是满足结合律的，所以我们有：

不妨将k[1]..k[j]划分的更好一点？

其中(k[1],k[2]...k[j])2表示将n表示成二进制数后每一位的数字。上面这个公式同时满足这样一个性质：

结合这两者我们可以得到一个算法：
1. 先计算出所有的{a^1, a^2, a^4 ... a^(2^j)}，因为该数列满足递推公式，时间复杂度为O(logN)
2. 将指数n二进制化，再利用公式将对应的a^j相乘计算出a^n，时间复杂度仍然为O(logN)
则总的时间复杂度为O(logN)
这种算法因为能够在很短时间内求出幂，我们称之为“快速幂”算法。

#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <vector>
#include <algorithm>
#define MOD 19999997
#define MAX 100000000

using namespace std;

/*
int fibo(int x)
{
	if (x < 3) return (1);
	return (fibo(x - 1) + fibo(x - 2));
}
*/

int to2(long long n, vector <int>& k)
{
	while (n)
	{
		k.push_back(n % 2);
		n /= 2;
	}
	reverse(k.begin(), k.end());//其实这一步是脱裤子打屁
	return 0;
}

long long quick_fibo (long long n)
{
	std::vector<int> v;
	to2(n, v);
	int j = v.size();
	long long a, b, c, d;
	long long x, y, z, w;
	a = d = 1; c = b = 0;
	x = 0; y = z = w = 1;
	for (int i = 0; i < j; i++)
	{
		if(v[j - i - 1])
		{
			long long tmp_a = ((a % MOD * x % MOD) % MOD + (b % MOD * z % MOD) % MOD) % MOD;
			long long tmp_b = ((a % MOD * y % MOD) % MOD + (b % MOD * w % MOD) % MOD) % MOD;
			long long tmp_c = ((c % MOD * x % MOD) % MOD + (d % MOD * z % MOD) % MOD) % MOD;
			long long tmp_d = ((c % MOD * y % MOD) % MOD + (d % MOD * w % MOD) % MOD) % MOD;
			a = tmp_a;
			b = tmp_b;
			c = tmp_c;
			d = tmp_d;
		}
		long long tmp_x = ((x % MOD * x % MOD) % MOD + (y % MOD * z % MOD) % MOD) % MOD;
		long long tmp_y = ((x % MOD * y % MOD) % MOD + (y % MOD * w % MOD) % MOD) % MOD;
		long long tmp_z = ((x % MOD * z % MOD) % MOD + (w % MOD * z % MOD) % MOD) % MOD;
		long long tmp_w = ((w % MOD * w % MOD) % MOD + (z % MOD * y % MOD) % MOD) % MOD;
		x = tmp_x;
		y = tmp_y;
		z = tmp_z;
		w = tmp_w;
	}
	return d;
}

int main(void)
{
	long long n;
	scanf("%lld", &n);
	cout << quick_fibo(n) << endl;
	return 0;
}