在網站上一直沒有找到有關於快速冪演算法的一個詳細的描述和解釋，這裡，我給出快速冪演算法的完整解釋，用的是C語言，不同語言的讀者只好換個位啦，畢竟讀C的人較多~

所謂的快速冪，實際上是快速冪取模的縮寫，簡單的說，就是快速的求一個冪式的模(餘)。在程式設計過程中，經常要去求一些大數對於某個數的餘數，為了得到更快、計算範圍更大的演算法，產生了快速冪取模演算法。[有讀者反映在講快速冪部分時有點含糊，所以在這裡對本文進行了修改，作了更詳細的補充，爭取讓更多的讀者一目瞭然]

我們先從簡單的例子入手：求= 幾。

演算法1.首先直接地來設計這個演算法：

int ans = 1;

for(int i = 1;i<=b;i++)

{

ans = ans * a;

}

ans = ans % c;

這個演算法的時間複雜度體現在for迴圈中，為O（b）.這個演算法存在著明顯的問題，如果a和b過大，很容易就會溢位。

那麼，我們先來看看第一個改進方案：在講這個方案之前，要先有這樣一個公式：

.這個公式大家在離散數學或者數論當中應該學過，不過這裡為了方便大家的閱讀，還是給出證明：

引理1：

上面公式為下面公式的引理，即積的取餘等於取餘的積的取餘。

證明了以上的公式以後，我們可以先讓a關於c取餘，這樣可以大大減少a的大小，

於是不用思考的進行了改進：

演算法2：

int ans = 1;

a = a % c; //加上這一句

for

(int i = 1;i<=b;i++)

{

ans = ans * a;

}

ans = ans % c;

聰明的讀者應該可以想到，既然某個因子取餘之後相乘再取餘保持餘數不變，那麼新算得的ans也可以進行取餘，所以得到比較良好的改進版本。

演算法3：

int ans = 1;

a = a % c; //加上這一句

for(int i = 1;i<=b;i++)

{

ans = (ans * a) % c;//這裡再取了一次餘

}

ans = ans % c;

這個演算法在時間複雜度上沒有改進，仍為O(b)，不過已經好很多的，但是在c過大的條件下，還是很有可能超時，所以，我們推出以下的快速冪演算法

。

快速冪演算法依賴於以下明顯的公式，我就不證明了。

有了上述兩個公式後，我們可以得出以下的結論：

1.如果b是偶數，我們可以記k = a2 mod c，那麼求(k)b/2 mod c就可以了。

2.如果b是奇數，我們也可以記k = a2 mod c，那麼求

((k)b/2 mod c × a ) mo d c =((k)b/2 mod c * a) mod c 就可以了。

那麼我們可以得到以下演算法：

演算法4：

int ans = 1;

a = a % c;

if(b%2==1)

ans = (ans * a) mod c; //如果是奇數，要多求一步，可以提前算到ans中

k = (a*a) % c; //我們取a2而不是a

for(int i = 1;i<=b/2;i++)

{

ans = (ans * k) % c;

}

ans = ans % c;

我們可以看到，我們把時間複雜度變成了O(b/2).當然，這樣子治標不治本。但我們可以看到，當我們令k = (a * a) mod c時，狀態已經發生了變化，我們所要求的最終結果即為(k)b/2 mod c而不是原來的ab mod c，所以我們發現這個過程是可以迭代下去的。當然，對於奇數的情形會多出一項a mod c，所以為了完成迭代，當b是奇數時，我們通過

ans = (ans * a) % c;來彌補多出來的這一項，此時剩餘的部分就可以進行迭代了。

形如上式的迭代下去後，當b=0時，所有的因子都已經相乘，演算法結束。於是便可以在O（log b）的時間內完成了。於是，有了最終的演算法：快速冪演算法。

演算法5：快速冪演算法

int ans = 1;

a = a % c;

while(b>0)

{

if(b % 2 == 1)

ans = (ans * a) % c;

b = b/2;

a = (a * a) % c;

}

將上述的程式碼結構化，也就是寫成函式：

int PowerMod(int a, int b, int c)

{

int ans = 1;

a = a % c;

while(b>0)

{

if(b % 2 = = 1)

ans = (ans * a) % c;

b = b/2;

a = (a * a) % c;

}

return ans;

}

本演算法的時間複雜度為O（logb），能在幾乎所有的程式設計（競賽）過程中通過，是目前最常用的演算法之一。

以下內容僅供參考：

擴充套件：有關於快速冪的演算法的推導，還可以從另一個角度來想。

=？ 求解這個問題，我們也可以從進位制轉換來考慮：

將10進位制的b轉化成2進位制的表示式：

那麼，實際上，.

所以

注意此處的要麼為0，要麼為1，如果某一項，那麼這一項就是1，這個對應了上面演算法過程中b是偶數的情況，為1對應了b是奇數的情況[不要搞反了，讀者自己好好分析，可以聯絡10進位制轉2進位制的方法]，我們從依次乘到。對於每一項的計算，計算後一項的結果時用前一項的結果的平方取餘。對於要求的結果而言，為時ans不用把它乘起來，[因為這一項值為1]，為1項時要乘以此項再取餘。這個演算法和上面的演算法在本質上是一樣的，讀者可以自行分析，這裡我說不多說了，希望本文有助於讀者掌握快速冪演算法的知識點，當然，要真正的掌握，不多練習是不行的。

矩陣 快速求冪 blog1

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

 }

改進版（取模）a^b mod n ：

int modexp(int a, int b, int n)

{

   int t = a, ret = 1;

   while(b != 0)

  {

     if(b % 2 == 1) ret *= t % n;

     t = t * t % n;

     b /= 2;

  } 

  return ret;

}

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

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

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

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

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

1.多重揹包問題

2.樹狀陣列

3.狀態壓縮DP

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

快速冪或者矩陣快速冪在算大指數次方時是很高效的，他的基本原理是二進位制，下面的A可以是一個數也可以是一個矩陣（本文特指方陣），若是數就是快速冪演算法，若是矩陣就是矩陣快速冪演算法，用c++只需把矩陣設成一個類就可以，然後過載一下乘法就可以，注意為矩陣時則ANS=1，應該是ANS=E，E是單位矩陣，即主對角線是1其餘的部分都是0的特殊方陣了。

 舉個例子若你要算A^7你會怎麼算一般你會用O(N)的演算法A^7=A*A*A*A*A*A*A也許你覺得這並不慢但是若要你算A^10000000000000000呢，是不是會覺得O(N)的演算法也太慢了吧這不得算死我啊，計算機也不想算了，因為有更高效的演算法我們把A的指數寫成二進位制，這樣就有了

A^7=A^111(2)  現在我們可以這麼算 令ANS=1;MULTI=A  ,N=7

while(N)

{

if(N%2==1) //亦可以寫成N&1 或N%2

ANS*=MULTI;

MULTI*=MULTI;

N/=2;//c++中可以寫成 N>>=1;直接用位運算更快

}

寫出上面的程式碼的執行過程就是

ANS=1;MULTI=A; 

N=7 ;N%2=1;   ANS*=MULTI; 所以ANS=A;  MULTI*=MULTI; 所以MULTI=A^2

然後 N/=2;N=3; N%2=1; ANS*=MULTI; 所以 ANS=A*A^2=A^3 ; 又MULTI*=MULTI; 所以MULTI=A^4

然後N/=2;N=1;N%2=1;ANS*=MULTI; 所以 ANS=A*A^2*A^4=A^7;又MULTI*=MULTI; 所以MULTI=A^8

然後N/=2;N=0;演算法結束 是不是很巧妙呢，實際上用的乘法次數是 6次你可能覺得，那個A^7=A*A*A*A*A*A*A，不也是用了6次乘法嗎有什麼區別？

那是因為這個演算法是log2(n)   (表示以2為底n的對數) 的複雜度，還有一個係數，大約是2 實際上計算次數就是 2*log2(n) 而普通的連乘計算的複雜度是n 乘法計算次數是n-1

這樣在n很小時差別不大，但隨著n的增長差距會迅速擴大，例如 n=1024時 普通方法得計算1023次乘法，但快速冪最多（因為當上面的程式執行時N的中間結果為偶數那麼  ANS*=MULTI，將不會被執行，故實際的計算次數要小於 2*log2(n)）只算2*log2(n) =20次乘法，是不是很快！！！！！！！！！！

但是為什麼呢？好像還有點不懂。。

實際上A^7=A^1*A^2*A^4這樣每次計算乘法乘的因子都是遞增的，而且還是指數遞增，還有這些因子是可以遞推產生的就是可以利用上次的計算每次平方就可以了，這中其實是使用的二進位制的思想，因為任意一個數都可以，表示成二進位制，故 A^N以定可以寫成

A^(一個二進位制數如101010)=A^(100000)*A^(00000)*A(1000)*A^(000)*A^(10)*A^(0)=A^(2^5)*A^(2^3)*A^(2^1)

而我們的MULTI 其實是一個數列 A^1，A^2，A^4，A^8，A^16........即A^(2^0),A^(2^1),A^(2^2)，A^(2^3),A^(2^4).......................注意到他的指數都是二進位制的位權（不知道是不是這個名詞，就像十進位制的位權是 1 10 100 1000 10000，一樣如1243=1*1000+2*100+3*10+4*1；而二進位制的1011 是 1*2^(3)+0*2^(2)+1*2^(1)+1*2^(0) 這樣是不是應該理解位權了呢？）實際上任何一個A^N都可以寫成這個數列的某些項的乘積，因為N始終都可以表示成二進位制，而把N表示成二進位制後如果某項為1則說明需要乘上MULTI 否則不用乘上MULTI

於是就有了上面的程式碼，，，，哎怎麼感覺我說的還是很不清楚呢?那就沒辦法

下面附上程式碼，另外一般要用快速冪的題都要取模，一般保證模的平方不超數的範圍， 因為指數太大的數是會爆掉int 和long long 的

#include<iostream>

using namespace std;

#define mod 1000000007

long long quick_pow(int n,int base)

//n是指數 base是底 即計算的是base^n 當然結果是取模了的

{

    long long ans=1;//預設ans大於等於1因為不能算負指數

    long long multi=base;

    while(n)

    {

        if(n%2) ans*=multi;

        ans%=mod;//由於數太大一般要取模

        n/=2;

        multi*=multi;

        multi%=mod;

    }

return ans;

}

int main()

{

    int n,base;

    while(cin>>n>>base)

        cout<<quick_pow(n,base)<<endl;

    return 0;

}

可能你會問了這個演算法有什麼用呢？其實用的更多是使用矩陣快速冪，算遞推式，注意是遞推式，簡單的如斐波那契數列的第一億項的結果模上10000000後是多少你還能用遞推式去，逐項遞推嗎？當然不能，這裡就可以發揮矩陣快速冪的神威了，那斐波那契數列和矩陣快速冪能有一毛錢的關係？答案是有而且很大

斐波那契的定義是f(1)=f(2)=1; 然後f(n)=f(n-1)+f(n-2) (n>=2) 我們也可以這樣定義f(1)=f(2)=1; [f(n),f(n-1)]=[f(n-1),f(n-2)][1,1,1,0],其中[1,1,1,0] 是一個2*2的矩陣 上面一行是1,1，下面一行是1,0，這樣就可以化簡了寫成[f(n),f(n-1)]=[f(2),f(1)]*[1,1,1,0]^(n-2)

化簡一下

這樣就可以用矩陣快速冪，快速的推出斐波那契數列的第一億項的值了（當然是取模的值了）是不是很神奇，類似的遞推式也可以，化成這種形式，用矩陣快速冪進行計算

下面附一個矩陣快速冪的程式碼，當然所有矩陣都是要模的

# include<cstdio>

# include<cstring>

using namespace std;

#define NUM 50

int MAXN,n,mod;

struct Matrix//矩陣的類

{

int a[NUM][NUM];

void init()           //將其初始化為單位矩陣

{

memset(a,0,sizeof(a));

for(int i=0;i<MAXN;i++)

a[i][i]=1;

}

} A;

Matrix mul(Matrix a,Matrix b)  //(a*b)%mod  矩陣乘法

{

Matrix ans;

for(int i=0;i<MAXN;i++)

for(int j=0;j<MAXN;j++)

{

ans.a[i][j]=0;

for(int k=0;k<MAXN;k++)

ans.a[i][j]+=a.a[i][k]*b.a[k][j];

ans.a[i][j]%=mod;

}

return ans;

}

Matrix add(Matrix a,Matrix b)  //(a+b)%mod  //矩陣加法

{

int i,j,k;

Matrix ans;

for(i=0;i<MAXN;i++)

for(j=0;j<MAXN;j++)

{

ans.a[i][j]=a.a[i][j]+b.a[i][j];

ans.a[i][j]%=mod;

}

return ans;

}

Matrix pow(Matrix a,int n)    //(a^n)%mod  //矩陣快速冪

{

Matrix ans;

ans.init();

while(n)

{

if(n%2)//n&1

ans=mul(ans,a);

n/=2;

a=mul(a,a);

}

return ans;

}

Matrix sum(Matrix a,int n)  //(a+a^2+a^3....+a^n)%mod// 矩陣的冪和

{

int m;

Matrix ans,pre;

if(n==1)

return a;

m=n/2;

pre=sum(a,m);                      //[1,n/2]

ans=add(pre,mul(pre,pow(a,m)));   //ans=[1,n/2]+a^(n/2)*[1,n/2]

if(n&1)

ans=add(ans,pow(a,n));          //ans=ans+a^n

return ans;

}

void output(Matrix a)//輸出矩陣

{

for(int i=0;i<MAXN;i++)

for(int j=0;j<MAXN;j++)

printf("%d%c",a.a[i][j],j==MAXN-1?'\n':' ');

}

int main()

{

freopen("in.txt","r",stdin);

Matrix ans;

scanf("%d%d%d",&MAXN,&n,&mod);

for(int i=0;i<MAXN;i++)

for(int j=0;j<MAXN;j++)

{

scanf("%d",&A.a[i][j]);

A.a[i][j]%=mod;

}

ans=sum(A,n);

output(ans);

return 0;

}

快速冪乘：