/* Copyright (C) 2016  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. */

#include "pari.h"
#include "paripriv.h"
/*******************************************************************/
/*                  LOGARITHMIC CLASS GROUP                        */
/*******************************************************************/
/* min(v, v(Log_p Norm_{F_\p/Q_p}(x))) */
static long
vlognorm(GEN nf, GEN T, GEN x, GEN p, long v)
{
  GEN a = nf_to_scalar_or_alg(nf, x);
  GEN N = RgXQ_norm(a, T);
  if (typ(N) != t_PADIC) N = cvtop(N, p, v);
  return minss(v, valp( Qp_log(N) ));
}
/* K number field, pr a maximal ideal, let K_pr be the attached local
 * field, K_pr = Q_p[X] / (T), T irreducible. Return \tilde{e}(K_pr/Q_p) */
static long
etilde(GEN nf, GEN pr, GEN T)
{
  GEN gp = pr_get_p(pr);
  ulong e = pr_get_e(pr);
  long v, voo, vmin, p, k;

  if (!u_pval(e, gp))
  {
    v = u_pval(pr_get_f(pr), gp);
    return itou( mului(e, powiu(gp, v)) );
  }
  nf = checknf(nf);
  p = itou(gp);
  k = e / (p-1) + 1;
  /* log Norm_{F_P/Q_p} (1 + P^k) = Tr(P^k) = p^[(k + v(Diff))/ e] Z_p */
  voo = (k + idealval(nf, nf_get_diff(nf), pr)) / e;
  vmin = vlognorm(nf, T, pr_get_gen(pr), gp, voo);
  if (k > 1)
  {
    GEN U = idealprincipalunits(nf, pr, k);
    GEN gen = abgrp_get_gen(U), cyc = abgrp_get_cyc(U);
    long i, l = lg(cyc);
    for (i = 1; i < l; i++)
    {
      if (voo - Z_lval(gel(cyc,i), p) >= vmin) break;
      vmin = vlognorm(nf, T, gel(gen,i), gp, vmin);
    }
  }
  v = u_lval(degpol(T), p) + (p == 2UL? 2 : 1) - vmin;
  (void)u_lvalrem(e, p, &e);
  return e * upowuu(p,v);
}
static long
ftilde_from_e(GEN pr, long e) { return pr_get_e(pr) * pr_get_f(pr) / e; }
static long
ftilde(GEN K, GEN pr, GEN T) { return ftilde_from_e(pr, etilde(K,pr, T)); }

static long
get_ZpX_index(GEN K, GEN pr, GEN T)
{
  GEN p, pi;
  long j, l = lg(T);
  if (l == 2) return 1;
  p = pr_get_p(pr); pi = nf_to_scalar_or_alg(K, pr_get_gen(pr));
  for (j = 1; j < l; j++)
  {
    GEN t = gel(T,j);
    if (t && gvaluation(RgXQ_norm(pi, t), p)) return j;
  }
  return 0;
}

/* Given a number field K and a prime p, return
 * S = places of K above p [primedec]
 * R = corresponding p-adic factors of K.pol (mod p^k), in the same order */
static GEN
padicfact(GEN K, GEN S, long k)
{
    GEN R, p = pr_get_p(gel(S,1));
  GEN T = gel(factorpadic(nf_get_pol(K), p, k), 1);
  long l, i;
  S = idealprimedec(K, p);
  R = cgetg_copy(S, &l);
  for (i = 1; i < l; i++)
  {
    long j = get_ZpX_index(K, gel(S,i), T);
    gel(R,i) = gel(T,j);
    gel(T,j) = NULL;
  }
  return R;
}

/* K a bnf, compute Cl'(K) = ell-Sylow of Cl(K) / (places above ell).
 * Return [D, u, R0, U0, ordS]
 * - D: cyclic factors for Cl'(K)
 * - u: generators of cyclic factors (all coprime to ell)
 * - R0: subgroup isprincipal(<S>) (divides K.cyc)
 * - U0: generators of R0 are of the form S . U0
 * - ordS[i] = order of S[i] in CL(K)  */
static GEN
CL_prime(GEN K, GEN ell, GEN Sell)
{
  GEN g, ordS, R0, U0, U, D, u, cyc = bnf_get_cyc(K);
  long i, l, lD, lS = lg(Sell);

  g = leafcopy(bnf_get_gen(K));
  l = lg(g);
  for (i = 1; i < l; i++)
  {
    GEN A = gel(g,i), a = gcoeff(A,1,1);
    long v = Z_pvalrem(a, ell, &a);
    if (v) gel(g,i) = hnfmodid(A, a); /* make coprime to ell */
  }
  R0 = cgetg(lS, t_MAT);
  ordS = cgetg(lS, t_VEC);
  for (i = 1; i < lS; i++)
  {
    gel(R0,i) = isprincipal(K, gel(Sell,i));
    gel(ordS,i) = charorder(cyc, gel(R0,i)); /* order of Sell[i] */
  }
  R0 = shallowconcat(R0, diagonal_shallow(cyc));
  /* R0 = subgroup generated by S in Cl(K) [ divides diagonal(K.cyc) ]*/
  R0 = ZM_hnfall(R0, &U0, 2); /* [S | cyc] * U0 = R0 in HNF */
  D = ZM_snfall(R0, &U,NULL);
  D = RgM_diagonal_shallow(D);
  lD = lg(D);
  u = ZM_inv(U, NULL); settyp(u, t_VEC);
  for (i = 1; i < lD; i++) gel(u,i) = idealfactorback(K,g,gel(u,i),1);
  setlg(U0, l);
  U0 = rowslice(U0,1,lS-1); /* restrict to 'S' part */
  return mkvec5(D, u, R0, U0, ordS);
}

static GEN
ell1(GEN ell) { return equaliu(ell,2)? utoipos(5): addiu(ell,1); }

/* log N_{F_P/Q_p}(x) / deg_F P */
static GEN
vtilde_i(GEN K, GEN x, GEN T, GEN deg, GEN ell, long prec)
{
  GEN L, cx;
  if (typ(x) != t_POL) x = nf_to_scalar_or_alg(K, x);
  if (typ(x) != t_POL) { cx = x; L = gen_0; }
  else
  {
    GEN N;
    x = Q_primitive_part(x,&cx);
    N = RgXQ_norm(x, T);
    L = Qp_log(cvtop(N,ell,prec));
  }
  if (cx)
  {
    Q_pvalrem(cx, ell, &cx);
    if (!isint1(cx))
      L = gadd(L, gmulsg(degpol(T), Qp_log(cvtop(cx,ell,prec))));
  }
  return gdiv(L, deg);
}
static GEN
vtilde(GEN K, GEN x, GEN T, GEN deg, GEN ell, long prec)
{
  GEN G, E, vG;
  long i, l;
  if (typ(x) != t_MAT) return vtilde_i(K,x,T,deg,ell,prec);
  G = gel(x,1); vG = cgetg_copy(G, &l);
  E = gel(x,2);
  for (i = 1; i < l; i++) gel(vG, i) = vtilde_i(K, gel(G,i),T,deg,ell,prec);
  return RgV_dotproduct(E, vG);
}

/* v[i] = deg S[i] mod p^prec */
static GEN
get_vdegS(GEN Ftilde, GEN ell, long prec)
{
  long i, l = lg(Ftilde);
  GEN v = cgetg(l, t_VEC), degell = Qp_log( cvtop(ell1(ell), ell, prec) );
  for (i = 1; i < l; i++) gel(v,i) = gmulsg(Ftilde[i], degell);
  return v;
}
/* K a bnf. Compute kernel \tilde{Cl}_K(ell); return cyclic factors.
 * Set *pM to (vtilde_S[i](US[j]))_{i,j} */
static GEN
CL_tilde(GEN K, GEN US, GEN ell, GEN T, GEN Ftilde, GEN *pM, long prec)
{
  GEN D, M, ellk, vdegS;
  long i, j, imin, vmin, k, lD, l = lg(T), lU = lg(US);

  *pM = cgetg(1, t_MAT);
  if (l == 2) return cgetg(1, t_VEC); /* p = P^e: \tilde{Cl}(l) = (1) */
  vdegS = get_vdegS(Ftilde, ell, prec);
  imin = 1; vmin = l; /* upper bound */
  for (i = 1; i < l; i++)
  {
    long v = z_pval(Ftilde[i], ell);
    if (v < vmin) { vmin = v; imin = i; }
  }
  M = cgetg(lU, t_MAT);
  for (j = 1; j < lU; j++)
  {
    GEN c = cgetg(l, t_COL), a = gel(US,j);
    for (i = 1; i < l; i++)
      gel(c,i) = vtilde(K, a, gel(T,i), gel(vdegS,i), ell, prec);
    gel(M,j) = c;
  }
  k = padicprec(M, ell); ellk = powiu(ell, k);
  *pM = M = gmod(M, ellk);
  M = rowsplice(M, imin);
  l--;
  if (l == 1) return cgetg(1, t_VEC);
  M = ZM_hnfmodid(M, ellk);
  D = matsnf0(M, 4); lD = lg(D);
  if (lD > 1 && Z_pval(gel(D,1), ell) >= k) return NULL;
  return D;
}

/* [L:K] = ell^k; return 1 if L/K is locally cyclotomic at ell, 0 otherwise */
long
rnfislocalcyclo(GEN rnf)
{
  pari_sp av = avma;
  GEN K, L, S, SK, TK, SLs, SL2, TL, ell;
  ulong ll;
  long i, j, k, lk, lSK;
  checkrnf(rnf);
  lk = rnf_get_degree(rnf);
  if (lk == 1) return 1;
  k = uisprimepower(lk, &ll);
  if (!k) pari_err_IMPL("rnfislocalcyclo for non-l-extensions");
  ell = utoi(ll);
  K = rnf_get_nf(rnf);
  L = rnf_build_nfabs(rnf, nf_get_prec(K));
  S = rnfidealprimedec(rnf, ell);
  SK  = gel(S,1);
  SLs = gel(S,2);
  SL2 = shallowconcat1(SLs);
  TK = padicfact(K, SK, 100); lSK = lg(SK);
  TL = padicfact(L, SL2, 100);
  for (i = 1; i < lSK; i++)
  {
    long eK = etilde(K, gel(SK,i), gel(TK,i));
    GEN SL = gel(SLs,i);
    long lSL = lg(SL);
    for (j = 1; j < lSL; j++)
    {
      long iS = gen_search(SL2, gel(SL,j), 0, (void*)&cmp_prime_over_p,
                &cmp_nodata);
      long eL = etilde(L, gel(SL,j), gel(TL,iS));
      if (dvdui(eL/eK, ell)) { avma = av; return 0; }
    }
  };
  avma = av; return 1;
}

#if 0
/* Return 1 if L/Q is locally cyclotomic at ell */
static int
islocalcycloQ(GEN L, GEN ell)
{
  GEN SL = idealprimedec(L,ell), TL;
  long i, lSL = lg(SL);
  TL = padicfact(L,  SL, 100);
  for (i = 1; i < lSL; i++)
  {
    long eL = etilde(L, gel(SL,i), gel(TL,i));
    if (dvdui(eL,ell)) return 0;
  }
  return 1;
}
#endif

/* true nf, pr a prid */
static long
nfislocalpower_i(GEN nf, GEN pr, GEN a, GEN n)
{
  long v, e, t;
  GEN p, G, L;
  a = nf_to_scalar_or_basis(nf,a);
  if (!signe(n)) return isint1(a);
  v = nfvalrem(nf, a, pr, &a); if (!dvdsi(v, n)) return 0;
  p = pr_get_p(pr);
  v = Z_pvalrem(n, p, &n);
  if (!equali1(n))
  {
    GEN T, modpr = zk_to_Fq_init(nf, &pr, &T, &p);
    GEN ap = nf_to_Fq(nf, a, modpr);
    if (!Fq_ispower(ap, n, T, p)) return 0;
  }
  if (!v) return 1;
  e = pr_get_e(pr);
  if (v == 1) /* optimal formula */
    t = itos( divii(mului(e,p), subiu(p,1)) ) + 1;
  else /* straight Hensel */
    t = 2 * e * v + 1;
  G = Idealstarprk(nf, pr, t, nf_INIT);
  L = ideallog(nf, a, G);
  return (ZV_equal0(L) || ZV_pval(L, p) >= v);
}
long
nfislocalpower(GEN nf, GEN pr, GEN a, GEN n)
{
  pari_sp av = avma;
  long r;
  if (typ(n) != t_INT) pari_err_TYPE("nfislocalpower",n);
  nf = checknf(nf);
  checkprid(pr);
  r = nfislocalpower_i(nf, pr, a, n);
  avma = av; return r;
}

/* v_ell(  exponent(D) ) */
static long
ellexpo(GEN D, GEN ell) { return lg(D) == 1? 0: Z_pval(gel(D,1), ell); }

static GEN
ellsylow(GEN cyc, GEN ell)
{
  long i, l;
  GEN d = cgetg_copy(cyc, &l);
  for (i = 1; i < l; i++)
  {
    GEN c = gel(cyc,i), a;
    if (!Z_pvalrem(c, ell, &a)) break;
    gel(d,i) = diviiexact(c, a);
  }
  setlg(d, i); return d;
}

static long
vnorm_x(GEN nf, GEN x, GEN ell)
{
  x = nf_to_scalar_or_alg(nf,x);
  if (typ(x) != t_POL) return 0;
  x = Q_primpart(x);
  return Q_pval(nfnorm(nf,x), ell);
}
static long
vtilde_prec_x(GEN nf, GEN x, GEN ell)
{
  long i, l, v;
  GEN G;
  if (typ(x) != t_MAT) return vnorm_x(nf,x,ell);
  G = gel(x,1); l = lg(G); v = 0;
  for (i = 1; i < l; i++) v = maxss(v, vnorm_x(nf,gel(G,i),ell));
  return v;
}
/* upper bound for \delta(vec): estimate loss of accuracy when evaluating
 * \tilde{v} on the vec[i] */
static long
vtilde_prec(GEN nf, GEN vec, GEN ell)
{
  long v0 = 0, i, l = lg(vec);
  for (i = 1; i < l; i++)
    v0 = maxss(v0, vtilde_prec_x(nf, gel(vec,i), ell));
  return 3 + v0 + z_pval(nf_get_degree(nf), ell);
}

static GEN
bnflog_i(GEN bnf, GEN ell)
{
  long prec0, prec;
  GEN nf, US, vdegS, S, T, M, CLp, CLt, Ftilde, vtG, ellk;
  GEN D, Ap, cycAp, bnfS;
  long i, j, lS, lvAp;

  checkbnf(bnf);
  nf = checknf(bnf);
  S = idealprimedec(nf, ell);
  bnfS = bnfsunit0(bnf, S, nf_GENMAT, LOWDEFAULTPREC); /* S-units */
  US = leafcopy(gel(bnfS,1));
  prec0 = maxss(30, vtilde_prec(nf, US, ell));
  US = shallowconcat(bnf_get_fu(bnf), US);
  settyp(US, t_COL);
  T = padicfact(nf, S, prec0);
  lS = lg(S); Ftilde = cgetg(lS, t_VECSMALL);
  for (j = 1; j < lS; j++) Ftilde[j] = ftilde(nf, gel(S,j), gel(T,j));
  CLp = CL_prime(bnf, ell, S);
  cycAp = gel(CLp,1);
  Ap = gel(CLp,2);
  for(;;)
  {
    CLt = CL_tilde(nf, US, ell, T, Ftilde, &vtG, prec0);
    if (CLt) break;
    prec0 <<= 1;
    T = padicfact(nf, S, prec0);
  }
  prec = ellexpo(cycAp, ell) + ellexpo(CLt,ell) + 1;
  if (prec == 1) return mkvec3(cgetg(1,t_VEC), cgetg(1,t_VEC), cgetg(1,t_VEC));

  vdegS = get_vdegS(Ftilde, ell, prec0);
  ellk = powiu(ell, prec);
  lvAp = lg(Ap);
  if (lvAp > 1)
  {
    GEN Kcyc = bnf_get_cyc(bnf);
    GEN C = zeromatcopy(lvAp-1, lS-1);
    GEN Rell = gel(CLp,3), Uell = gel(CLp,4), ordS = gel(CLp,5);
    for (i = 1; i < lvAp; i++)
    {
      GEN a, b, bi, A = gel(Ap,i), d = gel(cycAp,i);
      bi = isprincipal(bnf, A);
      a = vecmodii(ZC_Z_mul(bi,d), Kcyc);
      /* a in subgroup generated by S = Rell; hence b integral */
      b = hnf_invimage(Rell, a);
      b = vecmodii(ZM_ZC_mul(Uell, ZC_neg(b)), ordS);
      A = mkvec2(A, trivial_fact());
      A = idealpowred(nf, A, d);
      /* find a principal representative of A_i^cycA_i up to elements of S */
      a = isprincipalfact(bnf,gel(A,1),S,b,nf_GENMAT|nf_FORCE);
      if (!gequal0(gel(a,1))) pari_err_BUG("bnflog");
      a = famat_mul_shallow(gel(A,2), gel(a,2)); /* principal part */
      if (lg(a) == 1) continue;
      for (j = 1; j < lS; j++)
        gcoeff(C,i,j) = vtilde(nf, a, gel(T,j), gel(vdegS,j), ell, prec0);
    }
    C = gmod(gneg(C),ellk);
    C = shallowtrans(C);
    M = mkmat2(mkcol2(diagonal_shallow(cycAp), C), mkcol2(gen_0, vtG));
    M = shallowmatconcat(M); /* relation matrix */
  }
  else
    M = vtG;
  M = ZM_hnfmodid(M, ellk);
  D = matsnf0(M, 4);
  if (lg(D) == 1 || !dvdii(gel(D,1), ellk))
    pari_err_BUG("bnflog [missing Z_l component]");
  D = vecslice(D,2,lg(D)-1);
  return mkvec3(D, CLt, ellsylow(cycAp, ell));
}
GEN
bnflog(GEN bnf, GEN ell)
{
  pari_sp av = avma;
  return gerepilecopy(av, bnflog_i(bnf, ell));
}

GEN
bnflogef(GEN nf, GEN pr)
{
  pari_sp av = avma;
  long e, f, ef;
  GEN p;
  checkprid(pr); p = pr_get_p(pr);
  nf = checknf(nf);
  e = pr_get_e(pr);
  f = pr_get_f(pr); ef = e*f;
  if (u_pval(ef, p))
  {
    GEN T = gel(factorpadic(nf_get_pol(nf), p, 100), 1);
    long j = get_ZpX_index(nf, pr, T);
    e = etilde(nf, pr, gel(T,j));
    f = ef / e;
  }
  avma = av; return mkvec2s(e,f);
}

GEN
bnflogdegree(GEN nf, GEN A, GEN ell)
{
  pari_sp av = avma;
  GEN AZ, A0Z, NA0;
  long vAZ;

  if (typ(ell) != t_INT) pari_err_TYPE("bnflogdegree", ell);
  nf = checknf(nf);
  A = idealhnf(nf, A);
  AZ = gcoeff(A,1,1);
  vAZ = Z_pvalrem(AZ, ell, &A0Z);
  if (is_pm1(A0Z))
    NA0 = gen_1;
  else
    (void)Z_pvalrem(idealnorm(nf,A), ell, &NA0);
  if (vAZ)
  {
    GEN Aell = ZM_hnfmodid(A, powiu(ell,vAZ));
    GEN S = idealprimedec(nf, ell), T;
    long l, i, s = 0;
    T = padicfact(nf, S, 100);
    l = lg(S);
    for (i = 1; i < l; i++)
    {
      GEN P = gel(S,i);
      long v = idealval(nf, Aell, P);
      if (v) s += v * ftilde(nf, P, gel(T,i));
    }
    if (s) NA0 = gmul(NA0, gpowgs(ell1(ell), s));
  }
  return gerepileupto(av, NA0);
}