Improvement of FPPR method to solve ECDLP

Solving the elliptic curve discrete logarithm problem (ECDLP) by using Gröbner basis has recently appeared as a new threat to the security of elliptic curve cryptography and pairing-based cryptosystems. At Eurocrypt 2012, Faugère, Perret, Petit and Renault proposed a new method (FPPR method) using a multivariable polynomial system to solve ECDLP over finite fields of characteristic 2. At Asiacrypt 2012, Petit and Quisquater showed that this method may beat generic algorithms for extension degrees larger than about 2000. In this paper, we propose a variant of FPPR method that practically reduces the computation time and memory required. Our variant is based on the idea of symmetrization. This idea already provided practical improvements in several previous works for composite-degree extension fields, but its application to prime-degree extension fields has been more challenging. To exploit symmetries in an efficient way in that case, we specialize the definition of factor basis used in FPPR method to replace the original polynomial system by a new and simpler one. We provide theoretical and experimental evidence that our method is faster and requires less memory than FPPR method when the extension degree is large enough.


Introduction
In the last two decades, elliptic curves have become increasingly important.In 2009, the American National Security Agency (NSA) to advocate the use of elliptic curves for public key cryptography [14] which are based on the hardness of elliptic curve discrete logarithm problem (ECDLP) or other hardness problem on elliptic curves.Elliptic curves used in practice are defined either over a prime field F p or over a binary field F 2 n .Like any other discrete logarithm problem, ECDLP can be solved with generic algorithms such as Baby-step Giant-step algorithm, Pollard's ρ method and their variants [1,16,17,19].These algorithms can be parallelized very efficiently, but the parallel versions still have an exponential complexity in the size of the parameters.Better algorithms based on the index calculus framework have long been known for discrete logarithm problems over multiplicative groups of finite fields or hyperelliptic curves, but generic algorithms have remained the best algorithms for solving ECDLP until recently.
A key step of an index calculus algorithm for solving ECDLP is to solve the point decomposition problem.
In 2004, Semaev introduced the summation polynomials (also known as Semaev's polynomials) to solve this problem.Solving Semaev's polynomials is not a trivial task in general, in particular if K is a prime field.For extension fields K = F q n , Gaudry and Diem [2,9] independently proposed to define V as the subfield F q and to apply a Weil descent to further reduce the resolution of Semaev's polynomials to the resolution of a polynomial system of equations over F q .Diem generalized these ideas by defining V as a vector subspace of F q n [3].Using generic complexity bounds on the resolution of polynomial systems, these authors provided attacks that can beat generic algorithms and can even have subexponential complexity for specific families of curves [2].At Eurocrypt 2012, Faugère, Perret, Petit and Renault re-analized Diem's attack [3] in the case F 2 n (denoted as FPPR method in this work), and showed that the systems arising from the Weil descent on Semaev's polynomials are much easier to solve than generic systems [7].Later at Asiacrypt 2012, Petit and Quisquater provided heuristic evidence that ECDLP is subexponential for that very important family of curves, and would beat generic algorithms when n is larger than about 2000 [15].In 2013, Shantz and Teske provided further experimental results using the so-called "delta method" with smaller factor basis to solve the FPPR system [7,20].
Even though these recent results suggest that ECDLP is weaker than previously expected for binary curves, the attacks are still far from being practical.This is mainly due to the large memory and time required to solve the polynomial systems arising from the Weil descent in practice.In particular, the experimental results presented in [15] for primes n were limited to n = 17.In order to validate the heuristic assumptions taken in Petit and Quisquater's analysis and to estimate the exact security level of binary elliptic curves in practice, experiments on larger parameters are definitely required.
In this paper, we focus on Diem's version of index calculus for ECDLP over a binary field of prime extension degree n [3,7,15].In that case, the Weil descent is performed on a vector space that is not a subfield of F 2 n , and the resulting polynomial system cannot be re-written in terms of symmetric variables only.We therefore introduce a different method to take advantage of symmetries even in the prime degree extension case.While Shantz and Teske use the same multivariate system as FPPR method [7,20], in this work we re-write the system with both symmetric and non-symmetric variables.The total number of variables is increased compared to [7,15], but we limit this increase as much as possible thanks to an appropriate choice of the vector space V .On the other hand, the use of symmetric variables in our system allows reducing the degrees of the equations significantly.Our experimental results show that our systems can be solved faster than the original systems of [7,15] as long as n is large enough.

Notations.
In this work, we are interested in solving the elliptic curve discrete logarithm problem on a curve E defined over a finite field F 2 n , where n is a prime number.We denote by E α,β the elliptic curve over F 2 n defined by the equation y 2 + xy = x 3 + αx 2 + β.For a given point P ∈ E, we use x(P) and y(P) to indicate the x-coordinate and y-coordinate of P respectively.From now on, we use the specific symbols P, Q and k for the parameters and solution of the ECDLP: P ∈ E, Q ∈ P , and k is the smallest non-negative integer such that Q = [k] P. We assume that the order of P is prime here.We identify the field , where h is an irreducible polynomial of degree n.Any element e ∈ F 2 n can then be represented as poly(e) For any set S, we use the symbol #S to mean the order of S. We denote the degree of regularity as D reg , which is the maximum degree appearing when solving the multivariate polynomial system with Gröbner basis routine.
Outline.The remaining of this paper is organized as follows.In Section 2, we recall previous index calculus algorithms for ECDLP, in particular FPPR method attack on binary elliptic curves and previous work exploiting the symmetry of Semaev's polynomials when the extension degree is composite.In Section 3, we describe our variant of FPPR method taking advantage of the symmetries even when the extension degree is prime.In Section 4, we provide experimental results supporting our method with respect to FPPR original attack.Finally in Section 5, we conclude the paper and we introduce further work.
Remark.This is a full version of the paper [10] published at the 8th International Workshop on Security (IWSEC 2013), held at Okinawa, Japan.

The index calculus method
For a given point P ∈ E α,β , let Q be a point in P .The index calculus method can be adapted to elliptic curves to compute the discrete logarithm of Q with respect to P.
As shown in Algorithm 1, we first select a factor base F ⊂ E α,β and we perform a relation search expressed as the loop between the line 3 and 7 of Algorithm 1.This part is currently the efficiency bottleneck of the algorithm.For each step in the loop, we compute R := [a] P+[b] Q for random integers a and b and we apply the Decompose function on R to find all tuples (sol m ) of m elements P j ∈ F such that P j 1 + P j 2 + • • • + P j m + R = O.Note that we may obtain several decompositions for each point R.In the line 6, the AddRelationToMatrix function encodes every decomposition of a point R into a row vector of the matrix M.More precisely, the first #F columns of M correspond to the elements of F, the last two columns correspond to P and Q, and the coefficients corresponding to these points are encoded in the matrix.In the line 8, the ReducedRowEchelonForm function reduces M into a row echelon form.When the rank of M reaches #F +1, the last row of the reduced M is of the form (0, A straightforward method to implement the Decompose function would be to exhaustively compute the sums of all m-tuples of points in F and to compare these sums to R. However, this method would not be efficient enough.

Semaev's polynomials
Semaev's polynomials [18] allow replacing the complicated addition law involved in the point decomposition problem by a somewhat simpler polynomial equation over Definition 1.The m-th Semaev's polynomial s m for E α,β is defined as follows: The polynomial s m is symmetric and has degree 2 m−2 with respect to each variable.Definition 1 provides a straightforward method to compute it.In practice, computing large Semaev's polynomials may not be a trivial task, even if the symmetry of the polynomials can be used to accelerate it [12].Semaev's polynomials have the following property: Proposition 1.We have s m (x 1 , x 2 , . . ., x m ) = 0 if and only if there exist y j ∈ F 2 n such that P j = (x j , y j ) ∈ E α,β and P 1 + P 2 + . . .
In his seminal paper [18], Semaev proposed to choose the factor base F in Algorithm 1 as where V is some subset of the base field of the curve.According to Proposition 1, finding a decomposition of and then finding the corresponding points A straightforward Decompose function using Semaev's polynomials is described in Algorithm 2.
In this algorithm, Semaev's polynomials are solved by a naive exhaustive search method.Since every x-coordinate corresponds to at most two points on the elliptic curve E α,β , each solution of s m+1 (x 1 , x 2 , . . ., x m , x(R)) = 0 may correspond to up to 2 m possible solutions in E α,β .These potential solutions are tested in the line 5 of Algorithm 2. As such, Algorithm 2 still involves some exhaustive search and can clearly not solve ECDLP faster than generic algorithms.

Algorithm 2 Decompose function with s m+1
Input: R = Output: sol m contains the decomposition elements of R w.r.t.F V

FPPR method
At Eurocrypt 2012, following similar approaches by Gaudry [9] and Diem [2,3], FPPR method provided V with the structure of a vector space, to reduce the resolution of Semaev's polynomial to a system of multivariate polynomial equations.They then solved this system using Gröbner basis algorithms [7].More precisely, FPPR method suggested to fix V as a random vector subspace of F 2 n /F 2 with dimension n .If {v 1 , . . ., v n } is a basis of this vector space, the resolution of Semaev's polynomial is then reduced to a polynomial system as follows.For any fixed P ∈ F V , we can write x(P ) as where c ∈ F 2 are known elements.Similarly, we can write all the variables x j ∈ V in s m+1 | x m+1 =x(R) as where c j, are binary variables and r ∈ F 2 are known.Using these equations to substitute the variables x j in s m+1 , we obtain an equation We have s m+1 | x m+1 =x(R) = 0 if and only if each binary coefficient polynomial f is equal to 0. Solving Semaev's polynomial s m+1 is now equivalent to solving the binary multivariable polynomial system The Decompose function using this system is described in Algorithm 3.

Algorithm 3 Decompose function with binary multivariable polynomial system (FPPR) [7]
Output: sol m contains the decomposition elements of R w.r.t.F V We first substitute x m+1 with x(R) in s m+1 .The Trans-FromSemaevToBinaryWithSym function transforms the equation s m+1 | x m+1 =x(R) = 0 into system f 1 , f 2 , . . ., f m as described above.To solve this system, we compute its Gröbner basis with respect to a lexicographic ordering using an algorithm such as F 4 or F 5 algorithm [4,5].A Gröbner basis of the system we solved here always contains some univariate polynomial (the polynomial 1 when there is no solution) with lexicographic ordering, and the solutions of f 1 , f 2 , . . ., f m can be obtained from the roots of this polynomial.However, since it is much more efficient to compute a Gröbner basis for a graded-reversed lexicographic order than for a lexicographic ordering, a Gröbner basis of f 1 , f 2 , . . ., f m is first computed for a graded-reverse lexicographic ordering and then transformed into a Gröbner basis for a lexicographic ordering using FGLM algorithm [6].
After getting the solutions of f 1 , f 2 , . . ., f m , we find the corresponding solutions over E α,β .As before, this requires to check whether P 1 + P 2 + . . .+ P m + R = O for all the potential solutions in the line 6 of Algorithm 3.
Although FPPR approach provides a systematic way to solve Semaev's polynomials, their algorithm is still not practical.Petit and Quisquater estimated that the method could beat generic algorithms for extension degrees n larger than about 2000 [15].This number is much larger than the parameter n = 160 that is currently used in applications.In fact, the degrees of the equations in f 1 , f 2 , . . ., f m grow quadratically with m, and the number of monomial terms in the equations is exponential in this degree.In practice, the sole computation of the Semaev's polynomial s m+1 seems to be a challenging task for m larger than 7.Because of the large computation costs (both in time and memory), no experimental result has been provided in [7] for n larger than 20.
In this work, we provide a variant of FPPR method that practically improves its complexity.Our method exploits the symmetry of Semaev's polynomials to reduce both the degree of the equations and the number of monomial terms appearing during the computation of a Gröbner basis of the system f 1 , f 2 , . . ., f m .

Use of symmetries in previous works
The symmetry of Semaev's polynomials has been exploited in previous works, but always for finite fields F p n with composite extension degrees n.The approach was already described by Gaudry [9] as a mean to accelerate the Gröbner basis computations.The symmetry of Semaev's polynomials has also been used by Joux and Vitse's to establish new ECDLP records for composite extension degree fields [12,13].Extra symmetries resulting from the existence of a rational 2-torsion point have also been exploited by Faugère et al. for twisted Edward curves and twisted Jacobi curves [8].In all these approaches, exploiting the symmetries of the system allows reducing the degrees of the equations and the number of monomials involved in the Gröbner basis computation, hence it reduces both the time and the memory costs.
To exploit the symmetry in ECDLP index calculus algorithms, we first rewrite Semaev's polynomial s m+1 with the elementary symmetric polynomials.Definition 2. Let x 1 , x 2 , . . ., x m be m variables, then the elementary symmetric polynomials are defined as Any symmetric polynomial can be written as an algebraic combination of these elementary symmetric polynomials.We denote the symmetrized version of Semaev's polynomial s m by s m .For example for the curve E α,β in characteristic 2, we have where x 3 is supposed to be fixed to some x(R).The elementary symmetric polynomials are The symmetrized version of s 3 is therefore Since x 3 is fixed and the squaring is a linear operation over F 2 , we see that symmetrization leads to a much simpler polynomial.
Let us now assume that n is a composite number with a non-trivial factor n .In this case, we can fix the vector space V as the subfield F p n of F p n .We note that all arithmetic operations are closed on the elements of V for this special choice.In particular, we have Let now {v 1 , v 2 , . . ., v n/n } be a basis of F p n /F p n .We can write where r ∈ F p n are known and the variables d j,0 are defined over F p n .These relations can be substituted in the equation s m+1 | x m+1 =x(R) = 0 to obtain a system of n/n equations in the m variables d j,0 only.Since the total degree and the degree of s m with respect to each symmetric variable σ i are lower than those of s m with respect to all non-symmetric variables x i , the degrees of the equations in the resulting system are also lower and the system is easier to solve.As long as n/n ≈ m, the system has a reasonable chance to have a solution.
Given a solution (σ 1 , . . ., σ m ) for this system, we can recover all possible corresponding values for the variables x 1 , . . ., x m (if there is any) by solving the system given in Definition 2, or equivalently by solving the symmetric polynomial equation Note that the existence of a non-trivial factor of n and the special choice for V are crucial here.Indeed, they allow building a new system that only involves symmetric variables and that is significantly simpler to solve than the previous one.

Using symmetries with prime extension degrees
When n is prime, the only subfield of F 2 n is F 2 , but choosing V = F 2 would imply to choose m = n, hence to work with Semaev's polynomial s n+1 which would not be practical when n is large.In Diem's and FPPR attacks [3,7], the set V is therefore a generic vector subspace of F 2 n /F 2 with dimension n .In that case, Implication (2) does not hold, but we now show how to nevertheless take advantage of symmetries in Semaev's polynomials.

A new system with both symmetric and non-symmetric variables
Let n be an arbitrary integer (possibly prime) and let V be a vector subspace of F 2 n /F 2 with dimension n .Let {v 1 , . . ., v n } be a basis of V .We can write where c j, with 1 ≤ j ≤ m and 1 ≤ ≤ n are variables but r , 1 ≤ ≤ n are known elements in F 2 .
Like in the composite extension degree case, we can use the elementary symmetric polynomials to write Semaev's polynomial s m+1 as a polynomial s m+1 in the variables σ j only.However since V is not a field anymore, constraining x j in V does not constrain σ j in V anymore.Since σ j ∈ F 2 n , we can however write where d j, with 1 ≤ j ≤ m and 1 ≤ ≤ n are binary variables.Using these equations, we can substitute σ j in s m+1 to obtain where f 1 , f 2 , . . ., f n are polynomials in the binary variables d j, .Applying a Weil descent on the symmetrized Semaev's polynomial equation s m+1 = 0, we therefore obtain a polynomial system f 1 = f 2 = . . .= f n = 0 in the mn binary variables d j, .The variables d j, must also satisfy certain constraints provided by System (1).More precisely, substituting both the x j and the σ j variables for binary variables in the equation we obtain where g j, are polynomials in the mn binary variables c j, only.In other words, applying a Weil descent on each equation of System (1), we obtain mn new equations d j, = g j, in the mn + mn binary variables c j, and d j, .The resulting system has mn + n equations in mn + mn binary variables.As before, the system is expected to have solutions if mn ≈ n, and it can then be solved using a Gröbner basis algorithm.
In comparison with the FPPR [7], the number of variables is multiplied by a factor roughly (m + 1).However, the degrees of our equations are also decreased thanks to the symmetrization, and this may decrease the degree of regularity of the system.In order to compare the time and memory complexities of both approaches, let D FPPR and D Ours be the degrees of regularity of the corresponding systems.The time and memory costs are respectively roughly #var 2D reg and #var 3D reg .Assuming that neither D FPPR nor D Ours depends on n (as suggested by Petit and Quisquater's experiments [15]), that D Ours < D FPPR (thanks to the use of symmetric variables) and that m is small enough, then the extra (m+1) factors in the number of variables will be a small price to pay for large enough parameters.In practice, experiments are limited to very small n and m values.For these small parameters, we could not observe any significant advantage of this variant with respect to FPPR.However, the complexity can be improved even further in practice with a clever choice of vector space.

A special vector space
In the prime degree extension case, V cannot be a subfield, hence the symmetric variables σ j are not restricted to V .This led us to introduce mn variables d j, instead of m variables only in the composite extension degree case.However, we point out that some vector spaces may be "closer to a subfield" than other ones.In particular if V is generated by the basis {1, ω, ω 2 , . . ., ω n −1 }, then we have More generally, we can write Applying a Weil descent on s m+1 | x m+1 =x(R) and each equation of System (1) as before, we obtain a new polynomial system + m variables.When m is large and mn ≈ n, the number of variables is decreased by a factor 2 if we use our special choice of vector space instead of a random one.For m = 4 and n ≈ 4n , the number of variables is reduced from about 5n to about 7n/2.For m = 3 and n ≈ 3n , the number of variables is reduced from about 4n to about 3n thanks to our special choice for V .In practice, this improvement turns out to be significant.
Table 1 is the comparison of different strategies used in the decomposition algorithm.Note that the degree of regularity is decreased from 7 to 4 when m = 3 by rewriting s m+1 to s m+1 with the symmetric function.It is difficult to estimate how many degrees of regularity are reduced for m other than 3 so far since we don't have enough experimental results due to the large polynomial system and the little resource.Our experimental results in section 4 implies the heuristic "D Ours < D FPPR " will be true for any m as long as s m+1 had simpler structure and smaller degree than s m+1 .The lack of the the data of degree of regularity for m > 3 makes the difficulty of the prediction of the degree of regularity in terms of m.This makes the complexity analysis following the step of [15] impossible even for a restricted model.If we make a model for a fixed m = 3, then the algorithm become more likely an exhaustive search instead of a sub-exponential algorithm.We will leave the estimation of the degree of regularity as a future work.

New decomposition algorithm
Our new algorithm for the decomposition problem is therefore using a new multivariate polynomial system by adopting the symmetric function and the special vector space V described above, denoted as ThisWork.The only difference between FPPR and ThisWork comes from a different TransFromSemaevToBinary function in the line 1 of Algorithm 3.Although the system solved in ThisWork contains more variables and equations than the system solved in FPPR, the degrees of the equations are smaller and they involve less monomial terms.We now describe our experimental results.

Experimental results
To validate our analysis and experimentally compare our method with FPPR, we implemented both algorithms in Magma.All our experiments were conducted on a CPU with four AMD Opteron Processor 6276 with 16 cores, running at 2.3 GHz with a L3 cache of 16 MB.The Operating System was LinuxMint 14 with 512GB memory.The programming platform was Magma V2.18-9 in its 64-bit version.Gröbner basis were computed with the Groeb-nerBasis function of Magma.Our implementations of FPPR and ThisWork share the same program, except the different TransFromSemaevToBinary function at line 1 of Algorithm 3. We first focus on the relation search, then we describe experimental results for a whole ECDLP computation.

Relation search
The relation search is the core of both FPPR and our variant.In our experiments, we considered a fixed randomly chosen curve E α,β , a fixed ECDLP with respect to P, and a fixed m = 3 for all values of the parameters n and n .For random integers a and b, we used both FPPR and ThisWork to find factor basis elements P j ∈ F V such that We focused on = 3 (fourth Semaev's polynomial) in our experiments.Indeed, there is no hope to solve ECDLP faster than with generic algorithms using m = 2 because of the linear algebra stage at the end of the index calculus algorithm a .On the other hand, the method appears unpractical for m = 4 even for very small values of n because of the exponential increase with m of the degrees in Semaev's polynomials.
The experimental results are given in Tables 2 and 3.For most values of the parameters n and n , the experiment was repeated 200 times and average values are presented in the table.For large values n = 6, the experiment was only repeated 3 times due to the long execution time.
We noticed that the time required to solve one system varied significantly depending on whether it had solutions or not.Tables 2 and 3 therefore present results for each case in separate columns.The table contains the following information: D reg is degree of regularity; t trans and t groe are respectively the time (in seconds) needed to transform the polynomial s m+1 into a binary system and to compute a Gröbner basis of this system; mem is the memory required by the experiment (in MB).
The experiments show that the degrees of regularity of the systems occurring during the relation search are decreased from values between 6 and 7 in FPPR to values between 3 and 4 in our method.This is particularly important since the complexity of Gröebner basis algorithms is exponential in this degree.As noticed in Section 3, this huge advantage of our method comes at the cost of a significant increase in the number of variables, which itself tends to increase the complexity of Gröbner basis algorithms.However, while our method may require more memory and time for small parameters (n, n ), it becomes more efficient than FPPR when the parameters increase.We remark that although the time required to solve the system may be larger with our method than with FPPR method for small parameters, the time required to build this system is always smaller.This is due to the much simpler structure of s m+1 compared to s m+1 (lower degrees and less monomial terms).Our method seems to work particularly well compared to FPPR when there is no solution for the system, which will happen most of the times when solving an ECDLP instance.

Whole ECDLP computation
In a next step, we also implemented the whole ECDLP algorithm with the two strategies FPPR and ThisWork.For the specified n, we ran the whole attack using m = 3 and several values for n .The orders of the curves we picked in our experiments are shown in Table 4 together with the experimental results for the best value of n , which turned out to be 3 in all cases.Timings provided in the table are in seconds.Table 4 clearly shows that ThisWork is more efficient than FPPR.It may look strange that n = 3 leads to optimal timings at first sight.Indeed, the ECDLP attacks described above use mn ≈ n and a constant value for n leads to a method close to exhaustive search.However, this is consistent with the observation already made in [7,15] that exhaustive search is more efficient than index calculus for small parameters.Table 5 also shows that while increasing n increases the probability to have solutions, it also increases the complexity of the Gröebner basis algorithm.This increase turns out to be significant for small parameters.

Conclusion and future work
In this paper, we proposed a variant of FPPR attack on the binary elliptic curve discrete logarithm problem (ECDLP).Our variant takes advantage of the symmetry of Semaev's polynomials to compute relations more efficiently.While symmetries had also been exploited in similar ECDLP algorithms for curves defined over finite fields with composite extension degrees, our method is the first one in the case of extension fields with prime extension degrees, which is the most interesting case for applications.
At Asiacrypt 2012, Petit and Quisquater estimated that FPPR method would beat generic discrete logarithm  algorithms for any extension degree larger than roughly 2000.We provided heuristic arguments and experimental data showing that our method reduces both the time and the memory required to compute a relation in FPPR, unless the parameters are very small.Our results therefore imply that Petit and Quisquater's bound can be lowered a little.
Our work raises several interesting questions.On a theoretical side, it would be interesting to prove that the degrees of regularity of the systems appearing in the relation search will not rise when n increases.It would also be interesting to provide a more precise analysis of our method and to precisely estimate for which values of the parameters it will become better than FPPR.
On a practical side, it would be interesting to improve the resolution of the systems even further.One idea in that direction is pre-computation of the invariant of this algorithm such as the transformation and the Gröbner basis of part of the system.In fact, even the resolution of the system could potentially be improved using special Gröebner basis algorithms such as F 4 trace [4,11].
Using Gröbner basis algorithms to solve ECDLP is a very recent idea.We expect that the index calculus algorithms that have recently appeared in the literature will be subject to further theoretical improvements and practical optimizations in a close future.

Endnote
a In fact, even m = 3 would require a double large prime variant of the index calculus algorithm described above in order to beat generic discrete logarithm algorithms [9].

Table 4 Comparison of ECDLP (m = 3, n' = 3)
*It is the product symbol which denoted the order of the EC group.