===================================================================
RCS file: /home/cvs/OpenXM_contrib/pari-2.2/src/basemath/Attic/alglin1.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/alglin1.c	2001/10/02 11:17:00	1.1.1.1
+++ OpenXM_contrib/pari-2.2/src/basemath/Attic/alglin1.c	2002/09/11 07:26:48	1.2
@@ -1,4 +1,4 @@
-/* $Id: alglin1.c,v 1.1.1.1 2001/10/02 11:17:00 noro Exp $
+/* $Id: alglin1.c,v 1.2 2002/09/11 07:26:48 noro Exp $
 
 Copyright (C) 2000  The PARI group.
 
@@ -23,9 +23,9 @@ Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 
 /* for linear algebra mod p */
 #ifdef LONG_IS_64BIT
-#  define MASK (0x7fffffff00000000UL)
+#  define MASK (0xffffffff00000000UL)
 #else
-#  define MASK (0x7fff0000UL)
+#  define MASK (0xffff0000UL)
 #endif
 
 /* 2p^2 < 2^BITS_IN_LONG */
@@ -36,7 +36,7 @@ Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 #endif
 #define OK_ULONG(p) (lgefint(p) == 3 && u_OK_ULONG(p[2]))
 extern GEN ZM_init_CRT(GEN Hp, ulong p);
-extern ulong u_invmod(ulong x, ulong p);
+extern int ZM_incremental_CRT(GEN H, GEN Hp, GEN q, GEN qp, ulong p);
 
 /*******************************************************************/
 /*                                                                 */
@@ -75,7 +75,6 @@ gtrans_i(GEN x)
   return y;
 }
 
-
 GEN
 gtrans(GEN x)
 {
@@ -135,7 +134,8 @@ GEN
 concatsp3(GEN x, GEN y, GEN z)
 {
   long i, lx = lg(x), ly = lg(y), lz = lg(z);
-  GEN r = cgetg(lx+ly+lz-2, t_MAT), t = r;
+  GEN t, r = cgetg(lx+ly+lz-2, t_MAT);
+  t = r;
   for (i=1; i<lx; i++) *++t = *++x;
   for (i=1; i<ly; i++) *++t = *++y;
   for (i=1; i<lz; i++) *++t = *++z;
@@ -149,6 +149,8 @@ vconcat(GEN A, GEN B)
   long la,ha,hb,hc,i,j;
   GEN M,a,b,c;
 
+  if (!A) return B;
+  if (!B) return A;
   la = lg(A); if (la==1) return A;
   ha = lg(A[1]); M = cgetg(la,t_MAT);
   hb = lg(B[1]); hc = ha+hb-1;
@@ -162,11 +164,26 @@ vconcat(GEN A, GEN B)
 }
 
 GEN
+_veccopy(GEN x) { GEN v = cgetg(2, t_VEC); v[1] = lcopy(x); return v; }
+GEN
 _vec(GEN x) { GEN v = cgetg(2, t_VEC); v[1] = (long)x; return v; }
 GEN
 _col(GEN x) { GEN v = cgetg(2, t_COL); v[1] = (long)x; return v; }
 
+/* routines for naive growarrays */
 GEN
+cget1(long l, long t)
+{
+  GEN z = new_chunk(l);
+  z[0] = evaltyp(t) | evallg(1); return z;
+}
+void
+appendL(GEN x, GEN t)
+{
+  long l = lg(x); x[l] = (long)t; setlg(x, l+1);
+}
+
+GEN
 concatsp(GEN x, GEN y)
 {
   long tx=typ(x),ty=typ(y),lx=lg(x),ly=lg(y),i;
@@ -283,7 +300,7 @@ concat(GEN x, GEN y)
 
   if (!y)
   {
-    ulong av = avma;
+    gpmem_t av = avma;
     if (tx == t_LIST)
       { lx = lgef(x); i = 2; }
     else if (tx == t_VEC)
@@ -420,7 +437,7 @@ get_range(char *s, long *a, long *b, long *cmpl, long 
 
   *a = 1; *b = max;
   if (*s == '^') { *cmpl = 1; s++; } else *cmpl = 0;
-  if (*s == 0) return 0;
+  if (!*s) return 0;
   if (*s != '.')
   {
     *a = str_to_long(s, &s);
@@ -461,7 +478,27 @@ rowextract_i(GEN A, long x1, long x2)
   return B;
 }
 
+/* A[x0,] */
 GEN
+row(GEN A, long x0)
+{
+  long i, lB = lg(A);
+  GEN B  = cgetg(lB, t_VEC);
+  for (i=1; i<lB; i++) B[i] = coeff(A, x0, i);
+  return B;
+}
+
+/* A[x0, x1..x2] */
+GEN
+row_i(GEN A, long x0, long x1, long x2)
+{
+  long i, lB = x2 - x1 + 2;
+  GEN B  = cgetg(lB, t_VEC);
+  for (i=x1; i<=x2; i++) B[i] = coeff(A, x0, i);
+  return B;
+}
+
+GEN
 vecextract_p(GEN A, GEN p)
 {
   long i,lB = lg(p);
@@ -499,7 +536,8 @@ rowselect_p(GEN A, GEN B, GEN p, long init)
 GEN
 extract(GEN x, GEN l)
 {
-  long av,i,j, tl = typ(l), tx = typ(x), lx = lg(x);
+  gpmem_t av;
+  long i,j, tl = typ(l), tx = typ(x), lx = lg(x);
   GEN y;
 
   if (! is_matvec_t(tx)) err(typeer,"extract");
@@ -585,7 +623,7 @@ extract(GEN x, GEN l)
 GEN
 matextract(GEN x, GEN l1, GEN l2)
 {
-  long av = avma, tetpil;
+  gpmem_t av = avma, tetpil;
 
   if (typ(x)!=t_MAT) err(typeer,"matextract");
   x = extract(gtrans(extract(x,l2)),l1); tetpil=avma;
@@ -599,6 +637,19 @@ extract0(GEN x, GEN l1, GEN l2)
   return matextract(x,l1,l2);
 }
 
+/* v[a] + ... + v[b] */
+GEN
+sum(GEN v, long a, long b)
+{
+  GEN p;
+  long i;
+  if (a > b) return gzero;
+  if (b > lg(v)-1) err(talker,"incorrect length in sum");
+  p = (GEN)v[a];
+  for (i=a+1; i<=b; i++) p = gadd(p, (GEN)v[i]);
+  return p;
+}
+
 /*******************************************************************/
 /*                                                                 */
 /*                     SCALAR-MATRIX OPERATIONS                    */
@@ -668,15 +719,16 @@ GEN
 zeromat(long m, long n)
 {
   GEN y = cgetg(n+1,t_MAT);
-  long i; for (i=1; i<=n; i++) y[i]=(long)zerocol(m);
+  GEN v = zerocol(m);
+  long i; for (i=1; i<=n; i++) y[i]=(long)v;
   return y;
 }
 
 GEN
-gscalcol(GEN x, long n) 
-{ 
-  GEN y=gscalcol_proto(gzero,gzero,n); 
-  if (n) y[1]=lcopy(x); 
+gscalcol(GEN x, long n)
+{
+  GEN y=gscalcol_proto(gzero,gzero,n);
+  if (n) y[1]=lcopy(x);
   return y;
 }
 
@@ -847,22 +899,42 @@ mattodiagonal(GEN m)
 GEN
 gaddmat(GEN x, GEN y)
 {
-  long ly,dy,i,j;
-  GEN z;
+  long l,d,i,j;
+  GEN z, cz,cy;
 
-  ly=lg(y); if (ly==1) return cgetg(1,t_MAT);
-  dy=lg(y[1]);
-  if (typ(y)!=t_MAT || ly!=dy) err(mattype1,"gaddmat");
-  z=cgetg(ly,t_MAT);
-  for (i=1; i<ly; i++)
+  l = lg(y); if (l==1) return cgetg(1,t_MAT);
+  d = lg(y[1]);
+  if (typ(y)!=t_MAT || l!=d) err(mattype1,"gaddmat");
+  z=cgetg(l,t_MAT);
+  for (i=1; i<l; i++)
   {
-    z[i]=lgetg(dy,t_COL);
-    for (j=1; j<dy; j++)
-      coeff(z,j,i) = i==j? ladd(x,gcoeff(y,j,i)): lcopy(gcoeff(y,j,i));
+    cz = cgetg(d,t_COL); z[i] = (long)cz; cy = (GEN)y[i];
+    for (j=1; j<d; j++)
+      cz[j] = i==j? ladd(x,(GEN)cy[j]): lcopy((GEN)cy[j]);
   }
   return z;
 }
 
+/* same, no copy */
+GEN
+gaddmat_i(GEN x, GEN y)
+{
+  long l,d,i,j;
+  GEN z, cz,cy;
+
+  l = lg(y); if (l==1) return cgetg(1,t_MAT);
+  d = lg(y[1]);
+  if (typ(y)!=t_MAT || l!=d) err(mattype1,"gaddmat");
+  z = cgetg(l,t_MAT);
+  for (i=1; i<l; i++)
+  {
+    cz = cgetg(d,t_COL); z[i] = (long)cz; cy = (GEN)y[i];
+    for (j=1; j<d; j++)
+      cz[j] = i==j? ladd(x,(GEN)cy[j]): cy[j];
+  }
+  return z;
+}
+
 /*******************************************************************/
 /*                                                                 */
 /*                       Solve A*X=B (Gauss pivot)                 */
@@ -871,47 +943,45 @@ gaddmat(GEN x, GEN y)
 #define swap(x,y) { long _t=x; x=y; y=_t; }
 
 /* Assume x is a non-empty matrix. Return 0 if maximal pivot should not be
- * used, and the matrix precision (min real precision of coeffs) otherwise.
- */
-static long
-matprec(GEN x)
+ * used, 1 otherwise */
+static int
+use_maximal_pivot(GEN x)
 {
-  long tx,i,j,l, k = VERYBIGINT, lx = lg(x), ly = lg(x[1]);
+  int res = 0;
+  long tx,i,j, lx = lg(x), ly = lg(x[1]);
   GEN p1;
   for (i=1; i<lx; i++)
     for (j=1; j<ly; j++)
     {
       p1 = gmael(x,i,j); tx = typ(p1);
       if (!is_scalar_t(tx)) return 0;
-      l = precision(p1); if (l && l<k) k = l;
+      if (precision(p1)) res = 1;
     }
-  return (k==VERYBIGINT)? 0: k;
+  return res;
 }
 
-/* As above, returning 1 if the precision would be non-zero, 0 otherwise */
-static long
-use_maximal_pivot(GEN x)
+static GEN
+col_to_mat(GEN b)
 {
-  long tx,i,j, lx = lg(x), ly = lg(x[1]);
-  GEN p1;
-  for (i=1; i<lx; i++)
-    for (j=1; j<ly; j++)
-    {
-      p1 = gmael(x,i,j); tx = typ(p1);
-      if (!is_scalar_t(tx)) return 0;
-      if (precision(p1)) return 1;
-    }
-  return 0;
+  GEN B = cgetg(2, t_MAT);
+  B[1] = (long)b; return B;
 }
 
 static GEN
-check_b(GEN b, long nbli)
+check_b(GEN b, long nbli, int *iscol)
 {
   GEN col;
+  *iscol = (b && typ(b) == t_COL);
   if (!b) return idmat(nbli);
-  col = (typ(b) == t_MAT)? (GEN)b[1]: b;
+  b = dummycopy(b);
+  if (*iscol) { col = b; b = col_to_mat(b); }
+  else
+  {
+    if (lg(b) == 1) err(consister,"gauss");
+    col = (GEN)b[1];
+  }
   if (nbli != lg(col)-1) err(talker,"incompatible matrix dimensions in gauss");
-  return dummycopy(b);
+  return b;
 }
 
 /* C = A^(-1)(tB) where A, B, C are integral, A is upper triangular, t t_INT */
@@ -925,21 +995,77 @@ gauss_triangle_i(GEN A, GEN B, GEN t)
   for (k=1; k<=n; k++)
   {
     GEN u = cgetg(n+1, t_COL), b = (GEN)B[k];
-    ulong av = avma;
+    gpmem_t av = avma;
     c[k] = (long)u; m = mulii((GEN)b[n],t);
-    u[n] = (long)gerepileuptoint(av, divii(m, gcoeff(A,n,n)));
+    u[n] = lpileuptoint(av, divii(m, gcoeff(A,n,n)));
     for (i=n-1; i>0; i--)
     {
       av = avma; m = mulii(negi((GEN)b[i]),t);
       for (j=i+1; j<=n; j++)
         m = addii(m, mulii(gcoeff(A,i,j),(GEN) u[j]));
-      u[i] = (long)gerepileuptoint(av, divii(negi(m), gcoeff(A,i,i)));
+      u[i] = lpileuptoint(av, divii(negi(m), gcoeff(A,i,i)));
     }
   }
   return c;
 }
 
+/* A, B integral upper HNF. A^(-1) B integral ? */
+int
+hnfdivide(GEN A, GEN B)
+{
+  gpmem_t av = avma;
+  long n = lg(A)-1, i,j,k;
+  GEN u, b, m, r;
+
+  if (!n) return 1;
+  u = cgetg(n+1, t_COL);
+  for (k=1; k<=n; k++)
+  {
+    b = (GEN)B[k];
+    m = (GEN)b[k];
+    u[k] = ldvmdii(m, gcoeff(A,k,k), &r);
+    if (r != gzero) { avma = av; return 0; }
+    for (i=k-1; i>0; i--)
+    {
+      m = negi((GEN)b[i]);
+      for (j=i+1; j<=k; j++)
+        m = addii(m, mulii(gcoeff(A,i,j),(GEN) u[j]));
+      m = dvmdii(m, gcoeff(A,i,i), &r);
+      if (r != gzero) { avma = av; return 0; }
+      u[i] = lnegi(m);
+    }
+  }
+  avma = av; return 1;
+}
+
+/* A upper HNF, b integral vector. Return A^(-1) b if integral,
+ * NULL otherwise. */
 GEN
+hnf_invimage(GEN A, GEN b)
+{
+  gpmem_t av = avma, av2;
+  long n = lg(A)-1, i,j;
+  GEN u, m, r;
+
+  if (!n) return NULL;
+  u = cgetg(n+1, t_COL);
+  m = (GEN)b[n];
+  u[n] = ldvmdii(m, gcoeff(A,n,n), &r);
+  if (r != gzero) { avma = av; return NULL; }
+  for (i=n-1; i>0; i--)
+  {
+    av2 = avma;
+    m = negi((GEN)b[i]);
+    for (j=i+1; j<=n; j++)
+      m = addii(m, mulii(gcoeff(A,i,j),(GEN) u[j]));
+    m = dvmdii(m, gcoeff(A,i,i), &r);
+    if (r != gzero) { avma = av; return NULL; }
+    u[i] = lpileuptoint(av2, negi(m));
+  }
+  return u;
+}
+
+GEN
 gauss_get_col(GEN a, GEN b, GEN p, long li)
 {
   GEN m, u=cgetg(li+1,t_COL);
@@ -948,10 +1074,11 @@ gauss_get_col(GEN a, GEN b, GEN p, long li)
   u[li] = ldiv((GEN)b[li], p);
   for (i=li-1; i>0; i--)
   {
+    gpmem_t av = avma;
     m = gneg_i((GEN)b[i]);
     for (j=i+1; j<=li; j++)
       m = gadd(m, gmul(gcoeff(a,i,j), (GEN)u[j]));
-    u[i] = ldiv(gneg_i(m), gcoeff(a,i,i));
+    u[i] = lpileupto(av, gdiv(gneg_i(m), gcoeff(a,i,i)));
   }
   return u;
 }
@@ -965,38 +1092,79 @@ Fp_gauss_get_col(GEN a, GEN b, GEN piv, long li, GEN p
   u[li] = lresii(mulii((GEN)b[li], mpinvmod(piv,p)), p);
   for (i=li-1; i>0; i--)
   {
+    gpmem_t av = avma;
     m = (GEN)b[i];
     for (j=i+1; j<=li; j++)
       m = subii(m, mulii(gcoeff(a,i,j), (GEN)u[j]));
     m = resii(m, p);
-    u[i] = lresii(mulii(m, mpinvmod(gcoeff(a,i,i), p)), p);
+    u[i] = lpileuptoint(av, resii(mulii(m, mpinvmod(gcoeff(a,i,i), p)), p));
   }
   return u;
 }
+GEN
+Fq_gauss_get_col(GEN a, GEN b, GEN piv, long li, GEN T, GEN p)
+{
+  GEN m, u=cgetg(li+1,t_COL);
+  long i,j;
 
+  u[li] = (long)FpXQ_mul((GEN)b[li], FpXQ_inv(piv,T,p), T,p);
+  for (i=li-1; i>0; i--)
+  {
+    gpmem_t av = avma;
+    m = (GEN)b[i];
+    for (j=i+1; j<=li; j++)
+      m = gsub(m, gmul(gcoeff(a,i,j), (GEN)u[j]));
+    m = FpX_res(m, T,p);
+    u[i] = lpileupto(av, FpXQ_mul(m, FpXQ_inv(gcoeff(a,i,i), T,p), T,p));
+  }
+  return u;
+}
+
 /* assume -p < a < p, return 1/a mod p */
 static long
 u_Fp_inv(long a, long p)
 {
   if (a < 0) a = p + a; /* pb with ulongs < 0 */
-  return u_invmod(a,p);
+  return (long)u_invmod((ulong)a,(ulong)p);
 }
 
+/* assume 0 <= a[i,i] < p */
 GEN
-u_Fp_gauss_get_col(GEN a, GEN b, long piv, long li, long p)
+u_Fp_gauss_get_col(GEN a, GEN b, ulong piv, long li, ulong p)
 {
-  GEN u=cgetg(li+1,t_VECSMALL);
-  long i,j, m;
+  GEN u = cgetg(li+1,t_VECSMALL);
+  ulong m;
+  long i,j;
 
-  u[li] = (b[li] *  u_Fp_inv(piv,p)) % p;
-  if (u[li] < 0) u[li] += p;
-  for (i=li-1; i>0; i--)
+  m = b[li] % p;
+  if (u_OK_ULONG(p))
   {
-    m = b[i];
-    for (j=i+1; j<=li; j++) { m -= coeff(a,i,j) * u[j]; if (m & MASK) m %= p; }
-    u[i] = ((m%p) * u_Fp_inv(coeff(a,i,i), p)) % p;
-    if (u[i] < 0) u[i] += p;
+    u[li] = (m * u_Fp_inv(piv,p)) % p;
+    for (i=li-1; i>0; i--)
+    {
+      m = p - b[i]%p;
+      for (j=i+1; j<=li; j++) {
+        if (m & HIGHBIT) m %= p;
+        m += coeff(a,i,j) * u[j];
+      }
+      m %= p;
+      if (m) m = ((p-m) * u_invmod((ulong)coeff(a,i,i), p)) % p; 
+      u[i] = m;
+    }
   }
+  else
+  {
+    u[li] = mulssmod(m, u_Fp_inv(piv,p), p);
+    for (i=li-1; i>0; i--)
+    {
+      m = p - b[i]%p;
+      for (j=i+1; j<=li; j++)
+        m += mulssmod(coeff(a,i,j), u[j], p);
+      m %= p;
+      if (m) m = mulssmod(p-m, u_invmod((ulong)coeff(a,i,i), p), p);
+      u[i] = m;
+    }
+  }
   return u;
 }
 
@@ -1015,16 +1183,37 @@ _Fp_addmul(GEN b, long k, long i, GEN m, GEN p)
   b[k] = laddii((GEN)b[k], mulii(m, (GEN)b[i]));
 }
 
+/* same, reduce mod (T,p) */
 static void
-_u_Fp_addmul(GEN b, long k, long i, long m, long p)
+_Fq_addmul(GEN b, long k, long i, GEN m, GEN T, GEN p)
 {
-  long a;
-  if (b[i] & MASK) b[i] %= p;
-  a = b[k] + m * b[i];
-  if (a & MASK) a %= p;
-  b[k] = a;
+  b[i] = (long)FpX_res((GEN)b[i], T,p);
+  b[k] = (long)ladd((GEN)b[k], gmul(m, (GEN)b[i]));
 }
 
+/* assume m < p && u_OK_ULONG(p) && (! (b[i] & MASK)) */
+static void
+_u_Fp_addmul_OK(GEN b, long k, long i, ulong m, ulong p)
+{
+  b[k] += m * b[i];
+  if (b[k] & MASK) b[k] %= p;
+}
+/* assume m < p */
+static void
+_u_Fp_addmul(GEN b, long k, long i, ulong m, ulong p)
+{
+  b[i] %= p;
+  b[k] += mulssmod(m, b[i], p);
+  if (b[k] & MASK) b[k] %= p;
+}
+/* same m = 1 */
+static void
+_u_Fp_add(GEN b, long k, long i, ulong p)
+{
+  b[k] += b[i];
+  if (b[k] & MASK) b[k] %= p;
+}
+
 /* Gaussan Elimination. Compute a^(-1)*b
  * b is a matrix or column vector, NULL meaning: take the identity matrix
  * If a and b are empty, the result is the empty matrix.
@@ -1037,7 +1226,9 @@ _u_Fp_addmul(GEN b, long k, long i, long m, long p)
 GEN
 gauss_intern(GEN a, GEN b)
 {
-  long inexact,iscol,i,j,k,av,lim,li,bco, aco = lg(a)-1;
+  gpmem_t av, lim;
+  long i,j,k,li,bco, aco = lg(a)-1;
+  int inexact, iscol;
   GEN p,m,u;
 
   if (typ(a)!=t_MAT) err(mattype1,"gauss");
@@ -1048,15 +1239,13 @@ gauss_intern(GEN a, GEN b)
     if (DEBUGLEVEL) err(warner,"in Gauss lg(a)=1 lg(b)=%ld", b?1:-1);
     return cgetg(1,t_MAT);
   }
-  av=avma; lim=stack_lim(av,1);
+  av = avma; lim = stack_lim(av,1);
   li = lg(a[1])-1;
   if (li != aco && (li < aco || b)) err(mattype1,"gauss");
   a = dummycopy(a);
-  b = check_b(b,li); bco = lg(b)-1;
+  b = check_b(b,li, &iscol); bco = lg(b)-1;
   inexact = use_maximal_pivot(a);
-  iscol   = (typ(b)==t_COL);
-  if(DEBUGLEVEL>4)
-    fprintferr("Entering gauss with inexact=%ld iscol=%ld\n",inexact,iscol);
+  if(DEBUGLEVEL>4) fprintferr("Entering gauss with inexact=%ld\n",inexact);
 
   p = NULL; /* gcc -Wall */
   for (i=1; i<=aco; i++)
@@ -1083,51 +1272,32 @@ gauss_intern(GEN a, GEN b)
     if (k != i)
     {
       for (j=i; j<=aco; j++) swap(coeff(a,i,j), coeff(a,k,j));
-      if (iscol) { swap(b[i],b[k]); }
-      else
-        for (j=1; j<=bco; j++) swap(coeff(b,i,j), coeff(b,k,j));
+      for (j=1; j<=bco; j++) swap(coeff(b,i,j), coeff(b,k,j));
       p = gcoeff(a,i,i);
     }
     if (i == aco) break;
 
     for (k=i+1; k<=li; k++)
     {
-      m=gcoeff(a,k,i);
+      m = gcoeff(a,k,i);
       if (!gcmp0(m))
       {
 	m = gneg_i(gdiv(m,p));
 	for (j=i+1; j<=aco; j++) _addmul((GEN)a[j],k,i,m);
-	if (iscol) _addmul(b,k,i,m);
-        else 
-          for (j=1; j<=bco; j++) _addmul((GEN)b[j],k,i,m); 
+        for (j=1;   j<=bco; j++) _addmul((GEN)b[j],k,i,m);
       }
     }
     if (low_stack(lim, stack_lim(av,1)))
     {
-      GEN *gptr[2]; gptr[0]=&a; gptr[1]=&b;
       if(DEBUGMEM>1) err(warnmem,"gauss. i=%ld",i);
-      gerepilemany(av,gptr,2);
+      gerepileall(av,2, &a,&b);
     }
   }
 
   if(DEBUGLEVEL>4) fprintferr("Solving the triangular system\n");
-  if (iscol) u = gauss_get_col(a,b,p,aco);
-  else
-  {
-    long av1 = avma;
-    lim = stack_lim(av1,1); u=cgetg(bco+1,t_MAT);
-    for (j=2; j<=bco; j++) u[j] = zero; /* dummy */
-    for (j=1; j<=bco; j++)
-    {
-      u[j] = (long)gauss_get_col(a,(GEN)b[j],p,aco);
-      if (low_stack(lim, stack_lim(av1,1)))
-      {
-        if(DEBUGMEM>1) err(warnmem,"gauss[2]. j=%ld", j);
-        u = gerepilecopy(av1, u);
-      }
-    }
-  }
-  return gerepilecopy(av, u);
+  u = cgetg(bco+1,t_MAT);
+  for (j=1; j<=bco; j++) u[j] = (long)gauss_get_col(a,(GEN)b[j],p,aco);
+  return gerepilecopy(av, iscol? (GEN)u[1]: u);
 }
 
 GEN
@@ -1138,29 +1308,32 @@ gauss(GEN a, GEN b)
   return z;
 }
 
-GEN
-u_FpM_gauss(GEN a, GEN b, long p)
+/* destroy a, b */
+static GEN
+u_FpM_gauss_sp(GEN a, GEN b, ulong p)
 {
-  long piv,m,iscol,i,j,k,li,bco, aco = lg(a)-1;
+  long iscol,i,j,k,li,bco, aco = lg(a)-1;
+  ulong piv, m;
   GEN u;
 
   if (!aco) return cgetg(1,t_MAT);
   li = lg(a[1])-1;
   bco = lg(b)-1;
-  iscol = (typ(b)==t_COL);
+  iscol = (typ(b)!=t_MAT);
+  if (iscol) b = col_to_mat(b);
   piv = 0; /* gcc -Wall */
   for (i=1; i<=aco; i++)
   {
     /* k is the line where we find the pivot */
-    coeff(a,i,i) = coeff(a,i,i) % p;
-    piv = coeff(a,i,i); k = i;
+    piv = coeff(a,i,i) % p;
+    coeff(a,i,i) = piv; k = i;
     if (!piv)
     {
       for (k++; k <= li; k++)
       {
         coeff(a,k,i) %= p;
         if (coeff(a,k,i)) break;
-      } 
+      }
       if (k>li) return NULL;
     }
 
@@ -1168,41 +1341,65 @@ u_FpM_gauss(GEN a, GEN b, long p)
     if (k != i)
     {
       for (j=i; j<=aco; j++) swap(coeff(a,i,j), coeff(a,k,j));
-      if (iscol) { swap(b[i],b[k]); }
-      else
-        for (j=1; j<=bco; j++) swap(coeff(b,i,j), coeff(b,k,j));
+      for (j=1; j<=bco; j++) swap(coeff(b,i,j), coeff(b,k,j));
       piv = coeff(a,i,i);
     }
     if (i == aco) break;
 
     for (k=i+1; k<=li; k++)
     {
-      coeff(a,k,i) %= p;
-      m = coeff(a,k,i); coeff(a,k,i) = 0;
+      m = coeff(a,k,i) % p; coeff(a,k,i) = 0;
       if (m)
       {
-	m = - (m * u_Fp_inv(piv,p)) % p;
-	for (j=i+1; j<=aco; j++) _u_Fp_addmul((GEN)a[j],k,i,m, p);
-	if (iscol) _u_Fp_addmul(b,k,i,m, p);
-        else 
-          for (j=1; j<=bco; j++) _u_Fp_addmul((GEN)b[j],k,i,m, p); 
+	m = mulssmod(m, u_invmod(piv,p), p);
+        m = p - m;
+        if (u_OK_ULONG(p))
+        {
+          for (j=i+1; j<=aco; j++) _u_Fp_addmul_OK((GEN)a[j],k,i,m, p);
+          for (j=1;   j<=bco; j++) _u_Fp_addmul_OK((GEN)b[j],k,i,m, p);
+        }
+        else
+        {
+          for (j=i+1; j<=aco; j++) _u_Fp_addmul((GEN)a[j],k,i,m, p);
+          for (j=1;   j<=bco; j++) _u_Fp_addmul((GEN)b[j],k,i,m, p);
+        }
       }
     }
   }
-  if (iscol) u = u_Fp_gauss_get_col(a,b,piv,aco,p);
-  else
-  {
-    u=cgetg(bco+1,t_MAT);
-    for (j=1; j<=bco; j++)
-      u[j] = (long)u_Fp_gauss_get_col(a,(GEN)b[j],piv,aco,p);
-  }
-  return u;
+  u = cgetg(bco+1,t_MAT);
+  for (j=1; j<=bco; j++) u[j] = (long)u_Fp_gauss_get_col(a,(GEN)b[j],piv,aco,p);
+  return iscol? (GEN)u[1]: u;
 }
 
+static GEN
+u_idmat(long n)
+{
+  GEN y = cgetg(n+1,t_MAT);
+  long i;
+  if (n < 0) err(talker,"negative size in u_idmat");
+  for (i=1; i<=n; i++) { y[i] = (long)vecsmall_const(n, 0); coeff(y, i,i) = 1; }
+  return y;
+}
+
 GEN
+u_FpM_gauss(GEN a, GEN b, ulong p) {
+  return u_FpM_gauss_sp(dummycopy(a), dummycopy(b), p);
+}
+static GEN
+u_FpM_inv_sp(GEN a, ulong p) {
+  return u_FpM_gauss_sp(a, u_idmat(lg(a)-1), p);
+}
+GEN
+u_FpM_inv(GEN a, ulong p) {
+  return u_FpM_inv_sp(dummycopy(a), p);
+}
+
+GEN
 FpM_gauss(GEN a, GEN b, GEN p)
 {
-  long iscol,i,j,k,av,lim,li,bco, aco = lg(a)-1;
+  gpmem_t av,lim;
+  long i,j,k,li,bco, aco = lg(a)-1;
+  int iscol;
   GEN piv,m,u;
 
   if (typ(a)!=t_MAT) err(mattype1,"gauss");
@@ -1215,19 +1412,19 @@ FpM_gauss(GEN a, GEN b, GEN p)
   }
   li = lg(a[1])-1;
   if (li != aco && (li < aco || b)) err(mattype1,"gauss");
-  b = check_b(b,li); av = avma;
+  b = check_b(b,li, &iscol); av = avma;
   if (OK_ULONG(p))
   {
     ulong pp=p[2];
     a = u_Fp_FpM(a, pp);
     b = u_Fp_FpM(b, pp);
-    u = u_FpM_gauss(a,b, pp);
-    return gerepileupto(av, small_to_mat(u));
+    u = u_FpM_gauss_sp(a,b, pp);
+    u = iscol? small_to_col((GEN)u[1]): small_to_mat(u);
+    return gerepileupto(av, u);
   }
   lim = stack_lim(av,1);
   a = dummycopy(a);
   bco = lg(b)-1;
-  iscol = (typ(b)==t_COL);
   piv = NULL; /* gcc -Wall */
   for (i=1; i<=aco; i++)
   {
@@ -1240,7 +1437,7 @@ FpM_gauss(GEN a, GEN b, GEN p)
       {
         coeff(a,k,i) = lresii(gcoeff(a,k,i), p);
         if (signe(coeff(a,k,i))) break;
-      } 
+      }
       if (k>li) return NULL;
     }
 
@@ -1248,9 +1445,7 @@ FpM_gauss(GEN a, GEN b, GEN p)
     if (k != i)
     {
       for (j=i; j<=aco; j++) swap(coeff(a,i,j), coeff(a,k,j));
-      if (iscol) { swap(b[i],b[k]); }
-      else
-        for (j=1; j<=bco; j++) swap(coeff(b,i,j), coeff(b,k,j));
+      for (j=1; j<=bco; j++) swap(coeff(b,i,j), coeff(b,k,j));
       piv = gcoeff(a,i,i);
     }
     if (i == aco) break;
@@ -1264,9 +1459,7 @@ FpM_gauss(GEN a, GEN b, GEN p)
 	m = mulii(m, mpinvmod(piv,p));
         m = resii(negi(m), p);
 	for (j=i+1; j<=aco; j++) _Fp_addmul((GEN)a[j],k,i,m, p);
-	if (iscol) _Fp_addmul(b,k,i,m, p);
-        else 
-          for (j=1; j<=bco; j++) _Fp_addmul((GEN)b[j],k,i,m, p); 
+        for (j=1  ; j<=bco; j++) _Fp_addmul((GEN)b[j],k,i,m, p);
       }
     }
     if (low_stack(lim, stack_lim(av,1)))
@@ -1278,42 +1471,89 @@ FpM_gauss(GEN a, GEN b, GEN p)
   }
 
   if(DEBUGLEVEL>4) fprintferr("Solving the triangular system\n");
-  if (iscol) u = Fp_gauss_get_col(a,b,piv,aco,p);
-  else
+  u = cgetg(bco+1,t_MAT);
+  for (j=1; j<=bco; j++) u[j] = (long)Fp_gauss_get_col(a,(GEN)b[j],piv,aco,p);
+  return gerepilecopy(av, iscol? (GEN)u[1]: u);
+}
+GEN
+FqM_gauss(GEN a, GEN b, GEN T, GEN p)
+{
+  gpmem_t av,lim;
+  long i,j,k,li,bco, aco = lg(a)-1;
+  int iscol;
+  GEN piv,m,u;
+
+  if (!T) return FpM_gauss(a,b,p);
+
+  if (typ(a)!=t_MAT) err(mattype1,"gauss");
+  if (b && typ(b)!=t_COL && typ(b)!=t_MAT) err(typeer,"gauss");
+  if (!aco)
   {
-    long av1 = avma;
-    lim = stack_lim(av1,1); u=cgetg(bco+1,t_MAT);
-    for (j=2; j<=bco; j++) u[j] = zero; /* dummy */
-    for (j=1; j<=bco; j++)
+    if (b && lg(b)!=1) err(consister,"gauss");
+    if (DEBUGLEVEL) err(warner,"in Gauss lg(a)=1 lg(b)=%ld", b?1:-1);
+    return cgetg(1,t_MAT);
+  }
+  li = lg(a[1])-1;
+  if (li != aco && (li < aco || b)) err(mattype1,"gauss");
+  b = check_b(b,li,&iscol); av = avma;
+  lim = stack_lim(av,1);
+  a = dummycopy(a);
+  bco = lg(b)-1;
+  piv = NULL; /* gcc -Wall */
+  for (i=1; i<=aco; i++)
+  {
+    /* k is the line where we find the pivot */
+    coeff(a,i,i) = (long)FpX_res(gcoeff(a,i,i), T,p);
+    piv = gcoeff(a,i,i); k = i;
+    if (!signe(piv))
     {
-      u[j] = (long)Fp_gauss_get_col(a,(GEN)b[j],piv,aco,p);
-      if (low_stack(lim, stack_lim(av1,1)))
+      for (k++; k <= li; k++)
       {
-        if(DEBUGMEM>1) err(warnmem,"gauss[2]. j=%ld", j);
-        u = gerepilecopy(av1, u);
+        coeff(a,k,i) = (long)FpX_res(gcoeff(a,k,i), T,p);
+        if (signe(coeff(a,k,i))) break;
       }
+      if (k>li) return NULL;
     }
+
+    /* if (k!=i), exchange the lines s.t. k = i */
+    if (k != i)
+    {
+      for (j=i; j<=aco; j++) swap(coeff(a,i,j), coeff(a,k,j));
+      for (j=1; j<=bco; j++) swap(coeff(b,i,j), coeff(b,k,j));
+      piv = gcoeff(a,i,i);
+    }
+    if (i == aco) break;
+
+    for (k=i+1; k<=li; k++)
+    {
+      coeff(a,k,i) = (long)FpX_res(gcoeff(a,k,i), T,p);
+      m = gcoeff(a,k,i); coeff(a,k,i) = zero;
+      if (signe(m))
+      {
+	m = FpXQ_mul(m, FpXQ_inv(piv,T,p),T,p);
+        m = resii(negi(m), p);
+	for (j=i+1; j<=aco; j++) _Fq_addmul((GEN)a[j],k,i,m, T,p);
+        for (j=1;   j<=bco; j++) _Fq_addmul((GEN)b[j],k,i,m, T,p);
+      }
+    }
+    if (low_stack(lim, stack_lim(av,1)))
+    {
+      GEN *gptr[2]; gptr[0]=&a; gptr[1]=&b;
+      if(DEBUGMEM>1) err(warnmem,"FpM_gauss. i=%ld",i);
+      gerepilemany(av,gptr,2);
+    }
   }
-  return gerepilecopy(av, u);
+
+  if(DEBUGLEVEL>4) fprintferr("Solving the triangular system\n");
+  u = cgetg(bco+1,t_MAT);
+  for (j=1; j<=bco; j++) u[j] = (long)Fq_gauss_get_col(a,(GEN)b[j],piv,aco,T,p);
+  return gerepilecopy(av, iscol? (GEN)u[1]: u);
 }
 
+
 GEN
 FpM_inv(GEN x, GEN p) { return FpM_gauss(x, NULL, p); }
 
-static GEN
-u_idmat(long n)
-{
-  GEN y = cgetg(n+1,t_MAT);
-  long i,j;
-  if (n < 0) err(talker,"negative size in u_idmat");
-  for (i=1; i<=n; i++)
-  {
-    y[i]=lgetg(n+1,t_VECSMALL);
-    for (j=1; j<=n; j++) coeff(y,j,i) = (i==j)? 1: 0;
-  }
-  return y;
-}
-
 /* set y = x * y */
 #define HIGHWORD(a) ((a) >> BITS_IN_HALFULONG)
 static GEN
@@ -1332,11 +1572,12 @@ u_FpM_Fp_mul_ip(GEN y, long x, long p)
 }
 
 /* M integral, dM such that M' = dM M^-1 is integral [e.g det(M)]. Return M' */
-GEN 
+GEN
 ZM_inv(GEN M, GEN dM)
 {
-  GEN  ID,Hp,q,H;
-  ulong p,dMp, av2, av = avma, lim = stack_lim(av,1);
+  gpmem_t av2, av = avma, lim = stack_lim(av,1);
+  GEN Hp,q,H;
+  ulong p,dMp;
   byteptr d = diffptr;
   long lM = lg(M), stable = 0;
 
@@ -1346,16 +1587,14 @@ ZM_inv(GEN M, GEN dM)
   av2 = avma;
   H = NULL;
   d += 3000; /* 27449 = prime(3000) */
-  for(p = 27449; ; p+= *d++)
+  for(p = 27449; ; )
   {
-    if (!*d) err(primer1);
     dMp = umodiu(dM,p);
-    if (!dMp) continue;
-    ID = u_idmat(lM-1);
-    Hp = u_FpM_gauss(u_Fp_FpM(M,p), ID, p);
+    if (!dMp) goto repeat;
+    Hp = u_FpM_inv_sp(u_Fp_FpM(M,p), p);
     if (dMp != 1) Hp = u_FpM_Fp_mul_ip(Hp, dMp, p);
 
-    if (!H) 
+    if (!H)
     {
       H = ZM_init_CRT(Hp, p);
       q = utoi(p);
@@ -1368,27 +1607,28 @@ ZM_inv(GEN M, GEN dM)
     }
     if (DEBUGLEVEL>5) msgtimer("inverse mod %ld (stable=%ld)", p,stable);
     if (stable && isscalarmat(gmul(M, H), dM)) break; /* DONE */
-    
+
     if (low_stack(lim, stack_lim(av,2)))
     {
       GEN *gptr[2]; gptr[0] = &H; gptr[1] = &q;
       if (DEBUGMEM>1) err(warnmem,"ZM_inv");
       gerepilemany(av2,gptr, 2);
     }
- 
+   repeat:
+    NEXT_PRIME_VIADIFF_CHECK(p,d);
   }
   if (DEBUGLEVEL>5) msgtimer("ZM_inv done");
   return gerepilecopy(av, H);
 }
 
 /* same as above, M rational */
-GEN 
+GEN
 QM_inv(GEN M, GEN dM)
 {
-  ulong av = avma;
-  GEN cM = content(M);
-  if (is_pm1(cM)) { avma = av; return ZM_inv(M,dM); }
-  return gerepileupto(av, ZM_inv(gdiv(M,cM), gdiv(dM,cM)));
+  gpmem_t av = avma;
+  GEN cM, pM = Q_primitive_part(M, &cM);
+  if (!cM) return ZM_inv(pM,dM);
+  return gerepileupto(av, ZM_inv(pM, gdiv(dM,cM)));
 }
 
 /* x a matrix with integer coefficients. Return a multiple of the determinant
@@ -1397,8 +1637,9 @@ QM_inv(GEN M, GEN dM)
 GEN
 detint(GEN x)
 {
+  gpmem_t av = avma, av1, lim;
   GEN pass,c,v,det1,piv,pivprec,vi,p1;
-  long i,j,k,rg,n,m,m1,av=avma,av1,lim;
+  long i,j,k,rg,n,m,m1;
 
   if (typ(x)!=t_MAT) err(typeer,"detint");
   n=lg(x)-1; if (!n) return gun;
@@ -1462,10 +1703,11 @@ detint(GEN x)
 }
 
 static void
-gerepile_gauss_ker(GEN x, long m, long n, long k, long t, ulong av)
+gerepile_gauss_ker(GEN x, long k, long t, gpmem_t av)
 {
-  ulong tetpil = avma,A;
-  long dec,u,i;
+  gpmem_t tetpil = avma, A;
+  long u,i, n = lg(x)-1, m = n? lg(x[1])-1: 0;
+  size_t dec;
 
   if (DEBUGMEM > 1) err(warnmem,"gauss_pivot_ker. k=%ld, n=%ld",k,n);
   for (u=t+1; u<=m; u++) copyifstack(coeff(x,u,k), coeff(x,u,k));
@@ -1475,22 +1717,23 @@ gerepile_gauss_ker(GEN x, long m, long n, long k, long
   (void)gerepile(av,tetpil,NULL); dec = av-tetpil;
   for (u=t+1; u<=m; u++)
   {
-    A=coeff(x,u,k);
+    A=(gpmem_t)coeff(x,u,k);
     if (A<av && A>=bot) coeff(x,u,k)+=dec;
   }
   for (i=k+1; i<=n; i++)
     for (u=1; u<=m; u++)
     {
-      A=coeff(x,u,i);
+      A=(gpmem_t)coeff(x,u,i);
       if (A<av && A>=bot) coeff(x,u,i)+=dec;
     }
 }
 
 static void
-gerepile_gauss_FpM_ker(GEN x, GEN p, long m, long n, long k, long t, ulong av)
+gerepile_gauss_FpM_ker(GEN x, GEN p, long k, long t, gpmem_t av)
 {
-  ulong tetpil = avma,A;
-  long dec,u,i;
+  gpmem_t tetpil = avma, A;
+  long u,i, n = lg(x)-1, m = n? lg(x[1])-1: 0;
+  size_t dec;
 
   if (DEBUGMEM > 1) err(warnmem,"gauss_pivot_ker. k=%ld, n=%ld",k,n);
   for (u=t+1; u<=m; u++)
@@ -1502,13 +1745,13 @@ gerepile_gauss_FpM_ker(GEN x, GEN p, long m, long n, l
   (void)gerepile(av,tetpil,NULL); dec = av-tetpil;
   for (u=t+1; u<=m; u++)
   {
-    A=coeff(x,u,k);
+    A=(gpmem_t)coeff(x,u,k);
     if (A<av && A>=bot) coeff(x,u,k)+=dec;
   }
   for (i=k+1; i<=n; i++)
     for (u=1; u<=m; u++)
     {
-      A=coeff(x,u,i);
+      A=(gpmem_t)coeff(x,u,i);
       if (A<av && A>=bot) coeff(x,u,i)+=dec;
     }
 }
@@ -1516,10 +1759,11 @@ gerepile_gauss_FpM_ker(GEN x, GEN p, long m, long n, l
 /* special gerepile for huge matrices */
 
 static void
-gerepile_gauss(GEN x,long m,long n,long k,long t,ulong av, long j, GEN c)
+gerepile_gauss(GEN x,long k,long t,gpmem_t av, long j, GEN c)
 {
-  ulong tetpil = avma,A;
-  long dec,u,i;
+  gpmem_t tetpil = avma, A;
+  long u,i, n = lg(x)-1, m = n? lg(x[1])-1: 0;
+  size_t dec;
 
   if (DEBUGMEM > 1) err(warnmem,"gauss_pivot. k=%ld, n=%ld",k,n);
   for (u=t+1; u<=m; u++)
@@ -1532,14 +1776,14 @@ gerepile_gauss(GEN x,long m,long n,long k,long t,ulong
   for (u=t+1; u<=m; u++)
     if (u==j || !c[u])
     {
-      A=coeff(x,u,k);
+      A=(gpmem_t)coeff(x,u,k);
       if (A<av && A>=bot) coeff(x,u,k)+=dec;
     }
   for (u=1; u<=m; u++)
     if (u==j || !c[u])
       for (i=k+1; i<=n; i++)
       {
-        A=coeff(x,u,i);
+        A=(gpmem_t)coeff(x,u,i);
         if (A<av && A>=bot) coeff(x,u,i)+=dec;
       }
 }
@@ -1550,13 +1794,65 @@ gerepile_gauss(GEN x,long m,long n,long k,long t,ulong
 /*          return n - rk(x) linearly independant vectors          */
 /*                                                                 */
 /*******************************************************************/
+static long
+gauss_get_pivot_NZ(GEN x, GEN x0/*unused*/, GEN c, long i0)
+{
+  long i,lx = lg(x);
+  (void)x0;
+  if (c)
+    for (i=i0; i<lx; i++)
+    {
+      if (c[i]) continue; /* already a pivot in line i */
+      if (!gcmp0((GEN)x[i])) break;
+    }
+  else
+    for (i=i0; i<lx; i++)
+      if (!gcmp0((GEN)x[i])) break;
+  return i;
 
+}
+
+/* x ~ 0 compared to reference y */
+int
+approx_0(GEN x, GEN y)
+{
+  long tx = typ(x);
+  if (tx == t_COMPLEX)
+    return approx_0((GEN)x[1], y) && approx_0((GEN)x[2], y);
+  return gcmp0(x) ||
+         (tx == t_REAL && gexpo(y) - gexpo(x) > bit_accuracy(lg(x)));
+}
+
+static long
+gauss_get_pivot_max(GEN x, GEN x0, GEN c, long i0)
+{
+  long i,e, k = i0, ex = - (long)HIGHEXPOBIT, lx = lg(x);
+  GEN p,r;
+  if (c)
+    for (i=i0; i<lx; i++)
+    {
+      if (c[i]) continue;
+      e = gexpo((GEN)x[i]);
+      if (e > ex) { ex=e; k=i; }
+    }
+  else
+    for (i=i0; i<lx; i++)
+    {
+      e = gexpo((GEN)x[i]);
+      if (e > ex) { ex=e; k=i; }
+    }
+  p = (GEN)x[k];
+  r = (GEN)x0[k]; if (isexactzero(r)) r = x0;
+  return approx_0(p, r)? i: k;
+}
+
 /* x has INTEGER coefficients */
 GEN
 keri(GEN x)
 {
-  GEN c,d,y,p,pp;
-  long i,j,k,r,t,n,m,av,av0,tetpil,lim;
+  gpmem_t av, av0, tetpil, lim;
+  GEN c,l,y,p,pp;
+  long i,j,k,r,t,n,m;
 
   if (typ(x)!=t_MAT) err(typeer,"keri");
   n=lg(x)-1; if (!n) return cgetg(1,t_MAT);
@@ -1564,42 +1860,40 @@ keri(GEN x)
   av0=avma; m=lg(x[1])-1; r=0;
   pp=cgetg(n+1,t_COL);
   x=dummycopy(x); p=gun;
-  c=new_chunk(m+1); for (k=1; k<=m; k++) c[k]=0;
-  d=new_chunk(n+1); av=avma; lim=stack_lim(av,1);
+  c=cgetg(m+1, t_VECSMALL); for (k=1; k<=m; k++) c[k]=0;
+  l=cgetg(n+1, t_VECSMALL);
+  av = avma; lim = stack_lim(av,1);
   for (k=1; k<=n; k++)
   {
-    j=1;
-    while (j<=m && (c[j] || !signe(gcoeff(x,j,k))) ) j++;
-    if (j>m)
+    j = 1;
+    while ( j <= m && (c[j] || !signe(gcoeff(x,j,k))) ) j++;
+    if (j > m)
     {
-      r++; d[k]=0;
+      r++; l[k] = 0;
       for(j=1; j<k; j++)
-	if (d[j]) coeff(x,d[j],k) = lclone(gcoeff(x,d[j],k));
+	if (l[j]) coeff(x,l[j],k) = lclone(gcoeff(x,l[j],k));
       pp[k]=lclone(p);
     }
     else
     {
       GEN p0 = p;
-      long av1;
+      c[j]=k; l[k]=j; p = gcoeff(x,j,k);
 
-      c[j]=k; d[k]=j; p = gcoeff(x,j,k);
-
       for (t=1; t<=m; t++)
 	if (t!=j)
 	{
-	  GEN q=gcoeff(x,t,k), p1,p2;
+	  GEN q=gcoeff(x,t,k), p1;
 	  for (i=k+1; i<=n; i++)
 	  {
-	    av1=avma; (void)new_chunk(lgefint(p0));
-	    p1=mulii(q,gcoeff(x,j,i));
-	    p2=mulii(p,gcoeff(x,t,i));
-	    p1=subii(p2,p1); avma=av1;
-	    coeff(x,t,i) = ldivii(p1,p0);
+	    gpmem_t av1 = avma;
+	    p1 = subii(mulii(p,gcoeff(x,t,i)), mulii(q,gcoeff(x,j,i)));
+            p1 = gerepileuptoint(av1, diviiexact(p1,p0));
+	    coeff(x,t,i) = (long)p1;
 	  }
 	  if (low_stack(lim, stack_lim(av,1)))
           {
             p1 = gclone(p);
-            gerepile_gauss_ker(x,m,n,k,t,av);
+            gerepile_gauss_ker(x,k,t,av);
             p = gcopy(p1); gunclone(p1);
           }
 	}
@@ -1612,11 +1906,11 @@ keri(GEN x)
   for (j=k=1; j<=r; j++,k++)
   {
     p = cgetg(n+1, t_COL);
-    y[j]=(long)p; while (d[k]) k++;
+    y[j]=(long)p; while (l[k]) k++;
     for (i=1; i<k; i++)
-      if (d[i])
+      if (l[i])
       {
-	c=gcoeff(x,d[i],k);
+	c=gcoeff(x,l[i],k);
 	p[i] = (long) forcecopy(c); gunclone(c);
       }
       else
@@ -1628,51 +1922,54 @@ keri(GEN x)
 }
 
 GEN
-deplin(GEN x)
+deplin(GEN x0)
 {
-  long i,j,k,t,nc,nl, av=avma;
-  GEN y,q,c,l,d;
+  gpmem_t av = avma;
+  long i,j,k,t,nc,nl;
+  GEN x,y,piv,q,c,l,*d,*ck,*cj;
 
-  if (typ(x) != t_MAT) err(typeer,"deplin");
-  nc=lg(x)-1; if (!nc) err(talker,"empty matrix in deplin");
-  nl=lg(x[1])-1;
-  c=new_chunk(nl+1);
-  l=new_chunk(nc+1);
-  d=cgetg(nl+1,t_VEC);
-  for (i=1; i<=nl; i++) { d[i]=un; c[i]=0; }
-  k=1; t=1;
-  while (t<=nl && k<=nc)
+  t = typ(x0);
+  if (t == t_MAT) x = dummycopy(x0);
+  else
   {
+    if (t != t_VEC) err(typeer,"deplin");
+    x = gtomat(x0);
+  }
+  nc = lg(x)-1; if (!nc) err(talker,"empty matrix in deplin");
+  nl = lg(x[1])-1;
+  d = (GEN*)cgetg(nl+1,t_VEC); /* pivot list */
+  c = cgetg(nl+1, t_VECSMALL);
+  l = cgetg(nc+1, t_VECSMALL); /* not initialized */
+  for (i=1; i<=nl; i++) { d[i] = gun; c[i] = 0; }
+  ck = NULL; /* gcc -Wall */
+  for (k=1; k<=nc; k++)
+  {
+    ck = (GEN*)x[k];
     for (j=1; j<k; j++)
-     for (i=1; i<=nl; i++)
-      if (i!=l[j])
-       coeff(x,i,k)=lsub(gmul((GEN) d[j],gcoeff(x,i,k)),
-                         gmul(gcoeff(x,i,j),gcoeff(x,l[j],k)));
-    t=1;
-    while ( t<=nl && (c[t] || gcmp0(gcoeff(x,t,k))) ) t++;
-    if (t<=nl)
     {
-      d[k]=coeff(x,t,k);
-      c[t]=k; l[k++]=t;
+      cj = (GEN*)x[j]; piv = d[j]; q = gneg(ck[l[j]]);
+      for (i=1; i<=nl; i++)
+        if (i != l[j]) ck[i] = gadd(gmul(piv, ck[i]), gmul(q, cj[i]));
     }
+
+    i = gauss_get_pivot_NZ((GEN)ck, NULL, c, 1);
+    if (i > nl) break;
+
+    d[k] = ck[i];
+    c[i] = k; l[k] = i; /* pivot d[k] in x[i,k] */
   }
-  if (k>nc)
+  if (k > nc) { avma = av; return zerocol(nc); }
+  if (k == 1) { avma = av; return gscalcol_i(gun,nc); }
+  y = cgetg(nc+1,t_COL);
+  y[1] = (long)ck[ l[1] ];
+  for (q=d[1],j=2; j<k; j++)
   {
-    avma=av; y=cgetg(nc+1,t_COL);
-    for (j=1; j<=nc; j++) y[j]=zero;
-    return y;
+    y[j] = lmul(ck[ l[j] ], q);
+    q = gmul(q, d[j]);
   }
-  y=cgetg(nc+1,t_COL);
-  y[1]=(k>1)? coeff(x,l[1],k): un;
-  for (q=gun,j=2; j<k; j++)
-  {
-    q=gmul(q,(GEN) d[j-1]);
-    y[j]=lmul(gcoeff(x,l[j],k),q);
-  }
-  if (k>1) y[k]=lneg(gmul(q,(GEN) d[k-1]));
-  for (j=k+1; j<=nc; j++) y[j]=zero;
-  d=content(y); t=avma;
-  return gerepile(av,t,gdiv(y,d));
+  y[j] = lneg(q);
+  for (j++; j<=nc; j++) y[j] = zero;
+  return gerepileupto(av, gdiv(y,content(y)));
 }
 
 /*******************************************************************/
@@ -1681,78 +1978,6 @@ deplin(GEN x)
 /*           (kernel, image, complementary image, rank)            */
 /*                                                                 */
 /*******************************************************************/
-static long gauss_ex;
-static int (*gauss_is_zero)(GEN);
-
-static int
-real0(GEN x)
-{
-  return gcmp0(x) || (gexpo(x) < gauss_ex);
-}
-
-static void
-gauss_get_prec(GEN x, long prec)
-{
-  long pr = matprec(x);
-
-  if (!pr) { gauss_is_zero = &gcmp0; return; }
-  if (pr > prec) prec = pr;
-
-  gauss_ex = - (long)(0.85 * bit_accuracy(prec));
-  gauss_is_zero = &real0;
-}
-
-static long
-gauss_get_pivot_NZ(GEN x, GEN x0/* unused */, GEN c, long i0)
-{
-  long i,lx = lg(x);
-  if (c)
-    for (i=i0; i<lx; i++)
-    {
-      if (c[i]) continue;
-      if (!gcmp0((GEN)x[i])) break;
-    }
-  else
-    for (i=i0; i<lx; i++)
-      if (!gcmp0((GEN)x[i])) break;
-  return i;
-
-}
-
-/* ~ 0 compared to reference y */
-static int
-approx_0(GEN x, GEN y)
-{
-  long tx = typ(x);
-  if (tx == t_COMPLEX)
-    return approx_0((GEN)x[1], y) && approx_0((GEN)x[2], y);
-  return gcmp0(x) ||
-         (tx == t_REAL && gexpo(y) - gexpo(x) > bit_accuracy(lg(x)));
-}
-
-static long
-gauss_get_pivot_max(GEN x, GEN x0, GEN c, long i0)
-{
-  long i,e, k = i0, ex = -HIGHEXPOBIT, lx = lg(x);
-  GEN p,r;
-  if (c)
-    for (i=i0; i<lx; i++)
-    {
-      if (c[i]) continue;
-      e = gexpo((GEN)x[i]);
-      if (e > ex) { ex=e; k=i; }
-    }
-  else
-    for (i=i0; i<lx; i++)
-    {
-      e = gexpo((GEN)x[i]);
-      if (e > ex) { ex=e; k=i; }
-    }
-  p = (GEN)x[k];
-  r = (GEN)x0[k]; if (isexactzero(r)) r = x0;
-  return approx_0(p, r)? i: k;
-}
-
 /* return the transform of x under a standard Gauss pivot. r = dim ker(x).
  * d[k] contains the index of the first non-zero pivot in column k
  * If a != NULL, use x - a Id instead (for eigen)
@@ -1761,7 +1986,8 @@ static GEN
 gauss_pivot_ker(GEN x0, GEN a, GEN *dd, long *rr)
 {
   GEN x,c,d,p,mun;
-  long i,j,k,r,t,n,m,av,lim;
+  gpmem_t av, lim;
+  long i,j,k,r,t,n,m;
   long (*get_pivot)(GEN,GEN,GEN,long);
 
   if (typ(x0)!=t_MAT) err(typeer,"gauss_pivot");
@@ -1801,7 +2027,7 @@ gauss_pivot_ker(GEN x0, GEN a, GEN *dd, long *rr)
 	  for (i=k+1; i<=n; i++)
 	    coeff(x,t,i) = ladd(gcoeff(x,t,i),gmul(p,gcoeff(x,j,i)));
           if (low_stack(lim, stack_lim(av,1)))
-            gerepile_gauss_ker(x,m,n,k,t,av);
+            gerepile_gauss_ker(x,k,t,av);
 	}
     }
   }
@@ -1815,7 +2041,8 @@ static void
 gauss_pivot(GEN x0, GEN *dd, long *rr)
 {
   GEN x,c,d,d0,mun,p;
-  long i,j,k,r,t,n,m,av,lim;
+  long i, j, k, r, t, n, m;
+  gpmem_t av, lim;
   long (*get_pivot)(GEN,GEN,GEN,long);
 
   if (typ(x0)!=t_MAT) err(typeer,"gauss_pivot");
@@ -1856,7 +2083,7 @@ gauss_pivot(GEN x0, GEN *dd, long *rr)
           for (i=k+1; i<=n; i++)
             coeff(x,t,i) = ladd(gcoeff(x,t,i), gmul(p,gcoeff(x,j,i)));
           if (low_stack(lim, stack_lim(av,1)))
-            gerepile_gauss(x,m,n,k,t,av,j,c);
+            gerepile_gauss(x,k,t,av,j,c);
         }
       for (i=k; i<=n; i++) coeff(x,j,i) = zero; /* dummy */
     }
@@ -1868,8 +2095,9 @@ gauss_pivot(GEN x0, GEN *dd, long *rr)
 static GEN
 ker0(GEN x, GEN a)
 {
+  gpmem_t av = avma, tetpil;
   GEN d,y;
-  long i,j,k,r,n, av = avma, tetpil;
+  long i,j,k,r,n;
 
   x = gauss_pivot_ker(x,a,&d,&r);
   if (!r) { avma=av; return cgetg(1,t_MAT); }
@@ -1907,8 +2135,9 @@ matker0(GEN x,long flag)
 GEN
 image(GEN x)
 {
+  gpmem_t av = avma;
   GEN d,y;
-  long j,k,r, av = avma;
+  long j,k,r;
 
   gauss_pivot(x,&d,&r);
 
@@ -1926,8 +2155,9 @@ image(GEN x)
 GEN
 imagecompl(GEN x)
 {
+  gpmem_t av = avma;
   GEN d,y;
-  long j,i,r,av = avma;
+  long j,i,r;
 
   gauss_pivot(x,&d,&r);
   avma=av; y=cgetg(r+1,t_VEC);
@@ -1940,10 +2170,11 @@ imagecompl(GEN x)
 GEN
 imagecomplspec(GEN x, long *nlze)
 {
+  gpmem_t av = avma;
   GEN d,y;
-  long i,j,k,l,r,av = avma;
+  long i,j,k,l,r;
 
-  x = gtrans(x); l = lg(x);
+  x = gtrans_i(x); l = lg(x);
   gauss_pivot(x,&d,&r);
   avma=av; y = cgetg(l,t_VECSMALL);
   for (i=j=1, k=r+1; i<l; i++)
@@ -1955,7 +2186,8 @@ imagecomplspec(GEN x, long *nlze)
 static GEN
 sinverseimage(GEN mat, GEN y)
 {
-  long av=avma,tetpil,i, nbcol = lg(mat);
+  gpmem_t av=avma,tetpil;
+  long i, nbcol = lg(mat);
   GEN col,p1 = cgetg(nbcol+1,t_MAT);
 
   if (nbcol==1) return NULL;
@@ -1977,7 +2209,8 @@ sinverseimage(GEN mat, GEN y)
 GEN
 inverseimage(GEN m,GEN v)
 {
-  long av=avma,j,lv,tv=typ(v);
+  gpmem_t av=avma;
+  long j,lv,tv=typ(v);
   GEN y,p1;
 
   if (typ(m)!=t_MAT) err(typeer,"inverseimage");
@@ -1999,47 +2232,92 @@ inverseimage(GEN m,GEN v)
   return y;
 }
 
-/* x is an n x k matrix, rank(x) = k <= n. Return an invertible n x n matrix
- * whose first k columns are given by x. If rank(x)<k, the result may be wrong
- */
+/* i-th vector in the standard basis */
 GEN
-suppl_intern(GEN x, GEN myid)
+_ei(long n, long i)
 {
-  long av = avma, lx = lg(x), n,i,j;
-  GEN y,p1;
-  stackzone *zone;
+  GEN e = zerocol(n); e[i] = un; return e;
+}
 
-  if (typ(x) != t_MAT) err(typeer,"suppl");
-  if (lx==1) err(talker,"empty matrix in suppl");
-  n=lg(x[1]); if (lx>n) err(suppler2);
-  if (lx == n) return gcopy(x);
+/* NB: d freed */
+static GEN
+get_suppl(GEN x, GEN d, long r)
+{
+  gpmem_t av;
+  GEN y,c;
+  long j,k,rx,n;
 
-  zone  = switch_stack(NULL, n*n);
-  switch_stack(zone,1);
-  y = myid? dummycopy(myid): idmat(n-1);
-  switch_stack(zone,0);
-  gauss_get_prec(x,0);
-  for (i=1; i<lx; i++)
-  {
-    p1 = gauss(y,(GEN)x[i]); j=i;
-    while (j<n && gauss_is_zero((GEN)p1[j])) j++;
-    if (j>=n) err(suppler2);
-    p1=(GEN)y[i]; y[i]=x[i]; if (i!=j) y[j]=(long)p1;
-  }
-  avma = av; y = gcopy(y);
-  free(zone); return y;
+  rx = lg(x)-1;
+  if (!rx) err(talker,"empty matrix in suppl");
+  n = lg(x[1])-1;
+  if (rx == n && r == 0) { free(d); return gcopy(x); }
+  y = cgetg(n+1, t_MAT);
+  av = avma;
+  c = vecsmall_const(n,0);
+  k = 1;
+  /* c = lines containing pivots (could get it from gauss_pivot, but cheap)
+   * In theory r = 0 and d[j] > 0 for all j, but why take chances? */
+  for (j=1; j<=rx; j++)
+    if (d[j])
+    {
+      c[ d[j] ] = 1;
+      y[k++] = x[j];
+    }
+  for (j=1; j<=n; j++)
+    if (!c[j]) y[k++] = j;
+  avma = av;
+
+  rx -= r;
+  for (j=1; j<=rx; j++)
+    y[j] = lcopy((GEN)y[j]);
+  for (   ; j<=n; j++)
+    y[j] = (long)_ei(n, y[j]);
+  free(d); return y;
 }
 
+/* x is an n x k matrix, rank(x) = k <= n. Return an invertible n x n matrix
+ * whose first k columns are given by x. If rank(x) < k, undefined result. */
 GEN
 suppl(GEN x)
 {
-  return suppl_intern(x,NULL);
+  gpmem_t av = avma;
+  GEN d;
+  long r;
+
+  gauss_pivot(x,&d,&r);
+  avma = av; return get_suppl(x,d,r);
 }
 
+static void FpM_gauss_pivot(GEN x, GEN p, GEN *dd, long *rr);
+static void FqM_gauss_pivot(GEN x, GEN T, GEN p, GEN *dd, long *rr);
+
 GEN
+FpM_suppl(GEN x, GEN p)
+{
+  gpmem_t av = avma;
+  GEN d;
+  long r;
+
+  FpM_gauss_pivot(x,p, &d,&r);
+  avma = av; return get_suppl(x,d,r);
+}
+GEN
+FqM_suppl(GEN x, GEN T, GEN p)
+{
+  gpmem_t av = avma;
+  GEN d;
+  long r;
+
+  if (!T) return FpM_suppl(x,p);
+  FqM_gauss_pivot(x,T,p, &d,&r);
+  avma = av; return get_suppl(x,d,r);
+}
+
+GEN
 image2(GEN x)
 {
-  long av=avma,tetpil,k,n,i;
+  gpmem_t av=avma,tetpil;
+  long k,n,i;
   GEN p1,p2;
 
   if (typ(x)!=t_MAT) err(typeer,"image2");
@@ -2068,7 +2346,8 @@ matimage0(GEN x,long flag)
 long
 rank(GEN x)
 {
-  long av = avma, r;
+  gpmem_t av = avma;
+  long r;
   GEN d;
 
   gauss_pivot(x,&d,&r);
@@ -2078,13 +2357,12 @@ rank(GEN x)
   return lg(x)-1 - r;
 }
 
-static void FpM_gauss_pivot(GEN x, GEN p, GEN *dd, long *rr);
-
 /* if p != NULL, assume x integral and compute rank over Fp */
 static GEN
 indexrank0(GEN x, GEN p, int small)
 {
-  long av = avma, i,j,n,r;
+  gpmem_t av = avma;
+  long i,j,n,r;
   GEN res,d,p1,p2;
 
   /* yield r = dim ker(x) */
@@ -2108,13 +2386,13 @@ indexrank0(GEN x, GEN p, int small)
   return res;
 }
 
-GEN 
+GEN
 indexrank(GEN x) { return indexrank0(x,NULL,0); }
 
-GEN 
+GEN
 sindexrank(GEN x) { return indexrank0(x,NULL,1); }
 
-GEN 
+GEN
 FpM_sindexrank(GEN x, GEN p) { return indexrank0(x,p,1); }
 
 /*******************************************************************/
@@ -2122,7 +2400,7 @@ FpM_sindexrank(GEN x, GEN p) { return indexrank0(x,p,1
 /*                    LINEAR ALGEBRA MODULO P                      */
 /*                                                                 */
 /*******************************************************************/
-  
+
 /*If p is NULL no reduction is performed.*/
 GEN
 FpM_mul(GEN x, GEN y, GEN p)
@@ -2130,7 +2408,7 @@ FpM_mul(GEN x, GEN y, GEN p)
   long i,j,l,lx=lg(x), ly=lg(y);
   GEN z;
   if (ly==1) return cgetg(ly,t_MAT);
-  if (lx != lg(y[1])) err(operi,"* [mod p]",t_MAT,t_MAT);
+  if (lx != lg(y[1])) err(operi,"* [mod p]",x,y);
   z=cgetg(ly,t_MAT);
   if (lx==1)
   {
@@ -2142,8 +2420,8 @@ FpM_mul(GEN x, GEN y, GEN p)
   {
     z[j] = lgetg(l,t_COL);
     for (i=1; i<l; i++)
-    { 
-      ulong av;
+    {
+      gpmem_t av;
       GEN p1,p2;
       int k;
       p1=gzero; av=avma;
@@ -2162,45 +2440,37 @@ FpM_mul(GEN x, GEN y, GEN p)
 GEN
 FpM_FpV_mul(GEN x, GEN y, GEN p)
 {
-  long i,l,lx=lg(x), ly=lg(y);
+  long i,k,l,lx=lg(x), ly=lg(y);
   GEN z;
-  if (lx != ly) err(operi,"* [mod p]",t_MAT,t_COL);
-  z=cgetg(ly,t_COL);
+  if (lx != ly) err(operi,"* [mod p]",x,y);
   if (lx==1) return cgetg(1,t_COL);
-  l=lg(x[1]);
+  l = lg(x[1]);
+  z = cgetg(l,t_COL);
   for (i=1; i<l; i++)
-  { 
-    ulong av;
-    GEN p1,p2;
-    int k;
-    p1=gzero; av=avma;
+  {
+    gpmem_t av = avma;
+    GEN p1 = gzero;
     for (k=1; k<lx; k++)
-    {
-      p2=mulii(gcoeff(x,i,k),(GEN)y[k]);
-      p1=addii(p1,p2);
-    }
-    z[i]=lpileupto(av,p?modii(p1,p):p1);
+      p1 = addii(p1, mulii(gcoeff(x,i,k),(GEN)y[k]));
+    z[i] = lpileupto(av,p?modii(p1,p):p1);
   }
   return z;
 }
 
+/* in place, destroy x */
 static GEN
-u_FpM_ker(GEN x, GEN pp, long nontriv)
+u_FpM_ker_sp(GEN x, ulong p, long deplin)
 {
   GEN y,c,d;
-  long a,i,j,k,r,t,n,m,av0,tetpil, p = pp[2], piv;
+  long i,j,k,r,t,n;
+  ulong a, piv, m;
 
   n = lg(x)-1;
-  m=lg(x[1])-1; r=0; av0 = avma;
-  x = dummycopy(x);
-  for (i=1; i<=n; i++)
-  {
-    GEN p1 = (GEN)x[i];
-    for (j=1; j<=m; j++) p1[j] = itos((GEN)p1[j]);
-  }
-  c=new_chunk(m+1); for (k=1; k<=m; k++) c[k]=0;
-  d=new_chunk(n+1);
-  a=0; /* for gcc -Wall */
+  m=lg(x[1])-1; r=0;
+
+  c = new_chunk(m+1); for (k=1; k<=m; k++) c[k] = 0;
+  d = new_chunk(n+1);
+  a = 0; /* for gcc -Wall */
   for (k=1; k<=n; k++)
   {
     for (j=1; j<=m; j++)
@@ -2209,62 +2479,82 @@ u_FpM_ker(GEN x, GEN pp, long nontriv)
         a = coeff(x,j,k) % p;
         if (a) break;
       }
-    if (j>m)
+    if (j > m)
     {
-      if (nontriv) { avma=av0; return NULL; }
-      r++; d[k]=0;
+      if (deplin) {
+        c = cgetg(n+1, t_VECSMALL);
+        for (i=1; i<k; i++) c[i] = coeff(x,d[i],k) % p;
+        c[k] = 1; for (i=k+1; i<=n; i++) c[i] = 0;
+        return c;
+      }
+      r++; d[k] = 0;
     }
     else
     {
-      c[j]=k; d[k]=j;
-      piv = - u_Fp_inv(a, p);
-      coeff(x,j,k) = -1;
+      c[j] = k; d[k] = j;
+      piv = p - u_invmod(a, p); /* -1/a */
+      coeff(x,j,k) = p-1;
       for (i=k+1; i<=n; i++)
 	coeff(x,j,i) = (piv * coeff(x,j,i)) % p;
       for (t=1; t<=m; t++)
-	if (t!=j)
-	{
-	  piv = coeff(x,t,k) % p;
-          if (piv)
-          {
-            coeff(x,t,k) = 0;
+      {
+	if (t==j) continue;
+
+        piv = coeff(x,t,k) % p;
+        if (!piv) continue;
+
+        coeff(x,t,k) = 0;
+        if (piv == 1)
+          for (i=k+1; i<=n; i++) _u_Fp_add((GEN)x[i],t,j,p);
+        else
+        {
+          if (u_OK_ULONG(p))
+            for (i=k+1; i<=n; i++) _u_Fp_addmul_OK((GEN)x[i],t,j,piv,p);
+          else
             for (i=k+1; i<=n; i++) _u_Fp_addmul((GEN)x[i],t,j,piv,p);
-          }
-	}
+        }
+      }
     }
   }
+  if (deplin) return NULL;
 
-  tetpil=avma; y=cgetg(r+1,t_MAT);
+  y = cgetg(r+1, t_MAT);
   for (j=k=1; j<=r; j++,k++)
   {
-    GEN c = cgetg(n+1,t_COL);
+    GEN c = cgetg(n+1, t_VECSMALL);
 
-    y[j]=(long)c; while (d[k]) k++;
+    y[j] = (long)c; while (d[k]) k++;
     for (i=1; i<k; i++)
       if (d[i])
-      {
-        long a = coeff(x,d[i],k) % p;
-        if (a < 0) a += p;
-	c[i] = lstoi(a);
-      }
+        c[i] = coeff(x,d[i],k) % p;
       else
-	c[i] = zero;
-    c[k]=un; for (i=k+1; i<=n; i++) c[i]=zero;
+	c[i] = 0;
+    c[k] = 1; for (i=k+1; i<=n; i++) c[i] = 0;
   }
-  return gerepile(av0,tetpil,y);
+  return y;
 }
 
+/* assume x has integer entries */
 static GEN
-FpM_ker_i(GEN x, GEN p, long nontriv)
+FpM_ker_i(GEN x, GEN p, long deplin)
 {
+  gpmem_t av0 = avma, av,lim,tetpil;
   GEN y,c,d,piv,mun;
-  long i,j,k,r,t,n,m,av0,av,lim,tetpil;
+  long i,j,k,r,t,n,m;
 
   if (typ(x)!=t_MAT) err(typeer,"FpM_ker");
   n=lg(x)-1; if (!n) return cgetg(1,t_MAT);
-  if (OK_ULONG(p)) return u_FpM_ker(x, p, nontriv);
+  if (OK_ULONG(p))
+  {
+    ulong pp = p[2];
+    y = u_Fp_FpM(x, pp);
+    y = u_FpM_ker_sp(y, pp, deplin);
+    if (!y) return y;
+    y = deplin? small_to_col(y): small_to_mat(y);
+    return gerepileupto(av0, y);
+  }
 
-  m=lg(x[1])-1; r=0; av0 = avma;
+  m=lg(x[1])-1; r=0;
   x=dummycopy(x); mun=negi(gun);
   c=new_chunk(m+1); for (k=1; k<=m; k++) c[k]=0;
   d=new_chunk(n+1);
@@ -2279,7 +2569,12 @@ FpM_ker_i(GEN x, GEN p, long nontriv)
       }
     if (j>m)
     {
-      if (nontriv) { avma = av0; return NULL; }
+      if (deplin) {
+        c = cgetg(n+1, t_COL);
+        for (i=1; i<k; i++) c[i] = lmodii(gcoeff(x,d[i],k), p);
+        c[k] = un; for (i=k+1; i<=n; i++) c[i] = zero;
+        return gerepileupto(av0, c);
+      }
       r++; d[k]=0;
       for(j=1; j<k; j++)
         if (d[j]) coeff(x,d[j],k) = lclone(gcoeff(x,d[j],k));
@@ -2291,20 +2586,21 @@ FpM_ker_i(GEN x, GEN p, long nontriv)
       for (i=k+1; i<=n; i++)
 	coeff(x,j,i) = lmodii(mulii(piv,gcoeff(x,j,i)), p);
       for (t=1; t<=m; t++)
-	if (t!=j)
-	{
-	  piv = modii(gcoeff(x,t,k), p);
-          if (signe(piv))
-          {
-            coeff(x,t,k)=zero;
-            for (i=k+1; i<=n; i++)
-              coeff(x,t,i) = laddii(gcoeff(x,t,i),mulii(piv,gcoeff(x,j,i)));
-            if (low_stack(lim, stack_lim(av,1)))
-              gerepile_gauss_FpM_ker(x,p,m,n,k,t,av);
-          }
-	}
+      {
+	if (t==j) continue;
+
+        piv = modii(gcoeff(x,t,k), p);
+        if (!signe(piv)) continue;
+
+        coeff(x,t,k)=zero;
+        for (i=k+1; i<=n; i++)
+          coeff(x,t,i) = laddii(gcoeff(x,t,i),mulii(piv,gcoeff(x,j,i)));
+        if (low_stack(lim, stack_lim(av,1)))
+          gerepile_gauss_FpM_ker(x,p,k,t,av);
+      }
     }
   }
+  if (deplin) { avma = av0; return NULL; }
 
   tetpil=avma; y=cgetg(r+1,t_MAT);
   for (j=k=1; j<=r; j++,k++)
@@ -2325,13 +2621,37 @@ FpM_ker_i(GEN x, GEN p, long nontriv)
   return gerepile(av0,tetpil,y);
 }
 
+GEN
+FpM_intersect(GEN x, GEN y, GEN p)
+{
+  gpmem_t av = avma;
+  long j, lx = lg(x);
+  GEN z;
+
+  if (lx==1 || lg(y)==1) return cgetg(1,t_MAT);
+  z = FpM_ker(concatsp(x,y), p);
+  for (j=lg(z)-1; j; j--) setlg(z[j],lx);
+  return gerepileupto(av, FpM_mul(x,z,p));
+}
+
+GEN
+FpM_ker(GEN x, GEN p) { return FpM_ker_i(x, p, 0); }
+GEN
+FpM_deplin(GEN x, GEN p) { return FpM_ker_i(x, p, 1); }
+/* not memory clean */
+GEN
+u_FpM_ker(GEN x, ulong p) { return u_FpM_ker_sp(dummycopy(x), p, 0); }
+GEN
+u_FpM_deplin(GEN x, ulong p) { return u_FpM_ker_sp(dummycopy(x), p, 1); }
+
 static void
 FpM_gauss_pivot(GEN x, GEN p, GEN *dd, long *rr)
 {
+  gpmem_t av,lim;
   GEN c,d,piv;
-  long i,j,k,r,t,n,m,av,lim;
+  long i,j,k,r,t,n,m;
 
-  if (!p) return gauss_pivot(x,dd,rr);
+  if (!p) { gauss_pivot(x,dd,rr); return; }
   if (typ(x)!=t_MAT) err(typeer,"FpM_gauss_pivot");
   n=lg(x)-1; if (!n) { *dd=NULL; *rr=0; return; }
 
@@ -2363,7 +2683,7 @@ FpM_gauss_pivot(GEN x, GEN p, GEN *dd, long *rr)
             for (i=k+1; i<=n; i++)
               coeff(x,t,i) = laddii(gcoeff(x,t,i), mulii(piv,gcoeff(x,j,i)));
             if (low_stack(lim, stack_lim(av,1)))
-              gerepile_gauss(x,m,n,k,t,av,j,c);
+              gerepile_gauss(x,k,t,av,j,c);
           }
         }
       for (i=k; i<=n; i++) coeff(x,j,i) = zero; /* dummy */
@@ -2371,24 +2691,59 @@ FpM_gauss_pivot(GEN x, GEN p, GEN *dd, long *rr)
   }
   *dd=d; *rr=r;
 }
-
-GEN
-FpM_ker(GEN x, GEN p)
+static void
+FqM_gauss_pivot(GEN x, GEN T, GEN p, GEN *dd, long *rr)
 {
-  return FpM_ker_i(x, p, 0);
-}
+  gpmem_t av,lim;
+  GEN c,d,piv;
+  long i,j,k,r,t,n,m;
 
-int
-ker_trivial_mod_p(GEN x, GEN p)
-{
-  return FpM_ker_i(x, p, 1)==NULL;
+  if (typ(x)!=t_MAT) err(typeer,"FqM_gauss_pivot");
+  n=lg(x)-1; if (!n) { *dd=NULL; *rr=0; return; }
+
+  m=lg(x[1])-1; r=0;
+  x=dummycopy(x);
+  c=new_chunk(m+1); for (k=1; k<=m; k++) c[k]=0;
+  d=(GEN)gpmalloc((n+1)*sizeof(long)); av=avma; lim=stack_lim(av,1);
+  for (k=1; k<=n; k++)
+  {
+    for (j=1; j<=m; j++)
+      if (!c[j])
+      {
+        coeff(x,j,k) = (long)FpX_res(gcoeff(x,j,k), T,p);
+        if (signe(coeff(x,j,k))) break;
+      }
+    if (j>m) { r++; d[k]=0; }
+    else
+    {
+      c[j]=k; d[k]=j; piv = gneg(FpXQ_inv(gcoeff(x,j,k), T,p));
+      for (i=k+1; i<=n; i++)
+	coeff(x,j,i) = (long)Fq_mul(piv,gcoeff(x,j,i), T, p);
+      for (t=1; t<=m; t++)
+        if (!c[t]) /* no pivot on that line yet */
+        {
+          piv=gcoeff(x,t,k);
+          if (signe(piv))
+          {
+            coeff(x,t,k)=zero;
+            for (i=k+1; i<=n; i++)
+              coeff(x,t,i) = ladd(gcoeff(x,t,i), gmul(piv,gcoeff(x,j,i)));
+            if (low_stack(lim, stack_lim(av,1)))
+              gerepile_gauss(x,k,t,av,j,c);
+          }
+        }
+      for (i=k; i<=n; i++) coeff(x,j,i) = zero; /* dummy */
+    }
+  }
+  *dd=d; *rr=r;
 }
 
 GEN
 FpM_image(GEN x, GEN p)
 {
+  gpmem_t av = avma;
   GEN d,y;
-  long j,k,r, av = avma;
+  long j,k,r;
 
   FpM_gauss_pivot(x,p,&d,&r);
 
@@ -2406,7 +2761,8 @@ FpM_image(GEN x, GEN p)
 static GEN
 sFpM_invimage(GEN mat, GEN y, GEN p)
 {
-  long av=avma,i, nbcol = lg(mat);
+  gpmem_t av=avma;
+  long i, nbcol = lg(mat);
   GEN col,p1 = cgetg(nbcol+1,t_MAT),res;
 
   if (nbcol==1) return NULL;
@@ -2434,7 +2790,8 @@ sFpM_invimage(GEN mat, GEN y, GEN p)
 GEN
 FpM_invimage(GEN m, GEN v, GEN p)
 {
-  long av=avma,j,lv,tv=typ(v);
+  gpmem_t av=avma;
+  long j,lv,tv=typ(v);
   GEN y,p1;
 
   if (typ(m)!=t_MAT) err(typeer,"inverseimage");
@@ -2455,7 +2812,7 @@ FpM_invimage(GEN m, GEN v, GEN p)
   }
   return y;
 }
-/************************************************************** 
+/**************************************************************
  Rather stupid implementation of linear algebra in finite fields
  **************************************************************/
 static GEN
@@ -2467,7 +2824,7 @@ Fq_add(GEN x, GEN y, GEN T/*unused*/, GEN p)
     case 1: return FpX_Fp_add(x,y,p);
     case 2: return FpX_Fp_add(y,x,p);
     case 3: return FpX_add(x,y,p);
-  } 
+  }
   return NULL;
 }
 #if 0
@@ -2475,7 +2832,7 @@ Fq_add(GEN x, GEN y, GEN T/*unused*/, GEN p)
  * For consistency we write the code there.
  * To avoid having to remove static status, we rewrite it in polarit1.c
  */
-static GEN 
+static GEN
 Fq_neg(GEN x, GEN T, GEN p)
 {
   switch(typ(x)==t_POL)
@@ -2487,7 +2844,7 @@ Fq_neg(GEN x, GEN T, GEN p)
 }
 #endif
 
-static GEN
+GEN
 Fq_mul(GEN x, GEN y, GEN T, GEN p)
 {
   switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
@@ -2514,7 +2871,7 @@ Fq_neg_inv(GEN x, GEN T, GEN p)
 static GEN
 Fq_res(GEN x, GEN T, GEN p)
 {
-  ulong ltop=avma;
+  gpmem_t ltop=avma;
   switch(typ(x)==t_POL)
   {
     case 0: return modii(x,p);
@@ -2524,11 +2881,11 @@ Fq_res(GEN x, GEN T, GEN p)
 }
 
 static void
-Fq_gerepile_gauss_ker(GEN x, GEN T, GEN p, long m, long n, long k, long t,
-                      ulong av)
+Fq_gerepile_gauss_ker(GEN x, GEN T, GEN p, long k, long t, gpmem_t av)
 {
-  ulong tetpil = avma,A;
-  long dec,u,i;
+  gpmem_t tetpil = avma,A;
+  long u,i, n = lg(x)-1, m = n? lg(x[1])-1: 0;
+  size_t dec;
 
   if (DEBUGMEM > 1) err(warnmem,"gauss_pivot_ker. k=%ld, n=%ld",k,n);
   for (u=t+1; u<=m; u++)
@@ -2540,23 +2897,26 @@ Fq_gerepile_gauss_ker(GEN x, GEN T, GEN p, long m, lon
   (void)gerepile(av,tetpil,NULL); dec = av-tetpil;
   for (u=t+1; u<=m; u++)
   {
-    A=coeff(x,u,k);
+    A=(gpmem_t)coeff(x,u,k);
     if (A<av && A>=bot) coeff(x,u,k)+=dec;
   }
   for (i=k+1; i<=n; i++)
     for (u=1; u<=m; u++)
     {
-      A=coeff(x,u,i);
+      A=(gpmem_t)coeff(x,u,i);
       if (A<av && A>=bot) coeff(x,u,i)+=dec;
     }
 }
 
 static GEN
-FqM_ker_i(GEN x, GEN T, GEN p, long nontriv)
+FqM_ker_i(GEN x, GEN T, GEN p, long deplin)
 {
+  gpmem_t av0,av,lim,tetpil;
   GEN y,c,d,piv,mun;
-  long i,j,k,r,t,n,m,av0,av,lim,tetpil;
+  long i,j,k,r,t,n,m;
 
+  if (!T) return FpM_ker_i(x,p,deplin);
+
   if (typ(x)!=t_MAT) err(typeer,"FpM_ker");
   n=lg(x)-1; if (!n) return cgetg(1,t_MAT);
 
@@ -2575,7 +2935,12 @@ FqM_ker_i(GEN x, GEN T, GEN p, long nontriv)
       }
     if (j>m)
     {
-      if (nontriv) { avma = av0; return NULL; }
+      if (deplin) {
+        c = cgetg(n+1, t_COL);
+        for (i=1; i<k; i++) c[i] = (long)Fq_res(gcoeff(x,d[i],k), T, p);
+        c[k] = un; for (i=k+1; i<=n; i++) c[i] = zero;
+        return gerepileupto(av0, c); 
+      }
       r++; d[k]=0;
       for(j=1; j<k; j++)
         if (d[j]) coeff(x,d[j],k) = lclone(gcoeff(x,d[j],k));
@@ -2587,20 +2952,21 @@ FqM_ker_i(GEN x, GEN T, GEN p, long nontriv)
       for (i=k+1; i<=n; i++)
 	coeff(x,j,i) = (long) Fq_mul(piv,gcoeff(x,j,i), T, p);
       for (t=1; t<=m; t++)
-	if (t!=j)
-	{
-	  piv = Fq_res(gcoeff(x,t,k), T, p);
-          if (signe(piv))
-          {
-            coeff(x,t,k)=zero;
-            for (i=k+1; i<=n; i++)
-              coeff(x,t,i) = (long)Fq_add(gcoeff(x,t,i),Fq_mul(piv,gcoeff(x,j,i), T, p), T, p);
-            if (low_stack(lim, stack_lim(av,1)))
-              Fq_gerepile_gauss_ker(x,T,p,m,n,k,t,av);
-          }
-	}
+      {
+        if (t==j) continue;
+
+        piv = Fq_res(gcoeff(x,t,k), T, p);
+        if (!signe(piv)) continue;
+
+        coeff(x,t,k)=zero;
+        for (i=k+1; i<=n; i++)
+          coeff(x,t,i) = (long)Fq_add(gcoeff(x,t,i),Fq_mul(piv,gcoeff(x,j,i), T, p), T, p);
+        if (low_stack(lim, stack_lim(av,1)))
+          Fq_gerepile_gauss_ker(x,T,p,k,t,av);
+      }
     }
   }
+  if (deplin) { avma = av0; return NULL; }
 
   tetpil=avma; y=cgetg(r+1,t_MAT);
   for (j=k=1; j<=r; j++,k++)
@@ -2638,7 +3004,7 @@ eigen(GEN x, long prec)
 {
   GEN y,rr,p,ssesp,r1,r2,r3;
   long e,i,k,l,ly,ex, n = lg(x);
-  ulong av = avma;
+  gpmem_t av = avma;
 
   if (typ(x)!=t_MAT) err(typeer,"eigen");
   if (n != 1 && n != lg(x[1])) err(mattype1,"eigen");
@@ -2696,9 +3062,10 @@ det0(GEN a,long flag)
 
 /* Exact types: choose the first non-zero pivot. Otherwise: maximal pivot */
 static GEN
-det_simple_gauss(GEN a, long inexact)
+det_simple_gauss(GEN a, int inexact)
 {
-  long i,j,k,av,av1,s, nbco = lg(a)-1;
+  gpmem_t av, av1;
+  long i,j,k,s, nbco = lg(a)-1;
   GEN x,p;
 
   av=avma; s=1; x=gun; a=dummycopy(a);
@@ -2745,63 +3112,14 @@ det_simple_gauss(GEN a, long inexact)
   av1=avma; return gerepile(av,av1,gmul(x,gcoeff(a,nbco,nbco)));
 }
 
-/* a has integer entries, N = P^n */
 GEN
-det_mod_P_n(GEN a, GEN N, GEN P)
-{
-  long va,i,j,k,s, av = avma, nbco = lg(a)-1;
-  GEN x,p;
-
-  s=1; va=0; x=gun; a=dummycopy(a);
-  p = NULL; /* for gcc -Wall */
-  for (i=1; i<nbco; i++)
-  {
-    long fl = 0;
-    for(;;)
-    {
-      for (k=i; k<=nbco; k++)
-      {
-        p=gcoeff(a,i,k);
-        if (signe(p))
-        {
-          fl = 1;
-          if (resii(p,P) != gzero) break;
-        }
-      }
-      if (k <= nbco) break;
-      va++; N = divii(N, P);
-      if (!fl || is_pm1(N)) { avma=av; return gzero; }
-      for (k=i; k<=nbco; k++) coeff(a,i,k) = ldivii(gcoeff(a,i,k), P);
-    }
-    if (k != i) { swap(a[i],a[k]); s = -s; }
-
-    x = gmul(x,p); p = mpinvmod(p,N);
-    for (k=i+1; k<=nbco; k++)
-    {
-      GEN m = resii(gcoeff(a,i,k), N);
-      coeff(a,i,k) = zero;
-      if (signe(m))
-      {
-	m = negi(resii(mulii(m,p), N));
-	for (j=i+1; j<=nbco; j++)
-	  coeff(a,j,k) = lresii(addii(gcoeff(a,j,k),
-                                mulii(m,gcoeff(a,j,i))), N);
-      }
-    }
-  }
-  if (s<0) x = negi(x);
-  x = resii(mulii(x,gcoeff(a,nbco,nbco)), N);
-  return gerepileuptoint(av, mulii(x, gpowgs(P,va)));
-}
-
-GEN
 det2(GEN a)
 {
   long nbco = lg(a)-1;
   if (typ(a)!=t_MAT) err(mattype1,"det2");
   if (!nbco) return gun;
   if (nbco != lg(a[1])-1) err(mattype1,"det2");
-  return det_simple_gauss(a,use_maximal_pivot(a));
+  return det_simple_gauss(a, use_maximal_pivot(a));
 }
 
 static GEN
@@ -2818,7 +3136,7 @@ mydiv(GEN x, GEN y)
 GEN
 det(GEN a)
 {
-  ulong av, lim;
+  gpmem_t av, lim;
   long nbco = lg(a)-1,i,j,k,s;
   GEN p,pprec;
 
@@ -2826,9 +3144,9 @@ det(GEN a)
   if (!nbco) return gun;
   if (nbco != lg(a[1])-1) err(mattype1,"det");
   if (use_maximal_pivot(a)) return det_simple_gauss(a,1);
-  if (DEBUGLEVEL > 7) timer2();
+  if (DEBUGLEVEL > 7) (void)timer2();
 
-  av = avma; lim = stack_lim(av,2); 
+  av = avma; lim = stack_lim(av,2);
   a = dummycopy(a); s = 1;
   for (pprec=gun,i=1; i<nbco; i++,pprec=p)
   {
@@ -2892,7 +3210,8 @@ det(GEN a)
 static GEN
 gaussmoduloall(GEN M, GEN D, GEN Y, GEN *ptu1)
 {
-  long n,m,i,j,lM,av=avma,tetpil;
+  gpmem_t av = avma, tetpil;
+  long n,m,i,j,lM;
   GEN p1,delta,H,U,u1,u2,x;
 
   if (typ(M)!=t_MAT) err(typeer,"gaussmodulo");
@@ -2950,7 +3269,7 @@ gaussmoduloall(GEN M, GEN D, GEN Y, GEN *ptu1)
 GEN
 matsolvemod0(GEN M, GEN D, GEN Y, long flag)
 {
-  long av;
+  gpmem_t av;
   GEN p1,y;
 
   if (!flag) return gaussmoduloall(M,D,Y,NULL);