===================================================================
RCS file: /home/cvs/OpenXM_contrib/pari-2.2/src/basemath/Attic/polarit1.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/polarit1.c	2001/10/02 11:17:04	1.1
+++ OpenXM_contrib/pari-2.2/src/basemath/Attic/polarit1.c	2002/09/11 07:26:51	1.2
@@ -1,4 +1,4 @@
-/* $Id: polarit1.c,v 1.1 2001/10/02 11:17:04 noro Exp $
+/* $Id: polarit1.c,v 1.2 2002/09/11 07:26:51 noro Exp $
 
 Copyright (C) 2000  The PARI group.
 
@@ -21,14 +21,17 @@ Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 /***********************************************************************/
 #include "pari.h"
 extern GEN get_bas_den(GEN bas);
-extern GEN get_mul_table(GEN x,GEN bas,GEN invbas,GEN *T);
+extern GEN get_mul_table(GEN x,GEN bas,GEN invbas);
 extern GEN pol_to_monic(GEN pol, GEN *lead);
 
+#define lswap(x,y) { long _t=x; x=y; y=_t; }
+#define swap(x,y) { GEN _t=x; x=y; y=_t; }
+
 /* see splitgen() for how to use these two */
 GEN
 setloop(GEN a)
 {
-  a=icopy(a); new_chunk(2); /* dummy to get two cells of extra space */
+  a=icopy(a); (void)new_chunk(2); /* dummy to get two cells of extra space */
   return a;
 }
 
@@ -86,14 +89,14 @@ incloop(GEN a)
 int
 gdivise(GEN x, GEN y)
 {
-  long av=avma;
+  gpmem_t av=avma;
   x=gmod(x,y); avma=av; return gcmp0(x);
 }
 
 int
 poldivis(GEN x, GEN y, GEN *z)
 {
-  long av = avma;
+  gpmem_t av = avma;
   GEN p1 = poldivres(x,y,ONLY_DIVIDES);
   if (p1) { *z = p1; return 1; }
   avma=av; return 0;
@@ -108,21 +111,16 @@ poldivis(GEN x, GEN y, GEN *z)
  *   if z = ONLY_REM  return remainder, otherwise return quotient
  *   if z != NULL set *z to remainder
  *   *z is the last object on stack (and thus can be disposed of with cgiv
- *   instead of gerepile)
- */
+ *   instead of gerepile) */
 GEN
 poldivres(GEN x, GEN y, GEN *pr)
 {
-  ulong avy,av,av1;
+  gpmem_t avy, av, av1;
   long ty=typ(y),tx,vx,vy,dx,dy,dz,i,j,sx,lrem;
-  int remainder;
   GEN z,p1,rem,y_lead,mod;
   GEN (*f)(GEN,GEN);
 
-  if (pr == ONLY_DIVIDES_EXACT)
-    { f = gdivexact; pr = ONLY_DIVIDES; } 
-  else
-    f = gdiv;
+  f = gdiv;
   if (is_scalar_t(ty))
   {
     if (pr == ONLY_REM) return gzero;
@@ -149,9 +147,10 @@ poldivres(GEN x, GEN y, GEN *pr)
     }
     return f(x,y);
   }
-  if (!signe(y)) err(talker,"euclidean division by zero (poldivres)");
 
-  dy=degpol(y); y_lead = (GEN)y[dy+2];
+  if (!signe(y)) err(talker,"division by zero in poldivrem");
+  dy = degpol(y);
+  y_lead = (GEN)y[dy+2];
   if (gcmp0(y_lead)) /* normalize denominator if leading term is 0 */
   {
     err(warner,"normalizing a polynomial with 0 leading term");
@@ -170,7 +169,7 @@ poldivres(GEN x, GEN y, GEN *pr)
     }
     return f(x, constant_term(y));
   }
-  dx=degpol(x);
+  dx = degpol(x);
   if (vx>vy || dx<dy)
   {
     if (pr)
@@ -181,11 +180,14 @@ poldivres(GEN x, GEN y, GEN *pr)
     }
     return zeropol(vy);
   }
-  dz=dx-dy; av=avma; /* to avoid gsub's later on */
+
+  /* x,y in R[X], y non constant */
+  dz = dx-dy; av = avma;
   p1 = new_chunk(dy+3);
   for (i=2; i<dy+3; i++)
   {
     GEN p2 = (GEN)y[i];
+    /* gneg to avoid gsub's later on */
     p1[i] = isexactzero(p2)? 0: (long)gneg_i(p2);
   }
   y = p1;
@@ -198,30 +200,34 @@ poldivres(GEN x, GEN y, GEN *pr)
     default: if (gcmp1(y_lead)) y_lead = NULL;
       mod = NULL;
   }
-  avy=avma; z=cgetg(dz+3,t_POL);
-  z[1]=evalsigne(1) | evallgef(dz+3) | evalvarn(vx);
+  avy=avma;
+  z = cgetg(dz+3,t_POL);
+  z[1] = evalsigne(1) | evallgef(dz+3) | evalvarn(vx);
   x += 2; y += 2; z += 2;
 
-  p1 = (GEN)x[dx]; remainder = (pr == ONLY_REM);
-  z[dz]=y_lead? (long)f(p1,y_lead): lcopy(p1);
+  p1 = (GEN)x[dx];
+  z[dz] = y_lead? (long)f(p1,y_lead): lcopy(p1);
   for (i=dx-1; i>=dy; i--)
   {
     av1=avma; p1=(GEN)x[i];
     for (j=i-dy+1; j<=i && j<=dz; j++)
-      if (y[i-j]) p1 = gadd(p1, gmul((GEN)z[j],(GEN)y[i-j]));
+      if (y[i-j] && z[j] != zero) p1 = gadd(p1, gmul((GEN)z[j],(GEN)y[i-j]));
     if (y_lead) p1 = f(p1,y_lead);
-    if (!remainder) p1 = avma==av1? gcopy(p1): gerepileupto(av1,p1);
+
+    if (isexactzero(p1)) { avma=av1; p1 = gzero; }
+    else
+      p1 = avma==av1? gcopy(p1): gerepileupto(av1,p1);
     z[i-dy] = (long)p1;
   }
   if (!pr) return gerepileupto(av,z-2);
 
-  rem = (GEN)avma; av1 = (long)new_chunk(dx+3);
+  rem = (GEN)avma; av1 = (gpmem_t)new_chunk(dx+3);
   for (sx=0; ; i--)
   {
     p1 = (GEN)x[i];
     /* we always enter this loop at least once */
     for (j=0; j<=i && j<=dz; j++)
-      if (y[i-j]) p1 = gadd(p1, gmul((GEN)z[j],(GEN)y[i-j]));
+      if (y[i-j] && z[j] != zero) p1 = gadd(p1, gmul((GEN)z[j],(GEN)y[i-j]));
     if (mod && avma==av1) p1 = gmul(p1,mod);
     if (!gcmp0(p1)) { sx = 1; break; } /* remainder is non-zero */
     if (!isinexactreal(p1) && !isexactzero(p1)) break;
@@ -231,12 +237,12 @@ poldivres(GEN x, GEN y, GEN *pr)
   if (pr == ONLY_DIVIDES)
   {
     if (sx) { avma=av; return NULL; }
-    avma = (long)rem;
+    avma = (gpmem_t)rem;
     return gerepileupto(av,z-2);
   }
   lrem=i+3; rem -= lrem;
-  if (avma==av1) { avma = (long)rem; p1 = gcopy(p1); }
-  else p1 = gerepileupto((long)rem,p1);
+  if (avma==av1) { avma = (gpmem_t)rem; p1 = gcopy(p1); }
+  else p1 = gerepileupto((gpmem_t)rem,p1);
   rem[0]=evaltyp(t_POL) | evallg(lrem);
   rem[1]=evalsigne(1) | evalvarn(vx) | evallgef(lrem);
   rem += 2;
@@ -245,13 +251,13 @@ poldivres(GEN x, GEN y, GEN *pr)
   {
     av1=avma; p1 = (GEN)x[i];
     for (j=0; j<=i && j<=dz; j++)
-      if (y[i-j]) p1 = gadd(p1, gmul((GEN)z[j],(GEN)y[i-j]));
+      if (y[i-j] && z[j] != zero) p1 = gadd(p1, gmul((GEN)z[j],(GEN)y[i-j]));
     if (mod && avma==av1) p1 = gmul(p1,mod);
     rem[i]=avma==av1? lcopy(p1):lpileupto(av1,p1);
   }
   rem -= 2;
-  if (!sx) normalizepol_i(rem, lrem);
-  if (remainder) return gerepileupto(av,rem);
+  if (!sx) (void)normalizepol_i(rem, lrem);
+  if (pr == ONLY_REM) return gerepileupto(av,rem);
   z -= 2;
   {
     GEN *gptr[2]; gptr[0]=&z; gptr[1]=&rem;
@@ -286,7 +292,7 @@ factmod_init(GEN *F, GEN pp, long *p)
   f = gmul(f, mod(gun,pp));
   if (!signe(f)) err(zeropoler,"factmod");
   f = lift_intern(f); d = lgef(f);
-  for (i=2; i <d; i++) 
+  for (i=2; i <d; i++)
     if (typ(f[i])!=t_INT) err(impl,"factormod for general polynomials");
   *F = f; return d-3;
 }
@@ -352,7 +358,8 @@ GEN
 rootmod2(GEN f, GEN pp)
 {
   GEN g,y,ss,q,r, x_minus_s;
-  long p,av = avma,av1,d,i,nbrac;
+  long p, d, i, nbrac;
+  gpmem_t av = avma, av1;
 
   if (!(d = factmod_init(&f, pp, &p))) { avma=av; return cgetg(1,t_COL); }
   if (!p) err(talker,"prime too big in rootmod2");
@@ -389,15 +396,47 @@ rootmod2(GEN f, GEN pp)
   free(y); return g;
 }
 
+/* assume x reduced mod p, monic, squarefree. Return one root, or NULL if
+ * irreducible */
+static GEN
+quadsolvemod(GEN x, GEN p, int unknown)
+{
+  GEN b = (GEN)x[3], c = (GEN)x[2];
+  GEN s,u,D;
+
+  if (egalii(p, gdeux))
+  {
+    if (signe(c)) return NULL;
+    return gun;
+  }
+  D = subii(sqri(b), shifti(c,2));
+  D = resii(D,p);
+  if (unknown && kronecker(D,p) == -1) return NULL;
+
+  s = mpsqrtmod(D,p);
+  if (!s) err(talker,"not a prime in quadsolvemod");
+  u = addis(shifti(p,-1), 1); /* = 1/2 */
+  return modii(mulii(u, subii(s,b)), p);
+}
+
+static GEN
+otherroot(GEN x, GEN r, GEN p)
+{
+  GEN s = addii((GEN)x[3], r);
+  s = subii(p, s); if (signe(s) < 0) s = addii(s,p);
+  return s;
+}
+
 /* by splitting */
 GEN
 rootmod(GEN f, GEN p)
 {
-  long av = avma,tetpil,n,i,j,la,lb;
+  long n, i, j, la, lb;
+  gpmem_t av = avma, tetpil;
   GEN y,pol,a,b,q,pol0;
 
   if (!factmod_init(&f, p, &i)) { avma=av; return cgetg(1,t_COL); }
-  i = p[lgefint(p)-1];
+  i = modBIL(p);
   if ((i & 1) == 0) { avma = av; return root_mod_even(f,i); }
   i=2; while (!signe(f[i])) i++;
   if (i == 2) j = 1;
@@ -437,27 +476,28 @@ rootmod(GEN f, GEN p)
   if (la) y[j+lb] = (long)FpX_normalize(a,p);
   pol = gadd(polx[varn(f)], gun); pol0 = constant_term(pol);
   while (j<=n)
-  {
-    a=(GEN)y[j]; la=degpol(a);
+  { /* cf FpX_split_berlekamp */
+    a = (GEN)y[j]; la = degpol(a);
     if (la==1)
       y[j++] = lsubii(p, constant_term(a));
     else if (la==2)
     {
-      GEN d = subii(sqri((GEN)a[3]), shifti((GEN)a[2],2));
-      GEN e = mpsqrtmod(d,p), u = addis(q, 1); /* u = 1/2 */
-      y[j++] = lmodii(mulii(u,subii(e,(GEN)a[3])), p);
-      y[j++] = lmodii(mulii(u,negi(addii(e,(GEN)a[3]))), p);
+      GEN r = quadsolvemod(a, p, 0);
+      y[j++] = (long)r;
+      y[j++] = (long)otherroot(a,r, p);
     }
     else for (pol0[2]=1; ; pol0[2]++)
     {
       b = ZX_s_add(FpXQ_pow(pol,q, a,p), -1); /* pol^(p-1)/2 - 1 */
       b = FpX_gcd(a,b, p); lb = degpol(b);
-      if (lb && lb<la)
+      if (lb && lb < la)
       {
         b = FpX_normalize(b, p);
         y[j+lb] = (long)FpX_div(a,b, p);
         y[j]    = (long)b; break;
       }
+      if (pol0[2] == 100 && !BSW_psp(p))
+        err(talker, "not a prime in polrootsmod");
     }
   }
   tetpil = avma; y = gerepile(av,tetpil,sort(y));
@@ -483,28 +523,36 @@ rootmod0(GEN f, GEN p, long flag)
 /*                     FACTORISATION MODULO p                      */
 /*                                                                 */
 /*******************************************************************/
+#define FqX_div(x,y,T,p) FpXQX_divres((x),(y),(T),(p),NULL)
+#define FqX_rem(x,y,T,p) FpXQX_divres((x),(y),(T),(p),ONLY_REM)
+#define FqX_red FpXQX_red
+#define FqX_sqr FpXQX_sqr
 static GEN spec_FpXQ_pow(GEN x, GEN p, GEN S);
+extern GEN pol_to_vec(GEN x, long N);
+extern GEN FqXQ_pow(GEN x, GEN n, GEN S, GEN T, GEN p);
+extern GEN FpXQX_from_Kronecker(GEN z, GEN pol, GEN p);
+extern GEN FpXQX_safegcd(GEN P, GEN Q, GEN T, GEN p);
+extern GEN u_normalizepol(GEN x, long lx);
+extern GEN Fq_pow(GEN x, GEN n, GEN pol, GEN p);
 /* Functions giving information on the factorisation. */
 
 /* u in Z[X], return kernel of (Frob - Id) over Fp[X] / u */
 static GEN
-Berlekamp_ker(GEN u, GEN p)
+FpM_Berlekamp_ker(GEN u, GEN p)
 {
-  long i,j,d,N = degpol(u);
-  GEN vker,v,w,Q,p1,p2;
-  if (DEBUGLEVEL > 7) timer2();
+  long j,N = degpol(u);
+  GEN vker,v,w,Q,p1;
+  if (DEBUGLEVEL > 7) (void)timer2();
   Q = cgetg(N+1,t_MAT); Q[1] = (long)zerocol(N);
   w = v = FpXQ_pow(polx[varn(u)],p,u,p);
   for (j=2; j<=N; j++)
   {
-    Q[j] = lgetg(N+1,t_COL); p1 = (GEN)Q[j];
-    d = lgef(w)-1; p2 = w+1;
-    for (i=1; i<d ; i++) p1[i] = p2[i];
-    for (   ; i<=N; i++) p1[i] = zero;
+    p1 = pol_to_vec(w, N);
     p1[j] = laddis((GEN)p1[j], -1);
+    Q[j] = (long)p1;
     if (j < N)
     {
-      ulong av = avma;
+      gpmem_t av = avma;
       w = gerepileupto(av, FpX_res(gmul(w,v), u, p));
     }
   }
@@ -514,11 +562,36 @@ Berlekamp_ker(GEN u, GEN p)
   return vker;
 }
 
+static GEN
+FqM_Berlekamp_ker(GEN u, GEN T, GEN q, GEN p)
+{
+  long j,N = degpol(u);
+  GEN vker,v,w,Q,p1;
+  if (DEBUGLEVEL > 7) (void)timer2();
+  Q = cgetg(N+1,t_MAT); Q[1] = (long)zerocol(N);
+  w = v = FqXQ_pow(polx[varn(u)], q, u, T, p);
+  for (j=2; j<=N; j++)
+  {
+    p1 = pol_to_vec(w, N);
+    p1[j] = laddgs((GEN)p1[j], -1);
+    Q[j] = (long)p1;
+    if (j < N)
+    {
+      gpmem_t av = avma;
+      w = gerepileupto(av, FpXQX_divres(FpXQX_mul(w,v, T,p), u,T,p,ONLY_REM));
+    }
+  }
+  if (DEBUGLEVEL > 7) msgtimer("frobenius");
+  vker = FqM_ker(Q,T,p);
+  if (DEBUGLEVEL > 7) msgtimer("kernel");
+  return vker;
+}
+
 /* f in ZZ[X] and p a prime number. */
 long
 FpX_is_squarefree(GEN f, GEN p)
 {
-  long av = avma;
+  gpmem_t av = avma;
   GEN z;
   z = FpX_gcd(f,derivpol(f),p);
   avma = av;
@@ -533,7 +606,8 @@ FpX_is_squarefree(GEN f, GEN p)
 long
 FpX_nbroots(GEN f, GEN p)
 {
-  long av = avma, n=lgef(f);
+  long n=lgef(f);
+  gpmem_t av = avma;
   GEN z;
   if (n <= 4) return n-3;
   f = FpX_red(f, p);
@@ -545,10 +619,11 @@ FpX_nbroots(GEN f, GEN p)
 long
 FpX_is_totally_split(GEN f, GEN p)
 {
-  long av = avma, n=lgef(f);
+  long n=degpol(f);
+  gpmem_t av = avma;
   GEN z;
-  if (n <= 4) return 1;
-  if (!is_bigint(p) && n-3 > p[2]) return 0;
+  if (n <= 1) return 1;
+  if (!is_bigint(p) && n > p[2]) return 0;
   f = FpX_red(f, p);
   z = FpXQ_pow(polx[varn(f)], p, f, p);
   avma = av; return lgef(z)==4 && gcmp1((GEN)z[3]) && !signe(z[2]);
@@ -559,8 +634,8 @@ FpX_is_totally_split(GEN f, GEN p)
 long
 FpX_nbfact(GEN u, GEN p)
 {
-  ulong av = avma;
-  GEN vker = Berlekamp_ker(u,p);
+  gpmem_t av = avma;
+  GEN vker = FpM_Berlekamp_ker(u,p);
   avma = av; return lg(vker)-1;
 }
 
@@ -576,11 +651,11 @@ static GEN gsmul(GEN a,GEN b){return FpX_mul(a,b,modul
 GEN
 FpV_roots_to_pol(GEN V, GEN p, long v)
 {
-  ulong ltop=avma;
+  gpmem_t ltop=avma;
   long i;
   GEN g=cgetg(lg(V),t_VEC);
   for(i=1;i<lg(V);i++)
-    g[i]=(long)deg1pol(gun,negi((GEN)V[i]),v);
+    g[i]=(long)deg1pol(gun,modii(negi((GEN)V[i]),p),v);
   modulo=p;
   g=divide_conquer_prod(g,&gsmul);
   return gerepileupto(ltop,g);
@@ -595,24 +670,6 @@ trivfact(void)
   y[2]=lgetg(1,t_COL); return y;
 }
 
-static void
-fqunclone(GEN x, GEN a, GEN p)
-{
-  long i,j,lx = lgef(x);
-  for (i=2; i<lx; i++)
-  {
-    GEN p1 = (GEN)x[i];
-    if (typ(p1) == t_POLMOD) { p1[1] = (long)a; p1 = (GEN)p1[2]; }
-    if (typ(p1) == t_INTMOD) p1[1] = (long)p;
-    else /* t_POL */
-      for (j = lgef(p1)-1; j > 1; j--)
-      {
-        GEN p2 = (GEN)p1[j];
-        if (typ(p2) == t_INTMOD) p2[1] = (long)p;
-      }
-  }
-}
-
 static GEN
 try_pow(GEN w0, GEN pol, GEN p, GEN q, long r)
 {
@@ -639,7 +696,8 @@ try_pow(GEN w0, GEN pol, GEN p, GEN q, long r)
 static void
 split(long m, GEN *t, long d, GEN p, GEN q, long r, GEN S)
 {
-  long ps,l,v,dv,av0,av;
+  long ps, l, v, dv;
+  gpmem_t av0, av;
   GEN w,w0;
 
   dv=degpol(*t); if (dv==d) return;
@@ -648,7 +706,7 @@ split(long m, GEN *t, long d, GEN p, GEN q, long r, GE
   {
     if (ps==2)
     {
-      w0=w=gpuigs(polx[v],m-1); m+=2;
+      w0 = w = FpXQ_pow(polx[v], stoi(m-1), *t, gdeux); m += 2;
       for (l=1; l<d; l++)
         w = gadd(w0, spec_FpXQ_pow(w, p, S));
     }
@@ -672,7 +730,8 @@ split(long m, GEN *t, long d, GEN p, GEN q, long r, GE
 static void
 splitgen(GEN m, GEN *t, long d, GEN  p, GEN q, long r)
 {
-  long l,v,dv,av;
+  long l, v, dv;
+  gpmem_t av;
   GEN w;
 
   dv=degpol(*t); if (dv==d) return;
@@ -705,21 +764,22 @@ init_pow_p_mod_pT(GEN p, GEN T)
   S[1] = (long)FpXQ_pow(polx[v], p, T, p);
   /* use as many squarings as possible */
   for (i=2; i < n; i+=2)
-  {     
+  {
     p1 = gsqr((GEN)S[i>>1]);
     S[i]   = (long)FpX_res(p1, T, p);
     if (i == n-1) break;
     p1 = gmul((GEN)S[i], (GEN)S[1]);
     S[i+1] = (long)FpX_res(p1, T, p);
-  }       
+  }
   return S;
-} 
+}
 
 /* compute x^p, S is as above */
 static GEN
 spec_FpXQ_pow(GEN x, GEN p, GEN S)
 {
-  long av = avma, lim = stack_lim(av,1), i,dx = degpol(x);
+  long i, dx = degpol(x);
+  gpmem_t av = avma, lim = stack_lim(av, 1);
   GEN x0 = x+2, z;
   z = (GEN)x0[0];
   if (dx < 0) err(talker, "zero polynomial in FpXQ_pow. %Z not prime", p);
@@ -746,7 +806,8 @@ spec_FpXQ_pow(GEN x, GEN p, GEN S)
 static GEN
 factcantor0(GEN f, GEN pp, long flag)
 {
-  long i,j,k,d,e,vf,p,nbfact,tetpil,av = avma;
+  long i, j, k, d, e, vf, p, nbfact;
+  gpmem_t tetpil, av = avma;
   GEN ex,y,f2,g,g1,u,v,pd,q;
   GEN *t;
 
@@ -861,68 +922,170 @@ simplefactmod(GEN f, GEN p)
   return factcantor0(f,p,1);
 }
 
-/* vector of polynomials (in v) whose coeffs are given by the columns of x */
 GEN
-mat_to_vecpol(GEN x, long v)
+col_to_ff(GEN x, long v)
 {
-  long i,j, lx = lg(x), lcol = lg(x[1]);
-  GEN y = cgetg(lx, t_VEC);
+  long i, k = lg(x);
+  GEN p;
 
-  for (j=1; j<lx; j++)
+  while (--k && gcmp0((GEN)x[k]));
+  if (k <= 1) return k? (GEN)x[1]: gzero;
+  i = k+2; p = cgetg(i,t_POL);
+  p[1] = evalsigne(1) | evallgef(i) | evalvarn(v);
+  x--; for (k=2; k<i; k++) p[k] = x[k];
+  return p;
+}
+
+GEN
+vec_to_pol(GEN x, long v)
+{
+  long i, k = lg(x);
+  GEN p;
+
+  while (--k && gcmp0((GEN)x[k]));
+  if (!k) return zeropol(v);
+  i = k+2; p = cgetg(i,t_POL);
+  p[1] = evalsigne(1) | evallgef(i) | evalvarn(v);
+  x--; for (k=2; k<i; k++) p[k] = x[k];
+  return p;
+}
+
+GEN
+u_vec_to_pol(GEN x)
+{
+  long i, k = lg(x);
+  GEN p;
+
+  while (--k && !x[k]);
+  if (!k) return u_zeropol();
+  i = k+2; p = cgetg(i,t_POL);
+  p[1] = evalsigne(1) | evallgef(i) | evalvarn(0);
+  x--; for (k=2; k<i; k++) p[k] = x[k];
+  return p;
+}
+
+#if 0
+static GEN
+u_FpM_FpV_mul(GEN x, GEN y, ulong p)
+{
+  long i,k,l,lx=lg(x), ly=lg(y);
+  GEN z;
+  if (lx != ly) err(operi,"* [mod p]",x,y);
+  if (lx==1) return cgetg(1,t_VECSMALL);
+  l = lg(x[1]);
+  z = cgetg(l,t_VECSMALL);
+  for (i=1; i<l; i++)
   {
-    GEN p1, col = (GEN)x[j];
-    long k = lcol;
+    ulong p1 = 0;
+    for (k=1; k<lx; k++)
+    {
+      p1 += coeff(x,i,k) * y[k];
+      if (p1 & HIGHBIT) p1 %= p; 
+    }
+    z[i] = p1 % p;
+  }
+  return z;
+}
+#endif
 
-    while (k-- && gcmp0((GEN)col[k]));
-    i=k+2; p1=cgetg(i,t_POL);
-    p1[1] = evalsigne(1) | evallgef(i) | evalvarn(v);
-    col--; for (k=2; k<i; k++) p1[k] = col[k];
-    y[j] = (long)p1;
+/* u_vec_to_pol(u_FpM_FpV_mul(x, u_pol_to_vec(y), p)) */
+GEN
+u_FpM_FpX_mul(GEN x, GEN y, ulong p)
+{
+  long i,k,l, ly=lgef(y)-1;
+  GEN z;
+  if (ly==1) return u_zeropol();
+  l = lg(x[1]);
+  y++;
+  z = vecsmall_const(l,0) + 1;
+  for (k=1; k<ly; k++)
+  {
+    GEN c;
+    if (!y[k]) continue;
+    c = (GEN)x[k];
+    if (y[k] == 1)
+      for (i=1; i<l; i++)
+      {
+        z[i] += c[i];
+        if (z[i] & HIGHBIT) z[i] %= p; 
+      }
+    else
+      for (i=1; i<l; i++)
+      {
+        z[i] += c[i] * y[k];
+        if (z[i] & HIGHBIT) z[i] %= p; 
+      }
   }
+  for (i=1; i<l; i++) z[i] %= p;
+  while (--l && !z[l]);
+  if (!l) return u_zeropol();
+  *z-- = evalsigne(1) | evallgef(l+2) | evalvarn(0);
+  return z;
+}
+
+/* return the (N-dimensional) vector of coeffs of p */
+GEN
+pol_to_vec(GEN x, long N)
+{
+  long i, l;
+  GEN z = cgetg(N+1,t_COL);
+  if (typ(x) != t_POL)
+  {
+    z[1] = (long)x;
+    for (i=2; i<=N; i++) z[i]=zero;
+    return z;
+  }
+  l = lgef(x)-1; x++;
+  for (i=1; i<l ; i++) z[i]=x[i];
+  for (   ; i<=N; i++) z[i]=zero;
+  return z;
+}
+
+/* return the (N-dimensional) vector of coeffs of p */
+GEN
+u_pol_to_vec(GEN x, long N)
+{
+  long i, l;
+  GEN z = cgetg(N+1,t_VECSMALL);
+  if (typ(x) != t_VECSMALL) err(typeer,"u_pol_to_vec");
+  l = lgef(x)-1; x++;
+  for (i=1; i<l ; i++) z[i]=x[i];
+  for (   ; i<=N; i++) z[i]=0;
+  return z;
+}
+
+/* vector of polynomials (in v) whose coeffs are given by the columns of x */
+GEN
+mat_to_vecpol(GEN x, long v)
+{
+  long j, lx = lg(x);
+  GEN y = cgetg(lx, t_VEC);
+  for (j=1; j<lx; j++) y[j] = (long)vec_to_pol((GEN)x[j], v);
   return y;
 }
 
-/* matrix whose entries are given by the coeffs of the polynomials in 
+/* matrix whose entries are given by the coeffs of the polynomials in
  * vector v (considered as degree n-1 polynomials) */
 GEN
 vecpol_to_mat(GEN v, long n)
 {
-  long i,j,d,N = lg(v);
-  GEN p1,w, y = cgetg(N, t_MAT);
+  long j, N = lg(v);
+  GEN y = cgetg(N, t_MAT);
   if (typ(v) != t_VEC) err(typeer,"vecpol_to_mat");
-  n++;
-  for (j=1; j<N; j++)
-  {
-    p1 = cgetg(n,t_COL); y[j] = (long)p1;
-    w = (GEN)v[j];
-    if (typ(w) != t_POL) { p1[1] = (long)w; i=2; }
-    else
-    {
-      d=lgef(w)-1; w++;
-      for (i=1; i<d; i++) p1[i] = w[i];
-    }
-    for ( ; i<n; i++) p1[i] = zero;
-  }
+  for (j=1; j<N; j++) y[j] = (long)pol_to_vec((GEN)v[j], n);
   return y;
 }
+
 /* polynomial (in v) of polynomials (in w) whose coeffs are given by the columns of x */
 GEN
 mat_to_polpol(GEN x, long v,long w)
 {
-  long i,j, lx = lg(x), lcol = lg(x[1]);
+  long j, lx = lg(x);
   GEN y = cgetg(lx+1, t_POL);
   y[1]=evalsigne(1) | evallgef(lx+1) | evalvarn(v);
   y++;
-  for (j=1; j<lx; j++)
-  {
-    GEN p1, col = (GEN)x[j];
-    long k;
-    i=lcol+1; p1=cgetg(i,t_POL);
-    p1[1] = evalsigne(1) | evallgef(i) | evalvarn(w);
-    col--; for (k=2; k<i; k++) p1[k] = col[k];
-    y[j] = (long)normalizepol_i(p1,i);
-  }
-  return normalizepol_i(--y,lx+1);
+  for (j=1; j<lx; j++) y[j] = (long)vec_to_pol((GEN)x[j], w);
+  return normalizepol_i(--y, lx+1);
 }
 
 /* matrix whose entries are given by the coeffs of the polynomial v in
@@ -930,113 +1093,156 @@ mat_to_polpol(GEN x, long v,long w)
 GEN
 polpol_to_mat(GEN v, long n)
 {
-  long i,j,d,N = lgef(v)-1;
-  GEN p1,w, y = cgetg(N, t_MAT);
+  long j, N = lgef(v)-1;
+  GEN y = cgetg(N, t_MAT);
   if (typ(v) != t_POL) err(typeer,"polpol_to_mat");
-  n++;v++;
-  for (j=1; j<N; j++)
+  v++;
+  for (j=1; j<N; j++) y[j] = (long)pol_to_vec((GEN)v[j], n);
+  return y;
+}
+
+/* P(X,Y) --> P(Y,X), n-1 is the degree in Y */
+GEN
+swap_polpol(GEN x, long n, long w)
+{
+  long j, lx = lgef(x), ly = n+3;
+  long v=varn(x);
+  GEN y = cgetg(ly, t_POL);
+  y[1]=evalsigne(1) | evallgef(ly) | evalvarn(v);
+  for (j=2; j<ly; j++)
   {
-    p1 = cgetg(n,t_COL); y[j] = (long)p1;
-    w = (GEN)v[j];
-    if (typ(w) != t_POL) { p1[1] = (long)w; i=2; }
-    else
-    {
-      d=lgef(w)-1; w++;
-      for (i=1; i<d; i++) p1[i] = w[i];
-    }
-    for ( ; i<n; i++) p1[i] = zero;
+    long k;
+    GEN p1=cgetg(lx,t_POL);
+    p1[1] = evalsigne(1) | evallgef(lx) | evalvarn(w);
+    for (k=2; k<lx; k++)
+      if( j<lgef(x[k]))
+        p1[k] = mael(x,k,j);
+      else
+        p1[k] = zero;
+    y[j] = (long)normalizepol_i(p1,lx);
   }
-  return y;
+  return normalizepol_i(y,ly);
 }
 
-
 /* set x <-- x + c*y mod p */
 static void
-split_berlekamp_addmul(GEN x, GEN y, long c, long p)
+u_FpX_addmul(GEN x, GEN y, long c, long p)
 {
   long i,lx,ly,l;
   if (!c) return;
-  lx = lgef(x); ly = lgef(y); l = min(lx,ly);
+  lx = lgef(x);
+  ly = lgef(y); l = min(lx,ly);
   if (p & ~MAXHALFULONG)
   {
-    for (i=2; i<l;  i++) x[i] = ((ulong)x[i]+ (ulong)mulssmod(c,y[i],p)) % p; 
-    for (   ; i<ly; i++) x[i] = mulssmod(c,y[i],p);
+    for (i=2; i<l;  i++) x[i] = ((ulong)x[i]+ (ulong)mulssmod(c,y[i],p)) % p;
+    if (l == lx)
+      for (   ; i<ly; i++) x[i] = mulssmod(c,y[i],p);
   }
   else
   {
-    for (i=2; i<l;  i++) x[i] = ((ulong)x[i] + (ulong)(c*y[i])) % p; 
-    for (   ; i<ly; i++) x[i] = (c*y[i]) % p; 
+    for (i=2; i<l;  i++) x[i] = ((ulong)x[i] + (ulong)(c*y[i])) % p;
+    if (l == lx)
+      for (   ; i<ly; i++) x[i] = (c*y[i]) % p;
   }
-  do i--; while (i>1 && !x[i]);
-  if (i==1) setsigne(x,0); else { setsigne(x,1); setlgef(x,i+1); }
+  (void)u_normalizepol(x,i);
 }
 
+static long
+small_rand(long p)
+{
+  return (p==2)? ((mymyrand() & 0x1000) == 0): mymyrand() % p;
+}
+
+GEN
+FpX_rand(long d1, long v, GEN p)
+{
+  long i, d = d1+2;
+  GEN y;
+  y = cgetg(d,t_POL); y[1] = evalsigne(1) | evalvarn(v);
+  for (i=2; i<d; i++) y[i] = (long)genrand(p);
+  (void)normalizepol_i(y,d); return y;
+}
+
+/* return a random polynomial in F_q[v], degree < d1 */
+GEN
+FqX_rand(long d1, long v, GEN T, GEN p)
+{
+  long i, d = d1+2, k = degpol(T), w = varn(T);
+  GEN y;
+  y = cgetg(d,t_POL); y[1] = evalsigne(1) | evalvarn(v);
+  for (i=2; i<d; i++) y[i] = (long)FpX_rand(k, w, p);
+  (void)normalizepol_i(y,d); return y;
+}
+
+#define set_irred(i) { if ((i)>ir) swap(t[i],t[ir]); ir++;}
+
 long
-split_berlekamp(GEN *t, GEN pp, GEN pps2)
+FpX_split_berlekamp(GEN *t, GEN p)
 {
-  GEN u = *t, p1, p2, vker,pol;
-  long av,N = degpol(u), d,i,kk,l1,l2,p, vu = varn(u);
-  ulong av0 = avma;
-  
-  vker = Berlekamp_ker(u,pp);
+  GEN u = *t, a,b,po2,vker,pol;
+  long N = degpol(u), d, i, ir, L, la, lb, ps, vu = varn(u);
+  gpmem_t av0 = avma;
+
+  vker = FpM_Berlekamp_ker(u,p);
   vker = mat_to_vecpol(vker,vu);
   d = lg(vker)-1;
-  p = is_bigint(pp)? 0: pp[2];
-  if (p)
+  ps = is_bigint(p)? 0: p[2];
+  if (ps)
   {
-    avma = av0; p1 = cgetg(d+1, t_VEC); /* hack: hidden gerepile */
-    for (i=1; i<=d; i++) p1[i] = (long)pol_to_small((GEN)vker[i]);
-    vker = p1;
+    avma = av0; a = cgetg(d+1, t_VEC); /* hack: hidden gerepile */
+    for (i=1; i<=d; i++) a[i] = (long)pol_to_small((GEN)vker[i]);
+    vker = a;
   }
+  po2 = shifti(p, -1); /* (p-1) / 2 */
   pol = cgetg(N+3,t_POL);
-
-  for (kk=1; kk<d; )
+  ir = 0;
+  /* t[i] irreducible for i <= ir, still to be treated for i < L */
+  for (L=1; L<d; )
   {
     GEN polt;
-    if (p)
+    if (ps)
     {
-      if (p==2)
-      {
-        pol[2] = ((mymyrand() & 0x1000) == 0);
-        pol[1] = evallgef(pol[2]? 3: 2);
-        for (i=2; i<=d; i++)
-          split_berlekamp_addmul(pol,(GEN)vker[i],(mymyrand()&0x1000)?0:1, p);
-      }
-      else
-      {
-        pol[2] = mymyrand()%p; /* vker[1] = 1 */
-        pol[1] = evallgef(pol[2]? 3: 2);
-        for (i=2; i<=d; i++)
-          split_berlekamp_addmul(pol,(GEN)vker[i],mymyrand()%p, p);
-      }
+      pol[2] = small_rand(ps); /* vker[1] = 1 */
+      pol[1] = evallgef(pol[2]? 3: 2);
+      for (i=2; i<=d; i++)
+        u_FpX_addmul(pol,(GEN)vker[i],small_rand(ps), ps);
       polt = small_to_pol(pol,vu);
     }
     else
     {
-      pol[2] = (long)genrand(pp);
+      pol[2] = (long)genrand(p);
       pol[1] = evallgef(signe(pol[2])? 3: 2) | evalvarn(vu);
       for (i=2; i<=d; i++)
-        pol = gadd(pol, gmul((GEN)vker[i], genrand(pp)));
-      polt = FpX_red(pol,pp);
+        pol = gadd(pol, gmul((GEN)vker[i], genrand(p)));
+      polt = FpX_red(pol,p);
     }
-    for (i=1; i<=kk && kk<d; i++)
+    for (i=ir; i<L && L<d; i++)
     {
-      p1=t[i-1]; l1=degpol(p1);
-      if (l1>1)
+      a = t[i]; la = degpol(a);
+      if (la == 1) { set_irred(i); }
+      else if (la == 2)
       {
-        av = avma; p2 = FpX_res(polt, p1, pp);
-        if (lgef(p2) <= 3) { avma=av; continue; }
-        p2 = FpXQ_pow(p2,pps2, p1,pp);
-        if (!signe(p2)) err(talker,"%Z not a prime in split_berlekamp",pp);
-        p2 = ZX_s_add(p2, -1);
-        p2 = FpX_gcd(p1,p2, pp); l2=degpol(p2);
-        if (l2>0 && l2<l1)
+        GEN r = quadsolvemod(a,p,1);
+        if (r)
         {
-          p2 = FpX_normalize(p2, pp);
-          t[i-1] = p2; kk++;
-          t[kk-1] = FpX_div(p1,p2,pp);
-          if (DEBUGLEVEL > 7) msgtimer("new factor");
+          t[i] = deg1pol_i(gun, subii(p,r), vu); r = otherroot(a,r,p);
+          t[L] = deg1pol_i(gun, subii(p,r), vu); L++;
         }
+        set_irred(i);
+      }
+      else
+      {
+        gpmem_t av = avma;
+        b = FpX_res(polt, a, p);
+        if (degpol(b) == 0) { avma=av; continue; }
+        b = ZX_s_add(FpXQ_pow(b,po2, a,p), -1);
+        b = FpX_gcd(a,b, p); lb = degpol(b);
+        if (lb && lb < la)
+        {
+          b = FpX_normalize(b, p);
+          t[L] = FpX_div(a,b,p);
+          t[i]= b; L++;
+        }
         else avma = av;
       }
     }
@@ -1044,10 +1250,81 @@ split_berlekamp(GEN *t, GEN pp, GEN pps2)
   return d;
 }
 
+static GEN
+FqX_deriv(GEN f, GEN T, GEN p)
+{
+  return FqX_red(derivpol(f), T, p);
+}
 GEN
+FqX_gcd(GEN P, GEN Q, GEN T, GEN p)
+{
+  GEN g = T? FpXQX_safegcd(P,Q,T,p): FpX_gcd(P,Q,p);
+  if (!g) err(talker,"factmod9: %Z is reducible mod p!", T);
+  return g;
+}
+long
+FqX_is_squarefree(GEN P, GEN T, GEN p)
+{
+  gpmem_t av = avma;
+  GEN z = FqX_gcd(P, derivpol(P), T, p);
+  avma = av;
+  return degpol(z)==0;
+
+}
+
+long
+FqX_split_berlekamp(GEN *t, GEN q, GEN T, GEN p)
+{
+  GEN u = *t, a,b,qo2,vker,pol;
+  long N = degpol(u), vu = varn(u), vT = varn(T), dT = degpol(T);
+  long d, i, ir, L, la, lb;
+
+  vker = FqM_Berlekamp_ker(u,T,q,p);
+  vker = mat_to_vecpol(vker,vu);
+  d = lg(vker)-1;
+  qo2 = shifti(q, -1); /* (q-1) / 2 */
+  pol = cgetg(N+3,t_POL);
+  ir = 0;
+  /* t[i] irreducible for i < ir, still to be treated for i < L */
+  for (L=1; L<d; )
+  {
+    GEN polt;
+    pol[2] = (long)FpX_rand(dT,vT,p);
+    pol[1] = evallgef(signe(pol[2])? 3: 2) | evalvarn(vu);
+    for (i=2; i<=d; i++)
+      pol = gadd(pol, gmul((GEN)vker[i], FpX_rand(dT,vT,p)));
+    polt = FqX_red(pol,T,p);
+    for (i=ir; i<L && L<d; i++)
+    {
+      a = t[i]; la = degpol(a);
+      if (la == 1) { set_irred(i); }
+      else
+      {
+        gpmem_t av = avma;
+        b = FqX_rem(polt, a, T,p);
+        if (!degpol(b)) { avma=av; continue; }
+        b = FqXQ_pow(b,qo2, a,T,p);
+        if (!degpol(b)) { avma=av; continue; }
+        b[2] = ladd((GEN)b[2], gun);
+        b = FqX_gcd(a,b, T,p); lb = degpol(b);
+        if (lb && lb < la)
+        {
+          b = FpXQX_normalize(b, T,p);
+          t[L] = FqX_div(a,b,T,p);
+          t[i]= b; L++;
+        }
+        else avma = av;
+      }
+    }
+  }
+  return d;
+}
+
+GEN
 factmod0(GEN f, GEN pp)
 {
-  long i,j,k,e,p,N,nbfact,av = avma,tetpil,d;
+  long i, j, k, e, p, N, nbfact, d;
+  gpmem_t av = avma, tetpil;
   GEN pps2,ex,y,f2,p1,g1,u, *t;
 
   if (!(d = factmod_init(&f, pp, &p))) { avma=av; return trivfact(); }
@@ -1065,17 +1342,17 @@ factmod0(GEN f, GEN pp)
     {
       k++; if (p && !(k%p)) { k++; f2 = FpX_div(f2,g1,pp); }
       p1 = FpX_gcd(f2,g1, pp); u = g1; g1 = p1;
-      if (lgef(p1)!=3)
+      if (degpol(p1))
       {
         u = FpX_div( u,p1,pp);
         f2= FpX_div(f2,p1,pp);
       }
+      /* u is square-free (product of irred. of multiplicity e * k) */
       N = degpol(u);
       if (N)
       {
-        /* here u is square-free (product of irred. of multiplicity e * k) */
         t[nbfact] = FpX_normalize(u,pp);
-        d = (N==1)? 1: split_berlekamp(t+nbfact, pp, pps2);
+        d = (N==1)? 1: FpX_split_berlekamp(t+nbfact, pp);
         for (j=0; j<d; j++) ex[nbfact+j] = e*k;
         nbfact += d;
       }
@@ -1097,7 +1374,7 @@ factmod0(GEN f, GEN pp)
 GEN
 factmod(GEN f, GEN pp)
 {
-  long tetpil,av=avma;
+  gpmem_t tetpil, av=avma;
   long nbfact;
   long j;
   GEN y,u,v;
@@ -1135,7 +1412,14 @@ static GEN
 padic_pol_to_int(GEN f)
 {
   long i, l = lgef(f);
-  f = gdiv(f,content(f));
+  GEN c = content(f);
+  if (gcmp0(c)) /*  O(p^n) can occur */
+  {
+    if (typ(c) != t_PADIC) err(typeer,"padic_pol_to_int");
+    f = gdiv(f, gpowgs((GEN)c[2], valp(c)));
+  }
+  else
+    f = gdiv(f,c);
   for (i=2; i<l; i++)
     switch(typ(f[i]))
     {
@@ -1146,12 +1430,13 @@ padic_pol_to_int(GEN f)
   return f;
 }
 
-/* return invlead * (x + O(pr)), x in Z or Z_p, pr = p^r */
+/* return invlead * (x + O(pr)), x in Z or Z_p, pr = p^r. May return gzero */
 static GEN
 int_to_padic(GEN x, GEN p, GEN pr, long r, GEN invlead)
 {
   GEN p1,y;
-  long v,sx, av = avma;
+  long v, sx;
+  gpmem_t av = avma;
 
   if (typ(x) == t_PADIC)
   {
@@ -1171,16 +1456,25 @@ int_to_padic(GEN x, GEN p, GEN pr, long r, GEN invlead
   {
     y[4] = lmodii(p1,pr); r -= v;
   }
-  else 
+  else
   {
     y[4] = zero; v = r; r = 0;
   }
-  y[3] = (long)pr; 
+  y[3] = (long)pr;
   y[2] = (long)p;
   y[1] = evalprecp(r)|evalvalp(v);
   return invlead? gerepileupto(av, gmul(invlead,y)): y;
 }
 
+/* as above. always return a t_PADIC */
+static GEN
+strict_int_to_padic(GEN x, GEN p, GEN pr, long r, GEN invlead)
+{
+  GEN y = int_to_padic(x, p, pr, r, invlead);
+  if (y == gzero) y = ggrandocp(p,r);
+  return y;
+}
+
 /* return (x + O(p^r)) normalized (multiply by a unit such that leading coeff
  * is a power of p), x in Z[X] (or Z_p[X]) */
 static GEN
@@ -1192,7 +1486,7 @@ pol_to_padic(GEN x, GEN pr, GEN p, long r)
   if (gcmp1(lead)) lead = NULL;
   else
   {
-    long av = avma;
+    gpmem_t av = avma;
     v = ggval(lead,p);
     if (v) lead = gdiv(lead, gpowgs(p,v));
     lead = int_to_padic(lead,p,pr,r,NULL);
@@ -1206,139 +1500,169 @@ pol_to_padic(GEN x, GEN pr, GEN p, long r)
 static GEN
 padic_trivfact(GEN x, GEN p, long r)
 {
-  GEN p1, y = cgetg(3,t_MAT);
-  p1=cgetg(2,t_COL); y[1]=(long)p1; p = icopy(p);
-  p1[1]=(long)pol_to_padic(x,gpowgs(p,r),p,r);
-  p1=cgetg(2,t_COL); y[2]=(long)p1;
-  p1[1]=un; return y;
+  GEN y = cgetg(3,t_MAT);
+  p = icopy(p);
+  y[1]=(long)_col(pol_to_padic(x,gpowgs(p,r),p,r));
+  y[2]=(long)_col(gun); return y;
 }
 
-/* a etant un p-adique, retourne le vecteur des racines p-adiques de f
- * congrues a a modulo p dans le cas ou on suppose f(a) congru a 0 modulo p
- * (ou a 4 si p=2).
- */
+/* Assume x reduced mod p. q = p^e (simplified version of int_to_padic).
+ * If p = 2, is defined (and reduced) mod 4 [from rootmod] */
 GEN
-apprgen(GEN f, GEN a)
+Fp_to_Zp(GEN x, GEN p, GEN q, long e)
 {
-  GEN fp,p1,p,P,pro,x,x2,u,ip;
-  long av=avma,tetpil,v,Ps,i,j,k,lu,n,fl2;
+  GEN y = cgetg(5, t_PADIC);
+  if (egalii(p, x)) /* implies p = x = 2 */
+  {
+    x = gun; q = shifti(q, -1);
+    y[1] = evalprecp(e-1) | evalvalp(1);
+  }
+  else if (!signe(x)) y[1] = evalprecp(0) | evalvalp(e);
+  else                y[1] = evalprecp(e) | evalvalp(0);
+  y[2] = (long)p;
+  y[3] = (long)q;
+  y[4] = (long)x; return y;
+}
 
-  if (typ(f)!=t_POL) err(notpoler,"apprgen");
-  if (gcmp0(f)) err(zeropoler,"apprgen");
-  if (typ(a) != t_PADIC) err(rootper1);
+/* f != 0 in Z[X], primitive, a t_PADIC s.t f(a) = 0 mod p */
+static GEN
+apprgen_i(GEN f, GEN a)
+{
+  GEN fp,u,p,q,P,res,a0,rac;
+  long prec,i,j,k;
+
+  prec = gcmp0(a)? valp(a): precp(a);
+  if (prec <= 1) return _vec(a);
+  fp = derivpol(f); u = ggcd(f,fp);
+  if (degpol(u) > 0) { f = gdeuc(f,u); fp = derivpol(f); }
+  p = (GEN)a[2];
+  P = egalii(p,gdeux)? stoi(4): p;
+  a0= gmod(a, P);
+#if 0 /* assumption */
+  if (!gcmp0(FpX_eval(f,a0,P))) err(rootper2);
+#endif
+  if (!gcmp0(FpX_eval(fp,a0,p))) /* simple zero */
+  {
+    res = rootpadiclift(f, a0, p, prec);
+    res = strict_int_to_padic(res,p,gpowgs(p,prec),prec,NULL);
+    return _vec(res);
+  }
+
+  f = poleval(f, gadd(a, gmul(P,polx[varn(f)])));
   f = padic_pol_to_int(f);
-  fp=derivpol(f); p1=ggcd(f,fp);
-  if (lgef(p1)>3) { f=gdeuc(f,p1); fp=derivpol(f); }
-  p=(GEN)a[2]; p1=poleval(f,a);
-  v=ggval(p1,p); if (v <= 0) err(rootper2);
-  fl2=egalii(p,gdeux);
-  if (fl2 && v==1) err(rootper2);
-  v=ggval(poleval(fp,a),p);
-  if (!v) /* simple zero */
+  f = gdiv(f, gpowgs(p,ggval(f, p)));
+
+  res = cgetg(degpol(f)+1,t_VEC);
+  q = gpowgs(p,prec);
+  rac = lift_intern(rootmod(f, P));
+  for (j=i=1; i<lg(rac); i++)
   {
-    while (!gcmp0(p1))
-    {
-      a = gsub(a,gdiv(p1,poleval(fp,a)));
-      p1 = poleval(f,a);
-    }
-    tetpil=avma; pro=cgetg(2,t_VEC); pro[1]=lcopy(a);
-    return gerepile(av,tetpil,pro);
+    u = apprgen_i(f, Fp_to_Zp((GEN)rac[i], p,q,prec));
+    for (k=1; k<lg(u); k++) res[j++] = ladd(a, gmul(P,(GEN)u[k]));
   }
-  n=degpol(f); pro=cgetg(n+1,t_VEC);
+  setlg(res,j); return res;
+}
 
-  if (is_bigint(p)) err(impl,"apprgen for p>=2^31");
-  x = ggrandocp(p, valp(a) | precp(a)); 
-  if (fl2)
+/* a t_PADIC, return vector of p-adic roots of f equal to a (mod p)
+ * We assume 1) f(a) = 0 mod p (mod 4 if p=2).
+ *           2) leading coeff prime to p. */
+GEN
+apprgen(GEN f, GEN a)
+{
+  gpmem_t av = avma;
+  if (typ(f) != t_POL) err(notpoler,"apprgen");
+  if (gcmp0(f)) err(zeropoler,"apprgen");
+  if (typ(a) != t_PADIC) err(rootper1);
+  return gerepilecopy(av, apprgen_i(padic_pol_to_int(f), a));
+}
+
+static GEN
+rootpadic_i(GEN f, GEN p, long prec)
+{
+  GEN fp,y,z,q,rac;
+  long lx,i,j,k;
+
+  fp = derivpol(f); z = ggcd(f,fp);
+  if (degpol(z) > 0) { f = gdeuc(f,z); fp = derivpol(f); }
+  rac = rootmod(f, (egalii(p,gdeux) && prec>=2)? stoi(4): p);
+  rac = lift_intern(rac); lx = lg(rac);
+  if (prec==1)
   {
-    x2=ggrandocp(p,2); P = stoi(4);
+    y = cgetg(lx,t_COL);
+    for (i=1; i<lx; i++)
+      y[i] = (long)Fp_to_Zp((GEN)rac[i],p,p,1);
+    return y;
   }
-  else
+  y = cgetg(degpol(f)+1,t_COL);
+  q = gpowgs(p,prec);
+  for (j=i=1; i<lx; i++)
   {
-    x2=ggrandocp(p,1); P = p;
+    z = apprgen_i(f, Fp_to_Zp((GEN)rac[i], p,q,prec));
+    for (k=1; k<lg(z); k++,j++) y[j] = z[k];
   }
-  f = poleval(f, gadd(a,gmul(P,polx[varn(f)])));
-  if (!gcmp0(f)) f = gdiv(f,gpuigs(p,ggval(f, p)));
-  Ps = itos(P);
-  for (j=0,i=0; i<Ps; i++)
+  setlg(y,j); return y;
+}
+
+/* reverse x in place */
+static void
+polreverse(GEN x)
+{
+  long i, j;
+  if (typ(x) != t_POL) err(typeer,"polreverse");
+  for (i=2, j=lgef(x)-1; i<j; i++, j--) lswap(x[i], x[j]);
+  (void)normalizepol(x);
+}
+
+static GEN
+pnormalize(GEN f, GEN p, long prec, long n, GEN *plead, long *pprec, int *prev)
+{
+  *plead = leading_term(f);
+  *pprec = prec;
+  *prev = 0;
+  if (!gcmp1(*plead))
   {
-    ip=stoi(i);
-    if (gcmp0(poleval(f,gadd(ip,x2))))
+    long val = ggval(*plead,p), val1 = ggval(constant_term(f),p);
+    if (val1 < val)
     {
-      u=apprgen(f,gadd(x,ip)); lu=lg(u);
-      for (k=1; k<lu; k++)
-      {
-        j++; pro[j]=ladd(a,gmul(P,(GEN)u[k]));
-      }
+      *prev = 1; polreverse(f);
+     /* take care of loss of precision from leading coeff of factor
+      * (whose valuation is <= val) */
+      *pprec += val;
+      val = val1;
     }
+    *pprec += val * n;
   }
-  setlg(pro,j+1); return gerepilecopy(av,pro);
+  return pol_to_monic(f, plead);
 }
 
-/* Retourne le vecteur des racines p-adiques de f en precision r */
+/* return p-adic roots of f, precision prec */
 GEN
-rootpadic(GEN f, GEN p, long r)
+rootpadic(GEN f, GEN p, long prec)
 {
-  GEN fp,y,z,p1,pr,rac;
-  long lx,i,j,k,n,av=avma,tetpil,fl2;
+  gpmem_t av = avma;
+  GEN lead,y;
+  long PREC,i,k;
+  int reverse;
 
   if (typ(f)!=t_POL) err(notpoler,"rootpadic");
   if (gcmp0(f)) err(zeropoler,"rootpadic");
-  if (r<=0) err(rootper4);
+  if (prec <= 0) err(rootper4);
   f = padic_pol_to_int(f);
-  fp=derivpol(f); p1=ggcd(f,fp);
-  if (lgef(p1)>3) { f=gdeuc(f,p1); fp=derivpol(f); }
-  fl2=egalii(p,gdeux); rac=(fl2 && r>=2)? rootmod(f,stoi(4)): rootmod(f,p);
-  lx=lg(rac); p=gclone(p);
-  if (r==1)
-  {
-    tetpil=avma; y=cgetg(lx,t_COL);
-    for (i=1; i<lx; i++)
-    {
-      z=cgetg(5,t_PADIC); y[i]=(long)z;
-      z[1] = evalprecp(1)|evalvalp(0);
-      z[2] = z[3] = (long)p;
-      z[4] = lcopy(gmael(rac,i,2));
-    }
-    return gerepile(av,tetpil,y);
-  }
-  n=degpol(f); y=cgetg(n+1,t_COL);
-  j=0; pr = NULL;
-  z = cgetg(5,t_PADIC);
-  z[2] = (long)p;
-  for (i=1; i<lx; i++)
-  {
-    p1 = gmael(rac,i,2);
-    if (signe(p1))
-    {
-      if (!fl2 || mod2(p1))
-      {
-        z[1] = evalvalp(0)|evalprecp(r);
-        z[4] = (long)p1;
-      }
-      else
-      {
-        z[1] = evalvalp(1)|evalprecp(r);
-        z[4] = un;
-      }
-      if (!pr) pr=gpuigs(p,r);
-      z[3] = (long)pr;
-    }
-    else
-    {
-      z[1] = evalvalp(r);
-      z[3] = un;
-      z[4] = (long)p1;
-    }
-    p1 = apprgen(f,z);
-    for (k=1; k<lg(p1); k++) y[++j]=p1[k];
-  }
-  setlg(y,j+1); return gerepilecopy(av,y);
+  f = pnormalize(f, p, prec, 1, &lead, &PREC, &reverse);
+
+  y = rootpadic_i(f,p,PREC);
+  k = lg(y);
+  if (lead)
+    for (i=1; i<k; i++) y[i] = ldiv((GEN)y[i], lead);
+  if (reverse)
+    for (i=1; i<k; i++) y[i] = linv((GEN)y[i]);
+  return gerepilecopy(av, y);
 }
 /*************************************************************************/
 /*                             rootpadicfast                             */
 /*************************************************************************/
 
-/*lift accelerator. The author of the idea is unknown.*/
+/* lift accelerator */
 long hensel_lift_accel(long n, long *pmask)
 {
   long a,j;
@@ -1365,7 +1689,7 @@ long hensel_lift_accel(long n, long *pmask)
 
 /* STANDARD USE
    There exists p a prime number and a>0 such that
-   q=p^a    
+   q=p^a
    f in ZZ[X], with leading term prime to p.
    S must be a simple root mod p.
 
@@ -1375,7 +1699,7 @@ long hensel_lift_accel(long n, long *pmask)
 GEN
 rootpadiclift(GEN T, GEN S, GEN p, long e)
 {
-  ulong   ltop=avma;
+  gpmem_t ltop=avma;
   long    x;
   GEN     qold, q, qm1;
   GEN     W, Tr, Sr, Wr = gzero;
@@ -1387,7 +1711,7 @@ rootpadiclift(GEN T, GEN S, GEN p, long e)
   W=FpX_eval(deriv(Tr, x),S,q);
   W=mpinvmod(W,q);
   for(i=0;i<nb;i++)
-  {  
+  {
     qm1 = (mask&(1<<i))?sqri(qm1):mulii(qm1, q);
     q   =  mulii(qm1, p);
     Tr = FpX_red(T,q);
@@ -1424,7 +1748,7 @@ rootpadicliftroots(GEN f, GEN S, GEN q, long e)
     y[n-1]=(long) rootpadiclift(f, (GEN) S[n-1], q, e);
   else/* distinct-->totally split-->use trace trick */
   {
-    ulong av=avma;
+    gpmem_t av=avma;
     GEN z;
     z=(GEN)f[lgef(f)-2];/*-trace(roots)*/
     for(i=1; i<n-1;i++)
@@ -1441,12 +1765,12 @@ rootpadicliftroots(GEN f, GEN S, GEN q, long e)
  f must have no multiple roots mod p.
 
  return p-adics roots of f with prec pr, as integers (implicitly mod pr)
- 
+
 */
 GEN
 rootpadicfast(GEN f, GEN p, long e)
 {
-  ulong ltop=avma;
+  gpmem_t ltop=avma;
   GEN S,y;
   S=lift(rootmod(f,p));/*no multiplicity*/
   if (lg(S)==1)/*no roots*/
@@ -1461,13 +1785,13 @@ rootpadicfast(GEN f, GEN p, long e)
   return y;
 }
 /* Same as rootpadiclift for the polynomial X^n-a,
- * but here, n can be much larger. 
+ * but here, n can be much larger.
  * TODO: generalize to sparse polynomials.
  */
 GEN
 padicsqrtnlift(GEN a, GEN n, GEN S, GEN p, long e)
 {
-  ulong   ltop=avma;
+  gpmem_t ltop=avma;
   GEN     qold, q, qm1;
   GEN     W, Sr, Wr = gzero;
   long    i, nb, mask;
@@ -1476,7 +1800,7 @@ padicsqrtnlift(GEN a, GEN n, GEN S, GEN p, long e)
   W    = modii(mulii(n,powmodulo(S,subii(n,gun),q)),q);
   W    = mpinvmod(W,q);
   for(i=0;i<nb;i++)
-  {  
+  {
     qm1 = (mask&(1<<i))?sqri(qm1):mulii(qm1, q);
     q   =  mulii(qm1, p);
     Sr = S;
@@ -1512,29 +1836,29 @@ getprec(GEN x, long prec, GEN *p)
   return prec;
 }
 
-/* a appartenant a une extension finie de Q_p, retourne le vecteur des
- * racines de f congrues a a modulo p dans le cas ou on suppose f(a) congru a
- * 0 modulo p (ou a 4 si p=2).
- */
+/* a belongs to finite extension of Q_p, return all roots of f equal to a
+ * mod p. We assume f(a) = 0 (mod p) [mod 4 if p=2] */
 GEN
 apprgen9(GEN f, GEN a)
 {
-  GEN fp,p1,p,pro,x,x2,u,ip,t,vecg;
-  long av=avma,tetpil,v,ps_1,i,j,k,lu,n,prec,d,va,fl2;
+  GEN fp,p1,p,pro,x,x2,u,ip,T,vecg;
+  long v, ps_1, i, j, k, n, prec, d, va, fl2;
+  gpmem_t av=avma, tetpil;
 
   if (typ(f)!=t_POL) err(notpoler,"apprgen9");
   if (gcmp0(f)) err(zeropoler,"apprgen9");
   if (typ(a)==t_PADIC) return apprgen(f,a);
   if (typ(a)!=t_POLMOD || typ(a[2])!=t_POL) err(rootper1);
-  fp=derivpol(f); p1=ggcd(f,fp);
-  if (lgef(p1)>3) { f=gdeuc(f,p1); fp=derivpol(f); }
-  t=(GEN)a[1];
+
+  fp = derivpol(f); u = ggcd(f,fp);
+  if (degpol(u) > 0) { f = gdeuc(f,u); fp = derivpol(f); }
+  T = (GEN)a[1];
   prec = getprec((GEN)a[2], BIGINT, &p);
-  prec = getprec(t, prec, &p);
+  prec = getprec(T, prec, &p);
   if (prec==BIGINT) err(rootper1);
 
-  p1=poleval(f,a); v=ggval(lift_intern(p1),p); if (v<=0) err(rootper2);
-  fl2=egalii(p,gdeux);
+  p1 = poleval(f,a); v = ggval(lift_intern(p1),p); if (v<=0) err(rootper2);
+  fl2 = egalii(p,gdeux);
   if (fl2 && v==1) err(rootper2);
   v=ggval(lift_intern(poleval(fp,a)), p);
   if (!v)
@@ -1547,10 +1871,10 @@ apprgen9(GEN f, GEN a)
     tetpil=avma; pro=cgetg(2,t_COL); pro[1]=lcopy(a);
     return gerepile(av,tetpil,pro);
   }
-  n=degpol(f); pro=cgetg(n+1,t_COL); j=0;
+  n=degpol(f); pro=cgetg(n+1,t_COL);
 
   if (is_bigint(p)) err(impl,"apprgen9 for p>=2^31");
-  x=gmodulcp(ggrandocp(p,prec), t);
+  x=gmodulcp(ggrandocp(p,prec), T);
   if (fl2)
   {
     ps_1=3; x2=ggrandocp(p,2); p=stoi(4);
@@ -1560,21 +1884,19 @@ apprgen9(GEN f, GEN a)
     ps_1=itos(p)-1; x2=ggrandocp(p,1);
   }
   f = poleval(f,gadd(a,gmul(p,polx[varn(f)])));
-  if (!gcmp0(f)) f=gdiv(f,gpuigs(p,ggval(f,p)));
-  d=degpol(t); vecg=cgetg(d+1,t_COL);
+  f = gdiv(f,gpowgs(p,ggval(f,p)));
+
+  d=degpol(T); vecg=cgetg(d+1,t_COL);
   for (i=1; i<=d; i++)
     vecg[i] = (long)setloop(gzero);
-  va=varn(t);
+  va=varn(T); j = 1;
   for(;;) /* loop through F_q */
   {
-    ip=gmodulcp(gtopoly(vecg,va),t);
+    ip=gmodulcp(gtopoly(vecg,va),T);
     if (gcmp0(poleval(f,gadd(ip,x2))))
     {
-      u=apprgen9(f,gadd(ip,x)); lu=lg(u);
-      for (k=1; k<lu; k++)
-      {
-        j++; pro[j]=ladd(a,gmul(p,(GEN)u[k]));
-      }
+      u = apprgen9(f,gadd(ip,x));
+      for (k=1; k<lg(u); k++) pro[j++] = ladd(a, gmul(p,(GEN)u[k]));
     }
     for (i=d; i; i--)
     {
@@ -1584,12 +1906,17 @@ apprgen9(GEN f, GEN a)
     }
     if (!i) break;
   }
-  setlg(pro,j+1); return gerepilecopy(av,pro);
+  setlg(pro,j); return gerepilecopy(av,pro);
 }
 
-/*****************************************/
-/*  Factorisation p-adique d'un polynome */
-/*****************************************/
+/*******************************************************************/
+/*                                                                 */
+/*                      FACTORIZATION in Zp[X]                     */
+/*                                                                 */
+/*******************************************************************/
+extern GEN ZX_squff(GEN f, GEN *ex);
+extern GEN fact_from_DDF(GEN fa, GEN e, long n);
+
 int
 cmp_padic(GEN x, GEN y)
 {
@@ -1603,52 +1930,59 @@ cmp_padic(GEN x, GEN y)
   return cmpii((GEN)x[4], (GEN)y[4]);
 }
 
-/* factorise le polynome T=nf[1] dans Zp avec la precision pr */
+/*******************************/
+/*   Using Buchmann--Lenstra   */
+/*******************************/
+
+/* factor T = nf[1] in Zp to precision k */
 static GEN
-padicff2(GEN nf,GEN p,long pr)
+padicff2(GEN nf,GEN p,long k)
 {
-  long N=degpol(nf[1]),i,j,d,l;
-  GEN mat,V,D,fa,p1,pk,dec_p,pke,a,theta;
+  GEN z, mat, D, U,fa, pk, dec_p, Ui, mulx;
+  long i, l;
 
-  pk=gpuigs(p,pr); dec_p=primedec(nf,p);
-  l=lg(dec_p); fa=cgetg(l,t_COL);
+  mulx = eltmul_get_table(nf, gmael(nf,8,2)); /* mul by (x mod T) */
+  dec_p = primedec(nf,p);
+  l = lg(dec_p); fa = cgetg(l,t_COL);
+  D = NULL; /* -Wall */
   for (i=1; i<l; i++)
   {
-    p1 = (GEN)dec_p[i];
-    pke = idealpows(nf,p1, pr * itos((GEN)p1[3]));
-    p1=smith2(pke); V=(GEN)p1[3]; D=(GEN)p1[1];
-    for (d=1; d<=N; d++)
-      if (! egalii(gcoeff(V,d,d),pk)) break;
-    a=ginv(D); theta=gmael(nf,8,2); mat=cgetg(d,t_MAT);
-    for (j=1; j<d; j++)
-    {
-      p1 = gmul(D, element_mul(nf,theta,(GEN)a[j]));
-      setlg(p1,d); mat[j]=(long)p1;
-    }
-    fa[i]=(long)caradj(mat,0,NULL);
+    GEN P = (GEN)dec_p[i];
+    long e = itos((GEN)P[3]), ef = e * itos((GEN)P[4]);
+    D = smithall(idealpows(nf,P, k*e), &U, NULL);
+    Ui= ginv(U); setlg(Ui, ef+1); /* cf smithrel */
+    U = rowextract_i(U, 1, ef);
+    mat = gmul(U, gmul(mulx, Ui));
+    fa[i] = (long)caradj(mat,0,NULL);
   }
-  a = cgetg(l,t_COL); pk = icopy(pk);
+  pk = gcoeff(D,1,1); /* = p^k */
+  z = cgetg(l,t_COL); pk = icopy(pk);
   for (i=1; i<l; i++)
-    a[i] = (long)pol_to_padic((GEN)fa[i],pk,p,pr);
-  return a;
+    z[i] = (long)pol_to_padic((GEN)fa[i],pk,p,k);
+  return z;
 }
 
 static GEN
 padicff(GEN x,GEN p,long pr)
 {
-  GEN q,basden,bas,invbas,mul,dx,nf,mat;
-  long n=degpol(x),av=avma;
+  GEN q,basden,bas,invbas,mul,dx,dK,nf,fa,g,e;
+  long n=degpol(x);
+  gpmem_t av=avma;
 
-  nf=cgetg(10,t_VEC); nf[1]=(long)x; dx=discsr(x);
-  mat=cgetg(3,t_MAT); mat[1]=lgetg(3,t_COL); mat[2]=lgetg(3,t_COL);
-  coeff(mat,1,1)=(long)p; coeff(mat,1,2)=lstoi(pvaluation(dx,p,&q));
-  coeff(mat,2,1)=(long)q; coeff(mat,2,2)=un;
-  bas=allbase4(x,(long)mat,(GEN*)(nf+3),NULL);
-  if (!carrecomplet(divii(dx,(GEN)nf[3]),(GEN*)(nf+4)))
-    err(bugparier,"factorpadic2 (incorrect discriminant)");
+  nf=cgetg(10,t_VEC); nf[1]=(long)x;
+  fa = cgetg(3,t_MAT);
+  g = cgetg(3,t_COL); fa[1] = (long)g;
+  e = cgetg(3,t_COL); fa[2] = (long)e;
+  dx = ZX_disc(x);
+  g[1] = (long)p; e[1] = lstoi(pvaluation(absi(dx),p,&q));
+  g[2] = (long)q; e[2] = un;
+
+  bas = nfbasis(x, &dK, 0, fa);
+  nf[3] = (long)dK;
+  nf[4] = divise( diviiexact(dx, dK), p )? (long)p: un;
   basden = get_bas_den(bas);
   invbas = QM_inv(vecpol_to_mat(bas,n), gun);
-  mul = get_mul_table(x,basden,invbas,NULL);
+  mul = get_mul_table(x,basden,invbas);
   nf[7]=(long)bas;
   nf[8]=(long)invbas;
   nf[9]=(long)mul; nf[2]=nf[5]=nf[6]=zero;
@@ -1656,93 +1990,57 @@ padicff(GEN x,GEN p,long pr)
 }
 
 GEN
-factorpadic2(GEN x, GEN p, long r)
+factorpadic2(GEN f, GEN p, long prec)
 {
-  long av=avma,av2,k,i,j,i1,f,nbfac;
-  GEN res,p1,p2,y,d,a,ap,t,v,w;
-  GEN *fa;
+  gpmem_t av = avma;
+  GEN fa,ex,y;
+  long n,i,l;
 
-  if (typ(x)!=t_POL) err(notpoler,"factorpadic");
-  if (gcmp0(x)) err(zeropoler,"factorpadic");
-  if (r<=0) err(rootper4);
+  if (typ(f)!=t_POL) err(notpoler,"factorpadic");
+  if (gcmp0(f)) err(zeropoler,"factorpadic");
+  if (prec<=0) err(rootper4);
 
-  if (lgef(x)==3) return trivfact();
- if (!gcmp1(leading_term(x)))
+  n = degpol(f);
+  if (n==0) return trivfact();
+  if (n==1) return padic_trivfact(f,p,prec);
+  if (!gcmp1(leading_term(f)))
     err(impl,"factorpadic2 for non-monic polynomial");
-  if (lgef(x)==4) return padic_trivfact(x,p,r);
-  y=cgetg(3,t_MAT);
-  fa = (GEN*)new_chunk(lgef(x)-2);
-  d=content(x); a=gdiv(x,d);
-  ap=derivpol(a); t=ggcd(a,ap); v=gdeuc(a,t);
-  w=gdeuc(ap,t); j=0; f=1; nbfac=0;
-  while (f)
+  f = padic_pol_to_int(f);
+
+  fa = ZX_squff(f, &ex);
+  l = lg(fa); n = 0;
+  for (i=1; i<l; i++)
   {
-    j++; w=gsub(w,derivpol(v)); f=signe(w);
-    if (f) { res=ggcd(v,w); v=gdeuc(v,res); w=gdeuc(w,res); }
-    else res=v;
-    fa[j]=(lgef(res)>3) ? padicff(res,p,r) : cgetg(1,t_COL);
-    nbfac += (lg(fa[j])-1);
+    fa[i] = (long)padicff((GEN)fa[i],p,prec);
+    n += lg(fa[i])-1;
   }
-  av2=avma; y=cgetg(3,t_MAT);
-  p1=cgetg(nbfac+1,t_COL); y[1]=(long)p1;
-  p2=cgetg(nbfac+1,t_COL); y[2]=(long)p2;
-  for (i=1,k=0; i<=j; i++)
-    for (i1=1; i1<lg(fa[i]); i1++)
-    {
-      p1[++k]=lcopy((GEN)fa[i][i1]); p2[k]=lstoi(i);
-    }
-  y = gerepile(av,av2,y);
-  sort_factor(y, cmp_padic); return y;
+
+  y = fact_from_DDF(fa,ex,n);
+  return gerepileupto(av, sort_factor(y, cmp_padic));
 }
 
-/*******************************************************************/
-/*                                                                 */
-/*                FACTORISATION P-adique avec ROUND 4              */
-/*                                                                 */
-/*******************************************************************/
+/***********************/
+/*   Using ROUND 4     */
+/***********************/
 extern GEN Decomp(GEN p,GEN f,long mf,GEN theta,GEN chi,GEN nu,long r);
 extern GEN nilord(GEN p, GEN fx, long mf, GEN gx, long flag);
 extern GEN hensel_lift_fact(GEN pol, GEN Q, GEN T, GEN p, GEN pev, long e);
 
-static GEN
-squarefree(GEN f, GEN *ex)
+static int
+expo_is_squarefree(GEN e)
 {
-  GEN T,V,W,A,B;
-  long n,i,k;
-
-  T=ggcd(derivpol(f),f); V=gdeuc(f,T);
-  n=lgef(f)-2; A=cgetg(n,t_COL); B=cgetg(n,t_COL);
-  k=1; i=1;
-  do
-  {
-    W=ggcd(T,V); T=gdeuc(T,W);
-    if (lgef(V) != lgef(W))
-    {
-      A[i]=ldeuc(V,W); B[i]=k; i++;
-    }
-    k++; V=W;
-  }
-  while (lgef(V)>3);
-  setlg(A,i); *ex=B; return A;
+  long i, l = lg(e);
+  for (i=1; i<l; i++)
+    if (e[i] != 1) return 0;
+  return 1;
 }
 
-#define swap(x,y) { long _t=x; x=y; y=_t; }
-
-/* reverse x in place */
-static void
-polreverse(GEN x)
-{
-  long i, j;
-  if (typ(x) != t_POL) err(typeer,"polreverse");
-  for (i=2, j=lgef(x)-1; i<j; i++, j--) swap(x[i], x[j]);
-  (void)normalizepol(x);
-}
-
 GEN
 factorpadic4(GEN f,GEN p,long prec)
 {
-  GEN w,g,poly,fx,y,p1,p2,ex,pols,exps,ppow,lead;
-  long v=varn(f),n=degpol(f),av,tetpil,mfx,i,k,j,r,pr;
+  gpmem_t av = avma;
+  GEN w,g,poly,y,p1,p2,ex,pols,exps,ppow,lead;
+  long v=varn(f),n=degpol(f),mfx,i,k,j,r,pr;
   int reverse = 0;
 
   if (typ(f)!=t_POL) err(notpoler,"factorpadic");
@@ -1751,35 +2049,21 @@ factorpadic4(GEN f,GEN p,long prec)
   if (prec<=0) err(rootper4);
 
   if (n==0) return trivfact();
-  av=avma; f = padic_pol_to_int(f);
-  if (n==1) return gerepileupto(av, padic_trivfact(f,p,prec));
-  lead = leading_term(f); pr = prec;
-  if (!gcmp1(lead))
-  {
-    long val = ggval(lead,p), val1 = ggval(constant_term(f),p);
-    if (val1 < val)
-    {
-      reverse = 1; polreverse(f);
-     /* take care of loss of precision from leading coeff of factor
-      * (whose valuation is <= val) */
-      pr += val;
-      val = val1;
-    }
-    pr += val * (n-1);
-  }
-  f = pol_to_monic(f, &lead);
+  if (n==1) return padic_trivfact(f,p,prec);
+  f = padic_pol_to_int(f);
+  f = pnormalize(f, p, prec, n-1, &lead, &pr, &reverse);
 
-  poly=squarefree(f,&ex);
-  pols=cgetg(n+1,t_COL);
-  exps=cgetg(n+1,t_COL); n = lg(poly);
-  for (j=1,i=1; i<n; i++)
+  poly = ZX_squff(f,&ex);
+  pols = cgetg(n+1,t_COL);
+  exps = cgetg(n+1,t_COL); n = lg(poly);
+  for (j=i=1; i<n; i++)
   {
-    long av1 = avma;
-    fx=(GEN)poly[i]; mfx=ggval(discsr(fx),p);
-    w = (GEN)factmod(fx,p)[1];
-    if (!mfx)
+    gpmem_t av1 = avma;
+    GEN fx = (GEN)poly[i], fa = factmod0(fx,p);
+    w = (GEN)fa[1];
+    if (expo_is_squarefree((GEN)fa[2]))
     { /* no repeated factors: Hensel lift */
-      p1 = hensel_lift_fact(fx, lift_intern(w), NULL, p, gpowgs(p,pr), pr);
+      p1 = hensel_lift_fact(fx, w, NULL, p, gpowgs(p,pr), pr);
       p2 = stoi(ex[i]);
       for (k=1; k<lg(p1); k++,j++)
       {
@@ -1789,48 +2073,42 @@ factorpadic4(GEN f,GEN p,long prec)
       continue;
     }
     /* use Round 4 */
+    mfx = ggval(ZX_disc(fx),p);
     r = lg(w)-1;
-    g = lift_intern((GEN)w[r]);
+    g = (GEN)w[r];
     p2 = (r == 1)? nilord(p,fx,mfx,g,pr)
                  : Decomp(p,fx,mfx,polx[v],fx,g, (pr<=mfx)? mfx+1: pr);
     if (p2)
     {
-      p2 = gerepileupto(av1,p2);
+      p2 = gerepilecopy(av1,p2);
       p1 = (GEN)p2[1];
       p2 = (GEN)p2[2];
       for (k=1; k<lg(p1); k++,j++)
       {
-        pols[j]=p1[k];
-        exps[j]=lmulis((GEN)p2[k],ex[i]);
+        pols[j] = p1[k];
+        exps[j] = lmulis((GEN)p2[k],ex[i]);
       }
     }
     else
     {
-      avma=av1;
-      pols[j]=(long)fx;
-      exps[j]=lstoi(ex[i]); j++;
+      avma = av1;
+      pols[j] = (long)fx;
+      exps[j] = lstoi(ex[i]); j++;
     }
   }
   if (lead)
-  {
-    p1 = gmul(polx[v],lead);
     for (i=1; i<j; i++)
-    {
-      p2 = poleval((GEN)pols[i], p1);
-      pols[i] = ldiv(p2, content(p2));
-    }
-  }
-
-  tetpil=avma; y=cgetg(3,t_MAT);
-  p1 = cgetg(j,t_COL); ppow = gpowgs(p,prec); p = icopy(p);
+      pols[i] = (long)primpart( unscale_pol((GEN)pols[i], lead) );
+  y = cgetg(3,t_MAT);
+  p1 = cgetg(j,t_COL); p = icopy(p); ppow = gpowgs(p,prec);
   for (i=1; i<j; i++)
   {
     if (reverse) polreverse((GEN)pols[i]);
     p1[i] = (long)pol_to_padic((GEN)pols[i],ppow,p,prec);
   }
-  y[1]=(long)p1; setlg(exps,j);
-  y[2]=lcopy(exps); y = gerepile(av,tetpil,y);
-  sort_factor(y, cmp_padic); return y;
+  y[1] = (long)p1; setlg(exps,j);
+  y[2] = lcopy(exps);
+  return gerepileupto(av, sort_factor(y, cmp_padic));
 }
 
 GEN
@@ -1851,146 +2129,116 @@ factorpadic0(GEN f,GEN p,long r,long flag)
 /*                                                                 */
 /*******************************************************************/
 extern GEN to_Kronecker(GEN P, GEN Q);
-extern GEN from_Kronecker(GEN z, GEN pol);
-static GEN spec_Fq_pow_mod_pol(GEN x, GEN p, GEN a, GEN S);
+static GEN spec_Fq_pow_mod_pol(GEN x, GEN S, GEN T, GEN p);
 
 static GEN
-to_fq(GEN x, GEN a, GEN p)
+to_Fq(GEN x, GEN T, GEN p)
 {
-  long i,lx = lgef(x);
-  GEN z = cgetg(3,t_POLMOD), pol = cgetg(lx,t_POL);
-  pol[1] = x[1];
-  if (lx == 2) setsigne(pol, 0);
+  long i,lx, tx = typ(x);
+  GEN z = cgetg(3,t_POLMOD), y;
+
+  if (tx == t_INT)
+    y = mod(x,p);
   else
-    for (i=2; i<lx; i++) pol[i] = (long)mod((GEN)x[i], p);
-  /* assume deg(pol) < deg(a) */
-  z[1] = (long)a;
-  z[2] = (long)pol; return z;
+  {
+    if (tx != t_POL) err(typeer,"to_Fq");
+    lx = lgef(x);
+    y = cgetg(lx,t_POL);
+    y[1] = x[1];
+    if (lx == 2) setsigne(y, 0);
+    else
+      for (i=2; i<lx; i++) y[i] = (long)mod((GEN)x[i], p);
+  }
+  /* assume deg(y) < deg(T) */
+  z[1] = (long)T;
+  z[2] = (long)y; return z;
 }
 
-/* x POLMOD over Fq, return lift(x^n) */
 static GEN
-Kronecker_powmod(GEN x, GEN mod, GEN n)
+to_Fq_pol(GEN x, GEN T, GEN p)
 {
-  long lim,av,av0 = avma, i,j,m,v = varn(x);
-  GEN y, p1, p = NULL, pol = NULL;
-
-  for (i=lgef(mod)-1; i>1; i--)
-  {
-    p1 = (GEN)mod[i];
-    if (typ(p1) == t_POLMOD) { pol = (GEN)p1[1] ; break; }
-  }
-  if (!pol) err(talker,"need POLMOD coeffs in Kronecker_powmod");
-  for (i=lgef(pol)-1; i>1; i--)
-  {
-    p1 = (GEN)pol[i];
-    if (typ(p1) == t_INTMOD) { p = (GEN)p1[1] ; break; }
-  }
-  if (!p) err(talker,"need Fq coeffs in Kronecker_powmod");
-  x = lift_intern(to_Kronecker(x,pol));
-
-  /* adapted from powgi */
-  av=avma; lim=stack_lim(av,1);
-  p1 = n+2; m = *p1;
-
-  y=x; j=1+bfffo(m); m<<=j; j = BITS_IN_LONG-j;
-  for (i=lgefint(n)-2;;)
-  {
-    for (; j; m<<=1,j--)
-    {
-      y = gsqr(y);
-      y = from_Kronecker(FpX(y,p), pol);
-      setvarn(y, v);
-      y = gres(y, mod);
-      y = lift_intern(to_Kronecker(y,pol));
-
-      if (m<0)
-      {
-        y = gmul(y,x);
-        y = from_Kronecker(FpX(y,p), pol);
-        setvarn(y, v);
-        y = gres(y, mod);
-        y = lift_intern(to_Kronecker(y,pol));
-      }
-      if (low_stack(lim, stack_lim(av,1)))
-      {
-        if(DEBUGMEM>1) err(warnmem,"Kronecker_powmod");
-        y = gerepilecopy(av, y);
-      }
-    }
-    if (--i == 0) break;
-    m = *++p1, j = BITS_IN_LONG;
-  }
-  y = from_Kronecker(FpX(y,p),pol);
-  setvarn(y, v); return gerepileupto(av0, y);
+  long i, lx, tx = typ(x);
+  if (tx != t_POL) err(typeer,"to_Fq_pol");
+  lx = lgef(x);
+  for (i=2; i<lx; i++) x[i] = (long)to_Fq((GEN)x[i],T,p);
+  return x;
 }
 
-GEN
-FpX_rand(long d1, long v, GEN p)
+/* assume T a clone */
+static GEN
+copy_then_free(GEN T)
 {
-  long i, d = d1+2;
-  GEN y;
-  y = cgetg(d,t_POL); y[1] = evalsigne(1) | evalvarn(v);
-  for (i=2; i<d; i++) y[i] = (long)genrand(p);
-  normalizepol_i(y,d); return y;
+  GEN t = forcecopy(T); gunclone(T); return t;
 }
 
-/* return a random polynomial in F_q[v], degree < d1 */
 static GEN
-FqX_rand(long d1, long v, GEN p, GEN a)
+to_Fq_fact(GEN t, GEN ex, long nbf, int sort, GEN unfq, gpmem_t av)
 {
-  long i,j, d = d1+2, k = lgef(a)-1;
-  GEN y,t;
+  GEN T = (GEN)unfq[1]/*clone*/, y = cgetg(3,t_MAT), u,v,p;
+  long j,k, l = lg(t);
 
-  y = cgetg(d,t_POL); y[1] = evalsigne(1) | evalvarn(v);
-  t = cgetg(k,t_POL); t[1] = a[1];
-  for (i=2; i<d; i++)
-  {
-    for (j=2; j<k; j++) t[j] = (long)genrand(p);
-    normalizepol_i(t,k); y[i]=(long)to_fq(t,a,p);
-  }
-  normalizepol_i(y,d); return y;
+  u = cgetg(nbf,t_COL); y[1] = (long)u;
+  v = cgetg(nbf,t_COL); y[2] = (long)v;
+  for (j=1,k=0; j<l; j++)
+    if (ex[j])
+    {
+      k++;
+      u[k] = (long)simplify((GEN)t[j]); /* may contain pols of degree 0 */
+      v[k] = lstoi(ex[j]);
+    }
+  y = gerepileupto(av, y);
+  if (sort) y = sort_factor(y,cmp_pol);
+  T = copy_then_free(T);
+  p = (GEN)leading_term(T)[1];
+  u = (GEN)y[1];
+  for (j=1; j<nbf; j++) u[j] = (long)to_Fq_pol((GEN)u[j], T,p);
+  return y;
 }
 
 /* split into r factors of degree d */
 static void
-split9(GEN *t, long d, GEN p, GEN q, GEN a, GEN S)
+FqX_split(GEN *t, long d, GEN q, GEN S, GEN T, GEN p)
 {
-  long l,v,av,is2,cnt, dt = degpol(*t), da = degpol(a);
+  long l, v, is2, cnt, dt = degpol(*t), dT = degpol(T);
+  gpmem_t av;
   GEN w,w0;
 
   if (dt == d) return;
   v = varn(*t);
-  if (DEBUGLEVEL > 6) timer2();
+  if (DEBUGLEVEL > 6) (void)timer2();
   av = avma; is2 = egalii(p, gdeux);
   for(cnt = 1;;cnt++)
   { /* splits *t with probability ~ 1 - 2^(1-r) */
-    w = w0 = FqX_rand(dt,v, p,a);
+    w = w0 = FqX_rand(dt,v, T,p);
     for (l=1; l<d; l++) /* sum_{0<i<d} w^(q^i), result in (F_q)^r */
-      w = gadd(w0, spec_Fq_pow_mod_pol(w, p, a, S));
+      w = gadd(w0, spec_Fq_pow_mod_pol(w, S, T, p));
+    w = FqX_red(w, T,p);
     if (is2)
     {
       w0 = w;
-      for (l=1; l<da; l++) /* sum_{0<i<k} w^(2^i), result in (F_2)^r */
-        w = gadd(w0, gres(gsqr(w), *t));
+      for (l=1; l<dT; l++) /* sum_{0<i<k} w^(2^i), result in (F_2)^r */
+      {
+        w = FqX_rem(FqX_sqr(w,T,p), *t, T,p);
+        w = FqX_red(gadd(w0,w), NULL,p);
+      }
     }
     else
     {
-      w = Kronecker_powmod(w, *t, shifti(q,-1));
+      w = FqXQ_pow(w, shifti(q,-1), *t, T, p);
       /* w in {-1,0,1}^r */
-      if (lgef(w) == 3) continue;
+      if (!degpol(w)) continue;
       w[2] = ladd((GEN)w[2], gun);
     }
-    w = ggcd(*t,w); l = degpol(w);
+    w = FqX_gcd(*t,w, T,p); l = degpol(w);
     if (l && l != dt) break;
     avma = av;
   }
   w = gerepileupto(av,w);
-  if (DEBUGLEVEL > 6) 
-    fprintferr("[split9] time for splitting: %ld (%ld trials)\n",timer2(),cnt);
-  l /= d; t[l]=gdeuc(*t,w); *t=w;
-  split9(t+l,d,p,q,a,S);
-  split9(t  ,d,p,q,a,S);
+  if (DEBUGLEVEL > 6)
+    fprintferr("[FqX_split] splitting time: %ld (%ld trials)\n",timer2(),cnt);
+  l /= d; t[l] = FqX_div(*t,w, T,p); *t = w;
+  FqX_split(t+l,d,q,S,T,p);
+  FqX_split(t  ,d,q,S,T,p);
 }
 
 /* to "compare" (real) scalars and t_INTMODs */
@@ -2019,52 +2267,52 @@ cmp_pol(GEN x, GEN y)
   return 0;
 }
 
-/* assume n > 1, X a POLMOD over Fq */
-/* return S = [ X^q, X^2q, ... X^(n-1)q ] mod T (in Fq[X]) in Kronecker form */
+/* assume n > 1, X over FqX */
+/* return S = [ X^q, X^2q, ... X^(n-1)q ] mod u (in Fq[X]) in Kronecker form */
 static GEN
-init_pow_q_mod_pT(GEN Xmod, GEN q, GEN a, GEN T)
+init_pow_q_mod_pT(GEN X, GEN q, GEN u, GEN T, GEN p)
 {
-  long i, n = degpol(T);
+  long i, n = degpol(u);
   GEN p1, S = cgetg(n, t_VEC);
 
-  S[1] = (long)Kronecker_powmod((GEN)Xmod[2], (GEN)Xmod[1], q);
+  S[1] = (long)FqXQ_pow(X, q, u, T, p);
 #if 1 /* use as many squarings as possible */
   for (i=2; i < n; i+=2)
-  {     
+  {
     p1 = gsqr((GEN)S[i>>1]);
-    S[i]   = lres(p1, T);
+    S[i]   = (long)FqX_rem(p1, u, T,p);
     if (i == n-1) break;
     p1 = gmul((GEN)S[i], (GEN)S[1]);
-    S[i+1] = lres(p1, T);
-  }       
+    S[i+1] = (long)FqX_rem(p1, u, T,p);
+  }
 #else
   for (i=2; i < n; i++)
-  { 
+  {
     p1 = gmul((GEN)S[i-1], (GEN)S[1]);
-    S[i] = lres(p1, T);
-  } 
+    S[i] = (long)FqX_rem(p1, u, T,p);
+  }
 #endif
   for (i=1; i < n; i++)
-    S[i] = (long)lift_intern(to_Kronecker((GEN)S[i], a));
+    S[i] = (long)to_Kronecker((GEN)S[i], T);
   return S;
 }
 
 /* compute x^q, S is as above */
 static GEN
-spec_Fq_pow_mod_pol(GEN x, GEN p, GEN a, GEN S)
+spec_Fq_pow_mod_pol(GEN x, GEN S, GEN T, GEN p)
 {
-  long av = avma, lim = stack_lim(av,1), i,dx = degpol(x);
+  long i, dx = degpol(x);
+  gpmem_t av = avma, lim = stack_lim(av, 1);
   GEN x0 = x+2, z,c;
 
-  c = (GEN)x0[0];
-  z = lift_intern(lift(c));
+  z = c = (GEN)x0[0];
   for (i = 1; i <= dx; i++)
   {
     GEN d;
     c = (GEN)x0[i];
     if (gcmp0(c)) continue;
     d = (GEN)S[i];
-    if (!gcmp1(c)) d = gmul(lift_intern(lift(c)),d);
+    if (!gcmp1(c)) d = gmul(c,d);
     z = gadd(z, d);
     if (low_stack(lim, stack_lim(av,1)))
     {
@@ -2072,125 +2320,108 @@ spec_Fq_pow_mod_pol(GEN x, GEN p, GEN a, GEN S)
       z = gerepileupto(av, z);
     }
   }
-  z = FpX(z, p);
-  z = from_Kronecker(z, a);
+  z = FpXQX_from_Kronecker(z, T, p);
   setvarn(z, varn(x)); return gerepileupto(av, z);
 }
 
-static long isabsolutepol(GEN f, GEN pp, GEN a)
+static long
+isabsolutepol(GEN f)
 {
-  int i,res=1;
-  GEN c;
-  for(i=2; i<lg(f); i++)
+  int i, l = lgef(f);
+  for(i=2; i<l; i++)
   {
-    c = (GEN) f[i];
-    switch(typ(c))
-    {
-    case t_INT: /* OK*/
-      break;
-    case t_INTMOD:
-      if (gcmp((GEN)c[1],pp))
-	err(typeer,"factmod9");
-      break;
-    case t_POLMOD:
-      if (gcmp((GEN)c[1],a))
-	err(typeer,"factmod9");
-      isabsolutepol((GEN)c[1],pp,gzero);
-      isabsolutepol((GEN)c[2],pp,gzero);
-      if (degpol(c[1])>0)
-	res = 0;
-      break;
-    case t_POL:
-      isabsolutepol(c,pp,gzero);
-      if (degpol(c)>0)
-	res = 0;
-      break;
-    default:
-      err(typeer,"factmod9");
-    }
+    GEN c = (GEN)f[i];
+    if (typ(c) == t_POL && degpol(c) > 0) return 0;
   }
-  return res;
+  return 1;
 }
 
 GEN
-factmod9(GEN f, GEN pp, GEN a)
+factmod9(GEN f, GEN p, GEN T)
 {
-  long av = avma, tetpil,p,i,j,k,d,e,vf,va,nbfact,nbf,pk;
-  GEN S,ex,y,f2,f3,df1,df2,g,g1,xmod,u,v,qqd,qq,unfp,unfq, *t;
+  gpmem_t av = avma;
+  long pg, i, j, k, d, e, N, vf, va, nbfact, nbf, pk;
+  GEN S,ex,f2,f3,df1,df2,g1,u,v,q,unfp,unfq, *t;
   GEN frobinv,X;
 
-  if (typ(a)!=t_POL || typ(f)!=t_POL || gcmp0(a)) err(typeer,"factmod9");
-  vf=varn(f); va=varn(a);
-  if (va<=vf) err(talker,"polynomial variable must be of higher priority than finite field\nvariable in factorff");
-  if (isabsolutepol(f, pp, a))
+  if (typ(T)!=t_POL || typ(f)!=t_POL || gcmp0(T)) err(typeer,"factmod9");
+  vf = varn(f);
+  va = varn(T);
+  if (va <= vf)
+    err(talker,"polynomial variable must have higher priority in factorff");
+  unfp = gmodulsg(1,p); T = gmul(unfp,T);
+  unfq = gmodulo(gmul(unfp,polun[va]), T);
+  f = gmul(unfq,f);
+  if (!signe(f)) err(zeropoler,"factmod9");
+  d = degpol(f); if (!d) { avma = av; return trivfact(); }
+
+  f = simplify(lift_intern(lift_intern(f)));
+  T = lift_intern(T);
+  if (isabsolutepol(f))
   {
-    GEN z= Fp_factor_rel0(simplify(lift(lift(f))), pp, lift(a));
-    GEN t=(GEN)z[1],ex=(GEN)z[2];
-    unfp = gmodulsg(1,pp);
-    unfq = gmodulcp(gmul(unfp, polun[va]), gmul(unfp,a));
-    nbfact=lg(t);
-    tetpil=avma;
-    y=cgetg(3,t_MAT);
-    u=cgetg(nbfact,t_COL); y[1]=(long)u;
-    v=cgetg(nbfact,t_COL); y[2]=(long)v;
-    for (j=1; j<nbfact; j++)
-    {
-      u[j] = lmul((GEN)t[j],unfq);
-      v[j] = lstoi(ex[j]);
-    }
-    return gerepile(av,tetpil,y);
+    GEN z = Fp_factor_rel0(f, p, T);
+    return to_Fq_fact((GEN)z[1],(GEN)z[2],lg(z[1]), 0, unfq,av);
   }
-  p = is_bigint(pp)? 0: itos(pp);
-  unfp=gmodulsg(1,pp); a=gmul(unfp,a);
-  unfq=gmodulo(gmul(unfp,polun[va]), a); a = (GEN)unfq[1];
-  f = gmul(unfq,f); if (!signe(f)) err(zeropoler,"factmod9");
-  d = degpol(f); if (!d) { avma=av; gunclone(a); return trivfact(); }
 
+  pg = is_bigint(p)? 0: itos(p);
   S = df2  = NULL; /* gcc -Wall */
-  pp = gmael(a,2,1); /* out of the stack */
   t = (GEN*)cgetg(d+1,t_VEC); ex = new_chunk(d+1);
 
-  frobinv = gpowgs(pp, lgef(a)-4);
-  xmod = cgetg(3,t_POLMOD);
-  X = gmul(polx[vf],unfq);
-  xmod[2] = (long)X;
-  qq=gpuigs(pp,degpol(a));
+  frobinv = gpowgs(p, degpol(T)-1);
+
+  X = polx[vf];
+  q = gpowgs(p, degpol(T));
   e = nbfact = 1;
-  pk=1; df1=derivpol(f); f3=NULL;
+  pk = 1;
+  f3 = NULL;
+  df1 = FqX_deriv(f, T, p);
   for(;;)
   {
-    long du,dg;
     while (gcmp0(df1))
-    { /* needs d >= pp: p = 0 can't happen  */
-      pk *= p; e=pk;
-      j=(degpol(f))/p+3; setlg(f,j); setlgef(f,j);
-      for (i=2; i<j; i++) f[i] = (long)powgi((GEN)f[p*(i-2)+2], frobinv);
-      df1=derivpol(f); f3=NULL;
+    { /* needs d >= p: pg = 0 can't happen  */
+      pk *= pg; e = pk;
+      j = (degpol(f))/pg+3; setlg(f,j); setlgef(f,j);
+      for (i=2; i<j; i++)
+        f[i] = (long)Fq_pow((GEN)f[pg*(i-2)+2], frobinv, T,p);
+      df1 = FqX_deriv(f, T, p); f3 = NULL;
     }
-    f2 = f3? f3: ggcd(f,df1);
-    if (lgef(f2)==3) u = f;
+    f2 = f3? f3: FqX_gcd(f,df1, T,p);
+    if (!degpol(f2)) u = f;
     else
     {
-      g1=gdeuc(f,f2); df2=derivpol(f2);
-      if (gcmp0(df2)) { u=g1; f3=f2; }
+      g1 = FqX_div(f,f2, T,p); df2 = FqX_deriv(f2, T,p);
+      if (gcmp0(df2)) { u = g1; f3 = f2; }
       else
       {
-        f3=ggcd(f2,df2);
-        if (lgef(f3)==3) u=gdeuc(g1,f2);
+        f3 = FqX_gcd(f2,df2, T,p);
+        if (!degpol(f3)) u = FqX_div(g1,f2, T,p);
         else
-          u=gdeuc(g1,gdeuc(f2,f3));
+          u = FqX_div(g1, FqX_div(f2,f3, T,p), T,p);
       }
     }
     /* u is square-free (product of irreducibles of multiplicity e) */
-    qqd=gun; xmod[1]=(long)u;
-    
-    du = degpol(u); v = X;
-    if (du > 1) S = init_pow_q_mod_pT(xmod, qq, a, u);
-    for (d=1; d <= du>>1; d++)
+
+#if 0
+    N = degpol(u);
+    if (N)
     {
-      qqd=mulii(qqd,qq);
-      v = spec_Fq_pow_mod_pol(v, pp, a, S);
-      g = ggcd(gsub(v,X),u);
+      t[nbfact] = FpXQX_normalize(u, T,p);
+      d = (N==1)? 1: FqX_split_berlekamp(t+nbfact, q, T, p);
+      for (j=0; j<d; j++) ex[nbfact+j] = e;
+      nbfact += d;
+    }
+#else
+{
+    GEN qqd = gun, g;
+    long dg;
+
+    N = degpol(u); v = X;
+    if (N > 1) S = init_pow_q_mod_pT(X, q, u, T, p);
+    for (d=1; d <= N>>1; d++)
+    {
+      qqd = mulii(qqd,q);
+      v = spec_Fq_pow_mod_pol(v, S, T, p);
+      g = FqX_gcd(gsub(v,X),u, T,p);
       dg = degpol(g);
       if (dg <= 0) continue;
 
@@ -2198,22 +2429,24 @@ factmod9(GEN f, GEN pp, GEN a)
       j = nbfact+dg/d;
 
       t[nbfact] = g;
-      split9(t+nbfact,d,pp,qq,a,S);
-      for (; nbfact<j; nbfact++) ex[nbfact]=e;
-      du -= dg;
-      u = gdeuc(u,g);
-      v = gres(v,u); xmod[1] = (long)u;
+      FqX_split(t+nbfact,d,q,S,T,p);
+      for (; nbfact<j; nbfact++) ex[nbfact] = e;
+      N -= dg;
+      u = FqX_div(u,g, T,p);
+      v = FqX_rem(v,u, T,p);
     }
-    if (du) { t[nbfact]=u; ex[nbfact++]=e; }
-    if (lgef(f2) == 3) break;
+    if (N) { t[nbfact] = u; ex[nbfact++] = e; }
+}
+#endif
+    if (!degpol(f2)) break;
 
-    f=f2; df1=df2; e += pk;
+    f = f2; df1 = df2; e += pk;
   }
 
-  nbf=nbfact; tetpil=avma; y=cgetg(3,t_MAT);
+  nbf = nbfact; setlg(t, nbfact);
   for (j=1; j<nbfact; j++)
   {
-    t[j]=gdiv((GEN)t[j],leading_term(t[j]));
+    t[j] = FpXQX_normalize((GEN)t[j], T,p);
     for (k=1; k<j; k++)
       if (ex[k] && gegal(t[j],t[k]))
       {
@@ -2221,25 +2454,7 @@ factmod9(GEN f, GEN pp, GEN a)
         nbf--; break;
       }
   }
-  u=cgetg(nbf,t_COL); y[1]=(long)u;
-  v=cgetg(nbf,t_COL); y[2]=(long)v;
-  for (j=1,k=0; j<nbfact; j++)
-    if (ex[j])
-    {
-      k++;
-      u[k]=(long)t[j];
-      v[k]=lstoi(ex[j]);
-    }
-  y = gerepile(av,tetpil,y);
-  u=(GEN)y[1];
-  { /* put a back on the stack */
-    GEN tokill = a;
-    a = forcecopy(a);
-    gunclone(tokill);
-  }
-  pp = (GEN)leading_term(a)[1];
-  for (j=1; j<nbf; j++) fqunclone((GEN)u[j], a, pp);
-  (void)sort_factor(y, cmp_pol); return y;
+  return to_Fq_fact((GEN)t,ex,nbf, 1, unfq,av);
 }
 /* See also: Isomorphisms between finite field and relative
  * factorization in polarit3.c */
@@ -2258,7 +2473,8 @@ GEN zrhqr(GEN a,long PREC);
 GEN
 rootsold(GEN x, long l)
 {
-  long av1=avma,i,j,f,g,gg,fr,deg,l0,l1,l2,l3,l4,ln;
+  long i, j, f, g, gg, fr, deg, ln;
+  gpmem_t av=avma, av0, av1, av2, av3;
   long exc,expmin,m,deg0,k,ti,h,ii,e,e1,emax,v;
   GEN y,xc,xd0,xd,xdabs,p1,p2,p3,p4,p5,p6,p7,p8;
   GEN p9,p10,p11,p12,p14,p15,pa,pax,pb,pp,pq,ps;
@@ -2283,7 +2499,7 @@ rootsold(GEN x, long l)
       if (gsigne(p2)>0) g=0;
     } else if (ti != t_INT && ti != t_INTMOD && !is_frac_t(ti)) g=0;
   }
-  l1=avma; p2=cgetg(3,t_COMPLEX);
+  av1=avma; p2=cgetg(3,t_COMPLEX);
   p2[1]=lmppi(DEFAULTPREC); p2[2]=ldivrs((GEN)p2[1],10);
   p11=cgetg(4,t_POL); p11[1]=evalsigne(1)+evallgef(4);
   setvarn(p11,v); p11[3]=un;
@@ -2317,22 +2533,22 @@ rootsold(GEN x, long l)
       deg=degpol(ps);
       if (deg)
       {
-        l3=avma; e=gexpo((GEN)ps[deg+2]); emax=e;
+        av3=avma; e=gexpo((GEN)ps[deg+2]); emax=e;
         for (i=2; i<deg+2; i++)
         {
           p3=(GEN)(ps[i]);
           e1=gexpo(p3); if (e1>emax) emax=e1;
         }
-        e=emax-e; if (e<0) e=0; avma=l3; if (ps!=pax) xd0=deriv(ps,v);
+        e=emax-e; if (e<0) e=0; avma=av3; if (ps!=pax) xd0=deriv(ps,v);
         xdabs=cgetg(deg+2,t_POL); xdabs[1]=xd0[1];
         for (i=2; i<deg+2; i++)
         {
-          l3=avma; p3=(GEN)xd0[i];
+          av3=avma; p3=(GEN)xd0[i];
           p4=gabs(greal(p3),l);
-          p5=gabs(gimag(p3),l); l4=avma;
-          xdabs[i]=lpile(l3,l4,gadd(p4,p5));
+          p5=gabs(gimag(p3),l);
+          xdabs[i]=lpileupto(av3, gadd(p4,p5));
         }
-        l0=avma; xc=gcopy(ps); xd=gcopy(xd0); l2=avma;
+        av0=avma; xc=gcopy(ps); xd=gcopy(xd0); av2=avma;
         for (i=1; i<=deg; i++)
         {
           if (i==deg)
@@ -2346,10 +2562,10 @@ rootsold(GEN x, long l)
             while (exc >= -20)
             {
               p7 = gneg_i(gdiv(p4, poleval(xd,p3)));
-              l3 = avma;
+              av3 = avma;
               if (gcmp0(p5)) exc = -32;
               else exc = expo(gnorm(p7))-expo(gnorm(p3));
-              avma = l3;
+              avma = av3;
               for (j=1; j<=10; j++)
               {
                 p8=gadd(p3,p7); p9=poleval(xc,p8); p10=gnorm(p9);
@@ -2358,14 +2574,14 @@ rootsold(GEN x, long l)
                   GEN *gptr[3];
                   p3=p8; p4=p9; p5=p10;
                   gptr[0]=&p3; gptr[1]=&p4; gptr[2]=&p5;
-                  gerepilemanysp(l2,l3,gptr,3);
+                  gerepilemanysp(av2,av3,gptr,3);
                   break;
                 }
-                gshiftz(p7,-2,p7); avma=l3;
+                gshiftz(p7,-2,p7); avma=av3;
               }
               if (j > 10)
               {
-                avma=av1;
+                avma=av;
                 if (DEBUGLEVEL)
                 {
                   fprintferr("too many iterations in rootsold(): ");
@@ -2375,7 +2591,7 @@ rootsold(GEN x, long l)
               }
             }
             p1=(GEN)y[k+m*i]; setlg(p1[1],3); setlg(p1[2],3); gaffect(p3,p1);
-            avma=l2; p14=(GEN)p1[1]; p15=(GEN)p1[2];
+            avma=av2; p14=(GEN)p1[1]; p15=(GEN)p1[2];
             for (ln=4; ln<=l; ln=(ln<<1)-2)
             {
               setlg(p14,ln); setlg(p15,ln);
@@ -2384,7 +2600,7 @@ rootsold(GEN x, long l)
               p4=poleval(xc,p1);
               p5=poleval(xd,p1); p6=gneg_i(gdiv(p4,p5));
               settyp(p14,t_REAL); settyp(p15,t_REAL);
-              gaffect(gadd(p1,p6),p1); avma=l2;
+              gaffect(gadd(p1,p6),p1); avma=av2;
             }
           }
           setlg(p14,l); setlg(p15,l);
@@ -2404,7 +2620,7 @@ rootsold(GEN x, long l)
           p6=gdiv(p4,poleval(xdabs,gabs(p7,l)));
           if (gexpo(p6)>=expmin)
           {
-            avma=av1;
+            avma=av;
             if (DEBUGLEVEL)
             {
               fprintferr("internal error in rootsold(): using roots2()\n");
@@ -2412,13 +2628,13 @@ rootsold(GEN x, long l)
             }
             return roots2(x,l);
           }
-          avma=l2;
+          avma=av2;
           if (expo(p1[2])<expmin && g)
           {
             gaffect(gzero,(GEN)p1[2]);
             for (j=1; j<m; j++) gaffect(p1,(GEN)y[k+(i-1)*m+j]);
             p11[2]=lneg((GEN)p1[1]);
-            l4=avma; xc=gerepile(l0,l4,gdeuc(xc,p11));
+            xc=gerepileupto(av0, gdeuc(xc,p11));
           }
           else
           {
@@ -2428,16 +2644,16 @@ rootsold(GEN x, long l)
               p1=gconj(p1);
               for (j=1; j<=m; j++) gaffect(p1,(GEN)y[k+i*m+j]);
               i++;
-              p12[2]=lnorm(p1); p12[3]=lmulsg(-2,(GEN)p1[1]); l4=avma;
-              xc=gerepile(l0,l4,gdeuc(xc,p12));
+              p12[2]=lnorm(p1); p12[3]=lmulsg(-2,(GEN)p1[1]);
+              xc=gerepileupto(av0, gdeuc(xc,p12));
             }
             else
             {
-              p11[2]=lneg(p1); l4=avma;
-              xc=gerepile(l0,l4,gdeuc(xc,p11));
+              p11[2]=lneg(p1);
+              xc=gerepileupto(av0, gdeuc(xc,p11));
             }
           }
-          xd=deriv(xc,v); l2=avma;
+          xd=deriv(xc,v); av2=avma;
         }
         k += deg*m;
       }
@@ -2445,7 +2661,7 @@ rootsold(GEN x, long l)
     }
     while (k!=deg0);
   }
-  avma=l1;
+  avma=av1;
   for (j=2; j<=deg0; j++)
   {
     p1 = (GEN)y[j];
@@ -2466,9 +2682,10 @@ rootsold(GEN x, long l)
 GEN
 roots2(GEN pol,long PREC)
 {
-  ulong av = avma;
+  gpmem_t av = avma;
   long N,flagexactpol,flagrealpol,flagrealrac,ti,i,j;
-  long nbpol,k,av1,multiqol,deg,nbroot,fr,f;
+  long nbpol, k, multiqol, deg, nbroot, fr, f;
+  gpmem_t av1;
   GEN p1,p2,rr,EPS,qol,qolbis,x,b,c,*ad,v,tabqol;
 
   if (typ(pol)!=t_POL) err(typeer,"roots2");
@@ -2481,7 +2698,7 @@ roots2(GEN pol,long PREC)
     p2=gneg_i(gdiv((GEN)pol[2],p1));
     return gerepilecopy(av,p2);
   }
-  EPS=realun(3); setexpo(EPS, 12 - bit_accuracy(PREC));
+  EPS = real2n(12 - bit_accuracy(PREC), 3);
   flagrealpol=1; flagexactpol=1;
   for (i=2; i<=N+2; i++)
   {
@@ -2611,19 +2828,23 @@ roots2(GEN pol,long PREC)
 static GEN
 laguer(GEN pol,long N,GEN y0,GEN EPS,long PREC)
 {
-  long av = avma, av1,MAXIT,iter,i,j;
+  long MAXIT, iter, j;
+  gpmem_t av = avma, av1;
   GEN rac,erre,I,x,abx,abp,abm,dx,x1,b,d,f,g,h,sq,gp,gm,g2,*ffrac;
 
-  MAXIT=MR*MT; rac=cgetg(3,t_COMPLEX);
-  rac[1]=lgetr(PREC); rac[2]=lgetr(PREC);
+  MAXIT = MR*MT; rac = cgetc(PREC);
   av1 = avma;
   I=cgetg(3,t_COMPLEX); I[1]=un; I[2]=un;
-  ffrac=(GEN*)new_chunk(MR+1); for (i=0; i<=MR; i++) ffrac[i]=cgetr(PREC);
-  affrr(dbltor(0.0), ffrac[0]); affrr(dbltor(0.5), ffrac[1]);
-  affrr(dbltor(0.25),ffrac[2]); affrr(dbltor(0.75),ffrac[3]);
-  affrr(dbltor(0.13),ffrac[4]); affrr(dbltor(0.38),ffrac[5]);
-  affrr(dbltor(0.62),ffrac[6]); affrr(dbltor(0.88),ffrac[7]);
-  affrr(dbltor(1.0),ffrac[8]);
+  ffrac = (GEN*)new_chunk(MR+1);
+  ffrac[0] = dbltor(0.0);
+  ffrac[1] = dbltor(0.5);
+  ffrac[2] = dbltor(0.25);
+  ffrac[3] = dbltor(0.75);
+  ffrac[4] = dbltor(0.13);
+  ffrac[5] = dbltor(0.38);
+  ffrac[6] = dbltor(0.62);
+  ffrac[7] = dbltor(0.88);
+  ffrac[8] = dbltor(1.0);
   x=y0;
   for (iter=1; iter<=MAXIT; iter++)
   {
@@ -2631,11 +2852,12 @@ laguer(GEN pol,long N,GEN y0,GEN EPS,long PREC)
     d=gzero; f=gzero; abx=QuickNormL1(x,PREC);
     for (j=N-1; j>=0; j--)
     {
-      f=gadd(gmul(x,f),d); d=gadd(gmul(x,d),b);
-      b=gadd(gmul(x,b),(GEN)pol[j+2]);
-      erre=gadd(QuickNormL1(b,PREC),gmul(abx,erre));
+      f = gadd(gmul(x,f),d);
+      d = gadd(gmul(x,d),b);
+      b = gadd(gmul(x,b), (GEN)pol[j+2]);
+      erre = gadd(QuickNormL1(b,PREC), gmul(abx,erre));
     }
-    erre=gmul(erre,EPS);
+    erre = gmul(erre,EPS);
     if (gcmp(QuickNormL1(b,PREC),erre)<=0)
     {
       gaffect(x,rac); avma = av1; return rac;
@@ -2674,7 +2896,7 @@ laguer(GEN pol,long N,GEN y0,GEN EPS,long PREC)
 static GEN
 balanc(GEN x)
 {
-  ulong av = avma;
+  gpmem_t av = avma;
   long last,i,j, sqrdx = (RADIX<<1), n = lg(x);
   GEN r,c,cofgen,a;
 
@@ -2728,7 +2950,7 @@ hqr(GEN mat) /* find all the eigenvalues of the matrix
 
   anorm=fabs(a[1][1]);
   for (i=2; i<=n; i++) for (j=(i-1); j<=n; j++) anorm+=fabs(a[i][j]);
-  nn=n; t=0.0;
+  nn=n; t=0.;
   if (DEBUGLEVEL>3) { fprintferr("* Finding eigenvalues\n"); flusherr(); }
   while (nn>=1)
   {
@@ -2737,11 +2959,11 @@ hqr(GEN mat) /* find all the eigenvalues of the matrix
     {
       for (l=nn; l>=2; l--)
       {
-        s=fabs(a[l-1][l-1])+fabs(a[l][l]); if (s==0.0) s=anorm;
+        s=fabs(a[l-1][l-1])+fabs(a[l][l]); if (s==0.) s=anorm;
         if ((double)(fabs(a[l][l-1])+s)==s) break;
       }
       x=a[nn][nn];
-      if (l==nn){ wr[nn]=x+t; wi[nn--]=0.0; }
+      if (l==nn){ wr[nn]=x+t; wi[nn--]=0.; }
       else
       {
         y=a[nn-1][nn-1];
@@ -2749,18 +2971,18 @@ hqr(GEN mat) /* find all the eigenvalues of the matrix
         if (l == nn-1)
         {
           p=0.5*(y-x); q=p*p+w; z=sqrt(fabs(q)); x+=t;
-          if (q>=0.0 || fabs(q)<=eps)
+          if (q>=0. || fabs(q)<=eps)
           {
             z=p+SIGN(z,p); wr[nn-1]=wr[nn]=x+z;
             if (fabs(z)>eps) wr[nn]=x-w/z;
-            wi[nn-1]=wi[nn]=0.0;
+            wi[nn-1]=wi[nn]=0.;
           }
           else{ wr[nn-1]=wr[nn]=x+p; wi[nn-1]=-(wi[nn]=z); }
           nn-=2;
         }
         else
         {
-          p = q = r = 0.0; /* for lint */
+          p = q = r = 0.; /* for lint */
           if (its==30) err(talker,"too many iterations in hqr");
           if (its==10 || its==20)
           {
@@ -2780,18 +3002,18 @@ hqr(GEN mat) /* find all the eigenvalues of the matrix
             v=fabs(p)*(fabs(a[m-1][m-1])+fabs(z)+fabs(a[m+1][m+1]));
             if ((double)(u+v)==v) break;
           }
-          for (i=m+2; i<=nn; i++){ a[i][i-2]=0.0; if (i!=(m+2)) a[i][i-3]=0.0; }
+          for (i=m+2; i<=nn; i++){ a[i][i-2]=0.; if (i!=(m+2)) a[i][i-3]=0.; }
           for (k=m; k<=nn-1; k++)
           {
             if (k!=m)
             {
               p=a[k][k-1]; q=a[k+1][k-1];
-              r = (k != nn-1)? a[k+2][k-1]: 0.0;
+              r = (k != nn-1)? a[k+2][k-1]: 0.;
               x = fabs(p)+fabs(q)+fabs(r);
-              if (x != 0.0) { p/=x; q/=x; r/=x; }
+              if (x != 0.) { p/=x; q/=x; r/=x; }
             }
             s = SIGN(sqrt(p*p+q*q+r*r),p);
-            if (s == 0.0) continue;
+            if (s == 0.) continue;
 
             if (k==m)
               { if (l!=m) a[k][k-1] = -a[k][k-1]; }
@@ -2819,10 +3041,10 @@ hqr(GEN mat) /* find all the eigenvalues of the matrix
   }
   for (j=2; j<=n; j++) /* ordering the roots */
   {
-    x=wr[j]; y=wi[j]; if (y) flj=1; else flj=0;
+    x=wr[j]; y=wi[j]; if (y != 0.) flj=1; else flj=0;
     for (k=j-1; k>=1; k--)
     {
-      if (wi[k]) flk=1; else flk=0;
+      if (wi[k] != 0.) flk=1; else flk=0;
       if (flk<flj) break;
       if (!flk && !flj && wr[k]<=x) break;
       if (flk&&flj&& wr[k]<x) break;
@@ -2835,7 +3057,7 @@ hqr(GEN mat) /* find all the eigenvalues of the matrix
   for (i=1; i<=n; i++) free(a[i]); free(a); eig=cgetg(n+1,t_COL);
   for (i=1; i<=n; i++)
   {
-    if (wi[i])
+    if (wi[i] != 0.)
     {
       GEN p1 = cgetg(3,t_COMPLEX);
       eig[i]=(long)p1;
@@ -2854,14 +3076,14 @@ hqr(GEN mat) /* find all the eigenvalues of the matrix
 GEN
 zrhqr(GEN a,long prec)
 {
-  ulong av = avma;
+  gpmem_t av = avma;
   long i,j,prec2, n = degpol(a), ex = -bit_accuracy(prec);
   GEN aa,b,p1,rt,rr,hess,x,dx,y,newval,oldval;
 
-  hess = cgetg(n+1,t_MAT); 
+  hess = cgetg(n+1,t_MAT);
   for (j=1; j<=n; j++)
   {
-    p1 = cgetg(n+1,t_COL); hess[j] = (long)p1; 
+    p1 = cgetg(n+1,t_COL); hess[j] = (long)p1;
     p1[1] = lneg(gdiv((GEN)a[n-j+2],(GEN)a[n+2]));
     for (i=2; i<=n; i++) p1[i] = (i==(j+1))? un: zero;
   }