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