思路：和整數快速冪一樣，唯一不同的就是存放結果的矩陣初始值為單位矩陣，通過過載運算子*後，程式碼可以大大簡化。另外需要注意的是取模問題，我把模M放在了全域性變數，這樣省卻一些麻煩，可以根據自身需要調整，這個無傷大雅。

程式碼示例：

#include<iostream>  
#include<cstdio>  
#include<cstring>  
using namespace std;  
const int maxn = 105;   
const int M = 1e9+7;  
struct Matrix{  
    int n,m;  
    int v[maxn][maxn];  
    Matrix(int n,int m):n(n),m(m){}  
    void init(){                    //初始化矩陣   
        memset(v,0,sizeof v);  
    }  
    Matrix operator* (const Matrix B) const {  
        Matrix C(n,B.m);            //用來存放答案   
        C.init();  
        for(int i = 0;i < n;i++)  
        for(int j = 0;j < B.m;j++)  
        for(int k = 0;k < m;k++)       
            C.v[i][j] = (C.v[i][j] + v[i][k]*B.v[k][j] % M) % M;  
        return C;  
    }  
    void print(){//輸出該矩陣，用來測試   
        for(int i = 0;i < n;i++){  
        for(int j = 0;j < m;j++)  
            cout << v[i][j] << " ";  
        cout << endl;  
        }  
    }  
};  
Matrix q_pow(Matrix& A,int b,int M){  
    Matrix C(A.n,A.m);  
    C.init();  
    for(int i = 0;i < C.n;i++)   C.v[i][i] = 1;  
    while(b){  
        if(b&1) C = C*A;  
        A = A*A;  
        b >>= 1;  
    }  
    return C;  
}