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 17184670 : pari_rint(double a)
22 : {
23 : #ifdef HAS_RINT
24 17184670 : 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 191958 : lll_trivial(GEN x, long flag)
45 : {
46 191958 : if (lg(x) == 1)
47 : { /* dim x = 0 */
48 15443 : 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 176515 : if (gequal0(gel(x,1)))
53 : {
54 91 : if (flag & LLL_KER) return matid(1);
55 91 : if (flag & (LLL_IM|LLL_INPLACE)) return cgetg(1,t_MAT);
56 28 : retmkvec2(matid(1), cgetg(1,t_MAT));
57 : }
58 176419 : if (flag & LLL_INPLACE) return gcopy(x);
59 73191 : if (flag & LLL_KER) return cgetg(1,t_MAT);
60 73191 : if (flag & LLL_IM) return matid(1);
61 28 : 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 3256173 : vectail_inplace(GEN x, long k)
67 : {
68 3256173 : if (!k) return x;
69 57552 : x[k] = ((ulong)x[0] & ~LGBITS) | evallg(lg(x) - k);
70 57553 : return x + k;
71 : }
72 :
73 : /* k = dim Kernel */
74 : static GEN
75 3329196 : lll_finish(GEN h, long k, long flag)
76 : {
77 : GEN g;
78 3329196 : if (!(flag & (LLL_IM|LLL_KER|LLL_ALL|LLL_INPLACE))) return h;
79 3256199 : if (flag & (LLL_IM|LLL_INPLACE)) return vectail_inplace(h, k);
80 96 : if (flag & LLL_KER) { setlg(h,k+1); return h; }
81 68 : 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 1926715 : mulshift(GEN y, GEN z, long e)
88 : {
89 1926715 : long ly = lgefint(y), lz;
90 : pari_sp av;
91 : GEN t;
92 1926715 : if (ly == 2) return gen_0;
93 204589 : lz = lgefint(z);
94 204589 : av = avma; (void)new_chunk(ly+lz+nbits2lg(e)); /* HACK */
95 204589 : t = mulii(z, y);
96 204589 : 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 6817078 : submulshift(GEN x, GEN y, GEN z, long e)
102 : {
103 6817078 : long lx = lgefint(x), ly, lz;
104 : pari_sp av;
105 : GEN t;
106 6817078 : if (!e) return submulii(x, y, z);
107 6786430 : if (lx == 2) { t = mulshift(y, z, e); togglesign(t); return t; }
108 5034992 : ly = lgefint(y);
109 5034992 : if (ly == 2) return icopy(x);
110 4556509 : lz = lgefint(z);
111 4556509 : av = avma; (void)new_chunk(lx+ly+lz+nbits2lg(e)); /* HACK */
112 4556509 : t = shifti(mulii(z, y), e);
113 4556509 : set_avma(av); return subii(x, t);
114 : }
115 : static void
116 211509099 : subzi(GEN *a, GEN b)
117 : {
118 211509099 : pari_sp av = avma;
119 211509099 : b = subii(*a, b);
120 211508951 : if (lgefint(b)<=lg(*a) && isonstack(*a)) { affii(b,*a); set_avma(av); }
121 10964882 : else *a = b;
122 211509060 : }
123 :
124 : static void
125 209919765 : addzi(GEN *a, GEN b)
126 : {
127 209919765 : pari_sp av = avma;
128 209919765 : b = addii(*a, b);
129 209919682 : if (lgefint(b)<=lg(*a) && isonstack(*a)) { affii(b,*a); set_avma(av); }
130 10743504 : else *a = b;
131 209919798 : }
132 :
133 : /* x - u*y * 2^e */
134 : INLINE GEN
135 20099554 : submuliu2n(GEN x, GEN y, ulong u, long e)
136 : {
137 : pari_sp av;
138 20099554 : long ly = lgefint(y);
139 20099554 : if (ly == 2) return x;
140 13881631 : av = avma;
141 13881631 : (void)new_chunk(3+ly+lgefint(x)+nbits2lg(e)); /* HACK */
142 13882681 : y = shifti(mului(u,y), e);
143 13881940 : set_avma(av); return subii(x, y);
144 : }
145 : /* *x -= u*y * 2^e */
146 : INLINE void
147 3451235 : submulzu2n(GEN *x, GEN y, ulong u, long e)
148 : {
149 : pari_sp av;
150 3451235 : long ly = lgefint(y);
151 3451235 : if (ly == 2) return;
152 2655728 : av = avma;
153 2655728 : (void)new_chunk(3+ly+lgefint(*x)+nbits2lg(e)); /* HACK */
154 2655728 : y = shifti(mului(u,y), e);
155 2655728 : set_avma(av); return subzi(x, y);
156 : }
157 :
158 : /* x + u*y * 2^e */
159 : INLINE GEN
160 20585373 : addmuliu2n(GEN x, GEN y, ulong u, long e)
161 : {
162 : pari_sp av;
163 20585373 : long ly = lgefint(y);
164 20585373 : if (ly == 2) return x;
165 14143170 : av = avma;
166 14143170 : (void)new_chunk(3+ly+lgefint(x)+nbits2lg(e)); /* HACK */
167 14144185 : y = shifti(mului(u,y), e);
168 14143451 : set_avma(av); return addii(x, y);
169 : }
170 :
171 : /* *x += u*y * 2^e */
172 : INLINE void
173 3258874 : addmulzu2n(GEN *x, GEN y, ulong u, long e)
174 : {
175 : pari_sp av;
176 3258874 : long ly = lgefint(y);
177 3258874 : if (ly == 2) return;
178 2530508 : av = avma;
179 2530508 : (void)new_chunk(3+ly+lgefint(*x)+nbits2lg(e)); /* HACK */
180 2530508 : y = shifti(mului(u,y), e);
181 2530508 : set_avma(av); return addzi(x, y);
182 : }
183 :
184 : /* n < 10; gerepileall supporting &NULL arguments. Maybe rename and export ? */
185 : INLINE void
186 9425 : gc_lll(pari_sp av, int n, ...)
187 : {
188 : int i, j;
189 : GEN *gptr[10];
190 : size_t s;
191 9425 : va_list a; va_start(a, n);
192 28275 : for (i=j=0; i<n; i++)
193 : {
194 18850 : GEN *x = va_arg(a,GEN*);
195 18850 : if (*x) { gptr[j++] = x; *x = (GEN)copy_bin(*x); }
196 : }
197 9425 : va_end(a); set_avma(av);
198 21602 : for (--j; j>=0; j--) *gptr[j] = bin_copy((GENbin*)*gptr[j]);
199 9425 : s = pari_mainstack->top - pari_mainstack->bot;
200 : /* size of saved objects ~ stacksize / 4 => overflow */
201 9425 : if (av - avma > (s >> 2))
202 : {
203 0 : size_t t = avma - pari_mainstack->bot;
204 0 : av = avma; new_chunk((s + t) / sizeof(long)); set_avma(av); /* double */
205 : }
206 9425 : }
207 :
208 : /********************************************************************/
209 : /** **/
210 : /** FPLLL (adapted from D. Stehle's code) **/
211 : /** **/
212 : /********************************************************************/
213 : /* Babai* and fplll* are a conversion to libpari API and data types
214 : of fplll-1.3 by Damien Stehle'.
215 :
216 : Copyright 2005, 2006 Damien Stehle'.
217 :
218 : This program is free software; you can redistribute it and/or modify it
219 : under the terms of the GNU General Public License as published by the
220 : Free Software Foundation; either version 2 of the License, or (at your
221 : option) any later version.
222 :
223 : This program implements ideas from the paper "Floating-point LLL Revisited",
224 : by Phong Nguyen and Damien Stehle', in the Proceedings of Eurocrypt'2005,
225 : Springer-Verlag; and was partly inspired by Shoup's NTL library:
226 : http://www.shoup.net/ntl/ */
227 :
228 : /* x t_REAL, |x| >= 1/2. Test whether |x| <= 3/2 */
229 : static int
230 1024349 : absrsmall2(GEN x)
231 : {
232 1024349 : long e = expo(x), l, i;
233 1024349 : if (e < 0) return 1;
234 325585 : if (e > 0 || (ulong)x[2] > (3UL << (BITS_IN_LONG-2))) return 0;
235 : /* line above assumes l > 2. OK since x != 0 */
236 167104 : l = lg(x); for (i = 3; i < l; i++) if (x[i]) return 0;
237 152436 : return 1;
238 : }
239 : /* x t_REAL; test whether |x| <= 1/2 */
240 : static int
241 2249244 : absrsmall(GEN x)
242 : {
243 : long e, l, i;
244 2249244 : if (!signe(x)) return 1;
245 2237422 : e = expo(x); if (e < -1) return 1;
246 1030721 : if (e > -1 || (ulong)x[2] > HIGHBIT) return 0;
247 7730 : l = lg(x); for (i = 3; i < l; i++) if (x[i]) return 0;
248 6372 : return 1;
249 : }
250 :
251 : static void
252 46427626 : rotate(GEN A, long k2, long k)
253 : {
254 : long i;
255 46427626 : GEN B = gel(A,k2);
256 143446722 : for (i = k2; i > k; i--) gel(A,i) = gel(A,i-1);
257 46427626 : gel(A,k) = B;
258 46427626 : }
259 :
260 : /************************* FAST version (double) ************************/
261 : #define dmael(x,i,j) ((x)[i][j])
262 : #define del(x,i) ((x)[i])
263 :
264 : static double *
265 42238881 : cget_dblvec(long d)
266 42238881 : { return (double*) stack_malloc_align(d*sizeof(double), sizeof(double)); }
267 :
268 : static double **
269 12986223 : cget_dblmat(long d) { return (double **) cgetg(d, t_VECSMALL); }
270 :
271 : static double
272 254262802 : itodbl_exp(GEN x, long *e)
273 : {
274 254262802 : pari_sp av = avma;
275 254262802 : GEN r = itor(x,DEFAULTPREC);
276 254231042 : *e = expo(r); setexpo(r,0);
277 254227712 : return gc_double(av, rtodbl(r));
278 : }
279 :
280 : static double
281 178606900 : dbldotproduct(double *x, double *y, long n)
282 : {
283 : long i;
284 178606900 : double sum = del(x,1) * del(y,1);
285 2877191432 : for (i=2; i<=n; i++) sum += del(x,i) * del(y,i);
286 178606900 : return sum;
287 : }
288 :
289 : static double
290 3608776 : dbldotsquare(double *x, long n)
291 : {
292 : long i;
293 3608776 : double sum = del(x,1) * del(x,1);
294 9433299 : for (i=2; i<=n; i++) sum += del(x,i) * del(x,i);
295 3608776 : return sum;
296 : }
297 :
298 : static long
299 33522812 : set_line(double *appv, GEN v, long n)
300 : {
301 33522812 : long i, maxexp = 0;
302 33522812 : pari_sp av = avma;
303 33522812 : GEN e = cgetg(n+1, t_VECSMALL);
304 287772459 : for (i = 1; i <= n; i++)
305 : {
306 254261535 : del(appv,i) = itodbl_exp(gel(v,i), e+i);
307 254248326 : if (e[i] > maxexp) maxexp = e[i];
308 : }
309 287850545 : for (i = 1; i <= n; i++) del(appv,i) = ldexp(del(appv,i), e[i]-maxexp);
310 33510924 : set_avma(av); return maxexp;
311 : }
312 :
313 : static void
314 51163542 : dblrotate(double **A, long k2, long k)
315 : {
316 : long i;
317 51163542 : double *B = del(A,k2);
318 158447476 : for (i = k2; i > k; i--) del(A,i) = del(A,i-1);
319 51163542 : del(A,k) = B;
320 51163542 : }
321 : /* update G[kappa][i] from appB */
322 : static void
323 30117569 : setG_fast(double **appB, long n, double **G, long kappa, long a, long b)
324 : { long i;
325 144851854 : for (i = a; i <= b; i++)
326 114734555 : dmael(G,kappa,i) = dbldotproduct(del(appB,kappa), del(appB,i), n);
327 30117299 : }
328 : /* update G[i][kappa] from appB */
329 : static void
330 24251753 : setG2_fast(double **appB, long n, double **G, long kappa, long a, long b)
331 : { long i;
332 88128071 : for (i = a; i <= b; i++)
333 63876403 : dmael(G,i,kappa) = dbldotproduct(del(appB,kappa), del(appB,i), n);
334 24251668 : }
335 : const long EX0 = -2; /* uninitialized; any value less than expo(0.51) = -1 */
336 :
337 : #ifdef LONG_IS_64BIT
338 : typedef long s64;
339 : #define addmuliu64_inplace addmuliu_inplace
340 : #define submuliu64_inplace submuliu_inplace
341 : #define submuliu642n submuliu2n
342 : #define addmuliu642n addmuliu2n
343 : #else
344 : typedef long long s64;
345 : typedef unsigned long long u64;
346 :
347 : INLINE GEN
348 33873301 : u64toi(u64 x)
349 : {
350 : GEN y;
351 : ulong h;
352 33873301 : if (!x) return gen_0;
353 33873301 : h = x>>32;
354 33873301 : if (!h) return utoipos(x);
355 3987690 : y = cgetipos(4);
356 3987690 : *int_LSW(y) = x&0xFFFFFFFF;
357 3987690 : *int_MSW(y) = x>>32;
358 3987690 : return y;
359 : }
360 :
361 : INLINE GEN
362 3195525 : u64toineg(u64 x)
363 : {
364 : GEN y;
365 : ulong h;
366 3195525 : if (!x) return gen_0;
367 3195525 : h = x>>32;
368 3195525 : if (!h) return utoineg(x);
369 3195525 : y = cgetineg(4);
370 3195525 : *int_LSW(y) = x&0xFFFFFFFF;
371 3195525 : *int_MSW(y) = x>>32;
372 3195525 : return y;
373 : }
374 : INLINE GEN
375 15253786 : addmuliu64_inplace(GEN x, GEN y, u64 u) { return addmulii(x, y, u64toi(u)); }
376 :
377 : INLINE GEN
378 15316291 : submuliu64_inplace(GEN x, GEN y, u64 u) { return submulii(x, y, u64toi(u)); }
379 :
380 : INLINE GEN
381 3195525 : addmuliu642n(GEN x, GEN y, u64 u, long e) { return submulshift(x, y, u64toineg(u), e); }
382 :
383 : INLINE GEN
384 3303224 : submuliu642n(GEN x, GEN y, u64 u, long e) { return submulshift(x, y, u64toi(u), e); }
385 :
386 : #endif
387 :
388 : /* Babai's Nearest Plane algorithm (iterative); see Babai() */
389 : static int
390 41081326 : Babai_fast(pari_sp av, long kappa, GEN *pB, GEN *pU, double **mu, double **r,
391 : double *s, double **appB, GEN expoB, double **G,
392 : long a, long zeros, long maxG, double eta)
393 : {
394 41081326 : GEN B = *pB, U = *pU;
395 41081326 : const long n = nbrows(B), d = U ? lg(U)-1: 0;
396 41081127 : long k, aa = (a > zeros)? a : zeros+1;
397 41081127 : long emaxmu = EX0, emax2mu = EX0;
398 : s64 xx;
399 41081127 : int did_something = 0;
400 : /* N.B: we set d = 0 (resp. n = 0) to avoid updating U (resp. B) */
401 :
402 24588596 : for (;;) {
403 65669723 : int go_on = 0;
404 65669723 : long i, j, emax3mu = emax2mu;
405 :
406 65669723 : if (gc_needed(av,2))
407 : {
408 1521 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[1], a=%ld", aa);
409 1521 : gc_lll(av,2,&B,&U);
410 : }
411 : /* Step2: compute the GSO for stage kappa */
412 65668888 : emax2mu = emaxmu; emaxmu = EX0;
413 269159926 : for (j=aa; j<kappa; j++)
414 : {
415 203491892 : double g = dmael(G,kappa,j);
416 1100248365 : for (k = zeros+1; k < j; k++) g -= dmael(mu,j,k) * dmael(r,kappa,k);
417 203491892 : dmael(r,kappa,j) = g;
418 203491892 : dmael(mu,kappa,j) = dmael(r,kappa,j) / dmael(r,j,j);
419 203491892 : emaxmu = maxss(emaxmu, expoB[kappa]-expoB[j]);
420 : }
421 : /* maxmu doesn't decrease fast enough */
422 65668034 : if (emax3mu != EX0 && emax3mu <= emax2mu + 5) {*pB = B; *pU = U; return 1;}
423 :
424 262072203 : for (j=kappa-1; j>zeros; j--)
425 : {
426 220996385 : double tmp = fabs(ldexp (dmael(mu,kappa,j), expoB[kappa]-expoB[j]));
427 220996385 : if (tmp>eta) { go_on = 1; break; }
428 : }
429 :
430 : /* Step3--5: compute the X_j's */
431 65664503 : if (go_on)
432 122952281 : for (j=kappa-1; j>zeros; j--)
433 : { /* The code below seemingly handles U = NULL, but in this case d = 0 */
434 98372377 : int e = expoB[j] - expoB[kappa];
435 98372377 : double tmp = ldexp(dmael(mu,kappa,j), -e), atmp = fabs(tmp);
436 : /* tmp = Inf is allowed */
437 98372377 : if (atmp <= .5) continue; /* size-reduced */
438 54314714 : if (gc_needed(av,2))
439 : {
440 4279 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[2], a=%ld, j=%ld", aa,j);
441 4279 : gc_lll(av,2,&B,&U);
442 : }
443 54316683 : did_something = 1;
444 : /* we consider separately the case |X| = 1 */
445 54316683 : if (atmp <= 1.5)
446 : {
447 36134798 : if (dmael(mu,kappa,j) > 0) { /* in this case, X = 1 */
448 94912123 : for (k=zeros+1; k<j; k++)
449 76354026 : dmael(mu,kappa,k) -= ldexp(dmael(mu,j,k), e);
450 342736652 : for (i=1; i<=n; i++)
451 324183390 : gmael(B,kappa,i) = subii(gmael(B,kappa,i), gmael(B,j,i));
452 113185711 : for (i=1; i<=d; i++)
453 94632283 : gmael(U,kappa,i) = subii(gmael(U,kappa,i), gmael(U,j,i));
454 : } else { /* otherwise X = -1 */
455 93272326 : for (k=zeros+1; k<j; k++)
456 75695625 : dmael(mu,kappa,k) += ldexp(dmael(mu,j,k), e);
457 337773779 : for (i=1; i<=n; i++)
458 320202268 : gmael(B,kappa,i) = addii(gmael(B,kappa,i), gmael(B,j,i));
459 108659978 : for (i=1; i<=d; i++)
460 91088315 : gmael(U,kappa,i) = addii(gmael(U,kappa,i), gmael(U,j,i));
461 : }
462 36125091 : continue;
463 : }
464 : /* we have |X| >= 2 */
465 18181885 : if (atmp < 9007199254740992.)
466 : {
467 15741167 : tmp = pari_rint(tmp);
468 38347965 : for (k=zeros+1; k<j; k++)
469 22606803 : dmael(mu,kappa,k) -= ldexp(tmp * dmael(mu,j,k), e);
470 15741162 : xx = (s64) tmp;
471 15741162 : if (xx > 0) /* = xx */
472 : {
473 80178977 : for (i=1; i<=n; i++)
474 72277238 : gmael(B,kappa,i) = submuliu64_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
475 44941686 : for (i=1; i<=d; i++)
476 37039785 : gmael(U,kappa,i) = submuliu64_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
477 : }
478 : else /* = -xx */
479 : {
480 79660851 : for (i=1; i<=n; i++)
481 71822364 : gmael(B,kappa,i) = addmuliu64_inplace(gmael(B,kappa,i), gmael(B,j,i), -xx);
482 44281641 : for (i=1; i<=d; i++)
483 36443075 : gmael(U,kappa,i) = addmuliu64_inplace(gmael(U,kappa,i), gmael(U,j,i), -xx);
484 : }
485 : }
486 : else
487 : {
488 : int E;
489 2440718 : xx = (s64) ldexp(frexp(dmael(mu,kappa,j), &E), 53);
490 2440718 : E -= e + 53;
491 2440718 : if (E <= 0)
492 : {
493 0 : xx = xx << -E;
494 0 : for (k=zeros+1; k<j; k++)
495 0 : dmael(mu,kappa,k) -= ldexp(((double)xx) * dmael(mu,j,k), e);
496 0 : if (xx > 0) /* = xx */
497 : {
498 0 : for (i=1; i<=n; i++)
499 0 : gmael(B,kappa,i) = submuliu64_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
500 0 : for (i=1; i<=d; i++)
501 0 : gmael(U,kappa,i) = submuliu64_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
502 : }
503 : else /* = -xx */
504 : {
505 0 : for (i=1; i<=n; i++)
506 0 : gmael(B,kappa,i) = addmuliu64_inplace(gmael(B,kappa,i), gmael(B,j,i), -xx);
507 0 : for (i=1; i<=d; i++)
508 0 : gmael(U,kappa,i) = addmuliu64_inplace(gmael(U,kappa,i), gmael(U,j,i), -xx);
509 : }
510 : } else
511 : {
512 16082302 : for (k=zeros+1; k<j; k++)
513 13641584 : dmael(mu,kappa,k) -= ldexp(((double)xx) * dmael(mu,j,k), E + e);
514 2440718 : if (xx > 0) /* = xx */
515 : {
516 22970848 : for (i=1; i<=n; i++)
517 21753881 : gmael(B,kappa,i) = submuliu642n(gmael(B,kappa,i), gmael(B,j,i), xx, E);
518 2112723 : for (i=1; i<=d; i++)
519 895756 : gmael(U,kappa,i) = submuliu642n(gmael(U,kappa,i), gmael(U,j,i), xx, E);
520 : }
521 : else /* = -xx */
522 : {
523 23414147 : for (i=1; i<=n; i++)
524 22190756 : gmael(B,kappa,i) = addmuliu642n(gmael(B,kappa,i), gmael(B,j,i), -xx, E);
525 2122694 : for (i=1; i<=d; i++)
526 899303 : gmael(U,kappa,i) = addmuliu642n(gmael(U,kappa,i), gmael(U,j,i), -xx, E);
527 : }
528 : }
529 : }
530 : }
531 65655705 : if (!go_on) break; /* Anything happened? */
532 24578340 : expoB[kappa] = set_line(del(appB,kappa), gel(B,kappa), n);
533 24588286 : setG_fast(appB, n, G, kappa, zeros+1, kappa-1);
534 24588596 : aa = zeros+1;
535 : }
536 41077365 : if (did_something) setG2_fast(appB, n, G, kappa, kappa, maxG);
537 :
538 41077702 : del(s,zeros+1) = dmael(G,kappa,kappa);
539 : /* the last s[kappa-1]=r[kappa][kappa] is computed only if kappa increases */
540 179805316 : for (k=zeros+1; k<=kappa-2; k++)
541 138727614 : del(s,k+1) = del(s,k) - dmael(mu,kappa,k)*dmael(r,kappa,k);
542 41077702 : *pB = B; *pU = U; return 0;
543 : }
544 :
545 : static void
546 18788861 : update_alpha(GEN alpha, long kappa, long kappa2, long kappamax)
547 : {
548 : long i;
549 59957608 : for (i = kappa; i < kappa2; i++)
550 41168747 : if (kappa <= alpha[i]) alpha[i] = kappa;
551 59957668 : for (i = kappa2; i > kappa; i--) alpha[i] = alpha[i-1];
552 49195826 : for (i = kappa2+1; i <= kappamax; i++)
553 30406965 : if (kappa < alpha[i]) alpha[i] = kappa;
554 18788861 : alpha[kappa] = kappa;
555 18788861 : }
556 : static void
557 1734202 : rotateG(GEN G, long kappa2, long kappa, long maxG, GEN Gtmp)
558 : {
559 : long i, j;
560 20879724 : for (i=1; i<=kappa2; i++) gel(Gtmp,i) = gmael(G,kappa2,i);
561 10446183 : for ( ; i<=maxG; i++) gel(Gtmp,i) = gmael(G,i,kappa2);
562 7141248 : for (i=kappa2; i>kappa; i--)
563 : {
564 37921853 : for (j=1; j<kappa; j++) gmael(G,i,j) = gmael(G,i-1,j);
565 5407046 : gmael(G,i,kappa) = gel(Gtmp,i-1);
566 25948538 : for (j=kappa+1; j<=i; j++) gmael(G,i,j) = gmael(G,i-1,j-1);
567 26940949 : for (j=kappa2+1; j<=maxG; j++) gmael(G,j,i) = gmael(G,j,i-1);
568 : }
569 13738476 : for (i=1; i<kappa; i++) gmael(G,kappa,i) = gel(Gtmp,i);
570 1734202 : gmael(G,kappa,kappa) = gel(Gtmp,kappa2);
571 10446183 : for (i=kappa2+1; i<=maxG; i++) gmael(G,i,kappa) = gel(Gtmp,i);
572 1734202 : }
573 : static void
574 17054580 : rotateG_fast(double **G, long kappa2, long kappa, long maxG, double *Gtmp)
575 : {
576 : long i, j;
577 103939611 : for (i=1; i<=kappa2; i++) del(Gtmp,i) = dmael(G,kappa2,i);
578 39734749 : for ( ; i<=maxG; i++) del(Gtmp,i) = dmael(G,i,kappa2);
579 52816149 : for (i=kappa2; i>kappa; i--)
580 : {
581 119435737 : for (j=1; j<kappa; j++) dmael(G,i,j) = dmael(G,i-1,j);
582 35761569 : dmael(G,i,kappa) = del(Gtmp,i-1);
583 128786832 : for (j=kappa+1; j<=i; j++) dmael(G,i,j) = dmael(G,i-1,j-1);
584 87068271 : for (j=kappa2+1; j<=maxG; j++) dmael(G,j,i) = dmael(G,j,i-1);
585 : }
586 51123721 : for (i=1; i<kappa; i++) dmael(G,kappa,i) = del(Gtmp,i);
587 17054580 : dmael(G,kappa,kappa) = del(Gtmp,kappa2);
588 39734788 : for (i=kappa2+1; i<=maxG; i++) dmael(G,i,kappa) = del(Gtmp,i);
589 17054580 : }
590 :
591 : /* LLL-reduces (B,U) in place [apply base change transforms to B and U].
592 : * Gram matrix, and GSO performed on matrices of 'double'.
593 : * If (keepfirst), never swap with first vector.
594 : * Return -1 on failure, else zeros = dim Kernel (>= 0) */
595 : static long
596 3246572 : fplll_fast(GEN *pB, GEN *pU, double delta, double eta, long keepfirst)
597 : {
598 : pari_sp av;
599 : long kappa, kappa2, d, n, i, j, zeros, kappamax, maxG;
600 : double **mu, **r, *s, tmp, *Gtmp, **G, **appB;
601 3246572 : GEN alpha, expoB, B = *pB, U;
602 3246572 : long cnt = 0;
603 :
604 3246572 : d = lg(B)-1;
605 3246572 : n = nbrows(B);
606 3246569 : U = *pU; /* NULL if inplace */
607 :
608 3246569 : G = cget_dblmat(d+1);
609 3246570 : appB = cget_dblmat(d+1);
610 3246568 : mu = cget_dblmat(d+1);
611 3246565 : r = cget_dblmat(d+1);
612 3246568 : s = cget_dblvec(d+1);
613 12183178 : for (j = 1; j <= d; j++)
614 : {
615 8936593 : del(mu,j) = cget_dblvec(d+1);
616 8936604 : del(r,j) = cget_dblvec(d+1);
617 8936610 : del(appB,j) = cget_dblvec(n+1);
618 8936607 : del(G,j) = cget_dblvec(d+1);
619 45451473 : for (i=1; i<=d; i++) dmael(G,j,i) = 0.;
620 : }
621 3246585 : expoB = cgetg(d+1, t_VECSMALL);
622 12183071 : for (i=1; i<=d; i++) expoB[i] = set_line(del(appB,i), gel(B,i), n);
623 3246502 : Gtmp = cget_dblvec(d+1);
624 3246567 : alpha = cgetg(d+1, t_VECSMALL);
625 3246559 : av = avma;
626 :
627 : /* Step2: Initializing the main loop */
628 3246559 : kappamax = 1;
629 3246559 : i = 1;
630 3246559 : maxG = d; /* later updated to kappamax */
631 :
632 : do {
633 3401101 : dmael(G,i,i) = dbldotsquare(del(appB,i),n);
634 3401103 : } while (dmael(G,i,i) <= 0 && (++i <=d));
635 3246561 : zeros = i-1; /* all vectors B[i] with i <= zeros are zero vectors */
636 3246561 : kappa = i;
637 3246561 : if (zeros < d) dmael(r,zeros+1,zeros+1) = dmael(G,zeros+1,zeros+1);
638 12028646 : for (i=zeros+1; i<=d; i++) alpha[i]=1;
639 44324437 : while (++kappa <= d)
640 : {
641 41081392 : if (kappa > kappamax)
642 : {
643 5529353 : if (DEBUGLEVEL>=4) err_printf("K%ld ",kappa);
644 5529353 : maxG = kappamax = kappa;
645 5529353 : setG_fast(appB, n, G, kappa, zeros+1, kappa);
646 : }
647 : /* Step3: Call to the Babai algorithm, mu,r,s updated in place */
648 41081401 : if (Babai_fast(av, kappa, &B,&U, mu,r,s, appB, expoB, G, alpha[kappa],
649 3531 : zeros, maxG, eta)) { *pB=B; *pU=U; return -1; }
650 :
651 41077642 : tmp = ldexp(r[kappa-1][kappa-1] * delta, 2*(expoB[kappa-1]-expoB[kappa]));
652 41077642 : if ((keepfirst && kappa == 2) || tmp <= del(s,kappa-1))
653 : { /* Step4: Success of Lovasz's condition */
654 24022980 : alpha[kappa] = kappa;
655 24022980 : tmp = dmael(mu,kappa,kappa-1) * dmael(r,kappa,kappa-1);
656 24022980 : dmael(r,kappa,kappa) = del(s,kappa-1)- tmp;
657 24022980 : continue;
658 : }
659 : /* Step5: Find the right insertion index kappa, kappa2 = initial kappa */
660 17054662 : if (DEBUGLEVEL>=4 && kappa==kappamax && del(s,kappa-1)!=0)
661 0 : if (++cnt > 20) { cnt = 0; err_printf("(%ld) ", 2*expoB[1] + dblexpo(del(s,1))); }
662 17054645 : kappa2 = kappa;
663 : do {
664 35761678 : kappa--;
665 35761678 : if (kappa<zeros+2 + (keepfirst ? 1: 0)) break;
666 25447016 : tmp = dmael(r,kappa-1,kappa-1) * delta;
667 25447016 : tmp = ldexp(tmp, 2*(expoB[kappa-1]-expoB[kappa2]));
668 25447016 : } while (del(s,kappa-1) <= tmp);
669 17054645 : update_alpha(alpha, kappa, kappa2, kappamax);
670 :
671 : /* Step6: Update the mu's and r's */
672 17054687 : dblrotate(mu,kappa2,kappa);
673 17054642 : dblrotate(r,kappa2,kappa);
674 17054605 : dmael(r,kappa,kappa) = del(s,kappa);
675 :
676 : /* Step7: Update B, appB, U, G */
677 17054605 : rotate(B,kappa2,kappa);
678 17054596 : dblrotate(appB,kappa2,kappa);
679 17054566 : if (U) rotate(U,kappa2,kappa);
680 17054564 : rotate(expoB,kappa2,kappa);
681 17054575 : rotateG_fast(G,kappa2,kappa, maxG, Gtmp);
682 :
683 : /* Step8: Prepare the next loop iteration */
684 17054896 : if (kappa == zeros+1 && dmael(G,kappa,kappa)<= 0)
685 : {
686 207677 : zeros++; kappa++;
687 207677 : dmael(G,kappa,kappa) = dbldotsquare(del(appB,kappa),n);
688 207677 : dmael(r,kappa,kappa) = dmael(G,kappa,kappa);
689 : }
690 : }
691 3243045 : *pB = B; *pU = U; return zeros;
692 : }
693 :
694 : /***************** HEURISTIC version (reduced precision) ****************/
695 : static GEN
696 329886 : realsqrdotproduct(GEN x)
697 : {
698 329886 : long i, l = lg(x);
699 329886 : GEN z = sqrr(gel(x,1));
700 4024724 : for (i=2; i<l; i++) z = addrr(z, sqrr(gel(x,i)));
701 329886 : return z;
702 : }
703 : /* x, y non-empty vector of t_REALs, same length */
704 : static GEN
705 3656008 : realdotproduct(GEN x, GEN y)
706 : {
707 : long i, l;
708 : GEN z;
709 3656008 : if (x == y) return realsqrdotproduct(x);
710 3326122 : l = lg(x); z = mulrr(gel(x,1),gel(y,1));
711 52976165 : for (i=2; i<l; i++) z = addrr(z, mulrr(gel(x,i), gel(y,i)));
712 3326122 : return z;
713 : }
714 : static void
715 337385 : setG_heuristic(GEN appB, GEN G, long kappa, long a, long b)
716 337385 : { pari_sp av = avma;
717 : long i;
718 2634186 : for (i = a; i <= b; i++)
719 2296801 : affrr(realdotproduct(gel(appB,kappa),gel(appB,i)), gmael(G,kappa,i));
720 337385 : set_avma(av);
721 337385 : }
722 : static void
723 318663 : setG2_heuristic(GEN appB, GEN G, long kappa, long a, long b)
724 318663 : { pari_sp av = avma;
725 : long i;
726 1677870 : for (i = a; i <= b; i++)
727 1359207 : affrr(realdotproduct(gel(appB,kappa),gel(appB,i)), gmael(G,i,kappa));
728 318663 : set_avma(av);
729 318663 : }
730 :
731 : /* approximate t_REAL x as m * 2^e, where |m| < 2^bit */
732 : static GEN
733 20439 : truncexpo(GEN x, long bit, long *e)
734 : {
735 20439 : *e = expo(x) + 1 - bit;
736 20439 : if (*e >= 0) return mantissa2nr(x, 0);
737 1787 : *e = 0; return roundr_safe(x);
738 : }
739 : /* Babai's Nearest Plane algorithm (iterative); see Babai() */
740 : static int
741 655382 : Babai_heuristic(pari_sp av, long kappa, GEN *pB, GEN *pU, GEN mu, GEN r, GEN s,
742 : GEN appB, GEN G, long a, long zeros, long maxG,
743 : GEN eta, long prec)
744 : {
745 655382 : GEN B = *pB, U = *pU;
746 655382 : const long n = nbrows(B), d = U ? lg(U)-1: 0, bit = prec2nbits(prec);
747 655382 : long k, aa = (a > zeros)? a : zeros+1;
748 655382 : int did_something = 0;
749 655382 : long emaxmu = EX0, emax2mu = EX0;
750 : /* N.B: we set d = 0 (resp. n = 0) to avoid updating U (resp. B) */
751 :
752 326162 : for (;;) {
753 981544 : int go_on = 0;
754 981544 : long i, j, emax3mu = emax2mu;
755 :
756 981544 : if (gc_needed(av,2))
757 : {
758 6 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[1], a=%ld", aa);
759 6 : gc_lll(av,2,&B,&U);
760 : }
761 : /* Step2: compute the GSO for stage kappa */
762 981544 : emax2mu = emaxmu; emaxmu = EX0;
763 5446427 : for (j=aa; j<kappa; j++)
764 : {
765 4464883 : pari_sp btop = avma;
766 4464883 : GEN g = gmael(G,kappa,j);
767 27941574 : for (k = zeros+1; k<j; k++)
768 23476691 : g = subrr(g, mulrr(gmael(mu,j,k), gmael(r,kappa,k)));
769 4464883 : affrr(g, gmael(r,kappa,j));
770 4464883 : affrr(divrr(gmael(r,kappa,j), gmael(r,j,j)), gmael(mu,kappa,j));
771 4464883 : emaxmu = maxss(emaxmu, expo(gmael(mu,kappa,j)));
772 4464883 : set_avma(btop);
773 : }
774 981544 : if (emax3mu != EX0 && emax3mu <= emax2mu + 5)
775 1177 : { *pB = B; *pU = U; return 1; }
776 :
777 5932720 : for (j=kappa-1; j>zeros; j--)
778 5278515 : if (abscmprr(gmael(mu,kappa,j), eta) > 0) { go_on = 1; break; }
779 :
780 : /* Step3--5: compute the X_j's */
781 980367 : if (go_on)
782 2575406 : for (j=kappa-1; j>zeros; j--)
783 : { /* The code below seemingly handles U = NULL, but in this case d = 0 */
784 : pari_sp btop;
785 2249244 : GEN tmp = gmael(mu,kappa,j);
786 2249244 : if (absrsmall(tmp)) continue; /* size-reduced */
787 :
788 1024349 : if (gc_needed(av,2))
789 : {
790 88 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[2], a=%ld, j=%ld", aa,j);
791 88 : gc_lll(av,2,&B,&U);
792 : }
793 1024349 : btop = avma; did_something = 1;
794 : /* we consider separately the case |X| = 1 */
795 1024349 : if (absrsmall2(tmp))
796 : {
797 851200 : if (signe(tmp) > 0) { /* in this case, X = 1 */
798 2381987 : for (k=zeros+1; k<j; k++)
799 1956009 : affrr(subrr(gmael(mu,kappa,k), gmael(mu,j,k)), gmael(mu,kappa,k));
800 425978 : set_avma(btop);
801 7408988 : for (i=1; i<=n; i++)
802 6983010 : gmael(B,kappa,i) = subii(gmael(B,kappa,i), gmael(B,j,i));
803 6194282 : for (i=1; i<=d; i++)
804 5768304 : gmael(U,kappa,i) = subii(gmael(U,kappa,i), gmael(U,j,i));
805 : } else { /* otherwise X = -1 */
806 2378459 : for (k=zeros+1; k<j; k++)
807 1953237 : affrr(addrr(gmael(mu,kappa,k), gmael(mu,j,k)), gmael(mu,kappa,k));
808 425222 : set_avma(btop);
809 7388388 : for (i=1; i<=n; i++)
810 6963166 : gmael(B,kappa,i) = addii(gmael(B,kappa,i), gmael(B,j,i));
811 6169988 : for (i=1; i<=d; i++)
812 5744766 : gmael(U,kappa,i) = addii(gmael(U,kappa,i),gmael(U,j,i));
813 : }
814 851200 : continue;
815 : }
816 : /* we have |X| >= 2 */
817 173149 : if (expo(tmp) < BITS_IN_LONG)
818 : {
819 152710 : ulong xx = roundr_safe(tmp)[2]; /* X fits in an ulong */
820 152710 : if (signe(tmp) > 0) /* = xx */
821 : {
822 251334 : for (k=zeros+1; k<j; k++)
823 175832 : affrr(subrr(gmael(mu,kappa,k), mulur(xx, gmael(mu,j,k))),
824 175832 : gmael(mu,kappa,k));
825 75502 : set_avma(btop);
826 912585 : for (i=1; i<=n; i++)
827 837083 : gmael(B,kappa,i) = submuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
828 722073 : for (i=1; i<=d; i++)
829 646571 : gmael(U,kappa,i) = submuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
830 : }
831 : else /* = -xx */
832 : {
833 259324 : for (k=zeros+1; k<j; k++)
834 182116 : affrr(addrr(gmael(mu,kappa,k), mulur(xx, gmael(mu,j,k))),
835 182116 : gmael(mu,kappa,k));
836 77208 : set_avma(btop);
837 941899 : for (i=1; i<=n; i++)
838 864691 : gmael(B,kappa,i) = addmuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
839 743053 : for (i=1; i<=d; i++)
840 665845 : gmael(U,kappa,i) = addmuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
841 : }
842 : }
843 : else
844 : {
845 : long e;
846 20439 : GEN X = truncexpo(tmp, bit, &e); /* tmp ~ X * 2^e */
847 20439 : btop = avma;
848 87782 : for (k=zeros+1; k<j; k++)
849 : {
850 67343 : GEN x = mulir(X, gmael(mu,j,k));
851 67343 : if (e) shiftr_inplace(x, e);
852 67343 : affrr(subrr(gmael(mu,kappa,k), x), gmael(mu,kappa,k));
853 : }
854 20439 : set_avma(btop);
855 263178 : for (i=1; i<=n; i++)
856 242739 : gmael(B,kappa,i) = submulshift(gmael(B,kappa,i), gmael(B,j,i), X, e);
857 96029 : for (i=1; i<=d; i++)
858 75590 : gmael(U,kappa,i) = submulshift(gmael(U,kappa,i), gmael(U,j,i), X, e);
859 : }
860 : }
861 980367 : if (!go_on) break; /* Anything happened? */
862 4326279 : for (i=1 ; i<=n; i++) affir(gmael(B,kappa,i), gmael(appB,kappa,i));
863 326162 : setG_heuristic(appB, G, kappa, zeros+1, kappa-1);
864 326162 : aa = zeros+1;
865 : }
866 654205 : if (did_something) setG2_heuristic(appB, G, kappa, kappa, maxG);
867 654205 : affrr(gmael(G,kappa,kappa), gel(s,zeros+1));
868 : /* the last s[kappa-1]=r[kappa][kappa] is computed only if kappa increases */
869 654205 : av = avma;
870 4521945 : for (k=zeros+1; k<=kappa-2; k++)
871 3867740 : affrr(subrr(gel(s,k), mulrr(gmael(mu,kappa,k), gmael(r,kappa,k))),
872 3867740 : gel(s,k+1));
873 654205 : *pB = B; *pU = U; return gc_bool(av, 0);
874 : }
875 :
876 : static GEN
877 17991 : ZC_to_RC(GEN x, long prec)
878 175874 : { pari_APPLY_type(t_COL,itor(gel(x,i),prec)) }
879 :
880 : static GEN
881 3531 : ZM_to_RM(GEN x, long prec)
882 21522 : { pari_APPLY_same(ZC_to_RC(gel(x,i),prec)) }
883 :
884 : /* LLL-reduces (B,U) in place [apply base change transforms to B and U].
885 : * Gram matrix made of t_REAL at precision prec2, performe GSO at prec.
886 : * If (keepfirst), never swap with first vector.
887 : * Return -1 on failure, else zeros = dim Kernel (>= 0) */
888 : static long
889 3531 : fplll_heuristic(GEN *pB, GEN *pU, double DELTA, double ETA, long keepfirst,
890 : long prec, long prec2)
891 : {
892 : pari_sp av, av2;
893 : long kappa, kappa2, d, i, j, zeros, kappamax, maxG;
894 3531 : GEN mu, r, s, tmp, Gtmp, alpha, G, appB, B = *pB, U;
895 3531 : GEN delta = dbltor(DELTA), eta = dbltor(ETA);
896 3531 : long cnt = 0;
897 :
898 3531 : d = lg(B)-1;
899 3531 : U = *pU; /* NULL if inplace */
900 :
901 3531 : G = cgetg(d+1, t_MAT);
902 3531 : mu = cgetg(d+1, t_MAT);
903 3531 : r = cgetg(d+1, t_MAT);
904 3531 : s = cgetg(d+1, t_VEC);
905 3531 : appB = ZM_to_RM(B, prec2);
906 21522 : for (j = 1; j <= d; j++)
907 : {
908 17991 : GEN M = cgetg(d+1, t_COL), R = cgetg(d+1, t_COL), S = cgetg(d+1, t_COL);
909 17991 : gel(mu,j)= M;
910 17991 : gel(r,j) = R;
911 17991 : gel(G,j) = S;
912 17991 : gel(s,j) = cgetr(prec);
913 167182 : for (i = 1; i <= d; i++)
914 : {
915 149191 : gel(R,i) = cgetr(prec);
916 149191 : gel(M,i) = cgetr(prec);
917 149191 : gel(S,i) = cgetr(prec2);
918 : }
919 : }
920 3531 : Gtmp = cgetg(d+1, t_VEC);
921 3531 : alpha = cgetg(d+1, t_VECSMALL);
922 3531 : av = avma;
923 :
924 : /* Step2: Initializing the main loop */
925 3531 : kappamax = 1;
926 3531 : i = 1;
927 3531 : maxG = d; /* later updated to kappamax */
928 :
929 : do {
930 3532 : affrr(RgV_dotsquare(gel(appB,i)), gmael(G,i,i));
931 3532 : } while (signe(gmael(G,i,i)) == 0 && (++i <=d));
932 3531 : zeros = i-1; /* all vectors B[i] with i <= zeros are zero vectors */
933 3531 : kappa = i;
934 3531 : if (zeros < d) affrr(gmael(G,zeros+1,zeros+1), gmael(r,zeros+1,zeros+1));
935 21521 : for (i=zeros+1; i<=d; i++) alpha[i]=1;
936 :
937 657736 : while (++kappa <= d)
938 : {
939 655382 : if (kappa > kappamax)
940 : {
941 11223 : if (DEBUGLEVEL>=4) err_printf("K%ld ",kappa);
942 11223 : maxG = kappamax = kappa;
943 11223 : setG_heuristic(appB, G, kappa, zeros+1, kappa);
944 : }
945 : /* Step3: Call to the Babai algorithm, mu,r,s updated in place */
946 655382 : if (Babai_heuristic(av, kappa, &B,&U, mu,r,s, appB, G, alpha[kappa], zeros,
947 1177 : maxG, eta, prec)) { *pB = B; *pU = U; return -1; }
948 654205 : av2 = avma;
949 1308302 : if ((keepfirst && kappa == 2) ||
950 654097 : cmprr(mulrr(gmael(r,kappa-1,kappa-1), delta), gel(s,kappa-1)) <= 0)
951 : { /* Step4: Success of Lovasz's condition */
952 440182 : alpha[kappa] = kappa;
953 440182 : tmp = mulrr(gmael(mu,kappa,kappa-1), gmael(r,kappa,kappa-1));
954 440182 : affrr(subrr(gel(s,kappa-1), tmp), gmael(r,kappa,kappa));
955 440182 : set_avma(av2); continue;
956 : }
957 : /* Step5: Find the right insertion index kappa, kappa2 = initial kappa */
958 214023 : if (DEBUGLEVEL>=4 && kappa==kappamax && signe(gel(s,kappa-1)))
959 0 : if (++cnt > 20) { cnt = 0; err_printf("(%ld) ", expo(gel(s,1))); }
960 214023 : kappa2 = kappa;
961 : do {
962 644194 : kappa--;
963 644194 : if (kappa < zeros+2 + (keepfirst ? 1: 0)) break;
964 614548 : tmp = mulrr(gmael(r,kappa-1,kappa-1), delta);
965 614548 : } while (cmprr(gel(s,kappa-1), tmp) <= 0 );
966 214023 : set_avma(av2);
967 214023 : update_alpha(alpha, kappa, kappa2, kappamax);
968 :
969 : /* Step6: Update the mu's and r's */
970 214023 : rotate(mu,kappa2,kappa);
971 214023 : rotate(r,kappa2,kappa);
972 214023 : affrr(gel(s,kappa), gmael(r,kappa,kappa));
973 :
974 : /* Step7: Update B, appB, U, G */
975 214023 : rotate(B,kappa2,kappa);
976 214023 : rotate(appB,kappa2,kappa);
977 214023 : if (U) rotate(U,kappa2,kappa);
978 214023 : rotateG(G,kappa2,kappa, maxG, Gtmp);
979 :
980 : /* Step8: Prepare the next loop iteration */
981 214023 : if (kappa == zeros+1 && !signe(gmael(G,kappa,kappa)))
982 : {
983 2 : zeros++; kappa++;
984 2 : affrr(RgV_dotsquare(gel(appB,kappa)), gmael(G,kappa,kappa));
985 2 : affrr(gmael(G,kappa,kappa), gmael(r,kappa,kappa));
986 : }
987 : }
988 2354 : *pB=B; *pU=U; return zeros;
989 : }
990 :
991 : /************************* PROVED version (t_INT) ***********************/
992 : /* dpe inspired by dpe.h by Patrick Pelissier, Paul Zimmermann
993 : * https://gforge.inria.fr/projects/dpe/
994 : */
995 :
996 : typedef struct
997 : {
998 : double d; /* significand */
999 : long e; /* exponent */
1000 : } dpe_t;
1001 :
1002 : #define Dmael(x,i,j) (&((x)[i][j]))
1003 : #define Del(x,i) (&((x)[i]))
1004 :
1005 : static void
1006 3040358 : dperotate(dpe_t **A, long k2, long k)
1007 : {
1008 : long i;
1009 3040358 : dpe_t *B = A[k2];
1010 12566062 : for (i = k2; i > k; i--) A[i] = A[i-1];
1011 3040358 : A[k] = B;
1012 3040358 : }
1013 :
1014 : static void
1015 833001992 : dpe_normalize0(dpe_t *x)
1016 : {
1017 : int e;
1018 833001992 : x->d = frexp(x->d, &e);
1019 833001992 : x->e += e;
1020 833001992 : }
1021 :
1022 : static void
1023 413237074 : dpe_normalize(dpe_t *x)
1024 : {
1025 413237074 : if (x->d == 0.0)
1026 2260279 : x->e = -LONG_MAX;
1027 : else
1028 410976795 : dpe_normalize0(x);
1029 413237142 : }
1030 :
1031 : static GEN
1032 93248 : dpetor(dpe_t *x)
1033 : {
1034 93248 : GEN r = dbltor(x->d);
1035 93248 : if (signe(r)==0) return r;
1036 93248 : setexpo(r, x->e-1);
1037 93248 : return r;
1038 : }
1039 :
1040 : static void
1041 93260932 : affdpe(dpe_t *y, dpe_t *x)
1042 : {
1043 93260932 : x->d = y->d;
1044 93260932 : x->e = y->e;
1045 93260932 : }
1046 :
1047 : static void
1048 65328654 : affidpe(GEN y, dpe_t *x)
1049 : {
1050 65328654 : pari_sp av = avma;
1051 65328654 : GEN r = itor(y, DEFAULTPREC);
1052 65327878 : x->e = expo(r)+1;
1053 65327878 : setexpo(r,-1);
1054 65327788 : x->d = rtodbl(r);
1055 65328110 : set_avma(av);
1056 65328178 : }
1057 :
1058 : static void
1059 6658256 : affdbldpe(double y, dpe_t *x)
1060 : {
1061 6658256 : x->d = (double)y;
1062 6658256 : x->e = 0;
1063 6658256 : dpe_normalize(x);
1064 6658258 : }
1065 :
1066 : static void
1067 395660729 : dpe_mulz(dpe_t *x, dpe_t *y, dpe_t *z)
1068 : {
1069 395660729 : z->d = x->d * y->d;
1070 395660729 : if (z->d == 0.0)
1071 22584710 : z->e = -LONG_MAX;
1072 : else
1073 : {
1074 373076019 : z->e = x->e + y->e;
1075 373076019 : dpe_normalize0(z);
1076 : }
1077 395661009 : }
1078 :
1079 : static void
1080 51310668 : dpe_divz(dpe_t *x, dpe_t *y, dpe_t *z)
1081 : {
1082 51310668 : z->d = x->d / y->d;
1083 51310668 : if (z->d == 0.0)
1084 2361067 : z->e = -LONG_MAX;
1085 : else
1086 : {
1087 48949601 : z->e = x->e - y->e;
1088 48949601 : dpe_normalize0(z);
1089 : }
1090 51310734 : }
1091 :
1092 : static void
1093 539807 : dpe_negz(dpe_t *y, dpe_t *x)
1094 : {
1095 539807 : x->d = - y->d;
1096 539807 : x->e = y->e;
1097 539807 : }
1098 :
1099 : static void
1100 30177587 : dpe_addz(dpe_t *y, dpe_t *z, dpe_t *x)
1101 : {
1102 30177587 : if (y->e > z->e + 53)
1103 1060602 : affdpe(y, x);
1104 29116985 : else if (z->e > y->e + 53)
1105 274619 : affdpe(z, x);
1106 : else
1107 : {
1108 28842366 : long d = y->e - z->e;
1109 :
1110 28842366 : if (d >= 0)
1111 : {
1112 21414458 : x->d = y->d + ldexp(z->d, -d);
1113 21414458 : x->e = y->e;
1114 : }
1115 : else
1116 : {
1117 7427908 : x->d = z->d + ldexp(y->d, d);
1118 7427908 : x->e = z->e;
1119 : }
1120 28842366 : dpe_normalize(x);
1121 : }
1122 30177594 : }
1123 : static void
1124 411416015 : dpe_subz(dpe_t *y, dpe_t *z, dpe_t *x)
1125 : {
1126 411416015 : if (y->e > z->e + 53)
1127 39094972 : affdpe(y, x);
1128 372321043 : else if (z->e > y->e + 53)
1129 539807 : dpe_negz(z, x);
1130 : else
1131 : {
1132 371781236 : long d = y->e - z->e;
1133 :
1134 371781236 : if (d >= 0)
1135 : {
1136 339380038 : x->d = y->d - ldexp(z->d, -d);
1137 339380038 : x->e = y->e;
1138 : }
1139 : else
1140 : {
1141 32401198 : x->d = ldexp(y->d, d) - z->d;
1142 32401198 : x->e = z->e;
1143 : }
1144 371781236 : dpe_normalize(x);
1145 : }
1146 411416364 : }
1147 :
1148 : static void
1149 5954957 : dpe_muluz(dpe_t *y, ulong t, dpe_t *x)
1150 : {
1151 5954957 : x->d = y->d * (double)t;
1152 5954957 : x->e = y->e;
1153 5954957 : dpe_normalize(x);
1154 5954956 : }
1155 :
1156 : static void
1157 2317688 : dpe_addmuluz(dpe_t *y, dpe_t *z, ulong t, dpe_t *x)
1158 : {
1159 : dpe_t tmp;
1160 2317688 : dpe_muluz(z, t, &tmp);
1161 2317688 : dpe_addz(y, &tmp, x);
1162 2317688 : }
1163 :
1164 : static void
1165 2446799 : dpe_submuluz(dpe_t *y, dpe_t *z, ulong t, dpe_t *x)
1166 : {
1167 : dpe_t tmp;
1168 2446799 : dpe_muluz(z, t, &tmp);
1169 2446799 : dpe_subz(y, &tmp, x);
1170 2446799 : }
1171 :
1172 : static void
1173 380968643 : dpe_submulz(dpe_t *y, dpe_t *z, dpe_t *t, dpe_t *x)
1174 : {
1175 : dpe_t tmp;
1176 380968643 : dpe_mulz(z, t, &tmp);
1177 380968664 : dpe_subz(y, &tmp, x);
1178 380968768 : }
1179 :
1180 : static int
1181 14692202 : dpe_cmp(dpe_t *x, dpe_t *y)
1182 : {
1183 14692202 : int sx = x->d < 0. ? -1: x->d > 0.;
1184 14692202 : int sy = y->d < 0. ? -1: y->d > 0.;
1185 14692202 : int d = sx - sy;
1186 :
1187 14692202 : if (d != 0)
1188 142858 : return d;
1189 14549344 : else if (x->e > y->e)
1190 2340772 : return (sx > 0) ? 1 : -1;
1191 12208572 : else if (y->e > x->e)
1192 4814384 : return (sx > 0) ? -1 : 1;
1193 : else
1194 7394188 : return (x->d < y->d) ? -1 : (x->d > y->d);
1195 : }
1196 :
1197 : static int
1198 63395780 : dpe_abscmp(dpe_t *x, dpe_t *y)
1199 : {
1200 63395780 : if (x->e > y->e)
1201 450541 : return 1;
1202 62945239 : else if (y->e > x->e)
1203 59961275 : return -1;
1204 : else
1205 2983964 : return (fabs(x->d) < fabs(y->d)) ? -1 : (fabs(x->d) > fabs(y->d));
1206 : }
1207 :
1208 : static int
1209 9684531 : dpe_abssmall(dpe_t *x)
1210 : {
1211 9684531 : return (x->e <= 0) || (x->e == 1 && fabs(x->d) <= .75);
1212 : }
1213 :
1214 : static int
1215 14692197 : dpe_cmpmul(dpe_t *x, dpe_t *y, dpe_t *z)
1216 : {
1217 : dpe_t t;
1218 14692197 : dpe_mulz(x,y,&t);
1219 14692200 : return dpe_cmp(&t, z);
1220 : }
1221 :
1222 : static dpe_t *
1223 21840130 : cget_dpevec(long d)
1224 21840130 : { return (dpe_t*) stack_malloc_align(d*sizeof(dpe_t), sizeof(dpe_t)); }
1225 :
1226 : static dpe_t **
1227 6658271 : cget_dpemat(long d) { return (dpe_t **) cgetg(d, t_VECSMALL); }
1228 :
1229 : static GEN
1230 19966 : dpeM_diagonal_shallow(dpe_t **m, long d)
1231 : {
1232 : long i;
1233 19966 : GEN y = cgetg(d+1,t_VEC);
1234 113214 : for (i=1; i<=d; i++) gel(y, i) = dpetor(Dmael(m,i,i));
1235 19966 : return y;
1236 : }
1237 :
1238 : static void
1239 9684509 : affii_or_copy_gc(pari_sp av, GEN x, GEN *y)
1240 : {
1241 9684509 : long l = lg(*y);
1242 9684509 : if (lgefint(x) <= l && isonstack(*y))
1243 : {
1244 9684081 : affii(x,*y);
1245 9684086 : set_avma(av);
1246 : }
1247 : else
1248 433 : *y = gerepileuptoint(av, x);
1249 9684520 : }
1250 :
1251 : /* *x -= u*y */
1252 : INLINE void
1253 26431384 : submulziu(GEN *x, GEN y, ulong u)
1254 : {
1255 : pari_sp av;
1256 26431384 : long ly = lgefint(y);
1257 26431384 : if (ly == 2) return;
1258 19558641 : av = avma;
1259 19558641 : (void)new_chunk(3+ly+lgefint(*x)); /* HACK */
1260 19558718 : y = mului(u,y);
1261 19558674 : set_avma(av); subzi(x, y);
1262 : }
1263 :
1264 : /* *x += u*y */
1265 : INLINE void
1266 24696100 : addmulziu(GEN *x, GEN y, ulong u)
1267 : {
1268 : pari_sp av;
1269 24696100 : long ly = lgefint(y);
1270 24696100 : if (ly == 2) return;
1271 18810341 : av = avma;
1272 18810341 : (void)new_chunk(3+ly+lgefint(*x)); /* HACK */
1273 18810386 : y = mului(u,y);
1274 18810351 : set_avma(av); addzi(x, y);
1275 : }
1276 :
1277 : /************************** PROVED version (dpe) *************************/
1278 :
1279 : /* Babai's Nearest Plane algorithm (iterative).
1280 : * Size-reduces b_kappa using mu_{i,j} and r_{i,j} for j<=i <kappa
1281 : * Update B[,kappa]; compute mu_{kappa,j}, r_{kappa,j} for j<=kappa and s[kappa]
1282 : * mu, r, s updated in place (affrr). Return 1 on failure, else 0. */
1283 : static int
1284 10286482 : Babai_dpe(pari_sp av, long kappa, GEN *pG, GEN *pB, GEN *pU, dpe_t **mu, dpe_t **r, dpe_t *s,
1285 : long a, long zeros, long maxG, dpe_t *eta)
1286 : {
1287 10286482 : GEN G = *pG, B = *pB, U = *pU, ztmp;
1288 10286482 : long k, d, n, aa = a > zeros? a: zeros+1;
1289 10286482 : long emaxmu = EX0, emax2mu = EX0;
1290 : /* N.B: we set d = 0 (resp. n = 0) to avoid updating U (resp. B) */
1291 10286482 : d = U? lg(U)-1: 0;
1292 10286482 : n = B? nbrows(B): 0;
1293 2170836 : for (;;) {
1294 12457354 : int go_on = 0;
1295 12457354 : long i, j, emax3mu = emax2mu;
1296 :
1297 12457354 : if (gc_needed(av,2))
1298 : {
1299 0 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[1], a=%ld", aa);
1300 0 : gc_lll(av,3,&G,&B,&U);
1301 : }
1302 : /* Step2: compute the GSO for stage kappa */
1303 12457323 : emax2mu = emaxmu; emaxmu = EX0;
1304 63768035 : for (j=aa; j<kappa; j++)
1305 : {
1306 : dpe_t g;
1307 51310702 : affidpe(gmael(G,kappa,j), &g);
1308 375974441 : for (k = zeros+1; k < j; k++)
1309 324663793 : dpe_submulz(&g, Dmael(mu,j,k), Dmael(r,kappa,k), &g);
1310 51310648 : affdpe(&g, Dmael(r,kappa,j));
1311 51310667 : dpe_divz(Dmael(r,kappa,j), Dmael(r,j,j), Dmael(mu,kappa,j));
1312 51310671 : emaxmu = maxss(emaxmu, Dmael(mu,kappa,j)->e);
1313 : }
1314 12457333 : if (emax3mu != EX0 && emax3mu <= emax2mu + 5) /* precision too low */
1315 0 : { *pG = G; *pB = B; *pU = U; return 1; }
1316 :
1317 73682302 : for (j=kappa-1; j>zeros; j--)
1318 63395786 : if (dpe_abscmp(Dmael(mu,kappa,j), eta) > 0) { go_on = 1; break; }
1319 :
1320 : /* Step3--5: compute the X_j's */
1321 12457333 : if (go_on)
1322 23470252 : for (j=kappa-1; j>zeros; j--)
1323 : {
1324 : pari_sp btop;
1325 21299420 : dpe_t *tmp = Dmael(mu,kappa,j);
1326 21299420 : if (tmp->e < 0) continue; /* (essentially) size-reduced */
1327 :
1328 9684530 : if (gc_needed(av,2))
1329 : {
1330 0 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[2], a=%ld, j=%ld", aa,j);
1331 0 : gc_lll(av,3,&G,&B,&U);
1332 : }
1333 : /* we consider separately the case |X| = 1 */
1334 9684530 : if (dpe_abssmall(tmp))
1335 : {
1336 8241037 : if (tmp->d > 0) { /* in this case, X = 1 */
1337 31513736 : for (k=zeros+1; k<j; k++)
1338 27386178 : dpe_subz(Dmael(mu,kappa,k), Dmael(mu,j,k), Dmael(mu,kappa,k));
1339 42750098 : for (i=1; i<=n; i++)
1340 38622540 : subzi(&gmael(B,kappa,i), gmael(B,j,i));
1341 68834844 : for (i=1; i<=d; i++)
1342 64707308 : subzi(&gmael(U,kappa,i), gmael(U,j,i));
1343 4127536 : btop = avma;
1344 4127536 : ztmp = subii(gmael(G,j,j), shifti(gmael(G,kappa,j), 1));
1345 4127551 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
1346 4127553 : affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
1347 35983082 : for (i=1; i<=j; i++)
1348 31855526 : subzi(&gmael(G,kappa,i), gmael(G,j,i));
1349 34572353 : for (i=j+1; i<kappa; i++)
1350 30444799 : subzi(&gmael(G,kappa,i), gmael(G,i,j));
1351 27792725 : for (i=kappa+1; i<=maxG; i++)
1352 23665172 : subzi(&gmael(G,i,kappa), gmael(G,i,j));
1353 : } else { /* otherwise X = -1 */
1354 31397351 : for (k=zeros+1; k<j; k++)
1355 27283867 : dpe_addz(Dmael(mu,kappa,k), Dmael(mu,j,k), Dmael(mu,kappa,k));
1356 42696668 : for (i=1; i<=n; i++)
1357 38583184 : addzi(&gmael(B,kappa,i),gmael(B,j,i));
1358 68485756 : for (i=1; i<=d; i++)
1359 64372291 : addzi(&gmael(U,kappa,i),gmael(U,j,i));
1360 4113465 : btop = avma;
1361 4113465 : ztmp = addii(gmael(G,j,j), shifti(gmael(G,kappa,j), 1));
1362 4113480 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
1363 4113479 : affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
1364 35787895 : for (i=1; i<=j; i++)
1365 31674413 : addzi(&gmael(G,kappa,i), gmael(G,j,i));
1366 34458589 : for (i=j+1; i<kappa; i++)
1367 30345106 : addzi(&gmael(G,kappa,i), gmael(G,i,j));
1368 27718097 : for (i=kappa+1; i<=maxG; i++)
1369 23604615 : addzi(&gmael(G,i,kappa), gmael(G,i,j));
1370 : }
1371 8241035 : continue;
1372 : }
1373 : /* we have |X| >= 2 */
1374 1443498 : if (tmp->e < BITS_IN_LONG-1)
1375 : {
1376 1268229 : if (tmp->d > 0)
1377 : {
1378 662302 : ulong xx = (ulong) pari_rint(ldexp(tmp->d, tmp->e)); /* X fits in an ulong */
1379 3109101 : for (k=zeros+1; k<j; k++)
1380 2446799 : dpe_submuluz(Dmael(mu,kappa,k), Dmael(mu,j,k), xx, Dmael(mu,kappa,k));
1381 4358711 : for (i=1; i<=n; i++)
1382 3696409 : submulziu(&gmael(B,kappa,i), gmael(B,j,i), xx);
1383 11195362 : for (i=1; i<=d; i++)
1384 10533060 : submulziu(&gmael(U,kappa,i), gmael(U,j,i), xx);
1385 662302 : btop = avma;
1386 662302 : ztmp = submuliu2n(mulii(gmael(G,j,j), sqru(xx)), gmael(G,kappa,j), xx, 1);
1387 662301 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
1388 662300 : affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
1389 4179472 : for (i=1; i<=j; i++)
1390 3517170 : submulziu(&gmael(G,kappa,i), gmael(G,j,i), xx);
1391 6602070 : for (i=j+1; i<kappa; i++)
1392 5939767 : submulziu(&gmael(G,kappa,i), gmael(G,i,j), xx);
1393 3407373 : for (i=kappa+1; i<=maxG; i++)
1394 2745069 : submulziu(&gmael(G,i,kappa), gmael(G,i,j), xx);
1395 : }
1396 : else
1397 : {
1398 605927 : ulong xx = (ulong) pari_rint(ldexp(-tmp->d, tmp->e)); /* X fits in an ulong */
1399 2923616 : for (k=zeros+1; k<j; k++)
1400 2317688 : dpe_addmuluz(Dmael(mu,kappa,k), Dmael(mu,j,k), xx, Dmael(mu,kappa,k));
1401 4297868 : for (i=1; i<=n; i++)
1402 3691940 : addmulziu(&gmael(B,kappa,i), gmael(B,j,i), xx);
1403 10260959 : for (i=1; i<=d; i++)
1404 9655030 : addmulziu(&gmael(U,kappa,i), gmael(U,j,i), xx);
1405 605929 : btop = avma;
1406 605929 : ztmp = addmuliu2n(mulii(gmael(G,j,j), sqru(xx)), gmael(G,kappa,j), xx, 1);
1407 605928 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
1408 605927 : affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
1409 3775519 : for (i=1; i<=j; i++)
1410 3169591 : addmulziu(&gmael(G,kappa,i), gmael(G,j,i), xx);
1411 6394449 : for (i=j+1; i<kappa; i++)
1412 5788520 : addmulziu(&gmael(G,kappa,i), gmael(G,i,j), xx);
1413 2997003 : for (i=kappa+1; i<=maxG; i++)
1414 2391074 : addmulziu(&gmael(G,i,kappa), gmael(G,i,j), xx);
1415 : }
1416 : }
1417 : else
1418 : {
1419 175269 : long e = tmp->e - BITS_IN_LONG + 1;
1420 175269 : if (tmp->d > 0)
1421 : {
1422 90612 : ulong xx = (ulong) pari_rint(ldexp(tmp->d, BITS_IN_LONG - 1));
1423 705052 : for (k=zeros+1; k<j; k++)
1424 : {
1425 : dpe_t x;
1426 614440 : dpe_muluz(Dmael(mu,j,k), xx, &x);
1427 614440 : x.e += e;
1428 614440 : dpe_subz(Dmael(mu,kappa,k), &x, Dmael(mu,kappa,k));
1429 : }
1430 1740992 : for (i=1; i<=n; i++)
1431 1650380 : submulzu2n(&gmael(B,kappa,i), gmael(B,j,i), xx, e);
1432 459717 : for (i=1; i<=d; i++)
1433 369105 : submulzu2n(&gmael(U,kappa,i), gmael(U,j,i), xx, e);
1434 90612 : btop = avma;
1435 90612 : ztmp = submuliu2n(mulshift(gmael(G,j,j), sqru(xx), 2*e),
1436 90612 : gmael(G,kappa,j), xx, e+1);
1437 90612 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
1438 90612 : affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
1439 795857 : for (i=1; i<=j; i++)
1440 705245 : submulzu2n(&gmael(G,kappa,i), gmael(G,j,i), xx, e);
1441 757140 : for ( ; i<kappa; i++)
1442 666528 : submulzu2n(&gmael(G,kappa,i), gmael(G,i,j), xx, e);
1443 150589 : for (i=kappa+1; i<=maxG; i++)
1444 59977 : submulzu2n(&gmael(G,i,kappa), gmael(G,i,j), xx, e);
1445 : } else
1446 : {
1447 84657 : ulong xx = (ulong) pari_rint(ldexp(-tmp->d, BITS_IN_LONG - 1));
1448 660698 : for (k=zeros+1; k<j; k++)
1449 : {
1450 : dpe_t x;
1451 576033 : dpe_muluz(Dmael(mu,j,k), xx, &x);
1452 576033 : x.e += e;
1453 576033 : dpe_addz(Dmael(mu,kappa,k), &x, Dmael(mu,kappa,k));
1454 : }
1455 1714144 : for (i=1; i<=n; i++)
1456 1629479 : addmulzu2n(&gmael(B,kappa,i), gmael(B,j,i), xx, e);
1457 373223 : for (i=1; i<=d; i++)
1458 288558 : addmulzu2n(&gmael(U,kappa,i), gmael(U,j,i), xx, e);
1459 84665 : btop = avma;
1460 84665 : ztmp = addmuliu2n(mulshift(gmael(G,j,j), sqru(xx), 2*e),
1461 84665 : gmael(G,kappa,j), xx, e+1);
1462 84665 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
1463 84665 : affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
1464 745524 : for (i=1; i<=j; i++)
1465 660859 : addmulzu2n(&gmael(G,kappa,i), gmael(G,j,i), xx, e);
1466 724460 : for ( ; i<kappa; i++)
1467 639795 : addmulzu2n(&gmael(G,kappa,i), gmael(G,i,j), xx, e);
1468 124848 : for (i=kappa+1; i<=maxG; i++)
1469 40183 : addmulzu2n(&gmael(G,i,kappa), gmael(G,i,j), xx, e);
1470 : }
1471 : }
1472 : }
1473 12457348 : if (!go_on) break; /* Anything happened? */
1474 2170836 : aa = zeros+1;
1475 : }
1476 :
1477 10286512 : affidpe(gmael(G,kappa,kappa), Del(s,zeros+1));
1478 : /* the last s[kappa-1]=r[kappa][kappa] is computed only if kappa increases */
1479 57825199 : for (k=zeros+1; k<=kappa-2; k++)
1480 47538701 : dpe_submulz(Del(s,k), Dmael(mu,kappa,k), Dmael(r,kappa,k), Del(s,k+1));
1481 10286498 : *pG = G; *pB = B; *pU = U; return 0;
1482 : }
1483 :
1484 : /* G integral Gram matrix, LLL-reduces (G,B,U) in place [apply base change
1485 : * transforms to B and U]. If (keepfirst), never swap with first vector.
1486 : * If G = NULL, we compute the Gram matrix incrementally.
1487 : * Return -1 on failure, else zeros = dim Kernel (>= 0) */
1488 : static long
1489 3329129 : fplll_dpe(GEN *pG, GEN *pB, GEN *pU, GEN *pr, double DELTA, double ETA,
1490 : long keepfirst)
1491 : {
1492 : pari_sp av;
1493 3329129 : GEN Gtmp, alpha, G = *pG, B = *pB, U = *pU;
1494 3329129 : long d, maxG, kappa, kappa2, i, j, zeros, kappamax, incgram = !G, cnt = 0;
1495 : dpe_t delta, eta, **mu, **r, *s;
1496 3329129 : affdbldpe(DELTA,&delta);
1497 3329131 : affdbldpe(ETA,&eta);
1498 :
1499 3329150 : if (incgram)
1500 : { /* incremental Gram matrix */
1501 3246592 : maxG = 2; d = lg(B)-1;
1502 3246592 : G = zeromatcopy(d, d);
1503 : }
1504 : else
1505 82558 : maxG = d = lg(G)-1;
1506 :
1507 3329149 : mu = cget_dpemat(d+1);
1508 3329133 : r = cget_dpemat(d+1);
1509 3329137 : s = cget_dpevec(d+1);
1510 12584702 : for (j = 1; j <= d; j++)
1511 : {
1512 9255553 : mu[j]= cget_dpevec(d+1);
1513 9255553 : r[j] = cget_dpevec(d+1);
1514 : }
1515 3329149 : Gtmp = cgetg(d+1, t_VEC);
1516 3329142 : alpha = cgetg(d+1, t_VECSMALL);
1517 3329141 : av = avma;
1518 :
1519 : /* Step2: Initializing the main loop */
1520 3329141 : kappamax = 1;
1521 3329141 : i = 1;
1522 : do {
1523 3696637 : if (incgram) gmael(G,i,i) = ZV_dotsquare(gel(B,i));
1524 3696552 : affidpe(gmael(G,i,i), Dmael(r,i,i));
1525 3696571 : } while (!signe(gmael(G,i,i)) && ++i <= d);
1526 3329075 : zeros = i-1; /* all basis vectors b_i with i <= zeros are zero vectors */
1527 3329075 : kappa = i;
1528 12216960 : for (i=zeros+1; i<=d; i++) alpha[i]=1;
1529 :
1530 13615619 : while (++kappa <= d)
1531 : {
1532 10286501 : if (kappa > kappamax)
1533 : {
1534 5558823 : if (DEBUGLEVEL>=4) err_printf("K%ld ",kappa);
1535 5558823 : kappamax = kappa;
1536 5558823 : if (incgram)
1537 : {
1538 21147530 : for (i=zeros+1; i<=kappa; i++)
1539 15819985 : gmael(G,kappa,i) = ZV_dotproduct(gel(B,kappa), gel(B,i));
1540 5327545 : maxG = kappamax;
1541 : }
1542 : }
1543 : /* Step3: Call to the Babai algorithm, mu,r,s updated in place */
1544 10286240 : if (Babai_dpe(av, kappa, &G,&B,&U, mu,r,s, alpha[kappa], zeros, maxG, &eta))
1545 0 : { *pG = incgram? NULL: G; *pB = B; *pU = U; return -1; }
1546 20475109 : if ((keepfirst && kappa == 2) ||
1547 10188592 : dpe_cmpmul(Dmael(r,kappa-1,kappa-1), &delta, Del(s,kappa-1)) <= 0)
1548 : { /* Step4: Success of Lovasz's condition */
1549 8766339 : alpha[kappa] = kappa;
1550 8766339 : dpe_submulz(Del(s,kappa-1), Dmael(mu,kappa,kappa-1), Dmael(r,kappa,kappa-1), Dmael(r,kappa,kappa));
1551 8766362 : continue;
1552 : }
1553 : /* Step5: Find the right insertion index kappa, kappa2 = initial kappa */
1554 1520178 : if (DEBUGLEVEL>=4 && kappa==kappamax && Del(s,kappa-1)->d)
1555 0 : if (++cnt > 20) { cnt = 0; err_printf("(%ld) ", Del(s,1)->e-1); }
1556 1520178 : kappa2 = kappa;
1557 : do {
1558 4762851 : kappa--;
1559 4762851 : if (kappa < zeros+2 + (keepfirst ? 1: 0)) break;
1560 4503597 : } while (dpe_cmpmul(Dmael(r,kappa-1,kappa-1), &delta, Del(s,kappa-1)) >= 0);
1561 1520178 : update_alpha(alpha, kappa, kappa2, kappamax);
1562 :
1563 : /* Step6: Update the mu's and r's */
1564 1520179 : dperotate(mu, kappa2, kappa);
1565 1520179 : dperotate(r, kappa2, kappa);
1566 1520179 : affdpe(Del(s,kappa), Dmael(r,kappa,kappa));
1567 :
1568 : /* Step7: Update G, B, U */
1569 1520179 : if (U) rotate(U, kappa2, kappa);
1570 1520179 : if (B) rotate(B, kappa2, kappa);
1571 1520179 : rotateG(G,kappa2,kappa, maxG, Gtmp);
1572 :
1573 : /* Step8: Prepare the next loop iteration */
1574 1520179 : if (kappa == zeros+1 && !signe(gmael(G,kappa,kappa)))
1575 : {
1576 35176 : zeros++; kappa++;
1577 35176 : affidpe(gmael(G,kappa,kappa), Dmael(r,kappa,kappa));
1578 : }
1579 : }
1580 3329118 : if (pr) *pr = dpeM_diagonal_shallow(r,d);
1581 3329118 : *pG = G; *pB = B; *pU = U; return zeros; /* success */
1582 : }
1583 :
1584 :
1585 : /************************** PROVED version (t_INT) *************************/
1586 :
1587 : /* Babai's Nearest Plane algorithm (iterative).
1588 : * Size-reduces b_kappa using mu_{i,j} and r_{i,j} for j<=i <kappa
1589 : * Update B[,kappa]; compute mu_{kappa,j}, r_{kappa,j} for j<=kappa and s[kappa]
1590 : * mu, r, s updated in place (affrr). Return 1 on failure, else 0. */
1591 : static int
1592 0 : Babai(pari_sp av, long kappa, GEN *pG, GEN *pB, GEN *pU, GEN mu, GEN r, GEN s,
1593 : long a, long zeros, long maxG, GEN eta, long prec)
1594 : {
1595 0 : GEN G = *pG, B = *pB, U = *pU, ztmp;
1596 0 : long k, aa = a > zeros? a: zeros+1;
1597 0 : const long n = B? nbrows(B): 0, d = U ? lg(U)-1: 0, bit = prec2nbits(prec);
1598 0 : long emaxmu = EX0, emax2mu = EX0;
1599 : /* N.B: we set d = 0 (resp. n = 0) to avoid updating U (resp. B) */
1600 :
1601 0 : for (;;) {
1602 0 : int go_on = 0;
1603 0 : long i, j, emax3mu = emax2mu;
1604 :
1605 0 : if (gc_needed(av,2))
1606 : {
1607 0 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[1], a=%ld", aa);
1608 0 : gc_lll(av,3,&G,&B,&U);
1609 : }
1610 : /* Step2: compute the GSO for stage kappa */
1611 0 : emax2mu = emaxmu; emaxmu = EX0;
1612 0 : for (j=aa; j<kappa; j++)
1613 : {
1614 0 : pari_sp btop = avma;
1615 0 : GEN g = gmael(G,kappa,j);
1616 0 : for (k = zeros+1; k < j; k++)
1617 0 : g = mpsub(g, mulrr(gmael(mu,j,k), gmael(r,kappa,k)));
1618 0 : mpaff(g, gmael(r,kappa,j));
1619 0 : affrr(divrr(gmael(r,kappa,j), gmael(r,j,j)), gmael(mu,kappa,j));
1620 0 : emaxmu = maxss(emaxmu, expo(gmael(mu,kappa,j)));
1621 0 : set_avma(btop);
1622 : }
1623 0 : if (emax3mu != EX0 && emax3mu <= emax2mu + 5) /* precision too low */
1624 0 : { *pG = G; *pB = B; *pU = U; return 1; }
1625 :
1626 0 : for (j=kappa-1; j>zeros; j--)
1627 0 : if (abscmprr(gmael(mu,kappa,j), eta) > 0) { go_on = 1; break; }
1628 :
1629 : /* Step3--5: compute the X_j's */
1630 0 : if (go_on)
1631 0 : for (j=kappa-1; j>zeros; j--)
1632 : {
1633 : pari_sp btop;
1634 0 : GEN tmp = gmael(mu,kappa,j);
1635 0 : if (absrsmall(tmp)) continue; /* size-reduced */
1636 :
1637 0 : if (gc_needed(av,2))
1638 : {
1639 0 : if(DEBUGMEM>1) pari_warn(warnmem,"Babai[2], a=%ld, j=%ld", aa,j);
1640 0 : gc_lll(av,3,&G,&B,&U);
1641 : }
1642 0 : btop = avma;
1643 : /* we consider separately the case |X| = 1 */
1644 0 : if (absrsmall2(tmp))
1645 : {
1646 0 : if (signe(tmp) > 0) { /* in this case, X = 1 */
1647 0 : for (k=zeros+1; k<j; k++)
1648 0 : affrr(subrr(gmael(mu,kappa,k), gmael(mu,j,k)), gmael(mu,kappa,k));
1649 0 : set_avma(btop);
1650 0 : for (i=1; i<=n; i++)
1651 0 : gmael(B,kappa,i) = subii(gmael(B,kappa,i), gmael(B,j,i));
1652 0 : for (i=1; i<=d; i++)
1653 0 : gmael(U,kappa,i) = subii(gmael(U,kappa,i), gmael(U,j,i));
1654 0 : btop = avma;
1655 0 : ztmp = subii(gmael(G,j,j), shifti(gmael(G,kappa,j), 1));
1656 0 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
1657 0 : gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
1658 0 : for (i=1; i<=j; i++)
1659 0 : gmael(G,kappa,i) = subii(gmael(G,kappa,i), gmael(G,j,i));
1660 0 : for (i=j+1; i<kappa; i++)
1661 0 : gmael(G,kappa,i) = subii(gmael(G,kappa,i), gmael(G,i,j));
1662 0 : for (i=kappa+1; i<=maxG; i++)
1663 0 : gmael(G,i,kappa) = subii(gmael(G,i,kappa), gmael(G,i,j));
1664 : } else { /* otherwise X = -1 */
1665 0 : for (k=zeros+1; k<j; k++)
1666 0 : affrr(addrr(gmael(mu,kappa,k), gmael(mu,j,k)), gmael(mu,kappa,k));
1667 0 : set_avma(btop);
1668 0 : for (i=1; i<=n; i++)
1669 0 : gmael(B,kappa,i) = addii(gmael(B,kappa,i),gmael(B,j,i));
1670 0 : for (i=1; i<=d; i++)
1671 0 : gmael(U,kappa,i) = addii(gmael(U,kappa,i),gmael(U,j,i));
1672 0 : btop = avma;
1673 0 : ztmp = addii(gmael(G,j,j), shifti(gmael(G,kappa,j), 1));
1674 0 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
1675 0 : gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
1676 0 : for (i=1; i<=j; i++)
1677 0 : gmael(G,kappa,i) = addii(gmael(G,kappa,i), gmael(G,j,i));
1678 0 : for (i=j+1; i<kappa; i++)
1679 0 : gmael(G,kappa,i) = addii(gmael(G,kappa,i), gmael(G,i,j));
1680 0 : for (i=kappa+1; i<=maxG; i++)
1681 0 : gmael(G,i,kappa) = addii(gmael(G,i,kappa), gmael(G,i,j));
1682 : }
1683 0 : continue;
1684 : }
1685 : /* we have |X| >= 2 */
1686 0 : if (expo(tmp) < BITS_IN_LONG)
1687 : {
1688 0 : ulong xx = roundr_safe(tmp)[2]; /* X fits in an ulong */
1689 0 : if (signe(tmp) > 0) /* = xx */
1690 : {
1691 0 : for (k=zeros+1; k<j; k++)
1692 0 : affrr(subrr(gmael(mu,kappa,k), mulur(xx, gmael(mu,j,k))),
1693 0 : gmael(mu,kappa,k));
1694 0 : set_avma(btop);
1695 0 : for (i=1; i<=n; i++)
1696 0 : gmael(B,kappa,i) = submuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
1697 0 : for (i=1; i<=d; i++)
1698 0 : gmael(U,kappa,i) = submuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
1699 0 : btop = avma;
1700 0 : ztmp = submuliu2n(mulii(gmael(G,j,j), sqru(xx)), gmael(G,kappa,j), xx, 1);
1701 0 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
1702 0 : gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
1703 0 : for (i=1; i<=j; i++)
1704 0 : gmael(G,kappa,i) = submuliu_inplace(gmael(G,kappa,i), gmael(G,j,i), xx);
1705 0 : for (i=j+1; i<kappa; i++)
1706 0 : gmael(G,kappa,i) = submuliu_inplace(gmael(G,kappa,i), gmael(G,i,j), xx);
1707 0 : for (i=kappa+1; i<=maxG; i++)
1708 0 : gmael(G,i,kappa) = submuliu_inplace(gmael(G,i,kappa), gmael(G,i,j), xx);
1709 : }
1710 : else /* = -xx */
1711 : {
1712 0 : for (k=zeros+1; k<j; k++)
1713 0 : affrr(addrr(gmael(mu,kappa,k), mulur(xx, gmael(mu,j,k))),
1714 0 : gmael(mu,kappa,k));
1715 0 : set_avma(btop);
1716 0 : for (i=1; i<=n; i++)
1717 0 : gmael(B,kappa,i) = addmuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
1718 0 : for (i=1; i<=d; i++)
1719 0 : gmael(U,kappa,i) = addmuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
1720 0 : btop = avma;
1721 0 : ztmp = addmuliu2n(mulii(gmael(G,j,j), sqru(xx)), gmael(G,kappa,j), xx, 1);
1722 0 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
1723 0 : gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
1724 0 : for (i=1; i<=j; i++)
1725 0 : gmael(G,kappa,i) = addmuliu_inplace(gmael(G,kappa,i), gmael(G,j,i), xx);
1726 0 : for (i=j+1; i<kappa; i++)
1727 0 : gmael(G,kappa,i) = addmuliu_inplace(gmael(G,kappa,i), gmael(G,i,j), xx);
1728 0 : for (i=kappa+1; i<=maxG; i++)
1729 0 : gmael(G,i,kappa) = addmuliu_inplace(gmael(G,i,kappa), gmael(G,i,j), xx);
1730 : }
1731 : }
1732 : else
1733 : {
1734 : long e;
1735 0 : GEN X = truncexpo(tmp, bit, &e); /* tmp ~ X * 2^e */
1736 0 : btop = avma;
1737 0 : for (k=zeros+1; k<j; k++)
1738 : {
1739 0 : GEN x = mulir(X, gmael(mu,j,k));
1740 0 : if (e) shiftr_inplace(x, e);
1741 0 : affrr(subrr(gmael(mu,kappa,k), x), gmael(mu,kappa,k));
1742 : }
1743 0 : set_avma(btop);
1744 0 : for (i=1; i<=n; i++)
1745 0 : gmael(B,kappa,i) = submulshift(gmael(B,kappa,i), gmael(B,j,i), X, e);
1746 0 : for (i=1; i<=d; i++)
1747 0 : gmael(U,kappa,i) = submulshift(gmael(U,kappa,i), gmael(U,j,i), X, e);
1748 0 : btop = avma;
1749 0 : ztmp = submulshift(mulshift(gmael(G,j,j), sqri(X), 2*e),
1750 0 : gmael(G,kappa,j), X, e+1);
1751 0 : ztmp = addii(gmael(G,kappa,kappa), ztmp);
1752 0 : gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
1753 0 : for (i=1; i<=j; i++)
1754 0 : gmael(G,kappa,i) = submulshift(gmael(G,kappa,i), gmael(G,j,i), X, e);
1755 0 : for ( ; i<kappa; i++)
1756 0 : gmael(G,kappa,i) = submulshift(gmael(G,kappa,i), gmael(G,i,j), X, e);
1757 0 : for (i=kappa+1; i<=maxG; i++)
1758 0 : gmael(G,i,kappa) = submulshift(gmael(G,i,kappa), gmael(G,i,j), X, e);
1759 : }
1760 : }
1761 0 : if (!go_on) break; /* Anything happened? */
1762 0 : aa = zeros+1;
1763 : }
1764 :
1765 0 : affir(gmael(G,kappa,kappa), gel(s,zeros+1));
1766 : /* the last s[kappa-1]=r[kappa][kappa] is computed only if kappa increases */
1767 0 : av = avma;
1768 0 : for (k=zeros+1; k<=kappa-2; k++)
1769 0 : affrr(subrr(gel(s,k), mulrr(gmael(mu,kappa,k), gmael(r,kappa,k))),
1770 0 : gel(s,k+1));
1771 0 : *pG = G; *pB = B; *pU = U; return gc_bool(av, 0);
1772 : }
1773 :
1774 : /* G integral Gram matrix, LLL-reduces (G,B,U) in place [apply base change
1775 : * transforms to B and U]. If (keepfirst), never swap with first vector.
1776 : * If G = NULL, we compute the Gram matrix incrementally.
1777 : * Return -1 on failure, else zeros = dim Kernel (>= 0) */
1778 : static long
1779 0 : fplll(GEN *pG, GEN *pB, GEN *pU, GEN *pr, double DELTA, double ETA,
1780 : long keepfirst, long prec)
1781 : {
1782 : pari_sp av, av2;
1783 0 : GEN mu, r, s, tmp, Gtmp, alpha, G = *pG, B = *pB, U = *pU;
1784 0 : GEN delta = dbltor(DELTA), eta = dbltor(ETA);
1785 0 : long d, maxG, kappa, kappa2, i, j, zeros, kappamax, incgram = !G, cnt = 0;
1786 :
1787 0 : if (incgram)
1788 : { /* incremental Gram matrix */
1789 0 : maxG = 2; d = lg(B)-1;
1790 0 : G = zeromatcopy(d, d);
1791 : }
1792 : else
1793 0 : maxG = d = lg(G)-1;
1794 :
1795 0 : mu = cgetg(d+1, t_MAT);
1796 0 : r = cgetg(d+1, t_MAT);
1797 0 : s = cgetg(d+1, t_VEC);
1798 0 : for (j = 1; j <= d; j++)
1799 : {
1800 0 : GEN M = cgetg(d+1, t_COL), R = cgetg(d+1, t_COL);
1801 0 : gel(mu,j)= M;
1802 0 : gel(r,j) = R;
1803 0 : gel(s,j) = cgetr(prec);
1804 0 : for (i = 1; i <= d; i++)
1805 : {
1806 0 : gel(R,i) = cgetr(prec);
1807 0 : gel(M,i) = cgetr(prec);
1808 : }
1809 : }
1810 0 : Gtmp = cgetg(d+1, t_VEC);
1811 0 : alpha = cgetg(d+1, t_VECSMALL);
1812 0 : av = avma;
1813 :
1814 : /* Step2: Initializing the main loop */
1815 0 : kappamax = 1;
1816 0 : i = 1;
1817 : do {
1818 0 : if (incgram) gmael(G,i,i) = ZV_dotsquare(gel(B,i));
1819 0 : affir(gmael(G,i,i), gmael(r,i,i));
1820 0 : } while (!signe(gmael(G,i,i)) && ++i <= d);
1821 0 : zeros = i-1; /* all basis vectors b_i with i <= zeros are zero vectors */
1822 0 : kappa = i;
1823 0 : for (i=zeros+1; i<=d; i++) alpha[i]=1;
1824 :
1825 0 : while (++kappa <= d)
1826 : {
1827 0 : if (kappa > kappamax)
1828 : {
1829 0 : if (DEBUGLEVEL>=4) err_printf("K%ld ",kappa);
1830 0 : kappamax = kappa;
1831 0 : if (incgram)
1832 : {
1833 0 : for (i=zeros+1; i<=kappa; i++)
1834 0 : gmael(G,kappa,i) = ZV_dotproduct(gel(B,kappa), gel(B,i));
1835 0 : maxG = kappamax;
1836 : }
1837 : }
1838 : /* Step3: Call to the Babai algorithm, mu,r,s updated in place */
1839 0 : if (Babai(av, kappa, &G,&B,&U, mu,r,s, alpha[kappa], zeros, maxG, eta, prec))
1840 0 : { *pG = incgram? NULL: G; *pB = B; *pU = U; return -1; }
1841 0 : av2 = avma;
1842 0 : if ((keepfirst && kappa == 2) ||
1843 0 : cmprr(mulrr(gmael(r,kappa-1,kappa-1), delta), gel(s,kappa-1)) <= 0)
1844 : { /* Step4: Success of Lovasz's condition */
1845 0 : alpha[kappa] = kappa;
1846 0 : tmp = mulrr(gmael(mu,kappa,kappa-1), gmael(r,kappa,kappa-1));
1847 0 : affrr(subrr(gel(s,kappa-1), tmp), gmael(r,kappa,kappa));
1848 0 : set_avma(av2); continue;
1849 : }
1850 : /* Step5: Find the right insertion index kappa, kappa2 = initial kappa */
1851 0 : if (DEBUGLEVEL>=4 && kappa==kappamax && signe(gel(s,kappa-1)))
1852 0 : if (++cnt > 20) { cnt = 0; err_printf("(%ld) ", expo(gel(s,1))); }
1853 0 : kappa2 = kappa;
1854 : do {
1855 0 : kappa--;
1856 0 : if (kappa < zeros+2 + (keepfirst ? 1: 0)) break;
1857 0 : tmp = mulrr(gmael(r,kappa-1,kappa-1), delta);
1858 0 : } while (cmprr(gel(s,kappa-1), tmp) <= 0);
1859 0 : set_avma(av2);
1860 0 : update_alpha(alpha, kappa, kappa2, kappamax);
1861 :
1862 : /* Step6: Update the mu's and r's */
1863 0 : rotate(mu, kappa2, kappa);
1864 0 : rotate(r, kappa2, kappa);
1865 0 : affrr(gel(s,kappa), gmael(r,kappa,kappa));
1866 :
1867 : /* Step7: Update G, B, U */
1868 0 : if (U) rotate(U, kappa2, kappa);
1869 0 : if (B) rotate(B, kappa2, kappa);
1870 0 : rotateG(G,kappa2,kappa, maxG, Gtmp);
1871 :
1872 : /* Step8: Prepare the next loop iteration */
1873 0 : if (kappa == zeros+1 && !signe(gmael(G,kappa,kappa)))
1874 : {
1875 0 : zeros++; kappa++;
1876 0 : affir(gmael(G,kappa,kappa), gmael(r,kappa,kappa));
1877 : }
1878 : }
1879 0 : if (pr) *pr = RgM_diagonal_shallow(r);
1880 0 : *pG = G; *pB = B; *pU = U; return zeros; /* success */
1881 : }
1882 :
1883 : /* Assume x a ZM, if pN != NULL, set it to Gram-Schmidt (squared) norms
1884 : * The following modes are supported:
1885 : * - flag & LLL_INPLACE: x a lattice basis, return x*U
1886 : * - flag & LLL_GRAM: x a Gram matrix / else x a lattice basis; return
1887 : * LLL base change matrix U [LLL_IM]
1888 : * kernel basis [LLL_KER, nonreduced]
1889 : * both [LLL_ALL] */
1890 : GEN
1891 3452156 : ZM_lll_norms(GEN x, double DELTA, long flag, GEN *pN)
1892 : {
1893 3452156 : pari_sp av = avma;
1894 3452156 : const double ETA = 0.51;
1895 3452156 : const long keepfirst = flag & LLL_KEEP_FIRST;
1896 3452156 : long p, zeros = -1, n = lg(x)-1;
1897 : GEN G, B, U;
1898 : pari_timer T;
1899 3452156 : if (n <= 1) return lll_trivial(x, flag);
1900 3344251 : if (nbrows(x)==0)
1901 : {
1902 15190 : if (flag & LLL_KER) return matid(n);
1903 15190 : if (flag & (LLL_INPLACE|LLL_IM)) return cgetg(1,t_MAT);
1904 0 : retmkvec2(matid(n), cgetg(1,t_MAT));
1905 : }
1906 3329104 : x = gcopy(x);
1907 3329133 : if (flag & LLL_GRAM)
1908 82548 : { G = x; B = NULL; U = matid(n); }
1909 : else
1910 3246585 : { G = NULL; B = x; U = (flag & LLL_INPLACE)? NULL: matid(n); }
1911 3329137 : if(DEBUGLEVEL>=4) timer_start(&T);
1912 3329137 : if (!(flag & LLL_GRAM))
1913 : {
1914 : long t;
1915 3246576 : if(DEBUGLEVEL>=4)
1916 0 : err_printf("Entering L^2 (double): LLL-parameters (%.3f,%.3f)\n",
1917 : DELTA,ETA);
1918 3246576 : zeros = fplll_fast(&B, &U, DELTA, ETA, keepfirst);
1919 3246576 : if (DEBUGLEVEL>=4) timer_printf(&T, zeros < 0? "LLL (failed)": "LLL");
1920 3250106 : for (p = DEFAULTPREC, t = 0; zeros < 0 && t < 1 ; p += EXTRAPREC64, t++)
1921 : {
1922 3531 : if (DEBUGLEVEL>=4)
1923 0 : err_printf("Entering L^2 (heuristic): LLL-parameters (%.3f,%.3f), prec = %d/%d\n", DELTA, ETA, p, p);
1924 3531 : zeros = fplll_heuristic(&B, &U, DELTA, ETA, keepfirst, p, p);
1925 3531 : gc_lll(av, 2, &B, &U);
1926 3531 : if (DEBUGLEVEL>=4) timer_printf(&T, zeros < 0? "LLL (failed)": "LLL");
1927 : }
1928 : }
1929 3329136 : if(DEBUGLEVEL>=4)
1930 0 : err_printf("Entering L^2 (dpe): LLL-parameters (%.3f,%.3f)\n", DELTA,ETA);
1931 3329136 : zeros = fplll_dpe(&G, &B, &U, pN, DELTA, ETA, keepfirst);
1932 3329115 : if (DEBUGLEVEL>=4) timer_printf(&T, zeros < 0? "LLL (failed)": "LLL");
1933 3329128 : if (zeros < 0)
1934 0 : for (p = DEFAULTPREC;; p += EXTRAPREC64)
1935 : {
1936 0 : if (DEBUGLEVEL>=4)
1937 0 : err_printf("Entering L^2: LLL-parameters (%.3f,%.3f), prec = %d\n",
1938 : DELTA,ETA, p);
1939 0 : zeros = fplll(&G, &B, &U, pN, DELTA, ETA, keepfirst, p);
1940 0 : if (DEBUGLEVEL>=4) timer_printf(&T, zeros < 0? "LLL (failed)": "LLL");
1941 0 : if (zeros >= 0) break;
1942 0 : gc_lll(av, 3, &G, &B, &U);
1943 : }
1944 3329128 : return lll_finish(U? U: B, zeros, flag);
1945 : }
1946 :
1947 : /********************************************************************/
1948 : /** **/
1949 : /** LLL OVER K[X] **/
1950 : /** **/
1951 : /********************************************************************/
1952 : static int
1953 504 : pslg(GEN x)
1954 : {
1955 : long tx;
1956 504 : if (gequal0(x)) return 2;
1957 448 : tx = typ(x); return is_scalar_t(tx)? 3: lg(x);
1958 : }
1959 :
1960 : static int
1961 196 : REDgen(long k, long l, GEN h, GEN L, GEN B)
1962 : {
1963 196 : GEN q, u = gcoeff(L,k,l);
1964 : long i;
1965 :
1966 196 : if (pslg(u) < pslg(B)) return 0;
1967 :
1968 140 : q = gneg(gdeuc(u,B));
1969 140 : gel(h,k) = gadd(gel(h,k), gmul(q,gel(h,l)));
1970 140 : for (i=1; i<l; i++) gcoeff(L,k,i) = gadd(gcoeff(L,k,i), gmul(q,gcoeff(L,l,i)));
1971 140 : gcoeff(L,k,l) = gadd(gcoeff(L,k,l), gmul(q,B)); return 1;
1972 : }
1973 :
1974 : static int
1975 196 : do_SWAPgen(GEN h, GEN L, GEN B, long k, GEN fl, int *flc)
1976 : {
1977 : GEN p1, la, la2, Bk;
1978 : long ps1, ps2, i, j, lx;
1979 :
1980 196 : if (!fl[k-1]) return 0;
1981 :
1982 140 : la = gcoeff(L,k,k-1); la2 = gsqr(la);
1983 140 : Bk = gel(B,k);
1984 140 : if (fl[k])
1985 : {
1986 56 : GEN q = gadd(la2, gmul(gel(B,k-1),gel(B,k+1)));
1987 56 : ps1 = pslg(gsqr(Bk));
1988 56 : ps2 = pslg(q);
1989 56 : if (ps1 <= ps2 && (ps1 < ps2 || !*flc)) return 0;
1990 28 : *flc = (ps1 != ps2);
1991 28 : gel(B,k) = gdiv(q, Bk);
1992 : }
1993 :
1994 112 : swap(gel(h,k-1), gel(h,k)); lx = lg(L);
1995 112 : for (j=1; j<k-1; j++) swap(gcoeff(L,k-1,j), gcoeff(L,k,j));
1996 112 : if (fl[k])
1997 : {
1998 28 : for (i=k+1; i<lx; i++)
1999 : {
2000 0 : GEN t = gcoeff(L,i,k);
2001 0 : p1 = gsub(gmul(gel(B,k+1),gcoeff(L,i,k-1)), gmul(la,t));
2002 0 : gcoeff(L,i,k) = gdiv(p1, Bk);
2003 0 : p1 = gadd(gmul(la,gcoeff(L,i,k-1)), gmul(gel(B,k-1),t));
2004 0 : gcoeff(L,i,k-1) = gdiv(p1, Bk);
2005 : }
2006 : }
2007 84 : else if (!gequal0(la))
2008 : {
2009 28 : p1 = gdiv(la2, Bk);
2010 28 : gel(B,k+1) = gel(B,k) = p1;
2011 28 : for (i=k+2; i<=lx; i++) gel(B,i) = gdiv(gmul(p1,gel(B,i)),Bk);
2012 28 : for (i=k+1; i<lx; i++)
2013 0 : gcoeff(L,i,k-1) = gdiv(gmul(la,gcoeff(L,i,k-1)), Bk);
2014 28 : for (j=k+1; j<lx-1; j++)
2015 0 : for (i=j+1; i<lx; i++)
2016 0 : gcoeff(L,i,j) = gdiv(gmul(p1,gcoeff(L,i,j)), Bk);
2017 : }
2018 : else
2019 : {
2020 56 : gcoeff(L,k,k-1) = gen_0;
2021 56 : for (i=k+1; i<lx; i++)
2022 : {
2023 0 : gcoeff(L,i,k) = gcoeff(L,i,k-1);
2024 0 : gcoeff(L,i,k-1) = gen_0;
2025 : }
2026 56 : gel(B,k) = gel(B,k-1); fl[k] = 1; fl[k-1] = 0;
2027 : }
2028 112 : return 1;
2029 : }
2030 :
2031 : static void
2032 168 : incrementalGSgen(GEN x, GEN L, GEN B, long k, GEN fl)
2033 : {
2034 168 : GEN u = NULL; /* gcc -Wall */
2035 : long i, j;
2036 420 : for (j = 1; j <= k; j++)
2037 252 : if (j==k || fl[j])
2038 : {
2039 252 : u = gcoeff(x,k,j);
2040 252 : if (!is_extscalar_t(typ(u))) pari_err_TYPE("incrementalGSgen",u);
2041 336 : for (i=1; i<j; i++)
2042 84 : if (fl[i])
2043 : {
2044 84 : u = gsub(gmul(gel(B,i+1),u), gmul(gcoeff(L,k,i),gcoeff(L,j,i)));
2045 84 : u = gdiv(u, gel(B,i));
2046 : }
2047 252 : gcoeff(L,k,j) = u;
2048 : }
2049 168 : if (gequal0(u)) gel(B,k+1) = gel(B,k);
2050 : else
2051 : {
2052 112 : gel(B,k+1) = gcoeff(L,k,k); gcoeff(L,k,k) = gen_1; fl[k] = 1;
2053 : }
2054 168 : }
2055 :
2056 : static GEN
2057 168 : lllgramallgen(GEN x, long flag)
2058 : {
2059 168 : long lx = lg(x), i, j, k, l, n;
2060 : pari_sp av;
2061 : GEN B, L, h, fl;
2062 : int flc;
2063 :
2064 168 : n = lx-1; if (n<=1) return lll_trivial(x,flag);
2065 84 : if (lgcols(x) != lx) pari_err_DIM("lllgramallgen");
2066 :
2067 84 : fl = cgetg(lx, t_VECSMALL);
2068 :
2069 84 : av = avma;
2070 84 : B = scalarcol_shallow(gen_1, lx);
2071 84 : L = cgetg(lx,t_MAT);
2072 252 : for (j=1; j<lx; j++) { gel(L,j) = zerocol(n); fl[j] = 0; }
2073 :
2074 84 : h = matid(n);
2075 252 : for (i=1; i<lx; i++)
2076 168 : incrementalGSgen(x, L, B, i, fl);
2077 84 : flc = 0;
2078 84 : for(k=2;;)
2079 : {
2080 196 : if (REDgen(k, k-1, h, L, gel(B,k))) flc = 1;
2081 196 : if (do_SWAPgen(h, L, B, k, fl, &flc)) { if (k > 2) k--; }
2082 : else
2083 : {
2084 84 : for (l=k-2; l>=1; l--)
2085 0 : if (REDgen(k, l, h, L, gel(B,l+1))) flc = 1;
2086 84 : if (++k > n) break;
2087 : }
2088 112 : if (gc_needed(av,1))
2089 : {
2090 0 : if(DEBUGMEM>1) pari_warn(warnmem,"lllgramallgen");
2091 0 : gerepileall(av,3,&B,&L,&h);
2092 : }
2093 : }
2094 140 : k=1; while (k<lx && !fl[k]) k++;
2095 84 : return lll_finish(h,k-1,flag);
2096 : }
2097 :
2098 : static int
2099 57688 : RgM_square(GEN x) { long l = lg(x); return l == 1 || l == lgcols(x); }
2100 : static GEN
2101 168 : lllallgen(GEN x, long flag)
2102 : {
2103 168 : pari_sp av = avma;
2104 168 : if (!(flag & LLL_GRAM)) x = gram_matrix(x);
2105 84 : else if (!RgM_square(x)) pari_err_DIM("qflllgram");
2106 168 : return gerepilecopy(av, lllgramallgen(x, flag));
2107 : }
2108 : GEN
2109 42 : lllgen(GEN x) { return lllallgen(x, LLL_IM); }
2110 : GEN
2111 42 : lllkerimgen(GEN x) { return lllallgen(x, LLL_ALL); }
2112 : GEN
2113 42 : lllgramgen(GEN x) { return lllallgen(x, LLL_IM|LLL_GRAM); }
2114 : GEN
2115 42 : lllgramkerimgen(GEN x) { return lllallgen(x, LLL_ALL|LLL_GRAM); }
2116 :
2117 : static GEN
2118 56139 : lllall(GEN x, long flag)
2119 56139 : { pari_sp av = avma; return gerepilecopy(av, ZM_lll(x, LLLDFT, flag)); }
2120 : GEN
2121 14683 : lllint(GEN x) { return lllall(x, LLL_IM); }
2122 : GEN
2123 35 : lllkerim(GEN x) { return lllall(x, LLL_ALL); }
2124 : GEN
2125 41379 : lllgramint(GEN x)
2126 41379 : { if (!RgM_square(x)) pari_err_DIM("qflllgram");
2127 41379 : return lllall(x, LLL_IM | LLL_GRAM); }
2128 : GEN
2129 35 : lllgramkerim(GEN x)
2130 35 : { if (!RgM_square(x)) pari_err_DIM("qflllgram");
2131 35 : return lllall(x, LLL_ALL | LLL_GRAM); }
2132 :
2133 : GEN
2134 1255019 : lllfp(GEN x, double D, long flag)
2135 : {
2136 1255019 : long n = lg(x)-1;
2137 1255019 : pari_sp av = avma;
2138 : GEN h;
2139 1255019 : if (n <= 1) return lll_trivial(x,flag);
2140 1171050 : if ((flag & LLL_GRAM) && !RgM_square(x)) pari_err_DIM("qflllgram");
2141 1171036 : h = ZM_lll(RgM_rescale_to_int(x), D, flag);
2142 1171008 : return gerepilecopy(av, h);
2143 : }
2144 :
2145 : GEN
2146 15967 : lllgram(GEN x) { return lllfp(x,LLLDFT,LLL_GRAM|LLL_IM); }
2147 : GEN
2148 1165545 : lll(GEN x) { return lllfp(x,LLLDFT,LLL_IM); }
2149 :
2150 : GEN
2151 301 : qflll0(GEN x, long flag)
2152 : {
2153 301 : if (typ(x) != t_MAT) pari_err_TYPE("qflll",x);
2154 301 : switch(flag)
2155 : {
2156 49 : case 0: return lll(x);
2157 63 : case 1: RgM_check_ZM(x,"qflll"); return lllint(x);
2158 49 : case 2: RgM_check_ZM(x,"qflll"); return lllintpartial(x);
2159 7 : case 3: RgM_check_ZM(x,"qflll"); return lllall(x, LLL_INPLACE);
2160 49 : case 4: RgM_check_ZM(x,"qflll"); return lllkerim(x);
2161 42 : case 5: return lllkerimgen(x);
2162 42 : case 8: return lllgen(x);
2163 0 : default: pari_err_FLAG("qflll");
2164 : }
2165 : return NULL; /* LCOV_EXCL_LINE */
2166 : }
2167 :
2168 : GEN
2169 245 : qflllgram0(GEN x, long flag)
2170 : {
2171 245 : if (typ(x) != t_MAT) pari_err_TYPE("qflllgram",x);
2172 245 : switch(flag)
2173 : {
2174 63 : case 0: return lllgram(x);
2175 49 : case 1: RgM_check_ZM(x,"qflllgram"); return lllgramint(x);
2176 49 : case 4: RgM_check_ZM(x,"qflllgram"); return lllgramkerim(x);
2177 42 : case 5: return lllgramkerimgen(x);
2178 42 : case 8: return lllgramgen(x);
2179 0 : default: pari_err_FLAG("qflllgram");
2180 : }
2181 : return NULL; /* LCOV_EXCL_LINE */
2182 : }
2183 :
2184 : /********************************************************************/
2185 : /** **/
2186 : /** INTEGRAL KERNEL (LLL REDUCED) **/
2187 : /** **/
2188 : /********************************************************************/
2189 : static GEN
2190 70 : kerint0(GEN M)
2191 : {
2192 : /* return ZM_lll(M, LLLDFT, LLL_KER); */
2193 70 : GEN U, H = ZM_hnflll(M,&U,1);
2194 70 : long d = lg(M)-lg(H);
2195 70 : if (!d) return cgetg(1, t_MAT);
2196 70 : return ZM_lll(vecslice(U,1,d), LLLDFT, LLL_INPLACE);
2197 : }
2198 : GEN
2199 42 : kerint(GEN M)
2200 : {
2201 42 : pari_sp av = avma;
2202 42 : return gerepilecopy(av, kerint0(M));
2203 : }
2204 : /* OBSOLETE: use kerint */
2205 : GEN
2206 28 : matkerint0(GEN M, long flag)
2207 : {
2208 28 : pari_sp av = avma;
2209 28 : if (typ(M) != t_MAT) pari_err_TYPE("matkerint",M);
2210 28 : M = Q_primpart(M);
2211 28 : RgM_check_ZM(M, "kerint");
2212 28 : switch(flag)
2213 : {
2214 28 : case 0:
2215 28 : case 1: return gerepilecopy(av, kerint0(M));
2216 0 : default: pari_err_FLAG("matkerint");
2217 : }
2218 : return NULL; /* LCOV_EXCL_LINE */
2219 : }
|
