1. 程式人生 > >淺談矩陣快速冪的那些事

淺談矩陣快速冪的那些事

矩陣的快速冪是用來高效地計算矩陣的高次方的。將樸素的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

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

程式碼如下:

#include <cstdlib>
#include <cstring>
#include <cstdio>
#include <iostream> 
using namespace std;
int N;
struct matrix
{
       int a[3][3];
}origin,res;
matrix multiply(matrix x,matrix y)
{
       matrix temp;
       memset(temp.a,0,sizeof(temp.a));
       for(int i=0;i<3;i++)
       {
               for(int j=0;j<3;j++)
               {
                       for(int k=0;k<3;k++)
                       {
                               temp.a[i][j]+=x.a[i][k]*y.a[k][j];
                       }
               }
       }
       return temp;
}
void init()
{
     printf("隨機陣列如下:\n");
     for(int i=0;i<3;i++)
     {
             for(int j=0;j<3;j++)
             {
                     origin.a[i][j]=rand()%10;
                     printf("%8d",origin.a[i][j]);
             }
             printf("\n");
     }//
     printf("\n");
     memset(res.a,0,sizeof(res.a));
     res.a[0][0]=res.a[1][1]=res.a[2][2]=1;                  //將res.a初始化為單位矩陣 
}
void calc(int n)
{
     while(n)
     {
             if(n&1)
                    res=multiply(res,origin);
             n>>=1;
             origin=multiply(origin,origin);
     }
     printf("%d次冪結果如下:\n",n);
     for(int i=0;i<3;i++)
     {
             for(int j=0;j<3;j++)
                     printf("%8d",res.a[i][j]);
             printf("\n");
     }
     printf("\n");
}
int main()
{
    while(scanf("%d",&N)
    {
            init();
            calc(N);
    }
    return 0;
}

三維的矩陣快速冪,適用於所有的快速冪,這裡還有二維的矩陣,作用類似三維。但是在單位矩陣賦值的時候,必須賦成一樣的數,造成的結果是隻可以計算2的倍數的快速冪。(為什麼呢?在二維的計算中,矩陣運算的基本演算法造成了這一現象)

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<iostream>
#include<cmath>
#include<string>
#include<queue>
#include<map>
using namespace std;
const int MOD=10000;
struct matrix
{
    int m[2][2];
}base,one;
matrix calculate(matrix a,matrix b)
{
    matrix ans;
    for(int i=0;i<2;++i)
    for(int j=0;j<2;++j)
    {
        ans.m[i][j]=0;//初始賦值0 
        for(int k=0;k<2;++k)
        ans.m[i][j]=(ans.m[i][j]+a.m[i][k]*b.m[k][j])%MOD;
    }
    return ans;
}//矩陣乘法運算函式 
int fastpow(int n)
{
    one.m[0][0]=one.m[0][1]=one.m[1][0]=1;
    one.m[1][1]=0;
    base.m[0][0]=base.m[1][1]=1;  
    base.m[0][1]=base.m[1][0]=0;//base初始化為單位矩陣 
    while(n!=0)
    {
        if(n&1)//末位是否為一 
        {
            base=calculate(base,one);
        }//實現base*=t;其中要先把base賦值給ans,然後用base=ans*t 
        one=calculate(one,one);
        n=n/2;
    }
    return base.m[0][1];
}//求矩陣one的n次冪 
int main()
{
    int n;
    while(true)
    {   
        scanf("%d",&n);
        if (n==-1) break;
        else printf("%d\n",fastpow(n));
    }
    return 0;
}

這個快速冪還有些問題,以後更新吧。
2017.4.13