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 method uses the bits of the exponent to determine which powers are computed.
This example shows how to compute formula_2 using this method.
The exponent, 13, is 1101 in binary. The bits are used in left to right order.
The exponent has 4 bits, so there are 4 iterations.
First, initialize the result to 1: formula_3.
If we write formula_12 in binary as formula_13, then this is equivalent to defining a sequence formula_14 by letting formula_15 and then defining formula_16 for formula_17, where formula_18 will equal formula_19.
This may be implemented as the following recursive algorithm: 
 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);
Although not tail-recursive, this algorithm may be rewritten into a tail recursive algorithm by introducing an auxiliary function:
 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 = 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).
A tail-recursive variant may also be constructed using a pair of accumulators instead of an auxiliary function as seen in the F# example below. The accumulators a1 and a2 can be thought of as storing the values formula_20 and formula_21 where i and j are initialized to 1 and 0 respectively. In the even case i is doubled, and in the odd case j is increased by i. The final result is formula_22 where formula_23.
let exp_by_squaring x n =
 let rec _exp x n' a1 a2 =
 if n' = 0 then 1
 elif n' = 1 then a1*a2
 elif n'%2 = 0 then _exp x (n'/2) (a1*a1) a2
 else _exp x (n'-1) a1 (a1*a2)
 _exp x n x 1
The iterative version of the algorithm also uses a bounded auxiliary space, and is given by
 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
Computational complexity.
A brief analysis shows that such an algorithm uses formula_24 squarings and at most formula_24 multiplications, where formula_26 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""n" is given by
2"k"-ary method.
This algorithm calculates the value of "xn" 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:
 y := 1; i := l - 1
 while i ≥ 0 do
 (s, u) := f(ni)
 for j := 1 to k - s do
 y := y2 
 y := y * xu
 for j := 1 to s do
 y := y2
 i := i - 1
 return y
For optimal efficiency, "k" should be the smallest integer satisfying
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 and calculate 1, x3, x6, x12, x24, x48, x49, x98, x99, x198, x199, x398. 
But, we can also compute 1, x3, x6, x12, x24, x48, x96, x192, x199, x398, which saves one multiplication and amounts to evaluating (110 001 110)2
Here is the general algorithm:
Algorithm:
Algorithm:
 y := 1; i := l - 1
 while i > -1 do
 if ni = 0 then
 y := y2' i := i - 1
 else
 while ns = 0 do
 s := s + 1
 for h := 1 to i - s + 1 do
 y := y2
 u := (ni, ni-1, ..., ns)2
 y := y * xu
 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"0)2 with "n"k−1 = 1, we can compute "xn" as follows:
 x1 = x; x2 = x2
 for i = k - 2 to 0 do
 If ni = 0 then
 x2 = x1 * x2; x1 = x12
 else
 x1 = x1 * x2; x2 = x22
 return x1
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 "xn" 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
where formula_32 for all formula_33.
Let .
Then the algorithm uses the equality
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 ni = j then
 u = u × xbi
 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 .
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 formula_19 in group , where is a natural integer, whose algorithm is given below, is using the following equality recursively:
where formula_37.
In other words, a Euclidean division of the exponent by is used to return a quotient and a rest .
Given the base element in group , and the exponent formula_12 written as in Yao's method, the element formula_19 is calculated using formula_40 precomputed values formula_41 and then the algorithm below.
 Begin loop 
 
 Break loop 
 
 End loop;
The algorithm first finds the largest value among the and then the supremum within the set of .
Then it raises to the power , multiplies this value with , and then assigns the result of this computation and the value modulo .
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
The method works in every semigroup and is often used to compute powers of matrices.
For example, the evaluation of
would take a very long time and much storage space if the naïve method were used: compute 13789722341, 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 formula_42 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.
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 is "fast" or has been precomputed. For example, when computing , the binary method requires multiplications and squarings. However, one could perform squarings to get and then multiply by to obtain .
To this end we define the signed-digit representation of an integer in radix as 
"Signed binary representation" corresponds to the particular choice and formula_44. It is denoted by formula_45. There are several methods for computing this representation. The representation is not unique. For example, take : two distinct signed-binary representations are given by formula_46 and formula_47, where formula_48 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 formula_49 and denoted by formula_50. For example, the NAF representation of 478 is formula_51. This representation always has minimal Hamming weight. A simple algorithm to compute the NAF representation of a given integer formula_52 with formula_53 is the following:
 for to do
 
Another algorithm by Koyama and Tsuruoka does not require the condition that formula_54; it still minimizes the Hamming weight.
Alternatives and generalizations.
Exponentiation by squaring can be viewed as a suboptimal addition-chain exponentiation algorithm: it computes the exponent by 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.