===================================================================
RCS file: /home/cvs/OpenXM_contrib/pari-2.2/src/basemath/Attic/base4.c,v
retrieving revision 1.1
retrieving revision 1.2
diff -u -p -r1.1 -r1.2
--- OpenXM_contrib/pari-2.2/src/basemath/Attic/base4.c	2001/10/02 11:17:02	1.1
+++ OpenXM_contrib/pari-2.2/src/basemath/Attic/base4.c	2002/09/11 07:26:49	1.2
@@ -1,4 +1,4 @@
-/* $Id: base4.c,v 1.1 2001/10/02 11:17:02 noro Exp $
+/* $Id: base4.c,v 1.2 2002/09/11 07:26:49 noro Exp $
 
 Copyright (C) 2000  The PARI group.
 
@@ -22,13 +22,18 @@ Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 #include "pari.h"
 #include "parinf.h"
 
-#define principalideal_aux(nf,x) (principalideal0((nf),(x),0))
+extern GEN makeprimetoideal(GEN nf,GEN UV,GEN uv,GEN x);
+extern GEN gauss_triangle_i(GEN A, GEN B,GEN t);
+extern GEN hnf_invimage(GEN A, GEN b);
+extern int hnfdivide(GEN A, GEN B);
+extern GEN colreducemodHNF(GEN x, GEN y, GEN *Q);
+extern GEN zinternallog_pk(GEN nf,GEN a0,GEN y,GEN pr,GEN prk,GEN list,GEN *psigne);
+extern GEN special_anti_uniformizer(GEN nf, GEN pr);
+extern GEN set_sign_mod_idele(GEN nf, GEN x, GEN y, GEN idele, GEN sarch);
+extern long int_elt_val(GEN nf, GEN x, GEN p, GEN b, GEN *newx);
 
-extern GEN element_muli(GEN nf, GEN x, GEN y);
-extern GEN colreducemodmat(GEN x, GEN y, GEN *Q);
+static GEN mat_ideal_two_elt2(GEN nf, GEN x, GEN a);
 
-static GEN nfbezout(GEN nf, GEN a, GEN b, GEN ida, GEN idb, GEN *u, GEN *v, GEN *w, GEN *di);
-
 /*******************************************************************/
 /*                                                                 */
 /*                     IDEAL OPERATIONS                            */
@@ -47,8 +52,7 @@ static GEN nfbezout(GEN nf, GEN a, GEN b, GEN ida, GEN
  * (depends on the chosen generator a). All subroutines work with either
  * ideles or ideals (an omitted V is assumed to be 0).
  *
- * All the output ideals will be in HNF form.
- */
+ * All ideals are output in HNF form. */
 
 /* types and conversions */
 
@@ -85,38 +89,18 @@ idealtyp(GEN *ideal, GEN *arch)
   *ideal = x; return t;
 }
 
-/* Assume ideal in HNF form */
-long
-ideal_is_zk(GEN ideal,long N)
-{
-  long i,j, lx = lg(ideal);
-
-  if (typ(ideal) != t_MAT || lx==1) return 0;
-  N++; if (lx != N || lg(ideal[1]) != N) return 0;
-  for (i=1; i<N; i++)
-  {
-    if (!gcmp1(gcoeff(ideal,i,i))) return 0;
-    for (j=i+1; j<N; j++)
-      if (!gcmp0(gcoeff(ideal,i,j))) return 0;
-  }
-  return 1;
-}
-
 static GEN
 prime_to_ideal_aux(GEN nf, GEN vp)
 {
-  GEN m,el;
-  long i, N = degpol(nf[1]);
-
-  m = cgetg(N+1,t_MAT); el = (GEN)vp[2];
-  for (i=1; i<=N; i++) m[i] = (long) element_mulid(nf,el,i);
-  return hnfmodid(m,(GEN)vp[1]);
+  GEN m = eltmul_get_table(nf, (GEN)vp[2]);
+  return hnfmodid(m, (GEN)vp[1]);
 }
 
+/* vp = [a,x,...], a in Z. Return (a,x)  [HACK: vp need not be prime] */
 GEN
 prime_to_ideal(GEN nf, GEN vp)
 {
-  long av=avma;
+  gpmem_t av=avma;
   if (typ(vp) == t_INT) return gscalmat(vp, degpol(nf[1]));
   return gerepileupto(av, prime_to_ideal_aux(nf,vp));
 }
@@ -126,12 +110,12 @@ static GEN
 idealmat_to_hnf(GEN nf, GEN x)
 {
   long rx,i,j,N;
-  GEN m,dx;
+  GEN m, cx;
 
   N=degpol(nf[1]); rx=lg(x)-1;
   if (!rx) return gscalmat(gzero,N);
 
-  dx=denom(x); if (gcmp1(dx)) dx = NULL; else x=gmul(dx,x);
+  x = Q_primitive_part(x, &cx);
   if (rx >= N) m = x;
   else
   {
@@ -140,8 +124,8 @@ idealmat_to_hnf(GEN nf, GEN x)
       for (j=1; j<=N; j++)
         m[(i-1)*N + j] = (long) element_mulid(nf,(GEN)x[i],j);
   }
-  x = hnfmod(m,detint(m));
-  return dx? gdiv(x,dx): x;
+  x = hnfmod(m, detint(m));
+  return cx? gmul(x,cx): x;
 }
 
 int
@@ -156,18 +140,31 @@ ishnfall(GEN x)
   }
   return (gsigne(gcoeff(x,1,1)) > 0);
 }
+/* same x is known to be integral */
+int
+Z_ishnfall(GEN x)
+{
+  long i,j, lx = lg(x);
+  for (i=2; i<lx; i++)
+  {
+    if (signe(gcoeff(x,i,i)) <= 0) return 0;
+    for (j=1; j<i; j++)
+      if (signe(gcoeff(x,i,j))) return 0;
+  }
+  return (signe(gcoeff(x,1,1)) > 0);
+}
 
 GEN
 idealhermite_aux(GEN nf, GEN x)
 {
   long N,tx,lx;
-  GEN z;
+  GEN z, cx;
 
   tx = idealtyp(&x,&z);
   if (tx == id_PRIME) return prime_to_ideal_aux(nf,x);
   if (tx == id_PRINCIPAL)
   {
-    x = principalideal(nf,x);
+    x = eltmul_get_table(nf, x);
     return idealmat_to_hnf(nf,x);
   }
   N=degpol(nf[1]); lx = lg(x);
@@ -175,127 +172,212 @@ idealhermite_aux(GEN nf, GEN x)
 
   if (lx == N+1 && ishnfall(x)) return x;
   if (lx <= N) return idealmat_to_hnf(nf,x);
-  z=denom(x); if (gcmp1(z)) z=NULL; else x = gmul(z,x);
+  x = Q_primitive_part(x, &cx);
   x = hnfmod(x,detint(x));
-  return z? gdiv(x,z): x;
+  return cx? gmul(x, cx): x;
 }
 
 GEN
 idealhermite(GEN nf, GEN x)
 {
-  long av=avma;
+  gpmem_t av=avma;
   GEN p1;
   nf = checknf(nf); p1 = idealhermite_aux(nf,x);
   if (p1==x || p1==(GEN)x[1]) return gcopy(p1);
   return gerepileupto(av,p1);
 }
 
-static GEN
-principalideal0(GEN nf, GEN x, long copy)
+GEN
+principalideal(GEN nf, GEN x)
 {
   GEN z = cgetg(2,t_MAT);
+
+  nf = checknf(nf);
   switch(typ(x))
   {
     case t_INT: case t_FRAC: case t_FRACN:
-      if (copy) x = gcopy(x);
-      x = gscalcol_i(x, degpol(nf[1])); break;
+      x = gscalcol(x, degpol(nf[1])); break;
 
     case t_POLMOD:
       x = checknfelt_mod(nf,x,"principalideal");
       /* fall through */
     case t_POL:
-      x = copy? algtobasis(nf,x): algtobasis_intern(nf,x);
+      x = algtobasis(nf,x);
       break;
 
     case t_MAT:
       if (lg(x)!=2) err(typeer,"principalideal");
       x = (GEN)x[1];
     case t_COL:
-      if (lg(x)==lgef(nf[1])-2)
-      {
-        if (copy) x = gcopy(x);
-        break;
-      }
+      if (lg(x)-1==degpol(nf[1])) { x = gcopy(x); break; }
     default: err(typeer,"principalideal");
   }
   z[1]=(long)x; return z;
 }
 
-GEN
-principalideal(GEN nf, GEN x)
+static GEN
+mylog(GEN x, long prec)
 {
-  nf = checknf(nf); return principalideal0(nf,x,1);
+  if (gcmp0(x)) err(precer,"get_arch");
+  return glog(x,prec);
 }
 
+GEN get_arch(GEN nf,GEN x,long prec);
+
 static GEN
-mylog(GEN x, long prec)
+famat_get_arch(GEN nf, GEN x, long prec)
 {
-  if (gcmp0(x))
-    err(precer,"get_arch");
-  return glog(x,prec);
+  GEN A, a, g = (GEN)x[1], e = (GEN)x[2];
+  long i, l = lg(e);
+
+  if (l <= 1)
+    return get_arch(nf, gun, prec);
+  A = NULL; /* -Wall */
+  for (i=1; i<l; i++)
+  {
+    a = get_arch(nf, (GEN)g[i], prec);
+    a = gmul((GEN)e[i], a);
+    A = (i == 1)? a: gadd(A, a);
+  }
+  return A;
 }
 
-/* for internal use */
+static GEN
+scalar_get_arch(long R1, long RU, GEN x, long prec)
+{
+  GEN v = cgetg(RU+1,t_VEC);
+  GEN p1 = glog(x, prec);
+  long i;
+
+  for (i=1; i<=R1; i++) v[i] = (long)p1;
+  if (i <= RU)
+  {
+    p1 = gmul2n(p1,1);
+    for (   ; i<=RU; i++) v[i] = (long)p1;
+  }
+  return v;
+}
+
+/* For internal use. Get archimedean components: [e_i log( sigma_i(x) )],
+ * with e_i = 1 (resp 2.) for i <= R1 (resp. > R1) */
 GEN
 get_arch(GEN nf,GEN x,long prec)
 {
-  long i,R1,RU;
-  GEN v,p1,p2;
+  long i, RU, R1 = nf_get_r1(nf);
+  GEN v;
 
-  R1=itos(gmael(nf,2,1)); RU = R1+itos(gmael(nf,2,2));
-  if (typ(x)!=t_COL) x = algtobasis_intern(nf,x);
+  RU = lg(nf[6]) - 1;
+  switch(typ(x))
+  {
+    case t_MAT: return famat_get_arch(nf,x,prec);
+    case t_POLMOD:
+    case t_POL: x = algtobasis_i(nf,x);   /* fall through */
+    case t_COL: if (!isnfscalar(x)) break;
+      x = (GEN)x[1]; /* fall through */
+    default: /* rational number */
+      return scalar_get_arch(R1, RU, x, prec);
+  }
+
   v = cgetg(RU+1,t_VEC);
   if (isnfscalar(x)) /* rational number */
   {
-    p1 = glog((GEN)x[1],prec);
-    p2 = (RU > R1)? gmul2n(p1,1): NULL;
-    for (i=1; i<=R1; i++) v[i]=(long)p1;
-    for (   ; i<=RU; i++) v[i]=(long)p2;
+    GEN p1 = glog((GEN)x[1], prec);
+
+    for (i=1; i<=R1; i++) v[i] = (long)p1;
+    if (i <= RU)
+    {
+      p1 = gmul2n(p1,1);
+      for (   ; i<=RU; i++) v[i] = (long)p1;
+    }
   }
-  else
+  x = gmul(gmael(nf,5,1),x);
+  for (i=1; i<=R1; i++) v[i] = (long)mylog((GEN)x[i],prec);
+  for (   ; i<=RU; i++) v[i] = lmul2n(mylog((GEN)x[i],prec),1);
+  return v;
+}
+
+GEN get_arch_real(GEN nf,GEN x,GEN *emb,long prec);
+
+static GEN
+famat_get_arch_real(GEN nf,GEN x,GEN *emb,long prec)
+{
+  GEN A, T, a, t, g = (GEN)x[1], e = (GEN)x[2];
+  long i, l = lg(e);
+
+  if (l <= 1)
+    return get_arch_real(nf, gun, emb, prec);
+  A = T = NULL; /* -Wall */
+  for (i=1; i<l; i++)
   {
-    x = gmul(gmael(nf,5,1),x);
-    for (i=1; i<=R1; i++) v[i] = (long)mylog((GEN)x[i],prec);
-    for (   ; i<=RU; i++) v[i] = lmul2n(mylog((GEN)x[i],prec),1);
+    a = get_arch_real(nf, (GEN)g[i], &t, prec);
+    if (!a) return NULL;
+    a = gmul((GEN)e[i], a);
+    t = vecpow(t, (GEN)e[i]);
+    if (i == 1) { A = a;          T = t; }
+    else        { A = gadd(A, a); T = vecmul(T, t); }
   }
-  return v;
+  *emb = T; return A;
 }
 
+static GEN
+scalar_get_arch_real(long R1, long RU, GEN u, GEN *emb, long prec)
+{
+  GEN v, x, p1;
+  long i, s;
+
+  s = gsigne(u);
+  if (!s) err(talker,"0 in get_arch_real");
+  x = cgetg(RU+1, t_COL);
+  for (i=1; i<=RU; i++) x[i] = (long)u;
+
+  v = cgetg(RU+1, t_COL);
+  p1 = (s > 0)? glog(u,prec): gzero;
+  for (i=1; i<=R1; i++) v[i] = (long)p1;
+  if (i <= RU)
+  {
+    p1 = gmul2n(p1,1);
+    for (   ; i<=RU; i++) v[i] = (long)p1;
+  }
+  *emb = x; return v;
+}
+
+static int
+low_prec(GEN x)
+{
+  return gcmp0(x) || (typ(x) == t_REAL && lg(x) == 3);
+}
+
 /* as above but return NULL if precision problem, and set *emb to the
  * embeddings of x */
 GEN
-get_arch_real(GEN nf,GEN x,GEN *emb,long prec)
+get_arch_real(GEN nf, GEN x, GEN *emb, long prec)
 {
-  long i,R1,RU;
-  GEN v,p1,p2;
+  long i, RU, R1 = nf_get_r1(nf);
+  GEN v, t;
 
-  R1=itos(gmael(nf,2,1)); RU = R1+itos(gmael(nf,2,2));
-  if (typ(x)!=t_COL) x = algtobasis_intern(nf,x);
+  RU = lg(nf[6])-1;
+  switch(typ(x))
+  {
+    case t_MAT: return famat_get_arch_real(nf,x,emb,prec);
+
+    case t_POLMOD:
+    case t_POL: x = algtobasis_i(nf,x);   /* fall through */
+    case t_COL: if (!isnfscalar(x)) break;
+      x = (GEN)x[1]; /* fall through */
+    default: /* rational number */
+      return scalar_get_arch_real(R1, RU, x, emb, prec);
+  }
   v = cgetg(RU+1,t_COL);
-  if (isnfscalar(x)) /* rational number */
+  x = gmul(gmael(nf,5,1), x);
+  for (i=1; i<=R1; i++)
   {
-    GEN u = (GEN)x[1];
-    i = signe(u);
-    if (!i) err(talker,"0 in get_arch_real");
-    p1= (i > 0)? glog(u,prec): gzero;
-    p2 = (RU > R1)? gmul2n(p1,1): NULL;
-    for (i=1; i<=R1; i++) v[i]=(long)p1;
-    for (   ; i<=RU; i++) v[i]=(long)p2;
+    t = gabs((GEN)x[i],prec); if (low_prec(t)) return NULL;
+    v[i] = llog(t,prec);
   }
-  else
+  for (   ; i<=RU; i++)
   {
-    GEN t;
-    x = gmul(gmael(nf,5,1),x);
-    for (i=1; i<=R1; i++)
-    {
-      t = gabs((GEN)x[i],prec); if (gcmp0(t)) return NULL;
-      v[i] = llog(t,prec);
-    }
-    for (   ; i<=RU; i++)
-    {
-      t = gnorm((GEN)x[i]); if (gcmp0(t)) return NULL;
-      v[i] = llog(t,prec);
-    }
+    t = gnorm((GEN)x[i]); if (low_prec(t)) return NULL;
+    v[i] = llog(t,prec);
   }
   *emb = x; return v;
 }
@@ -303,13 +385,12 @@ get_arch_real(GEN nf,GEN x,GEN *emb,long prec)
 GEN
 principalidele(GEN nf, GEN x, long prec)
 {
-  GEN p1,y = cgetg(3,t_VEC);
-  long av;
+  GEN p1, y = cgetg(3,t_VEC);
+  gpmem_t av;
 
-  nf = checknf(nf);
-  p1 = principalideal0(nf,x,1);
+  p1 = principalideal(nf,x);
   y[1] = (long)p1;
-  av =avma; p1 = get_arch(nf,(GEN)p1[1],prec);
+  av = avma; p1 = get_arch(nf,(GEN)p1[1],prec);
   y[2] = lpileupto(av,p1); return y;
 }
 
@@ -355,26 +436,17 @@ idealaddtoone0(GEN nf, GEN arg1, GEN arg2)
   return idealaddtoone(nf,arg1,arg2);
 }
 
-static GEN
-two_to_hnf(GEN nf, GEN a, GEN b)
-{
-  a = principalideal_aux(nf,a);
-  b = principalideal_aux(nf,b);
-  a = concatsp(a,b);
-  if (lgef(nf[1])==5) /* quadratic field: a has to be turned into idealmat */
-    a = idealmul(nf,idmat(2),a);
-  return idealmat_to_hnf(nf, a);
-}
-
 GEN
 idealhnf0(GEN nf, GEN a, GEN b)
 {
-  long av;
+  gpmem_t av;
+  GEN x;
   if (!b) return idealhermite(nf,a);
 
   /* HNF of aZ_K+bZ_K */
-  av = avma; nf=checknf(nf);
-  return gerepileupto(av, two_to_hnf(nf,a,b));
+  av = avma; nf = checknf(nf);
+  x = concatsp(eltmul_get_table(nf,a), eltmul_get_table(nf,b));
+  return gerepileupto(av, idealmat_to_hnf(nf, x));
 }
 
 GEN
@@ -386,59 +458,11 @@ idealhermite2(GEN nf, GEN a, GEN b)
 static int
 ok_elt(GEN x, GEN xZ, GEN y)
 {
-  ulong av = avma;
+  gpmem_t av = avma;
   int r = gegal(x, hnfmodid(y, xZ));
   avma = av; return r;
 }
 
-static void
-setprec(GEN x, long prec)
-{
-  long i,j, n=lg(x);
-  for (i=1;i<n;i++)
-  {
-    GEN p2,p1 = (GEN)x[i];
-    for (j=1;j<n;j++)
-    {
-      p2 = (GEN)p1[j];
-      if (typ(p2) == t_REAL) setlg(p2, prec);
-    }
-  }
-}
-
-/* find a basis of x whose elements have small norm
- * M a bound for the size of coeffs of x */
-GEN
-ideal_better_basis(GEN nf, GEN x, GEN M)
-{
-  GEN a,b;
-  long nfprec = nfgetprec(nf);
-  long prec = DEFAULTPREC + (expi(M) >> TWOPOTBITS_IN_LONG);
-
-  if (typ(nf[5]) != t_VEC) return x;
-  if ((prec<<1) < nfprec) prec = (prec+nfprec) >> 1;
-  a = qf_base_change(gmael(nf,5,3),x,1);
-  setprec(a,prec);
-  b = lllgramintern(a,4,1,prec);
-  if (!b)
-  {
-    if (DEBUGLEVEL)
-      err(warner, "precision too low in ideal_better_basis (1)");
-    if (nfprec > prec)
-    {
-      setprec(a,nfprec);
-      b = lllgramintern(a,4,1,nfprec);
-    }
-  }
-  if (!b)
-  {
-    if (DEBUGLEVEL)
-      err(warner, "precision too low in ideal_better_basis (2)");
-    b = lllint(x); /* better than nothing */
-  }
-  return gmul(x, b);
-}
-
 static GEN
 addmul_col(GEN a, long s, GEN b)
 {
@@ -468,53 +492,58 @@ addmul_mat(GEN a, long s, GEN b)
 static GEN
 mat_ideal_two_elt(GEN nf, GEN x)
 {
-  GEN y,a,beta,cx,xZ,mul, pol = (GEN)nf[1];
-  long i,j,lm, N = degpol(pol);
-  ulong av,tetpil;
+  GEN y,a,beta,cx,xZ,mul;
+  long i,lm, N = degpol(nf[1]);
+  gpmem_t av0,av,tetpil;
 
-  y=cgetg(3,t_VEC); av=avma;
-  if (lg(x[1])!=N+1) err(typeer,"ideal_two_elt");
+  y = cgetg(3,t_VEC); av = avma;
+  if (lg(x[1]) != N+1) err(typeer,"ideal_two_elt");
   if (N == 2)
   {
     y[1] = lcopy(gcoeff(x,1,1));
     y[2] = lcopy((GEN)x[2]); return y;
   }
 
-  cx = content(x); if (!gcmp1(cx)) x = gdiv(x,cx);
+  x = Q_primitive_part(x, &cx); if (!cx) cx = gun;
   if (lg(x) != N+1) x = idealhermite_aux(nf,x);
+
   xZ = gcoeff(x,1,1);
   if (gcmp1(xZ))
   {
     y[1] = lpilecopy(av,cx);
-    y[2] = (long)gscalcol(cx,N); return y;
+    y[2] = (long)gscalcol(cx, N); return y;
   }
+  av0 = avma;
   a = NULL; /* gcc -Wall */
   beta= cgetg(N+1, t_VEC);
   mul = cgetg(N+1, t_VEC); lm = 1; /* = lg(mul) */
   /* look for a in x such that a O/xZ = x O/xZ */
   for (i=2; i<=N; i++)
   {
-    GEN t, y = cgetg(N+1,t_MAT);
-    a = (GEN)x[i];
-    for (j=1; j<=N; j++) y[j] = (long)element_mulid(nf,a,j);
-    /* columns of mul[i] = canonical generators for x[i] O/xZ as Z-module */
-    t = gmod(y, xZ);
-    if (gcmp0(t)) continue;
-    if (ok_elt(x,xZ, t)) break;
+    gpmem_t av1 = avma;
+    GEN t, y = eltmul_get_table(nf, (GEN)x[i]);
+    t = FpM_red(y, xZ);
+    if (gcmp0(t)) { avma = av1; continue; }
+    if (ok_elt(x,xZ, t)) { a = (GEN)x[i]; break; }
     beta[lm]= x[i];
+    /* mul[i] = { canonical generators for x[i] O/xZ as Z-module } */
     mul[lm] = (long)t; lm++;
   }
-  if (i>N)
+  if (i > N)
   {
     GEN z = cgetg(lm, t_VECSMALL);
-    ulong av1, c = 0;
+    gpmem_t av1;
+    ulong c = 0;
 
     setlg(mul, lm);
     setlg(beta,lm);
-    if (DEBUGLEVEL>3) fprintferr("ideal_two_elt, hard case: ");
+    if (DEBUGLEVEL>3) fprintferr("ideal_two_elt, hard case:\n");
     for(av1=avma;;avma=av1)
     {
-      if (DEBUGLEVEL>3) fprintferr("%ld ", ++c);
+      if (++c == 100) {
+        if (DEBUGLEVEL>3) fprintferr("using approximation theorem\n");
+        a = mat_ideal_two_elt2(nf, x, xZ); goto END;
+      }
       for (a=NULL,i=1; i<lm; i++)
       {
         long t = (mymyrand() >> (BITS_IN_RANDOM-5)) - 7; /* in [-7,8] */
@@ -526,19 +555,19 @@ mat_ideal_two_elt(GEN nf, GEN x)
     }
     for (a=NULL,i=1; i<lm; i++)
       a = addmul_col(a, z[i], (GEN)beta[i]);
-    if (DEBUGLEVEL>3) fprintferr("\n");
   }
-  a = centermod(a, xZ); tetpil=avma;
+END:
+  a = centermod(a, xZ);
+  tetpil = avma;
   y[1] = lmul(xZ,cx);
   y[2] = lmul(a, cx);
   gerepilemanyvec(av,tetpil,y+1,2); return y;
 }
 
-/* Etant donne un ideal ix, ressort un vecteur [a,alpha] a deux composantes
- * tel que a soit rationnel et ix=aZ_K+alpha Z_K, alpha etant un vecteur
- * colonne sur la base d'entiers. On peut avoir a=0 ou alpha=0, mais on ne
- * cherche pas a determiner si ix est principal.
- */
+/* Given an ideal x, returns [a,alpha] such that a is in Q,
+ * x = a Z_K + alpha Z_K, alpha in K^*
+ * a = 0 or alpha = 0 are possible, but do not try to determine whether
+ * x is principal. */
 GEN
 ideal_two_elt(GEN nf, GEN x)
 {
@@ -604,30 +633,33 @@ factor_norm(GEN x)
 GEN
 idealfactor(GEN nf, GEN x)
 {
-  long av,tx, tetpil,i,j,k,lf,lc,N,l,v,vc,e;
+  gpmem_t av,tetpil;
+  long tx,i,j,k,lf,lc,N,l,v,vc,e;
   GEN f,f1,f2,c1,c2,y1,y2,y,p1,p2,cx,P;
 
   tx = idealtyp(&x,&y);
   if (tx == id_PRIME)
   {
     y=cgetg(3,t_MAT);
-    y[1]=lgetg(2,t_COL); mael(y,1,1)=lcopy(x);
-    y[2]=lgetg(2,t_COL); mael(y,2,1)=un; return y;
+    y[1]=(long)_col(gcopy(x));
+    y[2]=(long)_col(gun); return y;
   }
-  nf=checknf(nf); av=avma;
-  if (tx == id_PRINCIPAL) x = principalideal_aux(nf,x);
-
-  N=degpol(nf[1]); if (lg(x) != N+1) x = idealmat_to_hnf(nf,x);
+  av = avma;
+  nf = checknf(nf);
+  N = degpol(nf[1]);
+  if (tx == id_PRINCIPAL) x = idealhermite_aux(nf, x);
+  else
+    if (lg(x) != N+1) x = idealmat_to_hnf(nf,x);
   if (lg(x)==1) err(talker,"zero ideal in idealfactor");
-  cx = content(x);
-  if (gcmp1(cx))
+  x = Q_primitive_part(x, &cx);
+  if (!cx)
   {
     c1 = c2 = NULL; /* gcc -Wall */
-    lc=1;
+    lc = 1;
   }
   else
   {
-    f = factor(cx); x = gdiv(x,cx);
+    f = factor(cx);
     c1 = (GEN)f[1];
     c2 = (GEN)f[2]; lc = lg(c1);
   }
@@ -641,7 +673,7 @@ idealfactor(GEN nf, GEN x)
   {
     p1 = primedec(nf,(GEN)f1[i]);
     l = f2[i]; /* = v_p(Nx) */
-    vc = ggval(cx,(GEN)f1[i]);
+    vc = cx? ggval(cx,(GEN)f1[i]): 0;
     for (j=1; j<lg(p1); j++)
     {
       P = (GEN)p1[j]; e = itos((GEN)P[3]);
@@ -685,14 +717,15 @@ idealfactor(GEN nf, GEN x)
 long
 idealval(GEN nf, GEN ix, GEN P)
 {
-  long N,v,vd,w,av=avma,av1,lim,e,f,i,j,k, tx = typ(ix);
+  gpmem_t av=avma,av1,lim;
+  long N,v,vd,w,e,f,i,j,k, tx = typ(ix);
   GEN mul,mat,a,x,y,r,bp,p,pk,cx;
 
   nf=checknf(nf); checkprimeid(P);
   if (is_extscalar_t(tx) || tx==t_COL) return element_val(nf,ix,P);
   p=(GEN)P[1]; N=degpol(nf[1]);
   tx = idealtyp(&ix,&a);
-  cx = content(ix); if (!gcmp1(cx)) ix = gdiv(ix,cx);
+  ix = Q_primitive_part(ix, &cx);
   if (tx != id_MAT)
     ix = idealhermite_aux(nf,ix);
   else
@@ -707,7 +740,7 @@ idealval(GEN nf, GEN ix, GEN P)
   /* 0 <= ceil[v_P(ix) / e] <= v_p(ix \cap Z) --> v_P <= e * v_p */
   j = k * e;
   v = min(i,j); /* v_P(ix) <= v */
-  vd = ggval(cx,p) * e;
+  vd = cx? ggval(cx,p) * e: 0;
   if (!v) { avma = av; return vd; }
 
   mul = cgetg(N+1,t_MAT); bp=(GEN)P[5];
@@ -762,7 +795,8 @@ idealval(GEN nf, GEN ix, GEN P)
 GEN
 idealadd(GEN nf, GEN x, GEN y)
 {
-  long av=avma,N,tx,ty;
+  gpmem_t av=avma;
+  long N,tx,ty;
   GEN z,p1,dx,dy,dz;
   int modid;
 
@@ -774,24 +808,27 @@ idealadd(GEN nf, GEN x, GEN y)
   if (ty != id_MAT || lg(y)!=N+1) y = idealhermite_aux(nf,y);
   if (lg(x) == 1) return gerepileupto(av,y);
   if (lg(y) == 1) return gerepileupto(av,x); /* check for 0 ideal */
-  dx=denom(x);
-  dy=denom(y); dz=mulii(dx,dy);
-  if (gcmp1(dz)) dz = NULL; else { x=gmul(x,dz); y=gmul(y,dz); }
+  dx = denom(x);
+  dy = denom(y); dz = mulii(dx,dy);
+  if (gcmp1(dz)) dz = NULL; else {
+    x = Q_muli_to_int(x, dz);
+    y = Q_muli_to_int(y, dz);
+  }
   if (isnfscalar((GEN)x[1]) && isnfscalar((GEN)y[1]))
   {
-    p1 = mppgcd(gcoeff(x,1,1),gcoeff(y,1,1));
+    p1 = mppgcd(gcoeff(x,1,1), gcoeff(y,1,1));
     modid = 1;
   }
   else
   {
-    p1 = mppgcd(detint(x),detint(y));
+    p1 = mppgcd(detint(x), detint(y));
     modid = 0;
   }
   if (gcmp1(p1))
   {
     long i,j;
     if (!dz) { avma=av; return idmat(N); }
-    avma = (long)dz; dz = gerepileupto((long)z, ginv(dz));
+    avma = (gpmem_t)dz; dz = gerepileupto((gpmem_t)z, ginv(dz));
     for (i=1; i<=N; i++)
     {
       z[i]=lgetg(N+1,t_COL);
@@ -806,34 +843,36 @@ idealadd(GEN nf, GEN x, GEN y)
   return gerepileupto(av,z);
 }
 
+/* Assume x,y integral non zero (presumably coprime, which is checked).
+ * Return a in x such that 1-a in y. If (x + y) \cap Z is 1, then addone_aux
+ * is more efficient. */
+GEN
+addone_aux2(GEN x, GEN y)
+{
+  GEN U, H;
+  long i, l;
+
+  H = hnfall_i(concatsp(x,y), &U, 1);
+  l = lg(H);
+  for (i=1; i<l; i++)
+    if (!gcmp1(gcoeff(H,i,i)))
+      err(talker,"ideals don't sum to Z_K in idealaddtoone");
+  U = (GEN)U[l]; setlg(U, lg(x));
+  return gmul(x, U);
+}
+
+/* assume x,y integral, non zero in HNF */
 static GEN
-get_p1(GEN nf, GEN x, GEN y,long fl)
+addone_aux(GEN x, GEN y)
 {
-  GEN u,v,v1,v2,v3,v4;
-  long i,j,N;
+  GEN a, b, u, v;
 
-  switch(fl)
-  {
-    case 1:
-      v1 = gcoeff(x,1,1);
-      v2 = gcoeff(y,1,1);
-      if (typ(v1)!=t_INT || typ(v2)!=t_INT)
-        err(talker,"ideals don't sum to Z_K in idealaddtoone");
-      if (gcmp1(bezout(v1,v2,&u,&v)))
-        return gmul(u,(GEN)x[1]);
-    default:
-      v=hnfperm(concatsp(x,y));
-      v1=(GEN)v[1]; v2=(GEN)v[2]; v3=(GEN)v[3];
-      j=0; N = degpol(nf[1]);
-      for (i=1; i<=N; i++)
-      {
-        if (!gcmp1(gcoeff(v1,i,i)))
-          err(talker,"ideals don't sum to Z_K in idealaddtoone");
-        if (gcmp1((GEN)v3[i])) j=i;
-      }
-      v4=(GEN)v2[N+j]; setlg(v4,N+1);
-      return gmul(x,v4);
-  }
+  a = gcoeff(x,1,1);
+  b = gcoeff(y,1,1);
+  if (typ(a) != t_INT || typ(b) != t_INT)
+    err(talker,"ideals don't sum to Z_K in idealaddtoone");
+  if (gcmp1(bezout(a,b,&u,&v))) return gmul(u,(GEN)x[1]);
+  return addone_aux2(x, y);
 }
 
 GEN
@@ -849,15 +888,15 @@ idealaddtoone_i(GEN nf, GEN x, GEN y)
     fprintferr(" y = %Z\n",y);
   }
   t = idealtyp(&x,&p1);
-  if (t != id_MAT || lg(x) > 1 || lg(x) != lg(x[1]))
+  if (t != id_MAT || lg(x) == 1 || lg(x) != lg(x[1]))
     xh = idealhermite_aux(nf,x);
   else
-    { xh=x; fl = isnfscalar((GEN)x[1]); }
+    { xh = x; fl = isnfscalar((GEN)x[1]); }
   t = idealtyp(&y,&p1);
   if (t != id_MAT || lg(y) == 1 || lg(y) != lg(y[1]))
     yh = idealhermite_aux(nf,y);
   else
-    { yh=y; if (fl) fl = isnfscalar((GEN)y[1]); }
+    { yh = y; if (fl) fl = isnfscalar((GEN)y[1]); }
   if (lg(xh) == 1)
   {
     if (lg(yh) == 1 || !gcmp1(gabs(gcoeff(yh,1,1),0)))
@@ -872,7 +911,8 @@ idealaddtoone_i(GEN nf, GEN x, GEN y)
     return algtobasis(nf, gneg(p1));
   }
 
-  p1 = get_p1(nf,xh,yh,fl);
+  if (fl) p1 = addone_aux (xh,yh); /* xh[1], yh[1] scalar */
+  else    p1 = addone_aux2(xh,yh);
   p1 = element_reduce(nf,p1, idealmullll(nf,x,y));
   if (DEBUGLEVEL>4 && !gcmp0(p1))
     fprintferr(" leaving idealaddtoone: %Z\n",p1);
@@ -909,7 +949,7 @@ ideleaddone_aux(GEN nf,GEN x,GEN ideal)
   return nba? p3: gcopy(p3);
 }
 
-static GEN
+GEN
 unnf_minus_x(GEN x)
 {
   long i, N = lg(x);
@@ -923,11 +963,23 @@ unnf_minus_x(GEN x)
 static GEN
 addone(GEN f(GEN,GEN,GEN), GEN nf, GEN x, GEN y)
 {
-  GEN z = cgetg(3,t_VEC);
-  long av=avma;
+  GEN z = cgetg(3,t_VEC), a;
+  gpmem_t av = avma;
+  nf = checknf(nf);
+  a = gerepileupto(av, f(nf,x,y));
+  z[1]=(long)a;
+  z[2]=(long)unnf_minus_x(a); return z;
+}
 
-  nf=checknf(nf); x = gerepileupto(av, f(nf,x,y));
-  z[1]=(long)x; z[2]=(long)unnf_minus_x(x); return z;
+/* assume x,y HNF, non-zero */
+static GEN
+addone_nored(GEN x, GEN y)
+{
+  GEN z = cgetg(3,t_VEC), a;
+  gpmem_t av = avma;
+  a = gerepileupto(av, addone_aux(x,y));
+  z[1] = (long)a;
+  z[2] = (long)unnf_minus_x(a); return z;
 }
 
 GEN
@@ -942,51 +994,17 @@ ideleaddone(GEN nf,GEN x,GEN idele)
   return addone(ideleaddone_aux,nf,x,idele);
 }
 
-/* return integral x = 0 mod p/pr^e, (x,pr) = 1.
- * Don't reduce mod p here: caller may need result mod pr^k */
-GEN
-special_anti_uniformizer(GEN nf, GEN pr)
-{
-  GEN p = (GEN)pr[1], e = (GEN)pr[3];
-  return gdivexact(element_pow(nf,(GEN)pr[5],e), gpuigs(p,itos(e)-1));
-}
-
-GEN
-nfmodprinit(GEN nf, GEN pr)
-{
-  long av;
-  GEN p,p1,prhall;
-
-  nf = checknf(nf); checkprimeid(pr);
-  prhall = cgetg(3,t_VEC);
-  prhall[1] = (long) prime_to_ideal(nf,pr);
-
-  av = avma; p = (GEN)pr[1];
-  p1 = cgetg(2,t_MAT);
-  p1[1] = (long)gmod(special_anti_uniformizer(nf, pr), p);
-  p1 = hnfmodid(idealhermite_aux(nf,p1), p);
-  p1 = idealaddtoone_i(nf,pr,p1);
-
-  /* p1 = 1 mod pr, p1 = 0 mod q^{e_q} for all other primes q | p */
-  prhall[2] = lpileupto(av, unnf_minus_x(p1)); return prhall;
-}
-
 /* given an element x in Z_K and an integral ideal y with x, y coprime,
    outputs an element inverse of x modulo y */
 GEN
 element_invmodideal(GEN nf, GEN x, GEN y)
 {
-  long av=avma,N,i, fl = 1;
+  gpmem_t av=avma;
+  long N,i, fl = 1;
   GEN v,p1,xh,yh;
 
   nf=checknf(nf); N=degpol(nf[1]);
-  if (ideal_is_zk(y,N)) return zerocol(N);
-  if (DEBUGLEVEL>4)
-  {
-    fprintferr(" entree dans element_invmodideal() :\n");
-    fprintferr(" x = "); outerr(x);
-    fprintferr(" y = "); outerr(y);
-  }
+  if (gcmp1(gcoeff(y,1,1))) return zerocol(N);
   i = lg(y);
   if (typ(y)!=t_MAT || i==1 || i != lg(y[1])) yh=idealhermite_aux(nf,y);
   else
@@ -998,18 +1016,18 @@ element_invmodideal(GEN nf, GEN x, GEN y)
     default: err(typeer,"element_invmodideal");
       return NULL; /* not reached */
   }
-  p1 = get_p1(nf,xh,yh,fl);
+  if (fl) p1 = addone_aux (xh,yh);
+  else    p1 = addone_aux2(xh,yh);
   p1 = element_div(nf,p1,x);
   v = gerepileupto(av, nfreducemodideal(nf,p1,y));
-  if (DEBUGLEVEL>2)
-    { fprintferr(" sortie de element_invmodideal : v = "); outerr(v); }
   return v;
 }
 
 GEN
 idealaddmultoone(GEN nf, GEN list)
 {
-  long av=avma,tetpil,N,i,i1,j,k;
+  gpmem_t av=avma,tetpil;
+  long N,i,i1,j,k;
   GEN z,v,v1,v2,v3,p1;
 
   nf=checknf(nf); N=degpol(nf[1]);
@@ -1053,17 +1071,17 @@ idealaddmultoone(GEN nf, GEN list)
 /* multiplication */
 
 /* x integral ideal (without archimedean component) in HNF form
- * [a,alpha,n] corresponds to the ideal aZ_K+alpha Z_K (a is a
+ * y = [a,alpha] corresponds to the ideal aZ_K+alpha Z_K (a is a
  * rational integer). Multiply them
  */
 static GEN
-idealmulspec(GEN nf, GEN x, GEN a, GEN alpha)
+idealmulspec(GEN nf, GEN x, GEN y)
 {
   long i, N=lg(x)-1;
-  GEN m, mod;
+  GEN m, mod, a = (GEN)y[1], alpha = (GEN)y[2];
 
   if (isnfscalar(alpha))
-    return gmul(mppgcd(a,(GEN)alpha[1]),x);
+    return gmul(mppgcd(a, (GEN)alpha[1]),x);
   mod = mulii(a, gcoeff(x,1,1));
   m = cgetg((N<<1)+1,t_MAT);
   for (i=1; i<=N; i++) m[i]=(long)element_muli(nf,alpha,(GEN)x[i]);
@@ -1074,16 +1092,27 @@ idealmulspec(GEN nf, GEN x, GEN a, GEN alpha)
 /* x ideal (matrix form,maximal rank), vp prime ideal (primedec). Output the
  * product. Can be used for arbitrary vp of the form [p,a,e,f,b], IF vp
  * =pZ_K+aZ_K, p is an integer, and norm(vp) = p^f; e and b are not used.
- * For internal use.
- */
+ * For internal use. */
 GEN
 idealmulprime(GEN nf, GEN x, GEN vp)
 {
-  GEN denx = denom(x);
+  GEN cx;
+  x = Q_primitive_part(x, &cx);
+  x = idealmulspec(nf,x,vp);
+  return cx? gmul(x,cx): x;
+}
 
-  if (gcmp1(denx)) denx = NULL; else x = gmul(denx,x);
-  x = idealmulspec(nf,x, (GEN)vp[1], (GEN)vp[2]);
-  return denx? gdiv(x,denx): x;
+GEN arch_mul(GEN x, GEN y);
+GEN vecpow(GEN x, GEN n);
+GEN vecmul(GEN x, GEN y);
+
+static GEN
+mul(GEN nf, GEN x, GEN y)
+{
+  GEN yZ = gcoeff(y,1,1);
+  if (is_pm1(yZ)) return gcopy(x);
+  y = mat_ideal_two_elt(nf,y);
+  return idealmulspec(nf, x, y);
 }
 
 /* Assume ix and iy are integral in HNF form (or ideles of the same form).
@@ -1099,63 +1128,89 @@ idealmulh(GEN nf, GEN ix, GEN iy)
   if (typ(iy)==t_VEC && typ(iy[1]) != t_INT) {f+=2; y=(GEN)iy[1];} else y=iy;
   res = f? cgetg(3,t_VEC): NULL;
 
-  if (typ(y) != t_VEC) y = ideal_two_elt(nf,y);
-  y = idealmulspec(nf,x,(GEN)y[1],(GEN)y[2]);
+  if (typ(y) == t_VEC)
+    y = idealmulspec(nf,x,y);
+  else
+  {
+    GEN xZ = gcoeff(x,1,1);
+    GEN yZ = gcoeff(x,1,1);
+    y = (cmpii(xZ, yZ) < 0)? mul(nf,y,x): mul(nf,x,y);
+  }
   if (!f) return y;
 
-  res[1]=(long)y;
-  if (f==3) y = gadd((GEN)ix[2],(GEN)iy[2]);
+  res[1] = (long)y;
+  if (f==3) y = arch_mul((GEN)ix[2], (GEN)iy[2]);
   else
   {
     y = (f==2)? (GEN)iy[2]: (GEN)ix[2];
     y = gcopy(y);
   }
-  res[2]=(long)y; return res;
+  res[2] = (long)y; return res;
 }
 
+GEN
+mul_content(GEN cx, GEN cy)
+{
+  if (!cx) return cy;
+  if (!cy) return cx;
+  return gmul(cx,cy);
+}
+
 /* x and y are ideals in matrix form */
 static GEN
 idealmat_mul(GEN nf, GEN x, GEN y)
 {
-  long i,j, rx=lg(x)-1, ry=lg(y)-1;
-  GEN dx,dy,m;
+  long i,j, rx = lg(x)-1, ry = lg(y)-1;
+  GEN cx, cy, m;
 
-  dx=denom(x); if (!gcmp1(dx)) x=gmul(dx,x);
-  dy=denom(y); if (!gcmp1(dy)) y=gmul(dy,y);
-  dx = mulii(dx,dy);
-  if (rx<=2 || ry<=2)
+  x = Q_primitive_part(x, &cx);
+  y = Q_primitive_part(y, &cy); cx = mul_content(cx,cy);
+  if (rx <= 2 || ry <= 2)
   {
-    m=cgetg(rx*ry+1,t_MAT);
+    m = cgetg(rx*ry+1,t_MAT);
     for (i=1; i<=rx; i++)
       for (j=1; j<=ry; j++)
-        m[(i-1)*ry+j]=(long)element_muli(nf,(GEN)x[i],(GEN)y[j]);
+        m[(i-1)*ry+j] = (long)element_muli(nf,(GEN)x[i],(GEN)y[j]);
     if (isnfscalar((GEN)x[1]) && isnfscalar((GEN)y[1]))
     {
       GEN p1 = mulii(gcoeff(x,1,1),gcoeff(y,1,1));
       y = hnfmodid(m, p1);
     }
     else
-      y=hnfmod(m, detint(m));
+      y = hnfmod(m, detint(m));
   }
   else
   {
-    x=idealmat_to_hnf(nf,x);
-    y=idealmat_to_hnf(nf,y); y=idealmulh(nf,x,y);
+    if (!Z_ishnfall(x)) x = idealmat_to_hnf(nf,x);
+    if (!Z_ishnfall(y)) y = idealmat_to_hnf(nf,y);
+    y = idealmulh(nf,x,y);
   }
-  return gcmp1(dx)? y: gdiv(y,dx);
+  return cx? gmul(y,cx): y;
 }
 
-#if 0
-/* y is principal */
-static GEN
-add_arch(GEN nf, GEN ax, GEN y)
+/* operations on elements in factored form */
+
+GEN
+to_famat(GEN g, GEN e)
 {
-  long tetpil, av=avma, prec=precision(ax);
+  GEN h = cgetg(3,t_MAT);
+  if (typ(g) != t_COL) { g = dummycopy(g); settyp(g, t_COL); }
+  if (typ(e) != t_COL) { e = dummycopy(e); settyp(e, t_COL); }
+  h[1] = (long)g;
+  h[2] = (long)e; return h;
+}
 
-  y = get_arch(nf,y,prec); tetpil=avma;
-  return gerepile(av,tetpil,gadd(ax,y));
+GEN
+to_famat_all(GEN x, GEN y) { return to_famat(_col(x), _col(y)); }
+
+static GEN
+append(GEN v, GEN x)
+{
+  long i, l = lg(v);
+  GEN w = cgetg(l+1, typ(v));
+  for (i=1; i<l; i++) w[i] = lcopy((GEN)v[i]);
+  w[i] = lcopy(x); return w;
 }
-#endif
 
 /* add x^1 to factorisation f */
 static GEN
@@ -1169,14 +1224,14 @@ famat_add(GEN f, GEN x)
   }
   else
   {
-    h[1] = (long)concat((GEN)f[1], x);
+    h[1] = (long)append((GEN)f[1], x); /* x may be a t_COL */
     h[2] = (long)concat((GEN)f[2], gun);
   }
   return h;
 }
 
 /* cf merge_factor_i */
-static GEN
+GEN
 famat_mul(GEN f, GEN g)
 {
   GEN h;
@@ -1199,7 +1254,7 @@ famat_sqr(GEN f)
   h[2] = lmul2n((GEN)f[2],1);
   return h;
 }
-static GEN
+GEN
 famat_inv(GEN f)
 {
   GEN h;
@@ -1209,11 +1264,12 @@ famat_inv(GEN f)
   h[2] = lneg((GEN)f[2]);
   return h;
 }
-static GEN
+GEN
 famat_pow(GEN f, GEN n)
 {
   GEN h;
   if (lg(f) == 1) return cgetg(1,t_MAT);
+  if (typ(f) != t_MAT) return to_famat_all(f,n);
   h = cgetg(3,t_MAT);
   h[1] = lcopy((GEN)f[1]);
   h[2] = lmul((GEN)f[2],n);
@@ -1235,18 +1291,60 @@ famat_to_nf(GEN nf, GEN f)
   return t;
 }
 
-GEN
-to_famat(GEN g, GEN e)
+/* "compare" two nf elt. Goal is to quickly sort for uniqueness of
+ * representation, not uniqueness of represented element ! */
+static int
+elt_cmp(GEN x, GEN y)
 {
-  GEN h = cgetg(3,t_MAT);
-  h[1] = (long)g;
-  h[2] = (long)e; return h;
+  long tx = typ(x), ty = typ(y);
+  if (ty == tx)
+    return (tx == t_POL || tx == t_POLMOD)? cmp_pol(x,y): lexcmp(x,y);
+  return tx - ty;
 }
+static int
+elt_egal(GEN x, GEN y)
+{
+  if (typ(x) == typ(y)) return gegal(x,y);
+  return 0;
+}
 
 GEN
-to_famat_all(GEN x, GEN y) { return to_famat(_col(x), _col(y)); }
+famat_reduce(GEN fa)
+{
+  GEN E, F, G, L, g, e;
+  long i, k, l;
 
-/* assume (num(g[i]), id) = 1 and for all i. return prod g[i]^e[i] mod id */
+  if (lg(fa) == 1) return fa;
+  g = (GEN)fa[1]; l = lg(g);
+  e = (GEN)fa[2];
+  L = gen_sort(g, cmp_IND|cmp_C, &elt_cmp);
+  G = cgetg(l, t_COL);
+  E = cgetg(l, t_COL);
+  /* merge */
+  for (k=i=1; i<l; i++,k++) 
+  {
+    G[k] = g[L[i]];
+    E[k] = e[L[i]];
+    if (k > 1 && elt_egal((GEN)G[k], (GEN)G[k-1]))
+    {
+      E[k-1] = laddii((GEN)E[k], (GEN)E[k-1]);
+      k--;
+    }
+  }
+  /* kill 0 exponents */
+  l = k;
+  for (k=i=1; i<l; i++)
+    if (!gcmp0((GEN)E[i]))
+    {
+      G[k] = G[i];
+      E[k] = E[i]; k++;
+    }
+  F = cgetg(3, t_MAT);
+  setlg(G, k); F[1] = (long)G;
+  setlg(E, k); F[2] = (long)E; return F;
+}
+
+/* assume (num(g[i]), id) = 1 for all i. Return prod g[i]^e[i] mod id */
 GEN
 famat_to_nf_modideal_coprime(GEN nf, GEN g, GEN e, GEN id)
 {
@@ -1268,54 +1366,53 @@ famat_to_nf_modideal_coprime(GEN nf, GEN g, GEN e, GEN
   return t? t: gscalcol(gun, lg(id)-1);
 }
 
-/* assume prh has degree 1 and coprime to numerator(x) */
-GEN
-nf_to_Fp_simple(GEN x, GEN prh)
+/* assume pr has degree 1 and coprime to numerator(x) */
+static GEN
+nf_to_Fp_simple(GEN x, GEN modpr, GEN p)
 {
-  GEN ch = denom(x), p = gcoeff(prh,1,1);
-  if (!is_pm1(ch))
-  {
-    x = gmul(gmul(x,ch), mpinvmod(ch,p));
-  }
-  ch = colreducemodmat(gmod(x, p), prh, NULL);
-  return (GEN)ch[1]; /* in Fp^* */
+  GEN c, r = zk_to_ff(Q_primitive_part(x, &c), modpr);
+  if (c) r = gmod(gmul(r, c), p);
+  return r;
 }
 
-GEN
-famat_to_Fp_simple(GEN g, GEN e, GEN prh)
+static GEN
+famat_to_Fp_simple(GEN nf, GEN x, GEN modpr, GEN p)
 {
-  GEN t = gun, h,n, p = gcoeff(prh,1,1), q = subis(p,1);
-  long i, lx = lg(g);
-  for (i=1; i<lx; i++)
+  GEN h, n, t = gun, g = (GEN)x[1], e = (GEN)x[2], q = subis(p,1);
+  long i, l = lg(g);
+
+  for (i=1; i<l; i++)
   {
     n = (GEN)e[i]; n = modii(n,q);
     if (!signe(n)) continue;
-    h = nf_to_Fp_simple((GEN)g[i], prh);
+
+    h = (GEN)g[i];
+    switch(typ(h))
+    {
+      case t_POL: case t_POLMOD: h = algtobasis(nf, h);  /* fall through */
+      case t_COL: h = nf_to_Fp_simple(h, modpr, p); break;
+      default: h = gmod(h, p);
+    }
     t = mulii(t, powmodulo(h, n, p)); /* not worth reducing */
   }
   return modii(t, p);
 }
 
-/* cf famat_to_nf_modideal_coprime, but id is a prime of degree 1 (=prh) */
+/* cf famat_to_nf_modideal_coprime, but id is a prime of degree 1 (=pr) */
 GEN
-to_Fp_simple(GEN x, GEN prh)
+to_Fp_simple(GEN nf, GEN x, GEN pr)
 {
+  GEN T, p, modpr = zk_to_ff_init(nf, &pr, &T, &p);
   switch(typ(x))
   {
-    case t_COL: return nf_to_Fp_simple(x,prh);
-    case t_MAT: return famat_to_Fp_simple((GEN)x[1],(GEN)x[2],prh);
+    case t_COL: return nf_to_Fp_simple(x,modpr,p);
+    case t_MAT: return famat_to_Fp_simple(nf,x,modpr,p);
     default: err(impl,"generic conversion to finite field");
   }
   return NULL;
 }
 
-extern GEN zinternallog_pk(GEN nf,GEN a0,GEN y,GEN pr,GEN prk,GEN list,GEN *psigne);
-extern GEN colreducemodmat(GEN x, GEN y, GEN *Q);
-extern GEN special_anti_uniformizer(GEN nf, GEN pr);
-extern GEN set_sign_mod_idele(GEN nf, GEN x, GEN y, GEN idele, GEN sarch);
-extern long int_elt_val(GEN nf, GEN x, GEN p, GEN b, GEN *newx, long v);
-
-/* Compute t = prod g[i]^e[i] mod pr^n, assuming (t, pr) = 1.
+/* Compute t = prod g[i]^e[i] mod pr^k, assuming (t, pr) = 1.
  * Method: modify each g[i] so that it becomes coprime to pr :
  *  x / (p^k u) --> x * (b/p)^v_pr(x) / z^k u, where z = b^e/p^(e-1)
  * b/p = vp^(-1) times something prime to p; both numerator and denominator
@@ -1324,54 +1421,56 @@ extern long int_elt_val(GEN nf, GEN x, GEN p, GEN b, G
  * EX = exponent of (O_K / pr^k)^* used to reduce the product in case the
  * e[i] are large */
 static GEN
-famat_makecoprime(GEN nf, GEN g, GEN e, GEN pr, GEN prn, GEN EX)
+famat_makecoprime(GEN nf, GEN g, GEN e, GEN pr, GEN prk, GEN EX)
 {
-  long i,k, l = lg(g), N = degpol(nf[1]);
-  GEN prnZ,cx,x,u,z, zpow = gzero, p = (GEN)pr[1], b = (GEN)pr[5];
-  GEN mul = cgetg(N+1,t_MAT);
+  long i, l = lg(g);
+  GEN prkZ,cx,x,u, zpow = gzero, p = (GEN)pr[1], b = (GEN)pr[5];
+  GEN mul = eltmul_get_table(nf, b);
   GEN newg = cgetg(l+1, t_VEC); /* room for z */
 
-  prnZ = gcoeff(prn, 1,1);
-  z = gmod(special_anti_uniformizer(nf, pr), prnZ);
-  for (i=1; i<=N; i++) mul[i] = (long)element_mulid(nf,b,i);
+  prkZ = gcoeff(prk, 1,1);
   for (i=1; i < l; i++)
   {
     x = (GEN)g[i];
     if (typ(x) != t_COL) x = algtobasis(nf, x);
-    cx = denom(x); x = gmul(x,cx);
-    k = pvaluation(cx, p, &u);
-    if (!gcmp1(u)) /* could avoid the inversion, but prnZ is small--> cheap */
-      x = gmul(x, mpinvmod(u, prnZ));
-    if (k)
-      zpow = addii(zpow, mulsi(k, (GEN)e[i]));
-    (void)int_elt_val(nf, x, p, mul, &x, VERYBIGINT);
-    newg[i] = (long)colreducemodmat(x, prn, NULL);
+    x = Q_remove_denom(x, &cx);
+    if (cx)
+    {
+      long k = pvaluation(cx, p, &u);
+      if (!gcmp1(u)) /* could avoid the inversion, but prkZ is small--> cheap */
+        x = gmul(x, mpinvmod(u, prkZ));
+      if (k)
+        zpow = addii(zpow, mulsi(k, (GEN)e[i]));
+    }
+    (void)int_elt_val(nf, x, p, mul, &x);
+    newg[i] = (long)colreducemodHNF(x, prk, NULL);
   }
   if (zpow == gzero) setlg(newg, l);
   else
   {
-    newg[i] = (long)z;
+    newg[i] = (long)FpV_red(special_anti_uniformizer(nf, pr), prkZ);
     e = concatsp(e, negi(zpow));
   }
   e = gmod(e, EX);
-  return famat_to_nf_modideal_coprime(nf, newg, e, prn);
+  return famat_to_nf_modideal_coprime(nf, newg, e, prk);
 }
 
 GEN
 famat_ideallog(GEN nf, GEN g, GEN e, GEN bid)
 {
-  ulong av = avma;
+  gpmem_t av = avma;
   GEN vp = gmael(bid, 3,1), ep = gmael(bid, 3,2), arch = gmael(bid,1,2);
   GEN cyc = gmael(bid,2,2), list_set = (GEN)bid[4], U = (GEN)bid[5];
   GEN p1,y0,x,y, psigne;
-  long i;
+  long i, l;
   if (lg(cyc) == 1) return cgetg(1,t_COL);
   y0 = y = cgetg(lg(U), t_COL);
   psigne = zsigne(nf, to_famat(g,e), arch);
-  for (i=1; i<lg(vp); i++)
+  l = lg(vp);
+  for (i=1; i < l; i++)
   {
     GEN pr = (GEN)vp[i], prk;
-    prk = idealpow(nf, pr, (GEN)ep[i]);
+    prk = (l==2)? gmael(bid,1,1): idealpow(nf, pr, (GEN)ep[i]);
     /* TODO: FIX group exponent (should be mod prk, not f !) */
     x = famat_makecoprime(nf, g, e, pr, prk, (GEN)cyc[1]);
     y = zinternallog_pk(nf, x, y, pr, prk, (GEN)list_set[i], &psigne);
@@ -1407,6 +1506,29 @@ famat_to_nf_modidele(GEN nf, GEN g, GEN e, GEN bid)
 }
 
 GEN
+famat_to_arch(GEN nf, GEN fa, long prec)
+{
+  GEN g,e, y = NULL;
+  long i,l;
+
+  if (typ(fa) != t_MAT) return get_arch(nf, fa, prec);
+  if (lg(fa) == 1) return zerovec(lg(nf[6])-1);
+  g = (GEN)fa[1];
+  e = (GEN)fa[2]; l = lg(e);
+  for (i=1; i<l; i++)
+  {
+    GEN t, x = (GEN)g[i];
+    if (typ(x) != t_COL) x = algtobasis(nf,x);
+    x = Q_primpart(x);
+    /* multiplicative arch would be better (save logs), but exponents overflow
+     * [ could keep track of expo separately, but not worth it ] */
+    t = gmul(get_arch(nf,x,prec), (GEN)e[i]);
+    y = y? gadd(y,t): t;
+  }
+  return y;
+}
+
+GEN
 vecmul(GEN x, GEN y)
 {
   long i,lx, tx = typ(x);
@@ -1429,6 +1551,17 @@ vecinv(GEN x)
 }
 
 GEN
+vecpow(GEN x, GEN n)
+{
+  long i,lx, tx = typ(x);
+  GEN z;
+  if (is_scalar_t(tx)) return powgi(x,n);
+  lx = lg(x); z = cgetg(lx, tx);
+  for (i=1; i<lx; i++) z[i] = (long)vecpow((GEN)x[i], n);
+  return z;
+}
+
+GEN
 vecdiv(GEN x, GEN y) { return vecmul(x, vecinv(y)); }
 
 /* x,y assumed to be of the same type, either
@@ -1450,8 +1583,9 @@ GEN
 arch_inv(GEN x) {
   switch (typ(x)) {
     case t_POLMOD: return ginv(x);
+    case t_COL:    return vecinv(x);
     case t_MAT:    return famat_inv(x);
-    default:       return gneg(x); /* t_COL, t_VEC */
+    default:       return gneg(x); /* t_VEC */
   }
 }
 
@@ -1459,6 +1593,7 @@ GEN
 arch_pow(GEN x, GEN n) {
   switch (typ(x)) {
     case t_POLMOD: return powgi(x,n);
+    case t_COL:    return vecpow(x,n);
     case t_MAT:    return famat_pow(x,n);
     default:       return gmul(n,x);
   }
@@ -1468,7 +1603,8 @@ arch_pow(GEN x, GEN n) {
 GEN
 idealmul(GEN nf, GEN x, GEN y)
 {
-  long tx,ty,av,f;
+  gpmem_t av;
+  long tx,ty,f;
   GEN res,ax,ay,p1;
 
   tx = idealtyp(&x,&ax);
@@ -1489,11 +1625,15 @@ idealmul(GEN nf, GEN x, GEN y)
           p1 = idealhermite_aux(nf, element_mul(nf,x,y));
           break;
         case id_PRIME:
-          p1 = gmul((GEN)y[1],x);
-          p1 = two_to_hnf(nf,p1, element_mul(nf,(GEN)y[2],x));
+        {
+          GEN mx = eltmul_get_table(nf, x);
+          GEN mpi= eltmul_get_table(nf, (GEN)y[2]);
+          p1 = concatsp(gmul(mx,(GEN)y[1]), gmul(mx,mpi));
+          p1 = idealmat_to_hnf(nf, p1);
           break;
+        }
         default: /* id_MAT */
-          p1 = idealmat_mul(nf,y, principalideal_aux(nf,x));
+          p1 = idealmat_to_hnf(nf, gmul(eltmul_get_table(nf,x), y));
       }break;
 
     case id_PRIME:
@@ -1518,7 +1658,7 @@ idealmul(GEN nf, GEN x, GEN y)
 GEN
 idealnorm(GEN nf, GEN x)
 {
-  long av = avma,tetpil;
+  gpmem_t av = avma,tetpil;
   GEN y;
 
   nf = checknf(nf);
@@ -1536,7 +1676,6 @@ idealnorm(GEN nf, GEN x)
 }
 
 /* inverse */
-extern GEN gauss_triangle_i(GEN A, GEN B,GEN t);
 
 /* rewritten from original code by P.M & M.H.
  *
@@ -1547,9 +1686,9 @@ extern GEN gauss_triangle_i(GEN A, GEN B,GEN t);
 static GEN
 hnfideal_inv(GEN nf, GEN I)
 {
-  GEN J, dI = denom(I), IZ,dual;
+  GEN J, dI, IZ,dual;
 
-  if (gcmp1(dI)) dI = NULL; else I = gmul(I,dI);
+  I = Q_remove_denom(I, &dI);
   if (lg(I)==1) err(talker, "cannot invert zero ideal");
   IZ = gcoeff(I,1,1); /* I \cap Z */
   if (!signe(IZ)) err(talker, "cannot invert zero ideal");
@@ -1576,7 +1715,8 @@ GEN
 idealinv(GEN nf, GEN x)
 {
   GEN res,ax;
-  long av=avma, tx = idealtyp(&x,&ax);
+  gpmem_t av=avma;
+  long tx = idealtyp(&x,&ax);
 
   res = ax? cgetg(3,t_VEC): NULL;
   nf=checknf(nf); av=avma;
@@ -1595,7 +1735,7 @@ idealinv(GEN nf, GEN x)
           case t_COL: x = gmul((GEN)nf[7],x); break;
           case t_POLMOD: x = (GEN)x[2]; break;
         }
-        x = ginvmod(x,(GEN)nf[1]);
+        x = QX_invmod(x,(GEN)nf[1]);
       }
       x = idealhermite_aux(nf,x); break;
     case id_PRIME:
@@ -1617,18 +1757,31 @@ idealpowprime_spec(GEN nf, GEN vp, GEN n, GEN *d)
   if (s < 0) n = negi(n);
   /* now n > 0 */
   x = dummycopy(vp);
-  n1 = dvmdii(n, (GEN)x[3], &r);
-  if (signe(r)) n1 = addis(n1,1); /* n1 = ceil(n/e) */
-  x[1] = (long)powgi((GEN)x[1],n1);
-  if (s < 0)
+  if (is_pm1(n)) /* n = 1 special cased for efficiency */
   {
-    x[2] = ldiv(element_pow(nf,(GEN)x[5],n), powgi((GEN)vp[1],subii(n,n1)));
-    *d = (GEN)x[1];
+    if (s < 0)
+    {
+      x[2] = x[5];
+      *d = (GEN)x[1];
+    }
+    else
+      *d = NULL;
   }
   else
   {
-    x[2] = (long)element_pow(nf,(GEN)x[2],n);
-    *d = NULL;
+    n1 = dvmdii(n, (GEN)x[3], &r);
+    if (signe(r)) n1 = addis(n1,1); /* n1 = ceil(n/e) */
+    x[1] = (long)powgi((GEN)x[1],n1);
+    if (s < 0)
+    {
+      x[2] = ldiv(element_pow(nf,(GEN)x[5],n), powgi((GEN)vp[1],subii(n,n1)));
+      *d = (GEN)x[1];
+    }
+    else
+    {
+      x[2] = (long)element_pow(nf,(GEN)x[2],n);
+      *d = NULL;
+    }
   }
   return x;
 }
@@ -1647,30 +1800,45 @@ idealpowprime(GEN nf, GEN vp, GEN n)
   return x;
 }
 
-/* x * vp^n */
+/* x * vp^n. Assume x in HNF (possibly non-integral) */
 GEN
 idealmulpowprime(GEN nf, GEN x, GEN vp, GEN n)
 {
-  GEN denx,y,d;
+  GEN cx,y,dx;
 
   if (!signe(n)) return x;
   nf = checknf(nf);
-  y = idealpowprime_spec(nf, vp, n, &d);
-  denx = denom(x);
-  if (gcmp1(denx)) denx = d; else
+
+  /* inert, special cased for efficiency */
+  if (itos((GEN)vp[4]) == degpol(nf[1]))
+    return gmul(x, powgi((GEN)vp[1], n));
+
+  y = idealpowprime_spec(nf, vp, n, &dx);
+  x = Q_primitive_part(x, &cx);
+  if (cx && dx)
   {
-    x = gmul(denx,x);
-    if (d) denx = mulii(d,denx);
+    cx = gdiv(cx, dx);
+    if (typ(cx) != t_FRAC) dx = NULL;
+    else { dx = (GEN)cx[2]; cx = (GEN)cx[1]; }
+    if (is_pm1(cx)) cx = NULL;
   }
-  x = idealmulspec(nf,x, (GEN)y[1], (GEN)y[2]);
-  return denx? gdiv(x,denx): x;
+  x = idealmulspec(nf,x,y);
+  if (cx) x = gmul(x,cx);
+  if (dx) x = gdiv(x,dx);
+  return x;
 }
+GEN
+idealdivpowprime(GEN nf, GEN x, GEN vp, GEN n)
+{
+  return idealmulpowprime(nf,x,vp, negi(n));
+}
 
 /* raise the ideal x to the power n (in Z) */
 GEN
 idealpow(GEN nf, GEN x, GEN n)
 {
-  long tx,N,av,s,i;
+  gpmem_t av;
+  long tx,N,s;
   GEN res,ax,m,cx,n1,a,alpha;
 
   if (typ(n) != t_INT) err(talker,"non-integral exponent in idealpow");
@@ -1696,12 +1864,11 @@ idealpow(GEN nf, GEN x, GEN n)
       default:
         n1 = (s<0)? negi(n): n;
 
-        cx = content(x); if (gcmp1(cx)) cx = NULL; else x = gdiv(x,cx);
+        x = Q_primitive_part(x, &cx);
         a=ideal_two_elt(nf,x); alpha=(GEN)a[2]; a=(GEN)a[1];
-        m = cgetg(N+1,t_MAT); a = powgi(a,n1);
         alpha = element_pow(nf,alpha,n1);
-        for (i=1; i<=N; i++) m[i]=(long)element_mulid(nf,alpha,i);
-        x = hnfmodid(m, a);
+        m = eltmul_get_table(nf, alpha);
+        x = hnfmodid(m, powgi(a,n1));
         if (s<0) x = hnfideal_inv(nf,x);
         if (cx) x = gmul(x, powgi(cx,n));
     }
@@ -1717,59 +1884,65 @@ idealpow(GEN nf, GEN x, GEN n)
 GEN
 idealpows(GEN nf, GEN ideal, long e)
 {
-  long court[] = {evaltyp(t_INT) | m_evallg(3),0,0};
+  long court[] = {evaltyp(t_INT) | _evallg(3),0,0};
   affsi(e,court); return idealpow(nf,ideal,court);
 }
 
-GEN
-init_idele(GEN nf)
+static GEN
+_idealmulred(GEN nf, GEN x, GEN y, long prec)
 {
-  GEN x = cgetg(3,t_VEC);
-  long RU;
-  nf = checknf(nf); RU = lg(nf[6])-1;
-  x[2] = (long)zerovec(RU); return x;
+  return ideallllred(nf,idealmul(nf,x,y), NULL, prec);
 }
 
+typedef struct {
+  GEN nf;
+  long prec;
+} idealred_muldata;
+
+static GEN
+_mul(void *data, GEN x, GEN y)
+{
+  idealred_muldata *D = (idealred_muldata *)data;
+  return _idealmulred(D->nf,x,y,D->prec);
+}
+static GEN
+_sqr(void *data, GEN x)
+{
+  return _mul(data,x,x);
+}
+
 /* compute x^n (x ideal, n integer), reducing along the way */
 GEN
 idealpowred(GEN nf, GEN x, GEN n, long prec)
 {
-  long i,j,m,av=avma, s = signe(n);
-  GEN y, p1;
+  gpmem_t av=avma;
+  long s = signe(n);
+  idealred_muldata D;
+  GEN y;
 
   if (typ(n) != t_INT) err(talker,"non-integral exponent in idealpowred");
-  if (signe(n) == 0) return idealpow(nf,x,n);
-  p1 = n+2; m = *p1;
-  y = x; j=1+bfffo(m); m<<=j; j = BITS_IN_LONG-j;
-  for (i=lgefint(n)-2;;)
-  {
-    for (; j; m<<=1,j--)
-    {
-      y = idealmul(nf,y,y);
-      if (m < 0) y = idealmul(nf,y,x);
-      y = ideallllred(nf,y,NULL,prec);
-    }
-    if (--i == 0) break;
-    m = *++p1; j = BITS_IN_LONG;
-  }
+  if (s == 0) return idealpow(nf,x,n);
+  D.nf  = nf;
+  D.prec= prec;
+  y = leftright_pow(x, n, (void*)&D, &_sqr, &_mul);
+
   if (s < 0) y = idealinv(nf,y);
-  if (y == x) y = ideallllred(nf,x,NULL,prec);
+  if (s < 0 || is_pm1(n)) y = ideallllred(nf,y,NULL,prec);
   return gerepileupto(av,y);
 }
 
 GEN
 idealmulred(GEN nf, GEN x, GEN y, long prec)
 {
-  long av=avma;
-  x = idealmul(nf,x,y);
-  return gerepileupto(av, ideallllred(nf,x,NULL,prec));
+  gpmem_t av = avma;
+  return gerepileupto(av, _idealmulred(nf,x,y,prec));
 }
 
 long
 isideal(GEN nf,GEN x)
 {
-  long N,av,i,j,k,tx=typ(x),lx;
-  GEN p1,minv;
+  gpmem_t av;
+  long N,i,j,tx=typ(x),lx;
 
   nf=checknf(nf); lx=lg(x);
   if (tx==t_VEC && lx==3) { x=(GEN)x[1]; tx=typ(x); lx=lg(x); }
@@ -1779,26 +1952,25 @@ isideal(GEN nf,GEN x)
   if (typ(x)==t_VEC) return (lx==6);
   if (typ(x)!=t_MAT) return 0;
   if (lx == 1) return 1;
-  N=lgef(nf[1])-2; if (lg(x[1]) != N) return 0;
+  N = degpol(nf[1]); if (lg(x[1])-1 != N) return 0;
 
-  av=avma;
-  if (lx != N) x = idealmat_to_hnf(nf,x);
-  x = gdiv(x,content(x)); minv=ginv(x);
+  av = avma;
+  if (lx-1 != N) x = idealmat_to_hnf(nf,x);
+  x = primpart(x);
 
-  for (i=1; i<N; i++)
-    for (j=1; j<N; j++)
-    {
-      p1=gmul(minv, element_mulid(nf,(GEN)x[i],j));
-      for (k=1; k<N; k++)
-	if (typ(p1[k])!=t_INT) { avma=av; return 0; }
-    }
+  for (i=1; i<=N; i++)
+    for (j=2; j<=N; j++)
+      if (! hnf_invimage(x, element_mulid(nf,(GEN)x[i],j)))
+      {
+        avma = av; return 0;
+      }
   avma=av; return 1;
 }
 
 GEN
 idealdiv(GEN nf, GEN x, GEN y)
 {
-  long av=avma,tetpil;
+  gpmem_t av=avma,tetpil;
   GEN z=idealinv(nf,y);
 
   tetpil=avma; return gerepile(av,tetpil,idealmul(nf,x,z));
@@ -1826,7 +1998,7 @@ idealdiv(GEN nf, GEN x, GEN y)
 GEN
 idealdivexact(GEN nf, GEN x0, GEN y0)
 {
-  ulong av = avma;
+  gpmem_t av = avma;
   GEN x,y,Nx,Ny,Nz, cy = content(y0);
 
   nf = checknf(nf);
@@ -1860,7 +2032,8 @@ idealdivexact(GEN nf, GEN x0, GEN y0)
 GEN
 idealintersect(GEN nf, GEN x, GEN y)
 {
-  long av=avma,lz,i,N;
+  gpmem_t av=avma;
+  long lz,i,N;
   GEN z,dx,dy;
 
   nf=checknf(nf); N=degpol(nf[1]);
@@ -1877,55 +2050,75 @@ idealintersect(GEN nf, GEN x, GEN y)
   return gerepileupto(av,z);
 }
 
-static GEN
-computet2twist(GEN nf, GEN vdir)
+/*******************************************************************/
+/*                                                                 */
+/*                      T2-IDEAL REDUCTION                         */
+/*                                                                 */
+/*******************************************************************/
+
+GEN
+computeGtwist(GEN nf, GEN vdir)
 {
-  long j, ru = lg(nf[6]);
-  GEN p1,MC, mat = (GEN)nf[5];
+  long i, j, k, l, lG, v, r1, r2;
+  GEN p1, G = gmael(nf,5,2);
 
-  if (!vdir) return (GEN)mat[3];
-  MC=(GEN)mat[2]; p1=cgetg(ru,t_MAT);
-  for (j=1; j<ru; j++)
+  if (!vdir) return G;
+  l = lg(vdir); lG = lg(G);
+  p1 = dummycopy(G);
+  nf_get_sign(nf, &r1, &r2);
+  for (i=1; i<l; i++)
   {
-    GEN v = (GEN)vdir[j];
-    if (gcmp0(v))
-      p1[j] = MC[j];
-    else if (typ(v) == t_INT)
-      p1[j] = lmul2n((GEN)MC[j],itos(v)<<1);
-    else
-      p1[j] = lmul((GEN)MC[j],gpui(stoi(4),v,0));
+    v = vdir[i]; if (!v) continue;
+    if (i <= r1) {
+      for (j=1; j<lG; j++) coeff(p1,i,j) = lmul2n(gcoeff(p1,i,j), v);
+    } else {
+      k = (i<<1) - r1;
+      for (j=1; j<lG; j++)
+      {
+        coeff(p1,k-1,j) = lmul2n(gcoeff(p1,k-1,j), v);
+        coeff(p1,k  ,j) = lmul2n(gcoeff(p1,k  ,j), v);
+      }
+    }
   }
-  return mulmat_real(p1,(GEN)mat[1]);
+  return p1;
 }
 
 static GEN
 chk_vdir(GEN nf, GEN vdir)
 {
+  long l, i, t;
+  GEN v;
   if (!vdir || gcmp0(vdir)) return NULL;
-  if (typ(vdir)!=t_VEC || lg(vdir) != lg(nf[6])) err(idealer5);
-  return vdir;
+  l = lg(vdir);
+  if (l != lg(nf[6])) err(talker, "incorrect vector length in idealred");
+  t = typ(vdir);
+  if (t == t_VECSMALL) return vdir;
+  if (t != t_VEC) err(talker, "not a vector in idealred");
+  v = cgetg(l, t_VECSMALL);
+  for (i=1; i<l; i++) v[i] = itos(gceil((GEN)vdir[i]));
+  return v;
 }
 
 /* assume I in NxN matrix form (not necessarily HNF) */
 static GEN
-ideallllred_elt_i(GEN *ptnf, GEN I, GEN vdir, long *ptprec)
+idealred_elt_i(GEN *ptnf, GEN I, GEN vdir, long *ptprec)
 {
-  GEN T2, u, y, nf = *ptnf;
+  GEN G, u, y, nf = *ptnf;
   long i, e, prec = *ptprec;
 
+  e = gexpo(I)>>TWOPOTBITS_IN_LONG;
+  if (e < 0) e = 0;
   for (i=1; ; i++)
   {
-    T2 = computet2twist(nf,vdir);
-    y = qf_base_change(T2,I,1);
-    e = 1 + (gexpo(y)>>TWOPOTBITS_IN_LONG);
-    if (e < 0) e = 0;
-    u = lllgramintern(y,100,1, e + prec);
+    G = computeGtwist(nf,vdir);
+    y = gmul(G, I);
+    u = lllintern(y, 100, 1, prec);
     if (u) break;
 
     if (i == MAXITERPOL) err(accurer,"ideallllred");
     prec = (prec<<1)-2;
     if (DEBUGLEVEL) err(warnprec,"ideallllred",prec);
-    nf = nfnewprec(nf, (e>>1)+prec);
+    nf = nfnewprec(nf, e + prec);
   }
   *ptprec = prec;
   *ptnf = nf;
@@ -1933,16 +2126,21 @@ ideallllred_elt_i(GEN *ptnf, GEN I, GEN vdir, long *pt
 }
 
 GEN
-ideallllred_elt(GEN nf, GEN I)
+ideallllred_elt(GEN nf, GEN I, GEN vdir)
 {
   long prec = DEFAULTPREC;
-  return ideallllred_elt_i(&nf, I, NULL, &prec);
+  return idealred_elt_i(&nf, I, vdir, &prec);
 }
+GEN
+idealred_elt(GEN nf, GEN I)
+{
+  return ideallllred_elt(nf, I, NULL);
+}
 
 GEN
 ideallllred(GEN nf, GEN I, GEN vdir, long prec)
 {
-  ulong av = avma;
+  gpmem_t av = avma;
   long N,i,nfprec;
   GEN J,I0,Ired,res,aI,y,x,Nx,b,c1,c,pol;
 
@@ -1959,12 +2157,12 @@ ideallllred(GEN nf, GEN I, GEN vdir, long prec)
   if (DEBUGLEVEL>5) msgtimer("entering idealllred");
   I0 = I;
   if (typ(I) != id_MAT || lg(I) != N+1) I = idealhermite_aux(nf,I);
-  c1 = content(I); if (gcmp1(c1)) c1 = NULL; else I = gdiv(I,c1);
+  I = primitive_part(I, &c1);
   if (2 * expi(gcoeff(I,1,1)) >= bit_accuracy(nfprec))
-    Ired = gmul(I, lllintpartial(I));
+    Ired = lllint_ip(I,4);
   else
     Ired = I;
-  y = ideallllred_elt_i(&nf, Ired, chk_vdir(nf,vdir), &prec);
+  y = idealred_elt_i(&nf, Ired, chk_vdir(nf,vdir), &prec);
 
   if (isnfscalar(y))
   { /* already reduced */
@@ -1974,16 +2172,16 @@ ideallllred(GEN nf, GEN I, GEN vdir, long prec)
   if (DEBUGLEVEL>5) msgtimer("LLL reduction");
 
   x = gmul((GEN)nf[7], y); Nx = subres(pol,x);
-  b = gmul(Nx, ginvmod(x,pol));
-  b = algtobasis_intern(nf,b);
+  b = gmul(Nx, QX_invmod(x,pol));
+  b = algtobasis_i(nf,b);
   J = cgetg(N+1,t_MAT); /* = I Nx / x integral */
   for (i=1; i<=N; i++)
     J[i] = (long)element_muli(nf,b,(GEN)Ired[i]);
-  c = content(J); if (!gcmp1(c)) J = gdiv(J,c);
+  J = Q_primitive_part(J, &c);
  /* c = content (I Nx / x) = Nx / den(I/x) --> d = den(I/x) = Nx / c
   * J = (d I / x); I[1,1] = I \cap Z --> d I[1,1] belongs to J and Z */
   if (isnfscalar((GEN)I[1]))
-    b = mulii(gcoeff(I,1,1), divii(Nx, c));
+    b = mulii(gcoeff(I,1,1), c? diviiexact(Nx, c): Nx);
   else
     b = detint(J);
   I = hnfmodid(J,b);
@@ -1994,15 +2192,24 @@ END:
 
   switch(typ(aI))
   {
-    case t_POLMOD: case t_MAT: /* compute y, I0 = J y */
+    case t_POLMOD: /* compute y, I0 = J y */
       if (!Nx) y = c1;
       else
       {
-        if (c1) c = gmul(c,c1);
-        y = gmul(x, gdiv(c,Nx));
+        c = mul_content(c,c1);
+        y = c? gmul(x, gdiv(c,Nx)): gdiv(x, Nx);
       }
       break;
 
+    case t_MAT: /* compute y, I0 = J y */
+      if (!Nx) y = c1;
+      else
+      {
+        c = mul_content(c,c1);
+        y = c? gmul(y, gdiv(c,Nx)): gdiv(y, Nx);
+      }
+      break;
+
     case t_COL:
       if (y) y = vecinv(gmul(gmael(nf,5,1), y));
       break;
@@ -2020,8 +2227,9 @@ END:
 GEN
 minideal(GEN nf, GEN x, GEN vdir, long prec)
 {
-  long av = avma, N, tx;
-  GEN p1,y;
+  gpmem_t av = avma;
+  long N, tx;
+  GEN y;
 
   nf = checknf(nf);
   vdir = chk_vdir(nf,vdir);
@@ -2030,311 +2238,381 @@ minideal(GEN nf, GEN x, GEN vdir, long prec)
   if (tx == id_PRINCIPAL) return gcopy(x);
   if (tx != id_MAT || lg(x) != N+1) x = idealhermite_aux(nf,x);
 
-  p1 = computet2twist(nf,vdir);
-  y = qf_base_change(p1,x,0);
-  y = gmul(x, (GEN)lllgram(y,prec)[1]);
+  y = gmul(computeGtwist(nf,vdir), x);
+  y = gmul(x, (GEN)lll(y,prec)[1]);
   return gerepileupto(av, principalidele(nf,y,prec));
 }
+
+/*******************************************************************/
+/*                                                                 */
+/*                   APPROXIMATION THEOREM                         */
+/*                                                                 */
+/*******************************************************************/
+
+/* write x = x1 x2, x2 maximal s.t. (x2,f) = 1, return x2 */
 static GEN
-appr_reduce(GEN s, GEN y, long N)
+coprime_part(GEN x, GEN f)
 {
-  GEN p1,u,z = cgetg(N+2,t_MAT);
-  long i;
-
-  s=gmod(s,gcoeff(y,1,1)); y=gmul(y,lllint(y));
-  for (i=1; i<=N; i++) z[i]=y[i]; z[N+1]=(long)s;
-  u=(GEN)ker(z)[1]; p1 = denom(u);
-  if (!gcmp1(p1)) u=gmul(u,p1);
-  p1=(GEN)u[N+1]; setlg(u,N+1);
-  for (i=1; i<=N; i++) u[i]=lround(gdiv((GEN)u[i],p1));
-  return gadd(s, gmul(y,u));
+  for (;;)
+  {
+    f = mppgcd(x, f); if (is_pm1(f)) break;
+    x = diviiexact(x, f);
+  }
+  return x;
 }
 
-/* Given a fractional ideal x (if fl=0) or a prime ideal factorization
- * with possibly zero or negative exponents (if fl=1), gives a b such that
- * v_p(b)=v_p(x) for all prime ideals p dividing x (or in the factorization)
- * and v_p(b)>=0 for all other p, using the (standard) proof given in GTM 138.
- * Certainly not the most efficient, but sure.
- */
-GEN
-idealappr0(GEN nf, GEN x, long fl)
+/* x t_INT, f ideal. Write x = x1 x2, sqf(x1) | f, (x2,f) = 1. Return x2 */
+static GEN
+nf_coprime_part(GEN nf, GEN x, GEN *listpr)
 {
-  long av=avma,tetpil,i,j,k,l,N,r,r2;
-  GEN fact,fact2,list,ep,ep1,ep2,y,z,v,p1,p2,p3,p4,s,pr,alpha,beta,den;
+  long v, j, lp = lg(listpr), N = degpol(nf[1]);
+  GEN x1, x2, ex, p, pr;
 
-  if (DEBUGLEVEL>4)
+#if 0 /*1) via many gcds. Expensive ! */
+  GEN f = idealprodprime(nf, (GEN)listpr);
+  f = hnfmodid(f, x); /* first gcd is less expensive since x in Z */
+  x = gscalmat(x, N);
+  for (;;)
   {
-    fprintferr(" entree dans idealappr0() :\n");
-    fprintferr(" x = "); outerr(x);
+    if (gcmp1(gcoeff(f,1,1))) break;
+    x = idealdivexact(nf, x, f);
+    f = hnfmodid(concatsp(f,x), gcoeff(x,1,1)); /* gcd(f,x) */
   }
-  nf=checknf(nf); N=degpol(nf[1]);
-  if (fl)
+  x2 = x;
+#else /*2) from prime decomposition */
+  x1 = NULL;
+  for (j=1; j<lp; j++)
   {
-    if (typ(x)!=t_MAT || lg(x)!=3)
-      err(talker,"not a prime ideal factorization in idealappr0");
-    fact=x; list=(GEN)fact[1]; ep=(GEN)fact[2]; r=lg(list);
-    if (r==1) return gscalcol_i(gun,N);
-    for (i=1; i<r; i++)
-      if (signe(ep[i]) < 0) break;
-    if (i < r)
-    {
-      ep1=cgetg(r,t_COL);
-      for (i=1; i<r; i++)
-        ep1[i] = (signe(ep[i])>=0)? zero: lnegi((GEN)ep[i]);
-      fact[2]=(long)ep1; beta=idealappr0(nf,fact,1);
-      fact2=idealfactor(nf,beta);
-      p1=(GEN)fact2[1]; r2=lg(p1);
-      ep2=(GEN)fact2[2]; l=r+r2-1;
-      z=cgetg(l,t_VEC); for (i=1; i<r; i++) z[i]=list[i];
-      ep1=cgetg(l,t_VEC);
-      for (i=1; i<r; i++)
-        ep1[i] = (signe(ep[i])<=0)? zero: licopy((GEN)ep[i]);
-      j=r-1;
-      for (i=1; i<r2; i++)
-      {
-        p3=(GEN)p1[i]; k=1;
-        while (k<r &&
-          (    !gegal((GEN)p3[1],gmael(list,k,1))
-            || !element_val(nf,(GEN)p3[2],(GEN)list[k]) )) k++;
-        if (k==r) { j++; z[j]=(long)p3; ep1[j]=ep2[i]; }
-      }
-      fact=cgetg(3,t_MAT);
-      fact[1]=(long)z; setlg(z,j+1);
-      fact[2]=(long)ep1; setlg(ep1,j+1);
-      alpha=idealappr0(nf,fact,1); tetpil=avma;
-      if (DEBUGLEVEL>2)
-      {
-        fprintferr(" alpha = "); outerr(alpha);
-        fprintferr(" beta = "); outerr(beta);
-      }
-      return gerepile(av,tetpil,element_div(nf,alpha,beta));
-    }	
-    y=idmat(N);
-    for (i=1; i<r; i++)
-    {
-      pr=(GEN)list[i];
-      if (signe(ep[i]))
-      {
-        p4=addsi(1,(GEN)ep[i]); p1=powgi((GEN)pr[1],p4);
-	if (cmpis((GEN)pr[4],N))
-	{
-	  p2=cgetg(3,t_MAT);
-          p2[1]=(long)gscalcol_i(p1, N);
-	  p2[2]=(long)element_pow(nf,(GEN)pr[2],p4);
-          y=idealmat_mul(nf,y,p2);
-	}
-	else y=gmul(p1,y);
-      }
-      else y=idealmulprime(nf,y,pr);
-    }
+    pr = listpr[j]; p = (GEN)pr[1];
+    v = ggval(x, p); if (!v) continue;
+
+    ex = mulsi(v, (GEN)pr[3]); /* = v_pr(x) > 0 */
+    x1 = x1? idealmulpowprime(nf, x1, pr, ex)
+           : idealpow(nf, pr, ex);
   }
+  x = gscalmat(x, N);
+  x2 = x1? idealdivexact(nf, x, x1): x;
+#endif
+  return x2;
+}
+
+/* assume L is a list of prime ideals. Return the product */
+GEN
+idealprodprime(GEN nf, GEN L)
+{
+  long l = lg(L), i;
+  GEN z;
+
+  if (l == 1) return idmat(degpol(nf[1]));
+  z = prime_to_ideal(nf, (GEN)L[1]);
+  for (i=2; i<l; i++) z = idealmulprime(nf,z, (GEN)L[i]);
+  return z;
+}
+
+/* assume L is a list of prime ideals. Return prod L[i]^e[i] */
+GEN
+factorbackprime(GEN nf, GEN L, GEN e)
+{
+  long l = lg(L), i;
+  GEN z;
+
+  if (l == 1) return idmat(degpol(nf[1]));
+  z = idealpow(nf, (GEN)L[1], (GEN)e[1]);
+  for (i=2; i<l; i++)
+    if (signe(e[i])) z = idealmulpowprime(nf,z, (GEN)L[i],(GEN)e[i]);
+  return z;
+}
+
+/* compute anti-uniformizer for pr, coprime to f outside of pr, integral
+ * outside of p below pr.
+ * sqf = product or primes dividing f, NULL if f a prime power*/
+GEN
+anti_unif_mod_f(GEN nf, GEN pr, GEN sqf)
+{
+  GEN U, V, UV, uv, cx, t, p = (GEN)pr[1];
+  if (!sqf) return gdiv((GEN)pr[5], p);
   else
   {
-    den=denom(x); if (gcmp1(den)) den=NULL; else x=gmul(den,x);
-    x=idealhermite_aux(nf,x);
-    fact=idealfactor(nf,x);
-    list=(GEN)fact[1]; ep=(GEN)fact[2]; r=lg(list);
-    if (r==1) { avma=av; return gscalcol_i(gun,N); }
-    if (den)
-    {
-      fact2=idealfactor(nf,den);
-      p1=(GEN)fact2[1]; r2=lg(p1);
-      l=r+r2-1;
-      z=cgetg(l,t_COL);   for (i=1; i<r; i++) z[i]=list[i];
-      ep1=cgetg(l,t_COL); for (i=1; i<r; i++) ep1[i]=ep[i];
-      j=r-1;
-      for (i=1; i<r2; i++)
-      {
-	p3=(GEN)p1[i]; k=1;
-	while (k<r && !gegal((GEN)list[k],p3)) k++;
-	if (k==r){ j++; z[j]=(long)p3; ep1[j]=zero; }
-      }
-      fact=cgetg(3,t_MAT);
-      fact[1]=(long)z; setlg(z,j+1);
-      fact[2]=(long)ep1; setlg(ep1,j+1);
-      alpha=idealappr0(nf,fact,1);
-      if (DEBUGLEVEL>2) { fprintferr(" alpha = "); outerr(alpha); }
-      tetpil=avma; return gerepile(av,tetpil,gdiv(alpha,den));
-    }
-    y=x; for (i=1; i<r; i++) y=idealmulprime(nf,y,(GEN)list[i]);
+    GEN d,d1,d2, sqfZ = gcoeff(sqf,1,1); /* = (sqf \cap Z) = (V \cap Z) * p */
+    U = idealpow(nf,pr,gdeux);
+    V = idealdivpowprime(nf,sqf,pr,gun);
+    uv = addone_nored(U, V);
+    UV = idealmul(nf,sqf,pr);
+    t = (GEN)pr[2];
+    t = makeprimetoideal(nf, UV, uv, t); /* cf unif_mod_f */
+
+    t = element_inv(nf, t);
+    t = primitive_part(t, &cx); /* p | denom(cx) */
+    d = denom(cx);
+    d2 = coprime_part(d, sqfZ);
+    d1 = diviiexact(d, d2); /* p | d1 */
+    cx = gmod(gmul(cx,d1), mulii(d1, gcoeff(V,1,1)));
+    t = gdiv(colreducemodHNF(gmul(cx,t), gmul(d1,V), NULL), d1);
+    return t; /* v_pr[i](t) = -1, v_pr[j](t) = 0 if i != j */
   }
+}
 
-  z=cgetg(r,t_VEC);
+/* compute integral uniformizer for pr, coprime to f outside of pr.
+ * sqf = product or primes dividing f, NULL if f a prime power*/
+GEN
+unif_mod_f(GEN nf, GEN pr, GEN sqf)
+{
+  GEN U, V, UV, uv, t = (GEN)pr[2];
+  if (!sqf) return t;
+  else
+  {
+    U = idealpow(nf,pr,gdeux);
+    V = idealdivpowprime(nf,sqf,pr,gun);
+    uv = addone_nored(U, V);
+    UV = idealmulprime(nf,sqf,pr);
+    return makeprimetoideal(nf, UV, uv, t);
+  }
+}
+
+static GEN
+appr_reduce(GEN s, GEN y)
+{
+  long i, N = lg(y)-1;
+  GEN d, u, z = cgetg(N+2,t_MAT);
+
+  s = gmod(s, gcoeff(y,1,1));
+  y = lllint_ip(y,100);
+  for (i=1; i<=N; i++) z[i] = y[i];
+  z[N+1] = (long)s;
+  u = (GEN)ker(z)[1];
+  u = Q_remove_denom(u, NULL);
+  d = (GEN)u[N+1]; setlg(u,N+1);
+  for (i=1; i<=N; i++) u[i] = (long)diviiround((GEN)u[i],d);
+  return gadd(s, gmul(y,u));
+}
+
+/* L0 in K^*. 
+ * If ptd1 == NULL, assume (L0,f) = 1
+ *   return L integral, L0 = L mod f 
+ *
+ * Otherwise, assume v_pr(L0) <= 0 for all pr | f and set *ptd1 = d1
+ *   return L integral, L0 = L/d1 mod f, and such that 
+ *   if (L*I,f) = 1 for some integral I, then d1 | L*I  */
+GEN
+make_integral(GEN nf, GEN L0, GEN f, GEN *listpr, GEN *ptd1)
+{
+  GEN fZ, t, L, D2, d1, d2, d;
+
+  if (ptd1) *ptd1 = NULL;
+  L = Q_remove_denom(L0, &d);
+  if (!d) return L0;
+  
+  /* L0 = L / d, L integral */
+  fZ = gcoeff(f,1,1);
+  /* Kill denom part coprime to fZ */
+  d2 = coprime_part(d, fZ);
+  t = mpinvmod(d2, fZ); if (!is_pm1(t)) L = gmul(L,t);
+  if (egalii(d, d2)) return L;
+
+  d1 = diviiexact(d, d2);
+  /* L0 = (L / d1) mod f. d1 not coprime to f
+   * write (d1) = D1 D2, D2 minimal, (D2,f) = 1. */
+  D2 = nf_coprime_part(nf, d1, listpr);
+  t = idealaddtoone_i(nf, D2, f); /* in D2, 1 mod f */
+  L = element_muli(nf,t,L);
+
+  /* if (L0, f) = 1, then L in D1 ==> in D1 D2 = (d1) */
+  if (!ptd1) return Q_div_to_int(L, d1); /* exact division */
+
+  *ptd1 = d1; return L;
+}
+
+void
+check_listpr(GEN x)
+{
+  long l = lg(x), i;
+  for (i=1; i<l; i++) checkprimeid((GEN)x[i]);
+}
+
+/* Given a prime ideal factorization with possibly zero or negative
+ * exponents, gives b such that v_p(b) = v_p(x) for all prime ideals pr | x
+ * and v_pr(b)> = 0 for all other pr.
+ * For optimal performance, all [anti-]uniformizers should be precomputed,
+ * but no support for this yet.
+ * No garbage collecting */
+GEN
+idealapprfact_i(GEN nf, GEN x)
+{
+  gpmem_t av = avma;
+  GEN tau, pi, z, d, sqf, sqfsafe, list, e, e2;
+  long s, flag, i, r, N;
+
+  nf = checknf(nf);
+  if (typ(x) != t_MAT || lg(x) != 3)
+    err(talker,"not a factorization in idealapprfact");
+  check_listpr((GEN)x[1]);
+  if (DEBUGLEVEL>4) {
+    fprintferr(" entering idealapprfact():\n");
+    fprintferr(" x = %Z\n", x);
+  }
+  list= (GEN)x[1];
+  e   = (GEN)x[2]; r = lg(e); N = degpol(nf[1]);
+  if (r==1) return gscalcol_i(gun, N);
+
+  sqf = (r == 2)? NULL: idealprodprime(nf, list);
+  sqfsafe = NULL;
+  z = gun; flag = 0;
   for (i=1; i<r; i++)
   {
-    pr=(GEN)list[i]; p4=addsi(1, (GEN)ep[i]); p1=powgi((GEN)pr[1],p4);
-    if (cmpis((GEN)pr[4],N))
+    s = signe(e[i]);
+    if (s < 0)
     {
-      p2=cgetg(3,t_MAT);
-      p2[1]=(long)gscalcol_i(p1,N);
-      p2[2]=(long)element_pow(nf,(GEN)pr[5],p4);
-      z[i]=ldiv(idealmat_mul(nf,y,p2),p1);
+      flag = 1;
+      if (!sqfsafe) sqfsafe = sqf? sqf: prime_to_ideal(nf, (GEN)list[1]);
+      tau = anti_unif_mod_f(nf, (GEN)list[i], sqf);
+      tau = make_integral(nf, tau, sqfsafe, (GEN*)list, &d);
+      tau = gdiv(tau, d);
+      z = element_mul(nf, z, element_pow(nf, tau, negi((GEN)e[i])));
     }
-    else z[i]=ldiv(y,p1);
+    else if (s > 0)
+    {
+      pi = unif_mod_f(nf, (GEN)list[i], sqf);
+      z = element_mul(nf, z, element_pow(nf, pi, (GEN)e[i]));
+    }
   }
-  v=idealaddmultoone(nf,z);
-  s=cgetg(N+1,t_COL); for (i=1; i<=N; i++) s[i]=zero;
-  for (i=1; i<r; i++)
+  if (z == gun) { avma = av; return gscalcol_i(gun, N); }
+  e2 = cgetg(r, t_VEC);
+  for (i=1; i<r; i++) e2[i] = laddis((GEN)e[i], 1);
+  x = factorbackprime(nf, list,e2);
+  if (flag) /* denominator */
   {
-    pr=(GEN)list[i];
-    if (signe(ep[i]))
-      s=gadd(s,element_mul(nf,(GEN)v[i],element_pow(nf,(GEN)pr[2],(GEN)ep[i])));
-    else s=gadd(s,(GEN)v[i]);
+    z = Q_remove_denom(z, &d);
+    x = Q_muli_to_int(x, d);
   }
-  p3 = appr_reduce(s,y,N);
-  if (DEBUGLEVEL>2)
-    { fprintferr(" sortie de idealappr0 p3 = "); outerr(p3); }
-  return gerepileupto(av,p3);
+  else d = NULL;
+  z = appr_reduce(z, x);
+  if (d) z = gdiv(z, d);
+  return z;
 }
 
+GEN
+idealapprfact(GEN nf, GEN x) {
+  gpmem_t av = avma;
+  return gerepileupto(av, idealapprfact_i(nf, x));
+}
+
+GEN
+idealappr(GEN nf, GEN x) {
+  gpmem_t av = avma;
+  return gerepileupto(av, idealapprfact_i(nf, idealfactor(nf, x)));
+}
+
+GEN
+idealappr0(GEN nf, GEN x, long fl) {
+  return fl? idealapprfact(nf, x): idealappr(nf, x);
+}
+
 /* Given a prime ideal factorization x with possibly zero or negative exponents,
- * and a vector y of elements of nf, gives a b such that
- * v_p(b-y_p)>=v_p(x) for all prime ideals p in the ideal factorization
+ * and a vector w of elements of nf, gives a b such that
+ * v_p(b-w_p)>=v_p(x) for all prime ideals p in the ideal factorization
  * and v_p(b)>=0 for all other p, using the (standard) proof given in GTM 138.
- * Certainly not the most efficient, but sure.
- */
+ * Certainly not the most efficient, but sure. */
 GEN
-idealchinese(GEN nf, GEN x, GEN y)
+idealchinese(GEN nf, GEN x, GEN w)
 {
-  long ty=typ(y),av=avma,i,j,k,l,N,r,r2;
-  GEN fact,fact2,list,ep,ep1,ep2,z,t,v,p1,p2,p3,p4,s,pr,den;
+  gpmem_t av = avma;
+  long ty = typ(w), i,N,r;
+  GEN L,e,z,t,y,v,s,den;
 
-  if (DEBUGLEVEL>4)
-  {
-    fprintferr(" entree dans idealchinese() :\n");
-    fprintferr(" x = "); outerr(x);
-    fprintferr(" y = "); outerr(y);
-  }
   nf=checknf(nf); N=degpol(nf[1]);
-  if (typ(x)!=t_MAT ||(lg(x)!=3))
+  if (typ(x) != t_MAT || lg(x) != 3)
     err(talker,"not a prime ideal factorization in idealchinese");
-  fact=x; list=(GEN)fact[1]; ep=(GEN)fact[2]; r=lg(list);
-  if (!is_vec_t(ty) || lg(y)!=r)
+  L = (GEN)x[1]; r = lg(L);
+  e = (GEN)x[2];
+  if (!is_vec_t(ty) || lg(w) != r)
     err(talker,"not a suitable vector of elements in idealchinese");
   if (r==1) return gscalcol_i(gun,N);
 
-  den=denom(y);
+  den = denom(w);
   if (!gcmp1(den))
   {
-    fact2=idealfactor(nf,den);
-    p1=(GEN)fact2[1]; r2=lg(p1);
-    ep2=(GEN)fact2[2]; l=r+r2-1;
-    z=cgetg(l,t_VEC); for (i=1; i<r; i++) z[i]=list[i];
-    ep1=cgetg(l,t_VEC); for (i=1; i<r; i++) ep1[i]=ep[i];
-    j=r-1;
-    for (i=1; i<r2; i++)
-    {
-      p3=(GEN)p1[i]; k=1;
-      while (k<r && !gegal((GEN)list[k],p3)) k++;
-      if (k==r) { j++; z[j]=(long)p3; ep1[j]=ep2[i]; }
-      else ep1[k]=ladd((GEN)ep1[k],(GEN)ep2[i]);
-    }
-    r=j+1; setlg(z,r); setlg(ep1,r); list=z; ep=ep1;
+    GEN fa = famat_reduce(famat_mul(x, idealfactor(nf,den)));
+    L = (GEN)fa[1]; r = lg(L);
+    e = (GEN)fa[2];
   }
   for (i=1; i<r; i++)
-    if (signe(ep[i])<0) ep[i] = zero;
-  t=idmat(N);
+    if (signe(e[i]) < 0) e[i] = zero;
+
+  y = factorbackprime(nf, L, e);
+  z = cgetg(r,t_VEC);
   for (i=1; i<r; i++)
+    z[i] = (long)idealdivpowprime(nf,y, (GEN)L[i], (GEN)e[i]);
+  v = idealaddmultoone(nf,z); s = NULL;
+  for (i=1; i<r; i++)
   {
-    pr=(GEN)list[i]; p4=(GEN)ep[i];
-    if (signe(p4))
-    {
-      if (cmpis((GEN)pr[4],N))
-      {
-	p2=cgetg(3,t_MAT);
-        p2[1]=(long)gscalcol_i(powgi((GEN)pr[1],p4), N);
-	p2[2]=(long)element_pow(nf,(GEN)pr[2],p4);
-        t=idealmat_mul(nf,t,p2);
-      }
-      else t=gmul(powgi((GEN)pr[1],p4),t);
-    }
+    t = element_mul(nf, (GEN)v[i], (GEN)w[i]);
+    s = s? gadd(s,t): t;
   }
-  z=cgetg(r,t_VEC);
+  return gerepileupto(av, appr_reduce(s,y));
+}
+
+static GEN
+mat_ideal_two_elt2(GEN nf, GEN x, GEN a)
+{
+  GEN L, e, fact = idealfactor(nf,a);
+  long i, v, r;
+  L = (GEN)fact[1]; r = lg(L);
+  e = (GEN)fact[2];
   for (i=1; i<r; i++)
   {
-    pr=(GEN)list[i]; p4=(GEN)ep[i];
-    if (cmpis((GEN)pr[4],N))
-    {
-      p2=cgetg(3,t_MAT); p1=powgi((GEN)pr[1],p4);
-      p2[1]=(long)gscalcol_i(p1,N);
-      p2[2]=(long)element_pow(nf,(GEN)pr[5],p4);
-      z[i]=ldiv(idealmat_mul(nf,t,p2),p1);
-    }
-    else z[i]=ldiv(t,powgi((GEN)pr[1],p4));
+    v = idealval(nf,x,(GEN)L[i]);
+    e[i] = lstoi(v);
   }
-  v=idealaddmultoone(nf,z);
-  s=cgetg(N+1,t_COL); for (i=1; i<=N; i++) s[i]=zero;
-  for (i=1; i<r; i++)
-    s = gadd(s,element_mul(nf,(GEN)v[i],(GEN)y[i]));
-
-  p3 = appr_reduce(s,t,N);
-  if (DEBUGLEVEL>2)
-    { fprintferr(" sortie de idealchinese() : p3 = "); outerr(p3); }
-  return gerepileupto(av,p3);
+  return centermod(idealapprfact_i(nf,fact), gcoeff(x,1,1));
 }
 
-GEN
-idealappr(GEN nf, GEN x) { return idealappr0(nf,x,0); }
-
-GEN
-idealapprfact(GEN nf, GEN x) { return idealappr0(nf,x,1); }
-
 /* Given an integral ideal x and a in x, gives a b such that
- * x=aZ_K+bZ_K using a different algorithm than ideal_two_elt
- */
+ * x = aZ_K + bZ_K using the approximation theorem */
 GEN
 ideal_two_elt2(GEN nf, GEN x, GEN a)
 {
-  long ta=typ(a), av=avma,tetpil,i,r;
-  GEN con,ep,b,list,fact;
+  gpmem_t av = avma;
+  GEN cx, b;
 
   nf = checknf(nf);
-  if (!is_extscalar_t(ta) && typ(a)!=t_COL)
-    err(typeer,"ideal_two_elt2");
+  if (typ(a) != t_COL) a = algtobasis(nf, a);
   x = idealhermite_aux(nf,x);
   if (gcmp0(x))
   {
     if (!gcmp0(a)) err(talker,"element not in ideal in ideal_two_elt2");
     avma=av; return gcopy(a);
   }
-  con = content(x);
-  if (gcmp1(con)) con = NULL; else { x = gdiv(x,con); a = gdiv(a,con); }
-  a = principalideal(nf,a);
-  if (!gcmp1(denom(gauss(x,a))))
+  x = Q_primitive_part(x, &cx);
+  if (cx) a = gdiv(a, cx);
+  if (!hnf_invimage(x, a))
     err(talker,"element does not belong to ideal in ideal_two_elt2");
 
-  fact=idealfactor(nf,a); list=(GEN)fact[1];
-  r=lg(list); ep = (GEN)fact[2];
-  for (i=1; i<r; i++) ep[i]=lstoi(idealval(nf,x,(GEN)list[i]));
-  b = centermod(idealappr0(nf,fact,1), gcoeff(x,1,1));
-  tetpil=avma; b = con? gmul(b,con): gcopy(b);
-  return gerepile(av,tetpil,b);
+  b = mat_ideal_two_elt2(nf, x, a);
+  b = cx? gmul(b,cx): gcopy(b);
+  return gerepileupto(av, b);
 }
 
-/* Given 2 integral ideals x and y in a number field nf gives a beta
- * belonging to nf such that beta.x is an integral ideal coprime to y
- */
+/* Given 2 integral ideals x and y in nf, returns a beta in nf such that
+ * beta * x is an integral ideal coprime to y */
 GEN
 idealcoprime(GEN nf, GEN x, GEN y)
 {
-  long av=avma,tetpil,i,r;
-  GEN fact,list,p2,ep;
+  gpmem_t av = avma;
+  long i, r;
+  GEN list, ep, fact = idealfactor(nf,y);
 
-  if (DEBUGLEVEL>4)
-  {
-    fprintferr(" entree dans idealcoprime() :\n");
-    fprintferr(" x = "); outerr(x);
-    fprintferr(" y = "); outerr(y);
-  }
-  fact=idealfactor(nf,y); list=(GEN)fact[1];
-  r=lg(list); ep = (GEN)fact[2];
-  for (i=1; i<r; i++) ep[i]=lstoi(-idealval(nf,x,(GEN)list[i]));
-  tetpil=avma; p2=idealappr0(nf,fact,1);
-  if (DEBUGLEVEL>4)
-    { fprintferr(" sortie de idealcoprime() : p2 = "); outerr(p2); }
-  return gerepile(av,tetpil,p2);
+  list = (GEN)fact[1];
+  ep   = (GEN)fact[2]; r = lg(ep);
+  for (i=1; i<r; i++) ep[i] = lstoi(-idealval(nf,x,(GEN)list[i]));
+  return gerepileupto(av, idealapprfact_i(nf,fact));
 }
 
+/*******************************************************************/
+/*                                                                 */
+/*                  LINEAR ALGEBRA OVER Z_K  (HNF,SNF)             */
+/*                                                                 */
+/*******************************************************************/
+
 /* returns the first index i<=n such that x=v[i] if it exits, 0 otherwise */
 long
 isinvector(GEN v, GEN x, long n)
@@ -2347,44 +2625,19 @@ isinvector(GEN v, GEN x, long n)
 }
 
 GEN
-elt_mul_get_table(GEN nf, GEN x)
-{
-  long i,lx = lg(x);
-  GEN mul=cgetg(lx,t_MAT);
-
-  /* assume w_1 = 1 */
-  mul[1]=(long)x;
-  for (i=2; i<lx; i++) mul[i] = (long)element_mulid(nf,x,i);
-  return mul;
-}
-
-GEN
-elt_mul_table(GEN mul, GEN z)
-{
-  long av = avma, i, lx = lg(mul);
-  GEN p1 = gmul((GEN)z[1], (GEN)mul[1]);
-
-  for (i=2; i<lx; i++)
-    if (!gcmp0((GEN)z[i])) p1 = gadd(p1, gmul((GEN)z[i], (GEN)mul[i]));
-  return gerepileupto(av, p1);
-}
-
-GEN
 element_mulvec(GEN nf, GEN x, GEN v)
 {
-  long lv=lg(v),i;
+  long lv = lg(v), i;
   GEN y = cgetg(lv,t_COL);
 
   if (typ(x) == t_COL)
   {
-    GEN mul = elt_mul_get_table(nf,x);
-    for (i=1; i<lv; i++)
-      y[i] = (long)elt_mul_table(mul,(GEN)v[i]);
+    GEN mul = eltmul_get_table(nf,x);
+    for (i=1; i<lv; i++) y[i] = lmul(mul,(GEN)v[i]);
   }
   else
   { /* scalar */
-    for (i=1; i<lv; i++)
-      y[i] = lmul(x, (GEN)v[i]);
+    for (i=1; i<lv; i++) y[i] = lmul(x, (GEN)v[i]);
   }
   return y;
 }
@@ -2393,34 +2646,31 @@ static GEN
 element_mulvecrow(GEN nf, GEN x, GEN m, long i, long lim)
 {
   long lv,j;
-  GEN y, mul = elt_mul_get_table(nf,x);
+  GEN y, mul = eltmul_get_table(nf,x);
 
   lv=min(lg(m),lim+1); y=cgetg(lv,t_VEC);
-  for (j=1; j<lv; j++)
-    y[j] = (long)elt_mul_table(mul,gcoeff(m,i,j));
+  for (j=1; j<lv; j++) y[j] = lmul(mul,gcoeff(m,i,j));
   return y;
 }
 
 /* Given an element x and an ideal in matrix form (not necessarily HNF),
- * gives an r such that x-r is in ideal and r is small. No checks
- */
+ * gives an r such that x-r is in ideal and r is small. No checks */
 GEN
 element_reduce(GEN nf, GEN x, GEN ideal)
 {
-  long tx=typ(x),av=avma,tetpil,N,i;
-  GEN p1,u;
+  long tx = typ(x), N, i;
+  gpmem_t av=avma;
+  GEN u,d;
 
   if (is_extscalar_t(tx))
-    x = algtobasis_intern(checknf(nf), x);
+    x = algtobasis_i(checknf(nf), x);
   N = lg(x);
   if (typ(ideal) != t_MAT || lg(ideal) != N) err(typeer,"element_reduce");
-  p1=cgetg(N+1,t_MAT);
-  for (i=1; i<N; i++) p1[i]=ideal[i];
-  p1[N]=(long)x; u=(GEN)ker(p1)[1];
-  p1=(GEN)u[N]; setlg(u,N);
-  for (i=1; i<N; i++) u[i]=lround(gdiv((GEN)u[i],p1));
-  u=gmul(ideal,u); tetpil=avma;
-  return gerepile(av,tetpil,gadd(u,x));
+  u = ker( concatsp(ideal,x) ); /* do NOT use deplin. Much, much slower -- KB */
+  u = (GEN)u[1];
+  d = (GEN)u[N]; setlg(u,N);
+  for (i=1; i<N; i++) u[i] = (long)gdivround((GEN)u[i],d);
+  return gerepileupto(av, gadd(x, gmul(ideal,u)));
 }
 
 /* A torsion-free module M over Z_K will be given by a row vector [A,I] with
@@ -2430,98 +2680,153 @@ element_reduce(GEN nf, GEN x, GEN ideal)
  * that [A,I] is a pseudo-basis if k=n
  */
 
-/* Given a torsion-free module x as above outputs a pseudo-basis for x in
- * Hermite Normal Form
- */
+#define swap(x,y) { long _t=x; x=y; y=_t; }
+
+/* Given a torsion-free module x as above outputs a pseudo-basis for x in HNF */
 GEN
 nfhermite(GEN nf, GEN x)
 {
-  long av0 = avma, av,lim,i,j,def,k,m;
+  long i, j, def, k, m;
+  gpmem_t av0 = avma, av, lim;
   GEN p1,p2,y,A,I,J;
 
-  nf=checknf(nf);
+  nf = checknf(nf);
   if (typ(x)!=t_VEC || lg(x)!=3) err(talker,"not a module in nfhermite");
-  A=(GEN)x[1]; I=(GEN)x[2]; k=lg(A)-1;
+  A = (GEN)x[1];
+  I = (GEN)x[2]; k = lg(A)-1;
   if (typ(A)!=t_MAT) err(talker,"not a matrix in nfhermite");
   if (typ(I)!=t_VEC || lg(I)!=k+1)
     err(talker,"not a correct ideal list in nfhermite");
   if (!k) err(talker,"not a matrix of maximal rank in nfhermite");
-  m=lg(A[1])-1;
-  if (k<m) err(talker,"not a matrix of maximal rank in nfhermite");
+  m = lg(A[1])-1;
+  if (k < m) err(talker,"not a matrix of maximal rank in nfhermite");
 
   av = avma; lim = stack_lim(av, 1);
-  p1 = cgetg(k+1,t_MAT); for (j=1; j<=k; j++) p1[j]=A[j];
-  A = p1; I = dummycopy(I);
-  J = cgetg(k+1,t_VEC);
-  for (j=1; j<=k; j++)
-  {
-    if (typ(I[j])!=t_MAT) I[j]=(long)idealhermite_aux(nf,(GEN)I[j]);
-    J[j] = zero;
-  }
-
-  def = k+1;
+  A = matalgtobasis(nf,A);
+  I = dummycopy(I);
+  J = zerovec(k); def = k+1;
   for (i=m; i>=1; i--)
   {
     GEN den,p4,p5,p6,u,v,newid, invnewid = NULL;
 
     def--; j=def; while (j>=1 && gcmp0(gcoeff(A,i,j))) j--;
     if (!j) err(talker,"not a matrix of maximal rank in nfhermite");
-    if (j==def) j--;
-    else
-    {
-      p1=(GEN)A[j]; A[j]=A[def]; A[def]=(long)p1;
-      p1=(GEN)I[j]; I[j]=I[def]; I[def]=(long)p1;
-    }
-    p1=gcoeff(A,i,def); p2=element_inv(nf,p1);
-    A[def]=(long)element_mulvec(nf,p2,(GEN)A[def]);
-    I[def]=(long)idealmul(nf,p1,(GEN)I[def]);
+    if (j==def) j--; else { swap(A[j], A[def]); swap(I[j], I[def]); }
+    p1 = gcoeff(A,i,def);
+    p2 = element_inv(nf, p1);
+    A[def] = (long)element_mulvec(nf,p2,(GEN)A[def]);
+    I[def] = (long)idealmul(nf,p1,(GEN)I[def]);
     for (  ; j; j--)
     {
-      p1=gcoeff(A,i,j);
+      p1 = gcoeff(A,i,j);
       if (gcmp0(p1)) continue;
 
-      p2=idealmul(nf,p1,(GEN)I[j]);
+      p2 = idealmul(nf,p1,(GEN)I[j]);
       newid = idealadd(nf,p2,(GEN)I[def]);
+
       invnewid = hnfideal_inv(nf,newid);
       p4 = idealmul(nf, p2,        invnewid);
       p5 = idealmul(nf,(GEN)I[def],invnewid);
       y = idealaddtoone(nf,p4,p5);
-      u=element_div(nf,(GEN)y[1],p1); v=(GEN)y[2];
-      p6=gsub((GEN)A[j],element_mulvec(nf,p1,(GEN)A[def]));
-      A[def]=ladd(element_mulvec(nf,u,(GEN)A[j]),
-                  element_mulvec(nf,v,(GEN)A[def]));
-      A[j]=(long)p6;
-      I[j]=(long)idealmul(nf,idealmul(nf,(GEN)I[j],(GEN)I[def]),invnewid);
-      I[def]=(long)newid; den=denom((GEN)I[j]);
+      u = (GEN)y[1]; u = element_div(nf, u, p1);
+      v = (GEN)y[2];
+      p6 = gsub((GEN)A[j], element_mulvec(nf, p1,(GEN)A[def]));
+      A[def] = ladd(element_mulvec(nf, u, (GEN)A[j]),
+                    element_mulvec(nf, v, (GEN)A[def]));
+      A[j] = (long)p6;
+      I[j] = (long)idealmul(nf,idealmul(nf,(GEN)I[j],(GEN)I[def]),invnewid);
+      I[def] = (long)newid; den = denom((GEN)I[j]);
       if (!gcmp1(den))
       {
-        I[j]=lmul(den,(GEN)I[j]);
-        A[j]=ldiv((GEN)A[j],den);
+        I[j] = lmul(den,(GEN)I[j]);
+        A[j] = ldiv((GEN)A[j],den);
       }
     }
     if (!invnewid) invnewid = hnfideal_inv(nf,(GEN)I[def]);
-    p1=(GEN)I[def]; J[def]=(long)invnewid;
+    J[def] = (long)invnewid;
+    p1 = (GEN)I[def];
     for (j=def+1; j<=k; j++)
     {
-      p2 = gsub(element_reduce(nf,gcoeff(A,i,j),idealmul(nf,p1,(GEN)J[j])),
-                gcoeff(A,i,j));
-      A[j] = ladd((GEN)A[j],element_mulvec(nf,p2,(GEN)A[def]));
+      GEN c = gcoeff(A,i,j);
+      p2 = gsub(element_reduce(nf,c, idealmul(nf,p1,(GEN)J[j])), c);
+      A[j] = ladd((GEN)A[j], element_mulvec(nf, p2,(GEN)A[def]));
     }
     if (low_stack(lim, stack_lim(av1,1)))
     {
-      GEN *gptr[3];
+      GEN *gptr[3]; gptr[0]=&A; gptr[1]=&I; gptr[2]=&J;
       if(DEBUGMEM>1) err(warnmem,"nfhermite, i = %ld", i);
-      gptr[0]=&A; gptr[1]=&I; gptr[2]=&J; gerepilemany(av,gptr,3);
+      gerepilemany(av,gptr,3);
     }
   }
   y = cgetg(3,t_VEC);
-  p1 = cgetg(m+1,t_MAT); y[1] = (long)p1;
-  p2 = cgetg(m+1,t_VEC); y[2] = (long)p2;
-  for (j=1; j<=m; j++) p1[j] = lcopy((GEN)A[j+k-m]);
-  for (j=1; j<=m; j++) p2[j] = lcopy((GEN)I[j+k-m]);
-  return gerepileupto(av0, y);
+  A += k-m; A[0] = evaltyp(t_MAT)|evallg(m+1); y[1] = (long)A;
+  I += k-m; I[0] = evaltyp(t_VEC)|evallg(m+1); y[2] = (long)I;
+  return gerepilecopy(av0, y);
 }
 
+static GEN
+idealmulelt(GEN nf, GEN elt, GEN x)
+{
+  long t = typ(elt);
+  GEN z;
+  if (t == t_POL || t == t_POLMOD) elt = algtobasis(nf,elt);
+  if (isnfscalar(elt)) elt = (GEN)elt[1];
+  z = element_mulvec(nf, elt, x);
+  settyp(z, t_MAT); return z;
+}
+
+static GEN
+zero_nfbezout(GEN nf,GEN b, GEN A,GEN B,GEN *u,GEN *v,GEN *w,GEN *di)
+{
+  gpmem_t av, tetpil;
+  GEN pab,d;
+
+  d=idealmulelt(nf,b,B); *di=idealinv(nf,idealmat_to_hnf(nf,d));
+  av=avma; pab=idealmul(nf,A,B); tetpil=avma;
+  *w=gerepile(av,tetpil, idealmul(nf,pab,*di));
+  *v=element_inv(nf,b);
+  *u=gzero; return d;
+}
+
+/* Given elements a,b and ideals A, B, outputs d = a.A+b.B and gives
+ * di=d^-1, w=A.B.di, u, v such that au+bv=1 and u in A.di, v in
+ * B.di. Assume A, B non-zero, but a or b can be zero (not both)
+ */
+static GEN
+nfbezout(GEN nf,GEN a,GEN b, GEN A,GEN B, GEN *u,GEN *v,GEN *w,GEN *di)
+{
+  GEN pab,pu,pv,pw,uv,d,dinv,pa,pb,pa1,pb1, *gptr[5];
+  gpmem_t av, tetpil;
+
+  if (gcmp0(a))
+  {
+    if (gcmp0(b)) err(talker,"both elements zero in nfbezout");
+    return zero_nfbezout(nf,b,A,B,u,v,w,di);
+  }
+  if (gcmp0(b))
+    return zero_nfbezout(nf,a,B,A,v,u,w,di);
+
+  av = avma;
+  pa = idealmulelt(nf,a,A);
+  pb = idealmulelt(nf,b,B);
+
+  d=idealadd(nf,pa,pb); dinv=idealinv(nf,d);
+  pa1=idealmullll(nf,pa,dinv);
+  pb1=idealmullll(nf,pb,dinv);
+  uv=idealaddtoone(nf,pa1,pb1);
+  pab=idealmul(nf,A,B); tetpil=avma;
+
+  pu=element_div(nf,(GEN)uv[1],a);
+  pv=element_div(nf,(GEN)uv[2],b);
+  d=gcopy(d); dinv=gcopy(dinv);
+  pw=idealmul(nf,pab,dinv);
+
+  *u=pu; *v=pv; *w=pw; *di=dinv;
+  gptr[0]=u; gptr[1]=v; gptr[2]=w; gptr[3]=di;
+  gptr[4]=&d; gerepilemanysp(av,tetpil,gptr,5);
+  return d;
+}
+
 /* A torsion module M over Z_K will be given by a row vector [A,I,J] with
  * three components. I=[b_1,...,b_n] is a row vector of k fractional ideals
  * given in HNF, J=[a_1,...,a_n] is a row vector of n fractional ideals in
@@ -2536,7 +2841,8 @@ nfhermite(GEN nf, GEN x)
 GEN
 nfsmith(GEN nf, GEN x)
 {
-  long av,tetpil,i,j,k,l,lim,c,n,m,N;
+  long i, j, k, l, c, n, m, N;
+  gpmem_t av, tetpil, lim;
   GEN p1,p2,p3,p4,z,b,u,v,w,d,dinv,unnf,A,I,J;
 
   nf=checknf(nf); N=degpol(nf[1]);
@@ -2562,82 +2868,75 @@ nfsmith(GEN nf, GEN x)
   {
     do
     {
-      c=0;
+      c = 0;
       for (j=i-1; j>=1; j--)
       {
-	p1=gcoeff(A,i,j);
-	if (!gcmp0(p1))
-	{
-	  p2=gcoeff(A,i,i);
-	  d=nfbezout(nf,p2,p1,(GEN)J[i],(GEN)J[j],&u,&v,&w,&dinv);
-	  if (!gcmp0(u))
-	  {
-	    if (!gcmp0(v))
-	      b=gadd(element_mulvec(nf,u,(GEN)A[i]),
-	             element_mulvec(nf,v,(GEN)A[j]));
-	    else b=element_mulvec(nf,u,(GEN)A[i]);
-	  }
-	  else b=element_mulvec(nf,v,(GEN)A[j]);
-	  A[j]=lsub(element_mulvec(nf,p2,(GEN)A[j]),
-	            element_mulvec(nf,p1,(GEN)A[i]));
-	  A[i]=(long)b; J[j]=(long)w; J[i]=(long)d;
-	}
+	p1=gcoeff(A,i,j); if (gcmp0(p1)) continue;
+
+        p2 = gcoeff(A,i,i);
+        d = nfbezout(nf,p2,p1,(GEN)J[i],(GEN)J[j],&u,&v,&w,&dinv);
+        if (gcmp0(u)) b = element_mulvec(nf,v,(GEN)A[j]);
+        else
+        {
+          if (gcmp0(v)) b = element_mulvec(nf,u,(GEN)A[i]);
+          else
+            b = gadd(element_mulvec(nf,u,(GEN)A[i]),
+                     element_mulvec(nf,v,(GEN)A[j]));
+        }
+        A[j] = lsub(element_mulvec(nf,p2,(GEN)A[j]),
+                    element_mulvec(nf,p1,(GEN)A[i]));
+        A[i] = (long)b; J[j] = (long)w; J[i] = (long)d;
       }
       for (j=i-1; j>=1; j--)
       {
-	p1=gcoeff(A,j,i);
-	if (!gcmp0(p1))
-	{
-	  p2=gcoeff(A,i,i);
-	  d=nfbezout(nf,p2,p1,(GEN)I[i],(GEN)I[j],&u,&v,&w,&dinv);
-	  if (gcmp0(u))
-	    b=element_mulvecrow(nf,v,A,j,i);
-	  else
-	  {
-	    if (gcmp0(v))
-	      b=element_mulvecrow(nf,u,A,i,i);
-	    else
-	      b=gadd(element_mulvecrow(nf,u,A,i,i),
-	             element_mulvecrow(nf,v,A,j,i));
-	  }
-	  p3=gsub(element_mulvecrow(nf,p2,A,j,i),
-	          element_mulvecrow(nf,p1,A,i,i));
-	  for (k=1; k<=i; k++) { coeff(A,j,k)=p3[k]; coeff(A,i,k)=b[k]; }
-	  I[j]=(long)w; I[i]=(long)d; c++;
-	}
+	p1 = gcoeff(A,j,i); if (gcmp0(p1)) continue;
+
+        p2 = gcoeff(A,i,i);
+        d = nfbezout(nf,p2,p1,(GEN)I[i],(GEN)I[j],&u,&v,&w,&dinv);
+        if (gcmp0(u)) b = element_mulvecrow(nf,v,A,j,i);
+        else
+        {
+          if (gcmp0(v))
+            b = element_mulvecrow(nf,u,A,i,i);
+          else
+            b = gadd(element_mulvecrow(nf,u,A,i,i),
+                     element_mulvecrow(nf,v,A,j,i));
+        }
+        p3 = gsub(element_mulvecrow(nf,p2,A,j,i),
+                  element_mulvecrow(nf,p1,A,i,i));
+        for (k=1; k<=i; k++) { coeff(A,j,k) = p3[k]; coeff(A,i,k) = b[k]; }
+        I[j] = (long)w; I[i] = (long)d; c = 1;
       }
-      if (!c)
-      {
-	b=gcoeff(A,i,i); if (gcmp0(b)) break;
+      if (c) continue;
 
-	b=idealmul(nf,b,idealmul(nf,(GEN)J[i],(GEN)I[i]));
-	for (k=1; k<i; k++)
-	  for (l=1; l<i; l++)
-	  {
-	    p3 = gcoeff(A,k,l);
-	    if (!gcmp0(p3))
-            {
-              p4 = idealmul(nf,p3,idealmul(nf,(GEN)J[l],(GEN)I[k]));
-	      if (!gegal(idealadd(nf,b,p4), b))
-              {
-                b=idealdiv(nf,(GEN)I[k],(GEN)I[i]);
-                p4=gauss(idealdiv(nf,(GEN)J[i],idealmul(nf,p3,(GEN)J[l])),b);
-                l=1; while (l<=N && gcmp1(denom((GEN)p4[l]))) l++;
-                if (l>N) err(talker,"bug2 in nfsmith");
-                p3=element_mulvecrow(nf,(GEN)b[l],A,k,i);
-                for (l=1; l<=i; l++)
-                  coeff(A,i,l) = ladd(gcoeff(A,i,l),(GEN)p3[l]);
+      b=gcoeff(A,i,i); if (gcmp0(b)) break;
 
-                k = l = i; c = 1;
-              }
-            }
-	  }
-      }
+      b=idealmul(nf,b,idealmul(nf,(GEN)J[i],(GEN)I[i]));
+      for (k=1; k<i; k++)
+        for (l=1; l<i; l++)
+        {
+          p3 = gcoeff(A,k,l);
+          if (gcmp0(p3)) continue;
+
+          p4 = idealmul(nf,p3,idealmul(nf,(GEN)J[l],(GEN)I[k]));
+          if (hnfdivide(b, p4)) continue;
+
+          b = idealdiv(nf,(GEN)I[k],(GEN)I[i]);
+          p1 = idealdiv(nf,(GEN)J[i], idealmul(nf,p3,(GEN)J[l]));
+          p4 = gauss(p4, b);
+          l=1; while (l<=N && gcmp1(denom((GEN)p4[l]))) l++;
+          if (l>N) err(talker,"bug2 in nfsmith");
+          p3 = element_mulvecrow(nf,(GEN)b[l],A,k,i);
+          for (l=1; l<=i; l++)
+            coeff(A,i,l) = ladd(gcoeff(A,i,l),(GEN)p3[l]);
+
+          k = l = i; c = 1;
+        }
       if (low_stack(lim, stack_lim(av,1)))
       {
-        GEN *gptr[3];
+        GEN *gptr[3]; gptr[0]=&A; gptr[1]=&I; gptr[2]=&J;
 	if(DEBUGMEM>1) err(warnmem,"nfsmith");
-        gptr[0]=&A; gptr[1]=&I; gptr[2]=&J; gerepilemany(av,gptr,3);
+        gerepilemany(av,gptr,3);
       }
     }
     while (c);
@@ -2652,260 +2951,108 @@ nfsmith(GEN nf, GEN x)
   return gerepile(av,tetpil,z);
 }
 
-/*******************************************************************/
-/*                                                                 */
-/*          ALGEBRE LINEAIRE DANS LES CORPS DE NOMBRES             */
-/*                                                                 */
-/*******************************************************************/
-
 #define trivlift(x) ((typ(x)==t_POLMOD)? (GEN)x[2]: lift_intern(x))
 
 GEN
-element_mulmodpr2(GEN nf, GEN x, GEN y, GEN prhall)
+element_mulmodpr(GEN nf, GEN x, GEN y, GEN modpr)
 {
-  long av=avma;
+  gpmem_t av=avma;
   GEN p1;
 
-  nf=checknf(nf); checkprhall(prhall);
+  nf=checknf(nf); checkmodpr(modpr);
   p1 = element_mul(nf,x,y);
-  return gerepileupto(av,nfreducemodpr(nf,p1,prhall));
+  return gerepileupto(av,nfreducemodpr(nf,p1,modpr));
 }
 
-/* On ne peut PAS definir ca comme les autres par
- * #define element_divmodpr() nfreducemodpr(element_div())
- * car le element_div ne marche pas en general
- */
 GEN
-element_divmodpr(GEN nf, GEN x, GEN y, GEN prhall)
+element_divmodpr(GEN nf, GEN x, GEN y, GEN modpr)
 {
-  long av=avma;
+  gpmem_t av=avma;
   GEN p1;
 
-  nf=checknf(nf); checkprhall(prhall);
+  nf=checknf(nf); checkmodpr(modpr);
   p1=lift_intern(gdiv(gmodulcp(gmul((GEN)nf[7],trivlift(x)), (GEN)nf[1]),
                       gmodulcp(gmul((GEN)nf[7],trivlift(y)), (GEN)nf[1])));
-  p1=algtobasis_intern(nf,p1);
-  return gerepileupto(av,nfreducemodpr(nf,p1,prhall));
+  p1=algtobasis_i(nf,p1);
+  return gerepileupto(av,nfreducemodpr(nf,p1,modpr));
 }
 
 GEN
-element_invmodpr(GEN nf, GEN y, GEN prhall)
+element_invmodpr(GEN nf, GEN y, GEN modpr)
 {
-  long av=avma;
+  gpmem_t av=avma;
   GEN p1;
 
-  p1=ginvmod(gmul((GEN)nf[7],trivlift(y)), (GEN)nf[1]);
-  p1=algtobasis_intern(nf,p1);
-  return gerepileupto(av,nfreducemodpr(nf,p1,prhall));
+  p1=QX_invmod(gmul((GEN)nf[7],trivlift(y)), (GEN)nf[1]);
+  p1=algtobasis_i(nf,p1);
+  return gerepileupto(av,nfreducemodpr(nf,p1,modpr));
 }
 
 GEN
-element_powmodpr(GEN nf,GEN x,GEN k,GEN prhall)
+element_powmodpr(GEN nf,GEN x,GEN k,GEN pr)
 {
-  long av=avma,N,s;
-  GEN y,z;
+  gpmem_t av=avma;
+  GEN z,T,p,modpr;
 
-  nf=checknf(nf); checkprhall(prhall);
-  N=degpol(nf[1]);
-  s=signe(k); k=(s>=0)?k:negi(k);
-  z=x; y = gscalcol_i(gun,N);
-  for(;;)
-  {
-    if (mpodd(k)) y=element_mulmodpr(nf,z,y,prhall);
-    k=shifti(k,-1);
-    if (signe(k)) z=element_sqrmodpr(nf,z,prhall);
-    else
-    {
-      cgiv(k); if (s<0) y = element_invmodpr(nf,y,prhall);
-      return gerepileupto(av,y);
-    }
-  }
+  nf = checknf(nf);
+  modpr = nf_to_ff_init(nf,&pr,&T,&p);
+  z = nf_to_ff(nf,x,modpr);
+  z = FpXQ_pow(z,k,T,p);
+  z = ff_to_nf(z,modpr);
+  if (typ(z) != t_COL) z = algtobasis(nf, z);
+  return gerepilecopy(av, z);
 }
 
-/* x est une matrice dont les coefficients sont des vecteurs dans la base
- * d'entiers modulo un ideal premier prhall, sous forme reduite modulo prhall.
- */
 GEN
-nfkermodpr(GEN nf, GEN x, GEN prhall)
+nfkermodpr(GEN nf, GEN x, GEN pr)
 {
-  ulong av0,av,av1,lim;
-  long i,j,k,r,t,n,m,N;
-  GEN c,d,y,unnf,munnf,zeromodp,zeronf,p,pp,prh;
+  gpmem_t av = avma;
+  GEN T,p,modpr;
 
-  nf=checknf(nf); checkprhall(prhall);
+  nf = checknf(nf);
   if (typ(x)!=t_MAT) err(typeer,"nfkermodpr");
-  n=lg(x)-1; if (!n) return cgetg(1,t_MAT);
-  prh=(GEN)prhall[1]; av0=avma;
-  N=degpol(nf[1]); pp=gcoeff(prh,1,1);
-
-  zeromodp=gmodulsg(0,pp);
-  unnf=cgetg(N+1,t_COL); unnf[1]=(long)gmodulsg(1,pp);
-  zeronf=cgetg(N+1,t_COL); zeronf[1]=(long)zeromodp;
-
-  av=avma; munnf=cgetg(N+1,t_COL); munnf[1]=(long)gmodulsg(-1,pp);
-  for (i=2; i<=N; i++)
-    zeronf[i] = munnf[i] = unnf[i] = (long)zeromodp;
-
-  m=lg(x[1])-1; x=dummycopy(x); r=0;
-  c=new_chunk(m+1); for (k=1; k<=m; k++) c[k]=0;
-  d=new_chunk(n+1); av1=avma; lim=stack_lim(av1,1);
-  for (k=1; k<=n; k++)
-  {
-    j=1;
-    while (j<=m && (c[j] || gcmp0(gcoeff(x,j,k)))) j++;
-    if (j > m) { r++; d[k]=0; continue; }
-
-      p=element_divmodpr(nf,munnf,gcoeff(x,j,k),prhall);
-      c[j]=k; d[k]=j; coeff(x,j,k)=(long)munnf;
-      for (i=k+1; i<=n; i++)
-	coeff(x,j,i)=(long)element_mulmodpr(nf,p,gcoeff(x,j,i),prhall);
-      for (t=1; t<=m; t++)
-	if (t!=j)
-	{
-        p = gcoeff(x,t,k); if (gcmp0(p)) continue;
-        coeff(x,t,k) = (long)zeronf;
-	  for (i=k+1; i<=n; i++)
-            coeff(x,t,i)=ladd(gcoeff(x,t,i),
-	                      element_mulmodpr(nf,p,gcoeff(x,j,i),prhall));
-          if (low_stack(lim, stack_lim(av1,1)))
-          {
-            if (DEBUGMEM>1) err(warnmem,"nfkermodpr, k = %ld / %ld",k,n);
-            x=gerepilecopy(av1,x);
-          }
-	}
-  }
-  if (!r) { avma=av0; return cgetg(1,t_MAT); }
-  av1=avma; y=cgetg(r+1,t_MAT);
-  for (j=k=1; j<=r; j++,k++)
-  {
-    p=cgetg(n+1,t_COL); y[j]=(long)p; while (d[k]) k++;
-    for (i=1; i<k; i++) p[i]=d[i]? lcopy(gcoeff(x,d[i],k)): (long)zeronf;
-    p[k]=(long)unnf; for (i=k+1; i<=n; i++) p[i]=(long)zeronf;
-  }
-  return gerepile(av,av1,y);
+  modpr = nf_to_ff_init(nf, &pr,&T,&p);
+  x = modprM(lift(x), nf, modpr);
+  return gerepilecopy(av, modprM_lift(FqM_ker(x,T,p), modpr));
 }
 
 /* a.x=b ou b est un vecteur */
 GEN
-nfsolvemodpr(GEN nf, GEN a, GEN b, GEN prhall)
+nfsolvemodpr(GEN nf, GEN a, GEN b, GEN pr)
 {
-  ulong av = avma;
-  long nbli,nbco,i,j,k;
-  GEN aa,x,p,m,u;
+  gpmem_t av = avma;
+  GEN T,p,modpr;
 
-  nf=checknf(nf); checkprhall(prhall);
+  nf = checknf(nf);
   if (typ(a)!=t_MAT) err(typeer,"nfsolvemodpr");
-  nbco=lg(a)-1; nbli=lg(a[1])-1;
-  if (typ(b)!=t_COL) err(typeer,"nfsolvemodpr");
-  if (lg(b)!=nbco+1) err(mattype1,"nfsolvemodpr");
-  x=cgetg(nbli+1,t_COL);
-  for (j=1; j<=nbco; j++) x[j]=b[j];
-  aa=cgetg(nbco+1,t_MAT);
-  for (j=1; j<=nbco; j++)
-  {
-    aa[j]=lgetg(nbli+1,t_COL);
-    for (i=1; i<=nbli; i++) coeff(aa,i,j)=coeff(a,i,j);
-  }
-  for (i=1; i<nbli; i++)
-  {
-    p=gcoeff(aa,i,i); k=i;
-    if (gcmp0(p))
-    {
-      k=i+1; while (k<=nbli && gcmp0(gcoeff(aa,k,i))) k++;
-      if (k>nbco) err(matinv1);
-      for (j=i; j<=nbco; j++)
-      {
-	u=gcoeff(aa,i,j); coeff(aa,i,j)=coeff(aa,k,j);
-	coeff(aa,k,j)=(long)u;
-      }
-      u=(GEN)x[i]; x[i]=x[k]; x[k]=(long)u;
-      p=gcoeff(aa,i,i);
-    }
-    for (k=i+1; k<=nbli; k++)
-    {
-      m=gcoeff(aa,k,i);
-      if (!gcmp0(m))
-      {
-	m=element_divmodpr(nf,m,p,prhall);
-	for (j=i+1; j<=nbco; j++)
-	  coeff(aa,k,j)=lsub(gcoeff(aa,k,j),
-	                     element_mulmodpr(nf,m,gcoeff(aa,i,j),prhall));
-	x[k]=lsub((GEN)x[k],element_mulmodpr(nf,m,(GEN)x[i],prhall));
-      }
-    }
-  }
-  /* Resolution systeme triangularise */
-  p=gcoeff(aa,nbli,nbco); if (gcmp0(p)) err(matinv1);
-
-  x[nbli]=(long)element_divmodpr(nf,(GEN)x[nbli],p,prhall);
-  for (i=nbli-1; i>0; i--)
-  {
-    m=(GEN)x[i];
-    for (j=i+1; j<=nbco; j++)
-      m=gsub(m,element_mulmodpr(nf,gcoeff(aa,i,j),(GEN)x[j],prhall));
-    x[i]=(long)element_divmodpr(nf,m,gcoeff(aa,i,i),prhall);
-  }
-  return gerepilecopy(av,x);
+  modpr = nf_to_ff_init(nf, &pr,&T,&p);
+  a = modprM(lift(a), nf, modpr);
+  b = modprM(lift(b), nf, modpr);
+  return gerepilecopy(av, modprM_lift(FqM_gauss(a,b,T,p), modpr));
 }
 
+#if 0
 GEN
-nfsuppl(GEN nf, GEN x, long n, GEN prhall)
+nfsuppl(GEN nf, GEN x, GEN pr)
 {
-  long av=avma,av2,k,s,t,N,lx=lg(x);
-  GEN y,p1,p2,p,unmodp,zeromodp,unnf,zeronf,prh;
-  stackzone *zone;
+  gpmem_t av = avma;
+  GEN T,p,modpr;
 
-  k=lx-1; if (k>n) err(suppler2);
-  if (k && lg(x[1])!=n+1) err(talker,"incorrect dimension in nfsupl");
-  N=degpol(nf[1]); prh=(GEN)prhall[1]; p=gcoeff(prh,1,1);
-
-  zone  = switch_stack(NULL, 2*(3 + 2*lg(p) + N+1) + (n+3)*(n+1));
-  switch_stack(zone,1);
-  unmodp=gmodulsg(1,p); zeromodp=gmodulsg(0,p);
-  unnf=gscalcol_proto(unmodp,zeromodp,N);
-  zeronf=gscalcol_proto(zeromodp,zeromodp,N);
-  y = idmat_intern(n,unnf,zeronf);
-  switch_stack(zone,0); av2=avma;
-
-  for (s=1; s<=k; s++)
-  {
-    p1=nfsolvemodpr(nf,y,(GEN)x[s],prhall); t=s;
-    while (t<=n && gcmp0((GEN)p1[t])) t++;
-    avma=av2; if (t>n) err(suppler2);
-    p2=(GEN)y[s]; y[s]=x[s]; if (s!=t) y[t]=(long)p2;
-  }
-  avma=av; y=gcopy(y);
-  free(zone); return y;
+  nf = checknf(nf);
+  if (typ(x)!=t_MAT) err(typeer,"nfkermodpr");
+  modpr = nf_to_ff_init(nf, &pr,&T,&p);
+  x = modprM(lift(x), nf, modpr);
+  return gerepilecopy(av, modprM_lift(FqM_suppl(x,T,p), modpr));
 }
+#endif
 
-/* Given two fractional ideals a and b, gives x in a, y in b, z in b^-1,
-   t in a^-1 such that xt-yz=1. In the present version, z is in Z. */
-GEN
-nfidealdet1(GEN nf, GEN a, GEN b)
-{
-  long av=avma;
-  GEN x,p1,res,da,db;
-
-  a = idealinv(nf,a);
-  da = denom(a); if (!gcmp1(da)) a = gmul(da,a);
-  db = denom(b); if (!gcmp1(db)) b = gmul(db,b);
-  x = idealcoprime(nf,a,b);
-  p1 = idealaddtoone(nf, idealmul(nf,x,a), b);
-
-  res = cgetg(5,t_VEC);
-  res[1] = lmul(x,da);
-  res[2] = ldiv((GEN)p1[2],db);
-  res[3] = lnegi(db);
-  res[4] = (long) element_div(nf,(GEN)p1[1],(GEN)res[1]);
-  return gerepileupto(av,res);
-}
-
 /* Given a pseudo-basis pseudo, outputs a multiple of its ideal determinant */
 GEN
 nfdetint(GEN nf,GEN pseudo)
 {
   GEN pass,c,v,det1,piv,pivprec,vi,p1,x,I,unnf,zeronf,id,idprod;
-  long i,j,k,rg,n,n1,m,m1,av=avma,av1,tetpil,lim,cm=0,N;
+  long i, j, k, rg, n, n1, m, m1, cm=0, N;
+  gpmem_t av=avma, av1, tetpil, lim;
 
   nf=checknf(nf); N=degpol(nf[1]);
   if (typ(pseudo)!=t_VEC || lg(pseudo)!=3)
@@ -2999,58 +3146,6 @@ nfcleanmod(GEN nf, GEN x, long lim, GEN detmat)
     x[i]=(long)element_reduce(nf,(GEN)x[i],detmat);
 }
 
-static GEN
-zero_nfbezout(GEN nf,GEN b, GEN A,GEN B,GEN *u,GEN *v,GEN *w,GEN *di)
-{
-  long av, tetpil;
-  GEN pab,d;
-
-  d=idealmulelt(nf,b,B); *di=idealinv(nf,idealmat_to_hnf(nf,d));
-  av=avma; pab=idealmul(nf,A,B); tetpil=avma;
-  *w=gerepile(av,tetpil, idealmul(nf,pab,*di));
-  *v=element_inv(nf,b);
-  *u=gzero; return d;
-}
-
-/* Given elements a,b and ideals A, B, outputs d = a.A+b.B and gives
- * di=d^-1, w=A.B.di, u, v such that au+bv=1 and u in A.di, v in
- * B.di. Assume A, B non-zero, but a or b can be zero (not both)
- */
-static GEN
-nfbezout(GEN nf,GEN a,GEN b, GEN A,GEN B, GEN *u,GEN *v,GEN *w,GEN *di)
-{
-  GEN pab,pu,pv,pw,uv,d,dinv,pa,pb,pa1,pb1, *gptr[5];
-  long av,tetpil;
-
-  if (gcmp0(a))
-  {
-    if (gcmp0(b)) err(talker,"both elements zero in nfbezout");
-    return zero_nfbezout(nf,b,A,B,u,v,w,di);
-  }
-  if (gcmp0(b))
-    return zero_nfbezout(nf,a,B,A,v,u,w,di);
-
-  av = avma;
-  pa=idealmulelt(nf,a,A);
-  pb=idealmulelt(nf,b,B);
-
-  d=idealadd(nf,pa,pb); dinv=idealinv(nf,d);
-  pa1=idealmullll(nf,pa,dinv);
-  pb1=idealmullll(nf,pb,dinv);
-  uv=idealaddtoone(nf,pa1,pb1);
-  pab=idealmul(nf,A,B); tetpil=avma;
-
-  pu=element_div(nf,(GEN)uv[1],a);
-  pv=element_div(nf,(GEN)uv[2],b);
-  d=gcopy(d); dinv=gcopy(dinv);
-  pw=idealmul(nf,pab,dinv);
-
-  *u=pu; *v=pv; *w=pw; *di=dinv;
-  gptr[0]=u; gptr[1]=v; gptr[2]=w; gptr[3]=di;
-  gptr[4]=&d; gerepilemanysp(av,tetpil,gptr,5);
-  return d;
-}
-
 /* A usage interne. Pas de verifs ni gestion de pile */
 GEN
 idealoplll(GEN op(GEN,GEN,GEN), GEN nf, GEN x, GEN y)
@@ -3062,44 +3157,34 @@ idealoplll(GEN op(GEN,GEN,GEN), GEN nf, GEN x, GEN y)
   return den? gdiv(z,den): z;
 }
 
-/* A usage interne. Pas de verifs ni gestion de pile */
 GEN
-idealmulelt(GEN nf, GEN elt, GEN x)
-{
-  long t = typ(elt);
-  GEN z;
-  if (t == t_POL || t == t_POLMOD) elt = algtobasis(nf,elt);
-  if (isnfscalar(elt)) elt = (GEN)elt[1];
-  z = element_mulvec(nf, elt, x);
-  settyp(z, t_MAT); return z;
-}
-
-GEN
 nfhermitemod(GEN nf, GEN pseudo, GEN detmat)
 {
-  long av0=avma,li,co,av,tetpil,i,j,jm1,def,ldef,lim,N;
-  GEN b,q,w,p1,p2,d,u,v,den,x,I,J,dinv,unnf,wh;
+  long li, co, i, j, jm1, def, ldef, N;
+  gpmem_t av0=avma, av, lim;
+  GEN b,q,w,p1,d,u,v,den,x,I,J,dinv,wh,unnf;
 
   nf=checknf(nf); N=degpol(nf[1]);
   if (typ(pseudo)!=t_VEC || lg(pseudo)!=3)
     err(talker,"not a module in nfhermitemod");
-  x=(GEN)pseudo[1]; I=(GEN)pseudo[2];
-  if (typ(x)!=t_MAT) err(talker,"not a matrix in nfhermitemod");
-  co=lg(x);
-  if (typ(I)!=t_VEC || lg(I)!=co)
+  x = (GEN)pseudo[1];
+  I = (GEN)pseudo[2];
+  if (typ(x) != t_MAT) err(talker,"not a matrix in nfhermitemod");
+  co = lg(x);
+  if (typ(I) != t_VEC || lg(I) != co)
     err(talker,"not a correct ideal list in nfhermitemod");
   if (co==1) return cgetg(1,t_MAT);
 
-  li=lg(x[1]); x=dummycopy(x); I=dummycopy(I);
-  unnf=gscalcol_i(gun,N);
-  for (j=1; j<co; j++)
-    if (typ(I[j])!=t_MAT) I[j]=(long)idealhermite_aux(nf,(GEN)I[j]);
+  li= lg(x[1]);
+  x = matalgtobasis(nf, x);
+  I = dummycopy(I);
+  unnf = gscalcol(gun,N);
 
-  den=denom(detmat); if (!gcmp1(den)) detmat=gmul(den,detmat);
-  detmat=gmul(detmat,lllintpartial(detmat));
+  den = denom(detmat); if (!gcmp1(den)) detmat = gmul(den,detmat);
+  detmat = lllint_ip(detmat,100);
 
-  av=avma; lim=stack_lim(av,1);
-  def=co; ldef=(li>co)?li-co+1:1;
+  av = avma; lim = stack_lim(av,1);
+  def = co; ldef = (li>co)? li-co+1: 1;
   for (i=li-1; i>=ldef; i--)
   {
     def--; j=def-1; while (j && gcmp0(gcoeff(x,i,j))) j--;
@@ -3124,39 +3209,40 @@ nfhermitemod(GEN nf, GEN pseudo, GEN detmat)
     }
     if (low_stack(lim, stack_lim(av,1)))
     {
-      GEN *gptr[2];
+      GEN *gptr[2]; gptr[0]=&x; gptr[1]=&I;
       if(DEBUGMEM>1) err(warnmem,"[1]: nfhermitemod");
-      gptr[0]=&x; gptr[1]=&I; gerepilemany(av,gptr,2);
+      gerepilemany(av,gptr,2);
     }
   }
-  b=detmat; wh=cgetg(li,t_MAT); def--;
+  b = detmat; wh = cgetg(li,t_MAT); def--;
   for (i=li-1; i>=1; i--)
   {
     d = nfbezout(nf,gcoeff(x,i,i+def),unnf,(GEN)I[i+def],b,&u,&v,&w,&dinv);
     p1 = element_mulvec(nf,u,(GEN)x[i+def]);
     nfcleanmod(nf,p1,i,idealmullll(nf,b,dinv));
-    wh[i]=(long)p1; coeff(wh,i,i)=(long)unnf; I[i+def]=(long)d;
-    if (i>1) b=idealmul(nf,b,dinv);
+    wh[i] = (long)p1; coeff(wh,i,i) = (long)unnf;
+    I[i+def] = (long)d;
+    if (i>1) b = idealmul(nf,b,dinv);
   }
   J=cgetg(li,t_VEC); J[1]=zero;
-  for (j=2; j<li; j++) J[j]=(long)idealinv(nf,(GEN)I[j+def]);
+  for (j=2; j<li; j++) J[j] = (long)idealinv(nf,(GEN)I[j+def]);
   for (i=li-2; i>=1; i--)
   {
     for (j=i+1; j<li; j++)
     {
-      q=idealmul(nf,(GEN)I[i+def],(GEN)J[j]);
-      p1=gsub(element_reduce(nf,gcoeff(wh,i,j),q),gcoeff(wh,i,j));
-      wh[j]=(long)gadd((GEN)wh[j],element_mulvec(nf,p1,(GEN)wh[i]));
+      q = idealmul(nf,(GEN)I[i+def],(GEN)J[j]);
+      p1 = gsub(element_reduce(nf,gcoeff(wh,i,j),q), gcoeff(wh,i,j));
+      wh[j] = ladd((GEN)wh[j],element_mulvec(nf,p1,(GEN)wh[i]));
     }
     if (low_stack(lim, stack_lim(av,1)))
     {
-      GEN *gptr[3];
+      GEN *gptr[3]; gptr[0]=&wh; gptr[1]=&I; gptr[2]=&J;
       if(DEBUGMEM>1) err(warnmem,"[2]: nfhermitemod");
-      gptr[0]=&wh; gptr[1]=&I; gptr[2]=&J; gerepilemany(av,gptr,3);
+      gerepilemany(av,gptr,3);
     }
   }
-  tetpil=avma; p1=cgetg(3,t_VEC); p1[1]=lcopy(wh);
-  p2=cgetg(li,t_VEC); p1[2]=(long)p2;
-  for (j=1; j<li; j++) p2[j]=lcopy((GEN)I[j+def]);
-  return gerepile(av0,tetpil,p1);
+  I += def; I[0] = evaltyp(t_VEC)|evallg(li);
+  p1=cgetg(3,t_VEC);
+  p1[1] = (long)wh;
+  p1[2] = (long)I; return gerepilecopy(av0, p1);
 }