OpenXM_contrib/pari-2.2/src/basemath/bibli1.c - diff

Return to bibli1.c CVS log

Up to [local] / OpenXM_contrib / pari-2.2 / src / basemath

Diff for /OpenXM_contrib/pari-2.2/src/basemath/Attic/bibli1.c between version 1.1 and 1.2

-version 1.1, 2001/10/02 11:17:02
+version 1.2, 2002/09/11 07:26:49
 Line 21  Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 Line 21  Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 Line 21  Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  #include "pari.h"
  #include "parinf.h"
  extern GEN ZV_lincomb(GEN u, GEN v, GEN X, GEN Y);
- extern GEN roots_to_pol_r1r2(GEN a, long r1, long v);
+ extern int addcolumntomatrix(GEN V,GEN INVP,GEN L);
- extern GEN makebasis(GEN nf,GEN pol);
+ extern long expodb(double x);
- extern GEN caractducos(GEN p, GEN x, int v);
+ /* default quality ratio for LLL: 99/100 */
+ #define LLLDFT 100
  /* scalar product x.x */
  GEN
  sqscal(GEN x)
  {
-   long i,av,lx;
+   long i, lx;
+   gpmem_t av;
    GEN z;
    lx = lg(x);
    if (lx == 1) return gzero;
-Line 44  sqscal(GEN x)
+Line 47  sqscal(GEN x)
 Line 44  sqscal(GEN x)
 Line 47  sqscal(GEN x)
  GEN
  gscal(GEN x,GEN y)
  {
-   long i,av,lx;
+   long i, lx;
+   gpmem_t av;
    GEN z;
    if (x == y) return sqscal(x);
    lx = lg(x);
-Line 59  gscal(GEN x,GEN y)
+Line 63  gscal(GEN x,GEN y)
 Line 59  gscal(GEN x,GEN y)
 Line 63  gscal(GEN x,GEN y)
  static GEN
  sqscali(GEN x)
  {
-   long i,av,lx;
+   long i, lx;
+   gpmem_t av;
    GEN z;
    lx = lg(x);
    if (lx == 1) return gzero;
-Line 73  sqscali(GEN x)
+Line 78  sqscali(GEN x)
 Line 73  sqscali(GEN x)
 Line 78  sqscali(GEN x)
  static GEN
  gscali(GEN x,GEN y)
  {
-   long i,av,lx;
+   long i, lx;
+   gpmem_t av;
    GEN z;
    if (x == y) return sqscali(x);
    lx = lg(x);
-Line 85  gscali(GEN x,GEN y)
+Line 91  gscali(GEN x,GEN y)
 Line 85  gscali(GEN x,GEN y)
 Line 91  gscali(GEN x,GEN y)
    return gerepileuptoint(av,z);
  }
+ /********************************************************************/
+ /**             QR Factorization via Householder matrices          **/
+ /********************************************************************/
+ /* zero x[1..k-1], fill mu */
+ static int
+ FindApplyQ(GEN x, GEN mu, GEN B, long k, GEN Q, long prec)
+ {
+   long i, lx = lg(x)-1, lv = lx - (k-1);
+   GEN v, p1, beta, Nx, x2, x1, xd = x + (k-1);
+   x1 = (GEN)xd[1];
+   x2 = gsqr(x1);
+   if (k < lx)
+   {
+     for (i=2; i<=lv; i++) x2 = mpadd(x2, gsqr((GEN)xd[i]));
+     Nx = gsqrt(x2, prec);
+     if (signe(x1) < 0) setsigne(Nx, -1);
+     v = cgetg(lv+1, t_VEC);
+     v[1] = lmpadd(x1, Nx);
+     for (i=2; i<=lv; i++) v[i] = xd[i];
+     if (gcmp0(x2)) return 0;
+     if (gcmp0(x1))
+       beta = mpmul(x2, realun(prec)); /* make sure typ(beta) != t_INT */
+     else
+       beta = mpadd(x2, mpmul(Nx,x1));
+     beta = ginv(beta);
+     p1 = cgetg(3,t_VEC); Q[k] = (long)p1;
+     p1[1] = (long)beta;
+     p1[2] = (long)v;
+     coeff(mu,k,k) = lmpneg(Nx);
+   }
+   else
+     coeff(mu,k,k) = x[k];
+   if (B)
+   {
+     B[k] = (long)x2;
+     for (i=1; i<k; i++) coeff(mu,k,i) = x[i];
+   }
+   else
+     for (i=1; i<k; i++) coeff(mu,i,k) = x[i];
+   return (typ(x2) != t_REAL || lg(x2) >  3 || expo(x2) < 10);
+ }
+ static void
+ ApplyQ(GEN Q, GEN r)
+ {
+   GEN s, rd, beta = (GEN)Q[1], v = (GEN)Q[2];
+   long i, l = lg(v), lr = lg(r);
+   rd = r + (lr - l);
+   s = mpmul((GEN)v[1], (GEN)rd[1]);
+   for (i=2; i<l; i++) s = mpadd(s, mpmul((GEN)v[i] ,(GEN)rd[i]));
+   s = mpneg(mpmul(beta, s));
+   for (i=1; i<l; i++) rd[i] = lmpadd((GEN)rd[i], mpmul(s, (GEN)v[i]));
+ }
  static GEN
- lllall_trivial(GEN x, long n, long flag)
+ ApplyAllQ(GEN Q, GEN r0, long k)
  {
-   GEN y;
+   gpmem_t av = avma;
-   if (!n)
+   GEN r = dummycopy(r0);
+   long j;
+   for (j=1; j<k; j++) ApplyQ((GEN)Q[j], r);
+   return gerepilecopy(av, r);
+ }
+ /* compute B[k] = | x[k] |^2, update mu(k, 1..k-1) using Householder matrices
+  * (Q = Householder(x[1..k-1]) in factored form) */
+ static int
+ incrementalQ(GEN x, GEN L, GEN B, GEN Q, long k, long prec)
+ {
+   GEN r = ApplyAllQ(Q, (GEN)x[k], k);
+   return FindApplyQ(r, L, B, k, Q, prec);
+ }
+ /* Q vector of Householder matrices orthogonalizing x[1..j0].
+  * Q[i] = 0 means not computed yet */
+ static int
+ Householder_get_mu(GEN x, GEN L, GEN B, long k, GEN Q, long prec)
+ {
+   GEN Nx, invNx, m;
+   long i, j, j0;
+   if (!Q) Q = zerovec(k);
+   for (j=1; j<=k; j++)
+     if (typ(Q[j]) == t_INT) break;
+   j0 = j;
+   for (   ; j<=k; j++)
+     if (! incrementalQ(x, L, B, Q, j, prec)) return 0;
+   for (j=1; j<=k; j++)
    {
-     if (flag != lll_ALL) return cgetg(1,t_MAT);
+     m = (GEN)L[j]; Nx = (GEN)m[j]; /* should set m[j] = un; but need it later */
+     if (j == k) break;
+     invNx = ginv(Nx);
+     for (i=max(j0, j+1); i<=k; i++) m[i] = lmpmul(invNx, (GEN)m[i]);
+   }
+   return 1;
+ }
+ GEN
+ sqred1_from_QR(GEN x, long prec)
+ {
+   long j, k = lg(x)-1;
+   GEN L, B = zerovec(k);
+   L = cgetg(k+1, t_MAT);
+   for (j=1; j<=k; j++) L[j] = (long)zerocol(k);
+   if (!Householder_get_mu(x, L, B, k, NULL, prec)) return NULL;
+   for (j=1; j<=k; j++) coeff(L,j,j) = B[j];
+   return gtrans_i(L);
+ }
+ GEN
+ R_from_QR(GEN x, long prec)
+ {
+   long j, k = lg(x)-1;
+   GEN L, B = zerovec(k), Q = cgetg(k+1, t_VEC);
+   L = cgetg(k+1, t_MAT);
+   for (j=1; j<=k; j++) L[j] = (long)zerocol(k);
+   for (j=1; j<=k; j++)
+     if (!incrementalQ(x, L, B, Q, j, prec)) return NULL;
+   return gtrans_i(L);
+ }
+ /********************************************************************/
+ /**             QR Factorization via Gram-Schmidt                  **/
+ /********************************************************************/
+ /* compute B[k] = | x[k] |^2, update mu(k, 1..k-1).
+  * Classical Gram-Schmidt (unstable!) */
+ static int
+ incrementalGS(GEN x, GEN mu, GEN B, long k)
+ {
+   GEN s,A = cgetg(k+1, t_COL); /* scratch space */
+   long i, j;
+   gpmem_t av;
+   A[1] = coeff(x,k,1);
+   for(j=1;j<k;)
+   {
+     coeff(mu,k,j) = lmpdiv((GEN)A[j], (GEN)B[j]);
+     j++; av = avma;
+     /* A[j] <-- x[k,j] - sum_{i<j} mu[j,i] A[i] */
+     s = mpmul(gcoeff(mu,j,1),(GEN)A[1]);
+     for (i=2; i<j; i++) s = mpadd(s, mpmul(gcoeff(mu,j,i),(GEN)A[i]));
+     s = mpneg(s); A[j] = (long)gerepileuptoleaf(av, mpadd(gcoeff(x,k,j), s));
+   }
+   B[k] = A[k]; return (signe((GEN)B[k]) > 0);
+ }
+ #if 0
+ /* return Gram-Schmidt orthogonal basis (f) associated to (e), B is the
+  * vector of the (f_i . f_i)*/
+ GEN
+ gram_schmidt(GEN e, GEN *ptB)
+ {
+   long i,j,lx = lg(e);
+   GEN f = dummycopy(e), B, iB;
+   B = cgetg(lx, t_VEC);
+   iB= cgetg(lx, t_VEC);
+   for (i=1;i<lx;i++)
+   {
+     GEN p1 = NULL;
+     gpmem_t av;
+     B[i] = (long)sqscal((GEN)f[i]);
+     iB[i]= linv((GEN)B[i]); av = avma;
+     for (j=1; j<i; j++)
+     {
+       GEN mu = gmul(gscal((GEN)e[i],(GEN)f[j]), (GEN)iB[j]);
+       GEN p2 = gmul(mu, (GEN)f[j]);
+       p1 = p1? gadd(p1,p2): p2;
+     }
+     p1 = p1? gerepileupto(av, gsub((GEN)e[i], p1)): (GEN)e[i];
+     f[i] = (long)p1;
+   }
+   *ptB = B; return f;
+ }
+ #endif
+ /********************************************************************/
+ /**                                                                **/
+ /**                          LLL ALGORITHM                         **/
+ /**                                                                **/
+ /********************************************************************/
+ #define swap(x,y) { long _t=x; x=y; y=_t; }
+ #define gswap(x,y) { GEN _t=x; x=y; y=_t; }
+ static GEN
+ lll_trivial(GEN x, long flag)
+ {
+   GEN y;
+   if (lg(x) == 1)
+   { /* dim x = 0 */
+     if (! (flag & lll_ALL)) return cgetg(1,t_MAT);
      y=cgetg(3,t_VEC);
      y[1]=lgetg(1,t_MAT);
      y[2]=lgetg(1,t_MAT); return y;
    }
-   /* here n = 1 */
+   /* here dim = 1 */
    if (gcmp0((GEN)x[1]))
    {
-     switch(flag ^ lll_GRAM)
+     switch(flag & (~lll_GRAM))
      {
        case lll_KER: return idmat(1);
        case lll_IM : return cgetg(1,t_MAT);
-Line 120  lllall_trivial(GEN x, long n, long flag)
+Line 315  lllall_trivial(GEN x, long n, long flag)
 Line 120  lllall_trivial(GEN x, long n, long flag)
 Line 315  lllall_trivial(GEN x, long n, long flag)
  }
  static GEN
- lllgramall_finish(GEN h,GEN fl,long flag,long n)
+ lll_finish(GEN h,GEN fl,long flag)
  {
-   long k;
+   long i, k, l = lg(fl);
-   GEN y;
+   GEN y, g;
-   k=1; while (k<=n && !fl[k]) k++;
+   k=1; while (k<l && !fl[k]) k++;
-   switch(flag)
+   switch(flag & (~lll_GRAM))
    {
      case lll_KER: setlg(h,k);
-       y = gcopy(h); break;
+       return h;
-     case lll_IM: h += k-1; h[0] = evaltyp(t_MAT)| evallg(n-k+2);
+     case lll_IM: h += k-1; h[0] = evaltyp(t_MAT) | evallg(l-k+1);
-       y = gcopy(h); break;
+       return h;
-     default: setlg(h,k); y=cgetg(3,t_VEC);
-       y[1] = lcopy(h);
-       h += k-1; h[0] = evaltyp(t_MAT)| evallg(n-k+2);
-       y[2] = lcopy(h);
-       break;
    }
-   return y;
+   y = cgetg(3,t_VEC);
+   g = cgetg(k, t_MAT); for (i=1; i<k; i++) g[i] = h[i];
+   y[1] = (long)g;
+   h += k-1; h[0] = evaltyp(t_MAT) | evallg(l-k+1);
+   y[2] = (long)h; return y;
  }
- /********************************************************************/
+ /* h[,k] += q * h[,l] */
- /**                                                                **/
+ static void
- /**                          LLL with content                      **/
+ Zupdate_col(long k, long l, GEN q, long K, GEN h)
- /**                                                                **/
- /********************************************************************/
- /* real Gram matrix has coeffs X[i,j] = x[i,j]*veccon[i]*veccon[j] */
- static GEN
- lllgramintwithcontent(GEN x, GEN veccon, long flag)
  {
-   long av0=avma,av,tetpil,lx=lg(x),i,j,k,l,n,lim,kmax;
+   GEN *hl, *hk;
-   GEN u,u2,B,lam,q,r,h,la,bb,p1,p2,p3,p4,fl,corr,corr2,newcon;
+   long i;
-   GEN *gptr[8];
-   if (typ(x) != t_MAT) err(typeer,"lllgramintwithcontent");
+   if (!h) return;
-   n=lx-1; if (n<=1) return lllall_trivial(x,n,flag);
+   hl = (GEN*)h[l]; hk = (GEN*)h[k];
-   if (lg((GEN)x[1])!=lx) err(mattype1,"lllgramintwithcontent");
+   if (is_pm1(q))
-   fl = new_chunk(lx);
+   {
+     if (signe(q) > 0)
+       for (i=1;i<=K;i++) { if (signe(hl[i])) hk[i] = addii(hk[i],hl[i]); }
+     else
+       for (i=1;i<=K;i++) { if (signe(hl[i])) hk[i] = subii(hk[i],hl[i]); }
+   } else
+       for (i=1;i<=K;i++) if (signe(hl[i])) hk[i] = addii(hk[i],mulii(q,hl[i]));
+ }
-   av=avma; lim=stack_lim(av,1);
+ /* L[k,] += q * L[l,], l < k */
-   x=dummycopy(x); veccon=dummycopy(veccon);
+ static void
-   B=cgetg(lx+1,t_COL); B[1]=un;
+ Zupdate_row(long k, long l, GEN q, GEN L, GEN B)
-   /* B[i+1]=B_i; le vrai B_i est B_i*prod(1,j=1,i,veccon[j]^2) */
+ {
+   long i;
-   for (i=1; i<lx; i++) { B[i+1]=zero; fl[i]=0; }
+   if (is_pm1(q))
-   lam=cgetg(lx,t_MAT);
-   for (j=1; j<lx; j++)
-   { p2=cgetg(lx,t_COL); lam[j]=(long)p2; for (i=1; i<lx; i++) p2[i]=zero; }
- /* le vrai lam[i,j] est
-    lam[i,j]*veccon[i]*veccon[j]*prod(1,l=1,j-1,veccon[l]^2) */
-   k=2; h=idmat(n); kmax=1;
-   u=gcoeff(x,1,1); if (typ(u)!=t_INT) err(lllger4);
-   if (signe(u)) { B[2]=(long)u; coeff(lam,1,1)=un; fl[1]=1; }
-   else { B[2]=un; fl[1]=0; }
-   if (DEBUGLEVEL>5) { fprintferr("k = %ld",k); flusherr(); }
-   for(;;)
    {
-     if (k>kmax)
+     if (signe(q) > 0)
      {
-       kmax=k;
+       for (i=1;i<l; i++) coeff(L,k,i) = laddii(gcoeff(L,k,i),gcoeff(L,l,i));
-       for (j=1; j<=k; j++)
+       coeff(L,k,l) = laddii(gcoeff(L,k,l), B);
-       {
+     } else {
-         if (j==k || fl[j])
+       for (i=1;i<l; i++) coeff(L,k,i) = lsubii(gcoeff(L,k,i),gcoeff(L,l,i));
-         {
+       coeff(L,k,l) = laddii(gcoeff(L,k,l), negi(B));
-           u=gcoeff(x,k,j); if (typ(u)!=t_INT) err(lllger4);
-           for (i=1; i<j; i++)
-             if (fl[i])
-               u=divii(subii(mulii((GEN)B[i+1],u),mulii(gcoeff(lam,k,i),gcoeff(lam,j,i))),(GEN)B[i]);
-           if (j<k) coeff(lam,k,j)=(long)u;
-           else
-           {
-             if (signe(u)) { B[k+1]=(long)u; coeff(lam,k,k)=un; fl[k]=1; }
-             else { B[k+1]=B[k]; fl[k]=0; }
-           }
-         }
-       }
-       if (low_stack(lim, stack_lim(av,1)))
-       {
-         if(DEBUGMEM>1) err(warnmem,"[1]: lllgramintwithcontent");
-         gptr[0]=&B; gptr[1]=&lam; gptr[2]=&h;
-         gptr[3]=&x; gptr[4]=&veccon; gerepilemany(av,gptr,5);
-       }
      }
-     u=shifti(mulii(gcoeff(lam,k,k-1),(GEN)veccon[k]),1);
+   } else {
-     u2=mulii((GEN)B[k],(GEN)veccon[k-1]);
+     for(i=1;i<l;i++)  coeff(L,k,i)=laddii(gcoeff(L,k,i),mulii(q,gcoeff(L,l,i)));
-     if (cmpii(absi(u),u2)>0)
+     coeff(L,k,l) = laddii(gcoeff(L,k,l), mulii(q,B));
-     {
-       q=dvmdii(addii(u,u2),shifti(u2,1),&r);
-       if (signe(r)<0) q=addsi(-1,q);
-       h[k]=lsub((GEN)h[k],gmul(q,(GEN)h[k-1]));
-       newcon=mppgcd((GEN)veccon[k],(GEN)veccon[k-1]);
-       corr=divii((GEN)veccon[k],newcon); veccon[k]=(long)newcon;
-       if(!gcmp1(corr))
-       {
-         corr2=sqri(corr);
-         for (j=1; j<=n; j++)
-           coeff(x,j,k)=coeff(x,k,j)=lmulii(corr,gcoeff(x,k,j));
-         coeff(x,k,k)=lmulii(corr,gcoeff(x,k,k));
-         for (j=k; j<=kmax; j++) B[j+1]=lmulii(corr2,(GEN)B[j+1]);
-         for (i=1; i<k; i++) coeff(lam,k,i)=lmulii(corr,gcoeff(lam,k,i));
-         for (i=k+1; i<=kmax; i++)
-         {
-           coeff(lam,i,k)=lmulii(corr,gcoeff(lam,i,k));
-           for (j=k+1; j<i; j++)
-             coeff(lam,i,j)=lmulii(corr2,gcoeff(lam,i,j));
-         }
-       }
-       r=negi(mulii(q,divii((GEN)veccon[k-1],(GEN)veccon[k])));
-       p1=gcoeff(x,k,k-1);
-       for (j=1; j<=n; j++)
-         coeff(x,j,k)=coeff(x,k,j)=laddii(gcoeff(x,j,k),mulii(r,(j==k)?p1:gcoeff(x,j,k-1)));
-       coeff(x,k,k)=laddii(gcoeff(x,k,k),mulii(r,gcoeff(x,k-1,k)));
-       coeff(lam,k,k-1)=laddii(gcoeff(lam,k,k-1),mulii(r,(GEN)B[k]));
-       for (i=1; i<k-1; i++)
-         coeff(lam,k,i)=laddii(gcoeff(lam,k,i),mulii(r,gcoeff(lam,k-1,i)));
-     }
-     if (low_stack(lim, stack_lim(av,1)))
-     {
-       if(DEBUGMEM>1) err(warnmem,"[2]: lllgramintwithcontent");
-       gptr[0]=&B; gptr[1]=&lam; gptr[2]=&h;
-       gptr[3]=&x; gptr[4]=&veccon; gerepilemany(av,gptr,5);
-     }
-     p3=mulii((GEN)B[k-1],(GEN)B[k+1]);la=gcoeff(lam,k,k-1);p4=mulii(la,la);
-     p2=mulsi(100,mulii(mulii((GEN)veccon[k],(GEN)veccon[k]),addii(p3,p4)));
-     p1=mulii((GEN)veccon[k-1],(GEN)B[k]);p1=mulsi(99,mulii(p1,p1));
-     if (fl[k-1] && (cmpii(p1,p2)>0 || !fl[k]))
-     {
-       if (DEBUGLEVEL>=4 && k==n)
-         { fprintferr(" (%ld)", expi(p1)-expi(p2)); flusherr(); }
-       p1=(GEN)h[k-1]; h[k-1]=h[k]; h[k]=(long)p1;
-       p1=(GEN)x[k-1]; x[k-1]=x[k]; x[k]=(long)p1;
-       for (j=1; j<=n; j++)
-       { p1=gcoeff(x,k-1,j); coeff(x,k-1,j)=coeff(x,k,j); coeff(x,k,j)=(long)p1; }
-       p1=(GEN)veccon[k-1]; veccon[k-1]=veccon[k]; veccon[k]=(long)p1;
-       for (j=1; j<=k-2; j++)
-       { p1=gcoeff(lam,k-1,j); coeff(lam,k-1,j)=coeff(lam,k,j); coeff(lam,k,j)=(long)p1; }
-       if (fl[k])
-       {
-         for (i=k+1; i<=kmax; i++)
-         {
-           bb=gcoeff(lam,i,k);
-           coeff(lam,i,k)=ldivii(subii(mulii((GEN)B[k+1],gcoeff(lam,i,k-1)),mulii(la,bb)),(GEN)B[k]);
-           coeff(lam,i,k-1)=ldivii(addii(mulii(la,gcoeff(lam,i,k-1)),mulii((GEN)B[k-1],bb)),(GEN)B[k]);
-           if (low_stack(lim, stack_lim(av,1)))
-           {
-             if(DEBUGMEM>1) err(warnmem,"[3]: lllgramintwithcontent");
-             gptr[0]=&B; gptr[1]=&lam; gptr[2]=&h;
-             gptr[3]=&x; gptr[4]=&la; gptr[5]=&veccon; gptr[6]=&p3;
-             gptr[7]=&p4; gerepilemany(av,gptr,8);
-           }
-         }
-         B[k]=ldivii(addii(p3,p4),(GEN)B[k]);
-       }
-       else
-       {
-         if (signe(la))
-         {
-           p2=(GEN)B[k]; p1=divii(p4,p2);
-           for (i=k+1; i<=kmax; i++)
-             coeff(lam,i,k-1)=ldivii(mulii(la,gcoeff(lam,i,k-1)),p2);
-           for (j=k+1; j<kmax; j++)
-           {
-             for (i=j+1; i<=kmax; i++)
-               coeff(lam,i,j)=ldivii(mulii(p1,gcoeff(lam,i,j)),p2);
-             if (low_stack(lim, stack_lim(av,1)))
-             {
-               if(DEBUGMEM>1) err(warnmem,"[4]: lllgramintwithcontent");
-               gptr[0]=&B; gptr[1]=&lam; gptr[2]=&h;
-               gptr[3]=&x; gptr[4]=&la; gptr[5]=&veccon; gptr[6]=&p1;
-               gptr[7]=&p2; gerepilemany(av,gptr,8);
-             }
-           }
-           B[k+1]=B[k]=(long)p1;
-           for (i=k+2; i<=lx; i++)
-             B[i]=ldivii(mulii(p1,(GEN)B[i]),p2);
-         }
-         else
-         {
-           coeff(lam,k,k-1)=zero;
-           for (i=k+1; i<=kmax; i++)
-           {
-             coeff(lam,i,k)=coeff(lam,i,k-1);
-             coeff(lam,i,k-1)=zero;
-           }
-           B[k]=B[k-1]; fl[k]=1; fl[k-1]=0;
-         }
-         if (low_stack(lim, stack_lim(av,1)))
-         {
-           if(DEBUGMEM>1) err(warnmem,"[5]: lllgramintwithcontent");
-           gptr[0]=&B; gptr[1]=&lam; gptr[2]=&h;
-           gptr[3]=&x; gptr[4]=&veccon;
-           gerepilemany(av,gptr,5);
-         }
-       }
-       if (k>2) k--;
-       if (DEBUGLEVEL>5) { fprintferr(" %ld",k); flusherr(); }
-     }
-     else
-     {
-       for (l=k-2; l>=1; l--)
-       {
-         u=shifti(mulii(gcoeff(lam,k,l),(GEN)veccon[k]),1);
-         u2=mulii((GEN)B[l+1],(GEN)veccon[l]);
-         if (cmpii(absi(u),u2)>0)
-         {
-           q=dvmdii(addii(u,u2),shifti(u2,1),&r);
-           if (signe(r)<0) q=addsi(-1,q);
-           h[k]=lsub((GEN)h[k],gmul(q,(GEN)h[l]));
-           newcon=mppgcd((GEN)veccon[k],(GEN)veccon[l]);
-           corr=divii((GEN)veccon[k],newcon); veccon[k]=(long)newcon;
-           if(!gcmp1(corr))
-           {
-             corr2=sqri(corr);
-             for (j=1; j<=n; j++)
-               coeff(x,j,k)=coeff(x,k,j)=lmulii(corr,gcoeff(x,k,j));
-             coeff(x,k,k)=lmulii(corr,gcoeff(x,k,k));
-             for (j=k; j<=kmax; j++) B[j+1]=lmulii(corr2,(GEN)B[j+1]);
-             for (i=1; i<k; i++) coeff(lam,k,i)=lmulii(corr,gcoeff(lam,k,i));
-             for (i=k+1; i<=kmax; i++)
-             {
-               coeff(lam,i,k)=lmulii(corr,gcoeff(lam,i,k));
-               for (j=k+1; j<i; j++)
-                 coeff(lam,i,j)=lmulii(corr2,gcoeff(lam,i,j));
-             }
-           }
-           r=negi(mulii(q,divii((GEN)veccon[l],(GEN)veccon[k])));
-           p1=gcoeff(x,k,l);
-           for (j=1; j<=n; j++)
-             coeff(x,j,k)=coeff(x,k,j)=laddii(gcoeff(x,j,k),mulii(r,(j==k)?p1:gcoeff(x,j,l)));
-           coeff(x,k,k)=laddii(gcoeff(x,k,k),mulii(r,gcoeff(x,l,k)));
-           coeff(lam,k,l)=laddii(gcoeff(lam,k,l),mulii(r,(GEN)B[l+1]));
-           for (i=1; i<l; i++)
-             coeff(lam,k,i)=laddii(gcoeff(lam,k,i),mulii(r,gcoeff(lam,l,i)));
-         }
-         if (low_stack(lim, stack_lim(av,1)))
-         {
-           if(DEBUGMEM>1) err(warnmem,"[6]: lllgramintwithcontent");
-           gptr[0]=&B; gptr[1]=&lam; gptr[2]=&h;
-           gptr[3]=&x; gptr[4]=&veccon; gerepilemany(av,gptr,5);
-         }
-       }
-       k++; if (DEBUGLEVEL>5) { fprintferr(" %ld",k); flusherr(); }
-       if (k>n)
-       {
-         if (DEBUGLEVEL>5) { fprintferr("\n"); flusherr(); }
-         tetpil=avma;
-         return gerepile(av0,tetpil,lllgramall_finish(h,fl,flag,n));
-       }
-     }
-     if (low_stack(lim, stack_lim(av,1)))
-     {
-       if(DEBUGMEM>1) err(warnmem,"[7]: lllgramintwithcontent");
-       gptr[0]=&B; gptr[1]=&lam; gptr[2]=&h;
-       gptr[3]=&x; gptr[4]=&veccon; gerepilemany(av,gptr,5);
-     }
    }
  }
- static GEN
+ static void
- lllintwithcontent(GEN x)
+ update_row(long k, long l, GEN q, GEN L)
  {
-   long lx=lg(x),i,j,av,tetpil;
+   long i;
-   GEN veccon,con,xred,g;
+   if (is_pm1(q))
-   if (typ(x) != t_MAT) err(typeer,"lllintwithcontent");
-   if (lx==1) return cgetg(1,t_MAT);
-   av=avma; veccon=cgetg(lx,t_VEC); g=cgetg(lx,t_MAT); xred=cgetg(lx,t_MAT);
-   for (j=1; j<lx; j++)
    {
-     g[j]=lgetg(lx,t_COL); con=content((GEN)x[j]);
+     if (signe(q) > 0)
-     xred[j]=ldiv((GEN)x[j],con); veccon[j]=(long)con;
+     {
+       for (i=1;i<l; i++) coeff(L,k,i) = lmpadd(gcoeff(L,k,i),gcoeff(L,l,i));
+     } else {
+       for (i=1;i<l; i++) coeff(L,k,i) = lmpsub(gcoeff(L,k,i),gcoeff(L,l,i));
+     }
+   } else {
+     for(i=1;i<l;i++)  coeff(L,k,i)=lmpadd(gcoeff(L,k,i),mpmul(q,gcoeff(L,l,i)));
    }
-   for (i=1; i<lx; i++)
+   coeff(L,k,l) = lmpadd(gcoeff(L,k,l), q);
-     for (j=1; j<=i; j++)
-       coeff(g,i,j)=coeff(g,j,i)=(long)gscal((GEN)xred[i],(GEN)xred[j]);
-   tetpil=avma;
-   return gerepile(av,tetpil,lllgramintwithcontent(g,veccon,2));
  }
- /********************************************************************/
- /**                                                                **/
- /**                          LLL ALGORITHM                         **/
- /**                                                                **/
- /********************************************************************/
- #define swap(x,y) { long _t=x; x=y; y=_t; }
- #define gswap(x,y) { GEN _t=x; x=y; y=_t; }
  static void
- REDI(GEN x, GEN h, GEN L, GEN B, long K, long k, long l)
+ ZRED_gram(long k, long l, GEN x, GEN h, GEN L, GEN B, long K)
  {
    long i,lx;
    GEN q = truedvmdii(addii(B,shifti(gcoeff(L,k,l),1)), shifti(B,1), NULL);
-   GEN *hk,*hl,xk,xl;
+   GEN xk,xl;
    if (!signe(q)) return;
-   q = negi(q); lx = lg(x);
+   q = negi(q);
-   xl = (GEN)x[l]; hl = (GEN*)h[l];
+   xl = (GEN) x[l]; xk = (GEN) x[k];
-   xk = (GEN)x[k]; hk = (GEN*)h[k];
+   lx = lg(xl);
    if (is_pm1(q))
    {
      if (signe(q) > 0)
      {
-       for (i=1;i<=K;i++) hk[i]=addii(hk[i],hl[i]);
        xk[k] = laddii((GEN)xk[k], (GEN)xl[k]);
-       for (i=1;i<lx;i++) coeff(x,k,i)=xk[i]=laddii((GEN)xk[i], (GEN)xl[i]);
+       for (i=1;i<lx;i++) coeff(x,k,i) = xk[i] = laddii((GEN)xk[i], (GEN)xl[i]);
-       for (i=1;i<l; i++) coeff(L,k,i)=laddii(gcoeff(L,k,i),gcoeff(L,l,i));
-       q = B;
      } else {
-       for (i=1;i<=K;i++) hk[i]=subii(hk[i],hl[i]);
        xk[k] = lsubii((GEN)xk[k], (GEN)xl[k]);
-       for (i=1;i<lx;i++) coeff(x,k,i)=xk[i]=lsubii((GEN)xk[i], (GEN)xl[i]);
+       for (i=1;i<lx;i++) coeff(x,k,i) = xk[i] = lsubii((GEN)xk[i], (GEN)xl[i]);
-       for (i=1;i<l; i++) coeff(L,k,i)=lsubii(gcoeff(L,k,i),gcoeff(L,l,i));
-       q = negi(B);
      }
-   } else {
+   } else { /* h[,k] += q* h[,l]. x[,k] += q * x[,l]. L[k,] += q* L[l,] */
-     for(i=1;i<=K;i++) hk[i]=addii(hk[i],mulii(q,hl[i]));
      xk[k] = laddii((GEN)xk[k], mulii(q,(GEN)xl[k]));
      for(i=1;i<lx;i++) coeff(x,k,i)=xk[i]=laddii((GEN)xk[i],mulii(q,(GEN)xl[i]));
-     for(i=1;i<l;i++)  coeff(L,k,i)=laddii(gcoeff(L,k,i),mulii(q,gcoeff(L,l,i)));
-     q = mulii(q,B);
    }
-   coeff(L,k,l) = laddii(gcoeff(L,k,l), q);
+   Zupdate_row(k,l,q,L,B);
+   Zupdate_col (k,l,q,K,h);
  }
- /* b[k] <-- b[k] - round(L[k,l]) b[l], only b[1] ... b[K] modified so far
-  * assume l < k and update x = Gram(b), L = Gram Schmidt coeffs. */
  static void
- RED(GEN x, GEN h, GEN L, long K, long k, long l)
+ ZRED(long k, long l, GEN x, GEN h, GEN L, GEN B, long K)
  {
-   long e,i,lx;
+   GEN q = truedvmdii(addii(B,shifti(gcoeff(L,k,l),1)), shifti(B,1), NULL);
-   GEN q = grndtoi(gcoeff(L,k,l),&e);
-   GEN *hk,*hl,xk,xl;
-   if (DEBUGLEVEL>8)
-     fprintferr("error bits when rounding in lllgram: %ld\n",e);
    if (!signe(q)) return;
+   q = negi(q);
+   Zupdate_row(k,l,q,L,B);
+   Zupdate_col (k,l,q,K,h);
+   x[k] = (long)ZV_lincomb(gun, q, (GEN)x[k], (GEN)x[l]);
+ }
+ static GEN
+ round_safe(GEN q)
+ {
+   long e;
+   if (typ(q) == t_INT) return q;
+   if (expo(q) > 30)
+   {
+     q = grndtoi(q, &e);
+     if (e > 0) return NULL;
+   } else
+     q = ground(q);
+   return q;
+ }
+ /* b[k] <-- b[k] - round(L[k,l]) b[l], only b[1] ... b[K] modified so far
+  * assume l < k and update x = Gram(b), L = Gram Schmidt coeffs. */
+ static int
+ RED_gram(long k, long l, GEN x, GEN h, GEN L, long K)
+ {
+   long i,lx;
+   GEN q = round_safe(gcoeff(L,k,l));
+   GEN xk, xl;
+   if (!q) return 0;
+   if (!signe(q)) return 1;
    q = negi(q); lx = lg(x);
-   xl = (GEN)x[l]; hl = (GEN*)h[l];
+   xk = (GEN)x[k]; xl = (GEN)x[l];
-   xk = (GEN)x[k]; hk = (GEN*)h[k];
    if (is_pm1(q))
    {
      if (signe(q) > 0)
      {
-       for (i=1;i<=K;i++) hk[i]=addii(hk[i],hl[i]);
+       xk[k] = lmpadd((GEN)xk[k], (GEN)xl[k]);
-       xk[k] = ladd((GEN)xk[k], (GEN)xl[k]);
+       for (i=1;i<lx;i++) coeff(x,k,i)=xk[i]=lmpadd((GEN)xk[i], (GEN)xl[i]);
-       for (i=1;i<lx;i++) coeff(x,k,i)=xk[i]=ladd((GEN)xk[i], (GEN)xl[i]);
-       for (i=1;i<l; i++) coeff(L,k,i)=ladd(gcoeff(L,k,i),gcoeff(L,l,i));
      } else {
-       for (i=1;i<=K;i++) hk[i]=subii(hk[i],hl[i]);
+       xk[k] = lmpsub((GEN)xk[k], (GEN)xl[k]);
-       xk[k] = lsub((GEN)xk[k], (GEN)xl[k]);
+       for (i=1;i<lx;i++) coeff(x,k,i)=xk[i]=lmpsub((GEN)xk[i], (GEN)xl[i]);
-       for (i=1;i<lx;i++) coeff(x,k,i)=xk[i]=lsub((GEN)xk[i], (GEN)xl[i]);
-       for (i=1;i<l; i++) coeff(L,k,i)=lsub(gcoeff(L,k,i),gcoeff(L,l,i));
      }
    } else {
-     for (i=1;i<=K;i++) hk[i]=addii(hk[i],mulii(q,hl[i]));
+     xk[k] = lmpadd((GEN)xk[k], mpmul(q,(GEN)xl[k]));
-     xk[k] = ladd((GEN)xk[k], gmul(q,(GEN)xl[k]));
+     for (i=1;i<lx;i++)
-     for (i=1;i<lx;i++) coeff(x,k,i)=xk[i]=ladd((GEN)xk[i], gmul(q,(GEN)xl[i]));
+       coeff(x,k,i)=xk[i]=lmpadd((GEN)xk[i], mpmul(q,(GEN)xl[i]));
-     for (i=1;i<l; i++) coeff(L,k,i)=ladd(gcoeff(L,k,i),gmul(q,gcoeff(L,l,i)));
    }
-   coeff(L,k,l) = ladd(gcoeff(L,k,l),q);
+   update_row(k,l,q,L);
+   Zupdate_col(k,l,q,K,h); return 1;
  }
+ static void
+ update_col(long k, long l, GEN q, GEN x)
+ {
+   long i,lx;
+   GEN xk = (GEN)x[k], xl = (GEN)x[l];
+   lx = lg(xk);
+   if (is_pm1(q))
+   {
+     if (signe(q) > 0)
+     {
+       for (i=1;i<lx;i++) xk[i] = lmpadd((GEN)xk[i], (GEN)xl[i]);
+     } else {
+       for (i=1;i<lx;i++) xk[i] = lmpsub((GEN)xk[i], (GEN)xl[i]);
+     }
+   } else {
+     for (i=1;i<lx;i++) xk[i] = lmpadd((GEN)xk[i], mpmul(q,(GEN)xl[i]));
+   }
+ }
  static int
- do_SWAPI(GEN x, GEN h, GEN L, GEN B, long kmax, long k, long alpha, GEN fl)
+ RED(long k, long l, GEN x, GEN h, GEN L, long K)
  {
-   GEN la,la2,p1,p2,Bk;
+   GEN q = round_safe(gcoeff(L,k,l));
-   long av,i,j,lx;
+   if (!q) return 0;
+   if (!signe(q)) return 1;
+   q = negi(q);
+   update_col(k,l,q,x);
+   update_row(k,l,q,L);
+   Zupdate_col(k,l,q,K,h); return 1;
+ }
+ static int
+ do_ZSWAP(GEN x, GEN h, GEN L, GEN B, long kmax, long k, long D, GEN fl,
+          int gram)
+ {
+   GEN la,la2,p1,Bk;
+   long i, j, lx;
+   gpmem_t av;
    if (!fl[k-1]) return 0;
-   lx = lg(x); av = avma;
+   av = avma;
    la = gcoeff(L,k,k-1); la2 = sqri(la);
-   p2 = addii(la2, mulii((GEN)B[k-1],(GEN)B[k+1]));
    Bk = (GEN)B[k];
-   if (fl[k] && cmpii(mulsi(alpha-1,sqri(Bk)), mulsi(alpha,p2)) <= 0)
+   if (fl[k])
-     { avma = av; return 0; }
+   {
+     GEN q;
-   /* SWAPI(k-1,k) */
+     if (!D) return 0; /* only swap non-kernel + kernel vector */
-   if (DEBUGLEVEL>3 && k==kmax) {
+     q = addii(la2, mulii((GEN)B[k-1],(GEN)B[k+1]));
-     fprintferr(" (%ld)", expi(mulsi(alpha-1,sqri(Bk)))
+     if (cmpii(mulsi(D-1,sqri(Bk)), mulsi(D,q)) <= 0) {
-                        - expi(mulsi(alpha,p2))); flusherr();
+       avma = av; return 0;
+     }
+     B[k] = (long)diviiexact(q, Bk);
    }
-   swap(h[k-1], h[k]);
+   /* ZSWAP(k-1,k) */
-   swap(x[k-1], x[k]);
+   if (DEBUGLEVEL>3 && k==kmax)
-   for (j=1; j< lx; j++) swap(coeff(x,k-1,j), coeff(x,k,j));
+   { /* output diagnostics associated to re-normalized rational quantities */
+     gpmem_t av1 = avma;
+     GEN d = mulii((GEN)B[k-1],(GEN)B[k+1]);
+     p1 = subii(mulsi(D-1, sqri(Bk)), mulsi(D, la2));
+     fprintferr(" (%ld)", expi(p1) - expi(mulsi(D, d)));
+     avma = av1;
+   }
+   if (h) swap(h[k-1], h[k]);
+   swap(x[k-1], x[k]); lx = lg(x);
+   if (gram)
+     for (j=1; j < lx; j++) swap(coeff(x,k-1,j), coeff(x,k,j));
    for (j=1; j<k-1; j++) swap(coeff(L,k-1,j), coeff(L,k,j))
    if (fl[k])
    {
-Line 514  do_SWAPI(GEN x, GEN h, GEN L, GEN B, long kmax, long k
+Line 553  do_SWAPI(GEN x, GEN h, GEN L, GEN B, long kmax, long k
 Line 514  do_SWAPI(GEN x, GEN h, GEN L, GEN B, long kmax, long k
 Line 553  do_SWAPI(GEN x, GEN h, GEN L, GEN B, long kmax, long k
      {
        GEN t = gcoeff(L,i,k);
        p1 = subii(mulii((GEN)B[k+1],gcoeff(L,i,k-1)), mulii(la,t));
-       p1 = divii(p1, Bk);
+       p1 = diviiexact(p1, Bk);
-       av = avma = coeff(L,i,k) = (long)icopy_av(p1,(GEN)av);
+       coeff(L,i,k) = (long)icopy_av(p1,(GEN)av);
+       av = avma = (gpmem_t)coeff(L,i,k);
        p1 = addii(mulii(la,gcoeff(L,i,k-1)), mulii((GEN)B[k-1],t));
-       p1 = divii(p1, Bk);
+       p1 = diviiexact(p1, Bk);
-       av = avma = coeff(L,i,k-1) = (long)icopy_av(p1,(GEN)av);
+       coeff(L,i,k-1) = (long)icopy_av(p1,(GEN)av);
+       av = avma = (gpmem_t)coeff(L,i,k-1);
      }
-     B[k] = ldivii(p2,Bk);
    }
+   else if (signe(la))
+   {
+     p1 = diviiexact(la2, Bk);
+     B[k+1] = B[k] = (long)p1;
+     for (i=k+2; i<=lx; i++) B[i] = (long)diviiexact(mulii(p1,(GEN)B[i]), Bk);
+     for (i=k+1; i<=kmax; i++)
+       coeff(L,i,k-1) = (long)diviiexact(mulii(la,gcoeff(L,i,k-1)), Bk);
+     for (j=k+1; j<kmax; j++)
+       for (i=j+1; i<=kmax; i++)
+         coeff(L,i,j) = (long)diviiexact(mulii(p1,gcoeff(L,i,j)), Bk);
+   }
    else
    {
-     if (signe(la))
+     for (i=k+1; i<=kmax; i++)
      {
-       p1 = divii(la2, Bk);
+       coeff(L,i,k) = coeff(L,i,k-1);
-       B[k+1] = B[k] = (long)p1;
+       coeff(L,i,k-1) = zero;
-       for (i=k+2; i<=lx; i++) B[i] = ldivii(mulii(p1,(GEN)B[i]), Bk);
-       for (i=k+1; i<=kmax; i++)
-         coeff(L,i,k-1) = ldivii(mulii(la,gcoeff(L,i,k-1)), Bk);
-       for (j=k+1; j<kmax; j++)
-         for (i=j+1; i<=kmax; i++)
-           coeff(L,i,j) = ldivii(mulii(p1,gcoeff(L,i,j)), Bk);
      }
-     else
+     B[k] = B[k-1]; fl[k] = 1; fl[k-1] = 0;
-     {
-       for (i=k+1; i<=kmax; i++)
-       {
-         coeff(L,i,k) = coeff(L,i,k-1);
-         coeff(L,i,k-1) = zero;
-       }
-       B[k] = B[k-1]; fl[k] = 1; fl[k-1] = 0;
-     }
    }
    return 1;
  }
  static int
- do_SWAP(GEN x, GEN h, GEN L, GEN B, long kmax, long k, GEN QR)
+ do_SWAP(GEN x, GEN h, GEN L, GEN B, long kmax, long k, GEN delta, int gram)
  {
    GEN la,la2, BK,BB,q;
-   long av,i,j,lx;
+   long i, j, lx;
+   gpmem_t av;
-   lx = lg(x); av = avma;
+   av = avma;
    la = gcoeff(L,k,k-1); la2 = gsqr(la);
-   q = gmul((GEN)B[k-1], gsub(QR,la2));
+   q = mpmul((GEN)B[k-1], mpsub(delta,la2));
-   if (gcmp(q, (GEN)B[k]) <= 0) { avma = av; return 0; }
+   if (mpcmp(q, (GEN)B[k]) <= 0) { avma = av; return 0; }
    /* SWAP(k-1,k) */
    if (DEBUGLEVEL>3 && k==kmax) {
      fprintferr(" (%ld)", gexpo(q) - gexpo((GEN)B[k])); flusherr();
    }
-   BB = gadd((GEN)B[k], gmul((GEN)B[k-1],la2));
+   BB = mpadd((GEN)B[k], mpmul((GEN)B[k-1],la2));
-   if (gcmp0(BB)) { B[k] = 0; return 1; } /* precision pb */
+   if (!signe(BB)) { B[k] = 0; return 1; } /* precision pb */
-   coeff(L,k,k-1) = ldiv(gmul(la,(GEN)B[k-1]), BB);
+   coeff(L,k,k-1) = lmpdiv(mpmul(la,(GEN)B[k-1]), BB);
-   BK = gdiv((GEN)B[k], BB);
+   BK = mpdiv((GEN)B[k], BB);
-   B[k] = lmul((GEN)B[k-1], BK);
+   B[k] = lmpmul((GEN)B[k-1], BK);
    B[k-1] = (long)BB;
    swap(h[k-1],h[k]);
-   swap(x[k-1],x[k]);
+   swap(x[k-1],x[k]); lx = lg(x);
-   for (j=1; j<lx ; j++) swap(coeff(x,k-1,j), coeff(x,k,j));
+   if (gram)
+     for (j=1; j < lx; j++) swap(coeff(x,k-1,j), coeff(x,k,j));
    for (j=1; j<k-1; j++) swap(coeff(L,k-1,j), coeff(L,k,j))
    for (i=k+1; i<=kmax; i++)
    {
      GEN t = gcoeff(L,i,k);
-     coeff(L,i,k) = lsub(gcoeff(L,i,k-1), gmul(la,t));
+     coeff(L,i,k) = lmpsub(gcoeff(L,i,k-1), mpmul(la,t));
-     coeff(L,i,k-1) = ladd(gmul(gcoeff(L,k,k-1),gcoeff(L,i,k-1)), gmul(BK,t));
+     coeff(L,i,k-1) = lmpadd(t, mpmul(gcoeff(L,k,k-1), gcoeff(L,i,k)));
-    /*              = ladd(t, gmul(gcoeff(L,k,k-1), gcoeff(L,i,k)));
-     * would save one multiplication, but loses precision faster... */
    }
    return 1;
  }
+ static void
+ ZincrementalGS(GEN x, GEN L, GEN B, long k, GEN fl, int gram)
+ {
+   GEN u = NULL; /* gcc -Wall */
+   long i, j, s;
+   for (j=1; j<=k; j++)
+     if (j==k || fl[j])
+     {
+       gpmem_t av = avma;
+       u = gram? gcoeff(x,k,j): gscali((GEN)x[k], (GEN)x[j]);
+       for (i=1; i<j; i++)
+         if (fl[i])
+         {
+           u = subii(mulii((GEN)B[i+1], u), mulii(gcoeff(L,k,i), gcoeff(L,j,i)));
+           u = diviiexact(u, (GEN)B[i]);
+         }
+       coeff(L,k,j) = lpileuptoint(av, u);
+     }
+   s = signe(u);
+   if (s == 0) B[k+1] = B[k];
+   else
+   {
+     if (s < 0) err(lllger3);
+     B[k+1] = coeff(L,k,k); coeff(L,k,k) = un; fl[k] = 1;
+   }
+ }
- /* x integer matrix */
+ /* x integer matrix. Beware: this function can return NULL */
  GEN
- lllgramall(GEN x, long alpha, long flag)
+ lllint_marked(long MARKED, GEN x, long D, int gram,
+               GEN *pth, GEN *ptfl, GEN *ptB)
  {
-   long av0=avma,av,tetpil,lim,lx=lg(x),i,j,k,l,n,s,kmax;
+   long lx = lg(x), hx, i, j, k, l, n, kmax;
-   GEN u,B,L,h,fl, *gptr[4];
+   gpmem_t av, lim;
+   GEN B,L,h,fl;
-   if (typ(x) != t_MAT) err(typeer,"lllgramall");
+   if (typ(x) != t_MAT) err(typeer,"lllint");
-   n=lx-1; if (n<=1) return lllall_trivial(x,n,flag);
+   n = lx-1; if (n <= 1) return NULL;
-   if (lg((GEN)x[1])!=lx) err(mattype1,"lllgramall");
+   hx = lg(x[1]);
+   if (gram && hx != lx) err(mattype1,"lllint");
    fl = cgetg(lx,t_VECSMALL);
-   av=avma; lim=stack_lim(av,1); x=dummycopy(x);
+   av = avma; lim = stack_lim(av,1); x = dummycopy(x);
-   B=gscalcol(gun, lx);
+   B = gscalcol(gun, lx);
-   L=cgetg(lx,t_MAT);
+   L = cgetg(lx,t_MAT);
    for (j=1; j<lx; j++)
    {
-     for (i=1; i<lx; i++)
+     for (i=1; i<hx; i++)
-       if (typ(gcoeff(x,i,j))!=t_INT) err(lllger4);
+       if (typ(gcoeff(x,i,j)) != t_INT) err(lllger4);
      fl[j] = 0; L[j] = (long)zerocol(n);
    }
-   k=2; h=idmat(n); kmax=1;
+   h = pth? idmat(n): NULL;
-   u=gcoeff(x,1,1); s= signe(u);
+   ZincrementalGS(x, L, B, 1, fl, gram);
-   if (s == 0) B[2]=un;
+   kmax = 1;
-   else
+   if (DEBUGLEVEL>5) fprintferr("k = ");
+   for (k=2;;)
    {
-     if (s < 0) err(lllger3);
-     B[2]=(long)u; coeff(L,1,1)=un; fl[1]=1;
-   }
-   if (DEBUGLEVEL>5) fprintferr("k =");
-   for(;;)
-   {
      if (k > kmax)
      {
-       if (DEBUGLEVEL>3) {fprintferr(" K%ld",k);flusherr();}
+       if (DEBUGLEVEL>3) fprintferr("K%ld ",k);
+       ZincrementalGS(x, L, B, k, fl, gram);
        kmax = k;
-       for (j=1; j<=k; j++)
+     }
-         if (j==k || fl[j])
+     if (k != MARKED)
-         {
-           long av1 = avma;
-           u = gcoeff(x,k,j);
-           for (i=1; i<j; i++) if (fl[i])
-             u = divii(subii(mulii((GEN)B[i+1],u),
-                             mulii(gcoeff(L,k,i),gcoeff(L,j,i))),
-                       (GEN)B[i]);
-           u = gerepileuptoint(av1, u);
-           if (j<k) coeff(L,k,j)=(long)u;
-           else
-           {
-             s = signe(u);
-             if (s < 0) err(lllger3);
-             if (s)
-               { B[k+1]=(long)u; coeff(L,k,k)=un; fl[k]=1; }
-             else
-               { B[k+1]=B[k]; fl[k]=0; }
-           }
-         }
-     } else if (DEBUGLEVEL>5) {fprintferr(" %ld",k); flusherr();}
-     REDI(x,h,L,(GEN)B[k],kmax,k,k-1);
-     if (do_SWAPI(x,h,L,B,kmax,k,alpha,fl))
      {
-       if (k>2) k--;
+       if (!gram) ZRED(k,k-1, x,h,L,(GEN)B[k],kmax);
+       else  ZRED_gram(k,k-1, x,h,L,(GEN)B[k],kmax);
      }
+     if (do_ZSWAP(x,h,L,B,kmax,k,D,fl,gram))
+     {
+       if      (MARKED == k)   MARKED = k-1;
+       else if (MARKED == k-1) MARKED = k;
+       if (k > 2) k--;
+     }
      else
      {
-       for (l=k-2; l; l--)
+       if (k != MARKED)
-       {
+         for (l=k-2; l; l--)
-         REDI(x,h,L,(GEN)B[l+1],kmax,k,l);
-         if (low_stack(lim, stack_lim(av,1)))
          {
-           if(DEBUGMEM>1) err(warnmem,"lllgramall[1]");
+           if (!gram) ZRED(k,l, x,h,L,(GEN)B[l+1],kmax);
-           gptr[0]=&B; gptr[1]=&L; gptr[2]=&h; gptr[3]=&x;
+           else  ZRED_gram(k,l, x,h,L,(GEN)B[l+1],kmax);
-           gerepilemany(av,gptr,4);
+           if (low_stack(lim, stack_lim(av,1)))
+           {
+             if(DEBUGMEM>1) err(warnmem,"lllint[1], kmax = %ld", kmax);
+             gerepileall(av,h?4:3,&B,&L,&x,&h);
+           }
          }
-       }
        if (++k > n) break;
      }
      if (low_stack(lim, stack_lim(av,1)))
      {
-       if(DEBUGMEM>1) err(warnmem,"lllgramall[2]");
+       if(DEBUGMEM>1) err(warnmem,"lllint[2], kmax = %ld", kmax);
-       gptr[0]=&B; gptr[1]=&L; gptr[2]=&h; gptr[3]=&x;
+       gerepileall(av,h?4:3,&B,&L,&x,&h);
-       gerepilemany(av,gptr,4);
      }
    }
    if (DEBUGLEVEL>3) fprintferr("\n");
-   tetpil=avma; return gerepile(av0,tetpil,lllgramall_finish(h,fl,flag,n));
+   if (ptB)  *ptB  = B;
+   if (ptfl) *ptfl = fl;
+   if (pth)  *pth = h;
+   if (MARKED && MARKED != 1)
+   {
+     if (h) { swap(h[1], h[MARKED]); } else swap(x[1], x[MARKED]);
+   }
+   return h? h: x;
  }
- static GEN
+ /* Beware: this function can return NULL (dim x <= 1) */
- lllgramall0(GEN x, long flag) { return lllgramall(x,100,flag); }
+ GEN
+ lllint_i(GEN x, long D, int gram, GEN *pth, GEN *ptfl, GEN *ptB)
+ {
+   return lllint_marked(0, x,D,gram,pth,ptfl,ptB);
+ }
+ /* return x * lllint(x). No garbage collection */
+ GEN
+ lllint_ip(GEN x, long D)
+ {
+   GEN h = lllint_i(x, D, 0, NULL, NULL, NULL);
+   if (!h) h = x;
+   return h;
+ }
+ GEN
+ lllall(GEN x, long D, int gram, long flag)
+ {
+   gpmem_t av = avma;
+   GEN fl, junk, h = lllint_i(x, D, gram, &junk, &fl, NULL);
+   if (!h) return lll_trivial(x,flag);
+   return gerepilecopy(av, lll_finish(h,fl,flag));
+ }
+ GEN
+ lllint(GEN x) { return lllall(x,LLLDFT,0, lll_IM); }
+ GEN
+ lllkerim(GEN x) { return lllall(x,LLLDFT,0, lll_ALL); }
+ GEN
+ lllgramint(GEN x) { return lllall(x,LLLDFT,1, lll_IM | lll_GRAM); }
+ GEN
+ lllgramkerim(GEN x) { return lllall(x,LLLDFT,1, lll_ALL | lll_GRAM); }
  static int
  pslg(GEN x)
  {
    long tx;
    if (gcmp0(x)) return 2;
-   tx=typ(x); return is_scalar_t(tx)? 3: lgef(x);
+   tx = typ(x); return is_scalar_t(tx)? 3: lgef(x);
  }
  static GEN
- lllgramallgen(GEN x, long flag)
+ to_MP(GEN x, long prec)
  {
-   long av0=avma,av,tetpil,lx=lg(x),tu,i,j,k,l,n,lim;
+   GEN y = cgetr(prec); gaffect(x, y); return y;
-   GEN u,B,lam,q,cq,h,la,bb,p1,p2,p3,p4,fl;
+ }
-   int ps1,ps2,flc;
+ static GEN
+ col_to_MP(GEN x, long prec)
+ {
+   long j, l = lg(x);
+   GEN y = cgetg(l, t_COL);
+   for (j=1; j<l; j++) y[j] = (long)to_MP((GEN)x[j], prec);
+   return y;
+ }
-   if (typ(x) != t_MAT) err(typeer,"lllgramallgen");
+ static GEN
-   n=lx-1; if (n<=1) return lllall_trivial(x,n,flag);
+ to_mp(GEN x, long prec)
-   if (lg((GEN)x[1])!=lx) err(mattype1,"lllgramallgen");
+ {
+   GEN y;
+   if (typ(x) == t_INT) return x;
+   y = cgetr(prec); gaffect(x, y); return y;
+ }
+ static GEN
+ col_to_mp(GEN x, long prec)
+ {
+   long j, l = lg(x);
+   GEN y = cgetg(l, t_COL);
+   for (j=1; j<l; j++) y[j] = (long)to_mp((GEN)x[j], prec);
+   return y;
+ }
+ static GEN
+ mat_to_mp(GEN x, long prec)
+ {
+   long j, l = lg(x);
+   GEN y = cgetg(l, t_MAT);
+   for (j=1; j<l; j++) y[j] = (long)col_to_mp((GEN)x[j], prec);
+   return y;
+ }
-   fl = new_chunk(lx);
+ static int
+ REDgen(long k, long l, GEN h, GEN L, GEN B)
+ {
+   GEN u, q, cq;
+   long i;
-   av=avma; lim=stack_lim(av,1);
+   u = gcoeff(L,k,l);
-   B=cgetg(lx+1,t_COL);
+   if (pslg(u) < pslg((GEN)B[l+1])) return 0;
-   B[1]=un; lam=cgetg(lx,t_MAT);
-   for (j=1; j<lx; j++) lam[j]=lgetg(lx,t_COL);
+   q = gdeuc(u,(GEN)B[l+1]);
-   for (i=1; i<lx; i++)
+   cq = ginv(content(q));
-     for (j=1; j<=i; j++)
+   q = gneg(gmul(q,cq));
+   h[k] = ladd(gmul(cq,(GEN)h[k]), gmul(q,(GEN)h[l]));
+   coeff(L,k,l) = ladd(gmul(cq,gcoeff(L,k,l)), gmul(q,(GEN)B[l+1]));
+   for (i=1; i<l; i++)
+     coeff(L,k,i) = ladd(gmul(cq,gcoeff(L,k,i)), gmul(q,gcoeff(L,l,i)));
+   return 1;
+ }
+ static int
+ do_SWAPgen(GEN h, GEN L, GEN B, long k, GEN fl, int *flc)
+ {
+   GEN p1, la, la2, Bk;
+   long ps1, ps2, i, j, lx;
+   if (!fl[k-1]) return 0;
+   la = gcoeff(L,k,k-1); la2 = gsqr(la);
+   Bk = (GEN)B[k];
+   if (fl[k])
+   {
+     GEN q = gadd(la2, gmul((GEN)B[k-1],(GEN)B[k+1]));
+     ps1 = pslg(gsqr(Bk));
+     ps2 = pslg(q);
+     if (ps1 <= ps2 && (ps1 < ps2 || !*flc)) return 0;
+     *flc = (ps1 != ps2);
+     B[k] = ldiv(q, Bk);
+   }
+   swap(h[k-1], h[k]); lx = lg(L);
+   for (j=1; j<=k-2; j++) swap(coeff(L,k-1,j), coeff(L,k,j));
+   if (fl[k])
+   {
+     for (i=k+1; i<lx; i++)
      {
-       if (j<i && !fl[j]) coeff(lam,i,j)=coeff(lam,j,i)=zero;
+       GEN t = gcoeff(L,i,k);
-       else
+       p1 = gsub(gmul((GEN)B[k+1],gcoeff(L,i,k-1)), gmul(la,t));
-       {
+       coeff(L,i,k)   = ldiv(p1, Bk);
-         u=gcoeff(x,i,j); tu=typ(u);
+       p1 = gadd(gmul(la,gcoeff(L,i,k-1)), gmul((GEN)B[k-1],t));
-         if (! is_scalar_t(tu) && tu != t_POL) err(lllger4);
+       coeff(L,i,k-1) = ldiv(p1, Bk);
-         for (k=1; k<j; k++)
-           if (fl[k])
-             u=gdiv(gsub( gmul((GEN)B[k+1],u),
-                          gmul(gcoeff(lam,i,k),gcoeff(lam,j,k)) ),
-                    (GEN)B[k]);
-         if (j<i) { coeff(lam,i,j)=(long)u; coeff(lam,j,i)=zero; }
-         else
-         {
-           if (!gcmp0(u)) { B[i+1]=(long)u; coeff(lam,i,i)=un; fl[i]=1; }
-           else { B[i+1]=B[i]; coeff(lam,i,i)=zero; fl[i]=0; }
-         }
-       }
      }
-   k=2; h=idmat(n); flc=0;
+   }
-   for(;;)
+   else if (!gcmp0(la))
    {
-     u=gcoeff(lam,k,k-1);
+     p1 = gdiv(la2, Bk);
-     if (pslg(u) >= pslg((GEN)B[k]))
+     B[k+1] = B[k] = (long)p1;
+     for (i=k+2; i<=lx; i++) B[i] = ldiv(gmul(p1,(GEN)B[i]),Bk);
+     for (i=k+1; i<lx; i++)
+       coeff(L,i,k-1) = ldiv(gmul(la,gcoeff(L,i,k-1)), Bk);
+     for (j=k+1; j<lx-1; j++)
+       for (i=j+1; i<lx; i++)
+         coeff(L,i,j) = ldiv(gmul(p1,gcoeff(L,i,j)), Bk);
+   }
+   else
+   {
+     coeff(L,k,k-1) = zero;
+     for (i=k+1; i<lx; i++)
      {
-       q=gdeuc(u,(GEN)B[k]); cq=gdivsg(1,content(q)); q=gmul(q,cq); flc=1;
+       coeff(L,i,k) = coeff(L,i,k-1);
-       h[k]=lsub(gmul(cq,(GEN)h[k]),gmul(q,(GEN)h[k-1]));
+       coeff(L,i,k-1) = zero;
-       coeff(lam,k,k-1)=lsub(gmul(cq,gcoeff(lam,k,k-1)),gmul(q,(GEN)B[k]));
-       for (i=1; i<k-1; i++)
-         coeff(lam,k,i)=lsub(gmul(cq,gcoeff(lam,k,i)),gmul(q,gcoeff(lam,k-1,i)));
      }
-     ps1 = pslg(gsqr((GEN)B[k]));
+     B[k] = B[k-1]; fl[k] = 1; fl[k-1] = 0;
-     p3 = gmul((GEN)B[k-1],(GEN)B[k+1]);
+   }
-     la=gcoeff(lam,k,k-1); p4 = gmul(la,gcoeff(lam,k,k-1));
+   return 1;
-     ps2=pslg(gadd(p3,p4));
+ }
-     if (fl[k-1] && (ps1>ps2 || (ps1==ps2 && flc) || !fl[k]))
+ static void
+ incrementalGSgen(GEN x, GEN L, GEN B, long k, GEN fl)
+ {
+   GEN u = NULL; /* gcc -Wall */
+   long i, j, tu;
+   for (j=1; j<=k; j++)
+     if (j==k || fl[j])
      {
-       flc = (ps1!=ps2);
+       u = gcoeff(x,k,j); tu = typ(u);
-       swap(h[k-1],h[k]);
+       if (! is_scalar_t(tu) && tu != t_POL) err(lllger4);
-       for (j=1; j<=k-2; j++) swap(coeff(lam,k-1,j), coeff(lam,k,j));
+       for (i=1; i<j; i++)
-       if (fl[k])
+         if (fl[i])
-       {
+         {
-         for (i=k+1; i<=n; i++)
+           u = gsub(gmul((GEN)B[i+1],u), gmul(gcoeff(L,k,i),gcoeff(L,j,i)));
-         {
+           u = gdiv(u, (GEN)B[i]);
-           bb=gcoeff(lam,i,k);
+         }
-           coeff(lam,i,k)=ldiv(gsub(gmul((GEN)B[k+1],gcoeff(lam,i,k-1)),gmul(la,bb)),(GEN)B[k]);
+       coeff(L,k,j) = (long)u;
-           coeff(lam,i,k-1)=ldiv(gadd(gmul(la,gcoeff(lam,i,k-1)),gmul((GEN)B[k-1],bb)),(GEN)B[k]);
-          }
-         B[k]=ldiv(gadd(p3,p4),(GEN)B[k]);
-       }
-       else
-       {
-         if (!gcmp0(la))
-         {
-           p2=(GEN)B[k]; p1=gdiv(p4,p2);
-           for (i=k+1; i<lx; i++)
-             coeff(lam,i,k-1)=ldiv(gmul(la,gcoeff(lam,i,k-1)),p2);
-           for (j=k+1; j<lx-1; j++)
-             for (i=j+1; i<lx; i++)
-               coeff(lam,i,j)=ldiv(gmul(p1,gcoeff(lam,i,j)),p2);
-           B[k+1]=B[k]=(long)p1;
-           for (i=k+2; i<=lx; i++)
-             B[i]=ldiv(gmul(p1,(GEN)B[i]),p2);
-          }
-         else
-         {
-           coeff(lam,k,k-1)=zero;
-           for (i=k+1; i<lx; i++)
-           { coeff(lam,i,k)=coeff(lam,i,k-1); coeff(lam,i,k-1)=zero; }
-           B[k]=B[k-1]; fl[k]=1; fl[k-1]=0;
-          }
-       }
-       if (k>2) k--;
      }
+   if (gcmp0(u)) B[k+1] = B[k];
+   else
+   {
+     B[k+1] = coeff(L,k,k); coeff(L,k,k) = un; fl[k] = 1;
+   }
+ }
+ static GEN
+ lllgramallgen(GEN x, long flag)
+ {
+   long lx = lg(x), i, j, k, l, n;
+   gpmem_t av0 = avma, av, lim;
+   GEN B, L, h, fl;
+   int flc;
+   if (typ(x) != t_MAT) err(typeer,"lllgramallgen");
+   n = lx-1; if (n<=1) return lll_trivial(x,flag);
+   if (lg(x[1]) != lx) err(mattype1,"lllgramallgen");
+   fl = cgetg(lx, t_VECSMALL);
+   av = avma; lim = stack_lim(av,1);
+   B = gscalcol(gun, lx);
+   L = cgetg(lx,t_MAT);
+   for (j=1; j<lx; j++) { L[j] = (long)zerocol(n); fl[j] = 0; }
+   h = idmat(n);
+   for (i=1; i<lx; i++)
+     incrementalGSgen(x, L, B, i, fl);
+   flc = 0;
+   for(k=2;;)
+   {
+     if (REDgen(k, k-1, h, L, B)) flc = 1;
+     if (do_SWAPgen(h, L, B, k, fl, &flc)) { if (k > 2) k--; }
      else
      {
        for (l=k-2; l>=1; l--)
-       {
+         if (REDgen(k, l, h, L, B)) flc = 1;
-         u=gcoeff(lam,k,l);
-         if (pslg(u)>=pslg((GEN)B[l+1]))
-         {
-           q=gdeuc(u,(GEN)B[l+1]); cq=gdivsg(1,content(q));
-           q=gmul(q,cq); flc=1;
-           h[k]=lsub(gmul(cq,(GEN)h[k]),gmul(q,(GEN)h[l]));
-           coeff(lam,k,l)=lsub(gmul(cq,gcoeff(lam,k,l)),gmul(q,(GEN)B[l+1]));
-           for (i=1; i<l; i++)
-             coeff(lam,k,i)=lsub(gmul(cq,gcoeff(lam,k,i)),gmul(q,gcoeff(lam,l,i)));
-         }
-       }
        if (++k > n) break;
      }
      if (low_stack(lim, stack_lim(av,1)))
      {
-       GEN *gptr[4];
        if(DEBUGMEM>1) err(warnmem,"lllgramallgen");
-       gptr[0]=&B; gptr[1]=&lam; gptr[2]=&h;
+       gerepileall(av,3,&B,&L,&h);
-       gerepilemany(av,gptr,3);
      }
    }
-   tetpil=avma;
+   return gerepilecopy(av0, lll_finish(h,fl,flag));
-   return gerepile(av0,tetpil,lllgramall_finish(h,fl,flag,n));
  }
- /* compute B[k], update mu(k,1..k-1) */
+ static GEN
+ _mul2n(GEN x, long e) { return e? gmul2n(x, e): x; }
+ /* E = scaling exponents
+  * F = backward scaling exponents
+  * X = lattice basis, Xs = associated scaled basis
+  * Q = vector of Householder transformations
+  * h = base change matrix (from X0)
+  * R = from QR factorization of Xs[1..k-1] */
  static int
- get_Gram_Schmidt(GEN x, GEN mu, GEN B, long k)
+ HRS(int MARKED, long k, int prim, long kmax, GEN X, GEN Xs, GEN h, GEN R,
+     GEN Q, GEN E, GEN F)
  {
-   GEN s,A = cgetg(k+1, t_COL); /* scratch space */
+   long e, i, N = lg(X[k]);
-   long av,i,j;
+   const long prec = MEDDEFAULTPREC; /* 128 bits */
-   A[1] = coeff(x,k,1);
+   GEN q, tau2, rk;
-   for(j=1;j<k;)
+   int rounds = 0, overf;
+   E[k] = prim? E[k-1]: 0;
+   F[k] = 0;
+   Xs[k] = E[k]? lmul2n((GEN)X[k], E[k]): (long)dummycopy((GEN)X[k]);
+   rounds = 0;
+   if (k == MARKED) goto DONE; /* no size reduction/scaling */
+ UP:
+   rk = ApplyAllQ(Q, col_to_MP((GEN)Xs[k], prec), k);
+   overf = 0;
+   for (i = k-1; i > 0; i--)
+   { /* size seduction of Xs[k] */
+     GEN q2, mu, muint;
+     long ex;
+     gpmem_t av = avma;
+     mu = mpdiv((GEN)rk[i], gcoeff(R,i,i));
+     if (!signe(mu)) { avma = av; continue; }
+     mu = _mul2n(mu, -F[i]); /*backward scaling*/
+     ex = gexpo(mu);
+     if (ex > 10) {
+       overf = 1; /* large size reduction */
+       muint = ex > 30? ceil_safe(mu): ground(mu);
+     }
+     else if (fabs(gtodouble(mu)) > 0.51)
+       muint = ground(mu);
+     else { avma = av; continue; }
+     q = gerepileuptoint(av, negi(muint));
+     q2= _mul2n(q, F[i]); /*backward scaling*/
+     Zupdate_col(k, i, q2, N-1, Xs);
+     rk = gadd(rk, gmul(q2, (GEN)R[i]));
+     /* if in HRS', propagate the last SR on X (for SWAP test) */
+     if (prim && (i == k-1)) {
+       Zupdate_col(k, i, q, kmax, h);
+        update_col(k, i, q, X);
+     }
+   }
+   /* If large size-reduction performed, restart from exact data */
+   if (overf)
    {
-     coeff(mu,k,j) = ldiv((GEN)A[j],(GEN)B[j]);
+     if (++rounds > 100) return 0; /* detect infinite loop */
-     j++; av = avma;
+     goto UP;
-     /* A[j] <-- x[k,j] - sum_{i<j} mu[j,i] A[i] */
-     s = gmul(gcoeff(mu,j,1),(GEN)A[1]);
-     for (i=2; i<j; i++) s = gadd(s, gmul(gcoeff(mu,j,i),(GEN)A[i]));
-     s = gneg(s); A[j] = lpileupto(av, gadd(gcoeff(x,k,j),s));
    }
-   B[k] = A[k]; return (gsigne((GEN)B[k]) > 0);
+   if (prim || k == 1) goto DONE; /* DONE */
+   rounds = 0;
+   /* rescale Xs[k] ? */
+   tau2 = signe(rk[k])? gsqr((GEN)rk[k]): gzero;
+   for (i = k+1; i < N; i++)
+     if (signe(rk[i])) tau2 = mpadd(tau2, gsqr((GEN)rk[i]));
+   q = mpdiv(gsqr(gcoeff(R,1,1)), tau2);
+   e = gexpo(q)/2 + F[1]; /* ~ log_2 (||Xs[1]|| / tau), backward scaling*/
+   if (e > 0)
+   { /* rescale */
+     if (e > 30) e = 30;
+     Xs[k] = lmul2n((GEN)Xs[k], e);
+     E[k] += e; goto UP; /* one more round */
+   }
+   else if (e < 0)
+   { /* enable backward scaling */
+     for (i = 1; i < k; i++) F[i] -= e;
+   }
+ DONE:
+   rk = ApplyAllQ(Q, col_to_MP((GEN)Xs[k], prec), k);
+   (void)FindApplyQ(rk, R, NULL, k, Q, prec);
+   return 1;
  }
- /* x = Gram(b_i). If precision problems return NULL if flag=1, error otherwise.
+ static GEN
-  * Quality ratio = Q = (alpha-1)/alpha. Suggested value: alpha = 100
+ rescale_to_int(GEN x)
-  */
+ {
+   long e, prec = gprecision(x);
+   GEN y;
+   if (!prec) return Q_primpart(x);
+   y = gmul2n(x, bit_accuracy(prec) - gexpo(x));
+   return gcvtoi(y, &e);
+ }
  GEN
- lllgramintern(GEN x, long alpha, long flag, long prec)
+ lll_scaled(int MARKED, GEN X0, long D)
  {
-   GEN xinit,L,h,B,L1,QR;
+   GEN delta, X, Xs, h, R, Q, E, F;
-   long retry = 2, av = avma,lim,l,i,j,k,k1,lx=lg(x),n,kmax,KMAX;
+   long j, kmax = 1, k, N = lg(X0);
-   long last_prec;
+   gpmem_t lim, av, av0 = avma;
+   int retry = 0;
-   if (typ(x) != t_MAT) err(typeer,"lllgram");
+   delta = stor(D-1, DEFAULTPREC);
-   n=lx-1; if (n && lg((GEN)x[1])!=lx) err(mattype1,"lllgram");
+   delta = divrs(delta,D);
-   if (n<=1) return idmat(n);
+   E  = vecsmall_const(N-1, 0);
-   xinit = x; lim = stack_lim(av,1);
+   F  = vecsmall_const(N-1, 0);
-   for (k=2,j=1; j<lx; j++)
+   av = avma; lim = stack_lim(av, 1);
+   h = idmat(N-1);
+   X = rescale_to_int(X0);
+ PRECPB:
+   k = 1;
+   if (retry++) return _vec(h);
+   Q  = zerovec(N-1);
+   Xs = zeromat(N-1, N-1);
+   R  = cgetg(N, t_MAT);
+   for (j=1; j<N; j++) R[j] = (long)zerocol(N-1);
+   while (k < N)
    {
-     GEN p1=(GEN)x[j];
+     gpmem_t av1;
-     for (i=1; i<=j; i++)
+     if (k == 1) {
-       if (typ(p1[i]) == t_REAL) { l = lg((GEN)p1[i]); if (l>k) k=l; }
+       HRS(MARKED, 1, 0, kmax, X, Xs, h, R, Q, E, F);
+       k = 2;
+     }
+     if (k > kmax) {
+       kmax = k;
+       if (DEBUGLEVEL>3) {fprintferr(" K%ld",k);flusherr();}
+     }
+     if (!HRS(MARKED, k, 1, kmax, X, Xs, h, R, Q, E, F)) goto PRECPB;
+     av1 = avma;
+     if (mpcmp( mpmul(delta, gsqr(gcoeff(R, k-1,k-1))),
+                mpadd(gsqr(gcoeff(R,k-1,k)), gsqr(gcoeff(R, k,k)))) > 0)
+     {
+       if (DEBUGLEVEL>3 && k == kmax)
+       {
+         GEN q = mpsub( mpmul(delta, gsqr(gcoeff(R, k-1,k-1))),
+                        gsqr(gcoeff(R,k-1,k)));
+         fprintferr(" (%ld)", gexpo(q) - gexpo(gsqr(gcoeff(R, k,k))));
+       }
+       swap(X[k], X[k-1]);
+       swap(h[k], h[k-1]);
+       if      (MARKED == k)   MARKED = k-1;
+       else if (MARKED == k-1) MARKED = k;
+       avma = av1; k--;
+     }
+     else
+     {
+       avma = av1;
+       if (!HRS(MARKED, k, 0, kmax, X, Xs, h, R, Q, E, F)) goto PRECPB;
+       k++;
+     }
+     if (low_stack(lim, stack_lim(av,1)))
+     {
+       if(DEBUGMEM>1) err(warnmem,"lllfp[1]");
+       gerepileall(av,5,&X,&Xs, &R,&h,&Q);
+     }
    }
+   return gerepilecopy(av0, h);
+ }
+ /* x = Gram(b_i). If precision problems return NULL if flag=1, error otherwise.
+  * Quality ratio = delta = (D-1)/D. Suggested value: D = 100
+  *
+  * If MARKED != 0 make sure e[MARKED] is the first vector of the output basis
+  * (which may then not be LLL-reduced) */
+ GEN
+ lllfp_marked(int MARKED, GEN x, long D, long flag, long prec, int gram)
+ {
+   GEN xinit,L,h,B,L1,delta, Q;
+   long retry = 2, lx = lg(x), hx, l, i, j, k, k1, n, kmax, KMAX;
+   gpmem_t av0 = avma, av, lim;
+   if (typ(x) != t_MAT) err(typeer,"lllfp");
+   n = lx-1; if (n <= 1) return idmat(n);
+ #if 0 /* doesn't work yet */
+   return lll_scaled(MARKED, x, D);
+ #endif
+   hx = lg(x[1]);
+   if (gram && hx != lx) err(mattype1,"lllfp");
+   delta = stor(D-1, DEFAULTPREC);
+   delta = divrs(delta,D);
+   xinit = x;
+   av = avma; lim = stack_lim(av,1);
+   if (gram) {
+     for (k=2,j=1; j<lx; j++)
+     {
+       GEN p1=(GEN)x[j];
+       for (i=1; i<=j; i++)
+         if (typ(p1[i]) == t_REAL) { l = lg(p1[i]); if (l>k) k=l; }
+     }
+   } else {
+     for (k=2,j=1; j<lx; j++)
+     {
+       GEN p1=(GEN)x[j];
+       for (i=1; i<hx; i++)
+         if (typ(p1[i]) == t_REAL) { l = lg(p1[i]); if (l>k) k=l; }
+     }
+   }
    if (k == 2)
    {
-     if (!prec) return lllgramint(x);
+     if (!prec) return gram? lllgramint(x): lllint(x);
      x = gmul(x, realun(prec+1));
    }
    else
    {
      if (prec < k) prec = k;
-     x = gprec_w(x, prec+1);
+     x = mat_to_mp(x, prec+1);
    }
   /* kmax = maximum column index attained during this run
    * KMAX = same over all runs (after PRECPB) */
    kmax = KMAX = 1;
+   Q = zerovec(n);
  PRECPB:
    switch(retry--)
    {
      case 2: h = idmat(n); break; /* entry */
      case 1:
-       if (kmax > 2)
+       if (DEBUGLEVEL > 3) fprintferr("\n");
+       if (flag == 2) return _vec(h);
+       if (gram && kmax > 2)
        { /* some progress but precision loss, try again */
          prec = (prec<<1)-2; kmax = 1;
-         if (DEBUGLEVEL > 3) fprintferr("\n");
+         if (DEBUGLEVEL) err(warnprec,"lllfp",prec);
-         if (DEBUGLEVEL) err(warnprec,"lllgramintern",prec);
+         x = gprec_w(xinit,prec);
-         x = qf_base_change(gprec_w(xinit,prec),h,1);
+         if (gram) x = qf_base_change(x,h,1); else x = gmul(x,h);
-         {
+         gerepileall(av,2,&h,&x);
-           GEN *gsav[2]; gsav[0]=&h; gsav[1]=&x;
+         if (DEBUGLEVEL) err(warner,"lllfp starting over");
-           gerepilemany(av, gsav, 2);
-         }
-         if (DEBUGLEVEL) err(warner,"lllgramintern starting over");
          break;
        } /* fall through */
      case 0: /* give up */
        if (DEBUGLEVEL > 3) fprintferr("\n");
-       if (DEBUGLEVEL) err(warner,"lllgramintern giving up");
+       if (DEBUGLEVEL) err(warner,"lllfp giving up");
        if (flag) { avma=av; return NULL; }
-       if (DEBUGLEVEL) outerr(xinit);
        err(lllger3);
    }
-   last_prec = prec;
-   QR = cgetr(prec+1); affsr(alpha-1,QR);
-   QR = divrs(QR,alpha);
    L=cgetg(lx,t_MAT);
    B=cgetg(lx,t_COL);
    for (j=1; j<lx; j++) { L[j] = (long)zerocol(n); B[j] = zero; }
-   k=2; B[1]=coeff(x,1,1);
+   if (gram && !incrementalGS(x, L, B, 1))
-   if (gcmp0((GEN)B[1]))
    {
-     if (flag) return NULL;
+     if (!flag) return NULL;
      err(lllger3);
    }
    if (DEBUGLEVEL>5) fprintferr("k =");
-   for(;;)
+   for(k=2;;)
    {
-     if (k>kmax)
+     if (k > kmax)
      {
        kmax = k; if (KMAX < kmax) KMAX = kmax;
        if (DEBUGLEVEL>3) {fprintferr(" K%ld",k);flusherr();}
-       if (!get_Gram_Schmidt(x,L,B,k)) goto PRECPB;
+       if (gram) j = incrementalGS(x, L, B, k);
+       else      j = Householder_get_mu(x, L, B, k, Q, prec);
+       if (!j) goto PRECPB;
      }
      else if (DEBUGLEVEL>5) fprintferr(" %ld",k);
      L1 = gcoeff(L,k,k-1);
-     if (typ(L1) == t_REAL &&
+     if (typ(L1) == t_REAL && expo(L1) + 20 > bit_accuracy(lg(L1)))
-         (2*gexpo(L1) > bit_accuracy(lg(L1)) || 2*lg(L1) < last_prec))
      {
-       last_prec = lg(L1);
        if (DEBUGLEVEL>3)
-         fprintferr("\nRecomputing Gram-Schmidt, kmax = %ld, prec was %ld\n",
+         fprintferr("\nRecomputing Gram-Schmidt, kmax = %ld\n", kmax);
-                    kmax,last_prec);
+       if (!gram) goto PRECPB;
        for (k1=1; k1<=kmax; k1++)
-         if (!get_Gram_Schmidt(x,L,B,k1)) goto PRECPB;
+         if (!incrementalGS(x, L, B, k1)) goto PRECPB;
      }
-     RED(x,h,L,KMAX,k,k-1);
+     if (k != MARKED)
-     if (do_SWAP(x,h,L,B,kmax,k,QR))
      {
+       if (!gram) j = RED(k,k-1, x,h,L,KMAX);
+       else  j = RED_gram(k,k-1, x,h,L,KMAX);
+       if (!j) goto PRECPB;
+     }
+     if (do_SWAP(x,h,L,B,kmax,k,delta,gram))
+     {
+       if      (MARKED == k)   MARKED = k-1;
+       else if (MARKED == k-1) MARKED = k;
        if (!B[k]) goto PRECPB;
+       Q[k] = Q[k-1] = zero;
        if (k>2) k--;
      }
      else
      {
-       for (l=k-2; l; l--) RED(x,h,L,KMAX,k,l);
+       if (k != MARKED)
-       if (++k > n) break;
+         for (l=k-2; l; l--)
+         {
+           if (!gram) j = RED(k,l, x,h,L,KMAX);
+           else  j = RED_gram(k,l, x,h,L,KMAX);
+           if (!j) goto PRECPB;
+           if (low_stack(lim, stack_lim(av,1)))
+           {
+             if(DEBUGMEM>1) err(warnmem,"lllfp[1], kmax = %ld", kmax);
+             gerepileall(av,5,&B,&L,&h,&x,&Q);
+           }
+         }
+       if (++k > n)
+       {
+         if (!gram && Q[n-1] == zero)
+         {
+           if (DEBUGLEVEL>3) fprintferr("\nChecking LLL basis\n");
+           j = Householder_get_mu(gmul(xinit,h), L, B, n, Q, prec);
+           if (!j) goto PRECPB;
+           k = 2; continue;
+         }
+         break;
+       }
      }
      if (low_stack(lim, stack_lim(av,1)))
      {
-       GEN *gptr[5]; gptr[0]=&B; gptr[1]=&L; gptr[2]=&h; gptr[3]=&x; gptr[4]=&QR;
+       if(DEBUGMEM>1) err(warnmem,"lllfp[2], kmax = %ld", kmax);
-       if(DEBUGMEM>1)
+       gerepileall(av,5,&B,&L,&h,&x,&Q);
-       {
-         if (DEBUGLEVEL > 3) fprintferr("\n");
-         err(warnmem,"lllgram");
-       }
-       gerepilemany(av,gptr,5);
      }
    }
    if (DEBUGLEVEL>3) fprintferr("\n");
-   return gerepilecopy(av, h);
+   if (MARKED && MARKED != 1) swap(h[1], h[MARKED]);
+   return gerepilecopy(av0, h);
  }
- static GEN
+ GEN
- lllgram_noerr(GEN x,long prec) { return lllgramintern(x,100,1,prec); }
+ lllgramintern(GEN x, long D, long flag, long prec)
+ {
+   return lllfp_marked(0, x,D,flag,prec,1);
+ }
  GEN
- lllgram(GEN x,long prec) { return lllgramintern(x,100,0,prec); }
+ lllintern(GEN x, long D, long flag, long prec)
+ {
+   return lllfp_marked(0, x,D,flag,prec,0);
+ }
  GEN
+ lllgram(GEN x,long prec) { return lllgramintern(x,LLLDFT,0,prec); }
+ GEN
+ lll(GEN x,long prec) { return lllintern(x,LLLDFT,0,prec); }
+ GEN
  qflll0(GEN x, long flag, long prec)
  {
    switch(flag)
-Line 960  qflll0(GEN x, long flag, long prec)
+Line 1308  qflll0(GEN x, long flag, long prec)
 Line 960  qflll0(GEN x, long flag, long prec)
 Line 1308  qflll0(GEN x, long flag, long prec)
      case 0: return lll(x,prec);
      case 1: return lllint(x);
      case 2: return lllintpartial(x);
-     case 3: return lllrat(x);
      case 4: return lllkerim(x);
      case 5: return lllkerimgen(x);
-     case 7: return lll1(x,prec);
      case 8: return lllgen(x);
-     case 9: return lllintwithcontent(x);
      default: err(flagerr,"qflll");
    }
    return NULL; /* not reached */
-Line 980  qflllgram0(GEN x, long flag, long prec)
+Line 1325  qflllgram0(GEN x, long flag, long prec)
 Line 980  qflllgram0(GEN x, long flag, long prec)
 Line 1325  qflllgram0(GEN x, long flag, long prec)
      case 1: return lllgramint(x);
      case 4: return lllgramkerim(x);
      case 5: return lllgramkerimgen(x);
-     case 7: return lllgram1(x,prec);
      case 8: return lllgramgen(x);
      default: err(flagerr,"qflllgram");
    }
    return NULL; /* not reached */
  }
- /* x est la matrice d'une base b_i; retourne la matrice u (entiere
+ GEN
-  * unimodulaire) d'une base LLL-reduite c_i en fonction des b_i (la base
+ gram_matrix(GEN b)
-  * reduite est c=b*u).
-  */
- static GEN
- lll_proto(GEN x, GEN f(GEN, long), long prec)
  {
-   long lx=lg(x),i,j,av,av1;
+   long i,j, lx = lg(b);
    GEN g;
+   if (typ(b) != t_MAT) err(typeer,"gram");
-   if (typ(x) != t_MAT) err(typeer,"lll_proto");
+   g = cgetg(lx,t_MAT);
-   if (lx==1) return cgetg(1,t_MAT);
-   av=avma; g=cgetg(lx,t_MAT);
-   for (j=1; j<lx; j++) g[j]=lgetg(lx,t_COL);
    for (i=1; i<lx; i++)
+   {
+     g[i] = lgetg(lx,t_COL);
      for (j=1; j<=i; j++)
-       coeff(g,i,j)=coeff(g,j,i)=(long)gscal((GEN)x[i],(GEN)x[j]);
+       coeff(g,i,j) = coeff(g,j,i) = (long)gscal((GEN)b[i],(GEN)b[j]);
-   av1=avma; x = f(g,prec);
+   }
-   if (!x) { avma=av; return NULL; }
+   return g;
-   return gerepile(av,av1,x);
  }
  GEN
- lllintern(GEN x,long flag,long prec)
+ lllgen(GEN x) {
- {
+   gpmem_t av = avma;
-   return lll_proto(x,flag? lllgram_noerr: lllgram,prec);
+   return gerepileupto(av, lllgramgen(gram_matrix(x)));
  }
  GEN
- lll(GEN x,long prec) { return lll_proto(x,lllgram,prec); }
+ lllkerimgen(GEN x) {
+   gpmem_t av = avma;
+   return gerepileupto(av, lllgramallgen(gram_matrix(x),lll_ALL));
+ }
  GEN
- lll1(GEN x,long prec) { return lll_proto(x,lllgram1,prec); }
- GEN
- lllrat(GEN x) { return lll_proto(x,lllgram,lll_ALL); }
- GEN
- lllint(GEN x) { return lll_proto(x,lllgramall0,lll_IM); }
- GEN
- lllgen(GEN x) { return lll_proto(x,lllgramallgen,lll_IM); }
- GEN
  lllgramgen(GEN x) { return lllgramallgen(x,lll_IM); }
  GEN
  lllgramkerimgen(GEN x) { return lllgramallgen(x,lll_ALL); }
- static GEN
- lllkerim_proto(GEN x, GEN f(GEN,long))
- {
-   long lx=lg(x), i,j,av,av1;
-   GEN g;
-   if (typ(x) != t_MAT) err(typeer,"lllkerim_proto");
-   if (lx==1)
-   {
-     g=cgetg(3,t_VEC);
-     g[1]=lgetg(1,t_MAT);
-     g[2]=lgetg(1,t_MAT); return g;
-   }
-   if (lg((GEN)x[1])==1)
-   {
-     g=cgetg(3,t_VEC);
-     g[1]=(long)idmat(lx-1);
-     g[2]=lgetg(1,t_MAT); return g;
-   }
-   av=avma; g=cgetg(lx,t_MAT);
-   for (j=1; j<lx; j++) g[j]=lgetg(lx,t_COL);
-   for (i=1; i<lx; i++)
-     for (j=1; j<=i; j++)
-       coeff(g,i,j)=coeff(g,j,i)=(long)gscal((GEN)x[i],(GEN)x[j]);
-   av1=avma; return gerepile(av,av1,f(g,lll_ALL));
- }
- GEN
- lllkerim(GEN x) { return lllkerim_proto(x,lllgramall0); }
- GEN
- lllkerimgen(GEN x) { return lllkerim_proto(x,lllgramallgen); }
- /* x est ici la matrice de GRAM des bi.  */
- GEN
- lllgram1(GEN x, long prec)
- {
-   GEN mu,u,B,BB,BK,p,q,r,cst,unreel,sv,mu1,mu2;
-   long av,tetpil,lim,l,i,j,k,lx=lg(x),n,e;
-   if (typ(x) != t_MAT) err(typeer,"lllgram1");
-   if (lg((GEN)x[1])!=lx) err(mattype1,"lllgram1"); n=lx-1;
-   if (n<=1) return idmat(n);
-   cst=gdivgs(stoi(99),100); /* LLL proposent 0.75 */
-   if (prec)
-   {
-     unreel = realun(prec+1);
-     x = gmul(x,unreel);
-     cst = gmul(cst,unreel);
-   }
-   av=avma; lim=stack_lim(av,1);
-   mu=gtrans(sqred(x)); B=cgetg(lx,t_COL);
-   for (i=1,l=0; i<=n; i++)
-   {
-     if (gsigne((GEN)(B[i]=coeff(mu,i,i)))>0) l++;
-     coeff(mu,i,i)=un;
-   }
-   if (l<n) err(lllger3);
-   u=idmat(n); k=2;
-   do
-   {
-     if (!gcmp0(r=grndtoi(gcoeff(mu,k,k-1),&e)))
-     {
-       u[k]=lsub((GEN)u[k],gmul(r,(GEN)u[k-1]));
-       for (j=1; j<k-1; j++)
-         coeff(mu,k,j)=lsub(gcoeff(mu,k,j),gmul(r,gcoeff(mu,k-1,j)));
-       mu1=(GEN)(coeff(mu,k,k-1)=lsub(gcoeff(mu,k,k-1),r));
-     }
-     else mu1=gcoeff(mu,k,k-1);
-     q=gmul((GEN)B[k-1],gsub(cst,mu2=gsqr(mu1)));
-     if (gcmp(q,(GEN)B[k])>0)
-     {
-       BB=gadd((GEN)B[k],gmul((GEN)B[k-1],mu2));
-       coeff(mu,k,k-1)=ldiv(gmul(mu1,(GEN)B[k-1]),BB);
-       B[k]=lmul((GEN)B[k-1],BK=gdiv((GEN)B[k],BB));
-       B[k-1]=(long)BB;
-       swap(u[k-1],u[k]);
-       for (j=1; j<=k-2; j++) swap(coeff(mu,k-1,j), coeff(mu,k,j));
-       for (i=k+1; i<=n; i++)
-       {
-         p=gcoeff(mu,i,k);
-         coeff(mu,i,k)=lsub(gcoeff(mu,i,k-1),gmul(mu1,p));
-         coeff(mu,i,k-1)=ladd(gmul(BK,p),gmul(gcoeff(mu,k,k-1),gcoeff(mu,i,k-1)));
-       }
-       if (k>2) k--;
-     }
-     else
-     {
-       for (l=k-2; l; l--)
-       {
-         if (!gcmp0(r=grndtoi(gcoeff(mu,k,l),&e)))
-         {
-           u[k]=lsub((GEN)u[k],gmul(r,(GEN)u[l]));
-           for (j=1; j<l; j++)
-             coeff(mu,k,j)=lsub(gcoeff(mu,k,j),gmul(r,gcoeff(mu,l,j)));
-           coeff(mu,k,l)=lsub(gcoeff(mu,k,l),r);
-          }
-       }
-       k++;
-     }
-     if (low_stack(lim, stack_lim(av,1)))
-     {
-       if(DEBUGMEM>1) err(warnmem,"lllgram1");
-       tetpil=avma;
-       sv=cgetg(4,t_VEC);
-       sv[1]=lcopy(B); sv[2]=lcopy(u); sv[3]=lcopy(mu);
-       sv=gerepile(av,tetpil,sv);
-       B=(GEN)sv[1]; u=(GEN)sv[2]; mu=(GEN)sv[3];
-     }
-   }
-   while (k<=n);
-   return gerepilecopy(av,u);
- }
- GEN
- lllgramint(GEN x)
- {
-   return lllgramall0(x,lll_IM);
- }
- GEN
- lllgramkerim(GEN x)
- {
-   return lllgramall0(x,lll_ALL);
- }
  /*  This routine is functionally similar to lllint when all = 0.
   *
   *    The input is an integer matrix mat (not necessarily square) whose
-Line 1192  GEN
+Line 1394  GEN
 Line 1192  GEN
 Line 1394  GEN
  lllintpartialall(GEN m, long flag)
  {
    const long ncol = lg(m)-1;
-   const long ltop1 = avma;
+   const gpmem_t ltop1 = avma;
    long ncolnz;
-   GEN tm1, tm2, mid, *gptr[4];
+   GEN tm1, tm2, mid;
    if (typ(m) != t_MAT) err(typeer,"lllintpartialall");
    if (ncol <= 1) return idmat(ncol);
-Line 1217  lllintpartialall(GEN m, long flag)
+Line 1419  lllintpartialall(GEN m, long flag)
 Line 1217  lllintpartialall(GEN m, long flag)
 Line 1419  lllintpartialall(GEN m, long flag)
   */
      while (npass2 < 2 || progress)
      {
-       GEN dot12new,q = gdivround(dot12, dot22);
+       GEN dot12new,q = diviiround(dot12, dot22);
        npass2++; progress = signe(q);
        if (progress)
-Line 1244  lllintpartialall(GEN m, long flag)
+Line 1446  lllintpartialall(GEN m, long flag)
 Line 1244  lllintpartialall(GEN m, long flag)
 Line 1446  lllintpartialall(GEN m, long flag)
        /* Occasionally (including final pass) do garbage collection.  */
        if (npass2 % 8 == 0 || !progress)
-       {
+         gerepileall(ltop1, 4, &dot11,&dot12,&dot22,&tm);
-         gptr[0] = &dot11; gptr[1] = &dot12;
-         gptr[2] = &dot22; gptr[3] = &tm;
-         gerepilemany(ltop1, gptr, 4);
-       }
      } /* while npass2 < 2 || progress */
      {
-Line 1282  lllintpartialall(GEN m, long flag)
+Line 1480  lllintpartialall(GEN m, long flag)
 Line 1282  lllintpartialall(GEN m, long flag)
 Line 1480  lllintpartialall(GEN m, long flag)
          GEN q1neg = subii(mulii(dot12, dot2i), mulii(dot22, dot1i));
          GEN q2neg = subii(mulii(dot12, dot1i), mulii(dot11, dot2i));
-         q1neg = gdivround(q1neg, det12);
+         q1neg = diviiround(q1neg, det12);
-         q2neg = gdivround(q2neg, det12);
+         q2neg = diviiround(q2neg, det12);
          coeff(tm1, 1, icol) = ladd(gmul(q1neg, gcoeff(tm,1,1)),
                                     gmul(q2neg, gcoeff(tm,1,2)));
          coeff(tm1, 2, icol) = ladd(gmul(q1neg, gcoeff(tm,2,1)),
-Line 1298  lllintpartialall(GEN m, long flag)
+Line 1496  lllintpartialall(GEN m, long flag)
 Line 1298  lllintpartialall(GEN m, long flag)
 Line 1496  lllintpartialall(GEN m, long flag)
      fprintferr("tm1 = %Z",tm1);
      fprintferr("mid = %Z",mid);
    }
-   gptr[0] = &tm1; gptr[1] = &mid;
+   gerepileall(ltop1,2, &tm1, &mid);
-   gerepilemany(ltop1, gptr, 2);
    {
     /* For each pair of column vectors v and w in mid * tm2,
      * try to replace (v, w) by (v, v - q*w) for some q.
      * We compute all inner products and check them repeatedly.
      */
-     const long ltop3 = avma; /* Excludes region with tm1 and mid */
+     const gpmem_t ltop3 = avma; /* Excludes region with tm1 and mid */
-     long icol, lim, reductions, npass = 0;
+     const gpmem_t lim = stack_lim(ltop3,1);
+     long icol, reductions, npass = 0;
      GEN dotprd = cgetg(ncol+1, t_MAT);
      tm2 = idmat(ncol);
-Line 1319  lllintpartialall(GEN m, long flag)
+Line 1517  lllintpartialall(GEN m, long flag)
 Line 1319  lllintpartialall(GEN m, long flag)
 Line 1517  lllintpartialall(GEN m, long flag)
          coeff(dotprd,jcol,icol) = coeff(dotprd,icol,jcol) =
            (long)gscali((GEN)mid[icol], (GEN)mid[jcol]);
      } /* for icol */
-     lim = stack_lim(ltop3,1);
      for(;;)
      {
        reductions = 0;
-Line 1330  lllintpartialall(GEN m, long flag)
+Line 1527  lllintpartialall(GEN m, long flag)
 Line 1330  lllintpartialall(GEN m, long flag)
 Line 1527  lllintpartialall(GEN m, long flag)
          for (ijdif=1; ijdif < ncol; ijdif++)
          {
-           const long previous_avma = avma;
+           const gpmem_t previous_avma = avma;
            jcol = (icol + ijdif - 1) % ncol; jcol++; /* Hack for NeXTgcc 2.5.8 */
            k1 = (cmpii(gcoeff(dotprd,icol,icol),
-Line 1338  lllintpartialall(GEN m, long flag)
+Line 1535  lllintpartialall(GEN m, long flag)
 Line 1338  lllintpartialall(GEN m, long flag)
 Line 1535  lllintpartialall(GEN m, long flag)
                  ? icol : jcol;  /* index of larger column */
            k2 = icol + jcol - k1;        /* index of smaller column */
            codi = gcoeff(dotprd,k2,k2);
-           q = gcmp0(codi)? gzero: gdivround(gcoeff(dotprd,k1,k2), codi);
+           q = gcmp0(codi)? gzero: diviiround(gcoeff(dotprd,k1,k2), codi);
            /* Try to subtract a multiple of column k2 from column k1.  */
            if (gcmp0(q)) avma = previous_avma;
-Line 1360  lllintpartialall(GEN m, long flag)
+Line 1557  lllintpartialall(GEN m, long flag)
 Line 1360  lllintpartialall(GEN m, long flag)
 Line 1557  lllintpartialall(GEN m, long flag)
          if (low_stack(lim, stack_lim(ltop3,1)))
          {
            if(DEBUGMEM>1) err(warnmem,"lllintpartialall");
-           gptr[0] = &dotprd; gptr[1] = &tm2;
+           gerepileall(ltop3, 2, &dotprd,&tm2);
-           gerepilemany(ltop3, gptr, 2);
          }
        } /* for icol */
        if (!reductions) break;
-Line 1423  lllintpartial(GEN mat)
+Line 1619  lllintpartial(GEN mat)
 Line 1423  lllintpartial(GEN mat)
 Line 1619  lllintpartial(GEN mat)
  /********************************************************************/
  /**                                                                **/
- /**                   LINEAR & ALGEBRAIC DEPENDANCE                **/
+ /**                   LINEAR & ALGEBRAIC DEPENDENCE                **/
  /**                                                                **/
  /********************************************************************/
+ GEN pslq(GEN x, long prec);
+ GEN pslqL2(GEN x, long prec);
  GEN
  lindep0(GEN x,long bit,long prec)
  {
-   if (!bit) return lindep(x,prec);
+   switch (bit)
-   if (bit>0) return lindep2(x,bit);
+   {
-   return deplin(x);
+     case 0: return pslq(x,prec);
+     case -1: return lindep(x,prec);
+     case -2: return deplin(x);
+     case -3: return pslqL2(x,prec);
+     default: return lindep2(x,labs(bit));
+   }
  }
  static int
-Line 1446  real_indep(GEN re, GEN im, long bitprec)
+Line 1650  real_indep(GEN re, GEN im, long bitprec)
 Line 1446  real_indep(GEN re, GEN im, long bitprec)
 Line 1650  real_indep(GEN re, GEN im, long bitprec)
  GEN
  lindep2(GEN x, long bit)
  {
-   long tx=typ(x),lx=lg(x),ly,i,j,e, av = avma;
+   long tx=typ(x), lx=lg(x), ly, i, j, e;
+   gpmem_t av = avma;
    GEN re,im,p1,p2;
    if (! is_vec_t(tx)) err(typeer,"lindep2");
-Line 1467  lindep2(GEN x, long bit)
+Line 1672  lindep2(GEN x, long bit)
 Line 1467  lindep2(GEN x, long bit)
 Line 1672  lindep2(GEN x, long bit)
      p1[lx]           = lcvtoi(gshift((GEN)re[i],bit),&e);
      if (im) p1[lx+1] = lcvtoi(gshift((GEN)im[i],bit),&e);
    }
-   p1 = (GEN)gmul(p2,lllint(p2))[1];
+   p1 = (GEN)lllint_ip(p2,100)[1];
    p1[0] = evaltyp(t_VEC) | evallg(lx);
    return gerepilecopy(av, p1);
  }
-Line 1478  lindep(GEN x, long prec)
+Line 1683  lindep(GEN x, long prec)
 Line 1478  lindep(GEN x, long prec)
 Line 1683  lindep(GEN x, long prec)
  {
    GEN *b,*be,*bn,**m,qzer;
    GEN c1,c2,c3,px,py,pxy,re,im,p1,p2,r,f,em;
-   long i,j,fl,i1, lx = lg(x), tx = typ(x), n = lx-1;
+   long i,j,fl,k, lx = lg(x), tx = typ(x), n = lx-1;
-   long av = avma, lim = stack_lim(av,1), av0,av1,tetpil;
+   gpmem_t av = avma, lim = stack_lim(av,1), av0, av1;
    const long EXP = - bit_accuracy(prec) + 2*n;
    if (! is_vec_t(tx)) err(typeer,"lindep");
-   if (lx<=2) return cgetg(1,t_VEC);
+   if (n <= 1) return cgetg(1,t_VEC);
    x = gmul(x, realun(prec));
-   re=greal(x); im=gimag(x);
+   re = greal(x);
+   im = gimag(x);
    /* independant over R ? */
-   if (lx == 3 && real_indep(re,im,bit_accuracy(prec)))
+   if (n == 2 && real_indep(re,im,bit_accuracy(prec)))
      { avma = av; return cgetg(1, t_VEC); }
+   if (EXP > -10) err(precer,"lindep");
-   qzer = new_chunk(lx);
+   qzer = cgetg(lx, t_VECSMALL);
    b = (GEN*)idmat(n);
    be= (GEN*)new_chunk(lx);
    bn= (GEN*)new_chunk(lx);
    m = (GEN**)new_chunk(lx);
    for (i=1; i<=n; i++)
    {
-     bn[i]=cgetr(prec+1);
+     bn[i] = cgetr(prec+1);
-     be[i]=cgetg(lx,t_COL);
+     be[i] = cgetg(lx,t_COL);
-     m[i] = (GEN*)new_chunk(lx);
+     m[i]  = (GEN*)new_chunk(lx);
-     for (j=1; j<i ; j++) m[i][j]=cgetr(prec+1);
+     for (j=1; j<i ; j++) m[i][j] = cgetr(prec+1);
-     for (j=1; j<=n; j++) be[i][j]=lgetr(prec+1);
+     for (j=1; j<=n; j++) be[i][j]= lgetr(prec+1);
    }
-   px=sqscal(re);
+   px = sqscal(re);
-   py=sqscal(im); pxy=gscal(re,im);
+   py = sqscal(im); pxy = gscal(re,im);
-   p1=mpsub(mpmul(px,py),gsqr(pxy));
+   p1 = mpsub(mpmul(px,py), gsqr(pxy));
-   if (quazero(re)) { re=im; px=py; fl=1; } else fl=quazero(p1);
+   if (quazero(px)) { re = im; px = py; fl = 1; } else fl = quazero(p1);
    av0 = av1 = avma;
    for (i=1; i<=n; i++)
    {
      p2 = gscal(b[i],re);
-     if (fl) p2=gmul(gdiv(p2,px),re);
+     if (fl) p2 = gmul(gdiv(p2,px),re);
      else
      {
        GEN p5,p6,p7;
-       p5=gscal(b[i],im);
+       p5 = gscal(b[i],im);
-       p6=gdiv(gsub(gmul(py,p2),gmul(pxy,p5)),p1);
+       p6 = gdiv(gsub(gmul(py,p2),gmul(pxy,p5)), p1);
-       p7=gdiv(gsub(gmul(px,p5),gmul(pxy,p2)),p1);
+       p7 = gdiv(gsub(gmul(px,p5),gmul(pxy,p2)), p1);
-       p2=gadd(gmul(p6,re),gmul(p7,im));
+       p2 = gadd(gmul(p6,re), gmul(p7,im));
      }
-     if (tx!=t_COL) p2=gtrans(p2);
+     if (tx != t_COL) p2 = gtrans(p2);
-     p2=gsub(b[i],p2);
+     p2 = gsub(b[i],p2);
      for (j=1; j<i; j++)
-       if (qzer[j]) affrr(bn[j],m[i][j]);
+       if (qzer[j]) affrr(bn[j], m[i][j]);
        else
        {
-         gdivz(gscal(b[i],be[j]),bn[j],m[i][j]);
+         gdivz(gscal(b[i],be[j]),bn[j], m[i][j]);
-         p2=gsub(p2,gmul(m[i][j],be[j]));
+         p2 = gsub(p2, gmul(m[i][j],be[j]));
        }
-     gaffect(p2,be[i]); affrr(sqscal(be[i]),bn[i]);
+     gaffect(p2,          be[i]);
-     qzer[i]=quazero(bn[i]); avma=av1;
+     affrr(sqscal(be[i]), bn[i]);
+     qzer[i] = quazero(bn[i]); avma = av1;
    }
    while (qzer[n])
    {
      long e;
-     if (DEBUGLEVEL>9)
+     if (DEBUGLEVEL>8)
      {
        fprintferr("qzer[%ld]=%ld\n",n,qzer[n]);
-       for (i1=1; i1<=n; i1++)
+       for (k=1; k<=n; k++)
-         for (i=1; i<i1; i++) output(m[i1][i]);
+         for (i=1; i<k; i++) output(m[k][i]);
      }
      em=bn[1]; j=1;
      for (i=2; i<n; i++)
-Line 1547  lindep(GEN x, long prec)
+Line 1755  lindep(GEN x, long prec)
 Line 1547  lindep(GEN x, long prec)
 Line 1755  lindep(GEN x, long prec)
        p1=shiftr(bn[i],i);
        if (cmprr(p1,em)>0) { em=p1; j=i; }
      }
-     i=j; i1=i+1;
+     i=j; k=i+1;
-     avma = av1; r = grndtoi(m[i1][i], &e);
+     avma = av1; r = grndtoi(m[k][i], &e);
      if (e >= 0) err(precer,"lindep");
      r  = negi(r);
-     p1 = ZV_lincomb(gun,r, b[i1],b[i]);
+     p1 = ZV_lincomb(gun,r, b[k],b[i]);
+     b[k] = b[i];
+     b[i]  = p1;
      av1 = avma;
-     b[i1]=b[i]; b[i]=p1; f=addir(r,m[i1][i]);
+     f = addir(r,m[k][i]);
      for (j=1; j<i; j++)
        if (!qzer[j])
        {
-         p1=mpadd(m[i1][j],mulir(r,m[i][j]));
+         p1 = mpadd(m[k][j], mulir(r,m[i][j]));
-         affrr(m[i][j],m[i1][j]); mpaff(p1,m[i][j]);
+         affrr(m[i][j],m[k][j]);
+         mpaff(p1,m[i][j]);
        }
-     c1=addrr(bn[i1],mulrr(gsqr(f),bn[i]));
+     c1 = addrr(bn[k], mulrr(gsqr(f),bn[i]));
      if (!quazero(c1))
      {
-       c2=divrr(mulrr(bn[i],f),c1); affrr(c2,m[i1][i]);
+       c2 = divrr(mulrr(bn[i],f),c1); affrr(c2, m[k][i]);
-       c3=divrr(bn[i1],c1); mulrrz(c3,bn[i],bn[i1]);
+       c3 = divrr(bn[k],c1);
-       affrr(c1,bn[i]); qzer[i1]=quazero(bn[i1]); qzer[i]=0;
+       mulrrz(c3,bn[i], bn[k]); qzer[k] = quazero(bn[k]);
+       affrr(c1,        bn[i]); qzer[i] = 0;
        for (j=i+2; j<=n; j++)
        {
-         p1=addrr(mulrr(m[j][i1],c3),mulrr(m[j][i],c2));
+         p1 = addrr(mulrr(m[j][k],c3), mulrr(m[j][i],c2));
-         subrrz(m[j][i],mulrr(f,m[j][i1]), m[j][i1]);
+         subrrz(m[j][i],mulrr(f,m[j][k]), m[j][k]);
-         affrr(p1,m[j][i]);
+         affrr(p1, m[j][i]);
        }
      }
      else
      {
-       qzer[i1]=qzer[i]; qzer[i]=1;
+       affrr(bn[i], bn[k]); qzer[k] = qzer[i];
-       affrr(bn[i],bn[i1]); affrr(c1,bn[i]);
+       affrr(c1,    bn[i]); qzer[i] = 1;
-       for (j=i+2; j<=n; j++) affrr(m[j][i],m[j][i1]);
+       for (j=i+2; j<=n; j++) affrr(m[j][i], m[j][k]);
      }
      if (low_stack(lim, stack_lim(av,1)))
      {
-Line 1586  lindep(GEN x, long prec)
+Line 1798  lindep(GEN x, long prec)
 Line 1586  lindep(GEN x, long prec)
 Line 1798  lindep(GEN x, long prec)
        av1 = avma;
      }
    }
-   p1=cgetg(lx,t_COL); p1[n]=un; for (i=1; i<n; i++) p1[i]=zero;
+   p1 = cgetg(lx,t_COL); p1[n] = un; for (i=1; i<n; i++) p1[i] = zero;
    p2 = (GEN)b; p2[0] = evaltyp(t_MAT) | evallg(lx);
-   p2=gauss(gtrans(p2),p1); tetpil=avma;
+   p2 = gauss(gtrans_i(p2),p1);
-   return gerepile(av,tetpil,gtrans(p2));
+   return gerepileupto(av, gtrans(p2));
  }
- /* x is a vector of elts of a p-adic field */
+ /* PSLQ Programs */
- GEN
- plindep(GEN x)
+ typedef struct {
+   long vmind, t12, t1234, reda, fin;
+   long ct;
+ } pslq_timer;
+ /* WARNING: for double ** matrices, A[i] is the i-th ROW of A */
+ typedef struct {
+   double *y, **H, **A, **B;
+   double *W; /* scratch space */
+   long n;
+   pslq_timer *T;
+ } pslqL2_M;
+ typedef struct {
+   GEN y, H, A, B;
+   long n, EXP;
+   int flreal;
+   pslq_timer *T;
+ } pslq_M;
+ void
+ init_dalloc()
+ { /* correct alignment for dalloc */
+   avma -= avma % sizeof(double);
+   if (avma < bot) err(errpile);
+ }
+ double *
+ dalloc(size_t n)
  {
-   long av = avma,i,j, prec = VERYBIGINT, lx = lg(x)-1, ly,v;
+   if (avma - bot < n) err(errpile);
-   GEN p = NULL, pn,p1,m,a;
+   avma -= n; return (double*)avma;
+ }
-   if (lx < 2) return cgetg(1,t_VEC);
+ static double
-   for (i=1; i<=lx; i++)
+ conjd(double x) { return x; }
+ static double
+ sqrd(double a) { return a*a; }
+ static void
+ redall(pslq_M *M, long i, long jsup)
+ {
+   long j, k, n = M->n;
+   GEN t,b;
+   GEN A = M->A, B = M->B, H = M->H, y = M->y;
+   const GEN p1 = (GEN)B[i];
+   for (j=jsup; j>=1; j--)
    {
-     p1 = (GEN)x[i];
+ /* FIXME: use grndtoi */
-     if (typ(p1) != t_PADIC) continue;
+     t = ground(gdiv(gcoeff(H,i,j), gcoeff(H,j,j)));
+     if (gcmp0(t)) continue;
-     j = precp(p1); if (j < prec) prec = j;
+     b = (GEN)B[j];
-     if (!p) p = (GEN)p1[2];
+     y[j] = ladd((GEN)y[j], gmul(t,(GEN)y[i]));
-     else if (!egalii(p, (GEN)p1[2]))
+     for (k=1; k<=j; k++)
-       err(talker,"inconsistent primes in plindep");
+       coeff(H,i,k) = lsub(gcoeff(H,i,k), gmul(t,gcoeff(H,j,k)));
+     for (k=1; k<=n; k++)
+     {
+       coeff(A,i,k) = lsub(gcoeff(A,i,k), gmul(t,gcoeff(A,j,k)));
+       b[k] = ladd((GEN)b[k], gmul(t,(GEN)p1[k]));
+     }
    }
-   if (!p) err(talker,"not a p-adic vector in plindep");
+ }
-   v = ggval(x,p); pn = gpowgs(p,prec);
-   if (v != 0) x = gmul(x, gpowgs(p, -v));
-   x = lift_intern(gmul(x, gmodulcp(gun, pn)));
-   ly = 2*lx - 1;
+ static void
-   m = cgetg(ly+1,t_MAT);
+ redallbar(pslqL2_M *Mbar, long i, long jsup)
-   for (j=1; j<=ly; j++)
+ {
+   long j, k, n = Mbar->n;
+   double t;
+   double *hi = Mbar->H[i], *ai = Mbar->A[i], *hj, *aj;
+ #if 0
+ fprintferr("%ld:\n==\n",i);
+ #endif
+   for (j=jsup; j>=1; j--)
    {
-     p1 = cgetg(lx+1, t_COL); m[j] = (long)p1;
+     hj = Mbar->H[j];
-     for (i=1; i<=lx; i++) p1[i] = zero;
+     t = floor(0.5 + hi[j] / hj[j]);
+     if (!t) continue;
+ #if 0
+ fprintferr("%15.15e ",t);
+ #endif
+     aj = Mbar->A[j];
+     Mbar->y[j] += t * Mbar->y[i];
+     for (k=1; k<=j; k++) hi[k] -= t * hj[k];
+     for (k=1; k<=n; k++) {
+       ai[k]         -= t * aj[k];
+       Mbar->B[k][j] += t * Mbar->B[k][i];
+     }
+ #if 0
+ fprintferr("  %ld:\n",j); dprintmat(Mbar->H,n,n-1);
+ #endif
    }
-   a = negi((GEN)x[1]);
+ }
-   for (i=1; i<lx; i++)
+ static long
+ vecabsminind(GEN v)
+ {
+   long l = lg(v), m = 1, i;
+   GEN t, la = mpabs((GEN)v[1]);
+   for (i=2; i<l; i++)
    {
-     coeff(m,1+i,i) = (long)a;
+     t = mpabs((GEN)v[i]);
-     coeff(m,1  ,i) = x[i+1];
+     if (mpcmp(t,la) < 0) { la = t; m = i; }
    }
-   for (i=1; i<=lx; i++) coeff(m,i,lx-1+i) = (long)pn;
+   return m;
-   p1 = lllint(m);
-   return gerepileupto(av, gmul(m, (GEN)p1[1]));
  }
- GEN
+ static long
- algdep0(GEN x, long n, long bit, long prec)
+ vecmaxind(GEN v)
  {
-   long tx=typ(x),av,i,k;
+   long l = lg(v), m = 1, i;
-   GEN y,p1;
+   GEN t, la = (GEN)v[1];
+   for (i=2; i<l; i++)
+   {
+     t = (GEN)v[i];
+     if (mpcmp(t,la) > 0) { la = t; m = i; }
+   }
+   return m;
+ }
-   if (! is_scalar_t(tx)) err(typeer,"algdep0");
+ static long
-   if (tx==t_POLMOD) { y=forcecopy((GEN)x[1]); setvarn(y,0); return y; }
+ vecmaxindbar(double *v, long n)
-   if (gcmp0(x)) return gzero;
+ {
-   if (!n) return gun;
+   long m = 1, i;
+   double la = v[1];
+   for (i=2; i<=n; i++)
+     if (v[i] > la) { la = v[i]; m = i; }
+   return m;
+ }
-   av=avma; p1=cgetg(n+2,t_COL);
+ static GEN
-   p1[1] = un;
+ maxnorml2(pslq_M *M)
-   p1[2] = (long)x; /* n >= 1 */
+ {
-   for (i=3; i<=n+1; i++) p1[i]=lmul((GEN)p1[i-1],x);
+   long n = M->n, i, j;
-   if (typ(x) == t_PADIC)
+   GEN ma = gzero, s;
-     p1 = plindep(p1);
-   else
-     p1 = bit? lindep2(p1,bit): lindep(p1,prec);
-   if (lg(p1) < 2)
-     err(talker,"higher degree than expected in algdep");
-   y=cgetg(n+3,t_POL);
+   for (i=1; i<=n; i++)
-   y[1] = evalsigne(1) | evalvarn(0);
+   {
-   k=1; while (gcmp0((GEN)p1[k])) k++;
+     s = gzero;
-   for (i=0; i<=n+1-k; i++) y[i+2] = p1[k+i];
+     for (j=1; j<n; j++) s = gadd(s, gnorm(gcoeff(M->H,i,j)));
-   normalizepol_i(y, n+4-k);
+     ma = gmax(ma, s);
-   y = (gsigne(leading_term(y)) > 0)? gcopy(y): gneg(y);
+   }
-   return gerepileupto(av,y);
+   return gsqrt(ma, DEFAULTPREC);
  }
- GEN
+ static void
- algdep2(GEN x, long n, long bit)
+ init_timer(pslq_timer *T)
  {
-   return algdep0(x,n,bit,0);
+   T->vmind = T->t12 = T->t1234 = T->reda = T->fin = T->ct = 0;
  }
- GEN
+ static int
- algdep(GEN x, long n, long prec)
+ init_pslq(pslq_M *M, GEN x, long prec)
  {
-   return algdep0(x,n,0,prec);
+   long lx = lg(x), n = lx-1, i, j, k;
- }
+   GEN s1, s, sinv;
- /********************************************************************/
+   M->EXP = - bit_accuracy(prec) + 2*n;
- /**                                                                **/
+   M->flreal = (gexpo(gimag(x)) < M->EXP);
- /**                   INTEGRAL KERNEL (LLL REDUCED)                **/
+   if (!M->flreal)
- /**                                                                **/
+     return 0; /* FIXME */
- /********************************************************************/
+   else
+     x = greal(x);
- GEN
+   if (DEBUGLEVEL>=3) { (void)timer(); init_timer(M->T); }
- matkerint0(GEN x, long flag)
+   x = gmul(x, realun(prec)); settyp(x,t_VEC);
- {
+   M->n = n;
-   switch(flag)
+   M->A = idmat(n);
+   M->B = idmat(n);
+   s1 = cgetg(lx,t_VEC); s1[n] = lnorm((GEN)x[n]);
+   s  = cgetg(lx,t_VEC); s[n]  = (long)gabs((GEN)x[n],prec);
+   for (k=n-1; k>=1; k--)
    {
-     case 0: return kerint(x);
+     s1[k] = ladd((GEN)s1[k+1], gnorm((GEN)x[k]));
-     case 1: return kerint1(x);
+     s[k]  = (long)gsqrt((GEN)s1[k], prec);
-     case 2: return kerint2(x);
-     default: err(flagerr,"matkerint");
    }
-   return NULL; /* not reached */
+   sinv = ginv((GEN)s[1]);
+   s    = gmul(sinv,s);
+   M->y = gmul(sinv, x);
+   M->H = cgetg(n,t_MAT);
+   for (j=1; j<n; j++)
+   {
+     GEN d, c = cgetg(lx,t_COL);
+     M->H[j] = (long)c;
+     for (i=1; i<j; i++) c[i] = zero;
+     c[j] = ldiv((GEN)s[j+1],(GEN)s[j]);
+     d = gneg( gdiv((GEN)M->y[j], gmul((GEN)s[j],(GEN)s[j+1]) ));
+     for (i=j+1; i<=n; i++) c[i] = lmul(gconj((GEN)M->y[i]), d);
+   }
+   for (i=2; i<=n; i++) redall(M, i, i-1);
+   return 1;
  }
- GEN
+ static void
- kerint1(GEN x)
+ SWAP(pslq_M *M, long m)
  {
-   long av=avma,tetpil;
+   long j, n = M->n;
-   GEN p1,p2;
+   swap(M->y[m], M->y[m+1]);
+   swap(M->B[m], M->B[m+1]);
-   p1=matrixqz3(ker(x)); p2=lllint(p1); tetpil=avma;
+   for (j=1; j<=n; j++) swap(coeff(M->A,m,j), coeff(M->A,m+1,j));
-   return gerepile(av,tetpil,gmul(p1,p2));
+   for (j=1; j<n;  j++) swap(coeff(M->H,m,j), coeff(M->H,m+1,j));
  }
- GEN
+ #define dswap(x,y) { double _t=x; x=y; y=_t; }
- kerint2(GEN x)
+ #define ddswap(x,y) { double* _t=x; x=y; y=_t; }
- {
-   long lx=lg(x), i,j,av,av1;
-   GEN g,p1;
-   if (typ(x) != t_MAT) err(typeer,"kerint2");
+ static void
-   av=avma; g=cgetg(lx,t_MAT);
+ SWAPbar(pslqL2_M *M, long m)
-   for (j=1; j<lx; j++) g[j]=lgetg(lx,t_COL);
+ {
-   for (i=1; i<lx; i++)
+   long j, n = M->n;
-     for (j=1; j<=i; j++)
+   dswap(M->y[m], M->y[m+1]);
-       coeff(g,i,j) = coeff(g,j,i) = (long)gscal((GEN)x[i],(GEN)x[j]);
+   ddswap(M->A[m], M->A[m+1]);
-   g=lllgramall0(g,lll_KER); p1=lllint(g);
+   ddswap(M->H[m], M->H[m+1]);
-   av1=avma; return gerepile(av,av1,gmul(g,p1));
+   for (j=1; j<=n; j++) dswap(M->B[j][m], M->B[j][m+1]);
  }
  static GEN
- lllall0(GEN x, long flag)
+ one_step_gen(pslq_M *M, GEN tabga, long prec)
  {
-   long av0=avma,av,tetpil,lx=lg(x),i,j,k,l,n,lim,kmax;
+   GEN H = M->H, p1, t0, tinv, t1,t2,t3,t4;
-   GEN u,B,L,q,r,h,la,p1,p2,p4,fl;
+   long n = M->n, i, m;
-   if (typ(x) != t_MAT) err(typeer,"lllall0");
+   p1 = cgetg(n,t_VEC);
-   n=lx-1; if (n<=1) return lllall_trivial(x,n, flag | lll_GRAM);
+   for (i=1; i<n; i++) p1[i] = lmul((GEN)tabga[i], gabs(gcoeff(H,i,i),prec));
-   fl = new_chunk(lx);
+   m = vecmaxind(p1);
+   if (DEBUGLEVEL>=4) M->T->vmind += timer();
-   av=avma; lim=stack_lim(av,1); x=dummycopy(x);
+   SWAP(M, m);
-   B=gscalcol(gun, lx);
+   if (m <= n-2)
-   L=cgetg(lx,t_MAT);
-   for (k=lg(x[1]),j=1; j<lx; j++)
    {
-     for (i=1; i<k; i++)
+     t0 = gadd(gnorm(gcoeff(H,m,m)), gnorm(gcoeff(H,m,m+1)));
-       if (typ(gcoeff(x,i,j))!=t_INT) err(lllger4);
+     tinv = ginv(gsqrt(t0, prec));
-     fl[j] = 0; L[j] = (long)zerocol(n);
+     t1 = gmul(tinv, gcoeff(H,m,m));
-   }
+     t2 = gmul(tinv, gcoeff(H,m,m+1));
-   k=2; h=idmat(n); kmax=1;
+     if (DEBUGLEVEL>=4) M->T->t12 += timer();
-   u=sqscali((GEN)x[1]);
+     for (i=m; i<=n; i++)
-   if (signe(u)) { B[2]=(long)u; coeff(L,1,1)=un; fl[1]=1; } else B[2]=un;
-   for(;;)
-   {
-     if (k > kmax)
      {
-       kmax = k;
+       t3 = gcoeff(H,i,m);
-       for (j=1; j<=k; j++)
+       t4 = gcoeff(H,i,m+1);
-       {
+       if (M->flreal)
-         if (j==k || fl[j])
+         coeff(H,i,m) = ladd(gmul(t1,t3), gmul(t2,t4));
-         {
-           long av1 = avma;
-           u=gscali((GEN)x[k],(GEN)x[j]);
-           for (i=1; i<j; i++)
-             if (fl[i])
-               u = divii(subii(mulii((GEN)B[i+1],u),
-                               mulii(gcoeff(L,k,i),gcoeff(L,j,i))),
-                         (GEN)B[i]);
-           u = gerepileuptoint(av1, u);
-           if (j<k) coeff(L,k,j)=(long)u;
-           else
-           {
-             if (signe(u)) { B[k+1]=(long)u; coeff(L,k,k)=un; fl[k]=1; }
-             else { B[k+1]=B[k]; fl[k]=0; }
-           }
-         }
-       }
-     }
-     if (fl[k-1] && !fl[k])
-     {
-       u = shifti(gcoeff(L,k,k-1),1);
-       if (absi_cmp(u, (GEN)B[k]) > 0)
-       {
-         q = truedvmdii(addii(u,(GEN)B[k]),shifti((GEN)B[k],1), NULL);
-         r = negi(q);
-         h[k] = (long)ZV_lincomb(gun,r, (GEN)h[k],(GEN)h[k-1]);
-         x[k] = (long)ZV_lincomb(gun,r, (GEN)x[k],(GEN)x[k-1]);
-         coeff(L,k,k-1)=laddii(gcoeff(L,k,k-1),mulii(r,(GEN)B[k]));
-         for (i=1; i<k-1; i++)
-           coeff(L,k,i)=laddii(gcoeff(L,k,i),mulii(r,gcoeff(L,k-1,i)));
-       }
-       la=gcoeff(L,k,k-1); p4=sqri(la);
-       swap(h[k-1], h[k]);
-       swap(x[k-1], x[k]);
-       for (j=1; j<k-1; j++) swap(coeff(L,k-1,j), coeff(L,k,j));
-       if (signe(la))
-       {
-         p2=(GEN)B[k]; p1=divii(p4,p2);
-         B[k+1]=B[k]=(long)p1;
-         for (i=k+1; i<=kmax; i++)
-           coeff(L,i,k-1)=ldivii(mulii(la,gcoeff(L,i,k-1)),p2);
-         for (j=k+1; j<kmax; j++)
-           for (i=j+1; i<=kmax; i++)
-             coeff(L,i,j)=ldivii(mulii((GEN)p1,gcoeff(L,i,j)),p2);
-         for (i=k+2; i<=kmax+1; i++)
-           B[i]=ldivii(mulii((GEN)p1,(GEN)B[i]),p2);
-       }
        else
-       {
+         coeff(H,i,m) = ladd(gmul(gconj(t1),t3), gmul(gconj(t2),t4));
-         for (i=k+1; i<=kmax; i++)
+       coeff(H,i,m+1) = lsub(gmul(t1,t4), gmul(t2,t3));
-         { coeff(L,i,k)=coeff(L,i,k-1); coeff(L,i,k-1)=zero; }
-         B[k]=B[k-1]; fl[k]=1; fl[k-1]=0;
-       }
-       if (k>2) k--;
      }
-     else
+     if (DEBUGLEVEL>=4) M->T->t1234 += timer();
-     {
+   }
-       for (l=k-1; l>=1; l--)
+   for (i=1; i<=n-1; i++)
-       {
+     if (gexpo(gcoeff(H,i,i)) <= M->EXP) {
-         u = shifti(gcoeff(L,k,l),1);
+       m = vecabsminind(M->y); return (GEN)M->B[m];
-         if (absi_cmp(u,(GEN)B[l+1]) > 0)
-         {
-           q = truedvmdii(addii(u,(GEN)B[l+1]),shifti((GEN)B[l+1],1), NULL);
-           r = negi(q);
-           h[k] = (long)ZV_lincomb(gun,r,(GEN)h[k],(GEN)h[l]);
-           x[k] = (long)ZV_lincomb(gun,r,(GEN)x[k],(GEN)x[l]);
-           coeff(L,k,l)=laddii(gcoeff(L,k,l),mulii(r,(GEN)B[l+1]));
-           for (i=1; i<l; i++)
-             coeff(L,k,i)=laddii(gcoeff(L,k,i),mulii(r,gcoeff(L,l,i)));
-          }
-       }
-       if (++k > n) break;
      }
-     if (low_stack(lim, stack_lim(av,1)))
+   for (i=m+1; i<=n; i++) redall(M, i, min(i-1,m+1));
+   if (DEBUGLEVEL>=4) M->T->reda += timer();
+   if (gexpo(M->A) >= -M->EXP) return ginv(maxnorml2(M));
+   m = vecabsminind(M->y);
+   if (gexpo((GEN)M->y[m]) <= M->EXP) return (GEN)M->B[m];
+   if (DEBUGLEVEL>=3)
+   {
+     if (DEBUGLEVEL>=4) M->T->fin += timer();
+     M->T->ct++;
+     if ((M->T->ct&0xff) == 0)
      {
-       GEN *gptr[4];
+       if (DEBUGLEVEL == 3)
-       if(DEBUGMEM>1) err(warnmem,"lllall0");
+         fprintferr("time for ct = %ld : %ld\n",M->T->ct,timer());
-       gptr[0]=&B; gptr[1]=&L; gptr[2]=&h;
+       else
-       gptr[3]=&x; gerepilemany(av,gptr,4);
+         fprintferr("time [max,t12,loop,reds,fin] = [%ld, %ld, %ld, %ld, %ld]\n",
+                    M->T->vmind, M->T->t12, M->T->t1234, M->T->reda, M->T->fin);
      }
    }
-   tetpil=avma;
+   return NULL; /* nothing interesting */
-   return gerepile(av0,tetpil,lllgramall_finish(h,fl,flag,n));
  }
- GEN
+ static GEN
- kerint(GEN x)
+ get_tabga(int flreal, long n, long prec)
  {
-   long av=avma,av1;
+   GEN ga = mpsqrt( flreal? divrs(stor(4, prec), 3): stor(2, prec) );
-   GEN g,p1;
+   GEN tabga = cgetg(n,t_VEC);
+   long i;
-   g=lllall0(x,lll_KER); if (lg(g)==1) return g;
+   tabga[1] = (long)ga;
-   p1=lllint(g); av1=avma;
+   for (i = 2; i < n; i++) tabga[i] = lmul((GEN)tabga[i-1],ga);
-   return gerepile(av,av1,gmul(g,p1));
+   return tabga;
  }
- /********************************************************************/
+ GEN
- /**                                                                **/
+ pslq(GEN x, long prec)
- /**                        POLRED & CO.                            **/
- /**                                                                **/
- /********************************************************************/
- /* remove duplicate polynomials in y, updating a (same indexes), in place */
- static long
- remove_duplicates(GEN y, GEN a)
  {
-   long k,i, nv = lg(y), av = avma;
+   GEN tabga, p1;
-   GEN z;
+   long tx = typ(x);
+   gpmem_t av0 = avma, lim = stack_lim(av0,1), av;
+   pslq_M M;
+   pslq_timer T;
-   if (nv < 2) return nv;
+   if (! is_vec_t(tx)) err(typeer,"pslq");
-   z = new_chunk(3);
+   if (lg(x) <= 2) return cgetg(1,t_VEC);
-   z[1] = (long)y;
-   z[2] = (long)a; (void)sort_factor(z, gcmp);
+   M.T = &T; init_pslq(&M, x, prec);
-   for  (k=1, i=2; i<nv; i++)
+   if (!M.flreal) return lindep(x, prec); /* FIXME */
-     if (!gegal((GEN)y[i], (GEN)y[k]))
+   tabga = get_tabga(M.flreal, M.n, prec);
+   av = avma;
+   if (DEBUGLEVEL>=3) printf("Initialization time = %ld\n",timer());
+   for (;;)
+   {
+     if ((p1 = one_step_gen(&M, tabga, prec)))
+       return gerepilecopy(av0, p1);
+     if (low_stack(lim, stack_lim(av,1)))
      {
-       k++;
+       if(DEBUGMEM>1) err(warnmem,"pslq");
-       a[k] = a[i];
+       gerepileall(av,4,&M.y,&M.H,&M.A,&M.B);
-       y[k] = y[i];
      }
-   nv = k+1; setlg(a,nv); setlg(y,nv);
+   }
-   avma = av; return nv;
  }
- /* in place; make sure second non-zero coeff is negative (choose x or -x) */
+ /* W de longueur n-1 */
- int
- canon_pol(GEN z)
+ static double
+ dnorml2(double *W, long n, long row)
  {
-   long i,s;
+   long i;
+   double s = 0.;
-   for (i=lgef(z)-2; i>=2; i-=2)
+   for (i=row; i<n; i++) s += W[i]*W[i];
+   return s;
+ }
+ /* Hbar *= Pbar */
+ static void
+ dmatmul(pslqL2_M *Mbar, double **Pbar, long row)
+ {
+   const long n = Mbar->n; /* > row */
+   long i, j, k;
+   double s, *W = Mbar->W, **H = Mbar->H;
+   for (i = row; i <= n; i++)
    {
-     s = signe(z[i]);
+     for (j = row; j < n; j++)
-     if (s > 0)
      {
-       for (; i>=2; i-=2) z[i]=lnegi((GEN)z[i]);
+       k = row; s = H[i][k] * Pbar[k][j];
-       return -1;
+       for ( ; k < n; k++) s += H[i][k] * Pbar[k][j];
+       W[j] = s;
      }
-     if (s) return 1;
+     for (j = row; j < n; j++) H[i][j] = W[j];
    }
-   return 0;
  }
- static GEN
+ /* compute n-1 x n-1 matrix Pbar */
- pols_for_polred(GEN x, GEN base, GEN LLLbase, GEN *pta,
+ static void
-                 int (*check)(GEN, GEN), GEN arg)
+ dmakep(pslqL2_M *Mbar, double **Pbar, long row)
  {
-   long i, v = varn(x), n = lg(base);
+   long i, j, n = Mbar->n;
-   GEN p1,p2,p3,y, a = cgetg(n,t_VEC);
+   double pro, nc, *C = Mbar->H[row], *W = Mbar->W;
-   for (i=1; i<n; i++) a[i] = lmul(base,(GEN)LLLbase[i]);
+   nc = sqrt(dnorml2(C,n,row));
-   y=cgetg(n,t_VEC);
+   W[row] = (C[row] < 0) ? C[row] - nc : C[row] + nc;
-   for (i=1; i<n; i++)
+   for (i=row; i<n; i++) W[i] = C[i];
+   pro = -2.0 / dnorml2(W, n, row);
+       /* must have dnorml2(W,n,row) = 2*nc*(nc+fabs(C[1])) */
+   for (j=row; j<n; j++)
    {
-     if (DEBUGLEVEL > 2) { fprintferr("i = %ld\n",i); flusherr(); }
+     for (i=j+1; i<n; i++)
-     p1=(GEN)a[i];
+       Pbar[j][i] = Pbar[i][j] = pro * W[i] * W[j];
-     p1 = primitive_part(p1, &p3);
+     Pbar[j][j] = 1.0 + pro * W[j] * W[j];
-     p1 = caractducos(x,p1,v);
-     if (p3) p1 = ZX_rescale_pol(p1,p3);
-     p2 = modulargcd(derivpol(p1),p1);
-     p3 = leading_term(p2); if (!gcmp1(p3)) p2=gdiv(p2,p3);
-     p1 = gdiv(p1,p2);
-     if (canon_pol(p1) < 0 && pta) a[i] = (long)gneg_i((GEN)a[i]);
-     y[i] = (long)p1; if (DEBUGLEVEL > 3) outerr(p1);
-     if (check && check(arg, p1)) return p1;
    }
-   if (check) return NULL; /* no suitable polynomial found */
-   (void)remove_duplicates(y,a);
-   if (pta) *pta = a;
-   return y;
  }
- GEN
+ static void
- nf_get_T2(GEN base, GEN polr)
+ dLQdec(pslqL2_M *Mbar, double **Pbar)
  {
-   long i,j, n = lg(base);
+   long row, j, n = Mbar->n;
-   GEN p1,p2=cgetg(n,t_MAT);
-   for (i=1; i<n; i++)
+   for (row=1; row<n; row++)
    {
-     p1=cgetg(n,t_COL); p2[i]=(long)p1;
+     dmakep(Mbar, Pbar, row);
-     for (j=1; j<n; j++)
+     dmatmul(Mbar, Pbar, row);
-       p1[j] = (long) poleval((GEN)base[i],(GEN)polr[j]);
+     for (j=row+1; j<n; j++) Mbar->H[row][j] = 0.;
    }
-   return mulmat_real(gconj(gtrans(p2)),p2);
  }
- /* compute Tr(w_i w_j) */
+ void
- static GEN
+ dprintvec(double *V, long m)
- nf_get_T(GEN x, GEN w)
  {
-   long i,j,k, n = degpol(x);
+   long i;
-   GEN p1,p2,p3;
+   fprintferr("[");
-   GEN ptrace = cgetg(n+2,t_VEC);
+   for (i=1; i<m; i++) fprintferr("%15.15e, ",V[i]);
-   GEN den = cgetg(n+1,t_VEC);
+   fprintferr("%15.15e]\n",V[m]); pariflush();
-   GEN T = cgetg(n+1,t_MAT);
+ }
-   ptrace[2]=lstoi(n);
+ void
-   for (k=2; k<=n; k++)
+ dprintmat(double **M, long r, long c)
-   { /* cf polsym */
+ {
-     GEN y = x + (n-k+1);
+   long i, j;
-     p1 = mulsi(k-1,(GEN)y[2]);
+   fprintferr("[");
-     for (i=3; i<=k; i++)
+   for (i=1; i<r; i++)
-       p1 = addii(p1,mulii((GEN)y[i],(GEN)ptrace[i]));
-     ptrace[i] = lnegi(p1);
-   }
-   w = dummycopy(w);
-   for (i=1; i<=n; i++)
    {
-     den[i] = (long)denom(content((GEN)w[i]));
+     for (j=1; j<c; j++) fprintferr("%15.15e, ",M[i][j]);
-     w[i] = lmul((GEN)w[i],(GEN)den[i]);
+     fprintferr("%15.15e;\n ",M[i][c]);
    }
+   for (j=1; j<c; j++) fprintferr("%15.15e, ",M[r][j]);
-   for (i=1; i<=n; i++)
+   fprintferr("%15.15e]\n",M[r][c]); pariflush();
-   {
-     p1=cgetg(n+1,t_COL); T[i]=(long)p1;
-     for (j=1; j<i ; j++) p1[j] = coeff(T,i,j);
-     for (   ; j<=n; j++)
-     { /* cf quicktrace */
-       p2 = gres(gmul((GEN)w[i],(GEN)w[j]),x);
-       p3 = gzero;
-       for (k=lgef(p2)-1; k>1; k--)
-         p3 = addii(p3, mulii((GEN)p2[k],(GEN)ptrace[k]));
-       p1[j]=(long)divii(p3, mulii((GEN)den[i],(GEN)den[j]));
-     }
-   }
-   return T;
  }
- /* Return the base change matrix giving the an LLL-reduced basis for the
+ static long
-  * maximal order of the nf given by x. Expressed in terms of the standard
+ initializedoubles(pslqL2_M *Mbar, pslq_M *M, long prec)
-  * HNF basis (as polynomials) given in base (ignored if x is an nf)
-  */
- GEN
- LLL_nfbasis(GEN *ptx, GEN polr, GEN base, long prec)
  {
-   GEN T2,p1, x = *ptx;
+   long i, j, n = Mbar->n;
-   int totally_real,n,i;
+   GEN ypro;
+   gpmem_t av = avma;
-   if (typ(x) != t_POL)
+   ypro = gdiv(M->y, vecmax(gabs(M->y,prec)));
+   for (i=1; i<=n; i++)
    {
-     p1=checknf(x); *ptx=x=(GEN)p1[1];
+     if (gexpo((GEN)ypro[i]) < -0x3ff) return 0;
-     base=(GEN)p1[7]; n=degpol(x);
+     Mbar->y[i] = rtodbl((GEN)ypro[i]);
-     totally_real = !signe(gmael(p1,2,2));
-     T2=gmael(p1,5,3); if (totally_real) T2 = ground(T2);
    }
-   else
+   avma = av;
-   {
+   for (j=1; j<=n; j++)
-     n=degpol(x); totally_real = (!prec || sturm(x)==n);
+     for (i=1; i<=n; i++)
-     if (typ(base) != t_VEC || lg(base)-1 != n)
+     {
-       err(talker,"incorrect Zk basis in LLL_nfbasis");
+       if (i==j) Mbar->A[i][j] = Mbar->B[i][j] = 1.;
-     if (!totally_real)
+       else      Mbar->A[i][j] = Mbar->B[i][j] = 0.;
-       T2 = nf_get_T2(base,polr? polr: roots(x,prec));
+       if (j < n)
-     else
+       {
-       T2 = nf_get_T(x, base);
+         GEN h = gcoeff(M->H,i,j);
-   }
+         if (!gcmp0(h) && labs(gexpo(h)) > 0x3ff) return 0;
-   if (totally_real) return lllgramint(T2);
+         Mbar->H[i][j] = rtodbl(h);
-   for (i=1; ; i++)
+       }
-   {
+     }
-     if ((p1 = lllgramintern(T2,100,1,prec))) return p1;
+   return 1;
-     if (i == MAXITERPOL) err(accurer,"allpolred");
-     prec=(prec<<1)-2;
-     if (DEBUGLEVEL) err(warnprec,"allpolred",prec);
-     T2=nf_get_T2(base,roots(x,prec));
-   }
  }
- /* x can be a polynomial, but also an nf or a bnf */
+ /* T(arget) := S(ource) */
- /* if check != NULL, return the first polynomial pol
+ static void
-    found such that check(arg, pol) != 0; return NULL
+ storeprecdoubles(pslqL2_M *T, pslqL2_M *S)
-    if there are no such pol */
- static GEN
- allpolred0(GEN x, GEN *pta, long code, long prec,
-            int (*check)(GEN,GEN), GEN arg)
  {
-   GEN y,p1, base = NULL, polr = NULL;
+   long n = T->n, i, j;
-   long av = avma;
-   if (typ(x) == t_POL)
+   for (i=1; i<=n; i++)
    {
-     if (!signe(x)) return gcopy(x);
+     for (j=1; j<n; j++)
-     check_pol_int(x,"polred");
-     if (!gcmp1(leading_term(x)))
-       err(impl,"allpolred for nonmonic polynomials");
-     base = allbase4(x,code,&p1,NULL); /* p1 is junk */
-   }
-   else
-   {
-     long i = lg(x);
-     if (typ(x) == t_VEC && i<=4 && i>=3 && typ(x[1])==t_POL)
-     { /* polynomial + integer basis */
-       base=(GEN)x[2]; x=(GEN)x[1];
-     }
-     else
      {
-       x = checknf(x);
+       T->H[i][j] = S->H[i][j];
-       base=(GEN)x[7]; x=(GEN)x[1];
+       T->A[i][j] = S->A[i][j];
+       T->B[i][j] = S->B[i][j];
      }
+     T->A[i][n] = S->A[i][n];
+     T->B[i][n] = S->B[i][n];
+     T->y[i] = S->y[i];
    }
-   p1 = LLL_nfbasis(&x,polr,base,prec);
-   y = pols_for_polred(x,base,p1,pta,check,arg);
-   if (check)
-   {
-     if (y) return gerepileupto(av, y);
-     avma = av; return NULL;
-   }
-   if (pta)
-   {
-     GEN *gptr[2]; gptr[0]=&y; gptr[1]=pta;
-     gerepilemany(av,gptr,2); return y;
-   }
-   return gerepileupto(av,y);
  }
- static GEN
+ static long
- allpolred(GEN x, GEN *pta, long code, long prec)
+ checkentries(pslqL2_M *Mbar)
  {
-   return allpolred0(x,pta,code,prec,NULL,NULL);
+   long n = Mbar->n, i, j;
- }
+   double *p1, *p2;
- GEN
+   for (i=1; i<=n; i++)
- polred0(GEN x,long flag, GEN p, long prec)
- {
-   GEN y;
-   long smll = (flag & 1);
-   if (p && !gcmp0(p)) smll=(long) p; /* factored polred */
-   if (flag & 2) /* polred2 */
    {
-     y=cgetg(3,t_MAT);
+     if (expodb(Mbar->y[i]) < -46) return 0;
-     y[2]=(long)allpolred(x,(GEN*)(y+1),smll,prec);
+     p1 = Mbar->A[i];
-     return y;
+     p2 = Mbar->B[i];
+     for (j=1; j<=n; j++)
+       if (expodb(p1[j]) > 43 || expodb(p2[j]) > 43) return 0;
    }
-   return allpolred(x,NULL,smll,prec);
+   return 1;
  }
- GEN
+ static long
- ordred(GEN x, long prec)
+ applybar(pslq_M *M, pslqL2_M *Mbar, GEN Abargen, GEN Bbargen)
  {
-   GEN base,y;
+   long n = Mbar->n, i, j;
-   long n=degpol(x),i,av=avma,v = varn(x);
+   double *p1, *p2;
-   if (typ(x) != t_POL) err(typeer,"ordred");
-   if (!signe(x)) return gcopy(x);
-   if (!gcmp1((GEN)x[n+2])) err(impl,"ordred for nonmonic polynomials");
-   base = cgetg(n+1,t_VEC); /* power basis */
    for (i=1; i<=n; i++)
-     base[i] = lpowgs(polx[v],i-1);
+   {
-   y = cgetg(3,t_VEC);
+     p1 = Mbar->A[i];
-   y[1] = (long)x;
+     p2 = Mbar->B[i];
-   y[2] = (long)base;
+     for (j=1; j<=n; j++)
-   return gerepileupto(av, allpolred(y,NULL,0,prec));
+     {
+       if (expodb(p1[j]) >= 52 || expodb(p2[j]) >= 52) return 0;
+       coeff(Abargen,i,j) = (long)mpent(dbltor(p1[j]));
+       coeff(Bbargen,i,j) = (long)mpent(dbltor(p2[j]));
+     }
+   }
+   M->y = gmul(M->y, Bbargen);
+   M->B = gmul(M->B, Bbargen);
+   M->A = gmul(Abargen, M->A);
+   M->H = gmul(Abargen, M->H); return 1;
  }
- GEN
+ static GEN
- sum(GEN v, long a, long b)
+ checkend(pslq_M *M, long prec)
  {
-   GEN p;
+   long i, m, n = M->n;
-   long i;
-   if (a > b) return gzero;
-   p = (GEN)v[a];
-   for (i=a+1; i<=b; i++) p = gadd(p, (GEN)v[i]);
-   return p;
- }
- GEN
+   for (i=1; i<=n-1; i++)
- T2_from_embed_norm(GEN x, long r1)
+     if (gexpo(gcoeff(M->H,i,i)) <= M->EXP)
- {
+     {
-   GEN p = sum(x, 1, r1);
+       m = vecabsminind(M->y);
-   GEN q = sum(x, r1+1, lg(x)-1);
+       return (GEN)M->B[m];
-   if (q != gzero) p = gadd(p, gmul2n(q,1));
+     }
-   return p;
+   if (gexpo(M->A) >= -M->EXP)
+     return ginv( maxnorml2(M) );
+   m = vecabsminind(M->y);
+   if (gexpo((GEN)M->y[m]) <= M->EXP) return (GEN)M->B[m];
+   return NULL;
  }
  GEN
- T2_from_embed(GEN x, long r1)
+ pslqL2(GEN x, long prec)
  {
-   return T2_from_embed_norm(gnorm(x), r1);
+   GEN Abargen, Bbargen, tabga, p1;
- }
+   long lx = lg(x), tx = typ(x), n = lx-1, i, m, ctpro, flreal, flit;
+   gpmem_t av0 = avma, lim = stack_lim(av0,1), av;
+   double *tabgabar, gabar, tinvbar, t1bar, t2bar, t3bar, t4bar;
+   double **Pbar, **Abar, **Bbar, **Hbar, *ybar;
+   pslqL2_M Mbar, Mbarst;
+   pslq_M M;
+   pslq_timer T;
- /* return T2 norm of the polynomial defining nf */
+   if (! is_vec_t(tx)) err(typeer,"pslq");
- static GEN
+   if (n <= 1) return cgetg(1,t_COL);
- get_Bnf(GEN nf)
+   M.T = &T; init_pslq(&M, x, prec);
- {
+   if (!M.flreal) return lindep(x, prec);
-   return T2_from_embed((GEN)nf[6], nf_get_r1(nf));
- }
- /* characteristic pol of x */
+   flreal = M.flreal;
- static GEN
+   tabga = get_tabga(flreal, n, prec);
- get_polchar(GEN data, GEN x)
+   Abargen = idmat(n);
- {
+   Bbargen = idmat(n);
-   GEN basis_embed = (GEN)data[1];
-   GEN g = gmul(basis_embed,x);
-   return ground(roots_to_pol_r1r2(g, data[0], 0));
- }
- /* return a defining polynomial for Q(x) */
+   Mbarst.n = Mbar.n = n;
- static GEN
+   Mbar.A = Abar = (double**)new_chunk(n+1);
- get_polmin(GEN data, GEN x)
+   Mbar.B = Bbar = (double**)new_chunk(n+1);
- {
+   Mbar.H = Hbar = (double**)new_chunk(n+1);
-   GEN g = get_polchar(data,x);
+   Mbarst.A = (double**)new_chunk(n+1);
-   GEN h = modulargcd(g,derivpol(g));
+   Mbarst.B = (double**)new_chunk(n+1);
-   if (lgef(h) > 3) g = gdivexact(g,h);
+   Mbarst.H = (double**)new_chunk(n+1);
-   return g;
+   Pbar   = (double**)new_chunk(n);
- }
- /* does x generate the correct field ? */
+   tabgabar = dalloc((n+1)*sizeof(double));
- static GEN
+   Mbar.y = ybar = dalloc((n+1)*sizeof(double));
- chk_gen(GEN data, GEN x)
+   Mbarst.y = dalloc((n+1)*sizeof(double));
- {
-   long av = avma;
-   GEN g = get_polchar(data,x);
-   if (lgef(modulargcd(g,derivpol(g))) > 3) { avma=av; return NULL; }
-   if (DEBUGLEVEL>3) fprintferr("  generator: %Z\n",g);
-   return g;
- }
- /* mat = base change matrix, gram = LLL-reduced T2 matrix */
+   Mbar.W = dalloc((n+1)*sizeof(double));
- static GEN
+   for (i=1; i< n; i++)  Pbar[i] = dalloc((n+1)*sizeof(double));
- chk_gen_init(FP_chk_fun *chk, GEN nf, GEN gram, GEN mat, long *ptprec)
+   for (i=1; i<=n; i++)  Abar[i] = dalloc((n+1)*sizeof(double));
- {
+   for (i=1; i<=n; i++)  Bbar[i] = dalloc((n+1)*sizeof(double));
-   GEN P,bound,prev,x,B,data, M = gmael(nf,5,1);
+   for (i=1; i<=n; i++)  Hbar[i] = dalloc(n*sizeof(double));
-   long N = lg(nf[7]), n = N-1,i,prec,prec2;
+   for (i=1; i<=n; i++) Mbarst.A[i] = dalloc((n+1)*sizeof(double));
-   int skipfirst = 1; /* [1,0...0] --> x rational */
+   for (i=1; i<=n; i++) Mbarst.B[i] = dalloc((n+1)*sizeof(double));
+   for (i=1; i<=n; i++) Mbarst.H[i] = dalloc(n*sizeof(double));
- /* data[0] = r1
+   gabar = gtodouble((GEN)tabga[1]); tabgabar[1] = gabar;
-  * data[1] = embeddings of LLL-reduced Zk basis
+   for (i=2; i<n; i++) tabgabar[i] = tabgabar[i-1]*gabar;
-  * data[2] = LLL reduced Zk basis (in M_n(Z)) */
-   data = new_chunk(3);
-   data[0] = itos(gmael(nf,2,1));
-   data[1] = lmul(M, mat);
-   data[2] = lmul((GEN)nf[7],mat);
-   chk->data = data;
-   x = cgetg(N,t_COL);
+   av = avma;
-   bound = get_Bnf(nf); prev = NULL;
+   if (DEBUGLEVEL>=3) printf("Initialization time = %ld\n",timer());
-   for (i=1; i<N; i++) x[i]=zero;
+ RESTART:
-   for (i=2; i<N; i++)
+   flit = initializedoubles(&Mbar, &M, prec);
+   storeprecdoubles(&Mbarst, &Mbar);
+   if (flit) dLQdec(&Mbar, Pbar);
+   ctpro = 0;
+   for (;;)
    {
-     x[i] = un;
+     if (low_stack(lim, stack_lim(av,1)))
-     P = get_polmin(data,x);
-     x[i] = zero;
-     if (degpol(P) == n)
      {
-       B = gcoeff(gram,i,i);
+       if(DEBUGMEM>1) err(warnmem,"pslq");
-       if (gcmp(B,bound) < 0) bound = B ;
+       gerepileall(av,4,&M.y,&M.H,&M.A,&M.B);
      }
-     else
+     if (flit)
      {
-       if (DEBUGLEVEL>2) fprintferr("chk_gen_init: subfield %Z\n",P);
+       ctpro++;
-       if (skipfirst == i-1)
+       for (i=1; i<n; i++) Mbar.W[i] = tabgabar[i]*fabs(Hbar[i][i]);
+       m = vecmaxindbar(Mbar.W, n-1);
+       SWAPbar(&Mbar, m);
+       if (m <= n-2)
        {
-         if (prev && !gegal(prev,P))
+         tinvbar = 1.0 / sqrt(sqrd(Hbar[m][m]) + sqrd(Hbar[m][m+1]));
+         t1bar = tinvbar*Hbar[m][m];
+         t2bar = tinvbar*Hbar[m][m+1];
+         if (DEBUGLEVEL>=4) T.t12 += timer();
+         for (i=m; i<=n; i++)
+         {
+           t3bar = Hbar[i][m];
+           t4bar = Hbar[i][m+1];
+           if (flreal)
+             Hbar[i][m] = t1bar*t3bar + t2bar*t4bar;
+           else
+             Hbar[i][m] = conjd(t1bar)*t3bar + conjd(t2bar)*t4bar;
+           Hbar[i][m+1] = t1bar*t4bar - t2bar*t3bar;
+         }
+         if (DEBUGLEVEL>=4) T.t1234 += timer();
+       }
+       flit = checkentries(&Mbar);
+       if (flit)
+       {
+         storeprecdoubles(&Mbarst, &Mbar);
+         for (i=m+1; i<=n; i++) redallbar(&Mbar, i, min(i-1,m+1));
+       }
+       else
+       {
+         if (applybar(&M, &Mbar, Abargen,Bbargen))
+         {
+           if ( (p1 = checkend(&M,prec)) ) return gerepilecopy(av0, p1);
+           goto RESTART;
+         }
+         else
          {
-           if (degpol(prev) * degpol(P) > 32) continue; /* too expensive */
+           if (ctpro == 1) goto DOGEN;
-           P = (GEN)compositum(prev,P)[1];
+           storeprecdoubles(&Mbar, &Mbarst); /* restore */
-           if (degpol(P) == n) continue;
+           if (! applybar(&M, &Mbar, Abargen,Bbargen)) err(bugparier,"pslqL2");
-           if (DEBUGLEVEL>2 && lgef(P)>lgef(prev))
+           if ( (p1 = checkend(&M, prec)) ) return gerepilecopy(av0, p1);
-             fprintferr("chk_gen_init: subfield %Z\n",P);
+           goto RESTART;
          }
-         skipfirst++; prev = P;
        }
      }
+     else
+     {
+ DOGEN:
+       if ((p1 = one_step_gen(&M, tabga, prec)))
+         return gerepilecopy(av, p1);
+     }
    }
-   /* x_1,...,x_skipfirst generate a strict subfield [unless n=skipfirst=1] */
-   chk->skipfirst = skipfirst;
-   if (DEBUGLEVEL>2) fprintferr("chk_gen_init: skipfirst = %ld\n",skipfirst);
-   /* should be gexpo( [max_k C^n_k (bound/k) ^ k] ^ (1/2) ) */
-   prec2 = (1 + (((gexpo(bound)*n)/2) >> TWOPOTBITS_IN_LONG));
-   if (prec2 < 0) prec2 = 0;
-   prec = 3 + prec2;
-   prec2= nfgetprec(nf);
-   if (DEBUGLEVEL)
-     fprintferr("chk_gen_init: estimated prec = %ld (initially %ld)\n",
-                 prec, prec2);
-   if (prec > prec2) return NULL;
-   if (prec < prec2) data[1] = (long)gprec_w((GEN)data[1], prec);
-   *ptprec = prec; return bound;
  }
- static GEN
+ /* x is a vector of elts of a p-adic field */
- chk_gen_post(GEN data, GEN res)
+ GEN
+ plindep(GEN x)
  {
-   GEN basis = (GEN)data[2];
+   long i, j, prec = VERYBIGINT, lx = lg(x)-1, ly, v;
-   GEN p1 = (GEN)res[2];
+   gpmem_t av = avma;
-   long i, l = lg(p1);
+   GEN p = NULL, pn,p1,m,a;
-   for (i=1; i<l; i++) p1[i] = lmul(basis, (GEN)p1[i]);
-   return res;
- }
- /* no garbage collecting, done in polredabs0 */
+   if (lx < 2) return cgetg(1,t_VEC);
- static GEN
+   for (i=1; i<=lx; i++)
- findmindisc(GEN nf, GEN y, GEN a, GEN phimax, long flun)
- {
-   long i,k, c = lg(y);
-   GEN v,dmin,z,beta,discs = cgetg(c,t_VEC);
-   for (i=1; i<c; i++) discs[i] = labsi(discsr((GEN)y[i]));
-   v = sindexsort(discs);
-   k = v[1]; dmin = (GEN)discs[k]; z = (GEN)y[k]; beta = (GEN)a[k];
-   for (i=2; i<c; i++)
    {
-     k = v[i];
+     p1 = (GEN)x[i];
-     if (!egalii((GEN)discs[k],dmin)) break;
+     if (typ(p1) != t_PADIC) continue;
-     if (gpolcomp((GEN)y[k],z) < 0) { z = (GEN)y[k]; beta = (GEN)a[k]; }
+     j = precp(p1); if (j < prec) prec = j;
+     if (!p) p = (GEN)p1[2];
+     else if (!egalii(p, (GEN)p1[2]))
+       err(talker,"inconsistent primes in plindep");
    }
-   if (flun & nf_RAW)
+   if (!p) err(talker,"not a p-adic vector in plindep");
-   {
+   v = ggval(x,p); pn = gpowgs(p,prec);
-     y=cgetg(3,t_VEC);
+   if (v != 0) x = gmul(x, gpowgs(p, -v));
-     y[1]=lcopy(z);
+   x = lift_intern(gmul(x, gmodulcp(gun, pn)));
-     y[2]=lcopy(beta);
-   }
-   else if (phimax)
-   {
-     y=cgetg(3,t_VEC);
-     y[1]=lcopy(z);
-     beta=polymodrecip(gmodulcp(beta,(GEN)nf[1]));
-     y[2]=(long)poleval(phimax,beta);
-   }
-   else y = gcopy(z);
-   return y;
- }
- /* no garbage collecting, done in polredabs0 */
+   ly = 2*lx - 1;
- static GEN
+   m = cgetg(ly+1,t_MAT);
- storeallpols(GEN nf, GEN z, GEN a, GEN phimax, long flun)
+   for (j=1; j<=ly; j++)
- {
-   GEN p1,y,beta;
-   if (flun & nf_RAW)
    {
-     long i, c = lg(z);
+     p1 = cgetg(lx+1, t_COL); m[j] = (long)p1;
-     y=cgetg(c,t_VEC);
+     for (i=1; i<=lx; i++) p1[i] = zero;
-     for (i=1; i<c; i++)
-     {
-       p1=cgetg(3,t_VEC); y[i]=(long)p1;
-       p1[1]=lcopy((GEN)z[i]);
-       p1[2]=lcopy((GEN)a[i]);
-     }
    }
-   else if (phimax)
+   a = negi((GEN)x[1]);
+   for (i=1; i<lx; i++)
    {
-     long i, c = lg(z);
+     coeff(m,1+i,i) = (long)a;
-     beta = new_chunk(c);
+     coeff(m,1  ,i) = x[i+1];
-     for (i=1; i<c; i++)
-       beta[i] = (long)polymodrecip(gmodulcp((GEN)a[i],(GEN)nf[1]));
-     y=cgetg(c,t_VEC);
-     for (i=1; i<c; i++)
-     {
-       p1=cgetg(3,t_VEC); y[i]=(long)p1;
-       p1[1]=lcopy((GEN)z[i]);
-       p1[2]=(long)poleval(phimax,(GEN)beta[i]);
-     }
    }
-   else y = gcopy(z);
+   for (i=1; i<=lx; i++) coeff(m,i,lx-1+i) = (long)pn;
-   return y;
+   p1 = lllint(m);
+   return gerepileupto(av, gmul(m, (GEN)p1[1]));
  }
  GEN
- polredabs0(GEN x, long flun, long prec)
+ algdep0(GEN x, long n, long bit, long prec)
  {
-   long i,nv, av = avma;
+   long tx=typ(x), i, k;
-   GEN nf,v,y,a,phimax;
+   gpmem_t av;
-   GEN (*storepols)(GEN, GEN, GEN, GEN, long);
+   GEN y,p1;
-   FP_chk_fun *chk;
-   chk = (FP_chk_fun*)new_chunk(sizeof(FP_chk_fun));
+   if (! is_scalar_t(tx)) err(typeer,"algdep0");
-   chk->f         = &chk_gen;
+   if (tx==t_POLMOD) { y=forcecopy((GEN)x[1]); setvarn(y,0); return y; }
-   chk->f_init    = &chk_gen_init;
+   if (gcmp0(x)) return gzero;
-   chk->f_post    = &chk_gen_post;
+   if (!n) return gun;
-   if ((ulong)flun >= 16) err(flagerr,"polredabs");
+   av=avma; p1=cgetg(n+2,t_COL);
-   nf = initalgall0(x,nf_SMALL|nf_REGULAR,prec);
+   p1[1] = un;
-   if (lg(nf) == 3)
+   p1[2] = (long)x; /* n >= 1 */
-   {
+   for (i=3; i<=n+1; i++) p1[i]=lmul((GEN)p1[i-1],x);
-     phimax = lift_to_pol((GEN)nf[2]);
+   if (typ(x) == t_PADIC)
-     nf = (GEN)nf[1];
+     p1 = plindep(p1);
-   }
    else
-     phimax = (flun & nf_ORIG)? polx[0]: (GEN)NULL;
-   prec = nfgetprec(nf);
-   x = (GEN)nf[1];
-   if (lgef(x) == 4)
    {
-     y = _vec(polx[varn(x)]);
+     switch(bit)
-     a = _vec(gsub((GEN)y[1], (GEN)x[2]));
-   }
-   else
-   {
-     for (i=1; ; i++)
      {
-       v = fincke_pohst(nf,NULL,5000,3,prec, chk);
+       case 0: p1 = pslq(p1,prec); break;
-       if (v) break;
+       case -1: p1 = lindep(p1,prec); break;
-       if (i==MAXITERPOL) err(accurer,"polredabs0");
+       case -2: p1 = deplin(p1); break;
-       prec = (prec<<1)-2; nf = nfnewprec(nf,prec);
+       default: p1 = lindep2(p1,bit);
-       if (DEBUGLEVEL) err(warnprec,"polredabs0",prec);
      }
-     a = (GEN)v[2];
-     y = (GEN)v[1];
    }
-   nv = lg(a);
+   if ((!bit) && (typ(p1) == t_REAL))
-   for (i=1; i<nv; i++)
-     if (canon_pol((GEN)y[i]) < 0 && phimax)
-       a[i] = (long) gneg_i((GEN)a[i]);
-   nv = remove_duplicates(y,a);
-   if (DEBUGLEVEL)
-     { fprintferr("%ld minimal vectors found.\n",nv-1); flusherr(); }
-   if (nv >= 10000) flun &= (~nf_ALL); /* should not happen */
-   storepols = (flun & nf_ALL)? storeallpols: findmindisc;
-   if (DEBUGLEVEL) fprintferr("\n");
-   if (nv==1)
    {
-     y = _vec(x);
+     y = gcopy(p1); return gerepileupto(av,y);
-     a = _vec(polx[varn(x)]);
    }
-   if (varn(y[1]) != varn(x))
+   if (lg(p1) < 2)
-   {
+     err(talker,"higher degree than expected in algdep");
-     long vx = varn(x);
+   y=cgetg(n+3,t_POL);
-     for (i=1; i<nv; i++) setvarn(y[i], vx);
+   y[1] = evalsigne(1) | evalvarn(0);
-   }
+   k=1; while (gcmp0((GEN)p1[k])) k++;
-   return gerepileupto(av, storepols(nf,y,a,phimax,flun));
+   for (i=0; i<=n+1-k; i++) y[i+2] = p1[k+i];
+   (void)normalizepol_i(y, n+4-k);
+   y = (gsigne(leading_term(y)) > 0)? gcopy(y): gneg(y);
+   return gerepileupto(av,y);
  }
  GEN
- polredabsall(GEN x, long flun, long prec)
+ algdep2(GEN x, long n, long bit)
  {
-   return polredabs0(x, flun | nf_ALL, prec);
+   return algdep0(x,n,bit,0);
  }
  GEN
- polredabs(GEN x, long prec)
+ algdep(GEN x, long n, long prec)
  {
-   return polredabs0(x,nf_REGULAR,prec);
+   return algdep0(x,n,0,prec);
  }
- GEN
+ /********************************************************************/
- polredabs2(GEN x, long prec)
+ /**                                                                **/
- {
+ /**                   INTEGRAL KERNEL (LLL REDUCED)                **/
-   return polredabs0(x,nf_ORIG,prec);
+ /**                                                                **/
- }
+ /********************************************************************/
  GEN
- polredabsnored(GEN x, long prec)
+ matkerint0(GEN x, long flag)
  {
-   return polredabs0(x,nf_NORED,prec);
+   switch(flag)
+   {
+     case 0: return kerint(x);
+     case 1: return kerint1(x);
+     default: err(flagerr,"matkerint");
+   }
+   return NULL; /* not reached */
  }
  GEN
- polred(GEN x, long prec)
+ kerint1(GEN x)
  {
-   return allpolred(x,(GEN*)0,0,prec);
+   gpmem_t av = avma;
+   return gerepilecopy(av, lllint_ip(matrixqz3(ker(x)), 100));
  }
- GEN
- polredfirstpol(GEN x, long prec, int (*check)(GEN,GEN), GEN arg)
- {
-   return allpolred0(x,(GEN*)0,0,prec,check,arg);
- }
  GEN
- smallpolred(GEN x, long prec)
+ kerint(GEN x)
  {
-   return allpolred(x,(GEN*)0,1,prec);
+   gpmem_t av = avma;
+   GEN fl, junk, h = lllint_i(x, 0, 0, &junk, &fl, NULL);
+   if (h) h = lll_finish(h,fl, lll_KER);
+   else   h = lll_trivial(x, lll_KER);
+   if (lg(h)==1) { avma = av; return cgetg(1, t_MAT); }
+   return gerepilecopy(av, lllint_ip(h, 100));
  }
- GEN
- factoredpolred(GEN x, GEN p, long prec)
- {
-   return allpolred(x,(GEN*)0,(long)p,prec);
- }
- GEN
- polred2(GEN x, long prec)
- {
-   GEN y=cgetg(3,t_MAT);
-   y[2]= (long) allpolred(x,(GEN*)(y+1),0,prec);
-   return y;
- }
- GEN
- smallpolred2(GEN x, long prec)
- {
-   GEN y=cgetg(3,t_MAT);
-   y[2]= (long) allpolred(x,(GEN*)(y+1),1,prec);
-   return y;
- }
- GEN
- factoredpolred2(GEN x, GEN p, long prec)
- {
-   GEN y=cgetg(3,t_MAT);
-   y[2]= (long) allpolred(x,(GEN*)(y+1),(long)p,prec);
-   return y;
- }
- /* relative polredabs. Returns
-  * flag = 0: relative polynomial
-  * flag = 1: relative polynomial + element
-  * flag = 2: absolute polynomial */
- GEN
- rnfpolredabs(GEN nf, GEN relpol, long flag, long prec)
- {
-   GEN p1,bpol,rnf,elt,pol;
-   long va, av = avma;
-   if (typ(relpol)!=t_POL) err(typeer,"rnfpolredabs");
-   nf=checknf(nf); va = varn(relpol);
-   if (DEBUGLEVEL>1) timer2();
-   p1 = makebasis(nf, unifpol(nf,relpol,1));
-   rnf = (GEN)p1[3];
-   if (DEBUGLEVEL>1)
-   {
-     msgtimer("absolute basis");
-     fprintferr("original absolute generator: %Z\n",p1[1]);
-   }
-   p1 = polredabs0(p1, nf_RAW | nf_NORED, prec);
-   bpol = (GEN)p1[1];
-   if (DEBUGLEVEL>1) fprintferr("reduced absolute generator: %Z\n",bpol);
-   if (flag==2) return gerepileupto(av,bpol);
-   elt = rnfelementabstorel(rnf,(GEN)p1[2]);
-   p1=cgetg(3,t_VEC);
-   pol = rnfcharpoly(nf,relpol,elt,va);
-   if (!flag) p1 = pol;
-   else
-   {
-     p1[1]=(long)pol;
-     p1[2]=(long)polymodrecip(elt);
-   }
-   return gerepileupto(av,p1);
- }
  /********************************************************************/
  /**                                                                **/
  /**                              MINIM                             **/
  /**                                                                **/
  /********************************************************************/
- int addcolumntomatrix(GEN V,GEN INVP,GEN L);
  /* x is a non-empty matrix, y is a compatible VECSMALL (dimension > 0). */
  GEN
  gmul_mat_smallvec(GEN x, GEN y)
  {
-   long c = lg(x), l = lg(x[1]), av,i,j;
+   long c = lg(x), l = lg(x[1]), i, j;
+   gpmem_t av;
    GEN z=cgetg(l,t_COL), s;
    for (i=1; i<l; i++)
-Line 2520  gmul_mat_smallvec(GEN x, GEN y)
+Line 2599  gmul_mat_smallvec(GEN x, GEN y)
 Line 2520  gmul_mat_smallvec(GEN x, GEN y)
 Line 2599  gmul_mat_smallvec(GEN x, GEN y)
  GEN
  gmul_mati_smallvec(GEN x, GEN y)
  {
-   long c = lg(x), l = lg(x[1]), av,i,j;
+   long c = lg(x), l = lg(x[1]), i, j;
+   gpmem_t av;
    GEN z=cgetg(l,t_COL), s;
    for (i=1; i<l; i++)
-Line 2528  gmul_mati_smallvec(GEN x, GEN y)
+Line 2608  gmul_mati_smallvec(GEN x, GEN y)
 Line 2528  gmul_mati_smallvec(GEN x, GEN y)
 Line 2608  gmul_mati_smallvec(GEN x, GEN y)
      av = avma; s = mulis(gcoeff(x,i,1),y[1]);
      for (j=2; j<c; j++)
        if (y[j]) s = addii(s, mulis(gcoeff(x,i,j),y[j]));
-     z[i]=(long)gerepileuptoint(av,s);
+     z[i] = lpileuptoint(av,s);
    }
    return z;
  }
-Line 2537  void
+Line 2617  void
 Line 2537  void
 Line 2617  void
  minim_alloc(long n, double ***q, GEN *x, double **y,  double **z, double **v)
  {
    long i, s;
-   double **Q;
    *x = cgetg(n, t_VECSMALL);
    *q = (double**) new_chunk(n);
+   s = n * sizeof(double);
-   /* correct alignment for the following */
+   init_dalloc();
-   s = avma % sizeof(double); avma -= s;
+   *y = dalloc(s);
-   if (avma<bot) err(errpile);
+   *z = dalloc(s);
+   *v = dalloc(s);
-   s = (n * sizeof(double))/sizeof(long);
+   for (i=1; i<n; i++) (*q)[i] = dalloc(s);
-   *y = (double*) new_chunk(s);
-   *z = (double*) new_chunk(s);
-   *v = (double*) new_chunk(s); Q = *q;
-   for (i=1; i<n; i++) Q[i] = (double*) new_chunk(s);
  }
  /* Minimal vectors for the integral definite quadratic form: a.
-Line 2568  static GEN
+Line 2643  static GEN
 Line 2568  static GEN
 Line 2643  static GEN
  minim00(GEN a, GEN BORNE, GEN STOCKMAX, long flag)
  {
    GEN x,res,p1,u,r,liste,gnorme,invp,V, *gptr[2];
-   long n = lg(a), av0 = avma, av1,av,tetpil,lim, i,j,k,s,maxrank;
+   long n = lg(a), i, j, k, s, maxrank;
+   gpmem_t av0 = avma, av1, av, tetpil, lim;
    double p,maxnorm,BOUND,*v,*y,*z,**q, eps = 0.000001;
    maxrank = 0; res = V = invp = NULL; /* gcc -Wall */
-Line 2683  minim00(GEN a, GEN BORNE, GEN STOCKMAX, long flag)
+Line 2759  minim00(GEN a, GEN BORNE, GEN STOCKMAX, long flag)
 Line 2683  minim00(GEN a, GEN BORNE, GEN STOCKMAX, long flag)
 Line 2759  minim00(GEN a, GEN BORNE, GEN STOCKMAX, long flag)
      }
      else
      {
-       long av2 = avma;
+       gpmem_t av2 = avma;
        gnorme = ground(dbltor(p));
        if (gcmp(gnorme,BORNE) >= 0) avma = av2;
        else
-Line 2707  minim00(GEN a, GEN BORNE, GEN STOCKMAX, long flag)
+Line 2783  minim00(GEN a, GEN BORNE, GEN STOCKMAX, long flag)
 Line 2707  minim00(GEN a, GEN BORNE, GEN STOCKMAX, long flag)
 Line 2783  minim00(GEN a, GEN BORNE, GEN STOCKMAX, long flag)
        case min_PERF:
        {
-         long av2=avma, I=1;
+         long I=1;
+         gpmem_t av2=avma;
          for (i=1; i<=n; i++)
            for (j=i; j<=n; j++,I++) V[I] = x[i]*x[j];
-Line 2788  perf(GEN a)
+Line 2865  perf(GEN a)
 Line 2788  perf(GEN a)
 Line 2865  perf(GEN a)
  /* general program for positive definit quadratic forms (real coeffs).
   * One needs BORNE != 0; LLL reduction done in fincke_pohst().
-  * If (flag & 2) stop as soon as stockmax is reached.
+  * If (noer) return NULL on precision problems (no error).
-  * If (flag & 1) return NULL on precision problems (no error).
   * If (check != NULL consider only vectors passing the check [ assumes
   *   stockmax > 0 and we only want the smallest possible vectors ] */
  static GEN
- smallvectors(GEN a, GEN BORNE, long stockmax, long flag, FP_chk_fun *CHECK)
+ smallvectors(GEN a, GEN BORNE, long stockmax, long noer, FP_chk_fun *CHECK)
  {
-   long av,av1,lim,N,n,i,j,k,s,epsbit,prec, checkcnt = 1;
+   long N, n, i, j, k, s, epsbit, prec, checkcnt = 1;
+   gpmem_t av, av1, lim;
    GEN u,S,x,y,z,v,q,norme1,normax1,borne1,borne2,eps,p1,alpha,norms,dummy;
-   GEN (*check)(GEN,GEN) = CHECK? CHECK->f: NULL;
+   GEN (*check)(void *,GEN) = CHECK? CHECK->f: NULL;
-   GEN data = CHECK? CHECK->data: NULL;
+   void *data = CHECK? CHECK->data: NULL;
    int skipfirst = CHECK? CHECK->skipfirst: 0;
    if (DEBUGLEVEL)
      fprintferr("smallvectors looking for norm <= %Z\n",gprec_w(BORNE,3));
-   q = sqred1intern(a,flag&1);
+   q = sqred1intern(a, noer);
    if (q == NULL) return NULL;
    if (DEBUGLEVEL>5) fprintferr("q = %Z",q);
    prec = gprecision(q);
    epsbit = bit_accuracy(prec) >> 1;
-   eps = realun(prec); setexpo(eps,-epsbit);
+   eps = real2n(-epsbit, prec);
    alpha = dbltor(0.95);
    normax1 = gzero;
    borne1= gadd(BORNE,eps);
-Line 2870  smallvectors(GEN a, GEN BORNE, long stockmax, long fla
+Line 2947  smallvectors(GEN a, GEN BORNE, long stockmax, long fla
 Line 2870  smallvectors(GEN a, GEN BORNE, long stockmax, long fla
 Line 2947  smallvectors(GEN a, GEN BORNE, long stockmax, long fla
        if (low_stack(lim, stack_lim(av,1)))
        {
-         GEN *gptr[7];
          int cnt = 4;
          if(DEBUGMEM>1) err(warnmem,"smallvectors");
-         gptr[0]=&x; gptr[1]=&y; gptr[2]=&z; gptr[3]=&normax1;
          if (stockmax)
          { /* initialize to dummy values */
            GEN T = S;
-Line 2882  smallvectors(GEN a, GEN BORNE, long stockmax, long fla
+Line 2957  smallvectors(GEN a, GEN BORNE, long stockmax, long fla
 Line 2882  smallvectors(GEN a, GEN BORNE, long stockmax, long fla
 Line 2957  smallvectors(GEN a, GEN BORNE, long stockmax, long fla
          }
          if (check)
          {
-           cnt+=3; gptr[4]=&borne1; gptr[5]=&borne2; gptr[6]=&norms;
+           cnt+=3;
            for (i=s+1; i<=stockmax; i++) norms[i]=(long)dummy;
          }
-         gerepilemany(av,gptr,cnt);
+         gerepileall(av,cnt,&x,&y,&z,&normax1,&borne1,&borne2,&norms);
        }
        if (DEBUGLEVEL>3)
        {
-Line 2922  smallvectors(GEN a, GEN BORNE, long stockmax, long fla
+Line 2997  smallvectors(GEN a, GEN BORNE, long stockmax, long fla
 Line 2922  smallvectors(GEN a, GEN BORNE, long stockmax, long fla
 Line 2997  smallvectors(GEN a, GEN BORNE, long stockmax, long fla
      {
        if (check) norms[s] = (long)norme1;
        S[s] = (long)dummycopy(x);
-       if (s == stockmax && (flag&2) && check)
+       if (s == stockmax)
        {
-         long av1 = avma;
+         gpmem_t av2 = avma;
-         GEN per = sindexsort(norms);
+         GEN per;
+         if (!check) goto END;
+         per = sindexsort(norms);
          if (DEBUGLEVEL) fprintferr("sorting...\n");
          for (i=1; i<=s; i++)
          {
            long k = per[i];
            if (check(data,(GEN)S[k]))
            {
-             S[1] = S[k]; avma = av1;
+             S[1] = S[k]; avma = av2;
              borne1 = mpadd(norme1,eps);
              borne2 = mpmul(borne1,alpha);
              s = 1; checkcnt = 0; break;
-Line 2981  END:
+Line 3059  END:
 Line 2981  END:
 Line 3059  END:
    u[3] = (long)S; return u;
  }
- /* solve x~.a.x <= bound, a > 0. If check is non-NULL keep x only if check(x).
+ /* solve q(x) = x~.a.x <= bound, a > 0.
-  * flag & 1, return NULL in case of precision problems (sqred1 or lllgram)
+  * If check is non-NULL keep x only if check(x).
-  *   raise an error otherwise.
+  * If noer, return NULL in case of precision problems, raise an error otherwise.
-  * flag & 2, return as soon as stockmax vectors found.
+  * If a is a vector, assume a[1] is the LLL-reduced Cholesky form of q */
-  * If a is a number field, use its T2 matrix
+ static GEN
-  */
+ _fincke_pohst(GEN a,GEN B0,long stockmax,long noer, long PREC, FP_chk_fun *CHECK)
- GEN
- fincke_pohst(GEN a,GEN B0,long stockmax,long flag, long PREC, FP_chk_fun *CHECK)
  {
-   VOLATILE long pr,av=avma,i,j,n;
+   gpmem_t av = avma;
-   VOLATILE GEN B,nf,r,rinvtrans,u,v,res,z,vnorm,sperm,perm,uperm,gram;
+   long i,j,l, round = 0;
-   VOLATILE GEN bound = B0;
+   GEN B,r,rinvtrans,u,v,res,z,vnorm,sperm,perm,uperm,gram, bound = B0;
-   void *catcherr = NULL;
-   long prec = PREC;
-   if (DEBUGLEVEL>2) { fprintferr("entering fincke_pohst\n"); flusherr(); }
+   if (DEBUGLEVEL>2) fprintferr("entering fincke_pohst\n");
-   if (typ(a) == t_VEC) { nf = checknf(a); a = gmael(nf,5,3); } else nf = NULL;
+   if (typ(a) == t_VEC)
-   pr = gprecision(a);
-   if (pr) prec = pr; else a = gmul(a,realun(prec));
-   if (DEBUGLEVEL>2) fprintferr("first LLL: prec = %ld\n", prec);
-   if (nf && !signe(gmael(nf,2,2))) /* totally real */
    {
-     GEN T = nf_get_T((GEN)nf[1], (GEN)nf[7]);
+     r = (GEN)a[1];
-     u = lllgramint(T);
+     u = NULL;
-     prec += 2 * ((gexpo(u) + gexpo((GEN)nf[1])) >> TWOPOTBITS_IN_LONG);
-     nf = nfnewprec(nf, prec);
-     r = qf_base_change(T,u,1);
-     i = PREC + (gexpo(r) >> TWOPOTBITS_IN_LONG);
-     if (i < prec) prec = i;
-     r = gmul(r, realun(prec));
    }
    else
    {
-     u = lllgramintern(a,4,flag&1, (prec<<1)-2);
+     long prec = PREC;
-     if (!u) goto PRECPB;
+     l = lg(a);
+     if (l == 1)
+     {
+       if (CHECK) err(talker, "dimension 0 in fincke_pohst");
+       z = cgetg(4,t_VEC);
+       z[1] = z[2] = zero;
+       z[3] = lgetg(1,t_MAT); return z;
+     }
+     i = gprecision(a);
+     if (i) prec = i; else { a = gmul(a,realun(prec)); round = 1; }
+     if (DEBUGLEVEL>2) fprintferr("first LLL: prec = %ld\n", prec);
+     u = lllgramintern(a,4,noer, (prec<<1)-2);
+     if (!u) return NULL;
      r = qf_base_change(a,u,1);
+     r = sqred1intern(r,noer);
+     if (!r) return NULL;
+     for (i=1; i<l; i++)
+     {
+       GEN s = gsqrt(gcoeff(r,i,i), prec);
+       coeff(r,i,i) = (long)s;
+       for (j=i+1; j<l; j++) coeff(r,i,j) = lmul(s, gcoeff(r,i,j));
+     }
    }
-   r = sqred1intern(r,flag&1);
+   /* now r~ * r = a in LLL basis */
-   if (!r) goto PRECPB;
-   n = lg(a);
-   if (n == 1)
-   {
-     if (CHECK) err(talker, "dimension 0 in fincke_pohst");
-     avma = av; z=cgetg(4,t_VEC);
-     z[1]=z[2]=zero;
-     z[3]=lgetg(1,t_MAT); return z;
-   }
-   for (i=1; i<n; i++)
-   {
-     GEN p1 = gsqrt(gcoeff(r,i,i), prec);
-     coeff(r,i,i)=(long)p1;
-     for (j=i+1; j<n; j++)
-       coeff(r,i,j) = lmul(p1, gcoeff(r,i,j));
-   }
-   /* now r~ * r = a in approximate LLL basis */
    rinvtrans = gtrans(invmat(r));
+   if (DEBUGLEVEL>2)
+       fprintferr("final LLL: prec = %ld\n", gprecision(rinvtrans));
+   v = lllintern(rinvtrans, 100, noer, 0);
+   if (!v) return NULL;
+   rinvtrans = gmul(rinvtrans, v);
-   v = NULL;
+   v = ZM_inv(gtrans_i(v),gun);
-   for (i=1; i<6; i++) /* try to get close to a genuine LLL basis */
+   r = gmul(r,v);
-   {
+   u = u? gmul(u,v): v;
-     GEN p1;
-     if (DEBUGLEVEL>2)
-       fprintferr("final LLLs: prec = %ld, precision(rinvtrans) = %ld\n",
-                   prec,gprecision(rinvtrans));
-     p1 = lllintern(rinvtrans,flag&1, (gprecision(rinvtrans)<<1)-2);
-     if (!p1) goto PRECPB;
-     if (ishnfall(p1)) break; /* upper triangular */
-     if (v) v = gmul(v,p1); else v = p1;
-     rinvtrans = gmul(rinvtrans,p1);
-   }
-   if (i == 6) goto PRECPB; /* diverges... */
-   if (v)
+   l = lg(r);
-   {
+   vnorm = cgetg(l,t_VEC);
-     GEN u2 = ZM_inv(gtrans(v),gun);
+   for (j=1; j<l; j++) vnorm[j] = lnorml2((GEN)rinvtrans[j]);
-     r = gmul(r,u2); /* r = LLL basis now */
+   sperm = cgetg(l,t_MAT);
-     u = gmul(u,u2);
+   uperm = cgetg(l,t_MAT); perm = sindexsort(vnorm);
-   }
+   for (i=1; i<l; i++) { uperm[l-i] = u[perm[i]]; sperm[l-i] = r[perm[i]]; }
-   vnorm = cgetg(n,t_VEC);
-   for (j=1; j<n; j++) vnorm[j] = lnorml2((GEN)rinvtrans[j]);
-   sperm = cgetg(n,t_MAT);
-   uperm = cgetg(n,t_MAT); perm = sindexsort(vnorm);
-   for (i=1; i<n; i++) { uperm[n-i] = u[perm[i]]; sperm[n-i] = r[perm[i]]; }
    gram = gram_matrix(sperm);
-   B = gcoeff(gram,n-1,n-1);
+   B = gcoeff(gram,l-1,l-1);
-   if (gexpo(B) >= bit_accuracy(lg(B)-2)) goto PRECPB;
+   if (gexpo(B) >= bit_accuracy(lg(B)-2)) return NULL;
-   { /* catch precision problems (truncation error) */
+   res = NULL;
-     jmp_buf env;
+   CATCH(precer) { }
-     if (setjmp(env)) goto PRECPB;
+   TRY {
-     catcherr = err_catch(precer, env, NULL);
+     if (CHECK && CHECK->f_init)
-   }
-   if (CHECK && CHECK->f_init)
-   {
-     bound = CHECK->f_init(CHECK,nf,gram,uperm, &prec);
-     if (!bound) goto PRECPB;
-   }
-   if (!bound)
-   { /* set bound */
-     GEN x = cgetg(n,t_COL);
-     if (nf) bound = get_Bnf(nf);
-     for (i=2; i<n; i++) x[i]=zero;
-     for (i=1; i<n; i++)
      {
-       x[i] = un; B = gcoeff(gram,i,i);
+       bound = CHECK->f_init(CHECK,gram,uperm);
-       if (!bound || mpcmp(B,bound) < 0) bound = B;
+       if (!bound) err(precer,"fincke_pohst");
-       x[i] = zero;
      }
-   }
+     if (!bound) bound = gcoeff(gram,1,1);
-   if (DEBUGLEVEL>2) {fprintferr("entering smallvectors\n"); flusherr();}
+     if (DEBUGLEVEL>2) fprintferr("entering smallvectors\n");
-   for (i=1; i<n; i++)
+     for (i=1; i<l; i++)
+     {
+       res = smallvectors(gram, bound, stockmax,noer,CHECK);
+       if (!res) err(precer,"fincke_pohst");
+       if (!CHECK || bound || lg(res[2]) > 1) break;
+       if (DEBUGLEVEL>2) fprintferr("  i = %ld failed\n",i);
+     }
+   } ENDCATCH;
+   if (DEBUGLEVEL>2) fprintferr("leaving fincke_pohst\n");
+   if (CHECK || !res) return res;
+   z = cgetg(4,t_VEC);
+   z[1] = lcopy((GEN)res[1]);
+   z[2] = round? lround((GEN)res[2]): lcopy((GEN)res[2]);
+   z[3] = lmul(uperm, (GEN)res[3]); return gerepileupto(av,z);
+ }
+ GEN
+ fincke_pohst(GEN a,GEN B0,long stockmax,long flag, long PREC, FP_chk_fun *CHECK)
+ {
+   gpmem_t av = avma;
+   GEN z = _fincke_pohst(a,B0,stockmax,flag & 1,PREC,CHECK);
+   if (!z)
    {
-     res = smallvectors(gram, bound? bound: gcoeff(gram,i,i),
+     if (!(flag & 1)) err(precer,"fincke_pohst");
-                        stockmax,flag,CHECK);
+     avma = av;
-     if (!res) goto PRECPB;
-     if (!CHECK || bound || lg(res[2]) > 1) break;
-     if (DEBUGLEVEL>2) fprintferr("  i = %ld failed\n",i);
    }
-   err_leave(&catcherr); catcherr = NULL;
+   return z;
-   if (CHECK)
-   {
-     if (CHECK->f_post) res = CHECK->f_post(CHECK->data, res);
-     return res;
-   }
-   if (DEBUGLEVEL>2) {fprintferr("leaving fincke_pohst\n"); flusherr();}
-   z=cgetg(4,t_VEC);
-   z[1]=lcopy((GEN)res[1]);
-   z[2]=pr? lcopy((GEN)res[2]) : lround((GEN)res[2]);
-   z[3]=lmul(uperm, (GEN)res[3]); return gerepileupto(av,z);
- PRECPB:
-   if (catcherr) err_leave(&catcherr);
-   if (!(flag & 1))
-     err(talker,"not a positive definite form in fincke_pohst");
-   avma = av; return NULL;
  }

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>