Jean-Charles Faugère
email: jean-charles.faugere@inria.fr

Ayoub Otmani
email: ayoub.otmani@inria.fr

Ludovic Perret
email: ludovic.perret@lip6.fr

Jean-Pierre Tillich
email: jean-pierre.tillich@inria.fr

Algebraic Cryptanalysis of McEliece Variants with Compact Keys

In this paper we propose a new approach to investigate the security of the McEliece cryptosystem. We recall that this cryptosystem relies on the use of error-correcting codes. Since its invention thirty years ago, no efficient attack had been devised that managed to recover the private key. We prove that the private key of the cryptosystem satisfies a system of bi-homogeneous polynomial equations. This property is due to the particular class of codes considered which are alternant codes. We have used these highly structured algebraic equations to mount an efficient key-recovery attack against two recent variants of the McEliece cryptosystems that aim at reducing public key sizes. These two compact variants of McEliece managed to propose keys with less than 20,000 bits. To do so, they proposed to use quasi-cyclic or dyadic structures. An implementation of our algebraic attack in the computer algebra system MAGMA allows to find the secret-key in a negligible time (less than one second) for almost all the proposed challenges. For instance, a private key designed for a 256-bit security has been found in 0.06 seconds with about 2 17.8 operations.

Introduction

Alternative cryptography. Despite the fact that several hard problems have been proposed as a foundation for public-key primitives, those effectively used are essentially classical problems coming from number theory: integer factorization (e.g. in RSA) and discrete logarithm (e.g. in Diffie-Hellman key-exchange). It is well-known that, although polynomial-time algorithms for those problems have not yet been found, they are not safe from a theoretic breakthrough that would endanger the security of the corresponding schemes. For instance, the emergence of a new computer model such as quantum computers would make schemes based on these classical number theory problems totally insecure.

The lack of diversity in public key cryptography has been identified as a major concern in the field of information security.A good illustration of the potential damage of such lack of diversity is hash zoo. The portfolio of hash functions used so far in practice was mainly restricted to the same type of functions which are now almost all broken. Although the American National Institute of Standards and Technology (NIST) issued an international call 4 to design a new standard hash function, the cryptography community will remain in a fuzzy situation until 2012 (date of final decision). One of the main issues in public key cryptography is to identify hard problems that are not based on the classical ones coming from number theory. However, few emerged until now as a viable alternative. As pointed in [START_REF] Avanzi | Lightweight asymmetric cryptography and alternatives to rsa, ecrypt european network of excellence in cryptology[END_REF], promising candidates include: the problem of solving multivariate equations over a finite field, the problem of finding a short vector in a lattice and the problem of decoding a linear code. Those problems known for being NP-hard are not concerned with the quantum computer threat.

McEliece cryptosystem. Among those problems, code-based cryptosystems seem to offer the most promising alternative. McEliece public key cryptosystem [START_REF] Mceliece | A Public-Key System Based on Algebraic Coding Theory[END_REF] is one of the oldest public-key cryptosystems. Its security relies on the difficulty of decoding a linear code. The main advantage of this system is to have very fast encryption and decryption functions. Depending on how the parameters are chosen for a fixed security level, this cryptosystem is about five times faster for encryption and about 10 to 100 times faster for decryption than RSA [START_REF] Biswas | McEliece cryptosystem implementation: Theory and practice[END_REF]. Furthermore, it has withstood many attacking attempts. After more than thirty years now, it still belongs to the very few public key cryptosystems which remain unbroken. Following McEliece's pioneering work, several different public key cryptosystems based on the intractability of decoding a linear code have been proposed [START_REF] Niederreiter | A public-key cryptosystem based on shift register sequences[END_REF][START_REF] Gabidulin | Ideals over a non-commutative ring and their applications to cryptography[END_REF][START_REF] Sidelnikov | A public-key cryptosytem based on Reed-Muller codes[END_REF][START_REF] Janwa | McEliece public key cryptosystems using algebraic-geometric codes[END_REF][START_REF] Berger | Designing an efficient and secure public-key cryptosystem based on reducible rank codes[END_REF][START_REF] Berger | How to mask the structure of codes for a cryptographic use[END_REF][START_REF] Baldi | Cryptanalysis of a new instance of McEliece cryptosystem based on QC-LDPC codes[END_REF][START_REF] Baldi | A new analysis of the McEliece cryptosystem based on QC-LDPC codes[END_REF][START_REF] Berger | Reducing key length of the McEliece cryptosystem[END_REF][START_REF] Misoczki | Compact McEliece keys from Goppa codes[END_REF]. The original McEliece cryptosystem relies on Goppa codes whereas its variants suggested to use different codes. The Sidelnikov system [START_REF] Sidelnikov | A public-key cryptosytem based on Reed-Muller codes[END_REF] used Reed-Muller codes, the Janwa-Moreno system proposed to take algebraic geometric codes [START_REF] Janwa | McEliece public key cryptosystems using algebraic-geometric codes[END_REF] and the Gabidulin-Paramonov-Tretjakov (GPT) cryptosystem considered Gabidulin codes devised for the rank-metric. LDPC codes have also been repeatedly suggested for this use. Niederreiter is the first in [START_REF] Niederreiter | A public-key cryptosystem based on shift register sequences[END_REF] to bring in a significant modification of the McEliece system. However his suggestion to use Generalized Reed-Solomon codes turned out to be an insecure solution [START_REF] Sidelnikov | On the insecurity of cryptosystems based on generalized Reed-Solomon codes[END_REF]. Many of these schemes were broken [START_REF] Sidelnikov | On the insecurity of cryptosystems based on generalized Reed-Solomon codes[END_REF][START_REF] Gibson | Severely denting the Gabidulin version of the McEliece public key cryptosystem[END_REF][START_REF] Minder | Cryptanalysis of the Sidelnikov cryptosystem[END_REF][START_REF] Otmani | Cryptanalysis of McEliece cryptosystem based on quasi-cyclic ldpc codes[END_REF][START_REF] Faure | Cryptanalysis of the McEliece cryptosystem over hyperelliptic curves[END_REF][START_REF] Overbeck | Structural attacks for public key cryptosystems based on Gabidulin codes[END_REF][START_REF] Wieschebrink | Cryptanalysis of the Niederreiter public key scheme based on GRS subcodes[END_REF]. All these attacks result in a total break of the system (the secret key, or an equivalent secret key is recovered from the knowledge of the public key). However, the original McEliece remains unbroken. The fact that the best known attacks are still exponential speaks for itself [START_REF] Stern | A method for finding codewords of small weight[END_REF][START_REF] Bernstein | Attacking and defending the McEliece cryptosystem[END_REF][START_REF] Bernstein | Explicit bounds for generic decoding algorithms for code-based cryptography[END_REF][START_REF] Finiasz | Security bounds for the design of code-based cryptosystems[END_REF]. Despite its impressive resistance against a variety of attacks and its fast encryption and decryption, McEliece cryptosystem has not stood up to RSA for practical applications. This is most likely due to the large size of the public key which is between several hundred thousand and several million bits. To overcome this limitation, a new trend initiated in [START_REF] Gaborit | Shorter keys for code based cryptography[END_REF] manages to decrease the key size by choosing the public matrix defining the code in a particular form; for instance with a quasi-cyclic structure [START_REF] Gaborit | Shorter keys for code based cryptography[END_REF]. This enables to decrease significantly the key size. The same idea was used in [START_REF] Baldi | Cryptanalysis of a new instance of McEliece cryptosystem based on QC-LDPC codes[END_REF] with LDPC codes. Both schemes were broken in [START_REF] Otmani | Cryptanalysis of McEliece cryptosystem based on quasi-cyclic ldpc codes[END_REF]. It should be noted that both proposals did not use the Goppa codes of the McEliece cryptosystem, and the attacks have no impact on its security. This work was then followed by two independent proposals [START_REF] Berger | Reducing key length of the McEliece cryptosystem[END_REF][START_REF] Misoczki | Compact McEliece keys from Goppa codes[END_REF] that are based on the same kind of idea of using quasi-cyclic [START_REF] Berger | Reducing key length of the McEliece cryptosystem[END_REF] or dyadic structure [START_REF] Misoczki | Compact McEliece keys from Goppa codes[END_REF]. Both use the same type of codes called the alternant code family which contains Goppa codes. Actually the codes used in [START_REF] Misoczki | Compact McEliece keys from Goppa codes[END_REF] are Goppa codes. This approach is quite attractive because it results in a drastic improvement of the public key size. In [START_REF] Berger | Reducing key length of the McEliece cryptosystem[END_REF], the size ranges between 8, 000 and 20, 000 bits, whereas it lies between 4, 000 and 20, 000 bits in [START_REF] Misoczki | Compact McEliece keys from Goppa codes[END_REF]. Until now, these new proposals seem to be immune against the attack suggested in [START_REF] Otmani | Cryptanalysis of McEliece cryptosystem based on quasi-cyclic ldpc codes[END_REF].

Our contribution. In this paper we show that both schemes have a serious flaw that can be exploited to recover the private keys. We present an algebraic cryptanalysis 5of the quasi-cyclic and dyadic schemes [START_REF] Berger | Reducing key length of the McEliece cryptosystem[END_REF][START_REF] Misoczki | Compact McEliece keys from Goppa codes[END_REF]. Algebraic cryptanalysis is a general framework that permits to assess the security of theoretically all cryptographic schemes. So far, such type of attacks has been applied successfully against several multivariate schemes and stream ciphers. To our knowledge, it is the first time that such an approach is used against code-based cryptosystems. The basic principle of this cryptanalysis is to associate to a cryptographic primitive a set of algebraic equations. The system of equations is constructed in such a way to have a correspondence between the solutions of this system, and a secret information of the cryptographic primitive (for instance the secret key of an encryption scheme). In McEliece, the algebraic system that we have to solve for recovering the secret-key has the following very specific structure:

g i,0 Y 0 X j 0 + • • • + g i,n-1 Y n-1 X j n-1 = 0 i ∈ {0, . . . , k -1}, j ∈ {0, . . . , r -1} (1)
where the unknowns are the X i 's and the Y i 's and the g i, j 's are known coefficients with 0 ≤ i ≤ k -1, 0 ≤ j ≤ n -1 that belong to a certain field F q with q = 2 s . We look for solutions of this system in a certain extension field F q m . Here k is an integer which is at least equal to nrm. By denoting X def = (X 0 , . . . , X n-1 ) and Y def = (Y 0 , . . . ,Y n-1 ) we will refer to such an algebraic system by McE k,n,r (X, Y). The total number of equations is rk. The number of unknowns 2n and the maximum degree r -1 of the equations can be extremely high when cryptographic parameters are considered (e.g. n = 1024 and r -1 = 49). Thus it is not clear whether an algebraic attack can be mounted efficiently in general. However, in the case of the tweaked McEliece schemes we have either quasi-cyclic or dyadic [START_REF] Berger | Reducing key length of the McEliece cryptosystem[END_REF][START_REF] Misoczki | Compact McEliece keys from Goppa codes[END_REF]. It turns out that is possible to make use of this structure in orderto reduce considerably the number of unknowns in the algebraic system [START_REF] Adams | An Introduction to Gröbner Bases[END_REF]. Moreover, it also appears that by using only the linear equations involving the Y i 's gives a linear space of solutions which is of small dimension. We will explain in Section 4 and Section 5 respectively how to solve the underlying algebraic systems. We will see how the public-key structure (quasi-cyclic or dyadic) impacts on the difficulty of solving the algebraic system (1). In particular, the structure induces an imbalance between the X and Y variables. From a practical point of view, we have been able to recover the secretkey via Gröbner bases computations in a negligible time (less than one second) for most of the parameters proposed in [START_REF] Berger | Reducing key length of the McEliece cryptosystem[END_REF][START_REF] Misoczki | Compact McEliece keys from Goppa codes[END_REF]. Before that, we briefly recall in the next section the McEliece scheme and we explain in Section 3 how we derive the algebraic system (1).

McEliece Public-Key Cryptosystem

We recall here how the McEliece public-key cryptosystem is defined.

Secret key: the triplet (S, G s , P) of matrices defined over a finite field F q over q elements, with q being a power of two, that is q = 2 s . G s is a full rank matrix of size k × n, with k < n, S is of size k × k and is invertible, and P is permutation matrix of size n × n. Moreover G s defines a code (which is the set of all possible uG s with u ranging over F k q ) which has a decoding algorithm which can correct in polynomial time a set of errors of weight at most t. This means that it can recover in polynomial time u from the knowledge of uG s + e for all possible e ∈ F n q of Hamming weight at most t. Public key: the matrix product G = SG s P.

Encryption: A plaintext u ∈ F k q is encrypted by choosing a random vector e in F n q of weight at most t. The corresponding ciphertext is c = uG + e. Decryption: c ′ = cP -1 is computed from the ciphertext c. Notice that c ′ = (uSG s P + e)P -1 = uSG s + eP -1 and that eP -1 is of Hamming weight at most t. Therefore the aforementioned decoding algorithm can recover in polynomial time uS. This vector is multiplied by S -1 to obtain the plaintext u.

This describes the general scheme suggested by McEliece. From now on, we say that G is the public generator matrix and the vector space C spanned by its rows is the public code i.e.

C def = uG | u ∈ F k q .
What is generally referred to as the McEliece cryptosystem is this scheme together with a particular choice of the code, which consists in taking a binary Goppa code. This class of codes belongs to a more general class of codes, namely the alternant code family ([24, Chap. 12, p. 365]). The main feature of this last class of codes is the fact that they can be decoded in polynomial time.

Algebraic Approach

Setting up the algebraic system. In this part, we explain more precisely how we construct the algebraic system described in [START_REF] Adams | An Introduction to Gröbner Bases[END_REF]. As explained in the previous section, the McEliece cryptosystem relies on Goppa codes which belong to the class of alternant codes and inherit from this an efficient decoding algorithm. We will describe this class of codes in more details since both cryptosystems that we cryptanalyze use such codes. It is convenient to describe them through a parity-check matrix. This is an r × n matrix H defined over an extension F q m of the field over which the code is defined, which is such that

{uG s | u ∈ F k q } = {c ∈ F n q | Hc T = 0}. (2)
r satisfies in this case the condition r ≥ n-k m . In the case of alternant codes, there exists a parity-check matrix with a very special form related to Vandermonde matrices. More precisely there exist two vectors x = (x 0 , . . . , x n-1 ) and y = (y 0 , . . . , y n-1 ) in F n q m such that V r (x, y) is a parity-check matrix , with

V r (x, y) def =      y 0 • • • y n-1 y 0 x 0 • • • y n-1 x n-1 . . . . . . y 0 x r-1 0 • • • y n-1 x r-1 n-1      . (3)
We use the following notation in what follows Definition 1. The alternant code A r (x, y) of order r over F q associated to x = (x 0 , . . . , x n-1 ) where the x i 's are different elements of F q m and y = (y 0 , . . . , y n-1 ) where the y i 's are nonzero elements of F q m is defined by

A r (x, y) = {c ∈ F n q | V r (x, y)c T = 0}.
It should be noted that the public code in the McEliece scheme is also an alternant code.

We denote here by the public code, the set of vectors of the form

{uG | u ∈ F k q } = {cSG s P | c ∈ F k q }.
This is simple consequence of the fact that the set

{uSG s P | u ∈ F k q } is obtained from the secret code {uG s | u ∈ F k q }
by permuting coordinates in it with the help of P, since multiplying by an invertible matrix S of size k × k leaves the code globally invariant.

The key feature of an alternant code is the following fact Fact 1. There exists a polynomial time algorithm decoding an alternant code once a parity-check matrix H of the form H = V r (x, y) is given.

In other words, it is possible to break the McEliece scheme, if it is possible to find x * and y * in F n q m such that

{xG | x ∈ F n q } = {y ∈ F n q | V r (x * , y * )y T = 0}. (4)
From the knowledge of this matrix V r (x * , y * ), it is possible to decode the public code, that is to say to recover u from uG + e. Finding such a matrix clearly amounts to find a matrix V r (x * , y * ) such that V r (x * , y * )G T = 0. By bringing in 2n variables X 0 , . . . , X n-1 and Y 0 , . . . ,Y n-1 where X i corresponds to x * i and Y i to y * i we see that this is equivalent to solve the following system:

g i,0 Y 0 X j 0 + • • • + g i,n-1 Y n-1 X j n-1 = 0 i ∈ {0, . . . , k -1}, j ∈ {0, . . . , r -1} (5)
where the g i, j 's are the entries of the known matrix G with 0

≤ i ≤ k -1 and 0 ≤ j ≤ r -1.
The cryptosystems proposed in [START_REF] Berger | Reducing key length of the McEliece cryptosystem[END_REF][START_REF] Misoczki | Compact McEliece keys from Goppa codes[END_REF] follow the McEliece scheme [START_REF] Mceliece | A Public-Key System Based on Algebraic Coding Theory[END_REF] with the additional goal to design a public-key cryptosystem with very small key sizes. They both require to identify alternant codes having a property that allows matrices to be represented by very few rows. In the case of [START_REF] Berger | Reducing key length of the McEliece cryptosystem[END_REF] circulant matrices are chosen whereas the scheme [START_REF] Misoczki | Compact McEliece keys from Goppa codes[END_REF] focuses on dyadic matrices. These two families have in common the fact the matrices are completely described from the first row. The public generator matrix G in these schemes is a block matrix where each block is circulant in [START_REF] Berger | Reducing key length of the McEliece cryptosystem[END_REF] and dyadic in [START_REF] Misoczki | Compact McEliece keys from Goppa codes[END_REF].

We shall see that the algebraic approach previously described leads to a key-recovery in nearly all the parameters proposed in both schemes. The crucial point that makes the attack possible is the fact that we have a system with much less unknowns than in the case of the McEliece cryptosystem. This is due to both the particular structure of the matrices and their block form that describe the public alternant codes. We finally end this section with a simple remark which explains that we can basically set one of the Y i 's and two values of the X i 's to an arbitrary value in the algebraic system [START_REF] Adams | An Introduction to Gröbner Bases[END_REF].

Proposition 1 ([24, Chap. 10, p. 305]). Let (X 0 , . . . , X n-1 ), (Y 0 , . . . ,Y n-1 ) be a solution of ( 1) and a = 0, b,c = 0 be elements of F q m . Then (aX 0 + b, . . . , aX n-1 + b) and (cY 0 , . . . , cY n-1 ) is also a solution of [START_REF] Adams | An Introduction to Gröbner Bases[END_REF].

Solving the Algebraic System. In this part, we describe a general technique to solve the algebraic system McE k,n,r (X, Y) using Gröbner bases techniques [START_REF] Buchberger | Ein Algorithmus zum Auffinden der Basiselemente des Restklassenringes nach einem nulldimensionalen Polynomideal[END_REF][START_REF] Cox | Ideals, Varieties, and algorithms: an Introduction to Computational Algebraic Geometry and Commutative Algebra[END_REF][START_REF] Faugère | A new efficient algorithm for computing gröbner bases (f4)[END_REF][START_REF] Faugère | A new efficient algorithm for computing gröbner bases without reduction to zero : F5[END_REF]. Although the particular characteristics of the cryptosystems [START_REF] Berger | Reducing key length of the McEliece cryptosystem[END_REF][START_REF] Misoczki | Compact McEliece keys from Goppa codes[END_REF] studied here will further influence the shape of McE k,n,r (X, Y) (number of variables, number of equations, . . .), we have designed a special strategy for taking advantage as much as possible of the intrinsic structure. We have made an implementation of the strategy described here. We will present the experimental results, as well as the improvements which are possible due to the quasicyclic and dyadic structures, in Section 4 (quasi-cyclic case) and Section 5 (dyadic case). As a first general remark, it is readily seen that McE k,n,r (X, Y) is highly structured. For instance, it is very sparse as the only monomials occurring in the system are of the form Y i X j i , with 0 ≤ i ≤ k -1 and 0 ≤ j ≤ r -1. It can also be noticed that each block of k equations is bi-homogeneous, i.e. homogeneous if the variables of X (resp. Y) are considered alone. Note that such structure already appears in the cryptanalysis of the MinRank problem [START_REF] Faugère | Cryptanalysis of minrank[END_REF]. Due to the particular structure of the system, it makes sense to design a specific strategy for solving McE k,n,r (X, Y). A simple way for solving this system would consist in generating the equations and try to solve it directly with a generic algorithm (for instance, the Gröbner basis algorithm available in the MAGMA computer algebra software). This approach fails for most challenges proposed in [START_REF] Misoczki | Compact McEliece keys from Goppa codes[END_REF]. However, it is interesting to remark that this direct approach has been successful in practice for all challenges of [START_REF] Berger | Reducing key length of the McEliece cryptosystem[END_REF]. We only mention that it takes between few minutes to 24 hours of computation using a negligible amount of memory. As a comparaison, the improved strategy that we will describe now permits to solve (almost) all the challenges of [START_REF] Berger | Reducing key length of the McEliece cryptosystem[END_REF][START_REF] Misoczki | Compact McEliece keys from Goppa codes[END_REF] in few seconds using also a negligible amount of memory.

The first fundamental remark is that there are k linear equations in the n variables of the block Y in McE k,n,r (X, Y). This implies that all the variables of the block Y can be expressed in terms of d ≥ nk variables. From now on, we will always assume that the variables of the block Y ′ only refer to these d free variables. The first step is then to rewrite the system (1) only in function of the variables of X and Y ′ , i.e. the variables of Y \ Y ′ are substituted by linear combinations involving only variables of Y ′ . For the cryptosystems considered in this paper [START_REF] Berger | Reducing key length of the McEliece cryptosystem[END_REF][START_REF] Misoczki | Compact McEliece keys from Goppa codes[END_REF], the number of free variables d in Y ′ can be rather small (typically 1 or 2 for some challenges). Furthermore, the quasi-cyclic and dyadic structures provide additional linear equations in the variables of X and Y ′ which can be also used to rewrite/clean the system. In the sequel, we will denote by r (X,Y) by removing all the linear equations in X and Y.

McE k,n,r (X ′ , Y ′ ) the system obtained from McE k,n,
The second crucial point is that as soon as the the projection of the solutions on the variables Y ′ are known, the system (1) simplifies to:

g ′ i,0 X j 0 + • • • + g ′ i,n-1 X j n-1 = 0 i ∈ {0, . . . , k -1}, j ∈ {0, . . . , r -1} .
This system is readily solved by keeping only the equations in this system which correspond to powers of the X i 's which are powers of two. In other words we consider only the equations of the form g

′ i,0 X 2 j 0 + • • • + g ′ i,n-1 X 2 j
n-1 = 0 for j in {0, . . . , ⌊log 2 (r -1)⌋} and i in {0, . . . , k -1}. Hence, we obtain a quasi bi-linear system because the system is always defined over a field of characteristic two. Moreover, it has very few monomials per equation and can be easily solved in practice by computing a Gröbner basis.

The most difficult part of the computation is to find a projection of the solutions with respect to the variables of the block Y ′ . Notice that an exhaustive search on the d free variables of Y ′ leads to a practical attack for some of the challenges proposed in [START_REF] Berger | Reducing key length of the McEliece cryptosystem[END_REF][START_REF] Misoczki | Compact McEliece keys from Goppa codes[END_REF]. We will present below an even more efficient strategy to recover Y ′ . More formally, let I be the ideal generated by McE k,n,r (X ′ , Y ′ ) and V be the corresponding variety i.e. the set of solutions. The goal is to compute the projection of V , denoted by V ′ , on the variables of Y ′ . This can done by computing a Gröbner basis (w.r.t. a degree order) G deg of I ∩F q m [Y ′ ]. This is a classical problem in computer algebra which can be solved by using standard elimination techniques (for instance see [1, Chap. 2.3, p. 69] or [START_REF] Cox | Ideals, Varieties, and algorithms: an Introduction to Computational Algebraic Geometry and Commutative Algebra[END_REF]Chap. 3,p. 112]. In the appendix, we briefly recall basic facts about elimination theory. In our context, we have used a slightly modified version of F 4 [START_REF] Faugère | A new efficient algorithm for computing gröbner bases (f4)[END_REF] for computing a Gröbner basis G deg of I ∩ F q m [Y ′ ]. Roughly speaking, the idea is to adapt the algorithm for performing the Gröbner basis computation in

F q m [X ′ ][Y ′ ],
i.e. the set of polynomials in Y ′ whose coefficients are polynomials in F q m [X ′ ]. As for the usual F 4 , we process degree by degree. However, we consider only the degree of the polynomials w.r.t. the variables of X ′ . We stop the computation as soon as we have sufficiently many equations in Y ′ , in other words, as soon as we detect that V ′ has a finite number of solution. Below, we describe more precisely the procedure which allows to compute a Gröbner basis

G deg of I ∩ F q m [Y ′ ].
As already explained, this is the most difficult part. Furthermore, to be sure that the variety V ′ associated to I ∩ F q m [Y ′ ] has few solutions, we have to remove parasite solutions corresponding to X i = X j or to Y j = 0. A classical way to do that is to introduce new variables u i j and v i and add to

McE k,n,r (X ′ , Y ′ ) Algorithm 1: ComputeProjection Input : The system McE k,n,r (X ′ , Y ′ ) Output: A Gröbner basis G deg of I ∩ F q m [Y ′ ] Let F be the equations of degree 2 i , 1 ≤ i ≤ r -1 of McE k,n,r (X ′ , Y ′ )
Let F ′ be the system obtained from F by fixing "randomly" some variables of X ′ Compute a Gröbner basis G deg of I ∩ F q m [Y ′ ] using the tweaked version of F 4 Return G deg equations of the form:

u i j • (X i -X j ) + 1 = and v i •Y i + 1 = 0.
In practice, we have not added all theses equations; but only few of them namely 4 or 5. The reason is that we do not want to add too many new variables. In addition, including few of such equations already permits to remove trivial solutions. We also have to remove some degree of freedom in the algebraic system by fixing randomly few variables of X' as explained in Proposition 1. It is important to notice that since we are removing component of high dimension the new system is indeed much faster to solve. Finally, we have not considered all the equations of McE k,n,r (X ′ , Y ′ ) to compute G deg . This system being naturally over-defined, we can "safely" remove some equations. Typically, it makes sense to consider the smaller subset of equations such that V ′ is zerodimensional and for which we can efficiently compute G deg . The variety V ′ having few elements it is not difficult to recover this set from the knowledge of G deg .

Algebraic Cryptanalysis of the Quasi-Cyclic Variant

The scheme presented in [START_REF] Berger | Reducing key length of the McEliece cryptosystem[END_REF] suggests to use block matrices where each block is a circulant matrix. The public code C suggested in [START_REF] Berger | Reducing key length of the McEliece cryptosystem[END_REF] is defined on a field F q = F 2 s which is considered as a subfield of F q m for a certain integer m. Let α be a primitive element of F q m . Let ℓ and N 0 be such that q m -1 = N 0 ℓ and let β be an element of F q m of order ℓ. Although this is not explicitly stated in [START_REF] Berger | Reducing key length of the McEliece cryptosystem[END_REF], it is readily checked that C is defined from an r ×n parity-check matrix H over F q m which is the juxtaposition of n 0 (n = ℓn 0 ) matrices

H (0) , • • • , H (n 0 -1) of size r × ℓ. Each H (b) = (h (b) i, j ) with 0 ≤ b ≤ n 0 -1, 0 ≤ i ≤ r -1 and 0 ≤ j ≤ ℓ -1 is given by h (b) i, j = γ b β (d b + j)e α w b β d b + j i (6)
where γ b is a nonzero element of F q m , d b is an integer of {0, . . . , ℓ -1}, e is an integer of {0, . . . , ℓ -1} and the w b 's are distinct integers of {0, . . . , N 0 -1}. From this, it is now clear that C is an alternant code A r (x, y) of order r associated to x = (x 0 , . . . , x n-1 ) and y = (y 0 , . . . , y n-1 ) which satisfy for any j in {0, . . . , ℓ -1}

x bℓ+ j = α w b β d b + j (7)

y bℓ+ j = γ b β (d b + j)e , (8)
It can be checked (see [START_REF] Berger | Reducing key length of the McEliece cryptosystem[END_REF]) that C has a public generating matrix G which is block circulant of size k × n with k of the form k = k 0 ℓ for some integer k 0 (recall that k ≥ nrm).

We present now an algebraic attack against the quasi-cyclic variant proposed in [START_REF] Berger | Reducing key length of the McEliece cryptosystem[END_REF] that recovers x and y by setting up an algebraic system of the form McE k,n,r (X, Y) from the equation HG T = 0. This would also give a system with 2n unknowns. We can obtain a huge reduction of the number of unknowns by using Equations ( 7) and ( 8) which induce some linear relations between the x i 's and the y i 's. From these two equations we deduce that x bℓ+ j = x bℓ β j (9)

y bℓ+ j = y bℓ β je , (10)
for any j in {0, . . . , ℓ -1} and j in {0, . . . , n 0 -1}. Furthermore, since in the cases considered in [START_REF] Berger | Reducing key length of the McEliece cryptosystem[END_REF], e is small because it lies in the range {0, . . . , ℓ -1} and ℓ is less than 100, we may assume that:

Assumption 2. The secret integer e such that 0 ≤ e ≤ ℓ -1 is known.

This enables to simplify the description of the system McE k,n,r (X, Y). By setting up the unknown X b for x bℓ and the unknown Y b for y bℓ we obtain the following algebraic system.

Proposition 2. Let G = (g i, j ) be the k × n public generator matrix with k = k 0 ℓ and n = n 0 ℓ. For any 0 ≤ w ≤ r -1 and any 0 ≤ i ≤ k -1, the unknowns X 0 , . . . , X n 0 -1 and Y 0 , . . . ,Y n 0 -1 should satisfy:

n 0 -1 ∑ b=0 g ′ i,b,w Y b X w b = 0 where g ′ i,b,w def = ℓ-1 ∑ j=0 g i,bℓ+ j β j(e+w) . ( 11
)
Proof. We observe that

n-1 ∑ j=0 g i, j y j x w j = n 0 -1 ∑ b=0 ℓ-1 ∑ j=0 g i,bℓ+ j y bℓ+ j x w bℓ+ j = n 0 -1 ∑ b=0 ℓ-1 ∑ j=0 g i,bℓ+ j y bℓ x w bℓ β je β jw = n 0 -1 ∑ b=0 y bℓ x w bℓ ℓ-1 ∑ j=0 g i,bℓ+ j β je β jw
By setting X b for x bℓ and Y b for y bℓ we obtain the aforementioned system.

Theoretically by Proposition 1, we would be able to fix two variables, say X 0 and X 1 , and one variable Y j , for instance Y 0 , to arbitrary values as long as X 0 = X 1 and Y 0 = 0. However, if we do it, we then lose the linear relations between the x i 's given in [START_REF] Bernstein | Explicit bounds for generic decoding algorithms for code-based cryptography[END_REF]. Therefore we can only fix one X i and one Y i as stated in Lemma 1 that is straightforward to prove. Lemma 1. Let (X 0 , . . . , X n 0 -1 ), (Y 0 , . . . ,Y n 0 -1 ) be a solution of [START_REF] Buchberger | Ein Algorithmus zum Auffinden der Basiselemente des Restklassenringes nach einem nulldimensionalen Polynomideal[END_REF]. Then (aX 0 , . . . , aX n 0 -1 ) and (cY 0 , . . . , cY n 0 -1 ) is also a solution of [START_REF] Buchberger | Ein Algorithmus zum Auffinden der Basiselemente des Restklassenringes nach einem nulldimensionalen Polynomideal[END_REF] for any a = 0 and c = 0 of F q m . Hence, the total number of unknowns is 2(n 0 -1). Furthermore there are many redundant equations in Proposition 2. This comes from the block circulant form of G. From this form we know that g iℓ+u,bℓ+ j = g iℓ,bℓ+(( j-u) mod ℓ) for all u in {0, . . . , ℓ -1} and i in {0, . . . , k 0 -1}. We also have:

g ′ iℓ+u,b,w = ℓ-1 ∑ j=0 g iℓ+u,bℓ+ j β j(e+w) = ℓ-1 ∑ j=0 g iℓ,bℓ+(( j-u) mod ℓ) β j(e+w) = ℓ-1 ∑ j=0 g iℓ,bℓ+ j β j(e+w) β u(e+w) = g ′ iℓ,b,w β u(e+w)
We used the fact β ℓ(e+w) = 1. So for a given i, when u describes {0, . . . , ℓ -1}, the equations

∑ n 0 -1 b=0 g ′ iℓ+u,b,w Y b X w b = 0 are all equivalent.
This means that instead of having rk equations we have only rk ℓ = k 0 r algebraic equations. Proposition 3. The system [START_REF] Buchberger | Ein Algorithmus zum Auffinden der Basiselemente des Restklassenringes nach einem nulldimensionalen Polynomideal[END_REF] has (n 0 -1) unknowns Y i and (n 0 -1) unknowns X i . Furthermore, it has k 0 linear equations involving only the Y i 's and (r -1)k/ℓ = (r -1)k 0 polynomial equations involving the unknowns Y i X w i with w > 1. From now on, we will always assume that redundant equations have been removed and the variables X 0 and Y 0 are fixed. Finally, note that there are

d def = n 0 -1 -k 0 free variables for the Y i 's.
We now present experimental results obtained when solving the system described in (11) using the strategy described in Section 3. The experimental results have been obtained with several Xeon bi-processor 3.2 Ghz, with 16 Gb of Ram. The instances of our problem have been generated using the MAGMA software. We used the MAGMA version 2.15 for our computations. The F 5 [START_REF] Faugère | A new efficient algorithm for computing gröbner bases without reduction to zero : F5[END_REF] algorithm has been implemented in C language in the FGb software. We used this implementation for computing the first Gröbner basis (i.e. which is used in Algorithm 1). All the other computations are performed under MAGMA including factorizing some univariate polynomials and computing Gröbner bases using the F 4 [START_REF] Faugère | A new efficient algorithm for computing gröbner bases (f4)[END_REF] algorithm. The most important observation is that we have been able to solve all the challenges of [START_REF] Berger | Reducing key length of the McEliece cryptosystem[END_REF] in a negligible time because the dimension d = n 0 -1k 0 of the vector space solution for the Y i 's is very small. We also proposed a challenge QC 600 for showing the behaviour of our attack for high security levels.

Algebraic Cryptanalysis of the Dyadic Variant

The cryptosystem presented in [START_REF] Misoczki | Compact McEliece keys from Goppa codes[END_REF] considers particular alternant codes called quasidyadic Goppa codes. Goppa codes form an important subclass of alternant codes. Goppa codes are defined by means of a polynomial G(X) of degree ℓ with coefficients in F q m and for which the sequence x is assumed not to contain any root of G(X). The alternant code defined by the parity-check matrix V ℓ (x, y) with y i = G(x i ) -1 is called a Goppa code over F q and is denoted by G (x, G). A detailed description of the key generation is given in Appendix B. We only provide important facts that are useful for recovering the private key. We first state an important result that shows that G defines actually an alternant code. The proof is given in Appendix C. The last important fact to know about G is that it is a k × n matrix over F q such that n = n 0 ℓ and k = nmℓ where n 0 , ℓ are given integers.

Proposition 4. The code defined by the public generator matrix G is an alternant code

A ℓ (x, y) where for any 0 ≤ j ≤ n 0 -1 and 0 ≤ i, i ′ ≤ ℓ -1, we have the following equations: [START_REF] Cox | Ideals, Varieties, and algorithms: an Introduction to Computational Algebraic Geometry and Commutative Algebra[END_REF] where ⊕ is the bitwise exclusive-or on the binary representation of the indices.

   y jℓ+i = y jℓ x jℓ+i + x jℓ = x i + x 0 x jℓ+(i⊕i ′ ) = x jℓ+i + x jℓ+i ′ + x jℓ
The cryptanalysis of the system consists in defining n 0 unknowns Y 0 , . . . ,Y n 0 -1 that play the role of the y j 's and n unknowns X 0 , . . . , X n that represent the x j 's. We know specify the system of equations that we obtain directly from Proposition 4. Proposition 5. For any w, j, i and i ′ such that 0 ≤ w ≤ ℓ -1, 0 ≤ j ≤ n 0 -1 and 1 ≤ i, i ′ ≤ ℓ -1, we have:

                 n 0 -1 ∑ j=0 Y j ℓ-1 ∑ l=0 g i, jℓ+l X w jℓ+l = 0 X jℓ+i + X jℓ + X i + X 0 = 0 X jℓ+(i⊕i ′ ) + X jℓ+i + X jℓ+i ′ + X jℓ = 0 (13)
It is possible to simplify System (13) by observing, thanks to the third equation, that actually many variables X i 's can be expressed in function of some few variables, namely X 2 j with 0 ≤ j ≤ log 2 (ℓ -1) and X b with 0 ≤ b ≤ n 0 -1.

Corollary 1.

For any

1 ≤ i ≤ ℓ-1, if we write the binary decomposition of i = ∑ log 2 (ℓ-1) j=0 η j 2 j
then the following equation holds:

X i = X 0 + log 2 (ℓ-1) ∑ j=0 η j (X 2 j + X 0 ).
We are also able to provide the exact number of unknowns we can fix to arbitrary values.

Lemma 2. Let (X 0 , . . . , X n-1 ), (Y 0 , . . . ,Y n-1 ) be a solution of ( 13) and a = 0, b,c = 0 be elements of F q m . Then (aX 0 + b, . . . , aX n-1 + b) and (cY 0 , . . . , cY n-1 ) is also a solution of [START_REF] Faugère | A new efficient algorithm for computing gröbner bases (f4)[END_REF].

Proof. The only fact to prove is that (X 0 + b, . . . , X n-1 + b) is also a solution of the last two equations in [START_REF] Faugère | A new efficient algorithm for computing gröbner bases (f4)[END_REF]. It is readily checked since F q m is of charateristic two.

We can now completely give the effective number of equations after elimination of redundant equations. Proposition 6. The system [START_REF] Faugère | A new efficient algorithm for computing gröbner bases (f4)[END_REF] has n 0 -1 unknowns Y i and n 0 -2 + log 2 (ℓ) unknowns X i . Furthermore, it has n 0 -m linear equations involving only the Y i 's, and (ℓ-1)ℓ(n 0m) polynomial equations involving the unknowns Y i X w i with w > 1.

Proof. The number of variables Y j is (n 0 -1) since we can choose Y 0 = 1. As for the variables X j , we observe that they can all be expressed only in function of X 2 j and X iℓ with 0 ≤ j ≤ log 2 (ℓ -1) and 0 ≤ i ≤ n 0 -1. So the number of unknowns X j is log 2 (ℓ -1) + 1 + n 0 -2 since we can fix two different arbitrary values for two variables, say X 0 and X ℓ (Lemma 2). Using the fact that log 2 (ℓ-1) = log 2 (ℓ)-1 since ℓ is a power of 2, we get the claimed number of unknowns. Furthermore, because of the dyadicity of G, the equations obtained with w = 0 are identical when g describes all the rows of a dyadic block of G. This does not appear when w > 1. So we have k/ℓ = n 0m linear equations that involve the Y i 's and (ℓ -1)k = (ℓ -1)ℓ(n 0m) polynomial equations that contain variables of the form Y i X w i where w > 1.

We now present in Table 2 the experimental results we obtained when we solve the system described in (13) using the strategy described in Section 3. As previously, the experimental results have been obtained with several Xeon bi-processor 3.2Ghz, with 16 Gb of Ram. The instances of our problem have been generated using the MAGMA software. We used the MAGMA version 2.15 for our computations. The F 5 [START_REF] Faugère | A new efficient algorithm for computing gröbner bases without reduction to zero : F5[END_REF] algorithm has been implemented in C language in the FGb software. We used this implementation for computing the first Gröbner basis (i.e. which is used in Algorithm 1). All the other computations are performed under MAGMA including factorizing some univariate polynomials and computing Gröbner bases using the F 4 algorithm [START_REF] Faugère | A new efficient algorithm for computing gröbner bases (f4)[END_REF]. Table 2 also shows the impact of the degree extension m. Indeed, computation times indicate that the solutions are easy to find if m is small. This phenomenon is directly related to the size of the solution space of the variables Y i . We have seen in Section 3 that such variables satisfy a system of linear equations. From Proposition 6, the number of linear equations is k 0 = n 0m whereas the number of unknowns Y i is n 0 -1. Thus the dimension of the vector space solution for the Y i 's is m -1. We recall that once the Y i are solved, we find the solutions for the X i 's by solving a system of linear equations. We also give in Table 2 new parameters (Dyadic 256 and Dyadic 512 ) that are generated by "extrapolating" challenges given in [START_REF] Misoczki | Compact McEliece keys from Goppa codes[END_REF]. We observe that this solution space is manageable in practise as long as m < 16 because we did not succeed to find an efficient way to solve the challenges of [START_REF] Misoczki | Compact McEliece keys from Goppa codes[END_REF] when m = 16.

Conclusion

We described in this paper a new algebraic approach to assess the security of the McEliece cryptosystem. We showed that the private key of this scheme is a solution of a very structured system of bi-homogeneous polynomial equations in two sets of unknowns Y i and X i . The solutions belong to a finite field F q m whereas the coefficients of the system are in F q for some known integers m and q. This system comes from the particular structure of alternant codes used in McEliece. Indeed, the Goppa codes as proposed in [START_REF] Mceliece | A Public-Key System Based on Algebraic Coding Theory[END_REF] form a subfamily of alternant codes. Furthermore, the system is composed of two parts of equations: one part that consists of linear equations that involve only the unknowns Y i and a second part where the equations involve terms of the form Y i X j i .

We applied this approach to two new cryptosystems [START_REF] Misoczki | Compact McEliece keys from Goppa codes[END_REF][START_REF] Berger | Reducing key length of the McEliece cryptosystem[END_REF] that are variants of the McEliece scheme. Both aim at reducing the public keys by using very structured block matrices (cyclic matrices in [START_REF] Berger | Reducing key length of the McEliece cryptosystem[END_REF] and dyadic matrices in [START_REF] Misoczki | Compact McEliece keys from Goppa codes[END_REF]). We show that our new cryptanalytic point of view is very efficient for all the parameters proposed in [START_REF] Berger | Reducing key length of the McEliece cryptosystem[END_REF]. An implementation in MAGMA validates our attack and shows that the private key can be found in a negligible time. For the scheme [START_REF] Misoczki | Compact McEliece keys from Goppa codes[END_REF], we are also able to fully recover the private key in almost all cases. An implementation in MAGMA shows that this can be done in time comparable to [START_REF] Berger | Reducing key length of the McEliece cryptosystem[END_REF] as long as the dimension m of the solution vector space of the Y i 's is small.

Thanks to a very recent development [START_REF] Faugère | Gröbner bases of bihomogeneous ideals generated by polynomials of bidegree (1,1): Algorithms and complexity[END_REF] on the solving of bihomogeneous bilinear systems, it is very likely that the solving technique presented here can be replaced by a new version of F 5 dedicated to bi-linear systems. In our case, we can obtain a (quasi) bilinear system when consider the equations involving terms of the form Y i X 2 j i . Moreover, this will permit to precisely estimate the complexity of the attack presented in this paper and will provide a concrete criteria to evaluate the security of future compact McEliece's variants.

Finally, it would be interesting to study if this algebraic approach followed here can be improved in order to mount an attack on the original McEliece cryptosystem. In this case however, there are much more unknowns than in the cases considered here and there is much more freedom left on the Y i 's by looking at the linear equations involving only the Y i 's.

A Gröbner Basics.

The classical approach for computing a Gröbner basis of I ∩ F q m [Y ′ ] can be described as follows. A reader already familiar with polynomial system solving can skip this part. We have to choose a suitable ordering on the monomials (for a definition of such orders, see for instance [12, Chap. 2, p. 52]). In particular, we have to select an elimination ordering ([1, Chap. 2.3, p. 69]) on the blocks X ′ , Y ′ such that the variables occurring in X are greater that those of Y ′ (denoted by X ′ >> Y ′ ). According to [1, Theo. 2.3.4, Chap. 2.3, p. 69], this elimination ordering will permit to compute a Gröbner basis G deg of I ∩ F q m [Y ′ ] with respect to a degree order on the variables of Y ′ (i.e. this is the order induced when "removing" the variables of the block X ′ in the elimination ordering). In theory, to compute the variety V ′ associated to I ∩ F q m [Y ′ ], we have to perform a change of ordering on G deg to compute Gröbner basis G lex of I ∩F q m [Y ′ ]. If we assume that V ′ is zero-dimensional (i.e. has a finite number of solutions so that #V ′ < ∞), then an efficient tool to perform the change of ordering is the FGLM algorithm [START_REF] Faugère | Efficient computation of zero-dimensional gröbner bases by change of ordering[END_REF]. The complexity of computing G lex from G deg with FGLM is polynomial in the size of V ′ , i.e. O (#V ′ ) 3 . In our case, the size of V ′ is very small (< 10).

We have used a slightly modified version of F 4 [START_REF] Faugère | A new efficient algorithm for computing gröbner bases (f4)[END_REF] for computing a Gröbner basis

G deg of I ∩ F q m [Y ′ ].
The idea is to adapt the algorithm for performing the Gröbner basis computation in F q m [X ′ ][Y ′ ], i.e. the set of polynomials in Y ′ whose coefficients are polynomials in F q m [X ′ ]. As for the usual F 4 , we process degree by degree. However, we consider only the degree of the polynomials w.r.t. the variables of X ′ . We stop the computation as soon as we have sufficiently many equations in Y ′ (for instance, as soon as we detect that V ′ has a finite number of solution, i.e. of dimension zero ). The modified version is defined below.

(P d ), G) for h ∈ R do P := P ∪ {CritPair(h, g) | g ∈ G} G := G ∪ {h} return G ∩ F q m [Y ′ ]
For the definition of MATRIX REDUCTION, and CritPair, we refer to [START_REF] Faugère | A new efficient algorithm for computing gröbner bases (f4)[END_REF]. Briefly, the first function performs the usual polynomial reduction of Buchberger's algorithm [START_REF] Cox | Ideals, Varieties, and algorithms: an Introduction to Computational Algebraic Geometry and Commutative Algebra[END_REF] using linear algebra. The second function selects critical pairs with respect to a defined strategy.

B Description of the Variant based on Dyadic Goppa Codes

The cryptosystem presented in [START_REF] Misoczki | Compact McEliece keys from Goppa codes[END_REF] considers particular alternant codes called quasidyadic Goppa codes. Goppa codes form an important subclass of alternant codes. Goppa codes are defined by means of a polynomial G(X) of degree ℓ with coefficients in F q m and for which the sequence x is assumed not to contain any root of G(X). The alternant code defined by the parity-check matrix V ℓ (x, y) with y i = G(x i ) -1 is called a Goppa code over F q and is denoted by G (x, G). It has dimension nmℓ and minimum distance d ≥ ℓ + 1 [24, Chap. 12, p. 340]. In the special case where the roots z = (z 0 , . . . , z ℓ-1 ) of G(X) are distinct and all belong to F q m then G (x, G) admits a parity-check matrix C(z, x) in Cauchy form [24, p. 345]. The scheme in [START_REF] Misoczki | Compact McEliece keys from Goppa codes[END_REF] considers a Goppa code that admits a parity-check matrix that is both a Cauchy matrix and a block matrix where each block is dyadic. An ℓ × ℓ matrix ∆ = (∆ i, j ) with 0 ≤ i ≤ ℓ -1 and 0 ≤ j ≤ ℓ -1 is dyadic if and only if ∆ i, j = h i⊕ j where ⊕ is the bitwise exclusive-or on the binary representation of the indices and h = (h 0 , . . . , h ℓ-1 ) is the first row of ∆ . Let h = (h 0 , . . . , h N-1 ) be a vector of F N q m with ℓ ≤ N. We denote by ∆ ℓ (h) = (∆ i, j ) the ℓ × N matrix such that ∆ i, j = h i⊕ j . One can easily observe that ∆ ℓ (h) is the juxtaposition of N 0 dyadic matrices of size ℓ × ℓ when N = N 0 ℓ for some integer N 0 . Proposition 7 proved in [27, Theorem 2] characterizes dyadic Cauchy matrices. Proposition 7. A necessary and sufficient condition for ∆ ℓ (h) to be a Cauchy matrix C(z, x) is that F q m is of characteristic 2 and for any i, j in {0, . . . , N -1} we have:

1 h i⊕ j = 1 h j + 1 h i + 1 h 0 . ( 14
)
Furthermore, for any θ ∈ F q m and for any z * i = 1/h i + θ and x * j = 1/h j + 1/h 0 + θ , the Cauchy matrix C(z * , x * ) is equal to ∆ ℓ (h).

Indeed, the public generator matrix G is a k × n block matrix where each block is an ℓ × ℓ dyadic matrix with ℓ being a power of 2. The entries of G belong to F q and the integers k and n are chosen such that n = n 0 ℓ and k = nmℓ = ℓ(n 0m) where n 0 is some integer and m defines the extension F q m . The matrix G is obtained from a secret ℓ × n block parity-check matrix

H = ∆ ℓ ( f 0 )| • • • |∆ ℓ ( f n 0 -1 )
where each block ∆ ℓ ( f j ) is an ℓ × ℓ dyadic matrix and for any 0 ≤ j ≤ n 0 -1, f j is a vector of F ℓ q m such that f j = γ j h ω j ℓ⊕d j , h (ω j ℓ+1)⊕d j , . . . , h ((ωj+1)ℓ-1)⊕dj where h = (h 0 , . . . , h N-1 ) is a random vector of F N q m that satisfies Equation ( 14) and such that N = N 0 ℓ for some integer N 0 ≫ n 0 . The integers ω j , d j are chosen such that 0 ≤ ω j ≤ N 0 -1 and 0 ≤ d j ≤ ℓ -1. The coefficients γ j are non zero elements of F q m . Note that the integers ω j 's are different. The secret key consists then of the vectors h, ω = (ω 0 , . . . , ω n 0 -1 ), d = (d 0 , . . . , d n 0 -1 ) and γ = (γ 0 , . . . , γ n 0 -1 ).

C Proof of Proposition 4

Lemma 3. Let N = N 0 ℓ with ℓ = 2 e for some integer e and let h be a vector in F N q m that satisfies Equation [START_REF] Faugère | A new efficient algorithm for computing gröbner bases without reduction to zero : F5[END_REF]. Let G (a * , G) be the Goppa code such that a * = (a * 0 , . . . , a * N-1 ) is defined by a * j = 1/h j + 1/h 0 for 0 ≤ j ≤ N -1 and G(X) = ∏ ℓ-1 i=0 (Xz i ) with z i = 1/h i . Then for any i, j in {0, . . . , N -1} we have a * i⊕ j = a * i + a * j and for any 0 ≤ j ≤ N 0 -1 and 0 ≤ i, i ′ ≤ ℓ -1

G(a * jℓ+i ) -1 = ( j+1)ℓ-1 ∏ l= jℓ h l
Proof. The property that a * i⊕ j = a * i + a * j comes from Equation [START_REF] Faugère | A new efficient algorithm for computing gröbner bases without reduction to zero : F5[END_REF]. Furthermore, we have:

G(a * jℓ+i ) -1 = ℓ-1 ∏ ℓ=0 (z ℓ -a * jℓ+i ) -1 = ℓ-1 ∏ ℓ=0 (1/h ℓ + 1/h jℓ+i + 1/h 0 ) -1 = ℓ-1 ∏ ℓ=0 h jℓ+ℓ
which terminates the proof.

We remark in particular that we have G(a * jℓ+i ) = G(a * jℓ ) for any 0 ≤ j ≤ n 0 -1 and 0 ≤ i ≤ ℓ -1. The next lemma we give without proof shows that the action of a dyadic permutation can be simply characterized as a translation. We are now prepared to prove Proposition 4. Let (h, ω, d, γ) be the private key and let G be the public generator matrix. We shall see that a parity-check matrix for the code generated by G is V ℓ (a, λ ) with a jℓ+i = a * (ω j t+i)⊕d ℓ and λ jℓ+i = γ j G(a * ω j t ) -1 where a * and G(X) are defined as in Lemma 3. Indeed, we know that the code defined by the parity-check matrix ∆ ℓ (h) is also defined by the parity-check matrix V ℓ (a, λ ) where λ j = G(a j ) -1 for any 0 ≤ j ≤ N -1. Recall from Lemma 3 that G(a jℓ+i ) = G(a jℓ ) for any 0 ≤ j ≤ N 0 -1 and 0 ≤ i ≤ ℓ -1. The role of ω is to pick n 0 dyadic blocks from ∆ ℓ (h). These blocks correspond to the columns a * ω j ℓ , . . . , a * (ω j +1)ℓ-1 of V ℓ (a, λ ) when j describes {1, . . . , n 0 }. These columns are then multiplied by a dyadic permutation matrix ∆ ℓ (b d ℓ ) which leads to reorder the columns as a ω j ℓ⊕d j , . . . , a ((ωj+1)ℓ-1)⊕dj according to Lemma 4. Finally, each dyadic block is scaled by γ j which means that if we set λ jℓ+i = γ j G(a ω j ℓ ) -1 then V ℓ (a, λ ) is another parity-check matrix of the code generated by G. We are now going to show that for any 0 ≤ j ≤ n 0 -1 and 0 ≤ i, i ′ ≤ ℓ -1, we have the following equations:

   λ jℓ+i
= λ jℓ a jℓ+i + a jℓ = a i + a 0 a jℓ+i⊕i ′ = a jℓ+i + a jℓ+i ′ + a jℓ [START_REF] Faugère | Gröbner bases of bihomogeneous ideals generated by polynomials of bidegree (1,1): Algorithms and complexity[END_REF] It is clear from Lemma 3 that λ jℓ+i = λ jℓ . On the other hand, a jℓ+i = a (ω j ℓ+i)⊕d j = 1/h (ω j ℓ+i)⊕d j + 1/h 0 . From Equation ( 14) we thus have:

a jℓ+i = 1 h ω j ℓ+i + 1 h d j = 1 h ω j ℓ + 1 h i + 1 h 0 + 1 h d j = 1 h ω j ℓ⊕d j + 1 h i + 1 h 0 = a jℓ + 1 h i + 1 h 0 .
We observe in particular that a i + a 0 = 1/h i + 1/h 0 and since this quantity does not depend on ℓ, this is equivalent to say that a jℓ+i + a jℓ = a i + a 0 . Before proving the third equation, we can first see that a jℓ+i⊕i ′ + a jℓ = a i⊕i ′ + a 0 . So if we know that a i⊕i ′ = a i + a i ′ + a 0 then we would get a i⊕i ′ = a jℓ+i + a jℓ + a i ′ which finally implies a i⊕i ′ = a jℓ+i +a jℓ+i ′ +a 0 that leads to the expected result. Now we have a i⊕i ′ = a (ω 1 ℓ+i+i ′ )⊕d 1 = a ω 1 ℓ+i+i ′ + a d 1 = a ω 1 ℓ+i + a i ′ + a d 1 = a (ω 1 ℓ+i)⊕d 1 + a i ′ . Therefore we obtain:

a i⊕i ′ = a i + a i ′ + a ω 1 ℓ + a d 1 + a ω 1 ℓ + a d 1 = a i + a i ′ + a 0 .

Fig. 1 .

 1 Fig. 1. Algorithm F 4 (modified version)

Lemma 4 .

 4 Let t and d two integers such that0 ≤ d ≤ ℓ -1. For any vector v = (v 0 , . . . , v ℓ-1 ), we have: v × ∆ ℓ (b d ) = (v d , v 1⊕d , . . . , v (ℓ-1)⊕d )(15) where the vector b d = (b d,0 , . . . , b d,ℓ-1 ) is such that b d, j = 0 if j = d and b d,d = 1.

Table 1 .

 1 Cryptanalysis results for[START_REF] Berger | Reducing key length of the McEliece cryptosystem[END_REF] (m = 2)

	Challenge	q	ℓ n 0 d Security [5] Unknowns Equations Time (Operations, Memory)
	A 16 2 8 51 9 3 B 16 2 8 51 10 3 C 16 2 8 51 12 3 D 16 2 8 51 15 4	80 90 100 120	16 18 22 28	510 0.06 sec (2 18.9 op, 115 Meg) 612 0.03 sec (2 17.1 op, 116 Meg) 816 0.05 sec (2 16.2 op, 116 Meg) 1275 0.02 sec (2 14.7 op, 113 Meg)
	A 20 2 10 75 6 2 B 20 2 10 93 6 2 C 20 2 10 93 8 2	80 90 110	10 10 14	337 0.05 sec (2 15.8 op, 115 Meg) 418 0.05 sec (2 17.1 op, 115 Meg) 697 0.02 sec (2 14.5 op, 115 Meg)
	QC 600 2 8 255 15 3	600	28	6820 0.08 sec (2 16.6 op, 116 Meg)

Table 2 .

 2 Cryptanalysis results for[START_REF] Misoczki | Compact McEliece keys from Goppa codes[END_REF].

	Challenge q m	ℓ n 0 Security Unknowns Equations	Time (Operations, Memory)
	Table 2 2 2 8 64 56	128	115 193, 584 1, 776.3 sec (2 34.2 op, 360 Meg)
	Table 2 2 4 4 64 32	128	67 112, 924	0.50 sec (2 22.1 op, 118 Meg)
	Table 2 2 8 2 64 12	128	27	40, 330	0.03 sec (2 16.7 op, 35 Meg)
	Table 3 2 8 2 64 10	102	23	32, 264	0.03 sec (2 15.9 op, 113 Meg)
	Table 3 2 8 2 128 6	136	16	65, 028	0.02 sec (2 15.4 op, 113 Meg)
	Table 3 2 8 2 256 4	168	13 130, 562	0.11 sec (2 19.2 op, 113 Meg)
	Table 5 2 8 2 128 4	80	12	32, 514	0.06 sec (2 17.7 op, 35 Meg)
	Table 5 2 8 2 128 5	112	14	48, 771	0.02 sec (2 14.5 op, 35 Meg)
	Table 5 2 8 2 128 6	128	16	65, 028	0.01 sec (2 16.6 op, 35 Meg)
	Table 5 2 8 2 256 5	192	15 195, 843	0.05 sec (2 17.5 op, 35 Meg)
	Table 5 2 8 2 256 6	256	17 261, 124	0.06 sec (2 17.8 op, 35 Meg)
	Dyadic 256 2 4 4 128 32 Dyadic 512 2 8 2 512 6	256 512	68 455, 196 18 1, 046, 532	7.1 sec (2 26.1 op, 131 Meg) 0.15 sec (2 19.7 op, 38 Meg)

An independent and parallel work[START_REF] Gauthier Umana | Practical key recovery attacks on two McEliece variants[END_REF] took place that also proposed a cryptanalysis of these two schemes.


Laurent Boudin
email: laurent.boudin@upmc.fr

Julien Mathiaud
email: julien.mathiaud@cea.fr

A NUMERICAL SCHEME FOR THE ONE-DIMENSIONAL PRESSURELESS GASES SYSTEM

 for the continuous model.

Introduction

During the last two decades, there have been many contributions on the pressureless gases system, and it seems natural to tackle the question of its discretization. The pressureless gases system appears as a system of conservation laws on the mass and momentum. Hence, it is relevant to wonder if standard numerical schemes for conservation laws, like the upwind scheme, for instance, are fitted to this particular system. However, we emphasize that it is a degenerate hyperbolic system (the Jacobian is not diagonalizable).

Let us now recall the one-dimensional system describing a pressureless gas. Let T > 0. The gas density ̺(t, x) ≥ 0 and the momentum q(t, x) ∈ R satisfy the following equations in (0, T ) × R

∂ t ̺ + ∂ x (̺u) = 0, (1) ∂ t q + ∂ x (qu) = 0. (2)
One must define the velocity u(t, x) ∈ R as a quotient of q by ̺, but this may not be possible, since ̺ can be zero. We discuss this issue below, by recalling the notion of duality solutions [START_REF] Bouchut | One-dimensional transport equations with discontinuous coefficients[END_REF]. As already stated, each equation consists of a conservation law, (1) for mass and (2) for momentum. We obviously need initial conditions [START_REF] Berthon | A relaxation scheme for the approximation of the pressureless Euler equations[END_REF] ̺(0, •) = ̺ in , q(0, •) = q in , Date: May 5, 2011. in which the condition on the momentum can be replaced by an initial condition on the velocity u(0, •) = u in , and then written again as q(0, •) = ̺ in u in .

The previous system can be seen as a simplified model of the Euler equations, where the pressure has been set to zero. It can describe either cold plasmas or galaxies' dynamics [START_REF] Ya | Gravitational instability: An approximate theory for large density perturbations[END_REF]. This system (1)-( 2) and related problems (traffic models, magnetohydrodynamics, astrophysics, pressureless fluid equations...) have been widely studied, see, for instance, [START_REF] Bouchut | On zero pressure gas dynamics[END_REF][START_REF] Grenier | Existence globale pour le système des gaz sans pression[END_REF][START_REF] Rykov | Generalized variational principles, global weak solutions and behavior with random initial data for systems of conservation laws arising in adhesion particle dynamics[END_REF][START_REF] Brenier | Sticky particles and scalar conservation laws[END_REF][START_REF] Poupaud | Measure solutions to the linear multi-dimensional transport equation with non-smooth coefficients[END_REF][START_REF] Bouchut | Duality solutions for pressureless gases, monotone scalar conservation laws, and uniqueness[END_REF][START_REF] Boudin | A solution with bounded expansion rate to the model of viscous pressureless gases[END_REF][START_REF] Sever | An existence theorem in the large for zero-pressure gas dynamics[END_REF][START_REF] Berthelin | Existence and weak stability for a pressureless model with unilateral constraint[END_REF][START_REF] Poupaud | Diagonal defect measures, adhesion dynamics and Euler equation[END_REF][START_REF] Leveque | The dynamics of pressureless dust clouds and delta waves[END_REF][START_REF] Bouchut | Uniqueness and weak stability for multidimensional transport equations with one-sided Lipschitz coefficient[END_REF][START_REF] Gallagher | On pressureless gases driven by a strong inhomogeneous magnetic field[END_REF][START_REF] Coulombel | From gas dynamics to pressureless gas dynamics[END_REF][START_REF] Berthelin | A model for the formation and evolution of traffic jams[END_REF][START_REF] Nguyen | Pressureless Euler/Euler-Poisson systems via adhesion dynamics and scalar conservation laws[END_REF][START_REF] Brenier | A modified least action principle allowing mass concentrations for the early universe reconstruction problem[END_REF]. Those references use the same fluid point of view we choose here, or the kinetic one, involving the adhesion dynamics of the so-called sticky particles.

When one studies smooth solutions of the pressureless gases system, (2) can be replaced by the standard Burgers equation:

(4) ∂ t u + ∂ x u 2 2 = ∂ t u + u∂ x u = 0.
System (1)-( 2) is then uncoupled, since we obtain u from (4), and then ̺ from [START_REF] Berthelin | Existence and weak stability for a pressureless model with unilateral constraint[END_REF]. On the other hand, it is well-known that smooth initial data can result in mass concentration, for example, when the velocity does not increase. In that case, the velocity cannot satisfy (4) anymore.

In [START_REF] Bouchut | One-dimensional transport equations with discontinuous coefficients[END_REF], Bouchut and James introduced the notion of duality solution for one-dimensional transport equations and conservation laws. In [START_REF] Bouchut | Duality solutions for pressureless gases, monotone scalar conservation laws, and uniqueness[END_REF], they prove that this framework is fitted to the pressureless gases system. Let us briefly recall the results they obtain. Definition 1. A couple (̺, q) ∈ C(R + ; w*-M loc (R)) 2 , with ̺ ≥ 0, is a duality solution to (1)-(2), if there exists a bounded Borel function a and α ∈ L 1 loc (R * + ) such that ∂ x a ≤ α, q = a̺, in R * + × R, and, in the duality sense on (t 1 , t 2 ) × R, for any 0 < t 1 < t 2 ,

∂ t ̺ + ∂ x (̺a) = 0, ∂ t q + ∂ x (qa) = 0.
In that setting, u is defined ̺-almost everywhere, and we have u = a ̺-a.e. Bouchut and James prove that duality solutions are stable, and also entropic, i.e. the following inequality holds, in the distributional sense, (5)

∂ t (̺S(u)) + ∂ x (̺uS(u)) ≤ 0,
for any convex function S. Using those properties and the sticky particles dynamics, they obtain the following existence result.

Theorem 1. Let ̺ in , q in ∈ M loc (R), with ̺ in ≥ 0 and |q in | ≤ U ̺ in , U ≥ 0.
Then there exists a duality solution to (1)-( 3), and we have a ∞ ≤ U and

α(t) = 1/t.
As proven in [START_REF] Hoff | The sharp form of Oleȋnik's entropy condition in several space variables[END_REF], the one-sided Lipschitz (OSL) condition on the expansion rate ∂ x a ≤ 1/t, also known as the Oleinik entropy condition, is optimal for a convex scalar conservation law. In the proof of Theorem 1, it is clear that the standard convex entropy condition (5) is not enough, and the OSL condition is really required. Note that, when the solutions are smooth, this estimate can easily be proven, since the Burgers equation ( 4) lies in the class of convex scalar conservation laws [START_REF] Brenier | The discrete one-sided Lipschitz condition for convex scalar conservation laws[END_REF][START_REF] Leveque | Numerical methods for conservation laws[END_REF].

Eventually, Bouchut and James also obtain uniqueness when ̺ in is nonatomic (essentially meaning that ̺ in is smooth).

In this work, we also consider the viscous pressureless gases system. In this system, as explained in [START_REF] Boudin | A solution with bounded expansion rate to the model of viscous pressureless gases[END_REF], ( 2) is replaced by an equation on the velocity itself. Let us choose ε > 0. The gas density ̺(t, x) > 0 and the velocity u(t, x) ∈ R satisfy, in (0, T ) × R, Equation ( 1) and ( 6)

∂ t u + u∂ x u = ε ̺ ∂ 2 xx u,
with the same set of initial conditions (3). That writing imposes that ̺ remains nonnegative, which is true if ̺ in is also nonnegative, see [START_REF] Boudin | A solution with bounded expansion rate to the model of viscous pressureless gases[END_REF]. Note that ( 6) is equivalent, when ε is fixed, to

∂ t u + ∂ x u 2 2 = ε ̺ ∂ 2 xx u, (7)
if we take into account the smoothness of the viscous velocity given in [START_REF] Boudin | A solution with bounded expansion rate to the model of viscous pressureless gases[END_REF]. In fact, ( 6) or [START_REF] Bouchut | Uniqueness and weak stability for multidimensional transport equations with one-sided Lipschitz coefficient[END_REF] can also be rewritten as an equation on the momentum, with a viscosity term ε∂ 2 xx u on the right-hand side, ∂ t (̺u) + ∂ x (̺u 2 ) = ε∂ 2 xx u, which yields (2) when ε goes to 0. In [START_REF] Boudin | A solution with bounded expansion rate to the model of viscous pressureless gases[END_REF], the author proved the existence, in the sense of distributions, of solutions to the viscous system (1), ( 3) and [START_REF] Bouchut | Uniqueness and weak stability for multidimensional transport equations with one-sided Lipschitz coefficient[END_REF], and that the expansion rate is upper-bounded: ∂ x u ≤ A/(At + 1), when A = max(ess sup ∂ x u in , 0) is finite. He also obtains the convergence of the viscous solutions towards the duality solutions to the pressureless gases system when ε vanishes. More precisely, the following convergence result holds.

Theorem 2. Let (̺ in ε ), (u in ε ) such that, for any ε > 0,

̺ in ε > 0, ̺ in ε ∈ L ∞ (R), 1/̺ in ε L ∞ (R) ≤ Cε -1/4 , u in ε ∈ L 1 ∩ L ∞ (R), u in ε L ∞ (R) ≤ C, ∂ x u in ε L 1 ∩ L 2 (R), ess sup ∂ x u in ε ≤ Cε -1/2 . We assume that (̺ in ε ) ⇀ ̺ in and (̺ in ε u in ε ) ⇀ q in in w*-M loc (R).
Then, up to a subsequence, (̺ ε , ̺ ε u ε ), given by the solutions to (1) and (7), with initial datum (̺ in ε , ̺ in ε u in ε ), converges in C t (w*-M loc (R)) towards the duality solution (̺, q) of (1)-(3).

Both viscous and inviscid systems can also be studied in a periodic framework, i.e. we focus on the closed interval [0, 1] and impose that all the physical quantities have the same values at both x = 0 and x = 1, so that the solutions are 1-periodic. This work is dedicated to the numerical approximation of the pressureless gases system (1)- [START_REF] Berthelin | A model for the formation and evolution of traffic jams[END_REF], where the latter may be replaced by [START_REF] Bouchut | Uniqueness and weak stability for multidimensional transport equations with one-sided Lipschitz coefficient[END_REF]. For readability reasons, we choose the periodic framework.

There are two methods to get a priori relevant schemes. The first one is to use the natural kinetic framework which underlies the pressureless gas dynamics, with kinetic schemes, as in [START_REF] Bouchut | On zero pressure gas dynamics[END_REF][START_REF] Bouchut | Numerical approximations of pressureless and isothermal gas dynamics[END_REF], or with particle methods [START_REF] Chertock | A new sticky particle method for pressureless gas dynamics[END_REF]. The second one is related to the discretization of hyperbolic conservation laws. Gosse and James [START_REF] Gosse | Numerical approximations of one-dimensional linear conservation equations with discontinuous coefficients[END_REF] point out the relevance of two families of numerical schemes: the upwind schemes and the Lax-Friedrichs schemes. In [START_REF] Berthon | A relaxation scheme for the approximation of the pressureless Euler equations[END_REF], Berthon et al. investigate a relaxation scheme for the pressureless gases system in one and two-dimensional settings.

As we already pointed out, the key condition to obtain the duality solution is that the velocity expansion rate must be upper-bounded by 1/t. Brenier and Osher [START_REF] Brenier | The discrete one-sided Lipschitz condition for convex scalar conservation laws[END_REF] obtained the relevance of the OSL condition in a discrete framework for the convex scalar conservation laws. In this work, we first investigate the upwind scheme associated to (1)- [START_REF] Berthelin | A model for the formation and evolution of traffic jams[END_REF], and prove that it fails to ensure the OSL condition. Subsequently, we try the upwind diffusive scheme associated to (1) and [START_REF] Bouchut | Uniqueness and weak stability for multidimensional transport equations with one-sided Lipschitz coefficient[END_REF], and explain how we can obtain a good numerical approximation of the duality solution to the inviscid pressureless gases system using this scheme. We do not study the Lax-Friedrichs schemes described in [START_REF] Gosse | Numerical approximations of one-dimensional linear conservation equations with discontinuous coefficients[END_REF]. Indeed, the numerical dissipation induced by those first order schemes is too significant. Since it is then natural to use higher order schemes, we recover the same kind of problems as in the diffusive upwind scheme we here propose, involving second order terms.

In the remainder, let ∆t, ∆x > 0 such that N = T /∆t ∈ N * and I = 1/∆x ∈ N * , and set λ = ∆t/∆x. We respectively denote ̺ n i , q n i and u n i the approximate values of ̺, q and u at time n∆t ∈ [0, T ] and coordinate (i + 1/2)∆x ∈ [0, 1), for 0 ≤ n ≤ N and 0 ≤ i < I. Since we use a periodic framework, we define ̺ n i , q n i and u n i for any i ∈ Z, by

̺ n i+pI = ̺ n i , q n i+pI = q n i , u n i+pI = u n i , 0 ≤ i < I, p ∈ Z * .
For the sake of readability, in the previous notations, we may drop the time iteration index n and replace n + 1 by a prime symbol " ′ ". For instance, the velocity at time (n + 1)∆t and coordinate (i + 1/2)∆x can be written as u ′ i or u n+1 i . Apart from the density, momentum and velocity, the quantity of interest, which we name the numerical expansion rate, is, for each time and space indices n and i,

w n i := nλ(u n i+1 -u n i ).
Indeed, the OSL condition at time n∆t then reads max i w n i ≤ 1.

Upwind scheme

Let us first denote the positive and negative parts of a ∈ R a + = max(0, a), a -= min(0, a).

The upwind scheme writes, for any 0 ≤ i < I,

̺ ′ i = ̺ i -λ ̺ i (u i ) + -̺ i-1 (u i-1 ) + -λ ̺ i+1 (u i+1 ) --̺ i (u i ) -, q ′ i = q i -λ q i (u i ) + -q i-1 (u i-1 ) + -λ q i+1 (u i+1 ) --q i (u i ) -, u ′ i = q ′ i ̺ ′ i , if ̺ ′ i > 0,
and

u ′ i is not defined if ̺ ′ i = 0.
It is quite clear that the previous schemes on both ̺ and q are monotonic, if the standard Courant-Friedrichs-Lewy (CFL) condition λ max |u| ≤ 1 is satisfied. Hence, we only choose positive initial data to study positive velocities. The scheme then becomes, for any 0 ≤ i < I,

̺ ′ i = (1 -λu i )̺ i + λu i-1 ̺ i-1 , (8) q
′ i = (1 -λu i )q i + λu i-1 q i-1 , (9)
u ′ i = q ′ i ̺ ′ i . (10)
As we already stated, this last equality allows to define u ′ i only when ̺ ′ i > 0. This fits the mathematical setting of the pressureless gases system, since u can only be defined ̺-almost everywhere. Nevertheless, it is not satisfying from a numerical viewpoint, since the computations stop whenever the density becomes equal to 0. We can impose whichever value we want, for instance, u ′ i = 0, when ̺ ′ i = 0. Indeed, we do not care about the value of the velocity at a point where is there is no matter. But we must keep in mind not to use those artificial nil values of u ′ i to study the numerical expansion rate.

Thanks to [START_REF] Boudin | Convergence of a viscous numerical scheme towards the inviscid one-dimensional pressureless gases system[END_REF], we immediately have

̺ ′ i ̺ ′ i+1 (u ′ i+1 -u ′ i ) = ̺ ′ i q ′ i+1 -̺ ′ i+1 q ′ i , which implies ̺ ′ i ̺ ′ i+1 w ′ i (n + 1)λ = (1 -λu i+1 )̺ i+1 ̺ ′ i w i nλ + λu i-1 ̺ i-1 ̺ ′ i+1 w i-1 nλ .
Under the CFL condition λ max |u| ≤ 1, if (w i ) 0≤i<I are negative, and if (̺ ′ i ) 0≤i<I are nonnegative, it is clear that the quantities (w ′ i ) also remain negative. Unfortunately, if w j is nonnegative for a given j, the OSL condition w ′ i ≤ 1 for all i may not be satisfied.

Proposition 3. Let 0 < λ < 1 and U such that 0 < λU < 1, and choose an integer I > 2 + 1/λ. We consider the following initial data

(11) ̺ 0 i = 1, 0 ≤ i ≤ I -1, u 0 0 = U, u 0 i = 0, 1 ≤ i ≤ I -1.
Then the upwind scheme (8)- [START_REF] Boudin | Convergence of a viscous numerical scheme towards the inviscid one-dimensional pressureless gases system[END_REF] does not satisfy the OSL condition. More precisely, we have [START_REF] Brenier | Sticky particles and scalar conservation laws[END_REF] max

i w I-2 i > U.
The assumption on the Courant number 0 < λ < 1 is standard and is a natural consequence of the CFL condition λU ≤ 1 when U is large. The initial data can of course be defined without the discretization grid: we have ̺ in ≡ 1 and u in ≡ 0 except in 0 where u in (0) = U .

Proof. It is easy to simultaneously prove, by induction on the time index 0 ≤ j < I -1, that

̺ j i > 0, 0 ≤ i ≤ I -1, and u j 0 = U, u j j+1 > 0, u j i = 0, j + 2 ≤ i ≤ I -1.
We must emphasize that the nil values of u j i are obtained because q j i = 0 and ̺ j i > 0, and not because of the choice of nil velocity when the density equals 0.

Then we can write that Note that the first inequality is in fact an equality, but we do not need to prove it here.

w I-2 I-1 = (I -2)λ(u I-2 0 -u I-2 I-1 ) = (I -2)
Remark 1. As we already pointed out, the standard numerical version of the OSL condition reads max i w n i ≤ 1. It may have been relaxed into max i w n i ≤ K, where K is a nonnegative constant, which does not depend on the initial data. But [START_REF] Brenier | Sticky particles and scalar conservation laws[END_REF] implies that the quantity max i w I-2 i can be as large as we want, depending on the value of U . Proposition 3 means in particular that, if the space step ∆x is refined enough, the numerical OSL condition cannot be satisfied anymore, with initial data given by [START_REF] Brenier | A modified least action principle allowing mass concentrations for the early universe reconstruction problem[END_REF]. Moreover, we must point out that, whatever the final time is, one can find a discretization for which the upwind scheme cannot satisfy the OSL condition, because I does not depend on T .

The initial datum u in in the previous proposition is not smooth. Nevertheless, even with smooth (and periodic) initial data, the upwind scheme does not necessarily provide a solution satisfying the OSL condition, see 4.1.

Adding an artificial viscosity

As it was done in [START_REF] Boudin | A solution with bounded expansion rate to the model of viscous pressureless gases[END_REF], we now add a small viscosity term in (2) to obtain [START_REF] Bouchut | Uniqueness and weak stability for multidimensional transport equations with one-sided Lipschitz coefficient[END_REF], and we study the numerical approximation of ( 1) and [START_REF] Bouchut | Uniqueness and weak stability for multidimensional transport equations with one-sided Lipschitz coefficient[END_REF]. We still deal with arbitrary 1-perodic initial data

u in ≥ 0, ̺ in ≥ 0.
In what follows, we consider a fixed ε > 0, small enough. If necessary, we regularize both u in and ̺ in so that (keeping the same notations for both, even if they depend on ε) u in , ̺ in ∈ C 1 (R; R * + ), remain periodic and satisfy the assumptions of Theorem 2, which can be written as

(13) ̺ in (x) ≥ Cε 1/4 , u in (x) ≤ C, (u in ) ′ (x) ≤ C √ ε , ∀x ∈ [0, 1],
where C is a constant which does not depend on ε. The regularized ̺ in must lie in R * + , since the continuous diffusive model involves a division by ̺. In the following, we set

U = max [0,1] u in > 0, V = min [0,1] u in > 0, A = max(0, max [0,1] (u in ) ′ ) ≥ 0, R = min [0,1]
̺ in > 0.

The previous quantities can depend on ε, and must satisfy properties which come from (13), i.e.

(14) R ≥ Cε 1/4 , V ≤ U ≤ C, A ≤ C √ ε ,
where C does not depend on ε.

Then we consider ∆t, ∆x > 0, and set

λ = ∆t ∆x , σ = ∆t ∆x 2 .
In the remainder, we make the following assumptions on the time and space steps:

0 < ∆x ≤ 2V 1 + A , ( 15
) 0 < ∆t ≤ min 1 4A + 1 , 1 4U ∆x, R 4ε(1 + AT ) ∆x 2 . ( 16
)
In fact, [START_REF] Coulombel | From gas dynamics to pressureless gas dynamics[END_REF] and ( 16) are not so restrictive, since, eventually, ∆x and ∆t will go to 0, ε being fixed. From now on, even if we do not write the dependence on ε, we must keep in mind, in the numerical examples, that U , V , A and R can depend on ε and must satisfy [START_REF] Chertock | A new sticky particle method for pressureless gas dynamics[END_REF], at least for ε small enough. That dependence implies that, at most, ∆x is of order √ ε and ∆t of order 4 √ ε.

Note that it cannot prevent ∆t and ∆x from going to 0 while ε remains fixed.

With the same notations for quantities at times n∆t and (n + 1)∆t as in Section 2, we now focus on the following scheme, corresponding to the discretization of ( 1) and [START_REF] Bouchut | Uniqueness and weak stability for multidimensional transport equations with one-sided Lipschitz coefficient[END_REF].

u ′ i = u i -λ u i 2 2 - u i-1 2 2 + εσ ̺ i (u i-1 + u i+1 -2u i ), (17)
̺ ′ i = (1 -λu ′ i )̺ i + λu ′ i-1 ̺ i-1 . (18)
Note that ( 17) is obtained from [START_REF] Bouchut | Uniqueness and weak stability for multidimensional transport equations with one-sided Lipschitz coefficient[END_REF], which is written under a conservative form, as suggested in [START_REF] Brenier | The discrete one-sided Lipschitz condition for convex scalar conservation laws[END_REF].

If we choose u in ≡ 1, we can note that both upwind and diffusive schemes give u n i = 1 for any i and n, which is reassuring: in that case, and when ̺ remains nonnegative, the velocity satisfies the Burgers equation, which implies, at least formally, that u remains constant.

Remark 2. The velocity terms which appear in [START_REF] Gosse | Numerical approximations of one-dimensional linear conservation equations with discontinuous coefficients[END_REF] are the ones at time (n + 1)∆t. They must not be at time n∆t to ensure the lower bound on ̺, as we shall see in the proof of Theorem 4 below.

Numerical strategy. Let us here sum up the strategy used to build a relevant numerical solution to the pressureless gases system.

1. Consider 1-periodic initial data. 2. Fix ε > 0 small enough.

3. Regularize ̺ in , u in so that they become C 1 (R; R * + ) and satisfy (13). 4. Fix ∆x and ∆t satisfying ( 15)-( 16). 5. Use the numerical scheme ( 17)- [START_REF] Gosse | Numerical approximations of one-dimensional linear conservation equations with discontinuous coefficients[END_REF]. The previous strategy holds for two reasons. First, the following theorem states that the scheme ( 17)-( 18) is L ∞ -stable, consistent, monotonic, and that it satisfies the OSL condition. Consequently, (̺ n i ) and (u n i ) converge towards to ̺ and u, solutions to the viscous pressureless gases system when both ∆t and ∆x go to 0, ε being fixed. Second, thanks to Theorem 2, the scheme eventually provides a good approximation of a solution to the inviscid pressureless gases system, if one chooses ε small enough, and regularized initial data in C 1 (R; R * + ) close to the original ones and satisfying [START_REF] Brenier | The discrete one-sided Lipschitz condition for convex scalar conservation laws[END_REF]. The error between the diffusive numerical and the duality solutions is currently under study, see [START_REF] Boudin | Convergence of a viscous numerical scheme towards the inviscid one-dimensional pressureless gases system[END_REF].

Theorem 4. We assume that (15)-( 16) hold. Then we have, for any i and n ≥ 0,

V ≤ u n i ≤ U, (19)
u n i -u n i-1 ≤ A∆x 1 + An∆t , (20)
̺ n i ≥ R 1 + An∆t ≥ R 1 + AT > 0. (21)
Moreover, the discrete total mass is conserved, i.e., for any n ≥ 0, [START_REF] Leveque | The dynamics of pressureless dust clouds and delta waves[END_REF] i

̺ n i ∆x = i ̺ 0 i ∆x.
Finally, when ε > 0 is fixed, the scheme (17)-( 18) is consistent with (1) and [START_REF] Bouchut | Uniqueness and weak stability for multidimensional transport equations with one-sided Lipschitz coefficient[END_REF], is first order accurate in time and space, and is monotonic.

Equations ( 19) and ( 21) respectively correspond to the maximum principles on the velocity and the density, [START_REF] Hoff | The sharp form of Oleȋnik's entropy condition in several space variables[END_REF] stands for the discrete version of the OSL condition.

Remark 3. The assumptions [START_REF] Rykov | Generalized variational principles, global weak solutions and behavior with random initial data for systems of conservation laws arising in adhesion particle dynamics[END_REF] on ∆t ensure the stability of the scheme. More precisely, the second one is induced by the CFL condition and the third one is similar to standard stability conditions for explicit diffusive schemes. The first one is needed for the required properties of the scheme, as it will be detailed in the proof of Theorem 4.

Proof. We proceed by induction on n ∈ N, and first investigate the case when n = 0. Equations ( 19) and ( 21) are obviously satisfied by definitions of U , V and R, and thanks to [START_REF] Coulombel | From gas dynamics to pressureless gas dynamics[END_REF]. The fact that (20) holds comes from the fact that u in is smooth, and consequenly satisfies the intermediate values inequality.

In the remaining of the proof, we suppose that A > 0. The case when A = 0 can easily be treated. Let us assume that ( 19)-( 21) hold for a fixed n, and prove them for n + 1. We can rewrite Equation ( 17) as [START_REF] Nguyen | Pressureless Euler/Euler-Poisson systems via adhesion dynamics and scalar conservation laws[END_REF]

u ′ i = 1 -λ u i + u i-1 2 - 2εσ ̺ i u i + εσ ̺ i u i+1 + λ u i + u i-1 2 + εσ ̺ i u i-1 .
Under this form, u ′ i is a convex combination of u i-1 , u i and u i+1 , since the corresponding coefficients in [START_REF] Nguyen | Pressureless Euler/Euler-Poisson systems via adhesion dynamics and scalar conservation laws[END_REF] live in [0, 1] and their sum equals 1. Indeed, we clearly have, thanks to ( 16) and ( 21),

0 ≤ 2εσ ̺ i ≤ 1 2 ,
and, thanks to ( 16) and ( 19),

0 ≤ λ u i + u i-1 2 ≤ 1 4 .
Then it is easy to check that u ′ i satisfies [START_REF] Grenier | Existence globale pour le système des gaz sans pression[END_REF]. Let us now define, for any i,

δ i = u i+1 -u i - A∆x 1 + An∆t ,
which we know is negative, and prove that δ ′ i is also negative, for any i. Thanks to [START_REF] Nguyen | Pressureless Euler/Euler-Poisson systems via adhesion dynamics and scalar conservation laws[END_REF], we can write

u ′ i+1 -u ′ i = 1 - λ 2 (u i+1 + u i ) - εσ ̺ i - εσ ̺ i+1 (u i+1 -u i ) + εσ ̺ i+1 (u i+2 -u i+1 ) + λ 2 (u i + u i-1 ) + εσ ̺ i (u i -u i-1
).

Then we have

δ ′ i = εσ ̺ i+1 δ i+1 + λ 2 (u i + u i-1 ) + εσ ̺ i - Aλ∆x 2(1 + An∆t) δ i-1 + 1 - λ 2 (u i+1 + u i ) - εσ ̺ i - εσ ̺ i+1 - Aλ∆x 2(1 + An∆t) δ i + A∆x 1 + An∆t 1 - A∆t 1 + An∆t - A∆x 1 + A(n + 1)∆t .
The coefficient before δ i+1 is clearly positive. Let us check that the ones before δ i-1 and δ i are positive too. We have

Aλ∆x 2(1 + An∆t) ≤ λ 2 (u i + u i-1 ) and λ 2 (u i+1 + u i ) + εσ ̺ i + εσ ̺ i+1 + A∆t 2(1 + An∆t) ≤ 1,
because of ( 15)-( 21). Since the (δ i ) are all negative, we still have to prove that the remaining term is negative to get δ ′ i ≤ 0. After simplifying by A∆x, which has no influence on the sign, we write

1 + A(n -1)∆t (1 + An∆t) 2 - 1 1 + A(n + 1)∆t = -(A∆t) 2 (1 + An∆t) 2 (1 + A(n + 1)∆t) ,
which is clearly negative, and ensures that (20) holds for n + 1.

We now focus on the properties of ̺. We successively have, thanks to [START_REF] Leveque | Numerical methods for conservation laws[END_REF] for n and (20) for n + 1,

̺ ′ i ≥ 1 - A∆t 1 + A(n + 1)∆t R 1 + An∆t = R 1 + A(n + 1)∆t ,
which concludes the induction. Note that, as we pointed out in Remark 2, if [START_REF] Gosse | Numerical approximations of one-dimensional linear conservation equations with discontinuous coefficients[END_REF] only involved velocities at time n∆t, the previous inequality would not hold, and we would not get any maximum principle on 1/̺.

We easily notice that in the equality i

̺ ′ i ∆x = i ̺ i ∆x -λ∆x i ̺ i u ′ i + λ∆x i ̺ i-1 u ′ i-1 ,
the last two terms cancel, which ensures the discrete total mass conservation.

Finally, let us investigate some basic properties of the scheme ( 17)-( 18). The consistency is quite clear. Moreover, if we study u ′ i as a function of u i-1 , u i and u i+1 , we immediately have

∂u ′ i ∂u i-1 = λu i-1 + εσ ̺ i ≥ 0, ∂u ′ i ∂u i = 1 - 2εσ ̺ i -λu i ≥ 0, ∂u ′ i ∂u i+1 = εσ ̺ i ≥ 0,
which ensures the required property of monotonicity for [START_REF] Gallagher | On pressureless gases driven by a strong inhomogeneous magnetic field[END_REF], whereas it is clear for [START_REF] Gosse | Numerical approximations of one-dimensional linear conservation equations with discontinuous coefficients[END_REF].

That ends the proof of Theorem 4.

Remark 4. Let us check the behavior of the numerical total momentum. Indeed, in its continuous version [START_REF] Bouchut | Uniqueness and weak stability for multidimensional transport equations with one-sided Lipschitz coefficient[END_REF], the total momentum is conserved, since all the terms besides the time derivative of ̺u are partial derivatives in x.

Unfortunately, the scheme does not ensure the exact conservation of the total momentum. Nevertheless, we can write

i q ′ i = i ̺ i u ′ i + λ i ̺ i u ′ i (u ′ i+1 -u ′ i ),
which implies the following inequalities

[1 -λ(U -V )] i ̺ i u ′ i ≤ i q ′ i ≤ 1 + min 1 n + 1 , λ(U -V ) i ̺ i u ′ i .
Then we have to study the behavior of the quantity

i ̺ i u ′ i = i q i - λ 2 i ̺ i (u i 2 -u i-1 2 ),
for which we have

i q i -U min 1 n , λ(U -V ) i ̺ 0 i ≤ i ̺ i u ′ i ≤ i q i + λV (U -V ) i ̺ 0 i .
We eventually can write

i q ′ i ≥ [1 -λ(U -V )] i q i -U min 1 n , λ(U -V ) i ̺ 0 i , i q ′ i ≤ 1 + min 1 n + 1 , λ(U -V ) i q i + λV (U -V ) i ̺ 0 i ,
which is not really satisfactory. Nevertheless, since the time and space steps satisfy (16), we have

λ ≤ R 4ε(1 + AT ) ∆x,
which ensures that λ is small when both ∆x and ∆t go to 0, and ε > 0 is fixed. Of course, that will not prevent the numerical total momentum from varying, but, at least, from one time step to the next one, the variations have to remain small. It is interesting to note that, in the examples of the next section, the total momentum conservation almost holds, meaning that the previous estimates may be improved in some cases.

Numerical examples

As we already pointed out, a significant drawback of our scheme ( 17)-( 18) is that it does not ensure the exact conservation of the total momentum, since it involves a scheme on the velocity and not on the momentum. Moreover, initial data with vacuum need to be regularized since our scheme cannot stand nil values of ̺. In this section, apart from checking that the OSL condition is satisfied (or not, if studying the behavior of the upwind scheme), we shall also study the numerical total momentum.

Of course, we choose the time and space steps in the following tests such that the CFL condition is satisfied when using the upwind scheme, and ( 15)-( 16) when using the diffusive scheme. 4.1. Nil velocity almost everywhere. This test is the one described in Proposition 3 to prove that the OSL condition was eventually not satisfied by the upwind scheme. We choose ε = 10 -6 . The (regularized) initial data are given by ̺ in ≡ 1 and

u in (x) =              U + ε 2 + U -ε 2 cos πx √ ε if 0 ≤ x ≤ √ ε, ε if √ ε ≤ x ≤ 1 - √ ε, U + ε 2 - U -ε 2 cos π √ ε x -1 + √ ε if 1 - √ ε ≤ x ≤ 1,
We immediately check that min

u in = ε, max u in = U , max(u in ) ′ ≤ U π 2 √ ε
and min ̺ in = 1. We numerically choose U = 1. The space step is set to ∆x = 10 -4 on [0, 1], i.e. I = 10 4 , and the Courant number to λ = 0.25, so that ∆t = 2.5 10 -5 . We perform 100 iterations in time, i.e. T = 2.5 10 -3 s. Eventually, it is clear, on Figure 1, that the diffusive scheme is more efficient

Piecewise linear velocity.

There are other situations when the upwind scheme does not satisfy the OSL condition. For instance, let us consider the following set of initial data [START_REF] Poupaud | Diagonal defect measures, adhesion dynamics and Euler equation[END_REF] ̺ in (x) = 1,

u in (x) = 1 -x ≥ 0, ∀x ∈ [0, 1),
extended by 1-periodicity on R. In both tests, we choose T = 1.2 and ∆x = 10 -4 .

4.2.1. Using the upwind scheme. Using the upwind scheme implies choosing the Courant number λ so that the CFL condition holds. We set λ = 0.1, which ensures λ max u < 1. Then, on Figure 2, the positive part of the numerical expansion rate w is plotted on [0, 1]. It is then clear that there are some values of i such that w i > 1, and, in anticipation of the next paragraph, we must point out that, of course, choosing a lower Courant number does not have any effect on the behavior of the numerical expansion rate. 4.2.2. Using the diffusive scheme. We choose ε = 0.001. As explained in 3, the initial data must be regularized: both ̺ in and u in must be C 1 (R; R * + ), and u in is regularized near 0 in order to have a reasonable periodic agreement with the value in 1, and satisfy [START_REF] Brenier | The discrete one-sided Lipschitz condition for convex scalar conservation laws[END_REF]. Since ( 15)-( 16) must hold, it is possible to check that (λ = 0.01, ∆t = 10 -6 ) is a relevant choice. This time, the OSL condition is satisfied, as one can see on Figure 3 at x = 0.1, where the upwind scheme experiences trouble with the expansion rate for times smaller than 0.2.

Eventually, to investigate the total numerical momentum, on Figure 4, we show its behavior with respect to t, till T , and the result is quite convincing. On the same figure, we also show the total numerical mass, which is of course exactly conserved. Of course, the upwind scheme may often provide a numerical solution satisfying the OSL condition. It is then interesting to check the behavior of both upwind and diffusive schemes, which should be similar. We consider the following initial data for the density

̺ in (x) = 1, 0 ≤ x < 0.2, ̺ in (x) = 0.5, 0.2 ≤ x < 1,
and for the velocity

u in (x) = 0.5(1 -cos(10πx)), 0 ≤ x < 1,
extended by 1-periodicity on R. The final time is T = 2.

For the diffusive scheme, we pick ε = 10 -12 . Since min ̺ in = 0.5, the regularization of the initial density can be chosen not depending on ε. Since u in is already C 1 , we need no regularization, but we have to add a nonnegative term to ensure that min u in > 0, for instance, V = 0.032. And we note that max(u in ) ′ = 5π.

Then we choose ∆x = 0.002, and λ = 0.1 for both upwind and diffusive cases. The space step satisfies (15), as required, since ∆x ≤ 2V 5π . Let us get into some more details of the behavior of both schemes with respect to time. For small times, one can check on Figures 67that both schemes give very similar results for ̺, u and w. If we accurately study Figure 6b, we can see that the upwind scheme has very small variations with respect to the diffusive scheme near some points, which are in fact the jump points of the density, see Figure 8a It is important to note that the numerical expansion rates are still upperbounded by 1, for both schemes. The differences between the numerical solutions is consequently not related to the OSL condition. In fact, we believe that the diffusive scheme is more trustworthy. Indeed, the upwind scheme has natural numerical diffusion, which is responsible for the variations. This numerical diffusion seems to be fully avoided by the diffusive scheme: it is absorbed by the artificial viscosity inserted in the scheme, and its effect cannot numerically appear.

  λU, which, together with the choice of I, implies max i w I-2 i ≥ (I -2)λU > U.

Figure 1 .

 1 Figure 1. Positive part of the numerical expansion rate near 1 at final time T

Figure 2 .

 2 Figure2. "Upwind" plot of w + at t = 0.2 s with initial data[START_REF] Poupaud | Diagonal defect measures, adhesion dynamics and Euler equation[END_REF]

Figure 3 .

 3 Figure 3. "Diffusive" plot of w + at x = 0.1 with regularized initial data[START_REF] Poupaud | Diagonal defect measures, adhesion dynamics and Euler equation[END_REF]

Figure 4 . 3 .

 43 Figure 4. Numerical total mass and momentum

Figure 5 .

 5 Figure 5. (a) Numerical total mass and momentum computed with the diffusive scheme, (b) zoom on the total momentum First, we check on Figure 5a-b that the numerical total momentum is still well conserved by the diffusive scheme.

  .

Figure 6 .

 6 Figure 6. (a) Density at 0.04 s, (b) velocity at 0.2 s

Figure 7 .

 7 Figure 7. Numerical expansion rate at time 0.2 s

Figure 8 .Figure 9 .Figure 10 .

 8910 Figure 8. Density at times (a) 0.2 s, and (b) 1 s

Acknowledgements. L. Boudin wants to thank F. James and R. Eymard for the various discussions they shared about the topic of this article a few years ago.