矩陣快速冪 總結

剛做了一道矩陣快速冪的題，看了網上不少資料，決定整理一下，接下來再做的時候也可以參考。從網上各位大神那邊直接copy過來的

 

 

矩陣快速冪

矩陣的快速冪是用來高效地計算矩陣的高次方的。將朴素的o（n）的時間復雜度，降到log（n）。

這里先對原理（主要運用了矩陣乘法的結合律）做下簡單形象的介紹：

一般一個矩陣的n次方，我們會通過連乘n-1次來得到它的n次冪。

但做下簡單的改進就能減少連乘的次數，方法如下：

把n個矩陣進行兩兩分組，比如：A*A*A*A*A*A  =>  (A*A)*(A*A)*(A*A)

這樣變的好處是，你只需要計算一次A*A，然后將結果(A*A)連乘自己兩次就能得到A^6，即(A*A)^3=A^6。算一下發現這次一共乘了3次，少於原來的5次。

其實大家還可以取A^3作為一個基本單位。原理都一樣：利用矩陣乘法的結合律，來減少重復計算的次數。

以上都是取一個具體的數來作為最小單位的長度，這樣做雖然能夠改進效率，但缺陷也是很明顯的，取個極限的例子（可能有點不恰當，但基本能說明問題），當n無窮大的時候，你現在所取的長度其實和1沒什么區別。所以就需要我們找到一種與n增長速度”相適應“的”單位長度“，那這個長度到底怎么去取呢？？？這點是我們要思考的問題。

有了以上的知識，我們現在再來看看，到底怎么迅速地求得矩陣的N次冪。

既然要減少重復計算，那么就要充分利用現有的計算結果咯！~怎么充分利用計算結果呢？？？這里考慮二分的思想。。

大家首先要認識到這一點：任何一個整數N，都能用二進制來表示。。這點大家都應該知道，但其中的內涵真的很深很深（這點筆者感觸很深，在文章的最后，我將談談我對的感想）！！

計算機處理的是離散的信息，都是以0,1來作為信號的處理的。可想而知二進制在計算機上起着舉足輕重的地位。它能將模擬信號轉化成數字信號，將原來連續的實際模型，用一個離散的算法模型來解決。  好了，扯得有點多了，不過相信這寫對下面的講解還是有用的。

回頭看看矩陣的快速冪問題，我們是不是也能把它離散化呢？比如A^19  =>  （A^16）*（A^2）*（A^1），顯然采取這樣的方式計算時因子數將是log(n)級別的(原來的因子數是n)，不僅這樣，因子間也是存在某種聯系的，比如A^4能通過(A^2)*(A^2)得到，A^8又能通過(A^4)*(A^4)得到，這點也充分利用了現有的結果作為有利條件。下面舉個例子進行說明：

現在要求A^156,而156(10)=10011100(2) 

也就有A^156=>(A^4)*(A^8)*(A^16)*(A^128)  考慮到因子間的聯系，我們從二進制10011100中的最右端開始計算到最左端。細節就說到這，下面給核心代碼：

while(N)
{
    if(N&1)
        res=res*A;
    n>>=1;
    A=A*A;
}

里面的乘號，是矩陣乘的運算，res是結果矩陣。

第3行代碼每進行一次，二進制數就少了最后面的一個1。二進制數有多少個1就第3行代碼就執行多少次。

好吧，矩陣快速冪的講解就到這里吧。在文章我最后給出我實現快速冪的具體代碼（代碼以3*3的矩陣為例）。

現在我就說下我對二進制的感想吧：

我們在做很多”連續“的問題的時候都會用到二進制將他們離散簡化

1.多重背包問題

2.樹狀數組

3.狀態壓縮DP

……………還有很多。。。究其根本還是那句話：化連續為離散。。很多時候我們並不是為了解決一個問題而使用二進制，更多是時候是為了優化而使用它。所以如果你想讓你的程序更加能適應大數據的情況，那么學習學習二進制及其算法思想將會對你有很大幫助。

 

 

快速冪模板（整數+矩陣）

1//整數的快速冪 m^n  % k 的快速冪： 

long long  quickpow(long long   m , long long   n , long long   k){ 
    long long   ans = 1; 
    while(n){ 
        if(n&1)//如果n是奇數 
            ans = (ans * m) % k; 
        n = n >> 1;//位運算“右移1類似除1” 
        m = (m * m) % k; 
    } 
    return ans; 
} 

2矩陣快速冪： 
定義一個矩陣類，例如求（A^n）%mod 

class Matrix { 
public: 
 
     long long m[MAXN][MAXN]; 
//二維數組存放矩陣 
    Matrix(){} 
    //對數組的初始化 
    void init(long long  num[MAXN][MAXN]){ 
        for(int i = 0 ; i < MAXN ; i++){ 
            for(int j = 0 ; j < MAXN ; j++){ 
                m[i][j] = num[i][j]; 
           } 
       } 
    } 
    //重載矩陣的乘法運算 
 
    friend Matrix operator*(Matrix &m1 ,Matrix &m2) { 
        int i, j, k; 
        Matrix temp; 
        for (i = 0; i < MAXN; i++) { 
            for (j = 0; j < MAXN; j++) { 
                temp.m[i][j] = 0; 
                for(k = 0 ; k < MAXN ; k++) 
                   temp.m[i][j] += (m1.m[i][k] * m2.m[k][j])%mod 
                temp.m[i][j] %= mod; 
//注意每一步都進行取模 
           } 
        } 
        return temp; 
    } 
    //矩陣的快速冪 
 
    friend Matrix quickpow(Matrix &M , long long n){ 
        Matrix tempans; 
//初始化為單位矩陣 
        //初始化 
        for(int i = 0 ; i < MAXN ; i++){ 
            for(int j = 0 ; j < MAXN ; j++){ 
                if(i == j) 
                    tempans.m[i][j] = 1; 
                else 
                    tempans.m[i][j] = 0; 
            } 
        } 
        //快速冪（類似整數） 
        while(n){ 
            if(n & 1)    www.2cto.com
                tempans = tempans * M; 
//已經重載了* 
            n = n >> 1; 
            M = M * M; 
        } 
       return tempans; 
    } 
}; 
 
int main() { 
    Matrix A , ans; 
    long long T , n , k , sum; 
//數據類型為long long 
    long long num[MAXN][MAXN]; 
//輸入的數據存入數組 
    scanf("%lld" , &T); 
    while(T--){ 
        scanf("%lld%lld/n", &n , &k); 
        memset(num , 0 , sizeof(num)); 
        for(int i = 0 ; i < n ; i++){ 
            for(int j = 0 ; j < n ; j++) 
                scanf("%lld" , &num[i][j]); 
        } 
        A.init(num);//初始化A矩陣 
        ans = quickpow(A , k);//求出矩陣的快速冪 
    } 
} 

 

 

最后，還有剛哥整理的，鏈接

 

 芳姐的矩陣快速冪的模板

利用快速冪的思想 根據矩陣的結合律 可以遞歸二分求解 

struct Mat
{
    int mat[N][N];
};
int 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;//優化
            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;
}