1 | |
---|
2 | /**************************************************************************\ |
---|
3 | |
---|
4 | MODULE: LLL |
---|
5 | |
---|
6 | SUMMARY: |
---|
7 | |
---|
8 | Routines are provided for lattice basis reduction, including both |
---|
9 | exact-aritmetic variants (slow but sure) and floating-point variants |
---|
10 | (fast but only approximate). |
---|
11 | |
---|
12 | For an introduction to the basics of LLL reduction, see |
---|
13 | [H. Cohen, A Course in Computational Algebraic Number Theory, Springer, 1993]. |
---|
14 | |
---|
15 | The LLL algorithm was introduced in [A. K. Lenstra, H. W. Lenstra, and |
---|
16 | L. Lovasz, Math. Ann. 261 (1982), 515-534]. |
---|
17 | |
---|
18 | \**************************************************************************/ |
---|
19 | |
---|
20 | |
---|
21 | |
---|
22 | |
---|
23 | #include <NTL/mat_ZZ.h> |
---|
24 | |
---|
25 | |
---|
26 | |
---|
27 | /**************************************************************************\ |
---|
28 | |
---|
29 | Exact Arithmetic Variants |
---|
30 | |
---|
31 | \**************************************************************************/ |
---|
32 | |
---|
33 | |
---|
34 | |
---|
35 | |
---|
36 | long LLL(ZZ& det2, mat_ZZ& B, long verbose = 0); |
---|
37 | long LLL(ZZ& det2, mat_ZZ& B, mat_ZZ& U, long verbose = 0); |
---|
38 | |
---|
39 | long LLL(ZZ& det2, mat_ZZ& B, long a, long b, long verbose = 0); |
---|
40 | long LLL(ZZ& det2, mat_ZZ& B, mat_ZZ& U, long a, long b, long verbose = 0); |
---|
41 | |
---|
42 | |
---|
43 | // performs LLL reduction. |
---|
44 | |
---|
45 | // B is an m x n matrix, viewed as m rows of n-vectors. m may be less |
---|
46 | // than, equal to, or greater than n, and the rows need not be |
---|
47 | // linearly independent. B is transformed into an LLL-reduced basis, |
---|
48 | // and the return value is the rank r of B. The first m-r rows of B |
---|
49 | // are zero. |
---|
50 | |
---|
51 | // More specifically, elementary row transformations are performed on |
---|
52 | // B so that the non-zero rows of new-B form an LLL-reduced basis |
---|
53 | // for the lattice spanned by the rows of old-B. |
---|
54 | // The default reduction parameter is delta=3/4, which means |
---|
55 | // that the squared length of the first non-zero basis vector |
---|
56 | // is no more than 2^{r-1} times that of the shortest vector in |
---|
57 | // the lattice. |
---|
58 | |
---|
59 | // det2 is calculated as the *square* of the determinant |
---|
60 | // of the lattice---note that sqrt(det2) is in general an integer |
---|
61 | // only when r = n. |
---|
62 | |
---|
63 | // In the second version, U is set to the transformation matrix, so |
---|
64 | // that U is a unimodular m x m matrix with U * old-B = new-B. |
---|
65 | // Note that the first m-r rows of U form a basis (as a lattice) |
---|
66 | // for the kernel of old-B. |
---|
67 | |
---|
68 | // The third and fourth versions allow an arbitrary reduction |
---|
69 | // parameter delta=a/b, where 1/4 < a/b <= 1, where a and b are positive |
---|
70 | // integers. |
---|
71 | // For a basis reduced with parameter delta, the squared length |
---|
72 | // of the first non-zero basis vector is no more than |
---|
73 | // 1/(delta-1/4)^{r-1} times that of the shortest vector in the |
---|
74 | // lattice (see, e.g., the article by Schnorr and Euchner mentioned below). |
---|
75 | |
---|
76 | // The algorithm employed here is essentially the one in Cohen's book. |
---|
77 | |
---|
78 | |
---|
79 | // Some variations: |
---|
80 | |
---|
81 | long LLL_plus(vec_ZZ& D, mat_ZZ& B, long verbose = 0); |
---|
82 | long LLL_plus(vec_ZZ& D, mat_ZZ& B, mat_ZZ& U, long verbose = 0); |
---|
83 | |
---|
84 | long LLL_plus(vec_ZZ& D, mat_ZZ& B, long a, long b, long verbose = 0); |
---|
85 | long LLL_plus(vec_ZZ& D, mat_ZZ& B, mat_ZZ& U, long a, long b, |
---|
86 | long verbose = 0); |
---|
87 | |
---|
88 | // These are variations that return a bit more information about the |
---|
89 | // reduced basis. If r is the rank of B, then D is a vector of length |
---|
90 | // r+1, such that D[0] = 1, and for i = 1..r, D[i]/D[i-1] is equal to |
---|
91 | // the square of the length of the i-th vector of the Gram-Schmidt basis |
---|
92 | // corresponding to the (non-zero) rows of the LLL reduced basis B. |
---|
93 | // In particular, D[r] is equal to the value det2 computed by the |
---|
94 | // plain LLL routines. |
---|
95 | |
---|
96 | /**************************************************************************\ |
---|
97 | |
---|
98 | Computing Images and Kernels |
---|
99 | |
---|
100 | \**************************************************************************/ |
---|
101 | |
---|
102 | |
---|
103 | long image(ZZ& det2, mat_ZZ& B, long verbose = 0); |
---|
104 | long image(ZZ& det2, mat_ZZ& B, mat_ZZ& U, long verbose = 0); |
---|
105 | |
---|
106 | // This computes the image of B using a "cheap" version of the LLL: |
---|
107 | // it performs the usual "size reduction", but it only swaps |
---|
108 | // vectors when linear dependencies are found. |
---|
109 | // I haven't seen this described in the literature, but it works |
---|
110 | // fairly well in practice, and can also easily be shown |
---|
111 | // to run in a reasonable amount of time with reasonably bounded |
---|
112 | // numbers. |
---|
113 | |
---|
114 | // As in the above LLL routines, the return value is the rank r of B, and the |
---|
115 | // first m-r rows will be zero. U is a unimodular m x m matrix with |
---|
116 | // U * old-B = new-B. det2 has the same meaning as above. |
---|
117 | |
---|
118 | // Note that the first m-r rows of U form a basis (as a lattice) |
---|
119 | // for the kernel of old-B. |
---|
120 | // This is a reasonably practical algorithm for computing kernels. |
---|
121 | // One can also apply image() to the kernel to get somewhat |
---|
122 | // shorter basis vectors for the kernels (there are no linear |
---|
123 | // dependencies, but the size reduction may anyway help). |
---|
124 | // For even shorter kernel basis vectors, on can apply |
---|
125 | // LLL(). |
---|
126 | |
---|
127 | |
---|
128 | /**************************************************************************\ |
---|
129 | |
---|
130 | Finding a vector in a lattice |
---|
131 | |
---|
132 | \**************************************************************************/ |
---|
133 | |
---|
134 | long LatticeSolve(vec_ZZ& x, const mat_ZZ& A, const vec_ZZ& y, long reduce=0); |
---|
135 | |
---|
136 | // This tests if for given A and y, there exists x such that x*A = y; |
---|
137 | // if so, x is set to a solution, and the value 1 is returned; |
---|
138 | // otherwise, x is left unchanged, and the value 0 is returned. |
---|
139 | |
---|
140 | // The optional parameter reduce controls the 'quality' of the |
---|
141 | // solution vector; if the rows of A are linearly dependent, |
---|
142 | // there are many solutions, if there are any at all. |
---|
143 | // The value of reduce controls the amount of effort that goes |
---|
144 | // into finding a 'short' solution vector x. |
---|
145 | |
---|
146 | // reduce = 0: No particular effort is made to find a short solution. |
---|
147 | |
---|
148 | // reduce = 1: A simple 'size reduction' algorithm is run on kernel(A); |
---|
149 | // this is fast, and may yield somewhat shorter |
---|
150 | // solutions than the default, but not necessarily |
---|
151 | // very close at all to optimal. |
---|
152 | |
---|
153 | // reduce = 2: the LLL algorithm is run on kernel(A); |
---|
154 | // this may be significantly slower than the other options, |
---|
155 | // but yields solutions that are provably close to optimal. |
---|
156 | // More precisely, if kernel(A) has rank k, |
---|
157 | // then the squared length of the obtained solution |
---|
158 | // is no more than max(1, 2^(k-2)) times that of |
---|
159 | // the optimal solution. This makes use of slight |
---|
160 | // variation of Babai's approximately nearest vector algorithm. |
---|
161 | |
---|
162 | // Of course, if the the rows of A are linearly independent, then |
---|
163 | // the value of reduce is irrelevant: the solution, if it exists, |
---|
164 | // is unique. |
---|
165 | |
---|
166 | // Note that regardless of the value of reduce, the algorithm |
---|
167 | // runs in polynomial time, and hence the bit-length of the solution |
---|
168 | // vector is bounded by a polynomial in the bit-length of the inputs. |
---|
169 | |
---|
170 | |
---|
171 | |
---|
172 | |
---|
173 | /**************************************************************************\ |
---|
174 | |
---|
175 | Floating Point Variants |
---|
176 | |
---|
177 | There are a number of floating point LLL variants available: |
---|
178 | you can choose the precision, the orthogonalization strategy, |
---|
179 | and the reduction condition. |
---|
180 | |
---|
181 | The wide variety of choices may seem a bit bewildering. |
---|
182 | See below the discussion "How to choose?". |
---|
183 | |
---|
184 | *** Precision: |
---|
185 | |
---|
186 | FP -- double |
---|
187 | QP -- quad_float (quasi quadruple precision) |
---|
188 | this is useful when roundoff errors can cause problems |
---|
189 | XD -- xdouble (extended exponent doubles) |
---|
190 | this is useful when numbers get too big |
---|
191 | RR -- RR (arbitrary precision floating point) |
---|
192 | this is useful for large precision and magnitudes |
---|
193 | |
---|
194 | Generally speaking, the choice FP will be the fastest, |
---|
195 | but may be prone to roundoff errors and/or overflow. |
---|
196 | |
---|
197 | |
---|
198 | *** Orthogonalization Strategy: |
---|
199 | |
---|
200 | -- Classical Gramm-Schmidt Orthogonalization. |
---|
201 | This choice uses classical methods for computing |
---|
202 | the Gramm-Schmidt othogonalization. |
---|
203 | It is fast but prone to stability problems. |
---|
204 | This strategy was first proposed by Schnorr and Euchner |
---|
205 | [C. P. Schnorr and M. Euchner, Proc. Fundamentals of Computation Theory, |
---|
206 | LNCS 529, pp. 68-85, 1991]. |
---|
207 | The version implemented here is substantially different, improving |
---|
208 | both stability and performance. |
---|
209 | |
---|
210 | -- Givens Orthogonalization. |
---|
211 | This is a bit slower, but generally much more stable, |
---|
212 | and is really the preferred orthogonalization strategy. |
---|
213 | For a nice description of this, see Chapter 5 of |
---|
214 | [G. Golub and C. van Loan, Matrix Computations, 3rd edition, |
---|
215 | Johns Hopkins Univ. Press, 1996]. |
---|
216 | |
---|
217 | |
---|
218 | *** Reduction Condition: |
---|
219 | |
---|
220 | -- LLL: the classical LLL reduction condition. |
---|
221 | |
---|
222 | -- BKZ: Block Korkin-Zolotarev reduction. |
---|
223 | This is slower, but yields a higher-quality basis, |
---|
224 | i.e., one with shorter vectors. |
---|
225 | See the Schnorr-Euchner paper for a description of this. |
---|
226 | This basically generalizes the LLL reduction condition |
---|
227 | from blocks of size 2 to blocks of larger size. |
---|
228 | |
---|
229 | |
---|
230 | ************* Calling Syntax for LLL routines *************** |
---|
231 | |
---|
232 | long |
---|
233 | [G_]LLL_{FP,QP,XD,RR} (mat_ZZ& B, [ mat_ZZ& U, ] double delta = 0.99, |
---|
234 | long deep = 0, LLLCheckFct check = 0, long verbose = 0); |
---|
235 | |
---|
236 | * The [ ... ] notation indicates something optional, |
---|
237 | and the { ... } indicates something that is chosen from |
---|
238 | among several alternatives. |
---|
239 | |
---|
240 | * The return value is the rank of B (but see below if check != 0). |
---|
241 | |
---|
242 | * The optional prefix G_ indicates that Givens rotations are to be used; |
---|
243 | otherwise, classical Gramm-Schmidt is used. |
---|
244 | |
---|
245 | * The choice FP, QP, XD, RR determines the precision used. |
---|
246 | |
---|
247 | * If the optional parameter U is given, then U is computed |
---|
248 | as the transition matrix: |
---|
249 | |
---|
250 | U * old_B = new_B |
---|
251 | |
---|
252 | * The optional argument "delta" is the reduction parameter, and may |
---|
253 | be set so that 0.50 <= delta < 1. Setting it close to 1 yields |
---|
254 | shorter vectors, and also improves the stability, but increases the |
---|
255 | running time. Recommended value: delta = 0.99. |
---|
256 | |
---|
257 | * The optional parameter "deep" can be set to any positive integer, |
---|
258 | which allows "deep insertions" of row k into row i, provided i <= |
---|
259 | deep or k-i <= deep. Larger values of deep will usually yield |
---|
260 | shorter vectors, but the running increases exponentially. |
---|
261 | |
---|
262 | NOTE: use of "deep" is obsolete, and has been "deprecated". |
---|
263 | It is recommended to use BKZ_FP to achieve higher-quality reductions. |
---|
264 | Moreover, the Givens versions do not support "deep", and setting |
---|
265 | deep != 0 will raise an error in this case. |
---|
266 | |
---|
267 | * The optional parameter "check" is a function that is invoked after |
---|
268 | each size reduction with the current row as an argument. If this |
---|
269 | function returns a non-zero value, the LLL procedure is immediately |
---|
270 | terminated. Note that it is possible that some linear dependencies |
---|
271 | remain undiscovered, so that the calculated rank value is in fact |
---|
272 | too large. In any case, zero rows discovered by the algorithm |
---|
273 | will be placed at the beginning, as usual. |
---|
274 | |
---|
275 | The check argument (if not zero) should be a routine taking |
---|
276 | a const vec_ZZ& as an argument and return value of type long. |
---|
277 | LLLCheckFct is defined via a typedef as: |
---|
278 | |
---|
279 | typedef long (*LLLCheckFct)(const vec_ZZ&); |
---|
280 | |
---|
281 | See the file subset.c for an example of the use of this feature. |
---|
282 | |
---|
283 | * The optional parameter "verbose" can be set to see all kinds of fun |
---|
284 | things printed while the routine is executing. A status report is |
---|
285 | printed every once in a while, and the current basis is optionally |
---|
286 | dumped to a file. The behavior can be controlled with these global |
---|
287 | variables: |
---|
288 | |
---|
289 | extern char *LLLDumpFile; // file to dump basis, 0 => no dump; |
---|
290 | // initially 0 |
---|
291 | |
---|
292 | extern double LLLStatusInterval; // seconds between status reports |
---|
293 | // initially 900s = 15min |
---|
294 | |
---|
295 | |
---|
296 | |
---|
297 | |
---|
298 | ************* Calling Syntax for BKZ routines *************** |
---|
299 | |
---|
300 | long |
---|
301 | [G_]BKZ_{FP,QP,QP1,XD,RR} (mat_ZZ& B, [ mat_ZZ& U, ] double delta=0.99, |
---|
302 | long BlockSize=10, long prune=0, |
---|
303 | LLLCheckFct check = 0, long verbose = 0); |
---|
304 | |
---|
305 | These functions are equivalent to the LLL routines above, |
---|
306 | except that Block Korkin-Zolotarev reduction is applied. |
---|
307 | We describe here only the differences in the calling syntax. |
---|
308 | |
---|
309 | * The optional parameter "BlockSize" specifies the size of the blocks |
---|
310 | in the reduction. High values yield shorter vectors, but the |
---|
311 | running time increases exponentially with BlockSize. |
---|
312 | BlockSize should be between 2 and the number of rows of B. |
---|
313 | |
---|
314 | * The optional parameter "prune" can be set to any positive number to |
---|
315 | invoke the Volume Heuristic from [Schnorr and Horner, Eurocrypt |
---|
316 | '95]. This can significantly reduce the running time, and hence |
---|
317 | allow much bigger block size, but the quality of the reduction is |
---|
318 | of course not as good in general. Higher values of prune mean |
---|
319 | better quality, and slower running time. |
---|
320 | When prune == 0, pruning is disabled. |
---|
321 | Recommended usage: for BlockSize >= 30, set 10 <= prune <= 15. |
---|
322 | |
---|
323 | * The QP1 variant uses quad_float precision to compute Gramm-Schmidt, |
---|
324 | but uses double precision in the search phase of the block reduction |
---|
325 | algorithm. This seems adequate for most purposes, and is faster |
---|
326 | than QP, which uses quad_float precision uniformly throughout. |
---|
327 | |
---|
328 | |
---|
329 | ******************** How to choose? ********************* |
---|
330 | |
---|
331 | I think it is safe to say that nobody really understands |
---|
332 | how the LLL algorithm works. The theoretical analyses are a long way |
---|
333 | from describing what "really" happens in practice. Choosing the best |
---|
334 | variant for a certain application ultimately is a matter of trial |
---|
335 | and error. |
---|
336 | |
---|
337 | The first thing to try is LLL_FP. |
---|
338 | It is the fastest of the routines, and is adequate for many applications. |
---|
339 | |
---|
340 | If there are precision problems, you will most likely get |
---|
341 | a warning message, something like "warning--relaxing reduction". |
---|
342 | If there are overflow problems, you should get an error message |
---|
343 | saying that the numbers are too big. |
---|
344 | |
---|
345 | If either of these happens, the next thing to try is G_LLL_FP, |
---|
346 | which uses the somewhat slower, but more stable, Givens rotations. |
---|
347 | This approach also has the nice property that the numbers remain |
---|
348 | smaller, so there is less chance of an overflow. |
---|
349 | |
---|
350 | If you are still having precision problems with G_LLL_FP, |
---|
351 | try LLL_QP or G_LLL_QP, which uses quadratic precision. |
---|
352 | |
---|
353 | If you are still having overflow problems, try LLL_XD or G_LLL_XD. |
---|
354 | |
---|
355 | I haven't yet come across a case where one *really* needs the |
---|
356 | extra precision available in the RR variants. |
---|
357 | |
---|
358 | All of the above discussion applies to the BKZ variants as well. |
---|
359 | In addition, if you have a matrix with really big entries, you might try |
---|
360 | using G_LLL_FP or LLL_XD first to reduce the sizes of the numbers, |
---|
361 | before running one of the BKZ variants. |
---|
362 | |
---|
363 | Also, one shouldn't rule out using the "all integer" LLL routines. |
---|
364 | For some highly structured matrices, this is not necessarily |
---|
365 | much worse than some of the floating point versions, and can |
---|
366 | under certain circumstances even be better. |
---|
367 | |
---|
368 | |
---|
369 | ******************** Implementation notes ********************* |
---|
370 | |
---|
371 | For all the floating point variants, I use a "relaxed" size reduction |
---|
372 | condition. Normally in LLL one makes all |\mu_{i,j}| <= 1/2. |
---|
373 | However, this can easily lead to infinite loops in floating point arithemetic. |
---|
374 | So I use the condition |\mu_{i,j}| <= 1/2 + fudge, where fudge is |
---|
375 | a very small number. Even with this, one can fall into an infinite loop. |
---|
376 | To handle this situation, I added some logic that detects, at quite low cost, |
---|
377 | when an infinite loop has been entered. When that happens, fudge |
---|
378 | is replaced by fudge*2, and a warning message "relaxing reduction condition" |
---|
379 | is printed. We may do this relaxation several times. |
---|
380 | If fudge gets too big, we give up and abort, except that |
---|
381 | LLL_FP and BKZ_FP make one last attempt to recover: they try to compute the |
---|
382 | Gramm-Schmidt coefficients using RR and continue. As described above, |
---|
383 | if you run into these problems, which you'll see in the error/warning |
---|
384 | messages, it is more effective to use the QP and/or Givens variants. |
---|
385 | |
---|
386 | For the Gramm-Schmidt orthogonalization, lots of "bookeeping" is done |
---|
387 | to avoid computing the same thing twice. |
---|
388 | |
---|
389 | For the Givens orthogonalization, we cannot do so many bookeeping tricks. |
---|
390 | Instead, we "cache" a certain amount of information, which |
---|
391 | allows us to avoid computing certain things over and over again. |
---|
392 | |
---|
393 | There are many other hacks and tricks to speed things up even further. |
---|
394 | For example, if the matrix elements are small enough to fit in |
---|
395 | double precision floating point, the algorithms avoid almost |
---|
396 | all big integer arithmetic. This is done in a dynamic, on-line |
---|
397 | fashion, so even if the numbers start out big, whenever they |
---|
398 | get small, we automatically switch to floating point arithmetic. |
---|
399 | |
---|
400 | \**************************************************************************/ |
---|
401 | |
---|
402 | |
---|
403 | |
---|
404 | |
---|
405 | /**************************************************************************\ |
---|
406 | |
---|
407 | Other Stuff |
---|
408 | |
---|
409 | \**************************************************************************/ |
---|
410 | |
---|
411 | |
---|
412 | |
---|
413 | void ComputeGS(const mat_ZZ& B, mat_RR& mu, vec_RR& c); |
---|
414 | |
---|
415 | // Computes Gramm-Schmidt data for B. Assumes B is an m x n matrix of |
---|
416 | // rank m. Let if { B^*(i) } is the othogonal basis, then c(i) = |
---|
417 | // |B^*(i)|^2, and B^*(i) = B(i) - \sum_{j=1}^{i-1} mu(i,j) B^*(j). |
---|
418 | |
---|
419 | void NearVector(vec_ZZ& w, const mat_ZZ& B, const vec_ZZ& a); |
---|
420 | |
---|
421 | // Computes a vector w that is an approximation to the closest vector |
---|
422 | // in the lattice spanned by B to a, using the "closest plane" |
---|
423 | // algorithm from [Babai, Combinatorica 6:1-13, 1986]. B must be a |
---|
424 | // square matrix, and it is assumed that B is already LLL or BKZ |
---|
425 | // reduced (the better the reduction the better the approximation). |
---|
426 | // Note that arithmetic in RR is used with the current value of |
---|
427 | // RR::precision(). |
---|
428 | |
---|
429 | // NOTE: Both of these routines use classical Gramm-Schmidt |
---|
430 | // orthogonalization. |
---|
431 | |
---|
432 | |
---|