1. 程式人生 > >矩陣快速模冪 + 求斐波那契數列第n項(Fibonacci)

矩陣快速模冪 + 求斐波那契數列第n項(Fibonacci)

兩矩陣相乘,樸素演算法的複雜度是O(N^3)。如果求一次矩陣的M次冪,按樸素的寫法就是O(N^3*M)。既然是求冪,不免想到快速冪取模的演算法,前面有快速冪取模的介紹,a^b %m 的複雜度可以降到O(logb)。如果矩陣相乘是不是也可以實現O(N^3 * logM)的時間複雜度呢?答案是肯定的。

過載結構體運算子*和^表示矩陣乘法和矩陣求冪,取模的過程在 乘法中實現。

(Fibonacci)A = F(n - 1), B = F(N - 2),這樣使構造矩陣的n次冪乘以初始矩陣得到的結果就是

因為是2*2的據稱,所以一次相乘的時間複雜度是O(2^3),總的複雜度是O(logn * 2^3 + 2*2*1)。

#include <iostream>
#include <cstdio>
#include <string>
#include <cstring>
#include <fstream>
#include <algorithm>
#include <cmath>
#include <queue>
#include <stack>
#include <vector>
#include <map>
#include <set>
#include <iomanip>

using namespace std;
#define MOD 1000000007
#define N 2

struct Mat
{
    long long mat[N][N];
};

Mat operator * (Mat a, Mat b)
{
    Mat c;
    memset(c.mat, 0, sizeof(c.mat));
    int i, j, k;
    for(k = 0; k < N; ++k)
    {
        for(i = 0; i < N; ++i)
        {
            if(a.mat[i][k] <= 0)  continue;   //(針對ZOJ2853)剪枝,cpu運算乘法的效率並不是想像的那麼理想(加法的運算效率高於乘法,比如Strassen矩陣乘法)
            for(j = 0; j < N; ++j)
            {
                if(b.mat[k][j] <= 0)    continue;    //剪枝
                c.mat[i][j] = (c.mat[i][j] + (a.mat[i][k] * b.mat[k][j]) % MOD) % MOD;
            }
        }
    }
    return c;
}

Mat operator ^ (Mat a, int k)
{
    Mat c;
    int i, j;
    for(i = 0; i < N; ++i)
        for(j = 0; j < N; ++j)
            c.mat[i][j] = (i == j);    //初始化為單位矩陣

    for(; k; k >>= 1)
    {
        if(k&1) c = c*a;
        a = a*a;
    }
    return c;
}

int main()
{
    int a , b , n;
    while(cin >> a >> b >> n)
    {
        Mat m;
        m.mat[0][0] = 1;
        m.mat[0][1] = 1;
        m.mat[1][0] = 1;
        m.mat[1][1] = 0;
        if(n >= 1)
        {
            m = m ^ (n);
          /*  for(int i = 0; i < N; ++i)
            {
                for(int j = 0; j < N; ++j)
                    cout << m.mat[i][j] <<" ";
                cout << endl;
            }*/

            cout <<((m.mat[0][0] * a )% MOD+ (m.mat[1][0] * b )%MOD )%MOD << endl;
        }
        else cout << a << endl;
    }
    return 0 ;
}