模冪

模冪（modular exponentiation）是一種對模進行的冪運算，在電腦科學，尤其是公開密鑰加密方面有一定用途。

模冪運算是指求整數$b$的$e$次方$b^e$被正整數$m$所除得到的餘數$c$的過程，可用數學符號表示為$c=b^e \bmod m$。由$c$的定義可得$0\leq c <m$。

例如，給定$b=5$，$e=3$和$m=13$，$5^3=125$被13除得的餘數$c=8$。

指數$e$為負數時可使用擴展歐幾里得演算法找到$b$模除$m$的模逆元$d$來執行模冪運算，即：

$c=b^e \bmod m=d^{-e} \bmod m$，其中 $e<0$且$b\cdot d \equiv 1 \pmod{m}$。

即使在整數很大的情況下，上述模冪運算其實也是易於執行的。然而，計算模的離散對數（即在已知$b$，$c$和$m$時求出指數$e$）則較為困難。這種類似單向函數的表現使模冪運算可用於加密演算法。

## 直接演算法

計算模冪的最直接方法是直接算出$b^e$，再把結果模除$m$。假設已知$b=4$，$e=13$，以及$m=497$，要求$c$：

$$c\equiv 4^{13}\pmod{497}$$

可用計算器算得413結果為67,108,864，模除497，可得c等於445。

注意到$b$只有一位，$e$也只有兩位，但$b^e$的值卻有8位。

強加密時$b$通常至少有1024位。考慮$b=5\times 10^{76}$和$e=17$的情況，$b$的長度為77位，$e$的長度為2位，但是$b^e$的值以十進位表示卻已經有1304位。現代電腦雖然可以進行這種計算，但計算速度卻大大降低。

用這種演算法求模冪所需的時間取決於操作環境和處理器，時間複雜度為$\mathrm{O}(e)$。

## 空間優化

這種方法比第一種所需要的步驟更多，但所需記憶體空間和時間均大為減少，其原理為：

給定兩個整數$a$和$b$，以下兩個等式是等價的：

$$c \bmod m=(a\cdot b) \bmod m$$

$$c \bmod m=[(a \bmod m)\cdot (b \bmod m)]\bmod m$$

演算法如下：

- 令$c=1$，$e'=0$。
- $e'$自增1。
- 令$c=(b\cdot c) \bmod m$.
- 若$e'<e$，則返回第二步；否則，$c$即為$c\equiv b^e \pmod{m}$。

再以$b=4$，$e=13$，$m=497$為例說明，演算法第三步需要執行13次：

- $$e' = 1\cdot c=(1\cdot 4) \bmod 497=4 \bmod 497=4$$
- $$e' = 2\cdot c=(4\cdot 4) \bmod 497=16 \bmod 497=16$$
- $$e' = 3\cdot c=(16\cdot 4) \bmod 497=64 \bmod 497=64$$
- $$e' = 4\cdot c=(64\cdot 4) \bmod 497=256 \bmod 497=256$$
- $$e' = 5\cdot c=(256\cdot 4) \bmod 497=1024 \bmod 497=30$$
- $$e' = 6\cdot c=(30\cdot 4) \bmod 497=120 \bmod 497=120$$
- $$e' = 7\cdot c=(120\cdot 4) \bmod 497=480 \bmod 497=480$$
- $$e' = 8\cdot c=(480\cdot 4) \bmod 497=1920 \bmod 497=429$$
- $$e' = 9\cdot c=(429\cdot 4) \bmod 497=1716 \bmod 497=225$$
- $$e' = 10\cdot c=(225\cdot 4) \bmod 497=900 \bmod 497=403$$
- $$e' = 11\cdot c=(403\cdot 4) \bmod 497=1612 \bmod 497=121$$
- $$e' = 12\cdot c=(121\cdot 4) \bmod 497=484 \bmod 497=484$$
- $$e' = 13\cdot c=(484\cdot 4) \bmod 497=1936 \bmod 497=445$$

因此最終結果$c$為445，與第一種方法所求結果相等。

與第一種方法相同，這種演算法需要$\mathrm{O}(e)$的時間才能完成。但是，由於在計算過程中處理的數字比第一種演算法小得多，因此計算時間至少減少了$\mathrm{O}(e)$倍。

演算法偽代碼如下：

```
'''function''' modular_pow(b, e, m)
    '''if''' m = 1
        '''then''' '''return''' 0
    c := 1
    '''for''' e' = 0 '''to''' e-1
        c := (c * b) '''mod''' m
    '''return''' c
```

## 從右到左二位演算法

第三種方法結合了第二種演算法和平方求冪原理，使所需步驟大大減少，同時也與第二種方法一樣減少了記憶體占用量。

首先把$e$表示成二進位，即：

$$e = \sum_{i=0}^{n-1} a_i 2^i$$

此時$e$的長度為$n$位。對任意$i$（$0\leq i < n$），$a_i$可取0或1任一值。由定義有$a_{n-1}=1$。

$b^e$的值可寫作：

$$b^e = b^{\left( \sum_{i=0}^{n-1} a_i 2^i \right)} = \prod_{i=0}^{n-1} \left( b^{2^i} \right) ^ {a_i}$$

因此答案$c$即為：

$$c \equiv \prod_{i=0}^{n-1} \left( b^{2^i} \right) ^ {a_i}\ (\mbox{mod}\ m)$$

### 偽代碼

下述偽代碼基於布魯斯·施奈爾所著《應用密碼學》。其中base，exponent和modulus分別對應上式中的$b$，$e$和$m$。

```
'''function''' modular_pow(base, exponent, modulus)
    '''if''' modulus = 1 '''then''' '''return''' 0
    [[斷言 (程式)|Assert]] :: (modulus - 1) * (modulus - 1) does not overflow base
    result := 1
    base := base '''mod''' modulus
    '''while''' exponent > 0
        '''if''' (exponent '''mod''' 2 == 1):
           result := (result * base) '''mod''' modulus
        exponent := exponent >> 1
        base := (base * base) '''mod''' modulus
    '''return''' result
```

注意到在首次進入循環時，變量base等於$b$。在第三行代碼中重複執行平方運算，會確保在每次循環結束時，變量base等於$b^{2^i} \bmod m$，其中$i$是循環執行次數。

本例中，底數$b$的指數為$e = 13$。

指數用二進位表示為1101，有4位，故循環執行4次。

4位數字從右到左依次為1，0，1，1。

首先，初始化結果$R$為1，並將b的值保存在變量$x$中：

$R \leftarrow 1 \, ( = b^0) \text{且 } x \leftarrow b$.

第1步 第1位為1，故令$R \leftarrow R \cdot x \text{ }(= b^1)$。

令$x \leftarrow x^2 \text{ }(= b^2)$。

第2步 第2位為0，故不給R賦值。

令$x \leftarrow x^2 \text{ }(= b^4)$。

第3步 第3位為1，故令$R \leftarrow R \cdot x \text{ }(= b^5)$。

令$x \leftarrow x^2 \text{ }(= b^8)$。

第4步 第4位為1，故令$R \leftarrow R \cdot x \text{ }(= b^{13})$。

這是最後一步，所以不需要對x求平方。

綜上，$R$為$b^{13}$。

以下計算$b = 4$的$e = 13$次方對497求模的結果。

初始化：

$R \leftarrow 1 \, ( = b^0)$且$x \leftarrow b = 4$。

第1步 第1位為1，故令$R \leftarrow R \cdot 4 \pmod {497} \equiv 4 \pmod{497}$。

令$x \leftarrow x^2 \text{ }(= b^2) \equiv 4^2 \equiv 16 \pmod{497}$。

第2步 第2位為0，故不給R賦值。

令$x \leftarrow x^2 \text{ }(= b^4) \equiv 16^2 \pmod{497} \equiv 256 \pmod{497}$。

第3步 第3位為1，故令$R \leftarrow R \cdot x \text{ }(= b^5) \equiv 4 \cdot 256 \pmod{497} \equiv 30 \pmod{497}$。

令$x \leftarrow x^2 \text{ }(= b^8) \equiv 256^2 \pmod{497} \equiv 429 \pmod{497}$。

第4步 第4位為1，故令$R \leftarrow R \cdot x \text{ }(= b^{13}) \equiv 30 \cdot 429 \pmod{497} \equiv 445 \pmod{497}$。

綜上，$R$為$4^{13} \pmod{497} \equiv 445 \pmod{497}$，與先前演算法中所得結果相同。

該演算法時間複雜度為O(logexponent。指數exponent值較大時，這種演算法與前兩種O(exponent演算法相比具有明顯的速度優勢。例如，如果指數為220 = 1048576，此演算法只需執行20步，而非1,048,576步。

### Lua實現

```
'''function''' modPow(b,e,m)
        '''if''' m == 1 '''then'''
                '''return''' 0
        '''else'''
                '''local''' r = 1
                b = b % m
                '''while''' e > 0 '''do'''
                        '''if''' e % 2 == 1 '''then'''
                                r = (r*b) % m
                        '''end'''
                        e = e >> 1     --Lua 5.2或更早版本使用e = math.floor(e / 2)
                        b = (b^2) % m
                '''end'''
                '''return''' r
        '''end'''
'''end'''
```

## 軟體實現

鑑於模冪運算是電腦科學中的重要操作，並且已有高效演算法，所以許多程式語言和高精度整數庫都有執行模冪運算的函數：

- Python內置的`pow()`（求冪）函數
- .NET框架的`BigInteger`類的`[http://msdn.microsoft.com/en-us/library/system.numerics.biginteger.modpow%28v=vs.100%29.aspx#pow ModPow()]`方法
- Java的`java.math.BigInteger`類的java/math方法
- Perl的`Math::BigInt`模塊的`bmodpow()`方法
- Raku內置的`expmod`例程
- Go的`big.Int`類的`Exp()`（求冪）方法
- PHP的BC Math庫的`bcpowmod()`函數
- GNU多重精度運算庫（GMP）的`mpz_powm()`函數
- FileMaker Pro的`[http://www.briandunning.com/cf/1482 @PowerMod()]`函數（以1024位RSA加密為例）
- Ruby的`openssl`包的`OpenSSL::BN#mod_exp`方法
- HP Prime計算器的CAS.powmod函數

## 另見

- 蒙哥馬利算法，用於計算模很大時的餘數。

## 外部連結

- Paul Garrett, Fast Modular Exponentiation Java Applet