轉載原地址 http://blog.csdn.net/hikean/article/details/9749391
快速冪或者矩陣快速冪在算指數時是很高效的，他的基本原理是二進制，下面的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==) //亦可以寫成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;
}