L R García Carrillo

A Dzul

R Lozano

Hovering quad-rotor control: A comparison of nonlinear controllers using visual feedback

Keywords: Autonomous helicopter, nested saturations, backstepping, sliding modes, visual feedback, position stabilization

This article presents a comparison of three control techniques: Nested Saturations, Backstepping and Sliding Modes. The control objective consists of obtaining the best control strategy to stabilize the position of a quad-rotor unmanned aerial vehicle (UAV) when using visual feedback. We propose a visionbased method to measure translational speed as well as the UAV 3D position in a local frame. The three selected controllers were implemented and tested in real-time experiments. The obtained results demonstrate the performance of such methodologies applied to the quad-rotor system.

I. INTRODUCTION

With a wide range of applications in both civilian and military scenarios, automatic control of flying machines has attracted the attention of many researches in recent years. Operations of search and rescue, inspection and sensing of remote areas, hazardous material recovery, real-time forest fire monitoring, disaster relief support, surveillance of sensitive areas (borders, ports, oil pipelines), etc., represent some of the tasks where UAV's have demonstrated to be very powerful tools. UAV's can fly autonomously or semi-autonomously and, in addition, they are expendable or recoverable [START_REF] Reinhardt | Future Employment of UAVS: Issues of Jointness[END_REF].

Several earlier works have already proven the real-time embedded stabilization of a small quad-rotor [START_REF] Salazar | Embedded control system for a four rotor UAV[END_REF]. In addition, testing the performance of different controls over such systems is a subject already studied. PID and LQR controllers L.R. García Carrillo is with Heudiasyc, UMR CNRS 6599, Université de Technologie de Compiègne, Compiègne, FRANCE. e-mail: (lgarciac@hds.utc.fr).

A. Dzul is with Instituto Tecnologico de la Laguna, Torreon, Coahuila, MEXICO. e-mail: (dzul@faraday.itlalaguna.edu.mx). R. Lozano is with Heudiasyc, UMR CNRS 6599, FRANCE and LAFMIA UMI 3175, Cinvestav, MEXICO. e-mail: (rlozano@hds.utc.fr).

Manuscript received: April 2011.

are compared in [START_REF] Bouabdallah | PID vs LQ control techniques applied to an indoor micro quadrotor[END_REF], while in [START_REF] Bouabdallah | Backstepping and Sliding-mode Techniques Applied to an Indoor Micro Quadrotor[END_REF], the performance of a backstepping and a sliding modes controllers are tested. For the last two examples, experiments were performed over a quad-rotor platform, where 3 degrees of freedom are locked. Authors conclude that backstepping control technique is the most appropriate approach for their future works. In [START_REF] Altug | Control of a Quadrotor Helicopter Using Dual CameraVisual Feedback[END_REF], two control methods are studied over a quad-rotor platform equipped with visual feedback. These methods are based on feedback linearization and a backstepping-like control.

Vision systems have become a popular choice for obtaining information that can be used in the feedback control loop of autonomous vehicles. This data is usually combined with an Inertial Measurement Unit (IMU) to provide robust relative attitude information and allowing autonomous positioning and navigation. Also, note that an on-board vision system increases the performance of an UAV equipped with a global positioning system (GPS), which provides position information relative to an inertial frame, but fails in indoor or in noisy environments. On-board computer vision systems provide information which is obtained, for example, from the detection of landmarks. This information allows the UAV to estimate its position in a local frame. Once the UAV knows its position, a control strategy could be implemented in order to achieve a desired position. Several research works aiming at the control of UAV's, using a camera as a vision sensor, can be found in [START_REF] Saripalli | Visually-Guided Landing of an Unmanned Aerial Vehicle[END_REF], [START_REF] Salazar | Modeling and Real-Time Stabilization of an Aircraft Having Eight Rotors[END_REF], [START_REF] Rondon | Vision-based altitude, position and speed regulation of a quad-rotor rotorcraft[END_REF], [START_REF] Salazar | Vision Based Autonomous Hover of a Mini-rotorcraft[END_REF], [START_REF] Gomez | Vision Based Tracking for a Mini-Rotorcraft Using Vanishing Points[END_REF]. A vision algorithm for visual navigation and landing of a gas-powered radio-controlled model helicopter, equipped with a PC-104 is presented in [START_REF] Saripalli | Visually-Guided Landing of an Unmanned Aerial Vehicle[END_REF]. Concerning electric-powered vehicles, an approach based on optical flow techniques has been applied to the real-time stabilization of an Eight-Rotor UAV in [START_REF] Romero | Real-Time Stabilization of an Eight-Rotor UAV Using Optical Flow[END_REF]. Also, a two cameras system is presented in [START_REF] Rondon | Vision-Based Position Control of a Two-rotor VTOL miniUAV[END_REF], with the drawback that the cameras are not embedded on the UAV. A stereo vision system combined with inertial measurements and a laser range finder is proposed in [START_REF] Achtelik | Stereo vision and laser odometry for autonomous helicopters in GPS-denied indoor environments[END_REF], with the purpose of enabling relative localization and local navigation.

In this article, we are interested on experimentally evaluating the performance of three different control strategies for stabilizing a quad-rotor UAV. The aircraft, shown in Figure 1, is equipped with an embedded IMU and gyros, providing Euler angles and angular rates respectively. In addition, we have developed an algorithm for landmark detection and tracking, which estimates the UAV motion (relative position and linear velocity) with respect to a landing pad on the ground, using data provided by an onboard monocular camera [START_REF] Garcia Carrillo | Stabilization and trajectory tracking of a quad-rotor UAV using vision[END_REF].

The control objective of our experiments consists of validating the most effective controller for stabilizing the vehicle's position with respect to an artificial visual landmark on the ground. For our tests, we have chosen three controllers considered between the most commonly reported in the literature: nested saturations [START_REF] Castillo | Real-Time Stabilization and Tracking of a Four-Rotor Mini Rotorcraft[END_REF], [START_REF] Teel | Global stabilization and restricted tracking for multiple integrators with bounded controls[END_REF], backstepping [START_REF] Saberi | Global Stabilization of Partially Linear Composite Systems[END_REF], and sliding modes [START_REF] Khalil | Nonlinear Systems[END_REF]. Experimental results have shown that the nested saturation controllers offer a smoother UAV behavior than the other two controllers, which leads to less energy consumption.

The paper is organized as follows. Section II presents the dynamical equations of the quad-rotor UAV. In Section III we give details of the visual system setup. The control strategies for hovering flight are described in Section IV. The complete quad-rotor experimental platform is presented in Section V. In order to evaluate the effectiveness of the implemented controllers, experimental results are shown in Section VI. Some final concluding remarks are presented in Section VII.

II. QUAD-ROTOR DYNAMICAL MODEL

A quad-rotor helicopter is an underactuated dynamic vehicle composed of four input forces (the thrust provided by each propeller) and six output coordinates (fully spatial movements). It could be considered as a vertical takeoff and landing vehicle (VTOL), able to move omnidirectionally and with the ability to fly in hover. The quad-rotor behavior is controlled by varying the angular speed of the rotors. Each rotor produces a thrust and a torque, whose combination generates the main thrust, the yaw torque, the pitch torque, and the roll torque acting on the rotorcraft. In the quad-rotor, the front and rear rotors rotate counter-clockwise while the left and right rotors rotate clockwise, canceling gyroscopic effects and aerodynamic torques in stationary trimmed flight. Vertical motion is controlled by the collective throttle input, that is, the sum of the thrusts of each motor. As we can see in Figure 2, forward/backward motion is achieved by controlling the differential speed of the front and rear motors. This causes the quad-rotor to tilt around the corresponding axis, generating a pitch angle. The left/right motion of the vehicle is achieved by controlling the differential speed of the right and left motors, tilting around the corresponding axis and producing a roll angle. Finally, yaw movement is obtained by taking advantage of the two sets of rotors rotating in opposite direction. Thus, a yaw angular displacement is obtained by increasing (or decreasing) the speed of the front and rear motors while decreasing (or increasing) the speed of the lateral motors. This is done keeping the total thrust constant, then the altitude remains unchanged.

The quad-rotor representation used for this paper is shown in Figure 2. The dynamic model of this aircraft is obtained by representing the quad-rotor as a solid body evolving in 3D and subject to one force and 3 moments [START_REF] Etkin | Dynamics of Flight[END_REF] (2)

mz = u cos θ cos φ -mg (3) ψ = τψ (4) θ = τθ (5) φ = τφ (6)
where u is the main thrust directed out of the top of the aircraft, m is the mass of the quad-rotor, g denotes the gravity constant, x and y are coordinates in the horizontal plane, z is the vertical position, and τψ , τθ , and τφ are the yawing moment, pitching moment and rolling moment respectively, which are related to the generalized torques τ ψ , τ θ , τ φ .

III. VISUAL SYSTEM SETUP

Controlling the 3-dimensional position of an UAV depends on the knowledge of the (x, y, z) vehicle coordinates and ( ẋ, ẏ, ż) translational velocities with respect to a well-known reference frame. Such values are required data for the controller in order to generate the control inputs to stabilize the aircraft over a desired location. In order to fulfill this situation, a vision system can be implemented to provide the required position and velocities information. The vision system proposed consists of a calibrated camera onboard the UAV, a landing pad or artificial marker placed on ground, a vision algorithm running on a supervisory ground station PC, and a wireless link between the helicopter and the supervisory ground station. Figure 2 shows the proposed system which can be described as:

• quadrotor UAV: with a body fixed frame

(X h , Y h , Z h ),
assumed to be at its center of gravity. Z h represents the yaw axis, and pointing upwards. X h and Y h are the roll and pitch axis respectively.

• strapdown camera: pointing downwards, with a reference frame (X c , Y c , Z c ). When moving, the camera surveys the scene passing below the quadrotor. Since X c -Y c and X h -Y h are considered as parallel planes, then visual information collected by the camera can be used to stabilize the vehicle.

• landing pad: artificial landmark of known dimensions, formed by four circles of known coordinates, painted on high contrast background and placed underneath the rotorcraft. The coordinates frame (X lp , Y lp , Z lp ) represents the inertial reference frame. The planes formed by (X h -Y h ) and (X lp -Y lp ) are considered to be parallel because we assume that the rotorcraft is in hover flight over the landing pad.

A. Vision-based position estimation

In order to estimate the UAV position relative to the landing pad, the extrinsic parameters of the camera are computed at every image frame. This is achieved by implementing an homography estimation technique, which provides the (x, y, z) position and (ψ, θ, φ) orientation of the camera with respect to the artificial landmark in the image scene. The action of the homography can be expressed as [START_REF] Bradski | Learning OpenCV, Computer Vision with the OpenCV Library[END_REF]:

  x y 1   = sM r 1 r 2 r 3 t     X Y Z 1     (7)
where [x y 1] T represents the landing pad position in the camera image, s is a known scale factor, M ∈ R 3×3 represents the intrinsics parameters camera matrix, R = [r 1 r 2 r 3 ] ∈ R 3×3 are the extrinsics rotation parameters, t ∈ R 3×1 is the extrinsics translation parameters vector , [X Y Z 1] T is the real landing pad position, and H = sM [r 1 r 2 t] is the homography matrix. The homography matrix H is divided in two parts: the physical transformation (which locates the observed object plane) and the projection (the camera intrinsic matrix).

Rotation R is described by three angles and translation t is defined by three offsets; hence there are six unknowns for each view. The known planar object (the artificial landmark) provides eight equations, that is, the mapping of a rectangle into a quadrilateral can be described by four (x, y) image points. For every instant, when the aerial vehicle is in hovering, it is possible to compute the H homography matrix using the a priori knowledge of the position of the four centroids of the circles [START_REF] Horaud | An Analytic Solution for the Perspective 4-Point Problem[END_REF]. Using this estimated transformation matrix and the intrinsic camera matrix previously identified by an off-line calibration based on the method in [START_REF] Bouguet | Camera Calibration Toolbox for Matlab[END_REF], we are able to calculate the camera extrinsic parameters, and therefore we have the vehicle's (x, y, z) position with respect to the landing pad on the ground.

B. Translational velocities

An optical flow computation procedure is applied to compute the ( ẋ, ẏ, ż) translational velocities of the aerial vehicle with respect to the landing pad. In order to compute optical flow, we have implemented the Lucas-Kanade pyramidal algorithm [START_REF] Bouguet | Pyramidal Implementation of the Lucas Kanade Feature Tracker -Description of the algorithm[END_REF] in combination with a feature-detecting algorithm. This approach provides an accurate estimation of the motion field since it does not take into account the non landing pad areas, where the motion field cannot be accurately determined.

Consider the camera moving with respect to a rigid scene. The velocities and rotation rates of the camera in the inertial frame are expressed by (V x , V y , V z ) and (w x , w y , w z ) respectively. To accurately estimate the pseudo-speeds of the engine, we define a tracking zone surrounding the landing pad, in a way that the centroid of the zone and the center of the landing pad coincide. The most representative features over the zone are selected as features to track for. These features are usually the circles perimeter. Once this group of features has been identified, a tracking process is performed over the entire image. The optical flow is estimated based on the displacements of the tracked features.

The optical flow computed at point (x i k , y i k ) is composed of a translational and a rotational part as

OF x i k OF y i k = T OF k + R OF k (8)
with the translational part

T OF k = 1 z -f 0 x i k 0 -f y i k   V x V y V z   (9)
and the rotational part

R OF k = x i k y i k f -(f + (x i k ) 2 f ) y i k (f + (y i k ) 2 f ) - x i k y i k f -x i k   ω x ω y ω z   (10)
where OF d j i is the optical flow component in the coordinate j of the point p i , V k and ω k are the translation velocities and rotation rates, respectively, of the body in the coordinate k. Thus, the mean of the optical flow computed on all those points can be expressed as a function of the camera movement as follows

ŌF x = VOFx + K x VOFz + ROFx (11) ŌF y = VOFy + K y VOFz v + ROFy (12)
Using the results from [START_REF] Rondon | Optical flow-based controller for reactive and relative navigation dedicated to a four rotor rotorcraft[END_REF], the rotational optical flow is compensated and the pseudo-speeds ( VOFx , VOFy , VOFz ) are deduced. Since the camera system and the helicopter share the same movements, it can be said that the deduced pseudovelocities depend of the rotorcraft movement. Indeed, the camera is mounted onboard the quad-rotor and fixed in a way it has no freedom degree. Thus, it can be written that

VOFx = - f ẋ z ( 13
)
VOFy = - f ẏ z (14) VOFz = ż z (15)
where ( ẋ, ẏ, ż) is the speed vector of the rotorcraft center of gravity and z is the altitude. Thus, from these three equations the proposed optical flow vision system allows speed estimation of the rotorcraft up to a scale factor, when flying at constant altitude. Those estimations can be used to control the translational velocities of the rotorcraft.

If the landing pad is not successfully detected in the current image, the vision algorithm will fail. To overcome this situation, optical flow values can be used to estimate the position of the four circles centroids as

ρ k+1 x = ρ k x + ∆ T VOFx (16) ρ k+1 y = ρ k y + ∆ T VOFy (17)
where (ρ k x , ρ k y ) represents the circle's centroid position and ∆ T is the working frequency of the algorithm. This estimated centroids positions are used to estimate the homography at each time the detection of the landing pad fails.

IV. CONTROL STRATEGY

In this section, we describe three control strategies applied to the quad-rotor: nested saturations control method, backstepping approach, and sliding modes controller. All control laws are designed to stabilize the x, y, θ and φ states. The altitude z and the yaw angle ψ are stabilized by PD controllers.

A. Altitude and yaw control

The control of the vertical position (3) can be obtained by using the following control input [START_REF] Castillo | Real-Time Stabilization and Tracking of a Four-Rotor Mini Rotorcraft[END_REF]:

u = (r 1 + mg) 1 cos θ cos φ (18)
where

r 1 = -k vz ż -k pz e z (19)
with e z = z d -z as the z error position and z d as the desired altitude. k pz and k vz are positive constants. Thus, for the altitude dynamics, r 1 is a PD controller. In the case of the yaw angular position (4), we can apply

τψ = -k vψ ψ -k pψ e ψ ( 20
)
where e ψ = ψ d -ψ denotes the yaw error, ψ d represents the desired yaw angle, k pψ and k vψ denote the positive constants of a PD controller. Indeed, introducing equations ( 18) and ( 20) into set of equations ( 1)-( 4) and provided that cos θ cos φ = 0, we obtain

mẍ = (r 1 + mg) sin ψ tan φ cos θ + cos ψ tan θ (21) mÿ = (r 1 + mg) sin ψ tan θ - cos ψ tan φ cos θ (22) mz = -k vz ż -k pz e z ( 23
) ψ = -k vψ ψ -k pψ e ψ (24)
The control parameters k pψ , k vψ , k pz and k vz should be carefully chosen to ensure a stable well-damped response in the vertical and yaw axes [START_REF] Castillo | Real-Time Stabilization and Tracking of a Four-Rotor Mini Rotorcraft[END_REF]. From [START_REF] Bouguet | Camera Calibration Toolbox for Matlab[END_REF] and [START_REF] Bouguet | Pyramidal Implementation of the Lucas Kanade Feature Tracker -Description of the algorithm[END_REF] it follows that ψ → ψ d and z → z d .

B. Nested Saturations Control

Consider a system given by four integrators in cascade:

ẋ1 = αx 2 , ẋ2 = βx 3 , ẋ3 = γx 4 , ẋ4 = u (25)
where α, β, γ = 0 are constants. A nested saturations control input can be defined as [START_REF] Castillo | Real-Time Stabilization and Tracking of a Four-Rotor Mini Rotorcraft[END_REF], [START_REF] Teel | Global stabilization and restricted tracking for multiple integrators with bounded controls[END_REF]:

u = -σ b4 (k 4 z 4 + σ b3 (k 3 z 3 + σ b2 (k 2 z 2 + σ b1 (k 1 z 1 )))) ( 26
)
where z i , for i = 1 . . . 4, denotes a change of variables. k i > 0 are constants, and σ bi represent saturation functions defined as:

σ bi (s) =    -b i ; s < -b i s ; |s| ≤ b i b i ; s > b i (27)
where b i > 0 are constants denoting the bounds of the saturation functions. The z i 's are given by:

z 1 = x 4 + k 4 + k 3 + k 2 γ x 3 + k 3 k 4 + k 2 k 3 + k 2 k 4 βγ x 2 + k 2 k 3 k 4 αβγ x 1 (28)
z 2 = x 4 + k 4 + k 3 γ x 3 + k 3 k 4 βγ x 2 (29)
z 3 = x 4 + k 4 γ x 3 (30)
z 4 = x 4 (31)
Note that from ( 19) and ( 23) r 1 → 0. For a time T large enough, e z and e ψ are arbitrarily small, therefore, equations ( 21) and ( 22) reduce to

ẍ = g tan θ (32) ÿ = -g tan φ cos θ (33)
1) Control of the forward position and pitch angle: Consider the subsystem given by ( 5) and (32). Implementing a nonlinear control based on nested saturations allows in the limit a guarantee of arbitrary bounds for x, ẋ, θ and θ. To further simplify the analysis, we will impose a very small upper bound on |θ| in such a way that the difference tan(θ)-θ is arbitrarily small (θ ≈ tan(θ)). Therefore, the subsystem ( 5) and (32) becomes the following linearized system:

ẍ = gθ (34) θ = τθ (35)
which represents four integrators in cascade. Then, by using ( 25)-( 26) the controller is given by

τθ = -σ b4 θ + σ b3 θ + θ + σ b2 θ + 2θ + ẋ g +σ b1 θ + 3θ + 3 ẋ g + x g (36)
It is proved in [START_REF] Castillo | Real-Time Stabilization and Tracking of a Four-Rotor Mini Rotorcraft[END_REF] that θ, θ, x and ẋ converge to zero. To regulate x around a desired position, we rewrite the most inner term (associated to σ b1 ) as θ + 3 θ + 3 ẋ g + ex g , where e x is the position error, expressed as e x = x d -x. Here, x d represents the desired position reference for x.

2) Control of lateral position and roll angle: Consider the subsystem given by equations ( 6) and (33). Imposing a very small upper bound on |φ| in such a way that the difference tan(φ) -φ is arbitrarily small (φ ≈ tan(φ)), the subsystem ( 6) and (33) becomes

ÿ = -gφ (37) φ = τφ (38)
Using a similar procedure to the one proposed for the pitch control, we obtain

τφ = -σ b4 φ + σ b3 φ + φ + σ b2 φ + 2φ - ẏ g +σ b1 φ + 3φ -3 ẏ g - y g (39)
In order to regulate y around a desired position, we rewrite the most inner term (associated to σ b1 ) as φ + 3 φ -

3 ẏ g - ey g
, where e y is the position error, expressed as e y = y d -y. Here, y d represents the desired position reference for y.

C. Backstepping Control

The Backstepping technique provides a systematic method to obtain a control law from a chain of integrators. This methodology was introduced in [START_REF] Saberi | Global Stabilization of Partially Linear Composite Systems[END_REF].

1) Control of the forward position and pitch angle: Rewrite subsystem given by ( 34)-(35) as

ẋ1 = x 2 (40) ẋ2 = gθ 1 (41) θ1 = θ 2 (42) θ2 = τBθ (43)
where τBθ will define the final backstepping control input. In order to obtain this control input, consider each equation as a new subsystem, where the next state is taken as the input and it is defined as a virtual control to stabilize such a subsystem. Then, for our case, we begin with ẋ1 = x 2 (44)

ζ 1 = x 1 (45)
where x 2 represents the input, and ζ 1 the output. Let us propose a positive definite function V 1 = 1 2 x 2 1 , whose time derivative is given by

V1 = x 1 ẋ1 = x 1 x 2 (46)
and consider a virtual input

α 1 = (x 2 ) v = -k 1 x 1 ,
where k 1 is a positive constant. Then V1 = -k 1 x 2 1 . Now, let ζ 2 be the new output:

ζ 2 = x 2 -α 1 ( 47
)
The new subsystem that we are trying to stabilize is written as

ẋ1 = ζ 2 + α 1 (48) ζ2 = gθ 1 -α1 (49)
and let us propose a positive definite function

V 2 = V 1 + ζ 2 2 2 , then V2 = V1 + ζ 2 ζ2 = V1 + ζ 2 (gθ 1 -α1 ) (50) Define a virtual input α 2 = (gθ 1 ) = -k 2 ζ 2 + α1 , where k 2 is a positive constant. Then, V2 = V1 -k 2 ζ 2 2
. Now, let ζ 3 be a new output:

ζ 3 = gθ 1 -α 2 (51)
The new subsystem to be stabilized is written as

ẋ1 = ζ 2 + α 1 (52) ζ2 = gθ 1 -α1 (53) ζ3 = gθ 2 -α2 (54)
and let us propose the positive definite function

V 3 = V 2 + ζ 2 3 2 , then V3 = V2 + ζ 3 ζ3 = V2 + ζ 3 (gθ 2 -α2 ) (55) Define a virtual input α 3 = (gθ 2 ) = -k 3 ζ 3 + α2 , where k 3 is a positive constant. Then, V3 = V2 -k 3 ζ 2 3
. Let ζ 4 be the new output:

ζ 4 = gθ 2 -α 3 ( 56
)
The new subsystem to be stabilized is written as

ẋ1 = ζ 2 + α 1 (57) ζ2 = gθ 1 -α1 (58) ζ3 = gθ 2 -α2 (59) ζ4 = g θ2 -α3 (60)
and let us propose the Lyapunov candidate function

V 4 = V 3 + ζ 2 4 2 , then V4 = V3 + ζ 4 ζ4 = V3 + ζ 4 (gτ Bθ -α3 ) (61)
Let us propose the backstepping control input τBθ as

τBθ = 1 g (-k 4 ζ 4 + α3 ) (62)
where k 4 is a positive constant. Then

V4 = V3 -k 4 ζ 2 4 = -k 1 x 2 1 -k 2 ζ 2 2 -k 3 ζ 2 3 -k 4 ζ 2 4 (63)
With the proposed backstepping control input τBθ we have V4 < 0, and then the system (40)-( 43) is globally asymptotically stable. In order to express τBθ as a function of x 1 , x 2 , θ 1 and θ 2 , we need to rewrite ζ 4 and α3 as a function of such variables:

ζ 4 = gθ 2 + (k 1 + k 2 + k 3 )gθ 1 + (k 1 k 2 + k 1 k 3 + k 2 k 3 )x 2 +k 1 k 2 k 3 x 1 (64) α3 = -(k 1 + k 2 + k 3 )gθ 2 -(k 1 k 2 + k 1 k 3 + k 2 k 3 )gθ 1 -k 1 k 2 k 3 x 2 (65)
The final control input for the forward position and pitch angle is rewritten as

τBθ = - k 1 g x 1 - k 2 g x 2 -k 3 θ 1 -k 4 θ 2 (66)
where

k 1 = k 1 k 2 k 3 k 4 (67) k 2 = k 1 k 2 k 3 + k 1 k 2 k 4 + k 1 k 3 k 4 + k 2 k 3 k 4 (68) k 3 = k 1 k 2 + k 1 k 3 + k 1 k 4 + k 2 k 3 + k 2 k 4 + k 3 k 4 (69) k 4 = k 1 + k 2 + k 3 + k 4 (70)
2) Control of lateral position and roll angle: Rewrite the subsystem given by ( 37)-(38) as ẏ1 = y 2 , ẏ2 = -gφ 1 , φ1 = φ 2 and φ2 = τBφ . Using a similar procedure to the one proposed for the pitch control, the backstepping roll control input can be obtained as

τBφ = k 1 g y 1 + k 2 g y 2 -k 3 φ 1 -k 4 φ 2 (71)

D. Sliding Modes Control

Consider the system (25) of integrators in cascade, and rewrite this system in the form [START_REF] Khalil | Nonlinear Systems[END_REF]:

φ ε = f a (ϕ, ε) f b (ϕ, ε) + G(x)E(x)u + δ(t, x, u) (72) then, we have that ϕ = [x 1 x 2 x 3 ] T and ε = x 4 . Also f a (ϕ, ε) = αx 2 βx 3 γx 4 T and f b (ϕ, ε) = 0, G(x) = E(x) = 1, δ(t, x, u) = 0. Consider x 4 = -c 0 x 1 -c 1 x 2 - c 2 x 3 = φ(ϕ)
in order to stabilize the origin. The partial derivative of φ, with respect to ϕ, is given by

∂φ ∂ϕ = [-c 0 -c 1 -c 2 ] (73) then ∂φ ∂ϕ f a = -c 0 αx 2 -c 1 βx 3 -c 2 γx 4 (74)
Define a sliding surface

s = ε -φ(ϕ) = x 4 + c 0 x 1 + c 1 x 2 + c 2 x 3 = 0 (75)
whose time derivative is given by

ṡ = - ∂φ ∂ϕ f a (ϕ, ε) + u (76)
then, the control input can be expressed as

u = -c 0 αx 2 -c 1 βx 3 -c 2 γx 4 + v (77)
where ṡ = v and v = -ν(x) tanh s ε . By definition, we have that

ṡ = g(x)v + ∆(t, x, u) (78)
thus g(x) = 1 and ∆(t, x, u) = 0. Also, we need to satisfy

∆(t, x, u) g(x) ≤ ρ(x) + k 0 v (79) then 0 1 ≤ B ( 80
)
where B is a positive constant, thus ρ(x) = B and k 0 = 0. We have then

ν(x) ≥ ρ(x) 1 -k 0 ≥ B (81)
Now we can determine that v = -ν(x) tanh s ε or v = -B tanh s ε . Note that tanh (s/ε) is a smooth approximation of the function sign(s), which is used in order to reduce the chattering effect. ε is selected as a small constant. Finally, it is now possible to completely write the input signal equation as

u = -c 0 αx 2 -c 1 βx 3 -c 2 γx 4 -B tanh s ε (82)
For the forward position and pitch angle subsystem equations, given by ( 34)-( 35), and by analogy with the system (25), we obtain that α = 1, β = g and γ = 1, also, we have that x 1 = x, x 2 = ẋ, x 3 = θ and x 4 = θ. With this information, the sliding mode surface is given by

s = c 0 x + c 1 g ẋ + c 2 θ + θ (83)
and the forward position and pitch angle control input can be expressed as

τθ = -c 0 ẋ -c 1 θ -c 2 θ -B tanh s ε (84)
The constant terms c i 's should be carefully selected to obtain a stable output. In order to stabilize x in a position outside of the origin, we must place x 1 = e x . A similar procedure should be followed to generate the lateral position and roll angle control input.

V. QUAD-ROTOR SYSTEM Fig. 3: The quad-rotor experimental platform.

The vision algorithm and the controllers presented previously have been tested over a system composed by a four-rotor helicopter, a supervisory ground station, and a wireless video and data link. The whole system is shown in Figure3.

A. Aerial vehicle

The quad-rotor helicopter shown in Figure 1 was built using a group of off-the-shelf components. Some of its characteristics are resumed in Table I. Onboard electronics consist of two interconnected cards: the first board is the control unit, while the second one deals with the motor's speed controllers. Both cards are shown in Figure 4. The control unit card performs the essential tasks of sensing, communicating and computing the control law for stabilizing the UAV attitude during fly. The characteristics of this board can be summarized as follows.

• a Texas Instruments TMS320F2812 DSP module reads the signals of the embedded sensors and computes the control law for stabilizing the aircraft. Its working frequency is 500 Hz.

• an MIDG II INS/GPS IMU from Microbotics Inc measures the angular position of the rotorcraft at a frequency of 100 Hz.

• three ADXRS150 analog gyroscopes measure the angular rates at 500 Hz. We chosen analog rate measurement rather than IMU based measurements, since we can obtain a faster refresh of angular rates which, enables a better attitude stabilization of the UAV.

• battery voltage measurement circuit intended to provide the actual tension level of the supply battery. This information is used for several goals: perform a safety landing and turn-off before an unwanted discharge of tension (avoiding accidents).

• an XBee ZB ZigBee PRO Radio Modem is used to link the ground station and the aerial vehicle. This communication link can be used to introduce external control inputs, send the sensors information to the ground station, etc. The second board contains:

• Signal conditioning circuitry: In this stage, the motor's control signals are decoupled from the rest of the electronic systems. PWM signals are also filtered and conditioned.

B. Supervisory ground station

The supervosory ground station consists of a desktop PC, a flight simulator Cyborg-X joystick, an XBee ZB ZigBee PRO Radio Modem and a Diversity video receiver system. This ground station runs a supervisory control application allowing an user to send information to the helicopter and to chose between a manual control or an autonomous visionbased position hold. The supervisory ground station receives and saves data sent by the vehicle in order to debug and analyze the flight experiments. The control feedback between the supervisory ground station and the helicopter is performed at 30 Hz.

C. Vision system

The UAV vision system is shown in Figure 5. Real-time video is obtained by means of a high definition CTDM-5351 camera, with a resolution of 640 × 480 pixels. It is installed in the lower part of the helicopter and is placed pointing downwards. The camera is connected to a 200mW micro video and audio HF transmitter. Images from the vision system are recovered on ground by a 4-Antenna Diversity System Receiver. This receiver is connected to the supervisory ground station PC throughout a USB frame grabber. The frequency of the video transmission is performed at a rate of 30 Hz. Once received, real-time video is processed by a computer vision application developed specifically for computing the 3-dimensional position and the translational speed of the helicopter (see Section III). The computer vision application is programmed in Visual C++ and is based on OpenCV functions [START_REF]Open Computer Vision Lib[END_REF].

VI. EXPERIMENTAL RESULTS

The purpose of our experiments consists on validating the most appropriate control strategy for stabilizing the position of a quad-rotor UAV in hover flight. The quad-rotor is equipped with a vision system for estimating its relative 3-dimensional position (x, y, z) as well as its translational velocity ( ẋ, ẏ, ż), with respect to a landing pad on the ground. Three similar experiments were performed, which are explained next.

Once the UAV is located exactly on top of the landing pad, the operator uses the supervisory ground station to define the current vehicle's position as the desired (x d , y d ) position reference. The desired altitude z d is always fixed at 150cm, and the desired yaw angle ψ is fixed at 0 degrees. The parameter values used for the altitude and the yaw controller are: k pz = 0.68, k vz = 1.6, k pψ = 38, k pψ = 1350. The parameters used for the nested saturations controller are: σ b4 = 0.4700, σ b3 = 0.2349, σ b2 = 0.1174 and σ b1 = 0.0287. The parameters used for the backstepping controller are: k 1 = 0.002, k 2 = 0.004, k 3 = 9 and k 4 = 4. The parameters Figures 6789show the obtained behavior when applying the three controllers to the UAV. We can observe that all the controllers achieve hovering flight, however, smoother translational and angular behaviors are obtained when using the nested saturation controller. Table II and Table III show the mean and standard deviation values for the position and Euler angles signals respectively. Note that Table II and Table III were computed with only one experiment for each controller, considering that the UAV is in steady state response. We have also computed the Mean Square Errors for the Euler angles. The values are show in Table IV. This results shown that the Nested Saturations controller is the method that induces less angular corrections, which can be considered as less control inputs generated during flight, and consequently, less energy consumption.

Note in Table II and Table III that the values for the backstepping controller are closer to the desired reference values. However, if the important objective concerns the energy consumption, the nested saturations controller should be considered as the best option. A video of the experiments can be seen at http://www.youtube.com/watch?v=SQlSXruTnj0 VII. CONCLUSIONS Three control strategies were implemented and tested in a real-time application for a quad-rotor UAV by using visual feedback. The helicopter uses a combination of an IMU and three analog rate gyros to measure its angular dynamics. Using an homography estimation technique, an onboard monocular vision system was used to estimate the (x, y, z) 3-dimensional position of the aerial vehicle with respect to a landing pad. It was also implemented an optical flow technique for estimating the vehicle's ( ẋ, ẏ, ż) translational velocities. The control algorithms were implemented onboard to stabilize the UAV's 3-dimensional position and attitude. All control algorithms ensure that the Euler angles of the vehicle remain very close to the desired values. The experimental results show that the nested saturations control approach is the most appropriated strategy for our system, since it ensures a smoother vehicle behavior and reduces the energy consumption with respect to the other two controllers.

Fig. 1 :

 1 Fig. 1: The four-rotor aircraft experimental platform.

  . The position of the vehicle's center of gravity, with respect to the inertial frame, is denoted by ξ = [x y z] T ∈ R 3 , the three Euler angles (roll, pitch and yaw), which represent the orientation of the vehicle are expressed as η = [φ θ ψ] T ∈ R 3 . The full model rotorcraft dynamics is obtained from Euler-Lagrange equations [16]:

Fig. 2 :

 2 Fig. 2: Visual system setup.

Fig. 4 :

 4 Fig. 4: The electronics on board: (a) Electronic card for the DSP, analog rate gyros, IMU connections and wireless modem; (b) Signal conditioning board.

Fig. 5 :

 5 Fig. 5: The UAV vision system: (a) CTDM-5351 camera on board; (b) 4-Antenna Diversity System Receiver.

Fig. 6 :Fig. 7 :

 67 Fig. 6: Behavior of the X, Y and altitude error signals. The Backstepping, Nested Saturation and Sliding Modes are represented by solid line, dashed line and dotted line respectively.

Fig. 8 :

 8 Fig. 8: Behavior of the pitch angle.

Fig. 9 :

 9 Fig. 9: Behavior of the roll angle.

TABLE I :

 I Characteristics of the quad-rotor

	Parameter	Value
	Diameter between rotors	40 cm
	Weight	800 grams
	Autonomy	15 minutes
	Power	12 V, 2200 mAh Li-Po battery
	Motor	1200 Brushless

TABLE II :

 II Mean values of position and Euler angles.

	Parameter	Backstepping	Nested	Sliding
			saturations	modes
	Yaw angle	-0.0108 •	-0.2504 •	-0.0400 •
	Pitch angle	0.0903 •	-0.0706 •	0.2033 •
	Roll angle	-0.0315 •	0.4371 •	-0.0782 •
	e Z error	0.0176 cm	1.4550 cm	-6.8905 cm
	X position	-1.5662 cm	-2.7342 cm	4.3211 cm
	Y position	5.6848 cm	4.2451 cm	10.1918 cm

TABLE III :

 III Standard deviation of position and Euler angles.

	Parameter	Backstepping	Nested	Sliding
			saturations	modes
	Yaw angle	0.7968 •	1.2606 •	1.0260 •
	Pitch angle	1.1588 •	0.4473 •	0.9833 •
	Roll angle	1.1673 •	0.7635 •	1.0739 •
	e Z error	8.0868 cm	8.2162 cm	7.7071 cm
	X position	9.2075 cm	6.2887 cm	10.2763 cm
	Y position	7.1207 cm	4.1411 cm	9.2603 cm

ACKNOWLEDGEMENTS

The authors gratefully acknowledge the reviewer's comments.


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

Ludovic Perret
email: ludovic.perret@lip6.fr

Christophe Petit
email: christophe.petit@uclouvain.be

Guénaël Renault
email: guenael.renault@lip6.fr

Improving the Complexity of Index Calculus Algorithms in Elliptic Curves over Binary Fields

Keywords: Elliptic Curve Cryptography, Index Calculus, Polynomial System Solving. F.R.S

The goal of this paper is to further study the index calculus method that was first introduced by Semaev for solving the ECDLP and later developed by Gaudry and Diem. In particular, we focus on the step which consists in decomposing points of the curve with respect to an appropriately chosen factor basis. This part can be nicely reformulated as a purely algebraic problem consisting in finding solutions to a multivariate polynomial f (x1, . . . , xm) = 0 such that x1, . . . , xm all belong to some vector subspace of F2n /F2. Our main contribution is the identification of particular structures inherent to such polynomial systems and a dedicated method for tackling this problem. We solve it by means of Gröbner basis techniques and analyze its complexity using the multi-homogeneous structure of the equations. A direct consequence of our results is an index calculus algorithm solving ECDLP over any binary field F2n in time O(2 ω t ), with t ≈ n/2 (provided that a certain heuristic assumption holds). This has to be compared with Diem's [14] index calculus based approach for solving ECDLP over Fqn which has complexity exp O(n log(n) 1/2 ) for q = 2 and n a prime (but this holds without any heuristic assumption). We emphasize that the complexity obtained here is very conservative in comparison to experimental results. We hope the new ideas provided here may lead to efficient index calculus based methods for solving ECDLP in theory and practice.

Introduction

Elliptic curves were independently introduced to cryptography by Miller and Koblitz in 1985 [START_REF] Miller | Use of Elliptic Curves in Cryptography[END_REF][START_REF] Koblitz | Elliptic Curve Cryptosystems[END_REF]. The security of curve-based cryptosystems often relies on the difficulty of solving the well-known Discrete Logarithm Problem on Elliptic Curves (ECDLP). Given a finite cyclic group G = P and given an element Q in G, the discrete logarithm problem asks for an integer k such that Q = kP . For an elliptic curve E defined over a finite field K, the group G can chosen to be the set E(K) of rational points on E.

One of the main method for solving (EC)DLP is Index Calculus. This approach, which was first introduced by Kraitchik [START_REF] Kraitchik | Théorie des Nombres[END_REF] and later optimized by Adleman [START_REF] Adleman | A Subexponential Algorithm for the Discrete Logarithm Problem with Applications to Cryptography[END_REF], can be seen as a general framework [START_REF] Enge | A General Framework for Subexponential Discrete Logarithm Algorithms[END_REF]. To summarize, Index Calculus algorithms are composed of the following three steps:

1. Factor Basis definition. Identify a subset F = {π 1 , . . . , π s } ⊂ G. 2. Sieving step. Collect more than #F relations of the form aP + bQ = s i=1 e i π i where a, b are random integers. 3. Linear Algebra step. From these relations, construct a matrix with a nontrivial kernel. Then, find a non-trivial vector in this kernel and deduce k (a discrete logarithm) from such vector.

While the last step is independent of the choice of G, the efficiency of the first two steps relies on specific properties of the group. Depending on the ability to find a factor basis together with an algorithm for computing relations during the sieving step, the index calculus method may achieve a subexponential complexity. For instance, subexponential algorithms have been obtained for multiplicative groups of finite fields [START_REF] Adleman | A Subexponential Algorithm for Discrete Logarithms over the Rational Subgroup of the Jacobians of Large Genus Hyperelliptic Curves over Finite Fields[END_REF][START_REF] Adleman | The Function Field Sieve[END_REF][START_REF] Adleman | Function Field Sieve Method for Discrete Logarithms over Finite Fields[END_REF][START_REF] Joux | The Function Field Sieve in the Medium Prime Case[END_REF] and for Jacobian groups of hyperelleptic curves with large genus [START_REF] Adleman | A Subexponential Algorithm for Discrete Logarithms over the Rational Subgroup of the Jacobians of Large Genus Hyperelliptic Curves over Finite Fields[END_REF][START_REF] Gaudry | An Algorithm for Solving the Discrete Log Problem on Hyperelliptic Curves[END_REF][START_REF] Gaudry | A Double Large Prime Variation for Small Genus Hyperelliptic Index Calculus[END_REF].

Related Works and Contributions. Our results are related to the sieving step of the ECDLP index calculus method proposed by Semaev [START_REF] Semaev | Summation Polynomials and the Discrete Logarithm Problem on Elliptic Curves[END_REF] and later developed by Gaudry [START_REF] Gaudry | Index Calculus for Abelian Varieties of Small Simension and the Elliptic Curve Discrete Logarithm Problem[END_REF] and Diem [START_REF] Diem | On the Discrete Logarithm Problem in Elliptic Curves[END_REF][START_REF] Diem | On the Discrete Logarithm Problem in Elliptic Curves II[END_REF]. The main idea introduced by Semaev is the use of so-called summation polynomials. As soon as a factor basis is fixed, summation polynomials can be used for sieving elements of an elliptic curve E and thus an index calculus method follows. The main problem to get an efficient index calculus is to find a factor basis together with an efficient algorithm for solving summation polynomials. During the ECC conference following publication of Semaev's result, Gaudry and Diem independently proposed such solutions. More precisely, when E is defined over an extension F q n with n > 1 a composite integer, Gaudry [START_REF] Gaudry | Index Calculus for Abelian Varieties of Small Simension and the Elliptic Curve Discrete Logarithm Problem[END_REF] proposed to use the following set F = {(x, y) ∈ E(F q n ) | x ∈ F q } as a factor basis and provides an algorithm for solving the ECDLP with a complexity better than generic methods. Next, the problem of finding P 1 , . . . , P m in F such that R = P 1 + • • • + P m for a given R ∈ E is reduced to solve the equation S m+1 ((P 1 ) x , . . . , (P m ) x , R x ) = 0, where S r is the r-th Semaev's summation polynomial [START_REF] Semaev | Summation Polynomials and the Discrete Logarithm Problem on Elliptic Curves[END_REF] and (P ) x stands for the x-coordinate of P . Diem proposed a generalization of this approach by considering (in a simpler form here) the factor basis F V := {(x, y) ∈ E(K) | x ∈ V }, with V a vector subspace of F q n /F q . Thus, the main computational tool for sieving here is an algorithm solving efficiently a specific polynomial system. Recently, Diem presented in [START_REF] Diem | On the Discrete Logarithm Problem in Elliptic Curves II[END_REF] new complexity results for this generalization. He succeeds to prove -without any heuristic assumption -some subexponential complexity results for solving ECDLP. But, for q = 2 and n is a prime -the setting considered here -he has an index calculus algorithm of exponential com-

plexity e O n log(n) 1/2
. The polynomial systems occurring in [START_REF] Diem | On the Discrete Logarithm Problem in Elliptic Curves II[END_REF] are solved with a geometrical algorithm proposed by Rojas [START_REF] Rojas | Solving Degenerate Sparse Polynomial Systems Faster[END_REF]. Whilst such algorithm has a good complexity estimate, it is well known that its hard to implement it in practice.

In this work, we focus on the specific case q = 2 and n prime. We show that the polynomial systems occurring have a very specific structure. We provide a new (heuristic) algorithm taking advantage of the structure for the sieving step. We focus our study on the following point decomposition problem related to some vector space: Problem 1 (Point Decomposition Problem associated to a vector space V ). Let V be a vector space of

F 2 n /F 2 . Given a point R ∈ E(F 2 n ), the problem is to find -if any -m points P 1 , . . . , P m such that R = P 1 + • • • + P m with the additional constraint that (P i ) x ∈ V for all i ∈ {1, . . . , m}.
This problem can be reduced to a polynomial system solving problem using Semaev's summation polynomials (or using any other polynomial system modeling). Here, we show that the multi-homogeneous structure of the system constructed from Semaev's summation polynomials can be used to design a Gröbner based algorithm. To be more precise, we design an efficient algorithm for: Problem 2. Let t ≥ 1, V be a vector space of F 2 n /F 2 and f ∈ F 2 n [x 1 , . . . , x m ] be any multivariate polynomial of degree bounded by 2 t -1 in each variable. The problem is to find (z 1 , . . . , z m ) ∈ V m such that f (z 1 , . . . , z m ) = 0.

Since F 2 n is a vector space over F 2 , f can be rewritten (or deployed) as a polynomial system of m equations over F 2 and then can be solved using Gröbner bases algorithms. The prominent observation is to remark that this system is (affine) multi-homogeneous. While the complexity of solving bi-linear systems using Gröbner bases -that is to say polynomials of bi-degree (1, 1) -is now well understood [START_REF] Faugère | Gröbner Bases of Bihomogeneous Ideals Generated by Polynomials of Bidegree (1,1): Algorithms and Complexity[END_REF], the general case is not known. Consequently, we propose a simple ad-hoc algorithm to take advantage of the multihomogeneous structure. This is of independent interest in the more general context of computer algebra.

The main idea is to show that starting from the unique equation f = 0, we can generate many low-degree equations by deploying the equations mf = 0 over F 2 for a large number of appropriately chosen monomials m. Another main difference with unstructured polynomials is that -in some degree d -the number of monomials occurring in the polynomials coming from mf is much smaller than the number of all possible monomials in degree d. Indeed, due to the choice of m, the monomials occurring in the equations constructed from mf have still a multi-homogeneous structure and their degrees are well controlled.

As usual, to estimate the maximum degree D reached during the computation we study the number of equations minus the number of monomials. When this number is > 0, and under a reasonable linear independence assumption confirmed by our experimental results, the computation is finished. Using the structure of the polynomials, we prove that this degree D is much smaller than expected; assuming that a reasonable heuristic is true. It is worth noticing that although we describe our algorithm as a linearization method [START_REF] Macaulay | The Algebraic Theory of Modular Systems[END_REF], Gröbner basis algorithms like F 4 or F 5 [START_REF] Faugère | A New Efficient Algorithm for Computing Gröbner Basis (F4)[END_REF][START_REF] Faugère | A New Efficient Algorithm for Computing Gröbner Bases without Reduction to Zero (F5)[END_REF] can be advantageously used in practice to solve the corresponding polynomial systems. The algorithm presented in this paper, together with its complexity analysis, can thus be understood as a method to (heuristically) bound the complexity of computing the corresponding Gröbner basis similarly to the Macaulay's bound obtain by Lazard [START_REF] Lazard | Gröbner-Bases, Gaussian Elimination and Resolution of Systems of Algebraic Equations[END_REF] to bound the complexity of Gröbner bases in the general case. More precisely, we obtain: Theorem 2. Assuming some linear independence assumption ( The index calculus algorithm presented here has a better complexity than the one proposed recently by Diem [START_REF] Diem | On the Discrete Logarithm Problem in Elliptic Curves II[END_REF]. Moreover, we propose a novel approach for solving the sieving step. We consider that it is a major open challenge to further exploit the intrinsic algebraic structure of Problem 2 using Gröbner bases algorithms. The complexity obtained here for solving ECDLP is still exponential. We hope that the structures identified here can be further used to get a complexity better than generic algorithms (see preliminary experiments in Section 5.2) or to get a subexponential algorithm in a near future. Finally, we emphasize that the complexity analysis of our algorithm relies on a heuristic assumption on the rank of a linearized system. The validity of this assumption was experimentally checked (see Section 5.1). Whilst we pushed the experiments as far as possible, we pointed out that -due to the size of the systems involved -it is very difficult to verify experimentally the assumption for large parameters. Consequently, it is an open issue to prove Assumption 1.

Outline. The remaining of this paper is organized as follows. In Section 2, we detail our notations and we provide some background on Gröbner bases. In Section 3, we introduce and analyze our algorithm for solving Problem 2. In Section 4, we apply our new result to the ECDLP over binary fields. In Section 5, we present experimental results to give first evidences for Assumption 1. Section 6 concludes the paper and introduces future extensions of our work.

Preliminaries

In this section, we introduce definitions, notations and recall well known results concerning polynomial system solving.

Definition and Notation

Let F 2 be the finite field of cardinality 2. We will consider a degree n extension F 2 n of F 2 . We will often see F 2 n as an n dimensional vector space over F 2 . Let {θ 1 , . . . , θ n } be a basis of F 2 n over F 2 . We will use bold letters for elements, variables and polynomials over F 2 n and normal letters for elements, variables and polynomials over F 2 . If x 1 , . . . , x m are algebraic independent variables over a finite field K, we write R = K[x 1 , . . . , x m ] for the polynomial ring in those variables. Given a set of polynomials {f 1 , . . . , f } ∈ R, the ideal generated by this set will be denoted by

f 1 , . . . , f ⊂ R. We write Res xi (f 1 , f 2 ) for the resultant of f 1 ∈ R and f 2 ∈ R with respect to the variable x i . A power product, i.e. k i=1 x ei
i where e i ∈ N, is called a monomial. Finally, we introduce a structure which will be very useful in this paper.

Definition 1 ([24]). Let X 1 ∪ X 2 ∪• • • ∪X t = {x 1 , . . . , x m } be a partition of the variables set. We shall say that a polynomial f ∈ K[x 1 , . . . , x m ] = K[X 1 , . . . , X t ] is multi-homogeneous of multi-degree (d 1 , d 2 , . . . , d t ) if: ∀(α 1 , . . . , α t ) ∈ K t , f(α 1 X 1 , . . . , α t X t ) = α d1 1 • • • α dt f (X 1 , . . . , X t ). For all i, 1 ≤ i ≤ t, let X i = {x i,1 , . . . , x i,ni }. We shall say that f is affine multi-homogeneous if there exists f h ∈ K[X 1 , . .

. , X t ] a multi-homogeneous polynomial of same degree such that when one replaces (homogenization) variables

x i,ni by 1 we obtains f , i.e.: f (x 1,1 , . . . , x 1,n1-1 , . . . , x t,1 , . . . , x t,nt-1 ) is equal to f h i (x 1,1 , . . . , x 1,n1-1 , 1, . . . , x t,1 , . . . , x t,nt-1 ,

1). Finally, we shall say that f has a multi-homogeneous structure if it is multi-homogeneous or affine multihomogeneous. A system of equation has a multi-homogeneous structure if each equation has a multi-homogeneous structure (the equations can have different multi-degrees).

Given a number e = ∞ i=0 e i 2 i with e i ∈ {0, 1}, we define its Hamming weight as the number of non-zero elements in its binary expansion, i.e. W (e) := ∞ i=0 e i . We write n k for the number of choices of k elements among a set of n elements without repetition. We write O for the "big O" notation: given two functions f and g of n, we say that

f = O(g) if there exist N, c ∈ Z + such that n > N ⇒ f (n) ≤ cg(n).
The notation log stands for the binary logarithm. Finally, we write ω for the linear algebra constant. Depending on the algorithm used for linear algebra, we have 2.376 ≤ ω ≤ 3.

Gröbner Bases [10]

Recent methods such as Faugère's F 4 and F 5 [START_REF] Faugère | A New Efficient Algorithm for Computing Gröbner Basis (F4)[END_REF][START_REF] Faugère | A New Efficient Algorithm for Computing Gröbner Bases without Reduction to Zero (F5)[END_REF] algorithms reduce Gröbner basis computation to Gaussian eliminations on several matrices. The link between linear algebra and Gröbner bases has been established by Lazard [START_REF] Lazard | Gröbner-Bases, Gaussian Elimination and Resolution of Systems of Algebraic Equations[END_REF]. He showed that computing a Gröbner basis is equivalent to perform Gaussian elimination on the so-called Macaulay matrices as defined below:

Definition 2 (Macaulay Matrix [START_REF] Macaulay | The Algebraic Theory of Modular Systems[END_REF][START_REF] Macaulay | Some Properties of Enumeration in the Theory of Modular Systems[END_REF]).

c 1 i,j c 2 i,j . . . ⎛ ⎝ ⎞ ⎠ m 1 > m 2 > . . . . . . t i,j f i Let F = {f 1 , . . . , f } ⊂ R be a set of polynomials of degree ≤ d. Let B = {m 1 > m 2 > • • • } ⊂ R
i ∈ F . Rows of M d (F ) correspond to the coefficients vectors (c 1 i,j , c 2 i,j , . . .) of these polynomials t i,j f i = k c k i,j m k with respect to the basis B.
Precisely, Lazard [START_REF] Lazard | Gröbner-Bases, Gaussian Elimination and Resolution of Systems of Algebraic Equations[END_REF] proved the following fundamental result:

Theorem 1. Let F = {f 1 , . . . , f } ⊂ R.
There exists a positive integer D for which Gaussian elimination on all matrices M 1 (F ), M 2 (F ), . . . , M D (F ) computes a Gröbner basis of f 1 , . . . , f .

F 4 [START_REF] Faugère | A New Efficient Algorithm for Computing Gröbner Basis (F4)[END_REF] can be seen as another way to use linear algebra without knowing an a priori bound. It successively constructs and reduces matrices until a Gröbner basis is found. The same is true for F 5 when considered in "F 4 -style" [START_REF] Faugère | Solving Systems of Polynomial Equations with Symmetries using SAGBI-Gröbner bases[END_REF].

It is clear that an important parameter in Gröbner basis computation is the maximal degree D reached during the computation. This maximal degree is called the degree of regularity. However, it is a difficult problem to estimate a priori the degree of regularity. This degree is known and well mastered for specific families of systems called regular and semi-regular [START_REF] Bardet | Étude des Systèmes Algébriques Surdéterminés. Applications aux codes Correcteurs et à la Cryptographie[END_REF][START_REF] Bardet | Complexity of Gröbner Basis Computation for Semi-Regular Overdetermined Sequences over F2 with Solutions in F2[END_REF][START_REF] Bardet | Asymptotic Expansion of the Degree of Regularity for Semi-Regular Systems of Equations[END_REF]. It is classical to assume that the regularity of regular/semi-regular systems provides an extremely tight upper bound on the regularity of random system of equations (most of the times, we have equality). For example, the regularity degree of a regular sequence

f 1 , . . . , f ∈ R (with ≤ n) is D = 1 + i=1 (deg(f i ) -1)
. Ideals with special structures -typically arising in cryptography -may have a much lower regularity degree, hence a much better time complexity. This has permitted Gröbner bases techniques to successfully attack many cryptosystems (e.g. the Hidden Field Equation cryptosystem (HFE) [START_REF] Patarin | Hidden Fields Equations (HFE) and Isomorphisms of Polynomials (IP): Two New Families of Asymmetric Algorithms[END_REF][START_REF] Kipnis | Cryptanalysis of the HFE Public Key Cryptosystem by Relinearization[END_REF][START_REF] Faugère | Algebraic Cryptanalysis of Hidden Field Equation (HFE) Cryptosystems Using Gröbner Bases[END_REF][START_REF] Granboulan | Inverting HFE Is Quasipolynomial[END_REF], Multi-HFE [START_REF] Bettale | Cryptanalysis of HFE, multi-HFE and Variants for Odd and Even Characteristic[END_REF]). In many cases, the algebraic systems appearing in these applications were not generic and could be solved more efficiently than generic systems, sometimes using dedicated Gröbner basis algorithms.

In our case, the algebraic system considered in Problem 2 has a multihomogeneous structure (more precisely multi-linear ). Interestingly, the formal study of computing Gröbner basis of multi-homogeneous systems has been initiated by Faugère, Safey El Din and Spaenlehauer for bilinear systems [START_REF] Faugère | Gröbner Bases of Bihomogeneous Ideals Generated by Polynomials of Bidegree (1,1): Algorithms and Complexity[END_REF] leading already to new cryptanalytic results [START_REF] Faugère | Algebraic Cryptanalysis of McEliece Variants with Compact Keys[END_REF][START_REF] Faugère | Computing Loci of Rank Defects of Linear Matrices using Gröbner Bases and Applications to Cryptology[END_REF][START_REF] Faugère | Cryptanalysis of MinRank[END_REF]. However, the general case (more blocks, larger degrees) remains to be investigated.

As a consequence, we design in this paper a simple ad hoc algorithm to take advantage of the multi-homogeneous structure arising in Problem 2. Basically, the idea is to generate a submatrix of the Macaulay matrix at a degree D Lin which allows to linearize the system derived from Problem 2. To do so, we strongly exploit the multi-homogeneous structure of the equations. The fundamental remark is that many low-degree relations exist and can be explicitly predicted. Moreover, the number of columns in the Macaulay matrix is less than usual. This is due to the fact that all monomials occurring have still the multihomogeneous structure. Roughly, this allows to predict many zero columns for a suitably chosen subset of the rows. In section 3.2, we derive a bound on the degree D Lin which needs to be considered.

It is worth noticing that although we describe our algorithm as a linearization method, computing a Gröbner basis with F 4 or F 5 [START_REF] Faugère | A New Efficient Algorithm for Computing Gröbner Basis (F4)[END_REF][START_REF] Faugère | A New Efficient Algorithm for Computing Gröbner Bases without Reduction to Zero (F5)[END_REF] can be advantageously used in practice to solve the corresponding polynomial systems. The D Lin obtained has to be understood as an upper bound on the real regularity degree of the system. Any new result on multi-homogeneous systems would allow to improve the complexity of solving Problem 2 (partial results are presented in Section 5.2).

Solving Multivariate Polynomials with Linear Constraints

In this section, we describe the main result of this paper: an algorithm for solving Problem 2. Let V be a F 2 -vector subspace ⊂ F 2 n of dimension n and let f ∈ F 2 n [x 1 , . . . , x m ] be a multivariate polynomial with degree < 2 t in each variable. We want to solve f (x 1 , . . . , x m ) = 0 under the linear constraints x 1 , . . . , x m ∈ V .

To tackle this problem, we generalize the algorithm and analysis of [START_REF] Faugère | New Subexponential Algorithms for Factoring in SL(2, F n 2 )[END_REF] from the multi-linear case (t = 1) to arbitrary values of t. From now on, we assume that m • n ≈ n so that the problem has about one solution on average.

Modeling the Linear Constraints

Let {ν 1 , . . . , ν n } ⊂ F 2 n be a basis of V as a F 2 -vector space. Let y i,j be m • n variables defined by

x i = ν 1 y i,1 + ν 2 y i,2 + • • • + ν n y i,n .
We apply a Weil descent to Problem 2 (see e.g. [START_REF]Handbook of Elliptic and Hyperelliptic Curve Cryptography[END_REF]Chapter 7]). By replacing the variables x i in the polynomial f , we get a new polynomial

f V ∈ F 2 n [y 1,1 , . . . , y m,n ] with m • n variables.
The linear constraints on f are translated to Galoisian constraints by constraining the solutions of f V to F 2 . Using the field equations, f V is viewed more precisely as a a polynomial in the affine algebra A(F

2 n ) := F 2 n [y 1,1 , . . . , y m,n ]/ S fe ,
where S fe is the 0-dimensional ideal generated by the field equations:

S fe = {y 2 i,j -y i,j } 1≤j≤n 1≤i≤m . Problem 2 is then equivalent to compute a Gröbner basis of the ideal f V , S fe ⊂ F 2 n [y 1,1 , . . . , y m,n ].
It is generally more efficient to consider its resolution over F 2 . To do so, we consider A(F 2 n ) as a module over A(F 2 ) = F 2 [y 1,1 , . . . , y m,n ]/ S fe whose basis is {θ 1 , . . . , θ n }. We consider f V as an element of A(F 2 n ) and we deploy it as a A(F 2 )-linear combination of the basis {θ 1 , . . . , θ n }. Namely:

f V = [f V ] ↓ 1 θ 1 + [f V ] ↓ 2 θ 2 + • • • + [f V ] ↓ n θ n (1)
for some

[f V ] ↓ 1 , . . . , [f V ] ↓ n ∈ A(F 2
) that depend on f and the vector subspace V . Due to the linear independence of the θ 1 , . . . , θ n , Problem 2 is equivalent to solve:

S alg : [f V ] ↓ 1 = [f V ] ↓ 2 = • • • = [f V ] ↓ n = 0. (2)
In order to solve S alg , we will generate many new equations by deploying multiples of f ∈ F 2 n [x 1 , . . . , x m ] over the vector subspace V . The key point of this strategy is the existence of abnormally high number of low-degree equations arising as algebraic combinations of the equations in [START_REF] Adleman | The Function Field Sieve[END_REF].

From now on, we represent the classes of polynomials g V ∈ A(F 2 n ) and

[g V ] ↓ i ∈ A(F 2 ) corresponding to g ∈ F 2 n [x 1 , . . . ,
x m ] by their minimal elements, in other words by their normal forms modulo I fe (whose generators form a Gröbner basis). By abuse of notation, we use the same symbol for a class and its minimal representative in the underlying polynomial ring. When the context is not clear, we precise the algebra where the element is lying.

Low-Degree Equations

Let e 1 , . . . , e m ∈ N and let m = m i=1 x i ei be a monomial of F 2 n [x 1 , . . . , x m ]. Following the descent described in Section 3.1, we have

A(F 2 n ) (mf ) V = n k=1 [(mf ) V ] ↓ k θ k , with [(mf ) V ] ↓ k ∈ A(F 2 n ).
The equation f = 0 clearly implies mf = 0, hence [(mf ) V ] ↓ k = 0 for k = 1, . . . , n. We can then add these new equations to the polynomial system (2). The equations obtained in this way all share the same structure. More precisely, their minimal representatives, due to the normal form computation modulo I fe , are all affine multi-linear in F 2 n [y 1,1 , . . . , y m,n ]. Moreover, thanks to the evaluations done during the deployments, each block of variables X i = {y i,1 , . . . y i,n } naturally corresponds to the variable x i . From these structures, we deduce the following result.

Lemma 1. Let f ∈ F 2 n [x 1 , . . . , x m ] be a multivariate polynomial with degree < 2 t in each variable. Let e 1 , . . . , e m ∈ N and m = m i=1 x i ei be a monomial of F 2 n [x 1 , . . . , x m ]. There exist polynomials p j,k ∈ F 2 [y 1,1 , . . . , y m,n ] such that [(mf ) V ] ↓ k = n j=1 p j,k [f V ] ↓ j . Each polynomial p j,k has degree ≤ W (e i ) with re- spect to every block of variables X i = {y i,1 , . . . y i,n }, 1 ≤ i ≤ m. Moreover, each minimal representative of [(mf ) V ] ↓ k has degree ≤ max 0≤e i <2 t W (e i + e i ) w.r.t. each block of variables X i , 1 ≤ i ≤ m.
This lemma implies that the new equations (obtained from mf ) are algebraic combinations of the original ones (obtained from f ). In particular, they can a priori be recovered "in a hidden form" with any Gröbner basis algorithm at degree D apriori = mt + m j=1 W (e j ). The value D apriori is the degree that the equations [(mf ) V ] ↓ k should have a priori from the algebraic dependencies of Lemma 1. It is the sum of the degree of the deployments of f (at most mt) and the degree of each polynomial p j,k at most m j=1 W (e j ) . However, Lemma 1 also implies that the [(mf ) V ] ↓ k only have degree D actual = m j=1 (max 0≤e j <2 t W (e j + e j )). Thus:

D apriori -mt ≤ D actual ≤ D apriori .
Therefore, [(mf ) V ] ↓ k may have a degree drop as large as mt depending on the monomial m chosen. The existence such low-degree relations compared generic systems makes Gröbner basis algorithms faster in practice and allows a linearization strategy.

Following the general method of Macaulay [START_REF] Macaulay | The Algebraic Theory of Modular Systems[END_REF], we will linearize the polynomial system S alg ∪ {. . . , . . . , [(mf

) V ] ↓ 1 , [(mf ) V ] ↓ 2 , . . . , [(mf ) V ]
↓ n , . . . , . . .} using the low-degree equations identified in this section. The choice of the monomials m used to generate the equations are particularly important for the efficiency of the linearization strategy. In particular, the equations with the lowest degrees are the most interesting ones since they involve less monomial terms. Of course, this strategy requires that a substantial subset of all low-degree relations are linearly independent.

Linear Dependencies

In the previous section, we explained how low-degree relations can be produced. To be used in a linearization strategy, these equations must be linearly independent. In this section, we describe two sources of linear dependencies.

Frobenius Transforms. The first source of linear dependencies is due to the Frobenius endomorphism (as identified in [START_REF] Faugère | New Subexponential Algorithms for Factoring in SL(2, F n 2 )[END_REF]). Let {θ 1 , . . . , θ n } be a basis of F 2 n /F 2 . The set {θ 2 1 , . . . , θ 2 n } is another basis of

F 2 n /F 2 . Let a ij ∈ F 2 be such that θ 2 j = i a ij θ i . We have f V 2 = n j=1 f V 2 ↓ j θ j . However, f V 2 = n j=1 [f V ] ↓ j θ j 2 = n j=1 [f V ] ↓ j θ 2 j = n i=1 n j=1 a ij [f V ] ↓ j θ i . Thus, we obtain f V 2 ↓ i = n j=1 a ij [f V ] ↓ j . In other words, the polynomials f V 2 ↓ 1 , . . ., f V 2 ↓ n are linear combinations of [f V ] ↓ 1 , . . . , [f V ] ↓ n . Decomposing f V as a sum of monomials, we deduce that f V 2 ↓ i = m∈Mon(fV) [(mf ) V ] ↓ i = n j=1 a ij [f V ] ↓ j
. This provides a non trivial linear relation between some low-degree equations. More generally, we obtain similar relations if we replace f by m f in the above equation (for any monomial m ). These linear dependencies can be easily detected and prevented. Indeed, we can simply avoid every monomial m that is the leading term of (m f ) for some m .

Vector Dependencies. Another source of dependency, that we call vector dependencies, is induced by the vector space V . To illustrate the phenomena, consider the simplest polynomial g

= x 1 ∈ F 2 n [x 1 , . . . , x m ]. We have g V = ν 1 y 1,1 + ν 2 y 1,2 + • • • + ν n y 1,n ∈ A(F 2 n ). More generally, (g 2 i ) V = ν 2 i 1 y 1,1 + ν 2 i 2 y 1,2 + • + ν 2 i n y 1,n ∈ A(F 2 n ).
Since y 1,1 , . . . , y 1,n are linearly independent, we obtain for any k > n a non trivial linear dependency by considering different

g 2 i 's, i.e. ∃β 1 , . . . , β k ∈ F 2 n \ {(0, . . . , 0)} such that β 1 (g 2 i 1 ) V + β 2 (g 2 i 2 ) V + • • •+ β k (g 2 i k ) V = 0.
= m i=1 x i ei ∈ F 2 n [x 1 , . . . , x n ] with 1 ≤ e i < 2 n , such that all [(mf ) V ] ↓ k (1 ≤ k ≤ n) are in A(F 2 ) d . Finally, let E(d) := n • #Mon(d) and M (d) := #MLinMonB(d).
We are now ready to describe Algorithm 1, a simple linearization algorithm for solving Problem 2. The algorithm constructs a sub-matrix of the Macaulay (see Definition 2) matrix M d for system 2. We first gather

M (d) equations [(mf ) V ] ↓ i with m ∈ Mon(d).
By definition, all these equations are in A(F 2 ) d . Hence, they can be decomposed with respect to the basis MLinMonB(d). We then form the corresponding linear system S lin over F 2 , where each row corresponds to the coefficients involves in the equations and each column corresponds to an element in MLinMonB(d). We finally solve the linear system. This simple algorithm, that we call Sub-Macaulay, is not aimed to be optimal in practice but to derive complexity bounds.

The general linearization strategy and our analysis below rely on a heuristic assumption formalized below: Assumption 1. With a probability exponentially close to one, the equations generated by Algorithm 1 are linearly independent.

Particularly, the assumption states that the solutions of S lin are in one-to-one correspondence with the solutions of Problem 2.

Complexity Bounds for Solving Problem 2

We now derive an upper bound on the complexity of Algorithm 1. The main task is to estimate the values of M (d) (number of columns in S lin ) and E(d) (number of equations in S lin ). Due to the field equations, we only have multilinear monomials, i.e. variables can only have exponents 0 or 1. Therefore, the number of monomials of total degree d , 0 ≤ d ≤ d involving variables of the block X i is n d . For the m blocks, we get:

M (d) = d d =0 n d m . ( 3
)
Input: f ∈ F2n [x1, . . . , xm] of degree in each variable is bounded by 2 t -1, and V a F2-vector subspace ⊂ F2n of dimension n . Result: If not empty, the finite set {s1, . . .} ⊂ V m such that f (si) = 0. begin By definition, the degree of each variable x i occurring in the monomials of f may have any degree between 0 and 2 t -1. By Lemma 1, the degree of

(mf )V] ↓ 1 , . . . , [(mf )V] ↓ n ∈ F2[y1 1, . . . , ym nm ] be such that 6 (mf )V = n k=1 [(mf )V] ↓ k θ k ; S alg ← S alg ∪ ([(mf )V] ↓ 1 , . . . , [(mf )V] ↓ n ); Mon ← Mon ∪ [m];
[(mf ) V ] ↓ k with respect to X i is max 0≤e i <2 t W (e i + e i ) = W e i 2 t
+ t. Therefore, the number of exponents e i , 1 ≤ e i ≤ 2 n leading to degree d with respect to the block X i is

2 t n -t d -t . As a consequence 1 E(d) := n 2 tm d d =t n -t d -t m . ( 4
)
We derive the following asymptotic bound (proven in Appendix A) on the minimal value d that allows linearization.

Lemma 2. Let α be such that 1 -α < 1/2 < α < 1. Assuming that n = n α , m = n 1-α and t = m -1, then E(d) ≥ M (d) for d ≈ n α 2 when n is large enough.
In the next

log(M (d)) ≈ m log n α n β ≈ n (1-α) n β (α -β) log(n) = n (1-1 log(n) ) = n 2 .
4 Application to ECDLP over Binary Fields

Diem's Variant of Index Calculus

Let E be an elliptic curve defined over F 2 n by the equation

E : y 2 + xy = x 3 + x 2 + a 6 , for some a 6 ∈ F 2 n . (5)
Semaev' summation polynomials S r [START_REF] Semaev | Summation Polynomials and the Discrete Logarithm Problem on Elliptic Curves[END_REF] are multivariate polynomials with the following property: S r (x 1 , . . . , x r ) = 0 for some x 1 , . . . , x r ∈ F 2 n if and only if there exist y 1 , . . . , y r ∈ F 2 n such that (x i , y i ) ∈ E(F 2 n ) and (x 1 , y 1 ) + • • • + (x r , y r ) = P ∞ .

Proposition 1 ( [START_REF] Semaev | Summation Polynomials and the Discrete Logarithm Problem on Elliptic Curves[END_REF]). The summation polynomials of the elliptic curve (5) are recursively given by: S 2 (x 1 , x 2 ) := x 2 +x 1 , S 3 (x 1 , x 2 , x 3 ) := x 1 2 x 2 2 +x 1 2 x 3 2 + x 1 x 2 x 3 +x 2 2 x 3 2 +a 6 and for r ≥ 4 and any k, 1 ≤ k ≤ r-3, the r-th summation polynomial is S r (x 1 , . . . , x r ) := Res X S r-k (x 1 , . . . , x m-k-1 , X), S k+2 (x r-k , . . . , x r , X) .

Moreover, the polynomial S r is symmetric and has degree 2 r-2 in each variable x i as soon as r ≥ 2. The cost of constructing this polynomial is bounded by 2 O((r-1) 2 ) .

Following Diem [START_REF] Diem | On the Discrete Logarithm Problem in Elliptic Curves II[END_REF], we use summation polynomials in the sieving stage of an index calculus algorithm. Let V be a vector subspace of F 2 n /F 2 with a dimension n to be fixed later. We define the factor basis F V as

F V := {(x, y) ∈ E(F 2 n )|x ∈ V }.
Since the abscissas of points ∈ E are uniformly distributed in F q n [START_REF] Diem | On the Discrete Logarithm Problem in Elliptic Curves[END_REF][START_REF] Diem | On the Discrete Logarithm Problem in Elliptic Curves II[END_REF], we can assume that the set F V has size about 2 n . During the sieving stage, we compute about 2 n relations P ∞ = a i P + b i Q + Pj ∈FV e i j P j with P j ∈ F V for randomly chosen integer couples (a i , b i ). Each relation is obtained by solving an instance of the following problem, for some integer parameter m to be fixed later. Problem 3. Let a i , b i be fixed random integers and R = a

i P + b i Q. Find -if any -(x 1 , . . . , x m ) ∈ V m such that S m+1 (x 1 , . . . , x m , (R) x ) = 0.
Clearly, this problem is a particular instance of Problem 2.

A Linearization Strategy for Solving ECDLP over F 2 n

We now apply the analysis of Section 3 to Problem 3. Let α, 1/2 < α < 1. be a parameter that will be optimized later. We set n := n α and m := n 1-α as in Lemma 2. According to Proposition 1, the (m + 1)th Semaev's polynomial S m+1 can be computed in time O(2 t1 ), where

t 1 ≈ m 2 ≈ n 2(1-α) .
For each relation computed in the sieving stage, we generate and solve an instance of Problem 2 where f has degree 2 m-1 with respect to each of the m variables. According to Theorem 2, each instance can be solved in time O(2 ωτ ) and memory O(2 2τ ), where τ ≈ n/2. The probability that a point R can be written as a sum of m factor basis elements is 1 m! [START_REF] Diem | On the Discrete Logarithm Problem in Elliptic Curves II[END_REF]. Hence, we need m!2 n trials on average to obtain 2 n valid relations. Since log(m!) ≈ n (1-α) log n (1-α) , the total cost of the sieving stage is bounded by O(2 t2 ) where

t 2 ≈ n (1-α) log n (1-α) + n α + ω n 2 .
Finally, the last step of our algorithm consists in (sparse) linear algebra on a matrix of rank about 2 n with elements of size about n bits. The computation time of this part can be approximated by O(2 t3 ), where t 3 ≈ ω n α , ω being the sparse linear algebra constant. Finally, we obtain the following theorem. This has to be compared with the algorithm presented [START_REF] Diem | On the Discrete Logarithm Problem in Elliptic Curves II[END_REF] which is was so far the best algebraic index calculus based approach. For q = 2 and n prime, the algorithm of [START_REF] Diem | On the Discrete Logarithm Problem in Elliptic Curves II[END_REF] has complexity e O n log(n) 1/2 (but such complexity holds without any heuristic assumption). The complexity of our approach is also very close to O(2 n/2 ); the complexity of generic algorithms (e.g. Pollard's rho algorithm).

Experimental Results

Validation of the Heuristic Assumption

We have experimentally checked the validity of Assumption 1 using the computer algebra system Magma. During the sieving stage, most of polynomial systems generated have no solution (only 1/m! polynomial systems will produce a solution). Thus we mainly check the assumption for polynomial systems with no solution in the vector space V . In order to validate this assumption we proceeded in the following way. We chose many random polynomials f with degree 2 t < 2 m-1 in each of its m < 5 variables. The coefficients are in F 2 n with n a prime less than 40. Then, for each of these polynomials, we construct a large binary matrix M of size N × M . This matrix represents the deployed polynomials [(mf ) V ] ↓ 1 , . . . , [(mf ) V ] ↓ n for all monomials in Mon(d) with d is the smallest integer with M (d) < E(d) . We avoid the ones corresponding to possible linear dependencies as identified in Section 3.3. We want to demonstrate that a random square submatrix of size M × M of M is full rank. We recall that the probability that a random N × M boolean matrix has (2 r -2 j ). Hence P (100, 100, 100) is only 28.8% but P (105, 100, 100) = 96.9%. For this reason, we consider submatrices M of M of size (M + 5) × M . We check that the rank is M or M -1 or M -2; for a random Boolean matrix the probability is 99.999982%. We repeated the test 100 times and deduced an approximation of the success probability. In all our experiments we always obtain ≈ 100% of success. In particular, we validated the assumption for Problem 3 using Semaev's summation polynomials with m + 1 variables (m = 2, . . . , 4). These validations represent a huge computational effort since some of the matrices M encountered during the experiments had more than 10000 columns.

Gröbner Basis Computations

We performed actual Gröbner basis computation using the FGb software [START_REF] Faugère | FGb: A Library for Computing Gröbner Bases[END_REF] to compare the theoretical number of operations with a more realistic value. For instance, for n = 131, m = 2 we can solve the point decomposition problem in 2 74.5 operations using a variant of the hybrid method [START_REF] Bettale | Hybrid Approach for Solving Multivariate Systems over Finite Fields[END_REF]. We compare the cost of solving Semaev's equations using Algorithm 1 (complexity of solving a linear system of rank M (d)) and a Gröbner basis computation. This is the dominating part in the complexity of our index calculus based algorithm. Remark that the number of operations reported in this table are much below the theoretical estimate for Algorithm 1. We have a gain of (at least) a factor 2 in the exponent in favor of Gröbner computations. Further experiments suggest that a more advanced approach (i.e. Gröbner basis instead of linearization) can lead to an algorithm of better complexity. One can hope a gain of a factor m in the exponent. However, we are, for the moment, not able to provide theoretical evidences supporting such an improvement.

Conclusion and Perspectives

As a conclusion, we emphasize that the algorithm of Section 3 is of independent interest in the more general context of polynomial system solving. It shows that algebraic systems arising by deploying a multivariate polynomial equation from F 2 n to F 2 are easier to solve than generic systems. We have underlined the intrinsic structures which help for solving such systems. This is an open problem how to use such structure in an F 4 /F 5 based algorithm. We hope that such improvements may lead to a subexponential algorithm. Finally, the approach generalizes quite easily to other composite fields with "small" characteristics, resulting in similar algorithms with comparable asymptotic complexities. All these extensions will be discussed in an extended version. Finally, although the paper is mainly theoretical, we hope that it could be the building block towards the development of more efficient methods to solve the ECDLP problem.

1

  Let d be the smallest integer such that E(d) ≥ M (d); 2 Mon ← []; S alg ← [] // Empty lists of monomials and equations 3 for k = 1, . . . , E(d)/n do 4 Randomly pick a monomial m ∈ Mon(d) \ Mon; 5 Let [

7 end 8 Construct 9 If 11 Algorithm 1 .

 789111 S lin the linear system over F2 obtained by linearizing S alg ; S lin has solutions then solve S lin and return {s1, . . . , } else return no 10 solution end; end Sub-Macaulay

1 j=0( 2 N

 12 rank r is P (N, M, r) = 2 -N M r--2 j )

  see Assumption 1, p. 36), Problem 2 can asymptotically be solved in time O(2 ωτ ) and memory O(2 2τ ), where ω is the linear algebra constant and τ ≈ n/2. Under the same hypothesis, there exists an index calculus based algorithm solving ECDLP over F 2 n in the same time complexity.

  be the sorted set (w.r.t. a fixed monomial ordering) of degree ≤ d monomials. The set B is a basis of the vector space of degree ≤ d polynomials in R. The Macaulay matrix M d (F ) of degree d is defined as follows. We consider all the polynomials t i,j f i of degree ≤ d with t i,j ∈ B and f

  This simple example can be easily generalized to g = mf with m a monomial. Such linear dependency can be clearly prevented during the generation of equations [(mf ) V ] ↓ i 's by considering monomials m = m i=1 x i ei , with 0 ≤ e i < 2 n < 2 n . 3.4 Description of the Linearization Algorithm For any positive integer d, let MLinMonB(d) be the set of multi-linear monomials in F 2 [y 1,1 , . . . , y m,n ] of degrees ≤ d with respect to each block X i = {y i,1 , . . . , y i,n }. The image of MLinMonB(d) in A(F 2 ) is a basis of the vector subspace A(F 2 ) d composed of elements in A(F 2 ) with a minimal representative having degrees ≤ d with respect to each block X i . Let also Mon(d) be the set of monomials m

  table, we have computed the smallest d real such that E(d real ) ≥ M (d real ) for different values of n and α. Then, we compute β real = log n (d real ). According to Lemma 2, the theoretical value predicted for β = log n (d) is β theo = α -1 log(n) . The last column columns shows that β theo is extremely close to β real . Finally, we obtain an estimate on the complexity: Theorem 2. Let α, n , m, t be as in Lemma 2. Under Assumption 1, Problem 2 can asymptotically be solved in time O(2 ωτ ) and memory O(2 2τ ), where ω is the linear algebra constant and τ ≈ n 2 .

	n	d real	α	β real	α -1 log(n) -β real
	100000	1114 2/3 0.6093	-0.0029
	1000000	5102 2/3 0.617956	-0.0014
	10000000 23466 2/3 0.62435	-0.00069
	100000000 108353 2/3 0.62936	-0.00032
	1000000	1285 0.55 0.51815	-0.018
	10000000	4331 0.55 0.51951	-0.012
	100000000 14738 0.55 0.52105	-0.0089
	1000000000 50577 0.55 0.52266	-0.0061
	10000000000 174773 0.55 0.52425	-0.0043
	Proof. The above linearization algorithm reduces the problem to linear algebra
	on a matrix of rank M (d). According to Lemma 4 (Appendix A), we get that

Note that we assume that d ≥ t.

Acknowledgements. We thank the reviewers of Eurocrypt 2012 for their very careful evaluation of our paper and their meaningful comments. We are also indebted to Sylvie Baudine for her comments on a preliminary version of this paper. The third author would like to thank the PolSys team for hosting him in the beginning of 2011, as well as Kristin Lauter and Daniele Micciancio since this paper was partially written in San Diego.

IAP program P6/26 BCRYPT. This work was partly supported by the Commission of the European Communities through the ICT program under contract ICT-2007-216676 (ECRYPT-II) and by the French ANR under the Computer Algebra and Cryptography (CAC) project (ANR-09-JCJCJ-0064-01).

A Proof of Lemma 2

In order to show Lemma 2, We use the following well-known technical result (a proof is given in [START_REF] Faugère | New Subexponential Algorithms for Factoring in SL(2, F n 2 )[END_REF] for instance). Lemma 3. Let n be an integer and let δ, 0 < δ < 1/2 be a number such that δn ∈ N. Finally, let ν := δ 1-δ . Then, the following inequalities hold:

We can now come back to the proof of Lemma 2.

Proof (of Lemma 2). We recall that 1α < 1 2 < α < 1, and n = n α , m = n (1-α) , t = m-1. The number of equations in the system generated in Algorithm

where d = n β for n big enough.

A sum of binomials can be bounded from above by the last binomial occurring in the sum. We then get:

When n is large enough, we have δ < 1 3 . This leads to ν ≤ 1 2 and 1 1-ν ≤ 2. We are now in position to apply Lemma 3.

Hence, we want to find d such that n

We search now for the equality and consider the logarithm of each side of (6):

The following result is useful to derive an asymptotical equivalent of this equality:

For n large enough, we have:

Taking γ = 0, we deduce log

As a consequence (7) yields:

So that 1 ≈ (α-β) log(n). Thus, β ≈ α-1/ log(n) and d β ≈ n α-1 log(n) = n α /2. This concludes the proof of Lemma 2.