Fibonacci数列——快速矩阵幂
今天想起来之前一个oj题目,是求类似与斐波那契数列一个数列的第N位。那时候接触到一个算法叫快速矩阵幂。
在这里我就用斐波那契数列的列子来简单说明一下如何用快速矩阵幂来解决这个题目。
Fibonacci数列定义:\(F(0) = 1, F(1) = 1, F(2) = 1, F(3) = 2, ...... F(n) = F(n-1)+F(n-2)\)
首先说明一下,因为斐波那契数列增长速度非常迅速,得到的数字可能过大,因此我们将结果对10000007(\(10^7+7\))取余来进行对比。
最天真的做法是用递归来解决:
1 | long long fibNaive(long long n) |
不用说了,算法第一步就会介绍这个反例,来说明递归效率不一定会高(他的算法的运行时间随n的增长类似与Fibonacci数列的增长)。实际上这个做法到n = 40的时候就已经可以让人等的有点急了。
然后正常的做法是用简单的循环
用两个数来代表之前的两个值,求出新值后继续依次更新前两个值,直到得到正确的结果:
1 | long long fibNormal(long long n) |
这个算法时间复杂度是\(O(n)\)。\(O(n)\)已经是一个很好的复杂度了,那还有没有办法继续加快这个过程?? ## 快速矩阵幂 ## 观察斐波那契数列的生成过程,我们可以发现它们可以被写成下面的样子: \[ F(N) = F(N-1) + F(N-2) F(N-1) = F(N-1)+0*F(N-2) \] 上面的式子可以写成矩阵形式: \[ \left [ \begin{matrix} F(N)\\ F(N-1) \end{matrix} \right ] = \begin{bmatrix} 1&1\\ 1&0 \end{bmatrix} \begin{bmatrix} F(N-1)\\ F(N-2) \end{bmatrix} \] 不断重复上面过程,往后继续展开,我们可以得到: \[ \left [ \begin{matrix} F(N)\\ F(N-1) \end{matrix} \right ] = { \begin{bmatrix} 1&1\\ 1&0 \end{bmatrix} }^{n-1} \begin{bmatrix} F(1)\\ F(0) \end{bmatrix} \] 因此我们可以把重点放到怎么来求中间这个矩阵的幂。而快速矩阵幂的思想也很简单,就类似与对于数字的幂的求法一致。比如:\(X^9 = X^8 \cdot X\),而\(X^8 = (X^4)^2 = ((X^2)^2)^2\),因此需要3+1次乘法就可以算出来8次幂,容易看出来快速矩阵幂的时间复杂度是\(O(\log (n))\)。 因此利用快速矩阵幂,可以将斐波那契数列的求法进一步加速。
至于如何实现就是细节问题了,需要注意的时候是乘法取余数的时候,两方都小于10000007,乘积依然可能会超过int的范围(10000007*10000007),导致出错,因此我在这里选择long long类型,这样可以保证结果的正确性。
实现分为几步1:定义矩阵乘法: 1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17typedef vector<vector<long long>> Matrix;
long long d = 10000007;
Matrix m_mul(const Matrix &m,const Matrix &n)
{
//check(m,n);
Matrix result = vector<vector<long long>>(m.size(),vector<long long>(n[0].size()));
for(long long i = 0;i!=m.size();++i)
for(long long j = 0;j!=n[0].size();++j)
{
long long temp = 0;
for(long long k = 0;k!=n.size();++k )
temp = ((m[i][k]*n[k][j])%d + temp%d)%d;
result[i][j] = temp;
}
return result;
}1
2
3
4
5
6
7
8
9
10 Matrix help(const Matrix & m,long long n)
{
Matrix result;
if(n == 1)
return m;
else if(n == 2)
return m_mul(m,m);
result = help(m,n/2);
return m_mul(result,result);
}
1 | Matrix quickMatrixPower(const Matrix &m,long long n) |
最后用fib函数封装起来:
1 | long long fib(long long n) |
最后用main函数利用clock函数进行时间测试 1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23int main()
{
long long n,result;
double start,end;
while(cin>>n)
{
start = clock();
result = fib(n);
end = clock();
cout<<result<<" "<<end-start<<endl;
start = clock();
result = fibNormal(n);
end = clock();
cout<<result<<" "<<end-start<<endl;
start = clock();
if(n<45)
{
result = fibNaive(n);
end = clock();
cout<<result<<" "<<end-start<<endl;
}
}
}1
2
3
4
5
6
7
8
9
10输入:40
2334085 0
2334085 0
2334085 6956
输入:1000000//此时naive的算法已经无法求出来
9640841 0
9640841 19
输入:100000000
129680 0
129680 3295