ANALYSIS OF THE FLEXIBLE

WINDOW POWERING ALGORITHM

Henri Cohen

Abstract. We analyze precisely the parameters entering into the flexible window powering method, and apply this to the case of an elliptic curve over a large finite prime field.

1. The Flexible Window Powering Method

Let G be a group, let g ∈ G and let N be a positive integer. We want to compute efficiently g^{N}. Several issues enter into his problem, in particular the representation used for the elements of G and the operations between them, and the addition or addition-subtraction chain (if the inverse in G can be computed inexpensively) used for N. In this paper we will mainly be concerned with the latter problem.

If N is fixed and not too large, we can try to find the optimal addition or addition-subtraction chain for N. However, if N is large or even if N is not too large but is allowed to vary, this is not practical anymore and we must use reasonable approximations to the optimal chain.

A first quite reasonable idea is the use of the very simple but efficient binary algorithm (see [Knu] or [Coh0]). Denote by D the cost of a squaring in G, by A the cost of a multiplication in G, and let n = 1 + ⌊log_{2} N⌋ be the number of bits of N. The total cost of computing N · g by the binary algorithm is approximately equal to nD + (n/2)A. Note that we use the letters D (for doubling) and A (for addition) as if the group was written additively, because we reserve the use of the letters S (for squaring) and M (for multiplication) for more basic operations such as base field operations if for example G is an elliptic curve over a field.

It is well known that the binary algorithm can be improved by using a base 2^{e} for a suitable e (see for example [Coh0, Algorithm 1.2.4]). The idea is as follows: we write N in base 2^{e} and we precompute g^{1}, g^{3}, …, g^{2^{e}−1}. If

N = (((x_{v}2^{e} + x_{v−1})2^{e} + x_{v−2})2^{e} + ⋯ + x_{1})2^{e} + x_{0} ,

with 0 ≤ x_{i} < 2^{e}, we compute g^{N} by iterating operations of the form a ← a^{2^{e}}g^{x_{i}}. If x_{i} = 2^{k_{i}}y_{i} with y_{i} odd, this is done by writing a ← (a^{2^{e−k_{i}}} g^{y_{i}})^{2^{k_{i}}}, and since the g^{y_{i}} have been precomputed, this requires e squarings and 1 multiplication in G. Thus the total number of squarings in the computation of g^{n} is still approximately equal to n (it is in fact slightly less), but the total number of multiplications is now equal

Typeset by AMS-TEX

Key words and phrases. binary powering, addition chain, window method.

1991 AMS Classification 11Y40

to 2^{e−1} + v. Since v is approximately equal to n/e, we obtain a total computing cost of approximately nD + (n/e + 2^{e−1})A.

It is clear that this becomes rapidly better than the binary algorithm. In fact, choosing simply e = 3, we have n/3 + 4 < n/2 as soon as n ≥ 25. It is also not difficult to find the optimal value of e for given n by minimizing the function n/e + 2^{e−1}. A similar analysis is valid if we allow additions and subtractions.

The above method, consisting essentially in looking at the bits of N through a window of width e (some authors call w = e − 1 the width of the window, but we will stick to e here), is naturally called the window method. A perhaps surprising fact (which will be included in a future reprinting of [Coh0]), is that there is a very easy modification of the window method (which may be called the flexible window method) which leads to a better performance with essentially no extra work. This method is based on the following lemmas.

Lemma 1.1. Let e ≥ 1. There exist a unique representation of N as

N = ((((x_{v}2^{k_{v}} + x_{v−1})2^{k_{v−1}} + x_{v−2})2^{k_{v−2}} + ⋯ + x_{1})2^{k_{1}} + x_{0})2^{k_{0}} ,

satisfying 1 ≤ x_{i} ≤ 2^{e} − 1, x_{i} odd, k_{i} ≥ e for i > 0.

When subtractions are also allowed, this lemma is modified as follows.

Lemma 1.2. Let e ≥ 1. There exist a unique representation of N as

N = ((((x_{v}2^{k_{v}} + x_{v−1})2^{k_{v−1}} + x_{v−2})2^{k_{v−2}} + ⋯ + x_{1})2^{k_{1}} + x_{0})2^{k_{0}} ,

satisfying −2^{e−1} + 1 ≤ x_{i} ≤ 2^{e−1} − 1, x_{i} odd, x_{v} > 0, and k_{i} ≥ e for i > 0.

Both lemmas are easy to prove, and it is also easy to give algorithms for finding the x_{i} and k_{i}.

The point of Lemmas 1.1 and 1.2 is that the number of group squarings and the number of group multiplications during the precomputations will be the same, but on the other hand a cursory analysis (see for example [CoMiOn1]) shows that the average order v of group multiplications during the main computation will approximately be equal to n/(e + 1) instead of n/e. The purpose of this paper is to analyse precisely the behavior of this flexible window algorithm. The analysis is useful also to analyze further more marginal improvements, such as the computation of the initial quantity g^{N_{v−1}} with N_{v−1} = x_{v}2^{k_{v}} + x_{v−1}. Since the context in which this problem appeared was the case where G is the group of points of an elliptic curve over a finite field, for which the cost of subtraction is identical to the cost of addition, we will give the detailed analysis for the addition-subtraction algorithm based on Lemma 1.2. The analysis for the addition-only algorithm based on Lemma 1.1 is essentially identical, and we will only give the results.

2. Basic Mathematical Analysis

We fix a positive integer e (the window size), a positive integer n, and we let N range uniformly among the integers having exactly n bits, in other words 2^{n−1} ≤ N < 2^{n}. Let

N = ((((x_{v}2^{k_{v}} + x_{v−1})2^{k_{v−1}} + x_{v−2})2^{k_{v−2}} + ⋯ + x_{1})2^{k_{1}} + x_{0})2^{k_{0}}

be the unique representation of N given by Lemma 1, let d(N) = ∑_{0≤i≤v} k_{i} be the number of group squarings, and let a(N) = v be the number of group multiplications, so that the cost of the flexible window algorithm is D (for the computation of g^{2}) plus (2^{e−2} − 1)A (for the computation of g^{3}, g^{5}, …, g^{2^{e−1}−1}) plus d(N)D + a(N)A (for the main computation), hence in total a cost of

C(N) = (d(N) + 1)D + (a(N) + 2^{e−2} − 1)A .

We will analyze the quantities d(N) and a(N) (and similar quantities) in turn. We first need preliminary results on a special type of linear recurrence.

2.1. Analysis of the recurrence 2u(n + e) − u(n + e − 1) − u(n) = 0 The purpose of this section is to prove the following proposition.

Proposition 2.1. Let u(n) be a sequence satisfying the linear recursion 2u(n + e) − u(n + e − 1) − u(n) = 0 for n ≥ 1. Then there exists ρ > 1 such that as n → ∞

u(n) = \frac{1}{e + 1} (2u(e) + ∑_{i=1}^{e−1} u(i)) + O(ρ^{−n}) .

Proof. Let P_{e}(X) = 2X^{e}−X^{e−1}−1 be the characteristic polynomial of the recursion. Note that P_{e}(1) = 0. Denote by α_{1} = 1, α_{2}, …,α_{e} the roots of P_{e}(X). We have the following lemma.

Lemma 2.2. The α_{i} are distinct, and for 2 ≤ i ≤ e we have |α_{i}| < 1.

Proof. The roots of P′_{e}(X) are X = 0 and X = (e − 1)/(2e). Clearly P_{e}(0) ≠ 0 so 0 is not a root of P_{e}(X), and P_{e}((e − 1)/(2e)) = −(1/e)((e − 1)/(2e))^{e−1} − 1 < 0, so (e − 1)/(2e) is also not a root of P_{e}(X), hence all the roots of P_{e}(X) are simple and the α_{i} are thus distinct. Let us show that |α_{i} < 1 for i ≥ 2. Indeed, if |α_{i}| > 1, then |(2α_{i}−1)α_{i}^{e−1}| > (2|α_{i}|−1) > 1. But even if |α_{i}| = 1, then |(2α_{i}−1)α_{i}^{e−1}| = |2α_{i}−1| and if we set α_{i} = e^{iθ} this is equal to \sqrt{5 − 4 cos(θ)} which can be equal to 1 only if θ = 0, i.e. α_{i} = 1, proving the lemma. ⊓ ⊔

Resuming the proof of the proposition, we deduce that there exist constants A_{i} such that u(n) = A_{1} + ∑_{2≤i≤e} A_{i}α_{i}^{n}, hence by the lemma u(n) = A_{1} + O(ρ^{−n}) where ρ can be taken to be equal to the inverse of the largest modulus of the α_{i} for i > 1 (numerically we have ρ = \sqrt{2} = 1.414... for e = 3, ρ = 1.2157... for e = 4 and ρ = 1.1296... for e = 5).

It remains to compute the value of A_{1}. To do this, we must solve the system of equations ∑_{1≤j≤e} A_{j}α_{j}^{i} = u(i) for 0 ≤ i ≤ e − 1, where even though the value of u(0) may not be defined, we define it so as to satisfy the recursion, in other words we set u(0) = 2u(e) − u(e − 1). The matrix of this system is the Vandermonde matrix M = (m_{i,j}) with m_{i,j} = α_{j}^{i−1}. The inverse of the Vandermonde matrix is well known (see [Knu]) and is equal to M^{−1} = (n_{i,j}) with

n_{i,j} = \frac{[(α_{1} − X) ⋯ { \atop (α_{i} − X)} ⋯ (α_{e} − X)]_{j−1}}{(α_{1} − α_{i}) ⋯ { \atop (α_{i} − α_{i})} ⋯ (α_{e} − α_{i})}

where means that the corresponding term is omitted and [Q(X)]_{j−1} is the coefficient of X^{j−1} in the polynomial Q(X).

Thus in particular we have the formula

A_{1} = ∑_{j=1}^{e} \frac{[(α_{2} − X) ⋯ (α_{e} − X)]_{j−1}}{(α_{2} − α_{1}) ⋯ (α_{e} − α_{1})} b_{j−1} .

Set Q_{e}(X) = P_{e}(X)/(2(X −1)). Then Q_{e}(X) = ∏_{2≤j≤e}(X −α_{j}), and since α_{1} = 1, we have

A_{1} = ∑_{j=1}^{e} \frac{Q_{e}(X)_{j−1}}{Q_{e}(1)} u(j − 1) .

Now we easily check that

Q_{e}(X) = X^{e−1} + \frac{1}{2}(X^{e−2} + ⋯ + 1)

so that Q_{e}(X)_{j−1} = 1/2 for j ≤ e − 1 and Q_{e}(X)_{j−1} = 1 for j = e. In addition, Q_{e}(1) = (e + 1)/2. Hence finally we obtain

A_{1} = \frac{2}{e + 1} (u(e − 1) + \frac{1}{2} ∑_{j=1}^{e−1} u(j − 1)) = \frac{1}{e + 1} (u(e − 1) + ∑_{j=0}^{e−1} u(j)) .

Replacing the (possibly artificial) u(0) by 2u(e) − u(e − 1), we obtain the given formula. ⊓ ⊔

2.2. Analysis of d(N)

Set

f_{d}(n) = \frac{1}{2^{n−2}} ∑_{2^{n−1}≤N<2^{n}, N odd} d(N)

and

g_{d}(n) = \frac{1}{2^{n−1}} ∑_{2^{n−1}≤N<2^{n}} d(N) .

These are the averages of d(N) on n-bit odd numbers and all n-bit numbers respectively. Although we are only interested in the function g_{d}(n), the auxiliary function f_{d}(n) will be very useful.

First note the following. If 2^{n−1+e} ≤ N < 2^{n+e} is odd, we can write in a unique way N = q2^{e} + r with −2^{e−1} + 1 ≤ r ≤ 2^{e−1} − 1 odd, and we have either 2^{n−1} < q < 2^{n}, or q = 2^{n} and r < 0, or q = 2^{n−1} and r > 0. In addition, we have d(N) = d(q) + e. It follows that

2^{n+e−2}f_{d}(n + e) = 2^{e−1} ∑_{2^{n−1}<q<2^{n}}(e + d(q)) + 2^{e−2}(e + d(2^{n−1})) + 2^{e−2}(e + d(2^{n}))

= 2^{e−1}(2^{n−1}e + ∑_{2^{n−1}≤q<2^{n}} d(q)) + 2^{e−2}

= 2^{e+n−2}e + 2^{e−2} + 2^{e+n−2}g_{d}(n) , so that

f_{d}(n + e) = g_{d}(n) + e + 2^{−n} .

On the other hand, since an integer is either even or odd and d(2M) = d(M) + 1, we have

2^{n−2}f_{d}(n) = ∑_{2^{n−1}≤N<2^{n}} d(N) − ∑_{2^{n−2}≤N<2^{n−1}}(d(N) + 1)
= 2^{n−1}g_{d}(n) − 2^{n−2}g_{d}(n − 1) − 2^{n−2}

so that

f_{d}(n) = 2g_{d}(n) − g_{d}(n − 1) − 1 .

Replacing in the preceding recursion, we obtain the recursion involving only the function g_{d}(n), which is the one we are interested in:

2g_{d}(n + e) − g_{d}(n + e − 1) − g_{d}(n) = 1 + e + 2^{−n} .

This recurrence relation determines g_{d}(n) entirely as soon as we know the values of g_{d}(1), …,g_{d}(e). If N ≤ 2^{e−1} − 1 is odd, then d(N) = 0 hence f_{d}(n) = 0 for n ≤ e − 1. Since clearly g_{d}(1) = 0, the recursion 2g_{d}(n) − g_{d}(n − 1) − 1 = f_{d}(n) = 0 for 2 ≤ n ≤ e − 1 implies that g_{d}(n) = 1 − 2^{1−n} for 1 ≤ n ≤ e − 1. In addition, if 2^{e−1} < N < 2^{e} is odd, then N = 1 · 2^{e} + x_{0} with −2^{e−1} + 1 ≤ x_{0} ≤ −1 is the representation of Lemma 1.2, so that d(N) = e, hence f_{d}(e) = e and so 2g_{d}(e) = f_{d}(e) + g_{d}(e − 1) + 1 = e + 1 + 1 − 2^{2−e}, hence g_{d}(e) = (e + 2)/2 − 2^{1−e}.

We summarize what we have just proved in the following

Lemma 2.3. If g_{d}(n) is the average of the number d(N) of group squarings for n-bit numbers, then g_{d}(n) satisfies the recursion

2g_{d}(n + e) − g_{d}(n + e − 1) − g_{d}(n) = 1 + e + 2^{−n}

with the initial conditions g_{d}(n) = 1 − 2^{1−n} for 1 ≤ n ≤ e − 1 and g_{d}(e) = (e + 2)/2 − 2^{1−e}.

We now solve this linear recurrence explicitly. If we set h(n) = g_{d}(n) − n + 2^{−n}, then h(n) satisfies the homogeneous reccurence 2h(n+e)−h(n+e−1)−h(n) = 0 with the initial conditions h(n) = 1−n−2^{−n} for 1 ≤ n ≤ e−1 and h(e) = (2−e)/2−2^{−e}.

It follows from Proposition 2.1 that the solution h(n) to our recurrence relation satisfies h(n) = A + O(ρ^{−n}) for some ρ > 1, with

(e + 1)A = 2 − e − 2^{1−e} + ∑_{j=1}^{e−1}(1 − j − 2^{−j}) = −e(e − 1)/2

so that we have proved:

Proposition 2.4. Let g_{d}(n) be the average of the number d(N) of group squarings on n-bit numbers. There exists a real number ρ > 1 such that

g_{d}(n) = n − \frac{e(e − 1)}{2(e + 1)} + O(ρ^{−n}) .

2.3. Analysis of a(N)

The analysis will be very similar to that of d(N), so we can be quite brief. Set

f_{a}(n) = \frac{1}{2^{n−2}} ∑_{2^{n−1}≤N<2^{n}, N odd} a(N)

and

g_{a}(n) = \frac{1}{2^{n−1}} ∑_{2^{n−1}≤N<2^{n}} a(N) .

These are the averages of a(N) on n-bit odd numbers and all n-bit numbers respectively.

Since a(2M) = a(M) we immediately obtain the relation 2^{n−2}f_{a}(n) = 2^{n−1}g_{a}(n)− 2^{n−2}g_{a}(n − 1) so f_{a}(n) = 2g_{a}(n) − g_{a}(n − 1). On the other hand, we have a(q2^{e} + r) = a(q) + 1 and a(r) = 0 for −2^{e−1} + 1 ≤ r ≤ 2^{e−1} − 1, r odd. It follows by writing N = q2^{e} + r that

2^{n+e−2}f_{a}(n + e) = 2^{e−1} ∑_{2^{n−1}<q<2^{n}}(a(q) + 1) + 2^{e−2}(1 + a(2^{n−1})) + 2^{e−2}(1 + a(2^{n}))
= 2^{e−1}(2^{n−1} + 2^{n−1}g_{a}(n)) + 2^{e−2}(a(2^{n}) − a(2^{n−1}))
= 2^{e+n−2}(1 + g_{a}(n)) ,

so that f_{a}(n + e) = g_{a}(n) + 1. Combining with the other formula that we have obtained, this gives the recursion 2g_{a}(n + e) − g_{a}(n + e − 1) − g_{a}(n) = 1.

On the other hand, g_{a}(n) = 0 for 1 ≤ n ≤ e − 1, while f_{a}(e) = 1 so g_{a}(e) = 1/2. We have thus proved:

Lemma 2.5. If g_{a}(n) is the average of the number a(N) of group multiplications for n-bit numbers, then g_{a}(n) satisfies the recurrence relation 2g_{a}(n + e) − g_{a}(n + e − 1) − g_{a}(n) = 1 with the initial conditions g_{a}(n) = 0 for 1 ≤ n ≤ e − 1 and g_{a}(e) = 1/2.

To solve this recursion, we set h(n) = g_{a}(n) − n/(e + 1), so that h(n) satisfies the homogeneous recurrence relation 2h(n + e) − h(n + e − 1) − h(n) = 0 with the initial conditions h(n) = −n/(e + 1) for 1 ≤ n ≤ e − 1 and h(e) = 1/2 − e/(e + 1). Applying Proposition 2.1, this gives:

Proposition 2.6. Let g_{a}(n) be the average of the number a(N) of group multiplications on n-bit numbers. There exists a real number ρ > 1 such that

g_{a}(n) = \frac{n}{e + 1} − \frac{(e − 1)(e + 2)}{2(e + 1)^{2}} + O(ρ^{−n}) .

Corollary 2.7. The average cost of computing g^{N} for an n-bit number N using the addition-subtraction chain coming from Lemma 1.2 is approximately equal to

C(N) = (n + 1 − \frac{e(e − 1)}{2(e + 1)}) D + (2^{e−2} − 1 + \frac{n}{e + 1} − \frac{(e − 1)(e + 2)}{2(e + 1)^{2}}) A .

We will see below a more precise estimate in the case of elliptic curves.

2.3. Analysis of x_{v}(N)

Let z be an odd integer such that 1 ≤ z ≤ 2^{e−1}. We want to estimate the probability that the leading coefficient x_{v} is equal to z. Contrary to what intuition could suggest, this probability is not even close to being uniformly distributed among all values of z. In fact, we will find the well known phenomenon that the initial digit of a number is more often equal to 1 than to any other digit, in an appropriate sense.

For simplicity of notation, define δ_{z}(N) = 1 if x_{v}(N) = z and δ_{z}(N) = 0 otherwise. Set

f_{x,z}(n) = \frac{1}{2^{n−2}} ∑_{2^{n−1}≤N<2^{n}, N odd} δ_{z}(N)

and

g_{x,z}(n) = \frac{1}{2^{n−1}} ∑_{2^{n−1}≤N<2^{n}} δ_{z}(N) .

Since δ_{z}(2k) = δ_{z}(k) we immediately obtain the relation

2^{n−2}f_{x,z}(n) = 2^{n−1}g_{x,z}(n) − 2^{n−2}g_{x,z}(n − 1)

hence f_{x,z}(n) = 2g_{x,z}(n) − g_{x,z}(n − 1). On the other hand for q > 0 we have δ_{z}(q2^{e} + r) = δ_{z}(q) , hence it follows by writing N = q2^{e} + r that

2^{n+e−2}f_{x,z}(n + e) = 2^{e−1} ∑_{2^{n−1}<q<2^{n}} δ_{z}(q) + 2^{e−2} + 2^{e−2}

= 2^{e−1} ∑_{2^{n−1}≤q<2^{n}} δ_{z}(q)

= 2^{e+n−2}g_{x,z}(n) ,

so that f_{x,z}(n + e) = g_{x+z}(n). Combining with the other formula that we have obtained, this gives the recursion 2g_{x,z}(n + e) − g_{x,z}(n + e − 1) − g_{x,z}(n) = 0.

We now consider two cases.

Case 1: z = 1

Then we have δ_{z}(r) = 0 for 3 ≤ r < 2^{e−1}, r odd, δ_{z}(1) = 1 and δ_{z}(r) = 1 for 2^{e−1} < r < 2^{e}, r odd. Thus, f_{x,z}(1) = 1, f_{x,z}(n) = 0 for 2 ≤ n ≤ e − 1, f_{x,z}(e) = 1, hence g_{x,z}(n) = 1/2^{n−1} for 1 ≤ n ≤ e − 1, while g_{x,z}(e) = 1/2 + 1/2^{e−1}.

Using Proposition 2.1, we obtain that for z = 1 the average of δ_{z}(N) on n-bit numbers is equal to 3/(e + 1) + O(ρ^{−n}) for some ρ > 1.

Case 2: z > 1

Let m be the number of bits of z, so that 2^{m−1} ≤ z < 2^{m}. We have 2 ≤ m ≤ e−1. Then δ_{z}(r) = 0 for 1 ≤ r < 2^{e}, r odd, except for r = z. Thus, f_{x,z}(n) = 0

for 1 ≤ n ≤ e, n ≠ m, and f_{x,z}(m) = 1/2^{m−2}, hence g_{x,z}(n) = 0 for n < m − 1, g_{x,z}(n) = 1/2^{n−1} for m ≤ n ≤ e. Using Proposition 2.1, we obtain that for 2^{m−1} ≤ z < 2^{m}, 2 ≤ m ≤ e − 1, the average of δ_{z}(N) on n-bit numbers is equal to 1/(2^{m−2}(e + 1)) + O(ρ^{−n}) for some ρ > 1. To summarize, we can state:

Proposition 2.8. Let f_{x,z}(n) be the probability that the leading coefficient x_{v}(N) be equal to an odd nummber z for n-bit numbers. There exists a real number ρ > 1 such that

f_{x,z}(n) = \{ \begin{matrix}\frac{3}{e + 1} & if z = 1 \\ \frac{1}{2^{m−2}(e + 1)} & if z has exactly m bits. \end{matrix}

In particular, the probability that x_{v}(N) has 1 bit is equal to 3/(e + 1) while the probability that x_{v}(N) has exactly m bits is equal to 1/(e + 1) for m ≥ 2.

2.4. Analysis of s(N)

For a fine analysis of the window method, we must also know the average number of ±1 among the x_{i}(N). Thus, let us let

s(N) = ∑_{0≤i≤v, x_{i}(N)=±1} 1

be the number of ±1 in N.

f_{s}(n) = \frac{1}{2^{n−2}} ∑_{2^{n−1}≤N<2^{n}, N odd} s(N)

and

g_{s}(n) = \frac{1}{2^{n−1}} ∑_{2^{n−1}≤N<2^{n}} s(N) .

Since s(2k) = s(k) we immediately obtain the relation 2^{n−2}f_{s}(n) = 2^{n−1}g_{s}(n) − 2^{n−2}g_{s}(n − 1) so f_{s}(n) = 2g_{s}(n) − g_{s}(n − 1).

To simplify notation, set δ(r) = 1 if r = ±1, δ(r) = 0 otherwise. Then s(q2^{e} + r) = s(q) + δ(r) and s(r) = δ(r) for 1 ≤ r < 2^{e−1}, r odd, s(r) = 1 for 2^{e−1} < r ≤ 2^{e} − 3, r odd, and s(2^{e} − 1) = 2. It follows by writing N = q2^{e} + r that

2^{n+e−2}f_{s}(n + e) = ∑_{−2^{e−1}<r<2^{e−1}, odd} ∑_{2^{n−1}<q<2^{n}}(s(q) + δ(r))
+ ∑_{1≤r<2^{e−1}, odd}(1 + δ(r)) + ∑_{−2^{e−1}<r≤−1, odd}(1 + δ(r))

= ∑_{−2^{e−1}<r<2^{e−1}, odd} ∑_{2^{n−1}≤q<2^{n}}(s(q) + δ(r))

= 2^{e+n−2}g_{s}(n) + 2^{n} ,

so that f_{s}(n+e) = g_{s}(n)+1/2^{e−2}. Combining with the other formula that we have obtained, this gives the recursion 2g_{s}(n + e) − g_{s}(n + e − 1) − g_{s}(n) = 1/2^{e−2}.

On the other hand, f_{s}(1) = 2, f_{s}(n) = 0 for 2 ≤ n ≤ e−1, f_{s}(e) = 1+1/2^{e−2}, so g_{s}(n) = 1/2^{n−1} for 1 ≤ n ≤ e − 1 and g_{s}(e) = 1/2 + 1/2^{e−2}. We have thus proved:

Lemma 2.9. If g_{s}(N) is the average number of ±1 for n-bit numbers, then g_{s}(N) satisfies the recurrence relation 2g_{s}(n + e) − g_{s}(n + e − 1) − g_{s}(n) = 1/2^{e−2} with the initial conditions g_{s}(n) = 1/2^{n−1} for 1 ≤ n ≤ e − 1 and g_{s}(e) = 1/2 + 1/2^{e−2}.

To solve this recursion, we set h(n) = g_{s}(n) − n/(2^{e−2}(e + 1)), so that h(n) satisfies the homogeneous recurrence relation 2h(n + e) − h(n + e − 1) − h(n) = 0 with the initial conditions h(n) = 1/2^{n−1} − n/(2^{e−2}(e + 1)) for 1 ≤ n ≤ e − 1 and h(e) = 1/2 + 1/(2^{e−2}(e + 1)). Applying Proposition 2.1, we obtain finally:

Proposition 2.10. The average number of x_{i}(N) = ±1 for n bit numbers is equal to

\frac{n}{2^{e−2}(e + 1)} + \frac{3}{e + 1} − \frac{(e − 1)(e + 2)}{2^{e−1}(e + 1)^{2}} + O(ρ^{−n})

for some ρ > 1.

Since we have proved in the preceding section that the probability that x_{v}(N) = 1 is equal to 3/(e + 1) (x_{v}(N) = −1 is not possible), we obtain:

Corollary 2.11. The average number of x_{i}(N) = ±1 among the x_{i}(N) for 0 ≤ i ≤ v − 1 for n bit numbers is equal to

\frac{n}{2^{e−2}(e + 1)} − \frac{(e − 1)(e + 2)}{2^{e−1}(e + 1)^{2}} + O(ρ^{−n})

for some ρ > 1.

Since we may choose 2 among the 2^{e−1} possible digits, we see that this time intuition does not fool us, since the value given above is (up to O(ρ^{−n})) exactly equal to the average number of group multiplications given in Proposition 2.6 divided by 2^{e−2}.

2.5. The case of the addition-only window method

In this section, we give the corresponding results for the addition-only method using the representation of Lemma 1.1 instead of Lemma 1.2. The proofs are essentially identical to the addition/subtraction case so are omitted.

Proposition 2.12. In the case of the addition-only method, there exists ρ > 1 such that:

(1) The average number g_{d}(n) of group squarings is equal to

g_{d}(n) = n − \frac{e^{2} + e + 2}{2(e + 1)} + O(ρ^{−n}) .

(2) The average number of group multiplications is equal to

g_{a}(n) = \frac{n}{e + 1} − \frac{e(e + 3)}{2(e + 1)^{2}} + O(ρ^{−n}) .

(3) The probability f_{x,z}(n) that the leading coefficient x_{v}(N) is equal to an odd number z for n-bit numbers is equal to

f_{x,z}(n) = \frac{1}{2^{m−2}(e + 1)}

if z has exactly m bits, with 1 ≤ m ≤ e. In particular, the probability that x_{v}(N) has 1 bit is equal to 2/(e + 1) while the probability that x_{v}(N) has exactly m bits is equal to 1/(e + 1) for 2 ≤ m ≤ e.

(4) The average number of x_{i} = 1 for n bit numbers is equal to

\frac{n}{2^{e−1}(e + 1)} + \frac{2}{e + 1} − \frac{e(e + 3)}{2^{e}(e + 1)^{2}} + O(ρ^{−n}) .

3. Application to elliptic curve exponentiation

In [CoMiOn2], we have presented an improved method for elliptic curve exponentiation using the flexible window method together with two new ideas: the use of a new modified coordinate system J^{m}, and the use of mixed coordinates for a single computation. In the present section, we use the results of the preceding section together with the results of [CoMiOn2] to precisely analyse the method.

Let 𝔽_{p} be a finite prime field, where p is assumed to be large (for example at least 96 bits), and let E be an elliptic curve over 𝔽_{p} whose affine equation is

y^{2} = x^{3} + ax + b .

Apart from the choice of affine coordinates, denoted A, we also use the Jacobian coordinates (x, y, z) corresponding to the affine point (x/z^{2}, y/z^{3}) when z ≠ 0, denoted J, and the modified Jacobian coordinates (x, y, z, az^{4}) with the same meaning but an extra component, denoted J^{m} (in [CoMiOn2] we also considered ordinary projective coordinates and Chudnovsky Jacobian coordinates, but they are not useful in practice).

3.1. The ordinary Horner scheme

The strategy given in [CoMiOn2] is the following, assuming reasonable practical hypotheses about the respective times for inversion I, multiplication M and squaring S in the base field 𝔽_{p}.

Let P be a point on E given in affine coordinates. To simplify the analysis, we assume that we are in a generic situation so that no point encountered during the computation is the point at infinity. This is of course highly likely since the prime p has been assumed large.

To compute N · P for some large integer N using the base 2^{e} flexible window representation, we first precompute 2P in affine coordinates, then the points x · P in affine coordinates for 3 ≤ x ≤ 2^{e−1} − 1, x odd by 2^{e−2} − 1 addditions in affine coordinates. Since, with evident notations, we have t(A ← A + A) = I + 2M + S and t(A ← 2A) = I + 2M + 2S, the precomputations require a time T_{1} given by

T_{1} = S + 2^{e−2}(I + 2M + S) .

For the main computation, we use Horner’s scheme on the representation

N = ((((x_{v}2^{k_{v}} + x_{v−1})2^{k_{v−1}} + x_{v−2})2^{k_{v−2}} + ⋯ + x_{1})2^{k_{1}} + x_{0})2^{k_{0}} .

Thus, if we set

N_{j} = ((x_{v}2^{k_{v}} + x_{v−1})2^{k_{v−1}} + ⋯ + x_{j+1})2^{k_{j+1}} + x_{j} ,

and P_{j} = N_{j} · P, we initialize P_{v} ← x_{v} · P and we set for j = v − 1, …,j = 0,

P_{j} ← 2^{k_{j+1}} · P_{j+1} + x_{j} · P

and the final result is

N · P = 2^{k_{0}} · P_{0} .

The computation of P_{j} ← 2^{k_{j+1}} · P_{j+1} + x_{j} · P is performed as follows. By induction P_{j+1} will be in the modified Jacobian coordinates J^{m}. The k_{j+1} − 1 doublings will be done using these coordinates, while the last such doubling will use these coordinates but give the result in ordinary Jacobian coordinates. The following times are taken from [CoMiOn2]:

t(J^{m} ← 2J^{m}) = 4M + 4S ,
t(J ← 2J^{m}) = 3M + 4S ,
t(J^{m} ← 2A) = 3M + 4S ,
t(A ← 2J^{m}) = I + 6M + 5S ,
t(J^{m} ← J + A) = 9M + 5S ,
t(J^{m} ← A + A) = 5M + 4S ,
t(A ← J + A) = I + 11M + 4S .

Thus the computation of P_{j} ← 2^{k_{j+1}} · P_{j+1} + x_{j} · P for 1 ≤ j ≤ v − 2 requires a time t_{2,j} given by

t_{2,j} = k_{j+1}(4M + 4S) + 8M + 5S .

For j = 0, the time is given by the same formula of k_{0} > 0, but we can subtract M + 2S if k_{0} = 0 since we need the result only in Jacobian, but not modified Jacobian coordinates.

For j = v − 1, 2P has already been computed, in affine coordinates, so the computation time is

t_{2,v−1} = k_{v}(4M + 4S) + 3M + S .

Finally, assuming we want the result in affine coordinates, the final k_{0} doublings require the time for conversion from Jacobian to affine if k_{0} = 0, that is I +3M +S, and otherwise time

t_{2,−1} = k_{0}(4M + 4S) + I + 2M + S .

Thus, the time for Horner’s scheme is

T_{2} = ( ∑_{0≤j≤v} k_{j}) (4M + 4S) + I + v(8M + 5S) − 3M − (3 + 2δ(k_{0})S) ,

where δ(x) is the Kronecker symbol equal to 0 except when x = 0 for which it is equal to 1.

Using the notation used in the preceding section, this can be written

T_{2}(N) = d(N)(4M + 4S) + I + a(N)(8M + 5S) − 3M − (3 + 2δ(k_{0})S) .

Since the probability for k_{0} = 0 is equal to 1/2 (no need for a detailed analysis for that), we deduce from Section 2 that the average of T_{2}(N) over integers N having exactly n bits is approximately equal to

(n − \frac{e(e − 1)}{2(e + 1)}) (4M + 4S) + (\frac{n}{e + 1} − \frac{(e − 1)(e + 2)}{2(e + 1)^{2}}) (8M + 5S) + I − 3M − 4S .

3.2. Horner’s scheme with initialization

As already mentioned, it is possible to slightly improve on the first step of Horner’s scheme. For example, if x_{v} = 1 and k_{v} = e, instead of computing 2^{k_{v}} · P which requires k_{v} doublings, we can simply compute 2 · (P + (2^{e−1} − 1)P) which requires only one doubling and one addition, since (2^{e−1} − 1)P has been precomputed. This is especially attractive since by the results of Section 2 we know that x_{v} will be equal to 1 with high probability.

We could very precisely write the optimal addition chain to choose for the initialization step. Since the gain would be minimal (less than M), we will be content with the following approximation to optimality.

If 2^{m−1} ≤ x_{v} < 2^{m}, in other words if x_{v} has exactly m bits, with 1 ≤ m ≤ e − 1, we compute P_{v−1} ← x_{v}2^{k_{v}} ·P in ordinary Jacobian coordinates by using the identity

x_{v}2^{k_{v}} = 2^{k_{v}−(e−m)}((2^{e−1} − 1) + (x_{v}2^{e−m} + 1 − 2^{e−1})) .

Since 2^{m−1} ≤ x_{v} ≤ 2^{m} − 1 and m ≤ e − 1, we have

1 ≤ x_{v}2^{e−m} + 1 − 2^{e−1} ≤ 2^{e−1} − 2^{e−m} + 1 ≤ 2^{e−1} − 1

so the computation of x_{v}2^{k_{v}} using this method requires only k_{v} −(e−m) doublings and one addition, giving a gain of e − m doublings minus 1 addition.

More precisely, the usual method for computing x_{v}2^{k_{v}} · P giving the result in ordinary Jacobian coordinates requires time (k_{v} − 1)(4M + 4S) − 2M since 2P has been computed.

The above method requires time (k_{v} − e + m + 1)(4M + 4S). The improvement is thus equal to (e − m − 2)(4M + 4S) − 2M, so should be used when m ≤ e − 3. In particular, there is no possible gain when e ≤ 3.

It is easy to compute the average gain when e ≥ 4: by Proposition 2.8, we have x_{v} = 1 with probability 3/(e + 1) and x_{v} has m bits with probability 1/(e + 1) for 2 ≤ m ≤ e − 1. Since we use this method only for m ≤ e − 3, a small computation shows that the average gain is equal to

\frac{2}{e + 1}((e^{2} − 2e − 5)M + (e^{2} − e − 6)S) .

3.3. Using Montgomery’s trick

There is still another trick which can be used to speed up the computation. This is the use of a variant of Montgomery’s method for computing simultaneous inverses, applied to the precomputations which are necessary in the window method for elliptic curves. Recall that, if k inversions modulo p are to be performed simultaneously, using Montgomery’s method we may compute a single inverse modulo p at the cost of 3(k − 1) multiplications (see for example [Coh0, Algorithm 10.3.4]).

As explained above, we first compute 2P in affine coordinates. But now, we compute simultaneously 3P = 2P +P and 4P = 2P +2P (still in affine coordinates) using Montgomery’s idea, then simultaneously 5P = 4P + P, 7P = 4P + 3P, 8P = 4P + 4P, and so on until (2^{e−2} + 1)P = 2^{e−2}P + P, …, (2^{e−1} − 1)P = 2^{e−2}P + (2^{e−2} − 1)P and we do not compute 2^{e−1}P. Although we compute the apparently useless quantities 4P, 8P, …, 2^{e−2}P, a short calculation shows that the time for computing x · P for x odd, 3 ≤ x ≤ 2^{e−1} − 1 using this method is equal to

(e − 1)I + (5 · 2^{e−2} + 2e − 12)M + (2^{e−2} + 2e − 5)S

when e ≥ 3. For example, with e = 5 this gives 4I + 38M + 13S, compared with T_{1} = 8I + 16M + 9S using the ordinary method, and this is almost always faster since I is much larger than M.

References

[Coh0] H. Cohen, A Course in Computational Algebraic Number Theory, GTM 138, Springer-Verlag, Berlin, Heidelberg, New York, 1993.

[CoMiOn1] H. Cohen, A. Miyaji and T. Ono, Efficient elliptic curve exponentiation, Proceedings ICICS’97, Lecture Notes in Computer Science, vol. 1334, Springer-Verlag, 1997, pp. 282–290.

[CoMiOn2] H. Cohen, A. Miyaji and T. Ono, Efficient elliptic curve exponentiation using mixed coordinates, Proceedings ASIACRYPT’98, Lecture Notes in Computer Science, vol. 1514, Springer-Verlag, 1998, pp. 51–65.

[Knu] D. Knuth, The Art of Computer Programming, vol. II, third edition, Addison Wesley, 1997.

Laboratoire A2X, UPRESA 5465 du C.N.R.S., Université Bordeaux I, 351 Cours de la Libération, 33405 TALENCE Cedex, FRANCE