Paul Kirchner
email: paul.kirchner@irisa.fr

Thomas Espitau
email: t.espitau@gmail.com

Pierre-Alain Fouque
email: pierre-alain.fouque@irisa.fr

Towards faster polynomial-time lattice reduction

come

Introduction

Lattice reduction and cryptanalysis. Lattice reduction is of the utmost importance in public-key cryptanalysis, as testified, for instance, by the extensive survey of Joux and Stern [START_REF] Joux | Lattice reduction: A toolbox for the cryptanalyst[END_REF]. Indeed, many cryptographic problems are solved by constructing an appropriate lattice and retrieving one of its short vectors. Some standard examples include knapsack problems [START_REF] Lagarias | Knapsack public key cryptosystems and diophantine approximation[END_REF][START_REF] Lagarias | Solving low-density subset sum problems[END_REF][START_REF] Joux | Lattice reduction: A toolbox for the cryptanalyst[END_REF], breaking linear congruential generators [START_REF] Frieze | Linear congruential generators do not produce random sequences[END_REF][START_REF] Stern | Secret linear congruential generators are not cryptographically secure[END_REF]), Coppersmith attack [START_REF] Coppersmith | Small solutions to polynomial equations, and low exponent RSA vulnerabilities[END_REF] against RSA modulus by retrieving small roots of univariate polynomials over Z/N Z or bivariate polynomials over Z, or even attacks against the initial versions of the NTRU cryptosystem [START_REF] Coppersmith | Lattice attacks on NTRU[END_REF][START_REF] Gentry | Key recovery and message attacks on NTRU-composite[END_REF][START_REF] Albrecht | A subfield lattice attack on overstretched NTRU assumptions -cryptanalysis of some FHE and graded encoding schemes[END_REF]]. Yet, its field of applications extends way beyond cryptography, as lattice reduction is a cornerstone of many number theoretical algorithms, allowing factoring polynomials over Z[X] [START_REF] Lenstra | Factoring polynomials with rational coefficients[END_REF], finding integer relations [START_REF] Håstad | Polynomial time algorithms for finding integer relations among real numbers[END_REF], solving simultaneous diophantine approximation problems [START_REF] Lagarias | The computational complexity of simultaneous diophantine approximation problems[END_REF].

Essentially, lattice reduction means finding a short and nearly orthogonal basis to a lattice Λ (represented as a Z-basis). For many applications, finding a small non-zero lattice vector, i.e., solving the (approximate) Short Vector Problem (svp), shall suffice. Since the work of Minkowski, we know that there exists a vector with euclidean norm smaller than √ d(vol Λ)

1 d
, but the proof is not constructive. Nonetheless, the lll algorithm, introduced in 1982 by Lenstra, Lenstra, and Lovász [START_REF] Lenstra | Factoring polynomials with rational coefficients[END_REF] retrieves vector within an exponential factor to the shortest vector of a lattice of dimension d in time O d 6 B 3 where B is the bitsize of the input representation. One can also prove that the norm of the first vector of an lll-reduced basis is less than 4/3 d-1

2 (vol Λ)

1 d
. The approximation factor ( b 1 / vol Λ) 1/d is called the root Hermite factor (RHF), with b 1 a short vector. Later, Schnorr developed a hierarchy of algorithms to reach better RHF in β 1 2β in time 2 O(β) for large β [START_REF] Schnorr | A hierarchy of polynomial time lattice basis reduction algorithms[END_REF][START_REF] Schnorr | Block reduced lattice bases and successive minima[END_REF]. This family leads to a polynomialtime algorithm with a RHF 2 log log d log d [START_REF] Lagarias | Korkin-Zolotarev bases and successive minima of a lattice and its reciprocal lattice[END_REF]. Gama and Nguyen introduced the slide reduction to give an effective take on Mordell's inequality and further improve the RHF [START_REF] Gama | Finding short lattice vectors within Mordell's inequality[END_REF]. In an orthogonal direction, following Haståd and Lagarias, Seysen proposed a variant of lll aiming at simultaneously reduces the primal and dual basis [START_REF] Seysen | Simultaneous reduction of a lattice basis and its reciprocal basis[END_REF]. He defines a new reduceness measure which is closely related to the condition number of the matrix [START_REF] Maze | Some inequalities related to the Seysen measure of a lattice[END_REF].

Related work. The two most singular characteristics of lattices appearing in the cryptographic setting are their high dimension and the large bitsize of their matrix representation. As such, the reduction of cryptanalytically relevant lattices is a computationally intensive challenge. While the original lll implementation works with exact arithmetic on rational entries, Schnorr proposed in 1988 to replace it with floating-point arithmetic [START_REF] Schnorr | A more efficient algorithm for lattice basis reduction[END_REF], significantly improving its efficiency. Since 1996, Shoup maintains a heuristic yet very efficient version in the NTL library with fine control of the float-point precision. This code has been routinely used for more than a decade to break cryptographic schemes. Later, Nguyen and Stehlé precisely analyzed and decreased the asymptotic complexity to O d 5 (d + B)B in [START_REF] Nguyen | An LLL algorithm with quadratic complexity[END_REF], a.k.a. the quadratic lll or L 2 algorithm. This algorithm has been then implemented in fpLLL [START_REF] Albrecht | fplll-5.0, a floating-point LLL implementation[END_REF], which is the current stateof-the-art open-source implementation of lll. However, despite many theoretical improvements to reduce the complexity to quasi-linear in the bitsize using recursive local computation techniques [START_REF] Schönhage | Factorization of univariate integer polynomials by diophantine aproximation and an improved basis reduction algorithm[END_REF][START_REF] Koy | Segment LLL-reduction of lattice bases[END_REF][START_REF] Novocin | An LLL-reduction algorithm with quasilinear time complexity[END_REF][START_REF] Neumaier | Faster LLL-type Reduction of Lattice Bases[END_REF] and some attempts [START_REF] Buchmann | Reducing lattice bases by means of approximations[END_REF][START_REF] Saruchi | LLL reducing with the most significant bits[END_REF][START_REF] Bi | Rounding and chaining LLL: Finding faster small roots of univariate polynomial congruences[END_REF] to use only the most significant bits, the practical complexity of the best implementation available remains in O d 4 B 2 . As such, it struggles to reduce lattices with large entries in high dimensions. Consequently, cryptographers still assess their concrete parameters using L 2 as a reference for lll.

Thus, from a cryptanalytical standpoint, it is interesting to have a fast implementation of lattice reduction (with a controlled approximation factor) even though this algorithm might rely on some heuristics. Since lattice-based cryptography is becoming a strong contender for post-quantum cryptography and offers many interesting functionalities to cryptography, such as efficient Fully Homomorphic Encryption (fhe), new algorithms and implementations of lattice reduction have been designed to give better security estimate for latticebased cryptography. Some improvements mainly target the bkz algorithm since it allows to finely adjust the approximation factor [START_REF] Gama | Finding short lattice vectors within Mordell's inequality[END_REF][START_REF] Gama | Lattice enumeration using extreme pruning[END_REF][START_REF] Hanrot | Analyzing blockwise lattice algorithms using dynamical systems[END_REF][START_REF] Chen | BKZ 2.0: Better lattice security estimates[END_REF][START_REF] Micciancio | Practical, predictable lattice basis reduction[END_REF][START_REF] Albrecht | Faster enumeration-based lattice reduction: Root hermite factor k 1/(2k) time k k/8+o(k)[END_REF]. Others are heuristic and improve sieving technique for solving svp, use a reduction technique in the lattice dimension [START_REF] Albrecht | The General Sieve Kernel and New Records in Lattice Reduction[END_REF], exploit subfield structure and symmetries in structured lattices [START_REF] Pellet-Mary | Approx-SVP in ideal lattices with pre-processing[END_REF][START_REF] Kirchner | Fast reduction of algebraic lattices over cyclotomic fields[END_REF], or use the tensor core architecture of GPU [START_REF] Ducas | Advanced lattice sieving on GPUs, with tensor cores[END_REF]. Some of them with sieving SVP-oracle [START_REF] Albrecht | The General Sieve Kernel and New Records in Lattice Reduction[END_REF] are used to perform the security estimation of signatures and KEM, where the dimension are generally lying between 512 and 1024. However, fhe schemes over the integers use extremely large integers (several millions of bits) and high dimensional lattices (typically of a few thousand dimensions), but can be broken with high approximation factors. To deal with such settings, faster algorithms are required, in particular with complexity quasi-linear in the bitsize and not much more costly than matrix multiplication.

Our Contributions. To improve the running time of lattice reduction algorithms, we propose to exploit parallelism with many cores, make full use of computer's cache using block matrix implementation [START_REF] Golub | Matrix Computations[END_REF], and use a low precision while still controlling the approximation factor at the same time. Our implementation we describe allows reducing lattice in dimensions up to 2,000 with entries of up to millions of bits, which is intractable otherwise.

Our proposal of lattice reduction is a lll-type algorithm, i.e., using a sizereduction procedure jointly, together with many passes of a rank-2 reduction subprocess. The design rationale is to exploit fast block matrix operations and locality of operations. To do so, we use a block variant of the Cholesky factorization algorithm for computing the QR-decomposition [START_REF] Golub | Matrix Computations[END_REF]. We replace the sizereduction with a block variant of Seysen's size-reduction, which can be thought of as a rounded version of the multiplication by the inverse of the R factor of the QR-decomposition. To our knowledge, this algorithm has not been used since 1993. Contrary to the textbook lll, we do not swap vectors when the Lovász condition is not fulfilled, but we fully reduce the 2-dimensional corresponding projected sublattice, using Schönhage's algorithm. The global design is recursive, as was proposed before by Koy and Schnorr [START_REF] Koy | Segment LLL-reduction of lattice bases[END_REF] with Segmented lll and by Neumaier and Stehlé [START_REF] Neumaier | Faster LLL-type Reduction of Lattice Bases[END_REF]. However, in this work, we do not recurse on overlapping blocks but on separate ones; a technique proposed by Villard to achieve parallelism [START_REF] Villard | Parallel Lattice Basis Reduction[END_REF] with even and odd steps, also used recently in [START_REF] Kirchner | Fast reduction of algebraic lattices over cyclotomic fields[END_REF].

As all the computations are conducted in floating-point arithmetic, a systematic caveat concerns the precision required for computing the correct result. We claim and experimentally verify that on average, it decreases exponentially with the recursion depth as shown in section 4.1, allowing to reduce the overall complexity by a factor d. Additionally, we handle matrice multiplications in the Fourier domain to compute with large numbers. We conjecture and experimentally verify a complexity of approximately d ω • C/ log C, where ω is the exponent of matrix multiplication, and C is the logarithm of the condition number of the input matrix. We highlight that typically, the complexity of lattice reduction depends on the bitlength B of the input, instead of C. For cryptographic applications (Coppersmith and knapsack-type lattice amid others), C is close to B, while it can be up to d times larger in the worst case. It is well-known that a row-wise diagonal dominant matrice has a condition number bounded by a constant times the ratio beteen the largest diagonal entry and the smallest one, so the logarithm of the condition number will be close to the size of the entries.

Additionally, section 5 shows that one can reduce knapsack lattice in a time approximately equal to the reduction of a random lattice with a bitsize reduced by a factor of d. Such a phenomenon is already known for some algorithms like fplll [68, 1.5.3], noted [START_REF] Nguyen | LLL on the average[END_REF], and exploited [START_REF] Van Hoeij | Gradual sub-lattice reduction and a new complexity for factoring polynomials[END_REF]. We present a reduction between the two problems with this property. The idea is to iteratively double the number of columns reduced, and reduce the bitsize of the other ones. It has been implemented and tested.

The complexity of our algorithms can be analyzed in an arithmetic cost model with an analysis similar to [START_REF] Hanrot | Analyzing blockwise lattice algorithms using dynamical systems[END_REF] (sandpile model) for LLL with even and odd pass as in [START_REF] Kirchner | Algebraic and Euclidean Lattices: Optimal Lattice Reduction and Beyond[END_REF]. However, the specificity of our algorithm is to consider the precision. Without such attention, it would have been impossible to reduce high dimensional lattice with so many bits. In an exact arithmetic cost model, the complexity would have been comparable to previous algorithms, which is not the case in practice. Such heuristic algorithms is interesting, even without a full analysis, to assess the security of cryptographic instances. For instance, many fhe schemes over the integers base their security on the complexity of the best algorithm. However, a rigorous proof of the algorithm with the precision is highly technical in a numerical computational model and escape us so far. Consequently, we decided to present only all the ingredients of our implementation with its applications in this paper and postpone a proof for future work.

Regarding the applications, we first show in section 6 that our implementation is much faster than fplll with a factor between 30 and 45 on single-thread in all dimensions tractable by fplll. However, our implementation can exploit multicore processors and reduce lattices in much higher dimensions. Consequently, we run it on matrices of dimensions a few thousand and inputs of millions of bits, as reported in table 1. As a result, we attack many instances fhe over the integers to illustrate the efficiency of our code and evaluate its running time on large inputs. For these examples, the wall-clock time is six orders of magnitude smaller than the (estimated) cost of fplll. We broke knapsack instances from [START_REF] Coron | Practical multilinear maps over the integers[END_REF] in dimension 2,230 with 4.26 millions of bits in 22h with 18 cores, while the security level was evaluated at 2 62 . We also broke ntru instances with overstretched parameters proposed in [START_REF] Genise | Homomorphic encryption for finite automata[END_REF] in 5h (resp. 10 days) in dimension 2560 (resp. 3086) with 111 (resp. 883) bits and RHF 2 0.1105 (2 0.018 , equivalent to BKZ-25). In practice, at the bottom of the recursion tree, we use a small bkz to improve the approximation factor, whilst not altering too much the running time.

Background

Notations and conventions

The capitals Z, Q, R refer to the ring of integers, the field of rational and real. Given a real number x, its integral rounding denoted by x returns its closest integer. The logarithms are log for the binary one and ln for the natural one.

Matrix and norms. We denote by Q d×d the space of square matrices of size d over Q, GL d (Q) its group of invertible. We use bold fonts for matrices and denote the elementary matrix transformations by T i,j (λ) and D i (λ) for respectively the transvection (or shear mapping) and the dilatation of parameter λ. We use Diag(x 1 , . . . , x d ) to refers to a diagonal matrix of elements x 1 , . . . , x d . We generalize this definition to block matrices and overload it to the extraction of the diagonal of a given matrix. A triangular unipotent or unitriangular matrix is a triangular matrix with ones on the diagonal. We extend the product for any pair of matrices (A, B): for every matrix C with compatible size with A and B, we set: (A, B)•C = (AC, BC). We adopt the usual conventions for submatrix extraction: for any matrix M = (m i,j ) ∈ Q d×d and 1 u < v d, 1 w < x d, define the submatrix M[u : v, w : x] = (m i,j ) u i v,w j x , while M i refers to the i-th column of M. For a vector v (resp. matrix A = (a i,j ) 1 i,j d ), we denote by v (resp. A ) the Frobenius norm, i.e., A = 1 i,j d a 2 i,j . The condition number of an invertible matrix M measures how much the output value of the matrix can change for a small change in the input. It is defined as κ(M) = M M -1 and allows to compute the precision needed during the computation. We deal with block decomposition of matrices, with block of halfdimension. For matrices of odd dimension 2k + 1, the upper-left block to be of dimension k + 1 and the bottom-right one of dimension k.

Computational setting. We use the standard model in algorithmic theory, i.e., the word-RAM with unit cost and logarithmic size register (see [52, Section 2.2] for a comprehensive description). The number of bits in the register is w and the precision during the computation by p. All computations with rational/real values are conducted in floating-point, unless stated otherwise. For a non-negative integer d, we set ω(d) to be the exponent of matrix multiplication of d × d matrices. If the dimension d is clear from context we might omit it and write simply O(d ω ) for this complexity. We can assume that this exponent is not too close to 2, in particular ω(d) > 2 + 1/ log(d). Due to the conflict with Laudau's small omega notation, we use ω for the latter symbol.

Lattices and LLL reduction Definition 1 (Lattice

). A d-dimensional (real) lattice Λ ⊆ R d is the set of in- teger linear combinations d i=1 b i Z of some linearly independent vectors (b i ) 1 i d . The finite family (b 1 , . . . , b d ) ∈ Λ is called a basis of Λ.
Every basis has the same number of elements called the rank of the lattice. A measure of the density of the lattice is its (co)volume, defined to be the volume of the torus R d /Λ, which corresponds to the square root of the Gram-determinant of any basis (b 1 , . . . , b d ):

vol Λ = det( b i , b j ) 1 i,j d .
Two different bases of a lattice Λ are related by a unimodular transformation, i.e., a linear transformation represented by an element of GL d (Z), the set of d × d integer-valued matrices of determinant ±1. In essence, algorithms acting on lattice bases are sequences of unimodular transformations. Among these procedures, reduction algorithms are of the utmost importance. They aim at finding congenial classes of bases, which are quasi-orthogonal and with controlled norms. Fundamental constant associated to any rank d lattice Λ are its successive minima λ 1 , . . . , λ d . The ith minimum λ i (Λ) is the radius of the smallest sphere centered in the origin containing i linearly independent lattice vectors.

Orthogonalization, QR-decomposition. Let B = (b 1 , . . . , b d ) a family of linearly independent vectors. Let π i the orthogonal projection on (b 1 , . . . , b i-1 ) ⊥ , with the convention that π 1 = Id. The Gram-Schmidt orthogonalization process is an algorithmic method for orthogonalizing B while preserving the increasing chain of subspaces

( i j=1 b j R) 1 i d . It constructs the orthogonal set B * = (π 1 (b 1 ), . . . , π d (b d )).
For notational simplicity we refer generically to the orthogonalized vectors by b * i for π i (b i ). The computation of B * can be done inductively as follows: for all

1 i d, b * i = b i - i-1 j=1 bi,b * j b * j ,b * j b * j .
Collect the family B in a matrix also denoted by the same notation and set

R i,j = bj ,b * i b * i and Q = b * 1 b * 1 . . . b * d b * d
. Then, we have B = QR, with Q being an orthogonal matrix and R being upper triangular. This is the QR-decomposition of B. In the following, we work with the R part only, so that we present the computation of this matrix in the pseudo-code Orthogonalize below. We omit considerations on the required fp-precision here, to just focus on the core ideas of the algorithms. the input matrix B is of dimension d×d with B-bit entries. This process is called Size-reduction and corresponds to the following iterative algorithm3 Size-reduce.

Q i ← b i - b i ,Q j Q j ,Q j Q j end for end for return R = Q i ,b j Q i 1 i j d Algorithm 2 -Size-Reduce Input : Basis B, R part of qr-decomposition Output : Tranformation of SR basis U = Id d for i = 1 to d do for j = i -1 to 1 do (U, R) ← (U, R) • T i,j -R[i,j] R[i,i]

The LLL reduction algorithm

Lenstra, Lenstra, and Lovász [START_REF] Lenstra | Factoring polynomials with rational coefficients[END_REF] proposed a notion called lll-reduction and a polynomial-time algorithm that computes an lll-reduced basis from an arbitrary basis of the same lattice. Their reduction notion is formally defined as follows (presented directly in an algorithmic way with the QR-decomposition):

Definition 2 (LLL reduction). A basis B of a lattice, admitting the decomposition B = QR, is said to be δ-lll-reduced for 1/4 < δ 1, if the following two conditions are satisfied:

∀i < j, |R[i, j]| 1 2 |R[i, i]| (Size-Reduction condition) (1) ∀i, δ R j,j 0 2 R j,j+1 R j+1,j+1 2
(Lovász condition). (2)
The length of vectors and orthogonality defect is related to the parameter δ:

Proposition 1. Let 1/4 < δ 1 be an admissible lll parameter. Let (b 1 , . . . , b d ) a δ-lll reduced basis of rank-d lattice Λ. Then for any 1 k d: vol (b 1 , . . . , b k ) (δ -1/4) -(d-k)k 4 vol Λ k d .
In particular, we have that R i,i (δ -1/4) -1 R i+1,i+1 .

We recall that vol(Λ) = det(B) = d i=1 R i,i and the log-potential is defined as

Π(B) = d i=1 (d -i) log(R i,i
). For k = 1 and δ = 1, the Hermite approximation factor defined as b 1 / det(B) 1/d , will be (4/3) (d-1)/4 . To find a basis entailing the lll conditions, it suffices to iteratively modify it at any index violating one of these conditions. This process yields the simplest version of the lll algorithm. However, we choose to present a different take on this algorithm, closer to the algorithms we introduce later. The first remark is that for a given 1 j d -1, the lll-reduceness conditions correspond to saying that the basis

R j,j R j,j+1 0 R j+1,j+1
is Gauss-reduced. The global strategy given in algorithm 3 to reduce a lattice consists of iteratively applying a reduction procedure in rank 2 to projected sublattices, naturally using the Gauss reduction algorithm [START_REF] Hanrot | Analyzing blockwise lattice algorithms using dynamical systems[END_REF]. We start by reducing the sublattice spanned by b 1 , b 2 , then the projection onto the orthogonal subspace to b 1 sublattice spanned by b 2 , b 3 and so on. When we hit the end of the basis, this iteration restarts afresh until no more progress is achieved. We replace the outermost while loop by a for loop of a fixed number ρ of iterations. This parameter is set to be sufficiently large to ensure the reducedness of the output (using a dynamical system analysis à la [START_REF] Hanrot | Analyzing blockwise lattice algorithms using dynamical systems[END_REF] after O d 2 log B rounds, a vector within the lll quality bound is discovered).

We will use a slight generalization of the lll-reduction notion. In particular, a lll-reduced basis satisfying the Lovász conditions, is a Siegel reduced basis.

Definition 3 (Siegel reduction). The Siegel reduction problem consists in, given an integer matrix A of dimension d with

A , A -1
2 B , outputting a matrix AU with U a unimodular integer matrix such that with QR = AU the QR-decomposition, we have for all i: R i,i 2R i+1,i+1 .

Algorithm 3 -Reduction

Input

:

Initial basis B = (b 1 , . . . , b d ) Output : A δ-lll-reduced basis 1 while B is not lll-reduced do 2 R ← Orthogonalize(B) 3 U i ← Size-Reduce(R) 4 (B, R) ← (B, R) • U i 5 for j = 1 to d do 6 B ← R[j : j + 1, j : j + 1] 7 U ←Gauss(B ) 8 (U i , B) ← (U i , B) • Diag(Id j-1 , U , Id d-j-1 )
9 end for 10 end while 11 return ρ i=1 U i // ρ is the number of passes

Matrices Representation

A matrix A is represented as A 2 e where A is an integer matrix and e 0. The quantity log( A ) is the precision of the matrix. The standard algorithm for multiplying matrices with large entries consists in transforming the integers in A and B into polynomials of degree bounded by O p+w w (p is the precision and w the number of bits in registers), and computing their evaluations on roots of unity. The matrices of evaluations are then multiplied, and an inverse Fourier transform gives the product of the matrix of polynomials. Carries are then computed to obtain AB. Matrices can be multiplied quickly using the FFT:

Fast inversion of unitriangular matrices

We eventually conclude this preliminary section by introducing a natural recursive algorithm to invert unitriangular matrices-working with floating-point approximation. It is a direct application of the computation of Schur's complement in the case of a block triangular matrix, i.e., the observation that:

A C 0 D -1 = A -1 -A -1 CD -1 0 D -1 .
As both A and D are unitriangular, this inversion formula translates naturally in a recursive algorithm. Its base case corresponds to inverting a one dimensional unitriangular matrix, that is (1), which is its own inverse. The corresponding pseudo-code is given in Invert. Its complexity is easily analyzed to be asymptotically the cost of a matrix multiplication, as the dominant step of each recursive call is the computation of the complement -A -1 CD -1 .

Algorithm 4 -Invert

Input : A unitriangular matrix M Output : A fp-approximation of M -1 1 if dim(M) = 1 then 2 return (1)
3 end if We provide the precise analysis of this inversion. It can be extended to triangular matrices, and we consider that Invert also computes their inverse.

M -1 -M 2 √ d2 -p .
The case d = 1 is straightforward, so that we now deal with inductive case d > 1. Let E, δA and δD be matrices such that the top-right part of M is -A CD + E, A -1 = A + δA, and D -1 = D + δD. Consequently, we get:

M -1 -M = δA -A -1 ED -1 0 δD
. We can guarantee that E 2 -p -2p with a computation with intermediary bitsize O(p ). This leads to our intermediary result. Now let

M -1 = M + F, so M = (M(Id + M -1 F)) -1 = (Id + M -1 F) -1 M -1 and M -M -1 M -1 (Id + M -1 F) -1 -Id 2 -p .
The complexity comes from the matrix multiplication with words of size w.

Fast reduction of Euclidean lattices

This section is devoted to the description of our block recursive lattice reduction algorithm. In the following, let us fix a Euclidean lattice Λ of rank d, described by a basis collected in a rational matrix B in the canonical basis of R d . We generically denote by R the R-part of the QR-decomposition of this matrix. We recall that computations are conducted in floating-point arithmetic. However, for the sake of readability and ease of presentation, we defer the issue of the necessary precision to section 4.

We turn to a detailed breakdown of the essential parts of the algorithm. Each of the following subsections details and refers to the corresponding lines of algorithm 5, Reduce.

1 if d = 2 then return Schonhage(B) 2 for i = 1 to ρ do 3 R ← Block-Cholesky(B T B) 4 U i ← Size-Reduce(Diag(R) -1 • R) 5 (B, R) ← (B, R) • U i 6 for j = 1 + (i mod 2) to D/2 by step of 2 do 7 V 1 ← vol(R[(j -1)d + 1 : jd , (j -1)d + 1 : jd ]) 8 V 2 ← vol(R[jd + 1 : (j + 1)d -1, jd + 1 : (j + 1)d ]) 9 if V 1 2 2(1+ε)α(d ) 2 V 2 then 10 U ← Reduce(R[jd : (j + 2)d -1, jd : (j + 2)d -1]) 11 (U i , B) ← (U i , B) • Diag(Id jd , U , Id d-3jd ) 12 end if 13 end for 14 end for 15 return ρ i=1 U i //
The product is computed from the end

Base case: plane lattices [Line 1]

As in all variants of the lll algorithm, the base case of the reduction boils down to the two-dimensional case, usually handled by the celebrated Lagrange-Gauss reduction or some equivalent transformations. For instance, in the original lll algorithm, truncated steps of Lagrange-Gauss reduction are conducted on twodimensional projections of shape

π i (b i )Z ⊕ π i (b i+1 )Z.
For the sake of efficiency, we adapt Schönhage's algorithm [START_REF] Schönhage | Fast Reduction and Composition of Binary Quadratic Forms[END_REF], as in the algorithms of [START_REF] Kirchner | Fast reduction of algebraic lattices over cyclotomic fields[END_REF][START_REF] Hanrot | Analyzing blockwise lattice algorithms using dynamical systems[END_REF], to reduce these plane lattices. This algorithm is an extension to the bidimensional case of the so-called half-GCD algorithm [START_REF] Möller | On Schönhage's Algorithm and Subquadratic Integer GCD Computation[END_REF], likewise that Gauss' algorithm is a bidimensional generalization of the classical Euclid's gcd. The original algorithm of Schönhage only deals with the reduction of binary quadratic forms but can be straightforwardly adapted to reduce lattices, as well as returning the corresponding unimodular transformation matrix. In the following, we denote by Schonhage this modified procedure. Its complexity is quasilinear in the size of its input (which is to be compared with the quadratic complexity of the classical Gauss reduction).

Outer iteration [Line 2]

To reduce the lattice Λ, we adopt an iterative strategy to progressively modify the basis: for ρ > 0 steps, a reduction pass over the current basis is performed, ρ being a parameter set to optimize the complexity of the whole algorithm while still ensuring the reduceness of the basis. We defer the choice of this constant for the moment. This global iterative scheme is similar to the terminating variants of the bkz algorithm, for instance as in [START_REF] Hanrot | Terminating BKZ[END_REF] or [START_REF] Micciancio | Practical, predictable lattice basis reduction[END_REF], where a polynomial number of rounds is fixed to reduce the input.

Orthogonalization via Block-Cholesky decomposition [Line 6]

Gram-Schmidt Orthogonalization is a preliminary step of every lll-type algorithms, as it computes the so-called Gram-Schmidt vectors of the basis, which are ubiquitous in the definition of the reduction itself. On symmetric matrices as the Gram-matrix B T B of the basis, one computes the Cholesky factorization, which given a symmetric positive-definite matrix G, the factorization asserts the existence (and unicity) of an upper triangular matrix R such that

G = R T R which is the some R in the QR decomposition of B since G = B T B = R T R.
We use here a recursive block variant of the Cholesky factorization algorithm, allowing to compute a floating-point approximation of the matrix R, whose running time is heuristically the cost of a matrix multiplication. It relies heavily on the Invert procedure introduced in Section 2.5.

Remark 1. Block computations of decompositions seems to be folklore in numerical algebra (see, for instance, the complete monograph of Higham [START_REF] Higham | Accuracy and stability of numerical algorithms[END_REF] for multiple variants of block orthogonalization, such as modified Gram-Schmidt, Householder transformations, . . . ), but oddly, we were unable to find a proper reference to the block Cholesky factorization.

The decomposition is as follows, given as input a symmetric matrix G. We start by block splitting it (with blocks of half size):

G = A B B T C , where
A, C are also symmetric. Its Schur complement S = C -B T A -1 B is then also symmetric. Suppose that we know the factorization of the A and S in say:

A = R T A R A and S = R T S R S . Then, we set R = R A R -T A B 0 R S .
This matrix is indeed the Cholesky factorization of G, as ensured by the following computation:

R T R = R T A 0 B T R -1 A R T S • R A R -T A B 0 R S = R T A R A R T A R -T A B B T R -1 A R A B T R -1 A R -T A B + R T S R S = A B B T C , since B T R -1 A R -T A B + R T S R S = B T A -1 B + C -B T A -1 B = C by definition of the Schur complement.
This derivation yields a direct recursive algorithm, whose base case corresponds to the unidimensional instance, i.e., G = (g), admitting the trivial decomposition G = ( √ g) T ( √ g). This observation yields the procedure stated in pseudocode in algorithm 6 Block-Cholesky, computing a floating-point approximation of the Cholesky decomposition.

Algorithm 6 -Block-Cholesky

Input : A positive-definite symmetric matrix G Output : A fp-approx. of a triangular matrix R s.t.

R T R = G 1 if dim(G) = 1 then return √ G 2 A B B T C ← G // with blocks of half-dimension 3 R A ← Block-Cholesky(A) 4 R A ← Invert(R A ) 5 A ← R T A R A 6 R S ← Block-Cholesky C -B T A B 7 return R A R T A B 0 R S

Size-reduction [Line 4]

As in the lll algorithm, a size-reduction operation is conducted at each step of the reduction. It allows to control the size of the coefficients and ensures that the running time remains polynomial. However, in our case, we lean on a Seysen-like reduction to perform this operation [START_REF] Seysen | Simultaneous reduction of a lattice basis and its reciprocal basis[END_REF]. Our recursive procedure allows to sizereduce a unitriangular matrix (in our case, the matrix Diag(R) -1 R) in roughly the time of matrix multiplication.

Algorithm 7 -Size-Reduce Input : A unitriangular matrix T Output : An integer unitriangular matrix U, TU reduced

1 if dim(T) = 1 then return (1) 2 A C 0 D ← R // with half dimension 3 U 1 ← Size-Reduce(A) 4 U 2 ← Size-Reduce(D) 5 A ← Invert(AU 1 ) 6 W ← A CU 2 7 return U 1 -U 1 W 0 U 2
We start from the classical observation that the usual size-reduction process is a discretized version of the iterative Gram-Schmidt process (which is a way of computing the QRdecomposition of a matrix). Over the triangular matrix R, it corresponds to make iteratively the extra diagonal elements as close as possible to 0. However, instead of using an iterative process, we use a lattice reduction algorithm with block matrix operations.

Let us start with a unitriangular matrix R, split in block of half dimension: A C 0 D

. Assume for the moment that both unitriangular submatrices A and D are already size-reduced. Then, set

U = Id -A -1 C 0 Id ,
which is unimodular as its diagonal elements are all 1. Its action on R gives by elementary computation:

RU = A C -A A -1 C 0 D
and the top-right part is of the same magnitude as A. The inverse of RU is

A -1 -A -1 C -A -1 C D -1 0 D -1
and the top-right part is of the same magnitude as D -1 , ensuring that the norm of this block is controlled. The translation of this process in pseudocode yields algorithm 7 Size-reduce. Note that this algorithm is presented as yielding the transformation matrix instead of the reduced matrix, to be consistent with the presentation of Reduce (see proof in Appendix A).

Step reduction subroutine [Lines 3-13]

Fig. 1: Illustration of the parallel step reduction on the R-part of the QRdecomposition. Green 2 × 2 blocks are simultaneously reduced on odd steps and orange ones are reduced on even steps. This strategy is similar to [START_REF] Heckler | Complexity analysis of a parallel lattice basis reduction algorithm[END_REF].

From parallel design of lll... Let us now describe the step reduction pass, occuring once the size-reduction operation has been performed. As observed in section 2, the lll algorithm reduces lattice reduction to the reduction of rank two lattices (more precisely, iteratively reduce orthogonally projected rank-2 sublattices). A first idea would be to use the same paradigm here and pass over the current basis in a sequence of reduction of projected planar lattices. However, on the contrary to the standard lll or bkz-2 algorithms, remark that we are not forced to proceed progressively along the basis, but that we can reduce d/2 independent (nonoverlapping) rank-2 lattices at each step, namely the

(π 2i (b 2i Z ⊕ b 2i+1 Z)) 1 i d/2 and then, (π 2i+1 (b 2i+1 Z ⊕ b 2i+2 Z)) 0 i d/2 .
This design enables an efficient parallel implementation which reduces sublattices simultaneously, in the same way that the classical lll algorithm can be parallelized [START_REF] Villard | Parallel Lattice Basis Reduction[END_REF][START_REF] Heckler | Complexity analysis of a parallel lattice basis reduction algorithm[END_REF]. This technique can also be thought of as a parallelized bkz [START_REF] Micciancio | Practical, predictable lattice basis reduction[END_REF] or slide-reduction [START_REF] Aggarwal | Slide reduction, revisited -filling the gaps in SVP approximation[END_REF] with blocksize 2. ...to recursive block design A bottleneck with this strategy is that each round needs (at least) a matrix multiplication to be updated. Using a dynamical system analysis similar to [START_REF] Hanrot | Analyzing blockwise lattice algorithms using dynamical systems[END_REF], such a reduction would require ρ rounds to be a Ω(d 2 ) to ensure an lll approximation factor. This implies a dependency in the running time which would be at least quartic in the dimension d. However, one can notice that each round only makes local modifications on the basis. As a result, we propose to use a small number D of blocks, and let a round recursively reduces consecutive pairs of blocks of dimension d D . In this setting, the dynamical system analysis of [START_REF] Hanrot | Analyzing blockwise lattice algorithms using dynamical systems[END_REF] shows that a O D to enable the reduction of cross blocks. This step reduction is then restarted for ρ rounds as indicated in section 3.2.

On the volumetric Siegel condition. Remark the use of relaxation parameters ε, α > 0, acting on the approximation factor of the reduction. As an avatar of the relaxation factor δ of lll, they allow a slight tradeoff between the running time and the overall reduction quality. It is an equivalent of the Siegel condition between blocks: instead of recursively calling the reduction every time on the blocks R 2j ⊕ R 2j+1 , we only do it if the volume of the left block R 2j is sufficiently larger than the one of the right block R 2j+1 . We do not perform a recursive reduction if the slope between the blocks is already small enough.

In practice, these values are dependent on the depth of recursion to optimize the global running time. Section 4 addresses this technicality more thoroughly.

Complexity estimation and supporting experiments

We now turn to the fine-tuning of the implementation and describe some optimization tricks used. We backed up our choices by supporting experiments and eventually devise an empirical estimate of the bit-complexity of our algorithm.

Needed Precision

(max i R[i, i]/ min i R[i, i]).
As predicted by heuristic 1 the corresponding graph presents an exponential decay.

Since the implementation of the algorithm is done using floating-point arithmetic, we need to set a precision which is sufficient to handle the internal values during the computation. To do so, we set:

p = log max i R[i, i] min i R[i, i] ,
where the R[i, i] encodes the norm of the Gram-Schmidt vectors. As in floatingpoint variants of lll [START_REF] Schnorr | A more efficient algorithm for lattice basis reduction[END_REF][START_REF] Nguyen | An LLL algorithm with quadratic complexity[END_REF][START_REF] Koy | Segment LLL-reduction of lattice bases[END_REF][START_REF] Neumaier | Faster LLL-type Reduction of Lattice Bases[END_REF], it is straightforward that a O(p) is sufficient to handle the computation. However, the remaining question is the evolution of this quantity within the recursive calls. Indeed, as we have more and more recursive calls of the reduction algorithm on projected lattices of smaller dimensions, we would like to reduce them with a limited precision to get an overall faster reduction. The analysis of [START_REF] Neumaier | Faster LLL-type Reduction of Lattice Bases[END_REF] bounds the number of rounds, and reaches a complexity in d 3 C 1+o (1) with exact arithmetic (C is the log of the condition number), while the non-optimized algorithm [START_REF] Hanrot | Analyzing blockwise lattice algorithms using dynamical systems[END_REF] uses Ω(d 3 log C) local reductions. Consequently, to decrease the complexity of the reduction, we have to reduce the precision in the local operations. The justification of this fact comes from that in practice, the values of R[i, i] decrease roughly exponentially in i both in the input and the output matrices. To define our heuristic, we rely on the notion of slope of the basis, which is the opposite of the slope of the linear regression of the log of the norm the Gram-Schmidt vectors. Under the Geometric Series Assumption (GSA), this corresponds to the usual geometric decay factor. Heuristic 1 says that we reduce the slope, i.e. the logarithm potential

Π(B) = d i=1 (d -i) log(R i,i
). We consider that we have access to an oracle which reduces with a slope parameter of α, namely R[i, i]/R[i + 1, i + 1] ≈ 2 2α . The matrix returned will have a slope parameter of (1 + ε)α for 0 < ε < 1/2. Heuristic 1 If ρ is even, and ρ 2 (2D -1) D 3 , then the slope decreases exponentially quickly towards

(1 + ε)α, with rate 1 -O 1 D 2 . Remark 2.
For a smaller ρ, we would have several leaves in the recursion tree, which would be negligible compared to d 3 , making it unlikely to reduce the lattice by a significant amount. These values come from an heuristic analysis.

Figure 3 shows the evolution of the slope on a lattice of dimension 1024 where the phenomenon is observable. heuristic 1 has been tested on various types of lattices (Knapsacks, NTRU-like) in dimensions from 128 to 2048 without failing.

On the choice of the relaxation parameter ε and its relation to the global complexity

To finely tune our parameters, we need to estimate the decrease of the potential at each recursive call. Using heuristic 1 at any moment in the recursion, when called with a lattice of rank d and working precision p = log maxi Ri,i mini Ri,i , such that

d/2 i=1 R i,i > 2 (1+ε)αd 2 /2 d i=d/2+1 R i,i , (condition)
the output basis has a log-potential reduced by at least Ω(d 2 pε). Calling a recursive reduction only when the condition is fulfilled allows the callee to reduce the slope by a factor of roughly 1+ε. If this is actually done, the potential is reduced by Ω(ε d D 2 p ) where p is the precision used by the callee. The complexity of the callee, if not already in a leaf, and outside of its recursive calls is in

O D 2 (d/D) ω (p /w + 1) + (d/D) 2 p log p w .
Keeping only the first term and assuming p > w, we get that the complexity per unit reduction in potential should behave in

O D 2 (d/D) ω-2 w -1 ε -1 .
This suggests minimizing D, so that we set D = 2 and ρ = 6; and also we deduce that most of the complexity is at low-depth. While the global complexity is minimized for D = 2, considering a larger D leads to better running time when using multithreading (higher number of blocks can be treated in parallel).

If we write d i and ε i for their values at depth i, we obtain that the global approximation factor is the one at the leaf multiplied by exp( i ε i ). Also, the main term in the complexity is proportional to i d ω-2 i ε -1

i . Thus, we want ε i proportional to d 1-ω/2 i . If we want i ε i = Θ(δ), we get

ε i = δ(ω(d) -2)(d/d i ) 1-ω(di)/2 .
Summing the complexity at all depths, we see that the main term becomes:

O d ω C w(ω -2) 2 δ for any δ = O 1 ω -2 .

Using small-dimension fast enumeration in the leaves

Since almost all the complexity concentrates at low recursive depth, we can allocate more time in the leaves of the recursion tree to improve the quality of the reduction without altering much the global complexity. In practice, this means stopping the recursion before reaching rank-2 sublattices and using a stronger reduction process than lll on these (higher dimensional) leaves. Some instances of stronger algorithms are the bkz-type family, which are parameterized by a block size β, and have a complexity exponential in β [START_REF] Becker | New directions in nearest neighbor searching with applications to lattice sieving[END_REF]. This family includes Schnorr's original bkz algorithm, Terminated-bkz with less rounds [START_REF] Hanrot | Analyzing blockwise lattice algorithms using dynamical systems[END_REF], the self-dual bkz [START_REF] Micciancio | Practical, predictable lattice basis reduction[END_REF] or pressed-bkz [START_REF] Bai | Measuring, simulating and exploiting the head concavity phenomenon in BKZ[END_REF]-which is particularly good for low β. If the dimension at the leaf d l is significantly larger than β log β, the famous Geometric Series Assumption states that the Gram-Schmidt norms of the reduced basis are well approximated by a geometric series of rate 2 Θ( β log β ) . We can assume that the basis was already reduced with a constant slope 2α, so that the potential will overall decrease only by O d 3 . At each leaf, we can use a constant ε l-1 and thus expect the log-potential to decrease by at least

Ω(d 3 l log β β ). The number of calls is therefore O d 3 β d 3 l log β = O d 3 β 2
so we can choose any β smaller than Ω((ω -2) log C).

Complexity estimation

The sketch of analysis conducted previously let us conjecture that the complexity should have a dominant term in d ω C. We plot the single-thread running time on lattices with dimension d = 2n generated by the columns of the following matrix

qId n A 0 Id n
with A sampled uniformly modulo q, and C = log(q) ≈ n4 k-1 for k from 0 (green) to 3 (blue). The slope of the reduced matrix is 2α ≈ 0.065 (RHF= 2 0.032 = 1.02). To confirm this hypothesis, we perform a linear regression on the log/log data of the running time in function of the input dimension (ranging from 128 to 2048). The regression reveals a slope of 3.5, that is a complexity in O d 2.5 C as C is linear in d. Given the noise generated by the inherent complexity of the program, its libraries, and the complex processor architecture, this experiment seems to validate our conjectural complexity. Each line corresponds to experiments made with matrices with bitsize bounded by (dimension K) [from green (lower) line with K = 1/4 to blue (upper) line K = 16]. We propose the following complexity for our algorithm, using the small bkz-enumeration in the leaves.

Analysis 1 Let A be a matrix of dimension d with integer entries, with κ(A) 2 C such that C d/ log d. Reduce(A) returns the transformation matrix to a basis of AZ d having its first vector of norm bounded by

max √ d, 2 O(d(ω-2) log log C log C ) vol A 1 d
Further, the heuristic running time is

O d ω • C (ω -2) 2 log C + d 2 C log C + d 2 C (ω -2) 2 .
Remark 3. The values come from a heuristic analysis that we do not develop.

• In practice, the entire basis is reduced at the end of the algorithm (as lll algorithm gives a reduced basis with controlled decay of the Gram-Schmidt).

• When ω is bounded away from 2, and C is not extremely large (C = 2 o(d) ), the complexity simplifies to O d ω • C log C . • It is better to first reduce with a large δ (say min(log(C/d), 1 ω-2 )), and progressively reduce the slope by decreasing δ by a constant, so that the precision used is exponentially decreasing. For C > d2 1/(ω-2) , we obtain a heuristic complexity of:

O d ω • C (ω -2) log C + d 2 C log C .
• The dependency in the second term of the complexity (term in d 2 C log C) comes as a direct consequence of the complexity of the Schonhäge algorithm.

The implementation mixes multiple machine representation as it needs to manage efficiently both large and small matrices, with a large range in bitsizes. On the one hand, the "large matrices", e.g. with of dimension greater than 80 and of coefficient represented on few hundreds bits, are represented in the Fourier domain, that is to say by a collection of complex matrices, one for each evaluation point. The complex matrices are with double-precision floating-point coordinates. Large integers are transformed into polynomials, with between 14 and 16 bits per coefficient.

On the other hand, small matrices (dimension lower than 80) and with small bitsize are represented with an array of MPFR values [START_REF] Fousse | MPFR: A multiple-precision binary floating-point library with correct rounding[END_REF]. A reduction of small matrix with at most 300 bits is computed by repeatedly reducing matrices with at most 39 bits, which are in turned reduced using blocks of dimension 12. These matrices of dimension 12 and with at most 20 bits are reduced with the quadratic L 2 [START_REF] Nguyen | An LLL algorithm with quadratic complexity[END_REF] procedure.

Finally, matrices where p is small (around 30) and dimension up to 400 are treated in double precision, thanks to the use of the Householder QRdecomposition and the Seysen size-reduction.

Reduction of structured knapsack-like

In this section, we present a progressive strategy to provably speed-up the reduction of almost triangular matrices. Combined with the reduction of section 3, it gives a heuristic reduction process which estimated running time is essentially a O d ω-1 C log C + Cd log d . The general idea is that a knapsack-like matrix of dimension d and with log condition number C can be reduced as quickly as a matrix of dimension d and condition number 2 C/d . As this effect was already known for some algorithms like fplll [68, 1.5.3], noted [START_REF] Nguyen | LLL on the average[END_REF], and used in [START_REF] Van Hoeij | Gradual sub-lattice reduction and a new complexity for factoring polynomials[END_REF], we aim at giving a general framework to encompass this observation.

Setting

Definition 4 (Almost triangular matrix). A matrix B with d columns and O(d) rows is said to be (asymptotically) almost triangular if B i,j = 0 for any i O(j), with a uniform constant.

In order to analyze our strategy we also require the matrices to be well conditioned in the following sense: Definition 5 (Knapsack-like matrix). Let B ∈ Z d×d be an almost triangular matrix and set C d 2 such that λ k (C) 2 C/k , for all matrices C whose columns are a subset of those of B of dimension k. Set R to be the R-factor of the QR-decomposition of B. We say that B is C-knapsack-like if furthermore R -1 2 C/d and |R i,j | 2 C/i for all i, j.

Remark 4. The conditions detailed in the previous definition seems apparently strong but such matrices are actually widespread, as corresponding to generic instances of so-called knapsack problems or searching integer relations. In practice, one can easily computationally check that these matrices, as well as Hermite Normal Form matrices with decreasing round pivots verify the assumptions with a reasonably small B.

Iterative reduction strategy

Hypothesis 1 In all of the following suppose that we have access to a lattice reduction oracle red-Oracle, whose output is a transition matrix to a Siegelreduced and size-reduced basis. Its running time on a d × d matrix of condition number bounded by C is denoted by T (d, C).

The progressive reduction consists in reducing the first k = 2 i columns of B, for all successive powers of two until reaching d. At step 1 i log d , we use the-now reduced-first k vectors to size-reduce the remaining columns before concatenating them to the current basis and pursuing the reduction. Hence, the bitsize of the whole matrix is reduced for each i before being actively used in the lattice reduction oracle red-Oracle.

Formally, define inductively a family of matrices B i which represents the state of the matrix B computed in the i-th iteration.

Initialization: B 0 = B. Induction: Let i > 0, and suppose that B i is known. We start by reducing only the first k = 2 i vectors using red-Oracle of B i and denote by B i the result. Define Q i R i to be the QR-decomposition of B i [: 1, k]. Then, remark that for any x being a column of B i not in the span of

B i [: 1, k],
we can reduce its bitsize by replacing it by x -B j R -1 j Q T j x for increasing

1 k, i ← 1, 0 2 B 0 ← B 3 while k < d do 4 k ← 2k; i ← i + 1 5 B i [: 1, k] ← B i-1 [: 1, k]•red-Oracle(B i-1 [: 1, k]) 6 R i ← Block-Cholesky(B i [: 1, k] T B[: 1, k]) 7 Q i ← B i [: 1, k] • R -1 i 8 B i+1 [1 : k, k + 1 : d] ← B i [1 : k, k + 1 : d]- 9 B i • R -1 i • Q T i • B i [1 : k, k + 1 : d] 10 end while

Complexity analysis

We now present the complexity analysis of the algorithm presented, under the hypothesis made on the lattice reduction oracle. For readability, we defer the proof to the full version. The following lemma entails that the condition number of the input of the reduction oracle is sufficiently small.

Q i R i = B i [, 1 : 2 i ] the QR-decomposition of the matrix of 2 i first columns of B i . We get R i , R -1 i = 2 O(d+C2 -i ) for all i.
From this we have:

Theorem 3. Let B a rank d almost triangular matrix which is C-well condi- tioned, C d 2 . We can Siegel-reduce it in time O log d i=1 T (2 i , O C2 -i ) + d ω-1 ω -2 • C log C + dC log d .
Remark 5. • One can use such a procedure to quickly search a putative minimal polynomial; the knapsack-like condition is however not guaranteed. • The setting of theorem 3 includes both modular and integer knapsacks.

• Assuming algorithm of section 3 has heuristically the right properties (which is the case in all of our extensive experiments), the complexity of the reduction of knapsack like matrices then becomes:

O d ω-1 (ω -2) 2 • C log C + dC log C .

Applications

Lattice reduction algorithms have numerous applications in mathematics and computer science. We survey here the impact of the implementation our algorithm, starting with cryptanalysis. In particular, we can reduce lattices of dimension in the thousands and with millions of bits. We recall that the Gram-Schmidt norms in the output basis are expected to decrease geometrically with rate 2 2α so that the Hermite factor in dimension d is 2 αd . For all the presented experiments, we use an Intel CPU E5-2695 v4 with 18 cores running at 2.10GHz processors; and 768 GiB of RAM. Some SSD swap was slightly used in the largest computation. For comparison with older timings, we used a machine with an Intel i7-8650U with 4 cores at 1.9 GHz. The program was compiled with Intel's libraries and compiler with the standard -Ofast low-level optimization flag.

Comparison with state of the art

We start this section by a comparison with the state-of-the-art implementation of fplll. Its complexity is O d4 B 2 in the general case, and its heuristic complexity 4 is O d 2 B 2 for knapsack matrices, as reported in [68, 1.5.3]. When d 220, its practical efficiency drops sharply due to the need of multiprecision computations. The following table presents a running time comparison with fplll, in single-threaded mode, on classical types of lattices namely knapsack and NTRU matrices. On the all instances, our implementation is sensibly faster than fplll. x 2 i ≈ 2 γ . Hence, lattice reduction with root Hermite factor 2 α can be used to construct a non-zero vector y ∈ Z d such that y , | x, y | 2 γ/d+αd . Indeed, any vector in this lattice is of the form ( x, y , y), so that its squared norm is the sum of two contributions: y 2 + | x, y | 2 . It now suffices to remark that norm of a vector found by reduction is smaller than the normalized covolume 2 γ/d times the root Hermite factor 2 αd .

Type

By plugging back the definition of the (x i ), we have x, y = s a, y + e, y where a = (a 1 , . . . , a d ), e = (e 1 , . . . , e d ). Assuming 2 γ/d+αd 2 η-ρ / √ d, the Cauchy-Schwarz inequality implies that a, y = 0. This is enough to break the scheme; if the (d -1) first vectors of the basis have this length, then the last one must be proportional to a (and is ±a if the entries are coprime). The first d -1 first vectors are orthogonal to a and are independent, so the last one must be proportional to a since a is in the lattice orthogonal to these vectors. The optimal d -for maximizing α -is therefore close to γ/α, leading to the condition α (η-ρ) 2 4γ . A part of the original paper [START_REF] Van Dijk | Fully homomorphic encryption over the integers[END_REF] considers security against polynomial-time adversaries, so that they obtain the condition γ = ω(η 2 log λ) for a "security parameter" of λ. Another part of the original paper [70, Section 6.3], and almost all follow-ups [START_REF] Coron | Fully homomorphic encryption over the integers with shorter public keys[END_REF][START_REF] Coron | Public key compression and modulus switching for fully homomorphic encryption over the integers[END_REF][START_REF] Coron | Scale-invariant fully homomorphic encryption over the integers[END_REF][START_REF] Cheon | Batch fully homomorphic encryption over the integers[END_REF][START_REF] Coron | Practical multilinear maps over the integers[END_REF][START_REF] Coron | New multilinear maps over the integers[END_REF], consider however security against adversaries able to do 2 λ operations. However, the condition was copied without change, which possibly explains why a large α was chosen in several implementations 5 .

As the lattice reduction algorithm can easily reach α = 0.04, this means we can use a smaller d, close to γ η-ρ for many instances than the d in the table, which is the dimension where α is maximal. In the instance where γ = 1.02 • 10 6 and η -ρ = 376, we used d = 3600 so that α = 0.024 was needed, and we used a pressed-bkz-19 in the leaves. This choice was made due to memory concerns.

While the large problems are clearly quite difficult, even the largest instances of the table seem to be within range of (motivated!) academic attackers, with terabytes of SSD memory and perhaps around 2 65 flop.

Overstretched NTRU

It is well-known since the work of Albrecht et al. [START_REF] Albrecht | A subfield lattice attack on overstretched NTRU assumptions -cryptanalysis of some FHE and graded encoding schemes[END_REF], Cheon et al. [START_REF] Cheon | Cryptanalysis of middle lattice on the overstretched NTRU problem for general modulus polynomial[END_REF] and Kirchner and Fouque [START_REF] Kirchner | Revisiting lattice attacks on overstretched NTRU parameters[END_REF] that an NTRU scheme with a very large modulus q compared to the dimension of the lattice is prone to attacks. However, these cases often happen in NTRU-based Homomorphic encryption schemes such as of the secret key sublattice is of size about the product of the d smaller Gram-Schmidt, it is possible to recover the secret key.

The optimal d is around log(q) 4α , which corresponds to a volume close to 2 log(q) 2 /16α . The scheme of [START_REF] Genise | Homomorphic encryption for finite automata[END_REF] chooses entries in F, G as integer Gaussians of standard deviation σ = r/π where r is the dimension of their lattice. We can restrict the lattice reduction to the middle 2d square matrix, but the volume is conserved. A more precise estimate consists in using the volume of the sublattice, projected orthogonally to the first r vectors of the reduced basis [START_REF] Kirchner | Revisiting lattice attacks on overstretched NTRU parameters[END_REF]. We expect the i-th Gram-Schmidt norm of the projected basis to be around √ r + 1 -iσ, so that the volume can be computed with Stirling's formula r i=1 (r + 1 -i)σ to a SVP in dimension 101, and each call currently needs 2 11 core-seconds [START_REF] Albrecht | The General Sieve Kernel and New Records in Lattice Reduction[END_REF], this translates into a year of computation on our machine. Alternatively, each call can be computed in 2 2 seconds with a GPU [START_REF] Ducas | Advanced lattice sieving on GPUs, with tensor cores[END_REF]. A pressed-bkz of dimension 29 was used for the second one.

Miscellaneous

Integral relations Another use of lattice reduction is the discovery small linear integer relation between reals. It actually corresponds to the setting of section 6.2, where 2 η corresponds the norm of the relation, and γ the precision used to represent the reals. Then clearly, γ ≈ dη + d 2 α is enough to perform a search by reduction. In 2001, Bailey and Broadhurst believed [START_REF] Bailey | Parallel integer relation detection: Techniques and applications[END_REF] that their computation with γ ≈ 166000 and d = 110 was the largest performed. It took 44 hours, on 32 CPUs of a Cray-T3E (300 MHz). We report this takes 5 minutes on a laptop, or 600 times fewer cycles. As the task is identical (for large α) to breaking the integer homomorphic schemes, the running time for bigger examples can be found in the previous subsections.

Univariate polynomial factorization Yet another application is factoring univariate polynomials [START_REF] Van Hoeij | Factoring polynomials and the knapsack problem[END_REF][START_REF] Belabas | Factoring polynomials over global fields[END_REF] over the integers. The first step is to factor modulo some prime, and the number of factors n is the dimension of the modular vectorial knapsack we have to solve, namely we have to find very short vectors in the lattice generated by qId r A 0 Id n . The precision q, and number of coordinates r can essentially be freely chosen. For random polynomials, n is typically very small (e.g. logarithmic) and lattice reduction is not the bottleneck; but it can be as large as half the degree. Our choice is to take r small, say n/ log n, and then r log q ≈ αn 2 allows (heuristically) to obtain the last Gram-Schmidt norms larger than n 2 . Then, this restricts the solutions of the knapsack -known to be shorter than this -to the first few vectors of the reduced basis. At this point, one can recover the factors, and prove that they are irreducible. Taking n = 256, we get a solution in two minutes on one core of our laptop instead of ten with a 1 GHz Athlon [START_REF] Belabas | A relative van Hoeij algorithm over number fields[END_REF]; for n = 512 it takes 25 minutes instead of 500. For ω bounded away from 2, with α = O log log n log n the heuristic asymptotical complexity is O n ω+1 log log n log 2 n .

Conclusion and open questions

In this work, we introduced a recursive lattice reduction algorithm, whose heuristic complexity is equivalent to a few matrix multiplications. This algorithm and the heuristics used to complete the complexity analysis have been thoroughly tested and applied to reduce lattices of very large dimension. The implementation takes advantage of fast matrix multiplications, and fast Fourier transforms. This work raises several questions. First of all, the analysis we are making is so far heuristic and empirical. It is possible to get a provable result by mitigating the complexity, in particular it seems difficult to be able to formally prove the heuristic on the decrease of the needed precision, even though this fact is easily checkable in practice. Reaching a provable bound in d ω C is an open and interesting problem, and our algorithm is a first step in this direction. Remark 6. It is mandatory to have T well-conditioned if we want a U which is not much larger than T. This is also true for other variants of lll (including fplll): outputting the transition matrix may lead to a slow-down by a factor of d.

  for i = 1 to d do for j = i -1 to 1 do

  of a family of vectors. Let Λ be a rank d lattice given by a basis B = (b 1 , . . . , b d ), we might want to use the Gram-Schmidt process. However, since the quotients bi,b * j b * j ,b * j are not integral in general, the vectors b * i may not lie in Λ. The size-reduction process instead approximates the result of the Gram-Schmidt process by rounding to a nearest integer: each vector b i is replaced by b i -The whole process takes time O d 5 B 2 when

Theorem 1 .

 1 Given A and B two integer matrices of dimension d with log( A + B ) = p, the product AB can be computed in time O d ω p+w w + d 2 p w log 2 + p w .

  M // with dimension almost halved 5 A ← Invert(A) ; D ← Invert(D) // fp-approximations of A -1 , D -1 6 S ← -A CD // Computed in floating-point 7 return A S 0 D

Lemma 1 .

 1 Given an integral unitriangular matrix M of dimension d, with both M , M -1 2 p and p w + log(d), Invert returns a matrix M such that M -M -1 2 -p with a running time of O d ω p w + d 2 p . Proof. We set a working precision p = 1 + 3p + log d = O(p), and by induction on d, let us prove that

Algorithm 5 -

 5 Reduce Parameter : Relaxation factor α, ε > 0, number of rounds ρ, number of blocks D, d = d/D block size. Input : Basis B ∈ Z d×d of the lattice Λ Output : A unimodular transformation U ∈ Z d×d , UB reduced.

Theorem 2 .

 2 Given a d-dimensional unitriangular matrix T such that T and T -1 2 p and p w+log(d) 2 d, the algorithm Size-Reduce returns an integral unitriangular matrix U with U 2 O(p) such that TU , (TU) -1 d log d with a running time of O d ω p w + d 2 p .

Fig. 2 :

 2 Fig. 2: The block process on the Rpart of the basis. Green blocks are recursively reduced on odd steps and oranges one are reduced on even steps.

2 .

 2 2 log C bound on the number of iterations ρ is now adequate. Let us denote by R j the extracted submatrix (R a,b ) (j-1)d <a,b jd , with d = d/D. The lattice R j spanned by R j is the projection of Λ j = (j-1)d <a jd b a Z over the orthogonal space to the first (j -1)d vectors (b 1 , . . . , b (j-1)d -1 ). The step reduction subprocess simultaneously (and recursively) calls the reduction of all the shifted sublattices R 2j ⊕ R 2j+1 1 j< D Then the same is done on the sublattices R 2j+1 ⊕ R 2j+2 0 j< D 2

1 Fig. 3 :

 13 Fig. 3: Abscissa corresponds to the iteration time and ordinates corresponds to the value log(max i R[i, i]/ min i R[i, i]).As predicted by heuristic 1 the corresponding graph presents an exponential decay.

Fig. 4 :

 4 Fig. 4: Log-Log representation of the running time (in seconds) for increasing dimension, with constant on each line.

  Such a size-reduction can be computed on all the columns of B i [k + 1 : d] simultaneously using a single matrix multiplication and call C the corresponding vectors. Eventually set B i+1 to the concatenation B [: 1, k] | C . The corresponding pseudo-code is given in algorithm 8. Algorithm 8 -Reduction of Knapsack-like lattices Parameter : Reduction oracle red-Oracle Input : Matrix B ∈ Z d×d Output : A reduced basis of BZ

Lemma 2 .

 2 Let B a rank d almost triangular matrix which is C-knapsack-like. For any index 0 i log d , set B i to be the matrix computed by the execution of algorithm 8 on B. Denote by

  lwe) problem in dimension one. More precisely, given an integer secret |s| 2 η (typically a prime), this problem aims at retrieving s from given samples x i of the form a i s + e i where 0 a i 2 γ /|s| and |e i | 2 ρ are sampled uniformly and independently. The parameters verify γ η ρ. A natural lattice reduction attack consists in collecting d samples x = (x i ) 1 i d and building the matrix X = x 1 , . . . , x d Id d . The volume of the lattice X spanned by the columns of X is 1 +

	With Error (d
		i=1
		Dimension d	Bitsize B	fplll	This work
		128	100 000	88 min	6 min
	Knapsack	256	10 000	134 min	4 min
		384	10 000	388 min	13 min
		256	80	24 min	3 min
	NTRU	384	70	431 min	10 min
		512	70	1392 min	33 min
	6.2 Fully Homomorphic Encryption over the integers
	FHE scheme was first designed by Gentry [33] using number theoretical tools in
	2009. Soon after, an equivalent system was presented, using only integer arith-
	metic [70], and is based on a distant relative [17] of the celebrated Learning

Table 2 :

 2 can find the α required. The first one necessitates roughly 2 20 calls Experiments for overstretched NTRU problems. Dimension is the actual dimension of the problem, and effective dimension refers to the dimension required in practice to mount the attack.

		2	1/2	≈	√ √ rσ e	r	. Overall we obtain 2 log(q) 2 /16α ≈	r √ πe	r	,
	from which we Dimension r log q		α		Effective dimension Algorithm		Time
	1024	42	0.01274		1648	BKZ-101		-
	4096	111	0.01799		3086	BKZ-25	233 hours
	32768	883	0.1105		2560	LLL	263 min

We choose to present it using the matrix R and yielding the unimodular transformation matrix, for consistency with the description of our algorithms in section 3.

This estimation comes from the observation that dB swaps are performed, each taking d 2 operations on B/d bits. There are B/d reduction steps for each new vector, each takes d 2 operations on B bits.

Surprisingly, many log λ were "rounded up" to λ in the parameter choices.

Acknowledgement

This work was supported by the European Union PROMETHEUS project (Horizon 2020 Research and Innovation Program, grant 780701). Paul Kirchner is partly funded by the Direction Générale de l'Armement (Pôle de Recherche CYBER).

YASHE [START_REF] Bos | Improved Security for a Ring-Based Fully Homomorphic Encryption Scheme[END_REF] or LTV schemes. In 2019, a homomorphic scheme has been proposed by Genise et al. [START_REF] Genise | Homomorphic encryption for finite automata[END_REF] with a similar variant of this problem, hoping that the overstretched NTRU only works in algebraic setting using ring of polynomials. Some parameters proposed for performances evaluations have been broken in [START_REF] Lee | Lattice analysis on MiNTRU problem[END_REF] also showing that the assumptions used is flawed. Here, we break comparable parameters showing that the proposed parameters only achieve a low security level. In [START_REF] Pataki | Lattice determinants in reduced bases[END_REF], Pataki and Tural showed that the volume of any r-dimensional sublattice L of a lattice L is larger than the product of the r smaller Gram-Schmidt. Kirchner and Fouque combined this result with the fact that in any 2d-dimensional NTRU lattices, there is a sublattice of dimension d and volume roughly the size of secret-key to the power d, one can deduce that if the volume

A Proof of Theorem 2

Proof. We use a precision p = O p + log(d) 2 = O(p). We prove by induction on d that TU , (TU) -1 d log d . Initialization is clear, so that we now assume that d > 1. We have by direct computation that TU is

The top-right matrix is AU 1 ((AU 1 ) -1 CU 2 -W) and we have, by setting that A -(AU 1 ) -1 = δA:

The first term is bounded by 2 O(p) δA and the second by 2d/3. We choose the precision so that the first term is at most 1/3 and the result follows directly, as

The top-right matrix is ((AU


Bruno Denis

Jean-Jacques Lesage

Zulema Juarez Orozco

Jean-Jaques Lesage

Zulema Juárez-Orozco

Performance verification of discrete event systems using hybrid model-checking

Keywords: DES controller, hybrid plant model, model-based verification, modelchecking, linear hybrid automaton, HYTECH

come   L'archive ouverte pluridisciplinaire

INTRODUCTION

A physical system is not, in most cases, intrinsically either purely discrete or purely continuous; instead, it's the abstraction the control engineer undertakes to ensure automation specification is being met that lends one distinction or the other. In the area of manufacturing systems for example, many processes have mobility axes that enable products to circulate, establish their position, etc.; consequently, some of the major physical variables controlled are either displacements or speeds. Depending not only on the level of quality to be guaranteed for these controlled variables, but also on the aggressiveness or variability in the external environment (e.g. type and importance of perturbations, parameter variation interval for the law of movement) and on the relative weight of economic constraints, the control engineer is required to choose between a servo control or a discrete control for each of these axes. Such automation-related choices will yield the physical variables to be observed and controlled by the controller, with either continuous or discrete control abstraction. Once the requisite displacement axis positioning quality has been achieved and provided that the perturbations encountered remain tolerable, the control engineer will then select a discrete control for obvious cost reduction reasons. This discrete displacement control, despite often being able to accommodate mobile positioning quality requirements, still constitutes an abstraction, and as such necessarily a simplification of the associated physical variables. The logic control of a linear displacement between two extreme positions, observed by means of two limit switch sensors, clearly does not enable ascertaining the precise position of the mobile, nor the time elapsed to complete this displacement. Satisfying industrial system dependability requirements often necessitates conducting offline analyses, such as formal verification, before placing the automated system into operation (for further information on this topic, see the standard IEC 61508 entitled " Functional safety of electrical / electronic / programmable electronic safety-related systems "). These verifications, which are now frequently performed by means of model-checking (MacMillan, 1993), may be practiced by electing to incorporate or not a plant model.

In this paper, our efforts have focused on checking systems composed of a discrete controller coupled with a continuous (or partially-continuous) physical process whose entire set of observed and controlled variables constitute discrete abstractions of physical variables. To proceed, we will make use of modelchecking techniques by coupling the discrete controller model with a hybrid plant model; this setup will demonstrate that beyond the liveness and safety properties, it is indeed possible to check whether automated system performance is compatible with that stipulated in the specifications. We will also show that by parameterizing the hybrid process model, it becomes possible to use the modelchecker to evaluate the robustness of discrete control to perturbations encountered by the plant. This paper has been organized as follows. After having recalled the possibilities and limitations of both DES and HDS model-checking, we will introduce the expectations derived from a hybrid plant model for verifying a discrete controller. In order to illustrate our approach, the paper's second part will present the example of a positioning axis, which represents a component of a more complex assembly system. We will thus be able to show that use of a model-checker (such as HYTECH), by implementation of a hybrid process model, makes it possible to verify the expected performance of this DES. Furthermore, a sensitivity study conducted on model parameters will allow evaluating the robustness of discrete control when confronted with perturbations.

ACQUIRED KNOWLEGE AND LIMITATIONS FROM AUTOMATED SYSTEM VERIFICATION

DES verification.

Formal verification techniques stem from the field of computer science. Only recently have they been adapted and applied to DES verification and, more specifically, to model-checking [START_REF] Clarke | Automatic verification of finite state concurrent system using temporal logic[END_REF]. The general principle behind model-checking may be expressed as follows (see Figure 1). Let's start with a system that has been designed to verify an entire array of properties (logical correctness, dependability, liveness, etc.). The first task of model-checking consists of formalizing system behavior in the form of a finite state automaton: S, plus the properties to be verified within a temporal algebra such as CTL [START_REF] Emerson | Sometimes and Not Never revisited : on branching versus linear time temporal logic[END_REF]: ϕ. The model-checker then conducts a thorough analysis of the state space reachable by S, which serves either to prove that S ╞═ ϕ (this algebraic statement denotes that "the system model satisfies the set of properties ϕ") or, when such is not the case, to propose a counterexample that revokes those properties not verified by S. Moreover, a DES may be represented in a generic manner, as shown in Figure 2: a discrete controller acting in a closed loop on a plant. As part of a dependable controller design approach, the system being targeted for verification can thus be (according to [START_REF] Frey | Formal method in PLC programming[END_REF] either the controller on its own, presumed to be operating within an open loop on the plant (a non "model-based" verification), or the {controller + plant} assembly set interacting within a closed loop ("model-based" verification). The research work focusing on DES verification initially favored a non model-based approach [START_REF] Moon | Modeling programmable logic controllers for logic verification[END_REF]. The reachable state space of the controller model is thus to be built in the most permissive manner possible, i.e. such that the evolution of its inputs are in no way constrained by plant behavior. In this case, the safety properties capable of being demonstrated provide the basis for strong proof given that they can be demonstrated regardless of the evolution in controller inputs. On the other hand, a good number of liveness or accessibility properties cannot be demonstrated via a non model-based model-checking approach given the often fast-paced combinatory explosion of the reachable state space. One means for reducing this combinatory explosion, using realistic constraints that depict the interaction of plant behavior with controller behavior, is to conduct a model-based verification. [START_REF] Rausch | Formal verification of PLC programs[END_REF], (Machado, et al., 2003). Through reliance upon these results, we are now in a position to study the possibilities offered by the model-checking procedure in evaluating DES performance (in terms of quality of service provided by the automated system). To accomplish this task, it is necessary to undertake a model-based approach by selecting a hybrid plant model. The physical variables of the plant, and not their discrete abstraction by the controller, are what give the actual performance measures to be evaluated. The coupling of a hybrid plant model with a discrete controller model thus leads to a Hybrid Dynamic System (HDS) verification problem. We will now proceed by recalling the basic knowledge acquired and current limitations of HDS verification.

Hybrid systems verification.

While model design using hybrid automata is not exactly straightforward, their semantics is welladapted to the analysis of HDS behavior. We will then assume that HDS verification is tantamount to exploring the reachable state space of a hybrid automaton. When framed as such, the fundamental problem of HDS verification becomes one of computing the reachable state space of a hybrid automaton from an initial region. Not only is this computation prone to encountering a rapid combinatory explosion, but in the general case the problem is not-decidable, i.e. impossible to identify a computation algorithm with guaranteed convergence (Henzinger, et al., 1995a). The linear hybrid automata [START_REF] Alur | The algorithmic analysis of hybrid systems[END_REF]) constitute a category of hybrid automata that places verification within the domain of a decidable problem and for this reason, such automata have been selected for the purposes of our work. From an other hand, modular approaches are essential for modeling and verifying of industrial systems. This is why we felt that the HYTECH tool (Henzinger, et al., 1995b) was adapted to our needs, by virtue of allowing for the verification of models obtained by composition of linear hybrid automata. With the context and objectives of our research now established and both the category of hybrid automata and model-checker selected, we will turn our attention next to presenting our approach for verifying discrete control performance via the assessment of a pertinent example.

A CASE STUDY

3.1 Presentation of the production system.

The system we are targeting is one of assembling / disassembling bearings within gears; it comprises four workstations (see Figure 3a). The four workstations pertain respectively to the following functions: loading at the beginning of the line (Station 1); identification of the gear material and eventual presence of a bearing (Station 2); insertion or removal of a bearing in the gear (Station 3); and gear sorting, taking into account the material, in anticipation of their unloading (Station 4). This entire assembly line is controlled by a Programmable Logic Controller (PLC). Control of this assembly system is purely discrete; since this proves to be the most typical case for such systems, seeking more competitive costs leads to avoiding servo controls as much as possible in favor of logic controls. All continuous physical variables of the process being controlled (which for the most part consist of positions) are thus observed and controlled via their discrete abstractions. In all, 82 logical inputs and 50 logical outputs are managed by the PLC. The control model, written using Sequential Function Chart (SFC) syntax according to IEC 61131-3 standard has been verified with a timed model-checking technique [START_REF] Mokadem | Verification of a timed multitask system with UPPAAL[END_REF]. The objective of our study herein consists of verifying whether the quality of the positions obtained using such a discrete control is compatible with the precision required by the process. This concept of precision only takes on meaning with respect to the positioning of objects manipulated by either a displacement axis or a group of combined axes. It is thereby unnecessary to include the entire control model and the plant model in order to verify position quality, but rather just the subset that pertains to control of the axis (or group of axes) involved in the studied displacement. Our investigation will focus more heavily on the displacement axis of station 2 (Fig. 3b), with the only relevant information available for displacement control being given in Figure 4. The gear is initially located on the conveyor at the loading station (Load_pos). A motor rotation order in a clockwise direction (Go_to_test) implies the conveyor in displacement to the test position (Test_pos), a detailed view is provided in Figure 3c.

For the test to be successfully conducted, stopping precision must be ≤ 1.5 mm. The conveyor position between these two discrete positions is denoted X.

When the presence test of a bearing is completed, a subsequent motor rotation order in the clockwise direction positions the gear, where it gets unloaded (Unload_pos). A motor rotation in the counterclockwise direction then brings the conveyor into the loading position, where a new cycle can be launched. We are seeking to quantify the precision of gear positioning at the test station, enabled by this discrete control along the path Load_pos → Test_pos, so as to verify whether this precision is compatible with specification. Execution of this task will entail application of hybrid model-checking using the HYTECH tool (Henzinger, et al., 1995b).

Modelling.

The model of the entire {controller + displacement axis} set, which allows this positioning quality verification, has been laid out in Figure 5. The hybrid automaton is composed of two synchronized automata.

Q0 Q1 0 Tc ≥ 1 c T - = & ] Tcmax [Tcmin, : Tc Xtest) (X 0) (Tc = → < ∧ = Q1a Q1b Tcmax] [Tcmin, : Tc = → 0 Tc ≥ 1 c T - = & Tcmax] [Tcmin, : Tc = → 0 Tc ≥ 0 Tc = 1 c T - = & ] [ Vmax Vmin, X = & start Load_pos 0 X = & 0 X = 0 X = & 0 Do ≥ 1 o D - = & stop 0 Do = 0 X ≥ Domax] [Domin, : Do = → 0 : X = → 0 : X = → Go_to_ test

Test_pos

Displacement axis

Fig. 5. Hybrid automata of the axis control.

The automaton that model electromechanical axis contains three discrete locations: two correspond with conveyor stopping positions at the extremities (Load_pos and Test_pos), and the other location corresponds with the movement phase between these extreme positions (Go_to_test). The return to loading position, which does not enter into the scope of the present study, is modeled by means of a timer with variable duration D 0 . The two continuous state variables associated with this automaton are conveyor position X and timer duration D 0 . Upon initialization, the conveyor is placed into the loading position and hence X = 0. Both of the stopping situations are associated with the function dX/dt = 0 (denoted X & = 0); the movement location is associated with a continuous dynamics of the form X & = a, with a representing a rational constant arbitrarily chosen within an interval: [V min , V max ]. This arbitrary assignment of constant a enables abstracting the variability in conveyor displacement speed due to perturbations undergone by the displacement axis (mechanical friction, resistant torque on the motor axis, kinematic nonlinearities, etc.). Two synchronization events with the automaton used for modeling controller behavior, called start and stop respectively, are associated with the transitions that correspond to the conveyor movement phase. The automaton that models controller behavior is composed of four discrete locations and is initialized at location Q 0 . The start synchronization event with the axis automaton shows that the motor starts up during a PLC cycle change. The monitor of this PLC is of the "periodic" type, indicating that execution of the control program is to take place according to a constant-duration cycle (see Figure 6), which the user is required to setup in PLC. In practice, this cycle duration undergoes slight variations (on the order of 1%), called "jitters". To translate these operations into the hybrid controller automaton, the duration of each PLC processing cycle (Tc) is assigned an arbitrarily-chosen value within an interval: [Tc min , Tc max ]. Once location Q 1 has been entered, the automaton remains there until the axis has traveled the distance required to position it in alignment with the test station (at which point X = Xtest). Locations Q1 a and Q1 b serve to translate that this information is taken into account by the PLC during the subsequent input reading phase, followed by processing and then execution of the motor stop order (Figure 6). The stop synchronization event thereby triggers controller automaton re-initialization at the same time it activates the Test_pos situation of the displacement axis automaton. Figure 6 shows the synchronization of the two automata presented in Figure 5 by the stop event as well as the evolution in continuous variables Tc and X. The axis positioning error at the level of the test station thus results from the time lag between the moment the axis reaches the Xtest position and the moment the PLC is able to react at this event in the input reading phase, with the consequence being that the axis stop order may be issued either within the same PLC cycle or one cycle later. The experimental campaign conducted using HYTECH has revealed this finding and the following discussion is intended to provide greater detail.

Implementation and the experimental protocol.

All of our experiments have been carried out by constraining the variation of duration the conveyor's return to the loading position (D 0 ) over the interval [3000ms, 3500ms] (which thus yields Domin = 3000ms and Domax = 3500ms). Moreover, the abscissa (Xtest) of the conveyor at the test station is 80 mm. Given that these data remain fixed, the HYTECH code obtained from the hybrid automaton shown in Figure 5 is a parametric code, with intervals [Tcmin, Tcmax] (for the cycle time to be adjusted on the PLC) and [Vmin, Vmax] (for the translation speed of the conveyor subjected to perturbations) constituting the two parameters. For each experiment, numerical values are ascribed to these two intervals. The model-checker then computes the region the automaton is capable of reaching from its initial region. We will initially derive the intersection of this reachable region with the discrete location Test_pos so as to build the accessibility domain of the axis position, which is then obtained by projection of this sub-region with respect to the variable X. The set of positions reached by the conveyor thus assumes the form of an interval on X.

As an example, for Vmin = Vmax = 200mm/s, Tcmin = 11ms -0.5% and Tcmax = 11ms +0.5% , the conveyor stopping position in front of the test station is contained within the interval [83.182mm, 84.018mm], which represents a resultant position deviation of 0.836mm.

Position-stopping quality without perturbations.

For this initial study, the hypothesis has been adopted that the conveyor is not submitted to any perturbation. Its translation speed is thus presumed to be constant and yields: Vmin = Vmax = 200 mm/s. The changes in conveyor stopping position may then be observed as a function of the PLC cycle time value when such value varies by more or less than 0.5% around the period adjusted by the user (the socalled jitter phenomenon described in Section 3.2). Figure 7 reveals that the positioning error evolves in accordance with two phenomena that we will now differentiate. The first corresponds to a systematic error (positive slope line passing through the average error interval value), which can be corrected by an adapted adjustment of sensor position Test_pos. The second phenomenon is homogeneous to a random error (the fixed interval deviation on the positions reached for a given Tc), whose general amplitude shape increases with PLC cycle time. It can nonetheless be remarked that for certain Tc values (e.g. Tc = 9ms or Tc = 13.1ms), range of these "random" errors is narrower than that for proximate Tc values. The existence of these particular zones of smaller positioning error, while PLC cycle time is increasing, highlights a non-trivial temporal correlation phenomenon between the PLC cycle and occurrence of the test position arrival event depicted in Figure 8. When the temporal detection region of the X=Xtest event has been included within a single PLC cycle, conveyor stopping will systematically be controlled during the subsequent cycle output assignment phase (see Figure 8a). In contrast, while the temporal detection region of the X=Xtest event straddles two PLC cycles, stopping will intervene during the output assignment phase of one of the two subsequent cycles (Figure 8b). The range of the stopping zone will thus increase. This temporal dispersion in stop order execution is heavily correlated with the PLC cycle time jitter (Figure 8c). Though the ratio between the distance to be traversed and Tc is such that, despite the jitter on Tc, the number N of PLC cycles elapsed between the conveyor start order and detection of the test position (X=Xtest) is the same with each new conveyor displacement from X=0 to X=Xtest, the conveyor stop order will always be issued at the end of the same PLC cycle and the error range will always solely depend on the jitter and cycle number N (this configuration will be referred to as "Scenario 1" in the following discussion). On the other hand, should the ratio between the distance traveled and Tc be such that a variable number of cycles is necessary for conveyor displacement, the stop order may then be issued at the end of the various PLC cycles (Figure 8c). The error range is thus merely a function of Tc (this configuration will be called "Scenario 2"). Figure 9 clearly highlights this phenomenon. For Scenario 2, the stopping range is proportional to the cycle time value. Whereas for Scenario 1 positioning error is independent of cycle time (at least up to the limit of values selected for our experiment).

If specifications call for a maximum positioning error in Test_pos of 1.5mm, the PLC cycle period has to be set at a value of below 7.4ms. If the full PLC load (for execution of the overall assembly line program with its 82 inputs and 50 outputs) does not enable parameterizing such a low cycle time, it would then be necessary to target a Scenario 1 parameterization, e.g. 9.0ms or 13.1ms.

Positioning robustness.

It has already been demonstrated that the PLC cycle time may be judiciously chosen in order to guarantee the quality of conveyor positioning under a hypothesis of no translation speed perturbations. Figure 10 presents the results from experiments conducted in order to evaluate the influence of such perturbations on positioning quality for a given cycle time Tc=13.1ms (which remains contaminated with a jitter of ± 0.5%). This curve shape is the same for Tc = 9ms or Tc = 5.1ms. Two distinct zones appear on this curve. Above 2% speed perturbation, the positioning error is basically constant; this corresponds to Scenario 2, for which the error range is directly correlated with Tc without any major influence on the speed-related perturbation. Below this 2% level, the error increases linearly with perturbations since the system here is behaving according to Scenario 1, which is sensitive to both the jitter on Tc and perturbations on the conveyor translation speed. For a requirement that limits the positioning error to 1.5mm, it thus becomes necessary to ensure that speed perturbations remain less than 0.8% (in the case of a PLC adjustment to Tc=13.1ms). Should such not be the case, adopting a servo control in the conveyer position would need to be examined.

Assessment of the use of HYTECH.

The results provided above have necessitated several hundreds of model-checking experiments on the hybrid automaton shown in Figure 5 using the HYTECH tool, with different parameter sets being introduced for each experiment. The search for the reachable region of a model never took more than 5 minutes in most instances, 30 seconds were sufficient. Computation time did not therefore pose any problem for this study. The primary difficulty encountered had to do with the interruption of reachability computations due to exceeding the integer limit value (2 64 ) during the region computation step. Our case study did necessitate manipulating numerical values with quite varied scales. As an example, the conveyor position progressed over a total path length of 80mm, whereas the calculated stopping error at times amounted to just a few hundredths of a millimeter. For all of the distances manipulated, numerical values out to 4 significant digits thus had to be used. The same difficulty gets magnified when manipulating time since the duration of the conveyor position return lasts approximately 3s, while a ± 0.5% jitter on Tc = 1ms equals ± 5µs; seven significant digits therefore proved necessary. This numerical value coding is very detrimental when running HYTECH, which interrupts the computation of regions once the integer coding capacity has been exeeded, thereby making it impossible to draw any conclusion on the current verification. In order to avoid this predicament, the units for each model parameter set have been adjusted. For example, depending on the value of Tc, the time unit has been fixed at 1µs, 10µs or 100µs, which makes only 5% of the verifications inconclusive.

CONCLUSION AND OUTLOOK

In this paper, we have focused on the verification of systems composed of a discrete controller coupled with a continuous (or partially-continuous) physical process, whose full extent of observed and controlled variables are logical abstractions of physical variables. For that, we employed a model-checking technique that couples the timed discrete model of the controller with a hybrid plant model. By relying upon a case study, we showed the possibility of verifying whether or not system performance is compatible with the specifications. By parameterizing the hybrid system model, it becomes possible for the model-checker to evaluate discrete control robustness to plant perturbations. According to this approach, HYTECH has proven to be an efficient tool. It goes without saying that the size of the model being verified must remain compatible with the difficulties inherent in the wellknown combinatory explosion phenomenon typical of model-checking. For this reason, it is necessary to limit performance verification to just the {control + plant} subsystem involved in deriving the evaluated variables. Determining this subsystem being targeted by performance verification along with building the hybrid model used the abstraction are the current focus of a complementary research project.

  Fig. 1. Model-checking scheme.

  Fig. 2. A generic closed-loop DES.

Fig. 3 .

 3 Fig. 3. Assembly/disassembly line.

Fig. 4 .

 4 Fig. 4. Displacement axis of station 2.

Fig. 6 .

 6 Fig. 6. Synchronous evolutions of automata.

Fig. 7 .

 7 Fig. 7. Conveyor positioning error.

  Fig. 8. Correlation between Tc, Xtest event occurrences and stop order.

Fig. 9 .

 9 Fig. 9. Conveyor positioning error range.

PhD grant financed by CONACYT (Mexico