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 50516 : RgM_is_square_mat(GEN x) { long l = lg(x); return l == 1 || l == lgcols(x); }
22 :
23 : static long
24 4221788 : ZM_is_upper(GEN R)
25 : {
26 4221788 : long i,j, l = lg(R);
27 4221788 : if (l != lgcols(R)) return 0;
28 8160514 : for(i = 1; i < l; i++)
29 8851203 : for(j = 1; j < i; j++)
30 4572732 : if (signe(gcoeff(R,i,j))) return 0;
31 251338 : return 1;
32 : }
33 :
34 : static long
35 609771 : ZM_is_knapsack(GEN R)
36 : {
37 609771 : long i,j, l = lg(R);
38 609771 : if (l != lgcols(R)) return 0;
39 856805 : for(i = 2; i < l; i++)
40 3042117 : for(j = 1; j < l; j++)
41 2795083 : if ( i!=j && signe(gcoeff(R,i,j))) return 0;
42 94824 : return 1;
43 : }
44 :
45 : static long
46 1190687 : ZM_is_lower(GEN R)
47 : {
48 1190687 : long i,j, l = lg(R);
49 1190687 : if (l != lgcols(R)) return 0;
50 2077544 : for(i = 1; i < l; i++)
51 2394753 : for(j = 1; j < i; j++)
52 1294156 : if (signe(gcoeff(R,j,i))) return 0;
53 34653 : return 1;
54 : }
55 :
56 : static GEN
57 34654 : RgM_flip(GEN R)
58 : {
59 : GEN M;
60 : long i,j,l;
61 34654 : M = cgetg_copy(R, &l);
62 179309 : for(i = 1; i < l; i++)
63 : {
64 144655 : gel(M,i) = cgetg(l, t_COL);
65 890846 : for(j = 1; j < l; j++)
66 746191 : gmael(M,i,j) = gmael(R,l-i, l-j);
67 : }
68 34654 : 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 23076176 : mpabscmp(GEN x, GEN y)
89 : {
90 23076176 : return (typ(x)==t_INT) ? abscmpii(x,y) : abscmprr(x,y);
91 : }
92 :
93 : /****************************************************************************/
94 : /*** FLATTER ***/
95 : /****************************************************************************/
96 :
97 : /* Implementation of "FLATTER" algorithm based on
98 : <https://eprint.iacr.org/2023/237>
99 : Fast Practical Lattice Reduction through Iterated Compression
100 :
101 : Keegan Ryan, University of California, San Diego
102 : Nadia Heninger, University of California, San Diego
103 : BA20230925 */
104 :
105 : static long
106 5077836 : drop(GEN R)
107 : {
108 5077836 : long i, n = lg(R)-1;
109 5077836 : long s = 0, m = mpexpo(gcoeff(R, 1, 1));
110 20783955 : for (i = 2; i <= n; ++i)
111 : {
112 15706118 : if (mpabscmp(gcoeff(R, i, i), gcoeff(R, i - 1, i - 1)) >= 0)
113 : {
114 9591453 : s += m - mpexpo(gcoeff(R, i - 1, i - 1));
115 9591453 : m = mpexpo(gcoeff(R, i, i));
116 : }
117 : }
118 5077837 : s += m - mpexpo(gcoeff(R, n, n));
119 5077836 : return s;
120 : }
121 :
122 : static long
123 4181086 : gsisinv(GEN M)
124 : {
125 4181086 : long i, l = lg(M);
126 21590311 : for (i = 1; i < l; ++i)
127 17409629 : if (signe(gmael(M, i, i))==0)
128 404 : return 0;
129 4180682 : return 1;
130 : }
131 :
132 : INLINE long
133 6484430 : nbits2prec64(long n)
134 : {
135 6484430 : return nbits2prec(((n+63)>>6)<<6);
136 : }
137 :
138 : static GEN
139 4039 : RgM_Cholesky_dynprec(GEN M)
140 : {
141 4039 : pari_sp ltop = avma;
142 : GEN L;
143 4039 : long minprec = 3*(lg(M)-1) + 30, bitprec = minprec, prec;
144 : while (1)
145 4295 : {
146 : long mbitprec;
147 8334 : prec = nbits2prec64(bitprec);
148 8334 : L = RgM_Cholesky(RgM_gtofp(M, prec), prec);
149 8334 : if (!L)
150 : {
151 1801 : bitprec *= 2;
152 1801 : set_avma(ltop);
153 1801 : continue;
154 : }
155 6533 : mbitprec = minprec + 2*drop(L);
156 6533 : if (bitprec >= mbitprec)
157 4039 : break;
158 2494 : bitprec = mbitprec;
159 2494 : set_avma(ltop);
160 : }
161 4039 : return gerepilecopy(ltop, L);
162 : }
163 :
164 : static GEN
165 1455 : gramschmidt_upper(GEN M)
166 : {
167 1455 : long bitprec = 3*(lg(M)-1) + 30 + 2*drop(M);
168 1455 : return RgM_gtofp(M, nbits2prec64(bitprec));
169 : }
170 :
171 : static GEN
172 2708378 : gramschmidt_dynprec(GEN M)
173 : {
174 2708378 : pari_sp ltop = avma;
175 2708378 : long minprec = 3*(lg(M)-1) + 30, bitprec = minprec;
176 2708378 : if (ZM_is_upper(M)) return gramschmidt_upper(M);
177 : while (1)
178 2624721 : {
179 : GEN B, Q, L;
180 5331644 : long prec = nbits2prec64(bitprec), mbitprec;
181 5331634 : if (!QR_init(RgM_gtofp(M, prec), &B, &Q, &L, prec) || !gsisinv(L))
182 : {
183 1620013 : bitprec *= 2;
184 1620013 : set_avma(ltop);
185 1620018 : continue;
186 : }
187 3711622 : mbitprec = minprec + 2*drop(L);
188 3711622 : if (bitprec >= mbitprec)
189 2706919 : return gerepilecopy(ltop, shallowtrans(L));
190 1004703 : bitprec = mbitprec;
191 1004703 : set_avma(ltop);
192 : }
193 : }
194 : /* return -T1 * round(T1^-1*(R1^-1*R2)*T3) */
195 : static GEN
196 1354190 : sizered(GEN T1, GEN T3, GEN R1, GEN R2)
197 : {
198 1354190 : pari_sp ltop = avma;
199 : long e;
200 1354190 : return gerepileupto(ltop, gmul(ZM_neg(T1), grndtoi(gmul(ZM_inv(T1,NULL),
201 : gmul(RgM_mul(RgM_inv_upper(R1), R2), T3)), &e)));
202 : }
203 :
204 : static GEN
205 1354190 : flat(GEN M, long flag, GEN *pt_T, long *pt_s)
206 : {
207 1354190 : pari_sp ltop = avma;
208 : GEN R, R1, R2, R3, T1, T2, T3, T, S;
209 1354190 : long k = lg(M)-1, n = k>>1, n2 = k - n, m = n>>1;
210 1354190 : long keepfirst = flag & LLL_KEEP_FIRST, inplace = flag & LLL_INPLACE;
211 : /* for k = 3, we want n = 1; n2 = 2; m = 0 */
212 : /* for k = 5, n = 2; n2 = 3; m = 1 */
213 1354190 : R = gramschmidt_dynprec(M);
214 1354189 : R1 = matslice(R, 1, n, 1, n);
215 1354188 : R2 = matslice(R, 1, n, n + 1, k);
216 1354188 : R3 = matslice(R, n + 1, k, n + 1, k);
217 1354190 : T1 = lllfp(R1, 0.99, LLL_IM| LLL_UPPER| LLL_NOCERTIFY| (keepfirst ? LLL_KEEP_FIRST: 0));
218 1354189 : T3 = lllfp(R3, 0.99, LLL_IM| LLL_UPPER| LLL_NOCERTIFY);
219 1354190 : T2 = sizered(T1, T3, R1, R2);
220 1354190 : T = shallowmatconcat(mkmat22(T1,T2,gen_0,T3));
221 1354190 : M = ZM_mul(M, T);
222 1354189 : R = gramschmidt_dynprec(M);
223 1354190 : R3 = matslice(R, m + 1, m + n2, m + 1, m + n2);
224 1354190 : T3 = lllfp(R3, 0.99, LLL_IM| LLL_UPPER| LLL_NOCERTIFY);
225 2708380 : S = shallowmatconcat(diagonal(
226 575789 : m == 0 ? mkvec2(T3, matid(k - m - n2))
227 0 : : m+n2 == k ? mkvec2(matid(m), T3)
228 778401 : : mkvec3(matid(m), T3, matid(k - m - n2))));
229 1354190 : M = ZM_mul(M, S);
230 1354187 : if (!inplace) *pt_T = ZM_mul(T, S);
231 1354189 : *pt_s = drop(R);
232 1354189 : return gc_all(ltop, inplace ? 1: 2, &M, pt_T);
233 : }
234 :
235 : static GEN
236 627542 : ZM_flatter(GEN M, long flag)
237 : {
238 627542 : pari_sp ltop = avma, btop;
239 627542 : long i, n = lg(M)-1;
240 627542 : long s = -1;
241 : GEN T;
242 : pari_timer ti;
243 627542 : long lti = 1;
244 627542 : long inplace = flag & LLL_INPLACE, cert = !(flag & LLL_NOCERTIFY);
245 627542 : T = matid(n);
246 627542 : if (DEBUGLEVEL>=3)
247 : {
248 0 : timer_start(&ti);
249 0 : if (cert)
250 0 : err_printf("flatter dim = %ld size = %ld\n", n, ZM_max_expi(M));
251 : }
252 627542 : btop = avma;
253 627542 : for (i = 1;;i++)
254 726648 : {
255 : long t;
256 1354190 : GEN U, M2 = flat(M, flag, &U, &t);
257 1354189 : if (t==0) { s = t; break; }
258 770634 : if (s >= 0)
259 : {
260 440398 : if (s==t)
261 43986 : break;
262 396412 : if (s<t && i > 20)
263 : {
264 0 : if (DEBUGLEVEL >= 3) err_printf("BACK:%ld:%ld:%g\n", n, i, s);
265 0 : break;
266 : }
267 : }
268 726648 : if (DEBUGLEVEL>=3 && (cert || timer_get(&ti) > 1000))
269 : {
270 0 : if (i==lti)
271 0 : timer_printf(&ti, "FLATTER, dim %ld, steps %ld: \t slope=%0.10g", n, i, ((double)t)/n);
272 : else
273 0 : timer_printf(&ti, "FLATTER, dim %ld, steps %ld-%ld: \t slope=%0.10g", n, lti,i,((double)t)/n);
274 0 : lti = i+1;
275 : }
276 726648 : s = t;
277 726648 : M = M2;
278 726648 : if (!inplace) T = ZM_mul(T, U);
279 726648 : if (gc_needed(btop, 1))
280 0 : gerepileall(btop, 2, &M, &T);
281 : }
282 627541 : if (DEBUGLEVEL>=3 && (cert || timer_get(&ti) > 1000))
283 : {
284 0 : if (i==lti)
285 0 : timer_printf(&ti, "FLATTER, dim %ld, final \t slope=%0.10g", n, ((double)s)/n);
286 : else
287 0 : timer_printf(&ti, "FLATTER, dim %ld, steps %ld-final:\t slope=%0.10g", n, lti, ((double)s)/n);
288 : }
289 627541 : return gerepilecopy(ltop, inplace ? M: T);
290 : }
291 :
292 : static GEN
293 625559 : ZM_flatter_rank(GEN M, long rank, long flag)
294 : {
295 : pari_timer ti;
296 625559 : pari_sp ltop = avma;
297 : GEN T;
298 625559 : long i, n = lg(M)-1;
299 625559 : if (rank == n) return ZM_flatter(M, flag);
300 3764 : T = matid(n);
301 3764 : if (DEBUGLEVEL>=3) timer_start(&ti);
302 3764 : for (i = 1;; i++)
303 1981 : {
304 5745 : GEN S = ZM_flatter(vconcat(gshift(M,i),matid(n)), flag);
305 5745 : if (DEBUGLEVEL>=3)
306 0 : timer_printf(&ti,"FLATTERRANK step %ld: %ld",i,expi(gnorml2(S)));
307 5745 : if (ZM_isidentity(S)) break;
308 1981 : T = ZM_mul(T, S);
309 1981 : M = ZM_mul(M, S);
310 1981 : if (gc_needed(ltop, 1))
311 0 : gerepileall(ltop, 2, &M, &T);
312 : }
313 3764 : return gerepilecopy(ltop, T);
314 : }
315 :
316 : static GEN
317 4039 : flattergram_i(GEN M, long flag, long *pt_s)
318 : {
319 4039 : pari_sp ltop = avma;
320 : GEN R, T;
321 4039 : R = RgM_Cholesky_dynprec(M);
322 4039 : *pt_s = drop(R);
323 4039 : T = lllfp(R, 0.99, LLL_IM| LLL_UPPER| LLL_NOCERTIFY | (flag&LLL_KEEP_FIRST));
324 4039 : return gerepilecopy(ltop, T);
325 : }
326 :
327 : static GEN
328 1316 : ZM_flattergram(GEN M, long flag)
329 : {
330 1316 : pari_sp ltop = avma, btop;
331 : GEN T;
332 1316 : long i, n = lg(M)-1;
333 : pari_timer ti;
334 1316 : long s = -1;
335 1316 : if (DEBUGLEVEL>=3)
336 : {
337 0 : timer_start(&ti);
338 0 : err_printf("FLATTERGRAM dim = %ld size = %ld\n", n, expi(gnorml2(M)));
339 : }
340 1316 : btop = avma;
341 1316 : T = matid(n);
342 1316 : for (i = 1;; i++)
343 2723 : {
344 : long t;
345 4039 : GEN S = flattergram_i(M, flag, &t);
346 4039 : t = expi(gnorml2(S));
347 4039 : if (t==0) { s = t; break; }
348 4039 : if (s)
349 : {
350 4039 : double st = s-t;
351 4039 : if (st == 0)
352 1316 : break;
353 2723 : if (st < 0 && i > 20)
354 : {
355 0 : if (DEBUGLEVEL >= 3) err_printf("BACK:%ld:%ld:%0.10g\n", n, i, ((double)s)/n);
356 0 : break;
357 : }
358 : }
359 2723 : T = ZM_mul(T, S);
360 2723 : M = ZM_mul(shallowtrans(S), ZM_mul(M, S));
361 2723 : s = t;
362 2723 : if (DEBUGLEVEL >= 3)
363 0 : timer_printf(&ti, "FLATTERGRAM, dim %ld step %ld, slope=%0.10g", n, i, ((double)s)/n);
364 2723 : if (gc_needed(btop, 1))
365 0 : gerepileall(btop, 2, &M, &T);
366 : }
367 1316 : if (DEBUGLEVEL >= 3)
368 0 : timer_printf(&ti, "FLATTERGRAM, dim %ld, slope=%0.10g", n, ((double)s)/n);
369 1316 : return gerepilecopy(ltop, T);
370 : }
371 :
372 : static GEN
373 1316 : ZM_flattergram_rank(GEN M, long rank, long flag)
374 : {
375 : pari_timer ti;
376 1316 : pari_sp ltop = avma;
377 : GEN T;
378 1316 : long i, n = lg(M)-1;
379 1316 : if (rank == n) return ZM_flattergram(M, flag);
380 0 : T = matid(n);
381 0 : if (DEBUGLEVEL>=3) timer_start(&ti);
382 0 : for (i = 1;; i++)
383 0 : {
384 0 : GEN S = ZM_flattergram(RgM_Rg_add(gshift(M, i), gen_1), flag);
385 0 : if (DEBUGLEVEL>=3)
386 0 : timer_printf(&ti,"FLATTERGRAMRANK step %ld: %ld",i,expi(gnorml2(S)));
387 0 : if (ZM_isidentity(S)) break;
388 0 : T = ZM_mul(T, S);
389 0 : M = ZM_mul(shallowtrans(S), ZM_mul(M, S));
390 0 : if (gc_needed(ltop, 1))
391 0 : gerepileall(ltop, 2, &M, &T);
392 : }
393 0 : return gerepilecopy(ltop, T);
394 : }
395 :
396 : static double
397 10931485 : pari_rint(double a)
398 : {
399 : #ifdef HAS_RINT
400 10931485 : return rint(a);
401 : #else
402 : const double two_to_52 = 4.5035996273704960e+15;
403 : double fa = fabs(a);
404 : double r = two_to_52 + fa;
405 : if (fa >= two_to_52) {
406 : r = a;
407 : } else {
408 : r = r - two_to_52;
409 : if (a < 0) r = -r;
410 : }
411 : return r;
412 : #endif
413 : }
414 :
415 : /* default quality ratio for LLL */
416 : static const double LLLDFT = 0.99;
417 :
418 : /* assume flag & (LLL_KER|LLL_IM|LLL_ALL). LLL_INPLACE implies LLL_IM */
419 : static GEN
420 768600 : lll_trivial(GEN x, long flag)
421 : {
422 768600 : if (lg(x) == 1)
423 : { /* dim x = 0 */
424 15448 : if (! (flag & LLL_ALL)) return cgetg(1,t_MAT);
425 28 : retmkvec2(cgetg(1,t_MAT), cgetg(1,t_MAT));
426 : }
427 : /* dim x = 1 */
428 753152 : if (gequal0(gel(x,1)))
429 : {
430 91 : if (flag & LLL_KER) return matid(1);
431 91 : if (flag & (LLL_IM|LLL_INPLACE)) return cgetg(1,t_MAT);
432 28 : retmkvec2(matid(1), cgetg(1,t_MAT));
433 : }
434 753054 : if (flag & LLL_INPLACE) return gcopy(x);
435 649834 : if (flag & LLL_KER) return cgetg(1,t_MAT);
436 649834 : if (flag & LLL_IM) return matid(1);
437 27 : retmkvec2(cgetg(1,t_MAT), (flag & LLL_GRAM)? gcopy(x): matid(1));
438 : }
439 :
440 : /* vecslice(x,#x-k,#x) in place. Works for t_MAT, t_VEC/t_COL */
441 : static GEN
442 2066712 : vectail_inplace(GEN x, long k)
443 : {
444 2066712 : if (!k) return x;
445 57605 : x[k] = ((ulong)x[0] & ~LGBITS) | evallg(lg(x) - k);
446 57608 : return x + k;
447 : }
448 :
449 : /* k = dim Kernel */
450 : static GEN
451 2139790 : lll_finish(GEN h, long k, long flag)
452 : {
453 : GEN g;
454 2139790 : if (!(flag & (LLL_IM|LLL_KER|LLL_ALL|LLL_INPLACE))) return h;
455 2066739 : if (flag & (LLL_IM|LLL_INPLACE)) return vectail_inplace(h, k);
456 96 : if (flag & LLL_KER) { setlg(h,k+1); return h; }
457 68 : g = vecslice(h,1,k); /* done first: vectail_inplace kills h */
458 70 : return mkvec2(g, vectail_inplace(h, k));
459 : }
460 :
461 : /* y * z * 2^e, e >= 0; y,z t_INT */
462 : INLINE GEN
463 440248 : mulshift(GEN y, GEN z, long e)
464 : {
465 440248 : long ly = lgefint(y), lz;
466 : pari_sp av;
467 : GEN t;
468 440248 : if (ly == 2) return gen_0;
469 55256 : lz = lgefint(z);
470 55256 : av = avma; (void)new_chunk(ly+lz+nbits2lg(e)); /* HACK */
471 55256 : t = mulii(z, y);
472 55256 : set_avma(av); return shifti(t, e);
473 : }
474 :
475 : /* x - y * z * 2^e, e >= 0; x,y,z t_INT */
476 : INLINE GEN
477 2075086 : submulshift(GEN x, GEN y, GEN z, long e)
478 : {
479 2075086 : long lx = lgefint(x), ly, lz;
480 : pari_sp av;
481 : GEN t;
482 2075086 : if (!e) return submulii(x, y, z);
483 2055591 : if (lx == 2) { t = mulshift(y, z, e); togglesign(t); return t; }
484 1639513 : ly = lgefint(y);
485 1639513 : if (ly == 2) return icopy(x);
486 1184972 : lz = lgefint(z);
487 1184972 : av = avma; (void)new_chunk(lx+ly+lz+nbits2lg(e)); /* HACK */
488 1184972 : t = shifti(mulii(z, y), e);
489 1184972 : set_avma(av); return subii(x, t);
490 : }
491 : static void
492 18625350 : subzi(GEN *a, GEN b)
493 : {
494 18625350 : pari_sp av = avma;
495 18625350 : b = subii(*a, b);
496 18625317 : if (lgefint(b)<=lg(*a) && isonstack(*a)) { affii(b,*a); set_avma(av); }
497 2127896 : else *a = b;
498 18625352 : }
499 :
500 : static void
501 17910660 : addzi(GEN *a, GEN b)
502 : {
503 17910660 : pari_sp av = avma;
504 17910660 : b = addii(*a, b);
505 17910604 : if (lgefint(b)<=lg(*a) && isonstack(*a)) { affii(b,*a); set_avma(av); }
506 1909051 : else *a = b;
507 17910656 : }
508 :
509 : /* x - u*y * 2^e */
510 : INLINE GEN
511 4219432 : submuliu2n(GEN x, GEN y, ulong u, long e)
512 : {
513 : pari_sp av;
514 4219432 : long ly = lgefint(y);
515 4219432 : if (ly == 2) return x;
516 2863369 : av = avma;
517 2863369 : (void)new_chunk(3+ly+lgefint(x)+nbits2lg(e)); /* HACK */
518 2864443 : y = shifti(mului(u,y), e);
519 2863762 : set_avma(av); return subii(x, y);
520 : }
521 : /* *x -= u*y * 2^e */
522 : INLINE void
523 320022 : submulzu2n(GEN *x, GEN y, ulong u, long e)
524 : {
525 : pari_sp av;
526 320022 : long ly = lgefint(y);
527 320022 : if (ly == 2) return;
528 209921 : av = avma;
529 209921 : (void)new_chunk(3+ly+lgefint(*x)+nbits2lg(e)); /* HACK */
530 209921 : y = shifti(mului(u,y), e);
531 209921 : set_avma(av); return subzi(x, y);
532 : }
533 :
534 : /* x + u*y * 2^e */
535 : INLINE GEN
536 4099808 : addmuliu2n(GEN x, GEN y, ulong u, long e)
537 : {
538 : pari_sp av;
539 4099808 : long ly = lgefint(y);
540 4099808 : if (ly == 2) return x;
541 2798131 : av = avma;
542 2798131 : (void)new_chunk(3+ly+lgefint(x)+nbits2lg(e)); /* HACK */
543 2799046 : y = shifti(mului(u,y), e);
544 2798463 : set_avma(av); return addii(x, y);
545 : }
546 :
547 : /* *x += u*y * 2^e */
548 : INLINE void
549 326940 : addmulzu2n(GEN *x, GEN y, ulong u, long e)
550 : {
551 : pari_sp av;
552 326940 : long ly = lgefint(y);
553 326940 : if (ly == 2) return;
554 212692 : av = avma;
555 212692 : (void)new_chunk(3+ly+lgefint(*x)+nbits2lg(e)); /* HACK */
556 212692 : y = shifti(mului(u,y), e);
557 212692 : set_avma(av); return addzi(x, y);
558 : }
559 :
560 : /* n < 10; gerepileall supporting &NULL arguments. Maybe rename and export ? */
561 : INLINE void
562 4640 : gc_lll(pari_sp av, int n, ...)
563 : {
564 : int i, j;
565 : GEN *gptr[10];
566 : size_t s;
567 4640 : va_list a; va_start(a, n);
568 13920 : for (i=j=0; i<n; i++)
569 : {
570 9280 : GEN *x = va_arg(a,GEN*);
571 9280 : if (*x) { gptr[j++] = x; *x = (GEN)copy_bin(*x); }
572 : }
573 4640 : va_end(a); set_avma(av);
574 11539 : for (--j; j>=0; j--) *gptr[j] = bin_copy((GENbin*)*gptr[j]);
575 4640 : s = pari_mainstack->top - pari_mainstack->bot;
576 : /* size of saved objects ~ stacksize / 4 => overflow */
577 4640 : if (av - avma > (s >> 2))
578 : {
579 0 : size_t t = avma - pari_mainstack->bot;
580 0 : av = avma; new_chunk((s + t) / sizeof(long)); set_avma(av); /* double */
581 : }
582 4640 : }
583 :
584 : /********************************************************************/
585 : /** **/
586 : /** FPLLL (adapted from D. Stehle's code) **/
587 : /** **/
588 : /********************************************************************/
589 : /* Babai* and fplll* are a conversion to libpari API and data types
590 : of fplll-1.3 by Damien Stehle'.
591 :
592 : Copyright 2005, 2006 Damien Stehle'.
593 :
594 : This program is free software; you can redistribute it and/or modify it
595 : under the terms of the GNU General Public License as published by the
596 : Free Software Foundation; either version 2 of the License, or (at your
597 : option) any later version.
598 :
599 : This program implements ideas from the paper "Floating-point LLL Revisited",
600 : by Phong Nguyen and Damien Stehle', in the Proceedings of Eurocrypt'2005,
601 : Springer-Verlag; and was partly inspired by Shoup's NTL library:
602 : http://www.shoup.net/ntl/ */
603 :
604 : /* x t_REAL, |x| >= 1/2. Test whether |x| <= 3/2 */
605 : static int
606 469696 : absrsmall2(GEN x)
607 : {
608 469696 : long e = expo(x), l, i;
609 469696 : if (e < 0) return 1;
610 224807 : if (e > 0 || (ulong)x[2] > (3UL << (BITS_IN_LONG-2))) return 0;
611 : /* line above assumes l > 2. OK since x != 0 */
612 77329 : l = lg(x); for (i = 3; i < l; i++) if (x[i]) return 0;
613 66467 : return 1;
614 : }
615 : /* x t_REAL; test whether |x| <= 1/2 */
616 : static int
617 883924 : absrsmall(GEN x)
618 : {
619 : long e, l, i;
620 883924 : if (!signe(x)) return 1;
621 878279 : e = expo(x); if (e < -1) return 1;
622 476293 : if (e > -1 || (ulong)x[2] > HIGHBIT) return 0;
623 7528 : l = lg(x); for (i = 3; i < l; i++) if (x[i]) return 0;
624 6597 : return 1;
625 : }
626 :
627 : static void
628 31941614 : rotate(GEN A, long k2, long k)
629 : {
630 : long i;
631 31941614 : GEN B = gel(A,k2);
632 100992089 : for (i = k2; i > k; i--) gel(A,i) = gel(A,i-1);
633 31941614 : gel(A,k) = B;
634 31941614 : }
635 :
636 : /************************* FAST version (double) ************************/
637 : #define dmael(x,i,j) ((x)[i][j])
638 : #define del(x,i) ((x)[i])
639 :
640 : static double *
641 34677049 : cget_dblvec(long d)
642 34677049 : { return (double*) stack_malloc_align(d*sizeof(double), sizeof(double)); }
643 :
644 : static double **
645 8316981 : cget_dblmat(long d) { return (double **) cgetg(d, t_VECSMALL); }
646 :
647 : static double
648 159954959 : itodbl_exp(GEN x, long *e)
649 : {
650 159954959 : pari_sp av = avma;
651 159954959 : GEN r = itor(x,DEFAULTPREC);
652 159918607 : *e = expo(r); setexpo(r,0);
653 159915341 : return gc_double(av, rtodbl(r));
654 : }
655 :
656 : static double
657 116302375 : dbldotproduct(double *x, double *y, long n)
658 : {
659 : long i;
660 116302375 : double sum = del(x,1) * del(y,1);
661 1366193227 : for (i=2; i<=n; i++) sum += del(x,i) * del(y,i);
662 116302375 : return sum;
663 : }
664 :
665 : static double
666 2448893 : dbldotsquare(double *x, long n)
667 : {
668 : long i;
669 2448893 : double sum = del(x,1) * del(x,1);
670 8145565 : for (i=2; i<=n; i++) sum += del(x,i) * del(x,i);
671 2448893 : return sum;
672 : }
673 :
674 : static long
675 24728890 : set_line(double *appv, GEN v, long n)
676 : {
677 24728890 : long i, maxexp = 0;
678 24728890 : pari_sp av = avma;
679 24728890 : GEN e = cgetg(n+1, t_VECSMALL);
680 184671933 : for (i = 1; i <= n; i++)
681 : {
682 159953731 : del(appv,i) = itodbl_exp(gel(v,i), e+i);
683 159942325 : if (e[i] > maxexp) maxexp = e[i];
684 : }
685 184744462 : for (i = 1; i <= n; i++) del(appv,i) = ldexp(del(appv,i), e[i]-maxexp);
686 24718202 : set_avma(av); return maxexp;
687 : }
688 :
689 : static void
690 34590053 : dblrotate(double **A, long k2, long k)
691 : {
692 : long i;
693 34590053 : double *B = del(A,k2);
694 108466723 : for (i = k2; i > k; i--) del(A,i) = del(A,i-1);
695 34590053 : del(A,k) = B;
696 34590053 : }
697 : /* update G[kappa][i] from appB */
698 : static void
699 22490725 : setG_fast(double **appB, long n, double **G, long kappa, long a, long b)
700 : { long i;
701 100614357 : for (i = a; i <= b; i++)
702 78123805 : dmael(G,kappa,i) = dbldotproduct(del(appB,kappa), del(appB,i), n);
703 22490552 : }
704 : /* update G[i][kappa] from appB */
705 : static void
706 16876727 : setG2_fast(double **appB, long n, double **G, long kappa, long a, long b)
707 : { long i;
708 55058974 : for (i = a; i <= b; i++)
709 38182276 : dmael(G,i,kappa) = dbldotproduct(del(appB,kappa), del(appB,i), n);
710 16876698 : }
711 : const long EX0 = -2; /* uninitialized; any value less than expo(0.51) = -1 */
712 :
713 : #ifdef LONG_IS_64BIT
714 : typedef long s64;
715 : #define addmuliu64_inplace addmuliu_inplace
716 : #define submuliu64_inplace submuliu_inplace
717 : #define submuliu642n submuliu2n
718 : #define addmuliu642n addmuliu2n
719 : #else
720 : typedef long long s64;
721 : typedef unsigned long long u64;
722 :
723 : INLINE GEN
724 20062764 : u64toi(u64 x)
725 : {
726 : GEN y;
727 : ulong h;
728 20062764 : if (!x) return gen_0;
729 20062764 : h = x>>32;
730 20062764 : if (!h) return utoipos(x);
731 1149255 : y = cgetipos(4);
732 1149255 : *int_LSW(y) = x&0xFFFFFFFF;
733 1149255 : *int_MSW(y) = x>>32;
734 1149255 : return y;
735 : }
736 :
737 : INLINE GEN
738 665952 : u64toineg(u64 x)
739 : {
740 : GEN y;
741 : ulong h;
742 665952 : if (!x) return gen_0;
743 665952 : h = x>>32;
744 665952 : if (!h) return utoineg(x);
745 665952 : y = cgetineg(4);
746 665952 : *int_LSW(y) = x&0xFFFFFFFF;
747 665952 : *int_MSW(y) = x>>32;
748 665952 : return y;
749 : }
750 : INLINE GEN
751 9641729 : addmuliu64_inplace(GEN x, GEN y, u64 u) { return addmulii(x, y, u64toi(u)); }
752 :
753 : INLINE GEN
754 9732575 : submuliu64_inplace(GEN x, GEN y, u64 u) { return submulii(x, y, u64toi(u)); }
755 :
756 : INLINE GEN
757 665952 : addmuliu642n(GEN x, GEN y, u64 u, long e) { return submulshift(x, y, u64toineg(u), e); }
758 :
759 : INLINE GEN
760 688460 : submuliu642n(GEN x, GEN y, u64 u, long e) { return submulshift(x, y, u64toi(u), e); }
761 :
762 : #endif
763 :
764 : /* Babai's Nearest Plane algorithm (iterative); see Babai() */
765 : static int
766 29800525 : Babai_fast(pari_sp av, long kappa, GEN *pB, GEN *pU, double **mu, double **r,
767 : double *s, double **appB, GEN expoB, double **G,
768 : long a, long zeros, long maxG, double eta)
769 : {
770 29800525 : GEN B = *pB, U = *pU;
771 29800525 : const long n = nbrows(B), d = U ? lg(U)-1: 0;
772 29800349 : long k, aa = (a > zeros)? a : zeros+1;
773 29800349 : long emaxmu = EX0, emax2mu = EX0;
774 : s64 xx;
775 29800349 : int did_something = 0;
776 : /* N.B: we set d = 0 (resp. n = 0) to avoid updating U (resp. B) */
777 :
778 17101269 : for (;;) {
779 46901618 : int go_on = 0;
780 46901618 : long i, j, emax3mu = emax2mu;
781 :
782 46901618 : if (gc_needed(av,2))
783 : {
784 86 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[1], a=%ld", aa);
785 86 : gc_lll(av,2,&B,&U);
786 : }
787 : /* Step2: compute the GSO for stage kappa */
788 46900746 : emax2mu = emaxmu; emaxmu = EX0;
789 178474226 : for (j=aa; j<kappa; j++)
790 : {
791 131574894 : double g = dmael(G,kappa,j);
792 567063669 : for (k = zeros+1; k < j; k++) g -= dmael(mu,j,k) * dmael(r,kappa,k);
793 131574894 : dmael(r,kappa,j) = g;
794 131574894 : dmael(mu,kappa,j) = dmael(r,kappa,j) / dmael(r,j,j);
795 131574894 : emaxmu = maxss(emaxmu, expoB[kappa]-expoB[j]);
796 : }
797 : /* maxmu doesn't decrease fast enough */
798 46899332 : if (emax3mu != EX0 && emax3mu <= emax2mu + 5) {*pB = B; *pU = U; return 1;}
799 :
800 165627371 : for (j=kappa-1; j>zeros; j--)
801 : {
802 135833504 : double tmp = fabs(ldexp (dmael(mu,kappa,j), expoB[kappa]-expoB[j]));
803 135833504 : if (tmp>eta) { go_on = 1; break; }
804 : }
805 :
806 : /* Step3--5: compute the X_j's */
807 46895078 : if (go_on)
808 76859960 : for (j=kappa-1; j>zeros; j--)
809 : { /* The code below seemingly handles U = NULL, but in this case d = 0 */
810 59773375 : int e = expoB[j] - expoB[kappa];
811 59773375 : double tmp = ldexp(dmael(mu,kappa,j), -e), atmp = fabs(tmp);
812 : /* tmp = Inf is allowed */
813 59773375 : if (atmp <= .5) continue; /* size-reduced */
814 33702639 : if (gc_needed(av,2))
815 : {
816 241 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[2], a=%ld, j=%ld", aa,j);
817 241 : gc_lll(av,2,&B,&U);
818 : }
819 33704363 : did_something = 1;
820 : /* we consider separately the case |X| = 1 */
821 33704363 : if (atmp <= 1.5)
822 : {
823 22395959 : if (dmael(mu,kappa,j) > 0) { /* in this case, X = 1 */
824 46159086 : for (k=zeros+1; k<j; k++)
825 34730505 : dmael(mu,kappa,k) -= ldexp(dmael(mu,j,k), e);
826 155442529 : for (i=1; i<=n; i++)
827 144021629 : gmael(B,kappa,i) = subii(gmael(B,kappa,i), gmael(B,j,i));
828 102405256 : for (i=1; i<=d; i++)
829 90984147 : gmael(U,kappa,i) = subii(gmael(U,kappa,i), gmael(U,j,i));
830 : } else { /* otherwise X = -1 */
831 45328449 : for (k=zeros+1; k<j; k++)
832 34361071 : dmael(mu,kappa,k) += ldexp(dmael(mu,j,k), e);
833 152722299 : for (i=1; i<=n; i++)
834 141762870 : gmael(B,kappa,i) = addii(gmael(B,kappa,i), gmael(B,j,i));
835 99800182 : for (i=1; i<=d; i++)
836 88840642 : gmael(U,kappa,i) = addii(gmael(U,kappa,i), gmael(U,j,i));
837 : }
838 22380649 : continue;
839 : }
840 : /* we have |X| >= 2 */
841 11308404 : if (atmp < 9007199254740992.)
842 : {
843 10456343 : tmp = pari_rint(tmp);
844 24711201 : for (k=zeros+1; k<j; k++)
845 14254857 : dmael(mu,kappa,k) -= ldexp(tmp * dmael(mu,j,k), e);
846 10456344 : xx = (s64) tmp;
847 10456344 : if (xx > 0) /* = xx */
848 : {
849 46134099 : for (i=1; i<=n; i++)
850 40883280 : gmael(B,kappa,i) = submuliu64_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
851 33973324 : for (i=1; i<=d; i++)
852 28722452 : gmael(U,kappa,i) = submuliu64_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
853 : }
854 : else /* = -xx */
855 : {
856 45880253 : for (i=1; i<=n; i++)
857 40675586 : gmael(B,kappa,i) = addmuliu64_inplace(gmael(B,kappa,i), gmael(B,j,i), -xx);
858 33545857 : for (i=1; i<=d; i++)
859 28341257 : gmael(U,kappa,i) = addmuliu64_inplace(gmael(U,kappa,i), gmael(U,j,i), -xx);
860 : }
861 : }
862 : else
863 : {
864 : int E;
865 852061 : xx = (s64) ldexp(frexp(dmael(mu,kappa,j), &E), 53);
866 852061 : E -= e + 53;
867 852061 : if (E <= 0)
868 : {
869 0 : xx = xx << -E;
870 0 : for (k=zeros+1; k<j; k++)
871 0 : dmael(mu,kappa,k) -= ldexp(((double)xx) * dmael(mu,j,k), e);
872 0 : if (xx > 0) /* = xx */
873 : {
874 0 : for (i=1; i<=n; i++)
875 0 : gmael(B,kappa,i) = submuliu64_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
876 0 : for (i=1; i<=d; i++)
877 0 : gmael(U,kappa,i) = submuliu64_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
878 : }
879 : else /* = -xx */
880 : {
881 0 : for (i=1; i<=n; i++)
882 0 : gmael(B,kappa,i) = addmuliu64_inplace(gmael(B,kappa,i), gmael(B,j,i), -xx);
883 0 : for (i=1; i<=d; i++)
884 0 : gmael(U,kappa,i) = addmuliu64_inplace(gmael(U,kappa,i), gmael(U,j,i), -xx);
885 : }
886 : } else
887 : {
888 2775130 : for (k=zeros+1; k<j; k++)
889 1923069 : dmael(mu,kappa,k) -= ldexp(((double)xx) * dmael(mu,j,k), E + e);
890 852061 : if (xx > 0) /* = xx */
891 : {
892 3937643 : for (i=1; i<=n; i++)
893 3508155 : gmael(B,kappa,i) = submuliu642n(gmael(B,kappa,i), gmael(B,j,i), xx, E);
894 1569305 : for (i=1; i<=d; i++)
895 1139817 : gmael(U,kappa,i) = submuliu642n(gmael(U,kappa,i), gmael(U,j,i), xx, E);
896 : }
897 : else /* = -xx */
898 : {
899 3870304 : for (i=1; i<=n; i++)
900 3447906 : gmael(B,kappa,i) = addmuliu642n(gmael(B,kappa,i), gmael(B,j,i), -xx, E);
901 1524497 : for (i=1; i<=d; i++)
902 1102099 : gmael(U,kappa,i) = addmuliu642n(gmael(U,kappa,i), gmael(U,j,i), -xx, E);
903 : }
904 : }
905 : }
906 : }
907 46880446 : if (!go_on) break; /* Anything happened? */
908 17084485 : expoB[kappa] = set_line(del(appB,kappa), gel(B,kappa), n);
909 17100829 : setG_fast(appB, n, G, kappa, zeros+1, kappa-1);
910 17101269 : aa = zeros+1;
911 : }
912 29795961 : if (did_something) setG2_fast(appB, n, G, kappa, kappa, maxG);
913 :
914 29796358 : del(s,zeros+1) = dmael(G,kappa,kappa);
915 : /* the last s[kappa-1]=r[kappa][kappa] is computed only if kappa increases */
916 107758320 : for (k=zeros+1; k<=kappa-2; k++)
917 77961962 : del(s,k+1) = del(s,k) - dmael(mu,kappa,k)*dmael(r,kappa,k);
918 29796358 : *pB = B; *pU = U; return 0;
919 : }
920 :
921 : static void
922 11977471 : update_alpha(GEN alpha, long kappa, long kappa2, long kappamax)
923 : {
924 : long i;
925 37722061 : for (i = kappa; i < kappa2; i++)
926 25744590 : if (kappa <= alpha[i]) alpha[i] = kappa;
927 37722128 : for (i = kappa2; i > kappa; i--) alpha[i] = alpha[i-1];
928 22973254 : for (i = kappa2+1; i <= kappamax; i++)
929 10995783 : if (kappa < alpha[i]) alpha[i] = kappa;
930 11977471 : alpha[kappa] = kappa;
931 11977471 : }
932 : static void
933 447262 : rotateG(GEN G, long kappa2, long kappa, long maxG, GEN Gtmp)
934 : {
935 : long i, j;
936 3570201 : for (i=1; i<=kappa2; i++) gel(Gtmp,i) = gmael(G,kappa2,i);
937 1874447 : for ( ; i<=maxG; i++) gel(Gtmp,i) = gmael(G,i,kappa2);
938 1565822 : for (i=kappa2; i>kappa; i--)
939 : {
940 5648759 : for (j=1; j<kappa; j++) gmael(G,i,j) = gmael(G,i-1,j);
941 1118560 : gmael(G,i,kappa) = gel(Gtmp,i-1);
942 4082653 : for (j=kappa+1; j<=i; j++) gmael(G,i,j) = gmael(G,i-1,j-1);
943 4920736 : for (j=kappa2+1; j<=maxG; j++) gmael(G,j,i) = gmael(G,j,i-1);
944 : }
945 2004379 : for (i=1; i<kappa; i++) gmael(G,kappa,i) = gel(Gtmp,i);
946 447262 : gmael(G,kappa,kappa) = gel(Gtmp,kappa2);
947 1874447 : for (i=kappa2+1; i<=maxG; i++) gmael(G,i,kappa) = gel(Gtmp,i);
948 447262 : }
949 : static void
950 11530075 : rotateG_fast(double **G, long kappa2, long kappa, long maxG, double *Gtmp)
951 : {
952 : long i, j;
953 66673647 : for (i=1; i<=kappa2; i++) del(Gtmp,i) = dmael(G,kappa2,i);
954 21972471 : for ( ; i<=maxG; i++) del(Gtmp,i) = dmael(G,i,kappa2);
955 36155856 : for (i=kappa2; i>kappa; i--)
956 : {
957 69344271 : for (j=1; j<kappa; j++) dmael(G,i,j) = dmael(G,i-1,j);
958 24625781 : dmael(G,i,kappa) = del(Gtmp,i-1);
959 83325820 : for (j=kappa+1; j<=i; j++) dmael(G,i,j) = dmael(G,i-1,j-1);
960 46141549 : for (j=kappa2+1; j<=maxG; j++) dmael(G,j,i) = dmael(G,j,i-1);
961 : }
962 30518093 : for (i=1; i<kappa; i++) dmael(G,kappa,i) = del(Gtmp,i);
963 11530075 : dmael(G,kappa,kappa) = del(Gtmp,kappa2);
964 21972493 : for (i=kappa2+1; i<=maxG; i++) dmael(G,i,kappa) = del(Gtmp,i);
965 11530075 : }
966 :
967 : /* LLL-reduces (B,U) in place [apply base change transforms to B and U].
968 : * Gram matrix, and GSO performed on matrices of 'double'.
969 : * If (keepfirst), never swap with first vector.
970 : * Return -1 on failure, else zeros = dim Kernel (>= 0) */
971 : static long
972 2079256 : fplll_fast(GEN *pB, GEN *pU, double delta, double eta, long keepfirst)
973 : {
974 : pari_sp av;
975 : long kappa, kappa2, d, n, i, j, zeros, kappamax, maxG;
976 : double **mu, **r, *s, tmp, *Gtmp, **G, **appB;
977 2079256 : GEN alpha, expoB, B = *pB, U;
978 2079256 : long cnt = 0;
979 :
980 2079256 : d = lg(B)-1;
981 2079256 : n = nbrows(B);
982 2079255 : U = *pU; /* NULL if inplace */
983 :
984 2079255 : G = cget_dblmat(d+1);
985 2079254 : appB = cget_dblmat(d+1);
986 2079254 : mu = cget_dblmat(d+1);
987 2079254 : r = cget_dblmat(d+1);
988 2079255 : s = cget_dblvec(d+1);
989 9709002 : for (j = 1; j <= d; j++)
990 : {
991 7629744 : del(mu,j) = cget_dblvec(d+1);
992 7629748 : del(r,j) = cget_dblvec(d+1);
993 7629754 : del(appB,j) = cget_dblvec(n+1);
994 7629746 : del(G,j) = cget_dblvec(d+1);
995 47756580 : for (i=1; i<=d; i++) dmael(G,j,i) = 0.;
996 : }
997 2079258 : expoB = cgetg(d+1, t_VECSMALL);
998 9708948 : for (i=1; i<=d; i++) expoB[i] = set_line(del(appB,i), gel(B,i), n);
999 2079226 : Gtmp = cget_dblvec(d+1);
1000 2079252 : alpha = cgetg(d+1, t_VECSMALL);
1001 2079247 : av = avma;
1002 :
1003 : /* Step2: Initializing the main loop */
1004 2079247 : kappamax = 1;
1005 2079247 : i = 1;
1006 2079247 : maxG = d; /* later updated to kappamax */
1007 :
1008 : do {
1009 2235718 : dmael(G,i,i) = dbldotsquare(del(appB,i),n);
1010 2235723 : } while (dmael(G,i,i) <= 0 && (++i <=d));
1011 2079252 : zeros = i-1; /* all vectors B[i] with i <= zeros are zero vectors */
1012 2079252 : kappa = i;
1013 2079252 : if (zeros < d) dmael(r,zeros+1,zeros+1) = dmael(G,zeros+1,zeros+1);
1014 9552537 : for (i=zeros+1; i<=d; i++) alpha[i]=1;
1015 31875562 : while (++kappa <= d)
1016 : {
1017 29800555 : if (kappa > kappamax)
1018 : {
1019 5389897 : if (DEBUGLEVEL>=4) err_printf("K%ld ",kappa);
1020 5389897 : maxG = kappamax = kappa;
1021 5389897 : setG_fast(appB, n, G, kappa, zeros+1, kappa);
1022 : }
1023 : /* Step3: Call to the Babai algorithm, mu,r,s updated in place */
1024 29800562 : if (Babai_fast(av, kappa, &B,&U, mu,r,s, appB, expoB, G, alpha[kappa],
1025 4254 : zeros, maxG, eta)) { *pB=B; *pU=U; return -1; }
1026 :
1027 29796244 : tmp = ldexp(r[kappa-1][kappa-1] * delta, 2*(expoB[kappa-1]-expoB[kappa]));
1028 29796244 : if ((keepfirst && kappa == 2) || tmp <= del(s,kappa-1))
1029 : { /* Step4: Success of Lovasz's condition */
1030 18266036 : alpha[kappa] = kappa;
1031 18266036 : tmp = dmael(mu,kappa,kappa-1) * dmael(r,kappa,kappa-1);
1032 18266036 : dmael(r,kappa,kappa) = del(s,kappa-1)- tmp;
1033 18266036 : continue;
1034 : }
1035 : /* Step5: Find the right insertion index kappa, kappa2 = initial kappa */
1036 11530208 : if (DEBUGLEVEL>=4 && kappa==kappamax && del(s,kappa-1)!=0)
1037 0 : if (++cnt > 20) { cnt = 0; err_printf("(%ld) ", 2*expoB[1] + dblexpo(del(s,1))); }
1038 11530200 : kappa2 = kappa;
1039 : do {
1040 24625995 : kappa--;
1041 24625995 : if (kappa<zeros+2 + (keepfirst ? 1: 0)) break;
1042 17949592 : tmp = dmael(r,kappa-1,kappa-1) * delta;
1043 17949592 : tmp = ldexp(tmp, 2*(expoB[kappa-1]-expoB[kappa2]));
1044 17949592 : } while (del(s,kappa-1) <= tmp);
1045 11530200 : update_alpha(alpha, kappa, kappa2, kappamax);
1046 :
1047 : /* Step6: Update the mu's and r's */
1048 11530223 : dblrotate(mu,kappa2,kappa);
1049 11530182 : dblrotate(r,kappa2,kappa);
1050 11530134 : dmael(r,kappa,kappa) = del(s,kappa);
1051 :
1052 : /* Step7: Update B, appB, U, G */
1053 11530134 : rotate(B,kappa2,kappa);
1054 11530126 : dblrotate(appB,kappa2,kappa);
1055 11530100 : if (U) rotate(U,kappa2,kappa);
1056 11530099 : rotate(expoB,kappa2,kappa);
1057 11530083 : rotateG_fast(G,kappa2,kappa, maxG, Gtmp);
1058 :
1059 : /* Step8: Prepare the next loop iteration */
1060 11530274 : if (kappa == zeros+1 && dmael(G,kappa,kappa)<= 0)
1061 : {
1062 213175 : zeros++; kappa++;
1063 213175 : dmael(G,kappa,kappa) = dbldotsquare(del(appB,kappa),n);
1064 213175 : dmael(r,kappa,kappa) = dmael(G,kappa,kappa);
1065 : }
1066 : }
1067 2075007 : *pB = B; *pU = U; return zeros;
1068 : }
1069 :
1070 : /***************** HEURISTIC version (reduced precision) ****************/
1071 : static GEN
1072 186486 : realsqrdotproduct(GEN x)
1073 : {
1074 186486 : long i, l = lg(x);
1075 186486 : GEN z = sqrr(gel(x,1));
1076 1474993 : for (i=2; i<l; i++) z = addrr(z, sqrr(gel(x,i)));
1077 186486 : return z;
1078 : }
1079 : /* x, y non-empty vector of t_REALs, same length */
1080 : static GEN
1081 1397760 : realdotproduct(GEN x, GEN y)
1082 : {
1083 : long i, l;
1084 : GEN z;
1085 1397760 : if (x == y) return realsqrdotproduct(x);
1086 1211274 : l = lg(x); z = mulrr(gel(x,1),gel(y,1));
1087 14242421 : for (i=2; i<l; i++) z = addrr(z, mulrr(gel(x,i), gel(y,i)));
1088 1211274 : return z;
1089 : }
1090 : static void
1091 198575 : setG_heuristic(GEN appB, GEN G, long kappa, long a, long b)
1092 198575 : { pari_sp av = avma;
1093 : long i;
1094 1121435 : for (i = a; i <= b; i++)
1095 922860 : affrr(realdotproduct(gel(appB,kappa),gel(appB,i)), gmael(G,kappa,i));
1096 198575 : set_avma(av);
1097 198575 : }
1098 : static void
1099 175753 : setG2_heuristic(GEN appB, GEN G, long kappa, long a, long b)
1100 175753 : { pari_sp av = avma;
1101 : long i;
1102 650653 : for (i = a; i <= b; i++)
1103 474900 : affrr(realdotproduct(gel(appB,kappa),gel(appB,i)), gmael(G,i,kappa));
1104 175753 : set_avma(av);
1105 175753 : }
1106 :
1107 : /* approximate t_REAL x as m * 2^e, where |m| < 2^bit */
1108 : static GEN
1109 32828 : truncexpo(GEN x, long bit, long *e)
1110 : {
1111 32828 : *e = expo(x) + 1 - bit;
1112 32828 : if (*e >= 0) return mantissa2nr(x, 0);
1113 1156 : *e = 0; return roundr_safe(x);
1114 : }
1115 : /* Babai's Nearest Plane algorithm (iterative); see Babai() */
1116 : static int
1117 285708 : Babai_heuristic(pari_sp av, long kappa, GEN *pB, GEN *pU, GEN mu, GEN r, GEN s,
1118 : GEN appB, GEN G, long a, long zeros, long maxG,
1119 : GEN eta, long prec)
1120 : {
1121 285708 : GEN B = *pB, U = *pU;
1122 285708 : const long n = nbrows(B), d = U ? lg(U)-1: 0, bit = prec2nbits(prec);
1123 285708 : long k, aa = (a > zeros)? a : zeros+1;
1124 285708 : int did_something = 0;
1125 285708 : long emaxmu = EX0, emax2mu = EX0;
1126 : /* N.B: we set d = 0 (resp. n = 0) to avoid updating U (resp. B) */
1127 :
1128 187842 : for (;;) {
1129 473550 : int go_on = 0;
1130 473550 : long i, j, emax3mu = emax2mu;
1131 :
1132 473550 : if (gc_needed(av,2))
1133 : {
1134 25 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[1], a=%ld", aa);
1135 25 : gc_lll(av,2,&B,&U);
1136 : }
1137 : /* Step2: compute the GSO for stage kappa */
1138 473550 : emax2mu = emaxmu; emaxmu = EX0;
1139 2100814 : for (j=aa; j<kappa; j++)
1140 : {
1141 1627264 : pari_sp btop = avma;
1142 1627264 : GEN g = gmael(G,kappa,j);
1143 7831794 : for (k = zeros+1; k<j; k++)
1144 6204530 : g = subrr(g, mulrr(gmael(mu,j,k), gmael(r,kappa,k)));
1145 1627264 : affrr(g, gmael(r,kappa,j));
1146 1627264 : affrr(divrr(gmael(r,kappa,j), gmael(r,j,j)), gmael(mu,kappa,j));
1147 1627264 : emaxmu = maxss(emaxmu, expo(gmael(mu,kappa,j)));
1148 1627264 : set_avma(btop);
1149 : }
1150 473550 : if (emax3mu != EX0 && emax3mu <= emax2mu + 5)
1151 1334 : { *pB = B; *pU = U; return 1; }
1152 :
1153 2078167 : for (j=kappa-1; j>zeros; j--)
1154 1793793 : if (abscmprr(gmael(mu,kappa,j), eta) > 0) { go_on = 1; break; }
1155 :
1156 : /* Step3--5: compute the X_j's */
1157 472216 : if (go_on)
1158 1071766 : for (j=kappa-1; j>zeros; j--)
1159 : { /* The code below seemingly handles U = NULL, but in this case d = 0 */
1160 : pari_sp btop;
1161 883924 : GEN tmp = gmael(mu,kappa,j);
1162 883924 : if (absrsmall(tmp)) continue; /* size-reduced */
1163 :
1164 469696 : if (gc_needed(av,2))
1165 : {
1166 34 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[2], a=%ld, j=%ld", aa,j);
1167 34 : gc_lll(av,2,&B,&U);
1168 : }
1169 469696 : btop = avma; did_something = 1;
1170 : /* we consider separately the case |X| = 1 */
1171 469696 : if (absrsmall2(tmp))
1172 : {
1173 311356 : if (signe(tmp) > 0) { /* in this case, X = 1 */
1174 657803 : for (k=zeros+1; k<j; k++)
1175 502615 : affrr(subrr(gmael(mu,kappa,k), gmael(mu,j,k)), gmael(mu,kappa,k));
1176 155188 : set_avma(btop);
1177 1988579 : for (i=1; i<=n; i++)
1178 1833391 : gmael(B,kappa,i) = subii(gmael(B,kappa,i), gmael(B,j,i));
1179 1570446 : for (i=1; i<=d; i++)
1180 1415258 : gmael(U,kappa,i) = subii(gmael(U,kappa,i), gmael(U,j,i));
1181 : } else { /* otherwise X = -1 */
1182 665547 : for (k=zeros+1; k<j; k++)
1183 509379 : affrr(addrr(gmael(mu,kappa,k), gmael(mu,j,k)), gmael(mu,kappa,k));
1184 156168 : set_avma(btop);
1185 2005481 : for (i=1; i<=n; i++)
1186 1849313 : gmael(B,kappa,i) = addii(gmael(B,kappa,i), gmael(B,j,i));
1187 1574693 : for (i=1; i<=d; i++)
1188 1418525 : gmael(U,kappa,i) = addii(gmael(U,kappa,i),gmael(U,j,i));
1189 : }
1190 311356 : continue;
1191 : }
1192 : /* we have |X| >= 2 */
1193 158340 : if (expo(tmp) < BITS_IN_LONG)
1194 : {
1195 125512 : ulong xx = roundr_safe(tmp)[2]; /* X fits in an ulong */
1196 125512 : if (signe(tmp) > 0) /* = xx */
1197 : {
1198 152389 : for (k=zeros+1; k<j; k++)
1199 89272 : affrr(subrr(gmael(mu,kappa,k), mulur(xx, gmael(mu,j,k))),
1200 89272 : gmael(mu,kappa,k));
1201 63117 : set_avma(btop);
1202 463660 : for (i=1; i<=n; i++)
1203 400543 : gmael(B,kappa,i) = submuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
1204 359024 : for (i=1; i<=d; i++)
1205 295907 : gmael(U,kappa,i) = submuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
1206 : }
1207 : else /* = -xx */
1208 : {
1209 150016 : for (k=zeros+1; k<j; k++)
1210 87621 : affrr(addrr(gmael(mu,kappa,k), mulur(xx, gmael(mu,j,k))),
1211 87621 : gmael(mu,kappa,k));
1212 62395 : set_avma(btop);
1213 460046 : for (i=1; i<=n; i++)
1214 397651 : gmael(B,kappa,i) = addmuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
1215 349356 : for (i=1; i<=d; i++)
1216 286961 : gmael(U,kappa,i) = addmuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
1217 : }
1218 : }
1219 : else
1220 : {
1221 : long e;
1222 32828 : GEN X = truncexpo(tmp, bit, &e); /* tmp ~ X * 2^e */
1223 32828 : btop = avma;
1224 131412 : for (k=zeros+1; k<j; k++)
1225 : {
1226 98584 : GEN x = mulir(X, gmael(mu,j,k));
1227 98584 : if (e) shiftr_inplace(x, e);
1228 98584 : affrr(subrr(gmael(mu,kappa,k), x), gmael(mu,kappa,k));
1229 : }
1230 32828 : set_avma(btop);
1231 405295 : for (i=1; i<=n; i++)
1232 372467 : gmael(B,kappa,i) = submulshift(gmael(B,kappa,i), gmael(B,j,i), X, e);
1233 381035 : for (i=1; i<=d; i++)
1234 348207 : gmael(U,kappa,i) = submulshift(gmael(U,kappa,i), gmael(U,j,i), X, e);
1235 : }
1236 : }
1237 472216 : if (!go_on) break; /* Anything happened? */
1238 1681903 : for (i=1 ; i<=n; i++) affir(gmael(B,kappa,i), gmael(appB,kappa,i));
1239 187842 : setG_heuristic(appB, G, kappa, zeros+1, kappa-1);
1240 187842 : aa = zeros+1;
1241 : }
1242 284374 : if (did_something) setG2_heuristic(appB, G, kappa, kappa, maxG);
1243 284374 : affrr(gmael(G,kappa,kappa), gel(s,zeros+1));
1244 : /* the last s[kappa-1]=r[kappa][kappa] is computed only if kappa increases */
1245 284374 : av = avma;
1246 1458300 : for (k=zeros+1; k<=kappa-2; k++)
1247 1173926 : affrr(subrr(gel(s,k), mulrr(gmael(mu,kappa,k), gmael(r,kappa,k))),
1248 1173926 : gel(s,k+1));
1249 284374 : *pB = B; *pU = U; return gc_bool(av, 0);
1250 : }
1251 :
1252 : static GEN
1253 16931 : ZC_to_RC(GEN x, long prec)
1254 120832 : { pari_APPLY_type(t_COL,itor(gel(x,i),prec)) }
1255 :
1256 : static GEN
1257 4254 : ZM_to_RM(GEN x, long prec)
1258 21185 : { pari_APPLY_same(ZC_to_RC(gel(x,i),prec)) }
1259 :
1260 : /* LLL-reduces (B,U) in place [apply base change transforms to B and U].
1261 : * Gram matrix made of t_REAL at precision prec2, performe GSO at prec.
1262 : * If (keepfirst), never swap with first vector.
1263 : * Return -1 on failure, else zeros = dim Kernel (>= 0) */
1264 : static long
1265 4254 : fplll_heuristic(GEN *pB, GEN *pU, double DELTA, double ETA, long keepfirst,
1266 : long prec, long prec2)
1267 : {
1268 : pari_sp av, av2;
1269 : long kappa, kappa2, d, i, j, zeros, kappamax, maxG;
1270 4254 : GEN mu, r, s, tmp, Gtmp, alpha, G, appB, B = *pB, U;
1271 4254 : GEN delta = dbltor(DELTA), eta = dbltor(ETA);
1272 4254 : long cnt = 0;
1273 :
1274 4254 : d = lg(B)-1;
1275 4254 : U = *pU; /* NULL if inplace */
1276 :
1277 4254 : G = cgetg(d+1, t_MAT);
1278 4254 : mu = cgetg(d+1, t_MAT);
1279 4254 : r = cgetg(d+1, t_MAT);
1280 4254 : s = cgetg(d+1, t_VEC);
1281 4254 : appB = ZM_to_RM(B, prec2);
1282 21185 : for (j = 1; j <= d; j++)
1283 : {
1284 16931 : GEN M = cgetg(d+1, t_COL), R = cgetg(d+1, t_COL), S = cgetg(d+1, t_COL);
1285 16931 : gel(mu,j)= M;
1286 16931 : gel(r,j) = R;
1287 16931 : gel(G,j) = S;
1288 16931 : gel(s,j) = cgetr(prec);
1289 112508 : for (i = 1; i <= d; i++)
1290 : {
1291 95577 : gel(R,i) = cgetr(prec);
1292 95577 : gel(M,i) = cgetr(prec);
1293 95577 : gel(S,i) = cgetr(prec2);
1294 : }
1295 : }
1296 4254 : Gtmp = cgetg(d+1, t_VEC);
1297 4254 : alpha = cgetg(d+1, t_VECSMALL);
1298 4254 : av = avma;
1299 :
1300 : /* Step2: Initializing the main loop */
1301 4254 : kappamax = 1;
1302 4254 : i = 1;
1303 4254 : maxG = d; /* later updated to kappamax */
1304 :
1305 : do {
1306 4272 : affrr(RgV_dotsquare(gel(appB,i)), gmael(G,i,i));
1307 4272 : } while (signe(gmael(G,i,i)) == 0 && (++i <=d));
1308 4254 : zeros = i-1; /* all vectors B[i] with i <= zeros are zero vectors */
1309 4254 : kappa = i;
1310 4254 : if (zeros < d) affrr(gmael(G,zeros+1,zeros+1), gmael(r,zeros+1,zeros+1));
1311 21167 : for (i=zeros+1; i<=d; i++) alpha[i]=1;
1312 :
1313 288628 : while (++kappa <= d)
1314 : {
1315 285708 : if (kappa > kappamax)
1316 : {
1317 10733 : if (DEBUGLEVEL>=4) err_printf("K%ld ",kappa);
1318 10733 : maxG = kappamax = kappa;
1319 10733 : setG_heuristic(appB, G, kappa, zeros+1, kappa);
1320 : }
1321 : /* Step3: Call to the Babai algorithm, mu,r,s updated in place */
1322 285708 : if (Babai_heuristic(av, kappa, &B,&U, mu,r,s, appB, G, alpha[kappa], zeros,
1323 1334 : maxG, eta, prec)) { *pB = B; *pU = U; return -1; }
1324 284374 : av2 = avma;
1325 568662 : if ((keepfirst && kappa == 2) ||
1326 284288 : cmprr(mulrr(gmael(r,kappa-1,kappa-1), delta), gel(s,kappa-1)) <= 0)
1327 : { /* Step4: Success of Lovasz's condition */
1328 170045 : alpha[kappa] = kappa;
1329 170045 : tmp = mulrr(gmael(mu,kappa,kappa-1), gmael(r,kappa,kappa-1));
1330 170045 : affrr(subrr(gel(s,kappa-1), tmp), gmael(r,kappa,kappa));
1331 170045 : set_avma(av2); continue;
1332 : }
1333 : /* Step5: Find the right insertion index kappa, kappa2 = initial kappa */
1334 114329 : if (DEBUGLEVEL>=4 && kappa==kappamax && signe(gel(s,kappa-1)))
1335 0 : if (++cnt > 20) { cnt = 0; err_printf("(%ld) ", expo(gel(s,1))); }
1336 114329 : kappa2 = kappa;
1337 : do {
1338 274992 : kappa--;
1339 274992 : if (kappa < zeros+2 + (keepfirst ? 1: 0)) break;
1340 250815 : tmp = mulrr(gmael(r,kappa-1,kappa-1), delta);
1341 250815 : } while (cmprr(gel(s,kappa-1), tmp) <= 0 );
1342 114329 : set_avma(av2);
1343 114329 : update_alpha(alpha, kappa, kappa2, kappamax);
1344 :
1345 : /* Step6: Update the mu's and r's */
1346 114329 : rotate(mu,kappa2,kappa);
1347 114329 : rotate(r,kappa2,kappa);
1348 114329 : affrr(gel(s,kappa), gmael(r,kappa,kappa));
1349 :
1350 : /* Step7: Update B, appB, U, G */
1351 114329 : rotate(B,kappa2,kappa);
1352 114329 : rotate(appB,kappa2,kappa);
1353 114329 : if (U) rotate(U,kappa2,kappa);
1354 114329 : rotateG(G,kappa2,kappa, maxG, Gtmp);
1355 :
1356 : /* Step8: Prepare the next loop iteration */
1357 114329 : if (kappa == zeros+1 && !signe(gmael(G,kappa,kappa)))
1358 : {
1359 0 : zeros++; kappa++;
1360 0 : affrr(RgV_dotsquare(gel(appB,kappa)), gmael(G,kappa,kappa));
1361 0 : affrr(gmael(G,kappa,kappa), gmael(r,kappa,kappa));
1362 : }
1363 : }
1364 2920 : *pB=B; *pU=U; return zeros;
1365 : }
1366 :
1367 : /************************* PROVED version (t_INT) ***********************/
1368 : /* dpe inspired by dpe.h by Patrick Pelissier, Paul Zimmermann
1369 : * https://gforge.inria.fr/projects/dpe/
1370 : */
1371 :
1372 : typedef struct
1373 : {
1374 : double d; /* significand */
1375 : long e; /* exponent */
1376 : } dpe_t;
1377 :
1378 : #define Dmael(x,i,j) (&((x)[i][j]))
1379 : #define Del(x,i) (&((x)[i]))
1380 :
1381 : static void
1382 665866 : dperotate(dpe_t **A, long k2, long k)
1383 : {
1384 : long i;
1385 665866 : dpe_t *B = A[k2];
1386 2353002 : for (i = k2; i > k; i--) A[i] = A[i-1];
1387 665866 : A[k] = B;
1388 665866 : }
1389 :
1390 : static void
1391 108253158 : dpe_normalize0(dpe_t *x)
1392 : {
1393 : int e;
1394 108253158 : x->d = frexp(x->d, &e);
1395 108253158 : x->e += e;
1396 108253158 : }
1397 :
1398 : static void
1399 47943239 : dpe_normalize(dpe_t *x)
1400 : {
1401 47943239 : if (x->d == 0.0)
1402 500243 : x->e = -LONG_MAX;
1403 : else
1404 47442996 : dpe_normalize0(x);
1405 47943325 : }
1406 :
1407 : static GEN
1408 93410 : dpetor(dpe_t *x)
1409 : {
1410 93410 : GEN r = dbltor(x->d);
1411 93409 : if (signe(r)==0) return r;
1412 93409 : setexpo(r, x->e-1);
1413 93409 : return r;
1414 : }
1415 :
1416 : static void
1417 25742189 : affdpe(dpe_t *y, dpe_t *x)
1418 : {
1419 25742189 : x->d = y->d;
1420 25742189 : x->e = y->e;
1421 25742189 : }
1422 :
1423 : static void
1424 20655828 : affidpe(GEN y, dpe_t *x)
1425 : {
1426 20655828 : pari_sp av = avma;
1427 20655828 : GEN r = itor(y, DEFAULTPREC);
1428 20655155 : x->e = expo(r)+1;
1429 20655155 : setexpo(r,-1);
1430 20655073 : x->d = rtodbl(r);
1431 20655387 : set_avma(av);
1432 20655485 : }
1433 :
1434 : static void
1435 3149165 : affdbldpe(double y, dpe_t *x)
1436 : {
1437 3149165 : x->d = (double)y;
1438 3149165 : x->e = 0;
1439 3149165 : dpe_normalize(x);
1440 3149170 : }
1441 :
1442 : static void
1443 56850182 : dpe_mulz(dpe_t *x, dpe_t *y, dpe_t *z)
1444 : {
1445 56850182 : z->d = x->d * y->d;
1446 56850182 : if (z->d == 0.0)
1447 8170640 : z->e = -LONG_MAX;
1448 : else
1449 : {
1450 48679542 : z->e = x->e + y->e;
1451 48679542 : dpe_normalize0(z);
1452 : }
1453 56850369 : }
1454 :
1455 : static void
1456 14049802 : dpe_divz(dpe_t *x, dpe_t *y, dpe_t *z)
1457 : {
1458 14049802 : z->d = x->d / y->d;
1459 14049802 : if (z->d == 0.0)
1460 1918686 : z->e = -LONG_MAX;
1461 : else
1462 : {
1463 12131116 : z->e = x->e - y->e;
1464 12131116 : dpe_normalize0(z);
1465 : }
1466 14049842 : }
1467 :
1468 : static void
1469 242907 : dpe_negz(dpe_t *y, dpe_t *x)
1470 : {
1471 242907 : x->d = - y->d;
1472 242907 : x->e = y->e;
1473 242907 : }
1474 :
1475 : static void
1476 1952794 : dpe_addz(dpe_t *y, dpe_t *z, dpe_t *x)
1477 : {
1478 1952794 : if (y->e > z->e + 53)
1479 113340 : affdpe(y, x);
1480 1839454 : else if (z->e > y->e + 53)
1481 41842 : affdpe(z, x);
1482 : else
1483 : {
1484 1797612 : long d = y->e - z->e;
1485 :
1486 1797612 : if (d >= 0)
1487 : {
1488 1349030 : x->d = y->d + ldexp(z->d, -d);
1489 1349030 : x->e = y->e;
1490 : }
1491 : else
1492 : {
1493 448582 : x->d = z->d + ldexp(y->d, d);
1494 448582 : x->e = z->e;
1495 : }
1496 1797612 : dpe_normalize(x);
1497 : }
1498 1952793 : }
1499 : static void
1500 53644161 : dpe_subz(dpe_t *y, dpe_t *z, dpe_t *x)
1501 : {
1502 53644161 : if (y->e > z->e + 53)
1503 11204384 : affdpe(y, x);
1504 42439777 : else if (z->e > y->e + 53)
1505 242907 : dpe_negz(z, x);
1506 : else
1507 : {
1508 42196870 : long d = y->e - z->e;
1509 :
1510 42196870 : if (d >= 0)
1511 : {
1512 39423539 : x->d = y->d - ldexp(z->d, -d);
1513 39423539 : x->e = y->e;
1514 : }
1515 : else
1516 : {
1517 2773331 : x->d = ldexp(y->d, d) - z->d;
1518 2773331 : x->e = z->e;
1519 : }
1520 42196870 : dpe_normalize(x);
1521 : }
1522 53644390 : }
1523 :
1524 : static void
1525 799459 : dpe_muluz(dpe_t *y, ulong t, dpe_t *x)
1526 : {
1527 799459 : x->d = y->d * (double)t;
1528 799459 : x->e = y->e;
1529 799459 : dpe_normalize(x);
1530 799459 : }
1531 :
1532 : static void
1533 343410 : dpe_addmuluz(dpe_t *y, dpe_t *z, ulong t, dpe_t *x)
1534 : {
1535 : dpe_t tmp;
1536 343410 : dpe_muluz(z, t, &tmp);
1537 343410 : dpe_addz(y, &tmp, x);
1538 343410 : }
1539 :
1540 : static void
1541 408850 : dpe_submuluz(dpe_t *y, dpe_t *z, ulong t, dpe_t *x)
1542 : {
1543 : dpe_t tmp;
1544 408850 : dpe_muluz(z, t, &tmp);
1545 408850 : dpe_subz(y, &tmp, x);
1546 408850 : }
1547 :
1548 : static void
1549 51607565 : dpe_submulz(dpe_t *y, dpe_t *z, dpe_t *t, dpe_t *x)
1550 : {
1551 : dpe_t tmp;
1552 51607565 : dpe_mulz(z, t, &tmp);
1553 51607572 : dpe_subz(y, &tmp, x);
1554 51607605 : }
1555 :
1556 : static int
1557 5242720 : dpe_cmp(dpe_t *x, dpe_t *y)
1558 : {
1559 5242720 : int sx = x->d < 0. ? -1: x->d > 0.;
1560 5242720 : int sy = y->d < 0. ? -1: y->d > 0.;
1561 5242720 : int d = sx - sy;
1562 :
1563 5242720 : if (d != 0)
1564 141657 : return d;
1565 5101063 : else if (x->e > y->e)
1566 489061 : return (sx > 0) ? 1 : -1;
1567 4612002 : else if (y->e > x->e)
1568 2521880 : return (sx > 0) ? -1 : 1;
1569 : else
1570 2090122 : return (x->d < y->d) ? -1 : (x->d > y->d);
1571 : }
1572 :
1573 : static int
1574 14486056 : dpe_abscmp(dpe_t *x, dpe_t *y)
1575 : {
1576 14486056 : if (x->e > y->e)
1577 277378 : return 1;
1578 14208678 : else if (y->e > x->e)
1579 13354252 : return -1;
1580 : else
1581 854426 : return (fabs(x->d) < fabs(y->d)) ? -1 : (fabs(x->d) > fabs(y->d));
1582 : }
1583 :
1584 : static int
1585 1406695 : dpe_abssmall(dpe_t *x)
1586 : {
1587 1406695 : return (x->e <= 0) || (x->e == 1 && fabs(x->d) <= .75);
1588 : }
1589 :
1590 : static int
1591 5242722 : dpe_cmpmul(dpe_t *x, dpe_t *y, dpe_t *z)
1592 : {
1593 : dpe_t t;
1594 5242722 : dpe_mulz(x,y,&t);
1595 5242719 : return dpe_cmp(&t, z);
1596 : }
1597 :
1598 : static dpe_t *
1599 13100007 : cget_dpevec(long d)
1600 13100007 : { return (dpe_t*) stack_malloc_align(d*sizeof(dpe_t), sizeof(dpe_t)); }
1601 :
1602 : static dpe_t **
1603 3149154 : cget_dpemat(long d) { return (dpe_t **) cgetg(d, t_VECSMALL); }
1604 :
1605 : static GEN
1606 20028 : dpeM_diagonal_shallow(dpe_t **m, long d)
1607 : {
1608 : long i;
1609 20028 : GEN y = cgetg(d+1,t_VEC);
1610 113437 : for (i=1; i<=d; i++) gel(y, i) = dpetor(Dmael(m,i,i));
1611 20027 : return y;
1612 : }
1613 :
1614 : static void
1615 1406672 : affii_or_copy_gc(pari_sp av, GEN x, GEN *y)
1616 : {
1617 1406672 : long l = lg(*y);
1618 1406672 : if (lgefint(x) <= l && isonstack(*y))
1619 : {
1620 1406242 : affii(x,*y);
1621 1406249 : set_avma(av);
1622 : }
1623 : else
1624 432 : *y = gerepileuptoint(av, x);
1625 1406681 : }
1626 :
1627 : /* *x -= u*y */
1628 : INLINE void
1629 5908095 : submulziu(GEN *x, GEN y, ulong u)
1630 : {
1631 : pari_sp av;
1632 5908095 : long ly = lgefint(y);
1633 5908095 : if (ly == 2) return;
1634 3253814 : av = avma;
1635 3253814 : (void)new_chunk(3+ly+lgefint(*x)); /* HACK */
1636 3253871 : y = mului(u,y);
1637 3253863 : set_avma(av); subzi(x, y);
1638 : }
1639 :
1640 : /* *x += u*y */
1641 : INLINE void
1642 4600206 : addmulziu(GEN *x, GEN y, ulong u)
1643 : {
1644 : pari_sp av;
1645 4600206 : long ly = lgefint(y);
1646 4600206 : if (ly == 2) return;
1647 2777135 : av = avma;
1648 2777135 : (void)new_chunk(3+ly+lgefint(*x)); /* HACK */
1649 2777144 : y = mului(u,y);
1650 2777140 : set_avma(av); addzi(x, y);
1651 : }
1652 :
1653 : /************************** PROVED version (dpe) *************************/
1654 :
1655 : /* Babai's Nearest Plane algorithm (iterative).
1656 : * Size-reduces b_kappa using mu_{i,j} and r_{i,j} for j<=i <kappa
1657 : * Update B[,kappa]; compute mu_{kappa,j}, r_{kappa,j} for j<=kappa and s[kappa]
1658 : * mu, r, s updated in place (affrr). Return 1 on failure, else 0. */
1659 : static int
1660 4621562 : Babai_dpe(pari_sp av, long kappa, GEN *pG, GEN *pB, GEN *pU, dpe_t **mu, dpe_t **r, dpe_t *s,
1661 : long a, long zeros, long maxG, dpe_t *eta)
1662 : {
1663 4621562 : GEN G = *pG, B = *pB, U = *pU, ztmp;
1664 4621562 : long k, d, n, aa = a > zeros? a: zeros+1;
1665 4621562 : long emaxmu = EX0, emax2mu = EX0;
1666 : /* N.B: we set d = 0 (resp. n = 0) to avoid updating U (resp. B) */
1667 4621562 : d = U? lg(U)-1: 0;
1668 4621562 : n = B? nbrows(B): 0;
1669 532989 : for (;;) {
1670 5154585 : int go_on = 0;
1671 5154585 : long i, j, emax3mu = emax2mu;
1672 :
1673 5154585 : if (gc_needed(av,2))
1674 : {
1675 0 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[1], a=%ld", aa);
1676 0 : gc_lll(av,3,&G,&B,&U);
1677 : }
1678 : /* Step2: compute the GSO for stage kappa */
1679 5154574 : emax2mu = emaxmu; emaxmu = EX0;
1680 19204401 : for (j=aa; j<kappa; j++)
1681 : {
1682 : dpe_t g;
1683 14049836 : affidpe(gmael(G,kappa,j), &g);
1684 52455853 : for (k = zeros+1; k < j; k++)
1685 38406138 : dpe_submulz(&g, Dmael(mu,j,k), Dmael(r,kappa,k), &g);
1686 14049715 : affdpe(&g, Dmael(r,kappa,j));
1687 14049804 : dpe_divz(Dmael(r,kappa,j), Dmael(r,j,j), Dmael(mu,kappa,j));
1688 14049798 : emaxmu = maxss(emaxmu, Dmael(mu,kappa,j)->e);
1689 : }
1690 5154565 : if (emax3mu != EX0 && emax3mu <= emax2mu + 5) /* precision too low */
1691 0 : { *pG = G; *pB = B; *pU = U; return 1; }
1692 :
1693 19107634 : for (j=kappa-1; j>zeros; j--)
1694 14486058 : if (dpe_abscmp(Dmael(mu,kappa,j), eta) > 0) { go_on = 1; break; }
1695 :
1696 : /* Step3--5: compute the X_j's */
1697 5154561 : if (go_on)
1698 3079812 : for (j=kappa-1; j>zeros; j--)
1699 : {
1700 : pari_sp btop;
1701 2546824 : dpe_t *tmp = Dmael(mu,kappa,j);
1702 2546824 : if (tmp->e < 0) continue; /* (essentially) size-reduced */
1703 :
1704 1406696 : if (gc_needed(av,2))
1705 : {
1706 0 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[2], a=%ld, j=%ld", aa,j);
1707 0 : gc_lll(av,3,&G,&B,&U);
1708 : }
1709 : /* we consider separately the case |X| = 1 */
1710 1406696 : if (dpe_abssmall(tmp))
1711 : {
1712 931560 : if (tmp->d > 0) { /* in this case, X = 1 */
1713 2071054 : for (k=zeros+1; k<j; k++)
1714 1604357 : dpe_subz(Dmael(mu,kappa,k), Dmael(mu,j,k), Dmael(mu,kappa,k));
1715 3042279 : for (i=1; i<=n; i++)
1716 2575582 : subzi(&gmael(B,kappa,i), gmael(B,j,i));
1717 7002249 : for (i=1; i<=d; i++)
1718 6535558 : subzi(&gmael(U,kappa,i), gmael(U,j,i));
1719 466691 : btop = avma;
1720 466691 : ztmp = subii(gmael(G,j,j), shifti(gmael(G,kappa,j), 1));
1721 466697 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
1722 466695 : affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
1723 2878158 : for (i=1; i<=j; i++)
1724 2411462 : subzi(&gmael(G,kappa,i), gmael(G,j,i));
1725 2202727 : for (i=j+1; i<kappa; i++)
1726 1736032 : subzi(&gmael(G,kappa,i), gmael(G,i,j));
1727 2370051 : for (i=kappa+1; i<=maxG; i++)
1728 1903355 : subzi(&gmael(G,i,kappa), gmael(G,i,j));
1729 : } else { /* otherwise X = -1 */
1730 2050470 : for (k=zeros+1; k<j; k++)
1731 1585607 : dpe_addz(Dmael(mu,kappa,k), Dmael(mu,j,k), Dmael(mu,kappa,k));
1732 3039267 : for (i=1; i<=n; i++)
1733 2574404 : addzi(&gmael(B,kappa,i),gmael(B,j,i));
1734 6881449 : for (i=1; i<=d; i++)
1735 6416594 : addzi(&gmael(U,kappa,i),gmael(U,j,i));
1736 464855 : btop = avma;
1737 464855 : ztmp = addii(gmael(G,j,j), shifti(gmael(G,kappa,j), 1));
1738 464860 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
1739 464858 : affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
1740 2790967 : for (i=1; i<=j; i++)
1741 2326104 : addzi(&gmael(G,kappa,i), gmael(G,j,i));
1742 2209287 : for (i=j+1; i<kappa; i++)
1743 1744425 : addzi(&gmael(G,kappa,i), gmael(G,i,j));
1744 2324538 : for (i=kappa+1; i<=maxG; i++)
1745 1859675 : addzi(&gmael(G,i,kappa), gmael(G,i,j));
1746 : }
1747 931559 : continue;
1748 : }
1749 : /* we have |X| >= 2 */
1750 475135 : if (tmp->e < BITS_IN_LONG-1)
1751 : {
1752 450973 : if (tmp->d > 0)
1753 : {
1754 247632 : ulong xx = (ulong) pari_rint(ldexp(tmp->d, tmp->e)); /* X fits in an ulong */
1755 656482 : for (k=zeros+1; k<j; k++)
1756 408850 : dpe_submuluz(Dmael(mu,kappa,k), Dmael(mu,j,k), xx, Dmael(mu,kappa,k));
1757 734503 : for (i=1; i<=n; i++)
1758 486871 : submulziu(&gmael(B,kappa,i), gmael(B,j,i), xx);
1759 3096203 : for (i=1; i<=d; i++)
1760 2848576 : submulziu(&gmael(U,kappa,i), gmael(U,j,i), xx);
1761 247627 : btop = avma;
1762 247627 : ztmp = submuliu2n(mulii(gmael(G,j,j), sqru(xx)), gmael(G,kappa,j), xx, 1);
1763 247632 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
1764 247632 : affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
1765 1311672 : for (i=1; i<=j; i++)
1766 1064041 : submulziu(&gmael(G,kappa,i), gmael(G,j,i), xx);
1767 804807 : for (i=j+1; i<kappa; i++)
1768 557177 : submulziu(&gmael(G,kappa,i), gmael(G,i,j), xx);
1769 1199146 : for (i=kappa+1; i<=maxG; i++)
1770 951516 : submulziu(&gmael(G,i,kappa), gmael(G,i,j), xx);
1771 : }
1772 : else
1773 : {
1774 203341 : ulong xx = (ulong) pari_rint(ldexp(-tmp->d, tmp->e)); /* X fits in an ulong */
1775 546751 : for (k=zeros+1; k<j; k++)
1776 343410 : dpe_addmuluz(Dmael(mu,kappa,k), Dmael(mu,j,k), xx, Dmael(mu,kappa,k));
1777 703632 : for (i=1; i<=n; i++)
1778 500291 : addmulziu(&gmael(B,kappa,i), gmael(B,j,i), xx);
1779 2367790 : for (i=1; i<=d; i++)
1780 2164450 : addmulziu(&gmael(U,kappa,i), gmael(U,j,i), xx);
1781 203340 : btop = avma;
1782 203340 : ztmp = addmuliu2n(mulii(gmael(G,j,j), sqru(xx)), gmael(G,kappa,j), xx, 1);
1783 203340 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
1784 203340 : affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
1785 995595 : for (i=1; i<=j; i++)
1786 792254 : addmulziu(&gmael(G,kappa,i), gmael(G,j,i), xx);
1787 664672 : for (i=j+1; i<kappa; i++)
1788 461332 : addmulziu(&gmael(G,kappa,i), gmael(G,i,j), xx);
1789 885242 : for (i=kappa+1; i<=maxG; i++)
1790 681902 : addmulziu(&gmael(G,i,kappa), gmael(G,i,j), xx);
1791 : }
1792 : }
1793 : else
1794 : {
1795 24162 : long e = tmp->e - BITS_IN_LONG + 1;
1796 24162 : if (tmp->d > 0)
1797 : {
1798 12035 : ulong xx = (ulong) pari_rint(ldexp(tmp->d, BITS_IN_LONG - 1));
1799 35455 : for (k=zeros+1; k<j; k++)
1800 : {
1801 : dpe_t x;
1802 23420 : dpe_muluz(Dmael(mu,j,k), xx, &x);
1803 23420 : x.e += e;
1804 23420 : dpe_subz(Dmael(mu,kappa,k), &x, Dmael(mu,kappa,k));
1805 : }
1806 148590 : for (i=1; i<=n; i++)
1807 136555 : submulzu2n(&gmael(B,kappa,i), gmael(B,j,i), xx, e);
1808 111267 : for (i=1; i<=d; i++)
1809 99232 : submulzu2n(&gmael(U,kappa,i), gmael(U,j,i), xx, e);
1810 12035 : btop = avma;
1811 12035 : ztmp = submuliu2n(mulshift(gmael(G,j,j), sqru(xx), 2*e),
1812 12035 : gmael(G,kappa,j), xx, e+1);
1813 12035 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
1814 12035 : affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
1815 48096 : for (i=1; i<=j; i++)
1816 36061 : submulzu2n(&gmael(G,kappa,i), gmael(G,j,i), xx, e);
1817 58095 : for ( ; i<kappa; i++)
1818 46060 : submulzu2n(&gmael(G,kappa,i), gmael(G,i,j), xx, e);
1819 14149 : for (i=kappa+1; i<=maxG; i++)
1820 2114 : submulzu2n(&gmael(G,i,kappa), gmael(G,i,j), xx, e);
1821 : } else
1822 : {
1823 12127 : ulong xx = (ulong) pari_rint(ldexp(-tmp->d, BITS_IN_LONG - 1));
1824 35914 : for (k=zeros+1; k<j; k++)
1825 : {
1826 : dpe_t x;
1827 23779 : dpe_muluz(Dmael(mu,j,k), xx, &x);
1828 23779 : x.e += e;
1829 23779 : dpe_addz(Dmael(mu,kappa,k), &x, Dmael(mu,kappa,k));
1830 : }
1831 152427 : for (i=1; i<=n; i++)
1832 140292 : addmulzu2n(&gmael(B,kappa,i), gmael(B,j,i), xx, e);
1833 113043 : for (i=1; i<=d; i++)
1834 100908 : addmulzu2n(&gmael(U,kappa,i), gmael(U,j,i), xx, e);
1835 12135 : btop = avma;
1836 12135 : ztmp = addmuliu2n(mulshift(gmael(G,j,j), sqru(xx), 2*e),
1837 12135 : gmael(G,kappa,j), xx, e+1);
1838 12135 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
1839 12135 : affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
1840 48823 : for (i=1; i<=j; i++)
1841 36688 : addmulzu2n(&gmael(G,kappa,i), gmael(G,j,i), xx, e);
1842 59402 : for ( ; i<kappa; i++)
1843 47267 : addmulzu2n(&gmael(G,kappa,i), gmael(G,i,j), xx, e);
1844 13920 : for (i=kappa+1; i<=maxG; i++)
1845 1785 : addmulzu2n(&gmael(G,i,kappa), gmael(G,i,j), xx, e);
1846 : }
1847 : }
1848 : }
1849 5154564 : if (!go_on) break; /* Anything happened? */
1850 532989 : aa = zeros+1;
1851 : }
1852 :
1853 4621575 : affidpe(gmael(G,kappa,kappa), Del(s,zeros+1));
1854 : /* the last s[kappa-1]=r[kappa][kappa] is computed only if kappa increases */
1855 13534620 : for (k=zeros+1; k<=kappa-2; k++)
1856 8913051 : dpe_submulz(Del(s,k), Dmael(mu,kappa,k), Dmael(r,kappa,k), Del(s,k+1));
1857 4621569 : *pG = G; *pB = B; *pU = U; return 0;
1858 : }
1859 :
1860 : /* G integral Gram matrix, LLL-reduces (G,B,U) in place [apply base change
1861 : * transforms to B and U]. If (keepfirst), never swap with first vector.
1862 : * If G = NULL, we compute the Gram matrix incrementally.
1863 : * Return -1 on failure, else zeros = dim Kernel (>= 0) */
1864 : static long
1865 1574584 : fplll_dpe(GEN *pG, GEN *pB, GEN *pU, GEN *pr, double DELTA, double ETA,
1866 : long keepfirst)
1867 : {
1868 : pari_sp av;
1869 1574584 : GEN Gtmp, alpha, G = *pG, B = *pB, U = *pU;
1870 1574584 : long d, maxG, kappa, kappa2, i, j, zeros, kappamax, incgram = !G, cnt = 0;
1871 : dpe_t delta, eta, **mu, **r, *s;
1872 1574584 : affdbldpe(DELTA,&delta);
1873 1574587 : affdbldpe(ETA,&eta);
1874 :
1875 1574588 : if (incgram)
1876 : { /* incremental Gram matrix */
1877 1514126 : maxG = 2; d = lg(B)-1;
1878 1514126 : G = zeromatcopy(d, d);
1879 : }
1880 : else
1881 60462 : maxG = d = lg(G)-1;
1882 :
1883 1574585 : mu = cget_dpemat(d+1);
1884 1574576 : r = cget_dpemat(d+1);
1885 1574580 : s = cget_dpevec(d+1);
1886 7337326 : for (j = 1; j <= d; j++)
1887 : {
1888 5762747 : mu[j]= cget_dpevec(d+1);
1889 5762741 : r[j] = cget_dpevec(d+1);
1890 : }
1891 1574579 : Gtmp = cgetg(d+1, t_VEC);
1892 1574576 : alpha = cgetg(d+1, t_VECSMALL);
1893 1574576 : av = avma;
1894 :
1895 : /* Step2: Initializing the main loop */
1896 1574576 : kappamax = 1;
1897 1574576 : i = 1;
1898 : do {
1899 1949513 : if (incgram) gmael(G,i,i) = ZV_dotsquare(gel(B,i));
1900 1949469 : affidpe(gmael(G,i,i), Dmael(r,i,i));
1901 1949502 : } while (!signe(gmael(G,i,i)) && ++i <= d);
1902 1574565 : zeros = i-1; /* all basis vectors b_i with i <= zeros are zero vectors */
1903 1574565 : kappa = i;
1904 6962322 : for (i=zeros+1; i<=d; i++) alpha[i]=1;
1905 :
1906 6196159 : while (++kappa <= d)
1907 : {
1908 4621584 : if (kappa > kappamax)
1909 : {
1910 3813177 : if (DEBUGLEVEL>=4) err_printf("K%ld ",kappa);
1911 3813177 : kappamax = kappa;
1912 3813177 : if (incgram)
1913 : {
1914 15966531 : for (i=zeros+1; i<=kappa; i++)
1915 12355055 : gmael(G,kappa,i) = ZV_dotproduct(gel(B,kappa), gel(B,i));
1916 3611476 : maxG = kappamax;
1917 : }
1918 : }
1919 : /* Step3: Call to the Babai algorithm, mu,r,s updated in place */
1920 4621392 : if (Babai_dpe(av, kappa, &G,&B,&U, mu,r,s, alpha[kappa], zeros, maxG, &eta))
1921 0 : { *pG = incgram? NULL: G; *pB = B; *pU = U; return -1; }
1922 9145174 : if ((keepfirst && kappa == 2) ||
1923 4523577 : dpe_cmpmul(Dmael(r,kappa-1,kappa-1), &delta, Del(s,kappa-1)) <= 0)
1924 : { /* Step4: Success of Lovasz's condition */
1925 4288664 : alpha[kappa] = kappa;
1926 4288664 : dpe_submulz(Del(s,kappa-1), Dmael(mu,kappa,kappa-1), Dmael(r,kappa,kappa-1), Dmael(r,kappa,kappa));
1927 4288663 : continue;
1928 : }
1929 : /* Step5: Find the right insertion index kappa, kappa2 = initial kappa */
1930 332933 : if (DEBUGLEVEL>=4 && kappa==kappamax && Del(s,kappa-1)->d)
1931 0 : if (++cnt > 20) { cnt = 0; err_printf("(%ld) ", Del(s,1)->e-1); }
1932 332933 : kappa2 = kappa;
1933 : do {
1934 843568 : kappa--;
1935 843568 : if (kappa < zeros+2 + (keepfirst ? 1: 0)) break;
1936 719145 : } while (dpe_cmpmul(Dmael(r,kappa-1,kappa-1), &delta, Del(s,kappa-1)) >= 0);
1937 332933 : update_alpha(alpha, kappa, kappa2, kappamax);
1938 :
1939 : /* Step6: Update the mu's and r's */
1940 332933 : dperotate(mu, kappa2, kappa);
1941 332933 : dperotate(r, kappa2, kappa);
1942 332933 : affdpe(Del(s,kappa), Dmael(r,kappa,kappa));
1943 :
1944 : /* Step7: Update G, B, U */
1945 332933 : if (U) rotate(U, kappa2, kappa);
1946 332933 : if (B) rotate(B, kappa2, kappa);
1947 332933 : rotateG(G,kappa2,kappa, maxG, Gtmp);
1948 :
1949 : /* Step8: Prepare the next loop iteration */
1950 332933 : if (kappa == zeros+1 && !signe(gmael(G,kappa,kappa)))
1951 : {
1952 35161 : zeros++; kappa++;
1953 35161 : affidpe(gmael(G,kappa,kappa), Dmael(r,kappa,kappa));
1954 : }
1955 : }
1956 1574575 : if (pr) *pr = dpeM_diagonal_shallow(r,d);
1957 1574575 : *pG = G; *pB = B; *pU = U; return zeros; /* success */
1958 : }
1959 :
1960 :
1961 : /************************** PROVED version (t_INT) *************************/
1962 :
1963 : /* Babai's Nearest Plane algorithm (iterative).
1964 : * Size-reduces b_kappa using mu_{i,j} and r_{i,j} for j<=i <kappa
1965 : * Update B[,kappa]; compute mu_{kappa,j}, r_{kappa,j} for j<=kappa and s[kappa]
1966 : * mu, r, s updated in place (affrr). Return 1 on failure, else 0. */
1967 : static int
1968 0 : Babai(pari_sp av, long kappa, GEN *pG, GEN *pB, GEN *pU, GEN mu, GEN r, GEN s,
1969 : long a, long zeros, long maxG, GEN eta, long prec)
1970 : {
1971 0 : GEN G = *pG, B = *pB, U = *pU, ztmp;
1972 0 : long k, aa = a > zeros? a: zeros+1;
1973 0 : const long n = B? nbrows(B): 0, d = U ? lg(U)-1: 0, bit = prec2nbits(prec);
1974 0 : long emaxmu = EX0, emax2mu = EX0;
1975 : /* N.B: we set d = 0 (resp. n = 0) to avoid updating U (resp. B) */
1976 :
1977 0 : for (;;) {
1978 0 : int go_on = 0;
1979 0 : long i, j, emax3mu = emax2mu;
1980 :
1981 0 : if (gc_needed(av,2))
1982 : {
1983 0 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[1], a=%ld", aa);
1984 0 : gc_lll(av,3,&G,&B,&U);
1985 : }
1986 : /* Step2: compute the GSO for stage kappa */
1987 0 : emax2mu = emaxmu; emaxmu = EX0;
1988 0 : for (j=aa; j<kappa; j++)
1989 : {
1990 0 : pari_sp btop = avma;
1991 0 : GEN g = gmael(G,kappa,j);
1992 0 : for (k = zeros+1; k < j; k++)
1993 0 : g = mpsub(g, mulrr(gmael(mu,j,k), gmael(r,kappa,k)));
1994 0 : mpaff(g, gmael(r,kappa,j));
1995 0 : affrr(divrr(gmael(r,kappa,j), gmael(r,j,j)), gmael(mu,kappa,j));
1996 0 : emaxmu = maxss(emaxmu, expo(gmael(mu,kappa,j)));
1997 0 : set_avma(btop);
1998 : }
1999 0 : if (emax3mu != EX0 && emax3mu <= emax2mu + 5) /* precision too low */
2000 0 : { *pG = G; *pB = B; *pU = U; return 1; }
2001 :
2002 0 : for (j=kappa-1; j>zeros; j--)
2003 0 : if (abscmprr(gmael(mu,kappa,j), eta) > 0) { go_on = 1; break; }
2004 :
2005 : /* Step3--5: compute the X_j's */
2006 0 : if (go_on)
2007 0 : for (j=kappa-1; j>zeros; j--)
2008 : {
2009 : pari_sp btop;
2010 0 : GEN tmp = gmael(mu,kappa,j);
2011 0 : if (absrsmall(tmp)) continue; /* size-reduced */
2012 :
2013 0 : if (gc_needed(av,2))
2014 : {
2015 0 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[2], a=%ld, j=%ld", aa,j);
2016 0 : gc_lll(av,3,&G,&B,&U);
2017 : }
2018 0 : btop = avma;
2019 : /* we consider separately the case |X| = 1 */
2020 0 : if (absrsmall2(tmp))
2021 : {
2022 0 : if (signe(tmp) > 0) { /* in this case, X = 1 */
2023 0 : for (k=zeros+1; k<j; k++)
2024 0 : affrr(subrr(gmael(mu,kappa,k), gmael(mu,j,k)), gmael(mu,kappa,k));
2025 0 : set_avma(btop);
2026 0 : for (i=1; i<=n; i++)
2027 0 : gmael(B,kappa,i) = subii(gmael(B,kappa,i), gmael(B,j,i));
2028 0 : for (i=1; i<=d; i++)
2029 0 : gmael(U,kappa,i) = subii(gmael(U,kappa,i), gmael(U,j,i));
2030 0 : btop = avma;
2031 0 : ztmp = subii(gmael(G,j,j), shifti(gmael(G,kappa,j), 1));
2032 0 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
2033 0 : gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
2034 0 : for (i=1; i<=j; i++)
2035 0 : gmael(G,kappa,i) = subii(gmael(G,kappa,i), gmael(G,j,i));
2036 0 : for (i=j+1; i<kappa; i++)
2037 0 : gmael(G,kappa,i) = subii(gmael(G,kappa,i), gmael(G,i,j));
2038 0 : for (i=kappa+1; i<=maxG; i++)
2039 0 : gmael(G,i,kappa) = subii(gmael(G,i,kappa), gmael(G,i,j));
2040 : } else { /* otherwise X = -1 */
2041 0 : for (k=zeros+1; k<j; k++)
2042 0 : affrr(addrr(gmael(mu,kappa,k), gmael(mu,j,k)), gmael(mu,kappa,k));
2043 0 : set_avma(btop);
2044 0 : for (i=1; i<=n; i++)
2045 0 : gmael(B,kappa,i) = addii(gmael(B,kappa,i),gmael(B,j,i));
2046 0 : for (i=1; i<=d; i++)
2047 0 : gmael(U,kappa,i) = addii(gmael(U,kappa,i),gmael(U,j,i));
2048 0 : btop = avma;
2049 0 : ztmp = addii(gmael(G,j,j), shifti(gmael(G,kappa,j), 1));
2050 0 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
2051 0 : gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
2052 0 : for (i=1; i<=j; i++)
2053 0 : gmael(G,kappa,i) = addii(gmael(G,kappa,i), gmael(G,j,i));
2054 0 : for (i=j+1; i<kappa; i++)
2055 0 : gmael(G,kappa,i) = addii(gmael(G,kappa,i), gmael(G,i,j));
2056 0 : for (i=kappa+1; i<=maxG; i++)
2057 0 : gmael(G,i,kappa) = addii(gmael(G,i,kappa), gmael(G,i,j));
2058 : }
2059 0 : continue;
2060 : }
2061 : /* we have |X| >= 2 */
2062 0 : if (expo(tmp) < BITS_IN_LONG)
2063 : {
2064 0 : ulong xx = roundr_safe(tmp)[2]; /* X fits in an ulong */
2065 0 : if (signe(tmp) > 0) /* = xx */
2066 : {
2067 0 : for (k=zeros+1; k<j; k++)
2068 0 : affrr(subrr(gmael(mu,kappa,k), mulur(xx, gmael(mu,j,k))),
2069 0 : gmael(mu,kappa,k));
2070 0 : set_avma(btop);
2071 0 : for (i=1; i<=n; i++)
2072 0 : gmael(B,kappa,i) = submuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
2073 0 : for (i=1; i<=d; i++)
2074 0 : gmael(U,kappa,i) = submuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
2075 0 : btop = avma;
2076 0 : ztmp = submuliu2n(mulii(gmael(G,j,j), sqru(xx)), gmael(G,kappa,j), xx, 1);
2077 0 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
2078 0 : gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
2079 0 : for (i=1; i<=j; i++)
2080 0 : gmael(G,kappa,i) = submuliu_inplace(gmael(G,kappa,i), gmael(G,j,i), xx);
2081 0 : for (i=j+1; i<kappa; i++)
2082 0 : gmael(G,kappa,i) = submuliu_inplace(gmael(G,kappa,i), gmael(G,i,j), xx);
2083 0 : for (i=kappa+1; i<=maxG; i++)
2084 0 : gmael(G,i,kappa) = submuliu_inplace(gmael(G,i,kappa), gmael(G,i,j), xx);
2085 : }
2086 : else /* = -xx */
2087 : {
2088 0 : for (k=zeros+1; k<j; k++)
2089 0 : affrr(addrr(gmael(mu,kappa,k), mulur(xx, gmael(mu,j,k))),
2090 0 : gmael(mu,kappa,k));
2091 0 : set_avma(btop);
2092 0 : for (i=1; i<=n; i++)
2093 0 : gmael(B,kappa,i) = addmuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
2094 0 : for (i=1; i<=d; i++)
2095 0 : gmael(U,kappa,i) = addmuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
2096 0 : btop = avma;
2097 0 : ztmp = addmuliu2n(mulii(gmael(G,j,j), sqru(xx)), gmael(G,kappa,j), xx, 1);
2098 0 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
2099 0 : gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
2100 0 : for (i=1; i<=j; i++)
2101 0 : gmael(G,kappa,i) = addmuliu_inplace(gmael(G,kappa,i), gmael(G,j,i), xx);
2102 0 : for (i=j+1; i<kappa; i++)
2103 0 : gmael(G,kappa,i) = addmuliu_inplace(gmael(G,kappa,i), gmael(G,i,j), xx);
2104 0 : for (i=kappa+1; i<=maxG; i++)
2105 0 : gmael(G,i,kappa) = addmuliu_inplace(gmael(G,i,kappa), gmael(G,i,j), xx);
2106 : }
2107 : }
2108 : else
2109 : {
2110 : long e;
2111 0 : GEN X = truncexpo(tmp, bit, &e); /* tmp ~ X * 2^e */
2112 0 : btop = avma;
2113 0 : for (k=zeros+1; k<j; k++)
2114 : {
2115 0 : GEN x = mulir(X, gmael(mu,j,k));
2116 0 : if (e) shiftr_inplace(x, e);
2117 0 : affrr(subrr(gmael(mu,kappa,k), x), gmael(mu,kappa,k));
2118 : }
2119 0 : set_avma(btop);
2120 0 : for (i=1; i<=n; i++)
2121 0 : gmael(B,kappa,i) = submulshift(gmael(B,kappa,i), gmael(B,j,i), X, e);
2122 0 : for (i=1; i<=d; i++)
2123 0 : gmael(U,kappa,i) = submulshift(gmael(U,kappa,i), gmael(U,j,i), X, e);
2124 0 : btop = avma;
2125 0 : ztmp = submulshift(mulshift(gmael(G,j,j), sqri(X), 2*e),
2126 0 : gmael(G,kappa,j), X, e+1);
2127 0 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
2128 0 : gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
2129 0 : for (i=1; i<=j; i++)
2130 0 : gmael(G,kappa,i) = submulshift(gmael(G,kappa,i), gmael(G,j,i), X, e);
2131 0 : for ( ; i<kappa; i++)
2132 0 : gmael(G,kappa,i) = submulshift(gmael(G,kappa,i), gmael(G,i,j), X, e);
2133 0 : for (i=kappa+1; i<=maxG; i++)
2134 0 : gmael(G,i,kappa) = submulshift(gmael(G,i,kappa), gmael(G,i,j), X, e);
2135 : }
2136 : }
2137 0 : if (!go_on) break; /* Anything happened? */
2138 0 : aa = zeros+1;
2139 : }
2140 :
2141 0 : affir(gmael(G,kappa,kappa), gel(s,zeros+1));
2142 : /* the last s[kappa-1]=r[kappa][kappa] is computed only if kappa increases */
2143 0 : av = avma;
2144 0 : for (k=zeros+1; k<=kappa-2; k++)
2145 0 : affrr(subrr(gel(s,k), mulrr(gmael(mu,kappa,k), gmael(r,kappa,k))),
2146 0 : gel(s,k+1));
2147 0 : *pG = G; *pB = B; *pU = U; return gc_bool(av, 0);
2148 : }
2149 :
2150 : /* G integral Gram matrix, LLL-reduces (G,B,U) in place [apply base change
2151 : * transforms to B and U]. If (keepfirst), never swap with first vector.
2152 : * If G = NULL, we compute the Gram matrix incrementally.
2153 : * Return -1 on failure, else zeros = dim Kernel (>= 0) */
2154 : static long
2155 0 : fplll(GEN *pG, GEN *pB, GEN *pU, GEN *pr, double DELTA, double ETA,
2156 : long keepfirst, long prec)
2157 : {
2158 : pari_sp av, av2;
2159 0 : GEN mu, r, s, tmp, Gtmp, alpha, G = *pG, B = *pB, U = *pU;
2160 0 : GEN delta = dbltor(DELTA), eta = dbltor(ETA);
2161 0 : long d, maxG, kappa, kappa2, i, j, zeros, kappamax, incgram = !G, cnt = 0;
2162 :
2163 0 : if (incgram)
2164 : { /* incremental Gram matrix */
2165 0 : maxG = 2; d = lg(B)-1;
2166 0 : G = zeromatcopy(d, d);
2167 : }
2168 : else
2169 0 : maxG = d = lg(G)-1;
2170 :
2171 0 : mu = cgetg(d+1, t_MAT);
2172 0 : r = cgetg(d+1, t_MAT);
2173 0 : s = cgetg(d+1, t_VEC);
2174 0 : for (j = 1; j <= d; j++)
2175 : {
2176 0 : GEN M = cgetg(d+1, t_COL), R = cgetg(d+1, t_COL);
2177 0 : gel(mu,j)= M;
2178 0 : gel(r,j) = R;
2179 0 : gel(s,j) = cgetr(prec);
2180 0 : for (i = 1; i <= d; i++)
2181 : {
2182 0 : gel(R,i) = cgetr(prec);
2183 0 : gel(M,i) = cgetr(prec);
2184 : }
2185 : }
2186 0 : Gtmp = cgetg(d+1, t_VEC);
2187 0 : alpha = cgetg(d+1, t_VECSMALL);
2188 0 : av = avma;
2189 :
2190 : /* Step2: Initializing the main loop */
2191 0 : kappamax = 1;
2192 0 : i = 1;
2193 : do {
2194 0 : if (incgram) gmael(G,i,i) = ZV_dotsquare(gel(B,i));
2195 0 : affir(gmael(G,i,i), gmael(r,i,i));
2196 0 : } while (!signe(gmael(G,i,i)) && ++i <= d);
2197 0 : zeros = i-1; /* all basis vectors b_i with i <= zeros are zero vectors */
2198 0 : kappa = i;
2199 0 : for (i=zeros+1; i<=d; i++) alpha[i]=1;
2200 :
2201 0 : while (++kappa <= d)
2202 : {
2203 0 : if (kappa > kappamax)
2204 : {
2205 0 : if (DEBUGLEVEL>=4) err_printf("K%ld ",kappa);
2206 0 : kappamax = kappa;
2207 0 : if (incgram)
2208 : {
2209 0 : for (i=zeros+1; i<=kappa; i++)
2210 0 : gmael(G,kappa,i) = ZV_dotproduct(gel(B,kappa), gel(B,i));
2211 0 : maxG = kappamax;
2212 : }
2213 : }
2214 : /* Step3: Call to the Babai algorithm, mu,r,s updated in place */
2215 0 : if (Babai(av, kappa, &G,&B,&U, mu,r,s, alpha[kappa], zeros, maxG, eta, prec))
2216 0 : { *pG = incgram? NULL: G; *pB = B; *pU = U; return -1; }
2217 0 : av2 = avma;
2218 0 : if ((keepfirst && kappa == 2) ||
2219 0 : cmprr(mulrr(gmael(r,kappa-1,kappa-1), delta), gel(s,kappa-1)) <= 0)
2220 : { /* Step4: Success of Lovasz's condition */
2221 0 : alpha[kappa] = kappa;
2222 0 : tmp = mulrr(gmael(mu,kappa,kappa-1), gmael(r,kappa,kappa-1));
2223 0 : affrr(subrr(gel(s,kappa-1), tmp), gmael(r,kappa,kappa));
2224 0 : set_avma(av2); continue;
2225 : }
2226 : /* Step5: Find the right insertion index kappa, kappa2 = initial kappa */
2227 0 : if (DEBUGLEVEL>=4 && kappa==kappamax && signe(gel(s,kappa-1)))
2228 0 : if (++cnt > 20) { cnt = 0; err_printf("(%ld) ", expo(gel(s,1))); }
2229 0 : kappa2 = kappa;
2230 : do {
2231 0 : kappa--;
2232 0 : if (kappa < zeros+2 + (keepfirst ? 1: 0)) break;
2233 0 : tmp = mulrr(gmael(r,kappa-1,kappa-1), delta);
2234 0 : } while (cmprr(gel(s,kappa-1), tmp) <= 0);
2235 0 : set_avma(av2);
2236 0 : update_alpha(alpha, kappa, kappa2, kappamax);
2237 :
2238 : /* Step6: Update the mu's and r's */
2239 0 : rotate(mu, kappa2, kappa);
2240 0 : rotate(r, kappa2, kappa);
2241 0 : affrr(gel(s,kappa), gmael(r,kappa,kappa));
2242 :
2243 : /* Step7: Update G, B, U */
2244 0 : if (U) rotate(U, kappa2, kappa);
2245 0 : if (B) rotate(B, kappa2, kappa);
2246 0 : rotateG(G,kappa2,kappa, maxG, Gtmp);
2247 :
2248 : /* Step8: Prepare the next loop iteration */
2249 0 : if (kappa == zeros+1 && !signe(gmael(G,kappa,kappa)))
2250 : {
2251 0 : zeros++; kappa++;
2252 0 : affir(gmael(G,kappa,kappa), gmael(r,kappa,kappa));
2253 : }
2254 : }
2255 0 : if (pr) *pr = RgM_diagonal_shallow(r);
2256 0 : *pG = G; *pB = B; *pU = U; return zeros; /* success */
2257 : }
2258 :
2259 : /* do not support LLL_KER, LLL_ALL, LLL_KEEP_FIRST */
2260 : static GEN
2261 4683094 : ZM2_lll_norms(GEN x, long flag, GEN *pN)
2262 : {
2263 : GEN a,b,c,d;
2264 : GEN G, U;
2265 4683094 : if (flag & LLL_GRAM)
2266 15359 : G = x;
2267 : else
2268 4667735 : G = gram_matrix(x);
2269 4683081 : a = gcoeff(G,1,1); b = shifti(gcoeff(G,1,2),1); c = gcoeff(G,2,2);
2270 4683047 : d = qfb_disc3(a,b,c);
2271 4683042 : if (signe(d)>=0) return NULL;
2272 4682657 : G = redimagsl2(mkqfb(a,b,c,d),&U);
2273 4682713 : if (pN) (void) RgM_gram_schmidt(G, pN);
2274 4682712 : if (flag & LLL_INPLACE) return ZM2_mul(x,U);
2275 4682712 : return U;
2276 : }
2277 :
2278 : static void
2279 626875 : fplll_flatter(GEN *pG, GEN *pB, GEN *pU, long rank, long flag)
2280 : {
2281 626875 : if (!*pG)
2282 : {
2283 625559 : GEN T = ZM_flatter_rank(*pB, rank, flag);
2284 625560 : if (*pU)
2285 : {
2286 611269 : *pU = ZM_mul(*pU, T);
2287 611270 : *pB = ZM_mul(*pB, T);
2288 14291 : } else *pB = T;
2289 : } else
2290 : {
2291 1316 : GEN T, G = *pG;
2292 1316 : long i, j, l = lg(G);
2293 8692 : for (i = 1; i < l; i++)
2294 45781 : for(j = 1; j < i; j++)
2295 38405 : gmael(G,j,i) = gmael(G,i,j);
2296 1316 : T = ZM_flattergram_rank(G, rank, flag);
2297 1316 : if (*pU) *pU = ZM_mul(*pU, T);
2298 1316 : *pG = ZM_mul(shallowtrans(T), ZM_mul(*pG,T));
2299 : }
2300 626876 : }
2301 :
2302 : static long
2303 1127080 : spread(GEN R)
2304 : {
2305 1127080 : long i, n = lg(R)-1;
2306 1127080 : GEN m = gcoeff(R, 1, 1), M = m;
2307 4812130 : for (i = 2; i <= n; ++i)
2308 : {
2309 3685044 : if (mpabscmp(gcoeff(R, i, i), m) < 0)
2310 1290111 : m = gcoeff(R, i, i);
2311 3685048 : if (mpabscmp(gcoeff(R, i, i), M) > 0)
2312 857821 : M = gcoeff(R, i, i);
2313 : }
2314 1127086 : return mpexpo(M)-mpexpo(m);
2315 : }
2316 :
2317 : static GEN
2318 1093449 : get_gramschmidt(GEN M, long rank)
2319 : {
2320 : GEN B, Q, L;
2321 1093449 : long l = lg(M), bitprec = 3*(l-1) + 30;
2322 1093449 : long prec = nbits2prec64(bitprec);
2323 1093448 : if (rank < l-1) M = vconcat(gshift(M,1), matid(l-1));
2324 1093448 : if (!QR_init(RgM_gtofp(M, prec), &B, &Q, &L, prec) || !gsisinv(L)) return NULL;
2325 469066 : return L;
2326 : }
2327 :
2328 : static GEN
2329 49563 : get_gaussred(GEN M, long rank)
2330 : {
2331 49563 : pari_sp ltop = avma;
2332 49563 : long lM = lg(M);
2333 49563 : long bitprec = 3*(lM-1) + 30, prec = nbits2prec64(bitprec);
2334 : GEN R;
2335 49563 : if (rank < lM-1) M = RgM_Rg_add(gshift(M, 1), gen_1);
2336 49563 : R = RgM_Cholesky(RgM_gtofp(M, prec), prec);
2337 49563 : if (!R) return NULL;
2338 48247 : return gerepilecopy(ltop, R);
2339 : }
2340 :
2341 : /* Assume x a ZM, if pN != NULL, set it to Gram-Schmidt (squared) norms
2342 : * The following modes are supported:
2343 : * - flag & LLL_INPLACE: x a lattice basis, return x*U
2344 : * - flag & LLL_GRAM: x a Gram matrix / else x a lattice basis; return
2345 : * LLL base change matrix U [LLL_IM]
2346 : * kernel basis [LLL_KER, nonreduced]
2347 : * both [LLL_ALL] */
2348 : GEN
2349 6946232 : ZM_lll_norms(GEN x, double DELTA, long flag, GEN *pN)
2350 : {
2351 6946232 : pari_sp av = avma;
2352 6946232 : const double ETA = 0.51;
2353 6946232 : const long keepfirst = flag & LLL_KEEP_FIRST;
2354 6946232 : long p, zeros = -1, n = lg(x)-1, is_upper, is_lower, useflatter = 0, rank;
2355 6946232 : GEN G, B, U, L = NULL;
2356 : pari_timer T;
2357 6946232 : long thre[]={31783,34393,20894,22525,13533,1928,672,671,
2358 : 422,506,315,313,222,205,167,154,139,138,
2359 : 110,120,98,94,81,75,74,64,74,74,
2360 : 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,
2361 : 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};
2362 6946232 : long thsn[]={23280,30486,50077,44136,78724,15690,1801,1611,
2363 : 981,1359,978,1042,815,866,788,775,726,712,
2364 : 626,613,548,564,474,481,504,447,453,508,
2365 : 705,794,1008,946,767,898,886,763,842,757,
2366 : 725,774,639,655,705,627,635,704,511,613,
2367 : 583,595,568,640,541,640,567,540,577,584,
2368 : 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};
2369 6946232 : if (n <= 1) return lll_trivial(x, flag);
2370 6837522 : if (nbrows(x)==0)
2371 : {
2372 15190 : if (flag & LLL_KER) return matid(n);
2373 15190 : if (flag & (LLL_INPLACE|LLL_IM)) return cgetg(1,t_MAT);
2374 0 : retmkvec2(matid(n), cgetg(1,t_MAT));
2375 : }
2376 6822425 : if (n==2 && nbrows(x)==2 && (flag&LLL_IM) && !keepfirst)
2377 : {
2378 4683094 : U = ZM2_lll_norms(x, flag, pN);
2379 4683096 : if (U) return U;
2380 : }
2381 2139711 : rank = ZM_rank(x);
2382 2139707 : if (flag & LLL_GRAM)
2383 60461 : { G = x; B = NULL; U = matid(n); is_upper = 0; is_lower = 0; }
2384 : else
2385 : {
2386 2079246 : G = NULL; B = x; U = (flag & LLL_INPLACE)? NULL: matid(n);
2387 2079251 : is_upper = (flag & LLL_UPPER) || ZM_is_upper(B);
2388 2079250 : is_lower = !B || is_upper || keepfirst ? 0: ZM_is_lower(B);
2389 2079251 : if (is_lower) L = RgM_flip(B);
2390 : }
2391 2139712 : if(DEBUGLEVEL>=4) timer_start(&T);
2392 2139711 : if (n > 2)
2393 : {
2394 1752782 : GEN R = B ? is_upper ? B : is_lower ? L : get_gramschmidt(B, rank) : get_gaussred(G, rank);
2395 1752771 : if (R)
2396 : {
2397 1127082 : long spr = spread(R), sz = mpexpo(gsupnorm(R, DEFAULTPREC)), thr;
2398 1127091 : if (DEBUGLEVEL>=5) err_printf("LLL: dim %ld, size %ld, spread %ld\n",n, sz, spr);
2399 1127091 : if ((is_upper && ZM_is_knapsack(B)) || (is_lower && ZM_is_knapsack(L)))
2400 94824 : thr = thsn[minss(n-3,numberof(thsn)-1)];
2401 : else
2402 : {
2403 1032267 : thr = thre[minss(n-3,numberof(thre)-1)];
2404 1032265 : if (n>=10) sz = spr;
2405 : }
2406 1127089 : useflatter = sz >= thr;
2407 : } else
2408 625689 : useflatter = 1;
2409 386929 : } else useflatter = 0;
2410 2139707 : if (useflatter)
2411 : {
2412 626875 : if (is_lower)
2413 : {
2414 0 : fplll_flatter(&G, &L, &U, rank, flag | LLL_UPPER);
2415 0 : B = RgM_flop(L);
2416 0 : if (U) U = RgM_flop(U);
2417 : }
2418 : else
2419 626875 : fplll_flatter(&G, &B, &U, rank, flag | (is_upper? LLL_UPPER:0));
2420 626876 : if (DEBUGLEVEL>=4 && !(flag & LLL_NOCERTIFY))
2421 0 : timer_printf(&T, "FLATTER");
2422 : }
2423 2139708 : if (!(flag & LLL_GRAM))
2424 : {
2425 : long t;
2426 2079247 : B = gcopy(B);
2427 2079258 : if(DEBUGLEVEL>=4)
2428 0 : err_printf("Entering L^2 (double): dim %ld, LLL-parameters (%.3f,%.3f)\n",
2429 : n, DELTA,ETA);
2430 2079258 : zeros = fplll_fast(&B, &U, DELTA, ETA, keepfirst);
2431 2079261 : if (DEBUGLEVEL>=4) timer_printf(&T, zeros < 0? "LLL (failed)": "LLL");
2432 2083514 : for (p = DEFAULTPREC, t = 0; zeros < 0 && t < 1 ; p += EXTRAPREC64, t++)
2433 : {
2434 4254 : if (DEBUGLEVEL>=4)
2435 0 : err_printf("Entering L^2 (heuristic): LLL-parameters (%.3f,%.3f), prec = %d/%d\n", DELTA, ETA, p, p);
2436 4254 : zeros = fplll_heuristic(&B, &U, DELTA, ETA, keepfirst, p, p);
2437 4254 : gc_lll(av, 2, &B, &U);
2438 4254 : if (DEBUGLEVEL>=4) timer_printf(&T, zeros < 0? "LLL (failed)": "LLL");
2439 : }
2440 : } else
2441 60461 : G = gcopy(G);
2442 2139722 : if (zeros < 0 || !(flag & LLL_NOCERTIFY))
2443 : {
2444 1574585 : if(DEBUGLEVEL>=4)
2445 0 : err_printf("Entering L^2 (dpe): LLL-parameters (%.3f,%.3f)\n", DELTA,ETA);
2446 1574585 : zeros = fplll_dpe(&G, &B, &U, pN, DELTA, ETA, keepfirst);
2447 1574571 : if (DEBUGLEVEL>=4) timer_printf(&T, zeros < 0? "LLL (failed)": "LLL");
2448 1574574 : if (zeros < 0)
2449 0 : for (p = DEFAULTPREC;; p += EXTRAPREC64)
2450 : {
2451 0 : if (DEBUGLEVEL>=4)
2452 0 : err_printf("Entering L^2: LLL-parameters (%.3f,%.3f), prec = %d\n",
2453 : DELTA,ETA, p);
2454 0 : zeros = fplll(&G, &B, &U, pN, DELTA, ETA, keepfirst, p);
2455 0 : if (DEBUGLEVEL>=4) timer_printf(&T, zeros < 0? "LLL (failed)": "LLL");
2456 0 : if (zeros >= 0) break;
2457 0 : gc_lll(av, 3, &G, &B, &U);
2458 : }
2459 : }
2460 2139711 : return lll_finish(U? U: B, zeros, flag);
2461 : }
2462 :
2463 : /********************************************************************/
2464 : /** **/
2465 : /** LLL OVER K[X] **/
2466 : /** **/
2467 : /********************************************************************/
2468 : static int
2469 504 : pslg(GEN x)
2470 : {
2471 : long tx;
2472 504 : if (gequal0(x)) return 2;
2473 448 : tx = typ(x); return is_scalar_t(tx)? 3: lg(x);
2474 : }
2475 :
2476 : static int
2477 196 : REDgen(long k, long l, GEN h, GEN L, GEN B)
2478 : {
2479 196 : GEN q, u = gcoeff(L,k,l);
2480 : long i;
2481 :
2482 196 : if (pslg(u) < pslg(B)) return 0;
2483 :
2484 140 : q = gneg(gdeuc(u,B));
2485 140 : gel(h,k) = gadd(gel(h,k), gmul(q,gel(h,l)));
2486 140 : for (i=1; i<l; i++) gcoeff(L,k,i) = gadd(gcoeff(L,k,i), gmul(q,gcoeff(L,l,i)));
2487 140 : gcoeff(L,k,l) = gadd(gcoeff(L,k,l), gmul(q,B)); return 1;
2488 : }
2489 :
2490 : static int
2491 196 : do_SWAPgen(GEN h, GEN L, GEN B, long k, GEN fl, int *flc)
2492 : {
2493 : GEN p1, la, la2, Bk;
2494 : long ps1, ps2, i, j, lx;
2495 :
2496 196 : if (!fl[k-1]) return 0;
2497 :
2498 140 : la = gcoeff(L,k,k-1); la2 = gsqr(la);
2499 140 : Bk = gel(B,k);
2500 140 : if (fl[k])
2501 : {
2502 56 : GEN q = gadd(la2, gmul(gel(B,k-1),gel(B,k+1)));
2503 56 : ps1 = pslg(gsqr(Bk));
2504 56 : ps2 = pslg(q);
2505 56 : if (ps1 <= ps2 && (ps1 < ps2 || !*flc)) return 0;
2506 28 : *flc = (ps1 != ps2);
2507 28 : gel(B,k) = gdiv(q, Bk);
2508 : }
2509 :
2510 112 : swap(gel(h,k-1), gel(h,k)); lx = lg(L);
2511 112 : for (j=1; j<k-1; j++) swap(gcoeff(L,k-1,j), gcoeff(L,k,j));
2512 112 : if (fl[k])
2513 : {
2514 28 : for (i=k+1; i<lx; i++)
2515 : {
2516 0 : GEN t = gcoeff(L,i,k);
2517 0 : p1 = gsub(gmul(gel(B,k+1),gcoeff(L,i,k-1)), gmul(la,t));
2518 0 : gcoeff(L,i,k) = gdiv(p1, Bk);
2519 0 : p1 = gadd(gmul(la,gcoeff(L,i,k-1)), gmul(gel(B,k-1),t));
2520 0 : gcoeff(L,i,k-1) = gdiv(p1, Bk);
2521 : }
2522 : }
2523 84 : else if (!gequal0(la))
2524 : {
2525 28 : p1 = gdiv(la2, Bk);
2526 28 : gel(B,k+1) = gel(B,k) = p1;
2527 28 : for (i=k+2; i<=lx; i++) gel(B,i) = gdiv(gmul(p1,gel(B,i)),Bk);
2528 28 : for (i=k+1; i<lx; i++)
2529 0 : gcoeff(L,i,k-1) = gdiv(gmul(la,gcoeff(L,i,k-1)), Bk);
2530 28 : for (j=k+1; j<lx-1; j++)
2531 0 : for (i=j+1; i<lx; i++)
2532 0 : gcoeff(L,i,j) = gdiv(gmul(p1,gcoeff(L,i,j)), Bk);
2533 : }
2534 : else
2535 : {
2536 56 : gcoeff(L,k,k-1) = gen_0;
2537 56 : for (i=k+1; i<lx; i++)
2538 : {
2539 0 : gcoeff(L,i,k) = gcoeff(L,i,k-1);
2540 0 : gcoeff(L,i,k-1) = gen_0;
2541 : }
2542 56 : gel(B,k) = gel(B,k-1); fl[k] = 1; fl[k-1] = 0;
2543 : }
2544 112 : return 1;
2545 : }
2546 :
2547 : static void
2548 168 : incrementalGSgen(GEN x, GEN L, GEN B, long k, GEN fl)
2549 : {
2550 168 : GEN u = NULL; /* gcc -Wall */
2551 : long i, j;
2552 420 : for (j = 1; j <= k; j++)
2553 252 : if (j==k || fl[j])
2554 : {
2555 252 : u = gcoeff(x,k,j);
2556 252 : if (!is_extscalar_t(typ(u))) pari_err_TYPE("incrementalGSgen",u);
2557 336 : for (i=1; i<j; i++)
2558 84 : if (fl[i])
2559 : {
2560 84 : u = gsub(gmul(gel(B,i+1),u), gmul(gcoeff(L,k,i),gcoeff(L,j,i)));
2561 84 : u = gdiv(u, gel(B,i));
2562 : }
2563 252 : gcoeff(L,k,j) = u;
2564 : }
2565 168 : if (gequal0(u)) gel(B,k+1) = gel(B,k);
2566 : else
2567 : {
2568 112 : gel(B,k+1) = gcoeff(L,k,k); gcoeff(L,k,k) = gen_1; fl[k] = 1;
2569 : }
2570 168 : }
2571 :
2572 : static GEN
2573 168 : lllgramallgen(GEN x, long flag)
2574 : {
2575 168 : long lx = lg(x), i, j, k, l, n;
2576 : pari_sp av;
2577 : GEN B, L, h, fl;
2578 : int flc;
2579 :
2580 168 : n = lx-1; if (n<=1) return lll_trivial(x,flag);
2581 84 : if (lgcols(x) != lx) pari_err_DIM("lllgramallgen");
2582 :
2583 84 : fl = cgetg(lx, t_VECSMALL);
2584 :
2585 84 : av = avma;
2586 84 : B = scalarcol_shallow(gen_1, lx);
2587 84 : L = cgetg(lx,t_MAT);
2588 252 : for (j=1; j<lx; j++) { gel(L,j) = zerocol(n); fl[j] = 0; }
2589 :
2590 84 : h = matid(n);
2591 252 : for (i=1; i<lx; i++)
2592 168 : incrementalGSgen(x, L, B, i, fl);
2593 84 : flc = 0;
2594 84 : for(k=2;;)
2595 : {
2596 196 : if (REDgen(k, k-1, h, L, gel(B,k))) flc = 1;
2597 196 : if (do_SWAPgen(h, L, B, k, fl, &flc)) { if (k > 2) k--; }
2598 : else
2599 : {
2600 84 : for (l=k-2; l>=1; l--)
2601 0 : if (REDgen(k, l, h, L, gel(B,l+1))) flc = 1;
2602 84 : if (++k > n) break;
2603 : }
2604 112 : if (gc_needed(av,1))
2605 : {
2606 0 : if(DEBUGMEM>1) pari_warn(warnmem,"lllgramallgen");
2607 0 : gerepileall(av,3,&B,&L,&h);
2608 : }
2609 : }
2610 140 : k=1; while (k<lx && !fl[k]) k++;
2611 84 : return lll_finish(h,k-1,flag);
2612 : }
2613 :
2614 : static GEN
2615 168 : lllallgen(GEN x, long flag)
2616 : {
2617 168 : pari_sp av = avma;
2618 168 : if (!(flag & LLL_GRAM)) x = gram_matrix(x);
2619 84 : else if (!RgM_is_square_mat(x)) pari_err_DIM("qflllgram");
2620 168 : return gerepilecopy(av, lllgramallgen(x, flag));
2621 : }
2622 : GEN
2623 42 : lllgen(GEN x) { return lllallgen(x, LLL_IM); }
2624 : GEN
2625 42 : lllkerimgen(GEN x) { return lllallgen(x, LLL_ALL); }
2626 : GEN
2627 42 : lllgramgen(GEN x) { return lllallgen(x, LLL_IM|LLL_GRAM); }
2628 : GEN
2629 42 : lllgramkerimgen(GEN x) { return lllallgen(x, LLL_ALL|LLL_GRAM); }
2630 :
2631 : static GEN
2632 56109 : lllall(GEN x, long flag)
2633 56109 : { pari_sp av = avma; return gerepilecopy(av, ZM_lll(x, LLLDFT, flag)); }
2634 : GEN
2635 14693 : lllint(GEN x) { return lllall(x, LLL_IM); }
2636 : GEN
2637 35 : lllkerim(GEN x) { return lllall(x, LLL_ALL); }
2638 : GEN
2639 41339 : lllgramint(GEN x)
2640 41339 : { if (!RgM_is_square_mat(x)) pari_err_DIM("qflllgram");
2641 41339 : return lllall(x, LLL_IM | LLL_GRAM); }
2642 : GEN
2643 35 : lllgramkerim(GEN x)
2644 35 : { if (!RgM_is_square_mat(x)) pari_err_DIM("qflllgram");
2645 35 : return lllall(x, LLL_ALL | LLL_GRAM); }
2646 :
2647 : GEN
2648 5322387 : lllfp(GEN x, double D, long flag)
2649 : {
2650 5322387 : long n = lg(x)-1;
2651 5322387 : pari_sp av = avma;
2652 : GEN h;
2653 5322387 : if (n <= 1) return lll_trivial(x,flag);
2654 4662579 : if ((flag & LLL_GRAM) && !RgM_is_square_mat(x)) pari_err_DIM("qflllgram");
2655 4662565 : h = ZM_lll(RgM_rescale_to_int(x), D, flag);
2656 4662541 : return gerepilecopy(av, h);
2657 : }
2658 :
2659 : GEN
2660 8942 : lllgram(GEN x) { return lllfp(x,LLLDFT,LLL_GRAM|LLL_IM); }
2661 : GEN
2662 1173380 : lll(GEN x) { return lllfp(x,LLLDFT,LLL_IM); }
2663 :
2664 : GEN
2665 301 : qflll0(GEN x, long flag)
2666 : {
2667 301 : if (typ(x) != t_MAT) pari_err_TYPE("qflll",x);
2668 301 : switch(flag)
2669 : {
2670 49 : case 0: return lll(x);
2671 63 : case 1: RgM_check_ZM(x,"qflll"); return lllint(x);
2672 49 : case 2: RgM_check_ZM(x,"qflll"); return lllintpartial(x);
2673 7 : case 3: RgM_check_ZM(x,"qflll"); return lllall(x, LLL_INPLACE);
2674 49 : case 4: RgM_check_ZM(x,"qflll"); return lllkerim(x);
2675 42 : case 5: return lllkerimgen(x);
2676 42 : case 8: return lllgen(x);
2677 0 : default: pari_err_FLAG("qflll");
2678 : }
2679 : return NULL; /* LCOV_EXCL_LINE */
2680 : }
2681 :
2682 : GEN
2683 245 : qflllgram0(GEN x, long flag)
2684 : {
2685 245 : if (typ(x) != t_MAT) pari_err_TYPE("qflllgram",x);
2686 245 : switch(flag)
2687 : {
2688 63 : case 0: return lllgram(x);
2689 49 : case 1: RgM_check_ZM(x,"qflllgram"); return lllgramint(x);
2690 49 : case 4: RgM_check_ZM(x,"qflllgram"); return lllgramkerim(x);
2691 42 : case 5: return lllgramkerimgen(x);
2692 42 : case 8: return lllgramgen(x);
2693 0 : default: pari_err_FLAG("qflllgram");
2694 : }
2695 : return NULL; /* LCOV_EXCL_LINE */
2696 : }
2697 :
2698 : /********************************************************************/
2699 : /** **/
2700 : /** INTEGRAL KERNEL (LLL REDUCED) **/
2701 : /** **/
2702 : /********************************************************************/
2703 : static GEN
2704 70 : kerint0(GEN M)
2705 : {
2706 : /* return ZM_lll(M, LLLDFT, LLL_KER); */
2707 70 : GEN U, H = ZM_hnflll(M,&U,1);
2708 70 : long d = lg(M)-lg(H);
2709 70 : if (!d) return cgetg(1, t_MAT);
2710 70 : return ZM_lll(vecslice(U,1,d), LLLDFT, LLL_INPLACE);
2711 : }
2712 : GEN
2713 42 : kerint(GEN M)
2714 : {
2715 42 : pari_sp av = avma;
2716 42 : return gerepilecopy(av, kerint0(M));
2717 : }
2718 : /* OBSOLETE: use kerint */
2719 : GEN
2720 28 : matkerint0(GEN M, long flag)
2721 : {
2722 28 : pari_sp av = avma;
2723 28 : if (typ(M) != t_MAT) pari_err_TYPE("matkerint",M);
2724 28 : M = Q_primpart(M);
2725 28 : RgM_check_ZM(M, "kerint");
2726 28 : switch(flag)
2727 : {
2728 28 : case 0:
2729 28 : case 1: return gerepilecopy(av, kerint0(M));
2730 0 : default: pari_err_FLAG("matkerint");
2731 : }
2732 : return NULL; /* LCOV_EXCL_LINE */
2733 : }
|
