===================================================================
RCS file: /home/cvs/OpenXM_contrib/pari-2.2/src/basemath/Attic/alglin2.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/alglin2.c	2001/10/02 11:17:00	1.1.1.1
+++ OpenXM_contrib/pari-2.2/src/basemath/Attic/alglin2.c	2002/09/11 07:26:48	1.2
@@ -1,4 +1,4 @@
-/* $Id: alglin2.c,v 1.1.1.1 2001/10/02 11:17:00 noro Exp $
+/* $Id: alglin2.c,v 1.2 2002/09/11 07:26:48 noro Exp $
 
 Copyright (C) 2000  The PARI group.
 
@@ -43,7 +43,8 @@ charpoly0(GEN x, int v, long flag)
 static GEN
 caract2_i(GEN p, GEN x, int v, GEN (subres_f)(GEN,GEN,GEN*))
 {
-  long av = avma, d;
+  gpmem_t av = avma;
+  long d;
   GEN p1, p2 = leading_term(p);
 
   if (!signe(x)) return gpowgs(polx[v], degpol(p));
@@ -55,7 +56,7 @@ caract2_i(GEN p, GEN x, int v, GEN (subres_f)(GEN,GEN,
   else
     p1 = gsubst(p1,MAXVARN,polx[v]);
 
-  if (!gcmp1(p2) && (d=degpol(x)) > 0) p1 = gdiv(p1, gpuigs(p2,d));
+  if (!gcmp1(p2) && (d=degpol(x)) > 0) p1 = gdiv(p1, gpowgs(p2,d));
   return gerepileupto(av,p1);
 }
 
@@ -75,7 +76,8 @@ caractducos(GEN p, GEN x, int v)
 static GEN
 easychar(GEN x, int v, GEN *py)
 {
-  long l,tetpil,lx;
+  gpmem_t av;
+  long lx;
   GEN p1,p2;
 
   switch(typ(x))
@@ -96,11 +98,11 @@ easychar(GEN x, int v, GEN *py)
 
     case t_COMPLEX: case t_QUAD:
       if (py) err(typeer,"easychar");
-      p1=cgetg(5,t_POL);
-      p1[1]=evalsigne(1) | evallgef(5) | evalvarn(v);
-      p1[2]=lnorm(x); l=avma; p2=gtrace(x); tetpil=avma;
-      p1[3]=lpile(l,tetpil,gneg(p2));
-      p1[4]=un; return p1;
+      p1 = cgetg(5,t_POL);
+      p1[1] = evalsigne(1) | evallgef(5) | evalvarn(v);
+      p1[2] = lnorm(x); av = avma;
+      p1[3] = lpileupto(av, gneg(gtrace(x)));
+      p1[4] = un; return p1;
 
     case t_POLMOD:
       if (py) err(typeer,"easychar");
@@ -123,24 +125,26 @@ easychar(GEN x, int v, GEN *py)
 GEN
 caract(GEN x, int v)
 {
-  long n,k,l=avma,tetpil;
-  GEN  p1,p2,p3,p4,p5,p6;
+  long k, n;
+  gpmem_t av=avma;
+  GEN  p1,p2,p3,p4,x_k;
 
   if ((p1 = easychar(x,v,NULL))) return p1;
 
-  p1=gzero; p2=gun;
-  n=lg(x)-1; if (n&1) p2=gneg_i(p2);
-  p4=cgetg(3,t_RFRACN); p5=dummycopy(polx[v]); p4[2]=(long)p5;
-  p6=cgeti(3); p6[1] = evalsigne(-1) | evallgefint(3);
+  p1 = gzero; p2 = gun;
+  n = lg(x)-1; if (n&1) p2 = negi(p2);
+  x_k = dummycopy(polx[v]);
+  p4 = cgetg(3,t_RFRACN); p4[2] = (long)x_k;
   for (k=0; k<=n; k++)
   {
-    p3=det(gsub(gscalmat(stoi(k),n), x));
-    p4[1]=lmul(p3,p2); p6[2]=k;
-    p1=gadd(p4,p1); p5[2]=(long)p6;
-    if (k!=n) p2=gdivgs(gmulsg(k-n,p2),k+1);
+    p3 = det(gsub(gscalmat(stoi(k),n), x));
+    p4[1] = lmul(p3,p2); x_k[2] = lstoi(-k);
+    p1 = gadd(p4,p1);
+    if (k == n) break;
+
+    p2 = gdivgs(gmulsg(k-n,p2),k+1);
   }
-  p2=mpfact(n); tetpil=avma;
-  return gerepile(l,tetpil,gdiv((GEN) p1[1],p2));
+  return gerepileupto(av, gdiv((GEN)p1[1], mpfact(n)));
 }
 
 GEN
@@ -155,7 +159,8 @@ caradj0(GEN x, long v)
 GEN
 caradj(GEN x, long v, GEN *py)
 {
-  long i,j,k,l,av,tetpil;
+  gpmem_t av,tetpil;
+  long i,j,k,l;
   GEN p,y,t,*gptr[2];
 
   if ((p = easychar(x,v,py))) return p;
@@ -217,7 +222,7 @@ GEN
 adj(GEN x)
 {
   GEN y;
-  caradj(x,MAXVARN,&y); return y;
+  (void)caradj(x,MAXVARN,&y); return y;
 }
 
 /*******************************************************************/
@@ -225,13 +230,13 @@ adj(GEN x)
 /*                       HESSENBERG FORM                           */
 /*                                                                 */
 /*******************************************************************/
-#define swap(x,y) { long _t=x; x=y; y=_t; }
-#define gswap(x,y) { GEN _t=x; x=y; y=_t; }
+#define lswap(x,y) { long _t=x; x=y; y=_t; }
+#define swap(x,y) { GEN _t=x; x=y; y=_t; }
 
 GEN
 hess(GEN x)
 {
-  ulong av = avma;
+  gpmem_t av = avma;
   long lx=lg(x),m,i,j;
   GEN p,p1,p2;
 
@@ -246,8 +251,8 @@ hess(GEN x)
       p = gcoeff(x,i,m-1);
       if (gcmp0(p)) continue;
 
-      for (j=m-1; j<lx; j++) swap(coeff(x,i,j), coeff(x,m,j));
-      swap(x[i], x[m]); p = ginv(p);
+      for (j=m-1; j<lx; j++) lswap(coeff(x,i,j), coeff(x,m,j));
+      lswap(x[i], x[m]); p = ginv(p);
       for (i=m+1; i<lx; i++)
       {
         p1 = gcoeff(x,i,m-1);
@@ -268,7 +273,8 @@ hess(GEN x)
 GEN
 carhess(GEN x, long v)
 {
-  long av,tetpil,lx,r,i;
+  gpmem_t av,tetpil;
+  long lx,r,i;
   GEN *y,p1,p2,p3,p4;
 
   if ((p1 = easychar(x,v,NULL))) return p1;
@@ -299,7 +305,8 @@ carhess(GEN x, long v)
 GEN
 gnorm(GEN x)
 {
-  long l,lx,i,tetpil, tx=typ(x);
+  gpmem_t av;
+  long lx,i, tx=typ(x);
   GEN p1,p2,y;
 
   switch(tx)
@@ -313,27 +320,30 @@ gnorm(GEN x)
     case t_FRAC: case t_FRACN:
       return gsqr(x);
 
-    case t_COMPLEX:
-      l=avma; p1=gsqr((GEN) x[1]); p2=gsqr((GEN) x[2]);
-      tetpil=avma; return gerepile(l,tetpil,gadd(p1,p2));
+    case t_COMPLEX: av = avma;
+      return gerepileupto(av, gadd(gsqr((GEN)x[1]), gsqr((GEN)x[2])));
 
-    case t_QUAD:
-      l=avma; p1=(GEN)x[1];
-      p2=gmul((GEN) p1[2], gsqr((GEN) x[3]));
-      p1 = gcmp0((GEN) p1[3])? gsqr((GEN) x[2])
-                             : gmul((GEN) x[2], gadd((GEN) x[2],(GEN) x[3]));
-      tetpil=avma; return gerepile(l,tetpil,gadd(p1,p2));
+    case t_QUAD: av = avma;
+      p1 = (GEN)x[1];
+      p2 = gmul((GEN)p1[2], gsqr((GEN)x[3]));
+      p1 = gcmp0((GEN)p1[3])? gsqr((GEN)x[2])
+                            : gmul((GEN)x[2], gadd((GEN)x[2],(GEN)x[3]));
+      return gerepileupto(av, gadd(p1,p2));
 
-    case t_POL: case t_SER: case t_RFRAC: case t_RFRACN:
-      l=avma; p1=gmul(gconj(x),x); tetpil=avma;
-      return gerepile(l,tetpil,greal(p1));
+    case t_POL: case t_SER: case t_RFRAC: case t_RFRACN: av = avma;
+      return gerepileupto(av, greal(gmul(gconj(x),x)));
 
     case t_POLMOD:
-      p1=(GEN)x[1]; p2=leading_term(p1);
-      if (gcmp1(p2) || gcmp0((GEN) x[2])) return subres(p1,(GEN) x[2]);
-      l=avma; y=subres(p1,(GEN)x[2]); p1=gpuigs(p2,degpol(x[2]));
-      tetpil=avma; return gerepile(l,tetpil,gdiv(y,p1));
+    {
+      GEN T = (GEN)x[1], A = (GEN)x[2];
+      if (typ(A) != t_POL) return gpowgs(A, degpol(T));
 
+      y = subres(T, A); p1 = leading_term(T);
+      if (gcmp1(p1) || gcmp0(A)) return y;
+      av = avma; T = gpowgs(p1, degpol(A));
+      return gerepileupto(av, gdiv(y,T));
+    }
+
     case t_VEC: case t_COL: case t_MAT:
       lx=lg(x); y=cgetg(lx,tx);
       for (i=1; i<lx; i++) y[i]=lnorm((GEN) x[i]);
@@ -346,8 +356,9 @@ gnorm(GEN x)
 GEN
 gnorml2(GEN x)
 {
+  gpmem_t av,lim;
   GEN y;
-  long i,tx=typ(x),lx,av,lim;
+  long i,tx=typ(x),lx;
 
   if (! is_matvec_t(tx)) return gnorm(x);
   lx=lg(x); if (lx==1) return gzero;
@@ -368,7 +379,7 @@ gnorml2(GEN x)
 GEN
 QuickNormL2(GEN x, long prec)
 {
-  long av = avma;
+  gpmem_t av = avma;
   GEN y = gmul(x, realun(prec));
   if (typ(x) == t_POL)
     *++y = evaltyp(t_VEC) | evallg(lgef(x)-1);
@@ -378,7 +389,7 @@ QuickNormL2(GEN x, long prec)
 GEN
 gnorml1(GEN x,long prec)
 {
-  ulong av = avma;
+  gpmem_t av = avma;
   long lx,i;
   GEN s;
   switch(typ(x))
@@ -406,7 +417,7 @@ gnorml1(GEN x,long prec)
 GEN
 QuickNormL1(GEN x,long prec)
 {
-  ulong av = avma;
+  gpmem_t av = avma;
   long lx,i;
   GEN p1,p2,s;
   switch(typ(x))
@@ -497,7 +508,8 @@ gconj(GEN x)
 GEN
 conjvec(GEN x,long prec)
 {
-  long lx,s,av,tetpil,i,tx=typ(x);
+  gpmem_t av,tetpil;
+  long lx,s,i,tx=typ(x);
   GEN z,y,p1,p2,p;
 
   switch(tx)
@@ -511,9 +523,10 @@ conjvec(GEN x,long prec)
     case t_VEC: case t_COL:
       lx=lg(x); z=cgetg(lx,t_MAT);
       for (i=1; i<lx; i++) z[i]=(long)conjvec((GEN)x[i],prec);
-      s=lg(z[1]);
+      if (lx == 1) break;
+      s = lg(z[1]);
       for (i=2; i<lx; i++)
-	if (lg(z[i])!=s) err(talker,"incompatible field degrees in conjvec");
+        if (lg(z[i])!=s) err(talker,"incompatible field degrees in conjvec");
       break;
 
     case t_POLMOD:
@@ -573,15 +586,17 @@ assmat(GEN x)
       p1[j] = (j==(i+1))? un: zero;
   }
   p1=cgetg(lx,t_COL); y[i]=(long)p1;
-  if (gcmp1((GEN) x[lx+1]))
+  if (gcmp1((GEN)x[lx+1]))
     for (j=1; j<lx; j++)
       p1[j] = lneg((GEN)x[j+1]);
   else
   {
-    p2 = (GEN)x[lx+1]; gnegz(p2,p2);
+    gpmem_t av = avma;
+    p2 = gclone(gneg((GEN)x[lx+1]));
+    avma = av;
     for (j=1; j<lx; j++)
       p1[j] = ldiv((GEN)x[j+1],p2);
-    gnegz(p2,p2);
+    gunclone(p2);
   }
   return y;
 }
@@ -589,7 +604,8 @@ assmat(GEN x)
 GEN
 gtrace(GEN x)
 {
-  long i,l,n,tx=typ(x),lx,tetpil;
+  gpmem_t av, tetpil;
+  long i,n,tx=typ(x),lx;
   GEN y,p1,p2;
 
   switch(tx)
@@ -604,8 +620,8 @@ gtrace(GEN x)
       p1=(GEN)x[1];
       if (!gcmp0((GEN) p1[3]))
       {
-	l=avma; p2=gmul2n((GEN)x[2],1);
-	tetpil=avma; return gerepile(l,tetpil,gadd((GEN)x[3],p2));
+	av=avma; p2=gmul2n((GEN)x[2],1);
+	tetpil=avma; return gerepile(av,tetpil,gadd((GEN)x[3],p2));
       }
       return gmul2n((GEN)x[2],1);
 
@@ -620,11 +636,11 @@ gtrace(GEN x)
       return y;
 
     case t_POLMOD:
-      l=avma; n=(lgef(x[1])-4);
+      av=avma; n=(lgef(x[1])-4);
       p1=polsym((GEN)x[1],n); p2=gzero;
       for (i=0; i<=n; i++)
         p2=gadd(p2,gmul(truecoeff((GEN)x[2],i),(GEN)p1[i+1]));
-      return gerepileupto(l,p2);
+      return gerepileupto(av,p2);
 
     case t_RFRAC: case t_RFRACN:
       return gadd(x,gconj(x));
@@ -637,27 +653,26 @@ gtrace(GEN x)
     case t_MAT:
       lx=lg(x); if (lx==1) return gzero;/*now lx>=2*/
       if (lx!=lg(x[1])) err(mattype1,"gtrace");
-      l=avma; p1=gcoeff(x,1,1); if (lx==2) return gcopy(p1);
+      av=avma; p1=gcoeff(x,1,1); if (lx==2) return gcopy(p1);
       for (i=2; i<lx-1; i++)
 	p1=gadd(p1,gcoeff(x,i,i));
-      tetpil=avma; return gerepile(l,tetpil,gadd(p1,gcoeff(x,i,i)));
+      tetpil=avma; return gerepile(av,tetpil,gadd(p1,gcoeff(x,i,i)));
 
   }
   err(typeer,"gtrace");
   return NULL; /* not reached */
 }
 
-/* Gauss reduction for positive definite matrix a
+/* Cholesky Decomposition for positive definite matrix a
  * If a is not positive definite:
  *   if flag is zero: raise an error
- *   else: return NULL.
- */
+ *   else: return NULL. */
 GEN
 sqred1intern(GEN a,long flag)
 {
+  gpmem_t av = avma, lim=stack_lim(av,1);
   GEN b,p;
   long i,j,k, n = lg(a);
-  ulong av = avma, lim=stack_lim(av,1);
 
   if (typ(a)!=t_MAT) err(typeer,"sqred1");
   if (n == 1) return cgetg(1, t_MAT);
@@ -704,7 +719,7 @@ sqred1(GEN a)
 GEN
 sqred3(GEN a)
 {
-  ulong av = avma, lim = stack_lim(av,1);
+  gpmem_t av = avma, lim = stack_lim(av,1);
   long i,j,k,l, n = lg(a);
   GEN p1,b;
 
@@ -746,16 +761,17 @@ static GEN
 sqred2(GEN a, long no_signature)
 {
   GEN r,p,mun;
-  ulong av,av1,lim;
+  gpmem_t av,av1,lim;
   long n,i,j,k,l,sp,sn,t;
 
   if (typ(a)!=t_MAT) err(typeer,"sqred2");
-  n=lg(a); if (lg(a[1]) != n) err(mattype1,"sqred2");
+  n = lg(a); if (n > 1 && lg(a[1]) != n) err(mattype1,"sqred2");
 
-  av=avma; mun=negi(gun); r=new_chunk(n);
-  for (i=1; i<n; i++) r[i]=1;
-  av1=avma; lim=stack_lim(av1,1); a=dummycopy(a);
-  n--; t=n; sp=sn=0;
+  av = avma; mun = negi(gun);
+  r = vecsmall_const(n-1, 1);
+  av1= avma; lim = stack_lim(av1,1);
+  a = dummycopy(a);
+  n--; t = n; sp = sn = 0;
   while (t)
   {
     k=1; while (k<=n && (gcmp0(gcoeff(a,k,k)) || !r[k])) k++;
@@ -806,9 +822,10 @@ sqred2(GEN a, long no_signature)
       if (k>n) break;
     }
   }
-  if (no_signature) return gerepilecopy(av,a);
-  avma=av;
-  a=cgetg(3,t_VEC); a[1]=lstoi(sp); a[2]=lstoi(sn); return a;
+  if (no_signature) return gerepilecopy(av, a);
+  avma = av; a = cgetg(3,t_VEC);
+  a[1] = lstoi(sp);
+  a[2] = lstoi(sn); return a;
 }
 
 GEN
@@ -824,7 +841,8 @@ signat(GEN a) { return sqred2(a,0); }
 GEN
 jacobi(GEN a, long prec)
 {
-  long de,e,e1,e2,l,n,i,j,p,q,av1,av2;
+  gpmem_t av1, av2;
+  long de,e,e1,e2,l,n,i,j,p,q;
   GEN c,s,t,u,ja,lambda,r,unr,x,y,x1,y1;
 
   if (typ(a)!=t_MAT) err(mattype1,"jacobi");
@@ -840,9 +858,8 @@ jacobi(GEN a, long prec)
   r=cgetg(l,t_MAT); ja[2]=(long)r;
   for (j=1; j<=n; j++)
   {
-    r[j]=lgetg(l,t_COL);
-    for (i=1; i<=n; i++)
-      affsr(i==j, (GEN)(coeff(r,i,j)=lgetr(prec)));
+    r[j] = lgetg(l,t_COL);
+    for (i=1; i<=n; i++) coeff(r,i,j) = (long)stor(i==j, prec);
   }
   av1=avma;
 
@@ -940,63 +957,63 @@ matrixqz0(GEN x,GEN p)
   err(flagerr,"matrixqz"); return NULL; /* not reached */
 }
 
+static int
+ZV_isin(GEN x)
+{
+  long i,l = lg(x);
+  for (i=1; i<l; i++)
+    if (typ(x[i]) != t_INT) return 0;
+  return 1;
+}
+
 GEN
 matrixqz(GEN x, GEN p)
 {
-  ulong av = avma, av1, lim;
-  long i,j,j1,m,n,t,nfact;
-  GEN p1,p2,p3,unmodp;
+  gpmem_t av = avma, av1, lim;
+  long i,j,j1,m,n,nfact;
+  GEN p1,p2,p3;
 
-  if (typ(x)!=t_MAT) err(typeer,"matrixqz");
-  n=lg(x)-1; if (!n) return gcopy(x);
-  m=lg(x[1])-1;
+  if (typ(x) != t_MAT) err(typeer,"matrixqz");
+  n = lg(x)-1; if (!n) return gcopy(x);
+  m = lg(x[1])-1;
   if (n > m) err(talker,"more rows than columns in matrixqz");
   if (n==m)
   {
-    p1=det(x);
+    p1 = det(x);
     if (gcmp0(p1)) err(talker,"matrix of non-maximal rank in matrixqz");
-    avma=av; return idmat(n);
+    avma = av; return idmat(n);
   }
-  p1=cgetg(n+1,t_MAT);
+  /* m > n */
+  p1 = x; x = cgetg(n+1,t_MAT);
   for (j=1; j<=n; j++)
   {
-    p2=gun; p3=(GEN)x[j];
-    for (i=1; i<=m; i++)
-    {
-      t=typ(p3[i]);
-      if (t != t_INT && !is_frac_t(t))
-        err(talker,"not a rational or integral matrix in matrixqz");
-      p2=ggcd(p2,(GEN)p3[i]);
-    }
-    p1[j]=ldiv(p3,p2);
+    x[j] = (long)primpart((GEN)p1[j]);
+    if (!ZV_isin((GEN)p1[j])) err(talker, "not a rational matrix in matrixqz");
   }
-  x=p1; unmodp=cgetg(3,t_INTMOD); unmodp[2]=un;
+  /* x integral */
 
   if (gcmp0(p))
   {
-    p1=cgetg(n+1,t_MAT); p2=gtrans(x);
-    for (j=1; j<=n; j++) p1[j]=p2[j];
-    p3=det(p1); p1[n]=p2[n+1]; p3=ggcd(p3,det(p1));
-    if (!signe(p3))
+    p1 = gtrans(x); setlg(p1,n+1);
+    p2 = det(p1); p1[n] = p1[n+1]; p2 = ggcd(p2,det(p1));
+    if (!signe(p2))
       err(impl,"matrixqz when the first 2 dets are zero");
-    if (gcmp1(p3)) return gerepilecopy(av,x);
+    if (gcmp1(p2)) return gerepilecopy(av,x);
 
-    p3=factor(p3); p1=(GEN)p3[1]; nfact=lg(p1)-1;
+    p1 = (GEN)factor(p2)[1];
   }
-  else
-  {
-    p1=cgetg(2,t_VEC); p1[1]=(long)p; nfact=1;
-  }
-  av1=avma; lim=stack_lim(av1,1);
+  else p1 = _vec(p);
+  nfact = lg(p1)-1;
+  av1 = avma; lim = stack_lim(av1,1);
   for (i=1; i<=nfact; i++)
   {
-    p=(GEN)p1[i]; unmodp[1]=(long)p;
+    p = (GEN)p1[i];
     for(;;)
     {
-      p2=ker(gmul(unmodp,x));
+      p2 = FpM_ker(x, p);
       if (lg(p2)==1) break;
 
-      p2=centerlift(p2); p3=gdiv(gmul(x,p2),p);
+      p2 = centermod(p2,p); p3 = gdiv(gmul(x,p2), p);
       for (j=1; j<lg(p2); j++)
       {
 	j1=n; while (gcmp0(gcoeff(p2,j1,j))) j1--;
@@ -1005,7 +1022,7 @@ matrixqz(GEN x, GEN p)
       if (low_stack(lim, stack_lim(av1,1)))
       {
 	if (DEBUGMEM>1) err(warnmem,"matrixqz");
-	x=gerepilecopy(av1,x);
+	x = gerepilecopy(av1,x);
       }
     }
   }
@@ -1013,90 +1030,138 @@ matrixqz(GEN x, GEN p)
 }
 
 static GEN
-matrixqz_aux(GEN x, long m, long n)
+Z_V_mul(GEN u, GEN A)
 {
-  ulong av = avma, lim = stack_lim(av,1);
-  long i,j,k,in[2];
-  GEN p1;
+  if (gcmp1(u)) return A;
+  if (gcmp_1(u)) return gneg(A);
+  if (gcmp0(u)) return zerocol(lg(A)-1);
+  return gmul(u,A);
+}
 
-  for (i=1; i<=m; i++)
+static GEN
+QV_lincomb(GEN u, GEN v, GEN A, GEN B)
+{
+  if (!signe(u)) return Z_V_mul(v,B);
+  if (!signe(v)) return Z_V_mul(u,A);
+  return gadd(Z_V_mul(u,A), Z_V_mul(v,B));
+}
+
+/* cf ZV_elem */
+/* zero aj = Aij (!= 0)  using  ak = Aik (maybe 0), via linear combination of
+ * A[j] and A[k] of determinant 1. */
+static void
+QV_elem(GEN aj, GEN ak, GEN A, long j, long k)
+{
+  GEN p1,u,v,d, D;
+
+  if (gcmp0(ak)) { lswap(A[j],A[k]); return; }
+  D = mpppcm(denom(aj), denom(ak));
+  if (!is_pm1(D)) { aj = gmul(aj,D); ak = gmul(ak,D); }
+  d = bezout(aj,ak,&u,&v);
+  /* frequent special case (u,v) = (1,0) or (0,1) */
+  if (!signe(u))
+  { /* ak | aj */
+    p1 = negi(divii(aj,ak));
+    A[j]   = (long)QV_lincomb(gun, p1, (GEN)A[j], (GEN)A[k]);
+    return;
+  }
+  if (!signe(v))
+  { /* aj | ak */
+    p1 = negi(divii(ak,aj));
+    A[k]   = (long)QV_lincomb(gun, p1, (GEN)A[k], (GEN)A[j]);
+    lswap(A[j], A[k]);
+    return;
+  }
+
+  if (!is_pm1(d)) { aj = divii(aj,d); ak = divii(ak,d); }
+  p1 = (GEN)A[k]; aj = negi(aj);
+  A[k] = (long)QV_lincomb(u,v, (GEN)A[j],p1);
+  A[j] = (long)QV_lincomb(aj,ak, p1,(GEN)A[j]);
+}
+
+static GEN
+matrixqz_aux(GEN A)
+{
+  gpmem_t av = avma, lim = stack_lim(av,1);
+  long i,j,k,n,m;
+  GEN a;
+
+  n = lg(A); if (n == 1) return cgetg(1,t_MAT);
+  m = lg(A[1]);
+  for (i=1; i<m; i++)
   {
-    for(;;)
+    for (j=k=1; j<n; j++)
     {
-      long fl=0;
+      GEN a = gcoeff(A,i,j);
+      if (gcmp0(a)) continue;
 
-      for (j=1; j<=n; j++)
-	if (!gcmp0(gcoeff(x,i,j)))
-	  { in[fl]=j; fl++; if (fl==2) break; }
-      if (j>n) break;
-
-      j=(gcmp(gabs(gcoeff(x,i,in[0]),DEFAULTPREC),
-	      gabs(gcoeff(x,i,in[1]),DEFAULTPREC)) > 0)? in[1]: in[0];
-      p1=gcoeff(x,i,j);
-      for (k=1; k<=n; k++)
-	if (k!=j)
-	  x[k]=lsub((GEN)x[k],gmul(ground(gdiv(gcoeff(x,i,k),p1)),(GEN)x[j]));
+      k = j+1; if (k == n) k = 1;
+      /* zero a = Aij  using  b = Aik */
+      QV_elem(a, gcoeff(A,i,k), A, j,k);
     }
-
-    j=1; while (j<=n && gcmp0(gcoeff(x,i,j))) j++;
-    if (j<=n)
+    a = gcoeff(A,i,k);
+    if (!gcmp0(a))
     {
-      p1=denom(gcoeff(x,i,j));
-      if (!gcmp1(p1)) x[j]=lmul(p1,(GEN)x[j]);
+      a = denom(a);
+      if (!is_pm1(a)) A[k] = lmul((GEN)A[k], a);
     }
     if (low_stack(lim, stack_lim(av,1)))
     {
       if(DEBUGMEM>1) err(warnmem,"matrixqz_aux");
-      x=gerepilecopy(av,x);
+      A = gerepilecopy(av,A);
     }
   }
-  return hnf(x);
+  return m > 100? hnfall_i(A,NULL,1): hnf(A);
 }
 
 GEN
 matrixqz2(GEN x)
 {
-  long av = avma,m,n;
-
+  gpmem_t av = avma;
   if (typ(x)!=t_MAT) err(typeer,"matrixqz2");
-  n=lg(x)-1; if (!n) return gcopy(x);
-  m=lg(x[1])-1; x=dummycopy(x);
-  return gerepileupto(av, matrixqz_aux(x,m,n));
+  x = dummycopy(x);
+  return gerepileupto(av, matrixqz_aux(x));
 }
 
 GEN
 matrixqz3(GEN x)
 {
-  long av=avma,av1,j,j1,k,m,n,lim;
+  gpmem_t av = avma, av1, lim;
+  long j,j1,k,m,n;
   GEN c;
 
   if (typ(x)!=t_MAT) err(typeer,"matrixqz3");
-  n=lg(x)-1; if (!n) return gcopy(x);
-  m=lg(x[1])-1; x=dummycopy(x); c=new_chunk(n+1);
-  for (j=1; j<=n; j++) c[j]=0;
-  av1=avma; lim=stack_lim(av1,1);
-  for (k=1; k<=m; k++)
+  n = lg(x); if (n==1) return gcopy(x);
+  m = lg(x[1]); x = dummycopy(x);
+  c = cgetg(n,t_VECSMALL);
+  for (j=1; j<n; j++) c[j] = 0;
+  av1 = avma; lim = stack_lim(av1,1);
+  for (k=1; k<m; k++)
   {
-    j=1; while (j<=n && (c[j] || gcmp0(gcoeff(x,k,j)))) j++;
-    if (j<=n)
-    {
-      c[j]=k; x[j]=ldiv((GEN)x[j],gcoeff(x,k,j));
-      for (j1=1; j1<=n; j1++)
-	if (j1!=j) x[j1]=lsub((GEN)x[j1],gmul(gcoeff(x,k,j1),(GEN)x[j]));
-      if (low_stack(lim, stack_lim(av1,1)))
+    j=1; while (j<n && (c[j] || gcmp0(gcoeff(x,k,j)))) j++;
+    if (j==n) continue;
+   
+    c[j]=k; x[j]=ldiv((GEN)x[j],gcoeff(x,k,j));
+    for (j1=1; j1<n; j1++)
+      if (j1!=j)
       {
-	if(DEBUGMEM>1) err(warnmem,"matrixqz3");
-	x=gerepilecopy(av1,x);
+        GEN t = gcoeff(x,k,j1);
+        if (!gcmp0(t)) x[j1] = lsub((GEN)x[j1],gmul(t,(GEN)x[j]));
       }
+    if (low_stack(lim, stack_lim(av1,1)))
+    {
+      if(DEBUGMEM>1) err(warnmem,"matrixqz3");
+      x = gerepilecopy(av1,x);
     }
   }
-  return gerepileupto(av, matrixqz_aux(x,m,n));
+  return gerepileupto(av, matrixqz_aux(x));
 }
 
 GEN
 intersect(GEN x, GEN y)
 {
-  long av,tetpil,j, lx = lg(x);
+  gpmem_t av,tetpil;
+  long j, lx = lg(x);
   GEN z;
 
   if (typ(x)!=t_MAT || typ(y)!=t_MAT) err(typeer,"intersect");
@@ -1128,7 +1193,7 @@ mathnf0(GEN x, long flag)
 }
 
 static GEN
-init_hnf(GEN x, GEN *denx, long *co, long *li, long *av)
+init_hnf(GEN x, GEN *denx, long *co, long *li, gpmem_t *av)
 {
   if (typ(x) != t_MAT) err(typeer,"mathnf");
   *co=lg(x); if (*co==1) return NULL; /* empty matrix */
@@ -1136,7 +1201,7 @@ init_hnf(GEN x, GEN *denx, long *co, long *li, long *a
 
   if (gcmp1(*denx)) /* no denominator */
     { *denx = NULL; return dummycopy(x); }
-  return gmul(x,*denx);
+  return Q_muli_to_int(x, *denx);
 }
 
 GEN
@@ -1182,13 +1247,14 @@ ZV_Z_mul(GEN c, GEN X)
 GEN
 ZV_lincomb(GEN u, GEN v, GEN X, GEN Y)
 {
-  long av,i,lx,m;
+  gpmem_t av;
+  long i,lx,m;
   GEN p1,p2,A;
 
   if (!signe(u)) return ZV_Z_mul(v,Y);
   if (!signe(v)) return ZV_Z_mul(u,X);
   lx = lg(X); A = cgetg(lx,t_COL); m = lgefint(u)+lgefint(v)+4;
-  if (is_pm1(v)) { gswap(u,v); gswap(X,Y); }
+  if (is_pm1(v)) { swap(u,v); swap(X,Y); }
   if (is_pm1(u))
   {
     if (signe(u) > 0)
@@ -1243,7 +1309,8 @@ ZV_lincomb(GEN u, GEN v, GEN X, GEN Y)
 GEN
 hnf_special(GEN x, long remove)
 {
-  long av0, s,li,co,av,tetpil,i,j,k,def,ldef,lim;
+  gpmem_t av0,av,tetpil,lim;
+  long s,li,co,i,j,k,def,ldef;
   GEN p1,u,v,d,denx,a,b, x2,res;
 
   if (typ(x) != t_VEC || lg(x) !=3) err(typeer,"hnf_special");
@@ -1311,7 +1378,7 @@ hnf_special(GEN x, long remove)
     for (i=1,j=1; j<co; j++)
       if (!gcmp0((GEN)x[j])) 
       {
-        x[i] = x[j];
+        x[i]  = x[j];
         x2[i] = x2[j]; i++;
       }
     setlg(x,i);
@@ -1367,7 +1434,8 @@ hnffinal(GEN matgen,GEN perm,GEN* ptdep,GEN* ptB,GEN* 
 {
   GEN p1,p2,U,H,Hnew,Bnew,Cnew,diagH1;
   GEN B = *ptB, C = *ptC, dep = *ptdep, depnew;
-  long av,i,j,k,s,i1,j1,lim,zc;
+  gpmem_t av, lim;
+  long i,j,k,s,i1,j1,zc;
   long co = lg(C);
   long col = lg(matgen)-1;
   long lnz, nlze, lig;
@@ -1387,7 +1455,8 @@ hnffinal(GEN matgen,GEN perm,GEN* ptdep,GEN* ptB,GEN* 
     }
   }
   p1 = hnflll(matgen);
-  H = (GEN)p1[1]; /* lnz x lnz */
+  H = (GEN)p1[1]; /* lnz x lnz [disregarding initial 0 cols] */
+  H += (lg(H)-1 - lnz); H[0] = evaltyp(t_MAT) | evallg(lnz+1);
   U = (GEN)p1[2]; /* col x col */
   /* Only keep the part above the H (above the 0s is 0 since the dep rows
    * are dependant from the ones in matgen) */
@@ -1397,7 +1466,7 @@ hnffinal(GEN matgen,GEN perm,GEN* ptdep,GEN* ptB,GEN* 
   diagH1 = new_chunk(lnz+1); /* diagH1[i] = 0 iff H[i,i] != 1 (set later) */
 
   av = avma; lim = stack_lim(av,1);
-  Cnew = cgetg(co,t_MAT);
+  Cnew = cgetg(co, typ(C));
   setlg(C, col+1); p1 = gmul(C,U);
   for (j=1; j<=col; j++) Cnew[j] = p1[j];
   for (   ; j<co ; j++)  Cnew[j] = C[j];
@@ -1406,23 +1475,23 @@ hnffinal(GEN matgen,GEN perm,GEN* ptdep,GEN* ptB,GEN* 
   /* Clean up B using new H */
   for (s=0,i=lnz; i; i--)
   {
-    GEN h = gcoeff(H,i,i);
+    GEN Di = (GEN)dep[i], Hi = (GEN)H[i];
+    GEN h = (GEN)Hi[i]; /* H[i,i] */
     if ( (diagH1[i] = gcmp1(h)) ) { h = NULL; s++; }
     for (j=col+1; j<co; j++)
     {
       GEN z = (GEN)B[j-col];
-      p1 = (GEN)z[i+nlze]; if (h) p1 = gdivent(p1,h);
-      for (k=1; k<=nlze; k++)
-	z[k] = lsubii((GEN)z[k], mulii(p1, gcoeff(dep,k,i)));
-      for (   ; k<=lig; k++)
-	z[k] = lsubii((GEN)z[k], mulii(p1, gcoeff(H,k-nlze,i)));
+      p1 = (GEN)z[i+nlze]; if (!signe(p1)) continue;
+
+      if (h) p1 = gdivent(p1,h);
+      for (k=1; k<=nlze; k++) z[k] = lsubii((GEN)z[k], mulii(p1, (GEN)Di[k]));
+      for (   ; k<=lig;  k++) z[k] = lsubii((GEN)z[k], mulii(p1, (GEN)Hi[k-nlze]));
       Cnew[j] = lsub((GEN)Cnew[j], gmul(p1, (GEN)Cnew[i+zc]));
     }
     if (low_stack(lim, stack_lim(av,1)))
     {
-      GEN *gptr[2]; gptr[0]=&Cnew; gptr[1]=&B;
       if(DEBUGMEM>1) err(warnmem,"hnffinal, i = %ld",i);
-      gerepilemany(av,gptr,2);
+      gerepileall(av, 2, &Cnew, &B);
     }
   }
   p1 = cgetg(lnz+1,t_VEC); p2 = perm + nlze;
@@ -1452,16 +1521,16 @@ hnffinal(GEN matgen,GEN perm,GEN* ptdep,GEN* ptB,GEN* 
     GEN z = (GEN)H[j];
     if (diagH1[j])
     { /* hit exactly s times */
-      i1++; p1 = cgetg(lig+1,t_COL); Bnew[i1] = (long)p1;
-      C[i1+col] = Cnew[j+zc];
+      i1++; C[i1+col] = Cnew[j+zc];
+      p1 = cgetg(lig+1,t_COL); Bnew[i1] = (long)p1;
       for (i=1; i<=nlze; i++) p1[i] = coeff(dep,i,j);
       p1 += nlze;
     }
     else
     {
-      j1++; p1 = cgetg(lnz+1,t_COL); Hnew[j1] = (long)p1;
-      C[j1+zc] = Cnew[j+zc];
-      if (nlze) depnew[j1] = dep[j];
+      j1++; C[j1+zc] = Cnew[j+zc];
+      p1 = cgetg(lnz+1,t_COL); Hnew[j1] = (long)p1;
+      depnew[j1] = dep[j];
     }
     for (i=k=1; k<=lnz; i++)
       if (!diagH1[i]) p1[k++] = z[i];
@@ -1484,30 +1553,21 @@ hnffinal(GEN matgen,GEN perm,GEN* ptdep,GEN* ptB,GEN* 
       fprintferr("mit = %Z\nB = %Z\nC = %Z\n", Hnew, Bnew, C);
     }
   }
-  if (nlze) *ptdep = depnew;
+  *ptdep = depnew;
   *ptC = C;
   *ptB = Bnew; return Hnew;
 }
 
 /* for debugging */
 static void
-p_mat(long **mat, long *perm, long k0)
+p_mat(long **mat, GEN perm, long k)
 {
-  long av=avma, i,j;
-  GEN p1, matj, matgen;
-  long co = lg(mat);
-  long li = lg(perm);
-
+  gpmem_t av = avma;
+  perm = vecextract_i(perm, k+1, lg(perm)-1);
   fprintferr("Permutation: %Z\n",perm);
-  matgen = cgetg(co,t_MAT);
-  for (j=1; j<co; j++)
-  {
-    p1 = cgetg(li-k0,t_COL); matgen[j]=(long)p1;
-    p1 -= k0; matj = mat[j];
-    for (i=k0+1; i<li; i++) p1[i] = lstoi(matj[perm[i]]);
-  }
-  if (DEBUGLEVEL > 6) fprintferr("matgen = %Z\n",matgen);
-  avma=av;
+  if (DEBUGLEVEL > 6)
+    fprintferr("matgen = %Z\n", small_to_mat( rowextract_p((GEN)mat, perm) ));
+  avma = av;
 }
 
 static GEN
@@ -1524,7 +1584,7 @@ col_dup(long n, GEN col)
 **   C   = r x (co-1) matrix of GEN
 **   perm= permutation vector (length li-1), indexing the rows of mat: easier
 **     to maintain perm than to copy rows. For columns we can do it directly
-**     using e.g. swap(mat[i], mat[j])
+**     using e.g. lswap(mat[i], mat[j])
 **   k0 = integer. The k0 first lines of mat are dense, the others are sparse.
 **     Also if k0=0, mat is modified in place [from mathnfspec], otherwise
 **     it is left alone [from buchall]
@@ -1536,11 +1596,11 @@ col_dup(long n, GEN col)
 GEN
 hnfspec(long** mat0, GEN perm, GEN* ptdep, GEN* ptB, GEN* ptC, long k0)
 {
-  long av=avma,av2,*p,i,j,k,lk0,col,lig,*matj, **mat;
-  long n,s,t,lim,nlze,lnz,nr;
+  gpmem_t av=avma,av2,lim;
+  long *p,i,j,k,lk0,col,lig,*matj, **mat;
+  long n,s,t,nlze,lnz,nr;
   GEN p1,p2,matb,matbnew,vmax,matt,T,extramat;
   GEN B,H,dep,permpro;
-  GEN *gptr[4];
   long co = lg(mat0);
   long li = lg(perm); /* = lg(mat0[1]) */
   int updateT = 1;
@@ -1579,13 +1639,13 @@ hnfspec(long** mat0, GEN perm, GEN* ptdep, GEN* ptB, G
     switch( count(mat,perm[i],col,&n) )
     {
       case 0: /* move zero lines between k0+1 and lk0 */
-	lk0++; swap(perm[i], perm[lk0]);
+	lk0++; lswap(perm[i], perm[lk0]);
         i=lig; continue;
 
       case 1: /* move trivial generator between lig+1 and li */
-	swap(perm[i], perm[lig]);
-        swap(T[n], T[col]);
-	gswap(mat[n], mat[col]); p = mat[col];
+	lswap(perm[i], perm[lig]);
+        lswap(T[n], T[col]);
+	swap(mat[n], mat[col]); p = mat[col];
 	if (p[perm[lig]] < 0) /* = -1 */
 	{ /* convert relation -g = 0 to g = 0 */
 	  for (i=lk0+1; i<lig; i++) p[perm[i]] = -p[perm[i]];
@@ -1607,29 +1667,26 @@ hnfspec(long** mat0, GEN perm, GEN* ptdep, GEN* ptB, G
   }
 
 #define absmax(s,z) {long _z; _z = labs(z); if (_z > s) s = _z;}
-
-#if 0 /* TODO: check, and put back in */
   /* Get rid of all lines containing only 0 and ± 1, keeping track of column
    * operations in T. Leave the rows 1..lk0 alone [up to k0, coeff
-   * explosion, between k0+1 and lk0, row is 0]
-   */
+   * explosion, between k0+1 and lk0, row is 0] */
   s = 0;
   while (lig > lk0 && s < (HIGHBIT>>1))
   {
     for (i=lig; i>lk0; i--)
-      if (count(mat,perm[i],col,&n) >= 0) break;
+      if (count(mat,perm[i],col,&n) > 0) break;
     if (i == lk0) break;
 
     /* only 0, ±1 entries, at least 2 of them non-zero */
-    swap(perm[i], perm[lig]);
-    swap(T[n], T[col]); p1 = (GEN)T[col];
-    gswap(mat[n], mat[col]); p = mat[col];
+    lswap(perm[i], perm[lig]);
+    swap(mat[n], mat[col]); p = mat[col];
+    lswap(T[n], T[col]); p1 = (GEN)T[col];
     if (p[perm[lig]] < 0)
     {
       for (i=lk0+1; i<=lig; i++) p[perm[i]] = -p[perm[i]];
-      T[col] = lneg(p1);
+      p1 = gneg(p1); T[col] = (long)p1;
     }
-    for (j=1; j<n; j++)
+    for (j=1; j<col; j++)
     {
       matj = mat[j];
       if (! (t = matj[perm[lig]]) ) continue;
@@ -1653,7 +1710,6 @@ hnfspec(long** mat0, GEN perm, GEN* ptdep, GEN* ptB, G
       T = gerepilecopy(av2, T);
     }
   }
-#endif
   /* As above with lines containing a ±1 (no other assumption).
    * Stop when single precision becomes dangerous */
   for (j=1; j<=col; j++)
@@ -1668,10 +1724,10 @@ hnfspec(long** mat0, GEN perm, GEN* ptdep, GEN* ptB, G
       if ( (n = count2(mat,perm[i],col)) ) break;
     if (i == lk0) break;
 
-    swap(perm[i], perm[lig]);
-    swap(vmax[n], vmax[col]);
-    gswap(mat[n], mat[col]); p = mat[col];
-    swap(T[n], T[col]); p1 = (GEN)T[col];
+    lswap(vmax[n], vmax[col]);
+    lswap(perm[i], perm[lig]);
+    swap(mat[n], mat[col]); p = mat[col];
+    lswap(T[n], T[col]); p1 = (GEN)T[col];
     if (p[perm[lig]] < 0)
     {
       for (i=lk0+1; i<=lig; i++) p[perm[i]] = -p[perm[i]];
@@ -1709,7 +1765,7 @@ END2: /* clean up mat: remove everything to the right 
   if (DEBUGLEVEL>5)
   {
     fprintferr("    after phase2:\n");
-    p_mat(mat,perm,k0);
+    p_mat(mat,perm,lk0);
   }
   for (i=li-2; i>lig; i--)
   {
@@ -1748,10 +1804,10 @@ END2: /* clean up mat: remove everything to the right 
     {
       if(DEBUGMEM>1) err(warnmem,"hnfspec[3], i = %ld", i);
       for (j=1; j<co; j++) setlg(matb[j], i0+1); /* bottom can be forgotten */
-      gptr[0]=&T; gptr[1]=&matb; gerepilemany(av2,gptr,2);
+      gerepileall(av2, 2, &T, &matb);
     }
   }
-  gptr[0]=&T; gptr[1]=&matb; gerepilemany(av2,gptr,2);
+  gerepileall(av2, 2, &T, &matb);
   if (DEBUGLEVEL>5)
   {
     fprintferr("    matb cleaned up (using Id block)\n");
@@ -1829,10 +1885,7 @@ END2: /* clean up mat: remove everything to the right 
   *ptdep = dep;
   *ptB = B;
   H = hnffinal(matbnew,perm,ptdep,ptB,ptC);
-  gptr[0]=ptC;
-  gptr[1]=ptdep;
-  gptr[2]=ptB;
-  gptr[3]=&H; gerepilemany(av,gptr,4);
+  gerepileall(av, 4, ptC, ptdep, ptB, &H);
   if (DEBUGLEVEL)
     msgtimer("hnfspec [%ld x %ld] --> [%ld x %ld]",li-1,co-1, lig-1,col-1);
   return H;
@@ -1891,15 +1944,13 @@ GEN
 hnfadd(GEN H, GEN perm, GEN* ptdep, GEN* ptB, GEN* ptC, /* cf hnfspec */
        GEN extramat,GEN extraC)
 {
-  GEN p1,p2,p3,matb,extratop,Cnew,permpro;
-  GEN B=*ptB, C=*ptC, dep=*ptdep, *gptr[4];
-  long av = avma, i,j,lextra,colnew;
-  long li = lg(perm);
-  long co = lg(C);
-  long lB = lg(B);
-  long lig = li - lB;
-  long col = co - lB;
+  GEN matb, extratop, Cnew, permpro, B = *ptB, C = *ptC, dep = *ptdep;
+  gpmem_t av = avma;
+  long i;
   long lH = lg(H)-1;
+  long lB = lg(B)-1;
+  long li = lg(perm)-1, lig = li - lB;
+  long co = lg(C)-1,    col = co - lB;
   long nlze = lH? lg(dep[1])-1: lg(B[1])-1;
 
   if (DEBUGLEVEL>5)
@@ -1914,89 +1965,37 @@ hnfadd(GEN H, GEN perm, GEN* ptdep, GEN* ptB, GEN* ptC
   *       [--------|-----] lig
   *       [   0    | Id  ]
   *       [        |     ] li */
-  lextra = lg(extramat)-1;
-  extratop = cgetg(lextra+1,t_MAT); /* [1..lig] part (top) */
-  p2 = cgetg(lextra+1,t_MAT); /* bottom */
-  for (j=1; j<=lextra; j++)
-  {
-    GEN z = (GEN)extramat[j];
-    p1=cgetg(lig+1,t_COL); extratop[j] = (long)p1;
-    for (i=1; i<=lig; i++) p1[i] = z[i];
-    p1=cgetg(lB,t_COL); p2[j] = (long)p1;
-    p1 -= lig;
-    for (   ; i<li; i++) p1[i] = z[i];
-  }
-  if (li-1 != lig)
+  extratop = rowextract_i(extramat, 1, lig);
+  if (li != lig)
   { /* zero out bottom part, using the Id block */
-    GEN A = cgetg(lB,t_MAT);
-    for (j=1; j<lB; j++) A[j] = C[j+col];
-    extraC   = gsub(extraC,  gmul(A,p2));
-    extratop = gsub(extratop,gmul(B,p2));
+    GEN A = vecextract_i(C, col+1, co);
+    GEN c = rowextract_i(extramat, lig+1, li);
+    extraC   = gsub(extraC,   gmul(A, c));
+    extratop = gsub(extratop, gmul(B, c));
   }
 
-  colnew = lH + lextra;
-  extramat = cgetg(colnew+1,t_MAT);
-  Cnew = cgetg(lB+colnew,t_MAT);
-  for (j=1; j<=lextra; j++)
-  {
-    extramat[j] = extratop[j];
-    Cnew[j] = extraC[j];
-  }
-  for (   ; j<=colnew; j++)
-  {
-    p1 = cgetg(lig+1,t_COL); extramat[j] = (long)p1;
-    p2 = (GEN)dep[j-lextra]; for (i=1; i<=nlze; i++) p1[i] = p2[i];
-    p2 = (GEN)  H[j-lextra]; for (   ; i<=lig ; i++) p1[i] = p2[i-nlze];
-  }
-  for (j=lextra+1; j<lB+colnew; j++)
-    Cnew[j] = C[j-lextra+col-lH];
-  if (DEBUGLEVEL>5)
-  {
-    fprintferr("    1st phase done\n");
-    if (DEBUGLEVEL>6) fprintferr("extramat = %Z\n",extramat);
-  }
+  extramat = concatsp(extratop, vconcat(dep, H));
+  Cnew     = concatsp(extraC, vecextract_i(C, col-lH+1, co));
+  if (DEBUGLEVEL>5) fprintferr("    1st phase done\n");
   permpro = imagecomplspec(extramat, &nlze);
-  p1 = new_chunk(lig+1);
-  for (i=1; i<=lig; i++) p1[i] = perm[permpro[i]];
-  for (i=1; i<=lig; i++) perm[i] = p1[i];
+  extramat = rowextract_p(extramat, permpro);
+  *ptB     = rowextract_p(B,        permpro);
+  /* perm o= permpro */
+  permpro = vecextract_p(perm, permpro);
+  for (i=1; i<=lig; i++) perm[i] = permpro[i];
 
-  matb = cgetg(colnew+1,t_MAT);
-  dep = cgetg(colnew+1,t_MAT);
-  for (j=1; j<=colnew; j++)
-  {
-    GEN z = (GEN)extramat[j];
-    p1=cgetg(nlze+1,t_COL); dep[j]=(long)p1;
-    p2=cgetg(lig+1-nlze,t_COL); matb[j]=(long)p2;
-    p2 -= nlze;
-    for (i=1; i<=nlze; i++) p1[i] = z[permpro[i]];
-    for (   ; i<= lig; i++) p2[i] = z[permpro[i]];
-  }
-  p3 = cgetg(lB,t_MAT);
-  for (j=1; j<lB; j++)
-  {
-    p2 = (GEN)B[j];
-    p1 = cgetg(lig+1,t_COL); p3[j] = (long)p1;
-    for (i=1; i<=lig; i++) p1[i] = p2[permpro[i]];
-  }
-  B = p3;
+  *ptdep  = rowextract_i(extramat, 1, nlze);
+  matb    = rowextract_i(extramat, nlze+1, lig);
   if (DEBUGLEVEL>5) fprintferr("    2nd phase done\n");
-  *ptdep = dep;
-  *ptB = B;
   H = hnffinal(matb,perm,ptdep,ptB,&Cnew);
-  p1 = cgetg(co+lextra,t_MAT);
-  for (j=1; j <= col-lH; j++)   p1[j] = C[j];
-  C = Cnew - (col-lH);
-  for (   ; j < co+lextra; j++) p1[j] = C[j];
+  *ptC = concatsp(vecextract_i(C, 1, col-lH), Cnew);
 
-  gptr[0]=ptC; *ptC=p1;
-  gptr[1]=ptdep;
-  gptr[2]=ptB; 
-  gptr[3]=&H; gerepilemany(av,gptr,4);
   if (DEBUGLEVEL)
   {
-    if (DEBUGLEVEL>7) fprintferr("mit = %Z\nC = %Z\n",H,*ptC);
-    msgtimer("hnfadd (%ld)",lextra);
+    if (DEBUGLEVEL>7) fprintferr("H = %Z\nC = %Z\n",H,*ptC);
+    msgtimer("hnfadd (%ld)", lg(extramat)-1);
   }
+  gerepileall(av, 4, ptC, ptdep, ptB, &H);
   return H;
 }
 
@@ -2014,7 +2013,7 @@ ZV_elem(GEN aj, GEN ak, GEN A, GEN U, long j, long k)
 {
   GEN p1,u,v,d;
 
-  if (!signe(ak)) { swap(A[j],A[k]); if (U) swap(U[j],U[k]); return; }
+  if (!signe(ak)) { lswap(A[j],A[k]); if (U) lswap(U[j],U[k]); return; }
   d = bezout(aj,ak,&u,&v);
   /* frequent special case (u,v) = (1,0) or (0,1) */
   if (!signe(u))
@@ -2029,10 +2028,10 @@ ZV_elem(GEN aj, GEN ak, GEN A, GEN U, long j, long k)
   { /* aj | ak */
     p1 = negi(divii(ak,aj));
     A[k]   = (long)ZV_lincomb(gun, p1, (GEN)A[k], (GEN)A[j]);
-    swap(A[j], A[k]);
+    lswap(A[j], A[k]);
     if (U) {
       U[k] = (long)ZV_lincomb(gun, p1, (GEN)U[k], (GEN)U[j]);
-      swap(U[j], U[k]);
+      lswap(U[j], U[k]);
     }
     return;
   }
@@ -2077,7 +2076,8 @@ ZM_reduce(GEN A, GEN U, long i, long j0)
 GEN
 hnf0(GEN A, int remove)
 {
-  long av0 = avma, s,li,co,av,tetpil,i,j,k,def,ldef,lim;
+  gpmem_t av0 = avma, av, lim;
+  long s,li,co,i,j,k,def,ldef;
   GEN denx,a,p1;
 
   if (typ(A) == t_VEC) return hnf_special(A,remove);
@@ -2100,7 +2100,7 @@ hnf0(GEN A, int remove)
       if (low_stack(lim, stack_lim(av,1)))
       {
         if (DEBUGMEM>1) err(warnmem,"hnf[1]. i=%ld",i);
-        tetpil=avma; A=gerepile(av,tetpil,ZM_copy(A));
+        A = gerepilecopy(av, A);
       }
     }
     p1 = gcoeff(A,i,def); s = signe(p1);
@@ -2115,7 +2115,7 @@ hnf0(GEN A, int remove)
     if (low_stack(lim, stack_lim(av,1)))
     {
       if (DEBUGMEM>1) err(warnmem,"hnf[2]. i=%ld",i);
-      tetpil=avma; A=gerepile(av,tetpil,ZM_copy(A));
+      A = gerepilecopy(av, A);
     }
   }
   if (remove)
@@ -2124,9 +2124,8 @@ hnf0(GEN A, int remove)
       if (!gcmp0((GEN)A[j])) A[i++] = A[j];
     setlg(A,i);
   }
-  tetpil=avma;
   A = denx? gdiv(A,denx): ZM_copy(A);
-  return gerepile(av0,tetpil,A);
+  return gerepileupto(av0, A);
 }
 
 GEN
@@ -2137,7 +2136,7 @@ hnf(GEN x) { return hnf0(x,1); }
 static GEN
 allhnfmod(GEN x,GEN dm,int flag)
 {
-  ulong av,tetpil,lim;
+  gpmem_t av, lim;
   long li,co,i,j,k,def,ldef,ldm;
   GEN a,b,w,p1,p2,d,u,v;
 
@@ -2184,7 +2183,7 @@ allhnfmod(GEN x,GEN dm,int flag)
       if (low_stack(lim, stack_lim(av,1)))
       {
         if (DEBUGMEM>1) err(warnmem,"allhnfmod[1]. i=%ld",i);
-	tetpil=avma; x=gerepile(av,tetpil,ZM_copy(x));
+	x = gerepilecopy(av, x);
       }
     }
   w = cgetg(li,t_MAT); b = dm;
@@ -2214,12 +2213,11 @@ allhnfmod(GEN x,GEN dm,int flag)
       if (low_stack(lim, stack_lim(av,1)))
       {
         if (DEBUGMEM>1) err(warnmem,"allhnfmod[2]. i=%ld", i);
-        tetpil=avma; w=gerepile(av,tetpil,ZM_copy(w));
-        diag = gcoeff(w,i,i);
+        gerepileall(av, 2, &w, &dm); diag = gcoeff(w,i,i);
       }
     }
   }
-  tetpil=avma; return gerepile(av,tetpil,ZM_copy(w));
+  return gerepilecopy(av, w);
 }
 
 GEN
@@ -2238,7 +2236,7 @@ hnfmodid(GEN x, GEN p) { return allhnfmod(x,p,0); }
 static GEN
 mynegi(GEN x)
 {
-  static long mun[]={evaltyp(t_INT)|m_evallg(3),evalsigne(-1)|evallgefint(3),1};
+  static long mun[]={evaltyp(t_INT)|_evallg(3),evalsigne(-1)|evallgefint(3),1};
   long s = signe(x);
 
   if (!s) return gzero;
@@ -2247,22 +2245,6 @@ mynegi(GEN x)
   setsigne(x,-s); return x;
 }
 
-/* We assume y > 0 */
-static GEN
-divnearest(GEN x, GEN y)
-{
-  GEN r, q = dvmdii(x,y,&r);
-  long av = avma, s=signe(r), t;
-
-  if (!s) { cgiv(r); return q; }
-  if (s<0) r = mynegi(r);
-
-  y = shifti(y,-1); t = cmpii(r,y);
-  avma=av; cgiv(r);
-  if (t < 0 || (t == 0 && s > 0)) return q;
-  return addsi(s,q);
-}
-
 static void
 Minus(long j, GEN **lambda)
 {
@@ -2313,33 +2295,30 @@ findi(GEN M)
   return 0;
 }
 
+static long
+findi_normalize(GEN Aj, GEN B, long j, GEN **lambda)
+{
+  long r = findi(Aj);
+  if (r && signe(Aj[r]) < 0)
+  {
+    neg_col(Aj); if (B) neg_col((GEN)B[j]);
+    Minus(j,lambda);
+  }
+  return r;
+}
+
 static void
 reduce2(GEN A, GEN B, long k, long j, long *row, GEN **lambda, GEN *D)
 {
   GEN q;
   long i, row0, row1;
 
-  row0 = findi((GEN)A[j]);
-  if (row0 && signe(gcoeff(A,row0,j)) < 0)
-  {
-    neg_col((GEN)A[j]);
-    if (B) neg_col((GEN)B[j]);
-    Minus(j,lambda);
-  }
-
-  row1 = findi((GEN)A[k]);
-  if (row1 && signe(gcoeff(A,row1,k)) < 0)
-  {
-    neg_col((GEN)A[k]);
-    if (B) neg_col((GEN)B[k]);
-    Minus(k,lambda);
-  }
-
-  row[0]=row0; row[1]=row1;
+  row[0] = row0 = findi_normalize((GEN)A[j], B,j,lambda);
+  row[1] = row1 = findi_normalize((GEN)A[k], B,k,lambda);
   if (row0)
     q = truedvmdii(gcoeff(A,row0,k), gcoeff(A,row0,j), NULL);
   else if (absi_cmp(shifti(lambda[k][j], 1), D[j]) > 0)
-    q = divnearest(lambda[k][j], D[j]);
+    q = diviiround(lambda[k][j], D[j]);
   else
     return;
 
@@ -2371,17 +2350,15 @@ reduce2(GEN A, GEN B, long k, long j, long *row, GEN *
   }
 }
 
-#define SWAP(x,y) {GEN _t=y; y=x; x=_t;}
-
 static void
 hnfswap(GEN A, GEN B, long k, GEN **lambda, GEN *D)
 {
   GEN t,p1,p2;
   long i,j,n = lg(A);
 
-  swap(A[k],A[k-1]);
-  if (B) swap(B[k],B[k-1]);
-  for (j=k-2; j; j--) SWAP(lambda[k-1][j], lambda[k][j]);
+  lswap(A[k],A[k-1]);
+  if (B) lswap(B[k],B[k-1]);
+  for (j=k-2; j; j--) swap(lambda[k-1][j], lambda[k][j]);
   for (i=k+1; i<n; i++)
   {
     p1 = mulii(lambda[i][k-1], D[k]);
@@ -2399,33 +2376,7 @@ hnfswap(GEN A, GEN B, long k, GEN **lambda, GEN *D)
   D[k-1] = divii(addii(p1,p2), D[k-1]);
 }
 
-/* A[k] = 0,  A[nz] != 0,  k > nz,  A[j] = 0 for all j < nz.
- * "Deep insert" A[k] at nz */
-static void
-trivswap(GEN *A, long nz, long k, GEN **lambda, GEN *D)
-{
-  GEN p1;
-  long i,j,n = lg(A);
-
-  p1 = A[nz]; /* cycle A */
-  for (j = nz; j < k; j++) SWAP(A[j+1], p1); 
-  A[nz] = p1; /* = [0...0] */
-
-  p1 = D[nz]; /* cycle D */
-  for (j = nz; j < k; j++) SWAP(D[j+1], p1);
-  D[nz] = gun;
-
-  for (j=k-1; j>=nz; j--) /* cycle lambda */
-    for (i=k-1; i>=nz; i--) lambda[i+1][j+1] = lambda[i][j];
-  for (j=n-1; j>k; j--)
-    for (i=k-1; i>=nz; i--) 
-    {
-      lambda[i+1][j] = lambda[i][j];
-      lambda[j][i+1] = lambda[j][i];
-    }
-  for (i=1; i<n; i++) lambda[nz][i] = lambda[i][nz] = gzero;
-}
-
+/* reverse row order in matrix A */
 static GEN
 fix_rows(GEN A)
 {
@@ -2443,88 +2394,75 @@ fix_rows(GEN A)
 }
 
 GEN
-hnflll_i(GEN A, GEN *ptB)
+hnflll_i(GEN A, GEN *ptB, int remove)
 {
-  ulong av = avma, lim = stack_lim(av,3);
+  gpmem_t av = avma, lim = stack_lim(av,3);
   long m1 = 1, n1 = 1; /* alpha = m1/n1. Maybe 3/4 here ? */
   long row[2], do_swap,i,n,k;
-  long nzcol = 1; /* index of 1st (maybe) non-0 col [only updated when !B] */
   GEN z,B, **lambda, *D;
-  GEN *gptr[4];
 
   if (typ(A) != t_MAT) err(typeer,"hnflll");
   n = lg(A);
   A = ZM_copy(fix_rows(A));
   B = ptB? idmat(n-1): NULL;
-  D = (GEN*) cgetg(n+1, t_VEC); D++; /* hack: need a "sentinel" D[0] */
-  if (n == 2) /* handle trivial case: return negative diag coeff otherwise */
-  {
-    i = findi((GEN)A[1]);
-    if (i && signe(gcoeff(A,i,1)) < 0)
-    {
-      neg_col((GEN)A[1]);
-      if (B) neg_col((GEN)B[1]);
-    }
-  }
-  lambda = (GEN**) cgetg(n,t_MAT);
-  for (i=1; i<n; i++) { D[i] = gun; lambda[i] = (GEN*)zerocol(n-1); }
-  D[0] = gun; k = 2;
+  D = (GEN*)cgetg(n+1,t_VEC); lambda = (GEN**) cgetg(n,t_MAT);
+  D++;
+  for (i=0; i<n; i++) D[i] = gun;
+  for (i=1; i<n; i++) lambda[i] = (GEN*)zerocol(n-1);
+  k = 2;
   while (k < n)
   {
     reduce2(A,B,k,k-1,row,lambda,D);
-    if (!B)
-    { /* not interested in B. Eliminate 0 cols */
-      if (!row[0] || !row[1])
-      {
-        while (!findi((GEN)A[nzcol]) && nzcol < n) nzcol++;
-        /* A[nzcol] != 0,  A[i] = 0 for i < nzcol */
-        if (!row[0] && k-1 > nzcol)
-          {trivswap((GEN*)A,nzcol,k-1, lambda,D); nzcol++;}
-        if (!row[1] && k   > nzcol)
-          {trivswap((GEN*)A,nzcol,k  , lambda,D); nzcol++;}
-        if (k <= nzcol) k = nzcol+1;
-        continue;
-      }
-      do_swap = (row[0] && row[0] <= row[1]);
+    if (row[0])
+      do_swap = (!row[1] || row[0] <= row[1]);
+    else if (!row[1])
+    { /* row[0] == row[1] == 0 */
+      gpmem_t av1 = avma;
+      z = addii(mulii(D[k-2],D[k]), sqri(lambda[k][k-1]));
+      do_swap = (cmpii(mulsi(n1,z), mulsi(m1,sqri(D[k-1]))) < 0);
+      avma = av1;
     }
     else
-    {
-      if (row[0])
-        do_swap = (!row[1] || row[0] <= row[1]);
-      else if (row[1])
-        do_swap = 0;
-      else
-      { /* row[0] == row[1] == 0 */
-        ulong av1 = avma;
-        z = addii(mulii(D[k-2],D[k]), sqri(lambda[k][k-1]));
-        do_swap = (cmpii(mulsi(n1,z), mulsi(m1,sqri(D[k-1]))) < 0);
-        avma = av1;
-      }
-    }
+      do_swap = 0;
     if (do_swap)
     {
       hnfswap(A,B,k,lambda,D);
-      if (k > nzcol+1) k--;
+      if (k > 2) k--;
     }
     else
     {
-      for (i=k-2; i>=nzcol; i--) reduce2(A,B,k,i,row,lambda,D);
+      for (i=k-2; i; i--)
+      {
+        reduce2(A,B,k,i,row,lambda,D);
+        if (low_stack(lim, stack_lim(av,3)))
+        {
+          GEN b = (GEN)(D-1);
+          if (DEBUGMEM) err(warnmem,"hnflll (reducing), i = %ld",i);
+          gerepileall(av, B? 4: 3, &A, (GEN*)&lambda, &b, &B);
+          D = (GEN*)(b+1);
+        }
+      }
       k++;
     }
     if (low_stack(lim, stack_lim(av,3)))
     {
-      GEN a = (GEN)lambda, b = (GEN)(D-1); /* gcc -Wall */
-      gptr[0]=&A; gptr[1]=&a; gptr[2]=&b; gptr[3]=&B; 
-      if (DEBUGMEM) err(warnmem,"hnflll, k = %ld / %ld",k,n);
-      gerepilemany(av,gptr,B? 4: 3); lambda = (GEN**)a; D = (GEN*)(b+1);
+      GEN b = (GEN)(D-1);
+      if (DEBUGMEM) err(warnmem,"hnflll, k = %ld / %ld",k,n-1);
+      gerepileall(av, B? 4: 3, &A, (GEN*)&lambda, &b, &B);
+      D = (GEN*)(b+1);
     }
   }
-  for (i=nzcol; i<n; i++)
-    if (findi((GEN)A[i])) break;
-  i--; A += i; A[0] = evaltyp(t_MAT)|evallg(n-i);
+  /* handle trivial case: return negative diag coeff otherwise */
+  if (n == 2) (void)findi_normalize((GEN)A[1], B,1,lambda);
   A = fix_rows(A);
-  gptr[0] = &A; gptr[1] = &B;
-  gerepilemany(av, gptr, B? 2: 1);
+  if (remove)
+  {
+    for (i = 1; i < n; i++)
+      if (findi((GEN)A[i])) break;
+    i--;
+    A += i; A[0] = evaltyp(t_MAT) | evallg(n-i);
+  }
+  gerepileall(av, B? 2: 1, &A, &B);
   if (B) *ptB = B;
   return A;
 }
@@ -2533,7 +2471,7 @@ GEN
 hnflll(GEN A)
 {
   GEN B, z = cgetg(3, t_VEC);
-  z[1] = (long)hnflll_i(A, &B);
+  z[1] = (long)hnflll_i(A, &B, 0);
   z[2] = (long)B; return z;
 }
 
@@ -2546,21 +2484,21 @@ reduce1(GEN A, GEN B, long k, long j, GEN **lambda, GE
   long i;
 
   if (signe(A[j]))
-    q = divnearest((GEN)A[k],(GEN)A[j]);
+    q = diviiround((GEN)A[k],(GEN)A[j]);
   else if (absi_cmp(shifti(lambda[k][j], 1), D[j]) > 0)
-    q = divnearest(lambda[k][j], D[j]);
+    q = diviiround(lambda[k][j], D[j]);
   else
     return;
 
-  if (! gcmp0(q))
+  if (signe(q))
   {
+    GEN *Lk = lambda[k], *Lj = lambda[j];
     q = mynegi(q);
     A[k] = laddii((GEN)A[k], mulii(q,(GEN)A[j]));
     elt_col((GEN)B[k],(GEN)B[j],q);
-    lambda[k][j] = addii(lambda[k][j], mulii(q,D[j]));
+    Lk[j] = addii(Lk[j], mulii(q,D[j]));
     for (i=1; i<j; i++)
-      if (signe(lambda[j][i]))
-        lambda[k][i] = addii(lambda[k][i], mulii(q,lambda[j][i]));
+      if (signe(Lj[i])) Lk[i] = addii(Lk[i], mulii(q,Lj[i]));
   }
 }
 
@@ -2568,17 +2506,16 @@ GEN
 extendedgcd(GEN A)
 {
   long m1 = 1, n1 = 1; /* alpha = m1/n1. Maybe 3/4 here ? */
-  long av = avma, tetpil, do_swap,i,j,n,k;
-  GEN p1,p2,B, **lambda, *D;
+  gpmem_t av = avma;
+  long do_swap,i,n,k;
+  GEN z,B, **lambda, *D;
 
-  n = lg(A); B = idmat(n-1); A = ZM_copy(A);
-  D = (GEN*) cgeti(n); lambda = (GEN**) cgetg(n,t_MAT);
+  n = lg(A);
+  A = dummycopy(A);
+  B = idmat(n-1);
+  D = (GEN*)new_chunk(n); lambda = (GEN**) cgetg(n,t_MAT);
   for (i=0; i<n; i++) D[i] = gun;
-  for (i=1; i<n; i++)
-  {
-    lambda[i] = (GEN*) cgetg(n,t_COL);
-    for(j=1;j<n;j++) lambda[i][j] = gzero;
-  }
+  for (i=1; i<n; i++) lambda[i] = (GEN*)zerocol(n-1);
   k = 2;
   while (k < n)
   {
@@ -2586,19 +2523,17 @@ extendedgcd(GEN A)
     if (signe(A[k-1])) do_swap = 1;
     else if (!signe(A[k]))
     {
-      long av1=avma;
-      p1 = addii(mulii(D[k-2],D[k]), sqri(lambda[k][k-1]));
-      p1 = mulsi(n1,p1);
-      p2 = mulsi(m1,sqri(D[k-1]));
-      do_swap = (cmpii(p1,p2) < 0);
-      avma=av1;
+      gpmem_t av1 = avma;
+      z = addii(mulii(D[k-2],D[k]), sqri(lambda[k][k-1]));
+      do_swap = (cmpii(mulsi(n1,z), mulsi(m1,sqri(D[k-1]))) < 0);
+      avma = av1;
     }
     else do_swap = 0;
 
     if (do_swap)
     {
       hnfswap(A,B,k,lambda,D);
-      if (k>2) k--;
+      if (k > 2) k--;
     }
     else
     {
@@ -2606,16 +2541,15 @@ extendedgcd(GEN A)
       k++;
     }
   }
-  if (signe(A[n-1])<0)
+  if (signe(A[n-1]) < 0)
   {
     A[n-1] = (long)mynegi((GEN)A[n-1]);
     neg_col((GEN)B[n-1]);
   }
-  tetpil = avma;
-  p1 = cgetg(3,t_VEC);
-  p1[1]=lcopy((GEN)A[n-1]);
-  p1[2]=lcopy(B);
-  return gerepile(av,tetpil,p1);
+  z = cgetg(3,t_VEC);
+  z[1] = A[n-1];
+  z[2] = (long)B;
+  return gerepilecopy(av, z);
 }
 
 /* HNF with permutation. TODO: obsolete, replace with hnflll. */
@@ -2623,7 +2557,8 @@ GEN
 hnfperm(GEN A)
 {
   GEN U,c,l,perm,d,u,p,q,y,b;
-  long r,t,i,j,j1,k,m,n,av=avma,av1,tetpil,lim;
+  gpmem_t av=avma,av1,tetpil,lim;
+  long r,t,i,j,j1,k,m,n;
 
   if (typ(A) != t_MAT) err(typeer,"hnfperm");
   n = lg(A)-1;
@@ -2696,9 +2631,8 @@ hnfperm(GEN A)
     }
     if (low_stack(lim, stack_lim(av1,1)))
     {
-      GEN *gptr[2]; gptr[0]=&A; gptr[1]=&U;
       if (DEBUGMEM>1) err(warnmem,"hnfperm");
-      gerepilemany(av1,gptr,2);
+      gerepileall(av1, 2, &A, &U);
     }
   }
   if (r < m)
@@ -2735,8 +2669,8 @@ GEN
 hnfall_i(GEN A, GEN *ptB, long remove)
 {
   GEN B,c,h,x,p,a;
-  long m,n,r,i,j,k,li,av=avma,av1,tetpil,lim;
-  GEN *gptr[2]; 
+  gpmem_t av = avma, av1, lim;
+  long m,n,r,i,j,k,li;
 
   if (typ(A)!=t_MAT) err(typeer,"hnfall");
   n=lg(A)-1;
@@ -2766,11 +2700,8 @@ hnfall_i(GEN A, GEN *ptB, long remove)
         ZM_reduce(A,B, i,k);
         if (low_stack(lim, stack_lim(av1,1)))
         {
-          tetpil = avma;
-          A = ZM_copy(A); gptr[0]=&A;
-          if (B) { B = ZM_copy(B); gptr[1]=&B; }
           if (DEBUGMEM>1) err(warnmem,"hnfall[1], li = %ld", li);
-          gerepilemanysp(av1,tetpil,gptr,B? 2: 1);
+          gerepileall(av1, B? 2: 1, &A, &B);
         }	
       }
       x = gcoeff(A,li,j); if (signe(x)) break;
@@ -2780,8 +2711,8 @@ hnfall_i(GEN A, GEN *ptB, long remove)
     r--;
     if (j < r) /* A[j] != 0 */
     {
-      swap(A[j], A[r]);
-      if (B) swap(B[j], B[r]);
+      lswap(A[j], A[r]);
+      if (B) lswap(B[j], B[r]);
       h[j]=h[r]; h[r]=li; c[li]=r;
     }
     p = gcoeff(A,li,r);
@@ -2794,9 +2725,8 @@ hnfall_i(GEN A, GEN *ptB, long remove)
     ZM_reduce(A,B, li,r);
     if (low_stack(lim, stack_lim(av1,1)))
     {
-      GEN *gptr[2]; gptr[0]=&A; gptr[1]=&B;
       if (DEBUGMEM>1) err(warnmem,"hnfall[2], li = %ld", li);
-      gerepilemany(av1,gptr,B? 2: 1);
+      gerepileall(av1, B? 2: 1, &A, &B);
     }	
   }
   if (DEBUGLEVEL>5) fprintferr("\nhnfall, final phase: ");
@@ -2812,18 +2742,14 @@ hnfall_i(GEN A, GEN *ptB, long remove)
       ZM_reduce(A,B, i,k);
       if (low_stack(lim, stack_lim(av1,1)))
       {
-        tetpil = avma;
-        A = ZM_copy(A); gptr[0]=&A;
-        if (B) { B = ZM_copy(B); gptr[1]=&B; }
         if (DEBUGMEM>1) err(warnmem,"hnfall[3], j = %ld", j);
-        gerepilemanysp(av1,tetpil,gptr, B? 2: 1);
+        gerepileall(av1, B? 2: 1, &A, &B);
       }	
     }
   if (DEBUGLEVEL>5) fprintferr("\n");
   /* remove the first r columns */
   if (remove) { A += r; A[0] = evaltyp(t_MAT) | evallg(n-r+1); }
-  gptr[0] = &A; gptr[1] = &B;
-  gerepilemany(av, gptr, B? 2: 1);
+  gerepileall(av, B? 2: 1, &A, &B);
   if (B) *ptB = B;
   return A;
 }
@@ -2896,208 +2822,220 @@ trivsmith(long all)
   z[3]=lgetg(1,t_MAT); return z;
 }
 
-/* Return the smith normal form d of matrix x. If all != 0 return [d,u,v],
- * where d = u.x.v
- */
-static GEN
-smithall(GEN x, long all)
+static void
+snf_pile(gpmem_t av, GEN *x, GEN *U, GEN *V)
 {
-  long av,tetpil,i,j,k,l,c,fl,n,s1,s2,lim;
-  GEN p1,p2,p3,p4,z,b,u,v,d,ml,mr,mun,mdet,ys;
+  GEN *gptr[3];
+  int c = 1; gptr[0]=x;
+  if (*U) gptr[c++] = U;
+  if (*V) gptr[c++] = V;
+  gerepilemany(av,gptr,c);
+}
 
+/* Return the SNF D of matrix X. If ptU/ptV non-NULL set them to U/V
+ * to that D = UXV */
+GEN
+smithall(GEN x, GEN *ptU, GEN *ptV)
+{
+  gpmem_t av = avma, lim = stack_lim(av,1);
+  long i,j,k,l,c,n,s1,s2;
+  GEN p1,p2,p3,p4,b,u,v,d,U,V,mun,mdet,ys;
+
   if (typ(x)!=t_MAT) err(typeer,"smithall");
-  if (DEBUGLEVEL>=9) outerr(x);
-  n=lg(x)-1;
-  if (!n) return trivsmith(all);
+  if (DEBUGLEVEL>8) outerr(x);
+  n = lg(x)-1;
+  if (!n) {
+    if (ptU) *ptU = cgetg(1,t_MAT);
+    if (ptV) *ptV = cgetg(1,t_MAT);
+    return cgetg(1,t_MAT);
+  }
   if (lg(x[1]) != n+1) err(mattype1,"smithall");
   for (i=1; i<=n; i++)
     for (j=1; j<=n; j++)
       if (typ(coeff(x,i,j)) != t_INT)
         err(talker,"non integral matrix in smithall");
 
-  av = avma; lim = stack_lim(av,1);
-  x = dummycopy(x); mdet = detint(x);
-  if (ishnfall(x)) { if (all) { ml=idmat(n); mr=idmat(n); } }
+  U = ptU? gun: NULL;
+  V = ptV? gun: NULL;
+  x = dummycopy(x);
+  if (ishnfall(x))
+  {
+    mdet = dethnf_i(x);
+    if (V) V = idmat(n);
+  }
   else
   {
+    mdet = detint(x);
     if (signe(mdet))
     {
-      p1=hnfmod(x,mdet);
-      if (all) { ml=idmat(n); mr=gauss(x,p1); }
+      p1 = hnfmod(x,mdet);
+      if (V) V = gauss(x,p1);
     }
     else
     {
-      if (all)
-      {
-        p1 = hnfall0(x,0);
-        ml = idmat(n); mr = (GEN)p1[2]; p1 = (GEN)p1[1];
-      }
-      else
-        p1 = hnf0(x,0);
+      p1 = hnfall_i(x, ptV, 0);
+      if (ptV) V = *ptV;
     }
     x = p1;
   }
+  if (U) U = idmat(n);
+ 
   p1=cgetg(n+1,t_VEC); for (i=1; i<=n; i++) p1[i]=lnegi(gcoeff(x,i,i));
-  p2=sindexsort(p1); ys=cgetg(n+1,t_MAT);
+  p2=sindexsort(p1);
+  ys = cgetg(n+1,t_MAT);
   for (j=1; j<=n; j++)
   {
     p1=cgetg(n+1,t_COL); ys[j]=(long)p1;
     for (i=1; i<=n; i++) p1[i]=coeff(x,p2[i],p2[j]);
   }
   x = ys;
-  if (all)
+  if (U)
   {
-    p3=cgetg(n+1,t_MAT); p4=cgetg(n+1,t_MAT);
-    for (j=1; j<=n; j++) { p3[j]=ml[p2[j]]; p4[j]=mr[p2[j]]; }
-    ml=p3; mr=p4;
+    p1 = cgetg(n+1,t_MAT);
+    for (j=1; j<=n; j++) p1[j]=U[p2[j]];
+    U = p1;
   }
+  if (V)
+  {
+    p1 = cgetg(n+1,t_MAT);
+    for (j=1; j<=n; j++) p1[j]=V[p2[j]];
+    V = p1;
+  }
+ 
   if (signe(mdet))
   {
     p1 = hnfmod(x,mdet);
-    if (all) mr=gmul(mr,gauss(x,p1));
+    if (V) V = gmul(V, gauss(x,p1));
   }
   else
   {
-    if (all)
-    {
-      p1 = hnfall0(x,0);
-      mr = gmul(mr,(GEN)p1[2]); p1 = (GEN)p1[1];
-    }
-    else
-      p1 = hnf0(x,0);
+    p1 = hnfall_i(x,ptV,0);
+    if (ptV) V = gmul(V, *ptV);
   }
-  x=p1; mun = negi(gun);
+  x = p1; mun = negi(gun);
 
-  if (DEBUGLEVEL>7) {fprintferr("starting SNF loop");flusherr();}
-  for (i=n; i>=2; i--)
+  if (DEBUGLEVEL>7) fprintferr("starting SNF loop");
+  for (i=n; i>1; i--)
   {
-    if (DEBUGLEVEL>7) {fprintferr("\ni = %ld: ",i);flusherr();}
+    if (DEBUGLEVEL>7) fprintferr("\ni = %ld: ",i);
     for(;;)
     {
       c = 0;
       for (j=i-1; j>=1; j--)
       {
 	p1=gcoeff(x,i,j); s1 = signe(p1);
-	if (s1)
-	{
-	  p2=gcoeff(x,i,i);
-          if (!absi_cmp(p1,p2))
+	if (!s1) continue;
+       
+        p2=gcoeff(x,i,i);
+        if (!absi_cmp(p1,p2))
+        {
+          s2=signe(p2);
+          if (s1 == s2) { d=p1; u=gun; p4=gun; }
+          else
           {
-            s2=signe(p2);
-            if (s1 == s2) { d=p1; u=gun; p4=gun; }
-            else
-	    {
-              if (s2>0) { u = gun; p4 = mun; }
-              else      { u = mun; p4 = gun; }
-	      d=(s1>0)? p1: absi(p1);
-	    }
-            v = gzero; p3 = u;
+            if (s2>0) { u = gun; p4 = mun; }
+            else      { u = mun; p4 = gun; }
+            d=(s1>0)? p1: absi(p1);
           }
-          else { d=bezout(p2,p1,&u,&v); p3=divii(p2,d); p4=divii(p1,d); }
-	  for (k=1; k<=i; k++)
-	  {
-	    b=addii(mulii(u,gcoeff(x,k,i)),mulii(v,gcoeff(x,k,j)));
-	    coeff(x,k,j)=lsubii(mulii(p3,gcoeff(x,k,j)),
-	                        mulii(p4,gcoeff(x,k,i)));
-	    coeff(x,k,i)=(long)b;
-	  }
-	  if (all) update(u,v,p3,p4,(GEN*)(mr+i),(GEN*)(mr+j));
-          if (low_stack(lim, stack_lim(av,1)))
-	  {
-	    if (DEBUGMEM>1) err(warnmem,"[1]: smithall");
-	    if (all)
-	    {
-	      GEN *gptr[3]; gptr[0]=&x; gptr[1]=&ml; gptr[2]=&mr;
-	      gerepilemany(av,gptr,3);
-	    }
-	    else x=gerepileupto(av, ZM_copy(x));
-	  }
-	}
+          v = gzero; p3 = u;
+        }
+        else { d=bezout(p2,p1,&u,&v); p3=divii(p2,d); p4=divii(p1,d); }
+        for (k=1; k<=i; k++)
+        {
+          b=addii(mulii(u,gcoeff(x,k,i)),mulii(v,gcoeff(x,k,j)));
+          coeff(x,k,j)=lsubii(mulii(p3,gcoeff(x,k,j)),
+                              mulii(p4,gcoeff(x,k,i)));
+          coeff(x,k,i)=(long)b;
+        }
+        if (V) update(u,v,p3,p4,(GEN*)(V+i),(GEN*)(V+j));
+        if (low_stack(lim, stack_lim(av,1)))
+        {
+          if (DEBUGMEM>1) err(warnmem,"[1]: smithall");
+          snf_pile(av, &x,&U,&V);
+        }
       }
-      if (DEBUGLEVEL>=8) {fprintferr("; ");flusherr();}
+      if (DEBUGLEVEL>7) fprintferr("; ");
       for (j=i-1; j>=1; j--)
       {
 	p1=gcoeff(x,j,i); s1 = signe(p1);
-	if (s1)
-	{
-	  p2=gcoeff(x,i,i);
-	  if (!absi_cmp(p1,p2))
+	if (!s1) continue;
+       
+        p2=gcoeff(x,i,i);
+        if (!absi_cmp(p1,p2))
+        {
+          s2 = signe(p2);
+          if (s1 == s2) { d=p1; u=gun; p4=gun; }
+          else
           {
-            s2 = signe(p2);
-            if (s1 == s2) { d=p1; u=gun; p4=gun; }
-            else
-	    {
-              if (s2>0) { u = gun; p4 = mun; }
-              else      { u = mun; p4 = gun; }
-	      d=(s1>0)? p1: absi(p1);
-	    }
-            v = gzero; p3 = u;
+            if (s2>0) { u = gun; p4 = mun; }
+            else      { u = mun; p4 = gun; }
+            d=(s1>0)? p1: absi(p1);
           }
-          else { d=bezout(p2,p1,&u,&v); p3=divii(p2,d); p4=divii(p1,d); }
-	  for (k=1; k<=i; k++)
-	  {
-	    b=addii(mulii(u,gcoeff(x,i,k)),mulii(v,gcoeff(x,j,k)));
-	    coeff(x,j,k)=lsubii(mulii(p3,gcoeff(x,j,k)),
-	                        mulii(p4,gcoeff(x,i,k)));
-	    coeff(x,i,k)=(long)b;
-	  }
-	  if (all) update(u,v,p3,p4,(GEN*)(ml+i),(GEN*)(ml+j));
-	  c = 1;
-	}
+          v = gzero; p3 = u;
+        }
+        else { d=bezout(p2,p1,&u,&v); p3=divii(p2,d); p4=divii(p1,d); }
+        for (k=1; k<=i; k++)
+        {
+          b=addii(mulii(u,gcoeff(x,i,k)),mulii(v,gcoeff(x,j,k)));
+          coeff(x,j,k)=lsubii(mulii(p3,gcoeff(x,j,k)),
+                              mulii(p4,gcoeff(x,i,k)));
+          coeff(x,i,k)=(long)b;
+        }
+        if (U) update(u,v,p3,p4,(GEN*)(U+i),(GEN*)(U+j));
+        c = 1;
       }
       if (!c)
       {
-	b=gcoeff(x,i,i); fl=1;
-	if (signe(b))
-	{
-	  for (k=1; k<i && fl; k++)
-	    for (l=1; l<i && fl; l++)
-	      fl = (int)!signe(resii(gcoeff(x,k,l),b));
-          /* cast to (int) necessary for gcc-2.95 on sparcv9-64 (IS) */
-	  if (!fl)
-	  {
-	    k--;
-	    for (l=1; l<=i; l++)
-	      coeff(x,i,l)=laddii(gcoeff(x,i,l),gcoeff(x,k,l));
-	    if (all) ml[i]=ladd((GEN)ml[i],(GEN)ml[k]);
-	  }
-	}
-        if (fl) break;
+	b = gcoeff(x,i,i);
+	if (!signe(b)) break;
+       
+        for (k=1; k<i; k++)
+        {
+          for (l=1; l<i; l++)
+            if (signe(resii(gcoeff(x,k,l),b))) break;
+          if (l != i) break;
+        }
+        if (k == i) break;
+       
+        /* x[k,l] != 0 mod b */
+        for (l=1; l<=i; l++)
+          coeff(x,i,l) = laddii(gcoeff(x,i,l),gcoeff(x,k,l));
+        if (U) U[i] = ladd((GEN)U[i],(GEN)U[k]);
       }
       if (low_stack(lim, stack_lim(av,1)))
       {
-	if (DEBUGMEM>1) err(warnmem,"[2]: smithall");
-	if (all)
-	{
-	  GEN *gptr[3]; gptr[0]=&x; gptr[1]=&ml; gptr[2]=&mr;
-	  gerepilemany(av,gptr,3);
-	}
-	else x=gerepileupto(av,ZM_copy(x));
+        if (DEBUGMEM>1) err(warnmem,"[2]: smithall");
+        snf_pile(av, &x,&U,&V);
       }
     }
   }
-  if (DEBUGLEVEL>7) {fprintferr("\n");flusherr();}
-  if (all)
-  {
-    for (k=1; k<=n; k++)
-      if (signe(gcoeff(x,k,k))<0)
-        { mr[k]=lneg((GEN)mr[k]); coeff(x,k,k)=lnegi(gcoeff(x,k,k)); }
-    tetpil=avma; z=cgetg(4,t_VEC);
-    z[1]=ltrans(ml); z[2]=lcopy(mr); z[3]=lcopy(x);
-    return gerepile(av,tetpil,z);
-  }
-  tetpil=avma; z=cgetg(n+1,t_VEC); j=n;
-  for (k=n; k; k--)
-    if (signe(gcoeff(x,k,k))) z[j--]=labsi(gcoeff(x,k,k));
-  for (   ; j; j--) z[j]=zero;
-  return gerepile(av,tetpil,z);
+  if (DEBUGLEVEL>7) fprintferr("\n");
+  for (k=1; k<=n; k++)
+    if (signe(gcoeff(x,k,k)) < 0)
+    {
+      if (V) V[k]=lneg((GEN)V[k]);
+      coeff(x,k,k) = lnegi(gcoeff(x,k,k));
+    }
+  if (!U && !V)
+    return gerepileupto(av, mattodiagonal(x));
+   
+  if (U) U = gtrans_i(U);
+  snf_pile(av, &x,&U,&V);
+  if (ptU) *ptU = U;
+  if (ptV) *ptV = V;
+  return x;
 }
 
 GEN
-smith(GEN x) { return smithall(x,0); }
+smith(GEN x) { return smithall(x, NULL,NULL); }
 
 GEN
-smith2(GEN x) { return smithall(x,1); }
+smith2(GEN x)
+{
+  GEN z = cgetg(4, t_VEC);
+  z[3] = (long)smithall(x, (GEN*)(z+1),(GEN*)(z+2));
+  return z;
+}
 
 /* Assume z was computed by [g]smithall(). Remove the 1s on the diagonal */
 GEN
@@ -3136,7 +3074,8 @@ smithclean(GEN z)
 static GEN
 gsmithall(GEN x,long all)
 {
-  long av,tetpil,i,j,k,l,c,n, lim;
+  gpmem_t av,tetpil,lim;
+  long i,j,k,l,c,n;
   GEN p1,p2,p3,p4,z,b,u,v,d,ml,mr;
 
   if (typ(x)!=t_MAT) err(typeer,"gsmithall");
@@ -3254,13 +3193,13 @@ gsmithall(GEN x,long all)
 GEN
 matsnf0(GEN x,long flag)
 {
-  long av = avma;
+  gpmem_t av = avma;
   if (flag > 7) err(flagerr,"matsnf");
   if (typ(x) == t_VEC && flag & 4) return smithclean(x);
-  if (flag & 2) x = gsmithall(x,flag & 1);
-  else          x = smithall(x, flag & 1);
-  if (flag & 4) x = smithclean(x);
-  return gerepileupto(av, x);
+  if (flag & 2) x = flag&1 ? gsmith2(x): gsmith(x);
+  else          x = flag&1 ?  smith2(x):  smith(x);
+  if (flag & 4) x = gerepileupto(av, smithclean(x));
+  return x;
 }
 
 GEN