平方求冪

在數學和程式設計中，平方求冪（英語：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})$ 。

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

  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);

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

  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).

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

  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$ 求模

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 _{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(n_i)
    for j:=1 to k-s do
        y := y² 
    y := y*x^u
    for j:=1 to s do
        y := y²
    i := i-1
return y

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

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

的最小整數。^[1]

滑動窗口法

此方法是 $2^{k}$ 法的更高效的變體。例如，要計算398次冪，二進制展開為 (110 001 110)₂，採用長度為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)n₂ 的值

以下是一般演算法：

演算法：

輸入: 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_i=0 then
        y:=y²' i:=i-1
    else
        s:=max{i-k+1,0}
        while n_s=0 do
            s:=s+1 ^{[notes 1]}
        for h:=1 to i-s+1 do
            y:=y²
        u:=(n_i,n_i-1,....,n_s)₂
        y:=y*x^u
        i:=s-1
return y

蒙哥馬利階梯法

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

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

x₁=x; x₂=x²
for i=k-2 to 0 do
  If n_i=0 then
    x₂=x₁*x₂; x₁=x₁²
  else
    x₁=x₁*x₂; x₂=x₂²
return x₁

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

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

固定基底的冪

當基底固定而指數變化時，可以使用幾種方法來計算 $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 n_i=j then u=u*x^b_i
   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$ 個元素。^[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}$ 。

範例實現

通過2的冪進行計算

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

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

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

執行實例：計算 3¹⁰

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

Iteration 1
  n = 10 -> n is even
  x := x² = 3² = 9
  n := n / 2 = 5

Iteration 2
  n = 5 -> n is odd
      -> result := result * x = 1 * x = 1 * 3² = 9
         n := n - 1 = 4
  x := x² = 9² = 3⁴ = 81
  n := n / 2 = 2

Iteration 3
  n = 2 -> n is even
  x := x² = 81² = 3⁸ = 6561
  n := n / 2 = 1

Iteration 4
  n = 1 -> n is odd
      -> result := result * x = 3² * 3⁸ = 3¹⁰ = 9 * 6561 = 59049
         n := n - 1 = 0

return result

執行實例：計算 3¹⁰

result := 3
bin := "1010"

Iteration for digit 2:
  result := result² = 3² = 9
  1010_bin - Digit equals "0"

Iteration for digit 3:
  result := result² = (3²)² = 3⁴  = 81
  1010_bin - Digit equals "1" --> result := result*3 = (3²)²*3 = 3⁵  = 243

Iteration for digit 4:
  result := result² = ((3²)²*3)² = 3¹⁰  = 59049
  1010_bin - Digit equals "0"

return result

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

冪的乘積的計算

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

例子

式子 a⁷×b⁵ 可以分三步計算：

((a)²×a)²×a （計算 a⁷ 需要四次乘法）

((b)²)²×b （計算 b⁵ 需要三次乘法）

(a⁷)×(b⁵) （計算二者乘積需要一次乘法）

所以總共需要八次乘法。

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

((a×b)²×a)²×a×b

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

有數字的例子：

2⁷×3⁵ = ((2×3)²×2)²×2×3 = (6²×2)²×6 = 72²×6 = 31,104

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

使用變換

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

a⁷×b⁵ = a²×(ab)⁵ 其中 ab := a×b

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

a²×(ab)⁵ = ((ab)²×a)²×ab（四次乘法）

這個變換可以推廣成下面的方案：
對於計算 a^A×b^B×...×m^M×n^N
首先：定義 ab := a×b, abc = ab×c, ...
然後：計算變換後的表達式 a^A−B×ab^B−C×...×abc..m^M−N×abc..mn^N

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

例子

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

例子	a⁷×b⁵×c³	a⁵×b⁵×c³	a⁷×b⁴×c¹
分開計算	[((a)²×a)²×a] × [((b)²)²×b] × [(c)²×c] （11次乘法）	[((a)²)²×a] × [((b)²)²×b] × [(c)²×c] （10次乘法）	[((a)²×a)²×a] × [((b)²)²] × [c] （8次乘法）
同時計算	((a×b)²×a×c)²×a×b×c （8次乘法）	((a×b)²×c)²×a×b×c （7次乘法）	((a×b)²×a)²×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)²×abc （4次乘法 ⇒ 總共6次）	(ab×abc)²×abc （3次乘法 ⇒ 總共5次）	(a×ab)²×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$ 的有符號數字表示（英語：signed-digit representation）為

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表示的每個非零項都要計算乘法，因此感興趣的是找到非零項數量最少的有符號二進制表示，即具有最小漢明重量的表示。有一種方法是計算非相鄰形式（英語：non-adjacent form）（簡稱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$ 這樣的條件；它仍然可以讓漢明重量最小化。

替代方法及推廣

平方求冪可以看作是一個次優的加法鏈求冪（英語：addition-chain exponentiation）演算法：它通過由重複指數加倍（平方）和指數遞增（乘以 $x$ ）組成的加法鏈（英語：addition chain）來計算指數。更一般地，如果允許任何先前計算的指數相加（通過乘以 $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) 增長得更慢，所以這些演算法只能減小平方求冪的漸進複雜度的常數因子。

參見

註釋

^ 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 $x^{2^{k}-1}$ need be computed and only those specifically involved in the computation need be considered.

參考文獻

^ ^1.0 ^1.1 Cohen, H.; Frey, G. (編). Handbook of Elliptic and Hyperelliptic Curve Cryptography. Discrete Mathematics and Its Applications. Chapman & Hall/CRC. 2006. ISBN 9781584885184.
^ Montgomery, Peter L. Speeding the Pollard and Elliptic Curve Methods of Factorization (PDF). Math. Comput. 1987, 48 (177): 243–264 [2018-02-17]. （原始內容存檔 (PDF)於2018-01-27）.

[2] 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 $x^{2^{k}-1}$ need be computed and only those specifically involved in the computation need be considered.

[frey-1] 1.0 ^1.1 Cohen, H.; Frey, G. (編). Handbook of Elliptic and Hyperelliptic Curve Cryptography. Discrete Mathematics and Its Applications. Chapman & Hall/CRC. 2006. ISBN 9781584885184.

[ladder-3] Montgomery, Peter L. Speeding the Pollard and Elliptic Curve Methods of Factorization (PDF). Math. Comput. 1987, 48 (177): 243–264 [2018-02-17]. （原始內容存檔 (PDF)於2018-01-27）.

[1]

[notes 1]

[2]

平方求冪

目次

基本方法

計算複雜度

$2^{k}$ 法

滑動窗口法

蒙哥馬利階梯法

固定基底的冪

姚期智的方法

歐幾里得法

更多應用

範例實現

通過2的冪進行計算

執行實例：計算 3¹⁰

執行實例：計算 3¹⁰

冪的乘積的計算

例子

使用變換

例子

用有符號數字重新編碼

替代方法及推廣

參見

註釋

參考文獻

平方求冪

基本方法

計算複雜度

2 k {\displaystyle 2^{k}} 法

滑動窗口法

蒙哥馬利階梯法

固定基底的冪

姚期智的方法

歐幾里得法

更多應用

範例實現

通過2的冪進行計算

執行實例：計算 310

執行實例：計算 310

冪的乘積的計算

例子

使用變換

例子

用有符號數字重新編碼

替代方法及推廣

參見

註釋

參考文獻

$2^{k}$ 法

執行實例：計算 3¹⁰

執行實例：計算 3¹⁰