Exponentiation using Division Chains
Colin D. Walt er, Member, IEEE
Computation Department, U.M.I.S.T.,
PO Box 88, Sackville Street, Manchester M60 1QD, U.K.
C.Walter@umist.ac.uk

Abstract
Exponentiation may be performed faster than the
traditional square and multiply method by iteratively
reducing the exponent modulo numbers which, as
exponents themselves, require few multiplications. This
mainly includes those with few non-zero bits. For a
suitable choice of such divisors, the resulting mixedbasis representation of the exponent reduces the
expected number of non-squaring multiplications by overhalf at the cost of a single extra register. Pre-processing
effort depends entirely on the exponent and can be kept
down to the work saved in a single exponentiation.
Moreover, no pre-computed look-up tables are required,
so the method is especially applicable where space is at
a premium. In particular, it out performs the instance ofthe m-ary method which uses the same space. However,for 512-bit exponents, it beats every instance of the mary method, achieving well under 635 multiplications on
average. Both hardware and software implementationsof the RSA crypto-system can benefit from this algorithm.
Key Words Modular exponentiation, bit recoding, RSA
cryptosystem, addition chains, m-ary method, mixedbasis arithmetic, radix representation.

1

Introduction

Fast exponentiation is becoming increasingly
important with the widening use of encryption. Whereas
the most startling improvements in speed are achieved
through the use of dedicated hardware for multiplication,some small gains can also be made through a good
choice of algorithm for organising the order ofmultiplications at either the hardware or software level.
The expected number of multiplications (including
squarings) to compute Ae in the traditional way, usingsquare and multiply, is approximately 23 log 2 e [14],[15]. The pre-calculation or intermediate calculation of

http://www.co.umist.ac.uk/

auxiliary powers of A enables some unnecessary
repetition of work to be avoided. There are several
methods described below for reducing the average
number of operations to about 43 log 2 e using storage for
only one or two extra powers of A, but the variance of the
distribution is high: there is a very real chance of stillrequiring 23 log 2 e or even more operations. With morestorage, the coefficient43

can be further reduced [22],but the lower limit is still above 1 [11], [14] sincelog 2 e squarings are certainly necessary to generate anumber as large as Ae. However, fast registers arenormally very limited in number, so the use of too much
more storage is actually likely to retard theexponentiation: the cost of communication with slower
memory may easily outweigh any minimal reduction in
the coefficient.
This paper describes another technique which uses
only one register more than the minimum for the
standard square and multiply method. It achieves anaverage coefficient below 45 , and is consistently close to
its average. Some pre-computation is required, but only
enough to establish a suitable sequence of squarings andmultiplications. Thus the pre-computation need not be
repeated in a sequence of exponentiations using the same
exponent. A simple version of the technique achieves
5with a pre-computation effort equivalent to only a few
4
multiplications of integers the size of the exponent. So
the overall gain is small for a single exponentiation, butwhen the same key is used repeatedly as in the RSA
cryptosystem, the pre-computation cost can be ignored
and the full reduction claimed. The method is very
adaptable, with better results easily obtained from greater
effort or from allowing more register space. It also works
equally well for any group, so that the speed-ups work as
well for an elliptic curve group as for an integer residue
class group.

IEEE TRANSACTIONS ON COMPUTERS, Vol. 47, No. 7, July 1997, pp. 757-765.

2

one of m registers on each iteration has a critical path
length of order O(log 2 m) . Thus, without considerable
care, the complexity of the non-squaring multiplications
for this method might still be essentially O(log 2 e) .

Notation and Literature Review

On a sequential machine, any exponentiation by e
can be described by an addition chain a0, a1, a2, a3,..., an,where a0 = 1, an = e and, for each i > 0, ai = aj + ak
for some j, k < i [14]. The ith multiplication performedaaHalf the memory can be saved by pre-computing
only Ai for odd i < m [22]. The exponent is recoded asais
A i = A j×A k and exponentiation by e takes nmultiplications. Storage requirements for any chain canbe worked out easily from the sequence itself, although if
there is a choice of j and k for any i then the minimum
storage might not be clear. For given small exponents e,
the minimal number of multiplications can be found by asearch of all addition chains for e. As an NP-hardproblem [8], however, such a search is impractical for the
typical decryption keys e found in RSA cryptography.
Nevertheless, here we present a technique which is easily
tailored to the available time and space resources andultimately would yield a best addition chain.

e = r0 + m0 (r1 + m1 (r2 + m2 (...(rn −1 + mn −1rn )...))) (1)where mi = m is chosen whenever it makes ri odd, andotherwise mi = 2 is chosen with ri = 0. Then Ae isevaluated starting at the innermost bracketed expression:
A e = ((...(( A rn ) mn −1 A rn −1 ) mn − 2 ...) m1 A r1 ) m0 A r0

This requires m/2 registers besides the partial result,about log 2 e − 12 log 2 m + 1 squarings if m is a power of2, and about log 2 e / (log 2 m + 1) other multiplications
on average besides any used for the pre-computations.
(For i < n there are on average as many cases of mi = m
as mi = 2.) Taking m = 2 gives the binary method above
with coefficient 23 for log 2 e in the total number ofSuppose e = ∑nj = 0 r j 2 j is the binary representationof e. The standard method of square and multiply can be
performed by scanning the bits of e in either direction.
First, a Horner-style evaluationA e = ((...(( A rn ) 2 A rn −1 ) 2 ...) 2 A r1 ) 2 A r0

multiplications. Taking m = 4 means a third register,which holds A3, and this reduces the coefficient to 43 .

corresponds to an addition chain with ai+1 = 2ai (square)or ai+1 = ai + a0 (multiply). This requires just 2 storageaaregisters, containing A = A 0 and the partial result A i
respectively. Alternatively, and dually, rather thansquaring the partial result and multiplying in A as
required, A can be repeatedly squared and the resulting
power multiplied into the partial result when needed:
0

1

2

(2)

If the inverse A–1 exists and can be calculated
cheaply, then any sequence of 1s can be replaced by thesequence 10...0 1 as follows. Starting from the leastsignificant bit of e, whenever two adjacent 1s areencountered with no carry from lower down, the lower 1
is replaced by 1 and a carry of 1 generated to the upper
1. Similarly, for a carry of 1 into 10, the lower 0 is
replaced by 1 with a carry up to the 1. Otherwise carries
propagate up as usual. On average the carry is 1 in 50%
of cases, 2/3rds of these are into a following 1. Similarly,
a carry of 0 is to a digit 0 in 2/3rds of cases. Thus, after
this modified Booth recoding, a 0 digit occurs on averagein 2/3rds of all cases. So again about 43 log 2 e squaringsnA e = ( A 2 ) r0 ( A 2 ) r1 ( A 2 ) r2 ...( A 2 ) rniHere the first of the two registers now contains A2 for i =0,1,..., n. For natural number arithmetic this requires
larger registers than Horner’s method, but there is nodifference for finite rings, such as the residues of RSA
cryptosystems, or for real approximations. The numberof multiplications, excluding squarings, is one less thanthe Hamming weight of e, i.e. one less than the numberof non-zero bits in e; on average log 2 e / 2 .

or multiplications by A or A–1 are required, with one
multiplication for each non-zero digit [26]. Storage is 3
registers, one for each of A, A–1 and the partial result.
This method can be combined with the m-ary method [9],[10], [15]. Further similar methods include [12].

By expressing the exponent using the radix minstead of 2 and pre-computing the powers Ai for i = 1, 2,..., m–1 we obtain the m-ary method [14]. This follows
the Horner style evaluation above, requiring repeatedly
raising the partial result to the mth power and thenmultiplying by an Ai . It is usually convenient to pick m
as a power of 2. Then the number of squarings isroughly the same as before, but the number of othermultiplications excluding pre-computations reduces tolog m e(m − 1) / m on average. This is good for larger

The m-ary method has been developed in several
important ways. One of these uses the ideas of datacompression: common subsequences of bits are identified
and the corresponding powers of A computed once andstored for use each time the bit sequence is encountered.
In effect, this windowing technique corresponds toallowing the mi in (1) to be any power of 2 [1], [2], [16].
Another method employs vector addition chains so thatwork involved in computing several different powers of A
can be shared wherever possible [2],[18],[19],[20],[23].
Similar sharing is also considered in [17], [24], [25].

m, but storage requirements quickly become prohibitive
for RSA applications. Moreover, the hardware to select2Colin D. Walter

Exponentiation Using Division Chains

In many of these methods, the trade-off between
space and time is poor when e is static and A variable. In
hardware implementations, a speed-up of less than twofold, as here, should require less than twice the area.
Since the total hardware requirement for a digit parallel
modular multiplier is equivalent to only a small numberof registers [21], few of these techniques are likely to beof much value to hardware designers. However, we willproduce more than half of the maximum possible savings
by using only 1 more register than for the square andmultiply method.

3

Mixed Basis Representations.

In a mixed basis number representation, an exponent
e is written as a digit linear combinatione = r0 b0 + r1b1 +...+rn bn
of basis elements bi . Normally, bi is a multiple of bi–1,
with bi = mi in a standard radix m representation. The
divisibility property enables the nested representation of enecessary in the m-ary method. In the variation where
half the memory is saved, the basis elements are bi =m0m1...mi–1 = bi–1mi–1 with mi–1 = 2 or m. Dimitrov et al[6] consider the particular case of this representation
with m = 3, illustrating the value of m not being a powerof 2. There the digits ri are 0 or 1 and basis elements bi
are of the form 2j3k. This leads to 1.42361×log2e
multiplications on average. They suggest that each mi
could be any prime up to some bound m. Then the powers
Ai could be pre-computed for i < m and Ae calculated
following the m-ary method, that is, with Horner styleevaluation. In fact, the algorithm presented in [6] for m= 3 processes the digits in the opposite order; the cost is
the same for that case. If many of the possible bi th
powers of A are precomputed, then this latter order canlead to considerable time savings when a number of
exponentiations with fixed A and variable e are involved[4], [7]. However, this is a different problem and thespace cost may be large.
Here we retain the key divisibility property of thebasis elements, namely bi = bi–1mi–1, and develop
algorithms for making good choices of the mi . The
divisibility property enables the digits to be processed ineither order. Where space is expensive, using Horner’sorder of evaluation without precomputed powers leads toduplication of effort when digits are equal since then theri th power of A must be re-calculated. We concentrate on
scanning the digits in the reverse order: the bi th power ofA will be computed dynamically from the bi–1st and held
in a single register, just as the 2i th power is held in thecorresponding binary method. At the same time, theri–1bi–1st power will be formed and multiplied into thepartial result.

4

The Division Chain Method

The new means of reducing the number ofmultiplications used in exponentiation arises from theiterative application of a decomposition e = me’ + r
where r is usually the least non-negative residue of e
modulo m. At each repetition the divisor m is selected byreference to a pre-determined set of inexpensive pairs(m,r) and the powers Am and Ar are computed. Since Aesatisfies the relationship
Ae = ( Am )e’ Ar(3)r
it suffices to multiply A into a partial product registerand re-apply the process to the remaining problem ofraising Am to the power e’. Including all residues forsome divisor such as 2, 3 or 5 guarantees each step ispossible and that termination occurs. When e’ = 0 has
been processed the partial product register contains therequired output. So, if m0, m1, ..., mn is the list ofdivisors which this generates, and r0, r1, ..., rn are theassociated remainders then, as in (2),A e = A r0 + m0 ( r1 +m1 (...+ mn −1 ( rn −1 + mn −1rn )...))and evaluation is performed by processing the exponent
expression from left to right. Note that the divisors forthe algorithm can be chosen on the fly as the
exponentiation is performed or in advance whenever the
exponent is known.

We assume that it is known how best to calculate Am
and Ar. So divisor / residue pairs (m,r) must be selectedfrom a set for which this information is known. Then ateach step the cheapest such decomposition can be chosenfrom this fixed set of pairs. A sequence of pairs (m,r)
used to direct an exponentiation in this way will be calleda division chain by analogy with the addition chain
description. Although it is usually the case, we see laterthat the residue need not always be the least non-negative
one. However, if no alternative choice of residue is ever
permitted, then the sequence of divisors suffices todetermine the division chain.
When m = 2 at each step, this technique is just thestandard square and multiply method described above. If
r = 0 at each step it becomes the factor method [14]. So
the division chain method here generalises both of these.
Moreover, if mi = m always, it is the dual of the m-ary
method, processing digits in the opposite order. Bos andCoster [2] provide some heuristics for their
MakeSequence method of constructing addition chains.
This method shows how these can be realised.

IEEE TRANSACTIONS ON COMPUTERS, Vol. 47, No. 7, July 1997, pp. 757-765.
If F(e) is the minimum number of multiplications
used to compute Ae by any method, then thedecomposition (3) shows F(e) ≤ F(e’)+F(m)+F(r)+1.
Inequality will arise if either r is 0 or somemultiplications which are used to form Am may also be
used in constructing Ar. However, F(e) ≈ F(e’)+F(m)because e ≈ me’ and F(x) is always of order log(x). So
the main step of the algorithm can only be useful if Am
and Ar are both cheap to compute, with much of thiswork in common. In particular, this is the case when m= 2n+1 and r = 0, 2n–1+1 or 2i for some i ≤ n. Then theformation of Ar is a by-product of computing Am. In
general, useful pairs (m,r) normally have the property
that r lies in an addition chain for m of minimal lengthor, at worst, r is the sum of two or three members of sucha chain.
Consequently, if f(e) is the number ofmultiplications this method yields, and there aresufficient such pairs (m,r), then we will have f(e) ≈f(e’)+f(m) at each step and can expect to obtain areasonably efficient scheme for exponentiation.

with r chosen as the least non-negative residue modulo
m. Let f (e) = c log 2 e be the average number ofsquarings plus multiplications expected for a random
exponent e. Then f(e) can be expressed fairly accurately
as a sum of terms p×(v+f(e/m)), one for each case (m,r),where p is the probability of the case arising, m is theassociated divisor, and v is the number of multiplications
required to form Am and Ar and multiply Ar into theresult. For the case e ≡ 5 mod 9, we can use the addition
chain (1, 2, 4, 5, 9) to show that 4 multiplications yieldboth Am and Ar, giving v = 5. The other values of v areclear.
emod 18m

0,6,12
2,8,14
4,10,16
3
9
15
1
7
13
5
11
17

2
2
2
3
3
3
9
3
3
9
9
9

Rel
freq

4
10
10
 
6
 
 
9
 
 
9
 
Freq sums: 48

With suitable restrictions on the divisors m andresidues r, only three storage registers are needed,although more can be used. One register holds thepartial result which is the product of all the Ar fromprevious steps. The other two are used to form Am using
an addition chain that requires only two values to be
kept. When the components required for Ar are formed,
they are multiplied into the partial result register and so
do not interfere with the calculation of Am.

Distrib of new e mod 6
0
1
2
3
4
5
2
1
1
-

5
2
1
1

5
1
3
1
-

2
2
1
1

5
1
3
1
-

5
2
1
1

4

9

10

6

10

9

Table 1. Frequency counts for the Markov process.

349

As an illustration, we show how to compute A
starting with the divisor 17. First, A349 = (A17)20A9
requires computing B = A17 and A9. These are found in 5
multiplications using the addition chain (1, 2, 4, 8, 9,17), and A9 is placed in the partial product register. This
leaves B20 to be calculated and multiplied into theaccumulating product. Using the divisor 4 next, B20 =(B4)5. We obtain C = B4 with 2 squarings and must thencompute C5 for multiplying into the partial product.
Using the divisor 4 again, C5 = (C4)1C1. The C1 ismultiplied into the partial product, then C4 is computedwith two squarings and the result finally multiplied into
the partial product to yield A349. So exponentiation by
349 can be done this way with 11 multiplications rather
than the 13 required by the binary and A, A3 methods.

5

The values of the probabilities p are less obvious
since application of the method causes the residue classes
to be no longer uniformly distributed. Action depends
entirely on the residue of e modulo the lowest common
multiple of the divisors, namely 18, and not on any
previous exponent.
The sequence of residues for
successive exponents forms a Markov chain for which thestochastic transition matrix P = (pij) is easy to construct:
pij is the probability of obtaining the new exponent e’ ≡ j
mod 18 from an exponent e ≡ i mod 18. For example, e= 18k+7 uses m = 3, yielding e’ = 6k+2 ≡ 2, 8 or 14 mod18. So p7,2 = p7,8 = p7,14 = 1/3. The probabilities pi of
exponents which are i mod 18 eventually stabilise: therow vector p = (pi ) satisfies p = pP and can be calculated
easily. It turns out that pi only depends on i mod 6, with
p0 = p6 = p12 = 1/36, etc. The relative frequencies of the
18 classes are given in Table 1, which can readily be
used to check that these do indeed yield the equilibrium
values.

Calculation of the Coefficient

We will demonstrate the value of the method with asmall set of divisors. Suppose the strategy is as follows:Ife ≡ 0 mod 2
elseif e ≡ 0 mod 3
elseif e ≡ 1,2,5,8 mod 9else{e ≡ 4,7 mod 9}

So the 4 cases of the algorithm have probabilities 1/2,
1/8, 1/4 and 1/8 respectively. In consequence,f(e) = (1/2)(1+f(e/2)) + (1/8)(2+f(e/3))+ (1/4)(5+f(e/9)) + (1/8)(3+f(e/3))

then m = 2
then m = 3
then m = 9then m = 3

4

Colin D. Walter
Then using f(e/m) ≈ f(e)–f(m) ≈ f(e)–clog2m we obtain
0 = (1/2)(1–clog2 2) + (1/8)(2–clog23)+ (1/4)(5–2clog23) + (1/8)(3–clog23)
So c ≈ 1.4064 (cf. [6] with c ≈ 1.42361 for a slightly
simpler example).
Better values for the coefficient c arise from the useof more divisors. Consider:e ≡ 0 mod 33
then m = 33e ≡ 0 mod 17
then m = 17e ≡ 0 mod 2
then m = 2e ≡ 0 mod 5
then m = 5e ≡ 0 mod 3
then m = 3
e ≡ 1,2,4,8,16,17,32 mod 33
then m = 33
elseif e ≡ 1,2,4,8,9,16 mod 17
then m = 17else case e mod 90 of
1,17,77,29,43,47,67,83: m = 3
7,11,13,31,41,49,61,71,79,89: m = 2
19,23,37,53,59,73: m = 9end.

If
elseifelseifelseifelseifelseifAgain, choose r minimally. The more complex final caseis engineered to enable useful factors such as 2, 3, 5 or 9
to be picked up on the next iteration. As before theresidue classes modulo 90×11×17 are not equi-probable.
The transition matrix P has grown to around 228 entries
and made the direct solution of p = pP infeasible.
However, if i is the suitably scaled all 1s row vector thenp = limn→∞ iPn. In practice convergence is fairly rapid.
The Euclidean distance between successive values of p isalmost halved at each iteration, and so another decimal
place of c is established on every fourth iteration. After
just 7 iterations we obtain the correctly rounded c =1.3566. When the divisor 65 is included in a similar
fashion to 33, c is reduced to 1.343.
The time complexity of finding p is driven by the
number of non-zero entries in P. The matrix has size l×l
where l is the lowest common multiple of the moduli
used, and it has not far short of ln non-zero entries where
n is the number of different divisor/residue pairsconsidered. Thus, including more than a few small
divisors makes the vector p too big and the computation
too time consuming for a direct solution, and it forces astatistical approach to estimating c.
To achieve a coefficient noticeably less than the 43of the A, A3 method, take the following divisors:
2, 3, 5, 17, 33, 49, 65, 97
129, 257, 513, 1025

and allow residues which are either 0, or are in aminimal addition chain which uses just two memory
fields, or are the sum of two residues in such a chain.
(Details are contained in the Appendix.) At each
iteration simply choose the divisor m for which the ratio

Exponentiation Using Division Chainsv/log2m is least, where v is the number of multiplications
associated with the residue which occurs. Then, on
average, random 512-bit integers require 672.15
multiplications – almost 1.5% fewer than the average of
682 for the A, A3 method.
How much work is involved in choosing the division
chain for a given exponent and a given such strategy? A
small, variable but bounded number of modular divisions
are performed on the exponent at each step to decide thenext divisor. Each of these divisions requires at most asmall constant times more work to perform than thedivision which is actually chosen. Thus choosing thedivisor m at each step requires on average a constant, sayc', times more work than simply dividing the remaining
exponent by m. Overall, c' times more work is done tofind the division chain than is used in dividing theoriginal exponent by the chain of divisors. However,dividing the exponent by the chain of divisors requires
the same order of work as finding the product of two
integers the same size as the exponent. So generating thedivision chain has the same time complexity as forming
the product of two numbers with the size of the exponent.
In the RSA cryptosystem the exponent e fordecrypting is usually of the same length as the message A
requiring decryption. So choosing the division chain willthen be equivalent to the cost of a fixed number ofmultiplications, say c". The total work of exponentiation
will then be about c"+c log 2 e multiplications. Thus, fora sequence of exponentiations with the same key e, this ischeaper than the A, A3 method as long as c < 43 and theexponent e is large enough.

6

Choice and Ordering of Divisors
From the last section the formula for c isc=∑i pi vi

∑i pi log 2 mi

where the sum extends over cases i, which occur withprobability pi , have divisor mi and require vi
multiplications. Thus, if the relative probabilities ofthese cases remain essentially unchanged, it is worth
adding any new pair (m,r) with a ratio v/log2m which is
better than the current value of c. In particular, in
addition to using numbers with the form 2n+1 (2 bit
numbers), we might also consider divisors m ofHamming weight 3 (i.e. 3 non-zero bits) with residues r
which are 0 or are generated en route to m. Then aminimal addition chain using just two registers will
require n+2 multiplications where the highest power of 2
in m has order n. Useful examples therefore include 49
and 97, both with ratios under 45 for these residues.

IEEE TRANSACTIONS ON COMPUTERS, Vol. 47, No. 7, July 1997, pp. 757-765.
It is a matter of only a few minutes computing to
generate all optimal and near optimal addition chains for
all potentially useful divisors up to 10 or so bits in length,say (see Appendix), and hence establish the smallest
number of extra additions which will generate eachresidue. Here we should relax the apparent restriction r< m since, to minimise multiplications, it may benecessary to go beyond the least non-negative residue.
Thus 11 mod 13 is only obtained with a minimal numberof multiplications within the suggested memory
restrictions by allowing r = 12+12. However, if both r
and r+m are equally cheap to compute, it may makesense to select the one which makes the next exponent e’even so that the beneficial (2,0) could then be applied.

7

Algorithms for Large Exponents

Most of the well-known algorithms described in
Section 2, such as square and multiply, prescribe a single
course of action, as do the strategies described here so far
for division chains. In all cases the high variance means
frequent poor results. However, for division chainsconsiderable choice is possible. Fixing a specific orderfor trying divisors usually yields a sub-optimal method.
A better coefficient c is obtained when the best chain is
selected from all those generated by extending partiallyconstructed chains in every possible way. This also
reduces the variance and hence gives a good value morereliably. Unfortunately, it is not feasible for a large
exponent – there are too many combinations to evaluate
all possibilities. The search space must be reduced.

Using all m up to over 3×210, optimal division chainswere constructed sequentially for e up to over 220. The
frequencies of pairs (m,r) were recorded, and thisconfirmed the relative uselessness of pairs with a poor
ratio v/log2m unless m was small. It also showed thatcomposite m in this range are often superfluous sincefrequently no more multiplications are required thanwhen their factors are used as divisors instead.

Suppose the division chain (m0,r0), (m1,r1), (m2,r2),..., (mk–1,rk–1) reduces the exponentiation problem fromthe power e to the power e’. Then e = me’ + r where m =m0m1m2...mk–1 and (normally) r < m. The number ofmultiplications v associated with the chain is the sum of
the numbers vi associated with each pair (mi ,ri ), namely
the number of multiplications required to form the mi andri th powers and multiply the ri th power into the existing
partial product. If e’ is still large compared to m then e
is essentially reduced by a factor m for the cost of v
multiplications. Thus, given a number of such chains,
the best one to choose is the one for which v/log2m or,
better, v–clog2m is minimal. This process can be
repeated to reduce e’ in its turn. Eventually e’ becomes
small enough for the value of r to affect the choice of
chain noticeably (say e’ < m), and then a table might beused to complete the division chain optimally.

The obvious order in which to prioritise thedivisor/residue pairs for choice is in increasing cost per
bit, i.e. following the order of the ratios v/log2m.
However, larger divisors have a longer lasting effect thansmaller ones, and so they are better (resp. worse) choices
if they have similar but below (resp. above) average
ratios. To be more precise, if m is applied to the
exponent e, then we can expect v + clog2(e/m) = v –clog2(m) + clog2(e) multiplications on average. This isminimised if a divisor is chosen for which cost function
v–clog2m is minimal.

This yields an algorithm for large exponents. Fix asuitable length k for the division chain segments to beconsidered. The larger the choice of k, the longer thealgorithm takes to complete, but the better the result.
Now repeatedly perform the following:Thus, the best order for selecting pairs (m,r) should
be close to that determined by the values ofv – c log2m
not that of the ratios v/log2m. These criteria are called
the difference and ratio tests respectively. The former,
while slightly better, does need a good approximation toc in advance, whereas the latter requires no suchinformation and so could be used to generate theapproximation to c for the difference test.

i)generate all reasonably priced chains of k divisors
to reduce the current value of the exponent;ii) select the cheapest one under the chosen criterion;iii) apply this sub-chain to reduce the exponent.

The iterative procedure should terminate when thecurrent value of the exponent becomes less than theupper limit of a pre-calculated table of optimal chains forsmall exponents. Of course, if in the last few iterations kdivisors were to reduce the exponent to close to 0 so thatthe costing criterion becomes inaccurate, then theoffending chains can be curtailed earlier, say at the point
where the exponent falls into the range of the table.
Finally, with the exponent in the range of the table, thebest way of completing the chain can be looked up.

Using the better difference test explains the order inthe second example of Section 5 where (2,0) is not thebest first choice, as one might expect. Note that insituations such as the last case of that example, any
divisor that will be picked up automatically for the next
iteration needs to be taken into account when using either
cost function to assess the relative merits of different
divisors.

6Colin D. Walter
To cost this, suppose S is the sum over all divisors of
the probability of the divisor giving an acceptable
residue. So S will be at most the number of divisors
being used. For a better estimate, assuming a close touniform distribution on residues of exponents modulo
each divisor, S is approximately ∑i ni / mi where ni is
the number of residues associated with the divisor mi .
Then, to extend any division chain by one divisor, there
are roughly S possible choices (assuming successive
choices are mostly independent).
Next suppose that L is the result of averaging thelogarithms of the divisors used in an optimal sub-chain,weighted according to their frequencies, i.e. the average
number of bits by which a divisor in such a chain reduces
the exponent. Assume finally that T is the logarithm of
the maximum exponent in the look-up table. Then, for
exponents of N bits, the above algorithm will require
about (N–T)/Lk iterations, each of which generates about
Sk sub-chains. So the work involved in choosing adivision chain is roughly proportional toSk(N–T)D/Lk
where D is the effort required to perform a single
division of the exponent by a divisor. This shows that
selecting a division chain for very large exponents is notmuch more difficult to deal with than one for smaller
exponents; the work is proportional to the square of thenumber of bits.
For good performance, k must not be too small sincethen the sub-chains would not be particularly good on
average. Hence the best savings in time are made by
reducing S through a sensible choice of acceptable pairs.
An obvious choice is the set of pairs used most frequently
over a large range of optimal chains. Furthermore,
picking large k may be a waste of effort since the average
reduction in numbers of multiplications declines
exponentially with k (see Section 9).
Many variations in the algorithm are clearlypossible. For example, the subchains could be boundedby limiting the divisor product rather than by k, or k
might be varied when there is no good subchain of a
specific length at some point. Also, only the initial few
divisors from optimal sub-chains might be selected ateach iteration and the test for optimality (which contains
an approximation to c) might be varied dynamically. If a
good chain is still not obtained, choosing a different firstdivisor or using all initial subchains is also possible.

8

Test Results

We next consider test results from a specificimplementation of the algorithm in order to establish that

Exponentiation Using Division Chainsa coefficient of 45 is generally obtainable for exponentsof the size used in, say, the RSA cryptosystem. None of
the improvements suggested in the last paragraph were
included for the tabulated results.
First, for divisors m up to 28 the minimal addition
chains using only two registers were generated and atable was constructed of optimal division chains for
exponents up to over 220. In any subrange the method
was found to use fewer than 54 log 2 e multiplications andsquarings on average. In particular, between 215 and 216,
apart from a few short chains such as those for 215 itself
and 215+2a with a < 15, almost 90% require 19 or 20
multiplications, under 1.5% require 21, and none require
more. Hence the variation in numbers of multiplications
is very small, unlike for the binary or A, A3 methods.
Essentially no exponents require too many
multiplications, but there is an increasing tail of
exponents which can be dealt with using fewer
multiplications than expected. It is this that makes
searches for better division chains worthwhile.
Next, over ranges of at least 217 exponents above 220
and divisors up to 3×29, the average value of v/log2e wascalculated for optimal chains and found to be about 2%
under 45 . Indeed, successive ranges show this ratio has atendency to decrease as the exponents increase. This isbecause once the reduction operations have made theexponent less than the maximum divisor, the choice of
future divisors becomes progressively more limited. This
leads to a poorer ratio for smaller exponents.
The coefficient for the whole of a large exponent isessentially of the form (∑i vi ) (∑i log 2 mi ) where eachvi /log2mi comes from an optimal subchain. From theabove experimental data each such term should be lessthan 45 on average since each mi will fall within therange investigated (unless k is large). So the expected
coefficient ought to be better than 45 . Although a given
sub-chain being optimal at one step undoubtedly confers
some properties on the next exponent, the division
process seems to minimise most of these effects. So therelevant multiplicative properties of successive exponents
appear to be mostly independent for successive steps ofthe algorithm. Thus, individual steps resulting in a large
number of multiplications should not affect the rest of thealgorithm. Indeed, for larger exponents which require agreater number of steps, the variance in the distributionof the number of multiplications should be smaller. Thus
extreme cases will be evened out and this virtually
guarantees finding a sequence for e with under 54 log 2 emultiplications.

IEEE TRANSACTIONS ON COMPUTERS, Vol. 47, No. 7, July 1997, pp. 757-765.
which Sk £ 2620. So a complete chain is obtained with
an effort of roughly 2620×38 divisions of 256-bit
numbers by numbers averaging 6 or 7 bits.

To test 512-bit exponents and investigate how smalla coefficient c might be possible, 212 of the mostfrequently used divisor/residue pairs were chosen from
optimal chains of exponents up to over 220. They used
225 different divisors up to about 3×210, as listed in the
Appendix. Choosing k = 4, an average of under 635
multiplications was achieved for randomly generated
512-bit exponents, with a standard deviation of under 4(see Table 2). This clearly improves on the coefficient5. Assuming the distribution is close to normal, a total
4

As this is the equivalent of about 29 multiplicationsof pairs of 512-bit numbers, the effort is only justified in
RSA cryptography if the key is re-used a number oftimes. However, for k = 1 or 2 the effort is reduced toaround 5 or 14 such multiplications respectively. This
leads to nearly the same computational effort as the A, A3
method if searching for a multiplication scheme must be
done for every iteration.

of at most 640 multiplications or squarings can beeffectively guaranteed for over 91% of all 512-bit
exponents even before any of the enhancements describedin Section 7. In comparison, the square & multiply andA, A3 methods average about 23 × 511 = 766.5 and 2+
4
× 510
3

k

1

2

3

4

5

Av. No. Mults. 656.33 649.46 644.64 641.44 639.28= 682 multiplications with standard deviations

Stand. Dev.

4.0

3.6

3.6

4.0

4.7

of about 11.3 and 6.2 respectively. Interestingly, the best
choice of radix m for 512-bit exponents in the standard
m-ary method requires space for 32 powers to be held andstill only achieves a worse average of 635 multiplications(Table 2 of [15]).

Av. No. Divs.

139

157

160

164

169

1

k

2

3

Av. No. Mults. 650.36

641.7

Stand. Dev.

3.3

3.45

Table 4. Statistics for the 29 listed divisors.

Only a few divisors are actually required in order toachieve an average of under 45 . Twenty-nine divisors
with low Hamming weight were considered: the twelve
above together with 9, 41, 43, 81, 83, 161, 163, 193, 321,
323, 385, 641, 643, 769, 1281, 1283 and 1537. Here
numbers of the forms 5×2n+1 and 5×2n+3 have aminimal addition chain of length n+3 by using thesequence 1, 2, 3, 5, 10, 20, ..., 5×2n, 5×2n+1 or 5×2n+3.
About 29 of the most frequently occurring divisor/residue
pairs were selected. From Table 4, the choice k = 5 canbe seen to yield a constant of under 45 without having to4

637.33 634.67
3.4
12

Table 2. Multiplication numbers for 2

3.9(m,r) pairs.

For comparison, Table 3 contains results for the 12
divisors used in Section 5, namely 2, 3, 5, 17, 33, 49, 65,
97, 129, 257, 513 and 1025, for which the 165 mostfrequently occurring divisor/residue pairs were selected
and the difference test applied.
2

3

Av. No. Mults. 668.6

659.7

654.7

651.2 648.55

Stand. Dev.

4.75

4.35

3.8

3.54

3.6

Diff. Test c

1.30

1.29

1.278

1.27

1.265

Av. No. Divs.

168

181

184

188

192

k

1

4

use as many divisors as in Table 2. The work involved in
obtaining a division chain here can be estimated from thetable using S £ 7.12: about 7 and 27 512-bit
multiplications’ worth of effort for k = 1 and 2
respectively. Overall there is again little difference with
the A, A3 method unless the same RSA key e is re-used.

5

9

Asymptotic Results

Comparing the contents of Tables 2-4, it is clear that
adding further divisors brings diminishing returns.
Neglecting the case k = 1, the figures here (and further
results not reproduced) all suggest an exponentiallydeclining formula for the average number ofmultiplications that any choice of divisor/residue set will
require.
Good least squares approximations wereobtained as follows:

Table 3. Multiplication numbers for the 12 listed divisors.

The cost of this particular scheme for k = 5 and arandom exponent of N = 512 bits is easy to estimate. The
sum in Section 7 for approximating S shows there areabout 4.86 ways of extending any subchain by one more
pair, and, if 192 divisors are used on average, the typicaldivisor has N/192 bits. So L £ 2.7. Averaging over a few
sub-chains would also give L. If no look-up table is used,
there are around 192/k £ 38 iterations of the algorithm inTable 3: #mults £ 641.15 + 27.33 ×( 1.3656 )k–1
Table 4: #mults £ 634.68 + 21.61 ×( 1.4687 )k–1
Table 2: #mults £ 631.51 + 18.74 ×( 1.8074 )k–1

8

Colin D. Walter
Since the best divisors have already been included, thesegive a good indication that the method of division chainsas described here is unlikely to yield an average of, say,630 or fewer multiplications, although this could be
improved a little by some of the enhancements suggested
in Section 7.
The work involved in discovering an exponentiation
scheme using these methods is related to the S of Section
8. For Tables 3, 4 and 2 respectively, S is roughly 4.86,
7.12 and 19.50. Consequently, it is possible to estimate
the best divisor set to minimise the work necessary toachieve a given average number of multiplications. The
12 divisor case is best if the number of multiplications
need not be under 664.5. Then a sub-chain length k of 1should be chosen. The 29 divisor case becomes more
worthwhile once fewer than 664.5 multiplications are
required, again taking a subchain length of 1. However,the 225 divisor case would be chosen to reduce thenumber of multiplications below 650.5, choosing a
subchain length according to Table 2. Below 634, say, avariety of independent full length chains can begenerated from a number of different initial choices.
These should contain a spread of multiplication numbers,some perhaps being sufficiently below average for use.

10 Conclusion
A straightforward, cheap algorithm has been
described for reducing the number of multiplications
normally performed in an exponentiation. Almost no
extra memory is required: just one register more than forsquare and multiply. The average improvement is better
than established methods with the same memory
requirements, and exceeds that for any instance of thestandard m-ary method for 512 bit exponents.
Furthermore, the method is adaptable to a wide range ofspace and time resources, providing a variable search
space from which better evaluation orders can be found.
For any exponent e and very high probability of success,it has been shown how to find a scheme requiring under5
log 2 e multiplications or squarings when only three4
registers are available.

Acknowledgement: The author would like to thank theanonymous referee who corrected the argument in
Section 5 and kindly improved the presentation by
supplying a table similar to that given.

Exponentiation Using Division Chains

Bibliography[1]

I.E. Bocharova & B.D. Kudryashov, “Fast Exponentiation
In Cryptography”, Lecture Notes in Computer Science,Vol. 948, Springer-Verlag, 1995, pp. 146-157.

[2]

J. Bos & M. Coster, “Addition Chain Heuristics”, Proc.
Crypto ’89, Lecture Notes in Computer Science, Vol.
435, Springer-Verlag, 1990, pp. 400-407.

[3]E.F. Brickell, D.M. Gordon, K.S. McCurley and D.B.
Wilson, “Fast Exponentiation with Precomputation”,
Proc. Eurocrypt ’92, Lecture Notes in Computer Science,Vol. 658, Springer-Verlag, 1993, pp. 200-207.

[4]

C.Y. Chen, C.C. Chang & W.P. Yang, “Hybrid Method
for Modular Exponentiation with Precomputation”,
Electronics Letters, Vol. 32, No. 6, 1996, pp. 540-541.

[5]

L.S. Danilchenko, “Efficient Algorithms for Remainder
Computation and Exponentiation of Long Numbers”,
Cybernetics and Systems Analysis, Vol. 32, No. 3, 1996,
pp. 437- 441.

[6]

V. Dimitrov & T. Cooklev, “Two algorithms for Modular
Exponentiation Using Non-standard Arithmetics”, IEICE
Transactions on Fundamentals of Electronics,
Communications and Computer Sciences, Vol. E78-A,
No. 1, Jan 1995, pp. 82-87.

[7]

V.S. Dimitrov, G.A. Jullien and W.C. Miller, “Theory
and Applications for a Double-Base Number System”,
Proc. 13th IEEE Symposium on Computer Arithmetic,
6th-9th July 1997, Monterey, CA, IEEE, 1997, pp. 44-51.

[8]

P. Downey, B. Leong & R. Sethi, “Computing Sequenceswith Addition Chains”, SIAM J. Comput. vol. 10 No. 3,
1981, pp. 638-646.

[9]

Ö. Eg̃eciog̃lu & Ç. K. Koç, “Fast Modular
Exponentiation”, Comm., Control & Signal Processing,ed. E. Arikan, Elsevier Science, 1990, pp. 188-194.

[10] Ö. Eg̃eciog̃lu & Ç. K. Koç, “Exponentiation Using
Canonical Recoding”, Theoretical Computer Science,
1994, Vol. 129, No.2, pp. 407-417.
[11] P. Erdös, “Remarks on Number Theory III: On Addition
Chains”, Acta Arithmetica, Vol. 6, 1960, pp. 77-81.
[12] L. C. K. Hui & K.-Y. Lam, “Fast Square-and-Multiply
Exponentiation for RSA”, Electronics Letters, vol. 30, no.
17, Aug 1994, pp. 1396-1397.
[13] S. Kawamura, K. Takabayashi & A. Shimbo, “A Fast
Modular Exponentiation Algorithm”, IEICE Trans.
Comms, Electronics, Information and Systems, Vol E-74,
No. 8, Aug 1991, pp. 2136-2142.
[14] D. E. Knuth, “The Art of Computer Programming”, vol.
2, “Seminumerical Algorithms”, §4.6.3, 2nd Edition,
Addison-Wesley, 1981, pp. 441-466.
[15] Ç. K. Koç, “High Radix and Bit Recoding Techniques for
Modular Exponentiation”, International J. of Computer
Mathematics, vol. 40, 1991, No. 3-4, pp. 139-156.
[16] Ç.K. Koç, “Analysis of Sliding Window Techniques for
Exponentiation”, Computers & Mathematics with
Applications, Vol. 30, No.10, 1995, pp. 17-24.

IEEE TRANSACTIONS ON COMPUTERS, Vol. 47, No. 7, July 1997, pp. 757-765.
[17] D.C. Lou & C.C. Chang, “Fast Exponentiation Method
Obtained By Folding The Exponent In Half”, Electronics
Letters, Vol. 32, No.11, 1996, pp. 984-985.

multiplications on average for a 512-bit exponent.
Applying the difference test instead with 1.3 as theapproximation to c yields a much improved 668.63, as in
the table.

[18] J. Olivos, “On Vectorial Addition Chains”, Journal of
Algorithms, Vol. 2, No. 1, 1981, pp. 13-21.

Divisor 2+0: 0 ;
Divisor 3+0: 0 ;
Divisor 5+0: 0 ;
Divisor 17+0: 0 ;
Divisor 33+0: 0 ;

[19] P. de Rooij, “Efficient Exponentiation using
Precomputation and Vector Addition Chains”, Proc
Eurocrypt ’94, Lecture Notes in Computer Science, vol.
950, Springer-Verlag, 1995, pp. 389-399.
[20] Y. Tsuruoka & K. Koyama, “Fast Exponentiation
Algorithms
basedon
Batch-Processing
andPrecomputation”, IEICE Transactions on Fundamentals
of Electronics, Communications and Computer Sciences,Vol. E80-A, No.1, Jan. 1997, pp. 34-39.
[21] C. D. Walter, “Space/Time Trade-offs for Higher Radix
Modular Multiplication using Repeated Addition”, IEEE
Trans. Comp., Vol. 46 No. 2, 1997, pp. 139-141.

Min chain length = 1+1: 1
Min chain length = 2+1: 1 2Min chain length = 3+1: 1 2 3 4Min chain length = 5+1: 1 2 4 8 9 16 ;Min chain length = 6+1: 1 2 4 8 16 17 32 ;+2: 11 13+2: 19 25

Divisor 49
Min chain length = 7+0: 0+1: 2 3 4 6 8 12 16 17 24 25 32 33 48+2: 23
Divisor 65
Min chain length = 7+0: 0 ;
+1: 2 4 8 16 32 33 64+2: 24 37 49 56
Divisor 97
Min chain length = 8+0: 0+1: 2 3 4 6 8 12 16 24 32 33 48 49 64 65 96+2: 23 41 53 55 69
Divisor 129
Min chain length = 8+0: 0 ;
+1: 2 4 8 16 32 64 65 128+2: 67 73 81 96 97 192
Divisor 257
Min chain length = 9+0: 0 ;
+1: 2 4 8 16 32 64 128 129 256+2: 12 18 20 40 48 66 72 96 131 133 136 137
144 145 160 161 192 193+3: 139 147 149
Divisor 513
Min chain length = 10
+0: 0
+1: 2 4 8 16 32 64 128 256 257 512+2: 34 66 72 259 261 265 273 289 385+3: 269 277 281 293[22] Y. Yacobi, “Exponentiating Faster with Addition
Chains”, Advances in Cryptology - Eurocrypt 90, Lecture
Notes in Computer Science, vol. 473, Springer-Verlag,
1991, pp. 222-229.
[23] A. Yao, “On the Evaluation of Powers”, SIAM J.
Computing, vol. 5, 1976, pp. 100-103.
[24] S. M. Yen,
“Improved Common-Multiplicand
Multiplication and Fast Exponentiation by Exponent
Decomposition”, IEICE Trans. on Fundamentals ofElectronics Communications and Computer Sciences,Vol. E80-A, No.6, 1997, pp. 1160-1163.
[25] S.M. Yen & C.-S. Laih, “Common-Multiplicand
Multiplication and its Applications to Public Key
Cryptography”, Electronics Letters, Vol. 29, No.17, 1993,
pp. 1583-1584.
[26] C. N. Zhang, H. L. Martin, & D. Y. Y. Yun, “Parallel
Algorithms and Systolic Array Designs for RSA
Cryptosystem”, Proc. of the International Conference on
Systolic Arrays, K. Bromley, S.Y. Kung & E.
Swartzlander (Eds.), Computer Society Press, San Diego,
California, May 25-27, 1988, pp. 341-350.

Appendix

Divisor 1025
Min chain length = 11
+0: 0
+1: 1 2 4 8 16 32 64 128 256 512 513 1024
+2: 12 24 36 48 515 517 521 529 544 545 576
577 769 1152
+3: 523 531 547 549 561 579 581 585

This appendix lists the sets of divisors m used in theexamples of Sections 5 and 8, in some cases with i) the
length of a minimal addition chain used for
exponentiating to the power m (subject to the space
restrictions), and ii) lists of acceptable residues r
preceded by the number of extra multiplications required
to include the rth power in the partial result. The residue
listed is the best one to use, and so is not always the least
non-negative one.
First are the 12 divisors used in Table 3 and Section
5. Using the ratio test to order the different cases (andthe +3 lists not used in Section 5) achieves 671.66

Table 4, which achieved under 45 log 2 e multiplications,used the following 29 divisors in 519 pairs:
2
3
5
9 17 33 41 43 49 65
81 83 97 129 161 163 193 257 321 323
385 513 641 643 769 1025 1281 1283 1537

10

Colin D. Walter
The following 225 divisors were used in constructing
Table 2:2
3
5
7
9 11 13 14 17 19
23 26 31 33 37 41 43 49 50 59
61 65 66 67 73 77 83 97 98 107
109 113 121 129 131 133 137 145 149 163
168 193 194 197 199 211 227 229 233 257
259 261 265 281 289 290 292 293 296 298
321 371 373 385 386 389 397 409 437 449
481 483 513 515 517 521 529 540 545 546
552 553 560 561 562 577 578 580 581 584
586 588 608 641 642 643 644 646 648 656
673 677 683 730 739 769 770 773 775 780
781 785 793 801 803 809 812 813 846 848
852 869 884 888 897 899 900 901 902 904
905 906 912 914 920 928 944 961 964 966
968 972 976 984 1008 1025 1027 1029 1033 1041
1057 1072 1073 1097 1153 1154 1156 1157 1170 1200
1264 1281 1282 1283 1348 1352 1361 1364 1384 1396
1412 1415 1422 1424 1432 1444 1448 1456 1464 1468
1470 1477 1482 1488 1490 1500 1504 1520 1524 1528
1537 1538 1540 1541 1543 1544 1546 1550 1553 1556
1561 1564 1569 1572 1576 1585 1588 1592 1593 1601
1603 1604 1606 1608 1609 1610 1612 1613 1616 1617
1618 1624 1626 1648 1671

Exponentiation Using Division Chains