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 45834 : RgM_is_square_mat(GEN x) { long l = lg(x); return l == 1 || l == lgcols(x); }
22 :
23 : static long
24 4239494 : ZM_is_upper(GEN R)
25 : {
26 4239494 : long i,j, l = lg(R);
27 4239494 : if (l != lgcols(R)) return 0;
28 8195982 : for(i = 1; i < l; i++)
29 8908814 : for(j = 1; j < i; j++)
30 4602955 : if (signe(gcoeff(R,i,j))) return 0;
31 265375 : return 1;
32 : }
33 :
34 : static long
35 607498 : ZM_is_knapsack(GEN R)
36 : {
37 607498 : long i,j, l = lg(R);
38 607498 : if (l != lgcols(R)) return 0;
39 846687 : for(i = 2; i < l; i++)
40 2921459 : for(j = 1; j < l; j++)
41 2682270 : if ( i!=j && signe(gcoeff(R,i,j))) return 0;
42 92869 : return 1;
43 : }
44 :
45 : static long
46 1205606 : ZM_is_lower(GEN R)
47 : {
48 1205606 : long i,j, l = lg(R);
49 1205606 : if (l != lgcols(R)) return 0;
50 2090771 : for(i = 1; i < l; i++)
51 2420099 : for(j = 1; j < i; j++)
52 1311881 : if (signe(gcoeff(R,j,i))) return 0;
53 34953 : return 1;
54 : }
55 :
56 : static GEN
57 34953 : RgM_flip(GEN R)
58 : {
59 : GEN M;
60 : long i,j,l;
61 34953 : M = cgetg_copy(R, &l);
62 181994 : for(i = 1; i < l; i++)
63 : {
64 147041 : gel(M,i) = cgetg(l, t_COL);
65 921196 : for(j = 1; j < l; j++)
66 774155 : gmael(M,i,j) = gmael(R,l-i, l-j);
67 : }
68 34953 : 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 4108859 : mpabscmp(GEN x, GEN y)
89 : {
90 4108859 : return (typ(x)==t_INT) ? abscmpii(x,y) : abscmprr(x,y);
91 : }
92 :
93 : /****************************************************************************/
94 : /*** FLATTER ***/
95 : /****************************************************************************/
96 : /* Implementation of "FLATTER" algorithm based on
97 : * <https://eprint.iacr.org/2023/237>
98 : * Fast Practical Lattice Reduction through Iterated Compression
99 : *
100 : * Keegan Ryan, University of California, San Diego
101 : * Nadia Heninger, University of California, San Diego. BA20230925 */
102 : static long
103 1347760 : drop(GEN R)
104 : {
105 1347760 : long i, n = lg(R)-1;
106 1347760 : long s = 0, m = mpexpo(gcoeff(R, 1, 1));
107 5456619 : for (i = 2; i <= n; ++i)
108 : {
109 4108859 : if (mpabscmp(gcoeff(R, i, i), gcoeff(R, i - 1, i - 1)) >= 0)
110 : {
111 2786358 : s += m - mpexpo(gcoeff(R, i - 1, i - 1));
112 2786358 : m = mpexpo(gcoeff(R, i, i));
113 : }
114 : }
115 1347760 : s += m - mpexpo(gcoeff(R, n, n));
116 1347760 : return s;
117 : }
118 :
119 : static long
120 1347760 : potential(GEN R)
121 : {
122 1347760 : long i, n = lg(R)-1;
123 1347760 : long s = 0, mul = n-1;;
124 6804379 : for (i = 1; i <= n; i++, mul-=2) s += mul * mpexpo(gcoeff(R,i,i));
125 1347760 : return s;
126 : }
127 :
128 : /* U upper-triangular invertible:
129 : * Bound on the exponent of the condition number of U.
130 : * Algo 8.13 in Higham, Accuracy and stability of numercal algorithms. */
131 : static long
132 4729021 : condition_bound(GEN U, int lower)
133 : {
134 4729021 : long n = lg(U)-1, e, i, j;
135 : GEN y;
136 4729021 : pari_sp av = avma;
137 4729021 : y = cgetg(n+1, t_VECSMALL);
138 4729021 : e = y[n] = -gexpo(gcoeff(U,n,n));
139 18836195 : for (i=n-1; i>0; i--)
140 : {
141 14107174 : long s = 0;
142 50966101 : for (j=i+1; j<=n; j++)
143 36858927 : s = maxss(s, (lower? gexpo(gcoeff(U,j,i)): gexpo(gcoeff(U,i,j))) + y[j]);
144 14107174 : y[i] = s - gexpo(gcoeff(U,i,i));
145 14107174 : e = maxss(e, y[i]);
146 : }
147 4729021 : return gc_long(av, gexpo(U) + e);
148 : }
149 :
150 : INLINE long
151 7496219 : nbits2prec64(long n)
152 : {
153 7496219 : return nbits2prec(((n+63)>>6)<<6);
154 : }
155 :
156 : static long
157 5855635 : spread(GEN R)
158 : {
159 5855635 : long i, n = lg(R)-1, m = mpexpo(gcoeff(R, 1, 1)), M = m;
160 23619003 : for (i = 2; i <= n; ++i)
161 : {
162 17763368 : long e = mpexpo(gcoeff(R, i, i));
163 17763368 : if (e < m) m = e;
164 17763368 : if (e > M) M = e;
165 : }
166 5855635 : return M - m;
167 : }
168 :
169 : static long
170 4729021 : GS_extraprec(GEN L, int lower)
171 : {
172 4729021 : long C = condition_bound(L, lower), S = spread(L), n = lg(L)-1;
173 4729021 : return maxss(2*S+2*n, C-S-2*n); /* = 2*S + 2*n + maxss(0, C-3*S-4*n) */
174 : }
175 :
176 : static GEN
177 2988 : RgM_Cholesky_dynprec(GEN M)
178 : {
179 2988 : pari_sp ltop = avma;
180 : GEN L;
181 2988 : long minprec = lg(M) + 30, bitprec = minprec, prec;
182 : while (1)
183 4918 : {
184 : long mbitprec;
185 7906 : prec = nbits2prec64(bitprec);
186 7906 : L = RgM_Cholesky(RgM_gtofp(M, prec), prec); /* upper-triangular */
187 7906 : if (!L)
188 : {
189 1486 : bitprec *= 2;
190 1486 : set_avma(ltop);
191 1486 : continue;
192 : }
193 6420 : mbitprec = minprec + GS_extraprec(L, 0);
194 6420 : if (bitprec >= mbitprec)
195 2988 : break;
196 3432 : bitprec = maxss((4*bitprec)/3, mbitprec);
197 3432 : set_avma(ltop);
198 : }
199 2988 : return gc_GEN(ltop, L);
200 : }
201 :
202 : static GEN
203 1402 : gramschmidt_upper(GEN M)
204 : {
205 1402 : long bitprec = lg(M)-1 + 31 + GS_extraprec(M, 0);
206 1402 : return RgM_gtofp(M, nbits2prec64(bitprec));
207 : }
208 :
209 : static GEN
210 2695520 : gramschmidt_dynprec(GEN M)
211 : {
212 2695520 : pari_sp ltop = avma;
213 2695520 : long minprec = lg(M) + 30, bitprec = minprec;
214 2695520 : if (ZM_is_upper(M)) return gramschmidt_upper(M);
215 : while (1)
216 3648605 : {
217 : GEN B, Q, L;
218 6342723 : long prec = nbits2prec64(bitprec), mbitprec;
219 6342723 : if (!QR_init(RgM_gtofp(M, prec), &B, &Q, &L, prec))
220 : {
221 1621524 : bitprec *= 2;
222 1621524 : set_avma(ltop);
223 1621524 : continue;
224 : }
225 4721199 : mbitprec = minprec + GS_extraprec(L, 1);
226 4721199 : if (bitprec >= mbitprec)
227 2694118 : return gc_GEN(ltop, shallowtrans(L));
228 2027081 : bitprec = maxss((4*bitprec)/3, mbitprec);
229 2027081 : set_avma(ltop);
230 : }
231 : }
232 : /* return -T1 * round(T1^-1*(R1^-1*R2)*T3) */
233 : static GEN
234 1347760 : sizered(GEN T1, GEN T3, GEN R1, GEN R2)
235 : {
236 1347760 : pari_sp ltop = avma;
237 : long e;
238 1347760 : return gc_upto(ltop, ZM_mul(ZM_neg(T1), grndtoi(gmul(ZM_inv(T1,NULL),
239 : RgM_mul(RgM_mul(RgM_inv_upper(R1), R2), T3)), &e)));
240 : }
241 :
242 : static GEN
243 1347760 : flat(GEN M, long flag, GEN *pt_T, long *pt_s, long *pt_pot)
244 : {
245 1347760 : pari_sp ltop = avma;
246 : GEN R, R1, R2, R3, T1, T2, T3, T, S;
247 1347760 : long k = lg(M)-1, n = k>>1, n2 = k - n, m = n>>1;
248 1347760 : long keepfirst = flag & LLL_KEEP_FIRST, inplace = flag & LLL_INPLACE;
249 : /* for k = 3, we want n = 1; n2 = 2; m = 0 */
250 : /* for k = 5, n = 2; n2 = 3; m = 1 */
251 1347760 : R = gramschmidt_dynprec(M);
252 1347760 : R1 = matslice(R, 1, n, 1, n);
253 1347760 : R2 = matslice(R, 1, n, n + 1, k);
254 1347760 : R3 = matslice(R, n + 1, k, n + 1, k);
255 1347760 : T1 = lllfp(R1, 0.99, LLL_IM| LLL_UPPER| LLL_NOCERTIFY| (keepfirst ? LLL_KEEP_FIRST: 0));
256 1347760 : T3 = lllfp(R3, 0.99, LLL_IM| LLL_UPPER| LLL_NOCERTIFY);
257 1347760 : T2 = sizered(T1, T3, R1, R2);
258 1347760 : T = shallowmatconcat(mkmat22(T1,T2,gen_0,T3));
259 1347760 : M = ZM_mul(M, T);
260 1347760 : R = gramschmidt_dynprec(M);
261 1347760 : R3 = matslice(R, m + 1, m + n2, m + 1, m + n2);
262 1347760 : T3 = lllfp(R3, 0.99, LLL_IM| LLL_UPPER| LLL_NOCERTIFY);
263 2695520 : S = shallowmatconcat(diagonal(
264 577282 : m == 0 ? mkvec2(T3, matid(k - m - n2))
265 0 : : m+n2 == k ? mkvec2(matid(m), T3)
266 770478 : : mkvec3(matid(m), T3, matid(k - m - n2))));
267 1347760 : M = ZM_mul(M, S);
268 1347760 : if (!inplace) *pt_T = ZM_mul(T, S);
269 1347760 : *pt_s = drop(R);
270 1347760 : *pt_pot = potential(R);
271 1347760 : return gc_all(ltop, inplace ? 1: 2, &M, pt_T);
272 : }
273 :
274 : static void
275 0 : dbg_flatter(pari_timer *ti, long n, long i, long lti, double t, double pot2)
276 : {
277 0 : double s = t / n, p = pot2 / (n*(n+1));
278 : const char *str;
279 0 : if (i == -1)
280 0 : str = (i == lti)? "final"
281 0 : : stack_sprintf("steps %ld-final", lti);
282 : else
283 0 : str = (i == lti)? stack_sprintf("step %ld", i)
284 0 : : stack_sprintf("steps %ld-%ld", lti, i);
285 0 : timer_printf(ti, "FLATTER, dim %ld, %s: \t slope=%0.10g \t pot=%0.10g",
286 : n, str, s, p);
287 0 : }
288 :
289 : static GEN
290 627251 : ZM_flatter(GEN M, long flag)
291 : {
292 627251 : pari_sp av = avma;
293 627251 : long i, n = lg(M)-1, s = -1, lti = 1, pot = LONG_MAX;
294 627251 : GEN T = NULL;
295 : pari_timer ti;
296 627251 : long inplace = flag & LLL_INPLACE, cert = !(flag & LLL_NOCERTIFY);
297 :
298 627251 : if (DEBUGLEVEL>=3)
299 : {
300 0 : timer_start(&ti);
301 0 : if (cert) err_printf("FLATTER dim = %ld size = %ld\n", n, ZM_max_expi(M));
302 : }
303 627251 : for (i = 1;;i++)
304 720509 : {
305 : long t, pot2;
306 1347760 : GEN U, M2 = flat(M, flag, &U, &t, &pot2);
307 1347760 : if (t == 0) { s = t; break; }
308 764322 : if (s >= 0)
309 : {
310 437764 : if (s == t && pot>=pot2) break;
311 393951 : if (s < t && i > 20)
312 : {
313 0 : if (DEBUGLEVEL >= 3) err_printf("BACK:%ld:%ld:%g\n", n, i, s);
314 0 : break;
315 : }
316 : }
317 720509 : if (DEBUGLEVEL>=3 && (cert || timer_get(&ti) > 1000))
318 0 : dbg_flatter(&ti, n, i, lti, t, pot2);
319 720509 : s = t;
320 720509 : pot = pot2;
321 720509 : M = M2;
322 720509 : if (!inplace)
323 : {
324 692825 : T = T? ZM_mul(T, U): U;
325 692825 : if (gc_needed(av, 1)) (void)gc_all(av, 2, &M, &T);
326 : }
327 : else
328 27684 : if (gc_needed(av, 1)) M = gc_GEN(av, M);
329 : }
330 627251 : if (DEBUGLEVEL>=3 && (cert || timer_get(&ti) > 1000))
331 0 : dbg_flatter(&ti, n, -1, i == lti? -1: lti, s, pot);
332 627251 : if (!inplace)
333 : {
334 613236 : if (!T) return gc_NULL(av);
335 312711 : return gc_GEN(av, T);
336 : }
337 14015 : return gc_GEN(av, M);
338 : }
339 :
340 : static GEN
341 625237 : ZM_flatter_rank(GEN M, long rank, long flag)
342 : {
343 : pari_timer ti;
344 625237 : pari_sp av = avma;
345 625237 : GEN T = NULL;
346 625237 : long i, n = lg(M)-1, sm = LONG_MAX;
347 625237 : long inplace = flag & LLL_INPLACE;
348 :
349 625237 : if (rank == n) return ZM_flatter(M, flag);
350 3785 : if (DEBUGLEVEL>=3) timer_start(&ti);
351 3785 : for (i = 1;; i++)
352 2014 : {
353 5799 : GEN S = ZM_flatter(vconcat(gshift(M,i),matid(n)), flag);
354 : long s;
355 5799 : if (!S || (s = expi(gnorml2(S))) >= sm) break;
356 2014 : sm = s;
357 2014 : if (DEBUGLEVEL>=3) timer_printf(&ti,"FLATTERRANK step %ld: %ld",i,sm);
358 2014 : T = T? ZM_mul(T, S): S;
359 2014 : M = ZM_mul(M, S);
360 2014 : if (gc_needed(av, 1)) (void)gc_all(av, 2, &M, &T);
361 : }
362 3785 : if (!inplace)
363 : {
364 3778 : if (!T) { set_avma(av); return matid(n); }
365 1951 : return gc_GEN(av, T);
366 : }
367 7 : return gc_GEN(av, M);
368 : }
369 :
370 : static GEN
371 2988 : flattergram_i(GEN M, long flag)
372 : {
373 2988 : pari_sp av = avma;
374 2988 : GEN T, R = RgM_Cholesky_dynprec(M);
375 2988 : T = lllfp(R, 0.99, LLL_IM|LLL_UPPER|LLL_NOCERTIFY | (flag&LLL_KEEP_FIRST));
376 2988 : return gc_upto(av, T);
377 : }
378 :
379 : static void
380 0 : dbg_flattergram(pari_timer *t, long n, long i, long s)
381 0 : { timer_printf(t, "FLATTERGRAM, dim %ld step %ld, slope=%0.10g", n, i,
382 0 : ((double)s)/n); }
383 : /* return base change, NULL if identity */
384 : static GEN
385 968 : ZM_flattergram(GEN M, long flag)
386 : {
387 968 : pari_sp av = avma;
388 968 : GEN T = NULL;
389 968 : long i, n = lg(M)-1, s = -1;
390 :
391 : pari_timer ti;
392 968 : if (DEBUGLEVEL>=3)
393 : {
394 0 : timer_start(&ti);
395 0 : err_printf("FLATTERGRAM dim = %ld size = %ld\n", n, ZM_max_expi(M));
396 : }
397 968 : for (i = 1;; i++)
398 2020 : {
399 2988 : GEN S = flattergram_i(M, flag);
400 2988 : long t = expi(gnorml2(S));
401 2988 : if (t == 0) { s = t; break; }
402 2988 : if (s)
403 : {
404 2988 : double st = s - t;
405 2988 : if (st == 0) break;
406 2020 : if (st < 0 && i > 20)
407 : {
408 0 : if (DEBUGLEVEL >= 3)
409 0 : err_printf("BACK:%ld:%ld:%0.10g\n", n, i, ((double)s)/n);
410 0 : break;
411 : }
412 : }
413 2020 : T = T? ZM_mul(T, S): S;
414 2020 : M = qf_ZM_apply(M, S);
415 2020 : s = t;
416 2020 : if (DEBUGLEVEL >= 3) dbg_flattergram(&ti, n, i, s);
417 2020 : if (gc_needed(av, 1)) (void)gc_all(av, 2, &M, &T);
418 : }
419 968 : if (DEBUGLEVEL >= 3) dbg_flattergram(&ti, n, i, s);
420 968 : if (!T && ZM_isidentity(T)) return gc_NULL(av);
421 968 : return gc_GEN(av, T);
422 : }
423 :
424 : /* return base change, NULL if identity */
425 : static GEN
426 968 : ZM_flattergram_rank(GEN M, long rank, long flag)
427 : {
428 : pari_timer ti;
429 968 : pari_sp av = avma;
430 968 : GEN T = NULL;
431 968 : long i, n = lg(M)-1;
432 968 : if (rank == n) return ZM_flattergram(M, flag);
433 0 : if (DEBUGLEVEL>=3) timer_start(&ti);
434 0 : for (i = 1;; i++)
435 0 : {
436 0 : GEN S = ZM_flattergram(RgM_Rg_add(gshift(M, i), gen_1), flag);
437 0 : if (DEBUGLEVEL>=3)
438 0 : timer_printf(&ti,"FLATTERGRAMRANK step %ld: %ld",i,expi(gnorml2(S)));
439 0 : if (!S) break;
440 0 : T = T? ZM_mul(T, S): S;
441 0 : M = qf_ZM_apply(M, S);
442 0 : if (gc_needed(av, 1)) (void)gc_all(av, 2, &M, &T);
443 : }
444 0 : if (!T || ZM_isidentity(T)) return gc_NULL(av);
445 0 : return gc_GEN(av, T);
446 : }
447 :
448 : /* round to closest integer (as a double). If |a| >= 2^52, return it */
449 : static double
450 11656149 : pari_rint(double a)
451 : {
452 : #ifdef HAS_RINT
453 11656149 : return rint(a);
454 : #else
455 : const double pow2 = 4.5035996273704960e+15; /* 2^52 */
456 : double r, fa = fabs(a);
457 : if (fa >= pow2) return a;
458 : r = (pow2 + fa) - pow2;
459 : if (a < 0) r = -r;
460 : return r;
461 : #endif
462 : }
463 :
464 : /* default quality ratio for LLL */
465 : static const double LLLDFT = 0.99;
466 :
467 : /* assume flag & (LLL_KER|LLL_IM|LLL_ALL). LLL_INPLACE implies LLL_IM */
468 : static GEN
469 771847 : lll_trivial(GEN x, long flag)
470 : {
471 771847 : if (lg(x) == 1)
472 : { /* dim x = 0 */
473 15484 : if (! (flag & LLL_ALL)) return cgetg(1,t_MAT);
474 28 : retmkvec2(cgetg(1,t_MAT), cgetg(1,t_MAT));
475 : }
476 : /* dim x = 1 */
477 756363 : if (gequal0(gel(x,1)))
478 : {
479 151 : if (flag & LLL_KER) return matid(1);
480 151 : if (flag & (LLL_IM|LLL_INPLACE)) return cgetg(1,t_MAT);
481 28 : retmkvec2(matid(1), cgetg(1,t_MAT));
482 : }
483 756212 : if (flag & LLL_INPLACE) return gcopy(x);
484 652528 : if (flag & LLL_KER) return cgetg(1,t_MAT);
485 652528 : if (flag & LLL_IM) return matid(1);
486 28 : retmkvec2(cgetg(1,t_MAT), (flag & LLL_GRAM)? gcopy(x): matid(1));
487 : }
488 :
489 : /* vecslice(x,#x-k,#x) in place. Works for t_MAT, t_VEC/t_COL */
490 : static GEN
491 2093445 : vectail_inplace(GEN x, long k)
492 : {
493 2093445 : if (!k) return x;
494 58105 : x[k] = ((ulong)x[0] & ~LGBITS) | _evallg(lg(x) - k);
495 58105 : return x + k;
496 : }
497 :
498 : /* k = dim Kernel */
499 : static GEN
500 2168064 : lll_finish(GEN h, long k, long flag)
501 : {
502 : GEN g;
503 2168064 : if (!(flag & (LLL_IM|LLL_KER|LLL_ALL|LLL_INPLACE))) return h;
504 2093459 : if (flag & (LLL_IM|LLL_INPLACE)) return vectail_inplace(h, k);
505 84 : if (flag & LLL_KER) { setlg(h,k+1); return h; }
506 70 : g = vecslice(h,1,k); /* done first: vectail_inplace kills h */
507 70 : return mkvec2(g, vectail_inplace(h, k));
508 : }
509 :
510 : /* y * z * 2^e, e >= 0; y,z t_INT */
511 : INLINE GEN
512 938305 : mulshift(GEN y, GEN z, long e)
513 : {
514 938305 : long ly = lgefint(y), lz;
515 : pari_sp av;
516 : GEN t;
517 938305 : if (ly == 2) return gen_0;
518 456035 : lz = lgefint(z);
519 456035 : av = avma; (void)new_chunk(ly+lz+nbits2lg(e)); /* HACK */
520 456035 : t = mulii(z, y);
521 456035 : set_avma(av); return shifti(t, e);
522 : }
523 :
524 : /* x - y * z * 2^e, e >= 0; x,y,z t_INT */
525 : INLINE GEN
526 2074362 : submulshift(GEN x, GEN y, GEN z, long e)
527 : {
528 2074362 : long lx = lgefint(x), ly, lz;
529 : pari_sp av;
530 : GEN t;
531 2074362 : if (!e) return submulii(x, y, z);
532 2051874 : if (lx == 2) { t = mulshift(y, z, e); togglesign(t); return t; }
533 1538444 : ly = lgefint(y);
534 1538444 : if (ly == 2) return icopy(x);
535 1093910 : lz = lgefint(z);
536 1093910 : av = avma; (void)new_chunk(lx+ly+lz+nbits2lg(e)); /* HACK */
537 1093910 : t = shifti(mulii(z, y), e);
538 1093910 : set_avma(av); return subii(x, t);
539 : }
540 : static void
541 33271522 : subzi(GEN *a, GEN b)
542 : {
543 33271522 : pari_sp av = avma;
544 33271522 : b = subii(*a, b);
545 33271522 : if (lgefint(b)<=lg(*a) && isonstack(*a)) { affii(b,*a); set_avma(av); }
546 2429530 : else *a = b;
547 33271522 : }
548 :
549 : static void
550 32506340 : addzi(GEN *a, GEN b)
551 : {
552 32506340 : pari_sp av = avma;
553 32506340 : b = addii(*a, b);
554 32506340 : if (lgefint(b)<=lg(*a) && isonstack(*a)) { affii(b,*a); set_avma(av); }
555 2212338 : else *a = b;
556 32506340 : }
557 :
558 : /* x - u*y * 2^e */
559 : INLINE GEN
560 4742023 : submuliu2n(GEN x, GEN y, ulong u, long e)
561 : {
562 : pari_sp av;
563 4742023 : long ly = lgefint(y);
564 4742023 : if (ly == 2) return x;
565 3316703 : av = avma;
566 3316703 : (void)new_chunk(3+ly+lgefint(x)+nbits2lg(e)); /* HACK */
567 3316703 : y = shifti(mului(u,y), e);
568 3316703 : set_avma(av); return subii(x, y);
569 : }
570 : /* *x -= u*y * 2^e */
571 : INLINE void
572 16867897 : submulzu2n(GEN *x, GEN y, ulong u, long e)
573 : {
574 : pari_sp av;
575 16867897 : long ly = lgefint(y);
576 16867897 : if (ly == 2) return;
577 5891164 : av = avma;
578 5891164 : (void)new_chunk(3+ly+lgefint(*x)+nbits2lg(e)); /* HACK */
579 5891164 : y = shifti(mului(u,y), e);
580 5891164 : set_avma(av); return subzi(x, y);
581 : }
582 :
583 : /* x + u*y * 2^e */
584 : INLINE GEN
585 4666438 : addmuliu2n(GEN x, GEN y, ulong u, long e)
586 : {
587 : pari_sp av;
588 4666438 : long ly = lgefint(y);
589 4666438 : if (ly == 2) return x;
590 3276445 : av = avma;
591 3276445 : (void)new_chunk(3+ly+lgefint(x)+nbits2lg(e)); /* HACK */
592 3276445 : y = shifti(mului(u,y), e);
593 3276445 : set_avma(av); return addii(x, y);
594 : }
595 :
596 : /* *x += u*y * 2^e */
597 : INLINE void
598 17070592 : addmulzu2n(GEN *x, GEN y, ulong u, long e)
599 : {
600 : pari_sp av;
601 17070592 : long ly = lgefint(y);
602 17070592 : if (ly == 2) return;
603 5924878 : av = avma;
604 5924878 : (void)new_chunk(3+ly+lgefint(*x)+nbits2lg(e)); /* HACK */
605 5924878 : y = shifti(mului(u,y), e);
606 5924878 : set_avma(av); return addzi(x, y);
607 : }
608 :
609 : /* n < 10; (void)gc_all supporting &NULL arguments. Maybe rename and export ? */
610 : INLINE void
611 5460 : gc_lll(pari_sp av, int n, ...)
612 : {
613 : int i, j;
614 : GEN *gptr[10];
615 : size_t s;
616 5460 : va_list a; va_start(a, n);
617 16380 : for (i=j=0; i<n; i++)
618 : {
619 10920 : GEN *x = va_arg(a,GEN*);
620 10920 : if (*x) { gptr[j++] = x; *x = (GEN)copy_bin(*x); }
621 : }
622 5460 : va_end(a); set_avma(av);
623 13462 : for (--j; j>=0; j--) *gptr[j] = bin_copy((GENbin*)*gptr[j]);
624 5460 : s = pari_mainstack->top - pari_mainstack->bot;
625 : /* size of saved objects ~ stacksize / 4 => overflow */
626 5460 : if (av - avma > (s >> 2))
627 : {
628 0 : size_t t = avma - pari_mainstack->bot;
629 0 : av = avma; new_chunk((s + t) / sizeof(long)); set_avma(av); /* double */
630 : }
631 5460 : }
632 :
633 : /********************************************************************/
634 : /** **/
635 : /** FPLLL (adapted from D. Stehle's code) **/
636 : /** **/
637 : /********************************************************************/
638 : /* Babai* and fplll* are a conversion to libpari API and data types
639 : of fplll-1.3 by Damien Stehle'.
640 :
641 : Copyright 2005, 2006 Damien Stehle'.
642 :
643 : This program is free software; you can redistribute it and/or modify it
644 : under the terms of the GNU General Public License as published by the
645 : Free Software Foundation; either version 2 of the License, or (at your
646 : option) any later version.
647 :
648 : This program implements ideas from the paper "Floating-point LLL Revisited",
649 : by Phong Nguyen and Damien Stehle', in the Proceedings of Eurocrypt'2005,
650 : Springer-Verlag; and was partly inspired by Shoup's NTL library:
651 : http://www.shoup.net/ntl/ */
652 :
653 : /* x t_REAL, |x| >= 1/2. Test whether |x| <= 3/2 */
654 : static int
655 442325 : absrsmall2(GEN x)
656 : {
657 442325 : long e = expo(x), l, i;
658 442325 : if (e < 0) return 1;
659 230847 : if (e > 0 || (ulong)x[2] > (3UL << (BITS_IN_LONG-2))) return 0;
660 : /* line above assumes l > 2. OK since x != 0 */
661 79835 : l = lg(x); for (i = 3; i < l; i++) if (x[i]) return 0;
662 68365 : return 1;
663 : }
664 : /* x t_REAL; test whether |x| <= 1/2 */
665 : static int
666 761889 : absrsmall(GEN x)
667 : {
668 : long e, l, i;
669 761889 : if (!signe(x)) return 1;
670 755848 : e = expo(x); if (e < -1) return 1;
671 448631 : if (e > -1 || (ulong)x[2] > HIGHBIT) return 0;
672 7148 : l = lg(x); for (i = 3; i < l; i++) if (x[i]) return 0;
673 6306 : return 1;
674 : }
675 :
676 : static void
677 33254711 : rotate(GEN A, long k2, long k)
678 : {
679 : long i;
680 33254711 : GEN B = gel(A,k2);
681 107860127 : for (i = k2; i > k; i--) gel(A,i) = gel(A,i-1);
682 33254711 : gel(A,k) = B;
683 33254711 : }
684 :
685 : /************************* FAST version (double) ************************/
686 : #define dmael(x,i,j) ((x)[i][j])
687 : #define del(x,i) ((x)[i])
688 :
689 : static double *
690 35108508 : cget_dblvec(long d)
691 35108508 : { return (double*) stack_malloc_align(d*sizeof(double), sizeof(double)); }
692 :
693 : static double **
694 8429696 : cget_dblmat(long d) { return (double **) cgetg(d, t_VECSMALL); }
695 :
696 : static double
697 177075474 : itodbl_exp(GEN x, long *e)
698 : {
699 177075474 : pari_sp av = avma;
700 177075474 : GEN r = itor(x,DEFAULTPREC);
701 177075474 : *e = expo(r); setexpo(r,0);
702 177075474 : return gc_double(av, rtodbl(r));
703 : }
704 :
705 : static double
706 129209745 : dbldotproduct(double *x, double *y, long n)
707 : {
708 : long i;
709 129209745 : double sum = del(x,1) * del(y,1);
710 1671961631 : for (i=2; i<=n; i++) sum += del(x,i) * del(y,i);
711 129209745 : return sum;
712 : }
713 :
714 : static double
715 2484954 : dbldotsquare(double *x, long n)
716 : {
717 : long i;
718 2484954 : double sum = del(x,1) * del(x,1);
719 8250590 : for (i=2; i<=n; i++) sum += del(x,i) * del(x,i);
720 2484954 : return sum;
721 : }
722 :
723 : static long
724 25546661 : set_line(double *appv, GEN v, long n)
725 : {
726 25546661 : long i, maxexp = 0;
727 25546661 : pari_sp av = avma;
728 25546661 : GEN e = cgetg(n+1, t_VECSMALL);
729 202622135 : for (i = 1; i <= n; i++)
730 : {
731 177075474 : del(appv,i) = itodbl_exp(gel(v,i), e+i);
732 177075474 : if (e[i] > maxexp) maxexp = e[i];
733 : }
734 202622135 : for (i = 1; i <= n; i++) del(appv,i) = ldexp(del(appv,i), e[i]-maxexp);
735 25546661 : set_avma(av); return maxexp;
736 : }
737 :
738 : static void
739 35878554 : dblrotate(double **A, long k2, long k)
740 : {
741 : long i;
742 35878554 : double *B = del(A,k2);
743 115224945 : for (i = k2; i > k; i--) del(A,i) = del(A,i-1);
744 35878554 : del(A,k) = B;
745 35878554 : }
746 : /* update G[kappa][i] from appB */
747 : static void
748 23266896 : setG_fast(double **appB, long n, double **G, long kappa, long a, long b)
749 : { long i;
750 109217086 : for (i = a; i <= b; i++)
751 85950190 : dmael(G,kappa,i) = dbldotproduct(del(appB,kappa), del(appB,i), n);
752 23266896 : }
753 : /* update G[i][kappa] from appB */
754 : static void
755 17613082 : setG2_fast(double **appB, long n, double **G, long kappa, long a, long b)
756 : { long i;
757 60872637 : for (i = a; i <= b; i++)
758 43259555 : dmael(G,i,kappa) = dbldotproduct(del(appB,kappa), del(appB,i), n);
759 17613082 : }
760 : const long EX0 = -2; /* uninitialized; any value less than expo(0.51) = -1 */
761 :
762 : #ifdef LONG_IS_64BIT
763 : typedef long s64;
764 : #define addmuliu64_inplace addmuliu_inplace
765 : #define submuliu64_inplace submuliu_inplace
766 : #define submuliu642n submuliu2n
767 : #define addmuliu642n addmuliu2n
768 : #else
769 : typedef long long s64;
770 : typedef unsigned long long u64;
771 :
772 : INLINE GEN
773 22027403 : u64toi(u64 x)
774 : {
775 : GEN y;
776 : ulong h;
777 22027403 : if (!x) return gen_0;
778 22027403 : h = x>>32;
779 22027403 : if (!h) return utoipos(x);
780 1277225 : y = cgetipos(4);
781 1277225 : *int_LSW(y) = x&0xFFFFFFFF;
782 1277225 : *int_MSW(y) = x>>32;
783 1277225 : return y;
784 : }
785 :
786 : INLINE GEN
787 729530 : u64toineg(u64 x)
788 : {
789 : GEN y;
790 : ulong h;
791 729530 : if (!x) return gen_0;
792 729530 : h = x>>32;
793 729530 : if (!h) return utoineg(x);
794 729530 : y = cgetineg(4);
795 729530 : *int_LSW(y) = x&0xFFFFFFFF;
796 729530 : *int_MSW(y) = x>>32;
797 729530 : return y;
798 : }
799 : INLINE GEN
800 10613593 : addmuliu64_inplace(GEN x, GEN y, u64 u) { return addmulii(x, y, u64toi(u)); }
801 :
802 : INLINE GEN
803 10666126 : submuliu64_inplace(GEN x, GEN y, u64 u) { return submulii(x, y, u64toi(u)); }
804 :
805 : INLINE GEN
806 729530 : addmuliu642n(GEN x, GEN y, u64 u, long e) { return submulshift(x, y, u64toineg(u), e); }
807 :
808 : INLINE GEN
809 747684 : submuliu642n(GEN x, GEN y, u64 u, long e) { return submulshift(x, y, u64toi(u), e); }
810 :
811 : #endif
812 :
813 : /* Babai's Nearest Plane algorithm (iterative); see Babai() */
814 : static int
815 31678254 : Babai_fast(pari_sp av, long kappa, GEN *pB, GEN *pU, double **mu, double **r,
816 : double *s, double **appB, GEN expoB, double **G,
817 : long a, long zeros, long maxG, double eta)
818 : {
819 31678254 : GEN B = *pB, U = *pU;
820 31678254 : const long n = nbrows(B), d = U ? lg(U)-1: 0;
821 31678254 : long k, aa = (a > zeros)? a : zeros+1;
822 31678254 : long emaxmu = EX0, emax2mu = EX0;
823 : s64 xx;
824 31678254 : int did_something = 0;
825 : /* N.B: we set d = 0 (resp. n = 0) to avoid updating U (resp. B) */
826 :
827 17823246 : for (;;) {
828 49501500 : int go_on = 0;
829 49501500 : long i, j, emax3mu = emax2mu;
830 :
831 49501500 : if (gc_needed(av,2))
832 : {
833 235 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[1], a=%ld", aa);
834 235 : gc_lll(av,2,&B,&U);
835 : }
836 : /* Step2: compute the GSO for stage kappa */
837 49501500 : emax2mu = emaxmu; emaxmu = EX0;
838 196796624 : for (j=aa; j<kappa; j++)
839 : {
840 147295124 : double g = dmael(G,kappa,j);
841 682323156 : for (k = zeros+1; k < j; k++) g -= dmael(mu,j,k) * dmael(r,kappa,k);
842 147295124 : dmael(r,kappa,j) = g;
843 147295124 : dmael(mu,kappa,j) = dmael(r,kappa,j) / dmael(r,j,j);
844 147295124 : emaxmu = maxss(emaxmu, expoB[kappa]-expoB[j]);
845 : }
846 : /* maxmu doesn't decrease fast enough */
847 49501500 : if (emax3mu != EX0 && emax3mu <= emax2mu + 5) {*pB = B; *pU = U; return 1;}
848 :
849 186322828 : for (j=kappa-1; j>zeros; j--)
850 : {
851 154649280 : double tmp = fabs(ldexp (dmael(mu,kappa,j), expoB[kappa]-expoB[j]));
852 154649280 : if (tmp>eta) { go_on = 1; break; }
853 : }
854 :
855 : /* Step3--5: compute the X_j's */
856 49496794 : if (go_on)
857 85135798 : for (j=kappa-1; j>zeros; j--)
858 : { /* The code below seemingly handles U = NULL, but in this case d = 0 */
859 67312552 : int e = expoB[j] - expoB[kappa];
860 67312552 : double tmp = ldexp(dmael(mu,kappa,j), -e), atmp = fabs(tmp);
861 : /* tmp = Inf is allowed */
862 67312552 : if (atmp <= .5) continue; /* size-reduced */
863 36940453 : if (gc_needed(av,2))
864 : {
865 473 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[2], a=%ld, j=%ld", aa,j);
866 473 : gc_lll(av,2,&B,&U);
867 : }
868 36940453 : did_something = 1;
869 : /* we consider separately the case |X| = 1 */
870 36940453 : if (atmp <= 1.5)
871 : {
872 25471874 : if (dmael(mu,kappa,j) > 0) { /* in this case, X = 1 */
873 55924508 : for (k=zeros+1; k<j; k++)
874 42957757 : dmael(mu,kappa,k) -= ldexp(dmael(mu,j,k), e);
875 192277002 : for (i=1; i<=n; i++)
876 179310251 : gmael(B,kappa,i) = subii(gmael(B,kappa,i), gmael(B,j,i));
877 137843034 : for (i=1; i<=d; i++)
878 124876283 : gmael(U,kappa,i) = subii(gmael(U,kappa,i), gmael(U,j,i));
879 : } else { /* otherwise X = -1 */
880 55105826 : for (k=zeros+1; k<j; k++)
881 42600703 : dmael(mu,kappa,k) += ldexp(dmael(mu,j,k), e);
882 189622434 : for (i=1; i<=n; i++)
883 177117311 : gmael(B,kappa,i) = addii(gmael(B,kappa,i), gmael(B,j,i));
884 135098515 : for (i=1; i<=d; i++)
885 122593392 : gmael(U,kappa,i) = addii(gmael(U,kappa,i), gmael(U,j,i));
886 : }
887 25471874 : continue;
888 : }
889 : /* we have |X| >= 2 */
890 11468579 : if (atmp < 9007199254740992.)
891 : {
892 10608188 : tmp = pari_rint(tmp);
893 26390058 : for (k=zeros+1; k<j; k++)
894 15781870 : dmael(mu,kappa,k) -= ldexp(tmp * dmael(mu,j,k), e);
895 10608188 : xx = (s64) tmp;
896 10608188 : if (xx > 0) /* = xx */
897 : {
898 50244907 : for (i=1; i<=n; i++)
899 44910709 : gmael(B,kappa,i) = submuliu64_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
900 36930013 : for (i=1; i<=d; i++)
901 31595815 : gmael(U,kappa,i) = submuliu64_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
902 : }
903 : else /* = -xx */
904 : {
905 49948300 : for (i=1; i<=n; i++)
906 44674310 : gmael(B,kappa,i) = addmuliu64_inplace(gmael(B,kappa,i), gmael(B,j,i), -xx);
907 36560249 : for (i=1; i<=d; i++)
908 31286259 : gmael(U,kappa,i) = addmuliu64_inplace(gmael(U,kappa,i), gmael(U,j,i), -xx);
909 : }
910 : }
911 : else
912 : {
913 : int E;
914 860391 : xx = (s64) ldexp(frexp(dmael(mu,kappa,j), &E), 53);
915 860391 : E -= e + 53;
916 860391 : if (E <= 0)
917 : {
918 0 : xx = xx << -E;
919 0 : for (k=zeros+1; k<j; k++)
920 0 : dmael(mu,kappa,k) -= ldexp(((double)xx) * dmael(mu,j,k), e);
921 0 : if (xx > 0) /* = xx */
922 : {
923 0 : for (i=1; i<=n; i++)
924 0 : gmael(B,kappa,i) = submuliu64_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
925 0 : for (i=1; i<=d; i++)
926 0 : gmael(U,kappa,i) = submuliu64_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
927 : }
928 : else /* = -xx */
929 : {
930 0 : for (i=1; i<=n; i++)
931 0 : gmael(B,kappa,i) = addmuliu64_inplace(gmael(B,kappa,i), gmael(B,j,i), -xx);
932 0 : for (i=1; i<=d; i++)
933 0 : gmael(U,kappa,i) = addmuliu64_inplace(gmael(U,kappa,i), gmael(U,j,i), -xx);
934 : }
935 : } else
936 : {
937 2953572 : for (k=zeros+1; k<j; k++)
938 2093181 : dmael(mu,kappa,k) -= ldexp(((double)xx) * dmael(mu,j,k), E + e);
939 860391 : if (xx > 0) /* = xx */
940 : {
941 4186422 : for (i=1; i<=n; i++)
942 3753652 : gmael(B,kappa,i) = submuliu642n(gmael(B,kappa,i), gmael(B,j,i), xx, E);
943 1622207 : for (i=1; i<=d; i++)
944 1189437 : gmael(U,kappa,i) = submuliu642n(gmael(U,kappa,i), gmael(U,j,i), xx, E);
945 : }
946 : else /* = -xx */
947 : {
948 4142983 : for (i=1; i<=n; i++)
949 3715362 : gmael(B,kappa,i) = addmuliu642n(gmael(B,kappa,i), gmael(B,j,i), -xx, E);
950 1606884 : for (i=1; i<=d; i++)
951 1179263 : gmael(U,kappa,i) = addmuliu642n(gmael(U,kappa,i), gmael(U,j,i), -xx, E);
952 : }
953 : }
954 : }
955 : }
956 49496794 : if (!go_on) break; /* Anything happened? */
957 17823246 : expoB[kappa] = set_line(del(appB,kappa), gel(B,kappa), n);
958 17823246 : setG_fast(appB, n, G, kappa, zeros+1, kappa-1);
959 17823246 : aa = zeros+1;
960 : }
961 31673548 : if (did_something) setG2_fast(appB, n, G, kappa, kappa, maxG);
962 :
963 31673548 : del(s,zeros+1) = dmael(G,kappa,kappa);
964 : /* the last s[kappa-1]=r[kappa][kappa] is computed only if kappa increases */
965 124325468 : for (k=zeros+1; k<=kappa-2; k++)
966 92651920 : del(s,k+1) = del(s,k) - dmael(mu,kappa,k)*dmael(r,kappa,k);
967 31673548 : *pB = B; *pU = U; return 0;
968 : }
969 :
970 : static void
971 12439465 : update_alpha(GEN alpha, long kappa, long kappa2, long kappamax)
972 : {
973 : long i;
974 40112597 : for (i = kappa; i < kappa2; i++)
975 27673132 : if (kappa <= alpha[i]) alpha[i] = kappa;
976 40112597 : for (i = kappa2; i > kappa; i--) alpha[i] = alpha[i-1];
977 26387437 : for (i = kappa2+1; i <= kappamax; i++)
978 13947972 : if (kappa < alpha[i]) alpha[i] = kappa;
979 12439465 : alpha[kappa] = kappa;
980 12439465 : }
981 : static void
982 479947 : rotateG(GEN G, long kappa2, long kappa, long maxG, GEN Gtmp)
983 : {
984 : long i, j;
985 3863536 : for (i=1; i<=kappa2; i++) gel(Gtmp,i) = gmael(G,kappa2,i);
986 1935078 : for ( ; i<=maxG; i++) gel(Gtmp,i) = gmael(G,i,kappa2);
987 1704282 : for (i=kappa2; i>kappa; i--)
988 : {
989 6071968 : for (j=1; j<kappa; j++) gmael(G,i,j) = gmael(G,i-1,j);
990 1224335 : gmael(G,i,kappa) = gel(Gtmp,i-1);
991 4529719 : for (j=kappa+1; j<=i; j++) gmael(G,i,j) = gmael(G,i-1,j-1);
992 5084479 : for (j=kappa2+1; j<=maxG; j++) gmael(G,j,i) = gmael(G,j,i-1);
993 : }
994 2159254 : for (i=1; i<kappa; i++) gmael(G,kappa,i) = gel(Gtmp,i);
995 479947 : gmael(G,kappa,kappa) = gel(Gtmp,kappa2);
996 1935078 : for (i=kappa2+1; i<=maxG; i++) gmael(G,i,kappa) = gel(Gtmp,i);
997 479947 : }
998 : static void
999 11959518 : rotateG_fast(double **G, long kappa2, long kappa, long maxG, double *Gtmp)
1000 : {
1001 : long i, j;
1002 72530350 : for (i=1; i<=kappa2; i++) del(Gtmp,i) = dmael(G,kappa2,i);
1003 25314761 : for ( ; i<=maxG; i++) del(Gtmp,i) = dmael(G,i,kappa2);
1004 38408315 : for (i=kappa2; i>kappa; i--)
1005 : {
1006 79371538 : for (j=1; j<kappa; j++) dmael(G,i,j) = dmael(G,i-1,j);
1007 26448797 : dmael(G,i,kappa) = del(Gtmp,i-1);
1008 92853919 : for (j=kappa+1; j<=i; j++) dmael(G,i,j) = dmael(G,i-1,j-1);
1009 54999334 : for (j=kappa2+1; j<=maxG; j++) dmael(G,j,i) = dmael(G,j,i-1);
1010 : }
1011 34122035 : for (i=1; i<kappa; i++) dmael(G,kappa,i) = del(Gtmp,i);
1012 11959518 : dmael(G,kappa,kappa) = del(Gtmp,kappa2);
1013 25314761 : for (i=kappa2+1; i<=maxG; i++) dmael(G,i,kappa) = del(Gtmp,i);
1014 11959518 : }
1015 :
1016 : /* LLL-reduces (B,U) in place [apply base change transforms to B and U].
1017 : * Gram matrix, and GSO performed on matrices of 'double'.
1018 : * If (keepfirst), never swap with first vector.
1019 : * Return -1 on failure, else zeros = dim Kernel (>= 0) */
1020 : static long
1021 2107424 : fplll_fast(GEN *pB, GEN *pU, double delta, double eta, long keepfirst)
1022 : {
1023 : pari_sp av;
1024 : long kappa, kappa2, d, n, i, j, zeros, kappamax, maxG;
1025 : double **mu, **r, *s, tmp, *Gtmp, **G, **appB;
1026 2107424 : GEN alpha, expoB, B = *pB, U;
1027 2107424 : long cnt = 0;
1028 :
1029 2107424 : d = lg(B)-1;
1030 2107424 : n = nbrows(B);
1031 2107424 : U = *pU; /* NULL if inplace */
1032 :
1033 2107424 : G = cget_dblmat(d+1);
1034 2107424 : appB = cget_dblmat(d+1);
1035 2107424 : mu = cget_dblmat(d+1);
1036 2107424 : r = cget_dblmat(d+1);
1037 2107424 : s = cget_dblvec(d+1);
1038 9830839 : for (j = 1; j <= d; j++)
1039 : {
1040 7723415 : del(mu,j) = cget_dblvec(d+1);
1041 7723415 : del(r,j) = cget_dblvec(d+1);
1042 7723415 : del(appB,j) = cget_dblvec(n+1);
1043 7723415 : del(G,j) = cget_dblvec(d+1);
1044 47986034 : for (i=1; i<=d; i++) dmael(G,j,i) = 0.;
1045 : }
1046 2107424 : expoB = cgetg(d+1, t_VECSMALL);
1047 9830839 : for (i=1; i<=d; i++) expoB[i] = set_line(del(appB,i), gel(B,i), n);
1048 2107424 : Gtmp = cget_dblvec(d+1);
1049 2107424 : alpha = cgetg(d+1, t_VECSMALL);
1050 2107424 : av = avma;
1051 :
1052 : /* Step2: Initializing the main loop */
1053 2107424 : kappamax = 1;
1054 2107424 : i = 1;
1055 2107424 : maxG = d; /* later updated to kappamax */
1056 :
1057 : do {
1058 2272799 : dmael(G,i,i) = dbldotsquare(del(appB,i),n);
1059 2272799 : } while (dmael(G,i,i) <= 0 && (++i <=d));
1060 2107424 : zeros = i-1; /* all vectors B[i] with i <= zeros are zero vectors */
1061 2107424 : kappa = i;
1062 2107424 : if (zeros < d) dmael(r,zeros+1,zeros+1) = dmael(G,zeros+1,zeros+1);
1063 9665457 : for (i=zeros+1; i<=d; i++) alpha[i]=1;
1064 33780972 : while (++kappa <= d)
1065 : {
1066 31678254 : if (kappa > kappamax)
1067 : {
1068 5443650 : if (DEBUGLEVEL>=4) err_printf("K%ld ",kappa);
1069 5443650 : maxG = kappamax = kappa;
1070 5443650 : setG_fast(appB, n, G, kappa, zeros+1, kappa);
1071 : }
1072 : /* Step3: Call to the Babai algorithm, mu,r,s updated in place */
1073 31678254 : if (Babai_fast(av, kappa, &B,&U, mu,r,s, appB, expoB, G, alpha[kappa],
1074 4706 : zeros, maxG, eta)) { *pB=B; *pU=U; return -1; }
1075 :
1076 31673548 : tmp = ldexp(r[kappa-1][kappa-1] * delta, 2*(expoB[kappa-1]-expoB[kappa]));
1077 31673548 : if ((keepfirst && kappa == 2) || tmp <= del(s,kappa-1))
1078 : { /* Step4: Success of Lovasz's condition */
1079 19714030 : alpha[kappa] = kappa;
1080 19714030 : tmp = dmael(mu,kappa,kappa-1) * dmael(r,kappa,kappa-1);
1081 19714030 : dmael(r,kappa,kappa) = del(s,kappa-1)- tmp;
1082 19714030 : continue;
1083 : }
1084 : /* Step5: Find the right insertion index kappa, kappa2 = initial kappa */
1085 11959518 : if (DEBUGLEVEL>=4 && kappa==kappamax && del(s,kappa-1)!=0)
1086 0 : if (++cnt > 20) { cnt = 0; err_printf("(%ld) ", 2*expoB[1] + dblexpo(del(s,1))); }
1087 11959518 : kappa2 = kappa;
1088 : do {
1089 26448797 : kappa--;
1090 26448797 : if (kappa<zeros+2 + (keepfirst ? 1: 0)) break;
1091 19828605 : tmp = dmael(r,kappa-1,kappa-1) * delta;
1092 19828605 : tmp = ldexp(tmp, 2*(expoB[kappa-1]-expoB[kappa2]));
1093 19828605 : } while (del(s,kappa-1) <= tmp);
1094 11959518 : update_alpha(alpha, kappa, kappa2, kappamax);
1095 :
1096 : /* Step6: Update the mu's and r's */
1097 11959518 : dblrotate(mu,kappa2,kappa);
1098 11959518 : dblrotate(r,kappa2,kappa);
1099 11959518 : dmael(r,kappa,kappa) = del(s,kappa);
1100 :
1101 : /* Step7: Update B, appB, U, G */
1102 11959518 : rotate(B,kappa2,kappa);
1103 11959518 : dblrotate(appB,kappa2,kappa);
1104 11959518 : if (U) rotate(U,kappa2,kappa);
1105 11959518 : rotate(expoB,kappa2,kappa);
1106 11959518 : rotateG_fast(G,kappa2,kappa, maxG, Gtmp);
1107 :
1108 : /* Step8: Prepare the next loop iteration */
1109 11959518 : if (kappa == zeros+1 && dmael(G,kappa,kappa)<= 0)
1110 : {
1111 212155 : zeros++; kappa++;
1112 212155 : dmael(G,kappa,kappa) = dbldotsquare(del(appB,kappa),n);
1113 212155 : dmael(r,kappa,kappa) = dmael(G,kappa,kappa);
1114 : }
1115 : }
1116 2102718 : *pB = B; *pU = U; return zeros;
1117 : }
1118 :
1119 : /***************** HEURISTIC version (reduced precision) ****************/
1120 : static GEN
1121 207826 : realsqrdotproduct(GEN x)
1122 : {
1123 207826 : long i, l = lg(x);
1124 207826 : GEN z = sqrr(gel(x,1));
1125 1463083 : for (i=2; i<l; i++) z = addrr(z, sqrr(gel(x,i)));
1126 207826 : return z;
1127 : }
1128 : /* x, y non-empty vector of t_REALs, same length */
1129 : static GEN
1130 1289785 : realdotproduct(GEN x, GEN y)
1131 : {
1132 : long i, l;
1133 : GEN z;
1134 1289785 : if (x == y) return realsqrdotproduct(x);
1135 1081959 : l = lg(x); z = mulrr(gel(x,1),gel(y,1));
1136 10682294 : for (i=2; i<l; i++) z = addrr(z, mulrr(gel(x,i), gel(y,i)));
1137 1081959 : return z;
1138 : }
1139 : static void
1140 217970 : setG_heuristic(GEN appB, GEN G, long kappa, long a, long b)
1141 217970 : { pari_sp av = avma;
1142 : long i;
1143 1030069 : for (i = a; i <= b; i++)
1144 812099 : affrr(realdotproduct(gel(appB,kappa),gel(appB,i)), gmael(G,kappa,i));
1145 217970 : set_avma(av);
1146 217970 : }
1147 : static void
1148 195115 : setG2_heuristic(GEN appB, GEN G, long kappa, long a, long b)
1149 195115 : { pari_sp av = avma;
1150 : long i;
1151 672801 : for (i = a; i <= b; i++)
1152 477686 : affrr(realdotproduct(gel(appB,kappa),gel(appB,i)), gmael(G,i,kappa));
1153 195115 : set_avma(av);
1154 195115 : }
1155 :
1156 : /* approximate t_REAL x as m * 2^e, where |m| < 2^bit */
1157 : static GEN
1158 24469 : truncexpo(GEN x, long bit, long *e)
1159 : {
1160 24469 : *e = expo(x) + 1 - bit;
1161 24469 : if (*e >= 0) return mantissa2nr(x, 0);
1162 1266 : *e = 0; return roundr_safe(x);
1163 : }
1164 : /* Babai's Nearest Plane algorithm (iterative); see Babai() */
1165 : static int
1166 302144 : Babai_heuristic(pari_sp av, long kappa, GEN *pB, GEN *pU, GEN mu, GEN r, GEN s,
1167 : GEN appB, GEN G, long a, long zeros, long maxG,
1168 : GEN eta, long prec)
1169 : {
1170 302144 : GEN B = *pB, U = *pU;
1171 302144 : const long n = nbrows(B), d = U ? lg(U)-1: 0, bit = prec2nbits(prec);
1172 302144 : long k, aa = (a > zeros)? a : zeros+1;
1173 302144 : int did_something = 0;
1174 302144 : long emaxmu = EX0, emax2mu = EX0;
1175 : /* N.B: we set d = 0 (resp. n = 0) to avoid updating U (resp. B) */
1176 :
1177 205259 : for (;;) {
1178 507403 : int go_on = 0;
1179 507403 : long i, j, emax3mu = emax2mu;
1180 :
1181 507403 : if (gc_needed(av,2))
1182 : {
1183 36 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[1], a=%ld", aa);
1184 36 : gc_lll(av,2,&B,&U);
1185 : }
1186 : /* Step2: compute the GSO for stage kappa */
1187 507403 : emax2mu = emaxmu; emaxmu = EX0;
1188 1993648 : for (j=aa; j<kappa; j++)
1189 : {
1190 1486245 : pari_sp btop = avma;
1191 1486245 : GEN g = gmael(G,kappa,j);
1192 5028204 : for (k = zeros+1; k<j; k++)
1193 3541959 : g = subrr(g, mulrr(gmael(mu,j,k), gmael(r,kappa,k)));
1194 1486245 : affrr(g, gmael(r,kappa,j));
1195 1486245 : affrr(divrr(gmael(r,kappa,j), gmael(r,j,j)), gmael(mu,kappa,j));
1196 1486245 : emaxmu = maxss(emaxmu, expo(gmael(mu,kappa,j)));
1197 1486245 : set_avma(btop);
1198 : }
1199 507403 : if (emax3mu != EX0 && emax3mu <= emax2mu + 5)
1200 1741 : { *pB = B; *pU = U; return 1; }
1201 :
1202 1745586 : for (j=kappa-1; j>zeros; j--)
1203 1445183 : if (abscmprr(gmael(mu,kappa,j), eta) > 0) { go_on = 1; break; }
1204 :
1205 : /* Step3--5: compute the X_j's */
1206 505662 : if (go_on)
1207 967148 : for (j=kappa-1; j>zeros; j--)
1208 : { /* The code below seemingly handles U = NULL, but in this case d = 0 */
1209 : pari_sp btop;
1210 761889 : GEN tmp = gmael(mu,kappa,j);
1211 761889 : if (absrsmall(tmp)) continue; /* size-reduced */
1212 :
1213 442325 : if (gc_needed(av,2))
1214 : {
1215 10 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[2], a=%ld, j=%ld", aa,j);
1216 10 : gc_lll(av,2,&B,&U);
1217 : }
1218 442325 : btop = avma; did_something = 1;
1219 : /* we consider separately the case |X| = 1 */
1220 442325 : if (absrsmall2(tmp))
1221 : {
1222 279843 : if (signe(tmp) > 0) { /* in this case, X = 1 */
1223 418132 : for (k=zeros+1; k<j; k++)
1224 279018 : affrr(subrr(gmael(mu,kappa,k), gmael(mu,j,k)), gmael(mu,kappa,k));
1225 139114 : set_avma(btop);
1226 1359820 : for (i=1; i<=n; i++)
1227 1220706 : gmael(B,kappa,i) = subii(gmael(B,kappa,i), gmael(B,j,i));
1228 855692 : for (i=1; i<=d; i++)
1229 716578 : gmael(U,kappa,i) = subii(gmael(U,kappa,i), gmael(U,j,i));
1230 : } else { /* otherwise X = -1 */
1231 426109 : for (k=zeros+1; k<j; k++)
1232 285380 : affrr(addrr(gmael(mu,kappa,k), gmael(mu,j,k)), gmael(mu,kappa,k));
1233 140729 : set_avma(btop);
1234 1384891 : for (i=1; i<=n; i++)
1235 1244162 : gmael(B,kappa,i) = addii(gmael(B,kappa,i), gmael(B,j,i));
1236 859555 : for (i=1; i<=d; i++)
1237 718826 : gmael(U,kappa,i) = addii(gmael(U,kappa,i),gmael(U,j,i));
1238 : }
1239 279843 : continue;
1240 : }
1241 : /* we have |X| >= 2 */
1242 162482 : if (expo(tmp) < BITS_IN_LONG)
1243 : {
1244 138013 : ulong xx = roundr_safe(tmp)[2]; /* X fits in an ulong */
1245 138013 : if (signe(tmp) > 0) /* = xx */
1246 : {
1247 169140 : for (k=zeros+1; k<j; k++)
1248 99711 : affrr(subrr(gmael(mu,kappa,k), mulur(xx, gmael(mu,j,k))),
1249 99711 : gmael(mu,kappa,k));
1250 69429 : set_avma(btop);
1251 564786 : for (i=1; i<=n; i++)
1252 495357 : gmael(B,kappa,i) = submuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
1253 330898 : for (i=1; i<=d; i++)
1254 261469 : gmael(U,kappa,i) = submuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
1255 : }
1256 : else /* = -xx */
1257 : {
1258 168053 : for (k=zeros+1; k<j; k++)
1259 99469 : affrr(addrr(gmael(mu,kappa,k), mulur(xx, gmael(mu,j,k))),
1260 99469 : gmael(mu,kappa,k));
1261 68584 : set_avma(btop);
1262 568434 : for (i=1; i<=n; i++)
1263 499850 : gmael(B,kappa,i) = addmuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
1264 315911 : for (i=1; i<=d; i++)
1265 247327 : gmael(U,kappa,i) = addmuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
1266 : }
1267 : }
1268 : else
1269 : {
1270 : long e;
1271 24469 : GEN X = truncexpo(tmp, bit, &e); /* tmp ~ X * 2^e */
1272 24469 : btop = avma;
1273 111005 : for (k=zeros+1; k<j; k++)
1274 : {
1275 86536 : GEN x = mulir(X, gmael(mu,j,k));
1276 86536 : if (e) shiftr_inplace(x, e);
1277 86536 : affrr(subrr(gmael(mu,kappa,k), x), gmael(mu,kappa,k));
1278 : }
1279 24469 : set_avma(btop);
1280 557869 : for (i=1; i<=n; i++)
1281 533400 : gmael(B,kappa,i) = submulshift(gmael(B,kappa,i), gmael(B,j,i), X, e);
1282 88217 : for (i=1; i<=d; i++)
1283 63748 : gmael(U,kappa,i) = submulshift(gmael(U,kappa,i), gmael(U,j,i), X, e);
1284 : }
1285 : }
1286 505662 : if (!go_on) break; /* Anything happened? */
1287 1638364 : for (i=1 ; i<=n; i++) affir(gmael(B,kappa,i), gmael(appB,kappa,i));
1288 205259 : setG_heuristic(appB, G, kappa, zeros+1, kappa-1);
1289 205259 : aa = zeros+1;
1290 : }
1291 300403 : if (did_something) setG2_heuristic(appB, G, kappa, kappa, maxG);
1292 300403 : affrr(gmael(G,kappa,kappa), gel(s,zeros+1));
1293 : /* the last s[kappa-1]=r[kappa][kappa] is computed only if kappa increases */
1294 300403 : av = avma;
1295 1103203 : for (k=zeros+1; k<=kappa-2; k++)
1296 802800 : affrr(subrr(gel(s,k), mulrr(gmael(mu,kappa,k), gmael(r,kappa,k))),
1297 802800 : gel(s,k+1));
1298 300403 : *pB = B; *pU = U; return gc_bool(av, 0);
1299 : }
1300 :
1301 : static GEN
1302 22443 : ZC_to_RC(GEN x, long prec)
1303 322633 : { pari_APPLY_type(t_COL,itor(gel(x,i),prec)) }
1304 :
1305 : static GEN
1306 4706 : ZM_to_RM(GEN x, long prec)
1307 27149 : { pari_APPLY_same(ZC_to_RC(gel(x,i),prec)) }
1308 :
1309 : /* LLL-reduces (B,U) in place [apply base change transforms to B and U].
1310 : * Gram matrix made of t_REAL at precision prec2, performe GSO at prec.
1311 : * If (keepfirst), never swap with first vector.
1312 : * Return -1 on failure, else zeros = dim Kernel (>= 0) */
1313 : static long
1314 4706 : fplll_heuristic(GEN *pB, GEN *pU, double DELTA, double ETA, long keepfirst,
1315 : long prec, long prec2)
1316 : {
1317 : pari_sp av, av2;
1318 : long kappa, kappa2, d, i, j, zeros, kappamax, maxG;
1319 4706 : GEN mu, r, s, tmp, Gtmp, alpha, G, appB, B = *pB, U;
1320 4706 : GEN delta = dbltor(DELTA), eta = dbltor(ETA);
1321 4706 : long cnt = 0;
1322 :
1323 4706 : d = lg(B)-1;
1324 4706 : U = *pU; /* NULL if inplace */
1325 :
1326 4706 : G = cgetg(d+1, t_MAT);
1327 4706 : mu = cgetg(d+1, t_MAT);
1328 4706 : r = cgetg(d+1, t_MAT);
1329 4706 : s = cgetg(d+1, t_VEC);
1330 4706 : appB = ZM_to_RM(B, prec2);
1331 27149 : for (j = 1; j <= d; j++)
1332 : {
1333 22443 : GEN M = cgetg(d+1, t_COL), R = cgetg(d+1, t_COL), S = cgetg(d+1, t_COL);
1334 22443 : long l = nbits2lg(prec), l2 = nbits2lg(prec2);
1335 22443 : gel(mu,j)= M;
1336 22443 : gel(r,j) = R;
1337 22443 : gel(G,j) = S;
1338 22443 : gel(s,j) = cgetg(l, t_REAL);
1339 261950 : for (i = 1; i <= d; i++)
1340 : {
1341 239507 : gel(R,i) = cgetg(l, t_REAL);
1342 239507 : gel(M,i) = cgetg(l, t_REAL);
1343 239507 : gel(S,i) = cgetg(l2, t_REAL);
1344 : }
1345 : }
1346 4706 : Gtmp = cgetg(d+1, t_VEC);
1347 4706 : alpha = cgetg(d+1, t_VECSMALL);
1348 4706 : av = avma;
1349 :
1350 : /* Step2: Initializing the main loop */
1351 4706 : kappamax = 1;
1352 4706 : i = 1;
1353 4706 : maxG = d; /* later updated to kappamax */
1354 :
1355 : do {
1356 4709 : affrr(RgV_dotsquare(gel(appB,i)), gmael(G,i,i));
1357 4709 : } while (signe(gmael(G,i,i)) == 0 && (++i <=d));
1358 4706 : zeros = i-1; /* all vectors B[i] with i <= zeros are zero vectors */
1359 4706 : kappa = i;
1360 4706 : if (zeros < d) affrr(gmael(G,zeros+1,zeros+1), gmael(r,zeros+1,zeros+1));
1361 27146 : for (i=zeros+1; i<=d; i++) alpha[i]=1;
1362 :
1363 305109 : while (++kappa <= d)
1364 : {
1365 302144 : if (kappa > kappamax)
1366 : {
1367 12711 : if (DEBUGLEVEL>=4) err_printf("K%ld ",kappa);
1368 12711 : maxG = kappamax = kappa;
1369 12711 : setG_heuristic(appB, G, kappa, zeros+1, kappa);
1370 : }
1371 : /* Step3: Call to the Babai algorithm, mu,r,s updated in place */
1372 302144 : if (Babai_heuristic(av, kappa, &B,&U, mu,r,s, appB, G, alpha[kappa], zeros,
1373 1741 : maxG, eta, prec)) { *pB = B; *pU = U; return -1; }
1374 300403 : av2 = avma;
1375 600698 : if ((keepfirst && kappa == 2) ||
1376 300295 : cmprr(mulrr(gmael(r,kappa-1,kappa-1), delta), gel(s,kappa-1)) <= 0)
1377 : { /* Step4: Success of Lovasz's condition */
1378 179718 : alpha[kappa] = kappa;
1379 179718 : tmp = mulrr(gmael(mu,kappa,kappa-1), gmael(r,kappa,kappa-1));
1380 179718 : affrr(subrr(gel(s,kappa-1), tmp), gmael(r,kappa,kappa));
1381 179718 : set_avma(av2); continue;
1382 : }
1383 : /* Step5: Find the right insertion index kappa, kappa2 = initial kappa */
1384 120685 : if (DEBUGLEVEL>=4 && kappa==kappamax && signe(gel(s,kappa-1)))
1385 0 : if (++cnt > 20) { cnt = 0; err_printf("(%ld) ", expo(gel(s,1))); }
1386 120685 : kappa2 = kappa;
1387 : do {
1388 289498 : kappa--;
1389 289498 : if (kappa < zeros+2 + (keepfirst ? 1: 0)) break;
1390 259513 : tmp = mulrr(gmael(r,kappa-1,kappa-1), delta);
1391 259513 : } while (cmprr(gel(s,kappa-1), tmp) <= 0 );
1392 120685 : set_avma(av2);
1393 120685 : update_alpha(alpha, kappa, kappa2, kappamax);
1394 :
1395 : /* Step6: Update the mu's and r's */
1396 120685 : rotate(mu,kappa2,kappa);
1397 120685 : rotate(r,kappa2,kappa);
1398 120685 : affrr(gel(s,kappa), gmael(r,kappa,kappa));
1399 :
1400 : /* Step7: Update B, appB, U, G */
1401 120685 : rotate(B,kappa2,kappa);
1402 120685 : rotate(appB,kappa2,kappa);
1403 120685 : if (U) rotate(U,kappa2,kappa);
1404 120685 : rotateG(G,kappa2,kappa, maxG, Gtmp);
1405 :
1406 : /* Step8: Prepare the next loop iteration */
1407 120685 : if (kappa == zeros+1 && !signe(gmael(G,kappa,kappa)))
1408 : {
1409 7 : zeros++; kappa++;
1410 7 : affrr(RgV_dotsquare(gel(appB,kappa)), gmael(G,kappa,kappa));
1411 7 : affrr(gmael(G,kappa,kappa), gmael(r,kappa,kappa));
1412 : }
1413 : }
1414 2965 : *pB=B; *pU=U; return zeros;
1415 : }
1416 :
1417 : /************************* PROVED version (t_INT) ***********************/
1418 : /* dpe inspired by dpe.h by Patrick Pelissier, Paul Zimmermann
1419 : * https://gforge.inria.fr/projects/dpe/
1420 : */
1421 :
1422 : typedef struct
1423 : {
1424 : double d; /* significand */
1425 : long e; /* exponent */
1426 : } dpe_t;
1427 :
1428 : #define Dmael(x,i,j) (&((x)[i][j]))
1429 : #define Del(x,i) (&((x)[i]))
1430 :
1431 : static void
1432 718524 : dperotate(dpe_t **A, long k2, long k)
1433 : {
1434 : long i;
1435 718524 : dpe_t *B = A[k2];
1436 2588198 : for (i = k2; i > k; i--) A[i] = A[i-1];
1437 718524 : A[k] = B;
1438 718524 : }
1439 :
1440 : static void
1441 153382620 : dpe_normalize0(dpe_t *x)
1442 : {
1443 : int e;
1444 153382620 : x->d = frexp(x->d, &e);
1445 153382620 : x->e += e;
1446 153382620 : }
1447 :
1448 : static void
1449 77494763 : dpe_normalize(dpe_t *x)
1450 : {
1451 77494763 : if (x->d == 0.0)
1452 2211105 : x->e = -LONG_MAX;
1453 : else
1454 75283658 : dpe_normalize0(x);
1455 77494763 : }
1456 :
1457 : static GEN
1458 26670 : dpetor(dpe_t *x)
1459 : {
1460 26670 : GEN r = dbltor(x->d);
1461 26670 : if (signe(r)==0) return r;
1462 26607 : setexpo(r, x->e-1);
1463 26607 : return r;
1464 : }
1465 :
1466 : static void
1467 33168250 : affdpe(dpe_t *y, dpe_t *x)
1468 : {
1469 33168250 : x->d = y->d;
1470 33168250 : x->e = y->e;
1471 33168250 : }
1472 :
1473 : static void
1474 22440119 : affidpe(GEN y, dpe_t *x)
1475 : {
1476 22440119 : pari_sp av = avma;
1477 22440119 : GEN r = itor(y, DEFAULTPREC);
1478 22440119 : x->e = expo(r)+1;
1479 22440119 : setexpo(r,-1);
1480 22440119 : x->d = rtodbl(r);
1481 22440119 : set_avma(av);
1482 22440119 : }
1483 :
1484 : static void
1485 3210504 : affdbldpe(double y, dpe_t *x)
1486 : {
1487 3210504 : x->d = (double)y;
1488 3210504 : x->e = 0;
1489 3210504 : dpe_normalize(x);
1490 3210504 : }
1491 :
1492 : static void
1493 75175245 : dpe_mulz(dpe_t *x, dpe_t *y, dpe_t *z)
1494 : {
1495 75175245 : z->d = x->d * y->d;
1496 75175245 : if (z->d == 0.0)
1497 10768251 : z->e = -LONG_MAX;
1498 : else
1499 : {
1500 64406994 : z->e = x->e + y->e;
1501 64406994 : dpe_normalize0(z);
1502 : }
1503 75175245 : }
1504 :
1505 : static void
1506 15649100 : dpe_divz(dpe_t *x, dpe_t *y, dpe_t *z)
1507 : {
1508 15649100 : z->d = x->d / y->d;
1509 15649100 : if (z->d == 0.0)
1510 1957132 : z->e = -LONG_MAX;
1511 : else
1512 : {
1513 13691968 : z->e = x->e - y->e;
1514 13691968 : dpe_normalize0(z);
1515 : }
1516 15649100 : }
1517 :
1518 : static void
1519 366745 : dpe_negz(dpe_t *y, dpe_t *x)
1520 : {
1521 366745 : x->d = - y->d;
1522 366745 : x->e = y->e;
1523 366745 : }
1524 :
1525 : static void
1526 6685227 : dpe_addz(dpe_t *y, dpe_t *z, dpe_t *x)
1527 : {
1528 6685227 : if (y->e > z->e + 53)
1529 985781 : affdpe(y, x);
1530 5699446 : else if (z->e > y->e + 53)
1531 92121 : affdpe(z, x);
1532 : else
1533 : {
1534 5607325 : long d = y->e - z->e;
1535 :
1536 5607325 : if (d >= 0)
1537 : {
1538 4541004 : x->d = y->d + ldexp(z->d, -d);
1539 4541004 : x->e = y->e;
1540 : }
1541 : else
1542 : {
1543 1066321 : x->d = z->d + ldexp(y->d, d);
1544 1066321 : x->e = z->e;
1545 : }
1546 5607325 : dpe_normalize(x);
1547 : }
1548 6685227 : }
1549 : static void
1550 76421246 : dpe_subz(dpe_t *y, dpe_t *z, dpe_t *x)
1551 : {
1552 76421246 : if (y->e > z->e + 53)
1553 16081986 : affdpe(y, x);
1554 60339260 : else if (z->e > y->e + 53)
1555 366745 : dpe_negz(z, x);
1556 : else
1557 : {
1558 59972515 : long d = y->e - z->e;
1559 :
1560 59972515 : if (d >= 0)
1561 : {
1562 55732011 : x->d = y->d - ldexp(z->d, -d);
1563 55732011 : x->e = y->e;
1564 : }
1565 : else
1566 : {
1567 4240504 : x->d = ldexp(y->d, d) - z->d;
1568 4240504 : x->e = z->e;
1569 : }
1570 59972515 : dpe_normalize(x);
1571 : }
1572 76421246 : }
1573 :
1574 : static void
1575 8704419 : dpe_muluz(dpe_t *y, ulong t, dpe_t *x)
1576 : {
1577 8704419 : x->d = y->d * (double)t;
1578 8704419 : x->e = y->e;
1579 8704419 : dpe_normalize(x);
1580 8704419 : }
1581 :
1582 : static void
1583 1345057 : dpe_addmuluz(dpe_t *y, dpe_t *z, ulong t, dpe_t *x)
1584 : {
1585 : dpe_t tmp;
1586 1345057 : dpe_muluz(z, t, &tmp);
1587 1345057 : dpe_addz(y, &tmp, x);
1588 1345057 : }
1589 :
1590 : static void
1591 1433659 : dpe_submuluz(dpe_t *y, dpe_t *z, ulong t, dpe_t *x)
1592 : {
1593 : dpe_t tmp;
1594 1433659 : dpe_muluz(z, t, &tmp);
1595 1433659 : dpe_subz(y, &tmp, x);
1596 1433659 : }
1597 :
1598 : static void
1599 69695441 : dpe_submulz(dpe_t *y, dpe_t *z, dpe_t *t, dpe_t *x)
1600 : {
1601 : dpe_t tmp;
1602 69695441 : dpe_mulz(z, t, &tmp);
1603 69695441 : dpe_subz(y, &tmp, x);
1604 69695441 : }
1605 :
1606 : static int
1607 5479804 : dpe_cmp(dpe_t *x, dpe_t *y)
1608 : {
1609 5479804 : int sx = x->d < 0. ? -1: x->d > 0.;
1610 5479804 : int sy = y->d < 0. ? -1: y->d > 0.;
1611 5479804 : int d = sx - sy;
1612 :
1613 5479804 : if (d != 0)
1614 142999 : return d;
1615 5336805 : else if (x->e > y->e)
1616 550741 : return (sx > 0) ? 1 : -1;
1617 4786064 : else if (y->e > x->e)
1618 2604063 : return (sx > 0) ? -1 : 1;
1619 : else
1620 2182001 : return (x->d < y->d) ? -1 : (x->d > y->d);
1621 : }
1622 :
1623 : static int
1624 15812775 : dpe_abscmp(dpe_t *x, dpe_t *y)
1625 : {
1626 15812775 : if (x->e > y->e)
1627 312064 : return 1;
1628 15500711 : else if (y->e > x->e)
1629 14574496 : return -1;
1630 : else
1631 926215 : return (fabs(x->d) < fabs(y->d)) ? -1 : (fabs(x->d) > fabs(y->d));
1632 : }
1633 :
1634 : static int
1635 2183624 : dpe_abssmall(dpe_t *x)
1636 : {
1637 2183624 : return (x->e <= 0) || (x->e == 1 && fabs(x->d) <= .75);
1638 : }
1639 :
1640 : static int
1641 5479804 : dpe_cmpmul(dpe_t *x, dpe_t *y, dpe_t *z)
1642 : {
1643 : dpe_t t;
1644 5479804 : dpe_mulz(x,y,&t);
1645 5479804 : return dpe_cmp(&t, z);
1646 : }
1647 :
1648 : static dpe_t *
1649 13317616 : cget_dpevec(long d)
1650 13317616 : { return (dpe_t*) stack_malloc_align(d*sizeof(dpe_t), sizeof(dpe_t)); }
1651 :
1652 : static dpe_t **
1653 3210504 : cget_dpemat(long d) { return (dpe_t **) cgetg(d, t_VECSMALL); }
1654 :
1655 : static GEN
1656 1736 : dpeM_diagonal_shallow(dpe_t **m, long d)
1657 : {
1658 : long i;
1659 1736 : GEN y = cgetg(d+1,t_VEC);
1660 28406 : for (i=1; i<=d; i++) gel(y, i) = dpetor(Dmael(m,i,i));
1661 1736 : return y;
1662 : }
1663 :
1664 : static void
1665 2183624 : affii_or_copy_gc(pari_sp av, GEN x, GEN *y)
1666 : {
1667 2183624 : long l = lg(*y);
1668 2183624 : if (lgefint(x) <= l && isonstack(*y))
1669 : {
1670 2183612 : affii(x,*y);
1671 2183612 : set_avma(av);
1672 : }
1673 : else
1674 12 : *y = gc_INT(av, x);
1675 2183624 : }
1676 :
1677 : /* *x -= u*y */
1678 : INLINE void
1679 12229399 : submulziu(GEN *x, GEN y, ulong u)
1680 : {
1681 : pari_sp av;
1682 12229399 : long ly = lgefint(y);
1683 12229399 : if (ly == 2) return;
1684 6293883 : av = avma;
1685 6293883 : (void)new_chunk(3+ly+lgefint(*x)); /* HACK */
1686 6293883 : y = mului(u,y);
1687 6293883 : set_avma(av); subzi(x, y);
1688 : }
1689 :
1690 : /* *x += u*y */
1691 : INLINE void
1692 10826608 : addmulziu(GEN *x, GEN y, ulong u)
1693 : {
1694 : pari_sp av;
1695 10826608 : long ly = lgefint(y);
1696 10826608 : if (ly == 2) return;
1697 5767632 : av = avma;
1698 5767632 : (void)new_chunk(3+ly+lgefint(*x)); /* HACK */
1699 5767632 : y = mului(u,y);
1700 5767632 : set_avma(av); addzi(x, y);
1701 : }
1702 :
1703 : /************************** PROVED version (dpe) *************************/
1704 :
1705 : /* Babai's Nearest Plane algorithm (iterative).
1706 : * Size-reduces b_kappa using mu_{i,j} and r_{i,j} for j<=i <kappa
1707 : * Update B[,kappa]; compute mu_{kappa,j}, r_{kappa,j} for j<=kappa and s[kappa]
1708 : * mu, r, s updated in place (affrr). Return 1 on failure, else 0. */
1709 : static int
1710 4767714 : Babai_dpe(pari_sp av, long kappa, GEN *pG, GEN *pB, GEN *pU, dpe_t **mu, dpe_t **r, dpe_t *s,
1711 : long a, long zeros, long maxG, dpe_t *eta)
1712 : {
1713 4767714 : GEN G = *pG, B = *pB, U = *pU, ztmp;
1714 4767714 : long k, d, n, aa = a > zeros? a: zeros+1;
1715 4767714 : long emaxmu = EX0, emax2mu = EX0;
1716 : /* N.B: we set d = 0 (resp. n = 0) to avoid updating U (resp. B) */
1717 4767714 : d = U? lg(U)-1: 0;
1718 4767714 : n = B? nbrows(B): 0;
1719 603686 : for (;;) {
1720 5371400 : int go_on = 0;
1721 5371400 : long i, j, emax3mu = emax2mu;
1722 :
1723 5371400 : if (gc_needed(av,2))
1724 : {
1725 0 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[1], a=%ld", aa);
1726 0 : gc_lll(av,3,&G,&B,&U);
1727 : }
1728 : /* Step2: compute the GSO for stage kappa */
1729 5371400 : emax2mu = emaxmu; emaxmu = EX0;
1730 21020500 : for (j=aa; j<kappa; j++)
1731 : {
1732 : dpe_t g;
1733 15649100 : affidpe(gmael(G,kappa,j), &g);
1734 71000758 : for (k = zeros+1; k < j; k++)
1735 55351658 : dpe_submulz(&g, Dmael(mu,j,k), Dmael(r,kappa,k), &g);
1736 15649100 : affdpe(&g, Dmael(r,kappa,j));
1737 15649100 : dpe_divz(Dmael(r,kappa,j), Dmael(r,j,j), Dmael(mu,kappa,j));
1738 15649100 : emaxmu = maxss(emaxmu, Dmael(mu,kappa,j)->e);
1739 : }
1740 5371400 : if (emax3mu != EX0 && emax3mu <= emax2mu + 5) /* precision too low */
1741 0 : { *pG = G; *pB = B; *pU = U; return 1; }
1742 :
1743 20580489 : for (j=kappa-1; j>zeros; j--)
1744 15812775 : if (dpe_abscmp(Dmael(mu,kappa,j), eta) > 0) { go_on = 1; break; }
1745 :
1746 : /* Step3--5: compute the X_j's */
1747 5371400 : if (go_on)
1748 4228640 : for (j=kappa-1; j>zeros; j--)
1749 : {
1750 : pari_sp btop;
1751 3624954 : dpe_t *tmp = Dmael(mu,kappa,j);
1752 3624954 : if (tmp->e < 0) continue; /* (essentially) size-reduced */
1753 :
1754 2183624 : if (gc_needed(av,2))
1755 : {
1756 0 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[2], a=%ld, j=%ld", aa,j);
1757 0 : gc_lll(av,3,&G,&B,&U);
1758 : }
1759 : /* we consider separately the case |X| = 1 */
1760 2183624 : if (dpe_abssmall(tmp))
1761 : {
1762 1135663 : if (tmp->d > 0) { /* in this case, X = 1 */
1763 2927361 : for (k=zeros+1; k<j; k++)
1764 2358410 : dpe_subz(Dmael(mu,kappa,k), Dmael(mu,j,k), Dmael(mu,kappa,k));
1765 6368504 : for (i=1; i<=n; i++)
1766 5799553 : subzi(&gmael(B,kappa,i), gmael(B,j,i));
1767 7700073 : for (i=1; i<=d; i++)
1768 7131122 : subzi(&gmael(U,kappa,i), gmael(U,j,i));
1769 568951 : btop = avma;
1770 568951 : ztmp = subii(gmael(G,j,j), shifti(gmael(G,kappa,j), 1));
1771 568951 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
1772 568951 : affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
1773 3842751 : for (i=1; i<=j; i++)
1774 3273800 : subzi(&gmael(G,kappa,i), gmael(G,j,i));
1775 3197684 : for (i=j+1; i<kappa; i++)
1776 2628733 : subzi(&gmael(G,kappa,i), gmael(G,i,j));
1777 2822218 : for (i=kappa+1; i<=maxG; i++)
1778 2253267 : subzi(&gmael(G,i,kappa), gmael(G,i,j));
1779 : } else { /* otherwise X = -1 */
1780 2914915 : for (k=zeros+1; k<j; k++)
1781 2348203 : dpe_addz(Dmael(mu,kappa,k), Dmael(mu,j,k), Dmael(mu,kappa,k));
1782 6352825 : for (i=1; i<=n; i++)
1783 5786113 : addzi(&gmael(B,kappa,i),gmael(B,j,i));
1784 7574480 : for (i=1; i<=d; i++)
1785 7007768 : addzi(&gmael(U,kappa,i),gmael(U,j,i));
1786 566712 : btop = avma;
1787 566712 : ztmp = addii(gmael(G,j,j), shifti(gmael(G,kappa,j), 1));
1788 566712 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
1789 566712 : affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
1790 3763207 : for (i=1; i<=j; i++)
1791 3196495 : addzi(&gmael(G,kappa,i), gmael(G,j,i));
1792 3184909 : for (i=j+1; i<kappa; i++)
1793 2618197 : addzi(&gmael(G,kappa,i), gmael(G,i,j));
1794 2771969 : for (i=kappa+1; i<=maxG; i++)
1795 2205257 : addzi(&gmael(G,i,kappa), gmael(G,i,j));
1796 : }
1797 1135663 : continue;
1798 : }
1799 : /* we have |X| >= 2 */
1800 1047961 : if (tmp->e < BITS_IN_LONG-1)
1801 : {
1802 623086 : if (tmp->d > 0)
1803 : {
1804 335416 : ulong xx = (ulong) pari_rint(ldexp(tmp->d, tmp->e)); /* X fits in an ulong */
1805 1769075 : for (k=zeros+1; k<j; k++)
1806 1433659 : dpe_submuluz(Dmael(mu,kappa,k), Dmael(mu,j,k), xx, Dmael(mu,kappa,k));
1807 4765145 : for (i=1; i<=n; i++)
1808 4429729 : submulziu(&gmael(B,kappa,i), gmael(B,j,i), xx);
1809 3253134 : for (i=1; i<=d; i++)
1810 2917718 : submulziu(&gmael(U,kappa,i), gmael(U,j,i), xx);
1811 335416 : btop = avma;
1812 335416 : ztmp = submuliu2n(mulii(gmael(G,j,j), sqru(xx)), gmael(G,kappa,j), xx, 1);
1813 335416 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
1814 335416 : affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
1815 2514827 : for (i=1; i<=j; i++)
1816 2179411 : submulziu(&gmael(G,kappa,i), gmael(G,j,i), xx);
1817 2039934 : for (i=j+1; i<kappa; i++)
1818 1704518 : submulziu(&gmael(G,kappa,i), gmael(G,i,j), xx);
1819 1333439 : for (i=kappa+1; i<=maxG; i++)
1820 998023 : submulziu(&gmael(G,i,kappa), gmael(G,i,j), xx);
1821 : }
1822 : else
1823 : {
1824 287670 : ulong xx = (ulong) pari_rint(ldexp(-tmp->d, tmp->e)); /* X fits in an ulong */
1825 1632727 : for (k=zeros+1; k<j; k++)
1826 1345057 : dpe_addmuluz(Dmael(mu,kappa,k), Dmael(mu,j,k), xx, Dmael(mu,kappa,k));
1827 4688090 : for (i=1; i<=n; i++)
1828 4400420 : addmulziu(&gmael(B,kappa,i), gmael(B,j,i), xx);
1829 2503532 : for (i=1; i<=d; i++)
1830 2215862 : addmulziu(&gmael(U,kappa,i), gmael(U,j,i), xx);
1831 287670 : btop = avma;
1832 287670 : ztmp = addmuliu2n(mulii(gmael(G,j,j), sqru(xx)), gmael(G,kappa,j), xx, 1);
1833 287670 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
1834 287670 : affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
1835 2168749 : for (i=1; i<=j; i++)
1836 1881079 : addmulziu(&gmael(G,kappa,i), gmael(G,j,i), xx);
1837 1893326 : for (i=j+1; i<kappa; i++)
1838 1605656 : addmulziu(&gmael(G,kappa,i), gmael(G,i,j), xx);
1839 1011261 : for (i=kappa+1; i<=maxG; i++)
1840 723591 : addmulziu(&gmael(G,i,kappa), gmael(G,i,j), xx);
1841 : }
1842 : }
1843 : else
1844 : {
1845 424875 : long e = tmp->e - BITS_IN_LONG + 1;
1846 424875 : if (tmp->d > 0)
1847 : {
1848 211202 : ulong xx = (ulong) pari_rint(ldexp(tmp->d, BITS_IN_LONG - 1));
1849 3144938 : for (k=zeros+1; k<j; k++)
1850 : {
1851 : dpe_t x;
1852 2933736 : dpe_muluz(Dmael(mu,j,k), xx, &x);
1853 2933736 : x.e += e;
1854 2933736 : dpe_subz(Dmael(mu,kappa,k), &x, Dmael(mu,kappa,k));
1855 : }
1856 10752252 : for (i=1; i<=n; i++)
1857 10541050 : submulzu2n(&gmael(B,kappa,i), gmael(B,j,i), xx, e);
1858 319409 : for (i=1; i<=d; i++)
1859 108207 : submulzu2n(&gmael(U,kappa,i), gmael(U,j,i), xx, e);
1860 211202 : btop = avma;
1861 211202 : ztmp = submuliu2n(mulshift(gmael(G,j,j), sqru(xx), 2*e),
1862 211202 : gmael(G,kappa,j), xx, e+1);
1863 211202 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
1864 211202 : affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
1865 3358154 : for (i=1; i<=j; i++)
1866 3146952 : submulzu2n(&gmael(G,kappa,i), gmael(G,j,i), xx, e);
1867 3280865 : for ( ; i<kappa; i++)
1868 3069663 : submulzu2n(&gmael(G,kappa,i), gmael(G,i,j), xx, e);
1869 213227 : for (i=kappa+1; i<=maxG; i++)
1870 2025 : submulzu2n(&gmael(G,i,kappa), gmael(G,i,j), xx, e);
1871 : } else
1872 : {
1873 213673 : ulong xx = (ulong) pari_rint(ldexp(-tmp->d, BITS_IN_LONG - 1));
1874 3205640 : for (k=zeros+1; k<j; k++)
1875 : {
1876 : dpe_t x;
1877 2991967 : dpe_muluz(Dmael(mu,j,k), xx, &x);
1878 2991967 : x.e += e;
1879 2991967 : dpe_addz(Dmael(mu,kappa,k), &x, Dmael(mu,kappa,k));
1880 : }
1881 10891725 : for (i=1; i<=n; i++)
1882 10678052 : addmulzu2n(&gmael(B,kappa,i), gmael(B,j,i), xx, e);
1883 322287 : for (i=1; i<=d; i++)
1884 108614 : addmulzu2n(&gmael(U,kappa,i), gmael(U,j,i), xx, e);
1885 213673 : btop = avma;
1886 213673 : ztmp = addmuliu2n(mulshift(gmael(G,j,j), sqru(xx), 2*e),
1887 213673 : gmael(G,kappa,j), xx, e+1);
1888 213673 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
1889 213673 : affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
1890 3421132 : for (i=1; i<=j; i++)
1891 3207459 : addmulzu2n(&gmael(G,kappa,i), gmael(G,j,i), xx, e);
1892 3288274 : for ( ; i<kappa; i++)
1893 3074601 : addmulzu2n(&gmael(G,kappa,i), gmael(G,i,j), xx, e);
1894 215539 : for (i=kappa+1; i<=maxG; i++)
1895 1866 : addmulzu2n(&gmael(G,i,kappa), gmael(G,i,j), xx, e);
1896 : }
1897 : }
1898 : }
1899 5371400 : if (!go_on) break; /* Anything happened? */
1900 603686 : aa = zeros+1;
1901 : }
1902 :
1903 4767714 : affidpe(gmael(G,kappa,kappa), Del(s,zeros+1));
1904 : /* the last s[kappa-1]=r[kappa][kappa] is computed only if kappa increases */
1905 14703045 : for (k=zeros+1; k<=kappa-2; k++)
1906 9935331 : dpe_submulz(Del(s,k), Dmael(mu,kappa,k), Dmael(r,kappa,k), Del(s,k+1));
1907 4767714 : *pG = G; *pB = B; *pU = U; return 0;
1908 : }
1909 :
1910 : /* G integral Gram matrix, LLL-reduces (G,B,U) in place [apply base change
1911 : * transforms to B and U]. If (keepfirst), never swap with first vector.
1912 : * If G = NULL, we compute the Gram matrix incrementally.
1913 : * Return -1 on failure, else zeros = dim Kernel (>= 0) */
1914 : static long
1915 1605252 : fplll_dpe(GEN *pG, GEN *pB, GEN *pU, GEN *pr, double DELTA, double ETA,
1916 : long keepfirst)
1917 : {
1918 : pari_sp av;
1919 1605252 : GEN Gtmp, alpha, G = *pG, B = *pB, U = *pU;
1920 1605252 : long d, maxG, kappa, kappa2, i, j, zeros, kappamax, incgram = !G, cnt = 0;
1921 : dpe_t delta, eta, **mu, **r, *s;
1922 1605252 : affdbldpe(DELTA,&delta);
1923 1605252 : affdbldpe(ETA,&eta);
1924 :
1925 1605252 : if (incgram)
1926 : { /* incremental Gram matrix */
1927 1544696 : maxG = 2; d = lg(B)-1;
1928 1544696 : G = zeromatcopy(d, d);
1929 : }
1930 : else
1931 60556 : maxG = d = lg(G)-1;
1932 :
1933 1605252 : mu = cget_dpemat(d+1);
1934 1605252 : r = cget_dpemat(d+1);
1935 1605252 : s = cget_dpevec(d+1);
1936 7461434 : for (j = 1; j <= d; j++)
1937 : {
1938 5856182 : mu[j]= cget_dpevec(d+1);
1939 5856182 : r[j] = cget_dpevec(d+1);
1940 : }
1941 1605252 : Gtmp = cgetg(d+1, t_VEC);
1942 1605252 : alpha = cgetg(d+1, t_VECSMALL);
1943 1605252 : av = avma;
1944 :
1945 : /* Step2: Initializing the main loop */
1946 1605252 : kappamax = 1;
1947 1605252 : i = 1;
1948 : do {
1949 1988109 : if (incgram) gmael(G,i,i) = ZV_dotsquare(gel(B,i));
1950 1988109 : affidpe(gmael(G,i,i), Dmael(r,i,i));
1951 1988109 : } while (!signe(gmael(G,i,i)) && ++i <= d);
1952 1605252 : zeros = i-1; /* all basis vectors b_i with i <= zeros are zero vectors */
1953 1605252 : kappa = i;
1954 7078570 : for (i=zeros+1; i<=d; i++) alpha[i]=1;
1955 :
1956 6372966 : while (++kappa <= d)
1957 : {
1958 4767714 : if (kappa > kappamax)
1959 : {
1960 3868073 : if (DEBUGLEVEL>=4) err_printf("K%ld ",kappa);
1961 3868073 : kappamax = kappa;
1962 3868073 : if (incgram)
1963 : {
1964 16227078 : for (i=zeros+1; i<=kappa; i++)
1965 12559526 : gmael(G,kappa,i) = ZV_dotproduct(gel(B,kappa), gel(B,i));
1966 3667552 : maxG = kappamax;
1967 : }
1968 : }
1969 : /* Step3: Call to the Babai algorithm, mu,r,s updated in place */
1970 4767714 : if (Babai_dpe(av, kappa, &G,&B,&U, mu,r,s, alpha[kappa], zeros, maxG, &eta))
1971 0 : { *pG = incgram? NULL: G; *pB = B; *pU = U; return -1; }
1972 9435211 : if ((keepfirst && kappa == 2) ||
1973 4667497 : dpe_cmpmul(Dmael(r,kappa-1,kappa-1), &delta, Del(s,kappa-1)) <= 0)
1974 : { /* Step4: Success of Lovasz's condition */
1975 4408452 : alpha[kappa] = kappa;
1976 4408452 : dpe_submulz(Del(s,kappa-1), Dmael(mu,kappa,kappa-1), Dmael(r,kappa,kappa-1), Dmael(r,kappa,kappa));
1977 4408452 : continue;
1978 : }
1979 : /* Step5: Find the right insertion index kappa, kappa2 = initial kappa */
1980 359262 : if (DEBUGLEVEL>=4 && kappa==kappamax && Del(s,kappa-1)->d)
1981 0 : if (++cnt > 20) { cnt = 0; err_printf("(%ld) ", Del(s,1)->e-1); }
1982 359262 : kappa2 = kappa;
1983 : do {
1984 934837 : kappa--;
1985 934837 : if (kappa < zeros+2 + (keepfirst ? 1: 0)) break;
1986 812307 : } while (dpe_cmpmul(Dmael(r,kappa-1,kappa-1), &delta, Del(s,kappa-1)) >= 0);
1987 359262 : update_alpha(alpha, kappa, kappa2, kappamax);
1988 :
1989 : /* Step6: Update the mu's and r's */
1990 359262 : dperotate(mu, kappa2, kappa);
1991 359262 : dperotate(r, kappa2, kappa);
1992 359262 : affdpe(Del(s,kappa), Dmael(r,kappa,kappa));
1993 :
1994 : /* Step7: Update G, B, U */
1995 359262 : if (U) rotate(U, kappa2, kappa);
1996 359262 : if (B) rotate(B, kappa2, kappa);
1997 359262 : rotateG(G,kappa2,kappa, maxG, Gtmp);
1998 :
1999 : /* Step8: Prepare the next loop iteration */
2000 359262 : if (kappa == zeros+1 && !signe(gmael(G,kappa,kappa)))
2001 : {
2002 35196 : zeros++; kappa++;
2003 35196 : affidpe(gmael(G,kappa,kappa), Dmael(r,kappa,kappa));
2004 : }
2005 : }
2006 1605252 : if (pr) *pr = dpeM_diagonal_shallow(r,d);
2007 1605252 : *pG = G; *pB = B; *pU = U; return zeros; /* success */
2008 : }
2009 :
2010 :
2011 : /************************** PROVED version (t_INT) *************************/
2012 :
2013 : /* Babai's Nearest Plane algorithm (iterative).
2014 : * Size-reduces b_kappa using mu_{i,j} and r_{i,j} for j<=i <kappa
2015 : * Update B[,kappa]; compute mu_{kappa,j}, r_{kappa,j} for j<=kappa and s[kappa]
2016 : * mu, r, s updated in place (affrr). Return 1 on failure, else 0. */
2017 : static int
2018 0 : Babai(pari_sp av, long kappa, GEN *pG, GEN *pB, GEN *pU, GEN mu, GEN r, GEN s,
2019 : long a, long zeros, long maxG, GEN eta, long prec)
2020 : {
2021 0 : GEN G = *pG, B = *pB, U = *pU, ztmp;
2022 0 : long k, aa = a > zeros? a: zeros+1;
2023 0 : const long n = B? nbrows(B): 0, d = U ? lg(U)-1: 0, bit = prec2nbits(prec);
2024 0 : long emaxmu = EX0, emax2mu = EX0;
2025 : /* N.B: we set d = 0 (resp. n = 0) to avoid updating U (resp. B) */
2026 :
2027 0 : for (;;) {
2028 0 : int go_on = 0;
2029 0 : long i, j, emax3mu = emax2mu;
2030 :
2031 0 : if (gc_needed(av,2))
2032 : {
2033 0 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[1], a=%ld", aa);
2034 0 : gc_lll(av,3,&G,&B,&U);
2035 : }
2036 : /* Step2: compute the GSO for stage kappa */
2037 0 : emax2mu = emaxmu; emaxmu = EX0;
2038 0 : for (j=aa; j<kappa; j++)
2039 : {
2040 0 : pari_sp btop = avma;
2041 0 : GEN g = gmael(G,kappa,j);
2042 0 : k = zeros + 1;
2043 0 : if (k >= j)
2044 0 : affir(g, gmael(r,kappa,j));
2045 : else
2046 : {
2047 0 : g = subir(g, mulrr(gmael(mu,j,k), gmael(r,kappa,k)));
2048 0 : for (k++; k < j; k++)
2049 0 : g = subrr(g, mulrr(gmael(mu,j,k), gmael(r,kappa,k)));
2050 0 : affrr(g, gmael(r,kappa,j));
2051 : }
2052 0 : affrr(divrr(gmael(r,kappa,j), gmael(r,j,j)), gmael(mu,kappa,j));
2053 0 : emaxmu = maxss(emaxmu, expo(gmael(mu,kappa,j)));
2054 0 : set_avma(btop);
2055 : }
2056 0 : if (emax3mu != EX0 && emax3mu <= emax2mu + 5) /* precision too low */
2057 0 : { *pG = G; *pB = B; *pU = U; return 1; }
2058 :
2059 0 : for (j=kappa-1; j>zeros; j--)
2060 0 : if (abscmprr(gmael(mu,kappa,j), eta) > 0) { go_on = 1; break; }
2061 :
2062 : /* Step3--5: compute the X_j's */
2063 0 : if (go_on)
2064 0 : for (j=kappa-1; j>zeros; j--)
2065 : {
2066 : pari_sp btop;
2067 0 : GEN tmp = gmael(mu,kappa,j);
2068 0 : if (absrsmall(tmp)) continue; /* size-reduced */
2069 :
2070 0 : if (gc_needed(av,2))
2071 : {
2072 0 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[2], a=%ld, j=%ld", aa,j);
2073 0 : gc_lll(av,3,&G,&B,&U);
2074 : }
2075 0 : btop = avma;
2076 : /* we consider separately the case |X| = 1 */
2077 0 : if (absrsmall2(tmp))
2078 : {
2079 0 : if (signe(tmp) > 0) { /* in this case, X = 1 */
2080 0 : for (k=zeros+1; k<j; k++)
2081 0 : affrr(subrr(gmael(mu,kappa,k), gmael(mu,j,k)), gmael(mu,kappa,k));
2082 0 : set_avma(btop);
2083 0 : for (i=1; i<=n; i++)
2084 0 : gmael(B,kappa,i) = subii(gmael(B,kappa,i), gmael(B,j,i));
2085 0 : for (i=1; i<=d; i++)
2086 0 : gmael(U,kappa,i) = subii(gmael(U,kappa,i), gmael(U,j,i));
2087 0 : btop = avma;
2088 0 : ztmp = subii(gmael(G,j,j), shifti(gmael(G,kappa,j), 1));
2089 0 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
2090 0 : gmael(G,kappa,kappa) = gc_INT(btop, ztmp);
2091 0 : for (i=1; i<=j; i++)
2092 0 : gmael(G,kappa,i) = subii(gmael(G,kappa,i), gmael(G,j,i));
2093 0 : for (i=j+1; i<kappa; i++)
2094 0 : gmael(G,kappa,i) = subii(gmael(G,kappa,i), gmael(G,i,j));
2095 0 : for (i=kappa+1; i<=maxG; i++)
2096 0 : gmael(G,i,kappa) = subii(gmael(G,i,kappa), gmael(G,i,j));
2097 : } else { /* otherwise X = -1 */
2098 0 : for (k=zeros+1; k<j; k++)
2099 0 : affrr(addrr(gmael(mu,kappa,k), gmael(mu,j,k)), gmael(mu,kappa,k));
2100 0 : set_avma(btop);
2101 0 : for (i=1; i<=n; i++)
2102 0 : gmael(B,kappa,i) = addii(gmael(B,kappa,i),gmael(B,j,i));
2103 0 : for (i=1; i<=d; i++)
2104 0 : gmael(U,kappa,i) = addii(gmael(U,kappa,i),gmael(U,j,i));
2105 0 : btop = avma;
2106 0 : ztmp = addii(gmael(G,j,j), shifti(gmael(G,kappa,j), 1));
2107 0 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
2108 0 : gmael(G,kappa,kappa) = gc_INT(btop, ztmp);
2109 0 : for (i=1; i<=j; i++)
2110 0 : gmael(G,kappa,i) = addii(gmael(G,kappa,i), gmael(G,j,i));
2111 0 : for (i=j+1; i<kappa; i++)
2112 0 : gmael(G,kappa,i) = addii(gmael(G,kappa,i), gmael(G,i,j));
2113 0 : for (i=kappa+1; i<=maxG; i++)
2114 0 : gmael(G,i,kappa) = addii(gmael(G,i,kappa), gmael(G,i,j));
2115 : }
2116 0 : continue;
2117 : }
2118 : /* we have |X| >= 2 */
2119 0 : if (expo(tmp) < BITS_IN_LONG)
2120 : {
2121 0 : ulong xx = roundr_safe(tmp)[2]; /* X fits in an ulong */
2122 0 : if (signe(tmp) > 0) /* = xx */
2123 : {
2124 0 : for (k=zeros+1; k<j; k++)
2125 0 : affrr(subrr(gmael(mu,kappa,k), mulur(xx, gmael(mu,j,k))),
2126 0 : gmael(mu,kappa,k));
2127 0 : set_avma(btop);
2128 0 : for (i=1; i<=n; i++)
2129 0 : gmael(B,kappa,i) = submuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
2130 0 : for (i=1; i<=d; i++)
2131 0 : gmael(U,kappa,i) = submuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
2132 0 : btop = avma;
2133 0 : ztmp = submuliu2n(mulii(gmael(G,j,j), sqru(xx)), gmael(G,kappa,j), xx, 1);
2134 0 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
2135 0 : gmael(G,kappa,kappa) = gc_INT(btop, ztmp);
2136 0 : for (i=1; i<=j; i++)
2137 0 : gmael(G,kappa,i) = submuliu_inplace(gmael(G,kappa,i), gmael(G,j,i), xx);
2138 0 : for (i=j+1; i<kappa; i++)
2139 0 : gmael(G,kappa,i) = submuliu_inplace(gmael(G,kappa,i), gmael(G,i,j), xx);
2140 0 : for (i=kappa+1; i<=maxG; i++)
2141 0 : gmael(G,i,kappa) = submuliu_inplace(gmael(G,i,kappa), gmael(G,i,j), xx);
2142 : }
2143 : else /* = -xx */
2144 : {
2145 0 : for (k=zeros+1; k<j; k++)
2146 0 : affrr(addrr(gmael(mu,kappa,k), mulur(xx, gmael(mu,j,k))),
2147 0 : gmael(mu,kappa,k));
2148 0 : set_avma(btop);
2149 0 : for (i=1; i<=n; i++)
2150 0 : gmael(B,kappa,i) = addmuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
2151 0 : for (i=1; i<=d; i++)
2152 0 : gmael(U,kappa,i) = addmuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
2153 0 : btop = avma;
2154 0 : ztmp = addmuliu2n(mulii(gmael(G,j,j), sqru(xx)), gmael(G,kappa,j), xx, 1);
2155 0 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
2156 0 : gmael(G,kappa,kappa) = gc_INT(btop, ztmp);
2157 0 : for (i=1; i<=j; i++)
2158 0 : gmael(G,kappa,i) = addmuliu_inplace(gmael(G,kappa,i), gmael(G,j,i), xx);
2159 0 : for (i=j+1; i<kappa; i++)
2160 0 : gmael(G,kappa,i) = addmuliu_inplace(gmael(G,kappa,i), gmael(G,i,j), xx);
2161 0 : for (i=kappa+1; i<=maxG; i++)
2162 0 : gmael(G,i,kappa) = addmuliu_inplace(gmael(G,i,kappa), gmael(G,i,j), xx);
2163 : }
2164 : }
2165 : else
2166 : {
2167 : long e;
2168 0 : GEN X = truncexpo(tmp, bit, &e); /* tmp ~ X * 2^e */
2169 0 : btop = avma;
2170 0 : for (k=zeros+1; k<j; k++)
2171 : {
2172 0 : GEN x = mulir(X, gmael(mu,j,k));
2173 0 : if (e) shiftr_inplace(x, e);
2174 0 : affrr(subrr(gmael(mu,kappa,k), x), gmael(mu,kappa,k));
2175 : }
2176 0 : set_avma(btop);
2177 0 : for (i=1; i<=n; i++)
2178 0 : gmael(B,kappa,i) = submulshift(gmael(B,kappa,i), gmael(B,j,i), X, e);
2179 0 : for (i=1; i<=d; i++)
2180 0 : gmael(U,kappa,i) = submulshift(gmael(U,kappa,i), gmael(U,j,i), X, e);
2181 0 : btop = avma;
2182 0 : ztmp = submulshift(mulshift(gmael(G,j,j), sqri(X), 2*e),
2183 0 : gmael(G,kappa,j), X, e+1);
2184 0 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
2185 0 : gmael(G,kappa,kappa) = gc_INT(btop, ztmp);
2186 0 : for (i=1; i<=j; i++)
2187 0 : gmael(G,kappa,i) = submulshift(gmael(G,kappa,i), gmael(G,j,i), X, e);
2188 0 : for ( ; i<kappa; i++)
2189 0 : gmael(G,kappa,i) = submulshift(gmael(G,kappa,i), gmael(G,i,j), X, e);
2190 0 : for (i=kappa+1; i<=maxG; i++)
2191 0 : gmael(G,i,kappa) = submulshift(gmael(G,i,kappa), gmael(G,i,j), X, e);
2192 : }
2193 : }
2194 0 : if (!go_on) break; /* Anything happened? */
2195 0 : aa = zeros+1;
2196 : }
2197 :
2198 0 : affir(gmael(G,kappa,kappa), gel(s,zeros+1));
2199 : /* the last s[kappa-1]=r[kappa][kappa] is computed only if kappa increases */
2200 0 : av = avma;
2201 0 : for (k=zeros+1; k<=kappa-2; k++)
2202 0 : affrr(subrr(gel(s,k), mulrr(gmael(mu,kappa,k), gmael(r,kappa,k))),
2203 0 : gel(s,k+1));
2204 0 : *pG = G; *pB = B; *pU = U; return gc_bool(av, 0);
2205 : }
2206 :
2207 : /* G integral Gram matrix, LLL-reduces (G,B,U) in place [apply base change
2208 : * transforms to B and U]. If (keepfirst), never swap with first vector.
2209 : * If G = NULL, we compute the Gram matrix incrementally.
2210 : * Return -1 on failure, else zeros = dim Kernel (>= 0) */
2211 : static long
2212 0 : fplll(GEN *pG, GEN *pB, GEN *pU, GEN *pr, double DELTA, double ETA,
2213 : long keepfirst, long prec)
2214 : {
2215 : pari_sp av, av2;
2216 0 : GEN mu, r, s, tmp, Gtmp, alpha, G = *pG, B = *pB, U = *pU;
2217 0 : GEN delta = dbltor(DELTA), eta = dbltor(ETA);
2218 0 : long d, maxG, kappa, kappa2, i, j, zeros, kappamax, incgram = !G, cnt = 0;
2219 :
2220 0 : if (incgram)
2221 : { /* incremental Gram matrix */
2222 0 : maxG = 2; d = lg(B)-1;
2223 0 : G = zeromatcopy(d, d);
2224 : }
2225 : else
2226 0 : maxG = d = lg(G)-1;
2227 :
2228 0 : mu = cgetg(d+1, t_MAT);
2229 0 : r = cgetg(d+1, t_MAT);
2230 0 : s = cgetg(d+1, t_VEC);
2231 0 : for (j = 1; j <= d; j++)
2232 : {
2233 0 : GEN M = cgetg(d+1, t_COL), R = cgetg(d+1, t_COL);
2234 0 : long l = nbits2lg(prec);
2235 0 : gel(mu,j)= M;
2236 0 : gel(r,j) = R;
2237 0 : gel(s,j) = cgetg(l, t_REAL);
2238 0 : for (i = 1; i <= d; i++)
2239 : {
2240 0 : gel(R,i) = cgetg(l, t_REAL);
2241 0 : gel(M,i) = cgetg(l, t_REAL);
2242 : }
2243 : }
2244 0 : Gtmp = cgetg(d+1, t_VEC);
2245 0 : alpha = cgetg(d+1, t_VECSMALL);
2246 0 : av = avma;
2247 :
2248 : /* Step2: Initializing the main loop */
2249 0 : kappamax = 1;
2250 0 : i = 1;
2251 : do {
2252 0 : if (incgram) gmael(G,i,i) = ZV_dotsquare(gel(B,i));
2253 0 : affir(gmael(G,i,i), gmael(r,i,i));
2254 0 : } while (!signe(gmael(G,i,i)) && ++i <= d);
2255 0 : zeros = i-1; /* all basis vectors b_i with i <= zeros are zero vectors */
2256 0 : kappa = i;
2257 0 : for (i=zeros+1; i<=d; i++) alpha[i]=1;
2258 :
2259 0 : while (++kappa <= d)
2260 : {
2261 0 : if (kappa > kappamax)
2262 : {
2263 0 : if (DEBUGLEVEL>=4) err_printf("K%ld ",kappa);
2264 0 : kappamax = kappa;
2265 0 : if (incgram)
2266 : {
2267 0 : for (i=zeros+1; i<=kappa; i++)
2268 0 : gmael(G,kappa,i) = ZV_dotproduct(gel(B,kappa), gel(B,i));
2269 0 : maxG = kappamax;
2270 : }
2271 : }
2272 : /* Step3: Call to the Babai algorithm, mu,r,s updated in place */
2273 0 : if (Babai(av, kappa, &G,&B,&U, mu,r,s, alpha[kappa], zeros, maxG, eta, prec))
2274 0 : { *pG = incgram? NULL: G; *pB = B; *pU = U; return -1; }
2275 0 : av2 = avma;
2276 0 : if ((keepfirst && kappa == 2) ||
2277 0 : cmprr(mulrr(gmael(r,kappa-1,kappa-1), delta), gel(s,kappa-1)) <= 0)
2278 : { /* Step4: Success of Lovasz's condition */
2279 0 : alpha[kappa] = kappa;
2280 0 : tmp = mulrr(gmael(mu,kappa,kappa-1), gmael(r,kappa,kappa-1));
2281 0 : affrr(subrr(gel(s,kappa-1), tmp), gmael(r,kappa,kappa));
2282 0 : set_avma(av2); continue;
2283 : }
2284 : /* Step5: Find the right insertion index kappa, kappa2 = initial kappa */
2285 0 : if (DEBUGLEVEL>=4 && kappa==kappamax && signe(gel(s,kappa-1)))
2286 0 : if (++cnt > 20) { cnt = 0; err_printf("(%ld) ", expo(gel(s,1))); }
2287 0 : kappa2 = kappa;
2288 : do {
2289 0 : kappa--;
2290 0 : if (kappa < zeros+2 + (keepfirst ? 1: 0)) break;
2291 0 : tmp = mulrr(gmael(r,kappa-1,kappa-1), delta);
2292 0 : } while (cmprr(gel(s,kappa-1), tmp) <= 0);
2293 0 : set_avma(av2);
2294 0 : update_alpha(alpha, kappa, kappa2, kappamax);
2295 :
2296 : /* Step6: Update the mu's and r's */
2297 0 : rotate(mu, kappa2, kappa);
2298 0 : rotate(r, kappa2, kappa);
2299 0 : affrr(gel(s,kappa), gmael(r,kappa,kappa));
2300 :
2301 : /* Step7: Update G, B, U */
2302 0 : if (U) rotate(U, kappa2, kappa);
2303 0 : if (B) rotate(B, kappa2, kappa);
2304 0 : rotateG(G,kappa2,kappa, maxG, Gtmp);
2305 :
2306 : /* Step8: Prepare the next loop iteration */
2307 0 : if (kappa == zeros+1 && !signe(gmael(G,kappa,kappa)))
2308 : {
2309 0 : zeros++; kappa++;
2310 0 : affir(gmael(G,kappa,kappa), gmael(r,kappa,kappa));
2311 : }
2312 : }
2313 0 : if (pr) *pr = RgM_diagonal_shallow(r);
2314 0 : *pG = G; *pB = B; *pU = U; return zeros; /* success */
2315 : }
2316 :
2317 : /* do not support LLL_KER, LLL_ALL, LLL_KEEP_FIRST */
2318 : static GEN
2319 4851223 : ZM2_lll_norms(GEN x, long flag, GEN *pN)
2320 : {
2321 : GEN a,b,c,d;
2322 : GEN G, U;
2323 4851223 : if (flag & LLL_GRAM)
2324 7356 : G = x;
2325 : else
2326 4843867 : G = gram_matrix(x);
2327 4851223 : a = gcoeff(G,1,1); b = shifti(gcoeff(G,1,2),1); c = gcoeff(G,2,2);
2328 4851223 : d = qfb_disc3(a,b,c);
2329 4851223 : if (signe(d)>=0) return NULL;
2330 4850838 : G = redimagsl2(mkqfb(a,b,c,d),&U);
2331 4850838 : if (pN) (void) RgM_gram_schmidt(G, pN);
2332 4850838 : if (flag & LLL_INPLACE) return ZM2_mul(x,U);
2333 4850838 : return U;
2334 : }
2335 :
2336 : static void
2337 626205 : fplll_flatter(GEN *pG, GEN *pB, GEN *pU, long rank, long flag)
2338 : {
2339 626205 : if (!*pG)
2340 : {
2341 625237 : GEN T = ZM_flatter_rank(*pB, rank, flag);
2342 625237 : if (T)
2343 : {
2344 328484 : if (*pU)
2345 : {
2346 314476 : *pU = ZM_mul(*pU, T);
2347 314476 : *pB = ZM_mul(*pB, T);
2348 : }
2349 14008 : else *pB = T;
2350 : }
2351 : }
2352 : else
2353 : {
2354 968 : GEN T, G = *pG;
2355 968 : long i, j, l = lg(G);
2356 7634 : for (i = 1; i < l; i++)
2357 56193 : for(j = 1; j < i; j++) gmael(G,j,i) = gmael(G,i,j);
2358 968 : T = ZM_flattergram_rank(G, rank, flag);
2359 968 : if (T)
2360 : {
2361 968 : if (*pU) *pU = ZM_mul(*pU, T);
2362 968 : *pG = qf_ZM_apply(*pG, T);
2363 : }
2364 : }
2365 626205 : }
2366 :
2367 : static GEN
2368 1099641 : get_gramschmidt(GEN M, long rank)
2369 : {
2370 : GEN B, Q, L;
2371 1099641 : long r = lg(M)-1, prec = nbits2prec64(3*r + 30);
2372 1099641 : if (rank < r) M = vconcat(gshift(M,1), matid(r));
2373 1099641 : if (!QR_init(RgM_gtofp(M, prec), &B, &Q, &L, prec)) return NULL;
2374 475537 : return L;
2375 : }
2376 :
2377 : static GEN
2378 44547 : get_cholesky(GEN M, long rank)
2379 : {
2380 44547 : long r = lg(M)-1, prec = nbits2prec64(3*r + 30);
2381 44547 : if (rank < r) M = RgM_Rg_add(gshift(M, 1), gen_1);
2382 44547 : return RgM_Cholesky(RgM_gtofp(M, prec), prec);
2383 : }
2384 :
2385 : static long
2386 92869 : thsn(long n)
2387 : {
2388 92869 : long T[]={23280,30486,50077,44136,78724,15690,1801,1611,
2389 : 981,1359,978,1042,815,866,788,775,726,712,
2390 : 626,613,548,564,474,481,504,447,453,508,
2391 : 705,794,1008,946,767,898,886,763,842,757,
2392 : 725,774,639,655,705,627,635,704,511,613,
2393 : 583,595,568,640,541,640,567,540,577,584,
2394 : 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};
2395 92869 : return T[minss(n-3,numberof(T)-1)];
2396 : }
2397 : static long
2398 1033745 : thre(long n)
2399 : {
2400 1033745 : long T[]={31783,34393,20894,22525,13533,1928,672,671,
2401 : 422,506,315,313,222,205,167,154,139,138,
2402 : 110,120,98,94,81,75,74,64,74,74,
2403 : 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,
2404 : 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};
2405 1033745 : return T[minss(n-3,numberof(T)-1)];
2406 : }
2407 :
2408 : /* Assume x a ZM, if pN != NULL, set it to Gram-Schmidt (squared) norms
2409 : * The following modes are supported:
2410 : * - flag & LLL_INPLACE: x a lattice basis, return x*U
2411 : * - flag & LLL_GRAM: x a Gram matrix / else x a lattice basis; return
2412 : * LLL base change matrix U [LLL_IM]
2413 : * kernel basis [LLL_KER, nonreduced]
2414 : * both [LLL_ALL] */
2415 : GEN
2416 7144037 : ZM_lll_norms(GEN x, double DELTA, long flag, GEN *pN)
2417 : {
2418 7144037 : pari_sp av = avma;
2419 7144037 : const double ETA = 0.51;
2420 7144037 : const long keepfirst = flag & LLL_KEEP_FIRST;
2421 7144037 : long p, zeros = -1, n = lg(x)-1, is_upper, is_lower, useflatter, rank;
2422 7144037 : GEN G, B, U, L = NULL;
2423 : pari_timer T;
2424 7144037 : if (n <= 1) return lll_trivial(x, flag);
2425 7033969 : if (nbrows(x)==0)
2426 : {
2427 15151 : if (flag & LLL_KER) return matid(n);
2428 15151 : if (flag & (LLL_INPLACE|LLL_IM)) return cgetg(1,t_MAT);
2429 0 : retmkvec2(matid(n), cgetg(1,t_MAT));
2430 : }
2431 7018818 : if (n==2 && nbrows(x)==2 && (flag&LLL_IM) && !keepfirst)
2432 : {
2433 4851223 : U = ZM2_lll_norms(x, flag, pN);
2434 4851223 : if (U) return U;
2435 : }
2436 2167980 : if (flag & LLL_GRAM)
2437 60556 : { G = x; B = NULL; U = matid(n); is_upper = 0; is_lower = 0; }
2438 : else
2439 : {
2440 2107424 : G = NULL; B = x; U = (flag & LLL_INPLACE)? NULL: matid(n);
2441 2107424 : is_upper = (flag & LLL_UPPER) || ZM_is_upper(B);
2442 2107424 : is_lower = !B || is_upper || keepfirst ? 0: ZM_is_lower(B);
2443 2107424 : if (is_lower) L = RgM_flip(B);
2444 : }
2445 2167980 : rank = useflatter = 0;
2446 2167980 : if (n > 2 && !(flag&LLL_NOFLATTER))
2447 : {
2448 1751686 : pari_sp av2 = avma;
2449 : GEN R;
2450 1751686 : rank = ZM_rank(x);
2451 1707139 : R = B ? (is_upper ? B : (is_lower ? L : get_gramschmidt(B, rank)))
2452 3458825 : : get_cholesky(G, rank);
2453 1751686 : if (R)
2454 : {
2455 1126614 : long spr = spread(R), sz = gexpo(R), thr;
2456 1126614 : if (DEBUGLEVEL>=5)
2457 0 : err_printf("LLL: dim %ld, size %ld, spread %ld\n",n, sz, spr);
2458 1126614 : if ((is_upper && ZM_is_knapsack(B)) || (is_lower && ZM_is_knapsack(L)))
2459 92869 : thr = thsn(n);
2460 : else
2461 : {
2462 1033745 : thr = thre(n);
2463 1033745 : if (n >= 10) sz = spr;
2464 : }
2465 1126614 : useflatter = sz >= thr;
2466 : } else
2467 625072 : useflatter = 1;
2468 1751686 : set_avma(av2);
2469 : }
2470 2167980 : if(DEBUGLEVEL>=4) timer_start(&T);
2471 2167980 : if (useflatter)
2472 : {
2473 626205 : if (is_lower)
2474 : {
2475 0 : fplll_flatter(&G, &L, &U, rank, flag | LLL_UPPER);
2476 0 : B = RgM_flop(L);
2477 0 : if (U) U = RgM_flop(U);
2478 : }
2479 : else
2480 626205 : fplll_flatter(&G, &B, &U, rank, flag | (is_upper? LLL_UPPER:0));
2481 626205 : if (DEBUGLEVEL>=4 && !(flag & LLL_NOCERTIFY))
2482 0 : timer_printf(&T, "FLATTER");
2483 : }
2484 2167980 : if (!(flag & LLL_GRAM))
2485 : {
2486 : long t;
2487 2107424 : long heu_max = n<100 ? 1: 2; /* need better tuning */
2488 2107424 : B = gcopy(B);
2489 2107424 : if(DEBUGLEVEL>=4)
2490 0 : err_printf("Entering L^2 (double): dim %ld, LLL-parameters (%.3f,%.3f)\n",
2491 : n, DELTA,ETA);
2492 2107424 : zeros = fplll_fast(&B, &U, DELTA, ETA, keepfirst);
2493 2107424 : if (DEBUGLEVEL>=4) timer_printf(&T, zeros < 0? "LLL (failed)": "LLL");
2494 2112130 : for (p = DEFAULTPREC, t = 0; zeros < 0 && t < heu_max ; p += EXTRAPREC64, t++)
2495 : {
2496 4706 : if (DEBUGLEVEL>=4)
2497 0 : err_printf("Entering L^2 (heuristic): LLL-parameters (%.3f,%.3f), prec = %d/%d\n", DELTA, ETA, p, p);
2498 4706 : zeros = fplll_heuristic(&B, &U, DELTA, ETA, keepfirst, p, p);
2499 4706 : gc_lll(av, 2, &B, &U);
2500 4706 : if (DEBUGLEVEL>=4) timer_printf(&T, zeros < 0? "LLL (failed)": "LLL");
2501 : }
2502 : } else
2503 60556 : G = gcopy(G);
2504 2167980 : if (zeros < 0 || !(flag & LLL_NOCERTIFY))
2505 : {
2506 1605252 : if(DEBUGLEVEL>=4)
2507 0 : err_printf("Entering L^2 (dpe): LLL-parameters (%.3f,%.3f)\n", DELTA,ETA);
2508 1605252 : zeros = fplll_dpe(&G, &B, &U, pN, DELTA, ETA, keepfirst);
2509 1605252 : if (DEBUGLEVEL>=4) timer_printf(&T, zeros < 0? "LLL (failed)": "LLL");
2510 1605252 : if (zeros < 0)
2511 0 : for (p = DEFAULTPREC;; p += EXTRAPREC64)
2512 : {
2513 0 : if (DEBUGLEVEL>=4)
2514 0 : err_printf("Entering L^2: LLL-parameters (%.3f,%.3f), prec = %d\n",
2515 : DELTA,ETA, p);
2516 0 : zeros = fplll(&G, &B, &U, pN, DELTA, ETA, keepfirst, p);
2517 0 : if (DEBUGLEVEL>=4) timer_printf(&T, zeros < 0? "LLL (failed)": "LLL");
2518 0 : if (zeros >= 0) break;
2519 0 : gc_lll(av, 3, &G, &B, &U);
2520 : }
2521 : }
2522 2167980 : return lll_finish(U? U: B, zeros, flag);
2523 : }
2524 :
2525 : /********************************************************************/
2526 : /** **/
2527 : /** LLL OVER K[X] **/
2528 : /** **/
2529 : /********************************************************************/
2530 : static int
2531 504 : pslg(GEN x)
2532 : {
2533 : long tx;
2534 504 : if (gequal0(x)) return 2;
2535 448 : tx = typ(x); return is_scalar_t(tx)? 3: lg(x);
2536 : }
2537 :
2538 : static int
2539 196 : REDgen(long k, long l, GEN h, GEN L, GEN B)
2540 : {
2541 196 : GEN q, u = gcoeff(L,k,l);
2542 : long i;
2543 :
2544 196 : if (pslg(u) < pslg(B)) return 0;
2545 :
2546 140 : q = gneg(gdeuc(u,B));
2547 140 : gel(h,k) = gadd(gel(h,k), gmul(q,gel(h,l)));
2548 140 : for (i=1; i<l; i++) gcoeff(L,k,i) = gadd(gcoeff(L,k,i), gmul(q,gcoeff(L,l,i)));
2549 140 : gcoeff(L,k,l) = gadd(gcoeff(L,k,l), gmul(q,B)); return 1;
2550 : }
2551 :
2552 : static int
2553 196 : do_SWAPgen(GEN h, GEN L, GEN B, long k, GEN fl, int *flc)
2554 : {
2555 : GEN p1, la, la2, Bk;
2556 : long ps1, ps2, i, j, lx;
2557 :
2558 196 : if (!fl[k-1]) return 0;
2559 :
2560 140 : la = gcoeff(L,k,k-1); la2 = gsqr(la);
2561 140 : Bk = gel(B,k);
2562 140 : if (fl[k])
2563 : {
2564 56 : GEN q = gadd(la2, gmul(gel(B,k-1),gel(B,k+1)));
2565 56 : ps1 = pslg(gsqr(Bk));
2566 56 : ps2 = pslg(q);
2567 56 : if (ps1 <= ps2 && (ps1 < ps2 || !*flc)) return 0;
2568 28 : *flc = (ps1 != ps2);
2569 28 : gel(B,k) = gdiv(q, Bk);
2570 : }
2571 :
2572 112 : swap(gel(h,k-1), gel(h,k)); lx = lg(L);
2573 112 : for (j=1; j<k-1; j++) swap(gcoeff(L,k-1,j), gcoeff(L,k,j));
2574 112 : if (fl[k])
2575 : {
2576 28 : for (i=k+1; i<lx; i++)
2577 : {
2578 0 : GEN t = gcoeff(L,i,k);
2579 0 : p1 = gsub(gmul(gel(B,k+1),gcoeff(L,i,k-1)), gmul(la,t));
2580 0 : gcoeff(L,i,k) = gdiv(p1, Bk);
2581 0 : p1 = gadd(gmul(la,gcoeff(L,i,k-1)), gmul(gel(B,k-1),t));
2582 0 : gcoeff(L,i,k-1) = gdiv(p1, Bk);
2583 : }
2584 : }
2585 84 : else if (!gequal0(la))
2586 : {
2587 28 : p1 = gdiv(la2, Bk);
2588 28 : gel(B,k+1) = gel(B,k) = p1;
2589 28 : for (i=k+2; i<=lx; i++) gel(B,i) = gdiv(gmul(p1,gel(B,i)),Bk);
2590 28 : for (i=k+1; i<lx; i++)
2591 0 : gcoeff(L,i,k-1) = gdiv(gmul(la,gcoeff(L,i,k-1)), Bk);
2592 28 : for (j=k+1; j<lx-1; j++)
2593 0 : for (i=j+1; i<lx; i++)
2594 0 : gcoeff(L,i,j) = gdiv(gmul(p1,gcoeff(L,i,j)), Bk);
2595 : }
2596 : else
2597 : {
2598 56 : gcoeff(L,k,k-1) = gen_0;
2599 56 : for (i=k+1; i<lx; i++)
2600 : {
2601 0 : gcoeff(L,i,k) = gcoeff(L,i,k-1);
2602 0 : gcoeff(L,i,k-1) = gen_0;
2603 : }
2604 56 : gel(B,k) = gel(B,k-1); fl[k] = 1; fl[k-1] = 0;
2605 : }
2606 112 : return 1;
2607 : }
2608 :
2609 : static void
2610 168 : incrementalGSgen(GEN x, GEN L, GEN B, long k, GEN fl)
2611 : {
2612 168 : GEN u = NULL; /* gcc -Wall */
2613 : long i, j;
2614 420 : for (j = 1; j <= k; j++)
2615 252 : if (j==k || fl[j])
2616 : {
2617 252 : u = gcoeff(x,k,j);
2618 252 : if (!is_extscalar_t(typ(u))) pari_err_TYPE("incrementalGSgen",u);
2619 336 : for (i=1; i<j; i++)
2620 84 : if (fl[i])
2621 : {
2622 84 : u = gsub(gmul(gel(B,i+1),u), gmul(gcoeff(L,k,i),gcoeff(L,j,i)));
2623 84 : u = gdiv(u, gel(B,i));
2624 : }
2625 252 : gcoeff(L,k,j) = u;
2626 : }
2627 168 : if (gequal0(u)) gel(B,k+1) = gel(B,k);
2628 : else
2629 : {
2630 112 : gel(B,k+1) = gcoeff(L,k,k); gcoeff(L,k,k) = gen_1; fl[k] = 1;
2631 : }
2632 168 : }
2633 :
2634 : static GEN
2635 168 : lllgramallgen(GEN x, long flag)
2636 : {
2637 168 : long lx = lg(x), i, j, k, l, n;
2638 : pari_sp av;
2639 : GEN B, L, h, fl;
2640 : int flc;
2641 :
2642 168 : n = lx-1; if (n<=1) return lll_trivial(x,flag);
2643 84 : if (lgcols(x) != lx) pari_err_DIM("lllgramallgen");
2644 :
2645 84 : fl = cgetg(lx, t_VECSMALL);
2646 :
2647 84 : av = avma;
2648 84 : B = scalarcol_shallow(gen_1, lx);
2649 84 : L = cgetg(lx,t_MAT);
2650 252 : for (j=1; j<lx; j++) { gel(L,j) = zerocol(n); fl[j] = 0; }
2651 :
2652 84 : h = matid(n);
2653 252 : for (i=1; i<lx; i++)
2654 168 : incrementalGSgen(x, L, B, i, fl);
2655 84 : flc = 0;
2656 84 : for(k=2;;)
2657 : {
2658 196 : if (REDgen(k, k-1, h, L, gel(B,k))) flc = 1;
2659 196 : if (do_SWAPgen(h, L, B, k, fl, &flc)) { if (k > 2) k--; }
2660 : else
2661 : {
2662 84 : for (l=k-2; l>=1; l--)
2663 0 : if (REDgen(k, l, h, L, gel(B,l+1))) flc = 1;
2664 84 : if (++k > n) break;
2665 : }
2666 112 : if (gc_needed(av,1))
2667 : {
2668 0 : if(DEBUGMEM>1) pari_warn(warnmem,"lllgramallgen");
2669 0 : (void)gc_all(av,3,&B,&L,&h);
2670 : }
2671 : }
2672 140 : k=1; while (k<lx && !fl[k]) k++;
2673 84 : return lll_finish(h,k-1,flag);
2674 : }
2675 :
2676 : static GEN
2677 168 : lllallgen(GEN x, long flag)
2678 : {
2679 168 : pari_sp av = avma;
2680 168 : if (!(flag & LLL_GRAM)) x = gram_matrix(x);
2681 84 : else if (!RgM_is_square_mat(x)) pari_err_DIM("qflllgram");
2682 168 : return gc_GEN(av, lllgramallgen(x, flag));
2683 : }
2684 : GEN
2685 42 : lllgen(GEN x) { return lllallgen(x, LLL_IM); }
2686 : GEN
2687 42 : lllkerimgen(GEN x) { return lllallgen(x, LLL_ALL); }
2688 : GEN
2689 42 : lllgramgen(GEN x) { return lllallgen(x, LLL_IM|LLL_GRAM); }
2690 : GEN
2691 42 : lllgramkerimgen(GEN x) { return lllallgen(x, LLL_ALL|LLL_GRAM); }
2692 :
2693 : static GEN
2694 36701 : lllall(GEN x, long flag)
2695 36701 : { pari_sp av = avma; return gc_GEN(av, ZM_lll(x, LLLDFT, flag)); }
2696 : GEN
2697 183 : lllint(GEN x) { return lllall(x, LLL_IM); }
2698 : GEN
2699 35 : lllkerim(GEN x) { return lllall(x, LLL_ALL); }
2700 : GEN
2701 36441 : lllgramint(GEN x)
2702 36441 : { if (!RgM_is_square_mat(x)) pari_err_DIM("qflllgram");
2703 36441 : return lllall(x, LLL_IM | LLL_GRAM); }
2704 : GEN
2705 35 : lllgramkerim(GEN x)
2706 35 : { if (!RgM_is_square_mat(x)) pari_err_DIM("qflllgram");
2707 35 : return lllall(x, LLL_ALL | LLL_GRAM); }
2708 :
2709 : GEN
2710 5375534 : lllfp(GEN x, double D, long flag)
2711 : {
2712 5375534 : long n = lg(x)-1;
2713 5375534 : pari_sp av = avma;
2714 : GEN h;
2715 5375534 : if (n <= 1) return lll_trivial(x,flag);
2716 4713839 : if (flag & LLL_GRAM)
2717 : {
2718 9274 : if (!RgM_is_square_mat(x)) pari_err_DIM("qflllgram");
2719 9260 : if (isinexact(x))
2720 : {
2721 9169 : x = RgM_Cholesky(x, gprecision(x));
2722 9169 : if (!x) return NULL;
2723 9169 : flag &= ~LLL_GRAM;
2724 : }
2725 : }
2726 4713825 : h = ZM_lll(RgM_rescale_to_int(x), D, flag);
2727 4713769 : return gc_GEN(av, h);
2728 : }
2729 :
2730 : GEN
2731 9093 : lllgram(GEN x) { return lllfp(x,LLLDFT,LLL_GRAM|LLL_IM); }
2732 : GEN
2733 1244995 : lll(GEN x) { return lllfp(x,LLLDFT,LLL_IM); }
2734 :
2735 : static GEN
2736 63 : qflllgram(GEN x)
2737 : {
2738 63 : GEN T = lllgram(x);
2739 42 : if (!T) pari_err_PREC("qflllgram");
2740 42 : return T;
2741 : }
2742 :
2743 : GEN
2744 301 : qflll0(GEN x, long flag)
2745 : {
2746 301 : if (typ(x) != t_MAT) pari_err_TYPE("qflll",x);
2747 301 : switch(flag)
2748 : {
2749 49 : case 0: return lll(x);
2750 63 : case 1: return lllfp(x, LLLDFT, LLL_IM | LLL_NOFLATTER);
2751 49 : case 2: RgM_check_ZM(x,"qflll"); return lllintpartial(x);
2752 7 : case 3: RgM_check_ZM(x,"qflll"); return lllall(x, LLL_INPLACE);
2753 49 : case 4: RgM_check_ZM(x,"qflll"); return lllkerim(x);
2754 42 : case 5: return lllkerimgen(x);
2755 42 : case 8: return lllgen(x);
2756 0 : default: pari_err_FLAG("qflll");
2757 : }
2758 : return NULL; /* LCOV_EXCL_LINE */
2759 : }
2760 :
2761 : GEN
2762 245 : qflllgram0(GEN x, long flag)
2763 : {
2764 245 : if (typ(x) != t_MAT) pari_err_TYPE("qflllgram",x);
2765 245 : switch(flag)
2766 : {
2767 63 : case 0: return qflllgram(x);
2768 49 : case 1: return lllfp(x, LLLDFT, LLL_IM | LLL_GRAM | LLL_NOFLATTER);
2769 49 : case 4: RgM_check_ZM(x,"qflllgram"); return lllgramkerim(x);
2770 42 : case 5: return lllgramkerimgen(x);
2771 42 : case 8: return lllgramgen(x);
2772 0 : default: pari_err_FLAG("qflllgram");
2773 : }
2774 : return NULL; /* LCOV_EXCL_LINE */
2775 : }
2776 :
2777 : /********************************************************************/
2778 : /** **/
2779 : /** INTEGRAL KERNEL (LLL REDUCED) **/
2780 : /** **/
2781 : /********************************************************************/
2782 : static GEN
2783 56 : kerint0(GEN M)
2784 : {
2785 : /* return ZM_lll(M, LLLDFT, LLL_KER); */
2786 56 : GEN U, H = ZM_hnflll(M,&U,1);
2787 56 : long d = lg(M)-lg(H);
2788 56 : if (!d) return cgetg(1, t_MAT);
2789 56 : return ZM_lll(vecslice(U,1,d), LLLDFT, LLL_INPLACE);
2790 : }
2791 : GEN
2792 28 : kerint(GEN M)
2793 : {
2794 28 : pari_sp av = avma;
2795 28 : return gc_GEN(av, kerint0(M));
2796 : }
2797 : /* OBSOLETE: use kerint */
2798 : GEN
2799 28 : matkerint0(GEN M, long flag)
2800 : {
2801 28 : pari_sp av = avma;
2802 28 : if (typ(M) != t_MAT) pari_err_TYPE("matkerint",M);
2803 28 : M = Q_primpart(M);
2804 28 : RgM_check_ZM(M, "kerint");
2805 28 : switch(flag)
2806 : {
2807 28 : case 0:
2808 28 : case 1: return gc_GEN(av, kerint0(M));
2809 0 : default: pari_err_FLAG("matkerint");
2810 : }
2811 : return NULL; /* LCOV_EXCL_LINE */
2812 : }
|
