平方求幂

在数学和程序设计中，平方求幂（英语：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) 增长得更慢，所以这些算法只能减小平方求幂的渐进复杂度的常数因子。

参见

模幂运算
向量加法链（英语：Vectorial addition chain）
蒙哥马利算法
非相邻形式（英语：Non-adjacent form）
加法链（英语：Addition chain）

注释

^ 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¹⁰