===================================================================
RCS file: /home/cvs/OpenXM_contrib/pari-2.2/src/basemath/Attic/base2.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/base2.c	2001/10/02 11:17:02	1.1
+++ OpenXM_contrib/pari-2.2/src/basemath/Attic/base2.c	2002/09/11 07:26:49	1.2
@@ -1,4 +1,4 @@
-/* $Id: base2.c,v 1.1 2001/10/02 11:17:02 noro Exp $
+/* $Id: base2.c,v 1.2 2002/09/11 07:26:49 noro Exp $
 
 Copyright (C) 2000  The PARI group.
 
@@ -19,42 +19,62 @@ Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 /*                                                                 */
 /*******************************************************************/
 #include "pari.h"
+#include "parinf.h"
 
-extern GEN caractducos(GEN p, GEN x, int v);
-extern GEN element_muli(GEN nf, GEN x, GEN y);
-extern GEN element_mulid(GEN nf, GEN x, long i);
-extern GEN eleval(GEN f,GEN h,GEN a);
-extern GEN ideal_better_basis(GEN nf, GEN x, GEN M);
-extern long int_elt_val(GEN nf, GEN x, GEN p, GEN bp, GEN *t, long v);
+#define RXQX_rem(x,y,T) RXQX_divrem((x),(y),(T),ONLY_REM)
+#define FpX_rem(x,y,p) FpX_divres((x),(y),(p),ONLY_REM)
+extern GEN addone_aux2(GEN x, GEN y);
+extern GEN addshiftw(GEN x, GEN y, long d);
+extern GEN gmul_mat_smallvec(GEN x, GEN y);
+extern GEN hnf_invimage(GEN A, GEN b);
+extern GEN norm_by_embed(long r1, GEN x);
+extern GEN ZX_resultant_all(GEN A, GEN B, GEN dB, ulong bound);
+extern GEN polsym_gen(GEN P, GEN y0, long n, GEN T, GEN N);
+extern GEN FqX_gcd(GEN P, GEN Q, GEN T, GEN p);
+extern GEN FqX_factor(GEN x, GEN T, GEN p);
+extern GEN DDF2(GEN x, long hint);
+extern GEN eltabstorel(GEN x, GEN T, GEN pol, GEN k);
+extern GEN element_powid_mod_p(GEN nf, long I, GEN n, GEN p);
+extern GEN FpVQX_red(GEN z, GEN T, GEN p);
+extern GEN Fp_factor_irred(GEN P,GEN l, GEN Q);
+extern GEN RXQX_divrem(GEN x, GEN y, GEN T, GEN *pr);
+extern GEN RXQX_mul(GEN x, GEN y, GEN T);
+extern GEN ZY_ZXY_resultant_all(GEN A, GEN B0, long *lambda, GEN *LPRS);
+extern GEN _ei(long n, long i);
+extern GEN col_to_ff(GEN x, long v);
+extern GEN element_mulidid(GEN nf, long i, long j);
+extern GEN eltmulid_get_table(GEN nf, long i);
+extern GEN idealaddtoone_i(GEN nf, GEN x, GEN y);
 extern GEN mat_to_vecpol(GEN x, long v);
-extern GEN nfidealdet1(GEN nf, GEN a, GEN b);
-extern GEN nfsuppl(GEN nf, GEN x, long n, GEN prhall);
+extern GEN merge_factor_i(GEN f, GEN g);
+extern GEN mulmat_pol(GEN A, GEN x);
+extern GEN nfgcd(GEN P, GEN Q, GEN nf, GEN den);
+extern GEN pidealprimeinv(GEN nf, GEN x);
 extern GEN pol_to_monic(GEN pol, GEN *lead);
-extern GEN pol_to_vec(GEN x, long N);
 extern GEN quicktrace(GEN x, GEN sym);
-extern GEN respm(GEN f1,GEN f2,GEN pm);
+extern GEN sqr_by_tab(GEN tab, GEN x);
+extern GEN to_polmod(GEN x, GEN mod);
+extern GEN unnf_minus_x(GEN x);
+extern int isrational(GEN x);
+extern long int_elt_val(GEN nf, GEN x, GEN p, GEN bp, GEN *t);
 
+/* FIXME: backward compatibility. Should use the proper nf_* equivalents */
+#define compat_PARTIAL 1
+#define compat_ROUND2  2
+
 static void
-allbase_check_args(GEN f, long code, GEN *y, GEN *ptw1, GEN *ptw2)
+allbase_check_args(GEN f, long flag, GEN *dx, GEN *ptw)
 {
-  GEN w;
+  GEN w = *ptw;
+
+  if (DEBUGLEVEL) (void)timer2();
   if (typ(f)!=t_POL) err(notpoler,"allbase");
   if (degpol(f) <= 0) err(constpoler,"allbase");
-  if (DEBUGLEVEL) timer2();
-  switch(code)
-  {
-    case 0: case 1:
-      *y = ZX_disc(f);
-      if (!signe(*y)) err(talker,"reducible polynomial in allbase");
-      w = auxdecomp(absi(*y),1-code);
-      break;
-    default:
-      w = (GEN)code;
-      *y = factorback(w, NULL);
-  }
+
+  *dx = w? factorback(w, NULL): ZX_disc(f);
+  if (!signe(*dx)) err(talker,"reducible polynomial in allbase");
+  if (!w) *ptw = auxdecomp(absi(*dx), (flag & nf_PARTIALFACT)? 0: 1);
   if (DEBUGLEVEL) msgtimer("disc. factorisation");
-  *ptw1 = (GEN)w[1];
-  *ptw2 = (GEN)w[2];
 }
 
 /*******************************************************************/
@@ -62,39 +82,6 @@ allbase_check_args(GEN f, long code, GEN *y, GEN *ptw1
 /*                            ROUND 2                              */
 /*                                                                 */
 /*******************************************************************/
-/*  Normalized quotient and remainder ( -1/2 |y| < r = x-q*y <= 1/2 |y| ) */
-static GEN
-rquot(GEN x, GEN y)
-{
-  long av=avma,av1;
-  GEN u,v,w,p;
-
-  u=absi(y); v=shifti(x,1); w=shifti(y,1);
-  if (cmpii(u,v)>0) p=subii(v,u);
-  else p=addsi(-1,addii(u,v));
-  av1=avma; return gerepile(av,av1,divii(p,w));
-}
-
-/* space needed lx + 2*ly */
-static GEN
-rrmdr(GEN x, GEN y)
-{
-  long av=avma,tetpil,k;
-  GEN r,ys2;
-
-  if (!signe(x)) return gzero;
-  r = resii(x,y); tetpil = avma;
-  ys2 = shifti(y,-1);
-  k = absi_cmp(r, ys2);
-  if (k>0 || (k==0 && signe(r)>0))
-  {
-    avma = tetpil;
-    if (signe(y) == signe(r)) r = subii(r,y); else r = addii(r,y);
-    return gerepile(av,tetpil,r);
-  }
-  avma = tetpil; return r;
-}
-
 /* companion matrix of unitary polynomial x */
 static GEN
 companion(GEN x) /* cf assmat */
@@ -117,7 +104,8 @@ companion(GEN x) /* cf assmat */
 static GEN
 mulmati(GEN x, GEN y)
 {
-  long n = lg(x),i,j,k,av;
+  gpmem_t av;
+  long n = lg(x),i,j,k;
   GEN z = cgetg(n,t_MAT),p1,p2;
 
   for (j=1; j<n; j++)
@@ -138,24 +126,26 @@ mulmati(GEN x, GEN y)
 }
 
 static GEN
-powmati(GEN x, long m)
-{
-  long av=avma,j;
-  GEN y = x;
+_mulmati(void *data /*ignored*/, GEN x, GEN y) {
+  (void)data; return mulmati(x,y);
+}
+static GEN
+_sqrmati(void *data /*ignored*/, GEN x) {
+  (void)data; return mulmati(x,x);
+}
 
-  j=1+bfffo(m); m<<=j; j = BITS_IN_LONG-j;
-  for (; j; m<<=1,j--)
-  {
-    y=mulmati(y,y);
-    if (m<0) y=mulmati(y,x);
-  }
+static GEN
+powmati(GEN x, GEN n)
+{
+  gpmem_t av = avma;
+  GEN y = leftright_pow(x, n, NULL, &_sqrmati, &_mulmati);
   return gerepileupto(av,y);
 }
 
 static GEN
 rtran(GEN v, GEN w, GEN q)
 {
-  long av,tetpil;
+  gpmem_t av,tetpil;
   GEN p1;
 
   if (signe(q))
@@ -169,21 +159,19 @@ rtran(GEN v, GEN w, GEN q)
 /* return (v - qw) mod m (only compute entries k0,..,n)
  * v and w are expected to have entries smaller than m */
 static GEN
-mtran(GEN v, GEN w, GEN q, GEN m, long k0)
+mtran(GEN v, GEN w, GEN q, GEN m, GEN mo2, long k0)
 {
-  long k,l;
+  long k;
   GEN p1;
 
   if (signe(q))
-  {
-    l = lgefint(m) << 2;
-    for (k=lg(v)-1; k>= k0; k--)
+    for (k=lg(v)-1; k >= k0; k--)
     {
-      long av = avma; (void)new_chunk(l);
+      gpmem_t av = avma;
       p1 = subii((GEN)v[k], mulii(q,(GEN)w[k]));
-      avma = av; v[k]=(long)rrmdr(p1, m);
+      p1 = centermodii(p1, m, mo2);
+      v[k] = lpileuptoint(av, p1);
     }
-  }
   return v;
 }
 
@@ -236,24 +224,24 @@ static void
 rowred(GEN a, GEN rmod)
 {
   long j,k,pro, c = lg(a), r = lg(a[1]);
-  long av=avma, lim=stack_lim(av,1);
-  GEN q;
+  gpmem_t av=avma, lim=stack_lim(av,1);
+  GEN q, rmodo2 = shifti(rmod,-1);
 
   for (j=1; j<r; j++)
   {
     for (k=j+1; k<c; k++)
       while (signe(gcoeff(a,j,k)))
       {
-	q=rquot(gcoeff(a,j,j),gcoeff(a,j,k));
-	pro=(long)mtran((GEN)a[j],(GEN)a[k],q,rmod, j);
+	q=diviiround(gcoeff(a,j,j),gcoeff(a,j,k));
+	pro=(long)mtran((GEN)a[j],(GEN)a[k],q,rmod,rmodo2, j);
 	a[j]=a[k]; a[k]=pro;
       }
     if (signe(gcoeff(a,j,j)) < 0)
       for (k=j; k<r; k++) coeff(a,k,j)=lnegi(gcoeff(a,k,j));
     for (k=1; k<j; k++)
     {
-      q=rquot(gcoeff(a,j,k),gcoeff(a,j,j));
-      a[k]=(long)mtran((GEN)a[k],(GEN)a[j],q,rmod, k);
+      q=diviiround(gcoeff(a,j,k),gcoeff(a,j,j));
+      a[k]=(long)mtran((GEN)a[k],(GEN)a[j],q,rmod,rmodo2, k);
     }
     if (low_stack(lim, stack_lim(av,1)))
     {
@@ -268,19 +256,19 @@ rowred(GEN a, GEN rmod)
 }
 
 /* Calcule d/x  ou  d est entier et x matrice triangulaire inferieure
- * entiere dont les coeff diagonaux divisent d (resultat entier).
- */
+ * entiere dont les coeff diagonaux divisent d (resultat entier). */
 static GEN
-matinv(GEN x, GEN d, long n)
+matinv(GEN x, GEN d)
 {
-  long i,j,k,av,av1;
+  gpmem_t av,av1;
+  long i,j,k, n = lg(x[1]); /* Warning: lg(x) from ordmax is bogus */
   GEN y,h;
 
-  y=idmat(n);
-  for (i=1; i<=n; i++)
-    coeff(y,i,i)=ldivii(d,gcoeff(x,i,i));
+  y = idmat(n-1);
+  for (i=1; i<n; i++)
+    coeff(y,i,i) = (long)diviiexact(d,gcoeff(x,i,i));
   av=avma;
-  for (i=2; i<=n; i++)
+  for (i=2; i<n; i++)
     for (j=i-1; j; j--)
     {
       for (h=gzero,k=j+1; k<=i; k++)
@@ -298,14 +286,16 @@ matinv(GEN x, GEN d, long n)
 static GEN
 ordmax(GEN *cf, GEN p, long epsilon, GEN *ptdelta)
 {
-  long sp,hard_case_exponent,i,n=lg(cf)-1,av=avma, av2,limit;
-  GEN T,T2,Tn,m,v,delta, *w;
+  long sp,i,n=lg(cf)-1;
+  gpmem_t av=avma, av2,limit;
+  GEN T,T2,Tn,m,v,delta,hard_case_exponent, *w;
   const GEN pp = sqri(p);
+  const GEN ppo2 = shifti(pp,-1);
   const long pps = (2*expi(pp)+2<BITS_IN_LONG)? pp[2]: 0;
 
   if (cmpis(p,n) > 0)
   {
-    hard_case_exponent = 0;
+    hard_case_exponent = NULL;
     sp = 0; /* gcc -Wall */
   }
   else
@@ -313,26 +303,28 @@ ordmax(GEN *cf, GEN p, long epsilon, GEN *ptdelta)
     long k;
     k = sp = itos(p);
     i=1; while (k < n) { k *= sp; i++; }
-    hard_case_exponent = i;
+    hard_case_exponent = stoi(i);
   }
   T=cgetg(n+1,t_MAT); for (i=1; i<=n; i++) T[i]=lgetg(n+1,t_COL);
   T2=cgetg(2*n+1,t_MAT); for (i=1; i<=2*n; i++) T2[i]=lgetg(n+1,t_COL);
   Tn=cgetg(n*n+1,t_MAT); for (i=1; i<=n*n; i++) Tn[i]=lgetg(n+1,t_COL);
   v = new_chunk(n+1);
-  w =  (GEN*)new_chunk(n+1);
+  w = (GEN*)new_chunk(n+1);
 
   av2 = avma; limit = stack_lim(av2,1);
   delta=gun; m=idmat(n);
 
   for(;;)
   {
-    long j,k,h, av0 = avma;
+    long j, k, h;
+    gpmem_t av0 = avma;
     GEN t,b,jp,hh,index,p1, dd = sqri(delta), ppdd = mulii(dd,pp);
+    GEN ppddo2 = shifti(ppdd,-1);
 
     if (DEBUGLEVEL > 3)
       fprintferr("ROUND2: epsilon = %ld\tavma = %ld\n",epsilon,avma);
 
-    b=matinv(m,delta,n);
+    b=matinv(m,delta);
     for (i=1; i<=n; i++)
     {
       for (j=1; j<=n; j++)
@@ -344,12 +336,12 @@ ordmax(GEN *cf, GEN p, long epsilon, GEN *ptdelta)
 	    GEN p2 = mulii(gcoeff(m,i,h),gcoeff(cf[h],j,k));
             if (p2!=gzero) p1 = addii(p1,p2);
           }
-          coeff(T,j,k) = (long)rrmdr(p1, ppdd);
+          coeff(T,j,k) = (long)centermodii(p1, ppdd, ppddo2);
         }
       p1 = mulmati(m, mulmati(T,b));
       for (j=1; j<=n; j++)
 	for (k=1; k<=n; k++)
-	  coeff(p1,j,k)=(long)rrmdr(divii(gcoeff(p1,j,k),dd),pp);
+	  coeff(p1,j,k)=(long)centermodii(divii(gcoeff(p1,j,k),dd),pp,ppo2);
       w[i] = p1;
     }
 
@@ -381,14 +373,14 @@ ordmax(GEN *cf, GEN p, long epsilon, GEN *ptdelta)
       for (i=1; i<=n; i++)
 	for (j=1; j<=n; j++)
 	{
-          long av1 = avma;
+          gpmem_t av1 = avma;
           p1 = gzero;
 	  for (k=1; k<=n; k++)
 	    for (h=1; h<=n; h++)
 	    {
 	      const GEN r=modii(gcoeff(w[i],k,h),p);
 	      const GEN s=modii(gcoeff(w[j],h,k),p);
-              const GEN p2 = mulii(r,s); 
+              const GEN p2 = mulii(r,s);
 	      if (p2!=gzero) p1 = addii(p1,p2);
 	    }
 	  coeff(T,i,j) = lpileupto(av1,p1);
@@ -405,7 +397,7 @@ ordmax(GEN *cf, GEN p, long epsilon, GEN *ptdelta)
           coeff(T2,j,i)=(i==j)? ps: 0;
           coeff(T2,j,n+i)=smodis(gcoeff(t,i,j),ps);
         }
-      rowred_long(T2,pps); 
+      rowred_long(T2,pps);
     }
     else
     {
@@ -417,12 +409,12 @@ ordmax(GEN *cf, GEN p, long epsilon, GEN *ptdelta)
         }
       rowred(T2,pp);
     }
-    jp=matinv(T2,p,n);
+    jp=matinv(T2,p);
     if (pps)
     {
       for (k=1; k<=n; k++)
       {
-        long av1=avma;
+        gpmem_t av1=avma;
         t = mulmati(mulmati(jp,w[k]), T2);
         for (h=i=1; i<=n; i++)
           for (j=1; j<=n; j++)
@@ -447,7 +439,7 @@ ordmax(GEN *cf, GEN p, long epsilon, GEN *ptdelta)
       index = mulii(index,gcoeff(Tn,i,i));
     if (gcmp1(index)) break;
 
-    m = mulmati(matinv(Tn,index,n), m);
+    m = mulmati(matinv(Tn,index), m);
     hh = delta = mulii(index,delta);
     for (i=1; i<=n; i++)
       for (j=1; j<=n; j++)
@@ -486,13 +478,17 @@ ordmax(GEN *cf, GEN p, long epsilon, GEN *ptdelta)
  *
  *  2) discriminant of K (in *y).
  */
-GEN
-allbase(GEN f, long code, GEN *y)
+static GEN
+allbase2(GEN f, int flag, GEN *dx, GEN *dK, GEN *ptw)
 {
-  GEN w1,w2,a,pro,at,bt,b,da,db,q, *cf,*gptr[2];
-  long av=avma,tetpil,n,h,j,i,k,r,s,t,v,mf;
+  GEN w,w1,w2,a,pro,at,bt,b,da,db,q, *cf,*gptr[2];
+  gpmem_t av=avma,tetpil;
+  long n,h,j,i,k,r,s,t,v,mf;
 
-  allbase_check_args(f,code,y, &w1,&w2);
+  w = ptw? *ptw: NULL;
+  allbase_check_args(f,flag,dx, &w);
+  w1 = (GEN)w[1];
+  w2 = (GEN)w[2];
   v = varn(f); n = degpol(f); h = lg(w1)-1;
   cf = (GEN*)cgetg(n+1,t_VEC);
   cf[2]=companion(f);
@@ -501,7 +497,7 @@ allbase(GEN f, long code, GEN *y)
   a=idmat(n); da=gun;
   for (i=1; i<=h; i++)
   {
-    long av1 = avma;
+    gpmem_t av1 = avma;
     mf=itos((GEN)w2[i]); if (mf==1) continue;
     if (DEBUGLEVEL) fprintferr("Treating p^k = %Z^%ld\n",w1[i],mf);
 
@@ -513,11 +509,11 @@ allbase(GEN f, long code, GEN *y)
       for (s=r; s; s--)
         while (signe(gcoeff(bt,s,r)))
         {
-          q=rquot(gcoeff(at,s,s),gcoeff(bt,s,r));
+          q=diviiround(gcoeff(at,s,s),gcoeff(bt,s,r));
           pro=rtran((GEN)at[s],(GEN)bt[r],q);
           for (t=s-1; t; t--)
           {
-            q=rquot(gcoeff(at,t,s),gcoeff(at,t,t));
+            q=diviiround(gcoeff(at,t,s),gcoeff(at,t,t));
             pro=rtran(pro,(GEN)at[t],q);
           }
           at[s]=bt[r]; bt[r]=(long)pro;
@@ -528,7 +524,7 @@ allbase(GEN f, long code, GEN *y)
       {
         while (signe(gcoeff(at,j,k)))
         {
-          q=rquot(gcoeff(at,j,j),gcoeff(at,j,k));
+          q=diviiround(gcoeff(at,j,j),gcoeff(at,j,k));
           pro=rtran((GEN)at[j],(GEN)at[k],q);
           at[j]=at[k]; at[k]=(long)pro;
         }
@@ -537,7 +533,7 @@ allbase(GEN f, long code, GEN *y)
         for (k=1; k<=j; k++) coeff(at,k,j)=lnegi(gcoeff(at,k,j));
       for (k=j+1; k<=n; k++)
       {
-        q=rquot(gcoeff(at,j,k),gcoeff(at,j,j));
+        q=diviiround(gcoeff(at,j,k),gcoeff(at,j,j));
         at[k]=(long)rtran((GEN)at[k],(GEN)at[j],q);
       }
     }
@@ -549,42 +545,33 @@ allbase(GEN f, long code, GEN *y)
       }
     tetpil=avma; a=gtrans(at);
     {
-      GEN *gptr[2]; 
+      GEN *gptr[2];
       da = icopy(da); gptr[0]=&a; gptr[1]=&da;
       gerepilemanysp(av1,tetpil,gptr,2);
     }
   }
+  *dK = *dx;
   for (j=1; j<=n; j++)
-    *y = divii(mulii(*y,sqri(gcoeff(a,j,j))), sqri(da));
-  tetpil=avma; *y=icopy(*y);
+    *dK = diviiexact(mulii(*dK,sqri(gcoeff(a,j,j))), sqri(da));
+  tetpil=avma; *dK = icopy(*dK);
   at=cgetg(n+1,t_VEC); v=varn(f);
   for (k=1; k<=n; k++)
   {
     q=cgetg(k+2,t_POL); at[k]=(long)q;
     q[1] = evalsigne(1) | evallgef(2+k) | evalvarn(v);
-    for (j=1; j<=k; j++) q[j+1]=ldiv(gcoeff(a,k,j),da);
+    for (j=1; j<=k; j++) q[j+1] = ldiv(gcoeff(a,k,j),da);
   }
-  gptr[0]=&at; gptr[1]=y;
+  gptr[0] = &at; gptr[1] = dK;
   gerepilemanysp(av,tetpil,gptr,2);
   return at;
 }
 
 GEN
-base2(GEN x, GEN *y)
-{
-  return allbase(x,0,y);
-}
+base2(GEN x, GEN *pdK) { return nfbasis(x, pdK, compat_ROUND2, NULL); }
 
 GEN
-discf2(GEN x)
-{
-  GEN y;
-  long av=avma,tetpil;
+discf2(GEN x) { return nfdiscf0(x, compat_ROUND2, NULL); }
 
-  allbase(x,0,&y); tetpil=avma;
-  return gerepile(av,tetpil,icopy(y));
-}
-
 /*******************************************************************/
 /*                                                                 */
 /*                            ROUND 4                              */
@@ -595,7 +582,6 @@ GEN nilord(GEN p,GEN fx,long mf,GEN gx,long flag);
 GEN Decomp(GEN p,GEN f,long mf,GEN theta,GEN chi,GEN nu,long r);
 static GEN dbasis(GEN p, GEN f, long mf, GEN alpha, GEN U);
 static GEN maxord(GEN p,GEN f,long mf);
-static GEN nbasis(GEN ibas,GEN pd);
 static GEN testb2(GEN p,GEN fa,long Fa,GEN theta,GEN pmf,long Ft,GEN ns);
 static GEN testc2(GEN p,GEN fa,GEN pmr,GEN pmf,GEN alph2,
 		  long Ea,GEN thet2,long Et,GEN ns);
@@ -609,107 +595,158 @@ fnz(GEN x,long j)
   return 1;
 }
 
-/* retourne la base, dans y le discf et dans ptw la factorisation (peut
- etre partielle) de discf */
+/* return list u[i], 2 by 2 coprime with the same prime divisors as ab */
+static GEN
+get_coprimes(GEN a, GEN b)
+{
+  long i, k = 1;
+  GEN u = cgetg(3, t_VEC);
+  u[1] = (long)a;
+  u[2] = (long)b;
+  /* u1,..., uk 2 by 2 coprime */
+  while (k+1 < lg(u))
+  {
+    GEN d, c = (GEN)u[k+1];
+    if (is_pm1(c)) { k++; continue; }
+    for (i=1; i<=k; i++)
+    {
+      if (is_pm1(u[i])) continue;
+      d = mppgcd(c, (GEN)u[i]);
+      if (d == gun) continue;
+      c = diviiexact(c, d);
+      u[i] = (long)diviiexact((GEN)u[i], d);
+      u = concatsp(u, d);
+    }
+    u[++k] = (long)c;
+  }
+  for (i = k = 1; i < lg(u); i++)
+    if (!is_pm1(u[i])) u[k++] = u[i];
+  setlg(u, k); return u;
+}
+
+/* return integer basis. Set dK = disc(K), dx = disc(f), w (possibly partial)
+ * factorization of dK. *ptw can be set by the caller, in which case it is
+ * taken to be the factorization of disc(f), then overwritten
+ * [no consistency check] */
 GEN
-allbase4(GEN f,long code, GEN *y, GEN *ptw)
+allbase(GEN f, int flag, GEN *dx, GEN *dK, GEN *index, GEN *ptw)
 {
-  GEN w,w1,w2,a,da,b,db,bas,q,p1,*gptr[3];
-  long v,n,mf,h,lfa,i,j,k,l,tetpil,av = avma;
+  GEN w, w1, w2, a, da, p1, ordmax;
+  long n, lw, i, j, k, l;
 
-  allbase_check_args(f,code,y, &w1,&w2);
-  v = varn(f); n = degpol(f); h = lg(w1)-1;
-  a = NULL; /* gcc -Wall */
-  da= NULL;
-  for (i=1; i<=h; i++)
+  if (flag & nf_ROUND2) return allbase2(f,flag,dx,dK,ptw);
+  w = ptw? *ptw: NULL;
+  allbase_check_args(f, flag, dx, &w);
+  w1 = (GEN)w[1];
+  w2 = (GEN)w[2];
+  n = degpol(f); lw = lg(w1);
+  ordmax = cgetg(1, t_VEC);
+  /* "complete" factorization first */
+  for (i=1; i<lw; i++)
   {
-    mf=itos((GEN)w2[i]); if (mf == 1) continue;
-    if (DEBUGLEVEL) fprintferr("Treating p^k = %Z^%ld\n",w1[i],mf);
+    long mf = itos((GEN)w2[i]);
+    if (mf == 1) { ordmax = concatsp(ordmax, gun); continue; }
 
-    b = maxord((GEN)w1[i],f,mf); db = gun;
+    CATCH(invmoder) { /* caught false prime, update factorization */
+      GEN x = (GEN)global_err_data;
+      GEN p = mppgcd((GEN)x[1], (GEN)x[2]);
+      GEN N, u;
+      if (DEBUGLEVEL) err(warner,"impossible inverse: %Z", x);
+
+      u = get_coprimes(p, diviiexact((GEN)x[1],p));
+      l = lg(u);
+      /* no small factors, but often a prime power */
+      for (k = 1; k < l; k++) u[k] = coeff(auxdecomp((GEN)u[k], 2),1,1);
+
+      w1[i] = u[1];
+      w1 = concatsp(w1, vecextract_i(u, 2, l-1));
+      N = *dx;
+      w2[i] = lstoi(pvaluation(N, (GEN)w1[i], &N));
+      k  = lw;
+      lw = lg(w1);
+      for ( ; k < lw; k++) w2[k] = lstoi(pvaluation(N, (GEN)w1[k], &N));
+    } RETRY {
+      if (DEBUGLEVEL) fprintferr("Treating p^k = %Z^%ld\n",w1[i],mf);
+      ordmax = concatsp(ordmax, _vec( maxord((GEN)w1[i],f,mf) ));
+    } ENDCATCH;
+  }
+
+  a = NULL; /* gcc -Wall */
+  da = NULL;
+  for (i=1; i<lw; i++)
+  {
+    GEN db, b = (GEN)ordmax[i];
+    if (b == gun) continue;
+    db = gun;
     for (j=1; j<=n; j++)
     {
       p1 = denom(gcoeff(b,j,j));
       if (cmpii(p1,db) > 0) db = p1;
     }
-    if (db != gun)
-    { /* db = denom(diag(b)), (da,db) = 1 */
-      b = gmul(b,db);
-      if (!da) { da=db; a=b; }
-      else
+    if (db == gun) continue;
+
+    /* db = denom(diag(b)), (da,db) = 1 */
+    b = Q_muli_to_int(b,db);
+    if (!da) { da = db; a = b; }
+    else
+    {
+      j=1; while (j<=n && fnz((GEN)a[j],j) && fnz((GEN)b[j],j)) j++;
+      b = gmul(da,b);
+      a = gmul(db,a); da = mulii(da,db);
+      k = j-1; p1 = cgetg(2*n-k+1,t_MAT);
+      for (j=1; j<=k; j++)
       {
-        j=1; while (j<=n && fnz((GEN)a[j],j) && fnz((GEN)b[j],j)) j++;
-        b = gmul(da,b); a = gmul(db,a);
-        k=j-1; p1=cgetg(2*n-k+1,t_MAT);
-        for (j=1; j<=k; j++)
-        {
-          p1[j] = a[j];
-          coeff(p1,j,j) = lmppgcd(gcoeff(a,j,j),gcoeff(b,j,j));
-        }
-        for (  ; j<=n;     j++) p1[j] = a[j];
-        for (  ; j<=2*n-k; j++) p1[j] = b[j+k-n];
-        da = mulii(da,db); a = hnfmodid(p1, da);
+        p1[j] = a[j];
+        coeff(p1,j,j) = lmppgcd(gcoeff(a,j,j),gcoeff(b,j,j));
       }
+      for (  ; j<=n;     j++) p1[j] = a[j];
+      for (  ; j<=2*n-k; j++) p1[j] = b[j+k-n];
+      a = hnfmodid(p1, da);
     }
-    if (DEBUGLEVEL>5)
-      fprintferr("Result for prime %Z is:\n%Z\n",w1[i],b);
+    if (DEBUGLEVEL>5) fprintferr("Result for prime %Z is:\n%Z\n",w1[i],b);
   }
+  *dK = *dx;
   if (da)
   {
-    for (j=1; j<=n; j++)
-      *y = mulii(divii(*y,sqri(da)),sqri(gcoeff(a,j,j)));
+    *index = diviiexact(da, gcoeff(a,1,1));
+    for (j=2; j<=n; j++) *index = mulii(*index, diviiexact(da, gcoeff(a,j,j)));
+    *dK = diviiexact(*dK, sqri(*index));
     for (j=n-1; j; j--)
       if (cmpis(gcoeff(a,j,j),2) > 0)
       {
-        p1=shifti(gcoeff(a,j,j),-1);
+        p1 = shifti(gcoeff(a,j,j),-1);
         for (k=j+1; k<=n; k++)
           if (cmpii(gcoeff(a,j,k),p1) > 0)
             for (l=1; l<=j; l++)
-              coeff(a,l,k)=lsubii(gcoeff(a,l,k),gcoeff(a,l,j));
+              coeff(a,l,k) = lsubii(gcoeff(a,l,k),gcoeff(a,l,j));
       }
+    a = gdiv(a, da);
   }
-  lfa = 0;
-  if (ptw)
+  else
   {
-    for (j=1; j<=h; j++)
-    {
-      k=ggval(*y,(GEN)w1[j]);
-      if (k) { lfa++; w1[lfa]=w1[j]; w2[lfa]=k; }
-    }
+    *index = gun;
+    a = idmat(n);
   }
-  tetpil=avma; *y=icopy(*y);
-  bas=cgetg(n+1,t_VEC); v=varn(f);
-  for (k=1; k<=n; k++)
-  {
-    q=cgetg(k+2,t_POL); bas[k]=(long)q;
-    q[1] = evalsigne(1) | evallgef(k+2) | evalvarn(v);
-    if (da)
-      for (j=1; j<=k; j++) q[j+1]=ldiv(gcoeff(a,j,k),da);
-    else
-    {
-      for (j=2; j<=k; j++) q[j]=zero;
-      q[j]=un;
-    }
-  }
+
   if (ptw)
   {
-    *ptw=w=cgetg(3,t_MAT);
-    w[1]=lgetg(lfa+1,t_COL);
-    w[2]=lgetg(lfa+1,t_COL);
-    for (j=1; j<=lfa; j++)
+    long lfa = 1;
+    GEN W1, W2, D = *dK;
+
+    w = cgetg(3,t_MAT);
+    W1 = cgetg(lw, t_COL); w[1] = (long)W1;
+    W2 = cgetg(lw, t_COL); w[2] = (long)W2;
+    for (j=1; j<lw; j++)
     {
-      coeff(w,j,1)=(long)icopy((GEN)w1[j]);
-      coeff(w,j,2)=lstoi(w2[j]);
+      k = pvaluation(D, (GEN)w1[j], &D);
+      if (k) { W1[lfa] = w1[j]; W2[lfa] = lstoi(k); lfa++; }
     }
-    gptr[2]=ptw;
+    setlg(W1, lfa);
+    setlg(W2, lfa); *ptw = w;
   }
-  gptr[0]=&bas; gptr[1]=y;
-  gerepilemanysp(av,tetpil,gptr, ptw?3:2);
-  return bas;
+  return mat_to_vecpol(a, varn(f));
 }
 
-extern GEN merge_factor_i(GEN f, GEN g);
-
 static GEN
 update_fact(GEN x, GEN f)
 {
@@ -723,127 +760,86 @@ update_fact(GEN x, GEN f)
   for (i=1; i<l; i++)
   {
     k = pvaluation(d, (GEN)p[i], &d);
-    if (k) { q[iq] = p[i]; e[iq] = lstoi(k); iq++; } 
+    if (k) { q[iq] = p[i]; e[iq] = lstoi(k); iq++; }
   }
   setlg(q,iq); setlg(e,iq);
   return merge_factor_i(decomp(d), g);
 }
 
-/* if y is non-NULL, it receives the discriminant
- * return basis if (ret_basis != 0), discriminant otherwise
- */
 static GEN
-nfbasis00(GEN x0, long flag, GEN p, long ret_basis, GEN *y)
+unscale_vecpol(GEN v, GEN h)
 {
-  GEN x, disc, basis, lead;
-  GEN *gptr[2];
-  long k, tetpil, av = avma, l = lgef(x0), smll;
+  long i, l;
+  GEN w;
+  if (!h) return v;
+  l = lg(v); w = cgetg(l, typ(v));
+  for (i=1; i<l; i++) w[i] = (long)unscale_pol((GEN)v[i], h);
+  return w;
+}
 
-  if (typ(x0)!=t_POL) err(typeer,"nfbasis00");
-  if (l<=3) err(zeropoler,"nfbasis00");
-  for (k=2; k<l; k++)
-    if (typ(x0[k])!=t_INT) err(talker,"polynomial not in Z[X] in nfbasis");
+/* FIXME: have to deal with compatibility flags */
+static void
+_nfbasis(GEN x0, long flag, GEN fa, GEN *pbas, GEN *pdK)
+{
+  GEN x, dx, dK, basis, lead, index;
+  long fl = 0;
 
-  x = pol_to_monic(x0,&lead);
+  if (typ(x0)!=t_POL) err(typeer,"nfbasis");
+  if (!degpol(x0)) err(zeropoler,"nfbasis");
+  check_pol_int(x0, "nfbasis");
 
-  if (!p || gcmp0(p))
-    smll = (flag & 1); /* small basis */
-  else
-  {
-    if (lead) p = update_fact(x,p);
-    smll = (long) p;   /* factored basis */
-  }
-
-  if (flag & 2)
-    basis = allbase(x,smll,&disc); /* round 2 */
-  else
-    basis = allbase4(x,smll,&disc,NULL); /* round 4 */
-
-  if (!ret_basis) return gerepilecopy(av,disc);
-
-  tetpil=avma;
-  if (!lead) basis = gcopy(basis);
-  else
-  {
-    long v = varn(x);
-    GEN pol = gmul(polx[v],lead);
-
-    tetpil = avma; basis = gsubst(basis,v,pol);
-  }
-  if (!y)
-    return gerepile(av,tetpil,basis);
-
-  *y = gcopy(disc);
-  gptr[0]=&basis; gptr[1]=y;
-  gerepilemanysp(av,tetpil,gptr,2);
-  return basis;
+  x = pol_to_monic(x0, &lead);
+  if (fa && gcmp0(fa)) fa = NULL; /* compatibility. NULL is the proper arg */
+  if (fa && lead) fa = update_fact(x, fa);
+  if (flag & compat_PARTIAL) fl |= nf_PARTIALFACT;
+  if (flag & compat_ROUND2)  fl |= nf_ROUND2;
+  basis = allbase(x, fl, &dx, &dK, &index, &fa);
+  if (pbas) *pbas = unscale_vecpol(basis, lead);
+  if (pdK)  *pdK = dK;
 }
 
 GEN
-nfbasis(GEN x, GEN *y, long flag, GEN p)
+nfbasis(GEN x, GEN *pdK, long flag, GEN fa)
 {
-  return nfbasis00(x,flag,p,1,y);
+  gpmem_t av = avma;
+  GEN bas; _nfbasis(x, flag, fa, &bas, pdK);
+  gerepileall(av, 2, &bas, pdK); return bas;
 }
 
 GEN
-nfbasis0(GEN x, long flag, GEN p)
+nfbasis0(GEN x, long flag, GEN fa)
 {
-  return nfbasis00(x,flag,p,1,NULL);
+  gpmem_t av = avma;
+  GEN bas; _nfbasis(x, flag, fa, &bas, NULL);
+  return gerepilecopy(av, bas);
 }
 
 GEN
-nfdiscf0(GEN x, long flag, GEN p)
+nfdiscf0(GEN x, long flag, GEN fa)
 {
-  return nfbasis00(x,flag,p,0,&p);
+  gpmem_t av = avma;
+  GEN dK; _nfbasis(x, flag, fa, NULL, &dK);
+  return gerepilecopy(av, dK);
 }
 
 GEN
-base(GEN x, GEN *y)
-{
-  return allbase4(x,0,y,NULL);
-}
+base(GEN x, GEN *pdK) { return nfbasis(x, pdK, 0, NULL); }
 
 GEN
-smallbase(GEN x, GEN *y)
-{
-  return allbase4(x,1,y,NULL);
-}
+smallbase(GEN x, GEN *pdK) { return nfbasis(x, pdK, compat_PARTIAL, NULL); }
 
 GEN
-factoredbase(GEN x, GEN p, GEN *y)
-{
-  return allbase4(x,(long)p,y,NULL);
-}
+factoredbase(GEN x, GEN fa, GEN *pdK) { return nfbasis(x, pdK, 0, fa); }
 
 GEN
-discf(GEN x)
-{
-  GEN y;
-  long av=avma,tetpil;
+discf(GEN x) { return nfdiscf0(x, 0, NULL); }
 
-  allbase4(x,0,&y,NULL); tetpil=avma;
-  return gerepile(av,tetpil,icopy(y));
-}
-
 GEN
-smalldiscf(GEN x)
-{
-  GEN y;
-  long av=avma,tetpil;
+smalldiscf(GEN x) { return nfdiscf0(x, nf_PARTIALFACT, NULL); }
 
-  allbase4(x,1,&y,NULL); tetpil=avma;
-  return gerepile(av,tetpil,icopy(y));
-}
-
 GEN
-factoreddiscf(GEN x, GEN p)
-{
-  GEN y;
-  long av=avma,tetpil;
+factoreddiscf(GEN x, GEN fa) { return nfdiscf0(x, 0, fa); }
 
-  allbase4(x,(long)p,&y,NULL); tetpil=avma;
-  return gerepile(av,tetpil,icopy(y));
-}
 
 /* return U if Z[alpha] is not maximal or 2*dU < m-1; else return NULL */
 static GEN
@@ -852,19 +848,16 @@ dedek(GEN f, long mf, GEN p,GEN g)
   GEN k,h;
   long dk;
 
-  if (DEBUGLEVEL>=3)
-  {
-    fprintferr("  entering dedek ");
-    if (DEBUGLEVEL>5)
-      fprintferr("with parameters p=%Z,\n  f=%Z",p,f);
-    fprintferr("\n");
-  }
   h = FpX_div(f,g,p);
   k = gdivexact(gadd(f, gneg_i(gmul(g,h))), p);
   k = FpX_gcd(k, FpX_gcd(g,h, p), p);
 
   dk = degpol(k);
-  if (DEBUGLEVEL>=3) fprintferr("  gcd has degree %ld\n", dk);
+  if (DEBUGLEVEL>2)
+  {
+    fprintferr("  dedek: gcd has degree %ld\n", dk);
+    if (DEBUGLEVEL>5) fprintferr("initial parameters p=%Z,\n  f=%Z\n",p,f);
+  }
   if (2*dk >= mf-1) return FpX_div(f,k,p);
   return dk? (GEN)NULL: f;
 }
@@ -873,7 +866,8 @@ dedek(GEN f, long mf, GEN p,GEN g)
 static GEN
 maxord(GEN p,GEN f,long mf)
 {
-  long j,r, av = avma, flw = (cmpsi(degpol(f),p) < 0);
+  const gpmem_t av = avma;
+  long j,r, flw = (cmpsi(degpol(f),p) < 0);
   GEN w,g,h,res;
 
   if (flw)
@@ -902,11 +896,8 @@ maxord(GEN p,GEN f,long mf)
 static GEN
 polmodiaux(GEN x, GEN y, GEN ys2)
 {
-  if (typ(x)!=t_INT)
-    x = mulii((GEN)x[1], mpinvmod((GEN)x[2],y));
-  x = modii(x,y); 
-  if (cmpii(x,ys2) > 0) x = subii(x,y);
-  return x;
+  if (typ(x)!=t_INT) x = mulii((GEN)x[1], mpinvmod((GEN)x[2],y));
+  return centermodii(x,y,ys2);
 }
 
 /* x polynomial with integer or rational coeff. Reduce them mod y IN PLACE */
@@ -931,79 +922,138 @@ polmodi_keep(GEN x, GEN y)
 }
 
 static GEN
+shiftpol(GEN x, long v)
+{
+  x = addshiftw(x, zeropol(v), 1);
+  setvarn(x,v); return x;
+}
+
+/* Sylvester's matrix, mod p^m (assumes f1 monic) */
+static GEN
+sylpm(GEN f1,GEN f2,GEN pm)
+{
+  long n,j,v=varn(f1);
+  GEN a,h;
+
+  n=degpol(f1); a=cgetg(n+1,t_MAT);
+  h = FpX_res(f2,f1,pm);
+  for (j=1;; j++)
+  {
+    a[j] = (long)pol_to_vec(h,n);
+    if (j == n) break;
+    h = FpX_res(shiftpol(h,v),f1,pm);
+  }
+  return hnfmodid(a,pm);
+}
+
+/* polynomial gcd mod p^m (assumes f1 monic) */
+GEN
+gcdpm(GEN f1,GEN f2,GEN pm)
+{
+  gpmem_t av=avma,tetpil;
+  long n,c,v=varn(f1);
+  GEN a,col;
+
+  n=degpol(f1); a=sylpm(f1,f2,pm);
+  for (c=1; c<=n; c++)
+    if (signe(resii(gcoeff(a,c,c),pm))) break;
+  if (c > n) { avma=av; return zeropol(v); }
+
+  col = gdiv((GEN)a[c], gcoeff(a,c,c)); tetpil=avma;
+  return gerepile(av,tetpil, gtopolyrev(col,v));
+}
+
+/* reduced resultant mod p^m (assumes x monic) */
+GEN
+respm(GEN x,GEN y,GEN pm)
+{
+  gpmem_t av = avma;
+  GEN p1 = sylpm(x,y,pm);
+
+  p1 = gcoeff(p1,1,1);
+  if (egalii(p1,pm)) { avma = av; return gzero; }
+  return gerepileuptoint(av, icopy(p1));
+}
+
+static GEN
 dbasis(GEN p, GEN f, long mf, GEN alpha, GEN U)
 {
-  long n=degpol(f),dU,c;
+  long n=degpol(f),dU,i;
   GEN b,ha,pd,pdp;
 
   if (n == 1) return gscalmat(gun, 1);
-  if (DEBUGLEVEL>=3)
+  if (DEBUGLEVEL>5)
   {
-    fprintferr("  entering Dedekind Basis ");
-    if (DEBUGLEVEL>5)
-    {
-      fprintferr("with parameters p=%Z\n",p);
-      fprintferr("  f = %Z,\n  alpha = %Z",f,alpha);
-    }
-    fprintferr("\n");
+    fprintferr("  entering Dedekind Basis with parameters p=%Z\n",p);
+    fprintferr("  f = %Z,\n  alpha = %Z\n",f,alpha);
   }
-  ha = pd = gpuigs(p,mf/2); pdp = mulii(pd,p);
+  ha = pd = gpowgs(p,mf/2); pdp = mulii(pd,p);
   dU = typ(U)==t_POL? degpol(U): 0;
   b = cgetg(n,t_MAT); /* Z[a] + U/p Z[a] is maximal */
   /* skip first column = gscalcol(pd,n) */
-  for (c=1; c<n; c++)
+  for (i=1; i<n; i++)
   {
-    if (c == dU)
+    if (i == dU)
     {
-      ha = gdiv(gmul(pd,eleval(f,U,alpha)),p);
+      ha = gdiv(gmul(pd,RX_RXQ_compo(U,alpha,f)),p);
       ha = polmodi(ha,pdp);
     }
     else
     {
-      GEN p2, mod;
+      GEN d, mod;
       ha = gmul(ha,alpha);
-      p2 = content(ha); /* to cancel denominator */
-      if (gcmp1(p2)) { p2 = NULL; mod = pdp; }
-      else
-      {
-        ha = gdiv(ha,p2);
-        if (typ(p2)==t_INT)
-          mod = divii(pdp, mppgcd(pdp,p2));
-        else
-          mod = mulii(pdp, (GEN)p2[2]); /* p2 = a / p^e */
-      }
+      ha = Q_remove_denom(ha, &d);
+      mod = d? mulii(pdp,d): pdp;
       ha = FpX_res(ha, f, mod);
-      if (p2) ha = gmul(ha,p2);
+      if (d) ha = gdivexact(ha,d);
     }
-    b[c] = (long)pol_to_vec(ha,n);
+    b[i] = (long)pol_to_vec(ha,n);
   }
   b = hnfmodid(b,pd);
   if (DEBUGLEVEL>5) fprintferr("  new order: %Z\n",b);
-  return gdiv(b,pd);
+  return gdiv(b, pd);
 }
 
 static GEN
 get_partial_order_as_pols(GEN p, GEN f)
 {
-  long i,j, n = degpol(f), vf = varn(f);
-  GEN b,ib,h,col;
+  GEN b = maxord(p,f, ggval(ZX_disc(f),p));
+  return mat_to_vecpol(b, varn(f));
+}
 
-  b = maxord(p,f, ggval(ZX_disc(f),p));
-  ib = cgetg(n+1,t_VEC);
-  for (i=1; i<=n; i++)
+long
+FpX_val(GEN x0, GEN t, GEN p, GEN *py)
+{
+  long k;
+  GEN r, y, x = x0;
+
+  for (k=0; ; k++)
   {
-    h=cgetg(i+2,t_POL); ib[i]=(long)h; col=(GEN)b[i];
-    h[1]=evalsigne(1)|evallgef(i+2)|evalvarn(vf);
-    for (j=1;j<=i;j++) h[j+1]=col[j];
+    y = FpX_divres(x, t, p, &r);
+    if (signe(r)) break;
+    x = y;
   }
-  return ib;
+  *py = x; return k;
 }
 
+/* e in Qp, f i Zp. Return r = e mod (f, pk) */
+static GEN
+QpX_mod(GEN e, GEN f, GEN pk)
+{
+  GEN mod, d;
+  e = Q_remove_denom(e, &d);
+  mod = d? mulii(pk,d): pk;
+  e = FpX_res(e, centermod(f, mod), mod);
+  e = FpX_center(e, mod);
+  if (d) e = gdiv(e, d);
+  return e;
+}
+
 /* if flag != 0, factorization to precision r (maximal order otherwise) */
 GEN
 Decomp(GEN p,GEN f,long mf,GEN theta,GEN chi,GEN nu,long flag)
 {
-  GEN res,pr,pk,ph,pdr,unmodp,b1,b2,b3,a1,e,f1,f2;
+  GEN fred,res,pr,pk,ph,pdr,b1,b2,a,e,f1,f2;
 
   if (DEBUGLEVEL>2)
   {
@@ -1016,77 +1066,65 @@ Decomp(GEN p,GEN f,long mf,GEN theta,GEN chi,GEN nu,lo
     }
     fprintferr("\n");
   }
-  unmodp = gmodulsg(1,p);
-  b1=lift_intern(gmul(chi,unmodp));
-  a1=gun; b2=gun;
-  b3=lift_intern(gmul(nu,unmodp));
-  while (degpol(b3) > 0)
-  {
-    GEN p1;
-    b1 = FpX_div(b1,b3, p);
-    b2 = FpX_red(gmul(b2,b3), p);
-    b3 = FpX_extgcd(b2,b1, p, &a1,&p1); /* p1 = junk */
-    p1 = leading_term(b3);
-    if (!gcmp1(p1))
-    { /* FpX_extgcd does not return normalized gcd */
-      p1 = mpinvmod(p1,p);
-      b3 = gmul(b3,p1);
-      a1 = gmul(a1,p1);
-    }
-  }
-  pdr = respm(f,derivpol(f),gpuigs(p,mf+1));
-  e = eleval(f,FpX_red(gmul(a1,b2), p),theta);
+
+  (void)FpX_val(chi, nu, p, &b1); /* nu irreducible mod p */
+  b2 = FpX_div(chi, b1, p);
+  a = FpX_mul(FpXQ_inv(b2, b1, p), b2, p);
+  pdr = respm(f, derivpol(f), gpowgs(p,mf+1));
+  e = RX_RXQ_compo(a, theta, f);
   e = gdiv(polmodi(gmul(pdr,e), mulii(pdr,p)),pdr);
 
   pr = flag? gpowgs(p,flag): mulii(p,sqri(pdr));
-  pk=p; ph=mulii(pdr,pr);
+  pk = p; ph = mulii(pdr,pr);
   /* E(t) - e(t) belongs to p^k Op, which is contained in p^(k-df)*Zp[xi] */
   while (cmpii(pk,ph) < 0)
-  {
+  { /* e <-- e^2(3-2e) mod p^2k */
+    pk = sqri(pk);
     e = gmul(gsqr(e), gsubsg(3,gmul2n(e,1)));
-    e = gres(e,f); pk = sqri(pk);
-    e = gdiv(polmodi(gmul(pdr,e), mulii(pdr,pk)), pdr);
+    e = QpX_mod(e, f, pk);
   }
-  f1 = gcdpm(f,gmul(pdr,gsubsg(1,e)), ph);
-  f1 = FpX_res(f1,f, pr);
-  f2 = FpX_res(FpX_div(f,f1, pr), f, pr);
+  fred = centermod(f, ph);
+  f1 = gcdpm(fred, gmul(pdr,gsubsg(1,e)), ph);
+  fred = centermod(fred, pr); /* pr | ph */
+  f1 = centermod(f1, pr);
+  f2 = FpX_div(fred,f1, pr);
+  f2 = FpX_center(f2, pr);
 
-  if (DEBUGLEVEL>2)
+  if (DEBUGLEVEL>5)
   {
     fprintferr("  leaving Decomp");
-    if (DEBUGLEVEL>5)
-      fprintferr(" with parameters: f1 = %Z\nf2 = %Z\ne = %Z\n", f1,f2,e);
-    fprintferr("\n");
+    fprintferr(" with parameters: f1 = %Z\nf2 = %Z\ne = %Z\n\n", f1,f2,e);
   }
 
   if (flag)
   {
-    b1=factorpadic4(f1,p,flag);
-    b2=factorpadic4(f2,p,flag); res=cgetg(3,t_MAT);
-    res[1]=lconcat((GEN)b1[1],(GEN)b2[1]);
-    res[2]=lconcat((GEN)b1[2],(GEN)b2[2]); return res;
+    b1 = factorpadic4(f1,p,flag);
+    b2 = factorpadic4(f2,p,flag); res = cgetg(3,t_MAT);
+    res[1] = (long)concatsp((GEN)b1[1],(GEN)b2[1]);
+    res[2] = (long)concatsp((GEN)b1[2],(GEN)b2[2]); return res;
   }
   else
   {
-    GEN ib1,ib2;
-    long n1,n2,i;
-    ib1 = get_partial_order_as_pols(p,f1); n1=lg(ib1)-1;
-    ib2 = get_partial_order_as_pols(p,f2); n2=lg(ib2)-1;
-    res=cgetg(n1+n2+1,t_VEC);
+    GEN ib1, ib2;
+    long n, n1, n2, i;
+    ib1 = get_partial_order_as_pols(p,f1); n1 = lg(ib1)-1;
+    ib2 = get_partial_order_as_pols(p,f2); n2 = lg(ib2)-1;
+    n = n1+n2;
+    res = cgetg(n+1, t_VEC);
     for (i=1; i<=n1; i++)
-      res[i]=(long)polmodi(gmod(gmul(gmul(pdr,(GEN)ib1[i]),e),f), pdr);
-    e=gsubsg(1,e); ib2 -= n1;
-    for (   ; i<=n1+n2; i++)
-      res[i]=(long)polmodi(gmod(gmul(gmul(pdr,(GEN)ib2[i]),e),f), pdr);
-    return nbasis(res,pdr);
+      res[i] = (long)QpX_mod(gmul(gmul(pdr,(GEN)ib1[i]),e), f, pdr);
+    e = gsubsg(1,e); ib2 -= n1;
+    for (   ; i<=n; i++)
+      res[i] = (long)QpX_mod(gmul(gmul(pdr,(GEN)ib2[i]),e), f, pdr);
+    res = vecpol_to_mat(res, n);
+    return gdiv(hnfmodid(res,pdr), pdr); /* normalized integral basis */
   }
 }
 
 /* minimum extension valuation: res[0]/res[1] (both are longs) */
-static long *
-vstar(GEN p,GEN h)
+static void
+vstar(GEN p,GEN h, long *L, long *E)
 {
-  static long res[2];
   long m,first,j,k,v,w;
 
   m=degpol(h); first=1; k=1; v=0;
@@ -1098,7 +1136,8 @@ vstar(GEN p,GEN h)
       first=0;
     }
   m = cgcd(v,k);
-  res[0]=v/m; res[1]=k/m; return res;
+  *L = v/m;
+  *E = k/m;
 }
 
 /* reduce the element elt modulo rd, taking first care of the denominators */
@@ -1107,7 +1146,7 @@ redelt(GEN elt, GEN rd, GEN pd)
 {
   GEN den, nelt, nrd, relt;
 
-  den  = ggcd(denom(content(elt)), pd);
+  den  = ggcd(Q_denom(elt), pd);
   nelt = gmul(den, elt);
   nrd  = gmul(den, rd);
 
@@ -1123,18 +1162,19 @@ redelt(GEN elt, GEN rd, GEN pd)
 GEN
 polsymmodpp(GEN g, GEN pp)
 {
-  long av1, av2, d = degpol(g), i, k;
+  gpmem_t av1, av2;
+  long d = degpol(g), i, k;
   GEN s , y;
 
-  y = cgetg(d + 1, t_COL); 
+  y = cgetg(d + 1, t_COL);
   y[1] = lstoi(d);
   for (k = 1; k < d; k++)
   {
-    av1 = avma; 
+    av1 = avma;
     s = centermod(gmulsg(k, polcoeff0(g,d-k,-1)), pp);
     for (i = 1; i < k; i++)
       s = gadd(s, gmul((GEN)y[k-i+1], polcoeff0(g,d-i,-1)));
-    av2 = avma; 
+    av2 = avma;
     y[k + 1] = lpile(av1, av2, centermod(gneg(s), pp));
   }
 
@@ -1152,26 +1192,27 @@ manage_cache(GEN chi, GEN pp, GEN ns)
   {
     if (DEBUGLEVEL > 4)
       fprintferr("newtonsums: result too large to fit in cache\n");
-    return polsymmodpp(chi, pp); 
+    return polsymmodpp(chi, pp);
   }
 
   if (!signe((GEN)ns[1]))
   {
-    ns2 = polsymmodpp(chi, pp); 
-    for (j = 1; j <= n; j++) 
+    ns2 = polsymmodpp(chi, pp);
+    for (j = 1; j <= n; j++)
       affii((GEN)ns2[j], (GEN)ns[j]);
   }
-  
+
   return ns;
 }
 
-/* compute the Newton sums modulo pp of the characteristic 
+/* compute the Newton sums modulo pp of the characteristic
    polynomial of a(x) mod g(x) */
 static GEN
 newtonsums(GEN a, GEN chi, GEN pp, GEN ns)
 {
   GEN va, pa, s, ns2;
-  long j, k, n = degpol(chi), av2, lim;
+  long j, k, n = degpol(chi);
+  gpmem_t av2, lim;
 
   ns2 = manage_cache(chi, pp, ns);
 
@@ -1184,17 +1225,17 @@ newtonsums(GEN a, GEN chi, GEN pp, GEN ns)
   for (j = 1; j <= n; j++)
   {
     pa = gmul(pa, a);
-    if (pp) pa = polmodi(pa, pp);
+    pa = polmodi(pa, pp);
     pa = gmod(pa, chi);
-    if (pp) pa = polmodi(pa, pp);
+    pa = polmodi(pa, pp);
 
     s  = gzero;
 
     for (k = 0; k <= n-1; k++)
       s = addii(s, mulii(polcoeff0(pa, k, -1), (GEN)ns2[k+1]));
-    
-    if (pp) va[j] = (long)centermod(s, pp);
-    
+
+    va[j] = (long)centermod(s, pp);
+
     if (low_stack(lim, stack_lim(av2, 1)))
     {
       GEN *gptr[2];
@@ -1213,7 +1254,8 @@ static GEN
 newtoncharpoly(GEN a, GEN chi, GEN pp, GEN ns)
 {
   GEN v, c, s, t;
-  long n = degpol(chi), j, k, vn = varn(chi), av = avma, av2, lim;
+  long n = degpol(chi), j, k, vn = varn(chi);
+  gpmem_t av = avma, av2, lim;
 
   v = newtonsums(a, chi, pp, ns);
   av2 = avma;
@@ -1240,67 +1282,61 @@ newtoncharpoly(GEN a, GEN chi, GEN pp, GEN ns)
       c = gerepilecopy(av2, c);
     }
   }
- 
+
   k = (n%2)? 1: 2;
-  for (  ; k <= n+1; k += 2) 
+  for (  ; k <= n+1; k += 2)
     c[k] = lneg((GEN)c[k]);
 
   return gerepileupto(av, gtopoly(c, vn));
 }
 
-static GEN 
+/* return v_p(n!) */
+static long
+val_fact(long n, GEN pp)
+{
+  long p, q, v;
+  if (is_bigint(pp)) return 0;
+  q = p = itos(pp); v = 0;
+  do { v += n/q; q *= p; } while (n >= q);
+  return v;
+}
+
+static GEN
 mycaract(GEN f, GEN beta, GEN p, GEN pp, GEN ns)
 {
-  GEN p1, p2, chi, chi2, npp;
-  long j, a, v = varn(f), n = degpol(f);
+  GEN p1, chi, npp;
+  long v = varn(f), n = degpol(f);
 
   if (gcmp0(beta)) return zeropol(v);
 
-  p1 = content(beta);
-  if (gcmp1(p1)) p1 = NULL; 
-  else beta = gdiv(beta, p1);
-
+  beta = Q_primitive_part(beta,&p1);
   if (!pp)
-    chi = caractducos(f, beta, v);
+    chi = ZX_caract(f, beta, v);
   else
   {
-    a = 0;
-    for (j = 1; j <= n; j++) /* compute the extra precision needed */
-      a += ggval(stoi(j), p);
-    npp = mulii(pp, gpowgs(p, a));
-    if (p1) npp = gmul(npp, gpowgs(denom(p1), n));
+    npp = mulii(pp, gpowgs(p, val_fact(n, p)));
+    if (p1) npp = mulii(npp, gpowgs(denom(p1), n));
 
     chi = newtoncharpoly(beta, f, npp, ns);
   }
 
-  if (p1)
-  {
-    chi2 = cgetg(n+3, t_POL);
-    chi2[1] = chi[1];
-    p2 = gun;
-    for (j = n+2; j >= 2; j--)
-    {
-      chi2[j] = lmul((GEN)chi[j], p2);
-      p2 = gmul(p2, p1);
-    }
-  }
-  else 
-    chi2 = chi;
-  
-  if (!pp) return chi2;
+  if (p1) chi = rescale_pol(chi,p1);
 
+  if (!pp) return chi;
+
   /* this can happen only if gamma is incorrect (not an integer) */
-  if (divise(denom(content(chi2)), p)) return NULL;
-  
-  return redelt(chi2, pp, pp);
+  if (divise(Q_denom(chi), p)) return NULL;
+
+  return redelt(chi, pp, pp);
 }
 
 /* Factor characteristic polynomial of beta mod p */
 static GEN
 factcp(GEN p, GEN f, GEN beta, GEN pp, GEN ns)
 {
-  long av,l;
+  gpmem_t av;
   GEN chi,nu, b = cgetg(4,t_VEC);
+  long l;
 
   chi=mycaract(f,beta,p,pp,ns);
   av=avma; nu=(GEN)factmod(chi,p)[1]; l=lg(nu)-1;
@@ -1314,39 +1350,33 @@ factcp(GEN p, GEN f, GEN beta, GEN pp, GEN ns)
 static GEN
 getprime(GEN p, GEN chi, GEN phi, GEN chip, GEN nup, long *Lp, long *Ep)
 {
-  long v = varn(chi), L, E, r, s; 
-  GEN chin, pip, pp, vn;
+  long v = varn(chi), L, E, r, s;
+  GEN chin, pip, pp;
 
   if (gegal(nup, polx[v]))
     chin = chip;
   else
     chin = mycaract(chip, nup, p, NULL, NULL);
-  
-  vn = vstar(p, chin);
-  L  = vn[0]; 
-  E  = vn[1];
-  
-  cbezout(L, -E, &r, &s);
 
-  if (r <= 0) 
-  { 
-    long q = (-r) / E; 
-    q++; 
-    r += q*E; 
-    s += q*L; 
+  vstar(p, chin, &L, &E);
+  (void)cbezout(L, -E, &r, &s);
+  if (r <= 0)
+  {
+    long q = 1 + ((-r) / E);
+    r += q*E;
+    s += q*L;
   }
 
-  pip = eleval(chi, nup, phi);
-  pip = lift_intern(gpuigs(gmodulcp(pip, chi), r));
-  pp  = gpuigs(p, s); 
+  pip = RX_RXQ_compo(nup, phi, chi);
+  pip = lift_intern(gpowgs(gmodulcp(pip, chi), r));
+  pp  = gpowgs(p, s);
 
-  *Lp = L; 
-  *Ep = E;
-  return gdiv(pip, pp);
+  *Lp = L;
+  *Ep = E; return gdiv(pip, pp);
 }
 
 static GEN
-update_alpha(GEN p, GEN fx, GEN alph, GEN chi, GEN pmr, GEN pmf, long mf, 
+update_alpha(GEN p, GEN fx, GEN alph, GEN chi, GEN pmr, GEN pmf, long mf,
 	     GEN ns)
 {
   long l, v = varn(fx);
@@ -1366,18 +1396,19 @@ update_alpha(GEN p, GEN fx, GEN alph, GEN chi, GEN pmr
   for (;;)
   {
     if (signe(pdr)) break;
-    if (!nalph) nalph = gadd(alph, gmul(p, polx[v])); 
-    /* nchi was too reduced at this point */
-    nchi = mycaract(fx, nalph, NULL, NULL, ns); 
-    pdr = respm(nchi, derivpol(nchi), pmf);
+    if (!nalph) nalph = gadd(alph, gmul(p, polx[v]));
+    /* nchi was too reduced at this point; try a larger precision */
+    pmf  = sqri(pmf);
+    nchi = mycaract(fx, nalph, p, pmf, ns);
+    pdr  = respm(nchi, derivpol(nchi), pmf);
     if (signe(pdr)) break;
     if (DEBUGLEVEL >= 6)
       fprintferr("  non separable polynomial in update_alpha!\n");
-    /* at this point, we assume that chi is not square-free */ 
+    /* at this point, we assume that chi is not square-free */
     nalph = gadd(nalph, gmul(p, polx[v]));
     w = factcp(p, fx, nalph, NULL, ns);
-    nchi = (GEN)w[1]; 
-    nnu  = (GEN)w[2]; 
+    nchi = (GEN)w[1];
+    nnu  = (GEN)w[2];
     l    = itos((GEN)w[3]);
     if (l > 1) return Decomp(p, fx, mf, nalph, nchi, nnu, 0);
     pdr = respm(nchi, derivpol(nchi), pmr);
@@ -1389,8 +1420,7 @@ update_alpha(GEN p, GEN fx, GEN alph, GEN chi, GEN pmr
   {
     npmr  = mulii(sqri(pdr), p);
     nchi  = polmodi(nchi, npmr);
-    if (!nalph) nalph = redelt(alph, npmr, pmf);
-    else nalph = redelt(nalph, npmr, pmf);
+    nalph = redelt(nalph? nalph: alph, npmr, pmf);
   }
 
   affii(gzero, (GEN)ns[1]); /* kill cache again (contains data for fx) */
@@ -1404,31 +1434,26 @@ update_alpha(GEN p, GEN fx, GEN alph, GEN chi, GEN pmr
   return rep;
 }
 
-extern GEN Fp_factor_irred(GEN P,GEN l, GEN Q);
-
-/* flag != 0 iff we're looking for the p-adic factorization, 
+/* flag != 0 iff we're looking for the p-adic factorization,
    in which case it is the p-adic precision we want */
 GEN
 nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
 {
-  long Fa, La, Ea, oE, Fg, eq, er, v = varn(fx), i, nv, Le, Ee, N, l, vn;
+  long L, E, Fa, La, Ea, oE, Fg, eq, er, v = varn(fx), i, nv, Le, Ee, N, l, vn;
   GEN p1, alph, chi, nu, w, phi, pmf, pdr, pmr, kapp, pie, chib, ns;
-  GEN gamm, chig, nug, delt, beta, eta, chie, nue, pia, vb, opa;
+  GEN gamm, chig, nug, delt, beta, eta, chie, nue, pia, opa;
 
-  if (DEBUGLEVEL >= 3)
+  if (DEBUGLEVEL>2)
   {
-    if (flag)
-      fprintferr("  entering Nilord2 (factorization)");
-    else
-      fprintferr("  entering Nilord2 (basis/discriminant)");
-    if (DEBUGLEVEL >= 5)
+    fprintferr("  entering Nilord");
+    if (DEBUGLEVEL>4)
     {
       fprintferr(" with parameters: p = %Z, expo = %ld\n", p, mf);
       fprintferr("  fx = %Z, gx = %Z", fx, gx);
     }
     fprintferr("\n");
   }
-  
+
   pmf = gpowgs(p, mf + 1);
   pdr = respm(fx, derivpol(fx), pmf);
   pmr = mulii(sqri(pdr), p);
@@ -1436,7 +1461,7 @@ nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
   chi = polmodi_keep(fx, pmr);
 
   alph = polx[v];
-  nu = gx; 
+  nu = gx;
   N  = degpol(fx);
   oE = 0;
   opa = NULL;
@@ -1452,7 +1477,7 @@ nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
 
   for(;;)
   {
-    /* kappa need to be recomputed */
+    /* kappa needs to be recomputed */
     kapp = NULL;
     Fa   = degpol(nu);
     /* the prime element in Zp[alpha] */
@@ -1471,16 +1496,16 @@ nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
       pia  = getprime(p, chi, polx[v], chi, nu, &La, &Ea);
       pia  = redelt(pia, pmr, pmf);
     }
-    
-    oE = Ea; opa = eleval(fx, pia, alph);
 
+    oE = Ea; opa = RX_RXQ_compo(pia, alph, fx);
+
     if (DEBUGLEVEL >= 5)
       fprintferr("  Fa = %ld and Ea = %ld \n", Fa, Ea);
-   
+
     /* we change alpha such that nu = pia */
     if (La > 1)
     {
-      alph = gadd(alph, eleval(fx, pia, alph));
+      alph = gadd(alph, RX_RXQ_compo(pia, alph, fx));
 
       w = update_alpha(p, fx, alph, NULL, pmr, pmf, mf, ns);
       alph = (GEN)w[1];
@@ -1490,7 +1515,7 @@ nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
     }
 
     /* if Ea*Fa == N then O = Zp[alpha] */
-    if (Ea*Fa == N) 
+    if (Ea*Fa == N)
     {
       if (flag) return NULL;
 
@@ -1515,15 +1540,15 @@ nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
 	eq = (long)(vn / N);
 	er = (long)(vn*Ea/N - eq*Ea);
       }
-      else  
+      else
       {
 	chib = mycaract(chi, beta, NULL, NULL, ns);
-	vb = vstar(p, chib);
-	eq = (long)(vb[0] / vb[1]);
-	er = (long)(vb[0]*Ea / vb[1] - eq*Ea);
+	vstar(p, chib, &L, &E);
+	eq = (long)(L / E);
+	er = (long)(L*Ea / E - eq*Ea);
       }
-      
-      /* eq and er are such that gamma = beta.p^-eq.nu^-er is a unit */ 
+
+      /* eq and er are such that gamma = beta.p^-eq.nu^-er is a unit */
       if (eq) gamm = gdiv(beta, gpowgs(p, eq));
       else gamm = beta;
 
@@ -1532,7 +1557,7 @@ nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
 	/* kappa = nu^-1 in Zp[alpha] */
 	if (!kapp)
 	{
-	  kapp = ginvmod(nu, chi);
+	  kapp = QX_invmod(nu, chi);
 	  kapp = redelt(kapp, pmr, pmf);
 	  kapp = gmodulcp(kapp, chi);
 	}
@@ -1544,7 +1569,7 @@ nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
 	fprintferr("  gamma = %Z\n", gamm);
 
       if (er || !chib)
-	chig = mycaract(chi, gamm, p, pmf, ns); 
+	chig = mycaract(chi, gamm, p, pmf, ns);
       else
       {
 	chig = poleval(chib, gmul(polx[v], gpowgs(p, eq)));
@@ -1552,23 +1577,24 @@ nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
 	chig = polmodi(chig, pmf);
       }
 
-      if (!chig || !gcmp1(denom(content(chig))))
+      if (!chig || !gcmp1(Q_denom(chig)))
       {
-	/* the valuation of beta was wrong... This also means 
-	   that chi_gamma has more than one factor modulo p   */
-	if (!chig) chig = mycaract(chi, gamm, p, NULL, NULL);
+	/* Valuation of beta was wrong. This means that
+	   gamma fails the v*-test */
+	chib = mycaract(chi, beta, p, NULL, ns);
+	vstar(p, chib, &L, &E);
+	eq = (long)(-L / E);
+	er = (long)(-L*Ea / E - eq*Ea);
 
-	vb = vstar(p, chig);
-	eq = (long)(-vb[0] / vb[1]);
-	er = (long)(-vb[0]*Ea / vb[1] - eq*Ea);
-	if (eq) gamm = gmul(gamm, gpowgs(p, eq));
-	if (er) 
+	if (eq) gamm = gmul(beta, gpowgs(p, eq));
+	else gamm = beta;
+	if (er)
 	{
 	  gamm = gmul(gamm, gpowgs(nu, er));
 	  gamm = gmod(gamm, chi);
 	  gamm = redelt(gamm, p, pmr);
 	}
-	if (eq || er) chig = mycaract(chi, gamm, p, pmf, ns); 
+	chig = mycaract(chi, gamm, p, pmf, ns);
       }
 
       nug  = (GEN)factmod(chig, p)[1];
@@ -1578,7 +1604,7 @@ nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
       if (l > 1)
       {
 	/* there are at least 2 factors mod. p => chi can be split */
-	phi  = eleval(fx, gamm, alph);
+	phi  = RX_RXQ_compo(gamm, alph, fx);
 	phi  = redelt(phi, p, pmf);
 	if (flag) mf += 3;
         return Decomp(p, fx, mf, phi, chig, nug, flag);
@@ -1594,7 +1620,7 @@ nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
 	if (gcmp1((GEN)w[1]))
 	{
 	  /* there are at least 2 factors mod. p => chi can be split */
-	  phi = eleval(fx, (GEN)w[2], alph);
+	  phi = RX_RXQ_compo((GEN)w[2], alph, fx);
 	  phi = redelt(phi, p, pmf);
           if (flag) mf += 3;
           return Decomp(p, fx, mf, phi, (GEN)w[3], (GEN)w[4], flag);
@@ -1613,10 +1639,10 @@ nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
       for (i = 1;; i++)
       {
 	if (i >= lg(w))
-          err(talker, "bug in nilord (no suitable root), is p = %Z a prime?", 
+          err(talker, "bug in nilord (no suitable root), is p = %Z a prime?",
 	      (long)p);
         delt = gneg_i(gsubst(gcoeff(w, 2, i), nv, polx[v]));
-        eta  = gsub(gamm, delt);	  
+        eta  = gsub(gamm, delt);	
         if (typ(delt) == t_INT)
         {
           chie = poleval(chig, gadd(polx[v], delt));
@@ -1625,28 +1651,28 @@ nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
           nue  = lift((GEN)nue[l]);
         }
         else
-        { 
+        {
           p1   = factcp(p, chi, eta, pmf, ns);
           chie = (GEN)p1[1];
-          nue  = (GEN)p1[2]; 
+          nue  = (GEN)p1[2];
           l    = itos((GEN)p1[3]);
         }
         if (l > 1)
         {
           /* there are at least 2 factors mod. p => chi can be split */
-          delete_var();      
-          phi = eleval(fx, eta, alph);
+          (void)delete_var();
+          phi = RX_RXQ_compo(eta, alph, fx);
           phi = redelt(phi, p, pmf);
           if (flag) mf += 3;
           return Decomp(p, fx, mf, phi, chie, nue, flag);
         }
-        
+
         /* if vp(eta) = vp(gamma - delta) > 0 */
         if (gegal(nue, polx[v])) break;
       }
-      delete_var();    
+      (void)delete_var();
 
-      if (!signe(modii((GEN)chie[2], pmr))) 
+      if (!signe(modii(constant_term(chie), pmr)))
 	chie = mycaract(chi, eta, p, pmf, ns);
 	
       pie = getprime(p, chi, eta, chie, nue, &Le, &Ee);
@@ -1660,21 +1686,21 @@ nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
 	if (gcmp1((GEN)w[1]))
 	{
 	  /* there are at least 2 factors mod. p => chi can be split */
-	  phi = eleval(fx, (GEN)w[2], alph);
+	  phi = RX_RXQ_compo((GEN)w[2], alph, fx);
 	  phi = redelt(phi, p, pmf);
           if (flag) mf += 3;
           return Decomp(p, fx, mf, phi, (GEN)w[3], (GEN)w[4], flag);
 	}
 	break;
       }
-      
+
       if (eq) delt = gmul(delt, gpowgs(p, eq));
       if (er) delt = gmul(delt, gpowgs(nu, er));
       beta = gsub(beta, delt);
     }
 
     /* we replace alpha by a new alpha with a larger F or E */
-    alph = eleval(fx, (GEN)w[2], alph);
+    alph = RX_RXQ_compo((GEN)w[2], alph, fx);
     chi  = (GEN)w[3];
     nu   = (GEN)w[4];
 
@@ -1702,8 +1728,8 @@ nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
 }
 
 /* Returns [1,phi,chi,nu] if phi non-primary
- *         [2,phi,chi,nu] with F_phi = lcm (F_alpha, F_theta) 
- *                         and E_phi = E_alpha 
+ *         [2,phi,chi,nu] with F_phi = lcm (F_alpha, F_theta)
+ *                         and E_phi = E_alpha
  */
 static GEN
 testb2(GEN p, GEN fa, long Fa, GEN theta, GEN pmf, long Ft, GEN ns)
@@ -1711,9 +1737,9 @@ testb2(GEN p, GEN fa, long Fa, GEN theta, GEN pmf, lon
   long m, Dat, t, v = varn(fa);
   GEN b, w, phi, h;
 
-  Dat = clcm(Fa, Ft) + 3; 
+  Dat = clcm(Fa, Ft) + 3;
   b = cgetg(5, t_VEC);
-  m = p[2]; 
+  m = p[2];
   if (degpol(p) > 0 || m < 0) m = 0;
 
   for (t = 1;; t++)
@@ -1725,10 +1751,10 @@ testb2(GEN p, GEN fa, long Fa, GEN theta, GEN pmf, lon
     if (h[2] > 1) { b[1] = un; break; }
     if (lgef(w[2]) == Dat) { b[1] = deux; break; }
   }
-  
-  b[2] = (long)phi; 
-  b[3] = w[1]; 
-  b[4] = w[2]; 
+
+  b[2] = (long)phi;
+  b[3] = w[1];
+  b[4] = w[2];
   return b;
 }
 
@@ -1736,7 +1762,7 @@ testb2(GEN p, GEN fa, long Fa, GEN theta, GEN pmf, lon
  *         [2, phi, chi, nu] if E_phi = lcm (E_alpha, E_theta)
  */
 static GEN
-testc2(GEN p, GEN fa, GEN pmr, GEN pmf, GEN alph2, long Ea, GEN thet2, 
+testc2(GEN p, GEN fa, GEN pmr, GEN pmf, GEN alph2, long Ea, GEN thet2,
        long Et, GEN ns)
 {
   GEN b, c1, c2, c3, psi, phi, w, h;
@@ -1744,13 +1770,13 @@ testc2(GEN p, GEN fa, GEN pmr, GEN pmf, GEN alph2, lon
 
   b=cgetg(5, t_VEC);
 
-  cbezout(Ea, Et, &r, &s); t = 0;
+  (void)cbezout(Ea, Et, &r, &s); t = 0;
   while (r < 0) { r = r + Et; t++; }
   while (s < 0) { s = s + Ea; t++; }
 
-  c1 = lift_intern(gpuigs(gmodulcp(alph2, fa), s));
-  c2 = lift_intern(gpuigs(gmodulcp(thet2, fa), r));
-  c3 = gdiv(gmod(gmul(c1, c2), fa), gpuigs(p, t));
+  c1 = lift_intern(gpowgs(gmodulcp(alph2, fa), s));
+  c2 = lift_intern(gpowgs(gmodulcp(thet2, fa), r));
+  c3 = gdiv(gmod(gmul(c1, c2), fa), gpowgs(p, t));
 
   psi = redelt(c3, pmr, pmr);
   phi = gadd(polx[v], psi);
@@ -1758,549 +1784,815 @@ testc2(GEN p, GEN fa, GEN pmr, GEN pmf, GEN alph2, lon
   w = factcp(p, fa, phi, pmf, ns);
   h = (GEN)w[3];
   b[1] = (h[2] > 1)? un: deux;
-  b[2] = (long)phi; 
-  b[3] = w[1]; 
-  b[4] = w[2]; 
+  b[2] = (long)phi;
+  b[3] = w[1];
+  b[4] = w[2];
   return b;
 }
 
-/* evaluate h(a) mod f */
-GEN
-eleval(GEN f,GEN h,GEN a)
+/*******************************************************************/
+/*                                                                 */
+/*    2-ELT REPRESENTATION FOR PRIME IDEALS (dividing index)       */
+/*                                                                 */
+/*******************************************************************/
+/* beta = generators of prime ideal (vectors). Return "small" random elt in P */
+static GEN
+random_elt_in_P(GEN beta, long small)
 {
-  long n,av,tetpil;
-  GEN y;
+  long z, i, j, la, lbeta = lg(beta);
+  GEN a = NULL, B;
 
-  if (typ(h) != t_POL) return gcopy(h);
-  av = tetpil = avma;
-  n=lgef(h)-1; y=(GEN)h[n];
-  for (n--; n>=2; n--)
+  /* compute sum random(-7..8) * beta[i] */
+  if (small)
   {
-    y = gadd(gmul(y,a),(GEN)h[n]);
-    tetpil=avma; y = gmod(y,f);
+    la = lg(beta[1]);
+    if (small > 7) small = 0;
+    for (i=1; i<lbeta; i++)
+    { /* Warning: change definition of 'small' if you change this */
+      z = mymyrand() >> (BITS_IN_RANDOM-5); /* in [0,15] */
+      if (z >= 9) z -= 7;
+      if (small) z %= small;
+      if (!z) continue;
+      B = (GEN)beta[i];
+      if (a)
+        for (j = 1; j < la; j++) a[j] += z * B[j];
+      else {
+        a = cgetg(la, t_VECSMALL);
+        for (j = 1; j <la; j++) a[j] = z * B[j];
+      }
+    }
   }
-  return gerepile(av,tetpil,y);
+  else
+    for (i=1; i<lbeta; i++)
+    {
+      z = mymyrand() >> (BITS_IN_RANDOM-5);
+      if (z >= 9) z -= 7;
+      if (!z) continue;
+      B = gmulsg(z, (GEN)beta[i]);
+      a = a? gadd(a, B): B;
+    }
+  return a;
 }
 
-GEN addshiftw(GEN x, GEN y, long d);
+/* routines to check uniformizer algebraically (as t_POL) */
 
-static GEN
-shiftpol(GEN x, long v)
+/* a t_POL, D t_INT (or NULL if 1), q = p^(f+1). a/D in pr | p, norm(pr) = pf.
+ * Return 1 if (a/D,p) = pr, and 0 otherwise */
+static int
+is_uniformizer(GEN D, GEN a, GEN T, GEN q)
 {
-  x = addshiftw(x, zeropol(v), 1);
-  setvarn(x,v); return x;
+  GEN N = ZX_resultant_all(T, a, D, 0); /* norm(a) */
+  return (resii(N, q) != gzero);
 }
 
-/* Sylvester's matrix, mod p^m (assumes f1 monic) */
+/* a/D or a/D + p uniformizer ? */
 static GEN
-sylpm(GEN f1,GEN f2,GEN pm)
+prime_check_elt(GEN D, GEN Dp, GEN a, GEN T, GEN q, int ramif)
 {
-  long n,j,v=varn(f1);
-  GEN a,h;
+  if (is_uniformizer(D,a,T,q)) return a;
+  /* can't succeed if e > 1 */
+  if (!ramif && is_uniformizer(D,gadd(a,Dp),T,q)) return a;
+  return NULL;
+}
 
-  n=degpol(f1); a=cgetg(n+1,t_MAT);
-  h = FpX_res(f2,f1,pm);
-  for (j=1;; j++)
+/* P given by an Fp basis */
+static GEN
+compute_pr2(GEN nf, GEN P, GEN p, int *ramif)
+{
+  long i, j, m = lg(P)-1;
+  GEN mul, P2 = gmul(p, P);
+  for (i=1; i<=m; i++)
   {
-    a[j] = (long)pol_to_vec(h,n);
-    if (j == n) break;
-    h = FpX_res(shiftpol(h,v),f1,pm);
+    mul = eltmul_get_table(nf, (GEN)P[i]);
+    for (j=1; j<=i; j++)
+      P2 = concatsp(P2, gmul(mul, (GEN)P[j]));
   }
-  return hnfmodid(a,pm);
+  P2 = hnfmodid(P2, sqri(p)); /* HNF of pr^2 */
+  *ramif = egalii(gcoeff(P2,1,1), p); /* pr/p ramified ? */
+  return P2;
 }
 
-/* polynomial gcd mod p^m (assumes f1 monic) */
-GEN
-gcdpm(GEN f1,GEN f2,GEN pm)
+static GEN
+random_unif_loop_pol(GEN nf, GEN P, GEN D, GEN Dp, GEN beta, GEN pol,
+                     GEN p, GEN q)
 {
-  long n,c,v=varn(f1),av=avma,tetpil;
-  GEN a,col;
+  long small, i, c = 0, m = lg(P)-1, N = degpol(nf[1]), keep = getrand();
+  int ramif;
+  GEN a, P2;
+  gpmem_t av;
 
-  n=degpol(f1); a=sylpm(f1,f2,pm);
-  for (c=1; c<=n; c++)
-    if (signe(resii(gcoeff(a,c,c),pm))) break;
-  if (c > n) { avma=av; return zeropol(v); }
+  for(i=1; i<=m; i++)
+    if ((a = prime_check_elt(D,Dp,(GEN)beta[i],pol,q,0))) return a;
 
-  col = gdiv((GEN)a[c], gcoeff(a,c,c)); tetpil=avma;
-  return gerepile(av,tetpil, gtopolyrev(col,v));
+  (void)setrand(1);
+  if (DEBUGLEVEL) fprintferr("uniformizer_loop, hard case: ");
+  P2 = compute_pr2(nf, P, p, &ramif);
+
+  a = mulis(Dp, 8*degpol(nf[1]));
+  if (is_bigint(a)) small = 0;
+  else
+  {
+    ulong mod = itou(Dp);
+    small = itos(p);
+    for (i=1; i<=m; i++)
+      beta[i] = (long)u_Fp_FpV(pol_to_vec((GEN)beta[i], N), mod);
+  }
+
+  for(av = avma; ;avma = av)
+  {
+    if (DEBUGLEVEL && (++c & 0x3f) == 0) fprintferr("%d ", c);
+    a = random_elt_in_P(beta, small);
+    if (!a) continue;
+    if (small) a = vec_to_pol(small_to_vec(a), varn(pol));
+    a = centermod(a, Dp);
+    if ((a = prime_check_elt(D,Dp,a,pol,q,ramif)))
+    {
+      if (DEBUGLEVEL) fprintferr("\n");
+      (void)setrand(keep); return a;
+    }
+  }
 }
 
-/* reduced resultant mod p^m (assumes x monic) */
-GEN
-respm(GEN x,GEN y,GEN pm)
-{
-  long av = avma;
-  GEN p1 = sylpm(x,y,pm);
+/* routines to check uniformizer numerically (as t_COL) */
 
-  p1 = gcoeff(p1,1,1);
-  if (egalii(p1,pm)) { avma = av; return gzero; }
-  return gerepileuptoint(av, icopy(p1));
+static int
+vec_is_uniformizer(GEN a, GEN q, long r1)
+{
+  long e;
+  GEN N = grndtoi( norm_by_embed(r1, a), &e );
+  if (e > -5) err(precer, "vec_is_uniformizer");
+  return (resii(N, q) != gzero);
 }
 
-/* Normalized integral basis */
 static GEN
-nbasis(GEN ibas,GEN pd)
+prime_check_eltvec(GEN a, GEN M, GEN p, GEN q, long r1, long small, int ramif)
 {
-  long k, n = lg(ibas)-1;
-  GEN a = cgetg(n+1,t_MAT);
-  for (k=1; k<=n; k++) a[k] = (long)pol_to_vec((GEN)ibas[k],n);
-  return gdiv(hnfmodid(a,pd), pd);
+  GEN ar;
+  long i,l;
+  a = centermod(a, p);
+  if (small) ar = gmul_mat_smallvec(M, a);
+  else       ar = gmul(M, a);
+
+  if (vec_is_uniformizer(ar, q, r1)) return small? small_to_col(a): a;
+  /* can't succeed if e > 1 */
+  if (ramif) return NULL;
+  l = lg(ar);
+  for (i=1; i<l; i++) ar[i] = ladd((GEN)ar[i], p);
+  if (!vec_is_uniformizer(ar, q, r1)) return NULL;
+  if (small) a = small_to_col(a);
+  a[1] = laddii((GEN)a[1],p); return a;
 }
 
-/*******************************************************************/
-/*                                                                 */
-/*                   BUCHMANN-LENSTRA ALGORITHM                    */
-/*                                                                 */
-/*******************************************************************/
-static GEN lens(GEN nf,GEN p,GEN a);
-GEN element_powid_mod_p(GEN nf, long I, GEN n, GEN p);
-
-/* return a Z basis of Z_K's p-radical, modfrob = x--> x^p-x */
 static GEN
-pradical(GEN nf, GEN p, GEN *modfrob)
+random_unif_loop_vec(GEN nf, GEN P, GEN p, GEN q)
 {
-  long i,N = degpol(nf[1]);
-  GEN p1,m,frob,rad;
+  long small, r1, i, c = 0, m = lg(P)-1, keep = getrand();
+  int ramif;
+  GEN a, P2, beta, M = gmael(nf,5,1);
+  gpmem_t av;
 
-  frob = cgetg(N+1,t_MAT);
-  for (i=1; i<=N; i++)
-    frob[i] = (long) element_powid_mod_p(nf,i,p, p);
-
-  /* p1 = smallest power of p st p^k >= N */
-  p1=p; while (cmpis(p1,N)<0) p1=mulii(p1,p);
-  if (p1==p) m = frob;
+  r1 = nf_get_r1(nf);
+  a = mulis(p, 8*degpol(nf[1]));
+  if (is_bigint(a)) { small = 0; beta = P; }
   else
   {
-    m=cgetg(N+1,t_MAT); p1 = divii(p1,p);
-    for (i=1; i<=N; i++)
-      m[i]=(long)element_pow_mod_p(nf,(GEN)frob[i],p1, p);
+    small = itos(p); beta = cgetg(m+1, t_MAT);
+    for (i=1; i<=m; i++)
+      beta[i] = (long)u_Fp_FpV((GEN)P[i], itou(p));
   }
-  rad = FpM_ker(m, p);
-  for (i=1; i<=N; i++)
-    coeff(frob,i,i) = lsubis(gcoeff(frob,i,i), 1);
-  *modfrob = frob; return rad;
-}
+  for(i=1; i<=m; i++)
+    if ((a = prime_check_eltvec((GEN)beta[i],M,p,q,r1,small,0))) return a;
 
-static GEN
-project(GEN algebre, GEN x, long k, long kbar)
-{
-  x = inverseimage(algebre,x);
-  x += k; x[0] = evaltyp(t_COL) | evallg(kbar+1);
-  return x;
-}
+  (void)setrand(1);
+  if (DEBUGLEVEL) fprintferr("uniformizer_loop, hard case: ");
+  P2 = compute_pr2(nf, P, p, &ramif);
 
-/* Calcule le polynome minimal de alpha dans algebre (coeffs dans Z) */
-static GEN
-pol_min(GEN alpha,GEN nf,GEN algebre,long kbar,GEN p)
-{
-  long av=avma,tetpil,i,N,k;
-  GEN p1,puiss;
-
-  N = lg(nf[1])-3; puiss=cgetg(N+2,t_MAT);
-  k = N-kbar; p1=alpha;
-  puiss[1] = (long)gscalcol_i(gun,kbar);
-  for (i=2; i<=N+1; i++)
+  for(av = avma; ;avma = av)
   {
-    if (i>2) p1 = element_mul(nf,p1,alpha);
-    puiss[i] = (long) project(algebre,p1,k,kbar);
+    if (DEBUGLEVEL && (++c & 0x3f) == 0) fprintferr("%d ", c);
+    a = random_elt_in_P(beta, small);
+    if (!a) continue;
+    if ((a = prime_check_eltvec(a,M,p,q,r1,small,0)))
+    {
+      if (DEBUGLEVEL) fprintferr("\n");
+      (void)setrand(keep); return a;
+    }
   }
-  puiss = lift_intern(puiss);
-  p1 = (GEN)FpM_ker(puiss, p)[1]; tetpil=avma;
-  return gerepile(av,tetpil,gtopolyrev(p1,0));
 }
 
-/* Evalue le polynome pol en alpha,element de nf */
+/* L = Fp-bases for prime ideals of O_K/p. Return a p-uniformizer for
+ * L[i] using the approximation theorem */
 static GEN
-eval_pol(GEN nf,GEN pol,GEN alpha,GEN algebre,GEN algebre1)
+uniformizer_appr(GEN nf, GEN L, long i, GEN p)
 {
-  long av=avma,tetpil,i,kbar,k, lx = lgef(pol)-1, N = degpol(nf[1]);
-  GEN res;
+  long j, l;
+  GEN inter = NULL, P = (GEN)L[i], P2, Q, u, v, pi;
+  int ramif;
 
-  kbar = lg(algebre1)-1; k = N-kbar;
-  res = gscalcol_i((GEN)pol[lx], N);
-  for (i=2; i<lx; i++)
+  l = lg(L)-1;
+  for (j=l; j; j--)
   {
-    res = element_mul(nf,alpha,res);
-    res[1] = ladd((GEN)res[1],(GEN)pol[i]);
+    if (j == i) continue;
+    Q = (GEN)L[j]; if (typ(Q) != t_MAT) break;
+    inter = inter? FpM_intersect(inter, Q, p): Q;
   }
-  res = project(algebre,res,k,kbar); tetpil=avma;
-  return gerepile(av,tetpil,gmul(algebre1,res));
+  if (!inter)
+  {
+    setlg(L, l);
+    inter = idealprodprime(nf, L);
+  }
+  else
+  {
+    inter = hnfmodid(inter, p);
+    for (   ; j; j--)
+    {
+      Q = (GEN)L[j];
+      inter = idealmul(nf, inter, Q);
+    }
+  }
+
+  /* inter = prod L[j], j != i */
+  P2 = compute_pr2(nf, P, p, &ramif);
+  u = addone_aux2(P2, inter);
+  v = unnf_minus_x(u);
+  if (!ramif) pi = gmul(p, v);
+  else
+  {
+    l = lg(P)-1;
+    for (j=l; j; j--)
+      if (!hnf_invimage(P2, (GEN)P[j])) break;
+    pi = element_muli(nf,(GEN)P[j],v);
+  }
+  pi = gadd(pi, u);
+  pi =  centermod(pi, p);
+  if (!ramif && hnf_invimage(P2, pi)) pi[1] = laddii((GEN)pi[1], p);
+  return pi;
 }
 
+/* Input: an ideal mod p, P != Z_K (Fp-basis, in matrix form)
+ * Output: x such that P = (p,x) */
 static GEN
-kerlens2(GEN x, GEN p)
+uniformizer(GEN nf, GEN P, GEN p)
 {
-  long i,j,k,t,nbc,nbl,av,av1;
-  GEN a,c,l,d,y,q;
+  GEN q, D, Dp, w, beta, a, T = (GEN)nf[1];
+  long i, f, N = degpol(T), m = lg(P)-1;
 
-  av=avma; a=gmul(x,gmodulsg(1,p));
-  nbl=nbc=lg(x)-1;
-  c=new_chunk(nbl+1); for (i=1; i<=nbl; i++) c[i]=0;
-  l=new_chunk(nbc+1);
-  d=new_chunk(nbc+1);
-  k = t = 1;
-  while (t<=nbl && k<=nbc)
+  if (!m) return gscalcol_i(p,N);
+
+  /* we want v_p(Norm(beta)) = p^f, f = N-m */
+  f = N-m;
+  P = centermod(P, p);
+  q = mulii(gpowgs(p,f), p);
+
+  if (typ(nf[5]) == t_VEC) /* dummy nf from padicff */
   {
-    for (j=1; j<k; j++)
-      for (i=1; i<=nbl; i++)
-	if (i!=l[j])
-	  coeff(a,i,k) = lsub(gmul((GEN)d[j],gcoeff(a,i,k)),
-	                      gmul(gcoeff(a,l[j],k),gcoeff(a,i,j)));
-    t=1; while (t<=nbl && (c[t] || gcmp0(gcoeff(a,t,k)))) t++;
-    if (t<=nbl) { d[k]=coeff(a,t,k); c[t]=k; l[k]=t; k++; }
+    GEN M = gmael(nf,5,1);
+    long ex, prec = gprecision(M);
+    ex = gexpo(M) + gexpo(mulsi(8 * N, p));
+    /* enough prec to use norm_by_embed */
+    if (N * ex <= bit_accuracy(prec))
+      return random_unif_loop_vec(nf, P, p, q);
   }
-  if (k>nbc) err(bugparier,"kerlens2");
-  y=cgetg(nbc+1,t_COL);
-  y[1]=(k>1)?coeff(a,l[1],k):un;
-  for (q=gun,j=2; j<k; j++)
+
+  w = Q_remove_denom((GEN)nf[7], &D);
+  if (D)
   {
-    q=gmul(q,(GEN)d[j-1]);
-    y[j]=lmul(gcoeff(a,l[j],k),q);
+    long v = pvaluation(D, p, &a);
+    D = gpowgs(p, v);
+    Dp = mulii(D, p);
+  } else {
+    w = dummycopy(w);
+    Dp = p;
   }
-  if (k>1) y[k]=lneg(gmul(q,(GEN)d[k-1]));
-  for (j=k+1; j<=nbc; j++) y[j]=zero;
-  av1=avma; return gerepile(av,av1,lift(y));
+  for (i=1; i<=N; i++)
+    w[i] = (long)FpX_red((GEN)w[i], Dp); /* w[i] / D defined mod p */
+  beta = gmul(w, P);
+  a = random_unif_loop_pol(nf,P,D,Dp,beta,T,p,q);
+  a = algtobasis_i(nf,a);
+  if (D) a = gdivexact(a,D);
+  a = centermod(a, p);
+  if (!is_uniformizer(D,gmul(w,a), T,q)) a[1] = laddii((GEN)a[1],p);
+  return a;
 }
 
+/* Assuming P = (p,u) prime, return tau such that p Z + tau Z = p P^(-1)*/
 static GEN
-kerlens(GEN x, GEN pgen)
+anti_uniformizer(GEN nf, GEN p, GEN u)
 {
-  long av = avma, i,j,k,t,nbc,nbl,p,q,*c,*l,*d,**a;
-  GEN y;
+  gpmem_t av = avma;
+  GEN mat = eltmul_get_table(nf, u);
+  return gerepileupto(av, FpM_deplin(mat,p));
+}
 
-  if (cmpis(pgen, MAXHALFULONG>>1) > 0)
-    return kerlens2(x,pgen);
-  /* ici p <= (MAXHALFULONG>>1) ==> long du C */
-  p=itos(pgen); nbl=nbc=lg(x)-1;
-  a=(long**)new_chunk(nbc+1);
-  for (j=1; j<=nbc; j++)
+/*******************************************************************/
+/*                                                                 */
+/*                   BUCHMANN-LENSTRA ALGORITHM                    */
+/*                                                                 */
+/*******************************************************************/
+
+/* pr = (p,u) of ramification index e */
+GEN
+apply_kummer(GEN nf,GEN u,GEN e,GEN p)
+{
+  GEN t, T = (GEN)nf[1], pr = cgetg(6,t_VEC);
+  long f = degpol(u), N = degpol(T);
+
+  pr[1] = (long)p;
+  pr[3] = (long)e;
+  pr[4] = lstoi(f);
+  if (f == N) /* inert */
   {
-    c=a[j]=new_chunk(nbl+1);
-    for (i=1; i<=nbl; i++) c[i]=smodis(gcoeff(x,i,j),p);
+    pr[2] = (long)gscalcol_i(p,N);
+    pr[5] = (long)gscalcol_i(gun,N);
   }
-  c=new_chunk(nbl+1); for (i=1; i<=nbl; i++) c[i]=0;
-  l=new_chunk(nbc+1);
-  d=new_chunk(nbc+1);
-  k = t = 1;
-  while (t<=nbl && k<=nbc)
-  {
-    for (j=1; j<k; j++)
-      for (i=1; i<=nbl; i++)
-	if (i!=l[j])
-          a[k][i] = (d[j]*a[k][i] - a[j][i]*a[k][l[j]]) % p;
-    t=1; while (t<=nbl && (c[t] || !a[k][t])) t++;
-    if (t<=nbl) { d[k]=a[k][t]; c[t]=k; l[k++]=t; }
+  else
+  { /* make sure v_pr(u) = 1 (automatic if e>1) */
+    if (is_pm1(e) && !is_uniformizer(NULL,u, T, gpowgs(p,f+1)))
+      u[2] = laddii((GEN)u[2], p);
+    t = algtobasis_i(nf, FpX_div(T,u,p));
+    pr[2] = (long)algtobasis_i(nf,u);
+    pr[5] = (long)centermod(t, p);
   }
-  if (k>nbc) err(bugparier,"kerlens");
-  avma=av; y=cgetg(nbc+1,t_COL);
-  t=(k>1) ? a[k][l[1]]:1;
-  y[1]=(t>0)? lstoi(t):lstoi(t+p);
-  for (q=1,j=2; j<k; j++)
-  {
-    q = (q*d[j-1]) % p;
-    t = (a[k][l[j]]*q) % p;
-    y[j] = (t>0) ? lstoi(t) : lstoi(t+p);
-  }
-  if (k>1)
-  {
-    t = (q*d[k-1]) % p;
-    y[k] = (t>0) ? lstoi(p-t) : lstoi(-t);
-  }
-  for (j=k+1; j<=nbc; j++) y[j]=zero;
-  return y;
+  return pr;
 }
 
-/* Calcule la constante de lenstra de l'ideal p.Z_K+a.Z_K ou a est un
-vecteur sur la base d'entiers */
+/* return a Z basis of Z_K's p-radical, phi = x--> x^p-x */
 static GEN
-lens(GEN nf, GEN p, GEN a)
+pradical(GEN nf, GEN p, GEN *phi)
 {
-  long av=avma,tetpil,N=degpol(nf[1]),j;
-  GEN mat=cgetg(N+1,t_MAT);
-  for (j=1; j<=N; j++) mat[j]=(long)element_mulid(nf,a,j);
-  tetpil=avma; return gerepile(av,tetpil,kerlens(mat,p));
-}
+  long i,N = degpol(nf[1]);
+  GEN q,m,frob,rad;
 
-extern GEN det_mod_P_n(GEN a, GEN N, GEN P);
-GEN sylvestermatrix_i(GEN x, GEN y);
+  /* matrix of Frob: x->x^p over Z_K/p */
+  frob = cgetg(N+1,t_MAT);
+  for (i=1; i<=N; i++)
+    frob[i] = (long)element_powid_mod_p(nf,i,p, p);
 
-/* check if p^va doesnt divide norm x (or norm(x+p)) */
-#if 0
-/* compute norm mod p^whatneeded using Sylvester's matrix */
-/* looks slower than the new subresultant. Have to re-check this */
+  m = frob; q = p;
+  while (cmpis(q,N) < 0) { q = mulii(q,p); m = FpM_mul(m, frob, p); }
+  rad = FpM_ker(m, p); /* m = Frob^k, s.t p^k >= N */
+  for (i=1; i<=N; i++)
+    coeff(frob,i,i) = lsubis(gcoeff(frob,i,i), 1);
+  *phi = frob; return rad;
+}
+
+/* return powers of a: a^0, ... , a^d,  d = dim A */
 static GEN
-prime_check_elt(GEN a, GEN pol, GEN p, GEN pf)
+get_powers(GEN mul, GEN p)
 {
-  GEN M,mod,x, c = denom(content(a));
-  long v = pvaluation(c, p, &x); /* x is junk */
+  long i, d = lg(mul[1]);
+  GEN z, pow = cgetg(d+2,t_MAT), P = pow+1;
 
-  mod = mulii(pf, gpowgs(p, degpol(pol)*v + 1));
-
-  x = FpX_red(gmul(a,c), mod);
-  M = sylvestermatrix_i(pol,x);
-  if (det_mod_P_n(M,mod,p) == gzero)
+  P[0] = (long)gscalcol_i(gun, d-1);
+  z = (GEN)mul[1];
+  for (i=1; i<=d; i++)
   {
-    x[2] = ladd((GEN)x[2], mulii(p,c));
-    M = sylvestermatrix_i(pol,x);
-    if (det_mod_P_n(M,mod,p) == gzero) return NULL;
-    a[2] = ladd((GEN)a[2], p);
+    P[i] = (long)z; /* a^i */
+    if (i!=d) z = FpM_FpV_mul(mul, z, p);
   }
-  return a;
+  return pow;
 }
-#else
-/* use subres to compute norm */
+
+/* minimal polynomial of a in A (dim A = d).
+ * mul = multiplication table by a in A */
 static GEN
-prime_check_elt(GEN a, GEN pol, GEN p, GEN pf)
+pol_min(GEN mul, GEN p)
 {
-  GEN norme=subres(pol,a);
-  if (resii(divii(norme,pf),p) != gzero) return a;
-  /* Note: a+p can't succeed if e > 1, can we know this at this point ? */
-  a=gadd(a,p); norme=subres(pol,a);
-  if (resii(divii(norme,pf),p) != gzero) return a;
-  return NULL;
+  gpmem_t av = avma;
+  GEN z, pow = get_powers(mul, p);
+  z = FpM_deplin(pow, p);
+  return gerepileupto(av, gtopolyrev(z,0));
 }
-#endif
 
-#if 0
 static GEN
-prime_two_elt_loop(GEN beta, GEN pol, GEN p, GEN pf)
+get_pr(GEN nf, GEN L, long i, GEN p, int appr)
 {
-  long av, m = lg(beta)-1;
-  int i,j,K, *x = (int*)new_chunk(m+1);
-  GEN a;
+  GEN pr, u, t, P = (GEN)L[i];
+  long e, f;
+  gpmem_t av;
 
-  K = 1; av = avma;
-  for(;;)
-  { /* x runs through strictly increasing sequences of length K,
-     * 1 <= x[i] <= m */
-nextK:
-    if (DEBUGLEVEL) fprintferr("K = %d\n", K);
-    for (i=1; i<=K; i++) x[i] = i;
-    for(;;)
-    {
-      if (DEBUGLEVEL > 1)
-      {
-        for (i=1; i<=K; i++) fprintferr("%d ",x[i]);
-        fprintferr("\n"); flusherr();
-      }
-      a = (GEN)beta[x[1]];
-      for (i=2; i<=K; i++) a = gadd(a, (GEN)beta[x[i]]);
-      if ((a = prime_check_elt(a,pol,p,pf))) return a;
-      avma = av;
+  if (typ(P) == t_VEC) return P; /* already done (Kummer) */
 
-      /* start: i = K+1; */
-      do
-      {
-        if (--i == 0)
-        {
-          if (++K > m) return NULL; /* fail */
-          goto nextK;
-        }
-        x[i]++;
-      } while (x[i] > m - K + i);
-      for (j=i; j<K; j++) x[j+1] = x[j]+1;
-    }
-  }
+  av = avma;
+  if (appr) u = uniformizer_appr(nf, L, i, p);
+  else      u = uniformizer(nf, P, p);
+  u = gerepilecopy(av, u);
+  t = anti_uniformizer(nf,p,u);
+  av = avma;
+  e = 1 + int_elt_val(nf,t,p,t,NULL);
+  f = degpol(nf[1]) - (lg(P)-1);
+  avma = av;
+  pr = cgetg(6,t_VEC);
+  pr[1] = (long)p;
+  pr[2] = (long)u;
+  pr[3] = lstoi(e);
+  pr[4] = lstoi(f);
+  pr[5] = (long)t; return pr;
 }
-#endif
 
-static GEN
-random_prime_two_elt_loop(GEN beta, GEN pol, GEN p, GEN pf)
+static int
+use_appr(GEN L, GEN pp, long N)
 {
-  long av = avma, z,i, m = lg(beta)-1;
-  long keep = getrand();
-  int c = 0;
-  GEN a;
+  long i, f, l = lg(L);
+  double prod = 1., NP;
+  GEN P;
+  gpmem_t av = avma;
 
-  for(i=1; i<=m; i++)
-    if ((a = prime_check_elt((GEN)beta[i],pol,p,pf))) return a;
-  (void)setrand(1);
-  if (DEBUGLEVEL) fprintferr("prime_two_elt_loop, hard case: ");
-  for(;;avma=av)
+  for (i=1; i<l; i++)
   {
-    if (DEBUGLEVEL) fprintferr("%d ", ++c);
-    a = gzero;
-    for (i=1; i<=m; i++)
-    {
-      z = mymyrand() >> (BITS_IN_RANDOM-5); /* in [0,15] */
-      if (z >= 9) z -= 7;
-      a = gadd(a,gmulsg(z,(GEN)beta[i]));
-    }
-    if ((a = prime_check_elt(a,pol,p,pf)))
-    {
-      if (DEBUGLEVEL) fprintferr("\n");
-      (void)setrand(keep); return a;
-    }
+    P = (GEN)L[i];
+    f = typ(P)==t_VEC? itos((GEN)P[4]): N - (lg(P)-1);
+    NP = gtodouble(gpowgs(pp, f));
+    prod *= (1 - 1./NP);
   }
+  avma = av;
+  return (prod * N * 10 < 1);
 }
 
-/* Input: an ideal mod p (!= Z_K)
- * Output: a 2-elt representation [p, x] */
+/* prime ideal decomposition of p */
 static GEN
-prime_two_elt(GEN nf, GEN p, GEN ideal)
+_primedec(GEN nf, GEN p)
 {
-  GEN beta,a,pf, pol = (GEN)nf[1];
-  long f, N=degpol(pol), m=lg(ideal)-1;
-  ulong av;
+  long i,k,c,iL,N;
+  GEN ex,F,L,Lpr,Ip,H,phi,mat1,T,f,g,h,p1,UN;
+  int appr;
 
-  if (!m) return gscalcol_i(p,N);
+  if (DEBUGLEVEL>2) (void)timer2();
+  nf = checknf(nf); T = (GEN)nf[1];
+  F = factmod(T,p);
+  ex = (GEN)F[2];
+  F  = (GEN)F[1]; F = centerlift(F);
+  if (DEBUGLEVEL>5) msgtimer("factmod");
 
-  /* we want v_p(Norm(beta)) = p^f, f = N-m */
-  av = avma; f = N-m; pf = gpuigs(p,f);
-  ideal = centerlift(ideal);
-  ideal = concatsp(gscalcol(p,N), ideal);
-  ideal = ideal_better_basis(nf, ideal, p);
-  beta = gmul((GEN)nf[7], ideal);
+  k = lg(F);
+  if (signe(modii((GEN)nf[4],p))) /* p doesn't divide index */
+  {
+    L = cgetg(k,t_VEC);
+    for (i=1; i<k; i++)
+      L[i] = (long)apply_kummer(nf,(GEN)F[i],(GEN)ex[i],p);
+    if (DEBUGLEVEL>5) msgtimer("simple primedec");
+    return L;
+  }
 
-#if 0
-  a = prime_two_elt_loop(beta,pol,p,pf);
-  if (!a) err(bugparier, "prime_two_elt (failed)");
-#else
-  a = random_prime_two_elt_loop(beta,pol,p,pf);
-#endif
+  g = (GEN)F[1];
+  for (i=2; i<k; i++) g = FpX_mul(g,(GEN)F[i], p);
+  h = FpX_div(T,g,p);
+  f = FpX_red(gdivexact(gsub(gmul(g,h), T), p), p);
 
-  a = centermod(algtobasis_intern(nf,a), p);
-  if (resii(divii(subres(gmul((GEN)nf[7],a),pol),pf),p) == gzero)
-    a[1] = laddii((GEN)a[1],p);
-  return gerepilecopy(av,a);
-}
+  N = degpol(T);
+  L = cgetg(N+1,t_VEC); iL = 1;
+  for (i=1; i<k; i++)
+    if (is_pm1(ex[i]) || signe(FpX_res(f,(GEN)F[i],p)))
+      L[iL++] = (long)apply_kummer(nf,(GEN)F[i],(GEN)ex[i],p);
+    else /* F[i] | (f,g,h), happens at least once by Dedekind criterion */
+      ex[i] = 0;
+  if (DEBUGLEVEL>2) msgtimer("%ld Kummer factors", iL-1);
 
-static GEN
-apply_kummer(GEN nf,GEN pol,GEN e,GEN p,long N,GEN *beta)
-{
-  GEN T,p1, res = cgetg(6,t_VEC);
-  long f = degpol(pol);
+  /* phi matrix of x -> x^p - x in algebra Z_K/p */
+  Ip = pradical(nf,p,&phi);
+  if (DEBUGLEVEL>2) msgtimer("pradical");
 
-  res[1]=(long)p;
-  res[3]=(long)e;
-  res[4]=lstoi(f);
-  if (f == N) /* inert */
-  {
-    res[2]=(long)gscalcol_i(p,N);
-    res[5]=(long)gscalcol_i(gun,N);
+  /* split etale algebra Z_K / (p,Ip) */
+  h = cgetg(N+1,t_VEC);
+  if (iL > 1)
+  { /* split off Kummer factors */
+    GEN mulbeta, beta = NULL;
+    for (i=1; i<k; i++)
+      if (!ex[i]) beta = beta? FpX_mul(beta, (GEN)F[i], p): (GEN)F[i];
+    beta = FpV_red(algtobasis_i(nf,beta), p);
+
+    mulbeta = FpM_red(eltmul_get_table(nf, beta), p);
+    p1 = concatsp(mulbeta, Ip);
+    /* Fp-base of ideal (Ip, beta) in ZK/p */
+    h[1] = (long)FpM_image(p1, p);
   }
   else
-  {
-    T = (GEN) nf[1];
-    if (ggval(subres(pol,T),p) > f)
-      pol[2] = laddii((GEN)pol[2],p);
-    res[2] = (long) algtobasis_intern(nf,pol);
+    h[1] = (long)Ip;
 
-    p1 = FpX_div(T,pol,p);
-    res[5] = (long) centermod(algtobasis_intern(nf,p1), p);
+  UN = gscalcol(gun, N);
+  for (c=1; c; c--)
+  { /* Let A:= (Z_K/p) / Ip; try to split A2 := A / Im H ~ Im M2
+       H * ? + M2 * Mi2 = Id_N ==> M2 * Mi2 projector A --> A2 */
+    GEN M, Mi, M2, Mi2, phi2;
+    long dim;
 
-    if (beta)
-      *beta = *beta? FpX_div(*beta,pol,p): p1;
+    H = (GEN)h[c]; k = lg(H)-1;
+    M   = FpM_suppl(concatsp(H,UN), p);
+    Mi  = FpM_inv(M, p);
+    M2  = vecextract_i(M, k+1,N); /* M = (H|M2) invertible */
+    Mi2 = rowextract_i(Mi,k+1,N);
+    /* FIXME: FpM_mul(,M2) could be done with vecextract_p */
+    phi2 = FpM_mul(Mi2, FpM_mul(phi,M2, p), p);
+    mat1 = FpM_ker(phi2, p);
+    dim = lg(mat1)-1; /* A2 product of 'dim' fields */
+    if (dim > 1)
+    { /* phi2 v = 0 <==> a = M2 v in Ker phi */
+      GEN I,R,r,a,mula,mul2, v = (GEN)mat1[2];
+      long n;
+
+      a = FpM_FpV_mul(M2,v, p);
+      mula = FpM_red(eltmul_get_table(nf, a), p);
+      mul2 = FpM_mul(Mi2, FpM_mul(mula,M2, p), p);
+      R = rootmod(pol_min(mul2,p), p); /* totally split mod p */
+
+      n = lg(R)-1;
+      for (i=1; i<=n; i++)
+      {
+        r = lift_intern((GEN)R[i]);
+        I = gaddmat_i(negi(r), mula);
+	h[c++] = (long)FpM_image(concatsp(H, I), p);
+      }
+      if (n == dim)
+        for (i=1; i<=n; i++) { H = (GEN)h[--c]; L[iL++] = (long)H; }
+    }
+    else /* A2 field ==> H maximal, f = N-k = dim(A2) */
+      L[iL++] = (long)H;
   }
-  return res;
+  setlg(L, iL);
+  Lpr = cgetg(iL, t_VEC);
+  if (DEBUGLEVEL>2) msgtimer("splitting %ld factors",iL-1);
+  appr = use_appr(L,p,N);
+  if (DEBUGLEVEL>2 && appr) fprintferr("using the approximation theorem\n");
+  for (i=1; i<iL; i++) Lpr[i] = (long)get_pr(nf, L, i, p, appr);
+  if (DEBUGLEVEL>2) msgtimer("finding uniformizers");
+  return Lpr;
 }
 
-/* prime ideal decomposition of p sorted by increasing residual degree */
 GEN
 primedec(GEN nf, GEN p)
 {
-  long av=avma,tetpil,i,j,k,kbar,np,c,indice,N,lp;
-  GEN ex,f,list,ip,elth,h;
-  GEN modfrob,algebre,algebre1,b,mat1,T;
-  GEN alpha,beta,p1,p2,unmodp,zmodp,idmodp;
+  gpmem_t av = avma;
+  return gerepileupto(av, gen_sort(_primedec(nf,p), 0, cmp_prime_over_p));
+}
 
-  if (DEBUGLEVEL>=3) timer2();
-  nf=checknf(nf); T=(GEN)nf[1]; N=degpol(T);
-  f=factmod(T,p); ex=(GEN)f[2];
-  f=centerlift((GEN)f[1]); np=lg(f);
-  if (DEBUGLEVEL>=6) msgtimer("factmod");
+/* return [Fp[x]: Fp] */
+static long
+ffdegree(GEN x, GEN frob, GEN p)
+{
+  gpmem_t av = avma;
+  long d, f = lg(frob)-1;
+  GEN y = x;
 
-  if (signe(modii((GEN)nf[4],p))) /* p doesn't divide index */
+  for (d=1; d < f; d++)
   {
-    list=cgetg(np,t_VEC);
-    for (i=1; i<np; i++)
-      list[i]=(long)apply_kummer(nf,(GEN)f[i],(GEN)ex[i],p,N, NULL);
-    if (DEBUGLEVEL>=6) msgtimer("simple primedec");
-    p1=stoi(4); tetpil=avma;
-    return gerepile(av,tetpil,vecsort(list,p1));
+    y = FpM_FpV_mul(frob, y, p);
+    if (gegal(y, x)) break;
   }
+  avma = av; return d;
+}
 
-  p1 = (GEN)f[1];
-  for (i=2; i<np; i++)
-    p1 = FpX_red(gmul(p1, (GEN)f[i]), p);
-  p1 = FpX_red(gdiv(gadd(gmul(p1, FpX_div(T,p1,p)), gneg(T)), p), p);
-  list = cgetg(N+1,t_VEC);
-  indice=1; beta=NULL;
-  for (i=1; i<np; i++) /* e = 1 or f[i] does not divide p1 (mod p) */
-    if (is_pm1(ex[i]) || signe(FpX_res(p1,(GEN)f[i],p)))
-      list[indice++] = (long)apply_kummer(nf,(GEN)f[i],(GEN)ex[i],p,N,&beta);
-  if (DEBUGLEVEL>=3) msgtimer("unramified factors");
+static GEN
+lift_to_zk(GEN v, GEN c, long N)
+{
+  GEN w = zerocol(N);
+  long i, l = lg(c);
+  for (i=1; i<l; i++) w[c[i]] = v[i];
+  return w;
+}
 
-  /* modfrob = modified Frobenius: x -> x^p - x mod p */
-  ip = pradical(nf,p,&modfrob);
-  if (DEBUGLEVEL>=3) msgtimer("pradical");
+/* 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), gpowgs(p,itos(e)-1));
+}
 
-  if (beta)
+/* return t = 1 mod pr, t = 0 mod p / pr^e(pr/p) */
+static GEN
+anti_uniformizer2(GEN nf, GEN pr)
+{
+  GEN p = (GEN)pr[1], z;
+  z = gmod(special_anti_uniformizer(nf, pr), p);
+  z = eltmul_get_table(nf, z);
+  z = hnfmodid(z, p);
+  z = idealaddtoone_i(nf,pr,z);
+  return unnf_minus_x(z);
+}
+
+#define mpr_TAU 1
+#define mpr_FFP 2
+#define mpr_PR  3
+#define mpr_T   4
+#define mpr_NFP 5
+#define SMALLMODPR 4
+#define LARGEMODPR 6
+static GEN
+modpr_TAU(GEN modpr)
+{
+  GEN tau = (GEN)modpr[mpr_TAU];
+  if (typ(tau) == t_INT && signe(tau) == 0) tau = NULL;
+  return tau;
+}
+
+/* prh = HNF matrix, which is identity but for the first line. Return a
+ * projector to ZK / prh ~ Z/prh[1,1] */
+GEN
+dim1proj(GEN prh)
+{
+  long i, N = lg(prh)-1;
+  GEN ffproj = cgetg(N+1, t_VEC);
+  GEN x, q = gcoeff(prh,1,1);
+  ffproj[1] = un;
+  for (i=2; i<=N; i++)
   {
-    beta = algtobasis_intern(nf,beta);
-    lp=lg(ip)-1; p1=cgetg(2*lp+N+1,t_MAT);
-    for (i=1; i<=N; i++) p1[i]=(long)element_mulid(nf,beta,i);
-    for (   ; i<=N+lp; i++)
-    {
-      p2 = (GEN) ip[i-N];
-      p1[i+lp] = (long) p2;
-      p1[i] = ldiv(element_mul(nf,lift(p2),beta),p);
-    }
-    ip = FpM_image(p1, p);
-    if (lg(ip)>N) err(bugparier,"primedec (bad pradical)");
+    x = gcoeff(prh,1,i);
+    if (signe(x)) x = subii(q,x);
+    ffproj[i] = (long)x;
   }
-  unmodp=gmodulsg(1,p); zmodp=gmodulsg(0,p);
-  idmodp = idmat_intern(N,unmodp,zmodp);
-  ip = gmul(ip, unmodp);
+  return ffproj;
+}
 
-  h=cgetg(N+1,t_VEC); h[1]=(long)ip;
-  for (c=1; c; c--)
+static GEN
+modprinit(GEN nf, GEN pr, int zk)
+{
+  gpmem_t av = avma;
+  GEN res, tau, mul, x, p, T, pow, ffproj, nfproj, prh, c, gf;
+  long N, i, k, f;
+
+  nf = checknf(nf); checkprimeid(pr);
+  gf = (GEN)pr[4];
+  f = itos(gf);
+  N = degpol(nf[1]);
+  prh = prime_to_ideal(nf, pr);
+  tau = zk? gzero: anti_uniformizer2(nf, pr);
+  p = (GEN)pr[1];
+
+  if (f == 1)
   {
-    elth=(GEN)h[c]; k=lg(elth)-1; kbar=N-k;
-    p1 = concatsp(elth,(GEN)idmodp[1]);
-    algebre = suppl_intern(p1,idmodp);
-    algebre1 = cgetg(kbar+1,t_MAT);
-    for (i=1; i<=kbar; i++) algebre1[i]=algebre[i+k];
-    b = gmul(modfrob,algebre1);
-    for (i=1;i<=kbar;i++)
-      b[i] = (long) project(algebre,(GEN) b[i],k,kbar);
-    mat1 = FpM_ker(lift_intern(b), p);
-    if (lg(mat1)>2)
+    res = cgetg(SMALLMODPR, t_COL);
+    res[mpr_TAU] = (long)tau;
+    res[mpr_FFP] = (long)dim1proj(prh);
+    res[mpr_PR ] = (long)pr; return gerepilecopy(av, res);
+  }
+
+  c = cgetg(f+1, t_VECSMALL);
+  ffproj = cgetg(N+1, t_MAT);
+  for (k=i=1; i<=N; i++)
+  {
+    x = gcoeff(prh, i,i);
+    if (!is_pm1(x)) { c[k] = i; ffproj[i] = (long)_ei(N, i); k++; }
+    else
+      ffproj[i] = lneg((GEN)prh[i]);
+  }
+  ffproj = rowextract_p(ffproj, c);
+  if (! divise((GEN)nf[4], p))
+  {
+    GEN basis, proj_modT;
+    if (N == f) T = (GEN)nf[1]; /* pr inert */
+    else
+      T = primpart(gmul((GEN)nf[7], (GEN)pr[2]));
+
+    basis = (GEN)nf[7];
+    proj_modT = cgetg(f+1, t_MAT);
+    for (i=1; i<=f; i++)
     {
-      GEN mat2 = cgetg(k+N+1,t_MAT);
-      for (i=1; i<=k; i++) mat2[i]=elth[i];
-      alpha=gmul(algebre1,(GEN)mat1[2]);
-      p1 = pol_min(alpha,nf,algebre,kbar,p);
-      p1 = (GEN)factmod(p1,p)[1];
-      for (i=1; i<lg(p1); i++)
+      GEN cx, w = Q_primitive_part((GEN)basis[c[i]], &cx);
+      w = FpX_rem(w, T, p);
+      if (cx)
       {
-	beta = eval_pol(nf,(GEN)p1[i],alpha,algebre,algebre1);
-        beta = lift_intern(beta);
-	for (j=1; j<=N; j++)
-	  mat2[k+j] = (long)FpV(element_mulid(nf,beta,j), p);
-	h[c] = (long) image(mat2); c++;
+        cx = gmod(cx, p);
+        w = FpX_Fp_mul(w, cx, p);
       }
+      proj_modT[i] = (long)pol_to_vec(w, f); /* w_i mod (T,p) */
     }
-    else
+    ffproj = FpM_mul(proj_modT,ffproj,p);
+
+    res = cgetg(SMALLMODPR+1, t_COL);
+    res[mpr_TAU] = (long)tau;
+    res[mpr_FFP] = (long)ffproj;
+    res[mpr_PR ] = (long)pr;
+    res[mpr_T  ] = (long)T; return gerepilecopy(av, res);
+  }
+
+  if (isprime(gf))
+  {
+    mul = eltmulid_get_table(nf, c[2]);
+    mul = vecextract_p(mul, c);
+  }
+  else
+  {
+    GEN v, u, u2, frob;
+    long deg,deg1,deg2;
+
+    /* matrix of Frob: x->x^p over Z_K/pr = < w[c1], ..., w[cf] > over Fp */
+    frob = cgetg(f+1, t_MAT);
+    for (i=1; i<=f; i++)
     {
-      long av1; p1 = cgetg(6,t_VEC);
-      list[indice++] = (long)p1;
-      p1[1]=(long)p; p1[4]=lstoi(kbar);
-      p1[2]=(long)prime_two_elt(nf,p,elth);
-      p1[5]=(long)lens(nf,p,(GEN)p1[2]);
-      av1=avma;
-      i = int_elt_val(nf,(GEN)p1[5],p,(GEN)p1[5],NULL,N-1);
-      avma=av1;
-      p1[3]=lstoi(i+1);
+      x = element_powid_mod_p(nf,c[i],p, p);
+      frob[i] = (long)FpM_FpV_mul(ffproj, x, p);
     }
-    if (DEBUGLEVEL>=3) msgtimer("h[%ld]",c);
+    u = _ei(f,2); k = 2;
+    deg1 = ffdegree(u, frob, p);
+    while (deg1 < f)
+    {
+      k++; u2 = _ei(f, k);
+      deg2 = ffdegree(u2, frob, p);
+      deg = clcm(deg1,deg2);
+      if (deg == deg1) continue;
+      if (deg == deg2) { deg1 = deg2; u = u2; continue; }
+      u = gadd(u, u2);
+      while (ffdegree(u, frob, p) < deg) u = gadd(u, u2);
+      deg1 = deg;
+    }
+    v = lift_to_zk(u,c,N);
+
+    mul = cgetg(f+1,t_MAT);
+    mul[1] = (long)v; /* assume w_1 = 1 */
+    for (i=2; i<=f; i++) mul[i] = (long)element_mulid(nf,v,c[i]);
   }
-  setlg(list, indice); tetpil=avma;
-  return gerepile(av,tetpil,gen_sort(list,0,cmp_prime_over_p));
+
+  /* Z_K/pr = Fp(v), mul = mul by v */
+  mul = FpM_red(mul, p);
+  mul = FpM_mul(ffproj, mul, p);
+
+  pow = get_powers(mul, p);
+  T = gtopolyrev(FpM_deplin(pow, p), varn(nf[1]));
+  nfproj = cgetg(f+1, t_MAT);
+  for (i=1; i<=f; i++) nfproj[i] = (long)lift_to_zk((GEN)pow[i], c, N);
+  nfproj = gmul((GEN)nf[7], nfproj);
+
+  setlg(pow, f+1);
+  ffproj = FpM_mul(FpM_inv(pow, p), ffproj, p);
+
+  res = cgetg(LARGEMODPR, t_COL);
+  res[mpr_TAU] = (long)tau;
+  res[mpr_FFP] = (long)ffproj;
+  res[mpr_PR ] = (long)pr;
+  res[mpr_T  ] = (long)T;
+  res[mpr_NFP] = (long)nfproj; return gerepilecopy(av, res);
 }
 
+GEN
+nfmodprinit(GEN nf, GEN pr) { return modprinit(nf, pr, 0); }
+GEN
+zkmodprinit(GEN nf, GEN pr) { return modprinit(nf, pr, 1); }
+
+void
+checkmodpr(GEN modpr)
+{
+  if (typ(modpr) != t_COL || lg(modpr) < SMALLMODPR)
+    err(talker,"incorrect modpr format");
+  checkprimeid((GEN)modpr[mpr_PR]);
+}
+
+
+static GEN
+to_ff_init(GEN nf, GEN *pr, GEN *T, GEN *p, int zk)
+{
+  GEN modpr = (typ(*pr) == t_COL)? *pr: modprinit(nf, *pr, zk);
+  *T = lg(modpr)==SMALLMODPR? NULL: (GEN)modpr[mpr_T];
+  *pr = (GEN)modpr[mpr_PR];
+  *p = (GEN)(*pr)[1]; return modpr;
+}
+GEN
+nf_to_ff_init(GEN nf, GEN *pr, GEN *T, GEN *p) {
+  GEN modpr = to_ff_init(nf,pr,T,p,0);
+  GEN tau = modpr_TAU(modpr);
+  if (!tau) modpr[mpr_TAU] = (long)anti_uniformizer2(nf, *pr);
+  return modpr;
+}
+GEN
+zk_to_ff_init(GEN nf, GEN *pr, GEN *T, GEN *p) {
+  return to_ff_init(nf,pr,T,p,1);
+}
+
+/* assume x in 'basis' form (t_COL) */
+GEN
+zk_to_ff(GEN x, GEN modpr)
+{
+  GEN pr = (GEN)modpr[mpr_PR];
+  GEN p = (GEN)pr[1];
+  GEN y = gmul((GEN)modpr[mpr_FFP], x);
+  if (lg(modpr) == SMALLMODPR) return modii(y,p);
+  y = FpV_red(y, p);
+  return col_to_ff(y, varn(modpr[mpr_T]));
+}
+
 /* REDUCTION Modulo a prime ideal */
 
+/* assume x in t_COL form, v_pr(x) >= 0 */
+static GEN
+kill_denom(GEN x, GEN nf, GEN p, GEN modpr)
+{
+  GEN cx, den = denom(x);
+  long v;
+  if (gcmp1(den)) return x;
+
+  v = ggval(den,p);
+  if (v) 
+  {
+    GEN tau = modpr_TAU(modpr);
+    if (!tau) err(talker,"modpr initialized for integers only!");
+    x = element_mul(nf,x, element_pow(nf,tau,stoi(v)));
+  }
+  x = Q_primitive_part(x, &cx);
+  x = FpV_red(x, p);
+  if (cx) x = FpV_red(gmul(gmod(cx, p), x), p);
+  return x;
+}
+
 /* x integral, reduce mod prh in HNF */
 GEN
 nfreducemodpr_i(GEN x, GEN prh)
@@ -2323,35 +2615,116 @@ nfreducemodpr_i(GEN x, GEN prh)
   x[1] = lresii((GEN)x[1], p); return x;
 }
 
-/* for internal use */
 GEN
-nfreducemodpr(GEN nf, GEN x, GEN prhall)
+nfreducemodpr(GEN nf, GEN x, GEN modpr)
 {
-  long i,v;
-  GEN p,prh,den;
+  gpmem_t av = avma;
+  long i;
+  GEN pr, p;
 
+  checkmodpr(modpr);
+  pr = (GEN)modpr[mpr_PR];
+  p = (GEN)pr[1];
+  if (typ(x) != t_COL) x = algtobasis(nf,x);
   for (i=lg(x)-1; i>0; i--)
-    if (typ(x[i]) == t_INTMOD) { x=lift_intern(x); break; }
-  prh=(GEN)prhall[1]; p=gcoeff(prh,1,1);
-  den=denom(x);
-  if (!gcmp1(den))
+    if (typ(x[i]) == t_INTMOD) { x = lift(x); break; }
+  x = kill_denom(x, nf, p, modpr);
+  x = ff_to_nf(zk_to_ff(x,modpr), modpr);
+  if (typ(x) != t_COL) x = algtobasis(nf, x);
+  return gerepileupto(av, FpV(x, p));
+}
+
+GEN
+nf_to_ff(GEN nf, GEN x, GEN modpr)
+{
+  gpmem_t av = avma;
+  GEN pr = (GEN)modpr[mpr_PR];
+  GEN p = (GEN)pr[1];
+  long t = typ(x);
+
+  if (t == t_POLMOD) { x = (GEN)x[2]; t = typ(x); }
+  switch(t)
   {
-    v=ggval(den,p);
-    if (v) x=element_mul(nf,x,element_pow(nf,(GEN)prhall[2],stoi(v)));
-    x = gmod(x,p);
+    case t_INT: return modii(x, p);
+    case t_FRAC: return gmod(x, p);
+    case t_POL: x = algtobasis(nf, x); break;
+    case t_COL: break;
+    default: err(typeer,"nf_to_ff");
   }
-  return FpV(nfreducemodpr_i(x, prh), p);
+  x = kill_denom(x, nf, p, modpr);
+  return gerepilecopy(av, zk_to_ff(x, modpr));
 }
 
-/* public function */
 GEN
-nfreducemodpr2(GEN nf, GEN x, GEN prhall)
+ff_to_nf(GEN x, GEN modpr)
 {
-  long av = avma; checkprhall(prhall);
-  if (typ(x) != t_COL) x = algtobasis(nf,x);
-  return gerepileupto(av, nfreducemodpr(nf,x,prhall));
+  if (lg(modpr) < LARGEMODPR) return x;
+  return mulmat_pol((GEN)modpr[mpr_NFP], x);
 }
+GEN
+modprM_lift(GEN x, GEN modpr)
+{
+  long i,j,h,l = lg(x);
+  GEN y = cgetg(l, t_MAT);
+  if (l == 1) return y;
 
+  h = lg(x[1]);
+  for (j=1; j<l; j++)
+  {
+    GEN p1 = cgetg(h,t_COL); y[j] = (long)p1;
+    for (i=1; i<h; i++) p1[i]=(long)ff_to_nf(gcoeff(x,i,j), modpr);
+  }
+  return y;
+}
+GEN
+modprX_lift(GEN x, GEN modpr)
+{
+  long i, l;
+  GEN z;
+
+  if (typ(x)!=t_POL) return gcopy(x); /* scalar */
+  l = lgef(x);
+  z = cgetg(l, t_POL); z[1] = x[1];
+  for (i=2; i<l; i++) z[i] = (long)ff_to_nf((GEN)x[i], modpr);
+  return z;
+}
+
+/* reduce the coefficients of pol modulo modpr */
+GEN
+modprX(GEN x, GEN nf,GEN modpr)
+{
+  long i, l;
+  GEN z;
+
+  if (typ(x)!=t_POL) return nf_to_ff(nf,x,modpr);
+  l = lgef(x);
+  z = cgetg(l,t_POL); z[1] = x[1];
+  for (i=2; i<l; i++) z[i] = (long)nf_to_ff(nf,(GEN)x[i],modpr);
+  return normalizepol(z);
+}
+GEN
+modprV(GEN z, GEN nf,GEN modpr)
+{
+  long i,l = lg(z);
+  GEN x;
+  x = cgetg(l,typ(z));
+  for (i=1; i<l; i++) x[i] = (long)nf_to_ff(nf,(GEN)z[i], modpr);
+  return x;
+}
+/* assume z a t_VEC/t_COL/t_MAT */
+GEN
+modprM(GEN z, GEN nf,GEN modpr)
+{
+  long i,l = lg(z), m;
+  GEN x;
+
+  if (typ(z) != t_MAT) return modprV(z,nf,modpr);
+  x = cgetg(l,t_MAT); if (l==1) return x;
+  m = lg((GEN)z[1]);
+  for (i=1; i<l; i++) x[i] = (long)modprV((GEN)z[i],nf,modpr);
+  return x;
+}
+
 /*                        relative ROUND 2
  *
  * input: nf = base field K
@@ -2363,8 +2736,7 @@ nfreducemodpr2(GEN nf, GEN x, GEN prhall)
  */
 
 /* given MODULES x and y by their pseudo-bases in HNF, gives a
- * pseudo-basis of the module generated by x and y. For internal use.
- */
+ * pseudo-basis of the module generated by x and y. */
 static GEN
 rnfjoinmodules(GEN nf, GEN x, GEN y)
 {
@@ -2372,116 +2744,63 @@ rnfjoinmodules(GEN nf, GEN x, GEN y)
   GEN p1,p2,z,Hx,Hy,Ix,Iy;
 
   if (x == NULL) return y;
-  Hx=(GEN)x[1]; lx=lg(Hx); Ix=(GEN)x[2];
-  Hy=(GEN)y[1]; ly=lg(Hy); Iy=(GEN)y[2];
+  Hx = (GEN)x[1]; lx = lg(Hx); Ix = (GEN)x[2];
+  Hy = (GEN)y[1]; ly = lg(Hy); Iy = (GEN)y[2];
   i = lx+ly-1;
   z = (GEN)gpmalloc(sizeof(long*)*(3+2*i));
   *z = evaltyp(t_VEC)|evallg(3);
   p1 =  z+3; z[1]=(long)p1; *p1 = evaltyp(t_MAT)|evallg(i);
   p2 = p1+i; z[2]=(long)p2; *p2 = evaltyp(t_VEC)|evallg(i);
 
-  for (i=1; i<lx; i++) { p1[i]=Hx[i]; p2[i]=Ix[i]; }
-  for (   ; i<lx+ly-1; i++) { p1[i]=Hy[i-lx+1]; p2[i]=Iy[i-lx+1]; }
+  for (i=1; i<lx; i++) { p1[i] = Hx[i]; p2[i] = Ix[i]; }
+  p1 += lx-1;
+  p2 += lx-1;
+  for (i=1; i<ly; i++) { p1[i] = Hy[i]; p2[i] = Iy[i]; }
   x = nfhermite(nf,z); free(z); return x;
 }
 
-/* a usage interne, pas de gestion de pile : x et y sont des vecteurs dont
- * les coefficients sont les composantes sur nf[7]; avec reduction mod pr sauf
- * si prhall=NULL
- */
+typedef struct {
+  GEN nf, multab, modpr,T,p;
+  long h;
+} rnfeltmod_muldata;
+
 static GEN
-rnfelement_mulidmod(GEN nf, GEN multab, GEN unnf, GEN x, long h, GEN prhall)
+_mul(void *data, GEN x, GEN y/* base; ignored */)
 {
-  long j,k,N;
-  GEN p1,c,v,s,znf;
-
-  if (h==1) return gcopy(x);
-  N = lg(x)-1; multab += (h-1)*N;
-  x = lift(x); v = cgetg(N+1,t_COL);
-  znf = gscalcol_i(gzero,lg(unnf)-1);
-  for (k=1; k<=N; k++)
-  {
-    s = gzero;
-    for (j=1; j<=N; j++)
-      if (!gcmp0(p1 = (GEN)x[j]) && !gcmp0(c = gcoeff(multab,k,j)))
-      {
-        if (!gegal(c,unnf)) p1 = element_mul(nf,p1,c);
-        s = gadd(s,p1);
-      }
-    if (s == gzero) s = znf;
-    else
-      if (prhall) s = nfreducemodpr(nf,s,prhall);
-    v[k] = (long)s;
-  }
-  return v;
+  rnfeltmod_muldata *D = (rnfeltmod_muldata *) data;
+  GEN z = x? element_mulid(D->multab,x,D->h)
+           : element_mulidid(D->multab,D->h,D->h);
+  (void)y;
+  return FpVQX_red(z,D->T,D->p);
 }
-
-/* a usage interne, pas de gestion de pile : x est un vecteur dont
- * les coefficients sont les composantes sur nf[7]
- */
 static GEN
-rnfelement_sqrmod(GEN nf, GEN multab, GEN unnf, GEN x, GEN prhall)
+_sqr(void *data, GEN x)
 {
-  long i,j,k,n;
-  GEN p1,c,z,s;
-
-  n=lg(x)-1; x=lift(x); z=cgetg(n+1,t_COL);
-  for (k=1; k<=n; k++)
-  {
-    if (k == 1)
-      s = element_sqr(nf,(GEN)x[1]);
-    else
-      s = gmul2n(element_mul(nf,(GEN)x[1],(GEN)x[k]), 1);
-    for (i=2; i<=n; i++)
-    {
-      c = gcoeff(multab,k,(i-1)*n+i);
-      if (!gcmp0(c))
-      {
-	p1=element_sqr(nf,(GEN)x[i]);
-	if (!gegal(c,unnf)) p1 = element_mul(nf,p1,c);
-        s = gadd(s,p1);
-      }
-      for (j=i+1; j<=n; j++)
-      {
-	c = gcoeff(multab,k,(i-1)*n+j);
-	if (!gcmp0(c))
-	{
-	  p1=gmul2n(element_mul(nf,(GEN)x[i],(GEN)x[j]),1);
-	  if (!gegal(c,unnf)) p1 = element_mul(nf,p1,c);
-          s = gadd(s,p1);
-	}
-      }
-    }
-    if (prhall) s = nfreducemodpr(nf,s,prhall);
-    z[k]=(long)s;
-  }
-  return z;
+  rnfeltmod_muldata *D = (rnfeltmod_muldata *) data;
+  GEN z = x? sqr_by_tab(D->multab,x)
+           : element_mulidid(D->multab,D->h,D->h);
+  return FpVQX_red(z,D->T,D->p);
 }
 
-/* Compute x^n mod pr in the extension, assume n >= 0 [cf puissii()]*/
+/* Compute x^n mod pr in the extension, assume n >= 0 */
 static GEN
-rnfelementid_powmod(GEN nf, GEN multab, GEN matId, long h, GEN n, GEN prhall)
+rnfelementid_powmod(GEN multab, long h, GEN n, GEN T, GEN p)
 {
-  ulong av = avma;
-  long i,k,m;
-  GEN y,p1, unrnf=(GEN)matId[1], unnf=(GEN)unrnf[1];
+  gpmem_t av = avma;
+  GEN y;
+  rnfeltmod_muldata D;
 
-  if (!signe(n)) return unrnf;
-  y = (GEN)matId[h]; p1 = n+2; m = *p1;
-  k = 1+bfffo(m); m<<=k; k = BITS_IN_LONG-k;
-  for (i=lgefint(n)-2;;)
-  {
-    for (; k; m<<=1,k--)
-    {
-      y = rnfelement_sqrmod(nf,multab,unnf,y,prhall);
-      if (m < 0) y = rnfelement_mulidmod(nf,multab,unnf,y,h,prhall);
-    }
-    if (--i == 0) break;
-    m = *++p1; k = BITS_IN_LONG;
-  }
+  if (!signe(n)) return gun;
+
+  D.multab = multab;
+  D.h = h;
+  D.T = T;
+  D.p = p;
+  y = leftright_pow(NULL, n, (void*)&D, &_sqr, &_mul);
   return gerepilecopy(av, y);
 }
 
+#if 0
 GEN
 mymod(GEN x, GEN p)
 {
@@ -2501,158 +2820,251 @@ mymod(GEN x, GEN p)
   for (i=1; i<lx; i++) y[i] = (long)mymod((GEN)x[i],p);
   return y;
 }
+#endif
 
+#define FqX_div(x,y,T,p) FpXQX_divres((x),(y),(T),(p),NULL)
+#define FqX_mul FpXQX_mul
+
+/* relative Dedekind criterion over nf, applied to the order defined by a
+ * root of irreducible polynomial P, modulo the prime ideal pr. Returns
+ * [flag,basis,val], where basis is a pseudo-basis of the enlarged order,
+ * flag is 1 iff this order is pr-maximal, and val is the valuation at pr of
+ * the order discriminant */
 static GEN
-rnfordmax(GEN nf, GEN pol, GEN pr, GEN unnf, GEN id, GEN matId)
+rnfdedekind_i(GEN nf,GEN P,GEN pr, GEN disc)
 {
-  long av=avma,tetpil,av1,lim,i,j,k,n,v1,v2,vpol,m,cmpt,sep;
-  GEN p,q,q1,prhall,A,Aa,Aaa,A1,I,R,p1,p2,p3,multab,multabmod,Aainv;
+  long vt, r, d, n, m, i, j;
+  gpmem_t av = avma;
+  GEN Prd,p1,p2,p,tau,g,matid;
+  GEN modpr,res,h,k,base,nfT,T,gzk,hzk;
+
+  if (!disc) disc = discsr(P);
+  P = lift(P);
+  nf = checknf(nf); nfT = (GEN)nf[1];
+  res = cgetg(4,t_VEC);
+  modpr = nf_to_ff_init(nf,&pr, &T, &p);
+  tau = gmul((GEN)nf[7], (GEN)pr[5]);
+  n = degpol(nfT);
+  m = degpol(P);
+
+  Prd = modprX(P, nf, modpr);
+  p1 = (GEN)FqX_factor(Prd,T,p)[1];
+  r = lg(p1); if (r < 2) err(constpoler,"rnfdedekind");
+  g = (GEN)p1[1];
+  for (i=2; i<r; i++) g = FqX_mul(g, (GEN)p1[i], T, p);
+  h = FqX_div(Prd,g, T, p);
+  gzk = modprX_lift(g, modpr);
+  hzk = modprX_lift(h, modpr);
+
+  k = gsub(P, RXQX_mul(gzk,hzk, nfT));
+  k = gdiv(RXQX_mul(tau,k,nfT), p);
+  k = modprX(k, nf, modpr);
+
+  p2 = FqX_gcd(g,h,  T,p);
+  k  = FqX_gcd(p2,k, T,p); d = degpol(k);  /* <= m */
+
+  vt = idealval(nf,disc,pr) - 2*d;
+  res[3] = lstoi(vt);
+  res[1] = (!d || vt<=1)? un: zero;
+
+  base = cgetg(3,t_VEC);
+  p1 = cgetg(m+d+1,t_MAT); base[1] = (long)p1;
+  p2 = cgetg(m+d+1,t_VEC); base[2] = (long)p2;
+ /* if d > 0, base[2] temporarily multiplied by p, for the final nfhermitemod
+  * [ which requires integral ideals ] */
+  matid = gscalmat(d? p: gun, n);
+  for (j=1; j<=m; j++)
+  {
+    p1[j] = (long)_ei(m, j);
+    p2[j] = (long)matid;
+  }
+  if (d)
+  {
+    GEN prinvp = pidealprimeinv(nf,pr); /* again multiplied by p */
+    GEN X = polx[varn(P)], pal;
+
+    pal = FqX_div(Prd,k, T,p);
+    pal = modprX_lift(pal, modpr);
+    for (   ; j<=m+d; j++)
+    {
+      p1[j] = (long)pol_to_vec(pal,m);
+      p2[j] = (long)prinvp;
+      pal = RXQX_rem(RXQX_mul(pal,X,nfT),P,nfT);
+    }
+    /* the modulus is integral */
+    base = nfhermitemod(nf,base, gmul(gpowgs(p, m-d),
+				      idealpows(nf, prinvp, d)));
+    base[2] = ldiv((GEN)base[2], p); /* cancel the factor p */
+  }
+  res[2] = (long)base; return gerepilecopy(av, res);
+}
+
+GEN
+rnfdedekind(GEN nf,GEN P0,GEN pr)
+{
+  return rnfdedekind_i(nf,P0,pr,NULL);
+}
+
+static GEN
+rnfordmax(GEN nf, GEN pol, GEN pr, GEN disc)
+{
+  gpmem_t av=avma,av1,lim;
+  long i,j,k,n,v1,v2,vpol,m,cmpt,sep;
+  GEN p,T,q,q1,modpr,A,Aa,Aaa,A1,I,R,p1,p2,p3,multab,multabmod,Aainv;
   GEN pip,baseIp,baseOp,alpha,matprod,alphainv,matC,matG,vecpro,matH;
-  GEN neworder,H,Hid,alphalistinv,alphalist,prhinv;
+  GEN neworder,H,Hid,alphalistinv,alphalist,prhinv,nfT,id,rnfId;
 
   if (DEBUGLEVEL>1) fprintferr(" treating %Z\n",pr);
-  prhall=nfmodprinit(nf,pr);
-  q=cgetg(3,t_VEC); q[1]=(long)pr;q[2]=(long)prhall;
-  p1=rnfdedekind(nf,pol,q);
+  modpr = nf_to_ff_init(nf,&pr,&T,&p);
+  p1 = rnfdedekind_i(nf,pol,modpr,disc);
   if (gcmp1((GEN)p1[1])) return gerepilecopy(av,(GEN)p1[2]);
 
-  sep=itos((GEN)p1[3]);
-  A=gmael(p1,2,1);
-  I=gmael(p1,2,2);
+  sep = itos((GEN)p1[3]);
+  A = gmael(p1,2,1);
+  I = gmael(p1,2,2);
 
-  n=degpol(pol); vpol=varn(pol);
-  p=(GEN)pr[1]; q=powgi(p,(GEN)pr[4]); pip=(GEN)pr[2];
-  q1=q; while (cmpis(q1,n)<0) q1=mulii(q1,q);
+  pip = basistoalg(nf, (GEN)pr[2]);
+  nfT = (GEN)nf[1];
+  n = degpol(pol); vpol = varn(pol);
+  q = T? gpowgs(p,degpol(T)): p;
+  q1 = q; while (cmpis(q1,n) < 0) q1 = mulii(q1,q);
+  rnfId = idmat(degpol(pol));
+  id    = idmat(degpol(nfT));
 
-  multab=cgetg(n*n+1,t_MAT);
-  for (j=1; j<=n*n; j++) multab[j]=lgetg(n+1,t_COL);
-  prhinv = idealinv(nf,(GEN)prhall[1]);
-  alphalistinv=cgetg(n+1,t_VEC);
-  alphalist=cgetg(n+1,t_VEC);
-  A1=cgetg(n+1,t_MAT);
-  av1=avma; lim=stack_lim(av1,1);
+  multab = cgetg(n*n+1,t_MAT);
+  for (j=1; j<=n*n; j++) multab[j] = lgetg(n+1,t_COL);
+  prhinv = idealinv(nf, pr);
+  alphalistinv = cgetg(n+1,t_VEC);
+  alphalist    = cgetg(n+1,t_VEC);
+  A1 = cgetg(n+1,t_MAT);
+  av1 = avma; lim = stack_lim(av1,1);
   for(cmpt=1; ; cmpt++)
   {
-    if (DEBUGLEVEL>1)
+    if (DEBUGLEVEL>1) fprintferr("    %ld%s pass\n",cmpt,eng_ord(cmpt));
+    for (j=1; j<=n; j++)
     {
-      fprintferr("    %ld%s pass\n",cmpt,eng_ord(cmpt));
-      flusherr();
-    }
-    for (i=1; i<=n; i++)
-    {
-      if (gegal((GEN)I[i],id)) alphalist[i]=alphalistinv[i]=(long)unnf;
+      if (gegal((GEN)I[j],id))
+      {
+        alphalist[j] = alphalistinv[j] = un;
+        A1[j] = A[j]; continue;
+      }
+
+      p1 = ideal_two_elt(nf,(GEN)I[j]);
+      if (gcmp0((GEN)p1[1]))
+        alpha = (GEN)p1[2];
       else
       {
-	p1=ideal_two_elt(nf,(GEN)I[i]);
-	v1=gcmp0((GEN)p1[1])? EXP220
-                            : element_val(nf,(GEN)p1[1],pr);
-	v2=element_val(nf,(GEN)p1[2],pr);
-	if (v1>v2) p2=(GEN)p1[2]; else p2=(GEN)p1[1];
-	alphalist[i]=(long)p2;
-        alphalistinv[i]=(long)element_inv(nf,p2);
+        v1 = element_val(nf,(GEN)p1[1],pr);
+        v2 = element_val(nf,(GEN)p1[2],pr);
+        alpha = (v1>v2)? (GEN)p1[2]: (GEN)p1[1];
       }
-    }
-    for (j=1; j<=n; j++)
-    {
-      p1=cgetg(n+1,t_COL); A1[j]=(long)p1;
+      alphalist[j]    = (long)alpha;
+      alphalistinv[j] = (long)element_inv(nf,alpha);
+
+      p1 = cgetg(n+1,t_COL); A1[j] = (long)p1;
       for (i=1; i<=n; i++)
-	p1[i] = (long)element_mul(nf,gcoeff(A,i,j),(GEN)alphalist[j]);
+        p1[i] = (long)element_mul(nf,gcoeff(A,i,j),alpha);
     }
-    Aa=basistoalg(nf,A1); Aainv=lift_intern(ginv(Aa));
-    Aaa=mat_to_vecpol(Aa,vpol);
-    for (i=1; i<=n; i++) for (j=i; j<=n; j++)
-    {
-      long tp;
-      p1 = lift_intern(gres(gmul((GEN)Aaa[i],(GEN)Aaa[j]), pol));
-      tp = typ(p1);
-      if (is_scalar_t(tp) || (tp==t_POL && varn(p1)>vpol))
-	p2 = gmul(p1, (GEN)Aainv[1]);
-      else
-        p2 = gmul(Aainv, pol_to_vec(p1, n));
-
-      p3 = algtobasis(nf,p2);
-      for (k=1; k<=n; k++)
+    Aa = basistoalg(nf,A1);
+    Aainv = lift_intern(ginv(Aa));
+    Aaa = lift_intern(mat_to_vecpol(Aa,vpol));
+    for (i=1; i<=n; i++)
+      for (j=i; j<=n; j++)
       {
-	coeff(multab,k,(i-1)*n+j) = p3[k];
-	coeff(multab,k,(j-1)*n+i) = p3[k];
+        long tp;
+        p1 = RXQX_rem(gmul((GEN)Aaa[i],(GEN)Aaa[j]), pol, nfT);
+        tp = typ(p1);
+        if (is_scalar_t(tp) || (tp==t_POL && varn(p1)>vpol))
+          p2 = gmul(p1, (GEN)Aainv[1]);
+        else
+          p2 = mulmat_pol(Aainv, p1);
+        for (k=1; k<=n; k++)
+        {
+          GEN c = gres((GEN)p2[k], nfT);
+          coeff(multab,k,(i-1)*n+j) = (long)c;
+          coeff(multab,k,(j-1)*n+i) = (long)c;
+        }
       }
-    }
-    R=cgetg(n+1,t_MAT); multabmod = mymod(multab,p);
-    R[1] = matId[1];
+    multabmod = modprM(multab,nf,modpr);
+    R = cgetg(n+1,t_MAT);
+    R[1] = rnfId[1];
     for (j=2; j<=n; j++)
-      R[j] = (long) rnfelementid_powmod(nf,multabmod,matId, j,q1,prhall);
-    baseIp = nfkermodpr(nf,R,prhall);
-    baseOp = lift_intern(nfsuppl(nf,baseIp,n,prhall));
-    alpha=cgetg(n+1,t_MAT);
-    for (j=1; j<lg(baseIp); j++) alpha[j]=baseOp[j];
+      R[j] = (long) rnfelementid_powmod(multabmod,j,q1,T,p);
+    baseIp = FqM_ker(R,T,p);
+    baseOp = lg(baseIp)==1? rnfId: FqM_suppl(baseIp,T,p);
+    alpha = modprM_lift(baseOp, modpr);
+    for (j=1; j<lg(baseIp); j++)
+    {
+      p1 = (GEN)alpha[j];
+      for (i=1; i<=n; i++) p1[i] = (long)to_polmod((GEN)p1[i], nfT);
+    }
     for (   ; j<=n; j++)
     {
-      p1=cgetg(n+1,t_COL); alpha[j]=(long)p1;
-      for (i=1; i<=n; i++)
-	p1[i]=(long)element_mul(nf,pip,gcoeff(baseOp,i,j));
+      p1 = (GEN)alpha[j];
+      for (i=1; i<=n; i++) p1[i] = lmul(pip, (GEN)p1[i]);
     }
-    matprod=cgetg(n+1,t_MAT);
+    alphainv = lift_intern(ginv(alpha));
+    alpha = lift_intern(alpha);
+
+    matprod = cgetg(n+1,t_MAT);
     for (j=1; j<=n; j++)
     {
-      p1=cgetg(n+1,t_COL); matprod[j]=(long)p1;
+      p1 = cgetg(n+1,t_COL); matprod[j] = (long)p1;
       for (i=1; i<=n; i++)
-      {
-        p2 = rnfelement_mulidmod(nf,multab,unnf, (GEN)alpha[i],j, NULL);
-        for (k=1; k<=n; k++)
-          p2[k] = lmul((GEN)nf[7], (GEN)p2[k]);
-        p1[i] = (long)p2;
-      }
+        p1[i] = (long)gmod(element_mulid(multab, (GEN)alpha[i],j), nfT);
     }
-    alphainv = lift_intern(ginv(basistoalg(nf,alpha)));
     matC = cgetg(n+1,t_MAT);
     for (j=1; j<=n; j++)
     {
-      p1=cgetg(n*n+1,t_COL); matC[j]=(long)p1;
+      p1 = cgetg(n*n+1,t_COL); matC[j] = (long)p1;
       for (i=1; i<=n; i++)
       {
 	p2 = gmul(alphainv, gcoeff(matprod,i,j));
 	for (k=1; k<=n; k++)
-	  p1[(i-1)*n+k]=(long)nfreducemodpr(nf,algtobasis(nf,(GEN)p2[k]),prhall);
+        {
+          GEN c = gres((GEN)p2[k], nfT);
+          p1[(i-1)*n+k] = (long)nf_to_ff(nf,c,modpr);
+        }
       }
     }
-    matG=nfkermodpr(nf,matC,prhall); m=lg(matG)-1;
-    vecpro=cgetg(3,t_VEC);
-    p1=cgetg(n+m+1,t_MAT); vecpro[1]=(long)p1;
-    p2=cgetg(n+m+1,t_VEC); vecpro[2]=(long)p2;
-    for (j=1; j<=m; j++)
-    {
-      p1[j] = llift((GEN)matG[j]);
-      p2[j] = (long)prhinv;
-    }
-    p1 += m;
+    matG = FqM_ker(matC,T,p); m = lg(matG)-1;
+    matG = modprM_lift(matG, modpr);
+    vecpro = cgetg(3,t_VEC);
+    vecpro[1] = (long)concatsp(matG,rnfId);
+
+    p2 = cgetg(n+m+1,t_VEC);
+    vecpro[2] = (long)p2;
+    for (j=1; j<=m; j++) p2[j] = (long)prhinv;
     p2 += m;
     for (j=1; j<=n; j++)
-    {
-      p1[j] = matId[j];
       p2[j] = (long)idealmul(nf,(GEN)I[j],(GEN)alphalistinv[j]);
-    }
-    matH=nfhermite(nf,vecpro);
-    p1=algtobasis(nf,gmul(Aa,basistoalg(nf,(GEN)matH[1])));
-    p2=(GEN)matH[2];
 
-    tetpil=avma; neworder=cgetg(3,t_VEC);
-    H=cgetg(n+1,t_MAT); Hid=cgetg(n+1,t_VEC);
+    matH = nfhermite(nf,vecpro);
+    p1 = algtobasis(nf, gmul(Aa, basistoalg(nf,(GEN)matH[1])));
+    p2 = (GEN)matH[2];
+
+    neworder = cgetg(3,t_VEC);
+    H   = cgetg(n+1,t_MAT);
+    Hid = cgetg(n+1,t_VEC);
     for (j=1; j<=n; j++)
     {
+      alphainv = (GEN)alphalistinv[j];
+      if (alphainv == gun) { H[j] = p1[j]; Hid[j] = p2[j]; continue; }
+
       p3=cgetg(n+1,t_COL); H[j]=(long)p3;
       for (i=1; i<=n; i++)
-	p3[i]=(long)element_mul(nf,gcoeff(p1,i,j),(GEN)alphalistinv[j]);
-      Hid[j]=(long)idealmul(nf,(GEN)p2[j],(GEN)alphalist[j]);
+	p3[i] = (long)element_mul(nf,gcoeff(p1,i,j),alphainv);
+      Hid[j] = (long)idealmul(nf,(GEN)p2[j],(GEN)alphalist[j]);
     }
-    if (DEBUGLEVEL>3)
-      { fprintferr(" new order:\n"); outerr(H); outerr(Hid); }
+    if (DEBUGLEVEL>3) { fprintferr(" new order:\n"); outerr(H); outerr(Hid); }
     if (sep == 2 || gegal(I,Hid))
     {
-      neworder[1]=(long)H; neworder[2]=(long)Hid;
-      return gerepile(av,tetpil,neworder);
+      neworder[1] = (long)H;
+      neworder[2] = (long)Hid;
+      return gerepilecopy(av, neworder);
     }
 
-    A=H; I=Hid;
+    A = H; I = Hid;
     if (low_stack(lim, stack_lim(av1,1)) || (cmpt & 3) == 0)
     {
       GEN *gptr[2]; gptr[0]=&A; gptr[1]=&I;
@@ -2665,7 +3077,7 @@ rnfordmax(GEN nf, GEN pol, GEN pr, GEN unnf, GEN id, G
 static void
 check_pol(GEN x, long v)
 {
-  long i,tx, lx = lg(x);
+  long i,tx, lx = lgef(x);
   if (varn(x) != v)
     err(talker,"incorrect variable in rnf function");
   for (i=2; i<lx; i++)
@@ -2680,7 +3092,7 @@ GEN
 fix_relative_pol(GEN nf, GEN x, int chk_lead)
 {
   GEN xnf = (typ(nf) == t_POL)? nf: (GEN)nf[1];
-  long i, vnf = varn(xnf), lx = lg(x);
+  long i, vnf = varn(xnf), lx = lgef(x);
   if (typ(x) != t_POL || varn(x) >= vnf)
     err(talker,"incorrect polynomial in rnf function");
   x = dummycopy(x);
@@ -2703,14 +3115,19 @@ fix_relative_pol(GEN nf, GEN x, int chk_lead)
 static GEN
 rnfround2all(GEN nf, GEN pol, long all)
 {
-  long av=avma,tetpil,i,j,n,N,nbidp,vpol,*ep;
-  GEN p1,p2,p3,p4,polnf,list,unnf,id,matId,I,W,pseudo,y,discpol,d,D,sym;
+  gpmem_t av = avma;
+  long i,j,n,N,nbidp,vpol,*ep;
+  GEN A,p1,p2,nfT,list,id,I,W,pseudo,y,d,D,sym,unnf,disc;
 
-  nf=checknf(nf); polnf=(GEN)nf[1]; vpol=varn(pol);
+  nf=checknf(nf); nfT=(GEN)nf[1]; vpol=varn(pol);
   pol = fix_relative_pol(nf,pol,1);
-  N=degpol(polnf); n=degpol(pol); discpol=discsr(pol);
-  list=idealfactor(nf,discpol); ep=(long*)list[2]; list=(GEN)list[1];
-  nbidp=lg(list)-1; for(i=1;i<=nbidp;i++) ep[i]=itos((GEN)ep[i]);
+  N = degpol(nfT);
+  n = degpol(pol);
+  disc = discsr(pol); pol = lift(pol);
+  list = idealfactor(nf, disc);
+  ep  = (long*)list[2];
+  list= (GEN)  list[1];
+  nbidp=lg(list)-1; for(i=1;i<=nbidp;i++) ep[i] = itos((GEN)ep[i]);
   if (DEBUGLEVEL>1)
   {
     fprintferr("Ideals to consider:\n");
@@ -2718,51 +3135,55 @@ rnfround2all(GEN nf, GEN pol, long all)
       if (ep[i]>1) fprintferr("%Z^%ld\n",list[i],ep[i]);
     flusherr();
   }
-  id=idmat(N); unnf=gscalcol_i(gun,N);
-  matId=idmat_intern(n,unnf, gscalcol_i(gzero,N));
+  id = idmat(N); unnf = (GEN)id[1];
   pseudo = NULL;
   for (i=1; i<=nbidp; i++)
     if (ep[i] > 1)
     {
-      y=rnfordmax(nf,pol,(GEN)list[i],unnf,id,matId);
+      y = rnfordmax(nf, pol, (GEN)list[i], disc);
       pseudo = rnfjoinmodules(nf,pseudo,y);
     }
-  if (!pseudo)
+  if (pseudo)
   {
-    I=cgetg(n+1,t_VEC); for (i=1; i<=n; i++) I[i]=(long)id;
-    pseudo=cgetg(3,t_VEC); pseudo[1]=(long)matId; pseudo[2]=(long)I;
+    A = (GEN)pseudo[1];
+    I = (GEN)pseudo[2];
   }
-  W=gmodulcp(mat_to_vecpol(basistoalg(nf,(GEN)pseudo[1]),vpol),pol);
-  p2=cgetg(n+1,t_MAT); for (j=1; j<=n; j++) p2[j]=lgetg(n+1,t_COL);
-  sym=polsym(pol,n-1);
+  else
+  {
+    I = cgetg(n+1,t_VEC); for (i=1; i<=n; i++) I[i] = (long)id;
+    A = idmat_intern(n, unnf, zerocol(N));
+  }
+  W = mat_to_vecpol(lift_intern(basistoalg(nf,A)), vpol);
+  sym = polsym_gen(pol, NULL, n-1, nfT, NULL);
+  p2 = cgetg(n+1,t_MAT); for (j=1; j<=n; j++) p2[j] = lgetg(n+1,t_COL);
   for (j=1; j<=n; j++)
     for (i=j; i<=n; i++)
     {
-      p1 = lift_intern(gmul((GEN)W[i],(GEN)W[j]));
-      coeff(p2,j,i)=coeff(p2,i,j)=(long)quicktrace(p1,sym);
+      p1 = RXQX_mul((GEN)W[i],(GEN)W[j], nfT);
+      p1 = RXQX_rem(p1, pol, nfT);
+      p1 = simplify_i(quicktrace(p1,sym));
+      coeff(p2,j,i)=coeff(p2,i,j) = (long)p1;
     }
-  d = algtobasis_intern(nf,det(p2));
+  d = algtobasis_i(nf, det(p2));
 
-  I=(GEN)pseudo[2];
-  i=1; while (i<=n && gegal((GEN)I[i],id)) i++;
-  if (i>n) D=id;
+  i=1; while (i<=n && gegal((GEN)I[i], id)) i++;
+  if (i > n) D = id;
   else
   {
     D = (GEN)I[i];
     for (i++; i<=n; i++)
-      if (!gegal((GEN)I[i],id)) D = idealmul(nf,D,(GEN)I[i]);
+      if (!gegal((GEN)I[i], id)) D = idealmul(nf,D,(GEN)I[i]);
     D = idealpow(nf,D,gdeux);
   }
-  p4=gun; p3=auxdecomp(content(d),0);
-  for (i=1; i<lg(p3[1]); i++)
-    p4 = gmul(p4, gpuigs(gcoeff(p3,i,1), itos(gcoeff(p3,i,2))>>1));
-  p4 = gsqr(p4); tetpil=avma;
+  p2 = core2partial(content(d));
+  p2 = sqri((GEN)p2[2]);
+
   i = all? 2: 0;
-  p1=cgetg(3 + i,t_VEC);
-  if (i) { p1[1]=lcopy((GEN)pseudo[1]); p1[2]=lcopy(I); }
+  p1=cgetg(3 + i, t_VEC);
+  if (i) { p1[1] = lcopy(A); p1[2] = lcopy(I); }
   p1[1+i] = (long)idealmul(nf,D,d);
-  p1[2+i] = ldiv(d,p4);
-  return gerepile(av,tetpil,p1);
+  p1[2+i] = (long)gdiv(d, p2);
+  return gerepileupto(av,p1);
 }
 
 GEN
@@ -2777,54 +3198,61 @@ rnfdiscf(GEN nf, GEN pol)
   return rnfround2all(nf,pol,0);
 }
 
-/* given bnf as output by buchinit and a pseudo-basis of an order
- * in HNF [A,I] (or [A,I,D,d] it does not matter), tries to simplify the
- * HNF as much as possible. The resulting matrix will be upper triangular
- * but the diagonal coefficients will not be equal to 1. The ideals
- * are guaranteed to be integral and primitive.
- */
 GEN
+gen_if_principal(GEN bnf, GEN x)
+{
+  gpmem_t av = avma;
+  GEN z = isprincipalall(bnf,x, nf_GEN_IF_PRINCIPAL | nf_FORCE);
+  if (typ(z) == t_INT) { avma = av; return NULL; }
+  return z;
+}
+
+/* given bnf and a pseudo-basis of an order in HNF [A,I] (or [A,I,D,d] it
+ * does not matter), tries to simplify the HNF as much as possible. The
+ * resulting matrix will be upper triangular but the diagonal coefficients
+ * will not be equal to 1. The ideals are guaranteed to be integral and
+ * primitive. */
+GEN
 rnfsimplifybasis(GEN bnf, GEN order)
 {
-  long av=avma,tetpil,j,N,n;
+  gpmem_t av = avma;
+  long j, n, l;
   GEN p1,id,Az,Iz,nf,A,I;
 
-  bnf = checkbnf(bnf);
+  bnf = checkbnf(bnf); nf = (GEN)bnf[7];
   if (typ(order)!=t_VEC || lg(order)<3)
     err(talker,"not a pseudo-basis in nfsimplifybasis");
-  A=(GEN)order[1]; I=(GEN)order[2]; n=lg(A)-1; nf=(GEN)bnf[7];
-  N=degpol(nf[1]); id=idmat(N); Iz=cgetg(n+1,t_VEC); Az=cgetg(n+1,t_MAT);
+  A = (GEN)order[1];
+  I = (GEN)order[2]; n = lg(A)-1;
+  id = idmat(degpol(nf[1]));
+  Iz = cgetg(n+1,t_VEC);
+  Az = cgetg(n+1,t_MAT);
   for (j=1; j<=n; j++)
   {
-    if (gegal((GEN)I[j],id)) { Iz[j]=(long)id; Az[j]=A[j]; }
-    else
+    if (gegal((GEN)I[j],id)) { Iz[j]=(long)id; Az[j]=A[j]; continue; }
+
+    Iz[j] = (long)Q_primitive_part((GEN)I[j], &p1);
+    Az[j] = p1? lmul((GEN)A[j],p1): A[j];
+    if (p1 && gegal((GEN)Iz[j], id)) continue;
+
+    p1 = gen_if_principal(bnf, (GEN)Iz[j]);
+    if (p1)
     {
-      p1=content((GEN)I[j]);
-      if (!gcmp1(p1))
-      {
-	Iz[j]=(long)gdiv((GEN)I[j],p1); Az[j]=lmul((GEN)A[j],p1);
-      }
-      else Az[j]=A[j];
-      if (!gegal((GEN)Iz[j],id))
-      {
-	p1=isprincipalgen(bnf,(GEN)Iz[j]);
-	if (gcmp0((GEN)p1[1]))
-	{
-	  p1=(GEN)p1[2]; Iz[j]=(long)id;
-	  Az[j]=(long)element_mulvec(nf,p1,(GEN)Az[j]);
-	}
-      }
+      Iz[j] = (long)id;
+      Az[j] = (long)element_mulvec(nf,p1,(GEN)Az[j]);
     }
   }
-  tetpil=avma; p1=cgetg(lg(order),t_VEC); p1[1]=lcopy(Az); p1[2]=lcopy(Iz);
-  for (j=3; j<lg(order); j++) p1[j]=lcopy((GEN)order[j]);
-  return gerepile(av,tetpil,p1);
+  l = lg(order); p1 = cgetg(l, t_VEC);
+  p1[1] = (long)Az;
+  p1[2] = (long)Iz; for (j=3; j<l; j++) p1[j] = order[j];
+  return gerepilecopy(av, p1);
 }
 
 GEN
 rnfdet2(GEN nf, GEN A, GEN I)
 {
-  long av,tetpil,i;
+  gpmem_t av,tetpil;
+  long i;
   GEN p1;
 
   nf=checknf(nf); av = tetpil = avma;
@@ -2847,179 +3275,192 @@ rnfdet0(GEN nf, GEN x, GEN y)
   return y? rnfdet2(nf,x,y): rnfdet(nf,x);
 }
 
-/* given a pseudo-basis of an order in HNF [A,I] (or [A,I,D,d] it does
- * not matter), gives an nxn matrix (not in HNF) of a pseudo-basis and
- * an ideal vector [id,id,...,id,I] such that order=nf[7]^(n-1)xI.
- * Since it uses the approximation theorem, can be long.
- */
+/* 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. */
+static GEN
+nfidealdet1(GEN nf, GEN a, GEN b)
+{
+  gpmem_t av = avma;
+  GEN x,p1,res,u,v,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);
+  u = (GEN)p1[1];
+  v = (GEN)p1[2];
+
+  res = cgetg(5,t_VEC);
+  res[1] = (long)gmul(x,da);
+  res[2] = (long)gdiv(v,db);
+  res[3] = lnegi(db);
+  res[4] = (long)element_div(nf, u, (GEN)res[1]);
+  return gerepilecopy(av,res);
+}
+
+static GEN
+get_order(GEN nf, GEN O, char *s)
+{
+  if (typ(O) == t_POL)
+    return rnfpseudobasis(nf, O);
+  if (typ(O)!=t_VEC || lg(O) < 3)
+    err(talker,"not a pseudo-matrix in %s", s);
+  return O;
+}
+
+/* given a pseudo-basis of an order in HNF [A,I] (or [A,I,D,d]), gives an
+ * n x n matrix (not in HNF) of a pseudo-basis and an ideal vector
+ * [id,id,...,id,I] such that order = nf[7]^(n-1) x I.
+ * Uses the approximation theorem ==> slow. */
 GEN
 rnfsteinitz(GEN nf, GEN order)
 {
-  long av=avma,tetpil,i,n;
+  gpmem_t av = avma;
+  long i,n,l;
   GEN Id,A,I,p1,a,b;
 
   nf = checknf(nf);
   Id = idmat(degpol(nf[1]));
-  if (typ(order)==t_POL) order=rnfpseudobasis(nf,order);
-  if (typ(order)!=t_VEC || lg(order)<3)
-    err(talker,"not a pseudo-matrix in rnfsteinitz");
-  A=dummycopy((GEN)order[1]);
-  I=dummycopy((GEN)order[2]); n=lg(A)-1;
+  order = get_order(nf, order, "rnfsteinitz");
+  A = matalgtobasis(nf, (GEN)order[1]);
+  I = dummycopy((GEN)order[2]); n=lg(A)-1;
   if (typ(A) != t_MAT || typ(I) != t_VEC || lg(I) != n+1)
     err(typeer,"rnfsteinitz");
   for (i=1; i<n; i++)
   {
-    a = (GEN)I[i];
-    if (!gegal(a,Id))
+    GEN c1,c2;
+    a = (GEN)I[i]; if (gegal(a,Id)) continue;
+
+    c1 = (GEN)A[i];
+    c2 = (GEN)A[i+1];
+    b = (GEN)I[i+1];
+    if (gegal(b,Id))
     {
-      GEN c1 = (GEN)A[i];
-      GEN c2 = (GEN)A[i+1];
-      b = (GEN)I[i+1];
-      if (gegal(b,Id))
-      {
-        A[i]  = (long)c2;
-        A[i+1]= lneg(c1);
-	I[i]  = (long)b;
-        I[i+1]= (long)a;
-      }
-      else
-      {
-	p1 = nfidealdet1(nf,a,b);
-	A[i]  = ladd(element_mulvec(nf,(GEN)p1[1], c1),
-		     element_mulvec(nf,(GEN)p1[2], c2));
-	A[i+1]= ladd(element_mulvec(nf,(GEN)p1[3], c1),
-	             element_mulvec(nf,(GEN)p1[4], c2));
-	I[i]  =(long)Id;
-        I[i+1]=(long)idealmul(nf,a,b); p1 = content((GEN)I[i+1]);
-	if (!gcmp1(p1))
-	{
-	  I[i+1] = ldiv((GEN)I[i+1],p1);
-	  A[i+1] = lmul(p1,(GEN)A[i+1]);
-	}
-      }
+      A[i]  = (long)c2;
+      A[i+1]= lneg(c1);
+      I[i]  = (long)b;
+      I[i+1]= (long)a;
     }
+    else
+    {
+      p1 = nfidealdet1(nf,a,b);
+      A[i]  = ladd(element_mulvec(nf, (GEN)p1[1], c1),
+                   element_mulvec(nf, (GEN)p1[2], c2));
+      A[i+1]= ladd(element_mulvec(nf, (GEN)p1[3], c1),
+                   element_mulvec(nf, (GEN)p1[4], c2));
+      I[i]  = (long)Id;
+      I[i+1]= (long)Q_primitive_part(idealmul(nf,a,b), &p1);
+      if (p1) A[i+1] = (long)element_mulvec(nf, p1,(GEN)A[i+1]);
+    }
   }
-  tetpil=avma; p1=cgetg(lg(order),t_VEC);
-  p1[1]=lcopy(A); p1[2]=lcopy(I);
-  for (i=3; i<lg(order); i++) p1[i]=lcopy((GEN)order[i]);
-  return gerepile(av,tetpil,p1);
+  l = lg(order);
+  p1 = cgetg(l,t_VEC);
+  p1[1]=(long)A;
+  p1[2]=(long)I; for (i=3; i<l; i++) p1[i]=order[i];
+  return gerepilecopy(av, p1);
 }
 
-/* Given bnf as output by buchinit and either an order as output by
- * rnfpseudobasis or a polynomial, and outputs a basis if it is free,
- * an n+1-generating set if it is not
- */
+/* Given bnf and either an order as output by rnfpseudobasis or a polynomial,
+ * and outputs a basis if it is free, an n+1-generating set if it is not */
 GEN
 rnfbasis(GEN bnf, GEN order)
 {
-  ulong av = avma;
-  long j,N,n;
-  GEN nf,A,I,classe,p1,p2,id;
+  gpmem_t av = avma;
+  long j, n;
+  GEN nf, A, I, cl, col, a, id;
 
-  bnf = checkbnf(bnf);
-  nf=(GEN)bnf[7]; N=degpol(nf[1]); id=idmat(N);
-  if (typ(order)==t_POL) order=rnfpseudobasis(nf,order);
-  if (typ(order)!=t_VEC || lg(order)<3)
-    err(talker,"not a pseudo-matrix in rnfbasis");
-  A=(GEN)order[1]; I=(GEN)order[2]; n=lg(A)-1;
+  bnf = checkbnf(bnf); nf = (GEN)bnf[7];
+  id = idmat(degpol(nf[1]));
+  order = get_order(nf, order, "rnfbasis");
+  I = (GEN)order[2]; n = lg(I)-1;
   j=1; while (j<n && gegal((GEN)I[j],id)) j++;
-  if (j<n) order=rnfsteinitz(nf,order);
-  A=(GEN)order[1]; I=(GEN)order[2]; classe=(GEN)I[n];
-  p1=isprincipalgen(bnf,classe);
-  if (gcmp0((GEN)p1[1]))
+  if (j<n)
   {
-    p2=cgetg(n+1,t_MAT);
-    p2[n]=(long)element_mulvec(nf,(GEN)p1[2],(GEN)A[n]);
+    order = rnfsteinitz(nf,order);
+    I = (GEN)order[2];
   }
-  else
+  A = (GEN)order[1];
+  col= (GEN)A[n]; A = vecextract_i(A, 1, n-1);
+  cl = (GEN)I[n];
+  a = gen_if_principal(bnf, cl);
+  if (!a)
   {
-    p1=ideal_two_elt(nf,classe);
-    p2=cgetg(n+2,t_MAT);
-    p2[n]=lmul((GEN)p1[1],(GEN)A[n]);
-    p2[n+1]=(long)element_mulvec(nf,(GEN)p1[2],(GEN)A[n]);
+    GEN p1 = ideal_two_elt(nf, cl);
+    A = concatsp(A, gmul((GEN)p1[1], col));
+    a = (GEN)p1[2];
   }
-  for (j=1; j<n; j++) p2[j]=A[j];
-  return gerepilecopy(av,p2);
+  A = concatsp(A, element_mulvec(nf, a, col));
+  return gerepilecopy(av, A);
 }
 
-/* Given bnf as output by buchinit and either an order as output by
- * rnfpseudobasis or a polynomial, and outputs a basis (not pseudo)
- * in Hermite Normal Form if it exists, zero if not
+/* Given bnf and either an order as output by rnfpseudobasis or a polynomial,
+ * and outputs a basis (not pseudo) in Hermite Normal Form if it exists, zero
+ * if not
  */
 GEN
 rnfhermitebasis(GEN bnf, GEN order)
 {
-  ulong av = avma;
-  long j,N,n;
-  GEN nf,A,I,p1,id;
+  gpmem_t av = avma;
+  long j, n;
+  GEN nf, A, I, a, id;
 
-  bnf = checkbnf(bnf); nf=(GEN)bnf[7];
-  N=degpol(nf[1]); id=idmat(N);
-  if (typ(order)==t_POL)
-  {
-    order=rnfpseudobasis(nf,order);
-    A=(GEN)order[1];
-  }
-  else
-  {
-    if (typ(order)!=t_VEC || lg(order)<3)
-      err(talker,"not a pseudo-matrix in rnfbasis");
-    A=gcopy((GEN)order[1]);
-  }
-  I=(GEN)order[2]; n=lg(A)-1;
+  bnf = checkbnf(bnf); nf = (GEN)bnf[7];
+  id = idmat(degpol(nf[1]));
+  order = get_order(nf, order, "rnfbasis");
+  A = (GEN)order[1]; A = dummycopy(A);
+  I = (GEN)order[2]; n = lg(A)-1;
   for (j=1; j<=n; j++)
   {
-    if (!gegal((GEN)I[j],id))
-    {
-      p1=isprincipalgen(bnf,(GEN)I[j]);
-      if (gcmp0((GEN)p1[1]))
-	A[j]=(long)element_mulvec(nf,(GEN)p1[2],(GEN)A[j]);
-      else { avma=av; return gzero; }
-    }
+    if (gegal((GEN)I[j],id)) continue;
+
+    a = gen_if_principal(bnf, (GEN)I[j]);
+    if (!a) { avma = av; return gzero; }
+    A[j] = (long)element_mulvec(nf, a, (GEN)A[j]);
   }
   return gerepilecopy(av,A);
 }
 
-long
-rnfisfree(GEN bnf, GEN order)
+static long
+_rnfisfree(GEN bnf, GEN order)
 {
-  long av=avma,n,N,j;
-  GEN nf,p1,id,I;
+  long n, j;
+  GEN nf, p1, id, I;
 
   bnf = checkbnf(bnf);
   if (gcmp1(gmael3(bnf,8,1,1))) return 1;
 
-  nf=(GEN)bnf[7]; N=degpol(nf[1]); id=idmat(N);
-  if (typ(order)==t_POL) order=rnfpseudobasis(nf,order);
-  if (typ(order)!=t_VEC || lg(order)<3)
-    err(talker,"not a pseudo-matrix in rnfisfree");
-
-  I=(GEN)order[2]; n=lg(I)-1;
+  nf = (GEN)bnf[7]; id = idmat(degpol(nf[1]));
+  order = get_order(nf, order, "rnfisfree");
+  I = (GEN)order[2]; n = lg(I)-1;
   j=1; while (j<=n && gegal((GEN)I[j],id)) j++;
-  if (j>n) { avma=av; return 1; }
+  if (j>n) return 1;
 
-  p1=(GEN)I[j];
+  p1 = (GEN)I[j];
   for (j++; j<=n; j++)
-    if (!gegal((GEN)I[j],id)) p1=idealmul(nf,p1,(GEN)I[j]);
-  j = gcmp0(isprincipal(bnf,p1));
-  avma=av; return j;
+    if (!gegal((GEN)I[j],id)) p1 = idealmul(nf,p1,(GEN)I[j]);
+  return gcmp0( isprincipal(bnf,p1) );
 }
 
+long
+rnfisfree(GEN bnf, GEN order)
+{
+  gpmem_t av = avma;
+  long n = _rnfisfree(bnf, order);
+  avma = av; return n;
+}
+
 /**********************************************************************/
 /**								     **/
 /**		      COMPOSITUM OF TWO NUMBER FIELDS                **/
 /**								     **/
 /**********************************************************************/
-extern GEN ZY_ZXY_resultant_all(GEN A, GEN B0, long *lambda, GEN *LPRS);
-extern GEN squff2(GEN x, long klim, long hint);
-extern GEN to_polmod(GEN x, GEN mod);
-
-/* modular version. TODO: check that compositum2 is not slower */
+/* modular version */
 GEN
 polcompositum0(GEN A, GEN B, long flall)
 {
-  ulong av = avma;
-  long v,k;
+  gpmem_t av = avma;
+  long v, k;
   GEN C, LPRS;
 
   if (typ(A)!=t_POL || typ(B)!=t_POL) err(typeer,"polcompositum0");
@@ -3034,18 +3475,19 @@ polcompositum0(GEN A, GEN B, long flall)
   if (!ZX_is_squarefree(B)) err(talker,"compositum: %Z not separable", B);
 
   k = 1; C = ZY_ZXY_resultant_all(A, B, &k, flall? &LPRS: NULL);
-  C = squff2(C,0,0); /* C = Res_Y (A, B(X + kY)) guaranteed squarefree */
+  C = DDF2(C, 0); /* C = Res_Y (A, B(X + kY)) guaranteed squarefree */
+  C = sort_vecpol(C);
   if (flall)
   {
     long i,l = lg(C);
     GEN w,a,b; /* a,b,c root of A,B,C = compositum, c = b - k a */
     for (i=1; i<l; i++)
     { /* invmod possibly very costly */
-      a = gmul((GEN)LPRS[1], ZX_invmod((GEN)LPRS[2], (GEN)C[i]));
+      a = gmul((GEN)LPRS[1], QX_invmod((GEN)LPRS[2], (GEN)C[i]));
       a = gneg_i(gmod(a, (GEN)C[i]));
       b = gadd(polx[v], gmulsg(k,a));
       w = cgetg(5,t_VEC); /* [C, a, b, n ] */
-      w[1] = C[i]; 
+      w[1] = C[i];
       w[2] = (long)to_polmod(a, (GEN)w[1]);
       w[3] = (long)to_polmod(b, (GEN)w[1]);
       w[4] = lstoi(-k); C[i] = (long)w;
@@ -3066,13 +3508,10 @@ compositum2(GEN pol1,GEN pol2)
   return polcompositum0(pol1,pol2,1);
 }
 
-extern int isrational(GEN x);
-extern GEN nfgcd(GEN P, GEN Q, GEN nf, GEN den);
-
 int
 nfissquarefree(GEN nf, GEN x)
 {
-  ulong av = avma;
+  gpmem_t av = avma;
   GEN g, y = derivpol(x);
   if (isrational(x))
     g = modulargcd(x, y);
@@ -3082,40 +3521,50 @@ nfissquarefree(GEN nf, GEN x)
 }
 
 GEN
-rnfequation0(GEN nf, GEN B, long flall)
+_rnfequation(GEN A, GEN B, long *pk, GEN *pLPRS)
 {
-  ulong av = avma;
-  long v,vpol,k,lA,lB;
-  GEN cC,A,C,LPRS;
+  long k, lA, lB;
+  GEN nf, C;
 
-  if (typ(nf)==t_POL) A=nf; else { nf=checknf(nf); A=(GEN)nf[1]; }
-  B = fix_relative_pol(nf,B,1);
-  v   = varn(A); lA = lgef(A);
-  vpol= varn(B); lB = lgef(B);
+  A = get_nfpol(A, &nf);
+  B = fix_relative_pol(A,B,1);
+  lA = lgef(A);
+  lB = lgef(B);
   if (lA<=3 || lB<=3) err(constpoler,"rnfequation");
 
   check_pol_int(A,"rnfequation");
-  B = lift_intern(B); B = gdiv(B, content(B));
+  B = primpart(lift_intern(B));
   for (k=2; k<lB; k++)
     if (lgef(B[k]) >= lA) B[k] = lres((GEN)B[k],A);
 
   if (!nfissquarefree(A,B))
     err(talker,"not k separable relative equation in rnfequation");
 
-  k = 0; C = ZY_ZXY_resultant_all(A, B, &k, flall? &LPRS: NULL);
+  *pk = 0; C = ZY_ZXY_resultant_all(A, B, pk, pLPRS);
   if (gsigne(leadingcoeff(C)) < 0) C = gneg_i(C);
-  C = primitive_part(C, &cC);
+  *pk = -*pk; return primpart(C);
+}
+
+GEN
+rnfequation0(GEN A, GEN B, long flall)
+{
+  gpmem_t av = avma;
+  long k;
+  GEN LPRS, nf, C;
+
+  A = get_nfpol(A, &nf);
+  C = _rnfequation(A, B, &k, flall? &LPRS: NULL);
   if (flall)
   {
-    GEN w,a,b; /* a,b,c root of A,B,C = compositum, c = b - k a */
+    GEN w,a,b; /* a,b,c root of A,B,C = compositum, c = b + k a */
     /* invmod possibly very costly */
-    a = gmul((GEN)LPRS[1], ZX_invmod((GEN)LPRS[2], C));
+    a = gmul((GEN)LPRS[1], QX_invmod((GEN)LPRS[2], C));
     a = gneg_i(gmod(a, C));
-    b = gadd(polx[v], gmulsg(k,a));
+    b = gadd(polx[varn(A)], gmulsg(k,a));
     w = cgetg(4,t_VEC); /* [C, a, n ] */
-    w[1] = (long)C; 
+    w[1] = (long)C;
     w[2] = (long)to_polmod(a, (GEN)w[1]);
-    w[3] = lstoi(-k); C = w;
+    w[3] = lstoi(k); C = w;
   }
   return gerepilecopy(av, C);
 }
@@ -3186,9 +3635,8 @@ static GEN
 rnfdiv(GEN x, GEN y)
 {
   long i, ru = lg(x);
-  GEN z;
+  GEN z = cgetg(ru,t_COL);
 
-  z=cgetg(ru,t_COL);
   for (i=1; i<ru; i++) z[i]=(long)gdiv((GEN)x[i],(GEN)y[i]);
   return z;
 }
@@ -3197,9 +3645,8 @@ static GEN
 rnfmul(GEN x, GEN y)
 {
   long i, ru = lg(x);
-  GEN z;
+  GEN z = cgetg(ru,t_COL);
 
-  z=cgetg(ru,t_COL);
   for (i=1; i<ru; i++) z[i]=(long)gmul((GEN)x[i],(GEN)y[i]);
   return z;
 }
@@ -3208,47 +3655,46 @@ static GEN
 rnfvecmul(GEN x, GEN v)
 {
   long i, lx = lg(v);
-  GEN y;
+  GEN y = cgetg(lx,typ(v));
 
-  y=cgetg(lx,typ(v));
   for (i=1; i<lx; i++) y[i]=(long)rnfmul(x,(GEN)v[i]);
   return y;
 }
 
 static GEN
-allonge(GEN v, long N)
+allonge(GEN v, long l)
 {
-  long r,r2,i;
-  GEN y;
+  long r, r2, i;
+  GEN y = cgetg(l,t_COL);
 
-  r=lg(v)-1; r2=N-r;
-  y=cgetg(N+1,t_COL);
-  for (i=1; i<=r; i++) y[i]=v[i];
-  for ( ; i<=N; i++) y[i]=(long)gconj((GEN)v[i-r2]);
+  r = lg(v); r2 = l-r;
+  for (i=1; i < r; i++) y[i] = v[i];
+  for (   ; i < l; i++) y[i] = lconj((GEN)v[i-r2]);
   return y;
 }
 
 static GEN
 findmin(GEN nf, GEN ideal, GEN muf,long prec)
 {
-  long av=avma,N,tetpil,i;
-  GEN m,y;
+  gpmem_t av = avma;
+  long i, l;
+  GEN m,y, G = gmael(nf,5,2);
 
-  m = qf_base_change(gmael(nf,5,3), ideal, 0); /* nf[5][3] = T2 */
-  m = lllgramintern(m,4,1,prec);
+  m = gmul(G, ideal);
+  m = lllintern(m,4,1,prec);
   if (!m)
   {
-    m = lllint(ideal);
-    m = qf_base_change(gmael(nf,5,3), gmul(ideal,m), 0);
-    m = lllgramintern(m,4,1,prec);
+    ideal = lllint_ip(ideal,4);
+    m = gmul(G, ideal);
+    m = lllintern(m,4,1,prec);
     if (!m) err(precer,"rnflllgram");
   }
-  ideal=gmul(ideal,m);
-  N=lg(ideal)-1; y=cgetg(N+1,t_MAT);
-  for (i=1; i<=N; i++)
-    y[i] = (long) allonge(nftocomplex(nf,(GEN)ideal[i]),N);
-  m=ground(greal(gauss(y,allonge(muf,N))));
-  tetpil=avma; return gerepile(av,tetpil,gmul(ideal,m));
+  ideal = gmul(ideal,m);
+  l = lg(ideal); y = cgetg(l,t_MAT);
+  for (i=1; i<l; i++)
+    y[i] = (long)allonge(nftocomplex(nf,(GEN)ideal[i]),l);
+  m = ground(greal(gauss(y, allonge(muf,l))));
+  return gerepileupto(av, gmul(ideal,m));
 }
 
 #define swap(x,y) { long _t=x; x=y; y=_t; }
@@ -3260,7 +3706,8 @@ findmin(GEN nf, GEN ideal, GEN muf,long prec)
 GEN
 rnflllgram(GEN nf, GEN pol, GEN order,long prec)
 {
-  long av=avma,tetpil,i,j,k,l,kk,kmax,r1,ru,lx,vnf;
+  gpmem_t av=avma,tetpil;
+  long i,j,k,l,kk,kmax,r1,ru,lx,vnf;
   GEN p1,p2,M,I,U,ronf,poll,unro,roorder,powreorder,mth,s,MC,MPOL,MCS;
   GEN B,mu,Bf,temp,ideal,x,xc,xpol,muf,mufc,muno,y,z,Ikk_inv;
 
@@ -3431,117 +3878,174 @@ rnflllgram(GEN nf, GEN pol, GEN order,long prec)
 GEN
 rnfpolred(GEN nf, GEN pol, long prec)
 {
-  ulong av = avma;
-  long i,j,k,n,N, vpol = varn(pol);
-  GEN id,id2,newid,newor,p1,p2,al,newpol,w,z;
-  GEN bnf,zk,newideals,ideals,order,neworder;
+  gpmem_t av = avma;
+  long i, j, n, v = varn(pol);
+  GEN id, al, w, I, O, bnf;
 
   if (typ(pol)!=t_POL) err(typeer,"rnfpolred");
-  if (typ(nf)!=t_VEC) err(idealer1);
-  switch(lg(nf))
-  {
-    case 10: bnf = NULL; break;
-    case 11: bnf = nf; nf = checknf((GEN)nf[7]); break;
-    default: err(idealer1);
-      return NULL; /* not reached */
-  }
-  if (degpol(pol) <= 1) 
-  {
-    w=cgetg(2,t_VEC);
-    w[1]=lpolx[vpol]; return w;
-  }
-  id=rnfpseudobasis(nf,pol); N=degpol(nf[1]);
+  bnf = nf; nf = checknf(bnf);
+  bnf = (nf == bnf)? NULL: checkbnf(bnf);
+  if (degpol(pol) <= 1) return _vec(polx[v]);
+
+  id = rnfpseudobasis(nf,pol);
   if (bnf && gcmp1(gmael3(bnf,8,1,1))) /* if bnf is principal */
   {
-    ideals=(GEN)id[2]; n=lg(ideals)-1; order=(GEN)id[1];
-    newideals=cgetg(n+1,t_VEC); neworder=cgetg(n+1,t_MAT);
-    zk=idmat(N);
+    GEN newI, newO, zk = idmat(degpol(nf[1]));
+    O = (GEN)id[1];
+    I = (GEN)id[2]; n = lg(I)-1;
+    newI = cgetg(n+1,t_VEC);
+    newO = cgetg(n+1,t_MAT);
     for (j=1; j<=n; j++)
     {
-      newideals[j]=(long)zk; p1=cgetg(n+1,t_COL); neworder[j]=(long)p1;
-      p2=(GEN)order[j];
-      al=(GEN)isprincipalgen(bnf,(GEN)ideals[j])[2];
-      for (k=1; k<=n; k++)
-	p1[k]=(long)element_mul(nf,(GEN)p2[k],al);
+      newI[j] = (long)zk; al = gen_if_principal(bnf,(GEN)I[j]);
+      newO[j] = (long)element_mulvec(nf, al, (GEN)O[j]);
     }
-    id=cgetg(3,t_VEC); id[1]=(long)neworder; id[2]=(long)newideals;
+    id = cgetg(3,t_VEC);
+    id[1] = (long)newO;
+    id[2] = (long)newI;
   }
-  id2=rnflllgram(nf,pol,id,prec);
-  z=(GEN)id2[1]; newid=(GEN)z[2]; newor=(GEN)z[1];
-  n=lg(newor)-1; w=cgetg(n+1,t_VEC);
+  
+  id = (GEN)rnflllgram(nf,pol,id,prec)[1];
+  O = (GEN)id[1];
+  I = (GEN)id[2]; n = lg(I)-1;
+  w = cgetg(n+1,t_VEC);
+  pol = lift(pol);
   for (j=1; j<=n; j++)
   {
-    p1=(GEN)newid[j]; al=gmul(gcoeff(p1,1,1),(GEN)newor[j]);
-    p1=basistoalg(nf,(GEN)al[n]);
+    GEN p1, newpol;
+
+    p1 = (GEN)I[j]; al = gmul(gcoeff(p1,1,1),(GEN)O[j]);
+    p1 = basistoalg(nf,(GEN)al[n]);
     for (i=n-1; i; i--)
-      p1=gadd(basistoalg(nf,(GEN)al[i]),gmul(polx[vpol],p1));
-    newpol=gtopoly(gmodulcp(gtovec(caract2(lift(pol),lift(p1),vpol)),
-                            (GEN) nf[1]), vpol);
+      p1 = gadd(basistoalg(nf,(GEN)al[i]), gmul(polx[v],p1));
+    p1 = gmodulcp(gtovec(caract2(pol,lift(p1),v)), (GEN)nf[1]);
+    newpol = gtopoly(p1, v);
+
     p1 = ggcd(newpol, derivpol(newpol));
-    if (degpol(p1)>0)
+    if (degpol(p1) > 0)
     {
-      newpol=gdiv(newpol,p1);
-      newpol=gdiv(newpol,leading_term(newpol));
+      newpol = gdiv(newpol, p1);
+      newpol = gdiv(newpol, leading_term(newpol));
     }
-    w[j]=(long)newpol;
-    if (DEBUGLEVEL>=4) outerr(newpol);
+    w[j] = (long)newpol;
   }
   return gerepilecopy(av,w);
 }
 
-extern GEN vecpol_to_mat(GEN v, long n);
-
 /* given a relative polynomial pol over nf, compute a pseudo-basis for the
  * extension, then an absolute basis */
-GEN
-makebasis(GEN nf,GEN pol)
+static GEN
+makebasis(GEN nf, GEN pol, GEN rnfeq)
 {
-  GEN elts,ids,polabs,plg,B,bs,p1,p2,a,den,vbs,vbspro,vpro,rnf;
-  long av=avma,n,N,m,i,j, v = varn(pol);
+  GEN elts,ids,polabs,plg,plg0,B,bs,p1,den,vbs,d,vpro;
+  gpmem_t av = avma;
+  long n,N,m,i,j,k, v = varn(pol);
 
-  p1 = rnfequation2(nf,pol);
-  polabs= (GEN)p1[1];
-  plg   = (GEN)p1[2];
-  a     = (GEN)p1[3];
-  rnf = cgetg(12,t_VEC);
-  for (i=2;i<=9;i++) rnf[i]=zero;
-  rnf[1] =(long)pol;
-  rnf[10]=(long)nf; p2=cgetg(4,t_VEC);
-  rnf[11]=(long)p2; p2[1]=p2[2]=zero; p2[3]=(long)a;
-  if (signe(a))
-    pol = gsubst(pol,v,gsub(polx[v],
-                            gmul(a,gmodulcp(polx[varn(nf[1])],(GEN)nf[1]))));
-  p1=rnfpseudobasis(nf,pol);
+  polabs= (GEN)rnfeq[1];
+  plg   = (GEN)rnfeq[2]; plg = lift_intern(plg);
+  p1 = rnfpseudobasis(nf,pol);
   elts= (GEN)p1[1];
   ids = (GEN)p1[2];
-  if (DEBUGLEVEL>1) { fprintferr("relative basis computed\n"); flusherr(); }
-  N=degpol(pol); n=degpol((GEN)nf[1]); m=n*N;
-  den = denom(content(lift(plg)));
-  vbs = cgetg(n+1,t_VEC);
-  vbs[1] = un; 
-  vbs[2] = (long)plg; vbspro = gmul(den,plg);
-  for(i=3;i<=n;i++)
-    vbs[i] = ldiv(gmul((GEN)vbs[i-1],vbspro),den);
+  if (DEBUGLEVEL>1) fprintferr("relative basis computed\n");
+  N = degpol(pol);
+  n = degpol(nf[1]); m = n*N;
+
+  plg0= Q_remove_denom(plg, &den); /* plg = plg0/den */
+  /* nf = K = Q(a), vbs[i+1] = a^i as an elt of L = Q[X] / polabs */
+  vbs = RXQ_powers(plg0, polabs, n-1);
+  if (den)
+  { /* restore denominators */
+    vbs[2] = (long)plg; d = den;
+    for (i=3; i<=n; i++) { d = mulii(d,den); vbs[i] = ldiv((GEN)vbs[i], d); }
+  }
+
+  /* bs = integer basis of K, as elements of L */
   bs = gmul(vbs, vecpol_to_mat((GEN)nf[7],n));
 
-  vpro=cgetg(N+1,t_VEC);
-  for (i=1;i<=N;i++)
+  vpro = cgetg(N+1,t_VEC);
+  for (i=1; i<=N; i++) vpro[i] = lpowgs(polx[v], i-1);
+  vpro = gmul(vpro, elts);
+  B = cgetg(m+1, t_MAT);
+  for(i=k=1; i<=N; i++)
   {
-    p1=cgetg(3,t_POLMOD);
-    p1[1]=(long)polabs;
-    p1[2]=lpuigs(polx[v],i-1); vpro[i]=(long)p1;
+    GEN w = element_mulvec(nf, (GEN)vpro[i], (GEN)ids[i]);
+    for(j=1; j<=n; j++)
+    {
+      p1 = gres(gmul(bs, (GEN)w[j]), polabs);
+      B[k++] = (long)pol_to_vec(p1, m);
+    }
   }
-  vpro=gmul(vpro,elts); B = cgetg(m+1, t_MAT);
-  for(i=1;i<=N;i++)
-    for(j=1;j<=n;j++)
+  B = Q_remove_denom(B, &den);
+  if (den) { B = hnfmodid(B, den); B = gdiv(B, den); } else B = idmat(m);
+  p1 = cgetg(3,t_VEC);
+  p1[1] = (long)polabs;
+  p1[2] = (long)B; return gerepilecopy(av, p1);
+}
+
+/* relative polredabs. Returns relative polynomial by default (flag = 0)
+ * flag & nf_ORIG: + element (base change)
+ * flag & nf_ADDZK: + integer basis
+ * flag & nf_ABSOLUTE: absolute polynomial */
+GEN
+rnfpolredabs(GEN nf, GEN relpol, long flag)
+{
+  GEN red, bas, z, elt, POL, pol, T, a;
+  long v, fl;
+  gpmem_t av = avma;
+
+  if (typ(relpol)!=t_POL) err(typeer,"rnfpolredabs");
+  nf = checknf(nf); v = varn(relpol);
+  if (DEBUGLEVEL>1) (void)timer2();
+  relpol = unifpol(nf,relpol,1);
+  T = (GEN)nf[1];
+  if ((flag & nf_ADDZK) && (flag != (nf_ADDZK|nf_ABSOLUTE)))
+    err(impl,"this combination of flags in rnfpolredabs");
+  fl = (flag & nf_ADDZK)? nf_ORIG: nf_RAW;
+  if (flag & nf_PARTIALFACT)
+  {
+    long sa;
+    bas = NULL; /* -Wall */
+    POL = _rnfequation(nf, relpol, &sa, NULL);
+    red = polredabs0(POL, fl | nf_PARTIALFACT);
+    a = stoi(sa);
+  }
+  else
+  {
+    GEN rel, eq = rnfequation2(nf,relpol);
+    POL = (GEN)eq[1];
+    a   = (GEN)eq[3];
+    rel = poleval(relpol, gsub(polx[v],
+                               gmul(a, gmodulcp(polx[varn(T)],T))));
+    bas = makebasis(nf, rel, eq);
+    if (DEBUGLEVEL>1)
     {
-      p1 = gmul(bs, element_mul(nf,(GEN)vpro[i],gmael(ids,i,j)));
-      B[(i-1)*n+j] = (long)pol_to_vec(lift_intern(p1), m);
+      msgtimer("absolute basis");
+      fprintferr("original absolute generator: %Z\n", POL);
     }
-  p1 = denom(B); B = gmul(B,p1);
-  B = hnfmodid(B, p1); B = gdiv(B,p1);
-  p1=cgetg(4,t_VEC);
-  p1[1]=(long)polabs;
-  p1[2]=(long)B;
-  p1[3]=(long)rnf; return gerepilecopy(av, p1);
+    red = polredabs0(bas, fl);
+  }
+  pol = (GEN)red[1];
+  if (DEBUGLEVEL>1) fprintferr("reduced absolute generator: %Z\n",pol);
+  if (flag & nf_ABSOLUTE)
+  {
+    if (flag & nf_ADDZK)
+    {
+      GEN t = (GEN)red[2], B = (GEN)bas[2];
+      GEN v = RXQ_powers(lift_intern(t), pol, degpol(pol)-1);
+      z = cgetg(3, t_VEC);
+      z[1] = (long)pol;
+      z[2] = lmul(v, B);
+      return gerepilecopy(av, z);
+    }
+    return gerepilecopy(av, pol);
+  }
+
+  elt = eltabstorel((GEN)red[2], T, relpol, a);
+
+  pol = rnfcharpoly(nf,relpol,elt,v);
+  if (!(flag & nf_ORIG)) return gerepileupto(av, pol);
+  z = cgetg(3,t_VEC);
+  z[1] = (long)pol;
+  z[2] = (long)to_polmod(modreverse_i((GEN)elt[2],(GEN)elt[1]), pol);
+  return gerepilecopy(av, z);
 }