===================================================================
RCS file: /home/cvs/OpenXM_contrib/pari-2.2/src/basemath/Attic/galconj.c,v
retrieving revision 1.2
retrieving revision 1.3
diff -u -p -r1.2 -r1.3
--- OpenXM_contrib/pari-2.2/src/basemath/Attic/galconj.c	2002/07/25 08:06:07	1.2
+++ OpenXM_contrib/pari-2.2/src/basemath/Attic/galconj.c	2002/09/11 07:26:50	1.3
@@ -1,4 +1,4 @@
-/* $Id: galconj.c,v 1.2 2002/07/25 08:06:07 noro Exp $
+/* $Id: galconj.c,v 1.3 2002/09/11 07:26:50 noro Exp $
 
 Copyright (C) 2000  The PARI group.
 
@@ -13,18 +13,26 @@ Check the License for details. You should have receive
 with the package; see the file 'COPYING'. If not, write to the Free Software
 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */
 
+#include "pari.h"
+extern GEN bezout_lift_fact(GEN T, GEN Tmod, GEN p, long e);
+extern GEN respm(GEN x,GEN y,GEN pm);
+extern GEN ZX_disc_all(GEN,long);
+extern GEN polratlift(GEN P, GEN mod, GEN amax, GEN bmax, GEN denom);
+extern GEN znstar_hnf_elts(GEN Z, GEN H);
+extern long ZX_get_prec(GEN x);
+
 /*************************************************************************/
 /**									**/
 /**                           GALOIS CONJUGATES        		        **/
 /**									**/
 /*************************************************************************/
-#include "pari.h"
-#define mylogint(x,y) logint((x),(y),NULL)
+
 GEN
 galoisconj(GEN nf)
 {
   GEN     x, y, z;
-  long    i, lz, v, av = avma;
+  long i, lz, v;
+  gpmem_t av = avma;
   nf = checknf(nf);
   x = (GEN) nf[1];
   v = varn(x);
@@ -51,7 +59,8 @@ galoisconj(GEN nf)
 GEN
 galoisconj2pol(GEN x, long nbmax, long prec)
 {
-  long    av = avma, i, n, v, nbauto;
+  long i, n, v, nbauto;
+  gpmem_t av = avma;
   GEN     y, w, polr, p1, p2;
   n = degpol(x);
   if (n <= 0)
@@ -92,7 +101,8 @@ galoisconj2pol(GEN x, long nbmax, long prec)
 GEN
 galoisconj2(GEN nf, long nbmax, long prec)
 {
-  long    av = avma, i, j, n, r1, ru, nbauto;
+  long i, j, n, r1, ru, nbauto;
+  gpmem_t av = avma;
   GEN     x, y, w, polr, p1, p2;
   if (typ(nf) == t_POL)
     return galoisconj2pol(nf, nbmax, prec);
@@ -156,37 +166,60 @@ galoisconj2(GEN nf, long nbmax, long prec)
   7: memory
   9: complete detail
 */
-/* retourne la permutation identite */
+
+/* DP = multiple of disc(P) or NULL
+ * Return a multiple of the denominator of an algebraic integer (in Q[X]/(P))
+ * when expressed in terms of the power basis */
 GEN
-permidentity(long l)
+indexpartial(GEN P, GEN DP)
 {
-  GEN     perm;
-  int     i;
-  perm = cgetg(l + 1, t_VECSMALL);
-  for (i = 1; i <= l; i++)
-    perm[i] = i;
-  return perm;
-}
+  gpmem_t av = avma;
+  long i, nb;
+  GEN fa, p1, res = gun, dP;
+  pari_timer T;
 
-GEN vandermondeinverseprepold(GEN T, GEN L)
-{
-  int     i, n = lg(L);
-  GEN     V, dT;
-  dT = derivpol(T);
-  V = cgetg(n, t_VEC);
-  for (i = 1; i < n; i++)
-    V[i] = (long) poleval(dT, (GEN) L[i]);
-  return V;
+  dP = derivpol(P);
+  if(DEBUGLEVEL>=5) (void)TIMER(&T);
+  if(!DP) DP = ZX_disc(P);
+  DP = mpabs(DP);
+  if(DEBUGLEVEL>=5) msgTIMER(&T,"IndexPartial: discriminant");
+  fa = auxdecomp(DP, 0);
+  if(DEBUGLEVEL>=5) msgTIMER(&T,"IndexPartial: factorization");
+  nb = lg(fa[1]);
+  for (i = 1; i < nb; i++)
+  {
+    GEN p=gmael(fa,1,i);
+    GEN e=gmael(fa,2,i);
+    p1 = powgi(p,shifti(e,-1));
+    if ( i==nb-1 )
+    {
+      if ( mod2(e) && !isprime(p) )
+	p1 = mulii(p1,p);
+    }
+    else if ( cmpis(e,4)>=0 )
+    {
+      if(DEBUGLEVEL>=5) fprintferr("IndexPartial: factor %Z ",p1);
+      p1 = mppgcd(p1, respm(P,dP,p1));
+      if(DEBUGLEVEL>=5) 
+      {
+        fprintferr("--> %Z : ",p1);
+        msgTIMER(&T,"");
+      }
+    }
+    res=mulii(res,p1);
+  }
+  return gerepileupto(av,res);
 }
 
-GEN vandermondeinverseprep(GEN T, GEN L)
+GEN
+vandermondeinverseprep(GEN L)
 {
   int     i, j, n = lg(L);
   GEN     V;
   V = cgetg(n, t_VEC);
   for (i = 1; i < n; i++)
   {
-    ulong ltop=avma;
+    gpmem_t ltop=avma;
     GEN W=cgetg(n,t_VEC);
     for (j = 1; j < n; j++)
       if (i==j)
@@ -202,23 +235,19 @@ GEN vandermondeinverseprep(GEN T, GEN L)
 GEN
 vandermondeinverse(GEN L, GEN T, GEN den, GEN prep)
 {
-  ulong   ltop = avma, lbot;
-  int     i, j, n = lg(L);
+  gpmem_t ltop = avma;
+  int     i, n = lg(L)-1;
   long    x = varn(T);
   GEN     M, P;
   if (!prep)
-    prep=vandermondeinverseprep(T,L);
-  M = cgetg(n, t_MAT);
-  for (i = 1; i < n; i++)
+    prep = vandermondeinverseprep(L);
+  M = cgetg(n+1, t_MAT);
+  for (i = 1; i <= n; i++)
   {
-    M[i] = lgetg(n, t_COL);
     P = gdiv(gdeuc(T, gsub(polx[x], (GEN) L[i])), (GEN) prep[i]);
-    for (j = 1; j < n; j++)
-      ((GEN *) M)[i][j] = P[1 + j];
+    M[i] = (long)pol_to_vec(P,n);
   }
-  lbot = avma;
-  M = gmul(den, M);
-  return gerepile(ltop, lbot, M);
+  return gerepileupto(ltop, gmul(den, M));
 }
 
 /* Calcule les bornes sur les coefficients a chercher */
@@ -234,88 +263,111 @@ struct galois_borne
 };
 
 
-GEN indexpartial(GEN,GEN);
-GEN ZX_disc_all(GEN,long);
-
 GEN
-galoisborne(GEN T, GEN dn, struct galois_borne *gb, long ppp)
+initgaloisborne(GEN T, GEN dn, long prec, GEN *ptL, GEN *ptprep, GEN *ptdis)
 {
-  ulong   ltop = avma, av2;
-  GEN     borne, borneroots, borneabs;
-  int     i, j;
-  int     n;
-  GEN     L, M, z, prep, den;
-  long    prec;
-  n = degpol(T);
-  prec = 1;
-  for (i = 2; i < lgef(T); i++)
-    if (lg(T[i]) > prec)
-      prec = lg(T[i]);
-  /*prec++;*/
-  if (DEBUGLEVEL>=4) gentimer(3);
+  int i, n = degpol(T);
+  GEN L, z, prep, den;
+  pari_timer ti;
+
+  if (DEBUGLEVEL>=4) (void)TIMER(&ti);
   L = roots(T, prec);
-  if (DEBUGLEVEL>=4) genmsgtimer(3,"roots");
+  if (DEBUGLEVEL>=4) msgTIMER(&ti,"roots");
   for (i = 1; i <= n; i++)
   {
     z = (GEN) L[i];
-    if (signe(z[2]))
-      break;
+    if (signe(z[2])) break;
     L[i] = z[1];
   }
-  if (DEBUGLEVEL>=4) gentimer(3);
-  prep=vandermondeinverseprep(T, L);
+  if (DEBUGLEVEL>=4) (void)TIMER(&ti);
+  prep = vandermondeinverseprep(L);
   if (!dn)
   {
-    GEN res = divide_conquer_prod(gabs(prep,prec), mpmul);
-    GEN dis;
+    GEN dis, res = divide_conquer_prod(gabs(prep,prec), mpmul);
     disable_dbg(0);
-    dis = ZX_disc_all(T, 1+mylogint(res,gdeux));
+    dis = ZX_disc_all(T, 1+logint(res,gdeux,NULL));
     disable_dbg(-1);
-    den = gclone(indexpartial(T,dis));
+    den = indexpartial(T,dis);
+    if (ptdis) *ptdis = dis;
   }
   else
-    den=dn;
-  M = vandermondeinverse(L, gmul(T, realun(prec)), den, prep);
-  if (DEBUGLEVEL>=4) genmsgtimer(3,"vandermondeinverse");
-  borne = realzero(prec);
-  for (i = 1; i <= n; i++)
+    den = dn;
+  if (ptprep) *ptprep = prep;
+  *ptL = L; return den;
+}
+
+/* ||| M ||| with respect to || x ||_oo. Assume M square t_MAT */
+GEN
+matrixnorm(GEN M, long prec)
+{
+  long i,j, n = lg(M);
+  GEN B = realzero(prec);
+
+  for (i = 1; i < n; i++)
   {
-    z = gzero;
-    for (j = 1; j <= n; j++)
+    GEN z = gabs(gcoeff(M,i,1), prec);
+    for (j = 2; j < n; j++)
       z = gadd(z, gabs(gcoeff(M,i,j), prec));
-    if (gcmp(z, borne) > 0)
-      borne = z;
+    if (gcmp(z, B) > 0) B = z;
   }
-  borneroots = realzero(prec);
-  for (i = 1; i <= n; i++)
+  return B;
+}
+
+/* L a t_VEC/t_COL, return ||L||_oo */
+GEN
+supnorm(GEN L, long prec)
+{
+  long i, n = lg(L);
+  GEN z, B;
+
+  if (n == 1) return realzero(prec);
+  B = gabs((GEN)L[1], prec);
+  for (i = 2; i < n; i++)
   {
-    z = gabs((GEN) L[i], prec);
-    if (gcmp(z, borneroots) > 0)
-      borneroots = z;
+    z = gabs((GEN)L[i], prec);
+    if (gcmp(z, B) > 0) B = z;
   }
+  return B;
+}
+
+GEN
+galoisborne(GEN T, GEN dn, struct galois_borne *gb, long ppp)
+{
+  gpmem_t ltop = avma, av2;
+  GEN     borne, borneroots, borneabs;
+  int     n;
+  GEN     L, M, prep, den;
+  long    prec;
+  pari_timer ti;
+
+  prec = ZX_get_prec(T);
+  den = initgaloisborne(T,dn,prec, &L,&prep,NULL);
+  if (!dn) den = gclone(den);
+  if (DEBUGLEVEL>=4) TIMERstart(&ti);
+  M = vandermondeinverse(L, gmul(T, realun(prec)), den, prep);
+  if (DEBUGLEVEL>=4) msgTIMER(&ti,"vandermondeinverse");
+  borne = matrixnorm(M, prec);
+  borneroots = supnorm(L, prec);
+  n = degpol(T);
   borneabs = addsr(1, gmulsg(n, gpowgs(borneroots, n/ppp)));
   /*if (ppp == 1)
     borneabs = addsr(1, gmulsg(n, gpowgs(borneabs, 2)));*/
   borneroots = addsr(1, gmul(borne, borneroots));
   av2 = avma;
   /*We use d-1 test, so we must overlift to 2^BITS_IN_LONG*/
-  gb->valsol = mylogint(gmul2n(borneroots,2+BITS_IN_LONG), gb->l);
-  gb->valabs = mylogint(gmul2n(borneabs,2), gb->l);
-  gb->valabs = max(gb->valsol,gb->valabs);
+  gb->valsol = logint(gmul2n(borneroots,2+BITS_IN_LONG), gb->l,NULL);
+  gb->valabs = logint(gmul2n(borneabs,2), gb->l,NULL);
+  gb->valabs = max(gb->valsol, gb->valabs);
   if (DEBUGLEVEL >= 4)
     fprintferr("GaloisConj:val1=%ld val2=%ld\n", gb->valsol, gb->valabs);
   avma = av2;
-  gb->bornesol = gerepileupto(ltop, ceil_safe(borneroots));
+  gb->bornesol = gerepileupto(ltop, ceil_safe(mulrs(borneroots,2)));
   if (DEBUGLEVEL >= 9)
     fprintferr("GaloisConj: Bound %Z\n",borneroots);
   gb->ladicsol = gpowgs(gb->l, gb->valsol);
   gb->ladicabs = gpowgs(gb->l, gb->valabs);
   gb->lbornesol = subii(gb->ladicsol,gb->bornesol);
-  if (!dn)
-  {
-    dn=forcecopy(den);
-    gunclone(den);
-  }
+  if (!dn) { dn = icopy(den); gunclone(den); }
   return dn;
 }
 
@@ -334,7 +386,7 @@ struct galois_lift
 /* Initialise la structure galois_lift */
 GEN makeLden(GEN L,GEN den, struct galois_borne *gb)
 {
-  ulong ltop=avma;
+  gpmem_t ltop=avma;
   long i,l=lg(L);
   GEN Lden;
   Lden=cgetg(l,t_VEC);
@@ -347,7 +399,7 @@ GEN makeLden(GEN L,GEN den, struct galois_borne *gb)
 void
 initlift(GEN T, GEN den, GEN p, GEN L, GEN Lden, struct galois_borne *gb, struct galois_lift *gl)
 {
-  ulong   ltop, lbot;
+  gpmem_t ltop, lbot;
   gl->gb=gb;
   gl->T = T;
   gl->den = den;
@@ -355,109 +407,21 @@ initlift(GEN T, GEN den, GEN p, GEN L, GEN Lden, struc
   gl->L = L;
   gl->Lden = Lden;
   ltop = avma;
-  gl->e = mylogint(gmul2n(gb->bornesol, 2+BITS_IN_LONG),p);
+  gl->e = logint(gmul2n(gb->bornesol, 2+BITS_IN_LONG),p,NULL);
   gl->e = max(2,gl->e);
   lbot = avma;
   gl->Q = gpowgs(p, gl->e);
   gl->Q = gerepile(ltop, lbot, gl->Q);
   gl->TQ = FpX_red(T,gl->Q);
 }
-#if 0
-/*
- * T est le polynome \sum t_i*X^i S est un Polmod T 
- * Retourne un vecteur Spow
- * verifiant la condition: i>=1 et t_i!=0 ==> Spow[i]=S^(i-1)*t_i 
- * Essaie d'etre efficace sur les polynomes lacunaires
- */
-GEN
-compoTS(GEN T, GEN S, GEN Q, GEN mod)
-{
-  GEN     Spow;
-  int     i;
-  Spow = cgetg(lgef(T) - 2, t_VEC);
-  for (i = 3; i < lg(Spow); i++)
-    Spow[i] = 0;
-  Spow[1] = (long) polun[varn(S)];
-  Spow[2] = (long) S;
-  for (i = 2; i < lg(Spow) - 1; i++)
-  {
-    if (signe((GEN) T[i + 3]))
-      for (;;)
-      {
-	int     k, l;
-	for (k = 1; k <= (i >> 1); k++)
-	  if (Spow[k + 1] && Spow[i - k + 1])
-	    break;
-	if ((k << 1) < i)
-	{
-	  Spow[i + 1] = (long) FpXQ_mul((GEN) Spow[k + 1], (GEN) Spow[i - k + 1],Q,mod);
-	  break;
-	}
-	else if ((k << 1) == i)	
-	{
-	  Spow[i + 1] = (long) FpXQ_sqr((GEN) Spow[k + 1],Q,mod);
-	  break;
-	}
-	for (k = i - 1; k > 0; k--)
-	  if (Spow[k + 1])
-	    break;
-	if ((k << 1) < i)
-	{
-	  Spow[(k << 1) + 1] = (long) FpXQ_sqr((GEN) Spow[k + 1],Q,mod);
-	  continue;
-	}
-	for (l = i - k; l > 0; l--)
-	  if (Spow[l + 1])
-	    break;
-	if (Spow[i - l - k + 1])
-	  Spow[i - k + 1] = (long) FpXQ_mul((GEN) Spow[i - l - k + 1], (GEN) Spow[l + 1],Q,mod);
-	else
-	  Spow[l + k + 1] = (long) FpXQ_mul((GEN) Spow[k + 1], (GEN) Spow[l + 1],Q,mod);
-      }
-  }
-  for (i = 1; i < lg(Spow); i++)
-    if (signe((GEN) T[i + 2]))
-      Spow[i] = (long) FpX_Fp_mul((GEN) Spow[i], (GEN) T[i + 2],mod);
-  return Spow;
-}
 
 /*
- * Calcule T(S) a l'aide du vecteur Spow
- */
-static GEN
-calcTS(GEN Spow, GEN P, GEN S, GEN Q, GEN mod)
-{
-  GEN     res = gzero;
-  int     i;
-  for (i = 1; i < lg(Spow); i++)
-    if (signe((GEN) P[i + 2]))
-      res = FpX_add(res, (GEN) Spow[i],NULL);
-  res = FpXQ_mul(res,S,Q,mod);
-  res=FpX_Fp_add(res,(GEN)P[2],mod);
-  return res;
-}
-
-/*
- * Calcule T'(S) a l'aide du vecteur Spow
- */
-static GEN
-calcderivTS(GEN Spow, GEN P, GEN mod)
-{
-  GEN     res = gzero;
-  int     i;
-  for (i = 1; i < lg(Spow); i++)
-    if (signe((GEN) P[i + 2]))
-      res = FpX_add(res, FpX_Fp_mul((GEN) Spow[i], stoi(i),mod),NULL);
-  return FpX_red(res,mod);
-}
-#endif
-/*
  * Verifie que f est une solution presque surement et calcule sa permutation
  */
 static int
 poltopermtest(GEN f, struct galois_lift *gl, GEN pf)
 {
-  ulong   ltop;
+  gpmem_t ltop;
   GEN     fx, fp;
   long     i, j,ll;
   for (i = 2; i< lgef(f); i++)
@@ -466,6 +430,8 @@ poltopermtest(GEN f, struct galois_lift *gl, GEN pf)
     {
       if (DEBUGLEVEL>=4)
 	fprintferr("GaloisConj: Solution too large, discard it.\n");
+      if (DEBUGLEVEL>=8)
+	fprintferr("f=%Z\n borne=%Z\n l-borne=%Z\n",f,gl->gb->bornesol,gl->gb->lbornesol);
       return 0;
     }
   ll=lg(gl->L);
@@ -492,8 +458,6 @@ poltopermtest(GEN f, struct galois_lift *gl, GEN pf)
   return 1;
 }
 
-GEN polratlift(GEN P, GEN mod, GEN amax, GEN bmax, GEN denom);
-
 /*
  * Soit P un polynome de \ZZ[X] , p un nombre premier , S\in\FF_p[X]/(Q) tel
  * que P(S)=0 [p,Q] Relever S en S_0 tel que P(S_0)=0 [p^e,Q]
@@ -501,7 +465,7 @@ GEN polratlift(GEN P, GEN mod, GEN amax, GEN bmax, GEN
  */
 static long monoratlift(GEN S, GEN q, GEN qm1old,struct galois_lift *gl, GEN frob)
 {
-  ulong ltop=avma;
+  gpmem_t ltop=avma;
   GEN tlift = polratlift(S,q,qm1old,qm1old,gl->den);
   if (tlift)
   {
@@ -525,7 +489,7 @@ static long monoratlift(GEN S, GEN q, GEN qm1old,struc
   GEN
 monomorphismratlift(GEN P, GEN S, struct galois_lift *gl, GEN frob)
 {
-  ulong   ltop, lbot;
+  gpmem_t ltop, lbot;
   long rt;
   GEN     Q=gl->T, p=gl->p;
   long    e=gl->e, level=1;
@@ -534,8 +498,7 @@ monomorphismratlift(GEN P, GEN S, struct galois_lift *
   GEN     W, Pr, Qr, Sr, Wr = gzero, Prold, Qrold, Spow;
   long    i,nb,mask;
   GEN    *gptr[2];
-  if (DEBUGLEVEL == 1)
-    timer2();
+  if (DEBUGLEVEL == 1) (void)timer2();
   x = varn(P);
   rt = brent_kung_optpow(degpol(Q),1);
   q = p; qm1 = gun; /*during the run, we have p*qm1=q*/
@@ -554,7 +517,7 @@ monomorphismratlift(GEN P, GEN S, struct galois_lift *
     if (DEBUGLEVEL>=2)
     {
       level=(level<<1)-((mask>>i)&1);
-      timer2();
+      (void)timer2();
     }
     qm1 = (mask&(1<<i))?sqri(qm1):mulii(qm1, q);
     q   =  mulii(qm1, p);
@@ -562,19 +525,16 @@ monomorphismratlift(GEN P, GEN S, struct galois_lift *
     Qr = (P==Q)?Pr:FpX_red(Q, q);/*A little speed up for automorphismlift*/
     ltop = avma;
     Sr = S;
-    /*Spow = compoTS(Pr, Sr, Qr, q);*/
     Spow = FpXQ_powers(Sr, rt, Qr, q);
 
     if (i)
     {
-      /*W = FpXQ_mul(Wr, calcderivTS(Spow, Prold,qold), Qrold, qold);*/
       W = FpXQ_mul(Wr, FpX_FpXQV_compo(deriv(Pr,-1),FpXV_red(Spow,qold),Qrold,qold), Qrold, qold);
       W = FpX_neg(W, qold);
       W = FpX_Fp_add(W, gdeux, qold);
       W = FpXQ_mul(Wr, W, Qrold, qold);
     }
     Wr = W;
-    /*S = FpXQ_mul(Wr, calcTS(Spow, Pr, Sr, Qr, q),Qr,q);*/
     S = FpXQ_mul(Wr, FpX_FpXQV_compo(Pr, Spow, Qr, q),Qr,q);
     S = FpX_sub(Sr, S, NULL);
     lbot = avma;
@@ -621,11 +581,10 @@ struct galois_testlift
   GEN     C;
   GEN     Cd;
 };
-GEN bezout_lift_fact(GEN T, GEN Tmod, GEN p, long e);
 static GEN
 galoisdolift(struct galois_lift *gl, GEN frob)
 {
-  ulong ltop=avma;
+  gpmem_t ltop=avma;
   long v = varn(gl->T);
   GEN Tp = FpX_red(gl->T, gl->p);
   GEN S = FpXQ_pow(polx[v],gl->p, Tp,gl->p);
@@ -647,11 +606,11 @@ inittestlift( GEN plift, GEN Tmod, struct galois_lift 
   gt->pauto[2] = (long) gcopy(plift);
   if (gt->f > 2)
   {
-    ulong ltop=avma;
+    gpmem_t ltop=avma;
     int i;
     long nautpow=brent_kung_optpow(gt->n-1,gt->f-2);
     GEN autpow;
-    if (DEBUGLEVEL >= 1) timer2();
+    if (DEBUGLEVEL >= 1) (void)timer2();
     autpow = FpXQ_powers(plift,nautpow,gl->TQ,gl->Q);
     for (i = 3; i <= gt->f; i++)
       gt->pauto[i] = (long) FpX_FpXQV_compo((GEN)gt->pauto[i-1],autpow,gl->TQ,gl->Q);
@@ -663,7 +622,7 @@ inittestlift( GEN plift, GEN Tmod, struct galois_lift 
 
 long intheadlong(GEN x, GEN mod)
 {
-  ulong ltop=avma;
+  gpmem_t ltop=avma;
   GEN r;
   int s;
   long res;
@@ -698,15 +657,15 @@ long
 frobeniusliftall(GEN sg, long el, GEN *psi, struct galois_lift *gl,
 		 struct galois_testlift *gt, GEN frob)
 {
-  ulong   ltop = avma, av, ltop2;
+  gpmem_t av, ltop2, ltop = avma;
   long    d, z, m, c, n, ord;
   int     i, j, k;
   GEN     pf, u, v;
-  GEN     C, Cd;
-  long 	  Z, Zv;
-  long    stop=0;
-  GEN NN,NQ,NR;
-  long N1,N2,R1,Ni;
+  GEN     C, Cd, sgi, cache;
+  long    Z, c_idx=gt->g-1;
+  long    stop=0,hop=0;
+  GEN     NN,NQ,NR;
+  long    N1,N2,R1,Ni;
   m = gt->g;
   ord = gt->f;
   n = lg(gl->L) - 1;
@@ -722,8 +681,8 @@ frobeniusliftall(GEN sg, long el, GEN *psi, struct gal
   NQ=dvmdis(NN,N1,&NR);
   if (cmpis(NQ,1000000000)>0)
   {
-    err(warner,"Combinatorics too hard : would need %Z tests!\n \
-	I will skip it,but it may induce galoisinit to loop",NN);
+    err(warner,"Combinatorics too hard : would need %Z tests!\n"
+	"I will skip it,but it may induce galoisinit to loop",NN);
     avma = ltop;
     *psi = NULL;
     return 0;
@@ -732,55 +691,65 @@ frobeniusliftall(GEN sg, long el, GEN *psi, struct gal
   if (DEBUGLEVEL>=4)
   {
     stop=N1/20;
-    timer2();
+    (void)timer2();
   }
   avma = ltop2;
   C=gt->C;
   Cd=gt->Cd;
   v = FpX_Fp_mul(FpXQ_mul((GEN)gt->pauto[1+el%ord], (GEN)
 	gt->bezoutcoeff[m],gl->TQ,gl->Q),gl->den,gl->Q);
-  Zv=polheadlong(v,0,gl->Q);
+  sgi=cgetg(lg(sg),t_VECSMALL);
+  for(i=1;i<lg(sgi);i++)
+    sgi[i]=(el*sg[i])%ord + 1;
+  cache=cgetg(m+1,t_VECSMALL);
+  cache[m]=polheadlong(v,1,gl->Q);
+  Z=polheadlong(v,2,gl->Q);
   for (i = 1; i < m; i++)
     pf[i] = 1 + i / d;
   av = avma;
   for (Ni = 0, i = 0; ;i++)
   {
-    Z=Zv;
-    for (j = 1; j < m; j++)
+    for (j = c_idx ; j > 0; j--)
     {
       long h;
-      h=(el*sg[pf[j]])%ord + 1;
+      h=sgi[pf[j]];
       if (!mael(C,h,j))
       {
-	ulong av3=avma;
+	gpmem_t av3=avma;
 	GEN r;
 	r=FpX_Fp_mul(FpXQ_mul((GEN) gt->pauto[h], 
 	      (GEN) gt->bezoutcoeff[j],gl->TQ,gl->Q),gl->den,gl->Q);
 	mael(C,h,j)  = lclone(r);
-	mael(Cd,h,j) = polheadlong(r,0,gl->Q);
+	mael(Cd,h,j) = polheadlong(r,1,gl->Q);
 	avma=av3;
       }
-      Z += mael(Cd,h,j);
+      cache[j]=cache[j+1]+mael(Cd,h,j);
     }
-    if (labs(Z)<=n )
+    if (labs(cache[1])<=n )
     {
-      Z=polheadlong(v,1,gl->Q);
+      long ZZ=Z;
       for (j = 1; j < m; j++)
-	Z += polheadlong(gmael(C,(el*sg[pf[j]])%ord + 1,j),1,gl->Q);
-      if (labs(Z)<=n )
+	ZZ += polheadlong(gmael(C,sgi[pf[j]],j),2,gl->Q);
+      if (labs(ZZ)<=n )
       {
 	u = v;
 	for (j = 1; j < m; j++)
-	  u = FpX_add(u, gmael(C,1+(el*sg[pf[j]])%ord,j),NULL);
+	  u = FpX_add(u, gmael(C,sgi[pf[j]],j),NULL);
 	u = FpX_center(FpX_red(u, gl->Q), gl->Q);
 	if (poltopermtest(u, gl, frob))
 	{
 	  if (DEBUGLEVEL >= 4 )
+	  {
 	    msgtimer("");
+	    fprintferr("GaloisConj: %d hops on %Z tests\n",hop,addis(mulss(Ni,N1),i));
+	  }
 	  avma = ltop2;
 	  return 1;
 	}
+	else if (DEBUGLEVEL >= 4 )
+	  fprintferr("M");
       }
+      else hop++;
     }
     if (DEBUGLEVEL >= 4 && i==stop)
     {
@@ -799,7 +768,7 @@ frobeniusliftall(GEN sg, long el, GEN *psi, struct gal
       if (DEBUGLEVEL>=4)
       {
 	stop=N1/20;
-	timer2();
+	(void)timer2();
       }
     }
     for (j = 2; j < m && pf[j - 1] >= pf[j]; j++);
@@ -813,27 +782,14 @@ frobeniusliftall(GEN sg, long el, GEN *psi, struct gal
     z = pf[j];
     pf[j] = pf[k];
     pf[k] = z;
+    c_idx=j;
   }
+  if (DEBUGLEVEL>=4)
+    fprintferr("GaloisConj: not found, %d hops \n",hop);
   avma = ltop;
   *psi = NULL;
   return 0;
 }
-/*
- * applique une permutation t a un vecteur s sans duplication
- * Propre si s est un VECSMALL
- */
-static GEN
-permapply(GEN s, GEN t)
-{
-  GEN     u;
-  int     i;
-  if (lg(s) < lg(t))
-    err(talker, "First permutation shorter than second in permapply");
-  u = cgetg(lg(s), typ(s));
-  for (i = 1; i < lg(s); i++)
-    u[i] = s[t[i]];
-  return u;
-}
 
 /*
  * alloue une ardoise pour n entiers de longueurs pour les test de
@@ -870,7 +826,7 @@ struct galois_test
 static GEN
 Vmatrix(long n, struct galois_test *td)
 {
-  ulong   ltop = avma, lbot;
+  gpmem_t lbot, ltop = avma;
   GEN     V;
   long    i;
   V = cgetg(lg(td->L), t_VEC);
@@ -888,7 +844,7 @@ Vmatrix(long n, struct galois_test *td)
 static void
 inittest(GEN L, GEN M, GEN borne, GEN ladic, struct galois_test *td)
 {
-  ulong   ltop;
+  gpmem_t ltop;
   long    i;
   int     n = lg(L) - 1;
   if (DEBUGLEVEL >= 8)
@@ -909,7 +865,7 @@ inittest(GEN L, GEN M, GEN borne, GEN ladic, struct ga
   ltop = avma;
   td->PV[td->ordre[n]] = (long) gclone(Vmatrix(td->ordre[n], td));
   avma = ltop;
-  td->TM = gtrans(M);
+  td->TM = gtrans_i(M);
   settyp(td->TM, t_VEC);
   for (i = 1; i < lg(td->TM); i++)
     settyp(td->TM[i], t_VEC);
@@ -942,7 +898,7 @@ static long
 padicisint(GEN P, struct galois_test *td)
 {
   long    r;
-  ulong   ltop = avma;
+  gpmem_t ltop = avma;
   GEN     U;
   /*if (typ(P) != t_INT)    err(typeer, "padicisint");*/
   U = modii(P, td->ladic);
@@ -957,13 +913,13 @@ padicisint(GEN P, struct galois_test *td)
 static long
 verifietest(GEN pf, struct galois_test *td)
 {
-  ulong   av = avma;
+  gpmem_t av = avma;
   GEN     P, V;
   int     i, j;
   int     n = lg(td->L) - 1;
   if (DEBUGLEVEL >= 8)
     fprintferr("GaloisConj:Entree Verifie Test\n");
-  P = permapply(td->L, pf);
+  P = perm_mul(td->L, pf);
   for (i = 1; i < n; i++)
   {
     long    ord;
@@ -1028,7 +984,7 @@ GEN
 testpermutation(GEN F, GEN B, long s, long t, long cut,
 		struct galois_test *td)
 {
-  ulong   av, avm = avma;
+  gpmem_t av, avm = avma;
   int     a, b, c, d, n;
   GEN     pf, x, ar, y, *G;
   int     p1, p2, p3, p4, p5, p6;
@@ -1036,8 +992,7 @@ testpermutation(GEN F, GEN B, long s, long t, long cut
   long    i, j, cx, hop = 0, start = 0;
   GEN     W, NN, NQ, NR;
   long    V;
-  if (DEBUGLEVEL >= 1)
-    timer2();
+  if (DEBUGLEVEL >= 1) (void)timer2();
   a = lg(F) - 1;
   b = lg(F[1]) - 1;
   c = lg(B) - 1;
@@ -1071,7 +1026,7 @@ testpermutation(GEN F, GEN B, long s, long t, long cut
   {
     avma=avm;
     err(warner,"Combinatorics too hard : would need %Z tests!\n I'll skip it but you will get a partial result...",NN);
-    return permidentity(n);
+    return perm_identity(n);
   }
   N2=itos(NQ); R1=itos(NR);
   for (l2 = 0; l2 <= N2; l2++)
@@ -1195,106 +1150,20 @@ testpermutation(GEN F, GEN B, long s, long t, long cut
   avma = avm;
   return gzero;
 }
-/* Compute generators for the subgroup of (Z/nZ)* given in HNF. 
- * I apologize for the following spec:
- * If zns=znstar(2) then
- * zn2=gtovecsmall((GEN)zns[2]);
- * zn3=lift((GEN)zns[3]);
- * gen and ord : VECSMALL of length lg(zn3).
- * the result is in gen. 
- * ord contains the relatives orders of the generators.
- */
-GEN hnftogeneratorslist(long n, GEN zn2, GEN zn3, GEN lss, GEN gen, GEN ord)
-{
-  ulong ltop=avma;
-  int j,h;
-  GEN m=stoi(n);
-  for (j = 1; j < lg(gen); j++)
-  {
-    gen[j] = 1;
-    for (h = 1; h < lg(lss); h++)
-      gen[j] = mulssmod(gen[j], itos(powmodulo((GEN)zn3[h],gmael(lss,j,h),m)),n);
-    ord[j] = zn2[j] / itos(gmael(lss,j,j));
-  }
-  avma=ltop;
-  return gen;
-}
-/*in place sort. Return V for convenience only*/
-GEN sortvecsmall(GEN V)
-{
-  long i,j,k,l,m;
-  for(l=1;l<lg(V);l<<=1)
-    for(j=1;(j>>1)<lg(V);j+=l<<1)
-      for(k=j+l,i=j; i<k && k<j+(l<<1) && k<lg(V);)
-	if (V[i]>V[k])
-	{
-	  long z=V[k];
-	  for (m=k;m>i;m--)
-	    V[m]=V[m-1];
-	  V[i]=z;
-	  k++;
-	}
-	else
-	  i++;
-  return V ;
-}
-GEN uniqvecsmall(GEN V)
-{
-  ulong ltop=avma;
-  GEN W;
-  long i,j;
-  if ( lg(V) == 1 ) return gcopy(V);
-  W=cgetg(lg(V),t_VECSMALL);
-  W[1]=V[1];
-  for(i=2,j=1;i<lg(V);i++)
-    if (V[i]!=W[j])
-      W[++j]=V[i];
-  setlg(W,j+1);
-  return gerepileupto(ltop,W);
-}
 int pari_compare_lg(GEN x, GEN y)
 {
   return lg(x)-lg(y);
 }
-/* Compute the elements of the subgroup of (Z/nZ)* given in HNF. 
- * I apologize for the following spec:
- * If zns=znstar(2) then
- * zn2=gtovecsmall((GEN)zns[2]);
- * zn3=lift((GEN)zns[3]);
- * card=cardinal of the subgroup(i.e phi(n)/det(lss))
- */
-GEN
-hnftoelementslist(long n, GEN zn2, GEN zn3, GEN lss, long card)
-{
-  ulong   ltop;
-  GEN     sg, gen, ord;
-  int     k, j, l;
-  sg = cgetg(1 + card, t_VECSMALL);
-  ltop=avma;
-  gen = cgetg(lg(zn3), t_VECSMALL);
-  ord = cgetg(lg(zn3), t_VECSMALL);
-  sg[1] = 1;
-  hnftogeneratorslist(n,zn2,zn3,lss,gen,ord);
-  for (j = 1, l = 1; j < lg(gen); j++)
-  {
-    int     c = l * (ord[j] - 1);
-    for (k = 1; k <= c; k++)	/* I like it */
-      sg[++l] = mulssmod(sg[k], gen[j], n);
-  }
-  sortvecsmall(sg);
-  avma=ltop;
-  return sg;
-}
 
 /*
- * Calcule la liste des sous groupes de \ZZ/m\ZZ d'ordre divisant p et
+ * Calcule la liste des sous groupes de (\ZZ/m\ZZ)^* d'ordre divisant p et
  * retourne la liste de leurs elements
  */
 GEN
 listsousgroupes(long m, long p)
 {
-  ulong   ltop = avma;
-  GEN     zns, lss, sg, res, zn2, zn3;
+  gpmem_t ltop = avma;
+  GEN     zn, zns, lss, sg, res;
   int     k, card, i, phi;
   if (m == 2)
   {
@@ -1304,89 +1173,29 @@ listsousgroupes(long m, long p)
     sg[1] = 1;
     return res;
   }
-  zns = znstar(stoi(m));
-  phi = itos((GEN) zns[1]);
-  zn2 = gtovecsmall((GEN)zns[2]);
-  zn3 = lift((GEN)zns[3]);
-  lss = subgrouplist((GEN) zns[2], 0);
+  zn = znstar(stoi(m));
+  phi = itos((GEN) zn[1]);
+  zns = znstar_small(zn);
+  lss = subgrouplist((GEN) zn[2], NULL);
   res = cgetg(lg(lss), t_VEC);
   for (k = 1, i = lg(lss) - 1; i >= 1; i--)
   {
-    long    av;
+    gpmem_t av;
     av = avma;
     card = phi / itos(det((GEN) lss[i]));
     avma = av;
     if (p % card == 0)
-      res[k++] = (long) hnftoelementslist(m,zn2,zn3,(GEN)lss[i],card);
+      res[k++] = (long) znstar_hnf_elts(zns,(GEN)lss[i]);
   }
   setlg(res,k);
   res=gen_sort(res,0,&pari_compare_lg);
   return gerepileupto(ltop,res);
 }
 
-/* Compute the orbits decomposition of a permutation
- * or of a set of permutations given by a vector.
- * v must not be an empty vector, becaude there is no way then
- * to tell the length of the permutation.
- */
-
 GEN
-permorbite(GEN v)
-{
-  ulong   ltop = avma, lbot;
-  int     i, j, k, l, m, n, o, p, flag;
-  GEN     bit, cycle, cy, u;
-  if (typ(v) == t_VECSMALL)
-  {
-    u = cgetg(2, t_VEC);
-    u[1] = (long) v;
-    v = u;
-  }
-  n = lg(v[1]);
-  cycle = cgetg(n, t_VEC);
-  bit = cgetg(n, t_VECSMALL);
-  for (i = 1; i < n; i++)
-    bit[i] = 0;
-  for (k = 1, l = 1; k < n;)
-  {
-    for (j = 1; bit[j]; j++);
-    cy = cgetg(n, t_VECSMALL);
-    m = 1;
-    cy[m++] = j;
-    bit[j] = 1;
-    k++;
-    do
-    {
-      flag = 0;
-      for (o = 1; o < lg(v); o++)
-      {
-	for (p = 1; p < m; p++)	/* m varie! */
-	{
-	  j = ((long **) v)[o][cy[p]];
-	  if (!bit[j])
-	  {
-	    flag = 1;
-	    bit[j] = 1;
-	    k++;
-	    cy[m++] = j;
-	  }
-	}
-      }
-    }
-    while (flag);
-    setlg(cy, m);
-    cycle[l++] = (long) cy;
-  }
-  setlg(cycle, l);
-  lbot = avma;
-  cycle = gcopy(cycle);
-  return gerepile(ltop, lbot, cycle);
-}
-
-GEN
 fixedfieldpolsigma(GEN sigma, GEN p, GEN Tp, GEN sym, GEN deg, long g)
 {
-  ulong ltop=avma;
+  gpmem_t ltop=avma;
   long i, j, npows;
   GEN  s, f, pows;
   sigma=lift(gmul(sigma,gmodulsg(1,p)));
@@ -1418,14 +1227,14 @@ fixedfieldfactmod(GEN Sp, GEN p, GEN Tmod)
   long l=lg(Tmod);
   GEN F=cgetg(l,t_VEC);
   for(i=1;i<l;i++)
-    F[i]=(long)FpXQ_minpoly(Sp, (GEN) Tmod[i],p);
+    F[i]=(long)FpXQ_minpoly(FpX_res(Sp,(GEN) Tmod[i],p), (GEN) Tmod[i],p);
   return F;
 }
 
 GEN
 fixedfieldnewtonsumaut(GEN sigma, GEN p, GEN Tp, GEN e, long g)
 {
-  ulong ltop=avma;
+  gpmem_t ltop=avma;
   long i;
   GEN s,f,V;
   long rt;
@@ -1453,7 +1262,7 @@ fixedfieldnewtonsum(GEN O, GEN L, GEN mod, GEN e)
   PL=cgetg(lg(O), t_COL);
   for(i=1; i<=f; i++)
   {
-    ulong ltop=avma; 
+    gpmem_t ltop=avma;
     GEN s=gzero;
     for(j=1; j<=g; j++)
       s=addii(s,powmodulo((GEN)L[mael(O,i,j)],e,mod));
@@ -1465,7 +1274,7 @@ fixedfieldnewtonsum(GEN O, GEN L, GEN mod, GEN e)
 GEN
 fixedfieldpol(GEN O, GEN L, GEN mod, GEN sym, GEN deg)
 {
-  ulong ltop=avma;
+  gpmem_t ltop=avma;
   long i;
   GEN S=gzero;
   for(i=1;i<lg(sym);i++)
@@ -1531,7 +1340,7 @@ fixedfieldsurmer(GEN O, GEN L, GEN mod, GEN l, GEN p, 
     long k;
     if (fixedfieldtest(V))
     {
-      ulong av=avma;
+      gpmem_t av=avma;
       GEN s=fixedfieldpol(O,L,mod,S,deg);
       GEN P=FpV_roots_to_pol(s,mod,v);
       P=FpX_center(P,mod);
@@ -1567,7 +1376,7 @@ fixedfieldsurmer(GEN O, GEN L, GEN mod, GEN l, GEN p, 
 GEN
 fixedfieldsympol(GEN O, GEN L, GEN mod, GEN l, GEN p, GEN S, GEN deg, long v)
 {
-  ulong ltop=avma;
+  gpmem_t ltop=avma;
   GEN V=NULL;
   GEN LN=cgetg(lg(L),t_VEC);
   GEN Ll=FpV_red(L,l);
@@ -1624,7 +1433,7 @@ fixedfieldinclusion(GEN O, GEN PL)
 GEN
 vandermondeinversemod(GEN L, GEN T, GEN den, GEN mod)
 {
-  ulong   av;
+  gpmem_t av;
   int     i, j, n = lg(L);
   long    x = varn(T);
   GEN     M, P, Tp;
@@ -1651,7 +1460,8 @@ vandermondeinversemod(GEN L, GEN T, GEN den, GEN mod)
 static GEN
 vectopol(GEN v, GEN M, GEN den , GEN mod, long x)
 {
-  long n = lg(v),i,k,av;
+  long n = lg(v), i, k;
+  gpmem_t av;
   GEN z = cgetg(n+1,t_POL),p1,p2,mod2;
   av=avma;
   mod2=gclone(shifti(mod,-1));/*clone*/
@@ -1677,7 +1487,8 @@ vectopol(GEN v, GEN M, GEN den , GEN mod, long x)
 static GEN
 permtopol(GEN p, GEN L, GEN M, GEN den, GEN mod, long x)
 {
-  long n = lg(L),i,k,av;
+  long n = lg(L), i, k;
+  gpmem_t av;
   GEN z = cgetg(n+1,t_POL),p1,p2,mod2;
   if (lg(p) != n) err(talker,"incorrect permutation in permtopol");
   av=avma;
@@ -1714,116 +1525,7 @@ galoisgrouptopol( GEN res, GEN L, GEN M, GEN den, GEN 
   }
   return aut;
 }
-/* No more used and ugly.
- * However adding corepartial to GP seem a good idea.
- * 
- * factorise partiellement n sous la forme n=d*u*f^2 ou d est un
- * discriminant fondamental et u n'est pas un carre parfait et
- * retourne u*f 
- * Note: Parfois d n'est pas un discriminant (congru a 3 mod 4).
- * Cela se produit si u est congrue a 3 mod 4.
- */
-GEN
-corediscpartial(GEN n)
-{
-  ulong   av = avma;
-  long i,r,s;
-  GEN fa, p1, p2, p3, res = gun, res2 = gun, res3=gun;
-  /*d=res,f=res2,u=res3*/
-  if (gcmp1(n))
-    return gun;
-  fa = auxdecomp(n, 0);
-  p1 = (GEN) fa[1];
-  p2 = (GEN) fa[2];
-  for (i = 1; i < lg(p1) - 1; i++)
-  {
-    p3 = (GEN) p2[i];
-    if (mod2(p3))
-      res = mulii(res, (GEN) p1[i]);
-    if (!gcmp1(p3))
-      res2 = mulii(res2, gpow((GEN) p1[i], shifti(p3, -1), 0));
-  }
-  p3 = (GEN) p2[i];
-  if (mod2(p3))		/* impaire: verifions */
-  {
-    if (!gcmp1(p3))
-      res2 = mulii(res2, gpow((GEN) p1[i], shifti(p3, -1), 0));
-    if (isprime((GEN) p1[i]))
-      res = mulii(res, (GEN) p1[i]);
-    else
-      res3 = (GEN)p1[i];
-  }
-  else			/* paire:OK */
-    res2 = mulii(res2, gpow((GEN) p1[i], shifti(p3, -1), 0));
-  r = mod4(res);
-  if (signe(res) < 0)
-    r = 4 - r;
-  s = mod4(res3);/*res2 est >0*/
-  if (r == 3 && s!=3)
-    /*d est n'est pas un discriminant mais res2 l'est: corrige*/
-    res2 = gmul2n((GEN) res2, -1);
-  return gerepileupto(av,gmul(res2,res3));
-}
-GEN respm(GEN x,GEN y,GEN pm);
-GEN ZX_disc(GEN x);
 
-GEN
-indexpartial(GEN P, GEN DP)
-{
-  ulong   av = avma;
-  long i, nb;
-  GEN fa, p1, res = gun, dP;
-  dP = derivpol(P);
-  if(DEBUGLEVEL>=5) gentimer(3);
-  if(!DP) DP = ZX_disc(P);
-  DP = mpabs(DP);
-  if(DEBUGLEVEL>=5) genmsgtimer(3,"IndexPartial: discriminant");
-  fa = auxdecomp(DP, 0);
-  if(DEBUGLEVEL>=5) genmsgtimer(3,"IndexPartial: factorization");
-  nb = lg(fa[1]);
-  for (i = 1; i < nb; i++)
-  {
-    GEN p=gmael(fa,1,i);
-    GEN e=gmael(fa,2,i);
-    p1 = powgi(p,shifti(e,-1));
-    if ( i==nb-1 )
-    {
-      if ( mod2(e) && !isprime(p) )
-	p1 = mulii(p1,p);
-    }
-    else if ( cmpis(e,4)>=0 )
-    {
-      if(DEBUGLEVEL>=5) fprintferr("IndexPartial: factor %Z ",p1);
-      p1 = mppgcd(p1, respm(P,dP,p1));
-      if(DEBUGLEVEL>=5) 
-      {
-        fprintferr("--> %Z : ",p1);
-        genmsgtimer(3,"");
-      }
-    }
-    res=mulii(res,p1);
-  }
-  return gerepileupto(av,res);
-}
-
-/* Calcule la puissance exp d'une permutation decompose en cycle */
-GEN
-permcyclepow(GEN O, long exp)
-{
-  int     j, k, n;
-  GEN     p;
-  for (n = 0, j = 1; j < lg(O); j++)
-    n += lg(O[j]) - 1;
-  p = cgetg(n + 1, t_VECSMALL);
-  for (j = 1; j < lg(O); j++)
-  {
-    n = lg(O[j]) - 1;
-    for (k = 1; k <= n; k++)
-      p[mael(O,j,k)] = mael(O,j,1 + (k + exp - 1) % n);
-  }
-  return p;
-}
-
 /*
  * Casse l'orbite O en sous-orbite d'ordre premier correspondant a des
  * puissance de l'element
@@ -1831,22 +1533,19 @@ permcyclepow(GEN O, long exp)
 GEN
 splitorbite(GEN O)
 {
-  ulong   ltop = avma, lbot;
+  gpmem_t lbot, ltop = avma;
   int     i, n;
-  GEN     F, fc, res;
+  GEN     fc, res;
   n = lg(O[1]) - 1;
-  F = factor(stoi(n));
-  fc = cgetg(lg(F[1]), t_VECSMALL);
-  for (i = 1; i < lg(fc); i++)
-    fc[i] = itos(powgi(gmael(F,1,i), gmael(F,2,i)));
+  fc = decomp_primary_small(n);
   lbot = avma;
   res = cgetg(3, t_VEC);
   res[1] = lgetg(lg(fc), t_VEC);
   res[2] = lgetg(lg(fc), t_VECSMALL);
   for (i = 1; i < lg(fc); i++)
   {
-    ((GEN **) res)[1][lg(fc) - i] = permcyclepow(O, n / fc[i]);
-    ((GEN *) res)[2][lg(fc) - i] = fc[i];
+    mael(res,1,lg(fc) - i) = (long)cyc_powtoperm(O, n / fc[i]);
+    mael(res,2,lg(fc) - i) = fc[i];
   }
   if ( DEBUGLEVEL>=4 )
     fprintferr("GaloisConj:splitorbite: %Z\n",res);
@@ -1863,6 +1562,7 @@ splitorbite(GEN O)
  * modulo l ppp:  plus petit diviseur premier du degre de T. primepointer:
  * permet de calculer les premiers suivant p.
  */
+enum ga_code {ga_all_normal=1,ga_ext_2=2,ga_non_wss=4};
 struct galois_analysis
 {
   long    p;
@@ -1871,26 +1571,28 @@ struct galois_analysis
   long    l;
   long    ppp;
   long    p4;
-  enum {ga_all_normal=1,ga_ext_2=2,ga_non_wss=4} group;
+  enum ga_code group;
   byteptr primepointer;
 };
 void
 galoisanalysis(GEN T, struct galois_analysis *ga, long calcul_l, long karma_type)
 {
-  ulong ltop=avma;
+  gpmem_t ltop=avma;
   long n,p;
   long i;
-  enum {k_amoeba=0,k_snake=1,k_fish=2,k_bird=4,k_rodent=6,k_dog=8,k_human=9,k_cat=12} karma;
+  enum k_code {k_amoeba=0,k_snake=1,k_fish=2,k_bird=4,k_rodent=6,k_dog=8,k_human=9,k_cat=12} karma;
   long group,linf;
   /*TODO: complete the table to at least 200*/
   const int prim_nonss_orders[]={36,48,56,60,72,75,80,96,108,120,132,0};
   GEN F,Fp,Fe,Fpe,O;
-  long np;
+  long min_prime,np;
   long order,phi_order;
   long plift,nbmax,nbtest,deg;
   byteptr primepointer,pp;
-  if (DEBUGLEVEL >= 1)
-    timer2();
+
+  if (!ZX_is_squarefree(T))
+    err(talker, "Polynomial not squarefree in galoisinit");
+  if (DEBUGLEVEL >= 1) (void)timer2();
   n = degpol(T);
   if (!karma_type) karma_type=n;
   else err(warner,"entering black magic computation");
@@ -1939,6 +1641,7 @@ galoisanalysis(GEN T, struct galois_analysis *ga, long
       break;
   }
   /*Now, we study the orders of the Frobenius elements*/
+  min_prime=n*max((long)(BITS_IN_LONG-bfffo(n)-4),2);
   plift = 0;
   nbmax = 8+(n>>1);
   nbtest = 0; karma = k_amoeba;
@@ -1950,16 +1653,13 @@ galoisanalysis(GEN T, struct galois_analysis *ga, long
         || ((group&ga_non_wss) && order == Fp[np]))
          && (nbtest < 3 * nbmax || !(group&ga_non_wss)) ;)
   {
-    ulong   av;
-    long    prime_incr;
+    gpmem_t av;
     GEN     ip,FS,p1;
     long o,norm_o=1;
-    prime_incr = *primepointer++;
-    if (!prime_incr)
-      err(primer1);
-    p += prime_incr;
+
+    NEXT_PRIME_VIADIFF_CHECK(p,primepointer);
     /*discard small primes*/
-    if (p <= (n << 1))
+    if (p <= min_prime)
       continue;
     ip=stoi(p);
     if (!FpX_is_squarefree(T,ip))
@@ -2011,7 +1711,7 @@ galoisanalysis(GEN T, struct galois_analysis *ga, long
 	  deg = norm_o;
 	  order = o;
 	  plift = p;
-	  karma=cgcd(p-1,karma_type);
+	  karma=(enum k_code)cgcd(p-1,karma_type);
 	  pp = primepointer;
 	  group |= ga_all_normal;
 	}
@@ -2021,7 +1721,7 @@ galoisanalysis(GEN T, struct galois_analysis *ga, long
       {
 	order = o;
 	plift = p;
-	karma=cgcd(p-1,karma_type);
+	karma=(enum k_code)cgcd(p-1,karma_type);
 	pp = primepointer;
       }
     }
@@ -2042,18 +1742,15 @@ galoisanalysis(GEN T, struct galois_analysis *ga, long
   linf=n;
   if (calcul_l && O[1]<=linf)
   {
-    ulong   av;
-    long    prime_incr;
+    gpmem_t av;
     long    l=0;
     /*we need a totally splited prime l*/
     av = avma;
     while (l == 0)
     {
       long nb;
-      prime_incr = *primepointer++;
-      if (!prime_incr)
-	err(primer1);
-      p += prime_incr;
+
+      NEXT_PRIME_VIADIFF_CHECK(p,primepointer);
       if (p<=linf) continue;
       nb=FpX_nbroots(T,stoi(p));
       if (nb == n)
@@ -2072,7 +1769,7 @@ galoisanalysis(GEN T, struct galois_analysis *ga, long
     O[1]=l;
   }  
   ga->p = plift;
-  ga->group = group;
+  ga->group = (enum ga_code)group;
   ga->deg = deg;
   ga->ord = order;
   ga->l = O[1];
@@ -2091,7 +1788,7 @@ galoisanalysis(GEN T, struct galois_analysis *ga, long
 GEN
 a4galoisgen(GEN T, struct galois_test *td)
 {
-  ulong   ltop = avma, av, av2;
+  gpmem_t ltop = avma, av, av2;
   long    i, j, k;
   long    n;
   long    N, hop = 0;
@@ -2354,7 +2051,7 @@ a4galoisgen(GEN T, struct galois_test *td)
   orb[2] = (long) pfu;
   if (DEBUGLEVEL >= 4)
     fprintferr("A4GaloisConj:orb=%Z \n", orb);
-  O = (long**)permorbite(orb);
+  O = (long**) vecperm_orbits(orb, 12);
   if (DEBUGLEVEL >= 4)
     fprintferr("A4GaloisConj:O=%Z \n", O);
   av2 = avma;
@@ -2437,11 +2134,10 @@ s4makelift(GEN u, struct galois_lift *gl, GEN liftpow)
 static long
 s4test(GEN u, GEN liftpow, struct galois_lift *gl, GEN phi)
 {
-  ulong ltop=avma;
+  gpmem_t ltop=avma;
   GEN     res;
   int     bl,i,d = lgef(u)-2;
-  if (DEBUGLEVEL >= 6)
-    timer2();
+  if (DEBUGLEVEL >= 6) (void)timer2();
   if ( !d ) return 0;
   res=(GEN)u[2];
   for (i = 1; i < d; i++)
@@ -2473,7 +2169,7 @@ s4test(GEN u, GEN liftpow, struct galois_lift *gl, GEN
 static GEN
 s4releveauto(GEN misom,GEN Tmod,GEN Tp,GEN p,long a1,long a2,long a3,long a4,long a5,long a6)
 {
-  ulong ltop=avma;
+  gpmem_t ltop=avma;
   GEN u1,u2,u3,u4,u5;
   GEN pu1,pu2,pu3,pu4;
   pu1=FpX_mul( (GEN) Tmod[a2], (GEN) Tmod[a1],p);
@@ -2494,7 +2190,7 @@ GEN
 s4galoisgen(struct galois_lift *gl)
 {
   struct galois_testlift gt;
-  ulong   ltop = avma, av, ltop2;
+  gpmem_t av, ltop2, ltop = avma;
   GEN     Tmod, isom, isominv, misom;
   int     i, j;
   GEN     sg;
@@ -2554,10 +2250,9 @@ s4galoisgen(struct galois_lift *gl)
   av = avma;
   for (i = 0; i < 3; i++)
   {
-    ulong   av2;
+    gpmem_t av2, avm1, avm2;
     GEN     u;
     int     j1, j2, j3;
-    ulong   avm1, avm2;
     GEN     u1, u2, u3;
     if (i)
     {
@@ -2625,7 +2320,7 @@ suites4:
   avma = av;
   for (j = 1; j <= 3; j++)
   {
-    ulong   av2;
+    gpmem_t av2;
     GEN     u;
     int     w, l;
     int     z;
@@ -2639,8 +2334,8 @@ suites4:
       av2 = avma;
       for (w = 0; w < 4; w += 2)
       {
+	gpmem_t av3;
 	GEN     uu;
-	long    av3;
 	pj[6] = (w + pj[3]) & 3;
 	uu =FpX_add(FpXQ_mul((GEN) bezoutcoeff[sg[5]],
 			  (GEN) pauto[(pj[6] & 3) + 1],TQ,Q),
@@ -2683,7 +2378,7 @@ suites4_2:
   {
     int     abc, abcdef;
     GEN     u;
-    ulong   av2;
+    gpmem_t av2;
     abc = (pj[1] + pj[2] + pj[3]) & 3;
     abcdef = (((abc + pj[4] + pj[5] - pj[6]) & 3) >> 1);
     u = s4releveauto(misom,Tmod,Tp,p,sg[1],sg[4],sg[2],sg[5],sg[3],sg[6]);
@@ -2744,19 +2439,19 @@ struct galois_frobenius
 static GEN
 galoisfindgroups(GEN lo, GEN sg, long f)
 {
-  ulong ltop=avma;
+  gpmem_t ltop=avma;
   GEN V,W;
   long i,j,k;
   V=cgetg(lg(lo),t_VEC);
   for(j=1,i=1;i<lg(lo);i++)
   {
-    ulong av=avma;
+    gpmem_t av=avma;
     GEN U;
     W=cgetg(lg(lo[i]),t_VECSMALL);
     for(k=1;k<lg(lo[i]);k++)
       W[k]=mael(lo,i,k)%f;
-    sortvecsmall(W); 
-    U=uniqvecsmall(W);
+    vecsmall_sort(W); 
+    U=vecsmall_uniq(W);
     if (gegal(U, sg))
     {
       cgiv(U);
@@ -2773,24 +2468,14 @@ galoisfindgroups(GEN lo, GEN sg, long f)
 static long
 galoisfrobeniustest(GEN aut, struct galois_lift *gl, GEN frob)
 {
-  ulong ltop=avma;
+  gpmem_t ltop = avma;
   GEN tlift = FpX_center(FpX_Fp_mul(aut,gl->den,gl->Q), gl->Q);
-  long res=poltopermtest(tlift, gl, frob);
-  avma=ltop;
+  long res = poltopermtest(tlift, gl, frob);
+  avma = ltop;
   return res;
 }
 
 static GEN
-galoismakepsiab(long g)
-{
-  GEN psi=cgetg(g+1,t_VECSMALL);
-  long i;
-  for (i = 1; i <= g; i++)
-    psi[i] = 1;
-  return psi;
-}
-
-static GEN
 galoismakepsi(long g, GEN sg, GEN pf)
 {
   GEN psi=cgetg(g+1,t_VECSMALL);
@@ -2802,10 +2487,10 @@ galoismakepsi(long g, GEN sg, GEN pf)
 }
 
 static GEN
-galoisfrobeniuslift(GEN T, GEN den, GEN L,  GEN Lden, long gmask,
+galoisfrobeniuslift(GEN T, GEN den, GEN L,  GEN Lden, 
     struct galois_frobenius *gf,  struct galois_borne *gb) 
 {
-  ulong ltop=avma,av2;
+  gpmem_t ltop=avma, av2;
   struct galois_testlift gt;
   struct galois_lift gl;
   GEN res;
@@ -2816,7 +2501,7 @@ galoisfrobeniuslift(GEN T, GEN den, GEN L,  GEN Lden, 
   if (DEBUGLEVEL >= 4)
     fprintferr("GaloisConj:p=%ld deg=%ld fp=%ld\n", gf->p, deg, gf->fp);
   res = cgetg(lg(L), t_VECSMALL);
-  gf->psi = galoismakepsiab(g);
+  gf->psi = vecsmall_const(g,1);
   av2=avma;
   initlift(T, den, ip, L, Lden, gb, &gl);
   aut = galoisdolift(&gl, res);
@@ -2842,11 +2527,11 @@ galoisfrobeniuslift(GEN T, GEN den, GEN L,  GEN Lden, 
   frob = cgetg(lg(L), t_VECSMALL);
   for(k=lg(Fp)-1;k>=1;k--)
   {
-    long btop=avma;
+    gpmem_t btop=avma;
     GEN psi=NULL,fres=NULL,sg;
     long el=gf->fp, dg=1, dgf=1;
     long e,pr;
-    sg=permidentity(1);
+    sg=perm_identity(1);
     for(e=1;e<=Fe[k];e++)
     {
       long l;
@@ -2858,7 +2543,7 @@ galoisfrobeniuslift(GEN T, GEN den, GEN L,  GEN Lden, 
       if (galoisfrobeniustest((GEN)gt.pauto[el+1],&gl,frob))
       {
 	dgf = dg; 
-	psi = galoismakepsiab(g);
+	psi = vecsmall_const(g,1);
 	fres= gcopy(frob);
 	continue;
       }
@@ -2891,7 +2576,7 @@ galoisfrobeniuslift(GEN T, GEN den, GEN L,  GEN Lden, 
     } 
     else
     {
-      GEN cp=permapply(res,fres);
+      GEN cp=perm_mul(res,fres);
       for(i=1;i<lg(res);i++) res[i]=cp[i];
       for(i=1;i<lg(psi);i++) gf->psi[i]=(dgf*gf->psi[i]+deg*psi[i])%pr;
     }
@@ -2913,7 +2598,7 @@ galoisfrobeniuslift(GEN T, GEN den, GEN L,  GEN Lden, 
   {
     /*We need to normalise result so that psi[g]=1*/
     long im=itos(mpinvmod(stoi(gf->psi[g]),stoi(gf->fp)));
-    GEN cp=permcyclepow(permorbite(res), im);
+    GEN cp=perm_pow(res, im);
     for(i=1;i<lg(res);i++) res[i]=cp[i];
     for(i=1;i<lg(gf->psi);i++) gf->psi[i]=mulssmod(im,gf->psi[i],gf->fp);
     avma=av2;
@@ -2923,11 +2608,11 @@ galoisfrobeniuslift(GEN T, GEN den, GEN L,  GEN Lden, 
 }
 
 static GEN
-galoisfindfrobenius(GEN T, GEN L, GEN M, GEN den, struct galois_frobenius *gf,
+galoisfindfrobenius(GEN T, GEN L, GEN den, struct galois_frobenius *gf,
     struct galois_borne *gb, const struct galois_analysis *ga)
 {
-  ulong ltop=avma,lbot;
-  long try=0;
+  gpmem_t lbot, ltop=avma;
+  long Try=0;
   long n = degpol(T), deg, gmask;
   byteptr primepointer = ga->primepointer;
   GEN Lden,frob;
@@ -2936,9 +2621,9 @@ galoisfindfrobenius(GEN T, GEN L, GEN M, GEN den, stru
   gmask=(ga->group&ga_ext_2)?3:1;
   for (;;)
   {
-    ulong   av = avma;
+    gpmem_t av = avma;
     long    isram;
-    long    i,c;
+    long    i;
     GEN     ip,Tmod;
     ip = stoi(gf->p);
     Tmod = lift((GEN) factmod(T, ip));
@@ -2959,7 +2644,7 @@ galoisfindfrobenius(GEN T, GEN L, GEN M, GEN den, stru
       gf->Tmod=gcopy((GEN)Tmod[1]);
       if ( ((gmask&1) && gf->fp % deg == 0) || ((gmask&2) && gf->fp % 2== 0) )
       {
-	frob=galoisfrobeniuslift(T, den, L, Lden, gmask, gf, gb);
+	frob=galoisfrobeniuslift(T, den, L, Lden, gf, gb);
 	if (frob)
 	{
 	  GEN *gptr[3];
@@ -2978,15 +2663,12 @@ galoisfindfrobenius(GEN T, GEN L, GEN M, GEN den, stru
 	  avma = ltop;
 	  return NULL;
 	}
-	try++;
-	if ( (ga->group&ga_non_wss) && try > n )
+	Try++;
+	if ( (ga->group&ga_non_wss) && Try > n )
 	  err(warner, "galoisconj _may_ hang up for this polynomial");
       }
     }
-    c = *primepointer++;
-    if (!c)
-      err(primer1);
-    gf->p += c;
+    NEXT_PRIME_VIADIFF_CHECK(gf->p, primepointer);
     if (DEBUGLEVEL >= 4)
       fprintferr("GaloisConj:next p=%ld\n", gf->p);
     avma = av;
@@ -2995,20 +2677,99 @@ galoisfindfrobenius(GEN T, GEN L, GEN M, GEN den, stru
 
 GEN
 galoisgen(GEN T, GEN L, GEN M, GEN den, struct galois_borne *gb,
+	  const struct galois_analysis *ga);
+static GEN
+galoisgenfixedfield(GEN Tp, GEN Pmod, GEN V, GEN ip, struct galois_borne *gb, GEN Pg)
+{
+  gpmem_t ltop=avma;
+  GEN     P, PL, Pden, PM, Pp, Pladicabs;
+  GEN     tau, PG;
+  long    g,gp;
+  long    x=varn(Tp);
+  P=(GEN)V[2];
+  PL=(GEN)V[1];
+  gp=lg(Pmod)-1;
+  Pp = FpX_red(P,ip);
+  if (DEBUGLEVEL>=6)
+    fprintferr("GaloisConj: Fixed field %Z\n",P);
+  if (degpol(P)==2)
+  {
+    PG=cgetg(3,t_VEC);
+    PG[1]=lgetg(2,t_VEC);
+    PG[2]=lgetg(2,t_VECSMALL);
+    mael(PG,1,1)=lgetg(3,t_VECSMALL);
+    mael(PG,2,1)=2;
+    mael3(PG,1,1,1)=2;
+    mael3(PG,1,1,2)=1;
+    tau=deg1pol(stoi(-1),negi((GEN)P[3]),x);
+    tau = lift(gmul(tau,gmodulcp(gun,ip)));
+    tau = FpX_FpXQ_compo((GEN) Pmod[gp], tau,Pp,ip);
+    tau = FpX_gcd(Pp, tau,ip);
+    tau = FpX_Fp_mul(tau,mpinvmod((GEN) tau[lgef(tau) - 1],ip),ip);
+    for (g = 1; g <= gp; g++)
+      if (gegal(tau, (GEN) Pmod[g]))
+	break;
+    if (g == lg(Pmod))
+      return NULL;
+    Pg[1]=g;
+  }
+  else
+  {
+    struct galois_analysis Pga;
+    struct galois_borne Pgb;
+    long j;
+    galoisanalysis(P, &Pga, 0, 0);
+    if (Pga.deg == 0)
+      return NULL;		/* Avoid computing the discriminant */
+    Pgb.l = gb->l;
+    Pden = galoisborne(P, NULL, &Pgb, Pga.ppp);
+    Pladicabs=Pgb.ladicabs;
+    if (Pgb.valabs > gb->valabs)
+    {
+      if (DEBUGLEVEL>=4)
+	fprintferr("GaloisConj:increase prec of p-adic roots of %ld.\n"
+	    ,Pgb.valabs-gb->valabs);
+      PL = rootpadicliftroots(P,PL,gb->l,Pgb.valabs);
+    }
+    PM = vandermondeinversemod(PL, P, Pden, Pgb.ladicabs);
+    PG = galoisgen(P, PL, PM, Pden, &Pgb, &Pga);
+    if (PG == gzero) return NULL;
+    for (j = 1; j < lg(PG[1]); j++)
+    {
+      gpmem_t btop=avma;
+      tau = permtopol(gmael(PG,1,j), PL, PM, Pden, Pladicabs, x);
+      tau = lift(gmul(tau,gmodulcp(gun,ip)));
+      tau = FpX_FpXQ_compo((GEN) Pmod[gp], tau,Pp,ip);
+      tau = FpX_gcd(Pp, tau,ip);
+      tau = FpX_Fp_mul(tau,mpinvmod((GEN) tau[lgef(tau) - 1],ip),ip);
+      for (g = 1; g < lg(Pmod); g++)
+	if (gegal(tau, (GEN) Pmod[g]))
+	  break;
+      if (g == lg(Pmod))
+	return NULL;
+      avma=btop;
+      Pg[j]=g;
+    }
+  }
+  return gerepilecopy(ltop,PG);
+}
+
+GEN
+galoisgen(GEN T, GEN L, GEN M, GEN den, struct galois_borne *gb,
 	  const struct galois_analysis *ga)
 {
   struct galois_test td;
   struct galois_frobenius gf;
-  ulong   ltop = avma, lbot, ltop2;
+  gpmem_t lbot, ltop2, ltop = avma;
   long    n, p, deg;
-  long    fp,gp;
+  long    fp;
   long    x;
   int     i, j;
   GEN     Lden, sigma;
   GEN     Tmod, res, pf = gzero, split, ip;
   GEN     frob;
   GEN     O;
-  GEN     P, PG, PL, Pden, PM, Pmod, Pp, Pladicabs;
+  GEN     PG, Pg;
   n = degpol(T);
   if (!ga->deg)
     return gzero;
@@ -3017,7 +2778,7 @@ galoisgen(GEN T, GEN L, GEN M, GEN den, struct galois_
     fprintferr("GaloisConj:denominator:%Z\n", den);
   if (n == 12 && ga->ord==3)	/* A4 is very probable,so test it first */
   {
-    long    av = avma;
+    gpmem_t av = avma;
     if (DEBUGLEVEL >= 4)
       fprintferr("GaloisConj:Testing A4 first\n");
     inittest(L, M, gb->bornesol, gb->ladicsol, &td);
@@ -3030,7 +2791,7 @@ galoisgen(GEN T, GEN L, GEN M, GEN den, struct galois_
   }
   if (n == 24 && ga->ord==3)	/* S4 is very probable,so test it first */
   {
-    long    av = avma;
+    gpmem_t av = avma;
     struct galois_lift gl;
     if (DEBUGLEVEL >= 4)
       fprintferr("GaloisConj:Testing S4 first\n");
@@ -3042,17 +2803,17 @@ galoisgen(GEN T, GEN L, GEN M, GEN den, struct galois_
       return gerepile(ltop, lbot, PG);
     avma = av;
   }
-  frob=galoisfindfrobenius(T, L, M, den, &gf, gb, ga);
+  frob=galoisfindfrobenius(T, L, den, &gf, gb, ga);
   if (!frob)
   {
     ltop=avma;
     return gzero;
   }
-  p=gf.p;deg=gf.deg;
+  p=gf.p;
   ip=stoi(p);
   Tmod=gf.Tmod;
-  fp=gf.fp; gp=n/fp;
-  O = permorbite(frob);
+  fp=gf.fp;
+  O = perm_cycles(frob);
   split = splitorbite(O);
   deg=lg(O[1])-1;
   sigma = permtopol(frob, L, M, den, gb->ladicabs, x);
@@ -3073,67 +2834,26 @@ galoisgen(GEN T, GEN L, GEN M, GEN den, struct galois_
   }
   if (DEBUGLEVEL >= 9)
     fprintferr("GaloisConj:Frobenius:%Z\n", sigma);
+  Pg=cgetg(lg(O),t_VECSMALL);
   {
-    GEN     V, Tp, Sp, sym, dg;
+    gpmem_t btop=avma;
+    GEN     V, sym, dg, Tp, Sp, Pmod;
     sym=cgetg(lg(L),t_VECSMALL);
     dg=cgetg(lg(L),t_VECSMALL);
     V = fixedfieldsympol(O, L, gb->ladicabs, gb->l, ip, sym, dg, x);
-    P=(GEN)V[2];
-    PL=(GEN)V[1];
     Tp = FpX_red(T,ip);
     Sp = fixedfieldpolsigma(sigma,ip,Tp,sym,dg,deg);
     Pmod = fixedfieldfactmod(Sp,ip,Tmod);
-    Pp = FpX_red(P,ip);
-    if (DEBUGLEVEL >= 4)
+    PG=galoisgenfixedfield(Tp, Pmod, V, ip, gb, Pg);
+    if (PG == NULL)
     {
-      fprintferr("GaloisConj:P=%Z   \n", P);
-      fprintferr("GaloisConj:psi=%Z\n", gf.psi);
-      fprintferr("GaloisConj:Sp=%Z\n", Sp);
-      fprintferr("GaloisConj:Pmod=%Z\n", Pmod);
-      fprintferr("GaloisConj:Tmod=%Z\n", Tmod);
+      avma = ltop;
+      return gzero;
     }
-    if ( n == (deg<<1))
-    {
-      PG=cgetg(3,t_VEC);
-      PG[1]=lgetg(2,t_VEC);
-      PG[2]=lgetg(2,t_VECSMALL);
-      mael(PG,1,1)=lgetg(3,t_VECSMALL);
-      mael(PG,2,1)=2;
-      mael3(PG,1,1,1)=2;
-      mael3(PG,1,1,2)=1;
-      Pden=NULL; PM=NULL; Pladicabs=NULL;/*make KB happy*/
-    }
-    else
-    {
-      struct galois_analysis Pga;
-      struct galois_borne Pgb;
-      galoisanalysis(P, &Pga, 0, 0);
-      if (Pga.deg == 0)
-      {
-	avma = ltop;
-	return gzero;		/* Avoid computing the discriminant */
-      }
-      Pgb.l = gb->l;
-      Pden = galoisborne(P, NULL, &Pgb, Pga.ppp);
-      Pladicabs=Pgb.ladicabs;
-      if (Pgb.valabs > gb->valabs)
-      {
-	if (DEBUGLEVEL>=4)
-	  fprintferr("GaloisConj:increase prec of p-adic roots of %ld.\n"
-	      ,Pgb.valabs-gb->valabs);
-	PL = rootpadicliftroots(P,PL,gb->l,Pgb.valabs);
-      }
-      PM = vandermondeinversemod(PL, P, Pden, Pgb.ladicabs);
-      PG = galoisgen(P, PL, PM, Pden, &Pgb, &Pga);
-    }
+    if (DEBUGLEVEL >= 4)
+      fprintferr("GaloisConj:Back to Earth:%Z\n", PG);
+    PG=gerepileupto(btop, PG);
   }
-  if (DEBUGLEVEL >= 4)
-    fprintferr("GaloisConj:Retour sur Terre:%Z\n", PG);
-  if (PG == gzero)
-  {
-    avma = ltop;
-    return gzero;
-  }
   inittest(L, M, gb->bornesol, gb->ladicsol, &td);
   lbot = avma;
   res = cgetg(3, t_VEC);
@@ -3149,41 +2869,16 @@ galoisgen(GEN T, GEN L, GEN M, GEN den, struct galois_
   ltop2 = avma;
   for (j = 1; j < lg(PG[1]); j++)
   {
-    GEN     B, tau;
-    long    t, g;
+    GEN     B;
+    long    t;
     long    w, sr, dss;
     if (DEBUGLEVEL >= 6)
       fprintferr("GaloisConj:G[%d]=%Z  d'ordre relatif %d\n", j,
 	  ((GEN **) PG)[1][j], ((long **) PG)[2][j]);
-    B = permorbite(((GEN **) PG)[1][j]);
+    B = perm_cycles(gmael(PG,1,j));
     if (DEBUGLEVEL >= 6)
       fprintferr("GaloisConj:B=%Z\n", B);
-    if (n == (deg<<1))
-      tau=deg1pol(stoi(-1),negi((GEN)P[3]),x);
-    else
-      tau = permtopol(gmael(PG,1,j), PL, PM, Pden, Pladicabs, x);
-    if (DEBUGLEVEL >= 9)
-      fprintferr("GaloisConj:Paut=%Z\n", tau);
-    /*tau not in Z[X]*/
-    tau = lift(gmul(tau,gmodulcp(gun,ip)));
-    tau = FpX_FpXQ_compo((GEN) Pmod[gp], tau,Pp,ip);
-    tau = FpX_gcd(Pp, tau,ip);
-    /*normalize gcd*/
-    tau = FpX_Fp_mul(tau,mpinvmod((GEN) tau[lgef(tau) - 1],ip),ip);
-    if (DEBUGLEVEL >= 6)
-      fprintferr("GaloisConj:tau=%Z\n", tau);
-    for (g = 1; g < lg(Pmod); g++)
-      if (gegal(tau, (GEN) Pmod[g]))
-	break;
-    if (DEBUGLEVEL >= 6)
-      fprintferr("GaloisConj:g=%ld\n", g);
-    if (g == lg(Pmod))
-    {
-      freetest(&td);
-      avma = ltop;
-      return gzero;
-    }
-    w = gf.psi[g];
+    w = gf.psi[Pg[j]];
     dss = deg / cgcd(deg, w - 1);
     sr = 1;
     for (i = 1; i < lg(B[1]) - 1; i++)
@@ -3226,7 +2921,7 @@ galoisgen(GEN T, GEN L, GEN M, GEN den, struct galois_
 GEN
 galoisconj4(GEN T, GEN den, long flag, long karma)
 {
-  ulong   ltop = avma;
+  gpmem_t ltop = avma;
   GEN     G, L, M, res, aut, grp=NULL;/*keep gcc happy on the wall*/
   struct galois_analysis ga;
   struct galois_borne gb;
@@ -3234,8 +2929,7 @@ galoisconj4(GEN T, GEN den, long flag, long karma)
   if (typ(T) != t_POL)
   {
     T = checknf(T);
-    if (!den)
-      den = gcoeff((GEN) T[8], degpol(T[1]), degpol(T[1]));
+    if (!den) den = Q_denom((GEN)T[7]);
     T = (GEN) T[1];
   }
   n = degpol(T);
@@ -3274,8 +2968,7 @@ galoisconj4(GEN T, GEN den, long flag, long karma)
     den = absi(den);
   }
   gb.l = stoi(ga.l);
-  if (DEBUGLEVEL >= 1)
-    timer2();
+  if (DEBUGLEVEL >= 1) (void)timer2();
   den = galoisborne(T, den, &gb, ga.ppp);
   if (DEBUGLEVEL >= 1)
     msgtimer("galoisborne()");
@@ -3300,8 +2993,7 @@ galoisconj4(GEN T, GEN den, long flag, long karma)
     avma = ltop;
     return gzero;
   }
-  if (DEBUGLEVEL >= 1)
-    timer2();
+  if (DEBUGLEVEL >= 1) (void)timer2();
   if (flag)
   {
     grp = cgetg(9, t_VEC);
@@ -3317,13 +3009,13 @@ galoisconj4(GEN T, GEN den, long flag, long karma)
     grp[8] = (long) gcopy((GEN) G[2]);
   }
   res = cgetg(n + 1, t_VEC);
-  res[1] = (long) permidentity(n);
+  res[1] = (long) perm_identity(n);
   k = 1;
   for (i = 1; i < lg(G[1]); i++)
   {
     int     c = k * (((long **) G)[2][i] - 1);
     for (j = 1; j <= c; j++)	/* I like it */
-      res[++k] = (long) permapply((GEN) res[j], ((GEN **) G)[1][i]);
+      res[++k] = (long) perm_mul((GEN) res[j], ((GEN **) G)[1][i]);
   }
   if (flag)
   {
@@ -3333,7 +3025,7 @@ galoisconj4(GEN T, GEN den, long flag, long karma)
   aut = galoisgrouptopol(res,L,M,den,gb.ladicsol, varn(T));
   if (DEBUGLEVEL >= 1)
     msgtimer("Calcul polynomes");
-  return gerepileupto(ltop, aut);
+  return gerepileupto(ltop, gen_sort(aut, 0, cmp_pol));
 }
 
 /* Calcule le nombre d'automorphisme de T avec forte probabilité */
@@ -3341,7 +3033,7 @@ galoisconj4(GEN T, GEN den, long flag, long karma)
 long
 numberofconjugates(GEN T, long pdepart)
 {
-  ulong   ltop = avma, ltop2;
+  gpmem_t ltop2, ltop = avma;
   long    n, p, nbmax, nbtest;
   long    card;
   byteptr primepointer;
@@ -3357,13 +3049,11 @@ numberofconjugates(GEN T, long pdepart)
   ltop2 = avma;
   for (p = 0, primepointer = diffptr; nbtest < nbmax && card > 1;)
   {
-    long    s, c;
+    long    s;
     long    isram;
     GEN     S;
-    c = *primepointer++;
-    if (!c)
-      err(primer1);
-    p += c;
+
+    NEXT_PRIME_VIADIFF_CHECK(p,primepointer);
     if (p < pdepart)
       continue;
     S = (GEN) simplefactmod(T, stoi(p));
@@ -3397,7 +3087,7 @@ numberofconjugates(GEN T, long pdepart)
 GEN
 galoisconj0(GEN nf, long flag, GEN d, long prec)
 {
-  ulong   ltop;
+  gpmem_t ltop;
   GEN     G, T;
   long    card;
   if (typ(nf) != t_POL)
@@ -3457,10 +3147,10 @@ galoisconj0(GEN nf, long flag, GEN d, long prec)
 long
 isomborne(GEN P, GEN den, GEN p)
 {
-  ulong ltop=avma;
+  gpmem_t ltop=avma;
   struct galois_borne gb;
   gb.l=p;
-  galoisborne(P,den,&gb,degpol(P));
+  (void)galoisborne(P,den,&gb,degpol(P));
   avma=ltop;
   return gb.valsol;
 }
@@ -3522,7 +3212,7 @@ galoiscosets(GEN O, GEN perm)
   long i,j,k,u;
   long o = lg(O)-1, f = lg(O[1])-1;
   GEN C = cgetg(lg(O),t_VECSMALL);
-  ulong av=avma;
+  gpmem_t av=avma;
   GEN RC=cgetg(lg(perm),t_VECSMALL);
   for(i=1;i<lg(RC);i++)
     RC[i]=0;
@@ -3543,7 +3233,7 @@ GEN
 fixedfieldfactor(GEN L, GEN O, GEN perm, GEN M, GEN den, GEN mod,
                  long x,long y)
 {
-  ulong   ltop=avma;
+  gpmem_t ltop=avma;
   GEN     F,G,V,res,cosets;
   int     i, j, k;
   F=cgetg(lg(O[1])+1,t_COL);
@@ -3564,7 +3254,7 @@ fixedfieldfactor(GEN L, GEN O, GEN perm, GEN M, GEN de
   res=cgetg(lg(O),t_VEC);
   for (i = 1; i < lg(O); i++)
   {
-    ulong av=avma;
+    gpmem_t av=avma;
     long ii,jj;
     GEN G=cgetg(lg(O),t_VEC);
     for (ii = 1; ii < lg(O); ii++)
@@ -3588,7 +3278,7 @@ fixedfieldfactor(GEN L, GEN O, GEN perm, GEN M, GEN de
 GEN
 galoisfixedfield(GEN gal, GEN perm, long flag, long y)
 {
-  ulong   ltop = avma, lbot;
+  gpmem_t lbot, ltop = avma;
   GEN     L, P, S, PL, O, res, mod;
   long    x, n;
   int     i;
@@ -3600,16 +3290,18 @@ galoisfixedfield(GEN gal, GEN perm, long flag, long y)
     err(flagerr, "galoisfixedfield");
   if (typ(perm) == t_VEC)
   {
-    if (lg(perm)==1)
-      perm=permidentity(n);
-    else
-      for (i = 1; i < lg(perm); i++)
-	if (typ(perm[i]) != t_VECSMALL || lg(perm[i])!=n+1)
-	  err(typeer, "galoisfixedfield");
+    for (i = 1; i < lg(perm); i++)
+      if (typ(perm[i]) != t_VECSMALL || lg(perm[i])!=n+1)
+        err(typeer, "galoisfixedfield");
+    O = vecperm_orbits(perm, n);
   }
   else if (typ(perm) != t_VECSMALL || lg(perm)!=n+1 )
+  {
     err(typeer, "galoisfixedfield");
-  O = permorbite(perm);
+    return NULL; /* not reached */
+  }
+  else
+    O = perm_cycles(perm);
   {
     GEN sym=cgetg(lg(L),t_VECSMALL);
     GEN dg=cgetg(lg(L),t_VECSMALL);
@@ -3640,12 +3332,12 @@ galoisfixedfield(GEN gal, GEN perm, long flag, long y)
       Pden = galoisborne(P, NULL, &Pgb, 1);
       if (Pgb.valabs > val)
       {
-	if (DEBUGLEVEL>=4)
-	  fprintferr("GaloisConj:increase prec of p-adic roots of %ld.\n"
-	      ,Pgb.valabs-val);
-	PL = rootpadicliftroots(P,PL,Pgb.l,Pgb.valabs);
-	L = rootpadicliftroots((GEN) gal[1],L,Pgb.l,Pgb.valabs);
-	mod = Pgb.ladicabs;
+        if (DEBUGLEVEL>=4)
+          fprintferr("GaloisConj:increase prec of p-adic roots of %ld.\n"
+              ,Pgb.valabs-val);
+        PL = rootpadicliftroots(P,PL,Pgb.l,Pgb.valabs);
+        L = rootpadicliftroots((GEN) gal[1],L,Pgb.l,Pgb.valabs);
+        mod = Pgb.ladicabs;
       }
     }
     PM = vandermondeinversemod(PL, P, Pden, mod);
@@ -3658,7 +3350,7 @@ galoisfixedfield(GEN gal, GEN perm, long flag, long y)
     res[1] = (long) gcopy(P);
     res[2] = (long) gmodulcp(S, (GEN) gal[1]);
     res[3] = (long) fixedfieldfactor(L,O,(GEN)gal[6],
-	PM,Pden,mod,x,y);
+        PM,Pden,mod,x,y);
     return gerepile(ltop, lbot, res);
   }
 }
@@ -3667,39 +3359,38 @@ galoisfixedfield(GEN gal, GEN perm, long flag, long y)
 long 
 galoistestabelian(GEN gal)
 {
-  ulong ltop=avma;
-  long i,j;
+  gpmem_t ltop=avma;
+  GEN G;
+  long test;
   gal = checkgal(gal);
-  for(i=2;i<lg(gal[7]);i++)
-    for(j=1;j<i;j++)
-    {
-      long test=gegal(permapply(gmael(gal,7,i),gmael(gal,7,j)),
-	  permapply(gmael(gal,7,j),gmael(gal,7,i)));
-      avma=ltop;
-      if (!test) return 0;
-    }
-  return 1;  
+  G=cgetg(3,t_VEC);
+  G[1]=gal[7];
+  G[2]=gal[8];
+  test=group_isabelian(G);
+  avma=ltop;
+  return test;
 }
 
 GEN galoisisabelian(GEN gal, long flag)
 {
   long i, j;
-  long test, n=lg(gal[7]);
+  long test, n;
   GEN M;
   gal = checkgal(gal);
   test=galoistestabelian(gal);
   if (!test) return gzero;
   if (flag==1)  return gun;
   if (flag) err(flagerr,"galoisisabelian");
+  n=lg(gal[7]);
   M=cgetg(n,t_MAT);
   for(i=1;i<n;i++)
   {
-    ulong btop;
+    gpmem_t btop;
     GEN P;
     long k;
     M[i]=lgetg(n,t_COL);
     btop=avma;
-    P=permcyclepow(permorbite(gmael(gal,7,i)),mael(gal,8,i));
+    P=perm_pow(gmael(gal,7,i),mael(gal,8,i));
     for(j=1;j<lg(gal[6]);j++)
       if (gegal(P,gmael(gal,6,j)))
 	  break;
@@ -3717,295 +3408,26 @@ GEN galoisisabelian(GEN gal, long flag)
   }
   return M;
 }
-/* Calcule les orbites d'un sous-groupe de Z/nZ donne par un
- * generateur ou d'un ensemble de generateur donne par un vecteur. 
- */
-GEN
-subgroupcoset(long n, GEN v)
+GEN galoissubgroups(GEN G)
 {
-  ulong   ltop = avma, lbot;
-  int     i, j, k = 1, l, m, o, p, flag;
-  GEN     bit, cycle, cy;
-  cycle = cgetg(n, t_VEC);
-  bit = cgetg(n, t_VECSMALL);
-  for (i = 1; i < n; i++)
-    if (cgcd(i,n)==1)
-      bit[i] = 0;
-    else
-    {
-      bit[i] = -1;
-      k++;
-    }
-  for (l = 1; k < n;)
-  {
-    for (j = 1; bit[j]; j++);
-    cy = cgetg(n, t_VECSMALL);
-    m = 1;
-    cy[m++] = j;
-    bit[j] = 1;
-    k++;
-    do
-    {
-      flag = 0;
-      for (o = 1; o < lg(v); o++)
-      {
-	for (p = 1; p < m; p++)	/* m varie! */
-	{
-	  j = mulssmod(v[o],cy[p],n);
-	  if (!bit[j])
-	  {
-	    flag = 1;
-	    bit[j] = 1;
-	    k++;
-	    cy[m++] = j;
-	  }
-	}
-      }
-    }
-    while (flag);
-    setlg(cy, m);
-    cycle[l++] = (long) cy;
-  }
-  setlg(cycle, l);
-  lbot = avma;
-  cycle = gcopy(cycle);
-  return gerepile(ltop, lbot, cycle);
+  gpmem_t ltop=avma;
+  GEN H;
+  G = checkgal(G);
+  H=cgetg(3,t_VEC);
+  H[1]=G[7];
+  H[2]=G[8];
+  return gerepileupto(ltop,group_subgroups(H));
 }
-/* Calcule les elements d'un sous-groupe H de Z/nZ donne par un
- * generateur ou d'un ensemble de generateur donne par un vecteur (v). 
- *
- * cy liste des elements   VECSMALL de longueur au moins card H.
- * bit bitmap des elements VECSMALL de longueur au moins n.
- * retourne le nombre d'elements+1
- */
-long
-sousgroupeelem(long n, GEN v, GEN cy, GEN bit)
-{
-  int     j, m, o, p, flag;
-  for(j=1;j<n;j++) 
-    bit[j]=0;
-  m = 1;
-  bit[m] = 1;
-  cy[m++] = 1;
-  do
-  {
-    flag = 0;
-    for (o = 1; o < lg(v); o++)
-    {
-      for (p = 1; p < m; p++)	/* m varie! */
-      {
-	j = mulssmod(v[o],cy[p],n);
-	if (!bit[j])
-	{
-	  flag = 1;
-	  bit[j] = 1;
-	  cy[m++] = j;
-	}
-      }
-    }
-  }
-  while (flag);
-  return m;
-}
-/* n,v comme precedemment.
- * Calcule le conducteur et retourne le nouveau groupe de congruence dans V
- * V doit etre un t_VECSMALL de taille n+1 au moins.
- */
-long znconductor(long n, GEN v, GEN V)
-{
-  ulong ltop;
-  int i,j;
-  long m;
-  GEN F,W;
-  W = cgetg(n, t_VECSMALL);
-  ltop=avma;
-  m = sousgroupeelem(n,v,V,W);
-  setlg(V,m);
-  if (DEBUGLEVEL>=6)
-    fprintferr("SubCyclo:elements:%Z\n",V);
-  F = factor(stoi(n));  
-  for(i=lg((GEN)F[1])-1;i>0;i--)
-  {
-    long p,e,q;
-    p=itos(gcoeff(F,i,1));
-    e=itos(gcoeff(F,i,2));
-    if (DEBUGLEVEL>=4)
-      fprintferr("SubCyclo:testing %ld^%ld\n",p,e);
-    while (e>=1)
-    {
-      int z = 1;
-      q=n/p;
-      for(j=1;j<p;j++)
-      {
-	z += q;
-	if (!W[z] && z%p) break;
-      }
-      if (j<p)
-      {
-	if (DEBUGLEVEL>=4)
-	  fprintferr("SubCyclo:%ld not found\n",z);
-	break;
-      }
-      e--;
-      n=q;
-      if (DEBUGLEVEL>=4)
-	fprintferr("SubCyclo:new conductor:%ld\n",n);
-      m=sousgroupeelem(n,v,V,W);
-      setlg(V,m); 
-      if (DEBUGLEVEL>=6)
-	fprintferr("SubCyclo:elements:%Z\n",V);
-    }
-  } 
-  if (DEBUGLEVEL>=6)
-    fprintferr("SubCyclo:conductor:%ld\n",n);
-  avma=ltop;
-  return n;
-}
 
-static GEN gscycloconductor(GEN g, long n, long flag)
+GEN galoissubfields(GEN G, long flag, long v)
 {
-  if (flag==2)
-  {
-    GEN V=cgetg(3,t_VEC);
-    V[1]=lcopy(g);
-    V[2]=lstoi(n);
-    return V;
-  }
-  return g;
-}
-static GEN 
-lift_check_modulus(GEN H, GEN n)
-{
-  long t=typ(H), l=lg(H);
+  gpmem_t ltop=avma;
   long i;
-  GEN V;
-  switch(t)
-  {
-    case t_INTMOD:
-      if (cmpii(n,(GEN)H[1]))
-	err(talker,"wrong modulus in galoissubcyclo");
-      H = (GEN)H[2];
-    case t_INT:
-      if (!is_pm1(mppgcd(H,n)))
-	err(talker,"generators must be prime to conductor in galoissubcyclo");
-      return modii(H,n);
-    case t_VEC: case t_COL:
-      V=cgetg(l,t);
-      for(i=1;i<l;i++)
-	V[i] = (long)lift_check_modulus((GEN)H[i],n);
-      return V;
-    case t_VECSMALL:
-      return H;
-  }
-  err(talker,"wrong type in galoissubcyclo");
-  return NULL;/*not reached*/
+  GEN L = galoissubgroups(G);
+  long l2 = lg(L);
+  GEN p3 = cgetg(l2, t_VEC);
+  for (i = 1; i < l2; ++i)
+    p3[i] = (long) galoisfixedfield(G, gmael(L,i,1), flag, v);
+  return gerepileupto(ltop,p3);
 }
 
-GEN 
-galoissubcyclo(long n, GEN H, GEN Z, long v, long flag)
-{
-  ulong ltop=avma,av;
-  GEN l,borne,le,powz,z,V;
-  long i;
-  long e,val;
-  long u,o,j;
-  GEN O,g;
-  if (flag<0 || flag>2) err(flagerr,"galoisubcyclo");
-  if ( v==-1 ) v=0;
-  if ( n<1 ) err(arither2);
-  if ( n>=VERYBIGINT) 
-    err(impl,"galoissubcyclo for huge conductor");    
-  if ( typ(H)==t_MAT )
-  {
-    GEN zn2, zn3, gen, ord;
-    if (lg(H) == 1 || lg(H) != lg(H[1]))
-      err(talker,"not a HNF matrix in galoissubcyclo");
-    if (!Z)
-      Z=znstar(stoi(n));
-    else if (typ(Z)!=t_VEC || lg(Z)!=4) 
-      err(talker,"Optionnal parameter must be as output by znstar in galoissubcyclo");
-    zn2 = gtovecsmall((GEN)Z[2]);
-    zn3 = lift((GEN)Z[3]);
-    if ( lg(zn2) != lg(H) || lg(zn3) != lg(H))
-      err(talker,"Matrix of wrong dimensions in galoissubcyclo");
-    gen = cgetg(lg(zn3), t_VECSMALL);
-    ord = cgetg(lg(zn3), t_VECSMALL);
-    hnftogeneratorslist(n,zn2,zn3,H,gen,ord);
-    H=gen;
-  }
-  else
-  {
-    H=lift_check_modulus(H,stoi(n));
-    H=gtovecsmall(H);
-    for (i=1;i<lg(H);i++)
-      if (H[i]<0)
-	H[i]=mulssmod(-H[i],n-1,n);
-    /*Should check components are prime to n, but it is costly*/
-  }
-  V = cgetg(n, t_VECSMALL);
-  if (DEBUGLEVEL >= 1)
-    timer2();
-  n = znconductor(n,H,V);
-  if (flag==1)  {avma=ltop;return stoi(n);}
-  if (DEBUGLEVEL >= 1)
-    msgtimer("znconductor.");
-  H = V;
-  O = subgroupcoset(n,H);
-  if (DEBUGLEVEL >= 1)
-    msgtimer("subgroupcoset.");
-  if (DEBUGLEVEL >= 6)
-    fprintferr("Subcyclo: orbit=%Z\n",O);
-  if (lg(O)==1 || lg(O[1])==2)
-  {
-    avma=ltop;
-    return gscycloconductor(cyclo(n,v),n,flag);
-  }
-  u=lg(O)-1;o=lg(O[1])-1;
-  if (DEBUGLEVEL >= 4)
-    fprintferr("Subcyclo: %ld orbits with %ld elements each\n",u,o);
-  l=stoi(n+1);e=1;
-  while(!isprime(l)) 
-  { 
-    l=addis(l,n);
-    e++;
-  }
-  if (DEBUGLEVEL >= 4)
-    fprintferr("Subcyclo: prime l=%Z\n",l);
-  av=avma;
-  /*Borne utilise': 
-    Vecmax(Vec((x+o)^u)=max{binome(u,i)*o^i ;1<=i<=u} 
-  */
-  i=u-(1+u)/(1+o);
-  borne=mulii(binome(stoi(u),i),gpowgs(stoi(o),i));
-  if (DEBUGLEVEL >= 4)
-    fprintferr("Subcyclo: borne=%Z\n",borne);
-  val=mylogint(shifti(borne,1),l);
-  avma=av;
-  if (DEBUGLEVEL >= 4)
-    fprintferr("Subcyclo: val=%ld\n",val);
-  le=gpowgs(l,val);
-  z=lift(gpowgs(gener(l),e));
-  z=padicsqrtnlift(gun,stoi(n),z,l,val);
-  if (DEBUGLEVEL >= 1)
-    msgtimer("padicsqrtnlift.");
-  powz = cgetg(n,t_VEC); powz[1] = (long)z;
-  for (i=2; i<n; i++) powz[i] = lmodii(mulii(z,(GEN)powz[i-1]),le);
-  if (DEBUGLEVEL >= 1)
-    msgtimer("computing roots.");  
-  g=cgetg(u+1,t_VEC);
-  for(i=1;i<=u;i++)
-  {
-    GEN s;
-    s=gzero;
-    for(j=1;j<=o;j++)
-      s=addii(s,(GEN)powz[mael(O,i,j)]);
-    g[i]=(long)modii(negi(s),le);
-  }
-  if (DEBUGLEVEL >= 1)
-    msgtimer("computing new roots."); 
-  g=FpV_roots_to_pol(g,le,v);
-  if (DEBUGLEVEL >= 1)
-    msgtimer("computing products."); 
-  g=FpX_center(g,le);
-  return gerepileupto(ltop,gscycloconductor(g,n,flag));
-}