Line data Source code
1 : /* Copyright (C) 2008 The PARI group.
2 :
3 : This file is part of the PARI/GP package.
4 :
5 : PARI/GP is free software; you can redistribute it and/or modify it under the
6 : terms of the GNU General Public License as published by the Free Software
7 : Foundation; either version 2 of the License, or (at your option) any later
8 : version. It is distributed in the hope that it will be useful, but WITHOUT
9 : ANY WARRANTY WHATSOEVER.
10 :
11 : Check the License for details. You should have received a copy of it, along
12 : with the package; see the file 'COPYING'. If not, write to the Free Software
13 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14 :
15 : #include "pari.h"
16 : #include "paripriv.h"
17 :
18 : #define DEBUGLEVEL DEBUGLEVEL_qflll
19 :
20 : static int
21 45828 : RgM_is_square_mat(GEN x) { long l = lg(x); return l == 1 || l == lgcols(x); }
22 :
23 : static long
24 4178489 : ZM_is_upper(GEN R)
25 : {
26 4178489 : long i,j, l = lg(R);
27 4178489 : if (l != lgcols(R)) return 0;
28 8072300 : for(i = 1; i < l; i++)
29 8717051 : for(j = 1; j < i; j++)
30 4482388 : if (signe(gcoeff(R,i,j))) return 0;
31 252061 : return 1;
32 : }
33 :
34 : static long
35 606085 : ZM_is_knapsack(GEN R)
36 : {
37 606085 : long i,j, l = lg(R);
38 606085 : if (l != lgcols(R)) return 0;
39 843323 : for(i = 2; i < l; i++)
40 2900702 : for(j = 1; j < l; j++)
41 2663464 : if ( i!=j && signe(gcoeff(R,i,j))) return 0;
42 92365 : return 1;
43 : }
44 :
45 : static long
46 1182517 : ZM_is_lower(GEN R)
47 : {
48 1182517 : long i,j, l = lg(R);
49 1182517 : if (l != lgcols(R)) return 0;
50 2058906 : for(i = 1; i < l; i++)
51 2383985 : for(j = 1; j < i; j++)
52 1291920 : if (signe(gcoeff(R,j,i))) return 0;
53 34904 : return 1;
54 : }
55 :
56 : static GEN
57 34904 : RgM_flip(GEN R)
58 : {
59 : GEN M;
60 : long i,j,l;
61 34904 : M = cgetg_copy(R, &l);
62 181322 : for(i = 1; i < l; i++)
63 : {
64 146418 : gel(M,i) = cgetg(l, t_COL);
65 910402 : for(j = 1; j < l; j++)
66 763984 : gmael(M,i,j) = gmael(R,l-i, l-j);
67 : }
68 34904 : return M;
69 : }
70 :
71 : static GEN
72 0 : RgM_flop(GEN R)
73 : {
74 : GEN M;
75 : long i,j,l;
76 0 : M = cgetg_copy(R, &l);
77 0 : for(i = 1; i < l; i++)
78 : {
79 0 : gel(M,i) = cgetg(l, t_COL);
80 0 : for(j = 1; j < l; j++)
81 0 : gmael(M,i,j) = gmael(R,i, l-j);
82 : }
83 0 : return M;
84 : }
85 :
86 : /* Assume x and y has same type! */
87 : INLINE int
88 4063890 : mpabscmp(GEN x, GEN y)
89 : {
90 4063890 : return (typ(x)==t_INT) ? abscmpii(x,y) : abscmprr(x,y);
91 : }
92 :
93 : /****************************************************************************/
94 : /*** FLATTER ***/
95 : /****************************************************************************/
96 : /* Implementation of "FLATTER" algorithm based on
97 : * <https://eprint.iacr.org/2023/237>
98 : * Fast Practical Lattice Reduction through Iterated Compression
99 : *
100 : * Keegan Ryan, University of California, San Diego
101 : * Nadia Heninger, University of California, San Diego. BA20230925 */
102 : static long
103 1335750 : drop(GEN R)
104 : {
105 1335750 : long i, n = lg(R)-1;
106 1335750 : long s = 0, m = mpexpo(gcoeff(R, 1, 1));
107 5399640 : for (i = 2; i <= n; ++i)
108 : {
109 4063890 : if (mpabscmp(gcoeff(R, i, i), gcoeff(R, i - 1, i - 1)) >= 0)
110 : {
111 2753404 : s += m - mpexpo(gcoeff(R, i - 1, i - 1));
112 2753404 : m = mpexpo(gcoeff(R, i, i));
113 : }
114 : }
115 1335750 : s += m - mpexpo(gcoeff(R, n, n));
116 1335749 : return s;
117 : }
118 :
119 : static long
120 1335750 : potential(GEN R)
121 : {
122 1335750 : long i, n = lg(R)-1;
123 1335750 : long s = 0, mul = n-1;;
124 6735389 : for (i = 1; i <= n; i++, mul-=2) s += mul * mpexpo(gcoeff(R,i,i));
125 1335750 : return s;
126 : }
127 :
128 : /* U upper-triangular invertible:
129 : * Bound on the exponent of the condition number of U.
130 : * Algo 8.13 in Higham, Accuracy and stability of numercal algorithms. */
131 : static long
132 4689745 : condition_bound(GEN U, int lower)
133 : {
134 4689745 : long n = lg(U)-1, e, i, j;
135 : GEN y;
136 4689745 : pari_sp av = avma;
137 4689745 : y = cgetg(n+1, t_VECSMALL);
138 4689746 : e = y[n] = -gexpo(gcoeff(U,n,n));
139 18641823 : for (i=n-1; i>0; i--)
140 : {
141 13952077 : long s = 0;
142 49602021 : for (j=i+1; j<=n; j++)
143 35649943 : s = maxss(s, (lower? gexpo(gcoeff(U,j,i)): gexpo(gcoeff(U,i,j))) + y[j]);
144 13952078 : y[i] = s - gexpo(gcoeff(U,i,i));
145 13952080 : e = maxss(e, y[i]);
146 : }
147 4689746 : return gc_long(av, gexpo(U) + e);
148 : }
149 :
150 : static long
151 5149784 : gsisinv(GEN M)
152 : {
153 5149784 : long i, l = lg(M);
154 25907880 : for (i = 1; i < l; ++i)
155 20758504 : if (! signe(gmael(M, i, i))) return 0;
156 5149376 : return 1;
157 : }
158 :
159 : INLINE long
160 7415245 : nbits2prec64(long n)
161 : {
162 7415245 : return nbits2prec(((n+63)>>6)<<6);
163 : }
164 :
165 : static long
166 5806484 : spread(GEN R)
167 : {
168 5806484 : long i, n = lg(R)-1, m = mpexpo(gcoeff(R, 1, 1)), M = m;
169 23364005 : for (i = 2; i <= n; ++i)
170 : {
171 17557519 : long e = mpexpo(gcoeff(R, i, i));
172 17557518 : if (e < m) m = e;
173 17557518 : if (e > M) M = e;
174 : }
175 5806486 : return M - m;
176 : }
177 :
178 : static long
179 4689744 : GS_extraprec(GEN L, int lower)
180 : {
181 4689744 : long C = condition_bound(L, lower), S = spread(L), n = lg(L)-1;
182 4689746 : return maxss(2*S+2*n, C-S-2*n); /* = 2*S + 2*n + maxss(0, C-3*S-4*n) */
183 : }
184 :
185 : static GEN
186 2967 : RgM_Cholesky_dynprec(GEN M)
187 : {
188 2967 : pari_sp ltop = avma;
189 : GEN L;
190 2967 : long minprec = lg(M) + 30, bitprec = minprec, prec;
191 : while (1)
192 4877 : {
193 : long mbitprec;
194 7844 : prec = nbits2prec64(bitprec);
195 7844 : L = RgM_Cholesky(RgM_gtofp(M, prec), prec); /* upper-triangular */
196 7844 : if (!L)
197 : {
198 1458 : bitprec *= 2;
199 1458 : set_avma(ltop);
200 1458 : continue;
201 : }
202 6386 : mbitprec = minprec + GS_extraprec(L, 0);
203 6386 : if (bitprec >= mbitprec)
204 2967 : break;
205 3419 : bitprec = maxss((4*bitprec)/3, mbitprec);
206 3419 : set_avma(ltop);
207 : }
208 2967 : return gerepilecopy(ltop, L);
209 : }
210 :
211 : static GEN
212 1336 : gramschmidt_upper(GEN M)
213 : {
214 1336 : long bitprec = lg(M)-1 + 31 + GS_extraprec(M, 0);
215 1336 : return RgM_gtofp(M, nbits2prec64(bitprec));
216 : }
217 :
218 : static GEN
219 2671497 : gramschmidt_dynprec(GEN M)
220 : {
221 2671497 : pari_sp ltop = avma;
222 2671497 : long minprec = lg(M) + 30, bitprec = minprec;
223 2671497 : if (ZM_is_upper(M)) return gramschmidt_upper(M);
224 : while (1)
225 3608992 : {
226 : GEN B, Q, L;
227 6279153 : long prec = nbits2prec64(bitprec), mbitprec;
228 6279150 : if (!QR_init(RgM_gtofp(M, prec), &B, &Q, &L, prec) || !gsisinv(L))
229 : {
230 1597125 : bitprec *= 2;
231 1597125 : set_avma(ltop);
232 1597128 : continue;
233 : }
234 4682020 : mbitprec = minprec + GS_extraprec(L, 1);
235 4682024 : if (bitprec >= mbitprec)
236 2670161 : return gerepilecopy(ltop, shallowtrans(L));
237 2011863 : bitprec = maxss((4*bitprec)/3, mbitprec);
238 2011864 : set_avma(ltop);
239 : }
240 : }
241 : /* return -T1 * round(T1^-1*(R1^-1*R2)*T3) */
242 : static GEN
243 1335750 : sizered(GEN T1, GEN T3, GEN R1, GEN R2)
244 : {
245 1335750 : pari_sp ltop = avma;
246 : long e;
247 1335750 : return gerepileupto(ltop, ZM_mul(ZM_neg(T1), grndtoi(gmul(ZM_inv(T1,NULL),
248 : RgM_mul(RgM_mul(RgM_inv_upper(R1), R2), T3)), &e)));
249 : }
250 :
251 : static GEN
252 1335750 : flat(GEN M, long flag, GEN *pt_T, long *pt_s, long *pt_pot)
253 : {
254 1335750 : pari_sp ltop = avma;
255 : GEN R, R1, R2, R3, T1, T2, T3, T, S;
256 1335750 : long k = lg(M)-1, n = k>>1, n2 = k - n, m = n>>1;
257 1335750 : long keepfirst = flag & LLL_KEEP_FIRST, inplace = flag & LLL_INPLACE;
258 : /* for k = 3, we want n = 1; n2 = 2; m = 0 */
259 : /* for k = 5, n = 2; n2 = 3; m = 1 */
260 1335750 : R = gramschmidt_dynprec(M);
261 1335750 : R1 = matslice(R, 1, n, 1, n);
262 1335749 : R2 = matslice(R, 1, n, n + 1, k);
263 1335749 : R3 = matslice(R, n + 1, k, n + 1, k);
264 1335750 : T1 = lllfp(R1, 0.99, LLL_IM| LLL_UPPER| LLL_NOCERTIFY| (keepfirst ? LLL_KEEP_FIRST: 0));
265 1335749 : T3 = lllfp(R3, 0.99, LLL_IM| LLL_UPPER| LLL_NOCERTIFY);
266 1335750 : T2 = sizered(T1, T3, R1, R2);
267 1335750 : T = shallowmatconcat(mkmat22(T1,T2,gen_0,T3));
268 1335750 : M = ZM_mul(M, T);
269 1335748 : R = gramschmidt_dynprec(M);
270 1335750 : R3 = matslice(R, m + 1, m + n2, m + 1, m + n2);
271 1335750 : T3 = lllfp(R3, 0.99, LLL_IM| LLL_UPPER| LLL_NOCERTIFY);
272 2671500 : S = shallowmatconcat(diagonal(
273 576070 : m == 0 ? mkvec2(T3, matid(k - m - n2))
274 0 : : m+n2 == k ? mkvec2(matid(m), T3)
275 759680 : : mkvec3(matid(m), T3, matid(k - m - n2))));
276 1335750 : M = ZM_mul(M, S);
277 1335747 : if (!inplace) *pt_T = ZM_mul(T, S);
278 1335749 : *pt_s = drop(R);
279 1335750 : *pt_pot = potential(R);
280 1335749 : return gc_all(ltop, inplace ? 1: 2, &M, pt_T);
281 : }
282 :
283 : static void
284 0 : dbg_flatter(pari_timer *ti, long n, long i, long lti, double t, double pot2)
285 : {
286 0 : double s = t / n, p = pot2 / (n*(n+1));
287 : const char *str;
288 0 : if (i == -1)
289 0 : str = (i == lti)? "final"
290 0 : : stack_sprintf("steps %ld-final", lti);
291 : else
292 0 : str = (i == lti)? stack_sprintf("step %ld", i)
293 0 : : stack_sprintf("steps %ld-%ld", lti, i);
294 0 : timer_printf(ti, "FLATTER, dim %ld, %s: \t slope=%0.10g \t pot=%0.10g",
295 : n, str, s, p);
296 0 : }
297 :
298 : static GEN
299 618384 : ZM_flatter(GEN M, long flag)
300 : {
301 618384 : pari_sp av = avma;
302 618384 : long i, n = lg(M)-1, s = -1, lti = 1, pot = LONG_MAX;
303 618384 : GEN T = NULL;
304 : pari_timer ti;
305 618384 : long inplace = flag & LLL_INPLACE, cert = !(flag & LLL_NOCERTIFY);
306 :
307 618384 : if (DEBUGLEVEL>=3)
308 : {
309 0 : timer_start(&ti);
310 0 : if (cert) err_printf("FLATTER dim = %ld size = %ld\n", n, ZM_max_expi(M));
311 : }
312 618386 : for (i = 1;;i++)
313 717364 : {
314 : long t, pot2;
315 1335750 : GEN U, M2 = flat(M, flag, &U, &t, &pot2);
316 1335749 : if (t == 0) { s = t; break; }
317 760643 : if (s >= 0)
318 : {
319 435561 : if (s == t && pot>=pot2) break;
320 392282 : if (s < t && i > 20)
321 : {
322 0 : if (DEBUGLEVEL >= 3) err_printf("BACK:%ld:%ld:%g\n", n, i, s);
323 0 : break;
324 : }
325 : }
326 717364 : if (DEBUGLEVEL>=3 && (cert || timer_get(&ti) > 1000))
327 0 : dbg_flatter(&ti, n, i, lti, t, pot2);
328 717364 : s = t;
329 717364 : pot = pot2;
330 717364 : M = M2;
331 717364 : if (!inplace)
332 : {
333 689732 : T = T? ZM_mul(T, U): U;
334 689732 : if (gc_needed(av, 1)) gerepileall(av, 2, &M, &T);
335 : }
336 : else
337 27632 : if (gc_needed(av, 1)) M = gerepilecopy(av, M);
338 : }
339 618385 : if (DEBUGLEVEL>=3 && (cert || timer_get(&ti) > 1000))
340 0 : dbg_flatter(&ti, n, -1, i == lti? -1: lti, s, pot);
341 618385 : if (!inplace)
342 : {
343 604377 : if (!T) return gc_NULL(av);
344 311242 : return gerepilecopy(av, T);
345 : }
346 14008 : return gerepilecopy(av, M);
347 : }
348 :
349 : static GEN
350 616370 : ZM_flatter_rank(GEN M, long rank, long flag)
351 : {
352 : pari_timer ti;
353 616370 : pari_sp av = avma;
354 616370 : GEN T = NULL;
355 616370 : long i, n = lg(M)-1, sm = LONG_MAX;
356 616370 : long inplace = flag & LLL_INPLACE;
357 :
358 616370 : if (rank == n) return ZM_flatter(M, flag);
359 3785 : if (DEBUGLEVEL>=3) timer_start(&ti);
360 3785 : for (i = 1;; i++)
361 2014 : {
362 5799 : GEN S = ZM_flatter(vconcat(gshift(M,i),matid(n)), flag);
363 : long s;
364 5799 : if (!S || (s = expi(gnorml2(S))) >= sm) break;
365 2014 : sm = s;
366 2014 : if (DEBUGLEVEL>=3) timer_printf(&ti,"FLATTERRANK step %ld: %ld",i,sm);
367 2014 : T = T? ZM_mul(T, S): S;
368 2014 : M = ZM_mul(M, S);
369 2014 : if (gc_needed(av, 1)) gerepileall(av, 2, &M, &T);
370 : }
371 3785 : if (!inplace)
372 : {
373 3778 : if (!T) { set_avma(av); return matid(n); }
374 1951 : return gerepilecopy(av, T);
375 : }
376 7 : return gerepilecopy(av, M);
377 : }
378 :
379 : static GEN
380 2967 : flattergram_i(GEN M, long flag)
381 : {
382 2967 : pari_sp av = avma;
383 2967 : GEN T, R = RgM_Cholesky_dynprec(M);
384 2967 : T = lllfp(R, 0.99, LLL_IM|LLL_UPPER|LLL_NOCERTIFY | (flag&LLL_KEEP_FIRST));
385 2967 : return gerepileupto(av, T);
386 : }
387 :
388 : static void
389 0 : dbg_flattergram(pari_timer *t, long n, long i, long s)
390 0 : { timer_printf(t, "FLATTERGRAM, dim %ld step %ld, slope=%0.10g", n, i,
391 0 : ((double)s)/n); }
392 : /* return base change, NULL if identity */
393 : static GEN
394 961 : ZM_flattergram(GEN M, long flag)
395 : {
396 961 : pari_sp av = avma;
397 961 : GEN T = NULL;
398 961 : long i, n = lg(M)-1, s = -1;
399 :
400 : pari_timer ti;
401 961 : if (DEBUGLEVEL>=3)
402 : {
403 0 : timer_start(&ti);
404 0 : err_printf("FLATTERGRAM dim = %ld size = %ld\n", n, ZM_max_expi(M));
405 : }
406 961 : for (i = 1;; i++)
407 2006 : {
408 2967 : GEN S = flattergram_i(M, flag);
409 2967 : long t = expi(gnorml2(S));
410 2967 : if (t == 0) { s = t; break; }
411 2967 : if (s)
412 : {
413 2967 : double st = s - t;
414 2967 : if (st == 0) break;
415 2006 : if (st < 0 && i > 20)
416 : {
417 0 : if (DEBUGLEVEL >= 3)
418 0 : err_printf("BACK:%ld:%ld:%0.10g\n", n, i, ((double)s)/n);
419 0 : break;
420 : }
421 : }
422 2006 : T = T? ZM_mul(T, S): S;
423 2006 : M = qf_ZM_apply(M, S);
424 2006 : s = t;
425 2006 : if (DEBUGLEVEL >= 3) dbg_flattergram(&ti, n, i, s);
426 2006 : if (gc_needed(av, 1)) gerepileall(av, 2, &M, &T);
427 : }
428 961 : if (DEBUGLEVEL >= 3) dbg_flattergram(&ti, n, i, s);
429 961 : if (!T && ZM_isidentity(T)) return gc_NULL(av);
430 961 : return gerepilecopy(av, T);
431 : }
432 :
433 : /* return base change, NULL if identity */
434 : static GEN
435 961 : ZM_flattergram_rank(GEN M, long rank, long flag)
436 : {
437 : pari_timer ti;
438 961 : pari_sp av = avma;
439 961 : GEN T = NULL;
440 961 : long i, n = lg(M)-1;
441 961 : if (rank == n) return ZM_flattergram(M, flag);
442 0 : if (DEBUGLEVEL>=3) timer_start(&ti);
443 0 : for (i = 1;; i++)
444 0 : {
445 0 : GEN S = ZM_flattergram(RgM_Rg_add(gshift(M, i), gen_1), flag);
446 0 : if (DEBUGLEVEL>=3)
447 0 : timer_printf(&ti,"FLATTERGRAMRANK step %ld: %ld",i,expi(gnorml2(S)));
448 0 : if (!S) break;
449 0 : T = T? ZM_mul(T, S): S;
450 0 : M = qf_ZM_apply(M, S);
451 0 : if (gc_needed(av, 1)) gerepileall(av, 2, &M, &T);
452 : }
453 0 : if (!T || ZM_isidentity(T)) return gc_NULL(av);
454 0 : return gerepilecopy(av, T);
455 : }
456 :
457 : /* round to closest integer (as a double). If |a| >= 2^52, return it */
458 : static double
459 10719455 : pari_rint(double a)
460 : {
461 : #ifdef HAS_RINT
462 10719455 : return rint(a);
463 : #else
464 : const double pow2 = 4.5035996273704960e+15; /* 2^52 */
465 : double r, fa = fabs(a);
466 : if (fa >= pow2) return a;
467 : r = (pow2 + fa) - pow2;
468 : if (a < 0) r = -r;
469 : return r;
470 : #endif
471 : }
472 :
473 : /* default quality ratio for LLL */
474 : static const double LLLDFT = 0.99;
475 :
476 : /* assume flag & (LLL_KER|LLL_IM|LLL_ALL). LLL_INPLACE implies LLL_IM */
477 : static GEN
478 769926 : lll_trivial(GEN x, long flag)
479 : {
480 769926 : if (lg(x) == 1)
481 : { /* dim x = 0 */
482 15498 : if (! (flag & LLL_ALL)) return cgetg(1,t_MAT);
483 28 : retmkvec2(cgetg(1,t_MAT), cgetg(1,t_MAT));
484 : }
485 : /* dim x = 1 */
486 754428 : if (gequal0(gel(x,1)))
487 : {
488 145 : if (flag & LLL_KER) return matid(1);
489 145 : if (flag & (LLL_IM|LLL_INPLACE)) return cgetg(1,t_MAT);
490 28 : retmkvec2(matid(1), cgetg(1,t_MAT));
491 : }
492 754277 : if (flag & LLL_INPLACE) return gcopy(x);
493 650893 : if (flag & LLL_KER) return cgetg(1,t_MAT);
494 650893 : if (flag & LLL_IM) return matid(1);
495 28 : retmkvec2(cgetg(1,t_MAT), (flag & LLL_GRAM)? gcopy(x): matid(1));
496 : }
497 :
498 : /* vecslice(x,#x-k,#x) in place. Works for t_MAT, t_VEC/t_COL */
499 : static GEN
500 2055005 : vectail_inplace(GEN x, long k)
501 : {
502 2055005 : if (!k) return x;
503 57766 : x[k] = ((ulong)x[0] & ~LGBITS) | _evallg(lg(x) - k);
504 57766 : return x + k;
505 : }
506 :
507 : /* k = dim Kernel */
508 : static GEN
509 2128994 : lll_finish(GEN h, long k, long flag)
510 : {
511 : GEN g;
512 2128994 : if (!(flag & (LLL_IM|LLL_KER|LLL_ALL|LLL_INPLACE))) return h;
513 2055032 : if (flag & (LLL_IM|LLL_INPLACE)) return vectail_inplace(h, k);
514 97 : if (flag & LLL_KER) { setlg(h,k+1); return h; }
515 69 : g = vecslice(h,1,k); /* done first: vectail_inplace kills h */
516 70 : return mkvec2(g, vectail_inplace(h, k));
517 : }
518 :
519 : /* y * z * 2^e, e >= 0; y,z t_INT */
520 : INLINE GEN
521 318718 : mulshift(GEN y, GEN z, long e)
522 : {
523 318718 : long ly = lgefint(y), lz;
524 : pari_sp av;
525 : GEN t;
526 318718 : if (ly == 2) return gen_0;
527 46676 : lz = lgefint(z);
528 46676 : av = avma; (void)new_chunk(ly+lz+nbits2lg(e)); /* HACK */
529 46676 : t = mulii(z, y);
530 46676 : set_avma(av); return shifti(t, e);
531 : }
532 :
533 : /* x - y * z * 2^e, e >= 0; x,y,z t_INT */
534 : INLINE GEN
535 1508442 : submulshift(GEN x, GEN y, GEN z, long e)
536 : {
537 1508442 : long lx = lgefint(x), ly, lz;
538 : pari_sp av;
539 : GEN t;
540 1508442 : if (!e) return submulii(x, y, z);
541 1498813 : if (lx == 2) { t = mulshift(y, z, e); togglesign(t); return t; }
542 1199097 : ly = lgefint(y);
543 1199097 : if (ly == 2) return icopy(x);
544 960955 : lz = lgefint(z);
545 960955 : av = avma; (void)new_chunk(lx+ly+lz+nbits2lg(e)); /* HACK */
546 960955 : t = shifti(mulii(z, y), e);
547 960955 : set_avma(av); return subii(x, t);
548 : }
549 : static void
550 18515846 : subzi(GEN *a, GEN b)
551 : {
552 18515846 : pari_sp av = avma;
553 18515846 : b = subii(*a, b);
554 18515862 : if (lgefint(b)<=lg(*a) && isonstack(*a)) { affii(b,*a); set_avma(av); }
555 2104159 : else *a = b;
556 18515911 : }
557 :
558 : static void
559 17772248 : addzi(GEN *a, GEN b)
560 : {
561 17772248 : pari_sp av = avma;
562 17772248 : b = addii(*a, b);
563 17772215 : if (lgefint(b)<=lg(*a) && isonstack(*a)) { affii(b,*a); set_avma(av); }
564 1876518 : else *a = b;
565 17772315 : }
566 :
567 : /* x - u*y * 2^e */
568 : INLINE GEN
569 4168395 : submuliu2n(GEN x, GEN y, ulong u, long e)
570 : {
571 : pari_sp av;
572 4168395 : long ly = lgefint(y);
573 4168395 : if (ly == 2) return x;
574 2858119 : av = avma;
575 2858119 : (void)new_chunk(3+ly+lgefint(x)+nbits2lg(e)); /* HACK */
576 2860188 : y = shifti(mului(u,y), e);
577 2858964 : set_avma(av); return subii(x, y);
578 : }
579 : /* *x -= u*y * 2^e */
580 : INLINE void
581 262957 : submulzu2n(GEN *x, GEN y, ulong u, long e)
582 : {
583 : pari_sp av;
584 262957 : long ly = lgefint(y);
585 262957 : if (ly == 2) return;
586 172642 : av = avma;
587 172642 : (void)new_chunk(3+ly+lgefint(*x)+nbits2lg(e)); /* HACK */
588 172642 : y = shifti(mului(u,y), e);
589 172642 : set_avma(av); return subzi(x, y);
590 : }
591 :
592 : /* x + u*y * 2^e */
593 : INLINE GEN
594 4076185 : addmuliu2n(GEN x, GEN y, ulong u, long e)
595 : {
596 : pari_sp av;
597 4076185 : long ly = lgefint(y);
598 4076185 : if (ly == 2) return x;
599 2798025 : av = avma;
600 2798025 : (void)new_chunk(3+ly+lgefint(x)+nbits2lg(e)); /* HACK */
601 2799738 : y = shifti(mului(u,y), e);
602 2798650 : set_avma(av); return addii(x, y);
603 : }
604 :
605 : /* *x += u*y * 2^e */
606 : INLINE void
607 270968 : addmulzu2n(GEN *x, GEN y, ulong u, long e)
608 : {
609 : pari_sp av;
610 270968 : long ly = lgefint(y);
611 270968 : if (ly == 2) return;
612 178036 : av = avma;
613 178036 : (void)new_chunk(3+ly+lgefint(*x)+nbits2lg(e)); /* HACK */
614 178036 : y = shifti(mului(u,y), e);
615 178036 : set_avma(av); return addzi(x, y);
616 : }
617 :
618 : /* n < 10; gerepileall supporting &NULL arguments. Maybe rename and export ? */
619 : INLINE void
620 4728 : gc_lll(pari_sp av, int n, ...)
621 : {
622 : int i, j;
623 : GEN *gptr[10];
624 : size_t s;
625 4728 : va_list a; va_start(a, n);
626 14184 : for (i=j=0; i<n; i++)
627 : {
628 9456 : GEN *x = va_arg(a,GEN*);
629 9456 : if (*x) { gptr[j++] = x; *x = (GEN)copy_bin(*x); }
630 : }
631 4728 : va_end(a); set_avma(av);
632 11782 : for (--j; j>=0; j--) *gptr[j] = bin_copy((GENbin*)*gptr[j]);
633 4728 : s = pari_mainstack->top - pari_mainstack->bot;
634 : /* size of saved objects ~ stacksize / 4 => overflow */
635 4728 : if (av - avma > (s >> 2))
636 : {
637 0 : size_t t = avma - pari_mainstack->bot;
638 0 : av = avma; new_chunk((s + t) / sizeof(long)); set_avma(av); /* double */
639 : }
640 4728 : }
641 :
642 : /********************************************************************/
643 : /** **/
644 : /** FPLLL (adapted from D. Stehle's code) **/
645 : /** **/
646 : /********************************************************************/
647 : /* Babai* and fplll* are a conversion to libpari API and data types
648 : of fplll-1.3 by Damien Stehle'.
649 :
650 : Copyright 2005, 2006 Damien Stehle'.
651 :
652 : This program is free software; you can redistribute it and/or modify it
653 : under the terms of the GNU General Public License as published by the
654 : Free Software Foundation; either version 2 of the License, or (at your
655 : option) any later version.
656 :
657 : This program implements ideas from the paper "Floating-point LLL Revisited",
658 : by Phong Nguyen and Damien Stehle', in the Proceedings of Eurocrypt'2005,
659 : Springer-Verlag; and was partly inspired by Shoup's NTL library:
660 : http://www.shoup.net/ntl/ */
661 :
662 : /* x t_REAL, |x| >= 1/2. Test whether |x| <= 3/2 */
663 : static int
664 397680 : absrsmall2(GEN x)
665 : {
666 397680 : long e = expo(x), l, i;
667 397680 : if (e < 0) return 1;
668 200040 : if (e > 0 || (ulong)x[2] > (3UL << (BITS_IN_LONG-2))) return 0;
669 : /* line above assumes l > 2. OK since x != 0 */
670 74094 : l = lg(x); for (i = 3; i < l; i++) if (x[i]) return 0;
671 63544 : return 1;
672 : }
673 : /* x t_REAL; test whether |x| <= 1/2 */
674 : static int
675 691417 : absrsmall(GEN x)
676 : {
677 : long e, l, i;
678 691417 : if (!signe(x)) return 1;
679 687251 : e = expo(x); if (e < -1) return 1;
680 403090 : if (e > -1 || (ulong)x[2] > HIGHBIT) return 0;
681 6119 : l = lg(x); for (i = 3; i < l; i++) if (x[i]) return 0;
682 5410 : return 1;
683 : }
684 :
685 : static void
686 31732232 : rotate(GEN A, long k2, long k)
687 : {
688 : long i;
689 31732232 : GEN B = gel(A,k2);
690 101396349 : for (i = k2; i > k; i--) gel(A,i) = gel(A,i-1);
691 31732232 : gel(A,k) = B;
692 31732232 : }
693 :
694 : /************************* FAST version (double) ************************/
695 : #define dmael(x,i,j) ((x)[i][j])
696 : #define del(x,i) ((x)[i])
697 :
698 : static double *
699 34481144 : cget_dblvec(long d)
700 34481144 : { return (double*) stack_malloc_align(d*sizeof(double), sizeof(double)); }
701 :
702 : static double **
703 8274572 : cget_dblmat(long d) { return (double **) cgetg(d, t_VECSMALL); }
704 :
705 : static double
706 161624524 : itodbl_exp(GEN x, long *e)
707 : {
708 161624524 : pari_sp av = avma;
709 161624524 : GEN r = itor(x,DEFAULTPREC);
710 161614605 : *e = expo(r); setexpo(r,0);
711 161610988 : return gc_double(av, rtodbl(r));
712 : }
713 :
714 : static double
715 117592379 : dbldotproduct(double *x, double *y, long n)
716 : {
717 : long i;
718 117592379 : double sum = del(x,1) * del(y,1);
719 1381246583 : for (i=2; i<=n; i++) sum += del(x,i) * del(y,i);
720 117592379 : return sum;
721 : }
722 :
723 : static double
724 2440341 : dbldotsquare(double *x, long n)
725 : {
726 : long i;
727 2440341 : double sum = del(x,1) * del(x,1);
728 8093449 : for (i=2; i<=n; i++) sum += del(x,i) * del(x,i);
729 2440341 : return sum;
730 : }
731 :
732 : static long
733 24695254 : set_line(double *appv, GEN v, long n)
734 : {
735 24695254 : long i, maxexp = 0;
736 24695254 : pari_sp av = avma;
737 24695254 : GEN e = cgetg(n+1, t_VECSMALL);
738 186305766 : for (i = 1; i <= n; i++)
739 : {
740 161622914 : del(appv,i) = itodbl_exp(gel(v,i), e+i);
741 161608937 : if (e[i] > maxexp) maxexp = e[i];
742 : }
743 186394028 : for (i = 1; i <= n; i++) del(appv,i) = ldexp(del(appv,i), e[i]-maxexp);
744 24682852 : set_avma(av); return maxexp;
745 : }
746 :
747 : static void
748 34416670 : dblrotate(double **A, long k2, long k)
749 : {
750 : long i;
751 34416670 : double *B = del(A,k2);
752 109020346 : for (i = k2; i > k; i--) del(A,i) = del(A,i-1);
753 34416670 : del(A,k) = B;
754 34416670 : }
755 : /* update G[kappa][i] from appB */
756 : static void
757 22459777 : setG_fast(double **appB, long n, double **G, long kappa, long a, long b)
758 : { long i;
759 101256466 : for (i = a; i <= b; i++)
760 78796757 : dmael(G,kappa,i) = dbldotproduct(del(appB,kappa), del(appB,i), n);
761 22459709 : }
762 : /* update G[i][kappa] from appB */
763 : static void
764 16904956 : setG2_fast(double **appB, long n, double **G, long kappa, long a, long b)
765 : { long i;
766 55704039 : for (i = a; i <= b; i++)
767 38799053 : dmael(G,i,kappa) = dbldotproduct(del(appB,kappa), del(appB,i), n);
768 16904986 : }
769 : const long EX0 = -2; /* uninitialized; any value less than expo(0.51) = -1 */
770 :
771 : #ifdef LONG_IS_64BIT
772 : typedef long s64;
773 : #define addmuliu64_inplace addmuliu_inplace
774 : #define submuliu64_inplace submuliu_inplace
775 : #define submuliu642n submuliu2n
776 : #define addmuliu642n addmuliu2n
777 : #else
778 : typedef long long s64;
779 : typedef unsigned long long u64;
780 :
781 : INLINE GEN
782 19651882 : u64toi(u64 x)
783 : {
784 : GEN y;
785 : ulong h;
786 19651882 : if (!x) return gen_0;
787 19651882 : h = x>>32;
788 19651882 : if (!h) return utoipos(x);
789 1141961 : y = cgetipos(4);
790 1141961 : *int_LSW(y) = x&0xFFFFFFFF;
791 1141961 : *int_MSW(y) = x>>32;
792 1141961 : return y;
793 : }
794 :
795 : INLINE GEN
796 672343 : u64toineg(u64 x)
797 : {
798 : GEN y;
799 : ulong h;
800 672343 : if (!x) return gen_0;
801 672343 : h = x>>32;
802 672343 : if (!h) return utoineg(x);
803 672343 : y = cgetineg(4);
804 672343 : *int_LSW(y) = x&0xFFFFFFFF;
805 672343 : *int_MSW(y) = x>>32;
806 672343 : return y;
807 : }
808 : INLINE GEN
809 9454199 : addmuliu64_inplace(GEN x, GEN y, u64 u) { return addmulii(x, y, u64toi(u)); }
810 :
811 : INLINE GEN
812 9509362 : submuliu64_inplace(GEN x, GEN y, u64 u) { return submulii(x, y, u64toi(u)); }
813 :
814 : INLINE GEN
815 672343 : addmuliu642n(GEN x, GEN y, u64 u, long e) { return submulshift(x, y, u64toineg(u), e); }
816 :
817 : INLINE GEN
818 688321 : submuliu642n(GEN x, GEN y, u64 u, long e) { return submulshift(x, y, u64toi(u), e); }
819 :
820 : #endif
821 :
822 : /* Babai's Nearest Plane algorithm (iterative); see Babai() */
823 : static int
824 30008234 : Babai_fast(pari_sp av, long kappa, GEN *pB, GEN *pU, double **mu, double **r,
825 : double *s, double **appB, GEN expoB, double **G,
826 : long a, long zeros, long maxG, double eta)
827 : {
828 30008234 : GEN B = *pB, U = *pU;
829 30008234 : const long n = nbrows(B), d = U ? lg(U)-1: 0;
830 30008123 : long k, aa = (a > zeros)? a : zeros+1;
831 30008123 : long emaxmu = EX0, emax2mu = EX0;
832 : s64 xx;
833 30008123 : int did_something = 0;
834 : /* N.B: we set d = 0 (resp. n = 0) to avoid updating U (resp. B) */
835 :
836 17111999 : for (;;) {
837 47120122 : int go_on = 0;
838 47120122 : long i, j, emax3mu = emax2mu;
839 :
840 47120122 : if (gc_needed(av,2))
841 : {
842 197 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[1], a=%ld", aa);
843 197 : gc_lll(av,2,&B,&U);
844 : }
845 : /* Step2: compute the GSO for stage kappa */
846 47119394 : emax2mu = emaxmu; emaxmu = EX0;
847 180589030 : for (j=aa; j<kappa; j++)
848 : {
849 133470496 : double g = dmael(G,kappa,j);
850 573275159 : for (k = zeros+1; k < j; k++) g -= dmael(mu,j,k) * dmael(r,kappa,k);
851 133470496 : dmael(r,kappa,j) = g;
852 133470496 : dmael(mu,kappa,j) = dmael(r,kappa,j) / dmael(r,j,j);
853 133470496 : emaxmu = maxss(emaxmu, expoB[kappa]-expoB[j]);
854 : }
855 : /* maxmu doesn't decrease fast enough */
856 47118534 : if (emax3mu != EX0 && emax3mu <= emax2mu + 5) {*pB = B; *pU = U; return 1;}
857 :
858 167656497 : for (j=kappa-1; j>zeros; j--)
859 : {
860 137654154 : double tmp = fabs(ldexp (dmael(mu,kappa,j), expoB[kappa]-expoB[j]));
861 137654154 : if (tmp>eta) { go_on = 1; break; }
862 : }
863 :
864 : /* Step3--5: compute the X_j's */
865 47114395 : if (go_on)
866 77690701 : for (j=kappa-1; j>zeros; j--)
867 : { /* The code below seemingly handles U = NULL, but in this case d = 0 */
868 60580056 : int e = expoB[j] - expoB[kappa];
869 60580056 : double tmp = ldexp(dmael(mu,kappa,j), -e), atmp = fabs(tmp);
870 : /* tmp = Inf is allowed */
871 60580056 : if (atmp <= .5) continue; /* size-reduced */
872 33958690 : if (gc_needed(av,2))
873 : {
874 346 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[2], a=%ld, j=%ld", aa,j);
875 346 : gc_lll(av,2,&B,&U);
876 : }
877 33960896 : did_something = 1;
878 : /* we consider separately the case |X| = 1 */
879 33960896 : if (atmp <= 1.5)
880 : {
881 22870710 : if (dmael(mu,kappa,j) > 0) { /* in this case, X = 1 */
882 46828987 : for (k=zeros+1; k<j; k++)
883 35156165 : dmael(mu,kappa,k) -= ldexp(dmael(mu,j,k), e);
884 157967247 : for (i=1; i<=n; i++)
885 146296124 : gmael(B,kappa,i) = subii(gmael(B,kappa,i), gmael(B,j,i));
886 104234806 : for (i=1; i<=d; i++)
887 92563495 : gmael(U,kappa,i) = subii(gmael(U,kappa,i), gmael(U,j,i));
888 : } else { /* otherwise X = -1 */
889 45998590 : for (k=zeros+1; k<j; k++)
890 34800702 : dmael(mu,kappa,k) += ldexp(dmael(mu,j,k), e);
891 155285754 : for (i=1; i<=n; i++)
892 144090057 : gmael(B,kappa,i) = addii(gmael(B,kappa,i), gmael(B,j,i));
893 101538127 : for (i=1; i<=d; i++)
894 90342280 : gmael(U,kappa,i) = addii(gmael(U,kappa,i), gmael(U,j,i));
895 : }
896 22867158 : continue;
897 : }
898 : /* we have |X| >= 2 */
899 11090186 : if (atmp < 9007199254740992.)
900 : {
901 10252293 : tmp = pari_rint(tmp);
902 24482714 : for (k=zeros+1; k<j; k++)
903 14230425 : dmael(mu,kappa,k) -= ldexp(tmp * dmael(mu,j,k), e);
904 10252289 : xx = (s64) tmp;
905 10252289 : if (xx > 0) /* = xx */
906 : {
907 46042709 : for (i=1; i<=n; i++)
908 40886722 : gmael(B,kappa,i) = submuliu64_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
909 33094782 : for (i=1; i<=d; i++)
910 27938784 : gmael(U,kappa,i) = submuliu64_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
911 : }
912 : else /* = -xx */
913 : {
914 45720819 : for (i=1; i<=n; i++)
915 40624651 : gmael(B,kappa,i) = addmuliu64_inplace(gmael(B,kappa,i), gmael(B,j,i), -xx);
916 32725433 : for (i=1; i<=d; i++)
917 27629235 : gmael(U,kappa,i) = addmuliu64_inplace(gmael(U,kappa,i), gmael(U,j,i), -xx);
918 : }
919 : }
920 : else
921 : {
922 : int E;
923 837893 : xx = (s64) ldexp(frexp(dmael(mu,kappa,j), &E), 53);
924 837893 : E -= e + 53;
925 837893 : if (E <= 0)
926 : {
927 0 : xx = xx << -E;
928 0 : for (k=zeros+1; k<j; k++)
929 0 : dmael(mu,kappa,k) -= ldexp(((double)xx) * dmael(mu,j,k), e);
930 0 : if (xx > 0) /* = xx */
931 : {
932 0 : for (i=1; i<=n; i++)
933 0 : gmael(B,kappa,i) = submuliu64_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
934 0 : for (i=1; i<=d; i++)
935 0 : gmael(U,kappa,i) = submuliu64_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
936 : }
937 : else /* = -xx */
938 : {
939 0 : for (i=1; i<=n; i++)
940 0 : gmael(B,kappa,i) = addmuliu64_inplace(gmael(B,kappa,i), gmael(B,j,i), -xx);
941 0 : for (i=1; i<=d; i++)
942 0 : gmael(U,kappa,i) = addmuliu64_inplace(gmael(U,kappa,i), gmael(U,j,i), -xx);
943 : }
944 : } else
945 : {
946 2788092 : for (k=zeros+1; k<j; k++)
947 1950199 : dmael(mu,kappa,k) -= ldexp(((double)xx) * dmael(mu,j,k), E + e);
948 837893 : if (xx > 0) /* = xx */
949 : {
950 3934383 : for (i=1; i<=n; i++)
951 3512622 : gmael(B,kappa,i) = submuliu642n(gmael(B,kappa,i), gmael(B,j,i), xx, E);
952 1509238 : for (i=1; i<=d; i++)
953 1087477 : gmael(U,kappa,i) = submuliu642n(gmael(U,kappa,i), gmael(U,j,i), xx, E);
954 : }
955 : else /* = -xx */
956 : {
957 3884768 : for (i=1; i<=n; i++)
958 3468603 : gmael(B,kappa,i) = addmuliu642n(gmael(B,kappa,i), gmael(B,j,i), -xx, E);
959 1485525 : for (i=1; i<=d; i++)
960 1069360 : gmael(U,kappa,i) = addmuliu642n(gmael(U,kappa,i), gmael(U,j,i), -xx, E);
961 : }
962 : }
963 : }
964 : }
965 47112985 : if (!go_on) break; /* Anything happened? */
966 17109225 : expoB[kappa] = set_line(del(appB,kappa), gel(B,kappa), n);
967 17111525 : setG_fast(appB, n, G, kappa, zeros+1, kappa-1);
968 17111999 : aa = zeros+1;
969 : }
970 30003760 : if (did_something) setG2_fast(appB, n, G, kappa, kappa, maxG);
971 :
972 30004068 : del(s,zeros+1) = dmael(G,kappa,kappa);
973 : /* the last s[kappa-1]=r[kappa][kappa] is computed only if kappa increases */
974 109281740 : for (k=zeros+1; k<=kappa-2; k++)
975 79277672 : del(s,k+1) = del(s,k) - dmael(mu,kappa,k)*dmael(r,kappa,k);
976 30004068 : *pB = B; *pU = U; return 0;
977 : }
978 :
979 : static void
980 11913844 : update_alpha(GEN alpha, long kappa, long kappa2, long kappamax)
981 : {
982 : long i;
983 37887426 : for (i = kappa; i < kappa2; i++)
984 25973582 : if (kappa <= alpha[i]) alpha[i] = kappa;
985 37887444 : for (i = kappa2; i > kappa; i--) alpha[i] = alpha[i-1];
986 23106229 : for (i = kappa2+1; i <= kappamax; i++)
987 11192385 : if (kappa < alpha[i]) alpha[i] = kappa;
988 11913844 : alpha[kappa] = kappa;
989 11913844 : }
990 : static void
991 441451 : rotateG(GEN G, long kappa2, long kappa, long maxG, GEN Gtmp)
992 : {
993 : long i, j;
994 3418583 : for (i=1; i<=kappa2; i++) gel(Gtmp,i) = gmael(G,kappa2,i);
995 1817157 : for ( ; i<=maxG; i++) gel(Gtmp,i) = gmael(G,i,kappa2);
996 1546732 : for (i=kappa2; i>kappa; i--)
997 : {
998 5249935 : for (j=1; j<kappa; j++) gmael(G,i,j) = gmael(G,i-1,j);
999 1105281 : gmael(G,i,kappa) = gel(Gtmp,i-1);
1000 3983763 : for (j=kappa+1; j<=i; j++) gmael(G,i,j) = gmael(G,i-1,j-1);
1001 4789895 : for (j=kappa2+1; j<=maxG; j++) gmael(G,j,i) = gmael(G,j,i-1);
1002 : }
1003 1871851 : for (i=1; i<kappa; i++) gmael(G,kappa,i) = gel(Gtmp,i);
1004 441451 : gmael(G,kappa,kappa) = gel(Gtmp,kappa2);
1005 1817157 : for (i=kappa2+1; i<=maxG; i++) gmael(G,i,kappa) = gel(Gtmp,i);
1006 441451 : }
1007 : static void
1008 11472263 : rotateG_fast(double **G, long kappa2, long kappa, long maxG, double *Gtmp)
1009 : {
1010 : long i, j;
1011 66668512 : for (i=1; i<=kappa2; i++) del(Gtmp,i) = dmael(G,kappa2,i);
1012 22151265 : for ( ; i<=maxG; i++) del(Gtmp,i) = dmael(G,i,kappa2);
1013 36340225 : for (i=kappa2; i>kappa; i--)
1014 : {
1015 69475187 : for (j=1; j<kappa; j++) dmael(G,i,j) = dmael(G,i-1,j);
1016 24867962 : dmael(G,i,kappa) = del(Gtmp,i-1);
1017 84736890 : for (j=kappa+1; j<=i; j++) dmael(G,i,j) = dmael(G,i-1,j-1);
1018 46853357 : for (j=kappa2+1; j<=maxG; j++) dmael(G,j,i) = dmael(G,j,i-1);
1019 : }
1020 30328448 : for (i=1; i<kappa; i++) dmael(G,kappa,i) = del(Gtmp,i);
1021 11472263 : dmael(G,kappa,kappa) = del(Gtmp,kappa2);
1022 22151317 : for (i=kappa2+1; i<=maxG; i++) dmael(G,i,kappa) = del(Gtmp,i);
1023 11472263 : }
1024 :
1025 : /* LLL-reduces (B,U) in place [apply base change transforms to B and U].
1026 : * Gram matrix, and GSO performed on matrices of 'double'.
1027 : * If (keepfirst), never swap with first vector.
1028 : * Return -1 on failure, else zeros = dim Kernel (>= 0) */
1029 : static long
1030 2068648 : fplll_fast(GEN *pB, GEN *pU, double delta, double eta, long keepfirst)
1031 : {
1032 : pari_sp av;
1033 : long kappa, kappa2, d, n, i, j, zeros, kappamax, maxG;
1034 : double **mu, **r, *s, tmp, *Gtmp, **G, **appB;
1035 2068648 : GEN alpha, expoB, B = *pB, U;
1036 2068648 : long cnt = 0;
1037 :
1038 2068648 : d = lg(B)-1;
1039 2068648 : n = nbrows(B);
1040 2068651 : U = *pU; /* NULL if inplace */
1041 :
1042 2068651 : G = cget_dblmat(d+1);
1043 2068650 : appB = cget_dblmat(d+1);
1044 2068646 : mu = cget_dblmat(d+1);
1045 2068644 : r = cget_dblmat(d+1);
1046 2068645 : s = cget_dblvec(d+1);
1047 9654667 : for (j = 1; j <= d; j++)
1048 : {
1049 7586019 : del(mu,j) = cget_dblvec(d+1);
1050 7586012 : del(r,j) = cget_dblvec(d+1);
1051 7586010 : del(appB,j) = cget_dblvec(n+1);
1052 7586008 : del(G,j) = cget_dblvec(d+1);
1053 46808743 : for (i=1; i<=d; i++) dmael(G,j,i) = 0.;
1054 : }
1055 2068648 : expoB = cgetg(d+1, t_VECSMALL);
1056 9654559 : for (i=1; i<=d; i++) expoB[i] = set_line(del(appB,i), gel(B,i), n);
1057 2068602 : Gtmp = cget_dblvec(d+1);
1058 2068642 : alpha = cgetg(d+1, t_VECSMALL);
1059 2068637 : av = avma;
1060 :
1061 : /* Step2: Initializing the main loop */
1062 2068637 : kappamax = 1;
1063 2068637 : i = 1;
1064 2068637 : maxG = d; /* later updated to kappamax */
1065 :
1066 : do {
1067 2234296 : dmael(G,i,i) = dbldotsquare(del(appB,i),n);
1068 2234300 : } while (dmael(G,i,i) <= 0 && (++i <=d));
1069 2068641 : zeros = i-1; /* all vectors B[i] with i <= zeros are zero vectors */
1070 2068641 : kappa = i;
1071 2068641 : if (zeros < d) dmael(r,zeros+1,zeros+1) = dmael(G,zeros+1,zeros+1);
1072 9489005 : for (i=zeros+1; i<=d; i++) alpha[i]=1;
1073 32072759 : while (++kappa <= d)
1074 : {
1075 30008243 : if (kappa > kappamax)
1076 : {
1077 5348170 : if (DEBUGLEVEL>=4) err_printf("K%ld ",kappa);
1078 5348170 : maxG = kappamax = kappa;
1079 5348170 : setG_fast(appB, n, G, kappa, zeros+1, kappa);
1080 : }
1081 : /* Step3: Call to the Babai algorithm, mu,r,s updated in place */
1082 30008241 : if (Babai_fast(av, kappa, &B,&U, mu,r,s, appB, expoB, G, alpha[kappa],
1083 4139 : zeros, maxG, eta)) { *pB=B; *pU=U; return -1; }
1084 :
1085 30004046 : tmp = ldexp(r[kappa-1][kappa-1] * delta, 2*(expoB[kappa-1]-expoB[kappa]));
1086 30004046 : if ((keepfirst && kappa == 2) || tmp <= del(s,kappa-1))
1087 : { /* Step4: Success of Lovasz's condition */
1088 18531662 : alpha[kappa] = kappa;
1089 18531662 : tmp = dmael(mu,kappa,kappa-1) * dmael(r,kappa,kappa-1);
1090 18531662 : dmael(r,kappa,kappa) = del(s,kappa-1)- tmp;
1091 18531662 : continue;
1092 : }
1093 : /* Step5: Find the right insertion index kappa, kappa2 = initial kappa */
1094 11472384 : if (DEBUGLEVEL>=4 && kappa==kappamax && del(s,kappa-1)!=0)
1095 0 : if (++cnt > 20) { cnt = 0; err_printf("(%ld) ", 2*expoB[1] + dblexpo(del(s,1))); }
1096 11472384 : kappa2 = kappa;
1097 : do {
1098 24868279 : kappa--;
1099 24868279 : if (kappa<zeros+2 + (keepfirst ? 1: 0)) break;
1100 18328348 : tmp = dmael(r,kappa-1,kappa-1) * delta;
1101 18328348 : tmp = ldexp(tmp, 2*(expoB[kappa-1]-expoB[kappa2]));
1102 18328348 : } while (del(s,kappa-1) <= tmp);
1103 11472384 : update_alpha(alpha, kappa, kappa2, kappamax);
1104 :
1105 : /* Step6: Update the mu's and r's */
1106 11472395 : dblrotate(mu,kappa2,kappa);
1107 11472378 : dblrotate(r,kappa2,kappa);
1108 11472327 : dmael(r,kappa,kappa) = del(s,kappa);
1109 :
1110 : /* Step7: Update B, appB, U, G */
1111 11472327 : rotate(B,kappa2,kappa);
1112 11472320 : dblrotate(appB,kappa2,kappa);
1113 11472288 : if (U) rotate(U,kappa2,kappa);
1114 11472285 : rotate(expoB,kappa2,kappa);
1115 11472267 : rotateG_fast(G,kappa2,kappa, maxG, Gtmp);
1116 :
1117 : /* Step8: Prepare the next loop iteration */
1118 11472456 : if (kappa == zeros+1 && dmael(G,kappa,kappa)<= 0)
1119 : {
1120 206044 : zeros++; kappa++;
1121 206044 : dmael(G,kappa,kappa) = dbldotsquare(del(appB,kappa),n);
1122 206044 : dmael(r,kappa,kappa) = dmael(G,kappa,kappa);
1123 : }
1124 : }
1125 2064516 : *pB = B; *pU = U; return zeros;
1126 : }
1127 :
1128 : /***************** HEURISTIC version (reduced precision) ****************/
1129 : static GEN
1130 195019 : realsqrdotproduct(GEN x)
1131 : {
1132 195019 : long i, l = lg(x);
1133 195019 : GEN z = sqrr(gel(x,1));
1134 1261611 : for (i=2; i<l; i++) z = addrr(z, sqrr(gel(x,i)));
1135 195019 : return z;
1136 : }
1137 : /* x, y non-empty vector of t_REALs, same length */
1138 : static GEN
1139 1174069 : realdotproduct(GEN x, GEN y)
1140 : {
1141 : long i, l;
1142 : GEN z;
1143 1174069 : if (x == y) return realsqrdotproduct(x);
1144 979050 : l = lg(x); z = mulrr(gel(x,1),gel(y,1));
1145 8498392 : for (i=2; i<l; i++) z = addrr(z, mulrr(gel(x,i), gel(y,i)));
1146 979050 : return z;
1147 : }
1148 : static void
1149 203188 : setG_heuristic(GEN appB, GEN G, long kappa, long a, long b)
1150 203188 : { pari_sp av = avma;
1151 : long i;
1152 925332 : for (i = a; i <= b; i++)
1153 722144 : affrr(realdotproduct(gel(appB,kappa),gel(appB,i)), gmael(G,kappa,i));
1154 203188 : set_avma(av);
1155 203188 : }
1156 : static void
1157 185331 : setG2_heuristic(GEN appB, GEN G, long kappa, long a, long b)
1158 185331 : { pari_sp av = avma;
1159 : long i;
1160 637256 : for (i = a; i <= b; i++)
1161 451925 : affrr(realdotproduct(gel(appB,kappa),gel(appB,i)), gmael(G,i,kappa));
1162 185331 : set_avma(av);
1163 185331 : }
1164 :
1165 : /* approximate t_REAL x as m * 2^e, where |m| < 2^bit */
1166 : static GEN
1167 12119 : truncexpo(GEN x, long bit, long *e)
1168 : {
1169 12119 : *e = expo(x) + 1 - bit;
1170 12119 : if (*e >= 0) return mantissa2nr(x, 0);
1171 928 : *e = 0; return roundr_safe(x);
1172 : }
1173 : /* Babai's Nearest Plane algorithm (iterative); see Babai() */
1174 : static int
1175 285744 : Babai_heuristic(pari_sp av, long kappa, GEN *pB, GEN *pU, GEN mu, GEN r, GEN s,
1176 : GEN appB, GEN G, long a, long zeros, long maxG,
1177 : GEN eta, long prec)
1178 : {
1179 285744 : GEN B = *pB, U = *pU;
1180 285744 : const long n = nbrows(B), d = U ? lg(U)-1: 0, bit = prec2nbits(prec);
1181 285744 : long k, aa = (a > zeros)? a : zeros+1;
1182 285744 : int did_something = 0;
1183 285744 : long emaxmu = EX0, emax2mu = EX0;
1184 : /* N.B: we set d = 0 (resp. n = 0) to avoid updating U (resp. B) */
1185 :
1186 193500 : for (;;) {
1187 479244 : int go_on = 0;
1188 479244 : long i, j, emax3mu = emax2mu;
1189 :
1190 479244 : if (gc_needed(av,2))
1191 : {
1192 27 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[1], a=%ld", aa);
1193 27 : gc_lll(av,2,&B,&U);
1194 : }
1195 : /* Step2: compute the GSO for stage kappa */
1196 479244 : emax2mu = emaxmu; emaxmu = EX0;
1197 1846581 : for (j=aa; j<kappa; j++)
1198 : {
1199 1367337 : pari_sp btop = avma;
1200 1367337 : GEN g = gmael(G,kappa,j);
1201 4443130 : for (k = zeros+1; k<j; k++)
1202 3075793 : g = subrr(g, mulrr(gmael(mu,j,k), gmael(r,kappa,k)));
1203 1367337 : affrr(g, gmael(r,kappa,j));
1204 1367337 : affrr(divrr(gmael(r,kappa,j), gmael(r,j,j)), gmael(mu,kappa,j));
1205 1367337 : emaxmu = maxss(emaxmu, expo(gmael(mu,kappa,j)));
1206 1367337 : set_avma(btop);
1207 : }
1208 479244 : if (emax3mu != EX0 && emax3mu <= emax2mu + 5)
1209 1327 : { *pB = B; *pU = U; return 1; }
1210 :
1211 1624954 : for (j=kappa-1; j>zeros; j--)
1212 1340537 : if (abscmprr(gmael(mu,kappa,j), eta) > 0) { go_on = 1; break; }
1213 :
1214 : /* Step3--5: compute the X_j's */
1215 477917 : if (go_on)
1216 884917 : for (j=kappa-1; j>zeros; j--)
1217 : { /* The code below seemingly handles U = NULL, but in this case d = 0 */
1218 : pari_sp btop;
1219 691417 : GEN tmp = gmael(mu,kappa,j);
1220 691417 : if (absrsmall(tmp)) continue; /* size-reduced */
1221 :
1222 397680 : if (gc_needed(av,2))
1223 : {
1224 19 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[2], a=%ld, j=%ld", aa,j);
1225 19 : gc_lll(av,2,&B,&U);
1226 : }
1227 397680 : btop = avma; did_something = 1;
1228 : /* we consider separately the case |X| = 1 */
1229 397680 : if (absrsmall2(tmp))
1230 : {
1231 261184 : if (signe(tmp) > 0) { /* in this case, X = 1 */
1232 384089 : for (k=zeros+1; k<j; k++)
1233 254017 : affrr(subrr(gmael(mu,kappa,k), gmael(mu,j,k)), gmael(mu,kappa,k));
1234 130072 : set_avma(btop);
1235 1233168 : for (i=1; i<=n; i++)
1236 1103096 : gmael(B,kappa,i) = subii(gmael(B,kappa,i), gmael(B,j,i));
1237 805717 : for (i=1; i<=d; i++)
1238 675645 : gmael(U,kappa,i) = subii(gmael(U,kappa,i), gmael(U,j,i));
1239 : } else { /* otherwise X = -1 */
1240 388969 : for (k=zeros+1; k<j; k++)
1241 257857 : affrr(addrr(gmael(mu,kappa,k), gmael(mu,j,k)), gmael(mu,kappa,k));
1242 131112 : set_avma(btop);
1243 1247296 : for (i=1; i<=n; i++)
1244 1116184 : gmael(B,kappa,i) = addii(gmael(B,kappa,i), gmael(B,j,i));
1245 807308 : for (i=1; i<=d; i++)
1246 676196 : gmael(U,kappa,i) = addii(gmael(U,kappa,i),gmael(U,j,i));
1247 : }
1248 261184 : continue;
1249 : }
1250 : /* we have |X| >= 2 */
1251 136496 : if (expo(tmp) < BITS_IN_LONG)
1252 : {
1253 124377 : ulong xx = roundr_safe(tmp)[2]; /* X fits in an ulong */
1254 124377 : if (signe(tmp) > 0) /* = xx */
1255 : {
1256 136522 : for (k=zeros+1; k<j; k++)
1257 73806 : affrr(subrr(gmael(mu,kappa,k), mulur(xx, gmael(mu,j,k))),
1258 73806 : gmael(mu,kappa,k));
1259 62716 : set_avma(btop);
1260 416034 : for (i=1; i<=n; i++)
1261 353318 : gmael(B,kappa,i) = submuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
1262 307410 : for (i=1; i<=d; i++)
1263 244694 : gmael(U,kappa,i) = submuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
1264 : }
1265 : else /* = -xx */
1266 : {
1267 133239 : for (k=zeros+1; k<j; k++)
1268 71578 : affrr(addrr(gmael(mu,kappa,k), mulur(xx, gmael(mu,j,k))),
1269 71578 : gmael(mu,kappa,k));
1270 61661 : set_avma(btop);
1271 409060 : for (i=1; i<=n; i++)
1272 347399 : gmael(B,kappa,i) = addmuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
1273 294178 : for (i=1; i<=d; i++)
1274 232517 : gmael(U,kappa,i) = addmuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
1275 : }
1276 : }
1277 : else
1278 : {
1279 : long e;
1280 12119 : GEN X = truncexpo(tmp, bit, &e); /* tmp ~ X * 2^e */
1281 12119 : btop = avma;
1282 29481 : for (k=zeros+1; k<j; k++)
1283 : {
1284 17362 : GEN x = mulir(X, gmael(mu,j,k));
1285 17362 : if (e) shiftr_inplace(x, e);
1286 17362 : affrr(subrr(gmael(mu,kappa,k), x), gmael(mu,kappa,k));
1287 : }
1288 12119 : set_avma(btop);
1289 99017 : for (i=1; i<=n; i++)
1290 86898 : gmael(B,kappa,i) = submulshift(gmael(B,kappa,i), gmael(B,j,i), X, e);
1291 72999 : for (i=1; i<=d; i++)
1292 60880 : gmael(U,kappa,i) = submulshift(gmael(U,kappa,i), gmael(U,j,i), X, e);
1293 : }
1294 : }
1295 477917 : if (!go_on) break; /* Anything happened? */
1296 1448651 : for (i=1 ; i<=n; i++) affir(gmael(B,kappa,i), gmael(appB,kappa,i));
1297 193500 : setG_heuristic(appB, G, kappa, zeros+1, kappa-1);
1298 193500 : aa = zeros+1;
1299 : }
1300 284417 : if (did_something) setG2_heuristic(appB, G, kappa, kappa, maxG);
1301 284417 : affrr(gmael(G,kappa,kappa), gel(s,zeros+1));
1302 : /* the last s[kappa-1]=r[kappa][kappa] is computed only if kappa increases */
1303 284417 : av = avma;
1304 1018482 : for (k=zeros+1; k<=kappa-2; k++)
1305 734065 : affrr(subrr(gel(s,k), mulrr(gmael(mu,kappa,k), gmael(r,kappa,k))),
1306 734065 : gel(s,k+1));
1307 284417 : *pB = B; *pU = U; return gc_bool(av, 0);
1308 : }
1309 :
1310 : static GEN
1311 15729 : ZC_to_RC(GEN x, long prec)
1312 103309 : { pari_APPLY_type(t_COL,itor(gel(x,i),prec)) }
1313 :
1314 : static GEN
1315 4139 : ZM_to_RM(GEN x, long prec)
1316 19868 : { pari_APPLY_same(ZC_to_RC(gel(x,i),prec)) }
1317 :
1318 : /* LLL-reduces (B,U) in place [apply base change transforms to B and U].
1319 : * Gram matrix made of t_REAL at precision prec2, performe GSO at prec.
1320 : * If (keepfirst), never swap with first vector.
1321 : * Return -1 on failure, else zeros = dim Kernel (>= 0) */
1322 : static long
1323 4139 : fplll_heuristic(GEN *pB, GEN *pU, double DELTA, double ETA, long keepfirst,
1324 : long prec, long prec2)
1325 : {
1326 : pari_sp av, av2;
1327 : long kappa, kappa2, d, i, j, zeros, kappamax, maxG;
1328 4139 : GEN mu, r, s, tmp, Gtmp, alpha, G, appB, B = *pB, U;
1329 4139 : GEN delta = dbltor(DELTA), eta = dbltor(ETA);
1330 4139 : long cnt = 0;
1331 :
1332 4139 : d = lg(B)-1;
1333 4139 : U = *pU; /* NULL if inplace */
1334 :
1335 4139 : G = cgetg(d+1, t_MAT);
1336 4139 : mu = cgetg(d+1, t_MAT);
1337 4139 : r = cgetg(d+1, t_MAT);
1338 4139 : s = cgetg(d+1, t_VEC);
1339 4139 : appB = ZM_to_RM(B, prec2);
1340 19868 : for (j = 1; j <= d; j++)
1341 : {
1342 15729 : GEN M = cgetg(d+1, t_COL), R = cgetg(d+1, t_COL), S = cgetg(d+1, t_COL);
1343 15729 : gel(mu,j)= M;
1344 15729 : gel(r,j) = R;
1345 15729 : gel(G,j) = S;
1346 15729 : gel(s,j) = cgetr(prec);
1347 95060 : for (i = 1; i <= d; i++)
1348 : {
1349 79331 : gel(R,i) = cgetr(prec);
1350 79331 : gel(M,i) = cgetr(prec);
1351 79331 : gel(S,i) = cgetr(prec2);
1352 : }
1353 : }
1354 4139 : Gtmp = cgetg(d+1, t_VEC);
1355 4139 : alpha = cgetg(d+1, t_VECSMALL);
1356 4139 : av = avma;
1357 :
1358 : /* Step2: Initializing the main loop */
1359 4139 : kappamax = 1;
1360 4139 : i = 1;
1361 4139 : maxG = d; /* later updated to kappamax */
1362 :
1363 : do {
1364 4142 : affrr(RgV_dotsquare(gel(appB,i)), gmael(G,i,i));
1365 4142 : } while (signe(gmael(G,i,i)) == 0 && (++i <=d));
1366 4139 : zeros = i-1; /* all vectors B[i] with i <= zeros are zero vectors */
1367 4139 : kappa = i;
1368 4139 : if (zeros < d) affrr(gmael(G,zeros+1,zeros+1), gmael(r,zeros+1,zeros+1));
1369 19865 : for (i=zeros+1; i<=d; i++) alpha[i]=1;
1370 :
1371 288556 : while (++kappa <= d)
1372 : {
1373 285744 : if (kappa > kappamax)
1374 : {
1375 9688 : if (DEBUGLEVEL>=4) err_printf("K%ld ",kappa);
1376 9688 : maxG = kappamax = kappa;
1377 9688 : setG_heuristic(appB, G, kappa, zeros+1, kappa);
1378 : }
1379 : /* Step3: Call to the Babai algorithm, mu,r,s updated in place */
1380 285744 : if (Babai_heuristic(av, kappa, &B,&U, mu,r,s, appB, G, alpha[kappa], zeros,
1381 1327 : maxG, eta, prec)) { *pB = B; *pU = U; return -1; }
1382 284417 : av2 = avma;
1383 568724 : if ((keepfirst && kappa == 2) ||
1384 284307 : cmprr(mulrr(gmael(r,kappa-1,kappa-1), delta), gel(s,kappa-1)) <= 0)
1385 : { /* Step4: Success of Lovasz's condition */
1386 168661 : alpha[kappa] = kappa;
1387 168661 : tmp = mulrr(gmael(mu,kappa,kappa-1), gmael(r,kappa,kappa-1));
1388 168661 : affrr(subrr(gel(s,kappa-1), tmp), gmael(r,kappa,kappa));
1389 168661 : set_avma(av2); continue;
1390 : }
1391 : /* Step5: Find the right insertion index kappa, kappa2 = initial kappa */
1392 115756 : if (DEBUGLEVEL>=4 && kappa==kappamax && signe(gel(s,kappa-1)))
1393 0 : if (++cnt > 20) { cnt = 0; err_printf("(%ld) ", expo(gel(s,1))); }
1394 115756 : kappa2 = kappa;
1395 : do {
1396 276081 : kappa--;
1397 276081 : if (kappa < zeros+2 + (keepfirst ? 1: 0)) break;
1398 246553 : tmp = mulrr(gmael(r,kappa-1,kappa-1), delta);
1399 246553 : } while (cmprr(gel(s,kappa-1), tmp) <= 0 );
1400 115756 : set_avma(av2);
1401 115756 : update_alpha(alpha, kappa, kappa2, kappamax);
1402 :
1403 : /* Step6: Update the mu's and r's */
1404 115756 : rotate(mu,kappa2,kappa);
1405 115756 : rotate(r,kappa2,kappa);
1406 115756 : affrr(gel(s,kappa), gmael(r,kappa,kappa));
1407 :
1408 : /* Step7: Update B, appB, U, G */
1409 115756 : rotate(B,kappa2,kappa);
1410 115756 : rotate(appB,kappa2,kappa);
1411 115756 : if (U) rotate(U,kappa2,kappa);
1412 115756 : rotateG(G,kappa2,kappa, maxG, Gtmp);
1413 :
1414 : /* Step8: Prepare the next loop iteration */
1415 115756 : if (kappa == zeros+1 && !signe(gmael(G,kappa,kappa)))
1416 : {
1417 0 : zeros++; kappa++;
1418 0 : affrr(RgV_dotsquare(gel(appB,kappa)), gmael(G,kappa,kappa));
1419 0 : affrr(gmael(G,kappa,kappa), gmael(r,kappa,kappa));
1420 : }
1421 : }
1422 2812 : *pB=B; *pU=U; return zeros;
1423 : }
1424 :
1425 : /************************* PROVED version (t_INT) ***********************/
1426 : /* dpe inspired by dpe.h by Patrick Pelissier, Paul Zimmermann
1427 : * https://gforge.inria.fr/projects/dpe/
1428 : */
1429 :
1430 : typedef struct
1431 : {
1432 : double d; /* significand */
1433 : long e; /* exponent */
1434 : } dpe_t;
1435 :
1436 : #define Dmael(x,i,j) (&((x)[i][j]))
1437 : #define Del(x,i) (&((x)[i]))
1438 :
1439 : static void
1440 651390 : dperotate(dpe_t **A, long k2, long k)
1441 : {
1442 : long i;
1443 651390 : dpe_t *B = A[k2];
1444 2309790 : for (i = k2; i > k; i--) A[i] = A[i-1];
1445 651390 : A[k] = B;
1446 651390 : }
1447 :
1448 : static void
1449 108014810 : dpe_normalize0(dpe_t *x)
1450 : {
1451 : int e;
1452 108014810 : x->d = frexp(x->d, &e);
1453 108014810 : x->e += e;
1454 108014810 : }
1455 :
1456 : static void
1457 47895330 : dpe_normalize(dpe_t *x)
1458 : {
1459 47895330 : if (x->d == 0.0)
1460 496220 : x->e = -LONG_MAX;
1461 : else
1462 47399110 : dpe_normalize0(x);
1463 47895385 : }
1464 :
1465 : static GEN
1466 93935 : dpetor(dpe_t *x)
1467 : {
1468 93935 : GEN r = dbltor(x->d);
1469 93935 : if (signe(r)==0) return r;
1470 93935 : setexpo(r, x->e-1);
1471 93935 : return r;
1472 : }
1473 :
1474 : static void
1475 25559531 : affdpe(dpe_t *y, dpe_t *x)
1476 : {
1477 25559531 : x->d = y->d;
1478 25559531 : x->e = y->e;
1479 25559531 : }
1480 :
1481 : static void
1482 20520494 : affidpe(GEN y, dpe_t *x)
1483 : {
1484 20520494 : pari_sp av = avma;
1485 20520494 : GEN r = itor(y, DEFAULTPREC);
1486 20520175 : x->e = expo(r)+1;
1487 20520175 : setexpo(r,-1);
1488 20520122 : x->d = rtodbl(r);
1489 20520061 : set_avma(av);
1490 20520024 : }
1491 :
1492 : static void
1493 3135908 : affdbldpe(double y, dpe_t *x)
1494 : {
1495 3135908 : x->d = (double)y;
1496 3135908 : x->e = 0;
1497 3135908 : dpe_normalize(x);
1498 3135907 : }
1499 :
1500 : static void
1501 56705913 : dpe_mulz(dpe_t *x, dpe_t *y, dpe_t *z)
1502 : {
1503 56705913 : z->d = x->d * y->d;
1504 56705913 : if (z->d == 0.0)
1505 8143955 : z->e = -LONG_MAX;
1506 : else
1507 : {
1508 48561958 : z->e = x->e + y->e;
1509 48561958 : dpe_normalize0(z);
1510 : }
1511 56706139 : }
1512 :
1513 : static void
1514 13954654 : dpe_divz(dpe_t *x, dpe_t *y, dpe_t *z)
1515 : {
1516 13954654 : z->d = x->d / y->d;
1517 13954654 : if (z->d == 0.0)
1518 1900461 : z->e = -LONG_MAX;
1519 : else
1520 : {
1521 12054193 : z->e = x->e - y->e;
1522 12054193 : dpe_normalize0(z);
1523 : }
1524 13954700 : }
1525 :
1526 : static void
1527 244136 : dpe_negz(dpe_t *y, dpe_t *x)
1528 : {
1529 244136 : x->d = - y->d;
1530 244136 : x->e = y->e;
1531 244136 : }
1532 :
1533 : static void
1534 1941878 : dpe_addz(dpe_t *y, dpe_t *z, dpe_t *x)
1535 : {
1536 1941878 : if (y->e > z->e + 53)
1537 112039 : affdpe(y, x);
1538 1829839 : else if (z->e > y->e + 53)
1539 41670 : affdpe(z, x);
1540 : else
1541 : {
1542 1788169 : long d = y->e - z->e;
1543 :
1544 1788169 : if (d >= 0)
1545 : {
1546 1344057 : x->d = y->d + ldexp(z->d, -d);
1547 1344057 : x->e = y->e;
1548 : }
1549 : else
1550 : {
1551 444112 : x->d = z->d + ldexp(y->d, d);
1552 444112 : x->e = z->e;
1553 : }
1554 1788169 : dpe_normalize(x);
1555 : }
1556 1941885 : }
1557 : static void
1558 53541777 : dpe_subz(dpe_t *y, dpe_t *z, dpe_t *x)
1559 : {
1560 53541777 : if (y->e > z->e + 53)
1561 11125585 : affdpe(y, x);
1562 42416192 : else if (z->e > y->e + 53)
1563 244136 : dpe_negz(z, x);
1564 : else
1565 : {
1566 42172056 : long d = y->e - z->e;
1567 :
1568 42172056 : if (d >= 0)
1569 : {
1570 39412877 : x->d = y->d - ldexp(z->d, -d);
1571 39412877 : x->e = y->e;
1572 : }
1573 : else
1574 : {
1575 2759179 : x->d = ldexp(y->d, d) - z->d;
1576 2759179 : x->e = z->e;
1577 : }
1578 42172056 : dpe_normalize(x);
1579 : }
1580 53542037 : }
1581 :
1582 : static void
1583 799034 : dpe_muluz(dpe_t *y, ulong t, dpe_t *x)
1584 : {
1585 799034 : x->d = y->d * (double)t;
1586 799034 : x->e = y->e;
1587 799034 : dpe_normalize(x);
1588 799034 : }
1589 :
1590 : static void
1591 342647 : dpe_addmuluz(dpe_t *y, dpe_t *z, ulong t, dpe_t *x)
1592 : {
1593 : dpe_t tmp;
1594 342647 : dpe_muluz(z, t, &tmp);
1595 342647 : dpe_addz(y, &tmp, x);
1596 342647 : }
1597 :
1598 : static void
1599 411904 : dpe_submuluz(dpe_t *y, dpe_t *z, ulong t, dpe_t *x)
1600 : {
1601 : dpe_t tmp;
1602 411904 : dpe_muluz(z, t, &tmp);
1603 411904 : dpe_subz(y, &tmp, x);
1604 411904 : }
1605 :
1606 : static void
1607 51511860 : dpe_submulz(dpe_t *y, dpe_t *z, dpe_t *t, dpe_t *x)
1608 : {
1609 : dpe_t tmp;
1610 51511860 : dpe_mulz(z, t, &tmp);
1611 51511874 : dpe_subz(y, &tmp, x);
1612 51511923 : }
1613 :
1614 : static int
1615 5194183 : dpe_cmp(dpe_t *x, dpe_t *y)
1616 : {
1617 5194183 : int sx = x->d < 0. ? -1: x->d > 0.;
1618 5194183 : int sy = y->d < 0. ? -1: y->d > 0.;
1619 5194183 : int d = sx - sy;
1620 :
1621 5194183 : if (d != 0)
1622 141619 : return d;
1623 5052564 : else if (x->e > y->e)
1624 481157 : return (sx > 0) ? 1 : -1;
1625 4571407 : else if (y->e > x->e)
1626 2501617 : return (sx > 0) ? -1 : 1;
1627 : else
1628 2069790 : return (x->d < y->d) ? -1 : (x->d > y->d);
1629 : }
1630 :
1631 : static int
1632 14395511 : dpe_abscmp(dpe_t *x, dpe_t *y)
1633 : {
1634 14395511 : if (x->e > y->e)
1635 271158 : return 1;
1636 14124353 : else if (y->e > x->e)
1637 13273984 : return -1;
1638 : else
1639 850369 : return (fabs(x->d) < fabs(y->d)) ? -1 : (fabs(x->d) > fabs(y->d));
1640 : }
1641 :
1642 : static int
1643 1387783 : dpe_abssmall(dpe_t *x)
1644 : {
1645 1387783 : return (x->e <= 0) || (x->e == 1 && fabs(x->d) <= .75);
1646 : }
1647 :
1648 : static int
1649 5194177 : dpe_cmpmul(dpe_t *x, dpe_t *y, dpe_t *z)
1650 : {
1651 : dpe_t t;
1652 5194177 : dpe_mulz(x,y,&t);
1653 5194182 : return dpe_cmp(&t, z);
1654 : }
1655 :
1656 : static dpe_t *
1657 13041752 : cget_dpevec(long d)
1658 13041752 : { return (dpe_t*) stack_malloc_align(d*sizeof(dpe_t), sizeof(dpe_t)); }
1659 :
1660 : static dpe_t **
1661 3135899 : cget_dpemat(long d) { return (dpe_t **) cgetg(d, t_VECSMALL); }
1662 :
1663 : static GEN
1664 20252 : dpeM_diagonal_shallow(dpe_t **m, long d)
1665 : {
1666 : long i;
1667 20252 : GEN y = cgetg(d+1,t_VEC);
1668 114187 : for (i=1; i<=d; i++) gel(y, i) = dpetor(Dmael(m,i,i));
1669 20252 : return y;
1670 : }
1671 :
1672 : static void
1673 1387761 : affii_or_copy_gc(pari_sp av, GEN x, GEN *y)
1674 : {
1675 1387761 : long l = lg(*y);
1676 1387761 : if (lgefint(x) <= l && isonstack(*y))
1677 : {
1678 1387745 : affii(x,*y);
1679 1387746 : set_avma(av);
1680 : }
1681 : else
1682 12 : *y = gerepileuptoint(av, x);
1683 1387762 : }
1684 :
1685 : /* *x -= u*y */
1686 : INLINE void
1687 5919970 : submulziu(GEN *x, GEN y, ulong u)
1688 : {
1689 : pari_sp av;
1690 5919970 : long ly = lgefint(y);
1691 5919970 : if (ly == 2) return;
1692 3251885 : av = avma;
1693 3251885 : (void)new_chunk(3+ly+lgefint(*x)); /* HACK */
1694 3251960 : y = mului(u,y);
1695 3251938 : set_avma(av); subzi(x, y);
1696 : }
1697 :
1698 : /* *x += u*y */
1699 : INLINE void
1700 4576862 : addmulziu(GEN *x, GEN y, ulong u)
1701 : {
1702 : pari_sp av;
1703 4576862 : long ly = lgefint(y);
1704 4576862 : if (ly == 2) return;
1705 2755323 : av = avma;
1706 2755323 : (void)new_chunk(3+ly+lgefint(*x)); /* HACK */
1707 2755362 : y = mului(u,y);
1708 2755356 : set_avma(av); addzi(x, y);
1709 : }
1710 :
1711 : /************************** PROVED version (dpe) *************************/
1712 :
1713 : /* Babai's Nearest Plane algorithm (iterative).
1714 : * Size-reduces b_kappa using mu_{i,j} and r_{i,j} for j<=i <kappa
1715 : * Update B[,kappa]; compute mu_{kappa,j}, r_{kappa,j} for j<=kappa and s[kappa]
1716 : * mu, r, s updated in place (affrr). Return 1 on failure, else 0. */
1717 : static int
1718 4585965 : Babai_dpe(pari_sp av, long kappa, GEN *pG, GEN *pB, GEN *pU, dpe_t **mu, dpe_t **r, dpe_t *s,
1719 : long a, long zeros, long maxG, dpe_t *eta)
1720 : {
1721 4585965 : GEN G = *pG, B = *pB, U = *pU, ztmp;
1722 4585965 : long k, d, n, aa = a > zeros? a: zeros+1;
1723 4585965 : long emaxmu = EX0, emax2mu = EX0;
1724 : /* N.B: we set d = 0 (resp. n = 0) to avoid updating U (resp. B) */
1725 4585965 : d = U? lg(U)-1: 0;
1726 4585965 : n = B? nbrows(B): 0;
1727 522162 : for (;;) {
1728 5108136 : int go_on = 0;
1729 5108136 : long i, j, emax3mu = emax2mu;
1730 :
1731 5108136 : if (gc_needed(av,2))
1732 : {
1733 0 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[1], a=%ld", aa);
1734 0 : gc_lll(av,3,&G,&B,&U);
1735 : }
1736 : /* Step2: compute the GSO for stage kappa */
1737 5108133 : emax2mu = emaxmu; emaxmu = EX0;
1738 19062820 : for (j=aa; j<kappa; j++)
1739 : {
1740 : dpe_t g;
1741 13954694 : affidpe(gmael(G,kappa,j), &g);
1742 52323843 : for (k = zeros+1; k < j; k++)
1743 38369241 : dpe_submulz(&g, Dmael(mu,j,k), Dmael(r,kappa,k), &g);
1744 13954602 : affdpe(&g, Dmael(r,kappa,j));
1745 13954656 : dpe_divz(Dmael(r,kappa,j), Dmael(r,j,j), Dmael(mu,kappa,j));
1746 13954659 : emaxmu = maxss(emaxmu, Dmael(mu,kappa,j)->e);
1747 : }
1748 5108126 : if (emax3mu != EX0 && emax3mu <= emax2mu + 5) /* precision too low */
1749 0 : { *pG = G; *pB = B; *pU = U; return 1; }
1750 :
1751 18981468 : for (j=kappa-1; j>zeros; j--)
1752 14395510 : if (dpe_abscmp(Dmael(mu,kappa,j), eta) > 0) { go_on = 1; break; }
1753 :
1754 : /* Step3--5: compute the X_j's */
1755 5108128 : if (go_on)
1756 3030748 : for (j=kappa-1; j>zeros; j--)
1757 : {
1758 : pari_sp btop;
1759 2508580 : dpe_t *tmp = Dmael(mu,kappa,j);
1760 2508580 : if (tmp->e < 0) continue; /* (essentially) size-reduced */
1761 :
1762 1387784 : if (gc_needed(av,2))
1763 : {
1764 0 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[2], a=%ld, j=%ld", aa,j);
1765 0 : gc_lll(av,3,&G,&B,&U);
1766 : }
1767 : /* we consider separately the case |X| = 1 */
1768 1387784 : if (dpe_abssmall(tmp))
1769 : {
1770 920626 : if (tmp->d > 0) { /* in this case, X = 1 */
1771 2057648 : for (k=zeros+1; k<j; k++)
1772 1596103 : dpe_subz(Dmael(mu,kappa,k), Dmael(mu,j,k), Dmael(mu,kappa,k));
1773 3016918 : for (i=1; i<=n; i++)
1774 2555373 : subzi(&gmael(B,kappa,i), gmael(B,j,i));
1775 6970707 : for (i=1; i<=d; i++)
1776 6509169 : subzi(&gmael(U,kappa,i), gmael(U,j,i));
1777 461538 : btop = avma;
1778 461538 : ztmp = subii(gmael(G,j,j), shifti(gmael(G,kappa,j), 1));
1779 461545 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
1780 461545 : affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
1781 2859391 : for (i=1; i<=j; i++)
1782 2397847 : subzi(&gmael(G,kappa,i), gmael(G,j,i));
1783 2191653 : for (i=j+1; i<kappa; i++)
1784 1730110 : subzi(&gmael(G,kappa,i), gmael(G,i,j));
1785 2360874 : for (i=kappa+1; i<=maxG; i++)
1786 1899333 : subzi(&gmael(G,i,kappa), gmael(G,i,j));
1787 : } else { /* otherwise X = -1 */
1788 2035772 : for (k=zeros+1; k<j; k++)
1789 1576688 : dpe_addz(Dmael(mu,kappa,k), Dmael(mu,j,k), Dmael(mu,kappa,k));
1790 3011876 : for (i=1; i<=n; i++)
1791 2552792 : addzi(&gmael(B,kappa,i),gmael(B,j,i));
1792 6843675 : for (i=1; i<=d; i++)
1793 6384607 : addzi(&gmael(U,kappa,i),gmael(U,j,i));
1794 459068 : btop = avma;
1795 459068 : ztmp = addii(gmael(G,j,j), shifti(gmael(G,kappa,j), 1));
1796 459080 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
1797 459080 : affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
1798 2770241 : for (i=1; i<=j; i++)
1799 2311158 : addzi(&gmael(G,kappa,i), gmael(G,j,i));
1800 2195487 : for (i=j+1; i<kappa; i++)
1801 1736406 : addzi(&gmael(G,kappa,i), gmael(G,i,j));
1802 2313488 : for (i=kappa+1; i<=maxG; i++)
1803 1854406 : addzi(&gmael(G,i,kappa), gmael(G,i,j));
1804 : }
1805 920623 : continue;
1806 : }
1807 : /* we have |X| >= 2 */
1808 467158 : if (tmp->e < BITS_IN_LONG-1)
1809 : {
1810 448160 : if (tmp->d > 0)
1811 : {
1812 247183 : ulong xx = (ulong) pari_rint(ldexp(tmp->d, tmp->e)); /* X fits in an ulong */
1813 659087 : for (k=zeros+1; k<j; k++)
1814 411904 : dpe_submuluz(Dmael(mu,kappa,k), Dmael(mu,j,k), xx, Dmael(mu,kappa,k));
1815 722230 : for (i=1; i<=n; i++)
1816 475047 : submulziu(&gmael(B,kappa,i), gmael(B,j,i), xx);
1817 3107316 : for (i=1; i<=d; i++)
1818 2860129 : submulziu(&gmael(U,kappa,i), gmael(U,j,i), xx);
1819 247187 : btop = avma;
1820 247187 : ztmp = submuliu2n(mulii(gmael(G,j,j), sqru(xx)), gmael(G,kappa,j), xx, 1);
1821 247182 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
1822 247183 : affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
1823 1313763 : for (i=1; i<=j; i++)
1824 1066583 : submulziu(&gmael(G,kappa,i), gmael(G,j,i), xx);
1825 807865 : for (i=j+1; i<kappa; i++)
1826 560685 : submulziu(&gmael(G,kappa,i), gmael(G,i,j), xx);
1827 1204790 : for (i=kappa+1; i<=maxG; i++)
1828 957611 : submulziu(&gmael(G,i,kappa), gmael(G,i,j), xx);
1829 : }
1830 : else
1831 : {
1832 200977 : ulong xx = (ulong) pari_rint(ldexp(-tmp->d, tmp->e)); /* X fits in an ulong */
1833 543626 : for (k=zeros+1; k<j; k++)
1834 342647 : dpe_addmuluz(Dmael(mu,kappa,k), Dmael(mu,j,k), xx, Dmael(mu,kappa,k));
1835 687415 : for (i=1; i<=n; i++)
1836 486436 : addmulziu(&gmael(B,kappa,i), gmael(B,j,i), xx);
1837 2359203 : for (i=1; i<=d; i++)
1838 2158225 : addmulziu(&gmael(U,kappa,i), gmael(U,j,i), xx);
1839 200978 : btop = avma;
1840 200978 : ztmp = addmuliu2n(mulii(gmael(G,j,j), sqru(xx)), gmael(G,kappa,j), xx, 1);
1841 200977 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
1842 200977 : affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
1843 989994 : for (i=1; i<=j; i++)
1844 789017 : addmulziu(&gmael(G,kappa,i), gmael(G,j,i), xx);
1845 662168 : for (i=j+1; i<kappa; i++)
1846 461190 : addmulziu(&gmael(G,kappa,i), gmael(G,i,j), xx);
1847 883029 : for (i=kappa+1; i<=maxG; i++)
1848 682051 : addmulziu(&gmael(G,i,kappa), gmael(G,i,j), xx);
1849 : }
1850 : }
1851 : else
1852 : {
1853 18998 : long e = tmp->e - BITS_IN_LONG + 1;
1854 18998 : if (tmp->d > 0)
1855 : {
1856 9396 : ulong xx = (ulong) pari_rint(ldexp(tmp->d, BITS_IN_LONG - 1));
1857 31333 : for (k=zeros+1; k<j; k++)
1858 : {
1859 : dpe_t x;
1860 21937 : dpe_muluz(Dmael(mu,j,k), xx, &x);
1861 21937 : x.e += e;
1862 21937 : dpe_subz(Dmael(mu,kappa,k), &x, Dmael(mu,kappa,k));
1863 : }
1864 124323 : for (i=1; i<=n; i++)
1865 114927 : submulzu2n(&gmael(B,kappa,i), gmael(B,j,i), xx, e);
1866 87036 : for (i=1; i<=d; i++)
1867 77640 : submulzu2n(&gmael(U,kappa,i), gmael(U,j,i), xx, e);
1868 9396 : btop = avma;
1869 9396 : ztmp = submuliu2n(mulshift(gmael(G,j,j), sqru(xx), 2*e),
1870 9396 : gmael(G,kappa,j), xx, e+1);
1871 9396 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
1872 9396 : affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
1873 40927 : for (i=1; i<=j; i++)
1874 31531 : submulzu2n(&gmael(G,kappa,i), gmael(G,j,i), xx, e);
1875 47343 : for ( ; i<kappa; i++)
1876 37947 : submulzu2n(&gmael(G,kappa,i), gmael(G,i,j), xx, e);
1877 10308 : for (i=kappa+1; i<=maxG; i++)
1878 912 : submulzu2n(&gmael(G,i,kappa), gmael(G,i,j), xx, e);
1879 : } else
1880 : {
1881 9602 : ulong xx = (ulong) pari_rint(ldexp(-tmp->d, BITS_IN_LONG - 1));
1882 32153 : for (k=zeros+1; k<j; k++)
1883 : {
1884 : dpe_t x;
1885 22547 : dpe_muluz(Dmael(mu,j,k), xx, &x);
1886 22547 : x.e += e;
1887 22547 : dpe_addz(Dmael(mu,kappa,k), &x, Dmael(mu,kappa,k));
1888 : }
1889 128490 : for (i=1; i<=n; i++)
1890 118884 : addmulzu2n(&gmael(B,kappa,i), gmael(B,j,i), xx, e);
1891 89277 : for (i=1; i<=d; i++)
1892 79671 : addmulzu2n(&gmael(U,kappa,i), gmael(U,j,i), xx, e);
1893 9606 : btop = avma;
1894 9606 : ztmp = addmuliu2n(mulshift(gmael(G,j,j), sqru(xx), 2*e),
1895 9606 : gmael(G,kappa,j), xx, e+1);
1896 9606 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
1897 9606 : affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
1898 41957 : for (i=1; i<=j; i++)
1899 32351 : addmulzu2n(&gmael(G,kappa,i), gmael(G,j,i), xx, e);
1900 48921 : for ( ; i<kappa; i++)
1901 39315 : addmulzu2n(&gmael(G,kappa,i), gmael(G,i,j), xx, e);
1902 10353 : for (i=kappa+1; i<=maxG; i++)
1903 747 : addmulzu2n(&gmael(G,i,kappa), gmael(G,i,j), xx, e);
1904 : }
1905 : }
1906 : }
1907 5108126 : if (!go_on) break; /* Anything happened? */
1908 522162 : aa = zeros+1;
1909 : }
1910 :
1911 4585964 : affidpe(gmael(G,kappa,kappa), Del(s,zeros+1));
1912 : /* the last s[kappa-1]=r[kappa][kappa] is computed only if kappa increases */
1913 13468539 : for (k=zeros+1; k<=kappa-2; k++)
1914 8882602 : dpe_submulz(Del(s,k), Dmael(mu,kappa,k), Dmael(r,kappa,k), Del(s,k+1));
1915 4585937 : *pG = G; *pB = B; *pU = U; return 0;
1916 : }
1917 :
1918 : /* G integral Gram matrix, LLL-reduces (G,B,U) in place [apply base change
1919 : * transforms to B and U]. If (keepfirst), never swap with first vector.
1920 : * If G = NULL, we compute the Gram matrix incrementally.
1921 : * Return -1 on failure, else zeros = dim Kernel (>= 0) */
1922 : static long
1923 1567955 : fplll_dpe(GEN *pG, GEN *pB, GEN *pU, GEN *pr, double DELTA, double ETA,
1924 : long keepfirst)
1925 : {
1926 : pari_sp av;
1927 1567955 : GEN Gtmp, alpha, G = *pG, B = *pB, U = *pU;
1928 1567955 : long d, maxG, kappa, kappa2, i, j, zeros, kappamax, incgram = !G, cnt = 0;
1929 : dpe_t delta, eta, **mu, **r, *s;
1930 1567955 : affdbldpe(DELTA,&delta);
1931 1567954 : affdbldpe(ETA,&eta);
1932 :
1933 1567954 : if (incgram)
1934 : { /* incremental Gram matrix */
1935 1507678 : maxG = 2; d = lg(B)-1;
1936 1507678 : G = zeromatcopy(d, d);
1937 : }
1938 : else
1939 60276 : maxG = d = lg(G)-1;
1940 :
1941 1567951 : mu = cget_dpemat(d+1);
1942 1567949 : r = cget_dpemat(d+1);
1943 1567948 : s = cget_dpevec(d+1);
1944 7304865 : for (j = 1; j <= d; j++)
1945 : {
1946 5736918 : mu[j]= cget_dpevec(d+1);
1947 5736914 : r[j] = cget_dpevec(d+1);
1948 : }
1949 1567947 : Gtmp = cgetg(d+1, t_VEC);
1950 1567946 : alpha = cgetg(d+1, t_VECSMALL);
1951 1567944 : av = avma;
1952 :
1953 : /* Step2: Initializing the main loop */
1954 1567944 : kappamax = 1;
1955 1567944 : i = 1;
1956 : do {
1957 1944935 : if (incgram) gmael(G,i,i) = ZV_dotsquare(gel(B,i));
1958 1944877 : affidpe(gmael(G,i,i), Dmael(r,i,i));
1959 1944895 : } while (!signe(gmael(G,i,i)) && ++i <= d);
1960 1567904 : zeros = i-1; /* all basis vectors b_i with i <= zeros are zero vectors */
1961 1567904 : kappa = i;
1962 6927706 : for (i=zeros+1; i<=d; i++) alpha[i]=1;
1963 :
1964 6153880 : while (++kappa <= d)
1965 : {
1966 4585945 : if (kappa > kappamax)
1967 : {
1968 3791909 : if (DEBUGLEVEL>=4) err_printf("K%ld ",kappa);
1969 3791909 : kappamax = kappa;
1970 3791909 : if (incgram)
1971 : {
1972 15911189 : for (i=zeros+1; i<=kappa; i++)
1973 12318798 : gmael(G,kappa,i) = ZV_dotproduct(gel(B,kappa), gel(B,i));
1974 3592391 : maxG = kappamax;
1975 : }
1976 : }
1977 : /* Step3: Call to the Babai algorithm, mu,r,s updated in place */
1978 4585744 : if (Babai_dpe(av, kappa, &G,&B,&U, mu,r,s, alpha[kappa], zeros, maxG, &eta))
1979 0 : { *pG = incgram? NULL: G; *pB = B; *pU = U; return -1; }
1980 9072616 : if ((keepfirst && kappa == 2) ||
1981 4486653 : dpe_cmpmul(Dmael(r,kappa-1,kappa-1), &delta, Del(s,kappa-1)) <= 0)
1982 : { /* Step4: Success of Lovasz's condition */
1983 4260268 : alpha[kappa] = kappa;
1984 4260268 : dpe_submulz(Del(s,kappa-1), Dmael(mu,kappa,kappa-1), Dmael(r,kappa,kappa-1), Dmael(r,kappa,kappa));
1985 4260274 : continue;
1986 : }
1987 : /* Step5: Find the right insertion index kappa, kappa2 = initial kappa */
1988 325695 : if (DEBUGLEVEL>=4 && kappa==kappamax && Del(s,kappa-1)->d)
1989 0 : if (++cnt > 20) { cnt = 0; err_printf("(%ld) ", Del(s,1)->e-1); }
1990 325695 : kappa2 = kappa;
1991 : do {
1992 829199 : kappa--;
1993 829199 : if (kappa < zeros+2 + (keepfirst ? 1: 0)) break;
1994 707517 : } while (dpe_cmpmul(Dmael(r,kappa-1,kappa-1), &delta, Del(s,kappa-1)) >= 0);
1995 325695 : update_alpha(alpha, kappa, kappa2, kappamax);
1996 :
1997 : /* Step6: Update the mu's and r's */
1998 325695 : dperotate(mu, kappa2, kappa);
1999 325695 : dperotate(r, kappa2, kappa);
2000 325695 : affdpe(Del(s,kappa), Dmael(r,kappa,kappa));
2001 :
2002 : /* Step7: Update G, B, U */
2003 325695 : if (U) rotate(U, kappa2, kappa);
2004 325695 : if (B) rotate(B, kappa2, kappa);
2005 325695 : rotateG(G,kappa2,kappa, maxG, Gtmp);
2006 :
2007 : /* Step8: Prepare the next loop iteration */
2008 325695 : if (kappa == zeros+1 && !signe(gmael(G,kappa,kappa)))
2009 : {
2010 35161 : zeros++; kappa++;
2011 35161 : affidpe(gmael(G,kappa,kappa), Dmael(r,kappa,kappa));
2012 : }
2013 : }
2014 1567935 : if (pr) *pr = dpeM_diagonal_shallow(r,d);
2015 1567935 : *pG = G; *pB = B; *pU = U; return zeros; /* success */
2016 : }
2017 :
2018 :
2019 : /************************** PROVED version (t_INT) *************************/
2020 :
2021 : /* Babai's Nearest Plane algorithm (iterative).
2022 : * Size-reduces b_kappa using mu_{i,j} and r_{i,j} for j<=i <kappa
2023 : * Update B[,kappa]; compute mu_{kappa,j}, r_{kappa,j} for j<=kappa and s[kappa]
2024 : * mu, r, s updated in place (affrr). Return 1 on failure, else 0. */
2025 : static int
2026 0 : Babai(pari_sp av, long kappa, GEN *pG, GEN *pB, GEN *pU, GEN mu, GEN r, GEN s,
2027 : long a, long zeros, long maxG, GEN eta, long prec)
2028 : {
2029 0 : GEN G = *pG, B = *pB, U = *pU, ztmp;
2030 0 : long k, aa = a > zeros? a: zeros+1;
2031 0 : const long n = B? nbrows(B): 0, d = U ? lg(U)-1: 0, bit = prec2nbits(prec);
2032 0 : long emaxmu = EX0, emax2mu = EX0;
2033 : /* N.B: we set d = 0 (resp. n = 0) to avoid updating U (resp. B) */
2034 :
2035 0 : for (;;) {
2036 0 : int go_on = 0;
2037 0 : long i, j, emax3mu = emax2mu;
2038 :
2039 0 : if (gc_needed(av,2))
2040 : {
2041 0 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[1], a=%ld", aa);
2042 0 : gc_lll(av,3,&G,&B,&U);
2043 : }
2044 : /* Step2: compute the GSO for stage kappa */
2045 0 : emax2mu = emaxmu; emaxmu = EX0;
2046 0 : for (j=aa; j<kappa; j++)
2047 : {
2048 0 : pari_sp btop = avma;
2049 0 : GEN g = gmael(G,kappa,j);
2050 0 : for (k = zeros+1; k < j; k++)
2051 0 : g = mpsub(g, mulrr(gmael(mu,j,k), gmael(r,kappa,k)));
2052 0 : mpaff(g, gmael(r,kappa,j));
2053 0 : affrr(divrr(gmael(r,kappa,j), gmael(r,j,j)), gmael(mu,kappa,j));
2054 0 : emaxmu = maxss(emaxmu, expo(gmael(mu,kappa,j)));
2055 0 : set_avma(btop);
2056 : }
2057 0 : if (emax3mu != EX0 && emax3mu <= emax2mu + 5) /* precision too low */
2058 0 : { *pG = G; *pB = B; *pU = U; return 1; }
2059 :
2060 0 : for (j=kappa-1; j>zeros; j--)
2061 0 : if (abscmprr(gmael(mu,kappa,j), eta) > 0) { go_on = 1; break; }
2062 :
2063 : /* Step3--5: compute the X_j's */
2064 0 : if (go_on)
2065 0 : for (j=kappa-1; j>zeros; j--)
2066 : {
2067 : pari_sp btop;
2068 0 : GEN tmp = gmael(mu,kappa,j);
2069 0 : if (absrsmall(tmp)) continue; /* size-reduced */
2070 :
2071 0 : if (gc_needed(av,2))
2072 : {
2073 0 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[2], a=%ld, j=%ld", aa,j);
2074 0 : gc_lll(av,3,&G,&B,&U);
2075 : }
2076 0 : btop = avma;
2077 : /* we consider separately the case |X| = 1 */
2078 0 : if (absrsmall2(tmp))
2079 : {
2080 0 : if (signe(tmp) > 0) { /* in this case, X = 1 */
2081 0 : for (k=zeros+1; k<j; k++)
2082 0 : affrr(subrr(gmael(mu,kappa,k), gmael(mu,j,k)), gmael(mu,kappa,k));
2083 0 : set_avma(btop);
2084 0 : for (i=1; i<=n; i++)
2085 0 : gmael(B,kappa,i) = subii(gmael(B,kappa,i), gmael(B,j,i));
2086 0 : for (i=1; i<=d; i++)
2087 0 : gmael(U,kappa,i) = subii(gmael(U,kappa,i), gmael(U,j,i));
2088 0 : btop = avma;
2089 0 : ztmp = subii(gmael(G,j,j), shifti(gmael(G,kappa,j), 1));
2090 0 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
2091 0 : gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
2092 0 : for (i=1; i<=j; i++)
2093 0 : gmael(G,kappa,i) = subii(gmael(G,kappa,i), gmael(G,j,i));
2094 0 : for (i=j+1; i<kappa; i++)
2095 0 : gmael(G,kappa,i) = subii(gmael(G,kappa,i), gmael(G,i,j));
2096 0 : for (i=kappa+1; i<=maxG; i++)
2097 0 : gmael(G,i,kappa) = subii(gmael(G,i,kappa), gmael(G,i,j));
2098 : } else { /* otherwise X = -1 */
2099 0 : for (k=zeros+1; k<j; k++)
2100 0 : affrr(addrr(gmael(mu,kappa,k), gmael(mu,j,k)), gmael(mu,kappa,k));
2101 0 : set_avma(btop);
2102 0 : for (i=1; i<=n; i++)
2103 0 : gmael(B,kappa,i) = addii(gmael(B,kappa,i),gmael(B,j,i));
2104 0 : for (i=1; i<=d; i++)
2105 0 : gmael(U,kappa,i) = addii(gmael(U,kappa,i),gmael(U,j,i));
2106 0 : btop = avma;
2107 0 : ztmp = addii(gmael(G,j,j), shifti(gmael(G,kappa,j), 1));
2108 0 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
2109 0 : gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
2110 0 : for (i=1; i<=j; i++)
2111 0 : gmael(G,kappa,i) = addii(gmael(G,kappa,i), gmael(G,j,i));
2112 0 : for (i=j+1; i<kappa; i++)
2113 0 : gmael(G,kappa,i) = addii(gmael(G,kappa,i), gmael(G,i,j));
2114 0 : for (i=kappa+1; i<=maxG; i++)
2115 0 : gmael(G,i,kappa) = addii(gmael(G,i,kappa), gmael(G,i,j));
2116 : }
2117 0 : continue;
2118 : }
2119 : /* we have |X| >= 2 */
2120 0 : if (expo(tmp) < BITS_IN_LONG)
2121 : {
2122 0 : ulong xx = roundr_safe(tmp)[2]; /* X fits in an ulong */
2123 0 : if (signe(tmp) > 0) /* = xx */
2124 : {
2125 0 : for (k=zeros+1; k<j; k++)
2126 0 : affrr(subrr(gmael(mu,kappa,k), mulur(xx, gmael(mu,j,k))),
2127 0 : gmael(mu,kappa,k));
2128 0 : set_avma(btop);
2129 0 : for (i=1; i<=n; i++)
2130 0 : gmael(B,kappa,i) = submuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
2131 0 : for (i=1; i<=d; i++)
2132 0 : gmael(U,kappa,i) = submuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
2133 0 : btop = avma;
2134 0 : ztmp = submuliu2n(mulii(gmael(G,j,j), sqru(xx)), gmael(G,kappa,j), xx, 1);
2135 0 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
2136 0 : gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
2137 0 : for (i=1; i<=j; i++)
2138 0 : gmael(G,kappa,i) = submuliu_inplace(gmael(G,kappa,i), gmael(G,j,i), xx);
2139 0 : for (i=j+1; i<kappa; i++)
2140 0 : gmael(G,kappa,i) = submuliu_inplace(gmael(G,kappa,i), gmael(G,i,j), xx);
2141 0 : for (i=kappa+1; i<=maxG; i++)
2142 0 : gmael(G,i,kappa) = submuliu_inplace(gmael(G,i,kappa), gmael(G,i,j), xx);
2143 : }
2144 : else /* = -xx */
2145 : {
2146 0 : for (k=zeros+1; k<j; k++)
2147 0 : affrr(addrr(gmael(mu,kappa,k), mulur(xx, gmael(mu,j,k))),
2148 0 : gmael(mu,kappa,k));
2149 0 : set_avma(btop);
2150 0 : for (i=1; i<=n; i++)
2151 0 : gmael(B,kappa,i) = addmuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
2152 0 : for (i=1; i<=d; i++)
2153 0 : gmael(U,kappa,i) = addmuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
2154 0 : btop = avma;
2155 0 : ztmp = addmuliu2n(mulii(gmael(G,j,j), sqru(xx)), gmael(G,kappa,j), xx, 1);
2156 0 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
2157 0 : gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
2158 0 : for (i=1; i<=j; i++)
2159 0 : gmael(G,kappa,i) = addmuliu_inplace(gmael(G,kappa,i), gmael(G,j,i), xx);
2160 0 : for (i=j+1; i<kappa; i++)
2161 0 : gmael(G,kappa,i) = addmuliu_inplace(gmael(G,kappa,i), gmael(G,i,j), xx);
2162 0 : for (i=kappa+1; i<=maxG; i++)
2163 0 : gmael(G,i,kappa) = addmuliu_inplace(gmael(G,i,kappa), gmael(G,i,j), xx);
2164 : }
2165 : }
2166 : else
2167 : {
2168 : long e;
2169 0 : GEN X = truncexpo(tmp, bit, &e); /* tmp ~ X * 2^e */
2170 0 : btop = avma;
2171 0 : for (k=zeros+1; k<j; k++)
2172 : {
2173 0 : GEN x = mulir(X, gmael(mu,j,k));
2174 0 : if (e) shiftr_inplace(x, e);
2175 0 : affrr(subrr(gmael(mu,kappa,k), x), gmael(mu,kappa,k));
2176 : }
2177 0 : set_avma(btop);
2178 0 : for (i=1; i<=n; i++)
2179 0 : gmael(B,kappa,i) = submulshift(gmael(B,kappa,i), gmael(B,j,i), X, e);
2180 0 : for (i=1; i<=d; i++)
2181 0 : gmael(U,kappa,i) = submulshift(gmael(U,kappa,i), gmael(U,j,i), X, e);
2182 0 : btop = avma;
2183 0 : ztmp = submulshift(mulshift(gmael(G,j,j), sqri(X), 2*e),
2184 0 : gmael(G,kappa,j), X, e+1);
2185 0 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
2186 0 : gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
2187 0 : for (i=1; i<=j; i++)
2188 0 : gmael(G,kappa,i) = submulshift(gmael(G,kappa,i), gmael(G,j,i), X, e);
2189 0 : for ( ; i<kappa; i++)
2190 0 : gmael(G,kappa,i) = submulshift(gmael(G,kappa,i), gmael(G,i,j), X, e);
2191 0 : for (i=kappa+1; i<=maxG; i++)
2192 0 : gmael(G,i,kappa) = submulshift(gmael(G,i,kappa), gmael(G,i,j), X, e);
2193 : }
2194 : }
2195 0 : if (!go_on) break; /* Anything happened? */
2196 0 : aa = zeros+1;
2197 : }
2198 :
2199 0 : affir(gmael(G,kappa,kappa), gel(s,zeros+1));
2200 : /* the last s[kappa-1]=r[kappa][kappa] is computed only if kappa increases */
2201 0 : av = avma;
2202 0 : for (k=zeros+1; k<=kappa-2; k++)
2203 0 : affrr(subrr(gel(s,k), mulrr(gmael(mu,kappa,k), gmael(r,kappa,k))),
2204 0 : gel(s,k+1));
2205 0 : *pG = G; *pB = B; *pU = U; return gc_bool(av, 0);
2206 : }
2207 :
2208 : /* G integral Gram matrix, LLL-reduces (G,B,U) in place [apply base change
2209 : * transforms to B and U]. If (keepfirst), never swap with first vector.
2210 : * If G = NULL, we compute the Gram matrix incrementally.
2211 : * Return -1 on failure, else zeros = dim Kernel (>= 0) */
2212 : static long
2213 0 : fplll(GEN *pG, GEN *pB, GEN *pU, GEN *pr, double DELTA, double ETA,
2214 : long keepfirst, long prec)
2215 : {
2216 : pari_sp av, av2;
2217 0 : GEN mu, r, s, tmp, Gtmp, alpha, G = *pG, B = *pB, U = *pU;
2218 0 : GEN delta = dbltor(DELTA), eta = dbltor(ETA);
2219 0 : long d, maxG, kappa, kappa2, i, j, zeros, kappamax, incgram = !G, cnt = 0;
2220 :
2221 0 : if (incgram)
2222 : { /* incremental Gram matrix */
2223 0 : maxG = 2; d = lg(B)-1;
2224 0 : G = zeromatcopy(d, d);
2225 : }
2226 : else
2227 0 : maxG = d = lg(G)-1;
2228 :
2229 0 : mu = cgetg(d+1, t_MAT);
2230 0 : r = cgetg(d+1, t_MAT);
2231 0 : s = cgetg(d+1, t_VEC);
2232 0 : for (j = 1; j <= d; j++)
2233 : {
2234 0 : GEN M = cgetg(d+1, t_COL), R = cgetg(d+1, t_COL);
2235 0 : gel(mu,j)= M;
2236 0 : gel(r,j) = R;
2237 0 : gel(s,j) = cgetr(prec);
2238 0 : for (i = 1; i <= d; i++)
2239 : {
2240 0 : gel(R,i) = cgetr(prec);
2241 0 : gel(M,i) = cgetr(prec);
2242 : }
2243 : }
2244 0 : Gtmp = cgetg(d+1, t_VEC);
2245 0 : alpha = cgetg(d+1, t_VECSMALL);
2246 0 : av = avma;
2247 :
2248 : /* Step2: Initializing the main loop */
2249 0 : kappamax = 1;
2250 0 : i = 1;
2251 : do {
2252 0 : if (incgram) gmael(G,i,i) = ZV_dotsquare(gel(B,i));
2253 0 : affir(gmael(G,i,i), gmael(r,i,i));
2254 0 : } while (!signe(gmael(G,i,i)) && ++i <= d);
2255 0 : zeros = i-1; /* all basis vectors b_i with i <= zeros are zero vectors */
2256 0 : kappa = i;
2257 0 : for (i=zeros+1; i<=d; i++) alpha[i]=1;
2258 :
2259 0 : while (++kappa <= d)
2260 : {
2261 0 : if (kappa > kappamax)
2262 : {
2263 0 : if (DEBUGLEVEL>=4) err_printf("K%ld ",kappa);
2264 0 : kappamax = kappa;
2265 0 : if (incgram)
2266 : {
2267 0 : for (i=zeros+1; i<=kappa; i++)
2268 0 : gmael(G,kappa,i) = ZV_dotproduct(gel(B,kappa), gel(B,i));
2269 0 : maxG = kappamax;
2270 : }
2271 : }
2272 : /* Step3: Call to the Babai algorithm, mu,r,s updated in place */
2273 0 : if (Babai(av, kappa, &G,&B,&U, mu,r,s, alpha[kappa], zeros, maxG, eta, prec))
2274 0 : { *pG = incgram? NULL: G; *pB = B; *pU = U; return -1; }
2275 0 : av2 = avma;
2276 0 : if ((keepfirst && kappa == 2) ||
2277 0 : cmprr(mulrr(gmael(r,kappa-1,kappa-1), delta), gel(s,kappa-1)) <= 0)
2278 : { /* Step4: Success of Lovasz's condition */
2279 0 : alpha[kappa] = kappa;
2280 0 : tmp = mulrr(gmael(mu,kappa,kappa-1), gmael(r,kappa,kappa-1));
2281 0 : affrr(subrr(gel(s,kappa-1), tmp), gmael(r,kappa,kappa));
2282 0 : set_avma(av2); continue;
2283 : }
2284 : /* Step5: Find the right insertion index kappa, kappa2 = initial kappa */
2285 0 : if (DEBUGLEVEL>=4 && kappa==kappamax && signe(gel(s,kappa-1)))
2286 0 : if (++cnt > 20) { cnt = 0; err_printf("(%ld) ", expo(gel(s,1))); }
2287 0 : kappa2 = kappa;
2288 : do {
2289 0 : kappa--;
2290 0 : if (kappa < zeros+2 + (keepfirst ? 1: 0)) break;
2291 0 : tmp = mulrr(gmael(r,kappa-1,kappa-1), delta);
2292 0 : } while (cmprr(gel(s,kappa-1), tmp) <= 0);
2293 0 : set_avma(av2);
2294 0 : update_alpha(alpha, kappa, kappa2, kappamax);
2295 :
2296 : /* Step6: Update the mu's and r's */
2297 0 : rotate(mu, kappa2, kappa);
2298 0 : rotate(r, kappa2, kappa);
2299 0 : affrr(gel(s,kappa), gmael(r,kappa,kappa));
2300 :
2301 : /* Step7: Update G, B, U */
2302 0 : if (U) rotate(U, kappa2, kappa);
2303 0 : if (B) rotate(B, kappa2, kappa);
2304 0 : rotateG(G,kappa2,kappa, maxG, Gtmp);
2305 :
2306 : /* Step8: Prepare the next loop iteration */
2307 0 : if (kappa == zeros+1 && !signe(gmael(G,kappa,kappa)))
2308 : {
2309 0 : zeros++; kappa++;
2310 0 : affir(gmael(G,kappa,kappa), gmael(r,kappa,kappa));
2311 : }
2312 : }
2313 0 : if (pr) *pr = RgM_diagonal_shallow(r);
2314 0 : *pG = G; *pB = B; *pU = U; return zeros; /* success */
2315 : }
2316 :
2317 : /* do not support LLL_KER, LLL_ALL, LLL_KEEP_FIRST */
2318 : static GEN
2319 4701008 : ZM2_lll_norms(GEN x, long flag, GEN *pN)
2320 : {
2321 : GEN a,b,c,d;
2322 : GEN G, U;
2323 4701008 : if (flag & LLL_GRAM)
2324 7355 : G = x;
2325 : else
2326 4693653 : G = gram_matrix(x);
2327 4700996 : a = gcoeff(G,1,1); b = shifti(gcoeff(G,1,2),1); c = gcoeff(G,2,2);
2328 4700953 : d = qfb_disc3(a,b,c);
2329 4700955 : if (signe(d)>=0) return NULL;
2330 4700570 : G = redimagsl2(mkqfb(a,b,c,d),&U);
2331 4700624 : if (pN) (void) RgM_gram_schmidt(G, pN);
2332 4700623 : if (flag & LLL_INPLACE) return ZM2_mul(x,U);
2333 4700623 : return U;
2334 : }
2335 :
2336 : static void
2337 617331 : fplll_flatter(GEN *pG, GEN *pB, GEN *pU, long rank, long flag)
2338 : {
2339 617331 : if (!*pG)
2340 : {
2341 616370 : GEN T = ZM_flatter_rank(*pB, rank, flag);
2342 616371 : if (T)
2343 : {
2344 327008 : if (*pU)
2345 : {
2346 313007 : *pU = ZM_mul(*pU, T);
2347 313007 : *pB = ZM_mul(*pB, T);
2348 : }
2349 14001 : else *pB = T;
2350 : }
2351 : }
2352 : else
2353 : {
2354 961 : GEN T, G = *pG;
2355 961 : long i, j, l = lg(G);
2356 7207 : for (i = 1; i < l; i++)
2357 43383 : for(j = 1; j < i; j++) gmael(G,j,i) = gmael(G,i,j);
2358 961 : T = ZM_flattergram_rank(G, rank, flag);
2359 961 : if (T)
2360 : {
2361 961 : if (*pU) *pU = ZM_mul(*pU, T);
2362 961 : *pG = qf_ZM_apply(*pG, T);
2363 : }
2364 : }
2365 617332 : }
2366 :
2367 : static GEN
2368 1082658 : get_gramschmidt(GEN M, long rank)
2369 : {
2370 : GEN B, Q, L;
2371 1082658 : long r = lg(M)-1, bitprec = 3*r + 30;
2372 1082658 : long prec = nbits2prec64(bitprec);
2373 1082658 : if (rank < r) M = vconcat(gshift(M,1), matid(r));
2374 1082659 : if (!QR_init(RgM_gtofp(M, prec), &B, &Q, &L, prec) || !gsisinv(L)) return NULL;
2375 467360 : return L;
2376 : }
2377 :
2378 : static GEN
2379 44260 : get_gaussred(GEN M, long rank)
2380 : {
2381 44260 : pari_sp ltop = avma;
2382 44260 : long r = lg(M)-1, bitprec = 3*r + 30, prec = nbits2prec64(bitprec);
2383 : GEN R;
2384 44260 : if (rank < r) M = RgM_Rg_add(gshift(M, 1), gen_1);
2385 44260 : R = RgM_Cholesky(RgM_gtofp(M, prec), prec);
2386 44260 : if (!R) return NULL;
2387 43299 : return gerepilecopy(ltop, R);
2388 : }
2389 :
2390 : /* Assume x a ZM, if pN != NULL, set it to Gram-Schmidt (squared) norms
2391 : * The following modes are supported:
2392 : * - flag & LLL_INPLACE: x a lattice basis, return x*U
2393 : * - flag & LLL_GRAM: x a Gram matrix / else x a lattice basis; return
2394 : * LLL base change matrix U [LLL_IM]
2395 : * kernel basis [LLL_KER, nonreduced]
2396 : * both [LLL_ALL] */
2397 : GEN
2398 6953999 : ZM_lll_norms(GEN x, double DELTA, long flag, GEN *pN)
2399 : {
2400 6953999 : pari_sp av = avma;
2401 6953999 : const double ETA = 0.51;
2402 6953999 : const long keepfirst = flag & LLL_KEEP_FIRST;
2403 6953999 : long p, zeros = -1, n = lg(x)-1, is_upper, is_lower, useflatter = 0, rank;
2404 6953999 : GEN G, B, U, L = NULL;
2405 : pari_timer T;
2406 6953999 : long thre[]={31783,34393,20894,22525,13533,1928,672,671,
2407 : 422,506,315,313,222,205,167,154,139,138,
2408 : 110,120,98,94,81,75,74,64,74,74,
2409 : 79,96,112,111,105,104,96,86,84,78,75,70,66,62,62,57,56,47,45,52,50,44,48,42,36,35,35,34,40,33,34,32,36,31,
2410 : 38,38,40,38,38,37,35,31,34,36,34,32,34,32,28,27,25,31,25,27,28,26,25,21,21,25,25,22,21,24,24,22,21,23,22,22,22,22,21,24,21,22,19,20,19,20,19,19,19,18,19,18,18,20,19,20,18,19,18,21,18,20,18,18};
2411 6953999 : long thsn[]={23280,30486,50077,44136,78724,15690,1801,1611,
2412 : 981,1359,978,1042,815,866,788,775,726,712,
2413 : 626,613,548,564,474,481,504,447,453,508,
2414 : 705,794,1008,946,767,898,886,763,842,757,
2415 : 725,774,639,655,705,627,635,704,511,613,
2416 : 583,595,568,640,541,640,567,540,577,584,
2417 : 546,509,526,572,637,746,772,743,743,742,800,708,832,768,707,692,692,768,696,635,709,694,768,719,655,569,590,644,685,623,627,720,633,636,602,635,575,631,642,647,632,656,573,511,688,640,528,616,511,559,601,620,635,688,608,768,658,582,644,704,555,673,600,601,641,661,601,670};
2418 6953999 : if (n <= 1) return lll_trivial(x, flag);
2419 6844555 : if (nbrows(x)==0)
2420 : {
2421 15171 : if (flag & LLL_KER) return matid(n);
2422 15171 : if (flag & (LLL_INPLACE|LLL_IM)) return cgetg(1,t_MAT);
2423 0 : retmkvec2(matid(n), cgetg(1,t_MAT));
2424 : }
2425 6829557 : if (n==2 && nbrows(x)==2 && (flag&LLL_IM) && !keepfirst)
2426 : {
2427 4701009 : U = ZM2_lll_norms(x, flag, pN);
2428 4701008 : if (U) return U;
2429 : }
2430 2128929 : if (flag & LLL_GRAM)
2431 60276 : { G = x; B = NULL; U = matid(n); is_upper = 0; is_lower = 0; }
2432 : else
2433 : {
2434 2068653 : G = NULL; B = x; U = (flag & LLL_INPLACE)? NULL: matid(n);
2435 2068652 : is_upper = (flag & LLL_UPPER) || ZM_is_upper(B);
2436 2068652 : is_lower = !B || is_upper || keepfirst ? 0: ZM_is_lower(B);
2437 2068650 : if (is_lower) L = RgM_flip(B);
2438 : }
2439 2128926 : rank = (flag&LLL_NOFLATTER) ? 0: ZM_rank(x);
2440 2128921 : if (n > 2 && !(flag&LLL_NOFLATTER))
2441 1733003 : {
2442 1688745 : GEN R = B ? (is_upper ? B : (is_lower ? L : get_gramschmidt(B, rank)))
2443 3421740 : : get_gaussred(G, rank);
2444 1732997 : if (R)
2445 : {
2446 1116744 : long spr = spread(R), sz = mpexpo(gsupnorm(R, DEFAULTPREC)), thr;
2447 1116749 : if (DEBUGLEVEL>=5) err_printf("LLL: dim %ld, size %ld, spread %ld\n",n, sz, spr);
2448 1116749 : if ((is_upper && ZM_is_knapsack(B)) || (is_lower && ZM_is_knapsack(L)))
2449 92364 : thr = thsn[minss(n-3,numberof(thsn)-1)];
2450 : else
2451 : {
2452 1024385 : thr = thre[minss(n-3,numberof(thre)-1)];
2453 1024385 : if (n>=10) sz = spr;
2454 : }
2455 1116750 : useflatter = sz >= thr;
2456 : } else
2457 616253 : useflatter = 1;
2458 395916 : } else useflatter = 0;
2459 2128919 : if(DEBUGLEVEL>=4) timer_start(&T);
2460 2128919 : if (useflatter)
2461 : {
2462 617331 : if (is_lower)
2463 : {
2464 0 : fplll_flatter(&G, &L, &U, rank, flag | LLL_UPPER);
2465 0 : B = RgM_flop(L);
2466 0 : if (U) U = RgM_flop(U);
2467 : }
2468 : else
2469 617331 : fplll_flatter(&G, &B, &U, rank, flag | (is_upper? LLL_UPPER:0));
2470 617332 : if (DEBUGLEVEL>=4 && !(flag & LLL_NOCERTIFY))
2471 0 : timer_printf(&T, "FLATTER");
2472 : }
2473 2128918 : if (!(flag & LLL_GRAM))
2474 : {
2475 : long t;
2476 2068642 : B = gcopy(B);
2477 2068645 : if(DEBUGLEVEL>=4)
2478 0 : err_printf("Entering L^2 (double): dim %ld, LLL-parameters (%.3f,%.3f)\n",
2479 : n, DELTA,ETA);
2480 2068645 : zeros = fplll_fast(&B, &U, DELTA, ETA, keepfirst);
2481 2068655 : if (DEBUGLEVEL>=4) timer_printf(&T, zeros < 0? "LLL (failed)": "LLL");
2482 2072794 : for (p = DEFAULTPREC, t = 0; zeros < 0 && t < 1 ; p += EXTRAPREC64, t++)
2483 : {
2484 4139 : if (DEBUGLEVEL>=4)
2485 0 : err_printf("Entering L^2 (heuristic): LLL-parameters (%.3f,%.3f), prec = %d/%d\n", DELTA, ETA, p, p);
2486 4139 : zeros = fplll_heuristic(&B, &U, DELTA, ETA, keepfirst, p, p);
2487 4139 : gc_lll(av, 2, &B, &U);
2488 4139 : if (DEBUGLEVEL>=4) timer_printf(&T, zeros < 0? "LLL (failed)": "LLL");
2489 : }
2490 : } else
2491 60276 : G = gcopy(G);
2492 2128931 : if (zeros < 0 || !(flag & LLL_NOCERTIFY))
2493 : {
2494 1567955 : if(DEBUGLEVEL>=4)
2495 0 : err_printf("Entering L^2 (dpe): LLL-parameters (%.3f,%.3f)\n", DELTA,ETA);
2496 1567955 : zeros = fplll_dpe(&G, &B, &U, pN, DELTA, ETA, keepfirst);
2497 1567935 : if (DEBUGLEVEL>=4) timer_printf(&T, zeros < 0? "LLL (failed)": "LLL");
2498 1567939 : if (zeros < 0)
2499 0 : for (p = DEFAULTPREC;; p += EXTRAPREC64)
2500 : {
2501 0 : if (DEBUGLEVEL>=4)
2502 0 : err_printf("Entering L^2: LLL-parameters (%.3f,%.3f), prec = %d\n",
2503 : DELTA,ETA, p);
2504 0 : zeros = fplll(&G, &B, &U, pN, DELTA, ETA, keepfirst, p);
2505 0 : if (DEBUGLEVEL>=4) timer_printf(&T, zeros < 0? "LLL (failed)": "LLL");
2506 0 : if (zeros >= 0) break;
2507 0 : gc_lll(av, 3, &G, &B, &U);
2508 : }
2509 : }
2510 2128915 : return lll_finish(U? U: B, zeros, flag);
2511 : }
2512 :
2513 : /********************************************************************/
2514 : /** **/
2515 : /** LLL OVER K[X] **/
2516 : /** **/
2517 : /********************************************************************/
2518 : static int
2519 504 : pslg(GEN x)
2520 : {
2521 : long tx;
2522 504 : if (gequal0(x)) return 2;
2523 448 : tx = typ(x); return is_scalar_t(tx)? 3: lg(x);
2524 : }
2525 :
2526 : static int
2527 196 : REDgen(long k, long l, GEN h, GEN L, GEN B)
2528 : {
2529 196 : GEN q, u = gcoeff(L,k,l);
2530 : long i;
2531 :
2532 196 : if (pslg(u) < pslg(B)) return 0;
2533 :
2534 140 : q = gneg(gdeuc(u,B));
2535 140 : gel(h,k) = gadd(gel(h,k), gmul(q,gel(h,l)));
2536 140 : for (i=1; i<l; i++) gcoeff(L,k,i) = gadd(gcoeff(L,k,i), gmul(q,gcoeff(L,l,i)));
2537 140 : gcoeff(L,k,l) = gadd(gcoeff(L,k,l), gmul(q,B)); return 1;
2538 : }
2539 :
2540 : static int
2541 196 : do_SWAPgen(GEN h, GEN L, GEN B, long k, GEN fl, int *flc)
2542 : {
2543 : GEN p1, la, la2, Bk;
2544 : long ps1, ps2, i, j, lx;
2545 :
2546 196 : if (!fl[k-1]) return 0;
2547 :
2548 140 : la = gcoeff(L,k,k-1); la2 = gsqr(la);
2549 140 : Bk = gel(B,k);
2550 140 : if (fl[k])
2551 : {
2552 56 : GEN q = gadd(la2, gmul(gel(B,k-1),gel(B,k+1)));
2553 56 : ps1 = pslg(gsqr(Bk));
2554 56 : ps2 = pslg(q);
2555 56 : if (ps1 <= ps2 && (ps1 < ps2 || !*flc)) return 0;
2556 28 : *flc = (ps1 != ps2);
2557 28 : gel(B,k) = gdiv(q, Bk);
2558 : }
2559 :
2560 112 : swap(gel(h,k-1), gel(h,k)); lx = lg(L);
2561 112 : for (j=1; j<k-1; j++) swap(gcoeff(L,k-1,j), gcoeff(L,k,j));
2562 112 : if (fl[k])
2563 : {
2564 28 : for (i=k+1; i<lx; i++)
2565 : {
2566 0 : GEN t = gcoeff(L,i,k);
2567 0 : p1 = gsub(gmul(gel(B,k+1),gcoeff(L,i,k-1)), gmul(la,t));
2568 0 : gcoeff(L,i,k) = gdiv(p1, Bk);
2569 0 : p1 = gadd(gmul(la,gcoeff(L,i,k-1)), gmul(gel(B,k-1),t));
2570 0 : gcoeff(L,i,k-1) = gdiv(p1, Bk);
2571 : }
2572 : }
2573 84 : else if (!gequal0(la))
2574 : {
2575 28 : p1 = gdiv(la2, Bk);
2576 28 : gel(B,k+1) = gel(B,k) = p1;
2577 28 : for (i=k+2; i<=lx; i++) gel(B,i) = gdiv(gmul(p1,gel(B,i)),Bk);
2578 28 : for (i=k+1; i<lx; i++)
2579 0 : gcoeff(L,i,k-1) = gdiv(gmul(la,gcoeff(L,i,k-1)), Bk);
2580 28 : for (j=k+1; j<lx-1; j++)
2581 0 : for (i=j+1; i<lx; i++)
2582 0 : gcoeff(L,i,j) = gdiv(gmul(p1,gcoeff(L,i,j)), Bk);
2583 : }
2584 : else
2585 : {
2586 56 : gcoeff(L,k,k-1) = gen_0;
2587 56 : for (i=k+1; i<lx; i++)
2588 : {
2589 0 : gcoeff(L,i,k) = gcoeff(L,i,k-1);
2590 0 : gcoeff(L,i,k-1) = gen_0;
2591 : }
2592 56 : gel(B,k) = gel(B,k-1); fl[k] = 1; fl[k-1] = 0;
2593 : }
2594 112 : return 1;
2595 : }
2596 :
2597 : static void
2598 168 : incrementalGSgen(GEN x, GEN L, GEN B, long k, GEN fl)
2599 : {
2600 168 : GEN u = NULL; /* gcc -Wall */
2601 : long i, j;
2602 420 : for (j = 1; j <= k; j++)
2603 252 : if (j==k || fl[j])
2604 : {
2605 252 : u = gcoeff(x,k,j);
2606 252 : if (!is_extscalar_t(typ(u))) pari_err_TYPE("incrementalGSgen",u);
2607 336 : for (i=1; i<j; i++)
2608 84 : if (fl[i])
2609 : {
2610 84 : u = gsub(gmul(gel(B,i+1),u), gmul(gcoeff(L,k,i),gcoeff(L,j,i)));
2611 84 : u = gdiv(u, gel(B,i));
2612 : }
2613 252 : gcoeff(L,k,j) = u;
2614 : }
2615 168 : if (gequal0(u)) gel(B,k+1) = gel(B,k);
2616 : else
2617 : {
2618 112 : gel(B,k+1) = gcoeff(L,k,k); gcoeff(L,k,k) = gen_1; fl[k] = 1;
2619 : }
2620 168 : }
2621 :
2622 : static GEN
2623 168 : lllgramallgen(GEN x, long flag)
2624 : {
2625 168 : long lx = lg(x), i, j, k, l, n;
2626 : pari_sp av;
2627 : GEN B, L, h, fl;
2628 : int flc;
2629 :
2630 168 : n = lx-1; if (n<=1) return lll_trivial(x,flag);
2631 84 : if (lgcols(x) != lx) pari_err_DIM("lllgramallgen");
2632 :
2633 84 : fl = cgetg(lx, t_VECSMALL);
2634 :
2635 84 : av = avma;
2636 84 : B = scalarcol_shallow(gen_1, lx);
2637 84 : L = cgetg(lx,t_MAT);
2638 252 : for (j=1; j<lx; j++) { gel(L,j) = zerocol(n); fl[j] = 0; }
2639 :
2640 84 : h = matid(n);
2641 252 : for (i=1; i<lx; i++)
2642 168 : incrementalGSgen(x, L, B, i, fl);
2643 84 : flc = 0;
2644 84 : for(k=2;;)
2645 : {
2646 196 : if (REDgen(k, k-1, h, L, gel(B,k))) flc = 1;
2647 196 : if (do_SWAPgen(h, L, B, k, fl, &flc)) { if (k > 2) k--; }
2648 : else
2649 : {
2650 84 : for (l=k-2; l>=1; l--)
2651 0 : if (REDgen(k, l, h, L, gel(B,l+1))) flc = 1;
2652 84 : if (++k > n) break;
2653 : }
2654 112 : if (gc_needed(av,1))
2655 : {
2656 0 : if(DEBUGMEM>1) pari_warn(warnmem,"lllgramallgen");
2657 0 : gerepileall(av,3,&B,&L,&h);
2658 : }
2659 : }
2660 140 : k=1; while (k<lx && !fl[k]) k++;
2661 84 : return lll_finish(h,k-1,flag);
2662 : }
2663 :
2664 : static GEN
2665 168 : lllallgen(GEN x, long flag)
2666 : {
2667 168 : pari_sp av = avma;
2668 168 : if (!(flag & LLL_GRAM)) x = gram_matrix(x);
2669 84 : else if (!RgM_is_square_mat(x)) pari_err_DIM("qflllgram");
2670 168 : return gerepilecopy(av, lllgramallgen(x, flag));
2671 : }
2672 : GEN
2673 42 : lllgen(GEN x) { return lllallgen(x, LLL_IM); }
2674 : GEN
2675 42 : lllkerimgen(GEN x) { return lllallgen(x, LLL_ALL); }
2676 : GEN
2677 42 : lllgramgen(GEN x) { return lllallgen(x, LLL_IM|LLL_GRAM); }
2678 : GEN
2679 42 : lllgramkerimgen(GEN x) { return lllallgen(x, LLL_ALL|LLL_GRAM); }
2680 :
2681 : static GEN
2682 36699 : lllall(GEN x, long flag)
2683 36699 : { pari_sp av = avma; return gerepilecopy(av, ZM_lll(x, LLLDFT, flag)); }
2684 : GEN
2685 183 : lllint(GEN x) { return lllall(x, LLL_IM); }
2686 : GEN
2687 35 : lllkerim(GEN x) { return lllall(x, LLL_ALL); }
2688 : GEN
2689 36439 : lllgramint(GEN x)
2690 36439 : { if (!RgM_is_square_mat(x)) pari_err_DIM("qflllgram");
2691 36439 : return lllall(x, LLL_IM | LLL_GRAM); }
2692 : GEN
2693 35 : lllgramkerim(GEN x)
2694 35 : { if (!RgM_is_square_mat(x)) pari_err_DIM("qflllgram");
2695 35 : return lllall(x, LLL_ALL | LLL_GRAM); }
2696 :
2697 : GEN
2698 5268849 : lllfp(GEN x, double D, long flag)
2699 : {
2700 5268849 : long n = lg(x)-1;
2701 5268849 : pari_sp av = avma;
2702 : GEN h;
2703 5268849 : if (n <= 1) return lll_trivial(x,flag);
2704 4608452 : if (flag & LLL_GRAM)
2705 : {
2706 9270 : if (!RgM_is_square_mat(x)) pari_err_DIM("qflllgram");
2707 9256 : if (isinexact(x))
2708 : {
2709 9165 : x = RgM_Cholesky(x, gprecision(x));
2710 9165 : if (!x) return NULL;
2711 9165 : flag &= ~LLL_GRAM;
2712 : }
2713 : }
2714 4608438 : h = ZM_lll(RgM_rescale_to_int(x), D, flag);
2715 4608390 : return gerepilecopy(av, h);
2716 : }
2717 :
2718 : GEN
2719 9089 : lllgram(GEN x) { return lllfp(x,LLLDFT,LLL_GRAM|LLL_IM); }
2720 : GEN
2721 1175019 : lll(GEN x) { return lllfp(x,LLLDFT,LLL_IM); }
2722 :
2723 : static GEN
2724 63 : qflllgram(GEN x)
2725 : {
2726 63 : GEN T = lllgram(x);
2727 42 : if (!T) pari_err_PREC("qflllgram");
2728 42 : return T;
2729 : }
2730 :
2731 : GEN
2732 301 : qflll0(GEN x, long flag)
2733 : {
2734 301 : if (typ(x) != t_MAT) pari_err_TYPE("qflll",x);
2735 301 : switch(flag)
2736 : {
2737 49 : case 0: return lll(x);
2738 63 : case 1: return lllfp(x, LLLDFT, LLL_IM | LLL_NOFLATTER);
2739 49 : case 2: RgM_check_ZM(x,"qflll"); return lllintpartial(x);
2740 7 : case 3: RgM_check_ZM(x,"qflll"); return lllall(x, LLL_INPLACE);
2741 49 : case 4: RgM_check_ZM(x,"qflll"); return lllkerim(x);
2742 42 : case 5: return lllkerimgen(x);
2743 42 : case 8: return lllgen(x);
2744 0 : default: pari_err_FLAG("qflll");
2745 : }
2746 : return NULL; /* LCOV_EXCL_LINE */
2747 : }
2748 :
2749 : GEN
2750 245 : qflllgram0(GEN x, long flag)
2751 : {
2752 245 : if (typ(x) != t_MAT) pari_err_TYPE("qflllgram",x);
2753 245 : switch(flag)
2754 : {
2755 63 : case 0: return qflllgram(x);
2756 49 : case 1: return lllfp(x, LLLDFT, LLL_IM | LLL_GRAM | LLL_NOFLATTER);
2757 49 : case 4: RgM_check_ZM(x,"qflllgram"); return lllgramkerim(x);
2758 42 : case 5: return lllgramkerimgen(x);
2759 42 : case 8: return lllgramgen(x);
2760 0 : default: pari_err_FLAG("qflllgram");
2761 : }
2762 : return NULL; /* LCOV_EXCL_LINE */
2763 : }
2764 :
2765 : /********************************************************************/
2766 : /** **/
2767 : /** INTEGRAL KERNEL (LLL REDUCED) **/
2768 : /** **/
2769 : /********************************************************************/
2770 : static GEN
2771 70 : kerint0(GEN M)
2772 : {
2773 : /* return ZM_lll(M, LLLDFT, LLL_KER); */
2774 70 : GEN U, H = ZM_hnflll(M,&U,1);
2775 70 : long d = lg(M)-lg(H);
2776 70 : if (!d) return cgetg(1, t_MAT);
2777 70 : return ZM_lll(vecslice(U,1,d), LLLDFT, LLL_INPLACE);
2778 : }
2779 : GEN
2780 42 : kerint(GEN M)
2781 : {
2782 42 : pari_sp av = avma;
2783 42 : return gerepilecopy(av, kerint0(M));
2784 : }
2785 : /* OBSOLETE: use kerint */
2786 : GEN
2787 28 : matkerint0(GEN M, long flag)
2788 : {
2789 28 : pari_sp av = avma;
2790 28 : if (typ(M) != t_MAT) pari_err_TYPE("matkerint",M);
2791 28 : M = Q_primpart(M);
2792 28 : RgM_check_ZM(M, "kerint");
2793 28 : switch(flag)
2794 : {
2795 28 : case 0:
2796 28 : case 1: return gerepilecopy(av, kerint0(M));
2797 0 : default: pari_err_FLAG("matkerint");
2798 : }
2799 : return NULL; /* LCOV_EXCL_LINE */
2800 : }
|
