快速冪取餘，冪取餘

快速冪取餘，冪取餘

求a^b mod c

 演算法1.

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

int  ans=1, i;for(i=1;i<=b;i++)    ans*=a;ans%=c;

這個演算法的時間複雜度體現在for迴圈中，為O（b）.

這個演算法存在著明顯的問題,如果a和b過大，很容易就會溢出。

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

a^b mod c=(a mod c)^b

引理：

(a *b) mod c =[ ( a mod c )* (b mod c) ] mod c ;

證明： 設 a mod c =d，b mod c= e;

       則：a=t*c + d ;  b=k*c + e ;

       (a*b)mod c = (t*c+d)(t*c+e)

                 = (tk c^2 + ( te+dk ) *c + d*e) mod c

                 =de mod c

即積的取餘等於取餘的積的取餘.

(a ^ b)mod c 由上述公式迭代即可得到 ( a mod c)^b.

證明了以上的公式以後，我們可以先讓a關於c取餘，這樣可以大大減少a的大小，於是不用思考的進行了改進：

 演算法2：

int ans = 1 , i ;   a = a % c; //加上這一句  for ( i = 1;i<=b;i++)      ans = ans * a;   ans = ans % c;  

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

演算法3：

int  ans = 1 ,i ;  a = a % c;for(int i = 1;i<=b;i++)     ans = (ans * a) % c; //這裡再取了一次餘 ans = ans % c; 

 

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

快速冪取餘依賴於以下公式：

 

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

演算法4： 

int  ans = 1 ，i ; a = a % c; if (b%2==1)  ans = (ans * a) mod c; //如果是奇數，要多求一步，                     可以提前算到 ans 中。k = (a*a) % c;  //我們取a^2 而不是a for( 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

而不是原來的a^b 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;  }  

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

long long  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;       //   b>>=1;        a = (a * a) % c;     }     return ans; }  

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

 

以下內容僅供參考： 

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

a^b%c 求解這個問題，我們也可以從二進位轉換來考慮： 

將10進位的b轉化成2

進位的運算式：

 

注意此處的an要麼為0，要麼為1，如果為0，那麼這一項就是1，這個對應了上面演算法過程中b是偶數的情況；

為1對應了b是奇數的情況.

pascal 快速冪：取餘

{輸入b,p,k的值,求b^p mod k的值。其中b,p,k*k為長整型數。
範例輸入：2 10 1000，輸出：24}
var b,p,k,t,ans:longint;
begin
readln(b,p,k);
ans:=1; t:=b;
while p>0 do
begin
if odd(p) then ans:=ans*t mod k;
p:=p div 2;
t:=t*t mod k;
end;
writeln(ans);
end.
 
下面是一個矩陣快速冪解斐波那契數列第n項對第m項取餘的程式，但是總是運行出錯，解

按照正常的邏輯是只要求a[2][2]={1,1,1,0}這個矩陣的n次方就可以得到斐波那契數列的第n項（即a[0][1]）的值。但是你忽略了一點，就是你在求a[0][1]，a[1][0]，a[1][1]的值的時候你的a[0][0]的值其實已經改變了，導致你求得的a[0][1]的值不正確，繼而影響到a[1][0]和a[1][1]的值。
因此，我們在遍曆這四個值的時候是不能改變其中任何一個的，只有改完之後才能去改變它的值。所以我們可以用幾個變數先將求得的新矩陣各值存起來，如下：
#include<stdio.h>
int a[2][2]={1,1,1,0},b[2][2]={1,1,1,0};//使用兩個二維數組表示快速冪演算法要用到的矩陣
int main()
{
int n,m,i,j,t,u;
int a1,a2,a3,a4;
while(scanf("%d %d",&n,&m)!=EOF)
{
if(m==-1 && n==-1)//當輸入的m.n值為-1時結束整個演算法
return 0;
/*以下是對於第n項斐波那契數的計算*/
if(n==0)
a[0][0]=0;
else if(n==1 || n==2)
a[0][0]=1;
else
{
for(i=3;i<=n;i++)
{
a1=a[0][0]*b[0][0]+a[0][1]*b[1][0];
a2=a[0][0]*b[0][1]+a[0][1]*b[1][1];
a3=a[1][0]*b[0][0]+a[1][1]*b[1][0];
a4=a[1][0]*b[0][1]+a[1][1]*b[1][1];
a[0][0]=a1;
a[0][1]=a2;
a[1][0]=a3;
a[1][1]=a4;
}
}
t=a[0][0];
a[0][0]=1;//重設矩陣
a[0][1]=1;
a[1][0]=1;
a[1][1]=0;
if(m==0)
a[0][0]=0;
else if(m==1 || m==2)
a[0][0]=1;
else
{
for(j=3;j<=m;j++)
{
a1=a[0][0]*b[0][0]+a[0][1]*b[1][0];
a2=a[0][0]*b[0][1]+a[0][1]*b[1][1];
a3=a[1][0]*b[0][0]+a[1][1]*b[1][0];
a4=a[1][0]*b[0][1]+a[1][1]*b[1][1];
a[0][0]=a1;
a[0][1]=a2;
a[1][0]=a3;
a[1][1]=a4;
}
}
u=a[0][0];
a[0][0]=1;//重設矩陣
a[0][1]=1;
a[1][0]=1;
a[1][1]=0;
t=t%u;
printf("%d\n",t);
}
return 0;
}

還有一點，就是你的矩陣相乘後兩項寫錯了，自己對比一下改過來吧。

這樣就能得到你想要的結果了。
...餘下全文>>