Annotation of OpenXM_contrib/pari/src/basemath/bibli1.c, Revision 1.1
1.1 ! maekawa 1: /********************************************************************/
! 2: /** **/
! 3: /** LLL Algorithm and close friends **/
! 4: /** **/
! 5: /********************************************************************/
! 6: /* $Id: bibli1.c,v 1.4 1999/09/24 10:51:24 karim Exp $ */
! 7: #include "pari.h"
! 8: #include "parinf.h"
! 9: GEN lincomb_integral(GEN u, GEN v, GEN X, GEN Y);
! 10:
! 11: /* scalar product of x with himself */
! 12: static GEN
! 13: sqscal(GEN x)
! 14: {
! 15: long i,av=avma,lx=lg(x);
! 16: GEN z=gzero;
! 17:
! 18: for (i=1; i<lx; i++)
! 19: z = gadd(z, gsqr((GEN)x[i]));
! 20: return gerepileupto(av,z);
! 21: }
! 22:
! 23: /* scalar product x . y */
! 24: static GEN
! 25: gscal(GEN x,GEN y)
! 26: {
! 27: long i,av=avma,lx=lg(x);
! 28: GEN z=gzero;
! 29:
! 30: for (i=1; i<lx; i++)
! 31: z = gadd(z, gmul((GEN)x[i],(GEN)y[i]));
! 32: return gerepileupto(av,z);
! 33: }
! 34:
! 35: static GEN
! 36: sqscali(GEN x)
! 37: {
! 38: GEN z = gzero;
! 39: long i,lx=lg(x);
! 40:
! 41: for (i=1; i<lx; i++)
! 42: z = addii(z,sqri((GEN)x[i]));
! 43: return z;
! 44: }
! 45:
! 46: static GEN
! 47: gscali(GEN x,GEN y)
! 48: {
! 49: GEN z = gzero;
! 50: long i,lx=lg(x);
! 51:
! 52: for (i=1; i<lx; i++)
! 53: z = addii(z, mulii((GEN)x[i],(GEN)y[i]));
! 54: return z;
! 55: }
! 56:
! 57: static GEN
! 58: lllall_trivial(GEN x, long n, long flag)
! 59: {
! 60: GEN y;
! 61: if (!n)
! 62: {
! 63: if (flag != lll_ALL) return cgetg(1,t_MAT);
! 64: y=cgetg(3,t_VEC);
! 65: y[1]=lgetg(1,t_MAT);
! 66: y[2]=lgetg(1,t_MAT); return y;
! 67: }
! 68: /* here n = 1 */
! 69: if (gcmp0((GEN)x[1]))
! 70: {
! 71: switch(flag)
! 72: {
! 73: case lll_KER: return idmat(1);
! 74: case lll_IM : return cgetg(1,t_MAT);
! 75: default: y=cgetg(3,t_VEC);
! 76: y[1]=(long)idmat(1);
! 77: y[2]=lgetg(1,t_MAT); return y;
! 78: }
! 79: }
! 80: if (flag & lll_GRAM) flag ^= lll_GRAM; else x = NULL;
! 81: switch (flag)
! 82: {
! 83: case lll_KER: return cgetg(1,t_MAT);
! 84: case lll_IM : return idmat(1);
! 85: default: y=cgetg(3,t_VEC);
! 86: y[1]=lgetg(1,t_MAT);
! 87: y[2]=x? lcopy(x): (long)idmat(1); return y;
! 88: }
! 89: }
! 90:
! 91: static GEN
! 92: lllgramall_finish(GEN h,GEN fl,long flag,long n)
! 93: {
! 94: long k;
! 95: GEN y;
! 96:
! 97: k=1; while (k<=n && !fl[k]) k++;
! 98: switch(flag)
! 99: {
! 100: case lll_KER: setlg(h,k);
! 101: y = gcopy(h); break;
! 102:
! 103: case lll_IM: h += k-1; h[0] = evaltyp(t_MAT)| evallg(n-k+2);
! 104: y = gcopy(h); break;
! 105:
! 106: default: setlg(h,k); y=cgetg(3,t_VEC);
! 107: y[1] = lcopy(h);
! 108: h += k-1; h[0] = evaltyp(t_MAT)| evallg(n-k+2);
! 109: y[2] = lcopy(h);
! 110: break;
! 111: }
! 112: return y;
! 113: }
! 114:
! 115: /********************************************************************/
! 116: /** **/
! 117: /** LLL with content **/
! 118: /** **/
! 119: /********************************************************************/
! 120:
! 121: /* real Gram matrix has coeffs X[i,j] = x[i,j]*veccon[i]*veccon[j] */
! 122: static GEN
! 123: lllgramintwithcontent(GEN x, GEN veccon, long flag)
! 124: {
! 125: long av0=avma,av,tetpil,lx=lg(x),i,j,k,l,n,lim,kmax;
! 126: GEN u,u2,B,lam,q,r,h,la,bb,p1,p2,p3,p4,fl,corr,corr2,newcon;
! 127: GEN *gptr[8];
! 128:
! 129: if (typ(x) != t_MAT) err(typeer,"lllgramintwithcontent");
! 130: n=lx-1; if (n<=1) return lllall_trivial(x,n,flag);
! 131: if (lg((GEN)x[1])!=lx) err(mattype1,"lllgramintwithcontent");
! 132: fl = new_chunk(lx);
! 133:
! 134: av=avma; lim=stack_lim(av,1);
! 135: x=dummycopy(x); veccon=dummycopy(veccon);
! 136: B=cgetg(lx+1,t_COL); B[1]=un;
! 137: /* B[i+1]=B_i; le vrai B_i est B_i*prod(1,j=1,i,veccon[j]^2) */
! 138:
! 139: for (i=1; i<lx; i++) { B[i+1]=zero; fl[i]=0; }
! 140: lam=cgetg(lx,t_MAT);
! 141: for (j=1; j<lx; j++)
! 142: { p2=cgetg(lx,t_COL); lam[j]=(long)p2; for (i=1; i<lx; i++) p2[i]=zero; }
! 143: /* le vrai lam[i,j] est
! 144: lam[i,j]*veccon[i]*veccon[j]*prod(1,l=1,j-1,veccon[l]^2) */
! 145: k=2; h=idmat(n); kmax=1;
! 146: u=gcoeff(x,1,1); if (typ(u)!=t_INT) err(lllger4);
! 147: if (signe(u)) { B[2]=(long)u; coeff(lam,1,1)=un; fl[1]=1; }
! 148: else { B[2]=un; fl[1]=0; }
! 149: if (DEBUGLEVEL>5) { fprintferr("k = %ld",k); flusherr(); }
! 150: for(;;)
! 151: {
! 152: if (k>kmax)
! 153: {
! 154: kmax=k;
! 155: for (j=1; j<=k; j++)
! 156: {
! 157: if (j==k || fl[j])
! 158: {
! 159: u=gcoeff(x,k,j); if (typ(u)!=t_INT) err(lllger4);
! 160: for (i=1; i<j; i++)
! 161: if (fl[i])
! 162: u=divii(subii(mulii((GEN)B[i+1],u),mulii(gcoeff(lam,k,i),gcoeff(lam,j,i))),(GEN)B[i]);
! 163: if (j<k) coeff(lam,k,j)=(long)u;
! 164: else
! 165: {
! 166: if (signe(u)) { B[k+1]=(long)u; coeff(lam,k,k)=un; fl[k]=1; }
! 167: else { B[k+1]=B[k]; fl[k]=0; }
! 168: }
! 169: }
! 170: }
! 171: if (low_stack(lim, stack_lim(av,1)))
! 172: {
! 173: if(DEBUGMEM>1) err(warnmem,"[1]: lllgramintwithcontent");
! 174: gptr[0]=&B; gptr[1]=&lam; gptr[2]=&h;
! 175: gptr[3]=&x; gptr[4]=&veccon; gerepilemany(av,gptr,5);
! 176: }
! 177: }
! 178: u=shifti(mulii(gcoeff(lam,k,k-1),(GEN)veccon[k]),1);
! 179: u2=mulii((GEN)B[k],(GEN)veccon[k-1]);
! 180: if (cmpii(absi(u),u2)>0)
! 181: {
! 182: q=dvmdii(addii(u,u2),shifti(u2,1),&r);
! 183: if (signe(r)<0) q=addsi(-1,q);
! 184: h[k]=lsub((GEN)h[k],gmul(q,(GEN)h[k-1]));
! 185: newcon=mppgcd((GEN)veccon[k],(GEN)veccon[k-1]);
! 186: corr=divii((GEN)veccon[k],newcon); veccon[k]=(long)newcon;
! 187: if(!gcmp1(corr))
! 188: {
! 189: corr2=sqri(corr);
! 190: for (j=1; j<=n; j++)
! 191: coeff(x,j,k)=coeff(x,k,j)=lmulii(corr,gcoeff(x,k,j));
! 192: coeff(x,k,k)=lmulii(corr,gcoeff(x,k,k));
! 193: for (j=k; j<=kmax; j++) B[j+1]=lmulii(corr2,(GEN)B[j+1]);
! 194: for (i=1; i<k; i++) coeff(lam,k,i)=lmulii(corr,gcoeff(lam,k,i));
! 195: for (i=k+1; i<=kmax; i++)
! 196: {
! 197: coeff(lam,i,k)=lmulii(corr,gcoeff(lam,i,k));
! 198: for (j=k+1; j<i; j++)
! 199: coeff(lam,i,j)=lmulii(corr2,gcoeff(lam,i,j));
! 200: }
! 201: }
! 202: r=negi(mulii(q,divii((GEN)veccon[k-1],(GEN)veccon[k])));
! 203: p1=gcoeff(x,k,k-1);
! 204: for (j=1; j<=n; j++)
! 205: coeff(x,j,k)=coeff(x,k,j)=laddii(gcoeff(x,j,k),mulii(r,(j==k)?p1:gcoeff(x,j,k-1)));
! 206: coeff(x,k,k)=laddii(gcoeff(x,k,k),mulii(r,gcoeff(x,k-1,k)));
! 207: coeff(lam,k,k-1)=laddii(gcoeff(lam,k,k-1),mulii(r,(GEN)B[k]));
! 208: for (i=1; i<k-1; i++)
! 209: coeff(lam,k,i)=laddii(gcoeff(lam,k,i),mulii(r,gcoeff(lam,k-1,i)));
! 210: }
! 211: if (low_stack(lim, stack_lim(av,1)))
! 212: {
! 213: if(DEBUGMEM>1) err(warnmem,"[2]: lllgramintwithcontent");
! 214: gptr[0]=&B; gptr[1]=&lam; gptr[2]=&h;
! 215: gptr[3]=&x; gptr[4]=&veccon; gerepilemany(av,gptr,5);
! 216: }
! 217: p3=mulii((GEN)B[k-1],(GEN)B[k+1]);la=gcoeff(lam,k,k-1);p4=mulii(la,la);
! 218: p2=mulsi(100,mulii(mulii((GEN)veccon[k],(GEN)veccon[k]),addii(p3,p4)));
! 219: p1=mulii((GEN)veccon[k-1],(GEN)B[k]);p1=mulsi(99,mulii(p1,p1));
! 220: if (fl[k-1] && (cmpii(p1,p2)>0 || !fl[k]))
! 221: {
! 222: if (DEBUGLEVEL>=4 && k==n)
! 223: { fprintferr(" (%ld)", expi(p1)-expi(p2)); flusherr(); }
! 224: p1=(GEN)h[k-1]; h[k-1]=h[k]; h[k]=(long)p1;
! 225: p1=(GEN)x[k-1]; x[k-1]=x[k]; x[k]=(long)p1;
! 226: for (j=1; j<=n; j++)
! 227: { p1=gcoeff(x,k-1,j); coeff(x,k-1,j)=coeff(x,k,j); coeff(x,k,j)=(long)p1; }
! 228: p1=(GEN)veccon[k-1]; veccon[k-1]=veccon[k]; veccon[k]=(long)p1;
! 229: for (j=1; j<=k-2; j++)
! 230: { p1=gcoeff(lam,k-1,j); coeff(lam,k-1,j)=coeff(lam,k,j); coeff(lam,k,j)=(long)p1; }
! 231: if (fl[k])
! 232: {
! 233: for (i=k+1; i<=kmax; i++)
! 234: {
! 235: bb=gcoeff(lam,i,k);
! 236: coeff(lam,i,k)=ldivii(subii(mulii((GEN)B[k+1],gcoeff(lam,i,k-1)),mulii(la,bb)),(GEN)B[k]);
! 237: coeff(lam,i,k-1)=ldivii(addii(mulii(la,gcoeff(lam,i,k-1)),mulii((GEN)B[k-1],bb)),(GEN)B[k]);
! 238: if (low_stack(lim, stack_lim(av,1)))
! 239: {
! 240: if(DEBUGMEM>1) err(warnmem,"[3]: lllgramintwithcontent");
! 241: gptr[0]=&B; gptr[1]=&lam; gptr[2]=&h;
! 242: gptr[3]=&x; gptr[4]=&la; gptr[5]=&veccon; gptr[6]=&p3;
! 243: gptr[7]=&p4; gerepilemany(av,gptr,8);
! 244: }
! 245: }
! 246: B[k]=ldivii(addii(p3,p4),(GEN)B[k]);
! 247: }
! 248: else
! 249: {
! 250: if (signe(la))
! 251: {
! 252: p2=(GEN)B[k]; p1=divii(p4,p2);
! 253: for (i=k+1; i<=kmax; i++)
! 254: coeff(lam,i,k-1)=ldivii(mulii(la,gcoeff(lam,i,k-1)),p2);
! 255: for (j=k+1; j<kmax; j++)
! 256: {
! 257: for (i=j+1; i<=kmax; i++)
! 258: coeff(lam,i,j)=ldivii(mulii(p1,gcoeff(lam,i,j)),p2);
! 259: if (low_stack(lim, stack_lim(av,1)))
! 260: {
! 261: if(DEBUGMEM>1) err(warnmem,"[4]: lllgramintwithcontent");
! 262: gptr[0]=&B; gptr[1]=&lam; gptr[2]=&h;
! 263: gptr[3]=&x; gptr[4]=&la; gptr[5]=&veccon; gptr[6]=&p1;
! 264: gptr[7]=&p2; gerepilemany(av,gptr,8);
! 265: }
! 266: }
! 267: B[k+1]=B[k]=(long)p1;
! 268: for (i=k+2; i<=lx; i++)
! 269: B[i]=ldivii(mulii(p1,(GEN)B[i]),p2);
! 270: }
! 271: else
! 272: {
! 273: coeff(lam,k,k-1)=zero;
! 274: for (i=k+1; i<=kmax; i++)
! 275: {
! 276: coeff(lam,i,k)=coeff(lam,i,k-1);
! 277: coeff(lam,i,k-1)=zero;
! 278: }
! 279: B[k]=B[k-1]; fl[k]=1; fl[k-1]=0;
! 280: }
! 281:
! 282: if (low_stack(lim, stack_lim(av,1)))
! 283: {
! 284: if(DEBUGMEM>1) err(warnmem,"[5]: lllgramintwithcontent");
! 285: gptr[0]=&B; gptr[1]=&lam; gptr[2]=&h;
! 286: gptr[3]=&x; gptr[4]=&veccon;
! 287: gerepilemany(av,gptr,5);
! 288: }
! 289: }
! 290: if (k>2) k--;
! 291: if (DEBUGLEVEL>5) { fprintferr(" %ld",k); flusherr(); }
! 292: }
! 293: else
! 294: {
! 295: for (l=k-2; l>=1; l--)
! 296: {
! 297: u=shifti(mulii(gcoeff(lam,k,l),(GEN)veccon[k]),1);
! 298: u2=mulii((GEN)B[l+1],(GEN)veccon[l]);
! 299: if (cmpii(absi(u),u2)>0)
! 300: {
! 301: q=dvmdii(addii(u,u2),shifti(u2,1),&r);
! 302: if (signe(r)<0) q=addsi(-1,q);
! 303: h[k]=lsub((GEN)h[k],gmul(q,(GEN)h[l]));
! 304: newcon=mppgcd((GEN)veccon[k],(GEN)veccon[l]);
! 305: corr=divii((GEN)veccon[k],newcon); veccon[k]=(long)newcon;
! 306: if(!gcmp1(corr))
! 307: {
! 308: corr2=sqri(corr);
! 309: for (j=1; j<=n; j++)
! 310: coeff(x,j,k)=coeff(x,k,j)=lmulii(corr,gcoeff(x,k,j));
! 311: coeff(x,k,k)=lmulii(corr,gcoeff(x,k,k));
! 312: for (j=k; j<=kmax; j++) B[j+1]=lmulii(corr2,(GEN)B[j+1]);
! 313: for (i=1; i<k; i++) coeff(lam,k,i)=lmulii(corr,gcoeff(lam,k,i));
! 314: for (i=k+1; i<=kmax; i++)
! 315: {
! 316: coeff(lam,i,k)=lmulii(corr,gcoeff(lam,i,k));
! 317: for (j=k+1; j<i; j++)
! 318: coeff(lam,i,j)=lmulii(corr2,gcoeff(lam,i,j));
! 319: }
! 320: }
! 321: r=negi(mulii(q,divii((GEN)veccon[l],(GEN)veccon[k])));
! 322: p1=gcoeff(x,k,l);
! 323: for (j=1; j<=n; j++)
! 324: coeff(x,j,k)=coeff(x,k,j)=laddii(gcoeff(x,j,k),mulii(r,(j==k)?p1:gcoeff(x,j,l)));
! 325: coeff(x,k,k)=laddii(gcoeff(x,k,k),mulii(r,gcoeff(x,l,k)));
! 326: coeff(lam,k,l)=laddii(gcoeff(lam,k,l),mulii(r,(GEN)B[l+1]));
! 327: for (i=1; i<l; i++)
! 328: coeff(lam,k,i)=laddii(gcoeff(lam,k,i),mulii(r,gcoeff(lam,l,i)));
! 329: }
! 330: if (low_stack(lim, stack_lim(av,1)))
! 331: {
! 332: if(DEBUGMEM>1) err(warnmem,"[6]: lllgramintwithcontent");
! 333: gptr[0]=&B; gptr[1]=&lam; gptr[2]=&h;
! 334: gptr[3]=&x; gptr[4]=&veccon; gerepilemany(av,gptr,5);
! 335: }
! 336: }
! 337: k++; if (DEBUGLEVEL>5) { fprintferr(" %ld",k); flusherr(); }
! 338: if (k>n)
! 339: {
! 340: if (DEBUGLEVEL>5) { fprintferr("\n"); flusherr(); }
! 341: tetpil=avma;
! 342: return gerepile(av0,tetpil,lllgramall_finish(h,fl,flag,n));
! 343: }
! 344: }
! 345: if (low_stack(lim, stack_lim(av,1)))
! 346: {
! 347: if(DEBUGMEM>1) err(warnmem,"[7]: lllgramintwithcontent");
! 348: gptr[0]=&B; gptr[1]=&lam; gptr[2]=&h;
! 349: gptr[3]=&x; gptr[4]=&veccon; gerepilemany(av,gptr,5);
! 350: }
! 351: }
! 352: }
! 353:
! 354: static GEN
! 355: lllintwithcontent(GEN x)
! 356: {
! 357: long lx=lg(x),i,j,av,tetpil;
! 358: GEN veccon,con,xred,g;
! 359:
! 360: if (typ(x) != t_MAT) err(typeer,"lllintwithcontent");
! 361: if (lx==1) return cgetg(1,t_MAT);
! 362: av=avma; veccon=cgetg(lx,t_VEC); g=cgetg(lx,t_MAT); xred=cgetg(lx,t_MAT);
! 363: for (j=1; j<lx; j++)
! 364: {
! 365: g[j]=lgetg(lx,t_COL); con=content((GEN)x[j]);
! 366: xred[j]=ldiv((GEN)x[j],con); veccon[j]=(long)con;
! 367: }
! 368: for (i=1; i<lx; i++)
! 369: for (j=1; j<=i; j++)
! 370: coeff(g,i,j)=coeff(g,j,i)=(long)gscal((GEN)xred[i],(GEN)xred[j]);
! 371: tetpil=avma;
! 372: return gerepile(av,tetpil,lllgramintwithcontent(g,veccon,2));
! 373: }
! 374:
! 375: /********************************************************************/
! 376: /** **/
! 377: /** LLL ALGORITHM **/
! 378: /** **/
! 379: /********************************************************************/
! 380: #define swap(x,y) { long _t=x; x=y; y=_t; }
! 381: #define gswap(x,y) { GEN _t=x; x=y; y=_t; }
! 382:
! 383: static void
! 384: lllupdate_int(GEN x, GEN h, GEN L, GEN B, long K, long k, long l)
! 385: {
! 386: long i,lx;
! 387: GEN r = truedvmdii(addii(B,shifti(gcoeff(L,k,l),1)), shifti(B,1), NULL);
! 388: GEN *hk,*hl,*xk,*xl;
! 389: if (!signe(r)) return;
! 390: r = negi(r); lx = lg(x);
! 391: hk = (GEN*)h[k]; hl = (GEN*)h[l];
! 392: xk = (GEN*)x[k]; xl = (GEN*)x[l];
! 393: if (is_pm1(r))
! 394: {
! 395: if (signe(r) > 0)
! 396: {
! 397: for (i=1;i<=K;i++) hk[i]=addii(hk[i],hl[i]);
! 398: for (i=1;i<lx;i++) xk[i]=addii(xk[i],xl[i]);
! 399: for (i=1;i<lx;i++) coeff(x,k,i)=laddii(gcoeff(x,k,i),gcoeff(x,l,i));
! 400: for (i=1;i<l; i++) coeff(L,k,i)=laddii(gcoeff(L,k,i),gcoeff(L,l,i));
! 401: } else {
! 402: for (i=1;i<=K;i++) hk[i]=subii(hk[i],hl[i]);
! 403: for (i=1;i<lx;i++) xk[i]=subii(xk[i],xl[i]);
! 404: for (i=1;i<lx;i++) coeff(x,k,i)=lsubii(gcoeff(x,k,i),gcoeff(x,l,i));
! 405: for (i=1;i<l; i++) coeff(L,k,i)=lsubii(gcoeff(L,k,i),gcoeff(L,l,i));
! 406: }
! 407: } else {
! 408: for(i=1;i<=K;i++) hk[i]=addii(hk[i],mulii(r,hl[i]));
! 409: for(i=1;i<lx;i++) xk[i]=addii(xk[i],mulii(r,xl[i]));
! 410: for(i=1;i<lx;i++) coeff(x,k,i)=laddii(gcoeff(x,k,i),mulii(r,gcoeff(x,l,i)));
! 411: for(i=1;i<l;i++) coeff(L,k,i)=laddii(gcoeff(L,k,i),mulii(r,gcoeff(L,l,i)));
! 412: }
! 413: coeff(L,k,l)=laddii(gcoeff(L,k,l),mulii(r,B));
! 414: }
! 415:
! 416: static void
! 417: lllupdate(GEN x, GEN h, GEN L, long K, long k, long l)
! 418: {
! 419: long e,i,lx;
! 420: GEN r = grndtoi(gcoeff(L,k,l),&e);
! 421: GEN *hk,*hl,*xk,*xl;
! 422: if (DEBUGLEVEL>8)
! 423: fprintferr("error bits when rounding in lllgram: %ld\n",e);
! 424: if (!signe(r)) return;
! 425: r = negi(r); lx = lg(x);
! 426: hk = (GEN*)h[k]; hl = (GEN*)h[l];
! 427: xk = (GEN*)x[k]; xl = (GEN*)x[l];
! 428: if (is_pm1(r))
! 429: {
! 430: if (signe(r) > 0)
! 431: {
! 432: for (i=1;i<=K;i++) hk[i]=addii(hk[i],hl[i]);
! 433: for (i=1;i<lx;i++) xk[i]=gadd(xk[i],xl[i]);
! 434: for (i=1;i<lx;i++) coeff(x,k,i)=ladd(gcoeff(x,k,i),gcoeff(x,l,i));
! 435: for (i=1;i<l; i++) coeff(L,k,i)=ladd(gcoeff(L,k,i),gcoeff(L,l,i));
! 436: } else {
! 437: for (i=1;i<=K;i++) hk[i]=subii(hk[i],hl[i]);
! 438: for (i=1;i<lx;i++) xk[i]=gsub(xk[i],xl[i]);
! 439: for (i=1;i<lx;i++) coeff(x,k,i)=lsub(gcoeff(x,k,i),gcoeff(x,l,i));
! 440: for (i=1;i<l; i++) coeff(L,k,i)=lsub(gcoeff(L,k,i),gcoeff(L,l,i));
! 441: }
! 442: } else {
! 443: for (i=1;i<=K;i++) hk[i]=addii(hk[i],mulii(r,hl[i]));
! 444: for (i=1;i<lx;i++) xk[i]=gadd(xk[i],gmul(r,xl[i]));
! 445: for (i=1;i<lx;i++) coeff(x,k,i)=ladd(gcoeff(x,k,i),gmul(r,gcoeff(x,l,i)));
! 446: for (i=1;i<l; i++) coeff(L,k,i)=ladd(gcoeff(L,k,i),gmul(r,gcoeff(L,l,i)));
! 447: }
! 448: coeff(L,k,l)=ladd(gcoeff(L,k,l),r);
! 449: }
! 450:
! 451: /* x integer matrix */
! 452: GEN
! 453: lllgramall(GEN x, long alpha, long flag)
! 454: {
! 455: long av0=avma,av,tetpil,lim,lx=lg(x),i,j,k,l,n,kmax;
! 456: GEN u,B,L,h,la,p1,p2,p3,p4,fl, *gptr[6];
! 457:
! 458: if (typ(x) != t_MAT) err(typeer,"lllgramall");
! 459: n=lx-1; if (n<=1) return lllall_trivial(x,n,flag);
! 460: if (lg((GEN)x[1])!=lx) err(mattype1,"lllgramall");
! 461: fl = new_chunk(lx);
! 462:
! 463: av=avma; lim=stack_lim(av,1); x=dummycopy(x);
! 464: B=gscalcol(gun, lx);
! 465: L=cgetg(lx,t_MAT);
! 466: for (j=1; j<lx; j++)
! 467: {
! 468: for (i=1; i<lx; i++)
! 469: if (typ(gcoeff(x,i,j))!=t_INT) err(lllger4);
! 470: fl[j] = 0; L[j] = (long)zerocol(n);
! 471: }
! 472: k=2; h=idmat(n); kmax=1;
! 473: u=gcoeff(x,1,1);
! 474: if (signe(u)) { B[2]=(long)u; coeff(L,1,1)=un; fl[1]=1; } else B[2]=un;
! 475: if (DEBUGLEVEL>5) fprintferr("k =");
! 476: for(;;)
! 477: {
! 478: if (k > kmax)
! 479: {
! 480: if (DEBUGLEVEL>3) {fprintferr(" K%ld",k);flusherr();}
! 481: kmax = k;
! 482: for (j=1; j<=k; j++)
! 483: if (j==k || fl[j])
! 484: {
! 485: long av1 = avma;
! 486: u = gcoeff(x,k,j);
! 487: for (i=1; i<j; i++) if (fl[i])
! 488: u = divii(subii(mulii((GEN)B[i+1],u),
! 489: mulii(gcoeff(L,k,i),gcoeff(L,j,i))),
! 490: (GEN)B[i]);
! 491: u = gerepileuptoint(av1, u);
! 492: if (j<k) coeff(L,k,j)=(long)u;
! 493: else
! 494: {
! 495: if (signe(u)) { B[k+1]=(long)u; coeff(L,k,k)=un; fl[k]=1; }
! 496: else { B[k+1]=B[k]; fl[k]=0; }
! 497: }
! 498: }
! 499: } else if (DEBUGLEVEL>5) {fprintferr(" %ld",k); flusherr();}
! 500: lllupdate_int(x,h,L,(GEN)B[k],kmax,k,k-1);
! 501: if (fl[k-1] &&
! 502: (cmpii(mulsi(alpha-1,sqri((GEN)B[k])),
! 503: mulsi(alpha,p3=addii(mulii((GEN)B[k-1],(GEN)B[k+1]),
! 504: p4=sqri(la=gcoeff(L,k,k-1))))) > 0
! 505: || !fl[k]))
! 506: {
! 507: if (DEBUGLEVEL>3 && k==n)
! 508: {
! 509: fprintferr(" (%ld)", expi(mulsi(alpha-1,sqri((GEN)B[k])))
! 510: - expi(mulsi(alpha,p3)));
! 511: flusherr();
! 512: }
! 513: swap(h[k-1], h[k]);
! 514: swap(x[k-1], x[k]);
! 515: for (j=1; j<=n; j++) swap(coeff(x,k-1,j), coeff(x,k,j));
! 516: for (j=1; j<k-1; j++) swap(coeff(L,k-1,j), coeff(L,k,j))
! 517: if (fl[k])
! 518: {
! 519: GEN Bk = (GEN)B[k];
! 520: long av1 = avma;
! 521: for (i=k+1; i<=kmax; i++)
! 522: {
! 523: GEN bb=gcoeff(L,i,k);
! 524: p1 = divii(subii(mulii((GEN)B[k+1],gcoeff(L,i,k-1)),
! 525: mulii(la,bb)), Bk);
! 526: av1=avma=coeff(L,i,k) = (long)icopy_av(p1,(GEN)av1);
! 527: p1 = divii(addii(mulii(la,gcoeff(L,i,k-1)),
! 528: mulii((GEN)B[k-1],bb)), Bk);
! 529: av1=avma=coeff(L,i,k-1) = (long)icopy_av(p1,(GEN)av1);
! 530: }
! 531: B[k] = ldivii(p3,Bk);
! 532: }
! 533: else
! 534: {
! 535: if (signe(la))
! 536: {
! 537: p2=(GEN)B[k]; p1=divii(p4,p2);
! 538: B[k+1]=B[k]=(long)p1;
! 539: for (i=k+2; i<=lx; i++)
! 540: B[i]=ldivii(mulii(p1,(GEN)B[i]),p2);
! 541: for (i=k+1; i<=kmax; i++)
! 542: coeff(L,i,k-1)=ldivii(mulii(la,gcoeff(L,i,k-1)),p2);
! 543: for (j=k+1; j<kmax; j++)
! 544: {
! 545: for (i=j+1; i<=kmax; i++)
! 546: coeff(L,i,j)=ldivii(mulii(p1,gcoeff(L,i,j)),p2);
! 547: if (low_stack(lim, stack_lim(av,1)))
! 548: {
! 549: if(DEBUGMEM>1) err(warnmem,"lllgramall[1]");
! 550: gptr[0]=&B; gptr[1]=&L; gptr[2]=&h;
! 551: gptr[3]=&x; gptr[4]=&p2;
! 552: gerepilemany(av,gptr,5); p1=(GEN)B[k];
! 553: }
! 554: }
! 555: }
! 556: else
! 557: {
! 558: for (i=k+1; i<=kmax; i++)
! 559: { coeff(L,i,k)=coeff(L,i,k-1); coeff(L,i,k-1)=zero; }
! 560: B[k]=B[k-1]; fl[k]=1; fl[k-1]=0;
! 561: }
! 562: }
! 563: if (k>2) k--;
! 564: }
! 565: else
! 566: {
! 567: for (l=k-2; l; l--)
! 568: {
! 569: lllupdate_int(x,h,L,(GEN)B[l+1],kmax,k,l);
! 570: if (low_stack(lim, stack_lim(av,1)))
! 571: {
! 572: if(DEBUGMEM>1) err(warnmem,"lllgramall[2]");
! 573: gptr[0]=&B; gptr[1]=&L; gptr[2]=&h;
! 574: gptr[3]=&x; gerepilemany(av,gptr,4);
! 575: }
! 576: }
! 577: if (++k > n) break;
! 578: }
! 579: if (low_stack(lim, stack_lim(av,1)))
! 580: {
! 581: if(DEBUGMEM>1) err(warnmem,"lllgramall[3]");
! 582: gptr[0]=&B; gptr[1]=&L; gptr[2]=&h;
! 583: gptr[3]=&x; gerepilemany(av,gptr,4);
! 584: }
! 585: }
! 586: if (DEBUGLEVEL>3) fprintferr("\n");
! 587: tetpil=avma; return gerepile(av0,tetpil,lllgramall_finish(h,fl,flag,n));
! 588: }
! 589:
! 590: static GEN
! 591: lllgramall0(GEN x, long flag) { return lllgramall(x,100,flag); }
! 592:
! 593: static int
! 594: pslg(GEN x)
! 595: {
! 596: long tx;
! 597: if (gcmp0(x)) return 2;
! 598: tx=typ(x); return is_scalar_t(tx)? 3: lgef(x);
! 599: }
! 600:
! 601: static GEN
! 602: lllgramallgen(GEN x, long flag)
! 603: {
! 604: long av0=avma,av,tetpil,lx=lg(x),tu,i,j,k,l,n,lim;
! 605: GEN u,B,lam,q,cq,h,la,bb,p1,p2,p3,p4,fl;
! 606: int ps1,ps2,flc;
! 607:
! 608: if (typ(x) != t_MAT) err(typeer,"lllgramallgen");
! 609: n=lx-1; if (n<=1) return lllall_trivial(x,n,flag);
! 610: if (lg((GEN)x[1])!=lx) err(mattype1,"lllgramallgen");
! 611:
! 612: fl = new_chunk(lx);
! 613:
! 614: av=avma; lim=stack_lim(av,1);
! 615: B=cgetg(lx+1,t_COL);
! 616: B[1]=un; lam=cgetg(lx,t_MAT);
! 617: for (j=1; j<lx; j++) lam[j]=lgetg(lx,t_COL);
! 618: for (i=1; i<lx; i++)
! 619: for (j=1; j<=i; j++)
! 620: {
! 621: if (j<i && !fl[j]) coeff(lam,i,j)=coeff(lam,j,i)=zero;
! 622: else
! 623: {
! 624: u=gcoeff(x,i,j); tu=typ(u);
! 625: if (! is_scalar_t(tu) && tu != t_POL) err(lllger4);
! 626: for (k=1; k<j; k++)
! 627: if (fl[k])
! 628: u=gdiv(gsub( gmul((GEN)B[k+1],u),
! 629: gmul(gcoeff(lam,i,k),gcoeff(lam,j,k)) ),
! 630: (GEN)B[k]);
! 631: if (j<i) { coeff(lam,i,j)=(long)u; coeff(lam,j,i)=zero; }
! 632: else
! 633: {
! 634: if (!gcmp0(u)) { B[i+1]=(long)u; coeff(lam,i,i)=un; fl[i]=1; }
! 635: else { B[i+1]=B[i]; coeff(lam,i,i)=zero; fl[i]=0; }
! 636: }
! 637: }
! 638: }
! 639: k=2; h=idmat(n); flc=0;
! 640: for(;;)
! 641: {
! 642: u=gcoeff(lam,k,k-1);
! 643: if (pslg(u) >= pslg((GEN)B[k]))
! 644: {
! 645: q=gdeuc(u,(GEN)B[k]); cq=gdivsg(1,content(q)); q=gmul(q,cq); flc=1;
! 646: h[k]=lsub(gmul(cq,(GEN)h[k]),gmul(q,(GEN)h[k-1]));
! 647: coeff(lam,k,k-1)=lsub(gmul(cq,gcoeff(lam,k,k-1)),gmul(q,(GEN)B[k]));
! 648: for (i=1; i<k-1; i++)
! 649: coeff(lam,k,i)=lsub(gmul(cq,gcoeff(lam,k,i)),gmul(q,gcoeff(lam,k-1,i)));
! 650: }
! 651: ps1 = pslg(gsqr((GEN)B[k]));
! 652: p3 = gmul((GEN)B[k-1],(GEN)B[k+1]);
! 653: la=gcoeff(lam,k,k-1); p4 = gmul(la,gcoeff(lam,k,k-1));
! 654: ps2=pslg(gadd(p3,p4));
! 655: if (fl[k-1] && (ps1>ps2 || (ps1==ps2 && flc) || !fl[k]))
! 656: {
! 657: flc = (ps1!=ps2);
! 658: swap(h[k-1],h[k]);
! 659: for (j=1; j<=k-2; j++) swap(coeff(lam,k-1,j), coeff(lam,k,j));
! 660: if (fl[k])
! 661: {
! 662: for (i=k+1; i<=n; i++)
! 663: {
! 664: bb=gcoeff(lam,i,k);
! 665: coeff(lam,i,k)=ldiv(gsub(gmul((GEN)B[k+1],gcoeff(lam,i,k-1)),gmul(la,bb)),(GEN)B[k]);
! 666: coeff(lam,i,k-1)=ldiv(gadd(gmul(la,gcoeff(lam,i,k-1)),gmul((GEN)B[k-1],bb)),(GEN)B[k]);
! 667: }
! 668: B[k]=ldiv(gadd(p3,p4),(GEN)B[k]);
! 669: }
! 670: else
! 671: {
! 672: if (!gcmp0(la))
! 673: {
! 674: p2=(GEN)B[k]; p1=gdiv(p4,p2);
! 675: for (i=k+1; i<lx; i++)
! 676: coeff(lam,i,k-1)=ldiv(gmul(la,gcoeff(lam,i,k-1)),p2);
! 677: for (j=k+1; j<lx-1; j++)
! 678: for (i=j+1; i<lx; i++)
! 679: coeff(lam,i,j)=ldiv(gmul(p1,gcoeff(lam,i,j)),p2);
! 680: B[k+1]=B[k]=(long)p1;
! 681: for (i=k+2; i<=lx; i++)
! 682: B[i]=ldiv(gmul(p1,(GEN)B[i]),p2);
! 683: }
! 684: else
! 685: {
! 686: coeff(lam,k,k-1)=zero;
! 687: for (i=k+1; i<lx; i++)
! 688: { coeff(lam,i,k)=coeff(lam,i,k-1); coeff(lam,i,k-1)=zero; }
! 689: B[k]=B[k-1]; fl[k]=1; fl[k-1]=0;
! 690: }
! 691: }
! 692: if (k>2) k--;
! 693: }
! 694: else
! 695: {
! 696: for (l=k-2; l>=1; l--)
! 697: {
! 698: u=gcoeff(lam,k,l);
! 699: if (pslg(u)>=pslg((GEN)B[l+1]))
! 700: {
! 701: q=gdeuc(u,(GEN)B[l+1]); cq=gdivsg(1,content(q));
! 702: q=gmul(q,cq); flc=1;
! 703: h[k]=lsub(gmul(cq,(GEN)h[k]),gmul(q,(GEN)h[l]));
! 704: coeff(lam,k,l)=lsub(gmul(cq,gcoeff(lam,k,l)),gmul(q,(GEN)B[l+1]));
! 705: for (i=1; i<l; i++)
! 706: coeff(lam,k,i)=lsub(gmul(cq,gcoeff(lam,k,i)),gmul(q,gcoeff(lam,l,i)));
! 707: }
! 708: }
! 709: if (++k > n) break;
! 710: }
! 711: if (low_stack(lim, stack_lim(av,1)))
! 712: {
! 713: GEN *gptr[4];
! 714: if(DEBUGMEM>1) err(warnmem,"lllgramallgen");
! 715: gptr[0]=&B; gptr[1]=&lam; gptr[2]=&h;
! 716: gerepilemany(av,gptr,3);
! 717: }
! 718: }
! 719: tetpil=avma;
! 720: return gerepile(av0,tetpil,lllgramall_finish(h,fl,flag,n));
! 721: }
! 722:
! 723: /* return x[k,j] - (mu.A)[j] */
! 724: #ifdef INLINE
! 725: INLINE
! 726: #endif
! 727: GEN
! 728: get_Aj(GEN x, GEN mu, GEN A, long j, long k)
! 729: {
! 730: long av,i;
! 731: GEN s;
! 732:
! 733: if (j==1) return gcopy(gcoeff(x,k,1));
! 734: av = avma; s = gmul(gcoeff(mu,j,1),(GEN)A[1]);
! 735: for (i=2; i<j; i++) s=gadd(s,gmul(gcoeff(mu,j,i),(GEN)A[i]));
! 736: s=gneg(s); return gerepileupto(av, gadd(gcoeff(x,k,j),s));
! 737: }
! 738:
! 739: /* x = Gram(b_i). If precision problems return NULL if flag=1, error otherwise.
! 740: * Quality ratio = (alpha-1)/alpha. Suggested value: alpha = 100
! 741: */
! 742: GEN
! 743: lllgramintern(GEN x, long alpha, long flag, long prec)
! 744: {
! 745: GEN xinit,L,h,A,B,L1,L2,q,cst;
! 746: long retry = 2, av = avma,tetpil,lim,l,i,j,k,k1,lx=lg(x),n,kmax;
! 747:
! 748: if (typ(x) != t_MAT) err(typeer,"lllgram");
! 749: n=lx-1; if (n && lg((GEN)x[1])!=lx) err(mattype1,"lllgram");
! 750: if (n<=1) return idmat(n);
! 751: lim = stack_lim(av,1); xinit=x;
! 752: for (k=2,j=1; j<lx; j++)
! 753: {
! 754: GEN p1=(GEN)x[j];
! 755: for (i=1; i<lx; i++) /* FIXME: lg <-> expo */
! 756: if (typ(p1[i]) == t_REAL) { l = lg((GEN)p1[i]); if (l>k) k=l; }
! 757: }
! 758: if (k == 2)
! 759: {
! 760: if (!prec) return lllgramint(x);
! 761: x = gmul(x, realun(prec+1));
! 762: }
! 763: else if (prec < k) prec = k;
! 764: h = idmat(n);
! 765: x = gprec_w(x, prec+1);
! 766:
! 767: LABLLLGRAM:
! 768: switch(retry--)
! 769: {
! 770: case 2: /* entry */ break;
! 771: case 1: /* failed already */
! 772: tetpil = avma; h = gcopy(h);
! 773: prec = (prec<<1)-2;
! 774: if (DEBUGLEVEL > 3) fprintferr("\n");
! 775: if (DEBUGLEVEL) err(warnprec,"lllgramintern",prec);
! 776: x = qf_base_change(gprec_w(xinit,prec),h,1);
! 777: {
! 778: GEN *gsav[2]; gsav[0]=&h; gsav[1]=&x;
! 779: gerepilemanysp(av, tetpil, gsav, 2);
! 780: }
! 781: if (DEBUGLEVEL) err(warner,"lllgramintern starting over");
! 782: break;
! 783: case 0: /* give up */
! 784: if (DEBUGLEVEL > 3) fprintferr("\n");
! 785: if (DEBUGLEVEL) err(warner,"lllgramintern giving up");
! 786: if (flag) { avma=av; return NULL; }
! 787: if (DEBUGLEVEL) outerr(xinit);
! 788: err(lllger3);
! 789: }
! 790: cst = cgetr(prec+1); affsr(alpha-1,cst);
! 791: cst = divrs(cst,alpha);
! 792:
! 793: L=cgetg(lx,t_MAT);
! 794: B=cgetg(lx,t_COL);
! 795: A=cgetg(lx,t_VEC);
! 796: for (j=1; j<lx; j++)
! 797: {
! 798: L[j] = (long)zerocol(n);
! 799: A[j] = B[j] = zero;
! 800: }
! 801: k=2; kmax=1; B[1]=coeff(x,1,1);
! 802: if (gcmp0((GEN)B[1]))
! 803: {
! 804: if (flag) return NULL;
! 805: err(lllger3);
! 806: }
! 807: if (DEBUGLEVEL>5) fprintferr("k =");
! 808: for(;;)
! 809: {
! 810: if (k>kmax)
! 811: {
! 812: if (DEBUGLEVEL>3) {fprintferr(" K%ld",k);flusherr();}
! 813: kmax=k;
! 814: for (j=1; j<k; j++)
! 815: {
! 816: A[j] = (long)get_Aj(x,L,A,j,k);
! 817: coeff(L,k,j) = ldiv((GEN)A[j],(GEN)B[j]);
! 818: }
! 819: B[k] = (long)get_Aj(x,L,A,k,k);
! 820: if (gsigne((GEN)B[k]) <= 0)
! 821: {
! 822: if (kmax == 2) retry = 0;
! 823: goto LABLLLGRAM;
! 824: }
! 825: }
! 826: else if (DEBUGLEVEL>5) fprintferr(" %ld",k);
! 827: L1=gcoeff(L,k,k-1);
! 828: if (DEBUGLEVEL>9)
! 829: {
! 830: fprintferr(" %ld", gexpo(L1) - bit_accuracy(lg(L1)));
! 831: for (i=1; i<lx; i++)
! 832: fprintferr("%ld: %Z\n",i,qfeval(x,(GEN)h[i]));
! 833: flusherr();
! 834: }
! 835: if (2*gexpo(L1) > bit_accuracy(lg(L1)))
! 836: {
! 837: if (DEBUGLEVEL>3)
! 838: {
! 839: fprintferr("\nRecomputing Gram-Schmidt, kmax = %ld\n",kmax);
! 840: if (DEBUGLEVEL>9) fprintferr("Old B = %Z\n",B);
! 841: }
! 842: for (k1=1; k1<=kmax; k1++)
! 843: {
! 844: for (j=1; j<k1; j++)
! 845: {
! 846: A[j] = (long)get_Aj(x,L,A,j,k1);
! 847: coeff(L,k1,j) = ldiv((GEN)A[j],(GEN)B[j]);
! 848: }
! 849: B[k1] = (long)get_Aj(x,L,A,k1,k1);
! 850: if (gsigne((GEN)B[k1]) <= 0) goto LABLLLGRAM;
! 851: }
! 852: if (DEBUGLEVEL>9) fprintferr("New B = %Z\n",B);
! 853: }
! 854: lllupdate(x,h,L,kmax,k,k-1);
! 855: L1 = gcoeff(L,k,k-1);
! 856: L2 = gsqr(L1);
! 857: q = gmul((GEN)B[k-1], gsub(cst,L2));
! 858: if (gcmp(q,(GEN)B[k]) > 0)
! 859: {
! 860: GEN BK,BB;
! 861: if (DEBUGLEVEL>3 && k==kmax)
! 862: { fprintferr(" (%ld)",gexpo(q)-gexpo((GEN)B[k])); flusherr(); }
! 863: BB = gadd((GEN)B[k], gmul((GEN)B[k-1],L2));
! 864: if (gcmp0(BB)) goto LABLLLGRAM;
! 865: coeff(L,k,k-1) = ldiv(gmul(L1,(GEN)B[k-1]), BB);
! 866: BK = gdiv((GEN)B[k],BB);
! 867: B[k] = lmul((GEN)B[k-1], BK);
! 868: B[k-1] = (long)BB;
! 869: swap(h[k-1],h[k]);
! 870: swap(x[k-1],x[k]);
! 871: for (j=1; j<=n; j++) swap(coeff(x,k-1,j), coeff(x,k,j));
! 872: for (j=1; j<k-1; j++) swap(coeff(L,k-1,j), coeff(L,k,j))
! 873: for (i=k+1; i<=kmax; i++)
! 874: {
! 875: GEN p=gcoeff(L,i,k);
! 876: coeff(L,i,k) = lsub(gcoeff(L,i,k-1),gmul(L1,p));
! 877: coeff(L,i,k-1)=ladd(gmul(BK,p), gmul(gcoeff(L,k,k-1),gcoeff(L,i,k-1)));
! 878: }
! 879: if (k>2) k--;
! 880: }
! 881: else
! 882: {
! 883: for (l=k-2; l; l--) lllupdate(x,h,L,kmax,k,l);
! 884: if (++k > n) break;
! 885: }
! 886: if (low_stack(lim, stack_lim(av,1)))
! 887: {
! 888: GEN *gptr[6];
! 889: if(DEBUGMEM>1)
! 890: {
! 891: if (DEBUGLEVEL > 3) fprintferr("\n");
! 892: err(warnmem,"lllgram");
! 893: }
! 894: gptr[0]=&B; gptr[1]=&L; gptr[2]=&h; gptr[3]=&A;
! 895: gptr[4]=&x; gptr[5]=&cst; gerepilemany(av,gptr,6);
! 896: }
! 897: }
! 898: if (DEBUGLEVEL>3) fprintferr("\n");
! 899: tetpil=avma; return gerepile(av,tetpil,gcopy(h));
! 900: }
! 901:
! 902: static GEN
! 903: lllgram_noerr(GEN x,long prec) { return lllgramintern(x,100,1,prec); }
! 904:
! 905: GEN
! 906: lllgram(GEN x,long prec) { return lllgramintern(x,100,0,prec); }
! 907:
! 908: GEN
! 909: qflll0(GEN x, long flag, long prec)
! 910: {
! 911: switch(flag)
! 912: {
! 913: case 0: return lll(x,prec);
! 914: case 1: return lllint(x);
! 915: case 2: return lllintpartial(x);
! 916: case 3: return lllrat(x);
! 917: case 4: return lllkerim(x);
! 918: case 5: return lllkerimgen(x);
! 919: case 7: return lll1(x,prec);
! 920: case 8: return lllgen(x);
! 921: case 9: return lllintwithcontent(x);
! 922: default: err(flagerr,"qflll");
! 923: }
! 924: return NULL; /* not reached */
! 925: }
! 926:
! 927: GEN
! 928: qflllgram0(GEN x, long flag, long prec)
! 929: {
! 930: switch(flag)
! 931: {
! 932: case 0: return lllgram(x,prec);
! 933: case 1: return lllgramint(x);
! 934: case 4: return lllgramkerim(x);
! 935: case 5: return lllgramkerimgen(x);
! 936: case 7: return lllgram1(x,prec);
! 937: case 8: return lllgramgen(x);
! 938: default: err(flagerr,"qflllgram");
! 939: }
! 940: return NULL; /* not reached */
! 941: }
! 942:
! 943: /* x est la matrice d'une base b_i; retourne la matrice u (entiere
! 944: * unimodulaire) d'une base LLL-reduite c_i en fonction des b_i (la base
! 945: * reduite est c=b*u).
! 946: */
! 947: static GEN
! 948: lll_proto(GEN x, GEN f(GEN, long), long prec)
! 949: {
! 950: long lx=lg(x),i,j,av,av1;
! 951: GEN g;
! 952:
! 953: if (typ(x) != t_MAT) err(typeer,"lll_proto");
! 954: if (lx==1) return cgetg(1,t_MAT);
! 955: av=avma; g=cgetg(lx,t_MAT);
! 956: for (j=1; j<lx; j++) g[j]=lgetg(lx,t_COL);
! 957: for (i=1; i<lx; i++)
! 958: for (j=1; j<=i; j++)
! 959: coeff(g,i,j)=coeff(g,j,i)=(long)gscal((GEN)x[i],(GEN)x[j]);
! 960: av1=avma; x = f(g,prec);
! 961: if (!x) { avma=av; return NULL; }
! 962: return gerepile(av,av1,x);
! 963: }
! 964:
! 965: GEN
! 966: lllintern(GEN x,long flag,long prec)
! 967: {
! 968: return lll_proto(x,flag? lllgram_noerr: lllgram,prec);
! 969: }
! 970:
! 971: GEN
! 972: lll(GEN x,long prec) { return lll_proto(x,lllgram,prec); }
! 973:
! 974: GEN
! 975: lll1(GEN x,long prec) { return lll_proto(x,lllgram1,prec); }
! 976:
! 977: GEN
! 978: lllrat(GEN x) { return lll_proto(x,lllgram,lll_ALL); }
! 979:
! 980: GEN
! 981: lllint(GEN x) { return lll_proto(x,lllgramall0,lll_IM); }
! 982:
! 983: GEN
! 984: lllgen(GEN x) { return lll_proto(x,lllgramallgen,lll_IM); }
! 985:
! 986: GEN
! 987: lllgramgen(GEN x) { return lllgramallgen(x,lll_IM); }
! 988:
! 989: GEN
! 990: lllgramkerimgen(GEN x) { return lllgramallgen(x,lll_ALL); }
! 991:
! 992: static GEN
! 993: lllkerim_proto(GEN x, GEN f(GEN,long))
! 994: {
! 995: long lx=lg(x), i,j,av,av1;
! 996: GEN g;
! 997:
! 998: if (typ(x) != t_MAT) err(typeer,"lllkerim_proto");
! 999: if (lx==1)
! 1000: {
! 1001: g=cgetg(3,t_VEC);
! 1002: g[1]=lgetg(1,t_MAT);
! 1003: g[2]=lgetg(1,t_MAT); return g;
! 1004: }
! 1005: if (lg((GEN)x[1])==1)
! 1006: {
! 1007: g=cgetg(3,t_VEC);
! 1008: g[1]=(long)idmat(lx-1);
! 1009: g[2]=lgetg(1,t_MAT); return g;
! 1010: }
! 1011: av=avma; g=cgetg(lx,t_MAT);
! 1012: for (j=1; j<lx; j++) g[j]=lgetg(lx,t_COL);
! 1013: for (i=1; i<lx; i++)
! 1014: for (j=1; j<=i; j++)
! 1015: coeff(g,i,j)=coeff(g,j,i)=(long)gscal((GEN)x[i],(GEN)x[j]);
! 1016: av1=avma; return gerepile(av,av1,f(g,lll_ALL));
! 1017: }
! 1018:
! 1019: GEN
! 1020: lllkerim(GEN x) { return lllkerim_proto(x,lllgramall0); }
! 1021:
! 1022: GEN
! 1023: lllkerimgen(GEN x) { return lllkerim_proto(x,lllgramallgen); }
! 1024:
! 1025: /* x est ici la matrice de GRAM des bi. */
! 1026: GEN
! 1027: lllgram1(GEN x, long prec)
! 1028: {
! 1029: GEN mu,u,B,BB,BK,p,q,r,cst,unreel,sv,mu1,mu2;
! 1030: long av,tetpil,lim,l,i,j,k,lx=lg(x),n,e;
! 1031:
! 1032: if (typ(x) != t_MAT) err(typeer,"lllgram1");
! 1033: if (lg((GEN)x[1])!=lx) err(mattype1,"lllgram1"); n=lx-1;
! 1034: if (n<=1) return idmat(n);
! 1035: cst=gdivgs(stoi(99),100); /* LLL proposent 0.75 */
! 1036: if (prec)
! 1037: {
! 1038: unreel = realun(prec+1);
! 1039: x = gmul(x,unreel);
! 1040: cst = gmul(cst,unreel);
! 1041: }
! 1042: av=avma; lim=stack_lim(av,1);
! 1043: mu=gtrans(sqred(x)); B=cgetg(lx,t_COL);
! 1044: for (i=1,l=0; i<=n; i++)
! 1045: {
! 1046: if (gsigne((GEN)(B[i]=coeff(mu,i,i)))>0) l++;
! 1047: coeff(mu,i,i)=un;
! 1048: }
! 1049: if (l<n) err(lllger3);
! 1050:
! 1051: u=idmat(n); k=2;
! 1052: do
! 1053: {
! 1054: if (!gcmp0(r=grndtoi(gcoeff(mu,k,k-1),&e)))
! 1055: {
! 1056: u[k]=lsub((GEN)u[k],gmul(r,(GEN)u[k-1]));
! 1057: for (j=1; j<k-1; j++)
! 1058: coeff(mu,k,j)=lsub(gcoeff(mu,k,j),gmul(r,gcoeff(mu,k-1,j)));
! 1059: mu1=(GEN)(coeff(mu,k,k-1)=lsub(gcoeff(mu,k,k-1),r));
! 1060: }
! 1061: else mu1=gcoeff(mu,k,k-1);
! 1062: q=gmul((GEN)B[k-1],gsub(cst,mu2=gsqr(mu1)));
! 1063: if (gcmp(q,(GEN)B[k])>0)
! 1064: {
! 1065: BB=gadd((GEN)B[k],gmul((GEN)B[k-1],mu2));
! 1066: coeff(mu,k,k-1)=ldiv(gmul(mu1,(GEN)B[k-1]),BB);
! 1067: B[k]=lmul((GEN)B[k-1],BK=gdiv((GEN)B[k],BB));
! 1068: B[k-1]=(long)BB;
! 1069: swap(u[k-1],u[k]);
! 1070: for (j=1; j<=k-2; j++) swap(coeff(mu,k-1,j), coeff(mu,k,j));
! 1071: for (i=k+1; i<=n; i++)
! 1072: {
! 1073: p=gcoeff(mu,i,k);
! 1074: coeff(mu,i,k)=lsub(gcoeff(mu,i,k-1),gmul(mu1,p));
! 1075: coeff(mu,i,k-1)=ladd(gmul(BK,p),gmul(gcoeff(mu,k,k-1),gcoeff(mu,i,k-1)));
! 1076: }
! 1077: if (k>2) k--;
! 1078: }
! 1079: else
! 1080: {
! 1081: for (l=k-2; l; l--)
! 1082: {
! 1083: if (!gcmp0(r=grndtoi(gcoeff(mu,k,l),&e)))
! 1084: {
! 1085: u[k]=lsub((GEN)u[k],gmul(r,(GEN)u[l]));
! 1086: for (j=1; j<l; j++)
! 1087: coeff(mu,k,j)=lsub(gcoeff(mu,k,j),gmul(r,gcoeff(mu,l,j)));
! 1088: coeff(mu,k,l)=lsub(gcoeff(mu,k,l),r);
! 1089: }
! 1090: }
! 1091: k++;
! 1092: }
! 1093: if (low_stack(lim, stack_lim(av,1)))
! 1094: {
! 1095: if(DEBUGMEM>1) err(warnmem,"lllgram1");
! 1096: tetpil=avma;
! 1097: sv=cgetg(4,t_VEC);
! 1098: sv[1]=lcopy(B); sv[2]=lcopy(u); sv[3]=lcopy(mu);
! 1099: sv=gerepile(av,tetpil,sv);
! 1100: B=(GEN)sv[1]; u=(GEN)sv[2]; mu=(GEN)sv[3];
! 1101: }
! 1102: }
! 1103: while (k<=n);
! 1104: tetpil=avma; return gerepile(av,tetpil,gcopy(u));
! 1105: }
! 1106:
! 1107: GEN
! 1108: lllgramint(GEN x)
! 1109: {
! 1110: return lllgramall0(x,lll_IM);
! 1111: }
! 1112:
! 1113: GEN
! 1114: lllgramkerim(GEN x)
! 1115: {
! 1116: return lllgramall0(x,lll_ALL);
! 1117: }
! 1118:
! 1119: /* This routine is functionally similar to lllint when all = 0.
! 1120: *
! 1121: * The input is an integer matrix mat (not necessarily square) whose
! 1122: * columns are linearly independent. It outputs another matrix T such that
! 1123: * mat * T is partially reduced. If all = 0, the size of mat * T is the
! 1124: * same as the size of mat. If all = 1 the number of columns of mat * T
! 1125: * is at most equal to its number of rows. A matrix M is said to be
! 1126: * -partially reduced- if | m1 +- m2 | >= |m1| for any two distinct
! 1127: * columns m1, m2, in M.
! 1128: *
! 1129: * This routine is designed to quickly reduce lattices in which one row
! 1130: * is huge compared to the other rows. For example, when searching for a
! 1131: * polynomial of degree 3 with root a mod p, the four input vectors might
! 1132: * be the coefficients of
! 1133: * X^3 - (a^3 mod p), X^2 - (a^2 mod p), X - (a mod p), p.
! 1134: * All four constant coefficients are O(p) and the rest are O(1). By the
! 1135: * pigeon-hole principle, the coefficients of the smallest vector in the
! 1136: * lattice are O(p^(1/4)), Hence significant reduction of vector lengths
! 1137: * can be anticipated.
! 1138: *
! 1139: * Peter Montgomery (July, 1994)
! 1140: *
! 1141: * If flag = 1 complete the reduction using lllint, otherwise return
! 1142: * partially reduced basis.
! 1143: */
! 1144: GEN
! 1145: lllintpartialall(GEN m, long flag)
! 1146: {
! 1147: const long ncol = lg(m)-1;
! 1148: const long ltop1 = avma;
! 1149: long ncolnz;
! 1150: GEN tm1, tm2, mid, *gptr[4];
! 1151:
! 1152: if (typ(m) != t_MAT) err(typeer,"lllintpartialall");
! 1153: if (ncol <= 1) return idmat(ncol);
! 1154:
! 1155: {
! 1156: GEN dot11 = sqscali((GEN)m[1]);
! 1157: GEN dot22 = sqscali((GEN)m[2]);
! 1158: GEN dot12 = gscali((GEN)m[1], (GEN)m[2]);
! 1159: GEN tm = idmat(2); /* For first two columns only */
! 1160:
! 1161: int progress = 0;
! 1162: long npass2 = 0;
! 1163:
! 1164: /* Try to row reduce the first two columns of m.
! 1165: * Our best result so far is (first two columns of m)*tm.
! 1166: *
! 1167: * Initially tm = 2 x 2 identity matrix.
! 1168: * The inner products of the reduced matrix are in
! 1169: * dot11, dot12, dot22.
! 1170: */
! 1171: while (npass2 < 2 || progress)
! 1172: {
! 1173: GEN dot12new,q = gdivround(dot12, dot22);
! 1174:
! 1175: npass2++; progress = signe(q);
! 1176: if (progress)
! 1177: {
! 1178: /* Conceptually replace (v1, v2) by (v1 - q*v2, v2),
! 1179: * where v1 and v2 represent the reduced basis for the
! 1180: * first two columns of the matrix.
! 1181: *
! 1182: * We do this by updating tm and the inner products.
! 1183: *
! 1184: * An improved algorithm would look only at the leading
! 1185: * digits of dot11, dot12, dot22. It would use
! 1186: * single-precision calculations as much as possible.
! 1187: */
! 1188: q = negi(q);
! 1189: dot12new = addii(dot12, mulii(q, dot22));
! 1190: dot11 = addii(dot11, mulii(q, addii(dot12, dot12new)));
! 1191: dot12 = dot12new;
! 1192: tm[1] = (long)lincomb_integral(gun,q, (GEN)tm[1],(GEN)tm[2]);
! 1193: }
! 1194:
! 1195: /* Interchange the output vectors v1 and v2. */
! 1196: gswap(dot11,dot22); swap(tm[1],tm[2]);
! 1197:
! 1198: /* Occasionally (including final pass) do garbage collection. */
! 1199: if (npass2 % 8 == 0 || !progress)
! 1200: {
! 1201: gptr[0] = &dot11; gptr[1] = &dot12;
! 1202: gptr[2] = &dot22; gptr[3] = &tm;
! 1203: gerepilemany(ltop1, gptr, 4);
! 1204: }
! 1205: } /* while npass2 < 2 || progress */
! 1206:
! 1207: {
! 1208: long icol;
! 1209: GEN det12 = subii(mulii(dot11, dot22), mulii(dot12, dot12));
! 1210:
! 1211: tm1 = idmat(ncol);
! 1212: mid = cgetg(ncol+1, t_MAT);
! 1213: for (icol = 1; icol <= 2; icol++)
! 1214: {
! 1215: coeff(tm1,1,icol) = coeff(tm,1,icol);
! 1216: coeff(tm1,2,icol) = coeff(tm,2,icol);
! 1217: mid[icol] = (long)lincomb_integral(
! 1218: gcoeff(tm,1,icol),gcoeff(tm,2,icol), (GEN)m[1],(GEN)m[2]);
! 1219: }
! 1220: for (icol = 3; icol <= ncol; icol++)
! 1221: {
! 1222: GEN curcol = (GEN)m[icol];
! 1223: GEN dot1i = gscali((GEN)mid[1], curcol);
! 1224: GEN dot2i = gscali((GEN)mid[2], curcol);
! 1225: /* Try to solve
! 1226: *
! 1227: * ( dot11 dot12 ) (q1) ( dot1i )
! 1228: * ( dot12 dot22 ) (q2) = ( dot2i )
! 1229: *
! 1230: * Round -q1 and -q2 to the nearest integer.
! 1231: * Then compute curcol - q1*mid[1] - q2*mid[2].
! 1232: * This will be approximately orthogonal to the
! 1233: * first two vectors in the new basis.
! 1234: */
! 1235: GEN q1neg = subii(mulii(dot12, dot2i), mulii(dot22, dot1i));
! 1236: GEN q2neg = subii(mulii(dot12, dot1i), mulii(dot11, dot2i));
! 1237:
! 1238: q1neg = gdivround(q1neg, det12);
! 1239: q2neg = gdivround(q2neg, det12);
! 1240: coeff(tm1, 1, icol) = ladd(gmul(q1neg, gcoeff(tm,1,1)),
! 1241: gmul(q2neg, gcoeff(tm,1,2)));
! 1242: coeff(tm1, 2, icol) = ladd(gmul(q1neg, gcoeff(tm,2,1)),
! 1243: gmul(q2neg, gcoeff(tm,2,2)));
! 1244: mid[icol] = ladd(curcol,
! 1245: lincomb_integral(q1neg,q2neg, (GEN)mid[1],(GEN)mid[2]));
! 1246: } /* for icol */
! 1247: } /* local block */
! 1248: }
! 1249: if (DEBUGLEVEL>4)
! 1250: {
! 1251: fprintferr("tm1 = "); outbeauterr(tm1);
! 1252: fprintferr("mid = "); outbeauterr(mid);
! 1253: }
! 1254: gptr[0] = &tm1; gptr[1] = ∣
! 1255: gerepilemany(ltop1, gptr, 2);
! 1256: {
! 1257: /* For each pair of column vectors v and w in mid * tm2,
! 1258: * try to replace (v, w) by (v, v - q*w) for some q.
! 1259: * We compute all inner products and check them repeatedly.
! 1260: */
! 1261: const long ltop3 = avma; /* Excludes region with tm1 and mid */
! 1262: long icol, lim, reductions, npass = 0;
! 1263: GEN dotprd = cgetg(ncol+1, t_MAT);
! 1264:
! 1265: tm2 = idmat(ncol);
! 1266: for (icol=1; icol <= ncol; icol++)
! 1267: {
! 1268: long jcol;
! 1269:
! 1270: dotprd[icol] = lgetg(ncol+1,t_COL);
! 1271: for (jcol=1; jcol <= icol; jcol++)
! 1272: coeff(dotprd,jcol,icol) = coeff(dotprd,icol,jcol) =
! 1273: (long)gscali((GEN)mid[icol], (GEN)mid[jcol]);
! 1274: } /* for icol */
! 1275: lim = stack_lim(ltop3,1);
! 1276: for(;;)
! 1277: {
! 1278: reductions = 0;
! 1279: for (icol=1; icol <= ncol; icol++)
! 1280: {
! 1281: long ijdif, jcol, k1, k2;
! 1282: GEN codi, q;
! 1283:
! 1284: for (ijdif=1; ijdif < ncol; ijdif++)
! 1285: {
! 1286: const long previous_avma = avma;
! 1287:
! 1288: jcol = (icol + ijdif - 1) % ncol; jcol++; /* Hack for NeXTgcc 2.5.8 */
! 1289: k1 = (cmpii(gcoeff(dotprd,icol,icol),
! 1290: gcoeff(dotprd,jcol,jcol) ) > 0)
! 1291: ? icol : jcol; /* index of larger column */
! 1292: k2 = icol + jcol - k1; /* index of smaller column */
! 1293: codi = gcoeff(dotprd,k2,k2);
! 1294: q = gcmp0(codi)? gzero: gdivround(gcoeff(dotprd,k1,k2), codi);
! 1295:
! 1296: /* Try to subtract a multiple of column k2 from column k1. */
! 1297: if (gcmp0(q)) avma = previous_avma;
! 1298: else
! 1299: {
! 1300: long dcol;
! 1301:
! 1302: reductions++; q = negi(q);
! 1303: tm2[k1]=(long)
! 1304: lincomb_integral(gun,q, (GEN)tm2[k1], (GEN)tm2[k2]);
! 1305: dotprd[k1]=(long)
! 1306: lincomb_integral(gun,q, (GEN)dotprd[k1], (GEN)dotprd[k2]);
! 1307: coeff(dotprd, k1, k1) = laddii(gcoeff(dotprd,k1,k1),
! 1308: mulii(q, gcoeff(dotprd,k2,k1)));
! 1309: for (dcol = 1; dcol <= ncol; dcol++)
! 1310: coeff(dotprd,k1,dcol) = coeff(dotprd,dcol,k1);
! 1311: } /* if q != 0 */
! 1312: } /* for ijdif */
! 1313: if (low_stack(lim, stack_lim(ltop3,1)))
! 1314: {
! 1315: if(DEBUGMEM>1) err(warnmem,"lllintpartialall");
! 1316: gptr[0] = &dotprd; gptr[1] = &tm2;
! 1317: gerepilemany(ltop3, gptr, 2);
! 1318: }
! 1319: } /* for icol */
! 1320: if (!reductions) break;
! 1321: if (DEBUGLEVEL>4)
! 1322: {
! 1323: GEN diag_prod = dbltor(1.0);
! 1324: for (icol = 1; icol <= ncol; icol++)
! 1325: diag_prod = gmul(diag_prod, gcoeff(dotprd,icol,icol));
! 1326: npass++;
! 1327: fprintferr("npass = %ld, red. last time = %ld, diag_prod = %Z\n\n",
! 1328: npass, reductions, diag_prod);
! 1329: }
! 1330: } /* for(;;)*/
! 1331:
! 1332: /* Sort columns so smallest comes first in m * tm1 * tm2.
! 1333: * Use insertion sort.
! 1334: */
! 1335: for (icol = 1; icol < ncol; icol++)
! 1336: {
! 1337: long jcol, s = icol;
! 1338:
! 1339: for (jcol = icol+1; jcol <= ncol; jcol++)
! 1340: if (cmpii(gcoeff(dotprd,s,s),gcoeff(dotprd,jcol,jcol)) > 0) s = jcol;
! 1341: if (icol != s)
! 1342: { /* Exchange with proper column */
! 1343: /* Only diagonal of dotprd is updated */
! 1344: swap(tm2[icol], tm2[s]);
! 1345: swap(coeff(dotprd,icol,icol), coeff(dotprd,s,s));
! 1346: }
! 1347: } /* for icol */
! 1348: icol=1;
! 1349: while (icol <= ncol && !signe(gcoeff(dotprd,icol,icol))) icol++;
! 1350: ncolnz = ncol - icol + 1;
! 1351: } /* local block */
! 1352:
! 1353: if (flag)
! 1354: {
! 1355: if (ncolnz == lg((GEN)m[1])-1)
! 1356: {
! 1357: tm2 += (ncol-ncolnz);
! 1358: tm2[0]=evaltyp(t_MAT)|evallg(ncolnz+1);
! 1359: }
! 1360: else
! 1361: {
! 1362: tm1 = gmul(tm1, tm2);
! 1363: tm2 = lllint(gmul(m, tm1));
! 1364: }
! 1365: }
! 1366: if (DEBUGLEVEL>4)
! 1367: { fprintferr("lllintpartial output = "); outbeauterr(gmul(tm1, tm2)); }
! 1368: return gerepileupto(ltop1, gmul(tm1, tm2));
! 1369: }
! 1370:
! 1371: GEN
! 1372: lllintpartial(GEN mat)
! 1373: {
! 1374: return lllintpartialall(mat,1);
! 1375: }
! 1376:
! 1377: /********************************************************************/
! 1378: /** **/
! 1379: /** LINEAR & ALGEBRAIC DEPENDANCE **/
! 1380: /** **/
! 1381: /********************************************************************/
! 1382:
! 1383: GEN
! 1384: lindep0(GEN x,long bit,long prec)
! 1385: {
! 1386: if (!bit) return lindep(x,prec);
! 1387: if (bit>0) return lindep2(x,bit);
! 1388: return deplin(x);
! 1389: }
! 1390:
! 1391: GEN
! 1392: lindep2(GEN x, long bit)
! 1393: {
! 1394: long tx=typ(x),lx=lg(x),ly,i,j,flag,e,tetpil, av = avma;
! 1395: GEN re,im,p1,p2;
! 1396:
! 1397: if (! is_vec_t(tx)) err(typeer,"lindep2");
! 1398: if (lx<=2) return cgetg(1,t_VEC);
! 1399: re=greal(x); im=gimag(x); flag = !gcmp0(im);
! 1400: if (lx == 3)
! 1401: { /* independant over R ? */
! 1402: if (gexpo(gsub(gmul((GEN)re[1],(GEN)im[2]),
! 1403: gmul((GEN)re[2],(GEN)im[1]))) > - bit)
! 1404: {
! 1405: avma = av; return cgetg(1, t_VEC);
! 1406: }
! 1407: }
! 1408: ly = flag? lx+2: lx+1;
! 1409: p2=cgetg(lx,t_MAT); bit = (long) (bit/L2SL10);
! 1410: for (i=1; i<lx; i++)
! 1411: {
! 1412: p1=cgetg(ly,t_COL); p2[i]=(long)p1;
! 1413: for (j=1; j<lx; j++) p1[j]=(i==j) ? un : zero;
! 1414: p1[lx]=lcvtoi(gshift((GEN)re[i],bit),&e);
! 1415: if (flag) p1[lx+1]=lcvtoi(gshift((GEN)im[i],bit),&e);
! 1416: }
! 1417: p1=gmul(p2,lllint(p2)); p1=(GEN)p1[1];
! 1418: p1[0]=evaltyp(t_VEC) | evallg(lx); tetpil=avma;
! 1419: return gerepile(av,tetpil,gcopy(p1));
! 1420: }
! 1421:
! 1422: #define quazero(x) (gcmp0(x) || (typ(x)==t_REAL && expo(x) < EXP))
! 1423: GEN
! 1424: lindep(GEN x, long prec)
! 1425: {
! 1426: GEN *b,*be,*bn,**m,qzer;
! 1427: GEN c1,c2,c3,px,py,pxy,re,im,p3,p4,r,f,em;
! 1428: long i,j,fl,i1, lx = lg(x), tx = typ(x), n = lx-1;
! 1429: long av = avma, lim = stack_lim(av,1), av0,av1,tetpil;
! 1430: const long EXP = - bit_accuracy(prec) + 2*n;
! 1431:
! 1432: if (! is_vec_t(tx)) err(typeer,"lindep");
! 1433: if (lx<=2) return cgetg(1,t_VEC);
! 1434: x = gmul(x, realun(prec));
! 1435: re=greal(x); im=gimag(x);
! 1436: if (lx == 3)
! 1437: { /* independant over R ? */
! 1438: if (gexpo(gsub(gmul((GEN)re[1],(GEN)im[2]),
! 1439: gmul((GEN)re[2],(GEN)im[1]))) > - bit_accuracy(prec))
! 1440: {
! 1441: avma = av; return cgetg(1, t_VEC);
! 1442: }
! 1443: }
! 1444: qzer = new_chunk(lx);
! 1445: b = (GEN*)idmat(n);
! 1446: be= (GEN*)new_chunk(lx);
! 1447: bn= (GEN*)new_chunk(lx);
! 1448: m = (GEN**)new_chunk(lx);
! 1449: for (i=1; i<=n; i++)
! 1450: {
! 1451: bn[i]=cgetr(prec+1);
! 1452: be[i]=cgetg(lx,t_COL);
! 1453: m[i] = (GEN*)new_chunk(lx);
! 1454: for (j=1; j<i ; j++) m[i][j]=cgetr(prec+1);
! 1455: for (j=1; j<=n; j++) be[i][j]=lgetr(prec+1);
! 1456: }
! 1457: px=sqscal(re);
! 1458: py=sqscal(im); pxy=gscal(re,im);
! 1459: p3=mpsub(mpmul(px,py),gsqr(pxy));
! 1460: if (quazero(re)) { re=im; px=py; fl=1; } else fl=quazero(p3);
! 1461: av0 = av1 = avma;
! 1462: for (i=1; i<=n; i++)
! 1463: {
! 1464: p4 = gscal(b[i],re);
! 1465: if (fl) p4=gmul(gdiv(p4,px),re);
! 1466: else
! 1467: {
! 1468: GEN p5,p6,p7;
! 1469: p5=gscal(b[i],im);
! 1470: p6=gdiv(gsub(gmul(py,p4),gmul(pxy,p5)),p3);
! 1471: p7=gdiv(gsub(gmul(px,p5),gmul(pxy,p4)),p3);
! 1472: p4=gadd(gmul(p6,re),gmul(p7,im));
! 1473: }
! 1474: if (tx!=t_COL) p4=gtrans(p4);
! 1475: p4=gsub(b[i],p4);
! 1476: for (j=1; j<i; j++)
! 1477: if (qzer[j]) affrr(bn[j],m[i][j]);
! 1478: else
! 1479: {
! 1480: gdivz(gscal(b[i],be[j]),bn[j],m[i][j]);
! 1481: p4=gsub(p4,gmul(m[i][j],be[j]));
! 1482: }
! 1483: gaffect(p4,be[i]); affrr(sqscal(be[i]),bn[i]);
! 1484: qzer[i]=quazero(bn[i]); avma=av1;
! 1485: }
! 1486: while (qzer[n])
! 1487: {
! 1488: long e;
! 1489: if (DEBUGLEVEL>9)
! 1490: {
! 1491: fprintferr("qzer[%ld]=%ld\n",n,qzer[n]);
! 1492: for (i1=1; i1<=n; i1++)
! 1493: for (i=1; i<i1; i++) output(m[i1][i]);
! 1494: }
! 1495: em=bn[1]; j=1;
! 1496: for (i=2; i<n; i++)
! 1497: {
! 1498: p3=shiftr(bn[i],i);
! 1499: if (cmprr(p3,em)>0) { em=p3; j=i; }
! 1500: }
! 1501: i=j; i1=i+1;
! 1502: avma = av1; r = grndtoi(m[i1][i], &e);
! 1503: if (e >= 0) err(talker,"precision too low in lindep");
! 1504: r = negi(r);
! 1505: p3 = lincomb_integral(gun,r, b[i1],b[i]);
! 1506: av1 = avma;
! 1507: b[i1]=b[i]; b[i]=p3; f=addir(r,m[i1][i]);
! 1508: for (j=1; j<i; j++)
! 1509: if (!qzer[j])
! 1510: {
! 1511: p3=mpadd(m[i1][j],mulir(r,m[i][j]));
! 1512: affrr(m[i][j],m[i1][j]); mpaff(p3,m[i][j]);
! 1513: }
! 1514: c1=addrr(bn[i1],mulrr(gsqr(f),bn[i]));
! 1515: if (!quazero(c1))
! 1516: {
! 1517: c2=divrr(mulrr(bn[i],f),c1); affrr(c2,m[i1][i]);
! 1518: c3=divrr(bn[i1],c1); mulrrz(c3,bn[i],bn[i1]);
! 1519: affrr(c1,bn[i]); qzer[i1]=quazero(bn[i1]); qzer[i]=0;
! 1520: for (j=i+2; j<=n; j++)
! 1521: {
! 1522: p3=addrr(mulrr(m[j][i1],c3),mulrr(m[j][i],c2));
! 1523: subrrz(m[j][i],mulrr(f,m[j][i1]), m[j][i1]);
! 1524: affrr(p3,m[j][i]);
! 1525: }
! 1526: }
! 1527: else
! 1528: {
! 1529: qzer[i1]=qzer[i]; qzer[i]=1;
! 1530: affrr(bn[i],bn[i1]); affrr(c1,bn[i]);
! 1531: for (j=i+2; j<=n; j++) affrr(m[j][i],m[j][i1]);
! 1532: }
! 1533: if (low_stack(lim, stack_lim(av,1)))
! 1534: {
! 1535: if(DEBUGMEM>1) err(warnmem,"lindep");
! 1536: b = (GEN*)gerepileupto(av0, gcopy((GEN)b));
! 1537: av1 = avma;
! 1538: }
! 1539: }
! 1540: p3=cgetg(lx,t_COL); p3[n]=un; for (i=1; i<n; i++) p3[i]=zero;
! 1541: p4 = (GEN)b; p4[0] = evaltyp(t_MAT) | evallg(lx);
! 1542: p4=gauss(gtrans(p4),p3); tetpil=avma;
! 1543: return gerepile(av,tetpil,gtrans(p4));
! 1544: }
! 1545:
! 1546: GEN
! 1547: algdep0(GEN x, long n, long bit, long prec)
! 1548: {
! 1549: long tx=typ(x),av,i,k;
! 1550: GEN y,p1;
! 1551:
! 1552: if (! is_scalar_t(tx)) err(typeer,"algdep0");
! 1553: if (tx==t_POLMOD) { y=forcecopy((GEN)x[1]); setvarn(y,0); return y; }
! 1554: if (gcmp0(x)) return gzero;
! 1555: if (!n) return gun;
! 1556:
! 1557: av=avma; p1=cgetg(n+2,t_COL); p1[1]=un;
! 1558: for (i=2; i<=n+1; i++) p1[i]=lmul((GEN)p1[i-1],x);
! 1559:
! 1560: p1 = bit? lindep2(p1,bit): lindep(p1,prec);
! 1561: if (lg(p1) < 2)
! 1562: err(talker,"higher degree than expected in algdep");
! 1563:
! 1564: y=cgetg(n+3,t_POL);
! 1565: y[1] = evalsigne(1) | evalvarn(0);
! 1566: k=1; while (gcmp0((GEN)p1[k])) k++;
! 1567: for (i=0; i<=n+1-k; i++) y[i+2] = p1[k+i];
! 1568: normalizepol_i(y, n+4-k);
! 1569: y = (gsigne(leading_term(y)) > 0)? gcopy(y): gneg(y);
! 1570: return gerepileupto(av,y);
! 1571: }
! 1572:
! 1573: GEN
! 1574: algdep2(GEN x, long n, long bit)
! 1575: {
! 1576: return algdep0(x,n,bit,0);
! 1577: }
! 1578:
! 1579: GEN
! 1580: algdep(GEN x, long n, long prec)
! 1581: {
! 1582: return algdep0(x,n,0,prec);
! 1583: }
! 1584:
! 1585: /********************************************************************/
! 1586: /** **/
! 1587: /** INTEGRAL KERNEL (LLL REDUCED) **/
! 1588: /** **/
! 1589: /********************************************************************/
! 1590:
! 1591: GEN
! 1592: matkerint0(GEN x, long flag)
! 1593: {
! 1594: switch(flag)
! 1595: {
! 1596: case 0: return kerint(x);
! 1597: case 1: return kerint1(x);
! 1598: case 2: return kerint2(x);
! 1599: default: err(flagerr,"matkerint");
! 1600: }
! 1601: return NULL; /* not reached */
! 1602: }
! 1603:
! 1604: GEN
! 1605: kerint1(GEN x)
! 1606: {
! 1607: long av=avma,tetpil;
! 1608: GEN p1,p2;
! 1609:
! 1610: p1=matrixqz3(ker(x)); p2=lllint(p1); tetpil=avma;
! 1611: return gerepile(av,tetpil,gmul(p1,p2));
! 1612: }
! 1613:
! 1614: GEN
! 1615: kerint2(GEN x)
! 1616: {
! 1617: long lx=lg(x), i,j,av,av1;
! 1618: GEN g,p1;
! 1619:
! 1620: if (typ(x) != t_MAT) err(typeer,"kerint2");
! 1621: av=avma; g=cgetg(lx,t_MAT);
! 1622: for (j=1; j<lx; j++) g[j]=lgetg(lx,t_COL);
! 1623: for (i=1; i<lx; i++)
! 1624: for (j=1; j<=i; j++)
! 1625: coeff(g,i,j) = coeff(g,j,i) = (long)gscal((GEN)x[i],(GEN)x[j]);
! 1626: g=lllgramall0(g,lll_KER); p1=lllint(g);
! 1627: av1=avma; return gerepile(av,av1,gmul(g,p1));
! 1628: }
! 1629:
! 1630: static GEN
! 1631: lllall0(GEN x, long flag)
! 1632: {
! 1633: long av0=avma,av,tetpil,lx=lg(x),i,j,k,l,n,lim,kmax;
! 1634: GEN u,B,L,q,r,h,la,p1,p2,p4,fl;
! 1635:
! 1636: if (typ(x) != t_MAT) err(typeer,"lllall0");
! 1637: n=lx-1; if (n<=1) return lllall_trivial(x,n, flag | lll_GRAM);
! 1638: fl = new_chunk(lx);
! 1639:
! 1640: av=avma; lim=stack_lim(av,1); x=dummycopy(x);
! 1641: B=gscalcol(gun, lx);
! 1642: L=cgetg(lx,t_MAT);
! 1643: for (k=lg(x[1]),j=1; j<lx; j++)
! 1644: {
! 1645: for (i=1; i<k; i++)
! 1646: if (typ(gcoeff(x,i,j))!=t_INT) err(lllger4);
! 1647: fl[j] = 0; L[j] = (long)zerocol(n);
! 1648: }
! 1649: k=2; h=idmat(n); kmax=1;
! 1650: u=sqscali((GEN)x[1]);
! 1651: if (signe(u)) { B[2]=(long)u; coeff(L,1,1)=un; fl[1]=1; } else B[2]=un;
! 1652: for(;;)
! 1653: {
! 1654: if (k > kmax)
! 1655: {
! 1656: kmax = k;
! 1657: for (j=1; j<=k; j++)
! 1658: {
! 1659: if (j==k || fl[j])
! 1660: {
! 1661: long av1 = avma;
! 1662: u=gscali((GEN)x[k],(GEN)x[j]);
! 1663: for (i=1; i<j; i++)
! 1664: if (fl[i])
! 1665: u = divii(subii(mulii((GEN)B[i+1],u),
! 1666: mulii(gcoeff(L,k,i),gcoeff(L,j,i))),
! 1667: (GEN)B[i]);
! 1668: u = gerepileuptoint(av1, u);
! 1669: if (j<k) coeff(L,k,j)=(long)u;
! 1670: else
! 1671: {
! 1672: if (signe(u)) { B[k+1]=(long)u; coeff(L,k,k)=un; fl[k]=1; }
! 1673: else { B[k+1]=B[k]; fl[k]=0; }
! 1674: }
! 1675: }
! 1676: }
! 1677: }
! 1678: if (fl[k-1] && !fl[k])
! 1679: {
! 1680: u = shifti(gcoeff(L,k,k-1),1);
! 1681: if (absi_cmp(u, (GEN)B[k]) > 0)
! 1682: {
! 1683: q = truedvmdii(addii(u,(GEN)B[k]),shifti((GEN)B[k],1), NULL);
! 1684: r = negi(q);
! 1685: h[k] = (long)lincomb_integral(gun,r, (GEN)h[k],(GEN)h[k-1]);
! 1686: x[k] = (long)lincomb_integral(gun,r, (GEN)x[k],(GEN)x[k-1]);
! 1687: coeff(L,k,k-1)=laddii(gcoeff(L,k,k-1),mulii(r,(GEN)B[k]));
! 1688: for (i=1; i<k-1; i++)
! 1689: coeff(L,k,i)=laddii(gcoeff(L,k,i),mulii(r,gcoeff(L,k-1,i)));
! 1690: }
! 1691: la=gcoeff(L,k,k-1); p4=sqri(la);
! 1692: swap(h[k-1], h[k]);
! 1693: swap(x[k-1], x[k]);
! 1694: for (j=1; j<k-1; j++) swap(coeff(L,k-1,j), coeff(L,k,j));
! 1695: if (signe(la))
! 1696: {
! 1697: p2=(GEN)B[k]; p1=divii(p4,p2);
! 1698: B[k+1]=B[k]=(long)p1;
! 1699: for (i=k+1; i<=kmax; i++)
! 1700: coeff(L,i,k-1)=ldivii(mulii(la,gcoeff(L,i,k-1)),p2);
! 1701: for (j=k+1; j<kmax; j++)
! 1702: for (i=j+1; i<=kmax; i++)
! 1703: coeff(L,i,j)=ldivii(mulii((GEN)p1,gcoeff(L,i,j)),p2);
! 1704: for (i=k+2; i<=kmax+1; i++)
! 1705: B[i]=ldivii(mulii((GEN)p1,(GEN)B[i]),p2);
! 1706: }
! 1707: else
! 1708: {
! 1709: for (i=k+1; i<=kmax; i++)
! 1710: { coeff(L,i,k)=coeff(L,i,k-1); coeff(L,i,k-1)=zero; }
! 1711: B[k]=B[k-1]; fl[k]=1; fl[k-1]=0;
! 1712: }
! 1713: if (k>2) k--;
! 1714: }
! 1715: else
! 1716: {
! 1717: for (l=k-1; l>=1; l--)
! 1718: {
! 1719: u = shifti(gcoeff(L,k,l),1);
! 1720: if (absi_cmp(u,(GEN)B[l+1]) > 0)
! 1721: {
! 1722: q = truedvmdii(addii(u,(GEN)B[l+1]),shifti((GEN)B[l+1],1), NULL);
! 1723: r = negi(q);
! 1724: h[k] = (long)lincomb_integral(gun,r,(GEN)h[k],(GEN)h[l]);
! 1725: x[k] = (long)lincomb_integral(gun,r,(GEN)x[k],(GEN)x[l]);
! 1726: coeff(L,k,l)=laddii(gcoeff(L,k,l),mulii(r,(GEN)B[l+1]));
! 1727: for (i=1; i<l; i++)
! 1728: coeff(L,k,i)=laddii(gcoeff(L,k,i),mulii(r,gcoeff(L,l,i)));
! 1729: }
! 1730: }
! 1731: if (++k > n) break;
! 1732: }
! 1733: if (low_stack(lim, stack_lim(av,1)))
! 1734: {
! 1735: GEN *gptr[4];
! 1736: if(DEBUGMEM>1) err(warnmem,"lllall0");
! 1737: gptr[0]=&B; gptr[1]=&L; gptr[2]=&h;
! 1738: gptr[3]=&x; gerepilemany(av,gptr,4);
! 1739: }
! 1740: }
! 1741: tetpil=avma;
! 1742: return gerepile(av0,tetpil,lllgramall_finish(h,fl,flag,n));
! 1743: }
! 1744:
! 1745: GEN
! 1746: kerint(GEN x)
! 1747: {
! 1748: long av=avma,av1;
! 1749: GEN g,p1;
! 1750:
! 1751: g=lllall0(x,lll_KER); if (lg(g)==1) return g;
! 1752: p1=lllint(g); av1=avma;
! 1753: return gerepile(av,av1,gmul(g,p1));
! 1754: }
! 1755:
! 1756: /********************************************************************/
! 1757: /** **/
! 1758: /** POLRED & CO. **/
! 1759: /** **/
! 1760: /********************************************************************/
! 1761: /* remove duplicate polynomials in y, updating a (same indexes), in place */
! 1762: static long
! 1763: remove_duplicates(GEN y, GEN a)
! 1764: {
! 1765: long k,i, nv = lg(y), av = avma;
! 1766: GEN z;
! 1767:
! 1768: if (nv < 2) return nv;
! 1769: z = new_chunk(3);
! 1770: z[1] = (long)y;
! 1771: z[2] = (long)a; (void)sort_factor(z, gcmp);
! 1772: for (k=1, i=2; i<nv; i++)
! 1773: if (!gegal((GEN)y[i], (GEN)y[k]))
! 1774: {
! 1775: k++;
! 1776: a[k] = a[i];
! 1777: y[k] = y[i];
! 1778: }
! 1779: nv = k+1; setlg(a,nv); setlg(y,nv);
! 1780: avma = av; return nv;
! 1781: }
! 1782:
! 1783: /* in place; make sure second non-zero coeff is negative (choose x or -x) */
! 1784: int
! 1785: canon_pol(GEN z)
! 1786: {
! 1787: long i,s;
! 1788:
! 1789: for (i=lgef(z)-2; i>=2; i-=2)
! 1790: {
! 1791: s = signe(z[i]);
! 1792: if (s > 0)
! 1793: {
! 1794: for (; i>=2; i-=2) z[i]=lnegi((GEN)z[i]);
! 1795: return -1;
! 1796: }
! 1797: if (s) return 1;
! 1798: }
! 1799: return 0;
! 1800: }
! 1801:
! 1802: static GEN
! 1803: pols_for_polred(GEN x, GEN base, GEN LLLbase, GEN *pta)
! 1804: {
! 1805: long i,j, v = varn(x), n = lg(base);
! 1806: GEN p1,p2,p3,y, a = cgetg(n,t_VEC);
! 1807:
! 1808: for (i=1; i<n; i++) a[i] = lmul(base,(GEN)LLLbase[i]);
! 1809: y=cgetg(n,t_VEC);
! 1810: for (i=1; i<n; i++)
! 1811: {
! 1812: if (DEBUGLEVEL>2) { fprintferr("i = %ld\n",i); flusherr(); }
! 1813: p1=(GEN)a[i]; p3=content(p1);
! 1814: if (gcmp1(p3)) p3 = NULL; else p1 = gdiv(p1,p3);
! 1815: p1 = caract2(x,p1,v);
! 1816: if (p3)
! 1817: for (p2=gun, j=lgef(p1)-2; j>=2; j--)
! 1818: {
! 1819: p2 = gmul(p2,p3); p1[j] = lmul((GEN)p1[j], p2);
! 1820: }
! 1821: p2=modulargcd(derivpol(p1),p1); p3=leading_term(p2);
! 1822: if (!gcmp1(p3)) p2=gdiv(p2,p3);
! 1823: p1 = gdiv(p1,p2);
! 1824: if (canon_pol(p1) < 0 && pta)
! 1825: a[i] = (long) gneg_i((GEN)a[i]);
! 1826: y[i] = (long)p1; if (DEBUGLEVEL>=4) outerr(p1);
! 1827: }
! 1828: (void)remove_duplicates(y,a);
! 1829: if (pta) *pta = a;
! 1830: return y;
! 1831: }
! 1832:
! 1833: GEN
! 1834: nf_get_T2(GEN base, GEN polr)
! 1835: {
! 1836: long i,j, n = lg(base);
! 1837: GEN p1,p2=cgetg(n,t_MAT);
! 1838:
! 1839: for (i=1; i<n; i++)
! 1840: {
! 1841: p1=cgetg(n,t_COL); p2[i]=(long)p1;
! 1842: for (j=1; j<n; j++)
! 1843: p1[j] = (long) poleval((GEN)base[i],(GEN)polr[j]);
! 1844: }
! 1845: return mulmat_real(gconj(gtrans(p2)),p2);
! 1846: }
! 1847:
! 1848: /* Return the base change matrix giving the an LLL-reduced basis for the
! 1849: * maximal order of the nf given by x. Expressed in terms of the standard
! 1850: * HNF basis (as polynomials) given in base (ignored if x is an nf)
! 1851: */
! 1852: GEN
! 1853: LLL_nfbasis(GEN *ptx, GEN polr, GEN base, long prec)
! 1854: {
! 1855: GEN T2,p1, x = *ptx;
! 1856: int totally_real,n,i,j;
! 1857:
! 1858: if (typ(x) != t_POL)
! 1859: {
! 1860: p1=checknf(x); *ptx=x=(GEN)p1[1];
! 1861: base=(GEN)p1[7]; n=lgef(x)-3;
! 1862: totally_real = !signe(gmael(p1,2,2));
! 1863: T2=gmael(p1,5,3); if (totally_real) T2 = ground(T2);
! 1864: }
! 1865: else
! 1866: {
! 1867: n=lgef(x)-3; totally_real = (!prec || sturm(x)==n);
! 1868: if (!totally_real)
! 1869: T2 = nf_get_T2(base,polr? polr: roots(x,prec));
! 1870: else
! 1871: { /* totally real */
! 1872: GEN ptrace=cgetg(n+2,t_VEC);
! 1873: long k;
! 1874:
! 1875: ptrace[2]=lstoi(n);
! 1876: for (k=2; k<=n; k++)
! 1877: { /* cf polsym */
! 1878: GEN y = x + (n-k+1);
! 1879: p1 = gmulsg(k-1,(GEN)y[2]);
! 1880: for (i=3; i<=k; i++)
! 1881: p1 = gadd(p1,gmul((GEN)y[i],(GEN)ptrace[i]));
! 1882: ptrace[i] = lneg(p1);
! 1883: }
! 1884: T2=cgetg(n+1,t_MAT);
! 1885: for (i=1; i<=n; i++)
! 1886: {
! 1887: p1=cgetg(n+1,t_COL); T2[i]=(long)p1;
! 1888: for (j=1; j<i ; j++) p1[j] = coeff(T2,i,j);
! 1889: for ( ; j<=n; j++)
! 1890: { /* cf quicktrace */
! 1891: GEN p2 = gres(gmul((GEN)base[i],(GEN)base[j]),x);
! 1892: GEN p3 = gzero;
! 1893: for (k=lgef(p2)-1; k>1; k--)
! 1894: p3 = gadd(p3, gmul((GEN)p2[k],(GEN)ptrace[k]));
! 1895: p1[j]=(long)p3;
! 1896: }
! 1897: }
! 1898: }
! 1899: }
! 1900: if (totally_real) return lllgramint(T2);
! 1901: for (i=1; ; i++)
! 1902: {
! 1903: if ((p1 = lllgramintern(T2,100,1,prec))) return p1;
! 1904: if (i == MAXITERPOL) err(accurer,"allpolred");
! 1905: prec=(prec<<1)-2;
! 1906: if (DEBUGLEVEL) err(warnprec,"allpolred",prec);
! 1907: T2=nf_get_T2(base,roots(x,prec));
! 1908: }
! 1909: }
! 1910:
! 1911: /* x can be a polynomial, but also an nf or a bnf */
! 1912: static GEN
! 1913: allpolred(GEN x, GEN *pta, long code, long prec)
! 1914: {
! 1915: GEN y,p1, base = NULL, polr = NULL;
! 1916: long av = avma;
! 1917:
! 1918: if (typ(x) == t_POL)
! 1919: {
! 1920: if (!signe(x)) return gcopy(x);
! 1921: check_pol_int(x);
! 1922: if (!gcmp1(leading_term(x)))
! 1923: err(impl,"allpolred for nonmonic polynomials");
! 1924: base = allbase4(x,code,&p1,NULL); /* p1 is junk */
! 1925: }
! 1926: else
! 1927: {
! 1928: long i = lg(x);
! 1929: if (typ(x) == t_VEC && i<=4 && i>=3 && typ(x[1])==t_POL)
! 1930: { /* polynomial + integer basis */
! 1931: base=(GEN)x[2]; x=(GEN)x[1];
! 1932: }
! 1933: }
! 1934: p1 = LLL_nfbasis(&x,polr,base,prec);
! 1935: y = pols_for_polred(x,base,p1,pta);
! 1936: if (pta)
! 1937: {
! 1938: GEN *gptr[2]; gptr[0]=&y; gptr[1]=pta;
! 1939: gerepilemany(av,gptr,2); return y;
! 1940: }
! 1941: return gerepileupto(av,y);
! 1942: }
! 1943:
! 1944: GEN
! 1945: polred0(GEN x,long flag, GEN p, long prec)
! 1946: {
! 1947: GEN y;
! 1948: long smll = (flag & 1);
! 1949:
! 1950: if (p && !gcmp0(p)) smll=(long) p; /* factored polred */
! 1951: if (flag & 2) /* polred2 */
! 1952: {
! 1953: y=cgetg(3,t_MAT);
! 1954: y[2]=(long)allpolred(x,(GEN*)(y+1),smll,prec);
! 1955: return y;
! 1956: }
! 1957: return allpolred(x,NULL,smll,prec);
! 1958: }
! 1959:
! 1960: GEN
! 1961: ordred(GEN x, long prec)
! 1962: {
! 1963: GEN p1,p2,p3,base;
! 1964: long n=lgef(x)-3,i,j,av=avma,v = varn(x),tetpil;
! 1965:
! 1966: if (typ(x) != t_POL) err(typeer,"ordred");
! 1967: if (!signe(x)) return gcopy(x);
! 1968: if (!gcmp1((GEN)x[n+2])) err(impl,"ordred for nonmonic polynomials");
! 1969:
! 1970: p2=roots(x,prec); p3=cgetg(n+1,t_MAT);
! 1971:
! 1972: base=cgetg(n+1,t_VEC); /* power basis */
! 1973: for (i=1; i<=n; i++)
! 1974: {
! 1975: base[i]=lpowgs(polx[v],i-1);
! 1976: p1=cgetg(n+1,t_COL); p3[i]=(long)p1;
! 1977: for (j=1; j<=n; j++)
! 1978: p1[j]=lpuigs((GEN)p2[j],i-1);
! 1979: }
! 1980: p2 = mulmat_real(gconj(gtrans(p3)),p3);
! 1981: p1 = pols_for_polred(x,base,lllgram(p2,prec), NULL);
! 1982:
! 1983: tetpil=avma; return gerepile(av,tetpil, gcopy(p1));
! 1984: }
! 1985:
! 1986: GEN roots_to_pol_r1r2(GEN a, long r1, long v);
! 1987: static GEN chk_basis_embed;
! 1988: long chk_r1, chk_v;
! 1989:
! 1990: static GEN
! 1991: init_chk(GEN nf, GEN mat, GEN bound)
! 1992: {
! 1993: GEN M = gmael(nf,5,1);
! 1994: long n = lg(nf[7])-1, prec,prec2;
! 1995:
! 1996: if (!bound)
! 1997: {
! 1998: chk_v = varn(nf[1]);
! 1999: chk_r1 = itos(gmael(nf,2,1));
! 2000: chk_basis_embed = gmul(M, mat);
! 2001: return gmul((GEN)nf[7],mat);
! 2002: }
! 2003:
! 2004: /* should be [max_k C^n_k (B/k) ^ k] ^ (1/2) */
! 2005: bound = gsqrt(gpowgs(bound, n), 3);
! 2006: prec2 = (1 + (gexpo(bound) >> TWOPOTBITS_IN_LONG));
! 2007: if (prec2 < 0) prec2 = 0;
! 2008: prec = 3 + prec2;
! 2009: prec2= (long)nfnewprec(nf,-1);
! 2010: if (DEBUGLEVEL)
! 2011: fprintferr("init_chk: estimated prec = %ld (initially %ld)\n",
! 2012: prec, prec2);
! 2013: if (prec > prec2) return NULL;
! 2014: if (prec < prec2) chk_basis_embed = gprec_w(chk_basis_embed, prec);
! 2015: return (GEN)prec;
! 2016: }
! 2017:
! 2018: static GEN
! 2019: checkgenerator(GEN x)
! 2020: {
! 2021: long l,i,av = avma;
! 2022: GEN g = gmul(chk_basis_embed,x);
! 2023: GEN r = gneg((GEN)g[1]);
! 2024: long epsbit = 5 - bit_accuracy(gprecision(r));
! 2025:
! 2026: l = lg(g)-1;
! 2027: for (i=2; i<=l; i++)
! 2028: if (gexpo(gadd((GEN)g[i], r)) < epsbit) { avma=av; return NULL; }
! 2029:
! 2030: g = roots_to_pol_r1r2(g, chk_r1, chk_v);
! 2031: g = ground(g);
! 2032: if (lgef(modulargcd(g,derivpol(g))) > 3) { avma=av; return NULL; }
! 2033: if (DEBUGLEVEL>3) fprintferr(" generator: %Z\n",g);
! 2034: return g;
! 2035: }
! 2036:
! 2037: /* no garbage collecting, done in polredabs0 */
! 2038: static GEN
! 2039: findmindisc(GEN nf, GEN y, GEN a, GEN phimax, long flun)
! 2040: {
! 2041: long i,k, c = lg(y);
! 2042: GEN v,dmin,z,beta,discs = cgetg(c,t_VEC);
! 2043:
! 2044: for (i=1; i<c; i++) discs[i] = labsi(discsr((GEN)y[i]));
! 2045: v = sindexsort(discs);
! 2046: k = v[1]; dmin = (GEN)discs[k]; z = (GEN)y[k]; beta = (GEN)a[k];
! 2047: for (i=2; i<c; i++)
! 2048: {
! 2049: k = v[i];
! 2050: if (!egalii((GEN)discs[k],dmin)) break;
! 2051: if (gpolcomp((GEN)y[k],z) < 0) { z = (GEN)y[k]; beta = (GEN)a[k]; }
! 2052: }
! 2053: if (flun & nf_RAW)
! 2054: {
! 2055: y=cgetg(3,t_VEC);
! 2056: y[1]=lcopy(z);
! 2057: y[2]=lcopy(beta);
! 2058: }
! 2059: else if (phimax)
! 2060: {
! 2061: y=cgetg(3,t_VEC);
! 2062: y[1]=lcopy(z);
! 2063: beta=polymodrecip(gmodulcp(beta,(GEN)nf[1]));
! 2064: y[2]=(long)poleval(phimax,beta);
! 2065: }
! 2066: else y = gcopy(z);
! 2067: return y;
! 2068: }
! 2069:
! 2070: /* no garbage collecting, done in polredabs0 */
! 2071: static GEN
! 2072: storeallpols(GEN nf, GEN z, GEN a, GEN phimax, long flun)
! 2073: {
! 2074: GEN p1,y,beta;
! 2075:
! 2076: if (flun & nf_RAW)
! 2077: {
! 2078: long i, c = lg(z);
! 2079: y=cgetg(c,t_VEC);
! 2080: for (i=1; i<c; i++)
! 2081: {
! 2082: p1=cgetg(3,t_VEC); y[i]=(long)p1;
! 2083: p1[1]=lcopy((GEN)z[i]);
! 2084: p1[2]=lcopy((GEN)a[i]);
! 2085: }
! 2086: }
! 2087: else if (phimax)
! 2088: {
! 2089: long i, c = lg(z);
! 2090: beta = new_chunk(c);
! 2091: for (i=1; i<c; i++)
! 2092: beta[i] = (long)polymodrecip(gmodulcp((GEN)a[i],(GEN)nf[1]));
! 2093:
! 2094: y=cgetg(c,t_VEC);
! 2095: for (i=1; i<c; i++)
! 2096: {
! 2097: p1=cgetg(3,t_VEC); y[i]=(long)p1;
! 2098: p1[1]=lcopy((GEN)z[i]);
! 2099: p1[2]=(long)poleval(phimax,(GEN)beta[i]);
! 2100: }
! 2101: }
! 2102: else y = gcopy(z);
! 2103: return y;
! 2104: }
! 2105:
! 2106: GEN
! 2107: polredabs0(GEN x, long flun, long prec)
! 2108: {
! 2109: long i,nv, av = avma;
! 2110: GEN nf,v,y,a,phimax;
! 2111: GEN (*storepols)(ANYARG);
! 2112:
! 2113: if ((ulong)flun >= 16) err(flagerr,"polredabs");
! 2114: nf = initalgall0(x,nf_SMALL|nf_REGULAR,prec);
! 2115: if (lg(nf) == 3)
! 2116: {
! 2117: phimax = lift_to_pol((GEN)nf[2]);
! 2118: nf = (GEN)nf[1];
! 2119: }
! 2120: else
! 2121: phimax = (flun & nf_ORIG)? polx[0]: (GEN)NULL;
! 2122: prec = (long)nfnewprec(nf,0);
! 2123:
! 2124: nv = lgef(nf[1])==4? 1: 2;
! 2125: for (i=1; ; i++)
! 2126: {
! 2127: v = fincke_pohst(nf,NULL,stoi(5000),3,prec, &checkgenerator);
! 2128: if (v) break;
! 2129: if (i==MAXITERPOL) err(accurer,"polredabs0");
! 2130: prec = (prec<<1)-2; nf = nfnewprec(nf,prec);
! 2131: if (DEBUGLEVEL) err(warnprec,"polredabs0",prec);
! 2132: }
! 2133: a = (GEN)v[2];
! 2134: y = (GEN)v[1];
! 2135: nv = lg(a);
! 2136: for (i=1; i<nv; i++)
! 2137: if (canon_pol((GEN)y[i]) < 0 && phimax)
! 2138: a[i] = (long) gneg_i((GEN)a[i]);
! 2139: nv = remove_duplicates(y,a);
! 2140:
! 2141: if (DEBUGLEVEL)
! 2142: { fprintferr("%ld minimal vectors found.\n",nv-1); flusherr(); }
! 2143: if (nv >= 10000) flun &= (~nf_ALL); /* should not happen */
! 2144: storepols = (flun & nf_ALL)? storeallpols: findmindisc;
! 2145:
! 2146: if (DEBUGLEVEL) fprintferr("\n");
! 2147: if (nv==1)
! 2148: {
! 2149: x=(GEN)nf[1];
! 2150: y = cgetg(2,t_VEC); y[1]=(long)x;
! 2151: a = cgetg(2,t_VEC); a[1]=(long)polx[varn(x)];
! 2152: }
! 2153: return gerepileupto(av, storepols(nf,y,a,phimax,flun));
! 2154: }
! 2155:
! 2156: GEN
! 2157: polredabsall(GEN x, long flun, long prec)
! 2158: {
! 2159: return polredabs0(x, flun | nf_ALL, prec);
! 2160: }
! 2161:
! 2162: GEN
! 2163: polredabs(GEN x, long prec)
! 2164: {
! 2165: return polredabs0(x,nf_REGULAR,prec);
! 2166: }
! 2167:
! 2168: GEN
! 2169: polredabs2(GEN x, long prec)
! 2170: {
! 2171: return polredabs0(x,nf_ORIG,prec);
! 2172: }
! 2173:
! 2174: GEN
! 2175: polredabsnored(GEN x, long prec)
! 2176: {
! 2177: return polredabs0(x,nf_NORED,prec);
! 2178: }
! 2179:
! 2180: GEN
! 2181: polred(GEN x, long prec)
! 2182: {
! 2183: return allpolred(x,(GEN*)0,0,prec);
! 2184: }
! 2185:
! 2186: GEN
! 2187: smallpolred(GEN x, long prec)
! 2188: {
! 2189: return allpolred(x,(GEN*)0,1,prec);
! 2190: }
! 2191:
! 2192: GEN
! 2193: factoredpolred(GEN x, GEN p, long prec)
! 2194: {
! 2195: return allpolred(x,(GEN*)0,(long)p,prec);
! 2196: }
! 2197:
! 2198: GEN
! 2199: polred2(GEN x, long prec)
! 2200: {
! 2201: GEN y=cgetg(3,t_MAT);
! 2202:
! 2203: y[2]= (long) allpolred(x,(GEN*)(y+1),0,prec);
! 2204: return y;
! 2205: }
! 2206:
! 2207: GEN
! 2208: smallpolred2(GEN x, long prec)
! 2209: {
! 2210: GEN y=cgetg(3,t_MAT);
! 2211:
! 2212: y[2]= (long) allpolred(x,(GEN*)(y+1),1,prec);
! 2213: return y;
! 2214: }
! 2215:
! 2216: GEN
! 2217: factoredpolred2(GEN x, GEN p, long prec)
! 2218: {
! 2219: GEN y=cgetg(3,t_MAT);
! 2220:
! 2221: y[2]= (long) allpolred(x,(GEN*)(y+1),(long)p,prec);
! 2222: return y;
! 2223: }
! 2224:
! 2225: GEN makebasis(GEN nf,GEN pol);
! 2226: /* relative polredabs. Returns
! 2227: * flag = 0: relative polynomial
! 2228: * flag = 1: relative polynomial + element
! 2229: * flag = 2: absolute polynomial */
! 2230: GEN
! 2231: rnfpolredabs(GEN nf, GEN relpol, long flag, long prec)
! 2232: {
! 2233: GEN p1,bpol,rnf,elt,pol;
! 2234: long va, av = avma;
! 2235:
! 2236: if (typ(relpol)!=t_POL) err(typeer,"rnfpolredabs");
! 2237: nf=checknf(nf); va = varn(relpol);
! 2238: if (DEBUGLEVEL>1) timer2();
! 2239: p1 = makebasis(nf, unifpol(nf,relpol,1));
! 2240: rnf = (GEN)p1[3];
! 2241: if (DEBUGLEVEL>1)
! 2242: {
! 2243: msgtimer("absolute basis");
! 2244: fprintferr("original absolute generator: %Z\n",p1[1]);
! 2245: }
! 2246: p1 = polredabs0(p1, nf_RAW | nf_NORED, prec);
! 2247: bpol = (GEN)p1[1];
! 2248: if (DEBUGLEVEL>1) fprintferr("reduced absolute generator: %Z\n",bpol);
! 2249: if (flag==2) return gerepileupto(av,bpol);
! 2250:
! 2251: elt = rnfelementabstorel(rnf,(GEN)p1[2]);
! 2252: p1=cgetg(3,t_VEC);
! 2253: #if 1
! 2254: pol = rnfcharpoly(nf,relpol,elt,va);
! 2255: #else
! 2256: {
! 2257: long i;
! 2258: p1=(GEN)nffactor(nf,bpol)[1];
! 2259: for (i=lg(p1)-1; i; i--)
! 2260: if(gcmp0(gsubst((GEN)p1[i],va,elt))) { pol=(GEN)p1[i]; break; }
! 2261: if (!i) err(bugparier,"rnfpolredabs (pol not found)");
! 2262: }
! 2263: #endif
! 2264:
! 2265: if (!flag) p1 = pol;
! 2266: else
! 2267: {
! 2268: p1[1]=(long)pol;
! 2269: p1[2]=(long)polymodrecip(elt);
! 2270: }
! 2271: return gerepileupto(av,p1);
! 2272: }
! 2273:
! 2274: /********************************************************************/
! 2275: /** **/
! 2276: /** MINIM **/
! 2277: /** **/
! 2278: /********************************************************************/
! 2279: long addcolumntomatrix(long *V,long n,long r,GEN *INVP,long *L);
! 2280: GEN gmul_mat_smallvec(GEN x, GEN y, long hx, long ly);
! 2281:
! 2282: /* Minimal vectors for the integral definite quadratic form: a.
! 2283: * Result u:
! 2284: * u[1]= Number of vectors of square norm <= BORNE
! 2285: * u[2]= maximum norm found
! 2286: * u[3]= list of vectors found (at most STOCKMAX)
! 2287: *
! 2288: * If BORNE = gzero: Minimal non-zero vectors.
! 2289: * flag = min_ALL, as above
! 2290: * flag = min_FIRST, exits when first suitable vector is found.
! 2291: * flag = min_PERF, only compute rank of the family of v.v~ (v min.)
! 2292: */
! 2293: static GEN
! 2294: minim00(GEN a, GEN BORNE, GEN STOCKMAX, long flag)
! 2295: {
! 2296: GEN res,p1,u,r,liste,gnorme,gnorme_max,invp,V, *gptr[2];
! 2297: long n = lg(a), av0 = avma, av1,av,tetpil,lim, i,j,k,s,maxrank,*x;
! 2298: double p,borne,*v,*y,*z,**q, eps = 0.000001;
! 2299:
! 2300: switch(flag)
! 2301: {
! 2302: case min_FIRST: res = cgetg(3,t_VEC); break;
! 2303: case min_ALL: res = cgetg(4,t_VEC);
! 2304: }
! 2305: av=avma;
! 2306:
! 2307: x = (long*) new_chunk(n);
! 2308: q = (double**) new_chunk(n);
! 2309:
! 2310: /* correct alignment for the following */
! 2311: s = avma % sizeof(double); avma -= s;
! 2312: if (avma<bot) err(errpile);
! 2313:
! 2314: s = (n * sizeof(double))/sizeof(long);
! 2315: y = (double*) new_chunk(s);
! 2316: z = (double*) new_chunk(s);
! 2317: v = (double*) new_chunk(s);
! 2318: for (j=1; j<n; j++) q[j] = (double*) new_chunk(s);
! 2319: av1=avma;
! 2320:
! 2321: u = lllgramint(a); a = qf_base_change(a,u,1);
! 2322: if (lg(a) != n)
! 2323: err(talker,"not a definite form in minim00");
! 2324:
! 2325: n--;
! 2326: a = gmul(a, realun(DEFAULTPREC)); r = sqred1(a);
! 2327: if (DEBUGLEVEL>4) { fprintferr("minim: r = "); outerr(r); }
! 2328: for (j=1; j<=n; j++)
! 2329: {
! 2330: v[j] = rtodbl(gcoeff(r,j,j));
! 2331: for (i=1; i<j; i++) q[i][j] = rtodbl(gcoeff(r,i,j));
! 2332: }
! 2333:
! 2334: if (flag==min_PERF || gcmp0(BORNE))
! 2335: {
! 2336: double c, b = rtodbl(gcoeff(a,1,1));
! 2337:
! 2338: for (i=2; i<=n; i++)
! 2339: { c=rtodbl(gcoeff(a,i,i)); if (c<b) b=c; }
! 2340: borne = b+eps;
! 2341: BORNE = ground(dbltor(borne));
! 2342: gnorme_max = NULL;
! 2343: }
! 2344: else
! 2345: {
! 2346: BORNE = gfloor(BORNE);
! 2347: borne = gtodouble(BORNE)+eps;
! 2348: gnorme_max = gzero;
! 2349: }
! 2350:
! 2351: switch(flag)
! 2352: {
! 2353: case min_ALL:
! 2354: maxrank=itos(STOCKMAX);
! 2355: liste = new_chunk(1+maxrank);
! 2356: break;
! 2357: case min_PERF:
! 2358: BORNE = gerepileupto(av1,BORNE);
! 2359: maxrank = (n*(n+1))>>1;
! 2360: liste = new_chunk(1+maxrank); V = new_chunk(1+maxrank);
! 2361: for (i=1; i<=maxrank; i++) liste[i]=0;
! 2362: }
! 2363:
! 2364: s=0; av1=avma; lim = stack_lim(av1,1);
! 2365: k = n; y[n] = z[n] = 0;
! 2366: x[n] = (long) sqrt(borne/v[n]+eps);
! 2367: if (flag == min_PERF) invp = idmat(maxrank);
! 2368: for(;;)
! 2369: {
! 2370: do
! 2371: {
! 2372: if (k>1)
! 2373: {
! 2374: k--; z[k]=0;
! 2375: for (j=k+1; j<=n; j++) z[k] += q[k][j]*x[j];
! 2376: p = x[k+1]+z[k+1];
! 2377: y[k] = y[k+1] + p*p*v[k+1];
! 2378: x[k] = (long) floor(sqrt((borne-y[k]+eps)/v[k])-z[k]);
! 2379: }
! 2380: for(;;)
! 2381: {
! 2382: p=x[k]+z[k];
! 2383: if (y[k] + p*p*v[k] <= borne+eps) break;
! 2384: k++; x[k]--;
! 2385: }
! 2386: }
! 2387: while (k>1);
! 2388: if (! x[1] && y[1]<=eps) break;
! 2389: p = x[1]+z[1]; gnorme = ground( dbltor(y[1]+ p*p*v[1]) );
! 2390: if (gnorme_max)
! 2391: { if (gcmp(gnorme,gnorme_max) > 0) gnorme_max=gnorme; }
! 2392: else
! 2393: {
! 2394: if (gcmp(gnorme,BORNE) < 0)
! 2395: {
! 2396: borne=gtodouble(gnorme); s=0;
! 2397: affii(gnorme,BORNE); avma=av1;
! 2398: if (flag == min_PERF) invp = idmat(maxrank);
! 2399: }
! 2400: }
! 2401:
! 2402: switch(flag)
! 2403: {
! 2404: case min_ALL:
! 2405: s++;
! 2406: if (s<=maxrank)
! 2407: {
! 2408: p1 = new_chunk(n+1); liste[s] = (long) p1;
! 2409: for (i=1; i<=n; i++) p1[i] = x[i];
! 2410: }
! 2411: break;
! 2412:
! 2413: case min_FIRST:
! 2414: if (! gnorme_max || gcmp(gnorme,BORNE)>0) break;
! 2415:
! 2416: tetpil=avma; gnorme = icopy(gnorme); r = gmul_mat_smallvec(r,x,n,n);
! 2417: gptr[0]=&gnorme; gptr[1]=&r; gerepilemanysp(av,tetpil,gptr,2);
! 2418: res[1]=(long)gnorme;
! 2419: res[2]=(long)r; return res;
! 2420:
! 2421: case min_PERF:
! 2422: {
! 2423: long av2=avma, I=1, newran;
! 2424:
! 2425: for (i=1; i<=n; i++)
! 2426: for (j=i; j<=n; j++,I++) V[I] = x[i]*x[j];
! 2427: newran = addcolumntomatrix(V,maxrank,s,&invp,liste);
! 2428: if (newran == s)
! 2429: {
! 2430: avma=av2;
! 2431: if (DEBUGLEVEL>1) { fprintferr("."); flusherr(); }
! 2432: }
! 2433: else
! 2434: {
! 2435: if (DEBUGLEVEL>1) { fprintferr("*"); flusherr(); }
! 2436: s = newran;
! 2437: if (s == maxrank)
! 2438: {
! 2439: if (DEBUGLEVEL>1) { fprintferr("\n"); flusherr(); }
! 2440: avma=av0; return stoi(s);
! 2441: }
! 2442:
! 2443: if (low_stack(lim, stack_lim(av1,1)))
! 2444: {
! 2445: if(DEBUGMEM>1) err(warnmem,"minim00");
! 2446: if (DEBUGLEVEL>1)
! 2447: {
! 2448: fprintferr("\ngerepile in qfperfection. rank>=%ld\n",s);
! 2449: flusherr();
! 2450: }
! 2451: tetpil=avma; invp = gerepile(av1,tetpil,gcopy(invp));
! 2452: }
! 2453: }
! 2454: }
! 2455: }
! 2456: x[1]--;
! 2457: }
! 2458: switch(flag)
! 2459: {
! 2460: case min_FIRST:
! 2461: avma=av0; return cgetg(1,t_VEC);
! 2462: case min_PERF:
! 2463: if (DEBUGLEVEL>1) { fprintferr("\n"); flusherr(); }
! 2464: avma=av0; return stoi(s);
! 2465: }
! 2466: k = min(s,maxrank);
! 2467:
! 2468: tetpil = avma; p1=cgetg(k+1,t_MAT);
! 2469: for (j=1; j<=k; j++)
! 2470: p1[j] = (long) gmul_mat_smallvec(u,(GEN)liste[j],n,n);
! 2471: liste = p1;
! 2472: r = gnorme_max? gnorme_max: BORNE;
! 2473:
! 2474: r=icopy(r); gptr[0]=&r; gptr[1]=&liste;
! 2475: gerepilemanysp(av,tetpil,gptr,2);
! 2476: res[1]=lstoi(s<<1);
! 2477: res[2]=(long)r;
! 2478: res[3]=(long)liste; return res;
! 2479: }
! 2480:
! 2481: GEN
! 2482: minim0(GEN a, GEN borne, GEN stockmax, long flag, long prec)
! 2483: {
! 2484: switch(flag)
! 2485: {
! 2486: case 0: return minim00(a,borne,stockmax,min_ALL);
! 2487: case 1: return minim00(a,borne,gzero ,min_FIRST);
! 2488: case 2: return fincke_pohst(a,borne,stockmax,0,prec,NULL);
! 2489: default: err(flagerr,"qfminim");
! 2490: }
! 2491: return NULL; /* not reached */
! 2492: }
! 2493:
! 2494: GEN
! 2495: minim(GEN a, GEN borne, GEN stockmax)
! 2496: {
! 2497: return minim00(a,borne,stockmax,min_ALL);
! 2498: }
! 2499:
! 2500: GEN
! 2501: minim2(GEN a, GEN borne, GEN stockmax)
! 2502: {
! 2503: return minim00(a,borne,stockmax,min_FIRST);
! 2504: }
! 2505:
! 2506: GEN
! 2507: perf(GEN a)
! 2508: {
! 2509: return minim00(a,gzero,gzero,min_PERF);
! 2510: }
! 2511:
! 2512: /* programme general pour les formes quadratiques definies positives
! 2513: * quelconques (a coeffs reels). On doit avoir BORNE != 0; on n'entre dans
! 2514: * cette fonction qu'a partir de fincke_pohst (la reduction LLL n'est donc
! 2515: * pas faite ici). Si flag >= 2, on s'arrete des que stockmax est atteint.
! 2516: * Si flag&1 == 1, pas d'erreur dans sqred1intern. */
! 2517: static GEN
! 2518: smallvectors(GEN a, GEN BORNE, GEN STOCKMAX, long flag, long prec,
! 2519: GEN (*check)(GEN))
! 2520: {
! 2521: long av = avma,av1,av2,lim,N,n,i,j,k,s,stockmax,epsbit,checkcnt = 0;
! 2522: GEN u,S,x,y,z,v,q,norme1,normax1,borne1,borne2,eps,p1,alpha,norms;
! 2523:
! 2524: N=lg(a); n=N-1; stockmax=itos(STOCKMAX);
! 2525: if (DEBUGLEVEL)
! 2526: fprintferr("smallvectors looking for norm <= %Z\n",gprec_w(BORNE,3));
! 2527:
! 2528: lim=stack_lim(av,1); q=sqred1intern(a,flag&1);
! 2529: if (q == NULL) { avma=av; return NULL; }
! 2530: if (DEBUGLEVEL>5) fprintferr("q = %Z",q);
! 2531: epsbit = bit_accuracy(prec) >> 1;
! 2532: eps=realun(prec); setexpo(eps,-epsbit);
! 2533: alpha = dbltor(0.95);
! 2534: normax1=gzero;
! 2535: borne1=gadd(BORNE,eps);
! 2536: borne2=mpmul(borne1,alpha);
! 2537: v=cgetg(N,t_VEC);
! 2538:
! 2539: av2=avma;
! 2540: S = cgetg(stockmax+1,t_MAT);
! 2541: norms = cgetg(stockmax+1,t_VEC);
! 2542: x=cgetg(N,t_COL);
! 2543: y=cgetg(N,t_COL);
! 2544: z=cgetg(N,t_COL);
! 2545: for (i=1; i<N; i++) { v[i] = coeff(q,i,i); x[i]=y[i]=z[i] = zero; }
! 2546:
! 2547: x[n] = lmpent(mpsqrt(gdiv(borne1,(GEN)v[n])));
! 2548: if (DEBUGLEVEL>3) { fprintferr("\nx[%ld] = %Z\n",n,x[n]); flusherr(); }
! 2549:
! 2550: s=0; k=n;
! 2551: for(;;)
! 2552: {
! 2553: do
! 2554: {
! 2555: int fl=0;
! 2556: if (k>1)
! 2557: {
! 2558: av1=avma; k--;
! 2559: p1 = mpmul(gcoeff(q,k,k+1),(GEN)x[k+1]);
! 2560: for (j=k+2; j<N; j++)
! 2561: p1 = mpadd(p1, mpmul(gcoeff(q,k,j),(GEN)x[j]));
! 2562: z[k] = lpileupto(av1,p1);
! 2563:
! 2564: av1=avma;
! 2565: p1 = gsqr(mpadd((GEN)x[k+1],(GEN)z[k+1]));
! 2566: y[k] = lpileupto(av1, mpadd((GEN)y[k+1], mpmul(p1,(GEN)v[k+1])));
! 2567:
! 2568: av1=avma;
! 2569: p1 = mpsub(borne1, (GEN)y[k]);
! 2570: if (signe(p1) < 0) { avma=av1; fl = 1; }
! 2571: else
! 2572: {
! 2573: p1 = mpadd(eps,mpsub(mpsqrt(gdiv(p1,(GEN)v[k])), (GEN)z[k]));
! 2574: avma=av1; x[k] = lmpent(p1); /* safe */
! 2575: }
! 2576: }
! 2577: else if (check)
! 2578: { /* don't waste time on the [x,0,...0] */
! 2579: for (i=k+1; i<N; i++)
! 2580: if (signe(x[i])) break;
! 2581: if (i == N && signe(x[1])) x[1] = un;
! 2582: }
! 2583: for(;;)
! 2584: {
! 2585: av1=avma;
! 2586: if (!fl)
! 2587: {
! 2588: p1 = mpmul((GEN)v[k], gsqr(mpadd((GEN)x[k],(GEN)z[k])));
! 2589: if (mpcmp(mpsub(mpadd(p1,(GEN)y[k]), borne1),
! 2590: gmul2n(mpabs(p1),-epsbit)) <= 0) { avma=av1; break; }
! 2591: }
! 2592: avma=av1; k++; fl=0;
! 2593: x[k]=laddis((GEN)x[k],-1);
! 2594: }
! 2595: if (low_stack(lim, stack_lim(av,1)))
! 2596: {
! 2597: GEN *gptr[8];
! 2598: int a = 4;
! 2599: if(DEBUGMEM>1) err(warnmem,"smallvectors");
! 2600: gptr[0]=&x; gptr[1]=&y; gptr[2]=&z; gptr[3]=&normax1;
! 2601: if (stockmax)
! 2602: {
! 2603: for (i=s+1; i<=stockmax; i++) S[i]=zero;
! 2604: gptr[4]=&S; a++;
! 2605: }
! 2606: if (check)
! 2607: {
! 2608: gptr[5]=&borne1; gptr[6]=&borne2; gptr[7]=&norms;
! 2609: for (i=s+1; i<=stockmax; i++) norms[i]=zero;
! 2610: a+=3;
! 2611: }
! 2612: gerepilemany(av2,gptr,a);
! 2613: }
! 2614: if (DEBUGLEVEL>3)
! 2615: {
! 2616: if (DEBUGLEVEL>5) fprintferr("%ld ",k);
! 2617: if (k==n) fprintferr("\nx[%ld] = %Z\n",n,x[n]);
! 2618: flusherr();
! 2619: }
! 2620: }
! 2621: while (k>1);
! 2622: if (!signe(x[1]) && gexpo((GEN)y[1]) <= -epsbit) break;
! 2623:
! 2624: av1=avma; p1 = gsqr(mpadd((GEN)x[1],(GEN)z[1]));
! 2625: norme1 = mpadd((GEN)y[1], mpmul(p1, (GEN)v[1]));
! 2626: if (mpcmp(norme1,borne1) > 0) avma=av1;
! 2627: else
! 2628: {
! 2629: norme1 = gerepileupto(av1,norme1);
! 2630: if (check)
! 2631: {
! 2632: if (checkcnt < 5 && mpcmp(norme1, borne2) < 0)
! 2633: {
! 2634: if (check(x))
! 2635: {
! 2636: borne1 = mpadd(norme1,eps);
! 2637: borne2 = mpmul(borne1,alpha);
! 2638: s = 0; checkcnt = 0;
! 2639: }
! 2640: else { checkcnt++ ; goto CONTINUE; }
! 2641: }
! 2642: }
! 2643: else if (mpcmp(norme1,normax1) > 0) normax1=norme1;
! 2644: if (++s <= stockmax)
! 2645: {
! 2646: norms[s] = (long)norme1;
! 2647: p1 = cgetg(N,t_COL); S[s] = (long)p1;
! 2648: for (i=1; i<N; i++) p1[i]=x[i];
! 2649: if (s == stockmax && (flag&2) && check)
! 2650: {
! 2651: long av1 = avma;
! 2652: GEN per = sindexsort(norms);
! 2653: if (DEBUGLEVEL) fprintferr("sorting...\n");
! 2654: for (i=1; i<=s; i++)
! 2655: {
! 2656: long k = per[i];
! 2657: if (check((GEN)S[k]))
! 2658: {
! 2659: S[1] = S[k]; avma = av1;
! 2660: borne1 = mpadd(norme1,eps);
! 2661: borne2 = mpmul(borne1,alpha);
! 2662: s = 1; checkcnt = 0; break;
! 2663: }
! 2664: }
! 2665: if (checkcnt) s = 0;
! 2666: }
! 2667: }
! 2668: }
! 2669: CONTINUE:
! 2670: x[k] = laddis((GEN)x[k],-1);
! 2671: }
! 2672: if (s<stockmax) stockmax = s;
! 2673: if (check)
! 2674: {
! 2675: GEN per, alph,pols,p;
! 2676: if (DEBUGLEVEL) fprintferr("final sort & check...\n");
! 2677: setlg(norms,s+1); per = sindexsort(norms);
! 2678: alph = cgetg(s+1,t_VEC);
! 2679: pols = cgetg(s+1,t_VEC);
! 2680: for (j=0,i=1; i<=s; i++)
! 2681: {
! 2682: long k = per[i];
! 2683: if (j && mpcmp((GEN)norms[k], borne1) > 0) break;
! 2684: if ((p = check((GEN)S[k])))
! 2685: {
! 2686: if (!j) borne1 = gadd((GEN)norms[k],eps);
! 2687: j++; pols[j]=(long)p; alph[j]=S[k];
! 2688: }
! 2689: }
! 2690: u = cgetg(3,t_VEC);
! 2691: setlg(pols,j+1); u[1] = (long)pols;
! 2692: setlg(alph,j+1); u[2] = (long)alph; return u;
! 2693: }
! 2694: u=cgetg(4,t_VEC);
! 2695: u[1]=lstoi(s<<1);
! 2696: u[2]=(long)normax1; setlg(S,stockmax+1);
! 2697: u[3]=(long)S; return u;
! 2698: }
! 2699:
! 2700: /* return T2 norm of the polynomial defining nf */
! 2701: static GEN
! 2702: get_Bnf(GEN nf)
! 2703: {
! 2704: GEN p = gzero, r = (GEN)nf[6];
! 2705: long i, r1 = itos(gmael(nf,2,1)), ru = lg(r)-1;
! 2706: for (i=ru; i>0; i--)
! 2707: {
! 2708: if (i == r1) p = gmul2n(p, 1);
! 2709: p = gadd(p, gnorm((GEN)r[i]));
! 2710: }
! 2711: if (i == r1) p = gmul2n(p, 1);
! 2712: return p;
! 2713: }
! 2714:
! 2715: /* solve x~.a.x <= borne, a > 0. If check is non-NULL keep x only if check(x).
! 2716: * flag = 1, return NULL in case of precision problems (sqred1 or lllgram)
! 2717: * raise an error otherwse.
! 2718: * flag = 2, return as soon as stockmax vectors found.
! 2719: * flag = 3, corresponds to 1+2 */
! 2720: GEN
! 2721: fincke_pohst(GEN a,GEN bound,GEN stockmax,long flag, long prec,
! 2722: GEN (*check)(GEN))
! 2723: {
! 2724: long pr,av=avma,i,j,n;
! 2725: GEN B,nf,r,rinvtrans,v,v1,u,s,res,z,vnorm,sperm,perm,uperm,basis,gram;
! 2726:
! 2727: if (DEBUGLEVEL>2) { fprintferr("entering fincke_pohst\n"); flusherr(); }
! 2728: if (typ(a) == t_VEC) { nf = a; a = gmael(nf,5,3); } else nf = NULL;
! 2729: pr = gprecision(a);
! 2730: if (pr) prec = pr; else a = gmul(a,realun(prec));
! 2731: if (DEBUGLEVEL>2) fprintferr("first LLL: prec = %ld\n", prec);
! 2732: v1 = lllgramintern(a,4,flag&1, (prec<<1)-2);
! 2733: if (v1 == NULL) goto PRECPB;
! 2734: r = qf_base_change(a,v1,1);
! 2735: r = sqred1intern(r,flag&1);
! 2736: if (r == NULL) goto PRECPB;
! 2737:
! 2738: n = lg(a);
! 2739: for (i=1; i<n; i++)
! 2740: {
! 2741: GEN p1 = gsqrt(gcoeff(r,i,i), prec);
! 2742: coeff(r,i,i)=(long)p1;
! 2743: for (j=i+1; j<n; j++)
! 2744: coeff(r,i,j) = lmpmul(p1, gcoeff(r,i,j));
! 2745: }
! 2746: /* now r~ * r = a in LLL basis */
! 2747: rinvtrans = gtrans(invmat(r));
! 2748: if (DEBUGLEVEL>2)
! 2749: fprintferr("final LLL: prec = %ld, precision(rinvtrans) = %ld\n",
! 2750: prec,gprecision(rinvtrans));
! 2751: v = lllintern(rinvtrans,flag&1, (gprecision(rinvtrans)<<1)-2);
! 2752: if (v == NULL) goto PRECPB;
! 2753: rinvtrans = gmul(rinvtrans,v);
! 2754:
! 2755: u = invmat(gtrans(v)); s = gmul(r,u);
! 2756: u = gmul(v1, u);
! 2757: vnorm=cgetg(n,t_VEC);
! 2758: for (j=1; j<n; j++) vnorm[j] = lnorml2((GEN)rinvtrans[j]);
! 2759: perm = sindexsort(vnorm);
! 2760: sperm = cgetg(n,t_MAT);
! 2761: uperm = cgetg(n,t_MAT);
! 2762: for (i=1; i<n; i++) { uperm[n-i] = u[perm[i]]; sperm[n-i] = s[perm[i]]; }
! 2763:
! 2764: gram = gram_matrix(sperm);
! 2765: B = gcoeff(gram,n-1,n-1);
! 2766: if (gexpo(B) >= bit_accuracy(lg(B)-2)) goto PRECPB;
! 2767:
! 2768: if (check && nf) basis = init_chk(nf,uperm,NULL);
! 2769: if (!bound)
! 2770: { /* polred */
! 2771: GEN x = cgetg(n,t_COL);
! 2772:
! 2773: if (nf) bound = get_Bnf(nf);
! 2774: for (j=1; j<n; j++) x[j]=zero;
! 2775: for (j=2; j<n; j++)
! 2776: {
! 2777: x[j]=un; B = gcoeff(gram,j,j);
! 2778: if (!bound || mpcmp(B,bound) < 0)
! 2779: if (!check || check(x)) bound = B;
! 2780: x[j]=zero;
! 2781: }
! 2782: }
! 2783:
! 2784: if (DEBUGLEVEL>2) {fprintferr("entering smallvectors\n"); flusherr();}
! 2785: if (check && nf)
! 2786: if (! (prec = (long)init_chk(nf,uperm,bound))) goto PRECPB;
! 2787: i = check? 2: 1; if (i == n) i--;
! 2788: for ( ; i<n; i++)
! 2789: {
! 2790: res = smallvectors(gram,bound? bound: gceil(gcoeff(gram,i,i)),
! 2791: stockmax,flag,prec,check);
! 2792: if (!res) goto PRECPB;
! 2793: if (!check || lg(res[2]) > 1) break;
! 2794: if (DEBUGLEVEL>2) fprintferr(" i = %ld failed\n",i);
! 2795: }
! 2796: if (check)
! 2797: {
! 2798: GEN p1 = (GEN)res[2];
! 2799: for (i=1; i<lg(p1); i++) p1[i] = lmul(basis, (GEN)p1[i]);
! 2800: return res;
! 2801: }
! 2802:
! 2803: if (DEBUGLEVEL>2) {fprintferr("leaving fincke_pohst\n"); flusherr();}
! 2804: z=cgetg(4,t_VEC);
! 2805: z[1]=lcopy((GEN)res[1]);
! 2806: z[2]=pr? lcopy((GEN)res[2]) : lround((GEN)res[2]);
! 2807: z[3]=lmul(uperm, (GEN)res[3]); return gerepileupto(av,z);
! 2808: PRECPB:
! 2809: if (!(flag & 1))
! 2810: err(talker,"not a positive definite form in fincke_pohst");
! 2811: avma = av; return NULL;
! 2812: }
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 / src / basemath