ńņš. 5 |

In Section VI.3.1 we showed that if (ak , bk )ā is the AGM sequence, then

k=0

j (Eak ,bk ) ā” j(Ek ) (mod 2k+1 ) ,

where Ek is the canonical lift of Ļ k (E). Expressing j(Eak ,bk ) as a function of

ak and bk shows that ak and bk must be correct modulo 2k+3 . Therefore, we

conclude that

amā’3 amā’3+n

Tr F ā” + 2n (mod 2m ) .

amā’3+n amā’3

This procedure is summarized in Algorithm VI.5.

120 VI. POINT COUNTING

Algorithm VI.5: AGM

An elliptic curve E : y 2 + xy = x3 + c over F2n with

INPUT:

j(E) ā F4 .

OUTPUT: The number of points on E(F2n ).

1. m ā n + 2. ā

2

2. Take r ā Zq with r ā” c (mod 2).

3. a ā 1 + 4r, b ā 1 ā’ 4r.

4. For i = 4 to m do:

ā

(a, b) ā (a + b)/2, ab (mod 2i ).

5.

6. a0 ā a.

7. For i = 0 to n ā’ 1 do:

ā

(a, b) ā (a + b)/2, ab (mod 2m ).

8.

9. t ā a0 /a (mod 2m ).

10. If t2 > 2n+2 then t ā t ā’ 2m .

11. Return 2n + 1 ā’ t.

The complexity of Algorithm VI.5 is determined by Step 7, which requires

O(n) square root computations to precision m n/2. Since each square root

computation takes O(1) multiplications at full precision, the complexity of

Algorithm VI.5 is O(n2Āµ+1 ) bit-operations. The space complexity clearly is

O(n2 ), since only O(1) elements of Zq /(2m Zq ) are required. Note that it is

possible to replace the For-loop in Step 7 of Algorithm VI.5 by one AGM

iteration and a norm computation. Indeed, the trace of Frobenius satisļ¬es

t ā” NQq /Qp (a0 /a1 ) (mod 2min(n,m) ). We refer to Section VI.5 for eļ¬cient norm

computation algorithms.

VI.3.4. Univariate AGM. Given the bivariate AGM sequence (ak , bk )ā k=0

with ak ā” bk ā” 1 (mod 4) and ak ā” bk (mod 8), we can easily deļ¬ne a

univariate AGM sequence (Ī»k )ā by Ī»k = bk /ak . This sequence corresponds

k=0

to the elliptic curves

EĪ»k : y 2 = x(x ā’ 1)(x ā’ Ī»2 ) .

k

ā

Since (ak+1 , bk+1 ) = ((ak + bk )/2, ak bk ), Ī»k+1 follows from Ī»k by

ā

2 Ī»k

Ī»k+1 = . (VI.12)

1 + Ī»k

In Section VI.3.1 we showed that for an ordinary curve E : y 2 + xy = x3 + c

with c ā Fā— , the AGM sequence can be initialized as (a0 , b0 ) = (1 + 4u, 1 ā’ 4u)

q

where u ā Zq satisļ¬es u2 ā” c (mod 2). However, these values cannot be

used to initialize the univariate AGM sequence, since (a0 , b0 ) is only correct

modulo 8, which would lead to Ī»0 ā” 1 (mod 8). However, this problem can

VI.4. GENERALIZED NEWTON ITERATION 121

be solved easily by computing (a1 , b1 ) ā” (1, 1 ā’ 8u2 ) (mod 16) and starting

from

Ī»1 ā” 1 + 8u2 ā” 1 + 8c (mod 16) ,

with c ā Zq and c ā” c (mod 2). Unlike the bivariate AGM sequence, the uni-

variate AGM sequence does converge, in the sense that Ī»k ā” Ī»k+n (mod 2k+3 ).

In Section VI.3.2 we showed that

bk+1 bk

ā” Ī£( ) (mod 2k+3 ) ,

ak+1 ak

and thus Ī»k+1 ā” Ī£(Ī»k ) (mod 2k+3 ). Gaudry [143] suggested substituting this

in equation (VI.12) which shows that Ī»k satisļ¬es

Ī£(Z)2 (1 + Z)2 ā’ 4Z ā” 0 (mod 2k+3 ) and Z ā” 1 + 8Ī£kā’1 (c) (mod 16) .

Deļ¬ne the polynomial Ī2 (X, Y ) = Y 2 (1 + X)2 ā’ 4X; then Ī»k is a solution

of Ī2 (X, Ī£(X)) ā” 0 (mod 2k+3 ). Recall that the j-invariant of the canonical

lift of E satisļ¬es a similar type of equation, i.e., Ī¦p (X, Ī£(X)) = 0 with

Ī¦p (X, Y ) the pth modular polynomial. The main diļ¬erence is that both

partial derivatives of Ī2 (X, Y ) vanish modulo 2. Therefore, we make the

change of variables X ā’ 1 + 8X and Y ā’ 1 + 8Y to obtain the modiļ¬ed

modular polynomial

Ī„2 (X, Y ) = (X + 2Y + 8XY )2 + Y + 4XY .

Let Ī³k be deļ¬ned by Ī»k = 1+8Ī³k ; then Ī³k satisļ¬es Ī„2 (X, Ī£(X)) ā” 0 (mod 2k )

and Ī³k ā” Ļ kā’1 (c) (mod 2). The partial derivatives of Ī„2 are

ā‚Ī„2

(X, Y ) = 2(X + 2Y + 8XY )(1 + 8Y ) + 4Y ,

ā‚X

ā‚Ī„2

(X, Y ) = (4(X + 2Y + 8XY ) + 1)(1 + 4X) ,

ā‚Y

which shows that ā‚Ī„2 ā” 0 (mod 2), but ā‚Ī„2 ā” 1 (mod 2).

ā‚X ā‚Y

The trace of Frobenius follows easily from Section VI.3.2 and is given by

1

TrF ā” NQq /Qp (mod 2min(n,k) ) . (VI.13)

1 + 4Ī³k

VI.4. Generalized Newton Iteration

In Sections VI.2.3 and VI.3.4 we solved equations of the form

Ī“(X, Ī£(X)) = 0 for Ī“(X, Y ) ā Zq [X, Y ]

without computing the Frobenius substitution Ī£. The aim of this section is

to show that in a well-chosen basis of Qq , the Frobenius substitution can be

computed eļ¬ciently, which can be used to devise faster algorithms to solve

the above equation Ī“(X, Ī£(X)) = 0.

122 VI. POINT COUNTING

VI.4.1. Fast Frobenius Substitution. The ļ¬rst construction admitting a

fast Frobenius substitution is due to Satoh, Skjernaa and Taguchi [287] and

works for all ļ¬nite ļ¬elds Fq Fp [t]/(f (t)). Let f (t) ā Zp [t] be the unique

polynomial deļ¬ned by

f (t)|(tq ā’ t) and f (t) ā” f (t) (mod p).

Then Qq can be represented as Qp [t]/(f (t)). The polynomial f (t) can be com-

puted modulo pm in O((nm)Āµ log m) time using an algorithm due to Harley;

a detailed description can be found in [333, 3.10.2].

As shown in Section VI.1, this implies that Ī£(t) = tp and the Frobenius

nā’1

substitution of an element Ī± = i=0 Ī±i ti can be computed eļ¬ciently as

nā’1 nā’1

i

Ī±i tip ,

Ī£(Ī±) = Ī£ Ī±i t =

i=0 i=0

where the result is reduced modulo f (t). Note that in general the polynomial

f (t) is dense, so a fast reduction algorithm should be used, e.g. the one

described in [141, Chapter 9].

Similarly, the inverse Frobenius substitution Ī£ā’1 can be computed as

pā’1

nā’1

ā’1 ā’1 i

Ī±pk+j tk

Ī£ (Ī±) = Ī£ Ī±i t = Cj (t),

i=0 j=0 0ā¤pk+j<n

nā’1

where Cj (t) = Ī£ā’1 (tj ) ā” tjp (mod f (t)). If we precompute Cj (t) for

ā’1

j = 0, . . . , pā’1, computing Ī£ (Ī±) for Ī± ā Zq only takes pā’1 multiplications

in Zq .

The second construction is due to Kim, Park, Cheon, Park, Kim and

Hahn [199] who proposed the use of ļ¬nite ļ¬elds with a Gaussian Normal

Basis (GNB) of small type. Such a basis can be lifted trivially to Zq and

allows eļ¬cient computation of arbitrary iterates of Frobenius. Recall that a

basis B of Qq /Qp is called normal if there exists an element Ī² ā Qq such that

B = {Ī(Ī²) | Ī ā Gal(Qq /Qp )}. A proof of the next proposition can be found

in [199].

Proposition VI.12. Let p be a prime and n, t positive integers such that

nt + 1 is a prime diļ¬erent from p. Let Ī³ be a primitive (nt + 1)-th root of

unity in some extension ļ¬eld of Qp . If gcd(nt/e, n) = 1, with e the order of

p modulo nt + 1, then for any primitive tth root of unity Ļ„ in Z/(nt + 1)Z

tā’1

i

Ī³Ļ„

Ī²=

i=0

is a normal element and [Qp (Ī²) : Qp ] = n. Such a basis is called a Gaussian

Normal Basis of type t.

VI.4. GENERALIZED NEWTON ITERATION 123

Kim et al. represent elements of Zq as elements of the ring

Zp [x]/(xnt+1 ā’ 1) .

Multiplication of two elements in Zq /(pm Zq ) therefore requires O((nmt)Āµ )

bit-operations, so one normally restricts t ā¤ 2.

For t = 1 we have Ī² = Ī³ and the minimal polynomial of Ī² therefore is

xn+1 ā’ 1

= xn + xnā’1 + Ā· Ā· Ā· + x + 1 .

f (x) =

xā’1

To speed up the Frobenius substitution, Kim et al. use a redundant represen-

nā’1

tation, i.e., they embed Zq in Zp [x]/(xn+1 ā’ 1) by mapping Ī± = i=0 Ī±i Ī² i to

k

nā’1

Ī±(x) = i=0 Ī±i xi + 0xn . Since Ī£k (Ī²) = Ī² p , we have

n n

ipk

k

Ī±j/pk (mod (n+1)) xj .

Ī£ (Ī±(x)) = Ī±i x = a0 +

i=0 j=1

Therefore, we can compute Ī£k (Ī±) by a simple permutation of the coeļ¬cients

of Ī±(x), which only requires O(n) bit-operations.

VI.4.2. Satohā“Skjernaaā“Taguchi Algorithm. Let x ā Zq be a root of

Ī“(X, Ī£(X)) = 0 for Ī“(X, Y ) ā Zq [X, Y ] and assume we have computed the

approximation xm ā” x (mod pm ). Deļ¬ne Ī“m = (x ā’ xm )/pm ; then the Taylor

series expansion around xm gives

0 = Ī“(x, Ī£(x)) = Ī“(xm + pm Ī“m , Ī£(xm + pm Ī“m ))

(VI.14)

ā” Ī“(xm , Ī£(xm )) + pm (Ī“m āx + Ī£(Ī“m )āy ) (mod p2m ) ,

with āx ā” ā‚X (xm , Ī£(xm )) (mod pm ) and āy ā” ā‚Y (xm , Ī£(xm )) (mod pm ).

ā‚Ī“ ā‚Ī“

Since Ī“(xm , Ī£(xm )) ā” 0 (mod pm ), we can divide this equation by pm and

obtain a relation for Ī“m modulo pm :

Ī“(xm , Ī£(xm ))

+ Ī“m āx + Ī£(Ī“m )āy ā” 0 (mod pm ) .

m

p

For simplicity we will assume that ordp (āy ) = 0, i.e., that āy is a unit in Zq

and that ordp (āx ) > 0. Reducing the above equation modulo p then gives

Ī“(xm , Ī£(xm ))

Ī“m ā” ā’

p

(mod p) .

pm āy

Taking the unique pth root Ī“m ā Fq leads to a better approximation of x given

by xm +pm Ī“m ā” x (mod pm+1 ). To avoid computing the pth root in the above

equation, Satoh, Skjernaa and Taguchi replace the equation Ī“(X, Ī£(X)) = 0

by Ī“(Ī£ā’1 (X), X) = 0, which clearly is equivalent. The new Ī“m then becomes

Ī“(Ī£ā’1 (xm ), xm )

Ī“m ā” ā’ m ā‚Ī“ ā’1 (mod p) .

p ā‚Y (Ī£ (xm ), xm )

124 VI. POINT COUNTING

Note that since Ī“(Ī£ā’1 (xm ), xm ) ā” 0 (mod pm ), we only need to compute the

inverse of ā‚Y (Ī£ā’1 (xm ), xm ) modulo p. Implementing these ideas naively leads

ā‚Ī“

to Algorithm VI.6, which can be used instead of Algorithm VI.4.

Algorithm VI.6: Satoh Skjernaa Taguchi Naive

Polynomial Ī“(X, Y ) ā Zq , element x0 ā Zq satisfying

INPUT:

Ī“(Ī£ā’1 (x0 ), x0 ) ā” 0 (mod p) and precision m.

OUTPUT: Element xm ā Zq , with Ī“(Ī£ā’1 (xm ), xm ) ā” 0 (mod pm ) and

xm ā” x0 (mod p).

ā’1

d ā ā‚Y (Ī£ā’1 (x0 ), x0 )

ā‚Ī“

(mod p).

1.

y ā x0 (mod p).

2.

For i = 2 to m do:

3.

x ā Ī£ā’1 (y) (mod pi ).

4.

y ā y ā’ dĪ“(x, y) (mod pi ).

5.

Return y.

6.

The main problem of Algorithm VI.6 is that in each iteration it recomputes

Ī“(x, y) although the values of x and y in Step i + 1 are very close to the ones

in Step i. Assume we have computed xW ā” x (mod pW ) for some W . Then,

for all s ā N,

Ī“ Ī£ā’1 (xsW +i ), xsW +i ā” Ī“ Ī£ā’1 (xsW ), xsW + ā (mod p(s+1)W ) , (VI.15)

with

ā‚Ī“ ā’1 ā‚Ī“ ā’1

(Ī£ (xsW ), xsW )Ī£ā’1 (Ī“) +

ā = psW (Ī£ (xsW ), xsW )Ī“ .

ā‚X ā‚Y

Furthermore, note that it is suļ¬cient to know the partial derivatives

ā‚Ī“ ā’1 ā‚Ī“ ā’1

(Ī£ (xsW ), xsW ) and (Ī£ (xsW ), xsW )

ā‚X ā‚Y

modulo pW only. Using (VI.15), we can thus compute Ī“(Ī£ā’1 (xsW +i ), xsW +i )

from Ī“(Ī£ā’1 (xsW ), xsW ) as long as i < W .

Algorithm VI.7 consists of two stages: in the ļ¬rst stage we compute

xW ā” x (mod pW ) using Algorithm VI.6 and in the second stage we use

equation (VI.15) to update Ī“(x, y).

VI.4. GENERALIZED NEWTON ITERATION 125

Algorithm VI.7: Satoh Skjernaa Taguchi

Polynomial Ī“(X, Y ) ā Zq , element x0 ā Zq satisfying

INPUT:

Ī“(Ī£ā’1 (x0 ), x0 ) ā” 0 (mod p) and precision m.

OUTPUT: Element xm ā Zq , with Ī“(Ī£ā’1 (xm ), xm ) ā” 0 (mod pm ) and

xm ā” x0 (mod p).

1. y ā Satoh Skjernaa Taguchi Naive(x0 , W ).

2. x ā Ī£ā’1 (y) (mod pW ).

3. āx ā ā‚X (x, y) (mod pW ).

ā‚Ī“

4. āy ā ā‚Y (x, y) (mod pW ).

ā‚Ī“

5. For s = 1 to (m ā’ 1)/W do:

x ā Ī£ā’1 (y) (mod p(s+1)W ).

6.

V ā Ī“(x, y) (mod p(s+1)W ).

7.

For i = 0 to W ā’ 1 do:

8.

Ī“y ā ā’ dpā’(sW +i) V (mod p).

9.

Ī“x ā Ī£ā’1 (Ī“y ) (mod pW ā’i ).

10.

y ā y + psW +i Ī“y (mod p(s+1)W ).

11.

V ā V + p(sW +i) (āx Ī“x + āy Ī“y ) (mod p(s+1)W ).

12.

13. Return y.

Satoh, Skjernaa and Taguchi prove that for W nĀµ/(1+Āµ) , Algorithm VI.7

runs in time O(nĀµ mĀµ+1/(1+Āµ) ). Although the value for W minimizes the time

complexity, in practise it is faster to take W a multiple of the CPU word size.

VI.4.3. Solving Artinā“Schreier Equations. Recall that if Fq is a ļ¬eld of

characteristic p, an equation of the form xp ā’x+c = 0 with c ā Fq is called an

Artinā“Schreier equation. Since Ī£(x) ā” xp (mod p), a natural generalization

to Zq is an equation of the form

a Ī£(x) + b x + c = 0 , (VI.16)

with a, b, c ā Zq . Let x ā Zq be a solution of Ī“(X, Ī£(X)) = 0 and let

xm ā” x (mod pm ). Deļ¬ne Ī“m = (x ā’ xm )/pm ; then equation (VI.14) shows

that Ī“m satisļ¬es an Artinā“Schreier equation modulo pm with

ā‚Ī“ ā‚Ī“ Ī“(xm , Ī£(xm ))

a= (xm , Ī£(xm )) , b = (xm , Ī£(xm )) , c = .

pm

ā‚Y ā‚X

Since we assumed that a is a unit in Zq , any solution Ī“ to this Artinā“Schreier

equation satisļ¬es Ī“ ā” Ī“m (mod pm ), and thus xm + pm Ī“ ā” x (mod p2m ).

Therefore, if we can solve Artinā“Schreier equations eļ¬ciently, we obtain an

algorithm to ļ¬nd a root of Ī“(X, Ī£(X)) = 0 that converges quadratically. This

strategy is formalized in Algorithm VI.8. Since the precision we compute with

doubles in every step, the complexity of Algorithm VI.8 is the same as the

complexity of the algorithm to solve an Artinā“Schreier equation.

126 VI. POINT COUNTING

Algorithm VI.8: Gen Newton Lift

Polynomial Ī“(X, Y ) ā Zq , element x0 ā Zq satisfying

INPUT:

Ī“(x0 , Ī£(x0 )) ā” 0 (mod p) and precision m.

OUTPUT: Element xm ā Zq , with Ī“(xm , Ī£(xm )) ā” 0 (mod pm ) and

xm ā” x0 (mod p).

1. If m = 1 then

x ā x0 .

2.

3. Else

mā m .

4. 2

x ā Gen Newton Lift(Ī“, x0 , m );

5.

y ā Ī£(x ) (mod pm ).

6.

V ā Ī“(x , y ) (mod pm ).

7.

āx ā ā‚X (x , y ) (mod pm ).

ā‚Ī“

8.

āy ā ā‚Y (x , y ) (mod pm ).

ā‚Ī“

9.

Ī“ ā Artin Schreier Root (āy , āx , V /pm , m ā’ m ).

10.

x ā x + pm Ī“ (mod pm ).

11.

12. Return x.

Lercierā“Lubicz Algorithm: In [218] Lercier and Lubicz described the ļ¬rst

quadratically convergent algorithm to solve Artinā“Schreier equations. Since

the algorithm uses iterated Frobenius substitutions, it is only eļ¬cient if Qq

admits a Gaussian Normal Basis of small type.

Dividing equation (VI.16) by a, we obtain

Ī£(x) ā’ Ī² x ā’ Ī³ = 0 ,

with Ī² = ā’b/a, Ī³ = ā’c/a ā Zq . Let Ī²1 = Ī², Ī³1 = Ī³; then Ī£(x) = Ī²1 x + Ī³1

and deļ¬ne Ī£k (x) = Ī²k x + Ī³k for k = 2, . . . , n by induction. Since Ī£n (x) = x

for all x ā Qq , we conclude that the unique solution to the above equation is

given by Ī³n /(1 ā’ Ī²n ). To compute the Ī²k , Ī³k ā Zq , Lercier and Lubicz use a

simple square and multiply algorithm based on the formula

Ī£k+l (x) = Ī£l (Ī²k x + Ī³k ) = Ī£l (Ī²k )(Ī²l x + Ī³l ) + Ī£l (Ī³k ) .

The above formula immediately leads to Algorithm VI.9, which returns Ī²k

and Ī³k modulo pm for any k.

The complexity of Algorithm VI.9 is determined by Steps 6 and 7, which

require O(1) multiplications and O(1) iterated Frobenius substitutions in

Zq /(pm Zq ). For ļ¬elds with a Gaussian Normal Basis of small type, such a

substitution takes O(n) bit-operations as shown in Section VI.4.1 and a mul-

tiplication takes O((nm)Āµ ) bit-operations. Since the algorithm needs O(log n)

recursive calls, the time and space complexity for ļ¬elds with Gaussian Normal

Basis of small type are O((nm)Āµ log n) and O(nm), respectively.

VI.4. GENERALIZED NEWTON ITERATION 127

Algorithm VI.9: Lercier Lubicz

INPUT: Elements Ī², Ī³ ā Zq , power k, precision m.

OUTPUT: Elements Ī²k , Ī³k ā Zq such that Ī£k (x) ā” Ī²k x + Ī³k (mod pm )

and Ī£(x) ā” Ī²x + Ī³ (mod pm ).

1. If k = 1 then

Ī²k ā Ī² (mod pm ), Ī³k ā Ī³ (mod pm ).

2.

3. Else

kā k.

4. 2

(Ī²k , Ī³k ) ā Lercier Lubicz(Ī², Ī³, k , m).

5.

Ī²k ā Ī²k Ī£k (Ī²k ) (mod pm ).

6.

Ī³k ā Ī³k Ī£k (Ī²k ) + Ī£k (Ī³k ) (mod pm ).

7.

If k ā” 1 (mod 2) then

8.

Ī²k ā Ī² Ī£(Ī²k ) (mod pm ).

9.

Ī³k ā Ī³ Ī£(Ī²k ) + Ī£(Ī³k ) (mod pm ).

10.

11. Return (Ī²k , Ī³k ).

Harleyā™s Algorithm: In an e-mail to the NMBRTHRY list [158], Harley an-

nounced an algorithm to count the number of points on an elliptic curve over

Fpn which runs in O(n2Āµ log n) time and uses O(n2 ) space for ļ¬xed p. Note

that this algorithm is protected by a U.S. patent [159] ļ¬led by Harley. The

key idea is a doubly recursive algorithm to solve the Artinā“Schreier equation

a Ī£(x) + b x + c ā” 0 (mod pm ) , (VI.17)

with a, b, c ā Zq , ordp (a) = 0 and ordp (b) > 0. Let m = m/2 and assume

we have an algorithm which returns a zero xm of (VI.17) to precision m .

Then we can use the same algorithm to ļ¬nd a solution xm to precision m.

Indeed, substituting xm = xm + pm ām in equation (VI.17) and dividing by

pm gives

a Ī£(xm ) + b xm + c

ā” 0 (mod pmā’m ) .

a Ī£(ām ) + b ām + m

p

This shows that ām itself also satisļ¬es an Artinā“Schreier equation of the

form (VI.17) but to a lower precision. Since m ā’ m ā¤ m we can use the

same algorithm to determine ām (mod pmā’m ) and therefore xm .

This immediately leads to a doubly recursive algorithm if we can solve the

base case, i.e., ļ¬nd a solution to (VI.17) modulo p. Since we assumed that

ordp (a) = 0 and ordp (b) > 0, we simply have to solve a xp + c ā” 0 (mod p).

128 VI. POINT COUNTING

Algorithm VI.10: Harley

INPUT: Elements a, b, c ā Zq , a ā ZĆ— , ordp (b) > 0, precision m.

q

OUTPUT: Element x ā Zq such that a Ī£(x) + b x + c ā” (mod pm ).

If m = 1 then

1.

x ā (ā’c/a)1/p (mod p).

2.

3. Else

m ā m ; M =mā’m.

4. 2

x ā Harley(a, b, c, m ).

5.

c ā aĪ£(x p)+bx +c (mod pM ).

6. m

ā ā Harley(a, b, c , M ).

7.

x ā x + pm ā (mod pm ).

8.

Return x.

9.

The pth root in Step 2 of Algorithm VI.10 should not be computed by

naively taking the pnā’1 th power. Instead, let Fq = Fp [Īø]; then

1/p pā’1

nā’1

i k

ai Īø = apk+j Īø Cj (Īø),

i=0 j=0 0ā¤pk+j<n

jpnā’1

j

. This shows that for z ā Fq , we can compute

with Cj (Īø) = (Īø )1/p = Īø

z with p ā’ 1 multiplications over Fq .

1/p

The complexity of Algorithm VI.10 is determined by the recursive calls in

Steps 5 and 7 and the Frobenius substitution in Step 6.

If we assume that Zq is represented as in Section VI.4.1, the Frobenius

substitution in Zq /(pm Zq ) can be computed using O((nm)Āµ ) bit-operations.

If T (m) is the running time of Algorithm VI.10 for precision m, then we have

T (m) ā¤ 2T ( m/2 ) + c(nm)Āµ ,

for some constant c. The above relation implies by induction that the time

complexity of Algorithm VI.10 is O((nm)Āµ log m). The space complexity is

clearly given by O(nm).

VI.5. Norm Computation

The ļ¬nal step in all point counting algorithms is to compute the norm

nā’1

NQq /Qp (Ī±) = i=0 Ī£i (Ī±) of a unit Ī± ā Zq . In this section we give an overview

of the existing algorithms to compute NQq /Qp (Ī±).

VI.5. NORM COMPUTATION 129

VI.5.1. Norm Computation I. Kedlaya [196] suggested a basic square

and multiply approach by computing

i

Ī±i+1 = Ī£2 (Ī±i ) Ī±i for i = 0, . . . , log2 n ,

with Ī±0 = Ī±, and to combine these to recover NQq /Qp (Ī±) = Ī£nā’1 (Ī±) Ā· Ā· Ā· Ī£(Ī±)Ī±.

Let n = l ni 2i , with ni ā {0, 1}, and nl = 1; then we can write

i=0

l

i+1 +2i+2 +Ā·Ā·Ā·+2l n

Ī£2 (Ī±i i ) ,

NQq /Qp (Ī±) =

i=0

where the sum 2i+1 +Ā· Ā· Ā·+2l is deļ¬ned zero for i ā„ l. This formula immediately

leads to Algorithm VI.11.

Algorithm VI.11: Norm I

INPUT: An element Ī± ā Zq with q = pn and a precision m.

OUTPUT: The norm NQq /Qp (Ī±) (mod pm ).

1. i ā n, j ā 0, r ā 1, s ā Ī±.

2. While i > 0 do:

j

If i ā” 1 (mod 2) then r ā Ī£2 (r) s (mod pm ).

3.

j

If i > 1 then s ā Ī£2 (s) s (mod pm ).

4.

j ā j + 1, i ā i/2 .

5.

6. Return r.

This algorithm is particularly attractive for p-adic ļ¬elds with Gaussian

Normal Basis of small type, due to eļ¬cient repeated Frobenius substitutions.

In this case the time complexity is determined by the O(log n) multiplications

in Zq /(pm Zq ), which amount to O((nm)Āµ log n) bit-operations and the space

complexity is O(nm).

VI.5.2. Norm Computation II. Satoh, Skjernaa and Taguchi [287] also

proposed a fast norm computation algorithm based on an analytic method.

Let exp(x) = ā xi /i! and log(x) = ā (ā’1)iā’1 (x ā’ 1)i /i be the p-adic

i=0 i=1

exponential and logarithmic function. An easy calculation shows that exp(x)

converges for ordp (x) > 1/(p ā’ 1) and log(x) converges for ordp (x ā’ 1) > 0.

Assume ļ¬rst that Ī± is close to unity, i.e., ordp (Ī± ā’ 1) > 1/(p ā’ 1); then

NQq /Qp (Ī±) = exp TrQq /Qp (log(Ī±)) ,

since Ī£ is continuous and both series converge. The dominant step in the

above formula is the evaluation of the logarithm. Using a naive algorithm,

this would take O(m) multiplications over Zq /(pm Zq ) or O(nĀµ mĀµ+1 ) bit-

operations.

k

Satoh, Skjernaa and Taguchi solve this problem by noting that Ī±p for

k

k ā N is very close to unity, i.e., ordp (Ī±p ā’ 1) > k + pā’1 . If Ī± ā Zq /(pm Zq ),

1

130 VI. POINT COUNTING

k

then Ī±p is well deļ¬ned in Zq /(pm+k Zq ) and can be computed with O(k)

multiplications in Zq /(pm+k Zq ). Furthermore, note that

k

log(Ī±) ā” pā’k (log(Ī±p ) (mod pm+k )) (mod pm )

k

and that log(ap ) (mod pm+k ) can be computed with O(m/k) multiplications

ā

over Zq /(pm+k Zq ). So if we take k m, then log(Ī±) (mod pm ) can be

computed in O(nĀµ mĀµ+0.5 ) bit-operations.

In characteristic 2, we can use an extra trick. Since ord2 (Ī± ā’ 1) > 1, we

have Ī± ā” 1 (mod 2Ī½ ) for Ī½ ā„ 2. Let z = Ī± ā’ 1 ā 2Ī½ Zq /(2m Zq ) and deļ¬ne

Ī³ = 2+z ā 2Ī½ā’1 Zq /(2mā’1 Zq ). Then Ī± = 1 + z = 1+Ī³ and thus

z

1ā’Ī³

ā

Ī³ 2jā’1

log(Ī±) = log(1 + z) = 2 .

2j ā’ 1

j=1

Note that all the denominators in the above formula are odd. Reducing this

equation modulo 2m therefore leads to

Ī³ 2jā’1

log(Ī±) ā” log(1 + z) ā” 2 (mod 2m ) .

2j ā’ 1

1ā¤(Ī½ā’1)(2jā’1)<mā’1

To compute the trace TrQq /Qp , Satoh, Skjernaa and Taguchi proceed as

nā’1

follows: deļ¬ne Ī±i ā Zp by log(Ī±) ā” i=0 Ī±i ti (mod pm ); then

nā’1

TrQq /Qp (log(Ī±)) ā” Ī±i TrQq /Qp (ti ) (mod pm )

i=0

and each TrQq /Qp (ti ) for i = 0, . . . , n ā’ 1 can be precomputed using Newtonā™s

formula:

iā’1

TrQq /Qp (tiā’j )fnā’j + ifnā’i ā” 0 (mod pm ) ,

i

TrQq /Qp (t ) +

j=1

with f (t) = n fi ti the deļ¬ning polynomial of Qq ā¼ Qp [t]/(f (t)).

=

i=0

Let ta = TrQq /Qp (log(a)) ā Zp . Then ordp (ta ) > 1/(p ā’ 1), so ta ā pZp for

p ā„ 3 and ta ā 4Z2 for p = 2. If we precompute exp(p) (mod pm ) for p ā„ 3

or exp(4) (mod 2m ) for p = 2, then

exp(ta ) ā” exp(p)ta /p (mod pm ) for p ā„ 3,

exp(ta ) ā” exp(4)ta /4 (mod 2m ) for p = 2,

and we can use a simple square and multiply algorithm to compute exp(ta ).

Therefore, if Ī± is close to unity, then NQq /Qp (Ī±) (mod pm ) can be computed

in O(nĀµ mĀµ+0.5 ) bit-operations and O(nm) space. Algorithm VI.12 computes

the norm of an element in 1 + 2Ī½ Zq , with q = 2n , assuming that exp(4) and

TrQq /Qp (ti ) for i = 0, . . . , n ā’ 1 are precomputed.

VI.5. NORM COMPUTATION 131

Algorithm VI.12: Norm II

INPUT: An element Ī± ā 1 + 2Ī½ Zq with Ī½ ā„ 2 and a precision m.

OUTPUT: The norm NQq /Qp (Ī±) (mod 2m ).

ā

1. s ā m/2 .

s

2. z ā Ī±2 ā’ 1 (mod 2m+s ).

3. w ā log(1 + z) (mod 2m+s ).

4. w ā 2ā’s w (mod 2m ).

5. u ā 2ā’Ī½ TrQq /Qp (w) (mod 2mā’Ī½ ).

6. Return exp(4)u (mod 2m ).

Equations (VI.5) and (VI.13) show that for p = 2 we can indeed assume

that ordp (Ī± ā’ 1) > 1/(p ā’ 1); however, for p ā„ 3 we also need to consider the

more general situation where Ī± ā Zq is not close to unity. Let Ī± ā Fq denote

the reduction of Ī± modulo p and Ī±t ā Zq the TeichmĀØller lift of Ī±, i.e., the

u

unique (q ā’ 1)th root of unity which reduces to Ī±. Consider the equality

ā’1

NQq /Qp (Ī±) = NQq /Qp (Ī±t )NQq /Qp (Ī±t Ī±) ;

ā’1

then ordp (Ī±t Ī± ā’ 1) ā„ 1 > 1/(p ā’ 1) since p is odd. Furthermore, note that

NQq /Qp (Ī±t ) is equal to the TeichmĀØller lift of NFq /Fp (Ī±). Satoh [286] pointed

u

out that the TeichmĀØller lift Ī±t can be computed eļ¬ciently as the solution of

u

and X ā” Ī± (mod p) .

Ī£(X) = X p

Note that we can apply the same algorithms as in Section VI.4.

VI.5.3. Norm Computation III. In an e-mail to the NMBRTHRY list [158],

Harley suggested an asymptotically fast norm computation algorithm based

on a formula from number theory that expresses the norm as a resultant.

The resultant itself can be computed using an adaptation of Moenckā™s fast

extended GCD algorithm [252].

Let Zq ā¼ Zp [t]/(f (t)) with f (t) ā Zp [t] a monic irreducible polynomial of

=

degree n. Let Īø ā Zq be a root of f (t); then f (t) splits completely over Zq

nā’1 nā’1

as f (t) = i=0 (x ā’ Ī£i (Īø)). For Ī± = i=0 Ī±i Īøi ā Zq , deļ¬ne the polynomial

nā’1

A(x) = i=0 Ī±i xi ā Zp [x]. By deļ¬nition of the norm and the resultant we

have

nā’1 nā’1

i

A(Ī£i (Īø)) = Res(f (x), A(x)) .

NQq /Qp (Ī±) = Ī£ (Ī±) =

i=0 i=0

The resultant Res(f (x), A(x)) can be computed in softly linear time using a

variant of Moenckā™s fast extended GCD algorithm [252] as described in detail

in Chapter 11 of the book by von zur Gathen and Gerhard [141]. The result

is an algorithm to compute NQq /Qp (Ī±) (mod pm ) in time O((nm)Āµ log n).

132 VI. POINT COUNTING

VI.6. Concluding Remarks

The p-adic approach, introduced by Satoh at the end of 1999, has caused a

fundamental dichotomy in point counting algorithms that essentially depends

on the characteristic of the ļ¬nite ļ¬eld. For elliptic curves over ļ¬nite ļ¬elds Fpn

with p moderately large, the l-adic approach, i.e., the SEA algorithm, remains

the only eļ¬cient solution and requires O(n2Āµ+2 log2Āµ+2 p) time. This stands in

stark contrast with the p-adic approach which only requires O(n2Āµ log n) time

for ļ¬xed p. The O-constant grows exponentially in p though, which implies

that p-adic algorithms are only eļ¬cient for small p. The huge eļ¬ciency gap

between the l-adic and p-adic approaches in small characteristic is illustrated

by the following examples: for an elliptic curve over F2163 (resp. F2509 ) the

SEA algorithm runs about two (resp. three) orders of magnitude slower than

the p-adic algorithms.

In this chapter we have only considered the advances in elliptic curve

point counting algorithms. To provide a starting point, we brieļ¬‚y review the

most important developments for higher genus curves, especially hyperelliptic

curves. The p-adic approach currently comes in two ļ¬‚avours:

Canonical Lift / AGM. This method computes the Serre-Tate canonical

lift of an ordinary abelian variety and recovers the zeta function from the

action of Frobenius on a basis for the holomorphic diļ¬erentials. The main

references are as follows: Mestre [247, 248] generalized the AGM algorithm

to ordinary hyperelliptic curves over F2n . An optimized version by Lercier

and Lubicz [219] runs in time O(2g n2Āµ ) for an ordinary genus g hyperelliptic

curve over F2n . Ritzenthaler [282] extended Mestreā™s algorithm to ordinary

genus 3 non-hyperelliptic curves.

p-adic Cohomology. These methods compute the action of a Frobenius

operator on p-adic cohomology groups. Diļ¬erent cohomology theories give

rise to diļ¬erent algorithms: Lauder and Wan [213] used p-adic analysis a la

`

Dwork [113] to show that the zeta function of any algebraic variety can be

computed in polynomial time. Lauder and Wan [214, 215] also devized an

algorithm based on Dwork cohomology to compute the zeta function of an

Artinā“Schreier curve in time O(g 2Āµ+3 n3Āµ ). Kedlaya [196] gave algorithms to

compute with Monskyā“Washnitzer cohomology and obtained an O(g 2Āµ+2 n3Āµ )

time algorithm to compute the zeta function of a hyperelliptic curve in odd

characteristic. Kedlayaā™s algorithm was subsequently extended as follows: to

characteristic 2 by Denef and Vercauteren [96, 97]; to superelliptic curves

by Gaudry and GĀØrel [144] and, ļ¬nally, to Cab curves by Denef and Ver-

u

cauteren [98]. A nice introduction to these p-adic cohomology theories can

be found in the PhD thesis of Gerkmann [148].

CHAPTER VII

Hyperelliptic Curves and the HCDLP

P. Gaudry

In 1989, hyperelliptic curves were proposed by Koblitz for use in cryptog-

raphy as an alternative to elliptic curves. Recent improvements on the algo-

rithms for computing the group law tend to prove that indeed cryptosystems

based on hyperelliptic curves can be competitive. However, another interest

of the study of hyperelliptic curves ā“ and also of more general curves ā“ is

that the process called the Weil descent relates some instances of the ECDLP

to instances of the hyperelliptic discrete logarithm problem (HCDLP). The

purpose of this chapter is to give eļ¬cient algorithms for hyperelliptic curves

that are useful for the next chapter.

VII.1. Generalities on Hyperelliptic Curves

In this section, we recall brieļ¬‚y the basic deļ¬nitions and properties of

hyperelliptic curves and their Jacobians. A complete elementary introduction

to hyperelliptic curves that includes all the proofs can be found in [243]. For

a more general treatment see [226].

VII.1.1. Equation of a Hyperelliptic Curve. Hyperelliptic curves can

be deļ¬ned in an abstract manner as follows: a hyperelliptic curve is a curve

which is a covering of degree 2 of the projective line. In a more concrete

setting, a hyperelliptic curve C of genus g is a curve with an equation of the

form

C : Y 2 + H(X)Y = F (X),

where F and H are polynomials with deg F = 2g +1 and deg H ā¤ g. Further-

more we ask the curve to be non-singular as an aļ¬ne curve, which translates

into the non-simultaneous vanishing of the partial derivatives of the equation.

This choice of equation is not always possible. In general, the polynomial

H(X) has degree at most g + 1 and F (X) has degree at most 2g + 2. The

condition for being able to reduce the degrees to g and 2g+1 is to have a ratio-

nal WeierstraĆ point, which is equivalent to saying that H(X) has a rational

root if the characteristic is two, or that 4F (X) + H(X)2 has a rational root

if the characteristic is odd. For this, we apply an invertible transformation

X = aX+b which sends this rational root to inļ¬nity. If the polynomials F

cX+d

and H have degrees 2g + 1 and at most g, we call the corresponding equation

an imaginary model of the curve. Otherwise, the equation is said to be a real

133

134 VII. HYPERELLIPTIC CURVES

model of the curve. In the following, we always consider an imaginary model.

Using the real model introduces some extra technicalities, since there are now

zero or two points at inļ¬nity.

In the case where the characteristic is odd, it is also possible to set H(X) =

0 without lost of generality. Indeed, we can apply the transformation Y =

Y + H(X)/2, which kills the crossed term. The division by two, prevents one

from doing the same in characteristic two.

A hyperelliptic curve comes with a natural involution called the hyperel-

liptic involution, which we denote by Ī¹. If P = (x, y) is a point of C, then Ī¹(P )

is the point (x, ā’y ā’ h(x)), which is also on C. The points which are ļ¬xed

under the action of Ī¹ are called the ramiļ¬cation points and they are actually

the WeierstraĆ points of C.

Just like for elliptic curves, we manipulate the curve via an aļ¬ne equation,

but we virtually always work with the projective model. For that we denote

by ā the point of C at inļ¬nity, which is unique in the imaginary model. We

call aļ¬ne points the points of C that are not ā.

VII.1.2. The Jacobian of a Hyperelliptic Curve. Let C be a hyperel-

liptic curve. Similarly to elliptic curves, C comes with a group associated to

it. However, the group law is no longer deļ¬ned over the curve itself but on

another object called the Jacobian of C and denoted by JC . First, we need to

consider the group of divisors on C: these are the ļ¬nite formal sums of points

of C over Z:

ni ā Z,

Div(C) = D= ni Pi ; ni = 0 for almost all i .

Pi āC

The degree of a divisor D = ni Pi is the sum of the coeļ¬cients deg D = ni

and the set of points Pi for which ni is diļ¬erent from 0 is called the support

of D. The set of divisors of degree 0 forms a subgroup of Div(C).

The function ļ¬eld K(C) of C is the set of rational functions on C. Formally,

K(C) is deļ¬ned as the following ļ¬eld of fractions:

K(C) = Frac K[X, Y ]/(Y 2 + H(X)Y ā’ F (X)) .

To each non-zero function Ļ• in K(C) we can associate a divisor div(Ļ•) =

Ī½Pi (Ļ•)Pi , where Ī½Pi (Ļ•) is an integer deļ¬ned as follows: if Ļ• has a zero at

Pi , then Ī½Pi (Ļ•) is the multiplicity of this zero; if Ļ• has a pole at Pi , Ī½Pi (Ļ•)

is the opposite of the order of this pole. Otherwise, Ī½Pi (Ļ•) = 0. A non-

zero function has only ļ¬nitely many zeroes and poles, therefore almost all

the Ī½Pi (Ļ•) are zero and the divisor div(Ļ•) is well deļ¬ned. Such a divisor is

called a principal divisor. Furthermore, it can be shown that the degree of a

principal divisor is always zero and that the set of the principal divisors is a

subgroup of the group of degree 0 divisors. We deļ¬ne the Jacobian JC of C

to be the quotient group of the degree 0 divisors by the principal divisors. A

VII.1. GENERALITIES ON HYPERELLIPTIC CURVES 135

consequence of Riemannā“Roch theorem is that every element of the Jacobian

can be represented by a divisor of the form

D = P1 + P2 + Ā· Ā· Ā· + Pr ā’ rā,

where, for all i, the point Pi is an aļ¬ne point of C and the integer r is less

than or equal to the genus. If furthermore we ask that Pi be diļ¬erent from

Ī¹(Pj ) for all pairs i = j, then this representation is unique. Such a divisor,

with Pi = Ī¹(Pj ) for i = j, is called reduced.

Until now, we did not pay attention to the ļ¬eld of deļ¬nition of the objects

we considered. Usually, the good way to say that an object is deļ¬ned over a

ļ¬eld K is to say that it is invariant under the action of any element of the

Galois group Gal(K/K). In the case of divisors, this deļ¬nition is indeed the

one we want in the sense that it makes it possible to deļ¬ne the elements of

the Jacobian deļ¬ned over K as a subgroup of the Jacobian deļ¬ned over K.

However, this implies in particular that a divisor deļ¬ned over K can include

some points Pi in its support which are not deļ¬ned over K. In that case all

the Galois conjugates of Pi appear within D with the same multiplicities as Pi .

From the practical point of view, it is annoying to have to work in extensions

of the base ļ¬elds, even when working with elements which are deļ¬ned over

the base ļ¬eld. The fact that a point always comes with its conjugates leads

to the following representation of reduced divisors.

Proposition VII.1. Let C be a hyperelliptic curve of genus g deļ¬ned over

a ļ¬eld K with equation Y 2 + H(X)Y = F (X). Then the elements of the

Jacobian of C that are deļ¬ned over K are in one-to-one correspondance

with the pairs of polynomials (a(X), b(X)) with coeļ¬cients in K, such that

deg(b) < deg(a) ā¤ g, the polynomial a is monic, and a divides b2 + bH ā’ F .

If a and b are two polynomials that satisfy these conditions, we denote by

div(a, b) the corresponding element of JC .

Let D = 1ā¤iā¤r Pi ā’ rā be a reduced divisor and denote by (xi , yi ) the

coordinates of Pi . Then the a-polynomial is given by a(X) = 1ā¤iā¤r (X ā’

xi ) and the b-polynomial is the only polynomial of degree less than r such

that b(xi ) = yi . In the case where the same point occurs several times, the

multiplicity must be taken into account to properly deļ¬ne b(X). From now

on, when we talk about a reduced divisor, we are actually thinking about this

representation with the two polynomials a and b.

The weight of a reduced divisor D = div(a, b) is the degree of the polyno-

mial a. We denote it by wt(D), and in the representation

Pi ā’ rā

D=

1ā¤iā¤r

the weight is given by r. By deļ¬nition, the weight of a reduced divisor is an

integer in [0, g]. The neutral element of JC is the unique weight zero divisor

div(1, 0).

136 VII. HYPERELLIPTIC CURVES

VII.1.3. Hyperelliptic Curves Deļ¬ned over a Finite Field. When the

base ļ¬eld K = Fq is ļ¬nite of cardinality q, the set of (a, b) pairs of polynomials

that we can form is ļ¬nite, and therefore the Jacobian is a ļ¬nite abelian

group. Similarly to the case of elliptic curves, the Frobenius endomorphism

Ļ• obtained by extending the map x ā’ xq to the Jacobian plays a central role

in the question of evaluating the cardinalities. Due to the theorem of Weil,

the characteristic polynomial of Ļ• is of the form

P (T ) = T 2g + a1 T 2gā’1 + Ā· Ā· Ā· + ag T g + qagā’1 T gā’1 + Ā· Ā· Ā· + q gā’1 a1 T + q g ,

where ai is an integer for all i. Furthermore, the cardinality of the Jacobian

ā

is given by P (1), and all the roots of P (T ) have absolute value q. There

is also a link between the ai and the cardinality of C over the extensions Fqi

for 1 ā¤ i ā¤ g. We skip the details that can be found for instance in [226] or

[319] and we give the bounds that we obtain:

Theorem VII.2. Let C be a hyperelliptic curve of genus g deļ¬ned over a

ļ¬nite ļ¬eld Fq with q elements. Then we have

ā ā

( q ā’ 1)2g ā¤ #JC (Fq ) ā¤ ( q + 1)2g

and ā

|#C(Fq ) ā’ (q + 1)| ā¤ 2g q.

This theorem shows that indeed the Jacobians of hyperelliptic curves de-

ļ¬ned over a ļ¬nite ļ¬eld are potentially good candidates for cryptosystems

based on the discrete logarithm problem in an abelian group. Indeed, to get

a group of size N , it is necessary to have a curve of genus g deļ¬ned over a

ļ¬nite ļ¬eld with about N 1/g elements. If the group law is not too complicated

compared to elliptic curves, then the higher number of operations can be

somewhat balanced by the fact that they take place in a ļ¬eld whose bit-size

is g times smaller.

VII.2. Algorithms for Computing the Group Law

In this section we consider the problem of computing eļ¬ectively the group

law in the Jacobian of a hyperelliptic curve. Having a good algorithm for

the group law is crucial when one wants to design a cryptosystem based on

hyperelliptic curves.

In all this section, we consider a hyperelliptic curve C of genus g deļ¬ned

over a ļ¬eld K (not necessarily ļ¬nite) by the equation Y 2 + H(X)Y = F (X).

VII.2.1. Algorithms for Small Genus. The ļ¬rst algorithm to compute

in JC is due to Cantor [59] and is a translation of Gaussā™s composition and

reduction algorithms for quadratic forms. Here we give another algorithm

described in [275] and which goes back to Lagrange. In general this variant

saves many operations in the reduction phase. However, for genus 2 this is

exactly Cantorā™s algorithm.

VII.2. ALGORITHMS FOR COMPUTING THE GROUP LAW 137

Algorithm VII.1: Lagrangeā™s Algorithm for Hyperelliptic Group Law

INPUT: Two reduced divisors D1 = div(a1 , b1 ) and D2 = div(a2 , b2 ).

OUTPUT: The reduced divisor D3 = div(a3 , b3 ) = D1 + D2 .

Composition.

1. Use two extended GCD computations, to get

d = gcd(a1 , a2 , b1 + b2 + H) = u1 a1 + u2 a2 + u3 (b1 + b2 + H).

2. a3 ā a1 a2 /d2

3. b3 ā b1 + (u1 a1 (b2 ā’ b1 ) + u3 (F ā’ b2 ā’ b1 H))/d (mod a3 )

1

4. If deg(a3 ) ā¤ g then return div(a3 , b3 ).

Reduction.

5. a3 ā a3 , Ė3 ā b3 .

Ė b

6. a3 ā (F ā’ b3 H ā’ b2 )/a3 .

3

7. q, b3 ā Quotrem(ā’b3 ā’ H, a3 ).

8. While deg(a3 ) > g do

t ā a3 + q(b3 ā’ Ė3 ).

Ė b

9.

Ė3 ā b3 , a3 = a3 , a3 ā t.

b Ė

10.

q, b3 ā Quotrem(ā’b3 ā’ H, a3 ).

11.

12. Return div(a3 , b3 ).

The important case where D1 = D2 must be treated separately; in-

deed some operations can be saved. For instance, we readily know that

gcd(a1 , a2 ) = a1 .

To analyze the run time of this algorithm, we evaluate how the number of

base ļ¬eld operations grows with the genus. The composition steps require one

to perform the Euclidean algorithm on polynomials of degree at most g and a

constant number of basic operations with polynomials of degree at most 2g.

This can be done at a cost of O(g 2 ) operations in the base ļ¬eld. The cost of

the reduction phase is less obvious. If the base ļ¬eld is large enough compared

to the genus, each division is expected to generate a remainder of degree

exactly one less than the modulus. The quotients q involved in the algorithm

are then always of degree one and one quickly see that the cost is again O(g 2 )

operations in the base ļ¬eld for the whole reduction process. It can be shown

[317] that even if the degrees of the quotients behave in an erratic way, the

heavy cost of a step where the quotient is large is amortized by the fact that

we have less steps to do and ļ¬nally the overall cost of reduction is always

O(g 2 ) base ļ¬eld operations.

VII.2.2. Asymptotically Fast Algorithms. In the late eighties, Shanks

[298] proposed a variant of Gaussā™s algorithm for composing positive deļ¬nite

quadratic forms. The main feature of this algorithm, called NUCOMP, is that

most of the operations can be done with integers half the size of the integers

that are required in the classical algorithm. This makes a great diļ¬erence if

138 VII. HYPERELLIPTIC CURVES

we are in the range where applying NUCOMP allows us to use machine word

sized integers instead of multiprecision integers.

In the classical algorithm, just like in Lagrangeā™s algorithm, the compo-

sition step produces large elements that are subsequently reduced. The idea

of the NUCOMP algorithm is to more or less merge the composition and the

reduction phases, in order to avoid those large elements. A key ingredient is

the partial GCD algorithm (also known as the Half GCD algorithm, which is

related to the continued fractions algorithm). This is essentially a Euclidean

algorithm that we stop in the middle. In NUCOMP we need a partial ex-

tended version of it, where only one of the coeļ¬cients of the representation

is computed.

In [182], it is shown that the NUCOMP algorithm can be adapted to

the case of the group law of a hyperelliptic curve. It suļ¬ces more or less to

replace integers by polynomials. We recall the Partial Extended Half GCD

algorithm here for completeness:

Algorithm VII.2: Partial Extended Half GCD (PartExtHGCD)

Two polynomials A and B of degree ā¤ n and a bound

INPUT:

m

OUTPUT: Four polynomials a, b, Ī±, Ī² such that

a ā” Ī±A (mod B) and b ā” Ī²A (mod B)

and deg a and deg b are bounded by m.

An integer z recording the number of operations

Ī± ā 1, Ī² ā 0, z ā 0, a ā A mod B, b ā B.

1.

While deg b > m and a = 0 do

2.

q, t ā Quotrem(b, a).

3.

b ā a, a ā t, t ā Ī² ā’ qĪ±, Ī² ā Ī±. Ī± ā t.

4.

z ā z + 1.

5.

If a = 0 then

6.

Return (A mod B, B, 1, 0, 0).

7.

If z is odd then

8.

Ī² ā ā’Ī², b ā ā’b.

9.

Return (a, b, Ī±, Ī², z).

10.

We are now ready to give the NUCOMP algorithm as described in [182].

We restrict ourselves to the case of odd characteristic where the equation of

the curve can be put in the form Y 2 = F (X). The minor modiļ¬cations to

handle the characteristic 2 case can be found in [182].

VII.2. ALGORITHMS FOR COMPUTING THE GROUP LAW 139

Algorithm VII.3: NUCOMP Algorithm

INPUT: Two reduced divisors D1 = div(u1 , v1 ) and D2 = div(u2 , v2 ).

OUTPUT: The reduced divisor D3 = D1 + D2 .

1. w1 ā (v1 2 ā’ F )/u1 , w2 ā (v2 2 ā’ F )/u2 .

2. If deg w2 < deg w1 then exchange D1 and D2 .

3. s ā v1 + v2 , m ā v2 ā’ v1 .

4. By Extended GCD, compute G0 , b, c such that

G0 = gcd(u2 , u1 ) = bu2 + cu1 .

5. If G0 |s then

By Extended GCD, compute G, x, y such that

6.

G = gcd(G0 , s) = xG0 + ys.

H ā G0 /G, By ā u1 /G, Cy ā u2 /G, Dy ā s/G.

7.

l ā y(bw1 + cw2 ) mod H, Bx ā b(m/H) + l(By /H).

8.

9. Else

G ā G0 , Ax ā G, Bx ā mb.

10.

By ā u1 /G, Cy ā u2 /G, Dy ā s/G.

11.

12. (bx , by , x, y, z) ā PartExtHGCD(Bx mod By , By , (g + 2)/2).

13. ax ā Gx, ay ā Gy.

14. If z = 0 then

u3 ā by Cy , Q1 ā Cy bx , v3 ā v2 ā’ Q1 .

15.

16. Else

cx ā (Cy bx ā’ mx)/By , Q1 ā by cx , Q2 ā Q1 + m.

17.

dx ā (Dy bx ā’ w2 x)/By , Q3 ā ydx , Q4 ā Q3 + Dy .

18.

dy ā Q4 /x, cy ā Q2 /bx .

19.

u3 ā by cy ā’ ay dy , v3 ā (G(Q3 + Q4 ) ā’ Q1 ā’ Q2 )/2 mod u3 .

20.

21. While deg(u3 ) > g do

u3 ā (v3 2 ā’ F )/u3 , v3 ā ā’v3 mod u3 .

22.

23. Return div(u3 , v3 ).

In the NUCOMP algorithm, if the input divisors are distinct, then it is

likely that G0 = 1, and then Steps 5ā“8 are skipped. On the other hand, in the

case of doubling a divisor, these steps are always performed whereas Step 4

is trivial. Hence it is interesting to have another version of NUCOMP that is

speciļ¬c to doubling. Shanks had such an algorithm which he called NUDPL

and it has been extended to our case in [182].

Steps 21ā“22 are there to adjust the result in the case where the divisor

div(u3 , v3 ) is not completely reduced at the end of the main part of the al-

gorithm. In general, at most one or two loops are required. Therefore the

NUCOMP algorithm requires a constant number of Euclidean algorithms on

polynomials of degree O(g) and a constant number of basic operations with

polynomials of degree O(g). Thus the complexity is O(g 2 ), again, just as for

140 VII. HYPERELLIPTIC CURVES

Lagrangeā™s algorithm. However, as we said before, the NUCOMP algorithm

works with polynomials of smaller degree (in practice rarely larger than g).

Up to now we have considered naive algorithms for polynomial multipli-

cations and GCDs. If we used an FFT-based algorithm, two polynomials

of degree O(n) can be multiplied in time O(n log n log log n) operations in

the base ļ¬eld. It is then possible to design a recursive GCD algorithm that

requires O(n log2 n log log n) operations [141]. Furthermore, the same algo-

rithm can also be used to computed the Partial Extended Half GCD in the

same run time (this is actually a building block for the fast recursive GCD

computation). Plugging those algorithms in NUCOMP and NUDUPL yields

an overall complexity of O(g log2 g log log g) which is asymptotically faster

than the O(g 2 ) complexity of Lagrangeā™s algorithm.

VII.2.3. Which Algorithm in Practice. When the genus is large, the

best algorithm is not the same as when the genus is small. In particular,

for very small genus, it is possible to write down the explicit formulae cor-

responding to the Cantor or Lagrange algorithm. Then some operations can

be saved by computing only the coeļ¬cients that are necessary. For instance,

the quotient of two polynomials costs less if it is known in advance that the

division is exact (see [258]). For genus 2 and 3, such optimized formulae have

been worked out by many contributors. We refer to [212] for a survey of the

current situation in genus 2, and to [277] for genus 3. In the case of genus 2,

the best formulae for adding (resp. doubling) have a cost of 25 multiplications

and 1 inversion (resp. 27 multiplications and 1 inversion).

When the genus is larger than 3, explicit formulae are not available and

the best algorithm is then Lagrangeā™s algorithm.

According to [182], the genus for which the NUCOMP and NUDUPL

algorithms are faster than Lagrangeā™s algorithm is around 10. This value de-

pends highly on the implementation of those two algorithms, on the particular

type of processor that is used, and on the size of the base ļ¬eld.

It is likely that the asymptotically fast algorithm for computing GCD will

never be useful for practical applications: the point at which it becomes faster

than the classical Euclidean algorithm is around genus 1000.

We ļ¬nally mention that an analogue of projective coordinates for elliptic

curves has been designed. This allows one to perform the group law without

any inversion, which can be crucial in some constrained environments such as

smart cards. In addition, in the spirit of [85], diļ¬erent types of coordinates

can be mixed to have a faster exponentiation algorithm [212].

VII.3. Classical Algorithms for HCDLP

The hyperelliptic curve discrete logarithm problem (HCDLP) is a gener-

alization of the elliptic curve discrete logarithm problem (ECDLP). As such

it comes as no surprise that every known attack on the ECDLP, for some

VII.3. CLASSICAL ALGORITHMS FOR HCDLP 141

particular instances, can be extended to an attack on the HCDLP, for the cor-

responding weak instances. We review all these attacks that we call classical,

in opposition to the index-calculus attacks which are speciļ¬c to hyperelliptic

curves. We assume that the reader is familiar with all attacks on ECDLP

that are listed in Chapter V of [ECC].

VII.3.1. The Simpliļ¬cation of Pohlig and Hellman. Let G denote an

abelian group of known order N whose factorization is also known. The algo-

rithm of Pohlig and Hellman reduces the computation of a discrete logarithm

in G to O(log N ) computations of discrete logarithms in subgroups of G of

order p where p is a prime factor of N . The method works for a generic group

and in particular can be applied to HCDLP. It means that solving a discrete

logarithm problem in the Jacobian of a hyperelliptic curve is not harder than

the discrete logarithm problem in the largest prime order subgroup of the

Jacobian.

VII.3.2. The MOV and Freyā“RĀØck Attacks. Let C be a hyperelliptic

u

curve of genus g over the ļ¬nite ļ¬eld Fq and let us consider a DL problem in a

subgroup of JC of prime order m. Denote by k the smallest integer such that

m divides q k ā’ 1, i.e., k is the order of q modulo m. Then there is an eļ¬cient

algorithm that converts the initial discrete logarithm problem into a discrete

logarithm problem in the multiplicative group of the ļ¬nite ļ¬eld Fqk . If k

is small enough, the existence of a subexponential algorithm for computing

discrete logarithms in ļ¬nite ļ¬elds can lead to an eļ¬cient way of solving the

initial problem in JC .

For more details on the method to transport the HCDLP to a classical

DLP we refer to Chapter IX. Indeed the key ingredients are the Weil and

Tate pairings which can, as we shall see in Chapter X, also be used in a

constructive way.

VII.3.3. The RĀØck Attack. The previous attack works only if the order m

u

of the subgroup is prime to the ļ¬eld characteristic (otherwise the parameter k

is not even deļ¬ned). In the case where m is precisely the ļ¬eld characteristic,

RĀØck [283] has extended the so-called anomalous curves attack (see [ECC,

u

Section V.3]) to hyperelliptic curves. He proved

Theorem VII.3. Let C be a hyperelliptic curve of genus g deļ¬ned over a

ļ¬nite ļ¬eld of characteristic p. The discrete logarithm problem in a subgroup

of JC of order p can be solved with at most O(log p) operations in the base

ļ¬eld and in the Jacobian.

In fact RĀØck proved his result in a more general setting where the curve

u

does not need to be hyperelliptic.

The principle is the following: let D be a degree 0 divisor whose class

in JC is of order p. Then, by deļ¬nition, there exists a function Ļ• such that

div(Ļ•) = p Ā· D. The holomorphic diļ¬erential dĻ•/Ļ• can be expressed using a

142 VII. HYPERELLIPTIC CURVES

local parameter t at the point at inļ¬nity ā: we can write dĻ•/Ļ• = ā‚Ļ•/ā‚t

dt.

Ļ•

We can then deļ¬ne the power series

ā

ā‚Ļ•/ā‚t

ai ti ,

=

Ļ• i=0

where the ai are scalars. The key point is that the function that maps the

class of D to the 2g ā’ 1 ļ¬rst coeļ¬cients (a0 , . . . , a2gā’2 ) is an homomorphism

of additive groups. The discrete logarithm problem is then transported into

a group where it is easily solvable.

The explicit computation of this homomorphism can be done eļ¬ciently

using a technique similar to the method used for computing the Weil pairing.

VII.3.4. Baby Step/Giant Step and Random Walks. The generic al-

gorithms, such as the Baby Step/Giant Step algorithm, as described in Sec-

tion V.4 of [ECC], can be applied to compute discrete logarithms in any

abelian group. In particular, it is straightforward to obtain an algorithm for

solving the HCDLP in a subgroup of order N of the Jacobian of a hyperelliptic

ā

curve in at most O( N ) group operations.

Also, the parallelized Pollard lambda method of van Oorschot and Wiener

ā

[266] directly applies to HCDLP, providing a O( N ) algorithm with a speed-

up by a factor of M if M machines are used and requiring a small amount of

memory.

Furthermore the presence of equivalence classes arising from the hyperel-

liptic involution as well as any other quickly computable automorphism can

be used just like for elliptic curves [111], [138].

VII.4. Smooth Divisors

In 1994, Adleman, DeMarrais and Huang [4] discovered a subexponential

algorithm for computing discrete logarithms in the Jacobian of hyperelliptic

curves of ālarge genus.ā Compared with points of elliptic curves, divisors in

the Jacobian of a hyperelliptic curve carry more structure. Indeed, a reduced

divisor is represented using two polynomials (a, b), and the factorization of

a as a polynomial in K[X] is somewhat compatible with the Jacobian group

law. This fact is the key stone for deļ¬ning a smoothness notion and then the

Adleman, DeMarrais and Huang index-calculus algorithm.

Lemma VII.4. Let C be a hyperelliptic curve over a ļ¬nite ļ¬eld Fq . Let D =

div(a, b) be a reduced divisor in JC (Fq ). Write a(X) = ai (X) as the factor-

ization in irreducible factors of a(X) in Fq [X]. Let bi (X) = b(X) (mod ai (X)).

Then Di = div(ai , bi ) is a reduced divisor and D = Di in JC .

With this result, a reduced divisor can be rewritten as the sum of reduced

divisors of smaller size, in the sense that the degree of their a-polynomial is

VII.4. SMOOTH DIVISORS 143

smaller. More precisely, we have wt(D) = wt(Di ). The divisors that can

not be rewritten are those for which the a-polynomial is irreducible.

Definition VII.5. A reduced divisor D = div(a, b) is said to be prime if the

polynomial a is irreducible.

Definition VII.6. Let B be an integer. A reduced divisor D = div(a, b) is

said to be B-smooth if the polynomial a is the product of polynomials of degree

at most B.

This deļ¬nition is very analogous to the smoothness notion that is used to

design index-calculus algorithms in the multiplicative group of a ļ¬nite ļ¬eld.

The next step is to study the proportion of smooth divisors.

We start with the easy case where B = 1, which will be of great impor-

tance in applications.

Proposition VII.7. Let C be a hyperelliptic curve of genus g over the ļ¬nite

ļ¬eld Fq . Then the proportion of 1-smooth reduced divisors in JC (Fq ) tends to

1/g! when g is ļ¬xed and q tends to inļ¬nity.

The proof relies on Theorem VII.2. The 1-smooth reduced divisors are

formed using at most g points of C(Fq ), and the order does not matter. There-

fore there are about q g /g! reduced divisors that are 1-smooth. This number

should be compared with #JC (Fq ) ā q g . The error terms given by the precise

bounds are an order of magnitude smaller, so that we obtain the announced

proportion.

If we do not ļ¬x B = 1, intuitively we expect the same kind of behaviour as

for polynomials over ļ¬nite ļ¬elds for which taking a subexponential smooth-

ness bound yields a subexponential proportion of smooth elements. This

heuristic has been made precise by Enge and Stein [118]. We recall the def-

inition of the so-called subexponential function parametrized by a real Ī± in

[0, 1] and a constant c:

LN (Ī±, c) = exp c(log N )Ī± (log log N )1ā’Ī± .

Theorem VII.8 (Enge and Stein). Let C be a hyperelliptic curve of genus

g deļ¬ned over the ļ¬nite ļ¬eld Fq . Let Ļ be a positive constant and let B =

logq Lqg ( 1 , Ļ) be a smoothness bound. Then the number of B-smooth reduced

2

divisors in JC (Fq ) is at least

1 1

q g Ā· Lq g , ā’ + o(1) .

2 2Ļ

It is worth noting that B is deļ¬ned as the ceiling of the usual expres-

sion. Indeed, as a degree, it has to be an integer. This looks inoļ¬ensive

but is actually the heart of the main problem that occurs below, namely the

subexponential behaviour can be observed only if the genus is large enough

compared to the size of the base ļ¬eld. If the genus is too small, the smooth-

ness bound would be smaller than 1, and we are back to the situation of the

previous proposition.

144 VII. HYPERELLIPTIC CURVES

VII.5. Index-Calculus Algorithm for Hyperelliptic Curves

The original algorithm by Adleman, DeMarrais and Huang is highly heuris-

tic and more theoretical than of practical use. Their idea was to work with

functions on the curve. This idea was pushed further by Flassenberg and

Paulus [121] who were able to design a sieving technique, and they did the

ļ¬rst implementation of this kind of algorithm. Later on, MĀØller, Stein and

u

Thiel [257] and Enge [116] managed to get rid of the heuristics, by working

with elements of the Jacobian instead of functions, so that Theorem VII.8

could be applied.

The method we shall now describe is the further improvement by Enge

and Gaudry [117]. We give a version which is close to what is really done

when it is implemented; there exists also a more complicated version that has

a heuristic-free complexity.

We ļ¬rst present a sketch of the algorithm, and then details about each

phase are given.

Algorithm VII.4: Hyperelliptic Index-Calculus Algorithm

A divisor D1 in JC (Fq ) with known order N = ord(D1 ),

INPUT:

and a divisor D2 ā D1

OUTPUT: An integer Ī» such that D2 = Ī»D1

1. Fix a smnoothness bound and construct the factor basis.

2. While not enough relations have been found do:

Pick a random element R = Ī±D1 + Ī²D2 .

3.

If R is smooth, record the corresponding relation.

4.

Solve the linear algebra system over Z/N Z.

5.

Return Ī».

6.

VII.5.1. Construction of the Factor Basis. First a smoothness bound

B is chosen. In the complexity analysis, we shall explain how to get the best

value for B. Then the factor basis F contains all the prime reduced divisors

of weight at most B:

F = {P ā JC : P is prime, wt(P ) ā¤ B}.

This set can be constructed in a naive way: for each monic polynomial of

degree at most B, check if it is irreducible and if it is the a-polynomial of

a reduced divisor. In that case, ļ¬nd all the compatible b-polynomials and

add the corresponding divisors to F. For convenience, we give names to the

elements of F:

F = {gi : i ā [1, #F]}.

VII.5. INDEX-CALCULUS ALGORITHM FOR HYPERELLIPTIC CURVES 145

VII.5.2. A Pseudo-Random Walk. Selecting a random element R =

Ī±D1 + Ī²D2 is costly: the values of Ī± and Ī² are randomly chosen in the

interval [1, N ] and then two scalar multiplications have to be done. Using

the binary powering method, the cost is O(log N ) group operations. We use

a pseudo-random walk instead, so that the construction of each new random

element costs just one group operation.

The pseudo-random walk is exactly the same as the one which is used in

[323] in discrete log algorithms based on the birthday paradox. For j from 1

to 20 we randomly choose aj and bj in [1, N ] and compute the āmultiplierā

Tj ā aj D1 + bj D2 . In our case where the group is traditionally written

additively, āsummandā would be a better name but we stick to the classical

name. The pseudo-random walk is then deļ¬ned as follows:

ā¢ R0 is given by Ī±0 D1 + Ī²0 D2 where Ī±0 and Ī²0 are randomly chosen in

[1, N ].

ā¢ Ri+1 is deļ¬ned by adding to Ri one of the multipliers Tj . The index

j is given by the value of a hash function evaluated at Ri . Then the

representation of Ri+1 in terms of D1 and D2 is deduced from the

corresponding represenation of Ri : we set Ī±i+1 = Ī±i + aj and Ī²i+i =

Ī²i + bj for the same j.

Practical experiments [323] suggest that by taking 20 multipliers the

pseudo-random walk behaves almost like a purely random walk. In our case,

it is not really necessary to take j deterministic in terms of Ri : picking a

random multiplier at each step would do the job. The determinism in the

pseudo-random walk was necessary in [323] to use distinguished points, but

it is not used in our algorithm.

VII.5.3. Collection of Relations. Each time we produce a new random

element R = Ī±D1 + Ī²D2 , we test it for smoothness. If R is not smooth, then

we continue the pseudo-random walk to get another element. The smoothness

test is done by factoring the a-polynomial of R. If all its irreducible factors

have degree at most B, then Lemma VII.4 allows us to write R as a sum of

elements of F. The kth smooth element Sk = Ī±k D1 + Ī²k D2 that is found is

stored in the kth column of a matrix M = (mik ) which has #F rows:

Sk = mik gi = Ī±k D1 + Ī²k D2 .

1ā¤iā¤#F

We also store the values of Ī±k and Ī²k for later use.

VII.5.4. Linear Algebra. After #F + 1 relations have been collected, the

matrix M has more columns than rows. As a consequence, there exists a

non-trivial vector (Ī³k ) in its kernel. Then by construction we have

Ī³k Sk = 0 = Ī³k Ī±k D1 + Ī³k Ī²k D2 .

146 VII. HYPERELLIPTIC CURVES

Now, with high probability, Ī³k Ī²k is diļ¬erent from 0 modulo N and we can

compute the wanted discrete logarithm

Ī» = ā’( Ī±k Ī³k )/( Ī²k Ī³k ) (mod N ).

To compute the vector (Ī³k ), we could use Gaussian elimination. However, we

can do better, because the matrix is very sparse. Indeed, in the kth column, a

non-zero entry at index i means that the ith element of F has a a-polynomial

which divides the a-polynomial of the reduced divisor Sk , the degree of which

is bounded by g, the genus of C. Hence there are at most g non-zero entries in

each column. For such sparse matrices, Lanczosā™s algorithm or Wiedemannā™s

algorithm are much faster than Gaussian elimination, though they are not

deterministic.

VII.6. Complexity Analysis

In our analysis, we denote by cJ the cost of one group operation in JC (Fq ),

by cF the cost of testing the B-smoothness of a polynomial of degree g over

Fq and by cN the cost of an operation in Z/N Z. With Lagrangeā™s algorithm

and naive multiplication algorithms in Fq , we obtain cJ = O(g 2 log2 q), and if

we use NUCOMP and asymptotically fast algorithms everywhere, we obtain

cJ = O(g log q Īµ(g, log q)), where Īµ(g, log q) is bounded by a polynomial in

log g and log log q. The B-smoothness test is the ļ¬rst phase of the Cantor-

Zassenhaus algorithm called Distinct-Degree-Factorization, and its costs is

cF = O(Bg 2 log3 q) if naive multiplication algorithms are used.

VII.6.1. Analysis When g is Large Enough. We write the smoothness

bound B in the form

1

B = logq Lqg ( , Ļ) ,

2

where Ļ is some positive real number. By Theorem VII.8, the number of

B-smooth reduced divisors in JC (Fq ) is at least

1 1

q g Ā· Lq g ,ā’ ,

2 2Ļ

so that the expected number of steps in the pseudo-random walk before ļ¬nd-

ing a B-smooth divisor is about Lqg 1 , 2Ļ . We need #F + 1 relations,

1

2

therefore the expected number of steps in the pseudo-random walk that are

required to ļ¬ll in the matrix is about

11

#F Ā· Lqg , .

2 2Ļ

The cardinality of the factor basis is bounded by a constant times the number

of polynomials of degree at most B. Hence #F is O(q B ) = O(Lqg ( 1 , Ļ)).

2

VII.6. COMPLEXITY ANALYSIS 147

Putting all of this together, we see that the algorithm requires

1 1

O Lq g ,Ļ +

2 2Ļ

steps in the pseudo-random walk. At each step, we need to perform one group

operation (which costs cJ ) and the factorization of the a-polynomial of the

divisor (which costs cF ).

Once the matrix M has been computed, we need to ļ¬nd a vector of the

kernel of M . This matrix has size #F and has at most g non-zero entries per

column. The cost of Lanczosā™s algorithm is then

1

O (#F)2 g = O g Lqg , 2Ļ

2

operations in Z/N Z.

Thus the total run time of the algorithm is

1 1 1

O (cJ + cF ) Ā· Lqg + (gcN ) Ā· Lqg

,Ļ + , 2Ļ .

2 2Ļ 2

To analyze this run time and ļ¬nd the optimal value for Ļ, we remark that

cJ , cF , cN and g are all polynomials in g and log q, whereas the Lqg (...) parts

are subexponential in those parameters. Therefore it makes sense to neglect

them when trying to balance both phases. We are then left with minimizing

the function

1 1 1 1 1

, 2Ļ ā Lqg

Lq g ,Ļ + + Lqg , max(Ļ + , 2Ļ) .

2 2Ļ 2 2 2Ļ

1

ā.

An easy calculation shows that the minimum is obtained for Ļ = With

2

this choice of Ļ, we obtain an overall complexity of

1ā

O (cJ + cF + gcN ) Ā· Lqg ,2 .

2

Before putting this result into a theorem, we need to make clear an abusive

simpliļ¬cation that we did along the way: when evaluating the size ā the factor

of

1

basis, we āforgotā the ceiling sign for B. However, if logq Lqg ( 2 , 2) is close

to 0, the fact that we round it to 1 introduces an exponential factor in the

complexity analysis. Therefore, what we did makes sense only if this value

is large enough so that the rounding has no eļ¬ect on the complexity. This

requirement is fullļ¬lled if and only if g is large enough compared to log q.

Theorem VII.9 (Engeā“Gaudry[117]). There exists an algorithm that can

ā

solve the HCDLP in time O Lqg 1 , 2 when q and g grow to inļ¬nity with

2

the condition that g stays ālarge enoughā compared to log q.

The complexity of the algorithm can be made more precise by assuming

that there exists a constant Ļ‘ such that g > Ļ‘ log q. Then the run time is in

148 VII. HYPERELLIPTIC CURVES

ā

O Lqg 1 , c(Ļ‘) , where c(Ļ‘) is an explicit function that tends to 2 when Ļ‘

2

tends to inļ¬nity.

VII.6.2. Complexity When g is Fixed. In this subsection, we analyze

the behaviour of Algorithm VII.4 when g is not large compared to log q, and in

particular when g is ļ¬xed and q tends to inļ¬nity. In that case, the theoretical

ā

optimal smoothness would be logq Lqg ( 1 , 2), which tends to 0. Therefore

2

we assume from now on that B = 1.

1

By Proposition VII.7, the probability for a divisor to be 1-smooth is g! .

Note that in our algorithm, we take random elements not in the whole Ja-

cobian but in the subgroup generated by D1 . Heuristically, as long as the

1

subgroup is large enough compared to g!, the smoothness probability of g!

still holds.

The expected time for ļ¬nding enough relations to ļ¬ll in the matrix is then

(cJ + cF )g!q, and the time for the linear algebra is O(cN gq 2 ). Putting this

together, we obtain

Theorem VII.10 (Gaudry[142]). There exists an algorithm that can solve

the HCDLP in large enough subgroups in time O((g 2 log3 q)g!q + (g 2 log q)q 2 ),

when q tends to inļ¬nity and g is ļ¬xed.

VII.6.3. Very Small Genera. In the case where the genus is ļ¬xed and very

small, Harley proposed an improvement to the previous technique. Taking

B = 1 is far from the theoretical optimal value for subexponentiality, but as a

degree it cannot be reduced. Still it is possible to keep in the basis only a small

q

proportion of the weight one reduced divisors. Fix #F = Ī± , where Ī± < q is to

be determined. A 1-smooth element has a probability 1/Ī±g to be represented

over this reduced basis, therefore the expected number of steps before ļ¬nding

a relation is Ī±g g!. On the other hand, the number of required relations is

also reduced: we need only q/Ī± + 1 columns in the matrix. Therefore the

expected number of steps in order to ļ¬ll in the matrix is qg!Ī±gā’1 . Also, the

linear algebra step is made easier by the fact that the factor basis is smaller,

the cost is O(q 2 /Ī±2 ). Balancing the two phases leads to an optimal value for

Ī± which is (q/g!)1/(g+1) . With this choice, we obtain an overall complexity of

(forgetting the coeļ¬cients which are polynomials in g and log q)

2g 2

O q g+1 Ā· (g!) g+1 .

The term in g! has a negligible contribution, especially because g is supposed

ńņš. 5 |