This articlemay beconfusing or unclear to readers. Please helpclarify the article. There is a discussion about this ontalk:Exponentiation by squaring § Least significant bit is first.(May 2022) (Learn how and when to remove this message)

This articleneeds additional citations forverification. Please helpimprove this article byadding citations to reliable sources. Unsourced material may be challenged and removed. Find sources: "Exponentiation by squaring" – news ·newspapers ·books ·scholar ·JSTOR(February 2018) (Learn how and when to remove this message)

Inmathematics andcomputer programming,exponentiating by squaring is a general method for fast computation of largepositive integer powers of anumber, or more generally of an element of asemigroup, like apolynomial or asquare matrix. Some variants are commonly referred to assquare-and-multiply algorithms orbinary exponentiation. These can be of quite general use, for example inmodular arithmetic or powering of matrices. For semigroups for whichadditive notation is commonly used, likeelliptic curves used incryptography, this method is also referred to asdouble-and-add.

The method is based on the observation that, for any integer${\displaystyle n>0}$, one has

If the exponentn 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,$${\displaystyle x^n={\left(\frac{1}{x}\right)}^{-n}.}$$

Together, these may be implemented directly as the followingrecursive algorithm:

Inputs: a real numberx; an integernOutput:xnfunction exp_by_squaring(x,n)isifn < 0thenreturn exp_by_squaring(1 /x, −n)else ifn = 0thenreturn 1else ifn is eventhenreturn exp_by_squaring(x ×x,n / 2)else ifn is oddthenreturnx × exp_by_squaring(x ×x, (n − 1) / 2)end function

In each recursive call, the least-significant digit of thebinary representation ofn is removed. It follows that the number of recursive calls is${\displaystyle \lceil {\log }_{2}n\rceil ,}$ the number ofbits of the binary representation ofn. 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 ofn. This logarithmic number of operations is to be compared with the trivial algorithm which requiresn − 1 multiplications.

This algorithm is nottail-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.

The variants described in this section are based on the formula

If one applies recursively this formula, by starting withy = 1, one gets eventually an exponent equal to0, and the desired result is then the left factor.

This may be implemented as a tail-recursive function:

Functionexp_by_squaring(x,n)returnexp_by_squaring2(1,x,n)Functionexp_by_squaring2(y,x,n)ifn<0thenreturnexp_by_squaring2(y,1/x,-n);elseifn=0thenreturny;elseifniseventhenreturnexp_by_squaring2(y,x*x,n/2);elseifnisoddthenreturnexp_by_squaring2(x*y,x*x,(n-1)/2).

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

Functionexp_by_squaring_iterative(x,n)ifn<0thenx:=1/x;n:=-n;ifn=0thenreturn1y:=1;whilen>1doifnisoddtheny:=x*y;n:=n-1;x:=x*x;n:=n/2;returnx*y

The correctness of the algorithm results from the fact that${\displaystyle y{x}^n}$ is invariant during the computation; it is${\displaystyle 1\cdot x^n=x^n}$ at the beginning; and it is${\displaystyle y{x}^1=xy}$ 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.

A brief analysis shows that such an algorithm uses${\displaystyle \lfloor {\log }_{2}n\rfloor }$ squarings and at most${\displaystyle \lfloor {\log }_{2}n\rfloor }$ multiplications, where${\displaystyle \lfloor \rfloor }$ denotes thefloor function. More precisely, the number of multiplications is one less than the number of ones present in thebinary expansion ofn. Forn 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 twod-digit numbers is implemented in O(d^k) operations for some fixedk, then the complexity of computingxⁿ is given by

${\displaystyle \sum _{i=0}^{O(\log n)}(2^{i}O(\log x){)}^k=O((n\log x{)}^k).}$

This algorithm calculates the value ofxⁿ after expanding the exponent in base 2^k. It was first proposed byBrauer in 1939. In the algorithm below we make use of the following functionf(0) = (k, 0) andf(m) = (s, u), wherem =u·2^s withu odd.

Algorithm:

Input

An elementx ofG, a parameterk > 0, a non-negative integern = (n_l−1,n_l−2, ...,n₀)_2^k and the precomputed values${\displaystyle x^3,x^5,...,x^{2^k-1}}$.

Output

The elementxⁿ inG

y := 1; i := l - 1while i ≥ 0 do    (s, u) := f(ni)for j := 1to k - sdo        y := y2    y := y * xufor j := 1to sdo        y := y2    i := i - 1return y

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

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 elementx ofG, a non negative integern=(n_l−1,n_l−2, ...,n₀)₂, a parameterk > 0 and the pre-computed values${\displaystyle x^3,x^5,...,x^{2^k-1}}$.

Output

The elementxⁿ ∈G.

Algorithm:

y := 1; i := l - 1while i > -1doif ni = 0then        y := y2        i := i - 1else        s := max{i - k + 1, 0}while ns = 0do            s := s + 1[notes 1]for h := 1to i - s + 1do            y := y2        u := (ni, ni-1, ..., ns)2        y := y * xu        i := s - 1return y

Many algorithms for exponentiation do not provide defence againstside-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 manypublic-key cryptosystems. A technique called "Montgomery's ladder"^[2] addresses this concern.

Given thebinary expansion of a positive, non-zero integern = (n_k−1...n₀)₂ withn_k−1 = 1, we can computexⁿ as follows:

x1 = x; x2 = x2for i = k - 2 to 0doif ni = 0then        x2 = x1 * x2; x1 = x12else        x1 = x1 * x2; x2 = x22return 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 cachetiming 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.^[3]

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

Yao's method is orthogonal to the2^k-ary method where the exponent is expanded in radixb = 2^k and the computation is as performed in the algorithm above. Letn,n_i,b, andb_i be integers.

Let the exponentn be written as

${\displaystyle n=\sum _{i=0}^{w-1}{n}_i{b}_i,}$

where for all${\displaystyle i\in [0,w-1]}$.

Letx_i =x^b_i.

Then the algorithm uses the equality

${\displaystyle x^n=\prod _{i=0}^{w-1}{x}_i^{n_i}=\prod _{j=1}^{h-1}[\prod _{n_i=j}{x}_i{]}^j.}$

Given the elementx ofG, and the exponentn written in the above form, along with the precomputed valuesx^b₀...x^b_w−1, the elementxⁿ is calculated using the algorithm below:

y = 1, u = 1, j = h - 1while j > 0dofor i = 0to w - 1doif ni = jthen            u = u × xbi    y = y × u    j = j - 1return y

If we seth = 2^k andb_i =hⁱ, then then_i values are simply the digits ofn in baseh. Yao's method collects inu first thosex_i that appear to the highest power⁠${\displaystyle h-1}$⁠; in the next round those with power⁠${\displaystyle h-2}$⁠ are collected inu as well etc. The variabley is multiplied⁠${\displaystyle h-1}$⁠ times with the initialu,⁠${\displaystyle h-2}$⁠ times with the next highest powers, and so on.The algorithm uses⁠${\displaystyle w+h-2}$⁠ multiplications, and⁠${\displaystyle w+1}$⁠ elements must be stored to computexⁿ.^[1]

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

This method for computing${\displaystyle x^n}$ in groupG, wheren is a natural integer, whose algorithm is given below, is using the following equality recursively:

${\displaystyle x_0^{n_0}\cdot x_1^{n_1}={\left(x_0\cdot x_1^q\right)}^{n_0}\cdot x_1^{n_1\mathrm{mod}{n}_0},}$

where${\displaystyle q=\lfloor \frac{n_1}{n_0}\rfloor }$.In other words, a Euclidean division of the exponentn₁ byn₀ is used to return a quotientq and a restn₁ modn₀.

Given the base elementx in groupG, and the exponent${\displaystyle n}$ written as in Yao's method, the element${\displaystyle x^n}$ is calculated using${\displaystyle l}$ precomputed values${\displaystyle x^{b_0},...,x^{b_{l_i}}}$ and then the algorithm below.

Begin loopFind${\displaystyle M\in [0,l-1]}$,such that${\displaystyle \mathrm{\forall }i\in [0,l-1],n_M\ge n_i}$.Find${\displaystyle N\in ([0,l-1]-M)}$,such that${\displaystyle \mathrm{\forall }i\in ([0,l-1]-M),n_N\ge n_i}$.Break loopif${\displaystyle n_N=0}$.Let${\displaystyle q=\lfloor n_M/n_N\rfloor }$,and thenlet${\displaystyle n_N=(n_M{modn}_N)}$.Compute recursively${\displaystyle x_M^q}$,and thenlet${\displaystyle x_N=x_N\cdot x_M^q}$.End loop;Return${\displaystyle x^n=x_M^{n_M}}$.

The algorithm first finds the largest value among then_i and then thesupremum within the set of{n_i \i ≠M }.Then it raisesx_M to the powerq, multiplies this value withx_N, and then assignsx_N the result of this computation andn_M the valuen_M modulon_N.

The approach also works withsemigroups that are not ofcharacteristic zero, for example allowing fast computation of largeexponents modulo a number. Especially incryptography, it is useful to compute powers in aring ofintegers moduloq. For example, the evaluation of

13789⁷²²³⁴¹ (mod 2345) = 2029

would take a very long time and much storage space if the naïve method of computing13789⁷²²³⁴¹ and then taking theremainder when divided by 2345 were used. Even using a more effective method will take a long time: square 13789, take the remainder when divided by 2345, multiply theresult by 13789, and so on.

Applying aboveexp-by-squaring algorithm, with "*" interpreted asx *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. Generally, any of these approaches will take fewer than2log₂(722340) ≤ 40 modular multiplications.

The approach can also be used to compute integer powers in agroup, using either of the rules

Power(x, −n) = Power(x⁻¹,n),

Power(x, −n) = (Power(x,n))⁻¹.

The approach also works innon-commutative semigroups and is often used to compute powers ofmatrices.

More generally, the approach works with positive integer exponents in everymagma for which thebinary operation ispower associative.

In certain computations it may be more efficient to allow negative coefficients and hence use the inverse of the base, provided inversion inG is "fast" or has been precomputed. For example, when computingx^2k−1, the binary method requiresk−1 multiplications andk−1 squarings. However, one could performk squarings to getx^2k and then multiply byx⁻¹ to obtainx^2k−1.

To this end we define thesigned-digit representation of an integern in radixb as

Signed binary representation corresponds to the particular choiceb = 2 and${\displaystyle n_i\in \{-1,0,1\}}$. It is denoted by${\displaystyle (n_{l-1}\dots n_0{)}_s}$. There are several methods for computing this representation. The representation is not unique. For example, taken = 478: two distinct signed-binary representations are given by${\displaystyle (10\overline{1}1100\overline{1}10{)}_s}$ and${\displaystyle (100\overline{1}1000\overline{1}0{)}_s}$, where${\displaystyle \overline{1}}$ is used to denote−1. Since the binary method computes a multiplication for every non-zero entry in the base-2 representation ofn, we are interested in finding the signed-binary representation with the smallest number of non-zero entries, that is, the one withminimalHamming weight. One method of doing this is to compute the representation innon-adjacent form, or NAF for short, which is one that satisfies${\displaystyle n_i{n}_{i+1}=0\text{ for all }i⩾0}$ and denoted by${\displaystyle (n_{l-1}\dots n_0{)}_{\text{NAF}}}$. For example, the NAF representation of 478 is${\displaystyle (1000\overline{1}000\overline{1}0{)}_{\text{NAF}}}$. This representation always has minimal Hamming weight. A simple algorithm to compute the NAF representation of a given integer${\displaystyle n=(n_l{n}_{l-1}\dots n_0{)}_2}$ with${\displaystyle n_l=n_{l-1}=0}$ is the following:

${\displaystyle c_0=0}$fori = 0 tol − 1 do${\displaystyle c_{i+1}=\lfloor \frac{1}{2}(c_i+n_i+n_{i+1})\rfloor }$${\displaystyle n_i^{\prime }=c_i+n_i-2c_{i+1}}$return${\displaystyle (n_{l-1}^{\prime }\dots n_0^{\prime }{)}_{\text{NAF}}}$

Another algorithm by Koyama and Tsuruoka does not require the condition that${\displaystyle n_i=n_{i+1}=0}$; it still minimizes the Hamming weight.

Exponentiation by squaring can be viewed as a suboptimaladdition-chain exponentiation algorithm: it computes the exponent by anaddition chain consisting of repeated exponent doublings (squarings) and/or incrementing exponents byone (multiplying byx) only. More generally, if one allowsany previously computed exponents to be summed (by multiplying those powers ofx), one can sometimes perform the exponentiation using fewer multiplications (but typically using more memory). The smallest power where this occurs is forn = 15:

${\displaystyle x^{15}=x\times (x\times [x\times x^2{]}^2{)}^2}$ (squaring, 6 multiplies),

${\displaystyle x^{15}=x^3\times ([x^3{]}^2{)}^2}$ (optimal addition chain, 5 multiplies ifx³ is re-used).

In general, finding theoptimal addition chain for a given exponent is a hard problem, for which no efficient algorithms are known, so optimal chains are typically used for small exponents only (e.g. incompilers where the chains for small powers have been pre-tabulated). However, there are a number ofheuristic 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Θ(logn), so these algorithms improve asymptotically upon exponentiation by squaring by only a constant factor at best.

• Modular exponentiation
• Vectorial addition chain
• Montgomery modular multiplication
• Non-adjacent form
• Addition chain

• Exponentials
• Computer arithmetic algorithms
• Computer arithmetic

• Articles with short description
• Short description is different from Wikidata
• Wikipedia articles needing clarification from May 2022
• All Wikipedia articles needing clarification
• Articles needing additional references from February 2018
• All articles needing additional references