矩阵快速幂入门 + 求Fibonacci数列

两矩阵相乘，朴素算法的复杂度是O(N^3)。如果求一次矩阵的M次幂，按朴素的写法就是O(N^3*M)。既然是求幂，不免想到快速幂取模的算法，a^b %m 的复杂度可以降到O(logb)。如果矩阵相乘是不是也可以实现O(N^3 * logM)的时间复杂度呢？答案是肯定的。

先定义矩阵数据结构：　　

struct Mat{
    int mat[MaxN][MaxN];
};

O(N^3)实现一次矩阵乘法（最基本的一次矩阵乘法，类似于生成元数据）

Mat operator * (Mat a, Mat b){
    Mat c;
    memset(c.mat,0,sizeof(c.mat));
    
    for(int k=0;k< Nsize;k++){
        for(int i=0;i< Nsize;i++){
            
            for(int j=0;j< Nsize;j++){
               c.mat[i][j]=( c.mat[i][j] + a.mat[i][k]*b.mat[k][j])% MOD;
            }
        
        }
    
    }
    return c;
}

下面介绍一种特殊的矩阵：单位矩阵

很明显的可以推知，任何矩阵乘以单位矩阵，其值不改变。

有了前边的介绍，就可以实现矩阵的快速连乘了。

Mat operator ^ (Mat a, int k){
    Mat c;
    for(int i=0;i< Nsize;i++)
        for(int j=0;j< Nsize;j++)
            c.mat[i][j]=(i==j); //初始化为单位矩阵
            
    for(;k;k>>=1){
        if(k&1) c=c*a;
        a=a*a;
    }
    return c;
}

举个例子:

　　求第n个Fibonacci数模M的值。如果这个n非常大的话，普通的递推时间复杂度为O(n)，这样的复杂度很有可能会挂掉。这里可以用矩阵做优化，复杂度可以降到O(logn * 2^3)

如图：

A = F(n - 1), B = F(N - 2)，这样使构造矩阵的n次幂乘以初始矩阵得到的结果就是。

因为是2*2的据称，所以一次相乘的时间复杂度是O(2^3)，总的复杂度是O(logn * 2^3 + 2*2*1)。

POJ 3070 Fibonacci 的实现代码如下，思路就是上面部分：

#include<cstdio>
#include<cstring>
#include<iostream>
#define MaxN 105
#define MOD 10000
using namespace std;
int Nsize=2;

struct Mat{
    int mat[MaxN][MaxN];
};


Mat operator * (Mat a, Mat b){
    Mat c;
    memset(c.mat,0,sizeof(c.mat));
    
    for(int k=0;k< Nsize;k++){
        for(int i=0;i< Nsize;i++){
            
            for(int j=0;j< Nsize;j++){
               c.mat[i][j]=( c.mat[i][j] + a.mat[i][k]*b.mat[k][j])% MOD;
            }
        
        }
    
    }
    return c;
}


Mat operator ^ (Mat a, int k){
    Mat c;
    for(int i=0;i< Nsize;i++)
        for(int j=0;j< Nsize;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 n;
    Mat a;
    a.mat[0][0]=a.mat[0][1]=a.mat[1][0]=1;
    a.mat[1][1]=0;
    
    while(scanf("%d",&n)){
        if(n==-1) {break;}
        
        Mat res;
        res.mat[0][0]=1;
        res.mat[1][0]=0;
        res=(a^n)*res;
        printf("%d\n",res.mat[1][0]);
    }
    return 0;
}