Fast Scalar Multiplication for ECC over GF(p)Using Division Chains
Colin D. Walter
Information Security Group, Royal Holloway, University of London,
Egham, Surrey, TW20 0EX, United Kingdom
Colin.Walter@rhul.ac.uk

Abstract. There have been many recent developments in formulae foreﬃcient composite elliptic curve operations of the form dP +Q for a small
integer d and points P and Q where the underlying ﬁeld is a prime ﬁeld.
To make best use of these in a scalar multiplication kP , it is necessary
to generate an eﬃcient “division chain” for the scalar where divisions ofk are by the values of d available through composite operations.
An algorithm-generating algorithm for this is presented that takes
into account the diﬀerent costs of using various representations for curve
points. This extends the applicability of methods presented by Longa& Gebotys at PKC 2009 to using speciﬁc characteristics of the target
device. It also enables the transfer of some scalar recoding computation
details to design time. An improved cost function also provides better
evaluation of alternatives in the relevant addition chain.
One result of these more general and improved methods includes aslight increase over the scalar multiplication speeds reported at PKC.
Furthermore, by the straightforward removal of rules for unusual cases,
some particularly concise yet eﬃcient presentations can be given foralgorithms in the target device.
Keywords: Scalar multiplication, multibase representation, addition chain,division chain, exponentiation, DBNS, elliptic curve cryptography.

1

Introduction

Exponentiation has been the subject of much study over the years. Classic summaries of the state of the art a decade ago are Knuth and Gordon [14,12]. The
aim of these methods is almost always minimising the operation count in an addition chain for the exponent. This is an NP-hard problem. More recent work hasconcentrated on wider optimisation issues, such as those imposed by constraints
on the size of any pre-computed table of values, the working space required,
the relative costs of squaring and multiplication, and reducing the side channel
leakage generated by a particular operation sequence.


The work described in this paper has been supported [in part] by the European
Commission through the ICT programme under contract ICT-2007-216676 ECRYPTII.

Y. Chung and M. Yung (Eds.): WISA 2010, LNCS 6513, pp. 61–75, 2010.
c Springer-Verlag Berlin Heidelberg 2010

62

C.D. Walter

Recently, eﬃcient composite elliptic curve operations of the form dP +Q havebeen developed for double-and-add (d = 2), triple-and-add (d = 3) and quintupling (d = 5) etc. over Fp [6,18,16]. By careful choice of curve formulae and aﬃnecoordinates for the table entries, these become faster than equivalent combinations of the normal point addition and doubling formulae, and so should lead tofaster scalar multiplication. With each application of the operation dP +Q to theaccumulating scalar multiple, the part of the key which has yet to be processedis reduced by a factor d using integer division. This is exactly the process forgenerating a division chain, and so division chains are the natural setting for
optimising scalar multiplication when these composite operations are available.
Division chains were ﬁrst deﬁned in [19] and [20] where the author uses themto derive eﬃcient addition chains in the context of limited table and working
space. The resulting chains are signiﬁcantly shorter than those for binary exponentiation when using just one more register. Without cheap inversion (as is thecase for RSA), exponentiation by k can be performed in a right to left direction
using an average of about 1.4 log2 (k) operations using divisors d = 2 and d = 3(op. cit. §5). This improves on the usual 32 log2 (k) operations using only the divisor d = 2 but is less eﬃcient than the 43 log2 (k) for sliding 2-bit windows1 which
also uses space equivalent to one more register than binary methods. Althoughwith more divisors the speed can be improved to under 54 log2 (k) operations using the same space, this requires much more pre-processing of the exponent andapplies only to the right to left direction for which the composite elliptic curve
operations dP +Q do not have more eﬃcient forms. Figures are not provided in[19] for the addition-subtraction chains relevant for elliptic curves, nor are they
provided for the left-to-right direction which is marginally less eﬃcient for spaceconstrained implementations but which is of relevance here where the composite
operations provide increased eﬃciency. One aim of the present work is to ﬁll thisgap.
For application to elliptic curves, it is the number of ﬁeld operations which
particularly determines the time eﬃciency of scalar point multiplication. This
number depends to a large extent on the form of the curve equation and thecoordinate representations chosen for the table entries and intermediate points
during processing2 [13]. The resulting cost of point doublings, additions, triplings
etc., rather than the length of the addition sub-chain, is the natural input tothe cost evaluation function from [19] which is used to compare competing alternative divisors. This explicit function enables the design time derivation of thealgorithm for choosing divisors at run time in a more reﬁned manner than in[17]. Moreover, it transfers all of the search for the best divisor to design time.
1

k is partitioned from right to left into “digit” subsequences 0, 01 and 11. Exponentiation is performed left to right using table entries for 01 and 11. Even digits 0 occur
with probability 12 whereas the odd digits 01 and 11 occur with probability 14 . So
an exponent bit has probability 13 of being associated with digit 0, and 23 with digit
01 or 11. Thus the average cost per exponent bit is asymptotically 13 · 11 + 23 · 32 = 43
operations.
2
Some speed-ups require additional space for four or more coordinates per point.

Fast Scalar Multiplication for ECC over GF(p) Using Division Chains

63

The less structured double base and multibase representations of Dimitrov
et al. [4,5,6] represent the scalar k as a signed sum of powers of two (or more)numbers. The extra freedom there allows the use of a slower, greedy algorithm
to obtain a very compact representation. However, compactness is not the aim
here and extra conditions are required to provide the structure required here
for an eﬃcient conversion to a fast addition sequence. These extra conditions
are those of a randomary representation [22]. Similar representations have beencreated more recently by Ciet, Longa et al. [3,15,17] but restricted to cases inwhich all digits are zero except those associated with divisor 2 (or its powers).
Here those restrictions are lifted to give a very general setting which encompasses all this previous work on the more structured multibase forms. A side
eﬀect of the generality is a modest increase in the speed records achieved in [17].
But the results are also applicable in wider contexts than the particular devices
targetted there. For example, the relatively cheap cost of tripling on DIK3 curves[10] makes use of the new availability of non-zero digits associated with base 3.
The main body of the paper describes the above in more detail: division
chains and their notation are covered initially; then a discussion on selecting
divisor and digit sets followed by details of the cost function for determining therecoding process, and results to evaluate a variety of parameter choices. Finally,
Appendix B contains an explicit and very straightforward example algorithm forgenerating an eﬃcient division chain in the case of divisor set {2,3}.

2

Addition Chains

The concept of addition chains and addition-subtraction chains needs extending
to the context here. Let OP be the set of “distinguished” operations which wewish to use in exponentiations schemes. For convenience only unary or binary
operations are considered. For the (additive) group G of interest, they combine
one or two multiples of a given g ∈ G into a higher multiple of g. So, for eachp ∈ OP there will be a pair (λ, μ) ∈ Z2 such that p(g, h) = λg+μh for all
g, h ∈ G. As special cases, (1, 1) yields g+h (an addition), (2, 0) yields 2g (adoubling) and (−1, 0) provides −g (an inversion). Thus OP can be viewed as asubset of Z2 . Recent research has provided double-and-add, triple-and-add andeven quintuple operations [18,16] to help populate OP. These are of the typep = (2, d) or p = (3, d) etc. where the d-th multiple of the initial point g is in the
pre-computed table. These operations lead to a more general addition chain:
Definition 1. For a set OP of linear binary operators on an additive group, an
OP-addition chain for k ∈ Z of length n is a sequence of quadruples (ai , bi , ki , pi )∈ Z3 ×OP, 0 ≤ i ≤ n, such that, for all i > 0,
– ki = λi ai +μi bi for pi = (λi , μi ) ∈ OP
– ai = ksi and bi = kti for some si and ti with 0 ≤ si < i and 0 ≤ ti < i– (a0 , b0 , k0 , p0 ) = (1, 0, 1, (1, 0))– kn = k.

64

3

C.D. Walter

Representations from Division Chains

Addition chains are too general for use in devices with limited memory. The
standard left-to-right binary exponentiation algorithm overwrites all previous
calculations with the latest value and so its operations can only access the mostrecent value and the initial input. Hence its addition chain can only contain
triples of the form (ki−1 , ki−1 , 2ki−1 ) (a squaring) or (ki−1 , 1, ki−1 +1) (a multiplication). Schemes with pre-computed tables are similarly restricted.
The natural extension of this restriction allows only quadruples of the form
(ki−1 , 1, λi ki−1 +μi 1, (λi , μi )) in an OP-addition chain, i.e. ki = λi ki−1 +μi . Consequently, just as the sequence of squares and multiplications completely determines the exponent k in binary exponentiation, so the sequence of operations(λi , μi ) also determines k. Given k, we can obtain such a sequence (in reverse
order) by iteratively choosing (λi , μi ) ∈ OP such that μi lies in the residue classof ki mod λi , and then performing the division ki−1 = (ki −μi )/λi . The process
starts with k = kn and ﬁnishes with k0 = 0. This deﬁnes a division chain for
an integer k, and it is speciﬁed by the list of pairs (λ, μ) corresponding to theconstituent operations. This is clear from the presentation of k as
k = (((μ1 λ2 + μ2 )λ3 + ... + μn−2 )λn−1 + μn−1 )λn + μn

(1)

Subscripts for addition chains are used in the opposite direction from those indigit representations of numbers. In order to obtain the customary notation forthe latter, let us re-number the pairs (λ, μ) by deﬁning (ri , di ) = (λn−i , μn−i ).
Then the division chain determines a randomary representation3 of k [19]:k = (((dn−1 rn−2 + dn−2 )rn−3 + ... + d2 )r1 + d1 )r0 + d0

(2)

in which the ri are called bases and the di are called digits. If all the divisors ri
were equal, then this would simply be a representation in radix r = ri . Digits
would be chosen in [0..r − 1] to yield the standard representation. Here the bases
will belong to some small set B ⊂ N such as {2,3,5} and the digits will lie in a
subset D ⊂ Z which may contain negative values and may depend on the base.
An example from [22] is23510 = (((((1)3 + 0)2 + 1)5 + 4)2 + 0)3 + 1.
Following standard subscript notation for specifying the base, this is written
23510 = 12 03 12 45 02 13 . Digit/base pairs are obtained using the usual change ofbase algorithm except that the base may be varied at each step: for i = 0, 1, ...
extract digit di by reducing modulo ri , and divide by ri . In the example of 23510 ,3

Developed from the words binary, ternary, quaternary etc. to indicate the apparent
random nature of base choices as used in, for example, [21]. This is a not a multi-base(or mixed base) representation in the original sense of that term, and so “randomary”
may 
be useful to distinguish the two forms. A double base (DBNS) example is writingk = i ±2bi 3ti , which has a diﬀerent structure.

Fast Scalar Multiplication for ECC over GF(p) Using Division Chains65

the ﬁrst base is 3 and 235 mod 3 = 1 is the ﬁrst digit (that of index 0). The
algorithm continues with input (235 − 1)/3 = 78 to determine the next digit.
Such representations determine a left-to-right exponentiation scheme using
Horner’s rule (which follows the bracketing above from the inside outwards). In
an additive group G the iterative step is
ki g = rn−i (ki−1 g)+(dn−i g)(3)using table entry dn−i g and starting at the most signiﬁcant, i.e. left, end.

4

Choosing the Base Set B

This section considers the choice of base set when the scalar k is used just once
and the full cost of creating the multibase representation must be counted. If acryptographic token were to use the same key throughout its lifetime then more
eﬀort could be spent on optimising the representation oﬄine before initialisation,and so diﬀerent rules would apply. Such a situation is not covered here.
Application of base r reduces the key by a factor of almost exactly r andso shortens k by close to log2 r bits. If this requires cr group operations whenused, then its approximate cost for comparing with other bases is cr / log2 r(when neglecting other beneﬁts). Apart from powers of 2 this ratio is lowest fornumbers of the form r = 2n ±1 as these consume very close to n scalar bits
for a cost of only n +1 or n +2 operations according to whether the digit iszero or represented in the pre-computed table. For large n this is going to bebetter than the average cn for the whole algorithm on a key of n bits, wherec is typically in the region of 1.2 to 1.25. With a table of T entries such bases
can be used in about 2T +1 out of r cases. k mod r must be computed every
time to check the usability of r and, for k in binary, the eﬀort involved in this is
proportional to n log(k). However, the saving to the scalar multiplication is less
than (c−1)n group operations per use, giving an average of (2T +1)(c−1)n 2−n
per calculation of k mod r . Hence, as n increases the re-coding cost will quickly
outweigh any beneﬁt obtained from its use. The same reasoning applies to anylarge non-2-power radix. Hence B should only contain small numbers.
Generation of addition-subtraction chains for all small integers shows that,with few exceptions, exponentiation by most bases r requires close to 1.25 log2 r
operations4 without counting the extra one for any point addition. With an aimof achieving an average of fewer than 1.25 operations per bit, such base choices
are not going to be of much use unless the associated digit is zero or there isforeseeable beneﬁt, such as a zero digit in the immediate future.
An investigation of small bases was carried out. It was noted ﬁrst that no baser below 29 required more than two registers for executing a shortest addition
4

The distribution is much more uniform than for NAF, and the 1.25 operations per
bit slowly decreases. NAF gives a chain with an asymptotic average of 43 log2 r operations. 27 is the ﬁrst number with an add/sub chain shorter than the NAF, and 427
is the ﬁrst number with an add/sub chain two operations shorter than the NAF. 31
is the ﬁrst that needs a subtraction for its shortest chain.

66

C.D. Walter

chain for r – the minimum possible for numbers other than powers of 2. So space
is not an issue. For the most part, only prime candidates need be considered forthe following reason. Suppose radix r = st is composite and it is to be used on
k with digit d. Then k would be reduced to (k−d)/st. However, application of
the base/digit pairs (s, 0), (t, d) has the same eﬀect. Consequently, the value of
including r as a divisor depends on whether or not a minimal addition chain forr is shorter than the sum of the chain lengths for any factorisation of r. As an
example, no powers of 2 above 2 itself need be included. Thus, the pair (4, 3) canbe replaced without cost by (2, 0) followed by (2, 3). However, 33 should be kept
in as it has a chain of length 6 which is less than the combined chain lengths for
3 and 11. The investigation revealed {33, 49, 63, 65, 77, 95, 121, 129, 133, 143, . . .}to be the sequence of composite numbers having shorter addition/subtraction
chains (by sequence length) than the sum of those of any factorisation.
In fact, shorter needs to be interpreted at a ﬁner level of detail. Powers of
3 then become more interesting: two applications of base 3 requires two point
doublings and two point additions but one application of base 9 can be achieved
with three doublings and one addition. This is the same number of point operations, but should be cheaper since a doubling ought to have a lower cost thana point addition. Similarly, 15 is cheaper than 3 and 5 separately.5
After omissions for the above reasons and because of poor ratios cr / log2 r,
there remains only the initial part of the sequence{2, 3, 5, 9, 15, 17, 31, 33, 47, 63, 65, 77, 97, 127, 129, . . .}.
In particular, 7, 11 and 13 have been removed. They have among the very highest
costs of any base, with cr / log2 r well above 43 . This high cost holds even when
the relative cost of a doubling is accounted for and has a value anywhere between
0.5 and 1.0 times the cost of an addition. As the ratio approaches 0.5 the popularbase 3 also becomes extremely expensive under this measure, whereas bases 5
and 9 are consistently good. 3 is better than 5 only once the ratio is above 0.87.
An exhaustive search was made for shortest divisor chains using choices fromthe above sequence to determine which bases were the most useful since any
resource-constrained device can only allow a few possibilities. Shortest was measured assuming point doubling to have a cost between half and one times thecost of an addition (see [17], Table 1). The digit set was allowed to be any subsetof {−7, ..., +7}. In all cases, the frequency of use was close to, in descendingorder, {2; 3; 5; 9, 7; 17, 13, 11; 33, 31, 23, 19; 65, 63, ...; 129, 127, ...} which is 2 concatenated with the sequences 2n+1 +1, ..., 2n +3 for n = 0, 1, 2, ...6 As expected,
the frequencies of bases r decrease as r increases. Consequently, it is reasonable
to consider only bases in the set {2, 3, 5, 9, 17}7 since, with such low frequencies,
more bases can only give negligible speed-up to a scalar multiplication. Furthermore, extra bases would require more code and more pre-computation. In fact,
5

The ﬁner detail of gains from using composite operations is treated below.
This listing assumes that bases r of interest have minimal addition chains of length1+log2 (r−2), which is true for the range of interest and the above restrictions.
7
In order of decreasing frequency when doubling is half the cost of addition.
6

Fast Scalar Multiplication for ECC over GF(p) Using Division Chains67

the decrease in frequencies is much greater when the cost of doubling approaches
that of point additions. In this case, it is hardly worth considering more thanjust {2, 3, 5, 17}.8 However, when doubling is just half the cost of adding, notethat 7, 11 and 13 are still possible options. Despite the heavy average cost per
bit by which these bases reduce k, they can lead to future values of k whichhave cheaper than average reductions and enable the extra cost to be recovered.
The search for optimal chains also revealed that odd bases rarely make use of
non-zero digits. This was the case even when the addition of a table element was
assigned just half the cost of the additions used in multiplying by a base. Thus
almost no extra speed is likely to be achieved if non-zero digits are allowed for
odd bases when composite curve operations are introduced.
This leads to the most sensible base choices for embedded cryptosystems being
B = {2, 3} or B = {2, 3, 5}. They were ﬁrst considered in [21,22] and more
recently by [18,17] where non-zero digits are only allowed for base 2. Larger sets
provide virtually no reasonable advantages and might only be considered out ofvanity for attaining a rather marginal increase in speed.

5Choosing the Digit Set D

Using the very rough cost function given by addition chain length in the previous
section it became clear that the fastest scalar multiplication methods would use
non-zero digits almost exclusively with the base 2. Consequently, digits which
are multiples of 2 can be eliminated as follows. Applying the base/digit pair (r, 0)
followed by (2, 2d) is equivalent to applying the (r, d) followed by (2, 0). So thetable need not contain even digit multiples of the initial point to cover this case.
Similarly, (2, d ) followed by (2, 2d) is equivalent to (2, d+d ) followed by (2, 0),
but the new expression saves one point addition. If d+d is not an acceptable
digit, (2, d+d ) could be implemented as a double and add (of multiple d) followedby an extra addition (of multiple d ), thereby again eliminating the need for any
even digit multiples in the table or any further adjustment to the representation.
Alternatively, d+d could be split into an acceptable digit and a carry up to thepreceding base/digit pair as in the constant base case, with this process repeateduntil the digit problem is resolved at no additional cost.
As in the case of NAF, and conﬁrmed by the search in §4, the choice of anon-zero digit is almost always going to be made to make the next one or moredigits zero. By taking the digit set D = {0, ±1, ±2, . . . , ±(2T −1)} for given tablesize T , it is possible to maximise the power of 2 achievable as a factor of thenext value of k during recoding. Because the composite operators work equallywell for any table entry, the desired operation set now has the form OP = B×D
where the two factors are now known.

6

Optimising the Representation

As just noted, a primary aim in the multibase algorithm is to select bases forwhich the digit is zero as this saves the point addition. When this is not possible,
8

In order of decreasing frequency when doubling has the same cost as addition.

68

C.D. Walter

the recoding strategy should make a base/digit choice to ensure it happens nexttime or as soon as possible thereafter. This requires knowing the key modulo a
power of each available base. Let π be a common multiple of the elements in thebase set B which contains only primes dividing available bases and in which it is
helpful to have a high power of two because of the frequency with which 2 will
be chosen as a base. Decisions will be based on ki mod π. Indeed, for arbitrary p
prime to all base choices, the value of ki mod p can clearly have no inﬂuence on
the choice of base. However, if π were picked large enough (e.g. larger than k) wewould have complete knowledge of k and so there would be enough information
to determine an optimal division chain.
Divisor choices are not independent. The choice of the ﬁrst divisor is partly
aﬀected by the best choice for the next, and less so by subsequent ones. After achain of four or more, simulations show that the eﬀect is minimal. Hence we wantto pick the ﬁrst divisor to minimise the average cost of the addition sub-chain
formed by the next four or so divisors. So let λ be the length of sub-chain to
be investigated. This would be the window size if the base were ﬁxed. Initially,let π be the least common multiple of the λth powers of each base in B. If the
pre-computed table has T elements and δ is minimal such that T ≤ 2δ then abase 2 digit can sometimes be chosen to guarantee divisibility by at least 2δ+2 .
(For example, with T = 3 and table entries for {±1, ±3, ±5}, we have δ = 2
and, when ki ≡5 mod 8, either of the digits −3, 5 might be used, but one will
make ki −di ≡0 mod 16.) In such cases, in eﬀect, we have a base of 2δ+2 . Thus
the base 2 digits cannot be chosen accurately towards the end of the window
unless π is augmented by another factor of at least 2δ+1 . Moreover, this example
shows that the cost beneﬁts of digit choices continue beyond the end of thewindow.
For each residue mod π, all possible sub-chains of length λ are generated.
The ﬁrst pair (r, d) is then chosen from the “best” sub-chain and assigned tobe the choice whenever ki has that value modulo π. Explicit code for this ispresented
 in the ﬁrst algorithm of Appendix A. Since the search takes a total ofO(π{T r∈B r1 }λ ) time, λ = 3 is easily achievable with reasonable resources.
In the next section, the device-speciﬁc deﬁnition of “best” is given properly.
However, the above association of pairs (r, d) with residues mod π is obtained
from an oﬀ-line search and yields a simple algorithm for programming into thetarget device, such as the second algorithm in Appendix A. Most divisor choices
are determined by very simple rules such as selecting the pairs (2, 0) and (3, 0)
respectively when the residue mod π is divisible by 2 or 3. So the table of data
can generally be converted to a much more compact set of simple rules. The fullset of rules mod π may contain many exceptional cases. They become less criticalto the performance of the algorithm as the number of residues aﬀected decreases,
and so some very small cases can be absorbed into larger ones with very little loss
in performance, thereby reducing code size. Appendix B contains two examplesof the recoding algorithm’s iterative step showing just how compact the resulting
code can be and yet still provide a very eﬃcient scheme.

Fast Scalar Multiplication for ECC over GF(p) Using Division Chains

6.1

69

The Detailed Cost Function

In §4 the cost per bit of applying base r is given as cr / log2 r where cr was initiallytaken as the length of the addition chain used for r in order to ascertain whichwere the best values to put in B. From now on, costs cr,d are taken at the accuracy required to diﬀerentiate successfully between competing choices for (r, d).
If reading and writing operations, ﬁeld additions and scalar ﬁeld multiplications
are relatively cheap, cr,d is given by the appropriately weighted sum of ﬁeld
multiplication and squaring counts. They depend on the form of the equationand the coordinate representation chosen for the elliptic curve and its points.
The reader is referred to [17], Table 1, for examples covering the specialised
eﬃcient doubling, tripling and quintupling operations of interest here.
A sequence of λ iterative steps of the type (3) can be composed into a single
operation of the same type, represented by a base/digit pair (r, d) where the newbase r is the product of those of the constituent pairs. The cost function cr,dcan be extended to this by adding the costs of the component pairs. However,the per bit cost cr,d / log2 r is no longer a good estimate of the relative value ofthe choice of either (r, d) or the ﬁrst pair in the window. The reduction achievedby any sub-chain is typically followed by an average reduction. Suppose c is theaverage cost per bit and an alternative sub-chain has cost cr ,d to achieve a
reduction by (r , d ). Removing a further factor of r/r from (ki −d )/r costs anaverage of c log2 (r/r ), giving a total of cr ,d +c log2 (r/r ) for reducing ki by afactor of r. So the ﬁrst choice is better on average if cr,d < cr ,d +c log2 (r/r ),i.e. if cr,d −c log2 r < cr ,d −c log2 r . Thus the “best” sub-chain is taken to be
the one minimising
cr,d − c log2 r(4)This is simply the extra work required above the average expected. Although useof (4) requires advance knowledge of the average cost per bit c, various values
can be tried in practice, of course, in order to converge on the smallest.
Finally, no choice of base/digit sequence for the window makes full use of
the value ki mod π. For a sub-chain equivalent to (r, d), cr,d does not take into
account properties of (ki −d)/r mod π/r. As noted above, the ﬁnal digit in thewindow is usually chosen, like the others, to give exact divisibility in the following
digit when that possible. So the sub-chain pair (r, d) in (4) should be augmented
to include pairs (r , 0) representing any remaining divisibility in ki −d mod π.
6.2

A Markov Process

An interesting consequence of specifying the divisor/digit pair to choose given
the residue of ki modulo π is that these residues become non-uniform after theﬁrst iteration and, after about half a dozen choices, have settled close to their
asymptotic probabilities. This is a Markov process which, in the limit, leads toa precise value for the work c per bit required when evaluating (4).
The process is described by a π×π transition probability matrix P = (pij )where pij is the probability that an input of i mod π will generate an outputof j mod π. If the pair (ri , di ) is associated with i mod π then the output j

70

C.D. Walter

will be one of the ri values (i−di +tπ)/ri mod π for 0 ≤ t ≤ ri −1. As these areall equally likely, each of these pij has the value ri −1 and for other values of j,(m)
pij is zero. The matrix P m contains entries pij which are the probabilities of
input i mod π resulting in an output of j mod π after m divisors are applied
to k using the algorithm. Averaging over all i gives the probability of residue
j mod π after m iterations. This converges very quickly to a steady state inwhich the probabilities typically lie in a range between 34 π −1 and 54 π −1 .
If pj is the asymptotic probability for residue class j mod π and cj the costof the pair (rj , dj ) associated with that class, then the average cost per bit isπ−1
j=0 pj cj(5)
c = π−1
j=0 pj log2 (rj )In fact, using a slightly more aggressive (smaller) value for c in (4) seems to yield
the best value for c in (5). Our overall objective is, of course, to minimise this.
The transition matrix is usually too big for practical calculations. Instead, aMonte Carlo simulation can be performed. A large number of random keys k is
generated, written in base π for convenience, and recoded to ﬁnd the probabilities. Since (4) may have failed to capture all the useful information in the associated residue modulo π, this process can be repeated on neighbouring schemes
to identify the local minimum for c: for each residue i mod π alternative pairs(ri , di ) can be tried to improve c. The same process can also be used to evaluate
potential simpliﬁcations to the scheme.

7

Test Results and Evaluation

Table 1 shows ﬁgures comparing the above algorithm with the best earlier results
as provided by the reﬁned NAF method of Longa and Gebotys in [17], Table
4. It illustrates the speed-ups for three sets of projective coordinate systems.
These are standard Jacobian coordinates (with a = −3), inverted Edwards coordinates [7,1] and extended coordinates on the Jacobi Quartic form [2,13]. The
cost per operator columns give the number of ﬁeld multiplications or equivalents (counting 0.8 for each squaring) which are used for doubling (D), tripling(T), quintupling (Q) and the addition of a table entry (A) within a doubleadd, triple-add etc. operation in the given coordinate system. The stored points
were the ﬁrst few of 1, 3, 5, 7,... times the initial point. The number of them isrecorded in the third column group, and the cost of initialising them is includedin the operator count totals. In order to make scaling to other key lengths easier,an initialisation optimisation [17,8,16] which saves 10 to 12 ﬁeld multiplications
has not been applied to either set of results. The ﬁnal column is the asymptotic
number of ﬁeld multiplications (or equivalents) per bit of k, calculated as in §6.2.
This is the value of c in (5).
The methodology here performs a systematic search at design time to determine the best run-time choice. This choice is wider than that in [17] because,
for example, 2 is not automatically the divisor when k is even (see Appendix B
for an example), bases which divide k exactly need not be chosen as divisors,

Fast Scalar Multiplication for ECC over GF(p) Using Division Chains

71

Table 1. Comparative Point Multiplication Speeds for 160-bit Scalar (cf [17], Table 4)Mult Cost/Op
B#Stored π
Op. Count Cost/Bit
D T Q APts
Sq = 0.8M
c
JQuart 6.0 11.2
8.4 {2,3}
81229.2
7.50
”
” ”
”
”
1
1288.48.10
InvEdw 6.2 12.2
8.8”
81277.1
7.79
Jacobian 7.0 12.6
10.2”
81445.1
8.91
JQuart 6.0 11.2
8.4 {2,3}
8
213 331207.1
7.32
”
” ”
”
”
4
213 331212.0
7.52
”
” ”
”
”
2
217 341253.8
7.84
”
” ”
”
”
1
217 341287.6
8.10
InvEd 6.2 12.2
8.8”
8213 331242.5
7.55
Jacobian 7.0 12.6
10.2”
8
213 331414.9
8.57
JQuart 6.0 11.2 17.2 8.4 {2,3,5}
81226.5
7.48
Jacobian 7.0 12.6 19.6 10.2”
81435.2
8.84
JQuart 6.0 11.2 17.2 8.4 {2,3,5}
8
29 32 52 1203.7
7.30
Jacobian 7.0 12.6 19.6 10.2”
8
29 32 52 1409.8
8.53

Method Coords
[17]
”
”
”
Here
”
”
”
”
”
[17]
”
Here
”

and divisors other than 2 can have non-zero digits associated with them. All thisautomatically leads to the possibility of faster scalar multiplication. However,
the tabulated results here represent a less deep search – a larger value for π isnecessary to uncover all the conditions listed in [17], Tables 2 and 3. Further
computation could have been used to obtain them. There is no clear cost function in [17] but an explicit cost function has been given above in (4). This is not
the obvious one (which is sub-optimal) but arguments have been presented tojustify the conclusion that it is better. This diﬀerence probably also contributes
to the greater speeds here, which are asymptotically 2.5% faster than [17].
The best value for money seems to occur with tables of 4 points, being under
1% slower than their 8 point counterparts. For the fastest times, observe that
the Jacobi Quartic coordinate representation [13] can be used to perform 160bit scalar multiplication with divisors 2 and 3 and a table of 8 pre-computed
points using, on average, the equivalent of fewer than 1200 ﬁeld multiplications
when the optimisations of [8,16] are applied. This can be achieved with verystraightforward code, which is given in the ﬁrst example of Appendix B.

8Conclusion

A detailed methodology has been presented for deriving compact, device-speciﬁc
algorithms which perform very eﬃcient scalar multiplication for elliptic curves
in modestly resource constrained environments when the scalar is an unseen
fresh, random, cryptographic key. The method is based on multibase representations of the scalar and careful construction of a cost function for evaluating
alternatives. A valuable output demonstrating the eﬃcacy of the approach isthe fastest known algorithm in certain contexts, such as when composite curve

72

C.D. Walter

operations are available. Another output is a couple of examples which illustrate
the compactness of resulting code for the iterated recoding step.
With an initial recoding cost equivalent to two or three ﬁeld multiplications,the exponentiation algorithm is eﬃcient and avoids an expensive run-time searchby resolving the search questions at design time. Part of the speed-up is derived
from a much fuller choice of base/digit pairs at each step than prior algorithms
and part from a new, explicit device-speciﬁc cost function for evaluating diﬀerent
options in the search for the most eﬃcient addition chain. That function can,and should moreover, use detailed knowledge of the times for various ﬁeld operations. Exactly the same methods can be applied with more extensive preparatorycomputation to generate still faster or compact algorithms than tabulated here.

References
1. Bernstein, D., Lange, T.: Analysis and Optimization of Elliptic-Curve Single-Scalar
Multiplication. Cryptology ePrint Archive, Report 2007/455, IACR 2007 (2007)
2. Billet, O., Joye, M.: The Jacobi Model of an Elliptic Curve and Side-Channel
Analysis. In: Fossorier, M.P.C., Høholdt, T., Poli, A. (eds.) AAECC 2003. LNCS,
vol. 2643, pp. 34–42. Springer, Heidelberg (2003)
3. Ciet, M., Joye, M., Lauter, K., Montgomery, P.: Trading Inversions for Multiplications in Elliptic Curve Cryptography. Designs, Codes and Cryptography 39(2),189–206 (2006)
4. Dimitrov, V., Cooklev, T.: Two Algorithms for Modular Exponentiation using
Non-Standard Arithmetics. IEICE Transactions on Fundamentals of Electronics,
Communications and Computer Sciences E78-A(1), 82–87 (1995)
5. Dimitrov, V.S., Jullien, G.A., Miller, W.C.: Theory and Applications for a DoubleBase Number System. In: Proc. 13th IEEE Symposium on Computer Arithmetic,
Monterey, July 6-9, pp. 44–51. IEEE, Los Alamitos (1997)
6. Dimitrov, V.S., Imbert, L., Mishra, P.K.: Eﬃcient and Secure Elliptic Curve Point
Multiplication using Double-Base Chains. In: Roy, B. (ed.) ASIACRYPT 2005.
LNCS, vol. 3788, pp. 59–78. Springer, Heidelberg (2005)
7. Edwards, H.: A Normal Form for Elliptic Curves. Bull. Amer. Math. Soc. 44, 393–
422 (2007)
8. Elmegaard-Fessel, L.: Eﬃcient Scalar Multiplication and Security against Power
Analysis in Cryptosystems based on the NIST Elliptic Curves over Prime Fields,
Masters Thesis, University of Copenhagen (2006)
9. Fouque, P.-A., Valette, F.: The Doubling Attack – Why upwards is better than
downwards. In: Walter, C.D., Koç, Ç.K., Paar, C. (eds.) CHES 2003. LNCS,
vol. 2779, pp. 269–280. Springer, Heidelberg (2003)
10. Doche, C., Icart, T., Kohel, D.R.: Eﬃcient Scalar Multiplication by Isogeny Decompositions. In: Yung, M., Dodis, Y., Kiayias, A., Malkin, T.G. (eds.) PKC 2006.
LNCS, vol. 3958, pp. 191–206. Springer, Heidelberg (2006)
11. Giessmann, E.-G.: Ein schneller Algorithmus zur Punktevervielfachung, der gegen
Seitenanalattacken resistent ist. In: Workshop über Theoretische und praktische
Aspekte von Kryptographie mit Elliptischen Kurven, Berlin (2001)
12. Gordon, D.M.: A Survey of Fast Exponentiation Algorithms. Journal of Algorithms 27, 129–146 (1998)
13. Hisil, H., Wong, K., Carter, G., Dawson, E.: Faster Group Operations on Elliptic
Curves. Cryptology ePrint Archive, Report 2007/441, IACR (2007)

Fast Scalar Multiplication for ECC over GF(p) Using Division Chains73

14. Knuth, D.E.: The Art of Computer Programming, 2nd edn. Seminumerical Algorithms, vol. 2, §4.6.3, pp. 441–466. Addison-Wesley, Reading (1981)
15. Longa, P.: Accelerating the Scalar Multiplication on Elliptic Curve Cryptosystems
over Prime Fields, Masters Thesis, University of Ottawa (2007)
16. Longa, P., Miri, A.: New Composite Operations and Precomputation Scheme for
Elliptic Curve Cryptosystems over Prime Fields. In: Cramer, R. (ed.) PKC 2008.
LNCS, vol. 4939, pp. 229–247. Springer, Heidelberg (2008)
17. Longa, P., Gebotys, C.: Fast Multibase Methods and Other Several Optimizations
for Elliptic Curve Scalar Multiplication. In: Jarecki, S., Tsudik, G. (eds.) PKC
2009. LNCS, vol. 5443, pp. 443–462. Springer, Heidelberg (2009)
18. Mishra, P.K., Dimitrov, V.: Eﬃcient Quintuple Formulas for Elliptic Curves andEﬃcient Scalar Multiplication using Multibase Number Representation. In: Garay,
J.A., Lenstra, A.K., Mambo, M., Peralta, R. (eds.) ISC 2007. LNCS, vol. 4779, pp.
390–406. Springer, Heidelberg (2007)
19. Walter, C.D.: Exponentiation using Division Chains. In: Proc. 13th IEEE Symposium on Computer Arithmetic, Monterey, CA, July 6-9, pp. 92–98. IEEE, Los
Alamitos (1997)
20. Walter, C.D.: Exponentiation using Division Chains. IEEE Transactions on Computers 47(7), 757–765 (1998)
21. Walter, C.D.: MIST: An Eﬃcient, Randomized Exponentiation Algorithm for Resisting Power Analysis. In: Preneel, B. (ed.) CT-RSA 2002. LNCS, vol. 2271, pp.
53–66. Springer, Heidelberg (2002)
22. Walter, C.D.: Some Security Aspects of the MIST Randomized Exponentiation
Algorithm. In: Kaliski Jr., B.S., Koç, Ç.K., Paar, C. (eds.) CHES 2002. LNCS,
vol. 2523, pp. 276–290. Springer, Heidelberg (2003)
23. Yao, A.C.-C.: On the Evaluation of Powers. SIAM J. Comput. 5(1), 100–103 (1976)

Appendix A. Pseudo-Code for the Recoding Scheme
The algorithm in Fig. 1 determines the base/digit pairs to associate with eachresidue i mod π of the key k during generation of the multibase representation.
Its inputs include a choice of base set B and digit set D suitable for the target
device, the window width λ and the value of π (the λth power of the lcm of basesin B times a power of 2 determined by the table size). The output is deposited in
two arrays, Base and Digit. The function Cost is that in (4), and it is describedin detail in §6.1.
The algorithm in Fig. 2 uses the arrays Base and Digit to recode the key k,storing the resulting multibase representation in an array MB. Usually the array
accesses in the loop body will be replaced by coded rules as in Appendix B.

Appendix B. Code for the Iterative Step
The very simple code of Figure 3 provides mod 26 32 rules for selecting the next
divisor/ residue pair (r, d) for k with base set B = {2, 3} and a precomputed
table of 8 points containing 1, 3, 5,..., 15 times the input point. Once the pair(r, d) is chosen, k is then reduced to (k−d)/r for the next application of therules. The rules were obtained using a sub-chain length of λ = 2 with c = 7.328

74

C.D. Walter
For i: 0 ≤ i < π do
{ BestCost ← MaxCost
For all r1 in B and all d1 in D do
If (i-d1) mod r1 = 0 then
{ i1 = (i-d1)/r1
For all r2 in B and all d2 in D do
If (i1-d2) mod r2 = 0 then
{ ...
If Cost((r1,d1),(r2,d2),...,(rλ,dλ)) < BestCost then{ BestCost ← Cost((r1,d1),(r2,d2),...,(rλ,dλ))Base[i] ← r1
Digit[i] ← d1}}}}
Fig. 1. Algorithm to determine Base/Digit pairs for each Residue mod πn ← 0While k = 0 do
{ i ← k mod π
MB[n].r ← Base[i]MB[n].d ← Digit[i]k ← (k-Digit[i])/Base[i]n ← n+1}
MB.len ← nFig. 2. Algorithm to recode k to a Division Chain Representation
If k = 0 mod 9 and k = 0 mod 4 thenr ← 3, d ← 0else if k = 0 mod 2 thenr ← 2, d ← 0else if k = 0 mod 3 and 18 < (k mod 64) < 46
and ((k mod 64)-32) = 0 mod 3 thenr ← 3, d ← 0else
r ← 2, d ← ((k+16) mod 32) - 16.
Fig. 3. A Choice of Mod 26 32 Recoding Rules for a Table Size of 8

in (4) and the Jacobi Quartic form costs of the composite operations as listed inTable 1. This generated the same value for c in (5). Increasing λ by 1 only saves
1 more multiplication per 160-bit scalar.
Similar conditions are given in [17], Tables 2 and 3. The main diﬀerence hereis that some multiples of 2 are assigned the divisor 3 rather than the 2 forcedby [17]. This enables the next residue to be modiﬁed by a large table entry to
make it much more highly divisible by 2. Note, however, that divisor 3 is notFast Scalar Multiplication for ECC over GF(p) Using Division Chains75

If k = 0 mod 9 and k = 0 mod 4
and (16 < (k mod 256) < 240) thenr ← 3, d ← 0else if k = 0 mod 2 thenr ← 2, d ← 0else if k = 0 mod 3 and 8 < (k mod 32) < 24
and ((k mod 32)-16) = 0 mod 3 thenr ← 3, d ← 0else
r ← 2, d ← ((k+8) mod 16) - 8.
Fig. 4. A Choice of Mod 28 32 Recoding Rules for a Table Size of 4

used with a non-zero digit although the generating algorithm could have allowed
it. Including the pre-computations for table entries, the total ﬁeld multiplication
count for 160-bit scalars is just under 1208 (at 1 squaring = 0.8 multiplications).
This was obtained by a Monte Carlo simulation. The optimisations referred toin [17,16] enable a further 10-12 multiplications to be saved during initialisation
using [8], §4.2.2, thereby taking the total to under 1200.
In the same context as the above except for having a pre-computed table offour instead of eight points, the code in Fig. 4 works with π = 28 32 and generates
a recoding which requires 1214.2 ﬁeld multiplications on average. This is over
22.7 multiplications better than code with no triplings.