===================================================================
RCS file: /home/cvs/OpenXM_contrib/pari-2.2/src/modules/Attic/kummer.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/kummer.c	2001/10/02 11:17:11	1.1
+++ OpenXM_contrib/pari-2.2/src/modules/Attic/kummer.c	2002/09/11 07:27:05	1.2
@@ -1,4 +1,4 @@
-/* $Id: kummer.c,v 1.1 2001/10/02 11:17:11 noro Exp $
+/* $Id: kummer.c,v 1.2 2002/09/11 07:27:05 noro Exp $
 
 Copyright (C) 2000  The PARI group.
 
@@ -19,1212 +19,1067 @@ Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 /*                                                                 */
 /*******************************************************************/
 #include "pari.h"
+#include "parinf.h"
+extern GEN gmul_mati_smallvec(GEN x, GEN y);
 extern GEN check_and_build_cycgen(GEN bnf);
 extern GEN get_arch_real(GEN nf,GEN x,GEN *emb,long prec);
-extern GEN vecmul(GEN x, GEN y);
-extern GEN vecinv(GEN x);
 extern GEN T2_from_embed_norm(GEN x, long r1);
-extern GEN pol_to_vec(GEN x, long N);
-extern GEN famat_to_nf(GEN nf, GEN f);
+extern GEN vconcat(GEN A, GEN B);
+extern long int_elt_val(GEN nf, GEN x, GEN p, GEN b, GEN *newx);
 
-static long rc,ell,degK,degKz,m,d,vnf,dv;
-static GEN matexpoteta1,nf,raycyc,polnf;
-static GEN bnfz,nfz,U,uu,gell,cyc,gencyc,vecalpha,R,g;
-static GEN listmod,listprSp,listbid,listunif,listellrank;
-static GEN listbidsup,listellranksup,vecw;
+extern GEN famat_inv(GEN f);
+extern GEN famat_pow(GEN f, GEN n);
+extern GEN famat_mul(GEN f, GEN g);
+extern GEN famat_reduce(GEN fa);
+extern GEN famat_ideallog(GEN nf, GEN g, GEN e, GEN bid);
+extern GEN to_famat(GEN g, GEN e);
 
-/* row vector B x matrix T : c_j=prod_i (b_i^T_ij) */
-static GEN
-groupproduct(GEN B, GEN T)
-{
-  long lB,lT,i,j;
-  GEN c,p1;
+typedef struct {
+  GEN x;  /* tau ( Mod(x, nf.pol) ) */
+  GEN zk; /* action of tau on nf.zk (as t_MAT) */
+} tau_s;
 
-  lB=lg(B)-1;
-  lT=lg(T)-1;
-  c=cgetg(lT+1,t_VEC);
-  for (j=1; j<=lT; j++)
-  {
-    p1=gun;
-    for (i=1; i<=lB; i++) p1=gmul(p1,gpui((GEN)B[i],gcoeff(T,i,j),0));
-    c[j]=(long)p1;
-  }
-  return c;
-}
+typedef struct {
+  GEN polnf, invexpoteta1;
+  tau_s *tau;
+  long m;
+} toK_s;
 
-static GEN
-grouppows(GEN B, long ex)
+long
+prank(GEN cyc, long ell)
 {
-  long lB = lg(B),j;
-  GEN c;
-
-  c = cgetg(lB,t_VEC);
-  for (j=1; j<lB; j++) c[j] = lpowgs((GEN)B[j], ex);
-  return c;
+  long i;
+  for (i=1; i<lg(cyc); i++)
+    if (smodis((GEN)cyc[i],ell)) break;
+  return i-1;
 }
 
+/* increment y, which runs through [0,ell-1]^(k-1). Return 0 when done. */
 static int
-ok_x(GEN X, GEN arch, GEN vecmunit2, GEN msign)
+increment(GEN y, long k, long ell)
 {
-  long i, l = lg(vecmunit2);
-  GEN p1;
-  for (i=1; i<l; i++)
+  long i = k, j;
+  do
   {
-    p1 = FpV_red(gmul((GEN)vecmunit2[i], X), gell);
-    if (gcmp0(p1)) return 0;
-  }
-  p1 = lift(gmul(msign,X)); settyp(p1,t_VEC);
-  return gegal(p1, arch);
+    if (--i == 0) return 0;
+    y[i]++;
+  } while (y[i] >= ell);
+  for (j = i+1; j < k; j++) y[j] = 0;
+  return 1;
 }
 
-static GEN
-get_listx(GEN arch,GEN msign,GEN munit,GEN vecmunit2,long ell,long lSp,long nbvunit)
+/* as above, y increasing (y[i] <= y[i+1]) */
+static int
+increment_inc(GEN y, long k, long ell)
 {
-  GEN kermat,p2,X,y, listx = cgetg(1,t_VEC);
-  long i,j,cmpt,lker;
-
-  kermat = FpM_ker(munit,gell); lker=lg(kermat)-1;
-  if (!lker) return listx;
-  y = cgetg(lker,t_VECSMALL);
-  for (i=1; i<lker; i++) y[i] = 0;
-  for(;;)
+  long i = k, j;
+  do
   {
-    p2 = cgetg(2,t_VEC);
-    X = (GEN)kermat[lker];
-    for (j=1; j<lker; j++) X = gadd(X, gmulsg(y[j],(GEN)kermat[j]));
-    X = FpV_red(X, gell);
-    cmpt = 0; for (j=1; j<=lSp; j++) if (gcmp0((GEN)X[j+nbvunit])) cmpt++;
-    if (!cmpt)
-    {
-      if (ok_x(X, arch, vecmunit2, msign))
-        { p2[1] = (long)X; listx = concatsp(listx,p2); }
-    }
-    if (!lSp)
-    {
-      j = 1; while (j<lker && !y[j]) j++;
-      if (j<lker && y[j] == 1)
-      {
-        X = gsub(X,(GEN)kermat[lker]);
-        if (ok_x(X, arch, vecmunit2, msign))
-          { p2[1] = (long)X; listx = concatsp(listx,p2); }
-      }
-    }
-    i = lker;
-    do
-    {
-      i--; if (!i) return listx;
-      if (i < lker-1) y[i+1] = 0;
-      y[i]++;
-    }
-    while (y[i] >= ell);
-  }
+    if (--i == 0) return 0;
+    y[i]++;
+  } while (y[i] >= ell);
+  for (j = i+1; j < k; j++) y[j] = y[i];
+  return 1;
 }
 
-static GEN
-reducealpha(GEN nf,GEN x,GEN gell)
-/* etant donne un nombre algebrique x du corps nf -- non necessairement
-   entier -- calcule un entier algebrique de la forme x*g^gell */
+static int
+ok_congruence(GEN X, GEN ell, long lW, GEN vecMsup)
 {
-  long tx=typ(x),i;
-  GEN den,fa,fac,ep,p1,y;
-
-  nf = checknf(nf);
-  if (tx==t_POL || tx==t_POLMOD) y = algtobasis(nf,x);
-  else { y = x; x = basistoalg(nf,y); }
-  den = denom(y);
-  if (gcmp1(den)) return x;
-  fa = decomp(den); fac = (GEN)fa[1];ep = (GEN)fa[2];
-  p1 = gun;
-  for (i=1; i<lg(fac); i++)
-    p1 = mulii(p1, powgi((GEN)fac[i], gceil(gdiv((GEN)ep[i],gell))));
-  return gmul(powgi(p1, gell), x);
+  long i, l;
+  if (gcmp0(X)) return 0;
+  l = lg(X);
+  for (i=lW; i<l; i++)
+    if (gcmp0((GEN)X[i])) return 0;
+  l = lg(vecMsup);
+  for (i=1; i<l; i++)
+    if (gcmp0(FpV_red(gmul((GEN)vecMsup[i],X), ell))) return 0;
+  return 1;
 }
 
-long
-ellrank(GEN cyc, long ell)
+static int
+ok_sign(GEN X, GEN msign, GEN arch)
 {
-  long i;
-  for (i=1; i<lg(cyc); i++)
-    if (smodis((GEN)cyc[i],ell)) break;
-  return i-1;
+  GEN p1 = lift(gmul(msign,X)); settyp(p1,t_VEC);
+  return gegal(p1, arch);
 }
 
+/* REDUCTION MOD ell-TH POWERS */
+
 static GEN
-no_sol(long all, int i)
+fix_be(GEN bnfz, GEN be, GEN u)
 {
-  if (!all) err(talker,"bug%d in kummer",i);
-  return cgetg(1,t_VEC);
+  GEN nf = checknf(bnfz), fu = gmael(bnfz,8,5);
+  return element_mul(nf, be, factorbackelt(fu, u, nf));
 }
 
-static void
-_append(GEN x, GEN t)
+static GEN
+logarch2arch(GEN x, long r1, long prec)
 {
-  long l = lg(x); x[l] = (long)t; setlg(x, l+1);
+  long i, lx = lg(x), tx = typ(x);
+  GEN y = cgetg(lx, tx);
+
+  if (tx == t_MAT)
+  {
+    for (i=1; i<lx; i++) y[i] = (long)logarch2arch((GEN)x[i], r1, prec);
+    return y;
+  }
+  for (i=1; i<=r1;i++) y[i] = lexp((GEN)x[i],prec);
+  for (   ; i<lx; i++) y[i] = lexp(gmul2n((GEN)x[i],-1),prec);
+  return y;
 }
 
-/* si all!=0, donne toutes les equations correspondant a un sousgroupe
-   de determinant all (i.e. de degre all) */
+/* multiply be by ell-th powers of units as to find small L2-norm for new be */
 static GEN
-rnfkummersimple(GEN bnr, GEN subgroup, long all)
+reducebetanaive(GEN bnfz, GEN be, GEN b, GEN ell)
 {
-  long r1,degnf,ell,j,i,l;
-  long nbgenclK,lSml2,lSl,lSp,rc,nbvunit;
-  long lastbid,nbcol,nblig,k,llistx,e,vp,vd;
+  long i,k,n,ru,r1, prec = nfgetprec(bnfz);
+  GEN z,p1,p2,nmax,c, nf = checknf(bnfz);
 
-  GEN bnf,nf,classgroup,cyclic,bid,ideal,arch,cycgen,gell,p1,p2,p3;
-  GEN cyclicK,genK,listgamma,listalpha,fa,prm,expom,Sm,Sml1,Sml2,Sl,ESml2;
-  GEN p,factell,Sp,vecbeta,matexpo,vunit,id,pr,vsml2,vecalpha0;
-  GEN cycpro,munit,vecmunit2,msign,archartif,listx,listal,listg,listgamma0;
-  GEN vecbeta0;
-
-  checkbnrgen(bnr);
-  bnf = (GEN)bnr[1];
-  nf = (GEN)bnf[7]; r1 = nf_get_r1(nf);
-  polnf = (GEN)nf[1]; vnf = varn(polnf); degnf = degpol(polnf);
-  if (vnf==0) err(talker,"main variable in kummer must not be x");
-
-  p1 = conductor(bnr,all ? gzero : subgroup,2);
-  bnr = (GEN)p1[2];
-  if (!all) subgroup = (GEN)p1[3];
-  classgroup = (GEN)bnr[5];
-  cyclic = (GEN)classgroup[2];
-  bid = (GEN)bnr[2];
-  ideal= gmael(bid,1,1);
-  arch = gmael(bid,1,2); /* this is the conductor */
-  if (!all && gcmp0(subgroup)) subgroup = diagonal(cyclic);
-  gell = all? stoi(all): det(subgroup);
-  ell = itos(gell);
-
-  cyclicK= gmael3(bnf,8,1,2); rc = ellrank(cyclicK,ell);
-  genK   = gmael3(bnf,8,1,3); nbgenclK = lg(genK)-1;
-  listgamma0=cgetg(nbgenclK+1,t_VEC);
-  listgamma=cgetg(nbgenclK+1,t_VEC);
-  vecalpha0=cgetg(rc+1,t_VEC);
-  listalpha=cgetg(rc+1,t_VEC);
-  cycgen = check_and_build_cycgen(bnf);
-  p1 = gmul(gell,ideal);
-  for (i=1; i<=rc; i++)
+  r1 = nf_get_r1(nf);
+  if (!b)
   {
-    p3 = basistoalg(nf,idealcoprime(nf,(GEN)genK[i],p1));
-    p2 = basistoalg(nf, famat_to_nf(nf, (GEN)cycgen[i]));
-    listgamma[i] = listgamma0[i] = linv(p3);
-    vecalpha0[i] = (long)p2;
-    listalpha[i] = lmul((GEN)vecalpha0[i], powgi(p3,(GEN)cyclicK[i]));
+    if (typ(be) != t_COL) be = algtobasis(nf, be);
+    b = gmul(gmael(nf,5,1), be);
   }
-  for (   ; i<=nbgenclK; i++)
+  n = max((itos(ell)>>1), 3);
+  z = cgetg(n+1, t_VEC);
+  c = gmul(greal((GEN)bnfz[3]), ell);
+  c = logarch2arch(c, r1, prec); /* = embeddings of fu^ell */
+  c = gprec_w(gnorm(c), DEFAULTPREC);
+  b = gprec_w(gnorm(b), DEFAULTPREC); /* need little precision */
+  z[1] = (long)concatsp(c, vecinv(c));
+  for (k=2; k<=n; k++) z[k] = (long) vecmul((GEN)z[1], (GEN)z[k-1]);
+  nmax = T2_from_embed_norm(b, r1);
+  ru = lg(c)-1; c = zerovec(ru);
+  for(;;)
   {
-    long k;
-    p3 = basistoalg(nf,idealcoprime(nf,(GEN)genK[i],p1));
-    p2 = basistoalg(nf, famat_to_nf(nf, (GEN)cycgen[i]));
-    k = itos(mpinvmod((GEN)cyclicK[i], gell));
-    p2 = gpowgs(p2,k);
-    listgamma0[i]= (long)p2;
-    listgamma[i] = lmul(p2, gpowgs(p3, k * itos((GEN)cyclicK[i]) - 1));
-  }
-  fa = (GEN)bid[3];
-  prm   = (GEN)fa[1];
-  expom = (GEN)fa[2]; l = lg(prm);
-  Sm  = cgetg(l,t_VEC); setlg(Sm,1);
-  Sml1= cgetg(l,t_VEC); setlg(Sml1,1);
-  Sml2= cgetg(l,t_VEC); setlg(Sml2,1);
-  Sl  = cgetg(l+degnf,t_VEC); setlg(Sl,1);
-  ESml2=cgetg(l,t_VECSMALL); setlg(ESml2,1);
-  for (i=1; i<l; i++)
-  {
-    pr = (GEN)prm[i]; p = (GEN)pr[1]; e = itos((GEN)pr[3]);
-    vp = itos((GEN)expom[i]);
-    if (!egalii(p,gell))
-    {
-      if (vp != 1) return no_sol(all,1);
-      _append(Sm, pr);
-    }
-    else
-    {
-      vd = (vp-1)*(ell-1)-ell*e;
-      if (vd > 0) return no_sol(all,4);
-      if (vd==0)
-        _append(Sml1, pr);
-      else
+    GEN B = NULL;
+    long besti = 0, bestk = 0;
+    for (k=1; k<=n; k++)
+      for (i=1; i<=ru; i++)
       {
-	if (vp==1) return no_sol(all,2);
-        _append(Sml2, pr);
-        _append(ESml2,(GEN)vp);
+        p1 = vecmul(b, gmael(z,k,i));    p2 = T2_from_embed_norm(p1,r1);
+        if (gcmp(p2,nmax) < 0) { B=p1; nmax=p2; besti=i; bestk = k; continue; }
+        p1 = vecmul(b, gmael(z,k,i+ru)); p2 = T2_from_embed_norm(p1,r1);
+        if (gcmp(p2,nmax) < 0) { B=p1; nmax=p2; besti=i; bestk =-k; }
       }
-    }
+    if (!B) break;
+    b = B; c[besti] = laddis((GEN)c[besti], bestk);
   }
-  factell = primedec(nf,gell); l = lg(factell);
-  for (i=1; i<l; i++)
+  if (DEBUGLEVEL) fprintferr("naive reduction mod U^l: unit exp. = %Z\n",c);
+  return fix_be(bnfz, be, gmul(ell,c));
+}
+
+static GEN
+reduce_mod_Qell(GEN bnfz, GEN be, GEN gell)
+{
+  GEN c, fa;
+  if (typ(be) != t_COL) be = algtobasis(bnfz, be);
+  be = primitive_part(be, &c);
+  if (c)
   {
-    pr = (GEN)factell[i];
-    if (!idealval(nf,ideal,pr)) _append(Sl, pr);
+    fa = factor(c);
+    fa[2] = (long)FpV_red((GEN)fa[2], gell);
+    c = factorback(fa, NULL);
+    be = gmul(be, c);
   }
-  lSml2=lg(Sml2)-1; lSl=lg(Sl)-1;
-  Sp = concatsp(Sm,Sml1); lSp = lg(Sp)-1;
-  vecbeta = cgetg(lSp+1,t_VEC); matexpo=cgetg(lSp+1,t_MAT);
-  vecbeta0= cgetg(lSp+1,t_VEC);
-  for (j=1; j<=lSp; j++)
+  return be;
+}
+
+/* return q, q^n r = x, v_pr(r) < n for all pr. Insist q is a genuine n-th
+ * root (i.e r = 1) if strict != 0. */
+static GEN
+idealsqrtn(GEN nf, GEN x, GEN gn, int strict)
+{
+  long i, l, n = itos(gn);
+  GEN fa, q, Ex, Pr;
+
+  fa = idealfactor(nf, x);
+  Pr = (GEN)fa[1]; l = lg(Pr);
+  Ex = (GEN)fa[2]; q = NULL;
+  for (i=1; i<l; i++)
   {
-    p1=isprincipalgenforce(bnf,(GEN)Sp[j]);
-    p2=basistoalg(nf,(GEN)p1[2]);
-    for (i=1; i<=rc; i++)
-      p2 = gmul(p2, powgi((GEN)listgamma[i],gmael(p1,1,i)));
-    p3 = p2;
-    for (   ; i<=nbgenclK; i++)
-    {
-      p2 = gmul(p2, powgi((GEN)listgamma[i],gmael(p1,1,i)));
-      p3 = gmul(p3, powgi((GEN)listgamma0[i],gmael(p1,1,i)));
-    }
-    matexpo[j]=(long)p1[1];
-    vecbeta[j]=(long)p2; /* attention, ceci sont les beta modifies */
-    vecbeta0[j]=(long)p3;
+    long ex = itos((GEN)Ex[i]);
+    GEN e = stoi(ex / n);
+    if (strict && ex % n) err(talker,"not an n-th power in idealsqrtn");
+    if (q) q = idealmulpowprime(nf, q, (GEN)Pr[i], e);
+    else   q = idealpow(nf, (GEN)Pr[i], e);
   }
-  listg = gmodulcp(concatsp(gmael(bnf,8,5),gmael3(bnf,8,4,2)),polnf);
-  vunit = concatsp(listg, listalpha);
-  nbvunit=lg(vunit)-1;
-  id=idmat(degnf);
-  for (i=1; i<=lSl; i++)
+  return q? q: gun;
+}
+
+static GEN
+reducebeta(GEN bnfz, GEN be, GEN ell)
+{
+  long j,ru, prec = nfgetprec(bnfz);
+  GEN emb,z,u,matunit, nf = checknf(bnfz);
+
+  if (DEBUGLEVEL>1) fprintferr("reducing beta = %Z\n",be);
+  /* reduce mod Q^ell */
+  be = reduce_mod_Qell(nf, be, ell);
+  /* reduce l-th root */
+  z = idealsqrtn(nf, be, ell, 0);
+  if (typ(z) == t_MAT && !gcmp1(gcoeff(z,1,1)))
   {
-    pr = (GEN)Sl[i]; e = itos((GEN)pr[3]);
-    id = idealmul(nf, id, idealpows(nf,pr,(ell*e)/(ell-1)));
+    z = idealred_elt(nf, z);
+    be = element_div(nf, be, element_pow(nf, z, ell));
+    /* make be integral */
+    be = reduce_mod_Qell(nf, be, ell);
   }
-  vsml2=cgetg(lSml2+1,t_VEC);
-  for (i=1; i<=lSml2; i++)
+  if (DEBUGLEVEL>1) fprintferr("beta reduced via ell-th root = %Z\n",be);
+
+  matunit = gmul(greal((GEN)bnfz[3]), ell); /* log. embeddings of fu^ell */
+  for (;;)
   {
-    pr = (GEN)Sml2[i]; e = itos((GEN)pr[3]);
-    p1=idealpows(nf,pr, (e*ell)/(ell-1) + 1 - ESml2[i]);
-    id = idealmul(nf,id,p1);
-    p2 = idealmul(nf,p1,pr);
-    p3 = zidealstarinitgen(nf,p2);
-    vsml2[i] = (long)p3;
-    cycpro = gmael(p3,2,2); l = lg(cycpro)-1;
-    if (! gdivise((GEN)cycpro[l],gell)) err(talker,"bug5 in kummer");
+    z = get_arch_real(nf, be, &emb, prec);
+    if (z) break;
+    prec = (prec-1)<<1;
+    if (DEBUGLEVEL) err(warnprec,"reducebeta",prec);
+    nf = nfnewprec(nf,prec);
   }
-  bid = zidealstarinitgen(nf,id);
-  lastbid = ellrank(gmael(bid,2,2), ell);
-  nbcol=nbvunit+lSp; nblig=lastbid+rc;
-  munit=cgetg(nbcol+1,t_MAT); vecmunit2=cgetg(lSml2+1,t_VEC);
-  msign=cgetg(nbcol+1,t_MAT);
-  for (k=1; k<=lSml2; k++) vecmunit2[k]=lgetg(nbcol+1,t_MAT);
-  archartif=cgetg(r1+1,t_VEC); for (j=1; j<=r1; j++) archartif[j]=un;
-  for (j=1; j<=nbvunit; j++)
+  z = concatsp(matunit, z);
+  u = lllintern(z, 100, 1, prec);
+  if (u)
   {
-    p1 = zideallog(nf,(GEN)vunit[j],bid);
-    p2 = cgetg(nblig+1,t_COL); munit[j]=(long)p2;
-    for (i=1; i<=lastbid; i++) p2[i]=p1[i];
-    for (   ; i<=nblig; i++) p2[i]=zero;
-    for (k=1; k<=lSml2; k++)
-      mael(vecmunit2,k,j) = (long)zideallog(nf,(GEN)vunit[j],(GEN)vsml2[k]);
-    msign[j] = (long)zsigne(nf,(GEN)vunit[j],archartif);
+    ru = lg(u);
+    for (j=1; j < ru; j++)
+      if (gcmp1(gcoeff(u,ru-1,j))) break;
+    if (j < ru)
+    {
+      u = (GEN)u[j]; /* coords on (fu^ell, be) of a small generator */
+      ru--; setlg(u, ru);
+      be = fix_be(bnfz, be, gmul(ell,u));
+    }
   }
-  for (j=1; j<=lSp; j++)
-  {
-    p1 = zideallog(nf,(GEN)vecbeta[j],bid);
-    p2 = cgetg(nblig+1,t_COL); munit[j+nbvunit]=(long)p2;
-    for (i=1; i<=lastbid; i++) p2[i]=p1[i];
-    for (   ; i<=nblig; i++) p2[i]=coeff(matexpo,i-lastbid,j);
-    for (k=1; k<=lSml2; k++)
-      mael(vecmunit2,k,j+nbvunit) = (long)zideallog(nf,(GEN)vecbeta[j],(GEN)vsml2[k]);
-    msign[j+nbvunit] = (long)zsigne(nf,(GEN)vecbeta[j],archartif);
-  }
-  listx = get_listx(arch,msign,munit,vecmunit2,ell,lSp,nbvunit);
-  llistx= lg(listx);
-  listal= cgetg(llistx,t_VEC);
-  listg = concatsp(listg, concatsp(vecalpha0,vecbeta0));
-  l = lg(listg);
-  for (i=1; i<llistx; i++)
-  {
-    p1 = gun; p2 = (GEN)listx[i];
-    for (j=1; j<l; j++)
-      p1 = gmul(p1, powgi((GEN)listg[j],(GEN)p2[j]));
-    listal[i] = (long)reducealpha(nf,p1,gell);
-  }
- /* A ce stade, tous les alpha dans la liste listal statisfont les
-  * congruences, noncongruences et signatures, et ne sont pas des puissances
-  * l-iemes donc x^l-alpha est irreductible, sa signature est correcte ainsi
-  * que son discriminant relatif. Reste a determiner son groupe de normes. */
-  if (DEBUGLEVEL) fprintferr("listalpha = %Z\n",listal);
-  p2 = cgetg(1,t_VEC);
-  for (i=1; i<llistx; i++)
-  {
-    p1 = gsub(gpuigs(polx[0],ell), (GEN)listal[i]);
-    if (all || gegal(rnfnormgroup(bnr,p1),subgroup)) p2 = concatsp(p2,p1);
-  }
-  if (all) return gcopy(p2);
-  switch(lg(p2)-1)
-  {
-    case 0: err(talker,"bug 6: no equation found in kummer");
-    case 1: break; /* OK */
-    default: fprintferr("equations = \n%Z",p2);
-      err(talker,"bug 7: more than one equation found in kummer");
-  }
-  return gcopy((GEN)p2[1]);
+  if (DEBUGLEVEL>1) fprintferr("beta LLL-reduced mod U^l = %Z\n",be);
+  return reducebetanaive(bnfz, be, NULL, ell);
 }
 
 static GEN
-tauofideal(GEN id)
-{
-  long j;
-  GEN p1,p2;
-
-  p1=gsubst(gmul((GEN)nfz[7],id),vnf,U);
-  p2=cgetg(lg(p1),t_MAT);
-  for (j=1; j<lg(p1); j++) p2[j]=(long)algtobasis(nfz,(GEN)p1[j]);
-  return p2;
+tauofalg(GEN x, GEN U) {
+  return gsubst(lift(x), varn(U[1]), U);
 }
 
-static GEN
-tauofprimeideal(GEN pr)
+static tau_s *
+get_tau(tau_s *tau, GEN nf, GEN U)
 {
-  GEN p1 = dummycopy(pr);
-
-  p1[2] = (long)algtobasis(nfz, gsubst(gmul((GEN)nfz[7],(GEN)pr[2]),vnf,U));
-  return gcoeff(idealfactor(nfz,prime_to_ideal(nfz,p1)),1,1);
+  GEN bas = (GEN)nf[7], Uzk;
+  long i, l = lg(bas);
+  Uzk = cgetg(l, t_MAT);
+  for (i=1; i<l; i++)
+    Uzk[i] = (long)algtobasis(nf, tauofalg((GEN)bas[i], U));
+  tau->zk = Uzk;
+  tau->x  = U; return tau;
 }
 
-static long
-isprimeidealconj(GEN pr1, GEN pr2)
-{
-  GEN pr=pr1;
+static GEN tauoffamat(GEN x, tau_s *tau);
 
-  do
+static GEN
+tauofelt(GEN x, tau_s *tau)
+{
+  switch(typ(x))
   {
-    if (gegal(pr,pr2)) return 1;
-    pr = tauofprimeideal(pr);
+    case t_COL: return gmul(tau->zk, x);
+    case t_MAT: return tauoffamat(x, tau);
+    default: return tauofalg(x, tau->x);
   }
-  while (!gegal(pr,pr1));
-  return 0;
 }
+static GEN
+tauofvec(GEN x, tau_s *tau)
+{
+  long i, l = lg(x);
+  GEN y = cgetg(l, typ(x));
 
-static long
-isconjinprimelist(GEN listpr, GEN pr2)
+  for (i=1; i<l; i++) y[i] = (long)tauofelt((GEN)x[i], tau);
+  return y;
+}
+/* [x, tau(x), ..., tau^m(x)] */
+static GEN
+powtau(GEN x, long m, tau_s *tau)
 {
-  long ll=lg(listpr)-1,i;
+  GEN y = cgetg(m+1, t_VEC);
+  long i;
+  y[1] = (long)x;
+  for (i=2; i<=m; i++) y[i] = (long)tauofelt((GEN)y[i-1], tau);
+  return y;
+}
 
-  for (i=1; i<=ll; i++)
-    if (isprimeidealconj((GEN)listpr[i],pr2)) return 1;
-  return 0;
+static GEN
+tauoffamat(GEN x, tau_s *tau)
+{
+  GEN y = cgetg(3, t_MAT);
+  y[1] = (long)tauofvec((GEN)x[1], tau);
+  y[2] = x[2]; return y;
 }
 
-static void
-computematexpoteta1(GEN A1, GEN R)
+/* TODO: wasteful to reduce to 2-elt form. Compute image directly */
+static GEN
+tauofideal(GEN nfz, GEN id, tau_s *tau)
 {
-  long j;
-  GEN Aj = polun[vnf];
-  matexpoteta1 = cgetg(degK+1,t_MAT);
-  for (j=1; j<=degK; j++)
+  GEN I = ideal_two_elt(nfz,id);
+  I[2] = (long)tauofelt((GEN)I[2], tau);
+  return prime_to_ideal(nfz,I);
+}
+
+static int
+isprimeidealconj(GEN nfz, GEN pr1, GEN pr2, tau_s *tau)
+{
+  GEN p = (GEN)pr1[1];
+  GEN x = (GEN)pr1[2];
+  GEN b1= (GEN)pr1[5];
+  GEN b2= (GEN)pr2[5];
+  if (!egalii(p, (GEN)pr2[1])
+   || !egalii((GEN)pr1[3], (GEN)pr2[3])
+   || !egalii((GEN)pr1[4], (GEN)pr2[4])) return 0;
+  if (gegal(x,(GEN)pr2[2])) return 1;
+  for(;;)
   {
-    matexpoteta1[j] = (long)pol_to_vec(Aj, degKz);
-    if (j<degK) Aj = gmod(gmul(Aj,A1), R);
+    if (int_elt_val(nfz,x,p,b2,NULL)) return 1;
+    x = FpV_red(tauofelt(x, tau), p);
+    if (int_elt_val(nfz,x,p,b1,NULL)) return 0;
   }
 }
 
-static GEN
-downtoK(GEN x)
+static int
+isconjinprimelist(GEN nfz, GEN S, GEN pr, tau_s *tau)
 {
-  long i;
-  GEN p2,p3;
+  long i, l;
+  
+  if (!tau) return 0;
+  l = lg(S);
+  for (i=1; i<l; i++)
+    if (isprimeidealconj(nfz, (GEN)S[i],pr,tau)) return 1;
+  return 0;
+}
 
-  p2 = inverseimage(matexpoteta1, pol_to_vec(lift(x), degKz));
-  if (lg(p2)==1) err(talker,"not an element of K in downtoK");
-  p3 = (GEN)p2[degK];
-  for (i=degK-1; i; i--) p3 = gadd((GEN)p2[i],gmul(polx[vnf],p3));
-  return gmodulcp(p3,polnf);
+static GEN
+downtoK(toK_s *T, GEN x)
+{
+  long degKz = lg(T->invexpoteta1) - 1;
+  GEN y = gmul(T->invexpoteta1, pol_to_vec(lift(x), degKz));
+  return gmodulcp(gtopolyrev(y,varn(T->polnf)), T->polnf);
 }
 
 static GEN
-tracetoK(GEN x)
+tracetoK(toK_s *T, GEN x)
 {
+  GEN p1 = x;
   long i;
-  GEN p1,p2;
-
-  p1=x; p2=x;
-  for (i=1; i<=m-1; i++)
-  {
-    p1=gsubst(lift(p1),vnf,U);
-    p2=gadd(p2,p1);
-  }
-  return downtoK(p2);
+  for (i=1; i < T->m; i++) p1 = gadd(x, tauofelt(p1,T->tau));
+  return downtoK(T, p1);
 }
 
 static GEN
-normtoK(GEN x)
+normtoK(toK_s *T, GEN x)
 {
+  GEN p1 = x;
   long i;
-  GEN p1,p2;
-
-  p1=x; p2=x;
-  for (i=1; i<=m-1; i++)
-  {
-    p1=gsubst(lift(p1),vnf,U);
-    p2=gmul(p2,p1);
-  }
-  return downtoK(p2);
+  for (i=1; i < T->m; i++) p1 = gmul(x, tauofelt(p1,T->tau));
+  return downtoK(T, p1);
 }
 
 static GEN
-computepolrel(void)
+no_sol(long all, int i)
 {
-  long i,j;
-  GEN p1,p2;
-
-  p1=gun; p2=gmodulcp(polx[vnf],R);
-  for (i=0; i<=m-1; i++)
-  {
-    p1=gmul(p1,gsub(polx[0],p2));
-    if (i<m-1) p2=gsubst(lift(p2),vnf,U);
-  }
-  for (j=2; j<=m+2; j++) p1[j]=(long)downtoK((GEN)p1[j]);
-  return p1;
+  if (!all) err(talker,"bug%d in kummer",i);
+  return cgetg(1,t_VEC);
 }
 
-/* alg. 5.2.15. with remark */
 static GEN
-isprincipalell(GEN id)
+get_gell(GEN bnr, GEN subgp, long all)
 {
-  long i, l = lg(vecalpha);
-  GEN y,logdisc,be;
-
-  y = isprincipalgenforce(bnfz,id);
-  logdisc = (GEN)y[1];
-  be = basistoalg(bnfz,(GEN)y[2]);
-  for (i=rc+1; i<l; i++)
-    be = gmul(be, powgi((GEN)vecalpha[i], modii(mulii((GEN)logdisc[i],(GEN)uu[i]),gell)));
-  y = cgetg(3,t_VEC);
-  y[1]=(long)logdisc; setlg(logdisc,rc+1);
-  y[2]=(long)be;
-  return y;
+  GEN gell;
+  if (all)        gell = stoi(all);
+  else if (subgp) gell = det(subgp);
+  else            gell = det(diagonal(gmael(bnr,5,2)));
+  if (typ(gell) != t_INT) err(arither1);
+  return gell;
 }
 
-/* alg. 5.3.11. */
-static GEN
-isvirtualunit(GEN v)
+typedef struct {
+  GEN Sm, Sml1, Sml2, Sl, ESml2;
+} primlist;
+
+static int
+build_list_Hecke(primlist *L, GEN nfz, GEN fa, GEN gothf, GEN gell, tau_s *tau)
 {
-  long llist,i,j,ex, l = lg(vecalpha);
-  GEN p1,listex,listpr,q,ga,eps,vecy,logdisc;
+  GEN listpr, listex, pr, p, factell;
+  long vd, vp, e, i, l, ell = itos(gell), degKz = degpol(nfz[1]);
 
-  p1=idealfactor(nfz,v);
-  listpr = (GEN)p1[1]; llist = lg(listpr)-1;
-  listex = (GEN)p1[2]; q = idmat(degKz);
-  for (i=1; i<=llist; i++)
+  if (!fa) fa = idealfactor(nfz, gothf);
+  listpr = (GEN)fa[1];
+  listex = (GEN)fa[2]; l = lg(listpr);
+  L->Sm  = cget1(l,t_VEC);
+  L->Sml1= cget1(l,t_VEC);
+  L->Sml2= cget1(l,t_VEC);
+  L->Sl  = cget1(l+degKz,t_VEC);
+  L->ESml2=cget1(l,t_VECSMALL);
+  for (i=1; i<l; i++)
   {
-    ex = itos((GEN)listex[i]);
-    if (ex%ell) err(talker,"not a virtual unit in isvirtualunit");
-    q = idealmul(nfz,q, idealpows(nfz,(GEN)listpr[i], ex/ell));
+    pr = (GEN)listpr[i]; p = (GEN)pr[1]; e = itos((GEN)pr[3]);
+    vp = itos((GEN)listex[i]);
+    if (!egalii(p,gell))
+    {
+      if (vp != 1) return 1;
+      if (!isconjinprimelist(nfz, L->Sm,pr,tau)) appendL(L->Sm,pr);
+    }
+    else
+    {
+      vd = (vp-1)*(ell-1)-ell*e;
+      if (vd > 0) return 4;
+      if (vd==0)
+      {
+	if (!isconjinprimelist(nfz, L->Sml1,pr,tau)) appendL(L->Sml1, pr);
+      }
+      else
+      {
+	if (vp==1) return 2;
+        if (!isconjinprimelist(nfz, L->Sml2,pr,tau))
+        {
+          appendL(L->Sml2, pr);
+          appendL(L->ESml2,(GEN)vp);
+        }
+      }
+    }
   }
-  /* q^ell = (v) */
-  p1 = isprincipalgenforce(bnfz,q);
-  logdisc = (GEN)p1[1];
-  ga = basistoalg(nfz,(GEN)p1[2]);
-  for (j=rc+1; j<l; j++)
-    ga = gmul(ga, powgi((GEN)vecalpha[j],divii((GEN)logdisc[j],(GEN)cyc[j])));
-  eps = gpuigs(ga,ell);
-  vecy = cgetg(rc+1,t_COL);
-  for (j=1; j<=rc; j++)
+  factell = primedec(nfz,gell); l = lg(factell);
+  for (i=1; i<l; i++)
   {
-    vecy[j] = (long)divii((GEN)logdisc[j], divii((GEN)cyc[j],gell));
-    eps = gmul(eps, powgi((GEN)vecalpha[j],(GEN)vecy[j]));
+    pr = (GEN)factell[i];
+    if (!idealval(nfz,gothf,pr))
+      if (!isconjinprimelist(nfz, L->Sl,pr,tau)) appendL(L->Sl, pr);
   }
-  eps = gdiv(v,eps);
-  p1 = cgetg(3,t_VEC);
-  p1[1] = (long)concatsp(vecy, lift(isunit(bnfz,eps)));
-  p1[2] = (long)ga;
-  return p1;
+  return 0; /* OK */
 }
 
 static GEN
-lifttokz(GEN id, GEN A1)
+logall(GEN nf, GEN vec, long lW, long mginv, long ell, GEN pr, long ex)
 {
-  long i,j;
-  GEN p1,p2,p3;
+  GEN m, M, bid = zidealstarinitgen(nf, idealpows(nf, pr, ex));
+  long ellrank, i, l = lg(vec);
 
-  p1=gsubst(gmul((GEN)nf[7],id),vnf,A1);
-  p2=gmodulcp((GEN)nfz[7],R);
-  p3=cgetg(degK*degKz+1,t_MAT);
-  for (i=1; i<=degK; i++)
-    for (j=1; j<=degKz; j++)
-      p3[(i-1)*degKz+j]=(long)algtobasis(nfz,gmul((GEN)p1[i],(GEN)p2[j]));
-  return hnfmod(p3,detint(p3));
-}
-
-static GEN
-steinitzaux(GEN id, GEN polrel)
-{
-  long i,j;
-  GEN p1,p2,p3,vecid,matid,pseudomat,pid;
-
-  p1=gsubst(gmul((GEN)nfz[7],id),vnf,polx[0]);
-  for (j=1; j<=degKz; j++)
-    p1[j]=(long)gmod((GEN)p1[j],polrel);
-  p2=cgetg(degKz+1,t_MAT);
-  for (j=1; j<=degKz; j++)
+  ellrank = prank(gmael(bid,2,2), ell);
+  M = cgetg(l,t_MAT);
+  for (i=1; i<l; i++)
   {
-    p3=cgetg(m+1,t_COL); p2[j]=(long)p3;
-    for (i=1; i<=m; i++) p3[i]=(long)algtobasis(nf,truecoeff((GEN)p1[j],i-1));
+    m = zideallog(nf, (GEN)vec[i], bid);
+    setlg(m, ellrank+1);
+    if (i < lW) m = gmulsg(mginv, m);
+    M[i] = (long)m;
   }
-  vecid=cgetg(degKz+1,t_VEC); matid=idmat(degK);
-  for (j=1; j<=degKz; j++) vecid[j]=(long)matid;
-  pseudomat=cgetg(3,t_VEC);
-  pseudomat[1]=(long)p2; pseudomat[2]=(long)vecid;
-  pid=(GEN)nfhermite(nf,pseudomat)[2];
-  p1=matid;
-  for (j=1; j<=m; j++) p1=idealmul(nf,p1,(GEN)pid[j]);
-  return p1;
+  return M;
 }
 
+/* compute the u_j (see remark 5.2.15.) */
 static GEN
-normrelz(GEN id, GEN polrel, GEN steinitzZk)
+get_u(GEN cyc, long rc, GEN gell)
 {
-  GEN p1 = steinitzaux(idealhermite(nfz, id), polrel);
-  return idealdiv(nf,p1,steinitzZk);
+  long i, l = lg(cyc);
+  GEN u = cgetg(l,t_VEC);
+  for (i=1; i<=rc; i++) u[i] = zero;
+  for (   ; i<  l; i++) u[i] = lmpinvmod((GEN)cyc[i], gell);
+  return u;
 }
 
+/* alg. 5.2.15. with remark */
 static GEN
-invimsubgroup(GEN bnrz, GEN bnr, GEN subgroup)
+isprincipalell(GEN bnfz, GEN id, GEN cycgen, GEN u, GEN gell, long rc)
 {
-  long l,j;
-  GEN g,Plog,raycycz,rayclgpz,genraycycz,U,polrel,steinitzZk;
+  long i, l = lg(cycgen);
+  GEN logdisc, b, y = quick_isprincipalgen(bnfz, id);
 
-  polrel = computepolrel();
-  steinitzZk = steinitzaux(idmat(degKz), polrel); 
-  rayclgpz = (GEN)bnrz[5];
-  raycycz   = (GEN)rayclgpz[2]; l=lg(raycycz);
-  genraycycz= (GEN)rayclgpz[3];
-  Plog = cgetg(l,t_MAT);
-  for (j=1; j<l; j++)
+  logdisc = gmod((GEN)y[1], gell);
+  b = (GEN)y[2];
+  for (i=rc+1; i<l; i++)
   {
-    g = normrelz((GEN)genraycycz[j],polrel,steinitzZk);
-    Plog[j] = (long)isprincipalray(bnr, g);
+    GEN e = modii(mulii((GEN)logdisc[i],(GEN)u[i]), gell);
+    if (signe(e)) b = famat_mul(b, famat_pow((GEN)cycgen[i], e));
   }
-  U = (GEN)hnfall(concatsp(Plog, subgroup))[2];
-  setlg(U, l); for (j=1; j<l; j++) setlg(U[j], l);
-  return hnfmod(concatsp(U, diagonal(raycycz)), (GEN)raycycz[1]);
+  y = cgetg(3,t_VEC);
+  y[1] = (long)logdisc; setlg(logdisc,rc+1);
+  y[2] = (long)b; return y;
 }
 
 static GEN
-ideallogaux(long i, GEN al)
+famat_factorback(GEN v, GEN e)
 {
-  long llogli,valal;
-  GEN p1;
-
-  valal = element_val(nfz,al,(GEN)listprSp[i]);
-  al = gmul(al,gpuigs((GEN)listunif[i],valal));
-  p1 = zideallog(nfz,al,(GEN)listbid[i]);
-  llogli = listellrank[i];
-  setlg(p1,llogli+1); return p1;
+  long i, l = lg(e);
+  GEN V = cgetg(1, t_MAT);
+  for (i=1; i<l; i++) 
+    if (signe(e[i])) V = famat_mul(V, famat_pow((GEN)v[i], (GEN)e[i]));
+  return V;
 }
 
 static GEN
-ideallogauxsup(long i, GEN al)
+compute_beta(GEN X, GEN vecWB, GEN ell, GEN bnfz)
 {
-  long llogli,valal;
-  GEN p1;
+  GEN BE, be;
+  BE = famat_reduce(famat_factorback(vecWB, X));
+  BE[2] = (long)centermod((GEN)BE[2], ell);
+  be = factorbackelt(BE, bnfz, NULL);
+  be = reducebeta(bnfz, be, ell);
+  if (DEBUGLEVEL>1) fprintferr("beta reduced = %Z\n",be);
+  return basistoalg(bnfz, be); /* FIXME */
+}
 
-  valal = element_val(nfz,al,(GEN)listprSp[i]);
-  al = gmul(al,gpuigs((GEN)listunif[i],valal));
-  p1 = zideallog(nfz,al,(GEN)listbidsup[i]);
-  llogli = listellranksup[i];
-  setlg(p1,llogli+1); return p1;
+static GEN
+get_Selmer(GEN bnf, GEN cycgen, long rc)
+{
+  GEN fu = check_units(bnf,"rnfkummer");
+  GEN tu = gmael3(bnf,8,4,2);
+  return concatsp(algtobasis(bnf,concatsp(fu,tu)), vecextract_i(cycgen,1,rc));
 }
 
+/* if all!=0, give all equations of degree 'all'. Assume bnr modulus is the
+ * conductor */
 static GEN
-vectau(GEN list)
+rnfkummersimple(GEN bnr, GEN subgroup, GEN gell, long all)
 {
-  long i,j,k,n;
-  GEN listz,listc,yz,yc,listfl,s, y = cgetg(3,t_VEC);
+  long ell, i, j, degK, dK;
+  long lSml2, lSl2, lSp, rc, lW;
+  long prec;
 
-  listz = (GEN)list[1];
-  listc = (GEN)list[2]; n = lg(listz);
-  yz = cgetg(n,t_VEC); y[1] = (long)yz;
-  yc = cgetg(n,t_VEC); y[2] = (long)yc;
-  listfl=cgetg(n,t_VECSMALL); for (i=1; i<n; i++) listfl[i] = 1;
-  k = 1;
-  for (i=1; i<n; i++)
+  GEN bnf,nf,bid,ideal,arch,cycgen;
+  GEN clgp,cyc;
+  GEN Sp,listprSp,matP;
+  GEN res,u,M,K,y,vecMsup,vecW,vecWB,vecBp,msign;
+  primlist L;
+
+  bnf = (GEN)bnr[1];
+  nf  = (GEN)bnf[7];
+  degK = degpol(nf[1]);
+  
+  bid = (GEN)bnr[2];
+  ideal= gmael(bid,1,1);
+  arch = gmael(bid,1,2); /* this is the conductor */
+  ell = itos(gell);
+  i = build_list_Hecke(&L, nf, (GEN)bid[3], ideal, gell, NULL);
+  if (i) return no_sol(all,i);
+
+  lSml2 = lg(L.Sml2)-1;
+  Sp = concatsp(L.Sm, L.Sml1); lSp = lg(Sp)-1;
+  listprSp = concatsp(L.Sml2, L.Sl); lSl2 = lg(listprSp)-1;
+
+  cycgen = check_and_build_cycgen(bnf);
+  clgp = gmael(bnf,8,1);
+  cyc = (GEN)clgp[2]; rc = prank(cyc, ell);
+
+  vecW = get_Selmer(bnf, cycgen, rc);
+  u = get_u(cyc, rc, gell);
+
+  vecBp = cgetg(lSp+1, t_VEC);
+  matP  = cgetg(lSp+1, t_MAT);
+  for (j=1; j<=lSp; j++)
   {
-    if (!listfl[i]) continue;
+    GEN e, a, L;
+    L = isprincipalell(bnf,(GEN)Sp[j], cycgen,u,gell,rc);
+    e = (GEN)L[1]; matP[j]  = (long)e;
+    a = (GEN)L[2]; vecBp[j] = (long)a;
+  }
+  vecWB = concatsp(vecW, vecBp);
 
-    yz[k] = listz[i];
-    s = (GEN)listc[i];
-    for (j=i+1; j<n; j++)
+  prec = DEFAULTPREC +
+    ((gexpo(vecWB) + gexpo(gmael(nf,5,1))) >> TWOPOTBYTES_IN_LONG);
+  if (nfgetprec(nf) < prec) nf = nfnewprec(nf, prec);
+  msign = zsigns(nf, vecWB);
+
+  vecMsup = cgetg(lSml2+1,t_VEC);
+  M = NULL;
+  for (i=1; i<=lSl2; i++)
+  {
+    GEN pr = (GEN)listprSp[i];
+    long e = itos((GEN)pr[3]), z = ell * (e / (ell-1));
+
+    if (i <= lSml2)
     {
-      if (listfl[j] && gegal((GEN)listz[j],(GEN)listz[i]))
-      {
-        s = gadd(s,(GEN)listc[j]);
-        listfl[j] = 0;
-      }
+      z += 1 - L.ESml2[i];
+      vecMsup[i] = (long)logall(nf, vecWB, 0,0, ell, pr,z+1);
     }
-    yc[k] = (long)s; k++;
+    M = vconcat(M, logall(nf, vecWB, 0,0, ell, pr,z));
   }
-  setlg(yz, k);
-  setlg(yc, k); return y;
-}
+  lW = lg(vecW);
+  M = vconcat(M, concatsp(zeromat(rc,lW-1), matP));
 
-static GEN
-subtau(GEN listx, GEN listy)
-{
-  GEN y = cgetg(3,t_VEC);
-  y[1] = (long)concatsp((GEN)listx[1], (GEN)listy[1]);
-  y[2] = (long)concatsp((GEN)listx[2], gneg_i((GEN)listy[2]));
-  return vectau(y);
+  K = FpM_ker(M, gell);
+  dK = lg(K)-1;
+  y = cgetg(dK+1,t_VECSMALL);
+  res = cgetg(1,t_VEC); /* in case all = 1 */
+  while (dK)
+  {
+    for (i=1; i<dK; i++) y[i] = 0;
+    y[i] = 1; /* y = [0,...,0,1,0,...,0], 1 at dK'th position */
+    do
+    {
+      GEN be, P, X = FpV_red(gmul_mati_smallvec(K, y), gell);
+      if (ok_congruence(X, gell, lW, vecMsup) && ok_sign(X, msign, arch)) 
+      {/* be satisfies all congruences, x^ell - be is irreducible, signature
+        * and relative discriminant are correct */
+        be = compute_beta(X, vecWB, gell, bnf);
+        P = gsub(gpowgs(polx[0],ell), be);
+        if (!all && gegal(rnfnormgroup(bnr,P),subgroup)) return P; /*DONE*/
+        res = concatsp(res, P);
+      }
+    } while (increment(y, dK, ell));
+    y[dK--] = 0;
+  }
+  if (all) return res;
+  return gzero;
 }
 
+/* alg. 5.3.11 (return only discrete log mod ell) */
 static GEN
-negtau(GEN list)
+isvirtualunit(GEN bnf, GEN v, GEN cycgen, GEN cyc, GEN gell, long rc)
 {
-  GEN y = cgetg(3,t_VEC);
-  y[1] = list[1];
-  y[2] = lneg((GEN)list[2]);
-  return y;
+  GEN L, b, eps, y, q, nf = checknf(bnf);
+  long i, l = lg(cycgen);
+
+  L = quick_isprincipalgen(bnf, idealsqrtn(nf, v, gell, 1));
+  q = (GEN)L[1];
+  if (gcmp0(q)) { eps = v; y = q; }
+  else
+  {
+    b = (GEN)L[2];
+    y = cgetg(l,t_COL);
+    for (i=1; i<l; i++)
+      y[i] = (long)diviiexact(mulii(gell,(GEN)q[i]), (GEN)cyc[i]);
+    eps = famat_mul(famat_factorback(cycgen, y), famat_pow(b, gell));
+    eps = famat_mul(famat_inv(eps), v);
+  }
+  setlg(y, rc+1);
+  b = isunit(bnf,eps);
+  if (lg(b) == 1) err(bugparier,"isvirtualunit");
+  return concatsp(lift_intern(b), y);
 }
 
+/* id a vector of elements in nfz = relative extension of nf by polrel,
+ * return the Steinitz element associated to the module generated by id */
 static GEN
-multau(GEN listx, GEN listy)
+steinitzaux(GEN nf, GEN id, GEN polrel)
 {
-  GEN lzx,lzy,lcx,lcy,lzz,lcz, y = cgetg(3,t_VEC);
-  long nx,ny,i,j,k;
+  long i, degKz = lg(id)-1, degK = degpol(nf[1]);
+  GEN V = dummycopy(id),vecid,matid,pseudomat,pid;
 
-  lzx=(GEN)listx[1]; lcx=(GEN)listx[2]; nx=lg(lzx)-1;
-  lzy=(GEN)listy[1]; lcy=(GEN)listy[2]; ny=lg(lzy)-1;
-  lzz = cgetg(nx*ny+1,t_VEC); y[1]=(long)lzz;
-  lcz = cgetg(nx*ny+1,t_VEC); y[2]=(long)lcz;
-  k = 0;
-  for (i=1; i<=nx; i++)
-    for (j=1; j<=ny; j++)
-    {
-      k++;
-      lzz[k] = ladd((GEN)lzx[i],(GEN)lzy[j]);
-      lcz[k] = lmul((GEN)lcx[i],(GEN)lcy[j]);
-    }
-  return vectau(y);
+  vecid = cgetg(degKz+1,t_VEC); matid = idmat(degK);
+  for (i=1; i<=degKz; i++) vecid[i] = (long)matid;
+  for (i=1; i<=degKz; i++)
+  {
+    GEN v = (GEN)V[i];
+    if (typ(v) == t_POL) { v = dummycopy(v); setvarn(v, 0); }
+    V[i] = (long)gmod(v, polrel);
+  }
+  pseudomat = cgetg(3,t_VEC);
+  pseudomat[1] = (long)vecpol_to_mat(V, degpol(polrel));
+  pseudomat[2] = (long)vecid;
+  pid = (GEN)nfhermite(nf,pseudomat)[2];
+  return factorback(pid, nf); /* product */
 }
 
 static GEN
-mulpoltau(GEN poltau, GEN list)
+polrelKzK(toK_s *T, GEN x)
 {
-  long i,j;
-  GEN y;
-
-  j = lg(poltau)-2;
-  y = cgetg(j+3,t_VEC);
-  y[1] = (long)negtau(multau(list,(GEN)poltau[1]));
-  for (i=2; i<=j+1; i++)
-    y[i] = (long)subtau((GEN)poltau[i-1],multau(list,(GEN)poltau[i]));
-  y[j+2] = poltau[j+1]; return y;
+  GEN P = roots_to_pol(powtau(x, T->m, T->tau), 0);
+  long i, l = lg(P);
+  for (i=2; i<l; i++) P[i] = (long)downtoK(T, (GEN)P[i]);
+  return P;
 }
 
-/* th. 5.3.5. and prop. 5.3.9. */
+/* N: Cl_m(Kz) --> Cl_m(K), lift subgroup from bnr to bnrz using Algo 4.1.11 */
 static GEN
-computepolrelbeta(GEN be)
+invimsubgroup(GEN bnrz, GEN bnr, GEN subgroup, toK_s *T)
 {
-  long i,a,b,j;
-  GEN e,u,u1,u2,u3,p1,p2,p3,p4,zet,be1,be2,listr,s,veczi,vecci,powtaubet;
+  long l, j;
+  GEN P,raycycz,rayclgpz,raygenz,U,polrel,steinitzZk;
+  GEN nf = checknf(bnr), nfz = checknf(bnrz), polz = (GEN)nfz[1];
 
-  switch (ell)
+  polrel = polrelKzK(T, polx[varn(polz)]);
+  steinitzZk = steinitzaux(nf, (GEN)nfz[7], polrel); 
+  rayclgpz = (GEN)bnrz[5];
+  raycycz = (GEN)rayclgpz[2]; l = lg(raycycz);
+  raygenz = (GEN)rayclgpz[3];
+  P = cgetg(l,t_MAT);
+  for (j=1; j<l; j++)
   {
-    case 2: err(talker,"you should not be here in rnfkummer !!"); break;
-    case 3: e=normtoK(be); u=tracetoK(be);
-      return gsub(gmul(polx[0],gsub(gsqr(polx[0]),gmulsg(3,e))),gmul(e,u));
-    case 5: e=normtoK(be);
-      if (d==2)
-      {
-	u=tracetoK(gpuigs(be,3));
-	p1=gadd(gmulsg(5,gsqr(e)), gmul(gsqr(polx[0]), gsub(gsqr(polx[0]),gmulsg(5,e))));
-	return gsub(gmul(polx[0],p1),gmul(e,u));
-      }
-      else
-      {
-	be1=gsubst(lift(be),vnf,U);
-        be2=gsubst(lift(be1),vnf,U);
-	u1=tracetoK(gmul(be,be1));
-        u2=tracetoK(gmul(gmul(be,be2),gsqr(be1)));
-	u3=tracetoK(gmul(gmul(gsqr(be),gpuigs(be1,3)),be2));
-	p1=gsub(gsqr(polx[0]),gmulsg(10,e));
-	p1=gsub(gmul(polx[0],p1),gmulsg(5,gmul(e,u1)));
-	p1=gadd(gmul(polx[0],p1),gmul(gmulsg(5,e),gsub(e,u2)));
-	p1=gsub(gmul(polx[0],p1),gmul(e,u3));
-	return p1;
-      }
-    default: p1=cgetg(2,t_VEC); p2=cgetg(3,t_VEC); p3=cgetg(2,t_VEC);
-      p4=cgetg(2,t_VEC); p3[1]=zero; p4[1]=un;
-      p2[1]=(long)p3; p2[2]=(long)p4; p1[1]=(long)p2;
-      zet=gmodulcp(polx[vnf], cyclo(ell,vnf));
-      listr=cgetg(m+1,t_VEC);
-      listr[1]=un;
-      for (i=2; i<=m; i++) listr[i]=(long)modii(mulii((GEN)listr[i-1],g),gell);
-      veczi=cgetg(m+1,t_VEC);
-      for (b=0; b<m; b++)
-      {
-	s=gzero;
-	for (a=0; a<m; a++)
-	  s=gadd(gmul(polx[0],s),modii(mulii((GEN)listr[b+1],(GEN)listr[a+1]),gell));
-	veczi[b+1]=(long)s;
-      }
-      for (j=0; j<ell; j++)
-      {
-	vecci=cgetg(m+1,t_VEC);
-	for (b=0; b<m; b++) vecci[b+1]=(long)gpui(zet,mulsi(j,(GEN)listr[b+1]),0);
-	p4=cgetg(3,t_VEC);
-	p4[1]=(long)veczi; p4[2]=(long)vecci;
-	p1=mulpoltau(p1,p4);
-      }
-      powtaubet=cgetg(m+1,t_VEC);
-      powtaubet[1]=(long)be;
-      for (i=2; i<=m; i++) powtaubet[i]=(long)gsubst(lift((GEN)powtaubet[i-1]),vnf,U);
-      err(impl,"difficult Kummer for ell>=7"); return gzero;
+    GEN g, id = idealhermite(nfz, (GEN)raygenz[j]);
+    g = steinitzaux(nf, gmul((GEN)nfz[7], id), polrel);
+    g = idealdiv(nf, g, steinitzZk); /* N_{Kz/K}(gen[j]) */
+    P[j] = (long)isprincipalray(bnr, g);
   }
-  return NULL; /* not reached */
+  U = (GEN)hnfall(concatsp(P, subgroup))[2];
+  setlg(U, l); for (j=1; j<l; j++) setlg(U[j], l);
+  return hnfmod(concatsp(U, diagonal(raycycz)), (GEN)raycycz[1]);
 }
 
+/* t(b1,...,b_{k-1}) [prop 5.3.9] */
 static GEN
-fix_be(GEN be, GEN u)
+compute_t(GEN b, GEN r, long m, long ell)
 {
-  GEN e,g, nf = checknf(bnfz), fu = gmael(bnfz,8,5);
-  long i,lu;
+  long z, s, a, i, k = lg(b);
+  GEN t = cgetg(m+1, t_COL);
 
-  lu = lg(u);
-  for (i=1; i<lu; i++)
+  for (a = 0; a < m; a++)
   {
-    if (!signe(u[i])) continue;
-    e = mulsi(ell, (GEN)u[i]);
-    g = gmodulcp((GEN)fu[i],(GEN)nf[1]);
-    be = gmul(be, powgi(g, e));
+    z = m-1 - a; s = r[1 + z] - ell;
+    for (i=1; i < k; i++) s += r[1 + ((z + b[i]) % m)];
+    t[1 + a] = lstoi(s / ell);
   }
-  return be;
+  return t;
 }
 
+/* Return multinomial(k-1; m1,...,m_{k-1}) where the mi are the
+ * multiplicities of the bi [ assume b1 <= ... <= b_{k-1} ] */
 static GEN
-logarch2arch(GEN x, long r1, long prec)
+get_multinomial(GEN b)
 {
-  long i, lx = lg(x), tx = typ(x);
-  GEN y = cgetg(lx, tx);
+  long i, k = lg(b), a, m;
+  GEN z = mpfact(k - 1);
 
-  if (tx == t_MAT)
+  a = b[1]; m = 1;
+  for (i = 2; i < k; i++)
   {
-    for (i=1; i<lx; i++) y[i] = (long)logarch2arch((GEN)x[i], r1, prec);
-    return y;
+    if (b[i] == a) m++;
+    else
+    {
+      if (m > 1) z = diviiexact(z, mpfact(m));
+      a = b[i]; m = 1; 
+    }
   }
-  for (i=1; i<=r1;i++) y[i] = lexp((GEN)x[i],prec);
-  for (   ; i<lx; i++) y[i] = lexp(gmul2n((GEN)x[i],-1),prec);
-  return y;
+  if (m > 1) z = diviiexact(z, mpfact(m));
+  return z;
 }
 
-/* multiply be by ell-th powers of units as to find small L2-norm for new be */
+/* r[b[1]] + ... + r[b[k-1]] + 1 = 0 mod ell ?*/
+static int
+b_suitable(GEN b, GEN r, long k, long ell)
+{
+  long i, s = 1;
+  for (i = 1; i < k; i++) s += r[ 1 + b[i] ];
+  return (s % ell) == 0;
+}
+
 static GEN
-reducebetanaive(GEN be, GEN b)
+pol_from_Newton(GEN S)
 {
-  long i,k,n,ru,r1, prec = nfgetprec(bnfz);
-  GEN z,p1,p2,nmax,c, nf = checknf(bnfz);
+  long i, k, l = lg(S);
+  GEN c = cgetg(l, t_VEC);
 
-  if (DEBUGLEVEL) fprintferr("reduce modulo (Z_K^*)^l\n");
-  r1 = nf_get_r1(nf);
-  if (!b) b = gmul(gmael(nf,5,1), algtobasis(nf,be));
-  n = max((ell>>1), 3);
-  z = cgetg(n+1, t_VEC);
-  c = gmulgs(greal((GEN)bnfz[3]), ell);
-  c = logarch2arch(c, r1, prec); /* = embeddings of fu^ell */
-  c = gprec_w(gnorm(c), DEFAULTPREC);
-  b = gprec_w(gnorm(b), DEFAULTPREC); /* need little precision */
-  z[1] = (long)concatsp(c, vecinv(c));
-  for (k=2; k<=n; k++) z[k] = (long) vecmul((GEN)z[1], (GEN)z[k-1]);
-  nmax = T2_from_embed_norm(b, r1);
-  ru = lg(c)-1; c = zerovec(ru);
-  for(;;)
+  c[1] = S[1];
+  for (k = 2; k < l; k++)
   {
-    GEN B = NULL;
-    long besti = 0, bestk = 0;
-    for (k=1; k<=n; k++)
-      for (i=1; i<=ru; i++)
-      {
-        p1 = vecmul(b, gmael(z,k,i));    p2 = T2_from_embed_norm(p1,r1);
-        if (gcmp(p2,nmax) < 0) { B=p1; nmax=p2; besti=i; bestk = k; continue; }
-        p1 = vecmul(b, gmael(z,k,i+ru)); p2 = T2_from_embed_norm(p1,r1);
-        if (gcmp(p2,nmax) < 0) { B=p1; nmax=p2; besti=i; bestk =-k; }
-      }
-    if (!B) break;
-    b = B; c[besti] = laddis((GEN)c[besti], bestk);
+    GEN s = (GEN)S[k];
+    for (i = 1; i < k; i++) s = gadd(s, gmul((GEN)S[i], (GEN)c[k-i]));
+    c[k] = ldivgs(s, -k);
   }
-  if (DEBUGLEVEL) fprintferr("unit exponents = %Z\n",c);
-  return fix_be(be,c);
+  return gadd(gpowgs(polx[0], l-1), gtopoly(c, 0));
 }
 
+/* th. 5.3.5. and prop. 5.3.9. */
 static GEN
-reducebeta(GEN be)
+compute_polrel(GEN nfz, toK_s *T, GEN be, long g, long ell)
 {
-  long j,ru, prec = nfgetprec(bnfz);
-  GEN emb,z,u,matunit, nf = checknf(bnfz);
-
-  if (gcmp1(gnorm(be))) return reducebetanaive(be,NULL);
-  matunit = gmulgs(greal((GEN)bnfz[3]), ell); /* log. embeddings of fu^ell */
-  for (;;)
+  long i, k, m = T->m;
+  GEN r, powtaubet, S, X = polx[0], e = normtoK(T,be);
+    
+  switch (ell)
+  { /* special-cased for efficiency */
+    GEN p1, u;
+    case 2: err(bugparier,"rnfkummer (-1 not in nf?)"); break;
+    case 3: u = tracetoK(T,be);
+      p1 = gsub(gsqr(X), gmulsg(3,e));
+      return gsub(gmul(X,p1), gmul(e,u));
+    case 5:
+      if (m == 2)
+      {
+	u = tracetoK(T, gpowgs(be,3));
+	p1 = gadd(gmulsg(5,gsqr(e)), gmul(gsqr(X), gsub(gsqr(X),gmulsg(5,e))));
+	return gsub(gmul(X,p1), gmul(e,u));
+      }
+      else
+      { /* m = 4 */
+        GEN be1, be2, u1, u2, u3;
+        be1 = tauofelt(be, T->tau);
+        be2 = tauofelt(be1,T->tau);
+        u1 = tracetoK(T, gmul(be,be1));
+        u2 = tracetoK(T, gmul(gmul(be,be2),gsqr(be1)));
+        u3 = tracetoK(T, gmul(gmul(gsqr(be),gpowgs(be1,3)),be2));
+        p1 = gsub(gsqr(X), gmulsg(10,e));
+        p1 = gsub(gmul(X,p1), gmulsg(5,gmul(e,u1)));
+        p1 = gadd(gmul(X,p1), gmul(gmulsg(5,e),gsub(e,u2)));
+        return gsub(gmul(X,p1), gmul(e,u3));
+      }
+  }
+  /* general case */
+  r = cgetg(m+1,t_VECSMALL); /* r[i+1] = g^i mod ell */
+  r[1] = 1;
+  for (i=2; i<=m; i++) r[i] = (r[i-1] * g) % ell;
+  powtaubet = powtau(be, m, T->tau);
+  S = cgetg(ell+1, t_VEC); /* Newton sums */
+  S[1] = zero;
+  for (k = 2; k <= ell; k++)
   {
-    z = get_arch_real(nf, be, &emb, prec);
-    if (z) break;
-    prec = (prec-1)<<1;
-    if (DEBUGLEVEL) err(warnprec,"reducebeta",prec);
-    nf = nfnewprec(nf,prec);
+    GEN z, g = gzero, b = vecsmall_const(k-1, 0);
+    do
+    {
+      if (! b_suitable(b, r, k, ell)) continue;
+      z = factorbackelt(powtaubet, compute_t(b, r, m, ell), nfz);
+      if (typ(z) == t_COL) z = basistoalg(nfz, z);
+      g = gadd(g, gmul(get_multinomial(b), z));
+    } while (increment_inc(b, k, m));
+    S[k] = lmul(gmulsg(ell, e), tracetoK(T,g));
   }
-  z = concatsp(matunit, z);
-  u = lllintern(z, 1, prec);
-  if (!u) return reducebetanaive(be,emb); /* shouldn't occur */
-  ru = lg(u);
-  for (j=1; j<ru; j++)
-    if (smodis(gcoeff(u,ru-1,j), ell)) break; /* prime to ell */
-  u = (GEN)u[j]; /* coords on (fu^ell, be) of a small generator */
-  ru--; setlg(u, ru);
-  be = powgi(be, (GEN)u[ru]);
-  return reducebetanaive(fix_be(be, u), NULL);
+  return pol_from_Newton(S);
 }
 
-/* cf. algo 5.3.18 */
+typedef struct {
+  GEN R; /* compositum(P,Q) */
+  GEN p; /* Mod(p,R) root of P */
+  GEN q; /* Mod(q,R) root of Q */
+  GEN k; /* Q[X]/R generated by q + k p */
+  GEN rev;
+} compo_s;
+
 static GEN
-testx(GEN bnf, GEN X, GEN module, GEN subgroup, GEN vecMsup)
+lifttoKz(GEN nfz, GEN nf, GEN id, compo_s *C)
 {
-  long i,v,l,lX;
-  GEN be,polrelbe,p1;
-
-  if (gcmp0(X)) return NULL;
-  lX = lg(X);
-  for (i=dv+1; i<lX; i++)
-    if (gcmp0((GEN)X[i])) return NULL;
-  l = lg(vecMsup);
-  for (i=1; i<l; i++)
-    if (gcmp0(FpV_red(gmul((GEN)vecMsup[i],X), gell))) return NULL;
-  be = gun;
-  for (i=1; i<lX; i++)
-    be = gmul(be, powgi((GEN)vecw[i], (GEN)X[i]));
-  if (DEBUGLEVEL>1) fprintferr("reducing beta = %Z\n",be);
-  be = reducebeta(be);
-  if (DEBUGLEVEL>1) fprintferr("beta reduced = %Z\n",be);
-  polrelbe = computepolrelbeta(be);
-  v = varn(polrelbe);
-  p1 = unifpol((GEN)bnf[7],polrelbe,0);
-  p1 = denom(gtovec(p1));
-  polrelbe = gsubst(polrelbe,v, gdiv(polx[v],p1));
-  polrelbe = gmul(polrelbe, gpowgs(p1, degpol(polrelbe)));
-  if (DEBUGLEVEL>1) fprintferr("polrelbe = %Z\n",polrelbe);
-  p1 = rnfconductor(bnf,polrelbe,0);
-  if (!gegal((GEN)p1[1],module) || !gegal((GEN)p1[3],subgroup)) return NULL;
-  return polrelbe;
+  GEN I = ideal_two_elt(nf,id);
+  GEN x = gmul((GEN)nf[7], (GEN)I[2]);
+  I[2] = (long)algtobasis(nfz, RX_RXQ_compo(x, C->p, C->R));
+  return prime_to_ideal(nfz,I);
 }
+  
+static void
+compositum_red(compo_s *C, GEN P, GEN Q)
+{
+  GEN a, z = (GEN)compositum2(P, Q)[1];
+  C->R = (GEN)z[1];
+  C->p = (GEN)z[2];
+  C->q = (GEN)z[3];
+  C->k = (GEN)z[4];
+  /* reduce R */
+  z = polredabs0(C->R, nf_ORIG|nf_PARTIALFACT);
+  C->R = (GEN)z[1];
+  if (DEBUGLEVEL>1) fprintferr("polred(compositum) = %Z\n",C->R);
+  a    = (GEN)z[2];
+  C->p = poleval(lift_intern(C->p), a);
+  C->q = poleval(lift_intern(C->q), a);
+  C->rev = modreverse_i((GEN)a[2], (GEN)a[1]);
+}
 
-GEN
-rnfkummer(GEN bnr, GEN subgroup, long all, long prec)
+static GEN
+_rnfkummer(GEN bnr, GEN subgroup, long all, long prec)
 {
-  long i,j,l,av=avma,e,vp,vd,dK,lSl,lSp,lSl2,lSml2,dc,ru,rv,nbcol;
-  GEN p1,p2,p3,p4,wk,pr;
-  GEN bnf,rayclgp,bid,ideal,cycgen,vselmer;
-  GEN kk,clgp,fununits,torsunit,vecB,vecC,Tc,Tv,P;
-  GEN Q,Qt,idealz,gothf,factgothf,listpr,listex,factell,p,vecnul;
-  GEN M,al,K,Msup,X,finalresult,y,module,A1,A2,vecMsup;
-  GEN listmodsup,vecalphap,vecbetap,mginv,matP,ESml2,Sp,Sm,Sml1,Sml2,Sl;
+  long ell, i, j, m, d, dK, dc, rc, ru, rv, g, mginv, degK, degKz, vnf;
+  long l, lSp, lSml2, lSl2, lW;
+  GEN polnf,bnf,nf,bnfz,nfz,bid,ideal,cycgen,gell,p1,p2,wk,U,vselmer;
+  GEN clgp,cyc,gen;
+  GEN Q,idealz,gothf;
+  GEN res,u,M,K,y,vecMsup,vecW,vecWA,vecWB,vecB,vecC,vecAp,vecBp;
+  GEN matP,Sp,listprSp,Tc,Tv,P;
+  primlist L;
+  toK_s T;
+  tau_s _tau, *tau;
+  compo_s COMPO;
 
   checkbnrgen(bnr);
-  wk = gmael4(bnr,1,8,4,1);
-  if (all) gell = stoi(all);
-  else
-  {
-    if (!gcmp0(subgroup)) gell = det(subgroup);
-    else gell = det(diagonal(gmael(bnr,5,2)));
-  }
-  if (gcmp1(gell)) { avma = av; return polx[varn(gmael3(bnr,1,7,1))]; }
-  if (!isprime(gell)) err(impl,"kummer for composite relative degree");
-  if (divise(wk,gell))
-    return gerepileupto(av,rnfkummersimple(bnr,subgroup,all));
-  if (all && gcmp0(subgroup))
-    err(talker,"kummer when zeta not in K requires a specific subgroup");
   bnf = (GEN)bnr[1];
-  nf = (GEN)bnf[7];
-  polnf = (GEN)nf[1]; vnf = varn(polnf); degK = degpol(polnf);
+  nf  = (GEN)bnf[7];
+  polnf = (GEN)nf[1]; vnf = varn(polnf);
   if (!vnf) err(talker,"main variable in kummer must not be x");
-      /* step 7 */
-  p1 = conductor(bnr,subgroup,2);
-      /* fin step 7 */
-  bnr = (GEN)p1[2]; 
+  wk = gmael3(bnf,8,4,1);
+  /* step 7 */
+  if (all) subgroup = NULL;
+  p1 = conductor(bnr, subgroup, 2);
+  bnr      = (GEN)p1[2]; 
   subgroup = (GEN)p1[3];
-  rayclgp = (GEN)bnr[5];
-  raycyc = (GEN)rayclgp[2];
-  bid = (GEN)bnr[2]; module=(GEN)bid[1];
-  ideal = (GEN)module[1];
-  if (gcmp0(subgroup)) subgroup = diagonal(raycyc);
+  gell = get_gell(bnr,subgroup,all);
+  if (gcmp1(gell)) return polx[vnf];
+  if (!isprime(gell)) err(impl,"kummer for composite relative degree");
+  if (divise(wk,gell)) return rnfkummersimple(bnr, subgroup, gell, all);
+
+  bid = (GEN)bnr[2];
+  ideal = gmael(bid,1,1);
   ell = itos(gell);
-      /* step 1 of alg 5.3.5. */
+  /* step 1 of alg 5.3.5. */
   if (DEBUGLEVEL>2) fprintferr("Step 1\n");
-  
-  p1 = (GEN)compositum2(polnf, cyclo(ell,vnf))[1];
-  R = (GEN)p1[1];
-  A1= (GEN)p1[2];
-  A2= (GEN)p1[3];
-  kk= (GEN)p1[4];
-  if (signe(leadingcoeff(R)) < 0)
-  {
-    R = gneg_i(R);
-    A1 = gmodulcp(lift(A1),R);
-    A2 = gmodulcp(lift(A2),R);
-  }
-      /* step 2 */
+  compositum_red(&COMPO, polnf, cyclo(ell,vnf));
+  /* step 2 */
   if (DEBUGLEVEL>2) fprintferr("Step 2\n");
-  degKz = degpol(R);
-  m = degKz/degK;
-  d = (ell-1)/m;
-  g = lift(gpuigs(gener(gell),d)); /* g has order m in all (Z/ell^k)^* */
-  if (gcmp1(powmodulo(g, stoi(m), stoi(ell*ell)))) g = addsi(ell,g);
-      /* step reduction of R */
-  if (degKz<=20)
-  {
-    GEN A3,A3rev;
-
-    if (DEBUGLEVEL>2) fprintferr("Step reduction\n");
-    p1 = polredabs2(R,prec);
-    if (DEBUGLEVEL>2) fprintferr("polredabs = %Z",p1[1]);
-    R = (GEN)p1[1];
-    A3= (GEN)p1[2];
-    A1 = poleval(lift(A1), A3);
-    A2 = poleval(lift(A2), A3);
-    A3rev= polymodrecip(A3);
-    U = gadd(powgi(A2,g), gmul(kk,A1));
-    U = poleval(lift(A3rev), U);
-  }
-  else U = gadd(powgi(A2,g), gmul(kk,A1));
-      /* step 3 */
-      /* one could factor disc(R) using th. 2.1.6. */
+  degK  = degpol(polnf);
+  degKz = degpol(COMPO.R);
+  m = degKz / degK;
+  d = (ell-1) / m;
+  g = powuumod(u_gener(ell), d, ell);
+  if (powuumod(g, m, ell*ell) == 1) g += ell; /* ord(g)=m in all (Z/ell^k)^* */
+  /* step 3 */
   if (DEBUGLEVEL>2) fprintferr("Step 3\n");
-  bnfz = bnfinit0(R,1,NULL,prec);
+  /* could factor disc(R) using th. 2.1.6. */
+  bnfz = bnfinit0(COMPO.R,1,NULL,prec);
+  cycgen = check_and_build_cycgen(bnfz);
   nfz = (GEN)bnfz[7];
   clgp = gmael(bnfz,8,1);
-  cyc   = (GEN)clgp[2]; rc = ellrank(cyc,ell);
-  gencyc= (GEN)clgp[3]; l = lg(cyc);
-  vecalpha=cgetg(l,t_VEC);
-  cycgen = check_and_build_cycgen(bnfz);
-  for (j=1; j<l; j++)
-    vecalpha[j] = (long)basistoalg(nfz, famat_to_nf(nfz, (GEN)cycgen[j]));
-      /* computation of the uu(j) (see remark 5.2.15.) */
-  uu = cgetg(l,t_VEC);
-  for (j=1; j<=rc; j++) uu[j] = zero;
-  for (   ; j<  l; j++) uu[j] = lmpinvmod((GEN)cyc[j], gell);
+  cyc = (GEN)clgp[2]; rc = prank(cyc,ell);
+  gen = (GEN)clgp[3]; l = lg(cyc);
+  u = get_u(cyc, rc, gell);
 
-  fununits = check_units(bnfz,"rnfkummer");
-  torsunit = gmael3(bnfz,8,4,2);
-  ru=(degKz>>1)-1;
-  rv=rc+ru+1;
-  vselmer = cgetg(rv+1,t_VEC);
-  for (j=1; j<=rc; j++) vselmer[j] = vecalpha[j];
-  for (   ; j< rv; j++) vselmer[j] = lmodulcp((GEN)fununits[j-rc],R);
-  vselmer[rv]=(long)gmodulcp((GEN)torsunit,R);
-    /* step 4 */
+  vselmer = get_Selmer(bnfz, cycgen, rc);
+  ru = (degKz>>1)-1;
+  rv = rc+ru+1;
+
+  /* compute action of tau */
+  U = gadd(gpowgs(COMPO.q, g), gmul(COMPO.k, COMPO.p));
+  U = poleval(COMPO.rev, U);
+  tau = get_tau(&_tau, nfz, U);
+
+  /* step 4 */
   if (DEBUGLEVEL>2) fprintferr("Step 4\n");
   vecB=cgetg(rc+1,t_VEC);
   Tc=cgetg(rc+1,t_MAT);
   for (j=1; j<=rc; j++)
   {
-    p1 = isprincipalell(tauofideal((GEN)gencyc[j]));
+    p1 = tauofideal(nfz, (GEN)gen[j], tau);
+    p1 = isprincipalell(bnfz, p1, cycgen,u,gell,rc);
     Tc[j]  = p1[1];
     vecB[j]= p1[2];
   }
-  p1=cgetg(m,t_VEC);
-  p1[1]=(long)idmat(rc);
-  for (j=2; j<=m-1; j++) p1[j]=lmul((GEN)p1[j-1],Tc);
-  p2=cgetg(rc+1,t_VEC);
-  for (j=1; j<=rc; j++) p2[j]=un;
-  p3=vecB;
+
+  vecC = cgetg(rc+1,t_VEC);
+  for (j=1; j<=rc; j++) vecC[j] = lgetg(1, t_MAT);
+  p1 = cgetg(m,t_VEC);
+  p1[1] = (long)idmat(rc);
+  for (j=2; j<=m-1; j++) p1[j] = lmul((GEN)p1[j-1],Tc);
+  p2 = vecB;
   for (j=1; j<=m-1; j++)
   {
-    p3 = gsubst(lift(p3),vnf,U);
-    p4 = groupproduct(grouppows(p3,(j*d)%ell),(GEN)p1[m-j]);
-    for (i=1; i<=rc; i++) p2[i] = lmul((GEN)p2[i],(GEN)p4[i]);
+    GEN T = FpM_red(gmulsg((j*d)%ell,(GEN)p1[m-j]), gell);
+    p2 = tauofvec(p2, tau);
+    for (i=1; i<=rc; i++)
+      vecC[i] = (long)famat_mul((GEN)vecC[i], famat_factorback(p2, (GEN)T[i]));
   }
-  vecC=p2;
-      /* step 5 */
+  for (i=1; i<=rc; i++) vecC[i] = (long)famat_reduce((GEN)vecC[i]);
+  /* step 5 */
   if (DEBUGLEVEL>2) fprintferr("Step 5\n");
   Tv = cgetg(rv+1,t_MAT);
   for (j=1; j<=rv; j++)
   {
-    Tv[j] = isvirtualunit(gsubst(lift((GEN)vselmer[j]),vnf,U))[1];
-    if (DEBUGLEVEL>2) fprintferr("   %ld\n",j);
+    p1 = tauofelt((GEN)vselmer[j], tau);
+    if (typ(p1) == t_MAT) /* famat */
+      p1 = factorbackelt((GEN)p1[1], FpV_red((GEN)p1[2],gell), nfz);
+    Tv[j] = (long)isvirtualunit(bnfz, p1, cycgen,cyc,gell,rc);
   }
-  P = FpM_ker(gsub(Tv, g), gell);
-  dv= lg(P)-1; vecw = cgetg(dv+1,t_VEC);
-  for (j=1; j<=dv; j++)
-  {
-    p1 = gun;
-    for (i=1; i<=rv; i++) p1 = gmul(p1, powgi((GEN)vselmer[i],gcoeff(P,i,j)));
-    vecw[j] = (long)p1;
-  }
-      /* step 6 */
+  P = FpM_ker(gsubgs(Tv, g), gell);
+  lW = lg(P); vecW = cgetg(lW,t_VEC);
+  for (j=1; j<lW; j++) vecW[j] = (long)famat_factorback(vselmer, (GEN)P[j]);
+  /* step 6 */
   if (DEBUGLEVEL>2) fprintferr("Step 6\n");
-  Q = FpM_ker(gsub(gtrans(Tc), g), gell);
-  Qt = gtrans(Q); dc = lg(Q)-1;
-      /* step 7 done above */
-      /* step 8 */
-  if (DEBUGLEVEL>2) fprintferr("Step 7 and 8\n");
-  idealz=lifttokz(ideal, A1);
-  computematexpoteta1(lift(A1), R);
-  if (!divise(idealnorm(nf,ideal),gell)) gothf = idealz;
+  Q = FpM_ker(gsubgs(gtrans(Tc), g), gell);
+  /* step 8 */
+  if (DEBUGLEVEL>2) fprintferr("Step 8\n");
+  p1 = RXQ_powers(lift_intern(COMPO.p), COMPO.R, degK-1);
+  p1 = vecpol_to_mat(p1, degKz);
+  T.invexpoteta1 = invmat(p1); /* left inverse */
+  T.polnf = polnf;
+  T.tau = tau;
+  T.m = m;
+ 
+  idealz = lifttoKz(nfz, nf, ideal, &COMPO);
+  if (smodis(idealnorm(nf,ideal), ell)) gothf = idealz;
   else
-  {
-    GEN bnrz, subgroupz;
-    bnrz = buchrayinitgen(bnfz,idealz);
-    subgroupz = invimsubgroup(bnrz,bnr,subgroup);
+  { /* ell | N(ideal) */
+    GEN bnrz = buchrayinitgen(bnfz, idealz);
+    GEN subgroupz = invimsubgroup(bnrz, bnr, subgroup, &T);
     gothf = conductor(bnrz,subgroupz,0);
   }
-      /* step 9 */
-  if (DEBUGLEVEL>2) fprintferr("Step 9\n");
-  factgothf=idealfactor(nfz,gothf);
-  listpr = (GEN)factgothf[1];
-  listex = (GEN)factgothf[2]; l = lg(listpr);
-      /* step 10 and step 11 */
-  if (DEBUGLEVEL>2) fprintferr("Step 10 and 11\n");
-  Sm  = cgetg(l,t_VEC); setlg(Sm,1);
-  Sml1= cgetg(l,t_VEC); setlg(Sml1,1);
-  Sml2= cgetg(l,t_VEC); setlg(Sml2,1);
-  Sl  = cgetg(l+degKz,t_VEC); setlg(Sl,1);
-  ESml2=cgetg(l,t_VECSMALL); setlg(ESml2,1);
-  for (i=1; i<l; i++)
-  {
-    pr = (GEN)listpr[i]; p = (GEN)pr[1]; e = itos((GEN)pr[3]);
-    vp = itos((GEN)listex[i]);
-    if (!egalii(p,gell))
-    {
-      if (vp != 1) { avma = av; return gzero; }
-      if (!isconjinprimelist(Sm,pr)) _append(Sm,pr);
-    }
-    else
-    {
-      vd = (vp-1)*(ell-1)-ell*e;
-      if (vd > 0) { avma = av; return gzero; }
-      if (vd==0)
-      {
-	if (!isconjinprimelist(Sml1,pr)) _append(Sml1, pr);
-      }
-      else
-      {
-	if (vp==1) { avma = av; return gzero; }
-        if (!isconjinprimelist(Sml2,pr))
-        {
-          _append(Sml2, pr);
-          _append(ESml2,(GEN)vp);
-        }
-      }
-    }
-  }
-  factell = primedec(nfz,gell); l = lg(factell);
-  for (i=1; i<l; i++)
-  {
-    pr = (GEN)factell[i];
-    if (!idealval(nfz,gothf,pr))
-      if (!isconjinprimelist(Sl,pr)) _append(Sl, pr);
-  }
-  lSml2=lg(Sml2)-1; lSl=lg(Sl)-1; lSl2=lSml2+lSl;
-  Sp=concatsp(Sm,Sml1); lSp=lg(Sp)-1;
-      /* step 12 */
+  /* step 9, 10, 11 */
+  if (DEBUGLEVEL>2) fprintferr("Step 9, 10 and 11\n");
+  i = build_list_Hecke(&L, nfz, NULL, gothf, gell, tau);
+  if (i) return no_sol(all,i);
+
+  lSml2 = lg(L.Sml2)-1;
+  Sp = concatsp(L.Sm, L.Sml1); lSp = lg(Sp)-1;
+  listprSp = concatsp(L.Sml2, L.Sl); lSl2 = lg(listprSp)-1;
+
+  /* step 12 */
   if (DEBUGLEVEL>2) fprintferr("Step 12\n");
-  vecbetap=cgetg(lSp+1,t_VEC);
-  vecalphap=cgetg(lSp+1,t_VEC);
-  matP=cgetg(lSp+1,t_MAT);
+  vecAp = cgetg(lSp+1, t_VEC);
+  vecBp = cgetg(lSp+1, t_VEC);
+  matP  = cgetg(lSp+1, t_MAT);
   for (j=1; j<=lSp; j++)
   {
-    p1=isprincipalell((GEN)Sp[j]);
-    matP[j]=p1[1];
-    p2=gun;
-    for (i=1; i<=rc; i++)
-      p2 = gmul(p2, powgi((GEN)vecC[i],gmael(p1,1,i)));
-    p3=gdiv((GEN)p1[2],p2); vecbetap[j]=(long)p3;
-    p2=gun;
+    GEN e, a, ap;
+    p1 = isprincipalell(bnfz, (GEN)Sp[j], cycgen,u,gell,rc);
+    e = (GEN)p1[1]; matP[j] = (long)e;
+    a = (GEN)p1[2];
+    p2 = famat_mul(famat_factorback(vecC, gneg(e)), a);
+    vecBp[j] = (long)p2;
+    ap = cgetg(1, t_MAT);
     for (i=0; i<m; i++)
     {
-      p2 = gmul(p2, powgi(p3, powmodulo(g,stoi(m-1-i),gell)));
-      if (i<m-1) p3=gsubst(lift(p3),vnf,U);
+      ap = famat_mul(ap, famat_pow(p2, utoi(powuumod(g,m-1-i,ell))));
+      if (i < m-1) p2 = tauofelt(p2, tau);
     }
-    vecalphap[j]=(long)p2;
+    vecAp[j] = (long)ap;
   }
-      /* step 13 */
+  /* step 13 */
   if (DEBUGLEVEL>2) fprintferr("Step 13\n");
-  nbcol=lSp+dv;
-  vecw=concatsp(vecw,vecbetap);
-  listmod=cgetg(lSl2+1,t_VEC);
-  listunif=cgetg(lSl2+1,t_VEC);
-  listprSp = concatsp(Sml2, Sl);
-  for (j=1; j<=lSl2; j++)
+  vecWA = concatsp(vecW, vecAp);
+  vecWB = concatsp(vecW, vecBp);
+
+  /* step 14, 15, and 17 */
+  if (DEBUGLEVEL>2) fprintferr("Step 14, 15 and 17\n");
+  mginv = (m * u_invmod(g,ell)) % ell;
+  vecMsup = cgetg(lSml2+1,t_VEC);
+  M = NULL;
+  for (i=1; i<=lSl2; i++)
   {
-    GEN pr = (GEN)listprSp[j];
+    GEN pr = (GEN)listprSp[i];
     long e = itos((GEN)pr[3]), z = ell * (e / (ell-1));
-    if (j <= lSml2) z += 1 - ESml2[j];
-    listmod[j]=(long)idealpows(nfz,pr, z);
-    listunif[j]=(long)basistoalg(nfz,gdiv((GEN)pr[5],(GEN)pr[1]));
+
+    if (i <= lSml2)
+    {
+      z += 1 - L.ESml2[i];
+      vecMsup[i] = (long)logall(nfz, vecWA,lW,mginv,ell, pr,z+1);
+    }
+    M = vconcat(M, logall(nfz, vecWA,lW,mginv,ell, pr,z));
   }
-      /* step 14 and step 15 */
-  if (DEBUGLEVEL>2) fprintferr("Step 14 and 15\n");
-  listbid=cgetg(lSl2+1,t_VEC);
-  listellrank = cgetg(lSl2+1,t_VECSMALL);
-  for (i=1; i<=lSl2; i++)
+  dc = lg(Q)-1;
+  if (dc)
   {
-    listbid[i]=(long)zidealstarinitgen(nfz,(GEN)listmod[i]);
-    listellrank[i] = ellrank(gmael3(listbid,i,2,2), ell);
+    GEN QtP = gmul(gtrans_i(Q), matP);
+    M = vconcat(M, concatsp(zeromat(dc,lW-1), QtP));
   }
-  mginv = modii(mulsi(m, mpinvmod(g,gell)), gell);
-  vecnul=cgetg(dc+1,t_COL); for (i=1; i<=dc; i++) vecnul[i]=zero;
-  M=cgetg(nbcol+1,t_MAT);
-  for (j=1; j<=dv; j++)
-  {
-    p1=cgetg(1,t_COL);
-    al=(GEN)vecw[j];
-    for (i=1; i<=lSl2; i++) p1 = concatsp(p1,ideallogaux(i,al));
-    p1=gmul(mginv,p1);
-    M[j]=(long)concatsp(p1,vecnul);
-  }
-  for (   ; j<=nbcol; j++)
-  {
-    p1=cgetg(1,t_COL);
-    al=(GEN)vecalphap[j-dv];
-    for (i=1; i<=lSl2; i++) p1 = concatsp(p1,ideallogaux(i,al));
-    M[j]=(long)concatsp(p1,gmul(Qt,(GEN)matP[j-dv]));
-  }
-      /* step 16 */
+  if (!M) M = zeromat(1, lSp + lW - 1);
+  /* step 16 */
   if (DEBUGLEVEL>2) fprintferr("Step 16\n");
   K = FpM_ker(M, gell);
-  dK= lg(K)-1; if (!dK) { avma=av; return gzero; }
-      /* step 17 */
-  if (DEBUGLEVEL>2) fprintferr("Step 17\n");
-  listmodsup=cgetg(lSml2+1,t_VEC);
-  listbidsup=cgetg(lSml2+1,t_VEC);
-  listellranksup=cgetg(lSml2+1,t_VECSMALL);
-  for (i=1; i<=lSml2; i++)
-  {
-    listmodsup[i]=(long)idealmul(nfz,(GEN)listprSp[i],(GEN)listmod[i]);
-    listbidsup[i]=(long)zidealstarinitgen(nfz,(GEN)listmodsup[i]);
-    listellranksup[i] = ellrank(gmael3(listbidsup,i,2,2), ell);
-  }
-  vecMsup=cgetg(lSml2+1,t_VEC);
-  for (i=1; i<=lSml2; i++)
-  {
-    Msup=cgetg(nbcol+1,t_MAT); vecMsup[i]=(long)Msup;
-    for (j=1; j<=dv; j++) Msup[j]=lmul(mginv,ideallogauxsup(i,(GEN)vecw[j]));
-    for (   ; j<=nbcol; j++) Msup[j]=(long)ideallogauxsup(i,(GEN)vecalphap[j-dv]);
-  }
-      /* step 18 */
+  /* step 18 & ff */
   if (DEBUGLEVEL>2) fprintferr("Step 18\n");
-  do
+  dK = lg(K)-1;
+  y = cgetg(dK+1,t_VECSMALL);
+  res = cgetg(1, t_VEC); /* in case all = 1 */
+  while (dK)
   {
-    y = cgetg(dK,t_VECSMALL);
     for (i=1; i<dK; i++) y[i] = 0;
-	/* step 19 */
-    for(;;)
-    {
-      X = (GEN)K[dK];
-      for (j=1; j<dK; j++) X = gadd(X, gmulsg(y[j],(GEN)K[j]));
-      X = FpV_red(X, gell);
-      finalresult = testx(bnf,X,module,subgroup,vecMsup);
-      if (finalresult) return gerepilecopy(av, finalresult);
-        /* step 20,21,22 */
-      i = dK;
-      do
+    y[i] = 1; /* y = [0,...,0,1,0,...,0], 1 at dK'th position */
+    do
+    { /* cf. algo 5.3.18 */
+      GEN be, P, X = FpV_red(gmul_mati_smallvec(K, y), gell);
+      if (ok_congruence(X, gell, lW, vecMsup))
       {
-        i--; if (!i) goto DECREASE;
-        if (i < dK-1) y[i+1] = 0;
-        y[i]++;
-      } while (y[i] >= ell);
-    }
-DECREASE:
-    dK--;
+        be = compute_beta(X, vecWB, gell, bnfz);
+        P = compute_polrel(nfz, &T, be, g, ell);
+        if (DEBUGLEVEL>1) fprintferr("polrel(beta) = %Z\n", P);
+        if (!all && gegal(subgroup, rnfnormgroup(bnr, P))) return P; /* DONE */
+        res = concatsp(res, P);
+      }
+    } while (increment(y, dK, ell));
+    y[dK--] = 0;
   }
-  while (dK);
-  avma = av; return gzero;
+  if (all) return res;
+  return gzero; /* FAIL */
+}
+
+GEN
+rnfkummer(GEN bnr, GEN subgroup, long all, long prec)
+{
+  gpmem_t av = avma;
+  return gerepilecopy(av, _rnfkummer(bnr, subgroup, all, prec));
 }