Pierrick Gaudry

Robert Harley

Counting Points on Hyperelliptic Curves over Finite Fields

We describe some algorithms for computing the cardinality of hyperelliptic curves and their Jacobians over nite elds. They include several methods for obtaining the result modulo small primes and prime powers, in particular an algorithm a la Schoof for genus 2 using Cantor's division polynomials. These are combined with a birthday paradox algorithm to calculate the cardinality. Our methods are practical and we give actual results computed using our current implementation. The Jacobian groups we handle are larger than those previously reported in the literature.

Introduction

In recent years there has been a surge of interest in algorithmic aspects of curves. When presented with any curve, a natural task is to compute the number of points on it with coordinates in some nite eld. When the nite eld is large this is generally di cult to do.

Ren e Schoof gave a polynomial time algorithm for counting points on elliptic curves i.e., those of genus 1, in his ground-breaking paper Sch85]. Subsequent improvements by Elkies and Atkin ( Sch95], Mor95], Elk98]) lowered the exponent to the point where e cient implementations became possible. After further improvements [START_REF] Cartier | Une nouvelle op eration sur les formes di erentielles[END_REF], Ler97]) several implementations of the Schoof-Elkies-Atkin algorithm were actually written and very large nite elds can now be handled in practice [START_REF] Manin | The Hasse-Witt matrix of an algebraic curve[END_REF], Ver99]).

For higher genus, signi cant theoretical progress was made by Pila who gave a polynomial time algorithm in Pil90] (see also HI98]). However to date these methods have not been developed as extensively as the elliptic case. As a rst step towards closing this gap it is fruitful to concentrate on low genus hyperelliptic curves, as these are a natural rst generalization of elliptic curves and techniques used in the elliptic case can be adapted. Such techniques include Schoof-like methods and several others which all contribute to a practical algorithm.

We mention two possible applications of the ability to count points on low genus hyperelliptic curves. An early theoretical application was the proof that primality testing is in probabilistic polynomial time, AH92]. A practical application results from the apparent di culty of computing discrete logarithms in the Jacobian groups of these curves. In low genus, no sub-exponential algorithms are currently known, except for some very thin sets of examples ( R uc99], FR94]) and hence the Jacobian group of a random curve is likely to be suitable for constructing cryptosystems Kob89]. To build such a cryptosystem, it is rst desirable to check that the group order has a large prime factor since otherwise the logarithm could be computed in small subgroups PH78].

We restrict ourselves to odd characteristic for simplicity. We will work with models of odd degree where arithmetic is analogous to that of imaginary quadratic elds. For the even degree alternative, which is similar to real quadratic elds, see the recent paper ST99] which describes a birthday paradox algorithm optimized using an analogue of Shanks' infrastructure.

Our contribution contains several complementary approaches to the problem of nding the size of Jacobian groups, all of which have been implemented. By combining these approaches we have been able to count larger groups than previously reported in the literature.

The rst approach is an e cient birthday paradox algorithm for hyperelliptic curves. We have lled in all the details required for a large-scale distributed implementation, although the basic idea has been known for 20 years. In our implementation we also use an optimized group operation for genus 2, in which we have reduced the number of eld operations required.

The time taken grows as a small power of the eld size and this algorithm, if used in isolation, would take a prohibitive amount of time to handle large groups such as those of cryptographic size. However our version of it can take advantage of prior information on the result modulo some integer. We elaborate various strategies for collecting as much of this information as possible.

We show that when the characteristic p is not too large, the result modulo p can be obtained surprisingly easily using the Cartier-Manin operator. It provides an elegant and self-contained method based on theoretical material proved in the 1960's.

To go further, we also extend Schoof's algorithm to genus 2 curves using Cantor's division polynomials. On the basis of previous outlines existing in the literature, but not directly implementable, we elaborated a practical algorithm and programmed it in Magma. For the case where the modulus is a power of 2, we are able to bypass computations with division polynomials and use a much faster technique based on formulae for halving in the Jacobian.

The combinations of these techniques has allowed us to count genus 2 groups with as many as 10 38 elements.

We would particularly like to thank Eric Schost of the GAGE laboratory at Ecole Polytechnique for helpful discussion concerning algebraic systems. Furthermore his assistance in computing Gr obner bases was invaluable and allowed us to compute group orders modulo larger powers of 2 than would otherwise have been possible.

We also thank Fran cois Morain for many constructive comments on this paper.

Prerequisites and Notations

We will take a concrete approach, concentrating on arithmetic and algorithmic aspects rather than more abstract geometric ones.

Let g be a positive integer and let F q be the nite eld of q = p n elements, where p is an odd prime. For our purposes, a hyperelliptic curve of genus g is the set of solutions (x; y) of the equation y 2 = f(x), where f(x) is a monic polynomial of degree 2g +1 with coe cients in F q and with distinct roots 1 . Note that the coordinates may be in the base eld F q or in an extension eld.

When a point P = (x P ; y P ) is on the curve C, its opposite is the point P = (x P ; y P ). A divisor 2 is a formal sum D = P i P i of points on C. Note that points may be repeated with some multiplicity in the sum. A semi-reduced divisor is a divisor with no two points opposite. Such a divisor with k points is said to have weight k. A reduced divisor is a semi-reduced divisor of weight k g.

The Jacobian, denoted J, is the set of reduced divisors. An important fact is that one can de ne an addition operation on reduced divisors which makes J into a group, whereas this is not possible on the curve itself directly. This group law is denoted by + and will be described in the next section.

A convenient representation of reduced (and semi-reduced) divisors, due to Mumford Mum84], uses a pair of polynomials hu(x); v(x)i. Here u(x) = Q i (x

x i ) and v(x) interpolates the points P i respecting multiplicities. More precisely v = 0 or deg v < deg u, and u divides f v 2 . We say that a semi-reduced divisor is de ned over a eld F when the coe cients of u and v are in F (even though the coordinates x i and y i may be in an extension eld) and write J=F for the set of such divisors.

Most reduced divisors have weight g. The set of those with strictly lower weight is called . A divisor of weight 1 i.e., with a single point P = (x P ; y P ), is represented by hu(x); v(x)i = hx x P ; y P i. The unique divisor of weight 0, O = hu(x); v(x)i = h1; 0i, is the neutral element of the addition law. Scalar multiplication by an integer l is denoted by: l]D = D + D + + D : (l times)

(1)

We say that D is an l{torsion divisor whenever l]D = O. The set of all l{torsion divisors, including those de ned over extension elds, is denoted by J l].

We concentrate particularly on genus-2 curves and in this case the divisors in J n have the form: D = hx 2 + u 1 x + u 2 ; v 0 x + v 1 i :

(2)

1 Strictly speaking, this is the a ne part of a smooth projective model. In genus 2 every curve is birationally equivalent to such a curve provided the base eld is large enough. 2 Strictly speaking, these are degree-0 divisors with the multiplicity of the point at in nity left implicit.

1 Group Law in the Jacobian

We will sketch the group law i.e., addition of reduced divisors, using the intuitive `sum of points' notation and then describe e cient formulae for computing the law in genus 2 using Mumford's representation.

The computation of D 1 + D 2 can be viewed, at a high level of abstraction, as the following three steps:

{ form a (non-reduced) divisor with all the points of D 1 and D 2 , { semi-reduce it by eliminating all pairs of opposite points, { reduce it completely.

The third step is the only one that presents any di culty3 . When we reach it we have a semi-reduced divisor with at most 2g points. If there are g or fewer then no reduction is necessary but if there are more than g we reduce by a higher-genus analogue of the well known chord-and-tangent operation for elliptic curves.

Reduction Step

Fix g = 2 and, for the moment, consider a semi-reduced divisor R with 3 distinct points. The reduction of R is as follows.

Let y = a(x) be the equation of the parabola (or perhaps line) interpolating the three points. The roots of f a 2 are the abscissae of the intersections between the parabola and the curve. This is a quintic polynomial so there are ve intersections (including multiplicities). We already know 3 of them, the points of R. Form a divisor S with the other two, and the result of the reduction is S.

In the more frequent case where R has 4 points, choose an interpolating cubic (or lower degree) polynomial a(x) instead. Then f a 2 has degree 5 or 6 and we know 4 intersections. Form S with the others and the result is S.

In cases where some points of R are repeated, the interpolation step is adjusted to ensure tangency to the curve with su cient multiplicity. Also, in genus g > 2 the reduction step may need to be repeated several times.

In practice it would be ine cient to compute the group law this way using the representation of divisors as sums of points, since the individual points may be de ned over extension elds. By using Mumford's notation we can work entirely in the eld of de nition of the divisors.

Group Law in Mumford's Notation

Cantor gave two forms of the group law using Mumford's notation in Can87]. One was a direct analogue of Gauss's reduction of binary quadratic forms of negative discriminant, the other an asymptotically fast algorithm for high genus making clever use of fast polynomial arithmetic.

We describe an e cient algorithm, carefully optimized to reduce the number of operations required. We nd that in genus 2 doubling a divisor or adding two divisors both take 30 multiplication operations and 2 inversions, in general. Note for comparison that optimized elliptic curve operations typically take 3 or 4 multiplications and 1 inversion.

Space limits us to a brief description of the genus 2 doubling operation. Let D = hu; vi. We cover the cases that may occur, in order of increasing complexity.

Simple case: If v = 0 the result is simply O.

Weight 1: Here u(x) = x x P and v(x) = y P . The result is h(x x P ) 2 ; ax + bi, where ax + b is the tangent line at P with a = f 0 (x P )=2y P and b = y P ax P .

Weight 2: Compute the resultant r = u 2 v 2 0 + v 2 1 u 1 v 0 v 1 of u and v.

Resultant 0: If r = 0 then u and v have a root in common i.e., D has a point with ordinate 0. Isolate the other point by x P = u 1 + v 1 =v 0 , y P = v 1 + v 0 x P and return to the weight 1 case above.

General case: Consider the fact that v is a square root of f modulo u. We can double the multiplicity of all points in D by using a Newton iteration to compute a square root modulo u 2 .

{ Newton iteration: set U = u 2 , and V = (v + f=v)=2 mod U, { Get `other' roots: set U = (f V 2 )=U, { Make U monic, { Reduce V modulo U, The result is hU; V i.

Several observations help to optimize calculation with these formulae: after the rst step, V v mod u; also the division by U in the second step is exact; not all coe cients of the polynomials are really needed; nally some multiplications can be avoided using Karatsuba's algorithm.

The general addition operation is similar to doubling although the Newton iteration is replaced by a little Chinese Remainder calculation and more cases need to be handled. Since the details are somewhat tedious, we give the resulting pseudo-code and sample C code at the following Web site: http://cristal.inria.fr/~harley/hyper/ 2 Frobenius Endomorphism

In this section we collect some useful results and quote them without proof. A starting point for the reader interested in pursuing this material is IR82] and the references therein.

We rst describe properties of the q{power Frobenius endomorphism (x) = x q . Note that it has no e ect on elements of F q but it becomes non-trivial in extension elds. This map extends naturally to points, by transforming their x and y coordinates. It extends further to divisors by acting point-wise.

Crucially, this latter action is equivalent to acting on each coe cient of the u and v polynomials in Mumford's notation. When a divisor is de ned over F q , may permute its points but it leaves the divisor as a whole invariant.

Characteristic Polynomial

The operator acts linearly and has a characteristic polynomial of degree 2g with integer coe cients. In genus 2 it is known to have the form:

(t) = t 4 s 1 t 3 + s 2 t 2 s 1 qt + q 2 ;

(3) so that ( ) is the identity map on all of J, in other words: 8P 2 J; 4 (P ) s 1 ] 3 (P ) + s 2 ] 2 (P ) s 1 q] (P ) + q 2 ]P = O : (4)

The so-called Riemann hypothesis for curves, on the roots of their zeta functions, was proved by Weil and implies that the complex roots of have absolute value p q. Hence, in genus 2 the following bounds apply: js 1 j 4 p q and js 2 j 6q.

Relations Between Frobenius and Cardinalities

The Frobenius is intimately related to the number of points on the curve and the number of divisors in J, over the base eld and its extensions.

First of all, knowledge of is equivalent to that of #C=F q i for 1 i g. In genus 2 the following formulae relate them:

#C=F q = q s 1 and #C=F q 2 = q 2 s 2 1 + 2s 2 :

(5) Furthermore #J=F q is completely determined by according to the formula #J=F q = (1). An important consequence is that the group order is constrained to a rather small interval, the Hasse-Weil interval:

( p q 1) 2g #J=F q ( p q + 1) 2g :

(6)

In the reverse direction, knowledge of #J=F q almost determines for q large enough. For instance in genus 2, (#J=F q ) q 2 1 = s 2 s 1 (q +1) and the bound on s 2 given above ensures that there are O(1) possibilities.

Birthday Paradox Algorithm

To compute the group order N = #J=F q exactly we search for it in the Hasse-Weil interval which has width w close to 4gq g 1=2 . The rst few coe cients s i of can be computed by exhaustively counting points on the curve over F qi . Doing so for i I reduces the search interval to width w = O(q g (I+1)=2 ) but costs O(q I ) (see Elk98]). In genus 2 this is not useful and one simply takes w = 2b4(q + 1) p qc.

Computing the Order of the Group

Assume for the moment that we know how to compute the order n of a randomly chosen divisor D in J=F q (from now on the term \divisor" always refers to a reduced divisor). Writing e for the group exponent, we have n j e and e j N and thus N is restricted to at most d(w + 1)=ne possibilities. Usually n w and so N is completely determined.

It is possible for n to be smaller, though. In such a case we could try several other randomly chosen divisors, taking n to be the least common multiple of their orders and stopping if n > w. After a few tries n will converge to e and if e > w the method terminates.

However in rare cases the exponent itself may be small, e w. It is known that J=F q is the product of at most 2g cyclic groups and thus e p q 1 and in fact this lower bound can be attained. It is possible to obtain further information by determining the orders of divisors in the Jacobian group of the quadratic twist curve, but even this may not be su cient. We do not yet have a completely satisfactory solution for such a rare case, however we mention that the Weil pairing may provide one.

Computing the Order of One Divisor

To determine the order n of an arbitrary divisor D we nd some multiple of n, factor it and search for its smallest factor d such that d]D = O.

There are certainly multiples of n in the search interval (since the group order is one such) and we can nd one of them using a birthday paradox algorithm, in particular a distributed version of Pollard's lambda method Pol78] with distinguished points. For a similar Pollard rho method see vOW99].

Since the width of the search interval is w, we expect to determine the multiple after O( p w) operations in the Jacobian. By using distinguished points and distributing the computation on M machines, this takes negligible space and O( p w(log q) 2 =M) time 4 .

The birthday paradox algorithm is as follows.

{ Choose some distinguishing characteristic. { Choose a hash function h that hashes divisors to the range 0::19, say. { Pick 20 random step lengths l i > 0 with average roughly M p w, { Precompute the 20 divisors D

i = l i ]D. { Precompute E = c]D.
Here c is the center of the search interval. In genus 2 c = q 2 + 6q + 1. The calculation then consists of many `chains' of iterations run on M client machines:

{ Pick a random r < w and compute R = r]D. { Pick a random bit b and if it is 1 set R to R + E. { While R is not distinguished, set r := r + l h(R) and R := R + D h(R) .
{ Store the distinguished R on a central server along with r and b.

The distinguishing feature must be chosen to occur with probability signicantly less than p w=M, say 50 times less. Thus each chain takes about p w=M=50 steps and has length about w=50.

Note that chains with b = 0 visit many pseudo-random divisors in the set S 1 = f r]D j 0 r < wg and a few with larger r. Chains with b = 1 visit many divisors in S 2 = fE + r]D j 0 r < wg and a few with larger r. However the choice of E guarantees that the intersection I = S 1 \ S 2 contains at least w=2 divisors.

Now after a total of O( p w) steps have been performed, O( p w) divisors have been visited and O( p w) of them are expected to be in I. Then the birthday paradox guarantees a signi cant chance that a useful collision occurs i.e., that the same divisor R is visited twice with di erent bits b. Shortly afterwards a useful collision of distinguished points is detected at the server, between R 0 and R 1 say.

Therefore r 0 c + r 1 modulo n and nally c + r 1 r 0 is the desired multiple of n.

Beyond the Birthday Paradox

To handle larger examples than is possible with the birthday paradox algorithm alone, we precompute the Jacobian order modulo some integer. If N is known modulo m then the search for a multiple of a divisor's order can be restricted to an arithmetic progression modulo m, rather than the entire search interval 5 . In this way the expected number of operations can be reduced by a factor p m.

The algorithm outlined above needs to be modi ed as follows (we can assume that m is much smaller than w since otherwise no birthday paradox algorithm would be required!).

{ Increase the frequency of the distinguishing characteristic by a factor p m. { The step lengths must be multiples of m chosen with average length M p wm. { Replace E with z]D where z is nearest c such that z N mod m.

To compute N modulo m with m as large as possible, we will rst compute it modulo small primes and prime powers using various techniques explained in the next few sections. Then the Chinese Remainder Theorem gives N modulo their product.

To date this use of local information has speeded up the birthday paradox algorithm by a signi cant factor in practice. It should be pointed out however that while the birthday paradox algorithm takes exponential time, the Schoof-like algorithm described below takes polynomial time. Hence it can be expected that for future calculations with very large Jacobians, the Schoof part will provide most of the information. 5 Note that we could also take advantage of partial information that restricted N to several arithmetic progressions modulo m.

Cartier-Manin Operator and Hasse-Witt Matrix

We propose a method for calculating the order of the Jacobian modulo the characteristic p of the base eld, by using the so-called Cartier-Manin operator and its concrete representation as the Hasse-Witt matrix (see Car57]). In the case of hyperelliptic curves, this g g matrix can be computed by a method given in Yui78] which generalizes the computation of the Hasse invariant for elliptic curves. Yui's result is as follows:

Theorem 1. Let y 2 = f(x) with deg f = 2g +1 be the equation of a genus g hyperelliptic curve. Denote by c i the coe cient of x i in the polynomial f(x) (p 1)=2 .

Then the Hasse-Witt matrix is given by A = (c ip j ) 1 i;j g :

In Man65], Manin relates it to the characteristic polynomial of the Frobenius modulo p. For a matrix A = (a ij ), let A (p) denote the elementwise p{th power i.e., (a p ij ). Then Manin proved the following result:

Theorem 2. Let C be a curve of genus g de ned over a nite eld F p n . Let A be the Hasse-Witt matrix of C, and let A = AA (p) A (p n 1 ) . Let (t) be the characteristic polynomial of the matrix A , and (t) the characteristic polynomial of the Frobenius endomorphism. Then (t) ( 1) g t g (t) mod p :

Now it is straightforward to compute (t) modulo the characteristic p and hence #J=F q mod p, provided that p is not too large (say at most 100000).

Note that this is a very e cient way to get information on the Jacobian order, particularly when p is moderately large. Such a situation can occur in practice with elds chosen, for implementation reasons, to be of the form F p n with p close to a power of 2 such as p = 2 8 5 or p = 2 16 15.

Algorithm a la Schoof

In this section we describe a polynomial time algorithm a la Schoof for computing the cardinality of J=F q in genus 2. This algorithm follows theoretical work of Pila Pil90] and Kampk otter Kam91]. We make extensive use of the division polynomials described by Cantor Can94].

Hyperelliptic Analogue of Schoof's Algorithm

The hyperelliptic analogue of Schoof's algorithm consists of computing modulo some small primes l by working in J l]. Once this has been done, modulo su ciently many primes (or prime powers), then can be recovered exactly by the Chinese Remainder Theorem. From the bounds on s i above, it su ces to consider l = O(log q). In practice we use a few small l, determine modulo their product, and use this information to optimize a birthday paradox search as described previously.

Let l be a prime power co-prime with the characteristic. Then the subgroup of l{torsion points has the structure J l] = (Z=lZ) 2g . Moreover, the Frobenius acts linearly on this subgroup and Tate's theorem Tat66] states that the characteristic polynomial of the induced endomorphism is precisely the characteristic polynomial of the Frobenius endomorphism on J with its coe cients reduced modulo l. Hence by computing the elements of J l] and the Frobenius action on them, we can get the characteristic polynomial modulo l.

The following lemma due to Kampk otter simpli es the problem.

Lemma 1. If l is an odd prime power, then the set Jn contains a Z=lZ{basis of J l]. Thus the Frobenius endomorphism on J l] is completely determined by its action on J l] n .

Let D = hx 2 + u 1 x + u 2 ; v 0 x + v 1 i be a divisor in J n , then the condition l]D = O can be expressed by a nite set of rational equations in u 1 ; u 2 ; v 0 ; v 1 . More precisely, there exists an ideal I l of the polynomial ring F q U 1 ; U 2 ; V 0 ; V 1 ] such that D lies in J l] n if and only if f(u 1 ; u 2 ; v 0 ; v 1 ) = 0 for all polynomials f in (a generating set of) the ideal I l . In Kam91], Kampk otter gives explicit formulae for multivariate polynomials generating I l .

From now on, we can represent a generic element of J l] n by the quotient ring F q U 1 ; U 2 ; V 0 ; V 1 ]=I l . The Frobenius action can be computed for this element and it is possible to nd its minimal polynomial by brute force. The characteristic polynomial is then easy to recover (at least in the case where l is a prime) and we are done. This method due to Pila and Kampk otter has polynomialtime complexity, however it involves arithmetic on ideals which requires timeconsuming computations of Gr obner bases. In the following we propose another method which avoids the use of ideals.

Cantor's Division Polynomials

In Can94], Cantor de ned division polynomials of hyperelliptic curves, generalizing the elliptic case, and gave an e cient recursion to build them. These polynomials are closely related to Kampk otter's ideal I l , but they allow a Schoof-like algorithm to work mostly with one instead of four variables. An approximate interpretation of the phenomenon is that the division polynomials lead to a representation of I l directly computed in a convenient form (almost a Gr obner basis for a lexicographical order).

Cantor's construction provides 6 sequences of polynomials d (l) 0 ; d (l) 1 ; d (l) 2 and e (l) 0 ; e (l) 1 ; e (l) 2 such that for divisors P = hx x P ; y P i of weight 1 in general position, we get l]P = * x 2 + d (l) 1 (x P ) d (l) 0 (x P )

x + d (l) 2 (x P ) d (l) 0 (x P )

; y P e (l) 1 (x P ) e (l) 0 (x P )

x + e (l) 2 (x P ) e (l) 0 (x P ) !+ :

The degrees of these division polynomials are d0 d1 d2 e0 e1 e2 2l 2 1 2l 2 2 2l 2 3 3l 2 2 3l 2 2 3l 2 3 By lemma 1 it is su cient to consider divisors D 6 2 . In order to multiply D = hu(x); v(x)i by l we express it as a sum of two divisors of weight 1 i.e., we write D = P 1 + P 2 . These divisors are given by P 1 = hx x 1 ; y 1 i and P 2 = hx x 2 ; y 2 i where x 1 and x 2 are the roots of u(x) and y i = v(x i ). Clearly l]D = l]P 1 + l]P 2 .

The divisor D is an l{torsion divisor if and only if l]P 1 and l]P 2 are opposite divisors. This last condition is converted into a condition on the polynomial representations l]P 1 = hu P1 (x); v P1 (x)i and l]P 2 = hu P2 (x); v P2 (x)i. Indeed two divisors are opposite if their u polynomials are equal and their v polynomials are opposite. Hence the elements of J l]n are characterized by a set of rational equations in the 4 indeterminates x 1 ; x 2 ; y 1 ; y 2 , two of them involving only the two indeterminates x 1 and x 2 . Thus we get an ideal similar to I l represented in a convenient form: we can eliminate x 2 with the two bivariate equations by computing some resultants, then we have a univariate polynomial in x 1 and for each root x 1 it is not di cult to recover the corresponding values of x 2 ; y 1 and y 2 .

Details of the Algorithm

Next we explain the computation of the characteristic polynomial modulo a xed prime power l. Here we will assume that l is odd (the even case discussed in the next section).

Building an Elimination Polynomial for x 1 . We rst compute Cantor's l{division polynomials. We refer to the original paper Can94] for the recursion formulae and the proof of the construction. This phase takes negligible time compared to what follows.

The second step is to eliminate x 2 in the two bivariate equations. The system looks like E 1 (x 1 ; x 2 ) = d 1 (x 1 )d 2 (x 2 ) d 1 (x 2 )d 2 (x 1 ) = 0 ; E 2 (x 1 ; x 2 ) = d 0 (x 1 )d 2 (x 2 ) d 0 (x 2 )d 2 (x 1 ) = 0 :

(10)

The polynomial (x 1 x 2 ) is clearly a common factor of E 1 and E 2 , and this factor is a parasite: it does not lead to a l{torsion divisor 6 . We throw away this factor and consider the new reduced system, still denoting the two equations by 6 If there is another common factor of E1 and E2, we have to throw it away. This occurs when a non trivial l{torsion divisor is in . The values for the degrees assume that we are in the generic case.

E 1 (x 1 ; x 2 ) and E 2 (x 1 ; x 2 ) . Then we eliminate x 2 by computing the following resultant R(x 1 ) = Res x2 (E 1 (x 1 ; x 2 ); E 2 (x 1 ; x 2 )) = 0 :

(11)

We can then note that R(x 1 ) is divisible by some high power of d 2 (x 1 ). Indeed, if d 2 (x 1 ) = 0 then the expressions E 1 and E 2 have common roots (at the roots of d 2 (x 2 )). The power of d 2 in R is = 2l 2 2. We assume that the base eld is large enough and we specialize the system at many distinct values for x 1 . Substituting i for x 1 , the system becomes two univariate polynomials in x 2 , for which we compute the resultant r i . With enough pairs ( i ; r i ) i.e., one more than a bound on the degree of R(x 1 ) = R(x 1 )=(d 2 (x 1 )) , we can recover R(x 1 ) by interpolation. Knowing the degrees of d 0 ; d 1 ; d 2 , it is easy to get deg R(x 1 ) = 4l 4 10l 2 + 6 :

(12)

Eliminating the Parasites (optional). As previously mentioned there are l 4 divisors of l{torsion and thus the degree of R(x 1 ) is too high by a factor 4. This means that there are still a lot of parasite factors, due to the fact that we only took conditions on the abscissae x 1 ; x 2 into account and imposed nothing on the ordinates y 1 ; y 2 . Two strategies can be used: we can decide to live with these parasites and go on to the next step or we can compute another resultant to eliminate them (and get a polynomial of degree l 4 1). The choice depends on the relative speeds of the resultant computation and the root-nding algorithm.

In order to eliminate the parasites we construct a third equation E 3 (x 1 ; x 2 ), coming from the fact that the ordinates of l]P 1 and l]P 2 are opposite. We write that the coe cients are opposite, ;

(13) and this system implies that E 3 (x 1 ; x 2 ) = e 1 (x 1 )e 2 (x 2 ) e 1 (x 2 )e 2 (x 1 ) = 0. Taking the resultant between E 1 and E 3 , we get a polynomial S(x 1 ) of degree 12l 4 30l 2 + 18 whose GCD with R(x 1 ) is of degree l 4 1 (in general, a few parasites may remain in rare cases). We still denote this GCD by R(x 1 ) for convenience.

Recovering the Result Modulo l. To nd the result we factor R(x 1 ) and, for each irreducible factor, we construct an extension of F q using this factor to get a root X 1 of R(x 1 ). Then we substitute this root into E 1 and E 2 and recover the corresponding root X 2 . Using the equation of the curve we get the ordinates Y 1 and Y 2 , which may be in a quadratic extension. We get the two divisors P 1 = hx X 1 ; Y 1 i and P 2 = hx X 2 ; Y 2 i and check whether l](P 1 + P 2 ) = O or l](P 1 P 2 ) = O. If neither holds, then we started from a parasite solution and we try another factor of R(x 1 ). In the favorable case we get an l{torsion divisor D with which we check the Frobenius equation. To do so we compute s 1 ] 3 (D) + qs 1 mod l] (D) ;

(14) for every s 1 2 0; l 1] and 4 (D) + s 2 ] 2 (D) + q 2 mod l]D ;

(15) for every s 2 2 0; l 1]. We only keep the pairs (s 1 ; s 2 ) for which these are equal.

If there is only one pair (s 1 ; s 2 ) left, or if there are several pairs all leading to the same value for the cardinality modulo l, then it is not necessary to continue with another factor. Thus it is usually not necessary to have a complete factorization of R(x 1 ) and the computation is faster if one starts with irreducible factors of smallest degree.

We summarize the above in the following:

Algorithm. Computation of #J=F q modulo l.

1. Compute R(x 1 ). 2. Find a factor of R(x 1 ) of smallest degree.

3. Build P 1 and P 2 with this factor. 4. Check if P 1 + P 2 or P 1 P 2 is an l{torsion divisor. If so call it D, else go back to step 2.

5. For each remaining pair (s 1 ; s 2 ), check the Frobenius equation for D.

6. Compute the set of possible values of #J=F q from the remaining values of (s 1 ; s 2 ). If there are several values left, go back to step 2. If there is just one, return it.

Complexity

We evaluate the cost of this algorithm by counting the number of operations in the base eld F q . We neglect all the log l factors, and denote by M(x) the number of eld operations required to multiply two polynomials of degree x. The rst step requires O(l 4 ) resultant computations, each of which can be done in M(l 2 ) operations, and the interpolation of a degree O(l 4 ) polynomial which can be done in M(l 4 ) operations. For the analysis of the remaining steps, we will denote by d the degree of the smallest factor of R(x 1 ) that allows us to conclude. We assume moreover that the most costly part of the factorization is the distinct degree factorization (which is the case if d is small and if the number of factors of degree d is not too large). Then the cost of nding the factor is O(d log(q))M(l 4 ). Thereafter the computation relies on manipulations of polynomials of degree d and the complexity is O(l + log(q))M(d), where l re ects the l possible values of s 1 and of s 2 and log(q) re ects the Frobenius computations. Hence the (heuristic) overall cost for the algorithm is O(l 4 )M(l 2 ) + O(d log q)M(l 4 ) + O(l 2 + log q)M(d) (16) operations in the base eld. Now we would like to obtain a complexity for the whole Schoof-like algorithm.

For that we will keep only the primes l for which d = O(l); this should occur heuristically with a xed probability (this is an analogue of `Elkies primes' for elliptic curves). Then we have to use a set of O(log q) primes l, each of them satisfying l = O(log q). Moreover we will assume fast polynomial arithmetic and thus M(x) = O(x) (ignoring logarithmic factors). Hence the cost of the algorithm is heuristically O(log 7 q) operations in F q . Each operation can be performed in O(log 2 q) bit operations using classical arithmetic and we get that the complexity of the Schoof-like algorithm is O(log 9 q).

Remark. This analysis is heuristic, but one could obtain a rigorous proof that the algorithm runs in polynomial time. The algorithm could also be made deterministic by avoiding polynomial factorizations. However in both cases the exponent would be higher than 9.

6 Lifting the 2-Power Torsion Divisors

In this section, we will show how to obtain some information on the #J=F q modulo small powers of 2. Factoring f gives some information immediately.

To go further we iterate a method for `halving' divisors in the Jacobian. This quickly leads to divisors de ned over large extensions, so that the run-time grows exponentially. In practice we can use this technique to obtain partial information modulo 256, say.

The divisors of order 1 or 2 are precisely the D = hu(x); 0i for which u(x) divides f(x) and is of degree at most g. When f has n irreducible factors, then it has 2 n factors altogether. Exactly half of them have degree at most g, since f is square-free of degree 2g + 1. Hence the number of such divisors is 2 n 1 , and 2 n 1 j #J=F q . Furthermore, when f is irreducible then the 2{part is trivial and #J=F q is odd.

Halving in the Jacobian

Let D = hu(x); v(x)i be a divisor di erent from O. We would like to nd a divisor such that 2] = D. Note that there are 2 2g solutions, any two of which di er by a 2{torsion divisor. In general, is de ned over an extension of the eld of de nition of D.

Writing = hũ(x); ṽ(x)i, we derive a rational expression for the divisor 2] using the formulae of section 1. Then equating this expression with D, we get a set of 2g polynomial equations in 2g indeterminates ũi and ṽi with 2g parameters u i and v i . There are g 2 such systems corresponding to the di erent possible weights of D and . We consider the most frequent case where D and are both of weight g.

The corresponding system has at most 2 2g solutions and these can be obtained by constructing a Gr obner basis for a lexicographical order, factoring the last polynomial in the basis and propagating the solution to the other polynomials.

All this can be done in time polynomial in log q provided that the divisor D we are dealing with is de ned over an extension of bounded degree of F q .

In order to speed up the computations in the case where D is de ned over a large extension, we can avoid repeated Gr obner-basis computations and instead compute a single generic Gr obner basis for the system, where the coe cients of D are parameters. As the halving is algebraic over F q (because the curve is de ned over F q ), the generic basis is also de ned over F q . After this computation we can halve any divisor D, even when de ned over a large extension, by plugging its coe cients into the generic basis to get the specialized one.

We are indebted to Eric Schost who kindly performed the construction of this generic Gr obner basis for the curves we studied Sch]. For his construction, he made use of the Kronecker package Lec99] written by Gr egoire Lecerf. This package behaves very well on these types of problem (lifting from specialized systems to generic ones), and it is likely that we would not have been able to do this lifting by using classical algorithms for Gr obner-basis computations.

Example. Let C be de ned by y 2 = x 5 + 1597 x 4 + 1041 x 3 + 5503 x 2 + 6101 x + 1887 ;

(17) over the nite eld F p with p = 10 17 + 3. We will search for all rational 2{power torsion divisors i.e., those de ned over F p . Two irreducible factors of f(x) have degree at most 2, they are f 1 = x + 28555025517563816 and f 2 = x + 74658844563359755 ; Thus there are three rational divisors of order two: P 1 = hf 1 ; 0i, P 2 = hf 2 ; 0i and P 1 + P 2 . The halving method applied to P 1 nds four rational divisors of order 4. They are hu; vi and hu; vi where: u = x 2 + 1571353025997967 x + 12198441063534328 v = 32227723250469108 x + 68133247565452990 and: u = x 2 + 70887725815800572 x + 94321182398888258 v = 42016761890161508 x + 3182371156137467 :

There are 16 solutions altogether but the others are in extension elds (the Gr obner bases are too large to include them here!) Applying the method to P 2 and to P 1 + P 2 nds no further rational 4{torsion divisors. By continuing in the same manner one nds 8 divisors of order 8, 16 of order 16, 32 of order 32 and no more. Thus the 2{part of the rational Jacobian is of the form (Z=2) (Z=32) and hence #J=F p 64 mod 128. This type of exhaustive search in the base eld determines the exact power of 2 dividing #J=F p . In the next section we show how to nd information modulo larger powers of 2.

6.2 Algorithm for Computing #J=F q mod 2 k Next we go into extension elds to nd some 2 k {torsion divisors and we substitute them into , the characteristic equation of the Frobenius endomorphism, to determine values of its coe cients modulo 2 k and hence the value of #J=F q mod 2 k , for increasing k.

Algorithm (for g = 2).

1. Factor f to nd a 2{torsion divisor. Halve it to get a 4{torsion divisor D 4 .

2. Find the pair (s 1 ; s 2 ) mod 4 for which (D 4 ) = O. Set k to 2.

3. Compute the generic Gr obner basis for halving (weight 2) divisors in the given Jacobian.

4. Build a 2 k+1 {torsion divisor D 2 k+1 by substituting the coe cients of D 2 k in the system, computing a root of the eliminating polynomial in an extension of minimal degree, and propagating it throughout the system.

5. For each pair (s 1 ; s 2 ) mod 2 k+1 compatible with the previously found pair modulo 2 k , plug D 2 k+1 into and nd the pair for which (D 2 k+1 ) = O.

6. Set k = k + 1, and go back to Step 4.

Note that this is an idealized description of the algorithm. In fact there will frequently be several pairs (s 1 ; s 2 ) remaining after checking the Frobenius equation for one 2 k {torsion divisor. We can eliminate false candidates by checking with other 2 k {torsion divisors. It can be costly to eliminate all of them when the required divisors are in large extensions; an alternative strategy is to continue and expect the false candidates to be eliminated later using 2 k+1 {torsion divisors.

In this algorithm, we could skip step 3 and compute speci c Gr obner bases at each time in step 4. However, the generic Gr obner basis is more e cient and allows one to perform one or two extra iterations for the same run-time.

Combining these Algorithms | Practical Results

We have implemented all these algorithms and tested their performance for real computation. Some of them were written in the C programming language, and others were implemented in the Magma computer algebra system BC97].

Prime Field

In the case where the curve is de ned over a prime eld F p , where p is a large prime, we use all the methods described in previous sections except for Cartier-Manin. We give some data for a `random' curve for which we computed the cardinality of the Jacobian. Let the curve C be de ned by y 2 = x 5 + 3141592653589793238 x 4 + 4626433832795028841 x 3 + 9716939937510582097 x 2 + 4944592307816406286 x + 2089986280348253421 ;

over the prime eld of order p = 10 19 + 51. The cardinality of its Jacobian is #J = 99999999982871020671452277000281660080 ;

(19) and the characteristic polynomial of the Frobenius has coe cients: s 1 = 1712898036 and s 2 = 11452277089352355350 : The rst step of this computation is to factor f(x). It has 3 irreducible factors, thus we already know that #J 0 mod 4.

The second step is to lift the 2{power torsion divisors. The computation of the generic halving Gr obner basis (done by E. Schost) took about one hour on an Alpha workstation. Then we lifted the divisors several times and checked the Frobenius equation. In the following table we give the degree of the extension where we found a 2 k {torsion divisor, and the information on #J that we got (timings on a Pentium 450). The run-time for l = 3 is surprisingly large in this table. For our curve, an unlucky event occurs, which becomes rare as l increases. Indeed, after testing the Frobenius equation for all the 3{torsion divisors several candidates (s 1 ; s 2 ) still remain, yielding several possibilities for #J mod 3. What this means is that the minimal polynomial of is not the characteristic polynomial. Each remaining candidate for (s 1 ; s 2 ) gives a multiple of the minimal polynomial. By taking their GCD we obtain the exact minimal polynomial, from which we can deduce the characteristic polynomial7 and #J mod 3.

In our case, there are 3 pairs left after testing all the 3{torsion points, leading to the following candidates for #J mod 3.

(s1; s2) mod 3 #J mod 3 (t) mod 3 (0; 2) 1 t 4 t 2 + 1 (1; 2) 2 t 4 t 3 t 2 t + 1 (2; 2) 0 t 4 + t 3 t 2 + t + 1

The third case is impossible because if #J 0 mod 3 then we would have found a rational 3{torsion divisor earlier. In order to decide between the two rst cases we determine the minimal polynomial, which is t 2 + 1 and thus the characteristic polynomial must be (t 2 + 1) 2 and nally #J 1 mod 3.

However to do this we have to build all the 3{torsion divisors. This explains why the running time is higher than for l = 5, where we found a rational 5{ torsion divisor and immediately deduced that #J = 0 mod 5. The nal step is the birthday paradox computation. The width of the Hasse-Weil interval is roughly 2:5 10 29 . The search space is reduced by a factor 2 8 3 5 7 11 13 = 3843840 leaving 6:6 10 22 candidates. The search was performed on ten Alpha workstations working in parallel and calculated 5 10 11 operations in the Jacobian. On a single 500 MHz workstation, this computation would have taken close to 50 days.

Non-prime Fields

Let C be a genus 2 curve de ned over F p n , where p is a small odd prime. We assume that C is not de ned over a small sub eld, for in that case it is easy to compute (t) using a theorem due to Weil.

Here the rst step is to use Cartier-Manin to get (t) mod p quickly and then continue as before, except that we avoid l = p in the Schoof part.

Examples: We did not try to build big examples, however we give two medium ones. For the rst, let the curve C be de ned by y 2 = x 5 + x 4 + x 3 + x 2 + tx + 1 ;

(20) over the nite eld F 3 30 = F 3 t]=(t 30 + t 1). The cardinality of its Jacobian is #J = 42391156018493425614913594804 :

(21) The second example illustrates the advantage given by Cartier-Manin in a favorable case where p = 2 16 15. Let the curve C be de ned by y 2 = x 5 + x 4 + x 3 + x 2 + x + t ;

(22) over the nite eld F p 4 = F p t]=(t 4 17). The Cartier-Manin computation gave us #J 58976 mod p in 17 minutes, and nishing using our other methods gave #J = 339659790214687297284652908385855015466 :

(23)

8 Perspectives for Further Research

The present paper reports on practical algorithms for counting points on hyperelliptic curves over large nite elds and on implementations for genus 2.

Although it is now possible to deal with almost cryptographic-size Jacobians, there is still a substantial amount of work to be done. Some improvements or generalizations seem to be accessible in the near future, whereas others are still quite vague. Among them we would like to mention:

{ Extension of the algorithm to even characteristic. This is only a matter of translating the formulae, in order to deal with an equation of the form y 2 + h(x) y = f(x). The Cartier-Manin part and the lifting of the 2{torsion should merge, giving an e cient way to compute the result modulo 2 k . For the Schoof-like part, the formulae of Cantor's division polynomials have to be adapted, which does not appear to be too di cult.

{ Extension of the Schoof-like algorithm to genus g > 2. The main di culty is that it does not appear possible to avoid manipulation of ideals.

{ More use could certainly be made of the Jacobian of the twist curve.

{ We believe that it may be possible to lift the curve to a local eld with residue eld F q and use Cartier-Manin to compute (t) modulo small powers of the characteristic. We do not yet know how to compute the lift, however.

{ A major improvement would be to elaborate a genus 2 version of the Elkies-Atkin approach for elliptic curves, which would lead to computations with polynomials of lower degree. We conjecture that it is possible to work with degrees reduced from O(l 4 ) to O(l 3 ). The rst task is to construct modular equations for Siegel modular forms, instead of classical ones. This requires a description of isogenies for each small prime degree, which can be given by lists of cosets under left actions of the symplectic group Sp 4 (Z) instead of the classical modular group SL 2 (Z). Starting points for studying the relevant forms and groups include Fre83] and Kli90]. This will be explained in more detail elsewhere Har]. All the above is the subject of active research.

In a more classical treatment the reduction would be described as choosing a representative for an equivalence class of degree 0 divisors modulo linear equivalence, where linearly equivalent divisors are those that di er by a principal divisor.

Using classical algorithms for eld arithmetic.

See Kam91] for more about this.


Julien Bénabès
email: julien.benabes@irccyn.ec-nantes.fr

Fouad Bennis

Emilie Poirson

Yannick Ravaut
email: yannick.ravaut@fr.thalesgroup.com

An Interactive-based approach to the Layout design optimization

Keywords: Layout problem, Interactive optimization, Genetic algorithm

Layout design plays an important role in the design and usability of many engineering products and systems. Because of the great complexity of most of industrial layout problems, the decision of the acceptable layout is a hard and critical task since the special layout can have a significant consequence on the global performances. Thus, in order to propose to the designer an optimal spatial arrangement in a reasonable time, this paper develops an interactive optimization strategy that is tested on the facilities layout problem of a shelter. This problem is innovative because it introduces the concept of space of accessibility.

INTRODUCTION

Layout problem is inherently a multidisciplinary task [START_REF] Giassi | Multidisciplinary design optimization and robust design approaches applied to concurrent design[END_REF]. It covers all the aspects of the product design life cycle from the conceptual to the detailed stage and makes necessary the collaboration between experts of technical and economical disciplines. In fact, layout design is usually formulated as an optimization problem: find the best arrangement (location and orientation) of components in a given available space satisfying geometrical and functional constraints. A non-overlapping constraint is basically a common geometrical constraint for all 3-Dimensional Layout problems, while alignment, orientation or gathering components refer to functional constraints. Because of the geometrical complexity, the 3-Dimensional layout problem optimisation is non-linear and NP-hard. It means that the problem is intrinsically harder than those that can be solved by a nondeterministic Turing machine in polynomial time. The objective and constraints evaluation is generally time consuming. It is essential to distinguish between Cutting and Packing (C&P) problems and 3-Dimensional Layout problems. In C&P problems, components are only geometrically related to each other, whereas in layout problems, components are geometrically and functionally related and connected. This difference leads different tools and methods to solve each class of problem being aware of the common nonoverlap constraints in the two problems. Typologies of C&P problems have been proposed [START_REF] Dyckhoff | A typology of cutting and packing problems[END_REF] but as far as we know, there is no general typology of layout problems. Amine Drira has described a tree representation of facility layout problems [START_REF] Drira | Facility layout problems: A survey[END_REF] that depends on design constraints and objectives of the location of facilities inside a plant. Actually, layout problems can be divided into several kinds of specific problems, which have their own solving method. C&P problems can be assimilated as a particular application of each specific problem. Layout problems can also be classified according to three criteria: the compactness of the problem, the number and type of design constraints and objectives and the geometrical complexity of the design components. For example, the container loading problem and the engine layout design are not in the same category of problem. It means that, for the two problems, the compactness is important but in the engine layout design, constraints and objectives are multiple (non-overlap and functional constraints, accessibility objective...) and the different parts of the engine have complex lines. The container loading problem does not have the same characteristics. The formulation of layout problems uses mono-objective or multi-objective optimisation. The designer can make an early decision by using an aggregation function in order to transform a multi-objective into a mono-objective one. This approach is only effective when all data and information on the aggregation are available or if the designer is familiar with the specific layout problem. In this paper, we use multi-objective optimization. The decision on the preferences between objective functions is delayed so that the designer can use the Pareto-front in order to select the most appropriate solution. In this approach, the designer has to simultaneously optimize two or more conflicting objectives subject to constraints. The general formulation of an optimization problem can be written as the equation 1, where m is the number of objective functions and n the number of design variables.

) ( 0 ) ( ) ,..., , ( ) ( ) ,..., , ( 2 1 2 1 x h x g f f f x F x x x m n x (1)
The designer has to compare two solutions represented by two vectors of objectives ) ,..., , ( ) (

2 1 mU U U f f f U F = and ) ,..., , ( ) ( 2 1 mV V V f f f V F =
where iU f is the ith component of the vector of objectives F for the design variable U. In fact, U dominates V (Pareto dominance) if U is as good as V for all the objectives and U is better than V for at least one objective. Mathematically, this can be formulated by:

⎪ ⎩ ⎪ ⎨ ⎧ < ∈ ∃ ≤ ∈ ∀ jV jU iV iU f f n j f f n i ] ,..., 1 [ ] ,..., 1 [ (2)
Multi-objective optimization searches for the set of nondominated points (assimilated to Pareto-optimal points in the next sections of this paper) in the objective space given by efficient solutions. Figure 1 represents the Pareto front for an optimization problem defined by two objectives ) 2 1 f f min , (min where U dominates V.

Figure 1: Pareto front of a multi-objective problem. One finds multiple search algorithms to solve layout optimization problems in two or three dimensions. Traditional optimization approaches for three dimensional layout problems are described by Jonathan Cagan [START_REF] Cagan | A survey of computational approaches to the three-dimensional layout problems[END_REF]. They use genetic algorithms [START_REF] Yi | Vehicle configuration design with a packing genetic algorithm[END_REF], simulated-annealing algorithms [START_REF] Szyckman | A simulated annealing-based approach to the three-dimensional packing[END_REF][START_REF] Szyckman | Constrained threedimensional component layout using simulated annealing[END_REF] or extended pattern search algorithms [START_REF] Su | An extended pattern search algorithm for three-dimensional component layout[END_REF]. Most search algorithms are developed for a specific problem and they provide an effective optimization strategy for this kind of problem. Therefore, they are not generic and can not be adapted to other layout problems. In this paper, the proposed method is based on a generic technique for solving layout problems. The design strategy uses a Genetic Algorithm coupled with a Separation Algorithm, in order to ensure a good diversity of solutions computed by the algorithm, and allows the designer to interact with the Pareto-optimal solutions. This paper is organized as follows: in section 2, the synopsis of the proposed optimization method is presented. In section 3, the proposed method is tested on the layout problem of facilities inside a shelter. Optimization problem formulation and results obtained by the method are described and analyzed. Sections 4 and 5 are dedicated to an outlook on future work and the conclusion.

PROBLEM SOLVING STRATEGY

The constraint space of design variables is highly composite because of the geometric complexity, the nonoverlap constraints and the relative location between the components. This property leads us to recognise that there is no choice but using stochastic or heuristic techniques for a 2-D or 3-D layout problem. This kind of technique makes it possible to explore efficiently the design space and avoid a local optimum. We also recognise that a multi-objective Genetic algorithm is suitable for this problem. Basically, a genetic algorithm search uses the mechanics of natural selection and natural genetics to evolve a population of initial solutions into a near-optima solution. The common idea behind these algorithms is the same: given a set of individuals from an initial population, the genetic algorithm uses basically 3 operators in order to create a set of new candidate solutions. This process can iterate until candidates with sufficient good fitness are obtained. The three operators are namely: selection, mutation and crossover or recombination. Mutation and recombination act to create the necessary diversity and novelty and selection acts as a force of pushing quality (fitness). Repeated selection from the same population would produce nothing better than multiple copies of the best individual originally in it. For improvement to be able to occur, some novelty must be introduced into the population between selection steps. The genetic operators modify the selected parents by manipulating their genotype. Since the genetic algorithm is based on stochastic operators and parameters, the progression of multiobjective layout optimisation is time consuming. It depends also on the number of design variables and the number of components and the types of design constraints. In order to improve the algorithm we introduce two new steps into the global process of the genetic algorithm: separation techniques and interaction with the designer. The initial random population of the GA is improved using the separation technique and interaction. In fact, the initial random population of the GA leads to a high number of overlap components thus the GA fails to find efficient candidates. The objective of the search algorithm proposed in this paper is to generate a uniformly distributed global Pareto front for layout optimization problems. Our strategy consists of initializing the multi-objective optimizer with a population of individuals which have been locally modified by a separation algorithm and designer interaction in order to reduce the violation of placement constraints. Thus, the strategy is based on three complementary approaches, which are clearly separated: 1. Firstly, the generation of a database of mixed designs that respect non-overlap constraints. 2. Secondly, the optimization of this database by considering all the design objectives (with a Genetic Algorithm). 3. At the end of the GA iteration process, the interactive choice of the suitable solution by the designer is conducted using the non dominated Pareto front generated by the GA.

Separation Algorithm

Several separation algorithms have been proposed [START_REF] Imamichi | A multi-sphere scheme for 2D and 3D packing problems, SLS[END_REF][START_REF] Imamichi | Designing algorithms with multi-sphere scheme[END_REF]. However, the key idea is always the same: given a configuration that doesn't satisfy location constraints, the objective of the separation algorithm is to minimize the non-respect of overlap between components and protrusion (overlap between components and the non allowed space).

For solving simple layout problems in two dimensions, the separation problem is formulated as an unconstrained minimization problem defined by:

⎩ ⎨ ⎧ ≠ ∈ = ∑ j i n j i ij A x F , ] ,..., 1 [ , ) ( min Algo) (Sep (3)
where ij A represents the intersection area between the components i and j. Consequently, it is possible to define the violation of placement constraints F as the total sum of intersection areas between the different elements which make up the layout design. For example, let us consider that all the items of the layout design are rectangles, the intersection area between two items is equal to:

)] 2 , 2 ( ) 2 , 2 ( , 0 [ )] 2 , 2 ( ) 2 , 2 ( , 0 [ j j i i j j i i j j i i j j i i ij L y L y L y L y l x l x l x l x A - - + + × - - + + = max - min max max - min max (4)
where ) , (

i i y x
are the coordinates of the geometric center of the rectangle i. i L and i l represent respectively the length and the width of the rectangle i.

f 1 f V f 2 U f 2 U f 1 V f 1 U V
The algorithm used to minimize F is based on the Broyden-Fletcher-Goldfarb-Shanno (BFGS) method. This algorithm computes a finite-difference approximation of the gradient and the hessian of the function F in order to locally modify the optimization variables and to minimize F. The algorithm stops after a certain number of iterations. in the container. It means that the algorithm searches the optimal configuration that reduces the violation of placement constraints F, which has been previously defined as the total sum of intersection areas between the square items and the container. Figure 2 shows simulations results, considering different values of the problem density.

Multi-objective optimization

The multi-objective optimizer is in charge of exploring efficiently the search space to propose trade-off solutions. The proposed method uses the Multi-Objective Genetic Algorithm (MOGA-II) [START_REF] Poles | MOGA-II, an improved multiobjective genetic algorithm[END_REF]. MOGA-II is an efficient Multi-Objective Genetic Algorithm (MOGA) [START_REF] Deb | Multi-objective genetic algorithms: problems difficulties and construction of test problems[END_REF] that uses a smart multi-search elitism. This new elitism operator is able to preserve some excellent solutions without bringing premature convergence to local-optimal frontiers. For simplicity, MOGA-II requires only very few user-provided parameters. Several other parameters are internally settled in order to provide robustness and efficiency to the optimizer. Three genetic operators are used to generate new solutions:

• Directional Cross-over: the cross-over operator is a method of recombination where parents produce offspring by sharing information. The aim of this operator is to obtain individuals with better characteristics while maintaining the diversity of the population. • Selection: this value gives the probability that design configurations are not changed during the evolution. In order to maintain a good diversity between points, this parameter should be kept small, • Mutation: this value gives the probability that a design configuration is randomly changed.

For the application studied in section 3, the genetic operators have been set as shown in

Interactive process

In general the development of an engineering object is considered as a single process involving multicriteria identification of the mathematical model followed by multicriteria optimization of the object design on the basis of this mathematical model. The process of statementsolution of engineering design problems without the interference of the design is impossible. In solving the design problem, the designer almost always has to correct either the mathematical model, the dimension of the vectors of design variables and criteria, the design variable ranges, and so on. This creative process of correcting an initial statement is natural when solving engineering problems. The direct participation of the designer in the construction of the feasible design and nonformal analysis are the essential stage of the search for the optimal design [START_REF] Serna | Some ways of implementation of interactive design[END_REF]. The simulation tools provide powerful solutions for planning and designing of complex mechanical systems. The problem with these is the representation and the interpretation of the results by the engineer. Important for the engineer is not only the value of the point but its variation and the information about the most favorable directions. The exploitation of the results is not obvious and the link with the performance value of the real phenomenon is not trivial. When one analyzes the communication between the operator and the computer, it can be perceived that the operator immersion in the digital model is very weak. It was pointed out in section 1, that the decision can be made earlier in the design optimisation process by the creation of an aggregation of the objective function or later using the non-dominated points generated by the genetic algorithm. The proposed strategy uses two interactive steps. The first one is the interactivity of the designer for the selection of feasible solution to present as an initial population of the genetic algorithm. The first interaction step is limited to the geometrical non-overlap of components. The GA algorithm can be stopped after a fixed number of iterations. Since all the non-dominated points are potentially good acceptable solutions of the given Layout problem, the designer has to explore the set of these points and select the best solution. It is well recognised by the expert of the optimisation that it is always very hard to express all the designer requirements using only the objective function. Several subjective functions are qualitative and can not easily be expressed explicitly by using a numerical value. In order to take into account this subjective aspect of the layout problem, selected solutions of the non-dominated front are presented to the designer and an interactive numerical environment is used to support the decision. Since the number of non-pareto points can be to very large, we also use a reduction approach that regroups all the equivalent solutions. Then for each presented solution, the designer can act on the environment and locally modify the position or orientation of the appropriate components. This step is different to the initial interaction since the designer has also the evaluation of all the objective function. The interactivity is not only limited to the geometrical visual evaluation but the value of the objectives function and there deviation are helpful at this step.

APPLICATION

Problem description

In this section we consider an application of the proposed strategy to find the optimal layout of facilities in a shelter.

Several components have to be arranged in the shelter including electric and energetic cabinets, desks and electrical boxes. The CAD model of the shelter is presented Figure 3. Here, the virtual components represent the spaces of accessibility of the cabinets and the desks. For example, the space of accessibility of the cabinet is the required space to insert some materials into the cabinet. These spaces are symbolized in Figure 4 by dotted rectangles. With this problem formulation, design constraints depend on the category of components. It means that overlap is allowed between two spaces of accessibility, taking account that operations of materials loading are sequentially made whereas overlap has to be minimized between two material components. Moreover, the space, represented by hatching in Figure 4, is the space below the air-conditioner where no cabinet can be placed. This space is also a virtual component that is fixed during the optimization process The dimensions of the shelter are 2150mm x 2740mm. The density of this configuration, without considering the spaces of accessibility of the different components, is equal to 50%. If the spaces of accessibility are considered, this density increases up to 90%.

Problem formulation

Problem formulation is a very important step of the optimization process. The optimization problem studied here is an under constrained multi-objective problem. Let us see how the variables, constraints and design objectives are defined. : the coordinates of each element (a continued variable along X axis and another one along Y axis) and the rotation angle (one discreet variable along Z axis). Consequently, the number of optimization variables for this problem is equal to 24 (= 8 items x 3 coordinates). Because of the rotation of each component, variables X and Y are bounded according to the following relation (for the variable i X for example):

) , ( ) , ( i i sh i i i L l l X L l min min - < < ( 5
)
where sh l represents the width of the shelter. Here, i l is the dimension of the component along X axis (it does not have to be confused as the width of the component i).

Design constraints

The design constraints of this layout problem are nonoverlap constraints. They are divided in four categories, according to the following classification:

• Non-overlap constraints between components (C1).

• Non-overlap constraints between components and spaces of accessibility (dotted rectangle represented in Figure 4) (C2). • Non-protrusion constraints between components, spaces of accessibility and the shelter (C3). • Non-overlap constraints between cabinets and the space below the air-conditioner (hatching represented in Figure 4) (C4). The rectangular shape of components simplifies the formulation of design constraints. Thus, the non-overlap constraint between the rectangles i and j is equal to the intersection area between the component i and j (in cm²). This area has been defined in equation 4. Actually, the objective of the separation algorithm is defined as:

4 3 2 1 C C C C F + + + = (6)

Design objectives

In collaboration with the company experts of this specific problem, we have considered for this optimization problem the five following design objectives:

• To minimize the distance between the center of gravity of components and the geometrical center of the shelter, in order to balance the masses inside the shelter (O1). • To maximize the distance between the cabinet 1 and the cabinets 2 and 3 and the electrical box 2, in order to limit interactions between energy and electric network (O2, O3, O4). • To minimize the distance between the electrical box 2 and one of the shelter's walls, in order to establish a connection with exterior (O5). The design objectives O2, O3, O4 and O5 are formulated by the distance between the centers of gravity of elements. For example, the distance ij d between the components i and j is equal to:

)² ( )² ( i j i j ij y y x x d - + - = (7)
where x and y are the coordinates of the centers of gravity of the items i and j. Let us consider the coordinates of the center of gravity of all the elements which are placed in the shelter. These coordinates are equal to:

∑ ∑ ∑ ∑ = = = = × = × = N i i N i i i gra N i i N i i i gra m m y Y m m x X 1 1 1 1 ) ( , ) ( (8)
where N is equal to the number of elements which have a mass: the cabinets, the desks, the electrical boxes and the air-conditioners. Then, by considering equation 7, the objective 1 (O1) is computed. More designer's knowledge could be incorporated in the layout problem formulation. It means for example, in the configuration design of the shelter studied in this paper, the design objective O5 can be deleted and the degree of freedom of the electrical box 2 can be reduced, in order to force it to displace along one of the walls of the shelter. The designer's contribution for formulation should simplify the search of feasible solutions by reducing the number of possible solutions.

Results and analysis

The resolution of this optimization problem has been firstly realized only with the multi-objective optimizer MOGA-II. The algorithm has been randomly initialized with a population of 240 designs. Most of these initial designs did not respect the non-overlap constraints because they have been randomly generated. Because of the great density of this layout problem, only one or two feasible configurations were computed for each simulation. A configuration is defined as: Design j is a new configuration if it differs from the design i by at least one of the following criteria: • One of the components of the layout has been displaced from at least Δ mm along one of the axis X or Y, ( Δ is set to 500 mm in this application).

• One of the components has been rotated.

• The minimum difference between the objective values of the two designs is bigger than a limit, for example, fixed at 10 cm. These results lead us to use the method proposed in this paper in order to generate, with only one optimization simulation, a set of well distributed Pareto-optimal designs. Thus, the results obtained for each step of the method are described here:

1. Separation algorithm and first interaction with the designer: the algorithm has been randomly initialized with designs that do not respect non-overlap constraints. Then, a set of 162 "feasible"' (it means that respect design constraints) designs have been computed. This population then has been completed with 78 individuals randomly generated in order to create the first population (240 individuals) and to guarantee the diversity of the genetic algorithm.

2. Multi-Objective Genetic Algorithm: the algorithm has searched optimal solutions by considering the design constraints and all the objectives of the problem. Then, after a hundred generations, a set of 172 feasible configurations have been computed. 41 of these solutions are Pareto-optimal designs. Figure 5 shows three of these Pareto-optimal configurations and the initial solution. It's important to mention that this initial configuration is an intuitive solution which has been generated only by considering geometric aspects.

Figure 5: Solution computed by the algorithm. The Pareto-optimal solutions, shown in Figure 5, point out a problem that is present in all the solutions computed by the algorithm. In fact, Figure 5 shows that, even though the design is a feasible configuration that respects all the design constraints, some facilities (for example facilities n°2, 4, 5, 8 in the first solution) are not accessible from the shelter's entry. It means that a design constraint or an objective is missing in the layout problem formulation, in order to characterize the accessibility to a facility from the shelter's entry. Consequently, the problem description has been changed and a free corridor, located in the middle of the shelter, has been added. In fact, this corridor is a space of living that can be considered as a fixed virtual component, where all the material components can not be placed. This space is going to resolve the problem of accessibility inside the shelter. Then, with the separation algorithm, a set of feasible designs have been computed and by interacting directly with them and by relaxing the design constraints (until 150 cm²), the designer has selected 78 different designs. This population then has been completed with 162 individuals randomly generated. Next, the genetic algorithm has generated 14 feasible configurations whose 7 are Pareto-optimal solutions. 3. Interactive decision making: the 7 Pareto-optimal designs do not dominate the initial solution. On the other hand, the initial solution does not dominate either the solutions computed by the proposed method. Actually, it means that the designer is the only person who can make the final design choices.

In order to make a decision, the method provides to the designer an interactive geometric and numeric visualization of the designs. The designer can explore the set of non-dominated solutions, compare their objective values and interact with them. Figure 6 represents the interface that allows the designer to visualize at the same time the layout design and the associated design constraints and objectives values.

The designer can also compare two solutions. In fact,

when the designer displaces one component or changes its direction, the design constraints and objectives are automatically actualized.

Figure 6: Interactive layout interface. Actually, among the 7 Pareto-optimal solutions computed by the algorithm, let us focus on the third solution. Locally changing the location of some components of this design improves its performances. Figure 7 shows this solution 3 and the solution that results from the modifications made by the designer.

RESTROSPECT AND PERSPECTIVE

This paper has introduced a new interactive optimization strategy for solving layout problems. The method can be divided into several steps, as shown in Figure 8. Firstly, a population of designs, randomly initialized, is optimized by the separation algorithm. Then, the designer interacts with the solutions computed by the algorithm and selects some individuals according to the design constraints. Secondly, the new population is optimized by the multi-objective optimizer by considering all the design objectives. Then, the designer can locally modify some computed designs in order to improve their objectives. Actually, our strategy has the innovative particularity to allow the designer to interact with the optimization process in order to improve the performances of the Pareto-optimal designs and to keep a good diversity in computed solutions. The application, which has been studied in this paper, has an innovative problem formulation because it introduces the concept of space of accessibility, which can be considered as a virtual component of the layout design. Moreover, this application emphasizes the problem of accessibility to a facility from the container's entry. This problem has been resolved by inserting the user's knowledge in the problem description. Here, the designer decided to insert a free corridor in the shelter in order to keep a space of accessibility to all the components. However, when the designer decides to insert his job knowledge in the problem formulation, he automatically influences the search of solutions realised by the algorithm. If he wants to find innovative solutions, he has to simplify the problem description and to adapt the design constraints and the objectives to his design preferences. For example, in the application studied in this paper, the idea is to consider the accessibility to the facilities as a design constraint or an objective. Next research works have to explore this new concept. Actually, this innovative optimization process proposed in this paper also suggests that the method could be improved according to the designer preferences:

• Qualitative fitness could be inserted into design objectives. When solving complex design problems, as a layout problem for example, the translation of some constraints and objectives into simple mathematical expressions can be very difficult. It means that these constraints and objectives could be replaced by a mark given by designers in order to characterize their designs. Then, this qualitative fitness could be considered by the algorithm as a design objective. Alexandra Melike Brintrup has already developed an interactive genetic algorithm-based framework for handling qualitative criteria in design optimization [START_REF] Brintrup | An interactive genetic algorithm-based framework for handling qualitative criteria in design optimization[END_REF]. Actually, it should improve the performances of Pareto-optimal solutions. • The designer could interact with design variables during the optimization process. Stopping the optimizer would allow the designer to firstly analyze a specific solution, secondly locally modify the design configuration and then decide to keep this modified design in the next generation of the genetic algorithm. We can find in [START_REF] Michalek | Interactive design optimization of architectural layouts[END_REF] a significant contribution to this concept applied to the design optimization of architectural layouts.

CONCLUSION

This article presents an innovative layout problem formulation including the concept of space of accessibility defined in section 3. It shows that problem formulation is a very important step in the optimization process because it has a great impact on computed solutions. Secondly, the hybridization of the separation algorithm and the multiobjective algorithm is a very efficient method to ensure a good diversity in a Pareto-optimal solutions set. Moreover, the strategy is designed to allow the interaction between the user and the optimization process in order to improve the performances of Pareto-optimal designs. Actually, for industrial experts, design optimization has great advantages. On the one hand, it allows the designers to explore more alternative solutions to their problem. This is a very good way to encourage the innovation. On the other hand, using design optimization User interaction 2 (to perform manually some solutions) lets the designer to easily make his design choices and justifying them with quantitative values related to his problem formulation.

2

 2

Figure 2 :

 2 Figure 2: Separation algorithm BFGS test. In order to understand the principle of the separation algorithm, let us consider a layout problem test. The dimensions of the square container are

  violation of placement = 0 cm² density = 75% , 142 iterations violation of placement = 0.65 cm² density = 100% , 389 iterations violation of placement = 1 m²

Figure 3 :

 3 Figure 3: Overall view of the shelter.The layout optimization of this shelter is a three dimensional optimization problem. However for the presentation, and fortunately because the cabinets are the full height of the shelter and prevent a superposition of elements, the model is simplified and conceptualized in

Figure 4 :

 4 Figure 4: Configuration model of the shelter in 2D.The formulation of this layout problem is innovative because the components can be classified in two categories: those which have a mass (material components) and those which no have mass (virtual components). Here, the virtual components represent the spaces of accessibility of the cabinets and the desks. For example, the space of accessibility of the cabinet is the required space to insert some materials into the cabinet. These spaces are symbolized in Figure4by dotted rectangles. With this problem formulation, design constraints depend on the category of components. It means that overlap is allowed between two spaces of accessibility, taking account that operations of materials loading are sequentially made whereas overlap has to be minimized between two material components. Moreover, the space, represented by hatching in Figure4, is the space below the air-conditioner where no cabinet can be placed. This space is also a virtual component that is fixed during the optimization process The dimensions of the shelter are 2150mm x 2740mm. The density of this configuration, without considering the spaces of accessibility of the different components, is equal to 50%. If the spaces of accessibility are considered, this density increases up to 90%.

Figure 7 :

 7 Figure 7: Local modification of the solution 3. Table3describes all the objectives values for the initial solution generated by the expert, the solution 3 and the solution 3 locally modified by the designer. Actually, due to the interaction of the designer with the final solutions, a new solution, better than the initial one, is created.Initial solution 3Solution 3 modified

Figure 8 :

 8 Figure 8: Schematic representation of the optimization strategy.

  individuals among P2)Multi-objective optimizer (to generate non-dominated solutions)

Table 1 .

 1 The number of individuals in the initial generation is equal to 240 because a rule of thumb suggests the possibly to accumulate an initial population of at least 16 design configurations and possibly more than 2 x Number of variables x Number of objectives = 2 x 24 x 5 = 240.

	Number of individuals in the initial	240
	generation
	Probability of Directional Cross-Over	0.1
	Probability of Selection	0.05
	Probability of Mutation	0.45
	Probability of Classical Cross-Over	0.4
	Elitism	enabled
	DNA String Mutation Ratio	0.5
	Number max of generations	100

Table 1 :

 1 Multi-Objective Genetic Algorithm parameters.

Table 3 :

 3 Industrial solution vs solution 3.

	O1 (minimize)	25.41 cm	8.67 cm	3.48 cm
	O2 (maximize) 240.58 cm	238.30 cm	240.58 cm
	O3 (maximize) 198.50 cm	177.80 cm	198.50 cm
	O4 (maximize) 165.80 cm	168.29 cm	172.51 cm
	05 (minimize)	0 cm	0 cm	0 cm

Acknowledgement

The authors would like to acknowledge THALES for the application study.