# 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

The method is based on the observation that, for a positive integer *n*, we have

This may be easily implemented as the following recursive algorithm:

Functionexp-by-squaring(x,n)ifn<0then returnexp-by-squaring(1/x, -n);else ifn=0then return1;else ifn=1then returnx;else ifnis eventhen returnexp-by-squaring(x^{2},n/2);else ifnis oddthen returnx*exp-by-squaring(x^{2}, (n-1)/2).

Moreover, the above algorithm can be rewritten to be tail-recursive and minimizes the stack-depth, making it efficient.

## Computational complexity

A brief analysis shows that such an algorithm uses O(log_{2}*n*) squarings and O(log_{2}*n*) multiplications. For n > 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*^{n} is given by:

## 2^{k}-ary method

This algorithm calculates the value of x^{n} 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
*n*= (*n*_{l−1},*n*_{l−2}, ...,*n*_{0})_{2k}and the precomputed values .

- Output
- The element x
^{n}in*G*

1. y := 1; i := l-1 2.whilei>=0 do 3. (s,u) := f(n_{i}) 4.forj:=1tok-sdo5. y := y^{2}6. y := y*x^{u}7.forj:=1tosdo8. y := y^{2}9. i := i-1 10.returny

For optimal efficiency, *k* should be the smallest integer satisfying ^{[1]}

## 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)_{2}, we take a window of length 3 using the 2^{k}-ary method algorithm we calculate 1,x^{3},x^{6},x^{12},x^{24},x^{48},x^{49},x^{98},x^{99},x^{198},x^{199},x^{398}.
But, we can also compute 1,x^{3},x^{6},x^{12},x^{24},x^{48},x^{96},x^{192},x^{199},
x^{398} which saves one multiplication and amounts to evaluating (110 001 110)n_{2}

Here is the general algorithm:

Algorithm:

- Input
- An element
*x*of*G*,a non negative integer*n*=(*n*_{l},*n*_{l-1},...,*n*_{0})_{2}, a parameter k>0 and the pre-computed values .

- Output
- The element
*x*∈^{n}*G*

Algorithm:

1. y := 1; i := l-1 2.whilei > -1do3.ifn_{i}=0theny:=y^{2}' i:=i-1 4.else5. s:=max{i-k+1,0} 6.whilen_{s}=0dos:=s+1^{[2]}7.forh:=1toi-s+1doy:=y^{2}8. u:=(n_{i},n_{i-1},....,n_{s})_{2}9. y:=y*x^{u}10. i:=s-1 11.returny

## 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^{[3]} addresses this concern.

Given the binary expansion of a positive, non-zero integer n=(n_{k-1}...n_{0})_{2} with n_{k-1}=1 we can compute x^{n} as follows:

x_{1}=x; x_{2}=x^{2}for i=k-2 to 0 do If n_{i}=0 then x_{2}=x_{1}*x_{2}; x_{1}=x_{1}^{2}else x_{1}=x_{1}*x_{2}; x_{2}=x_{2}^{2}return x_{1}

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.

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 you access different variables depending on the value of bits of the secret exponent.

## Fixed base exponent

There are several methods which can be employed to calculate x^{n} 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 2^{k}-ary method where the exponent is expanded in radix b=2^{k} and the computation is as performed in the algorithm above. Let "n", "n_{i}", "b", and "b_{i}" be integers.

Let the exponent "n" be written as

Let x_{i} = x^{bi}. Then the algorithm uses the equality

Given the element 'x' of G, and the exponent 'n' written in the above form, along with the precomputed values x^{b0}....x^{bw-1} the element x^{n} is calculated using the algorithm below.

`
`

- 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^{bi}

- if n
- y=y*u
- j=j-1

- for i=0 to w-1 do
- return y

`
If we set h=2`^{k} and b_{i} = h^{i} then the n_{i} 's are simply the digits of n in base h. Yao's method collects in u first those x_{i} which appear to the highest power h-1; in the next round those with power h-2 are collected in u as well etc. The variable y is multiplied h-1 times with the initial u, h-2 times with the next highest powers etc.
The algorithm uses w+h-2 multiplications and w+1 elements must be stored to compute x^{n} (see ^{[1]}).

### Euclidean method

The Euclidean method was first introduced in *Efficient exponentiation using precomputation and vector addition chains* by P.D Rooij.

This method for computing in group **G**, where is a natural integer, whose algorithm is given below, is using the following equality recursively:

- , where
- (in other words an Euclidian division of the exponent
*n*_{1}by*n*_{0}is used to return a quotient Template:Mvar and a rest ).

Given the base element in group **G**, and the exponent written as in Yao's method, the element is calculated using precomputed values and then the algorithm below.

Begin loopFind , such that ; Find , such that ;Break loopif ;Let, and thenlet; Compute recursively , and thenlet;End loop;Return.

The algorithm first finds the largest value amongst the *n*_{i} and then the supremum within the set of Template:( *n*_{i} \ *i* ≠ *M* Template:).
Then it raises *x*_{M} to the power Template:Mvar, multiplies this value with *x*_{N}, and then assigns *x*_{N} the result of this computation and *n*_{M} the value *n*_{M} modulo *n*_{N}.

## Further applications

The same idea allows fast computation of large exponents modulo a number. Especially in cryptography, it is useful to compute powers in a ring of integers modulo *q*. It can also be used to compute integer powers in a group, using the rule

- Power(
*x*, −*n*) = (Power(*x*,*n*))^{−1}.

The method works in every semigroup and is often used to compute powers of matrices,

For example, the evaluation of

- 13789
^{722341}(mod 2345)

would take a very long time and lots of storage space if the naïve method were used: compute 13789^{722341} then take the remainder when divided by 2345. Even using a more effective method will take a long time: square 13789, take the remainder when divided by 2345, multiply the result by 13789, and so on. This will take less than modular multiplications.

Applying above *exp-by-squaring* algorithm, with "*" interpreted as *x***y* = *xy* mod 2345 (that is a multiplication followed by a division with remainder) leads to only 27 multiplications and divisions of integers which may all be stored in a single machine word.

## Example implementations

### Computation by powers of 2

This is a non-recursive implementation of the above algorithm in Ruby.

In most statically typed languages, `result=1` must be replaced with code assigning an identity matrix of the same size as `x` to `result` to get a matrix exponentiating algorithm. In Ruby, thanks to coercion, `result` is automatically upgraded to the appropriate type, so this function works with matrices as well as with integers and floats. Note that n=n-1 is redundant when n=n/2 implicitly rounds towards zero, as lower level languages would do.
n[0] is the rightmost bit of the binary representation of n, so if it is 1, the number is odd, if it is zero, the number is even.

```
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
```

#### Runtime example: compute 3^{10}

parameter x = 3 parameter n = 10 result := 1Iteration 1n = 10 -> n is even x := x^{2}= 3^{2}= 9 n := n / 2 = 5Iteration 2n = 5 -> n is odd -> result := result * x = 1 * x = 1 * 3^{2}= 9 n := n - 1 = 4 x := x^{2}= 9^{2}= 3^{4}= 81 n := n / 2 = 2Iteration 3n = 2 -> n is even x := x^{2}= 81^{2}= 3^{8}= 6561 n := n / 2 = 1Iteration 4n = 1 -> n is odd -> result := result * x = 3^{2}* 3^{8}= 3^{10}= 9 * 6561 = 59049 n := n - 1 = 0 return result

#### Runtime example: compute 3^{10}

result := 3 bin := "1010"Iteration for digit 2:result := result^{2}= 3^{2}= 9 1010_{bin}- Digit equals "0"Iteration for digit 3:result := result^{2}= (3^{2})^{2}= 3^{4}= 81 1010_{bin}- Digit equals "1" --> result := result*3 = (3^{2})^{2}*3 = 3^{5}= 243Iteration for digit 4:result := result^{2}= ((3^{2})^{2}*3)^{2}= 3^{10}= 59049 1010_{bin}- Digit equals "0" return result

JavaScript-Demonstration: http://home.mnet-online.de/wzwz.de/temp/ebs/en.htm

### Calculation of products of powers

Exponentiation by squaring may also be used to calculate the product of 2 or more powers. If the underlying group or semigroup is commutative then it is often possible to reduce the number of multiplication by computing the product simultaneously.

#### Example

The formula a^{7}×b^{5} may be calculated within 3 steps:

- ((a)
^{2}×a)^{2}×a (four multiplications for calculating a^{7}) - ((b)
^{2})^{2}×b (three multiplications for calculating b^{5}) - (a
^{7})×(b^{5}) (one multiplication to calculate the product of the two)

so one gets eight multiplications in total.

A faster solution is to calculate both powers simultaneously:

- ((a×b)
^{2}×a)^{2}×a×b

which needs only 6 multiplications in total. Note that a×b is calculated twice, the result could be stored after the first calculation which reduces the count of multiplication to 5.

Example with numbers:

- 2
^{7}×3^{5}= ((2×3)^{2}×2)^{2}×2×3 = (6^{2}×2)^{2}×6 = 72^{2}×6 = 31,104

Calculating the powers simultaneously instead of calculating them separately always reduces the count of multiplications if at least two of the exponents are greater than 1.

#### Using transformation

The example above a^{7}×b^{5} may also be calculated with only 5
multiplications if the expression is transformed before calculation:

a^{7}×b^{5} = a^{2}×(ab)^{5} with ab := a×b

- ab := a×b (one multiplication)
- a
^{2}×(ab)^{5}= ((ab)^{2}×a)^{2}×ab (four multiplications)

Generalization of transformation shows the following scheme:

For calculating a^{A}×b^{B}×...×m^{M}×n^{N}

1st: define ab := a×b, abc = ab×c, ...

2nd: calculate the transformed expression a^{A−B}×ab^{B−C}×...×abc..m^{M−N}×abc..mn^{N}

Transformation before calculation often reduces the count of multiplications but in some cases it also increases the count (see the last one of the examples below), so it may be a good idea to check the count of multiplications before using the transformed expression for calculation.

#### Examples

For the following expressions the count of multiplications is shown for calculating each power separately, calculating them simultaneously without transformation and calculating them simultaneously after transformation.

Example: a^{7}×b^{5}×c^{3}

separate: [((a)^{2}×a)^{2}×a] × [((b)^{2})^{2}×b] × [(c)^{2}×c] ( **11** multiplications )

simultaneous: ((a×b)^{2}×a×c)^{2}×a×b×c ( **8** multiplications )

transformation: a := 2 ab := a×b abc := ab×c ( 2 multiplications )

calculation after that: (a×ab×abc)^{2}×abc ( 4 multiplications ⇒ **6** in total )

Example: a^{5}×b^{5}×c^{3}

separate: [((a)^{2})^{2}×a] × [((b)^{2})^{2}×b] × [(c)^{2}×c] ( **10** multiplications )

simultaneous: ((a×b)^{2}×c)^{2}×a×b×c ( **7** multiplications )

transformation: a := 2 ab := a×b abc := ab×c ( 2 multiplications )

calculation after that: (ab×abc)^{2}×abc ( 3 multiplications ⇒ **5** in total )

Example: a^{7}×b^{4}×c^{1}

separate: [((a)^{2}×a)^{2}×a] × [((b)^{2})^{2}] × [c] ( **8** multiplications )

simultaneous: ((a×b)^{2}×a)^{2}×a×c ( **6** multiplications )

transformation: a := 2 ab := a×b abc := ab×c ( 2 multiplications )

calculation after that: (a×ab)^{2}×a×ab×abc ( 5 multiplications ⇒ **7** in total )

## Signed-digit recoding

In certain computations it may be more efficient to allow negative coefficients and hence use the inverse of the base, provided inversion in G is ' fast' or has been precomputed. For example, when computing x^{2k−1} the binary method requires k−1 multiplications and k−1 squarings . However one could perform k squarings to get x^{2k} and then multiply by x^{−1} to obtain x^{2k−1}.

To this end we define the signed-digit representation of an integer in radix as

'Signed Binary Representation' corresponds to the particular choice and . It is denoted by . There are several methods for computing this representation. The representation is not unique, for example take . Two distinct signed-binary representations are given by and , where is used to denote . Since the binary method computes a multiplication for every non-zero entry in the base 2 representation of , we are interested in finding the signed-binary representation with the smallest number of non-zero entries, that is, the one with *minimal* Hamming weight. One method of doing this is to compute the representation in non-adjacent form, or NAF for short, which is one that satisfies and denoted by . For example the NAF representation of 478 is equal to . This representation always has minimal Hamming weight. A simple algorithm to compute the NAF representation of a given integer with is the following:

Another algorithm by Koyama and Tsuruoka does not require the condition that ; it still minimizes the Hamming weight.

## Alternatives and generalizations

{{#invoke:main|main}}
Exponentiation by squaring can be viewed as a suboptimal addition-chain exponentiation algorithm: it computes the exponent via an addition chain consisting of repeated exponent doublings (squarings) and/or incrementing exponents by *one* (multiplying by *x*) only. More generally, if one allows *any* previously computed exponents to be summed (by multiplying those powers of *x*), one can sometimes perform the exponentiation using fewer multiplications (but typically using more memory). The smallest power where this occurs is for *n*=15:

In general, finding the *optimal* addition chain for a given exponent is a hard problem, for which no efficient algorithms are known, so optimal chains are typically only used for small exponents (e.g. in compilers where the chains for small powers have been pre-tabulated). However, there are a number of heuristic algorithms that, while not being optimal, have fewer multiplications than exponentiation by squaring at the cost of additional bookkeeping work and memory usage. Regardless, the number of multiplications never grows more slowly than Θ(log *n*), so these algorithms only improve asymptotically upon exponentiation by squaring by a constant factor at best.

## See also

- Modular exponentiation
- Vectorial addition chain
- Montgomery reduction
- Non-adjacent form
- Addition chain

## Notes

- ↑
^{1.0}^{1.1}Cohen, H., Frey, G. (editors): Handbook of elliptic and hyperelliptic curve cryptography. Discrete Math.Appl., Chapman & Hall/CRC (2006) - ↑ 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 need be computed and only those specifically involved in the computation need be considered.
- ↑ Montgomery, P. L. "Speeding the Pollard and Elliptic Curve Methods of Factorization." Math. Comput. 48, 243-264, 1987.