===================================================================
RCS file: /home/cvs/OpenXM_contrib/pari-2.2/src/modules/Attic/subfield.c,v
retrieving revision 1.1
retrieving revision 1.2
diff -u -p -r1.1 -r1.2
--- OpenXM_contrib/pari-2.2/src/modules/Attic/subfield.c	2001/10/02 11:17:12	1.1
+++ OpenXM_contrib/pari-2.2/src/modules/Attic/subfield.c	2002/09/11 07:27:06	1.2
@@ -1,4 +1,4 @@
-/* $Id: subfield.c,v 1.1 2001/10/02 11:17:12 noro Exp $
+/* $Id: subfield.c,v 1.2 2002/09/11 07:27:06 noro Exp $
 
 Copyright (C) 2000  The PARI group.
 
@@ -24,6 +24,10 @@ extern GEN matrixpow(long n, long m, GEN y, GEN P,GEN 
 extern GEN Fp_factor_irred(GEN P,GEN l, GEN Q);
 extern GEN FpX_rand(long d1, long v, GEN p);
 extern GEN hensel_lift_fact(GEN pol, GEN Q, GEN T, GEN p, GEN pe, long e);
+extern GEN ZX_disc_all(GEN x, ulong bound);
+extern GEN indexpartial(GEN P, GEN DP);
+extern GEN initgaloisborne(GEN T, GEN dn, long prec, GEN *ptL, GEN *ptprep, GEN *ptdis);
+extern long ZX_get_prec(GEN x);
 
 static GEN print_block_system(long N,GEN Y,long d,GEN vbs,long maxl);
 
@@ -83,7 +87,7 @@ calc_block(long N,GEN Z,long d,GEN Y,GEN vbs,ulong max
       }
       if (dk == nn)
       {
-	long av=avma;
+	gpmem_t av=avma;
         int Z_equal_Zp = 1;
 
         for (j=1; j<lnon; j++) non[j]=0;
@@ -108,7 +112,7 @@ calc_block(long N,GEN Z,long d,GEN Y,GEN vbs,ulong max
           setlg(Zpp, u);
 	  vbs = calc_block(N,Zpp,d,Yp,vbs,maxl);
 	}
-        if (vbs && lg(vbs) > maxl) return vbs;
+        if (vbs && (ulong)lg(vbs) > maxl) return vbs;
         avma=av;
       }
     }
@@ -119,7 +123,8 @@ calc_block(long N,GEN Z,long d,GEN Y,GEN vbs,ulong max
 static GEN
 potential_block_systems(long N, long d, GEN n, ulong maxl)
 {
-  long av=avma,r,i,j,k;
+  long r, i, j, k;
+  gpmem_t av=avma;
   GEN p1,vbs,Z;
 
   r=lg(n); Z=cgetg(r,t_VEC);
@@ -134,9 +139,10 @@ potential_block_systems(long N, long d, GEN n, ulong m
 
 /* product of permutations. Put the result in perm1. */
 static void
-perm_mul(GEN perm1,GEN perm2)
+perm_mul_i(GEN perm1,GEN perm2)
 {
-  long av = avma,i, N = lg(perm1);
+  long i, N = lg(perm1);
+  gpmem_t av = avma;
   GEN perm=new_chunk(N);
   for (i=1; i<N; i++) perm[i]=perm1[perm2[i]];
   for (i=1; i<N; i++) perm1[i]=perm[i];
@@ -153,14 +159,14 @@ cycle_power_to_perm(GEN perm,GEN cy,long l)
   for (i=1; i<N; i++) perm[i] = i;
   if (lp)
   {
-    long av = avma;
+    gpmem_t av = avma;
     GEN p1 = new_chunk(N);
     b = cy[1];
     for (i=1; i<lcy; i++) b = (perm[b] = cy[i+1]);
     perm[b] = cy[1];
     for (i=1; i<N; i++) p1[i] = perm[i];
 
-    for (j=2; j<=lp; j++) perm_mul(perm,p1);
+    for (j=2; j<=lp; j++) perm_mul_i(perm,p1);
     avma = av;
   }
   return perm;
@@ -221,8 +227,7 @@ append_vbs(GEN vbs, GEN D)
   l = lg(vbs); maxl = vbs[-1];
   if (l == maxl)
   {
-    vbs = (GEN)gprealloc((void*)(vbs-1), (2 + (maxl<<1))*sizeof(GEN),
-                                         (2 + maxl)*sizeof(GEN));
+    vbs = (GEN)gprealloc((void*)(vbs-1), (2 + (maxl<<1))*sizeof(GEN));
     *vbs = maxl<<1; vbs++; setlg(vbs, 1);
   }
   if (DEBUGLEVEL>5) fprintferr("appending D = %Z\n",D);
@@ -256,7 +261,7 @@ print_block_system(long N,GEN Y,long d,GEN vbs,long ma
   if (DEBUGLEVEL>5) fprintferr("Y = %Z\n",Y);
   empty = cgetg(1,t_VEC);
   n = new_chunk(N+1);
-  D = cgetg(N+1, t_VEC); setlg(D,1);
+  D = cget1(N+1, t_VEC);
   t=new_chunk(r+1); k=new_chunk(r+1); Z=cgetg(r+1,t_VEC);
   for (ns=0,i=1; i<r; i++)
   {
@@ -299,7 +304,7 @@ print_block_system(long N,GEN Y,long d,GEN vbs,long ma
   perm = cgetg(N+1,t_VEC); i=ns;
   do
   {
-    long av = avma;
+    gpmem_t av = avma;
     if (DEBUGLEVEL>5)
     {
       for (l=1; l<=ns; l++)
@@ -313,7 +318,7 @@ print_block_system(long N,GEN Y,long d,GEN vbs,long ma
     for (u=1; u<=N; u++) perm[u]=u;
     for (u=1; u<=ns; u++)
       for (v=1; v<=t[u]; v++)
-	perm_mul(perm, cycle_power_to_perm(cyperm,gmael(Z,u,v),e[u][v]));
+	perm_mul_i(perm, cycle_power_to_perm(cyperm,gmael(Z,u,v),e[u][v]));
     vbs = append_vbs(vbs, im_block_by_perm(D,perm));
     if (lg(vbs) > maxl) return vbs;
     avma = av;
@@ -376,7 +381,7 @@ cand_for_subfields(GEN A,GEN DATA,GEN *ptlistdelta)
 {
   long N,m,i,j,d,lf;
   GEN M,T,pe,p,pol,fhk,g;
-  GEN _d_1_term,delta,listdelta,whichdelta,firstroot;
+  GEN d_1_term,delta,listdelta,whichdelta,firstroot;
 
   pol=(GEN)DATA[1];
   p = (GEN)DATA[2];
@@ -390,12 +395,12 @@ cand_for_subfields(GEN A,GEN DATA,GEN *ptlistdelta)
   lf = lg(firstroot);
   listdelta = cgetg(lf, t_VEC);
   whichdelta = cgetg(N+1, t_VECSMALL);
-  _d_1_term = gzero;
+  d_1_term = gzero;
   for (i=1; i<=m; i++)
   {
     GEN Ai = (GEN)A[i], p1 = gun;
     for (j=1; j<=d; j++)
-      p1 = FpXQ_mul(p1, (GEN)fhk[Ai[j]], T,pe);
+      p1 = Fq_mul(p1, (GEN)fhk[Ai[j]], T,pe);
     delta[i] = (long)p1;
     if (DEBUGLEVEL>2) fprintferr("delta[%ld] = %Z\n",i,p1);
     /* g = prod (X - delta[i])
@@ -403,10 +408,10 @@ cand_for_subfields(GEN A,GEN DATA,GEN *ptlistdelta)
     /* fk[k] belongs to block number whichdelta[k] */
     for (j=1; j<=d; j++) whichdelta[Ai[j]] = i;
     if (typ(p1) == t_POL) p1 = constant_term(p1);
-    _d_1_term = addii(_d_1_term, p1);
+    d_1_term = addii(d_1_term, p1);
   }
-  _d_1_term = centermod(_d_1_term, pe); /* Tr(g) */
-  if (absi_cmp(_d_1_term, (GEN)M[3]) > 0) return gdeux; /* d-1 test */
+  d_1_term = centermod(d_1_term, pe); /* Tr(g) */
+  if (absi_cmp(d_1_term, (GEN)M[3]) > 0) return gdeux; /* d-1 test */
   g = FqV_roots_to_pol(delta, T, pe, 0);
   g = centermod(polsimplify(g), pe); /* assume g in Z[X] */
   if (DEBUGLEVEL>2) fprintferr("pol. found = %Z\n",g);
@@ -497,6 +502,36 @@ chinese_retrieve_pol(GEN DATA, GEN listdelta)
   return FpX_res(FpX_red(S, p), FpX_red((GEN)DATA[1],p), p);
 }
 
+/* return P(X + c) using destructive Horner */
+static GEN
+TR_pol(GEN P, GEN c)
+{
+  GEN Q = dummycopy(P), *R;
+  long i,k,n;
+  
+  if (!signe(P)) return Q;
+  R = (GEN*)(Q+2); n = degpol(P);
+  if (gcmp1(c))
+  {
+    for (i=1; i<=n; i++)
+      for (k=n-i; k<n; k++)
+        R[k] = gadd(R[k], R[k+1]);
+  }
+  else if (gcmp_1(c))
+  {
+    for (i=1; i<=n; i++)
+      for (k=n-i; k<n; k++)
+        R[k] = gsub(R[k], R[k+1]);
+  }
+  else
+  {
+    for (i=1; i<=n; i++)
+      for (k=n-i; k<n; k++)
+        R[k] = gadd(R[k], gmul(c, R[k+1]));
+  }
+  return Q;
+}
+
 /* g in Z[X] potentially defines a subfield of Q[X]/f. It is a subfield iff A
  * (cf cand_for_subfields) was a block system; then there
  * exists h in Q[X] such that f | g o h. listdelta determines h s.t f | g o h
@@ -506,7 +541,7 @@ embedding_of_potential_subfields(GEN g,GEN DATA,GEN li
 {
   GEN TR,w0_Q,w0,w1_Q,w1,wpow,h0,gp,T,q2,q,p,ind,maxp,a;
   long rt;
-  ulong av;
+  gpmem_t av;
 
   T = (GEN)DATA[1]; rt = brent_kung_optpow(degpol(T), 2);
   p = (GEN)DATA[2];
@@ -535,7 +570,7 @@ embedding_of_potential_subfields(GEN g,GEN DATA,GEN li
     w1_Q = centermod(gmul(w1, resii(ind,q)), q);
     if (gegal(w1_Q, w0_Q) || cmpii(q,maxp) > 0)
     {
-      GEN G = gcmp1(ind)? g: ZX_rescale_pol(g,ind);
+      GEN G = gcmp1(ind)? g: rescale_pol(g,ind);
       if (gcmp0(poleval(G, gmodulcp(w1_Q,T)))) break;
     }
     if (cmpii(q, maxp) > 0)
@@ -560,39 +595,42 @@ embedding_of_potential_subfields(GEN g,GEN DATA,GEN li
     w0 = w1; w0_Q = w1_Q; p = q; q = q2;
   }
   TR = (GEN)DATA[14];
-  if (!gcmp0(TR)) w1_Q = poleval(w1_Q, gadd(polx[0], TR));
+  if (!gcmp0(TR)) w1_Q = TR_pol(w1_Q, TR);
   return gdiv(w1_Q,ind);
 }
 
 static GEN
 choose_prime(GEN pol,GEN dpol,long d,GEN *ptff,GEN *ptlistpotbl, long *ptlcm)
 {
-  ulong av;
-  long j,k,oldllist,llist,r,lcm,oldlcm,N,pp;
-  GEN p,listpotbl,oldlistpotbl,ff,oldff,n,oldn;
+  gpmem_t av;
   byteptr di=diffptr;
-
-  if (DEBUGLEVEL) timer2();
-  di++; p = stoi(2); N = degpol(pol);
-  while (p[2]<=N) p[2] += *di++;
+  long j,k,oldllist,llist,r,lcm,oldlcm,pp,minp, N = degpol(pol), m = N/d;
+  GEN p,listpotbl,oldlistpotbl,ff,oldff,n,oldn;
+ 
+  minp = N*(m-1)/2;
+  if (DEBUGLEVEL) (void)timer2();
+  di++; p = stoi(2);
+  while (p[2]<=minp)
+      NEXT_PRIME_VIADIFF(p[2], di);
   oldllist = 100000;
   oldlcm = 0;
   oldlistpotbl = oldff = oldn = NULL; pp = 0; /* gcc -Wall */
   av = avma;
-  for(k=1; k<11 || !oldn; k++,p[2]+= *di++,avma = av)
+  for(k=1; k<11 || !oldn; k++,avma = av)
   {
-    while (!smodis(dpol,p[2])) p[2] += *di++;
+    while (!smodis(dpol,p[2]))
+      NEXT_PRIME_VIADIFF(p[2], di);
     if (k > 50) err(talker,"sorry, too many block systems in nfsubfields");
-    ff=(GEN)factmod(pol,p)[1]; r=lg(ff)-1;
-    if (r == 1 || r == N) continue;
+    ff=(GEN)factmod0(pol,p)[1]; r=lg(ff)-1;
+    if (r == 1 || r == N) goto repeat;
 
     n = cgetg(r+1, t_VECSMALL);
     lcm = n[1] = degpol(ff[1]);
     for (j=2; j<=r; j++) { n[j] = degpol(ff[j]); lcm = clcm(lcm,n[j]); }
-    if (lcm < oldlcm) continue; /* false when oldlcm = 0 */
-    if (r >= BIL) { err(warner,"subfields: overflow in calc_block"); continue; }
+    if (lcm < oldlcm) goto repeat; /* false when oldlcm = 0 */
+    if (r >= BIL) { err(warner,"subfields: overflow in calc_block"); goto repeat; }
     if (DEBUGLEVEL) fprintferr("p = %ld,\tlcm = %ld,\torbits: %Z\n",p[2],lcm,n);
-    if (oldn && gegal(n,oldn)) continue;
+    if (oldn && gegal(n,oldn)) goto repeat;
 
     listpotbl = potential_block_systems(N,d,n, oldllist);
     if (!listpotbl) { oldlistpotbl = NULL; pp = p[2]; break; }
@@ -601,21 +639,16 @@ choose_prime(GEN pol,GEN dpol,long d,GEN *ptff,GEN *pt
     {
       if (DEBUGLEVEL) msgtimer("#pbs >= %ld [aborted]",oldllist);
       for (j=1; j<llist; j++) free((void*)listpotbl[j]);
-      free((void*)(listpotbl-1)); continue;
+      free((void*)(listpotbl-1)); goto repeat;
     }
     oldllist = llist; oldlistpotbl = listpotbl;
     pp = p[2]; oldff = ff; oldlcm = lcm; oldn = n;
     if (DEBUGLEVEL) msgtimer("#pbs = %ld",llist);
     av = avma;
+   repeat:
+    NEXT_PRIME_VIADIFF(p[2], di);
   }
-  if (DEBUGLEVEL)
-  {
-    fprintferr("Chosen prime: p = %ld\n",pp);
-    if (DEBUGLEVEL>2)
-      fprintferr("Potential block systems of size %ld: %Z\n", d,oldlistpotbl);
-    flusherr();
-  }
-  if (oldff) oldff = lift_intern(oldff);
+  if (DEBUGLEVEL) fprintferr("Chosen prime: p = %ld\n",pp);
   *ptlistpotbl=oldlistpotbl; *ptff=oldff; *ptlcm=oldlcm; return stoi(pp);
 }
 
@@ -625,7 +658,9 @@ bound_for_coeff(long m,GEN rr, GEN *maxroot)
   long i,r1, lrr=lg(rr);
   GEN p1,b1,b2,B,M, C = matpascal(m-1);
 
-  for (r1=0; typ(rr[r1+1]) == t_REAL; r1++) /* empty */;
+  for (r1=1; r1 < lrr; r1++)
+    if (typ(rr[r1]) != t_REAL) break;
+  r1--;
 
   rr = gabs(rr,0); *maxroot = vecmax(rr);
   for (i=1; i<lrr; i++)
@@ -689,7 +724,7 @@ interpol(GEN H, GEN T, GEN p)
   {
     d = constant_term(H[i]); /* -D[i] */
     p1 = FpXQX_mul(p1,gadd(X,d), T,p);
-    p2 = FpXQ_mul(p2, gadd(a,d), T,p);
+    p2 = Fq_mul(p2, gadd(a,d), T,p);
   }
   p2 = FpXQ_inv(p2, T,p);
   return FpXQX_FpXQ_mul(p1,p2, T,p);
@@ -732,8 +767,8 @@ struct poldata
 static GEN
 compute_data(GEN DATA, struct poldata PD, long d, GEN ff, GEN T,GEN p)
 {
-  long i,j,l,e,N;
-  GEN den,roo,pe,p1,p2,fk,fhk,MM,maxroot,pol,interp,bezoutC;
+  long i,j,l,e,N, lff = lg(ff);
+  GEN ffL,den,roo,pe,p1,p2,fk,fhk,MM,maxroot,pol,interp,bezoutC;    
 
   if (DEBUGLEVEL>1) { fprintferr("Entering compute_data()\n\n"); flusherr(); }
   pol = PD.pol; N = degpol(pol);
@@ -743,8 +778,9 @@ compute_data(GEN DATA, struct poldata PD, long d, GEN 
   {
     GEN Xm1 = gsub(polx[varn(pol)], gun);
     GEN TR = addis((GEN)DATA[14],1);
+    GEN mTR = negi(TR), mun = negi(gun);
     DATA[14] = (long)TR;
-    pol = poleval(pol, gsub(polx[varn(pol)], TR));
+    pol = TR_pol(pol, mTR);
     p1 = dummycopy(roo); l = lg(p1);
     for (i=1; i<l; i++) p1[i] = ladd(TR, (GEN)p1[i]);
     roo = p1;
@@ -758,29 +794,31 @@ compute_data(GEN DATA, struct poldata PD, long d, GEN 
     {
       if (degpol(interp[i]) > 0) /* do not turn polun[0] into gun */
       {
-        p1 = poleval((GEN)interp[i], Xm1);
+        p1 = TR_pol((GEN)interp[i], mun);
         interp[i] = (long)FpXQX_red(p1, NULL,p);
       }
       if (degpol(bezoutC[i]) > 0)
       {
-        p1 = poleval((GEN)bezoutC[i], Xm1);
+        p1 = TR_pol((GEN)bezoutC[i], mun);
         bezoutC[i] = (long)FpXQX_red(p1, NULL,p);
       }
     }
+    p1 = cgetg(lff, t_VEC);
+    for (i=1; i<lff; i++)
+      p1[i] = (long)FpX_red(TR_pol((GEN)ff[i], mTR), p);
+    ff = p1;
   }
   else
   {
-    GEN firstroot;    
-    long r = lg(ff);
+    GEN firstroot = cgetg(lff, t_VECSMALL);
     DATA = cgetg(15,t_VEC);
     DATA[2] = (long)p;
     DATA[4] = (long)T;
-    interp = cgetg(r, t_VEC);
-    firstroot = cgetg(r, t_VECSMALL);
+    interp = cgetg(lff, t_VEC);
     fk = cgetg(N+1,t_VEC);
-    for (l=1,j=1; j<r; j++)
+    for (l=1,j=1; j<lff; j++)
     { /* compute roots and fix ordering (Frobenius cycles) */
-      p1 = Fp_factor_irred((GEN)ff[j], p, (GEN)DATA[4]);
+      p1 = Fp_factor_irred((GEN)ff[j], p,T);
       interp[j] = (long)interpol(p1,T,p);
       firstroot[j] = l;
       for (i=1; i<lg(p1); i++,l++) fk[l] = p1[i];
@@ -799,7 +837,23 @@ compute_data(GEN DATA, struct poldata PD, long d, GEN 
   e = logint(shifti(vecmax(MM),20), p, &pe); /* overlift 2^20 [for d-1 test] */
   DATA[3] = (long)pe;
   DATA[6] = (long)roots_from_deg1(fk);
+
+#if 0
   fhk = hensel_lift_fact(pol,fk,(GEN)DATA[4],p,pe,e);
+#else
+  /* first lift in Zp to precision p^e */
+  ffL = hensel_lift_fact(pol,ff, NULL,p,pe,e);
+  fhk = NULL;
+  for (l=i=1; i<lff; i++)
+  { /* lift factorization of ff[i] in Qp[X] / T */
+    GEN L = (GEN)ffL[i];
+    long di = degpol(L);
+    p2 = cgetg(di+1, t_VEC);
+    for (j=1; j<=di; j++) p2[j] = fk[l++];
+    p1 = hensel_lift_fact(L,p2,T,p,pe,e);
+    fhk = fhk? concatsp(fhk, p1): p1;
+  }
+#endif
   DATA[5] = (long)roots_from_deg1(fhk);
 
   p1 = gmul(stoi(N), gsqrt(gpuigs(stoi(N-1),N-1),DEFAULTPREC));
@@ -831,8 +885,8 @@ _subfield(GEN g, GEN h)
 static GEN
 subfields_of_given_degree(struct poldata PD,long d)
 {
-  ulong av,av2;
-  long llist,i,nn;
+  gpmem_t av, av2;
+  long llist,i,nn,v;
   GEN listpotbl,ff,A,CSF,ESF,LSB,p,T,DATA,listdelta;
   GEN pol = PD.pol, dpol = PD.dis;
 
@@ -840,7 +894,8 @@ subfields_of_given_degree(struct poldata PD,long d)
   av = avma;
   p = choose_prime(pol,dpol,degpol(pol)/d,&ff,&listpotbl,&nn);
   if (!listpotbl) { avma=av; return cgetg(1,t_VEC); }
-  T = lift_intern(ffinit(p,nn, fetch_var()));
+  v = fetch_var(); name_var(v,"y");
+  T = lift_intern(ffinit(p,nn, v));
   DATA = NULL; LSB = cgetg(1,t_VEC); 
   i = 1; llist = lg(listpotbl);
 CHANGE: 
@@ -871,7 +926,7 @@ CHANGE: 
       LSB = gerepileupto(av2, concat(LSB, _subfield(CSF,ESF)));
     }
   }
-  delete_var();
+  (void)delete_var();
   for (i=1; i<llist; i++) free((void*)listpotbl[i]);
   free((void*)(listpotbl-1));
   if (DEBUGLEVEL) fprintferr("\nSubfields of degree %ld: %Z\n",d, LSB);
@@ -887,69 +942,33 @@ fix_var(GEN x, long v)
   return x;
 }
 
-extern GEN get_nfpol(GEN x, GEN *nf);
-extern GEN vandermondeinverseprep(GEN T, GEN L);
-extern GEN ZX_disc_all(GEN x, ulong bound);
-extern GEN indexpartial(GEN P, GEN DP);
-
 void
 subfields_poldata(GEN T, struct poldata *PD)
 {
-  int     i, n;
-  GEN     L, z, prep, den;
-  long    prec;
+  GEN  nf,L,dis;
 
-  GEN nf, dis;
   T = get_nfpol(T, &nf);
   PD->pol = dummycopy(T); /* may change variable number later */
   if (nf)
   {
-    PD->den = (GEN)nf[4];
+    PD->den = Q_denom((GEN)nf[7]);
     PD->roo = (GEN)nf[6];
     PD->dis = mulii(absi((GEN)nf[3]), sqri((GEN)nf[4]));
-    return;
   }
-
-  /* copy-paste from galconj.c:galoisborne. Merge as soon as possible */
-  /* start of copy-paste */
-  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);
-  L = roots(T, prec);
-  if (DEBUGLEVEL>=4) genmsgtimer(3,"roots");
-  for (i = 1; i <= n; i++)
+  else
   {
-    z = (GEN) L[i];
-    if (signe(z[2]))
-      break;
-    L[i] = z[1];
+    PD->den = initgaloisborne(T,NULL,ZX_get_prec(T), &L,NULL,&dis);
+    PD->roo = L;
+    PD->dis = absi(dis);
   }
-  if (DEBUGLEVEL>=4) gentimer(3);
-  prep=vandermondeinverseprep(T, L);
-  /* end of copy-paste */
-  {
-    GEN res = divide_conquer_prod(gabs(prep,prec), mpmul);
-    disable_dbg(0);
-    dis = ZX_disc_all(T, 1+logint(res,gdeux,NULL));
-    disable_dbg(-1);
-    den = indexpartial(T,dis);
-  }
-
-  PD->den = den;
-  PD->roo = L;
-  PD->dis = absi(dis);
 }
 
 GEN
 subfields(GEN nf,GEN d)
 {
-  ulong av = avma;
-  long di,N,v0;
-  GEN LSB,pol;
+  gpmem_t av = avma;
+  long di, N, v0;
+  GEN LSB, pol, G;
   struct poldata PD;
 
   pol = get_nfpol(nf, &nf); /* in order to treat trivial cases */
@@ -958,7 +977,25 @@ subfields(GEN nf,GEN d)
   if (di==1) return gerepilecopy(av, _subfield(polx[v0], pol));
   if (di < 1 || di > N || N % di) return cgetg(1,t_VEC);
 
-  subfields_poldata(nf, &PD);
+  /* much easier if nf is Galois (WSS) */
+  G = galoisconj4(nf? nf: pol, NULL, 1, 0);
+  if (typ(G) != t_INT)
+  { /* Bingo */
+    GEN L = galoissubgroups(G), F;
+    long k,i, l = lg(L), o = N/di;
+    F = cgetg(l, t_VEC);
+    k = 1;
+    for (i=1; i<l; i++)
+    {
+      GEN H = (GEN)L[i];
+      if (group_order(H) == o)
+        F[k++] = (long) lift_intern(galoisfixedfield(G, (GEN)H[1], 0, v0));
+    }
+    setlg(F, k);
+    return gerepilecopy(av, F);
+  }
+
+  subfields_poldata(nf? nf:pol, &PD);
   pol = PD.pol;
   setvarn(pol, 0);
   LSB = subfields_of_given_degree(PD, di);
@@ -968,11 +1005,19 @@ subfields(GEN nf,GEN d)
 static GEN
 subfieldsall(GEN nf)
 {
-  ulong av = avma;
-  long N,ld,i,v0;
-  GEN pol,dg,LSB,NLSB;
+  gpmem_t av = avma;
+  long N, ld, i, v0;
+  GEN G, pol, dg, LSB, NLSB;
   struct poldata PD;
 
+  /* much easier if nf is Galois (WSS) */
+  G = galoisconj4(nf, NULL, 1, 0);
+  if (typ(G) != t_INT)
+  {
+    pol = get_nfpol(nf, &nf);
+    return gerepilecopy(av, lift_intern( galoissubfields(G, 0, varn(pol)) ));
+  }
+
   subfields_poldata(nf, &PD);
   pol = PD.pol;
 
@@ -986,7 +1031,7 @@ subfieldsall(GEN nf)
     setvarn(pol, 0);
     for (i=2; i<ld; i++)
     {
-      ulong av1 = avma;
+      gpmem_t av1 = avma;
       NLSB = subfields_of_given_degree(PD, N / itos((GEN)dg[i]));
       if (lg(NLSB) > 1) LSB = concatsp(LSB,NLSB); else avma = av1;
     }
@@ -1002,20 +1047,3 @@ subfields0(GEN nf,GEN d)
   return d? subfields(nf,d): subfieldsall(nf);
 }
 
-/* irreducible (unitary) polynomial of degree n over Fp[v] */
-GEN
-ffinit(GEN p,long n,long v)
-{
-  long av = avma;
-  GEN pol;
-
-  if (n<=0) err(talker,"non positive degree in ffinit");
-  if (typ(p) != t_INT) err(typeer,"ffinit");
-  if (v<0) v = 0;
-  for(;; avma = av)
-  {
-    pol = gadd(gpowgs(polx[v],n), FpX_rand(n-1,v, p));
-    if (FpX_is_irred(pol, p)) break;
-  }
-  return gerepileupto(av, FpX(pol,p));
-}