pari-2.13.1.tar.gz: pari-2.13.1/src/basemath/bibli1.c

"Fossies" - the Fresh Open Source Software Archive

Member "pari-2.13.1/src/basemath/bibli1.c" (12 Dec 2020, 55776 Bytes) of package /linux/misc/pari-2.13.1.tar.gz:

As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "bibli1.c" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report:

2.13.0_vs_2.13.1.

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