平方求冪

在數學和程序設計中，平方求冪（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)$

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

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

```pascal
  Function exp_by_squaring(x, n)
    if n < 0  then return exp_by_squaring(1 / x, -n);
    else if n = 0  then return  1;
    else if n = 1  then return  x ;
    else if n is even  then return exp_by_squaring(x * x,  n / 2);
    else if n is odd  then return x * exp_by_squaring(x * x, (n - 1) / 2);
```

儘管不是尾調用，但是通過引入輔助函數，該演算法可以被重寫成尾遞歸演算法：

```pascal
  Function exp_by_squaring(x, n)
    exp_by_squaring2(1, x, n)
  Function exp_by_squaring2(y, x, n)
    if n < 0  then return exp_by_squaring2(y, 1 / x, - n);
    else if n = 0  then return  y;
    else if n = 1  then return  x * y;
    else if n is even  then return exp_by_squaring2(y, x * x,  n / 2);
    else if n is odd  then return exp_by_squaring2(x * y, x * x, (n - 1) / 2).
```

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

```pascal
  Function exp_by_squaring_iterative(x, n)
    if n < 0 then
      x := 1 / x;
      n := -n;
    if n = 0 then return 1
    y := 1;
    while n > 1 do
      if n is even then 
        x := x * x;
        n := n / 2;
      else
        y := x * y;
        x := x * x;
        n := (n – 1) / 2;
    return x * y
```

c++實現（非遞歸）

返回$p^n$ 對$Mod$求模

```c++
long long power(long long p, long long n)
{
	long long ans = 1;
	while (n)
	{
		if (n & 1) ans = (ans * p) % Mod;
		p = p * p % Mod;
		n >>= 1;
	}
	return ans;
}
```

## 計算複雜度

簡要分析表明此演算法用了 $\lfloor \log_2n\rfloor$ 次平方，以及至多 $\lfloor \log_2n\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(n<sub>i</sub>)
    '''for''' j:=1 '''to''' k-s '''do'''
        y := y<sup>2</sup> 
    y := y*x<sup>u</sup>
    '''for''' j:=1 '''to''' s '''do'''
        y := y<sup>2</sup>
    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''' n<sub>i</sub>=0 '''then'''
        y:=y<sup>2</sup>' i:=i-1
    '''else'''
        s:=max{i-k+1,0}
        '''while''' n<sub>s</sub>=0 '''do'''
            s:=s+1 <ref group=notes>In this line, the loop finds the longest string of length less than or equal to 'k' which ends in a non zero value. And not all odd powers of 2 up to <math style="vertical-align:baseline;">x^{2^k-1}</math> need be computed and only those specifically involved in the computation need be considered.</ref>
        '''for''' h:=1 '''to''' i-s+1 '''do'''
            y:=y<sup>2</sup>
        u:=(n<sub>i</sub>,n<sub>i-1</sub>,....,n<sub>s</sub>)<sub>2</sub>
        y:=y*x<sup>u</sup>
        i:=s-1
'''return''' y
```

## 蒙哥馬利階梯法

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

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

```
x<sub>1</sub>=x; x<sub>2</sub>=x<sup>2</sup>
for i=k-2 to 0 do
  If n<sub>i</sub>=0 then
    x<sub>2</sub>=x<sub>1</sub>*x<sub>2</sub>; x<sub>1</sub>=x<sub>1</sub><sup>2</sup>
  else
    x<sub>1</sub>=x<sub>1</sub>*x<sub>2</sub>; x<sub>2</sub>=x<sub>2</sub><sup>2</sup>
return x<sub>1</sub>
```

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

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

## 固定基底的冪

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

### 姚期智的方法

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

把指數 $n$ 寫成

$n = \sum_{i=0}^{w-1} n_ib_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 n<sub>i</sub>=j then u=u*x<sup>b<sub>i</sub></sup>
   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 <math>M \in \left[ 0, l - 1 \right]</math>, such that <math>\forall i \in \left[ 0, l - 1 \right], {n_M} \ge {n_i}</math>;
        Find <math>N \in \left( \left[ 0, l - 1 \right] - M \right)</math>, such that <math>\forall i \in \left( \left[ 0, l - 1 \right] - M \right), {n_N} \ge {n_i}</math>;
        '''Break loop''' if <math>{n_N} = 0</math>;
        '''Let''' <math>q = \left\lfloor {n_M} / {n_N} \right\rfloor</math>, and then '''let''' <math>{n_N} = \left( {n_M} \bmod {n_N} \right)</math>;
        Compute recursively <math>{x_M}^q</math>, and then '''let''' <math>{x_N} = {x_N} \cdot {x_M}^q</math>;
    '''End loop''';
    '''Return''' <math>x^n = {x_M}^{n_M}</math>.
```

該演算法首先在 $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 對那些語言來說，就是冗餘的步驟了。

n0是n的二進位表示的最右邊的位，所以如果它是1，則該數是奇數，如果它是零，則該數是偶數。它也是以2為模n的餘數。

```ruby
def power(x,n)
  result = 1
  while n.nonzero?
    if n[0].nonzero?
      result *= x
      n -= 1
    end
    x *= x
    n /= 2
  end
  return result
end
```

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

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

'''Iteration 1'''
  n = 10 -&gt; n is even
  x := x<sup>2</sup> = 3<sup>2</sup> = 9
  n := n / 2 = 5

'''Iteration 2'''
  n = 5 -&gt; n is odd
      -&gt; result := result * x = 1 * x = 1 * 3<sup>2</sup> = 9
         n := n - 1 = 4
  x := x<sup>2</sup> = 9<sup>2</sup> = 3<sup>4</sup> = 81
  n := n / 2 = 2

'''Iteration 3'''
  n = 2 -&gt; n is even
  x := x<sup>2</sup> = 81<sup>2</sup> = 3<sup>8</sup> = 6561
  n := n / 2 = 1

'''Iteration 4'''
  n = 1 -&gt; n is odd
      -&gt; result := result * x = 3<sup>2</sup> * 3<sup>8</sup> = 3<sup>10</sup> = 9 * 6561 = 59049
         n := n - 1 = 0

return result
```

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

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

'''Iteration for digit 2:'''
  result := result<sup>2</sup> = 3<sup>2</sup> = 9
  1'''0'''10<sub>bin</sub> - Digit equals "0"

'''Iteration for digit 3:'''
  result := result<sup>2</sup> = (3<sup>2</sup>)<sup>2</sup> = 3<sup>4</sup>  = 81
  10'''1'''0<sub>bin</sub> - Digit equals "1" --&gt;&nbsp;result := result*3 = (3<sup>2</sup>)<sup>2</sup>*3 = 3<sup>5</sup>  = 243

'''Iteration for digit 4:'''
  result := result<sup>2</sup> = ((3<sup>2</sup>)<sup>2</sup>*3)<sup>2</sup> = 3<sup>10</sup>  = 59049
  101'''0'''<sub>bin</sub> - 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

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

#### 例子

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

例子｜a7×b5×c3｜a5×b5×c3｜a7×b4×c1

分開計算｜((a)2×a)2×a × ((b)2)2×b × (c)2×c（11次乘法）｜((a)2)2×a × ((b)2)2×b × (c)2×c（10次乘法）｜((a)2×a)2×a × ((b)2)2 × c（8次乘法）

同時計算｜((a×b)2×a×c)2×a×b×c（8次乘法）｜((a×b)2×c)2×a×b×c（7次乘法）｜((a×b)2×a)2×a×c（6次乘法）

變換｜a := 2 ab := a×b abc := ab×c（2次乘法）｜a := 2 ab := a×b abc := ab×c（2次乘法）｜a := 2 ab := a×b abc := ab×c（2次乘法）

之後的計算｜(a×ab×abc)2×abc（4次乘法 ⇒ 總共6次）｜(ab×abc)2×abc（3次乘法 ⇒ 總共5次）｜(a×ab)2×a×ab×abc（5次乘法 ⇒ 總共7次）

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

在某些計算中，如果允許負係數（也就會需要用基底的倒數）的話，只要在 $\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_ib^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_in_{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_ln_{l-1}\dots n_0)_2$ 的NAF表示：

```
<math>c_0=0</math>
for 1=''i'' = 0 to ''l'' - 1 do
  <math>c_{i+1}=\left\lfloor\frac{1}{2}(c_i+n_i+n_{i+1})\right\rfloor</math>
  <math>n_i'=c_i+n_i-2c_{i+1}</math>
return <math>(n_{l-1}'\dots n_0')_{\text{NAF}}</math>
```

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) 增長得更慢，所以這些演算法只能減小平方求冪的漸進複雜度的常數因子。