Exponentiation by squaring

In mathematics and computer programming, exponentiating by squaring is a general method for fast computation of large positive integer powers of a number, or more generally of an element of a semigroup, like a polynomial or a square matrix. Some variants are commonly referred to as square-and-multiply algorithms or binary exponentiation. These can be of quite general use, for example in modular arithmetic or powering of matrices. For semigroups for which additive notation is commonly used, like elliptic curves used in cryptography, this method is also referred to as double-and-add.

Basic method

Recursive version

The method is based on the observation that, for any integer <math>n > 0</math>, one has

<math display="block"> x^n=

\begin{cases}

x \, (x^2)^{(n-1)/2} & \text{if } n \text{ is odd}, \\

(x^2)^{n/2} & \text{if } n \text{ is even}.

\end{cases}

</math>

If the exponent is zero, then the answer is 1. If the exponent is negative then we can reuse the previous formula by rewriting the value using a positive exponent. That is,

<math display="block">

x^n = \left(\frac{1}{x}\right)^{-n}.

</math>

Together, these may be implemented directly as the following recursive algorithm:

Inputs: a real number x; an integer n

Output: xⁿ

function exp_by_squaring(x, n) is

if n < 0 then

return exp_by_squaring(1 / x, −n)

else if n = 0 then

return 1

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)

end function

In each recursive call, the least-significant digit of the binary representation of is removed. It follows that the number of recursive calls is <math>\lceil \log_2 n\rceil,</math> the number of bits of the binary representation of . So this algorithm computes this number of squares and a lower number of multiplication, which is equal to the number of 1s in the binary representation of . This logarithmic number of operations is to be compared with the trivial algorithm which requires multiplications.

This algorithm is not tail-recursive. This implies that it requires an amount of auxiliary memory that is roughly proportional to the number of recursive calls, or perhaps higher if the amount of data per iteration is increasing.

The algorithms of the next section use a different approach, and the resulting algorithms needs the same number of operations, but use an auxiliary memory that is roughly the same as the memory required to store the result.

With constant auxiliary memory

The variants described in this section are based on the formula

<math display="block">

y x^n = \begin{cases}

yx\,(x^2)^{(n-1)/2} & \text{if } n \text{ is odd}, \\

y\,(x^2)^{n/2} & \text{if } n \text{ is even}.

\end{cases}

</math>

If one applies recursively this formula, by starting with , one gets eventually an exponent equal to , and the desired result is then the left factor.

This may be implemented as a tail-recursive function:

<syntaxhighlight lang="pascal">

Function exp_by_squaring(x, n)

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

</syntaxhighlight>

The iterative version of the algorithm also uses a bounded auxiliary space, and is given by

<syntaxhighlight lang="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 odd then

y := x * y;

n := n - 1;

x := x * x;

n := n / 2;

return x * y

</syntaxhighlight>

The correctness of the algorithm results from the fact that <math>yx^n</math> is invariant during the computation; it is <math>1 \cdot x^n = x^n</math> at the beginning; and it is <math>yx^1 = xy</math> at the end.

These algorithms use exactly the same number of operations as the algorithm of the preceding section, but the multiplications are done in a different order.

Computational complexity

A brief analysis shows that such an algorithm uses <math>\lfloor \log_2n\rfloor</math> squarings and at most <math>\lfloor \log_2n\rfloor</math> multiplications, where <math>\lfloor\;\rfloor</math> denotes the floor function. More precisely, the number of multiplications is one less than the number of ones present in the binary expansion of n. For n greater than about 4 this is computationally more efficient than naively multiplying the base with itself repeatedly.

Each squaring results in approximately double the number of digits of the previous, and so, if multiplication of two d-digit numbers is implemented in O(d^k) operations for some fixed k, then the complexity of computing xⁿ is given by

: <math>

\sum\limits_{i=0}^{O(\log n)} \big(2^i O(\log x)\big)^k = O\big((n \log x)^k\big).

</math>

2^k-ary method

This algorithm calculates the value of xⁿ after expanding the exponent in base 2^k. It was first proposed by Brauer in 1939. In the algorithm below we make use of the following function f(0) = (k, 0) and f(m) = (s, u), where m = u·2^s with u odd.

Algorithm:

;Input: An element x of G, a parameter k > 0, a non-negative integer and the precomputed values <math>x^3, x^5, ... , x^{2^k-1}</math>.

;Output: The element xⁿ in G

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

For optimal efficiency, k should be the smallest integer satisfying

: <math>\lg n < \frac{k(k + 1) \cdot 2^{2k{2^{k+1} - k - 2} + 1.</math>

Sliding-window method

This method is an efficient variant of the 2^k-ary method. For example, to calculate the exponent 398, which has binary expansion (110 001 110)₂, we take a window of length 3 using the 2^k-ary method algorithm and calculate 1, x³, x⁶, x¹², x²⁴, x⁴⁸, x⁴⁹, x⁹⁸, x⁹⁹, x¹⁹⁸, x¹⁹⁹, x³⁹⁸.

But, we can also compute 1, x³, x⁶, x¹², x²⁴, x⁴⁸, x⁹⁶, x¹⁹², x¹⁹⁹, x³⁹⁸, which saves one multiplication and amounts to evaluating (110 001 110)₂

Here is the general algorithm:

Algorithm:

;Input: An element x of G, a non negative integer , a parameter k > 0 and the pre-computed values <math>x^3, x^5, ... ,x^{2^k-1}</math>.

;Output: The element xⁿ &isin; G.

Algorithm:

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

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

Montgomery's ladder technique

Many algorithms for exponentiation do not provide defence against side-channel attacks. Namely, an attacker observing the sequence of squarings and multiplications can (partially) recover the exponent involved in the computation. This is a problem if the exponent should remain secret, as with many public-key cryptosystems. A technique called "Montgomery's ladder" addresses this concern.

Given the binary expansion of a positive, non-zero integer n = (n_k−1...n₀)₂ with n_k−1 = 1, we can compute xⁿ as follows:

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₁

The algorithm performs a fixed sequence of operations (up to log n): a multiplication and squaring takes place for each bit in the exponent, regardless of the bit's specific value. A similar algorithm for multiplication by doubling exists.

This specific implementation of Montgomery's ladder is not yet protected against cache timing attacks: memory access latencies might still be observable to an attacker, as different variables are accessed depending on the value of bits of the secret exponent. Modern cryptographic implementations use a "scatter" technique to make sure the processor always misses the faster cache.

Fixed-base exponent

There are several methods which can be employed to calculate xⁿ when the base is fixed and the exponent varies. As one can see, precomputations play a key role in these algorithms.

Yao's method

Yao's method is orthogonal to the -ary method where the exponent is expanded in radix and the computation is as performed in the algorithm above. Let , , , and be integers.

Let the exponent be written as

: <math> n = \sum_{i=0}^{w-1} n_i b_i,</math>

where <math>0 \leqslant n_i < h</math> for all <math>i \in [0, w-1]</math>.

Let .

Then the algorithm uses the equality

: <math>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.</math>

Given the element of , and the exponent written in the above form, along with the precomputed values , the element is calculated using the algorithm below:

y = 1, u = 1, 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

If we set and , then the values are simply the digits of in base . Yao's method collects in u first those that appear to the highest power ; in the next round those with power are collected in as well etc. The variable y is multiplied times with the initial , times with the next highest powers, and so on.

The algorithm uses multiplications, and elements must be stored to compute .