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 double
21 18781733 : pari_rint(double a)
22 : {
23 : #ifdef HAS_RINT
24 18781733 : return rint(a);
25 : #else
26 : const double two_to_52 = 4.5035996273704960e+15;
27 : double fa = fabs(a);
28 : double r = two_to_52 + fa;
29 : if (fa >= two_to_52) {
30 : r = a;
31 : } else {
32 : r = r - two_to_52;
33 : if (a < 0) r = -r;
34 : }
35 : return r;
36 : #endif
37 : }
38 :
39 : /* default quality ratio for LLL */
40 : static const double LLLDFT = 0.99;
41 :
42 : /* assume flag & (LLL_KER|LLL_IM|LLL_ALL). LLL_INPLACE implies LLL_IM */
43 : static GEN
44 192883 : lll_trivial(GEN x, long flag)
45 : {
46 192883 : if (lg(x) == 1)
47 : { /* dim x = 0 */
48 15401 : if (! (flag & LLL_ALL)) return cgetg(1,t_MAT);
49 28 : retmkvec2(cgetg(1,t_MAT), cgetg(1,t_MAT));
50 : }
51 : /* dim x = 1 */
52 177482 : if (gequal0(gel(x,1)))
53 : {
54 1386 : if (flag & LLL_KER) return matid(1);
55 1386 : if (flag & (LLL_IM|LLL_INPLACE)) return cgetg(1,t_MAT);
56 28 : retmkvec2(matid(1), cgetg(1,t_MAT));
57 : }
58 176094 : if (flag & LLL_INPLACE) return gcopy(x);
59 73005 : if (flag & LLL_KER) return cgetg(1,t_MAT);
60 73005 : if (flag & LLL_IM) return matid(1);
61 27 : retmkvec2(cgetg(1,t_MAT), (flag & LLL_GRAM)? gcopy(x): matid(1));
62 : }
63 :
64 : /* vecslice(x,#x-k,#x) in place. Works for t_MAT, t_VEC/t_COL */
65 : static GEN
66 2335986 : vectail_inplace(GEN x, long k)
67 : {
68 2335986 : if (!k) return x;
69 57258 : x[k] = ((ulong)x[0] & ~LGBITS) | evallg(lg(x) - k);
70 57262 : return x + k;
71 : }
72 :
73 : /* k = dim Kernel */
74 : static GEN
75 2409079 : lll_finish(GEN h, long k, long flag)
76 : {
77 : GEN g;
78 2409079 : if (!(flag & (LLL_IM|LLL_KER|LLL_ALL|LLL_INPLACE))) return h;
79 2336009 : if (flag & (LLL_IM|LLL_INPLACE)) return vectail_inplace(h, k);
80 90 : if (flag & LLL_KER) { setlg(h,k+1); return h; }
81 62 : g = vecslice(h,1,k); /* done first: vectail_inplace kills h */
82 70 : return mkvec2(g, vectail_inplace(h, k));
83 : }
84 :
85 : /* y * z * 2^e, e >= 0; y,z t_INT */
86 : INLINE GEN
87 1946899 : mulshift(GEN y, GEN z, long e)
88 : {
89 1946899 : long ly = lgefint(y), lz;
90 : pari_sp av;
91 : GEN t;
92 1946899 : if (ly == 2) return gen_0;
93 206968 : lz = lgefint(z);
94 206968 : av = avma; (void)new_chunk(ly+lz+nbits2lg(e)); /* HACK */
95 206968 : t = mulii(z, y);
96 206968 : set_avma(av); return shifti(t, e);
97 : }
98 :
99 : /* x - y * z * 2^e, e >= 0; x,y,z t_INT */
100 : INLINE GEN
101 6914675 : submulshift(GEN x, GEN y, GEN z, long e)
102 : {
103 6914675 : long lx = lgefint(x), ly, lz;
104 : pari_sp av;
105 : GEN t;
106 6914675 : if (!e) return submulii(x, y, z);
107 6879000 : if (lx == 2) { t = mulshift(y, z, e); togglesign(t); return t; }
108 5107627 : ly = lgefint(y);
109 5107627 : if (ly == 2) return icopy(x);
110 4629875 : lz = lgefint(z);
111 4629875 : av = avma; (void)new_chunk(lx+ly+lz+nbits2lg(e)); /* HACK */
112 4629875 : t = shifti(mulii(z, y), e);
113 4629875 : set_avma(av); return subii(x, t);
114 : }
115 :
116 : /* x - u*y * 2^e */
117 : INLINE GEN
118 23535325 : submuliu2n(GEN x, GEN y, ulong u, long e)
119 : {
120 : pari_sp av;
121 23535325 : long ly = lgefint(y);
122 23535325 : if (ly == 2) return x;
123 16550041 : av = avma;
124 16550041 : (void)new_chunk(3+ly+lgefint(x)+nbits2lg(e)); /* HACK */
125 16551733 : y = shifti(mului(u,y), e);
126 16550806 : set_avma(av); return subii(x, y);
127 : }
128 : /* x + u*y * 2^e */
129 : INLINE GEN
130 23802883 : addmuliu2n(GEN x, GEN y, ulong u, long e)
131 : {
132 : pari_sp av;
133 23802883 : long ly = lgefint(y);
134 23802883 : if (ly == 2) return x;
135 16666249 : av = avma;
136 16666249 : (void)new_chunk(3+ly+lgefint(x)+nbits2lg(e)); /* HACK */
137 16667836 : y = shifti(mului(u,y), e);
138 16667119 : set_avma(av); return addii(x, y);
139 : }
140 :
141 : /* n < 10; gerepileall supporting &NULL arguments. Maybe rename and export ? */
142 : INLINE void
143 7526 : gc_lll(pari_sp av, int n, ...)
144 : {
145 : int i, j;
146 : GEN *gptr[10];
147 : size_t s;
148 7526 : va_list a; va_start(a, n);
149 24200 : for (i=j=0; i<n; i++)
150 : {
151 16674 : GEN *x = va_arg(a,GEN*);
152 16674 : if (*x) { gptr[j++] = x; *x = (GEN)copy_bin(*x); }
153 : }
154 7526 : va_end(a); set_avma(av);
155 18894 : for (--j; j>=0; j--) *gptr[j] = bin_copy((GENbin*)*gptr[j]);
156 7526 : s = pari_mainstack->top - pari_mainstack->bot;
157 : /* size of saved objects ~ stacksize / 4 => overflow */
158 7526 : if (av - avma > (s >> 2))
159 : {
160 0 : size_t t = avma - pari_mainstack->bot;
161 0 : av = avma; new_chunk((s + t) / sizeof(long)); set_avma(av); /* double */
162 : }
163 7526 : }
164 :
165 : /********************************************************************/
166 : /** **/
167 : /** FPLLL (adapted from D. Stehle's code) **/
168 : /** **/
169 : /********************************************************************/
170 : /* Babai* and fplll* are a conversion to libpari API and data types
171 : of fplll-1.3 by Damien Stehle'.
172 :
173 : Copyright 2005, 2006 Damien Stehle'.
174 :
175 : This program is free software; you can redistribute it and/or modify it
176 : under the terms of the GNU General Public License as published by the
177 : Free Software Foundation; either version 2 of the License, or (at your
178 : option) any later version.
179 :
180 : This program implements ideas from the paper "Floating-point LLL Revisited",
181 : by Phong Nguyen and Damien Stehle', in the Proceedings of Eurocrypt'2005,
182 : Springer-Verlag; and was partly inspired by Shoup's NTL library:
183 : http://www.shoup.net/ntl/ */
184 :
185 : /* x t_REAL, |x| >= 1/2. Test whether |x| <= 3/2 */
186 : static int
187 1095308 : absrsmall2(GEN x)
188 : {
189 1095308 : long e = expo(x), l, i;
190 1095308 : if (e < 0) return 1;
191 344999 : if (e > 0 || (ulong)x[2] > (3UL << (BITS_IN_LONG-2))) return 0;
192 : /* line above assumes l > 2. OK since x != 0 */
193 178320 : l = lg(x); for (i = 3; i < l; i++) if (x[i]) return 0;
194 151059 : return 1;
195 : }
196 : /* x t_REAL; test whether |x| <= 1/2 */
197 : static int
198 2403848 : absrsmall(GEN x)
199 : {
200 : long e, l, i;
201 2403848 : if (!signe(x)) return 1;
202 2393033 : e = expo(x); if (e < -1) return 1;
203 1101635 : if (e > -1 || (ulong)x[2] > HIGHBIT) return 0;
204 7571 : l = lg(x); for (i = 3; i < l; i++) if (x[i]) return 0;
205 6327 : return 1;
206 : }
207 :
208 : static void
209 43794359 : rotate(GEN A, long k2, long k)
210 : {
211 : long i;
212 43794359 : GEN B = gel(A,k2);
213 140759355 : for (i = k2; i > k; i--) gel(A,i) = gel(A,i-1);
214 43794359 : gel(A,k) = B;
215 43794359 : }
216 :
217 : /************************* FAST version (double) ************************/
218 : #define dmael(x,i,j) ((x)[i][j])
219 : #define del(x,i) ((x)[i])
220 :
221 : static double *
222 37037420 : cget_dblvec(long d)
223 37037420 : { return (double*) stack_malloc_align(d*sizeof(double), sizeof(double)); }
224 :
225 : static double **
226 9307924 : cget_dblmat(long d) { return (double **) cgetg(d, t_VECSMALL); }
227 :
228 : static double
229 274753563 : itodbl_exp(GEN x, long *e)
230 : {
231 274753563 : pari_sp av = avma;
232 274753563 : GEN r = itor(x,DEFAULTPREC);
233 274685636 : *e = expo(r); setexpo(r,0);
234 274687389 : return gc_double(av, rtodbl(r));
235 : }
236 :
237 : static double
238 186155931 : dbldotproduct(double *x, double *y, long n)
239 : {
240 : long i;
241 186155931 : double sum = del(x,1) * del(y,1);
242 2969147867 : for (i=2; i<=n; i++) sum += del(x,i) * del(y,i);
243 186155931 : return sum;
244 : }
245 :
246 : static double
247 2674820 : dbldotsquare(double *x, long n)
248 : {
249 : long i;
250 2674820 : double sum = del(x,1) * del(x,1);
251 8606492 : for (i=2; i<=n; i++) sum += del(x,i) * del(x,i);
252 2674820 : return sum;
253 : }
254 :
255 : static long
256 32611705 : set_line(double *appv, GEN v, long n)
257 : {
258 32611705 : long i, maxexp = 0;
259 32611705 : pari_sp av = avma;
260 32611705 : GEN e = cgetg(n+1, t_VECSMALL);
261 307338846 : for (i = 1; i <= n; i++)
262 : {
263 274752809 : del(appv,i) = itodbl_exp(gel(v,i), e+i);
264 274725985 : if (e[i] > maxexp) maxexp = e[i];
265 : }
266 307468948 : for (i = 1; i <= n; i++) del(appv,i) = ldexp(del(appv,i), e[i]-maxexp);
267 32586037 : set_avma(av); return maxexp;
268 : }
269 :
270 : static void
271 47891975 : dblrotate(double **A, long k2, long k)
272 : {
273 : long i;
274 47891975 : double *B = del(A,k2);
275 153958996 : for (i = k2; i > k; i--) del(A,i) = del(A,i-1);
276 47891975 : del(A,k) = B;
277 47891975 : }
278 : /* update G[kappa][i] from appB */
279 : static void
280 30127160 : setG_fast(double **appB, long n, double **G, long kappa, long a, long b)
281 : { long i;
282 152203063 : for (i = a; i <= b; i++)
283 122076433 : dmael(G,kappa,i) = dbldotproduct(del(appB,kappa), del(appB,i), n);
284 30126630 : }
285 : /* update G[i][kappa] from appB */
286 : static void
287 24199319 : setG2_fast(double **appB, long n, double **G, long kappa, long a, long b)
288 : { long i;
289 88284639 : for (i = a; i <= b; i++)
290 64085570 : dmael(G,i,kappa) = dbldotproduct(del(appB,kappa), del(appB,i), n);
291 24199069 : }
292 : const long EX0 = -2; /* uninitialized; any value less than expo(0.51) = -1 */
293 :
294 : #ifdef LONG_IS_64BIT
295 : typedef long s64;
296 : #define addmuliu64_inplace addmuliu_inplace
297 : #define submuliu64_inplace submuliu_inplace
298 : #define submuliu642n submuliu2n
299 : #define addmuliu642n addmuliu2n
300 : #else
301 : typedef long long s64;
302 : typedef unsigned long long u64;
303 :
304 : INLINE GEN
305 42949151 : u64toi(u64 x)
306 : {
307 : GEN y;
308 : ulong h;
309 42949151 : if (!x) return gen_0;
310 42949151 : h = x>>32;
311 42949151 : if (!h) return utoipos(x);
312 4038329 : y = cgetipos(4);
313 4038329 : *int_LSW(y) = x&0xFFFFFFFF;
314 4038329 : *int_MSW(y) = x>>32;
315 4038329 : return y;
316 : }
317 :
318 : INLINE GEN
319 3309490 : u64toineg(u64 x)
320 : {
321 : GEN y;
322 : ulong h;
323 3309490 : if (!x) return gen_0;
324 3309490 : h = x>>32;
325 3309490 : if (!h) return utoineg(x);
326 3309490 : y = cgetineg(4);
327 3309490 : *int_LSW(y) = x&0xFFFFFFFF;
328 3309490 : *int_MSW(y) = x>>32;
329 3309490 : return y;
330 : }
331 : INLINE GEN
332 19770843 : addmuliu64_inplace(GEN x, GEN y, u64 u) { return addmulii(x, y, u64toi(u)); }
333 :
334 : INLINE GEN
335 19926645 : submuliu64_inplace(GEN x, GEN y, u64 u) { return submulii(x, y, u64toi(u)); }
336 :
337 : INLINE GEN
338 3309490 : addmuliu642n(GEN x, GEN y, u64 u, long e) { return submulshift(x, y, u64toineg(u), e); }
339 :
340 : INLINE GEN
341 3251663 : submuliu642n(GEN x, GEN y, u64 u, long e) { return submulshift(x, y, u64toi(u), e); }
342 :
343 : #endif
344 :
345 : /* Babai's Nearest Plane algorithm (iterative); see Babai() */
346 : static int
347 40768813 : Babai_fast(pari_sp av, long kappa, GEN *pB, GEN *pU, double **mu, double **r,
348 : double *s, double **appB, GEN expoB, double **G,
349 : long a, long zeros, long maxG, double eta)
350 : {
351 40768813 : GEN B = *pB, U = *pU;
352 40768813 : const long n = nbrows(B), d = U ? lg(U)-1: 0;
353 40768516 : long k, aa = (a > zeros)? a : zeros+1;
354 40768516 : long emaxmu = EX0, emax2mu = EX0;
355 : s64 xx;
356 40768516 : int did_something = 0;
357 : /* N.B: we set d = 0 (resp. n = 0) to avoid updating U (resp. B) */
358 :
359 24518052 : for (;;) {
360 65286568 : int go_on = 0;
361 65286568 : long i, j, emax3mu = emax2mu;
362 :
363 65286568 : if (gc_needed(av,2))
364 : {
365 479 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[1], a=%ld", aa);
366 479 : gc_lll(av,2,&B,&U);
367 : }
368 : /* Step2: compute the GSO for stage kappa */
369 65285140 : emax2mu = emaxmu; emaxmu = EX0;
370 276327412 : for (j=aa; j<kappa; j++)
371 : {
372 211043208 : double g = dmael(G,kappa,j);
373 1119687488 : for (k = zeros+1; k < j; k++) g -= dmael(mu,j,k) * dmael(r,kappa,k);
374 211043208 : dmael(r,kappa,j) = g;
375 211043208 : dmael(mu,kappa,j) = dmael(r,kappa,j) / dmael(r,j,j);
376 211043208 : emaxmu = maxss(emaxmu, expoB[kappa]-expoB[j]);
377 : }
378 : /* maxmu doesn't decrease fast enough */
379 65284204 : if (emax3mu != EX0 && emax3mu <= emax2mu + 5) {*pB = B; *pU = U; return 1;}
380 :
381 266441535 : for (j=kappa-1; j>zeros; j--)
382 : {
383 225679101 : double tmp = fabs(ldexp (dmael(mu,kappa,j), expoB[kappa]-expoB[j]));
384 225679101 : if (tmp>eta) { go_on = 1; break; }
385 : }
386 :
387 : /* Step3--5: compute the X_j's */
388 65280757 : if (go_on)
389 126705921 : for (j=kappa-1; j>zeros; j--)
390 : { /* The code below seemingly handles U = NULL, but in this case d = 0 */
391 102207925 : int e = expoB[j] - expoB[kappa];
392 102207925 : double tmp = ldexp(dmael(mu,kappa,j), -e), atmp = fabs(tmp);
393 : /* tmp = Inf is allowed */
394 102207925 : if (atmp <= .5) continue; /* size-reduced */
395 56909007 : if (gc_needed(av,2))
396 : {
397 1880 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[2], a=%ld, j=%ld", aa,j);
398 1880 : gc_lll(av,2,&B,&U);
399 : }
400 56911505 : did_something = 1;
401 : /* we consider separately the case |X| = 1 */
402 56911505 : if (atmp <= 1.5)
403 : {
404 37176639 : if (dmael(mu,kappa,j) > 0) { /* in this case, X = 1 */
405 95564236 : for (k=zeros+1; k<j; k++)
406 76445474 : dmael(mu,kappa,k) -= ldexp(dmael(mu,j,k), e);
407 347538831 : for (i=1; i<=n; i++)
408 328431338 : gmael(B,kappa,i) = subii(gmael(B,kappa,i), gmael(B,j,i));
409 120374972 : for (i=1; i<=d; i++)
410 101267246 : gmael(U,kappa,i) = subii(gmael(U,kappa,i), gmael(U,j,i));
411 : } else { /* otherwise X = -1 */
412 93967246 : for (k=zeros+1; k<j; k++)
413 75909369 : dmael(mu,kappa,k) += ldexp(dmael(mu,j,k), e);
414 341533678 : for (i=1; i<=n; i++)
415 323486036 : gmael(B,kappa,i) = addii(gmael(B,kappa,i), gmael(B,j,i));
416 114799660 : for (i=1; i<=d; i++)
417 96751880 : gmael(U,kappa,i) = addii(gmael(U,kappa,i), gmael(U,j,i));
418 : }
419 37155506 : continue;
420 : }
421 : /* we have |X| >= 2 */
422 19734866 : if (atmp < 9007199254740992.)
423 : {
424 17343643 : tmp = pari_rint(tmp);
425 44504243 : for (k=zeros+1; k<j; k++)
426 27160601 : dmael(mu,kappa,k) -= ldexp(tmp * dmael(mu,j,k), e);
427 17343642 : xx = (s64) tmp;
428 17343642 : if (xx > 0) /* = xx */
429 : {
430 94706934 : for (i=1; i<=n; i++)
431 85989449 : gmael(B,kappa,i) = submuliu64_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
432 62242490 : for (i=1; i<=d; i++)
433 53524873 : gmael(U,kappa,i) = submuliu64_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
434 : }
435 : else /* = -xx */
436 : {
437 93911514 : for (i=1; i<=n; i++)
438 85286889 : gmael(B,kappa,i) = addmuliu64_inplace(gmael(B,kappa,i), gmael(B,j,i), -xx);
439 61497995 : for (i=1; i<=d; i++)
440 52873217 : gmael(U,kappa,i) = addmuliu64_inplace(gmael(U,kappa,i), gmael(U,j,i), -xx);
441 : }
442 : }
443 : else
444 : {
445 : int E;
446 2391223 : xx = (s64) ldexp(frexp(dmael(mu,kappa,j), &E), 53);
447 2391223 : E -= e + 53;
448 2391223 : if (E <= 0)
449 : {
450 0 : xx = xx << -E;
451 0 : for (k=zeros+1; k<j; k++)
452 0 : dmael(mu,kappa,k) -= ldexp(((double)xx) * dmael(mu,j,k), e);
453 0 : if (xx > 0) /* = xx */
454 : {
455 0 : for (i=1; i<=n; i++)
456 0 : gmael(B,kappa,i) = submuliu64_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
457 0 : for (i=1; i<=d; i++)
458 0 : gmael(U,kappa,i) = submuliu64_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
459 : }
460 : else /* = -xx */
461 : {
462 0 : for (i=1; i<=n; i++)
463 0 : gmael(B,kappa,i) = addmuliu64_inplace(gmael(B,kappa,i), gmael(B,j,i), -xx);
464 0 : for (i=1; i<=d; i++)
465 0 : gmael(U,kappa,i) = addmuliu64_inplace(gmael(U,kappa,i), gmael(U,j,i), -xx);
466 : }
467 : } else
468 : {
469 15961730 : for (k=zeros+1; k<j; k++)
470 13570507 : dmael(mu,kappa,k) -= ldexp(((double)xx) * dmael(mu,j,k), E + e);
471 2391223 : if (xx > 0) /* = xx */
472 : {
473 22707755 : for (i=1; i<=n; i++)
474 21518924 : gmael(B,kappa,i) = submuliu642n(gmael(B,kappa,i), gmael(B,j,i), xx, E);
475 2265022 : for (i=1; i<=d; i++)
476 1076191 : gmael(U,kappa,i) = submuliu642n(gmael(U,kappa,i), gmael(U,j,i), xx, E);
477 : }
478 : else /* = -xx */
479 : {
480 23320863 : for (i=1; i<=n; i++)
481 22118830 : gmael(B,kappa,i) = addmuliu642n(gmael(B,kappa,i), gmael(B,j,i), -xx, E);
482 2258739 : for (i=1; i<=d; i++)
483 1056706 : gmael(U,kappa,i) = addmuliu642n(gmael(U,kappa,i), gmael(U,j,i), -xx, E);
484 : }
485 : }
486 : }
487 : }
488 65260515 : if (!go_on) break; /* Anything happened? */
489 24495513 : expoB[kappa] = set_line(del(appB,kappa), gel(B,kappa), n);
490 24517433 : setG_fast(appB, n, G, kappa, zeros+1, kappa-1);
491 24518052 : aa = zeros+1;
492 : }
493 40765002 : if (did_something) setG2_fast(appB, n, G, kappa, kappa, maxG);
494 :
495 40765471 : del(s,zeros+1) = dmael(G,kappa,kappa);
496 : /* the last s[kappa-1]=r[kappa][kappa] is computed only if kappa increases */
497 183334809 : for (k=zeros+1; k<=kappa-2; k++)
498 142569338 : del(s,k+1) = del(s,k) - dmael(mu,kappa,k)*dmael(r,kappa,k);
499 40765471 : *pB = B; *pU = U; return 0;
500 : }
501 :
502 : static void
503 17701252 : update_alpha(GEN alpha, long kappa, long kappa2, long kappamax)
504 : {
505 : long i;
506 58481975 : for (i = kappa; i < kappa2; i++)
507 40780723 : if (kappa <= alpha[i]) alpha[i] = kappa;
508 58482028 : for (i = kappa2; i > kappa; i--) alpha[i] = alpha[i-1];
509 48023003 : for (i = kappa2+1; i <= kappamax; i++)
510 30321751 : if (kappa < alpha[i]) alpha[i] = kappa;
511 17701252 : alpha[kappa] = kappa;
512 17701252 : }
513 : static void
514 1737043 : rotateG(GEN G, long kappa2, long kappa, long maxG, GEN Gtmp)
515 : {
516 : long i, j;
517 20963591 : for (i=1; i<=kappa2; i++) gel(Gtmp,i) = gmael(G,kappa2,i);
518 10456122 : for ( ; i<=maxG; i++) gel(Gtmp,i) = gmael(G,i,kappa2);
519 7161434 : for (i=kappa2; i>kappa; i--)
520 : {
521 38112995 : for (j=1; j<kappa; j++) gmael(G,i,j) = gmael(G,i-1,j);
522 5424391 : gmael(G,i,kappa) = gel(Gtmp,i-1);
523 26054763 : for (j=kappa+1; j<=i; j++) gmael(G,i,j) = gmael(G,i-1,j-1);
524 26994555 : for (j=kappa2+1; j<=maxG; j++) gmael(G,j,i) = gmael(G,j,i-1);
525 : }
526 13802157 : for (i=1; i<kappa; i++) gmael(G,kappa,i) = gel(Gtmp,i);
527 1737043 : gmael(G,kappa,kappa) = gel(Gtmp,kappa2);
528 10456122 : for (i=kappa2+1; i<=maxG; i++) gmael(G,i,kappa) = gel(Gtmp,i);
529 1737043 : }
530 : static void
531 15963995 : rotateG_fast(double **G, long kappa2, long kappa, long maxG, double *Gtmp)
532 : {
533 : long i, j;
534 101898432 : for (i=1; i<=kappa2; i++) del(Gtmp,i) = dmael(G,kappa2,i);
535 38551532 : for ( ; i<=maxG; i++) del(Gtmp,i) = dmael(G,i,kappa2);
536 51319961 : for (i=kappa2; i>kappa; i--)
537 : {
538 121033075 : for (j=1; j<kappa; j++) dmael(G,i,j) = dmael(G,i-1,j);
539 35355966 : dmael(G,i,kappa) = del(Gtmp,i-1);
540 128703744 : for (j=kappa+1; j<=i; j++) dmael(G,i,j) = dmael(G,i-1,j-1);
541 86313214 : for (j=kappa2+1; j<=maxG; j++) dmael(G,j,i) = dmael(G,j,i-1);
542 : }
543 50579098 : for (i=1; i<kappa; i++) dmael(G,kappa,i) = del(Gtmp,i);
544 15963995 : dmael(G,kappa,kappa) = del(Gtmp,kappa2);
545 38551602 : for (i=kappa2+1; i<=maxG; i++) dmael(G,i,kappa) = del(Gtmp,i);
546 15963995 : }
547 :
548 : /* LLL-reduces (B,U) in place [apply base change transforms to B and U].
549 : * Gram matrix, and GSO performed on matrices of 'double'.
550 : * If (keepfirst), never swap with first vector.
551 : * Return -1 on failure, else zeros = dim Kernel (>= 0) */
552 : static long
553 2327010 : fplll_fast(GEN *pB, GEN *pU, double delta, double eta, long keepfirst)
554 : {
555 : pari_sp av;
556 : long kappa, kappa2, d, n, i, j, zeros, kappamax, maxG;
557 : double **mu, **r, *s, tmp, *Gtmp, **G, **appB;
558 2327010 : GEN alpha, expoB, B = *pB, U;
559 2327010 : long cnt = 0;
560 :
561 2327010 : d = lg(B)-1;
562 2327010 : n = nbrows(B);
563 2327010 : U = *pU; /* NULL if inplace */
564 :
565 2327010 : G = cget_dblmat(d+1);
566 2327000 : appB = cget_dblmat(d+1);
567 2327000 : mu = cget_dblmat(d+1);
568 2326995 : r = cget_dblmat(d+1);
569 2326999 : s = cget_dblvec(d+1);
570 10423047 : for (j = 1; j <= d; j++)
571 : {
572 8096043 : del(mu,j) = cget_dblvec(d+1);
573 8096045 : del(r,j) = cget_dblvec(d+1);
574 8096068 : del(appB,j) = cget_dblvec(n+1);
575 8096048 : del(G,j) = cget_dblvec(d+1);
576 : }
577 2327004 : expoB = cgetg(d+1, t_VECSMALL);
578 10422920 : for (i=1; i<=d; i++) expoB[i] = set_line(del(appB,i), gel(B,i), n);
579 2326938 : Gtmp = cget_dblvec(d+1);
580 2326993 : alpha = cgetg(d+1, t_VECSMALL);
581 2326983 : av = avma;
582 :
583 : /* Step2: Initializing the main loop */
584 2326983 : kappamax = 1;
585 2326983 : i = 1;
586 2326983 : maxG = d; /* later updated to kappamax */
587 :
588 : do {
589 2480046 : dmael(G,i,i) = dbldotsquare(del(appB,i),n);
590 2480053 : } while (dmael(G,i,i) <= 0 && (++i <=d));
591 2326990 : zeros = i-1; /* all vectors B[i] with i <= zeros are zero vectors */
592 2326990 : kappa = i;
593 2326990 : if (zeros < d) dmael(r,zeros+1,zeros+1) = dmael(G,zeros+1,zeros+1);
594 10269983 : for (i=zeros+1; i<=d; i++) alpha[i]=1;
595 43092596 : while (++kappa <= d)
596 : {
597 40769016 : if (kappa > kappamax)
598 : {
599 5609776 : if (DEBUGLEVEL>=4) err_printf("K%ld ",kappa);
600 5609776 : maxG = kappamax = kappa;
601 5609776 : setG_fast(appB, n, G, kappa, zeros+1, kappa);
602 : }
603 : /* Step3: Call to the Babai algorithm, mu,r,s updated in place */
604 40769009 : if (Babai_fast(av, kappa, &B,&U, mu,r,s, appB, expoB, G, alpha[kappa],
605 3447 : zeros, maxG, eta)) { *pB=B; *pU=U; return -1; }
606 :
607 40765362 : tmp = ldexp(r[kappa-1][kappa-1] * delta, 2*(expoB[kappa-1]-expoB[kappa]));
608 40765362 : if ((keepfirst && kappa == 2) || tmp <= del(s,kappa-1))
609 : { /* Step4: Success of Lovasz's condition */
610 24801166 : alpha[kappa] = kappa;
611 24801166 : tmp = dmael(mu,kappa,kappa-1) * dmael(r,kappa,kappa-1);
612 24801166 : dmael(r,kappa,kappa) = del(s,kappa-1)- tmp;
613 24801166 : continue;
614 : }
615 : /* Step5: Find the right insertion index kappa, kappa2 = initial kappa */
616 15964196 : if (DEBUGLEVEL>=4 && kappa==kappamax && del(s,kappa-1)!=0)
617 0 : if (++cnt > 20) { cnt = 0; err_printf("(%ld) ", 2*expoB[1] + dblexpo(del(s,1))); }
618 15964191 : kappa2 = kappa;
619 : do {
620 35356262 : kappa--;
621 35356262 : if (kappa<zeros+2 + (keepfirst ? 1: 0)) break;
622 26541908 : tmp = dmael(r,kappa-1,kappa-1) * delta;
623 26541908 : tmp = ldexp(tmp, 2*(expoB[kappa-1]-expoB[kappa2]));
624 26541908 : } while (del(s,kappa-1) <= tmp);
625 15964191 : update_alpha(alpha, kappa, kappa2, kappamax);
626 :
627 : /* Step6: Update the mu's and r's */
628 15964241 : dblrotate(mu,kappa2,kappa);
629 15964199 : dblrotate(r,kappa2,kappa);
630 15964152 : dmael(r,kappa,kappa) = del(s,kappa);
631 :
632 : /* Step7: Update B, appB, U, G */
633 15964152 : rotate(B,kappa2,kappa);
634 15964096 : dblrotate(appB,kappa2,kappa);
635 15964057 : if (U) rotate(U,kappa2,kappa);
636 15964050 : rotate(expoB,kappa2,kappa);
637 15964005 : rotateG_fast(G,kappa2,kappa, maxG, Gtmp);
638 :
639 : /* Step8: Prepare the next loop iteration */
640 15964440 : if (kappa == zeros+1 && dmael(G,kappa,kappa)<= 0)
641 : {
642 194770 : zeros++; kappa++;
643 194770 : dmael(G,kappa,kappa) = dbldotsquare(del(appB,kappa),n);
644 194770 : dmael(r,kappa,kappa) = dmael(G,kappa,kappa);
645 : }
646 : }
647 2323580 : *pB = B; *pU = U; return zeros;
648 : }
649 :
650 : /***************** HEURISTIC version (reduced precision) ****************/
651 : static GEN
652 340849 : realsqrdotproduct(GEN x)
653 : {
654 340849 : long i, l = lg(x);
655 340849 : GEN z = sqrr(gel(x,1));
656 4236909 : for (i=2; i<l; i++) z = addrr(z, sqrr(gel(x,i)));
657 340849 : return z;
658 : }
659 : /* x, y non-empty vector of t_REALs, same length */
660 : static GEN
661 3855818 : realdotproduct(GEN x, GEN y)
662 : {
663 : long i, l;
664 : GEN z;
665 3855818 : if (x == y) return realsqrdotproduct(x);
666 3514969 : l = lg(x); z = mulrr(gel(x,1),gel(y,1));
667 57087233 : for (i=2; i<l; i++) z = addrr(z, mulrr(gel(x,i), gel(y,i)));
668 3514969 : return z;
669 : }
670 : static void
671 348202 : setG_heuristic(GEN appB, GEN G, long kappa, long a, long b)
672 348202 : { pari_sp av = avma;
673 : long i;
674 2785548 : for (i = a; i <= b; i++)
675 2437346 : affrr(realdotproduct(gel(appB,kappa),gel(appB,i)), gmael(G,kappa,i));
676 348202 : set_avma(av);
677 348202 : }
678 : static void
679 329719 : setG2_heuristic(GEN appB, GEN G, long kappa, long a, long b)
680 329719 : { pari_sp av = avma;
681 : long i;
682 1748191 : for (i = a; i <= b; i++)
683 1418472 : affrr(realdotproduct(gel(appB,kappa),gel(appB,i)), gmael(G,i,kappa));
684 329719 : set_avma(av);
685 329719 : }
686 :
687 : /* approximate t_REAL x as m * 2^e, where |m| < 2^bit */
688 : static GEN
689 21060 : truncexpo(GEN x, long bit, long *e)
690 : {
691 21060 : *e = expo(x) + 1 - bit;
692 21060 : if (*e >= 0) return mantissa2nr(x, 0);
693 1870 : *e = 0; return roundr_safe(x);
694 : }
695 : /* Babai's Nearest Plane algorithm (iterative); see Babai() */
696 : static int
697 677423 : Babai_heuristic(pari_sp av, long kappa, GEN *pB, GEN *pU, GEN mu, GEN r, GEN s,
698 : GEN appB, GEN G, long a, long zeros, long maxG,
699 : GEN eta, long prec)
700 : {
701 677423 : GEN B = *pB, U = *pU;
702 677423 : const long n = nbrows(B), d = U ? lg(U)-1: 0, bit = prec2nbits(prec);
703 677423 : long k, aa = (a > zeros)? a : zeros+1;
704 677423 : int did_something = 0;
705 677423 : long emaxmu = EX0, emax2mu = EX0;
706 : /* N.B: we set d = 0 (resp. n = 0) to avoid updating U (resp. B) */
707 :
708 337072 : for (;;) {
709 1014495 : int go_on = 0;
710 1014495 : long i, j, emax3mu = emax2mu;
711 :
712 1014495 : if (gc_needed(av,2))
713 : {
714 14 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[1], a=%ld", aa);
715 14 : gc_lll(av,2,&B,&U);
716 : }
717 : /* Step2: compute the GSO for stage kappa */
718 1014495 : emax2mu = emaxmu; emaxmu = EX0;
719 5723022 : for (j=aa; j<kappa; j++)
720 : {
721 4708527 : pari_sp btop = avma;
722 4708527 : GEN g = gmael(G,kappa,j);
723 30245965 : for (k = zeros+1; k<j; k++)
724 25537438 : g = subrr(g, mulrr(gmael(mu,j,k), gmael(r,kappa,k)));
725 4708527 : affrr(g, gmael(r,kappa,j));
726 4708527 : affrr(divrr(gmael(r,kappa,j), gmael(r,j,j)), gmael(mu,kappa,j));
727 4708527 : emaxmu = maxss(emaxmu, expo(gmael(mu,kappa,j)));
728 4708527 : set_avma(btop);
729 : }
730 1014495 : if (emax3mu != EX0 && emax3mu <= emax2mu + 5)
731 1146 : { *pB = B; *pU = U; return 1; }
732 :
733 6250095 : for (j=kappa-1; j>zeros; j--)
734 5573818 : if (abscmprr(gmael(mu,kappa,j), eta) > 0) { go_on = 1; break; }
735 :
736 : /* Step3--5: compute the X_j's */
737 1013349 : if (go_on)
738 2726189 : for (j=kappa-1; j>zeros; j--)
739 : { /* The code below seemingly handles U = NULL, but in this case d = 0 */
740 : pari_sp btop;
741 2389117 : GEN tmp = gmael(mu,kappa,j);
742 2389117 : if (absrsmall(tmp)) continue; /* size-reduced */
743 :
744 1086567 : if (gc_needed(av,2))
745 : {
746 84 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[2], a=%ld, j=%ld", aa,j);
747 84 : gc_lll(av,2,&B,&U);
748 : }
749 1086567 : btop = avma; did_something = 1;
750 : /* we consider separately the case |X| = 1 */
751 1086567 : if (absrsmall2(tmp))
752 : {
753 897647 : if (signe(tmp) > 0) { /* in this case, X = 1 */
754 2582780 : for (k=zeros+1; k<j; k++)
755 2134607 : affrr(subrr(gmael(mu,kappa,k), gmael(mu,j,k)), gmael(mu,kappa,k));
756 448173 : set_avma(btop);
757 7894422 : for (i=1; i<=n; i++)
758 7446249 : gmael(B,kappa,i) = subii(gmael(B,kappa,i), gmael(B,j,i));
759 6676484 : for (i=1; i<=d; i++)
760 6228311 : gmael(U,kappa,i) = subii(gmael(U,kappa,i), gmael(U,j,i));
761 : } else { /* otherwise X = -1 */
762 2586266 : for (k=zeros+1; k<j; k++)
763 2136792 : affrr(addrr(gmael(mu,kappa,k), gmael(mu,j,k)), gmael(mu,kappa,k));
764 449474 : set_avma(btop);
765 7921890 : for (i=1; i<=n; i++)
766 7472416 : gmael(B,kappa,i) = addii(gmael(B,kappa,i), gmael(B,j,i));
767 6700869 : for (i=1; i<=d; i++)
768 6251395 : gmael(U,kappa,i) = addii(gmael(U,kappa,i),gmael(U,j,i));
769 : }
770 897647 : continue;
771 : }
772 : /* we have |X| >= 2 */
773 188920 : if (expo(tmp) < BITS_IN_LONG)
774 : {
775 168380 : ulong xx = roundr_safe(tmp)[2]; /* X fits in an ulong */
776 168380 : if (signe(tmp) > 0) /* = xx */
777 : {
778 307368 : for (k=zeros+1; k<j; k++)
779 223569 : affrr(subrr(gmael(mu,kappa,k), mulur(xx, gmael(mu,j,k))),
780 223569 : gmael(mu,kappa,k));
781 83799 : set_avma(btop);
782 1110001 : for (i=1; i<=n; i++)
783 1026202 : gmael(B,kappa,i) = submuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
784 919889 : for (i=1; i<=d; i++)
785 836090 : gmael(U,kappa,i) = submuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
786 : }
787 : else /* = -xx */
788 : {
789 309539 : for (k=zeros+1; k<j; k++)
790 224958 : affrr(addrr(gmael(mu,kappa,k), mulur(xx, gmael(mu,j,k))),
791 224958 : gmael(mu,kappa,k));
792 84581 : set_avma(btop);
793 1113565 : for (i=1; i<=n; i++)
794 1028984 : gmael(B,kappa,i) = addmuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
795 915355 : for (i=1; i<=d; i++)
796 830774 : gmael(U,kappa,i) = addmuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
797 : }
798 : }
799 : else
800 : {
801 : long e;
802 20540 : GEN X = truncexpo(tmp, bit, &e); /* tmp ~ X * 2^e */
803 20540 : btop = avma;
804 88630 : for (k=zeros+1; k<j; k++)
805 : {
806 68090 : GEN x = mulir(X, gmael(mu,j,k));
807 68090 : if (e) shiftr_inplace(x, e);
808 68090 : affrr(subrr(gmael(mu,kappa,k), x), gmael(mu,kappa,k));
809 : }
810 20540 : set_avma(btop);
811 265641 : for (i=1; i<=n; i++)
812 245101 : gmael(B,kappa,i) = submulshift(gmael(B,kappa,i), gmael(B,j,i), X, e);
813 98271 : for (i=1; i<=d; i++)
814 77731 : gmael(U,kappa,i) = submulshift(gmael(U,kappa,i), gmael(U,j,i), X, e);
815 : }
816 : }
817 1013349 : if (!go_on) break; /* Anything happened? */
818 4546534 : for (i=1 ; i<=n; i++) affir(gmael(B,kappa,i), gmael(appB,kappa,i));
819 337072 : setG_heuristic(appB, G, kappa, zeros+1, kappa-1);
820 337072 : aa = zeros+1;
821 : }
822 676277 : if (did_something) setG2_heuristic(appB, G, kappa, kappa, maxG);
823 676277 : affrr(gmael(G,kappa,kappa), gel(s,zeros+1));
824 : /* the last s[kappa-1]=r[kappa][kappa] is computed only if kappa increases */
825 676277 : av = avma;
826 4790819 : for (k=zeros+1; k<=kappa-2; k++)
827 4114542 : affrr(subrr(gel(s,k), mulrr(gmael(mu,kappa,k), gmael(r,kappa,k))),
828 4114542 : gel(s,k+1));
829 676277 : *pB = B; *pU = U; return gc_bool(av, 0);
830 : }
831 :
832 : static GEN
833 17816 : ZC_to_RC(GEN x, long prec)
834 177243 : { pari_APPLY_type(t_COL,itor(gel(x,i),prec)) }
835 :
836 : static GEN
837 3447 : ZM_to_RM(GEN x, long prec)
838 21263 : { pari_APPLY_same(ZC_to_RC(gel(x,i),prec)) }
839 :
840 : /* LLL-reduces (B,U) in place [apply base change transforms to B and U].
841 : * Gram matrix made of t_REAL at precision prec2, performe GSO at prec.
842 : * If (keepfirst), never swap with first vector.
843 : * Return -1 on failure, else zeros = dim Kernel (>= 0) */
844 : static long
845 3447 : fplll_heuristic(GEN *pB, GEN *pU, double DELTA, double ETA, long keepfirst,
846 : long prec, long prec2)
847 : {
848 : pari_sp av, av2;
849 : long kappa, kappa2, d, i, j, zeros, kappamax, maxG;
850 3447 : GEN mu, r, s, tmp, Gtmp, alpha, G, appB, B = *pB, U;
851 3447 : GEN delta = dbltor(DELTA), eta = dbltor(ETA);
852 3447 : long cnt = 0;
853 :
854 3447 : d = lg(B)-1;
855 3447 : U = *pU; /* NULL if inplace */
856 :
857 3447 : G = cgetg(d+1, t_MAT);
858 3447 : mu = cgetg(d+1, t_MAT);
859 3447 : r = cgetg(d+1, t_MAT);
860 3447 : s = cgetg(d+1, t_VEC);
861 3447 : appB = ZM_to_RM(B, prec2);
862 21263 : for (j = 1; j <= d; j++)
863 : {
864 17816 : GEN M = cgetg(d+1, t_COL), R = cgetg(d+1, t_COL), S = cgetg(d+1, t_COL);
865 17816 : gel(mu,j)= M;
866 17816 : gel(r,j) = R;
867 17816 : gel(G,j) = S;
868 17816 : gel(s,j) = cgetr(prec);
869 168482 : for (i = 1; i <= d; i++)
870 : {
871 150666 : gel(R,i) = cgetr(prec);
872 150666 : gel(M,i) = cgetr(prec);
873 150666 : gel(S,i) = cgetr(prec2);
874 : }
875 : }
876 3447 : Gtmp = cgetg(d+1, t_VEC);
877 3447 : alpha = cgetg(d+1, t_VECSMALL);
878 3447 : av = avma;
879 :
880 : /* Step2: Initializing the main loop */
881 3447 : kappamax = 1;
882 3447 : i = 1;
883 3447 : maxG = d; /* later updated to kappamax */
884 :
885 : do {
886 3447 : affrr(RgV_dotsquare(gel(appB,i)), gmael(G,i,i));
887 3447 : } while (signe(gmael(G,i,i)) == 0 && (++i <=d));
888 3447 : zeros = i-1; /* all vectors B[i] with i <= zeros are zero vectors */
889 3447 : kappa = i;
890 3447 : if (zeros < d) affrr(gmael(G,zeros+1,zeros+1), gmael(r,zeros+1,zeros+1));
891 21263 : for (i=zeros+1; i<=d; i++) alpha[i]=1;
892 :
893 679724 : while (++kappa <= d)
894 : {
895 677423 : if (kappa > kappamax)
896 : {
897 11130 : if (DEBUGLEVEL>=4) err_printf("K%ld ",kappa);
898 11130 : maxG = kappamax = kappa;
899 11130 : setG_heuristic(appB, G, kappa, zeros+1, kappa);
900 : }
901 : /* Step3: Call to the Babai algorithm, mu,r,s updated in place */
902 677423 : if (Babai_heuristic(av, kappa, &B,&U, mu,r,s, appB, G, alpha[kappa], zeros,
903 1146 : maxG, eta, prec)) { *pB = B; *pU = U; return -1; }
904 676277 : av2 = avma;
905 1352462 : if ((keepfirst && kappa == 2) ||
906 676185 : cmprr(mulrr(gmael(r,kappa-1,kappa-1), delta), gel(s,kappa-1)) <= 0)
907 : { /* Step4: Success of Lovasz's condition */
908 456943 : alpha[kappa] = kappa;
909 456943 : tmp = mulrr(gmael(mu,kappa,kappa-1), gmael(r,kappa,kappa-1));
910 456943 : affrr(subrr(gel(s,kappa-1), tmp), gmael(r,kappa,kappa));
911 456943 : set_avma(av2); continue;
912 : }
913 : /* Step5: Find the right insertion index kappa, kappa2 = initial kappa */
914 219334 : if (DEBUGLEVEL>=4 && kappa==kappamax && signe(gel(s,kappa-1)))
915 0 : if (++cnt > 20) { cnt = 0; err_printf("(%ld) ", expo(gel(s,1))); }
916 219334 : kappa2 = kappa;
917 : do {
918 666320 : kappa--;
919 666320 : if (kappa < zeros+2 + (keepfirst ? 1: 0)) break;
920 636368 : tmp = mulrr(gmael(r,kappa-1,kappa-1), delta);
921 636368 : } while (cmprr(gel(s,kappa-1), tmp) <= 0 );
922 219334 : set_avma(av2);
923 219334 : update_alpha(alpha, kappa, kappa2, kappamax);
924 :
925 : /* Step6: Update the mu's and r's */
926 219334 : rotate(mu,kappa2,kappa);
927 219334 : rotate(r,kappa2,kappa);
928 219334 : affrr(gel(s,kappa), gmael(r,kappa,kappa));
929 :
930 : /* Step7: Update B, appB, U, G */
931 219334 : rotate(B,kappa2,kappa);
932 219334 : rotate(appB,kappa2,kappa);
933 219334 : if (U) rotate(U,kappa2,kappa);
934 219334 : rotateG(G,kappa2,kappa, maxG, Gtmp);
935 :
936 : /* Step8: Prepare the next loop iteration */
937 219334 : if (kappa == zeros+1 && !signe(gmael(G,kappa,kappa)))
938 : {
939 6 : zeros++; kappa++;
940 6 : affrr(RgV_dotsquare(gel(appB,kappa)), gmael(G,kappa,kappa));
941 6 : affrr(gmael(G,kappa,kappa), gmael(r,kappa,kappa));
942 : }
943 : }
944 2301 : *pB=B; *pU=U; return zeros;
945 : }
946 :
947 : /************************* PROVED version (t_INT) ***********************/
948 : /* dpe inspired by dpe.h by Patrick Pelissier, Paul Zimmermann
949 : * https://gforge.inria.fr/projects/dpe/
950 : */
951 :
952 : typedef struct
953 : {
954 : double d; /* significand */
955 : long e; /* exponent */
956 : } dpe_t;
957 :
958 : #define Dmael(x,i,j) (&((x)[i][j]))
959 : #define Del(x,i) (&((x)[i]))
960 :
961 : static void
962 3033930 : dperotate(dpe_t **A, long k2, long k)
963 : {
964 : long i;
965 3033930 : dpe_t *B = A[k2];
966 12545734 : for (i = k2; i > k; i--) A[i] = A[i-1];
967 3033930 : A[k] = B;
968 3033930 : }
969 :
970 : static void
971 859035435 : dpe_normalize0(dpe_t *x)
972 : {
973 : int e;
974 859035435 : x->d = frexp(x->d, &e);
975 859035435 : x->e += e;
976 859035435 : }
977 :
978 : static void
979 423805195 : dpe_normalize(dpe_t *x)
980 : {
981 423805195 : if (x->d == 0.0)
982 2301693 : x->e = -LONG_MAX;
983 : else
984 421503502 : dpe_normalize0(x);
985 423805335 : }
986 :
987 : static GEN
988 83121 : dpetor(dpe_t *x)
989 : {
990 83121 : GEN r = dbltor(x->d);
991 83121 : if (signe(r)==0) return r;
992 83121 : setexpo(r, x->e-1);
993 83121 : return r;
994 : }
995 :
996 : static void
997 96288774 : affdpe(dpe_t *y, dpe_t *x)
998 : {
999 96288774 : x->d = y->d;
1000 96288774 : x->e = y->e;
1001 96288774 : }
1002 :
1003 : static void
1004 67843691 : affidpe(GEN y, dpe_t *x)
1005 : {
1006 67843691 : pari_sp av = avma;
1007 67843691 : GEN r = itor(y, DEFAULTPREC);
1008 67842170 : x->e = expo(r)+1;
1009 67842170 : setexpo(r,-1);
1010 67842052 : x->d = rtodbl(r);
1011 67842831 : set_avma(av);
1012 67842996 : }
1013 :
1014 : static void
1015 4818148 : affdbldpe(double y, dpe_t *x)
1016 : {
1017 4818148 : x->d = (double)y;
1018 4818148 : x->e = 0;
1019 4818148 : dpe_normalize(x);
1020 4818150 : }
1021 :
1022 : static void
1023 408028197 : dpe_mulz(dpe_t *x, dpe_t *y, dpe_t *z)
1024 : {
1025 408028197 : z->d = x->d * y->d;
1026 408028197 : if (z->d == 0.0)
1027 22777805 : z->e = -LONG_MAX;
1028 : else
1029 : {
1030 385250392 : z->e = x->e + y->e;
1031 385250392 : dpe_normalize0(z);
1032 : }
1033 408028567 : }
1034 :
1035 : static void
1036 54674415 : dpe_divz(dpe_t *x, dpe_t *y, dpe_t *z)
1037 : {
1038 54674415 : z->d = x->d / y->d;
1039 54674415 : if (z->d == 0.0)
1040 2392122 : z->e = -LONG_MAX;
1041 : else
1042 : {
1043 52282293 : z->e = x->e - y->e;
1044 52282293 : dpe_normalize0(z);
1045 : }
1046 54674506 : }
1047 :
1048 : static void
1049 574136 : dpe_negz(dpe_t *y, dpe_t *x)
1050 : {
1051 574136 : x->d = - y->d;
1052 574136 : x->e = y->e;
1053 574136 : }
1054 :
1055 : static void
1056 30111819 : dpe_addz(dpe_t *y, dpe_t *z, dpe_t *x)
1057 : {
1058 30111819 : if (y->e > z->e + 53)
1059 1062345 : affdpe(y, x);
1060 29049474 : else if (z->e > y->e + 53)
1061 274741 : affdpe(z, x);
1062 : else
1063 : {
1064 28774733 : long d = y->e - z->e;
1065 :
1066 28774733 : if (d >= 0)
1067 : {
1068 21361901 : x->d = y->d + ldexp(z->d, -d);
1069 21361901 : x->e = y->e;
1070 : }
1071 : else
1072 : {
1073 7412832 : x->d = z->d + ldexp(y->d, d);
1074 7412832 : x->e = z->e;
1075 : }
1076 28774733 : dpe_normalize(x);
1077 : }
1078 30111826 : }
1079 : static void
1080 423634128 : dpe_subz(dpe_t *y, dpe_t *z, dpe_t *x)
1081 : {
1082 423634128 : if (y->e > z->e + 53)
1083 38760506 : affdpe(y, x);
1084 384873622 : else if (z->e > y->e + 53)
1085 574136 : dpe_negz(z, x);
1086 : else
1087 : {
1088 384299486 : long d = y->e - z->e;
1089 :
1090 384299486 : if (d >= 0)
1091 : {
1092 351260835 : x->d = y->d - ldexp(z->d, -d);
1093 351260835 : x->e = y->e;
1094 : }
1095 : else
1096 : {
1097 33038651 : x->d = ldexp(y->d, d) - z->d;
1098 33038651 : x->e = z->e;
1099 : }
1100 384299486 : dpe_normalize(x);
1101 : }
1102 423634549 : }
1103 :
1104 : static void
1105 5912638 : dpe_muluz(dpe_t *y, ulong t, dpe_t *x)
1106 : {
1107 5912638 : x->d = y->d * (double)t;
1108 5912638 : x->e = y->e;
1109 5912638 : dpe_normalize(x);
1110 5912638 : }
1111 :
1112 : static void
1113 2301191 : dpe_addmuluz(dpe_t *y, dpe_t *z, ulong t, dpe_t *x)
1114 : {
1115 : dpe_t tmp;
1116 2301191 : dpe_muluz(z, t, &tmp);
1117 2301191 : dpe_addz(y, &tmp, x);
1118 2301191 : }
1119 :
1120 : static void
1121 2430025 : dpe_submuluz(dpe_t *y, dpe_t *z, ulong t, dpe_t *x)
1122 : {
1123 : dpe_t tmp;
1124 2430025 : dpe_muluz(z, t, &tmp);
1125 2430025 : dpe_subz(y, &tmp, x);
1126 2430025 : }
1127 :
1128 : static void
1129 393255644 : dpe_submulz(dpe_t *y, dpe_t *z, dpe_t *t, dpe_t *x)
1130 : {
1131 : dpe_t tmp;
1132 393255644 : dpe_mulz(z, t, &tmp);
1133 393255656 : dpe_subz(y, &tmp, x);
1134 393255688 : }
1135 :
1136 : static int
1137 14772724 : dpe_cmp(dpe_t *x, dpe_t *y)
1138 : {
1139 14772724 : int sx = x->d < 0. ? -1: x->d > 0.;
1140 14772724 : int sy = y->d < 0. ? -1: y->d > 0.;
1141 14772724 : int d = sx - sy;
1142 :
1143 14772724 : if (d != 0)
1144 142833 : return d;
1145 14629891 : else if (x->e > y->e)
1146 2335819 : return (sx > 0) ? 1 : -1;
1147 12294072 : else if (y->e > x->e)
1148 4319417 : return (sx > 0) ? -1 : 1;
1149 : else
1150 7974655 : return (x->d < y->d) ? -1 : (x->d > y->d);
1151 : }
1152 :
1153 : static int
1154 66774454 : dpe_abscmp(dpe_t *x, dpe_t *y)
1155 : {
1156 66774454 : if (x->e > y->e)
1157 448895 : return 1;
1158 66325559 : else if (y->e > x->e)
1159 63279381 : return -1;
1160 : else
1161 3046178 : return (fabs(x->d) < fabs(y->d)) ? -1 : (fabs(x->d) > fabs(y->d));
1162 : }
1163 :
1164 : static int
1165 9664545 : dpe_abssmall(dpe_t *x)
1166 : {
1167 9664545 : return (x->e <= 0) || (x->e == 1 && fabs(x->d) <= .75);
1168 : }
1169 :
1170 : static int
1171 14772717 : dpe_cmpmul(dpe_t *x, dpe_t *y, dpe_t *z)
1172 : {
1173 : dpe_t t;
1174 14772717 : dpe_mulz(x,y,&t);
1175 14772718 : return dpe_cmp(&t, z);
1176 : }
1177 :
1178 : static dpe_t *
1179 19236844 : cget_dpevec(long d)
1180 19236844 : { return (dpe_t*) stack_malloc_align(d*sizeof(dpe_t), sizeof(dpe_t)); }
1181 :
1182 : static dpe_t **
1183 4818146 : cget_dpemat(long d) { return (dpe_t **) cgetg(d, t_VECSMALL); }
1184 :
1185 : static GEN
1186 17373 : dpeM_diagonal_shallow(dpe_t **m, long d)
1187 : {
1188 : long i;
1189 17373 : GEN y = cgetg(d+1,t_VEC);
1190 100494 : for (i=1; i<=d; i++) gel(y, i) = dpetor(Dmael(m,i,i));
1191 17373 : return y;
1192 : }
1193 :
1194 : /************************** PROVED version (dpe) *************************/
1195 :
1196 : /* Babai's Nearest Plane algorithm (iterative).
1197 : * Size-reduces b_kappa using mu_{i,j} and r_{i,j} for j<=i <kappa
1198 : * Update B[,kappa]; compute mu_{kappa,j}, r_{kappa,j} for j<=kappa and s[kappa]
1199 : * mu, r, s updated in place (affrr). Return 1 on failure, else 0. */
1200 : static int
1201 10372356 : Babai_dpe(pari_sp av, long kappa, GEN *pG, GEN *pB, GEN *pU, dpe_t **mu, dpe_t **r, dpe_t *s,
1202 : long a, long zeros, long maxG, dpe_t *eta)
1203 : {
1204 10372356 : GEN G = *pG, B = *pB, U = *pU, ztmp;
1205 10372356 : long k, d, n, aa = a > zeros? a: zeros+1;
1206 10372356 : long emaxmu = EX0, emax2mu = EX0;
1207 : /* N.B: we set d = 0 (resp. n = 0) to avoid updating U (resp. B) */
1208 10372356 : d = U? lg(U)-1: 0;
1209 10372356 : n = B? nbrows(B): 0;
1210 2165723 : for (;;) {
1211 12538124 : int go_on = 0;
1212 12538124 : long i, j, emax3mu = emax2mu;
1213 :
1214 12538124 : if (gc_needed(av,2))
1215 : {
1216 305 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[1], a=%ld", aa);
1217 305 : gc_lll(av,3,&G,&B,&U);
1218 : }
1219 : /* Step2: compute the GSO for stage kappa */
1220 12538087 : emax2mu = emaxmu; emaxmu = EX0;
1221 67212539 : for (j=aa; j<kappa; j++)
1222 : {
1223 : dpe_t g;
1224 54674455 : affidpe(gmael(G,kappa,j), &g);
1225 388233397 : for (k = zeros+1; k < j; k++)
1226 333559119 : dpe_submulz(&g, Dmael(mu,j,k), Dmael(r,kappa,k), &g);
1227 54674278 : affdpe(&g, Dmael(r,kappa,j));
1228 54674404 : dpe_divz(Dmael(r,kappa,j), Dmael(r,j,j), Dmael(mu,kappa,j));
1229 54674415 : emaxmu = maxss(emaxmu, Dmael(mu,kappa,j)->e);
1230 : }
1231 12538084 : if (emax3mu != EX0 && emax3mu <= emax2mu + 5) /* precision too low */
1232 2 : { *pG = G; *pB = B; *pU = U; return 1; }
1233 :
1234 77146818 : for (j=kappa-1; j>zeros; j--)
1235 66774456 : if (dpe_abscmp(Dmael(mu,kappa,j), eta) > 0) { go_on = 1; break; }
1236 :
1237 : /* Step3--5: compute the X_j's */
1238 12538080 : if (go_on)
1239 23430840 : for (j=kappa-1; j>zeros; j--)
1240 : {
1241 : pari_sp btop;
1242 21265113 : dpe_t *tmp = Dmael(mu,kappa,j);
1243 21265113 : if (tmp->e < 0) continue; /* (essentially) size-reduced */
1244 :
1245 9664546 : if (gc_needed(av,2))
1246 : {
1247 1317 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[2], a=%ld, j=%ld", aa,j);
1248 1317 : gc_lll(av,3,&G,&B,&U);
1249 : }
1250 : /* we consider separately the case |X| = 1 */
1251 9664546 : if (dpe_abssmall(tmp))
1252 : {
1253 8226480 : if (tmp->d > 0) { /* in this case, X = 1 */
1254 31459342 : for (k=zeros+1; k<j; k++)
1255 27338692 : dpe_subz(Dmael(mu,kappa,k), Dmael(mu,j,k), Dmael(mu,kappa,k));
1256 42578884 : for (i=1; i<=n; i++)
1257 38458234 : gmael(B,kappa,i) = subii(gmael(B,kappa,i), gmael(B,j,i));
1258 68769151 : for (i=1; i<=d; i++)
1259 64648515 : gmael(U,kappa,i) = subii(gmael(U,kappa,i), gmael(U,j,i));
1260 4120636 : btop = avma;
1261 4120636 : ztmp = subii(gmael(G,j,j), shifti(gmael(G,kappa,j), 1));
1262 4120650 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
1263 4120646 : gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
1264 35921505 : for (i=1; i<=j; i++)
1265 31800856 : gmael(G,kappa,i) = subii(gmael(G,kappa,i), gmael(G,j,i));
1266 34521402 : for (i=j+1; i<kappa; i++)
1267 30400753 : gmael(G,kappa,i) = subii(gmael(G,kappa,i), gmael(G,i,j));
1268 27765192 : for (i=kappa+1; i<=maxG; i++)
1269 23644546 : gmael(G,i,kappa) = subii(gmael(G,i,kappa), gmael(G,i,j));
1270 : } else { /* otherwise X = -1 */
1271 31344944 : for (k=zeros+1; k<j; k++)
1272 27239112 : dpe_addz(Dmael(mu,kappa,k), Dmael(mu,j,k), Dmael(mu,kappa,k));
1273 42535001 : for (i=1; i<=n; i++)
1274 38429169 : gmael(B,kappa,i) = addii(gmael(B,kappa,i),gmael(B,j,i));
1275 68415337 : for (i=1; i<=d; i++)
1276 64309524 : gmael(U,kappa,i) = addii(gmael(U,kappa,i),gmael(U,j,i));
1277 4105813 : btop = avma;
1278 4105813 : ztmp = addii(gmael(G,j,j), shifti(gmael(G,kappa,j), 1));
1279 4105829 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
1280 4105831 : gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
1281 35727592 : for (i=1; i<=j; i++)
1282 31621763 : gmael(G,kappa,i) = addii(gmael(G,kappa,i), gmael(G,j,i));
1283 34404424 : for (i=j+1; i<kappa; i++)
1284 30298600 : gmael(G,kappa,i) = addii(gmael(G,kappa,i), gmael(G,i,j));
1285 27688372 : for (i=kappa+1; i<=maxG; i++)
1286 23582549 : gmael(G,i,kappa) = addii(gmael(G,i,kappa), gmael(G,i,j));
1287 : }
1288 8226469 : continue;
1289 : }
1290 : /* we have |X| >= 2 */
1291 1438064 : if (tmp->e < BITS_IN_LONG-1)
1292 : {
1293 1263083 : long xx = (long) pari_rint(ldexp(tmp->d, tmp->e)); /* X fits in an long */
1294 1263084 : if (xx > 0) /* = xx */
1295 : {
1296 3090006 : for (k=zeros+1; k<j; k++)
1297 2430025 : dpe_submuluz(Dmael(mu,kappa,k), Dmael(mu,j,k), xx, Dmael(mu,kappa,k));
1298 4299041 : for (i=1; i<=n; i++)
1299 3639060 : gmael(B,kappa,i) = submuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
1300 11189367 : for (i=1; i<=d; i++)
1301 10529386 : gmael(U,kappa,i) = submuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
1302 659981 : btop = avma;
1303 659981 : ztmp = submuliu2n(mulii(gmael(G,j,j), sqru(xx)), gmael(G,kappa,j), xx, 1);
1304 659983 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
1305 659983 : gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
1306 4158343 : for (i=1; i<=j; i++)
1307 3498363 : gmael(G,kappa,i) = submuliu_inplace(gmael(G,kappa,i), gmael(G,j,i), xx);
1308 6575620 : for (i=j+1; i<kappa; i++)
1309 5915638 : gmael(G,kappa,i) = submuliu_inplace(gmael(G,kappa,i), gmael(G,i,j), xx);
1310 3401484 : for (i=kappa+1; i<=maxG; i++)
1311 2741501 : gmael(G,i,kappa) = submuliu_inplace(gmael(G,i,kappa), gmael(G,i,j), xx);
1312 : }
1313 : else /* = -xx */
1314 : {
1315 603103 : xx = -xx;
1316 2904294 : for (k=zeros+1; k<j; k++)
1317 2301191 : dpe_addmuluz(Dmael(mu,kappa,k), Dmael(mu,j,k), xx, Dmael(mu,kappa,k));
1318 4233219 : for (i=1; i<=n; i++)
1319 3630116 : gmael(B,kappa,i) = addmuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
1320 10255241 : for (i=1; i<=d; i++)
1321 9652142 : gmael(U,kappa,i) = addmuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
1322 603099 : btop = avma;
1323 603099 : ztmp = addmuliu2n(mulii(gmael(G,j,j), sqru(xx)), gmael(G,kappa,j), xx, 1);
1324 603102 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
1325 603102 : gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
1326 3753591 : for (i=1; i<=j; i++)
1327 3150491 : gmael(G,kappa,i) = addmuliu_inplace(gmael(G,kappa,i), gmael(G,j,i), xx);
1328 6367387 : for (i=j+1; i<kappa; i++)
1329 5764289 : gmael(G,kappa,i) = addmuliu_inplace(gmael(G,kappa,i), gmael(G,i,j), xx);
1330 2989996 : for (i=kappa+1; i<=maxG; i++)
1331 2386899 : gmael(G,i,kappa) = addmuliu_inplace(gmael(G,i,kappa), gmael(G,i,j), xx);
1332 : }
1333 : }
1334 : else
1335 : {
1336 174981 : long e = tmp->e - BITS_IN_LONG + 1;
1337 174981 : long xx = (long) pari_rint(ldexp(tmp->d, BITS_IN_LONG - 1));
1338 175006 : if (xx > 0)
1339 : {
1340 700388 : for (k=zeros+1; k<j; k++)
1341 : {
1342 : dpe_t x;
1343 609904 : dpe_muluz(Dmael(mu,j,k), xx, &x);
1344 609904 : x.e += e;
1345 609904 : dpe_subz(Dmael(mu,kappa,k), &x, Dmael(mu,kappa,k));
1346 : }
1347 1732528 : for (i=1; i<=n; i++)
1348 1642044 : gmael(B,kappa,i) = submuliu2n(gmael(B,kappa,i), gmael(B,j,i), xx, e);
1349 463822 : for (i=1; i<=d; i++)
1350 373338 : gmael(U,kappa,i) = submuliu2n(gmael(U,kappa,i), gmael(U,j,i), xx, e);
1351 90484 : btop = avma;
1352 90484 : ztmp = submuliu2n(mulshift(gmael(G,j,j), sqrs(xx), 2*e),
1353 90484 : gmael(G,kappa,j), xx, e+1);
1354 90484 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
1355 90484 : gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
1356 791428 : for (i=1; i<=j; i++)
1357 700944 : gmael(G,kappa,i) = submuliu2n(gmael(G,kappa,i), gmael(G,j,i), xx, e);
1358 751887 : for ( ; i<kappa; i++)
1359 661403 : gmael(G,kappa,i) = submuliu2n(gmael(G,kappa,i), gmael(G,i,j), xx, e);
1360 151612 : for (i=kappa+1; i<=maxG; i++)
1361 61128 : gmael(G,i,kappa) = submuliu2n(gmael(G,i,kappa), gmael(G,i,j), xx, e);
1362 : } else
1363 : {
1364 84522 : xx = -xx;
1365 656040 : for (k=zeros+1; k<j; k++)
1366 : {
1367 : dpe_t x;
1368 571518 : dpe_muluz(Dmael(mu,j,k), xx, &x);
1369 571518 : x.e += e;
1370 571518 : dpe_addz(Dmael(mu,kappa,k), &x, Dmael(mu,kappa,k));
1371 : }
1372 1704854 : for (i=1; i<=n; i++)
1373 1620332 : gmael(B,kappa,i) = addmuliu2n(gmael(B,kappa,i), gmael(B,j,i), xx, e);
1374 378247 : for (i=1; i<=d; i++)
1375 293725 : gmael(U,kappa,i) = addmuliu2n(gmael(U,kappa,i), gmael(U,j,i), xx, e);
1376 84522 : btop = avma;
1377 84522 : ztmp = addmuliu2n(mulshift(gmael(G,j,j), sqrs(xx), 2*e),
1378 84522 : gmael(G,kappa,j), xx, e+1);
1379 84522 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
1380 84522 : gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
1381 741027 : for (i=1; i<=j; i++)
1382 656505 : gmael(G,kappa,i) = addmuliu2n(gmael(G,kappa,i), gmael(G,j,i), xx, e);
1383 719172 : for ( ; i<kappa; i++)
1384 634650 : gmael(G,kappa,i) = addmuliu2n(gmael(G,kappa,i), gmael(G,i,j), xx, e);
1385 125875 : for (i=kappa+1; i<=maxG; i++)
1386 41353 : gmael(G,i,kappa) = addmuliu2n(gmael(G,i,kappa), gmael(G,i,j), xx, e);
1387 : }
1388 : }
1389 : }
1390 12538089 : if (!go_on) break; /* Anything happened? */
1391 2165723 : aa = zeros+1;
1392 : }
1393 :
1394 10372366 : affidpe(gmael(G,kappa,kappa), Del(s,zeros+1));
1395 : /* the last s[kappa-1]=r[kappa][kappa] is computed only if kappa increases */
1396 61213763 : for (k=zeros+1; k<=kappa-2; k++)
1397 50841424 : dpe_submulz(Del(s,k), Dmael(mu,kappa,k), Dmael(r,kappa,k), Del(s,k+1));
1398 10372339 : *pG = G; *pB = B; *pU = U; return 0;
1399 : }
1400 :
1401 : /* G integral Gram matrix, LLL-reduces (G,B,U) in place [apply base change
1402 : * transforms to B and U]. If (keepfirst), never swap with first vector.
1403 : * If G = NULL, we compute the Gram matrix incrementally.
1404 : * Return -1 on failure, else zeros = dim Kernel (>= 0) */
1405 : static long
1406 2409082 : fplll_dpe(GEN *pG, GEN *pB, GEN *pU, GEN *pr, double DELTA, double ETA,
1407 : long keepfirst)
1408 : {
1409 : pari_sp av;
1410 2409082 : GEN Gtmp, alpha, G = *pG, B = *pB, U = *pU;
1411 2409082 : long d, maxG, kappa, kappa2, i, j, zeros, kappamax, incgram = !G, cnt = 0;
1412 : dpe_t delta, eta, **mu, **r, *s;
1413 2409082 : affdbldpe(DELTA,&delta);
1414 2409078 : affdbldpe(ETA,&eta);
1415 :
1416 2409080 : if (incgram)
1417 : { /* incremental Gram matrix */
1418 2327020 : maxG = 2; d = lg(B)-1;
1419 2327020 : G = zeromatcopy(d, d);
1420 : }
1421 : else
1422 82060 : maxG = d = lg(G)-1;
1423 :
1424 2409084 : mu = cget_dpemat(d+1);
1425 2409078 : r = cget_dpemat(d+1);
1426 2409083 : s = cget_dpevec(d+1);
1427 10823045 : for (j = 1; j <= d; j++)
1428 : {
1429 8413960 : mu[j]= cget_dpevec(d+1);
1430 8413961 : r[j] = cget_dpevec(d+1);
1431 : }
1432 2409085 : Gtmp = cgetg(d+1, t_VEC);
1433 2409079 : alpha = cgetg(d+1, t_VECSMALL);
1434 2409068 : av = avma;
1435 :
1436 : /* Step2: Initializing the main loop */
1437 2409068 : kappamax = 1;
1438 2409068 : i = 1;
1439 : do {
1440 2762191 : if (incgram) gmael(G,i,i) = ZV_dotsquare(gel(B,i));
1441 2762097 : affidpe(gmael(G,i,i), Dmael(r,i,i));
1442 2762112 : } while (!signe(gmael(G,i,i)) && ++i <= d);
1443 2408989 : zeros = i-1; /* all basis vectors b_i with i <= zeros are zero vectors */
1444 2408989 : kappa = i;
1445 10469611 : for (i=zeros+1; i<=d; i++) alpha[i]=1;
1446 :
1447 12781383 : while (++kappa <= d)
1448 : {
1449 10372364 : if (kappa > kappamax)
1450 : {
1451 5651634 : if (DEBUGLEVEL>=4) err_printf("K%ld ",kappa);
1452 5651634 : kappamax = kappa;
1453 5651634 : if (incgram)
1454 : {
1455 24749675 : for (i=zeros+1; i<=kappa; i++)
1456 19328856 : gmael(G,kappa,i) = ZV_dotproduct(gel(B,kappa), gel(B,i));
1457 5420819 : maxG = kappamax;
1458 : }
1459 : }
1460 : /* Step3: Call to the Babai algorithm, mu,r,s updated in place */
1461 10372014 : if (Babai_dpe(av, kappa, &G,&B,&U, mu,r,s, alpha[kappa], zeros, maxG, &eta))
1462 2 : { *pG = incgram? NULL: G; *pB = B; *pU = U; return -1; }
1463 20646754 : if ((keepfirst && kappa == 2) ||
1464 10274377 : dpe_cmpmul(Dmael(r,kappa-1,kappa-1), &delta, Del(s,kappa-1)) <= 0)
1465 : { /* Step4: Success of Lovasz's condition */
1466 8855412 : alpha[kappa] = kappa;
1467 8855412 : dpe_submulz(Del(s,kappa-1), Dmael(mu,kappa,kappa-1), Dmael(r,kappa,kappa-1), Dmael(r,kappa,kappa));
1468 8855411 : continue;
1469 : }
1470 : /* Step5: Find the right insertion index kappa, kappa2 = initial kappa */
1471 1516965 : if (DEBUGLEVEL>=4 && kappa==kappamax && Del(s,kappa-1)->d)
1472 0 : if (++cnt > 20) { cnt = 0; err_printf("(%ld) ", Del(s,1)->e-1); }
1473 1516965 : kappa2 = kappa;
1474 : do {
1475 4755902 : kappa--;
1476 4755902 : if (kappa < zeros+2 + (keepfirst ? 1: 0)) break;
1477 4498324 : } while (dpe_cmpmul(Dmael(r,kappa-1,kappa-1), &delta, Del(s,kappa-1)) >= 0);
1478 1516965 : update_alpha(alpha, kappa, kappa2, kappamax);
1479 :
1480 : /* Step6: Update the mu's and r's */
1481 1516965 : dperotate(mu, kappa2, kappa);
1482 1516965 : dperotate(r, kappa2, kappa);
1483 1516965 : affdpe(Del(s,kappa), Dmael(r,kappa,kappa));
1484 :
1485 : /* Step7: Update G, B, U */
1486 1516965 : if (U) rotate(U, kappa2, kappa);
1487 1516965 : if (B) rotate(B, kappa2, kappa);
1488 1516965 : rotateG(G,kappa2,kappa, maxG, Gtmp);
1489 :
1490 : /* Step8: Prepare the next loop iteration */
1491 1516965 : if (kappa == zeros+1 && !signe(gmael(G,kappa,kappa)))
1492 : {
1493 35173 : zeros++; kappa++;
1494 35173 : affidpe(gmael(G,kappa,kappa), Dmael(r,kappa,kappa));
1495 : }
1496 : }
1497 2409019 : if (pr) *pr = dpeM_diagonal_shallow(r,d);
1498 2409019 : *pG = G; *pB = B; *pU = U; return zeros; /* success */
1499 : }
1500 :
1501 :
1502 : /************************** PROVED version (t_INT) *************************/
1503 :
1504 : /* Babai's Nearest Plane algorithm (iterative).
1505 : * Size-reduces b_kappa using mu_{i,j} and r_{i,j} for j<=i <kappa
1506 : * Update B[,kappa]; compute mu_{kappa,j}, r_{kappa,j} for j<=kappa and s[kappa]
1507 : * mu, r, s updated in place (affrr). Return 1 on failure, else 0. */
1508 : static int
1509 2208 : Babai(pari_sp av, long kappa, GEN *pG, GEN *pB, GEN *pU, GEN mu, GEN r, GEN s,
1510 : long a, long zeros, long maxG, GEN eta, long prec)
1511 : {
1512 2208 : GEN G = *pG, B = *pB, U = *pU, ztmp;
1513 2208 : long k, aa = a > zeros? a: zeros+1;
1514 2208 : const long n = B? nbrows(B): 0, d = U ? lg(U)-1: 0, bit = prec2nbits(prec);
1515 2208 : long emaxmu = EX0, emax2mu = EX0;
1516 : /* N.B: we set d = 0 (resp. n = 0) to avoid updating U (resp. B) */
1517 :
1518 1262 : for (;;) {
1519 3470 : int go_on = 0;
1520 3470 : long i, j, emax3mu = emax2mu;
1521 :
1522 3470 : if (gc_needed(av,2))
1523 : {
1524 0 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[1], a=%ld", aa);
1525 0 : gc_lll(av,3,&G,&B,&U);
1526 : }
1527 : /* Step2: compute the GSO for stage kappa */
1528 3470 : emax2mu = emaxmu; emaxmu = EX0;
1529 24881 : for (j=aa; j<kappa; j++)
1530 : {
1531 21411 : pari_sp btop = avma;
1532 21411 : GEN g = gmael(G,kappa,j);
1533 209529 : for (k = zeros+1; k < j; k++)
1534 188118 : g = mpsub(g, mulrr(gmael(mu,j,k), gmael(r,kappa,k)));
1535 21411 : mpaff(g, gmael(r,kappa,j));
1536 21411 : affrr(divrr(gmael(r,kappa,j), gmael(r,j,j)), gmael(mu,kappa,j));
1537 21411 : emaxmu = maxss(emaxmu, expo(gmael(mu,kappa,j)));
1538 21411 : set_avma(btop);
1539 : }
1540 3470 : if (emax3mu != EX0 && emax3mu <= emax2mu + 5) /* precision too low */
1541 0 : { *pG = G; *pB = B; *pU = U; return 1; }
1542 :
1543 25500 : for (j=kappa-1; j>zeros; j--)
1544 23292 : if (abscmprr(gmael(mu,kappa,j), eta) > 0) { go_on = 1; break; }
1545 :
1546 : /* Step3--5: compute the X_j's */
1547 3470 : if (go_on)
1548 15993 : for (j=kappa-1; j>zeros; j--)
1549 : {
1550 : pari_sp btop;
1551 14731 : GEN tmp = gmael(mu,kappa,j);
1552 14731 : if (absrsmall(tmp)) continue; /* size-reduced */
1553 :
1554 8741 : if (gc_needed(av,2))
1555 : {
1556 0 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[2], a=%ld, j=%ld", aa,j);
1557 0 : gc_lll(av,3,&G,&B,&U);
1558 : }
1559 8741 : btop = avma;
1560 : /* we consider separately the case |X| = 1 */
1561 8741 : if (absrsmall2(tmp))
1562 : {
1563 3721 : if (signe(tmp) > 0) { /* in this case, X = 1 */
1564 17772 : for (k=zeros+1; k<j; k++)
1565 15910 : affrr(subrr(gmael(mu,kappa,k), gmael(mu,j,k)), gmael(mu,kappa,k));
1566 1862 : set_avma(btop);
1567 45875 : for (i=1; i<=n; i++)
1568 44013 : gmael(B,kappa,i) = subii(gmael(B,kappa,i), gmael(B,j,i));
1569 7955 : for (i=1; i<=d; i++)
1570 6093 : gmael(U,kappa,i) = subii(gmael(U,kappa,i), gmael(U,j,i));
1571 1862 : btop = avma;
1572 1862 : ztmp = subii(gmael(G,j,j), shifti(gmael(G,kappa,j), 1));
1573 1862 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
1574 1862 : gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
1575 19634 : for (i=1; i<=j; i++)
1576 17772 : gmael(G,kappa,i) = subii(gmael(G,kappa,i), gmael(G,j,i));
1577 17730 : for (i=j+1; i<kappa; i++)
1578 15868 : gmael(G,kappa,i) = subii(gmael(G,kappa,i), gmael(G,i,j));
1579 7280 : for (i=kappa+1; i<=maxG; i++)
1580 5418 : gmael(G,i,kappa) = subii(gmael(G,i,kappa), gmael(G,i,j));
1581 : } else { /* otherwise X = -1 */
1582 18073 : for (k=zeros+1; k<j; k++)
1583 16214 : affrr(addrr(gmael(mu,kappa,k), gmael(mu,j,k)), gmael(mu,kappa,k));
1584 1859 : set_avma(btop);
1585 45822 : for (i=1; i<=n; i++)
1586 43963 : gmael(B,kappa,i) = addii(gmael(B,kappa,i),gmael(B,j,i));
1587 7934 : for (i=1; i<=d; i++)
1588 6075 : gmael(U,kappa,i) = addii(gmael(U,kappa,i),gmael(U,j,i));
1589 1859 : btop = avma;
1590 1859 : ztmp = addii(gmael(G,j,j), shifti(gmael(G,kappa,j), 1));
1591 1859 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
1592 1859 : gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
1593 19932 : for (i=1; i<=j; i++)
1594 18073 : gmael(G,kappa,i) = addii(gmael(G,kappa,i), gmael(G,j,i));
1595 17334 : for (i=j+1; i<kappa; i++)
1596 15475 : gmael(G,kappa,i) = addii(gmael(G,kappa,i), gmael(G,i,j));
1597 7296 : for (i=kappa+1; i<=maxG; i++)
1598 5437 : gmael(G,i,kappa) = addii(gmael(G,i,kappa), gmael(G,i,j));
1599 : }
1600 3721 : continue;
1601 : }
1602 : /* we have |X| >= 2 */
1603 5020 : if (expo(tmp) < BITS_IN_LONG)
1604 : {
1605 4500 : ulong xx = roundr_safe(tmp)[2]; /* X fits in an ulong */
1606 4500 : if (signe(tmp) > 0) /* = xx */
1607 : {
1608 25486 : for (k=zeros+1; k<j; k++)
1609 23274 : affrr(subrr(gmael(mu,kappa,k), mulur(xx, gmael(mu,j,k))),
1610 23274 : gmael(mu,kappa,k));
1611 2212 : set_avma(btop);
1612 68948 : for (i=1; i<=n; i++)
1613 66736 : gmael(B,kappa,i) = submuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
1614 3796 : for (i=1; i<=d; i++)
1615 1584 : gmael(U,kappa,i) = submuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
1616 2212 : btop = avma;
1617 2212 : ztmp = submuliu2n(mulii(gmael(G,j,j), sqru(xx)), gmael(G,kappa,j), xx, 1);
1618 2212 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
1619 2212 : gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
1620 27698 : for (i=1; i<=j; i++)
1621 25486 : gmael(G,kappa,i) = submuliu_inplace(gmael(G,kappa,i), gmael(G,j,i), xx);
1622 30761 : for (i=j+1; i<kappa; i++)
1623 28549 : gmael(G,kappa,i) = submuliu_inplace(gmael(G,kappa,i), gmael(G,i,j), xx);
1624 6073 : for (i=kappa+1; i<=maxG; i++)
1625 3861 : gmael(G,i,kappa) = submuliu_inplace(gmael(G,i,kappa), gmael(G,i,j), xx);
1626 : }
1627 : else /* = -xx */
1628 : {
1629 26514 : for (k=zeros+1; k<j; k++)
1630 24226 : affrr(addrr(gmael(mu,kappa,k), mulur(xx, gmael(mu,j,k))),
1631 24226 : gmael(mu,kappa,k));
1632 2288 : set_avma(btop);
1633 71318 : for (i=1; i<=n; i++)
1634 69030 : gmael(B,kappa,i) = addmuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
1635 3926 : for (i=1; i<=d; i++)
1636 1638 : gmael(U,kappa,i) = addmuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
1637 2288 : btop = avma;
1638 2288 : ztmp = addmuliu2n(mulii(gmael(G,j,j), sqru(xx)), gmael(G,kappa,j), xx, 1);
1639 2288 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
1640 2288 : gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
1641 28802 : for (i=1; i<=j; i++)
1642 26514 : gmael(G,kappa,i) = addmuliu_inplace(gmael(G,kappa,i), gmael(G,j,i), xx);
1643 31397 : for (i=j+1; i<kappa; i++)
1644 29109 : gmael(G,kappa,i) = addmuliu_inplace(gmael(G,kappa,i), gmael(G,i,j), xx);
1645 6391 : for (i=kappa+1; i<=maxG; i++)
1646 4103 : gmael(G,i,kappa) = addmuliu_inplace(gmael(G,i,kappa), gmael(G,i,j), xx);
1647 : }
1648 : }
1649 : else
1650 : {
1651 : long e;
1652 520 : GEN X = truncexpo(tmp, bit, &e); /* tmp ~ X * 2^e */
1653 520 : btop = avma;
1654 6849 : for (k=zeros+1; k<j; k++)
1655 : {
1656 6329 : GEN x = mulir(X, gmael(mu,j,k));
1657 6329 : if (e) shiftr_inplace(x, e);
1658 6329 : affrr(subrr(gmael(mu,kappa,k), x), gmael(mu,kappa,k));
1659 : }
1660 520 : set_avma(btop);
1661 16677 : for (i=1; i<=n; i++)
1662 16157 : gmael(B,kappa,i) = submulshift(gmael(B,kappa,i), gmael(B,j,i), X, e);
1663 709 : for (i=1; i<=d; i++)
1664 189 : gmael(U,kappa,i) = submulshift(gmael(U,kappa,i), gmael(U,j,i), X, e);
1665 520 : btop = avma;
1666 520 : ztmp = submulshift(mulshift(gmael(G,j,j), sqri(X), 2*e),
1667 520 : gmael(G,kappa,j), X, e+1);
1668 520 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
1669 520 : gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
1670 7369 : for (i=1; i<=j; i++)
1671 6849 : gmael(G,kappa,i) = submulshift(gmael(G,kappa,i), gmael(G,j,i), X, e);
1672 7435 : for ( ; i<kappa; i++)
1673 6915 : gmael(G,kappa,i) = submulshift(gmael(G,kappa,i), gmael(G,i,j), X, e);
1674 580 : for (i=kappa+1; i<=maxG; i++)
1675 60 : gmael(G,i,kappa) = submulshift(gmael(G,i,kappa), gmael(G,i,j), X, e);
1676 : }
1677 : }
1678 3470 : if (!go_on) break; /* Anything happened? */
1679 1262 : aa = zeros+1;
1680 : }
1681 :
1682 2208 : affir(gmael(G,kappa,kappa), gel(s,zeros+1));
1683 : /* the last s[kappa-1]=r[kappa][kappa] is computed only if kappa increases */
1684 2208 : av = avma;
1685 20786 : for (k=zeros+1; k<=kappa-2; k++)
1686 18578 : affrr(subrr(gel(s,k), mulrr(gmael(mu,kappa,k), gmael(r,kappa,k))),
1687 18578 : gel(s,k+1));
1688 2208 : *pG = G; *pB = B; *pU = U; return gc_bool(av, 0);
1689 : }
1690 :
1691 : /* G integral Gram matrix, LLL-reduces (G,B,U) in place [apply base change
1692 : * transforms to B and U]. If (keepfirst), never swap with first vector.
1693 : * If G = NULL, we compute the Gram matrix incrementally.
1694 : * Return -1 on failure, else zeros = dim Kernel (>= 0) */
1695 : static long
1696 2 : fplll(GEN *pG, GEN *pB, GEN *pU, GEN *pr, double DELTA, double ETA,
1697 : long keepfirst, long prec)
1698 : {
1699 : pari_sp av, av2;
1700 2 : GEN mu, r, s, tmp, Gtmp, alpha, G = *pG, B = *pB, U = *pU;
1701 2 : GEN delta = dbltor(DELTA), eta = dbltor(ETA);
1702 2 : long d, maxG, kappa, kappa2, i, j, zeros, kappamax, incgram = !G, cnt = 0;
1703 :
1704 2 : if (incgram)
1705 : { /* incremental Gram matrix */
1706 2 : maxG = 2; d = lg(B)-1;
1707 2 : G = zeromatcopy(d, d);
1708 : }
1709 : else
1710 0 : maxG = d = lg(G)-1;
1711 :
1712 2 : mu = cgetg(d+1, t_MAT);
1713 2 : r = cgetg(d+1, t_MAT);
1714 2 : s = cgetg(d+1, t_VEC);
1715 43 : for (j = 1; j <= d; j++)
1716 : {
1717 41 : GEN M = cgetg(d+1, t_COL), R = cgetg(d+1, t_COL);
1718 41 : gel(mu,j)= M;
1719 41 : gel(r,j) = R;
1720 41 : gel(s,j) = cgetr(prec);
1721 1146 : for (i = 1; i <= d; i++)
1722 : {
1723 1105 : gel(R,i) = cgetr(prec);
1724 1105 : gel(M,i) = cgetr(prec);
1725 : }
1726 : }
1727 2 : Gtmp = cgetg(d+1, t_VEC);
1728 2 : alpha = cgetg(d+1, t_VECSMALL);
1729 2 : av = avma;
1730 :
1731 : /* Step2: Initializing the main loop */
1732 2 : kappamax = 1;
1733 2 : i = 1;
1734 : do {
1735 2 : if (incgram) gmael(G,i,i) = ZV_dotsquare(gel(B,i));
1736 2 : affir(gmael(G,i,i), gmael(r,i,i));
1737 2 : } while (!signe(gmael(G,i,i)) && ++i <= d);
1738 2 : zeros = i-1; /* all basis vectors b_i with i <= zeros are zero vectors */
1739 2 : kappa = i;
1740 43 : for (i=zeros+1; i<=d; i++) alpha[i]=1;
1741 :
1742 2210 : while (++kappa <= d)
1743 : {
1744 2208 : if (kappa > kappamax)
1745 : {
1746 39 : if (DEBUGLEVEL>=4) err_printf("K%ld ",kappa);
1747 39 : kappamax = kappa;
1748 39 : if (incgram)
1749 : {
1750 610 : for (i=zeros+1; i<=kappa; i++)
1751 571 : gmael(G,kappa,i) = ZV_dotproduct(gel(B,kappa), gel(B,i));
1752 39 : maxG = kappamax;
1753 : }
1754 : }
1755 : /* Step3: Call to the Babai algorithm, mu,r,s updated in place */
1756 2208 : if (Babai(av, kappa, &G,&B,&U, mu,r,s, alpha[kappa], zeros, maxG, eta, prec))
1757 0 : { *pG = incgram? NULL: G; *pB = B; *pU = U; return -1; }
1758 2208 : av2 = avma;
1759 4416 : if ((keepfirst && kappa == 2) ||
1760 2208 : cmprr(mulrr(gmael(r,kappa-1,kappa-1), delta), gel(s,kappa-1)) <= 0)
1761 : { /* Step4: Success of Lovasz's condition */
1762 1464 : alpha[kappa] = kappa;
1763 1464 : tmp = mulrr(gmael(mu,kappa,kappa-1), gmael(r,kappa,kappa-1));
1764 1464 : affrr(subrr(gel(s,kappa-1), tmp), gmael(r,kappa,kappa));
1765 1464 : set_avma(av2); continue;
1766 : }
1767 : /* Step5: Find the right insertion index kappa, kappa2 = initial kappa */
1768 744 : if (DEBUGLEVEL>=4 && kappa==kappamax && signe(gel(s,kappa-1)))
1769 0 : if (++cnt > 20) { cnt = 0; err_printf("(%ld) ", expo(gel(s,1))); }
1770 744 : kappa2 = kappa;
1771 : do {
1772 2169 : kappa--;
1773 2169 : if (kappa < zeros+2 + (keepfirst ? 1: 0)) break;
1774 2033 : tmp = mulrr(gmael(r,kappa-1,kappa-1), delta);
1775 2033 : } while (cmprr(gel(s,kappa-1), tmp) <= 0);
1776 744 : set_avma(av2);
1777 744 : update_alpha(alpha, kappa, kappa2, kappamax);
1778 :
1779 : /* Step6: Update the mu's and r's */
1780 744 : rotate(mu, kappa2, kappa);
1781 744 : rotate(r, kappa2, kappa);
1782 744 : affrr(gel(s,kappa), gmael(r,kappa,kappa));
1783 :
1784 : /* Step7: Update G, B, U */
1785 744 : if (U) rotate(U, kappa2, kappa);
1786 744 : if (B) rotate(B, kappa2, kappa);
1787 744 : rotateG(G,kappa2,kappa, maxG, Gtmp);
1788 :
1789 : /* Step8: Prepare the next loop iteration */
1790 744 : if (kappa == zeros+1 && !signe(gmael(G,kappa,kappa)))
1791 : {
1792 0 : zeros++; kappa++;
1793 0 : affir(gmael(G,kappa,kappa), gmael(r,kappa,kappa));
1794 : }
1795 : }
1796 2 : if (pr) *pr = RgM_diagonal_shallow(r);
1797 2 : *pG = G; *pB = B; *pU = U; return zeros; /* success */
1798 : }
1799 :
1800 : /* Assume x a ZM, if pN != NULL, set it to Gram-Schmidt (squared) norms
1801 : * The following modes are supported:
1802 : * - flag & LLL_INPLACE: x a lattice basis, return x*U
1803 : * - flag & LLL_GRAM: x a Gram matrix / else x a lattice basis; return
1804 : * LLL base change matrix U [LLL_IM]
1805 : * kernel basis [LLL_KER, nonreduced]
1806 : * both [LLL_ALL] */
1807 : GEN
1808 2531887 : ZM_lll_norms(GEN x, double DELTA, long flag, GEN *pN)
1809 : {
1810 2531887 : pari_sp av = avma;
1811 2531887 : const double ETA = 0.51;
1812 2531887 : const long keepfirst = flag & LLL_KEEP_FIRST;
1813 2531887 : long p, zeros = -1, n = lg(x)-1;
1814 : GEN G, B, U;
1815 : pari_timer T;
1816 2531887 : if (n <= 1) return lll_trivial(x, flag);
1817 2422859 : if (nbrows(x)==0)
1818 : {
1819 13861 : if (flag & LLL_KER) return matid(n);
1820 13861 : if (flag & (LLL_INPLACE|LLL_IM)) return cgetg(1,t_MAT);
1821 0 : retmkvec2(matid(n), cgetg(1,t_MAT));
1822 : }
1823 2409045 : x = RgM_shallowcopy(x);
1824 2409056 : if (flag & LLL_GRAM)
1825 82040 : { G = x; B = NULL; U = matid(n); }
1826 : else
1827 2327016 : { G = NULL; B = x; U = (flag & LLL_INPLACE)? NULL: matid(n); }
1828 2409070 : if(DEBUGLEVEL>=4) timer_start(&T);
1829 2409070 : if (!(flag & LLL_GRAM))
1830 : {
1831 : long t;
1832 2327011 : if(DEBUGLEVEL>=4)
1833 0 : err_printf("Entering L^2 (double): LLL-parameters (%.3f,%.3f)\n",
1834 : DELTA,ETA);
1835 2327011 : zeros = fplll_fast(&B, &U, DELTA, ETA, keepfirst);
1836 2327027 : if (DEBUGLEVEL>=4) timer_printf(&T, zeros < 0? "LLL (failed)": "LLL");
1837 2330472 : for (p = DEFAULTPREC, t = 0; zeros < 0 && t < 1 ; p += EXTRAPREC64, t++)
1838 : {
1839 3447 : if (DEBUGLEVEL>=4)
1840 0 : err_printf("Entering L^2 (heuristic): LLL-parameters (%.3f,%.3f), prec = %d/%d\n", DELTA, ETA, p, p);
1841 3447 : zeros = fplll_heuristic(&B, &U, DELTA, ETA, keepfirst, p, p);
1842 3447 : gc_lll(av, 2, &B, &U);
1843 3447 : if (DEBUGLEVEL>=4) timer_printf(&T, zeros < 0? "LLL (failed)": "LLL");
1844 : }
1845 : }
1846 2409084 : if(DEBUGLEVEL>=4)
1847 0 : err_printf("Entering L^2 (dpe): LLL-parameters (%.3f,%.3f)\n", DELTA,ETA);
1848 2409084 : zeros = fplll_dpe(&G, &B, &U, pN, DELTA, ETA, keepfirst);
1849 2409022 : if (DEBUGLEVEL>=4) timer_printf(&T, zeros < 0? "LLL (failed)": "LLL");
1850 2409048 : if (zeros < 0)
1851 2 : for (p = DEFAULTPREC;; p += EXTRAPREC64)
1852 : {
1853 2 : if (DEBUGLEVEL>=4)
1854 0 : err_printf("Entering L^2: LLL-parameters (%.3f,%.3f), prec = %d\n",
1855 : DELTA,ETA, p);
1856 2 : zeros = fplll(&G, &B, &U, pN, DELTA, ETA, keepfirst, p);
1857 2 : if (DEBUGLEVEL>=4) timer_printf(&T, zeros < 0? "LLL (failed)": "LLL");
1858 2 : if (zeros >= 0) break;
1859 0 : gc_lll(av, 3, &G, &B, &U);
1860 : }
1861 2409048 : return lll_finish(U? U: B, zeros, flag);
1862 : }
1863 :
1864 : /********************************************************************/
1865 : /** **/
1866 : /** LLL OVER K[X] **/
1867 : /** **/
1868 : /********************************************************************/
1869 : static int
1870 378 : pslg(GEN x)
1871 : {
1872 : long tx;
1873 378 : if (gequal0(x)) return 2;
1874 336 : tx = typ(x); return is_scalar_t(tx)? 3: lg(x);
1875 : }
1876 :
1877 : static int
1878 147 : REDgen(long k, long l, GEN h, GEN L, GEN B)
1879 : {
1880 147 : GEN q, u = gcoeff(L,k,l);
1881 : long i;
1882 :
1883 147 : if (pslg(u) < pslg(B)) return 0;
1884 :
1885 105 : q = gneg(gdeuc(u,B));
1886 105 : gel(h,k) = gadd(gel(h,k), gmul(q,gel(h,l)));
1887 105 : for (i=1; i<l; i++) gcoeff(L,k,i) = gadd(gcoeff(L,k,i), gmul(q,gcoeff(L,l,i)));
1888 105 : gcoeff(L,k,l) = gadd(gcoeff(L,k,l), gmul(q,B)); return 1;
1889 : }
1890 :
1891 : static int
1892 147 : do_SWAPgen(GEN h, GEN L, GEN B, long k, GEN fl, int *flc)
1893 : {
1894 : GEN p1, la, la2, Bk;
1895 : long ps1, ps2, i, j, lx;
1896 :
1897 147 : if (!fl[k-1]) return 0;
1898 :
1899 105 : la = gcoeff(L,k,k-1); la2 = gsqr(la);
1900 105 : Bk = gel(B,k);
1901 105 : if (fl[k])
1902 : {
1903 42 : GEN q = gadd(la2, gmul(gel(B,k-1),gel(B,k+1)));
1904 42 : ps1 = pslg(gsqr(Bk));
1905 42 : ps2 = pslg(q);
1906 42 : if (ps1 <= ps2 && (ps1 < ps2 || !*flc)) return 0;
1907 21 : *flc = (ps1 != ps2);
1908 21 : gel(B,k) = gdiv(q, Bk);
1909 : }
1910 :
1911 84 : swap(gel(h,k-1), gel(h,k)); lx = lg(L);
1912 84 : for (j=1; j<k-1; j++) swap(gcoeff(L,k-1,j), gcoeff(L,k,j));
1913 84 : if (fl[k])
1914 : {
1915 21 : for (i=k+1; i<lx; i++)
1916 : {
1917 0 : GEN t = gcoeff(L,i,k);
1918 0 : p1 = gsub(gmul(gel(B,k+1),gcoeff(L,i,k-1)), gmul(la,t));
1919 0 : gcoeff(L,i,k) = gdiv(p1, Bk);
1920 0 : p1 = gadd(gmul(la,gcoeff(L,i,k-1)), gmul(gel(B,k-1),t));
1921 0 : gcoeff(L,i,k-1) = gdiv(p1, Bk);
1922 : }
1923 : }
1924 63 : else if (!gequal0(la))
1925 : {
1926 21 : p1 = gdiv(la2, Bk);
1927 21 : gel(B,k+1) = gel(B,k) = p1;
1928 21 : for (i=k+2; i<=lx; i++) gel(B,i) = gdiv(gmul(p1,gel(B,i)),Bk);
1929 21 : for (i=k+1; i<lx; i++)
1930 0 : gcoeff(L,i,k-1) = gdiv(gmul(la,gcoeff(L,i,k-1)), Bk);
1931 21 : for (j=k+1; j<lx-1; j++)
1932 0 : for (i=j+1; i<lx; i++)
1933 0 : gcoeff(L,i,j) = gdiv(gmul(p1,gcoeff(L,i,j)), Bk);
1934 : }
1935 : else
1936 : {
1937 42 : gcoeff(L,k,k-1) = gen_0;
1938 42 : for (i=k+1; i<lx; i++)
1939 : {
1940 0 : gcoeff(L,i,k) = gcoeff(L,i,k-1);
1941 0 : gcoeff(L,i,k-1) = gen_0;
1942 : }
1943 42 : gel(B,k) = gel(B,k-1); fl[k] = 1; fl[k-1] = 0;
1944 : }
1945 84 : return 1;
1946 : }
1947 :
1948 : static void
1949 126 : incrementalGSgen(GEN x, GEN L, GEN B, long k, GEN fl)
1950 : {
1951 126 : GEN u = NULL; /* gcc -Wall */
1952 : long i, j;
1953 315 : for (j = 1; j <= k; j++)
1954 189 : if (j==k || fl[j])
1955 : {
1956 189 : u = gcoeff(x,k,j);
1957 189 : if (!is_extscalar_t(typ(u))) pari_err_TYPE("incrementalGSgen",u);
1958 252 : for (i=1; i<j; i++)
1959 63 : if (fl[i])
1960 : {
1961 63 : u = gsub(gmul(gel(B,i+1),u), gmul(gcoeff(L,k,i),gcoeff(L,j,i)));
1962 63 : u = gdiv(u, gel(B,i));
1963 : }
1964 189 : gcoeff(L,k,j) = u;
1965 : }
1966 126 : if (gequal0(u)) gel(B,k+1) = gel(B,k);
1967 : else
1968 : {
1969 84 : gel(B,k+1) = gcoeff(L,k,k); gcoeff(L,k,k) = gen_1; fl[k] = 1;
1970 : }
1971 126 : }
1972 :
1973 : static GEN
1974 126 : lllgramallgen(GEN x, long flag)
1975 : {
1976 126 : long lx = lg(x), i, j, k, l, n;
1977 : pari_sp av;
1978 : GEN B, L, h, fl;
1979 : int flc;
1980 :
1981 126 : n = lx-1; if (n<=1) return lll_trivial(x,flag);
1982 63 : if (lgcols(x) != lx) pari_err_DIM("lllgramallgen");
1983 :
1984 63 : fl = cgetg(lx, t_VECSMALL);
1985 :
1986 63 : av = avma;
1987 63 : B = scalarcol_shallow(gen_1, lx);
1988 63 : L = cgetg(lx,t_MAT);
1989 189 : for (j=1; j<lx; j++) { gel(L,j) = zerocol(n); fl[j] = 0; }
1990 :
1991 63 : h = matid(n);
1992 189 : for (i=1; i<lx; i++)
1993 126 : incrementalGSgen(x, L, B, i, fl);
1994 63 : flc = 0;
1995 63 : for(k=2;;)
1996 : {
1997 147 : if (REDgen(k, k-1, h, L, gel(B,k))) flc = 1;
1998 147 : if (do_SWAPgen(h, L, B, k, fl, &flc)) { if (k > 2) k--; }
1999 : else
2000 : {
2001 63 : for (l=k-2; l>=1; l--)
2002 0 : if (REDgen(k, l, h, L, gel(B,l+1))) flc = 1;
2003 63 : if (++k > n) break;
2004 : }
2005 84 : if (gc_needed(av,1))
2006 : {
2007 0 : if(DEBUGMEM>1) pari_warn(warnmem,"lllgramallgen");
2008 0 : gerepileall(av,3,&B,&L,&h);
2009 : }
2010 : }
2011 105 : k=1; while (k<lx && !fl[k]) k++;
2012 63 : return lll_finish(h,k-1,flag);
2013 : }
2014 :
2015 : static int
2016 57163 : RgM_square(GEN x) { long l = lg(x); return l == 1 || l == lgcols(x); }
2017 : static GEN
2018 126 : lllallgen(GEN x, long flag)
2019 : {
2020 126 : pari_sp av = avma;
2021 126 : if (!(flag & LLL_GRAM)) x = gram_matrix(x);
2022 42 : else if (!RgM_square(x)) pari_err_DIM("qflllgram");
2023 126 : return gerepilecopy(av, lllgramallgen(x, flag));
2024 : }
2025 : GEN
2026 42 : lllgen(GEN x) { return lllallgen(x, LLL_IM); }
2027 : GEN
2028 42 : lllkerimgen(GEN x) { return lllallgen(x, LLL_ALL); }
2029 : GEN
2030 0 : lllgramgen(GEN x) { return lllallgen(x, LLL_IM|LLL_GRAM); }
2031 : GEN
2032 42 : lllgramkerimgen(GEN x) { return lllallgen(x, LLL_ALL|LLL_GRAM); }
2033 :
2034 : static GEN
2035 55922 : lllall(GEN x, long flag)
2036 55922 : { pari_sp av = avma; return gerepilecopy(av, ZM_lll(x, LLLDFT, flag)); }
2037 : GEN
2038 14662 : lllint(GEN x) { return lllall(x, LLL_IM); }
2039 : GEN
2040 35 : lllkerim(GEN x) { return lllall(x, LLL_ALL); }
2041 : GEN
2042 41183 : lllgramint(GEN x)
2043 41183 : { if (!RgM_square(x)) pari_err_DIM("qflllgram");
2044 41183 : return lllall(x, LLL_IM | LLL_GRAM); }
2045 : GEN
2046 35 : lllgramkerim(GEN x)
2047 35 : { if (!RgM_square(x)) pari_err_DIM("qflllgram");
2048 35 : return lllall(x, LLL_ALL | LLL_GRAM); }
2049 :
2050 : GEN
2051 645602 : lllfp(GEN x, double D, long flag)
2052 : {
2053 645602 : long n = lg(x)-1;
2054 645602 : pari_sp av = avma;
2055 : GEN h;
2056 645602 : if (n <= 1) return lll_trivial(x,flag);
2057 561811 : if ((flag & LLL_GRAM) && !RgM_square(x)) pari_err_DIM("qflllgram");
2058 561797 : h = ZM_lll(RgM_rescale_to_int(x), D, flag);
2059 561768 : return gerepilecopy(av, h);
2060 : }
2061 :
2062 : GEN
2063 15673 : lllgram(GEN x) { return lllfp(x,LLLDFT,LLL_GRAM|LLL_IM); }
2064 : GEN
2065 556349 : lll(GEN x) { return lllfp(x,LLLDFT,LLL_IM); }
2066 :
2067 : GEN
2068 301 : qflll0(GEN x, long flag)
2069 : {
2070 301 : if (typ(x) != t_MAT) pari_err_TYPE("qflll",x);
2071 301 : switch(flag)
2072 : {
2073 49 : case 0: return lll(x);
2074 63 : case 1: RgM_check_ZM(x,"qflll"); return lllint(x);
2075 49 : case 2: RgM_check_ZM(x,"qflll"); return lllintpartial(x);
2076 7 : case 3: RgM_check_ZM(x,"qflll"); return lllall(x, LLL_INPLACE);
2077 49 : case 4: RgM_check_ZM(x,"qflll"); return lllkerim(x);
2078 42 : case 5: return lllkerimgen(x);
2079 42 : case 8: return lllgen(x);
2080 0 : default: pari_err_FLAG("qflll");
2081 : }
2082 : return NULL; /* LCOV_EXCL_LINE */
2083 : }
2084 :
2085 : GEN
2086 203 : qflllgram0(GEN x, long flag)
2087 : {
2088 203 : if (typ(x) != t_MAT) pari_err_TYPE("qflllgram",x);
2089 203 : switch(flag)
2090 : {
2091 63 : case 0: return lllgram(x);
2092 49 : case 1: RgM_check_ZM(x,"qflllgram"); return lllgramint(x);
2093 49 : case 4: RgM_check_ZM(x,"qflllgram"); return lllgramkerim(x);
2094 42 : case 5: return lllgramkerimgen(x);
2095 0 : case 8: return lllgramgen(x);
2096 0 : default: pari_err_FLAG("qflllgram");
2097 : }
2098 : return NULL; /* LCOV_EXCL_LINE */
2099 : }
2100 :
2101 : /********************************************************************/
2102 : /** **/
2103 : /** INTEGRAL KERNEL (LLL REDUCED) **/
2104 : /** **/
2105 : /********************************************************************/
2106 : static GEN
2107 70 : kerint0(GEN M)
2108 : {
2109 : /* return ZM_lll(M, LLLDFT, LLL_KER); */
2110 70 : GEN U, H = ZM_hnflll(M,&U,1);
2111 70 : long d = lg(M)-lg(H);
2112 70 : if (!d) return cgetg(1, t_MAT);
2113 70 : return ZM_lll(vecslice(U,1,d), LLLDFT, LLL_INPLACE);
2114 : }
2115 : GEN
2116 42 : kerint(GEN M)
2117 : {
2118 42 : pari_sp av = avma;
2119 42 : return gerepilecopy(av, kerint0(M));
2120 : }
2121 : /* OBSOLETE: use kerint */
2122 : GEN
2123 28 : matkerint0(GEN M, long flag)
2124 : {
2125 28 : pari_sp av = avma;
2126 28 : if (typ(M) != t_MAT) pari_err_TYPE("matkerint",M);
2127 28 : M = Q_primpart(M);
2128 28 : RgM_check_ZM(M, "kerint");
2129 28 : switch(flag)
2130 : {
2131 28 : case 0:
2132 28 : case 1: return gerepilecopy(av, kerint0(M));
2133 0 : default: pari_err_FLAG("matkerint");
2134 : }
2135 : return NULL; /* LCOV_EXCL_LINE */
2136 : }
|
