Index: src/basemath/alglin1.c
===================================================================
RCS file: /home/cvs/pari/src/basemath/alglin1.c,v
retrieving revision 1.130
diff -u -r1.130 alglin1.c
--- src/basemath/alglin1.c	27 Aug 2003 15:25:09 -0000	1.130
+++ src/basemath/alglin1.c	26 Dec 2003 21:28:56 -0000
@@ -28,13 +28,6 @@
 #  define MASK (0xffff0000UL)
 #endif
 
-/* 2p^2 < 2^BITS_IN_LONG */
-#ifdef LONG_IS_64BIT
-#  define u_OK_ULONG(p) ((ulong)p <= 3037000493UL)
-#else
-#  define u_OK_ULONG(p) ((ulong)p <= 46337UL)
-#endif
-#define OK_ULONG(p) (lgefint(p) == 3 && u_OK_ULONG(p[2]))
 extern GEN ZM_init_CRT(GEN Hp, ulong p);
 extern int ZM_incremental_CRT(GEN H, GEN Hp, GEN q, GEN qp, ulong p);
 
@@ -1139,16 +1132,6 @@
   return u;
 }
 
-#if 0
-/* assume -p < a < p, return 1/a mod p */
-static long
-u_Fp_inv(long a, long p)
-{
-  if (a < 0) a = p + a; /* pb with ulongs < 0 */
-  return (long)invumod((ulong)a,(ulong)p);
-}
-#endif
-
 typedef ulong* uGEN;
 #define ucoeff(a,i,j)  ( ( (uGEN) ( (GEN) (a))[j]) [i])
 
@@ -1336,7 +1319,7 @@
 
 /* destroy a, b */
 static GEN
-u_FpM_gauss_sp(GEN a, GEN b, ulong p)
+Flm_gauss_sp(GEN a, GEN b, ulong p)
 {
   long iscol, i, j, k, li, bco, aco = lg(a)-1;
   ulong piv, invpiv, m;
@@ -1414,16 +1397,16 @@
 }
 
 GEN
-u_FpM_gauss(GEN a, GEN b, ulong p) {
-  return u_FpM_gauss_sp(dummycopy(a), dummycopy(b), p);
+Flm_gauss(GEN a, GEN b, ulong p) {
+  return Flm_gauss_sp(dummycopy(a), dummycopy(b), p);
 }
 static GEN
-u_FpM_inv_sp(GEN a, ulong p) {
-  return u_FpM_gauss_sp(a, u_idmat(lg(a)-1), p);
+Flm_inv_sp(GEN a, ulong p) {
+  return Flm_gauss_sp(a, u_idmat(lg(a)-1), p);
 }
 GEN
-u_FpM_inv(GEN a, ulong p) {
-  return u_FpM_inv_sp(dummycopy(a), p);
+Flm_inv(GEN a, ulong p) {
+  return Flm_inv_sp(dummycopy(a), p);
 }
 
 GEN
@@ -1448,9 +1431,9 @@
   if (OK_ULONG(p))
   {
     ulong pp = (ulong)p[2];
-    a = u_Fp_FpM(a, pp);
-    b = u_Fp_FpM(b, pp);
-    u = u_FpM_gauss_sp(a,b, pp);
+    a = ZM_Flm(a, pp);
+    b = ZM_Flm(b, pp);
+    u = Flm_gauss_sp(a,b, pp);
     u = iscol? small_to_col((GEN)u[1]): small_to_mat(u);
     return gerepileupto(av, u);
   }
@@ -1581,8 +1564,8 @@
 FpM_inv(GEN x, GEN p) { return FpM_gauss(x, NULL, p); }
 
 /* set y = x * y */
-static GEN
-u_FpM_Fp_mul_ip(GEN y, ulong x, ulong p)
+GEN
+Flm_Fl_mul_inplace(GEN y, ulong x, ulong p)
 {
   long i, j, m = lg(y[1]), l = lg(y);
   if (HIGHWORD(x | p))
@@ -1617,8 +1600,8 @@
     NEXT_PRIME_VIADIFF_CHECK(p,d);
     dMp = umodiu(dM,p);
     if (!dMp) continue;
-    Hp = u_FpM_inv_sp(u_Fp_FpM(M,p), p);
-    if (dMp != 1) Hp = u_FpM_Fp_mul_ip(Hp, dMp, p);
+    Hp = Flm_inv_sp(ZM_Flm(M,p), p);
+    if (dMp != 1) Hp = Flm_Fl_mul_inplace(Hp, dMp, p);
 
     if (!H)
     {
@@ -2507,9 +2490,33 @@
   return z;
 }
 
+/*FIXME: Does it assume OK_ULONG ?*/
+
+GEN
+Flm_Flv_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++)
+  {
+    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;
+}
+
 /* in place, destroy x */
-static GEN
-u_FpM_ker_sp(GEN x, ulong p, long deplin)
+GEN
+Flm_ker_sp(GEN x, ulong p, long deplin)
 {
   GEN y, c, d;
   long i, j, k, r, t, m, n;
@@ -2565,7 +2572,7 @@
   }
   if (deplin) return NULL;
 
-  y = cgetg(r+1, t_MAT);
+  y = cgetg(r+1, t_VEC);
   for (j=k=1; j<=r; j++,k++)
   {
     GEN C = cgetg(n+1, t_VECSMALL);
@@ -2594,8 +2601,8 @@
   if (OK_ULONG(p))
   {
     ulong pp = (ulong)p[2];
-    y = u_Fp_FpM(x, pp);
-    y = u_FpM_ker_sp(y, pp, deplin);
+    y = ZM_Flm(x, pp);
+    y = Flm_ker_sp(y, pp, deplin);
     if (!y) return y;
     y = deplin? small_to_col(y): small_to_mat(y);
     return gerepileupto(av0, y);
@@ -2687,9 +2694,9 @@
 FpM_deplin(GEN x, GEN p) { return FpM_ker_i(x, p, 1); }
 /* not memory clean */
 GEN
-u_FpM_ker(GEN x, ulong p) { return u_FpM_ker_sp(dummycopy(x), p, 0); }
+Flm_ker(GEN x, ulong p) { return Flm_ker_sp(gcopy(x), p, 0); }
 GEN
-u_FpM_deplin(GEN x, ulong p) { return u_FpM_ker_sp(dummycopy(x), p, 1); }
+Flm_deplin(GEN x, ulong p) { return Flm_ker_sp(gcopy(x), p, 1); }
 
 static void
 FpM_gauss_pivot(GEN x, GEN p, GEN *dd, long *rr)
@@ -2862,71 +2869,6 @@
 /**************************************************************
  Rather stupid implementation of linear algebra in finite fields
  **************************************************************/
-static GEN
-Fq_add(GEN x, GEN y, GEN T/*unused*/, GEN p)
-{
-  (void)T;
-  switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
-  {
-    case 0: return modii(addii(x,y),p);
-    case 1: return FpX_Fp_add(x,y,p);
-    case 2: return FpX_Fp_add(y,x,p);
-    case 3: return FpX_add(x,y,p);
-  }
-  return NULL;
-}
-#if 0
-/* this function is really for FpV_roots_to_pol in polarit1.c
- * For consistency we write the code there.
- * To avoid having to remove static status, we rewrite it in polarit1.c
- */
-static GEN
-Fq_neg(GEN x, GEN T, GEN p)
-{
-  switch(typ(x)==t_POL)
-  {
-    case 0: return signe(x)?subii(p,x):gzero;
-    case 1: return FpX_neg(x,p);
-  }
-  return NULL;
-}
-#endif
-
-GEN
-Fq_mul(GEN x, GEN y, GEN T, GEN p)
-{
-  switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
-  {
-    case 0: return modii(mulii(x,y),p);
-    case 1: return FpX_Fp_mul(x,y,p);
-    case 2: return FpX_Fp_mul(y,x,p);
-    case 3: return FpXQ_mul(x,y,T,p);
-  }
-  return NULL;
-}
-
-static GEN
-Fq_neg_inv(GEN x, GEN T, GEN p)
-{
-  switch(typ(x)==t_POL)
-  {
-    case 0: return mpinvmod(negi(x),p);
-    case 1: return FpXQ_inv(FpX_neg(x,p),T,p);
-  }
-  return NULL;
-}
-
-static GEN
-Fq_res(GEN x, GEN T, GEN p)
-{
-  pari_sp ltop=avma;
-  switch(typ(x)==t_POL)
-  {
-    case 0: return modii(x,p);
-    case 1: return gerepileupto(ltop,FpX_res(FpX_red(x,p),T,p));
-  }
-  return NULL;
-}
 
 static void
 Fq_gerepile_gauss_ker(GEN x, GEN T, GEN p, long k, long t, pari_sp av)
Index: src/basemath/alglin2.c
===================================================================
RCS file: /home/cvs/pari/src/basemath/alglin2.c,v
retrieving revision 1.111
diff -u -r1.111 alglin2.c
--- src/basemath/alglin2.c	19 Dec 2003 11:05:36 -0000	1.111
+++ src/basemath/alglin2.c	26 Dec 2003 21:29:01 -0000
@@ -2815,7 +2815,7 @@
   GEN U, perm, y = cgetg(4, t_VEC);
   y[1] = (long)hnfperm_i(A, &U, &perm);
   y[2] = (long)U;
-  y[3] = (long)small_to_vec(perm); return y;
+  y[3] = (long)vecsmall_vec(perm); return y;
 }
 
 /* Hermite Normal Form */
Index: src/basemath/arith1.c
===================================================================
RCS file: /home/cvs/pari/src/basemath/arith1.c,v
retrieving revision 1.122
diff -u -r1.122 arith1.c
--- src/basemath/arith1.c	19 Dec 2003 11:05:36 -0000	1.122
+++ src/basemath/arith1.c	26 Dec 2003 21:29:05 -0000
@@ -1345,6 +1345,15 @@
   return res;
 }
 
+GEN
+mpinvmodsafe(GEN a, GEN m)
+{
+  GEN res;
+  if (! invmod(a,m,&res))
+    return NULL;
+  return res;
+}
+
 /*********************************************************************/
 /**                                                                 **/
 /**                    MODULAR EXPONENTIATION                       **/
Index: src/basemath/base2.c
===================================================================
RCS file: /home/cvs/pari/src/basemath/base2.c,v
retrieving revision 1.234
diff -u -r1.234 base2.c
--- src/basemath/base2.c	19 Dec 2003 11:05:36 -0000	1.234
+++ src/basemath/base2.c	26 Dec 2003 21:29:20 -0000
@@ -21,7 +21,6 @@
 #include "parinf.h"
 
 #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 addshiftw(GEN x, GEN y, long d);
 extern GEN norm_by_embed(long r1, GEN x);
 extern GEN ZX_resultant_all(GEN A, GEN B, GEN dB, ulong bound);
@@ -31,7 +30,7 @@
 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 FqV_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);
@@ -2267,7 +2266,7 @@
     if (typ(w) != t_INT)
     {
       w = Q_primitive_part(w, &cx);
-      w = FpX_rem(w, T, p);
+      w = FpX_res(w, T, p);
       if (cx)
       {
         cx = gmod(cx, p);
@@ -2639,7 +2638,7 @@
   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);
+  return FqV_red(z,D->T,D->p);
 }
 static GEN
 _sqr(void *data, GEN x)
@@ -2647,7 +2646,7 @@
   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);
+  return FqV_red(z,D->T,D->p);
 }
 
 /* Compute W[h]^n mod pr in the extension, assume n >= 0 */
Index: src/basemath/base3.c
===================================================================
RCS file: /home/cvs/pari/src/basemath/base3.c,v
retrieving revision 1.113
diff -u -r1.113 base3.c
--- src/basemath/base3.c	19 Dec 2003 11:05:36 -0000	1.113
+++ src/basemath/base3.c	26 Dec 2003 21:29:22 -0000
@@ -21,8 +21,6 @@
 #include "pari.h"
 #include "parinf.h"
 extern GEN to_famat(GEN g, GEN e);
-extern GEN u_FpM_deplin(GEN x, ulong p);
-extern GEN u_FpM_inv(GEN a, ulong p);
 extern GEN famat_makecoprime(GEN nf, GEN g, GEN e, GEN pr, GEN prk, GEN EX);
 extern GEN famat_to_nf_modidele(GEN nf, GEN g, GEN e, GEN bid);
 extern GEN vconcat(GEN A, GEN B);
@@ -1559,7 +1557,7 @@
       c = (GEN)mat[lgmat];
       for (i=1; i<=nba; i++)
         c[i] = (gsigne((GEN)alpha[i]) < 0)? 1: 0;
-      if (u_FpM_deplin(mat, 2)) { avma = av1; continue; }
+      if (Flm_deplin(mat, 2)) { avma = av1; continue; }
 
       /* new vector indep. of previous ones */
       avma = av1; alpha = nfun;
@@ -1569,7 +1567,7 @@
       genarch[lgmat++] = (long)alpha;
       if (lgmat > nba)
       {
-        mat = small_to_mat( u_FpM_inv(mat, 2) );
+        mat = small_to_mat( Flm_inv(mat, 2) );
         gerepileall(av,2,&genarch,&mat);
         y[2] = (long)genarch;
         y[3] = (long)mat; return y;
Index: src/basemath/buch4.c
===================================================================
RCS file: /home/cvs/pari/src/basemath/buch4.c,v
retrieving revision 1.43
diff -u -r1.43 buch4.c
--- src/basemath/buch4.c	19 Dec 2003 11:05:36 -0000	1.43
+++ src/basemath/buch4.c	26 Dec 2003 21:29:24 -0000
@@ -414,7 +414,10 @@
   ordp= subis( p, 1);                 /* |F_p^*|        */
   modpr = nf_to_ff_init(nf, &pr,&T,&p);
   t = nf_to_ff(nf,t,modpr);
-  t = FpXQ_pow(t, diviiexact(ord, ordp), T,p); /* in F_p^* */
+  if (T)
+    t = FpXQ_pow(t, diviiexact(ord, ordp), T,p); /* in F_p^* */
+  else
+    t = powmodulo(t, diviiexact(ord, ordp), p);
   if (typ(t) == t_POL)
   {
     if (degpol(t)) err(bugparier,"nfhilbertp");
Index: src/basemath/gen3.c
===================================================================
RCS file: /home/cvs/pari/src/basemath/gen3.c,v
retrieving revision 1.111
diff -u -r1.111 gen3.c
--- src/basemath/gen3.c	19 Dec 2003 11:05:36 -0000	1.111
+++ src/basemath/gen3.c	26 Dec 2003 21:29:30 -0000
@@ -2519,7 +2519,7 @@
     for (i=1; i<=lx-2; i++) y[i]=lcopy((GEN)x[i]);
     return y;
   }
-  if (tx == t_VECSMALL) return small_to_vec(x);
+  if (tx == t_VECSMALL) return vecsmall_vec(x);
   if (!signe(x)) { y=cgetg(2,t_VEC); y[1]=zero; return y; }
   lx=lg(x); y=cgetg(lx-1,t_VEC); x++;
   for (i=1; i<=lx-2; i++) y[i]=lcopy((GEN)x[i]);
Index: src/basemath/perm.c
===================================================================
RCS file: /home/cvs/pari/src/basemath/perm.c,v
retrieving revision 1.19
diff -u -r1.19 perm.c
--- src/basemath/perm.c	19 Dec 2003 10:36:26 -0000	1.19
+++ src/basemath/perm.c	26 Dec 2003 21:29:33 -0000
@@ -21,6 +21,25 @@
 /**                                                                     **/
 /*************************************************************************/
 
+GEN
+vecsmall_vec(GEN z)
+{
+  long i, l = lg(z);
+  GEN x = cgetg(l,t_VEC);
+  for (i=1; i<l; i++) x[i] = lstoi(z[i]);
+  return x;
+}
+
+GEN
+vecsmall_copy(GEN x)
+{
+  long i, l = lg(x);
+  GEN z = cgetg(l, t_VECSMALL);
+  if (typ(x)!=t_VECSMALL) err(typeer,"vecsmall_copy");
+  for (i=1; i<l; i++) z[i] = x[i];
+  return z;
+}
+
 GEN vecsmall_const(long n, long c)
 {
   long i;
@@ -62,7 +81,7 @@
   pari_sp ltop=avma;
   GEN W;
   long i,j;
-  if ( lg(V) == 1 ) return gcopy(V);
+  if ( lg(V) == 1 ) return vecsmall_copy(V);
   W=cgetg(lg(V),t_VECSMALL);
   W[1]=V[1];
   for(i=2,j=1;i<lg(V);i++)
@@ -498,7 +517,7 @@
   long n = lg(H)-1;
   GEN L = cgetg(1+n*o, t_VEC);
   for(i=1; i<=n; i++)
-    L[i]=lcopy((GEN)H[i]);
+    L[i]= (long) vecsmall_copy((GEN)H[i]);
   for(k=n+1; k <= n*o; ++k)
     L[k] = (long) perm_mul((GEN) L[k-n], S);
   return L;
@@ -535,7 +554,7 @@
   GEN ord=(GEN) G[2];
   long card=group_order(G);
   res = cgetg(card + 1, t_VEC);
-  res[1] = lcopy(g);
+  res[1] = (long) vecsmall_copy(g);
   k = 1;
   for (i = 1; i < lg(gen); i++)
   {
@@ -556,7 +575,7 @@
   GEN ord=(GEN) G[2];
   long card=group_order(G);
   res = cgetg(card + 1, t_VEC);
-  res[1] = lcopy(g);
+  res[1] = (long) vecsmall_copy(g);
   k = 1;
   for (i = 1; i < lg(gen); i++)
   {
@@ -582,7 +601,7 @@
   GEN p2,p3,p4;
   p2 = cgetg(3, t_VEC);
   p3 = cgetg(2, t_VEC);
-  p3[1] = lcopy(g);
+  p3[1] = (long) vecsmall_copy(g);
   p4 = cgetg(2,t_VECSMALL);
   p4[1] = s;
   p2[1] = (long) p3;
@@ -597,8 +616,8 @@
   GEN H = cgetg(3, t_VEC);
   GEN p3,p4;
   p3 = cgetg(3, t_VEC);
-  p3[1] = lcopy((GEN)g1);
-  p3[2] = lcopy((GEN)g2);
+  p3[1] = (long) vecsmall_copy((GEN)g1);
+  p3[2] = (long) vecsmall_copy((GEN)g2);
   p4 = cgetg(3,t_VECSMALL);
   p4[1] = s1;
   p4[2] = s2;
@@ -642,7 +661,7 @@
   }
   p1 = cgetg(3,t_VEC);
   p1[1] = lcopy(p2);
-  p1[2]= (long) vecvecsmall_sort(p3);
+  p1[2] = (long) vecvecsmall_sort(p3);
   return gerepileupto(ltop,p1);
 }
 
@@ -940,7 +959,7 @@
   long i;
   long d=1;
   G[1]=lgetg(lg(v),t_VEC);
-  G[2]=lcopy(v);
+  G[2]=(long) vecsmall_copy(v);
   card=group_order(G);
   for(i=1;i<lg(v);i++)
   {
@@ -986,7 +1005,7 @@
   V=cgetg(l+1,t_VEC);
   for (i=1, j=1; i<=n ;i++)
     if (!bitvec_test(elts,i))
-      V[j++]=lcopy((GEN)S[i]);
+      V[j++]=(long) vecsmall_copy((GEN)S[i]);
   return gerepileupto(ltop,V);
 }
 
@@ -1055,7 +1074,7 @@
   {
     if (i > 1)
       s = concat(s, STRtoGENstr(", "));
-    s = concat(s, small_to_vec(gmael(G,1,i)));
+    s = concat(s, vecsmall_vec(gmael(G,1,i)));
   }
   s = concat(s, STRtoGENstr(">"));
   return gerepileupto(ltop,s);
Index: src/basemath/polarit1.c
===================================================================
RCS file: /home/cvs/pari/src/basemath/polarit1.c,v
retrieving revision 1.144
diff -u -r1.144 polarit1.c
--- src/basemath/polarit1.c	19 Dec 2003 11:05:36 -0000	1.144
+++ src/basemath/polarit1.c	26 Dec 2003 21:29:45 -0000
@@ -478,17 +478,15 @@
 #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. */
 
+/*TODO: merge with matrixpow*/
 /* u in Z[X], return kernel of (Frob - Id) over Fp[X] / u */
-static GEN
-FpM_Berlekamp_ker(GEN u, GEN p)
+GEN
+FpX_Berlekamp_ker(GEN u, GEN p)
 {
   long j,N = degpol(u);
   GEN v,w,Q,p1;
@@ -501,16 +499,17 @@
     Q[j] = (long)p1;
     if (j < N)
     {
-      pari_sp av = avma;
-      w = gerepileupto(av, FpX_res(gmul(w,v), u, p));
+      pari_sp av = avma; /*FpXQ_mul is not stack clean*/
+      w = gerepileupto(av, FpXQ_mul(w, v, u, p));
     }
   }
   return FpM_ker(Q,p);
 }
 
-static GEN
-FqM_Berlekamp_ker(GEN u, GEN T, GEN q, GEN p)
+GEN
+FqX_Berlekamp_ker(GEN u, GEN T, GEN q, GEN p)
 {
+  pari_sp ltop=avma;
   long j,N = degpol(u);
   GEN v,w,Q,p1;
   Q = cgetg(N+1,t_MAT); Q[1] = (long)zerocol(N);
@@ -526,7 +525,34 @@
       w = gerepileupto(av, FpXQX_divres(FpXQX_mul(w,v, T,p), u,T,p,ONLY_REM));
     }
   }
-  return FqM_ker(Q,T,p);
+  return gerepileupto(ltop,FqM_ker(Q,T,p));
+}
+
+GEN
+Flx_Berlekamp_ker(GEN u, ulong p)
+{
+  pari_sp ltop=avma;
+  long j,N = degpol(u);
+  GEN v,w,Q,p1;
+  pari_timer T;
+  TIMER(&T);
+  Q = cgetg(N+1,t_VEC); Q[1] = (long) vecsmall_const(N,0);
+  w = v = Flxq_pow(Flx_polx(u[1]),stoi(p),u,p);
+  for (j=2; j<=N; j++)
+  {
+    p1 = Flx_Flv_lg(w, N);
+    p1[j] = subuumod((ulong)p1[j], 1, p);
+    Q[j] = (long)p1;
+    if (j < N)
+    {
+      pari_sp av = avma; /*Flxq_mul is not stack clean*/
+      w = gerepileupto(av, Flxq_mul(w, v, u, p));
+    }
+  }
+  if(DEBUGLEVEL>=9) msgTIMER(&T,"Berlekamp matrix");
+  Q=Flm_ker_sp(Q,p,0);
+  if(DEBUGLEVEL>=9) msgTIMER(&T,"kernel");
+  return gerepileupto(ltop,Q);
 }
 
 /* f in ZZ[X] and p a prime number. */
@@ -541,37 +567,47 @@
 }
 
 /* product of terms of degree 1 in factorization of f */
-GEN
+static GEN
+FpX_split_part(GEN f, GEN p)
+{
+  long n = degpol(f);
+  GEN z, X = polx[varn(f)];
+  if (n <= 1) return f;
+  f = FpX_red(f, p);
+  z = FpXQ_pow(X, p, f, p);
+  z = FpX_sub(z, X, NULL);
+  return FpX_gcd(z,f,p);
+}
+
+static GEN
 FqX_split_part(GEN f, GEN T, GEN p)
 {
   long n = degpol(f);
   GEN z, X = polx[varn(f)];
   if (n <= 1) return f;
-  if (!T)
-  {
-    f = FpX_red(f, p);
-    z = FpXQ_pow(X, p, f, p);
-    z = FpX_sub(z, X, NULL);
-    return FpX_gcd(z,f,p);
-  } else {
-    f = FqX_red(f, T, p);
-    z = FqXQ_pow(X, gpowgs(p, degpol(T)), f, T, p);
-    z = gsub(z, X);
-    return FqX_gcd(z, f, T, p);
-  }
+  f = FqX_red(f, T, p);
+  z = FqXQ_pow(X, gpowgs(p, degpol(T)), f, T, p);
+  z = gsub(z, X);
+  return FqX_gcd(z, f, T, p);
 }
 
 /* Compute the number of roots in Fq without counting multiplicity
  * return -1 for 0 polynomial. lc(f) must be prime to p. */
 long
+FpX_nbroots(GEN f, GEN p) 
+{
+  pari_sp av = avma;
+  GEN z = FpX_split_part(f, p);
+  avma = av; return degpol(z);
+}
+
+long
 FqX_nbroots(GEN f, GEN T, GEN p)
 {
   pari_sp av = avma;
   GEN z = FqX_split_part(f, T, p);
   avma = av; return degpol(z);
 }
-long
-FpX_nbroots(GEN f, GEN p) { return FqX_nbroots(f, NULL, p); }
 
 long
 FpX_is_totally_split(GEN f, GEN p)
@@ -587,13 +623,14 @@
   return degpol(z) == 1 && gcmp1((GEN)z[3]) && !signe(z[2]); /* x^p = x ? */
 }
 
-/* u_vec_to_pol(u_FpM_FpV_mul(x, u_pol_to_vec(y), p)) */
+/*FIXME: what does that do ? */
 static GEN
-u_FpM_FpX_mul(GEN x, GEN y, ulong p)
+Flm_Flx_mul(GEN x, GEN y, ulong p)
 {
   long i,k,l, ly = lg(y)-1;
   GEN z;
-  if (ly==1) return u_zeropol();
+  long vs=y[1];
+  if (ly==1) return Flx_zero(vs);
   l = lg(x[1]);
   y++;
   z = vecsmall_const(l,0) + 1;
@@ -617,16 +654,12 @@
   }
   for (i=1; i<l; i++) z[i] %= p;
   while (--l && !z[l]);
-  if (!l) return u_zeropol();
-  *z-- = evalsigne(1) | evalvarn(0); return z;
+  if (!l) return Flx_zero(0);
+  *z-- = vs; return z;
 }
 
-extern GEN u_FpX_Kmul(GEN a, GEN b, ulong p, long na, long nb);
-#define u_FpX_mul(x,y, p) u_FpX_Kmul(x+2,y+2,p, lg(x)-2,lg(y)-2)
-#define u_FpX_div(x,y,p) u_FpX_divrem((x),(y),(p),NULL)
-
 static GEN
-u_FpM_Frobenius(GEN u, ulong p)
+Flx_Frobenius(GEN u, ulong p)
 {
   long j,N = degpol(u);
   GEN v,w,Q;
@@ -636,14 +669,14 @@
   Q = cgetg(N+1,t_MAT); 
   Q[1] = (long)vecsmall_const(N, 0);
   coeff(Q,1,1) = 1;
-  w = v = u_FpXQ_pow(pol_to_small(polx[varn(u)]), utoi(p), u, p);
+  w = v = Flxq_pow(Flx_polx(u[1]), utoi(p), u, p);
   for (j=2; j<=N; j++)
   {
-    Q[j] = (long)u_pol_to_vec(w, N);
+    Q[j] = (long)Flx_Flv_lg(w, N);
     if (j < N)
     {
       pari_sp av = avma;
-      w = gerepileupto(av, u_FpX_rem(u_FpX_mul(w, v, p), u, p));
+      w = gerepileupto(av, Flxq_mul(w, v, u, p));
     }
   }
   if (DEBUGLEVEL > 7) msgTIMER(&T, "frobenius");
@@ -654,17 +687,17 @@
  * u must be squarefree mod p.
  * leading term of u must be prime to p. */
 long
-u_FpX_nbfact(GEN z, long p)
+Flx_nbfact(GEN z, long p)
 {
   long lgg, nfacp = 0, d = 0, e = degpol(z);
-  GEN g, w, MP = u_FpM_Frobenius(z, p), PolX = pol_to_small(polx[0]);
+  GEN g, w, MP = Flx_Frobenius(z, p), PolX = Flx_polx(0);
 
   w = PolX;
   while (d < (e>>1))
   { /* here e = degpol(z) */
     d++;
-    w = u_FpM_FpX_mul(MP, w, p); /* w^p mod (z,p) */
-    g = u_FpX_gcd(z, u_FpX_sub(w, PolX, p), p);
+    w = Flm_Flx_mul(MP, w, p); /* w^p mod (z,p) */
+    g = Flx_gcd(z, Flx_sub(w, PolX, p), p);
     lgg = degpol(g);
     if (!lgg) continue;
 
@@ -673,8 +706,8 @@
     if (DEBUGLEVEL>5)
       fprintferr("   %3ld fact. of degree %3ld\n", lgg/d, d);
     if (!e) break;
-    z = u_FpX_div(z, g, p);
-    w = u_FpX_rem(w, z, p);
+    z = Flx_div(z, g, p);
+    w = Flx_res(w, z, p);
   }
   if (e)
   {
@@ -685,16 +718,16 @@
 }
 
 long
-u_FpX_nbroots(GEN f, long p)
+Flx_nbroots(GEN f, long p)
 {
   long n = degpol(f);
   pari_sp av = avma;
   GEN z, X;
   if (n <= 1) return n;
-  X = pol_to_small(polx[varn(f)]);
-  z = u_FpXQ_pow(X, utoi(p), f, p);
-  z = u_FpX_sub(z, X, p);
-  z = u_FpX_gcd(z, f, p);
+  X = Flx_polx(f[1]);
+  z = Flxq_pow(X, utoi(p), f, p);
+  z = Flx_sub(z, X, p);
+  z = Flx_gcd(z, f, p);
   avma = av; return degpol(z);
 }
 
@@ -702,7 +735,7 @@
 FpX_nbfact(GEN u, GEN p)
 {
   pari_sp av = avma;
-  GEN vker = FpM_Berlekamp_ker(u, p);
+  GEN vker = FpX_Berlekamp_ker(u, p);
   avma = av; return lg(vker)-1;
 }
 
@@ -712,7 +745,7 @@
   pari_sp av = avma;
   GEN vker;
   if (!T) return FpX_nbfact(u, p);
-  vker = FqM_Berlekamp_ker(u, T, gpowgs(p, degpol(T)), p);
+  vker = FqX_Berlekamp_ker(u, T, gpowgs(p, degpol(T)), p);
   avma = av; return lg(vker)-1;
 }
 
@@ -1021,44 +1054,6 @@
   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) | 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++)
-  {
-    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
-
 /* return the (N-dimensional) vector of coeffs of p */
 GEN
 pol_to_vec(GEN x, long N)
@@ -1077,19 +1072,6 @@
   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 = lg(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)
@@ -1160,29 +1142,22 @@
   return normalizepol_i(y,ly);
 }
 
-typedef ulong* uGEN;
-
 /* set x <-- x + c*y mod p */
+/* x is not required to be normalized.*/
 static void
-u_FpX_addmul(uGEN x, uGEN y, ulong c, ulong p)
+Flx_addmul_inplace(GEN gx, GEN gy, ulong c, ulong p)
 {
-  long i, lx, ly, l;
+  long i, lx, ly;
+  ulong *x=(ulong *)gx;
+  ulong *y=(ulong *)gy;
   if (!c) return;
-  lx = lg(x);
-  ly = lg(y); l = min(lx,ly);
-  if (p & ~MAXHALFULONG)
-  {
-    for (i=2; i<l;  i++) x[i] = (x[i] + muluumod(c,y[i],p)) % p;
-    if (l == lx)
-      for (   ; i<ly; i++) x[i] = muluumod(c,y[i],p);
-  }
+  lx = lg(gx);
+  ly = lg(gy);
+  if (lx<ly) err(bugparier,"lx<ly in Flx_addmul_inplace");
+  if (u_OK_ULONG(p))
+    for (i=2; i<ly;  i++) x[i] = (x[i] + c*y[i]) % p;
   else
-  {
-    for (i=2; i<l;  i++) x[i] = (x[i] + c*y[i]) % p;
-    if (l == lx)
-      for (   ; i<ly; i++) x[i] = (c*y[i]) % p;
-  }
-  (void)u_normalizepol((GEN)x,i);
+    for (i=2; i<ly;  i++) x[i] = adduumod(x[i], muluumod(c,y[i],p),p);
 }
 
 static long
@@ -1213,24 +1188,24 @@
 #define set_irred(i) { if ((i)>ir) swap(t[i],t[ir]); ir++;}
 
 long
-FpX_split_berlekamp(GEN *t, GEN p)
+FpX_split_Berlekamp(GEN *t, GEN p)
 {
-  GEN u = *t, a,b,po2,vker,pol;
-  long N = degpol(u), d, i, ir, L, la, lb, ps, vu = varn(u);
-  pari_sp av0 = avma;
-
-  vker = FpM_Berlekamp_ker(u,p);
-  vker = mat_to_vecpol(vker,vu);
-  d = lg(vker)-1;
-  ps = is_bigint(p)? 0: p[2];
+  GEN u = *t, a,b,po2,vker;
+  long d, i, ir, L, la, lb, vu = varn(u);
+  long l = lg(u);
+  ulong ps = is_bigint(p)? 0: p[2];
   if (ps)
   {
-    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;
+    vker = Flx_Berlekamp_ker(ZX_Flx(u,ps),ps);
+    vker = Flm_FlxV(vker, u[1]);
+  }
+  else
+  {
+    vker = FpX_Berlekamp_ker(u,p);
+    vker = mat_to_vecpol(vker,vu);
   }
+  d = lg(vker)-1;
   po2 = shifti(p, -1); /* (p-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; )
@@ -1238,19 +1213,20 @@
     GEN polt;
     if (ps)
     {
-      pol[2] = small_rand(ps); /* vker[1] = 1 */
-      setlg(pol, pol[2]? 3: 2);
+      GEN pol = vecsmall_const(l-2,0);
+      pol[1] = u[1];
+      pol[2] = small_rand(ps); /*Assume vker[1]=1*/
       for (i=2; i<=d; i++)
-        u_FpX_addmul((uGEN)pol, (uGEN)vker[i], (ulong)small_rand(ps), (ulong)ps);
-      polt = small_to_pol(pol,vu);
+        Flx_addmul_inplace(pol, (GEN) vker[i], (ulong)small_rand(ps), ps);
+      (void)Flx_renormalize((GEN)pol,l-1);
+      
+      polt = Flx_ZX(pol);
     }
     else
     {
-      pol[2] = (long)genrand(p);
-      setlg(pol, signe(pol[2])? 3: 2);
-      pol[1] = u[1];
+      GEN pol = scalarpol(genrand(p), vu);
       for (i=2; i<=d; i++)
-        pol = gadd(pol, gmul((GEN)vker[i], genrand(p)));
+        pol = FpX_add(pol, FpX_Fp_mul((GEN)vker[i], genrand(p),NULL),NULL);
       polt = FpX_red(pol,p);
     }
     for (i=ir; i<L && L<d; i++)
@@ -1271,7 +1247,7 @@
       {
         pari_sp av = avma;
         b = FpX_res(polt, a, p);
-        if (degpol(b) == 0) { avma=av; continue; }
+        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)
@@ -1287,11 +1263,30 @@
   return d;
 }
 
-static GEN
+GEN
+ZX_deriv(GEN x)
+{
+  long i,lx = lg(x)-1;
+  GEN y;
+
+  if (lx<3) return zeropol(varn(x));
+  y = cgetg(lx,t_POL);
+  for (i=2; i<lx ; i++) y[i] = lmulsi(i-1,(GEN)x[i+1]);
+  y[1] = x[1]; return y;
+}
+
+GEN
+FpX_deriv(GEN f, GEN p)
+{
+  return FpX_red(ZX_deriv(f), p);
+}
+
+GEN
 FqX_deriv(GEN f, GEN T, GEN p)
 {
-  return FqX_red(derivpol(f), T, p);
+  return FqX_red(ZX_deriv(f), T, p);
 }
+
 GEN
 FqX_gcd(GEN P, GEN Q, GEN T, GEN p)
 {
@@ -1310,13 +1305,13 @@
 }
 
 long
-FqX_split_berlekamp(GEN *t, GEN q, GEN T, GEN p)
+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 = FqX_Berlekamp_ker(u,T,q,p);
   vker = mat_to_vecpol(vker,vu);
   d = lg(vker)-1;
   qo2 = shifti(q, -1); /* (q-1) / 2 */
@@ -1393,7 +1388,7 @@
       if (N > 0)
       {
         t[nbfact] = FpX_normalize(u,pp);
-        d = (N==1)? 1: FpX_split_berlekamp(t+nbfact, pp);
+        d = (N==1)? 1: FpX_split_Berlekamp(t+nbfact, pp);
         for (j=0; j<d; j++) ex[nbfact+j] = e*k;
         nbfact += d;
       }
Index: src/basemath/polarit2.c
===================================================================
RCS file: /home/cvs/pari/src/basemath/polarit2.c,v
retrieving revision 1.292
diff -u -r1.292 polarit2.c
--- src/basemath/polarit2.c	19 Dec 2003 11:05:36 -0000	1.292
+++ src/basemath/polarit2.c	26 Dec 2003 21:29:57 -0000
@@ -38,7 +38,7 @@
 extern GEN to_polmod(GEN x, GEN mod);
 extern GEN vconcat(GEN Q1, GEN Q2);
 extern int approx_0(GEN x, GEN y);
-extern long FpX_split_berlekamp(GEN *t, GEN pp);
+extern long FpX_split_Berlekamp(GEN *t, GEN pp);
 extern long u_center(ulong u, ulong p, ulong ps2);
 extern void gerepilemanycoeffs2(pari_sp av, GEN x, int n, GEN y, int o);
 
@@ -342,9 +342,18 @@
     d = (GEN)d[2];
     if (!gcmp1(d))
     {
-      d = FpXQ_inv(d, T, p);
-      u = FpXQX_FpXQ_mul(u, d, T, p);
-      v = FpXQX_FpXQ_mul(v, d, T, p);
+      if (T)
+      {
+        d = FpXQ_inv(d, T, p);
+        u = FpXQX_FpXQ_mul(u, d, T, p);
+        v = FpXQX_FpXQ_mul(v, d, T, p);
+      }
+      else
+      {
+        d = mpinvmod(d, p);
+        u = FpX_Fp_mul(u, d, p);
+        v = FpX_Fp_mul(v, d, p);
+      }
     }
     W[j]   = (long)u;
     W[j+1] = (long)v;
@@ -1491,10 +1500,10 @@
   {
     NEXT_PRIME_VIADIFF_CHECK(p,pt);
     if (lead && !smodis(lead,p)) continue;
-    z = u_Fp_FpX(a, p);
-    if (!u_FpX_is_squarefree(z, p)) continue;
+    z = ZX_Flx(a, p);
+    if (!Flx_is_squarefree(z, p)) continue;
 
-    nfacp = fl? u_FpX_nbroots(z, p): u_FpX_nbfact(z, p);
+    nfacp = fl? Flx_nbroots(z, p): Flx_nbfact(z, p);
     if (DEBUGLEVEL>4)
       fprintferr("...tried prime %3ld (%-3ld %s). Time = %ld\n",
                   p, nfacp, fl?"roots": "factors", TIMER(&T2));
@@ -1516,7 +1525,7 @@
 
   famod = cgetg(nmax+1,t_COL);
   famod[1] = (long)ap;
-  if (nmax != FpX_split_berlekamp((GEN*)(famod+1), prime))
+  if (nmax != FpX_split_Berlekamp((GEN*)(famod+1), prime))
     err(bugparier,"DDF: wrong numbers of factors");
   if (DEBUGLEVEL>2)
   {
Index: src/basemath/polarit3.c
===================================================================
RCS file: /home/cvs/pari/src/basemath/polarit3.c,v
retrieving revision 1.177
diff -u -r1.177 polarit3.c
--- src/basemath/polarit3.c	19 Dec 2003 11:05:36 -0000	1.177
+++ src/basemath/polarit3.c	26 Dec 2003 21:30:07 -0000
@@ -25,14 +25,6 @@
 #define pswap(x,y) {GEN *_z=x; x=y; y=_z;}
 #define swap(x,y)  {GEN  _z=x; x=y; y=_z;}
 #define swapspec(x,y, nx,ny) {swap(x,y); lswap(nx,ny);}
-/* 2p^2 < 2^BITS_IN_LONG */
-#ifdef LONG_IS_64BIT
-#  define u_OK_ULONG(p) ((ulong)p <= 3037000493UL)
-#else
-#  define u_OK_ULONG(p) ((ulong)p <= 46337UL)
-#endif
-#define OK_ULONG(p) (lgefint(p) == 3 && u_OK_ULONG(p[2]))
-
 #define both_odd(x,y) ((x)&(y)&1)
 extern GEN caractducos(GEN p, GEN x, int v);
 extern double mylog2(GEN z);
@@ -77,347 +69,6 @@
 #  define MUL_LIMIT 10
 #endif
 
-/* multiplication in Fp[X], p small */
-
-GEN
-u_normalizepol(GEN x, long lx)
-{
-  long i;
-  for (i = lx-1; i>1; i--)
-    if (x[i]) break;
-  stackdummy(x + (i+1), lg(x) - (i+1));
-  setlg(x, i+1);
-  setsigne(x, (i > 1)); return x;
-}
-
-GEN
-u_FpX_deriv(GEN z, ulong p)
-{
-  long i,l = lg(z)-1;
-  GEN x;
-  if (l < 2) l = 2;
-  x = cgetg(l,t_VECSMALL); x[1] = z[1]; z++;
-  if (HIGHWORD(l | p))
-    for (i=2; i<l; i++) x[i] = (long)muluumod((ulong)i-1, (ulong)z[i], p);
-  else
-    for (i=2; i<l; i++) x[i] = ((i-1) * z[i]) % p;
-  return u_normalizepol(x,l);
-}
-
-static GEN
-u_FpX_gcd_i(GEN a, GEN b, ulong p)
-{
-  GEN c;
-  if (lg(b) > lg(a)) swap(a, b);
-  while (signe(b))
-  {
-    c = u_FpX_rem(a,b,p);
-    a = b; b = c;
-  }
-  return a;
-}
-
-GEN
-u_FpX_gcd(GEN a, GEN b, ulong p)
-{
-  pari_sp av = avma;
-  return gerepileupto(av, u_FpX_gcd_i(a,b,p));
-}
-
-int
-u_FpX_is_squarefree(GEN z, ulong p)
-{
-  pari_sp av = avma;
-  GEN d = u_FpX_gcd_i(z, u_FpX_deriv(z,p) , p);
-  avma = av; return (degpol(d) == 0);
-}
-
-static GEN
-u_FpX_addspec(GEN x, GEN y, ulong p, long lx, long ly)
-{
-  long i,lz;
-  GEN z;
-
-  if (ly>lx) swapspec(x,y, lx,ly);
-  lz = lx+2; z = cgetg(lz,t_VECSMALL) + 2;
-  for (i=0; i<ly; i++) z[i] = (long)adduumod((ulong)x[i], (ulong)y[i], p);
-  for (   ; i<lx; i++) z[i] = x[i];
-  z -= 2; z[1]=0; return u_normalizepol(z, lz);
-}
-
-INLINE ulong
-u_FpX_mullimb(GEN x, GEN y, ulong p, long a, long b)
-{
-  ulong p1 = 0;
-  long i;
-  for (i=a; i<b; i++)
-    if (y[i])
-    {
-      p1 += y[i] * x[-i];
-      if (p1 & HIGHBIT) p1 %= p;
-    }
-  return p1 % p;
-}
-
-ulong
-u_FpX_mullimb_safe(GEN x, GEN y, ulong p, long a, long b)
-{
-  ulong p1 = 0;
-  long i;
-  for (i=a; i<b; i++)
-    if (y[i])
-      p1 = adduumod(p1, muluumod(y[i],x[-i],p), p);
-  return p1;
-}
-
-/* assume nx >= ny > 0 */
-static GEN
-s_FpX_mulspec(GEN x, GEN y, ulong p, long nx, long ny)
-{
-  long i,lz,nz;
-  GEN z;
-
-  lz = nx+ny+1; nz = lz-2;
-  z = cgetg(lz, t_VECSMALL) + 2; /* x:y:z [i] = term of degree i */
-  if (u_OK_ULONG(p))
-  {
-    for (i=0; i<ny; i++)z[i] = (long)u_FpX_mullimb(x+i,y,p,0,i+1);
-    for (  ; i<nx; i++) z[i] = (long)u_FpX_mullimb(x+i,y,p,0,ny);
-    for (  ; i<nz; i++) z[i] = (long)u_FpX_mullimb(x+i,y,p,i-nx+1,ny);
-  }
-  else
-  {
-    for (i=0; i<ny; i++)z[i] = (long)u_FpX_mullimb_safe(x+i,y,p,0,i+1);
-    for (  ; i<nx; i++) z[i] = (long)u_FpX_mullimb_safe(x+i,y,p,0,ny);
-    for (  ; i<nz; i++) z[i] = (long)u_FpX_mullimb_safe(x+i,y,p,i-nx+1,ny);
-  }
-  z -= 2; z[1]=0; return u_normalizepol(z, lz);
-}
-
-/* return (x * X^d) + y. Assume d > 0, x > 0 and y >= 0 */
-GEN
-u_FpX_addshift(GEN x, GEN y, ulong p, long d)
-{
-  GEN xd,yd,zd = (GEN)avma;
-  long a,lz,ny = lg(y)-2, nx = lg(x)-2;
-
-  x += 2; y += 2; a = ny-d;
-  if (a <= 0)
-  {
-    lz = (a>nx)? ny+2: nx+d+2;
-    (void)new_chunk(lz); xd = x+nx; yd = y+ny;
-    while (xd > x) *--zd = *--xd;
-    x = zd + a;
-    while (zd > x) *--zd = 0;
-  }
-  else
-  {
-    xd = new_chunk(d); yd = y+d;
-    x = u_FpX_addspec(x,yd,p, nx,a);
-    lz = (a>nx)? ny+2: lg(x)+d;
-    x += 2; while (xd > x) *--zd = *--xd;
-  }
-  while (yd > y) *--zd = *--yd;
-  *--zd = evalsigne(1);
-  *--zd = evaltyp(t_VECSMALL) | evallg(lz); return zd;
-}
-
-#define u_FpX_mul(x,y, p) u_FpX_Kmul(x+2,y+2,p, lg(x)-2,lg(y)-2)
-#define u_FpX_sqr(x, p) u_FpX_Ksqr(x+2,p, lg(x)-2)
-#define u_FpX_add(x,y, p) u_FpX_addspec(x+2,y+2,p, lg(x)-2,lg(y)-2)
-
-GEN
-u_zeropol()
-{
-  GEN z = cgetg(2, t_VECSMALL);
-  z[1] = evalvarn(0); return z;
-}
-
-static GEN
-u_mallocpol(long d)
-{
-  GEN z = (GEN)gpmalloc((d+3) * sizeof(long));
-  z[0] = evaltyp(t_VECSMALL) | evallg(d+3);
-  z[1] = evalsigne((d >= 0)) | evalvarn(0);
-  return z;
-}
-GEN
-u_getpol(long d)
-{
-  GEN z = cgetg(d + 3, t_VECSMALL);
-  z[1] = evalsigne((d >= 0)) | evalvarn(0);
-  return z;
-}
-
-static GEN
-u_scalarpol(ulong x)
-{
-  GEN z;
-  if (!x) return u_zeropol();
-  z = u_getpol(0);
-  z[2] = (long)x; return z;
-}
-
-static GEN
-u_FpX_neg(GEN x, ulong p)
-{
-  long i, l = lg(x);
-  GEN z = cgetg(l, t_VECSMALL);
-  z[1] = x[1];
-  for (i=2; i<l; i++) z[i] = x[i]? (long)p - x[i]: 0;
-  return z;
-}
-
-static GEN
-u_FpX_neg_i(GEN x, ulong p)
-{
-  long i, l = lg(x);
-  for (i=2; i<l; i++)
-    if (x[i]) x[i] = (long)p - x[i];
-  return x;
-}
-
-/* shift polynomial + gerepile */
-static GEN
-u_FpX_shiftip(pari_sp av, GEN x, long v)
-{
-  long i, lx = lg(x), ly;
-  GEN y;
-  if (v <= 0 || !signe(x)) return gerepileupto(av, x);
-  avma = av; ly = lx + v;
-  x += lx; y = new_chunk(ly) + ly;
-  for (i = 2; i<lx; i++) *--y = *--x;
-  for (i = 0; i< v; i++) *--y = 0;
-  *--y = evalsigne(1);
-  *--y = evaltyp(t_VECSMALL) | evallg(ly); return y;
-}
-
-/* fast product (Karatsuba) of polynomials a,b. These are not real GENs, a+2,
- * b+2 were sent instead. na, nb = number of terms of a, b.
- * Only c, c0, c1, c2 are genuine GEN.
- */
-GEN
-u_FpX_Kmul(GEN a, GEN b, ulong p, long na, long nb)
-{
-  GEN a0,c,c0;
-  long n0, n0a, i, v = 0;
-  pari_sp av;
-
-  while (na && !a[0]) { a++; na--; v++; }
-  while (nb && !b[0]) { b++; nb--; v++; }
-  if (na < nb) swapspec(a,b, na,nb);
-  if (!nb) return u_zeropol();
-
-  av = avma;
-  if (nb < u_MUL_LIMIT)
-    return u_FpX_shiftip(av,s_FpX_mulspec(a,b,p,na,nb), v);
-  i=(na>>1); n0=na-i; na=i;
-  a0=a+n0; n0a=n0;
-  while (n0a && !a[n0a-1]) n0a--;
-
-  if (nb > n0)
-  {
-    GEN b0,c1,c2;
-    long n0b;
-
-    nb -= n0; b0 = b+n0; n0b = n0;
-    while (n0b && !b[n0b-1]) n0b--;
-    c =  u_FpX_Kmul(a,b,p,n0a,n0b);
-    c0 = u_FpX_Kmul(a0,b0,p,na,nb);
-
-    c2 = u_FpX_addspec(a0,a,p,na,n0a);
-    c1 = u_FpX_addspec(b0,b,p,nb,n0b);
-
-    c1 = u_FpX_mul(c1,c2,p);
-    c2 = u_FpX_add(c0,c,p);
-
-    c2 = u_FpX_neg_i(c2,p);
-    c2 = u_FpX_add(c1,c2,p);
-    c0 = u_FpX_addshift(c0,c2 ,p, n0);
-  }
-  else
-  {
-    c  = u_FpX_Kmul(a,b,p,n0a,nb);
-    c0 = u_FpX_Kmul(a0,b,p,na,nb);
-  }
-  c0 = u_FpX_addshift(c0,c,p,n0);
-  return u_FpX_shiftip(av,c0, v);
-}
-
-GEN
-u_FpX_sqrpol(GEN x, ulong p, long nx)
-{
-  long i,j,l,lz,nz;
-  ulong p1;
-  GEN z;
-
-  if (!nx) return u_zeropol();
-  lz = (nx << 1) + 1, nz = lz-2;
-  z = cgetg(lz,t_VECSMALL) + 2;
-  for (i=0; i<nx; i++)
-  {
-    p1 = 0; l = (i+1)>>1;
-    for (j=0; j<l; j++)
-    {
-      p1 += x[j] * x[i-j];
-      if (p1 & HIGHBIT) p1 %= p;
-    }
-    p1 <<= 1; if (p1 & HIGHBIT) p1 %= p;
-    if ((i&1) == 0)
-      p1 += x[i>>1] * x[i>>1];
-    z[i] = (long) (p1 % p);
-  }
-  for (  ; i<nz; i++)
-  {
-    p1 = 0; l = (i+1)>>1;
-    for (j=i-nx+1; j<l; j++)
-    {
-      p1 += x[j] * x[i-j];
-      if (p1 & HIGHBIT) p1 %= p;
-    }
-    p1 <<= 1; if (p1 & HIGHBIT) p1 %= p;
-    if ((i&1) == 0)
-      p1 += x[i>>1] * x[i>>1];
-    z[i] = (long) (p1 % p);
-  }
-  z -= 2; z[1]=0; return u_normalizepol(z,lz);
-}
-
-static GEN
-u_Fp_gmul2_1(GEN x, ulong p)
-{
-  long i, l = lg(x);
-  GEN z = cgetg(l, t_VECSMALL);
-  for (i=2; i<l; i++) z[i] = (long)adduumod(x[i], x[i], p);
-  z[1] = x[1]; return z;
-}
-
-GEN
-u_FpX_Ksqr(GEN a, ulong p, long na)
-{
-  GEN a0,c,c0,c1;
-  long n0, n0a, i, v = 0;
-  pari_sp av;
-
-  while (na && !a[0]) { a++; na--; v += 2; }
-  av = avma;
-  if (na < u_SQR_LIMIT) return u_FpX_shiftip(av, u_FpX_sqrpol(a,p,na), v);
-  i=(na>>1); n0=na-i; na=i;
-  a0=a+n0; n0a=n0;
-  while (n0a && !a[n0a-1]) n0a--;
-
-  c = u_FpX_Ksqr(a,p,n0a);
-  c0= u_FpX_Ksqr(a0,p,na);
-  if (p == 2) n0 *= 2;
-  else
-  {
-    c1 = u_Fp_gmul2_1(u_FpX_Kmul(a0,a,p, na,n0a), p);
-    c0 = u_FpX_addshift(c0,c1,p,n0);
-  }
-  c0 = u_FpX_addshift(c0,c,p,n0);
-  return u_FpX_shiftip(av,c0,v);
-}
-
 /* generic multiplication */
 
 static GEN
@@ -676,8 +327,8 @@
  *****************************************/
 
 /*********************************************************************
-This functions supposes polynomials already reduced.
-There are clean and memory efficient.
+These functions suppose polynomials to be already reduced.
+They are clean and memory efficient.
 **********************************************************************/
 
 GEN
@@ -778,22 +429,6 @@
   if (!p) return z;
   return FpX_red(z, p);
 }
-#define u_FpX_div(x,y,p) u_FpX_divrem((x),(y),(p),NULL)
-
-/* Product of y and x in Z/pZ[X]/(pol), as t_VECSMALL. Assume OK_ULONG(p) */
-static GEN
-u_FpXQ_mul(GEN y,GEN x,GEN pol,ulong p)
-{
-  GEN z = u_FpX_mul(y,x,p);
-  return u_FpX_rem(z,pol,p);
-}
-/* Square of y in Z/pZ[X]/(pol), as t_VECSMALL. Assume OK_ULONG(p) */
-static GEN
-u_FpXQ_sqr(GEN y,GEN pol,ulong p)
-{
-  GEN z = u_FpX_sqr(y,p);
-  return u_FpX_rem(z,pol,p);
-}
 
 /* Inverse of x in Z/pZ[X]/(pol) or NULL if inverse doesn't exist
  * return lift(1 / (x mod (p,pol))) */
@@ -802,10 +437,11 @@
 {
   GEN z, U, V;
 
-  if (typ(x) != t_POL) return mpinvmod(x, p);
+  if (typ(x) != t_POL) err(bugparier,"not a polynomial in FpXQ_invsafe");
   z = FpX_extgcd(x, T, p, &U, &V);
   if (degpol(z)) return NULL;
-  z = mpinvmod((GEN)z[2], p);
+  z = mpinvmodsafe((GEN)z[2], p);
+  if (!z) return NULL;
   return FpX_Fp_mul(U, z, p);
 }
 
@@ -932,7 +568,7 @@
   pari_sp av;
   GEN U;
 
-  if (!T) return mpinvmod(x,p);
+  if (!T) err(bugparier,"T=NULL in FpXQ_inv");
   av = avma;
   U = FpXQ_invsafe(x, T, p);
   if (!U) err(talker,"non invertible polynomial in FpXQ_inv");
@@ -1133,24 +769,8 @@
 typedef struct {
   GEN pol, p;
 } FpX_muldata;
-typedef struct {
-  GEN pol;
-  ulong p;
-} u_FpX_muldata;
 
 static GEN
-_u_sqr(void *data, GEN x)
-{
-  u_FpX_muldata *D = (u_FpX_muldata*)data;
-  return u_FpXQ_sqr(x, D->pol, D->p);
-}
-static GEN
-_u_mul(void *data, GEN x, GEN y)
-{
-  u_FpX_muldata *D = (u_FpX_muldata*)data;
-  return u_FpXQ_mul(x,y, D->pol, D->p);
-}
-static GEN
 _sqr(void *data, GEN x)
 {
   FpX_muldata *D = (FpX_muldata*)data;
@@ -1163,19 +783,6 @@
   return FpXQ_mul(x,y, D->pol, D->p);
 }
 
-/* assume n > 0 */
-GEN
-u_FpXQ_pow(GEN x, GEN n, GEN pol, ulong p)
-{
-  pari_sp av = avma;
-  u_FpX_muldata D;
-  GEN y;
-  D.pol = pol;
-  D.p   = p;
-  y = leftright_pow(x, n, (void*)&D, &_u_sqr, &_u_mul);
-  return gerepileupto(av, y);
-}
-
 /* x,pol in Z[X], p in Z, n in Z, compute lift(x^n mod (p, pol)) */
 GEN
 FpXQ_pow(GEN x, GEN n, GEN pol, GEN p)
@@ -1184,7 +791,7 @@
   FpX_muldata D;
   long vx;
   GEN y;
-  if (!pol) return powmodulo(x, n, p);
+  if (!pol) err(bugparier,"FpXQ_pow, pol is NULL");
   vx = varn(x);
   if (!signe(n)) return polun[vx];
   av = avma;
@@ -1195,13 +802,13 @@
   }
   else
     if (is_pm1(n)) return gcopy(x); /* n = 1 */
-  if (OK_ULONG(p))
+  if (!is_bigint(p))
   {
     ulong pp = p[2];
-    pol = u_Fp_FpX(pol, pp);
-    x   = u_Fp_FpX(x,   pp);
-    y = u_FpXQ_pow(x, n, pol, pp);
-    y = small_to_pol(y,vx);
+    pol = ZX_Flx(pol, pp);
+    x   = ZX_Flx(x,   pp);
+    y = Flxq_pow(x, n, pol, pp);
+    y = Flx_ZX(y);
   }
   else
   {
@@ -1213,16 +820,9 @@
   return gerepileupto(av, y);
 }
 
-GEN
-Fq_pow(GEN x, GEN n, GEN pol, GEN p)
-{
-  if (typ(x) == t_INT) return powmodulo(x,n,p);
-  return FpXQ_pow(x,n,pol,p);
-}
-
 /*******************************************************************/
 /*                                                                 */
-/*                             Fp[X][Y]                            */
+/*                             FpXX                                */
 /*                                                                 */
 /*******************************************************************/
 /*Polynomials whose coefficients are either polynomials or integers*/
@@ -1248,6 +848,26 @@
   return normalizepol_i(res,lg(res));
 }
 
+GEN 
+ZXX_Flxy(GEN B, ulong p, long vs)
+{
+  long lb=lg(B);
+  long i;
+  GEN b=cgetg(lb,t_POL);
+  b[1]=B[1];
+  for (i=2; i<lb; i++) 
+    switch (typ(B[i]))
+    {
+    case t_INT:  
+      b[i] = (long)Z_Flx((GEN)B[i], p, vs);
+      break;           
+    case t_POL:
+      b[i] = (long)ZX_Flx((GEN)B[i], p);
+      break;
+    }
+  return b;
+}
+
 /*******************************************************************/
 /*                                                                 */
 /*                             (Fp[X]/(Q))[Y]                      */
@@ -1277,21 +897,6 @@
 }
 
 GEN
-FpVQX_red(GEN z, GEN T, GEN p)
-{
-  GEN res;
-  int i, l = lg(z);
-  res=cgetg(l,typ(z));
-  for(i=1;i<l;i++)
-    if (typ(z[i]) == t_INT)
-      res[i] = lmodii((GEN)z[i],p);
-    else if (T)
-      res[i] = (long)FpX_res((GEN)z[i],T,p);
-    else
-      res[i] = (long)FpX_red((GEN)z[i],p);
-  return res;
-}
-GEN
 FpXQX_red(GEN z, GEN T, GEN p)
 {
   GEN res;
@@ -1378,7 +983,7 @@
   if (dg < 0) { swap(P, Q); dg = -dg; }
   for(;;)
   {
-    U = FpXQ_invsafe(leading_term(Q), T, p);
+    U = Fq_invsafe(leading_term(Q), T, p);
     if (!U) { avma = ltop; return NULL; }
     do /* set P := P % Q */
     {
@@ -1497,18 +1102,27 @@
 }
 /*******************************************************************/
 /*                                                                 */
-/*                             Fq[X]                               */
+/*                             Fq                                  */
 /*                                                                 */
 /*******************************************************************/
 
-/*Well nothing, for now. So we reuse FpXQX*/
-#define FqX_mul FpXQX_mul
-#define FqX_sqr FpXQX_sqr
-#define FqX_red FpXQX_red
-static GEN
-Fq_neg(GEN x, GEN T/*unused*/, GEN p)
+GEN
+Fq_add(GEN x, GEN y, GEN T/*unused*/, GEN p)
 {
   (void)T;
+  switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
+  {
+    case 0: return modii(addii(x,y),p);
+    case 1: return FpX_Fp_add(x,y,p);
+    case 2: return FpX_Fp_add(y,x,p);
+    case 3: return FpX_add(x,y,p);
+  }
+  return NULL;
+}
+
+GEN
+Fq_neg(GEN x, GEN T, GEN p)
+{
   switch(typ(x)==t_POL)
   {
     case 0: return signe(x)?subii(p,x):gzero;
@@ -1516,6 +1130,76 @@
   }
   return NULL;
 }
+
+GEN
+Fq_mul(GEN x, GEN y, GEN T, GEN p)
+{
+  switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
+  {
+    case 0: return modii(mulii(x,y),p);
+    case 1: return FpX_Fp_mul(x,y,p);
+    case 2: return FpX_Fp_mul(y,x,p);
+    case 3: return FpXQ_mul(x,y,T,p);
+  }
+  return NULL;
+}
+
+GEN
+Fq_neg_inv(GEN x, GEN T, GEN p)
+{
+  switch(typ(x)==t_POL)
+  {
+    case 0: return mpinvmod(negi(x),p);
+    case 1: return FpXQ_inv(FpX_neg(x,p),T,p);
+  }
+  return NULL;
+}
+
+GEN
+Fq_invsafe(GEN x, GEN pol, GEN p)
+{
+  if (typ(x) == t_INT) return mpinvmodsafe(x,p);
+  return FpXQ_invsafe(x,pol,p);
+}
+
+GEN
+Fq_inv(GEN x, GEN pol, GEN p)
+{
+  if (typ(x) == t_INT) return mpinvmod(x,p);
+  return FpXQ_inv(x,pol,p);
+}
+
+GEN
+Fq_pow(GEN x, GEN n, GEN pol, GEN p)
+{
+  if (typ(x) == t_INT) return powmodulo(x,n,p);
+  return FpXQ_pow(x,n,pol,p);
+}
+
+GEN
+Fq_res(GEN x, GEN T, GEN p)
+{
+  pari_sp ltop=avma;
+  switch(typ(x)==t_POL)
+  {
+    case 0: return modii(x,p);
+    case 1: return gerepileupto(ltop,FpX_res(FpX_red(x,p),T,p));
+  }
+  return NULL;
+}
+
+
+/*******************************************************************/
+/*                                                                 */
+/*                             Fq[X]                               */
+/*                                                                 */
+/*******************************************************************/
+
+/*Well nothing, for now. So we reuse FpXQX*/
+#define FqX_mul FpXQX_mul
+#define FqX_sqr FpXQX_sqr
+#define FqX_red FpXQX_red
+
 static GEN modulo,Tmodulo;
 static GEN fgmul(GEN a,GEN b){return FqX_mul(a,b,Tmodulo,modulo);}
 GEN 
@@ -1535,6 +1219,22 @@
   return gerepileupto(ltop, FqXV_mul(W, Tp, p));
 }
 
+GEN
+FqV_red(GEN z, GEN T, GEN p)
+{
+  GEN res;
+  int i, l = lg(z);
+  res=cgetg(l,typ(z));
+  for(i=1;i<l;i++)
+    if (typ(z[i]) == t_INT)
+      res[i] = lmodii((GEN)z[i],p);
+    else if (T)
+      res[i] = (long)FpX_res((GEN)z[i],T,p);
+    else
+      res[i] = (long)FpX_red((GEN)z[i],p);
+  return res;
+}
+
 /*******************************************************************/
 /*                                                                 */
 /*                       n-th ROOT in Fq                           */
@@ -1688,7 +1388,7 @@
 GEN
 Fp_inv_isom(GEN S,GEN T, GEN p)
 {
-  pari_sp lbot, ltop = avma;
+  pari_sp ltop = avma;
   GEN     M, V;
   int     n, i;
   long    x;
@@ -1700,9 +1400,8 @@
     V[i] = zero;
   V[2] = un;
   V = FpM_invimage(M,V,p);
-  lbot = avma;
-  V = gtopolyrev(V, x);
-  return gerepile(ltop, lbot, V);
+  V = vec_to_pol(V, x);
+  return gerepileupto(ltop, V);
 }
 
 /* Let M the matrix of the x^p Frobenius automorphism.
@@ -1777,7 +1476,7 @@
  * TODO: use left-right binary (tricky!)
  */
 GEN
-FpM_frobenius_pow(GEN M, long d, GEN T, GEN p)
+FpM_Frobenius_pow(GEN M, long d, GEN T, GEN p)
 {
   pari_sp ltop=avma;
   long i,l=degpol(T);
@@ -1992,7 +1691,7 @@
   n=degpol(P);
   m=n/d;
   if (DEBUGLEVEL>=4) (void)timer2();
-  M=FpM_frobenius_pow(MP,d,P,l);
+  M=FpM_Frobenius_pow(MP,d,P,l);
   if (DEBUGLEVEL>=4) msgtimer("Fp_factorgalois: frobenius power");
   Tl=gcopy(P); setvarn(Tl,w);
   V=cgetg(m+1,t_VEC);
@@ -2196,21 +1895,12 @@
 }
 
 GEN
-FpXQX_normalize(GEN z, GEN T, GEN p)
-{
-  GEN p1 = leading_term(z);
-  if (lg(z) == 2 || gcmp1(p1)) return z;
-  if (!T) return FpX_normalize(z,p);
-  return FpXQX_FpXQ_mul(z, FpXQ_inv(p1,T,p), T,p);
-}
-
-GEN
-small_to_vec(GEN z)
-{
-  long i, l = lg(z);
-  GEN x = cgetg(l,t_VEC);
-  for (i=1; i<l; i++) x[i] = lstoi(z[i]);
-  return x;
+FpXQX_normalize(GEN z, GEN T, GEN p)
+{
+  GEN p1 = leading_term(z);
+  if (lg(z) == 2 || gcmp1(p1)) return z;
+  if (!T) return FpX_normalize(z,p);
+  return FpXQX_FpXQ_mul(z, FpXQ_inv(p1,T,p), T,p);
 }
 
 GEN
@@ -2233,35 +1923,11 @@
 }
 
 GEN
-small_to_pol_i(GEN z, long l)
-{
-  long i;
-  GEN x = cgetg(l,t_POL);
-  for (i=2; i<l; i++) x[i] = lstoi(z[i]);
-  x[1] = z[1]; return x;
-}
-
-/* assume z[i] % p has been done */
-GEN
-small_to_pol(GEN z, long v)
-{
-  GEN x = small_to_pol_i(z, lg(z));
-  setvarn(x,v); return x;
-}
-/* same. In place */
-GEN
-small_to_pol_ip(GEN z, long v)
-{
-  long i, l = lg(z);
-  for (i=2; i<l; i++) z[i] = lstoi(z[i]);
-  settyp(z, t_POL); setvarn(z, v); return z;
-}
-
-GEN
 pol_to_small(GEN x)
 {
   long i, lx = lg(x);
-  GEN a = u_getpol(lx-3);
+  GEN a = cgetg(lx,t_VECSMALL);
+  a[1]=((ulong)x[1])&VARNBITS;
   for (i=2; i<lx; i++) a[i] = itos((GEN)x[i]);
   return a;
 }
@@ -2368,167 +2034,6 @@
 }
 #endif
 
-/* return x[0 .. dx] mod p as t_VECSMALL. Assume x a t_POL/VECSMALL/INT */
-GEN
-u_Fp_FpX(GEN x, ulong p)
-{
-  long i, lx;
-  GEN a;
-
-  switch (typ(x))
-  {
-    case t_VECSMALL: return x;
-    case t_INT: a = u_getpol(0);
-      a[2] = (long)umodiu(x, p); return a;
-  }
-  lx = lg(x); a = u_getpol(lx-3);
-  for (i=2; i<lx; i++) a[i] = (long)umodiu((GEN)x[i], p);
-  return u_normalizepol(a,lx);
-}
-
-GEN
-u_Fp_FpV(GEN x, ulong p)
-{
-  long i, n = lg(x);
-  GEN y = cgetg(n,t_VECSMALL);
-  for (i=1; i<n; i++) y[i] = (long)umodiu((GEN)x[i], p);
-  return y;
-}
-
-GEN
-u_Fp_FpM(GEN x, ulong p)
-{
-  long j,n = lg(x);
-  GEN y = cgetg(n,t_MAT);
-  if (n == 1) return y;
-  for (j=1; j<n; j++) y[j] = (long)u_Fp_FpV((GEN)x[j], p);
-  return y;
-}
-
-static GEN
-u_FpX_Fp_mul(GEN y, ulong x,ulong p)
-{
-  GEN z;
-  int i, l;
-  if (!x) return u_zeropol();
-  l = lg(y); z = u_getpol(l-3);
-  if (HIGHWORD(x | p))
-    for(i=2; i<l; i++) z[i] = (long)muluumod(y[i], x, p);
-  else
-    for(i=2; i<l; i++) z[i] = (y[i] * x) % p;
-  return z;
-}
-
-GEN
-u_FpX_normalize(GEN z, ulong p)
-{
-  long l = lg(z)-1;
-  ulong p1 = z[l]; /* leading term */
-  if (p1 == 1) return z;
-  return u_FpX_Fp_mul(z, invumod(p1,p), p);
-}
-
-static GEN
-u_copy(GEN x)
-{
-  long i, l = lg(x);
-  GEN z = u_getpol(l-3);
-  for (i=2; i<l; i++) z[i] = x[i];
-  return z;
-}
-
-/* as FpX_divres but working only on ulong types. ASSUME pr != ONLY_DIVIDES
- * if relevant, *pr is the last object on stack */
-GEN
-u_FpX_divrem(GEN x, GEN y, ulong p, GEN *pr)
-{
-  GEN z,q,c;
-  long dx,dy,dz,i,j;
-  ulong p1,inv;
-
-  if (pr == ONLY_REM) return u_FpX_rem(x, y, p);
-  dy = degpol(y);
-  if (!dy)
-  {
-    if (y[2] == 1UL)
-      q = u_copy(x);
-    else
-      q = u_FpX_Fp_mul(x, invumod(y[2], p), p);
-    if (pr) *pr = u_zeropol();
-    return q;
-  }
-  dx = degpol(x);
-  dz = dx-dy;
-  if (dz < 0)
-  {
-    q = u_zeropol();
-    if (pr) *pr = u_copy(x);
-    return q;
-  }
-  x += 2;
-  y += 2;
-  z = u_getpol(dz) + 2;
-  inv = (ulong)y[dy];
-  if (inv != 1UL) inv = invumod(inv,p);
-
-  if (u_OK_ULONG(p))
-  {
-    z[dz] = (long) (inv*x[dx]) % p;
-    for (i=dx-1; i>=dy; --i)
-    {
-      p1 = p - x[i]; /* compute -p1 instead of p1 (pb with ulongs otherwise) */
-      for (j=i-dy+1; j<=i && j<=dz; j++)
-      {
-        p1 += z[j]*y[i-j];
-        if (p1 & HIGHBIT) p1 %= p;
-      }
-      p1 %= p;
-      z[i-dy] = p1? (long) ((p - p1)*inv) % p: 0;
-    }
-  }
-  else
-  {
-    z[dz] = (long)muluumod(inv, x[dx], p);
-    for (i=dx-1; i>=dy; --i)
-    { /* compute -p1 instead of p1 (pb with ulongs otherwise) */
-      p1 = p - (ulong)x[i];
-      for (j=i-dy+1; j<=i && j<=dz; j++)
-        p1 = adduumod(p1, muluumod(z[j],y[i-j],p), p);
-      z[i-dy] = p1? (long)muluumod(p - p1, inv, p): 0;
-    }
-  }
-  q = u_normalizepol(z-2, dz+3);
-  if (!pr) return q;
-
-  c = u_getpol(dy) + 2;
-  if (u_OK_ULONG(p))
-  {
-    for (i=0; i<dy; i++)
-    {
-      p1 = (ulong)z[0]*y[i];
-      for (j=1; j<=i && j<=dz; j++)
-      {
-        p1 += (ulong)z[j]*y[i-j];
-        if (p1 & HIGHBIT) p1 %= p;
-      }
-      c[i] = (long)subuumod((ulong)x[i], p1%p, p);
-    }
-  }
-  else
-  {
-    for (i=0; i<dy; i++)
-    {
-      p1 = muluumod(z[0],y[i],p);
-      for (j=1; j<=i && j<=dz; j++)
-        p1 = adduumod(p1, muluumod(z[j],y[i-j],p), p);
-      c[i] = (long)subuumod((ulong)x[i], p1, p);
-    }
-  }
-  i=dy-1; while (i>=0 && !c[i]) i--;
-  c = u_normalizepol(c-2, i+3);
-  *pr = c; return q;
-}
-
 /*FIXME: Unify the following 3 divrem routines. Treat the case x,y (lifted) in
  * R[X], y non constant. Given: (lifted) [inv(), mul()], (delayed) red() in R */
 
@@ -2572,17 +2077,17 @@
   if (OK_ULONG(p))
   { /* assume ab != 0 mod p */
     ulong pp = (ulong)p[2];
-    GEN a = u_Fp_FpX(x, pp);
-    GEN b = u_Fp_FpX(y, pp);
-    z = u_FpX_divrem(a,b,pp, pr);
+    GEN a = ZX_Flx(x, pp);
+    GEN b = ZX_Flx(y, pp);
+    z = Flx_divres(a,b,pp, pr);
     avma = av0; /* HACK: assume pr last on stack, then z */
     setlg(z, lg(z)); z = dummycopy(z);
     if (pr && pr != ONLY_DIVIDES && pr != ONLY_REM)
     {
       setlg(*pr, lg(*pr)); *pr = dummycopy(*pr);
-      *pr = small_to_pol_ip(*pr,vx);
+      *pr = Flx_ZX_inplace(*pr);
     }
-    return small_to_pol_ip(z,vx);
+    return Flx_ZX_inplace(z);
   }
   lead = gcmp1(lead)? NULL: gclone(mpinvmod(lead,p));
   avma = av0;
@@ -2838,21 +2343,6 @@
   *pr = rem; return z-2;
 }
 
-static GEN
-FpX_gcd_long(GEN x, GEN y, GEN p)
-{
-  ulong pp = (ulong)p[2];
-  pari_sp av = avma;
-  GEN a,b;
-
-  (void)new_chunk((lg(x) + lg(y)) << 2); /* scratch space */
-  a = u_Fp_FpX(x, pp);
-  if (!signe(a)) { avma = av; return FpX_red(y,p); }
-  b = u_Fp_FpX(y, pp);
-  a = u_FpX_gcd_i(a,b, pp);
-  avma = av; return small_to_pol(a, varn(x));
-}
-
 /* x and y in Z[X], return lift(gcd(x mod p, y mod p)) */
 GEN
 FpX_gcd(GEN x, GEN y, GEN p)
@@ -2860,7 +2350,16 @@
   GEN a,b,c;
   pari_sp av0, av;
 
-  if (OK_ULONG(p)) return FpX_gcd_long(x,y,p);
+  if (OK_ULONG(p))
+  {
+    ulong pp=p[2];
+    av = avma;
+    (void)new_chunk((lg(x) + lg(y)) << 2); /* scratch space */
+    a = ZX_Flx(x, pp);
+    b = ZX_Flx(y, pp);
+    a = Flx_gcd_i(a,b, pp);
+    avma = av; return Flx_ZX(a);
+  }
   av0=avma;
   a = FpX_red(x, p); av = avma;
   b = FpX_red(y, p);
@@ -2892,30 +2391,9 @@
   avma = av; return gun;
 }
 
-GEN
-u_FpX_sub(GEN x, GEN y, ulong p)
-{
-  long i,lz,lx = lg(x), ly = lg(y);
-  GEN z;
-
-  if (ly <= lx)
-  {
-    lz = lx; z = cgetg(lz,t_VECSMALL);
-    for (i=2; i<ly; i++) z[i] = (long)subuumod((ulong)x[i],(ulong)y[i],p);
-    for (   ; i<lx; i++) z[i] = x[i];
-  }
-  else
-  {
-    lz = ly; z = cgetg(lz,t_VECSMALL);
-    for (i=2; i<lx; i++) z[i] = (long)subuumod((ulong)x[i],(ulong)y[i],p);
-    for (   ; i<ly; i++) z[i] = y[i]? (long)(p - y[i]): y[i];
-  }
-  z[1]=0; return u_normalizepol(z, lz);
-}
-
 /* list of u_FpX in gptr, return then as GEN */
 static void
-u_gerepilemany(pari_sp av, GEN* gptr[], long n, long v)
+u_gerepilemany(pari_sp av, GEN* gptr[], long n)
 {
   GEN *l = (GEN*)gpmalloc(n*sizeof(GEN));
   long i;
@@ -2926,45 +2404,25 @@
   /* copy them back, kill clones */
   for (--i; i>=0; i--)
   {
-    *(gptr[i]) = small_to_pol(l[i],v);
+    *(gptr[i]) = Flx_ZX(l[i]);
     gunclone(l[i]);
   }
   free(l);
 }
 
 static GEN
-u_FpX_extgcd(GEN a, GEN b, ulong p, GEN *ptu, GEN *ptv)
-{
-  GEN q,z,u,v, x = a, y = b;
-
-  u = u_zeropol();
-  v= u_scalarpol(1); /* v = 1 */
-  while (signe(y))
-  {
-    q = u_FpX_divrem(x,y,p,&z);
-    x = y; y = z; /* (x,y) = (y, x - q y) */
-    z = u_FpX_sub(u, u_FpX_mul(q,v, p), p);
-    u = v; v = z; /* (u,v) = (v, u - q v) */
-  }
-  z = u_FpX_sub(x, u_FpX_mul(b,u,p), p);
-  z = u_FpX_div(z,a,p);
-  *ptu = z;
-  *ptv = u; return x;
-}
-
-static GEN
 FpX_extgcd_long(GEN x, GEN y, GEN p, GEN *ptu, GEN *ptv)
 {
   ulong pp = (ulong)p[2];
   pari_sp av = avma;
   GEN a, b, d;
 
-  a = u_Fp_FpX(x, pp);
-  b = u_Fp_FpX(y, pp);
-  d = u_FpX_extgcd(a,b, pp, ptu,ptv);
+  a = ZX_Flx(x, pp);
+  b = ZX_Flx(y, pp);
+  d = Flx_extgcd(a,b, pp, ptu,ptv);
   {
     GEN *gptr[3]; gptr[0] = &d; gptr[1] = ptu; gptr[2] = ptv;
-    u_gerepilemany(av, gptr, 3, varn(x));
+    u_gerepilemany(av, gptr, 3);
   }
   return d;
 }
@@ -3067,7 +2525,7 @@
 
 /* return z = a mod q, b mod p (p,q) = 1. qinv = 1/q mod p */
 static GEN
-u_chrem_coprime(GEN a, ulong b, GEN q, ulong p, ulong qinv, GEN pq)
+u_chinese_coprime(GEN a, ulong b, GEN q, ulong p, ulong qinv, GEN pq)
 {
   ulong d, amod = umodiu(a, p);
   pari_sp av = avma;
@@ -3120,7 +2578,7 @@
   GEN h, lim = shifti(qp,-1);
   ulong qinv = invumod(umodiu(q,p), p);
   int stable = 1;
-  h = u_chrem_coprime(*H,Hp,q,p,qinv,qp);
+  h = u_chinese_coprime(*H,Hp,q,p,qinv,qp);
   if (h)
   {
     if (cmpii(h,lim) > 0) h = subii(h,qp);
@@ -3152,7 +2610,7 @@
   }
   for (i=2; i<lp; i++)
   {
-    h = u_chrem_coprime((GEN)H[i],Hp[i],q,p,qinv,qp);
+    h = u_chinese_coprime((GEN)H[i],Hp[i],q,p,qinv,qp);
     if (h)
     {
       if (cmpii(h,lim) > 0) h = subii(h,qp);
@@ -3161,7 +2619,7 @@
   }
   for (   ; i<l; i++)
   {
-    h = u_chrem_coprime((GEN)H[i],   0,q,p,qinv,qp);
+    h = u_chinese_coprime((GEN)H[i],   0,q,p,qinv,qp);
     if (h)
     {
       if (cmpii(h,lim) > 0) h = subii(h,qp);
@@ -3181,7 +2639,7 @@
   for (j=1; j<l; j++)
     for (i=1; i<m; i++)
     {
-      h = u_chrem_coprime(gcoeff(H,i,j), coeff(Hp,i,j),q,p,qinv,qp);
+      h = u_chinese_coprime(gcoeff(H,i,j), coeff(Hp,i,j),q,p,qinv,qp);
       if (h)
       {
         if (cmpii(h,lim) > 0) h = subii(h,qp);
@@ -3213,113 +2671,12 @@
   y[0] = evaltyp(t_POL) | evallg(l); return y;
 }
 
-/* separate from u_FpX_divrem for maximal speed. Implicitly malloc = 0  */
-GEN
-u_FpX_rem(GEN x, GEN y, ulong p)
-{
-  GEN z, c;
-  long dx,dy,dz,i,j;
-  ulong p1,inv;
-
-  dy = degpol(y); if (!dy) return u_zeropol();
-  dx = degpol(x);
-  dz = dx-dy; if (dz < 0) return u_copy(x);
-  x += 2;
-  y += 2;
-  z = u_mallocpol(dz) + 2;
-  inv = y[dy];
-  if (inv != 1UL) inv = invumod(inv,p);
-
-  c = u_getpol(dy) + 2;
-  if (u_OK_ULONG(p))
-  {
-    z[dz] = (inv*x[dx]) % p;
-    for (i=dx-1; i>=dy; --i)
-    {
-      p1 = p - x[i]; /* compute -p1 instead of p1 (pb with ulongs otherwise) */
-      for (j=i-dy+1; j<=i && j<=dz; j++)
-      {
-        p1 += z[j]*y[i-j];
-        if (p1 & HIGHBIT) p1 %= p;
-      }
-      p1 %= p;
-      z[i-dy] = p1? ((p - p1)*inv) % p: 0;
-    }
-    for (i=0; i<dy; i++)
-    {
-      p1 = z[0]*y[i];
-      for (j=1; j<=i && j<=dz; j++)
-      {
-        p1 += z[j]*y[i-j];
-        if (p1 & HIGHBIT) p1 %= p;
-      }
-      c[i] = (long)subuumod((ulong)x[i], p1%p, p);
-    }
-  }
-  else
-  {
-    z[dz] = muluumod(inv, x[dx], p);
-    for (i=dx-1; i>=dy; --i)
-    {
-      p1 = p - x[i]; /* compute -p1 instead of p1 (pb with ulongs otherwise) */
-      for (j=i-dy+1; j<=i && j<=dz; j++)
-        p1 = adduumod(p1, muluumod(z[j],y[i-j],p), p);
-      z[i-dy] = p1? muluumod(p - p1, inv, p): 0;
-    }
-    for (i=0; i<dy; i++)
-    {
-      p1 = muluumod(z[0],y[i],p);
-      for (j=1; j<=i && j<=dz; j++)
-        p1 = adduumod(p1, muluumod(z[j],y[i-j],p), p);
-      c[i] = (long)subuumod((ulong)x[i], p1, p);
-    }
-  }
-  i = dy-1; while (i>=0 && !c[i]) i--;
-  free(z-2); return u_normalizepol(c-2, i+3);
-}
-
-ulong
-u_FpX_resultant(GEN a, GEN b, ulong p)
-{
-  long da,db,dc,cnt;
-  ulong lb, res = 1UL;
-  pari_sp av;
-  GEN c;
-
-  if (!signe(a) || !signe(b)) return 0;
-  da = degpol(a);
-  db = degpol(b);
-  if (db > da)
-  {
-    swapspec(a,b, da,db);
-    if (both_odd(da,db)) res = p-res;
-  }
-  if (!da) return 1; /* = res * a[2] ^ db, since 0 <= db <= da = 0 */
-  cnt = 0; av = avma;
-  while (db)
-  {
-    lb = b[db+2];
-    c = u_FpX_rem(a,b, p);
-    a = b; b = c; dc = degpol(c);
-    if (dc < 0) { avma = av; return 0; }
-
-    if (both_odd(da,db)) res = p - res;
-    if (lb != 1) res = muluumod(res, powuumod(lb, da - dc, p), p);
-    if (++cnt == 4) { cnt = 0; avma = av; }
-    da = db; /* = degpol(a) */
-    db = dc; /* = degpol(b) */
-  }
-  avma = av; return muluumod(res, powuumod(b[2], da, p), p);
-}
-
 static GEN
 muliimod(GEN x, GEN y, GEN p)
 {
   return modii(mulii(x,y), p);
 }
 
-#define FpX_rem(x,y,p) FpX_divres((x),(y),(p),ONLY_REM)
-
 /* Res(A,B) = Res(B,R) * lc(B)^(a-r) * (-1)^(ab), with R=A%B, a=deg(A) ...*/
 GEN
 FpX_resultant(GEN a, GEN b, GEN p)
@@ -3341,7 +2698,7 @@
   while (db)
   {
     lb = (GEN)b[db+2];
-    c = FpX_rem(a,b, p);
+    c = FpX_res(a,b, p);
     a = b; b = c; dc = degpol(c);
     if (dc < 0) { avma = av; return 0; }
 
@@ -3359,50 +2716,6 @@
   return gerepileuptoint(av, res);
 }
 
-/* If resultant is 0, *ptU and *ptU are not set */
-ulong
-u_FpX_extresultant(GEN a, GEN b, ulong p, GEN *ptU, GEN *ptV)
-{
-  GEN z,q,u,v, x = a, y = b;
-  ulong lb, res = 1UL;
-  pari_sp av = avma;
-  long dx,dy,dz;
-
-  if (!signe(x) || !signe(y)) return 0;
-  dx = degpol(x);
-  dy = degpol(y);
-  if (dy > dx)
-  {
-    swap(x,y); lswap(dx,dy); pswap(ptU, ptV);
-    a = x; b = y;
-    if (both_odd(dx,dy)) res = p-res;
-  }
-  u = u_zeropol();
-  v = u_scalarpol(1); /* v = 1 */
-  while (dy)
-  { /* b u = x (a), b v = y (a) */
-    lb = y[dy+2];
-    q = u_FpX_divrem(x,y, p, &z);
-    x = y; y = z; /* (x,y) = (y, x - q y) */
-    dz = degpol(z); if (dz < 0) { avma = av; return 0; }
-    z = u_FpX_sub(u, u_FpX_mul(q,v, p), p);
-    u = v; v = z; /* (u,v) = (v, u - q v) */
-
-    if (both_odd(dx,dy)) res = p - res;
-    if (lb != 1) res = muluumod(res, powuumod(lb, dx-dz, p), p);
-    dx = dy; /* = degpol(x) */
-    dy = dz; /* = degpol(y) */
-  }
-  res = muluumod(res, powuumod(y[2], dx, p), p);
-  lb = muluumod(res, invumod(y[2],p), p);
-  v = gerepileupto(av, u_FpX_Fp_mul(v, lb, p));
-  av = avma;
-  u = u_FpX_sub(u_scalarpol(res), u_FpX_mul(b,v,p), p);
-  u = gerepileupto(av, u_FpX_div(u,a,p)); /* = (res - b v) / a */
-  *ptU = u;
-  *ptV = v; return res;
-}
-
 /* assuming the PRS finishes on a degree 1 polynomial C0 + C1X, with "generic"
  * degree sequence as given by dglist, set *Ci and return resultant(a,b) */
 ulong
@@ -3414,7 +2727,7 @@
   GEN c;
 
   if (C0) { *C0 = 1; *C1 = 0; }
-  if (!signe(a) || !signe(b)) return 0;
+  if (lg(a)-2==0 || lg(b)-2==0) return 0;
   da = degpol(a);
   db = degpol(b);
   if (db > da)
@@ -3428,7 +2741,7 @@
   while (db)
   {
     lb = b[db+2];
-    c = u_FpX_rem(a,b, p);
+    c = Flx_res(a,b, p);
     a = b; b = c; dc = degpol(c);
     if (dc < 0) { avma = av; return 0; }
 
@@ -3468,41 +2781,6 @@
   avma = av; return muluumod(res, lb, p);
 }
 
-static ulong global_pp;
-static GEN
-_u_FpX_mul(GEN a, GEN b)
-{
-  return u_FpX_mul(a,b, global_pp);
-}
-
-/* compute prod (x - a[i]) */
-GEN
-u_FpV_roots_to_pol(GEN a, ulong p)
-{
-  long i,k,lx = lg(a);
-  GEN p1,p2;
-  if (lx == 1) return polun[0];
-  p1 = cgetg(lx, t_VEC); global_pp = p;
-  for (k=1,i=1; i<lx-1; i+=2)
-  {
-    p2 = cgetg(5,t_VECSMALL); p1[k++] = (long)p2;
-    p2[2] = muluumod(a[i], a[i+1], p);
-    p2[3] = a[i] + a[i+1];
-    if ((ulong)p2[3] >= p) p2[3] -= p;
-    if (p2[3]) p2[3] = p - p2[3]; /* - (a[i] + a[i+1]) mod p */
-    p2[4] = 1; p2[1] = 0;
-  }
-  if (i < lx)
-  {
-    p2 = cgetg(4,t_POL); p1[k++] = (long)p2;
-    p2[1] = 0;
-    p2[2] = p - a[i];
-    p2[3] = 1;
-  }
-  setlg(p1, k); return divide_conquer_prod(p1, _u_FpX_mul);
-}
-
-
 /* u P(X) + v P(-X) */
 static GEN
 pol_comp(GEN P, GEN u, GEN v)
@@ -3519,23 +2797,6 @@
   y[1] = P[1]; return normalizepol_i(y,l);
 }
 
-static GEN
-u_pol_comp(GEN P, ulong u, ulong v, ulong p)
-{
-  long i, l = lg(P);
-  GEN y = u_getpol(l-3);
-  for (i=2; i<l; i++)
-  {
-    ulong t = P[i];
-    y[i] = (t == 0)? 0:
-                     (i&1)? muluumod(t, u + (p - v), p)
-                          : muluumod(t, u + v, p);
-  }
-  return u_normalizepol(y,l);
-}
-
-extern GEN roots_to_pol(GEN a, long v);
-
 GEN
 polint_triv(GEN xa, GEN ya)
 {
@@ -3568,72 +2829,6 @@
   return P? P: zeropol(0);
 }
 
-ulong
-u_FpX_eval(GEN x, ulong y, ulong p)
-{
-  ulong p1,r;
-  long j, i=lg(x)-1;
-  if (i<=2)
-    return (i==2)? x[2]: 0;
-  p1 = x[i];
-  /* specific attention to sparse polynomials (see poleval)*/
-  if (u_OK_ULONG(p))
-  {
-    for (i--; i>=2; i=j-1)
-    {
-      for (j=i; !x[j]; j--)
-        if (j==2)
-        {
-          if (i != j) y = powuumod(y, i-j+1, p);
-          return (p1 * y) % p;
-        }
-      r = (i==j)? y: powuumod(y, i-j+1, p);
-      p1 = ((p1*r) + x[j]) % p;
-    }
-  }
-  else
-  {
-    for (i--; i>=2; i=j-1)
-    {
-      for (j=i; !x[j]; j--)
-        if (j==2)
-        {
-          if (i != j) y = powuumod(y, i-j+1, p);
-          return muluumod(p1, y, p);
-        }
-      r = (i==j)? y: powuumod(y, i-j+1, p);
-      p1 = adduumod((ulong)x[j], muluumod(p1,r,p), p);
-    }
-  }
-  return p1;
-}
-
-static GEN
-u_FpX_div_by_X_x(GEN a, ulong x, ulong p)
-{
-  long l = lg(a), i;
-  GEN z = u_getpol(l-4), a0, z0;
-  a0 = a + l-1;
-  z0 = z + l-2; *z0 = *a0--;
-  if (u_OK_ULONG(p))
-  {
-    for (i=l-3; i>1; i--) /* z[i] = (a[i+1] + x*z[i+1]) % p */
-    {
-      ulong t = (*a0-- + x *  *z0--) % p;
-      *z0 = t;
-    }
-  }
-  else
-  {
-    for (i=l-3; i>1; i--)
-    {
-      ulong t = adduumod((ulong)*a0--, muluumod(x, *z0--, p), p);
-      *z0 = (long)t;
-    }
-  }
-  return z;
-}
-
 static GEN
 FpX_div_by_X_x(GEN a, GEN x, GEN p)
 {
@@ -3650,31 +2845,6 @@
   return z;
 }
 
-/* xa, ya = t_VECSMALL */
-GEN
-u_FpV_polint(GEN xa, GEN ya, ulong p)
-{
-  GEN T,dP, P = NULL, Q = u_FpV_roots_to_pol(xa, p);
-  long i, n = lg(xa);
-  ulong inv;
-  pari_sp av = avma;
-  for (i=1; i<n; i++)
-  {
-    if (!ya[i]) continue;
-    T = u_FpX_div_by_X_x(Q, xa[i], p);
-    inv = invumod(u_FpX_eval(T,xa[i], p), p);
-    if (i < n-1 && (ulong)(xa[i] + xa[i+1]) == p)
-    {
-      dP = u_pol_comp(T, muluumod(ya[i],inv,p), muluumod(ya[i+1],inv,p), p);
-      i++; /* x_i = -x_{i+1} */
-    }
-    else
-      dP = u_FpX_Fp_mul(T, muluumod(ya[i],inv,p), p);
-    P = P? u_FpX_add(P, dP, p): dP;
-  }
-  return P? gerepileupto(av, P): u_zeropol();
-}
-
 GEN
 FpV_polint(GEN xa, GEN ya, GEN p)
 {
@@ -3708,7 +2878,7 @@
 static void
 u_FpV_polint_all(GEN xa, GEN ya, GEN C0, GEN C1, ulong p)
 {
-  GEN T,Q = u_FpV_roots_to_pol(xa, p);
+  GEN T,Q = Flv_roots_to_pol(xa, p, 0);
   GEN dP  = NULL,  P = NULL;
   GEN dP0 = NULL, P0= NULL;
   GEN dP1 = NULL, P1= NULL;
@@ -3716,28 +2886,28 @@
   ulong inv;
   for (i=1; i<n; i++)
   {
-    T = u_FpX_div_by_X_x(Q, xa[i], p);
-    inv = invumod(u_FpX_eval(T,xa[i], p), p);
+    T = Flx_div_by_X_x(Q, xa[i], p);
+    inv = invumod(Flx_eval(T,xa[i], p), p);
 
     if (ya[i])
     {
-      dP = u_FpX_Fp_mul(T, muluumod(ya[i],inv,p), p);
-      P = P ? u_FpX_add(P , dP , p): dP;
+      dP = Flx_Fl_mul(T, muluumod(ya[i],inv,p), p);
+      P = P ? Flx_add(P , dP , p): dP;
     }
     if (C0[i])
     {
-      dP0= u_FpX_Fp_mul(T, muluumod(C0[i],inv,p), p);
-      P0= P0? u_FpX_add(P0, dP0, p): dP0;
+      dP0= Flx_Fl_mul(T, muluumod(C0[i],inv,p), p);
+      P0= P0? Flx_add(P0, dP0, p): dP0;
     }
     if (C1[i])
     {
-      dP1= u_FpX_Fp_mul(T, muluumod(C1[i],inv,p), p);
-      P1= P1? u_FpX_add(P1, dP1, p): dP1;
+      dP1= Flx_Fl_mul(T, muluumod(C1[i],inv,p), p);
+      P1= P1? Flx_add(P1, dP1, p): dP1;
     }
   }
-  ya[1] = (long) (P ? P : u_zeropol());
-  C0[1] = (long) (P0? P0: u_zeropol());
-  C1[1] = (long) (P1? P1: u_zeropol());
+  ya[1] = (long) (P ? P : Flx_zero(0));
+  C0[1] = (long) (P0? P0: Flx_zero(0));
+  C1[1] = (long) (P1? P1: Flx_zero(0));
 }
 
 /* b a vector of polynomials representing B in Fp[X][Y], evaluate at X = x,
@@ -3747,11 +2917,12 @@
 {
   GEN z;
   long i, lb = lg(b);
-  ulong leadz = u_FpX_eval((GEN)b[lb-1], x, p);
-  if (!leadz) return u_zeropol();
+  ulong leadz = Flx_eval(leading_term(b), x, p);
+  long vs=mael(b,2,1);
+  if (!leadz) return Flx_zero(vs);
 
-  z = u_getpol(lb-3);
-  for (i=2; i<lb-1; i++) z[i] = u_FpX_eval((GEN)b[i], x, p);
+  z = cgetg(lb, t_VECSMALL); z[1] = vs;
+  for (i=2; i<lb-1; i++) z[i] = Flx_eval((GEN)b[i], x, p);
   z[i] = leadz; return z;
 }
 
@@ -3761,9 +2932,10 @@
 {
   GEN z;
   long i, lb = lg(b);
-  z = u_getpol(lb-3);
-  for (i=2; i<lb; i++) z[i] = u_FpX_eval((GEN)b[i], x, p);
-  z = u_normalizepol(z, lb);
+  z = cgetg(lb,t_VECSMALL); z[1]=mael(b,2,1);
+
+  for (i=2; i<lb; i++) z[i] = Flx_eval((GEN)b[i], x, p);
+  z = Flx_renormalize(z, lb);
   *drop = lb - lg(z);
   return z;
 }
@@ -3812,7 +2984,7 @@
 {
   int drop;
   GEN ev = u_vec_FpX_eval_gen(b, n, p, &drop);
-  ulong r = u_FpX_resultant(a, ev, p);
+  ulong r = Flx_resultant(a, ev, p);
   if (drop && la != 1) r = muluumod(r, powuumod(la, drop,p),p);
   return r;
 }
@@ -3842,7 +3014,7 @@
   {
     x[++i] = 0;   y[i] = u_FpX_eval_resultant(a,b, x[i], p,la);
   }
-  return u_FpV_polint(x,y, p);
+  return Flv_polint(x,y, p, b[1]);
 }
 
 /* x^n mod p */
@@ -3855,22 +3027,6 @@
     return powuumod(x, (ulong)n, p);
 }
 
-/* x^n mod p, assume n > 0 */
-static GEN
-u_FpX_pow(GEN x, long n, ulong p)
-{
-  GEN y = u_scalarpol(1), z;
-  long m;
-  m = n;
-  z = x;
-  for (;;)
-  {
-    if (m&1) y = u_FpX_mul(y,z, p);
-    m >>= 1; if (!m) return y;
-    z = u_FpX_mul(z,z, p);
-  }
-}
-
 static GEN
 u_FpXX_pseudorem(GEN x, GEN y, ulong p)
 {
@@ -3884,13 +3040,13 @@
   av2 = avma; lim = stack_lim(av2,1);
   for (;;)
   {
-    x[0] = (long)u_FpX_neg((GEN)x[0], p); dp--;
+    x[0] = (long)Flx_neg((GEN)x[0], p); dp--;
     for (i=1; i<=dy; i++)
-      x[i] = (long)u_FpX_add( u_FpX_mul((GEN)y[0], (GEN)x[i], p),
-                              u_FpX_mul((GEN)x[0], (GEN)y[i], p), p );
+      x[i] = (long)Flx_add( Flx_mul((GEN)y[0], (GEN)x[i], p),
+                              Flx_mul((GEN)x[0], (GEN)y[i], p), p );
     for (   ; i<=dx; i++)
-      x[i] = (long)u_FpX_mul((GEN)y[0], (GEN)x[i], p);
-    do { x++; dx--; } while (dx >= 0 && !signe((GEN)x[0]));
+      x[i] = (long)Flx_mul((GEN)y[0], (GEN)x[i], p);
+    do { x++; dx--; } while (dx >= 0 && lg((GEN)x[0])==2);
     if (dx < dy) break;
     if (low_stack(lim,stack_lim(av2,1)))
     {
@@ -3898,16 +3054,16 @@
       gerepilemanycoeffs(av2,x,dx+1);
     }
   }
-  if (dx < 0) return u_zeropol();
+  if (dx < 0) return Flx_zero(0);
   lx = dx+3; x -= 2;
   x[0]=evaltyp(t_POL) | evallg(lx);
   x[1]=evalsigne(1) | evalvarn(vx);
   x = revpol(x) - 2;
   if (dp)
   { /* multiply by y[0]^dp   [beware dummy vars from FpY_FpXY_resultant] */
-    GEN t = u_FpX_pow((GEN)y[0], dp, p);
+    GEN t = Flx_pow((GEN)y[0], dp, p);
     for (i=2; i<lx; i++)
-      x[i] = (long)u_FpX_mul((GEN)x[i], t, p);
+      x[i] = (long)Flx_mul((GEN)x[i], t, p);
   }
   return gerepilecopy(av, x);
 }
@@ -3923,12 +3079,12 @@
     if (t == 1) return x;
     t = invumod(t, p);
     l = lg(x); z = cgetg(l, t_POL); z[1] = x[1];
-    for (i=2; i<l; i++) z[i] = (long)u_FpX_Fp_mul((GEN)x[i],t,p);
+    for (i=2; i<l; i++) z[i] = (long)Flx_Fl_mul((GEN)x[i],t,p);
   }
   else
   {
     l = lg(x); z = cgetg(l, t_POL); z[1] = x[1];
-    for (i=2; i<l; i++) z[i] = (long)u_FpX_div((GEN)x[i],y,p);
+    for (i=2; i<l; i++) z[i] = (long)Flx_div((GEN)x[i],y,p);
   }
   return z;
 }
@@ -3947,9 +3103,9 @@
     if (both_odd(dx, dy)) signh = -signh;
   }
   if (dy < 0) return gzero;
-  if (dy==0) return gerepileupto(av, u_FpX_pow((GEN)v[2],dx,p));
+  if (dy==0) return gerepileupto(av, Flx_pow((GEN)v[2],dx,p));
 
-  g = h = u_scalarpol(1); av2 = avma; lim = stack_lim(av2,1);
+  g = h = Fl_Flx(1,0); av2 = avma; lim = stack_lim(av2,1);
   for(;;)
   {
     r = u_FpXX_pseudorem(u,v,p); dr = lg(r);
@@ -3960,10 +3116,10 @@
     {
       case 0: break;
       case 1:
-        p1 = u_FpX_mul(h,p1, p); h = g; break;
+        p1 = Flx_mul(h,p1, p); h = g; break;
       default:
-        p1 = u_FpX_mul(u_FpX_pow(h,degq,p), p1, p);
-        h = u_FpX_div(u_FpX_pow(g,degq,p), u_FpX_pow(h,degq-1,p), p);
+        p1 = Flx_mul(Flx_pow(h,degq,p), p1, p);
+        h = Flx_div(Flx_pow(g,degq,p), Flx_pow(h,degq-1,p), p);
     }
     if (both_odd(du,dv)) signh = -signh;
     v = u_FpX_divexact(r, p1, p);
@@ -3975,8 +3131,8 @@
     }
   }
   z = (GEN)v[2];
-  if (dv > 1) z = u_FpX_div(u_FpX_pow(z,dv,p), u_FpX_pow(h,dv-1,p), p);
-  if (signh < 0) z = u_FpX_neg(z,p);
+  if (dv > 1) z = Flx_div(Flx_pow(z,dv,p), Flx_pow(h,dv-1,p), p);
+  if (signh < 0) z = Flx_neg(z,p);
   return gerepileupto(av, z);
 }
 
@@ -4009,16 +3165,16 @@
     ulong pp = (ulong)p[2];
     long l = lg(b);
 
-    a = u_Fp_FpX(a, pp);
+    a = ZX_Flx(a, pp);
     for (i=2; i<l; i++)
-      b[i] = (long)u_Fp_FpX((GEN)b[i], pp);
+      b[i] = (long)ZX_Flx((GEN)b[i], pp);
     if ((ulong)dres >= pp)
     {
       l = lg(a);
       a[0] = evaltyp(t_POL) | evallg(l);
       a[1] = evalsigne(1) | evalvarn(vY);
       for (i=2; i<l; i++)
-        a[i] = (long)u_scalarpol(a[i]);
+        a[i] = (long)Fl_Flx(a[i],evalvarn(vX));
       x = u_FpYX_subres(a, b, pp);
     }
     else
@@ -4027,7 +3183,7 @@
       y = cgetg(dres+2, t_VECSMALL);
       x = u_FpY_FpXY_resultant(a, b, pp, (ulong)dres, x, y);
     }
-    return small_to_pol(x, vX);
+    return Flx_ZX(x);
   }
  
   la = leading_term(a);
@@ -4131,7 +3287,7 @@
   while (p < (ulong)(dres<<1)) NEXT_PRIME_VIADIFF(p,d);
 
   H = H0 = H1 = NULL;
-  lb = lg(B); b = u_getpol(degpol(B));
+  lb = lg(B); 
   bound = ZY_ZXY_ResBound(A, B, dB);
   if (DEBUGLEVEL>4) fprintferr("bound for resultant coeffs: 2^%ld\n",bound);
   check_theta(bound);
@@ -4142,8 +3298,8 @@
     NEXT_PRIME_VIADIFF_CHECK(p,d);
     if (dB) { dp = smodis(dB, p); if (!dp) continue; }
 
-    a = u_Fp_FpX(A, p);
-    for (i=2; i<lb; i++) b[i] = (long)u_Fp_FpX((GEN)B[i], p);
+    a = ZX_Flx(A, p);
+    b = ZXX_Flxy(B, p, A[1]);
     if (LERS)
     {
       if (!b[lb-1] || degpol(a) < degA) continue; /* p | lc(A)lc(B) */
@@ -4184,9 +3340,9 @@
     else
       Hp = u_FpY_FpXY_resultant(a, b, p, (ulong)dres, x, y);
     if (!H && degpol(Hp) != dres) continue;
-    if (dp != 1) Hp = u_FpX_Fp_mul(Hp, powusmod(dp,-degA,p), p);
+    if (dp != 1) Hp = Flx_Fl_mul(Hp, powusmod(dp,-degA,p), p);
     if (checksqfree) {
-      if (!u_FpX_is_squarefree(Hp, p)) goto INIT;
+      if (!Flx_is_squarefree(Hp, p)) goto INIT;
       if (DEBUGLEVEL>4) fprintferr("Final lambda = %ld\n",*lambda);
       checksqfree = 0;
     }
@@ -4219,7 +3375,7 @@
     {
       GEN *gptr[4]; gptr[0] = &H; gptr[1] = &q; gptr[2] = &H0; gptr[3] = &H1;
       if (DEBUGMEM>1) err(warnmem,"ZY_ZXY_resultant");
-      gerepilemany(av2,gptr,LERS? 4: 2); b = u_getpol(degpol(B));
+      gerepilemany(av2,gptr,LERS? 4: 2); //b = u_getpol(degpol(B));
     }
   }
 END:
@@ -4339,9 +3495,9 @@
     NEXT_PRIME_VIADIFF_CHECK(p,d);
     if (dB) { dp = smodis(dB, p); if (!dp) continue; }
 
-    a = u_Fp_FpX(A, p);
-    b = u_Fp_FpX(B, p);
-    Hp= u_FpX_resultant(a, b, p);
+    a = ZX_Flx(A, p);
+    b = ZX_Flx(B, p);
+    Hp= Flx_resultant(a, b, p);
     if (dp != 1) Hp = muluumod(Hp, powusmod(dp, -degA, p), p);
     if (!H)
     {
@@ -4442,18 +3598,18 @@
     NEXT_PRIME_VIADIFF_CHECK(p,d);
     if (!umodiu(g,p)) continue;
 
-    a = u_Fp_FpX(A, p);
-    b = u_Fp_FpX(B, p); Hp = u_FpX_gcd_i(a,b, p);
+    a = ZX_Flx(A, p);
+    b = ZX_Flx(B, p); Hp = Flx_gcd_i(a,b, p);
     m = degpol(Hp);
     if (m == 0) { H = polun[varn(A0)]; break; } /* coprime. DONE */
     if (m > n) continue; /* p | Res(A/G, B/G). Discard */
 
     if (is_pm1(g)) /* make sure lead(H) = g mod p */
-      Hp = u_FpX_normalize(Hp, p);
+      Hp = Flx_normalize(Hp, p);
     else
     {
       ulong t = muluumod(umodiu(g, p), invumod(Hp[m+2],p), p);
-      Hp = u_FpX_Fp_mul(Hp, t, p);
+      Hp = Flx_Fl_mul(Hp, t, p);
     }
     if (m < n)
     { /* First time or degree drop [all previous p were as above; restart]. */
@@ -4505,10 +3661,10 @@
   for(;;)
   {
     NEXT_PRIME_VIADIFF_CHECK(p,d);
-    a = u_Fp_FpX(A, p);
-    b = u_Fp_FpX(B, p);
+    a = ZX_Flx(A, p);
+    b = ZX_Flx(B, p);
     /* if p | Res(A/G, B/G), discard */
-    if (!u_FpX_extresultant(b,a,p, &Vp,&Up)) continue;
+    if (!Flx_extresultant(b,a,p, &Vp,&Up)) continue;
 
     if (!U)
     { /* First time */
Index: src/headers/paricom.h
===================================================================
RCS file: /home/cvs/pari/src/headers/paricom.h,v
retrieving revision 1.43
diff -u -r1.43 paricom.h
--- src/headers/paricom.h	19 Dec 2003 11:05:36 -0000	1.43
+++ src/headers/paricom.h	26 Dec 2003 21:30:12 -0000
@@ -148,12 +148,17 @@
 #  define VERYBIGINT (9223372036854775807L) /* 2^63-1 */
 #  define EXP220 (1099511627776L)          /* 2^40   */
 #  define BIGINT (2147483647)              /* 2^31-1 */
+#  define u_OK_ULONG(p) ((ulong)p <= 3037000493UL)
 #else
 #  define VERYBIGINT (2147483647L) /* 2^31-1 */
 #  define EXP220 (1048576L)       /* 2^20   */
 #  define BIGINT (32767)          /* 2^15-1 */
+#  define u_OK_ULONG(p) ((ulong)p <= 46337UL)
 #endif
 
+/* 2p^2 < 2^BITS_IN_LONG */
+#define OK_ULONG(p) (lgefint(p) == 3 && u_OK_ULONG(p[2]))
+
 #ifndef HAS_EXP2
 #  ifdef __cplusplus
      inline double exp2(double x) {return exp(x*LOG2);}
@@ -321,6 +326,11 @@
 #define gres(x,y) (poldivres((x),(y),ONLY_REM))
 #define FpX_div(x,y,p) (FpX_divres((x),(y),(p), NULL))
 #define FpX_res(x,y,p) (FpX_divres((x),(y),(p), ONLY_REM))
+
+#define Flx_div(x,y,p) (Flx_divres((x),(y),(p), NULL))
+#define Flv_FpV(x)  (vecsmall_vec((x)))
+#define Flx_copy(x) (vecsmall_copy((x)))
+
 #define matpascal(n) matqpascal((n),NULL)
 #define sturm(x) (sturmpart((x),NULL,NULL))
 #define carreparfait(x) (carrecomplet((x),NULL))
Index: src/headers/paridecl.h
===================================================================
RCS file: /home/cvs/pari/src/headers/paridecl.h,v
retrieving revision 1.298
diff -u -r1.298 paridecl.h
--- src/headers/paridecl.h	19 Dec 2003 11:05:36 -0000	1.298
+++ src/headers/paridecl.h	26 Dec 2003 21:30:16 -0000
@@ -20,6 +20,10 @@
 /*******************************************************************/
 BEGINEXTERN
 /* alglin1.c */
+GEN     Flm_deplin(GEN x, ulong p);
+GEN     Flm_inv(GEN x, ulong p);
+GEN     Flm_ker(GEN x, ulong p);
+GEN     Flm_ker_sp(GEN x, ulong p, long deplin);
 GEN     FpV_FpL_mul(GEN x, GEN y, GEN p);
 GEN     FpM_FpV_mul(GEN x, GEN y, GEN p);
 GEN     FpM_deplin(GEN x, GEN p);
@@ -34,7 +38,6 @@
 GEN     FqM_gauss(GEN a, GEN b, GEN T, GEN p);
 GEN     FqM_ker(GEN x, GEN T, GEN p);
 GEN     FqM_suppl(GEN x, GEN T, GEN p);
-GEN     Fq_mul(GEN x, GEN y, GEN T, GEN p);
 GEN     QM_inv(GEN M, GEN dM);
 GEN     ZM_inv(GEN M, GEN dM);
 GEN     _col(GEN x);
@@ -222,6 +225,7 @@
 GEN     mpfact(long n);
 GEN     mpfactr(long n, long prec);
 GEN     mpinvmod(GEN a, GEN m);
+GEN     mpinvmodsafe(GEN a, GEN m);
 GEN     mppgcd(GEN a, GEN b);
 GEN     mpppcm(GEN a, GEN b);
 GEN     mpsqrtmod(GEN a, GEN p);
@@ -1044,6 +1048,10 @@
 GEN     zeropol(long v);
 GEN     zeroser(long v, long prec);
 
+/* groupid.c*/
+
+long    group_ident(GEN G, GEN S);
+
 /* ifactor1.c */
 
 long    BSW_psp(GEN N);
@@ -1052,10 +1060,6 @@
 GEN     plisprime(GEN N, long flag);
 GEN     precprime(GEN n);
 
-/* groupid.c*/
-
-long    group_ident(GEN G, GEN S);
-
 /* init.c */
 
 long    TIMER(pari_timer *T);
@@ -1267,13 +1271,15 @@
 GEN     quotient_group(GEN C, GEN G);
 GEN     vecperm_orbits(GEN v, long n);
 GEN     vecsmall_append(GEN V, long s);
-GEN     vecsmall_prepend(GEN V, long s);
 GEN     vecsmall_const(long n, long c);
+GEN     vecsmall_copy(GEN x);
 int     vecsmall_lexcmp(GEN x, GEN y);
 long    vecsmall_pack(GEN V, long base, long mod);
 int     vecsmall_prefixcmp(GEN x, GEN y);
+GEN     vecsmall_prepend(GEN V, long s);
 void    vecsmall_sort(GEN V);
 GEN     vecsmall_uniq(GEN V);
+GEN     vecsmall_vec(GEN z);
 
 /* polarit1.c */
 
@@ -1319,10 +1325,8 @@
 GEN     rootsold(GEN x, long l);
 GEN     simplefactmod(GEN f, GEN p);
 GEN     swap_polpol(GEN x, long n, long w);
-long    u_FpX_nbfact(GEN z, long p);
-long    u_FpX_nbroots(GEN f, long p);
-GEN     u_pol_to_vec(GEN x, long N);
-GEN     u_vec_to_pol(GEN x);
+long    Flx_nbfact(GEN z, long p);
+long    Flx_nbroots(GEN f, long p);
 GEN     vec_to_pol(GEN x, long v);
 GEN     vecpol_to_mat(GEN v, long n);
 
@@ -1439,11 +1443,20 @@
 void    Fp_intersect(long n,GEN P,GEN Q,GEN l,GEN *SP,GEN *SQ,GEN MA,GEN MB);
 GEN     Fp_inv_isom(GEN S,GEN Tp, GEN p);
 GEN     Fp_isom(GEN P,GEN Q,GEN l);
+GEN     Fq_inv(GEN x, GEN pol, GEN p);
+GEN     Fq_invsafe(GEN x, GEN pol, GEN p);
+GEN     Fq_add(GEN x, GEN y, GEN T/*unused*/, GEN p);
+GEN     Fq_mul(GEN x, GEN y, GEN T, GEN p);
+GEN     Fq_neg(GEN x, GEN T, GEN p);
+GEN     Fq_neg_inv(GEN x, GEN T, GEN p);
+GEN     Fq_pow(GEN x, GEN n, GEN pol, GEN p);
+GEN     Fq_res(GEN x, GEN T, GEN p);
 GEN     FqV_roots_to_pol(GEN V, GEN p, GEN Tp, long v);
+GEN     FqXQ_pow(GEN x, GEN n, GEN S, GEN T, GEN p);
 GEN     FqXV_mul(GEN V, GEN Tp, GEN p);
+GEN     QX_invmod(GEN A, GEN B);
 GEN     ZX_caract(GEN A, GEN B, long v);
 GEN     ZX_disc(GEN x);
-GEN     QX_invmod(GEN A, GEN B);
 int     ZX_is_squarefree(GEN x);
 GEN     ZX_resultant(GEN A, GEN B);
 GEN     ZX_QX_resultant(GEN A, GEN B);
@@ -1465,20 +1478,8 @@
 GEN     unscale_pol(GEN P, GEN h);
 GEN     stopoly(ulong m, ulong p, long v);
 GEN     stopoly_gen(GEN m, GEN p, long v);
-GEN     u_FpXQ_pow(GEN x, GEN n, GEN pol, ulong p);
-GEN     u_FpX_divrem(GEN x, GEN y, ulong p, GEN *pr);
-GEN     u_FpX_rem(GEN x, GEN y, ulong p);
-GEN     u_Fp_FpM(GEN x, ulong p);
-GEN     u_Fp_FpV(GEN x, ulong p);
-GEN     u_Fp_FpX(GEN x, ulong p);
-int     u_FpX_is_squarefree(GEN z, ulong p);
-GEN     u_FpX_normalize(GEN z, ulong p);
-GEN     u_FpX_sub(GEN x, GEN y, ulong p);
-GEN     u_FpX_gcd(GEN a, GEN b, ulong p);
-GEN     u_getpol(long d);
 int     u_pow(int p, int k);
 int     u_val(ulong n, ulong p);
-GEN     u_zeropol(void);
 
 /* rootpol.c */
 
@@ -1549,6 +1550,51 @@
 GEN     thue(GEN thueres, GEN rhs, GEN ne);
 GEN     thueinit(GEN poly, long flag, long prec);
 
+/* smalpol1.c */
+
+GEN     Fl_Flx(ulong x, long sv);
+GEN     Flm_FlxV(GEN x, long sv);
+GEN     Flv_polint(GEN xa, GEN ya, ulong p, long vs);
+GEN     Flv_roots_to_pol(GEN a, ulong p, long vs);
+GEN     Flx_Fl_mul(GEN y, ulong x, ulong p);
+GEN     Flx_Flv_lg(GEN x, long N);
+GEN     Flx_FpX_inplace(GEN z);
+GEN     Flx_add(GEN x, GEN y, ulong p);
+GEN     Flx_addshift(GEN x, GEN y, ulong p, long d);
+GEN     Flx_addspec(GEN x, GEN y, ulong p, long lx, long ly);
+GEN     Flx_deriv(GEN z, ulong p);
+GEN     Flx_div_by_X_x(GEN a, ulong x, ulong p);
+GEN     Flx_divres(GEN x, GEN y, ulong p, GEN *pr);
+ulong   Flx_eval(GEN x, ulong y, ulong p);
+GEN     Flx_extgcd(GEN a, GEN b, ulong p, GEN *ptu, GEN *ptv);
+ulong   Flx_extresultant(GEN a, GEN b, ulong p, GEN *ptU, GEN *ptV);
+GEN     Flx_gcd(GEN a, GEN b, ulong p);
+GEN     Flx_gcd_i(GEN a, GEN b, ulong p);
+GEN     Flx_get(long d);
+int     Flx_is_squarefree(GEN z, ulong p);
+GEN     Flx_mul(GEN x, GEN y, ulong p);
+GEN     Flx_neg(GEN x, ulong p);
+GEN     Flx_normalize(GEN z, ulong p);
+GEN     Flx_polx(long sv);
+GEN     Flx_pow(GEN x, long n, ulong p);
+GEN     Flx_renormalize(GEN x, long l);
+GEN     Flx_res(GEN x, GEN y, ulong p);
+ulong   Flx_resultant(GEN a, GEN b, ulong p);
+GEN     Flx_sqr(GEN x, ulong p);
+GEN     Flx_sub(GEN x, GEN y, ulong p);
+GEN     Flx_zero(long v);
+GEN     Flx_ZX(GEN z);
+GEN     Flx_ZX_inplace(GEN z);
+GEN     Flxq_inv(GEN x,GEN T,ulong p);
+GEN     Flxq_invsafe(GEN x, GEN T, ulong p);
+GEN     Flxq_mul(GEN y,GEN x,GEN pol,ulong p);
+GEN     Flxq_pow(GEN x, GEN n, GEN pol, ulong p);
+GEN     Flxq_sqr(GEN y,GEN pol,ulong p);
+GEN     Z_Flx(GEN x, ulong p, long vs);
+GEN     ZM_Flm(GEN x, ulong p);
+GEN     ZV_Flv(GEN x, ulong p);
+GEN     ZX_Flx(GEN x, ulong p);
+ 
 /* trans1.c */
 
 GEN     Pi2n(long n, long prec);
Index: src/language/es.c
===================================================================
RCS file: /home/cvs/pari/src/language/es.c,v
retrieving revision 1.134
diff -u -r1.134 es.c
--- src/language/es.c	19 Dec 2003 11:05:36 -0000	1.134
+++ src/language/es.c	26 Dec 2003 21:30:27 -0000
@@ -1786,6 +1786,7 @@
   if (typ(g) != t_MAT) { bruti(g,T,1); return; }
 
   r=lg(g); if (r==1 || lg(g[1])==1) { pariputs("[;]\n"); return; }
+  if (typ(g[1]) == t_VECSMALL) { bruti(g,T,1); return; }
   pariputc('\n'); l = lg(g[1]);
   for (i=1; i<l; i++)
   {
Index: src/modules/galois.c
===================================================================
RCS file: /home/cvs/pari/src/modules/galois.c,v
retrieving revision 1.40
diff -u -r1.40 galois.c
--- src/modules/galois.c	19 Dec 2003 11:05:36 -0000	1.40
+++ src/modules/galois.c	26 Dec 2003 21:30:35 -0000
@@ -730,8 +730,9 @@
   a = (GEN)BR->coef[l]; /* fill with random polynomial of degree <= l-1 */
   do
   {
+    a[1]=0;
     for (i=2; i < l+2; i++) a[i] = random_bits(3) + 1;
-    h = normalizepol( small_to_pol(a, 0) );
+    h = Flx_ZX(Flx_renormalize(a,l+2));
   } while (degpol(h) <= 0 || !ZX_is_squarefree(h));
   setvarn(h, v); k = 0;
   (void)ZX_caract_sqf(h, BR->p, &k, v);
@@ -2482,7 +2483,11 @@
     buildroot BR;
     long i;
     GEN z = cgetg(N + 1, t_VEC);
-    for (i = 1; i <= N; i++) z[i] = (long)u_getpol(i-1);
+    for (i = 1; i <= N; i++) 
+    {
+      z[i] = (long)cgetg(i+2,t_VECSMALL);
+      mael(z,i,1)=0;
+    }
     BR.coef = z;
     BR.p = pol;
     BR.pr = (gexpo( cauchy_bound(pol) ) >> TWOPOTBITS_IN_LONG) + prec;
Index: src/modules/nffactor.c
===================================================================
RCS file: /home/cvs/pari/src/modules/nffactor.c,v
retrieving revision 1.148
diff -u -r1.148 nffactor.c
--- src/modules/nffactor.c	19 Dec 2003 11:05:36 -0000	1.148
+++ src/modules/nffactor.c	26 Dec 2003 21:30:41 -0000
@@ -27,7 +27,7 @@
 extern GEN FqX_split_roots(GEN z, GEN T, GEN p, GEN pol);
 extern GEN FqX_split_all(GEN z, GEN T, GEN p);
 extern long FqX_split_by_degree(GEN *pz, GEN u, GEN q, GEN T, GEN p);
-extern long FpX_split_berlekamp(GEN *t, GEN p);
+extern long FpX_split_Berlekamp(GEN *t, GEN p);
 extern GEN get_proj_modT(GEN basis, GEN T, GEN p);
 extern void init_dalloc();
 extern double *dalloc(size_t n);
@@ -143,7 +143,7 @@
   if (!T)
   {
     rep = factmod0(x, p);
-    rep[2] = (long)small_to_vec((GEN)rep[2]);
+    rep[2] = (long)vecsmall_vec((GEN)rep[2]);
   }
   else
   {
@@ -1494,9 +1494,9 @@
     red = modprX(polbase, nf, modpr);
     if (!aT)
     { /* degree 1 */
-      red = u_Fp_FpX(red, pp);
-      if (!u_FpX_is_squarefree(red, pp)) { avma = av2; continue; }
-      anbf = fl? u_FpX_nbroots(red, pp): u_FpX_nbfact(red, pp);
+      red = ZX_Flx(red, pp);
+      if (!Flx_is_squarefree(red, pp)) { avma = av2; continue; }
+      anbf = fl? Flx_nbroots(red, pp): Flx_nbfact(red, pp);
     }
     else
     {
@@ -1567,7 +1567,7 @@
     {
       long d;
       rep = cgetg(dpol + 1, t_VEC); rep[1] = (long)polred;
-      d = FpX_split_berlekamp((GEN*)(rep + 1), L.p);
+      d = FpX_split_Berlekamp((GEN*)(rep + 1), L.p);
       setlg(rep, d + 1);
     }
     T.fact  = gerepilecopy(av2, sort_vecpol(rep, &cmp_pol));
Index: src/modules/thue.c
===================================================================
RCS file: /home/cvs/pari/src/modules/thue.c,v
retrieving revision 1.26
diff -u -r1.26 thue.c
--- src/modules/thue.c	6 Dec 2003 13:44:17 -0000	1.26
+++ src/modules/thue.c	26 Dec 2003 21:30:43 -0000
@@ -1038,7 +1038,7 @@
     if (!Nprimes) { sNx = 1; x = gun; }
     else
     {
-      x = isprincipalfact(bnf, Primes, small_to_vec(x), NULL,
+      x = isprincipalfact(bnf, Primes, vecsmall_vec(x), NULL,
                           nf_FORCE | nf_GEN_IF_PRINCIPAL);
       x = basistoalg(nf, x);
       sNx = signe(gnorm(x));
--- /dev/null	Sat Mar  1 19:52:54 2003
+++ src/basemath/smalpol1.c	Thu Dec 25 09:43:59 2003
@@ -0,0 +1,1091 @@
+/* $Id: polarit3.c,v 1.177 2003/12/19 11:05:36 kb Exp $
+
+Copyright (C) 2000  The PARI group.
+
+This file is part of the PARI/GP package.
+
+PARI/GP is free software; you can redistribute it and/or modify it under the
+terms of the GNU General Public License as published by the Free Software
+Foundation. It is distributed in the hope that it will be useful, but WITHOUT
+ANY WARRANTY WHATSOEVER.
+
+Check the License for details. You should have received a copy of it, along
+with the package; see the file 'COPYING'. If not, write to the Free Software
+Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */
+
+#include "pari.h"
+
+
+/*Not very fast arithmetic with polynomials with small coefficients.*/
+
+/***********************************************************************/
+/**                                                                   **/
+/**               Flx                                                 **/
+/**                                                                   **/
+/***********************************************************************/
+/* Flx objects are defined as follows:
+   They are t_VECSMALL:
+   x[0]: codeword
+   x[1]=variable number (signe is not stored).
+   x[2]=a_0 x[3]=a_1, etc.
+   With 0<=a_i<l
+   where l is a ulong.
+*/
+
+#define lswap(x,y) {long _z=x; x=y; y=_z;}
+#define pswap(x,y) {GEN *_z=x; x=y; y=_z;}
+#define swap(x,y)  {GEN  _z=x; x=y; y=_z;}
+#define swapspec(x,y, nx,ny) {swap(x,y); lswap(nx,ny);}
+
+long Flx_SQR_LIMIT = 6;
+long Flx_MUL_LIMIT = 10;
+
+#define both_odd(x,y) ((x)&(y)&1)
+
+/***********************************************************************/
+/**                                                                   **/
+/**          Conversion from Flx                                      **/
+/**                                                                   **/
+/***********************************************************************/
+
+/* Flx_ZX=zx_ZX*/
+GEN
+Flx_ZX(GEN z)
+{
+  long i;
+  long l=lg(z);
+  GEN x = cgetg(l,t_POL);
+  if (typ(z) != t_VECSMALL) err(bugparier,"Flx_ZX, not a VECSMALL");
+  for (i=2; i<l; i++) x[i] = lstoi(z[i]);
+  x[1] = evalsigne(l-2!=0)| z[1]; return x;
+}
+
+/* Flx_ZX_inplace=zx_ZX_inplace*/
+/* same. In place */
+GEN
+Flx_ZX_inplace(GEN z)
+{
+  long i, l = lg(z);
+  for (i=2; i<l; i++) z[i] = lstoi(z[i]);
+  settyp(z, t_POL); z[1]=evalsigne(l-2!=0)|z[1]; return z;
+}
+
+/*FIXME: should be zx_zv_lg*/
+
+GEN
+Flx_Flv_lg(GEN x, long N)
+{
+  long i, l;
+  GEN z = cgetg(N+1,t_VECSMALL);
+  if (typ(x) != t_VECSMALL) err(typeer,"Flx_Flv_lg");
+  l = lg(x)-1; x++;
+  for (i=1; i<l ; i++) z[i]=x[i];
+  for (   ; i<=N; i++) z[i]=0;
+  return z;
+}
+
+/*FIXME: should be zv_zx*/
+GEN
+Flv_Flx(GEN x, long vs)
+{
+  long i, l=lg(x)+1;
+  GEN z = cgetg(l,t_VECSMALL); z[1]=vs;
+  x--; 
+  for (i=2; i<l ; i++) z[i]=x[i];
+  return z;
+}
+
+/*FIXME: should be zm_zxV*/
+GEN
+Flm_FlxV(GEN x, long sv)
+{
+  long j, lx = lg(x);
+  GEN y = cgetg(lx, t_VEC);
+  for (j=1; j<lx; j++) y[j] = (long)Flv_Flx((GEN)x[j], sv);
+  return y;
+}
+
+/***********************************************************************/
+/**                                                                   **/
+/**          Conversion to Flx                                        **/
+/**                                                                   **/
+/***********************************************************************/
+
+/*sv is a evalvarn()*/
+/*FIXME: should be zero_Flx or something*/
+/*Flx_zero=zx_zero*/
+GEN
+Flx_zero(long sv)
+{
+  GEN x = cgetg(2, t_VECSMALL);
+  x[1] = sv; return x;
+}
+
+/* FIXME: should be polx_zx or something
+ * since it does not take a Flx as input.
+ *
+ * Flx_polx=zx_polx*/
+
+GEN
+Flx_polx(long sv)
+{
+  GEN z = cgetg(4, t_VECSMALL);
+  z[1] = sv;
+  z[2] = 0;
+  z[3] = 1;
+  return z;
+}
+
+
+/*FIXME: Should be z_zx_zero since modular reduction is not performed*/
+GEN
+Fl_Flx(ulong x, long sv)
+{
+  GEN z;
+  if (!x) return Flx_zero(sv);
+  z = cgetg(3, t_VECSMALL);
+  z[1] = sv;
+  z[2] = (long)x; return z;
+}
+
+/* Take an integer and return a scalar polynomial mod p,
+ * with evalvarn=vs */
+
+GEN
+Z_Flx(GEN x, ulong p, long vs)
+{
+  GEN a = cgetg(3, t_VECSMALL);
+  a[1] = vs;
+  a[2] = (long)umodiu(x, p); 
+  return a;
+}
+
+/* return x[0 .. dx] mod p as t_VECSMALL. Assume x a t_POL*/
+GEN
+ZX_Flx(GEN x, ulong p)
+{
+  long i;
+  long lx = lg(x); 
+  GEN a = cgetg(lx, t_VECSMALL);
+  a[1]=((ulong)x[1])&VARNBITS;
+  for (i=2; i<lx; i++) a[i] = (long)umodiu((GEN)x[i], p);
+  return Flx_renormalize(a,lx);
+}
+
+GEN
+ZV_Flv(GEN x, ulong p)
+{
+  long i, n = lg(x);
+  GEN y = cgetg(n,t_VECSMALL);
+  for (i=1; i<n; i++) y[i] = (long)umodiu((GEN)x[i], p);
+  return y;
+}
+
+GEN
+ZM_Flm(GEN x, ulong p)
+{
+  long j,n = lg(x);
+  GEN y = cgetg(n,t_MAT);
+  if (n == 1) return y;
+  for (j=1; j<n; j++) y[j] = (long)ZV_Flv((GEN)x[j], p);
+  return y;
+}
+
+/***********************************************************************/
+/**                                                                   **/
+/**          Basic operation on Flx                                   **/
+/**                                                                   **/
+/***********************************************************************/
+
+
+/*Similar to normalizepol, in place*/
+/*FIXME: should be zx_renormalize*/
+GEN
+Flx_renormalize(GEN /*in place*/ x, long lx)
+{
+  long i;
+  for (i = lx-1; i>1; i--)
+    if (x[i]) break;
+  stackdummy(x + (i+1), lg(x) - (i+1));
+  setlg(x, i+1);
+  return x;
+}
+
+GEN
+Flx_red(GEN z, ulong p)
+{
+  long i,l;
+  GEN x;
+  l = lg(z); x = cgetg(l, t_VECSMALL);
+  for (i=2; i<l; i++)
+    x[i] = z[i]%p;
+  x[1] = z[1]; return Flx_renormalize(x,l);
+}
+
+GEN
+Flx_addspec(GEN x, GEN y, ulong p, long lx, long ly)
+{
+  long i,lz;
+  GEN z;
+
+  if (ly>lx) swapspec(x,y, lx,ly);
+  lz = lx+2; z = cgetg(lz, t_VECSMALL) + 2;
+  for (i=0; i<ly; i++) z[i] = (long)adduumod((ulong)x[i], (ulong)y[i], p);
+  for (   ; i<lx; i++) z[i] = x[i];
+  z -= 2; return Flx_renormalize(z, lz);
+}
+
+GEN
+Flx_add(GEN x, GEN y, ulong p)
+{
+  long i,lz;
+  GEN z; 
+  long lx=lg(x);
+  long ly=lg(y);
+  if (ly>lx) swapspec(x,y, lx,ly);
+  lz = lx; z = cgetg(lz, t_VECSMALL); z[1]=x[1];
+  for (i=2; i<ly; i++) z[i] = (long)adduumod((ulong)x[i], (ulong)y[i], p);
+  for (   ; i<lx; i++) z[i] = x[i];
+  return Flx_renormalize(z, lz);
+}
+
+GEN
+Flx_sub(GEN x, GEN y, ulong p)
+{
+  long i,lz,lx = lg(x), ly = lg(y);
+  GEN z;
+
+  if (ly <= lx)
+  {
+    lz = lx; z = cgetg(lz, t_VECSMALL);
+    for (i=2; i<ly; i++) z[i] = (long)subuumod((ulong)x[i],(ulong)y[i],p);
+    for (   ; i<lx; i++) z[i] = x[i];
+  }
+  else
+  {
+    lz = ly; z = cgetg(lz, t_VECSMALL);
+    for (i=2; i<lx; i++) z[i] = (long)subuumod((ulong)x[i],(ulong)y[i],p);
+    for (   ; i<ly; i++) z[i] = y[i]? (long)(p - y[i]): y[i];
+  }
+  z[1]=x[1]; return Flx_renormalize(z, lz);
+}
+
+GEN
+Flx_neg(GEN x, ulong p)
+{
+  long i, l = lg(x);
+  GEN z = cgetg(l, t_VECSMALL);
+  z[1] = x[1];
+  for (i=2; i<l; i++) z[i] = x[i]? (long)p - x[i]: 0;
+  return z;
+}
+
+GEN
+Flx_neg_inplace(GEN x, ulong p)
+{
+  long i, l = lg(x);
+  for (i=2; i<l; i++)
+    if (x[i]) x[i] = (long)p - x[i];
+  return x;
+}
+
+GEN
+Flx_Fl_mul(GEN y, ulong x, ulong p)
+{
+  GEN z;
+  int i, l;
+  if (!x) return Flx_zero(y[1]);
+  l = lg(y); z = cgetg(l, t_VECSMALL); z[1]=y[1]; 
+  if (HIGHWORD(x | p))
+    for(i=2; i<l; i++) z[i] = (long)muluumod(y[i], x, p);
+  else
+    for(i=2; i<l; i++) z[i] = (y[i] * x) % p;
+  return z;
+}
+
+GEN
+Flx_normalize(GEN z, ulong p)
+{
+  long l = lg(z)-1;
+  ulong p1 = z[l]; /* leading term */
+  if (p1 == 1) return z;
+  return Flx_Fl_mul(z, invumod(p1,p), p);
+}
+
+/* return (x * X^d) + y. Assume d > 0, x > 0 and y >= 0 */
+GEN
+Flx_addshift(GEN x, GEN y, ulong p, long d)
+{
+  GEN xd,yd,zd = (GEN)avma;
+  long a,lz,ny = lg(y)-2, nx = lg(x)-2;
+  long vs = x[1];
+
+  x += 2; y += 2; a = ny-d;
+  if (a <= 0)
+  {
+    lz = (a>nx)? ny+2: nx+d+2;
+    (void)new_chunk(lz); xd = x+nx; yd = y+ny;
+    while (xd > x) *--zd = *--xd;
+    x = zd + a;
+    while (zd > x) *--zd = 0;
+  }
+  else
+  {
+    xd = new_chunk(d); yd = y+d;
+    x = Flx_addspec(x,yd,p, nx,a);
+    lz = (a>nx)? ny+2: lg(x)+d;
+    x += 2; while (xd > x) *--zd = *--xd;
+  }
+  while (yd > y) *--zd = *--yd;
+  *--zd = vs;
+  *--zd = evaltyp(t_VECSMALL) | evallg(lz); return zd;
+}
+
+/* shift polynomial + gerepile */
+/* Do not set evalvarn*/
+static GEN
+Flx_shiftip(pari_sp av, GEN x, long v)
+{
+  long i, lx = lg(x), ly;
+  GEN y;
+  if (v <= 0 || x-2==0) return gerepileupto(av, x);
+  avma = av; ly = lx + v;
+  x += lx; y = new_chunk(ly) + ly;
+  for (i = 2; i<lx; i++) *--y = *--x;
+  for (i = 0; i< v; i++) *--y = 0;
+  *--y = 0;
+  *--y = evaltyp(t_VECSMALL) | evallg(ly); return y;
+}
+
+/* Require OK_ULONG*/
+
+INLINE ulong
+Flx_mullimb_ok(GEN x, GEN y, ulong p, long a, long b)
+{
+  ulong p1 = 0;
+  long i;
+  for (i=a; i<b; i++)
+    if (y[i])
+    {
+      p1 += y[i] * x[-i];
+      if (p1 & HIGHBIT) p1 %= p;
+    }
+  return p1 % p;
+}
+
+INLINE ulong
+Flx_mullimb(GEN x, GEN y, ulong p, long a, long b)
+{
+  ulong p1 = 0;
+  long i;
+  for (i=a; i<b; i++)
+    if (y[i])
+      p1 = adduumod(p1, muluumod(y[i],x[-i],p), p);
+  return p1;
+}
+
+/* assume nx >= ny > 0 */
+GEN
+Flx_mulspec_basecase(GEN x, GEN y, ulong p, long nx, long ny)
+{
+  long i,lz,nz;
+  GEN z;
+
+  lz = nx+ny+1; nz = lz-2;
+  z = cgetg(lz, t_VECSMALL) + 2; /* x:y:z [i] = term of degree i */
+  if (u_OK_ULONG(p))
+  {
+    for (i=0; i<ny; i++)z[i] = (long)Flx_mullimb_ok(x+i,y,p,0,i+1);
+    for (  ; i<nx; i++) z[i] = (long)Flx_mullimb_ok(x+i,y,p,0,ny);
+    for (  ; i<nz; i++) z[i] = (long)Flx_mullimb_ok(x+i,y,p,i-nx+1,ny);
+  }
+  else
+  {
+    for (i=0; i<ny; i++)z[i] = (long)Flx_mullimb(x+i,y,p,0,i+1);
+    for (  ; i<nx; i++) z[i] = (long)Flx_mullimb(x+i,y,p,0,ny);
+    for (  ; i<nz; i++) z[i] = (long)Flx_mullimb(x+i,y,p,i-nx+1,ny);
+  }
+  z -= 2; z[1]=0; return Flx_renormalize(z, lz);
+}
+
+/* fast product (Karatsuba) of polynomials a,b. These are not real GENs, a+2,
+ * b+2 were sent instead. na, nb = number of terms of a, b.
+ * Only c, c0, c1, c2 are genuine GEN.
+ */
+GEN
+Flx_mulspec(GEN a, GEN b, ulong p, long na, long nb)
+{
+  GEN a0,c,c0;
+  long n0, n0a, i, v = 0;
+  pari_sp av;
+
+  while (na && !a[0]) { a++; na--; v++; }
+  while (nb && !b[0]) { b++; nb--; v++; }
+  if (na < nb) swapspec(a,b, na,nb);
+  if (!nb) return Flx_zero(0);
+
+  av = avma;
+  if (nb < Flx_MUL_LIMIT)
+    return Flx_shiftip(av,Flx_mulspec_basecase(a,b,p,na,nb), v);
+  i=(na>>1); n0=na-i; na=i;
+  a0=a+n0; n0a=n0;
+  while (n0a && !a[n0a-1]) n0a--;
+
+  if (nb > n0)
+  {
+    GEN b0,c1,c2;
+    long n0b;
+
+    nb -= n0; b0 = b+n0; n0b = n0;
+    while (n0b && !b[n0b-1]) n0b--;
+    c =  Flx_mulspec(a,b,p,n0a,n0b);
+    c0 = Flx_mulspec(a0,b0,p,na,nb);
+
+    c2 = Flx_addspec(a0,a,p,na,n0a);
+    c1 = Flx_addspec(b0,b,p,nb,n0b);
+
+    c1 = Flx_mul(c1,c2,p);
+    c2 = Flx_add(c0,c,p);
+
+    c2 = Flx_neg_inplace(c2,p);
+    c2 = Flx_add(c1,c2,p);
+    c0 = Flx_addshift(c0,c2 ,p, n0);
+  }
+  else
+  {
+    c  = Flx_mulspec(a,b,p,n0a,nb);
+    c0 = Flx_mulspec(a0,b,p,na,nb);
+  }
+  c0 = Flx_addshift(c0,c,p,n0);
+  return Flx_shiftip(av,c0, v);
+}
+
+GEN Flx_mul(GEN x, GEN y, ulong p)
+{
+ GEN z=Flx_mulspec(x+2,y+2,p, lg(x)-2,lg(y)-2);
+ z[1]=x[1];
+ return z;
+}
+/*FIXME: Assume OK_ULONG*/
+GEN
+Flx_sqrspec_basecase(GEN x, ulong p, long nx)
+{
+  long i,l,lz,nz;
+  ulong p1;
+  GEN z;
+
+  if (!nx) return Flx_zero(x[1]);
+  lz = (nx << 1) + 1, nz = lz-2;
+  z = cgetg(lz, t_VECSMALL) + 2;
+  for (i=0; i<nx; i++)
+  {
+    l = (i+1)>>1;
+    p1 = Flx_mullimb_ok(x,x,p,0,l);
+    p1 <<= 1; if (p1 & HIGHBIT) p1 %= p;
+    if ((i&1) == 0)
+      p1 += x[i>>1] * x[i>>1];
+    z[i] = (long) (p1 % p);
+  }
+  for (  ; i<nz; i++)
+  {
+    l = (i+1)>>1;
+    p1 = Flx_mullimb_ok(x,x,p,i-nx+1,l);
+    p1 <<= 1; if (p1 & HIGHBIT) p1 %= p;
+    if ((i&1) == 0)
+      p1 += x[i>>1] * x[i>>1];
+    z[i] = (long) (p1 % p);
+  }
+  z -= 2; z[1]=x[1]; return Flx_renormalize(z,lz);
+}
+
+GEN
+Flx_2_mul(GEN x, ulong p)
+{
+  long i, l = lg(x);
+  GEN z = cgetg(l, t_VECSMALL);
+  for (i=2; i<l; i++) z[i] = (long)adduumod(x[i], x[i], p);
+  z[1] = x[1]; return z;
+}
+
+GEN
+Flx_sqrspec(GEN a, ulong p, long na)
+{
+  GEN a0,c,c0,c1;
+  long n0, n0a, i, v = 0;
+  pari_sp av;
+
+  while (na && !a[0]) { a++; na--; v += 2; }
+  av = avma;
+  if (na == 0) return Flx_zero(0);
+  if (na < Flx_SQR_LIMIT) 
+    return Flx_shiftip(av, Flx_mulspec_basecase(a,a,p,na,na), v);
+  i=(na>>1); n0=na-i; na=i;
+  a0=a+n0; n0a=n0;
+  while (n0a && !a[n0a-1]) n0a--;
+
+  c = Flx_sqrspec(a,p,n0a);
+  c0= Flx_sqrspec(a0,p,na);
+  if (p == 2) n0 *= 2;
+  else
+  {
+    c1 = Flx_2_mul(Flx_mulspec(a0,a,p, na,n0a), p);
+    c0 = Flx_addshift(c0,c1,p,n0);
+  }
+  c0 = Flx_addshift(c0,c,p,n0);
+  return Flx_shiftip(av,c0,v);
+}
+
+GEN Flx_sqr(GEN x, ulong p)
+{
+  GEN z=Flx_sqrspec(x+2,p, lg(x)-2);
+  z[1]=x[1];
+  return z;
+}
+
+GEN
+Flx_pow(GEN x, long n, ulong p)
+{
+  GEN y = Fl_Flx(1,x[1]), z;
+  long m;
+  if (n == 0) return y;
+  m = n; z = x;
+  for (;;)
+  {
+    if (m&1) y = Flx_mul(y,z, p);
+    m >>= 1; if (!m) return y;
+    z = Flx_sqr(z, p);
+  }
+}
+
+/* separate from Flx_divres for maximal speed. */
+GEN
+Flx_res(GEN x, GEN y, ulong p)
+{ 
+  pari_sp av;
+  GEN z, c;
+  long dx,dy,dz,i,j;
+  ulong p1,inv;
+  long vs=x[1];
+
+  dy = degpol(y); if (!dy) return Flx_zero(x[1]);
+  dx = degpol(x);
+  dz = dx-dy; if (dz < 0) return Flx_copy(x);
+  x += 2; y += 2;
+  inv = y[dy];
+  if (inv != 1UL) inv = invumod(inv,p);
+
+  c = cgetg(dy+3, t_VECSMALL); c[1]=vs; c += 2; av=avma;
+  z = cgetg(dz+3, t_VECSMALL); z[1]=vs; z += 2;
+
+  if (u_OK_ULONG(p))
+  {
+    z[dz] = (inv*x[dx]) % p;
+    for (i=dx-1; i>=dy; --i)
+    {
+      p1 = p - x[i]; /* compute -p1 instead of p1 (pb with ulongs otherwise) */
+      for (j=i-dy+1; j<=i && j<=dz; j++)
+      {
+        p1 += z[j]*y[i-j];
+        if (p1 & HIGHBIT) p1 %= p;
+      }
+      p1 %= p;
+      z[i-dy] = p1? ((p - p1)*inv) % p: 0;
+    }
+    for (i=0; i<dy; i++)
+    {
+      p1 = z[0]*y[i];
+      for (j=1; j<=i && j<=dz; j++)
+      {
+        p1 += z[j]*y[i-j];
+        if (p1 & HIGHBIT) p1 %= p;
+      }
+      c[i] = (long)subuumod((ulong)x[i], p1%p, p);
+    }
+  }
+  else
+  {
+    z[dz] = muluumod(inv, x[dx], p);
+    for (i=dx-1; i>=dy; --i)
+    {
+      p1 = p - x[i]; /* compute -p1 instead of p1 (pb with ulongs otherwise) */
+      for (j=i-dy+1; j<=i && j<=dz; j++)
+        p1 = adduumod(p1, muluumod(z[j],y[i-j],p), p);
+      z[i-dy] = p1? muluumod(p - p1, inv, p): 0;
+    }
+    for (i=0; i<dy; i++)
+    {
+      p1 = muluumod(z[0],y[i],p);
+      for (j=1; j<=i && j<=dz; j++)
+        p1 = adduumod(p1, muluumod(z[j],y[i-j],p), p);
+      c[i] = (long)subuumod((ulong)x[i], p1, p);
+    }
+  }
+  i = dy-1; while (i>=0 && !c[i]) i--;
+  avma=av;
+  return Flx_renormalize(c-2, i+3);
+}
+
+/* as FpX_divres but working only on ulong types. ASSUME pr != ONLY_DIVIDES
+ * if relevant, *pr is the last object on stack */
+GEN
+Flx_divres(GEN x, GEN y, ulong p, GEN *pr)
+{
+  GEN z,q,c;
+  long dx,dy,dz,i,j;
+  ulong p1,inv;
+  long sv=x[1];
+
+  if (pr == ONLY_REM) return Flx_res(x, y, p);
+  dy = degpol(y);
+  if (!dy)
+  {
+    if (y[2] == 1UL)
+      q = Flx_copy(x);
+    else
+      q = Flx_Fl_mul(x, invumod(y[2], p), p);
+    if (pr) *pr = Flx_zero(sv);
+    return q;
+  }
+  dx = degpol(x);
+  dz = dx-dy;
+  if (dz < 0)
+  {
+    q = Flx_zero(sv);
+    if (pr) *pr = Flx_copy(x);
+    return q;
+  }
+  x += 2;
+  y += 2;
+  z = cgetg(dz + 3, t_VECSMALL); z[1] = sv; z += 2;
+  inv = (ulong)y[dy];
+  if (inv != 1UL) inv = invumod(inv,p);
+
+  if (u_OK_ULONG(p))
+  {
+    z[dz] = (long) (inv*x[dx]) % p;
+    for (i=dx-1; i>=dy; --i)
+    {
+      p1 = p - x[i]; /* compute -p1 instead of p1 (pb with ulongs otherwise) */
+      for (j=i-dy+1; j<=i && j<=dz; j++)
+      {
+        p1 += z[j]*y[i-j];
+        if (p1 & HIGHBIT) p1 %= p;
+      }
+      p1 %= p;
+      z[i-dy] = p1? (long) ((p - p1)*inv) % p: 0;
+    }
+  }
+  else
+  {
+    z[dz] = (long)muluumod(inv, x[dx], p);
+    for (i=dx-1; i>=dy; --i)
+    { /* compute -p1 instead of p1 (pb with ulongs otherwise) */
+      p1 = p - (ulong)x[i];
+      for (j=i-dy+1; j<=i && j<=dz; j++)
+        p1 = adduumod(p1, muluumod(z[j],y[i-j],p), p);
+      z[i-dy] = p1? (long)muluumod(p - p1, inv, p): 0;
+    }
+  }
+  q = Flx_renormalize(z-2, dz+3);
+  if (!pr) return q;
+
+  c = cgetg(dy + 3, t_VECSMALL); c[1] = sv; c += 2;
+  if (u_OK_ULONG(p))
+  {
+    for (i=0; i<dy; i++)
+    {
+      p1 = (ulong)z[0]*y[i];
+      for (j=1; j<=i && j<=dz; j++)
+      {
+        p1 += (ulong)z[j]*y[i-j];
+        if (p1 & HIGHBIT) p1 %= p;
+      }
+      c[i] = (long)subuumod((ulong)x[i], p1%p, p);
+    }
+  }
+  else
+  {
+    for (i=0; i<dy; i++)
+    {
+      p1 = muluumod(z[0],y[i],p);
+      for (j=1; j<=i && j<=dz; j++)
+        p1 = adduumod(p1, muluumod(z[j],y[i-j],p), p);
+      c[i] = (long)subuumod((ulong)x[i], p1, p);
+    }
+  }
+  i=dy-1; while (i>=0 && !c[i]) i--;
+  c = Flx_renormalize(c-2, i+3);
+  *pr = c; return q;
+}
+
+GEN
+Flx_deriv(GEN z, ulong p)
+{
+  long i,l = lg(z)-1;
+  GEN x;
+  if (l < 2) l = 2;
+  x = cgetg(l, t_VECSMALL); x[1] = z[1]; z++;
+  if (HIGHWORD(l | p))
+    for (i=2; i<l; i++) x[i] = (long)muluumod((ulong)i-1, (ulong)z[i], p);
+  else
+    for (i=2; i<l; i++) x[i] = ((i-1) * z[i]) % p;
+  return Flx_renormalize(x,l);
+}
+
+/*Do not garbage collect*/
+
+GEN
+Flx_gcd_i(GEN a, GEN b, ulong p)
+{
+  GEN c;
+  if (lg(b) > lg(a)) swap(a, b);
+  while (lg(b)-2)
+  {
+    c = Flx_res(a,b,p);
+    a = b; b = c;
+  }
+  return a;
+}
+
+GEN
+Flx_gcd(GEN a, GEN b, ulong p)
+{
+  pari_sp av = avma;
+  return gerepileupto(av, Flx_gcd_i(a,b,p));
+}
+
+int
+Flx_is_squarefree(GEN z, ulong p)
+{
+  pari_sp av = avma;
+  GEN d = Flx_gcd_i(z, Flx_deriv(z,p) , p);
+  long res= (degpol(d) == 0);
+  avma = av;
+  return res;
+}
+
+GEN
+Flx_extgcd(GEN a, GEN b, ulong p, GEN *ptu, GEN *ptv)
+{
+  GEN q,z,u,v, x = a, y = b;
+
+  u = Flx_zero(a[1]);
+  v = Fl_Flx(1,a[1]); /* v = 1 */
+  while (lg(y)-2)
+  {
+    q = Flx_divres(x,y,p,&z);
+    x = y; y = z; /* (x,y) = (y, x - q y) */
+    z = Flx_sub(u, Flx_mul(q,v, p), p);
+    u = v; v = z; /* (u,v) = (v, u - q v) */
+  }
+  z = Flx_sub(x, Flx_mul(b,u,p), p);
+  z = Flx_div(z,a,p);
+  *ptu = z;
+  *ptv = u; return x;
+}
+
+ulong
+Flx_resultant(GEN a, GEN b, ulong p)
+{
+  long da,db,dc,cnt;
+  ulong lb, res = 1UL;
+  pari_sp av;
+  GEN c;
+
+  if (lg(a)-2==0 || lg(b)-2==0) return 0;
+  da = degpol(a);
+  db = degpol(b);
+  if (db > da)
+  {
+    swapspec(a,b, da,db);
+    if (both_odd(da,db)) res = p-res;
+  }
+  if (!da) return 1; /* = res * a[2] ^ db, since 0 <= db <= da = 0 */
+  cnt = 0; av = avma;
+  while (db)
+  {
+    lb = b[db+2];
+    c = Flx_res(a,b, p);
+    a = b; b = c; dc = degpol(c);
+    if (dc < 0) { avma = av; return 0; }
+
+    if (both_odd(da,db)) res = p - res;
+    if (lb != 1) res = muluumod(res, powuumod(lb, da - dc, p), p);
+    if (++cnt == 4) { cnt = 0; avma = av; }
+    da = db; /* = degpol(a) */
+    db = dc; /* = degpol(b) */
+  }
+  avma = av; return muluumod(res, powuumod(b[2], da, p), p);
+}
+
+/* If resultant is 0, *ptU and *ptU are not set */
+ulong
+Flx_extresultant(GEN a, GEN b, ulong p, GEN *ptU, GEN *ptV)
+{
+  GEN z,q,u,v, x = a, y = b;
+  ulong lb, res = 1UL;
+  pari_sp av = avma;
+  long dx,dy,dz;
+  long vs=a[1];
+
+  if (!signe(x) || !signe(y)) return 0;
+  dx = degpol(x);
+  dy = degpol(y);
+  if (dy > dx)
+  {
+    swap(x,y); lswap(dx,dy); pswap(ptU, ptV);
+    a = x; b = y;
+    if (both_odd(dx,dy)) res = p-res;
+  }
+  u = Flx_zero(vs);
+  v = Fl_Flx(1,vs); /* v = 1 */
+  while (dy)
+  { /* b u = x (a), b v = y (a) */
+    lb = y[dy+2];
+    q = Flx_divres(x,y, p, &z);
+    x = y; y = z; /* (x,y) = (y, x - q y) */
+    dz = degpol(z); if (dz < 0) { avma = av; return 0; }
+    z = Flx_sub(u, Flx_mul(q,v, p), p);
+    u = v; v = z; /* (u,v) = (v, u - q v) */
+
+    if (both_odd(dx,dy)) res = p - res;
+    if (lb != 1) res = muluumod(res, powuumod(lb, dx-dz, p), p);
+    dx = dy; /* = degpol(x) */
+    dy = dz; /* = degpol(y) */
+  }
+  res = muluumod(res, powuumod(y[2], dx, p), p);
+  lb = muluumod(res, invumod(y[2],p), p);
+  v = gerepileupto(av, Flx_Fl_mul(v, lb, p));
+  av = avma;
+  u = Flx_sub(Fl_Flx(res,vs), Flx_mul(b,v,p), p);
+  u = gerepileupto(av, Flx_div(u,a,p)); /* = (res - b v) / a */
+  *ptU = u;
+  *ptV = v; return res;
+}
+
+ulong
+Flx_eval(GEN x, ulong y, ulong p)
+{
+  ulong p1,r;
+  long j, i=lg(x)-1;
+  if (i<=2)
+    return (i==2)? x[2]: 0;
+  p1 = x[i];
+  /* specific attention to sparse polynomials (see poleval)*/
+  if (u_OK_ULONG(p))
+  {
+    for (i--; i>=2; i=j-1)
+    {
+      for (j=i; !x[j]; j--)
+        if (j==2)
+        {
+          if (i != j) y = powuumod(y, i-j+1, p);
+          return (p1 * y) % p;
+        }
+      r = (i==j)? y: powuumod(y, i-j+1, p);
+      p1 = ((p1*r) + x[j]) % p;
+    }
+  }
+  else
+  {
+    for (i--; i>=2; i=j-1)
+    {
+      for (j=i; !x[j]; j--)
+        if (j==2)
+        {
+          if (i != j) y = powuumod(y, i-j+1, p);
+          return muluumod(p1, y, p);
+        }
+      r = (i==j)? y: powuumod(y, i-j+1, p);
+      p1 = adduumod((ulong)x[j], muluumod(p1,r,p), p);
+    }
+  }
+  return p1;
+}
+
+static ulong global_pp;
+static GEN
+_Flx_mul(GEN a, GEN b)
+{
+  return Flx_mul(a,b, global_pp);
+}
+
+/* compute prod (x - a[i]) */
+GEN
+Flv_roots_to_pol(GEN a, ulong p, long vs)
+{
+  long i,k,lx = lg(a);
+  GEN p1,p2;
+  if (lx == 1) return polun[0];
+  p1 = cgetg(lx, t_VEC); global_pp = p;
+  for (k=1,i=1; i<lx-1; i+=2)
+  {
+    p2 = cgetg(5,t_VECSMALL); p1[k++] = (long)p2;
+    p2[1] = vs;
+    p2[2] = muluumod(a[i], a[i+1], p);
+    p2[3] = a[i] + a[i+1];
+    if ((ulong)p2[3] >= p) p2[3] -= p;
+    if (p2[3]) p2[3] = p - p2[3]; /* - (a[i] + a[i+1]) mod p */
+    p2[4] = 1; p2[1] = 0;
+  }
+  if (i < lx)
+  {
+    p2 = cgetg(4,t_POL); p1[k++] = (long)p2;
+    p2[1] = 0;
+    p2[2] = p - a[i];
+    p2[3] = 1;
+  }
+  setlg(p1, k); return divide_conquer_prod(p1, _Flx_mul);
+}
+
+GEN
+Flx_div_by_X_x(GEN a, ulong x, ulong p)
+{
+  long l = lg(a), i;
+  GEN a0, z0;
+  GEN z = cgetg(l-1,t_POL);
+  z[1]=a[1];
+  a0 = a + l-1;
+  z0 = z + l-2; *z0 = *a0--;
+  if (u_OK_ULONG(p))
+  {
+    for (i=l-3; i>1; i--) /* z[i] = (a[i+1] + x*z[i+1]) % p */
+    {
+      ulong t = (*a0-- + x *  *z0--) % p;
+      *z0 = t;
+    }
+  }
+  else
+  {
+    for (i=l-3; i>1; i--)
+    {
+      ulong t = adduumod((ulong)*a0--, muluumod(x, *z0--, p), p);
+      *z0 = (long)t;
+    }
+  }
+  return z;
+}
+
+/* u P(X) + v P(-X) */
+GEN
+Flx_even_odd_comb(GEN P, ulong u, ulong v, ulong p)
+{
+  long i, l = lg(P);
+  GEN y = cgetg(l,t_VECSMALL);
+  y[1]=P[1];
+  for (i=2; i<l; i++)
+  {
+    ulong t = P[i];
+    y[i] = (t == 0)? 0:
+                     (i&1)? muluumod(t, u + (p - v), p)
+                          : muluumod(t, u + v, p);
+  }
+  return Flx_renormalize(y,l);
+}
+
+/* xa, ya = t_VECSMALL */
+GEN
+Flv_polint(GEN xa, GEN ya, ulong p, long vs)
+{
+  GEN T,dP, P = NULL, Q = Flv_roots_to_pol(xa, p, vs);
+  long i, n = lg(xa);
+  ulong inv;
+  pari_sp av = avma;
+  for (i=1; i<n; i++)
+  {
+    if (!ya[i]) continue;
+    T = Flx_div_by_X_x(Q, xa[i], p);
+    inv = invumod(Flx_eval(T,xa[i], p), p);
+    if (i < n-1 && (ulong)(xa[i] + xa[i+1]) == p)
+    {
+      dP = Flx_even_odd_comb(T, muluumod(ya[i],inv,p), muluumod(ya[i+1],inv,p), p);
+      i++; /* x_i = -x_{i+1} */
+    }
+    else
+      dP = Flx_Fl_mul(T, muluumod(ya[i],inv,p), p);
+    P = P? Flx_add(P, dP, p): dP;
+  }
+  return P? gerepileupto(av, P): Flx_zero(vs);
+}
+
+/***********************************************************************/
+/**                                                                   **/
+/**               Flxq                                                **/
+/**                                                                   **/
+/***********************************************************************/
+/* Flxq objects are defined as follows:
+   They are Flx modulo another Flx called q.
+*/
+
+/* Product of y and x in Z/pZ[X]/(pol), as t_VECSMALL. */
+GEN
+Flxq_mul(GEN y,GEN x,GEN pol,ulong p)
+{
+  GEN z = Flx_mul(y,x,p);
+  return Flx_res(z,pol,p);
+}
+
+/* Square of y in Z/pZ[X]/(pol), as t_VECSMALL. */
+GEN
+Flxq_sqr(GEN y,GEN pol,ulong p)
+{
+  GEN z = Flx_sqr(y,p);
+  return Flx_res(z,pol,p);
+}
+
+typedef struct {
+  GEN pol;
+  ulong p;
+} Flxq_muldata;
+
+
+static GEN
+_u_sqr(void *data, GEN x)
+{
+  Flxq_muldata *D = (Flxq_muldata*)data;
+  return Flxq_sqr(x, D->pol, D->p);
+}
+static GEN
+_u_mul(void *data, GEN x, GEN y)
+{
+  Flxq_muldata *D = (Flxq_muldata*)data;
+  return Flxq_mul(x,y, D->pol, D->p);
+}
+
+/* n-Power of x in Z/pZ[X]/(pol), as t_VECSMALL. */
+/* FIXME: assume n > 0 */
+GEN
+Flxq_pow(GEN x, GEN n, GEN pol, ulong p)
+{
+  pari_sp av = avma;
+  Flxq_muldata D;
+  GEN y;
+  D.pol = pol;
+  D.p   = p;
+  y = leftright_pow(x, n, (void*)&D, &_u_sqr, &_u_mul);
+  return gerepileupto(av, y);
+}
+
+/* Inverse of x in Z/pZ[X]/(pol) or NULL if inverse doesn't exist
+ * return lift(1 / (x mod (p,pol)))
+ * not stack clean.
+ * */
+GEN
+Flxq_invsafe(GEN x, GEN T, ulong p)
+{
+  GEN U, V;
+  GEN z = Flx_extgcd(x, T, p, &U, &V);
+  ulong iz;
+  if (degpol(z)) return NULL;
+  iz = invumod ((ulong)z[2], p);
+  return Flx_Fl_mul(U, iz, p);
+}
+
+GEN
+Flxq_inv(GEN x,GEN T,ulong p)
+{
+  pari_sp av=avma;
+  GEN U = Flxq_invsafe(x, T, p);
+  if (!U) err(talker,"non invertible polynomial in FpXQ_inv");
+  return gerepileupto(av, U);
+}