平方求幂

在數學和程式設計中，**平方求冪**（英語：exponentiating by squaring）或**快速冪**是快速計算一個數（或更一般地說，一個半群的元素，如多項式或方陣）的大正整數乘冪的一般方法。這些演算法可以非常通用，例如用在模算數或矩陣冪。對於通常使用加性表示法的半群，如密碼學中使用的橢圓曲線，這種方法也稱為 **double-and-add**。

## 基本方法

該方法是基於觀察到，對於正整數 $n$ ，可知

　
: $x^{n}={\begin{cases}x\,(x^{2})^{\frac {n-1}{2}},&{\mbox{if }}n{\mbox{ is odd}}\\(x^{2})^{\frac {n}{2}},&{\mbox{if }}n{\mbox{ is even}}.\end{cases}}$

該方法使用指數的位（二進位的位，即 bit，下文稱為「位」）來確定計算哪些冪。

此例顯示如何使用此方法計算 $x^{13}$ 。
冪指數 13 的二進位為 1101。這些位按照從左到右的順序使用。
指數有 4 位元，所以有 4 次迭代。

首先，將結果初始化為 1： $r\leftarrow 1\,(=x^{0})$

 1. $r\leftarrow r^{2}\,(=x^{0})$ ，第 1 位 = 1，所以計算 $r\leftarrow r\cdot x\,(=x^{1})$ 。
2. $r\leftarrow r^{2}\,(=x^{2})$ ，第 2 位 = 1，所以計算 $r\leftarrow r\cdot x\,(=x^{3})$ 。
3. $r\leftarrow r^{2}\,(=x^{6})$ ，第 3 位 = 0，所以這一步什麼都不做。
4. $r\leftarrow r^{2}\,(=x^{12})$ ，第 4 位元 = 1，所以計算 $r\leftarrow r\cdot x\,(=x^{13})$ 。

這可以按照下面的遞迴演算法來實現：

儘管不是尾呼叫，但是通過引入輔助函式，該演算法可以被重寫成尾遞迴演算法：

該演算法的迭代版本的輔助空間是有界的，代碼如下

c++ 實現（非遞迴）
返回 $p^{n}$ 對 $Mod$ 求模

## 計算複雜度

簡要分析表明此演算法用了 $\lfloor \log _{2}n\rfloor $ 次平方，以及至多 $\lfloor \log _{2}n\rfloor $ 次乘法，其中 $\lfloor \;\rfloor $ 表示向下取整函式。更確切地說，做乘法的次數比 $n$ 的二進位展開的次數要少一次。對於 $n$ 大於 4 左右的時候，這種演算法在計算上就已經比_天真地_將它與自身重複地相乘更高效了。

每次平方的結果大約是前一次結果的兩倍，因此，如果兩個 $d$ 位數的相乘的實現要進行 $\mathrm {O} (d^{k})$ 次計算（其中 $k$ 為一固定值），那麼計算 $x^{n}$ 的複雜度為：

　
: $\sum \limits _{i=0}^{O(\log(n))}(2^{i}O(\log(x)))^{k}=O((n\log(x))^{k})$

##  $2^{k}$ 法

此演算法先把指數展開成 $2^{k}$ 形式，然後再計算 $x^{n}$ 的值。它在 1939 年由 Brauer 首次提出。在下面的演算法中，使用以下函式 $f(0)=(k,0)$ 和 $f(m)=(s,u)$ ，其中 $m=u\cdot 2^{s}$ ， $u$ 為奇數。

演算法：

**輸入**
: G 的一個元素 $x$ ，參數 $k>0$ ，一個非零整數 $n=(n_{l-1},n_{l-2},\ldots n_{0})_{2}{^{k}}$ 以及預計算的值 $x^{3},x^{5},...,x^{2^{k}-1}$ 。

**輸出**
: _G_ 中的元素 $x^{n}$

```
y := 1; i := l-1
**while** i>=0 do
    (s,u) := f(ni)
    **for** j:=1 **to** k-s **do**
        y := y2 
    y := y*xu
    **for** j:=1 **to** s **do**
        y := y2
    i := i-1
**return** y
```

為了獲得最佳效率， _$k$_ 應該是滿足

　
: $\log(n)<{\frac {k(k+1)\cdot 2^{2k}}{2^{k+1}-k-2}}+1$

的最小整數。

## 滑動窗口法

此方法是 $2^{k}$ 法的更高效的變體。例如，要計算 398 次冪，二進位展開為 (110 001 110)2，採用長度為 3 的窗，使用 $2^{k}$ 法，需要計算 1, $x^{3},x^{6},x^{12},x^{24},x^{48},x^{49},x^{98},x^{99},x^{198},x^{199},x^{398}$ 。
但也可以計算 1, $x^{3},x^{6},x^{12},x^{24},x^{48},x^{49},x^{96},x^{192},x^{198},x^{199},x^{398}$ ，這就會省下一次乘法，相當於是計算 (110 001 110) n2 的值

以下是一般演算法：

演算法：

**輸入**
: _G_ 的元素 _$x$_ ，非負整數 $n=(n_{l-1},n_{l-2},\ldots n_{0})_{2}$ ，一個參數 $k>0$ ，以及預計算的值 $x^{3},x^{5},...,x^{2^{k}-1}$ 。

**輸出**
: 元素 $x^{n}\in G$

演算法：

```
y := 1; i := l-1
**while** i > -1 **do**
    **if** ni=0 **then**
        y:=y2' i:=i-1
    **else**
        s:=max{i-k+1,0}
        **while** ns=0 **do**
            s:=s+1 
        **for** h:=1 **to** i-s+1 **do**
            y:=y2
        u:=(ni,ni-1,....,ns)2
        y:=y*xu
        i:=s-1
**return** y
```

## 蒙哥馬利階梯法

求冪的許多演算法都不提供對旁路攻擊的防護。也就是說，監測到乘方運算的攻擊者可以（部分）還原所計算的指數。就如眾多公鑰加密系統中那樣，如果指數需要保密的話，這就是個問題了。一個叫做蒙哥馬利階梯的方法解決了這個問題。

給定一個非零正整數的二進位展開 $n=(n_{k-1}\ldots n_{0})_{2}$ （其中 $n_{k-1}=1$ ），可以以下面方式計算 $x^{n}$ ：

```
x1=x; x2=x2
for i=k-2 to 0 do
  If ni=0 then
    x2=x1*x2; x1=x12
  else
    x1=x1*x2; x2=x22
return x1
```

該演算法會執行一系列固定的操作（複雜度 $\log n$ ）：無論每一位的具體值如何，指數中的每一位都會進行乘法和平方。

蒙哥馬利階梯法的這種具體實現還無法抵禦快取時序攻擊：當根據秘密指數的位值存取不同的變數時，記憶體存取延遲仍可能被攻擊者觀察到。

## 固定基底的冪

當基底固定而指數變化時，可以使用幾種方法來計算 $x^{n}$ 。可以看出，預計算在這些演算法中起著關鍵作用。

### 姚期智的方法

姚期智的方法與 $2^{k}$ 法不同，是把指數以 $b=2^{k}$ 為基底展開，並按上面的演算法進行計算。令 $n$ , $n_{i}$ , $b$ , $b_{i}$ 為整數。

把指數 $n$ 寫成

　
: $n=\sum _{i=0}^{w-1}n_{i}b_{i}$ 其中對所有 $i\in [0,w-1]$ 都有 $0\leqslant n_{i}<h$

令 $x_{i}=x^{b_{i}}$ 。該演算法使用等式

　
: $x^{n}=\prod _{i=0}^{w-1}{x_{i}}^{n_{i}}=\prod _{j=1}^{h-1}{{\bigg [}\prod _{n_{i}=j}x_{i}{\bigg ]}}^{j}$

給定 $G$ 的元素 $x$ ，指數 $n$ 寫成上述形式，並且預先計算 $x^{b_{0}}\ldots x^{b_{w-1}}$ 的值，元素 $x^{n}$ 就可以用下面的演算法計算了。

```
 y=1,u=1 and j=h-1
 while j > 0 do
   for i=0 to w-1 do
     if ni=j then u=u*xbi
   y=y*u
   j=j-1
 return y
```

如果令 $h=2^{k}$ ， $b_{i}=h^{i}$ ，那麼這些 $n_{i}$ 就是 $n$ 以 $h$ 為基的每一位。姚期智的方法是把之前的那些 $x_{i}$ 收集到 $u$ 中，which appear to the highest power $h-1$ ；in the next round those with power $h-2$ are collected in $u$ as well etc. 變數 $y$ 被原始的 $u$ 乘了 $h-1$ 次，內第二高的指數乘了 $h-2$ 次……
該演算法計算 $x^{n}$ 要用 $w+h-2$ 次乘法，儲存 $w+1$ 個元素。

### 歐幾里得法

歐幾里德法首先在 P.D Rooij 的《使用預計算和向量加法鏈的高效求冪》（Efficient exponentiation using precomputation and vector addition chains）中介紹。

這種計算群 **G** 中 $x^{n}$ （ $n$ 為自然數）的方法是，遞迴地使用下面的等式：

　
: ${x_{0}}^{n_{0}}\cdot {x_{1}}^{n_{1}}={\left(x_{0}\cdot {x_{1}}^{q}\right)}^{n_{0}}\cdot {x_{1}}^{n_{1}\mod {n_{0}}}$ , where $q=\left\lfloor {\frac {n_{1}}{n_{0}}}\right\rfloor $

　
: （換句話說，用指數 $n_{1}$ 與 $n_{0}$ 的歐幾里得除法來返回商 $q$ 和餘數 $n_{1}{\bmod {n}}_{0}$ ）。

給定群 **G** 中的基底元素 $x$ ，把指數 $n$ 用姚期智的方法寫出來，然後就可以用預計算的 $l$ 個值 $x^{b_{0}},...,x^{b_{l_{i}}}$ 計算 $x^{n}$ 了。

```
    **Begin loop**   
        Find $M\in \left[0,l-1\right]$, such that $\forall i\in \left[0,l-1\right],{n_{M}}\geq {n_{i}}$;
        Find $N\in \left(\left[0,l-1\right]-M\right)$, such that $\forall i\in \left(\left[0,l-1\right]-M\right),{n_{N}}\geq {n_{i}}$;
        **Break loop** if ${n_{N}}=0$;
        **Let** $q=\left\lfloor {n_{M}}/{n_{N}}\right\rfloor $, and then **let** ${n_{N}}=\left({n_{M}}{\bmod {n_{N}}}\right)$;
        Compute recursively ${x_{M}}^{q}$, and then **let** ${x_{N}}={x_{N}}\cdot {x_{M}}^{q}$;
    **End loop**;
    **Return** $x^{n}={x_{M}}^{n_{M}}$.
```

該演算法首先在 $n_{i}$ 中找到最大值，在找到集合 $\{n_{i}\backslash i\neq M\}$ 中的最大值。
然後遞迴求 $x_{M}$ 的 $q$ 次冪，把這個值乘以 $x_{N}$ ，賦值給 $x_{N}$ ；把 $n_{M}$ 模 $n_{N}$ 的結果賦值給 $n_{M}$ 。

## 更多應用

這一思路還可用於計算大指數冪除以一個數的餘數。特別是在密碼學中，在整數除以 _q_ 的餘數的環中計算冪是很有用處的。在群中計算整數冪也是很有用的，使用法則

　
: Power(_x_, −_n_) = (Power(_x_, _n_))−1.

該方法對所有半群都適用，通常用於計算矩陣的冪。

例如，如果使用樸素的方法，計算

　
: 13789722341 (mod 2345)

將花費很長時間以及大量的儲存空間：計算 13789722341，並除以 2345 求余。即便使用更有效的方法也需要很長時間：求 13789 的平方，除以 2345 求余，將結果乘以 13789，再求余，一直往下計算。這會需要 $2\log _{2}(722340)\leq 40$ 次模乘。

把「\*」理解為 _x_\*_y_ = _xy_ mod 2345（即相乘後取模）後，應用上面的平方求冪演算法，只需要 27 次乘法和整數除法（所有這些都可以儲存在單個機器字中）就可以完成計算了。

## 範例實現

### 通過 2 的冪進行計算

這是用 Ruby 寫的上述演算法的非遞迴實現。

由於低階語言會將 n=n/2 隱式向 0 取整，n=n-1 對那些語言來說，就是冗餘的步驟了。
n [0] 是 n 的二進位表示的最右邊的位，所以如果它是 1，則該數是奇數，如果它是零，則該數是偶數。它也是以 2 為模 n 的餘數。

#### 執行實例：計算 310

```
parameter x =  3
parameter n = 10
result := 1

**Iteration 1**
  n = 10 -> n is even
  x := x2 = 32 = 9
  n := n / 2 = 5

**Iteration 2**
  n = 5 -> n is odd
      -> result := result * x = 1 * x = 1 * 32 = 9
         n := n - 1 = 4
  x := x2 = 92 = 34 = 81
  n := n / 2 = 2

**Iteration 3**
  n = 2 -> n is even
  x := x2 = 812 = 38 = 6561
  n := n / 2 = 1

**Iteration 4**
  n = 1 -> n is odd
      -> result := result * x = 32 * 38 = 310 = 9 * 6561 = 59049
         n := n - 1 = 0

return result
```

#### 執行實例：計算 310

```
result := 3
bin := "1010"

**Iteration for digit 2:**
  result := result2 = 32 = 9
  1**0**10bin - Digit equals "0"

**Iteration for digit 3:**
  result := result2 = (32)2 = 34  = 81
  10**1**0bin - Digit equals "1" --> result := result*3 = (32)2*3 = 35  = 243

**Iteration for digit 4:**
  result := result2 = ((32)2*3)2 = 310  = 59049
  101**0**bin - Digit equals "0"

return result
```

JavaScript-Demonstration: http://home.mnet-online.de/wzwz.de/temp/ebs/en.htm （頁面存檔備份，存於網際網路檔案館）

### 冪的乘積的計算

平方求冪也可用於計算 2 個或多個冪的乘積。
如果基礎群或半群是可交換的，那麼常常可以通過同時計算乘積來減少乘法的次數。

#### 例子

式子 a7×b5 可以分三步計算：

　
: ((a)2×a)2×a （計算 a7 需要四次乘法）

　
: ((b)2)2×b   （計算 b5 需要三次乘法）

　
: (a7)×(b5) （計算二者乘積需要一次乘法）

所以總共需要八次乘法。

更快的解法是同時計算這兩個冪

　
: ((a×b)2×a)2×a×b

總共只需要 6 次乘法。注意 a×b 計算了兩次；結果可以在第一次計算後儲存，這將乘法計數減少到 5 次。

有數字的例子：

　
: 27×35 = ((2×3)2×2)2×2×3 = (62×2)2×6 = 722×6 = 31,104

如果至少有兩個指數大於 1 的話，同時計算冪就會比單獨計算減少乘法次數。

#### 使用變換

如果表達式在計算前進行變換，上面的例子 a7×b5 也可以只用 5 次乘法就計算出來：

a7×b5 = a2×(ab)5 其中 ab := a×b

　
: ab := a×b（一次乘法）

　
: a2×(ab)5 = ((ab)2×a)2×ab（四次乘法）

這個變換可以推廣成下面的方案： 
對於計算 aA×bB×...×mM×nN 
首先：定義 ab := a×b, abc = ab×c, ... 
然後：計算變換後的表達式 aA−B×abB−C×...×abc..mM−N×abc..mnN

在計算之前進行變換通常會減少乘法計數，但在某些情況下也會增加計數（請參見下面最後一個範例），因此在使用變換後的表達式進行計算之前，最好檢查一下乘法的次數。

#### 例子

對於下面的表達式，表中顯示了分開計算每個冪，在不進行變換的情況下同時進行計算，以及在變換後同時進行計算的乘法次數。

## 用有符號數字重新編碼

在某些計算中，如果允許負係數（也就會需要用基底的倒數）的話，只要在 ${\boldsymbol {G}}$ 中計算倒數很快或者已經預先計算，求冪會更加高效。例如，當計算 $x^{2^{k}-1}$ 時，二進位方法需要 $k-1$ 次乘法和 $k-1$ 次平方。不過可以用 $k$ 次平方得到 $x^{2^{k}}$ ，然後乘以 $x^{-1}$ 得到 $x^{2^{k}-1}$ 。

為此，定義以 $b$ 為基數的整數 $n$ 的有符號數字表示為

　
: $n=\sum _{i=0}^{l-1}n_{i}b^{i}{\text{  with  }}|n_{i}|<b$

_有符號二進位表示_也就是選取 $b=2$ ， $n_{i}\in \{-1,0,1\}$ 的表示法。記為 $(n_{l-1}\dots n_{0})_{s}$ 。有多種計算這種表示的方法。該表示不是唯一的，例如，取 $n=478$ 。 $(10{\bar {1}}1100{\bar {1}}10)_{s}$ 和 $(100{\bar {1}}1000{\bar {1}}0)_{s}$ 給出了兩個不同的有符號二進位表示，其中 ${\bar {1}}$ 表示 -1。由於在二進位方法中， $n$ 的基 2 表示的每個非零項都要計算乘法，因此感興趣的是找到非零項數量最少的有符號二進位表示，即具有_最小_漢明重量的表示。有一種方法是計算非相鄰形式（簡稱 NAF）的有符號二進位表示，它滿足對所有 $i\geqslant 0$ ， $n_{i}n_{i+1}=0$ ，記為 $(n_{l-1}\dots n_{0})_{\text{NAF}}$ 。例如，478 的 NAF 表示為 $(1000{\bar {1}}000{\bar {1}}0)_{\text{NAF}}$ 。這種表示總是有最小的漢明重量。下面的簡單演算法可以計算 $n_{l}=n_{l-1}=0$ 的整數 $n=(n_{l}n_{l-1}\dots n_{0})_{2}$ 的 NAF 表示：

```
$c_{0}=0$
for _i_ = 0 to _l_ - 1 do
  $c_{i+1}=\left\lfloor {\frac {1}{2}}(c_{i}+n_{i}+n_{i+1})\right\rfloor $
  $n_{i}'=c_{i}+n_{i}-2c_{i+1}$
return $(n_{l-1}'\dots n_{0}')_{\text{NAF}}$
```

Koyama 和 Tsuruoka 的另一種演算法並不要求 $n_{i}=n_{i+1}=0$ 這樣的條件；它仍然可以讓漢明重量最小化。

## 替代方法及推廣

平方求冪可以看作是一個次優的加法鏈求冪演算法：它通過由重複指數加倍（平方）和指數遞增（乘以 $x$ ）組成的加法鏈來計算指數。更一般地，如果允許_任何_先前計算的指數相加（通過乘以 _$x$_ 的冪），有時可以讓求冪運算的乘法次數更少（但通常使用更多的記憶體）。 $n=15$ 時的最少次數：

　
: $a^{15}=x\times (x\times [x\times x^{2}]^{2})^{2}\!$ （平方求冪，6 次乘法）

　
: $a^{15}=x^{3}\times ([x^{3}]^{2})^{2}\!$ （最佳加法鏈，在復用 $x^{3}$ 的情況下只需要 5 次乘法）

一般來說，求給定指數的最佳加法鏈是一個難題，因為沒有已知的高效演算法，所以最佳鏈通常只用於小指數（比如，在編譯器中已經預先儲存了小指數冪的最佳鏈）。不過，有一些啟發式演算法，雖然不是最佳的，但是由於額外的簿記工作和記憶體使用量的增加而導致的乘法次數少於平方求冪。無論如何，乘法的次數永遠不會比 Θ(log _n_) 增長得更慢，所以這些演算法只能減小平方求冪的漸進複雜度的常數因子。

## 參見

* 模冪運算
* 向量加法鏈
* 蒙哥馬利演算法
* 非相鄰形式
* 加法鏈

## 注釋

## 參考文獻