快速冪演算法及其擴充

快速冪演算法

問題引入：求an(a,n∈N+) 
樸素演算法：令ans初始值為1，乘n次a得到an 
樸素演算法時間複雜度：O(n)

問題：如果n非常大，比如高達1015，怎麼辦？ 
思考：樸素演算法哪裡可以優化？

樸素演算法的特點是，連乘過程中底數始終為a，這很不聰明。考慮下例： 
a=2,n=15 
我們沒有必要乘15次2，注意到15=23+22+21+20，不妨將215拆分為223×222×221×220 
然後看一下n的二進位制形式：15=(1111)2.為什麼要看這個呢？再舉一個例子。

a=2,n=6 
26=222×221，而6=(110)2

觀察發現，當n的第i位（從低到高，i>=0）為1時，an的拆解表示式裡就要乘上一項a2i.

主要思想：對於a^b,當b很大時，容易超時，如果b&1=1，結果要累乘，每次迴圈b要右移1位，（底數a每次要累乘）例如當b=13時（1011）

於是我們有了樸素演算法改進的思路：逐位判斷冪次n的二進位制位是否為1，若是，給答案乘上一個a2i. 
改進後的演算法的虛擬碼描述如下：

public class MyPower {  
      
    public static void main(String[] args){  
           
        int a = 3;  
        int b = 13;  
        int m = power(a,b);  
        System.out.println(m);  
    }  
  
    private static int power(int a, int b) {  
        int result = 1;  
        while(b > 0){  
              
            if((b & 1) == 1)//b的某位上為1時才累乘  
                result *= a;  
              
            a *= a;//數學公式所得  
            b >>= 1;//右移1位  
        }  
        return result;  
    }  
}  

a, n, ans;
ans=1;
while n>0
    a = a*a;
    if n%2
        ans = ans*a;
    n = n/2;

• 1
• 2
• 3
• 4
• 5
• 6
• 7

• 1
• 2
• 3
• 4
• 5
• 6
• 7

這個演算法的時間複雜度是多少呢？很顯然，它取決於n的二進位制形式有多少位，因此T(n)=⌈log2(n)⌉=O(logn).這就是快速冪演算法。

其實快速冪演算法還可以遞迴實現，因為： 
當n為偶數時，an=(an/2)2 
當n為奇數時，an=(an/2)2×a 
邊界條件：當n=1，答案為a 
C語言描述如下：

int quick_power(int a, int n)
{
    if(n == 1) return a;
    int x = quick_power(a, n/2);
    long long ans = (long long)x*x;
    if(n%2) ans *= a;
    return (int)ans;
}

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8

程式碼中加入了防溢位處理，用快速冪演算法的時候比較容易犯的一個錯誤就是忘了考慮溢位，因此在使用的時候要看清楚資料範圍，估算一下答案上界。另外，快速冪演算法不推薦用遞迴實現，因為非遞迴版本不但程式碼也很簡潔，而且效率還更優。

擴充一：快速冪模M演算法

有時候所求冪的結果可能很大，於是問題要求對結果模上一個數M。我們只需要在原來演算法的基礎上運用一下模運算的性質即可。所謂模運算性質是指以下兩條： 
∑ni=1ai mod M=(∑ni=1ai mod M)mod M 
∏ni=1ai mod M=(∏ni=1ai mod M)mod M 
演算法非遞迴實現的虛擬碼描述為  

int runFermatPower(int a, int b, int m){  
          
        int result = 1;  
        while (b > 0) {  
            if ((b & 1) == 1)  
                result = (result * a);  
  
            a = (a * a) % m;  
            b >>= 1;  
        }  
        return result % m;  
}  
解釋：當b=1101（2）時，從第1位開始result累乘，a = （a*a）%m加上迴圈可以看成表示式（a%m)^2 % m....因為(a%m)^2 % m = a^2 % m ,所以我們可以把a^13%m看成((a^1) % m) ,((a^4) % m) ,((a^8) % m)的累乘，最後對m取模。


但累乘很容易造成溢位，所以我們可以把程式碼改成如下形式：

int runFermatPower(int a, int b, int m){  
          
        int result = 1;  
        while (b > 0) {  
            if ((b & 1) == 1)  
                result = (result * a) % m;  
  
            a = (a * a) % m;  
            b >>= 1;  
        }  
        return result;  
}  

//在累乘過程中 res=(res*a)%m就取模，避免溢位

a, n, ans, M;
ans=1;
while n>0
    a = a*a % M;
    if n%2
        ans = ans*a % M;
    n = n/2;

• 1
• 2
• 3
• 4
• 5
• 6
• 7

• 1
• 2
• 3
• 4
• 5
• 6
• 7

擴充二：矩陣快速冪演算法

快速冪演算法解決的是整數求冪的問題，而矩陣快速冪解決的是矩陣求冪問題，兩者沒有本質的區別。如果用C++實現，我們只要定義一個矩陣類，然後過載一下乘法運算子，原先的快速冪演算法幾乎不需要改變。 
矩陣快速冪常常用於線性遞推式的加速。以下僅舉一例。

快速求斐波那契數列第n項

對於這個問題，普通求法的複雜度是O(n)，現在我們用矩陣快速冪將它優化到O(logn). 
首先，將遞推式f(n)=f(n−1)+f(n−2)改寫成矩陣形式 

進而得到 

接下來，用矩陣加速演算法求出[1110]n−4，再做一次矩陣乘法，所得矩陣的(0,0)元素就是最終結果。

對於更一般的線性遞推式，構造其加速矩陣的方式超出了本文的論述範圍，在此省略。