Certified lattice reduction

TL;DR

This paper introduces an adaptive-precision LLL algorithm using Interval Arithmetic to certify lattice reduction of real quadratic forms, providing guarantees beyond traditional floating-point methods, with applications in algebraic number theory.

Contribution

It presents a novel adaptive-precision LLL algorithm that certifies lattice reduction of real quadratic forms using Interval Arithmetic, extending the reduction theory.

Findings

Certifies lattice reduction with guaranteed correctness.

Applicable to real quadratic forms and general lattices.

Demonstrates effectiveness in algebraic number theory applications.

Abstract

Quadratic form reduction and lattice reduction are fundamental tools in computational number theory and in computer science, especially in cryptography. The celebrated Lenstra-Lenstra-Lov\'asz reduction algorithm (so-called LLL) has been improved in many ways through the past decades and remains one of the central methods used for reducing integral lattice basis. In particular, its floating-point variants-where the rational arithmetic required by Gram-Schmidt orthogonalization is replaced by floating-point arithmetic-are now the fastest known. However, the systematic study of the reduction theory of real quadratic forms or, more generally, of real lattices is not widely represented in the literature. When the problem arises, the lattice is usually replaced by an integral approximation of (a multiple of) the original lattice, which is then reduced. While practically useful and proven in some special cases, this method doesn't offer any guarantee of success in general. In this work, we present an adaptive-precision version of a generalized LLL algorithm that covers this case in all generality. In particular, we replace floating-point arithmetic by Interval Arithmetic to certify the behavior of the algorithm. We conclude by giving a typical application of the result in algebraic number theory for the reduction of ideal lattices in number fields.

Equations

$$b_1\mathbf{R}\subset b_1\mathbf{R}\oplus b_2\mathbf{R}\subset \cdots \subset \underset{i=1}{\overset{r}{⨁}}{b}_i\mathbf{R}.$$

$${\pi }_1(b_1)$$

$$\forall 1<i\le r,{\pi }_i(b_i)$$

$$\operatorname{covol}\mathopen{}\mathclose{{}\left({S}}\right)=\prod_{i=1}^{r}\|\pi_{i}(b_{i})\|$$

$$\Vert x-y\Vert \ge {ϵ}_{\Lambda }.$$

$$\forall i<j,\quad\mathopen{}\mathclose{{}\left|\langle{{b_{j}}},{\pi_{i}(b_{i})}\rangle}\right|\leq\frac{1}{2}{\|\pi_{i}(b_{i})\|^{2}}\quad\textrm{(size-reduction condition)}$$

$$\forall i,\quad\delta\|\pi_{i}({b_{i}})\|^{2}\leq\mathopen{}\mathclose{{}\left(\|\pi_{i+1}(b_{i+1})\|^{2}+\frac{\langle{b_{i+1}},{\pi_{i}(b_{i})}\rangle}{\|\pi_{i}(b_{i})\|^{2}}}\right)\quad\textrm{(Lovász condition)}$$

$$\textrm{O}\mathopen{}\mathclose{{}\left(\frac{d^{5}\log^{2}\|B\|_{\textrm{max}}}{d+\log\|B\|_{\textrm{max}}}\mathcal{M}(d+\log\|B\|_{\textrm{max}})}\right).$$

$$\operatorname{covol}\mathopen{}\mathclose{{}\left({b_{1},\ldots,b_{k}}}\right)\leq\mathopen{}\mathclose{{}\left(\delta-\frac{1}{4}}\right)^{-\frac{(d-k)k}{4}}\operatorname{covol}\mathopen{}\mathclose{{}\left({\Lambda}}\right)^{\frac{k}{d}}.$$

$$\forall i<j,\quad\mathopen{}\mathclose{{}\left|\langle{{b_{j}}},{\pi_{i}(b_{i})}\rangle}\right|\leq\eta{\|\pi_{i}(b_{i})\|^{2}}\quad\textrm{(Approximate size-reduction condition)}$$

$$\textrm{O}\mathopen{}\mathclose{{}\left(d^{5}(d+\log(\|B\|_{\textrm{max}}))\log(\|B\|_{\textrm{max}})}\right).$$

$$G=M^T\cdot {B^{\ast }}^T\cdot B^{\ast }\cdot M=M^T\cdot D\cdot M$$

$${d\log\mathopen{}\mathclose{{}\left(\frac{(1+\eta)^{2}}{(\delta-\eta)^{2}}+\epsilon}\right)+o\mathopen{}\mathclose{{}\left(d}\right)}$$

$$\underline{x}_{n}=\mathopen{}\mathclose{{}\left[X_{n}-1,X_{n}+1}\right]$$

$$\underline{x}_{n}=\mathopen{}\mathclose{{}\left[X_{n}-\delta,X_{n}+\delta}\right]$$

$$\displaystyle\mathopen{}\mathclose{{}\left[\underline{x}^{-},\underline{x}^{+}}\right]+\mathopen{}\mathclose{{}\left[\underline{y}^{-},\underline{y}^{+}}\right]$$

$$\displaystyle\mathopen{}\mathclose{{}\left[\underline{x}^{-},\underline{x}^{+}}\right]-\mathopen{}\mathclose{{}\left[\underline{y}^{-},\underline{y}^{+}}\right]$$

$$\displaystyle\mathopen{}\mathclose{{}\left[\underline{x}^{-},\underline{x}^{+}}\right]\times\mathopen{}\mathclose{{}\left[\underline{y}^{-},\underline{y}^{+}}\right]$$

$$\displaystyle\mathopen{}\mathclose{{}\left[\underline{x}^{-},\underline{x}^{+}}\right]^{-1}$$

$${\underline{A}}_n=({\underline{a_{i,j}}}_n{)}_{(i,j)\in [1\cdots d{]}^2},$$

$$2^n{\lambda }_d(S)-2d\le {\lambda }_d(S^{\prime })\le 2^n{\lambda }_d(S)+2d,$$

$$⟨a,b⟩=A^T\cdot {\mathcal{G}}_{\gamma }\cdot B.$$

$$2^{n}\langle{a},{b}\rangle=A^{T}\cdot G\cdot B+A^{T}\cdot\Delta\cdot B,\leavevmode\nobreak\ \mbox{where}\leavevmode\nobreak\ |A^{T}\cdot\Delta\cdot B|\leq\mathopen{}\mathclose{{}\left(\sum_{i}|A_{i}|}\right)\mathopen{}\mathclose{{}\left(\sum_{i}|B_{i}|}\right).$$

$$\ell \ge 10+2{\log }_{2}d-{\log }_2\min (ϵ,\eta -1/2)+d(C+{\log }_2\rho )\text{mbox}where\text{leavevmode}\rho =\frac{(1+\eta {)}^2+ϵ}{\delta -{\eta }^2}.$$

$$\ell \ge T(d,\delta ,\eta )=10+2{\log }_{2}d-{\log }_2(\eta -1/2)+(\eta -1/2+{\log }_2\rho )d.$$

$$D(B)=\prod_{i=1}^{d}\operatorname{covol}\mathopen{}\mathclose{{}\left({B_{[1\ldots i]}}}\right).$$

$$D(B)^{2}\leq\mathopen{}\mathclose{{}\left(d^{2}\|\mathcal{G}_{\gamma}\|_{\textrm{max}}\|L\|_{\textrm{max}}^{2}}\right)^{d(d+1)/2}.$$

$$D(B{)}^2\ge {\lambda }_d({\mathcal{G}}_{\gamma }{)}^{d(d+1)/2}.$$

$$\tau\leq\textrm{O}\mathopen{}\mathclose{{}\left(d^{2}\mathopen{}\mathclose{{}\left(\log(\|L\|_{\textrm{max}}}\right)+\log\mathopen{}\mathclose{{}\left(\|\mathcal{G}_{\gamma}\|_{\textrm{max}}/\lambda_{d}(\mathcal{G}_{\gamma})}\right)+\log(d)}\right).$$

$$D^{\prime}(L)=\mathopen{}\mathclose{{}\left(\prod_{i=1}^{\dim{L}}d_{i}}\right)\cdot\mathopen{}\mathclose{{}\left(\prod_{i,b^{*}_{i}=0}2^{i}}\right).$$

$$\frac{\Vert L_i{\Vert }_1\Vert L_j{\Vert }_1}{|L_i^T{G}_n{L}_j|}.$$

$$er{r}_1\le \frac{\Vert L_i{\Vert }_1\Vert L_1{\Vert }_1}{|L_1^T{G}_n{L}_1|}\le \frac{{\max }_i\Vert L_i{\Vert }_1^2}{{\lambda }_d(G_n)}\le \frac{d{\max }_i\Vert L_i{\Vert }^2}{{\lambda }_d(G_n)}\le \frac{d{\max }_i\Vert b_i{\Vert }^2}{{\lambda }_d(G_n{)}^2}.$$

Full text

This work has been supported in part by the European Union as H2020 Programme under grant agreement number ERC-669891.

Thomas Espitau

Sorbonne Université

LIP 6, CNRS UMR 7606

Paris, France

Antoine Joux

Chaire de Cryptologie de la Fondation SU

Sorbonne Université, Institut de Mathématiques de Jussieu–Paris Rive Gauche

CNRS, INRIA, Univ Paris Diderot.

Campus Pierre et Marie Curie, F-75005, Paris, France

1. Introduction

In a general setting, a lattice $\Lambda$ is a free $\mathbf{Z}$-module of finite rank, endowed with a positive-definite bilinear form on its ambient space $\Lambda\otimes_{\mathbf{Z}}\mathbf{R}$, as presented for instance in [16]. In particular, this definition implies that $\Lambda$ is discrete in its ambient space for the topology induced by the scalar product. This formalism encompasses the well-known Euclidean lattices when taking the canonical scalar product of $\mathbf{R}^{d}$, but also lattices arising from ideals in rings of integers of number fields. The rank of the lattice $\Lambda$ is defined as the dimension of the vector space $\Lambda\otimes_{\mathbf{Z}}\mathbf{R}$. By definition of a finitely-generated free module, there exists a finite set of vectors $b_{1},\ldots,b_{\operatorname{rk}\Lambda}\in\Lambda$ such that $\Lambda=\bigoplus_{i=1}^{\operatorname{rk}\Lambda}b_{i}\mathbf{Z}$. Such a family is called a basis of the lattice and is not unique. In fact, as soon as $\operatorname{rk}\Lambda\geq 2$ there are infinitely many bases of $\Lambda$. Some among those have interesting properties, such as having reasonably small vectors and low orthogonality defects. They are informally called reduced bases and finding them is the goal of lattice reduction.

Numerous algorithms arising in algebraic number theory heavily rely on lattice reduction, for example, the computation of normal forms of integral matrices (see [10] for the Hermite Normal Form and [9] for the Smith Normal Form), class group computations in a number field [7, 2], or even the enumeration of points of small height near algebraic curves [6].

Even for lattices that use the canonical scalar product, there is a deep link with bilinear forms that clearly appears when considering the Gram matrix of a basis $\mathcal{B}=\{{b_{1}},\dotsc,{b_{d}}\}$, that is, the real symmetric matrix \mathcal{G}=\mathopen{}\mathclose{{}\left(\langle{b_{i}},{b_{j}}\rangle}\right)_{i,j}.

The study of these reduction problems is not recent and goes back to the works of Lagrange and Gauss. These early works were expressed in terms of reduction of quadratic form, more precisely integral binary quadratic forms111This can be viewed as the reduction of integral dimension-two lattices. and led to a method often called Gauss’ algorithm. This method can be seen as a 2-dimensional extension of the Euclid algorithm for computing the greatest common divisor of two integers. In 1850, Hermite proved a general upper bound on the length of the shortest vector in a lattice, given as a function of the dimension and of a very important invariant called the determinant, which is defined in Section 2.1. This bound involves the so-called Hermite constant and has recently been rephrased in algorithmic terms [20, Hermite’s Algorithms]. A century later, in 1982, Lenstra, Lenstra and Lovász designed the lll algorithm [14], with the polynomial factorization problem as an application, following the work of Lenstra on integer programming [15]. This algorithm constitutes a breakthrough in the history of lattice reduction algorithm, since it is the first to have a runtime polynomial in terms of the dimension. It was followed by many improvements lowering its complexity or improving the output’s quality.

Current implementations of lll often work with low precision approximations in order to greatly speed-up the computations. Indeed, the algorithm works surprinsingly well even with such reduced precisions, even if some care needs to be taken to avoid infinite loops. Moreover, once the result is obtained, it can verified efficiently as shown in [30].

We propose here an alternative strategy where we not only certify that the end-result is a reduced basis but also that the algorithm followed a valid computation path to reach it. This strongly deviates from other approaches that have been taken to obtain guaranteed lattice reduced basis. At first, this may seems irrelevant. After all, one might claim that a basis satisfying the end conditions of lll is what is desired and that the computation path doesn’t matter. However, as shown in [13] for Siegel-reduced bases, a reduced basis chosen uniformely at random behaves as the worst-case allowed by the final inequalities. By constrast, bases produced by the lll algorithm are usually much better than this worst-case. This argues in favor of trying to follow the algorithm defintion exactly to better understand the phenomenon. In particular, this option might be invaluable for experiments performed toward analyzing this gap.

The present article also relies on Interval Arithmetic, a representation of reals by intervals—whose endpoints are floating-point numbers—that contain them. Arithmetic operations, in particular the basic operations $+,-,\times,\div$ can be redefined in this context. The main interest of this representation lies in its certification property: if real numbers are represented by intervals, the interval resulting from the evaluation of an algebraic expression contains the exact value of the evaluated expression.

For some authors, Interval Arithmetic was introduced by R. Moore in 1962 in his Ph.D. thesis [18]. For others, it can be dated back to 1958, in an article of T. Sunaga [28] which describes an algebraic interpretation of the lattice of real intervals, or even sooner in 1931 as a proposal in the Ph.D. thesis [31] of R.C. Young at Cambridge. Its main asset—calculating directly on sets—is nowadays used to deterministically determine the global extrema of a continuous function [24] or localizing the zeroes of a function and (dis)proving their existence [11]. Another application of Interval Arithmetic is to be able to detect lack of precision at run-time of numerical algorithms, thanks to the guarantees it provides on computations. This can, in particular, be used to design adaptive-precision numerical algorithms.

In the present paper, we propose to transform and generalize the lll algorithm into an adaptive-precision version, which can reduce arbitrary lattices and follows a certified flow of execution. More precisely, it uses Interval Arithmetic to validate the size-reduction and exchange steps that occur within lll.

The interested reader may download an implementation of the algorithm from the webpage http://almacrypt.eu/outputs.php.

Organisation of the paper

In Section 2, we briefly introduce reduction theory and present the l2 variant of the lll algorithm. Section 3 aims at describing the basics of Interval Arithmetic used in Section 4 to handle the problem of representation of real lattices. The framework of this latter section is then used in Section 5 to derive a certified reduction algorithm for real lattices. Section 6 presents an application to algorithmic number theory.

Notations and conventions

General notations

As usual, the bold capitals $\mathbf{Z}$, $\mathbf{Q}$, $\mathbf{R}$ and $\mathbf{C}$ refer respectively to the ring of integers and the fields of rational, real and complex numbers. Given a real number $x$, the integral roundings floor, ceil and round to nearest integer are denoted respectively by $\lfloor x\rfloor,\lceil x\rceil,\lfloor x\rceil$. Note that the rounding operator is ambiguous when operating on half-integers. However, either choice when rounding is acceptable in lattice reduction algorithms. In fact, in this context, it is often enough to return an integer close to $x$, not necessarily the closest.

These operators are extended to operate on vectors and matrices by point-wise composition. The complex conjugation of $z\in\mathbf{C}$ is denoted by the usual bar $\bar{z}$ whereas the real and imaginary parts of a complex $z$ are indicated by respectively $\mathfrak{R}(z)$ and $\mathfrak{I}(z)$. All logarithms are taken in base $2$.

Matrices and norms

For a field $\mathbf{K}$, let us denote by $\mathbf{K}^{d\times d}$ the space of square matrices of dimension $d$ over $\mathbf{K}$, $\textrm{Gl}_{d}(\mathbf{K})$ its group of invertible elements and $\mathcal{S}_{d}(\mathbf{K})$ its subspace of symmetric matrices. For a complex matrix $A$, we write $A^{\dagger}$ for its conjugate transpose. For a vector $v$, we denote by $\|v\|_{\infty}$ its absolute (or infinity) norm, that is the maximum of the absolute value of its entries. We similarly define the matrix max norm $\|B\|_{\textrm{max}}=\max_{(i,j)\in[1\,\cdots\,d]^{2}}|B_{i,j}|$, for any matrix $B$.

Computational setting

The generic complexity model used in this work is the random-access machine (RAM) model and the computational cost is measured in bits operations. \mathcal{M}\mathopen{}\mathclose{{}\left({k}}\right) denotes the complexity of the multiplication of two integers of bit length at most $k$. It is also the cost of the multiplication of two floating-point numbers at precision $k$, since the cost of arithmetic over the exponents is negligible with regards to the cost of arithmetic over the mantissae.

2. Basics of Lattice Reduction

2.1. Orthogonalization

Let us fix an Euclidean space $(E,\langle{\cdot},{\cdot}\rangle)$, i.e. a real vector space $E$ together with a positive-definite bilinear form $\langle{\cdot},{\cdot}\rangle:E\times E\to\mathbf{R}$. As usual, two vectors $x,y\in E$ are said to be orthogonal—with respect to the form $\langle{\cdot},{\cdot}\rangle$—if $\langle{x},{y}\rangle=0$. More generally a family of vectors is orthogonal if its elements are pairwise orthogonal.

Now consider $S=(b_{1},\dots,b_{r}),$ a family of linearly independent vectors of $E$. The flag $\mathcal{F}_{S}$ associated to $S$ is the finite increasing chain of subspaces:

[TABLE]

The orthogonal complement $S^{\bot}$ is defined as the subspace $\{x\in E\leavevmode\nobreak\ |\leavevmode\nobreak\ \forall i,\leavevmode\nobreak\ \langle{x},{b_{i}}\rangle=0\}$. Denote by $\pi_{i}$ the orthogonal projection on $(b_{1},\dotsc,b_{i-1})^{\bot}$, with the convention that $\pi_{1}$ is the identity map. The Gram–Schmidt orthogonalization process—shorthanded as gso—is an algorithmic method for orthogonalizing $S$ while preserving its flag. It constructs the orthogonal set S^{*}=\mathopen{}\mathclose{{}\left(\pi_{1}(b_{1}),\ldots,\pi_{r}(b_{r})}\right). The computation $S^{*}$ can be done inductively as follows:

[TABLE]

Define the Gram matrix, associated to a family of vectors $S=(b_{1},\ldots,b_{r})$, as the symmetric matrix of scalar products: \mathcal{G}_{S}=\mathopen{}\mathclose{{}\left(\langle{b_{i}},{b_{j}}\rangle}\right)_{(i,j)\in[1\,\cdots\,r]^{2}}. The (co)volume of $S$, also called its determinant, is defined as the square root of the Gram determinant $\det\mathcal{G}_{S}$. It can be easily computed from the Gram-Schmidt vectors $S^{*}$ as:

[TABLE]

2.2. Lattices and reduction

Definition 2.1**.**

A (real) lattice $\Lambda$ is a finitely generated free $\mathbf{Z}$-module, endowed with a positive-definite bilinear form $\langle{\cdot},{\cdot}\rangle$ on its ambient space $\Lambda\otimes_{\mathbf{Z}}\mathbf{R}$.

By definition of the tensor product, there is a canonical injection that sends a vector $v$ to $v\otimes 1$ in the ambient space and preserves linear independence. Thus, the rank of $\Lambda$ as a $\mathbf{Z}$-module, is equal to the dimension of the vector space $\Lambda\otimes_{\mathbf{Z}}\mathbf{R}$.

Denoting by $d$ the rank of the lattice, a basis of $\Lambda$ is a family $b_{1},\ldots,b_{d}$ of elements of $\Lambda$ such that $\Lambda=\bigoplus_{i=1}^{d}b_{i}\mathbf{Z}$.

In the sequel, we identify $\Lambda$ with its canonical image $\Lambda\otimes 1$ and thus view the lattice as an additive subgroup of its ambient space $\Lambda\otimes_{\mathbf{Z}}\mathbf{R}$. When the context makes it clear, we may omit to write down the bilinear form associated to a lattice $\Lambda$. Throughout this section, $\|\cdot\|$ stands for the Euclidean norm induced by $\langle{\cdot},{\cdot}\rangle$, unless stated otherwise. As usual, any two bases $(b_{1},\ldots,b_{d})$ and $(b^{\prime}_{1},\ldots,b^{\prime}_{d})$ of $\Lambda$ are related by a unimodular transformation, i.e., a linear transformation represented by a $d\times d$ integer matrix of determinant $\pm 1$.

Lemma 2.2**.**

A lattice $\Lambda$ is discrete for the topology induced by the given norm on its ambient space. I.e., there exists a real $\epsilon_{\Lambda}>0$ such that for any pair $(x,y)$ of elements of $\Lambda$ with $x\neq y$ we have:

[TABLE]

The largest possible value for $\epsilon_{\Lambda}$ in the above inequality is equal to the norm of the shortest non-zero vector of $\Lambda$, which is traditionally called the first minimum or the minimum distance of the lattice and denoted by $\lambda_{1}(\Lambda)$.

Proof.

Let $\mathcal{B}=(b_{1},\ldots,b_{d})$ be a basis of $\Lambda$. Let $\mathcal{B}^{*}=(\pi_{1}(b_{1}),\ldots,\pi_{d}(b_{d}))$ be the orthogonal basis obtained by applying Gram-Schmidt orthogonalization to the canonical image of $\mathcal{B}$ in $\Lambda\otimes_{\mathbf{Z}}\mathbf{R}$. This orthogonalization is taken using as scalar product the given bilinear form.

Assume by contradiction that there exist pairs of distinct vectors with the norm of their difference arbitrarily small. Since the difference is also an element of $\Lambda$, there are non-zero elements of arbitrarily small norm in $\Lambda$. For any integer $i>0$, choose a vector $x_{i}$ in $\Lambda$ with $\|x_{i}\|^{2}\leq 2^{-i}.$ Decompose $x_{i}$ in the basis $\mathcal{B}^{*}$ as $x_{i}=\sum_{j=1}^{d}\chi_{i}^{(j)}\pi_{j}(b_{j}).$ For any pair of integers $i$, $j$ we see that $|\chi_{i}^{(j)}|^{2}\,\|\pi_{j}(b_{j})\|^{2}\leq\|x_{i}\|^{2}\leq 2^{-i}.$ As a consequence, each sequence $\chi^{(j)}$ converges to zero. Multiplying by the basis-change matrix, we see that the coordinates of $x_{i}$ in the basis $b_{1},\ldots,b_{d}$ also converge to zero. Since these coordinates are integral, the sequences are ultimately constant and $x$ is also ultimately constant (and null). This contradicts the choice of $x_{i}$ as a non-zero element. ∎

2.3. The LLL reduction algorithm

In 1982, Lenstra, Lenstra and Lovász [14] proposed a notion called lll reduction and a polynomial-time algorithm that computes an lll-reduced basis from arbitrary basis of the same lattice. Their reduction notion is formally defined as follows:

Definition 2.3** (lll reduction).**

A basis $\mathcal{B}=(b_{1},\ldots,b_{d})$ of a lattice is said to be $\delta$-lll-reduced for a parameter $1/4<\delta<1$, if the following conditions are satisfied:

[TABLE]

[TABLE]

In order to find a basis satisfying these conditions, it suffices to iteratively modify the current basis at any point where one of these conditions is violated. This yields the simplest version of the lll algorithm. As in [14], it is only defined for full-rank sublattice of $\mathbf{Z}^{d}$. It was remarked by Lovász and Scarf in [17] that the same algorithm also works with an arbitrary integral-valued scalar product. The method can be extended to deal with lattices described by a generating family rather than by a basis [23].

2.3.1. Decrease of the potential and complexity.

The algorithm can only terminate when the current lattice basis is lll-reduced. Moreover, as shown in [14], it terminates in polynomial time when $\delta<1$. Indeed, consider the (square of the) product of the covolumes of the flag associated with a basis: $\prod_{i=1}^{d}\|\pi(b_{i})\|^{2(d-i+1)},$ which is often called its potential. This value decreases by a factor at least $\delta^{-1}$ in each exchange step and is left unchanged by other operations. Indeed:

• •

The flag is not modified by any operation other than swaps.

• •

A swap between ${b_{k}}$ and ${b_{k-1}}$ only changes the sublattice spanned by the first $k-1$ vectors. The corresponding covolume $\prod_{i=1}^{k-1}\|\pi(b_{i})\|^{2}$ decreases by a factor at least $\delta^{-1}$ and so does the potential.

Since the total number of iterations can be bounded by twice the number of swaps plus the dimension of the lattice, this suffices to conclude that it is bounded by \textrm{O}\mathopen{}\mathclose{{}\left(d^{2}\log\|B\|_{\textrm{max}}}\right) where $B$ is the matrix of the initial basis.

As the cost of a loop iteration is of \textrm{O}\mathopen{}\mathclose{{}\left(d^{2}}\right) arithmetic operations on rational coefficients of length at most \textrm{O}\mathopen{}\mathclose{{}\left(d\log\|B\|_{\textrm{max}}}\right), the total cost in term of arithmetic operations is loosely bounded by \textrm{O}\mathopen{}\mathclose{{}\left(d^{6}\log^{3}\|B\|_{\textrm{max}}}\right). By being more precise in the majoration of the bit length of the integers appearing in lll, this analysis can be improved. Kaltofen in [12] bounds the complexity by

[TABLE]

2.3.2. A bound on the norm of reduced elements

Proposition 1**.**

Let $1/4<\delta<1$ be an admissible lll parameter. Let $(b_{1},\ldots,b_{d})$ be a $\delta$-lll reduced basis of rank-$d$ lattice $(\Lambda,\langle{\cdot},{\cdot}\rangle)$. Then for any $1\leq k\leq d$:

[TABLE]

Note that this is an easy generalization of the bound on the norm of $b_{1}$ which is given in most texts. It appears among other related inequalities in [22]. For completeness, a proof is given in Appendix.

2.3.3. Floating point representation

The total cost of the lll algorithm is dominated by the computation to handle arithmetic on rational values. A first idea of De Weger [5] to overcome this issue is to avoid the use of denominators by multiplying all the quantities by their common denominator. This is slightly more efficient in practice but doesn’t improve the asymptotics. Another idea is to remark that the norms of the rational values remain small and to try to use approximations instead of exact values. However, directly replacing rationals in the lll algorithm by floating-point approximations leads to severe drawbacks. The algorithm might not even terminate, and the output basis is not guaranteed to be lll-reduced.

The first provable floating-point version of the algorithm is due to Schnorr in [26], with complexity \textrm{O}\mathopen{}\mathclose{{}\left(d^{4}\log(\|B\|_{\textrm{max}})\mathcal{M}(d+\log\|B\|_{\textrm{max}})}\right). One of the key ingredients to achieve this reduction is to slightly relax the definition of the size-reduction, in order to compensate for the approximation errors introduced by the use of floating-point arithmetic. We call admissible any parameters $(\delta,\eta)$ satisfying $1/4<\delta<1$, and $1/2<\eta<\sqrt{\delta}$ and define:

Definition 2.4** ($(\delta,\eta)$-lll reduction).**

Let $(\delta,\eta)$ be admissible parameters. A basis $\mathcal{B}$ of a lattice is said to be $(\delta,\eta)$-LLL-reduced if the following condition is satisfied:

[TABLE]

together with the Lovász condition, which is kept unchanged from Definition 2.3.

Using naive multiplication, the cost of Schnorr’s algorithm is cubic in the size of the numbers, i.e. in $\log(\|B\|_{\textrm{max}})$. The introduction of approximate size reduction removes the need to know with extreme precision values close to half-integers. Instead, approximate size reduction of such values can be achieved by rounding either up or down in an arbitrary (possibly randomized) manner. In our pseudo-code, we use a function called $\eta$-Closest-Integer to achieve this rounding, returning an integer at distance at most $\eta$ of the function’s argument.

2.4. The $L^{2}$ algorithm

The l2 algorithm is a variant of Schnorr-Euchner version [27] of lll. By contrast with the original algorithm, l2 computes the gso coefficients on the fly as they are needed instead of doing a full orthogonalization at the start. It also uses a lazy size reduction inspired by the Cholesky factorization algorithm. These optimizations yield an improved lattice reduction with running time

[TABLE]

As usual in lattice reduction, while performing the Gram-Schmidt orthogonalization of $\mathcal{B}$, we also compute QR-decomposition of $B$ into $B^{*}\cdot M$ where $B^{*}$ is the matrix representing the \mathopen{}\mathclose{{}\left(\pi_{i}(b_{i})}\right)_{1\leq i\leq d}, and $M$ is the upper unitriangular matrix, whose coefficients with $j\geq i$ are $M_{i,j}=\frac{\langle{b_{j}},{\pi_{i}(b_{i})}\rangle}{\|\pi_{i}(b_{i})\|^{2}}$. Thus, the Gram matrix associated to the basis, i.e., $G=B^{T}B$ satisfies:

[TABLE]

where $D$ is a diagonal matrix whose entries are ${\|\pi_{i}(b_{i})\|^{2}}$. We denote by $R$ the matrix $D\cdot M$, and thus have $G=R^{T}\cdot M=M^{T}\cdot R.$

We give the pseudo-code of the Lazy Size-Reduction procedure as Algorithm 2 and of the l2 algorithm as Algorithm 3. Both use classical formulas relating $R$, $M$ and $B^{*}$ to perform the computations.

2.4.1. Precision required.

The precision required by the l2-Algorithm is

[TABLE]

bits for any $\epsilon>0$, i.e., almost linear in the dimension of the lattice. Moreover, as discussed in [21], it appears that—even though this bound can be shown to be sharp by specific examples—experiments indicate that the number of bits required on average is, in fact, lower.

This phenomenom is well-known and is often used in existing algorithms and softwares in the form of a compute-and-verify paradigm. For example, this is default strategy of the well-known FPLLL [29]. It relies on the fact that verifying that a lattice basis is indeed reduced is much less costly than the reduction itself, as shown in [30]. In addition, it is necessary to take several conservative measures in order to prevent the implementation to enter potentially infinite loops.

The approach we propose deviates from this paradigm. Instead of guaranteeing the end-result, we want to make sure that the whole computation follows the mathematical definition of the algorithm. With low-precision approximations, it is unclear how this could be done. However, interval-arithmetic offers a neat solution to achieve this goal.

3. Interval Arithmetic and its certification property

Interval arithmetic is a representation of reals by intervals that contain them. For instance, one can specify a value $x$ with an error $\epsilon$ by giving an interval of length $\epsilon$ containing $x$. For example, the constant $\pi$ can be represented with an error of $10^{-2}$ by the interval $[3.14,3.15]$. Interval arithmetic is crucial in the context of certified numerical computations, where reals can only be represented with finite precision. For more details, the interested reader can consult an extensive reference, such as [19].

In the following, we denote by $\underline{x}$ a closed interval $[\underline{x}^{-},\underline{x}^{+}]$. We define its diameter as the positive real $\underline{x}^{+}-\underline{x}^{-}$ and its center as the real $\frac{1}{2}(\underline{x}^{+}+\underline{x}^{-})$.

Given a real-valued function $f(x_{1},\ldots,x_{n})$ an interval-arithmetic realization of $f$ is an interval-valued function $F$ such that the interval $F(\underline{x_{1}},\ldots,\underline{x_{n}})$ contains all the values $f(x_{1},\ldots,x_{n})$ for $(x_{1},\ldots,x_{n})$ in $\underline{x_{1}}\times\cdots\times\underline{x_{n}}$.

If $F$ always returns the smallest possible interval, it is called a tight realization, otherwise it is called loose. In practice, tight realizations can only be achieved in very simple specific cases. However, even a loose realization can suffice to certify the correctness of a computation.

Another important property of interval arithmetic is that it can be used to compare numbers in a certified way, as long as the intervals that represent them are disjoint.

3.1. Some useful interval-arithmetic realizations

3.1.1. Integral representation of fixed length

A first convenient way to represent reals at finite precision is to use integers as an approximate representation.

Definition 3.1** (Integral representation of reals).**

Let $x\in\mathbf{R}$ be an arbitrary real number and $n\geq 0$ a non-negative integer. Define an integral representation at accuracy222 We use here the denomination of “accuracy” instead of “precision” to avoid confusions with the floating-point precision as defined in paragraph 3.1.3.  $n$ as an interval of diameter $2$:

[TABLE]

together with a guarantee that $2^{n}x$ belongs to $\underline{x}_{n}$.

This representation is very compact, since it only requires to store the center $X_{n}$ of the interval using $n+\lceil\log x\rceil$ bits. However, computing with this form of representation is not convenient. As a consequence, we only use it to represent immutable values and we convert to a different representation for computations. The reason for using the interval \mathopen{}\mathclose{{}\left[X_{n}-1,X_{n}+1}\right] of diameter $2$ rather than \mathopen{}\mathclose{{}\left[X_{n}-1/2,X_{n}+1/2}\right] (of diameter $1$) is that when $2^{n}x$ is very close to a half-integer, it remains possible to easily provide a valid value for $X_{n}$ without computing extraneous bits of the representation of $x$.

3.1.2. Fixed-point representations

In the context of lattice reduction, it is useful to compute linear combinations with exact integral coefficients. In order to do that with approximate values initially given by centered integral representation, it is possible to use a fixed-point representation.

Definition 3.2** (Fixed point representation of reals).**

Let $x\in\mathbf{R}$ be an arbitrary real number and $n\geq 0$ a non-negative integer. Define a fixed-point representation at accuracy $n$ of radius $\delta$ as an interval:

[TABLE]

together with a guarantee that $2^{n}x$ belongs to $\underline{x}_{n}$.

It is easy to add or subtract such intervals by doing the computation on the center and by adding the two radii. It is also easy to multiply by an exact integer by multiplying the center by the integer and the radius by its absolute value. Integral representations are a special case of fixed-point representations, with radius equal to $1$.

3.1.3. Floating-point representation.

Another way to handle real values is to use floating point representations of the two bounds of each interval. For example, if we denote by $\lfloor x\rfloor_{n}$ and $\lceil x\rceil_{n}$ respectively the largest floating-point number below $x$ and the lowest floating-point number above $x$ written with $n$ bits, the tightest floating-point representation of $x$ with $n$ bits of precision is the interval I_{n}(x)=\mathopen{}\mathclose{{}\left[\lfloor x\rfloor_{n},\lceil x\rceil_{n}}\right].

With such a representation, it becomes possible to create a realization of the elementary operations by using careful rounding when computing approximations of the bounds of the resulting interval, as shown in Figure 1. When speaking of the precision of such a representation, we simply refer to the common floating-point precision of the upper and lower bounds.

Once the elementary operations are available, they can be used to implement certified versions of any function that can classically be computed with floating point arithmetic.

4. Approximate lattices

The need to reduce lattices given by approximations, especially for number-theoretic applications as been known for long. In particular, Buchmann gives in [3] a bound on the required precision to achieve this goal by using a direct approximation of the input basis. However, this bound is computed in terms of a quantity called the defect that can be very large and also involves the first minimum of the lattice.

Using interval arithmetic, it becomes possible to get finer control on the precision required to perform the lattice reduction, even with approximate lattices.

4.1. Approximate representation of a positive-definite matrix

A matrix with real entries can easily be represented with the integral representation from Definition 3.1, using the same accuracy for all of its entries.

Definition 4.1** (Matrix integral representation).**

Let $A=(a_{i,j})_{i,j}\in\mathbf{R}^{d\times d}$ be an arbitrary real matrix of dimension $d$ and $n>0$ be a fixed positive integer. A matrix of intervals

[TABLE]

where each $\underline{a_{i,j}}_{n}$ is an integral representation of $a_{i,j}$ is said to integrally represent $A$ at accuracy $n$.

We may omit the subscript $n$ when the accuracy is clear from the context. Given a matrix $A$, and a matrix $B\in\underline{A}_{n}$, there exists a unique $d\times d$ matrix $\Delta$ with entries in $[-2,2]$ such that $B=2^{n}A+\Delta$.

In particular, we may apply this representation to symmetric matrices. In that case, we obtain the following useful lemma:

Lemma 4.2**.**

Let $S=(s_{i,j})_{i,j}\in\mathcal{S}_{d}(\mathbf{R})$ be a symmetric matrix of dimension $d$ and $\underline{S}_{n}$ an integral representation of $S$ at accuracy $n$. Then, for any symmetric matrix $S^{\prime}$ in $\underline{S}_{n}$, we have:

[TABLE]

where $\lambda_{d}(T)$ denotes the smallest eigenvalue of a $d$-dimensional symmetric matrix $T$.

Proof.

This is a direct consequence of Weyl’s inequalities for Hermitian matrices and of the relation $S^{\prime}=2^{n}S+\Delta$, where $\Delta$ is real symmetric with entries in $[-2,2]$. Note that the eigenvalues of $\Delta$ all belong to $[-2d,2d]$. ∎

4.2. Representation of lattices

In order to represent arbitrary lattices, we first need a description of their ambient space. We simply describe the ambient space $V$ of dimension $d$ by providing a basis $\gamma=(\gamma_{1},\ldots,\gamma_{d})$. Then, the scalar product $\langle{\cdot},{\cdot}\rangle$ on $V$ can be encoded by a Gram matrix \mathcal{G}_{\gamma}=\mathopen{}\mathclose{{}\left(\langle{\gamma_{i}},{\gamma_{j}}\rangle}\right)_{(i,j)\in[1\,\cdots\,d]^{2}}.

When the Gram matrix $\mathcal{G}_{\gamma}$ is integral, this already is a standard description of the lattice $\Gamma$ spanned by $\gamma$. This representation appears in particular in [4, Proposition 2.5.3]. We now extend this in order to represent bases and generating families of arbitrary sublattices of $\Gamma$. Let $\Lambda$ be a rank $r\leq d$ sublattice of $\Gamma$ given by a generating family \ell=\mathopen{}\mathclose{{}\left(\ell_{1},\ldots,\ell_{p}}\right). Since any vector in $\ell$ belongs to $\Gamma$, it can be expressed with integral coordinates in the basis $\gamma$. As a consequence, we can represent $\ell$ by a $p\times d$ integral matrix $L$. Moreover, the knowledge of $\mathcal{G}_{\gamma}$ allows us to easily compute the scalar product of any pair of vectors in $\Lambda$.

All this leads to the following definition:

Definition 4.3** (Approximate representation of a lattice).**

Let $\mathcal{G}_{\gamma}$ and $L$ be as above and $n$ be a non-negative integer. Denote by $G$ the matrix of centers of an integral representation $\underline{\mathcal{G}_{\gamma}}_{n}$at accuracy $n$ of the Gram matrix $\mathcal{G}_{\gamma}$. Then the pair $(G,L)\in\mathbf{Z}^{d\times d}\times\mathbf{Z}^{p\times d}$ of integral matrices is said to represent at accuracy $n$ the lattice $\Lambda$ in the basis $\gamma$ of $\Gamma$.

4.2.1. Computation of the inner product in Interval Arithmetic.

Let $a$ and $b$ be two vectors of $\Lambda$ described by their vectors $A$ and $B$ of coordinates in the basis $\gamma$. We know that:

[TABLE]

Thus:

[TABLE]

This directly gives an interval representation of $\langle{a},{b}\rangle$.

4.3. Lattice reduction of approximate lattices

Suppose now that the Gram matrix \mathcal{G}_{\gamma}=\mathopen{}\mathclose{{}\left(\langle{\gamma_{i}},{\gamma_{j}}\rangle}\right)_{(i,j)\in[1\,\cdots\,d]^{2}} representing the inner product of the ambient space $\Gamma\otimes_{\mathbf{Z}}\mathbf{R}$ in the basis $\gamma$ is given indirectly by an algorithm or an oracle $\mathcal{O}_{\gamma}$ that can compute each entry at any desired accuracy. We can restate the definition of a reduced basis in this framework as:

Definition 4.4** ($(\delta,\eta)$-lll reduction).**

Let $(\delta,\eta)$ be admissible lll parameters. Given an integral matrix $L\in\mathbf{Z}^{p\times d}$ which describes the vectors of a basis of a lattice $\Lambda$ in the basis $\gamma$, we say that $(\mathcal{G}_{\gamma},L)$ is a $(\delta,\eta)$-lll reduced basis of $\Lambda$ if and only if there exists an $n_{0}>0$ such that for any $n\geq n_{0}$ there exists a pair $({\underline{G}}_{n},L)$, where ${\underline{G}}_{n}$ is an integral representation of $\mathcal{G}_{\gamma}$ at accuracy $n$, which is a $(\delta,\eta)$-lll reduced basis.

The computational problem associated with reduction theory can then be written as:

Problem** (Lattice Reduction for approximate representation).**

Let $\delta,\eta$ be admissible lll parameters. Given as input an algorithm or oracle to compute $\mathcal{G}_{\gamma}$ at arbitrary precision and an integral matrix $L\in\mathbf{Z}^{p\times d}$ that describes the vectors of a generating family of a lattice $\Lambda$ in the basis $\gamma$: find a basis $L^{\prime}$ of $\Lambda$ such that $(\mathcal{G}_{\gamma},L^{\prime})$ is a $(\delta,\eta)$-lll reduced basis in the sense of Definition 4.4.

Note that using interval arithmetic it suffices to check the $(\delta,\eta)$-lll reduction condition at accuracy $n_{0}$ to be sure it holds at any larger accuracy. Indeed, an integral representation that satisfies the condition can be refined into a more precise integral representation by scaling up the integer representing the center by an adequate power of two. This refined representation continues to satisfy the condition.

4.3.1. Accuracy of representation and space complexity.

Let $({\underline{G}}_{n},L)$ be an integral representation of $\Lambda$, at accuracy $n$. Then, the magnitude of the entries of $G$ is $2^{n}$ times the magnitude of the entries of $\mathcal{G}_{\gamma}$. Thus, ${\underline{G}}_{n}$ can be encoded using \textrm{O}\mathopen{}\mathclose{{}\left(d^{2}(n+\log\|\mathcal{G}_{\gamma}\|_{\textrm{max}})}\right) bits.

5. Generalized LLL reduction with Interval Arithmetic

In this Section, we adapt lattice reduction algorithms to our setting. More precisely, we represent the information related to Gram-Schmidt vectors by interval arithmetic using a floating-point representation as described in Section 3.1.3. For the representation of the lattice itself, we consider two cases: either the underlying Gram matrix is integral, or it is given by an approximate integral representation as in Section 4.1. In the latter case, our algorithm also asks for representations with higher accuracy until it is sufficient to yield a reduced basis for the given lattice. The canonical case with the standard Euclidean scalar product is achieved by setting the Gram matrix to the (exact) identity matrix.

5.1. Interval Arithmetic $L^{2}$ reduction with fixed precision.

We first consider the simplified case where the lattice representation is fixed. It can be either exact or approximate with a given accuracy. In both cases, we fix a basis $\gamma=(\gamma_{1},\ldots,\gamma_{d})$ and a representation of a lattice $\Lambda$ in this basis. It is respectively an exact integral representation $(G,L)$ or an approximate representation $({\underline{G}}_{n},L)$ at accuracy $n$ of $(\mathcal{G}_{\gamma},L)$.

5.1.1. Using Interval Arithmetic in lll.

We now modify the l2 algorithm of [21] in a few relevant places to make use of interval arithmetic instead of floating-point arithmetic for the Gram-Schmidt-related values. Since the description of the lattice $\Lambda$ is already using intervals, it seems natural to use interval arithmetic in the lattice reduction algorithm. For completeness, when the input Gram matrix is exact, we make the updates to the Gram-Schmidt orthogonalized matrix used by lll explicit in the algorithm (except the simple displacements). This also emphasizes a subtle difference with the case of an approximate input Gram matrix. Indeed, in that case, we update the gso-values but recompute the errors rather than relying on the interval arithmetic to do it. This is important to gain a fine control on the error growth during updates.

In addition, when using the technique from [23] to be able to deal with lattices given by a generating family instead of a basis, we make a slightly different choice than in [21]. Instead of moving the zero vectors that are encountered during the computation during the reduction to the start of the basis, we simply remove them. Note that with an approximate matrix, if we discover a non-zero vector whose length is given by an interval containing [math], it is not possible to continue the computation. This means that the accuracy of the input is insufficient and we abort. The core modification with interval arithmetic appears while testing the Lovász condition. If it is not possible to decide whether the test is true or false because of interval overlap, we also abort due to lack of precision. To be more precise, when testing the Lovász condition, we also need to check that the corresponding $\mu$ coefficient is indeed smaller than $\eta$. The reason for this is that, when called with insufficient precision, the Lazy reduction routine may fail to ensure that property.

In addition, if a negative number occurs when computing the norm of a vector, it means that the given Gram matrix is not positive-definite and the algorithm returns an error accordingly.

5.1.2. Internal precision in the exact-input case.

For the classical l2 algorithm, Section 2.4.1 states that the precision that is needed for the computations only depends on the dimension of the lattice. It is natural to ask a similar question about the algorithm $\overline{\textrm{{ll}}}$: can the required internal accuracy be bounded independently of the entries appearing in the matrices $G$ and $L$. When $G$ is exact, i.e., integral, the adaptation is straightforward and we obtain the following result.

Theorem 5.1**.**

Let $(\delta,\eta)$ be admissible lll parameters. Let $c>\log\frac{(1+\eta)^{2}}{\delta-\eta^{2}}$ and let $(\Lambda,\langle{\cdot},{\cdot}\rangle)$ denote a rank-$d$ lattice, exactly described by the pair $(G,L)$. Let $B$ denotes the maximum entry in absolute value in $L^{T}GL$. Then, the $\overline{\textrm{{ll}}}$ of Figure 5 used with \ell=cd+o\mathopen{}\mathclose{{}\left(d}\right) outputs a $(\delta,\eta)$-lll-reduced basis in time \textrm{O}\mathopen{}\mathclose{{}\left(d^{3}\log{B}(d+\log{B})\mathcal{M}\mathopen{}\mathclose{{}\left({d}}\right)}\right). Furthermore, if $\tau$ denotes the number of main loop iterations, the running time is \textrm{O}\mathopen{}\mathclose{{}\left(d(\tau+d\log{dB})(d+\log{B})\mathcal{M}\mathopen{}\mathclose{{}\left({d}}\right)}\right).

In fact, the bound on $\ell$ is made explicit in [21]. More precisely, it states that for any arbitrary $C>0$ and an $\epsilon\in]0,1/2]$, it suffices to have:

[TABLE]

For example, choosing $C=\epsilon=\eta-1/2$ it suffices to have:

[TABLE]

When $\delta$ is close to $1$ and $\eta$ to $1/2$, the constant before $d$ becomes smaller than $1.6$.

5.1.3. Dealing with approximate inputs.

When dealing with lattices given in an approximate form, i.e., by a representation $({\underline{G}}_{n},L)$ at accuracy $n$ of $(\mathcal{G}_{\gamma},L)$, the analysis of the algorithms differs in three main places:

• •

When bounding the number of rounds $\tau$, we can no longer assume that the potential is an integer. As a consequence, in order to keep a polynomial bound on $\tau$, we need to provide a lower bound on the possible values of the potential, rather than rely on the trivial lower bound of $1$ for an integral-valued potential.

• •

Since the notion of lll-reduction is only well-defined for a positive definite $G$, we need to make sure that ${\underline{G}}_{n}$ is positive-definite during the algorithm. Otherwise, it should output an error; Algorithm 5 returns an error that ${\underline{G}}_{n}$ is incorrect whenever it encounters a vector with a negative norm.

• •

When ${\underline{G}}_{n}$ is approximate, the scalar products between lattice vectors can no longer be exactly computed. Thus, we need to able to make sure that the errors are small enough to be compatible with the inner precision used for Gram-Schmidt values. At first glance, this might seem easy. However, when using update formulas to avoid recomputation of scalar products, the estimates on errors provided by interval arithmetic can grow quite quickly. In fact,it would prevent the update strategy from working. The key insight is to remark that since the centers of the intervals are represented by integers, any computation on them is exact and we can use update formulas to compute them. However, it is essential to recompute the radii of the intervals, i.e., the errors, to prevent them from growing too quickly.

Number of rounds

Since interval arithmetic allows up to emulate exact computations as long as no failures are detected, we can analyze the number of rounds by assuming that all computations on non-integral values are done using an exact arithmetic oracle. In this context, the number of rounds can be studied by considering the potential as usual. Remember that the initial setting where lll operates on a basis the potential is defined as

[TABLE]

The key argument is that it decreases by a multiplicative factor whenever an exchange is performed.

However, in our context, the starting upper bound and the ending lower bound are different from the integer lattice setting. The initial upper bound needs to account from the presence of the positive definite matrix. So if the lattice is described by a pair $(\mathcal{G}_{\gamma},L)$ the upper bound becomes:

[TABLE]

More importantly, it is no longer possible to claim that the potential is an integer. Instead, we derive a lower bound by considering the smallest eigenvalue of $\mathcal{G}_{\gamma}$ and find:

[TABLE]

As a consequence, if we let $\tau$ denote the number of rounds of the algorithm, we can conclude that:

[TABLE]

When the lattice is given by a generating family $L$ rather than a basis $B$, we need a slightly different invariant. Following [21], we define $d_{i}$ to be the product of the first $i$ non-zero values $\|b^{*}_{j}\|$. Note that they are not necessarily consecutive, since zeroes may occur anywhere. We then let:

[TABLE]

This generalized potential is needed for the proof of Theorem 5.2. Note that, for lattices given by a basis, the two definitions coincide.

Necessary accuracy for the scalar products

In order to preserve the correctness of the algorithm when computing with internal precision $\ell$, we need to check that all conversions of scalar product values, using the calls to ConvertToFPinterval in Algorithms 4 and 5, have sufficient precision. For a pair of lattice elements, described by vectors $L_{i}$ and $L_{j}$, the relative precision on the value of their scalar product is:

[TABLE]

When the vectors are close to orthogonal with respect to the scalar product given by $G_{n}$, the error can be arbitrarily large. However, by carefully following the analysis of Theorem 3 in [21, Section 4.1], we can show that this Theorem remains true in our context. This suffices to ensure the correctness part of Theorem 5 of [21]. The first check is to verify that quantity called $err_{1}$ in the proof of the Theorem remains upper bounded by $2^{-\ell}$. Since the value is defined as the error on the scalar product of the vectors number $i$ and $1$ divided by the norm of the first vector, we have:

[TABLE]

Thus:

[TABLE]

As a consequence, it suffices to have:

[TABLE]

l2 with approximate inputs.

To complete the above properties on the number of rounds and necessary accuracy, it suffices to remark that the only additional line of code in the approximate l2 is the recomputation of interval radii on line 10. Since it suffices to know the $\ell$ high-order bits of the values, this recomputation can fully be done using arithmetic on $\ell$. Indeed, during the computations of $\|L_{i}\|_{1}$ no cancellation occurs. As a consequence, we get the following adaptation of Theorem 5.1. For completeness, we give here the case where the lattice is initially given by a generating family of $p$ vectors, has rank $d$ and lives in an ambient space of dimension $D$.

Theorem 5.2**.**

Let $(\delta,\eta)$ be such that $1/4<\delta<1$ and $1/2<\eta<\sqrt{\delta}$. Let $c>\log\frac{(1+\eta)^{2}}{\delta-\eta^{2}}$. Assume that we are given as input $(\Lambda,\langle{\cdot},{\cdot}\rangle)$ a rank-$d$ lattice $(\mathcal{G},L)$ described by $p\geq d$ generating vectors in a ambient space of dimension $D\geq d$. Further assume that it is approximately represented at accurary $N$ by the pair $(\underline{G_{N}},L)$ and let $B$ denote the maximum entry in absolute value in $L^{T}\mathcal{G}L$. Let \ell=cd+o\mathopen{}\mathclose{{}\left(d}\right) and

[TABLE]

Then, the $\overline{\textrm{{ll}}}$ of Figure 5 outputs a $(\delta,\eta)$-lll-reduced basis in time

[TABLE]

Furthermore, if $\tau$ denotes the number of main loop iterations, the running time is \textrm{O}\mathopen{}\mathclose{{}\left(DN\mathopen{}\mathclose{{}\left(dN+\tau}\right)\mathcal{M}\mathopen{}\mathclose{{}\left({d}}\right)}\right).

5.2. $L^{2}$ reduction with adaptive precision and accuracy.

5.2.1. Adaptive precision.

Since by construction the $\overline{\textrm{{ll}}}$ Algorithm can detect that the choice for internal precision $\ell$ is insufficient to correctly reduce the lattice $\Lambda$. The procedure can be wrapped in a loop that geometrically increases precision $\ell$ after each unsuccessful iteration. This yields an adaptive precision reduction algorithm adaptive-lll. Since the complexity of floating-point multiplication is superlinear, the use of a geometric precision growth guarantees that the total complexity of this lattice reduction is asymptotically dominated by its final iteration.333In practice, for lattices of rank few hundreds it appears nonetheless that the computational cost of the previous iterations lies between $20\%$ and $40\%$ of the total cost.

Moreover, the cost of operations in the floating-point realization of interval arithmetic is at most four times the cost of floating-point arithmetic at the same precision. Depending on the internal representation used, this constant can even be improved. As a consequence, for lattices that can be reduced with a low-enough precision, it can be faster to use interval arithmetic than floating-point arithmetic with the precision required by the bound from Section 2.4.1.

5.2.2. Adaptive accuracy.

We now turn to the setting of Section 4.3, where an algorithm or oracle $\mathcal{O}_{\gamma}$ can output an integral representation of the Gram matrix \mathcal{G}_{\gamma}=\mathopen{}\mathclose{{}\left(\langle{\gamma_{i}},{\gamma_{j}}\rangle}\right)_{(i,j)\in[1\,\cdots\,r]^{2}} at arbitrary accuracy $n$. In that context, we need to determine both the necessary accuracy and internal precision. When running Algorithm 5 with some given accuracy and precision, three outcomes are possible:

• •

Either the reduction terminates in which case the lattice is lll-reduced, which implies that both accuracy and precision are sufficient.

• •

The Lovász condition fails to be tested correctly, which indicates an insufficient precision. In that case, we need to test whether the precision is lower than theoretical bound $T(d,\delta,\eta)$ given after Theorem 5.1 or not. In the latter case, we know that the accuracy needs to be increased.

• •

The algorithm detects a non-zero vector whose norm is given by an interval containing 0. This directly indicates insufficient accuracy.

Depending on the result of Algorithm 5, we increase the precision or the accuracy and restart. The corresponding pseudo-code is given in Algorithm 6. Since the precision and accuracy both follow a geometric growth, the computation is dominated by its final iteration. In particular, we may use the complexity bound given by Theorem 5.2.

Note that when we increase the accuracy in Algorithm 6, we also reset the precision to its minimal value. This is a matter of preference that doesn’t affect the asymptotic complexity. In practice, it seems to be preferable.

It is important to note that we do need to precompute the eigenvalues of the Gram matrix, since Algorithm 6 automatically detects the needed accuracy.

5.3. Possible generalizations

The adaptative strategy we describe for lll can be generalized to other lattice reduction algorithm. In particular, enumeration algorithms are possible within our framework, which allows the implementation of the BKZ algorithm of [25].

It would be interesting to study a generalization to sieving techniques to adapt them to approximate lattices.

6. Application to Algebraic Number Theory

We now present a direct application of our lattice reduction strategy in algorithmic number theory. Namely, we consider some interesting lattices sitting inside number fields: ideal lattices.

6.1. Number fields, integers and ideal lattices

Number fields

A number field $\mathbf{K}$ is a finite-dimensional algebraic extension of $\mathbf{Q}$. It can be described as:

[TABLE]

where $P$ is a monic irreducible polynomial of degree $d$ in $\mathbf{Z}[X]$ and where $\alpha$ denotes the image of $X$ in the quotient.

Let \mathopen{}\mathclose{{}\left(\alpha_{1},\dotsc,\alpha_{d}}\right)\in\mathbf{C}^{d} denote the distinct complex roots of $P$. Then, there are $d$ distinct ring-embeddings of $\mathbf{K}$ in $\mathbf{C}$. We define the $i$-th embedding $\sigma_{i}:\mathbf{K}\to\mathbf{C}$ as the field homomorphism sending $\alpha$ to $\alpha_{i}$.

It is classical to distinguish embeddings induced by real roots, a.k.a., real embeddings from embeddings coming from (pairs of conjugate) complex roots, called complex embeddings. Those arising from complex roots called complex embeddings.

Assume that $P$ has $r_{1}$ real roots and $r_{2}$ pairs of conjugate complex roots, with $d=r_{1}+2r_{2}$. Since the embeddings corresponding to conjugate roots are related by conjugation on $\mathbf{C}$, we can either keep a single complex root in each pair or replace each pair by the real and imaginary part of the chosen root. This leads to the Archimedean embedding $\sigma$ defined as:

[TABLE]

This embedding allows us to define a real symmetric bilinear form on $\mathbf{K}$:

[TABLE]

The second equality explains the presence of the normalization factors $\sqrt{2}$ in the definition of $\sigma$. Note that the form is positive definite, thus endowing $\mathbf{K}$ with an Euclidean structure.

Integers

Any element $\gamma$ of $\mathbf{K}$ has a minimal polynomial, defined as the unique monic polynomial of least degree among all polynomials of $\mathbf{Q}[X]$ vanishing at $\gamma$. The algebraic number $\gamma$ is said to be integral if its minimal polynomial lies in $\mathbf{Z}[X]$. The set of all integers in $\mathbf{K}$ forms a ring, called the ring of integers of $\mathbf{K}$ and denoted $\mathfrak{o}_{\mathbf{K}}$. It is also a free $\mathbf{Z}$-module of rank $d$. A basis $(w_{1},\ldots,w_{d})$ of $\mathfrak{o}_{\mathbf{K}}$ (as a $\mathbf{Z}$-module) is called an integral basis of $\mathbf{K}$.

As a consequence, using the bilinear form $\langle{\cdot},{\cdot}\rangle_{\sigma}$, we can view $\mathfrak{o}_{\mathbf{K}}$ as a lattice.

Ideals

An ideal of $\mathfrak{o}_{\mathbf{K}}$ is defined as an $\mathfrak{o}_{\mathbf{K}}$-submodule of $\mathfrak{o}_{\mathbf{K}}$. In particular, it is a $\mathbf{Z}$-submodule of rank $d$. Every ideal $I$ can be described by a two-element representation, i.e. expressed as $I=\alpha\mathfrak{o}_{\mathbf{K}}+\beta\mathfrak{o}_{\mathbf{K}},$ with $\alpha$ and $\beta$ in $\mathfrak{o}_{\mathbf{K}}$. Alternatively, every ideal can also be described by a $\mathbf{Z}$-basis formed of $d$ elements.

6.2. Lattice reduction for ideals

With the above notations, we can directly use our lattice reduction algorithm to reduce an ideal lattice. More precisely, given an integral basis $(w_{1},\ldots,w_{d})$ and a two-element representation of $I$ by $\alpha$ and $\beta$, we proceed as follows:

1. (1)

Define the Gram matrix $\mathcal{G}_{w}$ with entries $\langle{w_{i}},{w_{j}}\rangle_{\sigma}$. It can be computed to any desired precision from approximations of the roots of $P$. The roots themselves can be computed, using, for example, the Gourdon-Schönhage algorithm [8]. 2. (2)

Let $L$ be the matrix formed of the (integral) coordinates of $(\alpha w_{1},\ldots,\alpha w_{d})$ and $(\beta w_{1},\ldots,\beta w_{d})$ in the basis $(w_{1},\ldots,w_{d}).$ 3. (3)

Directly apply Algorithm 6 to $(\mathcal{G}_{w},L).$

The same thing can be done, mutantis mudantis, for an ideal described by a $\mathbf{Z}$-basis.

A well-known special case

For some number fields, the Gram matrix is $\mathcal{G}_{w}$ is integral. In that case, the use of Algorithm 6 isn’t necessary and one can directly work with an exact lattice. This is described for the special case of reducing the full lattice corresponding to the ring of integers in [1, Section 4.2] for totally real fields. It can be generalized to CM-fields, since they satisfy the same essential property of having an integral Gram matrix. The same application is also discussed in [4, Section 4.4.2].

Non integral case

For the general case where the Gram matrix is real, [1] propose to multiply by $2^{e}$ and round to the closest integer. It also gives a bound on the necessary accuracy $e$ as the logarithm of (the inverse of) the smallest diagonal entry in the Cholesky decomposition of the Gram matrix. In some sense, this is similar to our approach. However, without any auxiliary information on this coefficient, it is proposed to continue increasing $e$ as long as it is deemed unsatisfactory.

By contrast, termination of our algorithm guarantees that lattice reduction is completed and that the output basis is lll-reduced.

Appendix A Proof of Proposition 1

We now show the more general statement for a $(\delta,\eta)$-lll reduced basis $(b_{1},\ldots,b_{d})$ of $(\Lambda,\langle{\cdot},{\cdot}\rangle)$. Namely that for any $1\leq k\leq d$ we have:

[TABLE]

Proof.

Using the Lovász condition at index $1\leq i<d$, we write:

[TABLE]

Thanks to the size-reduction condition, this implies:

[TABLE]

Let $K$ denote \mathopen{}\mathclose{{}\left(\delta-\eta^{2}}\right)^{-1/2} and $\ell_{i}$ be the norm of the vector $\pi_{i}(b_{i})$. Then, Equation (4) becomes:

[TABLE]

Recall that \operatorname{covol}\mathopen{}\mathclose{{}\left({b_{1},\ldots,b_{k}}}\right)=\prod_{i=1}^{k}\ell_{i}. This implies that for any $j>k$:

[TABLE]

Thus:

[TABLE]

∎