===================================================================
RCS file: /home/cvs/OpenXM_contrib/pari-2.2/src/basemath/Attic/polarit2.c,v
retrieving revision 1.1
retrieving revision 1.2
diff -u -p -r1.1 -r1.2
--- OpenXM_contrib/pari-2.2/src/basemath/Attic/polarit2.c	2001/10/02 11:17:05	1.1
+++ OpenXM_contrib/pari-2.2/src/basemath/Attic/polarit2.c	2002/09/11 07:26:52	1.2
@@ -1,4 +1,4 @@
-/* $Id: polarit2.c,v 1.1 2001/10/02 11:17:05 noro Exp $
+/* $Id: polarit2.c,v 1.2 2002/09/11 07:26:52 noro Exp $
 
 Copyright (C) 2000  The PARI group.
 
@@ -24,32 +24,90 @@ Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 #define swap(a,b) { GEN _x = a; a = b; b = _x; }
 #define lswap(a,b) { long _x = a; a = b; b = _x; }
 #define pswap(a,b) { GEN *_x = a; a = b; b = _x; }
-#define both_even(a,b) ((((a)|(b))&1) == 0)
+#define both_odd(a,b) ((a)&(b)&1)
+#define addshift(x,y) addshiftpol((x),(y),1)
 
-GEN gassoc_proto(GEN f(GEN,GEN),GEN,GEN);
+extern GEN FpXQX_safegcd(GEN P, GEN Q, GEN T, GEN p);
+extern GEN ZY_ZXY_resultant(GEN A, GEN B0, long *lambda);
+extern GEN addshiftpol(GEN x, GEN y, long d);
+extern GEN cauchy_bound(GEN p);
+extern GEN gassoc_proto(GEN f(GEN,GEN),GEN,GEN);
+extern GEN gauss_intern(GEN a, GEN b);
+extern GEN indexpartial(GEN P, GEN DP);
+extern GEN qf_disc(GEN x, GEN y, GEN z);
+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 u_center(ulong u, ulong p, ulong ps2);
+extern void gerepilemanycoeffs2(gpmem_t av, GEN x, int n, GEN y, int o);
 
+GEN matratlift(GEN M, GEN mod, GEN amax, GEN bmax, GEN denom);
+GEN nfgcd(GEN P, GEN Q, GEN nf, GEN den);
+
+/* compute Newton sums S_1(P), ... , S_n(P). S_k(P) = sum a_j^k, a_j root of P
+ * If N != NULL, assume p-adic roots and compute mod N [assume integer coeffs]
+ * If T != NULL, compute mod (T,N) [assume integer coeffs if N != NULL]
+ * If y0!= NULL, precomputed i-th powers, i=1..m, m = length(y0) */
 GEN
-polsym(GEN x, long n)
+polsym_gen(GEN P, GEN y0, long n, GEN T, GEN N)
 {
-  long av1,av2,dx=degpol(x),i,k;
-  GEN s,y,x_lead;
+  long dP=degpol(P), i, k, m;
+  gpmem_t av1, av2;
+  GEN s,y,P_lead;
 
   if (n<0) err(impl,"polsym of a negative n");
-  if (typ(x) != t_POL) err(typeer,"polsym");
-  if (!signe(x)) err(zeropoler,"polsym");
-  y=cgetg(n+2,t_COL); y[1]=lstoi(dx);
-  x_lead=(GEN)x[dx+2]; if (gcmp1(x_lead)) x_lead=NULL;
-  for (k=1; k<=n; k++)
+  if (typ(P) != t_POL) err(typeer,"polsym");
+  if (!signe(P)) err(zeropoler,"polsym");
+  y = cgetg(n+2,t_COL);
+  if (y0)
   {
-    av1=avma; s = (dx>=k)? gmulsg(k,(GEN)x[dx+2-k]): gzero;
-    for (i=1; i<k && i<=dx; i++)
-      s = gadd(s,gmul((GEN)y[k-i+1],(GEN)x[dx+2-i]));
-    if (x_lead) s = gdiv(s,x_lead);
-    av2=avma; y[k+1]=lpile(av1,av2,gneg(s));
+    if (typ(y0) != t_COL) err(typeer,"polsym_gen");
+    m = lg(y0)-1;
+    for (i=1; i<=m; i++) y[i] = lcopy((GEN)y0[i]);
   }
+  else
+  {
+    m = 1;
+    y[1] = lstoi(dP);
+  }
+  P += 2; /* strip codewords */
+
+  P_lead = (GEN)P[dP]; if (gcmp1(P_lead)) P_lead = NULL;
+  if (P_lead)
+  {
+    if (N) P_lead = FpXQ_inv(P_lead,T,N);
+    else if (T) P_lead = QX_invmod(P_lead,T);
+  }
+  for (k=m; k<=n; k++)
+  {
+    av1 = avma; s = (dP>=k)? gmulsg(k,(GEN)P[dP-k]): gzero;
+    for (i=1; i<k && i<=dP; i++)
+      s = gadd(s, gmul((GEN)y[k-i+1],(GEN)P[dP-i]));
+    if (N)
+    {
+      if (T) s = FpX_res(FpX_red(s,N), T, N);
+      else   s = modii(s, N);
+      if (P_lead) s = Fq_mul(s, P_lead, T, N);
+    }
+    else if (T)
+    {
+      s = gres(s, T);
+      if (P_lead) s = gres(gmul(s, P_lead), T);
+    }
+    else
+      if (P_lead) s = gdiv(s, P_lead);
+    av2 = avma; y[k+1] = lpile(av1,av2, gneg(s));
+  }
   return y;
 }
 
+GEN
+polsym(GEN x, long n)
+{
+  return polsym_gen(x, NULL, n, NULL,NULL);
+}
+
 static int (*polcmp_coeff_cmp)(GEN,GEN);
 
 /* assume x and y are polynomials in the same variable whose coeffs can be
@@ -60,10 +118,6 @@ polcmp(GEN x, GEN y)
 {
   long i, lx = lgef(x), ly = lgef(y);
   int fl;
-#if 0 /* for efficiency */
-  if (typ(x) != t_POL || typ(y) != t_POL)
-    err(talker,"not a polynomial in polcmp");
-#endif
   if (lx > ly) return  1;
   if (lx < ly) return -1;
   for (i=lx-1; i>1; i--)
@@ -78,7 +132,7 @@ sort_factor(GEN y, int (*cmp)(GEN,GEN))
 {
   int (*old)(GEN,GEN) = polcmp_coeff_cmp;
   polcmp_coeff_cmp = cmp;
-  (void)sort_factor_gen(y,polcmp);
+  (void)sort_factor_gen(y,&polcmp);
   polcmp_coeff_cmp = old; return y;
 }
 
@@ -86,7 +140,8 @@ sort_factor(GEN y, int (*cmp)(GEN,GEN))
 GEN
 sort_factor_gen(GEN y, int (*cmp)(GEN,GEN))
 {
-  long n,i, av = avma;
+  long n, i;
+  gpmem_t av = avma;
   GEN a,b,A,B,w;
   a = (GEN)y[1]; n = lg(a); A = new_chunk(n);
   b = (GEN)y[2];            B = new_chunk(n);
@@ -96,40 +151,74 @@ sort_factor_gen(GEN y, int (*cmp)(GEN,GEN))
   avma = av; return y;
 }
 
+GEN
+sort_vecpol(GEN a)
+{
+  long n, i;
+  gpmem_t av = avma;
+  GEN A,w;
+  n = lg(a); A = new_chunk(n);
+  w = gen_sort(a, cmp_IND | cmp_C, cmp_pol);
+  for (i=1; i<n; i++) A[i] = a[w[i]];
+  for (i=1; i<n; i++) a[i] = A[i];
+  avma = av; return a;
+}
+
+GEN
+centermodii(GEN x, GEN p, GEN po2)
+{
+  GEN y = modii(x,p);
+  if (cmpii(y,po2) > 0) return subii(y,p);
+  return y;
+}
+
+long
+s_centermod(long x, ulong pp, ulong pps2)
+{
+  long y = x % (long)pp;
+  if (y < 0) y += pp;
+  return u_center(y, pp,pps2);
+}
+
 /* for internal use */
 GEN
 centermod_i(GEN x, GEN p, GEN ps2)
 {
-  long av,i,lx;
-  GEN y,p1;
+  long i, lx;
+  gpmem_t av;
+  GEN y;
 
   if (!ps2) ps2 = shifti(p,-1);
   switch(typ(x))
   {
-    case t_INT:
-      y=modii(x,p);
-      if (cmpii(y,ps2)>0) return subii(y,p);
-      return y;
+    case t_INT: return centermodii(x,p,ps2);
 
-    case t_POL: lx=lgef(x);
-      y=cgetg(lx,t_POL); y[1]=x[1];
+    case t_POL: lx = lgef(x);
+      y = cgetg(lx,t_POL); y[1] = x[1];
       for (i=2; i<lx; i++)
       {
-	av=avma; p1=modii((GEN)x[i],p);
-	if (cmpii(p1,ps2)>0) p1=subii(p1,p);
-	y[i]=lpileupto(av,p1);
+	av = avma;
+	y[i] = lpileuptoint(av, centermodii((GEN)x[i],p,ps2));
       }
+      return normalizepol_i(y, lx);
+
+    case t_COL: lx = lg(x);
+      y = cgetg(lx,t_COL);
+      for (i=1; i<lx; i++) y[i] = (long)centermodii((GEN)x[i],p,ps2);
       return y;
 
-    case t_COL: lx=lg(x);
-      y=cgetg(lx,t_COL);
-      for (i=1; i<lx; i++)
-      {
-	p1=modii((GEN)x[i],p);
-	if (cmpii(p1,ps2)>0) p1=subii(p1,p);
-	y[i]=(long)p1;
-      }
+    case t_MAT: lx = lg(x);
+      y = cgetg(lx,t_MAT);
+      for (i=1; i<lx; i++) y[i] = (long)centermod_i((GEN)x[i],p,ps2);
       return y;
+    
+    case t_VECSMALL: lx = lg(x);
+    {
+      ulong pp = itou(p), pps2 = itou(ps2);
+      y = cgetg(lx,t_VECSMALL);
+      for (i=1; i<lx; i++) y[i] = s_centermod(x[i], pp, pps2);
+      return y;
+    }
   }
   return x;
 }
@@ -141,7 +230,8 @@ centermod(GEN x, GEN p) { return centermod_i(x,p,NULL)
 static GEN
 polidivis(GEN x, GEN y, GEN bound)
 {
-  long dx,dy,dz,i,j,av, vx = varn(x);
+  long dx, dy, dz, i, j, vx = varn(x);
+  gpmem_t av;
   GEN z,p1,y_lead;
 
   dy=degpol(y);
@@ -153,13 +243,13 @@ polidivis(GEN x, GEN y, GEN bound)
   if (gcmp1(y_lead)) y_lead = NULL;
 
   p1 = (GEN)x[dx];
-  z[dz]=y_lead? ldivii(p1,y_lead): licopy(p1);
+  z[dz]=y_lead? (long)diviiexact(p1,y_lead): licopy(p1);
   for (i=dx-1; i>=dy; i--)
   {
-    av=avma; p1=(GEN)x[i];
+    av = avma; p1 = (GEN)x[i];
     for (j=i-dy+1; j<=i && j<=dz; j++)
       p1 = subii(p1, mulii((GEN)z[j],(GEN)y[i-j]));
-    if (y_lead) { p1 = gdiv(p1,y_lead); if (typ(p1)!=t_INT) return NULL; }
+    if (y_lead) p1 = diviiexact(p1,y_lead);
     if (absi_cmp(p1, bound) > 0) return NULL;
     p1 = gerepileupto(av,p1);
     z[i-dy] = (long)p1;
@@ -176,57 +266,17 @@ polidivis(GEN x, GEN y, GEN bound)
     if (!i) break;
   }
   z -= 2;
-  z[1]=evalsigne(1) | evallgef(dz+3) | evalvarn(vx);
+  z[1] = evalsigne(1) | evallgef(dz+3) | evalvarn(vx);
   return z;
 }
 
-static long
-min_deg(long jmax,ulong tabbit[])
-{
-  long j, k = 1, j1 = 2;
-
-  for (j=jmax; j>=0; j--)
-  {
-    for (  ; k<=15; k++)
-    {
-      if (tabbit[j]&j1) return ((jmax-j)<<4)+k;
-      j1<<=1;
-    }
-    k = 0; j1 = 1;
-  }
-  return (jmax<<4)+15;
-}
-
-/* tabkbit is a bit vector (only lowest 32 bits of each word are used
- * on 64bit architecture): reading from right to left, bit i+1 is set iff
- * degree i is attainable from the factorisation mod p.
- *
- * record N modular factors of degree d. */
-static void
-record_factors(long N, long d, long jmax, ulong *tabkbit, ulong *tmp)
-{
-  long n,j, a = d>>4, b = d&15; /* d = 16 a + b */
-  ulong *tmp2 = tmp - a;
-
-  for (n = 1; n <= N; n++)
-  {
-    ulong pro, rem = 0;
-    for (j=jmax; j>=a; j--)
-    {
-      pro = tabkbit[j] << b;
-      tmp2[j] = (pro&0xffff) + rem; rem = pro >> 16;
-    }
-    for (j=jmax-a; j>=0; j--) tabkbit[j] |= tmp[j];
-  }
-}
-
 /***********************************************************************/
 /**                                                                   **/
 /**       QUADRATIC HENSEL LIFT (adapted from V. Shoup's NTL)         **/
 /**                                                                   **/
 /***********************************************************************/
 /* Setup for divide/conquer quadratic Hensel lift
- *  a = set of k t_POL in Z[X] corresponding to factors over Fp
+ *  a = set of k t_POL in Z[X] = factors over Fp (T=NULL) or Fp[Y]/(T)
  *  V = set of products of factors built as follows
  *  1) V[1..k] = initial a
  *  2) iterate:
@@ -238,57 +288,7 @@ record_factors(long N, long d, long jmax, ulong *tabkb
  * link[i] = -j if V[i] = a[j]
  *            j if V[i] = V[j] * V[j+1]
  * Arrays (link, V, W) pre-allocated for 2k - 2 elements */
-#if 0 /* restricted to Fp */
 static void
-BuildTree(GEN link, GEN V, GEN W, GEN a, GEN p)
-{
-  long k = lg(a)-1;
-  long i, j, s, minp, mind;
-
-  for (i=1; i<=k; i++) { V[i] = a[i]; link[i] = -i; }
-
-  for (j=1; j <= 2*k-5; j+=2,i++)
-  {
-    minp = j;
-    mind = degpol(V[j]);
-    for (s=j+1; s<i; s++)
-      if (degpol(V[s]) < mind) { minp = s; mind = degpol(V[s]); }
-
-    lswap(V[j], V[minp]);
-    lswap(link[j], link[minp]);
-
-    minp = j+1;
-    mind = degpol(V[j+1]);
-    for (s=j+2; s<i; s++)
-      if (degpol(V[s]) < mind) { minp = s; mind = degpol(V[s]); }
-
-    lswap(V[j+1], V[minp]);
-    lswap(link[j+1], link[minp]);
-
-    V[i] = (long)FpX_mul((GEN)V[j], (GEN)V[j+1], p);
-    link[i] = j;
-  }
-
-  for (j=1; j <= 2*k-3; j+=2)
-  {
-    GEN d, u, v;
-    d = FpX_extgcd((GEN)V[j], (GEN)V[j+1], p, &u, &v);
-    if (degpol(d) > 0) err(talker, "relatively prime polynomials expected");
-    d = (GEN)d[2];
-    if (!gcmp1(d))
-    {
-      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;
-  }
-}
-#endif
-
-/* same over Fp[Y] / (T) [copy-paste] */
-static void
 BuildTree(GEN link, GEN V, GEN W, GEN a, GEN T, GEN p)
 {
   long k = lg(a)-1;
@@ -340,7 +340,7 @@ BuildTree(GEN link, GEN V, GEN W, GEN a, GEN T, GEN p)
 static void
 HenselLift(GEN V, GEN W, long j, GEN f, GEN T, GEN pd, GEN p0, int noinv)
 {
-  ulong av = avma;
+  gpmem_t av = avma;
   long space = lgef(f) * (lgefint(pd) + lgefint(p0) - 2);
   GEN a2,b2,g,z,s,t;
   GEN a = (GEN)V[j], b = (GEN)V[j+1];
@@ -351,7 +351,7 @@ HenselLift(GEN V, GEN W, long j, GEN f, GEN T, GEN pd,
   (void)new_chunk(space); /* HACK */
   g = gadd(f, gneg_i(gmul(a,b)));
   if (T) g = FpXQX_red(g, T, mulii(p0,pd));
-  g = gdivexact(g, p0); g = FpXQX_red(g, NULL, pd);
+  g = gdivexact(g, p0); if (!T) g = FpXQX_red(g, NULL, pd);
   z = FpXQX_mul(v,g, T,pd);
   t = FpXQX_divres(z,a, T,pd, &s);
   t = gadd(gmul(u,g), gmul(t,b)); t = FpXQX_red(t, T, pd);
@@ -370,7 +370,7 @@ HenselLift(GEN V, GEN W, long j, GEN f, GEN T, GEN pd,
   g = gadd(gneg_i(g), gun);
 
   if (T) g = FpXQX_red(g, T, mulii(p0,pd));
-  g = gdivexact(g, p0); g = FpXQX_red(g, NULL, pd);
+  g = gdivexact(g, p0); if (!T) g = FpXQX_red(g, NULL, pd);
   z = FpXQX_mul(v,g, T,pd);
   t = FpXQX_divres(z,a, T,pd, &s);
   t = gadd(gmul(u,g), gmul(t,b)); t = FpXQX_red(t, T, pd);
@@ -422,7 +422,7 @@ MultiLift(GEN f, GEN a, GEN T, GEN p, long e0, int fla
   l = 0; E[++l] = e;
   while (e > 1) { e = (e+1)/2; E[++l] = e; }
 
-  if (DEBUGLEVEL > 3) timer2();
+  if (DEBUGLEVEL > 3) (void)timer2();
 
   if (flag != 2)
   {
@@ -528,14 +528,15 @@ GEN
 polhensellift(GEN pol, GEN fct, GEN p, long exp)
 {
   GEN p1, p2;
-  long av = avma, i, j, l;
+  long i, j, l;
+  gpmem_t av = avma;
 
   /* we check the arguments */
   if (typ(pol) != t_POL)
     err(talker, "not a polynomial in polhensellift");
   if ((typ(fct) != t_COL && typ(fct) != t_VEC) || (lg(fct) < 3))
     err(talker, "not a factorization in polhensellift");
-  if (typ(p) != t_INT || !isprime(p))
+  if (typ(p) != t_INT || !BSW_psp(p))
     err(talker, "not a prime number in polhensellift");
   if (exp < 1)
     err(talker, "not a positive exponent in polhensellift");
@@ -564,7 +565,7 @@ polhensellift(GEN pol, GEN fct, GEN p, long exp)
   {
     for (i = 1; i <= l; i++)
       for (j = 1; j < i; j++)
-        if (lgef(FpX_gcd((GEN)p1[i], (GEN)p1[j], p)) != 3)
+        if (degpol(FpX_gcd((GEN)p1[i], (GEN)p1[j], p)))
           err(talker, "polhensellift: factors %Z and %Z are not coprime",
                      p1[i], p1[j]);
   }
@@ -579,7 +580,8 @@ polhensellift(GEN pol, GEN fct, GEN p, long exp)
 static GEN
 two_factor_bound(GEN x)
 {
-  long av = avma, i,j, n = lgef(x) - 3;
+  long i, j, n = lgef(x) - 3;
+  gpmem_t av = avma;
   GEN *invbin, c, r = cgetr(3), z;
 
   x += 2; invbin = (GEN*)new_chunk(n+1);
@@ -604,86 +606,138 @@ two_factor_bound(GEN x)
 }
 #endif
 
+/* A | S ==> |a_i| <= binom(d-1, i-1) || S ||_2 + binom(d-1, i) lc(S) */
 static GEN
-uniform_Mignotte_bound(GEN x)
+Mignotte_bound(GEN S)
 {
-  long e, n = degpol(x);
-  GEN p1, N2 = mpsqrt(QuickNormL2(x,DEFAULTPREC));
+  long i, d = degpol(S);
+  GEN lS, C, binlS, bin, N2, p1;
+  
+  N2 = mpsqrt(QuickNormL2(S,DEFAULTPREC));
+  binlS = bin = vecbinome(d-1);
+  lS = leading_term(S);
+  if (!is_pm1(lS)) binlS = gmul(lS, bin);
 
-  p1 = grndtoi(gmul2n(N2, n), &e);
-  if (e>=0) p1 = addii(p1, shifti(gun, e));
-  return p1;
+  /* i = 0 */
+  C = (GEN)binlS[1];
+  /* i = d */
+  p1 = N2; if (gcmp(C, p1) < 0) C = p1;
+  for (i = 1; i < d; i++)
+  {
+    p1 = gadd(gmul((GEN)bin[i], N2), (GEN)binlS[i+1]);
+    if (gcmp(C, p1) < 0) C = p1;
+  }
+  return C;
 }
+/* A | S ==> |a_i|^2 <= 3^{3/2 + d} / (4 \pi d) [P]_2^2,
+ * where [P]_2 is Bombieri's 2-norm */
+static GEN
+Beauzamy_bound(GEN S)
+{
+  const long prec = DEFAULTPREC;
+  long i, d = degpol(S);
+  GEN bin, lS, s, C;
+  bin = vecbinome(d);
 
+  s = realzero(prec);
+  for (i=0; i<=d; i++)
+  {
+    GEN c = (GEN)S[i+2];
+    if (gcmp0(c)) continue;
+    /* s += P_i^2 / binomial(d,i) */
+    s = addrr(s, gdiv(itor(sqri(c), prec), (GEN)bin[i+1]));
+  }
+  /* s = [S]_2^2 */
+  C = gpow(stoi(3), dbltor(1.5 + d), prec); /* 3^{3/2 + d} */
+  C = gdiv(gmul(C, s), gmulsg(4*d, mppi(prec)));
+  lS = absi(leading_term(S));
+  return mulir(lS, mpsqrt(C));
+}
+
+static GEN
+factor_bound(GEN S)
+{
+  gpmem_t av = avma;
+  GEN a = Mignotte_bound(S);
+  GEN b = Beauzamy_bound(S);
+  if (DEBUGLEVEL>2)
+  {
+    fprintferr("Mignotte bound: %Z\n",a);
+    fprintferr("Beauzamy bound: %Z\n",b);
+  }
+  return gerepileupto(av, ceil_safe(gmin(a, b)));
+}
+
+#if 0
 /* all factors have coeffs less than the bound */
 static GEN
 all_factor_bound(GEN x)
 {
   long i, n = lgef(x) - 3;
   GEN t, z = gzero;
-  for (i=2; i<=n+2; i++) z = addii(z,sqri((GEN)x[i]));
+  for (i=2; i<=n+2; i++) z = addii(z, sqri((GEN)x[i]));
   t = absi((GEN)x[n+2]);
   z = addii(t, addsi(1, racine(z)));
   z = mulii(z, binome(stoi(n-1), n>>1));
   return shifti(mulii(t,z),1);
 }
+#endif
 
-/* x defined mod p^a, return mods ( (x - mods(x, p^b)) / p^b , p^(b-a) )  */
-static GEN
-TruncTrace(GEN x, GEN pb, GEN pa_b, GEN pa_bs2, GEN pbs2)
-{
-  GEN r, q = dvmdii(x, pb, &r);
-  if (cmpii(r,  pbs2) > 0) q = addis(q,1);
-  if (pa_bs2 && cmpii(q,pa_bs2) > 0) q = subii(q,pa_b);
-  return q;
-}
-
 /* Naive recombination of modular factors: combine up to maxK modular
  * factors, degree <= klim and divisible by hint
  *
  * target = polynomial we want to factor
  * famod = array of modular factors.  Product should be congruent to
  * target/lc(target) modulo p^a
- * All rational factors are bounded by p^b, b <= a and p^(b-a) < 2^(BIL-1) */
+ * For true factors: S1,S2 <= p^b, with b <= a and p^(b-a) < 2^31 */
 static GEN
-cmbf(GEN target, GEN famod, GEN p, long b, long a,
+cmbf(GEN pol, GEN famod, GEN bound, GEN p, long a, long b,
      long maxK, long klim,long hint)
 {
   long K = 1, cnt = 1, i,j,k, curdeg, lfamod = lg(famod)-1;
-  long spa_b, spa_bs2;
-  GEN lc, lctarget, pa = gpowgs(p,a), pas2 = shifti(pa,-1);
-  GEN trace    = cgetg(lfamod+1, t_VECSMALL);
+  long spa_b, spa_bs2, Sbound;
+  GEN lc, lcpol, pa = gpowgs(p,a), pas2 = shifti(pa,-1);
+  GEN trace1   = cgetg(lfamod+1, t_VECSMALL);
+  GEN trace2   = cgetg(lfamod+1, t_VECSMALL);
   GEN ind      = cgetg(lfamod+1, t_VECSMALL);
-  GEN degpol      = cgetg(lfamod+1, t_VECSMALL);
+  GEN degpol   = cgetg(lfamod+1, t_VECSMALL);
   GEN degsofar = cgetg(lfamod+1, t_VECSMALL);
   GEN listmod  = cgetg(lfamod+1, t_COL);
   GEN fa       = cgetg(lfamod+1, t_COL);
   GEN res = cgetg(3, t_VEC);
-  GEN bound = all_factor_bound(target);
 
   if (maxK < 0) maxK = lfamod-1;
 
-  lc = absi(leading_term(target));
-  lctarget = gmul(lc,target);
+  lc = absi(leading_term(pol));
+  if (is_pm1(lc)) lc = NULL;
+  lcpol = lc? gmul(lc,pol): pol;
 
   {
-    GEN pa_b,pa_bs2,pb,pbs2;
-    pa_b = gpowgs(p, a-b); /* make sure p^(a-b) < 2^(BIL-1) */
-    while (is_bigint(pa_b)) { b++; pa_b = divii(pa_b, p); }
+    GEN pa_b,pa_bs2,pb, lc2 = lc? sqri(lc): NULL;
+
+    pa_b = gpowgs(p, a-b); /* < 2^31 */
     pa_bs2 = shifti(pa_b,-1);
-    pb= gpowgs(p, b); pbs2 = shifti(pb,-1);
+    pb= gpowgs(p, b);
     for (i=1; i <= lfamod; i++)
     {
-      GEN T, p1 = (GEN)famod[i];
-      degpol[i] = degpol(p1);
-      if (!gcmp1(lc))
-        T = modii(mulii(lc, (GEN)p1[degpol[i]+1]), pa);
-      else
-        T = (GEN)p1[degpol[i]+1]; /* d-1 term */
-      trace[i] = itos( TruncTrace(T, pb,pa_b,NULL,pbs2) );
+      GEN T1,T2, P = (GEN)famod[i];
+      long d = degpol(P);
+
+      degpol[i] = d; P += 2;
+      T1 = (GEN)P[d-1];/* = - S_1 */
+      T2 = sqri(T1);
+      if (d > 1) T2 = subii(T2, shifti((GEN)P[d-2],1));
+      T2 = modii(T2, pa); /* = S_2 Newton sum */
+      if (lc)
+      {
+        T1 = modii(mulii(lc, T1), pa);
+        T2 = modii(mulii(lc2,T2), pa);
+      }
+      trace1[i] = itos(diviiround(T1, pb));
+      trace2[i] = itos(diviiround(T2, pb));
     }
-    spa_b   =   pa_b[2]; /* < 2^(BIL-1) */
-    spa_bs2 = pa_bs2[2]; /* < 2^(BIL-1) */
+    spa_b   =   pa_b[2]; /* < 2^31 */
+    spa_bs2 = pa_bs2[2]; /* < 2^31 */
   }
   degsofar[0] = 0; /* sentinel */
 
@@ -691,9 +745,10 @@ cmbf(GEN target, GEN famod, GEN p, long b, long a,
    * 1 <= ind[i] <= lfamod */
 nextK:
   if (K > maxK || 2*K > lfamod) goto END;
-  if (DEBUGLEVEL > 4)
+  if (DEBUGLEVEL > 3)
     fprintferr("\n### K = %d, %Z combinations\n", K,binome(stoi(lfamod), K));
   setlg(ind, K+1); ind[1] = 1;
+  Sbound = ((K+1)>>1);
   i = 1; curdeg = degpol[ind[1]];
   for(;;)
   { /* try all combinations of K factors */
@@ -704,73 +759,87 @@ nextK:
     }
     if (curdeg <= klim && curdeg % hint == 0) /* trial divide */
     {
-      GEN y, q, cont, list;
-      ulong av;
+      GEN y, q, list;
+      gpmem_t av;
       long t;
 
-      /* d - 1 test,  overflow is not a problem (correct mod 2^BIL) */
-      for (t=trace[ind[1]],i=2; i<=K; i++)
-        t = addssmod(t, trace[ind[i]], spa_b);
+      /* d - 1 test */
+      for (t=trace1[ind[1]],i=2; i<=K; i++)
+        t = addssmod(t, trace1[ind[i]], spa_b);
       if (t > spa_bs2) t -= spa_b;
-      if (labs(t) > ((K+1)>>1))
+      if (labs(t) > Sbound)
       {
         if (DEBUGLEVEL>6) fprintferr(".");
         goto NEXT;
       }
-      av = avma;
+      /* d - 2 test */
+      for (t=trace2[ind[1]],i=2; i<=K; i++)
+        t = addssmod(t, trace2[ind[i]], spa_b);
+      if (t > spa_bs2) t -= spa_b;
+      if (labs(t) > Sbound)
+      {
+        if (DEBUGLEVEL>6) fprintferr("|");
+        goto NEXT;
+      }
 
+      av = avma;
       /* check trailing coeff */
       y = lc;
       for (i=1; i<=K; i++)
-        y = centermod_i(mulii(y, gmael(famod,ind[i],2)), pa, pas2);
-      if (!signe(y) || resii((GEN)lctarget[2], y) != gzero)
       {
-        if (DEBUGLEVEL>6) fprintferr("T");
+        GEN q = constant_term((GEN)famod[ind[i]]);
+        if (y) q = mulii(y, q);
+        y = centermod_i(q, pa, pas2);
+      }
+      if (!signe(y) || resii((GEN)lcpol[2], y) != gzero)
+      {
+        if (DEBUGLEVEL>3) fprintferr("T");
         avma = av; goto NEXT;
       }
       y = lc; /* full computation */
       for (i=1; i<=K; i++)
-        y = centermod_i(gmul(y, (GEN)famod[ind[i]]), pa, pas2);
+      {
+        GEN q = (GEN)famod[ind[i]];
+        if (y) q = gmul(y, q);
+        y = centermod_i(q, pa, pas2);
+      }
 
       /* y is the candidate factor */
-      if (! (q = polidivis(lctarget,y,bound)) )
+      if (! (q = polidivis(lcpol,y,bound)) )
       {
-        if (DEBUGLEVEL>6) fprintferr("*");
+        if (DEBUGLEVEL>3) fprintferr("*");
         avma = av; goto NEXT;
       }
       /* found a factor */
-      cont = content(y);
-      if (signe(leading_term(y)) < 0) cont = negi(cont);
-      y = gdiv(y, cont);
-
       list = cgetg(K+1, t_VEC);
       listmod[cnt] = (long)list;
       for (i=1; i<=K; i++) list[i] = famod[ind[i]];
 
+      y = primpart(y);
       fa[cnt++] = (long)y;
-      /* fix up target */
-      target = gdiv(q, leading_term(y));
+      /* fix up pol */
+      pol = q;
+      if (lc) pol = gdivexact(pol, leading_term(y));
       for (i=j=k=1; i <= lfamod; i++)
       { /* remove used factors */
         if (j <= K && i == ind[j]) j++;
         else
         {
           famod[k] = famod[i];
-          trace[k] = trace[i];
+          trace1[k] = trace1[i];
+          trace2[k] = trace2[i];
           degpol[k] = degpol[i]; k++;
         }
       }
       lfamod -= K;
       if (lfamod < 2*K) goto END;
       i = 1; curdeg = degpol[ind[1]];
-      bound = all_factor_bound(target);
-      lc = absi(leading_term(target));
-      lctarget = gmul(lc,target);
+      bound = factor_bound(pol);
+      if (lc) lc = absi(leading_term(pol));
+      lcpol = lc? gmul(lc,pol): pol;
       if (DEBUGLEVEL > 2)
-      {
-        fprintferr("\n"); msgtimer("to find factor %Z",y);
-        fprintferr("remaining modular factor(s): %ld\n", lfamod);
-      }
+        fprintferr("\nfound factor %Z\nremaining modular factor(s): %ld\n",
+                   y, lfamod);
       continue;
     }
 
@@ -786,13 +855,13 @@ NEXT:
     }
   }
 END:
-  if (degpol(target) > 0)
+  if (degpol(pol) > 0)
   { /* leftover factor */
-    if (signe(leading_term(target)) < 0) target = gneg_i(target);
+    if (signe(leading_term(pol)) < 0) pol = gneg_i(pol);
 
     setlg(famod, lfamod+1);
     listmod[cnt] = (long)dummycopy(famod);
-    fa[cnt++] = (long)target;
+    fa[cnt++] = (long)pol;
   }
   if (DEBUGLEVEL>6) fprintferr("\n");
   setlg(listmod, cnt); setlg(fa, cnt);
@@ -822,7 +891,7 @@ factor_quad(GEN x, GEN res, long *ptcnt)
 long
 logint(GEN B, GEN y, GEN *ptq)
 {
-  ulong av = avma;
+  gpmem_t av = avma;
   long e,i,fl;
   GEN q,pow2, r = y;
 
@@ -858,78 +927,35 @@ END:
   if (ptq) *ptq = gerepileuptoint(av, icopy(r)); else avma = av;
   return e;
 }
-    
+
 /* recombination of modular factors: van Hoeij's algorithm */
 
-/* compute Newton sums of P (i-th powers of roots, i=1..n)
- * If N != NULL, assume p-adic roots and compute mod N [assume integer coeffs]
- * If y0!= NULL, precomputed i-th powers, i=1..m, m = length(y0) */
-GEN
-polsym_gen(GEN P, GEN y0, long n, GEN N)
-{
-  long av1,av2,dP=degpol(P),i,k,m;
-  GEN s,y,P_lead;
-
-  if (n<0) err(impl,"polsym of a negative n");
-  if (typ(P) != t_POL) err(typeer,"polsym");
-  if (!signe(P)) err(zeropoler,"polsym");
-  y = cgetg(n+2,t_COL);
-  if (y0)
-  {
-    if (typ(y0) != t_COL) err(typeer,"polsym_gen");
-    m = lg(y0)-1;
-    for (i=1; i<=m; i++) y[i] = lcopy((GEN)y0[i]);
-  }
-  else
-  {
-    m = 1;
-    y[1] = lstoi(dP);
-  }
-  P += 2; /* strip codewords */
-
-  P_lead=(GEN)P[dP]; if (gcmp1(P_lead)) P_lead=NULL;
-  if (N && P_lead)
-    if (!invmod(P_lead,N,&P_lead))
-      err(talker,"polsyn: non-invertible leading coeff: %Z", P_lead);
-  for (k=m; k<=n; k++)
-  {
-    av1=avma; s = (dP>=k)? gmulsg(k,(GEN)P[dP-k]): gzero;
-    for (i=1; i<k && i<=dP; i++)
-      s = gadd(s, gmul((GEN)y[k-i+1],(GEN)P[dP-i]));
-    if (N)
-    {
-      s = modii(s, N);
-      if (P_lead) { s = mulii(s,P_lead); s = modii(s,N); }
-    }
-    else
-      if (P_lead) s = gdiv(s,P_lead);
-    av2=avma; y[k+1]=lpile(av1,av2,gneg(s));
-  }
-  return y;
-}
-
-/* return integer y such that all roots of P are less than y */
+/* return integer y such that all |a| <= y if P(a) = 0 */
 static GEN
-root_bound(GEN P)
+root_bound(GEN P0)
 {
-  GEN P0 = dummycopy(P), lP = absi(leading_term(P)), x,y,z;
-  long k,d = degpol(P);
+  GEN Q = dummycopy(P0), lP = absi(leading_term(Q)), P,x,y,z;
+  long k, d = degpol(Q);
 
-  setlgef(P0, d+2); /* P = lP x^d + P0 */
-  P = P0+2; /* strip codewords */
+  setlgef(Q, d+2); /* P = lP x^d + Q */
+  P = Q+2; /* strip codewords */
   for (k=0; k<d; k++) P[k] = labsi((GEN)P[k]);
 
-  x = y = gun;
+  k = gexpo(cauchy_bound(P0)) - 5;
+  if (k < 0) k = 0;
+  x = y = shifti(gun, k);
   for (k=0; ; k++)
   {
-    if (cmpii(poleval(P0,y), mulii(lP, gpowgs(y, d))) < 0) break;
+    gpmem_t av = avma;
+    if (cmpii(poleval(Q,y), mulii(lP, gpowgs(y, d))) < 0) break;
+    avma = av;
     x = y; y = shifti(y,1);
   }
-  for(;;)
+  for(k=0; ; k++)
   {
     z = shifti(addii(x,y), -1);
-    if (egalii(x,z)) break;
-    if (cmpii(poleval(P0,z), mulii(lP, gpowgs(z, d))) < 0)
+    if (egalii(x,z) || k == 20) break;
+    if (cmpii(poleval(Q,z), mulii(lP, gpowgs(z, d))) < 0)
       y = z;
     else
       x = z;
@@ -937,44 +963,16 @@ root_bound(GEN P)
   return y;
 }
 
-extern GEN gscal(GEN x,GEN y);
-extern GEN sqscal(GEN x);
-
-/* return Gram-Schmidt orthogonal basis (f) associated to (e), B is the
- * vector of the (f_i . f_i)*/
-GEN
-gram_schmidt(GEN e, GEN *ptB)
+static GEN
+ZM_HNFimage(GEN x)
 {
-  long i,j,lx = lg(e);
-  GEN f = dummycopy(e), B, iB;
-
-  B = cgetg(lx, t_VEC);
-  iB= cgetg(lx, t_VEC);
-
-  for (i=1;i<lx;i++)
-  {
-    GEN p1 = NULL;
-    ulong av;
-    B[i] = (long)sqscal((GEN)f[i]);
-    iB[i]= linv((GEN)B[i]); av = avma;
-    for (j=1; j<i; j++)
-    {
-      GEN mu = gmul(gscal((GEN)e[i],(GEN)f[j]), (GEN)iB[j]);
-      GEN p2 = gmul(mu, (GEN)f[j]);
-      p1 = p1? gadd(p1,p2): p2;
-    }
-    p1 = p1? gerepileupto(av, gsub((GEN)e[i], p1)): (GEN)e[i];
-    f[i] = (long)p1;
-  }
-  *ptB = B; return f;
+  return (lg(x) > 50)? hnflll_i(x, NULL, 1): hnfall_i(x, NULL, 1);
 }
 
-extern GEN hnfall_i(GEN A, GEN *ptB, long remove);
-
 GEN
 special_pivot(GEN x)
 {
-  GEN t, H = hnfall_i(x, NULL, 1);
+  GEN t, H = ZM_HNFimage(x);
   long i,j, l = lg(H), h = lg(H[1]);
   for (i=1; i<h; i++)
   {
@@ -992,422 +990,480 @@ special_pivot(GEN x)
   return H;
 }
 
-/* x matrix: compute a bound for | \sum e_i x_i | ^ 2, e_i = 0,1 */
+/* B from lllint_i: return [ |b_i^*|^2, i ] */
+GEN
+GS_norms(GEN B, long prec)
+{
+  long i, l = lg(B);
+  GEN v = gmul(B, realun(prec));
+  l--; setlg(v, l);
+  for (i=1; i<l; i++)
+    v[i] = ldivrr((GEN)v[i+1], (GEN)v[i]);
+  return v;
+}
+
 static GEN
-my_norml2(GEN x)
+check_factors(GEN P, GEN M_L, GEN bound, GEN famod, GEN pa)
 {
-  long i,j, co = lg(x), li;
-  GEN S = gzero;
-  if (co == 1) return S;
-  li = lg(x[1]);
-  for (i=1; i<li; i++)
+  long i, j, r, n0;
+  GEN pol = P, lcpol, lc, list, piv, y, pas2;
+
+  piv = special_pivot(M_L);
+  if (!piv) return NULL;
+  if (DEBUGLEVEL>3) fprintferr("special_pivot output:\n%Z\n",piv);
+
+  pas2 = shifti(pa,-1);
+  r  = lg(piv)-1;
+  n0 = lg(piv[1])-1;
+  list = cgetg(r+1, t_COL);
+  lc = absi(leading_term(pol));
+  if (is_pm1(lc)) lc = NULL;
+  lcpol = lc? gmul(lc, pol): pol;
+  for (i=1; i<r; i++)
   {
-    GEN p = gzero;
-    GEN m = gzero;
-    for (j=1; j<co; j++)
+    GEN c = (GEN)piv[i];
+    if (DEBUGLEVEL) fprintferr("LLL_cmbf: checking factor %ld\n",i);
+
+    y = lc;
+    for (j=1; j<=n0; j++)
+      if (signe(c[j]))
+      {
+        GEN q = (GEN)famod[j];
+        if (y) q = gmul(y, q);
+        y = centermod_i(q, pa, pas2);
+      }
+
+    /* y is the candidate factor */
+    if (! (pol = polidivis(lcpol,y,bound)) ) return NULL;
+    y = primpart(y);
+    if (lc)
     {
-      GEN u = gcoeff(x,i,j);
-      long s = gsigne(u);
-      if      (s < 0) m = gadd(m,u);
-      else if (s > 0) p = gadd(p,u);
+      pol = gdivexact(pol, leading_term(y));
+      lc = absi(leading_term(pol));
     }
-    if (m != gzero) m = gneg(m);
-    if (gcmp(p,m) > 0) m = p;
+    lcpol = lc? gmul(lc, pol): pol;
+    list[i] = (long)y;
+  }
+  y = primpart(pol);
+  list[i] = (long)y; return list;
+}
 
-    S = gadd(S, gsqr(m));
+GEN
+LLL_check_progress(GEN Bnorm, long n0, GEN m, int final,
+                   pari_timer *T, long *ti_LLL)
+{
+  GEN B, norm, u;
+  long i, R;
+
+  if (DEBUGLEVEL>2) (void)TIMER(T);
+  u = lllint_i(m, final? 1000: 4, 0, NULL, NULL, &B);
+  if (DEBUGLEVEL>2) *ti_LLL += TIMER(T);
+  norm = GS_norms(B, DEFAULTPREC);
+  for (R=lg(m)-1; R > 0; R--)
+    if (cmprr((GEN)norm[R], Bnorm) < 0) break;
+  for (i=1; i<=R; i++) setlg(u[i], n0+1);
+  if (R <= 1)
+  {
+    if (!R) err(bugparier,"LLL_cmbf [no factor]");
+    return NULL; /* irreducible */
   }
-  return S;
+  setlg(u, R+1); return u;
 }
 
-/* each entry < 2^e */
-static long
-init_padic_prec(long e, int BitPerFactor, long n0, double LOGp2)
+static ulong
+next2pow(ulong a)
 {
-/* long b, goodb = (long) ((e - 0.175 * n0 * n0)  * LOGp2); */
-  long b, goodb = (long) (((e - BitPerFactor*n0)) * LOGp2);
-  b = (long) ((e -     32)*LOGp2); if (b < goodb) goodb = b;
-  /* b = (long) ((e - Max*32)*LOGp2); if (b > goodb) goodb = b; */
-  return goodb;
+  ulong b = 1;
+  while (b < a) b <<= 1;
+  return b;
 }
 
-extern GEN sindexrank(GEN x);
-extern GEN vconcat(GEN Q1, GEN Q2);
-extern GEN gauss_intern(GEN a, GEN b);
-
 /* Recombination phase of Berlekamp-Zassenhaus algorithm using a variant of
  * van Hoeij's knapsack
  *
- * P = monic squarefree in Z[X].
+ * P = squarefree in Z[X].
  * famod = array of (lifted) modular factors mod p^a
- * bound = Mignotte bound for the size of divisors of P (sor the sup norm)
- * previously recombined all set of factors with less than rec elts
- */
-GEN
+ * bound = Mignotte bound for the size of divisors of P (for the sup norm)
+ * previously recombined all set of factors with less than rec elts */
+static GEN
 LLL_cmbf(GEN P, GEN famod, GEN p, GEN pa, GEN bound, long a, long rec)
 {
-  long s = 2; /* # of traces added at each step */
-  long BitPerFactor = 3; /* nb bits in p^(a-b) / modular factor */
-  long i,j,C,r,S,tmax,n0,n,dP = degpol(P);
-  double logp = log(gtodouble(p)), LOGp2 = LOG2/logp;
-  double b0 = log((double)dP*2) / logp;
-  double k = gtodouble(glog(root_bound(P), DEFAULTPREC)) / logp;
-  GEN y, T, T2, TT, BL, m, u, norm, target, M, piv, list;
-  GEN run = realun(DEFAULTPREC);
-  ulong av,av2, lim;
-  int did_recompute_famod = 0;
+  const long N0 = 1; /* # of traces added at each step */
+  double BitPerFactor = 0.5; /* nb bits in p^(a-b) / modular factor */
+  long i,j,tmax,n0,C, dP = degpol(P);
+  double logp = log((double)itos(p)), LOGp2 = LOG2/logp;
+  double b0 = log((double)dP*2) / logp, logBr;
+  GEN lP, Br, Bnorm, Tra, T2, TT, CM_L, m, list, ZERO;
+  gpmem_t av, av2, lim;
+  long ti_LLL = 0, ti_CF  = 0;
+  pari_timer ti, ti2, TI;
 
-  n0 = n = r = lg(famod) - 1;
-  list = cgetg(n0+1, t_COL);
+  if (DEBUGLEVEL>2) (void)TIMER(&ti);
 
+  lP = absi(leading_term(P));
+  if (is_pm1(lP)) lP = NULL;
+  Br = root_bound(P);
+  if (lP) Br = gmul(lP, Br);
+  logBr = gtodouble(glog(Br, DEFAULTPREC)) / logp;
+
+  n0 = lg(famod) - 1;
+  C = (long)ceil( sqrt(N0 * n0 / 4.) ); /* > 1 */
+  Bnorm = dbltor(n0 * (C*C + N0*n0/4.) * 1.00001);
+  ZERO = zeromat(n0, N0);
+
   av = avma; lim = stack_lim(av, 1);
   TT = cgetg(n0+1, t_VEC);
-  T  = cgetg(n0+1, t_MAT);
+  Tra  = cgetg(n0+1, t_MAT);
   for (i=1; i<=n0; i++)
   {
-    TT[i] = 0;
-    T [i] = lgetg(s+1, t_COL);
+    TT[i]  = 0;
+    Tra[i] = lgetg(N0+1, t_COL);
   }
-  BL = idmat(n0);
-  /* tmax = current number of traces used (and computed so far)
-   * S = number of traces used at the round's end = tmax + s */
-  for(tmax = 0;; tmax = S)
+  CM_L = gscalsmat(C, n0);
+  /* tmax = current number of traces used (and computed so far) */
+  for (tmax = 0;; tmax += N0)
   {
-    GEN pas2, pa_b, BE;
-    long b, goodb;
-    double Nx;
+    long b, bmin, bgood, delta, tnew = tmax + N0, r = lg(CM_L)-1;
+    GEN M_L, q, CM_Lp, oldCM_L;
+    int first = 1;
+    
+    bmin = (long)ceil(b0 + tnew*logBr);
+    if (DEBUGLEVEL>2)
+      fprintferr("\nLLL_cmbf: %ld potential factors (tmax = %ld, bmin = %ld)\n",
+                 r, tmax, bmin);
 
-    if (DEBUGLEVEL>3)
-      fprintferr("LLL_cmbf: %ld potential factors (tmax = %ld)\n", r, tmax);
-    av2 = avma;
-    if (tmax > 0)
-    { /* bound small vector in terms of a modified L2 norm of a
-       * left inverse of BL */
-      GEN z = gauss_intern(BL,NULL); /* 1/BL */
-      if (!z) /* not maximal rank */
-      {
-        avma = av2;
-        BL = hnfall_i(BL,NULL,1);
-        r = lg(BL)-1; z = invmat(BL);
-        av2 = avma;
-      }
-      Nx = gtodouble(my_norml2(gmul(run, z)));
-      avma = av2;
-    }
-    else
-      Nx = (double)n0; /* first time: BL = id */
-    C = (long)sqrt(s*n0*n0/4. / Nx);
-    if (C == 0) C = 1; /* frequent after a few iterations */
-    M = dbltor((Nx * C*C + s*n0*n0/4.) * 1.00001);
-
-    S = tmax + s;
-    b = (long)ceil(b0 + S*k);
-    if (a <= b)
+    /* compute Newton sums (possibly relifting first) */
+    if (a <= bmin)
     {
-      a = (long)ceil(b + 3*s*k) + 1; /* enough for 3 more rounds */
+      a = (long)ceil(bmin + 3*N0*logBr) + 1; /* enough for 3 more rounds */
+      a = next2pow(a);
+
       pa = gpowgs(p,a);
       famod = hensel_lift_fact(P,famod,NULL,p,pa,a);
-      /* recompute old Newton sums to new precision */
-      for (i=1; i<=n0; i++)
-        TT[i] = (long)polsym_gen((GEN)famod[i], NULL, tmax, pa);
-      did_recompute_famod = 1;
+      for (i=1; i<=n0; i++) TT[i] = 0;
     }
     for (i=1; i<=n0; i++)
     {
-      GEN p1 = (GEN)T[i];
-      GEN p2 = polsym_gen((GEN)famod[i], (GEN)TT[i], S, pa);
+      GEN p1 = (GEN)Tra[i];
+      GEN p2 = polsym_gen((GEN)famod[i], (GEN)TT[i], tnew, NULL, pa);
       TT[i] = (long)p2;
       p2 += 1+tmax; /* ignore traces number 0...tmax */
-      for (j=1; j<=s; j++) p1[j] = p2[j];
+      for (j=1; j<=N0; j++) p1[j] = p2[j];
+      if (lP)
+      { /* make Newton sums integral */
+        GEN lPpow = gpowgs(lP, tmax);
+        for (j=1; j<=N0; j++)
+        {
+          lPpow = mulii(lPpow,lP);
+          p1[j] = lmulii((GEN)p1[j], lPpow);
+        }
+      }
     }
 
+    /* compute truncation parameter */
+    if (DEBUGLEVEL>2) { TIMERstart(&ti2); TIMERstart(&TI); }
+    oldCM_L = CM_L;
     av2 = avma;
-    T2 = gmod( gmul(T, BL), pa );
-    goodb = init_padic_prec(gexpo(T2), BitPerFactor, n0, LOGp2);
-    if (goodb > b) b = goodb;
-    if (DEBUGLEVEL>2)
-      fprintferr("LLL_cmbf: (a, b) = (%ld, %ld), expo(T) = %ld\n",
-                 a,b,gexpo(T2));
-    pa_b = gpowgs(p, a-b);
-    {
-      GEN pb = gpowgs(p, b);
-      GEN pa_bs2 = shifti(pa_b,-1);
-      GEN pbs2   = shifti(pb,-1);
-      for (i=1; i<=r; i++)
-      {
-        GEN p1 = (GEN)T2[i];
-        for (j=1; j<=s; j++)
-          p1[j] = (long)TruncTrace((GEN)p1[j],pb,pa_b,pa_bs2,pbs2);
-      }
+    delta = b = 0; /* -Wall */
+AGAIN:
+    M_L = gdivexact(CM_L, stoi(C));
+    T2 = centermod( gmul(Tra, M_L), pa );
+    if (first)
+    { /* initialize lattice, using few p-adic digits for traces */
+      double t = gexpo(T2) - max(32, BitPerFactor*r);
+      bgood = (long) (t * LOGp2);
+      b = max(bmin, bgood);
+      delta = a - b;
     }
-    if (gcmp0(T2)) { avma = av2; continue; }
+    else
+    { /* add more p-adic digits and continue reduction */
+      long b0 = gexpo(T2) * LOGp2;
+      if (b0 < b) b = b0;
+      b = max(b-delta, bmin);
+      if (b - delta/2 < bmin) b = bmin; /* near there. Go all the way */
+    }
 
-    BE = cgetg(s+1, t_MAT);
-    for (i=1; i<=s; i++)
+    q = gpowgs(p, b);
+    m = vconcat( CM_L, gdivround(T2, q) );
+    if (first)
     {
-      BE[i] = (long)zerocol(r+s);
-      coeff(BE, i+r, i) = (long)pa_b;
+      GEN P1 = gscalmat(gpowgs(p, a-b), N0);
+      first = 0;
+      m = concatsp( m, vconcat(ZERO, P1) );
+      /*     [ C M_L        0    ]
+       * m = [                   ]   square matrix
+       *     [  T2'  p^(a-b) I_s ]   T2' = Tra * M_L  truncated
+       */
     }
-    m = concatsp( vconcat( gscalmat(stoi(C), r), T2 ), BE );
-    /*     [  C        0      ]
-     * m = [                  ]   square matrix
-     *     [ T2   p^(a-b) I_S ]   T2 = T * BL  truncated
-     */
-    u = lllgramintern(gram_matrix(m), 4, 0, 0);
-    m = gmul(m,u);
-    (void)gram_schmidt(gmul(run,m), &norm);
-    for (i=r+s; i>0; i--)
-      if (cmprr((GEN)norm[i], M) < 0) break;
-    if (i > r)
-    { /* no progress */
-      avma = av2; BitPerFactor += 2;
-      if (DEBUGLEVEL>2)
-        fprintferr("LLL_cmbf: increasing BitPerFactor = %ld\n", BitPerFactor);
-#if 0
-      s++; for (i=1; i<=n; i++) T[i] = lgetg(s+1, t_COL);
-      if (DEBUGLEVEL>2) fprintferr("LLL_cmbf: increasing s = %ld\n", s);
-#endif
-      continue;
-    }
 
-    n = r; r = i;
-    if (r <= 1)
+    CM_L = LLL_check_progress(Bnorm, n0, m, b == bmin, /*dbg:*/ &ti, &ti_LLL);
+    if (DEBUGLEVEL>2)
+      fprintferr("LLL_cmbf: (a,b) =%4ld,%4ld; r =%3ld -->%3ld, time = %ld\n",
+                 a,b, lg(m)-1, CM_L? lg(CM_L)-1: 1, TIMER(&TI));
+    if (!CM_L) { list = _col(P); break; }
+    if (b > bmin) 
     {
-      if (r == 0) err(bugparier,"LLL_cmbf [no factor]");
-      if (DEBUGLEVEL>2) fprintferr("LLL_cmbf: 1 factor\n");
-      list[1] = (long)P; setlg(list,2); return list;
+      CM_L = gerepilecopy(av2, CM_L);
+      goto AGAIN;
     }
+    if (DEBUGLEVEL>2) msgTIMER(&ti2, "for this block of traces");
 
-    setlg(u, r+1);
-    for (i=1; i<=r; i++) setlg(u[i], n+1);
-    BL = gerepileupto(av2, gmul(BL, u));
-    if (low_stack(lim, stack_lim(av,1)))
+    i = lg(CM_L) - 1;
+    if (i == r && gegal(CM_L, oldCM_L))
     {
-      GEN *gptr[4]; gptr[0]=&BL; gptr[1]=&TT; gptr[2]=&T; gptr[3]=&famod;
-      if(DEBUGMEM>1) err(warnmem,"LLL_cmbf");
-      gerepilemany(av, gptr, did_recompute_famod? 4: 3);
+      CM_L = oldCM_L;
+      avma = av2; continue;
     }
-    if (r*rec >= n0) continue;
 
-    av2 = avma;
-    piv = special_pivot(BL);
-    if (DEBUGLEVEL) fprintferr("special_pivot output:\n%Z\n",piv);
-    if (!piv) { avma = av2; continue; }
-
-    pas2 = shifti(pa,-1); target = P;
-    for (i=1; i<=r; i++)
+    CM_Lp = FpM_image(CM_L, stoi(27449)); /* inexpensive test */
+    if (lg(CM_Lp) != lg(CM_L))
     {
-      GEN p1 = (GEN)piv[i];
-      if (DEBUGLEVEL) fprintferr("LLL_cmbf: checking factor %ld\n",i);
+      if (DEBUGLEVEL>2) fprintferr("LLL_cmbf: rank decrease\n");
+      CM_L = ZM_HNFimage(CM_L);
+    }
 
-      y = gun;
-      for (j=1; j<=n0; j++)
-        if (signe(p1[j]))
-          y = centermod_i(gmul(y, (GEN)famod[j]), pa, pas2);
-
-      /* y is the candidate factor */
-      if (! (target = polidivis(target,y,bound)) ) break;
-      list[i] = (long)y;
+    if (i <= r && i*rec < n0)
+    {
+      if (DEBUGLEVEL>2) (void)TIMER(&ti);
+      list = check_factors(P, M_L, bound, famod, pa);
+      if (DEBUGLEVEL>2) ti_CF += TIMER(&ti);
+      if (list) break;
+      CM_L = gerepilecopy(av2, CM_L);
     }
-    if (i > r)
+    if (low_stack(lim, stack_lim(av,1)))
     {
-      if (DEBUGLEVEL>2) fprintferr("LLL_cmbf: %ld factors\n", r);
-      setlg(list,r+1); return list;
+      if(DEBUGMEM>1) err(warnmem,"LLL_cmbf");
+      gerepileall(av, 5, &CM_L, &TT, &Tra, &famod, &pa);
     }
   }
+  if (DEBUGLEVEL>2)
+    fprintferr("* Time LLL: %ld\n* Time Check Factor: %ld\n",ti_LLL,ti_CF);
+  return list;
 }
 
-extern GEN primitive_pol_to_monic(GEN pol, GEN *ptlead);
-
-/* P(hx), in place. Assume P in Z[X], h in Z */
-void
-rescale_pol_i(GEN P, GEN h)
+/* Return P(h * x) */
+GEN
+unscale_pol(GEN P, GEN h)
 {
-  GEN hi = gun;
   long i, l = lgef(P);
+  GEN hi = gun, Q = cgetg(l, t_POL);
+  Q[1] = P[1];
+  Q[2] = lcopy((GEN)P[2]);
   for (i=3; i<l; i++)
   {
-    hi = mulii(hi,h);
-    P[i] = lmulii((GEN)P[i], hi);
+    hi = gmul(hi,h);
+    Q[i] = lmul((GEN)P[i], hi);
   }
+  return Q;
 }
 
-/* P(hx) in Fp[X], in place. Assume P in Z[X], h in Z */
-void
-FpX_rescale_i(GEN P, GEN h, GEN p)
+/* Return h^degpol(P) P(x / h) */
+GEN
+rescale_pol(GEN P, GEN h)
 {
-  GEN hi = gun;
   long i, l = lgef(P);
-  for (i=3; i<l; i++)
+  GEN Q = cgetg(l,t_POL), hi = gun;
+  Q[l-1] = P[l-1];
+  for (i=l-2; i>=2; i--)
   {
+    hi = gmul(hi,h);
+    Q[i] = lmul((GEN)P[i], hi);
+  }
+  Q[1] = P[1]; return Q;
+}
+
+/* as above over Fp[X] */
+GEN
+FpX_rescale(GEN P, GEN h, GEN p)
+{
+  long i, l = lgef(P);
+  GEN Q = cgetg(l,t_POL), hi = gun;
+  Q[l-1] = P[l-1];
+  for (i=l-2; i>=2; i--)
+  {
     hi = modii(mulii(hi,h), p);
-    P[i] = lmodii(mulii((GEN)P[i], hi), p);
+    Q[i] = lmodii(mulii((GEN)P[i], hi), p);
   }
+  Q[1] = P[1]; return Q;
 }
 
+/* Find a,b minimal such that A < q^a, B < q^b, 1 << q^(a-b) < 2^31 */
+int
+cmbf_precs(GEN q, GEN A, GEN B, long *pta, long *ptb, GEN *qa, GEN *qb)
+{
+  long a,b,amin,d = (long)(31 * LOG2/gtodouble(glog(q,DEFAULTPREC)) - 1e-5);
+  int fl = 0;
+
+  b = logint(B, q, qb);
+  amin = b + d;
+  if (gcmp(gpowgs(q, amin), A) <= 0)
+  {
+    a = logint(A, q, qa);
+    b = a - d; *qb = gpowgs(q, b);
+  }
+  else
+  { /* not enough room */
+    a = amin;  *qa = gpowgs(q, a);
+    fl = 1;
+  }
+  if (DEBUGLEVEL > 3) {
+    fprintferr("S_2   bound: %Z^%ld\n", q,b);
+    fprintferr("coeff bound: %Z^%ld\n", q,a);
+  }
+  *pta = a;
+  *ptb = b; return fl;
+}
+
 /* use van Hoeij's knapsack algorithm */
 static GEN
-combine_factors(GEN a, GEN famod, GEN p, long klim, long hint)
+combine_factors(GEN target, GEN famod, GEN p, long klim, long hint)
 {
-  GEN B = uniform_Mignotte_bound(a), res,y,lt,L,pe,pE,listmod,p1;
-  long i,E,e,l, maxK = 3, nft = lg(famod)-1;
+  GEN la, B, A, res, L, pa, pb, listmod;
+  long a,b, l, maxK = 3, nft = lg(famod)-1, n = degpol(target);
+  pari_timer T;
 
-  e = logint(B, p, &pe);
-  if (DEBUGLEVEL > 4) fprintferr("Mignotte bound: %Z^%ld\n", p,e);
+  A = factor_bound(target);
 
-  { /* overlift for the d-1 test */
-    int over = (int) (32 * LOG2 / log((double)itos(p)));
-    pE = mulii(pe, gpowgs(p,over)); E = e + over;
-  }
-  famod = hensel_lift_fact(a,famod,NULL,p,pE,E);
+  la = absi(leading_term(target));
+  B = mulsi(n, sqri(gmul(la, root_bound(target)))); /* = bound for S_2 */
 
+  (void)cmbf_precs(p, A, B, &a, &b, &pa, &pb);
+
+  if (DEBUGLEVEL>2) (void)TIMER(&T);
+  famod = hensel_lift_fact(target,famod,NULL,p,pa,a);
   if (nft < 11) maxK = -1; /* few modular factors: try all posibilities */
-  else
-  {
-    lt = leading_term(a);
-    if (!is_pm1(lt) && nft < 13) maxK = -1;
-  }
-  L = cmbf(a, famod, p, e, E, maxK, klim, hint);
+  if (DEBUGLEVEL>2) msgTIMER(&T, "Hensel lift (mod %Z^%ld)", p,a);
+  L = cmbf(target, famod, A, p, a, b, maxK, klim, hint);
+  if (DEBUGLEVEL>2) msgTIMER(&T, "Naive recombination");
 
   res     = (GEN)L[1];
-  listmod = (GEN)L[2]; l = lg(listmod);
-  famod = (GEN)listmod[l-1];
+  listmod = (GEN)L[2]; l = lg(listmod)-1;
+  famod = (GEN)listmod[l];
   if (maxK >= 0 && lg(famod)-1 > 2*maxK)
   {
-    a = (GEN)res[l-1];
-    lt = leading_term(a);
-    if (signe(lt) < 0) { a = gneg_i(a); lt = leading_term(a); }
+    if (l!=1) A = factor_bound((GEN)res[l]);
     if (DEBUGLEVEL > 4) fprintferr("last factor still to be checked\n");
-    if (gcmp1(lt))
-      lt = NULL;
-    else
-    {
-      GEN invlt, invLT;
-      if (DEBUGLEVEL > 4) fprintferr("making it monic\n");
-      a = primitive_pol_to_monic(a, &lt);
-      B = uniform_Mignotte_bound(a);
-      e = logint(B, p, &pe);
+    L = LLL_cmbf((GEN)res[l], famod, p, pa, A, a, maxK);
+    if (DEBUGLEVEL>2) msgTIMER(&T,"Knapsack");
+    /* remove last elt, possibly unfactored. Add all new ones. */
+    setlg(res, l); res = concatsp(res, L);
+  }
+  return res;
+}
 
-      invlt = mpinvmod(lt,p);
-      for (i = 1; i<lg(famod); i++)
-      { /* renormalize modular factors */
-        p1 = (GEN)famod[i]; FpX_rescale_i(p1, invlt, p);
-        invLT = powmodulo(lt, stoi(degpol(p1)), p);
-        famod[i] = (long)FpX_Fp_mul(p1, invLT, p);
-      }
-      famod = hensel_lift_fact(a,famod,NULL,p,pe,e);
-    }
-    setlg(res, l-1); /* remove last elt (possibly unfactored) */
-    L = LLL_cmbf(a, famod, p, pe, B, e, maxK);
-    if (lt)
+#define u_FpX_div(x,y,p) u_FpX_divrem((x),(y),(p),NULL)
+
+extern GEN u_FpX_Kmul(GEN a, GEN b, ulong p, long na, long nb);
+extern GEN u_pol_to_vec(GEN x, long N);
+extern GEN u_FpM_FpX_mul(GEN x, GEN y, ulong p);
+#define u_FpX_mul(x,y, p) u_FpX_Kmul(x+2,y+2,p, lgef(x)-2,lgef(y)-2)
+
+static GEN
+u_FpM_Frobenius(GEN u, ulong p)
+{
+  long j,N = degpol(u);
+  GEN v,w,Q;
+  if (DEBUGLEVEL > 7) (void)timer2();
+  Q = cgetg(N+1,t_MAT); 
+  Q[1] = (long)vecsmall_const(N, 0);
+  coeff(Q,1,1) = 1;
+  w = v = u_FpXQ_pow(u_Fp_FpX(polx[varn(u)], p), utoi(p), u, p);
+  for (j=2; j<=N; j++)
+  {
+    Q[j] = (long)u_pol_to_vec(w, N);
+    if (j < N)
     {
-      for (i=1; i<lg(L); i++)
-      {
-        y = (GEN)L[i]; rescale_pol_i(y, lt);
-        L[i] = (long)primitive_part(y, &p1);
-      }
+      gpmem_t av = avma;
+      w = gerepileupto(av, u_FpX_rem(u_FpX_mul(w, v, p), u, p));
     }
-    res = concatsp(res, L);
   }
-  return res;
+  if (DEBUGLEVEL > 7) msgtimer("frobenius");
+  return Q;
 }
 
-extern long split_berlekamp(GEN *t, GEN pp, GEN pps2);
-#define u_FpX_div(x,y,p) u_FpX_divrem((x),(y),(p),(0),NULL)
-
-/* assume degree(a) > 0, a(0) != 0, and a squarefree */
+/* Assume a squarefree, degree(a) > 0, a(0) != 0 */
 static GEN
-squff(GEN a, long hint)
+DDF(GEN a, long hint)
 {
-  GEN PolX,lead,res,prime,primes2,famod,y,g,z,w,*tabd,*tabdnew;
-  long da = degpol(a), va = varn(a);
-  long klim,chosenp,p,nfacp,lbit,j,d,e,np,nmax,nf,nft;
-  ulong *tabbit, *tabkbit, *tmp, av = avma;
-  byteptr pt=diffptr;
-  const int MAXNP = max(5, (long)sqrt(da));
+  GEN MP, PolX, lead, prime, famod, z, w;
+  const long da = degpol(a);
+  long chosenp, p, nfacp, d, e, np, nmax;
+  gpmem_t av = avma, av1;
+  byteptr pt = diffptr;
+  const int MAXNP = max(5, (long)sqrt((double)da));
+  long ti = 0;
+  pari_timer T;
 
+  if (DEBUGLEVEL>2) (void)TIMER(&T);
   if (hint <= 0) hint = 1;
-  if (DEBUGLEVEL > 2) timer2();
-  lbit = (da>>4)+1; nmax = da+1; klim = da>>1;
+  if (DEBUGLEVEL > 2) (void)timer2();
+  nmax = da+1;
   chosenp = 0;
-  tabd = NULL;
-  tabdnew = (GEN*)new_chunk(nmax);
-  tabbit  = (ulong*)new_chunk(lbit);
-  tabkbit = (ulong*)new_chunk(lbit);
-  tmp     = (ulong*)new_chunk(lbit);
   prime = icopy(gun);
-  lead = (GEN)a[da+2]; PolX = u_Fp_FpX(polx[0],0, 2);
-  for (p = np = 0; np < MAXNP; )
+  lead = (GEN)a[da+2]; if (gcmp1(lead)) lead = NULL;
+  PolX = u_Fp_FpX(polx[0], 2);
+  av1 = avma;
+  for (p = np = 0; np < MAXNP; avma = av1)
   {
-    ulong av0 = avma;
-    p += *pt++; if (!*pt) err(primer1);
-    if (!smodis(lead,p)) continue;
-    z = u_Fp_FpX(a,0, p);
-    if (!u_FpX_is_squarefree(z, p)) { avma = av0; continue ; }
+    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;
 
-    for (j=0; j<lbit-1; j++) tabkbit[j] = 0;
-    tabkbit[j] = 1;
+    MP = u_FpM_Frobenius(z, p);
+
     d = 0; e = da; nfacp = 0;
     w = PolX; prime[2] = p;
     while (d < (e>>1))
     {
       long lgg;
-      d++; w = u_FpXQ_pow(w, prime, z, p);
+      GEN g;
+      /* 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);
       lgg = degpol(g);
-      if (lgg == 0) g = NULL;
-      else
-      {
-	z = u_FpX_div(z, g, p); e = degpol(z);
-        w = u_FpX_rem(w, z, p);
-        lgg /= d; nfacp += lgg;
-        if (DEBUGLEVEL>5)
-          fprintferr("   %3ld factor%s of degree %3ld\n", lgg, lgg==1?"":"s",d);
-	record_factors(lgg, d, lbit-1, tabkbit, tmp);
-      }
-      tabdnew[d] = g;
+      if (!lgg) continue;
+
+      e -= lgg;
+      nfacp += lgg/d;
+      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);
     }
-    if (e > 0)
+    if (e)
     {
       if (DEBUGLEVEL>5) fprintferr("   %3ld factor of degree %3ld\n",1,e);
-      tabdnew[e] = z; nfacp++;
-      record_factors(1, e, lbit-1, tabkbit, tmp);
+      nfacp++;
     }
-
-    if (np)
-      for (j=0; j<lbit; j++) tabbit[j] &= tabkbit[j];
-    else
-      for (j=0; j<lbit; j++) tabbit[j] = tabkbit[j];
-    if (DEBUGLEVEL > 4)
+    if (DEBUGLEVEL>4)
       fprintferr("...tried prime %3ld (%-3ld factor%s). Time = %ld\n",
                   p, nfacp, nfacp==1?"": "s", timer2());
-    if (min_deg(lbit-1,tabbit) > klim) { avma = av; return _col(a); }
     if (nfacp < nmax)
     {
-      nmax = nfacp; tabd = tabdnew; chosenp = p;
-      for (j=d+1; j<e; j++) tabd[j] = NULL;
-      tabdnew = (GEN*)new_chunk(da);
+      if (nfacp == 1) { avma = av; return _col(a); }
+      nmax = nfacp; chosenp = p;
+      if (nmax < 5) break; /* very few factors. Enough */
     }
-    else avma = av0;
     np++;
   }
-  prime[2] = chosenp; primes2 = shifti(prime, -1);
-  nf = nmax; nft = 1;
-  y = cgetg(nf+1,t_COL); famod = cgetg(nf+1,t_COL);
-  for (d = 1; nft <= nf; d++)
+  prime[2] = chosenp;
+  famod = cgetg(nmax+1,t_COL);
+  famod[1] = (long)(lead? FpX_normalize(a, prime): FpX_red(a, prime));
+  if (nmax != FpX_split_berlekamp((GEN*)(famod+1),prime))
+    err(bugparier,"DDF: wrong numbers of factors");
+  if (DEBUGLEVEL>2)
   {
-    g = tabd[d]; 
-    if (g)
-    {
-      long n = degpol(g)/d;
-      famod[nft] = (long)small_to_pol(u_FpX_normalize(g, chosenp),va);
-      if (n > 1 && n != split_berlekamp((GEN*)(famod+nft),prime,primes2))
-        err(bugparier,"squff: wrong numbers of factors");
-      nft += n;
-    }
+    if (DEBUGLEVEL>4) msgtimer("splitting mod p = %ld",chosenp);
+    ti = TIMER(&T);
+    fprintferr("Time setup: %ld\n", ti);
   }
-  if (DEBUGLEVEL > 4) msgtimer("splitting mod p = %ld",chosenp);
-  res = combine_factors(a, famod, prime, da-1, hint);
-  return gerepilecopy(av, res);
+  z = combine_factors(a, famod, prime, da-1, hint);
+  if (DEBUGLEVEL>2)
+    fprintferr("Total Time: %ld\n===========\n", ti + TIMER(&T));
+  return gerepilecopy(av, z);
 }
 
 /* A(X^d) --> A(X) */
@@ -1442,10 +1498,11 @@ gdeflate(GEN x, long v, long d)
   long i, lx, tx = typ(x);
   GEN z;
   if (is_scalar_t(tx)) return gcopy(x);
+  if (d <= 0) err(talker,"need positive degree in gdeflate");
   if (tx == t_POL)
   {
     long vx = varn(x);
-    ulong av;
+    gpmem_t av;
     if (vx < v)
     {
       lx = lgef(x);
@@ -1499,13 +1556,15 @@ polinflate(GEN x0, long d)
   return y;
 }
 
+/* Distinct Degree Factorization (deflating first)
+ * Assume x squarefree, degree(x) > 0, x(0) != 0 */
 GEN
-squff2(GEN x, long hint)
+DDF2(GEN x, long hint)
 {
   GEN L;
   long m;
   x = poldeflate(x, &m);
-  L = squff(x, hint);
+  L = DDF(x, hint);
   if (m > 1)
   {
     GEN e, v, fa = decomp(stoi(m));
@@ -1525,69 +1584,87 @@ squff2(GEN x, long hint)
     {
       GEN L2 = cgetg(1,t_VEC);
       for (i=1; i < lg(L); i++)
-        L2 = concatsp(L2, squff(polinflate((GEN)L[i], v[k]), hint));
+        L2 = concatsp(L2, DDF(polinflate((GEN)L[i], v[k]), hint));
       L = L2;
     }
   }
   return L;
 }
 
-/* Factor x in Z[t]. Assume all factors have degree divisible by hint */
+/* SquareFree Factorization. f = prod P^e, all e distinct, in Z[X] (char 0
+ * would be enough, if modulargcd --> ggcd). Return (P), set *ex = (e) */
 GEN
-factpol(GEN x, long hint)
+ZX_squff(GEN f, GEN *ex)
 {
-  GEN fa,p1,d,t,v,w, y = cgetg(3,t_MAT);
-  long av=avma,av2,lx,vx,i,j,k,ex,nbfac,zval;
+  GEN T,V,W,P,e,cf;
+  long i,k,dW,n,val;
 
-  if (typ(x)!=t_POL) err(notpoler,"factpol");
-  if (!signe(x)) err(zeropoler,"factpol");
+  val = polvaluation(f, &f);
+  n = 1 + degpol(f); if (val) n++;
+  e = cgetg(n,t_VECSMALL);
+  P = cgetg(n,t_COL);
 
-  ex = nbfac = 0;
-  p1 = x+2; while (gcmp0((GEN)*p1)) p1++;
-  zval = p1 - (x + 2);
-  lx = lgef(x) - zval;
-  vx = varn(x);
-  if (zval)
+  cf = content(f); if (gsigne(leading_term(f)) < 0) cf = gneg_i(cf);
+  if (!gcmp1(cf)) f = gdiv(f,cf);
+
+  T = modulargcd(derivpol(f), f);
+  V = gdeuc(f,T);
+  for (k=i=1;; k++)
   {
-    x = cgetg(lx, t_POL); p1 -= 2;
-    for (i=2; i<lx; i++) x[i] = p1[i];
-    x[1] = evalsigne(1)|evalvarn(vx)|evallgef(lx);
-    nbfac++;
+    W = modulargcd(T,V); T = gdeuc(T,W); dW = degpol(W);
+    /* W = prod P^e, e > k; V = prod P^e, e >= k */
+    if (dW != degpol(V)) { P[i] = ldeuc(V,W); e[i] = k; i++; }
+    if (dW <= 0) break;
+    V = W;
   }
-  /* now x(0) != 0 */
-  if (lx==3) { fa = NULL;/* for lint */ goto END; }
-  p1 = cgetg(1,t_VEC); fa=cgetg(lx,t_VEC);
-  for (i=1; i<lx; i++) fa[i] = (long)p1;
-  d=content(x); if (gsigne(leading_term(x)) < 0) d = gneg_i(d);
-  if (!gcmp1(d)) x=gdiv(x,d);
-  if (lx==4) { nbfac++; ex++; fa[1] = (long)concatsp(p1,x); goto END; }
+  if (val) { P[i] = lpolx[varn(f)]; e[i] = val; i++;}
+  setlg(P,i);
+  setlg(e,i); *ex=e; return P;
+}
 
-  w=derivpol(x); t=modulargcd(x,w);
-  if (!gcmp1(t)) { x=gdeuc(x,t); w=gdeuc(w,t); }
-  k=1;
-  while (k)
+GEN
+fact_from_DDF(GEN fa, GEN e, long n)
+{
+  GEN v,w, y = cgetg(3, t_MAT);
+  long i,j,k, l = lg(fa);
+
+  v = cgetg(n+1,t_COL); y[1] = (long)v;
+  w = cgetg(n+1,t_COL); y[2] = (long)w;
+  for (k=i=1; i<l; i++)
   {
-    ex++; w=gadd(w, gneg_i(derivpol(x))); k=signe(w);
-    if (k) { t=modulargcd(x,w); x=gdeuc(x,t); w=gdeuc(w,t); } else t=x;
-    if (degpol(t) > 0)
+    GEN L = (GEN)fa[i], ex = stoi(e[i]);
+    long J = lg(L);
+    for (j=1; j<J; j++,k++)
     {
-      fa[ex] = (long)squff2(t,hint);
-      nbfac += lg(fa[ex])-1;
+      v[k] = lcopy((GEN)L[j]);
+      w[k] = (long)ex;
     }
   }
-END: av2=avma;
-  v=cgetg(nbfac+1,t_COL); y[1]=(long)v;
-  w=cgetg(nbfac+1,t_COL); y[2]=(long)w;
-  if (zval) { v[1]=lpolx[vx]; w[1]=lstoi(zval); k=1; } else k=0;
-  for (i=1; i<=ex; i++)
-    for (j=1; j<lg(fa[i]); j++)
-    {
-      k++; v[k]=lcopy(gmael(fa,i,j)); w[k]=lstoi(i);
-    }
-  gerepilemanyvec(av,av2,y+1,2);
-  return sort_factor(y, cmpii);
+  return y;
 }
 
+/* Factor x in Z[t]. Assume all factors have degree divisible by hint */
+GEN
+factpol(GEN x, long hint)
+{
+  gpmem_t av = avma;
+  GEN fa,ex,y;
+  long n,i,l;
+
+  if (typ(x)!=t_POL) err(notpoler,"factpol");
+  if (!signe(x)) err(zeropoler,"factpol");
+
+  fa = ZX_squff(x, &ex);
+  l = lg(fa); n = 0;
+  for (i=1; i<l; i++)
+  {
+    fa[i] = (long)DDF2((GEN)fa[i], hint);
+    n += lg(fa[i])-1;
+  }
+  y = fact_from_DDF(fa,ex,n);
+  return gerepileupto(av, sort_factor(y, cmpii));
+}
+
 /***********************************************************************/
 /**                                                                   **/
 /**                          FACTORIZATION                            **/
@@ -1630,12 +1707,7 @@ poltype(GEN x, GEN *ptp, GEN *ptpol, long *ptpa)
 	assign_or_fail((GEN)p1[1],p);
         t[3]=1; break;
       case t_COMPLEX:
-        if (!pcx)
-        {
-          pcx = cgetg(5,t_POL); /* x^2 + 1 */
-          pcx[1] = evalsigne(1)|evalvarn(0)|m_evallgef(5),
-          pcx[2]=pcx[4]=un; pcx[3]=zero;
-        }
+        if (!pcx) pcx = coefs_to_pol(3, gun,gzero,gun); /* x^2 + 1 */
 	for (j=1; j<=2; j++)
 	{
 	  p2 = (GEN)p1[j];
@@ -1756,9 +1828,8 @@ factor0(GEN x,long flag)
   return NULL; /* not reached */
 }
 
-/* assume f and g coprime integer factorizations */
 GEN
-merge_factor_i(GEN f, GEN g)
+concat_factor(GEN f, GEN g)
 {
   GEN h;
   if (lg(f) == 1) return g;
@@ -1766,13 +1837,33 @@ merge_factor_i(GEN f, GEN g)
   h = cgetg(3,t_MAT);
   h[1] = (long)concatsp((GEN)f[1], (GEN)g[1]);
   h[2] = (long)concatsp((GEN)f[2], (GEN)g[2]);
-  return sort_factor_gen(h, cmpii);
+  return h;
 }
 
+/* assume f and g coprime integer factorizations */
 GEN
+merge_factor_i(GEN f, GEN g)
+{
+  if (lg(f) == 1) return g;
+  if (lg(g) == 1) return f;
+  return sort_factor_gen(concat_factor(f,g), cmpii);
+}
+
+GEN
+deg1_from_roots(GEN L, long v)
+{
+  long i, l = lg(L);
+  GEN z = cgetg(l,t_COL);
+  for (i=1; i<l; i++)
+    z[i] = (long)deg1pol_i(gun, gneg((GEN)z[i]), v);
+  return z;
+}
+
+GEN
 factor(GEN x)
 {
-  long tx=typ(x),lx,av,tetpil,i,j,pa,v,r1;
+  long tx=typ(x), lx, i, j, pa, v, r1;
+  gpmem_t av, tetpil;
   GEN  y,p,p1,p2,p3,p4,p5,pol;
 
   if (is_matvec_t(tx))
@@ -1794,7 +1885,7 @@ factor(GEN x)
     case t_INT: return decomp(x);
 
     case t_FRACN:
-      x = gred(x); /* fall through */
+      return gerepileupto(av, factor(gred(x)));
     case t_FRAC:
       p1 = decomp((GEN)x[1]);
       p2 = decomp((GEN)x[2]); p2[2] = (long)gneg_i((GEN)p2[2]);
@@ -1810,9 +1901,7 @@ factor(GEN x)
 
 	case t_COMPLEX: y=cgetg(3,t_MAT); lx=lgef(x)-2; v=varn(x);
 	  av = avma; p1 = roots(x,pa); tetpil = avma;
-          p2=cgetg(lx,t_COL);
-	  for (i=1; i<lx; i++)
-            p2[i] = (long)deg1pol_i(gun, gneg((GEN)p1[i]), v);
+          p2 = deg1_from_roots(p1, v);
 	  y[1]=lpile(av,tetpil,p2);
 	  p3=cgetg(lx,t_COL); for (i=1; i<lx; i++) p3[i] = un;
           y[2]=(long)p3; return y;
@@ -1857,7 +1946,7 @@ factor(GEN x)
                 p2[2] = (long)deg1pol_i((GEN)p1[2], (GEN)p1[1], v);
             }
           }
-          killv = (avma != (ulong)pol);
+          killv = (avma != (gpmem_t)pol);
           if (killv) setvarn(pol, fetch_var());
           switch (typ2(tx))
           {
@@ -1869,7 +1958,7 @@ factor(GEN x)
           switch (typ1(tx))
           {
             case t_POLMOD:
-              if (killv) delete_var();
+              if (killv) (void)delete_var();
               return gerepileupto(av,p1);
             case t_COMPLEX: p5 = gi; break;
             case t_QUAD: p5=cgetg(4,t_QUAD);
@@ -1888,7 +1977,7 @@ factor(GEN x)
               if(typ(p4)==t_POLMOD) p3[j]=lsubst((GEN)p4[2],v,p5);
             }
           }
-          if (killv) delete_var();
+          if (killv) (void)delete_var();
           tetpil=avma; y=cgetg(3,t_MAT);
           y[1]=lcopy(p2);y[2]=lcopy((GEN)p1[2]);
           return gerepile(av,tetpil,y);
@@ -1896,7 +1985,7 @@ factor(GEN x)
       }
 
     case t_RFRACN:
-      x=gred_rfrac(x); /* fall through */
+      return gerepileupto(av, factor(gred_rfrac(x)));
     case t_RFRAC:
       p1=factor((GEN)x[1]);
       p2=factor((GEN)x[2]); p3=gneg_i((GEN)p2[2]);
@@ -1913,7 +2002,36 @@ factor(GEN x)
 #undef typs
 #undef tsh
 
+/* assume n != 0, t_INT. Compute x^|n| using left-right binary powering */
 GEN
+leftright_pow(GEN x, GEN n, void *data,
+             GEN (*sqr)(void*,GEN), GEN (*mul)(void*,GEN,GEN))
+{
+  GEN nd = n+2, y = x;
+  long i, m = *nd, j = 1+bfffo((ulong)m);
+  gpmem_t av = avma, lim = stack_lim(av, 1);
+
+  /* normalize, i.e set highest bit to 1 (we know m != 0) */
+  m<<=j; j = BITS_IN_LONG-j;
+  /* first bit is now implicit */
+  for (i=lgefint(n)-2;;)
+  {
+    for (; j; m<<=1,j--)
+    {
+      y = sqr(data,y);
+      if (m < 0) y = mul(data,y,x); /* first bit set: multiply by base */
+      if (low_stack(lim, stack_lim(av,1)))
+      {
+        if (DEBUGMEM>1) err(warnmem,"leftright_pow");
+        y = gerepilecopy(av, y);
+      }
+    }
+    if (--i == 0) return avma==av? gcopy(y): y;
+    m = *++nd; j = BITS_IN_LONG;
+  }
+}
+
+GEN
 divide_conquer_prod(GEN x, GEN (*mul)(GEN,GEN))
 {
   long i,k,lx = lg(x);
@@ -1936,67 +2054,96 @@ divide_conquer_prod(GEN x, GEN (*mul)(GEN,GEN))
 static GEN static_nf;
 
 static GEN
-myidealmulred(GEN x, GEN y) { return idealmulred(static_nf, x, y, 0); }
-
+idmulred(GEN x, GEN y) { return idealmulred(static_nf, x, y, 0); }
 static GEN
-myidealpowred(GEN x, GEN n) { return idealpowred(static_nf, x, n, 0); }
-
+idpowred(GEN x, GEN n) { return idealpowred(static_nf, x, n, 0); }
 static GEN
-myidealmul(GEN x, GEN y) { return idealmul(static_nf, x, y); }
-
+idmul(GEN x, GEN y) { return idealmul(static_nf, x, y); }
 static GEN
-myidealpow(GEN x, GEN n) { return idealpow(static_nf, x, n); }
+idpow(GEN x, GEN n) { return idealpow(static_nf, x, n); }
+static GEN
+eltmul(GEN x, GEN y) { return element_mul(static_nf, x, y); }
+static GEN
+eltpow(GEN x, GEN n) { return element_pow(static_nf, x, n); }
 
 GEN
-factorback_i(GEN fa, GEN nf, int red)
+_factorback(GEN fa, GEN e, GEN (*_mul)(GEN,GEN), GEN (*_pow)(GEN,GEN))
 {
-  long av=avma,k,l,lx,t=typ(fa);
-  GEN ex, x;
-  GEN (*_mul)(GEN,GEN);
-  GEN (*_pow)(GEN,GEN);
-  if (nf)
+  gpmem_t av = avma;
+  long k,l,lx,t = typ(fa);
+  GEN p,x;
+
+  if (e) /* supplied vector of exponents */
+    p = fa;
+  else /* genuine factorization */
   {
-    static_nf = nf;
-    if (red)
+    if (t != t_MAT || lg(fa) != 3)
     {
-      _mul = &myidealmulred;
-      _pow = &myidealpowred;
+      if (!is_vec_t(t)) err(talker,"not a factorisation in factorback");
+      return gerepileupto(av, divide_conquer_prod(fa, _mul));
     }
-    else
-    {
-      _mul = &myidealmul;
-      _pow = &myidealpow;
-    }
+    p = (GEN)fa[1];
+    e = (GEN)fa[2];
   }
-  else
+  lx = lg(p);
+  t = t_INT; /* dummy */
+  /* check whether e is an integral vector of correct length */
+  if (is_vec_t(typ(e)) && lx == lg(e))
   {
-    _mul = &gmul;
-    _pow = &powgi;
+    for (k=1; k<lx; k++)
+      if (typ(e[k]) != t_INT) break;
+    if (k == lx) t = t_MAT;
   }
-  if ( t == t_VEC || t == t_COL)
-    return gerepileupto(av, divide_conquer_prod(fa, _mul));
-  if (t!=t_MAT || lg(fa)!=3)
-    err(talker,"not a factorisation in factorback");
-  ex=(GEN)fa[2]; fa=(GEN)fa[1];
-  lx = lg(fa); if (lx == 1) return gun;
+  if (t != t_MAT) err(talker,"not a factorisation in factorback");
+  if (lx == 1) return gun;
   x = cgetg(lx,t_VEC);
   for (l=1,k=1; k<lx; k++)
-    if (signe(ex[k]))
-      x[l++] = (long)_pow((GEN)fa[k],(GEN)ex[k]);
+    if (signe(e[k]))
+      x[l++] = (long)_pow((GEN)p[k],(GEN)e[k]);
   setlg(x,l);
   return gerepileupto(av, divide_conquer_prod(x, _mul));
 }
 
 GEN
+factorback_i(GEN fa, GEN e, GEN nf, int red)
+{
+  if (!nf && e && lg(e) > 1 && typ(e[1]) != t_INT) { nf = e; e = NULL; }
+  if (!nf) return _factorback(fa, e, &gmul, &powgi);
+
+  static_nf = checknf(nf);
+  if (red)
+    return _factorback(fa, e, &idmulred, &idpowred);
+  else
+    return _factorback(fa, e, &idmul, &idpow);
+}
+
+GEN
+factorbackelt(GEN fa, GEN e, GEN nf)
+{
+  if (!nf && e && lg(e) > 1 && typ(e[1]) != t_INT) { nf = e; e = NULL; }
+  if (!nf) err(talker, "missing nf in factorbackelt");
+
+  static_nf = checknf(nf);
+  return _factorback(fa, e, &eltmul, &eltpow);
+}
+
+GEN
+factorback0(GEN fa, GEN e, GEN nf)
+{
+  return factorback_i(fa,e,nf,0);
+}
+
+GEN
 factorback(GEN fa, GEN nf)
 {
-  return factorback_i(fa,nf,0);
+  return factorback_i(fa,nf,NULL,0);
 }
 
 GEN
 gisirreducible(GEN x)
 {
-  long av=avma, tx = typ(x),l,i;
+  long tx = typ(x), l, i;
+  gpmem_t av=avma;
   GEN y;
 
   if (is_matvec_t(tx))
@@ -2035,7 +2182,7 @@ gcd0(GEN x, GEN y, long flag)
 static GEN
 triv_cont_gcd(GEN x, GEN y)
 {
-  long av = avma, tetpil;
+  gpmem_t av = avma, tetpil;
   GEN p1 = (typ(x)==t_COMPLEX)? ggcd((GEN)x[1],(GEN)x[2])
                               : ggcd((GEN)x[2],(GEN)x[3]);
   tetpil=avma; return gerepile(av,tetpil,ggcd(p1,y));
@@ -2044,7 +2191,7 @@ triv_cont_gcd(GEN x, GEN y)
 static GEN
 cont_gcd(GEN x, GEN y)
 {
-  long av = avma, tetpil;
+  gpmem_t av = avma, tetpil;
   GEN p1 = content(x);
 
   tetpil=avma; return gerepile(av,tetpil,ggcd(p1,y));
@@ -2056,19 +2203,105 @@ padic_gcd(GEN x, GEN y)
 {
   long v = ggval(x,(GEN)y[2]), w = valp(y);
   if (w < v) v = w;
-  return gpuigs((GEN)y[2], v);
+  return gpowgs((GEN)y[2], v);
 }
 
+/* x,y in Z[i], at least one of which is t_COMPLEX */
+static GEN
+gaussian_gcd(GEN x, GEN y)
+{
+  gpmem_t av = avma;
+  GEN dx = denom(x);
+  GEN dy = denom(y);
+  GEN den = mpppcm(dx, dy);
+
+  x = gmul(x, dx);
+  y = gmul(y, dy);
+  while (!gcmp0(y))
+  {
+    GEN z = gdiv(x,y);
+    GEN r0 = greal(z), r = gfloor(r0);
+    GEN i0 = gimag(z), i = gfloor(i0);
+    if (gcmp(gsub(r0,r), ghalf) > 0) r = addis(r,1);
+    if (gcmp(gsub(i0,i), ghalf) > 0) i = addis(i,1);
+    if (gcmp0(i)) z = r;
+    else
+    {
+      z = cgetg(3, t_COMPLEX);
+      z[1] = (long)r;
+      z[2] = (long)i;
+    }
+    z = gsub(x, gmul(z,y));
+    x = y; y = z;
+  }
+  if (signe(greal(x)) < 0) x = gneg(x);
+  if (signe(gimag(x)) < 0) x = gmul(x, gi);
+  if (typ(x) == t_COMPLEX)
+  {
+    if      (gcmp0((GEN)x[2])) x = (GEN)x[1];
+    else if (gcmp0((GEN)x[1])) x = (GEN)x[2];
+  }
+  return gerepileupto(av, gdiv(x, den));
+}
+
 #define fix_frac(z) if (signe(z[2])<0)\
 {\
   setsigne(z[1],-signe(z[1]));\
   setsigne(z[2],1);\
 }
 
+int
+isrational(GEN x)
+{
+  long i, t = typ(x);
+  if (t != t_POL) return is_rational_t(t);
+  for (i=lgef(x)-1; i>1; i--)
+  {
+    t = typ(x[i]);
+    if (!is_rational_t(t)) return 0;
+  }
+  return 1;
+}
+
+static int
+cx_isrational(GEN x)
+{
+  return (isrational((GEN)x[1]) && isrational((GEN)x[2]));
+}
+
+/* y == 0 */
+static GEN
+zero_gcd(GEN x, GEN y)
+{
+  if (typ(y) == t_PADIC)
+  {
+    GEN p = (GEN)y[2];
+    long v = ggval(x,p), w = valp(y);
+    if (w < v) return padiczero(p, w);
+    if (gcmp0(x)) return padiczero(p, v);
+    return gpowgs(p, v);
+  }
+  switch(typ(x))
+  {
+    case t_INT: case t_FRAC: case t_FRACN:
+      return gabs(x,0);
+
+    case t_COMPLEX:
+      if (cx_isrational(x)) return gaussian_gcd(x, gzero);
+      /* fall through */
+    case t_REAL:
+      return gun;
+
+    default:
+      return gcopy(x);
+  }
+}
+
 GEN
 ggcd(GEN x, GEN y)
 {
-  long l,av,tetpil,i,vx,vy, tx = typ(x), ty = typ(y);
+  long l, i, vx, vy, tx = typ(x), ty = typ(y);
+  gpmem_t av, tetpil;
   GEN p1,z;
 
   if (tx>ty) { swap(x,y); lswap(tx,ty); }
@@ -2078,9 +2311,9 @@ ggcd(GEN x, GEN y)
     for (i=1; i<l; i++) z[i]=lgcd(x,(GEN)y[i]);
     return z;
   }
-  if (is_noncalc_t(tx) || is_noncalc_t(ty)) err(operf,"g",tx,ty);
-  if (gcmp0(x)) return gcopy(y);
-  if (gcmp0(y)) return gcopy(x);
+  if (is_noncalc_t(tx) || is_noncalc_t(ty)) err(operf,"g",x,y);
+  if (gcmp0(x)) return zero_gcd(y, x);
+  if (gcmp0(y)) return zero_gcd(x, y);
   if (is_const_t(tx))
   {
     if (ty == t_FRACN)
@@ -2100,7 +2333,7 @@ ggcd(GEN x, GEN y)
 
       case t_INTMOD: z=cgetg(3,t_INTMOD);
         if (egalii((GEN)x[1],(GEN)y[1]))
-          { copyifstack(x[1],z[1]); }
+          copyifstack(x[1],z[1]);
         else
           z[1] = lmppgcd((GEN)x[1],(GEN)y[1]);
         if (gcmp1((GEN)z[1])) z[2] = zero;
@@ -2118,28 +2351,13 @@ ggcd(GEN x, GEN y)
         return z;
 
       case t_COMPLEX:
-        av=avma; p1=gdiv(x,y);
-        if (gcmp0((GEN)p1[2]))
-        {
-          p1 = (GEN)p1[1];
-          switch(typ(p1))
-          {
-            case t_INT:
-              avma=av; return gcopy(y);
-            case t_FRAC: case t_FRACN:
-              tetpil=avma; return gerepile(av,tetpil,gdiv(y,(GEN)p1[2]));
-            default: avma=av; return gun;
-          }
-        }
-        if (typ(p1[1])==t_INT && typ(p1[2])==t_INT) {avma=av; return gcopy(y);}
-        p1 = ginv(p1); avma=av;
-        if (typ(p1[1])==t_INT && typ(p1[2])==t_INT) return gcopy(x);
+        if (cx_isrational(x) && cx_isrational(y)) return gaussian_gcd(x,y);
         return triv_cont_gcd(y,x);
 
       case t_PADIC:
         if (!egalii((GEN)x[2],(GEN)y[2])) return gun;
         vx = valp(x);
-        vy = valp(y); return gpuigs((GEN)y[2], min(vy,vx));
+        vy = valp(y); return gpowgs((GEN)y[2], min(vy,vx));
 
       case t_QUAD:
         av=avma; p1=gdiv(x,y);
@@ -2171,7 +2389,9 @@ ggcd(GEN x, GEN y)
             z[1] = lmppgcd(x,(GEN)y[1]);
             z[2] = licopy((GEN)y[2]); return z;
 
-          case t_COMPLEX: case t_QUAD:
+          case t_COMPLEX:
+            if (cx_isrational(y)) return gaussian_gcd(x,y);
+          case t_QUAD:
             return triv_cont_gcd(y,x);
 
           case t_PADIC:
@@ -2185,7 +2405,7 @@ ggcd(GEN x, GEN y)
         {
           case t_FRAC:
             av = avma; p1=mppgcd((GEN)x[1],(GEN)y[2]); avma = av;
-            if (!is_pm1(p1)) err(operi,"g",tx,ty);
+            if (!is_pm1(p1)) err(operi,"g",x,y);
             return ggcd((GEN)y[1], x);
 
           case t_COMPLEX: case t_QUAD:
@@ -2198,7 +2418,9 @@ ggcd(GEN x, GEN y)
       case t_FRAC:
         switch(ty)
         {
-          case t_COMPLEX: case t_QUAD:
+          case t_COMPLEX:
+            if (cx_isrational(y)) return gaussian_gcd(x,y);
+          case t_QUAD:
             return triv_cont_gcd(y,x);
 
           case t_PADIC:
@@ -2226,7 +2448,7 @@ ggcd(GEN x, GEN y)
       {
 	case t_POLMOD: z=cgetg(3,t_POLMOD);
           if (gegal((GEN)x[1],(GEN)y[1]))
-	    { copyifstack(x[1],z[1]); }
+	    copyifstack(x[1],z[1]);
           else
             z[1] = lgcd((GEN)x[1],(GEN)y[1]);
 	  if (lgef(z[1])<=3) z[2]=zero;
@@ -2251,7 +2473,7 @@ ggcd(GEN x, GEN y)
 
 	case t_RFRAC:
 	  av = avma; p1=ggcd((GEN)x[1],(GEN)y[2]); avma = av;
-          if (!gcmp1(p1)) err(operi,"g",tx,ty);
+          if (!gcmp1(p1)) err(operi,"g",x,y);
 	  return ggcd((GEN)y[1],x);
 
 	case t_RFRACN:
@@ -2267,7 +2489,7 @@ ggcd(GEN x, GEN y)
 	  return srgcd(x,y);
 
 	case t_SER:
-	  return gpuigs(polx[vx],min(valp(y),gval(x,vx)));
+	  return gpowgs(polx[vx],min(valp(y),gval(x,vx)));
 
 	case t_RFRAC: case t_RFRACN: av=avma; z=cgetg(3,ty);
           z[1]=lgcd(x,(GEN)y[1]);
@@ -2279,21 +2501,21 @@ ggcd(GEN x, GEN y)
       switch(ty)
       {
 	case t_SER:
-	  return gpuigs(polx[vx],min(valp(x),valp(y)));
+	  return gpowgs(polx[vx],min(valp(x),valp(y)));
 
 	case t_RFRAC: case t_RFRACN:
-	  return gpuigs(polx[vx],min(valp(x),gval(y,vx)));
+	  return gpowgs(polx[vx],min(valp(x),gval(y,vx)));
       }
       break;
 
     case t_RFRAC: case t_RFRACN: z=cgetg(3,ty);
-      if (!is_rfrac_t(ty)) err(operf,"g",tx,ty);
+      if (!is_rfrac_t(ty)) err(operf,"g",x,y);
       p1 = gdiv((GEN)y[2], ggcd((GEN)x[2], (GEN)y[2]));
       tetpil = avma;
-      z[2] = lpile((long)z,tetpil,gmul(p1, (GEN)x[2]));
+      z[2] = lpile((gpmem_t)z,tetpil,gmul(p1, (GEN)x[2]));
       z[1] = lgcd((GEN)x[1], (GEN)y[1]); return z;
   }
-  err(operf,"g",tx,ty);
+  err(operf,"g",x,y);
   return NULL; /* not reached */
 }
 
@@ -2305,7 +2527,8 @@ GEN glcm0(GEN x, GEN y)
 GEN
 glcm(GEN x, GEN y)
 {
-  long av,tx,ty,i,l;
+  long tx, ty, i, l;
+  gpmem_t av;
   GEN p1,p2,z;
 
   ty = typ(y);
@@ -2340,22 +2563,36 @@ glcm(GEN x, GEN y)
   return gerepileupto(av,p2);
 }
 
+/* x + r ~ x ? Assume x,r are t_POL, deg(r) <= deg(x) */
+static int
+pol_approx0(GEN r, GEN x, int exact)
+{
+  long i, lx,lr;
+  if (exact) return gcmp0(r);
+  lx = lgef(x);
+  lr = lgef(r); if (lr < lx) lx = lr;
+  for (i=2; i<lx; i++)
+    if (!approx_0((GEN)r[i], (GEN)x[i])) return 0;
+  return 1;
+}
+
 static GEN
 polgcdnun(GEN x, GEN y)
 {
-  long av1, av = avma, lim = stack_lim(av,1);
-  GEN p1, yorig = y;
+  gpmem_t av1, av = avma, lim = stack_lim(av, 1);
+  GEN r, yorig = y;
+  int exact = !(isinexactreal(x) || isinexactreal(y));
 
   for(;;)
   {
-    av1 = avma; p1 = gres(x,y);
-    if (gcmp0(p1))
+    av1 = avma; r = gres(x,y);
+    if (pol_approx0(r, x, exact))
     {
       avma = av1;
-      if (lgef(y) == 3 && isinexactreal(y)) { avma = av; return gun; }
+      if (lgef(y) == 3 && !exact) { avma = av; return gun; }
       return (y==yorig)? gcopy(y): gerepileupto(av,y);
     }
-    x = y; y = p1;
+    x = y; y = r;
     if (low_stack(lim,stack_lim(av,1)))
     {
       GEN *gptr[2]; x = gcopy(x); gptr[0]=&x; gptr[1]=&y;
@@ -2365,39 +2602,32 @@ polgcdnun(GEN x, GEN y)
   }
 }
 
+static int issimplefield(GEN x);
+static int
+issimplepol(GEN x)
+{
+  long i, lx = lgef(x);
+  for (i=2; i<lx; i++)
+    if (issimplefield((GEN)x[i])) return 1;
+  return 0;
+}
+
 /* return 1 if coeff explosion is not possible */
 static int
 issimplefield(GEN x)
 {
-  long lx,i;
   switch(typ(x))
   {
     case t_REAL: case t_INTMOD: case t_PADIC: case t_SER:
       return 1;
-    case t_POL:
-      lx=lgef(x);
-      for (i=2; i<lx; i++)
-	if (issimplefield((GEN)x[i])) return 1;
-      return 0;
-    case t_COMPLEX: case t_POLMOD:
+    case t_COMPLEX:
       return issimplefield((GEN)x[1]) || issimplefield((GEN)x[2]);
+    case t_POLMOD: 
+      return issimplepol((GEN)x[1]) || issimplepol((GEN)x[2]);
   }
   return 0;
 }
 
-int
-isrational(GEN x)
-{
-  long i;
-  for (i=lgef(x)-1; i>1; i--)
-    switch(typ(x[i]))
-    {
-      case t_INT: case t_FRAC: break;
-      default: return 0;
-    }
-  return 1;
-}
-
 static int
 can_use_modular_gcd(GEN x, GEN *mod, long *v)
 {
@@ -2423,7 +2653,7 @@ can_use_modular_gcd(GEN x, GEN *mod, long *v)
         break;
       case t_POL:
         if (!isrational(p1)) return 0;
-        if (*v >= 0) 
+        if (*v >= 0)
         {
           if (*v != varn(p1)) return 0;
         } else *v = varn(p1);
@@ -2456,11 +2686,12 @@ isinexactfield(GEN x)
 static GEN
 gcdmonome(GEN x, GEN y)
 {
-  long tetpil,av=avma, dx=degpol(x), v=varn(x), e=gval(y,v);
+  long dx=degpol(x), v=varn(x), e=gval(y, v);
+  gpmem_t tetpil, av=avma;
   GEN p1,p2;
 
   if (dx < e) e = dx;
-  p1=ggcd(leading_term(x),content(y)); p2=gpuigs(polx[v],e);
+  p1=ggcd(leading_term(x),content(y)); p2=gpowgs(polx[v],e);
   tetpil=avma; return gerepile(av,tetpil,gmul(p1,p2));
 }
 
@@ -2473,7 +2704,7 @@ gcdmonome(GEN x, GEN y)
 static GEN
 polinvinexact(GEN x, GEN y)
 {
-  ulong av = avma;
+  gpmem_t av = avma;
   long i,dx=degpol(x),dy=degpol(y),lz=dx+dy;
   GEN v,z;
 
@@ -2490,7 +2721,8 @@ polinvinexact(GEN x, GEN y)
 static GEN
 polinvmod(GEN x, GEN y)
 {
-  long av,av1,tx,vx=varn(x),vy=varn(y);
+  long tx, vx=varn(x), vy=varn(y);
+  gpmem_t av, av1;
   GEN u,v,d,p1;
 
   while (vx!=vy)
@@ -2560,11 +2792,20 @@ GEN
 content(GEN x)
 {
   long lx,i,tx=typ(x);
-  ulong av,tetpil;
+  gpmem_t av, tetpil;
   GEN p1,p2;
 
   if (is_scalar_t(tx))
-    return tx==t_POLMOD? content((GEN)x[2]): gcopy(x);
+  {
+    switch (tx)
+    {
+      case t_INT: return absi(x);
+      case t_FRAC:
+      case t_FRACN: return gabs(x,0);
+      case t_POLMOD: return content((GEN)x[2]);
+    }
+    return gcopy(x);
+  }
   av = avma;
   switch(tx)
   {
@@ -2614,13 +2855,168 @@ content(GEN x)
 }
 
 GEN
-primitive_part(GEN x, GEN *c)
+primitive_part(GEN x, GEN *ptc)
 {
-  *c = content(x);
-  if (gcmp1(*c)) *c = NULL; else x = gdiv(x,*c);
+  gpmem_t av = avma;
+  GEN c = content(x);
+  if (gcmp1(c)) { avma = av; c = NULL; }
+  else if (!gcmp0(c)) x = gdiv(x,c);
+  if (ptc) *ptc = c;
   return x;
 }
 
+GEN
+Q_primitive_part(GEN x, GEN *ptc)
+{
+  gpmem_t av = avma;
+  GEN c = content(x);
+  if (gcmp1(c)) { avma = av; c = NULL; }
+  else if (!gcmp0(c)) x = Q_div_to_int(x,c);
+  if (ptc) *ptc = c;
+  return x;
+}
+
+GEN
+primpart(GEN x) { return primitive_part(x, NULL); }
+
+GEN
+Q_primpart(GEN x) { return Q_primitive_part(x, NULL); }
+
+/* NOT MEMORY CLEAN (because of t_FRAC).
+ * As denom(), but over Q. Treats polynomial as elts of Q[x1,...xn], instead
+ * of Q(x2,...,xn)[x1] */
+GEN
+Q_denom(GEN x)
+{
+  long i, l = lgef(x);
+  GEN d, D;
+  gpmem_t av;
+
+  switch(typ(x))
+  {
+    case t_INT: return gun;
+    case t_FRAC: return (GEN)x[2];
+
+    case t_VEC: case t_COL: case t_MAT:
+      l = lg(x); if (l == 1) return gun;
+      av = avma; d = Q_denom((GEN)x[1]);
+      for (i=2; i<l; i++)
+      {
+        D = Q_denom((GEN)x[i]);
+        if (D != gun) d = mpppcm(d, D);
+      }
+      return gerepileuptoint(av, d);
+
+    case t_POL:
+      l = lgef(x); if (l == 2) return gun;
+      av = avma; d = Q_denom((GEN)x[2]);
+      for (i=3; i<l; i++)
+      {
+        D = Q_denom((GEN)x[i]);
+        if (D != gun) d = mpppcm(d, D);
+      }
+      return gerepileuptoint(av, d);
+  }
+  err(typeer,"Q_denom");
+  return NULL; /* not reached */
+}
+
+GEN
+Q_remove_denom(GEN x, GEN *ptd)
+{
+  GEN d = Q_denom(x);
+  if (gcmp1(d)) d = NULL; else x = Q_muli_to_int(x,d);
+  if (ptd) *ptd = d;
+  return x;
+}
+
+/* return y = x * d, assuming x rational, and d,y integral */
+GEN
+Q_muli_to_int(GEN x, GEN d)
+{
+  long i, l, t = typ(x);
+  GEN y, xn, xd;
+  gpmem_t av;
+
+  if (typ(d) != t_INT) err(typeer,"Q_muli_to_int");
+  switch (t)
+  {
+    case t_VEC: case t_COL: case t_MAT:
+      l = lg(x); y = cgetg(l, t);
+      for (i=1; i<l; i++) y[i] = (long)Q_muli_to_int((GEN)x[i], d);
+      return y;
+
+    case t_POL: l = lgef(x);
+      y = cgetg(l, t_POL); y[1] = x[1];
+      for (i=2; i<l; i++) y[i] = (long)Q_muli_to_int((GEN)x[i], d);
+      return y;
+
+    case t_INT:
+      return mulii(x,d);
+
+    case t_FRAC:
+      xn = (GEN)x[1];
+      xd = (GEN)x[2]; av = avma;
+      y = mulii(xn, diviiexact(d, xd));
+      return gerepileuptoint(av, y);
+  }
+  err(typeer,"Q_muli_to_int");
+  return NULL; /* not reached */
+}
+
+/* return x * n/d. x: rational; d,n,result: integral. n = NULL represents 1 */
+GEN
+Q_divmuli_to_int(GEN x, GEN d, GEN n)
+{
+  long i, l, t = typ(x);
+  GEN y, xn, xd;
+  gpmem_t av;
+  
+  switch(t)
+  {
+    case t_VEC: case t_COL: case t_MAT:
+      l = lg(x); y = cgetg(l, t);
+      for (i=1; i<l; i++) y[i] = (long)Q_divmuli_to_int((GEN)x[i], d,n);
+      return y;
+
+    case t_POL: l = lgef(x);
+      y = cgetg(l, t_POL); y[1] = x[1];
+      for (i=2; i<l; i++) y[i] = (long)Q_divmuli_to_int((GEN)x[i], d,n);
+      return y;
+
+    case t_INT:
+      av = avma; y = diviiexact(x,d);
+      if (n) y = gerepileuptoint(av, mulii(y,n));
+      return y;
+
+    case t_FRAC:
+      xn = (GEN)x[1];
+      xd = (GEN)x[2]; av = avma;
+      y = mulii(diviiexact(xn, d), diviiexact(n, xd));
+      return gerepileuptoint(av, y);
+  }
+  err(typeer,"Q_divmuli_to_int");
+  return NULL; /* not reached */
+}
+
+/* return y = x / c, assuming x,c rational, and y integral */
+GEN
+Q_div_to_int(GEN x, GEN c)
+{
+  GEN d, n;
+  switch(typ(c))
+  {
+    case t_INT:
+      return Q_divmuli_to_int(x, c, NULL);
+    case t_FRAC:
+      n = (GEN)c[1];
+      d = (GEN)c[2]; if (gcmp1(n)) return Q_muli_to_int(x,d);
+      return Q_divmuli_to_int(x, n,d);
+  }
+  err(typeer,"Q_div_to_int");
+  return NULL; /* not reached */
+}
+
 /*******************************************************************/
 /*                                                                 */
 /*                           SUBRESULTANT                          */
@@ -2647,7 +3043,7 @@ gdivexact(GEN x, GEN y)
         case t_INTMOD:
         case t_POLMOD: return gmul(x,ginv(y));
         case t_POL:
-          if (varn(x)==varn(y)) return poldivres(x,y, ONLY_DIVIDES_EXACT);
+          if (varn(x)==varn(y)) return poldivres(x,y, ONLY_DIVIDES);
       }
       lx = lgef(x); z = cgetg(lx,t_POL);
       for (i=2; i<lx; i++) z[i]=(long)gdivexact((GEN)x[i],y);
@@ -2669,17 +3065,17 @@ init_resultant(GEN x, GEN y)
   tx = typ(x); ty = typ(y);
   if (is_scalar_t(tx) || is_scalar_t(ty))
   {
-    if (tx==t_POL) return gpuigs(y,degpol(x));
-    if (ty==t_POL) return gpuigs(x,degpol(y));
+    if (tx==t_POL) return gpowgs(y,degpol(x));
+    if (ty==t_POL) return gpowgs(x,degpol(y));
     return gun;
   }
   if (tx!=t_POL || ty!=t_POL) err(typeer,"subresall");
   if (varn(x)==varn(y)) return NULL;
-  return (varn(x)<varn(y))? gpuigs(y,degpol(x)): gpuigs(x,degpol(y));
+  return (varn(x)<varn(y))? gpowgs(y,degpol(x)): gpowgs(x,degpol(y));
 }
 
 /* return coefficients s.t x = x_0 X^n + ... + x_n */
-static GEN
+GEN
 revpol(GEN x)
 {
   long i,n = degpol(x);
@@ -2689,12 +3085,16 @@ revpol(GEN x)
   return y;
 }
 
-/* assume dx >= dy, y non constant */
-GEN
-pseudorem(GEN x, GEN y)
+/* assume dx >= dy, y non constant. */
+static GEN
+pseudorem_i(GEN x, GEN y, GEN mod)
 {
-  long av = avma, av2,lim, vx = varn(x),dx,dy,dz,i,lx, p;
+  long vx = varn(x), dx, dy, dz, i, lx, p;
+  gpmem_t av = avma, av2, lim;
+  GEN (*red)(GEN,GEN) = NULL;
 
+  if (mod) red = (typ(mod) == t_INT)? &FpX_red: &gmod;
+
   if (!signe(y)) err(talker,"euclidean division by zero (pseudorem)");
   (void)new_chunk(2);
   dx=degpol(x); x = revpol(x);
@@ -2704,9 +3104,15 @@ pseudorem(GEN x, GEN y)
   {
     x[0] = lneg((GEN)x[0]); p--;
     for (i=1; i<=dy; i++)
+    {
       x[i] = ladd(gmul((GEN)y[0], (GEN)x[i]), gmul((GEN)x[0],(GEN)y[i]));
+      if (red) x[i] = (long)red((GEN)x[i], mod);
+    }
     for (   ; i<=dx; i++)
+    {
       x[i] = lmul((GEN)y[0], (GEN)x[i]);
+      if (red) x[i] = (long)red((GEN)x[i], mod);
+    }
     do { x++; dx--; } while (dx >= 0 && gcmp0((GEN)x[0]));
     if (dx < dy) break;
     if (low_stack(lim,stack_lim(av2,1)))
@@ -2720,20 +3126,39 @@ pseudorem(GEN x, GEN y)
   x[0]=evaltyp(t_POL) | evallg(lx);
   x[1]=evalsigne(1) | evalvarn(vx) | evallgef(lx);
   x = revpol(x) - 2;
-  return gerepileupto(av, gmul(x, gpowgs((GEN)y[0], p)));
+  if (p)
+  { /* multiply by y[0]^p   [beware dummy vars from FpY_FpXY_resultant] */
+    GEN t = (GEN)y[0];
+    if (mod)
+    { /* assume p fairly small */
+      for (i=1; i<p; i++)
+        t = red(gmul(t, (GEN)y[0]), mod);
+    }
+    else
+      t = gpowgs(t, p);
+    for (i=2; i<lx; i++)
+    {
+      x[i] = lmul((GEN)x[i], t);
+      if (red) x[i] = (long)red((GEN)x[i], mod);
+    }
+    if (!red) return gerepileupto(av, x);
+  }
+  return gerepilecopy(av, x);
 }
 
-extern void gerepilemanycoeffs2(long av, GEN x, long n, GEN y, long o);
+GEN
+pseudorem(GEN x, GEN y) { return pseudorem_i(x,y, NULL); }
 
 /* assume dx >= dy, y non constant
  * Compute z,r s.t lc(y)^(dx-dy+1) x = z y + r */
 GEN
 pseudodiv(GEN x, GEN y, GEN *ptr)
 {
-  long av = avma, av2,lim, vx = varn(x),dx,dy,dz,i,iz,lx,lz,p;
+  long vx = varn(x), dx, dy, dz, i, iz, lx, lz, p;
+  gpmem_t av = avma, av2, lim;
   GEN z, r, ypow;
 
-  if (!signe(y)) err(talker,"euclidean division by zero (pseudorem)");
+  if (!signe(y)) err(talker,"euclidean division by zero (pseudodiv)");
   (void)new_chunk(2);
   dx=degpol(x); x = revpol(x);
   dy=degpol(y); y = revpol(y); dz=dx-dy; p = dz+1;
@@ -2751,7 +3176,7 @@ pseudodiv(GEN x, GEN y, GEN *ptr)
       x[i] = ladd(gmul((GEN)y[0], (GEN)x[i]), gmul((GEN)x[0],(GEN)y[i]));
     for (   ; i<=dx; i++)
       x[i] = lmul((GEN)y[0], (GEN)x[i]);
-    x++; dx--; 
+    x++; dx--;
     while (dx >= dy && gcmp0((GEN)x[0])) { x++; dx--; z[iz++] = zero; }
     if (dx < dy) break;
     if (low_stack(lim,stack_lim(av2,1)))
@@ -2761,7 +3186,7 @@ pseudodiv(GEN x, GEN y, GEN *ptr)
     }
   }
   while (dx >= 0 && gcmp0((GEN)x[0])) { x++; dx--; }
-  if (dx < 0) 
+  if (dx < 0)
     x = zeropol(vx);
   else
   {
@@ -2776,20 +3201,16 @@ pseudodiv(GEN x, GEN y, GEN *ptr)
   z[1] = evalsigne(1) | evalvarn(vx) | evallgef(lz);
   z = revpol(z) - 2;
   r = gmul(x, (GEN)ypow[p]);
-  {
-    GEN *gptr[2]; gptr[0] = &z; gptr[1] = &r;
-    gerepilemany(av,gptr,2);
-  }
+  gerepileall(av, 2, &z, &r);
   *ptr = r; return z;
 }
 
 /* Return resultant(u,v). If sol != NULL: set *sol to the last non-zero
- * polynomial in the prs IF the sequence was computed, and gzero otherwise
- */
+ * polynomial in the prs IF the sequence was computed, and gzero otherwise */
 GEN
 subresall(GEN u, GEN v, GEN *sol)
 {
-  ulong av,av2,lim;
+  gpmem_t av, av2, lim;
   long degq,dx,dy,du,dv,dr,signh;
   GEN z,g,h,r,p1,p2,cu,cv;
 
@@ -2797,13 +3218,13 @@ subresall(GEN u, GEN v, GEN *sol)
   if ((r = init_resultant(u,v))) return r;
 
   if (isinexactfield(u) || isinexactfield(v)) return resultant2(u,v);
-  dx=lgef(u); dy=lgef(v); signh=1;
-  if (dx<dy)
+  dx=degpol(u); dy=degpol(v); signh=1;
+  if (dx < dy)
   {
     swap(u,v); lswap(dx,dy);
-    if (both_even(dx, dy)) signh = -signh;
+    if (both_odd(dx, dy)) signh = -signh;
   }
-  if (dy==3) return gpowgs((GEN)v[2],dx-3);
+  if (dy==0) return gpowgs((GEN)v[2],dx);
   av = avma;
   u = primitive_part(u, &cu);
   v = primitive_part(v, &cv);
@@ -2813,10 +3234,10 @@ subresall(GEN u, GEN v, GEN *sol)
     r = pseudorem(u,v); dr = lgef(r);
     if (dr == 2)
     {
-      if (sol) { avma = (long)(r+2); *sol=gerepileupto(av,v); } else avma = av;
+      if (sol) { avma = (gpmem_t)(r+2); *sol=gerepileupto(av,v); } else avma = av;
       return gzero;
     }
-    du = lgef(u); dv = lgef(v); degq = du-dv;
+    du = degpol(u); dv = degpol(v); degq = du-dv;
     u = v; p1 = g; g = leading_term(u);
     switch(degq)
     {
@@ -2827,22 +3248,21 @@ subresall(GEN u, GEN v, GEN *sol)
         p1 = gmul(gpowgs(h,degq),p1);
         h = gdivexact(gpowgs(g,degq), gpowgs(h,degq-1));
     }
-    if (both_even(du,dv)) signh = -signh;
+    if (both_odd(du,dv)) signh = -signh;
     v = gdivexact(r,p1);
     if (dr==3) break;
     if (low_stack(lim,stack_lim(av2,1)))
     {
-      GEN *gptr[4]; gptr[0]=&u; gptr[1]=&v; gptr[2]=&g; gptr[3]=&h;
       if(DEBUGMEM>1) err(warnmem,"subresall, dr = %ld",dr);
-      gerepilemany(av2,gptr,4);
+      gerepileall(av2,4, &u, &v, &g, &h);
     }
   }
   z = (GEN)v[2];
-  if (dv > 4) z = gdivexact(gpowgs(z,dv-3), gpowgs(h,dv-4));
+  if (dv > 1) z = gdivexact(gpowgs(z,dv), gpowgs(h,dv-1));
   if (signh < 0) z = gneg(z); /* z = resultant(ppart(x), ppart(y)) */
   p2 = gun;
-  if (cu) p2 = gmul(p2, gpowgs(cu,dy-3));
-  if (cv) p2 = gmul(p2, gpowgs(cv,dx-3));
+  if (cu) p2 = gmul(p2, gpowgs(cu,dy));
+  if (cv) p2 = gmul(p2, gpowgs(cv,dx));
   z = gmul(z, p2);
 
   if (sol) u = gclone(u);
@@ -2854,17 +3274,16 @@ subresall(GEN u, GEN v, GEN *sol)
 static GEN
 scalar_res(GEN x, GEN y, GEN *U, GEN *V)
 {
-  long dx=lgef(x)-4;
-  *V=gpuigs(y,dx); *U=gzero; return gmul(y,*V);
+  *V=gpowgs(y,degpol(x)-1); *U=gzero; return gmul(y,*V);
 }
 
 /* compute U, V s.t Ux + Vy = resultant(x,y) */
 GEN
 subresext(GEN x, GEN y, GEN *U, GEN *V)
 {
-  ulong av,av2,tetpil,lim;
+  gpmem_t av, av2, tetpil, lim;
   long degq,tx,ty,dx,dy,du,dv,dr,signh;
-  GEN q,z,g,h,r,p1,p2,cu,cv,u,v,um1,uze,vze,xprim,yprim;
+  GEN q,z,g,h,r,p1,p2,cu,cv,u,v,um1,uze,vze,xprim,yprim, *gptr[3];
 
   if (gcmp0(x) || gcmp0(y)) { *U = *V = gzero; return gzero; }
   tx=typ(x); ty=typ(y);
@@ -2878,15 +3297,15 @@ subresext(GEN x, GEN y, GEN *U, GEN *V)
   if (varn(x) != varn(y))
     return (varn(x)<varn(y))? scalar_res(x,y,U,V)
                             : scalar_res(y,x,V,U);
-  dx = lgef(x); dy = lgef(y); signh = 1;
+  dx = degpol(x); dy = degpol(y); signh = 1;
   if (dx < dy)
   {
     pswap(U,V); lswap(dx,dy); swap(x,y);
-    if (both_even(dx, dy)) signh = -signh;
+    if (both_odd(dx, dy)) signh = -signh;
   }
-  if (dy == 3)
+  if (dy == 0)
   {
-    *V = gpowgs((GEN)y[2],dx-4);
+    *V = gpowgs((GEN)y[2],dx-1);
     *U = gzero; return gmul(*V,(GEN)y[2]);
   }
   av = avma;
@@ -2899,9 +3318,9 @@ subresext(GEN x, GEN y, GEN *U, GEN *V)
     q = pseudodiv(u,v, &r); dr = lgef(r);
     if (dr == 2) { *U = *V = gzero; avma = av; return gzero; }
 
-    du = lgef(u); dv = lgef(v); degq = du-dv;
+    du = degpol(u); dv = degpol(v); degq = du-dv;
     /* lead(v)^(degq + 1) * um1 - q * uze */
-    p1 = gsub(gmul(gpowgs((GEN)v[dv-1],degq+1),um1), gmul(q,uze));
+    p1 = gsub(gmul(gpowgs((GEN)v[dv+2],degq+1),um1), gmul(q,uze));
     um1 = uze; uze = p1;
     u = v; p1 = g; g = leading_term(u);
     switch(degq)
@@ -2912,42 +3331,42 @@ subresext(GEN x, GEN y, GEN *U, GEN *V)
         p1 = gmul(gpowgs(h,degq),p1);
         h = gdivexact(gpowgs(g,degq), gpowgs(h,degq-1));
     }
-    if (both_even(du, dv)) signh = -signh;
+    if (both_odd(du, dv)) signh = -signh;
     v  = gdivexact(r,p1);
     uze= gdivexact(uze,p1);
     if (dr==3) break;
     if (low_stack(lim,stack_lim(av2,1)))
     {
-      GEN *gptr[6]; gptr[0]=&u; gptr[1]=&v; gptr[2]=&g; gptr[3]=&h;
-      gptr[4]=&uze; gptr[5]=&um1;
       if(DEBUGMEM>1) err(warnmem,"subresext, dr = %ld",dr);
-      gerepilemany(av2,gptr,6);
+      gerepileall(av2,6, &u,&v,&g,&h,&uze,&um1);
     }
   }
   z = (GEN)v[2];
-  if (dv > 4) 
+  if (dv > 1)
   {
-    /* z = gdivexact(gpowgs(z,dv-3), gpowgs(h,dv-4)); */
-    p1 = gpowgs(gdiv(z,h),dv-4);
+    /* z = gdivexact(gpowgs(z,dv), gpowgs(h,dv-1)); */
+    p1 = gpowgs(gdiv(z,h),dv-1);
     z = gmul(z,p1); uze = gmul(uze,p1);
   }
   if (signh < 0) { z = gneg_i(z); uze = gneg_i(uze); }
   p1 = gadd(z, gneg(gmul(uze,xprim)));
-  if (!poldivis(p1,yprim,&vze)) err(bugparier,"subresext");
+  vze = poldivres(p1, yprim, &r);
+  if (!gcmp0(r)) err(warner,"inexact computation in subresext");
   /* uze ppart(x) + vze ppart(y) = z = resultant(ppart(x), ppart(y)), */
   p2 = gun;
-  if (cu) p2 = gmul(p2, gpowgs(cu,dy-3));
-  if (cv) p2 = gmul(p2, gpowgs(cv,dx-3));
+  if (cu) p2 = gmul(p2, gpowgs(cu,dy));
+  if (cv) p2 = gmul(p2, gpowgs(cv,dx));
   cu = cu? gdiv(p2,cu): p2;
   cv = cv? gdiv(p2,cv): p2;
 
   tetpil = avma;
-  z = gmul(z,p2); uze = gmul(uze,cu); vze = gmul(vze,cv);
-  {
-    GEN *gptr[3]; gptr[0]=&z; gptr[1]=&uze; gptr[2]=&vze;
-    gerepilemanysp(av,tetpil,gptr,3); 
-  }
-  *U = uze; *V = vze; return z;
+  z = gmul(z,p2);
+  uze = gmul(uze,cu);
+  vze = gmul(vze,cv);
+  gptr[0]=&z; gptr[1]=&uze; gptr[2]=&vze;
+  gerepilemanysp(av,tetpil,gptr,3);
+  *U = uze;
+  *V = vze; return z;
 }
 
 static GEN
@@ -2967,7 +3386,7 @@ zero_bezout(GEN y, GEN *U, GEN *V)
 GEN
 bezoutpol(GEN x, GEN y, GEN *U, GEN *V)
 {
-  ulong av,av2,tetpil,lim;
+  gpmem_t av, av2, tetpil, lim;
   long degq,tx,ty,dx,dy,du,dv,dr;
   GEN q,z,g,h,r,p1,cu,cv,u,v,um1,uze,vze,xprim,yprim, *gptr[3];
 
@@ -2984,12 +3403,12 @@ bezoutpol(GEN x, GEN y, GEN *U, GEN *V)
   if (varn(x)!=varn(y))
     return (varn(x)<varn(y))? scalar_bezout(x,y,U,V)
                             : scalar_bezout(y,x,V,U);
-  dx = lgef(x); dy = lgef(y);
+  dx = degpol(x); dy = degpol(y);
   if (dx < dy)
   {
     pswap(U,V); lswap(dx,dy); swap(x,y);
   }
-  if (dy==3) return scalar_bezout(x,y,U,V);
+  if (dy==0) return scalar_bezout(x,y,U,V);
 
   u = primitive_part(x, &cu); xprim = u;
   v = primitive_part(y, &cv); yprim = v;
@@ -3000,8 +3419,8 @@ bezoutpol(GEN x, GEN y, GEN *U, GEN *V)
     q = pseudodiv(u,v, &r); dr = lgef(r);
     if (dr == 2) break;
 
-    du = lgef(u); dv = lgef(v); degq = du-dv;
-    p1 = gsub(gmul(gpowgs((GEN)v[dv-1],degq+1),um1), gmul(q,uze));
+    du = degpol(u); dv = degpol(v); degq = du-dv;
+    p1 = gsub(gmul(gpowgs((GEN)v[dv+2],degq+1),um1), gmul(q,uze));
     um1 = uze; uze = p1;
     u = v; p1 = g; g  = leading_term(u);
     switch(degq)
@@ -3018,24 +3437,25 @@ bezoutpol(GEN x, GEN y, GEN *U, GEN *V)
     if (dr==3) break;
     if (low_stack(lim,stack_lim(av2,1)))
     {
-      GEN *gptr[6]; gptr[0]=&u; gptr[1]=&v; gptr[2]=&g; gptr[3]=&h;
-      gptr[4]=&uze; gptr[5]=&um1;
       if(DEBUGMEM>1) err(warnmem,"bezoutpol, dr = %ld",dr);
-      gerepilemany(av2,gptr,6);
+      gerepileall(av2,6,&u,&v,&g,&h,&uze,&um1);
     }
   }
-  if (!poldivis(gsub(v,gmul(uze,xprim)),yprim, &vze))
-    err(warner,"inexact computation in bezoutpol");
+  p1 = gsub(v,gmul(uze,xprim));
+  vze = poldivres(p1, yprim, &r);
+  if (!gcmp0(r)) err(warner,"inexact computation in bezoutpol");
   if (cu) uze = gdiv(uze,cu);
   if (cv) vze = gdiv(vze,cv);
-  p1 = ginv(content(v)); tetpil = avma;
-
+  p1 = ginv(content(v));
+  
+  tetpil = avma;
   uze = gmul(uze,p1);
   vze = gmul(vze,p1);
   z = gmul(v,p1);
   gptr[0]=&uze; gptr[1]=&vze; gptr[2]=&z;
   gerepilemanysp(av,tetpil,gptr,3);
-  *U = uze; *V = vze; return z;
+  *U = uze;
+  *V = vze; return z;
 }
 
 /*******************************************************************/
@@ -3075,14 +3495,12 @@ Lazard2(GEN F, GEN x, GEN y, long n)
   x = Lazard(x,y,n-1); return gdivexact(gmul(x,F),y);
 }
 
-extern GEN addshiftpol(GEN x, GEN y, long d);
-#define addshift(x,y) addshiftpol((x),(y),1)
-
 static GEN
 nextSousResultant(GEN P, GEN Q, GEN Z, GEN s)
 {
   GEN p0, q0, z0, H, A;
-  long p, q, j, av, lim, v = varn(P);
+  long p, q, j, v = varn(P);
+  gpmem_t av, lim;
 
   z0 = leading_term(Z);
   p = degpol(P); p0 = leading_term(P); P = reductum(P);
@@ -3099,9 +3517,8 @@ nextSousResultant(GEN P, GEN Q, GEN Z, GEN s)
     A = gadd(A,gmul((GEN)P[j+2],H));
     if (low_stack(lim,stack_lim(av,1)))
     {
-      GEN *gptr[2]; gptr[0]=&A; gptr[1]=&H;
       if(DEBUGMEM>1) err(warnmem,"nextSousResultant j = %ld/%ld",j,p);
-      gerepilemany(av,gptr,2);
+      gerepileall(av,2,&A,&H);
     }
   }
   P = normalizepol_i(P, q+2);
@@ -3115,7 +3532,7 @@ nextSousResultant(GEN P, GEN Q, GEN Z, GEN s)
 GEN
 resultantducos(GEN P, GEN Q)
 {
-  ulong av = avma, av2,lim;
+  gpmem_t av = avma, av2, lim;
   long dP,dQ,delta;
   GEN cP,cQ,Z,s;
 
@@ -3134,7 +3551,7 @@ resultantducos(GEN P, GEN Q)
   if (degpol(Q) > 0)
   {
     av2 = avma; lim = stack_lim(av2,1);
-    s = gpuigs(leading_term(Q),delta);
+    s = gpowgs(leading_term(Q),delta);
     Z = Q;
     Q = pseudorem(P, gneg(Q));
     P = Z;
@@ -3142,9 +3559,8 @@ resultantducos(GEN P, GEN Q)
     {
       if (low_stack(lim,stack_lim(av,1)))
       {
-        GEN *gptr[2]; gptr[0]=&P; gptr[1]=&Q;
         if(DEBUGMEM>1) err(warnmem,"resultantducos, degpol Q = %ld",degpol(Q));
-        gerepilemany(av2,gptr,2); s = leading_term(P);
+        gerepileall(av2,2,&P,&Q); s = leading_term(P);
       }
       delta = degpol(P) - degpol(Q);
       Z = Lazard2(Q, leading_term(Q), s, delta);
@@ -3215,7 +3631,7 @@ sylvestermatrix(GEN x, GEN y)
 GEN
 resultant2(GEN x, GEN y)
 {
-  long av;
+  gpmem_t av;
   GEN r;
   if ((r = init_resultant(x,y))) return r;
   av = avma; return gerepileupto(av,det(sylvestermatrix_i(x,y)));
@@ -3251,7 +3667,8 @@ fix_pol(GEN x, long v, long *mx)
 GEN
 polresultant0(GEN x, GEN y, long v, long flag)
 {
-  long av = avma, m = 0;
+  long m = 0;
+  gpmem_t av = avma;
 
   if (v >= 0)
   {
@@ -3274,13 +3691,11 @@ polresultant0(GEN x, GEN y, long v, long flag)
 /*                  GCD USING SUBRESULTANT                         */
 /*                                                                 */
 /*******************************************************************/
-extern GEN nfgcd(GEN P, GEN Q, GEN nf, GEN den);
-extern GEN to_polmod(GEN x, GEN mod);
-
 GEN
 srgcd(GEN x, GEN y)
 {
-  long av,av1,tetpil,tx=typ(x),ty=typ(y),dx,dy,vx,vmod,lim;
+  long tx=typ(x), ty=typ(y), dx, dy, vx, vmod;
+  gpmem_t av, av1, tetpil, lim;
   GEN d,g,h,p1,p2,u,v,mod,cx,cy;
 
   if (!signe(y)) return gcopy(x);
@@ -3309,7 +3724,7 @@ srgcd(GEN x, GEN y)
     if (vmod < 0) return modulargcd(x,y); /* Q[X] */
   }
 
-  if (issimplefield(x) || issimplefield(y)) x = polgcdnun(x,y);
+  if (issimplepol(x) || issimplepol(y)) x = polgcdnun(x,y);
   else
   {
     dx=lgef(x); dy=lgef(y);
@@ -3344,15 +3759,14 @@ srgcd(GEN x, GEN y)
           h = g = leading_term(u);
           break;
         default:
-          v = gdiv(r,gmul(gpuigs(h,degq),g));
+          v = gdiv(r,gmul(gpowgs(h,degq),g));
           g = leading_term(u);
-          h = gdiv(gpuigs(g,degq), gpuigs(h,degq-1));
+          h = gdiv(gpowgs(g,degq), gpowgs(h,degq-1));
       }
       if (low_stack(lim, stack_lim(av1,1)))
       {
-        GEN *gptr[4]; gptr[0]=&u; gptr[1]=&v; gptr[2]=&g; gptr[3]=&h;
         if(DEBUGMEM>1) err(warnmem,"srgcd");
-        gerepilemany(av1,gptr,4);
+        gerepileall(av1,4,&u,&v,&g,&h);
       }
     }
     p1 = content(v); if (!gcmp1(p1)) v = gdiv(v,p1);
@@ -3368,12 +3782,11 @@ srgcd(GEN x, GEN y)
   return gerepileupto(av,x);
 }
 
-extern GEN qf_disc(GEN x, GEN y, GEN z);
-
 GEN
 poldisc0(GEN x, long v)
 {
-  long av,tx=typ(x),i;
+  long tx=typ(x), i;
+  gpmem_t av;
   GEN z,p1,p2;
 
   switch(tx)
@@ -3415,7 +3828,8 @@ discsr(GEN x)
 GEN
 reduceddiscsmith(GEN pol)
 {
-  long av=avma,tetpil,i,j,n;
+  long i, j, n;
+  gpmem_t av=avma, tetpil;
   GEN polp,alpha,p1,m;
 
   if (typ(pol)!=t_POL) err(typeer,"reduceddiscsmith");
@@ -3446,7 +3860,8 @@ reduceddiscsmith(GEN pol)
 long
 sturmpart(GEN x, GEN a, GEN b)
 {
-  long av = avma, lim = stack_lim(av,1), sl,sr,s,t,r1;
+  long sl, sr, s, t, r1;
+  gpmem_t av = avma, lim = stack_lim(av, 1);
   GEN g,h,u,v;
 
   if (typ(x) != t_POL) err(typeer,"sturm");
@@ -3508,15 +3923,14 @@ sturmpart(GEN x, GEN a, GEN b)
       case 1:
         p1 = gmul(h,p1); h = g; break;
       default:
-        p1 = gmul(gpuigs(h,degq),p1);
-        h = gdivexact(gpuigs(g,degq), gpuigs(h,degq-1));
+        p1 = gmul(gpowgs(h,degq),p1);
+        h = gdivexact(gpowgs(g,degq), gpowgs(h,degq-1));
     }
     v = gdivexact(r,p1);
     if (low_stack(lim,stack_lim(av,1)))
     {
-      GEN *gptr[4]; gptr[0]=&u; gptr[1]=&v; gptr[2]=&g; gptr[3]=&h;
       if(DEBUGMEM>1) err(warnmem,"polsturm, dr = %ld",dr);
-      gerepilemany(av,gptr,4);
+      gerepileall(av,4,&u,&v,&g,&h);
     }
   }
 }
@@ -3530,33 +3944,35 @@ sturmpart(GEN x, GEN a, GEN b)
 GEN
 quadpoly0(GEN x, long v)
 {
-  long res,l,tetpil,i,sx, tx = typ(x);
+  long res, i, l, sx, tx = typ(x);
+  gpmem_t av;
   GEN y,p1;
 
   if (is_matvec_t(tx))
   {
-    l=lg(x); y=cgetg(l,tx);
-    for (i=1; i<l; i++) y[i]=(long)quadpoly0((GEN)x[i],v);
+    l = lg(x); y = cgetg(l,tx);
+    for (i=1; i<l; i++) y[i] = (long)quadpoly0((GEN)x[i],v);
     return y;
   }
-  if (tx!=t_INT) err(arither1);
+  if (tx != t_INT) err(arither1);
   if (v < 0) v = 0;
-  sx=signe(x);
+  sx = signe(x);
   if (!sx) err(talker,"zero discriminant in quadpoly");
-  y=cgetg(5,t_POL);
-  y[1]=evalsigne(1) | evalvarn(v) | evallgef(5); y[4]=un;
-  res=mod4(x); if (sx<0 && res) res=4-res;
-  if (res>1) err(funder2,"quadpoly");
+  res = mod4(x); if (sx < 0 && res) res = 4-res;
+  if (res > 1) err(funder2,"quadpoly");
 
-  l=avma; p1=shifti(x,-2); setsigne(p1,-signe(p1));
+  y = cgetg(5,t_POL);
+  y[1] = evalsigne(1) | evalvarn(v) | evallgef(5);
+
+  av = avma; p1 = shifti(x,-2); setsigne(p1,-signe(p1));
   y[2] = (long) p1;
   if (!res) y[3] = zero;
   else
   {
-    if (sx<0) { tetpil=avma; y[2] = lpile(l,tetpil,addsi(1,p1)); }
+    if (sx < 0) y[2] = lpileuptoint(av, addsi(1,p1));
     y[3] = lnegi(gun);
   }
-  return y;
+  y[4] = un; return y;
 }
 
 GEN
@@ -3609,7 +4025,7 @@ newtonpoly(GEN x, GEN p)
 {
   GEN y;
   long n,ind,a,b,c,u1,u2,r1,r2;
-  long *vval, num[] = {evaltyp(t_INT) | m_evallg(3), 0, 0};
+  long *vval, num[] = {evaltyp(t_INT) | _evallg(3), 0, 0};
 
   if (typ(x)!=t_POL) err(typeer,"newtonpoly");
   n=degpol(x); if (n<=0) { y=cgetg(1,t_VEC); return y; }
@@ -3636,17 +4052,12 @@ newtonpoly(GEN x, GEN p)
   free(vval); return y;
 }
 
-extern int cmp_pol(GEN x, GEN y);
-extern GEN ZY_ZXY_resultant(GEN A, GEN B0, long *lambda);
-GEN matratlift(GEN M, GEN mod, GEN amax, GEN bmax, GEN denom);
-GEN nfgcd(GEN P, GEN Q, GEN nf, GEN den);
-
 /* Factor polynomial a on the number field defined by polynomial t */
 GEN
 polfnf(GEN a, GEN t)
 {
-  ulong av = avma;
-  GEN x0,y,p1,p2,u,fa,n,unt,dent,alift;
+  gpmem_t av = avma;
+  GEN x0,y,p1,p2,u,G,fa,n,unt,dent,alift;
   long lx,i,k,e;
   int sqfree, tmonic;
 
@@ -3655,55 +4066,52 @@ polfnf(GEN a, GEN t)
   a = fix_relative_pol(t,a,0);
   alift = lift(a);
   p1 = content(alift); if (!gcmp1(p1)) { a = gdiv(a, p1); alift = lift(a); }
-  p1 = content(t); if (!gcmp1(t)) t = gdiv(t, p1);
+  t = primpart(t);
   tmonic = is_pm1(leading_term(t));
 
-  dent = ZX_disc(t); unt = gmodulsg(1,t);
-  u = nfgcd(alift,derivpol(alift), t, dent);
-  sqfree = gcmp1(u);
-  u = sqfree? alift: lift_intern(gdiv(a, gmul(unt,u)));
+  dent = indexpartial(t, NULL); unt = gmodulsg(1,t);
+  G = nfgcd(alift,derivpol(alift), t, dent);
+  sqfree = gcmp1(G);
+  u = sqfree? alift: lift_intern(gdiv(a, gmul(unt,G)));
   k = 0; n = ZY_ZXY_resultant(t, u, &k);
-  if (DEBUGLEVEL > 4) fprintferr("polfnf: choosing k = %ld\n",k);
+  if (DEBUGLEVEL>4) fprintferr("polfnf: choosing k = %ld\n",k);
+  if (!sqfree)
+  {
+    G = poleval(G, gadd(polx[varn(a)], gmulsg(k, polx[varn(t)])));
+    G = ZY_ZXY_resultant(t, G, NULL);
+  }
 
   /* n guaranteed to be squarefree */
-  fa = squff2(n,0); lx = lg(fa);
+  fa = DDF2(n,0); lx = lg(fa);
   y = cgetg(3,t_MAT);
   p1 = cgetg(lx,t_COL); y[1] = (long)p1;
   p2 = cgetg(lx,t_COL); y[2] = (long)p2;
+  if (lx == 2)
+  { /* P^k, k irreducible */
+    p1[1] = lmul(unt,u);
+    p2[1] = lstoi(degpol(a) / degpol(u));
+    return gerepilecopy(av, y);
+  }
   x0 = gadd(polx[varn(a)], gmulsg(-k, gmodulcp(polx[varn(t)], t)));
-  for (i=lx-1; i>1; i--)
+  for (i=lx-1; i>0; i--)
   {
-    GEN b, F = lift_intern(poleval((GEN)fa[i], x0));
-    if (!tmonic) F = gdiv(F, content(F));
+    GEN f = (GEN)fa[i], F = lift_intern(poleval(f, x0));
+    if (!tmonic) F = primpart(F);
     F = nfgcd(u, F, t, dent);
     if (typ(F) != t_POL || degpol(F) == 0)
       err(talker,"reducible modulus in factornf");
-    F = gmul(unt, F);
-    F = gdiv(F, leading_term(F));
-    u = lift_intern(gdeuc(u, F));
-    u = gdiv(u, content(u));
-    if (sqfree) e = 1;
-    else
+    e = 1;
+    if (!sqfree)
     {
-      e = 0; while (poldivis(a,F, &b)) { a = b; e++; }
-      if (degpol(a) == degpol(u)) sqfree = 1;
+      while (poldivis(G,f, &G)) e++;
+      if (degpol(G) == 0) sqfree = 1;
     }
-    p1[i] = (long)F;
+    p1[i] = ldiv(gmul(unt,F), leading_term(F));
     p2[i] = lstoi(e);
   }
-  u = gmul(unt, u);
-  u = gdiv(u, leading_term(u));
-  p1[1] = (long)u;
-  e = sqfree? 1: degpol(a)/degpol(u);
-  p2[1] = lstoi(e);
   return gerepilecopy(av, sort_factor(y, cmp_pol));
 }
 
-extern GEN FpXQX_safegcd(GEN P, GEN Q, GEN T, GEN p);
-extern GEN FpM(GEN z, GEN p);
-extern GEN polpol_to_mat(GEN v, long n);
-extern GEN mat_to_polpol(GEN x, long v, long w);
-
 static GEN
 to_frac(GEN a, GEN b)
 {
@@ -3736,7 +4144,7 @@ lift_to_frac(GEN t, GEN mod, GEN amax, GEN bmax, GEN d
 GEN
 matratlift(GEN M, GEN mod, GEN amax, GEN bmax, GEN denom)
 {
-  ulong ltop = avma;
+  gpmem_t ltop = avma;
   GEN N, a;
   long l2, l3;
   long i, j;
@@ -3759,7 +4167,7 @@ matratlift(GEN M, GEN mod, GEN amax, GEN bmax, GEN den
 GEN
 polratlift(GEN P, GEN mod, GEN amax, GEN bmax, GEN denom)
 {
-  ulong ltop = avma;
+  gpmem_t ltop = avma;
   GEN Q,a;
   long l2, j;
   if (typ(P)!=t_POL) err(typeer,"polratlift");
@@ -3793,13 +4201,11 @@ polratlift(GEN P, GEN mod, GEN amax, GEN bmax, GEN den
  * If not NULL, den must a a multiple of the denominator of the gcd,
  * for example the discriminant of nf.
  *
- * NOTE: if nf is not irreducible, nfgcd may loop forever, especially if the
- * gcd divides nf !
- */
+ * NOTE: if nf is not irreducible, nfgcd may loop forever, esp. if gcd | nf */
 GEN
 nfgcd(GEN P, GEN Q, GEN nf, GEN den)
 {
-  ulong ltop = avma;
+  gpmem_t ltop = avma;
   GEN sol, mod = NULL;
   long x = varn(P);
   long y = varn(nf);
@@ -3815,33 +4221,31 @@ nfgcd(GEN P, GEN Q, GEN nf, GEN den)
     den = mulii(den, mppgcd(ZX_resultant(lP, nf), ZX_resultant(lQ, nf)));
   /*Modular GCD*/
   {
-    ulong btop = avma;
-    ulong st_lim = stack_lim(btop, 1);
+    gpmem_t btop = avma, st_lim = stack_lim(btop, 1);
     long p;
     long dM=0, dR;
-    GEN M, dsol, dens;
+    GEN M, dsol;
     GEN R, ax, bo;
     byteptr primepointer;
-    for (p = 27449, primepointer = diffptr + 3000; ; p += *(primepointer++))
+    for (p = 27449, primepointer = diffptr + 3000; ; )
     {
       if (!*primepointer) err(primer1);
       if (!smodis(den, p))
-        continue;/*Discard primes dividing the disc(T) or leadingcoeff(PQ) */
-      if (DEBUGLEVEL>=5) fprintferr("nfgcd: p=%d\n",p);
+        goto repeat;/*Discard primes dividing disc(T) or leadingcoeff(PQ) */
+      if (DEBUGLEVEL>5) fprintferr("nfgcd: p=%d\n",p);
       if ((R = FpXQX_safegcd(P, Q, nf, stoi(p))) == NULL)
-        continue;/*Discard primes when modular gcd does not exist*/
+        goto repeat;/*Discard primes when modular gcd does not exist*/
       dR = degpol(R);
       if (dR == 0) return scalarpol(gun, x);
-      if (mod && dR > dM) continue; /* p divides Res(P/gcd, Q/gcd). Discard. */
+      if (mod && dR > dM) goto repeat; /* p divides Res(P/gcd, Q/gcd). Discard. */
 
       R = polpol_to_mat(R, d);
       /* previous primes divided Res(P/gcd, Q/gcd)? Discard them. */
-      if (!mod || dR < dM) { M = R; mod = stoi(p); dM = dR; continue; }
+      if (!mod || dR < dM) { M = R; mod = stoi(p); dM = dR; goto repeat; }
       if (low_stack(st_lim, stack_lim(btop, 1)))
       {
-	GEN *bptr[2]; bptr[0]=&M; bptr[1]=&mod;
 	if (DEBUGMEM>1) err(warnmem,"nfgcd");
-	gerepilemany(btop, bptr, 2);
+	gerepileall(btop, 2, &M, &mod);
       }
 
       ax = gmulgs(mpinvmod(stoi(p), mod), p);
@@ -3851,12 +4255,13 @@ nfgcd(GEN P, GEN Q, GEN nf, GEN den)
       /* I suspect it must be better to take amax > bmax*/
       bo = racine(shifti(mod, -1));
       if ((sol = matratlift(M, mod, bo, bo, den)) == NULL)
-        continue;
-      dens = denom(sol);
+        goto repeat;
       sol = mat_to_polpol(sol,x,y);
-      dsol = gmul(sol, gmodulcp(dens, nf));
-      if (gdivise(P, dsol) && gdivise(Q, dsol))
-        break;
+      dsol = primpart(sol);
+      if (gcmp0(pseudorem_i(P, dsol, nf))
+       && gcmp0(pseudorem_i(Q, dsol, nf))) break;
+    repeat:
+      NEXT_PRIME_VIADIFF_CHECK(p, primepointer);
     }
   }
   return gerepilecopy(ltop, sol);