Guillaume Barbu

Alberto Battistelllo

Guillaume Dabosville

Christophe Giraud

Guénaël Renault

Soline Renner

Rina Zeitoun

Alberto Battistello

• Guénaël

Renault 2⋆

Combined Attack on CRT-RSA. Why Public Verification Must Not Be Public?

Keywords: CRT-RSA, Combined Attacks, Fault Injection, Side-Channel Analysis, Coppersmith's methods

published or not. The documents may come

Introduction

Since the seminal work of Kocher published in 1996 [START_REF] Kocher | Timing Attacks on Implementations of Diffie-Hellman, RSA, DSS, and Other Systems[END_REF], Side-Channel Analysis (SCA) has raised a huge interest in both academic and industrial communities. This kind of attack is based on the fact that side-channel leakages of embedded devices contain information on the values manipulated inside the device. Therefore any sensitive variable carelessly used can be recovered by an attacker using SCA. Originally, time execution was used as side-channel leakage but the exploitation of power consumption and electromagnetic radiation became quickly the most efficient way to attack embedded cryptography [START_REF] Kocher | Differential Power Analysis[END_REF][START_REF] Quisquater | A New Tool for Non-intrusive Analysis of Smart Cards Based on Electro-magnetic Emissions, the SEMA and DEMA Methods[END_REF]. Over the years, many improvements have been made leading to very efficient attacks and very ingenious countermeasures [START_REF] Mangard | Power Analysis Attacks -Revealing the Secrets of Smartcards[END_REF].

In parallel to SCA, Fault Injection (FI) provides the attacker with another way to attack embedded devices. Such attacks aim at disturbing cryptographic computations and the analysis of corresponding faulty outputs allows the attacker to recover the secret key [START_REF] Giraud | A Survey on Fault Attacks[END_REF]. Shortly after the original publication focusing on RSA implementation [START_REF]Bellcore: New Threat Model Breaks Crypto Codes[END_REF], many other articles have been published to present FI attacks on various cryptosystems such as DES, ElGamal or DSA signature schemes [START_REF] Biham | Differential Fault Analysis of Secret Key Cryptosystem[END_REF][START_REF] Bao | Breaking Public Key Cryptosystems an Tamper Resistance Devices in the Presence of Transient Fault[END_REF]. As for SCA, FI has been deeply studied over the last decade [START_REF] Joye | Fault Analysis in Cryptography[END_REF] and the consequences of both attacks on the industry are huge since secure products must now be certified to prove their resistance against such threats.

Over the last few years, the cryptographic community has investigated the possibility of combining the two previous kinds of attacks. This has resulted in a new class of attacks called Combined Attacks (CA) that can defeat implementations which are meant to resist both SCA and FI. However, as far as we know only four CA have been published since their introduction in 2007, proving the difficulty to conceive such attacks [START_REF] Amiel | Passive and Active Combined Attacks -Combining Fault Attacks and Side Channel Analysis[END_REF][START_REF] Fang | To Infinity and Beyond: Combined Attack on ECC Using Points of Low Order[END_REF][START_REF] Roche | Combined Fault and Side-Channel Attack on Protected Implementations of AES[END_REF][START_REF] Robisson | Differential Behavioral Analysis[END_REF].

Nowadays, most embedded devices implement a large variety of cryptosystems to ensure the security of the sensitive assets they contain. As well as being the first practical public-key cryptosystem published, RSA [START_REF] Rivest | A Method for Obtaining Digital Signatures and Public-Key Cryptosystems[END_REF] has also been the most widely used for many years. In particular, the RSA using the Chinese Remainder Theorem (CRT), providing a speed-up factor of four compared to the original implementation, is available on most ID, banking and mobile smart cards. Obviously, this cryptosystem has been the main target of SCA and FI attackers leading to the development of efficient countermeasures having the smallest impact on both memory consumption and time execution due to the constraints of embedded environment.

In this article, we describe a new Combined Attack against a CRT-RSA implementation resistant to SCA by using blinding countermeasures and protected against FI by verifying the signature using the public exponent. Such an implementation is known to resist each and every kind of attack published so far. However, we demonstrate that when injecting a fault during the signature computation, a value depending on the message and on a multiple of a secret prime is manipulated in plain during the public verification. Therefore, we notice that one can use SCA to gain some information on such a sensitive value. The recovered information can then be used to factorize the RSA modulus and thus to reveal the whole private key. Besides, we exploit lattice techniques, and in particular Coppersmith's methods for finding small solutions to polynomial equations [START_REF] Coppersmith | Finding a small root of a univariate modular equation[END_REF][START_REF] Coppersmith | Small Solutions to Polynomial Equations, and Low Exponent RSA Vulnerabilities[END_REF], to significantly reduce the complexity of our CA.

The rest of this paper is organised as follows. In Section 2 we briefly recall some basics on CRT-RSA signature as well as the corresponding attacks and countermeasures. In Section 3 we describe our new CA on a CRT-RSA implementation that is known to resist both SCA and FI attacks. In Section 4 we present the results of our simulations which prove the efficiency of our new attack. We then improve its complexity by using lattice reduction techniques in Section 5. Finally, we suggest in Section 6 possible countermeasures having a negligible penalty on the performance of the algorithm.

Previous Works

In this section we briefly recall the RSA signature, in particular the CRT mode. Secondly we present the principal attacks on such an algorithm as well as the main countermeasures.

RSA on Embedded Systems

Since its introduction in 1978, the RSA cryptosystem has become one of the most used public-key cryptosystems, especially in electronic signature schemes [START_REF] Rivest | A Method for Obtaining Digital Signatures and Public-Key Cryptosystems[END_REF].

In the following we briefly recall how to compute the RSA signature in both standard and CRT modes.

Let N denote the public modulus being the product of two secret large prime integers p and q. Let d refer to the private exponent and e refer to the public exponent satisfying de = 1 mod ϕ(N ), where ϕ denotes Euler's totient function. The RSA signature of a message m ∈ Z N is then obtained by computing S = m d mod N . To verify the signature, one computes S e mod N and checks if the corresponding result is equal to m.

In embedded systems, most RSA implementations use the Chinese Remainder Theorem (CRT) which yields an expected speed-up factor of four [START_REF] Couvreur | Fast Decipherment Algorithm for RSA Public-Key Cryptosystem[END_REF]. Following the CRT-RSA algorithm, the signature generation is composed of two exponentiations S p = m dp mod p and S q = m dq mod q, where d p = d mod p -1 and d q = d mod q -1. The signature is then obtained by recombining S p and S q , which is usually done by using Garner's formula [START_REF] Garner | The Residue Number System[END_REF]:

S = CRT (S p , S q ) = S q + q(i q (S p -S q ) mod p) , (1)
where i q = q -1 mod p.

Attacks and Countermeasures

Side-Channel Analysis Side-Channel Analysis (SCA) has been introduced by the publication of the so-called timing attacks in 1996 [START_REF] Kocher | Timing Attacks on Implementations of Diffie-Hellman, RSA, DSS, and Other Systems[END_REF]. SCA exploits the dependency between the manipulated data or the executed instruction and the side-channel leakages which can be monitored during the algorithm execution.

Examples of such leakages are the power consumption or the electromagnetic radiation of the device. When only one measure is required to exploit sensitive information, the attack is called Simple Passive Analysis (SPA) [START_REF] Kocher | Differential Power Analysis[END_REF]. In the case of a straightforward Square-and-Multiply exponentiation, SPA consists for instance in observing if the squaring and multiplication operations have different patterns in the corresponding side-channel leakages [START_REF] Kocher | Differential Power Analysis[END_REF]. Hence, the secret exponent can be directly extracted from one measurement. In the literature, a common countermeasures consists in using a so-called regular algorithm which performs the same operation whatever the exponent bit value such as the Square-Always or Montgomery ladder algorithms [START_REF] Joye | The Montgomery Powering Ladder[END_REF][START_REF] Clavier | Square Always Exponentiation[END_REF]. Moreover, attacks based on side-channel leakages have evolved to a type of SCA called Differential Passive Analysis (DPA) [START_REF] Kocher | Differential Power Analysis[END_REF] which requires a large number of measurements. This type of attack applies a statistical treatment on the curves to recover information on the manipulated values.Nowadays a common statistic tool used to perform such a statistical treatment, is the Pearson correlation coefficient:

ρ k = cov(L, H) σ L σ H , (2)
where L is the set of curves and H depends on a known value m and on a guess of a small part of a secret k. Such an attack is called Correlation Power Analysis (CPA) [START_REF] Brier | Correlation Power Analysis with a Leakage Model[END_REF]. In the literature, many different CPAs have been published to attack the RSA cryptosystem [START_REF] Amiel | Power Analysis for Secret Recovering and Reverse Engineering of Public Key Algorithms[END_REF]. For instance in the CRT-mode, an attacker can mount a CPA to recover the private parameter q during the CRT-recombination, cf. Rel. [START_REF] Kocher | Timing Attacks on Implementations of Diffie-Hellman, RSA, DSS, and Other Systems[END_REF], by observing the leakage obtained during the manipulation of the value i q (S p -S q ) mod p and by making a guess on a few bits of q. Hence, an attacker can obtain the whole secret q by performing a CPA for each of its subpart (typically for each byte). However, half of q is sufficient to recover the rest of q by using Coppersmith's attack [START_REF] Coppersmith | Finding a small root of a bivariate integer equation; factoring with high bits known[END_REF]. Classical countermeasures to resist CPA consist in randomizing the modulus, the message and the exponent [START_REF] Amiel | Power Analysis for Secret Recovering and Reverse Engineering of Public Key Algorithms[END_REF].

Fault Injection RSA has been the first cryptosystem to succumb to Fault Injection [START_REF] Boneh | On the Importance of Checking Cryptographic Protocols for Faults[END_REF]. In the following, we describe such an attack in the CRT case. Assume that a fault is injected during the computation of S p leading to a faulty signature S. Since S ≡ S p mod p and S ≡ S q mod q, one can notice that S ≡ S mod q but S ≡ S mod p. Therefore, the secret parameter q can be easily recovered by computing the gcd of S -S and N . The rest of the private key can then be straightforwardly deduced. When it is not possible to sign the same message twice and if the message is known to the attacker, a variant of this attack consists in computing the gcd of S em and N to obtain the secret value q [START_REF] Lenstra | Memo on RSA Signature Generation in the Presence of Faults[END_REF].

Moreover, the effect of fault injections on CRT-RSA is not limited to the disturbance of S p or S q . Indeed, a fault injected in any part of the key parameters (i.e. p, q, d p , d q or i q ), in the message m at the beginning of either S p or S q computation, or even during the CRT-recombination can lead to a useful faulty signature.

In the literature, four different fault models are generally considered to define the attacker's capabilities [START_REF] Blömer | A New RSA-CRT Algorithm Secure against Bellcore Attacks[END_REF]:

the random fault model: the bits are changed to a uniformly distributed random value; -the bit-flip fault model: in that case, affected bits are flipped to their complementary value; -the stuck-at fault model: the fault sets the bits to 0 or to 1, depending on the underlying hardware; -the unknown constant fault model: the fault always sets the bits to the same unknown value.

Moreover, these faults do not necessarily modify a whole temporary result. Indeed, it is generally considered that the number of bits affected by the fault is linked to the CPU word-size which is generally 8, 16 or 32 bits.

The most natural way to counteract fault injection on RSA-type signature is to check the correctness of the signature S before outputting it [START_REF] Boneh | On the Importance of Checking Cryptographic Protocols for Faults[END_REF]. More precisely, the signature is returned iff S e mod N = m. Moreover, such a method requires very little overhead since the public exponent e is usually small in practice (typically 3, 17 or 2 16 + 1).

Other methods getting rid of e have also been proposed but they do not offer the same level of security and are generally slower than the public verification [START_REF] Giraud | An RSA Implementation Resistant to Fault Attacks and to Simple Power Analysis[END_REF][START_REF] Vigilant | RSA with CRT: A New Cost-Effective Solution to Thwart Fault Attacks[END_REF][START_REF] Rivain | Securing RSA against Fault Analysis by Double Addition Chain Exponentiation[END_REF].

Combined Attacks

The idea to combine SCA and FI appeared in 2007 when Amiel et al. proposed a so-called Combined Attack (CA) on an RSA implementation protected against FI and SPA [START_REF] Amiel | Passive and Active Combined Attacks -Combining Fault Attacks and Side Channel Analysis[END_REF]. They noticed that by setting to zero one of the temporary registers used in the Montgomery ladder, its structure becomes unbalanced, revealing the value of the secret exponent by SPA. Following this publication, three other papers have been published taking advantage of this new way of defeating embedded security. Two of them present a CA against a secured AES implementation [START_REF] Roche | Combined Fault and Side-Channel Attack on Protected Implementations of AES[END_REF][START_REF] Robisson | Differential Behavioral Analysis[END_REF]. The third one focuses on the elliptic curve scalar multiplication [START_REF] Fang | To Infinity and Beyond: Combined Attack on ECC Using Points of Low Order[END_REF].

Despite its theoretical effectiveness, the combination of SCA and FI is very difficult in practice, explaining the lack of practical experiments in the current literature.

Lattices Randomized RSA encoding schemes are usually considered to be resistant to traditional FI attacks since a part of the message is unknown to the attacker and varies for each signature computation. However this common assumption has to be mitigated regarding the works of [START_REF] Coron | Fault Attacks on RSA Signatures with Partially Unknown Messages[END_REF] and [START_REF] Coron | Fault Attacks against EMV Signatures[END_REF] which defeat two randomised RSA encoding schemes. These attacks use Coppersmiths method to solve a bivariate polynomial whose coefficients are built thanks to the generated faulty signatures.

More recently in [START_REF] Brier | Modulus Fault Attack against RSA-CRT Signatures[END_REF], the authors present an attack taking advantage of the disturbance of the public modulus. The generated faulty signatures allow them to build a lattice, which in turns leads to factorize the public modulus.

Although [START_REF] Coron | Fault Attacks on RSA Signatures with Partially Unknown Messages[END_REF], [START_REF] Coron | Fault Attacks against EMV Signatures[END_REF] and [START_REF] Brier | Modulus Fault Attack against RSA-CRT Signatures[END_REF] also apply lattice reduction techniques as carried out in this paper, the attack presented hereafter does not share the same context since we consider a secured implementation which never returns the faulty signature.

3 A New Combined Attack on CRT-RSA

Context and Principle

As stated in Section 2, several countermeasures have been developed to protect CRT-RSA embedded implementations against both SCA and FI. In the framework of this article, we consider an algorithm protected:

against SCA by using message and exponent blinding as suggested in [START_REF] Witteman | Defeating RSA Multiply-Always and Message Blinding Countermeasures[END_REF],

a regular exponentiation algorithm such as the Square Always [START_REF] Clavier | Square Always Exponentiation[END_REF] and a mask refreshing method along the exponentiation such as the one presented in [START_REF] Dupaquis | Redundant Modular Reduction Algorithms[END_REF]. Moreover, the blinding is kept all along the CRT-recombination. -against FI by verifying the signature using the public exponent e [START_REF] Boneh | On the Importance of Checking Cryptographic Protocols for Faults[END_REF]. In addition, we also use the approach presented in [START_REF] Dottax | On Second-Order Fault Analysis Resistance for CRT-RSA Implementations[END_REF] which mainly consists in checking the result of the verification twice to counteract double FI attacks.

Fig. 1 depicts the main steps of such an implementation where the k i 's are random values (typically of 64 bits) generated at each execution of the algorithm and S ′ p , S ′ q and S ′ represent the blinded version of S p , S q and S respectively.

In the following, we assume that the fault injected by the attacker follows either the bit-fault, the stuck-at or the unknown constant fault models (cf. Section 2.2). Moreover, we assume the attacker is able to choose which byte of the message is affected by the fault.

As mentioned in Section 2.2, injecting a fault during the signature computation leads to a faulty signature that allows the attacker to recover the private key. However in the implementation considered in this paper, the verification with the public exponent detects such a disturbance and the faulty signature is never revealed to the attacker. The main contribution of this paper is to show that in this case, an SCA can still allow the attacker to gain enough information on the faulty signature to recover the private key.

At first glance, it seems impossible to perform such an attack during the signature process due to the blinding countermeasure. However by observing Fig. 1, one may note that the faulty signature S remains blinded until the end of exponentiation with e modulo N . Therefore if we can express S e mod N in terms of the message m and of the private key then we can perform an SCA on this value. In the next section, we exhibit such a relation allowing us to mount a CA on an SCA-FI-resistant CRT-RSA implementation.

A Useful Relation

Proposition 1. If a fault ε is induced in m such that the faulty message m is equal to m + ε at the very beginning of the computation of S p then

S e ≡ m + εqi q mod N , (3)
where S corresponds to the faulty signature.

Proof. By definition of the CRT-RSA signature, we have:

S ≡ (m + ε) d mod p S ≡ m d mod q (4)
It comes then straightforwardly that:

S e ≡ m + ε mod p S e ≡ m mod q (5)
Finally, applying Gauss recombination to (5) leads to (3) since:

S e ≡ pi p m + qi q (m + ε) mod N (6) ≡ (pi p m + qi q m) + εqi q mod N (7) ≡ m + εqi q mod N , (8)
where i p = p -1 mod q. ⊓ ⊔

One may note that a similar relation holds if m is disturbed at the very beginning of S q computation due to the symmetrical roles of p and q in both branches of the CRT-RSA. For the sake of simplicity, we will use the case where S p computation is disturbed in the rest of this paper.

Recovering the Private Key

Following the attack's principle depicted in Section 3.1 and using Proposition 1, we will now present in detail the main steps of our attack.

Firstly, the attacker asks the embedded device to sign several messages m i through a CRT-RSA implemented as described in Section 3.1. For each signature, the computation of S q is performed correctly and a constant additive error ε is injected on the message m i at the beginning of each S p computation. Then during each signature verification, the attacker monitors the corresponding side-channel leakage L i which represents the manipulation of S e i mod N . From Proposition 1, we know that there exists a sensitive value k satisfying the relation S e i mod N = m i + k. Therefore, the attacker will perform a CPA to recover this sensitive value by computing ρ k (m i + k, L i ) for all the possible values of k (cf. Section 2.2).

Depending on the set {(m i , S e i mod N )} i , it follows from Rel. (3) that k will be equal either to εqi q mod N or to εqi q mod N -N . Therefore, the value k producing the strongest correlation at the end of the CPA will be one of these two values. Once k recovered, the attacker must then compute the gcd between k and N , which leads to the disclosure of q. From this value, the private key is straightforwardly computed.

Regarding the practicality of our fault model (i.e. a constant additive fault), one may note that by fixing a small part of the message (e.g. a byte), the disturbance of such a part in either the stuck-at, the bit-flip or the unknown constant fault model results in a constant additive error during the different signature computations. Therefore our fault model is definitely valid if the attacker can choose the messages to sign, or even if she can only have the knowledge of the messages and attack only those with a given common part.

Finally, one may note that it is not possible to perform a statistical attack targeting the full value of k at once due to its large size (i.e. ⌈log 2 (N )⌉ bits). However, one can attack each subpart of this value, for instance by attacking byte per byte starting with the least significant one in order to be able to propagate easily the carry. It is worth noticing that CPA only applies when the corresponding part of the message varies. Therefore, if the attacker fixes the MSB of the message, then the corresponding set of measurements can be used to recover the whole but last byte of k. In such a case, a brute force search can be used to recover the missing byte.

In the next section, we present simulations of our attack which prove the efficiency of our method and which are based on the attacker's capability to inject the same fault and on the noise of the side-channel measurements.

Experiments

The success of the attack presented in Section 3 relies on the ability of the attacker to both measure the side-channel leakage of the system during the signature verification and induce the same fault ε on the different manipulated messages.

In order to evaluate the effectiveness of this attack, we have experimented it on simulated curves of the side-channel leakage L, according to the following leakage model:

L(d) = HW (d) + N (µ, σ) (9)
with N (µ, σ) a Gaussian noise of mean µ and standard deviation σ, and HW (d) the Hamming weight function evaluated for the manipulated data d. In the framework of our experiments, we consider that the processor manipulates 8-bit words and we use three different levels of noise, namely σ = 0.1, 1 and 5.

As well as the side-channel leakage, the faults were also simulated by setting the most significant word of the message m to all-0 at the very beginning of the S p computation. These faults were induced with a given success rate r, varying in our different experiment campaigns (namely 50%, 10% and 1%).

Depending on the experimental settings, all the different words of the secret value will be equivalently correlated with the simulated curves. The graphs presented in Fig. 2 present the convergence of the correlation for each possible value k of one particular byte (the 5 th least-significant byte) of the secret depending on the number of side-channel measurements with different simulation settings σ and r.

As exposed in Fig. 2, the number of traces required to recover the secret value depends essentially on the fault injection success rate. This comes from the fact that every wrongly-faulted computation can be considered as noise in the scope of our statistical analysis. The number of curves required to retrieve the secret word grows as the fault injection success rate decreases and to a fewer extent as the noise of the side-channel leakage increases.

With regards to the results obtained when σ = 5 and r = 10%, which appear to be plausible values in practice, it took us 3.35 seconds to retrieve one byte of the secret value by performing the CPA on 15, 000 curves of 128 points each 4 . Assuming a genuine curve should be made of at least 5, 000 points, we can estimate the time required to practically perform the attack to about 1 minute 5 seconds per byte. That is to say, it takes about 2 hours 20 minutes to recover the complete secret value if we consider a 1024-bit RSA module.

For the sake of clarity, we restrained the experiments presented here to the case where the processor manipulates 8-bit words, and thus ε is an 8-bit error. The same experiments have been run for processor word-size up to 32 bits with success. Besides, about the same number of curves were necessary for the CPA to highlight the correct secret byte.

Section 5 shows how it is possible to considerably reduce the complexity of our attack thanks to the use of lattice techniques.

Reducing the Attack Complexity Using Coppersmith's Methods

This section aims at improving the attack complexity using Coppersmith's methods. It is in line with the problem of factorizing N knowing half part of prime p (or q), that was solved in [START_REF] Coppersmith | Small Solutions to Polynomial Equations, and Low Exponent RSA Vulnerabilities[END_REF]. With respect to our case, we highlight that if the CA presented in Section 3 provides about half of the secret εqi q mod N , then the other half part can be straightforwardly computed by solving a well-designed modular polynomial equation that we elaborate in the sequel. Besides, we deal with two cases (ε known and unknown), depending on the fault model that is considered.

Bringing Up the Original Problem to Solving a Modular Equation

Suppose we are given the t least significant bits (LSB) of the secret εqi q mod N . The latter value can be rewritten as follows:

εqi q ≡ 2 t x 0 + k mod N , (10)
where t and k are known values, and x 0 is the ⌈log 2 (N ) -t⌉-bit unknown integer that is to be recovered.

Lemma 1. The unknown secret part x 0 is solution of the polynomial P ε (x):

P ε (x) = x 2 + c (2 t+1 k -2 t ε) x + c (k 2 -kε) ≡ 0 mod N , (11)
where c = (2 2t ) -1 mod N , k, t, N are known, and ε is the induced fault.

Proof. The Bézout identity applied to our context yields that primes p and q interrelate with integers i p = p -1 mod q and i q = q -1 mod p by the following relation:

pi p + qi q ≡ 1 mod N . (12)
Multiplying ( 12) by ε leads to the relation εpi p +εqi q ≡ ε mod N , or equivalently to εpi p ≡ ε-εqi q mod N . Therefore, replacing εqi q using (10) allows us to deduce an equivalence for εpi p :

εpi p ≡ ε -2 t x 0 -k mod N . (13)
As N = pq, we then multiply (10) by [START_REF] Robisson | Differential Behavioral Analysis[END_REF], to get the relation:

εqi q • εpi p ≡ (2 t x 0 + k) • (ε -2 t x 0 -k) ≡ 0 mod N . (14)
Eventually, developing the right-hand side of ( 14), and multiplying it by c = (2 2t ) -1 mod N leads to the obtention of the monic polynomial P ε (x).

⊓ ⊔

The initial problem of retrieving the unknown part of εqi q mod N is thereby altered in solving the modular polynomial equation [START_REF] Fang | To Infinity and Beyond: Combined Attack on ECC Using Points of Low Order[END_REF]. In the sequel, we deal with two possible cases regarding ε, whether it is known to the adversary or not.

Case 1. The fault ε is known to the adversary This case corresponds to the bit-flip and stuck-at fault models (Section 2.2) since the message is known to the attacker and the fault location can be chosen. In both cases, since the fault ε is known, the problem is reduced to solving a univariate modular polynomial equation, cf. Rel. ((11)). This problem is known to be hard. However, when the integer solution x is small, Coppersmith showed [START_REF] Coppersmith | Finding a small root of a bivariate integer equation; factoring with high bits known[END_REF] that it can be retrieved using the well-known LLL algorithm. Accordingly, we induce the following proposition: Proposition 2. Given N = pq and the low order 1/2 log 2 (N ) bits of εqi q mod N and assuming ε is known, one can recover in time polynomial in (log 2 (N ), d) the factorization of N .

Proof. From Coppersmith's Theorem [START_REF] Coppersmith | Small Solutions to Polynomial Equations, and Low Exponent RSA Vulnerabilities[END_REF], we know that, given a monic polynomial P (x) of degree d, modulo an integer N of unknown factorization, and an upper bound X on the desired solution x 0 , one can find in polynomial time all integers x 0 such that

P (x 0 ) ≡ 0 mod N and |X| < N 1/d . (15)
In our case we have d = 2, and since x 0 is a ⌈log 2 (N ) -t⌉-bit integer, we know that |x 0 | < X = 2 ⌈log 2 (N )-t⌉ . Thus, the condition in (15) becomes 2 ⌈log 2 (N )-t⌉ < N 1/2 , i.e.

t > 1 2 log 2 (N ) . (16)
Therefore, knowing at least half part of the secret εqi q mod N allows to recover the whole secret. As previously done, computing gcd(εqi q mod N, N ) provides the factorization of N .

⊓ ⊔

Note that the method is deterministic, and as will be seen further (Table 1), it is reasonably fast.

Case 2. The fault ε is unknown to the adversary This case is met in the unknown constant fault model (Section 2.2). In such a case, one can consider the polynomial P ε (x) as a bivariate modular polynomial equation with unknown values x and ε. This specific scheme has also been studied by Coppersmith and includes an additional difficulty of algebraic dependency of vectors which induces the heuristic characteristic of the method [START_REF] Coppersmith | Finding a small root of a univariate modular equation[END_REF]. As depicted in Section 5.2, in our experiments nearly 100% of the tests verified the favorable property of independency. Accordingly, in this vast majority of cases, the following proposition holds: Proposition 3. Under an hypothesis of independency (see discussion above), given N = pq and the low order 1/2 log 2 (N ) + s bits of εqi q mod N , where s denotes the bitsize of ε, and assuming ε is unknown, one can recover in time polynomial in (log 2 (N ), d) the factorization of N .

Proof. Coppersmith's Theorem for the bivariate modular case [START_REF] Coppersmith | Finding a small root of a univariate modular equation[END_REF] notifies that given a polynomial P (x, ε) of total degree d, modulo an integer N of unknown factorization, and upper bounds X and E on the desired solutions x 0 , ε 0 , it may be possible (heuristic) to find in polynomial time all integer couples (x 0 , ε 0 ) such that P (x 0 , ε 0 ) ≡ 0 mod N and

|X • E| < N 1/d . (17)
In our case, we have d = 2 and E = 2 s . The integer x 0 is ⌈log 2 (N ) -t⌉-bit long, therefore we have X = 2 ⌈log 2 (N )-t⌉ . Thus, the condition in ( 17) becomes

2 ⌈log 2 (N )-t⌉ • 2 s < N 1/2 , i.e. t > 1 2 log 2 (N ) + s . (18)
This means that knowing s more bits of the secret εqi q mod N than before, would allow the recovering of the whole secret. ⊓ ⊔ Remark 1. The bound of success in Proposition 3 can actually be slightly improved using results of [START_REF] Blomer | A Tool Kit for Finding Small Roots of Bivariate Polynomials over the Integers[END_REF]. Indeed, Coppersmith's bound applies to polynomials whose monomials shape is rectangular, while in our case the monomial ε 2 does not appear in P (x, ε) which corresponds to what they called an extended rectangle in [START_REF] Blomer | A Tool Kit for Finding Small Roots of Bivariate Polynomials over the Integers[END_REF]. For the sake of simplicity, we only mentioned Coppersmith's bound since practical results are similar.

Results From Our Implementation

We have implemented this lattice-based improvement using Magma Software [START_REF] Bosma | The Magma algebra system. I. The user language[END_REF], with N a 1024-bit integer i.e. 128 bytes long, in the cases where ε is an 8bit known value (for Case 1) and a 32-bit unknown value (for Case 2). We chose Howgrave-Graham's method [START_REF] Howgrave-Graham | Finding small roots of univariate modular equations revisited[END_REF] for the univariate case, and its generalization by Jochemsz et al. [START_REF] Jochemsz | A strategy for finding roots of multivariate polynomials with new applications in attacking rsa variants[END_REF] for the bivariate case since both have the same bound of success as Coppersmith's method (sometimes even better for [START_REF] Jochemsz | A strategy for finding roots of multivariate polynomials with new applications in attacking rsa variants[END_REF]) and they are easier to implement. As we know, the theoretical bound given in Coppersmith's method is only asymptotic [START_REF] Coppersmith | Small Solutions to Polynomial Equations, and Low Exponent RSA Vulnerabilities[END_REF]. Thus, we report in Table 1 (for Case 1) and in Table 2 (for Case 2) the size t (in bytes) of the secret εqi q mod N that is known to the attacker before applying Coppersmith's method, the lattice dimension used to solve [START_REF] Fang | To Infinity and Beyond: Combined Attack on ECC Using Points of Low Order[END_REF] and finally the timings of our attack.

As depicted in Table 1, and combining these results with the experiments of Section 4, the best trade-off is to perform a CPA on the 66 first bytes, taking 66 × 1m05s = 1h11m30s, and to retrieve the 62 remaining bytes using lattices in 34.25s, bringing the total time up to 1 hour 12 minutes, instead of the previous 2 hours 20 minutes. In order to illustrate Case 2, we have chosen to rather show our results for ε being a 32-bit value, since when ε is 8-bit long, we obtained slightly better results by considering the 255 possible values of the variable ε together with their corresponding polynomials P ε (x), and by running the method on each of the polynomials until finding the solution x 0 that allows to factorize N. This

Countermeasures

In this section, we describe different countermeasures to protect an implementation against the CA presented in Section 3.

Blind Before Splitting

Our first proposition consists in avoiding the possibility to inject the same fault during several signature computations. To do so, we deport the blinding of the input message m before executing the two exponentiations modulo p and q:

m ′ = m + k 0 N mod k 1 N , (19)
with k 0 and k 1 two n-bit random values generated at each algorithm execution (n being typically 64). Hence S ′ p = m ′dp mod k 2 p and S ′ q = m ′dq mod k 3 q. This countermeasure prevents an attacker from injecting always the same error during the signature computation. Indeed if the fault is injected on m at the very beginning of one exponentiation, then the corresponding error cannot be fixed due to the blinding injected by Rel. [START_REF] Joye | The Montgomery Powering Ladder[END_REF].

Moreover, if the fault is injected when the message m is manipulated during [START_REF] Joye | The Montgomery Powering Ladder[END_REF], then the error ε impacts the computation of both S ′ p and S ′ q , leading to inexploitable faulty outputs.

Such a countermeasure induces a small overhead in terms of memory space since m ′ must be kept in memory during the first exponentiation but the execution time remains the same.

Verification Blinding

Our second countermeasure aims at annihilating the second hypothesis of our attack: a predictive variable is manipulated in plain during the verification. To do so, we inject a ⌈log 2 (N )⌉-bit random r before performing the final reduction with N , cf. Rel. [START_REF] Clavier | Square Always Exponentiation[END_REF]. Therefore, each and every variable manipulated during the verification is blinded.

(( S e + r -m) mod k 1 N ) mod N ? = r . (20)
One may note that the final comparison should be performed securely with regards to the attack described in [START_REF] Lomne | On the Need of Randomness in Fault Attack Countermeasures -Application to AES[END_REF] since information on εqi q could leak if such a comparison was performed through a substraction.

The cost of such a countermeasure is negligible since it mainly consists in generating a ⌈log 2 (N )⌉-bit random variable.

Conclusion

This paper introduces a new Combined Attack on CRT-RSA. Even if a secure implementation does not return the faulty signature when the computation is disturbed, we show how to combine FI with SCA during the verification process to obtain information on the faulty signature. Such information allows us to factorize the public modulus and thus to recover the whole private key. We also show that Coppersmith's methods to solve univariate and bivariate modular polynomial equations can be used to significantly reduce the complexity of our new attack. Finally, we provide simulations to confirm the efficiency of our method and we present two countermeasures which have a very small penalty on the performance of the algorithm.

Our main objective was to prove that stacking several countermeasures does not provide global security despite addressing each and every attack separately. Therefore, the main consequence of this paper is that fault injection countermeasures must also be designed to resist SCA and vice versa.

Fig. 1 .

 1 Fig. 1. Main steps of a CRT-RSA implementation secure against SCA and FI.

Fig. 2 .

 2 Fig.2. Convergence of the correlation for the 256 possible values ki for the secret (the correct one being depicted in black) depending on the number of side-channel measurements (×500) for different levels of noise σ and fault injection success rates r.

Table 1 .

 1 Size t required (in bytes) for the method to work and timings (Magma V2.17-1), as a function of the lattice dimension in Case 1 (ε known, being an 8-bit integer).

	Size t required 86 (bytes)	72	70	69	68	67	66	65	64
	Dimension of the lattice	3	9	11	15	17	23	37	73 Theoretical
	Time for LLL < 0.01 0.03 0.07 0.29 0.52 2.63 34.25 2587.7 bound (seconds)

Table 2 .

 2 Size t required (in bytes) for the method to work and timings (Magma V2.17-1), as a function of the lattice dimension in Case 2 (ε unknown, being a 32-bit integer).indeed leads to a best trade-off of 70 bytes required from the CPA and the 58 remaining bytes computed with lattices by performing 255 times the LLL algorithm in the worst case, for a total of 68 × 1m05s + 255 × 0.52s, i.e. 1 hour 16 minutes instead of 2 hours 20 minutes. Besides, this exhaustive search can be performed in parallel and it also has the advantage to be deterministic. However, when ε is 32-bit long, an exhaustive search becomes impractical and, as depicted in Table2, the best trade-off would be to perform a CPA on 72 bytes and to compute the 56 remaining bytes with lattices (even if heuristic, it worked in nearly 100% of the tests in practice), resulting in a total of 72×1m05s+30.22s, i.e. 1 hour 18 minutes instead of the previous 2 hours 20 minutes.

	Size t required 86 (bytes)	78	76	74	73	72	71	70	69
	Dimension of the lattice	5	12	22	35	51	70	117	201 Theoretical
	Time for LLL < 0.01 0.02 0.16 1.17 5.88 30.22 605.9 12071.1 bound (seconds)

The execution time given here and in Section

[START_REF] Giraud | A Survey on Fault Attacks[END_REF].2 have been obtained on a 32-bit CPU @3.2GHz

Acknowledgements. The authors would like to thank Emmanuelle Dottax for her useful comments on the preliminary version of this article. We would also like to thank Jean-Sébastien Coron and the anonymous reviewers of PKC'13 for their valuable comments and suggestions.

⋆ This work was partly supported by the French ANR under the Computer Algebra and Cryptography (CAC) project (ANR-09-JCJCJ-0064-01).


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

Jules Svartz
email: jules.svartz@lip6.fr

Solving Polynomial Systems Globally Invariant Under an Action of the Symmetric Group and Application to the Equilibria of N vortices in the Plane

Keywords: I, 1, 2 [Computing Methodologies]: Symbolic and Algebraic Ma-

We propose an efficient algorithm to solve polynomial systems of which equations are globally invariant under an action of the symmetric group S N acting on the variable x i with σ (x i ) = x σ (i) and the number of variables is a multiple of N. For instance, we can assume that swapping two variables (or two pairs of variables) in one equation gives rise to another equation of the system (perhaps changing the sign). The idea is to apply many times divided difference operators to the original system in order to obtain a new system of equations involving only the symmetric functions of a subset of the variables. The next step is to solve the system using Gröbner techniques; this is usually several order faster than computing the Gröbner basis of the original system since the number of solutions of the corresponding ideal, which is always finite has been divided by at least N!.

To illustrate the algorithm and to demonstrate its efficiency, we apply the method to a well known physical problem called equilibria positions of vortices. This problem has been studied for almost 150 years and goes back to works by von Helmholtz and Lord Kelvin. Assuming that all vortices have same vorticity, the problem can be reformulated as a system of polynomial equations invariant under an action of S N . Using numerical methods, physicists have been able to compute solutions up to N ≤ 7 but it was an open challenge to check whether the set of solution is complete. Direct naive approach of Gröbner bases techniques give rise to hard-to-solve polynomial system: for instance, when N = 5, it takes several days to compute the Gröbner basis and the number of solutions is 2060. By contrast, applying the new algorithm to the same problem gives rise to a system of 17 solutions that can be solved in less than 0.1 sec. Moreover, we are able to compute all equilibria when N ≤ 7.

INTRODUCTION

Solving general polynomial systems is a fundamental problem in Computer Algebra. However it must be emphasized that polynomial systems arising in applications are very often highly structured. For instance, in several algebraic problems coming from applications the solutions (the algebraic variety) is invariant under the action of a finite group: for example the Cyclic-n problem [START_REF] Björck | Functions of modulus 1 on Zn , whose Fourier transforms have constant modulus, and "cyclic n-roots[END_REF], in Cryptography the NTRU Cryptosystem [START_REF] Hoffstein | NTRU: a ring-based public key cryptosystem[END_REF] or the membrane inclusions curvature equations in biology [START_REF] Faugère | The membrane inclusions curvature equations[END_REF]. Hence an important subproblem is to solve efficiently such algebraic problems. We should consider two distinct cases. First, if all the equations are invariant under the action of the group, there are two ways to solve the system using the symmetries. In [START_REF] Colin | Solving a system of algebraic equations with symmetries[END_REF] Colin propose to use invariants [START_REF] Sturmfels | Algorithms in invariant theory[END_REF] to solve the system. This method is very efficient if the Hironaka Decomposition of the ring of invariants is simple, but for the Cyclic-n problem, for example, it seems better to use a second method based on SAGBI Gröbner Basis techniques [START_REF] Faugère | Solving systems of polynomial equations with symmetries using SAGBI-Gröbner bases[END_REF]. The second class of problems, which is probably the most important in practice, is to consider polynomial systems of which the set of solutions is globally invariant under the action of a finite group (this is the case of the biology example [START_REF] Faugère | The membrane inclusions curvature equations[END_REF]). The goal of the present paper is to propose an efficient method to solve such problems assuming that the group is the whole symmetric group. To illustrate the algorithm and to demonstrate its efficiency, we apply the method to a well known physical problem called equilibria positions of vortices.

The problem of finding and classifying all relative equilibria of N-point vortices in the plane is of long-standing interest. In the plane, attacks on the problem date back to the 1800s with the works of von Helmholtz [START_REF] Helmholtz | Über Integrale der hydrodynamischen Gleichungen, welche den Wirbelbewegungen entsprechen[END_REF] and later in the works by Thomson [START_REF] Lord | On vortex atoms[END_REF] (the later Lord Kelvin). A complete bibliography of papers on the subject can be found in [START_REF] Meleshko | A bibliography of vortex dynamics 1858-1956[END_REF] or [START_REF] Aref | Vortex crystals[END_REF]. Several families of equilibria have been found [START_REF] Aref | Vortex crystals[END_REF] and other solutions have been found numerically, see [START_REF] Dirksen | Close pairs of relative equilibria for identical point vortices[END_REF]. More generally, the problem of equilibria on manifolds with different potentials has been studied by Albouy [START_REF] Albouy | The symmetric central configurations of four equal masses[END_REF].

In the planar case, the problem is equivalent to solve the following algebraic system (in the following Z symbolizes the complex conjugate of z):

Z i = N ∑ j=1, j =i 1 z i -z j
. In this paper we describe a general algorithm and for each step we apply it to the equilibria of N-point vortices. The proposed algorithms is a 3 steps process: 1. We apply many times divided difference operators (see sec-tion 3) to the original system in order to obtain a new system of equations involving only invariant equations. For instance, the 4-vortex problem is equivalent to r 0 = s 1 = r 1 -6 = r 2 = 2 r 3 + 5 s 2 = 0 where r k = ∑ i Z i z k i and s k = ∑ i z k i is the Newton sum.

2. As explained in section 4, the second step is to eliminate all the variables but the z i . For that purpose, we require that the algebraic system fulfils the parameterization assumption (see definition 6). We derive a new system of equations involving only the symmetric functions of a subset of the variables. For instance, for the 4-vortex problem we obtain the symmetric system e 3 (e 2 2 +12e 4 ) 2 = e 2 (e 4 2 -16e 2 2 e 4 +9e 2 e 2 3 +48e 2 4 ) = 0. 3. The last step consists simply in solving the symmetric equations using standard Gröbner bases techniques.

The first step can be viewed as generalisation of [START_REF] Faugère | The membrane inclusions curvature equations[END_REF], we give algorithms to obtain invariant equations under the action of S N on several blocks of variables.

Applied to the vortex problem, our method has three advantages over previous method:

1. In theory, it is possible to solve directly the original equations. However, when N = 5, it take several days to compute the Gröbner basis and the number of solutions is 2060. By contrast, applying the new algorithm to the same problem give rise to a system of 17 solutions that can be solved in less than 0.1 sec. The case N = 7 can be completely solved in about 20 minutes.

2. We are sure to find all the solutions, so we give a certificate for the previous numerical solutions. For N ≥ 5, it is completely new.

3. Two distinct solutions could be so close, that 300 digits are needed to be sure that they are distinct, see [START_REF] Dirksen | Close pairs of relative equilibria for identical point vortices[END_REF] for example. With exact computations, the solutions appear to be distinct without further computations.

Since we are using only exact computations, our algorithms gives computational proofs of the solutions of the vortex problem.

VORTEX PROBLEM

Physical equations and first steps

We start with the equations of motion for the N body problem:

∂ 2 ∂t 2 r i = ∑ j =i m j U ′ (s i j )(r i -r j ) for i = 1, . . . , N (1)
where m i and r i , are, respectively the mass and the position vector (relative to the center of mass) for the i-th particle and s i j = |r ir j | 2 is the square of the distance between particle i and j; U(s) is the potential function such that U ′ (s) = s a for some real value. Without lost of generality we can assume that the center of mass is at the origin. Usually the potential is one of the two well known potential:

a potential Newton -3/2 U(r) = r -1 2 Vortex -1 U(r) = log(r)
In this paper we assume that: we are in the planar case, all the masses (vortices) are equal (that is to say m i = 1) and that the potential is the logarithmic one.

A central configuration is a configuration of bodies such that the acceleration vector for each body is a common scalar multiple of its position vector:

∂ 2 ∂t 2 r i = λ r i where λ is real 0 (2)
Central configurations are of interest for a variety of reasons: To every central configuration is corresponding a homothetic solution -a solution which retains its shape for all time, while expanding, contracting and rotating about the center of mass.

We identify the real plane R 2 with the complex plane C. As we will see, in the planar case, it is easier to work with complex positions z i = x i + ı y i = r i . Hence s i, j = |r ir j | 2 = (z iz j ) (z iz j ) where z is the complex conjugate of z. Combining (1) and ( 2) we obtain:

λ z i = ∑ j =i (z i -z j ) (z i -z j ) (z i -z j ) = ∑ j =i 1 zi -z j (3)
By summing over the i, we see that ∑ z i = 0 (the center of mass is at the origin), so e 1 , the first symmetric function of the z i , is equal to zero. Moreover, since λ > 0, we observe that the conjugate of equation ( 3) is equivalent to:

√ λ zi = N ∑ j=1, j =i 1 √ λ z i - √ λ z j
Therefore, we can suppose λ = 1 and recover the original solutions by multiplying the solutions of the case λ = 1 by √ λ . Knowing z i , it is easy to recover λ using the following property

Proposition 1 2λ ∑ N i=1 |z i | 2 = N (N -1) PROOF. λ ∑ N i=1 z i zi = ∑ N i=1 ∑ N j=1, j =i z i z i -z j = ∑ N i=1 ∑ N j=1, j =i (1 + z j z i -z j ) and thus λ ∑ N i=1 z i zi = N(N -1) -λ ∑ N j=1 z j z j
In summary, the central configuration problem is equivalent to:

zi = N ∑ j=1, j =i 1 z i -z j (E i )
2.2 Symmetries of the solutions

Action of S N

The permutation group S N acts on the variables {z 1 , . . . , z N } with σ (z i ) = z σ (i) , and σ sends the equation (E i ) on the equation E σ (i) . So, if (z 1 , . . . , z N ) is a solution of the problem, any of the N! N-uples (z σ (1) , . . . , z σ (N) ) is still a solution.

Action of O 2 (R)

The isometry group of R 2 can be identified to a transformation group on C generated by the rotations z → az with a a complex of modulus one, and the symmetry z → z. These transformations act on C N by (z 1 , . . . , z N ) → (az 1 , . . . , az N ) and (z 1 , . . . , z N ) → ( z1 , . . . , z N ). If we multiply (E i ) by ā or if we conjugate (E i ), we see that the set of solutions is invariant under this action.

Consequently, the set of solutions is invariant under the action of S N × O 2 (R). We will first focus on the action of S N to obtain invariant equations, and finally use the action of O 2 (R) to speed up the Gröbner Basis computation (see section 5).

Algebraic reformulation

Algebraically, it is impossible to separate ı and -ı, and therefore z and z. Thus, we introduce N new variables Z 1 , . . . , Z N that represent z1 , . . . , zN . The algebraic relations between these 2 N variables are :

Z i = ∑ j =i 1 z i -z j and z i = ∑ j =i 1 Z i -Z j (E i , Ēi )
To obtain polynomials, we multiply the equation E i by D i = ∏ j =i (z iz j ) to obtain the polynomial equation U i = 0 where

U i = Z i ∏ j =i (z i -z j ) -∑ j =i ∏ k =i, j (z i -z k ) ∈ Q[z 1 , . . . , z N , Z 1 , . . . , Z N ]
Remark 1 Observe that permuting z i and Z i for all i transforms the equation (E i ) in ( Ēi ), because of complex conjugation. Thus, for every relation in the ideal generated by the 2N equations

(E i , Ēi ) in Q[z 1 , . . . , z N , Z 1 , . . . , Z n ],
there is another one obtained by permuting z i and Z i .

The following lemma is useful to express the equations (E i , Ēi ) in a very compact way.

Lemma 1 All the solutions of the vortex problem satisfy the following rational equations:

Z i = 1 2 Q ′′ (z i ) Q ′ (z i ) (4)
where

Q(z) = ∏ N i=1 (z -z i ) = z N + e 2 z N-2 + • • • + (-1
) N e N (recall that e 1 = 0), where e i is the i-th elementary symmetric function of the z i .

PROOF. Let Q i (z)= Q(z) z-z i = ∏ j =i (z-z j ), then Q ′ i (z) Q i (z) = ∑ j =i 1 z-z j , so that Q ′ i (z i ) Q i (z i ) = Z i according to equation (E i ). But we can write Q(z) = (z -z i )Q i (z)
, and with two derivations, we obtain

Q ′ (z) = Q i (z) + (z -z i )Q ′ i (z) and Q ′′ (z) = 2Q ′ i (z) + (z -z i )Q ′′ i (z). Setting z = z i , we have Q i (z i ) = Q ′ (z i ) and Q ′′ (z i ) = 2 Q ′ i (z i
) so that we get the proof of the lemma.

S N acts on {1, . . . , N}, and therefore on {z 1 , . . . , z N , Z 1 , . . . , Z N }. The next step is to obtain equations depending only of the e i , the symmetric functions of the z i . In the next section we will see how to obtain equations of lower degree individually invariant under S N .

FROM INVARIANT SYSTEM TO INVARIANT EQUATIONS

In this section, we will generalize the previous situation with more than two blocks of variables. More precisely we assume that we have to solve the following polynomial system:

U i = 0 for i = 1, . . . , N
where each equation U i is a polynomial in A[Z , V ] where A is an integral domain (for instance, a polynomial ring with coefficients in K), Z is the set {z 1 , . . . , z N } and V is another set of variables. We assume that S N acts on Z ∪ V .

Invariant system under S N

We suppose that |V | is a multiple of N, and that S N acts on each block of N variables like it acts on Z such that σ

(z i ) = z σ (i) . In fact, the product S N × • • • × S N acts on the variables Z ∪ V . Let {(U i )} be a system of N equations, which is globally invariant under the subgroup of S N × • • • × S N of which elements are of the form σ × • • • × σ : that means that for every σ = σ × • • • × σ and every i, it exists j such that σ .U i = U j .
Remark 2 We could introduce the ring of invariants of this subgroup, but its Hironaka Decomposition [START_REF] Sturmfels | Algorithms in invariant theory[END_REF] is not simple enough to allow easy computations.

We want to obtain from the original set of equations {U i } a new set of equations {V i } which are individually invariant under the action of S N , which means that for all i and σ , σ .V i = V i . To this end, we will use divided differences.

Divided differences on one block

Here, we first assume that we only have one block of variables

Z = {z 1 , . . . , z N }, V = / 0 and N equations U i ∈ A[z 1 , . . . , z N ] such that σ (U i ) = U σ (i) for all σ in S N .
Definition 1 Given U 1 , . . . ,U N we define recursively the divided differences by:

[U i ] = U i for i = 1, . . . , N [U i 1 , • • • ,U i k ] = [Ui 1 ,••• ,U i k-1 ]-[Ui 1 ,••• ,U i k-2 ,U i k ] z i k-1 -z i k
for any given distinct integers {i 1 , . . . , i k } in {1, . . . , N}. PROOF. We prove the first part by induction on k ∈ {1 . . . , N}. For k = 1 it is obvious. Suppose the theorem is true for k -1, and let {i 1 , . . . , i k } be like in the statement.

Theorem 1 The divided difference [U i 1 , • • • ,U i k ] is
Let à = A[z 1 , . . . , z i k-1 -1 , z i k-1 +1 , . . . , z N ]. Since z i k-1 -z i k is monic as polynomial in Ã[z i k-1 ], we can perform the division of U i 1 , • • • ,U i k-1 -U i 1 , • • • ,U i k-2 ,U i k by z i k-1 -z i k ,

and by sending

z i k-1 on z i k , we see that the rest is equal to 0, so U i 1 , • • • ,U i k-2 ,U i k belongs to A[Z ]
. Moreover, by acting on the equality The two following lemmas will be useful later.

[U i 1 , • • • ,U i k ] (z i k-1 -z i k ) = U i 1 , • • • ,U i k-1 -U i 1 , • • • ,U i k-2 ,U
Lemma 2 For an univariate polynomial F(z) ∈ K[z] we have F(z 1 , . . . , z N ) = ∑ N i=1 F(z i ) Q ′ (z i ) where Q(z) = ∏ N i=1 (z -z i ). PROOF.
We prove this by induction on N. For N = 1, Q(z) = (zz 1 ) so the assertion is obvious. Suppose that we have proved this lemma for

N -1. Let U(z) = ∏ N-1 i=1 (z -z i ) and V (z) = ∏ N-2 i=1 (z -z i ) × (z -z N ). Hence, Q(z) = U(z)(z -z N ) = V (z)(z -z N-1 ) which implies that Q ′ (z) = U(z)+U ′ (z)(z-z N ) = V (z)+V ′ (z)(z- z N-1 ). Consequently, F(z 1 , . . . , z N ) = F(z 1 ) U ′ (z 1 ) + F(z N ) V ′ (z N ) +∑ N-1 i=2 F(z i ) U ′ (z i ) + F(z i ) V ′ (z i ) z N-1 -z N = ∑ N k=1 F(z k ) Q ′ (z k )
, and the lemma is proved. Definition 3 For Z = {z 1 , . . . , z N }, we define h k the k-th symmetric complete function as the sum of all monomials of degree k on the variables in Z . By extension, h k = 0 when k < 0 and h 0 = 1.

Lemma 3 For any k, if F(z) = z k , then F(z 1 , . . . , z N ) = h k-N+1 .
PROOF. We prove this by induction on N as well.

If N = 1, F(z 1 ) = z k 1 is the complete symmetric function of degree k on one variable z 1 . For N ≥ 2, F(z 1 , . . . , z N ) = F(z 1 ,...,z N-1 )-F(z 1 ,...,z N-2 ,z N ) z N-1 -z N By induction, F(z 1 , . . . , z N-1 )-F(z 1 , . . . , z N-2 , z N ) = ∑(z k-N+2-u N-1 -z k-N+2-u N
) × m, where the sum is all over the monomials m in z

1 , . . . , z N-2 of degree u ∈ {0, . . . , k -N + 2}. Writing z k-N+2-u N-1 - z k-N+2-u N = (z N-1 -z N ) ∑ m ′ ,
where the sum is all over the monomials m ′ in z N-1 , z N of degree k -N + 1u, we obtain exactly the complete symmetric function in z 1 , . . . , z N of degree k -N + 1.

Invariant equations

We explain here how to obtain invariant equations from divided differences in case of one block of variables.

Theorem 2 Let V i be ∑ P⊂{1,...,N},|P|=i

[U] P for all i ∈ {1, . . . , N}.

We obtain N equations invariant under S N , and the varieties associated respectively to {V i } and {U i } are the same, except maybe for points with two equal components.

PROOF. Any σ in S N realizes a permutation of the subsets of {1, . . . , N} with same cardinality, and also a permutation of the [U] P . Therefore, σ (V i ) = V i for all i in {1, . . . , N}. Suppose that α = (α 1 , . . . , α N ) is a common zero of the U i , without equal components. Then, we deduce easily that all the [U] P (α) are equal to zero, and also the V i (α). Conversely, if V N (α) = 0 then all the [U] P (α) with P of cardinality N -1 are equal, because V N can be written as

[U] P -[U] Q z k -z ℓ
where P and Q are two distinct subsets of cardinality N -1, z k = P\Q and z ℓ = Q\P. But their sum V N-1 (α) is equal to zero, so they are equal to zero. We can repeat it for N -2, N -3, . . . , 1 to deduce that U i (α) = 0 for all i.

Definition 4 We define the Reynolds operator [START_REF] Cox | Ideals, varieties, and algorithms[END_REF] on S N by

R : A[z 1 , . . . , z N ] -→ A[z 1 , . . . , z N ] P → 1 N! ∑ σ σ (P)
Remark 3 We just have to compute all the [U 1 , . . . ,U k ] for k in {1, . . . , N} to obtain the V i , because V i = i N R([U 1 , . . . ,U i ]). We deduce from this fact a simple algorithm to compute the set {V i }.

Let τ i, j be the transposition permuting i and j.

ComputeInvariantSystem Input: Variables {z 1 , . . . , z N } and the polynomials U i .

Output: The polynomials V 1 , . . . ,V N . 1. for k = 2 to N do [U 1 , . . . ,U k ] :=Quo([U 1 , . . . ,U k-1 ] -τ k-1,k ([U 1 , . . . ,U k-1 ]), z k-1 -z k ); end for. 2. return{ k N R([U 1 , . . . ,U k ]), k = 1 . . . N} where R is the Reynolds op- erator. Since K[z 1 , . . . , z N ] S N = K[e 1 , . . . , e N ],
we can reformulate the equations V i in terms of the symmetric functions of the z i . The goal is now to generalize this fact in the case of several blocks of variables.

Generalization to several blocks

Assume that the equations U i involve the set Z = {z 1 , . . . , z N } and another set of variables V and that S N acts on Z ∪V . Assume that the U i are equal to D i P i + Q i , where D i = ∏ j =i (z iz j ), Q i are polynomials in Z , and for all σ , σ (P i ) = P σ (i) and σ (Q i ) = Q σ (i) . The previous section corresponds to the case P i = 0, but when P i = 0 we can still apply divided differences in the same way, and construct [U i 1 , . . . ,U i k ] for given distinct integers; we obtain a similar theorem :

Theorem 3 (i) [U i 1 , . . . ,U i k ]
is a polynomial in Z and V which depends only on the set {i 1 , . . . , i k }. Moreover for any σ and any P, σ ([U] P ) = [U] σ (P) .

(ii) V i = ∑ |P|=i [U] P is invariant under the action of S N and the varieties associated to respectively V i and U i are the same, except maybe for points with two equal Z -components.

Remark 4

We can use the algorithm ComputeInvariantSystem to compute the V i .

Application to the vortex problem

We obtained this equation from the vortex problem :

U i = Z i ∏ j =i (z i -z j ) -∑ j =i ∏ k =i, j (z i -z k )
which can be written as

U i = D i P i + Q i , with P i = Z i and Q i = -∑ j =i ∏ k =i, j (z i -z k ). These polynomials verify σ (P i ) = P σ (i) and σ (Q i ) = Q σ (i) .
Example 1 For N = 3, it is easy to compute the invariant polynomials V 1 ,V 2 ,V 3 , and we obtain

V 1 = ∑ i Z i z 2 i -∑ j =i Z i z i z j + ∑ {i, j,k}={1,2,3} Z i z j z k V 2 = 2 ∑ i z i z i -∑ i = j Z i z j -9 V 3 = ∑ i Z i
Using the equations V 3 = 0 and V3 = z 1 + z 2 + z 3 = 0 in V 1 and V 2 , we can rewrite the system as

V 1 = 4 ∑ i Z i z 2 i V 2 = 3 ∑ i Z i z i -9 V 3 = 0 Definition 5
We introduce new invariants which will be very useful to reformulate the vortex problem in a very simple way:

r k = ∑ i Z i z k i and R k = ∑ i Z k i z i
It is also useful to define the Newton sums:

s k = ∑ i z k i and S k = ∑ i Z k i
Example 2 For N = 4, after reformulation, we obtain:

r 0 = s 1 = r 1 -6 = r 2 = 2 r 3 + 5 s 2 = 0
and the conjugate equations:

R 0 = S 1 = R 1 -6 = R 2 = 2 R 3 + 5 S 2 = 0
It turns out that, surprisingly, we can obtain a general and very simple expression of these equations for any N.

Theorem 4 (Invariant Equations) For any N and k, the solution of the vortex problem satisfies the following invariant equations:

2 r k = k-1 ∑ i=0 s i s k-1-i -k s k-1 with s 0 = N. (5)
To prove the theorem 4 we will first give a lemma.

Lemma 4 For N ∈ N ⋆ and k ∈ N, let e k , h k and s k be respectively the symmetric functions, the complete symmetric functions, and the Newton sum of degree k in the variables z 1 , . . . , z N . (With the convention that e 0 = h 0 = 1, s 0 = N and all these terms are zero when k < 0 and e k is zero for k > N). Then

N ∑ j=0 (-1) j e j (N -j)(N -j -1) h k-j-1 = k-1 ∑ i=0 s i s k-1-i -k s k-1
PROOF. We introduce the series

E(t) = ∑ e k t k = ∏ N i=1 (1 + t z i ), H(t) = ∑ h k t k = ∏ N i=1 1 1-t z i and S(t) = ∑ s k t k = ∑ N i=1 1 1-t z i . Let β k-1
and γ k-1 be the left and right terms of the equality we want to prove. We introduce β (t) and γ(t) the associated power series. Observe that (Nj)(Nj -1) = N 2 -N -2(N -1) j + j( j -1), so β (t) = (N 2 -N)E(-t)H(t) + 2(N -1)tE ′ (-t)H(t) + t 2 E ′′ (-t)H(t) and γ(t) = S 2 (t) -S(t) -tS ′ (t). But H(t)E(-t) = 1, and with two derivations, we have -E ′ (-t)H(t) + E(-t)H ′ (t) = 0 and E ′′ (-t)H(t) -2E ′ (-t)H ′ (t) + E(-t)H ′′ (t) = 0, so we can express β with H and its derivatives :

β (t) = (N 2 -N) + 2(N - 1) tH ′ (t) H(t) + 2 t 2 H ′2 (t) H 2 (t) - t 2 H ′′ (t)
H(t) . Since t H ′ (t) = H(t)(-N + S(t)), γ can be expressed with H and its derivatives too, and we find the same expression, that proves that β (t) = γ(t) and therefore

β k-1 = γ k-1 .
PROOF. We give now the proof of theorem 4. By lemma 1,

r k = ∑ N i=1 z k i Z i = ∑ z k i Q ′′ (z i ) 2Q ′ (z i ) = ∑ F(z i ) Q ′ (z i ) where F(z) = z k Q ′′ (z) 2
.

Writing

Q(z) = z N -e 1 z N-1 + e 2 z N-2 + • • • + (-1) N e N , we obtain F(z) = 1 2 ∑ N j=0 (-1) j (N -j)(N -j -1) e j z N-j+k-2 .
Using linearity and lemma 2, we have exactly 2r k = ∑ N j=0 (-1) j e j (Nj)(Nj -1) h k-j-1 . Using lemma 4, we obtain the theorem 4.

FROM TWO BLOCKS TO SYMMETRIC FUNCTIONS IN ONE BLOCK

General case under the rational parameterization assumption

We return now to the general case where each U i is an equation in A[Z , Z ] where Z = {z 1 , . . . , z N } and Z = {Z 1 , . . . , Z N }; in addition, we require that the algebraic system fulfils the following parameterization assumption: Definition 6 We say that the system U i = 0 is under parameterization assumption if for all i, Z i = R(z i ) where R(z) = N(z) M(z) ∈ A(z) with A = Q(e 1 , e 2 , . . . , e N ) an univariate rational function whose coefficients depend on the symmetric functions of the z i .

This hypothesis is well adapted to algebraic problems in the plane and more specifically to the vortex problem since lemma 1 can be reformulated: Proposition 1 The vortex problem satisfies the rational assumption since

Z i = 1 2 Q ′′ (z i ) Q ′ (z i ) with Q(z) = ∏ N i=1 (z -z i ) = z N + e 2 z N-2 + • • • + (-1) N e N .
For now, we assume that the system is under this assumption. We now describe an algorithm to obtain invariant equations under the action of S n . First, we apply the algorithm ComputeInvariantSystem to compute the invariant equations V i . Denote again by

Q(z) the polynomial ∏ i (z -z i ) = z N -e 1 z N-1 + • • • + (-1) N e N . There exist two polynomials B and C in K[e 1 , . . . , e N ][z] such that B Q + C M = R M , where R M is the resultant of Q and M with respect to the variable z. Then, R M Z i = R M N(z i ) M(z i ) = N(z i )C(z i ) since Q(z i ) = 0. More generally, R k M Z k i = (N k ×C k )(z i ) = (N k ×C k mod Q)(z i ) for all k ≥ 0. For each W ∈ {V i , Vi }, we substitute 1 R k M (N k ×C k mod Q)(z i ) to Z k i in each monomial of W .
Up to a multiplication by an appropriate power of R k M to obtain polynomials, we obtain equations involving only the variables z 1 , . . . , z N . These polynomials are invariant under S N , and can be expressed as equations in the e i .

These ideas give the following algorithm, we denote by ∂ Z P the total degree of P as polynomial of A[Z 1 , . . . , Z N ] with A = K[z 1 , . . . , z N ]. For all polynomial P in K[z 1 , . . . , z N ] S N , let Σ(P) be the expression of P as polynomial of K[e 1 , . . . , e N ].

ComputeSymmetricFunctionsSystem Input:

The invariant equations V i , Vi of variables Z = {z 1 , . . . , z N } and Z = {Z 1 , . . . , Z N } Output: A system of 2N equations of variables e i , the symmetric functions of the z i .

1. m := Max{∂ Z W, ,W ∈ {V i , Vi }}; L NC := [(NC) i mod Q, i ∈ 1..m]; 2. for W in {V i , Vi , i = 1..N} do d W := ∂ Z (W ); for U monomial of W do d U := ∂ Z (U); Z k i ← L NC [k](z i ) in U. end for. U ← R d W -d U M U in W . end for. 3. return {Σ(W )}
Example 3 For the vortex problem, according to (4), the denominator M is equal to Q ′ , and the relation that we have is BQ +CQ ′ = D, where D is the discriminant of Q(z) with respect to the variable z. Let P Z (z) be the polynomial of K[e 2 , . . . , e N ][z] equal to 1 2 Q ′′ C mod Q. Then, P Z sends z i on D Z i , and we can apply the previous algorithm to compute symmetric equations. From V k = 2r k -∑ k-1 i=0 s i s k-1-i + ks k-1 = 0, we obtain always 0, but not from Vk . In this specific case, instead of using previous algorithm, there is a faster way to compute the equations, explained in the subsection hereafter.

Invariant system of equations for the vortices problem

We introduce the two K[e 2 . . . , e N ]-modules morphisms :

S : K[e 2 , . . . , e N ][z] -→ K[e 2 , . . . , e N ] P(z) = ∑ a k z k -→ S(P) = ∑ a k s k H : K[e 2 , . . . , e N ][z] -→ K[e 2 , . . . , e N ] P(z) = ∑ a k z k -→ ∑ a k h k-N+1
Proposition 2 For any polynomial P in K[e 2 , . . . , e N ][z] we have S (P) = S (P mod Q) and H (P) = H (P mod Q). Moreover

S (P) = ∑ N i=1 P(z i ) and H (P) = ∑ N i=1 P(z i ) Q ′ (z i ) If deg(P) < N then H (P) = a N-1 .
From theorem 4, we know that 2

r k = ∑ k-1 i=0 s i s k-1-i -k s k-1 , therefore we have the conjugate equation: 2 R k = k-1 ∑ i=0 S i S k-1-i -k S k-1 (6)
One way to obtain directly symmetric equations is to compute: [START_REF] Faugère | The membrane inclusions curvature equations[END_REF] we obtain:

S k = ∑ N i=1 Z k i = 1 D k ∑ N i=1 P k Z (z i ) = 1 D k S (P k Z (z) mod Q) R k = ∑ N i=1 z i Z k i = 1 D k ∑ N i=1 z i P k Z (z i ) = 1 D k S (z P k Z (z) mod Q) Substituting these expressions in
Proposition 3 Given the Bézout relation B(z) Q(z) +C(z) Q ′ (z) = D,
for any N and k, the solution of the vortex problem satisfies the following symmetric equations:

2 1 D S (z P k Z ) = k-1 ∑ i=0 S (P i Z ) S (P k-1-i Z ) -k S (P k-1 Z )
where S (1) = N and

P Z = 1 2 Q ′′ C.
The drawback of the method is that high powers of the discriminant occur in the resultant equations. Instead of using the polynomial P Z to obtain equations with the e i 's, as explained in the previous section we will give a more efficient method using the H morphism.

Using the following lemma it is possible to compute R k and S k with powers of the discriminant divided by two:

Lemma 5 Given the Bézout relation B(z) Q(z) + C(z) Q ′ (z) = D, we have D Q ′′ Q ′2 (z k ) = A(z k )
where A(z) is the polynomial -(B(z) +C ′ (z)).

PROOF. By derivating the relation

B(z) Q(z) + C(z) Q ′ (z) = D, we obtain B ′ (z)Q(z) + (B ′ (z) +C(z))Q ′ (z) +C(z)Q ′′ (z) = 0. Mod- ulo Q(z), we have -(B + C ′ )Q ′ = CQ ′′ , so that -(B + C ′ )Q ′ C = -(B +C ′ )D = CQ ′′ where A(z) = -(B(z) +C ′ (z)). Modulo Q, we deduce that A = Q ′′ C 2 D = DQ ′′ Q ′2 .
Hence, with one power of A, there are two powers of Q ′ in the denominator, and only one power of D in the numerator. If we use the morphism H when k is odd, we obtain R k and S k . More exactly:

Homogeneity of the equations

= N(N - 1)⌊ k 2 ⌋ + 1 -k. PROOF. We started from 2R k = ∑ S i S k-1-i -kS k-1 . With Z i = ∑ k =i 1 z i -z k , we see that this equation is homogeneous in the z i of degree 1 -k. The discriminant D = ∏ i = j (z i -z j ) is homogeneous of degree 2 N 2 .
So, the previous equation is homogeneous in the z i with degree 2⌊ k 2 ⌋ N 2 + 1k. The symmetric function e k is homogeneous in the z i of degree k, that's why we took the degree d.

Recall that we have lost the equation 2r 1 = 2R 1 = N(N -1), but there is no surprise : we have seen that the set of solutions (the z i 's) is invariant under multiplication by a complex of modulus one. This implies that the algebraic variety with variables (e 2 , . . . , e N ) is invariant under the operation (e 2 , . . . , e N ) → (γ 2 e 2 , . . . , γ N e N ), with |γ| = 1. But an ideal associated to such a variety is homogeneous for the previous degree d : let P a polynomial in this ideal, and write P = ∑ u P u , with P u the homogenous part of degree u for the degree d, then if (e 2 , . . . , e N ) is a zero of P; we have ∑ u γ u P u (e 2 , . . . , e N ) = 0 for all γ of modulus 1. A non-zero univariate polynomial have only a finite number of roots, so this polynomial (in γ) is null and P u (e 2 , . . . , e N ) = 0 for all u. Another argument is that for k = 1, the degree of R 1 (in z i ) is zero, so we have to obtained from 2R 1 = N(N -1) an homogeneous equation of degree d equal to 0. There are solutions for all N to the vortex problem (see [START_REF] Aref | Vortex crystals[END_REF]), so this equation is 0.

Strategy to obtain a Gröbner Basis

Because of the homogeneity, we can suppose that any of the symmetric functions e i is equal to 1 or 0. If it is 0, we have again an homogeneous system, so we can suppose that another symmetric function e j is equal to 1 or 0, and so on. We have to add a new equation for being sure that all the z i are distinct : h × D = 1. (We can't solve it and remove the spurious solutions easily : for example, for N = 5, the system with e 2 = 1 without h × D = 1 is 1-dimensional.)

According to the benchmark, it seems that the fastest way to compute a Gröbner Basis is to separate the system into two parts, e 2 = 1 or e 2 = 0 and compute a Gröbner Basis with DRL order with h > e N > • • • > e 3 , and then perform a change of ordering from DRL to the lexicographic order with FGLM (see [START_REF] Faugère | Efficient Computation of Zero-dimensional Gröbner Bases by Change of Ordering[END_REF][START_REF] Faugère | Fast Algorithm for Change of Ordering of Zero-dimensional Gröbner Bases with Sparse Multiplication Matrices[END_REF]). For the component with e 2 = 0, we separate e 3 = 1 or e 3 = 0, and so on.

Then, we perform a Triangular Decomposition (see [START_REF] Lazard | A new method for solving algebraic systems of positive dimension[END_REF]) of each component.

Remark 6

To compute a Gröbner Basis, we supposed that e k = 1 for some k. But with this assumption, the solutions (z 1 , . . . , z N ) that we obtain are not solutions of the equations

(E i ) zi = ∑ j =i 1 z i -z j but of λ zi = ∑ j =i 1 z i -z j
for some λ > 0. Denote by (az 1 , . . . , az N ) the solutions of (E i ), where a can be supposed to be a positive real.

Then

2r 1 = 2R 1 = 2a 2 ∑ i |z i | 2 = N(N -1), and a = N(N-1) 2 ∑ |z i | 2 . The true value of e k is ∑ a z i 1 × • • • × a z i k = a k .
Example 5 With e 2 , e 3 or e 4 equal to 1, the solutions (z 1 , . . . , z 4 ) for N = 4 are drawn below :

In the case of the four aligned points, ∑ i |z i | 2 is equal to 2, so we have to perform a multiplication by √ 3 to obtain the solutions of

Removing spurious solutions

We can solve the system to obtain approximations of the e i and then approximations of the z i , but there are spurious solutions: we have to check that P Z (z i ) = D zi for each i to be sure that we have computed a true solution. Another way to perform it is to introduce two news variables x and z and add to the system the equations P Z (z) + D z = x and z N + e 2 z N-2 + • • • + (-1) N e N , with P Z the polynomial computed previously, which maps z i to DZ i . The next step is to perform a Gröbner elimination with lexicographical order z > e N > • • • > e k > x to obtain an univariate polynomial P ℜ in x. Then we isolate the real roots of this polynomial P ℜ using certificated methods.

Other symmetries

Suppose that we are in the component e 2 = 1, we have said that if e 2 , . . . , e N is a solution, then λ 2 e 2 , . . . , λ N e N too, for all λ of modulus 1. If λ = -1 (geometrically, we do a symmetry of center O), e 2 stays at 1, but e 3 is changed to -e 3 , so we can keep only half of the possible e 3 . (e 2 , . . . , e N ) → ( ē2 , . . . , ē N ) gives an other solution, so if e 3 is not real, we can suppose e 3 of imaginary part non negative. If e 3 is real and e 4 not, we can keep only the e 4 with imaginary part non negative, and so on.

Naive Approach

It is possible to solve directly the original system of 2N equations (E i , Ēi ) in z i and Z i . Because of invariance by multiplication by a scalar of modulus 1, we can suppose that z 1 is real, so we add the equation z 1 = Z 1 . This trick will give an ideal of dimension 0, if we suppose that z 1 = 0. We will split the ideal into two parts : in the first one, we add the equation z 1 × α = 1, and in the second one, we add z 1 = 0, and we can add z 2 = Z 2 . In each case, the ideal is zero dimensional, if we add the last equation ∏ i< j (z iz j )β = 1, for being sure that all the z i are distinct. We report in table 1 the following timings in Magma to compute the corresponding Gröbner basis (∞ means that we stopped the computation after five days):

It is possible to introduce the invariant ring of the subgroup of

Generating and solving the symmetric system

We have implemented the algorithm ComputeSymmetricFunc-tionsVorticesSystem in Maple and Magma to generate the symmetric system. We report in table 3 the timings to compute the systems depending only on the symmetric functions e i using Com-puteSymmetricFunctionsVorticesSystem algorithm with Magma (Intel Xeon 2.93 GHz with 128GB Ram). On the same computer, the times to compute a Gröbner Basis using Magma of the symmetric system and perform a triangular decomposition of each component (mostly for the component with e 2 = 1) are presented in table 4.

When N = 7 we use FGb [START_REF] Faugère | FGb: A Library for Computing Gröbner Bases[END_REF] to compute the corresponding Gröbner bases: it takes 144 sec to compute the system over F 65521 and about 20 minutes to compute a Gröbner basis and a triangular decomposition over Q. We postpone to section 7 the complete prime decomposition of the ideal corresponding to N = 7; using all the symmetries the problem admits 12 solutions; among them, 2 solutions are expressed as algebraic numbers of degree 82. For N = 8 the computation is still running but the most difficult part is already done (it takes 12 days to compute the first Gröbner basis). A web page was created to collect all the data: http://wwwsalsa.lip6.fr/~jcf/vortices/

  i k with any permutation σ , we deduce the second part again by induction on N. It is usual to introduce a special notation in the case of univariate polynomials: Definition 2 Let K be a field, and F(z) an univariate polynomial in K[z]. We define F(z 1 , . . . , z N ) =[F(z 1 ), . . . , F(z N )].

Proposition 2

 2 The equation we obtained in the previous section are homogeneous for the degree d = ∑ k k × ∂ e k , where ∂ e k is the degree in e k . More exactly, the k-th equation has degree d

Figure 1 :Figure 2 :

 12 Figure 1: Q(z) = z 4 + z 2 + 1 12

Figure 3 :

 3 Figure 3: Q(z) = z 4 + 1

  Magma 0.0s 0.0s 0.06s 70.6s 7649.6s Maple 0.0s 0.2s 0.9s 41.9s 2407.3s

  a polynomial in Z and depends only on the set {i 1 , . . . , i k }, so for any subset P = {i 1 , . . . , i k }, we set [U] P = [U i 1 , . . . ,U i k ]. Moreover, for any subset P of {1, . . . , N}, and for any σ in S N , σ ([U] P ) = [U] σ (P) .

Table 1 :

 1 Direct approach: Gröbner bases of the non symmetric systems with Magma.S N × S N , the elements of which are (σ , σ ). We report in table 2 the number of secondary invariants in the Hironaka decomposition or the number of fundamental invariants over Q, and the timings to compute them in Magma.

		3		4	5
	Q	0.02s 176.8s ∞
	F 65521 0.01s 0.2s	∞
			3	4	5	6	7
	Secondary Invariants		6	24	120	?	?
	Timings	0.0s 0.1s 225s	∞	∞
	Fundamental Invariants	9	14	20	27	?
	Timings	0.0s 0.1s 3.0s 400s ∞

Table 2 :

 2 Invariant Ring : Hironaka Decomposition and Fundamental Invariants with Magma.

Table 3 :

 3 Time to generate the symmetric systems with Maple or Magma.

Proposition 4

The expressions of S i and R i in terms of the symmetric functions of the z i 's are :

and all polynomials could be taken modulo Q.

Substituting these expressions in [START_REF] Faugère | The membrane inclusions curvature equations[END_REF] we obtain:

Theorem 5 (Symmetric Equations) Given the Bézout relation B(z) Q(z) + C(z) Q ′ (z) = D, for any N and k, the solution of the vortex problem satisfies the following symmetric equations:

where

PROOF. We substitute the expression of R 2k , R 2k+1 , S 2k , S 2k+1 given by proposition 4 into the equations:

This theorem gives a very efficient algorithm to compute a system involving only the e i , which solutions include all symmetric functions of the vortex problem. To simplify the description of the algorithm we introduce the following notation α i,k and β k , which depend only on the parity of i and k:

is even and k odd 1 otherwise

ComputeSymmetricFunctionsVorticesSystem Input: N, the polynomials Q, D and A = -B -C ′ , where B and C appear in the Bézout relation B Q +C Q ′ = D, and the two functions S and H Output: Symmetric polynomials in the e i 's.

Remark 5 The equation 2R 1 = N(N -1) gives always 0 = 0. We explain this fact in the next section. ). The next step is to replace z 3 by s 3 = 3e 3 , z 2 by s 2 = -2e 2 and z by p 1 = 0 so that 0 = S (z AQ ′′ ) = -16 e 3 (12 e 4 + e 2

2 ) 2 . In the same way, we compute the second equation

SOLVING THE EQUATIONS WITH THE SYMMETRIC FUNCTIONS

The goal of this section is to solve explicitly the symmetric equations by exact methods. To achieve this we use Gröbner bases computation.

The case N = 4

Interestingly enough, we can solve the vortex problem by hand when N = 4. Hence, we give the complete resolution of the case N = 4 without Gröbner Basis computation: The symmetric equations are given by [START_REF] Faugère | Fast Algorithm for Change of Ordering of Zero-dimensional Gröbner Bases with Sparse Multiplication Matrices[END_REF] and, in addition, we assume that the discriminant D = 16e 4 2 e 4 -4e 3 2 e 2 3 -128e 2 2 e 2 4 + 144e 2 e 2 3 e 4 -27e 4 3 + 256e 3 4 = 0, to ensures that the z i 's are distinct. Lemma 6 In equations [START_REF] Faugère | Fast Algorithm for Change of Ordering of Zero-dimensional Gröbner Bases with Sparse Multiplication Matrices[END_REF], if e 2 = 0, then e 3 = 0.

PROOF. We prove it by reduction to the absurd. If e 2 = 0 and e 3 = 0, the first equation gives e 4 = - 1 12 e 2 2 , and the second equation becomes 8e 3 2 + 27e 2 3 = 0, but (8e 3 2 + 27e 2 3 ) 2 is the discriminant of Q, up to a constant factor, which is a contradiction.

Then, if e 2 = 0, e 3 = 0, and the second equation becomes (e 2 2 -12e 4 )(e 2 2 -4e 4 ) = 0, but D becomes 16e 4 (e 2 2 -4e 4 ) 2 = 0, so e 4 = 1 12 e 2 2 . If e 2 = 0 then e 3 = 0 or e 4 = 0. We can conclude that: Proposition 5 When N = 4, there are three solutions to the vortex problem :

The indetermination on e 2 , e 3 or e 4 will be explained and solved in the next section as shown in the figures 1,2 and 3.

APPENDIX: SOLUTIONS FOR N=7