Marie-Claude Portmann
email: portmann@loria.fr

Mouloua Zerouk
email: mouloua@loria.fr

A window time negotiation approach at the scheduling level inside supply chains

Keywords: scheduling, supply chain, just in time, negotiations, decomposition approach, "pert cost" algorithm, meta-heuristics

We consider a supply chain, which consists in a network of enterprises with independent decision centers. The finished products (or sub products) of the assembly enterprise are produced using components and/or sub products supplied by other enterprises or by external suppliers. We assume we are at the scheduling level and each enterprise builds its own schedules associated with its own production centers. As an operation can be performed only when the production center has received the necessary components and/or sub products, the schedules are dependent. This induces negotiations between decision centers of enterprises, which can be expressed in term of penalty functions associated with soft and hard release dates and due dates. At each negotiation point, hard release dates and due dates are considered as imperative constraints, while soft release dates and due dates define soft intervals and induce earliness and tardiness penalties. Both soft and hard constraints can be modified during the negotiation process. A global solution is searched by an iterative decomposition approach including alternatively bilateral negotiations of the soft and hard constraints between the production decision centers and just in time scheduling, minimizing the local total sum of penalties, built by approximation approaches. We assume that production centers are flow shop (linear production and/or assembly line). To solve each local just-in-time scheduling problem, we propose an approximation approach based on meta-heuristics, which explores the set of feasible and infeasible solutions, in which a solution is described by the job order on each machine (permutation of integers) and is evaluated using a "pert cost" algorithm. The infeasible solutions are evaluated by a lower bound of the number of non verified imperative constraints and the feasible solutions are evaluated by the minimal sum of penalties corresponding to the considered order. A semi-decentralized control is suggested to assume the negotiation convergence.

Introduction

We consider a supply chain, which consists in a network of enterprises with independent decision centers. In consequence, we assume each decision center can prepare its own schedules and negotiate release dates and due dates with their suppliers (upstream level of the supply chain) and their customers (downstream level of the supply chain). The decision system could have been totally decentralized, but in this case, there is no guarantee of convergence of the global negotiation process. In consequence, we choose a semi-decentralized decision system in which a common memory is shared and global rules are agreed at the network level in order to help the negotiations to converge (eventually to infeasible solutions of the initial global just in time scheduling problem, by authorization of unlimited increase of the finished product tardiness in the worst cases). We consider the following production system. The finished products (or sub products) of the assembly enterprise are produced using components and/or sub products supplied by other enterprises or by external suppliers. We assume we are at the scheduling level and each enterprise builds its own schedules associated with its own production centers. We assume the transport between production centers is taken into account only by minimal time lags. As an operation on a product can be performed only when an enterprise has received the needed components and sub products, the schedules are in fact dependent. This induces negotiations between enterprises, which can be expressed in term of soft and hard release dates and due dates. Hard release dates and due dates can be changed in order to increase or decrease the set of feasible solutions of each local partner. Soft release dates and due dates and their associated earliness and tardiness penalties can be real or only virtual delays and penalties in order to lead the whole network towards a win-win global solution. A global solution is searched by an iterative decomposition approach including alternatively bilateral negotiations of the soft and hard delays between the production centers and just in time scheduling minimizing the total sum of penalties inside each decision or production center. A semi-decentralized global negotiation approach will be suggested in section 2. We assume that production centers are either flow shop (assembly line) or even job shop (this complicates only the just in time scheduling algorithms). To solve each local just-in-time scheduling problem associated with a given set of hard and soft release dates and due dates and a given set of penalty data, we propose an approximation approach based on meta-heuristics, which explores the set of feasible and infeasible solutions, in which a solution is described by the job order on each machine (permutation of integers) and is evaluated, using a polynomial "pert cost" algorithm. Section 3 presents the just in time scheduling problem and how the timing problem can be solved by a "pert cost" algorithm. Some experiments concerning the just in time scheduling problem are presented at the end of section 3.

Global negotiation model

Global problem description

In the framework of decentralized or semi-decentralized decision system inside supply chain, several models have been proposed at the medium term level [START_REF] Dudek | Negotiation-based collaborative planning between supply chains partners[END_REF], [START_REF] Hurtubise | Planning tools for managing the supply chain[END_REF]). This proposition is in continuation of our own research works at the medium term planning [START_REF] Ouzizi | A model for co-operative planning within a virtual enterprise[END_REF]. At the medium term level, volumes associated with family of finished products, sub products and components are forecasted and production and transportation phases are planned in order to take into account the limited capacities of the resources. Negotiations can be organized between the independent partners in order to take into account their limited capacities and aggregated linear models can be used in order to compute each local planning. At the medium term level, hard and soft cumulated curves of component arrival or product delivery are negotiated between the supply chain partners and complementary resources can be bought. At the scheduling level, the forecasted medium term volumes on product family have been converted into firm orders of very precisely definite products. Each order is associated with one or several jobs (if lot sizing or splitting is used). Hence a job corresponds either to only one product or to an inseparable series of products, for which the detail schedule parameters are known: production centers to be visited, main sub products to be waited for (release date on the first machine of a given shop or further on another assembly machine), one or more production process inside each center (generalized job shop or hybrid flow shop). We assume here that external due dates have been promised to the final customers, which do not belong to the supply chain, and external release dates have been agreed with the external suppliers, most of them being imperative, due for example to transportation delays. We had to find the global schedule of the considered supply chain, while minimizing the total production cost and the total external penalties: win-win strategy inside the supply chain.

Global negotiation process

The ideas used at the medium term level, i.e. hard and soft cumulated curves of component arrival or product delivery with earliness and tardiness penalties negotiated between the supply chain partners, could be extended to the scheduling level replacing hard and soft cumulated curves by hard and soft operation windows. The main idea of the global negotiation process is also an extension of the process proposed in [START_REF] Portmann | Méthodes de décompositions spatiales et temporelles en ordonnancement de la production[END_REF] and [START_REF] Chu | A splitting up approach to simplify job-shop scheduling problems[END_REF]. In these papers, the criterion considered was the sum of weighted tardiness (a regular criterion) and only active schedules associated with hard release dates were considered; in consequence, only hard release dates, but hard and soft due dates were used in this iterative splitting approximation approach, whose aim was to build a global schedule for a given production center, obtained by scheduling optimally cells and modifying hard release dates and hard and soft local due dates. This splitting approach gave very efficient solution methods. The extended similar process is the following one illustrated by a small example on one product.

Figure 1: external constraints

In figure 1, soft and hard due dates associated with a given finished product order due to an external customer are denoted by d E and D E . The assembly of this finished product order needs two types of sub products used at two successive operations of the assembly line and sub products, made by two different supply chain partners, need themselves components from external suppliers, whose soft and hard release dates are respectively r A , R A , r C and R C . If there exists only one decision center for the whole supply chain, then we have got only one just in time scheduling problem to be solved, in which production costs must be minimized while verifying the external hard release dates and due dates (if possible) and minimizing the total sum of earliness and tardiness external penalties. While we consider a semi-decentralized control of independent decision centers, each decision center computes its own just in time scheduling and negotiates hard and soft release dates with the other partners of the supply chain with the aim of finding a global solution as good as possible, which will be a trade-off between the interests of the partners. To organize the negotiations, we need initial soft and hard window constraints between the whole set of partners. They will be called "virtual"; because they are not imposed by the global just in time scheduling problem and they will change during the global negotiation process.

A B C D E R A r A R C r C d E D E

Virtual sub product soft window

Virtual sub product soft window Virtual product soft windows

d B = r D d B = r D
First, similarly to [START_REF] Portmann | Méthodes de décompositions spatiales et temporelles en ordonnancement de la production[END_REF], we use the soft release dates and due dates of each sub product, of the final product and their processing times and we call a sharing function (for example, soft window lengths are proportional to the corresponding sum of processing times) in order to define "virtual" windows for each manufacturing part (see figure 2). It can be remarked that release and due dates are associated to any operation, whose predecessors and/or successors are in another decision centers (minimal time lags can be added in the model for transportation duration).

Hard initial windows are also obtained by a given rule. For example, they can be very large at the beginning of the global process (each operation is left shifted on its own hard release date or on its previous left shifted operations for the hard release dates and each operation is right shifted on its own hard due date or on its following right shifted operations for the hard due dates) or the difference between the hard and soft window lengths of each sub product and/or product is computed by a function similar to the one use for the soft constraints (proportionality rule). Virtual penalties are attributed equally without any further information or individually increased when, for example, it is known that some machine of the production system is overloaded and upstream and downstream operations must respect as tight as possible the virtual given windows.

After this hard and soft windows attribution, the schedules of the various production centers become independent and an exact or approximation scheduling tool can be used in order to minimize the earliness and tardiness penalty sum while respecting the hard windows. Each schedule can be computed locally by the decision centers in real life applications. This initial schedule computation finishes the first step of the global process.

An analysis is made at the end of each step of this splitting approach. Case 1: If the external soft release dates and due dates are respected and any precedence constraint between two different production centers is respected, i.e. when an operation is late relatively to hard or soft "virtual" due dates, then we have the chance that the following operations are at least as late as it and when an operation is early relatively to its hard or soft "virtual" release date then the previous operations are sufficiently early for the global schedule to be feasible. It is the perfect case and the negotiation process can stop without any external negotiations. Case 2: If external hard release dates and due dates are respected, but not any soft external windows, and nevertheless the schedules of the various production centers are compatible, then two ways can be followed before the next step, either to negotiate new external soft release dates and/or due dates or to increase the penalty associated with the external earliness and tardiness hoping to get a better set of feasible local schedules. Case 3: If external hard windows are respected, but the set of schedules induce incompatibility between production centers, then the internal hard and soft windows must be modified in order to converge to a set of feasible schedules. In a semi-decentralized decision system approach, the order in which the windows are negotiated can be organized, going from supplier to customer (direct succession of schedules and windows modification) or from customer to supplier (retrograde succession of schedules and windows modification, knowing that the minimal hard release date must remain always positive). This organization is a way to assume the convergence of the global process (assuming that case 4 can enlarge the external hard windows to increase the set of feasible solutions). Case 4: If external hard windows are not respected, it is a very delicate situation. The explanation could be that the whole global problem has no feasible solution, but it is impossible to verify if the problem is too big to use an exact algorithm to solve the centralized global just in time scheduling problem. The only way to continue is to relax slightly the hard final due dates and to begin a new negotiation process between the partners. The convergence of the process can be assumed by the organization of the authorized negotiations (no loop in the negotiation process) and by the last assumption that hard final due dates can be relaxed. The efficiency of this global negotiation approach can only be verified by experimentation using simulation for modeling the human behavior inside the decision process.

Just in time scheduling problems

To complete this window negotiation process, we need a tool for solving the scheduling problem with hard and soft window inside each production center. We simplify here the scheduling problem assuming there are no production costs (such as inventory costs or overtime hours or sub contracting …) to be considered and we have only to minimize the sum of earliness and tardiness penalties.

If the local shop we consider is a flow shop and hard and soft release dates are associated to the first job operations, while hard and soft due dates are associated to the last job operations, then our problem can be denoted by F/ are respectively the hard and soft release dates and due dates of the first and last operations of each job and α r,i , β r,i , α d,i , β d,i are the earliness (α) or tardiness (β) penalty associated with the release date of the first operation ( r ) and the due date of the last operation ( d ). When we assume that release dates and/or due dates can be associated to any operations of the jobs, we can use a more general α/β/γ notation: J / tmin j , t* j , tmax j / Σα j E j +Σβ j T j , where t* j is the ideal position of the starting time of the operation j, tmin j is its minimal authorized value (hard constraint for earliness) and tmax j is its maximal authorized value (hard constraint for tardiness). We propose to use another approximation decomposition approach to solve this local problem. If we assume the operation sequence is given on each machine (without circuit for the job shop problem), then we have only to compute the optimal start time of each operation (i.e. where are put the idle times, because our criterion is not regular) and we use a meta-heuristic approach to explore the set of operation sequences on the machine in order to find a feasible solution as good as possible. We now present the timing scheduling problem when the sequence is given (by a chromosome) and the main feature of the meta-heuristic afterwards.

Timing scheduling problem (i.e. chromosome evaluation)

We have to solve the F seq /

i i i i d , d , r , r
/Σα r,i E r i + Σβ r,I T r i +Σα d,i E d i +Σβ d,I T d i scheduling problem or the J seq / tmin j , t* j , tmax j / Σα j E j +Σβ j T j , where "seq" means the sequences on the machine are given by a chromosome. For job shop problem, some chromosomes can be unacceptable whatever could be the hard windows because of circuit in the precedence graph, while this cannot happen for flow shop or single machine problem. We assume the meta-heuristic used is able to build only chromosome with no circuit in the precedence graph (we have not this problem in our first experiments because we consider only flow shops). When considering only the hard windows ( i i d , r

) or [tmin j , tmax j ], for fixed sequence on the resources, our problem is simply a "PERT" problem, which could be polynomially solved by computing the left shifted schedule (λ k,i ) or (λ j ) and the right shifted schedule (λ' k,i ) or (λ' j ) and verifying that for each operation k of each job i (denoted j) we have λ j ≤λ' j . We use the number of times this inequality is not verified to measure the infeasibility of the chromosomes relatively to the hard windows, but we keep this infeasible chromosome with very bad fitness as current solution of simulated annealing or inside the genetic algorithm population.

We now assume a chromosome is totally feasible (no circuit due to the sequences on the machine and there exists at least one solution verifying the hard windows) and we want to find each operation starting time of a feasible schedule in order to minimize the total sum of earliness and tardiness penalty. This problem can be modeled as a "pert cost" scheduling problem as follows. s is the source and t the sink of the "pert cost" extended graph CG, where each operation, assumed to be non preemptive can be represented by its beginning event. The duration of each operation is added to the value of the arc issued from this event. The operation sequence of the job and on the machine is represented by usual precedence constraints between the operations. We denote by H the greatest hard due date (or the greatest tmax j ), H is the length of the considered horizon. The figure 3 gives the sub graph representing earliness or tardiness penalty for any operation. OP j is the beginning of the operation j, t* j is its ideal position. If x i is equal to 0, then the operation j is ideally placed. If x i is positive, then the event is early and if it is negative then the event is late. We considered that the arc issued from s is compressed and induce a cost α j .x j when x j is positive and that the arc arriving in t is compressed and induce a cost -β j .x j when x j is negative. It is important to remark than the length of the path from s to t through BS j , OP j and ES j is exactly equal to H, the length of the horizon and the values x j will be chosen so that the maximum path in the complete precedence graph will also be equal to H, and in this case, any operation j is placed on the time axis by knowing the value of x j . t* j -x j must remain greater than tmin j and smaller than tmax j .

Figure 3: precedence constraints linked to earliness and tardiness penalties

When points of the graph CG are in series with the same set of successors or predecessors, we can delete the beginning of the point series and add the values of the delete arcs to the last arc kept. The data of a particular instance corresponding to a flow shop with two machines and earliness and tardiness cost relatively to the release dates for the first machine and the due dates for the second machine are given in table 1 and the simplified corresponding "pert cost" extended graph is given in the figure 4.

Table 1: precedence constraints linked to earliness and tardiness penalties

It is a particular "pert cost" problem in which you can only pay to decrease the duration of the arcs adjacent to the source s or to the sink t. You can solve it polynomially using a dual approach. But due to the particular shape of the graph, the primal algorithm (using alternatively "pert" algorithm and search of minimal cut on the critical sub graph), when using the maximal compression associated to each cut, is probably polynomial (we do not prove it formally until now). Any way, the minimal sum of earliness and tardiness penalty can be computed for each chromosome by using a "pert cost" algorithm and this sum is used to compute the fitness of each chromosome.

BS

Meta-heuristic approaches

We already used and published papers using meta-heuristic and more specifically genetic algorithms for which we analyze specifically the performance of the permutation cross-over operators for various criteria [START_REF] Djerid | How to Keep Good Schemata Using Cross-over Operators for Permutation Problems[END_REF], [START_REF] Portmann | Performances' study on crossover operators keeping good schemata for some scheduling problems[END_REF]. Nevertheless a rapid experimentation using the simple genetic algorithm (SGA) of [START_REF] Goldberg | Genetic Algorithms in Search, Optimization and Machine Learning[END_REF] provides us with deceiving results and our genetic algorithm needs to be improved. On another way, simulated annealing is very rapid to test and only to show meta-heuristic could be a good approach for our problem we use it also in the following results. The experiments were made on generated flow shop examples with 10 jobs and 5 machines. For each instance, three algorithms were used: random generation of a sequence (Random), simulated annealing (SA) and simple algorithm genetic (GA). To see the rapidity to find solution as good as possible, we compute the function BEST(Algo, k) which gives the value of the best found feasible solution after having evaluated k sequences. To eliminate the fact that hazard could be more or less favorable, we execute each algorithm 5 times on the same instance and we keep the mean value, which provides us with MEAN-BEST(Algo,k) and we present here in the figure 5 the average on four instances. Only 300 feasible solutions analyzed by "pert cost", where generated for each execution of the algorithm. Bad results obtained by the genetic algorithm can be due to the fact that the initial population is randomly built and it is possible, with only 300 sequences studied among 10!, SA or GA have not yet converged. We will continue our experimentation with longer computational time and with improvement of our very simple current meta-heuristic scheme.

First experimentation for the just in time scheduling problem

Conclusion

We propose here a new decomposition approach to organize the negotiation between various independent partners of a supply chain at the scheduling level using hard and soft windows as entity of negotiation. This decomposition can also be used as an approximation approach for building a just in time solution in presence of only one decision center and several production centers, as an extension of the method proposed in [START_REF] Chu | A splitting up approach to simplify job-shop scheduling problems[END_REF]. The interest of this last proposition is that it can be compared with other methods using experimentations, while the model with several decision centers and human decision models must be validated using simulation. Another important problem will be to convince the industrial decision makers and managers of the interest of the proposed models.

Figure 2

 2 Figure 2: initial soft window constraints

Figure 4 :

 4 Figure 4: simplified extended "pert cost" for a flow shop problem

Figure 5 :

 5 Figure 5: comparison of approximation approach


Aurore Guillevic
email: aurore.guillevic@inria.fr

FASTER INDIVIDUAL DISCRETE LOGARITHMS IN FINITE FIELDS OF COMPOSITE EXTENSION DEGREE

Keywords: Subject Classification. Primary 11T71: Cryptography Finite field, discrete logarithm, number field sieve, function field sieve, individual logarithm Waterloo alg. Algorithm

Computing discrete logarithms in finite fields is a main concern in cryptography. The best algorithms in large and medium characteristic fields (e.g., GF(p 2 ), GF(p 12 )) are the Number Field Sieve and its variants (special, high-degree, tower). The best algorithms in small characteristic finite fields (e.g., GF(3 6•509 )) are the Function Field Sieve, Joux's algorithm, and the quasipolynomial-time algorithm. The last step of this family of algorithms is the individual logarithm computation. It computes a smooth decomposition of a given target in two phases: an initial splitting, then a descent tree. While new improvements have been made to reduce the complexity of the dominating relation collection and linear algebra steps, resulting in a smaller factor basis (database of known logarithms of small elements), the last step remains at the same level of difficulty. Indeed, we have to find a smooth decomposition of a typically large element in the finite field. This work improves the initial splitting phase and applies to any nonprime finite field. It is very efficient when the extension degree is composite. It exploits the proper subfields, resulting in a much more smooth decomposition of the target. This leads to a new trade-off between the initial splitting step and the descent step in small characteristic. Moreover it reduces the width and the height of the subsequent descent tree.

Introduction

This work is interested in improving the last step of discrete logarithm (DL) computations in nonprime finite fields. The discrete logarithm instances that we target come from Diffie-Hellman (DH) [START_REF] Diffie | New directions in cryptography[END_REF] key-exchange, or from pairing-based cryptography. In the latter case, the security relies on the hardness of computing discrete logarithms in two groups: the group of points of a particular elliptic curve defined over a finite field, and a small extension of this finite field (in most of the cases of degree 2, 3, 4, 6, or 12).

The finite fields fall in three groups: small, medium and large characteristic finite fields, corresponding to the respective size of the characteristic p compared to the total size Q = p n of the finite field. This is formalized with the L notation:

(1.1) L Q [α, c] = e (c+o(1))(log Q) α (log log Q) 1-α , where Q = p n , α ∈ [0, 1], c = 0.
Small, medium and large characteristic correspond to α < 1/3, 1/3 < α < 2/3, and α > 2/3 respectively. The boundary cases are α = 1/3 and α = 2/3. In large characteristic, that is p = L Q [α, c] where α > 2/3, the Number Field Sieve (NFS) [START_REF] Daniel | Discrete logarithms in GF(p) using the number field sieve[END_REF][START_REF] Schirokauer | Discrete logarithms and local units[END_REF][START_REF] Joux | The function field sieve is quite special[END_REF] provides the best expected running time: in L Q [1/3, (64/9) 1/3 ≈ 1.923] and was used in the latest record computations in a 768-bit prime field [START_REF] Kleinjung | Computation of a 768-bit prime field discrete logarithm[END_REF]. Its special variant in expected running time L Q [1/3, (32/9) 1/3 ≈ 1.526] was used to break a 1024-bit trapdoored prime field [START_REF] Fried | A kilobit hidden SNFS discrete logarithm computation[END_REF]. In 2015 and 2016, the Tower-NFS construction of Schirokauer was revisited for prime fields [START_REF] Barbulescu | The tower number field sieve[END_REF], then Kim, Barbulescu and Jeong improved it for nonprime finite fields F p n where the extension degree n is composite [START_REF] Kim | Extended tower number field sieve: A new complexity for the medium prime case[END_REF][START_REF] Kim | Extended tower number field sieve with application to finite fields of arbitrary composite extension degree[END_REF], and used the name Extended TNFS algorithm. To avoid a confusion due to the profusion of names denoting variants of the same algorithm, in this paper we will use TNFS as a generic term to denote the family of all the variants of NFS that use a tower of number fields. Small characteristic means p = L Q [α, c] where α < 1/3. The first L[1/3] algorithm was proposed by Coppersmith, and generalized as the Function Field Sieve [START_REF] Adleman | The function field sieve[END_REF][START_REF] Adleman | Function field sieve method for discrete logarithms over finite fields[END_REF].

The NFS and FFS algorithms are made of four phases: polynomial selection (two polynomials are chosen), relation collection where relations between small elements are obtained, linear algebra (computing the kernel of a huge sparse matrix over an auxiliary large prime finite field), and individual discrete logarithm computation. In this work, we improve this last step. All the improvements of NFS, FFS, and related variants since the 90's decrease the size of the factor basis, that is, the database of known discrete logarithms of small elements obtained after the linear algebra step, small meaning an element represented by a polynomial of small degree (FFS), resp., an element whose pseudonorm is small (NFS). The effort required in the individual discrete logarithm step increases: one needs to find a decomposition of a given target into small elements, to be able to express its discrete logarithm in terms of already known logarithms of elements in the factor basis, while the factor basis has decreased at each major improvement. In characteristic 2 and 3 where the extension degree is composite, obtaining the discrete logarithms of the factor basis elements can be done in polynomial time. The individual discrete logarithm is the most costly part, in quasi-polynomial -time in the most favorable cases [START_REF] Barbulescu | A heuristic quasipolynomial algorithm for discrete logarithm in finite fields of small characteristic[END_REF][START_REF] Granger | Breaking '128-bit secure' supersingular binary curves -(or how to solve discrete logarithms in F 2 4•1223 and F 2 12•367[END_REF]. In practice, the record computations [START_REF] Göloglu | Discrete logarithms in GF(2 6 120), Number Theory list[END_REF][START_REF] Joux | [END_REF][START_REF] Adj | Computing discrete logarithms in F 3 6•137 and F 3 6•163 using Magma, Arithmetic of Finite Fields[END_REF][START_REF] Adj | Weakness of F 6 6•1429 and F 2 4•3041 for discrete logarithm cryptography[END_REF][START_REF] Joux | Discrete logarithm record in characteristic 3[END_REF][START_REF] Adj | 6 * 509[END_REF] implement hybrid algorithms made of Joux's L[1/4] algorithm [START_REF]A new index calculus algorithm with complexity L(1/4 + o(1)) in small characteristic[END_REF], and the individual discrete logarithm is computed with a continued fraction descent, then a classical descent, a QPA descent, and a Gröbner basis descent, or a powers-of-two descent algorithm (a.k.a. zig-zag descent) [START_REF] Granger | On the powers of 2[END_REF][START_REF]Discrete logarithms in the Jacobian of a genus 2 supersingular curve over GF(2 367 )[END_REF][START_REF] Granger | Discrete logarithms in GF(2[END_REF].

The heart of this paper relies on the following two observations. Firstly, to speed up the individual discrete logarithm phase, we start by speeding up the initial splitting step, and for that we compute a representation of a preimage of the given target of smaller degree, and/or whose coefficients are smaller. It will improve its smoothness probability. Secondly, to compute this preimage of smaller degree, we exploit the proper subfields of the finite field GF(p n ), and intensively use this key-ingredient: since we are computing discrete logarithms modulo (a prime divisor of) Φ n (p), we can freely multiply or divide the target by any element in a proper subfield without affecting its discrete logarithm modulo Φ n (p).

Organization of the paper. The background needed is presented as preliminaries in Section 2. We present our generic strategy to lower the degree of the polynomial representing a given element in GF(p n ) in Section 3. We apply it to characteristic two and three in Section 4. Preliminaries before the large characteristic case are given in Section 5. We apply our technique to medium and large characteristic finite fields, that is the NFS case and its tower variant in Section 6, and provide examples of cryptographic size in Section 7. Finally in Section 8 we present a more advanced strategy, to exploit several subfields at a time, and we apply it to F p 6 .

Preliminaries

2.1. Setting. In this paper, we are interested in nonprime finite fields GF(p n ), n > 1. To keep the same notation between small, medium, and large characteristic finite fields, we assume that the field F p n is defined by an extension of degree n 2 above an extension of degree n 1 , that is, F (p n 1 ) n 2 , and n = n 1 n 2 . The elements are of the form T = n2-1 i=0 n1-1 j=0 a ij y j x i , where the coefficients a ij are in F p , the coefficients

a i = n1-1 j=0 a ij y j are in F p n 1 = F p [y]/(h(y)), and F (p n 1 ) n 2 = F p n 1 [x]/(ψ(x))
, where h is a monic irreducible polynomial of F p [y] of degree n 1 and ψ is a monic irreducible polynomial of F p n 1 [x] of degree n 2 . In other words, T is represented as a polynomial of degree n 2 -1 in the variable x, and has coefficients a i ∈ F p n 1 . For the FFS and NFS algorithms, n 1 = 1 and n 2 = n; for finite fields from pairing constructions, n 2 > 1 is a strict divisor of n, and for the original version of TNFS, n 1 = n and n 2 = 1. Definition 2.1 (Smoothness). Let B be a positive integer. A polynomial is said to be B-smooth w.r.t. its degree if all its irreducible factors have a degree smaller than B. An integer is said to be B-smooth if all its prime divisors are less than B. An ideal in a number field is said to be B-smooth if it factors into prime ideals whose norms are bounded by B.

Definition 2.2 (Preimage). The preimage of an element a

= n2-1 i=0 n1-1 j=0 a ij y j x i ∈ F (p n 1
) n 2 will be, for the NFS and TNFS algorithms, the bivariate polynomial

n2-1 i=0 n1-1 j=0 a ′ ij y j x i ∈ Z[x, y],
where each coefficient a ′ ij is a lift in Z of the coefficient a ij in F p . It is a preimage for the reduction modulo (p, h, ψ), that we denote by ρ : Z[x, y] → F (p n 1 ) n 2 . In small characteristic, the preimage of a is a univariate polynomial in F p n 1 [x]. It is a preimage for the reduction modulo ψ, that we also denote by ρ :

F p n 1 [x] → F (p n 1 ) n 2 . Definition 2.3 (Pseudonorm). The integral pseudonorm w.r.t. a number field Q[x]/(f (x)) (f monic) of a polynomial T = deg f -1 i=0 a i x i of integer coefficients a i is computed as Res x (T (x), f (x)).
Since there is no chance for a preimage of a target T 0 to be B-smooth, the individual discrete logarithm is done in two steps: an initial splitting of the target, 1 and then a descent phase. 2 The initial splitting is an iterative process that tries many targets g t T 0 ∈ F * p n , where t is a known exponent (taken uniformly at random), until a B 1 -smooth decomposition of the preimage is found. Here smooth stands for a factorization into irreducible polynomials of F p n 1 [x] of degree at most B 1 in the small characteristic setting, resp., a pseudonorm that factors as an integer into a product of primes smaller than B 1 in the NFS (and TNFS) settings.

1 also called boot or smoothing step in large characteristic finite fields 2 in order to make no confusion with the mathematical descent, which is not involved in this process, we mention that in this step, the norm (with NFS) or the degree (with FFS) of the preimage decreases.

The second phase starts a recursive process for each element less than B 1 but greater than B 0 obtained after the initial splitting phase. Each of these mediumsized elements are processed until a complete decomposition over the factor basis is obtained. Each element obtained from the initial splitting is at the root of its descent tree. One finds a relation involving the original one and other ones whose degree, resp., pseudonorm, is strictly smaller than the degree, resp., pseudonorm, of the initial element at the root. These smaller elements form the new leaves of the descent tree. For each leaf, the process is repeated until all the leaves are elements in the factor basis. The discrete logarithm of an element output by the initial splitting can be computed by a tree traversal. This strategy is considered in [19, §6], [48, §7], [38, §3.5], [18, §4].

In small characteristic, the initial splitting step is known as the Waterloo 3 algorithm [START_REF] Blake | Computing logarithms in finite fields of characteristic two[END_REF][START_REF] Blake | Computing logarithms in GF(2 n )[END_REF]. It outputs T = U (x)/V (x) mod I(x), and U, V are two polynomials of degree ⌊(n 2 -1)/2⌋. It uses an Extended GCD computation. For prime fields, the continued fraction algorithm was already used with the Quadratic Sieve and Coppersmith-Odlyzko-Schroeppel algorithm. It expresses an integer N modulo p as a fraction N ≡ u/v mod p, and the numerator and denominator are of size about the square root of p. The generalization of this technique was used in [START_REF] Joux | The number field sieve in the medium prime case[END_REF]. As for the Waterloo algorithm, this technique provides a very good practical speed up but does not improve the asymptotic complexity.

this subfield tool was highlighted in [START_REF] Guillevic | Computing individual discrete logarithms faster in GF(p n ) with the NFS-DL algorithm[END_REF]; we will intensively use it.

Lemma 2.4 ([34, Lemma 1]). Let T ∈ F * p n , and let deg T < n. Let ℓ be a nontrivial prime divisor of Φ n (p). Let T ′ = u • T with u in a proper subfield of F p n . Then (2.1)
log T ′ ≡ log T mod Φ n (p) and in particular log T ′ ≡ log T mod ℓ .

3.

The heart of our strategy: representing elements in the cyclotomic subgroup of a nonprime finite field with less coefficients

In the FFS setting, n 1 = 1 and usually n 2 is prime and our technique cannot be helpful, but if n is not prime, our algorithm applies, and moreover in favorable cases Joux's L[1/4] algorithm and its variants can be used and our technique can provide a further notable speed-up in the descent. For the implementations in small characteristic, the factor basis is made of the irreducible polynomials of F p n 1 [x] of very small degree, e.g., of degrees 1, 2, 3, and 4 in [START_REF] Adj | Computing discrete logarithms in cryptographically-interesting characteristic-three finite fields[END_REF]. Our aim is to improve the smoothness probability of a preimage P ∈ F p n 1 [x] of a given target T ∈ F (p n 1 ) n 2 and for that we want to reduce the degree in x of the preimage P (as a lift of T in F p n 1 [x], P has degree at most n 2 -1 in x), while keeping the property log(ρ(P )) = log T mod ℓ , where ρ :

F p n 1 [x] → F (p n 1 ) n 2 is the reduction modulo ψ.
Let d denote the largest proper divisor of n, 1 < d < n (d might sometimes be equal to n 2 in the QPA setting). We will compute P in F p n 1 [x] of degree at most n 2d/n 1 in x (and coefficients in F p n 1 ) such that (3.1) P = uT (mod ψ), where u p d -1 = 1 .

3 the name comes from the authors' affiliation: the University of Waterloo, ON, Canada.

It means that we will cancel the d/n 1 -1 higher coefficients (in

F p n 1 ) of a preimage of T in F p n 1 [x].
There are two strategies: either handle coefficients in F p or in F p gcd(d,n 1 ) . We will consider the latter case. Let d ′ = d/ gcd(d, n 1 ) to simplify the notation, and let [1, U, . . . , U d ′ -1 ] be a polynomial basis of F p d ′ . Every product P = U i T satisfies (3.1). Define the d ′ × n 2 matrix L whose rows are made of the coefficients (in F p n 1 ) of U i T for 0 ≤ i ≤ d ′ -1:

L d ′ ×n2 =      T U T . . . U d ′ -1 T      ∈ M d ′ ,n (F p n 1 ) .
Then we compute a row-echelon form of this matrix by performing only F p gcd(n 1 ,d)linear operations over the rows, so that each row of the echeloned matrix is a F p gcd(n 1 ,d) -linear combination of the initial rows, that can be expressed as

P = d ′ -1 i=0 λ i U i T = uT, where λ i ∈ F p gcd(n 1 ,d) , U i ∈ F p d/ gcd(n 1 ,d)
so that P = uT with u p d -1 = 1. Assuming that the matrix is lower-triangular (the other option being an upper-triangular matrix), we take the first row of the matrix as the coefficients of a polynomial in F p n 1 [x] of degree at most4 n 2d/n 1 . This is formalized in Algorithm 1. We obtain the following Theorem 3.1.

Algorithm 1: Computing a representation by a polynomial of smaller degree Input: Finite field F p n represented as a tower We can now directly apply Algorithm 1 to improve the initial splitting algorithm in practice.

F (p n 1 ) n 2 = F p n 1 [x]/(ψ(x)) (one may have n 1 = 1), a proper divisor d of n (d | n, 1 < d < n), T ∈ F p n Output: P ∈ F p n 1 [x] a polynomial of degree ≤ n 2 -d/n 1 satisfying P mod ψ = uT , where u ∈ F p d d ′ = d/ gcd(n 1 , d) Compute a polynomial basis (1, U, U 2 , . . . , U d ′ -1 ) of the subfield F p d ′ Define L =      T U T . . . U d ′ -1 T      a d ′ × n 2 matrix of coefficients in F p n 1 M ← RowEchelonForm(L)

Application to small characteristic finite fields, and cryptographic-size examples

In all the examples of small characteristic finite fields from pairings, n is not prime, for instance n = 6 • 509. The notation in [START_REF]Weakness of F 3 6509 for discrete logarithm cryptography[END_REF] was n = lk, with the property p l ≈ k. With our notation, n 1 = l and n 2 = k. 4.1. Algorithm. We directly use Algorithm 1 as a subroutine of Algorithm 2. Then to improve it in practice, we list valuable modifications.

Remark 4.1. As was pointed out to us by F. Rodríguez-Henríquez [START_REF] Rodríguez-Henríquez | Another initial splitting in small characteristic finite fields[END_REF][START_REF] Adj | Computing discrete logarithms in cryptographically-interesting characteristic-three finite fields[END_REF], the elements of the form x i R(x) where R itself is of degree ≤ n 2d/n 1 are evenly interesting because the discrete logarithm of x i can be deduced from the discrete logarithm of x, which is known after linear algebra.

So we can increase the number of elements tested for B 1 -smoothness for each exponent t by a factor d ′ almost for free in the following way. We again run a Gaussian elimination algorithm on the matrix M but in the reverse side, for instance from row one to row d ′ and left to right if it was done from row d ′ to row one and right to left the first time. The matrix is in row-echelon form on the left-hand side and on the right-hand side (the upper right and lower left corners are filled with zeros). We obtain a matrix N of the form .

N =       * . . . *
The i-th row represents a polynomial P ′ i = x ei P i , where P i is of degree at most n 2d/n 1 , and e i ≈ (i -1) gcd(n 1 , d)/n 1 . Since x is in the factor basis (by construction, like all the degree one polynomials), its logarithm is known at this point (after the relation collection and linear algebra steps), hence the logarithm of any power x ei is known. It remains to compute the discrete logarithm of P i .

Algorithm 2: Initial splitting in small characteristic with the subfield technique Input: Finite field F p n of small characteristic (e.g., p = 2, 3), with a tower representation

F p n = F (p n 1 ) n 2 = F p n 1 [x]/(I(x)) (one may have n 1 = 1)
, generator g (of the order ℓ subgroup of the cyclotomic subgroup of

F p n ), target T 0 ∈ F (p n 1 ) n 2 , smoothness bound B 1 Output: t, P ∈ F p n 1 [x] a polynomial of degree ≤ n 2 -d/n 1 such that vlog g ρ(P ) = t + vlog g T 0 mod ℓ, and P (x) is B 1 -smooth (w.r.t. its degree in x) d ← the largest divisor of n, 1 < d < n d ′ ← d/ gcd(d, n 1 ) Compute U (x) ∈ F (p n 1 ) n 2 s.t. (1, U, U 2 , . . . , U d ′ -1 ) is a polynomial basis of the subfield F p d ′ repeat take t ∈ {1, . . . , ℓ -1} at random T ← g t T 0 in F (p n 1 ) n 2 Define L =      T U T . . . U d ′ -1 T      a d ′ × n 2 matrix of coefficients in F p n 1 M ← RowEchelonForm(L) (with F p gcd(d,n 1 ) -linear Gaussian elimination) P (x) ← the polynomial of lowest degree made of the first row of L until P (x) is B 1 -smooth return t, P (x)
In practice there are some technicalities: in the second Gaussian elimination, if the leading coefficient is zero, then two rows are swapped, and it cancels the previous Gaussian elimination (computed at the other end of the matrix) for that row. We end up with a matrix which is in row-echelon form on the right and almost row-echelon form on the left (or vice-versa). Since each set of subsequent n 1 / gcd(n 1 , d) rows produces polynomials of the same degree, swapping two rows from the same set will not change the degree in x of the polynomial. In average (this is what we observed in our experiments for F 3 6•509 and F 3 5•479 ), some rare polynomials will have a degree in x increased by one or two. This second Gaussian elimination increases the number of tests by a factor d ′ at a very cheap cost, since in fact it allows one to share the cost of computing the U i T and the two Gaussian eliminations over d ′ tests. Remark 4.2. If gcd(d, n 1 ) > 1 we can increase the number of rows by a small factor. We perform linear combinations of n 1 / gcd(d, n 1 ) subsequent rows (all giving a polynomial of same degree):

0≤j≤n1/ gcd(d,n1) µ j r j where µ ∈ F p gcd(d,n 1 ) , and it will result in new rows and new polynomials of same degree.

Remark 4.3. Other improvements are possible [START_REF] Rodríguez-Henríquez | Another initial splitting in small characteristic finite fields[END_REF][START_REF] Adj | Computing discrete logarithms in cryptographically-interesting characteristic-three finite fields[END_REF], for instance computing F p gcd(n 1 ,d) -linear combinations over a small number of rows corresponding to polynomials of almost the same degree. The resulting polynomial will have degree increased by one or two, which does not significantly affect its B 1 -smoothness probability in practice for cryptographic sizes. This technique allows one to produce many more candidates, at a very cheap cost of linear operations in

F p n 1 [x].
(d | n, 1 < d < n), and let d ′ = d/ gcd(d, n 1 ). Since d ′ | d | n = n 1 n 2 and gcd(d ′ , n 1 ) = 1, then d ′ | n 2 and d ′ ≤ n 2 . The computation of all the U i T of the matrix L costs at most d ′ n 2 2 multiplications in F p n 1 , since a schoolbook multiplication in F (p n 1 ) n 2 costs n 2 2 multiplications in F p n 1 . There are d ′ such multiplications. The complexity of a reduced row-echelon form computation of a (d ′ × n 2 )-matrix, d ′ ≤ n 2 , is less than O(d ′2 n 2 ) multiplications in F p n 1 [23, §13.4.2].
To simplify, we consider that the computation of the matrix L and of two Gaussian eliminations is done in time at most O(d ′ n 2

2 ). This cost is shared over d ′ polynomials P i to be tested for B 1 -smoothness. In this way, the complexity of computing a preimage P with our initial splitting algorithm is the same as in the Waterloo algorithm:

O(n 2
2 ), and moreover the smoothness probabilities are much higher for the targeted cryptographic cases coming from supersingular pairing-friendly curves. We also replace two B 1 -smoothness tests by only one, and that might save some time in practice (this saving disappears in the O notation). We present the theoretical costs in Tables 1 and2 from [START_REF] Flajolet | The complete analysis of a polynomial factorization algorithm over finite fields[END_REF]. XGCD stands for extended Euclidean algorithm, SQF stands for SQuare-free Factorization, DDF stands for Distinct Degree Factorization, and EDF stands for Equal Degree Factorization. All the polynomials to be factored are of degree smaller than n 2 ; we take n 2 as an upper bound to get the costs of Table 2.

Table 1. Costs for the initial splitting step. The preimage obtained with Algorithm 2 has degree d P ≤ n 2d/n 1 . The Waterloo algorithm [START_REF] Blake | Computing logarithms in finite fields of characteristic two[END_REF][START_REF] Blake | Computing logarithms in GF(2 n )[END_REF] produces two polynomials of degree

d P = ⌊n 2 /2⌋. Factorization cost Square-free (SQF) O(d 2 P ) Distinct degree (DDF) O(d 3 P log p n1 ) Equal degree (EDF) O(d 2 P log p n1 )
Table 2. Cost in multiplications in F p n 1 to compute one preimage to be tested for smoothness, in the initial splitting step.

Computation XGCD(T, To start, we recall some results on the smoothness probability of a polynomial of given degree.

I) matrix [U i T ] 0≤i≤d ′ -
Definition 4.4. Let N q (b; d) denote the number of monic polynomials over F q of degree d which are b-smooth. Let N q (b; d 1 , d 2 ) denote the number of coprime pairs of monic polynomials over F q of degrees d 1 and d 2 , respectively, which are b-smooth.

Let Pr q (b; d) denote the probability of a monic polynomial over F q of degree d to be b-smooth. Let Pr q (b; d 1 , d 2 ) denote the probability of two coprime monic polynomials over F q of degrees d 1 and d 2 to be both b-smooth. 2 (log n 2 ) 1-α b , and the degree of the polynomial to be tested for smoothness d = an 2 , where a ∈]0, 1[ and n 2 = log Q/ log p n1 , one obtains

Pr p n 1 (b, d) = L Q [1 -α b , -(1 -α b )a/γ] , where Q = p n1n2 . Theorem 4.5 ([22, Theorem 1]). Let δ > 0 be given. Then we have, uniformly for b, d 1 , d 2 → ∞ with d δ 1 ≤ b ≤ d 1-δ 1 and d δ 2 ≤ b ≤ d 1-δ 2 , N q (b; d 1 , d 2 ) ∼ 1 - 1 q N q (b; d 1 )N q (b; d 2 ) .
b, d 1 , d 2 → ∞ with d δ 1 ≤ b ≤ d 1-δ 1 and d δ 2 ≤ b ≤ d 1-δ 2 , Pr q (b; d 1 , d 2 ) ∼ 1 - 1 q Pr q (b; d 1 ) Pr q (b; d 2 ) .
We can now compare the Waterloo algorithm with this work. Assuming that B 1 = log p n 1 L p n [2/3, γ] for a certain γ, then the probability of a polynomial of degree an 2 , 0 < a < n 2 , to be B 1 -smooth is L p n [1/3, -a/(3γ)]. In the Waterloo algorithm, two polynomials of degree n 2 /2 should be B 1 -smooth at the same time, and the expected running time to find such a pair is

L p n [1/3, 1/(3γ)] (the square of L p n [1/3, 1(6γ)]). In our algorithm, a polynomial of degree ⌊n 2 -d/n 1 ⌋ = ⌊n 2 (1 -d/n)⌋ is tested for B 1 -smoothness, so finding a good one requires (4.2) L p n [1/3, a/(3γ)] tests, where a ≈ 1 -d/n ,
which is always faster than the Waterloo algorithm, for which a = 1. When n is even (this is always the case for finite fields of supersingular pairing-friendly curves), one can choose d = n/2, hence a = 1/2 and our algorithm has running time the square root of the running time of the Waterloo algorithm.

4.3.

Improving the record computation in GF(3 6•509 ). Adj, Menezes, Oliveira, and Rodríguez-Henríquez estimated in [START_REF]Weakness of F 3 6509 for discrete logarithm cryptography[END_REF] the cost to compute discrete logarithms in the 4841-bit finite field GF(3 6•509 ) and announced their record computation in July 2016 [START_REF] Adj | 6 * 509[END_REF]. The details of the computations are available in Adj's PhD thesis [START_REF]Logaritmo discreto en campos finitos de característica pequeña: atacando la criptogrfía basada en emparejamientos de tipo 1[END_REF] and the details for initial splitting and descent can be found in [START_REF] Andrés Canales-Martínez | Implementación eficiente de prueba de suavidad para polinomios[END_REF]. The elements are represented by polynomials of degree at most 508 whose coefficients are in F 3 6 . In this case n 1 = 6 and n 2 = 509. The initial splitting made with the Waterloo algorithm outputs two polynomials of degree 254. The probability that two independent and relatively prime polynomials of degree 254 over F 3 6 are simultaneously b-smooth is (1 -1/3 6 ) Pr 2 3 6 (254, b) [START_REF] Drmota | A rigorous proof of the Waterloo algorithm for the discrete logarithm problem[END_REF]. The term (1 -1/3 6 ) is negligible in practice for the values that we are considering. 4.3.1. Improvements. Our Algorithm 2 outputs one polynomial of degree 254, whose probability to be b-smooth is Pr 3 6 (n, b), i.e., the square root of the previous one. So we can take a much smaller b while reaching the same probability as before with the Waterloo algorithm. We list in Table 3a, p. 11, the values of b to obtain a probability between 2 -40 and 2 -20 . For instance, if we allow 2 30 trials, then we can set b = 28 with our algorithm, instead of b = 43 previously: we have Pr 2 3 6 (254, 43) = 2 -30.1 , and we only need to take b = 28 to get the same probability with this work: Pr 3 6 (254, 28) = 2 -29.6 . This will provide a good practical speed-up of the descent phase: much fewer elements need to be "reduced": this reduces the initial width of the tree, and they are of much smaller degree: this reduces the depth of the descent tree. 4.3.2. A 30-smooth initial splitting. The finite field is represented with n 1 = 6 and n 2 = 509, that is, as a first extension

F 3 6 = F p n 1 = F 3 [y]/(y 6 + 2y 4 + y 2 + 2y + 2), then a second extension F 3 6•509 = F 3 6 [x]/(I(x))
, where I(x) is the degree 509 irreducible factor of h 1 x q1h 0 , where q 1 = p n1 , h 1 = x 2 + y 424 x, and h 0 = y 316 x + y 135 . The generator is g = x + y 2 . As a proof of concept, we computed a 30-smooth initial splitting of the target T 0 = 508 i=0 (y ⌊π(3 6 ) i+1 ⌋ mod 3 6 )x i , with the parameters d = 3 × 509, d ′ = d/ gcd(d, n 1 ) = 509. Each trial g t T 0 produces d ′ = 509 polynomials to test for smoothness. We found that g 47233 T 0 = uvx 230 P , where u = 1 ∈ F 3 6 , v ∈ F 3 3•509 , and P is of degree 255 and 30-smooth. The equality (g 47233 T 0 )

p n -1 ℓ = (uvx 230 P ) p n -1 ℓ
is satisfied. The explicit value of P is available at https://members.loria.fr/AGuillevic/files/F3_6_509_30smooth.mag.txt.

The whole computation took less than 6 days (real time) on 48 cores Intel Xeon E5-2609 at 2.40GHz (274 core days, i.e., 0.75 core-years). This is obviously an overshot compared to the estimate of 2 26.6 , but this was done with a nonoptimized Magma implementation.

As a comparison, with the classical Waterloo algorithm, Adj et al. computed a 40-smooth initial splitting in 51.71 CPU (at 2.87GHz) years [2, Table 5.2, p. 87] and [START_REF] Adj | 6 * 509[END_REF]. They obtained irreducible polynomials of degree 40, 40, 39, 38, 37, and seven polynomials of degree between 22 and 35. They needed another 9.99 CPU years (at 2.66 GHz) to compute a classical descent from 40-smooth to 21-smooth polynomials. A complete comparison can be found in [START_REF] Adj | Computing discrete logarithms in cryptographically-interesting characteristic-three finite fields[END_REF] and [1]. In [START_REF] Adj | Computing discrete logarithms in cryptographically-interesting characteristic-three finite fields[END_REF], Adj et al. estimated that with our Algorithm 2 enriched as in Remarks 4.1 and 4.3, it is possible to compute discrete logarithms in F 3 6•709 at the same cost as in F 3 6•509 with the former Waterloo algorithm.

4.4.

Computing discrete logarithms in F 2 512 and F 2 1024 . In [28, §3.6] discrete logarithms in F 2 512 and F 2 1024 need to be computed modulo the full multiplicative group order 2 n -1. As pointed to us by R. Granger, our technique can be used to compute discrete logarithms in F 2 1024 . Our algorithm provides a decomposition of the target as the product uR where u is an element in the largest proper subfield F 2 512 , and P is an element of F 2 1024 of degree 512 instead of 1023. The discrete Table 3. Smoothness probabilities of polynomials over finite fields, comparison of the Waterloo algorithm and Algorithm 2. The values were computed with Odlyzko's induction formula [START_REF] Odlyzko | Discrete logarithms in finite fields and their cryptographic significance[END_REF] and Drmota and Panario's Theorem 4.5, as in [START_REF] Andrés Canales-Martínez | Implementación eficiente de prueba de suavidad para polinomios[END_REF].

(a) Probabilities for GF(3 6•509 ) two polynomials u(x), v(x) ∈ F 3 5 [x] of degree ⌊478/2⌋ = 239. Our Algorithm 1 outputs one polynomial of degree ⌊ 4 5 479⌋ = 383. This example is interesting because the smoothness probabilities are very close. We computed the exact values with Drmota-Panario's formulas, and give them in -49.34 ; and the cross-over point is for b = 24: in this case, we have Pr 2 3 5 (239, 24) = 2 -67.96 and Pr 3 5 (383, 24) = 2 -67. 59 , which is slightly larger. The probabilities would advise using the classical initial splitting with the Waterloo (extended GCD) algorithm. We remark that this algorithm would output two B 1 -smooth polynomials of degree (n 2 -1)/2. Each would factor into at least (n 2 -1)/(2B 1 ) irreducible polynomials of degree at most B 1 . Each such factor is sent as an input to the second step (descent step), that is, roughly n 2 B 1 factors. If we use Algorithm 2, the initial splitting will outputs one polynomial of degree 4/5n 2 = 383 that factors into at least 4/5n 2 /B 1 polynomials of degree at most B 1 , each of them sent as input to the second step, that is, the descent step is called 20% time less, and that would reduce the total width of the descent tree in the same proportion. Since the descent is the most costly part, and in particular, the memory size required is huge, this remark would need to be taken into consideration for a practical implementation.

As a proof of concept of our algorithm, we implemented in Magma our algorithm, took the same parameters, generator, and target as in [START_REF] Joux | Discrete logarithm record in characteristic 3[END_REF], and found a 50-smooth decomposition for the target given by the 471-th row of the matrix computed for g 23940 T 0 in 1239 core-hours (22.12 hours over 56 cores) on an Intel Xeon E5-2609 at 2.40GHz (compared to 5000 core-hours announced in [START_REF] Joux | Discrete logarithm record in characteristic 3[END_REF]).

The value can be found at https://members.loria.fr/AGuillevic/files/F3_5_479_50smooth.mag.txt. In our technique, we compute g t T 0 = uvR where u ∈ F 3 5 (this is the leading term of the polynomial), v ∈ F 3 479 , and R is 50-smooth. The discrete logarithm of u can be tabulated, however it remains quite hard to compute the discrete logarithm of v. Our technique is useful if it is easy (or not required) to compute discrete logarithms in the subfields.

Preliminaries before medium and large characteristic cases

In the first part of this paper, we were considering polynomials, and we wanted polynomials of smallest possible degree. Now we turn to the medium and large characteristic cases, where we do not have polynomials but ideals in number fields, and we want ideals of small norm. It requires testing whether large integers (norms) are smooth as fast as possible. We recall the results of Pomerance and Barbulescu on the early abort strategy. [START_REF] Pomerance | Analysis and comparison of some integer factoring algorithms, Computational methods in number theory[END_REF] introduced the Early Abort Strategy (EAS) to speed up the factorization of large integers, within Dixon's algorithm, the Morrison-Brillhart (continued fraction) algorithm, and the Schroeppel (linear sieve) and quadratic sieve, with several variations in the factorization sub-routine (trial-division, Pollard-Strassen method). The Early Abort Strategy provides an asymptotic improvement in the expected running time. Two versions are studied in [START_REF] Pomerance | Analysis and comparison of some integer factoring algorithms, Computational methods in number theory[END_REF]: one early-abort test, then many tests. In the relation collection step of the NFS algorithm, the partial factorization of the pseudonorms is done with ECM in time L Q [1/6] (Q = p n ), which is negligible compared to the total cost in L Q [1/3]. So Pomerance's EAS does not provide an asymptotic speedup, but a practical one. However, in the individual discrete logarithm computation, the initial splitting requires to find smooth integers (pseudonorms) of larger size: L Q [1]. This time the ECM cost is not negligible, and Pomerance's EAS matters. The speed-up was analyzed by Barbulescu in [START_REF] Barbulescu | Algorithmes de logarithmes discrets dans les corps finis[END_REF].

Pomerance's Early Abort Strategy. Pomerance in

Remark 5.1. Instead of the ECM test, it could be possible to use the hyperelliptic curve method test of H. Lenstra, Pila and Pomerance [START_REF] Lenstra | A hyperelliptic smoothness test, I[END_REF][START_REF]A hyperelliptic smoothness test. II[END_REF]. This was investigated for instance by Cosset [20,Chapter 4].

Pomerance's analysis is presented in the general framework of testing integers for smoothness. This is named smoothing problem in [10, Chapter 4]. In the individual discrete logarithm context, the numbers we want to test for smoothness are not integers in an interval, but pseudonorms, and their chances of being smooth do not exactly match the chances of random integers of the same size. However, we will make the usual heuristic assumption that for our asymptotic computations, the pseudonorms considered behave as integers of the same size. We give Pomerance's Early Abort Strategy with one test in Algorithm 3 and with k tests in Algorithm 4.

Algorithm 3: Pomerance's Early Abort Strategy (EAS)

Input: Integer m, smoothness bound B 1 , real numbers θ, b ∈]0, 1[ Output: B 1 -smooth decomposition of m, or ⊥ (m 0 , m 1 ) ← ECM (m, B θ 1 ) // cost: L B θ 1 [1/2, √ 2] // m 0 is a B θ 1 -smooth part of m // m 1 is the non-factorized part of m if m 1 ≤ m 1-b then (m 2 , m 3 ) ← ECM(m 1 , B 1 ) // cost: L B1 [1/2, √ 2] if m 3 = 1 then return B 1 -smooth decomposition m 1 , m 2 of m return ⊥
Writing the complexities as in Pomerance's paper, in terms of k early-abort tests, one obtains Theorems 5.2 and 5.3.

θ i < θ i+1 Output: B 1 -smooth decomposition of m, or ⊥ m i ← m i ← 0 S ← ∅ repeat (s i , m i+1 ) ← ECM (m i , B θi 1 ) // cost: L B θ i 1 [1/2, √ 2] // s i is a B θi 1 -smooth part of m i , m i+1 is not factorized S ← S ∪ s i m i ← m i+1 i ← i + 1 until (i > k) OR (m i = 1) OR (m i > m 1-i-1 j=0 bj ) if m i == 1 then return B 1 -smooth decomposition S of m return ⊥ the bound b i for 0 ≤ i ≤ k -1 on the remaining part m i in Algorithm 4 is b i = (2/3) 3(k-i) 19/(15 + 4(2/3) 3k ) ,
and the exponent θ i for 0 ≤ i ≤ k is

θ i = (4/9) k-i .
In Section 6.3, we will consider that pseudonorms behave in terms of smoothness like integers bounded by N e (instead of N ). We will need the following lemmas. We will mix Pomerance' strategy with our new initial splitting step to improve its running time.

LLL algorithm.

We recall an important property of the LLL algorithm [START_REF] Lenstra | Factoring polynomials with rational coefficients[END_REF] that we will widely use in this paper. Given a lattice L of Z n defined by a basis given by an n × n matrix L, and parameters 1 4 < δ < 1, 1 2 < η < √ δ, the LLL algorithm outputs a (η, δ)-reduced basis of the lattice. the coefficients of the first (shortest) vector are bounded by

(δ -η 2 ) n-1 4 det(L) 1/n .
In the remainder of this paper, we will simply denote by C this LLL approximation factor.

5.3. NFS and Tower variants.

5.3.1.

Settings. There exist many polynomial selection methods to initialize the NFS algorithm for large and medium characteristic finite fields. We give in Table 4 the properties of the polynomials that we need (degree and coefficient size) to deduce an upper bound of the pseudonorm, as in (5.3), and (5.4). Three polynomials define the NFS setting: ψ, f 0 , f 1 , where f 0 , f 1 are two polynomials of integer coefficients, irreducible over Q, of degree ≥ n, defining two nonisomorphic number fields, and whose GCD modulo p is an irreducible polynomial ψ of degree n, used to define the extension

Q K f0 K f1 deg f 0 ≥ n deg f 1 ≥ n (a) NFS number fields Q K h deg h = n 1 K f0 K f1 deg f 0 ≥ n 2 deg f 1 ≥ n 2 (b) Tower-NFS number fields
Z[x] Z[x]/(f 0 (x)) Z[x]/(f 1 (x)) F p n = F p [x]/(ϕ(x)) ρ 0 ρ 1 (a) NFS diagram for F p n R h = Z[y]/(h(y)) R h [x] R h [x]/(f 0 (x)) R h [x]/(f 1 (x)) F (p n 1 ) n 2 = (R h /pR h )[x]/(ψ(x)) ρ 0 ρ 1 (b) Tower-NFS diagram
F p n = F p [x]/(ψ(x)).
In a tower-NFS setting, one has n = n 1 n 2 , n 1 , n 2 = 1 and four polynomials are defined: h, ψ, f 0 , f 1 , where deg h = n 1 and h is irreducible modulo p, deg ψ = n 2 and ψ is irreducible modulo p, and gcd(f 0 mod (p, h), f 1 mod (p, h)) = ψ. It can be seen as a generalization of the NFS setting as follows: writing n = n 1 n 2 , one starts by defining a field extension F p n 1 = F p [y]/(h(y)) and then adapting any previously available polynomial selection designed for NFS in GF(p n2 ), using F p n 1 as the base field instead of F p . When gcd(n 1 , n 2 ) > 1, the polynomials f 0 , f 1 , resp., ψ, will have coefficients in Q[y]/(h(y)), resp., F p n 1 , instead of Q, resp., F p . Then one defines the second extension F p n 1 [x]/(ψ(x)) of degree n 2 = deg x ψ. Again, to cover all the cases, we consider F p n = F (p n 1 ) n 2 . The NFS case will correspond to n 1 = 1, n 2 = n and the original TNFS case to n 1 = n, n 2 = 1.

5.3.2.

Pseudonorm and upper bound. Let f be a monic irreducible polynomial over Q, and let K = Q[x]/(f (x)) be a number field. Write T ∈ K as a polynomial in x:

T = deg f -1 i=0 a i x i .
The norm is defined by a resultant computation:

(5.1)

Norm K/Q (T ) = Res(f, T ) .
In the NFS case, we will consider elements expressed as polynomials in x whose coefficients are integers. We define the pseudonorm as the resultant of the element with the given polynomial f :

T = deg f -1 i=0 a i x i , pseudonorm(T (x)) = Res(T (x), f (x)) .
We use Kalkbrener's bound [START_REF] Kalkbrener | An upper bound on the number of monomials in determinants of sparse matrices with symbolic entries[END_REF]Corollary 2] for an upper bound:

(5.2) | Res(f, T )| ≤ κ(deg f, deg T ) f deg T ∞ T deg f ∞ , where κ(n, m) = n+m n n+m-1 n
and f ∞ = max 0≤j≤deg f |f j | is the absolute value of the largest coefficient. An upper bound for κ(n, m) is (n + m)!. We will use the following bound in Section 6:

(5.3) Norm K f /Q (T ) ≤ (deg f + deg T )! f deg T ∞ T deg f ∞ .
In a Tower-NFS case, we nest two resultants:

T = deg f -1 i=0 deg h-1 j=0 a ij y j x i , pseudonorm(T (x, y)) = Res y (Res x (T (x), f (x)), h(y)) .
A bound is [45, §A Lemma 2]

(5.4)

|N K f /Q deg x P i=0 deg h-1 j=0 a ij α j h α i f | < a ij deg h deg f ∞ f deg x P deg h ∞ h (deg x P +deg f )(deg h-1) ∞ D(deg h, deg f ) , where a ij ∞ = max i,j |a ij | and D(d 1 , d 2 ) is a combinatorial term, D(d 1 , d 2 ) = ((2d 2 -1)(d 1 -1) + 1) d1/2 (d 1 + 1) (2d2-1)(d1-1)/2 ((2d 2 -1)!d 2d2 1 ) d1 .

Faster Initial Splitting with NFS and Tower variants for medium and large characteristic finite fields

We apply Algorithm 1 to the medium and large characteristic cases. For a general exposition, we assume that we are in a tower setting, where Q = p n = (p n1 ) n2 . The elements of F p n are represented as T = n1-1 i=0 n2-1 j=0 a i,j y j x i . the NFS setting corresponds to n 1 = 1, n 2 = n. When n is prime, the tower setting is n 1 = n, n 2 = 1 but our algorithm does not apply. Denote by h(y) the polynomial defining the field F p n 1 and by ψ the polynomial defining the degree n 2 extension F (p n 1 ) n 2 . Here we are not interested (only) in computing a preimage of degree as small as possible, but more generally one whose size of pseudonorm is as small as possible. According to the bounds (5.3), (5.4), we need to combine small coefficients a i,j (to reduce the contribution of a ij

deg h deg f ∞
) with a small degree in x (to reduce the contribution of f

deg x P deg h ∞
), and balance the two terms to find a pseudonorm of smaller size.

6.1. The algorithm. We start again with the same idea as in Algorithm 1: let d be the largest proper divisor of n, with 1 < d < n. 5 Assume we want to obtain a preimage P ∈ F p n 1 [x] of the target, of degree (nd)/n 1 ≤ deg P < deg f . We will use relations of the form P = uT (mod ψ), where u p d -1 = 1 as in (3.1). We use the relations

x i y j p = 0 (mod p, h, ψ) for 0 ≤ ij < d , P = U i T (mod p, h, ψ) ,
where {1, U, . . . , U d-1 } is a polynomial basis of F p d and where

x i y j ψ = 0 (mod p, h, ψ) for 0 ≤ j < n 1 , 0 ≤ i < deg(P ) -n 2 .
We define the lattice of these relations and we obtain a matrix

L n1(deg P +1)×n1(deg P +1) =                    p . . . p coeff(T ) coeff(U T ) . . . coeff(U d-1 T ) coeff(ψ) . . . coeff(x i (y j ψ mod h(y)))                   
We want to obtain a matrix in row-echelon form. The d first rows and the n 1 (deg Pn 2 ) last rows are in row-echelon form by construction. We compute Gaussian elimination to obtain a reduced row-echelon form for the rows U i T . We use F plinear combinations of these rows, and we allow divisions in F p so that the leading coefficient is one. We then obtain a square matrix of dimension n 1 (deg P +1) in rowechelon form. Now at this point we apply a lattice reduction algorithm such as LLL or BKZ to reduce the size of the coefficients of L. We obtain a matrix R whose first row has coefficients bounded by C LLL det(L) 1/(n1(deg P +1)) = p (n-d)/(n1(deg P +1)) .

6.2. Properties and pseudonorm size bound.

Proposition 6.1. The preimage P output by Algorithm 5 satisfies log g ρ(P ) ≡ log g g t T 0 = log g T 0 + t mod Φ n (p), where ρ : Z[x, y] → F (p n 1 ) n 2 was defined in Figure 2.

5 d = deg(h) = n 1
is the case studied independently in the preprint [START_REF] Zhu | Improvements on the individual logarithm step in extended tower number field sieve[END_REF]. Since an earlier version of this work was presented at Asiacrypt 2015 and ECC 2015, and the question of how to use larger subfields raised in discussions at these conferences, it is not surprising that other researchers though of similar ideas to improve individual discrete logarithms in the same time period.

Algorithm 5: Initial splitting, Tower-NFS setting Input: Finite field F p n , n = n 1 n 2 , monic irreducible polynomials h, ψ s.t.

F p n 1 = F p [y]/(h(y)), F (p n 1 ) n 2 = F p n 1 [x]/(ψ(x)), prime order subgroup ℓ | Φ n (p)
, generator g (of the order ℓ subgroup), target T 0 ∈ F p n , degree of the preimage deg P , polynomial f i , smoothness bound B 1 Output: t ∈ {1, . . . , ℓ -1}, P ∈ Z[x] s.t. log g ρ(P ) ≡ t + log g T 0 , and the pseudonorm Res y (Res

x (P, f i ), h) is B 1 -smooth d ← the largest divisor of n, 1 ≤ d < n Compute a polynomial basis (1, U, U 2 , . . . , U d-1 ) of the subfield F p d , where U satisfies U p d -1 = 1 ∈ F p n repeat take t ∈ {1, . . . , ℓ -1} uniformly at random T ← g t T 0 ∈ F p n L ←                    p . . . p coeff(T ) coeff(U T ) . . . coeff(U d-1 T ) coeff(ψ) . . . coeff(x i (y j ψ mod h(y)))                    Compute a F p -
U i T , 0 ≤ i ≤ d-1, i.e., an element d-1 i=0 λ i U i T , where λ i ∈ F p . We can factor T in the expression. Each element u j = d-1 i=0 λ i U i satisfies u p d -1 j = 1, i.e., is in F p d by construction. So each row represents an el- ement T j = u j T , where u p d -1 j = 1 (u j ∈ F p d ), so that log T j ≡ log T mod Φ n (p) by Lemma 2.4.
The second part of the proof uses the same argument: the short vector output by the LLL algorithm is a linear combination of the rows of the matrix N . Each row represents either 0 or a F p d -multiple T j of T , hence the short vector is also a F p d -multiple of T . We conclude thanks to Lemma 2.4, that log ρ(P ) ≡ log T mod Φ n (p). Proposition 6.2. The pseudonorm of P in Algorithm 5 has size

(6.1) | Res y (Res x (P, f i ), h)| = O Q (1-d n ) deg f i deg x P +1 f i n1 deg x P ∞ assuming that h ∞ = O(1).
Proof of Proposition 6.2. The matrix N computed in Algorithm 5 is a square matrix of (deg x P + 1)n 1 rows and columns, whose coefficients are in

F p . Its determi- nant is det N = p n-d = Q 1-d/n
. Using the LLL algorithm for the lattice reduction, the coefficients of the shortest vector P are bounded by

CQ (1-d/n)/((deg x P +1)n1) ,
where C is the LLL factor. We obtain the bound (6.1) according to the bound formula (5.4), and neglecting the combinatorial factor D(n 1 , deg f i ). Moreover in the Tower-NFS setting, the polynomial selection is designed such that h ∞ = O(1).

We finally obtain the following.

Theorem 6.3. Let GF(p n ) be a finite field, and let d be the largest divisor of n, d < n, and d = 1 if n is prime. Let n = n 1 n 2 and h, ψ, f i be given by a polynomial selection method. Let T ∈ F (p n 1 ) n 2 be an element which is not in a proper subfield of F p n . Then there exists a preimage

P ∈ Z[x, y] of T , of any degree (in x) between ⌊n 2 -d/n 1 ⌋ and deg f i -1, of coefficients bounded by O(Q (1-d n ) 1 (deg P +1)n 1
), and such that when P is mapped in F (p n 1 ) n 2 as ρ(P ), its discrete logarithm is equal to the discrete logarithm of T modulo Φ n (p) (and in particular modulo any prime divisor ℓ of Φ n (p)), that is, log ρ(P ) ≡ log T mod Φ n (p) . The degree of P in x and the polynomial f i can be chosen to minimize the resultant (pseudonorm):

min i min ⌊n2-d/n1⌋≤deg x P ≤deg fi-1 f i n1 deg x P ∞ Q (1-d n ) deg f i deg x P +1 .
We recall in Table 4 the degree and coefficient sizes of the polynomial selections published as of July 2017. Corollary 6.4. With the notation of Table 4 and the NFS setting corresponding to n 2 = n and n 1 = 1,

(1) For the polynomial selection methods where there is a side i such that f i ∞ = O(1) (GJL, Conjugation, Joux-Pierrot and Sarkar-Singh up to now), we do the initial splitting on this side and choose

deg x P = deg x f i -1 to obtain the smallest norm: | Res y (Res x (P, f i ), h)| = O Q 1-d n .
We obtain the same bound for NFS and its tower variants.

(2) When f i ∞ = Q 1/(2n) as for the JLSV 1 method, the bound is

Q (1-d n ) n 2 deg x P +1 + deg x P 2n 2 . When deg x P = deg x f i -1 = n 2 -1, one obtains Q 3 2 -d n -1 2n 2 .
In the NFS setting, n 2 = n, while in the tower setting, n 2 < n and the pseudonorm is slightly smaller.

(3) When f i ∞ = Q 1/(n1(D+1)) as for the JLSV 2 method, the lower bound is

Q deg x P D+1 +(1-d n ) n 2 deg x P +1
on the f 0 -side where deg f 0 = n 2 , and it is

Q deg x P D+1 +(1-d n ) D deg x P +1 on the f 1 -side, where deg f 1 = D ≥ n 2 .
According to the value of n, one can decide which value of deg x P will produce a smaller norm. Table 4. Properties: degree and coefficient size of the main polynomial selection methods for NFS-DL in F Q , where Q = p n .

We give a bound on the coefficient size of the polynomials with the notation f i ∞ = O(x). To lighten the notation, we only write x, without O(). In the Joux-Pierrot method, the prime p can be written p = p x (x 0 ), where p x is a polynomial of tiny coefficients and degree at least 2. This table takes into account the methods published until July 2017. method deg h [START_REF] Matyukhin | Effective version of the number field sieve for discrete logarithms in the field GF(p k )[END_REF][START_REF] Barbulescu | Algorithmes de logarithmes discrets dans les corps finis[END_REF][START_REF] Barbulescu | Improving NFS for the discrete logarithm problem in non-prime finite fields[END_REF] D + 1 D ≥ n log p Q 1/(D+1) Conjugation [START_REF] Barbulescu | Improving NFS for the discrete logarithm problem in non-prime finite fields[END_REF] 2n

deg f 0 deg f 1 f 0 ∞ f 1 ∞ NFS JLSV 1 [39] n n Q 1/2n Q 1/2n JLSV 2 [39] n D > n Q 1/(D+1) Q 1/(D+1) GJL
n log p Q 1/2n Joux-Pierrot [41] p = p x (x 0 ) n(deg p x ) n log p Q 1/(n deg px) Sarkar-Singh [59] n = n 1 n 2 , D ≥ n 2 (D + 1)n 1 Dn 1 log p Q 1/(n1(D+1)) Tower-NFS TNFS + base-m [14] n D 1 p 1/D p 1/D Tower-JLSV 1 n = n 1 n 2 n 1 n 2 n 2 Q 1/(2n) Q 1/(2n)
Tower-JLSV 2 n = n 1 n 2 [START_REF] Kim | Extended tower number field sieve: A new complexity for medium prime case[END_REF][START_REF] Kim | Extended tower number field sieve: A new complexity for the medium prime case[END_REF] n

1 n 2 D ≥ n 2 Q 1/(n1(D+1)) Q 1/(n1(D+1))
Tower-GJL n = n 1 n 2 [START_REF] Kim | Extended tower number field sieve: A new complexity for the medium prime case[END_REF]

n 1 D + 1 D ≥ n 2 log p Q 1/(n1(D+1))
Tower-Conjugation n = n 1 n 2 [START_REF] Barbulescu | An appendix for a recent paper of Kim[END_REF][START_REF] Kim | Extended tower number field sieve: A new complexity for the medium prime case[END_REF][START_REF] Kim | Extended tower number field sieve with application to finite fields of arbitrary composite extension degree[END_REF]]

n 1 2n 2 n 2 log p Q 1/(2n)
Tower-Joux-Pierrot n = n 1 n 2 , p = p x (x 0 ) [START_REF] Kim | Extended tower number field sieve: A new complexity for the medium prime case[END_REF][START_REF] Kim | Extended tower number field sieve with application to finite fields of arbitrary composite extension degree[END_REF]

n 1 n 2 (deg p x ) n 2 log p Q 1/(n deg px)
Tower-Sarkar-Singh n = n 1 n 2 n 3 , D ≥ n 3 [START_REF] Sarkar | A general polynomial selection method and new asymptotic complexities for the tower number field sieve algorithm[END_REF][START_REF]Tower number field sieve variant of a recent polynomial selection method[END_REF][START_REF]A generalisation of the conjugation method for polynomial selection for the extended tower number field sieve algorithm[END_REF]]

n 1 (D + 1)n 2 Dn 2 log p Q 1/(n1n2(D+1))
6.3. running time. To apply Lemma 5.4 to the initial splitting case, we make the usual heuristic assumption that the pseudonorms of the elements g t T 0 behave asymptotically like random integers of the same size. Their size is O(Q e ), so we replace N e by Q e . The basis {1, U, . . . , U d-1 } can be precomputed. The cost of computing the U i T for 0 ≤ i ≤ d -1 is at most dn 2 multiplications in F p with a schoolbook multiplication algorithm. We can roughly upper-bound it by O(n 3 ). The time needed to compute the reduced row-echelon form of a d × n matrix is in O(n 3 ) which is polynomial in n [START_REF] Dumas | Handbook of finite fields[END_REF]. These two complexities are asymptotically negligible compared to any L Q [α > 0]. We obtain the following.

Corollary 6.5. The running time of the initial splitting step with Algorithm 5 to find a B-smooth pseudonorm, where the pseudonorm has size O(Q e ) for a fixed real number e > 0 determined by the polynomial selection (Table 4, two right-most columns), is

(1) L Q [1/3, c = (3e) 1/3 ] with ECM to perform the smoothness tests;

(2) L Q [1/3, c = (3e) 1/3 (23/27) 2/3 ] with ECM and EAS;

(3) L Q [1/3, c = (3e) 1/3 ((15 + 4(2/3) 3k )/19) 2/3 ] with ECM and k-EAS. For each case, the lower bound was obtained for B = L Q [2/3, e/c]. Corollary 6.4 gives a bound on the size of the pseudonorms, from which we can deduce e to apply Corollary 6.5, and get the expected running time.

Examples

Example 7.1. Let p = ⌊10 25 π⌋ + 7926 = 31415926535897932384634359 be a 85-bit prime made of the first 26 decimals of π so that F p 6 is a 509-bit finite field. Moreover, Φ 6 (p) = p 2p + 1 is a 170-bit prime, we denote it by ℓ = 986960440108935861883947021513080740536833738706523. We want to compute discrete logarithms in the order-ℓ cyclotomic subgroup of F p 6 . The JLSV 1 method computes two polynomials f 0 , f 1 , where deg f 0 = deg f 1 = 6, and f i ∞ ≈ p 1/2 . In our example, we have log 2 f 0 ∞ = 44.67 and log 2 f 1 ∞ = 46.67 (and log 2 p/2 = 42.35):

f 0 = x 6 -11209975711932 x 5 -28024939279845 x 4 -20 x 3 +28024939279830 x 2 + 11209975711938 x + 1 f 1 = 5604994576830 x 6 + 20986447533158 x 5 -31608799819555 x 4 -112099891536600 x 3 -52466118832895 x 2 + 12643519927822 x +5604994576830.
Since f 0 is already of degree 6 and monic, it can define the extension F p 6 = F p [x]/(f 0 (x)). Let T 0 be our target in F p 6 whose coefficients are made of the decimals of π (starting at the 26-th decimal, since the first 25 ones were already used for p):

T 0 = 6427704988581508162162455 x 5 + 16240052432693899613177738 x 4 +4509390283780949909020139 x 3 + 3868374359445757647591444 x 2 +8209755913602112920808122 x + 3279502884197169399375105.

Let g = x + 3 be a generator of F p 6 . Let (1, U, U 2 ) be a polynomial basis of F p 3 considered as an implicit subfield of F p 6 , where U = g 1+p 3 = Norm F p 6 /F p 3 (g). We run Algorithm 5 and find that the fourth preimage of T = g 812630 T 0 gives a 61-smooth pseudonorm. We compute the reduced row-echelon form A common choice for the factor basis would be to set its smoothness bound to 30 or 32 bits. There are six prime ideals whose norm is larger than 30 bits, and that should be retreated to reach the factor basis. This initial splitting, testing all pseudonorms obtained for g i T 0 , i from 0 to 930000, that is, 5.58

=      T U T . . . U n/2-1 T      and compute M =       m 1,1 . . . m 1, n 2 -1 1
1 to be B 1 = L Q [2/3, ((1 -1 2n ) 2 /3) 1/3
]. The running time of Algorithm 5 will be

L Q [1/3, (3(1 -1 2n )) 1/3 ].
8. Optimal representation: monic polynomial of degree ϕ(n)

In Section 3, we exploited the largest proper subfield F p d of F p n to find an alternative representation of a given element T ∈ F p n , with n-d nonzero coefficients, and d -1 coefficients (in F p ) set to zero. The key ingredient was to compute an expression of the form P = uT , where P has d -1 coefficients set to zero, and u ∈ F p d , so that we have the equality (P/T ) (p n -1)/Φn(p) = 1. We can generalize this strategy: given an element T in the cyclotomic subgroup of F p n , of order Φ n (p), we would like to compute an element P ∈ F p n such that (P/T ) (p n -1)/Φn(p) = 1 and P has only ϕ n (p) = deg Φ n (x) non-zero coefficients in F p . To achieve that, we would like to compute an expression T = u 1 u 2 . . . u i P, where each u i is in a proper subfield F p d i of F p n .

Given an element T ∈ F p n such that T (p n -1)/Φn(p) = 1 (in other words, its order in the cyclotomic subgroup of F p n is not zero), we can sometimes compute an element P with ϕ(n) non-zero coefficients, where ϕ(n) is the Euler totient function, plus a monic leading term. Since in Algorithm 5 we do not need a one-to-one correspondence between the given elements of the cyclotomic subgroup on one hand, and their representation with only ϕ(n) non-zero non-one coefficients on the other hand, we can just solve a system of equations even if we do not expect a solution at all times. If no such compact representation is found, one picks a new t and tests for the next g t T 0 . To define the system we need to solve, we list all the distinct subfields F p d of F p n that are not themselves contained in another proper subfield, compute a polynomial basis for each of them, and allow a degree of freedom for the coefficients to be ϕ(d) for each subfield F p d . If we consider the system as a Gröbner basis computation, it becomes very costly even for F p 30 , where we need to handle nϕ(n) -1 = 21 variables. We give a numerical example for F p 6 .

What we do is different than what is done in XTR and CEILIDH compact representations. In the XTR cryptosystem [START_REF] Arjen | The XTR public key system[END_REF], the elements of the cyclotomic subgroup of F p 6 are represented with an optimal normal basis over F p 2 , also in normal basis representation. Only their trace over F p 2 is considered for representation, storage, and transmission. In [START_REF] Van Dijk | Asymptotically optimal communication for torusbased cryptography[END_REF][START_REF] Van Dijk | Practical cryptography in high dimensional tori[END_REF], the aim is to define a one-to-one correspondence between the elements in the torus of F p n and the set of coefficients (F p ) ϕ(n) . This optimal compression was achieved for n = 6 but not for n = 30. These techniques are not compatible with the representation of the elements in the NFS algorithm: one chooses a representation by choosing two polynomials f 0 , f 1 that define the two number fields involved in the algorithm. One cannot change the representation afterwards: the elements in the individual discrete logarithm phase should be represented in the same way as the elements of the factor basis. 8.1. Compressed representation of elements in the cyclotomic subgroup of F p 6 by a monic polynomial of degree 2. We consider the finite field F p 6 . We will use the two subfields F p 2 and F p 3 to cancel three coefficients. Let U ∈ F p 6 such that (1, U, U 2 ) is a basis of F p 3 ⊂ F p 6 . Let V ∈ F p 6 such that (1, V ) is a basis of F p 2 ⊂ F p 6 . We want to solve uvwT = (u 0 + u 1 U + u 2 U 2 )(v 0 + v 1 V )wT = P ,

where u = u 0 + u 1 U + u 2 U 2 ∈ F p 3 , v = v 0 + v 1 V ∈ F p 2 ,
w ∈ F p , and P ∈ F p 6 is represented by a monic polynomial in x of degree 2. To simplify, we set u 2 = v 1 = 1 so that we obtain equations where we can recursively eliminate the variables by computing resultants. We compute u, v, w such that uvwT = P , where P = a 0 + a 1 x + x 2 is monic of degree 2. We define the lattice

L =   p 0 0 0 p 0 a 0 a 1 1   .
The determinant of L is p 2 hence LLL(L) computes a short vector P of coefficient size bounded by Cp 2/3 , where C is the LLL approximation factor (we can take C ≈ 1 in this practical case). The pseudonorm of P will be in the JLSV 1 case | Res(P, f )| ≈ P 6 ∞ f 2 ∞ = p 5 = Q 5/6 . This is better than the bound Q 11/12 obtained with the cubic subfield cofactor method. This specific method can be generalized to specific cases of finite fields where reducing as much as possible the degree of the target is the best strategy, as in Example 8.1. This technique was implemented in [START_REF] Grémy | Computing discrete logarithms in F p 6[END_REF] for computing a new discrete logarithm record in F p 6 of 422 bits.

Example 8.1. We take the same finite field parameters as in Example 7.1, where F p 6 = F p [x]/(f (x)). g = x + 3 is a generator of F p 6 . (1, U, U 2 ) where U = g 1+p 3 is a basis of F p 3 and (1, V ) where V = g 1+p 2 +p 4 is a basis of F p 2 . We solve the system (u 0 + u 1 U + U 2 )(v 0 + V )T = P where u i , v i ∈ F p and P is monic and represented by a polynomial of degree 2 instead of 5. We ran Algorithm 5 with this modification on the same machine (Intel Xeon E5-2650 at 2.0GHz with hyperthreading turned on), from g 0 T 0 to g 90000 T 0 . On average, the set of I candidates g i T 0 led to six times more monic degree two polynomials P i . We found that the third polynomial output for T = g 60928 T 0 has a 64-bit-smooth pseudonorm. Testing the 90000 g i T 0 (that is, 2.7 • 10 5 pseudonorms) took 1.2 core-day: of 429 bits, which corresponds to the estimate log 2 Q 5/6 = 423 bits. We still have log g ρ(P ) ≡ log g T 0 + 60928 mod ℓ. The pseudonorm is 64-bit-smooth, and its factorization into prime ideals is The first three ideals are small enough to be in the factor basis, and eight ideals on side 0 remain to be descended.

Conclusion

The algorithms presented in this paper were implemented in Magma and used for cryptographic-size record computations. It was shown in [START_REF] Adj | Computing discrete logarithms in cryptographically-interesting characteristic-three finite fields[END_REF] that combined with a practical variant of Joux's algorithm, our Algorithm 2 allows to compute a discrete logarithm in the finite field F 3 6•709 at the same cost as in F 3 6•509 with the previous state of the art. The large characteristic variant (Algorithm 5) was used in [START_REF] Grémy | Computing discrete logarithms in F p 6[END_REF] for a 422-bit record computation in F p 6 . It would be interesting to be able to generalize it further, to be able to exploit at the same time several subfields, and provide a practical implementation of it for cryptographic sizes.

4. 2 .

 2 Complexity analysis. 4.2.1. Cost of computing one preimage P ∈ F p n 1 [x] in the initial splitting step. We use the notation of Algorithm 2: let d be the largest proper divisor of n

  Odlyzko gave the following estimation for Pr q (b; d) in [54, (4.5), p. 14]. (4.1) Pr q (b, d) -1 = exp (1 + o(1)) d b log e d b for d 1/100 ≤ b ≤ d 99/100 . Writing the smoothness bound degree b = log L Q [α b , c b ]/ log p n1 to match Odlyzko's convention b = c b n α b

Corollary 4 . 6 (

 46 [START_REF] Drmota | A rigorous proof of the Waterloo algorithm for the discrete logarithm problem[END_REF] Theorem 1]). Let δ > 0 be given. Then we have, uniformly for

Theorem 5 . 2 (

 52 [10, § 4.3]). The expected running time of the smoothing problem of an integer N with Pomerance's EAS and the ECM smoothness test is L N [1/3, c] where c = (23/3) 2/3 /3, the smoothness bound is B = L N [2/3, γ], where γ = 1/c, θ = 4/9, and b = 8/23.

Theorem 5 . 3 (

 53 [10, § 4.5 Th. 4.5.1]). The expected running time of the smoothing problem of an integer N with k tests of Pomerance's EAS and the ECM smoothness test is L N [1/3, c] where c = 3 1/3 ((15 + 4(2/3) 3k )/19) 2/3 , the smoothness bound is B = L N [2/3, γ], where γ = 1/c , Algorithm 4: Pomerance's Early Abort Strategy with k tests (k-EAS) Input: Integer m, smoothness bound B 1 , number of tests k ≥ 0, array of positive real numbers b = [b 0 , b 1 , . . . , b k ] where 0 < b i ≤ 1, and k i=0 b i = 1 array of positive real numbers θ = [θ 0 , . . . , θ k = 1] where

Lemma 5 . 4 (

 54 [18, §4.1],[34, Lemma 1] Running time of B-smooth decomposition of integers with ECM). Let N i be integers taken uniformly at random and bounded by N e , for a fixed real number e > 0. Write B = L N [α B , γ] the smoothness bound. Then the expected running time to obtain a B-smooth N i , using ECM for B-smooth tests, isL N [1/3, (3e) 1/3 ], obtained with B = L N [2/3, e/c = (e 2 /3) 1/3 ].Lemma 5.5 ([55, 10] Running time of B-smooth decomposition of integers with ECM and k-EAS). Let N i be integers taken uniformly at random and bounded by N e , for a fixed real number e > 0. Write B = L N [α B , γ] for the smoothness bound. Then the expected running time to obtain a B-smooth N i , using ECM for B-smooth tests and Pomerance's Early Abort Strategy with one test, is L N [1/3, c = (3e) 1/3 (23/27) 2/3 ], obtained with B = L N [2/3, e/c]. The expected running time with k-EAS is L N [1/3, c = (3e) 1/3 ((15 + 4(2/3) 3k )/19) 2/3 ] , with B = L N [2/3, e/c].

Figure 1 .

 1 Figure 1. Extensions of number fields for NFS and tower variants

Figure 2 .

 2 Figure 2. NFS and tower variant diagrams for F p n

  reduced row echelon form of the rows nd + 1 to n of L N ← LatticeReduction(L) P ← polynomial in Z[y, x] made of the shortest vector output by the LatticeReduction algorithm until Res y (Res x (P, f i ), h) is B 1 -smooth // ECM, ECM+EAS, or ECM+k-EAS return t, P , factorization of Res y (Res x (P, f i ), h) Proof of Proposition 6.1. Each row of the row-echelon matrix M represents a F plinear combination of the d elements

m 10 m 11 m 12 m 13 1 0 m 20 m 21 m 22 m 23 m 24 1 .

 1 30930778358987253373198053 m 01 = 16172276732961477886471865, m 02 = 251875570676859576731124 m 10 = 8981071706647180870633008, m 11 = 26297121233008662476505921 m 12 = 4999545867425989707589927, m 13 = 4380553940470247124926451 m 20 = 4787502941827866787698085, m 21 = 18855419729462744536987506 m 22 = 15450347628775338768673252, m 23 = 31092163492444411597011243 m 24 = 9824382756181109886988461.Then we reduce with the LLL algorithm the following lattice defined by the (6 × 6)matrix, where m ij stands for the coefficient at row i and column j of the above matrix M , and m i,3+i = 1: Each row of LLL(N ) gives us a preimage P ∈ Z[x] of short coefficients such that log 2 P ∞ ≈ log 2 p/2 = 42.34 bits and log ρ(P ) ≡ log T mod ℓ (in other words, (T /ρ(P ))

u

  = 12307232765040677532260293 + 18116887363761988927417497 U + U 2 v = 30422514788629575495025401 + V w = 21470888563719305004900851 P = uvwT = x 2 + 479190487430850236087613 x + 6943966382910680737931850 .We checked that (P/T ) p 6 -1 ℓ = 1, meaning that log g P = log g T = 60928 + log g T 0 . Then we reduce the lattice defined by the matrix polynomials of smaller coefficients, the third one beingR = 107301402613441938 x 2 -32014642452727111 x + 60125316588415598 whose pseudonorm is Res(R, f ) = 12474200655939339762647720853686893930822373172685245800138935320 22514918959041066623605301421497621878867497302294873400285994921

  , in F p n 1 [x], of degree n 2 -⌈d/n 1 ⌉ in x and coefficients in F p n 1 , and such that We use Algorithm 1 to compute P . The matrix has full rank since the U i s form a polynomial basis of F p d ′ . The linear combinations involve T and elements in F p d/ gcd(n 1 ,d) and F p gcd(n 1 ,d) that are in the proper subfield F p d by construction. The first row after Gaussian elimination will have at least d/n 1 -1 coefficients equal to zero at the right, and will represent a polynomial P of degree at most n 2d/n 1 , that satisfies P = uT (mod ψ) where u = λ i U i ∈ F p d , since in the process, T was multiplied only by elements whose images in F (p n 1 ) n 2 are in the subfield F p d . We have ρ(P ) = uT , u ∈ F p d , and the equality of logarithms follows by Lemma 2.4.

	P of T log(ρ(P )) = log T mod Φ n (p) .
	Proof.
	with only F p gcd(n 1 ,d) -linear combinations
	P (x) ← polynomial from the coefficients of the first row of L return P (x)
	p n Theorem 3.1. Let F p n be a finite field represented as a tower F (p n 1 ) n 2 . Let T ∈ F * be an element which is not in a proper subfield of F p n . Let d be the largest proper
	divisor of n, 1 < d < n (n is not prime). Assume that T is represented by a
	polynomial in F p n 1 [x] of degree larger than n 2 -d/n 1 . Then there exists a preimage

  Table 3b, p. 11. We obtain Pr 2 3 5 (239, 50) = 2 -20.96 (Waterloo) and Pr 3 5 (383, 50) = 2 -22.96 , i.e., our Algorithm 2 would be four times slower compared to Joux's and Pierrot's record; Pr 2 3 5 (239, 40) = 2 -30.61 and Pr 3 5 (383, 40) = 2 -32.34 ; Pr 2 3 5 (239, 30) = 2 -48.46 and Pr 3 5 (383, 30) = 2

  The fourth row has coefficients of at most 41.82 bits and givesP = 482165402365 x 5 + 3892831179802 x 4 + 2694050932529 x 3 +2325450478817 x 2 + 1117470283668 x + 3688595236671 .

	The pseudonorm of P w.r.t. f 0 is
	Res(P, f 0 ) =
	32601551184187978602887820222780280368556791213406352787959859478882009\
	89411710052105812763285379877699363515358275429392312189582741360186561
	of 471 bits, which is very close to log 309800481
	1402910243559283, x + 1034551157262971
	1587503571970639, x + 524543605465730
	36834399852305717, x + 24916507207930752
	242270403627311729, x + 170018299727614229
	1070632553963863603, x + 408232161861505290
	4305864084909925127, x + 3252872861595329896 .

2 Q 11/12 = 466 bits. Its factorization in prime ideals of K f0 is 3, x + 2

3

11, x + 5 17, x + 4 67, x + 44 2011, x + 463 501997, x + 18312 340575947, x + 27999767 506032577, x + 177467846 604579099, x +

  • 10 6 pseudonorms, with our Magma implementation, took 0.95 day on one node of 16 physical cores (32 virtual cores thanks to hyperthreading) Intel Xeon E5-2650 at 2.0GHz, that is, 15.2 core-days.Example 7.2 (A more general example with NFS). Assume that n is even and let T ∈ F p n . Compute a polynomial basis (1, U, U 2 , . . . , U n/2-1 ) of the subfield F p n/2 . Let

	L

  in reduced echelon form in the lower half. Moreover, if deg(f ) > n, then we add (deg fn -1) rows made of the coefficients of x i ψ where F p n = F p [x]/(ψ(x)), for 0 ≤ i < deg fn -1. Finally we apply the LLL algorithm to this matrix. The short vector gives us a preimage P whose pseudonorm is bounded by Q 1/2 , with a polynomial selection such that f ∞ = O(1) (such as conjugation or GJL). Applying Lemma 5.4, we set the bound B 1 to beB 1 = L Q [2/3, ((1/2) 2 /3) 1/3 ≈ 0.436].The running time of Algorithm 5 will be L q [1/3, (3/2) 1/3 ≈ 1.144]. We obtain preimages P whose pseudonorm is bounded by Q 1-1 2n with the JLSV 1 polynomial selection method as shown in Example 7.1. Applying Lemma 5.4, we set the bound B

		0	. . .	0 
	. . .	. . .	. . .	. . .	 
	. . .		. . .	  0 
	m n 2	. . .	m n 2 , n 2 -1	1

to be the reduced echelon form of L. Then we define the lower triangular matrix made of the n/2 × n/2 identity matrix with p on the diagonal in the upper left quarter, the n/2 × n/2 zero matrix in the upper right quarter, and the n/2 × n matrix M

  11, x + 8 23, x + 15 12239, x + 482 (small)

	1144616018827, x + 218590032699
	2682498999539, x + 1582479651452
	42175797334421, x + 14828919302862
	1195156519724071, x + 966160984838340
	13533793331200309, x + 12224259030902272
	92644276473186311, x + 5754482791048201
	101186915694167857, x + 42826432866764905
	20516170632026633467, x + 14633926248916275064 .

n 2d/n 1 is not necessarily an integer, meaning that the leading coefficient of the polynomial is some element in F p n 1 . Its degree in x is actually n 2 -⌈d/n 1 ⌉.

Acknowledgments

The author is grateful to Francisco Rodríguez-Henríquez, Frederik Vercauteren, Robert Granger and Thorsten Kleinjung, François Morain, Pierrick Gaudry, Laurent Grémy, Luca De Feo, and the other researchers who helped to improve this work. All these very fruitful discussions started at the ECC 2015 conference, the CATREL workshop and the Asiacrypt 2015 conference; in particular, the author would like to thank the anonymous reviewers of Asiacrypt 2015 who suggested a generalization.