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 45301 : RgM_is_square_mat(GEN x) { long l = lg(x); return l == 1 || l == lgcols(x); }
22 :
23 : static long
24 4170970 : ZM_is_upper(GEN R)
25 : {
26 4170970 : long i,j, l = lg(R);
27 4170970 : if (l != lgcols(R)) return 0;
28 8059936 : for(i = 1; i < l; i++)
29 8699384 : for(j = 1; j < i; j++)
30 4471365 : if (signe(gcoeff(R,i,j))) return 0;
31 252191 : return 1;
32 : }
33 :
34 : static long
35 604398 : ZM_is_knapsack(GEN R)
36 : {
37 604398 : long i,j, l = lg(R);
38 604398 : if (l != lgcols(R)) return 0;
39 841184 : for(i = 2; i < l; i++)
40 2892022 : for(j = 1; j < l; j++)
41 2655236 : if ( i!=j && signe(gcoeff(R,i,j))) return 0;
42 92314 : return 1;
43 : }
44 :
45 : static long
46 1179867 : ZM_is_lower(GEN R)
47 : {
48 1179867 : long i,j, l = lg(R);
49 1179867 : if (l != lgcols(R)) return 0;
50 2056815 : for(i = 1; i < l; i++)
51 2381233 : for(j = 1; j < i; j++)
52 1290334 : if (signe(gcoeff(R,j,i))) return 0;
53 34865 : return 1;
54 : }
55 :
56 : static GEN
57 34867 : RgM_flip(GEN R)
58 : {
59 : GEN M;
60 : long i,j,l;
61 34867 : M = cgetg_copy(R, &l);
62 181048 : for(i = 1; i < l; i++)
63 : {
64 146181 : gel(M,i) = cgetg(l, t_COL);
65 908698 : for(j = 1; j < l; j++)
66 762517 : gmael(M,i,j) = gmael(R,l-i, l-j);
67 : }
68 34867 : 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 4072826 : mpabscmp(GEN x, GEN y)
89 : {
90 4072826 : 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 1336638 : drop(GEN R)
104 : {
105 1336638 : long i, n = lg(R)-1;
106 1336638 : long s = 0, m = mpexpo(gcoeff(R, 1, 1));
107 5409464 : for (i = 2; i <= n; ++i)
108 : {
109 4072826 : if (mpabscmp(gcoeff(R, i, i), gcoeff(R, i - 1, i - 1)) >= 0)
110 : {
111 2757318 : s += m - mpexpo(gcoeff(R, i - 1, i - 1));
112 2757318 : m = mpexpo(gcoeff(R, i, i));
113 : }
114 : }
115 1336638 : s += m - mpexpo(gcoeff(R, n, n));
116 1336638 : return s;
117 : }
118 :
119 : static long
120 1333629 : potential(GEN R)
121 : {
122 1333629 : long i, n = lg(R)-1;
123 1333629 : long s = 0, mul = n-1;;
124 6722857 : for (i = 1; i <= n; i++, mul-=2) s += mul * mpexpo(gcoeff(R,i,i));
125 1333629 : 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 4684570 : condition_bound(GEN U, int lower)
133 : {
134 4684570 : long n = lg(U)-1, e, i, j;
135 : GEN y;
136 4684570 : pari_sp av = avma;
137 4684570 : y = cgetg(n+1, t_VECSMALL);
138 4684571 : e = y[n] = -gexpo(gcoeff(U,n,n));
139 18616012 : for (i=n-1; i>0; i--)
140 : {
141 13931442 : long s = 0;
142 49508222 : for (j=i+1; j<=n; j++)
143 35576803 : s = maxss(s, (lower? gexpo(gcoeff(U,j,i)): gexpo(gcoeff(U,i,j))) + y[j]);
144 13931419 : y[i] = s - gexpo(gcoeff(U,i,i));
145 13931442 : e = maxss(e, y[i]);
146 : }
147 4684570 : return gc_long(av, gexpo(U) + e);
148 : }
149 :
150 : static long
151 5143936 : gsisinv(GEN M)
152 : {
153 5143936 : long i, l = lg(M);
154 25864905 : for (i = 1; i < l; ++i)
155 20721369 : if (! signe(gmael(M, i, i))) return 0;
156 5143536 : return 1;
157 : }
158 :
159 : INLINE long
160 7403306 : nbits2prec64(long n)
161 : {
162 7403306 : return nbits2prec(((n+63)>>6)<<6);
163 : }
164 :
165 : static long
166 5798601 : spread(GEN R)
167 : {
168 5798601 : long i, n = lg(R)-1, m = mpexpo(gcoeff(R, 1, 1)), M = m;
169 23318274 : for (i = 2; i <= n; ++i)
170 : {
171 17519672 : long e = mpexpo(gcoeff(R, i, i));
172 17519670 : if (e < m) m = e;
173 17519670 : if (e > M) M = e;
174 : }
175 5798602 : return M - m;
176 : }
177 :
178 : static long
179 4684569 : GS_extraprec(GEN L, int lower)
180 : {
181 4684569 : long C = condition_bound(L, lower), S = spread(L), n = lg(L)-1;
182 4684571 : return maxss(2*S+2*n, C-S-2*n); /* = 2*S + 2*n + maxss(0, C-3*S-4*n) */
183 : }
184 :
185 : static GEN
186 3009 : RgM_Cholesky_dynprec(GEN M)
187 : {
188 3009 : pari_sp ltop = avma;
189 : GEN L;
190 3009 : long minprec = lg(M) + 30, bitprec = minprec, prec;
191 : while (1)
192 4933 : {
193 : long mbitprec;
194 7942 : prec = nbits2prec64(bitprec);
195 7942 : L = RgM_Cholesky(RgM_gtofp(M, prec), prec); /* upper-triangular */
196 7942 : if (!L)
197 : {
198 1472 : bitprec *= 2;
199 1472 : set_avma(ltop);
200 1472 : continue;
201 : }
202 6470 : mbitprec = minprec + GS_extraprec(L, 0);
203 6470 : if (bitprec >= mbitprec)
204 3009 : break;
205 3461 : bitprec = maxss((4*bitprec)/3, mbitprec);
206 3461 : set_avma(ltop);
207 : }
208 3009 : return gerepilecopy(ltop, L);
209 : }
210 :
211 : static GEN
212 1312 : gramschmidt_upper(GEN M)
213 : {
214 1312 : long bitprec = lg(M)-1 + 31 + GS_extraprec(M, 0);
215 1312 : return RgM_gtofp(M, nbits2prec64(bitprec));
216 : }
217 :
218 : static GEN
219 2667257 : gramschmidt_dynprec(GEN M)
220 : {
221 2667257 : pari_sp ltop = avma;
222 2667257 : long minprec = lg(M) + 30, bitprec = minprec;
223 2667257 : if (ZM_is_upper(M)) return gramschmidt_upper(M);
224 : while (1)
225 3603120 : {
226 : GEN B, Q, L;
227 6269065 : long prec = nbits2prec64(bitprec), mbitprec;
228 6269063 : if (!QR_init(RgM_gtofp(M, prec), &B, &Q, &L, prec) || !gsisinv(L))
229 : {
230 1592266 : bitprec *= 2;
231 1592266 : set_avma(ltop);
232 1592271 : continue;
233 : }
234 4676786 : mbitprec = minprec + GS_extraprec(L, 1);
235 4676789 : if (bitprec >= mbitprec)
236 2665943 : return gerepilecopy(ltop, shallowtrans(L));
237 2010846 : bitprec = maxss((4*bitprec)/3, mbitprec);
238 2010848 : set_avma(ltop);
239 : }
240 : }
241 : /* return -T1 * round(T1^-1*(R1^-1*R2)*T3) */
242 : static GEN
243 1333626 : sizered(GEN T1, GEN T3, GEN R1, GEN R2)
244 : {
245 1333626 : pari_sp ltop = avma;
246 : long e;
247 1333626 : return gerepileupto(ltop, ZM_mul(ZM_neg(T1), grndtoi(gmul(ZM_inv(T1,NULL),
248 : RgM_mul(RgM_mul(RgM_inv_upper(R1), R2), T3)), &e)));
249 : }
250 :
251 : static GEN
252 1333629 : flat(GEN M, long flag, GEN *pt_T, long *pt_s, long *pt_pot)
253 : {
254 1333629 : pari_sp ltop = avma;
255 : GEN R, R1, R2, R3, T1, T2, T3, T, S;
256 1333629 : long k = lg(M)-1, n = k>>1, n2 = k - n, m = n>>1;
257 1333629 : long keepfirst = flag & LLL_KEEP_FIRST, inplace = flag & LLL_INPLACE;
258 : /* for k = 3, we want n = 1; n2 = 2; m = 0 */
259 : /* for k = 5, n = 2; n2 = 3; m = 1 */
260 1333629 : R = gramschmidt_dynprec(M);
261 1333629 : R1 = matslice(R, 1, n, 1, n);
262 1333628 : R2 = matslice(R, 1, n, n + 1, k);
263 1333629 : R3 = matslice(R, n + 1, k, n + 1, k);
264 1333628 : T1 = lllfp(R1, 0.99, LLL_IM| LLL_UPPER| LLL_NOCERTIFY| (keepfirst ? LLL_KEEP_FIRST: 0));
265 1333629 : T3 = lllfp(R3, 0.99, LLL_IM| LLL_UPPER| LLL_NOCERTIFY);
266 1333629 : T2 = sizered(T1, T3, R1, R2);
267 1333629 : T = shallowmatconcat(mkmat22(T1,T2,gen_0,T3));
268 1333629 : M = ZM_mul(M, T);
269 1333628 : R = gramschmidt_dynprec(M);
270 1333629 : R3 = matslice(R, m + 1, m + n2, m + 1, m + n2);
271 1333629 : T3 = lllfp(R3, 0.99, LLL_IM| LLL_UPPER| LLL_NOCERTIFY);
272 2667258 : S = shallowmatconcat(diagonal(
273 575773 : m == 0 ? mkvec2(T3, matid(k - m - n2))
274 0 : : m+n2 == k ? mkvec2(matid(m), T3)
275 757856 : : mkvec3(matid(m), T3, matid(k - m - n2))));
276 1333629 : M = ZM_mul(M, S);
277 1333628 : if (!inplace) *pt_T = ZM_mul(T, S);
278 1333629 : *pt_s = drop(R);
279 1333629 : *pt_pot = potential(R);
280 1333629 : return gc_all(ltop, inplace ? 1: 2, &M, pt_T);
281 : }
282 :
283 : static GEN
284 617438 : ZM_flatter(GEN M, long flag)
285 : {
286 617438 : pari_sp ltop = avma, btop;
287 617438 : long i, n = lg(M)-1;
288 617438 : long s = -1, pot = LONG_MAX;
289 : GEN T;
290 : pari_timer ti;
291 617438 : long lti = 1;
292 617438 : long inplace = flag & LLL_INPLACE, cert = !(flag & LLL_NOCERTIFY);
293 617438 : T = matid(n);
294 617439 : if (DEBUGLEVEL>=3)
295 : {
296 0 : timer_start(&ti);
297 0 : if (cert)
298 0 : err_printf("flatter dim = %ld size = %ld\n", n, ZM_max_expi(M));
299 : }
300 617439 : btop = avma;
301 617439 : for (i = 1;;i++)
302 716190 : {
303 : long t, pot2;
304 1333629 : GEN U, M2 = flat(M, flag, &U, &t, &pot2);
305 1333629 : if (t==0) { s = t; break; }
306 759256 : if (s >= 0)
307 : {
308 434681 : if (s==t && pot>=pot2)
309 43066 : break;
310 391615 : if (s<t && i > 20)
311 : {
312 0 : if (DEBUGLEVEL >= 3) err_printf("BACK:%ld:%ld:%g\n", n, i, s);
313 0 : break;
314 : }
315 : }
316 716190 : if (DEBUGLEVEL>=3 && (cert || timer_get(&ti) > 1000))
317 : {
318 0 : if (i==lti)
319 0 : timer_printf(&ti, "FLATTER, dim %ld, step %ld: \t slope=%0.10g \t pot=%0.10g", n, i, ((double)t)/n, ((double)pot2)/(n*(n+1)));
320 : else
321 0 : timer_printf(&ti, "FLATTER, dim %ld, steps %ld-%ld: \t slope=%0.10g \t pot=%0.10g", n, lti,i, ((double)t)/n, ((double)pot2)/(n*(n+1)));
322 0 : lti = i+1;
323 : }
324 716190 : s = t;
325 716190 : pot = pot2;
326 716190 : M = M2;
327 716190 : if (!inplace) T = ZM_mul(T, U);
328 716190 : if (gc_needed(btop, 1))
329 0 : gerepileall(btop, 2, &M, &T);
330 : }
331 617439 : if (DEBUGLEVEL>=3 && (cert || timer_get(&ti) > 1000))
332 : {
333 0 : if (i==lti)
334 0 : timer_printf(&ti, "FLATTER, dim %ld, final \t slope=%0.10g \t pot=%0.10g", n, ((double)s)/n, ((double)pot)/(n*(n+1)));
335 : else
336 0 : timer_printf(&ti, "FLATTER, dim %ld, steps %ld-final:\t slope=%0.10g \t pot=%0.10g", n, lti, ((double)s)/n, ((double)pot)/(n*(n+1)));
337 : }
338 617439 : return gerepilecopy(ltop, inplace ? M: T);
339 : }
340 :
341 : static GEN
342 615437 : ZM_flatter_rank(GEN M, long rank, long flag)
343 : {
344 : pari_timer ti;
345 615437 : pari_sp ltop = avma;
346 : GEN T;
347 615437 : long i, n = lg(M)-1;
348 615437 : if (rank == n) return ZM_flatter(M, flag);
349 3771 : T = matid(n);
350 3771 : if (DEBUGLEVEL>=3) timer_start(&ti);
351 3771 : for (i = 1;; i++)
352 2001 : {
353 5772 : GEN S = ZM_flatter(vconcat(gshift(M,i),matid(n)), flag);
354 5772 : if (DEBUGLEVEL>=3)
355 0 : timer_printf(&ti,"FLATTERRANK step %ld: %ld",i,expi(gnorml2(S)));
356 5772 : if (ZM_isidentity(S)) break;
357 2001 : T = ZM_mul(T, S);
358 2001 : M = ZM_mul(M, S);
359 2001 : if (gc_needed(ltop, 1))
360 0 : gerepileall(ltop, 2, &M, &T);
361 : }
362 3771 : return gerepilecopy(ltop, T);
363 : }
364 :
365 : static GEN
366 3009 : flattergram_i(GEN M, long flag, long *pt_s)
367 : {
368 3009 : pari_sp ltop = avma;
369 : GEN R, T;
370 3009 : R = RgM_Cholesky_dynprec(M);
371 3009 : *pt_s = drop(R);
372 3009 : T = lllfp(R, 0.99, LLL_IM| LLL_UPPER| LLL_NOCERTIFY | (flag&LLL_KEEP_FIRST));
373 3009 : return gerepilecopy(ltop, T);
374 : }
375 :
376 : static GEN
377 975 : ZM_flattergram(GEN M, long flag)
378 : {
379 975 : pari_sp ltop = avma, btop;
380 : GEN T;
381 975 : long i, n = lg(M)-1;
382 : pari_timer ti;
383 975 : long s = -1;
384 975 : if (DEBUGLEVEL>=3)
385 : {
386 0 : timer_start(&ti);
387 0 : err_printf("FLATTERGRAM dim = %ld size = %ld\n", n, ZM_max_expi(M));
388 : }
389 975 : btop = avma;
390 975 : T = matid(n);
391 975 : for (i = 1;; i++)
392 2034 : {
393 : long t;
394 3009 : GEN S = flattergram_i(M, flag, &t);
395 3009 : t = expi(gnorml2(S));
396 3009 : if (t==0) { s = t; break; }
397 3009 : if (s)
398 : {
399 3009 : double st = s-t;
400 3009 : if (st == 0)
401 975 : break;
402 2034 : if (st < 0 && i > 20)
403 : {
404 0 : if (DEBUGLEVEL >= 3) err_printf("BACK:%ld:%ld:%0.10g\n", n, i, ((double)s)/n);
405 0 : break;
406 : }
407 : }
408 2034 : T = ZM_mul(T, S);
409 2034 : M = ZM_mul(shallowtrans(S), ZM_mul(M, S));
410 2034 : s = t;
411 2034 : if (DEBUGLEVEL >= 3)
412 0 : timer_printf(&ti, "FLATTERGRAM, dim %ld step %ld, slope=%0.10g", n, i, ((double)s)/n);
413 2034 : if (gc_needed(btop, 1))
414 0 : gerepileall(btop, 2, &M, &T);
415 : }
416 975 : if (DEBUGLEVEL >= 3)
417 0 : timer_printf(&ti, "FLATTERGRAM, dim %ld, slope=%0.10g", n, ((double)s)/n);
418 975 : return gerepilecopy(ltop, T);
419 : }
420 :
421 : static GEN
422 975 : ZM_flattergram_rank(GEN M, long rank, long flag)
423 : {
424 : pari_timer ti;
425 975 : pari_sp ltop = avma;
426 : GEN T;
427 975 : long i, n = lg(M)-1;
428 975 : if (rank == n) return ZM_flattergram(M, flag);
429 0 : T = matid(n);
430 0 : if (DEBUGLEVEL>=3) timer_start(&ti);
431 0 : for (i = 1;; i++)
432 0 : {
433 0 : GEN S = ZM_flattergram(RgM_Rg_add(gshift(M, i), gen_1), flag);
434 0 : if (DEBUGLEVEL>=3)
435 0 : timer_printf(&ti,"FLATTERGRAMRANK step %ld: %ld",i,expi(gnorml2(S)));
436 0 : if (ZM_isidentity(S)) break;
437 0 : T = ZM_mul(T, S);
438 0 : M = ZM_mul(shallowtrans(S), ZM_mul(M, S));
439 0 : if (gc_needed(ltop, 1))
440 0 : gerepileall(ltop, 2, &M, &T);
441 : }
442 0 : return gerepilecopy(ltop, T);
443 : }
444 :
445 : static double
446 10650713 : pari_rint(double a)
447 : {
448 : #ifdef HAS_RINT
449 10650713 : return rint(a);
450 : #else
451 : const double two_to_52 = 4.5035996273704960e+15;
452 : double fa = fabs(a);
453 : double r = two_to_52 + fa;
454 : if (fa >= two_to_52) {
455 : r = a;
456 : } else {
457 : r = r - two_to_52;
458 : if (a < 0) r = -r;
459 : }
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 768903 : lll_trivial(GEN x, long flag)
470 : {
471 768903 : if (lg(x) == 1)
472 : { /* dim x = 0 */
473 15456 : 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 753447 : if (gequal0(gel(x,1)))
478 : {
479 91 : if (flag & LLL_KER) return matid(1);
480 91 : if (flag & (LLL_IM|LLL_INPLACE)) return cgetg(1,t_MAT);
481 28 : retmkvec2(matid(1), cgetg(1,t_MAT));
482 : }
483 753355 : if (flag & LLL_INPLACE) return gcopy(x);
484 649995 : if (flag & LLL_KER) return cgetg(1,t_MAT);
485 649995 : if (flag & LLL_IM) return matid(1);
486 29 : 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 2049860 : vectail_inplace(GEN x, long k)
492 : {
493 2049860 : if (!k) return x;
494 57656 : x[k] = ((ulong)x[0] & ~LGBITS) | _evallg(lg(x) - k);
495 57656 : return x + k;
496 : }
497 :
498 : /* k = dim Kernel */
499 : static GEN
500 2123063 : lll_finish(GEN h, long k, long flag)
501 : {
502 : GEN g;
503 2123063 : if (!(flag & (LLL_IM|LLL_KER|LLL_ALL|LLL_INPLACE))) return h;
504 2049886 : if (flag & (LLL_IM|LLL_INPLACE)) return vectail_inplace(h, k);
505 97 : if (flag & LLL_KER) { setlg(h,k+1); return h; }
506 69 : 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 314194 : mulshift(GEN y, GEN z, long e)
513 : {
514 314194 : long ly = lgefint(y), lz;
515 : pari_sp av;
516 : GEN t;
517 314194 : if (ly == 2) return gen_0;
518 46123 : lz = lgefint(z);
519 46123 : av = avma; (void)new_chunk(ly+lz+nbits2lg(e)); /* HACK */
520 46123 : t = mulii(z, y);
521 46123 : 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 1482516 : submulshift(GEN x, GEN y, GEN z, long e)
527 : {
528 1482516 : long lx = lgefint(x), ly, lz;
529 : pari_sp av;
530 : GEN t;
531 1482516 : if (!e) return submulii(x, y, z);
532 1473343 : if (lx == 2) { t = mulshift(y, z, e); togglesign(t); return t; }
533 1178151 : ly = lgefint(y);
534 1178151 : if (ly == 2) return icopy(x);
535 944087 : lz = lgefint(z);
536 944087 : av = avma; (void)new_chunk(lx+ly+lz+nbits2lg(e)); /* HACK */
537 944087 : t = shifti(mulii(z, y), e);
538 944087 : set_avma(av); return subii(x, t);
539 : }
540 : static void
541 18508981 : subzi(GEN *a, GEN b)
542 : {
543 18508981 : pari_sp av = avma;
544 18508981 : b = subii(*a, b);
545 18508934 : if (lgefint(b)<=lg(*a) && isonstack(*a)) { affii(b,*a); set_avma(av); }
546 2102960 : else *a = b;
547 18508985 : }
548 :
549 : static void
550 17765234 : addzi(GEN *a, GEN b)
551 : {
552 17765234 : pari_sp av = avma;
553 17765234 : b = addii(*a, b);
554 17765215 : if (lgefint(b)<=lg(*a) && isonstack(*a)) { affii(b,*a); set_avma(av); }
555 1875272 : else *a = b;
556 17765264 : }
557 :
558 : /* x - u*y * 2^e */
559 : INLINE GEN
560 4138030 : submuliu2n(GEN x, GEN y, ulong u, long e)
561 : {
562 : pari_sp av;
563 4138030 : long ly = lgefint(y);
564 4138030 : if (ly == 2) return x;
565 2839464 : av = avma;
566 2839464 : (void)new_chunk(3+ly+lgefint(x)+nbits2lg(e)); /* HACK */
567 2841098 : y = shifti(mului(u,y), e);
568 2840076 : set_avma(av); return subii(x, y);
569 : }
570 : /* *x -= u*y * 2^e */
571 : INLINE void
572 262957 : submulzu2n(GEN *x, GEN y, ulong u, long e)
573 : {
574 : pari_sp av;
575 262957 : long ly = lgefint(y);
576 262957 : if (ly == 2) return;
577 172642 : av = avma;
578 172642 : (void)new_chunk(3+ly+lgefint(*x)+nbits2lg(e)); /* HACK */
579 172642 : y = shifti(mului(u,y), e);
580 172642 : set_avma(av); return subzi(x, y);
581 : }
582 :
583 : /* x + u*y * 2^e */
584 : INLINE GEN
585 4049644 : addmuliu2n(GEN x, GEN y, ulong u, long e)
586 : {
587 : pari_sp av;
588 4049644 : long ly = lgefint(y);
589 4049644 : if (ly == 2) return x;
590 2784453 : av = avma;
591 2784453 : (void)new_chunk(3+ly+lgefint(x)+nbits2lg(e)); /* HACK */
592 2785977 : y = shifti(mului(u,y), e);
593 2784990 : set_avma(av); return addii(x, y);
594 : }
595 :
596 : /* *x += u*y * 2^e */
597 : INLINE void
598 270968 : addmulzu2n(GEN *x, GEN y, ulong u, long e)
599 : {
600 : pari_sp av;
601 270968 : long ly = lgefint(y);
602 270968 : if (ly == 2) return;
603 178036 : av = avma;
604 178036 : (void)new_chunk(3+ly+lgefint(*x)+nbits2lg(e)); /* HACK */
605 178036 : y = shifti(mului(u,y), e);
606 178036 : set_avma(av); return addzi(x, y);
607 : }
608 :
609 : /* n < 10; gerepileall supporting &NULL arguments. Maybe rename and export ? */
610 : INLINE void
611 4669 : gc_lll(pari_sp av, int n, ...)
612 : {
613 : int i, j;
614 : GEN *gptr[10];
615 : size_t s;
616 4669 : va_list a; va_start(a, n);
617 14007 : for (i=j=0; i<n; i++)
618 : {
619 9338 : GEN *x = va_arg(a,GEN*);
620 9338 : if (*x) { gptr[j++] = x; *x = (GEN)copy_bin(*x); }
621 : }
622 4669 : va_end(a); set_avma(av);
623 11612 : for (--j; j>=0; j--) *gptr[j] = bin_copy((GENbin*)*gptr[j]);
624 4669 : s = pari_mainstack->top - pari_mainstack->bot;
625 : /* size of saved objects ~ stacksize / 4 => overflow */
626 4669 : 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 4669 : }
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 326730 : absrsmall2(GEN x)
656 : {
657 326730 : long e = expo(x), l, i;
658 326730 : if (e < 0) return 1;
659 174062 : if (e > 0 || (ulong)x[2] > (3UL << (BITS_IN_LONG-2))) return 0;
660 : /* line above assumes l > 2. OK since x != 0 */
661 58575 : l = lg(x); for (i = 3; i < l; i++) if (x[i]) return 0;
662 50620 : return 1;
663 : }
664 : /* x t_REAL; test whether |x| <= 1/2 */
665 : static int
666 563113 : absrsmall(GEN x)
667 : {
668 : long e, l, i;
669 563113 : if (!signe(x)) return 1;
670 558960 : e = expo(x); if (e < -1) return 1;
671 332140 : if (e > -1 || (ulong)x[2] > HIGHBIT) return 0;
672 6119 : l = lg(x); for (i = 3; i < l; i++) if (x[i]) return 0;
673 5410 : return 1;
674 : }
675 :
676 : static void
677 31500950 : rotate(GEN A, long k2, long k)
678 : {
679 : long i;
680 31500950 : GEN B = gel(A,k2);
681 100603768 : for (i = k2; i > k; i--) gel(A,i) = gel(A,i-1);
682 31500950 : gel(A,k) = B;
683 31500950 : }
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 34362630 : cget_dblvec(long d)
691 34362630 : { return (double*) stack_malloc_align(d*sizeof(double), sizeof(double)); }
692 :
693 : static double **
694 8252993 : cget_dblmat(long d) { return (double **) cgetg(d, t_VECSMALL); }
695 :
696 : static double
697 160891908 : itodbl_exp(GEN x, long *e)
698 : {
699 160891908 : pari_sp av = avma;
700 160891908 : GEN r = itor(x,DEFAULTPREC);
701 160880407 : *e = expo(r); setexpo(r,0);
702 160877453 : return gc_double(av, rtodbl(r));
703 : }
704 :
705 : static double
706 117055615 : dbldotproduct(double *x, double *y, long n)
707 : {
708 : long i;
709 117055615 : double sum = del(x,1) * del(y,1);
710 1376845812 : for (i=2; i<=n; i++) sum += del(x,i) * del(y,i);
711 117055615 : return sum;
712 : }
713 :
714 : static double
715 2426145 : dbldotsquare(double *x, long n)
716 : {
717 : long i;
718 2426145 : double sum = del(x,1) * del(x,1);
719 8064345 : for (i=2; i<=n; i++) sum += del(x,i) * del(x,i);
720 2426145 : return sum;
721 : }
722 :
723 : static long
724 24564570 : set_line(double *appv, GEN v, long n)
725 : {
726 24564570 : long i, maxexp = 0;
727 24564570 : pari_sp av = avma;
728 24564570 : GEN e = cgetg(n+1, t_VECSMALL);
729 185444162 : for (i = 1; i <= n; i++)
730 : {
731 160890270 : del(appv,i) = itodbl_exp(gel(v,i), e+i);
732 160878035 : if (e[i] > maxexp) maxexp = e[i];
733 : }
734 185506129 : for (i = 1; i <= n; i++) del(appv,i) = ldexp(del(appv,i), e[i]-maxexp);
735 24553892 : set_avma(av); return maxexp;
736 : }
737 :
738 : static void
739 34250613 : dblrotate(double **A, long k2, long k)
740 : {
741 : long i;
742 34250613 : double *B = del(A,k2);
743 108511127 : for (i = k2; i > k; i--) del(A,i) = del(A,i-1);
744 34250613 : del(A,k) = B;
745 34250613 : }
746 : /* update G[kappa][i] from appB */
747 : static void
748 22345522 : setG_fast(double **appB, long n, double **G, long kappa, long a, long b)
749 : { long i;
750 100816542 : for (i = a; i <= b; i++)
751 78471081 : dmael(G,kappa,i) = dbldotproduct(del(appB,kappa), del(appB,i), n);
752 22345461 : }
753 : /* update G[i][kappa] from appB */
754 : static void
755 16802481 : setG2_fast(double **appB, long n, double **G, long kappa, long a, long b)
756 : { long i;
757 55389590 : for (i = a; i <= b; i++)
758 38587092 : dmael(G,i,kappa) = dbldotproduct(del(appB,kappa), del(appB,i), n);
759 16802498 : }
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 19635274 : u64toi(u64 x)
774 : {
775 : GEN y;
776 : ulong h;
777 19635274 : if (!x) return gen_0;
778 19635274 : h = x>>32;
779 19635274 : if (!h) return utoipos(x);
780 1128778 : y = cgetipos(4);
781 1128778 : *int_LSW(y) = x&0xFFFFFFFF;
782 1128778 : *int_MSW(y) = x>>32;
783 1128778 : return y;
784 : }
785 :
786 : INLINE GEN
787 662927 : u64toineg(u64 x)
788 : {
789 : GEN y;
790 : ulong h;
791 662927 : if (!x) return gen_0;
792 662927 : h = x>>32;
793 662927 : if (!h) return utoineg(x);
794 662927 : y = cgetineg(4);
795 662927 : *int_LSW(y) = x&0xFFFFFFFF;
796 662927 : *int_MSW(y) = x>>32;
797 662927 : return y;
798 : }
799 : INLINE GEN
800 9453030 : addmuliu64_inplace(GEN x, GEN y, u64 u) { return addmulii(x, y, u64toi(u)); }
801 :
802 : INLINE GEN
803 9503081 : submuliu64_inplace(GEN x, GEN y, u64 u) { return submulii(x, y, u64toi(u)); }
804 :
805 : INLINE GEN
806 662927 : addmuliu642n(GEN x, GEN y, u64 u, long e) { return submulshift(x, y, u64toineg(u), e); }
807 :
808 : INLINE GEN
809 679163 : 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 29881047 : 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 29881047 : GEN B = *pB, U = *pU;
820 29881047 : const long n = nbrows(B), d = U ? lg(U)-1: 0;
821 29880962 : long k, aa = (a > zeros)? a : zeros+1;
822 29880962 : long emaxmu = EX0, emax2mu = EX0;
823 : s64 xx;
824 29880962 : int did_something = 0;
825 : /* N.B: we set d = 0 (resp. n = 0) to avoid updating U (resp. B) */
826 :
827 17008115 : for (;;) {
828 46889077 : int go_on = 0;
829 46889077 : long i, j, emax3mu = emax2mu;
830 :
831 46889077 : if (gc_needed(av,2))
832 : {
833 195 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[1], a=%ld", aa);
834 195 : gc_lll(av,2,&B,&U);
835 : }
836 : /* Step2: compute the GSO for stage kappa */
837 46888394 : emax2mu = emaxmu; emaxmu = EX0;
838 179795621 : for (j=aa; j<kappa; j++)
839 : {
840 132907653 : double g = dmael(G,kappa,j);
841 571721560 : for (k = zeros+1; k < j; k++) g -= dmael(mu,j,k) * dmael(r,kappa,k);
842 132907653 : dmael(r,kappa,j) = g;
843 132907653 : dmael(mu,kappa,j) = dmael(r,kappa,j) / dmael(r,j,j);
844 132907653 : emaxmu = maxss(emaxmu, expoB[kappa]-expoB[j]);
845 : }
846 : /* maxmu doesn't decrease fast enough */
847 46887968 : if (emax3mu != EX0 && emax3mu <= emax2mu + 5) {*pB = B; *pU = U; return 1;}
848 :
849 167024178 : for (j=kappa-1; j>zeros; j--)
850 : {
851 137148455 : double tmp = fabs(ldexp (dmael(mu,kappa,j), expoB[kappa]-expoB[j]));
852 137148455 : if (tmp>eta) { go_on = 1; break; }
853 : }
854 :
855 : /* Step3--5: compute the X_j's */
856 46883853 : if (go_on)
857 77298132 : for (j=kappa-1; j>zeros; j--)
858 : { /* The code below seemingly handles U = NULL, but in this case d = 0 */
859 60291716 : int e = expoB[j] - expoB[kappa];
860 60291716 : double tmp = ldexp(dmael(mu,kappa,j), -e), atmp = fabs(tmp);
861 : /* tmp = Inf is allowed */
862 60291716 : if (atmp <= .5) continue; /* size-reduced */
863 33786065 : if (gc_needed(av,2))
864 : {
865 348 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[2], a=%ld, j=%ld", aa,j);
866 348 : gc_lll(av,2,&B,&U);
867 : }
868 33788006 : did_something = 1;
869 : /* we consider separately the case |X| = 1 */
870 33788006 : if (atmp <= 1.5)
871 : {
872 22773138 : if (dmael(mu,kappa,j) > 0) { /* in this case, X = 1 */
873 46688227 : for (k=zeros+1; k<j; k++)
874 35065236 : dmael(mu,kappa,k) -= ldexp(dmael(mu,j,k), e);
875 157465742 : for (i=1; i<=n; i++)
876 145844617 : gmael(B,kappa,i) = subii(gmael(B,kappa,i), gmael(B,j,i));
877 104283929 : for (i=1; i<=d; i++)
878 92662706 : gmael(U,kappa,i) = subii(gmael(U,kappa,i), gmael(U,j,i));
879 : } else { /* otherwise X = -1 */
880 45859131 : for (k=zeros+1; k<j; k++)
881 34708984 : dmael(mu,kappa,k) += ldexp(dmael(mu,j,k), e);
882 154779312 : for (i=1; i<=n; i++)
883 143631199 : gmael(B,kappa,i) = addii(gmael(B,kappa,i), gmael(B,j,i));
884 101583599 : for (i=1; i<=d; i++)
885 90435402 : gmael(U,kappa,i) = addii(gmael(U,kappa,i), gmael(U,j,i));
886 : }
887 22769420 : continue;
888 : }
889 : /* we have |X| >= 2 */
890 11014868 : if (atmp < 9007199254740992.)
891 : {
892 10184704 : tmp = pari_rint(tmp);
893 24386096 : for (k=zeros+1; k<j; k++)
894 14201390 : dmael(mu,kappa,k) -= ldexp(tmp * dmael(mu,j,k), e);
895 10184706 : xx = (s64) tmp;
896 10184706 : if (xx > 0) /* = xx */
897 : {
898 45869743 : for (i=1; i<=n; i++)
899 40746890 : gmael(B,kappa,i) = submuliu64_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
900 33034072 : for (i=1; i<=d; i++)
901 27911185 : gmael(U,kappa,i) = submuliu64_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
902 : }
903 : else /* = -xx */
904 : {
905 45540125 : for (i=1; i<=n; i++)
906 40478348 : gmael(B,kappa,i) = addmuliu64_inplace(gmael(B,kappa,i), gmael(B,j,i), -xx);
907 32701070 : for (i=1; i<=d; i++)
908 27639262 : gmael(U,kappa,i) = addmuliu64_inplace(gmael(U,kappa,i), gmael(U,j,i), -xx);
909 : }
910 : }
911 : else
912 : {
913 : int E;
914 830164 : xx = (s64) ldexp(frexp(dmael(mu,kappa,j), &E), 53);
915 830164 : E -= e + 53;
916 830164 : 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 2756358 : for (k=zeros+1; k<j; k++)
938 1926194 : dmael(mu,kappa,k) -= ldexp(((double)xx) * dmael(mu,j,k), E + e);
939 830164 : if (xx > 0) /* = xx */
940 : {
941 3889076 : for (i=1; i<=n; i++)
942 3471261 : gmael(B,kappa,i) = submuliu642n(gmael(B,kappa,i), gmael(B,j,i), xx, E);
943 1507779 : for (i=1; i<=d; i++)
944 1089964 : gmael(U,kappa,i) = submuliu642n(gmael(U,kappa,i), gmael(U,j,i), xx, E);
945 : }
946 : else /* = -xx */
947 : {
948 3842256 : for (i=1; i<=n; i++)
949 3429837 : gmael(B,kappa,i) = addmuliu642n(gmael(B,kappa,i), gmael(B,j,i), -xx, E);
950 1485079 : for (i=1; i<=d; i++)
951 1072660 : gmael(U,kappa,i) = addmuliu642n(gmael(U,kappa,i), gmael(U,j,i), -xx, E);
952 : }
953 : }
954 : }
955 : }
956 46882137 : if (!go_on) break; /* Anything happened? */
957 17005323 : expoB[kappa] = set_line(del(appB,kappa), gel(B,kappa), n);
958 17007785 : setG_fast(appB, n, G, kappa, zeros+1, kappa-1);
959 17008115 : aa = zeros+1;
960 : }
961 29876814 : if (did_something) setG2_fast(appB, n, G, kappa, kappa, maxG);
962 :
963 29877058 : del(s,zeros+1) = dmael(G,kappa,kappa);
964 : /* the last s[kappa-1]=r[kappa][kappa] is computed only if kappa increases */
965 108927529 : for (k=zeros+1; k<=kappa-2; k++)
966 79050471 : del(s,k+1) = del(s,k) - dmael(mu,kappa,k)*dmael(r,kappa,k);
967 29877058 : *pB = B; *pU = U; return 0;
968 : }
969 :
970 : static void
971 11837457 : update_alpha(GEN alpha, long kappa, long kappa2, long kappamax)
972 : {
973 : long i;
974 37634602 : for (i = kappa; i < kappa2; i++)
975 25797145 : if (kappa <= alpha[i]) alpha[i] = kappa;
976 37634620 : for (i = kappa2; i > kappa; i--) alpha[i] = alpha[i-1];
977 22992888 : for (i = kappa2+1; i <= kappamax; i++)
978 11155431 : if (kappa < alpha[i]) alpha[i] = kappa;
979 11837457 : alpha[kappa] = kappa;
980 11837457 : }
981 : static void
982 420465 : rotateG(GEN G, long kappa2, long kappa, long maxG, GEN Gtmp)
983 : {
984 : long i, j;
985 3292694 : for (i=1; i<=kappa2; i++) gel(Gtmp,i) = gmael(G,kappa2,i);
986 1775030 : for ( ; i<=maxG; i++) gel(Gtmp,i) = gmael(G,i,kappa2);
987 1463791 : for (i=kappa2; i>kappa; i--)
988 : {
989 5149446 : for (j=1; j<kappa; j++) gmael(G,i,j) = gmael(G,i-1,j);
990 1043326 : gmael(G,i,kappa) = gel(Gtmp,i-1);
991 3771592 : for (j=kappa+1; j<=i; j++) gmael(G,i,j) = gmael(G,i-1,j-1);
992 4691063 : for (j=kappa2+1; j<=maxG; j++) gmael(G,j,i) = gmael(G,j,i-1);
993 : }
994 1828903 : for (i=1; i<kappa; i++) gmael(G,kappa,i) = gel(Gtmp,i);
995 420465 : gmael(G,kappa,kappa) = gel(Gtmp,kappa2);
996 1775030 : for (i=kappa2+1; i<=maxG; i++) gmael(G,i,kappa) = gel(Gtmp,i);
997 420465 : }
998 : static void
999 11416873 : rotateG_fast(double **G, long kappa2, long kappa, long maxG, double *Gtmp)
1000 : {
1001 : long i, j;
1002 66331677 : for (i=1; i<=kappa2; i++) del(Gtmp,i) = dmael(G,kappa2,i);
1003 22079940 : for ( ; i<=maxG; i++) del(Gtmp,i) = dmael(G,i,kappa2);
1004 36170378 : for (i=kappa2; i>kappa; i--)
1005 : {
1006 69342350 : for (j=1; j<kappa; j++) dmael(G,i,j) = dmael(G,i-1,j);
1007 24753505 : dmael(G,i,kappa) = del(Gtmp,i-1);
1008 84378460 : for (j=kappa+1; j<=i; j++) dmael(G,i,j) = dmael(G,i-1,j-1);
1009 46712791 : for (j=kappa2+1; j<=maxG; j++) dmael(G,j,i) = dmael(G,j,i-1);
1010 : }
1011 30161526 : for (i=1; i<kappa; i++) dmael(G,kappa,i) = del(Gtmp,i);
1012 11416873 : dmael(G,kappa,kappa) = del(Gtmp,kappa2);
1013 22079980 : for (i=kappa2+1; i<=maxG; i++) dmael(G,i,kappa) = del(Gtmp,i);
1014 11416873 : }
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 2063252 : 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 2063252 : GEN alpha, expoB, B = *pB, U;
1027 2063252 : long cnt = 0;
1028 :
1029 2063252 : d = lg(B)-1;
1030 2063252 : n = nbrows(B);
1031 2063255 : U = *pU; /* NULL if inplace */
1032 :
1033 2063255 : G = cget_dblmat(d+1);
1034 2063253 : appB = cget_dblmat(d+1);
1035 2063248 : mu = cget_dblmat(d+1);
1036 2063248 : r = cget_dblmat(d+1);
1037 2063249 : s = cget_dblvec(d+1);
1038 9622357 : for (j = 1; j <= d; j++)
1039 : {
1040 7559102 : del(mu,j) = cget_dblvec(d+1);
1041 7559100 : del(r,j) = cget_dblvec(d+1);
1042 7559095 : del(appB,j) = cget_dblvec(n+1);
1043 7559085 : del(G,j) = cget_dblvec(d+1);
1044 46559642 : for (i=1; i<=d; i++) dmael(G,j,i) = 0.;
1045 : }
1046 2063255 : expoB = cgetg(d+1, t_VECSMALL);
1047 9622223 : for (i=1; i<=d; i++) expoB[i] = set_line(del(appB,i), gel(B,i), n);
1048 2063204 : Gtmp = cget_dblvec(d+1);
1049 2063249 : alpha = cgetg(d+1, t_VECSMALL);
1050 2063248 : av = avma;
1051 :
1052 : /* Step2: Initializing the main loop */
1053 2063248 : kappamax = 1;
1054 2063248 : i = 1;
1055 2063248 : maxG = d; /* later updated to kappamax */
1056 :
1057 : do {
1058 2217948 : dmael(G,i,i) = dbldotsquare(del(appB,i),n);
1059 2217948 : } while (dmael(G,i,i) <= 0 && (++i <=d));
1060 2063248 : zeros = i-1; /* all vectors B[i] with i <= zeros are zero vectors */
1061 2063248 : kappa = i;
1062 2063248 : if (zeros < d) dmael(r,zeros+1,zeros+1) = dmael(G,zeros+1,zeros+1);
1063 9467660 : for (i=zeros+1; i<=d; i++) alpha[i]=1;
1064 31940243 : while (++kappa <= d)
1065 : {
1066 29881097 : if (kappa > kappamax)
1067 : {
1068 5337659 : if (DEBUGLEVEL>=4) err_printf("K%ld ",kappa);
1069 5337659 : maxG = kappamax = kappa;
1070 5337659 : setG_fast(appB, n, G, kappa, zeros+1, kappa);
1071 : }
1072 : /* Step3: Call to the Babai algorithm, mu,r,s updated in place */
1073 29881098 : if (Babai_fast(av, kappa, &B,&U, mu,r,s, appB, expoB, G, alpha[kappa],
1074 4115 : zeros, maxG, eta)) { *pB=B; *pU=U; return -1; }
1075 :
1076 29877028 : tmp = ldexp(r[kappa-1][kappa-1] * delta, 2*(expoB[kappa-1]-expoB[kappa]));
1077 29877028 : if ((keepfirst && kappa == 2) || tmp <= del(s,kappa-1))
1078 : { /* Step4: Success of Lovasz's condition */
1079 18460036 : alpha[kappa] = kappa;
1080 18460036 : tmp = dmael(mu,kappa,kappa-1) * dmael(r,kappa,kappa-1);
1081 18460036 : dmael(r,kappa,kappa) = del(s,kappa-1)- tmp;
1082 18460036 : continue;
1083 : }
1084 : /* Step5: Find the right insertion index kappa, kappa2 = initial kappa */
1085 11416992 : 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 11416990 : kappa2 = kappa;
1088 : do {
1089 24753817 : kappa--;
1090 24753817 : if (kappa<zeros+2 + (keepfirst ? 1: 0)) break;
1091 18251878 : tmp = dmael(r,kappa-1,kappa-1) * delta;
1092 18251878 : tmp = ldexp(tmp, 2*(expoB[kappa-1]-expoB[kappa2]));
1093 18251878 : } while (del(s,kappa-1) <= tmp);
1094 11416990 : update_alpha(alpha, kappa, kappa2, kappamax);
1095 :
1096 : /* Step6: Update the mu's and r's */
1097 11416997 : dblrotate(mu,kappa2,kappa);
1098 11416987 : dblrotate(r,kappa2,kappa);
1099 11416937 : dmael(r,kappa,kappa) = del(s,kappa);
1100 :
1101 : /* Step7: Update B, appB, U, G */
1102 11416937 : rotate(B,kappa2,kappa);
1103 11416930 : dblrotate(appB,kappa2,kappa);
1104 11416900 : if (U) rotate(U,kappa2,kappa);
1105 11416899 : rotate(expoB,kappa2,kappa);
1106 11416879 : rotateG_fast(G,kappa2,kappa, maxG, Gtmp);
1107 :
1108 : /* Step8: Prepare the next loop iteration */
1109 11416959 : if (kappa == zeros+1 && dmael(G,kappa,kappa)<= 0)
1110 : {
1111 208197 : zeros++; kappa++;
1112 208197 : dmael(G,kappa,kappa) = dbldotsquare(del(appB,kappa),n);
1113 208197 : dmael(r,kappa,kappa) = dmael(G,kappa,kappa);
1114 : }
1115 : }
1116 2059146 : *pB = B; *pU = U; return zeros;
1117 : }
1118 :
1119 : /***************** HEURISTIC version (reduced precision) ****************/
1120 : static GEN
1121 157807 : realsqrdotproduct(GEN x)
1122 : {
1123 157807 : long i, l = lg(x);
1124 157807 : GEN z = sqrr(gel(x,1));
1125 1030183 : for (i=2; i<l; i++) z = addrr(z, sqrr(gel(x,i)));
1126 157807 : return z;
1127 : }
1128 : /* x, y non-empty vector of t_REALs, same length */
1129 : static GEN
1130 945582 : realdotproduct(GEN x, GEN y)
1131 : {
1132 : long i, l;
1133 : GEN z;
1134 945582 : if (x == y) return realsqrdotproduct(x);
1135 787775 : l = lg(x); z = mulrr(gel(x,1),gel(y,1));
1136 7297201 : for (i=2; i<l; i++) z = addrr(z, mulrr(gel(x,i), gel(y,i)));
1137 787775 : return z;
1138 : }
1139 : static void
1140 165784 : setG_heuristic(GEN appB, GEN G, long kappa, long a, long b)
1141 165784 : { pari_sp av = avma;
1142 : long i;
1143 759230 : for (i = a; i <= b; i++)
1144 593446 : affrr(realdotproduct(gel(appB,kappa),gel(appB,i)), gmael(G,kappa,i));
1145 165784 : set_avma(av);
1146 165784 : }
1147 : static void
1148 148216 : setG2_heuristic(GEN appB, GEN G, long kappa, long a, long b)
1149 148216 : { pari_sp av = avma;
1150 : long i;
1151 500352 : for (i = a; i <= b; i++)
1152 352136 : affrr(realdotproduct(gel(appB,kappa),gel(appB,i)), gmael(G,i,kappa));
1153 148216 : set_avma(av);
1154 148216 : }
1155 :
1156 : /* approximate t_REAL x as m * 2^e, where |m| < 2^bit */
1157 : static GEN
1158 11507 : truncexpo(GEN x, long bit, long *e)
1159 : {
1160 11507 : *e = expo(x) + 1 - bit;
1161 11507 : if (*e >= 0) return mantissa2nr(x, 0);
1162 890 : *e = 0; return roundr_safe(x);
1163 : }
1164 : /* Babai's Nearest Plane algorithm (iterative); see Babai() */
1165 : static int
1166 225089 : 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 225089 : GEN B = *pB, U = *pU;
1171 225089 : const long n = nbrows(B), d = U ? lg(U)-1: 0, bit = prec2nbits(prec);
1172 225089 : long k, aa = (a > zeros)? a : zeros+1;
1173 225089 : int did_something = 0;
1174 225089 : long emaxmu = EX0, emax2mu = EX0;
1175 : /* N.B: we set d = 0 (resp. n = 0) to avoid updating U (resp. B) */
1176 :
1177 156193 : for (;;) {
1178 381282 : int go_on = 0;
1179 381282 : long i, j, emax3mu = emax2mu;
1180 :
1181 381282 : if (gc_needed(av,2))
1182 : {
1183 3 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[1], a=%ld", aa);
1184 3 : gc_lll(av,2,&B,&U);
1185 : }
1186 : /* Step2: compute the GSO for stage kappa */
1187 381282 : emax2mu = emaxmu; emaxmu = EX0;
1188 1471435 : for (j=aa; j<kappa; j++)
1189 : {
1190 1090153 : pari_sp btop = avma;
1191 1090153 : GEN g = gmael(G,kappa,j);
1192 3749567 : for (k = zeros+1; k<j; k++)
1193 2659414 : g = subrr(g, mulrr(gmael(mu,j,k), gmael(r,kappa,k)));
1194 1090153 : affrr(g, gmael(r,kappa,j));
1195 1090153 : affrr(divrr(gmael(r,kappa,j), gmael(r,j,j)), gmael(mu,kappa,j));
1196 1090153 : emaxmu = maxss(emaxmu, expo(gmael(mu,kappa,j)));
1197 1090153 : set_avma(btop);
1198 : }
1199 381282 : if (emax3mu != EX0 && emax3mu <= emax2mu + 5)
1200 1321 : { *pB = B; *pU = U; return 1; }
1201 :
1202 1311104 : for (j=kappa-1; j>zeros; j--)
1203 1087336 : if (abscmprr(gmael(mu,kappa,j), eta) > 0) { go_on = 1; break; }
1204 :
1205 : /* Step3--5: compute the X_j's */
1206 379961 : if (go_on)
1207 719306 : 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 563113 : GEN tmp = gmael(mu,kappa,j);
1211 563113 : if (absrsmall(tmp)) continue; /* size-reduced */
1212 :
1213 326730 : if (gc_needed(av,2))
1214 : {
1215 8 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[2], a=%ld, j=%ld", aa,j);
1216 8 : gc_lll(av,2,&B,&U);
1217 : }
1218 326730 : btop = avma; did_something = 1;
1219 : /* we consider separately the case |X| = 1 */
1220 326730 : if (absrsmall2(tmp))
1221 : {
1222 203288 : if (signe(tmp) > 0) { /* in this case, X = 1 */
1223 316514 : for (k=zeros+1; k<j; k++)
1224 215215 : affrr(subrr(gmael(mu,kappa,k), gmael(mu,j,k)), gmael(mu,kappa,k));
1225 101299 : set_avma(btop);
1226 1025873 : for (i=1; i<=n; i++)
1227 924574 : gmael(B,kappa,i) = subii(gmael(B,kappa,i), gmael(B,j,i));
1228 607494 : for (i=1; i<=d; i++)
1229 506195 : gmael(U,kappa,i) = subii(gmael(U,kappa,i), gmael(U,j,i));
1230 : } else { /* otherwise X = -1 */
1231 319674 : for (k=zeros+1; k<j; k++)
1232 217685 : affrr(addrr(gmael(mu,kappa,k), gmael(mu,j,k)), gmael(mu,kappa,k));
1233 101989 : set_avma(btop);
1234 1037281 : for (i=1; i<=n; i++)
1235 935292 : gmael(B,kappa,i) = addii(gmael(B,kappa,i), gmael(B,j,i));
1236 606521 : for (i=1; i<=d; i++)
1237 504532 : gmael(U,kappa,i) = addii(gmael(U,kappa,i),gmael(U,j,i));
1238 : }
1239 203288 : continue;
1240 : }
1241 : /* we have |X| >= 2 */
1242 123442 : if (expo(tmp) < BITS_IN_LONG)
1243 : {
1244 111935 : ulong xx = roundr_safe(tmp)[2]; /* X fits in an ulong */
1245 111935 : if (signe(tmp) > 0) /* = xx */
1246 : {
1247 123217 : for (k=zeros+1; k<j; k++)
1248 66747 : affrr(subrr(gmael(mu,kappa,k), mulur(xx, gmael(mu,j,k))),
1249 66747 : gmael(mu,kappa,k));
1250 56470 : set_avma(btop);
1251 369426 : for (i=1; i<=n; i++)
1252 312956 : gmael(B,kappa,i) = submuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
1253 264462 : for (i=1; i<=d; i++)
1254 207992 : gmael(U,kappa,i) = submuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
1255 : }
1256 : else /* = -xx */
1257 : {
1258 120489 : for (k=zeros+1; k<j; k++)
1259 65024 : affrr(addrr(gmael(mu,kappa,k), mulur(xx, gmael(mu,j,k))),
1260 65024 : gmael(mu,kappa,k));
1261 55465 : set_avma(btop);
1262 362758 : for (i=1; i<=n; i++)
1263 307293 : gmael(B,kappa,i) = addmuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
1264 251704 : for (i=1; i<=d; i++)
1265 196239 : gmael(U,kappa,i) = addmuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
1266 : }
1267 : }
1268 : else
1269 : {
1270 : long e;
1271 11507 : GEN X = truncexpo(tmp, bit, &e); /* tmp ~ X * 2^e */
1272 11507 : btop = avma;
1273 27784 : for (k=zeros+1; k<j; k++)
1274 : {
1275 16277 : GEN x = mulir(X, gmael(mu,j,k));
1276 16277 : if (e) shiftr_inplace(x, e);
1277 16277 : affrr(subrr(gmael(mu,kappa,k), x), gmael(mu,kappa,k));
1278 : }
1279 11507 : set_avma(btop);
1280 93973 : for (i=1; i<=n; i++)
1281 82466 : gmael(B,kappa,i) = submulshift(gmael(B,kappa,i), gmael(B,j,i), X, e);
1282 69467 : for (i=1; i<=d; i++)
1283 57960 : gmael(U,kappa,i) = submulshift(gmael(U,kappa,i), gmael(U,j,i), X, e);
1284 : }
1285 : }
1286 379961 : if (!go_on) break; /* Anything happened? */
1287 1179300 : for (i=1 ; i<=n; i++) affir(gmael(B,kappa,i), gmael(appB,kappa,i));
1288 156193 : setG_heuristic(appB, G, kappa, zeros+1, kappa-1);
1289 156193 : aa = zeros+1;
1290 : }
1291 223768 : if (did_something) setG2_heuristic(appB, G, kappa, kappa, maxG);
1292 223768 : 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 223768 : av = avma;
1295 832493 : for (k=zeros+1; k<=kappa-2; k++)
1296 608725 : affrr(subrr(gel(s,k), mulrr(gmael(mu,kappa,k), gmael(r,kappa,k))),
1297 608725 : gel(s,k+1));
1298 223768 : *pB = B; *pU = U; return gc_bool(av, 0);
1299 : }
1300 :
1301 : static GEN
1302 15608 : ZC_to_RC(GEN x, long prec)
1303 102455 : { pari_APPLY_type(t_COL,itor(gel(x,i),prec)) }
1304 :
1305 : static GEN
1306 4115 : ZM_to_RM(GEN x, long prec)
1307 19723 : { 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 4115 : 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 4115 : GEN mu, r, s, tmp, Gtmp, alpha, G, appB, B = *pB, U;
1320 4115 : GEN delta = dbltor(DELTA), eta = dbltor(ETA);
1321 4115 : long cnt = 0;
1322 :
1323 4115 : d = lg(B)-1;
1324 4115 : U = *pU; /* NULL if inplace */
1325 :
1326 4115 : G = cgetg(d+1, t_MAT);
1327 4115 : mu = cgetg(d+1, t_MAT);
1328 4115 : r = cgetg(d+1, t_MAT);
1329 4115 : s = cgetg(d+1, t_VEC);
1330 4115 : appB = ZM_to_RM(B, prec2);
1331 19723 : for (j = 1; j <= d; j++)
1332 : {
1333 15608 : GEN M = cgetg(d+1, t_COL), R = cgetg(d+1, t_COL), S = cgetg(d+1, t_COL);
1334 15608 : gel(mu,j)= M;
1335 15608 : gel(r,j) = R;
1336 15608 : gel(G,j) = S;
1337 15608 : gel(s,j) = cgetr(prec);
1338 94224 : for (i = 1; i <= d; i++)
1339 : {
1340 78616 : gel(R,i) = cgetr(prec);
1341 78616 : gel(M,i) = cgetr(prec);
1342 78616 : gel(S,i) = cgetr(prec2);
1343 : }
1344 : }
1345 4115 : Gtmp = cgetg(d+1, t_VEC);
1346 4115 : alpha = cgetg(d+1, t_VECSMALL);
1347 4115 : av = avma;
1348 :
1349 : /* Step2: Initializing the main loop */
1350 4115 : kappamax = 1;
1351 4115 : i = 1;
1352 4115 : maxG = d; /* later updated to kappamax */
1353 :
1354 : do {
1355 4118 : affrr(RgV_dotsquare(gel(appB,i)), gmael(G,i,i));
1356 4118 : } while (signe(gmael(G,i,i)) == 0 && (++i <=d));
1357 4115 : zeros = i-1; /* all vectors B[i] with i <= zeros are zero vectors */
1358 4115 : kappa = i;
1359 4115 : if (zeros < d) affrr(gmael(G,zeros+1,zeros+1), gmael(r,zeros+1,zeros+1));
1360 19720 : for (i=zeros+1; i<=d; i++) alpha[i]=1;
1361 :
1362 227883 : while (++kappa <= d)
1363 : {
1364 225089 : if (kappa > kappamax)
1365 : {
1366 9591 : if (DEBUGLEVEL>=4) err_printf("K%ld ",kappa);
1367 9591 : maxG = kappamax = kappa;
1368 9591 : setG_heuristic(appB, G, kappa, zeros+1, kappa);
1369 : }
1370 : /* Step3: Call to the Babai algorithm, mu,r,s updated in place */
1371 225089 : if (Babai_heuristic(av, kappa, &B,&U, mu,r,s, appB, G, alpha[kappa], zeros,
1372 1321 : maxG, eta, prec)) { *pB = B; *pU = U; return -1; }
1373 223768 : av2 = avma;
1374 447438 : if ((keepfirst && kappa == 2) ||
1375 223670 : cmprr(mulrr(gmael(r,kappa-1,kappa-1), delta), gel(s,kappa-1)) <= 0)
1376 : { /* Step4: Success of Lovasz's condition */
1377 128168 : alpha[kappa] = kappa;
1378 128168 : tmp = mulrr(gmael(mu,kappa,kappa-1), gmael(r,kappa,kappa-1));
1379 128168 : affrr(subrr(gel(s,kappa-1), tmp), gmael(r,kappa,kappa));
1380 128168 : set_avma(av2); continue;
1381 : }
1382 : /* Step5: Find the right insertion index kappa, kappa2 = initial kappa */
1383 95600 : if (DEBUGLEVEL>=4 && kappa==kappamax && signe(gel(s,kappa-1)))
1384 0 : if (++cnt > 20) { cnt = 0; err_printf("(%ld) ", expo(gel(s,1))); }
1385 95600 : kappa2 = kappa;
1386 : do {
1387 215523 : kappa--;
1388 215523 : if (kappa < zeros+2 + (keepfirst ? 1: 0)) break;
1389 193659 : tmp = mulrr(gmael(r,kappa-1,kappa-1), delta);
1390 193659 : } while (cmprr(gel(s,kappa-1), tmp) <= 0 );
1391 95600 : set_avma(av2);
1392 95600 : update_alpha(alpha, kappa, kappa2, kappamax);
1393 :
1394 : /* Step6: Update the mu's and r's */
1395 95600 : rotate(mu,kappa2,kappa);
1396 95600 : rotate(r,kappa2,kappa);
1397 95600 : affrr(gel(s,kappa), gmael(r,kappa,kappa));
1398 :
1399 : /* Step7: Update B, appB, U, G */
1400 95600 : rotate(B,kappa2,kappa);
1401 95600 : rotate(appB,kappa2,kappa);
1402 95600 : if (U) rotate(U,kappa2,kappa);
1403 95600 : rotateG(G,kappa2,kappa, maxG, Gtmp);
1404 :
1405 : /* Step8: Prepare the next loop iteration */
1406 95600 : if (kappa == zeros+1 && !signe(gmael(G,kappa,kappa)))
1407 : {
1408 0 : zeros++; kappa++;
1409 0 : affrr(RgV_dotsquare(gel(appB,kappa)), gmael(G,kappa,kappa));
1410 0 : affrr(gmael(G,kappa,kappa), gmael(r,kappa,kappa));
1411 : }
1412 : }
1413 2794 : *pB=B; *pU=U; return zeros;
1414 : }
1415 :
1416 : /************************* PROVED version (t_INT) ***********************/
1417 : /* dpe inspired by dpe.h by Patrick Pelissier, Paul Zimmermann
1418 : * https://gforge.inria.fr/projects/dpe/
1419 : */
1420 :
1421 : typedef struct
1422 : {
1423 : double d; /* significand */
1424 : long e; /* exponent */
1425 : } dpe_t;
1426 :
1427 : #define Dmael(x,i,j) (&((x)[i][j]))
1428 : #define Del(x,i) (&((x)[i]))
1429 :
1430 : static void
1431 649730 : dperotate(dpe_t **A, long k2, long k)
1432 : {
1433 : long i;
1434 649730 : dpe_t *B = A[k2];
1435 2305336 : for (i = k2; i > k; i--) A[i] = A[i-1];
1436 649730 : A[k] = B;
1437 649730 : }
1438 :
1439 : static void
1440 107903044 : dpe_normalize0(dpe_t *x)
1441 : {
1442 : int e;
1443 107903044 : x->d = frexp(x->d, &e);
1444 107903044 : x->e += e;
1445 107903044 : }
1446 :
1447 : static void
1448 47851622 : dpe_normalize(dpe_t *x)
1449 : {
1450 47851622 : if (x->d == 0.0)
1451 495389 : x->e = -LONG_MAX;
1452 : else
1453 47356233 : dpe_normalize0(x);
1454 47851680 : }
1455 :
1456 : static GEN
1457 93417 : dpetor(dpe_t *x)
1458 : {
1459 93417 : GEN r = dbltor(x->d);
1460 93416 : if (signe(r)==0) return r;
1461 93416 : setexpo(r, x->e-1);
1462 93416 : return r;
1463 : }
1464 :
1465 : static void
1466 25516342 : affdpe(dpe_t *y, dpe_t *x)
1467 : {
1468 25516342 : x->d = y->d;
1469 25516342 : x->e = y->e;
1470 25516342 : }
1471 :
1472 : static void
1473 20473771 : affidpe(GEN y, dpe_t *x)
1474 : {
1475 20473771 : pari_sp av = avma;
1476 20473771 : GEN r = itor(y, DEFAULTPREC);
1477 20473413 : x->e = expo(r)+1;
1478 20473413 : setexpo(r,-1);
1479 20473403 : x->d = rtodbl(r);
1480 20473345 : set_avma(av);
1481 20473308 : }
1482 :
1483 : static void
1484 3128274 : affdbldpe(double y, dpe_t *x)
1485 : {
1486 3128274 : x->d = (double)y;
1487 3128274 : x->e = 0;
1488 3128274 : dpe_normalize(x);
1489 3128274 : }
1490 :
1491 : static void
1492 56644219 : dpe_mulz(dpe_t *x, dpe_t *y, dpe_t *z)
1493 : {
1494 56644219 : z->d = x->d * y->d;
1495 56644219 : if (z->d == 0.0)
1496 8132563 : z->e = -LONG_MAX;
1497 : else
1498 : {
1499 48511656 : z->e = x->e + y->e;
1500 48511656 : dpe_normalize0(z);
1501 : }
1502 56644422 : }
1503 :
1504 : static void
1505 13930355 : dpe_divz(dpe_t *x, dpe_t *y, dpe_t *z)
1506 : {
1507 13930355 : z->d = x->d / y->d;
1508 13930355 : if (z->d == 0.0)
1509 1894837 : z->e = -LONG_MAX;
1510 : else
1511 : {
1512 12035518 : z->e = x->e - y->e;
1513 12035518 : dpe_normalize0(z);
1514 : }
1515 13930397 : }
1516 :
1517 : static void
1518 243913 : dpe_negz(dpe_t *y, dpe_t *x)
1519 : {
1520 243913 : x->d = - y->d;
1521 243913 : x->e = y->e;
1522 243913 : }
1523 :
1524 : static void
1525 1941014 : dpe_addz(dpe_t *y, dpe_t *z, dpe_t *x)
1526 : {
1527 1941014 : if (y->e > z->e + 53)
1528 111849 : affdpe(y, x);
1529 1829165 : else if (z->e > y->e + 53)
1530 41640 : affdpe(z, x);
1531 : else
1532 : {
1533 1787525 : long d = y->e - z->e;
1534 :
1535 1787525 : if (d >= 0)
1536 : {
1537 1343506 : x->d = y->d + ldexp(z->d, -d);
1538 1343506 : x->e = y->e;
1539 : }
1540 : else
1541 : {
1542 444019 : x->d = z->d + ldexp(y->d, d);
1543 444019 : x->e = z->e;
1544 : }
1545 1787525 : dpe_normalize(x);
1546 : }
1547 1941014 : }
1548 : static void
1549 53488916 : dpe_subz(dpe_t *y, dpe_t *z, dpe_t *x)
1550 : {
1551 53488916 : if (y->e > z->e + 53)
1552 11107787 : affdpe(y, x);
1553 42381129 : else if (z->e > y->e + 53)
1554 243913 : dpe_negz(z, x);
1555 : else
1556 : {
1557 42137216 : long d = y->e - z->e;
1558 :
1559 42137216 : if (d >= 0)
1560 : {
1561 39379711 : x->d = y->d - ldexp(z->d, -d);
1562 39379711 : x->e = y->e;
1563 : }
1564 : else
1565 : {
1566 2757505 : x->d = ldexp(y->d, d) - z->d;
1567 2757505 : x->e = z->e;
1568 : }
1569 42137216 : dpe_normalize(x);
1570 : }
1571 53489177 : }
1572 :
1573 : static void
1574 798420 : dpe_muluz(dpe_t *y, ulong t, dpe_t *x)
1575 : {
1576 798420 : x->d = y->d * (double)t;
1577 798420 : x->e = y->e;
1578 798420 : dpe_normalize(x);
1579 798420 : }
1580 :
1581 : static void
1582 342388 : dpe_addmuluz(dpe_t *y, dpe_t *z, ulong t, dpe_t *x)
1583 : {
1584 : dpe_t tmp;
1585 342388 : dpe_muluz(z, t, &tmp);
1586 342388 : dpe_addz(y, &tmp, x);
1587 342388 : }
1588 :
1589 : static void
1590 411548 : dpe_submuluz(dpe_t *y, dpe_t *z, ulong t, dpe_t *x)
1591 : {
1592 : dpe_t tmp;
1593 411548 : dpe_muluz(z, t, &tmp);
1594 411548 : dpe_subz(y, &tmp, x);
1595 411548 : }
1596 :
1597 : static void
1598 51459840 : dpe_submulz(dpe_t *y, dpe_t *z, dpe_t *t, dpe_t *x)
1599 : {
1600 : dpe_t tmp;
1601 51459840 : dpe_mulz(z, t, &tmp);
1602 51459838 : dpe_subz(y, &tmp, x);
1603 51459902 : }
1604 :
1605 : static int
1606 5184498 : dpe_cmp(dpe_t *x, dpe_t *y)
1607 : {
1608 5184498 : int sx = x->d < 0. ? -1: x->d > 0.;
1609 5184498 : int sy = y->d < 0. ? -1: y->d > 0.;
1610 5184498 : int d = sx - sy;
1611 :
1612 5184498 : if (d != 0)
1613 141619 : return d;
1614 5042879 : else if (x->e > y->e)
1615 479901 : return (sx > 0) ? 1 : -1;
1616 4562978 : else if (y->e > x->e)
1617 2496366 : return (sx > 0) ? -1 : 1;
1618 : else
1619 2066612 : return (x->d < y->d) ? -1 : (x->d > y->d);
1620 : }
1621 :
1622 : static int
1623 14371449 : dpe_abscmp(dpe_t *x, dpe_t *y)
1624 : {
1625 14371449 : if (x->e > y->e)
1626 270084 : return 1;
1627 14101365 : else if (y->e > x->e)
1628 13252637 : return -1;
1629 : else
1630 848728 : return (fabs(x->d) < fabs(y->d)) ? -1 : (fabs(x->d) > fabs(y->d));
1631 : }
1632 :
1633 : static int
1634 1385497 : dpe_abssmall(dpe_t *x)
1635 : {
1636 1385497 : return (x->e <= 0) || (x->e == 1 && fabs(x->d) <= .75);
1637 : }
1638 :
1639 : static int
1640 5184490 : dpe_cmpmul(dpe_t *x, dpe_t *y, dpe_t *z)
1641 : {
1642 : dpe_t t;
1643 5184490 : dpe_mulz(x,y,&t);
1644 5184493 : return dpe_cmp(&t, z);
1645 : }
1646 :
1647 : static dpe_t *
1648 12995903 : cget_dpevec(long d)
1649 12995903 : { return (dpe_t*) stack_malloc_align(d*sizeof(dpe_t), sizeof(dpe_t)); }
1650 :
1651 : static dpe_t **
1652 3128268 : cget_dpemat(long d) { return (dpe_t **) cgetg(d, t_VECSMALL); }
1653 :
1654 : static GEN
1655 20028 : dpeM_diagonal_shallow(dpe_t **m, long d)
1656 : {
1657 : long i;
1658 20028 : GEN y = cgetg(d+1,t_VEC);
1659 113444 : for (i=1; i<=d; i++) gel(y, i) = dpetor(Dmael(m,i,i));
1660 20027 : return y;
1661 : }
1662 :
1663 : static void
1664 1385480 : affii_or_copy_gc(pari_sp av, GEN x, GEN *y)
1665 : {
1666 1385480 : long l = lg(*y);
1667 1385480 : if (lgefint(x) <= l && isonstack(*y))
1668 : {
1669 1385465 : affii(x,*y);
1670 1385467 : set_avma(av);
1671 : }
1672 : else
1673 12 : *y = gerepileuptoint(av, x);
1674 1385480 : }
1675 :
1676 : /* *x -= u*y */
1677 : INLINE void
1678 5915847 : submulziu(GEN *x, GEN y, ulong u)
1679 : {
1680 : pari_sp av;
1681 5915847 : long ly = lgefint(y);
1682 5915847 : if (ly == 2) return;
1683 3249009 : av = avma;
1684 3249009 : (void)new_chunk(3+ly+lgefint(*x)); /* HACK */
1685 3249059 : y = mului(u,y);
1686 3249045 : set_avma(av); subzi(x, y);
1687 : }
1688 :
1689 : /* *x += u*y */
1690 : INLINE void
1691 4573674 : addmulziu(GEN *x, GEN y, ulong u)
1692 : {
1693 : pari_sp av;
1694 4573674 : long ly = lgefint(y);
1695 4573674 : if (ly == 2) return;
1696 2752962 : av = avma;
1697 2752962 : (void)new_chunk(3+ly+lgefint(*x)); /* HACK */
1698 2752998 : y = mului(u,y);
1699 2752985 : set_avma(av); addzi(x, y);
1700 : }
1701 :
1702 : /************************** PROVED version (dpe) *************************/
1703 :
1704 : /* Babai's Nearest Plane algorithm (iterative).
1705 : * Size-reduces b_kappa using mu_{i,j} and r_{i,j} for j<=i <kappa
1706 : * Update B[,kappa]; compute mu_{kappa,j}, r_{kappa,j} for j<=kappa and s[kappa]
1707 : * mu, r, s updated in place (affrr). Return 1 on failure, else 0. */
1708 : static int
1709 4576176 : Babai_dpe(pari_sp av, long kappa, GEN *pG, GEN *pB, GEN *pU, dpe_t **mu, dpe_t **r, dpe_t *s,
1710 : long a, long zeros, long maxG, dpe_t *eta)
1711 : {
1712 4576176 : GEN G = *pG, B = *pB, U = *pU, ztmp;
1713 4576176 : long k, d, n, aa = a > zeros? a: zeros+1;
1714 4576176 : long emaxmu = EX0, emax2mu = EX0;
1715 : /* N.B: we set d = 0 (resp. n = 0) to avoid updating U (resp. B) */
1716 4576176 : d = U? lg(U)-1: 0;
1717 4576176 : n = B? nbrows(B): 0;
1718 520747 : for (;;) {
1719 5096943 : int go_on = 0;
1720 5096943 : long i, j, emax3mu = emax2mu;
1721 :
1722 5096943 : if (gc_needed(av,2))
1723 : {
1724 0 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[1], a=%ld", aa);
1725 0 : gc_lll(av,3,&G,&B,&U);
1726 : }
1727 : /* Step2: compute the GSO for stage kappa */
1728 5096940 : emax2mu = emaxmu; emaxmu = EX0;
1729 19027321 : for (j=aa; j<kappa; j++)
1730 : {
1731 : dpe_t g;
1732 13930383 : affidpe(gmael(G,kappa,j), &g);
1733 52268859 : for (k = zeros+1; k < j; k++)
1734 38338559 : dpe_submulz(&g, Dmael(mu,j,k), Dmael(r,kappa,k), &g);
1735 13930300 : affdpe(&g, Dmael(r,kappa,j));
1736 13930353 : dpe_divz(Dmael(r,kappa,j), Dmael(r,j,j), Dmael(mu,kappa,j));
1737 13930348 : emaxmu = maxss(emaxmu, Dmael(mu,kappa,j)->e);
1738 : }
1739 5096938 : if (emax3mu != EX0 && emax3mu <= emax2mu + 5) /* precision too low */
1740 0 : { *pG = G; *pB = B; *pU = U; return 1; }
1741 :
1742 18947654 : for (j=kappa-1; j>zeros; j--)
1743 14371474 : if (dpe_abscmp(Dmael(mu,kappa,j), eta) > 0) { go_on = 1; break; }
1744 :
1745 : /* Step3--5: compute the X_j's */
1746 5096936 : if (go_on)
1747 3026309 : for (j=kappa-1; j>zeros; j--)
1748 : {
1749 : pari_sp btop;
1750 2505556 : dpe_t *tmp = Dmael(mu,kappa,j);
1751 2505556 : if (tmp->e < 0) continue; /* (essentially) size-reduced */
1752 :
1753 1385497 : if (gc_needed(av,2))
1754 : {
1755 0 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[2], a=%ld, j=%ld", aa,j);
1756 0 : gc_lll(av,3,&G,&B,&U);
1757 : }
1758 : /* we consider separately the case |X| = 1 */
1759 1385497 : if (dpe_abssmall(tmp))
1760 : {
1761 919494 : if (tmp->d > 0) { /* in this case, X = 1 */
1762 2056641 : for (k=zeros+1; k<j; k++)
1763 1595630 : dpe_subz(Dmael(mu,kappa,k), Dmael(mu,j,k), Dmael(mu,kappa,k));
1764 3016360 : for (i=1; i<=n; i++)
1765 2555349 : subzi(&gmael(B,kappa,i), gmael(B,j,i));
1766 6967971 : for (i=1; i<=d; i++)
1767 6506966 : subzi(&gmael(U,kappa,i), gmael(U,j,i));
1768 461005 : btop = avma;
1769 461005 : ztmp = subii(gmael(G,j,j), shifti(gmael(G,kappa,j), 1));
1770 461010 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
1771 461008 : affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
1772 2857843 : for (i=1; i<=j; i++)
1773 2396833 : subzi(&gmael(G,kappa,i), gmael(G,j,i));
1774 2190696 : for (i=j+1; i<kappa; i++)
1775 1729686 : subzi(&gmael(G,kappa,i), gmael(G,i,j));
1776 2360041 : for (i=kappa+1; i<=maxG; i++)
1777 1899032 : subzi(&gmael(G,i,kappa), gmael(G,i,j));
1778 : } else { /* otherwise X = -1 */
1779 2034562 : for (k=zeros+1; k<j; k++)
1780 1576079 : dpe_addz(Dmael(mu,kappa,k), Dmael(mu,j,k), Dmael(mu,kappa,k));
1781 3011275 : for (i=1; i<=n; i++)
1782 2552792 : addzi(&gmael(B,kappa,i),gmael(B,j,i));
1783 6840485 : for (i=1; i<=d; i++)
1784 6382012 : addzi(&gmael(U,kappa,i),gmael(U,j,i));
1785 458473 : btop = avma;
1786 458473 : ztmp = addii(gmael(G,j,j), shifti(gmael(G,kappa,j), 1));
1787 458479 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
1788 458481 : affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
1789 2768423 : for (i=1; i<=j; i++)
1790 2309942 : addzi(&gmael(G,kappa,i), gmael(G,j,i));
1791 2194342 : for (i=j+1; i<kappa; i++)
1792 1735862 : addzi(&gmael(G,kappa,i), gmael(G,i,j));
1793 2312636 : for (i=kappa+1; i<=maxG; i++)
1794 1854157 : addzi(&gmael(G,i,kappa), gmael(G,i,j));
1795 : }
1796 919488 : continue;
1797 : }
1798 : /* we have |X| >= 2 */
1799 466003 : if (tmp->e < BITS_IN_LONG-1)
1800 : {
1801 447008 : if (tmp->d > 0)
1802 : {
1803 246547 : ulong xx = (ulong) pari_rint(ldexp(tmp->d, tmp->e)); /* X fits in an ulong */
1804 658095 : for (k=zeros+1; k<j; k++)
1805 411548 : dpe_submuluz(Dmael(mu,kappa,k), Dmael(mu,j,k), xx, Dmael(mu,kappa,k));
1806 721378 : for (i=1; i<=n; i++)
1807 474831 : submulziu(&gmael(B,kappa,i), gmael(B,j,i), xx);
1808 3104472 : for (i=1; i<=d; i++)
1809 2857928 : submulziu(&gmael(U,kappa,i), gmael(U,j,i), xx);
1810 246544 : btop = avma;
1811 246544 : ztmp = submuliu2n(mulii(gmael(G,j,j), sqru(xx)), gmael(G,kappa,j), xx, 1);
1812 246548 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
1813 246547 : affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
1814 1312132 : for (i=1; i<=j; i++)
1815 1065589 : submulziu(&gmael(G,kappa,i), gmael(G,j,i), xx);
1816 806775 : for (i=j+1; i<kappa; i++)
1817 560230 : submulziu(&gmael(G,kappa,i), gmael(G,i,j), xx);
1818 1203873 : for (i=kappa+1; i<=maxG; i++)
1819 957328 : submulziu(&gmael(G,i,kappa), gmael(G,i,j), xx);
1820 : }
1821 : else
1822 : {
1823 200461 : ulong xx = (ulong) pari_rint(ldexp(-tmp->d, tmp->e)); /* X fits in an ulong */
1824 542849 : for (k=zeros+1; k<j; k++)
1825 342388 : dpe_addmuluz(Dmael(mu,kappa,k), Dmael(mu,j,k), xx, Dmael(mu,kappa,k));
1826 686753 : for (i=1; i<=n; i++)
1827 486292 : addmulziu(&gmael(B,kappa,i), gmael(B,j,i), xx);
1828 2356943 : for (i=1; i<=d; i++)
1829 2156485 : addmulziu(&gmael(U,kappa,i), gmael(U,j,i), xx);
1830 200458 : btop = avma;
1831 200458 : ztmp = addmuliu2n(mulii(gmael(G,j,j), sqru(xx)), gmael(G,kappa,j), xx, 1);
1832 200460 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
1833 200459 : affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
1834 988701 : for (i=1; i<=j; i++)
1835 788241 : addmulziu(&gmael(G,kappa,i), gmael(G,j,i), xx);
1836 661300 : for (i=j+1; i<kappa; i++)
1837 460841 : addmulziu(&gmael(G,kappa,i), gmael(G,i,j), xx);
1838 882322 : for (i=kappa+1; i<=maxG; i++)
1839 681863 : addmulziu(&gmael(G,i,kappa), gmael(G,i,j), xx);
1840 : }
1841 : }
1842 : else
1843 : {
1844 18995 : long e = tmp->e - BITS_IN_LONG + 1;
1845 18995 : if (tmp->d > 0)
1846 : {
1847 9396 : ulong xx = (ulong) pari_rint(ldexp(tmp->d, BITS_IN_LONG - 1));
1848 31333 : for (k=zeros+1; k<j; k++)
1849 : {
1850 : dpe_t x;
1851 21937 : dpe_muluz(Dmael(mu,j,k), xx, &x);
1852 21937 : x.e += e;
1853 21937 : dpe_subz(Dmael(mu,kappa,k), &x, Dmael(mu,kappa,k));
1854 : }
1855 124323 : for (i=1; i<=n; i++)
1856 114927 : submulzu2n(&gmael(B,kappa,i), gmael(B,j,i), xx, e);
1857 87036 : for (i=1; i<=d; i++)
1858 77640 : submulzu2n(&gmael(U,kappa,i), gmael(U,j,i), xx, e);
1859 9396 : btop = avma;
1860 9396 : ztmp = submuliu2n(mulshift(gmael(G,j,j), sqru(xx), 2*e),
1861 9396 : gmael(G,kappa,j), xx, e+1);
1862 9396 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
1863 9396 : affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
1864 40927 : for (i=1; i<=j; i++)
1865 31531 : submulzu2n(&gmael(G,kappa,i), gmael(G,j,i), xx, e);
1866 47343 : for ( ; i<kappa; i++)
1867 37947 : submulzu2n(&gmael(G,kappa,i), gmael(G,i,j), xx, e);
1868 10308 : for (i=kappa+1; i<=maxG; i++)
1869 912 : submulzu2n(&gmael(G,i,kappa), gmael(G,i,j), xx, e);
1870 : } else
1871 : {
1872 9599 : ulong xx = (ulong) pari_rint(ldexp(-tmp->d, BITS_IN_LONG - 1));
1873 32153 : for (k=zeros+1; k<j; k++)
1874 : {
1875 : dpe_t x;
1876 22547 : dpe_muluz(Dmael(mu,j,k), xx, &x);
1877 22547 : x.e += e;
1878 22547 : dpe_addz(Dmael(mu,kappa,k), &x, Dmael(mu,kappa,k));
1879 : }
1880 128490 : for (i=1; i<=n; i++)
1881 118884 : addmulzu2n(&gmael(B,kappa,i), gmael(B,j,i), xx, e);
1882 89277 : for (i=1; i<=d; i++)
1883 79671 : addmulzu2n(&gmael(U,kappa,i), gmael(U,j,i), xx, e);
1884 9606 : btop = avma;
1885 9606 : ztmp = addmuliu2n(mulshift(gmael(G,j,j), sqru(xx), 2*e),
1886 9606 : gmael(G,kappa,j), xx, e+1);
1887 9606 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
1888 9606 : affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
1889 41957 : for (i=1; i<=j; i++)
1890 32351 : addmulzu2n(&gmael(G,kappa,i), gmael(G,j,i), xx, e);
1891 48921 : for ( ; i<kappa; i++)
1892 39315 : addmulzu2n(&gmael(G,kappa,i), gmael(G,i,j), xx, e);
1893 10353 : for (i=kappa+1; i<=maxG; i++)
1894 747 : addmulzu2n(&gmael(G,i,kappa), gmael(G,i,j), xx, e);
1895 : }
1896 : }
1897 : }
1898 5096933 : if (!go_on) break; /* Anything happened? */
1899 520747 : aa = zeros+1;
1900 : }
1901 :
1902 4576186 : affidpe(gmael(G,kappa,kappa), Del(s,zeros+1));
1903 : /* the last s[kappa-1]=r[kappa][kappa] is computed only if kappa increases */
1904 13446356 : for (k=zeros+1; k<=kappa-2; k++)
1905 8870182 : dpe_submulz(Del(s,k), Dmael(mu,kappa,k), Dmael(r,kappa,k), Del(s,k+1));
1906 4576174 : *pG = G; *pB = B; *pU = U; return 0;
1907 : }
1908 :
1909 : /* G integral Gram matrix, LLL-reduces (G,B,U) in place [apply base change
1910 : * transforms to B and U]. If (keepfirst), never swap with first vector.
1911 : * If G = NULL, we compute the Gram matrix incrementally.
1912 : * Return -1 on failure, else zeros = dim Kernel (>= 0) */
1913 : static long
1914 1564138 : fplll_dpe(GEN *pG, GEN *pB, GEN *pU, GEN *pr, double DELTA, double ETA,
1915 : long keepfirst)
1916 : {
1917 : pari_sp av;
1918 1564138 : GEN Gtmp, alpha, G = *pG, B = *pB, U = *pU;
1919 1564138 : long d, maxG, kappa, kappa2, i, j, zeros, kappamax, incgram = !G, cnt = 0;
1920 : dpe_t delta, eta, **mu, **r, *s;
1921 1564138 : affdbldpe(DELTA,&delta);
1922 1564138 : affdbldpe(ETA,&eta);
1923 :
1924 1564139 : if (incgram)
1925 : { /* incremental Gram matrix */
1926 1504403 : maxG = 2; d = lg(B)-1;
1927 1504403 : G = zeromatcopy(d, d);
1928 : }
1929 : else
1930 59736 : maxG = d = lg(G)-1;
1931 :
1932 1564135 : mu = cget_dpemat(d+1);
1933 1564134 : r = cget_dpemat(d+1);
1934 1564129 : s = cget_dpevec(d+1);
1935 7280036 : for (j = 1; j <= d; j++)
1936 : {
1937 5715904 : mu[j]= cget_dpevec(d+1);
1938 5715902 : r[j] = cget_dpevec(d+1);
1939 : }
1940 1564132 : Gtmp = cgetg(d+1, t_VEC);
1941 1564129 : alpha = cgetg(d+1, t_VECSMALL);
1942 1564122 : av = avma;
1943 :
1944 : /* Step2: Initializing the main loop */
1945 1564122 : kappamax = 1;
1946 1564122 : i = 1;
1947 : do {
1948 1932298 : if (incgram) gmael(G,i,i) = ZV_dotsquare(gel(B,i));
1949 1932262 : affidpe(gmael(G,i,i), Dmael(r,i,i));
1950 1932263 : } while (!signe(gmael(G,i,i)) && ++i <= d);
1951 1564087 : zeros = i-1; /* all basis vectors b_i with i <= zeros are zero vectors */
1952 1564087 : kappa = i;
1953 6911701 : for (i=zeros+1; i<=d; i++) alpha[i]=1;
1954 :
1955 6140302 : while (++kappa <= d)
1956 : {
1957 4576181 : if (kappa > kappamax)
1958 : {
1959 3783543 : if (DEBUGLEVEL>=4) err_printf("K%ld ",kappa);
1960 3783543 : kappamax = kappa;
1961 3783543 : if (incgram)
1962 : {
1963 15878950 : for (i=zeros+1; i<=kappa; i++)
1964 12293935 : gmael(G,kappa,i) = ZV_dotproduct(gel(B,kappa), gel(B,i));
1965 3585015 : maxG = kappamax;
1966 : }
1967 : }
1968 : /* Step3: Call to the Babai algorithm, mu,r,s updated in place */
1969 4575946 : if (Babai_dpe(av, kappa, &G,&B,&U, mu,r,s, alpha[kappa], zeros, maxG, &eta))
1970 0 : { *pG = incgram? NULL: G; *pB = B; *pU = U; return -1; }
1971 9054066 : if ((keepfirst && kappa == 2) ||
1972 4477872 : dpe_cmpmul(Dmael(r,kappa-1,kappa-1), &delta, Del(s,kappa-1)) <= 0)
1973 : { /* Step4: Success of Lovasz's condition */
1974 4251330 : alpha[kappa] = kappa;
1975 4251330 : dpe_submulz(Del(s,kappa-1), Dmael(mu,kappa,kappa-1), Dmael(r,kappa,kappa-1), Dmael(r,kappa,kappa));
1976 4251336 : continue;
1977 : }
1978 : /* Step5: Find the right insertion index kappa, kappa2 = initial kappa */
1979 324864 : if (DEBUGLEVEL>=4 && kappa==kappamax && Del(s,kappa-1)->d)
1980 0 : if (++cnt > 20) { cnt = 0; err_printf("(%ld) ", Del(s,1)->e-1); }
1981 324864 : kappa2 = kappa;
1982 : do {
1983 827802 : kappa--;
1984 827802 : if (kappa < zeros+2 + (keepfirst ? 1: 0)) break;
1985 706609 : } while (dpe_cmpmul(Dmael(r,kappa-1,kappa-1), &delta, Del(s,kappa-1)) >= 0);
1986 324864 : update_alpha(alpha, kappa, kappa2, kappamax);
1987 :
1988 : /* Step6: Update the mu's and r's */
1989 324865 : dperotate(mu, kappa2, kappa);
1990 324865 : dperotate(r, kappa2, kappa);
1991 324865 : affdpe(Del(s,kappa), Dmael(r,kappa,kappa));
1992 :
1993 : /* Step7: Update G, B, U */
1994 324865 : if (U) rotate(U, kappa2, kappa);
1995 324865 : if (B) rotate(B, kappa2, kappa);
1996 324865 : rotateG(G,kappa2,kappa, maxG, Gtmp);
1997 :
1998 : /* Step8: Prepare the next loop iteration */
1999 324865 : if (kappa == zeros+1 && !signe(gmael(G,kappa,kappa)))
2000 : {
2001 35161 : zeros++; kappa++;
2002 35161 : affidpe(gmael(G,kappa,kappa), Dmael(r,kappa,kappa));
2003 : }
2004 : }
2005 1564121 : if (pr) *pr = dpeM_diagonal_shallow(r,d);
2006 1564121 : *pG = G; *pB = B; *pU = U; return zeros; /* success */
2007 : }
2008 :
2009 :
2010 : /************************** PROVED version (t_INT) *************************/
2011 :
2012 : /* Babai's Nearest Plane algorithm (iterative).
2013 : * Size-reduces b_kappa using mu_{i,j} and r_{i,j} for j<=i <kappa
2014 : * Update B[,kappa]; compute mu_{kappa,j}, r_{kappa,j} for j<=kappa and s[kappa]
2015 : * mu, r, s updated in place (affrr). Return 1 on failure, else 0. */
2016 : static int
2017 0 : Babai(pari_sp av, long kappa, GEN *pG, GEN *pB, GEN *pU, GEN mu, GEN r, GEN s,
2018 : long a, long zeros, long maxG, GEN eta, long prec)
2019 : {
2020 0 : GEN G = *pG, B = *pB, U = *pU, ztmp;
2021 0 : long k, aa = a > zeros? a: zeros+1;
2022 0 : const long n = B? nbrows(B): 0, d = U ? lg(U)-1: 0, bit = prec2nbits(prec);
2023 0 : long emaxmu = EX0, emax2mu = EX0;
2024 : /* N.B: we set d = 0 (resp. n = 0) to avoid updating U (resp. B) */
2025 :
2026 0 : for (;;) {
2027 0 : int go_on = 0;
2028 0 : long i, j, emax3mu = emax2mu;
2029 :
2030 0 : if (gc_needed(av,2))
2031 : {
2032 0 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[1], a=%ld", aa);
2033 0 : gc_lll(av,3,&G,&B,&U);
2034 : }
2035 : /* Step2: compute the GSO for stage kappa */
2036 0 : emax2mu = emaxmu; emaxmu = EX0;
2037 0 : for (j=aa; j<kappa; j++)
2038 : {
2039 0 : pari_sp btop = avma;
2040 0 : GEN g = gmael(G,kappa,j);
2041 0 : for (k = zeros+1; k < j; k++)
2042 0 : g = mpsub(g, mulrr(gmael(mu,j,k), gmael(r,kappa,k)));
2043 0 : mpaff(g, gmael(r,kappa,j));
2044 0 : affrr(divrr(gmael(r,kappa,j), gmael(r,j,j)), gmael(mu,kappa,j));
2045 0 : emaxmu = maxss(emaxmu, expo(gmael(mu,kappa,j)));
2046 0 : set_avma(btop);
2047 : }
2048 0 : if (emax3mu != EX0 && emax3mu <= emax2mu + 5) /* precision too low */
2049 0 : { *pG = G; *pB = B; *pU = U; return 1; }
2050 :
2051 0 : for (j=kappa-1; j>zeros; j--)
2052 0 : if (abscmprr(gmael(mu,kappa,j), eta) > 0) { go_on = 1; break; }
2053 :
2054 : /* Step3--5: compute the X_j's */
2055 0 : if (go_on)
2056 0 : for (j=kappa-1; j>zeros; j--)
2057 : {
2058 : pari_sp btop;
2059 0 : GEN tmp = gmael(mu,kappa,j);
2060 0 : if (absrsmall(tmp)) continue; /* size-reduced */
2061 :
2062 0 : if (gc_needed(av,2))
2063 : {
2064 0 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[2], a=%ld, j=%ld", aa,j);
2065 0 : gc_lll(av,3,&G,&B,&U);
2066 : }
2067 0 : btop = avma;
2068 : /* we consider separately the case |X| = 1 */
2069 0 : if (absrsmall2(tmp))
2070 : {
2071 0 : if (signe(tmp) > 0) { /* in this case, X = 1 */
2072 0 : for (k=zeros+1; k<j; k++)
2073 0 : affrr(subrr(gmael(mu,kappa,k), gmael(mu,j,k)), gmael(mu,kappa,k));
2074 0 : set_avma(btop);
2075 0 : for (i=1; i<=n; i++)
2076 0 : gmael(B,kappa,i) = subii(gmael(B,kappa,i), gmael(B,j,i));
2077 0 : for (i=1; i<=d; i++)
2078 0 : gmael(U,kappa,i) = subii(gmael(U,kappa,i), gmael(U,j,i));
2079 0 : btop = avma;
2080 0 : ztmp = subii(gmael(G,j,j), shifti(gmael(G,kappa,j), 1));
2081 0 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
2082 0 : gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
2083 0 : for (i=1; i<=j; i++)
2084 0 : gmael(G,kappa,i) = subii(gmael(G,kappa,i), gmael(G,j,i));
2085 0 : for (i=j+1; i<kappa; i++)
2086 0 : gmael(G,kappa,i) = subii(gmael(G,kappa,i), gmael(G,i,j));
2087 0 : for (i=kappa+1; i<=maxG; i++)
2088 0 : gmael(G,i,kappa) = subii(gmael(G,i,kappa), gmael(G,i,j));
2089 : } else { /* otherwise X = -1 */
2090 0 : for (k=zeros+1; k<j; k++)
2091 0 : affrr(addrr(gmael(mu,kappa,k), gmael(mu,j,k)), gmael(mu,kappa,k));
2092 0 : set_avma(btop);
2093 0 : for (i=1; i<=n; i++)
2094 0 : gmael(B,kappa,i) = addii(gmael(B,kappa,i),gmael(B,j,i));
2095 0 : for (i=1; i<=d; i++)
2096 0 : gmael(U,kappa,i) = addii(gmael(U,kappa,i),gmael(U,j,i));
2097 0 : btop = avma;
2098 0 : ztmp = addii(gmael(G,j,j), shifti(gmael(G,kappa,j), 1));
2099 0 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
2100 0 : gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
2101 0 : for (i=1; i<=j; i++)
2102 0 : gmael(G,kappa,i) = addii(gmael(G,kappa,i), gmael(G,j,i));
2103 0 : for (i=j+1; i<kappa; i++)
2104 0 : gmael(G,kappa,i) = addii(gmael(G,kappa,i), gmael(G,i,j));
2105 0 : for (i=kappa+1; i<=maxG; i++)
2106 0 : gmael(G,i,kappa) = addii(gmael(G,i,kappa), gmael(G,i,j));
2107 : }
2108 0 : continue;
2109 : }
2110 : /* we have |X| >= 2 */
2111 0 : if (expo(tmp) < BITS_IN_LONG)
2112 : {
2113 0 : ulong xx = roundr_safe(tmp)[2]; /* X fits in an ulong */
2114 0 : if (signe(tmp) > 0) /* = xx */
2115 : {
2116 0 : for (k=zeros+1; k<j; k++)
2117 0 : affrr(subrr(gmael(mu,kappa,k), mulur(xx, gmael(mu,j,k))),
2118 0 : gmael(mu,kappa,k));
2119 0 : set_avma(btop);
2120 0 : for (i=1; i<=n; i++)
2121 0 : gmael(B,kappa,i) = submuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
2122 0 : for (i=1; i<=d; i++)
2123 0 : gmael(U,kappa,i) = submuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
2124 0 : btop = avma;
2125 0 : ztmp = submuliu2n(mulii(gmael(G,j,j), sqru(xx)), gmael(G,kappa,j), xx, 1);
2126 0 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
2127 0 : gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
2128 0 : for (i=1; i<=j; i++)
2129 0 : gmael(G,kappa,i) = submuliu_inplace(gmael(G,kappa,i), gmael(G,j,i), xx);
2130 0 : for (i=j+1; i<kappa; i++)
2131 0 : gmael(G,kappa,i) = submuliu_inplace(gmael(G,kappa,i), gmael(G,i,j), xx);
2132 0 : for (i=kappa+1; i<=maxG; i++)
2133 0 : gmael(G,i,kappa) = submuliu_inplace(gmael(G,i,kappa), gmael(G,i,j), xx);
2134 : }
2135 : else /* = -xx */
2136 : {
2137 0 : for (k=zeros+1; k<j; k++)
2138 0 : affrr(addrr(gmael(mu,kappa,k), mulur(xx, gmael(mu,j,k))),
2139 0 : gmael(mu,kappa,k));
2140 0 : set_avma(btop);
2141 0 : for (i=1; i<=n; i++)
2142 0 : gmael(B,kappa,i) = addmuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
2143 0 : for (i=1; i<=d; i++)
2144 0 : gmael(U,kappa,i) = addmuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
2145 0 : btop = avma;
2146 0 : ztmp = addmuliu2n(mulii(gmael(G,j,j), sqru(xx)), gmael(G,kappa,j), xx, 1);
2147 0 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
2148 0 : gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
2149 0 : for (i=1; i<=j; i++)
2150 0 : gmael(G,kappa,i) = addmuliu_inplace(gmael(G,kappa,i), gmael(G,j,i), xx);
2151 0 : for (i=j+1; i<kappa; i++)
2152 0 : gmael(G,kappa,i) = addmuliu_inplace(gmael(G,kappa,i), gmael(G,i,j), xx);
2153 0 : for (i=kappa+1; i<=maxG; i++)
2154 0 : gmael(G,i,kappa) = addmuliu_inplace(gmael(G,i,kappa), gmael(G,i,j), xx);
2155 : }
2156 : }
2157 : else
2158 : {
2159 : long e;
2160 0 : GEN X = truncexpo(tmp, bit, &e); /* tmp ~ X * 2^e */
2161 0 : btop = avma;
2162 0 : for (k=zeros+1; k<j; k++)
2163 : {
2164 0 : GEN x = mulir(X, gmael(mu,j,k));
2165 0 : if (e) shiftr_inplace(x, e);
2166 0 : affrr(subrr(gmael(mu,kappa,k), x), gmael(mu,kappa,k));
2167 : }
2168 0 : set_avma(btop);
2169 0 : for (i=1; i<=n; i++)
2170 0 : gmael(B,kappa,i) = submulshift(gmael(B,kappa,i), gmael(B,j,i), X, e);
2171 0 : for (i=1; i<=d; i++)
2172 0 : gmael(U,kappa,i) = submulshift(gmael(U,kappa,i), gmael(U,j,i), X, e);
2173 0 : btop = avma;
2174 0 : ztmp = submulshift(mulshift(gmael(G,j,j), sqri(X), 2*e),
2175 0 : gmael(G,kappa,j), X, e+1);
2176 0 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
2177 0 : gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
2178 0 : for (i=1; i<=j; i++)
2179 0 : gmael(G,kappa,i) = submulshift(gmael(G,kappa,i), gmael(G,j,i), X, e);
2180 0 : for ( ; i<kappa; i++)
2181 0 : gmael(G,kappa,i) = submulshift(gmael(G,kappa,i), gmael(G,i,j), X, e);
2182 0 : for (i=kappa+1; i<=maxG; i++)
2183 0 : gmael(G,i,kappa) = submulshift(gmael(G,i,kappa), gmael(G,i,j), X, e);
2184 : }
2185 : }
2186 0 : if (!go_on) break; /* Anything happened? */
2187 0 : aa = zeros+1;
2188 : }
2189 :
2190 0 : affir(gmael(G,kappa,kappa), gel(s,zeros+1));
2191 : /* the last s[kappa-1]=r[kappa][kappa] is computed only if kappa increases */
2192 0 : av = avma;
2193 0 : for (k=zeros+1; k<=kappa-2; k++)
2194 0 : affrr(subrr(gel(s,k), mulrr(gmael(mu,kappa,k), gmael(r,kappa,k))),
2195 0 : gel(s,k+1));
2196 0 : *pG = G; *pB = B; *pU = U; return gc_bool(av, 0);
2197 : }
2198 :
2199 : /* G integral Gram matrix, LLL-reduces (G,B,U) in place [apply base change
2200 : * transforms to B and U]. If (keepfirst), never swap with first vector.
2201 : * If G = NULL, we compute the Gram matrix incrementally.
2202 : * Return -1 on failure, else zeros = dim Kernel (>= 0) */
2203 : static long
2204 0 : fplll(GEN *pG, GEN *pB, GEN *pU, GEN *pr, double DELTA, double ETA,
2205 : long keepfirst, long prec)
2206 : {
2207 : pari_sp av, av2;
2208 0 : GEN mu, r, s, tmp, Gtmp, alpha, G = *pG, B = *pB, U = *pU;
2209 0 : GEN delta = dbltor(DELTA), eta = dbltor(ETA);
2210 0 : long d, maxG, kappa, kappa2, i, j, zeros, kappamax, incgram = !G, cnt = 0;
2211 :
2212 0 : if (incgram)
2213 : { /* incremental Gram matrix */
2214 0 : maxG = 2; d = lg(B)-1;
2215 0 : G = zeromatcopy(d, d);
2216 : }
2217 : else
2218 0 : maxG = d = lg(G)-1;
2219 :
2220 0 : mu = cgetg(d+1, t_MAT);
2221 0 : r = cgetg(d+1, t_MAT);
2222 0 : s = cgetg(d+1, t_VEC);
2223 0 : for (j = 1; j <= d; j++)
2224 : {
2225 0 : GEN M = cgetg(d+1, t_COL), R = cgetg(d+1, t_COL);
2226 0 : gel(mu,j)= M;
2227 0 : gel(r,j) = R;
2228 0 : gel(s,j) = cgetr(prec);
2229 0 : for (i = 1; i <= d; i++)
2230 : {
2231 0 : gel(R,i) = cgetr(prec);
2232 0 : gel(M,i) = cgetr(prec);
2233 : }
2234 : }
2235 0 : Gtmp = cgetg(d+1, t_VEC);
2236 0 : alpha = cgetg(d+1, t_VECSMALL);
2237 0 : av = avma;
2238 :
2239 : /* Step2: Initializing the main loop */
2240 0 : kappamax = 1;
2241 0 : i = 1;
2242 : do {
2243 0 : if (incgram) gmael(G,i,i) = ZV_dotsquare(gel(B,i));
2244 0 : affir(gmael(G,i,i), gmael(r,i,i));
2245 0 : } while (!signe(gmael(G,i,i)) && ++i <= d);
2246 0 : zeros = i-1; /* all basis vectors b_i with i <= zeros are zero vectors */
2247 0 : kappa = i;
2248 0 : for (i=zeros+1; i<=d; i++) alpha[i]=1;
2249 :
2250 0 : while (++kappa <= d)
2251 : {
2252 0 : if (kappa > kappamax)
2253 : {
2254 0 : if (DEBUGLEVEL>=4) err_printf("K%ld ",kappa);
2255 0 : kappamax = kappa;
2256 0 : if (incgram)
2257 : {
2258 0 : for (i=zeros+1; i<=kappa; i++)
2259 0 : gmael(G,kappa,i) = ZV_dotproduct(gel(B,kappa), gel(B,i));
2260 0 : maxG = kappamax;
2261 : }
2262 : }
2263 : /* Step3: Call to the Babai algorithm, mu,r,s updated in place */
2264 0 : if (Babai(av, kappa, &G,&B,&U, mu,r,s, alpha[kappa], zeros, maxG, eta, prec))
2265 0 : { *pG = incgram? NULL: G; *pB = B; *pU = U; return -1; }
2266 0 : av2 = avma;
2267 0 : if ((keepfirst && kappa == 2) ||
2268 0 : cmprr(mulrr(gmael(r,kappa-1,kappa-1), delta), gel(s,kappa-1)) <= 0)
2269 : { /* Step4: Success of Lovasz's condition */
2270 0 : alpha[kappa] = kappa;
2271 0 : tmp = mulrr(gmael(mu,kappa,kappa-1), gmael(r,kappa,kappa-1));
2272 0 : affrr(subrr(gel(s,kappa-1), tmp), gmael(r,kappa,kappa));
2273 0 : set_avma(av2); continue;
2274 : }
2275 : /* Step5: Find the right insertion index kappa, kappa2 = initial kappa */
2276 0 : if (DEBUGLEVEL>=4 && kappa==kappamax && signe(gel(s,kappa-1)))
2277 0 : if (++cnt > 20) { cnt = 0; err_printf("(%ld) ", expo(gel(s,1))); }
2278 0 : kappa2 = kappa;
2279 : do {
2280 0 : kappa--;
2281 0 : if (kappa < zeros+2 + (keepfirst ? 1: 0)) break;
2282 0 : tmp = mulrr(gmael(r,kappa-1,kappa-1), delta);
2283 0 : } while (cmprr(gel(s,kappa-1), tmp) <= 0);
2284 0 : set_avma(av2);
2285 0 : update_alpha(alpha, kappa, kappa2, kappamax);
2286 :
2287 : /* Step6: Update the mu's and r's */
2288 0 : rotate(mu, kappa2, kappa);
2289 0 : rotate(r, kappa2, kappa);
2290 0 : affrr(gel(s,kappa), gmael(r,kappa,kappa));
2291 :
2292 : /* Step7: Update G, B, U */
2293 0 : if (U) rotate(U, kappa2, kappa);
2294 0 : if (B) rotate(B, kappa2, kappa);
2295 0 : rotateG(G,kappa2,kappa, maxG, Gtmp);
2296 :
2297 : /* Step8: Prepare the next loop iteration */
2298 0 : if (kappa == zeros+1 && !signe(gmael(G,kappa,kappa)))
2299 : {
2300 0 : zeros++; kappa++;
2301 0 : affir(gmael(G,kappa,kappa), gmael(r,kappa,kappa));
2302 : }
2303 : }
2304 0 : if (pr) *pr = RgM_diagonal_shallow(r);
2305 0 : *pG = G; *pB = B; *pU = U; return zeros; /* success */
2306 : }
2307 :
2308 : /* do not support LLL_KER, LLL_ALL, LLL_KEEP_FIRST */
2309 : static GEN
2310 4630824 : ZM2_lll_norms(GEN x, long flag, GEN *pN)
2311 : {
2312 : GEN a,b,c,d;
2313 : GEN G, U;
2314 4630824 : if (flag & LLL_GRAM)
2315 7271 : G = x;
2316 : else
2317 4623553 : G = gram_matrix(x);
2318 4630811 : a = gcoeff(G,1,1); b = shifti(gcoeff(G,1,2),1); c = gcoeff(G,2,2);
2319 4630793 : d = qfb_disc3(a,b,c);
2320 4630772 : if (signe(d)>=0) return NULL;
2321 4630387 : G = redimagsl2(mkqfb(a,b,c,d),&U);
2322 4630441 : if (pN) (void) RgM_gram_schmidt(G, pN);
2323 4630441 : if (flag & LLL_INPLACE) return ZM2_mul(x,U);
2324 4630441 : return U;
2325 : }
2326 :
2327 : static void
2328 616411 : fplll_flatter(GEN *pG, GEN *pB, GEN *pU, long rank, long flag)
2329 : {
2330 616411 : if (!*pG)
2331 : {
2332 615436 : GEN T = ZM_flatter_rank(*pB, rank, flag);
2333 615438 : if (*pU)
2334 : {
2335 601434 : *pU = ZM_mul(*pU, T);
2336 601434 : *pB = ZM_mul(*pB, T);
2337 14004 : } else *pB = T;
2338 : } else
2339 : {
2340 975 : GEN T, G = *pG;
2341 975 : long i, j, l = lg(G);
2342 7291 : for (i = 1; i < l; i++)
2343 43593 : for(j = 1; j < i; j++)
2344 37277 : gmael(G,j,i) = gmael(G,i,j);
2345 975 : T = ZM_flattergram_rank(G, rank, flag);
2346 975 : if (*pU) *pU = ZM_mul(*pU, T);
2347 975 : *pG = ZM_mul(shallowtrans(T), ZM_mul(*pG,T));
2348 : }
2349 616413 : }
2350 :
2351 : static GEN
2352 1081143 : get_gramschmidt(GEN M, long rank)
2353 : {
2354 : GEN B, Q, L;
2355 1081143 : long r = lg(M)-1, bitprec = 3*r + 30;
2356 1081143 : long prec = nbits2prec64(bitprec);
2357 1081143 : if (rank < r) M = vconcat(gshift(M,1), matid(r));
2358 1081143 : if (!QR_init(RgM_gtofp(M, prec), &B, &Q, &L, prec) || !gsisinv(L)) return NULL;
2359 466755 : return L;
2360 : }
2361 :
2362 : static GEN
2363 43854 : get_gaussred(GEN M, long rank)
2364 : {
2365 43854 : pari_sp ltop = avma;
2366 43854 : long r = lg(M)-1, bitprec = 3*r + 30, prec = nbits2prec64(bitprec);
2367 : GEN R;
2368 43854 : if (rank < r) M = RgM_Rg_add(gshift(M, 1), gen_1);
2369 43854 : R = RgM_Cholesky(RgM_gtofp(M, prec), prec);
2370 43853 : if (!R) return NULL;
2371 42878 : return gerepilecopy(ltop, R);
2372 : }
2373 :
2374 : /* Assume x a ZM, if pN != NULL, set it to Gram-Schmidt (squared) norms
2375 : * The following modes are supported:
2376 : * - flag & LLL_INPLACE: x a lattice basis, return x*U
2377 : * - flag & LLL_GRAM: x a Gram matrix / else x a lattice basis; return
2378 : * LLL base change matrix U [LLL_IM]
2379 : * kernel basis [LLL_KER, nonreduced]
2380 : * both [LLL_ALL] */
2381 : GEN
2382 6877410 : ZM_lll_norms(GEN x, double DELTA, long flag, GEN *pN)
2383 : {
2384 6877410 : pari_sp av = avma;
2385 6877410 : const double ETA = 0.51;
2386 6877410 : const long keepfirst = flag & LLL_KEEP_FIRST;
2387 6877410 : long p, zeros = -1, n = lg(x)-1, is_upper, is_lower, useflatter = 0, rank;
2388 6877410 : GEN G, B, U, L = NULL;
2389 : pari_timer T;
2390 6877410 : long thre[]={31783,34393,20894,22525,13533,1928,672,671,
2391 : 422,506,315,313,222,205,167,154,139,138,
2392 : 110,120,98,94,81,75,74,64,74,74,
2393 : 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,
2394 : 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};
2395 6877410 : long thsn[]={23280,30486,50077,44136,78724,15690,1801,1611,
2396 : 981,1359,978,1042,815,866,788,775,726,712,
2397 : 626,613,548,564,474,481,504,447,453,508,
2398 : 705,794,1008,946,767,898,886,763,842,757,
2399 : 725,774,639,655,705,627,635,704,511,613,
2400 : 583,595,568,640,541,640,567,540,577,584,
2401 : 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};
2402 6877410 : if (n <= 1) return lll_trivial(x, flag);
2403 6768473 : if (nbrows(x)==0)
2404 : {
2405 15197 : if (flag & LLL_KER) return matid(n);
2406 15197 : if (flag & (LLL_INPLACE|LLL_IM)) return cgetg(1,t_MAT);
2407 0 : retmkvec2(matid(n), cgetg(1,t_MAT));
2408 : }
2409 6753434 : if (n==2 && nbrows(x)==2 && (flag&LLL_IM) && !keepfirst)
2410 : {
2411 4630824 : U = ZM2_lll_norms(x, flag, pN);
2412 4630825 : if (U) return U;
2413 : }
2414 2122986 : if (flag & LLL_GRAM)
2415 59737 : { G = x; B = NULL; U = matid(n); is_upper = 0; is_lower = 0; }
2416 : else
2417 : {
2418 2063249 : G = NULL; B = x; U = (flag & LLL_INPLACE)? NULL: matid(n);
2419 2063258 : is_upper = (flag & LLL_UPPER) || ZM_is_upper(B);
2420 2063256 : is_lower = !B || is_upper || keepfirst ? 0: ZM_is_lower(B);
2421 2063253 : if (is_lower) L = RgM_flip(B);
2422 : }
2423 2122989 : rank = (flag&LLL_NOFLATTER) ? 0: ZM_rank(x);
2424 2122970 : if (n > 2 && !(flag&LLL_NOFLATTER))
2425 1729395 : {
2426 1685541 : GEN R = B ? (is_upper ? B : (is_lower ? L : get_gramschmidt(B, rank)))
2427 3414919 : : get_gaussred(G, rank);
2428 1729384 : if (R)
2429 : {
2430 1114032 : long spr = spread(R), sz = mpexpo(gsupnorm(R, DEFAULTPREC)), thr;
2431 1114043 : if (DEBUGLEVEL>=5) err_printf("LLL: dim %ld, size %ld, spread %ld\n",n, sz, spr);
2432 1114043 : if ((is_upper && ZM_is_knapsack(B)) || (is_lower && ZM_is_knapsack(L)))
2433 92314 : thr = thsn[minss(n-3,numberof(thsn)-1)];
2434 : else
2435 : {
2436 1021729 : thr = thre[minss(n-3,numberof(thre)-1)];
2437 1021729 : if (n>=10) sz = spr;
2438 : }
2439 1114043 : useflatter = sz >= thr;
2440 : } else
2441 615352 : useflatter = 1;
2442 393578 : } else useflatter = 0;
2443 2122973 : if(DEBUGLEVEL>=4) timer_start(&T);
2444 2122973 : if (useflatter)
2445 : {
2446 616411 : if (is_lower)
2447 : {
2448 0 : fplll_flatter(&G, &L, &U, rank, flag | LLL_UPPER);
2449 0 : B = RgM_flop(L);
2450 0 : if (U) U = RgM_flop(U);
2451 : }
2452 : else
2453 616411 : fplll_flatter(&G, &B, &U, rank, flag | (is_upper? LLL_UPPER:0));
2454 616413 : if (DEBUGLEVEL>=4 && !(flag & LLL_NOCERTIFY))
2455 0 : timer_printf(&T, "FLATTER");
2456 : }
2457 2122974 : if (!(flag & LLL_GRAM))
2458 : {
2459 : long t;
2460 2063238 : B = gcopy(B);
2461 2063250 : if(DEBUGLEVEL>=4)
2462 0 : err_printf("Entering L^2 (double): dim %ld, LLL-parameters (%.3f,%.3f)\n",
2463 : n, DELTA,ETA);
2464 2063250 : zeros = fplll_fast(&B, &U, DELTA, ETA, keepfirst);
2465 2063260 : if (DEBUGLEVEL>=4) timer_printf(&T, zeros < 0? "LLL (failed)": "LLL");
2466 2067375 : for (p = DEFAULTPREC, t = 0; zeros < 0 && t < 1 ; p += EXTRAPREC64, t++)
2467 : {
2468 4115 : if (DEBUGLEVEL>=4)
2469 0 : err_printf("Entering L^2 (heuristic): LLL-parameters (%.3f,%.3f), prec = %d/%d\n", DELTA, ETA, p, p);
2470 4115 : zeros = fplll_heuristic(&B, &U, DELTA, ETA, keepfirst, p, p);
2471 4115 : gc_lll(av, 2, &B, &U);
2472 4115 : if (DEBUGLEVEL>=4) timer_printf(&T, zeros < 0? "LLL (failed)": "LLL");
2473 : }
2474 : } else
2475 59736 : G = gcopy(G);
2476 2122996 : if (zeros < 0 || !(flag & LLL_NOCERTIFY))
2477 : {
2478 1564138 : if(DEBUGLEVEL>=4)
2479 0 : err_printf("Entering L^2 (dpe): LLL-parameters (%.3f,%.3f)\n", DELTA,ETA);
2480 1564138 : zeros = fplll_dpe(&G, &B, &U, pN, DELTA, ETA, keepfirst);
2481 1564120 : if (DEBUGLEVEL>=4) timer_printf(&T, zeros < 0? "LLL (failed)": "LLL");
2482 1564125 : if (zeros < 0)
2483 0 : for (p = DEFAULTPREC;; p += EXTRAPREC64)
2484 : {
2485 0 : if (DEBUGLEVEL>=4)
2486 0 : err_printf("Entering L^2: LLL-parameters (%.3f,%.3f), prec = %d\n",
2487 : DELTA,ETA, p);
2488 0 : zeros = fplll(&G, &B, &U, pN, DELTA, ETA, keepfirst, p);
2489 0 : if (DEBUGLEVEL>=4) timer_printf(&T, zeros < 0? "LLL (failed)": "LLL");
2490 0 : if (zeros >= 0) break;
2491 0 : gc_lll(av, 3, &G, &B, &U);
2492 : }
2493 : }
2494 2122983 : return lll_finish(U? U: B, zeros, flag);
2495 : }
2496 :
2497 : /********************************************************************/
2498 : /** **/
2499 : /** LLL OVER K[X] **/
2500 : /** **/
2501 : /********************************************************************/
2502 : static int
2503 504 : pslg(GEN x)
2504 : {
2505 : long tx;
2506 504 : if (gequal0(x)) return 2;
2507 448 : tx = typ(x); return is_scalar_t(tx)? 3: lg(x);
2508 : }
2509 :
2510 : static int
2511 196 : REDgen(long k, long l, GEN h, GEN L, GEN B)
2512 : {
2513 196 : GEN q, u = gcoeff(L,k,l);
2514 : long i;
2515 :
2516 196 : if (pslg(u) < pslg(B)) return 0;
2517 :
2518 140 : q = gneg(gdeuc(u,B));
2519 140 : gel(h,k) = gadd(gel(h,k), gmul(q,gel(h,l)));
2520 140 : for (i=1; i<l; i++) gcoeff(L,k,i) = gadd(gcoeff(L,k,i), gmul(q,gcoeff(L,l,i)));
2521 140 : gcoeff(L,k,l) = gadd(gcoeff(L,k,l), gmul(q,B)); return 1;
2522 : }
2523 :
2524 : static int
2525 196 : do_SWAPgen(GEN h, GEN L, GEN B, long k, GEN fl, int *flc)
2526 : {
2527 : GEN p1, la, la2, Bk;
2528 : long ps1, ps2, i, j, lx;
2529 :
2530 196 : if (!fl[k-1]) return 0;
2531 :
2532 140 : la = gcoeff(L,k,k-1); la2 = gsqr(la);
2533 140 : Bk = gel(B,k);
2534 140 : if (fl[k])
2535 : {
2536 56 : GEN q = gadd(la2, gmul(gel(B,k-1),gel(B,k+1)));
2537 56 : ps1 = pslg(gsqr(Bk));
2538 56 : ps2 = pslg(q);
2539 56 : if (ps1 <= ps2 && (ps1 < ps2 || !*flc)) return 0;
2540 28 : *flc = (ps1 != ps2);
2541 28 : gel(B,k) = gdiv(q, Bk);
2542 : }
2543 :
2544 112 : swap(gel(h,k-1), gel(h,k)); lx = lg(L);
2545 112 : for (j=1; j<k-1; j++) swap(gcoeff(L,k-1,j), gcoeff(L,k,j));
2546 112 : if (fl[k])
2547 : {
2548 28 : for (i=k+1; i<lx; i++)
2549 : {
2550 0 : GEN t = gcoeff(L,i,k);
2551 0 : p1 = gsub(gmul(gel(B,k+1),gcoeff(L,i,k-1)), gmul(la,t));
2552 0 : gcoeff(L,i,k) = gdiv(p1, Bk);
2553 0 : p1 = gadd(gmul(la,gcoeff(L,i,k-1)), gmul(gel(B,k-1),t));
2554 0 : gcoeff(L,i,k-1) = gdiv(p1, Bk);
2555 : }
2556 : }
2557 84 : else if (!gequal0(la))
2558 : {
2559 28 : p1 = gdiv(la2, Bk);
2560 28 : gel(B,k+1) = gel(B,k) = p1;
2561 28 : for (i=k+2; i<=lx; i++) gel(B,i) = gdiv(gmul(p1,gel(B,i)),Bk);
2562 28 : for (i=k+1; i<lx; i++)
2563 0 : gcoeff(L,i,k-1) = gdiv(gmul(la,gcoeff(L,i,k-1)), Bk);
2564 28 : for (j=k+1; j<lx-1; j++)
2565 0 : for (i=j+1; i<lx; i++)
2566 0 : gcoeff(L,i,j) = gdiv(gmul(p1,gcoeff(L,i,j)), Bk);
2567 : }
2568 : else
2569 : {
2570 56 : gcoeff(L,k,k-1) = gen_0;
2571 56 : for (i=k+1; i<lx; i++)
2572 : {
2573 0 : gcoeff(L,i,k) = gcoeff(L,i,k-1);
2574 0 : gcoeff(L,i,k-1) = gen_0;
2575 : }
2576 56 : gel(B,k) = gel(B,k-1); fl[k] = 1; fl[k-1] = 0;
2577 : }
2578 112 : return 1;
2579 : }
2580 :
2581 : static void
2582 168 : incrementalGSgen(GEN x, GEN L, GEN B, long k, GEN fl)
2583 : {
2584 168 : GEN u = NULL; /* gcc -Wall */
2585 : long i, j;
2586 420 : for (j = 1; j <= k; j++)
2587 252 : if (j==k || fl[j])
2588 : {
2589 252 : u = gcoeff(x,k,j);
2590 252 : if (!is_extscalar_t(typ(u))) pari_err_TYPE("incrementalGSgen",u);
2591 336 : for (i=1; i<j; i++)
2592 84 : if (fl[i])
2593 : {
2594 84 : u = gsub(gmul(gel(B,i+1),u), gmul(gcoeff(L,k,i),gcoeff(L,j,i)));
2595 84 : u = gdiv(u, gel(B,i));
2596 : }
2597 252 : gcoeff(L,k,j) = u;
2598 : }
2599 168 : if (gequal0(u)) gel(B,k+1) = gel(B,k);
2600 : else
2601 : {
2602 112 : gel(B,k+1) = gcoeff(L,k,k); gcoeff(L,k,k) = gen_1; fl[k] = 1;
2603 : }
2604 168 : }
2605 :
2606 : static GEN
2607 168 : lllgramallgen(GEN x, long flag)
2608 : {
2609 168 : long lx = lg(x), i, j, k, l, n;
2610 : pari_sp av;
2611 : GEN B, L, h, fl;
2612 : int flc;
2613 :
2614 168 : n = lx-1; if (n<=1) return lll_trivial(x,flag);
2615 84 : if (lgcols(x) != lx) pari_err_DIM("lllgramallgen");
2616 :
2617 84 : fl = cgetg(lx, t_VECSMALL);
2618 :
2619 84 : av = avma;
2620 84 : B = scalarcol_shallow(gen_1, lx);
2621 84 : L = cgetg(lx,t_MAT);
2622 252 : for (j=1; j<lx; j++) { gel(L,j) = zerocol(n); fl[j] = 0; }
2623 :
2624 84 : h = matid(n);
2625 252 : for (i=1; i<lx; i++)
2626 168 : incrementalGSgen(x, L, B, i, fl);
2627 84 : flc = 0;
2628 84 : for(k=2;;)
2629 : {
2630 196 : if (REDgen(k, k-1, h, L, gel(B,k))) flc = 1;
2631 196 : if (do_SWAPgen(h, L, B, k, fl, &flc)) { if (k > 2) k--; }
2632 : else
2633 : {
2634 84 : for (l=k-2; l>=1; l--)
2635 0 : if (REDgen(k, l, h, L, gel(B,l+1))) flc = 1;
2636 84 : if (++k > n) break;
2637 : }
2638 112 : if (gc_needed(av,1))
2639 : {
2640 0 : if(DEBUGMEM>1) pari_warn(warnmem,"lllgramallgen");
2641 0 : gerepileall(av,3,&B,&L,&h);
2642 : }
2643 : }
2644 140 : k=1; while (k<lx && !fl[k]) k++;
2645 84 : return lll_finish(h,k-1,flag);
2646 : }
2647 :
2648 : static GEN
2649 168 : lllallgen(GEN x, long flag)
2650 : {
2651 168 : pari_sp av = avma;
2652 168 : if (!(flag & LLL_GRAM)) x = gram_matrix(x);
2653 84 : else if (!RgM_is_square_mat(x)) pari_err_DIM("qflllgram");
2654 168 : return gerepilecopy(av, lllgramallgen(x, flag));
2655 : }
2656 : GEN
2657 42 : lllgen(GEN x) { return lllallgen(x, LLL_IM); }
2658 : GEN
2659 42 : lllkerimgen(GEN x) { return lllallgen(x, LLL_ALL); }
2660 : GEN
2661 42 : lllgramgen(GEN x) { return lllallgen(x, LLL_IM|LLL_GRAM); }
2662 : GEN
2663 42 : lllgramkerimgen(GEN x) { return lllallgen(x, LLL_ALL|LLL_GRAM); }
2664 :
2665 : static GEN
2666 50782 : lllall(GEN x, long flag)
2667 50782 : { pari_sp av = avma; return gerepilecopy(av, ZM_lll(x, LLLDFT, flag)); }
2668 : GEN
2669 14651 : lllint(GEN x) { return lllall(x, LLL_IM); }
2670 : GEN
2671 35 : lllkerim(GEN x) { return lllall(x, LLL_ALL); }
2672 : GEN
2673 36054 : lllgramint(GEN x)
2674 36054 : { if (!RgM_is_square_mat(x)) pari_err_DIM("qflllgram");
2675 36054 : return lllall(x, LLL_IM | LLL_GRAM); }
2676 : GEN
2677 35 : lllgramkerim(GEN x)
2678 35 : { if (!RgM_is_square_mat(x)) pari_err_DIM("qflllgram");
2679 35 : return lllall(x, LLL_ALL | LLL_GRAM); }
2680 :
2681 : GEN
2682 5261131 : lllfp(GEN x, double D, long flag)
2683 : {
2684 5261131 : long n = lg(x)-1;
2685 5261131 : pari_sp av = avma;
2686 : GEN h;
2687 5261131 : if (n <= 1) return lll_trivial(x,flag);
2688 4601249 : if (flag & LLL_GRAM)
2689 : {
2690 9128 : if (!RgM_is_square_mat(x)) pari_err_DIM("qflllgram");
2691 9114 : if (isinexact(x))
2692 : {
2693 9051 : x = RgM_Cholesky(x, gprecision(x));
2694 9051 : if (!x) return NULL;
2695 9051 : flag &= ~LLL_GRAM;
2696 : }
2697 : }
2698 4601235 : h = ZM_lll(RgM_rescale_to_int(x), D, flag);
2699 4601211 : return gerepilecopy(av, h);
2700 : }
2701 :
2702 : GEN
2703 8949 : lllgram(GEN x) { return lllfp(x,LLLDFT,LLL_GRAM|LLL_IM); }
2704 : GEN
2705 1174642 : lll(GEN x) { return lllfp(x,LLLDFT,LLL_IM); }
2706 :
2707 : static GEN
2708 63 : qflllgram(GEN x)
2709 : {
2710 63 : GEN T = lllgram(x);
2711 42 : if (!T) pari_err_PREC("qflllgram");
2712 42 : return T;
2713 : }
2714 :
2715 : GEN
2716 301 : qflll0(GEN x, long flag)
2717 : {
2718 301 : if (typ(x) != t_MAT) pari_err_TYPE("qflll",x);
2719 301 : switch(flag)
2720 : {
2721 49 : case 0: return lll(x);
2722 63 : case 1: RgM_check_ZM(x,"qflll"); return lllint(x);
2723 49 : case 2: RgM_check_ZM(x,"qflll"); return lllintpartial(x);
2724 7 : case 3: RgM_check_ZM(x,"qflll"); return lllall(x, LLL_INPLACE);
2725 49 : case 4: RgM_check_ZM(x,"qflll"); return lllkerim(x);
2726 42 : case 5: return lllkerimgen(x);
2727 42 : case 8: return lllgen(x);
2728 0 : default: pari_err_FLAG("qflll");
2729 : }
2730 : return NULL; /* LCOV_EXCL_LINE */
2731 : }
2732 :
2733 : GEN
2734 245 : qflllgram0(GEN x, long flag)
2735 : {
2736 245 : if (typ(x) != t_MAT) pari_err_TYPE("qflllgram",x);
2737 245 : switch(flag)
2738 : {
2739 63 : case 0: return qflllgram(x);
2740 49 : case 1: RgM_check_ZM(x,"qflllgram"); return lllgramint(x);
2741 49 : case 4: RgM_check_ZM(x,"qflllgram"); return lllgramkerim(x);
2742 42 : case 5: return lllgramkerimgen(x);
2743 42 : case 8: return lllgramgen(x);
2744 0 : default: pari_err_FLAG("qflllgram");
2745 : }
2746 : return NULL; /* LCOV_EXCL_LINE */
2747 : }
2748 :
2749 : /********************************************************************/
2750 : /** **/
2751 : /** INTEGRAL KERNEL (LLL REDUCED) **/
2752 : /** **/
2753 : /********************************************************************/
2754 : static GEN
2755 70 : kerint0(GEN M)
2756 : {
2757 : /* return ZM_lll(M, LLLDFT, LLL_KER); */
2758 70 : GEN U, H = ZM_hnflll(M,&U,1);
2759 70 : long d = lg(M)-lg(H);
2760 70 : if (!d) return cgetg(1, t_MAT);
2761 70 : return ZM_lll(vecslice(U,1,d), LLLDFT, LLL_INPLACE);
2762 : }
2763 : GEN
2764 42 : kerint(GEN M)
2765 : {
2766 42 : pari_sp av = avma;
2767 42 : return gerepilecopy(av, kerint0(M));
2768 : }
2769 : /* OBSOLETE: use kerint */
2770 : GEN
2771 28 : matkerint0(GEN M, long flag)
2772 : {
2773 28 : pari_sp av = avma;
2774 28 : if (typ(M) != t_MAT) pari_err_TYPE("matkerint",M);
2775 28 : M = Q_primpart(M);
2776 28 : RgM_check_ZM(M, "kerint");
2777 28 : switch(flag)
2778 : {
2779 28 : case 0:
2780 28 : case 1: return gerepilecopy(av, kerint0(M));
2781 0 : default: pari_err_FLAG("matkerint");
2782 : }
2783 : return NULL; /* LCOV_EXCL_LINE */
2784 : }
|
