/* Copyright (C) 2000  The PARI group.

This file is part of the PARI/GP package.

PARI/GP is free software; you can redistribute it and/or modify it under the
terms of the GNU General Public License as published by the Free Software
Foundation. It is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY WHATSOEVER.

Check the License for details. You should have received a copy of it, along
with the package; see the file 'COPYING'. If not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */

/*******************************************************************/
/*                                                                 */
/*                      KUMMER EXTENSIONS                          */
/*                                                                 */
/*******************************************************************/
#include "pari.h"
#include "paripriv.h"

typedef struct {
  GEN R; /* nf.pol */
  GEN x; /* tau ( Mod(x, R) ) */
  GEN zk;/* action of tau on nf.zk (as t_MAT) */
} tau_s;

typedef struct {
  GEN polnf, invexpoteta1, powg;
  tau_s *tau;
  long m;
} toK_s;

typedef struct {
  GEN R; /* ZX, compositum(P,Q) */
  GEN p; /* QX, Mod(p,R) root of P */
  GEN q; /* QX, Mod(q,R) root of Q */
  long k; /* Q[X]/R generated by q + k p */
  GEN rev;
} compo_s;

static long
prank(GEN cyc, long ell)
{
  long i;
  for (i=1; i<lg(cyc); i++)
    if (smodis(gel(cyc,i),ell)) break;
  return i-1;
}

/* increment y, which runs through [0,d-1]^(k-1). Return 0 when done. */
static int
increment(GEN y, long k, long d)
{
  long i = k, j;
  do
  {
    if (--i == 0) return 0;
    y[i]++;
  } while (y[i] >= d);
  for (j = i+1; j < k; j++) y[j] = 0;
  return 1;
}

static int
ok_congruence(GEN X, ulong ell, long lW, GEN vecMsup)
{
  long i, l;
  if (zv_equal0(X)) return 0;
  l = lg(X);
  for (i=lW; i<l; i++)
    if (X[i] == 0) return 0;
  l = lg(vecMsup);
  for (i=1; i<l; i++)
    if (zv_equal0(Flm_Flc_mul(gel(vecMsup,i),X, ell))) return 0;
  return 1;
}

static int
ok_sign(GEN X, GEN msign, GEN arch)
{
  return zv_equal(Flm_Flc_mul(msign, X, 2), arch);
}

/* REDUCTION MOD ell-TH POWERS */

#if 0
static GEN
logarch2arch(GEN x, long r1, long prec)
{
  long i, lx;
  GEN y = cgetg_copy(x, &lx);
  if (typ(x) == t_MAT)
  {
    for (i=1; i<lx; i++) gel(y,i) = logarch2arch(gel(x,i), r1, prec);
  }
  else
  {
    for (i=1; i<=r1;i++) gel(y,i) = gexp(gel(x,i),prec);
    for (   ; i<lx; i++) gel(y,i) = gexp(gmul2n(gel(x,i),-1),prec);
  }
  return y;
}
#endif

/* make be integral by multiplying by t in (Q^*)^ell */
static GEN
reduce_mod_Qell(GEN bnfz, GEN be, GEN gell)
{
  GEN c;
  be = nf_to_scalar_or_basis(bnfz, be);
  be = Q_primitive_part(be, &c);
  if (c)
  {
    GEN d, fa = factor(c);
    gel(fa,2) = FpC_red(gel(fa,2), gell);
    d = factorback(fa);
    be = typ(be) == t_INT? mulii(be,d): ZC_Z_mul(be, d);
  }
  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 = gel(fa,1); l = lg(Pr);
  Ex = gel(fa,2); q = NULL;
  for (i=1; i<l; i++)
  {
    long ex = itos(gel(Ex,i));
    GEN e = stoi(ex / n);
    if (strict && ex % n) pari_err_SQRTN("idealsqrtn", fa);
    if (q) q = idealmulpowprime(nf, q, gel(Pr,i), e);
    else   q = idealpow(nf, gel(Pr,i), e);
  }
  return q? q: gen_1;
}

static GEN
reducebeta(GEN bnfz, GEN b, GEN ell)
{
  GEN y, cb, nf = bnf_get_nf(bnfz);

  if (DEBUGLEVEL>1) err_printf("reducing beta = %Ps\n",b);
  b = reduce_mod_Qell(nf, b, ell);
  /* reduce l-th root */
  y = idealsqrtn(nf, b, ell, 0); /* (b) = y^ell I, I integral */
  if (typ(y) == t_MAT && !is_pm1(gcoeff(y,1,1)))
  {
    GEN T = idealred(nf, mkvec2(y, gen_1)), t = gel(T,2);
    /* (t)*T[1] = y, T[1] integral and small */
    if (gcmp(idealnorm(nf,t), gen_1) > 0)
      b = nfmul(nf, b, nfpow(nf, t, negi(ell)));
  }
  if (DEBUGLEVEL>1) err_printf("beta reduced via ell-th root = %Ps\n",b);
  b = Q_primitive_part(b, &cb);
  if (cb)
  {
    y = nfroots(nf, gsub(monomial(gen_1, itou(ell), fetch_var_higher()),
                         basistoalg(nf,b)));
    delete_var();
  }
  if (cb && lg(y) != 1) b = gen_1;
  else
  { /* log. embeddings of fu^ell */
    GEN fu = bnf_get_fu_nocheck(bnfz), logfu = bnf_get_logfu(bnfz);
    GEN elllogfu = RgM_Rg_mul(real_i(logfu), ell);
    long prec = nf_get_prec(nf);
    for (;;)
    {
      GEN emb, z = get_arch_real(nf, b, &emb, prec);
      if (z)
      {
        GEN ex = RgM_Babai(elllogfu, z);
        if (ex)
        {
          y = nffactorback(nf, fu, RgC_Rg_mul(ex,ell));
          b = nfdiv(nf, b, y); break;
        }
      }
      prec = precdbl(prec);
      if (DEBUGLEVEL) pari_warn(warnprec,"reducebeta",prec);
      nf = nfnewprec_shallow(nf,prec);
    }
  }
  if (cb) b = gmul(b, cb);
  if (DEBUGLEVEL>1) err_printf("beta LLL-reduced mod U^l = %Ps\n",b);
  return b;
}

/* FIXME: remove */
static GEN
tauofalg(GEN x, tau_s *tau) {
  long tx = typ(x);
  if (tx == t_POLMOD) { x = gel(x,2); tx = typ(x); }
  if (tx == t_POL) x = RgX_RgXQ_eval(x, tau->x, tau->R);
  return mkpolmod(x, tau->R);
}

/* compute Gal(K(\zeta_l)/K) */
static void
get_tau(tau_s *tau, GEN nf, compo_s *C, ulong g)
{
  GEN U;

  /* compute action of tau: q^g + kp */
  U = RgX_add(RgXQ_powu(C->q, g, C->R), RgX_muls(C->p, C->k));
  U = RgX_RgXQ_eval(C->rev, U, C->R);

  tau->x  = U;
  tau->R  = C->R;
  tau->zk = nfgaloismatrix(nf, U);
}

static GEN tauoffamat(GEN x, tau_s *tau);

static GEN
tauofelt(GEN x, tau_s *tau)
{
  switch(typ(x))
  {
    case t_COL: return RgM_RgC_mul(tau->zk, x);
    case t_MAT: return tauoffamat(x, tau);
    default: return tauofalg(x, tau);
  }
}
static GEN
tauofvec(GEN x, tau_s *tau)
{
  long i, l;
  GEN y = cgetg_copy(x, &l);
  for (i=1; i<l; i++) gel(y,i) = tauofelt(gel(x,i), tau);
  return y;
}
/* [x, tau(x), ..., tau^(m-1)(x)] */
static GEN
powtau(GEN x, long m, tau_s *tau)
{
  GEN y = cgetg(m+1, t_VEC);
  long i;
  gel(y,1) = x;
  for (i=2; i<=m; i++) gel(y,i) = tauofelt(gel(y,i-1), tau);
  return y;
}
/* x^lambda */
static GEN
lambdaofelt(GEN x, toK_s *T)
{
  tau_s *tau = T->tau;
  long i, m = T->m;
  GEN y = trivial_fact(), powg = T->powg; /* powg[i] = g^i */
  for (i=1; i<m; i++)
  {
    y = famat_mulpows_shallow(y, x, uel(powg,m-i+1));
    x = tauofelt(x, tau);
  }
  return famat_mul_shallow(y, x);
}
static GEN
lambdaofvec(GEN x, toK_s *T)
{
  long i, l;
  GEN y = cgetg_copy(x, &l);
  for (i=1; i<l; i++) gel(y,i) = lambdaofelt(gel(x,i), T);
  return y;
}

static GEN
tauoffamat(GEN x, tau_s *tau)
{
  return mkmat2(tauofvec(gel(x,1), tau), gel(x,2));
}

static GEN
tauofideal(GEN id, tau_s *tau)
{
  return ZM_hnfmodid(RgM_mul(tau->zk, id), gcoeff(id, 1,1));
}

static int
isprimeidealconj(GEN P, GEN Q, tau_s *tau)
{
  GEN p = pr_get_p(P);
  GEN x = pr_get_gen(P);
  if (!equalii(p, pr_get_p(Q))
   || pr_get_e(P) != pr_get_e(Q)
   || pr_get_f(P) != pr_get_f(Q)) return 0;
  if (ZV_equal(x, pr_get_gen(Q))) return 1;
  for(;;)
  {
    if (ZC_prdvd(x,Q)) return 1;
    x = FpC_red(tauofelt(x, tau), p);
    if (ZC_prdvd(x,P)) return 0;
  }
}

static int
isconjinprimelist(GEN S, GEN pr, tau_s *tau)
{
  long i, l;

  if (!tau) return 0;
  l = lg(S);
  for (i=1; i<l; i++)
    if (isprimeidealconj(gel(S,i),pr,tau)) return 1;
  return 0;
}

/* assume x in basistoalg form */
static GEN
downtoK(toK_s *T, GEN x)
{
  long degKz = lg(T->invexpoteta1) - 1;
  GEN y = gmul(T->invexpoteta1, Rg_to_RgC(lift_shallow(x), degKz));
  return gmodulo(gtopolyrev(y,varn(T->polnf)), T->polnf);
}

static GEN
no_sol(long all, long i)
{
  if (!all) pari_err_BUG(stack_sprintf("kummer [bug%ld]", i));
  return cgetg(1,t_VEC);
}

static GEN
get_gell(GEN bnr, GEN subgp, long all)
{
  GEN gell;
  if (all && all != -1) return utoipos(labs(all));
  if (!subgp) return ZV_prod(bnr_get_cyc(bnr));
  gell = det(subgp);
  if (typ(gell) != t_INT) pari_err_TYPE("rnfkummer",gell);
  return gell;
}

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)
{
  GEN listpr, listex, pr, factell;
  long vp, i, l, ell = itos(gell), degKz = nf_get_degree(nfz);

  if (!fa) fa = idealfactor(nfz, gothf);
  listpr = gel(fa,1);
  listex = gel(fa,2); l = lg(listpr);
  L->Sm  = vectrunc_init(l);
  L->Sml1= vectrunc_init(l);
  L->Sml2= vectrunc_init(l);
  L->Sl  = vectrunc_init(l+degKz);
  L->ESml2=vecsmalltrunc_init(l);
  for (i=1; i<l; i++)
  {
    pr = gel(listpr,i);
    vp = itos(gel(listex,i));
    if (!equalii(pr_get_p(pr), gell))
    {
      if (vp != 1) return 1;
      if (!isconjinprimelist(L->Sm,pr,tau)) vectrunc_append(L->Sm,pr);
    }
    else
    {
      long e = pr_get_e(pr), vd = (vp-1)*(ell-1)-ell*e;
      if (vd > 0) return 4;
      if (vd==0)
      {
        if (!isconjinprimelist(L->Sml1,pr,tau)) vectrunc_append(L->Sml1, pr);
      }
      else
      {
        if (vp==1) return 2;
        if (!isconjinprimelist(L->Sml2,pr,tau))
        {
          vectrunc_append(L->Sml2, pr);
          vecsmalltrunc_append(L->ESml2, vp);
        }
      }
    }
  }
  factell = idealprimedec(nfz,gell); l = lg(factell);
  for (i=1; i<l; i++)
  {
    pr = gel(factell,i);
    if (!idealval(nfz,gothf,pr) && !isconjinprimelist(L->Sl,pr,tau))
      vectrunc_append(L->Sl, pr);
  }
  return 0; /* OK */
}

/* Return a Flm */
static GEN
logall(GEN nf, GEN vec, long lW, long mginv, long ell, GEN pr, long ex)
{
  GEN m, M, sprk = zlog_pr_init(nf, pr, ex);
  long ellrank, i, l = lg(vec);

  ellrank = prank(gel(sprk,1), ell);
  M = cgetg(l,t_MAT);
  for (i=1; i<l; i++)
  {
    m = zlog_pr(nf, gel(vec,i), sprk);
    setlg(m, ellrank+1);
    if (i < lW) m = gmulsg(mginv, m);
    gel(M,i) = ZV_to_Flv(m, ell);
  }
  return M;
}

/* compute the u_j (see remark 5.2.15.) */
static GEN
get_u(GEN cyc, long rc, ulong ell)
{
  long i, l = lg(cyc);
  GEN u = cgetg(l,t_VECSMALL);
  for (i=1; i<=rc; i++) uel(u,i) = 0;
  for (   ; i<  l; i++) uel(u,i) = Fl_inv(uel(cyc,i), ell);
  return u;
}

/* alg. 5.2.15. with remark */
static GEN
isprincipalell(GEN bnfz, GEN id, GEN cycgen, GEN u, ulong ell, long rc)
{
  long i, l = lg(cycgen);
  GEN v, b, db, y = bnfisprincipal0(bnfz, id, nf_FORCE|nf_GENMAT);

  v = ZV_to_Flv(gel(y,1), ell);
  b = gel(y,2);
  if (typ(b) == t_COL)
  {
    b = Q_remove_denom(gel(y,2), &db);
    if (db) b = famat_mulpows_shallow(b, db, -1);
  }
  for (i=rc+1; i<l; i++)
  {
    ulong e = Fl_mul( uel(v,i), uel(u,i), ell);
    b = famat_mulpows_shallow(b, gel(cycgen,i), e);
  }
  setlg(v,rc+1); return mkvec2(v, b);
}

static GEN
famat_factorback(GEN v, GEN e)
{
  long i, l = lg(e);
  GEN V = trivial_fact();
  for (i=1; i<l; i++) V = famat_mulpow_shallow(V, gel(v,i), gel(e,i));
  return V;
}

static GEN
famat_factorbacks(GEN v, GEN e)
{
  long i, l = lg(e);
  GEN V = trivial_fact();
  for (i=1; i<l; i++) V = famat_mulpows_shallow(V, gel(v,i), uel(e,i));
  return V;
}

static GEN
compute_beta(GEN X, GEN vecWB, GEN ell, GEN bnfz)
{
  GEN BE, be;
  BE = famat_reduce(famat_factorbacks(vecWB, X));
  gel(BE,2) = centermod(gel(BE,2), ell);
  be = nffactorback(bnfz, BE, NULL);
  be = reducebeta(bnfz, be, ell);
  if (DEBUGLEVEL>1) err_printf("beta reduced = %Ps\n",be);
  return be;
}

static GEN
get_Selmer(GEN bnf, GEN cycgen, long rc)
{
  GEN U = bnf_build_units(bnf), tu = gel(U,1), fu = vecslice(U, 2, lg(U)-1);
  return shallowconcat(shallowconcat(fu,mkvec(tu)), vecslice(cycgen,1,rc));
}

GEN
lift_if_rational(GEN x)
{
  long lx, i;
  GEN y;

  switch(typ(x))
  {
    default: break;

    case t_POLMOD:
      y = gel(x,2);
      if (typ(y) == t_POL)
      {
        long d = degpol(y);
        if (d > 0) return x;
        return (d < 0)? gen_0: gel(y,2);
      }
      return y;

    case t_POL: lx = lg(x);
      for (i=2; i<lx; i++) gel(x,i) = lift_if_rational(gel(x,i));
      break;
    case t_VEC: case t_COL: case t_MAT: lx = lg(x);
      for (i=1; i<lx; i++) gel(x,i) = lift_if_rational(gel(x,i));
  }
  return x;
}

/* A column vector representing a subgroup of prime index */
static GEN
grptocol(GEN H)
{
  long i, j, l = lg(H);
  GEN col = cgetg(l, t_VECSMALL);
  for (i = 1; i < l; i++)
  {
    ulong ell = itou( gcoeff(H,i,i) );
    if (ell == 1) col[i] = 0; else { col[i] = ell-1; break; }
  }
  for (j=i; ++j < l; ) col[j] = itou( gcoeff(H,i,j) );
  return col;
}

/* Reorganize kernel basis so that the tests of ok_congruence can be ok
 * for y[ncyc]=1 and y[1..ncyc]=1 */
static GEN
fix_kernel(GEN K, GEN M, GEN vecMsup, long lW, long ell)
{
  pari_sp av = avma;
  long i, j, idx, ffree, dK = lg(K)-1;
  GEN Ki, Kidx = cgetg(dK+1, t_VECSMALL);

  /* First step: Gauss elimination on vectors lW...lg(M) */
  for (idx = lg(K), i=lg(M); --i >= lW; )
  {
    for (j=dK; j > 0; j--) if (coeff(K, i, j)) break;
    if (!j)
    { /* Do our best to ensure that K[dK,i] != 0 */
      if (coeff(K, i, dK)) continue;
      for (j = idx; j < dK; j++)
        if (coeff(K, i, j) && coeff(K, Kidx[j], dK) != ell - 1)
          Flv_add_inplace(gel(K,dK), gel(K,j), ell);
    }
    if (j != --idx) swap(gel(K, j), gel(K, idx));
    Kidx[idx] = i;
    if (coeff(K,i,idx) != 1)
      Flv_Fl_div_inplace(gel(K,idx), coeff(K,i,idx), ell);
    Ki = gel(K,idx);
    if (coeff(K,i,dK) != 1)
    {
      ulong t = Fl_sub(coeff(K,i,dK), 1, ell);
      Flv_sub_inplace(gel(K,dK), Flv_Fl_mul(Ki, t, ell), ell);
    }
    for (j = dK; --j > 0; )
    {
      if (j == idx) continue;
      if (coeff(K,i,j))
        Flv_sub_inplace(gel(K,j), Flv_Fl_mul(Ki, coeff(K,i,j), ell), ell);
    }
  }
  /* ffree = first vector that is not "free" for the scalar products */
  ffree = idx;
  /* Second step: for each hyperplane equation in vecMsup, do the same
   * thing as before. */
  for (i=1; i < lg(vecMsup); i++)
  {
    GEN Msup = gel(vecMsup,i);
    ulong dotprod;
    if (lgcols(Msup) != 2) continue;
    Msup = zm_row(Msup, 1);
    for (j=ffree; --j > 0; )
    {
      dotprod = Flv_dotproduct(Msup, gel(K,j), ell);
      if (dotprod)
      {
        if (j != --ffree) swap(gel(K, j), gel(K, ffree));
        if (dotprod != 1) Flv_Fl_div_inplace(gel(K, ffree), dotprod, ell);
        break;
      }
    }
    if (!j)
    { /* Do our best to ensure that vecMsup.K[dK] != 0 */
      if (Flv_dotproduct(Msup, gel(K,dK), ell) == 0)
      {
        for (j = ffree-1; j <= dK; j++)
          if (Flv_dotproduct(Msup, gel(K,j), ell)
              && coeff(K,Kidx[j],dK) != ell-1)
            Flv_add_inplace(gel(K,dK), gel(K,j), ell);
      }
      continue;
    }
    Ki = gel(K,ffree);
    dotprod = Flv_dotproduct(Msup, gel(K,dK), ell);
    if (dotprod != 1)
    {
      ulong t = Fl_sub(dotprod,1,ell);
      Flv_sub_inplace(gel(K,dK), Flv_Fl_mul(Ki,t,ell), ell);
    }
    for (j = dK; --j > 0; )
    {
      if (j == ffree) continue;
      dotprod = Flv_dotproduct(Msup, gel(K,j), ell);
      if (dotprod) Flv_sub_inplace(gel(K,j), Flv_Fl_mul(Ki,dotprod,ell), ell);
    }
  }
  if (ell == 2)
  {
    for (i = ffree, j = ffree-1; i <= dK && j; i++, j--)
    { swap(gel(K,i), gel(K,j)); }
  }
  /* Try to ensure that y = vec_ei(n, i) gives a good candidate */
  for (i = 1; i < dK; i++) Flv_add_inplace(gel(K,i), gel(K,dK), ell);
  return gerepilecopy(av, K);
}

static GEN
Flm_init(long m, long n)
{
  GEN M = cgetg(n+1, t_MAT);
  long i; for (i = 1; i <= n; i++) gel(M,i) = cgetg(m+1, t_VECSMALL);
  return M;
}
static void
Flv_fill(GEN v, GEN y)
{
  long i, l = lg(y);
  for (i = 1; i < l; i++) v[i] = y[i];
}

static GEN
get_badbnf(GEN bnf)
{
  long i, l;
  GEN bad = gen_1, gen = bnf_get_gen(bnf);
  l = lg(gen);
  for (i = 1; i < l; i++)
  {
    GEN g = gel(gen,i);
    bad = lcmii(bad, gcoeff(g,1,1));
  }
  return bad;
}
/* Let K base field, L/K described by bnr (conductor f) + H. Return a list of
 * primes coprime to f*ell of degree 1 in K whose images in Cl_f(K) generate H:
 * thus they all split in Lz/Kz; t in Kz is such that
 * t^(1/p) generates Lz => t is an ell-th power in k(pr) for all such primes.
 * Restrict to primes not dividing
 * - the index fz of the polynomial defining Kz, or
 * - the modulus, or
 * - ell, or
 * - a generator in bnf.gen or bnfz.gen */
static GEN
get_prlist(GEN bnr, GEN H, ulong ell, GEN bnfz)
{
  pari_sp av0 = avma;
  forprime_t T;
  ulong p;
  GEN L, nf, cyc, bad, cond, condZ, Hsofar;
  L = cgetg(1, t_VEC);
  cyc = bnr_get_cyc(bnr);
  nf = bnr_get_nf(bnr);

  cond = gel(bnr_get_mod(bnr), 1);
  condZ = gcoeff(cond,1,1);
  bad = get_badbnf(bnr_get_bnf(bnr));
  if (bnfz)
  {
    GEN badz = lcmii(get_badbnf(bnfz), nf_get_index(bnf_get_nf(bnfz)));
    bad = mulii(bad,badz);
  }
  bad = lcmii(muliu(condZ, ell), bad);
  /* restrict to primes not dividing bad */

  u_forprime_init(&T, 2, ULONG_MAX);
  Hsofar = cgetg(1, t_MAT);
  while ((p = u_forprime_next(&T)))
  {
    GEN LP;
    long i, l;
    if (p == ell || !umodiu(bad, p)) continue;
    LP = idealprimedec_limit_f(nf, utoipos(p), 1);
    l = lg(LP);
    for (i = 1; i < l; i++)
    {
      pari_sp av = avma;
      GEN M, P = gel(LP,i), v = bnrisprincipal(bnr, P, 0);
      if (!hnf_invimage(H, v)) { avma = av; continue; }
      M = shallowconcat(Hsofar, v);
      M = ZM_hnfmodid(M, cyc);
      if (ZM_equal(M, Hsofar)) continue;
      L = shallowconcat(L, mkvec(P));
      Hsofar = M;
      /* the primes in L generate H */
      if (ZM_equal(M, H)) return gerepilecopy(av0, L);
    }
  }
  pari_err_BUG("rnfkummer [get_prlist]");
  return NULL;
}
/*Lprz list of prime ideals in Kz that must split completely in Lz/Kz, vecWA
 * generators for the S-units used to build the Kummer generators. Return
 * matsmall M such that \prod WA[j]^x[j] ell-th power mod pr[i] iff
 * \sum M[i,j] x[j] = 0 (mod ell) */
static GEN
subgroup_info(GEN bnfz, GEN Lprz, long ell, GEN vecWA)
{
  GEN nfz = bnf_get_nf(bnfz), M, gell = utoipos(ell), Lell = mkvec(gell);
  long i, j, l = lg(vecWA), lz = lg(Lprz);
  M = cgetg(l, t_MAT);
  for (j=1; j<l; j++) gel(M,j) = cgetg(lz, t_VECSMALL);
  for (i=1; i < lz; i++)
  {
    GEN pr = gel(Lprz,i), EX = subiu(pr_norm(pr), 1);
    GEN N, g,T,p, prM = idealhnf(nfz, pr);
    GEN modpr = zk_to_Fq_init(nfz, &pr,&T,&p);
    long v = Z_lvalrem(divis(EX,ell), ell, &N) + 1; /* Norm(pr)-1 = N * ell^v */
    GEN ellv = powuu(ell, v);
    g = gener_Fq_local(T,p, Lell);
    g = Fq_pow(g,N, T,p); /* order ell^v */
    for (j=1; j < l; j++)
    {
      GEN logc, c = gel(vecWA,j);
      if (typ(c) == t_MAT) /* famat */
        c = famat_makecoprime(nfz, gel(c,1), gel(c,2), pr, prM, EX);
      c = nf_to_Fq(nfz, c, modpr);
      c = Fq_pow(c, N, T,p);
      logc = Fq_log(c, g, ellv, T,p);
      ucoeff(M, i,j) = umodiu(logc, ell);
    }
  }
  return M;
}

/* if all!=0, give all equations of degree 'all'. Assume bnr modulus is the
 * conductor */
static GEN
rnfkummersimple(GEN bnr, GEN subgroup, GEN gell, long all)
{
  long ell, i, j, degK, dK;
  long lSml2, lSl2, lSp, rc, lW;
  long prec;
  long rk=0, ncyc=0;
  GEN mat=NULL, matgrp=NULL, xell, be1 = NULL;
  long firstpass = all<0;

  GEN bnf,nf,bid,ideal,arch,cycgen;
  GEN cyc;
  GEN Sp,listprSp,matP;
  GEN res=NULL,u,M,K,y,vecMsup,vecW,vecWB,vecBp,msign;
  primlist L;

  bnf = bnr_get_bnf(bnr); (void)bnf_build_units(bnf);
  nf  = bnf_get_nf(bnf);
  degK = nf_get_degree(nf);

  bid = bnr_get_bid(bnr);
  ideal= bid_get_ideal(bid);
  arch = bid_get_arch(bid); /* this is the conductor */
  ell = itos(gell);
  i = build_list_Hecke(&L, nf, bid_get_fact2(bid), ideal, gell, NULL);
  if (i) return no_sol(all,i);

  lSml2 = lg(L.Sml2);
  Sp = shallowconcat(L.Sm, L.Sml1); lSp = lg(Sp);
  listprSp = shallowconcat(L.Sml2, L.Sl); lSl2 = lg(listprSp);

  cycgen = bnf_build_cycgen(bnf);
  cyc = bnf_get_cyc(bnf); rc = prank(cyc, ell);

  vecW = get_Selmer(bnf, cycgen, rc);
  u = get_u(ZV_to_Flv(cyc,ell), rc, ell);

  vecBp = cgetg(lSp, t_VEC);
  matP  = cgetg(lSp, t_MAT);
  for (j = 1; j < lSp; j++)
  {
    GEN L = isprincipalell(bnf,gel(Sp,j), cycgen,u,ell,rc);
    gel( matP,j) = gel(L,1);
    gel(vecBp,j) = gel(L,2);
  }
  vecWB = shallowconcat(vecW, vecBp);

  prec = DEFAULTPREC +
      nbits2extraprec(((degK-1) * (gexpo(vecWB) + gexpo(nf_get_M(nf)))));
  if (nf_get_prec(nf) < prec) nf = nfnewprec_shallow(nf, prec);
  msign = nfsign(nf, vecWB);
  arch = ZV_to_zv(arch);

  vecMsup = cgetg(lSml2,t_VEC);
  M = NULL;
  for (i = 1; i < lSl2; i++)
  {
    GEN pr = gel(listprSp,i);
    long e = pr_get_e(pr), z = ell * (e / (ell-1));

    if (i < lSml2)
    {
      z += 1 - L.ESml2[i];
      gel(vecMsup,i) = logall(nf, vecWB, 0,0, ell, pr,z+1);
    }
    M = vconcat(M, logall(nf, vecWB, 0,0, ell, pr,z));
  }
  lW = lg(vecW);
  M = vconcat(M, shallowconcat(zero_Flm(rc,lW-1), matP));
  if (!all)
  { /* primes landing in subgroup must be totally split */
    GEN Lpr = get_prlist(bnr, subgroup, ell, NULL);
    GEN M2 = subgroup_info(bnf, Lpr, ell, vecWB);
    M = vconcat(M, M2);
  }
  K = Flm_ker(M, ell);
  dK = lg(K)-1;

  if (all < 0)
    K = fix_kernel(K, M, vecMsup, lW, ell);

  y = cgetg(dK+1,t_VECSMALL);
  if (all) res = cgetg(1,t_VEC); /* in case all = 1 */
  if (all < 0)
  {
    ncyc = dK;
    mat = Flm_init(dK, ncyc);
    if (all == -1) matgrp = Flm_init(lg(bnr_get_cyc(bnr)), ncyc+1);
    rk = 0;
  }
  xell = pol_xn(ell, 0);
  do {
    dK = lg(K)-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
      {
        pari_sp av = avma;
        GEN be, P=NULL, X;
        if (all < 0)
        {
          Flv_fill(gel(mat, rk+1), y);
          setlg(mat, rk+2);
          if (Flm_rank(mat, ell) <= rk) continue;
        }
FOUND:  X = Flm_Flc_mul(K, y, ell);
        if (ok_congruence(X, ell, lW, vecMsup) && ok_sign(X, msign, arch))
        {/* be satisfies all congruences, x^ell - be is irreducible, signature
          * and relative discriminant are correct */
          if (all < 0) rk++;
          be = compute_beta(X, vecWB, gell, bnf);
          be = nf_to_scalar_or_alg(nf, be);
          if (typ(be) == t_POL) be = mkpolmod(be, nf_get_pol(nf));
          if (all == -1)
          {
            pari_sp av2 = avma;
            GEN Kgrp, colgrp = grptocol(rnfnormgroup(bnr, gsub(xell, be)));
            if (ell != 2)
            {
              if (rk == 1) be1 = be;
              else
              { /* Compute the pesky scalar */
                GEN K2, C = cgetg(4, t_MAT);
                gel(C,1) = gel(matgrp,1);
                gel(C,2) = colgrp;
                gel(C,3) = grptocol(rnfnormgroup(bnr, gsub(xell, gmul(be1,be))));
                K2 = Flm_ker(C, ell);
                if (lg(K2) != 2) pari_err_BUG("linear algebra");
                K2 = gel(K2,1);
                if (K2[1] != K2[2])
                  Flv_Fl_mul_inplace(colgrp, Fl_div(K2[2],K2[1],ell), ell);
              }
            }
            Flv_fill(gel(matgrp,rk), colgrp);
            setlg(matgrp, rk+1);
            Kgrp = Flm_ker(matgrp, ell);
            if (lg(Kgrp) == 2)
            {
              setlg(gel(Kgrp,1), rk+1);
              y = Flm_Flc_mul(mat, gel(Kgrp,1), ell);
              all = 0; goto FOUND;
            }
            avma = av2;
          }
          else
          {
            P = gsub(xell, be);
            if (all)
              res = shallowconcat(res, gerepileupto(av, P));
            else
            {
              if (ZM_equal(rnfnormgroup(bnr,P),subgroup)) return P; /*DONE*/
              avma = av; continue;
            }
          }
          if (all < 0 && rk == ncyc) return res;
          if (firstpass) break;
        }
        else avma = av;
      } while (increment(y, dK, ell));
      y[dK--] = 0;
    }
  } while (firstpass--);
  return all? res: gen_0;
}

/* alg. 5.3.11 (return only discrete log mod ell) */
static GEN
isvirtualunit(GEN bnf, GEN v, GEN cycgen, GEN cyc, GEN gell, long rc)
{
  GEN L, b, eps, y, q, nf = bnf_get_nf(bnf);
  GEN w = idealsqrtn(nf, v, gell, 1);
  long i, l = lg(cycgen);

  L = bnfisprincipal0(bnf, w, nf_GENMAT|nf_FORCE);
  q = gel(L,1);
  if (ZV_equal0(q)) { eps = v; y = q; }
  else
  {
    y = cgetg(l,t_COL);
    for (i=1; i<l; i++) gel(y,i) = diviiexact(mulii(gell,gel(q,i)), gel(cyc,i));
    eps = famat_mulpow_shallow(famat_factorback(cycgen,y), gel(L,2), gell);
    eps = famat_mul_shallow(famat_inv(eps), v);
  }
  setlg(y, rc+1);
  b = bnfisunit(bnf,eps);
  if (lg(b) == 1) pari_err_BUG("isvirtualunit");
  return shallowconcat(lift_shallow(b), y);
}

/* J a vector of elements in nfz = relative extension of nf by polrel,
 * return the Steinitz element attached to the module generated by J */
static GEN
Stelt(GEN nf, GEN J, GEN polrel)
{
  long i, l = lg(J), vx = varn(polrel);
  GEN A = cgetg(l, t_VEC), I = cgetg(l, t_VEC);
  for (i = 1; i < l; i++)
  {
    GEN v = gel(J,i);
    if (typ(v) == t_POL) { v = RgX_rem(v, polrel); setvarn(v,vx); }
    gel(A,i) = v;
    gel(I,i) = gen_1;
  }
  A = RgV_to_RgM(A, degpol(polrel));
  return idealprod(nf, gel(nfhnf(nf, mkvec2(A,I)),2));
}

static GEN
polrelKzK(toK_s *T, GEN x)
{
  GEN P = roots_to_pol(powtau(x, T->m, T->tau), 0);
  long i, l = lg(P);
  for (i=2; i<l; i++) gel(P,i) = downtoK(T, gel(P,i));
  return P;
}

/* N: Cl_m(Kz) --> Cl_m(K), lift subgroup from bnr to bnrz using Algo 4.1.11 */
static GEN
invimsubgroup(GEN bnrz, GEN bnr, GEN subgroup, toK_s *T)
{
  long l, j;
  GEN P, cyc, gen, U, polrel, StZk;
  GEN nf = bnr_get_nf(bnr), nfz = bnr_get_nf(bnrz);
  GEN polz = nf_get_pol(nfz), zkzD = nf_get_zkprimpart(nfz);

  polrel = polrelKzK(T, pol_x(varn(polz)));
  StZk = Stelt(nf, zkzD, polrel);
  cyc = bnr_get_cyc(bnrz); l = lg(cyc);
  gen = bnr_get_gen(bnrz);
  P = cgetg(l,t_MAT);
  for (j=1; j<l; j++)
  {
    GEN g, id = idealhnf_shallow(nfz, gel(gen,j));
    g = Stelt(nf, RgV_RgM_mul(zkzD, id), polrel);
    g = idealdiv(nf, g, StZk); /* N_{Kz/K}(gen[j]) */
    gel(P,j) = isprincipalray(bnr, g);
  }
  (void)ZM_hnfall_i(shallowconcat(P, subgroup), &U, 1);
  setlg(U, l); for (j=1; j<l; j++) setlg(U[j], l);
  return ZM_hnfmodid(U, cyc);
}

static GEN
pol_from_Newton(GEN S)
{
  long i, k, l = lg(S);
  GEN C = cgetg(l+1, t_VEC), c = C + 1;
  gel(c,0) = gen_1;
  gel(c,1) = gel(S,1); /* gen_0 in our case */
  for (k = 2; k < l; k++)
  {
    GEN s = gel(S,k);
    for (i = 2; i < k-1; i++) s = gadd(s, gmul(gel(S,i), gel(c,k-i)));
    gel(c,k) = gdivgs(s, -k);
  }
  return gtopoly(C, 0);
}

/* - mu_b = sum_{0 <= i < m} floor(r_b r_{d-1-i} / ell) tau^i */
static GEN
get_mmu(long b, GEN r, long ell)
{
  long i, m = lg(r)-1;
  GEN M = cgetg(m+1, t_VEC);
  for (i = 0; i < m; i++) gel(M,i+1) = stoi((r[b + 1] * r[m - i]) / ell);
  return M;
}

/* coeffs(x, a..b) in variable v >= varn(x) */
static GEN
split_pol(GEN x, long v, long a, long b)
{
  long i, l = degpol(x);
  GEN y = x + a, z;

  if (l < b) b = l;
  if (a > b || varn(x) != v) return pol_0(v);
  l = b-a + 3;
  z = cgetg(l, t_POL); z[1] = x[1];
  for (i = 2; i < l; i++) gel(z,i) = gel(y,i);
  return normalizepol_lg(z, l);
}

/* return (den_a * z) mod (v^ell - num_a/den_a), assuming deg(z) < 2*ell
 * allow either num/den to be NULL (= 1) */
static GEN
mod_Xell_a(GEN z, long v, long ell, GEN num_a, GEN den_a)
{
  GEN z1 = split_pol(z, v, ell, degpol(z));
  GEN z0 = split_pol(z, v, 0,   ell-1); /* z = v^ell z1 + z0*/
  if (den_a) z0 = gmul(den_a, z0);
  if (num_a) z1 = gmul(num_a, z1);
  return gadd(z0, z1);
}
static GEN
to_alg(GEN nfz, GEN c, long v)
{
  GEN z, D;
  if (typ(c) != t_COL) return c;
  z = gmul(nf_get_zkprimpart(nfz), c);
  if (typ(z) == t_POL) setvarn(z, v);
  D = nf_get_zkden(nfz);
  if (!equali1(D)) z = RgX_Rg_div(z, D);
  return z;
}

/* th. 5.3.5. and prop. 5.3.9. */
static GEN
compute_polrel(GEN nfz, toK_s *T, GEN be, long g, long ell)
{
  long i, k, m = T->m, vT = fetch_var(), vz = fetch_var();
  GEN r, powtaubet, S, p1, root, num_t, den_t, nfzpol, powtau_prim_invbe;
  GEN prim_Rk, C_Rk, prim_root, C_root, prim_invbe, C_invbe;
  pari_timer ti;

  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);
  if (DEBUGLEVEL>1) { err_printf("Computing Newton sums: "); timer_start(&ti); }
  prim_invbe = Q_primitive_part(nfinv(nfz, be), &C_invbe);
  powtau_prim_invbe = powtau(prim_invbe, m, T->tau);

  root = cgetg(ell + 2, t_POL);
  root[1] = evalsigne(1) | evalvarn(0);
  for (i = 0; i < ell; i++) gel(root,2+i) = gen_0;
  for (i = 0; i < m; i++)
  { /* compute (1/be) ^ (-mu) instead of be^mu [mu << 0].
     * 1/be = C_invbe * prim_invbe */
    GEN mmu = get_mmu(i, r, ell);
    /* p1 = prim_invbe ^ -mu */
    p1 = to_alg(nfz, nffactorback(nfz, powtau_prim_invbe, mmu), vz);
    if (C_invbe) p1 = gmul(p1, powgi(C_invbe, RgV_sumpart(mmu, m)));
    /* root += zeta_ell^{r_i} T^{r_i} be^mu_i */
    gel(root, 2 + r[i+1]) = monomial(p1, r[i+1], vT);
  }
  /* Other roots are as above with z_ell --> z_ell^j.
   * Treat all contents (C_*) and principal parts (prim_*) separately */
  prim_Rk = prim_root = Q_primitive_part(root, &C_root);
  C_Rk = C_root;

  r = vecsmall_reverse(r); /* theta^ell = be^( sum tau^a r_{d-1-a} ) */
  /* Compute modulo X^ell - 1, T^ell - t, nfzpol(vz) */
  p1 = to_alg(nfz, nffactorback(nfz, powtaubet, r), vz);
  num_t = Q_remove_denom(p1, &den_t);

  nfzpol = leafcopy(nf_get_pol(nfz));
  setvarn(nfzpol, vz);
  S = cgetg(ell+1, t_VEC); /* Newton sums */
  gel(S,1) = gen_0;
  for (k = 2; k <= ell; k++)
  { /* compute the k-th Newton sum */
    pari_sp av = avma;
    GEN z, D, Rk = gmul(prim_Rk, prim_root);
    C_Rk = mul_content(C_Rk, C_root);
    Rk = mod_Xell_a(Rk, 0, ell, NULL, NULL); /* mod X^ell - 1 */
    for (i = 2; i < lg(Rk); i++)
    {
      if (typ(gel(Rk,i)) != t_POL) continue;
      z = mod_Xell_a(gel(Rk,i), vT, ell, num_t,den_t); /* mod T^ell - t */
      gel(Rk,i) = RgXQX_red(z, nfzpol); /* mod nfz.pol */
    }
    if (den_t) C_Rk = mul_content(C_Rk, ginv(den_t));
    prim_Rk = Q_primitive_part(Rk, &D);
    C_Rk = mul_content(C_Rk, D); /* root^k = prim_Rk * C_Rk */

    /* Newton sum is ell * constant coeff (in X), which has degree 0 in T */
    z = polcoef_i(prim_Rk, 0, 0);
    z = polcoef_i(z      , 0,vT);
    z = downtoK(T, gmulgs(z, ell));
    if (C_Rk) z = gmul(z, C_Rk);
    gerepileall(av, C_Rk? 3: 2, &z, &prim_Rk, &C_Rk);
    if (DEBUGLEVEL>1) { err_printf("%ld(%ld) ", k, timer_delay(&ti)); err_flush(); }
    gel(S,k) = z;
  }
  if (DEBUGLEVEL>1) err_printf("\n");
  (void)delete_var();
  (void)delete_var(); return pol_from_Newton(S);
}

/* lift elt t in nf to nfz, algebraic form */
static GEN
lifttoKz(GEN nf, GEN t, compo_s *C)
{
  GEN x = nf_to_scalar_or_alg(nf, t);
  if (typ(x) != t_POL) return x;
  return RgX_RgXQ_eval(x, C->p, C->R);
}
/* lift ideal id in nf to nfz */
static GEN
ideallifttoKz(GEN nfz, GEN nf, GEN id, compo_s *C)
{
  GEN I = idealtwoelt(nf,id);
  GEN x = nf_to_scalar_or_alg(nf, gel(I,2));
  if (typ(x) != t_POL) return gel(I,1);
  gel(I,2) = algtobasis(nfz, RgX_RgXQ_eval(x, C->p, C->R));
  return idealhnf_two(nfz,I);
}
/* lift ideal pr in nf to ONE prime in nfz (the others are conjugate under tau
 * and bring no further information on e_1 W). Assume pr coprime to
 * index of both nf and nfz, and unramified in Kz/K (minor simplification) */
static GEN
prlifttoKz(GEN nfz, GEN nf, GEN pr, compo_s *C)
{
  GEN F, p = pr_get_p(pr), t = pr_get_gen(pr), T = nf_get_pol(nfz);
  if (nf_get_degree(nf) != 1)
  { /* restrict to primes above pr */
    t = Q_primpart( lifttoKz(nf,t,C) );
    T = FpX_gcd(FpX_red(T,p), FpX_red(t,p), p);
    T = FpX_normalize(T, p);
  }
  F = FpX_factor(T, p);
  return idealprimedec_kummer(nfz,gcoeff(F,1,1), pr_get_e(pr), p);
}
static GEN
get_przlist(GEN L, GEN nfz, GEN nf, compo_s *C)
{
  long i, l;
  GEN M = cgetg_copy(L, &l);
  for (i = 1; i < l; i++) gel(M,i) = prlifttoKz(nfz, nf, gel(L,i), C);
  return M;
}

static void
compositum_red(compo_s *C, GEN P, GEN Q)
{
  GEN p, q, a, z = gel(compositum2(P, Q),1);
  a = gel(z,1);
  p = gel(gel(z,2), 2);
  q = gel(gel(z,3), 2);
  C->k = itos( gel(z,4) );
  /* reduce R. FIXME: should be polredbest(a, 1), but breaks rnfkummer bench */
  z = polredabs0(a, nf_ORIG|nf_PARTIALFACT);
  C->R = gel(z,1);
  a = gel(gel(z,2), 2);
  C->p = RgX_RgXQ_eval(p, a, C->R);
  C->q = RgX_RgXQ_eval(q, a, C->R);
  C->rev = QXQ_reverse(a, C->R);
  if (DEBUGLEVEL>1) err_printf("polred(compositum) = %Ps\n",C->R);
}

/* replace P->C^(-deg P) P(xC) for the largest integer C such that coefficients
 * remain algebraic integers. Lift *rational* coefficients */
static void
nfX_Z_normalize(GEN nf, GEN P)
{
  long i, l;
  GEN C, Cj, PZ = cgetg_copy(P, &l);
  PZ[1] = P[1];
  for (i = 2; i < l; i++) /* minor variation on RgX_to_nfX (create PZ) */
  {
    GEN z = nf_to_scalar_or_basis(nf, gel(P,i));
    if (typ(z) == t_INT)
      gel(PZ,i) = gel(P,i) = z;
    else
      gel(PZ,i) = ZV_content(z);
  }
  (void)ZX_Z_normalize(PZ, &C);

  if (C == gen_1) return;
  Cj = C;
  for (i = l-2; i > 1; i--)
  {
    if (i != l-2) Cj = mulii(Cj, C);
    gel(P,i) = gdiv(gel(P,i), Cj);
  }
}

static GEN
_rnfkummer_step4(GEN bnfz, GEN gen, GEN cycgen, GEN u, ulong ell, long rc,
                 long d, long m, long g, tau_s *tau)
{
  long i, j;
  GEN vecB, vecC, Tc, Q;
  vecB=cgetg(rc+1,t_VEC);
  Tc=cgetg(rc+1,t_MAT);
  for (j=1; j<=rc; j++)
  {
    GEN p1 = tauofideal(gel(gen,j), tau);
    p1 = isprincipalell(bnfz, p1, cycgen,u,ell,rc);
    gel(Tc,j)  = gel(p1,1);
    gel(vecB,j)= gel(p1,2);
  }

  if (!rc) vecC = cgetg(1,t_VEC);
  else
  {
    GEN p1, p2;
    vecC = const_vec(rc, trivial_fact());
    p1 = Flm_powers(Tc, m-2, ell);
    p2 = vecB;
    for (j=1; j<=m-1; j++)
    {
      GEN z = Flm_Fl_mul(gel(p1,m-j), Fl_mul(j,d,ell), ell);
      p2 = tauofvec(p2, tau);
      for (i=1; i<=rc; i++)
        gel(vecC,i) = famat_mul_shallow(gel(vecC,i),
                                        famat_factorbacks(p2, gel(z,i)));
    }
    for (i=1; i<=rc; i++) gel(vecC,i) = famat_reduce(gel(vecC,i));
  }
  Q = Flm_ker(Flm_Fl_add(Flm_transpose(Tc), Fl_neg(g, ell), ell), ell);
  return mkvec2(vecC, Q);
}

static GEN
_rnfkummer_step5(GEN bnfz, GEN vselmer, GEN cycgen, GEN gell, long rc,
                 long rv, long g, tau_s *tau)
{
  GEN Tv, P, vecW;
  long j, lW;
  ulong ell = itou(gell);
  GEN cyc = bnf_get_cyc(bnfz);
  Tv = cgetg(rv+1,t_MAT);
  for (j=1; j<=rv; j++)
  {
    GEN p1 = tauofelt(gel(vselmer,j), tau);
    if (typ(p1) == t_MAT) /* famat */
      p1 = nffactorback(bnfz, gel(p1,1), FpC_red(gel(p1,2),gell));
    gel(Tv,j) = ZV_to_Flv(isvirtualunit(bnfz, p1, cycgen,cyc,gell,rc), ell);
  }
  P = Flm_ker(Flm_Fl_add(Tv, Fl_neg(g, ell), ell), ell);
  lW = lg(P);
  vecW = cgetg(lW,t_VEC);
  for (j=1; j<lW; j++) gel(vecW,j) = famat_factorbacks(vselmer, gel(P,j));
  return vecW;
}

static GEN
_rnfkummer_step18(toK_s *T, GEN bnr, GEN subgroup, GEN bnfz, GEN M,
     GEN vecWB, GEN vecMsup, ulong g, GEN gell, long lW, long all)
{
  GEN K, y, res = NULL, mat = NULL;
  long i, dK, ncyc = 0;
  ulong ell = itou(gell);
  GEN bnf = bnr_get_bnf(bnr);
  GEN nf  = bnf_get_nf(bnf);
  GEN polnf = nf_get_pol(nf);
  GEN nfz = bnf_get_nf(bnfz);
  long firstpass = all<0;
  long rk=0;
  K = Flm_ker(M, ell);
  if (all < 0)
    K = fix_kernel(K, M, vecMsup, lW, ell);
  if (DEBUGLEVEL>2) err_printf("Step 18\n");
  dK = lg(K)-1;
  y = cgetg(dK+1,t_VECSMALL);
  if (all) res = cgetg(1, t_VEC);
  if (all < 0) { ncyc = dK; rk = 0; mat = zero_Flm(lg(M)-1, ncyc); }

  do {
    dK = lg(K)-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
      { /* cf. algo 5.3.18 */
        GEN H, be, P, X = Flm_Flc_mul(K, y, ell);
        if (ok_congruence(X, ell, lW, vecMsup))
        {
          pari_sp av = avma;
          if (all < 0)
          {
            gel(mat, rk+1) = X;
            if (Flm_rank(mat,ell) <= rk) continue;
            rk++;
          }
          be = compute_beta(X, vecWB, gell, bnfz);
          P = compute_polrel(nfz, T, be, g, ell);
          nfX_Z_normalize(nf, P);
          if (DEBUGLEVEL>1) err_printf("polrel(beta) = %Ps\n", P);
          if (!all) {
            H = rnfnormgroup(bnr, P);
            if (ZM_equal(subgroup, H)) return P; /* DONE */
            avma = av; continue;
          } else {
            GEN P0 = Q_primpart(lift_shallow(P));
            GEN g = nfgcd(P0, RgX_deriv(P0), polnf, nf_get_index(nf));
            if (degpol(g)) continue;
            H = rnfnormgroup(bnr, P);
            if (!ZM_equal(subgroup,H) && !bnrisconductor(bnr,H)) continue;
          }
          P = gerepilecopy(av, P);
          res = shallowconcat(res, P);
          if (all < 0 && rk == ncyc) return res;
          if (firstpass) break;
        }
      } while (increment(y, dK, ell));
      y[dK--] = 0;
    }
  } while (firstpass--);
  return res;
}

static GEN
_rnfkummer(GEN bnr, GEN subgroup, long all, long prec)
{
  long i, j, m, d, dc, rc, ru, rv, mginv, degK, degKz, vnf;
  long lSp, lSml2, lSl2, lW;
  ulong g, ell;
  GEN polnf,bnf,nf,bnfz,nfz,bid,ideal,cycgen,gell,p1,vselmer;
  GEN cyc, gen, step4;
  GEN Q,idealz,gothf;
  GEN res=NULL,u,M,vecMsup,vecW,vecWA,vecWB,vecC,vecAp,vecBp;
  GEN matP, Sp, listprSp;
  primlist L;
  toK_s T;
  tau_s tau;
  compo_s COMPO;
  pari_timer t;

  if (DEBUGLEVEL) timer_start(&t);
  checkbnr(bnr);
  bnf = bnr_get_bnf(bnr);
  nf  = bnf_get_nf(bnf);
  polnf = nf_get_pol(nf); vnf = varn(polnf);
  if (!vnf) pari_err_PRIORITY("rnfkummer", polnf, "=", 0);
  /* step 7 */
  p1 = bnrconductor_i(bnr, subgroup, 2);
  if (DEBUGLEVEL) timer_printf(&t, "[rnfkummer] conductor");
  bnr      = gel(p1,2);
  subgroup = gel(p1,3);
  gell = get_gell(bnr,subgroup,all);
  ell = itou(gell);
  if (ell == 1) return pol_x(0);
  if (!uisprime(ell)) pari_err_IMPL("kummer for composite relative degree");
  if (all && all != -1 && umodiu(bnr_get_no(bnr), ell))
    return cgetg(1, t_VEC);
  if (bnf_get_tuN(bnf) % ell == 0)
    return rnfkummersimple(bnr, subgroup, gell, all);

  if (all == -1) all = 0;
  bid = bnr_get_bid(bnr);
  ideal = bid_get_ideal(bid);
  /* step 1 of alg 5.3.5. */
  if (DEBUGLEVEL>2) err_printf("Step 1\n");
  compositum_red(&COMPO, polnf, polcyclo(ell,vnf));
  /* step 2 */
  if (DEBUGLEVEL>2) err_printf("Step 2\n");
  if (DEBUGLEVEL) timer_printf(&t, "[rnfkummer] compositum");
  degK  = degpol(polnf);
  degKz = degpol(COMPO.R);
  m = degKz / degK;
  d = (ell-1) / m;
  g = Fl_powu(pgener_Fl(ell), d, ell);
  if (Fl_powu(g, m, ell*ell) == 1) g += ell;
  /* ord(g) = m in all (Z/ell^k)^* */
  /* step 3 */
  if (DEBUGLEVEL>2) err_printf("Step 3\n");
  /* could factor disc(R) using th. 2.1.6. */
  bnfz = Buchall(COMPO.R, nf_FORCE, maxss(prec,BIGDEFAULTPREC));
  if (DEBUGLEVEL) timer_printf(&t, "[rnfkummer] bnfinit(Kz)");
  cycgen = bnf_build_cycgen(bnfz);
  nfz = bnf_get_nf(bnfz);
  cyc = bnf_get_cyc(bnfz); rc = prank(cyc,ell);
  gen = bnf_get_gen(bnfz);
  u = get_u(ZV_to_Flv(cyc, ell), rc, ell);

  vselmer = get_Selmer(bnfz, cycgen, rc);
  if (DEBUGLEVEL) timer_printf(&t, "[rnfkummer] Selmer group");
  ru = (degKz>>1)-1;
  rv = rc+ru+1;
  get_tau(&tau, nfz, &COMPO, g);

  /* step 4 */
  if (DEBUGLEVEL>2) err_printf("Step 4\n");
  step4 = _rnfkummer_step4(bnfz, gen, cycgen, u, ell, rc, d, m, g, &tau);
  vecC = gel(step4,1);
  Q    = gel(step4,2);
  /* step 5 */
  if (DEBUGLEVEL>2) err_printf("Step 5\n");
  vecW = _rnfkummer_step5(bnfz, vselmer, cycgen, gell, rc, rv, g, &tau);
  lW = lg(vecW);
  /* step 8 */
  if (DEBUGLEVEL>2) err_printf("Step 8\n");
  p1 = RgXQ_matrix_pow(COMPO.p, degKz, degK, COMPO.R);
  T.invexpoteta1 = RgM_inv(p1); /* left inverse */
  T.polnf = polnf;
  T.tau = &tau;
  T.m = m;
  T.powg = Fl_powers(g, m, ell);

  idealz = ideallifttoKz(nfz, nf, ideal, &COMPO);
  if (umodiu(gcoeff(ideal,1,1), ell)) gothf = idealz;
  else
  { /* ell | N(ideal) */
    GEN bnrz = Buchray(bnfz, idealz, nf_INIT|nf_GEN);
    GEN subgroupz = invimsubgroup(bnrz, bnr, subgroup, &T);
    gothf = bnrconductor_i(bnrz,subgroupz,0);
  }
  /* step 9, 10, 11 */
  if (DEBUGLEVEL>2) err_printf("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);
  Sp = shallowconcat(L.Sm, L.Sml1); lSp = lg(Sp);
  listprSp = shallowconcat(L.Sml2, L.Sl); lSl2 = lg(listprSp);

  /* step 12 */
  if (DEBUGLEVEL>2) err_printf("Step 12\n");
  vecAp = cgetg(lSp, t_VEC);
  vecBp = cgetg(lSp, t_VEC);
  matP  = cgetg(lSp, t_MAT);

  for (j = 1; j < lSp; j++)
  {
    GEN e, a;
    p1 = isprincipalell(bnfz, gel(Sp,j), cycgen,u,ell,rc);
    e = gel(p1,1); gel(matP,j) = gel(p1, 1);
    a = gel(p1,2);
    gel(vecBp,j) = famat_mul_shallow(famat_factorbacks(vecC, zv_neg(e)), a);
  }
  vecAp = lambdaofvec(vecBp, &T);
  /* step 13 */
  if (DEBUGLEVEL>2) err_printf("Step 13\n");
  vecWA = shallowconcat(vecW, vecAp);
  vecWB = shallowconcat(vecW, vecBp);

  /* step 14, 15, and 17 */
  if (DEBUGLEVEL>2) err_printf("Step 14, 15 and 17\n");
  mginv = Fl_div(m, g, ell);
  vecMsup = cgetg(lSml2,t_VEC);
  M = NULL;
  for (i = 1; i < lSl2; i++)
  {
    GEN pr = gel(listprSp,i);
    long e = pr_get_e(pr), z = ell * (e / (ell-1));

    if (i < lSml2)
    {
      z += 1 - L.ESml2[i];
      gel(vecMsup,i) = logall(nfz, vecWA,lW,mginv,ell, pr,z+1);
    }
    M = vconcat(M, logall(nfz, vecWA,lW,mginv,ell, pr,z));
  }
  dc = lg(Q)-1;
  if (dc)
  {
    GEN QtP = Flm_mul(Flm_transpose(Q), matP, ell);
    M = vconcat(M, shallowconcat(zero_Flm(dc,lW-1), QtP));
  }
  if (!M) M = zero_Flm(1, lSp-1 + lW-1);

  if (!all)
  { /* primes landing in subgroup must be totally split */
    GEN lambdaWB = shallowconcat(lambdaofvec(vecW, &T), vecAp);/*vecWB^lambda*/
    GEN Lpr = get_prlist(bnr, subgroup, ell, bnfz);
    GEN Lprz= get_przlist(Lpr, nfz, nf, &COMPO);
    GEN M2 = subgroup_info(bnfz, Lprz, ell, lambdaWB);
    M = vconcat(M, M2);
  }
  if (DEBUGLEVEL>2) err_printf("Step 16\n");
  /* step 16 && 18 & ff */
  res = _rnfkummer_step18(&T,bnr,subgroup,bnfz, M, vecWB, vecMsup, g, gell, lW, all);
  return res? res: gen_0;
}

GEN
rnfkummer(GEN bnr, GEN subgroup, long all, long prec)
{
  pari_sp av = avma;
  return gerepilecopy(av, _rnfkummer(bnr, subgroup, all, prec));
}