/* 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;
  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 */

static GEN
fix_be(GEN bnfz, GEN be, GEN u)
{
  GEN nf = bnf_get_nf(bnfz), fu = bnf_get_fu_nocheck(bnfz);
  return nfmul(nf, be, nffactorback(nf, fu, u));
}

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;
}

/* multiply be by ell-th powers of units as to find small L2-norm for new be */
static GEN
reducebetanaive(GEN bnfz, GEN be, long ell, GEN elllogfu)
{
  GEN z, nmax, b, c, nf = bnf_get_nf(bnfz);
  const long r1 = nf_get_r1(nf), n = maxss(ell>>1,3);
  long i, k, ru;

  b = gmul(nf_get_M(nf), be);
  z = cgetg(n+1, t_VEC);
  c = logarch2arch(elllogfu, r1, nf_get_prec(nf)); /* embeddings of fu^ell */
  c = gprec_w(gnorm(c), DEFAULTPREC);
  b = gprec_w(gnorm(b), DEFAULTPREC); /* need little precision */
  gel(z,1) = shallowconcat(c, vecinv(c));
  for (k=2; k<=n; k++) gel(z,k) = vecmul(gel(z,1), gel(z,k-1));
  nmax = embednorm_T2(b, r1);
  ru = lg(c)-1; c = zerovec(ru);
  for(;;)
  {
    GEN B = NULL;
    long besti = 0, bestk = 0;
    for (k=1; k<=n; k++)
    {
      GEN zk = gel(z,k);
      for (i=1; i<=ru; i++)
      {
        GEN v, t;
        v = vecmul(b, gel(zk,i));    t = embednorm_T2(v,r1);
        if (gcmp(t,nmax) < 0) { B=v; nmax=t; besti=i; bestk = k; continue; }
        v = vecmul(b, gel(zk,i+ru)); t = embednorm_T2(v,r1);
        if (gcmp(t,nmax) < 0) { B=v; nmax=t; besti=i; bestk =-k; }
      }
    }
    if (!B) break;
    b = B; gel(c,besti) = addis(gel(c,besti), bestk);
  }
  if (DEBUGLEVEL) err_printf("naive reduction mod U^l: unit exp. = %Ps\n",c);
  return fix_be(bnfz, be, ZC_z_mul(c, 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. */
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 be, GEN ell)
{
  long j,ru, prec = nf_get_prec(bnfz);
  GEN emb, z, u, elllogfu, nf = bnf_get_nf(bnfz);

  if (DEBUGLEVEL>1) err_printf("reducing beta = %Ps\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 && !is_pm1(gcoeff(z,1,1)))
  {
    z = idealred_elt(nf, z);
    be = nfdiv(nf, be, nfpow(nf, z, ell));
    /* make be integral */
    be = reduce_mod_Qell(nf, be, ell);
  }
  if (DEBUGLEVEL>1) err_printf("beta reduced via ell-th root = %Ps\n",be);

  for (;;)
  {
    z = get_arch_real(nf, be, &emb, prec);
    if (z) break;
    prec = precdbl(prec);
    if (DEBUGLEVEL) pari_warn(warnprec,"reducebeta",prec);
    nf = nfnewprec_shallow(nf,prec);
  }
  /* log. embeddings of fu^ell */
  elllogfu = RgM_Rg_mul(real_i(bnf_get_logfu(bnfz)), ell);
  z = shallowconcat(elllogfu, z);
  u = lll(z);
  if (lg(u) == lg(z))
  {
    ru = lg(u);
    for (j=1; j < ru; j++)
      if (gequal1(gcoeff(u,ru-1,j))) break;
    if (j < ru)
    {
      u = gel(u,j); /* coords on (fu^ell, be) of a small generator */
      ru--; setlg(u, ru);
      be = fix_be(bnfz, be, ZC_Z_mul(u, ell));
    }
  }
  if (DEBUGLEVEL>1) err_printf("beta LLL-reduced mod U^l = %Ps\n",be);
  if (typ(be) == t_INT) return be;
  return reducebetanaive(bnfz, be, itos(ell), elllogfu);
}

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);
}

static tau_s *
get_tau(tau_s *tau, GEN nf, compo_s *C, long g)
{
  GEN bas = nf_get_zk(nf), U, Uzk;
  long i, l = lg(bas);

  /* 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;
  Uzk = cgetg(l, t_MAT);
  for (i=1; i<l; i++)
    gel(Uzk,i) = algtobasis(nf, tauofalg(gel(bas,i), tau));
  tau->zk = Uzk; return tau;
}

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 = lg(x);
  GEN y = cgetg(l, typ(x));

  for (i=1; i<l; i++) gel(y,i) = tauofelt(gel(x,i), tau);
  return y;
}
/* [x, tau(x), ..., tau^m(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;
}

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 nfz, 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(nfz,x,Q)) return 1;
    x = FpC_red(tauofelt(x, tau), p);
    if (ZC_prdvd(nfz,x,P)) return 0;
  }
}

static int
isconjinprimelist(GEN nfz, 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(nfz, 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_RgV(lift_intern(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(nfz, 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(nfz, L->Sml1,pr,tau)) vectrunc_append(L->Sml1, pr);
      }
      else
      {
        if (vp==1) return 2;
        if (!isconjinprimelist(nfz, 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))
      if (!isconjinprimelist(nfz, 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, bid = Idealstar(nf, idealpows(nf, pr, ex), nf_INIT);
  long ellrank, i, l = lg(vec);

  ellrank = prank(bid_get_cyc(bid), ell);
  M = cgetg(l,t_MAT);
  for (i=1; i<l; i++)
  {
    m = ideallog(nf, gel(vec,i), bid);
    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, GEN gell)
{
  long i, l = lg(cyc);
  GEN u = cgetg(l,t_VEC);
  for (i=1; i<=rc; i++) gel(u,i) = gen_0;
  for (   ; i<  l; i++) gel(u,i) = Fp_inv(gel(cyc,i), gell);
  return u;
}

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

  logdisc = FpC_red(gel(y,1), gell);
  b = gel(y,2);
  for (i=rc+1; i<l; i++)
  {
    GEN e = modii(mulii(gel(logdisc,i),gel(u,i)), gell);
    if (signe(e)) b = famat_mul(b, famat_pow(gel(cycgen,i), e));
  }
  setlg(logdisc,rc+1); return mkvec2(logdisc, b);
}

static GEN
famat_factorback(GEN v, GEN e)
{
  long i, l = lg(e);
  GEN V = cgetg(1, t_MAT);
  for (i=1; i<l; i++)
    if (signe(gel(e,i))) V = famat_mul(V, famat_pow(gel(v,i), gel(e,i)));
  return V;
}

static GEN
compute_beta(GEN X, GEN vecWB, GEN ell, GEN bnfz)
{
  GEN BE, be;
  BE = famat_reduce(famat_factorback(vecWB, zv_to_ZV(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 fu = bnf_get_fu_nocheck(bnf), tu = bnf_get_tuU(bnf);
  GEN units = matalgtobasis(bnf,shallowconcat(fu,tu));
  return shallowconcat(units, 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)
      Flc_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), Flc_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), Flc_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 = row_zm(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) Flc_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), Flc_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), Flc_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];
}

/* 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_get_fu(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, gel(bid,3), ideal, gell, NULL);
  if (i) return no_sol(all,i);

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

  cycgen = check_and_build_cycgen(bnf);
  cyc = bnf_get_cyc(bnf); 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++)
  {
    GEN L = isprincipalell(bnf,gel(Sp,j), cycgen,u,gell,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+1,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), ZM_to_Flm(matP, ell)));
  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 = monomial(gen_1, 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])
                  Flc_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 (gequal0(q)) { eps = v; y = q; }
  else
  {
    b = gel(L,2);
    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_mul(famat_factorback(cycgen, y), famat_pow(b, gell));
    eps = famat_mul(famat_inv(eps), v);
  }
  setlg(y, rc+1);
  b = bnfisunit(bnf,eps);
  if (lg(b) == 1) pari_err_BUG("isvirtualunit");
  return shallowconcat(lift_intern(b), y);
}

/* J a vector of elements in nfz = relative extension of nf by polrel,
 * return the Steinitz element associated to the module generated by J */
static GEN
Stelt(GEN nf, GEN J, GEN polrel)
{
  long i, l = lg(J);
  GEN A, I;

  A = cgetg(l, t_VEC);
  I = cgetg(l, t_VEC);
  for (i = 1; i < l; i++)
  {
    GEN v = gel(J,i);
    gel(A,i) = (typ(v) != t_POL)? v: RgX_rem(v, polrel);
    gel(I,i) = gen_1;
  }
  A = RgV_to_RgM(A, degpol(polrel));
  return prodid(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,raycycz,rayclgpz,raygenz,U,polrel,StZk;
  GEN nf = checknf(bnr), nfz = checknf(bnrz);
  GEN polz = nf_get_pol(nfz), zkz = nf_get_zk(nfz);

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

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;
}
/* theta^ell = be ^ ( sum tau^a r_{d-1-a} ) */
static GEN
get_reverse(GEN r)
{
  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[m - i]);
  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 v)
{
  GEN z;
  if (typ(v) != t_COL) return v;
  z = gmul(nf_get_zk(nfz), v);
  if (typ(z) == t_POL) setvarn(z, MAXVARN);
  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();
  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));
    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;

  /* Compute modulo X^ell - 1, T^ell - t, nfzpol(MAXVARN) */
  p1 = to_alg(nfz, nffactorback(nfz, powtaubet, get_reverse(r)));
  num_t = Q_remove_denom(p1, &den_t);

  nfzpol = leafcopy(nf_get_pol(nfz));
  setvarn(nfzpol, MAXVARN);
  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 = polcoeff_i(prim_Rk, 0, 0);
    z = polcoeff_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(); return pol_from_Newton(S);
}

static GEN
lifttoKz(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);
}

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(GEN bnr, GEN subgroup, long all, long prec)
{
  long ell, i, j, m, d, dK, dc, rc, ru, rv, g, mginv, degK, degKz, vnf;
  long lSp, lSml2, lSl2, lW;
  GEN polnf,bnf,nf,bnfz,nfz,bid,ideal,cycgen,gell,p1,p2,vselmer;
  GEN cyc,gen;
  GEN Q,idealz,gothf;
  GEN res=NULL,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;
  pari_timer t;
  long rk=0, ncyc=0;
  GEN mat=NULL;
  long firstpass = all<0;

  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(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 = itos(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 = (long)Fl_powu(pgener_Fl(ell), d, ell);
  if (Fl_powu((ulong)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 = check_and_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(cyc, rc, gell);

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

  /* step 4 */
  if (DEBUGLEVEL>2) err_printf("Step 4\n");
  vecB=cgetg(rc+1,t_VEC);
  Tc=cgetg(rc+1,t_MAT);
  for (j=1; j<=rc; j++)
  {
    p1 = tauofideal(gel(gen,j), tau);
    p1 = isprincipalell(bnfz, p1, cycgen,u,gell,rc);
    Tc[j]  = p1[1];
    vecB[j]= p1[2];
  }

  vecC = cgetg(rc+1,t_VEC);
  if (rc)
  {
    for (j=1; j<=rc; j++) gel(vecC,j) = cgetg(1, t_MAT);
    p1 = cgetg(m,t_VEC);
    gel(p1,1) = matid(rc);
    for (j=2; j<=m-1; j++) gel(p1,j) = gmul(gel(p1,j-1),Tc);
    p2 = vecB;
    for (j=1; j<=m-1; j++)
    {
      GEN z = FpM_red(gmulsg((j*d)%ell,gel(p1,m-j)), gell);
      p2 = tauofvec(p2, tau);
      for (i=1; i<=rc; i++)
        gel(vecC,i) = famat_mul(gel(vecC,i), famat_factorback(p2,gel(z,i)));
    }
    for (i=1; i<=rc; i++) gel(vecC,i) = famat_reduce(gel(vecC,i));
  }
  /* step 5 */
  if (DEBUGLEVEL>2) err_printf("Step 5\n");
  Tv = cgetg(rv+1,t_MAT);
  for (j=1; j<=rv; j++)
  {
    p1 = tauofelt(gel(vselmer,j), tau);
    if (typ(p1) == t_MAT) /* famat */
      p1 = nffactorback(nfz, gel(p1,1), FpC_red(gel(p1,2),gell));
    gel(Tv,j) = isvirtualunit(bnfz, p1, cycgen,cyc,gell,rc);
  }
  P = FpM_ker(RgM_Rg_add_shallow(Tv, stoi(-g)), gell);
  lW = lg(P); vecW = cgetg(lW,t_VEC);
  for (j=1; j<lW; j++) gel(vecW,j) = famat_factorback(vselmer, gel(P,j));
  /* step 6 */
  if (DEBUGLEVEL>2) err_printf("Step 6\n");
  Q = FpM_ker(RgM_Rg_add_shallow(shallowtrans(Tc), stoi(-g)), gell);
  /* 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;

  idealz = lifttoKz(nfz, nf, ideal, &COMPO);
  if (smodis(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(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)-1;
  Sp = shallowconcat(L.Sm, L.Sml1); lSp = lg(Sp)-1;
  listprSp = shallowconcat(L.Sml2, L.Sl); lSl2 = lg(listprSp)-1;

  /* step 12 */
  if (DEBUGLEVEL>2) err_printf("Step 12\n");
  vecAp = cgetg(lSp+1, t_VEC);
  vecBp = cgetg(lSp+1, t_VEC);
  matP  = cgetg(lSp+1, t_MAT);
  for (j=1; j<=lSp; j++)
  {
    GEN e, a, ap;
    p1 = isprincipalell(bnfz, gel(Sp,j), cycgen,u,gell,rc);
    e = gel(p1,1); gel(matP,j) = e;
    a = gel(p1,2);
    p2 = famat_mul(famat_factorback(vecC, gneg(e)), a);
    gel(vecBp,j) = p2;
    ap = cgetg(1, t_MAT);
    for (i=0; i<m; i++)
    {
      ap = famat_mul(ap, famat_pow(p2, utoi(Fl_powu(g,m-1-i,ell))));
      if (i < m-1) p2 = tauofelt(p2, tau);
    }
    gel(vecAp,j) = ap;
  }
  /* 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 = (m * Fl_inv(g,ell)) % ell;
  vecMsup = cgetg(lSml2+1,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 = gmul(shallowtrans(Q), matP);
    M = vconcat(M, shallowconcat(zero_Flm(dc,lW-1), ZM_to_Flm(QtP,ell)));
  }
  if (!M) M = zero_Flm(1, lSp + lW - 1);
  /* step 16 */
  if (DEBUGLEVEL>2) err_printf("Step 16\n");
  K = Flm_ker(M, ell);
  if (all < 0)
    K = fix_kernel(K, M, vecMsup, lW, ell);
  /* step 18 & ff */
  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 (DEBUGLEVEL) timer_printf(&t, "[rnfkummer] candidate list");
  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(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--);
  if (all) return res;
  return gen_0; /* FAIL */
}

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