===================================================================
RCS file: /home/cvs/OpenXM_contrib/pari-2.2/src/basemath/Attic/bibli1.c,v
retrieving revision 1.1.1.1
retrieving revision 1.2
diff -u -p -r1.1.1.1 -r1.2
--- OpenXM_contrib/pari-2.2/src/basemath/Attic/bibli1.c	2001/10/02 11:17:02	1.1.1.1
+++ OpenXM_contrib/pari-2.2/src/basemath/Attic/bibli1.c	2002/09/11 07:26:49	1.2
@@ -1,4 +1,4 @@
-/* $Id: bibli1.c,v 1.1.1.1 2001/10/02 11:17:02 noro Exp $
+/* $Id: bibli1.c,v 1.2 2002/09/11 07:26:49 noro Exp $
 
 Copyright (C) 2000  The PARI group.
 
@@ -21,15 +21,18 @@ 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 GEN makebasis(GEN nf,GEN pol);
-extern GEN caractducos(GEN p, GEN x, int v);
+extern int addcolumntomatrix(GEN V,GEN INVP,GEN L);
+extern long expodb(double x);
 
+/* 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;
@@ -44,7 +47,8 @@ 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);
@@ -59,7 +63,8 @@ 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;
@@ -73,7 +78,8 @@ 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);
@@ -85,21 +91,210 @@ 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;
-  if (!n)
+  gpmem_t av = avma;
+  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);
@@ -120,392 +315,236 @@ 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;
-  GEN y;
+  long i, k, l = lg(fl);
+  GEN y, g;
 
-  k=1; while (k<=n && !fl[k]) k++;
-  switch(flag)
+  k=1; while (k<l && !fl[k]) k++;
+  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);
-      y = gcopy(h); break;
-
-    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;
+    case lll_IM: h += k-1; h[0] = evaltyp(t_MAT) | evallg(l-k+1);
+      return h;
   }
-  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;
 }
 
-/********************************************************************/
-/**                                                                **/
-/**                          LLL with content                      **/
-/**                                                                **/
-/********************************************************************/
-
-/* real Gram matrix has coeffs X[i,j] = x[i,j]*veccon[i]*veccon[j] */
-static GEN
-lllgramintwithcontent(GEN x, GEN veccon, long flag)
+/* h[,k] += q * h[,l] */
+static void
+Zupdate_col(long k, long l, GEN q, long K, GEN h)
 {
-  long av0=avma,av,tetpil,lx=lg(x),i,j,k,l,n,lim,kmax;
-  GEN u,u2,B,lam,q,r,h,la,bb,p1,p2,p3,p4,fl,corr,corr2,newcon;
-  GEN *gptr[8];
+  GEN *hl, *hk;
+  long i;
 
-  if (typ(x) != t_MAT) err(typeer,"lllgramintwithcontent");
-  n=lx-1; if (n<=1) return lllall_trivial(x,n,flag);
-  if (lg((GEN)x[1])!=lx) err(mattype1,"lllgramintwithcontent");
-  fl = new_chunk(lx);
+  if (!h) return;
+  hl = (GEN*)h[l]; hk = (GEN*)h[k];
+  if (is_pm1(q))
+  {
+    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);
-  x=dummycopy(x); veccon=dummycopy(veccon);
-  B=cgetg(lx+1,t_COL); B[1]=un;
-  /* B[i+1]=B_i; le vrai B_i est B_i*prod(1,j=1,i,veccon[j]^2) */
-
-  for (i=1; i<lx; i++) { B[i+1]=zero; fl[i]=0; }
-  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(;;)
+/* L[k,] += q * L[l,], l < k */
+static void
+Zupdate_row(long k, long l, GEN q, GEN L, GEN B)
+{
+  long i;
+  if (is_pm1(q))
   {
-    if (k>kmax)
+    if (signe(q) > 0)
     {
-      kmax=k;
-      for (j=1; j<=k; j++)
-      {
-	if (j==k || fl[j])
-	{
-	  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);
-      }
+      for (i=1;i<l; i++) coeff(L,k,i) = laddii(gcoeff(L,k,i),gcoeff(L,l,i));
+      coeff(L,k,l) = laddii(gcoeff(L,k,l), B);
+    } else {
+      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=shifti(mulii(gcoeff(lam,k,k-1),(GEN)veccon[k]),1);
-    u2=mulii((GEN)B[k],(GEN)veccon[k-1]);
-    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[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);
-    }
+  } else {
+    for(i=1;i<l;i++)  coeff(L,k,i)=laddii(gcoeff(L,k,i),mulii(q,gcoeff(L,l,i)));
+    coeff(L,k,l) = laddii(gcoeff(L,k,l), mulii(q,B));
   }
 }
 
-static GEN
-lllintwithcontent(GEN x)
+static void
+update_row(long k, long l, GEN q, GEN L)
 {
-  long lx=lg(x),i,j,av,tetpil;
-  GEN veccon,con,xred,g;
-
-  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++)
+  long i;
+  if (is_pm1(q))
   {
-    g[j]=lgetg(lx,t_COL); con=content((GEN)x[j]);
-    xred[j]=ldiv((GEN)x[j],con); veccon[j]=(long)con;
+    if (signe(q) > 0)
+    {
+      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++)
-    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));
+  coeff(L,k,l) = lmpadd(gcoeff(L,k,l), q);
 }
 
-/********************************************************************/
-/**                                                                **/
-/**                          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);
-  xl = (GEN)x[l]; hl = (GEN*)h[l];
-  xk = (GEN)x[k]; hk = (GEN*)h[k];
+  q = negi(q);
+  xl = (GEN) x[l]; xk = (GEN) x[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<l; i++) coeff(L,k,i)=laddii(gcoeff(L,k,i),gcoeff(L,l,i));
-      q = B;
+      for (i=1;i<lx;i++) coeff(x,k,i) = xk[i] = laddii((GEN)xk[i], (GEN)xl[i]);
     } 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<l; i++) coeff(L,k,i)=lsubii(gcoeff(L,k,i),gcoeff(L,l,i));
-      q = negi(B);
+      for (i=1;i<lx;i++) coeff(x,k,i) = xk[i] = lsubii((GEN)xk[i], (GEN)xl[i]);
     }
-  } else {
-    for(i=1;i<=K;i++) hk[i]=addii(hk[i],mulii(q,hl[i]));
+  } else { /* h[,k] += q* h[,l]. x[,k] += q * x[,l]. L[k,] += q* L[l,] */
     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 = grndtoi(gcoeff(L,k,l),&e);
-  GEN *hk,*hl,xk,xl;
-  if (DEBUGLEVEL>8)
-    fprintferr("error bits when rounding in lllgram: %ld\n",e);
+  GEN q = truedvmdii(addii(B,shifti(gcoeff(L,k,l),1)), shifti(B,1), NULL);
   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]; hk = (GEN*)h[k];
+  xk = (GEN)x[k]; xl = (GEN)x[l];
   if (is_pm1(q))
   {
     if (signe(q) > 0)
     {
-      for (i=1;i<=K;i++) hk[i]=addii(hk[i],hl[i]);
-      xk[k] = ladd((GEN)xk[k], (GEN)xl[k]);
-      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));
+      xk[k] = lmpadd((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]);
     } else {
-      for (i=1;i<=K;i++) hk[i]=subii(hk[i],hl[i]);
-      xk[k] = lsub((GEN)xk[k], (GEN)xl[k]);
-      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));
+      xk[k] = lmpsub((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]);
     }
   } else {
-    for (i=1;i<=K;i++) hk[i]=addii(hk[i],mulii(q,hl[i]));
-    xk[k] = ladd((GEN)xk[k], gmul(q,(GEN)xl[k]));
-    for (i=1;i<lx;i++) coeff(x,k,i)=xk[i]=ladd((GEN)xk[i], gmul(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)));
+    xk[k] = lmpadd((GEN)xk[k], mpmul(q,(GEN)xl[k]));
+    for (i=1;i<lx;i++)
+      coeff(x,k,i)=xk[i]=lmpadd((GEN)xk[i], mpmul(q,(GEN)xl[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;
-  long av,i,j,lx;
+  GEN q = round_safe(gcoeff(L,k,l));
+  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)
-    { avma = av; return 0; }
-
-  /* SWAPI(k-1,k) */
-  if (DEBUGLEVEL>3 && k==kmax) {
-    fprintferr(" (%ld)", expi(mulsi(alpha-1,sqri(Bk)))
-                       - expi(mulsi(alpha,p2))); flusherr();
+  if (fl[k])
+  {
+    GEN q;
+    if (!D) return 0; /* only swap non-kernel + kernel vector */
+    q = addii(la2, mulii((GEN)B[k-1],(GEN)B[k+1]));
+    if (cmpii(mulsi(D-1,sqri(Bk)), mulsi(D,q)) <= 0) {
+      avma = av; return 0;
+    }
+    B[k] = (long)diviiexact(q, Bk);
   }
-  swap(h[k-1], h[k]);
-  swap(x[k-1], x[k]);
-  for (j=1; j< lx; j++) swap(coeff(x,k-1,j), coeff(x,k,j));
+  /* ZSWAP(k-1,k) */
+  if (DEBUGLEVEL>3 && k==kmax)
+  { /* 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])
   {
@@ -514,445 +553,754 @@ 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);
-      av = avma = coeff(L,i,k) = (long)icopy_av(p1,(GEN)av);
+      p1 = diviiexact(p1, Bk);
+      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);
-      av = avma = coeff(L,i,k-1) = (long)icopy_av(p1,(GEN)av);
+      p1 = diviiexact(p1, Bk);
+      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);
-      B[k+1] = B[k] = (long)p1;
-      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);
+      coeff(L,i,k) = coeff(L,i,k-1);
+      coeff(L,i,k-1) = zero;
     }
-    else
-    {
-      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;
-    }
+    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));
-  if (gcmp(q, (GEN)B[k]) <= 0) { avma = av; return 0; }
+  q = mpmul((GEN)B[k-1], mpsub(delta,la2));
+  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));
-  if (gcmp0(BB)) { B[k] = 0; return 1; } /* precision pb */
+  BB = mpadd((GEN)B[k], mpmul((GEN)B[k-1],la2));
+  if (!signe(BB)) { B[k] = 0; return 1; } /* precision pb */
 
-  coeff(L,k,k-1) = ldiv(gmul(la,(GEN)B[k-1]), BB);
-  BK = gdiv((GEN)B[k], BB);
-  B[k] = lmul((GEN)B[k-1], BK);
+  coeff(L,k,k-1) = lmpdiv(mpmul(la,(GEN)B[k-1]), BB);
+  BK = mpdiv((GEN)B[k], BB);
+  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]);
-  for (j=1; j<lx ; j++) swap(coeff(x,k-1,j), coeff(x,k,j));
+  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))
   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-1) = ladd(gmul(gcoeff(L,k,k-1),gcoeff(L,i,k-1)), gmul(BK,t));
-   /*              = ladd(t, gmul(gcoeff(L,k,k-1), gcoeff(L,i,k)));
-    * would save one multiplication, but loses precision faster... */
+    coeff(L,i,k) = lmpsub(gcoeff(L,i,k-1), mpmul(la,t));
+    coeff(L,i,k-1) = lmpadd(t, mpmul(gcoeff(L,k,k-1), gcoeff(L,i,k)));
   }
   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;
-  GEN u,B,L,h,fl, *gptr[4];
+  long lx = lg(x), hx, i, j, k, l, n, kmax;
+  gpmem_t av, lim;
+  GEN B,L,h,fl;
 
-  if (typ(x) != t_MAT) err(typeer,"lllgramall");
-  n=lx-1; if (n<=1) return lllall_trivial(x,n,flag);
-  if (lg((GEN)x[1])!=lx) err(mattype1,"lllgramall");
+  if (typ(x) != t_MAT) err(typeer,"lllint");
+  n = lx-1; if (n <= 1) return NULL;
+  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);
-  B=gscalcol(gun, lx);
-  L=cgetg(lx,t_MAT);
+  av = avma; lim = stack_lim(av,1); x = dummycopy(x);
+  B = gscalcol(gun, lx);
+  L = cgetg(lx,t_MAT);
   for (j=1; j<lx; j++)
   {
-    for (i=1; i<lx; i++)
-      if (typ(gcoeff(x,i,j))!=t_INT) err(lllger4);
+    for (i=1; i<hx; i++)
+      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;
-  u=gcoeff(x,1,1); s= signe(u);
-  if (s == 0) B[2]=un;
-  else
+  h = pth? idmat(n): NULL;
+  ZincrementalGS(x, L, B, 1, fl, gram);
+  kmax = 1;
+  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])
-	{
-          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 != MARKED)
     {
-      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--)
-      {
-        REDI(x,h,L,(GEN)B[l+1],kmax,k,l);
-        if (low_stack(lim, stack_lim(av,1)))
+      if (k != MARKED)
+        for (l=k-2; l; l--)
         {
-          if(DEBUGMEM>1) err(warnmem,"lllgramall[1]");
-          gptr[0]=&B; gptr[1]=&L; gptr[2]=&h; gptr[3]=&x;
-          gerepilemany(av,gptr,4);
+          if (!gram) ZRED(k,l, x,h,L,(GEN)B[l+1],kmax);
+          else  ZRED_gram(k,l, x,h,L,(GEN)B[l+1],kmax);
+          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]");
-      gptr[0]=&B; gptr[1]=&L; gptr[2]=&h; gptr[3]=&x;
-      gerepilemany(av,gptr,4);
+      if(DEBUGMEM>1) err(warnmem,"lllint[2], kmax = %ld", kmax);
+      gerepileall(av,h?4:3,&B,&L,&x,&h);
     }
   }
   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
-lllgramall0(GEN x, long flag) { return lllgramall(x,100,flag); }
+/* Beware: this function can return NULL (dim x <= 1) */
+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 u,B,lam,q,cq,h,la,bb,p1,p2,p3,p4,fl;
-  int ps1,ps2,flc;
+  GEN 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;
+}
 
-  if (typ(x) != t_MAT) err(typeer,"lllgramallgen");
-  n=lx-1; if (n<=1) return lllall_trivial(x,n,flag);
-  if (lg((GEN)x[1])!=lx) err(mattype1,"lllgramallgen");
+static GEN
+to_mp(GEN x, long prec)
+{
+  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);
-  B=cgetg(lx+1,t_COL);
-  B[1]=un; lam=cgetg(lx,t_MAT);
-  for (j=1; j<lx; j++) lam[j]=lgetg(lx,t_COL);
-  for (i=1; i<lx; i++)
-    for (j=1; j<=i; j++)
+  u = gcoeff(L,k,l);
+  if (pslg(u) < pslg((GEN)B[l+1])) return 0;
+
+  q = gdeuc(u,(GEN)B[l+1]);
+  cq = ginv(content(q));
+  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;
-      else
-      {
-	u=gcoeff(x,i,j); tu=typ(u);
-	if (! is_scalar_t(tu) && tu != t_POL) err(lllger4);
-	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; }
-	}
-      }
+      GEN t = gcoeff(L,i,k);
+      p1 = gsub(gmul((GEN)B[k+1],gcoeff(L,i,k-1)), gmul(la,t));
+      coeff(L,i,k)   = ldiv(p1, Bk);
+      p1 = gadd(gmul(la,gcoeff(L,i,k-1)), gmul((GEN)B[k-1],t));
+      coeff(L,i,k-1) = ldiv(p1, Bk);
     }
-  k=2; h=idmat(n); flc=0;
-  for(;;)
+  }
+  else if (!gcmp0(la))
   {
-    u=gcoeff(lam,k,k-1);
-    if (pslg(u) >= pslg((GEN)B[k]))
+    p1 = gdiv(la2, Bk);
+    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;
-      h[k]=lsub(gmul(cq,(GEN)h[k]),gmul(q,(GEN)h[k-1]));
-      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)));
+      coeff(L,i,k) = coeff(L,i,k-1);
+      coeff(L,i,k-1) = zero;
     }
-    ps1 = pslg(gsqr((GEN)B[k]));
-    p3 = gmul((GEN)B[k-1],(GEN)B[k+1]);
-    la=gcoeff(lam,k,k-1); p4 = gmul(la,gcoeff(lam,k,k-1));
-    ps2=pslg(gadd(p3,p4));
-    if (fl[k-1] && (ps1>ps2 || (ps1==ps2 && flc) || !fl[k]))
+    B[k] = B[k-1]; fl[k] = 1; fl[k-1] = 0;
+  }
+  return 1;
+}
+
+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);
-      swap(h[k-1],h[k]);
-      for (j=1; j<=k-2; j++) swap(coeff(lam,k-1,j), coeff(lam,k,j));
-      if (fl[k])
-      {
-	for (i=k+1; i<=n; 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(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--;
+      u = gcoeff(x,k,j); tu = typ(u);
+      if (! is_scalar_t(tu) && tu != t_POL) err(lllger4);
+      for (i=1; i<j; i++)
+        if (fl[i])
+        {
+          u = gsub(gmul((GEN)B[i+1],u), gmul(gcoeff(L,k,i),gcoeff(L,j,i)));
+          u = gdiv(u, (GEN)B[i]);
+        }
+      coeff(L,k,j) = (long)u;
     }
+  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--)
-      {
-	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 (REDgen(k, l, h, L, B)) flc = 1;
       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;
-      gerepilemany(av,gptr,3);
+      gerepileall(av,3,&B,&L,&h);
     }
   }
-  tetpil=avma;
-  return gerepile(av0,tetpil,lllgramall_finish(h,fl,flag,n));
+  return gerepilecopy(av0, lll_finish(h,fl,flag));
 }
 
-/* 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 av,i,j;
-  A[1] = coeff(x,k,1);
-  for(j=1;j<k;)
+  long e, i, N = lg(X[k]);
+  const long prec = MEDDEFAULTPREC; /* 128 bits */
+  GEN q, tau2, rk;
+  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]);
-    j++; av = avma;
-    /* 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));
+    if (++rounds > 100) return 0; /* detect infinite loop */
+    goto UP;
   }
-  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.
- * Quality ratio = Q = (alpha-1)/alpha. Suggested value: alpha = 100
- */
+static GEN
+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;
-  long retry = 2, av = avma,lim,l,i,j,k,k1,lx=lg(x),n,kmax,KMAX;
-  long last_prec;
+  GEN delta, X, Xs, h, R, Q, E, F;
+  long j, kmax = 1, k, N = lg(X0);
+  gpmem_t lim, av, av0 = avma;
+  int retry = 0;
 
-  if (typ(x) != t_MAT) err(typeer,"lllgram");
-  n=lx-1; if (n && lg((GEN)x[1])!=lx) err(mattype1,"lllgram");
-  if (n<=1) return idmat(n);
-  xinit = x; lim = stack_lim(av,1);
-  for (k=2,j=1; j<lx; j++)
+  delta = stor(D-1, DEFAULTPREC);
+  delta = divrs(delta,D);
+  E  = vecsmall_const(N-1, 0);
+  F  = vecsmall_const(N-1, 0);
+
+  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];
-    for (i=1; i<=j; i++)
-      if (typ(p1[i]) == t_REAL) { l = lg((GEN)p1[i]); if (l>k) k=l; }
+    gpmem_t av1;
+    if (k == 1) {
+      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 
+  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 = 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,"lllgramintern",prec);
-        x = qf_base_change(gprec_w(xinit,prec),h,1);
-        {
-          GEN *gsav[2]; gsav[0]=&h; gsav[1]=&x;
-          gerepilemany(av, gsav, 2);
-        }
-        if (DEBUGLEVEL) err(warner,"lllgramintern starting over");
+        if (DEBUGLEVEL) err(warnprec,"lllfp",prec);
+        x = gprec_w(xinit,prec);
+        if (gram) x = qf_base_change(x,h,1); else x = gmul(x,h);
+        gerepileall(av,2,&h,&x);
+        if (DEBUGLEVEL) err(warner,"lllfp 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 (gcmp0((GEN)B[1]))
+  if (gram && !incrementalGS(x, L, 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 &&
-        (2*gexpo(L1) > bit_accuracy(lg(L1)) || 2*lg(L1) < last_prec))
+    if (typ(L1) == t_REAL && expo(L1) + 20 > bit_accuracy(lg(L1)))
     {
-      last_prec = lg(L1);
       if (DEBUGLEVEL>3)
-	fprintferr("\nRecomputing Gram-Schmidt, kmax = %ld, prec was %ld\n",
-                   kmax,last_prec);
+	fprintferr("\nRecomputing Gram-Schmidt, kmax = %ld\n", kmax);
+      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 (do_SWAP(x,h,L,B,kmax,k,QR))
+    if (k != MARKED)
     {
+      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 > n) break;
+      if (k != MARKED)
+        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) 
-      {
-        if (DEBUGLEVEL > 3) fprintferr("\n");
-        err(warnmem,"lllgram");
-      }
-      gerepilemany(av,gptr,5);
+      if(DEBUGMEM>1) err(warnmem,"lllfp[2], kmax = %ld", kmax);
+      gerepileall(av,5,&B,&L,&h,&x,&Q);
     }
   }
   if (DEBUGLEVEL>3) fprintferr("\n");
-  return gerepilecopy(av, h);
+  if (MARKED && MARKED != 1) swap(h[1], h[MARKED]);
+  return gerepilecopy(av0, h);
 }
 
-static GEN
-lllgram_noerr(GEN x,long prec) { return lllgramintern(x,100,1,prec); }
+GEN
+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)
@@ -960,12 +1308,9 @@ 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 */
@@ -980,189 +1325,46 @@ 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
- * unimodulaire) d'une base LLL-reduite c_i en fonction des b_i (la base
- * reduite est c=b*u).
- */
-static GEN
-lll_proto(GEN x, GEN f(GEN, long), long prec)
+GEN
+gram_matrix(GEN b)
 {
-  long lx=lg(x),i,j,av,av1;
+  long i,j, lx = lg(b);
   GEN g;
-
-  if (typ(x) != t_MAT) err(typeer,"lll_proto");
-  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);
+  if (typ(b) != t_MAT) err(typeer,"gram");
+  g = cgetg(lx,t_MAT);
   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]);
-  av1=avma; x = f(g,prec);
-  if (!x) { avma=av; return NULL; }
-  return gerepile(av,av1,x);
+      coeff(g,i,j) = coeff(g,j,i) = (long)gscal((GEN)b[i],(GEN)b[j]);
+  }
+  return g;
 }
 
 GEN
-lllintern(GEN x,long flag,long prec)
-{
-  return lll_proto(x,flag? lllgram_noerr: lllgram,prec);
+lllgen(GEN x) {
+  gpmem_t av = avma;
+  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
@@ -1192,9 +1394,9 @@ 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);
@@ -1217,7 +1419,7 @@ 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)
@@ -1244,11 +1446,7 @@ lllintpartialall(GEN m, long flag)
 
       /* Occasionally (including final pass) do garbage collection.  */
       if (npass2 % 8 == 0 || !progress)
-      {
-        gptr[0] = &dot11; gptr[1] = &dot12;
-        gptr[2] = &dot22; gptr[3] = &tm;
-        gerepilemany(ltop1, gptr, 4);
-      }
+        gerepileall(ltop1, 4, &dot11,&dot12,&dot22,&tm);
     } /* while npass2 < 2 || progress */
 
     {
@@ -1282,8 +1480,8 @@ 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);
-        q2neg = gdivround(q2neg, det12);
+        q1neg = diviiround(q1neg, 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)),
@@ -1298,15 +1496,15 @@ lllintpartialall(GEN m, long flag)
     fprintferr("tm1 = %Z",tm1);
     fprintferr("mid = %Z",mid);
   }
-  gptr[0] = &tm1; gptr[1] = &mid;
-  gerepilemany(ltop1, gptr, 2);
+  gerepileall(ltop1,2, &tm1, &mid);
   {
    /* 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 */
-    long icol, lim, reductions, npass = 0;
+    const gpmem_t ltop3 = avma; /* Excludes region with tm1 and mid */
+    const gpmem_t lim = stack_lim(ltop3,1);
+    long icol, reductions, npass = 0;
     GEN dotprd = cgetg(ncol+1, t_MAT);
 
     tm2 = idmat(ncol);
@@ -1319,7 +1517,6 @@ 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;
@@ -1330,7 +1527,7 @@ 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),
@@ -1338,7 +1535,7 @@ 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;
@@ -1360,8 +1557,7 @@ 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;
-	  gerepilemany(ltop3, gptr, 2);
+	  gerepileall(ltop3, 2, &dotprd,&tm2);
         }
       } /* for icol */
       if (!reductions) break;
@@ -1423,16 +1619,24 @@ 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);
-  if (bit>0) return lindep2(x,bit);
-  return deplin(x);
+  switch (bit)
+  {
+    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
@@ -1446,7 +1650,8 @@ 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");
@@ -1467,7 +1672,7 @@ 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);
 }
@@ -1478,68 +1683,71 @@ 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 av = avma, lim = stack_lim(av,1), av0,av1,tetpil;
+  long i,j,fl,k, lx = lg(x), tx = typ(x), n = lx-1;
+  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);
-    be[i]=cgetg(lx,t_COL);
-    m[i] = (GEN*)new_chunk(lx);
-    for (j=1; j<i ; j++) m[i][j]=cgetr(prec+1);
-    for (j=1; j<=n; j++) be[i][j]=lgetr(prec+1);
+    bn[i] = cgetr(prec+1);
+    be[i] = cgetg(lx,t_COL);
+    m[i]  = (GEN*)new_chunk(lx);
+    for (j=1; j<i ; j++) m[i][j] = cgetr(prec+1);
+    for (j=1; j<=n; j++) be[i][j]= lgetr(prec+1);
   }
-  px=sqscal(re);
-  py=sqscal(im); pxy=gscal(re,im);
-  p1=mpsub(mpmul(px,py),gsqr(pxy));
-  if (quazero(re)) { re=im; px=py; fl=1; } else fl=quazero(p1);
+  px = sqscal(re);
+  py = sqscal(im); pxy = gscal(re,im);
+  p1 = mpsub(mpmul(px,py), gsqr(pxy));
+  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);
-      p6=gdiv(gsub(gmul(py,p2),gmul(pxy,p5)),p1);
-      p7=gdiv(gsub(gmul(px,p5),gmul(pxy,p2)),p1);
-      p2=gadd(gmul(p6,re),gmul(p7,im));
+      p5 = gscal(b[i],im);
+      p6 = gdiv(gsub(gmul(py,p2),gmul(pxy,p5)), p1);
+      p7 = gdiv(gsub(gmul(px,p5),gmul(pxy,p2)), p1);
+      p2 = gadd(gmul(p6,re), gmul(p7,im));
     }
-    if (tx!=t_COL) p2=gtrans(p2);
-    p2=gsub(b[i],p2);
+    if (tx != t_COL) p2 = gtrans(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]);
-        p2=gsub(p2,gmul(m[i][j],be[j]));
+        gdivz(gscal(b[i],be[j]),bn[j], m[i][j]);
+        p2 = gsub(p2, gmul(m[i][j],be[j]));
       }
-    gaffect(p2,be[i]); affrr(sqscal(be[i]),bn[i]);
-    qzer[i]=quazero(bn[i]); avma=av1;
+    gaffect(p2,          be[i]);
+    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 (i=1; i<i1; i++) output(m[i1][i]);
+      for (k=1; k<=n; k++)
+	for (i=1; i<k; i++) output(m[k][i]);
     }
     em=bn[1]; j=1;
     for (i=2; i<n; i++)
@@ -1547,37 +1755,41 @@ lindep(GEN x, long prec)
       p1=shiftr(bn[i],i);
       if (cmprr(p1,em)>0) { em=p1; j=i; }
     }
-    i=j; i1=i+1;
-    avma = av1; r = grndtoi(m[i1][i], &e);
+    i=j; k=i+1;
+    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]));
-        affrr(m[i][j],m[i1][j]); mpaff(p1,m[i][j]);
+        p1 = mpadd(m[k][j], mulir(r,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]);
-      c3=divrr(bn[i1],c1); mulrrz(c3,bn[i],bn[i1]);
-      affrr(c1,bn[i]); qzer[i1]=quazero(bn[i1]); qzer[i]=0;
+      c2 = divrr(mulrr(bn[i],f),c1); affrr(c2, m[k][i]);
+      c3 = divrr(bn[k],c1);
+      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));
-        subrrz(m[j][i],mulrr(f,m[j][i1]), m[j][i1]);
-        affrr(p1,m[j][i]);
+        p1 = addrr(mulrr(m[j][k],c3), mulrr(m[j][i],c2));
+        subrrz(m[j][i],mulrr(f,m[j][k]), m[j][k]);
+        affrr(p1, m[j][i]);
       }
     }
     else
     {
-      qzer[i1]=qzer[i]; qzer[i]=1;
-      affrr(bn[i],bn[i1]); affrr(c1,bn[i]);
-      for (j=i+2; j<=n; j++) affrr(m[j][i],m[j][i1]);
+      affrr(bn[i], bn[k]); qzer[k] = qzer[i];
+      affrr(c1,    bn[i]); qzer[i] = 1;
+      for (j=i+2; j<=n; j++) affrr(m[j][i], m[j][k]);
     }
     if (low_stack(lim, stack_lim(av,1)))
     {
@@ -1586,924 +1798,791 @@ 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;
-  return gerepile(av,tetpil,gtrans(p2));
+  p2 = gauss(gtrans_i(p2),p1);
+  return gerepileupto(av, gtrans(p2));
 }
 
-/* x is a vector of elts of a p-adic field */
-GEN
-plindep(GEN x)
+/* PSLQ Programs */
+
+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;
-  GEN p = NULL, pn,p1,m,a;
+  if (avma - bot < n) err(errpile);
+  avma -= n; return (double*)avma;
+}
 
-  if (lx < 2) return cgetg(1,t_VEC);
-  for (i=1; i<=lx; i++)
+static double
+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];
-    if (typ(p1) != t_PADIC) continue;
+/* FIXME: use grndtoi */
+    t = ground(gdiv(gcoeff(H,i,j), gcoeff(H,j,j)));
+    if (gcmp0(t)) continue;
 
-    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");
+    b = (GEN)B[j];
+    y[j] = ladd((GEN)y[j], gmul(t,(GEN)y[i]));
+    for (k=1; k<=j; k++)
+      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;
-  m = cgetg(ly+1,t_MAT);
-  for (j=1; j<=ly; j++)
+static void
+redallbar(pslqL2_M *Mbar, long i, long jsup)
+{
+  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;
-    for (i=1; i<=lx; i++) p1[i] = zero;
+    hj = Mbar->H[j];
+    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;
-    coeff(m,1  ,i) = x[i+1];
+    t = mpabs((GEN)v[i]);
+    if (mpcmp(t,la) < 0) { la = t; m = i; }
   }
-  for (i=1; i<=lx; i++) coeff(m,i,lx-1+i) = (long)pn;
-  p1 = lllint(m);
-  return gerepileupto(av, gmul(m, (GEN)p1[1]));
+  return m;
 }
 
-GEN
-algdep0(GEN x, long n, long bit, long prec)
+static long
+vecmaxind(GEN v)
 {
-  long tx=typ(x),av,i,k;
-  GEN y,p1;
+  long l = lg(v), m = 1, i;
+  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");
-  if (tx==t_POLMOD) { y=forcecopy((GEN)x[1]); setvarn(y,0); return y; }
-  if (gcmp0(x)) return gzero;
-  if (!n) return gun;
+static long
+vecmaxindbar(double *v, long n)
+{
+  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);
-  p1[1] = un;
-  p1[2] = (long)x; /* n >= 1 */
-  for (i=3; i<=n+1; i++) p1[i]=lmul((GEN)p1[i-1],x);
-  if (typ(x) == t_PADIC)
-    p1 = plindep(p1);
-  else
-    p1 = bit? lindep2(p1,bit): lindep(p1,prec);
-  if (lg(p1) < 2)
-    err(talker,"higher degree than expected in algdep");
+static GEN
+maxnorml2(pslq_M *M)
+{
+  long n = M->n, i, j;
+  GEN ma = gzero, s;
 
-  y=cgetg(n+3,t_POL);
-  y[1] = evalsigne(1) | evalvarn(0);
-  k=1; while (gcmp0((GEN)p1[k])) k++;
-  for (i=0; i<=n+1-k; i++) y[i+2] = p1[k+i];
-  normalizepol_i(y, n+4-k);
-  y = (gsigne(leading_term(y)) > 0)? gcopy(y): gneg(y);
-  return gerepileupto(av,y);
+  for (i=1; i<=n; i++)
+  {
+    s = gzero;
+    for (j=1; j<n; j++) s = gadd(s, gnorm(gcoeff(M->H,i,j)));
+    ma = gmax(ma, s);
+  }
+  return gsqrt(ma, DEFAULTPREC);
 }
 
-GEN
-algdep2(GEN x, long n, long bit)
+static void
+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
-algdep(GEN x, long n, long prec)
+static int
+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;
 
-/********************************************************************/
-/**                                                                **/
-/**                   INTEGRAL KERNEL (LLL REDUCED)                **/
-/**                                                                **/
-/********************************************************************/
+  M->EXP = - bit_accuracy(prec) + 2*n;
+  M->flreal = (gexpo(gimag(x)) < M->EXP);
+  if (!M->flreal)
+    return 0; /* FIXME */
+  else
+    x = greal(x);
 
-GEN
-matkerint0(GEN x, long flag)
-{
-  switch(flag)
+  if (DEBUGLEVEL>=3) { (void)timer(); init_timer(M->T); }
+  x = gmul(x, realun(prec)); settyp(x,t_VEC);
+  M->n = n;
+  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);
-    case 1: return kerint1(x);
-    case 2: return kerint2(x);
-    default: err(flagerr,"matkerint");
+    s1[k] = ladd((GEN)s1[k+1], gnorm((GEN)x[k]));
+    s[k]  = (long)gsqrt((GEN)s1[k], prec);
   }
-  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
-kerint1(GEN x)
+static void
+SWAP(pslq_M *M, long m)
 {
-  long av=avma,tetpil;
-  GEN p1,p2;
-
-  p1=matrixqz3(ker(x)); p2=lllint(p1); tetpil=avma;
-  return gerepile(av,tetpil,gmul(p1,p2));
+  long j, n = M->n;
+  swap(M->y[m], M->y[m+1]);
+  swap(M->B[m], M->B[m+1]);
+  for (j=1; j<=n; j++) swap(coeff(M->A,m,j), coeff(M->A,m+1,j));
+  for (j=1; j<n;  j++) swap(coeff(M->H,m,j), coeff(M->H,m+1,j));
 }
 
-GEN
-kerint2(GEN x)
-{
-  long lx=lg(x), i,j,av,av1;
-  GEN g,p1;
+#define dswap(x,y) { double _t=x; x=y; y=_t; }
+#define ddswap(x,y) { double* _t=x; x=y; y=_t; }
 
-  if (typ(x) != t_MAT) err(typeer,"kerint2");
-  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]);
-  g=lllgramall0(g,lll_KER); p1=lllint(g);
-  av1=avma; return gerepile(av,av1,gmul(g,p1));
+static void
+SWAPbar(pslqL2_M *M, long m)
+{
+  long j, n = M->n;
+  dswap(M->y[m], M->y[m+1]);
+  ddswap(M->A[m], M->A[m+1]);
+  ddswap(M->H[m], M->H[m+1]);
+  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 u,B,L,q,r,h,la,p1,p2,p4,fl;
+  GEN H = M->H, p1, t0, tinv, t1,t2,t3,t4;
+  long n = M->n, i, m;
 
-  if (typ(x) != t_MAT) err(typeer,"lllall0");
-  n=lx-1; if (n<=1) return lllall_trivial(x,n, flag | lll_GRAM);
-  fl = new_chunk(lx);
-
-  av=avma; lim=stack_lim(av,1); x=dummycopy(x);
-  B=gscalcol(gun, lx);
-  L=cgetg(lx,t_MAT);
-  for (k=lg(x[1]),j=1; j<lx; j++)
+  p1 = cgetg(n,t_VEC);
+  for (i=1; i<n; i++) p1[i] = lmul((GEN)tabga[i], gabs(gcoeff(H,i,i),prec));
+  m = vecmaxind(p1);
+  if (DEBUGLEVEL>=4) M->T->vmind += timer();
+  SWAP(M, m);
+  if (m <= n-2)
   {
-    for (i=1; i<k; i++)
-      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;
-  u=sqscali((GEN)x[1]);
-  if (signe(u)) { B[2]=(long)u; coeff(L,1,1)=un; fl[1]=1; } else B[2]=un;
-  for(;;)
-  {
-    if (k > kmax)
+    t0 = gadd(gnorm(gcoeff(H,m,m)), gnorm(gcoeff(H,m,m+1)));
+    tinv = ginv(gsqrt(t0, prec));
+    t1 = gmul(tinv, gcoeff(H,m,m));
+    t2 = gmul(tinv, gcoeff(H,m,m+1));
+    if (DEBUGLEVEL>=4) M->T->t12 += timer();
+    for (i=m; i<=n; i++)
     {
-      kmax = k;
-      for (j=1; j<=k; j++)
-      {
-	if (j==k || fl[j])
-	{
-          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);
-      }
+      t3 = gcoeff(H,i,m);
+      t4 = gcoeff(H,i,m+1);
+      if (M->flreal)
+        coeff(H,i,m) = ladd(gmul(t1,t3), gmul(t2,t4));
       else
-      {
-	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;
-      }
-      if (k>2) k--;
+        coeff(H,i,m) = ladd(gmul(gconj(t1),t3), gmul(gconj(t2),t4));
+      coeff(H,i,m+1) = lsub(gmul(t1,t4), gmul(t2,t3));
     }
-    else
-    {
-      for (l=k-1; l>=1; l--)
-      {
-        u = shifti(gcoeff(L,k,l),1);
-	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 (DEBUGLEVEL>=4) M->T->t1234 += timer();
+  }
+  for (i=1; i<=n-1; i++)
+    if (gexpo(gcoeff(H,i,i)) <= M->EXP) {
+      m = vecabsminind(M->y); return (GEN)M->B[m];
     }
-    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(DEBUGMEM>1) err(warnmem,"lllall0");
-      gptr[0]=&B; gptr[1]=&L; gptr[2]=&h;
-      gptr[3]=&x; gerepilemany(av,gptr,4);
+      if (DEBUGLEVEL == 3)
+        fprintferr("time for ct = %ld : %ld\n",M->T->ct,timer());
+      else
+        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 gerepile(av0,tetpil,lllgramall_finish(h,fl,flag,n));
+  return NULL; /* nothing interesting */
 }
 
-GEN
-kerint(GEN x)
+static GEN
+get_tabga(int flreal, long n, long prec)
 {
-  long av=avma,av1;
-  GEN g,p1;
-
-  g=lllall0(x,lll_KER); if (lg(g)==1) return g;
-  p1=lllint(g); av1=avma;
-  return gerepile(av,av1,gmul(g,p1));
+  GEN ga = mpsqrt( flreal? divrs(stor(4, prec), 3): stor(2, prec) );
+  GEN tabga = cgetg(n,t_VEC);
+  long i;
+  tabga[1] = (long)ga;
+  for (i = 2; i < n; i++) tabga[i] = lmul((GEN)tabga[i-1],ga);
+  return tabga;
 }
 
-/********************************************************************/
-/**                                                                **/
-/**                        POLRED & CO.                            **/
-/**                                                                **/
-/********************************************************************/
-/* remove duplicate polynomials in y, updating a (same indexes), in place */
-static long
-remove_duplicates(GEN y, GEN a)
+GEN
+pslq(GEN x, long prec)
 {
-  long k,i, nv = lg(y), av = avma;
-  GEN z;
+  GEN tabga, p1;
+  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;
-  z = new_chunk(3);
-  z[1] = (long)y;
-  z[2] = (long)a; (void)sort_factor(z, gcmp);
-  for  (k=1, i=2; i<nv; i++)
-    if (!gegal((GEN)y[i], (GEN)y[k]))
+  if (! is_vec_t(tx)) err(typeer,"pslq");
+  if (lg(x) <= 2) return cgetg(1,t_VEC);
+
+  M.T = &T; init_pslq(&M, x, prec);
+  if (!M.flreal) return lindep(x, prec); /* FIXME */
+
+  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++;
-      a[k] = a[i];
-      y[k] = y[i];
+      if(DEBUGMEM>1) err(warnmem,"pslq");
+      gerepileall(av,4,&M.y,&M.H,&M.A,&M.B);
     }
-  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) */
-int
-canon_pol(GEN z)
+/* W de longueur n-1 */
+
+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]);
-    if (s > 0)
+    for (j = row; j < n; j++)
     {
-      for (; i>=2; i-=2) z[i]=lnegi((GEN)z[i]);
-      return -1;
+      k = row; s = H[i][k] * Pbar[k][j];
+      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
-pols_for_polred(GEN x, GEN base, GEN LLLbase, GEN *pta, 
-		int (*check)(GEN, GEN), GEN arg)
+/* compute n-1 x n-1 matrix Pbar */
+static void
+dmakep(pslqL2_M *Mbar, double **Pbar, long row)
 {
-  long i, v = varn(x), n = lg(base);
-  GEN p1,p2,p3,y, a = cgetg(n,t_VEC);
+  long i, j, n = Mbar->n;
+  double pro, nc, *C = Mbar->H[row], *W = Mbar->W;
 
-  for (i=1; i<n; i++) a[i] = lmul(base,(GEN)LLLbase[i]);
-  y=cgetg(n,t_VEC);
-  for (i=1; i<n; i++)
+  nc = sqrt(dnorml2(C,n,row));
+  W[row] = (C[row] < 0) ? C[row] - nc : C[row] + nc;
+  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(); }
-    p1=(GEN)a[i];
-    p1 = primitive_part(p1, &p3);
-    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; 
+    for (i=j+1; i<n; i++)
+      Pbar[j][i] = Pbar[i][j] = pro * W[i] * W[j];
+    Pbar[j][j] = 1.0 + pro * W[j] * W[j];
   }
-  if (check) return NULL; /* no suitable polynomial found */
-  (void)remove_duplicates(y,a);
-  if (pta) *pta = a;
-  return y;
 }
 
-GEN
-nf_get_T2(GEN base, GEN polr)
+static void
+dLQdec(pslqL2_M *Mbar, double **Pbar)
 {
-  long i,j, n = lg(base);
-  GEN p1,p2=cgetg(n,t_MAT);
+  long row, j, n = Mbar->n;
 
-  for (i=1; i<n; i++)
+  for (row=1; row<n; row++)
   {
-    p1=cgetg(n,t_COL); p2[i]=(long)p1;
-    for (j=1; j<n; j++)
-      p1[j] = (long) poleval((GEN)base[i],(GEN)polr[j]);
+    dmakep(Mbar, Pbar, row);
+    dmatmul(Mbar, Pbar, row);
+    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) */
-static GEN
-nf_get_T(GEN x, GEN w)
+void
+dprintvec(double *V, long m)
 {
-  long i,j,k, n = degpol(x);
-  GEN p1,p2,p3;
-  GEN ptrace = cgetg(n+2,t_VEC);
-  GEN den = cgetg(n+1,t_VEC);
-  GEN T = cgetg(n+1,t_MAT);
+  long i;
+  fprintferr("[");
+  for (i=1; i<m; i++) fprintferr("%15.15e, ",V[i]);
+  fprintferr("%15.15e]\n",V[m]); pariflush();
+}
 
-  ptrace[2]=lstoi(n);
-  for (k=2; k<=n; k++)
-  { /* cf polsym */
-    GEN y = x + (n-k+1);
-    p1 = mulsi(k-1,(GEN)y[2]);
-    for (i=3; i<=k; i++)
-      p1 = addii(p1,mulii((GEN)y[i],(GEN)ptrace[i]));
-    ptrace[i] = lnegi(p1);
-  }
-  w = dummycopy(w);
-  for (i=1; i<=n; i++)
+void
+dprintmat(double **M, long r, long c)
+{
+  long i, j;
+  fprintferr("[");
+  for (i=1; i<r; i++)
   {
-    den[i] = (long)denom(content((GEN)w[i]));
-    w[i] = lmul((GEN)w[i],(GEN)den[i]);
+    for (j=1; j<c; j++) fprintferr("%15.15e, ",M[i][j]);
+    fprintferr("%15.15e;\n ",M[i][c]);
   }
-  
-  for (i=1; i<=n; i++)
-  {
-    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;
+  for (j=1; j<c; j++) fprintferr("%15.15e, ",M[r][j]);
+  fprintferr("%15.15e]\n",M[r][c]); pariflush();
 }
 
-/* Return the base change matrix giving the an LLL-reduced basis for the
- * maximal order of the nf given by x. Expressed in terms of the standard
- * HNF basis (as polynomials) given in base (ignored if x is an nf)
- */
-GEN
-LLL_nfbasis(GEN *ptx, GEN polr, GEN base, long prec)
+static long
+initializedoubles(pslqL2_M *Mbar, pslq_M *M, long prec)
 {
-  GEN T2,p1, x = *ptx;
-  int totally_real,n,i;
+  long i, j, n = Mbar->n;
+  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];
-    base=(GEN)p1[7]; n=degpol(x);
-    totally_real = !signe(gmael(p1,2,2));
-    T2=gmael(p1,5,3); if (totally_real) T2 = ground(T2);
+    if (gexpo((GEN)ypro[i]) < -0x3ff) return 0;
+    Mbar->y[i] = rtodbl((GEN)ypro[i]);
   }
-  else
-  {
-    n=degpol(x); totally_real = (!prec || sturm(x)==n);
-    if (typ(base) != t_VEC || lg(base)-1 != n)
-      err(talker,"incorrect Zk basis in LLL_nfbasis");
-    if (!totally_real)
-      T2 = nf_get_T2(base,polr? polr: roots(x,prec));
-    else
-      T2 = nf_get_T(x, base);
-  }
-  if (totally_real) return lllgramint(T2);
-  for (i=1; ; i++)
-  {
-    if ((p1 = lllgramintern(T2,100,1,prec))) return p1;
-    if (i == MAXITERPOL) err(accurer,"allpolred");
-    prec=(prec<<1)-2; 
-    if (DEBUGLEVEL) err(warnprec,"allpolred",prec);
-    T2=nf_get_T2(base,roots(x,prec));
-  }
+  avma = av;
+  for (j=1; j<=n; j++)
+    for (i=1; i<=n; i++)
+    {
+      if (i==j) Mbar->A[i][j] = Mbar->B[i][j] = 1.;
+      else      Mbar->A[i][j] = Mbar->B[i][j] = 0.;
+      if (j < n)
+      {
+        GEN h = gcoeff(M->H,i,j);
+        if (!gcmp0(h) && labs(gexpo(h)) > 0x3ff) return 0;
+        Mbar->H[i][j] = rtodbl(h);
+      }
+    }
+  return 1;
 }
 
-/* x can be a polynomial, but also an nf or a bnf */
-/* if check != NULL, return the first polynomial pol
-   found such that check(arg, pol) != 0; return NULL
-   if there are no such pol */
-static GEN
-allpolred0(GEN x, GEN *pta, long code, long prec, 
-	   int (*check)(GEN,GEN), GEN arg)
+/* T(arget) := S(ource) */
+static void
+storeprecdoubles(pslqL2_M *T, pslqL2_M *S)
 {
-  GEN y,p1, base = NULL, polr = NULL;
-  long av = avma;
+  long n = T->n, i, j;
 
-  if (typ(x) == t_POL)
+  for (i=1; i<=n; i++)
   {
-    if (!signe(x)) return gcopy(x);
-    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
+    for (j=1; j<n; j++)
     {
-      x = checknf(x);
-      base=(GEN)x[7]; x=(GEN)x[1];
+      T->H[i][j] = S->H[i][j];
+      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
-allpolred(GEN x, GEN *pta, long code, long prec)
+static long
+checkentries(pslqL2_M *Mbar)
 {
-  return allpolred0(x,pta,code,prec,NULL,NULL);
-}
+  long n = Mbar->n, i, j;
+  double *p1, *p2;
 
-GEN
-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 */
+  for (i=1; i<=n; i++)
   {
-    y=cgetg(3,t_MAT);
-    y[2]=(long)allpolred(x,(GEN*)(y+1),smll,prec);
-    return y;
+    if (expodb(Mbar->y[i]) < -46) return 0;
+    p1 = Mbar->A[i];
+    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
-ordred(GEN x, long prec)
+static long
+applybar(pslq_M *M, pslqL2_M *Mbar, GEN Abargen, GEN Bbargen)
 {
-  GEN base,y;
-  long n=degpol(x),i,av=avma,v = varn(x);
+  long n = Mbar->n, i, j;
+  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);
-  y[1] = (long)x;
-  y[2] = (long)base;
-  return gerepileupto(av, allpolred(y,NULL,0,prec));
+  {
+    p1 = Mbar->A[i];
+    p2 = Mbar->B[i];
+    for (j=1; j<=n; j++)
+    {
+      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
-sum(GEN v, long a, long b)
+static GEN
+checkend(pslq_M *M, long prec)
 {
-  GEN p;
-  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;
-}
+  long i, m, n = M->n;
 
-GEN
-T2_from_embed_norm(GEN x, long r1)
-{
-  GEN p = sum(x, 1, r1);
-  GEN q = sum(x, r1+1, lg(x)-1);
-  if (q != gzero) p = gadd(p, gmul2n(q,1));
-  return p;
+  for (i=1; i<=n-1; i++)
+    if (gexpo(gcoeff(M->H,i,i)) <= M->EXP)
+    {
+      m = vecabsminind(M->y);
+      return (GEN)M->B[m];
+    }
+  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 */
-static GEN 
-get_Bnf(GEN nf)
-{
-  return T2_from_embed((GEN)nf[6], nf_get_r1(nf));
-}
+  if (! is_vec_t(tx)) err(typeer,"pslq");
+  if (n <= 1) return cgetg(1,t_COL);
+  M.T = &T; init_pslq(&M, x, prec);
+  if (!M.flreal) return lindep(x, prec);
 
-/* characteristic pol of x */
-static GEN
-get_polchar(GEN data, GEN x)
-{
-  GEN basis_embed = (GEN)data[1];
-  GEN g = gmul(basis_embed,x);
-  return ground(roots_to_pol_r1r2(g, data[0], 0));
-}
+  flreal = M.flreal;
+  tabga = get_tabga(flreal, n, prec);
+  Abargen = idmat(n);
+  Bbargen = idmat(n);
 
-/* return a defining polynomial for Q(x) */
-static GEN
-get_polmin(GEN data, GEN x)
-{
-  GEN g = get_polchar(data,x);
-  GEN h = modulargcd(g,derivpol(g));
-  if (lgef(h) > 3) g = gdivexact(g,h);
-  return g;
-}
+  Mbarst.n = Mbar.n = n;
+  Mbar.A = Abar = (double**)new_chunk(n+1);
+  Mbar.B = Bbar = (double**)new_chunk(n+1);
+  Mbar.H = Hbar = (double**)new_chunk(n+1);
+  Mbarst.A = (double**)new_chunk(n+1);
+  Mbarst.B = (double**)new_chunk(n+1);
+  Mbarst.H = (double**)new_chunk(n+1);
+  Pbar   = (double**)new_chunk(n);
 
-/* does x generate the correct field ? */
-static GEN
-chk_gen(GEN data, GEN x)
-{
-  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;
-}
+  tabgabar = dalloc((n+1)*sizeof(double));
+  Mbar.y = ybar = dalloc((n+1)*sizeof(double));
+  Mbarst.y = dalloc((n+1)*sizeof(double));
 
-/* mat = base change matrix, gram = LLL-reduced T2 matrix */
-static GEN
-chk_gen_init(FP_chk_fun *chk, GEN nf, GEN gram, GEN mat, long *ptprec)
-{
-  GEN P,bound,prev,x,B,data, M = gmael(nf,5,1);
-  long N = lg(nf[7]), n = N-1,i,prec,prec2;
-  int skipfirst = 1; /* [1,0...0] --> x rational */
+  Mbar.W = dalloc((n+1)*sizeof(double));
+  for (i=1; i< n; i++)  Pbar[i] = dalloc((n+1)*sizeof(double));
+  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));
+  for (i=1; i<=n; i++)  Hbar[i] = dalloc(n*sizeof(double));
+  for (i=1; i<=n; i++) Mbarst.A[i] = dalloc((n+1)*sizeof(double));
+  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
- * data[1] = embeddings of LLL-reduced Zk basis
- * 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;
+  gabar = gtodouble((GEN)tabga[1]); tabgabar[1] = gabar;
+  for (i=2; i<n; i++) tabgabar[i] = tabgabar[i-1]*gabar;
 
-  x = cgetg(N,t_COL);
-  bound = get_Bnf(nf); prev = NULL;
-  for (i=1; i<N; i++) x[i]=zero;
-  for (i=2; i<N; i++)
+  av = avma;
+  if (DEBUGLEVEL>=3) printf("Initialization time = %ld\n",timer());
+RESTART:
+  flit = initializedoubles(&Mbar, &M, prec);
+  storeprecdoubles(&Mbarst, &Mbar);
+  if (flit) dLQdec(&Mbar, Pbar);
+  ctpro = 0;
+  for (;;)
   {
-    x[i] = un;
-    P = get_polmin(data,x);
-    x[i] = zero;
-    if (degpol(P) == n)
+    if (low_stack(lim, stack_lim(av,1)))
     {
-      B = gcoeff(gram,i,i);
-      if (gcmp(B,bound) < 0) bound = B ;
+      if(DEBUGMEM>1) err(warnmem,"pslq");
+      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);
-      if (skipfirst == i-1)
+      ctpro++;
+      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 */
-          P = (GEN)compositum(prev,P)[1];
-          if (degpol(P) == n) continue;
-          if (DEBUGLEVEL>2 && lgef(P)>lgef(prev))
-            fprintferr("chk_gen_init: subfield %Z\n",P);
+          if (ctpro == 1) goto DOGEN;
+          storeprecdoubles(&Mbar, &Mbarst); /* restore */
+          if (! applybar(&M, &Mbar, Abargen,Bbargen)) err(bugparier,"pslqL2");
+	  if ( (p1 = checkend(&M, prec)) ) return gerepilecopy(av0, p1);
+          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
-chk_gen_post(GEN data, GEN res)
+/* x is a vector of elts of a p-adic field */
+GEN
+plindep(GEN x)
 {
-  GEN basis = (GEN)data[2];
-  GEN p1 = (GEN)res[2];
-  long i, l = lg(p1);
-  for (i=1; i<l; i++) p1[i] = lmul(basis, (GEN)p1[i]);
-  return res;
-}
+  long i, j, prec = VERYBIGINT, lx = lg(x)-1, ly, v;
+  gpmem_t av = avma;
+  GEN p = NULL, pn,p1,m,a;
 
-/* no garbage collecting, done in polredabs0 */
-static GEN
-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++)
+  if (lx < 2) return cgetg(1,t_VEC);
+  for (i=1; i<=lx; i++)
   {
-    k = v[i];
-    if (!egalii((GEN)discs[k],dmin)) break;
-    if (gpolcomp((GEN)y[k],z) < 0) { z = (GEN)y[k]; beta = (GEN)a[k]; }
+    p1 = (GEN)x[i];
+    if (typ(p1) != t_PADIC) continue;
+
+    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)
-  {
-    y=cgetg(3,t_VEC);
-    y[1]=lcopy(z);
-    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;
-}
+  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)));
 
-/* no garbage collecting, done in polredabs0 */
-static GEN
-storeallpols(GEN nf, GEN z, GEN a, GEN phimax, long flun)
-{
-  GEN p1,y,beta;
-
-  if (flun & nf_RAW)
+  ly = 2*lx - 1;
+  m = cgetg(ly+1,t_MAT);
+  for (j=1; j<=ly; j++)
   {
-    long i, c = lg(z);
-    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]=lcopy((GEN)a[i]);
-    }
+    p1 = cgetg(lx+1, t_COL); m[j] = (long)p1;
+    for (i=1; i<=lx; i++) p1[i] = zero;
   }
-  else if (phimax)
+  a = negi((GEN)x[1]);
+  for (i=1; i<lx; i++)
   {
-    long i, c = lg(z);
-    beta = new_chunk(c);
-    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]);
-    }
+    coeff(m,1+i,i) = (long)a;
+    coeff(m,1  ,i) = x[i+1];
   }
-  else y = gcopy(z);
-  return y;
+  for (i=1; i<=lx; i++) coeff(m,i,lx-1+i) = (long)pn;
+  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;
-  GEN nf,v,y,a,phimax;
-  GEN (*storepols)(GEN, GEN, GEN, GEN, long);
-  FP_chk_fun *chk;
+  long tx=typ(x), i, k;
+  gpmem_t av;
+  GEN y,p1;
 
-  chk = (FP_chk_fun*)new_chunk(sizeof(FP_chk_fun));
-  chk->f         = &chk_gen;
-  chk->f_init    = &chk_gen_init;
-  chk->f_post    = &chk_gen_post;
+  if (! is_scalar_t(tx)) err(typeer,"algdep0");
+  if (tx==t_POLMOD) { y=forcecopy((GEN)x[1]); setvarn(y,0); return y; }
+  if (gcmp0(x)) return gzero;
+  if (!n) return gun;
 
-  if ((ulong)flun >= 16) err(flagerr,"polredabs");
-  nf = initalgall0(x,nf_SMALL|nf_REGULAR,prec);
-  if (lg(nf) == 3)
-  {
-    phimax = lift_to_pol((GEN)nf[2]);
-    nf = (GEN)nf[1];
-  }
+  av=avma; p1=cgetg(n+2,t_COL);
+  p1[1] = un;
+  p1[2] = (long)x; /* n >= 1 */
+  for (i=3; i<=n+1; i++) p1[i]=lmul((GEN)p1[i-1],x);
+  if (typ(x) == t_PADIC)
+    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)]);
-    a = _vec(gsub((GEN)y[1], (GEN)x[2]));
-  }
-  else
-  {
-    for (i=1; ; i++)
+    switch(bit)
     {
-      v = fincke_pohst(nf,NULL,5000,3,prec, chk);
-      if (v) break;
-      if (i==MAXITERPOL) err(accurer,"polredabs0");
-      prec = (prec<<1)-2; nf = nfnewprec(nf,prec);
-      if (DEBUGLEVEL) err(warnprec,"polredabs0",prec);
+      case 0: p1 = pslq(p1,prec); break;
+      case -1: p1 = lindep(p1,prec); break;
+      case -2: p1 = deplin(p1); break;
+      default: p1 = lindep2(p1,bit);
     }
-    a = (GEN)v[2];
-    y = (GEN)v[1];
   }
-  nv = lg(a);
-  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)
+  if ((!bit) && (typ(p1) == t_REAL))
   {
-    y = _vec(x);
-    a = _vec(polx[varn(x)]);
+    y = gcopy(p1); return gerepileupto(av,y);
   }
-  if (varn(y[1]) != varn(x))
-  {
-    long vx = varn(x);
-    for (i=1; i<nv; i++) setvarn(y[i], vx);
-  }
-  return gerepileupto(av, storepols(nf,y,a,phimax,flun));
+  if (lg(p1) < 2)
+    err(talker,"higher degree than expected in algdep");
+  y=cgetg(n+3,t_POL);
+  y[1] = evalsigne(1) | evalvarn(0);
+  k=1; while (gcmp0((GEN)p1[k])) k++;
+  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)
-{
-  return polredabs0(x,nf_ORIG,prec);
-}
+/********************************************************************/
+/**                                                                **/
+/**                   INTEGRAL KERNEL (LLL REDUCED)                **/
+/**                                                                **/
+/********************************************************************/
 
 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++)
@@ -2520,7 +2599,8 @@ 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++)
@@ -2528,7 +2608,7 @@ 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;
 }
@@ -2537,20 +2617,15 @@ 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);
-
-  /* correct alignment for the following */
-  s = avma % sizeof(double); avma -= s;
-  if (avma<bot) err(errpile);
-
-  s = (n * sizeof(double))/sizeof(long);
-  *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);
+  s = n * sizeof(double);
+  init_dalloc();
+  *y = dalloc(s);
+  *z = dalloc(s);
+  *v = dalloc(s);
+  for (i=1; i<n; i++) (*q)[i] = dalloc(s);
 }
 
 /* Minimal vectors for the integral definite quadratic form: a.
@@ -2568,7 +2643,8 @@ 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 */
@@ -2683,7 +2759,7 @@ 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
@@ -2707,7 +2783,8 @@ 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];
@@ -2788,28 +2865,28 @@ 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 (flag & 1) return NULL on precision problems (no error).
+ * If (noer) 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 data = CHECK? CHECK->data: NULL;
+  GEN (*check)(void *,GEN) = CHECK? CHECK->f: 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);
@@ -2870,10 +2947,8 @@ 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;
@@ -2882,10 +2957,10 @@ 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)
       {
@@ -2922,17 +2997,20 @@ 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;
-        GEN per = sindexsort(norms);
+        gpmem_t av2 = avma;
+        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;
@@ -2981,143 +3059,109 @@ 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).
- * flag & 1, return NULL in case of precision problems (sqred1 or lllgram)
- *   raise an error otherwise.
- * flag & 2, return as soon as stockmax vectors found.
- * If a is a number field, use its T2 matrix
- */
-GEN
-fincke_pohst(GEN a,GEN B0,long stockmax,long flag, long PREC, FP_chk_fun *CHECK)
+/* solve q(x) = x~.a.x <= bound, a > 0.
+ * If check is non-NULL keep x only if check(x).
+ * If noer, return NULL in case of precision problems, raise an error otherwise.
+ * If a is a vector, assume a[1] is the LLL-reduced Cholesky form of q */
+static GEN
+_fincke_pohst(GEN a,GEN B0,long stockmax,long noer, long PREC, FP_chk_fun *CHECK)
 {
-  VOLATILE long pr,av=avma,i,j,n;
-  VOLATILE GEN B,nf,r,rinvtrans,u,v,res,z,vnorm,sperm,perm,uperm,gram;
-  VOLATILE GEN bound = B0;
-  void *catcherr = NULL;
-  long prec = PREC;
+  gpmem_t av = avma;
+  long i,j,l, round = 0;
+  GEN B,r,rinvtrans,u,v,res,z,vnorm,sperm,perm,uperm,gram, bound = B0;
 
-  if (DEBUGLEVEL>2) { fprintferr("entering fincke_pohst\n"); flusherr(); }
-  if (typ(a) == t_VEC) { nf = checknf(a); a = gmael(nf,5,3); } else nf = NULL;
-  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 */
+  if (DEBUGLEVEL>2) fprintferr("entering fincke_pohst\n");
+  if (typ(a) == t_VEC)
   {
-    GEN T = nf_get_T((GEN)nf[1], (GEN)nf[7]);
-    u = lllgramint(T);
-    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));
+    r = (GEN)a[1];
+    u = NULL;
   }
   else
   {
-    u = lllgramintern(a,4,flag&1, (prec<<1)-2);
-    if (!u) goto PRECPB;
+    long prec = PREC;
+    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);
-  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 */
+  /* now r~ * r = a in 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;
-  for (i=1; i<6; i++) /* try to get close to a genuine LLL basis */
-  {
-    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... */
+  v = ZM_inv(gtrans_i(v),gun);
+  r = gmul(r,v);
+  u = u? gmul(u,v): v;
 
-  if (v)
-  {
-    GEN u2 = ZM_inv(gtrans(v),gun);
-    r = gmul(r,u2); /* r = LLL basis now */
-    u = gmul(u,u2);
-  }
+  l = lg(r);
+  vnorm = cgetg(l,t_VEC);
+  for (j=1; j<l; j++) vnorm[j] = lnorml2((GEN)rinvtrans[j]);
+  sperm = cgetg(l,t_MAT);
+  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);
-  if (gexpo(B) >= bit_accuracy(lg(B)-2)) goto PRECPB;
+  B = gcoeff(gram,l-1,l-1);
+  if (gexpo(B) >= bit_accuracy(lg(B)-2)) return NULL;
 
-  { /* catch precision problems (truncation error) */
-    jmp_buf env;
-    if (setjmp(env)) goto PRECPB;
-    catcherr = err_catch(precer, env, NULL);
-  }
-  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++)
+  res = NULL;
+  CATCH(precer) { }
+  TRY {
+    if (CHECK && CHECK->f_init)
     {
-      x[i] = un; B = gcoeff(gram,i,i);
-      if (!bound || mpcmp(B,bound) < 0) bound = B;
-      x[i] = zero;
+      bound = CHECK->f_init(CHECK,gram,uperm);
+      if (!bound) err(precer,"fincke_pohst");
     }
-  }
+    if (!bound) bound = gcoeff(gram,1,1);
 
-  if (DEBUGLEVEL>2) {fprintferr("entering smallvectors\n"); flusherr();}
-  for (i=1; i<n; i++)
+    if (DEBUGLEVEL>2) fprintferr("entering smallvectors\n");
+    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),
-                       stockmax,flag,CHECK);
-    if (!res) goto PRECPB;
-    if (!CHECK || bound || lg(res[2]) > 1) break;
-    if (DEBUGLEVEL>2) fprintferr("  i = %ld failed\n",i);
+    if (!(flag & 1)) err(precer,"fincke_pohst");
+    avma = av;
   }
-  err_leave(&catcherr); catcherr = NULL;
-  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;
+  return z;
 }