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