Keywords

1 Introduction

The index calculus algorithm originally denoted a technique to compute discrete logarithms modulo a prime number, but it now refers to a whole family of algorithms adapted to other finite fields and some algebraic curves. It includes the Number Field Sieve (NFS) [23], dedicated to logarithms in \(\mathbb {Z}_q\) and the algorithms of Gaudry [15] and Diem [8] for algebraic curves defined over \(\mathbb {F}_{q^n}\), where \(q=p^k\). Index calculus algorithms proceed in two main steps. The sieving (or point decomposition) step concentrates most of the number theory and algebraic geometry needed overall. By splitting random elements over a well-chosen factor base, it produces a large sparse matrix, the rows of which are “relations”. In a second phase, the matrix step produces “good” combinations of the relations by finding a non-trivial vector in the kernel of this matrix. This, in turn, enables the efficient computation of any discrete logarithm on the input domain. A crucial step of the index calculus on elliptic curves is to solve the point decomposition problem (pdp), by generating sufficiently many relations among suitable points on the curve. Using the so-called summation polynomials attached to the curve, this boils down to solving a system of polynomial equations whose solutions are the coordinates of points. The resulting algorithm has complexity \(O(q^{2-2/n})\), but this hides an exponential factor in n which comes from the hardness of solving the point decomposition problem.

Consequently, when q is large, \(n\ge 3\) is small and \(\log q>cm\) for some constant c, the Gaudry-Diem algorithm has a better asymptotic complexity than generic methods for solving the discrete logarithm problem and Gröbner basis algorithms have become a well-established technique [18] to solve these systems. Since a large number of instances of pdp needs to be solved, most of the research in the area has focused on improving the complexity of this step. Several simplifications such as symmetries and polynomials with lower degree obtained from the algebraic structure of the curve have been proposed [10].

When we consider elliptic curves defined over \(\mathbb {F}_{2^n}\) with n prime, solving the pdp system via Gröbner bases quickly becomes a bottleneck, and index calculus algorithms are slower than generic attacks, from a theoretical and a practical point of view. Moreover, it is not known how to define the factor base in order to exploit all the symmetries coming from the algebraic structure of the curve, without increasing the number of variables when solving pdp [36]. Finally, note that for random systems, pure Gröbner basis algorithms are both theoretically and practically slower than simpler methods, typically exhaustive search [6, 24], hybrid methods [2] and sat solvers. It is thus natural that we turn our attention towards combinatorics tools to solve the pdp in characteristic 2.

Until recent years, sat solvers have been proven to be a powerful tool in the cryptanalysis of symmetric schemes. They were successfully used for attacking secret key cryptosystems such as Bivium, Trivium, Grain, AES [16, 17, 22, 30, 31]. However, their use in public key cryptosystems has rarely been considered. A prominent example is the work of Galbraith and Gebregiyorgis [14], where they explore the possibility of replacing available Gröbner basis implementations with generic sat solvers (such as MiniSat), as a tool for solving the polynomial system for the pdp over binary curves. They observe experimentally that the use of sat solvers may potentially enable larger factor bases to be considered.

In this paper, we take important steps towards fully replacing Gröbner basis techniques for solving pdp with constraint programming ones. First, we model the point decomposition problem as a logical formula, with a reduced number of clauses, when compared to the model used in [14]. We compare different sat solvers and decide that the recently introduced WDSat solver [35] is most adapted to this problem and yields the fastest running times. Secondly, we propose a symmetry breaking technique and we implement it as an extension of this solver. We show that by using the extended solver, the proven worst-case complexity of solving a PDP is \(O(\frac{2^{ml}}{m!})\), where m is the number of points in the decomposition and l is the dimension of the vector space defining the factor base. This is to be compared against the Gröbner basis algorithm proposed in [11], whose runtime \(O(2^{\omega n/2})\) (with \(n\sim ml\) and \(\omega \) the linear algebra constant) is proven under heuristic assumptions.

We experimented with the index calculus attack on the discrete logarithm for elliptic curves over prime-degree binary extension fields. We obtain an important speedup in comparison with the best currently available implementation of Gröbner bases (F4 [11] in Magma [4]) and generic solvers [1, 31, 32]). Consequently, we were able to display results for a range of parameters l and n that were not feasible with previous approaches. In addition, our experiments show that Gröbner bases cannot compete with sat solvers techniques in terms of memory requirements. To illustrate, a system solved with the extended WDSat solver using only 17 MB of memory requires more than 200 GB when using the Gröbner basis method.

However, our experiments suggest that this improved pdp resolution does not render the index calculus attack faster than generic methods for solving the ECDLP in the case of prime-degree extension fields \(\mathbb {F}_{2^n}\).

This paper is organized as follows. Section 2 gives an overview of the index calculus algorithm on elliptic curves, introduces the pdp problem and briefly recalls algebraic and combinatorial techniques used in the literature to solve this problem. Section 3 details the logical models used in our experiments. Section 4 explains the symmetry breaking technique that we implemented in a sat solver. In Sect. 5 we give worst time complexity estimates for solving a pdp instance and derive the complexity of our sat-based index calculus algorithm. Finally, Sect. 6 presents benchmarks obtained with our implementation. We compare this against results obtained using Magma’s F4 implementation and several available best generic sat-solvers, such as MiniSat [32] and CryptoMiniSat [31].

2 An Overview of Index Calculus

In 2008 and 2009, Gaudry [15] and Diem [8] independently proposed a technique to perform the point decomposition step of the index calculus attack for elliptic curves over extension fields, using Semaev’s summation polynomials [27]. Since this paper focuses on binary elliptic curves, we introduce Semaev’s summation polynomials here directly for these curves.

Let \(\mathbb {F}_{2^n}\) be a finite field and E be an elliptic curve with j-invariant different from 0, defined by an equation

$$\begin{aligned} E: y^2+xy=x^3+ax^2+b, \end{aligned}$$
(1)

with \(a,b\in \mathbb {F}_{2^n}\). Using standard notation, we take \(\bar{\mathbb {F}}_{2^n}\) to be the algebraic closure of \(\mathbb {F}_{2^n}\) and \(E(\mathbb {F}_{2^n})\) (resp. \(E(\bar{\mathbb {F}}_{2^n})\)) to be the set of points on the elliptic curve defined over \(\mathbb {F}_{2^n}\) (resp. \(\bar{\mathbb {F}}_{2^n}\)). Let \(\mathcal {O}\) be the point at infinity on the elliptic curve. For \(m\in \mathbb {N}\), the m-th summation polynomial is a multivariate polynomial in \(\mathbb {F}_{2^n}[X_1,\ldots , X_m]\) with the property that, given points \(P_1, \ldots , P_m\in E(\bar{\mathbb {F}}_{2^n})\), then \(P_1 \pm \ldots \pm P_m=\mathcal {O}\) if and only if \(S_m(\mathbf{x} _{P_1}, \ldots , \mathbf{x} _{P_m})=0\). We have that

$$\begin{aligned}&S_2(X_1,X_2)=X_1+X_2,&\\&S_3(X_1,X_2,X_3)=X_1^2X_2^2+X_1^2X_3^2+X_1X_2X_3+X_2^2X_3^2+b,&\nonumber \end{aligned}$$
(2)

and for \(m \ge 4\) we have the following recursive formula:

$$\begin{aligned}&S_m(X_1,\ldots ,X_m)=&\\&Res_X(S_{m-k}(X_1,\ldots , X_{m-k-1},X),S_{k+2}(X_{m-k},\ldots , X_m,X)).&\nonumber \end{aligned}$$
(3)

The polynomial \(S_m\) is symmetric and has degree \(2^{m-2}\) in each of the variables. Let V be a vector subspace of \(\mathbb {F}_{2^n}/\mathbb {F}_2\), whose dimension l will be defined later. We define the factor basis \(\mathcal {B}\) to be:

$$\begin{aligned} \mathcal {B}=\{(\mathbf{x} ,\mathbf{y} )\in E(\mathbb {F}_{2^n})|\mathbf{x} \in V\}. \end{aligned}$$

Heuristically, we can easily see that the factor base has approximatively \(2^l\) elements. Given a point \(R\in E(\mathbb {F}_{2^n})\), the point decomposition problem is to find m points \(P_1, \ldots , P_m\in \mathcal {B}\) such that \(R=P_1 \pm \ldots \pm P_m\). Using Semaev’s polynomials, this problem is reduced to the one of solving a multivariate polynomial system.

Definition 1

Given \(s\ge 1\) and an l-dimensional vector subspace V of \(\mathbb {F}_{2^n}/\mathbb {F}_2\) and \(f\in ~\mathbb {F}_{2^n}[X_1,\ldots , X_m]\) any multivariate polynomial of degree bounded by s, find \((\mathbf{x} _1,\ldots , \mathbf{x} _m)\in ~V^m\) such that \(f(\mathbf{x} _1,\ldots ,\mathbf{x} _m)=0\).

Using the fact that \(\mathbb {F}_{2^n}\) is an n-dimensional vector space over \(\mathbb {F}_2\), the equation \(f(\mathbf{x} _1,\ldots ,\mathbf{x} _m)=0\) can be rewritten as a system of n equations over \(\mathbb {F}_2\), with ml variables. In the literature, this is called a Weil restriction [15] or Weil descent [26]. The probability of having a solution to this system depends on the ratio between n and l. Roughly, when \(n/l \sim m\) the system has a reasonable chance to have a solution.

Recent work on solving the decomposition problem has focused on using advanced methods for Gröbner basis computation such as Faugère’s \(F_4\) and \(F_5\) algorithms [11, 12]. This is a natural approach, given that similar techniques for small degree extension fields in characteristic \({>}2\) yielded index calculus algorithms which are faster than the generic attacks on the DLP.

A common technique when working with Semaev’s polynomials is to use a symmetrization process to further reduce the degree of the polynomials appearing in the pdp system. In short, since \(S_m\) is symmetric, we can rewrite it in terms of the elementary symmetric polynomials \(e_1=\sum _{1\le i_1 \le m}X_{i_1}\), \(e_2=\sum _{1\le i_1,i_2 \le m}X_{i_1}X_{i_2}\), \(\ldots \), \(e_m=\prod _{1\le i\le m} X_i\). We denote by \(S'_{m+1}\) the polynomial obtained after symmetrizing \(S_{m+1}\) in the first m variables, i.e. we have \(S'_{m+1}\in \mathbb {F}_{2^n}[e_1,\ldots , e_m,X_{m+1}].\)

In [36], the authors report on experiments carried on systems obtained using a careful choice of the vector space V and application of the symmetrization process. Using Magma’s \(F_4\) available implementation, we experimented with both the symmetric and the non-symmetric version for pdp systems and found, as in [36], that the symmetric version yields better results. Therefore, in order to set the notation, we detail this approach here.

Let t be a root of a defining polynomial of \(\mathbb {F}_{2^n}\) over \(\mathbb {F}_2\). Following [36], we choose the vector space V to be the l-dimensional subspace generated by \(1,t,t^2, \ldots , t^{l-1}\). Assuming that \(m(l-1) \le n\) we can write:

$$\begin{aligned} e_1= & {} d_{1,0}+\ldots +d_{1,l-1}t^{l-1} \nonumber \\ e_2= & {} d_{2,0}+\ldots +d_{2,2l-2}t^{2l-2} \\&\ldots&\nonumber \\ e_m= & {} d_{m,0}+\ldots +d_{m,m(l-1)}t^{m(l-1)}, \nonumber \end{aligned}$$
(4)

where the \(d_{i,j}\) with \(1\le i \le m\), \(0\le j\le i(l-1)\) are binary variables. After choosing \(\mathbf{x} _{m+1}\in \mathbb {F}_{2^n}\) and substituting \(e_1,\ldots , e_m\) as in Eq. (4), we get:

$$\begin{aligned} S'_{m+1}(e_1,\ldots , e_m,\mathbf{x} _{m+1})= & {} f_0+\ldots +f_{n-1}t^{n-1}, \end{aligned}$$

where \(f_{i}\), \(0\le i \le n-1\) are polynomials in the binary variables \(d_{i,j}\), \(1\le i \le m\), \(0\le j\le i(l-1)\). After a Weil descent, we obtain the following polynomial system

$$\begin{aligned} f_0=f_1=\ldots =f_{n-1}=0. \end{aligned}$$
(5)

One can see that with this approach, the number of variables is increased by a factor m, but the degrees of the polynomials in the system are significantly reduced. Further simplification of this system can be obtained if the elliptic curve has a rational point of order 2 or 4 [14]. Since this is a restriction, we did not implement this approach and used the system in Eq. (5) as the starting point for our sat model of the point decomposition problem.

2.1 Solving the Decomposition Problem Using SAT Solvers

Before presenting our approach for finding solutions to the pdp using sat solvers, we give preliminaries on the Satisfiability problem, its terminology and solving techniques. A sat solver is a special-purpose program to solve the sat problem. Using sat solvers as a cryptanalytic tool requires expressing the cryptographic problem as a Boolean formula in conjunctive normal form (cnf). The basic building block of a cnf formula is a literal, which is either a propositional variable or its negation. An or -clause is a non-exclusive disjunction (\(\vee \)) of literals \(x_1 \vee x_2 \vee \ldots \vee x_k\). A cnf formula is a unique or-clause or a conjunction (\(\wedge \)) of at least two or-clauses. An interpretation of a given propositional formula consists in assigning a truth value (\(\textsc {true} \)/\(\textsc {false} \)) to each of its variables. A cnf formula is said to be satisfiable if there exists at least one interpretation under which the formula is \(\textsc {true} \), and it is said to be unsatisfiable otherwise. The propositional satisfiability problem (sat) is the problem of determining whether a (usually cnf) formula is satisfiable.

In the remainder of this paper, we will refer to an or-clause simply by a clause, since cnf is the standard form used in sat solvers. A clause where the operation between literals is an exclusive or, will be referred to as a xor-clause. The use of the logical xor operator (\(\oplus \)) is common in cryptography. When working on cryptographic problems the cnf form can be extended to a cnf-xor form, which is a conjunction of both or-clauses and xor-clauses.

A straightforward method for solving the sat problem is to complete the truth table associated with the formula in question. This is equivalent to an exhaustive search method and thus impractical. Luckily, in some cases, a partial assignment on the set of variables can determine whether a clause is satisfiable. Assigning l, a literal from the partial assignment, to true will lead to:

  1. 1.

    Every clause containing l is removed (since the clause is satisfied).

  2. 2.

    In every clause that contains \(\lnot l\) this literal is deleted (since it can not contribute to the clause being satisfied).

The second rule above can lead to obtaining a clause composed of a single literal, called a unit clause. Since this is the only literal left that can satisfy the clause, it must be set to true and therefore propagated. The described method is called unit propagation. The reader can refer to [3] for more details.

A conflict occurs when it exists at least one clause with all literals assigned to false in the formula. If this case is a consequence of a direct assignment, or eventually of Unit Propagation, this has to be undone. This is commonly known as backtracking.

Example 1

For instance, these two atomic operations can be illustrated with the following example built of a set of 5 clauses numbered \(C_1\) to \(C_{5}\):

$$\begin{aligned}&C_1 : \lnot x_1 \vee x_2 \vee \lnot x_4 \\&C_2 : x_1 \vee x_3 \vee x_4 \\&C_3 : x_1 \vee \lnot x_3 \\&C_4 : x_1 \vee x_3 \\&C_5 : x_2 \vee x_4 \end{aligned}$$

Assigning the variable \(x_1\) to false leads the clause \(C_1\) to be satisfied by the literal \(x_1\). Another consequence is that the clauses \(C_2\), \(C_3\) and \(C_4\) cannot be satisfied by the literal \(x_1\). Hence, \(x_1\) can be deleted from these clauses. Then, \(C_3\) is a unit clause composed of the literal \(\lnot x_3\) and as a consequence, \(x_3\) has to be assigned to false. We say that the truth value of \(x_3\) is inferred through unit propagation.

When we set \(x_3\) to its inferred value false, we apply the second rule to clauses \(C_2\) and \(C_4\). As a consequence, clause \(C_4\) can not be satisfied by any of its literals. This constitutes a conflict and it invokes a backtracking procedure. The backtracking procedure consists in going back to the state that the formula was in before the last assumption was made. In our example, the last assumption was that \(x_1\) is false and thus, we go back to the initial state.

The basic backtracking search with unit propagation that we described composes the Davis-Putnam-Logemann-Loveland (dpll) algorithm [7], which is a state-of-the-art complete sat solving technique. dpll works by trying to assign a truth value to each variable in the cnf formula, recursively building a binary search tree of height equivalent (at worst) to the number of variables. After each variable assignment, the formula is simplified by unit propagation. If a conflict is met, a backtracking procedure is launched and the opposite truth value is assigned to the last assigned literal. If the opposite truth value results in conflict as well, we backtrack to an earlier assumption or conclude that the formula is unsatisfiable - when there are no earlier assumptions left. The number of conflicts is a good measure for the time complexity of a sat problem solved using a dpll -based solver. If the complete search tree is built, the worst-case complexity is \(O(2^v)\), where v is the number of variables in the formula. Figure 1 illustrates the binary search tree resulting from the resolution of Example 1.

Fig. 1.
figure 1

Binary search tree constructed with the dpll algorithm.

Full size image

A common variation of the dpll is the conflict-driven clause learning (cdcl) algorithm [29]. In this variation, each encountered conflict is described as a new clause which is learnt (added to the formula). State-of-the-art cdcl solvers, such as MiniSat [32] and Glucose [1], have been shown to be a powerful tool for solving cnf formulas. However, they are not equipped to handle xor-clauses and thus parity constraints have to be translated into cnf. Since handling cnf-clauses derived from xor constraints is not necessarily efficient, recent works have concentrated on coupling cdcl solvers with a xor-reasoning module. Furthermore, these techniques can be enhanced by Gaussian elimination, as in the works of Soos et al. (resulting in the CryptoMiniSat solver) [30, 31], Han and Jiang [17], Laitinen et al. [21, 22].

3 Model Description

This section gives in full detail the three models we used in our experiments: the algebraic one used by Yun-Ju et al. [36], the cnf model used by Galbraith and Gebregiyorgis [14] and the model we propose.

3.1 The Algebraic Model

Since the logical models are constructed starting from the algebraic one, we present first the model used when solving the pdp problem using Gröbner basis. The elementary symmetric polynomials \(e_i\) are written in terms of the \(d_{i,j}\) binary variables, as in Eq. (4). Similarly, since we look for a set of solutions \((\mathbf{x} _1,\ldots , \mathbf{x} _m)\in V^m\), the \(X_i\) variables are written formally as follows:

$$\begin{aligned} X_1=c_{1,0}+&\ldots +c_{1,l-1}t^{l-1}\\ X_2=c_{2,0}+&\ldots +c_{2,l-1}t^{l-1}\\&\ldots \\ X_m=c_{m,0}+&\ldots +c_{m,l-1}t^{l-1},\\ \end{aligned}$$

where \(c_{i,j}\), with \(1\le i \le m\), \(0\le j\le l-1\), are binary variables. Using Eq. (4), we derive the following equations:

$$\begin{aligned} d_{1,0}=c_{1,0} +&\ldots + c_{m,0} \nonumber \\ d_{1,1}=c_{1,1} +&\ldots + c_{m,1} \\&\ldots \nonumber \\ d_{m,m(l-1)}=c_{1,l} \cdot&\ldots \cdot c_{m,l}. \nonumber \end{aligned}$$
(6)

The remaining equations correspond to polynomials \(f_{i}\), \(0\le i\le n-1\), obtained via the Weil descent on \(S'_{m+1}\). Recall that these are polynomials in the binary variables \(d_{i,j}\). We now describe how we derive logical formulas from this system.

3.2 The CNF-XOR Model

When creating constraints from a boolean polynomial system, the multiplication of variables becomes a conjunction of literals and the sum of multiple terms becomes a xor-clause. From the two sets of equations in the algebraic model, we obtain two sets of xor-clauses, where the terms are single literals or conjunctions. To illustrate, the logical formula derived from Eq. (6) is as follows:

$$\begin{aligned}&\lnot d_{1,0} \oplus c_{1,0} \oplus \ldots \oplus c_{m,0} \nonumber \\&\lnot d_{1,1} \oplus c_{1,1} \oplus \ldots \oplus c_{m,1} \\&\ldots \nonumber \\&\lnot d_{m,m(l-1)} \oplus (c_{1,l} \wedge \ldots \wedge c_{m,l}). \nonumber \end{aligned}$$
(7)

sat solvers adapted for xor reasoning in the literature perform on xor clauses obtained by xoring single literals, and not conjunctions of several ones. To follow this paradigm, we have to transform the system above further. We substitute all conjunctions in a xor clause by a newly added variable. For example, let \(c'\) be the variable substituting a conjunction \((c_{i_1,j_1} \wedge c_{i_2,j_2} \wedge ... \wedge c_{i_k,j_k})\). We have \(c' \Leftrightarrow (c_{i_1,j_1} \wedge c_{i_2,j_2} \wedge ... \wedge c_{i_k,j_k})\), which rewrites as

$$\begin{aligned}&(c' \vee \lnot c_{i_1,j_1} \vee \lnot c_{i_2,j_2} \vee ... \vee \lnot c_{i_k,j_k})\, \wedge \nonumber \\&(\lnot c' \vee c_{i_1,j_1})\, \wedge \nonumber \\&(\lnot c' \vee c_{i_2,j_2})\, \wedge \\&\cdots \nonumber \\&(\lnot c' \vee c_{i_k,j_k})\, \nonumber \end{aligned}$$
(8)

For clarity, variables introduced by substitution of monomials containing exclusively the variables \(c_{i,j}\) will be denoted \(c'\) and clauses derived from these substitutions are said to be in the X-substitutions set of clauses. Similarly, substitutions of the monomials containing only the \(d_{i,j}\) variables are denoted by \(d'\) and the resulting set is referred to as the E-substitutions set of clauses.

After substituting conjunctions, we will refer to the set of clauses obtained from Eq. (7) as the E-X-relation set of clauses. Finally, the equations corresponding to polynomials \(f_{i}\), \(0\le i\le n-1\), are derived in the same manner and the resulting clauses will be referred to as the F set of clauses.

That concludes the four sets of clauses in our sat model. This model does not represent a cnf formula, since the E-X-relation set and the F set are made up of xor-clauses. Hence, it will be referred to as the cnf-xor model.

Proposition 1

Assigning all \(c_{i,j}\) variables, for \(1 \le i \le m\) and \({0\le j\le l-1}\), leads to the assignment of all variables in the cnf-xor model through unit propagation.

Proof

Let us examine the unit propagation process for each set of clauses separately.

  1. 1.

    Clauses in the X-substitutions set are obtained by transforming \(c' \Leftrightarrow (c_{i_1,j_1} \wedge c_{i_2,j_2} \wedge ... \wedge c_{i_k,j_k})\). We note that on the right of these equivalences there are only \(c_{i,j}\) variables and on the left, there is one single \(c'\) variable. The assignment of all of the \(c_{i,j}\) variables will yield the assignment of all variables on the left of the equivalences, i.e. all \(c'\) variables.

  2. 2.

    Clauses in the E-X-relations set are obtained by transforming the algebraic system in (6). We observe that on the right of the equations there are only \(c_{i,j}\) and \(c'\) variables and on the left there is one single \(d_{i,j}\) variable. When all \(c_{i,j}\) and all \(c'\) variables are assigned, all \(d_{i,j}\) variables will have their truth value assigned through unit propagation on the E-X-relation set.

  3. 3.

    Clauses in the E-substitutions set are obtained by transforming \(d' \Leftrightarrow (d_{i_1,j_1} \wedge d_{i_2,j_2} \wedge ... \wedge d_{i_k,j_k})\). Similarly as with the X-substitutions set, we have only \(d_{i,j}\) variables on the right of these equivalences and one single \(d'\) variable on the left. The assignment of all of the \(d_{i,j}\) variables will thus yield the assignment of all \(d'\) variables.

  4. 4.

    At this point, all variables in the parity constraints in the set F were assigned and we simply check whether the obtained interpretation satisfies the formula.

We conclude that variables in all four types of clauses of our CNF-XOR model were assigned through unit propagation.    \(\square \)

3.3 The CNF Model

Since most modern sat solvers read and process cnf formulas, we explain the classical technique for transforming a cnf-xor model to a cnf model. In fact, this is also the technique used in Magma’s available implementation for deriving a cnf model from a boolean polynomial system.

A xor-clause is said to be satisfied when it evaluates to true, i.e. when an odd number of literals in the clause are set to true and the rest are set to false. The cnf-encoding of a ternary xor-clause \((x_1 \oplus x_2 \oplus x_3)\) is

$$\begin{aligned} (x_1 \vee \lnot x_2 \vee \lnot x_3)\, \wedge \nonumber \\ (\lnot x_1 \vee x_2 \vee \lnot x_3)\, \wedge \\ (\lnot x_1 \vee \lnot x_2 \vee x_3)\, \wedge \nonumber \\ (x_1 \vee x_2 \vee x_3)\, \nonumber \end{aligned}$$
(9)

Similarly, a xor-clause of size k can be transformed into a conjunction of \(2^{k-1}\) or-clauses of size k. Since the number of introduced clauses grows exponentially with the size of the xor-clause, it is a good practice to cut up the xor-clause into manageable size clauses before proceeding with the transformation. To cut a xor-clause \((x_1 \oplus \ldots \oplus x_k)\) of size k in two, we introduce a new variable \(\varvec{x'}\) and we obtain the following two xor-clauses:

$$\begin{aligned}&(x_1 \oplus \ldots \oplus x_i \oplus \varvec{x'})\, \wedge \\&(x_{i+1} \oplus \ldots \oplus x_k \oplus \varvec{\lnot x'}). \end{aligned}$$

In our experiments with MiniSat in Sect. 6, we used a cnf model obtained after cutting into ternary xor-clauses, since any xor sat problem reduces in polynomial time to a 3-xor sat problem [3]. To the best of our knowledge, Magma’s implementation adopts a size 5 for xor clauses. The optimal size at which to cut the xor-clauses depends on the nature of the model and can be determined by running experiments using different values. Running these experiments was out of the scope of our work, as the WDSat solver does not use the cnf model.

We implemented all three models described in this section and we present Table 1 to serve as a comparison on the number of variables, equations and clauses. Values for the algebraic and cnf-xor model are exact, whereas those for the cnf model are averages obtained from experiments presented in Sect. 6. The value of m is always 3.

Table 1. The number of variables and equations/clauses for the three models.
Full size table

In 2014, Galbraith and Gebregiyorgis [14] used Magma’s implementation to compute the equivalent cnf logical formulas of the polynomial system resulting from the Weil descent of a pdp system and ran experiments using the general-purpose MiniSat solver to get solutions for these formulas. One can infer from Table 1 that the model they used has a significantly larger number of clauses and variables when compared to the cnf-xor model. This motivated our choice of the cnf-xor model for this work.

4 Breaking Symmetry

Since Semaev’s summation polynomials are symmetric, if \(\{\mathbf{x }_1, \ldots , \mathbf{x }_m\}\) is a solution, then all permutations of this set are solutions as well. These solutions are equivalent and finding more than one is of no use for the pdp. When a dpll -based sat solver is used (see Sect. 2.1), we observe redundancy in the binary search tree. Indeed, for \(m=3\) when a potential solution \(\{\mathbf{x }_1, \mathbf{x} _2, \mathbf{x} _3\}\) has been eliminated, \(\{\mathbf{x }_2, \mathbf{x} _1, \mathbf{x} _3\}\) does not need to be tried out. To avoid this redundancy, we establish the following constraint \(\mathbf{x} _1 \le \mathbf{x} _2 \le \ldots \le \mathbf{x} _m\), where \(\le \) is the lexicographic order on \(\{\textsc {false},\textsc {true} \}^l\) with false < true.

It would be tedious to add this constraint to the model itself, since this would imply adding new clauses and complexifying the sat model. Instead, we decided to add this constraint in the dpll algorithm using a tree-pruning-like technique. In a classical dpll implementation we try out both false and true for the truth value of a chosen variable. In our symmetry breaking variation of dpll, in some cases, the truth value of false will not be tried out as all potential solutions after this assignment would not satisfy the constraint \(\mathbf{x} _1 \le \mathbf{x} _2 \le \ldots \le \mathbf{x} _m\). Our variation of dpll is detailed in Algorithm 1 and the line numbers that distinguish it from a classical dpll algorithm are in bold. Note that one crucial difference between the two algorithms is the choice of a variable on line 4. While this choice is arbitrary in a classical dpll algorithm, in Algorithm 1 variables need to be chosen in the order from the leading bit of \(\mathbf{x} _1\) to the trailing bit of \(\mathbf{x} _m\). If this is not respected, our algorithm does not yield a correct answer.

figure a

Using the notation in Sect. 3, \(c_{i,j}\) corresponds to the jth bit of the ith \(\mathbf{x} \)-vector, where \(1 \le i \le m\) and \(0\le j\le l-1\). We recall from Proposition 1 that assigning all \(c_{i,j}\) variables in the cnf-xor model leads to the assignment of all variables through unit propagation. In Algorithm 1, we decide whether to try out the truth value of false for \(c_{i,j}\) or not by comparing two \(\mathbf{x} \)-vectors bit for bit, in the same way that we would compare binary numbers. When we are deciding on the truth value of \(c_{i,j}\) we have the following reasoning:

  • If \(c_{i-1,j}\) is false, we try to set \(c_{i,j}\) both to false and true (if false fails). When \(c_{i,j}\) is set to false, all of the potential \(\mathbf{x} _i\) solutions are greater than or equal to \(\mathbf{x} _{i-1}\), thus we continue with the same bit comparison on the next level. However, when \(c_{i,j}\) is set to true, all of the potential \(\mathbf{x} _i\) solutions are strictly greater than \(\mathbf{x} _{i-1}\) and we no longer do bit comparison on further levels.

  • If \(c_{i-1,j}\) is true, we only try out the truth value of true for \(c_{i,j}\) and we continue to do bit comparison since the potential \(\mathbf{x} _i\) solutions are still greater than or equal to \(\mathbf{x} _{i-1}\) at this point.

Lastly, we give further information which explains in full detail Algorithm 1. We use a flag denoted compare to instruct whether to do bit comparison at the current search tree level or not. On line 6 we reset the compare flag to true since \(c_{i,j}\), when \(j=0\), corresponds to a leading bit of the next \(\mathbf{x} \)-vector. Lastly, if-conditions on line 8 have to be checked in the specified order.

The assign procedure assigns the specified literal to true in a formula F, simplifies F and infers truth values for other literals. The backtrack procedure is used to undo all changes made to F after the last truth-value assignment. For more details on how these procedures are handled in the WDSat implementation, see [35].

5 Time Complexity Analysis

As we explained in Sect. 2, the time complexity of a sat problem in a dpll context is measured by the number of conflicts. This essentially corresponds to the number of leaves created in the binary search tree. The worst-case complexity of the algorithm is thus \(2^h\), where h is the height of the tree.

As per Proposition 1, we only reason on \(c_{i,j}\) variables from the cnf-xor model. Therefore, \(h=ml\) and the worst-case complexity for the pdp is \(2^{ml}\). Furthermore, using the symmetry breaking technique explained in Sect. 4, we optimize this complexity by a factor of m!. Indeed, out of the m! permutations of the solution set \(\{\mathbf{x }_1, \ldots , \mathbf{x} _m\}\), only one satisfies \(\mathbf{x} _1 \le \mathbf{x} _2 \le \ldots \le \mathbf{x} _m\) (neglecting the equality). This concludes that the worst-case number of conflicts reached for one pdp computation is

$$\begin{aligned} \frac{2^{ml}}{m!}. \end{aligned}$$
(10)

Going further in the time complexity analysis, we observe that to find one conflict we go through (in the worst case) all clauses in the model during unit propagation. Hence, the running time per conflict grows linearly with the number of clauses. First, let us count the number of clauses in the X-substitutions set. For every \(2\le d \le m\) there exist \(\left( {\begin{array}{c}m\\ d\end{array}}\right) \cdot l^d\) monomials of degree d given by products of variables \(c_{i,j}\), and they each yield \(d+1\) clauses (see Eq. (8)). In total, the number of clauses in the X-substitutions set is

$$\begin{aligned} (\sum _{d=2}^{m} \left( {\begin{array}{c}m\\ d\end{array}}\right) \cdot l^d)(d+1). \end{aligned}$$

Recall that degree one monomials are not substituted and thus do not produce new clauses. We can adapt this reasoning for the E-substitutions set as well.

The number of xor-clauses in the cnf-xor model is equivalent to the number of equations in the algebraic model. We have \(\frac{m(m+1)}{2}(l-1)+m\) in the E-X-relation set and n in the F set.

Remark 1

Using this analysis, we approximate the number of clauses, denoted by C, for \(m=3\), as all experiments presented in this paper are performed using the fourth summation polynomial.

$$\begin{aligned} C\approx & {} \left( {\begin{array}{c}3\\ 2\end{array}}\right) \cdot 3l^2 + \left( {\begin{array}{c}3\\ 3\end{array}}\right) \cdot 4l^3 + \left( \left( {\begin{array}{c}3\\ 2\end{array}}\right) \right) \cdot 3(3l-2)^2 + (6l-3) + n \approx \\ {}\approx & {} 4l^3+171l^2-210l+n+69. \nonumber \end{aligned}$$
(11)

In practice, many monomials have no occurrence in the system after the Weil descent. In fact, the value in Eq. (11) is a huge overestimate and exact values for \(l\in \{6,\ldots , 11\}\) are shown in Table 1.

Assuming that we take m small, we conclude that the number of clauses in our model is polynomial in l.

Let T be a constant representing the time to process one clause. The running time of the pdp is bounded by

$$\begin{aligned} C \cdot T \cdot 2^{ml}/m!. \end{aligned}$$

This allows us to establish the following result on the complexity of our SAT-based index calculus algorithm.

Theorem 1

The complexity of the index calculus algorithm for solving ECDLP on a curve defined over \(\mathbb {F}_{2^n}\), using a factor base given by a vector space of dimension l, is \(\tilde{O}(2^{n+l})\), where the \(\tilde{O}\) hides a polynomial factor in l.

Proof

In order to perform a whole ECDLP computation, one has to find \(2^l\) linearly independent relations. Following [9], the probability that a random point can be written as a sum of m factor basis elements is heuristically approximated by \(\frac{2^{ml}}{m!2^n}\). The time complexity for the full decomposition phase, using a DPLL-based solver coupled with the breaking symmetry technique is \(CT2^{n+l}\).   \(\square \)

This worst-case complexity is to be compared to the \(O(2^{\omega \frac{n}{2}+l})\) complexity of Faugère et al. [13]. Both approaches rely on the heuristic approximation of the probability that a random point can be decomposed in the factor base. However, we underline here that Faugère et al.’s proof of this result is based on an heuristic assumption on the Gröbner basis computation for pdp, while our analysis for the sat-based approach simply relies on the rigorously proved worst case for the dpll search tree (see Eq. (10)).

6 Experimental Results

We conducted experiments using \(S'_4\) on binary Koblitz elliptic curves [20] defined over \(\mathbb {F}_{2^n}\). We experimented with Gröbner bases and sat approaches. In [35], WDSat is reported to outperform the Gröbner basis methods, as well as all generic SAT solvers for this particular problem. First, we confirm this by experimenting with higher parameters and results are reported in Table 2. Secondly, we extend the WDSat solver with our symmetry breaking algorithm described in Sect. 4. Our symmetry breaking algorithm yields faster running times and we were able to perform experiments using greater parameters. Results are shown in Table 3. All tests were performed on a 2.40 GHz Intel Xeon E5-2640 processor. Our Weil descent implementation used to generate benchmarks is open source [34].

The Gröbner basis approach takes as input an algebraic model. We used the grevlex ordering, as this is considered to be optimal in the literature. The MiniSat solver processes a cnf model input, whereas CryptoMiniSat (CMS) and WDSat use the cnf-xor model. WDSat can also process directly an algebraic model in ANF form. Using the cnf-xor model is a huge advantage, as it has far fewer clauses and variables than the cnf model. Gaussian elimination can be beneficial for sat instances derived from cryptographic problems. However, it has been reported to yield slower running times for some instances, as performing the operation is very costly. For this reason, CryptoMiniSat and WDSat do not include Gaussian elimination by default, but the feature can be turned on explicitly. We experimented with both variants for both xor-able solvers.

With WDSat we set a custom order of branching variables, which allowed us to make use of unit propagation as explained in Proposition 1 and branch only on the \(c_{i,j}\) variables. CryptoMiniSat does not have this feature in the current version as the authors report that custom order of branching variables leads to slower running times in most cases. We added this feature to the source code of CryptoMiniSat and we ran tests both with a custom order as per Proposition 1 and with the order chosen by the solver.

Table 2 compares different approaches, showing results for optimal variants of each solving tool. Running times of all variants of CryptoMiniSat and WDSat are given in Appendix A. We experimented with different values of n for each l and we performed tests on 20 instances for each parameter size. Half of the instances have a solution and the other half do not. We show running time and memory averages on satisfiable and unsatisfiable instances separately since these values differ between the two cases. sat solvers stop as soon as they find a solution and if this is not the case they need to respond with certainty that a solution does not exist. Hence, running times of sat solvers are significantly slower when there is no solution. On the other hand, [36] indicates that the computational complexity of Gröbner bases is lower when a solution does not exist.

We set a timeout of 10 h and a memory limit of 200 GB for each run. Using MiniSat, we were not able to solve the highest parameter instances (\(l=8\)) within this time frame. On the other hand, Gröbner basis computations for these instances halted before timeout because of the memory limit. This data is in line with previous works. Indeed, [36] and [28] show experiments using the fourth summation polynomial with \(l=6\), whereas the highest parameter size achieved in [14] is \(l=8\).

Table 2 shows the average runtime in seconds, the average number of conflicts and the average memory use in MB. The WDSat solver allocates memory statically, according to predefined constant memory requirements. This explains why memory averages do not vary much between the different size parameters, or between satisfiable and unsatisfiable instances.

Table 2. Comparing different approaches for solving the pdp.
Full size table

Our experimental results show that performing Gaussian elimination on the system comes with a significant computational cost and yields a small decrease in the number of conflicts (see Table 4 in the Appendix). As this was the case for all instances derived from the Weil descent on \(S'_4\), we concluded that Gaussian elimination is not beneficial for this model. Choosing the WDSat variant without Gaussian elimination as optimal, we continued experiments for bigger size parameters using this variant coupled with the symmetry breaking technique. Table 3 shows results for \(l \in \{6,7,8,9,10,11\}\) and n sizes up to 89. All values are an average of 100 runs, as running times for satisfiable instances can vary remarkably. If we compare the number of conflicts for the first three values for l in this Table to that of the basic WDSat solver without the breaking symmetry extension in Table 2, we observe a speedup factor that rapidly approaches 6.Footnote 1 This confirms our claims in Sect. 5 that the symmetry breaking technique proposed in this paper yields a speedup by a factor of m!.

Comparing results for \(l=6\) and \(l=7\) in Table 3 with the equivalent results for the Gröbner basis method in Table 2, we observe that WDSat is up to 300 times faster than Gröbner bases for the cases where there is no solution and up to 1700 times faster for instances allowing a solution. This is a rough comparison, as the factor grows with parameters l and n.

Table 3. Experimental results using the complete WDSat solver. Running times are in seconds and memory use is in MB.
Full size table

Lastly, we experimented with the collision search [25] generic method, using the open source code at [33]. This implementation solves the discrete log problem in the case of prime field curves. We did not adapt the code for extension fields and the computation time for scalar multiplication on the curve might vary between the two cases. Even so, this allows for a rough comparison between the running times of generic methods and the work presented in this paper. In a uni-thread environment, a whole collision search computation for parameter \(n=59\) has an average runtime of 0.8 h on our platform. Computing \(2^l\) successful decompositions for parameters \(n=59\) and \(l=9\) would take more than 86 h according to results in Table 3. The estimated running time becomes considerably higher when we take into account unsuccessful decompositions as well. We conclude that for the case of prime-degree extension fields, even with the significant speedup that we achieved for the pdp, index calculus attacks are still not practical compared to the PCS generic method.

7 Conclusions and Future Work

Gröbner basis methods have been shown powerful in solving the pdp in the index calculus attack for elliptic curves defined over small degree extension fields in characteristic \({>}2\). In this paper, we argue that for finite fields in characteristic 2 a sat-based approach yields better results. We started by explaining that general-purpose sat solvers cannot yield considerably faster running times because the number of variables in a sat model is significantly larger than the number of variables in the algebraic model.

Our first contribution is to propose a pdp cnf-xor model with only ml core variables, whose assignment propagates all remaining variables in the model. To solve this model we use a sat solver dedicated to solving systems derived from a Weil descent. As our second contribution, we optimized the time complexity of this solver by a factor of m! using a symmetry breaking technique.

We presented experiments for the pdp on prime-degree extension fields in characteristic 2, using parameter sizes of up to \(l=11\) and \(n=89\). This presents a significant improvement over the current state-of-the-art, as experiments using \(l>8\) have never been shown before for this case. Moreover, memory is no longer a constraint for the pdp when the Gröbner basis computation is replaced with sat solving.

For technical reasons and lack of space, we were not able to provide here a complete comparison to other existing exhaustive search-based implementations, such as the libFes library [5] based on Bouillaguet et al.’s algorithm [6] and the Joux-Vitse hybrid algorithm [19]. For a more complete set of benchmarks, including experiments with Semaev’s polynomials for \(m > 3\), the interested reader is referred to the first author’s upcoming PhD thesis. It would also be interesting to test the performance of sat solvers on the simplified system obtained by considering the action of 2-torsion and 4-torsion points on the factor base, as in [14].