Annotation of OpenXM_contrib/pari-2.2/src/basemath/bibli1.c, Revision 1.2
1.2 ! noro 1: /* $Id: bibli1.c,v 1.143 2002/09/10 18:40:47 karim Exp $
1.1 noro 2:
3: Copyright (C) 2000 The PARI group.
4:
5: This file is part of the PARI/GP package.
6:
7: PARI/GP is free software; you can redistribute it and/or modify it under the
8: terms of the GNU General Public License as published by the Free Software
9: Foundation. It is distributed in the hope that it will be useful, but WITHOUT
10: ANY WARRANTY WHATSOEVER.
11:
12: Check the License for details. You should have received a copy of it, along
13: with the package; see the file 'COPYING'. If not, write to the Free Software
14: Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */
15:
16: /********************************************************************/
17: /** **/
18: /** LLL Algorithm and close friends **/
19: /** **/
20: /********************************************************************/
21: #include "pari.h"
22: #include "parinf.h"
23: extern GEN ZV_lincomb(GEN u, GEN v, GEN X, GEN Y);
1.2 ! noro 24: extern int addcolumntomatrix(GEN V,GEN INVP,GEN L);
! 25: extern long expodb(double x);
! 26:
! 27: /* default quality ratio for LLL: 99/100 */
! 28: #define LLLDFT 100
1.1 noro 29:
30: /* scalar product x.x */
31: GEN
32: sqscal(GEN x)
33: {
1.2 ! noro 34: long i, lx;
! 35: gpmem_t av;
1.1 noro 36: GEN z;
37: lx = lg(x);
38: if (lx == 1) return gzero;
39: av = avma;
40: z = gsqr((GEN)x[1]);
41: for (i=2; i<lx; i++)
42: z = gadd(z, gsqr((GEN)x[i]));
43: return gerepileupto(av,z);
44: }
45:
46: /* scalar product x.y */
47: GEN
48: gscal(GEN x,GEN y)
49: {
1.2 ! noro 50: long i, lx;
! 51: gpmem_t av;
1.1 noro 52: GEN z;
53: if (x == y) return sqscal(x);
54: lx = lg(x);
55: if (lx == 1) return gzero;
56: av = avma;
57: z = gmul((GEN)x[1],(GEN)y[1]);
58: for (i=2; i<lx; i++)
59: z = gadd(z, gmul((GEN)x[i],(GEN)y[i]));
60: return gerepileupto(av,z);
61: }
62:
63: static GEN
64: sqscali(GEN x)
65: {
1.2 ! noro 66: long i, lx;
! 67: gpmem_t av;
1.1 noro 68: GEN z;
69: lx = lg(x);
70: if (lx == 1) return gzero;
71: av = avma;
72: z = sqri((GEN)x[1]);
73: for (i=2; i<lx; i++)
74: z = addii(z, sqri((GEN)x[i]));
75: return gerepileuptoint(av,z);
76: }
77:
78: static GEN
79: gscali(GEN x,GEN y)
80: {
1.2 ! noro 81: long i, lx;
! 82: gpmem_t av;
1.1 noro 83: GEN z;
84: if (x == y) return sqscali(x);
85: lx = lg(x);
86: if (lx == 1) return gzero;
87: av = avma;
88: z = mulii((GEN)x[1],(GEN)y[1]);
89: for (i=2; i<lx; i++)
90: z = addii(z, mulii((GEN)x[i],(GEN)y[i]));
91: return gerepileuptoint(av,z);
92: }
93:
1.2 ! noro 94: /********************************************************************/
! 95: /** QR Factorization via Householder matrices **/
! 96: /********************************************************************/
! 97:
! 98: /* zero x[1..k-1], fill mu */
! 99: static int
! 100: FindApplyQ(GEN x, GEN mu, GEN B, long k, GEN Q, long prec)
! 101: {
! 102: long i, lx = lg(x)-1, lv = lx - (k-1);
! 103: GEN v, p1, beta, Nx, x2, x1, xd = x + (k-1);
! 104:
! 105: x1 = (GEN)xd[1];
! 106: x2 = gsqr(x1);
! 107: if (k < lx)
! 108: {
! 109: for (i=2; i<=lv; i++) x2 = mpadd(x2, gsqr((GEN)xd[i]));
! 110: Nx = gsqrt(x2, prec);
! 111: if (signe(x1) < 0) setsigne(Nx, -1);
! 112: v = cgetg(lv+1, t_VEC);
! 113: v[1] = lmpadd(x1, Nx);
! 114: for (i=2; i<=lv; i++) v[i] = xd[i];
! 115: if (gcmp0(x2)) return 0;
! 116:
! 117: if (gcmp0(x1))
! 118: beta = mpmul(x2, realun(prec)); /* make sure typ(beta) != t_INT */
! 119: else
! 120: beta = mpadd(x2, mpmul(Nx,x1));
! 121: beta = ginv(beta);
! 122: p1 = cgetg(3,t_VEC); Q[k] = (long)p1;
! 123: p1[1] = (long)beta;
! 124: p1[2] = (long)v;
! 125:
! 126: coeff(mu,k,k) = lmpneg(Nx);
! 127: }
! 128: else
! 129: coeff(mu,k,k) = x[k];
! 130: if (B)
! 131: {
! 132: B[k] = (long)x2;
! 133: for (i=1; i<k; i++) coeff(mu,k,i) = x[i];
! 134: }
! 135: else
! 136: for (i=1; i<k; i++) coeff(mu,i,k) = x[i];
! 137: return (typ(x2) != t_REAL || lg(x2) > 3 || expo(x2) < 10);
! 138: }
! 139:
! 140: static void
! 141: ApplyQ(GEN Q, GEN r)
! 142: {
! 143: GEN s, rd, beta = (GEN)Q[1], v = (GEN)Q[2];
! 144: long i, l = lg(v), lr = lg(r);
! 145:
! 146: rd = r + (lr - l);
! 147: s = mpmul((GEN)v[1], (GEN)rd[1]);
! 148: for (i=2; i<l; i++) s = mpadd(s, mpmul((GEN)v[i] ,(GEN)rd[i]));
! 149: s = mpneg(mpmul(beta, s));
! 150: for (i=1; i<l; i++) rd[i] = lmpadd((GEN)rd[i], mpmul(s, (GEN)v[i]));
! 151: }
! 152:
! 153: static GEN
! 154: ApplyAllQ(GEN Q, GEN r0, long k)
! 155: {
! 156: gpmem_t av = avma;
! 157: GEN r = dummycopy(r0);
! 158: long j;
! 159: for (j=1; j<k; j++) ApplyQ((GEN)Q[j], r);
! 160: return gerepilecopy(av, r);
! 161: }
! 162:
! 163: /* compute B[k] = | x[k] |^2, update mu(k, 1..k-1) using Householder matrices
! 164: * (Q = Householder(x[1..k-1]) in factored form) */
! 165: static int
! 166: incrementalQ(GEN x, GEN L, GEN B, GEN Q, long k, long prec)
! 167: {
! 168: GEN r = ApplyAllQ(Q, (GEN)x[k], k);
! 169: return FindApplyQ(r, L, B, k, Q, prec);
! 170: }
! 171:
! 172: /* Q vector of Householder matrices orthogonalizing x[1..j0].
! 173: * Q[i] = 0 means not computed yet */
! 174: static int
! 175: Householder_get_mu(GEN x, GEN L, GEN B, long k, GEN Q, long prec)
! 176: {
! 177: GEN Nx, invNx, m;
! 178: long i, j, j0;
! 179: if (!Q) Q = zerovec(k);
! 180: for (j=1; j<=k; j++)
! 181: if (typ(Q[j]) == t_INT) break;
! 182: j0 = j;
! 183: for ( ; j<=k; j++)
! 184: if (! incrementalQ(x, L, B, Q, j, prec)) return 0;
! 185: for (j=1; j<=k; j++)
! 186: {
! 187: m = (GEN)L[j]; Nx = (GEN)m[j]; /* should set m[j] = un; but need it later */
! 188: if (j == k) break;
! 189: invNx = ginv(Nx);
! 190: for (i=max(j0, j+1); i<=k; i++) m[i] = lmpmul(invNx, (GEN)m[i]);
! 191: }
! 192: return 1;
! 193: }
! 194:
! 195: GEN
! 196: sqred1_from_QR(GEN x, long prec)
! 197: {
! 198: long j, k = lg(x)-1;
! 199: GEN L, B = zerovec(k);
! 200: L = cgetg(k+1, t_MAT);
! 201: for (j=1; j<=k; j++) L[j] = (long)zerocol(k);
! 202: if (!Householder_get_mu(x, L, B, k, NULL, prec)) return NULL;
! 203: for (j=1; j<=k; j++) coeff(L,j,j) = B[j];
! 204: return gtrans_i(L);
! 205: }
! 206:
! 207: GEN
! 208: R_from_QR(GEN x, long prec)
! 209: {
! 210: long j, k = lg(x)-1;
! 211: GEN L, B = zerovec(k), Q = cgetg(k+1, t_VEC);
! 212: L = cgetg(k+1, t_MAT);
! 213: for (j=1; j<=k; j++) L[j] = (long)zerocol(k);
! 214: for (j=1; j<=k; j++)
! 215: if (!incrementalQ(x, L, B, Q, j, prec)) return NULL;
! 216: return gtrans_i(L);
! 217: }
! 218:
! 219: /********************************************************************/
! 220: /** QR Factorization via Gram-Schmidt **/
! 221: /********************************************************************/
! 222:
! 223: /* compute B[k] = | x[k] |^2, update mu(k, 1..k-1).
! 224: * Classical Gram-Schmidt (unstable!) */
! 225: static int
! 226: incrementalGS(GEN x, GEN mu, GEN B, long k)
! 227: {
! 228: GEN s,A = cgetg(k+1, t_COL); /* scratch space */
! 229: long i, j;
! 230: gpmem_t av;
! 231: A[1] = coeff(x,k,1);
! 232: for(j=1;j<k;)
! 233: {
! 234: coeff(mu,k,j) = lmpdiv((GEN)A[j], (GEN)B[j]);
! 235: j++; av = avma;
! 236: /* A[j] <-- x[k,j] - sum_{i<j} mu[j,i] A[i] */
! 237: s = mpmul(gcoeff(mu,j,1),(GEN)A[1]);
! 238: for (i=2; i<j; i++) s = mpadd(s, mpmul(gcoeff(mu,j,i),(GEN)A[i]));
! 239: s = mpneg(s); A[j] = (long)gerepileuptoleaf(av, mpadd(gcoeff(x,k,j), s));
! 240: }
! 241: B[k] = A[k]; return (signe((GEN)B[k]) > 0);
! 242: }
! 243:
! 244: #if 0
! 245: /* return Gram-Schmidt orthogonal basis (f) associated to (e), B is the
! 246: * vector of the (f_i . f_i)*/
! 247: GEN
! 248: gram_schmidt(GEN e, GEN *ptB)
! 249: {
! 250: long i,j,lx = lg(e);
! 251: GEN f = dummycopy(e), B, iB;
! 252:
! 253: B = cgetg(lx, t_VEC);
! 254: iB= cgetg(lx, t_VEC);
! 255:
! 256: for (i=1;i<lx;i++)
! 257: {
! 258: GEN p1 = NULL;
! 259: gpmem_t av;
! 260: B[i] = (long)sqscal((GEN)f[i]);
! 261: iB[i]= linv((GEN)B[i]); av = avma;
! 262: for (j=1; j<i; j++)
! 263: {
! 264: GEN mu = gmul(gscal((GEN)e[i],(GEN)f[j]), (GEN)iB[j]);
! 265: GEN p2 = gmul(mu, (GEN)f[j]);
! 266: p1 = p1? gadd(p1,p2): p2;
! 267: }
! 268: p1 = p1? gerepileupto(av, gsub((GEN)e[i], p1)): (GEN)e[i];
! 269: f[i] = (long)p1;
! 270: }
! 271: *ptB = B; return f;
! 272: }
! 273: #endif
! 274:
! 275: /********************************************************************/
! 276: /** **/
! 277: /** LLL ALGORITHM **/
! 278: /** **/
! 279: /********************************************************************/
! 280: #define swap(x,y) { long _t=x; x=y; y=_t; }
! 281: #define gswap(x,y) { GEN _t=x; x=y; y=_t; }
! 282:
1.1 noro 283: static GEN
1.2 ! noro 284: lll_trivial(GEN x, long flag)
1.1 noro 285: {
286: GEN y;
1.2 ! noro 287: if (lg(x) == 1)
! 288: { /* dim x = 0 */
! 289: if (! (flag & lll_ALL)) return cgetg(1,t_MAT);
1.1 noro 290: y=cgetg(3,t_VEC);
291: y[1]=lgetg(1,t_MAT);
292: y[2]=lgetg(1,t_MAT); return y;
293: }
1.2 ! noro 294: /* here dim = 1 */
1.1 noro 295: if (gcmp0((GEN)x[1]))
296: {
1.2 ! noro 297: switch(flag & (~lll_GRAM))
1.1 noro 298: {
299: case lll_KER: return idmat(1);
300: case lll_IM : return cgetg(1,t_MAT);
301: default: y=cgetg(3,t_VEC);
302: y[1]=(long)idmat(1);
303: y[2]=lgetg(1,t_MAT); return y;
304: }
305: }
306: if (flag & lll_GRAM) flag ^= lll_GRAM; else x = NULL;
307: switch (flag)
308: {
309: case lll_KER: return cgetg(1,t_MAT);
310: case lll_IM : return idmat(1);
311: default: y=cgetg(3,t_VEC);
312: y[1]=lgetg(1,t_MAT);
313: y[2]=x? lcopy(x): (long)idmat(1); return y;
314: }
315: }
316:
317: static GEN
1.2 ! noro 318: lll_finish(GEN h,GEN fl,long flag)
1.1 noro 319: {
1.2 ! noro 320: long i, k, l = lg(fl);
! 321: GEN y, g;
1.1 noro 322:
1.2 ! noro 323: k=1; while (k<l && !fl[k]) k++;
! 324: switch(flag & (~lll_GRAM))
1.1 noro 325: {
326: case lll_KER: setlg(h,k);
1.2 ! noro 327: return h;
1.1 noro 328:
1.2 ! noro 329: case lll_IM: h += k-1; h[0] = evaltyp(t_MAT) | evallg(l-k+1);
! 330: return h;
1.1 noro 331: }
1.2 ! noro 332: y = cgetg(3,t_VEC);
! 333: g = cgetg(k, t_MAT); for (i=1; i<k; i++) g[i] = h[i];
! 334: y[1] = (long)g;
! 335: h += k-1; h[0] = evaltyp(t_MAT) | evallg(l-k+1);
! 336: y[2] = (long)h; return y;
1.1 noro 337: }
338:
1.2 ! noro 339: /* h[,k] += q * h[,l] */
! 340: static void
! 341: Zupdate_col(long k, long l, GEN q, long K, GEN h)
1.1 noro 342: {
1.2 ! noro 343: GEN *hl, *hk;
! 344: long i;
1.1 noro 345:
1.2 ! noro 346: if (!h) return;
! 347: hl = (GEN*)h[l]; hk = (GEN*)h[k];
! 348: if (is_pm1(q))
1.1 noro 349: {
1.2 ! noro 350: if (signe(q) > 0)
! 351: for (i=1;i<=K;i++) { if (signe(hl[i])) hk[i] = addii(hk[i],hl[i]); }
! 352: else
! 353: for (i=1;i<=K;i++) { if (signe(hl[i])) hk[i] = subii(hk[i],hl[i]); }
! 354: } else
! 355: for (i=1;i<=K;i++) if (signe(hl[i])) hk[i] = addii(hk[i],mulii(q,hl[i]));
! 356: }
1.1 noro 357:
1.2 ! noro 358: /* L[k,] += q * L[l,], l < k */
! 359: static void
! 360: Zupdate_row(long k, long l, GEN q, GEN L, GEN B)
! 361: {
! 362: long i;
! 363: if (is_pm1(q))
! 364: {
! 365: if (signe(q) > 0)
1.1 noro 366: {
1.2 ! noro 367: for (i=1;i<l; i++) coeff(L,k,i) = laddii(gcoeff(L,k,i),gcoeff(L,l,i));
! 368: coeff(L,k,l) = laddii(gcoeff(L,k,l), B);
! 369: } else {
! 370: for (i=1;i<l; i++) coeff(L,k,i) = lsubii(gcoeff(L,k,i),gcoeff(L,l,i));
! 371: coeff(L,k,l) = laddii(gcoeff(L,k,l), negi(B));
1.1 noro 372: }
1.2 ! noro 373: } else {
! 374: for(i=1;i<l;i++) coeff(L,k,i)=laddii(gcoeff(L,k,i),mulii(q,gcoeff(L,l,i)));
! 375: coeff(L,k,l) = laddii(gcoeff(L,k,l), mulii(q,B));
1.1 noro 376: }
377: }
378:
1.2 ! noro 379: static void
! 380: update_row(long k, long l, GEN q, GEN L)
1.1 noro 381: {
1.2 ! noro 382: long i;
! 383: if (is_pm1(q))
1.1 noro 384: {
1.2 ! noro 385: if (signe(q) > 0)
! 386: {
! 387: for (i=1;i<l; i++) coeff(L,k,i) = lmpadd(gcoeff(L,k,i),gcoeff(L,l,i));
! 388: } else {
! 389: for (i=1;i<l; i++) coeff(L,k,i) = lmpsub(gcoeff(L,k,i),gcoeff(L,l,i));
! 390: }
! 391: } else {
! 392: for(i=1;i<l;i++) coeff(L,k,i)=lmpadd(gcoeff(L,k,i),mpmul(q,gcoeff(L,l,i)));
1.1 noro 393: }
1.2 ! noro 394: coeff(L,k,l) = lmpadd(gcoeff(L,k,l), q);
1.1 noro 395: }
396:
397: static void
1.2 ! noro 398: ZRED_gram(long k, long l, GEN x, GEN h, GEN L, GEN B, long K)
1.1 noro 399: {
400: long i,lx;
401: GEN q = truedvmdii(addii(B,shifti(gcoeff(L,k,l),1)), shifti(B,1), NULL);
1.2 ! noro 402: GEN xk,xl;
1.1 noro 403: if (!signe(q)) return;
1.2 ! noro 404: q = negi(q);
! 405: xl = (GEN) x[l]; xk = (GEN) x[k];
! 406: lx = lg(xl);
1.1 noro 407: if (is_pm1(q))
408: {
409: if (signe(q) > 0)
410: {
411: xk[k] = laddii((GEN)xk[k], (GEN)xl[k]);
1.2 ! noro 412: for (i=1;i<lx;i++) coeff(x,k,i) = xk[i] = laddii((GEN)xk[i], (GEN)xl[i]);
1.1 noro 413: } else {
414: xk[k] = lsubii((GEN)xk[k], (GEN)xl[k]);
1.2 ! noro 415: for (i=1;i<lx;i++) coeff(x,k,i) = xk[i] = lsubii((GEN)xk[i], (GEN)xl[i]);
1.1 noro 416: }
1.2 ! noro 417: } else { /* h[,k] += q* h[,l]. x[,k] += q * x[,l]. L[k,] += q* L[l,] */
1.1 noro 418: xk[k] = laddii((GEN)xk[k], mulii(q,(GEN)xl[k]));
419: for(i=1;i<lx;i++) coeff(x,k,i)=xk[i]=laddii((GEN)xk[i],mulii(q,(GEN)xl[i]));
420: }
1.2 ! noro 421: Zupdate_row(k,l,q,L,B);
! 422: Zupdate_col (k,l,q,K,h);
! 423: }
! 424:
! 425: static void
! 426: ZRED(long k, long l, GEN x, GEN h, GEN L, GEN B, long K)
! 427: {
! 428: GEN q = truedvmdii(addii(B,shifti(gcoeff(L,k,l),1)), shifti(B,1), NULL);
! 429: if (!signe(q)) return;
! 430: q = negi(q);
! 431: Zupdate_row(k,l,q,L,B);
! 432: Zupdate_col (k,l,q,K,h);
! 433: x[k] = (long)ZV_lincomb(gun, q, (GEN)x[k], (GEN)x[l]);
! 434: }
! 435:
! 436: static GEN
! 437: round_safe(GEN q)
! 438: {
! 439: long e;
! 440: if (typ(q) == t_INT) return q;
! 441: if (expo(q) > 30)
! 442: {
! 443: q = grndtoi(q, &e);
! 444: if (e > 0) return NULL;
! 445: } else
! 446: q = ground(q);
! 447: return q;
1.1 noro 448: }
449:
450: /* b[k] <-- b[k] - round(L[k,l]) b[l], only b[1] ... b[K] modified so far
451: * assume l < k and update x = Gram(b), L = Gram Schmidt coeffs. */
1.2 ! noro 452: static int
! 453: RED_gram(long k, long l, GEN x, GEN h, GEN L, long K)
! 454: {
! 455: long i,lx;
! 456: GEN q = round_safe(gcoeff(L,k,l));
! 457: GEN xk, xl;
! 458:
! 459: if (!q) return 0;
! 460: if (!signe(q)) return 1;
! 461: q = negi(q); lx = lg(x);
! 462: xk = (GEN)x[k]; xl = (GEN)x[l];
! 463: if (is_pm1(q))
! 464: {
! 465: if (signe(q) > 0)
! 466: {
! 467: xk[k] = lmpadd((GEN)xk[k], (GEN)xl[k]);
! 468: for (i=1;i<lx;i++) coeff(x,k,i)=xk[i]=lmpadd((GEN)xk[i], (GEN)xl[i]);
! 469: } else {
! 470: xk[k] = lmpsub((GEN)xk[k], (GEN)xl[k]);
! 471: for (i=1;i<lx;i++) coeff(x,k,i)=xk[i]=lmpsub((GEN)xk[i], (GEN)xl[i]);
! 472: }
! 473: } else {
! 474: xk[k] = lmpadd((GEN)xk[k], mpmul(q,(GEN)xl[k]));
! 475: for (i=1;i<lx;i++)
! 476: coeff(x,k,i)=xk[i]=lmpadd((GEN)xk[i], mpmul(q,(GEN)xl[i]));
! 477: }
! 478: update_row(k,l,q,L);
! 479: Zupdate_col(k,l,q,K,h); return 1;
! 480: }
! 481:
1.1 noro 482: static void
1.2 ! noro 483: update_col(long k, long l, GEN q, GEN x)
1.1 noro 484: {
1.2 ! noro 485: long i,lx;
! 486: GEN xk = (GEN)x[k], xl = (GEN)x[l];
! 487: lx = lg(xk);
1.1 noro 488: if (is_pm1(q))
489: {
490: if (signe(q) > 0)
491: {
1.2 ! noro 492: for (i=1;i<lx;i++) xk[i] = lmpadd((GEN)xk[i], (GEN)xl[i]);
1.1 noro 493: } else {
1.2 ! noro 494: for (i=1;i<lx;i++) xk[i] = lmpsub((GEN)xk[i], (GEN)xl[i]);
1.1 noro 495: }
496: } else {
1.2 ! noro 497: for (i=1;i<lx;i++) xk[i] = lmpadd((GEN)xk[i], mpmul(q,(GEN)xl[i]));
1.1 noro 498: }
499: }
500:
501: static int
1.2 ! noro 502: RED(long k, long l, GEN x, GEN h, GEN L, long K)
1.1 noro 503: {
1.2 ! noro 504: GEN q = round_safe(gcoeff(L,k,l));
! 505: if (!q) return 0;
! 506: if (!signe(q)) return 1;
! 507: q = negi(q);
! 508: update_col(k,l,q,x);
! 509: update_row(k,l,q,L);
! 510: Zupdate_col(k,l,q,K,h); return 1;
! 511: }
! 512:
! 513: static int
! 514: do_ZSWAP(GEN x, GEN h, GEN L, GEN B, long kmax, long k, long D, GEN fl,
! 515: int gram)
! 516: {
! 517: GEN la,la2,p1,Bk;
! 518: long i, j, lx;
! 519: gpmem_t av;
1.1 noro 520:
521: if (!fl[k-1]) return 0;
1.2 ! noro 522: av = avma;
1.1 noro 523: la = gcoeff(L,k,k-1); la2 = sqri(la);
524: Bk = (GEN)B[k];
1.2 ! noro 525: if (fl[k])
! 526: {
! 527: GEN q;
! 528: if (!D) return 0; /* only swap non-kernel + kernel vector */
! 529: q = addii(la2, mulii((GEN)B[k-1],(GEN)B[k+1]));
! 530: if (cmpii(mulsi(D-1,sqri(Bk)), mulsi(D,q)) <= 0) {
! 531: avma = av; return 0;
! 532: }
! 533: B[k] = (long)diviiexact(q, Bk);
! 534: }
! 535: /* ZSWAP(k-1,k) */
! 536: if (DEBUGLEVEL>3 && k==kmax)
! 537: { /* output diagnostics associated to re-normalized rational quantities */
! 538: gpmem_t av1 = avma;
! 539: GEN d = mulii((GEN)B[k-1],(GEN)B[k+1]);
! 540: p1 = subii(mulsi(D-1, sqri(Bk)), mulsi(D, la2));
! 541: fprintferr(" (%ld)", expi(p1) - expi(mulsi(D, d)));
! 542: avma = av1;
! 543: }
! 544: if (h) swap(h[k-1], h[k]);
! 545: swap(x[k-1], x[k]); lx = lg(x);
! 546: if (gram)
! 547: for (j=1; j < lx; j++) swap(coeff(x,k-1,j), coeff(x,k,j));
1.1 noro 548: for (j=1; j<k-1; j++) swap(coeff(L,k-1,j), coeff(L,k,j))
549: if (fl[k])
550: {
551: av = avma;
552: for (i=k+1; i<=kmax; i++)
553: {
554: GEN t = gcoeff(L,i,k);
555: p1 = subii(mulii((GEN)B[k+1],gcoeff(L,i,k-1)), mulii(la,t));
1.2 ! noro 556: p1 = diviiexact(p1, Bk);
! 557: coeff(L,i,k) = (long)icopy_av(p1,(GEN)av);
! 558: av = avma = (gpmem_t)coeff(L,i,k);
1.1 noro 559:
560: p1 = addii(mulii(la,gcoeff(L,i,k-1)), mulii((GEN)B[k-1],t));
1.2 ! noro 561: p1 = diviiexact(p1, Bk);
! 562: coeff(L,i,k-1) = (long)icopy_av(p1,(GEN)av);
! 563: av = avma = (gpmem_t)coeff(L,i,k-1);
1.1 noro 564: }
1.2 ! noro 565: }
! 566: else if (signe(la))
! 567: {
! 568: p1 = diviiexact(la2, Bk);
! 569: B[k+1] = B[k] = (long)p1;
! 570: for (i=k+2; i<=lx; i++) B[i] = (long)diviiexact(mulii(p1,(GEN)B[i]), Bk);
! 571: for (i=k+1; i<=kmax; i++)
! 572: coeff(L,i,k-1) = (long)diviiexact(mulii(la,gcoeff(L,i,k-1)), Bk);
! 573: for (j=k+1; j<kmax; j++)
! 574: for (i=j+1; i<=kmax; i++)
! 575: coeff(L,i,j) = (long)diviiexact(mulii(p1,gcoeff(L,i,j)), Bk);
1.1 noro 576: }
577: else
578: {
1.2 ! noro 579: for (i=k+1; i<=kmax; i++)
1.1 noro 580: {
1.2 ! noro 581: coeff(L,i,k) = coeff(L,i,k-1);
! 582: coeff(L,i,k-1) = zero;
1.1 noro 583: }
1.2 ! noro 584: B[k] = B[k-1]; fl[k] = 1; fl[k-1] = 0;
1.1 noro 585: }
586: return 1;
587: }
588:
589: static int
1.2 ! noro 590: do_SWAP(GEN x, GEN h, GEN L, GEN B, long kmax, long k, GEN delta, int gram)
1.1 noro 591: {
592: GEN la,la2, BK,BB,q;
1.2 ! noro 593: long i, j, lx;
! 594: gpmem_t av;
1.1 noro 595:
1.2 ! noro 596: av = avma;
1.1 noro 597: la = gcoeff(L,k,k-1); la2 = gsqr(la);
1.2 ! noro 598: q = mpmul((GEN)B[k-1], mpsub(delta,la2));
! 599: if (mpcmp(q, (GEN)B[k]) <= 0) { avma = av; return 0; }
1.1 noro 600:
601: /* SWAP(k-1,k) */
602: if (DEBUGLEVEL>3 && k==kmax) {
603: fprintferr(" (%ld)", gexpo(q) - gexpo((GEN)B[k])); flusherr();
604: }
1.2 ! noro 605: BB = mpadd((GEN)B[k], mpmul((GEN)B[k-1],la2));
! 606: if (!signe(BB)) { B[k] = 0; return 1; } /* precision pb */
1.1 noro 607:
1.2 ! noro 608: coeff(L,k,k-1) = lmpdiv(mpmul(la,(GEN)B[k-1]), BB);
! 609: BK = mpdiv((GEN)B[k], BB);
! 610: B[k] = lmpmul((GEN)B[k-1], BK);
1.1 noro 611: B[k-1] = (long)BB;
612:
613: swap(h[k-1],h[k]);
1.2 ! noro 614: swap(x[k-1],x[k]); lx = lg(x);
! 615: if (gram)
! 616: for (j=1; j < lx; j++) swap(coeff(x,k-1,j), coeff(x,k,j));
1.1 noro 617: for (j=1; j<k-1; j++) swap(coeff(L,k-1,j), coeff(L,k,j))
618: for (i=k+1; i<=kmax; i++)
619: {
620: GEN t = gcoeff(L,i,k);
1.2 ! noro 621: coeff(L,i,k) = lmpsub(gcoeff(L,i,k-1), mpmul(la,t));
! 622: coeff(L,i,k-1) = lmpadd(t, mpmul(gcoeff(L,k,k-1), gcoeff(L,i,k)));
1.1 noro 623: }
624: return 1;
625: }
1.2 ! noro 626: static void
! 627: ZincrementalGS(GEN x, GEN L, GEN B, long k, GEN fl, int gram)
! 628: {
! 629: GEN u = NULL; /* gcc -Wall */
! 630: long i, j, s;
! 631: for (j=1; j<=k; j++)
! 632: if (j==k || fl[j])
! 633: {
! 634: gpmem_t av = avma;
! 635: u = gram? gcoeff(x,k,j): gscali((GEN)x[k], (GEN)x[j]);
! 636: for (i=1; i<j; i++)
! 637: if (fl[i])
! 638: {
! 639: u = subii(mulii((GEN)B[i+1], u), mulii(gcoeff(L,k,i), gcoeff(L,j,i)));
! 640: u = diviiexact(u, (GEN)B[i]);
! 641: }
! 642: coeff(L,k,j) = lpileuptoint(av, u);
! 643: }
! 644: s = signe(u);
! 645: if (s == 0) B[k+1] = B[k];
! 646: else
! 647: {
! 648: if (s < 0) err(lllger3);
! 649: B[k+1] = coeff(L,k,k); coeff(L,k,k) = un; fl[k] = 1;
! 650: }
! 651: }
1.1 noro 652:
1.2 ! noro 653: /* x integer matrix. Beware: this function can return NULL */
1.1 noro 654: GEN
1.2 ! noro 655: lllint_marked(long MARKED, GEN x, long D, int gram,
! 656: GEN *pth, GEN *ptfl, GEN *ptB)
1.1 noro 657: {
1.2 ! noro 658: long lx = lg(x), hx, i, j, k, l, n, kmax;
! 659: gpmem_t av, lim;
! 660: GEN B,L,h,fl;
1.1 noro 661:
1.2 ! noro 662: if (typ(x) != t_MAT) err(typeer,"lllint");
! 663: n = lx-1; if (n <= 1) return NULL;
! 664: hx = lg(x[1]);
! 665: if (gram && hx != lx) err(mattype1,"lllint");
1.1 noro 666: fl = cgetg(lx,t_VECSMALL);
667:
1.2 ! noro 668: av = avma; lim = stack_lim(av,1); x = dummycopy(x);
! 669: B = gscalcol(gun, lx);
! 670: L = cgetg(lx,t_MAT);
1.1 noro 671: for (j=1; j<lx; j++)
672: {
1.2 ! noro 673: for (i=1; i<hx; i++)
! 674: if (typ(gcoeff(x,i,j)) != t_INT) err(lllger4);
1.1 noro 675: fl[j] = 0; L[j] = (long)zerocol(n);
676: }
1.2 ! noro 677: h = pth? idmat(n): NULL;
! 678: ZincrementalGS(x, L, B, 1, fl, gram);
! 679: kmax = 1;
! 680: if (DEBUGLEVEL>5) fprintferr("k = ");
! 681: for (k=2;;)
1.1 noro 682: {
683: if (k > kmax)
684: {
1.2 ! noro 685: if (DEBUGLEVEL>3) fprintferr("K%ld ",k);
! 686: ZincrementalGS(x, L, B, k, fl, gram);
1.1 noro 687: kmax = k;
1.2 ! noro 688: }
! 689: if (k != MARKED)
! 690: {
! 691: if (!gram) ZRED(k,k-1, x,h,L,(GEN)B[k],kmax);
! 692: else ZRED_gram(k,k-1, x,h,L,(GEN)B[k],kmax);
! 693: }
! 694: if (do_ZSWAP(x,h,L,B,kmax,k,D,fl,gram))
1.1 noro 695: {
1.2 ! noro 696: if (MARKED == k) MARKED = k-1;
! 697: else if (MARKED == k-1) MARKED = k;
! 698: if (k > 2) k--;
1.1 noro 699: }
700: else
701: {
1.2 ! noro 702: if (k != MARKED)
! 703: for (l=k-2; l; l--)
1.1 noro 704: {
1.2 ! noro 705: if (!gram) ZRED(k,l, x,h,L,(GEN)B[l+1],kmax);
! 706: else ZRED_gram(k,l, x,h,L,(GEN)B[l+1],kmax);
! 707: if (low_stack(lim, stack_lim(av,1)))
! 708: {
! 709: if(DEBUGMEM>1) err(warnmem,"lllint[1], kmax = %ld", kmax);
! 710: gerepileall(av,h?4:3,&B,&L,&x,&h);
! 711: }
1.1 noro 712: }
713: if (++k > n) break;
714: }
715: if (low_stack(lim, stack_lim(av,1)))
716: {
1.2 ! noro 717: if(DEBUGMEM>1) err(warnmem,"lllint[2], kmax = %ld", kmax);
! 718: gerepileall(av,h?4:3,&B,&L,&x,&h);
1.1 noro 719: }
720: }
721: if (DEBUGLEVEL>3) fprintferr("\n");
1.2 ! noro 722: if (ptB) *ptB = B;
! 723: if (ptfl) *ptfl = fl;
! 724: if (pth) *pth = h;
! 725: if (MARKED && MARKED != 1)
! 726: {
! 727: if (h) { swap(h[1], h[MARKED]); } else swap(x[1], x[MARKED]);
! 728: }
! 729: return h? h: x;
! 730: }
! 731:
! 732: /* Beware: this function can return NULL (dim x <= 1) */
! 733: GEN
! 734: lllint_i(GEN x, long D, int gram, GEN *pth, GEN *ptfl, GEN *ptB)
! 735: {
! 736: return lllint_marked(0, x,D,gram,pth,ptfl,ptB);
1.1 noro 737: }
738:
1.2 ! noro 739: /* return x * lllint(x). No garbage collection */
! 740: GEN
! 741: lllint_ip(GEN x, long D)
! 742: {
! 743: GEN h = lllint_i(x, D, 0, NULL, NULL, NULL);
! 744: if (!h) h = x;
! 745: return h;
! 746: }
! 747:
! 748: GEN
! 749: lllall(GEN x, long D, int gram, long flag)
! 750: {
! 751: gpmem_t av = avma;
! 752: GEN fl, junk, h = lllint_i(x, D, gram, &junk, &fl, NULL);
! 753: if (!h) return lll_trivial(x,flag);
! 754: return gerepilecopy(av, lll_finish(h,fl,flag));
! 755: }
! 756:
! 757: GEN
! 758: lllint(GEN x) { return lllall(x,LLLDFT,0, lll_IM); }
! 759:
! 760: GEN
! 761: lllkerim(GEN x) { return lllall(x,LLLDFT,0, lll_ALL); }
! 762:
! 763: GEN
! 764: lllgramint(GEN x) { return lllall(x,LLLDFT,1, lll_IM | lll_GRAM); }
! 765:
! 766: GEN
! 767: lllgramkerim(GEN x) { return lllall(x,LLLDFT,1, lll_ALL | lll_GRAM); }
1.1 noro 768:
769: static int
770: pslg(GEN x)
771: {
772: long tx;
773: if (gcmp0(x)) return 2;
1.2 ! noro 774: tx = typ(x); return is_scalar_t(tx)? 3: lgef(x);
! 775: }
! 776:
! 777: static GEN
! 778: to_MP(GEN x, long prec)
! 779: {
! 780: GEN y = cgetr(prec); gaffect(x, y); return y;
! 781: }
! 782: static GEN
! 783: col_to_MP(GEN x, long prec)
! 784: {
! 785: long j, l = lg(x);
! 786: GEN y = cgetg(l, t_COL);
! 787: for (j=1; j<l; j++) y[j] = (long)to_MP((GEN)x[j], prec);
! 788: return y;
! 789: }
! 790:
! 791: static GEN
! 792: to_mp(GEN x, long prec)
! 793: {
! 794: GEN y;
! 795: if (typ(x) == t_INT) return x;
! 796: y = cgetr(prec); gaffect(x, y); return y;
! 797: }
! 798: static GEN
! 799: col_to_mp(GEN x, long prec)
! 800: {
! 801: long j, l = lg(x);
! 802: GEN y = cgetg(l, t_COL);
! 803: for (j=1; j<l; j++) y[j] = (long)to_mp((GEN)x[j], prec);
! 804: return y;
! 805: }
! 806: static GEN
! 807: mat_to_mp(GEN x, long prec)
! 808: {
! 809: long j, l = lg(x);
! 810: GEN y = cgetg(l, t_MAT);
! 811: for (j=1; j<l; j++) y[j] = (long)col_to_mp((GEN)x[j], prec);
! 812: return y;
! 813: }
! 814:
! 815: static int
! 816: REDgen(long k, long l, GEN h, GEN L, GEN B)
! 817: {
! 818: GEN u, q, cq;
! 819: long i;
! 820:
! 821: u = gcoeff(L,k,l);
! 822: if (pslg(u) < pslg((GEN)B[l+1])) return 0;
! 823:
! 824: q = gdeuc(u,(GEN)B[l+1]);
! 825: cq = ginv(content(q));
! 826: q = gneg(gmul(q,cq));
! 827: h[k] = ladd(gmul(cq,(GEN)h[k]), gmul(q,(GEN)h[l]));
! 828: coeff(L,k,l) = ladd(gmul(cq,gcoeff(L,k,l)), gmul(q,(GEN)B[l+1]));
! 829: for (i=1; i<l; i++)
! 830: coeff(L,k,i) = ladd(gmul(cq,gcoeff(L,k,i)), gmul(q,gcoeff(L,l,i)));
! 831: return 1;
! 832: }
! 833:
! 834: static int
! 835: do_SWAPgen(GEN h, GEN L, GEN B, long k, GEN fl, int *flc)
! 836: {
! 837: GEN p1, la, la2, Bk;
! 838: long ps1, ps2, i, j, lx;
! 839:
! 840: if (!fl[k-1]) return 0;
! 841:
! 842: la = gcoeff(L,k,k-1); la2 = gsqr(la);
! 843: Bk = (GEN)B[k];
! 844: if (fl[k])
! 845: {
! 846: GEN q = gadd(la2, gmul((GEN)B[k-1],(GEN)B[k+1]));
! 847: ps1 = pslg(gsqr(Bk));
! 848: ps2 = pslg(q);
! 849: if (ps1 <= ps2 && (ps1 < ps2 || !*flc)) return 0;
! 850: *flc = (ps1 != ps2);
! 851: B[k] = ldiv(q, Bk);
! 852: }
! 853:
! 854: swap(h[k-1], h[k]); lx = lg(L);
! 855: for (j=1; j<=k-2; j++) swap(coeff(L,k-1,j), coeff(L,k,j));
! 856: if (fl[k])
! 857: {
! 858: for (i=k+1; i<lx; i++)
! 859: {
! 860: GEN t = gcoeff(L,i,k);
! 861: p1 = gsub(gmul((GEN)B[k+1],gcoeff(L,i,k-1)), gmul(la,t));
! 862: coeff(L,i,k) = ldiv(p1, Bk);
! 863: p1 = gadd(gmul(la,gcoeff(L,i,k-1)), gmul((GEN)B[k-1],t));
! 864: coeff(L,i,k-1) = ldiv(p1, Bk);
! 865: }
! 866: }
! 867: else if (!gcmp0(la))
! 868: {
! 869: p1 = gdiv(la2, Bk);
! 870: B[k+1] = B[k] = (long)p1;
! 871: for (i=k+2; i<=lx; i++) B[i] = ldiv(gmul(p1,(GEN)B[i]),Bk);
! 872: for (i=k+1; i<lx; i++)
! 873: coeff(L,i,k-1) = ldiv(gmul(la,gcoeff(L,i,k-1)), Bk);
! 874: for (j=k+1; j<lx-1; j++)
! 875: for (i=j+1; i<lx; i++)
! 876: coeff(L,i,j) = ldiv(gmul(p1,gcoeff(L,i,j)), Bk);
! 877: }
! 878: else
! 879: {
! 880: coeff(L,k,k-1) = zero;
! 881: for (i=k+1; i<lx; i++)
! 882: {
! 883: coeff(L,i,k) = coeff(L,i,k-1);
! 884: coeff(L,i,k-1) = zero;
! 885: }
! 886: B[k] = B[k-1]; fl[k] = 1; fl[k-1] = 0;
! 887: }
! 888: return 1;
! 889: }
! 890:
! 891: static void
! 892: incrementalGSgen(GEN x, GEN L, GEN B, long k, GEN fl)
! 893: {
! 894: GEN u = NULL; /* gcc -Wall */
! 895: long i, j, tu;
! 896: for (j=1; j<=k; j++)
! 897: if (j==k || fl[j])
! 898: {
! 899: u = gcoeff(x,k,j); tu = typ(u);
! 900: if (! is_scalar_t(tu) && tu != t_POL) err(lllger4);
! 901: for (i=1; i<j; i++)
! 902: if (fl[i])
! 903: {
! 904: u = gsub(gmul((GEN)B[i+1],u), gmul(gcoeff(L,k,i),gcoeff(L,j,i)));
! 905: u = gdiv(u, (GEN)B[i]);
! 906: }
! 907: coeff(L,k,j) = (long)u;
! 908: }
! 909: if (gcmp0(u)) B[k+1] = B[k];
! 910: else
! 911: {
! 912: B[k+1] = coeff(L,k,k); coeff(L,k,k) = un; fl[k] = 1;
! 913: }
1.1 noro 914: }
915:
916: static GEN
917: lllgramallgen(GEN x, long flag)
918: {
1.2 ! noro 919: long lx = lg(x), i, j, k, l, n;
! 920: gpmem_t av0 = avma, av, lim;
! 921: GEN B, L, h, fl;
! 922: int flc;
1.1 noro 923:
924: if (typ(x) != t_MAT) err(typeer,"lllgramallgen");
1.2 ! noro 925: n = lx-1; if (n<=1) return lll_trivial(x,flag);
! 926: if (lg(x[1]) != lx) err(mattype1,"lllgramallgen");
! 927:
! 928: fl = cgetg(lx, t_VECSMALL);
1.1 noro 929:
1.2 ! noro 930: av = avma; lim = stack_lim(av,1);
! 931: B = gscalcol(gun, lx);
! 932: L = cgetg(lx,t_MAT);
! 933: for (j=1; j<lx; j++) { L[j] = (long)zerocol(n); fl[j] = 0; }
1.1 noro 934:
1.2 ! noro 935: h = idmat(n);
1.1 noro 936: for (i=1; i<lx; i++)
1.2 ! noro 937: incrementalGSgen(x, L, B, i, fl);
! 938: flc = 0;
! 939: for(k=2;;)
1.1 noro 940: {
1.2 ! noro 941: if (REDgen(k, k-1, h, L, B)) flc = 1;
! 942: if (do_SWAPgen(h, L, B, k, fl, &flc)) { if (k > 2) k--; }
1.1 noro 943: else
944: {
945: for (l=k-2; l>=1; l--)
1.2 ! noro 946: if (REDgen(k, l, h, L, B)) flc = 1;
1.1 noro 947: if (++k > n) break;
948: }
949: if (low_stack(lim, stack_lim(av,1)))
950: {
951: if(DEBUGMEM>1) err(warnmem,"lllgramallgen");
1.2 ! noro 952: gerepileall(av,3,&B,&L,&h);
1.1 noro 953: }
954: }
1.2 ! noro 955: return gerepilecopy(av0, lll_finish(h,fl,flag));
1.1 noro 956: }
957:
1.2 ! noro 958: static GEN
! 959: _mul2n(GEN x, long e) { return e? gmul2n(x, e): x; }
! 960:
! 961: /* E = scaling exponents
! 962: * F = backward scaling exponents
! 963: * X = lattice basis, Xs = associated scaled basis
! 964: * Q = vector of Householder transformations
! 965: * h = base change matrix (from X0)
! 966: * R = from QR factorization of Xs[1..k-1] */
1.1 noro 967: static int
1.2 ! noro 968: HRS(int MARKED, long k, int prim, long kmax, GEN X, GEN Xs, GEN h, GEN R,
! 969: GEN Q, GEN E, GEN F)
! 970: {
! 971: long e, i, N = lg(X[k]);
! 972: const long prec = MEDDEFAULTPREC; /* 128 bits */
! 973: GEN q, tau2, rk;
! 974: int rounds = 0, overf;
! 975:
! 976: E[k] = prim? E[k-1]: 0;
! 977: F[k] = 0;
! 978: Xs[k] = E[k]? lmul2n((GEN)X[k], E[k]): (long)dummycopy((GEN)X[k]);
! 979: rounds = 0;
! 980: if (k == MARKED) goto DONE; /* no size reduction/scaling */
! 981:
! 982: UP:
! 983: rk = ApplyAllQ(Q, col_to_MP((GEN)Xs[k], prec), k);
! 984: overf = 0;
! 985: for (i = k-1; i > 0; i--)
! 986: { /* size seduction of Xs[k] */
! 987: GEN q2, mu, muint;
! 988: long ex;
! 989: gpmem_t av = avma;
! 990:
! 991: mu = mpdiv((GEN)rk[i], gcoeff(R,i,i));
! 992: if (!signe(mu)) { avma = av; continue; }
! 993: mu = _mul2n(mu, -F[i]); /*backward scaling*/
! 994: ex = gexpo(mu);
! 995: if (ex > 10) {
! 996: overf = 1; /* large size reduction */
! 997: muint = ex > 30? ceil_safe(mu): ground(mu);
! 998: }
! 999: else if (fabs(gtodouble(mu)) > 0.51)
! 1000: muint = ground(mu);
! 1001: else { avma = av; continue; }
! 1002:
! 1003: q = gerepileuptoint(av, negi(muint));
! 1004: q2= _mul2n(q, F[i]); /*backward scaling*/
! 1005: Zupdate_col(k, i, q2, N-1, Xs);
! 1006: rk = gadd(rk, gmul(q2, (GEN)R[i]));
! 1007: /* if in HRS', propagate the last SR on X (for SWAP test) */
! 1008: if (prim && (i == k-1)) {
! 1009: Zupdate_col(k, i, q, kmax, h);
! 1010: update_col(k, i, q, X);
! 1011: }
! 1012: }
! 1013: /* If large size-reduction performed, restart from exact data */
! 1014: if (overf)
! 1015: {
! 1016: if (++rounds > 100) return 0; /* detect infinite loop */
! 1017: goto UP;
! 1018: }
! 1019:
! 1020: if (prim || k == 1) goto DONE; /* DONE */
! 1021:
! 1022: rounds = 0;
! 1023: /* rescale Xs[k] ? */
! 1024: tau2 = signe(rk[k])? gsqr((GEN)rk[k]): gzero;
! 1025: for (i = k+1; i < N; i++)
! 1026: if (signe(rk[i])) tau2 = mpadd(tau2, gsqr((GEN)rk[i]));
! 1027: q = mpdiv(gsqr(gcoeff(R,1,1)), tau2);
! 1028: e = gexpo(q)/2 + F[1]; /* ~ log_2 (||Xs[1]|| / tau), backward scaling*/
! 1029: if (e > 0)
! 1030: { /* rescale */
! 1031: if (e > 30) e = 30;
! 1032: Xs[k] = lmul2n((GEN)Xs[k], e);
! 1033: E[k] += e; goto UP; /* one more round */
! 1034: }
! 1035: else if (e < 0)
! 1036: { /* enable backward scaling */
! 1037: for (i = 1; i < k; i++) F[i] -= e;
! 1038: }
! 1039:
! 1040: DONE:
! 1041: rk = ApplyAllQ(Q, col_to_MP((GEN)Xs[k], prec), k);
! 1042: (void)FindApplyQ(rk, R, NULL, k, Q, prec);
! 1043: return 1;
! 1044: }
! 1045:
! 1046: static GEN
! 1047: rescale_to_int(GEN x)
1.1 noro 1048: {
1.2 ! noro 1049: long e, prec = gprecision(x);
! 1050: GEN y;
! 1051: if (!prec) return Q_primpart(x);
! 1052:
! 1053: y = gmul2n(x, bit_accuracy(prec) - gexpo(x));
! 1054: return gcvtoi(y, &e);
! 1055: }
! 1056:
! 1057: GEN
! 1058: lll_scaled(int MARKED, GEN X0, long D)
! 1059: {
! 1060: GEN delta, X, Xs, h, R, Q, E, F;
! 1061: long j, kmax = 1, k, N = lg(X0);
! 1062: gpmem_t lim, av, av0 = avma;
! 1063: int retry = 0;
! 1064:
! 1065: delta = stor(D-1, DEFAULTPREC);
! 1066: delta = divrs(delta,D);
! 1067: E = vecsmall_const(N-1, 0);
! 1068: F = vecsmall_const(N-1, 0);
! 1069:
! 1070: av = avma; lim = stack_lim(av, 1);
! 1071: h = idmat(N-1);
! 1072: X = rescale_to_int(X0);
! 1073:
! 1074: PRECPB:
! 1075: k = 1;
! 1076: if (retry++) return _vec(h);
! 1077: Q = zerovec(N-1);
! 1078: Xs = zeromat(N-1, N-1);
! 1079: R = cgetg(N, t_MAT);
! 1080: for (j=1; j<N; j++) R[j] = (long)zerocol(N-1);
! 1081:
! 1082: while (k < N)
! 1083: {
! 1084: gpmem_t av1;
! 1085: if (k == 1) {
! 1086: HRS(MARKED, 1, 0, kmax, X, Xs, h, R, Q, E, F);
! 1087: k = 2;
! 1088: }
! 1089: if (k > kmax) {
! 1090: kmax = k;
! 1091: if (DEBUGLEVEL>3) {fprintferr(" K%ld",k);flusherr();}
! 1092: }
! 1093: if (!HRS(MARKED, k, 1, kmax, X, Xs, h, R, Q, E, F)) goto PRECPB;
! 1094:
! 1095: av1 = avma;
! 1096: if (mpcmp( mpmul(delta, gsqr(gcoeff(R, k-1,k-1))),
! 1097: mpadd(gsqr(gcoeff(R,k-1,k)), gsqr(gcoeff(R, k,k)))) > 0)
! 1098: {
! 1099: if (DEBUGLEVEL>3 && k == kmax)
! 1100: {
! 1101: GEN q = mpsub( mpmul(delta, gsqr(gcoeff(R, k-1,k-1))),
! 1102: gsqr(gcoeff(R,k-1,k)));
! 1103: fprintferr(" (%ld)", gexpo(q) - gexpo(gsqr(gcoeff(R, k,k))));
! 1104: }
! 1105: swap(X[k], X[k-1]);
! 1106: swap(h[k], h[k-1]);
! 1107: if (MARKED == k) MARKED = k-1;
! 1108: else if (MARKED == k-1) MARKED = k;
! 1109: avma = av1; k--;
! 1110: }
! 1111: else
! 1112: {
! 1113: avma = av1;
! 1114: if (!HRS(MARKED, k, 0, kmax, X, Xs, h, R, Q, E, F)) goto PRECPB;
! 1115: k++;
! 1116: }
! 1117: if (low_stack(lim, stack_lim(av,1)))
! 1118: {
! 1119: if(DEBUGMEM>1) err(warnmem,"lllfp[1]");
! 1120: gerepileall(av,5,&X,&Xs, &R,&h,&Q);
! 1121: }
1.1 noro 1122: }
1.2 ! noro 1123: return gerepilecopy(av0, h);
1.1 noro 1124: }
1125:
1126: /* x = Gram(b_i). If precision problems return NULL if flag=1, error otherwise.
1.2 ! noro 1127: * Quality ratio = delta = (D-1)/D. Suggested value: D = 100
! 1128: *
! 1129: * If MARKED != 0 make sure e[MARKED] is the first vector of the output basis
! 1130: * (which may then not be LLL-reduced) */
1.1 noro 1131: GEN
1.2 ! noro 1132: lllfp_marked(int MARKED, GEN x, long D, long flag, long prec, int gram)
1.1 noro 1133: {
1.2 ! noro 1134: GEN xinit,L,h,B,L1,delta, Q;
! 1135: long retry = 2, lx = lg(x), hx, l, i, j, k, k1, n, kmax, KMAX;
! 1136: gpmem_t av0 = avma, av, lim;
! 1137:
! 1138: if (typ(x) != t_MAT) err(typeer,"lllfp");
! 1139: n = lx-1; if (n <= 1) return idmat(n);
! 1140: #if 0 /* doesn't work yet */
! 1141: return lll_scaled(MARKED, x, D);
! 1142: #endif
! 1143:
! 1144: hx = lg(x[1]);
! 1145: if (gram && hx != lx) err(mattype1,"lllfp");
! 1146: delta = stor(D-1, DEFAULTPREC);
! 1147: delta = divrs(delta,D);
! 1148: xinit = x;
! 1149: av = avma; lim = stack_lim(av,1);
! 1150: if (gram) {
! 1151: for (k=2,j=1; j<lx; j++)
! 1152: {
! 1153: GEN p1=(GEN)x[j];
! 1154: for (i=1; i<=j; i++)
! 1155: if (typ(p1[i]) == t_REAL) { l = lg(p1[i]); if (l>k) k=l; }
! 1156: }
! 1157: } else {
! 1158: for (k=2,j=1; j<lx; j++)
! 1159: {
! 1160: GEN p1=(GEN)x[j];
! 1161: for (i=1; i<hx; i++)
! 1162: if (typ(p1[i]) == t_REAL) { l = lg(p1[i]); if (l>k) k=l; }
! 1163: }
1.1 noro 1164: }
1165: if (k == 2)
1166: {
1.2 ! noro 1167: if (!prec) return gram? lllgramint(x): lllint(x);
1.1 noro 1168: x = gmul(x, realun(prec+1));
1169: }
1.2 ! noro 1170: else
1.1 noro 1171: {
1172: if (prec < k) prec = k;
1.2 ! noro 1173: x = mat_to_mp(x, prec+1);
1.1 noro 1174: }
1.2 ! noro 1175: /* kmax = maximum column index attained during this run
1.1 noro 1176: * KMAX = same over all runs (after PRECPB) */
1177: kmax = KMAX = 1;
1178:
1.2 ! noro 1179: Q = zerovec(n);
! 1180:
1.1 noro 1181: PRECPB:
1182: switch(retry--)
1183: {
1184: case 2: h = idmat(n); break; /* entry */
1185: case 1:
1.2 ! noro 1186: if (DEBUGLEVEL > 3) fprintferr("\n");
! 1187: if (flag == 2) return _vec(h);
! 1188: if (gram && kmax > 2)
1.1 noro 1189: { /* some progress but precision loss, try again */
1190: prec = (prec<<1)-2; kmax = 1;
1.2 ! noro 1191: if (DEBUGLEVEL) err(warnprec,"lllfp",prec);
! 1192: x = gprec_w(xinit,prec);
! 1193: if (gram) x = qf_base_change(x,h,1); else x = gmul(x,h);
! 1194: gerepileall(av,2,&h,&x);
! 1195: if (DEBUGLEVEL) err(warner,"lllfp starting over");
1.1 noro 1196: break;
1197: } /* fall through */
1198: case 0: /* give up */
1199: if (DEBUGLEVEL > 3) fprintferr("\n");
1.2 ! noro 1200: if (DEBUGLEVEL) err(warner,"lllfp giving up");
1.1 noro 1201: if (flag) { avma=av; return NULL; }
1202: err(lllger3);
1203: }
1204:
1205: L=cgetg(lx,t_MAT);
1206: B=cgetg(lx,t_COL);
1207: for (j=1; j<lx; j++) { L[j] = (long)zerocol(n); B[j] = zero; }
1.2 ! noro 1208: if (gram && !incrementalGS(x, L, B, 1))
1.1 noro 1209: {
1.2 ! noro 1210: if (!flag) return NULL;
1.1 noro 1211: err(lllger3);
1212: }
1213: if (DEBUGLEVEL>5) fprintferr("k =");
1.2 ! noro 1214: for(k=2;;)
1.1 noro 1215: {
1.2 ! noro 1216: if (k > kmax)
1.1 noro 1217: {
1218: kmax = k; if (KMAX < kmax) KMAX = kmax;
1219: if (DEBUGLEVEL>3) {fprintferr(" K%ld",k);flusherr();}
1.2 ! noro 1220: if (gram) j = incrementalGS(x, L, B, k);
! 1221: else j = Householder_get_mu(x, L, B, k, Q, prec);
! 1222: if (!j) goto PRECPB;
1.1 noro 1223: }
1224: else if (DEBUGLEVEL>5) fprintferr(" %ld",k);
1225: L1 = gcoeff(L,k,k-1);
1.2 ! noro 1226: if (typ(L1) == t_REAL && expo(L1) + 20 > bit_accuracy(lg(L1)))
1.1 noro 1227: {
1228: if (DEBUGLEVEL>3)
1.2 ! noro 1229: fprintferr("\nRecomputing Gram-Schmidt, kmax = %ld\n", kmax);
! 1230: if (!gram) goto PRECPB;
1.1 noro 1231: for (k1=1; k1<=kmax; k1++)
1.2 ! noro 1232: if (!incrementalGS(x, L, B, k1)) goto PRECPB;
! 1233: }
! 1234: if (k != MARKED)
! 1235: {
! 1236: if (!gram) j = RED(k,k-1, x,h,L,KMAX);
! 1237: else j = RED_gram(k,k-1, x,h,L,KMAX);
! 1238: if (!j) goto PRECPB;
1.1 noro 1239: }
1.2 ! noro 1240: if (do_SWAP(x,h,L,B,kmax,k,delta,gram))
1.1 noro 1241: {
1.2 ! noro 1242: if (MARKED == k) MARKED = k-1;
! 1243: else if (MARKED == k-1) MARKED = k;
1.1 noro 1244: if (!B[k]) goto PRECPB;
1.2 ! noro 1245: Q[k] = Q[k-1] = zero;
1.1 noro 1246: if (k>2) k--;
1247: }
1248: else
1249: {
1.2 ! noro 1250: if (k != MARKED)
! 1251: for (l=k-2; l; l--)
! 1252: {
! 1253: if (!gram) j = RED(k,l, x,h,L,KMAX);
! 1254: else j = RED_gram(k,l, x,h,L,KMAX);
! 1255: if (!j) goto PRECPB;
! 1256: if (low_stack(lim, stack_lim(av,1)))
! 1257: {
! 1258: if(DEBUGMEM>1) err(warnmem,"lllfp[1], kmax = %ld", kmax);
! 1259: gerepileall(av,5,&B,&L,&h,&x,&Q);
! 1260: }
! 1261: }
! 1262: if (++k > n)
! 1263: {
! 1264: if (!gram && Q[n-1] == zero)
! 1265: {
! 1266: if (DEBUGLEVEL>3) fprintferr("\nChecking LLL basis\n");
! 1267: j = Householder_get_mu(gmul(xinit,h), L, B, n, Q, prec);
! 1268: if (!j) goto PRECPB;
! 1269: k = 2; continue;
! 1270: }
! 1271: break;
! 1272: }
1.1 noro 1273: }
1274: if (low_stack(lim, stack_lim(av,1)))
1275: {
1.2 ! noro 1276: if(DEBUGMEM>1) err(warnmem,"lllfp[2], kmax = %ld", kmax);
! 1277: gerepileall(av,5,&B,&L,&h,&x,&Q);
1.1 noro 1278: }
1279: }
1280: if (DEBUGLEVEL>3) fprintferr("\n");
1.2 ! noro 1281: if (MARKED && MARKED != 1) swap(h[1], h[MARKED]);
! 1282: return gerepilecopy(av0, h);
1.1 noro 1283: }
1284:
1.2 ! noro 1285: GEN
! 1286: lllgramintern(GEN x, long D, long flag, long prec)
! 1287: {
! 1288: return lllfp_marked(0, x,D,flag,prec,1);
! 1289: }
! 1290:
! 1291: GEN
! 1292: lllintern(GEN x, long D, long flag, long prec)
! 1293: {
! 1294: return lllfp_marked(0, x,D,flag,prec,0);
! 1295: }
! 1296:
! 1297: GEN
! 1298: lllgram(GEN x,long prec) { return lllgramintern(x,LLLDFT,0,prec); }
1.1 noro 1299:
1300: GEN
1.2 ! noro 1301: lll(GEN x,long prec) { return lllintern(x,LLLDFT,0,prec); }
1.1 noro 1302:
1303: GEN
1304: qflll0(GEN x, long flag, long prec)
1305: {
1306: switch(flag)
1307: {
1308: case 0: return lll(x,prec);
1309: case 1: return lllint(x);
1310: case 2: return lllintpartial(x);
1311: case 4: return lllkerim(x);
1312: case 5: return lllkerimgen(x);
1313: case 8: return lllgen(x);
1314: default: err(flagerr,"qflll");
1315: }
1316: return NULL; /* not reached */
1317: }
1318:
1319: GEN
1320: qflllgram0(GEN x, long flag, long prec)
1321: {
1322: switch(flag)
1323: {
1324: case 0: return lllgram(x,prec);
1325: case 1: return lllgramint(x);
1326: case 4: return lllgramkerim(x);
1327: case 5: return lllgramkerimgen(x);
1328: case 8: return lllgramgen(x);
1329: default: err(flagerr,"qflllgram");
1330: }
1331: return NULL; /* not reached */
1332: }
1333:
1334: GEN
1.2 ! noro 1335: gram_matrix(GEN b)
1.1 noro 1336: {
1.2 ! noro 1337: long i,j, lx = lg(b);
! 1338: GEN g;
! 1339: if (typ(b) != t_MAT) err(typeer,"gram");
! 1340: g = cgetg(lx,t_MAT);
! 1341: for (i=1; i<lx; i++)
1.1 noro 1342: {
1.2 ! noro 1343: g[i] = lgetg(lx,t_COL);
! 1344: for (j=1; j<=i; j++)
! 1345: coeff(g,i,j) = coeff(g,j,i) = (long)gscal((GEN)b[i],(GEN)b[j]);
1.1 noro 1346: }
1.2 ! noro 1347: return g;
! 1348: }
1.1 noro 1349:
1.2 ! noro 1350: GEN
! 1351: lllgen(GEN x) {
! 1352: gpmem_t av = avma;
! 1353: return gerepileupto(av, lllgramgen(gram_matrix(x)));
1.1 noro 1354: }
1355:
1356: GEN
1.2 ! noro 1357: lllkerimgen(GEN x) {
! 1358: gpmem_t av = avma;
! 1359: return gerepileupto(av, lllgramallgen(gram_matrix(x),lll_ALL));
1.1 noro 1360: }
1361:
1362: GEN
1.2 ! noro 1363: lllgramgen(GEN x) { return lllgramallgen(x,lll_IM); }
! 1364:
! 1365: GEN
! 1366: lllgramkerimgen(GEN x) { return lllgramallgen(x,lll_ALL); }
1.1 noro 1367:
1368: /* This routine is functionally similar to lllint when all = 0.
1369: *
1370: * The input is an integer matrix mat (not necessarily square) whose
1371: * columns are linearly independent. It outputs another matrix T such that
1372: * mat * T is partially reduced. If all = 0, the size of mat * T is the
1373: * same as the size of mat. If all = 1 the number of columns of mat * T
1374: * is at most equal to its number of rows. A matrix M is said to be
1375: * -partially reduced- if | m1 +- m2 | >= |m1| for any two distinct
1376: * columns m1, m2, in M.
1377: *
1378: * This routine is designed to quickly reduce lattices in which one row
1379: * is huge compared to the other rows. For example, when searching for a
1380: * polynomial of degree 3 with root a mod p, the four input vectors might
1381: * be the coefficients of
1382: * X^3 - (a^3 mod p), X^2 - (a^2 mod p), X - (a mod p), p.
1383: * All four constant coefficients are O(p) and the rest are O(1). By the
1384: * pigeon-hole principle, the coefficients of the smallest vector in the
1385: * lattice are O(p^(1/4)), Hence significant reduction of vector lengths
1386: * can be anticipated.
1387: *
1388: * Peter Montgomery (July, 1994)
1389: *
1390: * If flag = 1 complete the reduction using lllint, otherwise return
1391: * partially reduced basis.
1392: */
1393: GEN
1394: lllintpartialall(GEN m, long flag)
1395: {
1396: const long ncol = lg(m)-1;
1.2 ! noro 1397: const gpmem_t ltop1 = avma;
1.1 noro 1398: long ncolnz;
1.2 ! noro 1399: GEN tm1, tm2, mid;
1.1 noro 1400:
1401: if (typ(m) != t_MAT) err(typeer,"lllintpartialall");
1402: if (ncol <= 1) return idmat(ncol);
1403:
1404: {
1405: GEN dot11 = sqscali((GEN)m[1]);
1406: GEN dot22 = sqscali((GEN)m[2]);
1407: GEN dot12 = gscali((GEN)m[1], (GEN)m[2]);
1408: GEN tm = idmat(2); /* For first two columns only */
1409:
1410: int progress = 0;
1411: long npass2 = 0;
1412:
1413: /* Try to row reduce the first two columns of m.
1414: * Our best result so far is (first two columns of m)*tm.
1415: *
1416: * Initially tm = 2 x 2 identity matrix.
1417: * The inner products of the reduced matrix are in
1418: * dot11, dot12, dot22.
1419: */
1420: while (npass2 < 2 || progress)
1421: {
1.2 ! noro 1422: GEN dot12new,q = diviiround(dot12, dot22);
1.1 noro 1423:
1424: npass2++; progress = signe(q);
1425: if (progress)
1426: {
1427: /* Conceptually replace (v1, v2) by (v1 - q*v2, v2),
1428: * where v1 and v2 represent the reduced basis for the
1429: * first two columns of the matrix.
1430: *
1431: * We do this by updating tm and the inner products.
1432: *
1433: * An improved algorithm would look only at the leading
1434: * digits of dot11, dot12, dot22. It would use
1435: * single-precision calculations as much as possible.
1436: */
1437: q = negi(q);
1438: dot12new = addii(dot12, mulii(q, dot22));
1439: dot11 = addii(dot11, mulii(q, addii(dot12, dot12new)));
1440: dot12 = dot12new;
1441: tm[1] = (long)ZV_lincomb(gun,q, (GEN)tm[1],(GEN)tm[2]);
1442: }
1443:
1444: /* Interchange the output vectors v1 and v2. */
1445: gswap(dot11,dot22); swap(tm[1],tm[2]);
1446:
1447: /* Occasionally (including final pass) do garbage collection. */
1448: if (npass2 % 8 == 0 || !progress)
1.2 ! noro 1449: gerepileall(ltop1, 4, &dot11,&dot12,&dot22,&tm);
1.1 noro 1450: } /* while npass2 < 2 || progress */
1451:
1452: {
1453: long icol;
1454: GEN det12 = subii(mulii(dot11, dot22), mulii(dot12, dot12));
1455:
1456: tm1 = idmat(ncol);
1457: mid = cgetg(ncol+1, t_MAT);
1458: for (icol = 1; icol <= 2; icol++)
1459: {
1460: coeff(tm1,1,icol) = coeff(tm,1,icol);
1461: coeff(tm1,2,icol) = coeff(tm,2,icol);
1462: mid[icol] = (long)ZV_lincomb(
1463: gcoeff(tm,1,icol),gcoeff(tm,2,icol), (GEN)m[1],(GEN)m[2]);
1464: }
1465: for (icol = 3; icol <= ncol; icol++)
1466: {
1467: GEN curcol = (GEN)m[icol];
1468: GEN dot1i = gscali((GEN)mid[1], curcol);
1469: GEN dot2i = gscali((GEN)mid[2], curcol);
1470: /* Try to solve
1471: *
1472: * ( dot11 dot12 ) (q1) ( dot1i )
1473: * ( dot12 dot22 ) (q2) = ( dot2i )
1474: *
1475: * Round -q1 and -q2 to the nearest integer.
1476: * Then compute curcol - q1*mid[1] - q2*mid[2].
1477: * This will be approximately orthogonal to the
1478: * first two vectors in the new basis.
1479: */
1480: GEN q1neg = subii(mulii(dot12, dot2i), mulii(dot22, dot1i));
1481: GEN q2neg = subii(mulii(dot12, dot1i), mulii(dot11, dot2i));
1482:
1.2 ! noro 1483: q1neg = diviiround(q1neg, det12);
! 1484: q2neg = diviiround(q2neg, det12);
1.1 noro 1485: coeff(tm1, 1, icol) = ladd(gmul(q1neg, gcoeff(tm,1,1)),
1486: gmul(q2neg, gcoeff(tm,1,2)));
1487: coeff(tm1, 2, icol) = ladd(gmul(q1neg, gcoeff(tm,2,1)),
1488: gmul(q2neg, gcoeff(tm,2,2)));
1489: mid[icol] = ladd(curcol,
1490: ZV_lincomb(q1neg,q2neg, (GEN)mid[1],(GEN)mid[2]));
1491: } /* for icol */
1492: } /* local block */
1493: }
1494: if (DEBUGLEVEL>6)
1495: {
1496: fprintferr("tm1 = %Z",tm1);
1497: fprintferr("mid = %Z",mid);
1498: }
1.2 ! noro 1499: gerepileall(ltop1,2, &tm1, &mid);
1.1 noro 1500: {
1501: /* For each pair of column vectors v and w in mid * tm2,
1502: * try to replace (v, w) by (v, v - q*w) for some q.
1503: * We compute all inner products and check them repeatedly.
1504: */
1.2 ! noro 1505: const gpmem_t ltop3 = avma; /* Excludes region with tm1 and mid */
! 1506: const gpmem_t lim = stack_lim(ltop3,1);
! 1507: long icol, reductions, npass = 0;
1.1 noro 1508: GEN dotprd = cgetg(ncol+1, t_MAT);
1509:
1510: tm2 = idmat(ncol);
1511: for (icol=1; icol <= ncol; icol++)
1512: {
1513: long jcol;
1514:
1515: dotprd[icol] = lgetg(ncol+1,t_COL);
1516: for (jcol=1; jcol <= icol; jcol++)
1517: coeff(dotprd,jcol,icol) = coeff(dotprd,icol,jcol) =
1518: (long)gscali((GEN)mid[icol], (GEN)mid[jcol]);
1519: } /* for icol */
1520: for(;;)
1521: {
1522: reductions = 0;
1523: for (icol=1; icol <= ncol; icol++)
1524: {
1525: long ijdif, jcol, k1, k2;
1526: GEN codi, q;
1527:
1528: for (ijdif=1; ijdif < ncol; ijdif++)
1529: {
1.2 ! noro 1530: const gpmem_t previous_avma = avma;
1.1 noro 1531:
1532: jcol = (icol + ijdif - 1) % ncol; jcol++; /* Hack for NeXTgcc 2.5.8 */
1533: k1 = (cmpii(gcoeff(dotprd,icol,icol),
1534: gcoeff(dotprd,jcol,jcol) ) > 0)
1535: ? icol : jcol; /* index of larger column */
1536: k2 = icol + jcol - k1; /* index of smaller column */
1537: codi = gcoeff(dotprd,k2,k2);
1.2 ! noro 1538: q = gcmp0(codi)? gzero: diviiround(gcoeff(dotprd,k1,k2), codi);
1.1 noro 1539:
1540: /* Try to subtract a multiple of column k2 from column k1. */
1541: if (gcmp0(q)) avma = previous_avma;
1542: else
1543: {
1544: long dcol;
1545:
1546: reductions++; q = negi(q);
1547: tm2[k1]=(long)
1548: ZV_lincomb(gun,q, (GEN)tm2[k1], (GEN)tm2[k2]);
1549: dotprd[k1]=(long)
1550: ZV_lincomb(gun,q, (GEN)dotprd[k1], (GEN)dotprd[k2]);
1551: coeff(dotprd, k1, k1) = laddii(gcoeff(dotprd,k1,k1),
1552: mulii(q, gcoeff(dotprd,k2,k1)));
1553: for (dcol = 1; dcol <= ncol; dcol++)
1554: coeff(dotprd,k1,dcol) = coeff(dotprd,dcol,k1);
1555: } /* if q != 0 */
1556: } /* for ijdif */
1557: if (low_stack(lim, stack_lim(ltop3,1)))
1558: {
1559: if(DEBUGMEM>1) err(warnmem,"lllintpartialall");
1.2 ! noro 1560: gerepileall(ltop3, 2, &dotprd,&tm2);
1.1 noro 1561: }
1562: } /* for icol */
1563: if (!reductions) break;
1564: if (DEBUGLEVEL>6)
1565: {
1566: GEN diag_prod = dbltor(1.0);
1567: for (icol = 1; icol <= ncol; icol++)
1568: diag_prod = gmul(diag_prod, gcoeff(dotprd,icol,icol));
1569: npass++;
1570: fprintferr("npass = %ld, red. last time = %ld, diag_prod = %Z\n\n",
1571: npass, reductions, diag_prod);
1572: }
1573: } /* for(;;)*/
1574:
1575: /* Sort columns so smallest comes first in m * tm1 * tm2.
1576: * Use insertion sort.
1577: */
1578: for (icol = 1; icol < ncol; icol++)
1579: {
1580: long jcol, s = icol;
1581:
1582: for (jcol = icol+1; jcol <= ncol; jcol++)
1583: if (cmpii(gcoeff(dotprd,s,s),gcoeff(dotprd,jcol,jcol)) > 0) s = jcol;
1584: if (icol != s)
1585: { /* Exchange with proper column */
1586: /* Only diagonal of dotprd is updated */
1587: swap(tm2[icol], tm2[s]);
1588: swap(coeff(dotprd,icol,icol), coeff(dotprd,s,s));
1589: }
1590: } /* for icol */
1591: icol=1;
1592: while (icol <= ncol && !signe(gcoeff(dotprd,icol,icol))) icol++;
1593: ncolnz = ncol - icol + 1;
1594: } /* local block */
1595:
1596: if (flag)
1597: {
1598: if (ncolnz == lg((GEN)m[1])-1)
1599: {
1600: tm2 += (ncol-ncolnz);
1601: tm2[0]=evaltyp(t_MAT)|evallg(ncolnz+1);
1602: }
1603: else
1604: {
1605: tm1 = gmul(tm1, tm2);
1606: tm2 = lllint(gmul(m, tm1));
1607: }
1608: }
1609: if (DEBUGLEVEL>6)
1610: fprintferr("lllintpartial output = %Z", gmul(tm1, tm2));
1611: return gerepileupto(ltop1, gmul(tm1, tm2));
1612: }
1613:
1614: GEN
1615: lllintpartial(GEN mat)
1616: {
1617: return lllintpartialall(mat,1);
1618: }
1619:
1620: /********************************************************************/
1621: /** **/
1.2 ! noro 1622: /** LINEAR & ALGEBRAIC DEPENDENCE **/
1.1 noro 1623: /** **/
1624: /********************************************************************/
1625:
1.2 ! noro 1626: GEN pslq(GEN x, long prec);
! 1627: GEN pslqL2(GEN x, long prec);
! 1628:
1.1 noro 1629: GEN
1630: lindep0(GEN x,long bit,long prec)
1631: {
1.2 ! noro 1632: switch (bit)
! 1633: {
! 1634: case 0: return pslq(x,prec);
! 1635: case -1: return lindep(x,prec);
! 1636: case -2: return deplin(x);
! 1637: case -3: return pslqL2(x,prec);
! 1638: default: return lindep2(x,labs(bit));
! 1639: }
1.1 noro 1640: }
1641:
1642: static int
1643: real_indep(GEN re, GEN im, long bitprec)
1644: {
1645: GEN p1 = gsub(gmul((GEN)re[1],(GEN)im[2]),
1646: gmul((GEN)re[2],(GEN)im[1]));
1647: return (!gcmp0(p1) && gexpo(p1) > - bitprec);
1648: }
1649:
1650: GEN
1651: lindep2(GEN x, long bit)
1652: {
1.2 ! noro 1653: long tx=typ(x), lx=lg(x), ly, i, j, e;
! 1654: gpmem_t av = avma;
1.1 noro 1655: GEN re,im,p1,p2;
1656:
1657: if (! is_vec_t(tx)) err(typeer,"lindep2");
1658: if (lx<=2) return cgetg(1,t_VEC);
1659: re = greal(x);
1660: im = gimag(x); bit = (long) (bit/L2SL10);
1661: /* independant over R ? */
1662: if (lx == 3 && real_indep(re,im,bit))
1663: { avma = av; return cgetg(1, t_VEC); }
1664:
1665: if (gcmp0(im)) im = NULL;
1666: ly = im? lx+2: lx+1;
1667: p2=cgetg(lx,t_MAT);
1668: for (i=1; i<lx; i++)
1669: {
1670: p1 = cgetg(ly,t_COL); p2[i] = (long)p1;
1671: for (j=1; j<lx; j++) p1[j] = (i==j)? un: zero;
1672: p1[lx] = lcvtoi(gshift((GEN)re[i],bit),&e);
1673: if (im) p1[lx+1] = lcvtoi(gshift((GEN)im[i],bit),&e);
1674: }
1.2 ! noro 1675: p1 = (GEN)lllint_ip(p2,100)[1];
1.1 noro 1676: p1[0] = evaltyp(t_VEC) | evallg(lx);
1677: return gerepilecopy(av, p1);
1678: }
1679:
1680: #define quazero(x) (gcmp0(x) || (typ(x)==t_REAL && expo(x) < EXP))
1681: GEN
1682: lindep(GEN x, long prec)
1683: {
1684: GEN *b,*be,*bn,**m,qzer;
1685: GEN c1,c2,c3,px,py,pxy,re,im,p1,p2,r,f,em;
1.2 ! noro 1686: long i,j,fl,k, lx = lg(x), tx = typ(x), n = lx-1;
! 1687: gpmem_t av = avma, lim = stack_lim(av,1), av0, av1;
1.1 noro 1688: const long EXP = - bit_accuracy(prec) + 2*n;
1689:
1690: if (! is_vec_t(tx)) err(typeer,"lindep");
1.2 ! noro 1691: if (n <= 1) return cgetg(1,t_VEC);
1.1 noro 1692: x = gmul(x, realun(prec));
1.2 ! noro 1693: re = greal(x);
! 1694: im = gimag(x);
1.1 noro 1695: /* independant over R ? */
1.2 ! noro 1696: if (n == 2 && real_indep(re,im,bit_accuracy(prec)))
1.1 noro 1697: { avma = av; return cgetg(1, t_VEC); }
1.2 ! noro 1698: if (EXP > -10) err(precer,"lindep");
1.1 noro 1699:
1.2 ! noro 1700: qzer = cgetg(lx, t_VECSMALL);
1.1 noro 1701: b = (GEN*)idmat(n);
1702: be= (GEN*)new_chunk(lx);
1703: bn= (GEN*)new_chunk(lx);
1704: m = (GEN**)new_chunk(lx);
1705: for (i=1; i<=n; i++)
1706: {
1.2 ! noro 1707: bn[i] = cgetr(prec+1);
! 1708: be[i] = cgetg(lx,t_COL);
! 1709: m[i] = (GEN*)new_chunk(lx);
! 1710: for (j=1; j<i ; j++) m[i][j] = cgetr(prec+1);
! 1711: for (j=1; j<=n; j++) be[i][j]= lgetr(prec+1);
! 1712: }
! 1713: px = sqscal(re);
! 1714: py = sqscal(im); pxy = gscal(re,im);
! 1715: p1 = mpsub(mpmul(px,py), gsqr(pxy));
! 1716: if (quazero(px)) { re = im; px = py; fl = 1; } else fl = quazero(p1);
1.1 noro 1717: av0 = av1 = avma;
1718: for (i=1; i<=n; i++)
1719: {
1720: p2 = gscal(b[i],re);
1.2 ! noro 1721: if (fl) p2 = gmul(gdiv(p2,px),re);
1.1 noro 1722: else
1723: {
1724: GEN p5,p6,p7;
1.2 ! noro 1725: p5 = gscal(b[i],im);
! 1726: p6 = gdiv(gsub(gmul(py,p2),gmul(pxy,p5)), p1);
! 1727: p7 = gdiv(gsub(gmul(px,p5),gmul(pxy,p2)), p1);
! 1728: p2 = gadd(gmul(p6,re), gmul(p7,im));
1.1 noro 1729: }
1.2 ! noro 1730: if (tx != t_COL) p2 = gtrans(p2);
! 1731: p2 = gsub(b[i],p2);
1.1 noro 1732: for (j=1; j<i; j++)
1.2 ! noro 1733: if (qzer[j]) affrr(bn[j], m[i][j]);
1.1 noro 1734: else
1735: {
1.2 ! noro 1736: gdivz(gscal(b[i],be[j]),bn[j], m[i][j]);
! 1737: p2 = gsub(p2, gmul(m[i][j],be[j]));
1.1 noro 1738: }
1.2 ! noro 1739: gaffect(p2, be[i]);
! 1740: affrr(sqscal(be[i]), bn[i]);
! 1741: qzer[i] = quazero(bn[i]); avma = av1;
1.1 noro 1742: }
1743: while (qzer[n])
1744: {
1745: long e;
1.2 ! noro 1746: if (DEBUGLEVEL>8)
1.1 noro 1747: {
1748: fprintferr("qzer[%ld]=%ld\n",n,qzer[n]);
1.2 ! noro 1749: for (k=1; k<=n; k++)
! 1750: for (i=1; i<k; i++) output(m[k][i]);
1.1 noro 1751: }
1752: em=bn[1]; j=1;
1753: for (i=2; i<n; i++)
1754: {
1755: p1=shiftr(bn[i],i);
1756: if (cmprr(p1,em)>0) { em=p1; j=i; }
1757: }
1.2 ! noro 1758: i=j; k=i+1;
! 1759: avma = av1; r = grndtoi(m[k][i], &e);
1.1 noro 1760: if (e >= 0) err(precer,"lindep");
1761: r = negi(r);
1.2 ! noro 1762: p1 = ZV_lincomb(gun,r, b[k],b[i]);
! 1763: b[k] = b[i];
! 1764: b[i] = p1;
1.1 noro 1765: av1 = avma;
1.2 ! noro 1766: f = addir(r,m[k][i]);
1.1 noro 1767: for (j=1; j<i; j++)
1768: if (!qzer[j])
1769: {
1.2 ! noro 1770: p1 = mpadd(m[k][j], mulir(r,m[i][j]));
! 1771: affrr(m[i][j],m[k][j]);
! 1772: mpaff(p1,m[i][j]);
1.1 noro 1773: }
1.2 ! noro 1774: c1 = addrr(bn[k], mulrr(gsqr(f),bn[i]));
1.1 noro 1775: if (!quazero(c1))
1776: {
1.2 ! noro 1777: c2 = divrr(mulrr(bn[i],f),c1); affrr(c2, m[k][i]);
! 1778: c3 = divrr(bn[k],c1);
! 1779: mulrrz(c3,bn[i], bn[k]); qzer[k] = quazero(bn[k]);
! 1780: affrr(c1, bn[i]); qzer[i] = 0;
1.1 noro 1781: for (j=i+2; j<=n; j++)
1782: {
1.2 ! noro 1783: p1 = addrr(mulrr(m[j][k],c3), mulrr(m[j][i],c2));
! 1784: subrrz(m[j][i],mulrr(f,m[j][k]), m[j][k]);
! 1785: affrr(p1, m[j][i]);
1.1 noro 1786: }
1787: }
1788: else
1789: {
1.2 ! noro 1790: affrr(bn[i], bn[k]); qzer[k] = qzer[i];
! 1791: affrr(c1, bn[i]); qzer[i] = 1;
! 1792: for (j=i+2; j<=n; j++) affrr(m[j][i], m[j][k]);
1.1 noro 1793: }
1794: if (low_stack(lim, stack_lim(av,1)))
1795: {
1796: if(DEBUGMEM>1) err(warnmem,"lindep");
1797: b = (GEN*)gerepilecopy(av0, (GEN)b);
1798: av1 = avma;
1799: }
1800: }
1.2 ! noro 1801: p1 = cgetg(lx,t_COL); p1[n] = un; for (i=1; i<n; i++) p1[i] = zero;
1.1 noro 1802: p2 = (GEN)b; p2[0] = evaltyp(t_MAT) | evallg(lx);
1.2 ! noro 1803: p2 = gauss(gtrans_i(p2),p1);
! 1804: return gerepileupto(av, gtrans(p2));
1.1 noro 1805: }
1806:
1.2 ! noro 1807: /* PSLQ Programs */
1.1 noro 1808:
1.2 ! noro 1809: typedef struct {
! 1810: long vmind, t12, t1234, reda, fin;
! 1811: long ct;
! 1812: } pslq_timer;
! 1813:
! 1814: /* WARNING: for double ** matrices, A[i] is the i-th ROW of A */
! 1815: typedef struct {
! 1816: double *y, **H, **A, **B;
! 1817: double *W; /* scratch space */
! 1818: long n;
! 1819: pslq_timer *T;
! 1820: } pslqL2_M;
! 1821:
! 1822: typedef struct {
! 1823: GEN y, H, A, B;
! 1824: long n, EXP;
! 1825: int flreal;
! 1826: pslq_timer *T;
! 1827: } pslq_M;
1.1 noro 1828:
1.2 ! noro 1829: void
! 1830: init_dalloc()
! 1831: { /* correct alignment for dalloc */
! 1832: avma -= avma % sizeof(double);
! 1833: if (avma < bot) err(errpile);
1.1 noro 1834: }
1835:
1.2 ! noro 1836: double *
! 1837: dalloc(size_t n)
1.1 noro 1838: {
1.2 ! noro 1839: if (avma - bot < n) err(errpile);
! 1840: avma -= n; return (double*)avma;
! 1841: }
1.1 noro 1842:
1.2 ! noro 1843: static double
! 1844: conjd(double x) { return x; }
1.1 noro 1845:
1.2 ! noro 1846: static double
! 1847: sqrd(double a) { return a*a; }
1.1 noro 1848:
1.2 ! noro 1849: static void
! 1850: redall(pslq_M *M, long i, long jsup)
1.1 noro 1851: {
1.2 ! noro 1852: long j, k, n = M->n;
! 1853: GEN t,b;
! 1854: GEN A = M->A, B = M->B, H = M->H, y = M->y;
! 1855: const GEN p1 = (GEN)B[i];
! 1856:
! 1857: for (j=jsup; j>=1; j--)
! 1858: {
! 1859: /* FIXME: use grndtoi */
! 1860: t = ground(gdiv(gcoeff(H,i,j), gcoeff(H,j,j)));
! 1861: if (gcmp0(t)) continue;
! 1862:
! 1863: b = (GEN)B[j];
! 1864: y[j] = ladd((GEN)y[j], gmul(t,(GEN)y[i]));
! 1865: for (k=1; k<=j; k++)
! 1866: coeff(H,i,k) = lsub(gcoeff(H,i,k), gmul(t,gcoeff(H,j,k)));
! 1867: for (k=1; k<=n; k++)
! 1868: {
! 1869: coeff(A,i,k) = lsub(gcoeff(A,i,k), gmul(t,gcoeff(A,j,k)));
! 1870: b[k] = ladd((GEN)b[k], gmul(t,(GEN)p1[k]));
! 1871: }
! 1872: }
1.1 noro 1873: }
1874:
1.2 ! noro 1875: static void
! 1876: redallbar(pslqL2_M *Mbar, long i, long jsup)
1.1 noro 1877: {
1.2 ! noro 1878: long j, k, n = Mbar->n;
! 1879: double t;
! 1880: double *hi = Mbar->H[i], *ai = Mbar->A[i], *hj, *aj;
! 1881:
! 1882: #if 0
! 1883: fprintferr("%ld:\n==\n",i);
! 1884: #endif
! 1885: for (j=jsup; j>=1; j--)
! 1886: {
! 1887: hj = Mbar->H[j];
! 1888: t = floor(0.5 + hi[j] / hj[j]);
! 1889: if (!t) continue;
! 1890: #if 0
! 1891: fprintferr("%15.15e ",t);
! 1892: #endif
! 1893: aj = Mbar->A[j];
! 1894:
! 1895: Mbar->y[j] += t * Mbar->y[i];
! 1896: for (k=1; k<=j; k++) hi[k] -= t * hj[k];
! 1897: for (k=1; k<=n; k++) {
! 1898: ai[k] -= t * aj[k];
! 1899: Mbar->B[k][j] += t * Mbar->B[k][i];
! 1900: }
! 1901: #if 0
! 1902: fprintferr(" %ld:\n",j); dprintmat(Mbar->H,n,n-1);
! 1903: #endif
! 1904: }
1.1 noro 1905: }
1906:
1.2 ! noro 1907: static long
! 1908: vecabsminind(GEN v)
1.1 noro 1909: {
1.2 ! noro 1910: long l = lg(v), m = 1, i;
! 1911: GEN t, la = mpabs((GEN)v[1]);
! 1912: for (i=2; i<l; i++)
1.1 noro 1913: {
1.2 ! noro 1914: t = mpabs((GEN)v[i]);
! 1915: if (mpcmp(t,la) < 0) { la = t; m = i; }
1.1 noro 1916: }
1.2 ! noro 1917: return m;
1.1 noro 1918: }
1919:
1.2 ! noro 1920: static long
! 1921: vecmaxind(GEN v)
1.1 noro 1922: {
1.2 ! noro 1923: long l = lg(v), m = 1, i;
! 1924: GEN t, la = (GEN)v[1];
! 1925: for (i=2; i<l; i++)
! 1926: {
! 1927: t = (GEN)v[i];
! 1928: if (mpcmp(t,la) > 0) { la = t; m = i; }
! 1929: }
! 1930: return m;
1.1 noro 1931: }
1932:
1.2 ! noro 1933: static long
! 1934: vecmaxindbar(double *v, long n)
1.1 noro 1935: {
1.2 ! noro 1936: long m = 1, i;
! 1937: double la = v[1];
! 1938: for (i=2; i<=n; i++)
! 1939: if (v[i] > la) { la = v[i]; m = i; }
! 1940: return m;
1.1 noro 1941: }
1942:
1943: static GEN
1.2 ! noro 1944: maxnorml2(pslq_M *M)
1.1 noro 1945: {
1.2 ! noro 1946: long n = M->n, i, j;
! 1947: GEN ma = gzero, s;
1.1 noro 1948:
1.2 ! noro 1949: for (i=1; i<=n; i++)
1.1 noro 1950: {
1.2 ! noro 1951: s = gzero;
! 1952: for (j=1; j<n; j++) s = gadd(s, gnorm(gcoeff(M->H,i,j)));
! 1953: ma = gmax(ma, s);
1.1 noro 1954: }
1.2 ! noro 1955: return gsqrt(ma, DEFAULTPREC);
1.1 noro 1956: }
1957:
1.2 ! noro 1958: static void
! 1959: init_timer(pslq_timer *T)
1.1 noro 1960: {
1.2 ! noro 1961: T->vmind = T->t12 = T->t1234 = T->reda = T->fin = T->ct = 0;
1.1 noro 1962: }
1963:
1.2 ! noro 1964: static int
! 1965: init_pslq(pslq_M *M, GEN x, long prec)
1.1 noro 1966: {
1.2 ! noro 1967: long lx = lg(x), n = lx-1, i, j, k;
! 1968: GEN s1, s, sinv;
! 1969:
! 1970: M->EXP = - bit_accuracy(prec) + 2*n;
! 1971: M->flreal = (gexpo(gimag(x)) < M->EXP);
! 1972: if (!M->flreal)
! 1973: return 0; /* FIXME */
! 1974: else
! 1975: x = greal(x);
1.1 noro 1976:
1.2 ! noro 1977: if (DEBUGLEVEL>=3) { (void)timer(); init_timer(M->T); }
! 1978: x = gmul(x, realun(prec)); settyp(x,t_VEC);
! 1979: M->n = n;
! 1980: M->A = idmat(n);
! 1981: M->B = idmat(n);
! 1982: s1 = cgetg(lx,t_VEC); s1[n] = lnorm((GEN)x[n]);
! 1983: s = cgetg(lx,t_VEC); s[n] = (long)gabs((GEN)x[n],prec);
! 1984: for (k=n-1; k>=1; k--)
! 1985: {
! 1986: s1[k] = ladd((GEN)s1[k+1], gnorm((GEN)x[k]));
! 1987: s[k] = (long)gsqrt((GEN)s1[k], prec);
! 1988: }
! 1989: sinv = ginv((GEN)s[1]);
! 1990: s = gmul(sinv,s);
! 1991: M->y = gmul(sinv, x);
! 1992: M->H = cgetg(n,t_MAT);
! 1993: for (j=1; j<n; j++)
! 1994: {
! 1995: GEN d, c = cgetg(lx,t_COL);
! 1996:
! 1997: M->H[j] = (long)c;
! 1998: for (i=1; i<j; i++) c[i] = zero;
! 1999: c[j] = ldiv((GEN)s[j+1],(GEN)s[j]);
! 2000: d = gneg( gdiv((GEN)M->y[j], gmul((GEN)s[j],(GEN)s[j+1]) ));
! 2001: for (i=j+1; i<=n; i++) c[i] = lmul(gconj((GEN)M->y[i]), d);
! 2002: }
! 2003: for (i=2; i<=n; i++) redall(M, i, i-1);
! 2004: return 1;
1.1 noro 2005: }
2006:
1.2 ! noro 2007: static void
! 2008: SWAP(pslq_M *M, long m)
1.1 noro 2009: {
1.2 ! noro 2010: long j, n = M->n;
! 2011: swap(M->y[m], M->y[m+1]);
! 2012: swap(M->B[m], M->B[m+1]);
! 2013: for (j=1; j<=n; j++) swap(coeff(M->A,m,j), coeff(M->A,m+1,j));
! 2014: for (j=1; j<n; j++) swap(coeff(M->H,m,j), coeff(M->H,m+1,j));
! 2015: }
! 2016:
! 2017: #define dswap(x,y) { double _t=x; x=y; y=_t; }
! 2018: #define ddswap(x,y) { double* _t=x; x=y; y=_t; }
1.1 noro 2019:
1.2 ! noro 2020: static void
! 2021: SWAPbar(pslqL2_M *M, long m)
! 2022: {
! 2023: long j, n = M->n;
! 2024: dswap(M->y[m], M->y[m+1]);
! 2025: ddswap(M->A[m], M->A[m+1]);
! 2026: ddswap(M->H[m], M->H[m+1]);
! 2027: for (j=1; j<=n; j++) dswap(M->B[j][m], M->B[j][m+1]);
1.1 noro 2028: }
2029:
2030: static GEN
1.2 ! noro 2031: one_step_gen(pslq_M *M, GEN tabga, long prec)
! 2032: {
! 2033: GEN H = M->H, p1, t0, tinv, t1,t2,t3,t4;
! 2034: long n = M->n, i, m;
! 2035:
! 2036: p1 = cgetg(n,t_VEC);
! 2037: for (i=1; i<n; i++) p1[i] = lmul((GEN)tabga[i], gabs(gcoeff(H,i,i),prec));
! 2038: m = vecmaxind(p1);
! 2039: if (DEBUGLEVEL>=4) M->T->vmind += timer();
! 2040: SWAP(M, m);
! 2041: if (m <= n-2)
! 2042: {
! 2043: t0 = gadd(gnorm(gcoeff(H,m,m)), gnorm(gcoeff(H,m,m+1)));
! 2044: tinv = ginv(gsqrt(t0, prec));
! 2045: t1 = gmul(tinv, gcoeff(H,m,m));
! 2046: t2 = gmul(tinv, gcoeff(H,m,m+1));
! 2047: if (DEBUGLEVEL>=4) M->T->t12 += timer();
! 2048: for (i=m; i<=n; i++)
! 2049: {
! 2050: t3 = gcoeff(H,i,m);
! 2051: t4 = gcoeff(H,i,m+1);
! 2052: if (M->flreal)
! 2053: coeff(H,i,m) = ladd(gmul(t1,t3), gmul(t2,t4));
! 2054: else
! 2055: coeff(H,i,m) = ladd(gmul(gconj(t1),t3), gmul(gconj(t2),t4));
! 2056: coeff(H,i,m+1) = lsub(gmul(t1,t4), gmul(t2,t3));
! 2057: }
! 2058: if (DEBUGLEVEL>=4) M->T->t1234 += timer();
1.1 noro 2059: }
1.2 ! noro 2060: for (i=1; i<=n-1; i++)
! 2061: if (gexpo(gcoeff(H,i,i)) <= M->EXP) {
! 2062: m = vecabsminind(M->y); return (GEN)M->B[m];
! 2063: }
! 2064: for (i=m+1; i<=n; i++) redall(M, i, min(i-1,m+1));
1.1 noro 2065:
1.2 ! noro 2066: if (DEBUGLEVEL>=4) M->T->reda += timer();
! 2067: if (gexpo(M->A) >= -M->EXP) return ginv(maxnorml2(M));
! 2068: m = vecabsminind(M->y);
! 2069: if (gexpo((GEN)M->y[m]) <= M->EXP) return (GEN)M->B[m];
1.1 noro 2070:
1.2 ! noro 2071: if (DEBUGLEVEL>=3)
1.1 noro 2072: {
1.2 ! noro 2073: if (DEBUGLEVEL>=4) M->T->fin += timer();
! 2074: M->T->ct++;
! 2075: if ((M->T->ct&0xff) == 0)
! 2076: {
! 2077: if (DEBUGLEVEL == 3)
! 2078: fprintferr("time for ct = %ld : %ld\n",M->T->ct,timer());
! 2079: else
! 2080: fprintferr("time [max,t12,loop,reds,fin] = [%ld, %ld, %ld, %ld, %ld]\n",
! 2081: M->T->vmind, M->T->t12, M->T->t1234, M->T->reda, M->T->fin);
! 2082: }
1.1 noro 2083: }
1.2 ! noro 2084: return NULL; /* nothing interesting */
1.1 noro 2085: }
2086:
2087: static GEN
1.2 ! noro 2088: get_tabga(int flreal, long n, long prec)
1.1 noro 2089: {
1.2 ! noro 2090: GEN ga = mpsqrt( flreal? divrs(stor(4, prec), 3): stor(2, prec) );
! 2091: GEN tabga = cgetg(n,t_VEC);
! 2092: long i;
! 2093: tabga[1] = (long)ga;
! 2094: for (i = 2; i < n; i++) tabga[i] = lmul((GEN)tabga[i-1],ga);
! 2095: return tabga;
! 2096: }
1.1 noro 2097:
2098: GEN
1.2 ! noro 2099: pslq(GEN x, long prec)
1.1 noro 2100: {
1.2 ! noro 2101: GEN tabga, p1;
! 2102: long tx = typ(x);
! 2103: gpmem_t av0 = avma, lim = stack_lim(av0,1), av;
! 2104: pslq_M M;
! 2105: pslq_timer T;
1.1 noro 2106:
1.2 ! noro 2107: if (! is_vec_t(tx)) err(typeer,"pslq");
! 2108: if (lg(x) <= 2) return cgetg(1,t_VEC);
1.1 noro 2109:
1.2 ! noro 2110: M.T = &T; init_pslq(&M, x, prec);
! 2111: if (!M.flreal) return lindep(x, prec); /* FIXME */
1.1 noro 2112:
1.2 ! noro 2113: tabga = get_tabga(M.flreal, M.n, prec);
! 2114: av = avma;
! 2115: if (DEBUGLEVEL>=3) printf("Initialization time = %ld\n",timer());
! 2116: for (;;)
1.1 noro 2117: {
1.2 ! noro 2118: if ((p1 = one_step_gen(&M, tabga, prec)))
! 2119: return gerepilecopy(av0, p1);
! 2120:
! 2121: if (low_stack(lim, stack_lim(av,1)))
1.1 noro 2122: {
1.2 ! noro 2123: if(DEBUGMEM>1) err(warnmem,"pslq");
! 2124: gerepileall(av,4,&M.y,&M.H,&M.A,&M.B);
1.1 noro 2125: }
2126: }
2127: }
2128:
1.2 ! noro 2129: /* W de longueur n-1 */
! 2130:
! 2131: static double
! 2132: dnorml2(double *W, long n, long row)
1.1 noro 2133: {
1.2 ! noro 2134: long i;
! 2135: double s = 0.;
! 2136:
! 2137: for (i=row; i<n; i++) s += W[i]*W[i];
! 2138: return s;
1.1 noro 2139: }
2140:
1.2 ! noro 2141: /* Hbar *= Pbar */
! 2142: static void
! 2143: dmatmul(pslqL2_M *Mbar, double **Pbar, long row)
1.1 noro 2144: {
1.2 ! noro 2145: const long n = Mbar->n; /* > row */
! 2146: long i, j, k;
! 2147: double s, *W = Mbar->W, **H = Mbar->H;
1.1 noro 2148:
1.2 ! noro 2149: for (i = row; i <= n; i++)
1.1 noro 2150: {
1.2 ! noro 2151: for (j = row; j < n; j++)
! 2152: {
! 2153: k = row; s = H[i][k] * Pbar[k][j];
! 2154: for ( ; k < n; k++) s += H[i][k] * Pbar[k][j];
! 2155: W[j] = s;
! 2156: }
! 2157: for (j = row; j < n; j++) H[i][j] = W[j];
1.1 noro 2158: }
2159: }
2160:
1.2 ! noro 2161: /* compute n-1 x n-1 matrix Pbar */
! 2162: static void
! 2163: dmakep(pslqL2_M *Mbar, double **Pbar, long row)
1.1 noro 2164: {
1.2 ! noro 2165: long i, j, n = Mbar->n;
! 2166: double pro, nc, *C = Mbar->H[row], *W = Mbar->W;
1.1 noro 2167:
1.2 ! noro 2168: nc = sqrt(dnorml2(C,n,row));
! 2169: W[row] = (C[row] < 0) ? C[row] - nc : C[row] + nc;
! 2170: for (i=row; i<n; i++) W[i] = C[i];
! 2171: pro = -2.0 / dnorml2(W, n, row);
! 2172: /* must have dnorml2(W,n,row) = 2*nc*(nc+fabs(C[1])) */
! 2173: for (j=row; j<n; j++)
! 2174: {
! 2175: for (i=j+1; i<n; i++)
! 2176: Pbar[j][i] = Pbar[i][j] = pro * W[i] * W[j];
! 2177: Pbar[j][j] = 1.0 + pro * W[j] * W[j];
! 2178: }
! 2179: }
! 2180:
! 2181: static void
! 2182: dLQdec(pslqL2_M *Mbar, double **Pbar)
! 2183: {
! 2184: long row, j, n = Mbar->n;
1.1 noro 2185:
1.2 ! noro 2186: for (row=1; row<n; row++)
! 2187: {
! 2188: dmakep(Mbar, Pbar, row);
! 2189: dmatmul(Mbar, Pbar, row);
! 2190: for (j=row+1; j<n; j++) Mbar->H[row][j] = 0.;
! 2191: }
1.1 noro 2192: }
2193:
1.2 ! noro 2194: void
! 2195: dprintvec(double *V, long m)
1.1 noro 2196: {
2197: long i;
1.2 ! noro 2198: fprintferr("[");
! 2199: for (i=1; i<m; i++) fprintferr("%15.15e, ",V[i]);
! 2200: fprintferr("%15.15e]\n",V[m]); pariflush();
1.1 noro 2201: }
2202:
1.2 ! noro 2203: void
! 2204: dprintmat(double **M, long r, long c)
1.1 noro 2205: {
1.2 ! noro 2206: long i, j;
! 2207: fprintferr("[");
! 2208: for (i=1; i<r; i++)
! 2209: {
! 2210: for (j=1; j<c; j++) fprintferr("%15.15e, ",M[i][j]);
! 2211: fprintferr("%15.15e;\n ",M[i][c]);
! 2212: }
! 2213: for (j=1; j<c; j++) fprintferr("%15.15e, ",M[r][j]);
! 2214: fprintferr("%15.15e]\n",M[r][c]); pariflush();
1.1 noro 2215: }
2216:
1.2 ! noro 2217: static long
! 2218: initializedoubles(pslqL2_M *Mbar, pslq_M *M, long prec)
1.1 noro 2219: {
1.2 ! noro 2220: long i, j, n = Mbar->n;
! 2221: GEN ypro;
! 2222: gpmem_t av = avma;
! 2223:
! 2224: ypro = gdiv(M->y, vecmax(gabs(M->y,prec)));
! 2225: for (i=1; i<=n; i++)
! 2226: {
! 2227: if (gexpo((GEN)ypro[i]) < -0x3ff) return 0;
! 2228: Mbar->y[i] = rtodbl((GEN)ypro[i]);
! 2229: }
! 2230: avma = av;
! 2231: for (j=1; j<=n; j++)
! 2232: for (i=1; i<=n; i++)
! 2233: {
! 2234: if (i==j) Mbar->A[i][j] = Mbar->B[i][j] = 1.;
! 2235: else Mbar->A[i][j] = Mbar->B[i][j] = 0.;
! 2236: if (j < n)
! 2237: {
! 2238: GEN h = gcoeff(M->H,i,j);
! 2239: if (!gcmp0(h) && labs(gexpo(h)) > 0x3ff) return 0;
! 2240: Mbar->H[i][j] = rtodbl(h);
! 2241: }
! 2242: }
! 2243: return 1;
1.1 noro 2244: }
2245:
1.2 ! noro 2246: /* T(arget) := S(ource) */
! 2247: static void
! 2248: storeprecdoubles(pslqL2_M *T, pslqL2_M *S)
1.1 noro 2249: {
1.2 ! noro 2250: long n = T->n, i, j;
1.1 noro 2251:
1.2 ! noro 2252: for (i=1; i<=n; i++)
! 2253: {
! 2254: for (j=1; j<n; j++)
! 2255: {
! 2256: T->H[i][j] = S->H[i][j];
! 2257: T->A[i][j] = S->A[i][j];
! 2258: T->B[i][j] = S->B[i][j];
! 2259: }
! 2260: T->A[i][n] = S->A[i][n];
! 2261: T->B[i][n] = S->B[i][n];
! 2262: T->y[i] = S->y[i];
! 2263: }
1.1 noro 2264: }
2265:
1.2 ! noro 2266: static long
! 2267: checkentries(pslqL2_M *Mbar)
1.1 noro 2268: {
1.2 ! noro 2269: long n = Mbar->n, i, j;
! 2270: double *p1, *p2;
1.1 noro 2271:
1.2 ! noro 2272: for (i=1; i<=n; i++)
! 2273: {
! 2274: if (expodb(Mbar->y[i]) < -46) return 0;
! 2275: p1 = Mbar->A[i];
! 2276: p2 = Mbar->B[i];
! 2277: for (j=1; j<=n; j++)
! 2278: if (expodb(p1[j]) > 43 || expodb(p2[j]) > 43) return 0;
! 2279: }
! 2280: return 1;
1.1 noro 2281: }
2282:
1.2 ! noro 2283: static long
! 2284: applybar(pslq_M *M, pslqL2_M *Mbar, GEN Abargen, GEN Bbargen)
1.1 noro 2285: {
1.2 ! noro 2286: long n = Mbar->n, i, j;
! 2287: double *p1, *p2;
1.1 noro 2288:
1.2 ! noro 2289: for (i=1; i<=n; i++)
! 2290: {
! 2291: p1 = Mbar->A[i];
! 2292: p2 = Mbar->B[i];
! 2293: for (j=1; j<=n; j++)
1.1 noro 2294: {
1.2 ! noro 2295: if (expodb(p1[j]) >= 52 || expodb(p2[j]) >= 52) return 0;
! 2296: coeff(Abargen,i,j) = (long)mpent(dbltor(p1[j]));
! 2297: coeff(Bbargen,i,j) = (long)mpent(dbltor(p2[j]));
1.1 noro 2298: }
2299: }
1.2 ! noro 2300: M->y = gmul(M->y, Bbargen);
! 2301: M->B = gmul(M->B, Bbargen);
! 2302: M->A = gmul(Abargen, M->A);
! 2303: M->H = gmul(Abargen, M->H); return 1;
1.1 noro 2304: }
2305:
2306: static GEN
1.2 ! noro 2307: checkend(pslq_M *M, long prec)
1.1 noro 2308: {
1.2 ! noro 2309: long i, m, n = M->n;
! 2310:
! 2311: for (i=1; i<=n-1; i++)
! 2312: if (gexpo(gcoeff(M->H,i,i)) <= M->EXP)
! 2313: {
! 2314: m = vecabsminind(M->y);
! 2315: return (GEN)M->B[m];
! 2316: }
! 2317: if (gexpo(M->A) >= -M->EXP)
! 2318: return ginv( maxnorml2(M) );
! 2319: m = vecabsminind(M->y);
! 2320: if (gexpo((GEN)M->y[m]) <= M->EXP) return (GEN)M->B[m];
! 2321: return NULL;
1.1 noro 2322: }
2323:
1.2 ! noro 2324: GEN
! 2325: pslqL2(GEN x, long prec)
1.1 noro 2326: {
1.2 ! noro 2327: GEN Abargen, Bbargen, tabga, p1;
! 2328: long lx = lg(x), tx = typ(x), n = lx-1, i, m, ctpro, flreal, flit;
! 2329: gpmem_t av0 = avma, lim = stack_lim(av0,1), av;
! 2330: double *tabgabar, gabar, tinvbar, t1bar, t2bar, t3bar, t4bar;
! 2331: double **Pbar, **Abar, **Bbar, **Hbar, *ybar;
! 2332: pslqL2_M Mbar, Mbarst;
! 2333: pslq_M M;
! 2334: pslq_timer T;
! 2335:
! 2336: if (! is_vec_t(tx)) err(typeer,"pslq");
! 2337: if (n <= 1) return cgetg(1,t_COL);
! 2338: M.T = &T; init_pslq(&M, x, prec);
! 2339: if (!M.flreal) return lindep(x, prec);
! 2340:
! 2341: flreal = M.flreal;
! 2342: tabga = get_tabga(flreal, n, prec);
! 2343: Abargen = idmat(n);
! 2344: Bbargen = idmat(n);
! 2345:
! 2346: Mbarst.n = Mbar.n = n;
! 2347: Mbar.A = Abar = (double**)new_chunk(n+1);
! 2348: Mbar.B = Bbar = (double**)new_chunk(n+1);
! 2349: Mbar.H = Hbar = (double**)new_chunk(n+1);
! 2350: Mbarst.A = (double**)new_chunk(n+1);
! 2351: Mbarst.B = (double**)new_chunk(n+1);
! 2352: Mbarst.H = (double**)new_chunk(n+1);
! 2353: Pbar = (double**)new_chunk(n);
1.1 noro 2354:
1.2 ! noro 2355: tabgabar = dalloc((n+1)*sizeof(double));
! 2356: Mbar.y = ybar = dalloc((n+1)*sizeof(double));
! 2357: Mbarst.y = dalloc((n+1)*sizeof(double));
! 2358:
! 2359: Mbar.W = dalloc((n+1)*sizeof(double));
! 2360: for (i=1; i< n; i++) Pbar[i] = dalloc((n+1)*sizeof(double));
! 2361: for (i=1; i<=n; i++) Abar[i] = dalloc((n+1)*sizeof(double));
! 2362: for (i=1; i<=n; i++) Bbar[i] = dalloc((n+1)*sizeof(double));
! 2363: for (i=1; i<=n; i++) Hbar[i] = dalloc(n*sizeof(double));
! 2364: for (i=1; i<=n; i++) Mbarst.A[i] = dalloc((n+1)*sizeof(double));
! 2365: for (i=1; i<=n; i++) Mbarst.B[i] = dalloc((n+1)*sizeof(double));
! 2366: for (i=1; i<=n; i++) Mbarst.H[i] = dalloc(n*sizeof(double));
1.1 noro 2367:
1.2 ! noro 2368: gabar = gtodouble((GEN)tabga[1]); tabgabar[1] = gabar;
! 2369: for (i=2; i<n; i++) tabgabar[i] = tabgabar[i-1]*gabar;
1.1 noro 2370:
1.2 ! noro 2371: av = avma;
! 2372: if (DEBUGLEVEL>=3) printf("Initialization time = %ld\n",timer());
! 2373: RESTART:
! 2374: flit = initializedoubles(&Mbar, &M, prec);
! 2375: storeprecdoubles(&Mbarst, &Mbar);
! 2376: if (flit) dLQdec(&Mbar, Pbar);
! 2377: ctpro = 0;
! 2378: for (;;)
1.1 noro 2379: {
1.2 ! noro 2380: if (low_stack(lim, stack_lim(av,1)))
1.1 noro 2381: {
1.2 ! noro 2382: if(DEBUGMEM>1) err(warnmem,"pslq");
! 2383: gerepileall(av,4,&M.y,&M.H,&M.A,&M.B);
1.1 noro 2384: }
1.2 ! noro 2385: if (flit)
! 2386: {
! 2387: ctpro++;
! 2388: for (i=1; i<n; i++) Mbar.W[i] = tabgabar[i]*fabs(Hbar[i][i]);
! 2389: m = vecmaxindbar(Mbar.W, n-1);
! 2390: SWAPbar(&Mbar, m);
! 2391: if (m <= n-2)
! 2392: {
! 2393: tinvbar = 1.0 / sqrt(sqrd(Hbar[m][m]) + sqrd(Hbar[m][m+1]));
! 2394: t1bar = tinvbar*Hbar[m][m];
! 2395: t2bar = tinvbar*Hbar[m][m+1];
! 2396: if (DEBUGLEVEL>=4) T.t12 += timer();
! 2397: for (i=m; i<=n; i++)
! 2398: {
! 2399: t3bar = Hbar[i][m];
! 2400: t4bar = Hbar[i][m+1];
! 2401: if (flreal)
! 2402: Hbar[i][m] = t1bar*t3bar + t2bar*t4bar;
! 2403: else
! 2404: Hbar[i][m] = conjd(t1bar)*t3bar + conjd(t2bar)*t4bar;
! 2405: Hbar[i][m+1] = t1bar*t4bar - t2bar*t3bar;
! 2406: }
! 2407: if (DEBUGLEVEL>=4) T.t1234 += timer();
! 2408: }
1.1 noro 2409:
1.2 ! noro 2410: flit = checkentries(&Mbar);
! 2411: if (flit)
! 2412: {
! 2413: storeprecdoubles(&Mbarst, &Mbar);
! 2414: for (i=m+1; i<=n; i++) redallbar(&Mbar, i, min(i-1,m+1));
! 2415: }
! 2416: else
! 2417: {
! 2418: if (applybar(&M, &Mbar, Abargen,Bbargen))
! 2419: {
! 2420: if ( (p1 = checkend(&M,prec)) ) return gerepilecopy(av0, p1);
! 2421: goto RESTART;
! 2422: }
! 2423: else
! 2424: {
! 2425: if (ctpro == 1) goto DOGEN;
! 2426: storeprecdoubles(&Mbar, &Mbarst); /* restore */
! 2427: if (! applybar(&M, &Mbar, Abargen,Bbargen)) err(bugparier,"pslqL2");
! 2428: if ( (p1 = checkend(&M, prec)) ) return gerepilecopy(av0, p1);
! 2429: goto RESTART;
! 2430: }
! 2431: }
! 2432: }
! 2433: else
1.1 noro 2434: {
1.2 ! noro 2435: DOGEN:
! 2436: if ((p1 = one_step_gen(&M, tabga, prec)))
! 2437: return gerepilecopy(av, p1);
1.1 noro 2438: }
2439: }
2440: }
2441:
1.2 ! noro 2442: /* x is a vector of elts of a p-adic field */
1.1 noro 2443: GEN
1.2 ! noro 2444: plindep(GEN x)
1.1 noro 2445: {
1.2 ! noro 2446: long i, j, prec = VERYBIGINT, lx = lg(x)-1, ly, v;
! 2447: gpmem_t av = avma;
! 2448: GEN p = NULL, pn,p1,m,a;
! 2449:
! 2450: if (lx < 2) return cgetg(1,t_VEC);
! 2451: for (i=1; i<=lx; i++)
1.1 noro 2452: {
1.2 ! noro 2453: p1 = (GEN)x[i];
! 2454: if (typ(p1) != t_PADIC) continue;
1.1 noro 2455:
1.2 ! noro 2456: j = precp(p1); if (j < prec) prec = j;
! 2457: if (!p) p = (GEN)p1[2];
! 2458: else if (!egalii(p, (GEN)p1[2]))
! 2459: err(talker,"inconsistent primes in plindep");
1.1 noro 2460: }
1.2 ! noro 2461: if (!p) err(talker,"not a p-adic vector in plindep");
! 2462: v = ggval(x,p); pn = gpowgs(p,prec);
! 2463: if (v != 0) x = gmul(x, gpowgs(p, -v));
! 2464: x = lift_intern(gmul(x, gmodulcp(gun, pn)));
1.1 noro 2465:
1.2 ! noro 2466: ly = 2*lx - 1;
! 2467: m = cgetg(ly+1,t_MAT);
! 2468: for (j=1; j<=ly; j++)
1.1 noro 2469: {
1.2 ! noro 2470: p1 = cgetg(lx+1, t_COL); m[j] = (long)p1;
! 2471: for (i=1; i<=lx; i++) p1[i] = zero;
1.1 noro 2472: }
1.2 ! noro 2473: a = negi((GEN)x[1]);
! 2474: for (i=1; i<lx; i++)
1.1 noro 2475: {
1.2 ! noro 2476: coeff(m,1+i,i) = (long)a;
! 2477: coeff(m,1 ,i) = x[i+1];
1.1 noro 2478: }
1.2 ! noro 2479: for (i=1; i<=lx; i++) coeff(m,i,lx-1+i) = (long)pn;
! 2480: p1 = lllint(m);
! 2481: return gerepileupto(av, gmul(m, (GEN)p1[1]));
1.1 noro 2482: }
2483:
2484: GEN
1.2 ! noro 2485: algdep0(GEN x, long n, long bit, long prec)
1.1 noro 2486: {
1.2 ! noro 2487: long tx=typ(x), i, k;
! 2488: gpmem_t av;
! 2489: GEN y,p1;
1.1 noro 2490:
1.2 ! noro 2491: if (! is_scalar_t(tx)) err(typeer,"algdep0");
! 2492: if (tx==t_POLMOD) { y=forcecopy((GEN)x[1]); setvarn(y,0); return y; }
! 2493: if (gcmp0(x)) return gzero;
! 2494: if (!n) return gun;
1.1 noro 2495:
1.2 ! noro 2496: av=avma; p1=cgetg(n+2,t_COL);
! 2497: p1[1] = un;
! 2498: p1[2] = (long)x; /* n >= 1 */
! 2499: for (i=3; i<=n+1; i++) p1[i]=lmul((GEN)p1[i-1],x);
! 2500: if (typ(x) == t_PADIC)
! 2501: p1 = plindep(p1);
! 2502: else
! 2503: {
! 2504: switch(bit)
! 2505: {
! 2506: case 0: p1 = pslq(p1,prec); break;
! 2507: case -1: p1 = lindep(p1,prec); break;
! 2508: case -2: p1 = deplin(p1); break;
! 2509: default: p1 = lindep2(p1,bit);
! 2510: }
! 2511: }
! 2512: if ((!bit) && (typ(p1) == t_REAL))
! 2513: {
! 2514: y = gcopy(p1); return gerepileupto(av,y);
! 2515: }
! 2516: if (lg(p1) < 2)
! 2517: err(talker,"higher degree than expected in algdep");
! 2518: y=cgetg(n+3,t_POL);
! 2519: y[1] = evalsigne(1) | evalvarn(0);
! 2520: k=1; while (gcmp0((GEN)p1[k])) k++;
! 2521: for (i=0; i<=n+1-k; i++) y[i+2] = p1[k+i];
! 2522: (void)normalizepol_i(y, n+4-k);
! 2523: y = (gsigne(leading_term(y)) > 0)? gcopy(y): gneg(y);
! 2524: return gerepileupto(av,y);
1.1 noro 2525: }
2526:
2527: GEN
1.2 ! noro 2528: algdep2(GEN x, long n, long bit)
1.1 noro 2529: {
1.2 ! noro 2530: return algdep0(x,n,bit,0);
1.1 noro 2531: }
2532:
2533: GEN
1.2 ! noro 2534: algdep(GEN x, long n, long prec)
1.1 noro 2535: {
1.2 ! noro 2536: return algdep0(x,n,0,prec);
1.1 noro 2537: }
2538:
1.2 ! noro 2539: /********************************************************************/
! 2540: /** **/
! 2541: /** INTEGRAL KERNEL (LLL REDUCED) **/
! 2542: /** **/
! 2543: /********************************************************************/
1.1 noro 2544:
2545: GEN
1.2 ! noro 2546: matkerint0(GEN x, long flag)
1.1 noro 2547: {
1.2 ! noro 2548: switch(flag)
! 2549: {
! 2550: case 0: return kerint(x);
! 2551: case 1: return kerint1(x);
! 2552: default: err(flagerr,"matkerint");
! 2553: }
! 2554: return NULL; /* not reached */
1.1 noro 2555: }
2556:
2557: GEN
1.2 ! noro 2558: kerint1(GEN x)
1.1 noro 2559: {
1.2 ! noro 2560: gpmem_t av = avma;
! 2561: return gerepilecopy(av, lllint_ip(matrixqz3(ker(x)), 100));
1.1 noro 2562: }
2563:
2564: GEN
1.2 ! noro 2565: kerint(GEN x)
1.1 noro 2566: {
1.2 ! noro 2567: gpmem_t av = avma;
! 2568: GEN fl, junk, h = lllint_i(x, 0, 0, &junk, &fl, NULL);
! 2569: if (h) h = lll_finish(h,fl, lll_KER);
! 2570: else h = lll_trivial(x, lll_KER);
! 2571: if (lg(h)==1) { avma = av; return cgetg(1, t_MAT); }
! 2572: return gerepilecopy(av, lllint_ip(h, 100));
1.1 noro 2573: }
2574:
2575: /********************************************************************/
2576: /** **/
2577: /** MINIM **/
2578: /** **/
2579: /********************************************************************/
2580: /* x is a non-empty matrix, y is a compatible VECSMALL (dimension > 0). */
2581: GEN
2582: gmul_mat_smallvec(GEN x, GEN y)
2583: {
1.2 ! noro 2584: long c = lg(x), l = lg(x[1]), i, j;
! 2585: gpmem_t av;
1.1 noro 2586: GEN z=cgetg(l,t_COL), s;
2587:
2588: for (i=1; i<l; i++)
2589: {
2590: av = avma; s = gmulgs(gcoeff(x,i,1),y[1]);
2591: for (j=2; j<c; j++)
2592: if (y[j]) s = gadd(s, gmulgs(gcoeff(x,i,j),y[j]));
2593: z[i] = lpileupto(av,s);
2594: }
2595: return z;
2596: }
2597:
2598: /* same, x integral */
2599: GEN
2600: gmul_mati_smallvec(GEN x, GEN y)
2601: {
1.2 ! noro 2602: long c = lg(x), l = lg(x[1]), i, j;
! 2603: gpmem_t av;
1.1 noro 2604: GEN z=cgetg(l,t_COL), s;
2605:
2606: for (i=1; i<l; i++)
2607: {
2608: av = avma; s = mulis(gcoeff(x,i,1),y[1]);
2609: for (j=2; j<c; j++)
2610: if (y[j]) s = addii(s, mulis(gcoeff(x,i,j),y[j]));
1.2 ! noro 2611: z[i] = lpileuptoint(av,s);
1.1 noro 2612: }
2613: return z;
2614: }
2615:
2616: void
2617: minim_alloc(long n, double ***q, GEN *x, double **y, double **z, double **v)
2618: {
2619: long i, s;
2620:
2621: *x = cgetg(n, t_VECSMALL);
2622: *q = (double**) new_chunk(n);
1.2 ! noro 2623: s = n * sizeof(double);
! 2624: init_dalloc();
! 2625: *y = dalloc(s);
! 2626: *z = dalloc(s);
! 2627: *v = dalloc(s);
! 2628: for (i=1; i<n; i++) (*q)[i] = dalloc(s);
1.1 noro 2629: }
2630:
2631: /* Minimal vectors for the integral definite quadratic form: a.
2632: * Result u:
2633: * u[1]= Number of vectors of square norm <= BORNE
2634: * u[2]= maximum norm found
2635: * u[3]= list of vectors found (at most STOCKMAX)
2636: *
2637: * If BORNE = gzero: Minimal non-zero vectors.
2638: * flag = min_ALL, as above
2639: * flag = min_FIRST, exits when first suitable vector is found.
2640: * flag = min_PERF, only compute rank of the family of v.v~ (v min.)
2641: */
2642: static GEN
2643: minim00(GEN a, GEN BORNE, GEN STOCKMAX, long flag)
2644: {
2645: GEN x,res,p1,u,r,liste,gnorme,invp,V, *gptr[2];
1.2 ! noro 2646: long n = lg(a), i, j, k, s, maxrank;
! 2647: gpmem_t av0 = avma, av1, av, tetpil, lim;
1.1 noro 2648: double p,maxnorm,BOUND,*v,*y,*z,**q, eps = 0.000001;
2649:
2650: maxrank = 0; res = V = invp = NULL; /* gcc -Wall */
2651: switch(flag)
2652: {
2653: case min_FIRST:
2654: if (gcmp0(BORNE)) err(talker,"bound = 0 in minim2");
2655: res = cgetg(3,t_VEC); break;
2656: case min_ALL: res = cgetg(4,t_VEC); break;
2657: case min_PERF: break;
2658: default: err(talker, "incorrect flag in minim00");
2659: }
2660: if (n == 1)
2661: {
2662: switch(flag)
2663: {
2664: case min_FIRST: avma=av0; return cgetg(1,t_VEC);
2665: case min_PERF: avma=av0; return gzero;
2666: }
2667: res[1]=res[2]=zero;
2668: res[3]=lgetg(1,t_MAT); return res;
2669: }
2670:
2671: av = avma;
2672: minim_alloc(n, &q, &x, &y, &z, &v);
2673: av1 = avma;
2674:
2675: u = lllgramint(a);
2676: if (lg(u) != n)
2677: err(talker,"not a definite form in minim00");
2678: a = qf_base_change(a,u,1);
2679:
2680: n--;
2681: a = gmul(a, realun(DEFAULTPREC)); r = sqred1(a);
2682: if (DEBUGLEVEL>4) { fprintferr("minim: r = "); outerr(r); }
2683: for (j=1; j<=n; j++)
2684: {
2685: v[j] = rtodbl(gcoeff(r,j,j));
2686: for (i=1; i<j; i++) q[i][j] = rtodbl(gcoeff(r,i,j));
2687: }
2688:
2689: if (flag==min_PERF || gcmp0(BORNE))
2690: {
2691: double c, b = rtodbl(gcoeff(a,1,1));
2692:
2693: for (i=2; i<=n; i++)
2694: { c=rtodbl(gcoeff(a,i,i)); if (c<b) b=c; }
2695: BOUND = b+eps;
2696: BORNE = ground(dbltor(BOUND));
2697: maxnorm = -1.; /* don't update maxnorm */
2698: }
2699: else
2700: {
2701: BORNE = gfloor(BORNE);
2702: BOUND = gtodouble(BORNE)+eps;
2703: maxnorm = 0.;
2704: }
2705:
2706: switch(flag)
2707: {
2708: case min_ALL:
2709: maxrank=itos(STOCKMAX);
2710: liste = new_chunk(1+maxrank);
2711: break;
2712: case min_PERF:
2713: BORNE = gerepileupto(av1,BORNE);
2714: maxrank = (n*(n+1))>>1;
2715: liste = cgetg(1+maxrank, t_VECSMALL);
2716: V = cgetg(1+maxrank, t_VECSMALL);
2717: for (i=1; i<=maxrank; i++) liste[i]=0;
2718: }
2719:
2720: s=0; av1=avma; lim = stack_lim(av1,1);
2721: k = n; y[n] = z[n] = 0;
2722: x[n] = (long) sqrt(BOUND/v[n]);
2723: if (flag == min_PERF) invp = idmat(maxrank);
2724: for(;;x[1]--)
2725: {
2726: do
2727: {
2728: if (k>1)
2729: {
2730: long l = k-1;
2731: z[l]=0;
2732: for (j=k; j<=n; j++) z[l] += q[l][j]*x[j];
2733: p = x[k]+z[k];
2734: y[l] = y[k] + p*p*v[k];
2735: x[l] = (long) floor(sqrt((BOUND-y[l])/v[l])-z[l]);
2736: k = l;
2737: }
2738: for(;;)
2739: {
2740: p = x[k]+z[k];
2741: if (y[k] + p*p*v[k] <= BOUND) break;
2742: k++; x[k]--;
2743: }
2744: }
2745: while (k>1);
2746: if (! x[1] && y[1]<=eps) break;
2747: p = x[1]+z[1]; p = y[1] + p*p*v[1]; /* norm(x) */
2748: if (maxnorm >= 0)
2749: {
2750: if (flag == min_FIRST)
2751: {
2752: gnorme = dbltor(p);
2753: tetpil=avma; gnorme = ground(gnorme); r = gmul_mati_smallvec(u,x);
2754: gptr[0]=&gnorme; gptr[1]=&r; gerepilemanysp(av,tetpil,gptr,2);
2755: res[1]=(long)gnorme;
2756: res[2]=(long)r; return res;
2757: }
2758: if (p > maxnorm) maxnorm = p;
2759: }
2760: else
2761: {
1.2 ! noro 2762: gpmem_t av2 = avma;
1.1 noro 2763: gnorme = ground(dbltor(p));
2764: if (gcmp(gnorme,BORNE) >= 0) avma = av2;
2765: else
2766: {
2767: BOUND=gtodouble(gnorme)+eps; s=0;
2768: affii(gnorme,BORNE); avma=av1;
2769: if (flag == min_PERF) invp = idmat(maxrank);
2770: }
2771: }
2772: s++;
2773:
2774: switch(flag)
2775: {
2776: case min_ALL:
2777: if (s<=maxrank)
2778: {
2779: p1 = new_chunk(n+1); liste[s] = (long) p1;
2780: for (i=1; i<=n; i++) p1[i] = x[i];
2781: }
2782: break;
2783:
2784: case min_PERF:
2785: {
1.2 ! noro 2786: long I=1;
! 2787: gpmem_t av2=avma;
1.1 noro 2788:
2789: for (i=1; i<=n; i++)
2790: for (j=i; j<=n; j++,I++) V[I] = x[i]*x[j];
2791: if (! addcolumntomatrix(V,invp,liste))
2792: {
2793: if (DEBUGLEVEL>1) { fprintferr("."); flusherr(); }
2794: s--; avma=av2; continue;
2795: }
2796:
2797: if (DEBUGLEVEL>1) { fprintferr("*"); flusherr(); }
2798: if (s == maxrank)
2799: {
2800: if (DEBUGLEVEL>1) { fprintferr("\n"); flusherr(); }
2801: avma=av0; return stoi(s);
2802: }
2803:
2804: if (low_stack(lim, stack_lim(av1,1)))
2805: {
2806: if(DEBUGMEM>1) err(warnmem,"minim00, rank>=%ld",s);
2807: invp = gerepilecopy(av1, invp);
2808: }
2809: }
2810: }
2811: }
2812: switch(flag)
2813: {
2814: case min_FIRST:
2815: avma=av0; return cgetg(1,t_VEC);
2816: case min_PERF:
2817: if (DEBUGLEVEL>1) { fprintferr("\n"); flusherr(); }
2818: avma=av0; return stoi(s);
2819: }
2820: k = min(s,maxrank);
2821:
2822: tetpil = avma; p1=cgetg(k+1,t_MAT);
2823: for (j=1; j<=k; j++)
2824: p1[j] = (long) gmul_mati_smallvec(u,(GEN)liste[j]);
2825: liste = p1;
2826: r = (maxnorm >= 0) ? ground(dbltor(maxnorm)): BORNE;
2827:
2828: r=icopy(r); gptr[0]=&r; gptr[1]=&liste;
2829: gerepilemanysp(av,tetpil,gptr,2);
2830: res[1]=lstoi(s<<1);
2831: res[2]=(long)r;
2832: res[3]=(long)liste; return res;
2833: }
2834:
2835: GEN
2836: qfminim0(GEN a, GEN borne, GEN stockmax, long flag, long prec)
2837: {
2838: switch(flag)
2839: {
2840: case 0: return minim00(a,borne,stockmax,min_ALL);
2841: case 1: return minim00(a,borne,gzero ,min_FIRST);
2842: case 2: return fincke_pohst(a,borne,itos(stockmax),0,prec,NULL);
2843: default: err(flagerr,"qfminim");
2844: }
2845: return NULL; /* not reached */
2846: }
2847:
2848: GEN
2849: minim(GEN a, GEN borne, GEN stockmax)
2850: {
2851: return minim00(a,borne,stockmax,min_ALL);
2852: }
2853:
2854: GEN
2855: minim2(GEN a, GEN borne, GEN stockmax)
2856: {
2857: return minim00(a,borne,stockmax,min_FIRST);
2858: }
2859:
2860: GEN
2861: perf(GEN a)
2862: {
2863: return minim00(a,gzero,gzero,min_PERF);
2864: }
2865:
2866: /* general program for positive definit quadratic forms (real coeffs).
2867: * One needs BORNE != 0; LLL reduction done in fincke_pohst().
1.2 ! noro 2868: * If (noer) return NULL on precision problems (no error).
1.1 noro 2869: * If (check != NULL consider only vectors passing the check [ assumes
2870: * stockmax > 0 and we only want the smallest possible vectors ] */
2871: static GEN
1.2 ! noro 2872: smallvectors(GEN a, GEN BORNE, long stockmax, long noer, FP_chk_fun *CHECK)
1.1 noro 2873: {
1.2 ! noro 2874: long N, n, i, j, k, s, epsbit, prec, checkcnt = 1;
! 2875: gpmem_t av, av1, lim;
1.1 noro 2876: GEN u,S,x,y,z,v,q,norme1,normax1,borne1,borne2,eps,p1,alpha,norms,dummy;
1.2 ! noro 2877: GEN (*check)(void *,GEN) = CHECK? CHECK->f: NULL;
! 2878: void *data = CHECK? CHECK->data: NULL;
1.1 noro 2879: int skipfirst = CHECK? CHECK->skipfirst: 0;
2880:
2881: if (DEBUGLEVEL)
2882: fprintferr("smallvectors looking for norm <= %Z\n",gprec_w(BORNE,3));
2883:
1.2 ! noro 2884: q = sqred1intern(a, noer);
1.1 noro 2885: if (q == NULL) return NULL;
2886: if (DEBUGLEVEL>5) fprintferr("q = %Z",q);
2887: prec = gprecision(q);
2888: epsbit = bit_accuracy(prec) >> 1;
1.2 ! noro 2889: eps = real2n(-epsbit, prec);
1.1 noro 2890: alpha = dbltor(0.95);
2891: normax1 = gzero;
2892: borne1= gadd(BORNE,eps);
2893: borne2 = mpmul(borne1,alpha);
2894: N = lg(q); n = N-1;
2895: v = cgetg(N,t_VEC);
2896: dummy = cgetg(1,t_STR);
2897:
2898: av = avma; lim = stack_lim(av,1);
2899: if (check) norms = cgetg(stockmax+1,t_VEC);
2900: S = cgetg(stockmax+1,t_VEC);
2901: x = cgetg(N,t_COL);
2902: y = cgetg(N,t_COL);
2903: z = cgetg(N,t_COL);
2904: for (i=1; i<N; i++) { v[i] = coeff(q,i,i); x[i]=y[i]=z[i] = zero; }
2905:
2906: x[n] = lmpent(mpsqrt(gdiv(borne1,(GEN)v[n])));
2907: if (DEBUGLEVEL>3) { fprintferr("\nx[%ld] = %Z\n",n,x[n]); flusherr(); }
2908:
2909: s = 0; k = n;
2910: for(;; x[k] = laddis((GEN)x[k],-1)) /* main */
2911: {
2912: do
2913: {
2914: int fl = 0;
2915: if (k > 1)
2916: {
2917: av1=avma; k--; p1 = mpmul(gcoeff(q,k,k+1),(GEN)x[k+1]);
2918: for (j=k+2; j<N; j++)
2919: p1 = mpadd(p1, mpmul(gcoeff(q,k,j),(GEN)x[j]));
2920: z[k] = (long)gerepileuptoleaf(av1,p1);
2921:
2922: av1=avma; p1 = gsqr(mpadd((GEN)x[k+1],(GEN)z[k+1]));
2923: p1 = mpadd((GEN)y[k+1], mpmul(p1,(GEN)v[k+1]));
2924: y[k] = (long)gerepileuptoleaf(av1, p1);
2925:
2926: av1=avma; p1 = mpsub(borne1, (GEN)y[k]);
2927: if (signe(p1) < 0) { avma=av1; fl = 1; }
2928: else
2929: {
2930: p1 = mpadd(eps,mpsub(mpsqrt(gdiv(p1,(GEN)v[k])), (GEN)z[k]));
2931: x[k] = (long)gerepileuptoleaf(av1,mpent(p1));
2932: }
2933: /* reject the [x_1,...,x_skipfirst,0,...,0] */
2934: if (k <= skipfirst && !signe(y[skipfirst])) goto END;
2935: }
2936: for(;; x[k] = laddis((GEN)x[k],-1))
2937: {
2938: if (!fl)
2939: {
2940: av1 = avma; /* p1 >= 0 */
2941: p1 = mpmul((GEN)v[k], gsqr(mpadd((GEN)x[k],(GEN)z[k])));
2942: i = mpcmp(mpsub(mpadd(p1,(GEN)y[k]), borne1), gmul2n(p1,-epsbit));
2943: avma = av1; if (i <= 0) break;
2944: }
2945: k++; fl=0;
2946: }
2947:
2948: if (low_stack(lim, stack_lim(av,1)))
2949: {
2950: int cnt = 4;
2951: if(DEBUGMEM>1) err(warnmem,"smallvectors");
2952: if (stockmax)
2953: { /* initialize to dummy values */
2954: GEN T = S;
2955: for (i=s+1; i<=stockmax; i++) S[i]=(long)dummy;
2956: S = gclone(S); if (isclone(T)) gunclone(T);
2957: }
2958: if (check)
2959: {
1.2 ! noro 2960: cnt+=3;
1.1 noro 2961: for (i=s+1; i<=stockmax; i++) norms[i]=(long)dummy;
2962: }
1.2 ! noro 2963: gerepileall(av,cnt,&x,&y,&z,&normax1,&borne1,&borne2,&norms);
1.1 noro 2964: }
2965: if (DEBUGLEVEL>3)
2966: {
2967: if (DEBUGLEVEL>5) fprintferr("%ld ",k);
2968: if (k==n) fprintferr("\nx[%ld] = %Z\n",n,x[n]);
2969: flusherr();
2970: }
2971: }
2972: while (k > 1);
2973:
2974: /* x = 0: we're done */
2975: if (!signe(x[1]) && !signe(y[1])) goto END;
2976:
2977: av1=avma; p1 = gsqr(mpadd((GEN)x[1],(GEN)z[1]));
2978: norme1 = mpadd((GEN)y[1], mpmul(p1, (GEN)v[1]));
2979: if (mpcmp(norme1,borne1) > 0) { avma=av1; continue; /* main */ }
2980:
2981: norme1 = gerepileupto(av1,norme1);
2982: if (check)
2983: {
2984: if (checkcnt < 5 && mpcmp(norme1, borne2) < 0)
2985: {
2986: if (check(data,x))
2987: {
2988: borne1 = mpadd(norme1,eps);
2989: borne2 = mpmul(borne1,alpha);
2990: s = 0; checkcnt = 0;
2991: }
2992: else { checkcnt++ ; continue; /* main */}
2993: }
2994: }
2995: else if (mpcmp(norme1,normax1) > 0) normax1 = norme1;
2996: if (++s <= stockmax)
2997: {
2998: if (check) norms[s] = (long)norme1;
2999: S[s] = (long)dummycopy(x);
1.2 ! noro 3000: if (s == stockmax)
1.1 noro 3001: {
1.2 ! noro 3002: gpmem_t av2 = avma;
! 3003: GEN per;
! 3004:
! 3005: if (!check) goto END;
! 3006: per = sindexsort(norms);
1.1 noro 3007: if (DEBUGLEVEL) fprintferr("sorting...\n");
3008: for (i=1; i<=s; i++)
3009: {
3010: long k = per[i];
3011: if (check(data,(GEN)S[k]))
3012: {
1.2 ! noro 3013: S[1] = S[k]; avma = av2;
1.1 noro 3014: borne1 = mpadd(norme1,eps);
3015: borne2 = mpmul(borne1,alpha);
3016: s = 1; checkcnt = 0; break;
3017: }
3018: }
3019: if (checkcnt) s = 0;
3020: }
3021: }
3022: }
3023: END:
3024: if (s < stockmax) stockmax = s;
3025: if (check)
3026: {
3027: GEN per, alph,pols,p;
3028: if (DEBUGLEVEL) fprintferr("final sort & check...\n");
3029: setlg(norms,s+1); per = sindexsort(norms);
3030: alph = cgetg(s+1,t_VEC);
3031: pols = cgetg(s+1,t_VEC);
3032: for (j=0,i=1; i<=s; i++)
3033: {
3034: long k = per[i];
3035: if (j && mpcmp((GEN)norms[k], borne1) > 0) break;
3036: if ((p = check(data,(GEN)S[k])))
3037: {
3038: if (!j) borne1 = gadd((GEN)norms[k],eps);
3039: j++; pols[j]=(long)p; alph[j]=S[k];
3040: }
3041: }
3042: u = cgetg(3,t_VEC);
3043: setlg(pols,j+1); u[1] = (long)pols;
3044: setlg(alph,j+1); u[2] = (long)alph;
3045: if (isclone(S)) { u[2] = (long)forcecopy(alph); gunclone(S); }
3046: return u;
3047: }
3048: u = cgetg(4,t_VEC);
3049: u[1] = lstoi(s<<1);
3050: u[2] = (long)normax1;
3051: if (stockmax)
3052: {
3053: setlg(S,stockmax+1);
3054: settyp(S,t_MAT);
3055: if (isclone(S)) { p1 = S; S = forcecopy(S); gunclone(p1); }
3056: }
3057: else
3058: S = cgetg(1,t_MAT);
3059: u[3] = (long)S; return u;
3060: }
3061:
1.2 ! noro 3062: /* solve q(x) = x~.a.x <= bound, a > 0.
! 3063: * If check is non-NULL keep x only if check(x).
! 3064: * If noer, return NULL in case of precision problems, raise an error otherwise.
! 3065: * If a is a vector, assume a[1] is the LLL-reduced Cholesky form of q */
! 3066: static GEN
! 3067: _fincke_pohst(GEN a,GEN B0,long stockmax,long noer, long PREC, FP_chk_fun *CHECK)
! 3068: {
! 3069: gpmem_t av = avma;
! 3070: long i,j,l, round = 0;
! 3071: GEN B,r,rinvtrans,u,v,res,z,vnorm,sperm,perm,uperm,gram, bound = B0;
! 3072:
! 3073: if (DEBUGLEVEL>2) fprintferr("entering fincke_pohst\n");
! 3074: if (typ(a) == t_VEC)
! 3075: {
! 3076: r = (GEN)a[1];
! 3077: u = NULL;
1.1 noro 3078: }
3079: else
3080: {
1.2 ! noro 3081: long prec = PREC;
! 3082: l = lg(a);
! 3083: if (l == 1)
! 3084: {
! 3085: if (CHECK) err(talker, "dimension 0 in fincke_pohst");
! 3086: z = cgetg(4,t_VEC);
! 3087: z[1] = z[2] = zero;
! 3088: z[3] = lgetg(1,t_MAT); return z;
! 3089: }
! 3090: i = gprecision(a);
! 3091: if (i) prec = i; else { a = gmul(a,realun(prec)); round = 1; }
! 3092: if (DEBUGLEVEL>2) fprintferr("first LLL: prec = %ld\n", prec);
! 3093: u = lllgramintern(a,4,noer, (prec<<1)-2);
! 3094: if (!u) return NULL;
1.1 noro 3095: r = qf_base_change(a,u,1);
1.2 ! noro 3096: r = sqred1intern(r,noer);
! 3097: if (!r) return NULL;
! 3098: for (i=1; i<l; i++)
! 3099: {
! 3100: GEN s = gsqrt(gcoeff(r,i,i), prec);
! 3101: coeff(r,i,i) = (long)s;
! 3102: for (j=i+1; j<l; j++) coeff(r,i,j) = lmul(s, gcoeff(r,i,j));
! 3103: }
1.1 noro 3104: }
1.2 ! noro 3105: /* now r~ * r = a in LLL basis */
1.1 noro 3106: rinvtrans = gtrans(invmat(r));
1.2 ! noro 3107: if (DEBUGLEVEL>2)
! 3108: fprintferr("final LLL: prec = %ld\n", gprecision(rinvtrans));
! 3109: v = lllintern(rinvtrans, 100, noer, 0);
! 3110: if (!v) return NULL;
! 3111: rinvtrans = gmul(rinvtrans, v);
! 3112:
! 3113: v = ZM_inv(gtrans_i(v),gun);
! 3114: r = gmul(r,v);
! 3115: u = u? gmul(u,v): v;
! 3116:
! 3117: l = lg(r);
! 3118: vnorm = cgetg(l,t_VEC);
! 3119: for (j=1; j<l; j++) vnorm[j] = lnorml2((GEN)rinvtrans[j]);
! 3120: sperm = cgetg(l,t_MAT);
! 3121: uperm = cgetg(l,t_MAT); perm = sindexsort(vnorm);
! 3122: for (i=1; i<l; i++) { uperm[l-i] = u[perm[i]]; sperm[l-i] = r[perm[i]]; }
1.1 noro 3123:
3124: gram = gram_matrix(sperm);
1.2 ! noro 3125: B = gcoeff(gram,l-1,l-1);
! 3126: if (gexpo(B) >= bit_accuracy(lg(B)-2)) return NULL;
1.1 noro 3127:
1.2 ! noro 3128: res = NULL;
! 3129: CATCH(precer) { }
! 3130: TRY {
! 3131: if (CHECK && CHECK->f_init)
! 3132: {
! 3133: bound = CHECK->f_init(CHECK,gram,uperm);
! 3134: if (!bound) err(precer,"fincke_pohst");
! 3135: }
! 3136: if (!bound) bound = gcoeff(gram,1,1);
! 3137:
! 3138: if (DEBUGLEVEL>2) fprintferr("entering smallvectors\n");
! 3139: for (i=1; i<l; i++)
! 3140: {
! 3141: res = smallvectors(gram, bound, stockmax,noer,CHECK);
! 3142: if (!res) err(precer,"fincke_pohst");
! 3143: if (!CHECK || bound || lg(res[2]) > 1) break;
! 3144: if (DEBUGLEVEL>2) fprintferr(" i = %ld failed\n",i);
! 3145: }
! 3146: } ENDCATCH;
! 3147: if (DEBUGLEVEL>2) fprintferr("leaving fincke_pohst\n");
! 3148: if (CHECK || !res) return res;
! 3149:
! 3150: z = cgetg(4,t_VEC);
! 3151: z[1] = lcopy((GEN)res[1]);
! 3152: z[2] = round? lround((GEN)res[2]): lcopy((GEN)res[2]);
! 3153: z[3] = lmul(uperm, (GEN)res[3]); return gerepileupto(av,z);
! 3154: }
1.1 noro 3155:
1.2 ! noro 3156: GEN
! 3157: fincke_pohst(GEN a,GEN B0,long stockmax,long flag, long PREC, FP_chk_fun *CHECK)
! 3158: {
! 3159: gpmem_t av = avma;
! 3160: GEN z = _fincke_pohst(a,B0,stockmax,flag & 1,PREC,CHECK);
! 3161: if (!z)
1.1 noro 3162: {
1.2 ! noro 3163: if (!(flag & 1)) err(precer,"fincke_pohst");
! 3164: avma = av;
1.1 noro 3165: }
1.2 ! noro 3166: return z;
1.1 noro 3167: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to bibli1.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / pari-2.2 / src / basemath