Line data Source code
1 : /* Copyright (C) 2000 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 : /********************************************************************/
16 : /** **/
17 : /** LLL Algorithm and close friends **/
18 : /** **/
19 : /********************************************************************/
20 : #include "pari.h"
21 : #include "paripriv.h"
22 :
23 : #define DEBUGLEVEL DEBUGLEVEL_qf
24 :
25 : /********************************************************************/
26 : /** QR Factorization via Householder matrices **/
27 : /********************************************************************/
28 : static int
29 24775018 : no_prec_pb(GEN x)
30 : {
31 24699062 : return (typ(x) != t_REAL || realprec(x) > DEFAULTPREC
32 49474080 : || expo(x) < DEFAULTPREC>>1);
33 : }
34 : /* Find a Householder transformation which, applied to x[k..#x], zeroes
35 : * x[k+1..#x]; fill L = (mu_{i,j}). Return 0 if precision problem [obtained
36 : * a 0 vector], 1 otherwise */
37 : static int
38 24785920 : FindApplyQ(GEN x, GEN L, GEN B, long k, GEN Q, long prec)
39 : {
40 24785920 : long i, nx = lg(x)-1;
41 24785920 : GEN x2, x1, xd = x + (k-1);
42 :
43 24785920 : x1 = gel(xd,1);
44 24785920 : x2 = mpsqr(x1);
45 24785920 : if (k < nx)
46 : {
47 19543175 : long lv = nx - (k-1) + 1;
48 19543175 : GEN beta, Nx, v = cgetg(lv, t_VEC);
49 78001377 : for (i=2; i<lv; i++) {
50 58458202 : x2 = mpadd(x2, mpsqr(gel(xd,i)));
51 58458202 : gel(v,i) = gel(xd,i);
52 : }
53 19543175 : if (!signe(x2)) return 0;
54 19534951 : Nx = gsqrt(x2, prec); if (signe(x1) < 0) setsigne(Nx, -1);
55 19534951 : gel(v,1) = mpadd(x1, Nx);
56 :
57 19534951 : if (!signe(x1))
58 734109 : beta = gtofp(x2, prec); /* make sure typ(beta) != t_INT */
59 : else
60 18800842 : beta = mpadd(x2, mpmul(Nx,x1));
61 19534951 : gel(Q,k) = mkvec2(invr(beta), v); /* [t_REAL, vector of t_INT/t_REALs] */
62 :
63 19534951 : togglesign(Nx);
64 19534951 : gcoeff(L,k,k) = Nx; /* nonzero t_REAL */
65 : }
66 : else /* k = nx */
67 : {
68 5242745 : gcoeff(L,k,k) = x1; /* t_INT or t_REAL */
69 5242745 : if (!signe(x1)) return 0;
70 : }
71 24775018 : gel(B,k) = x2; /* t_INT or t_REAL */
72 71758475 : for (i=1; i<k; i++) gcoeff(L,k,i) = gel(x,i); /* t_INT or t_REAL */
73 24775018 : return no_prec_pb(x2);
74 : }
75 :
76 : /* apply Householder transformation Q = [beta,v] to r with t_INT/t_REAL
77 : * coefficients, in place: r -= ((0|v).r * beta) v */
78 : static void
79 47001261 : ApplyQ(GEN Q, GEN r)
80 : {
81 47001261 : GEN s, rd, beta = gel(Q,1), v = gel(Q,2);
82 47001261 : long i, l = lg(v), lr = lg(r);
83 :
84 47001261 : rd = r + (lr - l);
85 47001261 : s = mpmul(gel(v,1), gel(rd,1));
86 505949680 : for (i=2; i<l; i++) s = mpadd(s, mpmul(gel(v,i), gel(rd,i)));
87 47001261 : s = mpmul(beta, s);
88 552950941 : for (i=1; i<l; i++)
89 505949680 : if (signe(gel(v,i))) gel(rd,i) = mpsub(gel(rd,i), mpmul(s, gel(v,i)));
90 47001261 : }
91 : /* apply Q[1], ..., Q[j-1] to r */
92 : static GEN
93 17027474 : ApplyAllQ(GEN Q, GEN r, long j)
94 : {
95 17027474 : pari_sp av = avma;
96 : long i;
97 17027474 : r = leafcopy(r);
98 64028735 : for (i=1; i<j; i++) ApplyQ(gel(Q,i), r);
99 17027474 : return gc_GEN(av, r);
100 : }
101 :
102 : /* same, arbitrary coefficients [20% slower for t_REAL at DEFAULTPREC] */
103 : static void
104 22113 : RgC_ApplyQ(GEN Q, GEN r)
105 : {
106 22113 : GEN s, rd, beta = gel(Q,1), v = gel(Q,2);
107 22113 : long i, l = lg(v), lr = lg(r);
108 :
109 22113 : rd = r + (lr - l);
110 22113 : s = gmul(gel(v,1), gel(rd,1));
111 464373 : for (i=2; i<l; i++) s = gadd(s, gmul(gel(v,i), gel(rd,i)));
112 22113 : s = gmul(beta, s);
113 486486 : for (i=1; i<l; i++)
114 464373 : if (signe(gel(v,i))) gel(rd,i) = gsub(gel(rd,i), gmul(s, gel(v,i)));
115 22113 : }
116 : static GEN
117 567 : RgC_ApplyAllQ(GEN Q, GEN r, long j)
118 : {
119 567 : pari_sp av = avma;
120 : long i;
121 567 : r = leafcopy(r);
122 22680 : for (i=1; i<j; i++) RgC_ApplyQ(gel(Q,i), r);
123 567 : return gc_GEN(av, r);
124 : }
125 :
126 : int
127 21 : RgM_QR_init(GEN x, GEN *pB, GEN *pQ, GEN *pL, long prec)
128 : {
129 21 : x = RgM_gtomp(x, prec);
130 21 : return QR_init(x, pB, pQ, pL, prec);
131 : }
132 :
133 : static void
134 35 : check_householder(GEN Q)
135 : {
136 35 : long i, l = lg(Q);
137 35 : if (typ(Q) != t_VEC) pari_err_TYPE("mathouseholder", Q);
138 854 : for (i = 1; i < l; i++)
139 : {
140 826 : GEN q = gel(Q,i), v;
141 826 : if (typ(q) != t_VEC || lg(q) != 3) pari_err_TYPE("mathouseholder", Q);
142 826 : v = gel(q,2);
143 826 : if (typ(v) != t_VEC || lg(v)+i-2 != l) pari_err_TYPE("mathouseholder", Q);
144 : }
145 28 : }
146 :
147 : GEN
148 35 : mathouseholder(GEN Q, GEN x)
149 : {
150 35 : long l = lg(Q);
151 35 : check_householder(Q);
152 28 : switch(typ(x))
153 : {
154 14 : case t_MAT:
155 14 : if (lg(x) == 1) return cgetg(1, t_MAT);
156 14 : if (lgcols(x) != l+1) pari_err_TYPE("mathouseholder", x);
157 574 : pari_APPLY_same(RgC_ApplyAllQ(Q, gel(x,i), l));
158 7 : case t_COL:
159 7 : if (lg(x) == l+1) return RgC_ApplyAllQ(Q, x, l);
160 : }
161 7 : pari_err_TYPE("mathouseholder", x);
162 : return NULL; /* LCOV_EXCL_LINE */
163 : }
164 :
165 : GEN
166 35 : matqr(GEN x, long flag, long prec)
167 : {
168 35 : pari_sp av = avma;
169 : GEN B, Q, L;
170 35 : long n = lg(x)-1;
171 35 : if (typ(x) != t_MAT) pari_err_TYPE("matqr",x);
172 35 : if (!n)
173 : {
174 14 : if (!flag) retmkvec2(cgetg(1,t_MAT),cgetg(1,t_MAT));
175 7 : retmkvec2(cgetg(1,t_VEC),cgetg(1,t_MAT));
176 : }
177 21 : if (n != nbrows(x)) pari_err_DIM("matqr");
178 21 : if (!RgM_QR_init(x, &B,&Q,&L, prec)) pari_err_PREC("matqr");
179 21 : if (!flag) Q = shallowtrans(mathouseholder(Q, matid(n)));
180 21 : return gc_GEN(av, mkvec2(Q, shallowtrans(L)));
181 : }
182 :
183 : /* compute B = squared length of orthogonalized vectors x[k]^*,
184 : * Q = Householder transforms and L = mu_{i,j}. B[k] t_INT/t_REAL;
185 : * L[j,j] a t_REAL for j < #x */
186 : int
187 7758446 : QR_init(GEN x, GEN *pB, GEN *pQ, GEN *pL, long prec)
188 : {
189 7758446 : long j, k = lg(x)-1;
190 7758446 : GEN B = cgetg(k+1, t_VEC), Q = cgetg(k, t_VEC), L = zeromatcopy(k,k);
191 30298626 : for (j=1; j<=k; j++)
192 : {
193 24785920 : GEN r = gel(x,j);
194 24785920 : if (j > 1) r = ApplyAllQ(Q, r, j);
195 24785920 : if (!FindApplyQ(r, L, B, j, Q, prec)) return 0;
196 : }
197 5512706 : *pB = B; *pQ = Q; *pL = L; return 1;
198 : }
199 : /* x a square t_MAT with t_INT / t_REAL entries and maximal rank. Return
200 : * qfgaussred(x~*x) */
201 : GEN
202 301435 : gaussred_from_QR(GEN x, long prec)
203 : {
204 301435 : long j, k = lg(x)-1;
205 : GEN B, Q, L;
206 301435 : if (!QR_init(x, &B,&Q,&L, prec)) return NULL;
207 1075156 : for (j=1; j<k; j++)
208 : {
209 773721 : GEN m = gel(L,j), invNx = invr(gel(m,j));
210 : long i;
211 773721 : gel(m,j) = gel(B,j);
212 2995205 : for (i=j+1; i<=k; i++) gel(m,i) = mpmul(invNx, gel(m,i));
213 : }
214 301435 : gcoeff(L,j,j) = gel(B,j); /* t_INT or t_REAL */
215 301435 : return shallowtrans(L);
216 : }
217 : GEN
218 14280 : R_from_QR(GEN x, long prec)
219 : {
220 : GEN B, Q, L;
221 14280 : if (!QR_init(x, &B,&Q,&L, prec)) return NULL;
222 14266 : return shallowtrans(L);
223 : }
224 :
225 : /********************************************************************/
226 : /** QR Factorization via Gram-Schmidt **/
227 : /********************************************************************/
228 : /* return x + y, x possibly NULL (= initialized to gen_0) */
229 : static GEN
230 62737 : _add(GEN x, GEN y) { return x? gadd(x, y): y; }
231 :
232 : /* return Gram-Schmidt orthogonal basis (f) attached to (e), B is the
233 : * vector of the (f_i . f_i)*/
234 : GEN
235 56780 : RgM_gram_schmidt(GEN e, GEN *ptB)
236 : {
237 56780 : long i, j, lx = lg(e);
238 56780 : GEN f = RgM_shallowcopy(e), B = cgetg(lx, t_VEC), iB = cgetg(lx, t_VEC);
239 :
240 120645 : for (i = 1; i < lx; i++)
241 : {
242 63865 : pari_sp av = avma;
243 63865 : GEN c = NULL;
244 126602 : for (j = 1; j < i; j++)
245 : {
246 62737 : GEN mu = gmul(RgV_dotproduct(gel(e,i),gel(f,j)), gel(iB,j));
247 62737 : c = _add(c, gmul(mu, gel(f,j)));
248 : }
249 63865 : gel(f,i) = c? gc_upto(av, gsub(gel(e,i), c)): gel(e,i);
250 63865 : gel(B,i) = RgV_dotsquare(gel(f,i));
251 63865 : gel(iB,i) = ginv(gel(B,i));
252 : }
253 56780 : *ptB = B; return f;
254 : }
255 :
256 : /* B a Z-basis (which the caller should LLL-reduce for efficiency), t a vector.
257 : * Apply Babai's nearest plane algorithm to (B,t) */
258 : GEN
259 56780 : RgM_Babai(GEN B, GEN t)
260 : {
261 56780 : GEN C, N, G = RgM_gram_schmidt(B, &N), b = t;
262 56780 : long j, n = lg(B)-1;
263 :
264 56780 : C = cgetg(n+1,t_COL);
265 120645 : for (j = n; j > 0; j--)
266 : {
267 63865 : GEN c = gdiv( RgV_dotproduct(b, gel(G,j)), gel(N,j) );
268 : long e;
269 63865 : c = grndtoi(c,&e);
270 63865 : if (e >= 0) return NULL;
271 63865 : if (signe(c)) b = RgC_sub(b, RgC_Rg_mul(gel(B,j), c));
272 63865 : gel(C,j) = c;
273 : }
274 56780 : return C;
275 : }
276 :
277 : /********************************************************************/
278 : /** **/
279 : /** LLL ALGORITHM **/
280 : /** **/
281 : /********************************************************************/
282 : /* Def: a matrix M is said to be -partially reduced- if | m1 +- m2 | >= |m1|
283 : * for any two columns m1 != m2, in M.
284 : *
285 : * Input: an integer matrix mat whose columns are linearly independent. Find
286 : * another matrix T such that mat * T is partially reduced.
287 : *
288 : * Output: mat * T if flag = 0; T if flag != 0,
289 : *
290 : * This routine is designed to quickly reduce lattices in which one row
291 : * is huge compared to the other rows. For example, when searching for a
292 : * polynomial of degree 3 with root a mod N, the four input vectors might
293 : * be the coefficients of
294 : * X^3 - (a^3 mod N), X^2 - (a^2 mod N), X - (a mod N), N.
295 : * All four constant coefficients are O(p) and the rest are O(1). By the
296 : * pigeon-hole principle, the coefficients of the smallest vector in the
297 : * lattice are O(p^(1/4)), hence significant reduction of vector lengths
298 : * can be anticipated.
299 : *
300 : * An improved algorithm would look only at the leading digits of dot*. It
301 : * would use single-precision calculations as much as possible.
302 : *
303 : * Original code: Peter Montgomery (1994) */
304 : static GEN
305 35 : lllintpartialall(GEN m, long flag)
306 : {
307 35 : const long ncol = lg(m)-1;
308 35 : const pari_sp av = avma;
309 : GEN tm1, tm2, mid;
310 :
311 35 : if (ncol <= 1) return flag? matid(ncol): gcopy(m);
312 :
313 14 : tm1 = flag? matid(ncol): NULL;
314 : {
315 14 : const pari_sp av2 = avma;
316 14 : GEN dot11 = ZV_dotsquare(gel(m,1));
317 14 : GEN dot22 = ZV_dotsquare(gel(m,2));
318 14 : GEN dot12 = ZV_dotproduct(gel(m,1), gel(m,2));
319 14 : GEN tm = matid(2); /* For first two columns only */
320 :
321 14 : int progress = 0;
322 14 : long npass2 = 0;
323 :
324 : /* Row reduce the first two columns of m. Our best result so far is
325 : * (first two columns of m)*tm.
326 : *
327 : * Initially tm = 2 x 2 identity matrix.
328 : * Inner products of the reduced matrix are in dot11, dot12, dot22. */
329 49 : while (npass2 < 2 || progress)
330 : {
331 42 : GEN dot12new, q = diviiround(dot12, dot22);
332 :
333 35 : npass2++; progress = signe(q);
334 35 : if (progress)
335 : {/* Conceptually replace (v1, v2) by (v1 - q*v2, v2), where v1 and v2
336 : * represent the reduced basis for the first two columns of the matrix.
337 : * We do this by updating tm and the inner products. */
338 21 : togglesign(q);
339 21 : dot12new = addii(dot12, mulii(q, dot22));
340 21 : dot11 = addii(dot11, mulii(q, addii(dot12, dot12new)));
341 21 : dot12 = dot12new;
342 21 : ZC_lincomb1_inplace(gel(tm,1), gel(tm,2), q);
343 : }
344 :
345 : /* Interchange the output vectors v1 and v2. */
346 35 : swap(dot11,dot22);
347 35 : swap(gel(tm,1), gel(tm,2));
348 :
349 : /* Occasionally (including final pass) do garbage collection. */
350 35 : if ((npass2 & 0xff) == 0 || !progress)
351 14 : (void)gc_all(av2, 4, &dot11,&dot12,&dot22,&tm);
352 : } /* while npass2 < 2 || progress */
353 :
354 : {
355 : long i;
356 7 : GEN det12 = subii(mulii(dot11, dot22), sqri(dot12));
357 :
358 7 : mid = cgetg(ncol+1, t_MAT);
359 21 : for (i = 1; i <= 2; i++)
360 : {
361 14 : GEN tmi = gel(tm,i);
362 14 : if (tm1)
363 : {
364 14 : GEN tm1i = gel(tm1,i);
365 14 : gel(tm1i,1) = gel(tmi,1);
366 14 : gel(tm1i,2) = gel(tmi,2);
367 : }
368 14 : gel(mid,i) = ZC_lincomb(gel(tmi,1),gel(tmi,2), gel(m,1),gel(m,2));
369 : }
370 42 : for (i = 3; i <= ncol; i++)
371 : {
372 35 : GEN c = gel(m,i);
373 35 : GEN dot1i = ZV_dotproduct(gel(mid,1), c);
374 35 : GEN dot2i = ZV_dotproduct(gel(mid,2), c);
375 : /* ( dot11 dot12 ) (q1) ( dot1i )
376 : * ( dot12 dot22 ) (q2) = ( dot2i )
377 : *
378 : * Round -q1 and -q2 to nearest integer. Then compute
379 : * c - q1*mid[1] - q2*mid[2].
380 : * This will be approximately orthogonal to the first two vectors in
381 : * the new basis. */
382 35 : GEN q1neg = subii(mulii(dot12, dot2i), mulii(dot22, dot1i));
383 35 : GEN q2neg = subii(mulii(dot12, dot1i), mulii(dot11, dot2i));
384 :
385 35 : q1neg = diviiround(q1neg, det12);
386 35 : q2neg = diviiround(q2neg, det12);
387 35 : if (tm1)
388 : {
389 35 : gcoeff(tm1,1,i) = addii(mulii(q1neg, gcoeff(tm,1,1)),
390 35 : mulii(q2neg, gcoeff(tm,1,2)));
391 35 : gcoeff(tm1,2,i) = addii(mulii(q1neg, gcoeff(tm,2,1)),
392 35 : mulii(q2neg, gcoeff(tm,2,2)));
393 : }
394 35 : gel(mid,i) = ZC_add(c, ZC_lincomb(q1neg,q2neg, gel(mid,1),gel(mid,2)));
395 : } /* for i */
396 : } /* local block */
397 : }
398 7 : if (DEBUGLEVEL>6)
399 : {
400 0 : if (tm1) err_printf("tm1 = %Ps",tm1);
401 0 : err_printf("mid = %Ps",mid); /* = m * tm1 */
402 : }
403 7 : (void)gc_all(av, tm1? 2: 1, &mid, &tm1);
404 : {
405 : /* For each pair of column vectors v and w in mid * tm2,
406 : * try to replace (v, w) by (v, v - q*w) for some q.
407 : * We compute all inner products and check them repeatedly. */
408 7 : const pari_sp av3 = avma;
409 7 : long i, j, npass = 0, e = LONG_MAX;
410 7 : GEN dot = cgetg(ncol+1, t_MAT); /* scalar products */
411 :
412 7 : tm2 = matid(ncol);
413 56 : for (i=1; i <= ncol; i++)
414 : {
415 49 : gel(dot,i) = cgetg(ncol+1,t_COL);
416 245 : for (j=1; j <= i; j++)
417 196 : gcoeff(dot,j,i) = gcoeff(dot,i,j) = ZV_dotproduct(gel(mid,i),gel(mid,j));
418 : }
419 : for(;;)
420 35 : {
421 42 : long reductions = 0, olde = e;
422 336 : for (i=1; i <= ncol; i++)
423 : {
424 : long ijdif;
425 2058 : for (ijdif=1; ijdif < ncol; ijdif++)
426 : {
427 : long d, k1, k2;
428 : GEN codi, q;
429 :
430 1764 : j = i + ijdif; if (j > ncol) j -= ncol;
431 : /* let k1, resp. k2, index of larger, resp. smaller, column */
432 1764 : if (cmpii(gcoeff(dot,i,i), gcoeff(dot,j,j)) > 0) { k1 = i; k2 = j; }
433 1022 : else { k1 = j; k2 = i; }
434 1764 : codi = gcoeff(dot,k2,k2);
435 1764 : q = signe(codi)? diviiround(gcoeff(dot,k1,k2), codi): gen_0;
436 1764 : if (!signe(q)) continue;
437 :
438 : /* Try to subtract a multiple of column k2 from column k1. */
439 700 : reductions++; togglesign_safe(&q);
440 700 : ZC_lincomb1_inplace(gel(tm2,k1), gel(tm2,k2), q);
441 700 : ZC_lincomb1_inplace(gel(dot,k1), gel(dot,k2), q);
442 700 : gcoeff(dot,k1,k1) = addii(gcoeff(dot,k1,k1),
443 700 : mulii(q, gcoeff(dot,k2,k1)));
444 5600 : for (d = 1; d <= ncol; d++) gcoeff(dot,k1,d) = gcoeff(dot,d,k1);
445 : } /* for ijdif */
446 294 : if (gc_needed(av3,2))
447 : {
448 0 : if(DEBUGMEM>1) pari_warn(warnmem,"lllintpartialall");
449 0 : (void)gc_all(av3, 2, &dot,&tm2);
450 : }
451 : } /* for i */
452 42 : if (!reductions) break;
453 35 : e = 0;
454 280 : for (i = 1; i <= ncol; i++) e += expi( gcoeff(dot,i,i) );
455 35 : if (e == olde) break;
456 35 : if (DEBUGLEVEL>6)
457 : {
458 0 : npass++;
459 0 : err_printf("npass = %ld, red. last time = %ld, log_2(det) ~ %ld\n\n",
460 : npass, reductions, e);
461 : }
462 : } /* for(;;)*/
463 :
464 : /* Sort columns so smallest comes first in m * tm1 * tm2.
465 : * Use insertion sort. */
466 49 : for (i = 1; i < ncol; i++)
467 : {
468 42 : long j, s = i;
469 :
470 189 : for (j = i+1; j <= ncol; j++)
471 147 : if (cmpii(gcoeff(dot,s,s),gcoeff(dot,j,j)) > 0) s = j;
472 42 : if (i != s)
473 : { /* Exchange with proper column; only the diagonal of dot is updated */
474 28 : swap(gel(tm2,i), gel(tm2,s));
475 28 : swap(gcoeff(dot,i,i), gcoeff(dot,s,s));
476 : }
477 : }
478 : } /* local block */
479 7 : return gc_upto(av, ZM_mul(tm1? tm1: mid, tm2));
480 : }
481 :
482 : GEN
483 35 : lllintpartial(GEN mat) { return lllintpartialall(mat,1); }
484 :
485 : GEN
486 0 : lllintpartial_inplace(GEN mat) { return lllintpartialall(mat,0); }
487 :
488 : /********************************************************************/
489 : /** **/
490 : /** COPPERSMITH ALGORITHM **/
491 : /** Finding small roots of univariate equations. **/
492 : /** **/
493 : /********************************************************************/
494 :
495 : static int
496 882 : check(double b, double x, double rho, long d, long dim, long delta, long t)
497 : {
498 882 : double cond = delta * (d * (delta+1) - 2*b*dim + rho * (delta-1 + 2*t))
499 882 : + x*dim*(dim - 1);
500 882 : if (DEBUGLEVEL >= 4)
501 0 : err_printf("delta = %d, t = %d (%.1lf)\n", delta, t, cond);
502 882 : return (cond <= 0);
503 : }
504 :
505 : static void
506 21 : choose_params(GEN P, GEN N, GEN X, GEN B, long *pdelta, long *pt)
507 : {
508 21 : long d = degpol(P), dim;
509 21 : GEN P0 = leading_coeff(P);
510 21 : double logN = dbllog2(N), x, b, rho;
511 21 : x = dbllog2(X) / logN;
512 21 : b = B? dbllog2(B) / logN: 1.;
513 21 : if (x * d >= b * b) pari_err_OVERFLOW("zncoppersmith [bound too large]");
514 : /* TODO : remove P0 completely */
515 14 : rho = is_pm1(P0)? 0: dbllog2(P0) / logN;
516 :
517 : /* Enumerate (delta,t) by increasing lattice dimension */
518 14 : for(dim = d + 1;; dim++)
519 161 : {
520 : long delta, t; /* dim = d*delta + t in the loop */
521 1043 : for (delta = 0, t = dim; t >= 0; delta++, t -= d)
522 882 : if (check(b,x,rho,d,dim,delta,t)) { *pdelta = delta; *pt = t; return; }
523 : }
524 : }
525 :
526 : static int
527 14021 : sol_OK(GEN x, GEN N, GEN B)
528 14021 : { return B? (cmpii(gcdii(x,N),B) >= 0): dvdii(x,N); }
529 : /* deg(P) > 0, x >= 0. Find all j such that gcd(P(j), N) >= B, |j| <= x */
530 : static GEN
531 7 : do_exhaustive(GEN P, GEN N, long x, GEN B)
532 : {
533 7 : GEN Pe, Po, sol = vecsmalltrunc_init(2*x + 2);
534 : pari_sp av;
535 : long j;
536 7 : RgX_even_odd(P, &Pe,&Po); av = avma;
537 7 : if (sol_OK(gel(P,2), N,B)) vecsmalltrunc_append(sol, 0);
538 7007 : for (j = 1; j <= x; j++, set_avma(av))
539 : {
540 7000 : GEN j2 = sqru(j), E = FpX_eval(Pe,j2,N), O = FpX_eval(Po,j2,N);
541 7000 : if (sol_OK(addmuliu(E,O,j), N,B)) vecsmalltrunc_append(sol, j);
542 7000 : if (sol_OK(submuliu(E,O,j), N,B)) vecsmalltrunc_append(sol,-j);
543 : }
544 7 : vecsmall_sort(sol); return zv_to_ZV(sol);
545 : }
546 :
547 : /* General Coppersmith, look for a root x0 <= p, p >= B, p | N, |x0| <= X.
548 : * B = N coded as NULL */
549 : GEN
550 35 : zncoppersmith(GEN P, GEN N, GEN X, GEN B)
551 : {
552 : GEN Q, R, N0, M, sh, short_pol, *Xpowers, sol, nsp, cP, Z;
553 35 : long delta, i, j, row, d, l, t, dim, bnd = 10;
554 35 : const ulong X_SMALL = 1000;
555 35 : pari_sp av = avma;
556 :
557 35 : if (typ(P) != t_POL || !RgX_is_ZX(P)) pari_err_TYPE("zncoppersmith",P);
558 28 : if (typ(N) != t_INT) pari_err_TYPE("zncoppersmith",N);
559 28 : if (typ(X) != t_INT) {
560 7 : X = gfloor(X);
561 7 : if (typ(X) != t_INT) pari_err_TYPE("zncoppersmith",X);
562 : }
563 28 : if (signe(X) < 0) pari_err_DOMAIN("zncoppersmith", "X", "<", gen_0, X);
564 28 : P = FpX_red(P, N); d = degpol(P);
565 28 : if (d == 0) retgc_const(av, cgetg(1, t_VEC));
566 28 : if (d < 0) pari_err_ROOTS0("zncoppersmith");
567 28 : if (B && typ(B) != t_INT) B = gceil(B);
568 28 : if (abscmpiu(X, X_SMALL) <= 0)
569 7 : return gc_upto(av, do_exhaustive(P, N, itos(X), B));
570 :
571 21 : if (B && equalii(B,N)) B = NULL;
572 21 : if (B) bnd = 1; /* bnd-hack is only for the case B = N */
573 21 : cP = gel(P,d+2);
574 21 : if (!gequal1(cP))
575 : {
576 : GEN r, z;
577 14 : gel(P,d+2) = cP = bezout(cP, N, &z, &r);
578 35 : for (j = 0; j < d; j++) gel(P,j+2) = Fp_mul(gel(P,j+2), z, N);
579 14 : if (!is_pm1(cP))
580 : {
581 7 : P = Q_primitive_part(P, &cP);
582 7 : if (cP) { N = diviiexact(N,cP); B = gceil(gdiv(B, cP)); }
583 : }
584 : }
585 21 : if (DEBUGLEVEL >= 2) err_printf("Modified P: %Ps\n", P);
586 :
587 21 : choose_params(P, N, X, B, &delta, &t);
588 14 : if (DEBUGLEVEL >= 2)
589 0 : err_printf("Init: trying delta = %d, t = %d\n", delta, t);
590 : for(;;)
591 : {
592 14 : dim = d * delta + t;
593 : /* TODO: In case of failure do not recompute the full vector */
594 14 : Xpowers = (GEN*)new_chunk(dim + 1);
595 14 : Xpowers[0] = gen_1;
596 217 : for (j = 1; j <= dim; j++) Xpowers[j] = mulii(Xpowers[j-1], X);
597 :
598 : /* TODO: in case of failure, use the part of the matrix already computed */
599 14 : M = zeromatcopy(dim,dim);
600 :
601 : /* Rows of M correspond to the polynomials
602 : * N^delta, N^delta Xi, ... N^delta (Xi)^d-1,
603 : * N^(delta-1)P(Xi), N^(delta-1)XiP(Xi), ... N^(delta-1)P(Xi)(Xi)^d-1,
604 : * ...
605 : * P(Xi)^delta, XiP(Xi)^delta, ..., P(Xi)^delta(Xi)^t-1 */
606 42 : for (j = 1; j <= d; j++) gcoeff(M, j, j) = gel(Xpowers,j-1);
607 :
608 : /* P-part */
609 14 : if (delta) row = d + 1; else row = 0;
610 :
611 14 : Q = P;
612 70 : for (i = 1; i < delta; i++)
613 : {
614 182 : for (j = 0; j < d; j++,row++)
615 1239 : for (l = j + 1; l <= row; l++)
616 1113 : gcoeff(M, l, row) = mulii(Xpowers[l-1], gel(Q,l-j+1));
617 56 : Q = ZX_mul(Q, P);
618 : }
619 63 : for (j = 0; j < t; row++, j++)
620 490 : for (l = j + 1; l <= row; l++)
621 441 : gcoeff(M, l, row) = mulii(Xpowers[l-1], gel(Q,l-j+1));
622 :
623 : /* N-part */
624 14 : row = dim - t; N0 = N;
625 84 : while (row >= 1)
626 : {
627 224 : for (j = 0; j < d; j++,row--)
628 1421 : for (l = 1; l <= row; l++)
629 1267 : gcoeff(M, l, row) = mulii(gmael(M, row, l), N0);
630 70 : if (row >= 1) N0 = mulii(N0, N);
631 : }
632 : /* Z is the upper bound for the L^1 norm of the polynomial,
633 : ie. N^delta if B = N, B^delta otherwise */
634 14 : if (B) Z = powiu(B, delta); else Z = N0;
635 :
636 14 : if (DEBUGLEVEL >= 2)
637 : {
638 0 : if (DEBUGLEVEL >= 6) err_printf("Matrix to be reduced:\n%Ps\n", M);
639 0 : err_printf("Entering LLL\nbitsize bound: %ld\n", expi(Z));
640 0 : err_printf("expected shvector bitsize: %ld\n", expi(ZM_det_triangular(M))/dim);
641 : }
642 :
643 14 : sh = ZM_lll(M, 0.75, LLL_INPLACE);
644 : /* Take the first vector if it is non constant */
645 14 : short_pol = gel(sh,1);
646 14 : if (ZV_isscalar(short_pol)) short_pol = gel(sh, 2);
647 :
648 14 : nsp = gen_0;
649 217 : for (j = 1; j <= dim; j++) nsp = addii(nsp, absi_shallow(gel(short_pol,j)));
650 :
651 14 : if (DEBUGLEVEL >= 2)
652 : {
653 0 : err_printf("Candidate: %Ps\n", short_pol);
654 0 : err_printf("bitsize Norm: %ld\n", expi(nsp));
655 0 : err_printf("bitsize bound: %ld\n", expi(mului(bnd, Z)));
656 : }
657 14 : if (cmpii(nsp, mului(bnd, Z)) < 0) break; /* SUCCESS */
658 :
659 : /* Failed with the precomputed or supplied value */
660 0 : if (++t == d) { delta++; t = 1; }
661 0 : if (DEBUGLEVEL >= 2)
662 0 : err_printf("Increasing dim, delta = %d t = %d\n", delta, t);
663 : }
664 14 : bnd = itos(divii(nsp, Z)) + 1;
665 :
666 14 : while (!signe(gel(short_pol,dim))) dim--;
667 :
668 14 : R = cgetg(dim + 2, t_POL); R[1] = P[1];
669 217 : for (j = 1; j <= dim; j++)
670 203 : gel(R,j+1) = diviiexact(gel(short_pol,j), Xpowers[j-1]);
671 14 : gel(R,2) = subii(gel(R,2), mului(bnd - 1, N0));
672 :
673 14 : sol = cgetg(1, t_VEC);
674 84 : for (i = -bnd + 1; i < bnd; i++)
675 : {
676 70 : GEN r = nfrootsQ(R);
677 70 : if (DEBUGLEVEL >= 2) err_printf("Roots: %Ps\n", r);
678 91 : for (j = 1; j < lg(r); j++)
679 : {
680 21 : GEN z = gel(r,j);
681 21 : if (typ(z) == t_INT && sol_OK(FpX_eval(P,z,N), N,B))
682 14 : sol = shallowconcat(sol, z);
683 : }
684 70 : if (i < bnd) gel(R,2) = addii(gel(R,2), Z);
685 : }
686 14 : return gc_upto(av, ZV_sort_uniq(sol));
687 : }
688 :
689 : /********************************************************************/
690 : /** **/
691 : /** LINEAR & ALGEBRAIC DEPENDENCE **/
692 : /** **/
693 : /********************************************************************/
694 :
695 : static int
696 8123 : real_indep(GEN re, GEN im, long bit)
697 : {
698 8123 : GEN d = gsub(gmul(gel(re,1),gel(im,2)), gmul(gel(re,2),gel(im,1)));
699 8123 : return (!gequal0(d) && gexpo(d) > - bit);
700 : }
701 :
702 : GEN
703 15302 : lindepfull_bit(GEN x, long bit)
704 : {
705 15302 : long lx = lg(x), ly, i, j;
706 : GEN re, im, M;
707 :
708 15302 : if (! is_vec_t(typ(x))) pari_err_TYPE("lindep2",x);
709 15302 : if (lx <= 2)
710 : {
711 21 : if (lx == 2 && gequal0(x)) return matid(1);
712 14 : return NULL;
713 : }
714 15281 : re = real_i(x);
715 15281 : im = imag_i(x);
716 : /* independent over R ? */
717 15281 : if (lx == 3 && real_indep(re,im,bit)) return NULL;
718 15267 : if (gequal0(im)) im = NULL;
719 15267 : ly = im? lx+2: lx+1;
720 15267 : M = cgetg(lx,t_MAT);
721 60701 : for (i=1; i<lx; i++)
722 : {
723 45434 : GEN c = cgetg(ly,t_COL); gel(M,i) = c;
724 209394 : for (j=1; j<lx; j++) gel(c,j) = gen_0;
725 45434 : gel(c,i) = gen_1;
726 45434 : gel(c,lx) = gtrunc2n(gel(re,i), bit);
727 45434 : if (im) gel(c,lx+1) = gtrunc2n(gel(im,i), bit);
728 : }
729 15267 : return ZM_lll(M, 0.99, LLL_INPLACE);
730 : }
731 : GEN
732 3311 : lindep_bit(GEN x, long bit)
733 : {
734 3311 : pari_sp av = avma;
735 3311 : GEN v, M = lindepfull_bit(x,bit);
736 3311 : if (!M) retgc_const(av, cgetg(1, t_COL));
737 3283 : v = gel(M,1); setlg(v, lg(M));
738 3283 : return gc_GEN(av, v);
739 : }
740 : /* deprecated */
741 : GEN
742 112 : lindep2(GEN x, long dig)
743 : {
744 : long bit;
745 112 : if (dig < 0) pari_err_DOMAIN("lindep2", "accuracy", "<", gen_0, stoi(dig));
746 112 : if (dig) bit = (long) (dig/LOG10_2);
747 : else
748 : {
749 98 : bit = gprecision(x);
750 98 : if (!bit)
751 : {
752 35 : x = Q_primpart(x); /* left on stack */
753 35 : bit = 32 + gexpo(x);
754 : }
755 : else
756 63 : bit = (long)prec2nbits_mul(bit, 0.8);
757 : }
758 112 : return lindep_bit(x, bit);
759 : }
760 :
761 : /* x is a vector of elts of a p-adic field */
762 : GEN
763 28 : lindep_padic(GEN x)
764 : {
765 28 : long i, j, prec = LONG_MAX, nx = lg(x)-1, v;
766 28 : pari_sp av = avma;
767 28 : GEN p = NULL, pn, m, a;
768 :
769 28 : if (nx < 2) return cgetg(1,t_COL);
770 147 : for (i=1; i<=nx; i++)
771 : {
772 119 : GEN c = gel(x,i), q;
773 119 : if (typ(c) != t_PADIC) continue;
774 :
775 91 : j = precp(c); if (j < prec) prec = j;
776 91 : q = padic_p(c);
777 91 : if (!p) p = q; else if (!equalii(p, q)) pari_err_MODULUS("lindep_padic", p, q);
778 : }
779 28 : if (!p) pari_err_TYPE("lindep_padic [not a p-adic vector]",x);
780 28 : v = gvaluation(x,p); pn = powiu(p,prec);
781 28 : if (v) x = gmul(x, powis(p, -v));
782 28 : x = RgV_to_FpV(x, pn);
783 :
784 28 : a = negi(gel(x,1));
785 28 : m = cgetg(nx,t_MAT);
786 119 : for (i=1; i<nx; i++)
787 : {
788 91 : GEN c = zerocol(nx);
789 91 : gel(c,1+i) = a;
790 91 : gel(c,1) = gel(x,i+1);
791 91 : gel(m,i) = c;
792 : }
793 28 : m = ZM_lll(ZM_hnfmodid(m, pn), 0.99, LLL_INPLACE);
794 28 : return gc_GEN(av, gel(m,1));
795 : }
796 : /* x is a vector of t_POL/t_SER */
797 : GEN
798 77 : lindep_Xadic(GEN x)
799 : {
800 77 : long i, prec = LONG_MAX, deg = 0, lx = lg(x), vx, v;
801 77 : pari_sp av = avma;
802 : GEN m;
803 :
804 77 : if (lx == 1) return cgetg(1,t_COL);
805 77 : vx = gvar(x);
806 77 : if (gequal0(x)) return col_ei(lx-1,1);
807 70 : v = gvaluation(x, pol_x(vx));
808 70 : if (!v) x = shallowcopy(x);
809 0 : else if (v > 0) x = gdiv(x, pol_xn(v, vx));
810 0 : else x = gmul(x, pol_xn(-v, vx));
811 : /* all t_SER have valuation >= 0 */
812 735 : for (i=1; i<lx; i++)
813 : {
814 665 : GEN c = gel(x,i);
815 665 : if (gvar(c) != vx) { gel(x,i) = scalarpol_shallow(c, vx); continue; }
816 658 : switch(typ(c))
817 : {
818 231 : case t_POL: deg = maxss(deg, degpol(c)); break;
819 0 : case t_RFRAC: pari_err_TYPE("lindep_Xadic", c);
820 427 : case t_SER:
821 427 : prec = minss(prec, valser(c)+lg(c)-2);
822 427 : gel(x,i) = ser2rfrac_i(c);
823 : }
824 : }
825 70 : if (prec == LONG_MAX) prec = deg+1;
826 70 : m = RgXV_to_RgM(x, prec);
827 70 : return gc_upto(av, deplin(m));
828 : }
829 : static GEN
830 35 : vec_lindep(GEN x)
831 : {
832 35 : pari_sp av = avma;
833 35 : long i, l = lg(x); /* > 1 */
834 35 : long t = typ(gel(x,1)), h = lg(gel(x,1));
835 35 : GEN m = cgetg(l, t_MAT);
836 126 : for (i = 1; i < l; i++)
837 : {
838 98 : GEN y = gel(x,i);
839 98 : if (lg(y) != h || typ(y) != t) pari_err_TYPE("lindep",x);
840 91 : if (t != t_COL) y = shallowtrans(y); /* Sigh */
841 91 : gel(m,i) = y;
842 : }
843 28 : return gc_upto(av, deplin(m));
844 : }
845 :
846 : GEN
847 0 : lindep(GEN x) { return lindep2(x, 0); }
848 :
849 : GEN
850 434 : lindep0(GEN x,long bit)
851 : {
852 434 : long i, tx = typ(x);
853 434 : if (tx == t_MAT) return deplin(x);
854 147 : if (! is_vec_t(tx)) pari_err_TYPE("lindep",x);
855 441 : for (i = 1; i < lg(x); i++)
856 357 : switch(typ(gel(x,i)))
857 : {
858 7 : case t_PADIC: return lindep_padic(x);
859 21 : case t_POL:
860 : case t_RFRAC:
861 21 : case t_SER: return lindep_Xadic(x);
862 35 : case t_VEC:
863 35 : case t_COL: return vec_lindep(x);
864 : }
865 84 : return lindep2(x, bit);
866 : }
867 :
868 : GEN
869 77 : algdep0(GEN x, long n, long bit)
870 : {
871 77 : long tx = typ(x), i;
872 : pari_sp av;
873 : GEN y;
874 :
875 77 : if (! is_scalar_t(tx)) pari_err_TYPE("algdep0",x);
876 77 : if (tx == t_POLMOD)
877 : {
878 14 : av = avma; y = minpoly(x, 0);
879 14 : return (degpol(y) > n)? gc_const(av, gen_1): y;
880 : }
881 63 : if (gequal0(x)) return pol_x(0);
882 63 : if (n <= 0)
883 : {
884 14 : if (!n) return gen_1;
885 7 : pari_err_DOMAIN("algdep", "degree", "<", gen_0, stoi(n));
886 : }
887 :
888 49 : av = avma; y = cgetg(n+2,t_COL);
889 49 : gel(y,1) = gen_1;
890 49 : gel(y,2) = x; /* n >= 1 */
891 210 : for (i=3; i<=n+1; i++) gel(y,i) = gmul(gel(y,i-1),x);
892 49 : if (typ(x) == t_PADIC)
893 21 : y = lindep_padic(y);
894 : else
895 28 : y = lindep2(y, bit);
896 49 : if (lg(y) == 1) pari_err(e_DOMAIN,"algdep", "degree(x)",">", stoi(n), x);
897 49 : y = RgV_to_RgX(y, 0);
898 49 : if (signe(leading_coeff(y)) > 0) return gc_GEN(av, y);
899 14 : return gc_upto(av, ZX_neg(y));
900 : }
901 :
902 : GEN
903 0 : algdep(GEN x, long n)
904 : {
905 0 : return algdep0(x,n,0);
906 : }
907 :
908 : static GEN
909 56 : sertomat(GEN S, long p, long r, long vy)
910 : {
911 : long n, m;
912 56 : GEN v = cgetg(r*p+1, t_VEC); /* v[r*n+m+1] = s^n * y^m */
913 : /* n = 0 */
914 245 : for (m = 0; m < r; m++) gel(v, m+1) = pol_xn(m, vy);
915 175 : for(n=1; n < p; n++)
916 546 : for (m = 0; m < r; m++)
917 : {
918 427 : GEN c = gel(S,n);
919 427 : if (m)
920 : {
921 308 : c = shallowcopy(c);
922 308 : setvalser(c, valser(c) + m);
923 : }
924 427 : gel(v, r*n + m + 1) = c;
925 : }
926 56 : return v;
927 : }
928 :
929 : GEN
930 42 : seralgdep(GEN s, long p, long r)
931 : {
932 42 : pari_sp av = avma;
933 : long vy, i, n, prec;
934 : GEN S, v, D;
935 :
936 42 : if (typ(s) != t_SER) pari_err_TYPE("seralgdep",s);
937 42 : if (p <= 0) pari_err_DOMAIN("seralgdep", "p", "<=", gen_0, stoi(p));
938 42 : if (r < 0) pari_err_DOMAIN("seralgdep", "r", "<", gen_0, stoi(r));
939 42 : if (is_bigint(addiu(muluu(p, r), 1))) pari_err_OVERFLOW("seralgdep");
940 42 : vy = varn(s);
941 42 : if (!vy) pari_err_PRIORITY("seralgdep", s, ">", 0);
942 42 : r++; p++;
943 42 : prec = valser(s) + lg(s)-2;
944 42 : if (r > prec) r = prec;
945 42 : S = cgetg(p+1, t_VEC); gel(S, 1) = s;
946 119 : for (i = 2; i <= p; i++) gel(S,i) = gmul(gel(S,i-1), s);
947 42 : v = sertomat(S, p, r, vy);
948 42 : D = lindep_Xadic(v);
949 42 : if (lg(D) == 1) { set_avma(av); return gen_0; }
950 35 : v = cgetg(p+1, t_VEC);
951 133 : for (n = 0; n < p; n++)
952 98 : gel(v, n+1) = RgV_to_RgX(vecslice(D, r*n+1, r*n+r), vy);
953 35 : return gc_GEN(av, RgV_to_RgX(v, 0));
954 : }
955 :
956 : GEN
957 14 : serdiffdep(GEN s, long p, long r)
958 : {
959 14 : pari_sp av = avma;
960 : long vy, i, n, prec;
961 : GEN P, S, v, D;
962 :
963 14 : if (typ(s) != t_SER) pari_err_TYPE("serdiffdep",s);
964 14 : if (p <= 0) pari_err_DOMAIN("serdiffdep", "p", "<=", gen_0, stoi(p));
965 14 : if (r < 0) pari_err_DOMAIN("serdiffdep", "r", "<", gen_0, stoi(r));
966 14 : if (is_bigint(addiu(muluu(p, r), 1))) pari_err_OVERFLOW("serdiffdep");
967 14 : vy = varn(s);
968 14 : if (!vy) pari_err_PRIORITY("serdiffdep", s, ">", 0);
969 14 : r++; p++;
970 14 : prec = valser(s) + lg(s)-2;
971 14 : if (r > prec) r = prec;
972 14 : S = cgetg(p+1, t_VEC); gel(S, 1) = s;
973 56 : for (i = 2; i <= p; i++) gel(S,i) = derivser(gel(S,i-1));
974 14 : v = sertomat(S, p, r, vy);
975 14 : D = lindep_Xadic(v);
976 14 : if (lg(D) == 1) { set_avma(av); return gen_0; }
977 14 : P = RgV_to_RgX(vecslice(D, 1, r), vy);
978 14 : v = cgetg(p, t_VEC);
979 56 : for (n = 1; n < p; n++)
980 42 : gel(v, n) = RgV_to_RgX(vecslice(D, r*n+1, r*n+r), vy);
981 14 : return gc_GEN(av, mkvec2(RgV_to_RgX(v, 0), gneg(P)));
982 : }
983 :
984 : /* FIXME: could precompute ZM_lll attached to V[2..] */
985 : static GEN
986 11991 : lindepcx(GEN V, long bit)
987 : {
988 11991 : GEN Vr = real_i(V), Vi = imag_i(V);
989 11991 : long d = gexpo(Vr) - gexpo(Vi);
990 11991 : if (d < -bit) V = Vi;
991 11991 : else if (d > bit) V = Vr;
992 11991 : return lindepfull_bit(V, bit);
993 : }
994 : /* c floating point t_REAL or t_COMPLEX, T ZX, recognize in Q[x]/(T).
995 : * V helper vector (containing complex roots of T), MODIFIED */
996 : static GEN
997 11991 : cx_bestapprnf(GEN c, GEN T, GEN V, long bit)
998 : {
999 11991 : GEN M, a, v = NULL;
1000 : long i, l;
1001 11991 : gel(V,1) = gneg(c); M = lindepcx(V, bit);
1002 11991 : if (!M) pari_err(e_MISC, "cannot rationalize coeff in bestapprnf");
1003 11991 : l = lg(M); a = NULL;
1004 11991 : for (i = 1; i < l; i ++) { v = gel(M,i); a = gel(v,1); if (signe(a)) break; }
1005 11991 : v = RgC_Rg_div(vecslice(v, 2, lg(M)-1), a);
1006 11991 : if (!T) return gel(v,1);
1007 4830 : v = RgV_to_RgX(v, varn(T)); l = lg(v);
1008 4830 : if (l == 2) return gen_0;
1009 4165 : if (l == 3) return gel(v,2);
1010 3668 : return mkpolmod(v, T);
1011 : }
1012 : static GEN
1013 14784 : bestapprnf_i(GEN x, GEN T, GEN V, long bit)
1014 : {
1015 14784 : long i, l, tx = typ(x);
1016 : GEN z;
1017 14784 : switch (tx)
1018 : {
1019 1505 : case t_INT: case t_FRAC: return x;
1020 11991 : case t_REAL: case t_COMPLEX: return cx_bestapprnf(x, T, V, bit);
1021 0 : case t_POLMOD: if (RgX_equal(gel(x,1),T)) return x;
1022 0 : break;
1023 1288 : case t_POL: case t_SER: case t_VEC: case t_COL: case t_MAT:
1024 1288 : l = lg(x); z = cgetg(l, tx);
1025 1974 : for (i = 1; i < lontyp[tx]; i++) z[i] = x[i];
1026 13993 : for (; i < l; i++) gel(z,i) = bestapprnf_i(gel(x,i), T, V, bit);
1027 1288 : return z;
1028 : }
1029 0 : pari_err_TYPE("mfcxtoQ", x);
1030 : return NULL;/*LCOV_EXCL_LINE*/
1031 : }
1032 :
1033 : GEN
1034 2163 : bestapprnf(GEN x, GEN T, GEN roT, long prec)
1035 : {
1036 2163 : pari_sp av = avma;
1037 2163 : long tx = typ(x), dT = 1, bit;
1038 : GEN V;
1039 :
1040 2163 : if (T)
1041 : {
1042 1610 : if (typ(T) != t_POL) T = nf_get_pol(checknf(T));
1043 1610 : else if (!RgX_is_ZX(T)) pari_err_TYPE("bestapprnf", T);
1044 1610 : dT = degpol(T);
1045 : }
1046 2163 : if (is_rational_t(tx)) return gcopy(x);
1047 2079 : if (tx == t_POLMOD)
1048 : {
1049 0 : if (!T || !RgX_equal(T, gel(x,1))) pari_err_TYPE("bestapprnf",x);
1050 0 : return gcopy(x);
1051 : }
1052 :
1053 2079 : if (roT)
1054 : {
1055 644 : long l = gprecision(roT);
1056 644 : switch(typ(roT))
1057 : {
1058 644 : case t_INT: case t_FRAC: case t_REAL: case t_COMPLEX: break;
1059 0 : default: pari_err_TYPE("bestapprnf", roT);
1060 : }
1061 644 : if (prec < l) prec = l;
1062 : }
1063 1435 : else if (!T)
1064 525 : roT = gen_1;
1065 : else
1066 : {
1067 910 : long n = poliscyclo(T); /* cyclotomic is an important special case */
1068 910 : roT = n? rootsof1u_cx(n,prec): gel(QX_complex_roots(T,prec), 1);
1069 : }
1070 2079 : V = vec_prepend(gpowers(roT, dT-1), NULL);
1071 2079 : bit = prec2nbits_mul(prec, 0.8);
1072 2079 : return gc_GEN(av, bestapprnf_i(x, T, V, bit));
1073 : }
1074 :
1075 : /********************************************************************/
1076 : /** **/
1077 : /** MINIM **/
1078 : /** **/
1079 : /********************************************************************/
1080 : void
1081 124054 : minim_alloc(long n, double ***q, GEN *x, double **y, double **z, double **v)
1082 : {
1083 124054 : long i, s = n * sizeof(double);
1084 :
1085 124054 : *x = cgetg(n, t_VECSMALL);
1086 124054 : *q = (double**) new_chunk(n);
1087 124054 : *y = (double*) stack_malloc_align(s, sizeof(double));
1088 124054 : *z = (double*) stack_malloc_align(s, sizeof(double));
1089 124054 : *v = (double*) stack_malloc_align(s, sizeof(double));
1090 534974 : for (i=1; i<n; i++) (*q)[i] = (double*) stack_malloc_align(s, sizeof(double));
1091 124054 : }
1092 :
1093 : static void
1094 70 : cvp_alloc(long n, double **t, double **tpre)
1095 : {
1096 70 : long s = n * sizeof(double);
1097 70 : *t = (double*) stack_malloc_align(s, sizeof(double));
1098 70 : *tpre = (double*) stack_malloc_align(s, sizeof(double));
1099 70 : }
1100 :
1101 : static GEN
1102 5502 : ZC_canon(GEN V)
1103 : {
1104 5502 : long l = lg(V), j, s;
1105 11242 : for (j = 1; j < l; j++)
1106 11242 : if ((s = signe(gel(V,j)))) return s < 0? ZC_neg(V): V;
1107 0 : return V;
1108 : }
1109 : static GEN
1110 5502 : ZM_zc_mul_canon(GEN u, GEN x) { return ZC_canon(ZM_zc_mul(u,x)); }
1111 : static GEN
1112 240366 : ZM_zc_mul_canon_zm(GEN u, GEN x)
1113 : {
1114 240366 : pari_sp av = avma;
1115 240366 : GEN y = ZV_to_zv(ZM_zc_mul(u,x));
1116 240366 : zv_canon_inplace(y); return gc_upto(av, y);
1117 : }
1118 :
1119 : struct qfvec
1120 : {
1121 : GEN a, r, u;
1122 : };
1123 :
1124 : static void
1125 0 : err_minim(GEN a)
1126 : {
1127 0 : pari_err_DOMAIN("minim0","form","is not",
1128 : strtoGENstr("positive definite"),a);
1129 0 : }
1130 :
1131 : static GEN
1132 902 : minim_lll(GEN a, GEN *u)
1133 : {
1134 902 : *u = lllgramint(a);
1135 902 : if (lg(*u) != lg(a)) err_minim(a);
1136 902 : return qf_ZM_apply(a,*u);
1137 : }
1138 :
1139 : static void
1140 902 : forqfvec_init_dolll(struct qfvec *qv, GEN *pa, long dolll)
1141 : {
1142 902 : GEN r, u, a = *pa;
1143 902 : if (!dolll) u = NULL;
1144 : else
1145 : {
1146 860 : if (typ(a) != t_MAT || !RgM_is_ZM(a)) pari_err_TYPE("qfminim",a);
1147 860 : a = *pa = minim_lll(a, &u);
1148 : }
1149 902 : qv->a = RgM_gtofp(a, DEFAULTPREC);
1150 902 : r = qfgaussred_positive(qv->a);
1151 902 : if (!r)
1152 : {
1153 0 : r = qfgaussred_positive(a); /* exact computation */
1154 0 : if (!r) err_minim(a);
1155 0 : r = RgM_gtofp(r, DEFAULTPREC);
1156 : }
1157 902 : qv->r = r;
1158 902 : qv->u = u;
1159 902 : }
1160 :
1161 : static void
1162 42 : forqfvec_init(struct qfvec *qv, GEN a)
1163 42 : { forqfvec_init_dolll(qv, &a, 1); }
1164 :
1165 : static void
1166 42 : forqfvec_i(void *E, long (*fun)(void *, GEN, GEN, double), struct qfvec *qv, GEN BORNE)
1167 : {
1168 42 : GEN x, a = qv->a, r = qv->r, u = qv->u;
1169 42 : long n = lg(a)-1, i, j, k;
1170 : double p,BOUND,*v,*y,*z,**q;
1171 42 : const double eps = 1e-10;
1172 42 : if (!BORNE) BORNE = gen_0;
1173 : else
1174 : {
1175 28 : BORNE = gfloor(BORNE);
1176 28 : if (typ(BORNE) != t_INT) pari_err_TYPE("minim0",BORNE);
1177 35 : if (signe(BORNE) <= 0) return;
1178 : }
1179 35 : if (n == 0) return;
1180 28 : minim_alloc(n+1, &q, &x, &y, &z, &v);
1181 98 : for (j=1; j<=n; j++)
1182 : {
1183 70 : v[j] = rtodbl(gcoeff(r,j,j));
1184 133 : for (i=1; i<j; i++) q[i][j] = rtodbl(gcoeff(r,i,j));
1185 : }
1186 :
1187 28 : if (gequal0(BORNE))
1188 : {
1189 : double c;
1190 14 : p = rtodbl(gcoeff(a,1,1));
1191 42 : for (i=2; i<=n; i++) { c = rtodbl(gcoeff(a,i,i)); if (c < p) p = c; }
1192 14 : BORNE = roundr(dbltor(p));
1193 : }
1194 : else
1195 14 : p = gtodouble(BORNE);
1196 28 : BOUND = p * (1 + eps);
1197 28 : if (BOUND > (double)ULONG_MAX || (ulong)BOUND != (ulong)p)
1198 7 : pari_err_PREC("forqfvec");
1199 :
1200 21 : k = n; y[n] = z[n] = 0;
1201 21 : x[n] = (long)sqrt(BOUND/v[n]);
1202 56 : for(;;x[1]--)
1203 : {
1204 : do
1205 : {
1206 140 : if (k>1)
1207 : {
1208 84 : long l = k-1;
1209 84 : z[l] = 0;
1210 245 : for (j=k; j<=n; j++) z[l] += q[l][j]*x[j];
1211 84 : p = (double)x[k] + z[k];
1212 84 : y[l] = y[k] + p*p*v[k];
1213 84 : x[l] = (long)floor(sqrt((BOUND-y[l])/v[l])-z[l]);
1214 84 : k = l;
1215 : }
1216 : for(;;)
1217 : {
1218 189 : p = (double)x[k] + z[k];
1219 189 : if (y[k] + p*p*v[k] <= BOUND) break;
1220 49 : k++; x[k]--;
1221 : }
1222 140 : } while (k > 1);
1223 77 : if (! x[1] && y[1]<=eps) break;
1224 :
1225 56 : p = (double)x[1] + z[1]; p = y[1] + p*p*v[1]; /* norm(x) */
1226 56 : if (fun(E, u, x, p)) break;
1227 : }
1228 : }
1229 :
1230 : void
1231 0 : forqfvec(void *E, long (*fun)(void *, GEN, GEN, double), GEN a, GEN BORNE)
1232 : {
1233 0 : pari_sp av = avma;
1234 : struct qfvec qv;
1235 0 : forqfvec_init(&qv, a);
1236 0 : forqfvec_i(E, fun, &qv, BORNE);
1237 0 : set_avma(av);
1238 0 : }
1239 :
1240 : struct qfvecwrap
1241 : {
1242 : void *E;
1243 : long (*fun)(void *, GEN);
1244 : };
1245 :
1246 : static long
1247 56 : forqfvec_wrap(void *E, GEN u, GEN x, double d)
1248 : {
1249 56 : pari_sp av = avma;
1250 56 : struct qfvecwrap *W = (struct qfvecwrap *) E;
1251 : (void) d;
1252 56 : return gc_long(av, W->fun(W->E, ZM_zc_mul_canon(u, x)));
1253 : }
1254 :
1255 : void
1256 42 : forqfvec1(void *E, long (*fun)(void *, GEN), GEN a, GEN BORNE)
1257 : {
1258 42 : pari_sp av = avma;
1259 : struct qfvecwrap wr;
1260 : struct qfvec qv;
1261 42 : wr.E = E; wr.fun = fun;
1262 42 : forqfvec_init(&qv, a);
1263 42 : forqfvec_i((void*) &wr, forqfvec_wrap, &qv, BORNE);
1264 35 : set_avma(av);
1265 35 : }
1266 :
1267 : void
1268 42 : forqfvec0(GEN a, GEN BORNE, GEN code)
1269 42 : { EXPRVOID_WRAP(code, forqfvec1(EXPR_ARGVOID, a, BORNE)) }
1270 :
1271 : enum { min_ALL = 0, min_FIRST, min_VECSMALL, min_VECSMALL2 };
1272 :
1273 : static int
1274 923 : stockmax_init(const char *fun, GEN STOCKMAX, long *maxrank)
1275 : {
1276 923 : long r = 200;
1277 923 : if (!STOCKMAX) { *maxrank = 200; return 1; }
1278 511 : STOCKMAX = gfloor(STOCKMAX);
1279 511 : if (typ(STOCKMAX) != t_INT) pari_err_TYPE(fun, STOCKMAX);
1280 511 : r = itos(STOCKMAX);
1281 511 : if (r < 0)
1282 : {
1283 0 : char *e = stack_strcat(fun, "[negative number of vectors]");
1284 0 : pari_err_TYPE(e, STOCKMAX);
1285 : }
1286 511 : *maxrank = r; return 0;
1287 : }
1288 :
1289 : /* Minimal vectors for the integral definite quadratic form: a.
1290 : * Result u:
1291 : * u[1]= Number of vectors of square norm <= BORNE
1292 : * u[2]= maximum norm found
1293 : * u[3]= list of vectors found (at most STOCKMAX, unless NULL)
1294 : *
1295 : * If BORNE = NULL: Minimal nonzero vectors.
1296 : * flag = min_ALL, as above
1297 : * flag = min_FIRST, exits when first suitable vector is found.
1298 : * flag = min_VECSMALL, return a t_VECSMALL of (half) the number of vectors
1299 : * for each norm
1300 : * flag = min_VECSMALL2, same but count only vectors with even norm, and shift
1301 : * the answer */
1302 : static GEN
1303 847 : minim0_dolll(GEN a, GEN BORNE, GEN STOCKMAX, long flag, long dolll)
1304 : {
1305 : GEN x, u, r, L, gnorme;
1306 847 : long n = lg(a)-1, i, j, k, s, maxrank, sBORNE;
1307 847 : pari_sp av = avma, av1;
1308 : double p,maxnorm,BOUND,*v,*y,*z,**q;
1309 847 : const double eps = 1e-10;
1310 : int stockall;
1311 : struct qfvec qv;
1312 :
1313 847 : if (!BORNE)
1314 56 : sBORNE = 0;
1315 : else
1316 : {
1317 791 : BORNE = gfloor(BORNE);
1318 791 : if (typ(BORNE) != t_INT) pari_err_TYPE("minim0",BORNE);
1319 791 : if (is_bigint(BORNE)) pari_err_PREC( "qfminim");
1320 790 : sBORNE = itos(BORNE); set_avma(av);
1321 790 : if (sBORNE < 0) sBORNE = 0;
1322 : }
1323 846 : stockall = stockmax_init("minim0", STOCKMAX, &maxrank);
1324 :
1325 846 : switch(flag)
1326 : {
1327 462 : case min_VECSMALL:
1328 : case min_VECSMALL2:
1329 462 : if (sBORNE <= 0) return cgetg(1, t_VECSMALL);
1330 434 : L = zero_zv(sBORNE);
1331 434 : if (flag == min_VECSMALL2) sBORNE <<= 1;
1332 434 : if (n == 0) return L;
1333 434 : break;
1334 35 : case min_FIRST:
1335 35 : if (n == 0 || (!sBORNE && BORNE)) return cgetg(1,t_VEC);
1336 21 : L = NULL; /* gcc -Wall */
1337 21 : break;
1338 349 : case min_ALL:
1339 349 : if (n == 0 || (!sBORNE && BORNE))
1340 14 : retmkvec3(gen_0, gen_0, cgetg(1, t_MAT));
1341 335 : L = new_chunk(1+maxrank);
1342 335 : break;
1343 0 : default:
1344 0 : return NULL;
1345 : }
1346 790 : minim_alloc(n+1, &q, &x, &y, &z, &v);
1347 :
1348 790 : forqfvec_init_dolll(&qv, &a, dolll);
1349 790 : av1 = avma;
1350 790 : r = qv.r;
1351 790 : u = qv.u;
1352 5912 : for (j=1; j<=n; j++)
1353 : {
1354 5122 : v[j] = rtodbl(gcoeff(r,j,j));
1355 29579 : for (i=1; i<j; i++) q[i][j] = rtodbl(gcoeff(r,i,j)); /* |.| <= 1/2 */
1356 : }
1357 :
1358 790 : if (sBORNE) maxnorm = 0.;
1359 : else
1360 : {
1361 56 : GEN B = gcoeff(a,1,1);
1362 56 : long t = 1;
1363 616 : for (i=2; i<=n; i++)
1364 : {
1365 560 : GEN c = gcoeff(a,i,i);
1366 560 : if (cmpii(c, B) < 0) { B = c; t = i; }
1367 : }
1368 56 : if (flag == min_FIRST) return gc_GEN(av, mkvec2(B, gel(u,t)));
1369 49 : maxnorm = -1.; /* don't update maxnorm */
1370 49 : if (is_bigint(B)) return NULL;
1371 48 : sBORNE = itos(B);
1372 : }
1373 782 : BOUND = sBORNE * (1 + eps);
1374 782 : if ((long)BOUND != sBORNE) return NULL;
1375 :
1376 770 : s = 0;
1377 770 : k = n; y[n] = z[n] = 0;
1378 770 : x[n] = (long)sqrt(BOUND/v[n]);
1379 1223264 : for(;;x[1]--)
1380 : {
1381 : do
1382 : {
1383 2245614 : if (k>1)
1384 : {
1385 1022259 : long l = k-1;
1386 1022259 : z[l] = 0;
1387 11756080 : for (j=k; j<=n; j++) z[l] += q[l][j]*x[j];
1388 1022259 : p = (double)x[k] + z[k];
1389 1022259 : y[l] = y[k] + p*p*v[k];
1390 1022259 : x[l] = (long)floor(sqrt((BOUND-y[l])/v[l])-z[l]);
1391 1022259 : k = l;
1392 : }
1393 : for(;;)
1394 : {
1395 3263729 : p = (double)x[k] + z[k];
1396 3263729 : if (y[k] + p*p*v[k] <= BOUND) break;
1397 1018115 : k++; x[k]--;
1398 : }
1399 : }
1400 2245614 : while (k > 1);
1401 1224034 : if (! x[1] && y[1]<=eps) break;
1402 :
1403 1223271 : p = (double)x[1] + z[1];
1404 1223271 : p = y[1] + p*p*v[1]; /* norm(x) */
1405 1223271 : if (maxnorm >= 0)
1406 : {
1407 1220723 : if (p > maxnorm) maxnorm = p;
1408 : }
1409 : else
1410 : { /* maxnorm < 0 : only look for minimal vectors */
1411 2548 : pari_sp av2 = avma;
1412 2548 : gnorme = roundr(dbltor(p));
1413 2548 : if (cmpis(gnorme, sBORNE) >= 0) set_avma(av2);
1414 : else
1415 : {
1416 14 : sBORNE = itos(gnorme); set_avma(av1);
1417 14 : BOUND = sBORNE * (1+eps);
1418 14 : L = new_chunk(maxrank+1);
1419 14 : s = 0;
1420 : }
1421 : }
1422 1223271 : s++;
1423 :
1424 1223271 : switch(flag)
1425 : {
1426 7 : case min_FIRST:
1427 7 : if (dolll) x = ZM_zc_mul_canon(u,x);
1428 7 : return gc_GEN(av, mkvec2(roundr(dbltor(p)), x));
1429 :
1430 248241 : case min_ALL:
1431 248241 : if (s > maxrank && stockall) /* overflow */
1432 : {
1433 490 : long maxranknew = maxrank << 1;
1434 490 : GEN Lnew = new_chunk(1 + maxranknew);
1435 344890 : for (i=1; i<=maxrank; i++) Lnew[i] = L[i];
1436 490 : L = Lnew; maxrank = maxranknew;
1437 : }
1438 248241 : if (s<=maxrank) gel(L,s) = leafcopy(x);
1439 248241 : break;
1440 :
1441 39200 : case min_VECSMALL:
1442 39200 : { ulong norm = (ulong)(p + 0.5); L[norm]++; }
1443 39200 : break;
1444 :
1445 935823 : case min_VECSMALL2:
1446 935823 : { ulong norm = (ulong)(p + 0.5); if (!odd(norm)) L[norm>>1]++; }
1447 935823 : break;
1448 :
1449 : }
1450 : }
1451 763 : switch(flag)
1452 : {
1453 7 : case min_FIRST:
1454 7 : retgc_const(av, cgetg(1, t_VEC));
1455 434 : case min_VECSMALL:
1456 : case min_VECSMALL2:
1457 434 : set_avma((pari_sp)L); return L;
1458 : }
1459 322 : r = (maxnorm >= 0) ? roundr(dbltor(maxnorm)): stoi(sBORNE);
1460 322 : k = minss(s,maxrank);
1461 322 : L[0] = evaltyp(t_MAT) | evallg(k + 1);
1462 322 : if (dolll)
1463 246092 : for (j=1; j<=k; j++)
1464 245805 : gel(L,j) = dolll==1 ? ZM_zc_mul_canon(u, gel(L,j))
1465 245805 : : ZM_zc_mul_canon_zm(u, gel(L,j));
1466 322 : return gc_GEN(av, mkvec3(stoi(s<<1), r, L));
1467 : }
1468 :
1469 : /* Closest vectors for the integral definite quadratic form: a.
1470 : * Code bases on minim0_dolll
1471 : * Result u:
1472 : * u[1]= Number of closest vectors of square distance <= BORNE
1473 : * u[2]= maximum squared distance found
1474 : * u[3]= list of vectors found (at most STOCKMAX, unless NULL)
1475 : *
1476 : * If BORNE = NULL or <= 0.: returns closest vectors.
1477 : * flag = min_ALL, as above
1478 : * flag = min_FIRST, exits when first suitable vector is found.
1479 : */
1480 : static GEN
1481 91 : cvp0_dolll(GEN a, GEN target, GEN BORNE, GEN STOCKMAX, long flag, long dolll)
1482 : {
1483 : GEN x, u, r, L;
1484 91 : long n = lg(a)-1, i, j, k, s, maxrank;
1485 91 : pari_sp av = avma, av1;
1486 : double p,maxnorm,BOUND,*v,*y,*z,*tt,**q, *tpre, sBORNE;
1487 91 : const double eps = 1e-10;
1488 : int stockall;
1489 : struct qfvec qv;
1490 91 : int done = 0;
1491 :
1492 91 : if (!is_vec_t(typ(target))) pari_err_TYPE("cvp0",target);
1493 91 : if (n != lg(target)-1) pari_err_TYPE("cvp0 [different dimensions]",target);
1494 77 : if (!BORNE)
1495 0 : sBORNE = 0.;
1496 : else
1497 : {
1498 77 : if (!is_real_t(typ(BORNE))) pari_err_TYPE("cvp0",BORNE);
1499 77 : sBORNE = gtodouble(BORNE);
1500 77 : if (sBORNE < 0.) sBORNE = 0.;
1501 : }
1502 77 : stockall = stockmax_init("cvp0", STOCKMAX, &maxrank);
1503 :
1504 77 : L = (flag==min_ALL) ? new_chunk(1+maxrank) : NULL;
1505 77 : if (n == 0)
1506 : {
1507 7 : if (flag==min_ALL) retmkvec3(gen_0, gen_0, cgetg(1, t_MAT));
1508 0 : return cgetg(1,t_VEC);
1509 : }
1510 :
1511 70 : minim_alloc(n+1, &q, &x, &y, &z, &v);
1512 70 : cvp_alloc(n+1, &tt, &tpre);
1513 :
1514 70 : forqfvec_init_dolll(&qv, &a, dolll);
1515 70 : av1 = avma;
1516 70 : r = qv.r;
1517 70 : u = qv.u;
1518 392 : for (j=1; j<=n; j++)
1519 : {
1520 322 : v[j] = rtodbl(gcoeff(r,j,j));
1521 1729 : for (i=1; i<j; i++) q[i][j] = rtodbl(gcoeff(r,i,j)); /* |.| <= 1/2 */
1522 : }
1523 :
1524 70 : if (dolll)
1525 : {
1526 70 : GEN tv = RgM_RgC_mul(ZM_inv(u, NULL), target);
1527 392 : for (j=1; j<=n; j++) tt[j] = gtodouble(gel(tv, j));
1528 : } else
1529 0 : for (j=1; j<=n; j++) tt[j] = gtodouble(gel(target, j));
1530 : /* precompute contribution of tt to z[l] */
1531 392 : for(k=1; k <= n; k++)
1532 : {
1533 322 : tpre[k] = -tt[k];
1534 1729 : for(j=k+1; j<=n; j++) tpre[k] -= q[k][j] * tt[j];
1535 : }
1536 :
1537 70 : if (sBORNE) maxnorm = 0.;
1538 : else
1539 : {
1540 28 : GEN B = gcoeff(a,1,1);
1541 112 : for (i = 2; i <= n; i++) B = addii(B, gcoeff(a,i,i));
1542 28 : maxnorm = -1.; /* don't update maxnorm */
1543 28 : if (is_bigint(B)) return NULL;
1544 28 : sBORNE = 0.;
1545 140 : for(i=1; i<=n; i++) sBORNE += v[i];
1546 : }
1547 70 : BOUND = sBORNE * (1 + eps);
1548 :
1549 70 : s = 0;
1550 70 : k = n; y[n] = 0;
1551 70 : z[n] = tpre[n];
1552 70 : x[n] = (long)floor(sqrt(BOUND/v[n])-z[n]);
1553 889 : for(;;x[1]--)
1554 : {
1555 : do
1556 : {
1557 8582 : if (k>1)
1558 : {
1559 7665 : long l = k-1;
1560 7665 : z[l] = tpre[l];
1561 61488 : for (j=k; j<=n; j++) z[l] += q[l][j]*x[j];
1562 7665 : p = (double)x[k] + z[k];
1563 7665 : y[l] = y[k] + p*p*v[k];
1564 7665 : x[l] = (long)floor(sqrt((BOUND-y[l])/v[l])-z[l]);
1565 7665 : k = l;
1566 : }
1567 : for(;;)
1568 : {
1569 16247 : p = (double)x[k] + z[k];
1570 16247 : if (y[k] + p*p*v[k] <= BOUND) break;
1571 7735 : if (k >= n) { done = 1; break; }
1572 7665 : k++; x[k]--;
1573 : }
1574 : }
1575 8582 : while (k > 1 && !done);
1576 959 : if (done) break;
1577 :
1578 889 : p = (double)x[1] + z[1];
1579 889 : p = y[1] + p*p*v[1]; /* norm(x-target) */
1580 889 : if (maxnorm >= 0)
1581 : {
1582 175 : if (p > maxnorm) maxnorm = p;
1583 : }
1584 : else
1585 : { /* maxnorm < 0 : only look for closest vectors */
1586 714 : if (p * (1+10*eps) < sBORNE) {
1587 231 : sBORNE = p; set_avma(av1);
1588 231 : BOUND = sBORNE * (1+eps);
1589 231 : L = new_chunk(maxrank+1);
1590 231 : s = 0;
1591 : }
1592 : }
1593 889 : s++;
1594 :
1595 889 : switch(flag)
1596 : {
1597 0 : case min_FIRST:
1598 0 : if (dolll) x = ZM_zc_mul(u,x);
1599 0 : return gc_GEN(av, mkvec2(dbltor(p), x));
1600 :
1601 889 : case min_ALL:
1602 889 : if (s > maxrank && stockall) /* overflow */
1603 : {
1604 0 : long maxranknew = maxrank << 1;
1605 0 : GEN Lnew = new_chunk(1 + maxranknew);
1606 0 : for (i=1; i<=maxrank; i++) Lnew[i] = L[i];
1607 0 : L = Lnew; maxrank = maxranknew;
1608 : }
1609 889 : if (s<=maxrank) gel(L,s) = leafcopy(x);
1610 889 : break;
1611 : }
1612 : }
1613 70 : switch(flag)
1614 : {
1615 0 : case min_FIRST:
1616 0 : retgc_const(av, cgetg(1, t_VEC));
1617 : }
1618 70 : r = (maxnorm >= 0) ? dbltor(maxnorm): dbltor(sBORNE);
1619 70 : k = minss(s,maxrank);
1620 70 : L[0] = evaltyp(t_MAT) | evallg(k + 1);
1621 322 : for (j=1; j<=k; j++)
1622 252 : gel(L,j) = dolll==1 ? ZM_zc_mul(u, gel(L,j))
1623 252 : : zc_to_ZC(gel(L,j));
1624 70 : return gc_GEN(av, mkvec3(stoi(s), r, L));
1625 : }
1626 :
1627 : static GEN
1628 553 : minim0(GEN a, GEN BORNE, GEN STOCKMAX, long flag)
1629 : {
1630 553 : GEN v = minim0_dolll(a, BORNE, STOCKMAX, flag, 1);
1631 552 : if (!v) pari_err_PREC("qfminim");
1632 546 : return v;
1633 : }
1634 :
1635 : static GEN
1636 91 : cvp0(GEN a, GEN target, GEN BORNE, GEN STOCKMAX, long flag)
1637 : {
1638 91 : GEN v = cvp0_dolll(a, target, BORNE, STOCKMAX, flag, 1);
1639 77 : if (!v) pari_err_PREC("qfcvp");
1640 77 : return v;
1641 : }
1642 :
1643 : static GEN
1644 252 : minim0_zm(GEN a, GEN BORNE, GEN STOCKMAX, long flag)
1645 : {
1646 252 : GEN v = minim0_dolll(a, BORNE, STOCKMAX, flag, 2);
1647 252 : if (!v) pari_err_PREC("qfminim");
1648 252 : return v;
1649 : }
1650 :
1651 : GEN
1652 462 : qfrep0(GEN a, GEN borne, long flag)
1653 462 : { return minim0(a, borne, gen_0, (flag & 1)? min_VECSMALL2: min_VECSMALL); }
1654 :
1655 : GEN
1656 133 : qfminim0(GEN a, GEN borne, GEN stockmax, long flag, long prec)
1657 : {
1658 133 : switch(flag)
1659 : {
1660 49 : case 0: return minim0(a,borne,stockmax,min_ALL);
1661 35 : case 1: return minim0(a,borne,gen_0 ,min_FIRST);
1662 49 : case 2:
1663 : {
1664 49 : long maxnum = -1;
1665 49 : if (typ(a) != t_MAT) pari_err_TYPE("qfminim",a);
1666 49 : if (stockmax) {
1667 14 : if (typ(stockmax) != t_INT) pari_err_TYPE("qfminim",stockmax);
1668 14 : maxnum = itos(stockmax);
1669 : }
1670 49 : a = fincke_pohst(a,borne,maxnum,prec,NULL);
1671 42 : if (!a) pari_err_PREC("qfminim");
1672 42 : return a;
1673 : }
1674 0 : default: pari_err_FLAG("qfminim");
1675 : }
1676 : return NULL; /* LCOV_EXCL_LINE */
1677 : }
1678 :
1679 :
1680 : GEN
1681 91 : qfcvp0(GEN a, GEN target, GEN borne, GEN stockmax, long flag)
1682 : {
1683 91 : switch(flag)
1684 : {
1685 91 : case 0: return cvp0(a,target,borne,stockmax,min_ALL);
1686 0 : case 1: return cvp0(a,target,borne,gen_0 ,min_FIRST);
1687 : /* case 2:
1688 : TODO: more robust finke_pohst enumeration */
1689 0 : default: pari_err_FLAG("qfcvp");
1690 : }
1691 : return NULL; /* LCOV_EXCL_LINE */
1692 : }
1693 :
1694 : GEN
1695 7 : minim(GEN a, GEN borne, GEN stockmax)
1696 7 : { return minim0(a,borne,stockmax,min_ALL); }
1697 :
1698 : GEN
1699 252 : minim_zm(GEN a, GEN borne, GEN stockmax)
1700 252 : { return minim0_zm(a,borne,stockmax,min_ALL); }
1701 :
1702 : GEN
1703 42 : minim_raw(GEN a, GEN BORNE, GEN STOCKMAX)
1704 42 : { return minim0_dolll(a, BORNE, STOCKMAX, min_ALL, 0); }
1705 :
1706 : GEN
1707 0 : minim2(GEN a, GEN borne, GEN stockmax)
1708 0 : { return minim0(a,borne,stockmax,min_FIRST); }
1709 :
1710 : /* If V depends linearly from the columns of the matrix, return 0.
1711 : * Otherwise, update INVP and L and return 1. No GC. */
1712 : static int
1713 1652 : addcolumntomatrix(GEN V, GEN invp, GEN L)
1714 : {
1715 1652 : long i,j,k, n = lg(invp);
1716 1652 : GEN a = cgetg(n, t_COL), ak = NULL, mak;
1717 :
1718 84231 : for (k = 1; k < n; k++)
1719 83706 : if (!L[k])
1720 : {
1721 27902 : ak = RgMrow_zc_mul(invp, V, k);
1722 27902 : if (!gequal0(ak)) break;
1723 : }
1724 1652 : if (k == n) return 0;
1725 1127 : L[k] = 1;
1726 1127 : mak = gneg_i(ak);
1727 43253 : for (i=k+1; i<n; i++)
1728 42126 : gel(a,i) = gdiv(RgMrow_zc_mul(invp, V, i), mak);
1729 43883 : for (j=1; j<=k; j++)
1730 : {
1731 42756 : GEN c = gel(invp,j), ck = gel(c,k);
1732 42756 : if (gequal0(ck)) continue;
1733 8757 : gel(c,k) = gdiv(ck, ak);
1734 8757 : if (j==k)
1735 43253 : for (i=k+1; i<n; i++)
1736 42126 : gel(c,i) = gmul(gel(a,i), ck);
1737 : else
1738 184814 : for (i=k+1; i<n; i++)
1739 177184 : gel(c,i) = gadd(gel(c,i), gmul(gel(a,i), ck));
1740 : }
1741 1127 : return 1;
1742 : }
1743 :
1744 : GEN
1745 42 : qfperfection(GEN a)
1746 : {
1747 42 : pari_sp av = avma;
1748 : GEN u, L;
1749 42 : long r, s, k, l, n = lg(a)-1;
1750 :
1751 42 : if (!n) return gen_0;
1752 42 : if (typ(a) != t_MAT || !RgM_is_ZM(a)) pari_err_TYPE("qfperfection",a);
1753 42 : a = minim_lll(a, &u);
1754 42 : L = minim_raw(a,NULL,NULL);
1755 42 : r = (n*(n+1)) >> 1;
1756 42 : if (L)
1757 : {
1758 : GEN D, V, invp;
1759 35 : L = gel(L, 3); l = lg(L);
1760 35 : if (l == 2) { set_avma(av); return gen_1; }
1761 : /* |L[i]|^2 fits into a long for all i */
1762 21 : D = zero_zv(r);
1763 21 : V = cgetg(r+1, t_VECSMALL);
1764 21 : invp = matid(r);
1765 21 : s = 0;
1766 1659 : for (k = 1; k < l; k++)
1767 : {
1768 1652 : pari_sp av2 = avma;
1769 1652 : GEN x = gel(L,k);
1770 : long i, j, I;
1771 21098 : for (i = I = 1; i<=n; i++)
1772 145278 : for (j=i; j<=n; j++,I++) V[I] = x[i]*x[j];
1773 1652 : if (!addcolumntomatrix(V,invp,D)) set_avma(av2);
1774 1127 : else if (++s == r) break;
1775 : }
1776 : }
1777 : else
1778 : {
1779 : GEN M;
1780 7 : L = fincke_pohst(a,NULL,-1, DEFAULTPREC, NULL);
1781 7 : if (!L) pari_err_PREC("qfminim");
1782 7 : L = gel(L, 3); l = lg(L);
1783 7 : if (l == 2) { set_avma(av); return gen_1; }
1784 7 : M = cgetg(l, t_MAT);
1785 959 : for (k = 1; k < l; k++)
1786 : {
1787 952 : GEN x = gel(L,k), c = cgetg(r+1, t_COL);
1788 : long i, I, j;
1789 16184 : for (i = I = 1; i<=n; i++)
1790 144704 : for (j=i; j<=n; j++,I++) gel(c,I) = mulii(gel(x,i), gel(x,j));
1791 952 : gel(M,k) = c;
1792 : }
1793 7 : s = ZM_rank(M);
1794 : }
1795 28 : return gc_utoipos(av, s);
1796 : }
1797 :
1798 : static GEN
1799 140 : clonefill(GEN S, long s, long t)
1800 : { /* initialize to dummy values */
1801 140 : GEN T = S, dummy = cgetg(1, t_STR);
1802 : long i;
1803 308822 : for (i = s+1; i <= t; i++) gel(S,i) = dummy;
1804 140 : S = gclone(S); if (isclone(T)) gunclone(T);
1805 140 : return S;
1806 : }
1807 :
1808 : /* increment ZV x, by incrementing cell of index k. Initial value x0[k] was
1809 : * chosen to minimize qf(x) for given x0[1..k-1] and x0[k+1,..] = 0
1810 : * The last nonzero entry must be positive and goes through x0[k]+1,2,3,...
1811 : * Others entries go through: x0[k]+1,-1,2,-2,...*/
1812 : INLINE void
1813 2950008 : step(GEN x, GEN y, GEN inc, long k)
1814 : {
1815 2950008 : if (!signe(gel(y,k))) /* x[k+1..] = 0 */
1816 160807 : gel(x,k) = addiu(gel(x,k), 1); /* leading coeff > 0 */
1817 : else
1818 : {
1819 2789201 : long i = inc[k];
1820 2789201 : gel(x,k) = addis(gel(x,k), i),
1821 2789201 : inc[k] = (i > 0)? -1-i: 1-i;
1822 : }
1823 2950008 : }
1824 :
1825 : /* 1 if we are "sure" that x < y, up to few rounding errors, i.e.
1826 : * x < y - epsilon. More precisely :
1827 : * - sign(x - y) < 0
1828 : * - lgprec(x-y) > 3 || expo(x - y) - expo(x) > -24 */
1829 : static int
1830 1216062 : mplessthan(GEN x, GEN y)
1831 : {
1832 1216062 : pari_sp av = avma;
1833 1216062 : GEN z = mpsub(x, y);
1834 1216062 : set_avma(av);
1835 1216062 : if (typ(z) == t_INT) return (signe(z) < 0);
1836 1216062 : if (signe(z) >= 0) return 0;
1837 22371 : if (realprec(z) > LOWDEFAULTPREC) return 1;
1838 22371 : return ( expo(z) - mpexpo(x) > -24 );
1839 : }
1840 :
1841 : /* 1 if we are "sure" that x > y, up to few rounding errors, i.e.
1842 : * x > y + epsilon */
1843 : static int
1844 4616418 : mpgreaterthan(GEN x, GEN y)
1845 : {
1846 4616418 : pari_sp av = avma;
1847 4616418 : GEN z = mpsub(x, y);
1848 4616418 : set_avma(av);
1849 4616418 : if (typ(z) == t_INT) return (signe(z) > 0);
1850 4616418 : if (signe(z) <= 0) return 0;
1851 2690107 : if (realprec(z) > LOWDEFAULTPREC) return 1;
1852 476880 : return ( expo(z) - mpexpo(x) > -24 );
1853 : }
1854 :
1855 : /* x a t_INT, y t_INT or t_REAL */
1856 : INLINE GEN
1857 1228154 : mulimp(GEN x, GEN y)
1858 : {
1859 1228154 : if (typ(y) == t_INT) return mulii(x,y);
1860 1228154 : return signe(x) ? mulir(x,y): gen_0;
1861 : }
1862 : /* x + y*z, x,z two mp's, y a t_INT */
1863 : INLINE GEN
1864 13536473 : addmulimp(GEN x, GEN y, GEN z)
1865 : {
1866 13536473 : if (!signe(y)) return x;
1867 5830244 : if (typ(z) == t_INT) return mpadd(x, mulii(y, z));
1868 5830244 : return mpadd(x, mulir(y, z));
1869 : }
1870 :
1871 : /* yk + vk * (xk + zk)^2 */
1872 : static GEN
1873 5774835 : norm_aux(GEN xk, GEN yk, GEN zk, GEN vk)
1874 : {
1875 5774835 : GEN t = mpadd(xk, zk);
1876 5774835 : if (typ(t) == t_INT) { /* probably gen_0, avoid loss of accuracy */
1877 306205 : yk = addmulimp(yk, sqri(t), vk);
1878 : } else {
1879 5468630 : yk = mpadd(yk, mpmul(sqrr(t), vk));
1880 : }
1881 5774835 : return yk;
1882 : }
1883 : /* yk + vk * (xk + zk)^2 < B + epsilon */
1884 : static int
1885 4164187 : check_bound(GEN B, GEN xk, GEN yk, GEN zk, GEN vk)
1886 : {
1887 4164187 : pari_sp av = avma;
1888 4164187 : int f = mpgreaterthan(norm_aux(xk,yk,zk,vk), B);
1889 4164187 : return gc_bool(av, !f);
1890 : }
1891 :
1892 : /* q(k-th canonical basis vector), where q is given in Cholesky form
1893 : * q(x) = sum_{i = 1}^n q[i,i] (x[i] + sum_{j > i} q[i,j] x[j])^2.
1894 : * Namely q(e_k) = q[k,k] + sum_{i < k} q[i,i] q[i,k]^2
1895 : * Assume 1 <= k <= n. */
1896 : static GEN
1897 182 : cholesky_norm_ek(GEN q, long k)
1898 : {
1899 182 : GEN t = gcoeff(q,k,k);
1900 : long i;
1901 1484 : for (i = 1; i < k; i++) t = norm_aux(gen_0, t, gcoeff(q,i,k), gcoeff(q,i,i));
1902 182 : return t;
1903 : }
1904 :
1905 : /* q is the Cholesky decomposition of a quadratic form
1906 : * Enumerate vectors whose norm is less than BORNE (Algo 2.5.7),
1907 : * minimal vectors if BORNE = NULL (implies check = NULL).
1908 : * If (check != NULL) consider only vectors passing the check, and assumes
1909 : * we only want the smallest possible vectors */
1910 : static GEN
1911 14713 : smallvectors(GEN q, GEN BORNE, long maxnum, FP_chk_fun *CHECK)
1912 : {
1913 14713 : long N = lg(q), n = N-1, i, j, k, s, stockmax, checkcnt = 1;
1914 : pari_sp av, av1;
1915 : GEN inc, S, x, y, z, v, p1, alpha, norms;
1916 : GEN norme1, normax1, borne1, borne2;
1917 14713 : GEN (*check)(void *,GEN) = CHECK? CHECK->f: NULL;
1918 14713 : void *data = CHECK? CHECK->data: NULL;
1919 14713 : const long skipfirst = CHECK? CHECK->skipfirst: 0;
1920 14713 : const int stockall = (maxnum == -1);
1921 :
1922 14713 : alpha = dbltor(0.95);
1923 14713 : normax1 = gen_0;
1924 :
1925 14713 : v = cgetg(N,t_VEC);
1926 14713 : inc = const_vecsmall(n, 1);
1927 :
1928 14713 : av = avma;
1929 14713 : stockmax = stockall? 2000: maxnum;
1930 14713 : norms = cgetg(check?(stockmax+1): 1,t_VEC); /* unused if (!check) */
1931 14713 : S = cgetg(stockmax+1,t_VEC);
1932 14713 : x = cgetg(N,t_COL);
1933 14713 : y = cgetg(N,t_COL);
1934 14713 : z = cgetg(N,t_COL);
1935 97807 : for (i=1; i<N; i++) {
1936 83094 : gel(v,i) = gcoeff(q,i,i);
1937 83094 : gel(x,i) = gel(y,i) = gel(z,i) = gen_0;
1938 : }
1939 14713 : if (BORNE)
1940 : {
1941 14692 : borne1 = BORNE;
1942 14692 : if (gsigne(borne1) <= 0) retmkvec3(gen_0, gen_0, cgetg(1,t_MAT));
1943 14678 : if (typ(borne1) != t_REAL)
1944 : {
1945 : long prec;
1946 419 : prec = nbits2prec(gexpo(borne1) + 10);
1947 419 : borne1 = gtofp(borne1, maxss(prec, DEFAULTPREC));
1948 : }
1949 : }
1950 : else
1951 : {
1952 21 : borne1 = gcoeff(q,1,1);
1953 203 : for (i=2; i<N; i++)
1954 : {
1955 182 : GEN b = cholesky_norm_ek(q, i);
1956 182 : if (gcmp(b, borne1) < 0) borne1 = b;
1957 : }
1958 : /* borne1 = norm of smallest basis vector */
1959 : }
1960 14699 : borne2 = mulrr(borne1,alpha);
1961 14699 : if (DEBUGLEVEL>2)
1962 0 : err_printf("smallvectors looking for norm < %P.4G\n",borne1);
1963 14699 : s = 0; k = n;
1964 381909 : for(;; step(x,y,inc,k)) /* main */
1965 : { /* x (supposedly) small vector, ZV.
1966 : * For all t >= k, we have
1967 : * z[t] = sum_{j > t} q[t,j] * x[j]
1968 : * y[t] = sum_{i > t} q[i,i] * (x[i] + z[i])^2
1969 : * = 0 <=> x[i]=0 for all i>t */
1970 : do
1971 : {
1972 1610070 : int skip = 0;
1973 1610070 : if (k > 1)
1974 : {
1975 1228154 : long l = k-1;
1976 1228154 : av1 = avma;
1977 1228154 : p1 = mulimp(gel(x,k), gcoeff(q,l,k));
1978 14458422 : for (j=k+1; j<N; j++) p1 = addmulimp(p1, gel(x,j), gcoeff(q,l,j));
1979 1228154 : gel(z,l) = gc_leaf(av1,p1);
1980 :
1981 1228154 : av1 = avma;
1982 1228154 : p1 = norm_aux(gel(x,k), gel(y,k), gel(z,k), gel(v,k));
1983 1228154 : gel(y,l) = gc_leaf(av1, p1);
1984 : /* skip the [x_1,...,x_skipfirst,0,...,0] */
1985 1228154 : if ((l <= skipfirst && !signe(gel(y,skipfirst)))
1986 1228154 : || mplessthan(borne1, gel(y,l))) skip = 1;
1987 : else /* initial value, minimizing (x[l] + z[l])^2, hence qf(x) for
1988 : the given x[1..l-1] */
1989 1214172 : gel(x,l) = mpround( mpneg(gel(z,l)) );
1990 1228154 : k = l;
1991 : }
1992 1228154 : for(;; step(x,y,inc,k))
1993 : { /* at most 2n loops */
1994 2838224 : if (!skip)
1995 : {
1996 2824242 : if (check_bound(borne1, gel(x,k),gel(y,k),gel(z,k),gel(v,k))) break;
1997 1339945 : step(x,y,inc,k);
1998 1339945 : if (check_bound(borne1, gel(x,k),gel(y,k),gel(z,k),gel(v,k))) break;
1999 : }
2000 1242853 : skip = 0; inc[k] = 1;
2001 1242853 : if (++k > n) goto END;
2002 : }
2003 :
2004 1595371 : if (gc_needed(av,2))
2005 : {
2006 15 : if(DEBUGMEM>1) pari_warn(warnmem,"smallvectors");
2007 15 : if (stockmax) S = clonefill(S, s, stockmax);
2008 15 : if (check) {
2009 15 : GEN dummy = cgetg(1, t_STR);
2010 9629 : for (i=s+1; i<=stockmax; i++) gel(norms,i) = dummy;
2011 : }
2012 15 : (void)gc_all(av,7,&x,&y,&z,&normax1,&borne1,&borne2,&norms);
2013 : }
2014 : }
2015 1595371 : while (k > 1);
2016 381909 : if (!signe(gel(x,1)) && !signe(gel(y,1))) continue; /* exclude 0 */
2017 :
2018 381192 : av1 = avma;
2019 381192 : norme1 = norm_aux(gel(x,1),gel(y,1),gel(z,1),gel(v,1));
2020 381192 : if (mpgreaterthan(norme1,borne1)) { set_avma(av1); continue; /* main */ }
2021 :
2022 381192 : norme1 = gc_leaf(av1,norme1);
2023 381192 : if (check)
2024 : {
2025 312606 : if (checkcnt < 5 && mpcmp(norme1, borne2) < 0)
2026 : {
2027 4420 : if (!check(data,x)) { checkcnt++ ; continue; /* main */}
2028 474 : if (DEBUGLEVEL>4) err_printf("New bound: %Ps", norme1);
2029 474 : borne1 = norme1;
2030 474 : borne2 = mulrr(borne1, alpha);
2031 474 : s = 0; checkcnt = 0;
2032 : }
2033 : }
2034 : else
2035 : {
2036 68586 : if (!BORNE) /* find minimal vectors */
2037 : {
2038 1890 : if (mplessthan(norme1, borne1))
2039 : { /* strictly smaller vector than previously known */
2040 0 : borne1 = norme1; /* + epsilon */
2041 0 : s = 0;
2042 : }
2043 : }
2044 : else
2045 66696 : if (mpcmp(norme1,normax1) > 0) normax1 = norme1;
2046 : }
2047 377246 : if (++s > stockmax) continue; /* too many vectors: no longer remember */
2048 376315 : if (check) gel(norms,s) = norme1;
2049 376315 : gel(S,s) = leafcopy(x);
2050 376315 : if (s != stockmax) continue; /* still room, get next vector */
2051 :
2052 125 : if (check)
2053 : { /* overflow, eliminate vectors failing "check" */
2054 104 : pari_sp av2 = avma;
2055 : long imin, imax;
2056 104 : GEN per = indexsort(norms), S2 = cgetg(stockmax+1, t_VEC);
2057 104 : if (DEBUGLEVEL>2) err_printf("sorting... [%ld elts]\n",s);
2058 : /* let N be the minimal norm so far for x satisfying 'check'. Keep
2059 : * all elements of norm N */
2060 24639 : for (i = 1; i <= s; i++)
2061 : {
2062 24633 : long k = per[i];
2063 24633 : if (check(data,gel(S,k))) { borne1 = gel(norms,k); break; }
2064 : }
2065 104 : imin = i;
2066 21036 : for (; i <= s; i++)
2067 21016 : if (mpgreaterthan(gel(norms,per[i]), borne1)) break;
2068 104 : imax = i;
2069 21036 : for (i=imin, s=0; i < imax; i++) gel(S2,++s) = gel(S,per[i]);
2070 21036 : for (i = 1; i <= s; i++) gel(S,i) = gel(S2,i);
2071 104 : set_avma(av2);
2072 104 : if (s) { borne2 = mulrr(borne1, alpha); checkcnt = 0; }
2073 104 : if (!stockall) continue;
2074 104 : if (s > stockmax/2) stockmax <<= 1;
2075 104 : norms = cgetg(stockmax+1, t_VEC);
2076 21036 : for (i = 1; i <= s; i++) gel(norms,i) = borne1;
2077 : }
2078 : else
2079 : {
2080 21 : if (!stockall && BORNE) goto END;
2081 21 : if (!stockall) continue;
2082 21 : stockmax <<= 1;
2083 : }
2084 :
2085 : {
2086 125 : GEN Snew = clonefill(vec_lengthen(S,stockmax), s, stockmax);
2087 125 : if (isclone(S)) gunclone(S);
2088 125 : S = Snew;
2089 : }
2090 : }
2091 14699 : END:
2092 14699 : if (s < stockmax) stockmax = s;
2093 14699 : if (check)
2094 : {
2095 : GEN per, alph, pols, p;
2096 14671 : if (DEBUGLEVEL>2) err_printf("final sort & check...\n");
2097 14671 : setlg(norms,stockmax+1); per = indexsort(norms);
2098 14671 : alph = cgetg(stockmax+1,t_VEC);
2099 14671 : pols = cgetg(stockmax+1,t_VEC);
2100 84580 : for (j=0,i=1; i<=stockmax; i++)
2101 : {
2102 70171 : long t = per[i];
2103 70171 : GEN N = gel(norms,t);
2104 70171 : if (j && mpgreaterthan(N, borne1)) break;
2105 69909 : if ((p = check(data,gel(S,t))))
2106 : {
2107 55908 : if (!j) borne1 = N;
2108 55908 : j++;
2109 55908 : gel(pols,j) = p;
2110 55908 : gel(alph,j) = gel(S,t);
2111 : }
2112 : }
2113 14671 : setlg(pols,j+1);
2114 14671 : setlg(alph,j+1);
2115 14671 : if (stockmax && isclone(S)) { alph = gcopy(alph); gunclone(S); }
2116 14671 : return mkvec2(pols, alph);
2117 : }
2118 28 : if (stockmax)
2119 : {
2120 21 : setlg(S,stockmax+1);
2121 21 : settyp(S,t_MAT);
2122 21 : if (isclone(S)) { p1 = S; S = gcopy(S); gunclone(p1); }
2123 : }
2124 : else
2125 7 : S = cgetg(1,t_MAT);
2126 28 : return mkvec3(utoi(s<<1), borne1, S);
2127 : }
2128 :
2129 : static GEN
2130 14720 : fincke_pohst_i(GEN r, GEN u, GEN B, long stockmax, long prec, FP_chk_fun *CHECK)
2131 : {
2132 14720 : VOLATILE GEN bound = B, res = NULL;
2133 14720 : pari_CATCH(e_PREC) { }
2134 : pari_TRY {
2135 : GEN q;
2136 14720 : if (CHECK && CHECK->f_init) bound = CHECK->f_init(CHECK, r, u);
2137 14713 : q = gaussred_from_QR(r, prec);
2138 14713 : if (!q) return NULL;
2139 14713 : res = smallvectors(q, bound, stockmax, CHECK);
2140 14713 : } pari_ENDCATCH;
2141 14720 : return res;
2142 : }
2143 : /* solve q(x) = x~.a.x <= bound, a > 0.
2144 : * If check is non-NULL keep x only if check(x).
2145 : * If a is a vector, assume a[1] is the LLL-reduced Cholesky form of q */
2146 : GEN
2147 14734 : fincke_pohst(GEN a, GEN B0, long stockmax, long PREC, FP_chk_fun *CHECK)
2148 : {
2149 14734 : pari_sp av = avma;
2150 : long i, j, l;
2151 : GEN r, rinv, rinvtrans, u, v, res, z, vnorm, rperm, perm, uperm;
2152 :
2153 14734 : if (typ(a) == t_VEC)
2154 : {
2155 14266 : r = gel(a,1);
2156 14266 : u = NULL;
2157 : }
2158 : else
2159 : {
2160 468 : long prec = PREC;
2161 468 : l = lg(a);
2162 468 : if (l == 1)
2163 : {
2164 7 : if (CHECK) pari_err_TYPE("fincke_pohst [dimension 0]", a);
2165 7 : retmkvec3(gen_0, gen_0, cgetg(1,t_MAT));
2166 : }
2167 461 : u = lllfp(a, 0.75, LLL_GRAM | LLL_IM);
2168 454 : if (!u || lg(u) != lg(a)) return gc_NULL(av);
2169 454 : r = qf_RgM_apply(a,u);
2170 454 : i = gprecision(r);
2171 454 : if (i)
2172 412 : prec = i;
2173 : else {
2174 42 : prec = DEFAULTPREC + nbits2extraprec(gexpo(r));
2175 42 : if (prec < PREC) prec = PREC;
2176 : }
2177 454 : if (DEBUGLEVEL>2) err_printf("first LLL: prec = %ld\n", prec);
2178 454 : r = qfgaussred_positive(r);
2179 454 : if (!r) return gc_NULL(av);
2180 1984 : for (i=1; i<l; i++)
2181 : {
2182 1530 : GEN s = gsqrt(gcoeff(r,i,i), prec);
2183 1530 : gcoeff(r,i,i) = s;
2184 4236 : for (j=i+1; j<l; j++) gcoeff(r,i,j) = gmul(s, gcoeff(r,i,j));
2185 : }
2186 : }
2187 : /* now r~ * r = a in LLL basis */
2188 14720 : rinv = RgM_inv_upper(r);
2189 14720 : if (!rinv) return gc_NULL(av);
2190 14720 : rinvtrans = shallowtrans(rinv);
2191 14720 : if (DEBUGLEVEL>2)
2192 0 : err_printf("Fincke-Pohst, final LLL: prec = %ld\n", gprecision(rinvtrans));
2193 14720 : v = lll(rinvtrans);
2194 14720 : if (lg(v) != lg(rinvtrans)) return gc_NULL(av);
2195 :
2196 14720 : rinvtrans = RgM_mul(rinvtrans, v);
2197 14720 : v = ZM_inv(shallowtrans(v),NULL);
2198 14720 : r = RgM_mul(r,v);
2199 14720 : u = u? ZM_mul(u,v): v;
2200 :
2201 14720 : l = lg(r);
2202 14720 : vnorm = cgetg(l,t_VEC);
2203 97842 : for (j=1; j<l; j++) gel(vnorm,j) = gnorml2(gel(rinvtrans,j));
2204 14720 : rperm = cgetg(l,t_MAT);
2205 14720 : uperm = cgetg(l,t_MAT); perm = indexsort(vnorm);
2206 97842 : for (i=1; i<l; i++) { uperm[l-i] = u[perm[i]]; rperm[l-i] = r[perm[i]]; }
2207 14720 : res = fincke_pohst_i(rperm, uperm, B0, stockmax, gprecision(vnorm), CHECK);
2208 14720 : if (!res) return gc_NULL(av);
2209 14713 : if (CHECK)
2210 : {
2211 14671 : if (CHECK->f_post) res = CHECK->f_post(CHECK, res, u);
2212 14671 : return res;
2213 : }
2214 42 : z = cgetg(4,t_VEC);
2215 42 : gel(z,1) = gcopy(gel(res,1));
2216 42 : gel(z,2) = gcopy(gel(res,2));
2217 42 : gel(z,3) = ZM_mul(uperm, gel(res,3)); return gc_upto(av,z);
2218 : }
|
