Philippe Mahey

Said Oualibouch

Tao Pham Dinh

SAID OUALIBOUCHt Philippe Maheyt

Pham Dinh Tao

Proximal Decomposition on the Graph of a Maximal Monotone Operator

Keywords: proximal point algorithm, partial inverse, convex programming AMS subject classification. 90C25

published or not. The documents may come   L'archive ouverte pluridisciplinaire

1. Introduction. We consider in this paper the following constrained inclusion problem: let X be a finite dimensional vector space and A a subspace of X. Let us denote by B the orthogonal subspace of A, i.e., B = A.L. Let T be a maximal monotone operator on X and denote its graph by Gr(T), i.e., Gr(T) = {(x, y) E X x X l y E Tx }. Then, the problem is to find x E A and y E B such that y E Tx, which can be written:

(P) Find (x, y) EX x X such that (x, y) E A x B n Gr(T).

A typical situation, which is easily shown to give the form (P), is the problem of minimizing a convex lower semicontinuous function on a subspace. The particular applications we have in mind are the decomposition methods for separable convex programming. They have recently gained some new interest with the possibility of implementing them on massively parallel architectures to solve very large problems such as the ones that appear in network optimization or stochastic programming (see [START_REF] Dinh Tao | Parallel and Distributed Computation[END_REF]). There are many different ways to transform a separable convex program in the form (P), but the general idea is to represent the coupling between the subsystems by a subspace of the product space of the copies of the primal and dual variables.

We are aiming here at the application of the Proximal Point Algorithm (PPA) (cf. [START_REF] Rockafellar | Monotone operotors and the proximal point algorithm in convex progrom ming[END_REF]) to problem (P). In 1983, Spingarn [START_REF] Spingarn | Partial inverse of a monotone operotor[END_REF] proposed a generalization of PPA to solve (P) that was based on the notion of the Partial Inverse operator. If we denote by x A the orthogonal projection of x on a subspace A, the graph of the partial inverse operator T A is given by Gr(TA) = {(xA + YB . YA + XB) I y E Tx }.

Applying the PPA to this operator leads to the Partial Inverse Method (PIM) which we summarize here. ALGORITHM 1 (PIM). At iteration k, (xk,Yk) E A X B. Then, find (x�, y� ) such that Xk + Yk = x� + y� and � (YDA + (YDB E T((x�)A + � (x�)B).

Then, (xk+l , Yk+l ) = (( x� )A, (YD B ) .

The main problem that arises with this algorithm is the difficulty of performing the proximal step (1) when c =11 in most interesting situations including the decom position methods. When c = 1, then the proximal step is a proximal decomposition on the graph of T and the subspaces A and B only appear in the projection step. In §3 we present the resultant algorithm, indeed equivalent to PIM with c = 1. The convergence is proved without the need to consider the Partial Inverse operator. The iteration is now written in the following way. It appears that the modified proximal step is uniquely determined and corresponds to a proximal decomposition on the graph of .XT if .Xp, = 1. We recover then the scaled version of PIM proposed by Spingarn in [START_REF]Applications of the method of partial inverse to convex progromming: Decomposition[END_REF]. It is mentioned in [START_REF] Idrissi | Applications and numerical convergence of the partial inverse method[END_REF] that the performance of PIM is very sensitive to the scaling factor variations and we give an explanation of this fact, allowing its adjustment to an optimal value in the strongly monotone case.

In §4, we give some numerical results that confirm the accelerating properties of the scaling parameter.

2. The proximal decomposition on the graph of a maximal monotone operator. We recall here some known results on the "P rox" mapping (I + T)-1 associated to a maximal monotone operator T and focus on the properties of the decomposition on the graph of T. More details on that subject can be found in [START_REF] Brezis | Operateurs Maximaux Monotones[END_REF] and [START_REF] Goebel | Topics in Metric Fixed Point Theory[END_REF].

Let T be a maximal monotone operator on a Hilbert space X. The graph of T, denoted by Gr(T), is defined by Gr(T) = {(x,y ) EX x Xiy E Tx}.

Monotonicity implies that for all x, x' E X and for all y E Tx, for all y' E Tx', (y-y', x-x') 2: 0. As T is maximal, its graph is not properly contained in the graph of any other monotone operator.

If T is strongly monotone, then there exists a positive p such that Vx,x'EX and VyETx, Vy'ETx', (y-y',x-x')2: Pii x-x' ii2 .

We say that the operator T is Lipschitz with constant L if Vx, x' EX and Vy E Tx, Vy' E Tx', ii Y-y' il � L il x-x' ll•

For monotone operators that share both properties, we get the following explicit bounds:

(3) Pii x -x' ll � II Y-y' ll � L ll x-x' ii•
When T is a linear operator represented by a positive definite matrix T, the best estimates for p and L are, respectively, the smallest and the largest eigenvalues ofT. Of course, if T is maximal monotone, then for any >. > 0, >.T is maximal monotone and if, moreover, Tis strongly monotone with modulus p and Lipschitz with constant L, then >.T is strongly monotone with modulus >. p and Lipschitz with constant >.L.

The resolvent associated with maximal monotone operator T is defined by (I + T)-1. It is single-valued, defined on the whole space, and firmly nonexpansive, which means that, if we let U =( I+ T)-1 and V =I-U, then, [START_REF] Eckstein | On the Douglas-Rachford splitting method and the proximal point algorithm for maximal monotone operotors[END_REF] or equivalently

(5) Vx,x' EX, ! l Ux-Ux'ii2 + II V x-Vx' ii2 � ll x-x' ll2 II Ux-Ux' ii2 � (x-x', Ux-Ux').
Related interesting facts on this characteristic property of resolvents may be found in theses by Martinet [9] and Eckstein [START_REF] Eckstein | Splitting Methods for Monotone Operators with Applications to Parallel Opti mization[END_REF] (see also [START_REF] Goebel | Topics in Metric Fixed Point Theory[END_REF]). Indeed, resolvents and maximal firmly nonexpansive mappings coincide and, following [START_REF] Lawrence | On fixed points of nonexpansive piecewise isometric mappings[END_REF], one-to-one correspondences among these operators, maximal monotone, and maximal nonexpansive operators, may be stated. This fact is explored further in the appendix.

We introduce now the proximal decomposition on the graph of a maximal mono tone operator.

Given a maximal monotone operator T and a vector (x,y) EX x X, there exists a unique pair (u, v) EX x X called the proximal decomposition of (x, y) on the graph of T such that u+v=x+y and (u, v)EGr(T).

The unicity is a direct consequence of the maximality of T and we get u = (I+ T)-1 (x + y), v =( I+ r-1 )-1 (x + y).

3. The proximal decomposition algorithm. We return now to problem (P), which has been analyzed by Spingarn [START_REF]Applications of the method of partial inverse to convex progromming: Decomposition[END_REF]. Let T be a maximal monotone operator on X. Let A be a subspace and B its orthogonal subspace. The problem is to find

(x, y) EX x X such that (x, y) E Ax B n Gr(T).
This problem is a particular case of the general problem of finding a zero of the sum of two maximal monotone operators. The algorithms we are aiming at are splitting methods that alternate computations on each operator separately (see [START_REF] Lions | Splitting algorithms for the sum of two nonlinear operotors[END_REF]). Indeed, most large-scale optimization problems can be formulated as the problem of mini mizing a separable convex lower semicontinuous function on a very simple subspace which represents the coupling between the subsystems.

We propose then a generic algorithm that alternates a proximal decomposition on the graph of T with a projection on Ax B. Before going on with the analysis of the method, we observe that the other alternatives that come to mind to find a point in the intersection of two sets are not suitable.

1. We can use the classical successive projections method on the two sets. The problem is that Gr(T) is generally not convex in X x X.

2. We cannot use another proximal decomposition on A x B (which is indeed the graph of the maximal monotone operator 8XA , the subdifferential of the indicator function of the set A), because it would lead back to the original point! Indeed, if (x,y) E Ax B and (u,v) is the proximal decomposition of x +yon Gr(T), then x = ( u + v) A and y = ( u + v) B, which means that ( x, y) is the proximal decomposition of u + von the graph of 8XA •

The Algorithm PDG (proximal decomposition on the graph) is stated below.

ALGORITHM 2 (PDG). Let (x0,yo) E Ax B. k = 0. If (xk, Yk) E Gr(T ), then stop: (xk, Yk) is a solution of (P). Else compute the proximal decomposition (uk,vk) of Xk + Yk on the graph of T. If (uk,vk) E Ax B, stop: (uk,vk) is a solution of (P). Else, Xk+ l = (uk)A and Yk+ l = (vk)B • k= k+ l An iteration of the algorithm may be formally stated as

(x,y) E A x B � C(x,y) = x+ y = z EX� (u,v) =:Fz � PAxs(u, v) E A x B,
where L is isometric from X x X into X, :F is the proximal decomposition operator from X into X xX, and PAxB is the projection on Ax B. Let us denote the composed mapping by .J = PAx B o:FoC.

We verify now that any fixed point (x, y) of Algorithm PDG is a solution of (P). Indeed, (x, y) = PAxs(u, v) and (u, v) = :Fz with z = x + y. If (u, v) E Ax B, then (x, y) is a solution of (P). Else, we have

(u-x, v-y) E L = {(a, b) EX x X i a + b = 0}.
But, as (x,y) = PAxs(u, v), we can state

(u -x,v-y)EB x A.
A and B being orthogonal subspace, the unique intersection of L and B x A is (0, 0). Thus, (x,y) = (u,v) and (x,y) solves (P).

On the other hand, if (x, y) is a solution of (P), :F(x + y) = (x, y), and (x, y) E A x B, which means that (x,y) is a fixed point of Algorithm PDG.

The PDG Algorithm is a particular instance of Spingarn's Partial Inverse Method [START_REF] Spingarn | Partial inverse of a monotone operotor[END_REF]. Indeed, when c = 1, the proximal step on the Partial Inverse operator TA becomes: Find (x�, YD such that : Xk + Yk = x� + y� and (YD A + (YOB E T( (x�) A + (xDB), which means, of course, that (x�, YD is the proximal decomposition of (xk, Yk) on the graph of T. Thus, the convergence has been established by Spingarn who has used the properties of the PPA applied to the partial inverse operator. However, here we give a direct proof of this fact that does not use the concept of the Partial Inverse. The main interest is that we shall obtain as a corollary the numerical analysis of the scaled version of PDG in the strongly monotone case.

We prove first that the composed mapping :T associated with Algorithm PDG is firmly nonexpansive. It can easily be seen that the mapping U = C o:J o£-1 is indeed the proximal operator associated to the Partial Inverse of T, i. e., U = (I+ TA)-1. But, we do not use this fact to prove that :Tis firmly nonexpansive.

THEOREM 3.1. The mapping :T associated to Algorithm PDG is firmly nonex pansive if and only if T is monotone. Moreover , it is defined on the whole space Ax B if and only if T is maximal monotone.

Proof. Let (x,y) and (x',y') E Ax B, z, z' E C(x,y),C (x',y') respectively, i.e., z = x + y and z' = x' + y', (u,v) E :F(z) and (u',v') E :F(z'), i. e., u + v = z, u = (I+ T)-1z and u' + v' = z', u' = (I + T)-1z'. Finally, let (uA, VB ) and ( uA., v k ) E A x B be the respective projections of ( u, v) and ( u', v') on A x B.

It is clear that, as z E (I+T)u, dom(:F) = R(I+T), and dom(.J) = c-1 (dom(:F)) = {(zA, ZB) E A X B i z E dom(:F)}.

E B (x-x', (u -• u')A} = (z-z', (u-u')A} = ( (u + v-u'-v', (u-u') A} = ( (u-u') + (v-v'), (u-u')A} and (y-y', (v-v')B} = (z-z', (v-v')B} = ( (u-u') + (v-v'), (v-v')B}• Hence, inequality (6) becomes ((u-u') + (v-v'), (u-u')A} + ( (u-u') + (v-v'), (v-v')B} 2:: ll (u-u')A II 2 + ll (v-v')B II 2•
We can now use the orthogonal decomposition of u -u' and v -v' on the direct sum A$B to get

'v'(u, v), (u', v') E Gr(T), ( (u-u') A, (v-v') A) + ( (u-u')B, (v-v ')B) � 0.
Finally, remarking that

(u-u', v-v') = ( (u-u')A, (v -v')A) + ( (u-u')B, (v-v')B),
we can conclude that . [START_REF] Lawrence | On fixed points of nonexpansive piecewise isometric mappings[END_REF]

D

Assuming that (P) has a solution, the convergence of the algorithm follows directly from Opial's lemma (see [START_REF] Opial | Weak convergence of the J'Uccessive approximations for nonexpansive mappings in Banach spaces[END_REF]), which states that, if a fixed point exists, a firmly nonexpansive operator is asymptotically regular and generates a convergent sequence. This is the very same idea as used by Martinet in the original proof for the PPA [START_REF] Martinet | Algorithmes pour la Resolution de Problemes d' Optimisation et de Minimax[END_REF] and developed further by Rockafellar who included approximate computations of the proximal steps [START_REF] Rockafellar | Monotone operotors and the proximal point algorithm in convex progrom ming[END_REF] .

4. A scaled decomposition on the graph of T. We introduce now a scaled version of the decomposition on the graph of a maximal monotone operator.

DEFINITION 4.1. Let (x, y) EX x X, T be a maximal monotone operator and>. a positive number . Then, the scaled proximal decomposition of (x, y) on the graph of T is the unique ( u, v) such that u + >.v = x + >.y, (u, v) E Gr(T).

Again, the existence and unicity of that new decomposition is a consequence of T being maximal monotone. Indeed, if v E Tu, we can write

u + >.v E u + >.Tu => u = (I+ >.T)-1 (u + >.v ) = (I+ >.T)-1 (x + >.y) v = .x-1 (x + >.y-u).
Observe that we can also write the following inclusions using the inverse operator r-1 for a given positive J.L:

u E T-1v, v+ J.Lu E v+ J.LT-1v, then v =(I+ J.LT-1 )-1 (v + J.LU).
Now, if J.L satisfies J.L-l = >. , we get v + J.LU = J.L(u + >.v) and, using the fact that (J.LT)-1z = r-1(J.L-1z), we obtain v = >. -l u + J.LT-1 )-1 (u + >.v) =(I+ J.LT-1 )-1 (J.L x + y).

Resuming, the scaled decomposition on the graph of T can be defined by (7) u = (I + >.T) -1 (x + >.y) , v = (I + J. .t T-1 )-1 ( J. .t X + y) , which appears as a natural generalization of (1 ). But, in fact, only one scaling factor can be introduced to maintain the desired properties, this is why we must fix AJ. .t = 1.

We can now describe the iteration of a scaled version of Algorithm PDG .

ALGORITHM 3 (SPDG) . (xk,Yk) E Ax B. Compute the scaled decomposition of (xk, Yk) on the graph of T.

U k = (I + >.T)-1 (xk + >.y k) , V k = >.-l (Xk + AYk -U k) • If (uk, vk) E Ax B, then stop. Else, Xk+l = (uk)A and Yk +l = (vk)B •
Observe that the scaled proximal decomposition can be stated in the following way.

Let w = >.v and r = >.y . Then, if ( u, v) is the scaled proximal decomposition of (x, y) on the graph of T, (u, w) is the proximal decomposition of (x, r ) on the graph of >.T. Hence, from the preceding section, we know that the sequence {(xk, rk)} converges to a point in Ax B nG r(>. T) . This fact implies that the sequence {(xk, Yk) } converges to a solution of (P ) .

On the other hand, we can see that SPDG is equivalent to the scaled version of the Partial Inverse Method (with c = 1) described by Spingarn in [13, Algorithm 2, p. 208 ) for the minimization of a convex function on a subspace. It reduces, of course, to PDG, i.e., to PIM, when>.= 1. Again, as the decomposition on the graph of T is a proximal step, approximate rules for computations can be added as in [START_REF] Rockafellar | Monotone operotors and the proximal point algorithm in convex progrom ming[END_REF] to get an implementable algorithm. We prefer to omit these details to focus on the accelerating properties of the scaling parameter, which constitute the main contribution of the present work.

To analyze the influence of the scaling parameter on the speed ratio of convergence of SPDG , we consider now the case where Tis both strongly monotone and Lipschitz. From the composed nature of .:1>.. and using the relations (8 ) and ( 9), we deduce the following bounds: Observe that, as

ll (uA,WB)-(u A ,w � )ll2 � ll u-u' ll2 + llw-w'll2 �li z-z'll2-2(u-u',w-w'} � li z-z'll2 -2p ll u-u'll2 � ( 1 -(1 � ).. : £ ) 2 ) li z-z'll2 � (1-(1 ��£ ) 2 ) li (x,
L � p,r(>.. ) = 2>-.p 1 -(1 + ).. £ ) 2 < 1.
We easily deduce the theoretical optimal value for >.. : [START_REF] Opial | Weak convergence of the J'Uccessive approximations for nonexpansive mappings in Banach spaces[END_REF] "X= 1/L and r("X) = J 1 -; L .

When Tis a linear positive definite operator, we observe that bad conditioning implies a slowdown of the algorithm. The optimal value of the scaling parameter must be chosen very small if L = J.Lmax, the largest eigenvalue of the associated matrix, is very large. We may observe that the speed ratio obtained in Theorem 4. 2 is the same as the one given in [START_REF] Lions | Splitting algorithms for the sum of two nonlinear operotors[END_REF] for the Douglas-Rachford splitting algorithm. Indeed, the connection between that algorithm and the Partial Inverse Method has been established by Eckstein [3 ] and we give its precise meaning in the Appendix.

The influence of the Lipschitz constant on the number of iterations has been analyzed for quadratic convex functions that were minimized on a simple subspace. The sensitivity to that parameter is shown on the five graphics of Fig. 1 and2 for different values of L, p, and the dimension of the space. These results are shown in Table 1. The influence of the scaling parameter on the number of iterations is illustrated by comparing columns iter(":\) ( number of iterations when ).. = "X) and iter (1) (number of iterations when ).. = 1). The number of iterations corresponds to the implementation of Algorithm PDG associated with the graph of >.. T. We show below why it is faster than the straightforward application of SPDG even if the primal sequences {xk} coincide in both algorithms.

It is also interesting to analyze the behaviour of the sequence {(xk,Yk)} and to look for some values of the scaling parameter such that, that sequence is mapped by a contraction. To be more precise, let .1>.. and 1-l>.. be the maps associated with the sequences {(xk, rk)} and {(xk, Yk)}, respectively. Then, if D>.. is the mapping defined by D>.. (x, y) = (x, >.. y), t I I � : --------T--------------�- ---------------------�-- As (xk, Yk) = D;.1 (xk, rk), we already know that the sequence {(xk, Yk)} converges when {(xk, rk)} converges. Note that a direct proof of this fact seems rather hard to state. The reason is that 1i>. is not necessarily a contractive map for any >.. We study below the conditions on >. to get a contraction in the strongly monotone case. In the strongly monotone and Lipschitz cases, we already know that 1i>. is a contraction for >. = 1 . The next theorem shows that this remains true if >. lies in a specific interval containing one.

I I I I � I \ \ \ dimension = 1 0 �i� - '---- c +---�1�-
' I I � : dimension = 1 o L=1 p= l I I I I I I \,r ---- 0�-
� \ 17' • ------ '-- • • • • -- • • ••- - • - - •• ••- •• ••-••••-•• ••-••••-••- -••-•• 0 �--�--------------------------� � . \ \ � \ t dimension = 1 00 p=O,l10S3S L=0
•--••-••-• • -••-• •-••-••-••-••-••--••-••-••-••-••-••--••-••-••-• oT-----��----�----------------�� i FIG.
2. Number of itemtions for dim=lOO.

THEOREM 4.3. Suppose that T is strongly monotone with modulus p and Lipschitz with constant L. Then, if A E (1, p + y'1 + p 2 ), the mapping 'H>.. is a contraction.

Proof Again let u =( I+ AT)-1 (x + AY) and u' =(I+ AT)-1(x' + Ay'). We use successively the nonexpansiveness of the projection and the firmly nonexpansiveness of the resolvent to write II'H>.. ( x, y ) -'H>.. (x',y') ll k $\\ (u-u 1 ,A-1 (x-x' + A(y-y')-u + u')) ll k $ A-211X -x'll2 + IIY-y'll2 + A 2-� p -1 11u-u'll2• Using the Lipschitz property, we obtain Hence, a sufficient condition that ensures that 'H.>.. is a contraction is A � 1 and O(A) = A 2-2Ap-1 < 0. That condition does not depend on the Lipschitz constant (indeed, this happens because 0 < p < L). We observe now that 0(1) = -2 p < 0 and the desired interval must be : A E (1, p + y'1 + p 2 ).

0

The different behaviour of both sequences {(xk, Zk)} and {(xk, Yk)} is illustrated in Fig. 3. For a small A, the second sequence (which is the one that will yield a solution for the original problem (P )) converges much slower even if it presents a monotonic decrease toward the fixed point. We conclude with the following observations on the choice of the scaling pa rameter: if bad conditioning is due to a too-small p, then we must accelerate the convergence by choosing .X close to the optimal value 1/ L (if it is not too far from 1!). If bad conditioning is due to a too-large L, then we may choose >. close to 1 in [ l, p+ Jl + p 2 ).

Appendix. The relation between the partial inverse and the Douglas-Rachford splitting operator may be explained in the following way which is directly inspired by the work of Lawrence and Spingarn [START_REF] Lawrence | On fixed points of nonexpansive piecewise isometric mappings[END_REF]. It was later derived by Eckstein and Bertsekas [START_REF] Eckstein | On the Douglas-Rachford splitting method and the proximal point algorithm for maximal monotone operotors[END_REF].

We recall the one-to-one correspondences among maximal monotone operators, maximal nonexpansive, and proximal operators as described in [START_REF] Lawrence | On fixed points of nonexpansive piecewise isometric mappings[END_REF].

Let a : (x, y) �---> (x, 2y-x) be the one-to-one correspondence of the class of proximal operators onto the class of nonexpansive operators and let f3 : (x, y) �---> (x + y,x-y) be the one-to-one correspondence of the class of monotone operators onto the class of nonexpansive operators. Following [START_REF] Lawrence | On fixed points of nonexpansive piecewise isometric mappings[END_REF], let us define two types of composition operations.

Let p 1 *P2 = a-1 (a(p 1 ) oa(p2 )) be the proximal operator obtained by composing two proximal operators p1 and p2 through their associated respective nonexpansive images (which give indeed another nonexpansive operator when composed). Likewise, let T1 8 T2 = /3-1 (/3(T1) o /3(T2 )) be the monotone operator obtained by composing two monotone operators in the same way. A straightforward calculus shows that, if p1 and P2 are the resolvents of T1 and T2, respectively, then p = Pl * P 2 is the resolvent of T = T1 8 T2 . As observed in [START_REF] Lawrence | On fixed points of nonexpansive piecewise isometric mappings[END_REF], we have the following interpretation of the * operation :

P1 * P 2 = P1 ° (2 p2 -I) +I -P 2, which is the operator associated to the fixed point iteration of the Douglas-Rachford splitting method (see [START_REF] Lions | Splitting algorithms for the sum of two nonlinear operotors[END_REF]). Observe that the nonexpansive operator a(p 1 ) o a(p2 )) is the operator associated with the Peaceman-Rachford iteration.

On the other side, it is shown in [START_REF] Lawrence | On fixed points of nonexpansive piecewise isometric mappings[END_REF] that, when T1 is the subdifferential mapping of the indicator function of a subspace A, i.e., Gr(T1 ) = A x Al., then T1 8 T = TA, the Partial Inverse of T. Resuming these facts, we have the following proposition.

PROPOSITION. Let T1 and T2 be two maximal monotone operators on X. The Douglas-Rachford splitting operator p = p1 o ( 2p 2 -I) +I-p2 , where p 1 = (I + >.T 1 ) -1 and p2 = (I+ >.T2 )-1 , is a proximal operator, indeed p = (I+ T)-1 , where T = >.T1 8 >.T2. Moreover, if Gr(T 1 ) =A x Al. and T2 = T, then p = (I + (>.T ) A)-1 , the resolvent of the partial inverse of >.T. Then, the Douglas-Rachford iteration applied to problem (P) is the partial inverse method associated to >.T. SPDG is the corresponding algorithm defined in the product space X x X.

Observation . Clearly (I+ (>.T )A)-1 f. (I+ >.TA)-1 • This point is crucial because the computation can only be performed in the first expression (this is then the SPDG Algorithm) or in the second expression with >. = 1.

2 )

 2 Proximal decomposition. Find the unique (x� , YD such that x� + y � = Xk +Yk and (x� , YD E Gr(T) If (x� , YD E A x B, then stop. Else (xk+t. Yk+ 1 ) = (( x�) A , (YD B ) • The unique solution of the proximal decomposition step is given by (1) X�= (I + T)-1 (xk + Yk) , Y� =(I + T-1 )-1 (xk + Yk) • Of course, only one proximal calculus is needed as (I + r-1 )-1 = I-(I + T)-1 . We propose then a modified proximal decomposition algorithm by introducing scaling factors .X and p,. Indeed, problem (P) may be written in two ways : y E Tx � x + .Xy E (I + .XT )x, x E T-1 y � y + p,x E (I + p,T -1 )y , which induces the following fixed point iteration, a natural scaled version of (1). (Modified proximal decomposition. X�= (I + .XT ) -1 (xk + AYk), y� =(I + J.L T-1 ) -1 (Yk + p,xk). If (x� , y�) E A x B, then stop. Else (xk+t.Yk+ 1 ) = (( x� ) A ,(YD B ) •

( 6 )

 6 Now, :Tis firmly nonexpansive if and only if 't/ (x,y), (x',y') E dom(:T ) and 't/ (uA, VB ) E :T (x,y), (u A, ,v k ) E :T (x',y') ( (x, y)-(x', y'), (uA, VB )-(u A, , v k )} 2:: ll (uA, VB )-(u A, , v k ) ll2• But, we have ( (x,y)-(x',y'), .( uA,VB)-(u A, ,v k)} = (x-x', (u-u')A} + (y-y', ( v-v')B} and, as x, x' E A and y,y '

THEOREM 4 . 2 .

 42 When Tis strongly monotone with modulus p and Lipschitz with constant L, then the convergence of the sequence {(xk, rk) } genemted by SPDG with rk = >.y k is linear with speed mtio 2..\p 1 -( 1 + >. L ) 2 . Proof. If .J>.. is the composed operator associated to the monotone operator >.T, we define as in Theorem 3.1 (x, r) , (x ', r' ) E Ax B, z = x+ r, z' = x' +r', (u , w), (u', w') E G r (>.T) . Then, (uA,WB) E .J>.. (x,r) and (u�,wB) E .J>.. (x ' , r' ) . The strong monotonicity of >.T implies that {8) Vw E Tu, w' E Tu', (w-w', u-u') ;::: >.p ii u-u'll2 and, as z E (I + >.T)u and z' E {I + >.T)u' (9) liz-z'll $ (1 + L) ii u-u' ll•

  r)-(x',r')ll2 • Let (x*, r*) be the limit point of the sequence { (xk, rk) }. It is therefore a fixed point of the mapping :1> ... Applying the above inequality to the pairs (xk+ l . rk+ l ) and (x*, r*), we obtain the desired result: ll (xk +l> rk+ l )-(x*,r*) ll2 � ( 1 -(1 � ).. : £ ) 2 ) ll (xk,rk)-(x*,r*) ll2 .D

�FIG. 1 .

 1 FIG. 1. Number of iterations for dim=lO.

  ,584036 \ \ ; -• mru j ...... . .... -............ _________ •

  is firmly nonexpansive if and only if T is monotone. Moreover, as dom(.7) = { (x, y) E A x Bl x + y E dom(F)}, we obtain .7firmly nonexpansive } { Tmonotone T . 1 t

	dom(.7) =A x B	<===>	dom(F) =X

<===>

max t ma mono one.

TABLE 1

 1 Numerical tests for quadmtic problems. Comparison of both sequences

	1/L	.A Iter( .A) Iter(1) dim	tolerance
	1.712	2.05	17	34
	0.936	1.05	17	20
	0.202	0.26	18	29
	0.102	0.13	18	42
	0.051	0.07	18	60
	0.034	0.05	18	70
	0.010	0.02	18	92	100	0.01
	0.508	0.58	17	19
	0.171	0.12	18	35
	0.094	0.12	17	44
	0.049	0.06	18	60
	0.033	0.05	18	70
	0.020	0.03	18	82
	0.010	0.02	19	92
	1.712	2.04	16	32
	0.936	1.21	16	18
	0.202 0.241	17	27
	0.102 0.121	17	39
	0.051 0.061	17	54
	0.034 0.041	17	63
	0.010 0.021	17	75	10	0.001
	0.508	0.58	15	16
	0.171 0.211	16	29
	0.094 0.121	16	40
	0.049 0.061	17	54
	0.033 0.041	17	63
	0.020 0.031	17	72
	0.010 0.021	17	75
						p=O ,l L=12 l=O,OS

' Iteration FIG. 3.


Claire Delaplace

Pierre-Alain Fouque

Paul Kirchner

Alexander May

Solving ECDLP over F p with Pre-computation via Representation Technique

We present a new algorithm solving the elliptic curve discrete logarithm problem (ECDLP) with pre-computation. Our method is based on recent work from Delaplace and May (NuTMiC 2019) which aims to solve ECDLP over some quadratic field F p 2 using the so-called representation technique. We first revisit Delaplace and May's algorithm and show that with better routines, it runs in time p, improving on the p 1.314 time complexity given by the authors. Then, we show that this method can be converted to a p 4/5 -algorithm, assuming that we have performed a p 6/5 -pre-computation step in advance. Finally, we extend our method to other fields F p , where ≥ 2 is a constant, leading to a p 2 /(2 +1) -algorithm, after a p ( +1)/(2 +1) -pre-computation step. Although our method does not improve on previous ECDLP with precomputation algorithm (namely the one from Bernstein and Lange), we think it offers an elegant alternative, which can hopefully lead to further improvements, the idea of using the representation technique for ECDLP being very new.

Introduction

Elliptic curves are important tools of today's cryptography. In particular, the elliptic curve discrete logarithm problem (ECDLP) is one of the key challenges of current public key cryptography. Given an elliptic curve E defined over a finite field F q and two points P and Q on this curve, this problem consists in recovering some integer k such that kP = Q.

This problem has been widely studied by the community in the last fifty years (e.g. [Sha71,Sem98,Gau09,GG16,PKM16,GWZ17]). We refer the reader to [GG16] for a general overview of existing ECDLP algorithms.

For elliptic curves defined over F p where p is a prime, Pollard's Rho algorithm obtains the best time and memory complexity, requiring O √ p operations over F p , and constant storage. This method basically works by creating a sequence (R i ) i≥1 of elements of P with R 0 = O, the point at infinity of the elliptic curve, and R i+1 = g(R i ) for a well chosen function g, such that each R i satisfies R i = a i P + b i Q for some a i and b i in Z. When a collision appears between two elements R i and R j , we have (a i -a j )P = (b j -b i )Q. Assuming that b i = b j , this leads to k = ai-aj bj -bi being the solution. Another interesting algorithm due to Pollard is the so-called kangaroo method. The algorithm computes two sequences (P i ) (the "tamed kangaroo") and (R i ) (the "wild kangaroo") such that each P i satisfies P i = a i P for some a i in Z and each R i satisfies R i = Q -b i P . When a collision appears between two elements P i and R j of these sequences, we have a i P = Q -b i P leading to k = a i + b i being the solution. All in all this procedure uses O √ p operations over F p and requires to store O √ p elements of F p . Many variants exists in various settings (e.g. [VOW99,FJM14,BL12,BL13]). In this paper, we are particularly interested in the one of [BL12,BL13] proposed by Bernstein and Lange, which studies the case of solving (EC)DLP with pre-computation.

In this setting, we assume that the elliptic curve E is given along with a generator P of a (large) subgroup of the elliptic group. This procedures recovers k such that Q = kP , for any input point Q in the group G generated by P . Note that here, we assume that the point P is independent of the instance. If we want to find k such that Q = kP for a different generator P belonging to G, we can find k 1 and k 2 such that Q = k 1 P , P = k 2 P . The solution k would then be such that k • k 2 P = k 1 P . In other words k = k 1 k -1 2 mod r, where r is the (well known) order of P . This does not change the overall complexity of the algorithm.

The idea of [BL12,BL13] is to pre-compute p λ "tamed kangaroos" sequences with different starting points and terminate them when a so called distinguished point is reached (c.f. [VOW99]). For each sequence, only the first and the last points are stored. The algorithm then works by computing the "wild kangaroo" sequence (R i ), until a distinguished point D is reached. If this D is also the ending point of one of the "tamed kangaroos" sequence (P i ), then a collision can be found by re-evaluating (P i ) and (R i ). Otherwise the procedure is restarted with a different starting point for the "wild kangaroo" sequence. All in all this procedure requires to perform O p 1-λ operations over F p and to store O p 1-λ elements of F p . The best time and memory complexity given is O p 1/3 , after a pre-computation which costs O p 2/3 . For elliptic curves defined over F p for some constant ≤ 2 and p being a prime or a power of a prime, it is possible to solve the problem using index calculus type of algorithms. This particular case has been widely studied by the community (e.g. [GS99,GHS02b,GHS02a,Gau09] ) The idea is then to split a point into an -sum of points from a factor basis that is defined over the field F p . Relying on Weil's descent and Semaev's summation polynomials [Sem04], this factor basis can be computed via a Gröbner basis computation. This yields a relation between factor basis elements. If at least q -the size of the factor base -such relations are collected, ECDLP can be solved via linear algebra. In total, Gaudry's method [Gau09] runs in time O p 2-2 . For = 2 this is the same complexity as Pollard's Rho algorithm. These methods seem quite hard to transpose to the setting where p is small (e.g. p = 2) and is a large prime extension due to the fact that some ! factor is hidden in the O (). However the idea of using index calculus based method to improve algorithms for ECDLP in characteristic 2 has also been investigated (e.g. [FPPR12,PQ12,HPST13,GG14]).

Recently, Delaplace and May [DM19] proposed another method to solve ECDLP with elliptic curves defined over a field F p 2 . Their algorithm relies on the representation technique introduced by Howgrave-Graham and Joux [HGJ10] to solve the subset-sum problem. In this particular context, the representation technique can break below the square-root barrier [HGJ10,BCJ11]. Unfortunately, the results of [DM19] were not as successful. The idea behind [DM19] is that the solution k can be split as the sum k = k 1 +k 2 for k 1 and k 2 having a specific "good shape" in p different ways, each of these p representation (k 1 , k 2 ) of k leading to the solution. The idea is then to compute only a p -1 fraction of all possible k 1 P and Q -k 2 P (e.g. points (x, y) such that x ∈ F p ). There should still be a pair (k 1 , k 2 ) satisfying k 1 P = Q -k 2 P which survives. In order to do so the authors of [DM19] split again k 1 into the sum k 11 + k 12 and k 2 into k 21 + k 22 , for k ij of "good shape" coming up with a 4-list algorithm. Indeed the authors of [DM19] first build lists L 11 = {k 11 P }, L 12 = {k 12 P }, L 21 = {Q-k 21 } and L 22 = {-k 22 }. Then, they construct the lists L and L such that

L = {(x, y) ∈ L 11 +L 12 , x ∈ F p } and L = {(x, y) ∈ L 21 + L 22 , x ∈ F p },
and search for a collision between these two lists. In order to perform the second step, they define a "join" problem, which given two lists L 1 and L 2 of points of

E(F p 2 ) such that |L 1 ||L 2 | = p 3/2 consists in returning a list L of all (x, y) = P 1 + P 2 with P 1 ∈ L 1 P 2 ∈ L 2 such that x ∈ F p .
They show that if this problem can be solved in time T and in memory M , then one can also solve ECDLP over F p 2 in time T and memory M . Delaplace and May proposed a way to solve this problem, relying on multipoint-evaluation of bivariate polynomials. Using results from [NZ04] they obtained a O p 1.314 algorithm.

Our contributions. We extend the work of Delaplace and May in several ways. First we present a Õ p 1+ -algorithm for the join problem defined in [DM19], relying on Kedlaya and Umans multivariate multipoint evaluation algorithm [KU08] rather than Nüsken and Ziegler. This reduces the time complexity of Delaplace and May's algorithm.

Second, we show that with a few modifications, this algorithm can be turned to an algorithm solving ECDLP for an elliptic curve defined over F p 2 in time and memory Õ (p c+ ), for any constant > 0, and for any 4/5 ≤ c ≤ 1 after a Õ p 2-c+ -pre-computation step. The core idea is that the lists L 11 , L 12 , L 22 as well as the join L of L 11 and L 12 do not depend on the target Q and can then be pre-computed in advance. For values c smaller than 4/5, the same algorithm can still be used with a larger memory consumption as the storage of the list L will dominate the memory complexity. The question of reducing the memory consumption in this case without increasing the runtime remains as an open problem.

Finally we extend this method to other finite fields F p , for ≥ 2 a constant. We obtain a Õ (p c+ ) in time and space algorithm, for any constant > 0, and for any 2 /(2 + 1) ≤ c ≤ /2 after a Õ p -c+ -pre-computation step. Once again, we could potentially use smaller values of c, but the memory consumption would then increase due to the space required to store the join list L.

Although our new algorithm does not outperform the one from Bernstein and Lange, we think it offers a refreshing alternative. Furthermore, it gives us reason to think that the use of representation technique for ECDLP is worth further investigation. Hopefully the results presented here can still be improved since the idea of using representation technique for ECDLP is genuinely new.

Organisation of the paper. After recalling the important notions and properties in Section 2, we revisit and improve the algorithm from [DM19] in Section 3. Then we propose our new ECDLP algorithm with pre-computation for curves over F p 2 in Section 4. Finally we extend this result to other fields F p , for a constant ≥ 2 and discuss the limitations and possible improvements in Section 5.

Preliminaries

For any integer r, we denote by Z r the ring Z/rZ. For any prime or power of a prime q, we denote by F q the finite field with q elements. Let be a constant, every x ∈ F p can uniquely be expressed as x = x 0 + αx 1 + . . . α -1 x -1 where α is a symbolic root of an irreducible polynomial P of degree over F p . The tuple (1, α, . . . , α -1 ) is a basis of F p seen as an F p -vector space of dimension . For any x > 0, we denote by log(x) the logarithm of x in basis 2. Let (G, +) be an abelian group and let S 1 and S 2 be two sets of elements of G. We denote by S 1 + S 2 the set of all x ∈ G such that x = x 1 + x 2 with x 1 ∈ S 1 and x 2 ∈ S 2 . Let L = ( 0 , . . . , N -1 ) be a list. We denote by |L| = N the size of L and for all 0 ≤ i < N , we denote the (i+1)-th entry of L by L[i] = i . If L 1 = ( 0 , . . . , N -1 ) and L 2 = (t 0 , . . . , t M -1 ) are two lists whose entries belong to an abelian group, we denote by L 1 + L 2 the list ( 0 + t 0 , . . . , 0 + t M -1 , . . . , N -1 + t 0 , . . . , N -1 + t M -1 ) An empty list is represented by the symbol ⊥. Let f 0 , . . . , f N -1 be a set of polynomials in F p [X 1 , . . . , X m ]. We denote by f 0 , . . . , f N -1 the ideal spanned by f 0 , . . . , f N -1 .

Elliptic curves and discrete logarithm problem. Let p > 3 be a prime, and let q = p . Let E be an elliptic curve over F q . We denote by O the point at infinity of E, and by E(F p ) the group formed by all the points of E, called the elliptic group. Similar to [DM19] in this paper we represent elliptic curves using Weierstraß's equations. In other words, let f : F q → F q be such that f (x) = x 3 + ax + b for some a, b ∈ F q . We have E(F q ) = {(x, y) ∈ F 2 q : y 2 = f (x)} ∪ {O}. We denote by |E(F q )| the order of the elliptic group. By Hasse's theorem, we now that |E(F q )| = Θ(q). Let P = O be a point of E(F q ). We denote by P the subgroup of E(F q ) which is generated by P , that is P is the set of all Q such that there exists k ∈ Z satisfiying Q = kP . We denote by r the order of P , which is also the number of elements in P .

Definition 1 (ECDLP). Let F q be a finite field and let E be an elliptic curve over F q . Let G be a large subgroup (possibly equal to) of prime order of E(F q ), and let P be a generator of G. We denote by r = O (q) the order of G. Given a point Q in G, solving the discrete logarithm problem over G consists in recovering k ∈ Z r such that kP = Q.

Similar to [DM19] we also call the particular case where q = p 2 where p is a prime the p 2 -ECDLP problem.

The complexity analyses of the algorithms we present later rely on the following heuristic which was already stated in [DM19]. Some theoretical justification can be found in recent work by Kim and Tibouchi [KT19].

Heuristic 1. Let E be an elliptic curve over a finite field F q . The x-coordinate of the points P = O of E(F q ) are uniformly distributed in F q .

Complexities. We use the standard Landau notations to describe the complexities presented in this paper. All time complexities are given in number of elementary operations (negations, additions, multiplications, inverse) over F p , and the memory complexity is the number of elements of F p that we need to store. In particular, this means that for any constant , all elementary operations over F p and the sum of two points in E(F p ) can be performed in time O (1), and storing an element of F p or a point of E(F p ) requires O (1) memory. We also use the Õ (•) notation, to hide poly-logarithmic factors in p. In other words, Õ (X) = O (X(log p) c ), where c is a constant.

Representation technique. The representation technique was introduced in

[HGJ10] by Howgrave-Graham and Joux to solve the subset-sum problem. In this context, given a vector a ∈ Z n , and a target t ∈ Z, the goal is to recover a vector e ∈ {0, 1} n , so that the dot product between a and e satisfies a, e = t. Howgrave-Graham and Joux noticed that if e consists of exactly n/2 ones, then e can be represented as the sum of two vectors e 1 + e 2 such that e 1 and e 2 consist of n/4 ones in n/2 n/4 ways. Finding only one of these representations (e 1 , e 2 ) of e is all that is needed to solve the problem. As such, additional constraints can be enforced on the dot products a, e 1 and a, e 2 so that only one of these representations survives with high probability. This yields to the first subset-sum algorithm breaking below the square-root bound with a complexity of 2 0.337n . In [DM19], Delaplace and May investigated whether a similar method can be used to solve the p 2 -ECDLP problem. We detail and revisit their idea in Section 3. For now, we recall the definition of representation as it was given in

[DM19]. Definition 2. Let t > 0, r ∈ N and k ∈ Z r . Let S 1 , . . . , S t be subsets of N. A tuple (k 1 , . . . , k t ) ∈ S 1 × • • • × S t is a representation of k if k = k 1 + • • • + k t (mod r).
Useful results. Here, we recall some lemmas which we utilise to support our complexity claims. The first one is the Schwartz-Zippel Lemma, which estimates the probability that a random point in F m p is a zero of a given polynomial in m variables.

Lemma 1 (Schwartz-Zippel Lemma). Let F ∈ F p [X 1 , . . . , X m ] be a polynomial of total degree at most d. Let (x 1 , . . . , x m ) be a randomly chosen point in F m p , the probability that (x 1 , . . . x m ) is a zero of F satisfies

P[F(x 1 , . . . , x m ) = 0] ≤ d p .
We also rely on an algorithm due to Kedlaya and Umans [KU08] which evaluates simultaneously a polynomial

F ∈ F p [X 1 , . . . X m ] in N points of F m p .
The bit complexity of this algorithm was originally studied in [KU08]. These results have recently been refined and extended to other computational models (namely Turing machine) by Van der Hoeven and Lecerf [vdHL19]. In our complexity model however, the difference between the analysis of [KU08] and [vdHL19] disappears in the Õ ().

Lemma 2 ). Let F be a polynomial in m variables over a prime field F p , with individual degree less than d. For any fixed constant , it is possible to evaluate simultaneously

F in N points of F m p in time Õ (d m + N ) 1+ .
Finally we also rely on Weil-descent to transform a polynomial equations g(x 1 , . . . , x n ) = 0 over F p to a system of polynomial equations over F p .

Lemma 3 (Weil descent, stated as in [GG16] lemma 1). Let (1, α, . . . , α -1 ) be a basis of F p seen as a vector space over F p . Let g(x 1 , . . . , x m ) ∈ F p [x 1 , . . . , x m ]. We denote by y the tuple y = (y 1,0 , . . . , y 1, -1 , . . . , y m,0 , . . . , y m, -1 ). Then there exist unique polynomials g

i (y) ∈ F p [y], 0 ≤ i < such that f(y 1,0 + • • • + α -1 y 1, -1 , . . . , y m,0 + • • • + α -1 y m, -1 ) = -1 i=0 α i f i (y). Furthermore if f(a 1 , . . . , a m ) = 0, for some a 1 , . . . , a m ∈ F p , there exist b k,j ∈ F p , 1 ≤ k ≤ m, 0 ≤ j < , such that a k = -1 j=0 b k,j α j and for all i, f i (b) = 0, where b = (b 1,0 , . . . , b 1, -1 , . . . , b m,0 , . . . , b m, -1 ).

Revisiting Delaplace and May Algorithm

We focus here on solving the p 2 -ECDLP problem. We consider an elliptic curve E given by its Weierstraß equation, a large subgroup G of E(F p 2 ) of prime order r generated by a point P and a target point Q ∈ G such that Q = kP , for some unknown k ∈ Z r . Let (1, α) be a base of F p 2 seen as a F p -vector space.

Algorithm 1 DM19-ECDLP((a, b), P, Q, f)

Require: An elliptic curve E over F p 2 defined by the coefficients (a, b) of its Weierstraß's equation, a generator P of the group G, a target point Q and the polynomial f associated to E.

Ensure: k ∈ Zr such that kP = Q or ⊥. 1: L11 ← {k11P : k11 ∈ S11} 2: L12 ← {k12P : k12 ∈ S12} 3: L21 ← {Q -k21P : k21 ∈ S21} 4: L22 ← {-k22P : k22 ∈ S22} 5: L ← Join((a, b), L11, L12, f) 6: L ← Join((a, b), L21, L22, f) 7: for all R = Q -k2P ∈ L do 8: if R = k1P ∈ L then 9:
return k1 + k2 mod r 10: return ⊥

General Idea

In [DM19] Delaplace and May proposed a new algorithm using representation technique. The main idea is that the solution k can be represented by a pair

(k 1 , k 2 ) such that k 1 (resp. k 2 ) belongs to a certain well chosen subset S 1 (resp. S 2 ) of N, and (k 1 , k 2 ) satisfies k = k 1 + k 2 mod r.
The sets S 1 and S 2 are constructed so that there are at least r/p ≈ p representations (k 1 , k 2 ) of k. It is then possible to reduce the search space by computing only those which satisfy a certain criteria, which is verified by a p -1 fraction of all the points of E(F p 2 ). Their algorithm is recalled in Algorithm 1 and illustrated in Figure 1.

To summarise, the authors of [DM19] did the following. They computed the lists L and L such that

L = {k 1 P = (x, y) : k 1 ∈ S 1 , x ∈ F p }, L = {Q -k 2 P = (x , y ) : k 2 ∈ S 2 , x ∈ F p },
and searched for a collision k 1 P = Q -k 2 P between L and L . They argued that, as long as p is large enough, we can expect one collision to exist.

The difficulty lies in building the lists L and L efficiently. Delaplace and May proposed the following 4-list method.

P : k 12 ∈ S 12 } L 21 = {Q -k 21 P : k 21 ∈ S 21 } L 22 = {-k 22 P : k 22 ∈ S 22 } Join Join (x, y) x ∈ Fp (x , y ) x ∈ Fp (k 11 + k 12 )P = Q -(k 21 + k 22 )P L11 k11P L12 k12P L21 L22 Q -k21P -k22P √ p p √ p Fig. 1: [DM19] 4-list Algorithm.
The list L is then constructed as the join between L 11 and L 12 . This means that L is the sublist of L 11 +L 12 such that for all (x, y) in L, x ∈ F p . Accordingly, L is constructed as the join between L 21 and L 22 . The authors of [DM19] argued that any two points P 1 = (x 1 , y 1 ), P 2 = (x 2 , y 2 ) such that P 1 + P 2 = P 3 = (x 3 , y 3 ) with x 3 ∈ F p , satisfy a certain polynomial equation

f(u 1 , v 1 , u 2 , v 2 ) = 0,
where x i = u i + αv i for i ∈ {1, 2} and where f is a polynomial of constant degree obtained from Weierstraß's equation. We give more detail on how this polynomial is constructed in section 4.1. From here, they designed an algorithm for the join which is summarised in Algorithm 2. To compute the join of two lists, say L 1 and L 2 , the idea is to highlight a sublist of L 2 , such that all points (x 2 , y 2 ) satisfy

(x1,y1)∈L1 f(u 1 , v 1 , u 2 , v 2 ) = F(u 2 , v 2 ) = 0. (1)
The authors of [DM19] propose to evaluate simultaneously F in all (u 2 , v 2 ) such that (x 2 , y 2 ) ∈ L 2 in order to find this sublist. This step is described by the MultipointEvaluation procedure in Algorithm 2. Then, for each of this points

Algorithm 2 Join((a, b), f, L 1 , L 2 )
Require: An elliptic curve E over F p 2 defined by the coefficients a and b of its Weierstraß's equation, two lists of points L1, L2, the polynomial f associated to E.

Ensure: L = {(P0, k1 + k2) : P0 = P1 + P2 = (x, y), (Pi, ki) ∈ Li, i ∈ [1, 2], x ∈ Fp}. 1: L ← ⊥ 2: A ← {(u1, v1) : u1 + αv1 = x1, (x1, y1) ∈ L1} 3: B ← {(u2, v2) : u2 + αv2 = x2, (x2, y2) ∈ L2} 4: A f ← {fu 1 ,v 1 ∈ F 2 p [U, V ] : fu 1 ,v 1 (U, V ) = f(u1, v1, U, V )} 5: F ← A f fu 1 ,v 1 6: E ← MultipointEvaluation(F, B) 7: B ← {B[i] : E[i] = 0} 8: for all fu 1 ,v 1 ∈ A f do 9:
for all (u2, v2) ∈ B do 10:

if fu 1 ,v 1 (u2, v2) = 0 then 11: L ← L ∪ {P1 + P2 : Pi = (xi, yi) ∈ Li, xi = ui + αvi, i ∈ {1, 2}} 12: return L P 2 ,
they search all corresponding points P 1 ∈ L 1 such that (x, y) = P 1 + P 2 implies x ∈ F p . These points P 1 are recovered by exhaustive search.

We also have the following result from Delaplace and May.

Lemma 4 ([DM19] Theorem 1). If Algorithm 2 runs in time T and uses memory M then it is possible to solve p 2 -ECDLP in time mostly T using mostly M as memory.

Note that we do not actually need to store the list L . When a point R = Q -k 2 P is computed during the join procedure, we can directly check whether R ∈ L. If so the solution can be returned immediately.

Improving the Join Runtime

Delaplace and May used an algorithm from Nüsken and Ziegler [NZ04] for the MultipointEvaluation procedure. Given a bivariate polynomial F of individual degrees lower than some d, Nüsken and Ziegler's algorithm evaluates this polynomial simultaneously in N points in time O (N + d 2 )d ω2/2-1+ , for any fixed > 0 and where ω 2 is the linear algebra exponent of the multiplication of an n-by-n matrix by a n-by-n 2 one. The authors of [DM19] showed that this step dominates the whole runtime, and for input lists L 1 and L 2 of respective size √ p and p, this lead to a runtime T = O p 1.314 -replacing ω 2 by the best known bound (cf. [LG14]) -for the whole Join procedure, and thus for their p 2 -ECDLP algorithm.

Instead of Nüsken and Ziegler's results, we propose to use Kedlaya and Umans's algorithm whose complexity is recalled in Lemma 2. This leads to the following result. Proof. This proof is very similar to the one given in [DM19] Theorem 2. Let > 0 be a fixed constant. First we claim that steps 2 to 4 can be performed in time which is proportional to the input lists L 1 and L 2 , that is in O (|L 2 |) operation over F p . Using fast polynomial multiplication algorithms, it is possible to compute F in time Õ |L 1 | 2 . We now need to evaluate the bivariate polynomial F of individual degrees smaller than some

d = O (|L 1 |) in |L 2 | points of F 2 p . From Lemma 2, it can be done in time Õ (|L 1 | 2 + |L 2 |) 1+ , since |L 2 | = p c

New Complexity Claim

With this improved join procedure, we can reduce the complexity claim of [DM19].

We have the following result. Once again, we assume that we are given an elliptic curve E defined by its Weierstraß's equation

y 2 = f (x) = x 3 + ax + b, a subgroup G = P of prime order r = O p 2 , and a target point Q in G.
The idea is that we can pre-compute lists L 11 , L 12 , L 21 , L 22 such that

L 1j = {k 1j P : k 1j ∈ S 1j } (2)
L 2j = {-k 2j P : k 2j ∈ S 2j } i ∈ {1, 2},
for some well chosen subset S ij of N. It is also possible to pre-compute the join L between L 11 and L 12 , since it does not depend on Q.

The algorithm would then take as input L, L 21 and L 22 , the point Q and the polynomial f, and perform the following steps.

1. For all points P 21 ∈ L 21 replace P 21 by P 21 + Q. 2. Compute the join L between L 21 and L 22 . 3. Search for a collision between L and L .

In this section, we will show that the following result holds.

Theorem 2. For any δ > 0, Algorithm 5 solves p 2 -ECDLP in time Õ p 2c2+δ , using Õ max(p 3c1-1 , p 2c2+δ ) memory, after a pre-computation which takes time Õ p 2c1 .

The memory consumption of our algorithm is minimised when c 1 = 3 5 and c 2 = 2 5 . In this case, our algorithm runs in time Õ p 4/5+δ using Õ p 4/5+δ memory after a pre-computation step having runtime Õ p 6/5+δ .

Pre-computation

There are basically three steps during the pre-computation: (1) compute the polynomial f; (2) compute the lists L ij as given in Equation 2 for some well chosen sets S ij , i, j ∈ {1, 2}; (3) compute the join L between L 11 and L 12 . These steps are summarised in Algorithm 3.

Algorithm 3 Precomputation((a, b), P )

Require: An elliptic curve E over F p 2 defined by the coefficients (a, b) of its Weierstraß's equation y = f (x), a generator P of G. Ensure: Lists L, L21, L22, such that the L2j lists are defined as in Equation 2and L = {k1P = (x, y) : k1 ∈ S11 + S12, x ∈ Fp} and the polynomial f.

1: G ← ((X1-X2) 2 (X1+X2+X3)-f (X1)-f (X2)) 2 -4f (X1)f (X2) in F p 2 [X1, X2, X3] 2: g0, g1 ← Weil-Descent(G) in Fp[U1, V1, U2, V2, U3, V3] 3: g0 ← g0(U1, V1, U2, V2, U3, 0) in Fp[U1, V1, U2, V2, U3] 4: g1 ← g1(U1, V1, U2, V2, U3, 0) in Fp[U1, V1, U2, V2, U3] 5: f ← ResU 3 (g0, g1) 6: L11 ← {k11P : k11 ∈ S11} 7: L12 ← {k12P : k12 ∈ S12} 8: L21 ← {-k21P : k21 ∈ S21} 9: L22 ← {-k22P : k22 ∈ S22} 10: L ← Join((a, b), L11, L12, f) 11: return L, L21, L22, f
More sophisticated join. If we used the procedure given in Algorithm 2 to perform this last step, then the time complexity of the join would be dominated by the exhaustive search at the end. Indeed this search takes a time which is proportional to the product of |L 11 | and |L 12 |, where L 12 is the list of points in L 12 which pass the multipoint-evaluation test. We claim that, even in this case, it is still possible to perform the join in time Õ |L 12 | 1+ , assuming that |L 11 | 2 = |L 12 | = p 2c , for any constant c < 1. More precisely, for any constant γ such that 2 γ-1 -1 2 γ-1 < c ≤ 2 γ -1 2 γ , there is an algorithm computing the join of two lists L 1 and L 2 such that |L 1 | 2 = |L 2 | = p 2c . The Join procedure given in Algorithm 2 is a particular case of this procedure when γ = 1.

For the sake of simplicity, we present an algorithm which does the job only for 1/2 < c ≤ 3/4 in Algorithm 4.

The idea is to first build the list B of all points (u 2 , v 2 ) in B such that there is a polynomial f u1,v1 in A f such that f u1,v1 (u 2 , v 2 ) = 0. Then, for a well chosen parameter κ, we split

A f into p κ buckets A (i) f , such that |A (i) f | = p c-κ for all i. The we compute p κ sublists B (i) of B , such that for all (u 2 , v 2 ) ∈ B (i) , the product F i of all f u1,v1 ∈ A (i) f , satisfies F i (u 2 , v 2 ) = 0.
We re-iterate the procedure, with the A (i) f and B (i) lists, as long as required (if 2 γ-1 -1 2 γ-1 < c ≤ 2 γ -1 2 γ , then γ steps are to be expected).

We have following lemma.

Lemma 6. Let L 1 and L 2 be such that

|L 1 | 2 = |L 2 | = p 2c
, where c < 1 is a constant. For any fixed > 0, it is possible to compute the join between L 1 and

L 2 in time Õ |L 2 | 1+ . Algorithm 4 Join((a, b), f, L 1 , L 2 ) for |L 1 | ≤ p 3/4
Require: An elliptic curve E over F p 2 defined by the coefficients a and b of its Weierstraß's equation, two lists of points L1, L2, such that |L1| = p c |L2| = p 2c , 0 < c ≤ 3/4 the polynomial f associated to E.

Ensure: L = {(P0, k1 + k2) : P0 = P1 + P2 = (x, y), (Pi, ki) ∈ Li, i ∈ [1, 2], x ∈ Fp}. 1: L ← ⊥ 2: A ← {(u1, v1) : u1 + αv1 = x1, (x1, y1) ∈ L1} 3: B ← {(u2, v2) : u2 + αv2 = x2, (x2, y2) ∈ L2} 4: A f ← {fu 1 ,v 1 ∈ F 2 p [U, V ] : fu 1 ,v 1 (U, V ) = f(u1, v1, U, V )} 5: for all 1 ≤ i ≤ p κ 1 do 6: A (i) f ← A f [(i -1)p 2c-1 ..ip 2c-1 -1] 7: Fi = A (i) f fu 1 ,v 1 8: F ← i Fi 9: E ← MultipointEvaluation(F, B) 10: B ← {B[i] : E[i] = 0} 11: κ1 ← 1 -c 12: for all 1 ≤ i ≤ p κ 1 do 13: E (i) ← MultipointEvaluation(Fi, B) 14: B (i) ← {B [j] : E (i) [j] = 0} 15: for all fu 1 ,v 1 ∈ A (i) f do 16: for all (u2, v2) ∈ B (i) do 17: if fu 1 ,v 1 (u2, v2) = 0 then 18: L ← L ∪ {P1 + P2 : Pi = (xi, yi) ∈ Li, xi = ui + αvi, i ∈ {1, 2}} 19: return L
For simplicity we prove this lemma for the case 1/2 < c ≤ 3/4. As we will see in Section 4.3, this corresponds to the case where the memory consumption of our Algorithm is the lowest (c = 3/5).

Proof. It is already clear from Lemma 5 that Steps 1 to 10 of the Join procedure described in Algorithm 4 takes time Õ |L 2 | 1+ . The rest of the algorithm consists mostly in computing p κ1 polynomial multipoint evaluations and an exhaustive search at the end.

We first estimate the cost of the multipoint evaluations. Using Kedlaya-Umans algorithm, this would be

Õ p κ1 (|A (i) f | 2 + |B |) 1+ ,
for some constant > 0. Furthermore for all i, we have

|A (i) f | = p 2(c-κ1) = p 2(2c-1) , since κ 1 = 1 -c. We also have |B | = p 3c-1 . The reader can check that |B | ≥ |A (i) f |. It follows that the time complexity of the multipoint evaluation is Õ p κ1+(3c-1)(1+ ) = Õ p 2c+(3c-1) . We choose ≤ (3c -1) /2c so that this complexity is Õ |L 2 | 1+ .
The runtime of the exhaustive search at the end is given by

O p κ1 |A (i) f ||B (i) | .
Under heuristic 1, we can estimate |B (i) | thanks to Schwartz-Zippel lemma to be

|B (i) | = O |A (i) f ||B | p = O p 5c-3 .
It follows that the time complexity of this step is

O p κ1+5c-3+c-κ1 = Õ p 6c-3 . Having c < 3/4, implies that p 6c-3 < |L 2 |.
Using the same argument, the proof can be extended to all 2 γ-1 -1

2 γ-1 < c ≤ 2 γ -1 2 γ
, where γ is a constant positive integer, by choosing a sequence of parameters κ 0 , . . . , κ γ-1 , such that κ i = i(1 -c).

We are now ready to prove this following statement.

Lemma 7. For any constant , it is possible to perform the whole pre-computation in Õ |L 12 | 1+ operations over F p .

Proof. First, we claim that computing the polynomial f can be done in constant time. Then, using a similar argument than the one used by [DM19] to prove Lemma 4, we show that computing the join takes more time than creating the base lists.

Let us recall that any point

(x, y) of E, y = f (x). Let G be the polynomial in F p 2 [X 1 , X 2 , X 3 ] such that G(X 1 , X 2 , X 3 ) = ((X 1 -X 2 ) 2 (X 1 + X 2 + X 3 ) -f (X 1 ) -f (X 2 )) 2 -4f (X 1 )f (X 2 ).
If P 1 = (x 1 , y 1 ), P 2 = (x 2 , y 2 ) and P 3 = (x 3 , y 3 ) are such that P 1 + P 2 = P 3 , then G(x 1 , x 2 , x 3 ) = 0. Let (1, α) be a base of F p 2 seen as a F p -vector space of dimension 2. From Weil Descent, there exist polynomials g 0 and g

1 in F p [U 1 , V 1 , U 2 , V 2 , U 3 , V 3 ] such that G(U 1 + αV 1 , U 2 + αV 2 , U 3 + αV 3 ) = g 0 (U 1 , V 1 , U 2 , V 2 , U 3 , V 3 ) + αg 1 (U 1 , . . . , V 3 ).
G is a polynomial in 3 variables of constant individual degree (i.e. deg X1 (G) = deg X2 (G) = 6 and deg X3 (G) = 2), as such, computing g 0 and g 1 can be done in constant time, and the individual degrees of g 0 and g 1 are constants. As in [DM19], we consider polynomials

g i (U 1 , V 1 , U 2 , V 2 , U 3 ) = g i (U 1 , V 1 , U 2 , V 2 , U 3 , 0) for i ∈ {0, 1}.
and we compute the resultant f of g 0 and g 1 in U 3 . As both g i polynomials consist of a constant number of variables and have constant individual degrees, this can be performed in constant time.

As 1: Finally we need to compute the join between L 11 and L 12 . From Lemma 6, this can be performed in Õ (|L 12 |) 1+ for any chosen constant > 0.

L21 ← L21 + {Q} 2: L ← Join((a, b), L21, L22, f) 3: for all Q -k2P ∈ L do 4: if Q -k2P = k1P ∈ L

The Main Algorithm

Given Q and the precomputed lists L, L 21 , L 22 we want to recover k such that Q = kP . The method is summarised in Algorithm 5. Proof. The first step of the algorithm consists in replacing all points P 3 of L 21 by Q + P 3 , this can be done in O (|L 21 |) elementary operations over F p . For now on, we consider that L 21 is the list of all Q -k 21 P .

The second step consists in building the join L between L 21 and L 22 . For any constant > 0, using Kedlaya-Umans to perform the multi-evaluation step, it is possible to compute L in time

Õ |L 22 | 1+ , provided that |L 22 | = |L 21 | 2 .
Finally the last step consists in searching for a collision between the lists L and L . We can assume that L is sorted according to the x-coordinate of its points (sorting L could have been done during the pre-computation step in time O (|L|) under Heuristic 1). Under Heuristic 1, it is possible to check if a point R = (x, y), with x ∈ F p belongs to L in constant time. As such this whole step can be performed in time O (|L |). This gives us a total time complexity of

Õ |L 21 | + |L 22 | 1+ + |L | . We have |L 21 | 2 = |L 22 | as such |L 21 | < |L 22 |.
Furthermore, let us denote by L 22 the sublists of L 22 of all points (x 2 , y 2 ) with x 2 = u 2 + αv 2 such that there is

(x 1 , y 1 ) ∈ L 21 with x 2 = u 2 + αv 2 and f(u 1 , v 1 , u 2 , v 2 ) = 0. Under Heuristic 1,
we can show that the expected size of L is

E[|L |] = O |L 21 ||L 22 | p .
This comes from the fact that if the x-coordinates of the points of L 21 + L 22 are uniformly distributed, then a point (x, y) ∈ L 21 + L 22 will satisfy x ∈ F p with probability p -1 . Having

|L 21 | < p ensures that E[|L |] < |L 22 |.
Thus the whole time complexity of this algorithm is

T = Õ |L 22 | 1+ .
We now prove the claimed memory complexity. This algorithm requires to store the input lists L, L 21 and L 22 , as well as the list L . Some additional lists have to be stored during the join procedure, but their size is bounded by Õ |L 22 | 1+ . We have already argued that max(|L |, |L 21 |) ≤ |L 22 |. It follows that the memory required by this algorithm is

Õ max(|L|, |L 22 | 1+ ) ,
which concludes the proof.

Complexity Analysis

Choice of the S ij sets. Now we discuss how to define the sets S ij . Let 1 2 ≤ λ ≤ 1 be a parameter. For simplicity, we assume that p λ is an integer. If not, replacing p λ by its rounding to the closest integer, the argument still holds. We consider the following subsets of N.

S 1 = {k 1 < r : k 1 = s 1 p + s 0 , 0 ≤ s 0 < p λ } (3)
S 2 = {k 2 < r : k 2 = s 1 p + s 0 , 0 ≤ s 0 < p and p λ |s 0 } (4)

The shape of the elements from the S i sets is given in Figure 3a. The sets S ij are built so that

S i1 ∩ S i2 = {0}, S i1 + S i2 = S i and |S i1 | 2 = |S i2 | for i ∈ {1, 2}. More precisely let c 1 = 1 3 (λ + 1) and c 2 = 1 3 (2 -λ).
It is possible to check that the following holds:

1 2 ≤ c 1 ≤ λ ≤ 1 ≤ λ + c 2 ≤ 4 3 , since, 1 2 ≤ λ ≤ 1.
We construct the sets S ij as follows.

S 11 = {k 11 ∈ S 1 : k 11 < p c1 } S 12 = {k 12 ∈ S 1 : k 12 = s 1 p + s 0 , p c1 |s 0 } S 21 = {k 21 ∈ S 2 : k 21 < p λ+c2 } S 22 = {k 22 ∈ S 2 : k 22 = sp, p λ+c2-1 |s}
The shape of the elements from the S ij sets is given in Figure 3b. Lemma 9. The sets S i1 and S i2 satisfies S i1 ∩ S i2 = {0} S i1 + S i2 = S i , and

|S i1 | 2 = O (S i2 |), for i ∈ {1, 2}.
Proof. Let us first prove the relations between the sizes. By construction

|S 11 | = p c1 k 12 < r and 0 ≤ s 0 < p λ , |S 12 | = r p • p λ-c1 = O p 1+λ-c1 . Recall that c 1 = 1 3 (1 + c 1 ), it follows that |S 12 | = O p 2c1 .
We also have

|S 21 | = p p λ • p λ+c 2 p = p c2 , and |S 22 | = r p • 1 p c 2 +λ-1 = O p 2-λ-c2 . Recall that c 2 = 1 3 (2 -λ), it follows that |S 12 | = O p 2c2 .
We prove that the intersection is reduced to {0}. It is clear that 0 belongs to all sets. Let now be k 12 = s 1 p + s 0 p c1 an element of S 12 . We show that if k 12 ∈ S 11 , then k 12 = 0. Indeed k 12 ∈ S 11 means that k 12 < p c1 this implies that s 0 = s 1 = 0. The same argument can be used to show that S 21 ∩ S 22 = {0}.

Finally we need to show that S i1 + S i2 = S i . For any k 1 ∈ S 1 , k 1 = s 1 p + s 0 , where s 0 < p λ . Furthermore s 0 can be expressed as s 0 = s 01 p c1 + s 00 , for 0 ≤ s 00 < p c1 . Setting k 11 = s 1 p λ + s 01 p c1 ∈ S 11 and k 12 = s 00 ∈ S 12 , we have k 1 = k 11 + k 12 ∈ S 11 + S 12 . Conversely, if k 1 = k 11 + k 12 ∈ S 11 + S 12 , we have k 1 = s 1 p + s 0 + k 11 , such that s 0 = ap c1 < p λ and k 11 < p c1 . We have s 0 = s 0 + k 11 < p λ . As such k 1 ∈ S 1 . A similar argument can be used to show that S 21 + S 22 = S 2 .

We also have the following result, similar to [DM19] Lemma 1.

Lemma 10. There are at least r p + 1 representations of k as the sum

k 1 + k 2 with (k 1 , k 2 ) ∈ S 1 × S 2 . Proof. Let k ∈ N such that k (mod r) = k, 0 ≤ k < r. Let 0 ≤ k 0 < p
and a be the unique integers such that k = a p + k 0 , and let 0 ≤ c < p λ and b be the unique integers such that k 0 = b p λ + c . Let S 1 be the subset of S 1 such that S 1 = {k 1 = c + s 1 p}. We show that for any k 1 ∈ S 1 such that k 1 < k , there exist a k 2 ∈ S 2 such that k = k 1 + k 2 (mod r) = k (mod r).

For any

k 1 ∈ S 1 , k 1 ≤ k , there is some 0 ≤ s 1 ≤ a < r/p such that k 1 = s 1 p + c . Let 0 ≤ s 2 ≤ a be such that s 2 + s 1 = a , and let k 2 = b p λ + s 2 p. It is clear that k 2 ∈ S 2 , and k 1 + k 2 = c + b p λ + (s 1 + s 2 )p = k .
Similarly, if we consider k ∈ N such that k = k + r, we can find unique a , b , c such that k = a p + b p λ + c such that 0 ≤ c < p λ and 0 ≤ b p λ + c < p. We consider the sublist S 1 of S 1 of all k 1 such that k 1 = c + s 1 p. With the same argument, we can show that for all

k 1 = c + s 1 p λ ∈ S 1 such that k 1 > k , there is a k 2 such that k 2 = b p λ + s 2 p, with s 2 + s 1 = a such that k 1 + k 2 = k . Furthermore, since k 1 > k , k 2 = k + r -k 1 < r, which proves that k 2 is indeed in S 2 .
It follows that the number of representations (k 1 , k 2 ) ∈ S 1 × S 2 of k is at least equal to some N = N 1 + N 2 , where N 1 is the number of elements k 1 of S 1 such that k 1 < k and N 2 is the number of elements k 1 of S 1 such that k 1 > k . We have N 1 = a + 1 and N 2 ≈ r/p -a it follows that N ≈ r/p + 1, which concludes the proof.

Complexity. Using the results above, we are ready to prove Theorem 2.

Proof (of Theorem 2). For any δ > 0, we set = δ/(2c 2 ). With our definition of the S ij sets, we have

|L 11 | = p c1 , |L 12 | = O p 2c1 , |L 21 | = p c2 , |L 22 | = O p 2c2 .
Furthermore, under Heuristic 1, we have |L| = |L 11 ||L 12 |/p = p 3c1-1 . Combining this with the results of Lemma 8 and Lemma 7, we obtain that Algorithm 5 runs in time Õ p 2c2(1+ ) = Õ p 2c2+δ , using Õ max p 3c1-1 , p 2c2+δ , after a precomputation step which takes time Õ p 2c1+δ .

we have c 1 = 3 5 and c 2 = 2 5 , this gives us

3c 1 -1 = λ = 4 5 = 2c 2 .
In this case |L| = |L 2 |. For larger values of λ the size of L increases, and for smaller values of λ the size of L 2 increases, which makes λ = 4 5 the value for which the memory consumption of this algorithm is minimised.

Extending the Algorithm to Other Fields

Here we investigate an extension of this algorithm to fields F p , where ≥ 2 is a constant.

Let (1, α, . . . , α -1 ) be a basis of F p seen as an F p -vector space. We consider an elliptic curve E over F p defined by its Weierstraß's equation: y = f (x) = x 3 + ax + b, a subgroup G = P of prime order r = O p 2 and a target point

Q ∈ G.
Once again, the idea is to split the problem into four lists L 11 , L 12 , L 21 and L 22 such that

L 1j = {k 1j P : k 1j ∈ S 1j } L 2j = {-k 2j P : k 2j ∈ S 2j }.
For well chosen subsets S ij of N. These lists, which do not depend on Q can be pre-computed. We can also pre-compute the join L between L 11 and L 12 such that all points (x, y) ∈ L satisfy x = u 0 + αu 1 + • • • + α -2 u -2 . The main algorithm would then consist in replacing L 21 by L 21 + {Q}, compute the join L between L 21 and L 22 and search for a collision between L and L .

The pre-computation is given in Algorithm 6 and the main procedure in Algorithm 7.

In this section, we prove the following theorem.

Theorem 3. For any δ > 0, Algorithm 7 solves ECDLP, with base field F p in time Õ p c2+δ , using Õ max(p ( +1)c1-1 , p c2+δ ) memory, after a pre-computation which takes time Õ p 2 1 . The memory consumption of our algorithm is minimised when

c 1 = +1 2 +1
and c 2 = 2 +2 . In this case, our algorithm runs in time Õ p The idea is quite similar to the one that has been developed in the previous section. We simply give more details below.

Algorithm 6 Precomputation((a, b), P )

Require: An elliptic curve E over F p defined by the coefficients (a, b) of its Weierstraß's equation y = f (x), a generator P of G. Ensure: The lists L, L21, L22, and the polynomial f.

1: G ← ((X1-X2) 2 (X1+X2+X3)-f (X1)-f (X2)) 2 -4f (X1)f (X2) in F p [X1, X2, X3]
2: g0, g1, . . . , g -1 ← Weil-Descent(G) in Fp[U1,0, . . . U 1, -1 , U2,0, . . . , U 2, -1 , U3,0, . . . , U 3, -1 ] 3: for all 0 ≤ i < do 4: gi ← gi(U1,0, . . . , U 3, -2 , 0) in Fp[U1,0, . . . , U 3, -2 ] 5: f ← Elimination(g0, . . . , g -1 , U3,0, . . . , U 3, -2 )

6: L11 ← {k11P : k11 ∈ S11} 7: L12 ← {k12P : k12 ∈ S12} 8: L21 ← {-k21P : k21 ∈ S21} 9: L22 ← {-k22P : k22 ∈ S22} 10: L ← Join((a, b), L11, L12, f) 11: return L, L21, L22, f Algorithm 7 DM19-ECDLP((a, b), Q, L, L 21 , L 22 , f)
Require: An elliptic curve E over F p defined by the coefficients (a, b) of its Weierstraß's equation, a target point Q, the lists L, L21 and L22 the polynomial f associated to E.

Ensure: k ∈ Zr such that kP = Q or ⊥. 1: L21 ← L21 + {Q} 2: L ← Join((a, b), L21, L22, L, f) 3: for all Q -k2P ∈ L do 4: if Q -k2P = k1P ∈ L then 5: return k1 + k2

Analysis of the Algorithm

Computing the join. We generalise the result of lemma 6 to the case of an elliptic curve defined over F p .

Lemma 11. Let L 1 and L 2 be two lists of points of E(F p ), and let f ∈ F[U 1,0 , . . . , U 1, -1 , U 2,0 , . . . , U 2,1 ] be a constant degree polynomial. It is possible to compute the list L such that:

L = {P 1 + P 2 , P 1 = (x 1 , y 1 ) ∈ L 1 , P 2 = (x 2 , y 2 ) ∈ L 2 : f(u 1,0 , . . . u 2, -1 ) = 0 }, in time and space Õ |L 2 | 1+ provided that |L 1 | = |L 2 | = p c , where 0 < c < 1.
Once again, for simplicity we give the proof for 0 < c ≤ 1/2 and for 1/2 < c ≤ 3/4. We claim that the proof can be extended for any c such that 2 γ-1 -1 2 γ-1 < c ≤ 2 γ -1 2 γ , where γ is any chosen constant.

Proof. Let 0 < c ≤ 3/4.

The first step of the join computation is to build the lists A f and B, such that A f is the list of polynomials

f x1 in F p [U 2,0 , . . . , U 2, -1 ] with f x1 (U 2,0 , . . . , U 2, -1 ) = f(u 1,0 , . . . , u 1, -1 , U 2,0 , . . . , U 2, -1 ) x 1 = -1 i=0 u 1,i α i ,
for all (x 1 , y 1 ) ∈ L 1 , and

B = {(u 2,0 , . . . , u 2, -1 ) : x 2 = u 2,0 + • • • + α -1 u 2, -1 , (x 2 , y 2 ) ∈ L 2 }.
We claim that this can be performed in time and memory O (max

(|L 1 |, |L 2 |)) = O (|L 2 |).
The second step consists in computing the polynomial F = A f f (x1) . Using fast polynomial product this can be done in time and memory Õ |L 1 | .

In the third step we need to compute E = MultipointEvaluation(F, B). Using the Kedlaya-Umans Algorithm, for any > 0 this can be done in time and memory

Õ |L 1 | + |L 2 | 1+ = Õ |L 2 | 1+
. Afterward, we compute B the sublist of B of all points (u 2,0 , . . . u 2, -1 ) such that F(u 2,0 , . . . , u 2, -1 ) = 0. This can be done in time and memory O (|L 2 |).

Now, if 0 < c ≤ 1/2, the fourth and last step consists in finding all pairs (f x1 , (u 2,i ) -1 i=0 ) such that f x1 (u 2,0 , . . . , u 2, -1 ) = 0. This is performed by an exhaustive search which takes time Õ (|L 1 ||B |). We can estimate the size of B using Schwarzt-Zippel lemma. We have

|B | = O (|L 1 ||L 2 |/p) = O p ( +1)c-1 . It follows that |L 1 ||B | = O p ( +2)c-1 . And ( + 2)c -1 ≤ c when c ≤ 1/2. As such |L 1 ||B | If 1/2 < c ≤ 3/4
, we use a similar trick than in Section 4.1. We dispatch A f into p κ1 buckets A

(i) f such that |A (i) f | = p c-κ1
, where κ 1 = 1 -c. We compute all the E (i) lists such that E (i) = MultipointEvaluation(F i , B ) where

F i = A (i) f f x1 . Then we compute the list B (i) = {B [j] : E (i) = 0}.
(i) f | = p (c-κ1) = p (2c-1) , |B | = p ( +1)c-1 . It follows that |B | > |A (i) f | and thus the complexity of this step is Õ p κ1 |B | 1+ = Õ p c+(( +1)c-1) = Õ |L 2 | 1+ , for = (( + 1)c -1) /( c).
The last step consist in finding all pairs (f x1 , (u 2,i ) -1 i=0 ) such that f x1 (u 2,0 , . . . , u 2, -1 ) = 0. This is performed by an exhaustive search which takes time

Õ p κ1 |A (i) f ||B (i) | = Õ |L 1 ||B (i) | . We estimate |B (i) | with Schwartz- Zippel Lemma. |B (i) | = O p ( +3)c-3 . It follows that |L 1 ||B (i) | = O p ( +4)c-3 and ( +4)c-3 ≤ c when c ≤ 3/4. As such |L 1 ||B (i) | = O (|L 2 |), which conclude the proof for c ≤ 3/4.
This proof can be extended to 2 γ-1 -1

2 γ-1 < c ≤ 2 γ -1 2 γ
, where γ is any chosen constant, by repeating the same argument γ times, and choosing a sequence of parameter κ 0 , . . . , κ γ-1 such that κ i = i(1 -c).

Pre-computation step. We show the following result, regarding the complexity of the pre-computation step.

Lemma 12. For any constant > 0, the pre-computation step given in Algorithm 6 can be performed in O |L 12 | 1+ operations over F p , for a "good choice" of the S ij sets.

Proof. Similar to the proof of Lemma 7, we need three steps in this proof. The first step consists in showing that computing the polynomial f requires a constant number of operations over F p . Then we argue that building the lists L 11 , L 12 , L 21 and L 22 requires a time that is proportional to the size of the largest one, that is, following the same argument than in the proof of Lemma 7, O (|L 12 |) for a good choice of the sets S ij . The last step, finally, consists in showing that the join of L 11 and L 12 can be computed in time Õ |L 12 | 1+ . Assuming that |L 11 | = |L 12 | < p (this can be enforced by the choice of the S ij sets), this is a direct consequence of Lemma 11.

In the end, all we need to prove is that f can be computed in constant time. The argument is quite similar to the one given in the proof of Lemma 7. Since G is a polynomial of constant individual degree and is a constant as well, the polynomials g 0 , . . . , g -1 would be polynomial of constant degree individual degree as well, and would have a constant number of variables. Furthermore their computation can be done in constant time. Computing f would require to compute the elimination ideal of g 0 , . . . , g -1 in the variables {U 3,0 , . . . , U 3, -2 }. As long as the number of variables of each polynomials and the individual degree of these polynomial are all constants, this can be performed in constant time as well.

Main algorithm. We show the following result, regarding the complexity of the main algorithm. Choice of the S ij sets. We now discuss our choice of the S ij sets. Let -1 2 ≤ λ ≤ -1 be a parameter. We assume that p λ is an integer. If not the following stills holds by replacing p λ by its rounding to the closest integer. Let S 1 and S 2 be the following subsets of N.

S 1 = {k 1 < r : k 1 = s 1 p + s 0 , 0 ≤ s 0 < p λ } (5)

S 2 = {k 2 < r : k 2 = s 1 p + s 0 , 0 ≤ s 0 < p and p λ |s 0 }.

Similar to the choice that has been done in Section 4.3, the S ij sets are then built so that S i1 ∩ S i2 = {0}, S i1 + S i2 = S i and |S i1 | = |S i2 |. For instance, let c 1 = λ-1 -1 and let c 2 = -λ +1 . We have

0 < c 1 ≤ λ ≤ -1 ≤ c 2 + λ < .
Following the same kind of argument than the ones given in the proof of Lemma 9 reader can check that if we choose the S ij sets to be S 11 = {k 11 ∈ S 1 : k 11 < p c1 } S 12 = {k 12 ∈ S 1 : k 12 = s 1 p -1 + s 0 , p c1 |s 0 } S 21 = {k 21 ∈ S 2 : k 21 = sp -1 , p (c2-1)+λ+1 |s} S 22 = {k 22 ∈ S 2 : k 22 = s 1 p -1 + s 0 , s 1 < p (c2-1)+λ+1 , p λ |s 0 }, these conditions are satisfied. The shape of the elements of the S ij sets defined above is illustrated in Figure 4b.

Lemma 14. There are at least r p + 1 representations of k as the sum k 1 + k 2 with (k 1 , k 2 ) ∈ S 1 × S 2 .

The proof of this lemma is very similar to the one of Lemma 10, as such we do not detail it here.

Proof of Theorem 3. Combining all these results, we are ready to prove Theorem 3. Combining Lemma 13 and Lemma 12, we obtain that for any constant δ > 0, the time complexity of the main algorithm is Õ p c2+δ and the memory required is Õ max(p c2+δ , p ( +1)c1-1 ) after a pre-computation which requires Õ p c1+δ .

The memory consumption is reduced when c 2 = ( + 1)c 1 -1, Our choice of c 1 = +1 2 +1 and c 2 = 2 +2 , is derived from this equation and from the fact that c 1 + c 2 = 1.

Discussion

We now address some open questions and remarks regarding these algorithm and their possible improvements.

Lemma 12 to compute the join, we would need to consider polynomials f 1 , . . . , f , in O ( ) variables and compute a Gröbner basis to discard the dependency in the last variable x 3 . This would already cost more than we can afford during the pre-computation step. It could also be nice to extend the method to a field F p where p is a prime.

Lemma 5 .

 5 Fig. 2: Bits repartition of elements from S 11 , S 12 = S 22 , S 21 as vector of F log r 2 , the leftmost part being the least significant bits of the numbers. The white rectangles represent a portion of zeroes.

  for some constant c. Step 7 can be performed in time O (|L 2 |). Finally the last part of the algorithm consists in an exhaustive search which takes time Õ (|L 1 ||L 2 |), where L 2 is the sublist of L 2 of the elements which pass the multipoint evaluation test. Similar to [DM19], we rely on Schwartz-Zippel lemma and Heuristic 1 to argue that the expected value of |B | is E[|B |] ≤ d|L 2 | p = O p 3c/2-1 . Since c ≤ 1, we have 3c/2 -1 ≤ c/2 and thus |B | = O (|L 1 |). It follows that this last step can be performed in time O |L 1 | 2 . Taking into account that |L 1 | 2 = |L 2 | and combining all intermediary results, we get the claimed time complexity.

Theorem 1 .

 1 Using the algorithm from [KU08] as the MultipointEvaluation routine, for any constant > 0 Algorithm 1 solves p 2 -ECDLP in time and memory Õ p 1+ . Proof. In [DM19] Delaplace and May chose the sets S ij for i, j ∈ {1, 2} so that |L 11 | ≈ |L 21 | = O √ p and |L 12 | = |L 22 | = O (p), (see Figure 1). More precisely, they set S 11 = {k 11 ∈ N : k 11 < √ p } S 21 = {k 21 ∈ N : k 21 < p, ∃s ∈ N : k 21 = √ p s} S 12 = S 22 = {k 12 ∈ N : k 12 < r, ∃s ∈ N : k 12 = ps}, where √ p is the closest integer from √ p (see Figure 2 to see how the bits of these elements are dispatched). As direct consequence of Lemmas 4 and 5, for any constant > 0, Algorithm 1 solves p 2 -ECDLP in time and memory Õ max |L 12 | 1+δ , |L 22 | 1+ = Õ p 1+ . 4 A New p 2 -ECDLP Algorithm with Pre-computation In the section, we present a new algorithm for the p 2 -ECDLP problem, based on the work of [DM19].

  already argued in [DM19], building the lists L 12 , L 12 , L 21 , L 22 can be done in O (max i (|L i |)) operations over F p . In order to minimise the runtime of the join routine, we choose |L 11 | 2 = |L 12 | and |L 21 | 2 = |L 22 |. We would also like to Algorithm 5 DM19-ECDLP((a, b), Q, L, L 21 , L 22 , f) Require: An elliptic curve E over F p 2 defined by the coefficients (a, b) of its Weierstraß's equation, a target point Q, the lists L, L21 and L22 the polynomial f associated to E. Ensure: k ∈ Zr such that kP = Q or ⊥.

Lemma 8 .

 8 Given L, L 21 , L 22 as described in the previous section such that |L 21 | < p, for any fixed constant > 0, Algorithm 5 runs in time Õ (|L 22 |) and uses Õ max(|L|, |L 22 | 1+ ) elements of F p as storage.

  λ) log p log r -log p (a) Bits repartition of elements from S1 and S2, as vector of F log r 2 , the leftmost part being the least significant bits of the numbers. The white rectangles represent a portion of zeroes. λ) log p log r -(λ + c2) log p c2 log p λ log p (b) Bits repartition of elements from S11, S12, S21 and S22 as vector of F log r 2 , the leftmost part being the least significant bits of the numbers. The white rectangles represent a portion of zeroes.

Fig. 3 :

 3 Fig. 3: Bits repartition of the elements of S i and S ij for i ∈ {1, 2}, j ∈ {1, 2}.

  memory after a pre-computation step having runtime Õ p

  For any constant this step can be performed in time and memory O p κ1 |A (i) f | + |B | 1+ . Furthermore, we have |A

Lemma 13 .

 13 For any constant > 0, Algorithm 7 solves ECDLP for an elliptic curve defined over a finite field F p , in time Õ |L 21 | 1+ and memory Õ max |L|, |L 21 | 1+ . Proof. Replacing the points P 1 of L 21 by P 1 + Q can be done in O (|L 21 |) operations over F p . We have already argued that if |L 21 | = |L 22 | < p , for any constant > 0, the join between L 21 and L 22 can be performed in time Õ |L 22 | 1+ . From here, the proof is trivial.

Proof.

  Let the S ij sets be defined as above. This gives us |L 11 | = p c1 , |L 12 | = O p c1 , |L 21 | = O (p c2 ) and |L 22 | = p c2 . it follows that |L| = p ( +1)c1-1 .

  Consider that k 1 is the sum k 11 + k 12 for a unique pair (k 11 , k 12 ) ∈ S 11 × S 12 where S 11 and S 12 are well chosen subsets of S 1 such that S 11 + S 12 = S 1 . Accordingly, k 2 is split uniquely as k 2 = k 21 + k 22 , with (k 21 , k 22 ) ∈ S 21 × S 22 , where S 21 and S 22 are two subsets of S 2 satisfying S 21 + S 22 = S 2 . The algorithm then works by building base lists L 11 , L 12 , L 21 and L 22 such that L 11 = {k 11 P : k 12 ∈ S 11 } L 12 = {k 12

  From the Schwartz-Zippel Lemma, under Heuristic 1, we can expect |L 12 | to be O |L 12 | 3/2 /p , assuming that |L 11 | 2 = |L 12 |. Now, if |L 12 | ≥ p, |L 12 | ≥ √ p, and thus |L 11 ||L 12 | ≥ |L 12 |.

  then 5: return k1 + k2 compute the join between L 21 and L 22 faster than the one between L 11 and L 12 , as this join would be part of the actual algorithm. So it would be better if the lists L 21 and L 22 are smaller than L 11 and L 12 . It follows that L 12 is the largest list to be computed, and thus the runtime of this step is O (|L 12 |) elementary operations over F p .

Funded by DFG under Germany's Excellence Strategy -EXC 2092 CASA -390781972.

Reducing the memory consumption? As already mentioned, the main bottleneck of our algorithm is the large memory required to store L as soon as the parameter c 1 gets larger than +1 2 +1 . Note that this gets worse when increases. If we could find a way to reduce the storage of L, we should be able to extend the results given in Theorem 3 to a larger set of parameters, and get closer to Bernstein and Lange p /3 complexity, or possibly reach it for > 2. However doing this without increasing the time complexity of the algorithm nor the pre-computation does not appear to be an easy task.

Extending the algorithm to other fields? Another interesting open question is whether our algorithm can be extended to other fields. Here, we focused on fields F p for (small) constant and large (prime) p. Another cryptographically relevant case is the one where is a large prime and p a small prime (e.g. p = 2). It would be possible to extend the algorithms described here to this case with some changes to make it fit the model. However, similar to what is shown in