/*

packageaurel.gp

This is a special version of the abelianbnf package
  https://hal.inria.fr/hal-02961482
implementing algorithms from
  https://hal.inria.fr/hal-02497890
distributed only for the Atelier Pari/GP 2022.

*/


/*
  Sections:
    I)      Utilities
    II)     Accessors
    III)    Linear algebra
    IV)     Number fields utilities
    V)      Brauer class group
    VI)     Galois, relations and subfields computations
    VII)    Auxiliary functions for type 1
    VIII)   Main functions
    IX)     Certification
    X)      Example generation
*/

/*     I) Utilities         */

install(ZM_snf_group,"mGD&D&");

vprintab = 1; \\verbose print on/off

\\only remove powers of small primes
cheaprad(N) =
{
  my(M=N,D=1,k,m);
  forprime(p=2,10^5, \\useless to remove more: slow for only a few bits
    v = valuation(M,p);
    if(v>0, M \= p^v; D*=p);
  );
  k = ispower(M,,&m);
  if(k,m*D,M*D)
};

expo(x) = if(x==0,-oo,exponent(x)+1);

preceval0(c,ea,i) =
{
  if(i==0,return(0));
  if(c==0,return(0));
  i*ea + expo(c) + expo(i)
};
preceval(f,a) =
{
  my(ea=expo(a),d=poldegree(f));
  if(d<=0,return(0));
  expo(d) + vecmax([preceval0(polcoef(f,i),ea,i) | i <- [0..d]])
};


/*     II) Accessors         */

getexpo(abbnf) =
{
  my(cyc);
  cyc = getcyc(abbnf);
  if(cyc==[], 1, cyc[1]);
};

getcyc(abbnf) =
{
  if(abbnf[1]==2, abbnf[2]
  ,abbnf[1]==1, abbnf[2]
  ,/*else: 0*/ abbnf[2].cyc
  );
};

getpol(abbnf) =
{
  if(abbnf[1]==2, abbnf[4]
  ,abbnf[1]==1, abbnf[3]
  ,/*else: 0*/ abbnf[2].pol
  );
};

getdisc(abbnf) =
{
  if(abbnf[1]==2, error("not implemented");
  ,abbnf[1]==1, abbnf[18]
  ,/*else: 0*/ abbnf[2].disc
  );
};

abgal_gal(abgal) = abgal[1];
abgal_galgen(abgal) = abgal[2];
abgal_U(abgal) = abgal[3];
abgal_Ui(abgal) = abgal[4];
abgal_cyc(abgal) = abgal[5];


/*     III) Linear algebra         */

\\assume n is a prime power
kermodprimepow(M,n) =
{
  my(e,p,K,Kfin,pid);
  Kfin = matid(matsize(M)[2]);
  e = isprimepower(n,&p);
  pid = p*matid(matsize(M)[2]);
  for(i=1,e,
    K = liftint(matker(Mod(M,p)));
    K = matconcat([pid,K]);
    K = matimagemod(K,n);
    Kfin = (Kfin*K)%n;
    M = ((M*K)%n)\p;
  );
  Kfin = matimagemod(Kfin,n);
  Kfin
};

\\cardinal of subgroup of (Z/d)^m represented by Howell form H
cardmod(H,d) =
{
  my(m,n,card=1,i);
  [m,n] = matsize(H);
  if (m*n==0, return(1));
  i=m;
  forstep(j=n,1,-1,
    while(H[i,j]==0,i--);
    card *= d\H[i,j];
  );
  card
};

\\p-part / p'-part
splitcycsylow(cyc,p) =
{
  my(cycp,cycc);
  cycp = vector(#cyc,i,1);
  cycc = vector(#cyc,i,1);
  for(i=1,#cyc,
    cycp[i] = p^valuation(cyc[i],p);
    cycc[i] = cyc[i] \ cycp[i];
  );
  cycp = [d | d <- cycp, d!=1];
  cycc = [d | d <- cycc, d!=1];
  [cycp, cycc]
};

\\assumed to be coprime
mergecyc(cyc1,cyc2) =
{
  my(n);
  n = max(#cyc1,#cyc2);
  cyc1 = concat(cyc1,vector(n,i,1));
  cyc2 = concat(cyc2,vector(n,i,1));
  vector(n,i,cyc1[i]*cyc2[i])
};

cycsimplify(cyc,{d=0}) =
{
  [di | di <-apply(c->gcd(d,c),cyc), di!=1]
};


/*     IV) Number fields utilities         */

\\Minkowski bound
minkow(D,n,{r2=-1}) =
{
  if(r2<0,if(n%2,r2=0,r2=n\2));
  sqrt(abs(D))*(4/Pi)^r2*n!/n^n
};
nfminkow(nf) =
{
  floor(minkow(nf.disc,poldegree(nf.pol),nf.r2))
};

easynf(pol) =
{
  my(nf,cert,fa);
  if(vprintab,print("pol=",pol));
  fa = poldiscfactors(pol)[2];
  nf = nfinit([pol,fa[,1]]);
  cert = nfcertify(nf);
  if(#cert==0,return(nf));
  warning("easynf: computing nf the hard way!");
  /* Here, could be more clever, but rarely ever reached. */
  nfinit(pol)
};

nfhash3(pol) =
{
  my(fa);
  fa = nfdiscfactors([pol,10^5])[2];
  fa = matconcat([f | f <- Col(fa), f[1] <= 10^5]~);
  if(fa==[;],return(1));
  factorback(fa)
};
\\[degree, signature, abs(disc) at primes <= 10^5]
nfhash(pol) =
{
  my(n,s,d);
  n = poldegree(pol);
  s = polsturm(pol);
  d = nfhash3(pol);
  [n,s,d]
};

\\assume pol sqrfree mod p
\\return a vector of [p,1,alpha,f,pi] representing prime ideals:
\\  - [p,alpha] generates pr
\\  - f is the residue degree
\\  - pi has valuation 0 at pr and 1 at other prime divisors of p
myprimedec(pol,p) =
{
  my(fa,lifa,P,good);
  fa = factormod(pol,p)[,1]~;
  if(#fa==1, return([[p,p,poldegree(pol),1]]));
  lifa = liftint(fa);
  \\ensure generators have valuation 1 at corresponding primes
  while(1,
    P = (factorback(lifa)-pol)\p;
    good=1;
    for(i=1,#fa,
      if(Mod(P,fa[i])==0,good=0;break);
    );
    if(good,break);
    lifa = vector(#lifa,i,liftint(Mod(lifa[i]+p*liftint(random(fa[i])),p^2)));
  );
  vector(#fa,i,[p,1,lifa[i],poldegree(fa[i]),liftint(prod(j=1,#fa,if(j==i,1,fa[j])))]);
}

myidealnorm(polrel,idl) =
{
  my(n,k,a,na=[],d,m);
  [n,k,a] = idl;
  d = poldegree(polrel);
  a = Mod(a,polrel);
  m = n^(d*k+1);
  iferr(
    na = liftall(norm(Mod(a,m)));
  ,E,
    na = norm(a);
    na = liftall(Mod(na,m));
  );
  [n,d*k,na]
};

unitrank(nf) =
{
  nf.r1+nf.r2-1
};
colinfdeg(nf) =
{
  vectorv(nf.r1+nf.r2,i,if(i<=nf.r1,1,2));
};
preMi(bnf) =
{
  M = matconcat([real(bnf[3]),colinfdeg(bnf)]);
  M^(-1)
};
logisunit(Mi,logu) =
{
  my(X,rnd,e);
  X = Mi*logu;
  rnd = round(X,&e);
  if(e>-5, error("logisunit: too low precision! e=",e));
  rnd[1..-2]
};

mykeycx(z) = [abs(imag(z)),real(z),imag(z)];
myinfplaces(pol) =
{
  my(r1,r2,evalr,evalc,mycmpcx);
  r1 = polsturm(pol);
  r2 = (poldegree(pol)-r1)\2;
  if(vprintab,print("signature: ", [r1,r2], " unit rank: ", r1+r2-1));
  if(r1, evalr = polrootsreal(pol), evalr=[]~);
  if(r2,
    evalc = polroots(pol);
    evalc = vecsort(evalc,mykeycx);
    evalc = vector(r2,i,evalc[r1+2*i])~;
  ,\\else
    evalc=[]~
  );
  [r1,r2,concat(evalr,evalc)];
};

\\precision to evaluate h*R from Brauer relation
precbrauerhr(Lbnf,coefsbr) =
{
  vecmax(vector(#Lbnf,i,my(bnf=Lbnf[i][2]);
    expo(max(5,bnf.reg)*max(1,bnf.no/bnf.tu[1]))*2*abs(coefsbr[i])))
    + expo(#Lbnf) + 5;
};

\\precision to compute regulator
precdirectR0(matextru,r0,R00) =
{
  sum(i=1,r0,ceil(expo(norml2(matextru[,i]))/2))
    - exponent(R00) + 2*expo(r0) + 5;
};

whichprime(nfsub,decsub,prsub) =
{
  my(v,x=prsub[3]);
  for(i=1,#decsub,
    v = nfeltval(nfsub,x,decsub[i]);
    if(v>0,return(i));
  );
  0
};

\\compute valuations from subfields
valfromsub(bnfsub,x,isub,S,valdata=0) =
{
  my(data,vals,decsub,respr);
  concat(vector(#S,is,
    data = S[is];
    decsub = data[3][isub][1];
    respr = data[3][isub][2];
    if(valdata==0,
      vals = vector(#decsub,i,nfeltval(bnfsub,x,decsub[i]))
    ,\\else
      vals = vector(#decsub,i,valdata[1][vecsearch(valdata[2],decsub[i])])
    );
    vectorv(#respr,i,vals[respr[i]])
  ))
};

\\compute DL of residue modulo prime ideals from subfields
resfromsub(bnfsub,x,isub,T) =
{
  my(modpr,g,xmod,data,expo,n,q);
  vectorv(#T,it,
    data = T[it];
    q = data[1][1];
    g = data[2];
    [expo,n] = data[3];
    modpr = data[4][isub];
    xmod = nfmodpr(bnfsub,x,modpr);
    if(type(xmod)!="t_INT", xmod = polcoef(xmod.pol,0));
    xmod = Mod(xmod,q);
    znlog(xmod^expo,g,n);
  )
};

\\output: [S-units, matrix of their valuations, sorted S]
mysunits(bnf,S,den) =
{
  my(su,K);
  S = vecsort(S);
  su = bnfunits(bnf,S)[1][1..#S];
  K = matrix(#S,#su,i,j,nfeltval(bnf,su[j],S[i]));
  [su,K,S]
};

fetchbnf(pol,Lbnfdata) =
{
  my(n,s,d,L,isom,typ);
  n = poldegree(pol);
  L = [i | i <- [1..#Lbnfdata], Lbnfdata[i][1][1]==n];
  if(#L == 0, return(0));
  s = polsturm(pol);
  L = [i | i <- L, Lbnfdata[i][1][2]==s];
  if(#L == 0, return(0));
  d = nfhash3(pol);
  L = [i | i <- L, Lbnfdata[i][1][3]==d];
  if(#L == 0, return(0));
  foreach(L,i,
    if(getpol(Lbnfdata[i][2])==pol,return([Lbnfdata[i][2],variable(pol),variable(pol)]));
  );
  for(j=1,#L,
    my(i=L[j]);
    typ = Lbnfdata[i][2][1];
    if(typ==0,
      isom = nfisisom(Lbnfdata[i][2][2],pol)
    ,\\else
      isom = nfisisom(getpol(Lbnfdata[i][2]),pol)
    );
    if(isom, return([Lbnfdata[i][2],isom[1],lift(modreverse(Mod(isom[1],pol)))]))
  );
  0
};


/*     V) Brauer class group         */

\\class group at p
\\assume brauer is valid at p
Brauerclgp0(Lbnf,brauer,pk) =
{
  my(Lscyc,mult,cyc,m,co);
  Lscyc = apply(bnf->cycsimplify(getcyc(bnf),pk),Lbnf);
  mult = Map();
  for(i=1,#Lscyc,
    cyc = Lscyc[i];
    co = brauer[i];
    for(j=1,#cyc,
      if(!mapisdefined(mult,cyc[j],&m),
        mapput(mult,cyc[j],co)
      ,\\else
        mapput(mult,cyc[j],m+co)
      )
    );
  );
  mult = Vec(Mat(mult)~);
  mult = [m | m <- mult, m[2]!=0];
  if(#mult==0, [], vecsort(concat([vector(m[2],i,m[1]) | m <- mult, m[2]!=0]),,4))
};

\\class group at d, Lp list of primes corresponding to list of relations
Brauerclgp(Lbnf,Lp,brauer,d) =
{
  my(cyc,Lq,pk);
  Lq = factor(d/prod(i=1,#Lp,p=Lp[i];p^valuation(d,p)))[,1]~;
  cyc = [];
  for(i=1,#Lq,
    p = Lq[i];
    pk = p^valuation(d,p);
    cyc = mergecyc(cyc,Brauerclgp0(Lbnf,brauer[1],pk));
  );
  for(i=1,#Lp,
    p = Lp[i];
    pk = p^valuation(d,p);
    if(pk!=1,
      cyc = mergecyc(cyc,Brauerclgp0(Lbnf,brauer[i],pk));
    );
  );
  cyc
};


/*     VI) Galois, relations and subfields computations         */

\\find a suitable abelian subgroup of the Galois group
abgaloisanalysis(pol) =
{
  my(gal,galgen,U,Ui,cyc,rk,rkbest,rel,relbest,Lsub,sub);
  gal = galoisinit(pol);
  if(!gal, error("The field is not Galois."));
  rel = galoisisabelian(gal);
  if(rel,
    galgen = gal.gen;
  ,\\else: non-abelian
    Lsub = galoissubgroups(gal);
    cardbest=0;
    for(i=1,#Lsub,
      sub = Lsub[i];
      rel = galoisisabelian(sub);
      if(rel,
        rk = #matsnf(rel,4);
        if(rk>rkbest,
          rkbest = rk;
          relbest = rel;
          galgen = sub[1];
        )
      );
    );
    rel = relbest;
  );
  cyc = ZM_snf_group(rel,&U,&Ui);
  /* gal generators <--Ui--  --U--> snf generators */
  if(vprintab,print("galois group: ", cyc));
  [gal,galgen,U,Ui,cyc]
};

abgaloisfixedfield(abgal, K, flag, var) =
{
  my(gal,galgen,Ui,H);
  gal = abgal_gal(abgal);
  galgen = abgal_galgen(abgal);
  Ui = abgal_Ui(abgal);
  H = Ui*K;
  H = vector(#H[1,],i,factorback(galgen,H[,i]));
  galoisfixedfield(gal,H,flag,var)
};

\\[defining pol for subfield, embedding, relative pol]
getsubfield(abgal,K,{var='y},{reduce=0},{emb=1}) =
{
  my(res,redpol,isom,isomi,pol);
  if(emb,
    res = abgaloisfixedfield(abgal, K, 2, var);
    res = [res[1],res[2],res[3][1]*Mod(1,res[1])];
    if(reduce,
      [redpol,isomi] = polredbest(res[1],1);
      isom = lift(modreverse(isomi));
      pol = abgal_gal(abgal).pol;
      res = [redpol, Mod(subst(isom,variable(redpol),res[2]),pol),
        Mod(Pol(apply(c -> subst(lift(c),variable(c),isomi),Vec(res[3])),
        variable(res[3])),redpol)];
    );
  ,\\else
    res = abgaloisfixedfield(abgal, K, 1, var);
    res = subst(res,variable(res),var);
    if(reduce, res = polredbest(res));
    res = [res,0,0]
  );
  res
};

/*
  If p!=0, assume Gal=C*P with C cyclic and P p-Sylow, and write the coefficients
  of the relation coming from the p-Sylow.
*/
getsubfields(abgal,{var='y},{p=0},{reduce=0},{emb=1}) =
{
  my(sub=[], coef=0, r=0, v, valord, ord, C, coefbr = 0, K, cyc);
  cyc = abgal_cyc(abgal);
  if(p,
    r = #cyc;
    v = valuation(factorback(cyc),p);
    C = cyc[1]/p^valuation(cyc[1],p);
  ); 
  if(vprintab,print("generating list"));
  forvec(chi=vector(#cyc,i,[0,cyc[i]-1]),
    K = charker(cyc,chi);
    if(p,
      ord = charorder(cyc,chi);
      valord = valuation(ord,p);
      if(ord!=C*p^valord,next);
      if(valord==0,
        coefbr = (1-(p^r-1)/(p-1));
      ,\\chi_p != 0
        coefbr = if(chi%p!=0,1,1-p^(r-1));
      );
      coef = coefbr/p^(v-valord);
    );
    sub = concat(sub,[[K,[coef,coefbr]]]);
  );
  if(vprintab,print("eliminating duplicates #sub=", #sub));
  sub = Set(sub);
  if(vprintab,print("computing fields, #sub=", #sub));
  sub = vector(#sub,i,
    if(vprintab,print("i=",i,"/",#sub));
    concat(getsubfield(abgal,sub[i][1],var,reduce,emb),sub[i][2])
  ); 
  vecsort(sub,dat->[poldegree(dat[1]),poldisc(dat[1]),polsturm(dat[1])])
};

\\relation from Proposition 2.26
dualfunakurarelation(cyc) =
{
  my(sub=[],Lp,rp,p,H,ord,coef,coefbr,g,h,delta);
  g = vecprod(cyc);
  Lp = factor(cyc[1])[,1];
  rp = vector(#Lp,i,
    p = Lp[i];
    my(j=1);
    while(j<#cyc && cyc[j+1]%p==0, j++);
    j
  );
  forvec(chi=vector(#cyc,j,[0,cyc[j]-1]),
    H = charker(cyc,chi);
    ord = charorder(cyc,chi);
    h = g/ord;
    coefbr = 1;
    for(i=1,#Lp,
      p = Lp[i];
      if(ord%p==0,
        delta=1;
        for(j=1,#cyc,if(cyc[j]%p==0 && chi[j]%p, delta=0; break));
        if(delta, coefbr *= 1-p^(rp[i]-1));
      ,\\else
        coefbr *= -sum(j=1,rp[i]-1,p^j);
      );
    );
    if(coefbr,
      coef = coefbr/h;
      sub = concat(sub,[[H,[coef,coefbr]]]);
    );
  );
  if(vprintab,print("eliminating duplicates #sub=", #sub));
  sub = Set(sub);
  if(vprintab,print("#sub=", #sub));
  sub
};

\\G with two non-cyclic Sylows
\\return a list of subgroups and for each p, a Brauer relation over Z_p
\\relation from Theorem 2.27 (1)
abrelations(cyc) =
{
  my(Lp,p,cycp,cycc,Rels,rel,dat,H,coefbr,coefs);
  Rels = Map();
  Lp = factor(cyc[1])[,1]~;
  for(i=1,#Lp,
    p = Lp[i];
    [cycp,cycc] = splitcycsylow(cyc,p);
    rel = dualfunakurarelation(cycc);
    if(sum(j=1,#rel,rel[j][2])!=[1,1],error("incorrect relation! p=",p," sum=",sum(j=1,#rel,rel[j][2])));
    for(j=1,#rel,
      dat = rel[j];
      H = mathnfmodid(cycp[1]*matconcat([dat[1];vectorv(#cyc-#cycc)]),cyc);
      coefbr = dat[2][2];
      if(!mapisdefined(Rels,H,&coefs),coefs = vector(#Lp));
      coefs[i] = coefbr;
      mapput(Rels,H,coefs);
    );
  );
  Mat(Rels)
};

\\construct all subfields with Galois group C*P
\\with C cyclic and P a non-cyclic p-Sylow of G
getsubfields_CP(abgal,{var='y},{reduce=0},{emb=1}) =
{
  my(rels,cyc,perm);
  cyc = abgal_cyc(abgal);
  rels = abrelations(cyc);
  for(i=1,#rels[,1],
    if(vprintab,print("i=",i,"/",#rels[,1]));
    rels[i,1] = getsubfield(abgal,rels[i,1],var,reduce,emb);
  );
  perm = vecsort(rels[,1],dat->[poldegree(dat[1]),poldisc(dat[1]),polsturm(dat[1])],1);
  vecextract(rels,perm,[1,2])
};


/*     VII) Auxiliary functions for type 1         */

\\assume C*P case
myrootsofunity(pol,Lbnf,p) =
{
  my(w,kinf,ksup,bad,count=0,f,polcyc,nfcyc,q,n=poldegree(pol), disc);
  if(polsturm(pol)>0,return(2));
  w = lcm(vector(#Lbnf,i,Lbnf[i][2].tu[1])); \\from subfields: correct away from p
  kinf = valuation(w,p);
  if(n%(p-1), return(w));
  ksup = valuation(n,p)+1;
  if(ksup==kinf, return(w));
  disc = poldisc(pol);
  if(disc%p, return(w));
  bad = cheaprad(disc);
  w \= p^valuation(w,p);
  while(count<100,
    q = randomprime([10^3,10^5*expo(bad)]);
    if(!(bad%q), next);
    f = factormod(pol,q,1)[1,1];
    ksup = min(ksup, valuation(q^f-1,p));
    if(ksup==kinf, return(w*p^ksup));
    count++;
  );
  if(vprintab,print("kinf=",kinf," ksup=",ksup));
  while(ksup>kinf, \\probably only one iteration
    polcyc = polcyclo(p^ksup);
    nfcyc = nfinit(polcyc);
    if(nfisincl(nfcyc,pol), break);
    ksup--;
  );
  if(vprintab,print("k=",ksup));
  w*p^ksup;
};

maxspecialelt(pol) =
{
  my(s=2,nf);
  while(1,
    if(poldegree(pol)%2^(s-1),break);
    nf = nfinit(galoissubcyclo(2^(s+1),-1));
    if(!nfisincl(nf,pol), break);
    s++
  );
  s
};
\\detect Grunwald-Wang special case
isspecial(pol,Lbnf,p,w,r1,r2) =
{
  my(s,nf);
  if(p!=2, return(0));
  if(w%4==0, return(0));
  s = maxspecialelt(pol);
  if(r1>0, return(s));
  nf = nfinit(galoissubcyclo(2^(s+1),-Mod(5,2^(s+1))^(2^(s-2))));
  if(!nfisincl(nf,pol), return(s));
  0
};

\\identify restrictions of infinite places to subfield
resinfplaces(nf,r1,r2,ro2,emb) =
{
  my(g2,ro,resro2,perm,jbest,r,mul,d);
  g2 = variable(emb);
  emb = liftpol(emb);
  ro = nf.roots;
  resro2 = vector(#ro2, i, subst(emb,g2,ro2[i]));
  perm = vector(#ro2,i,
    r = resro2[i];
    if(imag(r)<0,r=conj(r));
    vecmin(vector(#ro,j,abs(ro[j]-r)),&jbest);
    jbest
  );
  d = poldegree(nf.pol);
  for(j=1,#ro,
    mul = #[i | i <- perm, i==j];
    if(mul != #ro2/d && mul != 2*#ro2/d,
      if(vprintab,print("#ro2: ", #ro2));
      if(vprintab,print("d=",d));
      if(vprintab,print("multiplicities: ", vector(#ro,j,#[i | i <- perm, i==j])));
      error("bug in resinfplaces [wrong multiplicities]")
    )
  );
  perm
};

\\log embedding in K from log embedding in the subfield:
unitmatrix(Lbnf,r1,r2,Lresinfpl) =
{
  Mat(concat(vector(#Lbnf,j,
    vector(unitrank(Lbnf[j][2]),j2,
      vectorv(r1+r2, i, if(Lresinfpl[j][i]<=Lbnf[j][2].r1 && i>r1,2,1)*real(Lbnf[j][2][3][Lresinfpl[j][i],j2]))
    )
  )))
};

increaseT(pol,Lsub,Lnf,T,prodT,nbT,prodS,bad,den,w) =
{
  my(boundq,q,fa,pr,data,prsub,decsub,nfsub,g,modu);
  boundq = 10^3*(1+expo(prodT)+expo(prodS)+expo(bad))*(nbT-#T)*poldegree(pol);
  modu = lcm(den,w);
  while(#T<nbT,
    q = 1+random(boundq)*modu;
    if(!(bad%q) || !(prodS%q) || !(prodT%q), next);
    if(!isprime(q), next);
    fa = factormod(pol,q,1);
    if(fa[1,1]>1, next);
    pr = myprimedec(pol,q)[1];
    sqrtn(Mod(1,q),den,&g);
    data = [pr,g,[(q-1)\den,den],vector(#Lnf,i,
      nfsub = Lnf[i][2];
      prsub = myidealnorm(Lsub[i][3],pr);
      decsub = idealprimedec(nfsub,q);
      prsub = decsub[whichprime(nfsub,decsub,prsub)];
      nfmodprinit(nfsub,prsub)
    )];
    T = concat(T,[data]);
    prodT *= q;
  );
  [T,prodT]
};

increaseS(pol,Lsub,Lbnf,S,prodS,qS,nbS,prodT,bad,disc) =
{
  my(fa,dec,data,nfsub,decsub,respr,prsub,pr);
  boundq = (expo(disc)^3+10^3)*(1+expo(prodT)+expo(prodS)+expo(bad))*(nbS-#S)*poldegree(pol);
  while(1,
    if(#S>=nbS,break);
    q = randomprime(boundq);
    if(!(bad%q) || !(prodT%q) || !(prodS%q), next);
    fa = factormod(pol,q,1);
    if(fa[1,1]>1, next);
    if(!isprime(q),next);
    if(vprintab,print1("    adding q=",q,"... "));
    dec = myprimedec(pol,q);
    data = [q,dec,vector(#Lbnf,i,
      nfsub = Lbnf[i][2];
      decsub = idealprimedec(nfsub,q);
      respr = vectorsmall(#dec,j,
        pr = dec[j];
        prsub = myidealnorm(Lsub[i][3],pr);
        whichprime(nfsub,decsub,prsub)
      );
      [decsub,respr]
    )];
    S = concat(S,[data]);
    qS = q;
    prodS *= q;
    if(vprintab,print("done adding q."));
  );
  [S,prodS,qS]
};

pruneS(S,relp) =
{
  my(m,n,keep=vector(#S),invpivot,i);
  if(vprintab,print("  pruning S"));
  [m,n] = matsize(relp);
  invpivot = vector(n);
  i=m;
  forstep(j=n,1,-1,
    while(i>0 && relp[i,j]==0, i--);
    if(i>0 && relp[i,j]==1, invpivot[j]=i);
  );
  invpivot = Set(invpivot);
  i=1;
  for(k=1,#S,
    for(l=1,#S[k][2],
      if(!vecsearch(invpivot,i),keep[k]=1);
      i++;
    );
  );
  if(vprintab,print("  keep: ", keep));
  [S[i] | i <- [1..#S], keep[i]]
};

\\fundamental units first, then root of unity, then S-units
valresSunits(S,T,Lbnf,subunits,unitcorresp,izp,den) =
{
  my(Mval,Mres,su,Ssub,bnfsub,uc,Stot,Lsu = vector(#Lbnf));
  Stot = sum(i=1,#S,#S[i][2]);
  if(vprintab,print1("  S-units of subfields..."));
  if(#S!=0,
    for(i=1,#Lbnf,
      bnfsub = Lbnf[i][2];
      Ssub = concat(vector(#S,j,S[j][3][i][1]));
      Lsu[i]=mysunits(bnfsub,Ssub,den)
    );
  );
  if(vprintab,print("done."));

  if(vprintab,print1("  valuations..."));
  Mval = matconcat([
    \\units: valuations 0
    matrix(Stot,matsize(subunits)[2]+1),
    \\sunits
    if(#S==0,[;],
      Mat(concat(vector(#Lbnf,i,
        bnfsub = Lbnf[i][2];
        su = Lsu[i];
        vector(#su[1],j,valfromsub(bnfsub,su[1][j],i,S,[su[2][,j],su[3]]))
      )))
    )
  ]);
  if(vprintab,print("done."));

  if(vprintab,print1("  residues and DLs..."));
  Lunits = [bnfunits(bnfsub[2])[1] | bnfsub <- Lbnf];
  Mres = matconcat([
    \\fundamental units
    Mat(vector(#unitcorresp,i,
      uc = unitcorresp[i];
      bnfsub = Lbnf[uc[1]][2];
      resfromsub(bnfsub,Lunits[uc[1]][uc[2]],uc[1],T)
    ))*subunits,
    \\root of unity
    resfromsub(Lbnf[izp][2],Lbnf[izp][2].tu[2],izp,T),
    \\sunits
    if(#S==0,[;],
      Mat(concat(vector(#Lbnf,i,
        bnfsub = Lbnf[i][2];
        su = Lsu[i];
        vector(#su[1],j,resfromsub(bnfsub,su[1][j],i,T))
      )))
    )
  ]);
  if(vprintab,print("done."));

  if(vprintab,print("  Stot: ", Stot));
  if(vprintab,print("  size Mval: ", matsize(Mval)));
  if(vprintab,print("  size Mres: ", matsize(Mres)));
  [Mval,Mres]
};

saturate(S,T,Lbnf,subunits,unitcorresp,izp,den,r1,r2,p,expobound) =
{
    my(Mval,Mres,Mrestot,Mresu,K,indexu,Mvalall);
    \\compute valuations and residues of S-units of subfields
    if(vprintab,print("  computing valuations and residues"));
    [Mval,Mres] = valresSunits(S,T,Lbnf,subunits,unitcorresp,izp,den^2);
    Mrestot = Mres;
    Mres = Mres%den;
    Mresu = Mres[,1..r1+r2];
    \\compute kernels for den-th powers
    if(vprintab,print("  Mresu: ", matsize(Mresu)));
    K = kermodprimepow(Mresu,den);
    if(vprintab,print("  kernel in units: ", matsize(K)));
    \\mod out root of unity component since already counted in w
    if(#K>0, K = matimagemod(K[1..r1+r2-1,],den));
    indexu = cardmod(K,den);
    if(vprintab,print("  indexu = ", indexu, " = ", p, "^", valuation(indexu,p)));
    K = kermodprimepow(matconcat([Mval;Mres]),den);
    if(vprintab,print("  kernel in S-units: ", matsize(K)));
    Mvalall = matconcat([Mval,(Mval*K)\den]);
    relp = matimagemod(Mvalall,expobound);
    if(vprintab,print("  relp: dims=", matsize(relp)));
    hp = expobound^matsize(Mval)[1]\cardmod(relp,expobound);
    if(vprintab,print("  hp = ", hp, " = ", p, "^", valuation(hp,p)));
    [indexu,relp,hp,Mrestot,Mval,K]
};


/*     VIII) Main functions         */

\\main function
abelianbnfinit(pol,{bad=1},{Lbnfdata=[]}) =
{
  my(gal,abgal,cycnc,p);
  if(vprintab,print("\npol=", pol, "\nstarting abelianbnfinit"));
  if(poldegree(pol)==1, return(abelianbnfinit0(pol)));
  abgal = abgaloisanalysis(pol);
  cyc = abgal_cyc(abgal);
  cycnc = vecprod(cyc[2..-1]);
  if(cycnc==1,
    return(abelianbnfinit0(pol))
  ,isprimepower(cycnc,&p),
    return(abelianbnfinit1(pol,abgal,p,bad,Lbnfdata))
  ,/*else: at least two non-cyclic Sylows*/
    return(abelianbnfinit2(pol,abgal,bad,Lbnfdata))
  );
};

\\denominator 1 relation
abelianbnfinit2(pol,abgal,bad0,Lbnfdata) =
{
  my(LsubC,LsubCP,Lbnf,bad,d,rel,abbnf,brauer,Lp,Lbnfdata2,polsub,k);
  if(vprintab,print("at least 2 non-cyclic Sylows"));
  LsubCP = getsubfields_CP(abgal,'y,0,0);
  Lp = factor(poldegree(pol))[,1]~;
  brauer = LsubCP[,2]~;
  brauer = vector(#Lp,i,vector(#brauer,j,brauer[j][i]));
  LsubCP = LsubCP[,1]~;
  LsubC = getsubfields(abgal,'z,,1,0);
  if(vprintab,print("computing bnf1"));
  Lbnfdata2 = vector(#LsubC);
  k=0;
  for(i=1,#LsubC,
    polsub = LsubC[i][1];
    if(poldegree(polsub)==1 || fetchbnf(polsub,Lbnfdata),next);
    if(vprintab,print1("\ni=",i,"/",#LsubC," "));
    abbnf = abelianbnfinit(polsub,Lbnfdata); 
    k++;
    Lbnfdata2[k] = [nfhash(polsub),abbnf];
  );
  Lbnfdata = vecsort(concat(Lbnfdata,Lbnfdata2[1..k]),1);
  Lbnf = vector(#LsubCP,i,
    if(vprintab,print1("\ni=",i,"/",#LsubCP," "));
    abbnf = abelianbnfinit(LsubCP[i][1],,Lbnfdata);
    abbnf
  );
  d = lcm(apply(getexpo,Lbnf));
  if(vprintab,print("\n", apply(getcyc,Lbnf)," d=", d));
  [2,Brauerclgp(Lbnf,Lp,brauer,d),Lbnf,pol,Lbnfdata,abgal_cyc(abgal)];
};

\\p-power denominator relation
abelianbnfinit1(pol,abgal,p,bad0,Lbnfdata) =
{
  my(Lsub,Lbnf,d,dc,bad,hc,coefs,den,coefsbr,w,hR,r1,r2,ro,Lresinfpl,units,
    Ubasis,normemb,normindices,R00,R0,check,nbT,nbS,qS,prodT,prodS,T,S,Mval,K,
    indexu,R1,hp,expobound,relp,subunits,denemb,unitcorresp,izp,snfc,snfp,
    cyctot,disc,f,spe,Up,Mrestot,abbnf,matextru,detubasis,detectoo=0,sep,rosub,
    bprec,emb,LMi,precu,prechr,r0,precR0,checksum);
  if(vprintab,print("non-cyclic galois group"));
  Lsub = getsubfields(abgal,'z,p,1);
  if(vprintab,print("number of subfields: ", #Lsub));
  coefs=[x[4] | x <- Lsub];
  checksum = sum(i=1,#coefs,coefs[i]);
  if(vprintab,print("coefs norm rel: ", coefs, " sum=", checksum, " (must be 1)"));
  if(checksum!=1, error("incorrect norm relation!"));
  den = denominator(coefs);
  coefs *= den;
  if(vprintab,print("den=", den));
  coefsbr = [x[5] | x <- Lsub];
  checksum = sum(i=1,#coefsbr,coefsbr[i]);
  if(vprintab,print("coefs brauer rel: ", coefsbr, " sum=", checksum, " (must be 1)"));
  if(checksum!=1, error("incorrect Brauer relation!"));

  Lbnf = vector(#Lsub,i,
    if(vprintab,print1("\ni=",i,"/",#Lsub," "));
    abbnf = fetchbnf(Lsub[i][1],Lbnfdata);
    if(abbnf,
      if(vprintab,print("found in list."));
      if(abbnf[2]!=x,
        Lsub[i][1] = getpol(abbnf[1]);
        Lsub[i][2] = subst(abbnf[2],variable(abbnf[2]),Lsub[i][2]);
        Lsub[i][3] = Mod(Pol(apply(c ->
          subst(lift(c),variable(c),Mod(abbnf[3],Lsub[i][1])), Vec(Lsub[i][3])),
          variable(Lsub[i][3])),Lsub[i][1]);
      );
      abbnf[1]
    ,\\else
      abbnf = abelianbnfinit0(Lsub[i][1]);
      abbnf
    )
  );

  disc = factorback(vector(#Lbnf,i,Lbnf[i][2].disc),coefsbr);
  f = poldisc(pol)/disc;
  check = issquare(f,&f);
  if(vprintab,printf("disc = %Ps\nf = %Ps\n", disc, f));
  if(!check, error("disc pol / disc is not square!"));

  d = lcm(apply(getexpo,Lbnf));
  dc = d/p^valuation(d,p);
  if(vprintab,print(apply(getcyc,Lbnf)," d=", d, " dc=", dc));
  bad = vector(#Lsub,i,poldisc(Lsub[i][1]));
  bad = concat([bad,[disc,f]]);
  bad = cheaprad(lcm(concat(apply(cheaprad,bad),bad0)));
  if(vprintab,print("bad=", bad));

  \\compute p'-part of class group
  snfc = Brauerclgp(Lbnf,[p],[coefsbr],dc);
  hc = vecprod(snfc);

  \\compute number of roots of unity
  w = myrootsofunity(pol,Lbnf,p);
  if(vprintab,print("w=", w));

  \\compute h*R using subfields
  prechr = precbrauerhr(Lbnf,coefsbr);
  if(prechr>bitprecision(Lbnf),
    if(vprintab,print("increasing prec! bitprec(Lbnf):",bitprecision(Lbnf),
      " prechr=",prechr));
    localbitprec(prechr);
    for(i=1,#Lbnf,Lbnf[i][2]=nfnewprec(Lbnf[i][2]));
  );
  hR = factorback(vector(#Lbnf,i,my(bnf=Lbnf[i][2]); bnf.no*bnf.reg/bnf.tu[1]),coefsbr)*w;
  if(vprintab,print("hR=", precision(hR,30)));
  [r1,r2,ro] = myinfplaces(pol);

  \\detect special case of GW
  spe = isspecial(pol,Lbnf,p,w,r1,r2);
  if(vprintab,print("special field? ", spe));
  if(spe && (den^2)%(2^(spe+1))==0,
    detectoo = 1;
    warning("Grunwald-Wang special field: infinite loop possible"));

  \\compute restrictions of infinite places
  bprec = bitprecision(ro);
  for(i=1,#Lbnf,
    rosub = Lbnf[i][2].nf.roots;
    if(#rosub==1,next);
    sep = vecmax([-expo(abs(rosub[i]-rosub[j])) | i <- [1..#rosub]; j <- [i+1..#rosub]]);
    sep = max(0,2+sep)+5;
    emb = liftpol(Lsub[i][2]);
    sep += vecmax([preceval(emb,r) | r <- ro]);
    if(sep>bprec,bprec=sep);
  );
  if(vprintab,print("resinfplaces: bprec=",bprec));
  localbitprec(bprec);
  [r1,r2,ro] = myinfplaces(pol);
  Lresinfpl = vector(#Lsub,i,resinfplaces(Lbnf[i][2].nf,r1,r2,ro,Lsub[i][2]));

  \\compute regulator of image of units from subfields
  units = unitmatrix(Lbnf,r1,r2,Lresinfpl);
  \\compute a basis of the group of subfield units
  LMi = vector(#Lbnf,i,preMi(Lbnf[i][2]));
  precu = vecmax(vector(#LMi,i,expo(norml2(LMi[i]))));
  precu += vecmax(abs(apply(exponent,units)));
  precu += 2*expo(r1+2*r2) + 5;
  if(precu>bitprecision(Lbnf),
    if(vprintab,print("increasing prec! bitprec(Lbnf):",bitprecision(Lbnf),
      " precu=",precu));
    localbitprec(precu);
    for(i=1,#Lbnf,Lbnf[i][2]=nfnewprec(Lbnf[i][2]));
    units = unitmatrix(Lbnf,r1,r2,Lresinfpl);
    LMi = vector(#Lbnf,i,preMi(Lbnf[i][2]));
  );
  normemb = Mat(vector(#units,j,
    \\norms to all subfields
    concat(vector(#Lsub,i,
      my(lognormu = vectorv(unitrank(Lbnf[i][2])+1));
      for(k=1,r1+r2,lognormu[Lresinfpl[i][k]]+=units[k,j]);
      logisunit(LMi[i],lognormu)
    ));
  ));
  \\compute independent subset of the units
  forprime(ell=100,oo,
    normindices = matindexrank(Mod(normemb,ell))[2];
    if(#normindices == r1+r2-1,
      if(vprintab,print("prime used for independent subset of units: ", ell));
    break);
  );
  if(vprintab,print("indices: ", normindices));
  matextru = vecextract(units,[1..r1+r2-1],normindices);
  r0 = matsize(matextru)[1];
  R00 = abs(matdet(matextru));
  precR0 = precdirectR0(matextru,r0,R00);
  if(precR0>bitprecision(Lbnf),
    if(vprintab,print("increasing prec! bitprec(Lbnf):",bitprecision(Lbnf),
      " precR0=",precR0));
    localbitprec(precR0);
    for(i=1,#Lbnf,Lbnf[i][2]=nfnewprec(Lbnf[i][2]));
    units = unitmatrix(Lbnf,r1,r2,Lresinfpl);
    matextru = vecextract(units,[1..r1+r2-1],normindices);
    R00 = abs(matdet(matextru));
  );

  normemb = matsolve(vecextract(normemb,normindices),normemb);
  denemb = denominator(normemb);
  if(vprintab,print("denemb: ", denemb));
  Ubasis = mathnf(denemb*normemb)/denemb;
  detubasis = matdet(Ubasis);
  if(vprintab,print("index: ", detubasis));
  R0 = R00*detubasis;
  if(vprintab,print("R0=", precision(R0,30))); \\regulator of units of subfields before saturation

  \\basis of units from subfields mod n, expressed in terms of the units of subfields
  matimagemod(denemb*normemb,den*p^valuation(denemb,p),&subunits);
  if(vprintab,print("normemb: ", matsize(normemb)));
  if(vprintab,print("subunits:", matsize(subunits)));
  unitcorresp = concat(vector(#Lsub,i,vector(unitrank(Lbnf[i][2]),j,[i,j])));

  \\subfield containing the root of unity of maximal p-power order
  vecmax(vector(#Lbnf,i,valuation(Lbnf[i][2].tu[1],p)),&izp);

  \\bound on the exponent of Cl(K)_p
  expobound=p^vecmax(vector(#Lbnf,i,if(Lbnf[i][2].no==1,0,valuation(Lbnf[i][2].cyc[1],p))));
  expobound *= den;
  expobound *= p;
  if(vprintab,print("expobound = ",expobound," = ",p,"^",valuation(expobound,p)));
  if(vprintab,print("den = ", den, " = ",p,"^",valuation(den,p)));

  \\initialise a T that hopefully defines an embedding of the Selmer group Sel_den(K)
  if(vprintab,print1("initialising T..."));
  T = [];
  nbT = r1+r2+10;
  [T,prodT] = increaseT(pol,Lsub,Lbnf,T,1,nbT,1,bad,den^2,w);
  if(vprintab,print("done."));

  \\initialise an S that hopefully generates Cl(K)_p
  S = [];
  nbS = 0;
  prodS = 1;
  qS = 1;

  check = oo;
  while(1,
    if(vprintab,print("\nnbT=", nbT));
    if(vprintab,print("nbS=", nbS, " total nb primes=", sum(i=1,#S,#S[i][2])));
    [indexu,relp,hp,Mrestot,Mval,K] = saturate(S,T,Lbnf,subunits,unitcorresp,
      izp,den,r1,r2,p,expobound);
    \\compute corresponding hR and check
    R1 = R0/indexu;
    check=hR/(R1*hc*hp);
    if(vprintab,print("  check: ",p,"^",round(log(check)/log(p))," (",precision(check,30),
    " ~ ", bestappr(check),")"));

    if(check<0.7, error("bug: non-integral check: ", check));
    if(check>1.4,
      if(detectoo>0,
        if(check<3, detectoo++);
        if(detectoo>3, warning("probable infinite loop detected."));
      );
      nbT = ceil(1.05*nbT);
      nbS = round(1.3*nbS);
      if(nbS==0, nbS=2);
      nbS += 1;
      S = pruneS(S,relp);
      if(vprintab,print("  increasing S... "));
      [S,prodS,qS] = increaseS(pol,Lsub,Lbnf,S,prodS,qS,nbS,prodT,bad,disc);
      if(vprintab,print1("  done.\n  increasing T... "));
      [T,prodT] = increaseT(pol,Lsub,Lbnf,T,prodT,nbT,prodS,bad,den^2,w);
      if(vprintab,print("done."));
    ,\\else
      break
    );
  );

  \\compute class group structure
  snfp = ZM_snf_group(relp,&Up);
  if(vprintab,print("snfp = ",snfp));
  cyctot = mergecyc(snfc,snfp);
  if(vprintab,print("cyctot = ",cyctot));

  [ 1,          \\1
    cyctot,     \\2
    pol,        \\3
    Lsub,       \\4
    Lbnf,       \\5
    snfc,       \\6
    0,          \\7
    dc,         \\8
    snfp,       \\9
    0,          \\10
    0,          \\11
    0,          \\12
    Up,         \\13
    coefs,      \\14
    S,          \\15
    T,          \\16
    den,        \\17
    disc,       \\18
    p,          \\19
    0,          \\20
    0,          \\21
    lcm([bad,prodS,prodT]), \\22
    coefsbr,    \\23
    [r1,r2],    \\24
    w,          \\25
    matextru,   \\26
    R00,        \\27
    detubasis,  \\28
    Lresinfpl,  \\29
    normindices,\\30
    hp*hc/indexu    \\31
    ]
};

\\base case: cyclic Galois group
abelianbnfinit0(pol) =
{
  my(nf,bnf,prec=default(realprecision));
  localprec(100);
  if(vprintab,print("base case, starting bnfinit"));
  nf = easynf(pol);
  if(vprintab,print("sig=",nf.sign," disc=",nf.disc," ~ 2^",exponent(nf.disc)));
  bnf = bnfinit(nf,1);
  if(vprintab,print("cyc=",bnf.cyc));
  if(vprintab,print1("increasing precision..."));
  localprec(max(100,prec));
  bnf = nfnewprec(bnf);
  if(vprintab,print("done."));
  [0,bnf];
};


/*     IX) Certification         */

fastcert = 0;

abelianbnfcertify(abbnf) =
{
  my(L=[]);
  if(abbnf[1]==2,
    L = abbnf[3];
    for(i=1,#L,
      if(vprintab,print("certify: i=",i,"/",#L));
      if(!abelianbnfcertify(L[i]),return(0))
    );
  ,abbnf[1]==1,
    L = abbnf[5];
    for(j=1,#L,
      if(vprintab,print("certify:   j=",j,"/",#L));
      if(!abelianbnfcertify(L[j]),return(0))
    );
  ,/*else: 0*/
    if(#Lp,
      return(bnfcertify(abbnf[2]))
    ,\\else
      return(bnfcertify(abbnf[2]))
    )
  );
  1
};


/*     X) Example generation    */

\\assume linearly disjoint
veccompositum(L) =
{
  my(pol1,pol2,n=#L,pol);
  if(n==1, return(L[1]));
  pol1 = veccompositum(L[1..n\2]);
  pol2 = veccompositum(L[n\2+1..-1]);
  pol = polcompositum(pol1,pol2,2);
  if(poldegree(pol)<30,pol=polredbest(pol));
  pol
};

multiquad(L,var='x) =
{
  veccompositum([var^2-a | a <- L])
};

kummer(n,L,var='x) =
{
  veccompositum(concat([var^n-a | a <- L],polcyclo(n,var)))
};

/*
multiquad(L,var='x) =
{
  my(pol1,pol2,n=#L,pol);
  if(n==1, return(var^2-L[1]));
  pol1 = multiquad(L[1..n\2],var);
  pol2 = multiquad(L[n\2+1..-1],var);
  pol = polcompositum(pol1,pol2,2);
  if(poldegree(pol)<30,pol=polredbest(pol));
  pol
};
*/