矩阵加速计算斐波那契数列

最新推荐文章于 2022-06-26 16:38:03 发布

Coder.AN

最新推荐文章于 2022-06-26 16:38:03 发布

阅读量1.6k

点赞数 2

分类专栏：数据结构与算法文章标签：算法矩阵

本文链接：https://blog.csdn.net/AnHongjun/article/details/112133081

版权

数据结构与算法专栏收录该内容

2 篇文章

订阅专栏

摘要

通过矩阵加速计算，将计算斐波那契数列的时间复杂度优化到 $O\left ( logn \right )$ ，空间复杂度 $O\left ( 1 \right )$ 。

一、斐波那契数列

斐波那契数列，又称黄金分割数列，是指这样一个数列：0，1，1，2，3，5，8，13，……

在数学上通常以如下递归方法定义： $F(0)=0, F(1)=1,F(n)=F(n-1)+F(n-2)(n\geq 2, n\in N+)$

二、传统计算方法

1.递归法

即使用最普通的递归算法进行计算，代码如下：

int Fib(int n){
    if(n == 0) return 0;
    if(n == 1) return 1;
    return Fib(n-1)+Fib(n-2);
}

该算法时间复杂度为 $O\left ( \left (\frac{2}{\sqrt{5} -1} \right ) ^{n}\right )$ ，时间复杂度推导可参考该篇博客点击跳转。可见当n较大时，时间复杂度将会非常大，满足不了时间要求。

2.记忆化递归（动态规划）

本质还是用递归的方法。但朴素递归存在重复计算的缺陷。如图是计算斐波那契数列的递归树。

可见，在递归过程中，很多量被重复计算过很多次。如果我们在计算的同时，将计算过的量存储起来，即加入“记忆化”，可以使计算变得相对高效。

int f[MaxN] = {0};

int Fib(int n){
    if(n == 0) return 0;
    if(n == 1) return 1;
    if(f[n]!=0) return f[n];
    return f[n] = Fib(n-1)+Fib(n-2);
}

由此算法时间复杂度可以优化到 $O\left ( n \right )$ ，但显然这种做法是以空间换时间，需要耗费相对较大的空间复杂度 $O\left ( n \right )$ 。虽然时间复杂度已经优化到线性，但面对较大的数据仍然略显吃力。

三、矩阵加速计算斐波那契数列

1.构造斐波那契矩阵

笔者假定读者已经熟悉矩阵乘法的计算规则，如果不是，请参考此处。

我们首先构造这样一个2x2的矩阵：

$\begin{bmatrix} F\left ( n+1 \right ) & F\left ( n \right )\\ F\left ( n \right )& F\left ( n-1 \right ) \end{bmatrix}$

再构造这样一个2x2的矩阵：

$\begin{bmatrix} 1&1 \\ 1&0 \end{bmatrix}$

如果我们把如上两个矩阵会得到什么效果？

$\begin{bmatrix} F\left ( n+1 \right ) & F\left ( n \right )\\ F\left ( n \right )& F\left ( n-1 \right ) \end{bmatrix} \times \begin{bmatrix} 1 &1 \\ 1&0 \end{bmatrix} = \begin{bmatrix} F\left ( n+1 \right ) +F\left ( n \right )&F\left ( n+1 \right ) \\ F\left ( n \right )+F\left ( n-1 \right )& F\left ( n \right ) \end{bmatrix} = \begin{bmatrix} F\left ( n+2 \right ) & F\left ( n+1 \right )\\ F\left ( n+1 \right ) & F\left ( n \right ) \end{bmatrix}$

由此我们可以得到一个递推公式，如果令：

$T\left ( n \right ) = \begin{bmatrix} F\left ( n+1 \right ) &F\left ( n \right ) \\ F\left ( n \right )& F\left ( n-1 \right ) \end{bmatrix}$

则：

$T\left ( n \right ) = T\left ( n-1 \right ) \times \begin{bmatrix} 1 &1 \\ 1 &0 \end{bmatrix} = T\left ( n-2 \right ) \times \begin{bmatrix} 1 &1 \\ 1&0 \end{bmatrix}^{2} =...... = \begin{bmatrix} 1&1 \\ 1& 0 \end{bmatrix}^{n}$

由此就出现了非常神奇的效果：

$\begin{bmatrix} 1&1 \\ 1&0 \end{bmatrix}^{1}= \begin{bmatrix} 1&1 \\ 1&0 \end{bmatrix}$ $F\left ( 1 \right )=1$

$\begin{bmatrix} 1&1 \\ 1&0 \end{bmatrix}^{2}= \begin{bmatrix} 2&1 \\ 1&1 \end{bmatrix}$ $F\left ( 2 \right )=1$

$\begin{bmatrix} 1&1 \\ 1&0 \end{bmatrix}^{3}= \begin{bmatrix} 3&2 \\ 2&1 \end{bmatrix}$ $F\left ( 3 \right )=2$

$\begin{bmatrix} 1&1 \\ 1&0 \end{bmatrix}^{4}= \begin{bmatrix} 5&3 \\ 3&2 \end{bmatrix}$ $F\left ( 4 \right )=3$

$\begin{bmatrix} 1&1 \\ 1&0 \end{bmatrix}^{5}= \begin{bmatrix} 8&5 \\ 5&3 \end{bmatrix}$ $F\left ( 5 \right )=5$

$\begin{bmatrix} 1&1 \\ 1&0 \end{bmatrix}^{6}= \begin{bmatrix} 13&8 \\ 8&5 \end{bmatrix}$ $F\left ( 6 \right )=8$

……

$\begin{bmatrix} 1&1 \\ 1&0 \end{bmatrix}^{n}= \begin{bmatrix} F\left ( n+1 \right )&F\left ( n \right ) \\ F\left ( n \right )&F\left ( n-1 \right ) \end{bmatrix}$ $F\left ( n \right )=F\left ( n \right )$

因此，如果我们要计算 $F\left ( n \right )$ ，我们只需要计算 $\begin{bmatrix} 1 & 1\\ 1&0 \end{bmatrix}^{n}$ ，然后取左下角（或右上角）的元素即可。

笔者将 $\begin{bmatrix} 1&1 \\ 1&0 \end{bmatrix}$ 称为单位斐波那契矩阵。

2.C++代码实现矩阵计算斐波那契数列

用结构体构造矩阵。

struct Matrix{
    int data[2][2];
}

构造一个成员初始化函数，方便对该数据类型进行初始化

struct Matrix{
    int data[2][2];
    Matrix(int data00, int data01, int data10, int data11){
        data[0][0] = data00; data[0][1] = data01;
        data[1][0] = data10; data[1][1] = data11;
    }
};

该成员函数构造好后，就可以用如下方法初始化Matrix类型变量了。

Matrix Fib = Matrix(1, 1, 1, 0);

对乘法运算符*重载，方便矩阵乘法运算。

struct Matrix{
    static const int MOD = 10007;
    int data[2][2];
    Matrix(int data00, int data01, int data10, int data11){
        data[0][0] = data00; data[0][1] = data01;
        data[1][0] = data10; data[1][1] = data11;
    }
    Matrix operator * (const Matrix &a) const{
        Matrix ans = Matrix(0, 0, 0, 0);
        for(int i = 0; i<2; ++i)
            for(int j = 0; j<2; ++j)
                for(int k = 0; k<2; ++k) 
                    ans.data[i][j] = (ans.data[i][j] + data[i][k]*a.data[k][j])%MOD;
        return ans;
    }
};

这样我们就可以像这样直接调用*运算符进行矩阵乘法运算了。

Matrix Fib = Matrix(1, 1, 1, 0);
Fib = Fib * Fib;
printf("%d\n", Fib.data[1][0]);  //F(2)

构造成员函数进行幂运算。

struct Matrix{
    static const int MOD = 10007;
    int data[2][2];
    Matrix(int data00, int data01, int data10, int data11){
        data[0][0] = data00; data[0][1] = data01;
        data[1][0] = data10; data[1][1] = data11;
    }
    Matrix operator * (const Matrix &a) const{
        Matrix ans = Matrix(0, 0, 0, 0);
        for(int i = 0; i<2; ++i)
            for(int j = 0; j<2; ++j)
                for(int k = 0; k<2; ++k) 
                    ans.data[i][j] = (ans.data[i][j] + data[i][k]*a.data[k][j])%MOD;
        return ans;
    }
    Matrix pow(int power){
        Matrix ans = (*this);
        for(int i = 2; i<=power; ++i)
            ans = ans * (*this);
        return ans;
    }
};

由此我们便能实现斐波那契数列的计算。

#include <cstdio>
using namespace std;

struct Matrix{
    static const int MOD = 10007;
    int data[2][2];
    Matrix(int data00, int data01, int data10, int data11){
        data[0][0] = data00; data[0][1] = data01;
        data[1][0] = data10; data[1][1] = data11;
    }
    Matrix operator * (const Matrix &a) const{
        Matrix ans = Matrix(0, 0, 0, 0);
        for(int i = 0; i<2; ++i)
            for(int j = 0; j<2; ++j)
                for(int k = 0; k<2; ++k) 
                    ans.data[i][j] = (ans.data[i][j] + data[i][k]*a.data[k][j])%MOD;
        return ans;
    }
    Matrix pow(int power){
        Matrix ans = (*this);
        for(int i = 2; i<=power; ++i)
            ans = ans * (*this);
        return ans;
    }
};

int main(){
    int n;
    scanf("%d", &n);
    Matrix ans = Matrix(1, 1, 1, 0);
    printf("%d\n", ans.pow(n).data[1][0]);
    return 0;
}

算法时间复杂度为 $O\left ( n \right )$ ，空间复杂度优化到 $O\left ( 1 \right )$ 。时间复杂度我们可以继续优化。

3.快速幂

关于快速幂算法可以参考这篇文章。

现在我们将快速幂推广到矩阵乘法。

由于n较大时运算结果可能非常大，在此我们对结果进行10007取余运算。

矩阵乘法作为一个代数系统，不难证明 $\begin{bmatrix} 1 &0 \\ 0& 1 \end{bmatrix}$ 是他的幺元。即 $\begin{bmatrix} a &b \\ c & d \end{bmatrix}^{0}} = \begin{bmatrix} 1& 0\\ 0&1 \end{bmatrix}$ 。

struct Matrix{
    static const int MOD = 10007;
    int data[2][2];
    Matrix(int data00, int data01, int data10, int data11){
        data[0][0] = data00; data[0][1] = data01;
        data[1][0] = data10; data[1][1] = data11;
    }
    Matrix operator * (const Matrix &a) const{
        Matrix ans = Matrix(0, 0, 0, 0);
        for(int i = 0; i<2; ++i)
            for(int j = 0; j<2; ++j)
                for(int k = 0; k<2; ++k) 
                    ans.data[i][j] = (ans.data[i][j] + data[i][k]*a.data[k][j])%MOD;
        return ans;
    }
    Matrix pow( int power){
        Matrix ans = Matrix(1, 0, 0, 1);
        Matrix base = (*this);
        while(power){
            if(power&1) ans = ans * base;
            base = base * base;
            power >>= 1;
        }
        return ans;
    }
};

由此我们将时间复杂度优化到 $O\left ( logn \right )$ ，空间复杂度为 $O\left ( 1 \right )$ ，可以说是一个非常优的算法了。

4.矩阵加速计算斐波那契数列全部代码

#include <cstdio>
using namespace std;

struct Matrix{
    static const int MOD = 10007;
    int data[2][2];
    Matrix(int data00, int data01, int data10, int data11){
        data[0][0] = data00; data[0][1] = data01;
        data[1][0] = data10; data[1][1] = data11;
    }
    Matrix operator * (const Matrix &a) const{
        Matrix ans = Matrix(0, 0, 0, 0);
        for(int i = 0; i<2; ++i)
            for(int j = 0; j<2; ++j)
                for(int k = 0; k<2; ++k) 
                    ans.data[i][j] = (ans.data[i][j] + data[i][k]*a.data[k][j])%MOD;
        return ans;
    }
    Matrix pow( int power){
        Matrix ans = Matrix(1, 0, 0, 1);
        Matrix base = (*this);
        while(power){
            if(power&1) ans = ans * base;
            base = base * base;
            power >>= 1;
        }
        return ans;
    }
};

int main(){
    int n;
    scanf("%d", &n);
    Matrix ans = Matrix(1, 1, 1, 0);
    printf("%d\n", ans.pow(n).data[1][0]);
    return 0;
}