Code coverage tests

This page documents the degree to which the PARI/GP source code is tested by our public test suite, distributed with the source distribution in directory src/test/. This is measured by the gcov utility; we then process gcov output using the lcov frond-end.

We test a few variants depending on Configure flags on the pari.math.u-bordeaux.fr machine (x86_64 architecture), and agregate them in the final report:

The target is to exceed 90% coverage for all mathematical modules (given that branches depending on DEBUGLEVEL or DEBUGMEM are not covered). This script is run to produce the results below.

LCOV - code coverage report
Current view: top level - basemath - kummer.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.12.1 lcov report (development 25685-2ac6d95b9c) Lines: 853 862 99.0 %
Date: 2020-08-06 06:06:08 Functions: 61 61 100.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : /*******************************************************************/
      15             : /*                                                                 */
      16             : /*                      KUMMER EXTENSIONS                          */
      17             : /*                                                                 */
      18             : /*******************************************************************/
      19             : #include "pari.h"
      20             : #include "paripriv.h"
      21             : 
      22             : typedef struct {
      23             :   GEN R; /* nf.pol */
      24             :   GEN x; /* tau ( Mod(x, R) ) */
      25             :   GEN zk;/* action of tau on nf.zk (as t_MAT) */
      26             : } tau_s;
      27             : 
      28             : typedef struct {
      29             :   GEN polnf, invexpoteta1, powg;
      30             :   tau_s *tau;
      31             :   long m;
      32             : } toK_s;
      33             : 
      34             : typedef struct {
      35             :   GEN R; /* ZX, compositum(P,Q) */
      36             :   GEN p; /* QX, Mod(p,R) root of P */
      37             :   GEN q; /* QX, Mod(q,R) root of Q */
      38             :   long k; /* Q[X]/R generated by q + k p */
      39             :   GEN rev;
      40             : } compo_s;
      41             : 
      42             : /* REDUCTION MOD ell-TH POWERS */
      43             : /* make b integral by multiplying by t in (Q^*)^ell */
      44             : static GEN
      45        1603 : reduce_mod_Qell(GEN nf, GEN b, GEN gell)
      46             : {
      47             :   GEN c;
      48        1603 :   b = nf_to_scalar_or_basis(nf, b);
      49        1603 :   b = Q_primitive_part(b, &c);
      50        1603 :   if (c)
      51             :   {
      52         952 :     GEN d, fa = Q_factor_limit(c, 1000000);
      53         952 :     gel(fa,2) = FpC_red(gel(fa,2), gell);
      54         952 :     d = factorback(fa);
      55         952 :     b = typ(b) == t_INT? mulii(b,d): ZC_Z_mul(b, d);
      56             :   }
      57        1603 :   return b;
      58             : }
      59             : 
      60             : static GEN
      61        1603 : reducebeta(GEN bnfz, GEN b, GEN gell)
      62             : {
      63        1603 :   GEN t, cb, fu, nf = bnf_get_nf(bnfz);
      64        1603 :   long ell = itou(gell);
      65             : 
      66        1603 :   if (DEBUGLEVEL>1) err_printf("reducing beta = %Ps\n",b);
      67        1603 :   b = reduce_mod_Qell(nf, b, gell);
      68        1603 :   t = idealredmodpower(nf, b, ell, 1000000);
      69        1603 :   if (!isint1(t)) b = nfmul(nf, b, nfpow_u(nf, t, ell));
      70        1603 :   if (DEBUGLEVEL>1) err_printf("beta reduced via ell-th root = %Ps\n",b);
      71        1603 :   b = Q_primitive_part(b, &cb);
      72        1603 :   if (cb)
      73             :   {
      74         488 :     GEN y = nfroots(nf, gsub(monomial(gen_1, ell, fetch_var_higher()),
      75             :                          basistoalg(nf,b)));
      76         488 :     if (lg(y) != 1) b = cb;
      77         488 :     delete_var();
      78             :   }
      79        1603 :   if (b != cb && (fu = bnf_build_cheapfu(bnfz)))
      80             :   { /* log. embeddings of fu^ell */
      81        1428 :     GEN elllogfu = RgM_Rg_mul(real_i(bnf_get_logfu(bnfz)), gell);
      82        1428 :     long prec = nf_get_prec(nf);
      83             :     for (;;)
      84          14 :     {
      85        1442 :       GEN z = nflogembed(nf, b, NULL, prec);
      86        1442 :       if (z)
      87             :       {
      88        1428 :         GEN y, ex = RgM_Babai(elllogfu, z);
      89        1428 :         if (ex)
      90             :         {
      91        1428 :           if (ZV_equal0(ex)) break;
      92         351 :           y = nffactorback(nf, fu, RgC_Rg_mul(ex,gell));
      93         351 :           b = nfdiv(nf, b, y); break;
      94             :         }
      95             :       }
      96          14 :       prec = precdbl(prec);
      97          14 :       if (DEBUGLEVEL) pari_warn(warnprec,"reducebeta",prec);
      98          14 :       nf = nfnewprec_shallow(nf,prec);
      99             :     }
     100        1428 :     if (cb) b = gmul(b, cb);
     101             :   }
     102        1603 :   if (DEBUGLEVEL>1) err_printf("beta LLL-reduced mod U^l = %Ps\n",b);
     103        1603 :   return b;
     104             : }
     105             : 
     106             : struct rnfkummer
     107             : {
     108             :   GEN bnfz, cycgenmod, u, vecC, tQ, vecW;
     109             :   ulong mgi, g, ell;
     110             :   long rc;
     111             :   compo_s COMPO;
     112             :   tau_s tau;
     113             :   toK_s T;
     114             : };
     115             : 
     116             : /* set kum->tau; compute Gal(K(\zeta_l)/K) */
     117             : static void
     118         553 : get_tau(struct rnfkummer *kum)
     119             : { /* compute action of tau: q^g + kp */
     120         553 :   compo_s *C = &kum->COMPO;
     121         553 :   GEN U = RgX_add(RgXQ_powu(C->q, kum->g, C->R), RgX_muls(C->p, C->k));
     122         553 :   kum->tau.x  = RgX_RgXQ_eval(C->rev, U, C->R);
     123         553 :   kum->tau.R  = C->R;
     124         553 :   kum->tau.zk = nfgaloismatrix(bnf_get_nf(kum->bnfz), kum->tau.x);
     125         553 : }
     126             : 
     127             : static GEN tauofvec(GEN x, tau_s *tau);
     128             : static GEN
     129       82581 : tauofelt(GEN x, tau_s *tau)
     130             : {
     131       82581 :   switch(typ(x))
     132             :   {
     133        4263 :     case t_INT: case t_FRAC: return x;
     134       75574 :     case t_COL: return RgM_RgC_mul(tau->zk, x);
     135        2744 :     case t_MAT: return mkmat2(tauofvec(gel(x,1), tau), gel(x,2));
     136             :     default: pari_err_TYPE("tauofelt",x); return NULL;/*LCOV_EXCL_LINE*/
     137             :   }
     138             : }
     139             : static GEN
     140        3143 : tauofvec(GEN x, tau_s *tau)
     141             : {
     142             :   long i, l;
     143        3143 :   GEN y = cgetg_copy(x, &l);
     144       78976 :   for (i=1; i<l; i++) gel(y,i) = tauofelt(gel(x,i), tau);
     145        3143 :   return y;
     146             : }
     147             : /* [x, tau(x), ..., tau^(m-1)(x)] */
     148             : static GEN
     149        1344 : powtau(GEN x, long m, tau_s *tau)
     150             : {
     151        1344 :   GEN y = cgetg(m+1, t_VEC);
     152             :   long i;
     153        1344 :   gel(y,1) = x;
     154        2352 :   for (i=2; i<=m; i++) gel(y,i) = tauofelt(gel(y,i-1), tau);
     155        1344 :   return y;
     156             : }
     157             : /* x^lambda */
     158             : static GEN
     159        1778 : lambdaofelt(GEN x, toK_s *T)
     160             : {
     161        1778 :   tau_s *tau = T->tau;
     162        1778 :   long i, m = T->m;
     163        1778 :   GEN y = trivial_fact(), powg = T->powg; /* powg[i] = g^i */
     164        4046 :   for (i=1; i<m; i++)
     165             :   {
     166        2268 :     y = famat_mulpows_shallow(y, x, uel(powg,m-i+1));
     167        2268 :     x = tauofelt(x, tau);
     168             :   }
     169        1778 :   return famat_mul_shallow(y, x);
     170             : }
     171             : static GEN
     172         672 : lambdaofvec(GEN x, toK_s *T)
     173             : {
     174             :   long i, l;
     175         672 :   GEN y = cgetg_copy(x, &l);
     176        2107 :   for (i=1; i<l; i++) gel(y,i) = lambdaofelt(gel(x,i), T);
     177         672 :   return y;
     178             : }
     179             : 
     180             : static GEN
     181         455 : tauofideal(GEN id, tau_s *tau)
     182             : {
     183         455 :   return ZM_hnfmodid(RgM_mul(tau->zk, id), gcoeff(id, 1,1));
     184             : }
     185             : 
     186             : static int
     187         938 : prconj(GEN P, GEN Q, tau_s *tau)
     188             : {
     189         938 :   GEN p = pr_get_p(P), x = pr_get_gen(P);
     190             :   for(;;)
     191             :   {
     192        2037 :     if (ZC_prdvd(x,Q)) return 1;
     193        1491 :     x = FpC_red(tauofelt(x, tau), p);
     194        1491 :     if (ZC_prdvd(x,P)) return 0;
     195             :   }
     196             : }
     197             : static int
     198        3843 : prconj_in_list(GEN S, GEN P, tau_s *tau)
     199             : {
     200             :   long i, l, e, f;
     201             :   GEN p, x;
     202        3843 :   if (!tau) return 0;
     203        1729 :   p = pr_get_p(P); x = pr_get_gen(P);
     204        1729 :   e = pr_get_e(P); f = pr_get_f(P); l = lg(S);
     205        2345 :   for (i = 1; i < l; i++)
     206             :   {
     207        1162 :     GEN Q = gel(S, i);
     208        1162 :     if (equalii(p, pr_get_p(Q)) && e == pr_get_e(Q) && f == pr_get_f(Q))
     209         938 :       if (ZV_equal(x, pr_get_gen(Q)) || prconj(gel(S,i), P, tau)) return 1;
     210             :   }
     211        1183 :   return 0;
     212             : }
     213             : 
     214             : /* >= ell */
     215             : static long
     216        2359 : get_z(GEN pr, long ell) { return ell * (pr_get_e(pr) / (ell-1)); }
     217             : /* zeta_ell in nfz */
     218             : static void
     219        1603 : list_Hecke(GEN *pSp, GEN *pvsprk, GEN nfz, GEN fa, GEN gell, tau_s *tau)
     220             : {
     221        1603 :   GEN P = gel(fa,1), E = gel(fa,2), faell, Sl, S, Sl1, Sl2, Vl, Vl2;
     222        1603 :   long i, l = lg(P), ell = gell[2];
     223             : 
     224        1603 :   S  = vectrunc_init(l);
     225        1603 :   Sl1= vectrunc_init(l);
     226        1603 :   Sl2= vectrunc_init(l);
     227        1603 :   Vl2= vectrunc_init(l);
     228        3549 :   for (i = 1; i < l; i++)
     229             :   {
     230        1946 :     GEN pr = gel(P,i);
     231        1946 :     if (!equaliu(pr_get_p(pr), ell))
     232        1673 :     { if (!prconj_in_list(S,pr,tau)) vectrunc_append(S,pr); }
     233             :     else
     234             :     { /* pr | ell */
     235         273 :       long a = get_z(pr, ell) + 1 - itou(gel(E,i));
     236         273 :       if (!a)
     237          28 :       { if (!prconj_in_list(Sl1,pr,tau)) vectrunc_append(Sl1, pr); }
     238         245 :       else if (a != 1 && !prconj_in_list(Sl2,pr,tau))
     239             :       {
     240          49 :         vectrunc_append(Sl2, pr);
     241          49 :         vectrunc_append(Vl2, log_prk_init(nfz, pr, a, gell));
     242             :       }
     243             :     }
     244             :   }
     245        1603 :   faell = idealprimedec(nfz, gell); l = lg(faell);
     246        1603 :   Vl = vectrunc_init(l);
     247        1603 :   Sl = vectrunc_init(l);
     248        3969 :   for (i = 1; i < l; i++)
     249             :   {
     250        2366 :     GEN pr = gel(faell,i);
     251        2366 :     if (!tablesearch(P, pr, cmp_prime_ideal) && !prconj_in_list(Sl, pr, tau))
     252             :     {
     253        2086 :       vectrunc_append(Sl, pr);
     254        2086 :       vectrunc_append(Vl, log_prk_init(nfz, pr, get_z(pr,ell), gell));
     255             :     }
     256             :   }
     257        1603 :   *pvsprk = shallowconcat(Vl2, Vl); /* divide ell */
     258        1603 :   *pSp = shallowconcat(S, Sl1);
     259        1603 : }
     260             : 
     261             : /* Return a Flm, sprk mod pr^k, pr | ell, k >= 2 */
     262             : static GEN
     263        2135 : logall(GEN nf, GEN v, long lW, long mgi, GEN gell, GEN sprk)
     264             : {
     265        2135 :   long i, ell = gell[2], l = lg(v);
     266        2135 :   GEN M = cgetg(l,t_MAT);
     267        9870 :   for (i = 1; i < l; i++)
     268             :   {
     269        7735 :     GEN c = log_prk(nf, gel(v,i), sprk, gell); /* ell-rank = #c */
     270        7735 :     c = ZV_to_Flv(c, ell);
     271        7735 :     if (i < lW) c = Flv_Fl_mul(c, mgi, ell);
     272        7735 :     gel(M,i) = c;
     273             :   }
     274        2135 :   return M;
     275             : }
     276             : static GEN
     277        1603 : matlogall(GEN nf, GEN v, long lW, long mgi, GEN gell, GEN vsprk)
     278             : {
     279        1603 :   GEN M = NULL;
     280        1603 :   long i, l = lg(vsprk);
     281        3738 :   for (i = 1; i < l; i++)
     282        2135 :     M = vconcat(M, logall(nf, v, lW, mgi, gell, gel(vsprk,i)));
     283        1603 :   return M;
     284             : }
     285             : 
     286             : /* id = (b) prod_{i <= rc} bnfz.gen[i]^v[i] (mod K^*)^ell,
     287             :  * - i <= rc: gen[i]^cyc[i] = (cycgenmod[i]); ell | cyc[i]
     288             :  * - i > rc: gen[i]^(u[i]*cyc[i]) = (cycgenmod[i]); u[i] cyc[i] = 1 mod ell */
     289             : static void
     290        1617 : isprincipalell(GEN bnfz, GEN id, GEN cycgenmod, ulong ell, long rc,
     291             :                GEN *pv, GEN *pb)
     292             : {
     293        1617 :   long i, l = lg(cycgenmod);
     294        1617 :   GEN y = bnfisprincipal0(bnfz, id, nf_FORCE|nf_GENMAT);
     295        1617 :   GEN v = ZV_to_Flv(gel(y,1), ell), b = gel(y,2);
     296        1862 :   for (i = rc+1; i < l; i++)
     297         245 :     b = famat_mulpows_shallow(b, gel(cycgenmod,i), v[i]);
     298        1617 :   setlg(v,rc+1); *pv = v; *pb = b;
     299        1617 : }
     300             : 
     301             : static GEN
     302        1603 : compute_beta(GEN X, GEN vecWB, GEN ell, GEN bnfz)
     303             : {
     304        1603 :   GEN be = famat_reduce(famatV_zv_factorback(vecWB, X));
     305        1603 :   if (typ(be) == t_MAT)
     306             :   {
     307        1561 :     gel(be,2) = centermod(gel(be,2), ell);
     308        1561 :     be = nffactorback(bnfz, be, NULL);
     309             :   }
     310        1603 :   be = reducebeta(bnfz, be, ell);
     311        1603 :   if (DEBUGLEVEL>1) err_printf("beta reduced = %Ps\n",be);
     312        1603 :   return be;
     313             : }
     314             : 
     315             : GEN
     316       56518 : lift_if_rational(GEN x)
     317             : {
     318             :   long lx, i;
     319             :   GEN y;
     320             : 
     321       56518 :   switch(typ(x))
     322             :   {
     323        8771 :     default: break;
     324             : 
     325       32550 :     case t_POLMOD:
     326       32550 :       y = gel(x,2);
     327       32550 :       if (typ(y) == t_POL)
     328             :       {
     329       12068 :         long d = degpol(y);
     330       12068 :         if (d > 0) return x;
     331        2247 :         return (d < 0)? gen_0: gel(y,2);
     332             :       }
     333       20482 :       return y;
     334             : 
     335        6923 :     case t_POL: lx = lg(x);
     336       29792 :       for (i=2; i<lx; i++) gel(x,i) = lift_if_rational(gel(x,i));
     337        6923 :       break;
     338        8274 :     case t_VEC: case t_COL: case t_MAT: lx = lg(x);
     339       36029 :       for (i=1; i<lx; i++) gel(x,i) = lift_if_rational(gel(x,i));
     340             :   }
     341       23968 :   return x;
     342             : }
     343             : 
     344             : /* lift elt t in nf to nfz, algebraic form */
     345             : static GEN
     346         665 : lifttoKz(GEN nf, GEN t, compo_s *C)
     347             : {
     348         665 :   GEN x = nf_to_scalar_or_alg(nf, t);
     349         665 :   if (typ(x) != t_POL) return x;
     350         665 :   return RgX_RgXQ_eval(x, C->p, C->R);
     351             : }
     352             : /* lift ideal id in nf to nfz */
     353             : static GEN
     354         672 : ideallifttoKz(GEN nfz, GEN nf, GEN id, compo_s *C)
     355             : {
     356         672 :   GEN I = idealtwoelt(nf,id);
     357         672 :   GEN x = nf_to_scalar_or_alg(nf, gel(I,2));
     358         672 :   if (typ(x) != t_POL) return gel(I,1);
     359         469 :   gel(I,2) = algtobasis(nfz, RgX_RgXQ_eval(x, C->p, C->R));
     360         469 :   return idealhnf_two(nfz,I);
     361             : }
     362             : 
     363             : static GEN
     364         693 : prlifttoKz_i(GEN nfz, GEN nf, GEN pr, compo_s *C)
     365             : {
     366         693 :   GEN p = pr_get_p(pr), T = nf_get_pol(nfz);
     367         693 :   if (nf_get_degree(nf) != 1)
     368             :   { /* restrict to primes above pr */
     369         665 :     GEN t = pr_get_gen(pr);
     370         665 :     t = Q_primpart( lifttoKz(nf,t,C) );
     371         665 :     T = FpX_gcd(FpX_red(T,p), FpX_red(t,p), p);
     372         665 :     T = FpX_normalize(T, p);
     373             :   }
     374         693 :   return gel(FpX_factor(T, p), 1);
     375             : }
     376             : /* lift ideal pr in nf to ONE prime in nfz (the others are conjugate under tau
     377             :  * and bring no further information on e_1 W). Assume pr coprime to
     378             :  * index of both nf and nfz, and unramified in Kz/K (minor simplification) */
     379             : static GEN
     380         245 : prlifttoKz(GEN nfz, GEN nf, GEN pr, compo_s *C)
     381             : {
     382         245 :   GEN P = prlifttoKz_i(nfz, nf, pr, C);
     383         245 :   return idealprimedec_kummer(nfz, gel(P,1), pr_get_e(pr), pr_get_p(pr));
     384             : }
     385             : static GEN
     386         448 : prlifttoKzall(GEN nfz, GEN nf, GEN pr, compo_s *C)
     387             : {
     388         448 :   GEN P = prlifttoKz_i(nfz, nf, pr, C), p = pr_get_p(pr), vP;
     389         448 :   long l = lg(P), e = pr_get_e(pr), i;
     390         448 :   vP = cgetg(l, t_VEC);
     391        1540 :   for (i = 1; i < l; i++)
     392        1092 :     gel(vP,i) = idealprimedec_kummer(nfz,gel(P,i), e, p);
     393         448 :   return vP;
     394             : }
     395             : 
     396             : static GEN
     397        2275 : get_badbnf(GEN bnf)
     398             : {
     399             :   long i, l;
     400        2275 :   GEN bad = gen_1, gen = bnf_get_gen(bnf);
     401        2275 :   l = lg(gen);
     402        4403 :   for (i = 1; i < l; i++)
     403             :   {
     404        2128 :     GEN g = gel(gen,i);
     405        2128 :     bad = lcmii(bad, gcoeff(g,1,1));
     406             :   }
     407        2275 :   return bad;
     408             : }
     409             : /* test whether H has index p */
     410             : static int
     411        3444 : H_is_good(GEN H, GEN p)
     412             : {
     413        3444 :   long l = lg(H), status = 0, i;
     414        8547 :   for (i = 1; i < l; i++)
     415        6825 :     if (equalii(gcoeff(H,i,i), p) && ++status > 1) return 0;
     416        1722 :   return status == 1;
     417             : }
     418             : static GEN
     419        1064 : bid_primes(GEN bid)
     420             : {
     421        1064 :   GEN v = leafcopy(gel(bid_get_fact(bid),1));
     422        1064 :   long i, l = lg(v);
     423        1911 :   for (i = 1; i < l; i++) gel(v,i) = pr_get_p(gel(v,i));
     424        1064 :   settyp(v, t_VEC); return v;
     425             : }
     426             : /* Let K base field, L/K described by bnr (conductor F) + H. Return a list of
     427             :  * primes coprime to f*ell of degree 1 in K whose images in Cl_f(K) together
     428             :  * with ell*Cl_f(K), generate H:
     429             :  * thus they all split in Lz/Kz; t in Kz is such that
     430             :  * t^(1/p) generates Lz => t is an ell-th power in k(pr) for all such primes.
     431             :  * Restrict to primes not dividing
     432             :  * - the index of the polynomial defining Kz,
     433             :  * - the modulus,
     434             :  * - ell,
     435             :  * - a generator in bnf.gen or bnfz.gen.
     436             :  * If ell | F and Kz != K, set compute the congruence group Hz over Kz
     437             :  * and set *pfa to the conductor factorization. */
     438             : static GEN
     439        1603 : get_prlist(GEN bnr, GEN H, GEN gell, GEN *pfa, struct rnfkummer *kum)
     440             : {
     441        1603 :   pari_sp av0 = avma;
     442        1603 :   GEN Hz = NULL, bnrz = NULL, cycz = NULL, nfz = NULL;
     443        1603 :   GEN cyc = bnr_get_cyc(bnr), nf = bnr_get_nf(bnr), F = gel(bnr_get_mod(bnr),1);
     444        1603 :   GEN bad, Hsofar, L = cgetg(1, t_VEC);
     445             :   forprime_t T;
     446        1603 :   ulong p, ell = gell[2];
     447        1603 :   int Ldone = 0;
     448             : 
     449        1603 :   bad = get_badbnf(bnr_get_bnf(bnr));
     450        1603 :   if (kum)
     451             :   {
     452         672 :     GEN bnfz = kum->bnfz, ideal = gel(bnr_get_mod(bnr), 1);
     453         672 :     GEN badz = lcmii(get_badbnf(bnfz), nf_get_index(bnf_get_nf(bnfz)));
     454         672 :     bad = lcmii(bad,badz);
     455         672 :     nfz = bnf_get_nf(bnfz);
     456         672 :     ideal = ideallifttoKz(nfz, nf, ideal, &kum->COMPO);
     457         672 :     *pfa = idealfactor_partial(nfz, ideal, bid_primes(bnr_get_bid(bnr)));
     458         672 :     if (dvdiu(idealdown(nf, ideal), ell))
     459             :     { /* ell | N(ideal), need Hz = Ker N: Cl_Kz(bothz) -> Cl_K(ideal)/H
     460             :        * to update conductor */
     461         133 :       bnrz = Buchraymod(bnfz, *pfa, nf_INIT, gell);
     462         133 :       cycz = bnr_get_cyc(bnrz);
     463         133 :       Hz = diagonal_shallow(ZV_snf_gcd(cycz, gell));
     464         133 :       if (H_is_good(Hz, gell))
     465             :       {
     466           0 :         *pfa = gel(bnrconductor_factored(bnrz, Hz), 2);
     467           0 :         gerepileall(av0, 2, &L, pfa); return L;
     468             :       }
     469             :     }
     470             :   }
     471        1603 :   bad = lcmii(gcoeff(F,1,1), bad);
     472        1603 :   cyc = ZV_snf_gcd(cyc, gell);
     473        1603 :   Hsofar = diagonal_shallow(cyc);
     474        1603 :   if (H_is_good(Hsofar, gell))
     475             :   {
     476         910 :     Ldone = 1;
     477         910 :     if (!Hz) { gerepileall(av0, pfa? 2: 1, &L, pfa); return L; }
     478             :   }
     479             :   /* restrict to primes not dividing bad and 1 mod ell */
     480         728 :   u_forprime_arith_init(&T, 2, ULONG_MAX, 1, ell);
     481        9555 :   while ((p = u_forprime_next(&T)))
     482             :   {
     483             :     GEN LP;
     484             :     long i, l;
     485        9555 :     if (!umodiu(bad, p)) continue;
     486        8539 :     LP = idealprimedec_limit_f(nf, utoipos(p), 1);
     487        8539 :     l = lg(LP);
     488       14077 :     for (i = 1; i < l; i++)
     489             :     {
     490        6266 :       pari_sp av = avma;
     491        6266 :       GEN M, P = gel(LP,i), v = bnrisprincipalmod(bnr, P, gell, 0);
     492        6266 :       if (!hnf_invimage(H, v)) { set_avma(av); continue; }
     493             :       /* P in H */
     494        1505 :       M = ZM_hnfmodid(shallowconcat(Hsofar, v), cyc);
     495        1505 :       if (Hz)
     496             :       { /* N_{Kz/K} P in H hence P in Hz */
     497         448 :         GEN vP = prlifttoKzall(nfz, nf, P, &kum->COMPO);
     498         448 :         long j, lv = lg(vP);
     499        1260 :         for (j = 1; j < lv; j++)
     500             :         {
     501         945 :           v = bnrisprincipalmod(bnrz, gel(vP,j), gell, 0);
     502         945 :           if (!ZV_equal0(v))
     503             :           {
     504         917 :             Hz = ZM_hnfmodid(shallowconcat(Hz,v), cycz);
     505         917 :             if (H_is_good(Hz, gell))
     506             :             {
     507         133 :               *pfa = gel(bnrconductor_factored(bnrz, Hz), 2);
     508         133 :               if (!Ldone) L = vec_append(L, gel(vP,1));
     509         133 :               gerepileall(av0, 2, &L, pfa); return L;
     510             :             }
     511             :           }
     512             :         }
     513         315 :         P = gel(vP,1);
     514             :       }
     515        1057 :       else if (kum) P = prlifttoKz(nfz, nf, P, &kum->COMPO);
     516        1372 :       if (Ldone || ZM_equal(M, Hsofar)) continue;
     517         791 :       L = vec_append(L, P);
     518         791 :       Hsofar = M;
     519         791 :       if (H_is_good(Hsofar, gell))
     520             :       {
     521         679 :         Ldone = 1;
     522         679 :         if (!Hz) { gerepileall(av0, pfa? 2: 1, &L, pfa); return L; }
     523             :       }
     524             :     }
     525             :   }
     526             :   pari_err_BUG("rnfkummer [get_prlist]"); return NULL;/*LCOV_EXCL_LINE*/
     527             : }
     528             : /*Lprz list of prime ideals in Kz that must split completely in Lz/Kz, vecWA
     529             :  * generators for the S-units used to build the Kummer generators. Return
     530             :  * matsmall M such that \prod WA[j]^x[j] ell-th power mod pr[i] iff
     531             :  * \sum M[i,j] x[j] = 0 (mod ell) */
     532             : static GEN
     533        1603 : subgroup_info(GEN bnfz, GEN Lprz, GEN gell, GEN vecWA)
     534             : {
     535        1603 :   GEN M, nfz = bnf_get_nf(bnfz), Lell = mkvec(gell);
     536        1603 :   long i, j, ell = gell[2], l = lg(vecWA), lz = lg(Lprz);
     537        1603 :   M = cgetg(l, t_MAT);
     538        7035 :   for (j=1; j<l; j++) gel(M,j) = cgetg(lz, t_VECSMALL);
     539        2408 :   for (i=1; i < lz; i++)
     540             :   {
     541         805 :     GEN pr = gel(Lprz,i), EX = subiu(pr_norm(pr), 1);
     542         805 :     GEN N, g,T,p, prM = idealhnf(nfz, pr);
     543         805 :     GEN modpr = zk_to_Fq_init(nfz, &pr,&T,&p);
     544         805 :     long v = Z_lvalrem(divis(EX,ell), ell, &N) + 1; /* Norm(pr)-1 = N * ell^v */
     545         805 :     GEN ellv = powuu(ell, v);
     546         805 :     g = gener_Fq_local(T,p, Lell);
     547         805 :     g = Fq_pow(g,N, T,p); /* order ell^v */
     548        4186 :     for (j=1; j < l; j++)
     549             :     {
     550        3381 :       GEN logc, c = gel(vecWA,j);
     551        3381 :       if (typ(c) == t_MAT) /* famat */
     552        2359 :         c = famat_makecoprime(nfz, gel(c,1), gel(c,2), pr, prM, EX);
     553        3381 :       c = nf_to_Fq(nfz, c, modpr);
     554        3381 :       c = Fq_pow(c, N, T,p);
     555        3381 :       logc = Fq_log(c, g, ellv, T,p);
     556        3381 :       ucoeff(M, i,j) = umodiu(logc, ell);
     557             :     }
     558             :   }
     559        1603 :   return M;
     560             : }
     561             : 
     562             : static GEN
     563        1484 : futu(GEN bnf)
     564             : {
     565        1484 :   GEN fu, tu, SUnits = bnf_get_sunits(bnf);
     566        1484 :   if (SUnits)
     567             :   {
     568         798 :     tu = nf_to_scalar_or_basis(bnf_get_nf(bnf), bnf_get_tuU(bnf));
     569         798 :     fu = bnf_compactfu(bnf);
     570             :   }
     571             :   else
     572             :   {
     573         686 :     GEN U = bnf_build_units(bnf);
     574         686 :     tu = gel(U,1); fu = vecslice(U, 2, lg(U)-1);
     575             :   }
     576        1484 :   return vec_append(fu, tu);
     577             : }
     578             : static GEN
     579        1484 : bnf_cycgenmod(GEN bnf, long ell, GEN *pselmer, long *prc)
     580             : {
     581        1484 :   GEN gen = bnf_get_gen(bnf), cyc = bnf_get_cyc(bnf), nf = bnf_get_nf(bnf);
     582        1484 :   GEN B, r = ZV_to_Flv(cyc, ell);
     583        1484 :   long i, rc, l = lg(gen);
     584        1484 :   B = cgetg(l, t_VEC);
     585        2541 :   for (i = 1; i < l && !r[i]; i++);
     586        1484 :   *prc = rc = i-1; /* ell-rank */
     587        2765 :   for (i = 1; i < l; i++)
     588             :   {
     589        1281 :     GEN G, q, c = gel(cyc,i), g = gel(gen,i);
     590        1281 :     if (i > rc && r[i] != 1) c  = muliu(c, Fl_inv(r[i], ell));
     591        1281 :     q = divis(c, ell); /* remainder = 0 (i <= rc) or 1 */
     592             :     /* compute (b) = g^c mod ell-th powers */
     593        1281 :     G = equali1(q)? g: idealpowred(nf, g, q); /* lose principal part */
     594        1281 :     G = idealpows(nf, G, ell);
     595        1281 :     if (i > rc) G = idealmul(nf, G, g);
     596        1281 :     gel(B,i) = gel(bnfisprincipal0(bnf, G, nf_GENMAT|nf_FORCE), 2);
     597             :   }
     598        1484 :   *pselmer = shallowconcat(futu(bnf), vecslice(B,1,rc));
     599        1484 :   return B;
     600             : }
     601             : 
     602             : static GEN
     603         931 : rnfkummersimple(GEN bnr, GEN H, long ell)
     604             : {
     605             :   long j, lSp, rc;
     606             :   GEN bnf, nf,bid, cycgenmod, Sp, vsprk, matP;
     607         931 :   GEN be, M, K, vecW, vecWB, vecBp, gell = utoipos(ell);
     608             :   /* primes landing in H must be totally split */
     609         931 :   GEN Lpr = get_prlist(bnr, H, gell, NULL, NULL);
     610             : 
     611         931 :   bnf = bnr_get_bnf(bnr); if (!bnf_get_sunits(bnf)) bnf_build_units(bnf);
     612         931 :   nf  = bnf_get_nf(bnf);
     613         931 :   bid = bnr_get_bid(bnr);
     614         931 :   list_Hecke(&Sp, &vsprk, nf, bid_get_fact2(bid), gell, NULL);
     615             : 
     616         931 :   cycgenmod = bnf_cycgenmod(bnf, ell, &vecW, &rc);
     617         931 :   lSp = lg(Sp);
     618         931 :   vecBp = cgetg(lSp, t_VEC);
     619         931 :   matP  = cgetg(lSp, t_MAT);
     620        1750 :   for (j = 1; j < lSp; j++)
     621         819 :     isprincipalell(bnf,gel(Sp,j), cycgenmod,ell,rc, &gel(matP,j),&gel(vecBp,j));
     622         931 :   vecWB = shallowconcat(vecW, vecBp);
     623             : 
     624         931 :   M = matlogall(nf, vecWB, 0, 0, gell, vsprk);
     625         931 :   M = vconcat(M, shallowconcat(zero_Flm(rc,lg(vecW)-1), matP));
     626         931 :   M = vconcat(M, subgroup_info(bnf, Lpr, gell, vecWB));
     627         931 :   K = Flm_ker(M, ell);
     628         931 :   if (ell == 2)
     629             :   {
     630         924 :     GEN msign = nfsign(nf, vecWB), y;
     631         924 :     GEN arch = ZV_to_zv(bid_get_arch(bid)); /* the conductor */
     632         924 :     msign = Flm_mul(msign, K, 2);
     633         924 :     y = Flm_ker(msign, 2);
     634         924 :     y = zv_equal0(arch)? gel(y,1): Flm_Flc_invimage(msign, arch, 2);
     635         924 :     K = Flm_Flc_mul(K, y, 2);
     636             :   }
     637             :   else
     638           7 :     K = gel(K,1);
     639         931 :   be = compute_beta(K, vecWB, gell, bnf);
     640         931 :   be = nf_to_scalar_or_alg(nf, be);
     641         931 :   if (typ(be) == t_POL) be = mkpolmod(be, nf_get_pol(nf));
     642         931 :   return gsub(pol_xn(ell, 0), be);
     643             : }
     644             : 
     645             : static ulong
     646       18844 : nf_to_logFl(GEN nf, GEN x, GEN modpr, ulong g, ulong q, ulong ell, ulong p)
     647             : {
     648       18844 :   x = nf_to_Fp_coprime(nf, x, modpr);
     649       18844 :   return Fl_log(Fl_powu(umodiu(x, p), q, p), g, ell, p);
     650             : }
     651             : static GEN
     652        5635 : nfV_to_logFlv(GEN nf, GEN x, GEN modpr, ulong g, ulong q, ulong ell, ulong p)
     653       24479 : { pari_APPLY_long(nf_to_logFl(nf, gel(x,i), modpr, g, q, ell, p)); }
     654             : 
     655             : /* Compute e_1 Cl(K)/Cl(K)^ell. If u = w^ell a virtual unit, compute
     656             :  * discrete log mod ell on units.gen + bnf.gen (efficient variant of algo
     657             :  * 5.3.11) by finding ru degree 1 primes Pj coprime to everything, and gj
     658             :  * in k(Pj)^* of order ell such that
     659             :  *      log_gj(u_i^((Pj.p - 1) / ell) mod Pj), j = 1..ru
     660             :  * has maximal F_ell rank ru then solve linear system */
     661             : static GEN
     662         553 : kervirtualunit(struct rnfkummer *kum, GEN vselmer)
     663             : {
     664         553 :   GEN bnf = kum->bnfz, cyc = bnf_get_cyc(bnf), nf = bnf_get_nf(bnf);
     665         553 :   GEN W, B, vy, vz, M, U1, U2, vtau, vell, SUnits = bnf_get_sunits(bnf);
     666         553 :   long i, j, r, l = lg(vselmer), rc = kum->rc, ru = l-1 - rc, ell = kum->ell;
     667         553 :   long LIMC = SUnits? itou(gel(SUnits,4)): 1;
     668             :   ulong p;
     669             :   forprime_t T;
     670             : 
     671         553 :   vtau = cgetg(l, t_VEC);
     672         553 :   vell = cgetg(l, t_VEC);
     673        2534 :   for (j = 1; j < l; j++)
     674             :   {
     675        1981 :     GEN t = gel(vselmer,j);
     676        1981 :     if (typ(t) == t_MAT)
     677             :     {
     678             :       GEN ct;
     679        1428 :       t = nffactorback(bnf, gel(t,1), ZV_to_Flv(gel(t,2), ell));
     680        1428 :       t = Q_primitive_part(t, &ct);
     681        1428 :       if (ct)
     682             :       {
     683         686 :         GEN F = Q_factor(ct);
     684         686 :         ct = factorback2(gel(F,1), ZV_to_Flv(gel(F,2), ell));
     685         686 :         t = (typ(t) == t_INT)? ct: ZC_Z_mul(t, ct);
     686             :       }
     687             :     }
     688        1981 :     gel(vell,j) = t; /* integral, not too far from primitive */
     689        1981 :     gel(vtau,j) = tauofelt(t, &kum->tau);
     690             :   }
     691         553 :   U1 = vecslice(vell, 1, ru); /* units */
     692         553 :   U2 = vecslice(vell, ru+1, ru+rc); /* cycgen (mod ell-th powers) */
     693         553 :   B = nf_get_index(nf); /* bad primes; from 1 to ru are LIMC-units */
     694        1008 :   for (i = 1; i <= rc; i++) B = mulii(B, nfnorm(nf, gel(U2,i)));
     695         553 :   if (LIMC > 1)
     696             :   {
     697         553 :     GEN F = absZ_factor_limit(B, LIMC), P = gel(F,1);
     698         553 :     long lP = lg(P);
     699         553 :     B = (lP > 1)? gel(P,lP-1): gen_1;
     700             :   }
     701         553 :   vy = cgetg(l, t_MAT);
     702        2079 :   for (j = 1; j <= ru; j++) gel(vy,j) = zero_Flv(rc); /* units */
     703        1008 :   for (     ; j < l; j++)
     704             :   {
     705         455 :     GEN y, w, u = gel(vtau, j); /* virtual unit */
     706         455 :     if (!idealispower(nf, u, ell, &w)) pari_err_BUG("kervirtualunit");
     707         455 :     y = isprincipal(bnf, w); setlg(y, rc+1);
     708         455 :     if (!ZV_equal0(y))
     709        1330 :       for (i = 1; i <= rc; i++)
     710         875 :         gel(y,i) = diviiexact(mului(ell,gel(y,i)), gel(cyc,i));
     711         455 :     gel(vy,j) = ZV_to_Flv(y, ell);
     712             :   }
     713         553 :   u_forprime_arith_init(&T, LIMC+1, ULONG_MAX, 1, ell);
     714         553 :   M = cgetg(ru+1, t_MAT); r = 1; setlg(M,2);
     715         553 :   vz = cgetg(ru+1, t_MAT);
     716        2219 :   while ((p = u_forprime_next(&T))) if (umodiu(B,p))
     717             :   {
     718        2177 :     GEN P = idealprimedec_limit_f(nf, utoipos(p), 1);
     719        2177 :     long nP = lg(P)-1;
     720        2177 :     ulong g = rootsof1_Fl(ell, p), q = p / ell; /* (p-1) / ell */
     721        4207 :     for (i = 1; i <= nP; i++)
     722             :     {
     723        2583 :       GEN modpr = zkmodprinit(nf, gel(P,i));
     724             :       GEN z, v2;
     725        2583 :       gel(M, r) = nfV_to_logFlv(nf, U1, modpr, g, q, ell, p); /* log futu */
     726        2583 :       if (Flm_rank(M, ell) < r) continue; /* discard */
     727             : 
     728        1526 :       v2 = nfV_to_logFlv(nf, U2, modpr, g, q, ell, p); /* log alpha[1..rc] */
     729        1526 :       gel(vz, r) = z = nfV_to_logFlv(nf, vtau, modpr, g, q, ell, p);
     730        2926 :       for (j = ru+1; j < l; j++)
     731        1400 :         uel(z,j) = Fl_sub(uel(z,j), Flv_dotproduct(v2, gel(vy,j), ell), ell);
     732        1526 :       if (r == ru) break;
     733         973 :       r++; setlg(M, r+1);
     734             :     }
     735        2177 :     if (i < nP) break;
     736             :   }
     737         553 :   if (r != ru) pari_err_BUG("kervirtualunit");
     738             :   /* Solve prod_k U[k]^x[j,k] = vtau[j] / prod_i alpha[i]^vy[j,i] mod (K^*)^ell
     739             :    * for 1 <= j <= #vtau. I.e. for a fixed j: M x[j] = vz[j] (mod ell) */
     740         553 :   M = Flm_inv(Flm_transpose(M), ell);
     741         553 :   vz = Flm_transpose(vz); /* now ru x #vtau */
     742        2534 :   for (j = 1; j < l; j++)
     743        1981 :     gel(vy,j) = shallowconcat(Flm_Flc_mul(M, gel(vz,j), ell), gel(vy,j));
     744         553 :   W = Flm_ker(Flm_Fl_sub(vy, kum->g, ell), ell); l = lg(W);
     745        1666 :   for (j = 1; j < l; j++)
     746        1113 :     gel(W,j) = famat_reduce(famatV_zv_factorback(vselmer, gel(W,j)));
     747         553 :   settyp(W, t_VEC); return W;
     748             : }
     749             : 
     750             : /* - mu_b = sum_{0 <= i < m} floor(r_b r_{m-1-i} / ell) tau^i.
     751             :  * Note that i is in fact restricted to i < m-1 */
     752             : static GEN
     753         840 : get_mmu(long b, GEN r, long ell)
     754             : {
     755         840 :   long i, m = lg(r)-1;
     756         840 :   GEN M = cgetg(m, t_VECSMALL);
     757        2156 :   for (i = 0; i < m-1; i++) M[i+1] = (r[b + 1] * r[m - i]) / ell;
     758         840 :   return M;
     759             : }
     760             : /* max_b zv_sum(mu_b) < m ell */
     761             : static long
     762         567 : max_smu(GEN r, long ell)
     763             : {
     764         567 :   long i, s = 0, z = vecsmall_max(r), l = lg(r);
     765        1288 :   for (i = 2; i < l; i++) s += (z * r[i]) / ell;
     766         567 :   return s;
     767             : }
     768             : 
     769             : /* coeffs(x, a..b) in variable 0 >= varn(x) */
     770             : static GEN
     771        3304 : split_pol(GEN x, long a, long b)
     772             : {
     773        3304 :   long i, l = degpol(x);
     774        3304 :   GEN y = x + a, z;
     775             : 
     776        3304 :   if (l < b) b = l;
     777        3304 :   if (a > b || varn(x) != 0) return pol_0(0);
     778        3304 :   l = b-a + 3;
     779        3304 :   z = cgetg(l, t_POL); z[1] = x[1];
     780       14056 :   for (i = 2; i < l; i++) gel(z,i) = gel(y,i);
     781        3304 :   return normalizepol_lg(z, l);
     782             : }
     783             : 
     784             : /* return (ad * z) mod (T^ell - an/ad), assuming deg_T(z) < 2*ell
     785             :  * allow ad to be NULL (= 1) */
     786             : static GEN
     787        1652 : mod_Xell_a(GEN z, long ell, GEN an, GEN ad, GEN T)
     788             : {
     789        1652 :   GEN z1 = split_pol(z, ell, degpol(z));
     790        1652 :   GEN z0 = split_pol(z, 0,   ell-1); /* z = v^ell z1 + z0*/
     791        1652 :   if (ad) z0 = ZXX_Z_mul(z0, ad);
     792        1652 :   return gadd(z0, ZXQX_ZXQ_mul(z1, an, T));
     793             : }
     794             : /* D*basistoalg(nfz, c), in variable v. Result is integral */
     795             : static GEN
     796        1512 : to_alg(GEN nfz, GEN c, GEN D)
     797             : {
     798        1512 :   if (typ(c) != t_COL) return D? mulii(D,c): c;
     799        1512 :   return RgV_dotproduct(nf_get_zkprimpart(nfz), c);
     800             : }
     801             : /* assume x in alg form */
     802             : static GEN
     803        1652 : downtoK(toK_s *T, GEN x)
     804             : {
     805        1652 :   if (typ(x) != t_POL) return x;
     806        1652 :   x = RgM_RgC_mul(T->invexpoteta1, RgX_to_RgC(x, lg(T->invexpoteta1) - 1));
     807        1652 :   return mkpolmod(RgV_to_RgX(x, varn(T->polnf)), T->polnf);
     808             : }
     809             : 
     810             : /* th. 5.3.5. and prop. 5.3.9. */
     811             : static GEN
     812         672 : compute_polrel(struct rnfkummer *kum, GEN be)
     813             : {
     814         672 :   toK_s *T = &kum->T;
     815         672 :   long i, k, MU = 0, ell = kum->ell, m = T->m;
     816         672 :   GEN r = Fl_powers(kum->g, m-1, ell); /* r[i+1] = g^i mod ell */
     817             :   GEN D, S, root, numa, powtau_Ninvbe, Ninvbe, Dinvbe;
     818         672 :   GEN C, prim_Rk, C_Rk, prim_root, C_root, mell = utoineg(ell);
     819         672 :   GEN nfz = bnf_get_nf(kum->bnfz), Tz = nf_get_pol(nfz), Dz = nf_get_zkden(nfz);
     820             :   pari_timer ti;
     821             : 
     822         672 :   if (DEBUGLEVEL>1) { err_printf("Computing Newton sums: "); timer_start(&ti); }
     823         672 :   if (equali1(Dz)) Dz = NULL;
     824         672 :   D = Dz;
     825         672 :   Ninvbe = Q_remove_denom(nfinv(nfz, be), &Dinvbe);
     826         672 :   powtau_Ninvbe = powtau(Ninvbe, m-1, T->tau);
     827         672 :   if (Dinvbe)
     828             :   {
     829         567 :     MU = max_smu(r, ell);
     830         567 :     D = mul_denom(Dz, powiu(Dinvbe, MU));
     831             :   }
     832             : 
     833         672 :   root = cgetg(ell + 2, t_POL); /* compute D*root, will correct at the end */
     834         672 :   root[1] = evalsigne(1) | evalvarn(0);
     835         672 :   gel(root,2) = gen_0;
     836         672 :   gel(root,3) = D? D: gen_1;
     837        1652 :   for (i = 2; i < ell; i++) gel(root,2+i) = gen_0;
     838        1512 :   for (i = 1; i < m; i++)
     839             :   { /* compute (1/be) ^ (-mu) instead of be^mu [mu < 0].
     840             :      * 1/be = Ninvbe / Dinvbe */
     841         840 :     GEN mmu = get_mmu(i, r, ell), t;
     842         840 :     t = to_alg(nfz, nffactorback(nfz, powtau_Ninvbe, mmu), Dz);/* Ninvbe^-mu */
     843         840 :     if (Dinvbe)
     844             :     {
     845         721 :       long a = MU - zv_sum(mmu);
     846         721 :       if (a) t = gmul(t, powiu(Dinvbe, a));
     847             :     }
     848         840 :     gel(root, 2 + r[i+1]) = t; /* root += D * (z_ell*T)^{r_i} be^mu_i */
     849             :   }
     850             :   /* Other roots are as above with z_ell -> z_ell^j.
     851             :    * Treat all contents (C_*) and principal parts (prim_*) separately */
     852         672 :   prim_root = Q_primitive_part(root, &C_root);
     853         672 :   C_root = div_content(C_root, D);
     854             : 
     855             :   /* theta^ell = be^( sum tau^a r_{d-1-a} ) = a = numa / Dz */
     856         672 :   numa = to_alg(nfz, nffactorback(nfz, powtau(be, m, T->tau),
     857             :                                   vecsmall_reverse(r)), Dz);
     858         672 :   if (DEBUGLEVEL>1) err_printf("root(%ld) ", timer_delay(&ti));
     859             : 
     860             :   /* Compute mod (X^ell - t, nfz.pol) */
     861         672 :   C_Rk = C_root; prim_Rk = prim_root;
     862         672 :   C = div_content(C_root, Dz);
     863         672 :   S = cgetg(ell+3, t_POL); /* Newton sums */
     864         672 :   S[1] = evalsigne(1) | evalvarn(0);
     865         672 :   gel(S,2) = gen_0;
     866        2324 :   for (k = 2; k <= ell; k++)
     867             :   { /* compute the k-th Newton sum; here C_Rk ~ C_root  */
     868        1652 :     pari_sp av = avma;
     869        1652 :     GEN z, C_z, d, Rk = ZXQX_mul(prim_Rk, prim_root, Tz);
     870        1652 :     Rk = mod_Xell_a(Rk, ell, numa, Dz, Tz); /* (mod X^ell - a, nfz.pol) */
     871        1652 :     prim_Rk = Q_primitive_part(Rk, &d); /* d C_root ~ 1 */
     872        1652 :     C_Rk = mul_content(C_Rk, mul_content(d, C));
     873             :     /* root^k = prim_Rk * C_Rk */
     874        1652 :     z = Q_primitive_part(gel(prim_Rk,2), &C_z); /* C_z ~ 1/C_root ~ 1/C_Rk */
     875        1652 :     z = downtoK(T, z);
     876        1652 :     C_z = mul_content(mul_content(C_z, C_Rk), mell);
     877        1652 :     z = gmul(z, C_z); /* C_z ~ 1 */
     878        1652 :     gerepileall(av, C_Rk? 3: 2, &z, &prim_Rk, &C_Rk);
     879        1652 :     if (DEBUGLEVEL>1) err_printf("%ld(%ld) ", k, timer_delay(&ti));
     880        1652 :     gel(S,k+1) = z; /* - Newton sum */
     881             :   }
     882         672 :   gel(S,ell+2) = gen_m1; if (DEBUGLEVEL>1) err_printf("\n");
     883         672 :   return RgX_recip(RgXn_expint(S,ell+1));
     884             : }
     885             : 
     886             : static void
     887         553 : compositum_red(compo_s *C, GEN P, GEN Q)
     888             : {
     889         553 :   GEN p, q, a, z = gel(compositum2(P, Q),1);
     890         553 :   a = gel(z,1);
     891         553 :   p = gel(gel(z,2), 2);
     892         553 :   q = gel(gel(z,3), 2);
     893         553 :   C->k = itos( gel(z,4) );
     894             :   /* reduce R. FIXME: should be polredbest(a, 1), but breaks rnfkummer bench */
     895         553 :   z = polredabs0(a, nf_ORIG|nf_PARTIALFACT);
     896         553 :   C->R = gel(z,1);
     897         553 :   a = gel(gel(z,2), 2);
     898         553 :   C->p = RgX_RgXQ_eval(p, a, C->R);
     899         553 :   C->q = RgX_RgXQ_eval(q, a, C->R);
     900         553 :   C->rev = QXQ_reverse(a, C->R);
     901         553 : }
     902             : 
     903             : /* replace P->C^(-deg P) P(xC) for the largest integer C such that coefficients
     904             :  * remain algebraic integers. Lift *rational* coefficients */
     905             : static void
     906         672 : nfX_Z_normalize(GEN nf, GEN P)
     907             : {
     908             :   long i, l;
     909         672 :   GEN C, Cj, PZ = cgetg_copy(P, &l);
     910         672 :   PZ[1] = P[1];
     911        3668 :   for (i = 2; i < l; i++) /* minor variation on RgX_to_nfX (create PZ) */
     912             :   {
     913        2996 :     GEN z = nf_to_scalar_or_basis(nf, gel(P,i));
     914        2996 :     if (typ(z) == t_INT)
     915        2043 :       gel(PZ,i) = gel(P,i) = z;
     916             :     else
     917         953 :       gel(PZ,i) = ZV_content(z);
     918             :   }
     919         672 :   (void)ZX_Z_normalize(PZ, &C);
     920             : 
     921         672 :   if (C == gen_1) return;
     922         200 :   Cj = C;
     923         828 :   for (i = l-2; i > 1; i--)
     924             :   {
     925         628 :     if (i != l-2) Cj = mulii(Cj, C);
     926         628 :     gel(P,i) = gdiv(gel(P,i), Cj);
     927             :   }
     928             : }
     929             : 
     930             : /* set kum->vecC, kum->tQ */
     931             : static void
     932         315 : _rnfkummer_step4(struct rnfkummer *kum, long d, long m)
     933             : {
     934         315 :   long i, j, rc = kum->rc;
     935         315 :   GEN Q, vT, vB, vC, vz, B = cgetg(rc+1,t_VEC), T = cgetg(rc+1,t_MAT);
     936         315 :   GEN gen = bnf_get_gen(kum->bnfz), cycgenmod = kum->cycgenmod;
     937         315 :   ulong ell = kum->ell;
     938             : 
     939         770 :   for (j = 1; j <= rc; j++)
     940             :   {
     941         455 :     GEN t = tauofideal(gel(gen,j), &kum->tau);
     942         455 :     isprincipalell(kum->bnfz, t, cycgenmod,ell,rc, &gel(T,j), &gel(B,j));
     943             :   }
     944         315 :   Q = Flm_ker(Flm_Fl_sub(Flm_transpose(T), kum->g, ell), ell);
     945         315 :   kum->tQ = lg(Q) == 1? NULL: Flm_transpose(Q);
     946         315 :   kum->vecC = vC = cgetg(rc+1, t_VEC);
     947             :   /* T = rc x rc matrix */
     948         315 :   vT = Flm_powers(T, m-2, ell);
     949         315 :   vB = cgetg(m, t_VEC);
     950         315 :   vz = cgetg(rc+1, t_VEC);
     951         770 :   for (i = 1; i <= rc; i++) gel(vz, i) = cgetg(m, t_VEC);
     952         714 :   for (j = 1; j < m; j++)
     953             :   {
     954         399 :     GEN Tj = Flm_Fl_mul(gel(vT,m-j), Fl_mul(j,d,ell), ell);
     955         399 :     gel(vB, j) = tauofvec(j == 1? B: gel(vB, j-1), &kum->tau);
     956         980 :     for (i = 1; i <= rc; i++) gmael(vz, i, j) = gel(Tj, i);
     957             :   }
     958         315 :   vB = shallowconcat1(vB);
     959         770 :   for (i = 1; i <= rc; i++)
     960             :   {
     961         455 :     GEN z = shallowconcat1(gel(vz,i));
     962         455 :     gel(vC,i) = famat_reduce(famatV_zv_factorback(vB, z));
     963             :   }
     964         315 : }
     965             : 
     966             : /* alg 5.3.5 */
     967             : static void
     968         553 : rnfkummer_init(struct rnfkummer *kum, GEN bnf, ulong ell, long prec)
     969             : {
     970         553 :   compo_s *COMPO = &kum->COMPO;
     971         553 :   toK_s *T = &kum->T;
     972         553 :   GEN nf  = bnf_get_nf(bnf), polnf = nf_get_pol(nf), vselmer, bnfz;
     973             :   long degK, degKz, m, d;
     974             :   ulong g;
     975             :   pari_timer ti;
     976         553 :   if (DEBUGLEVEL>2) err_printf("Step 1\n");
     977         553 :   if (DEBUGLEVEL) timer_start(&ti);
     978         553 :   compositum_red(COMPO, polnf, polcyclo(ell, varn(polnf)));
     979         553 :   if (DEBUGLEVEL)
     980             :   {
     981           0 :     timer_printf(&ti, "[rnfkummer] compositum");
     982           0 :     if (DEBUGLEVEL>1) err_printf("polred(compositum) = %Ps\n",COMPO->R);
     983             :   }
     984         553 :   if (DEBUGLEVEL>2) err_printf("Step 2\n");
     985         553 :   degK  = degpol(polnf);
     986         553 :   degKz = degpol(COMPO->R);
     987         553 :   m = degKz / degK; /* > 1 */
     988         553 :   d = (ell-1) / m;
     989         553 :   g = Fl_powu(pgener_Fl(ell), d, ell);
     990         553 :   if (Fl_powu(g, m, ell*ell) == 1) g += ell;
     991             :   /* ord(g) = m in all (Z/ell^k)^* */
     992         553 :   if (DEBUGLEVEL>2) err_printf("Step 3\n");
     993             :   /* could factor disc(R) using th. 2.1.6. */
     994         553 :   kum->bnfz = bnfz = Buchall(COMPO->R, nf_FORCE, maxss(prec,BIGDEFAULTPREC));
     995         553 :   if (DEBUGLEVEL) timer_printf(&ti, "[rnfkummer] bnfinit(Kz)");
     996         553 :   kum->cycgenmod = bnf_cycgenmod(bnfz, ell, &vselmer, &kum->rc);
     997         553 :   kum->ell = ell;
     998         553 :   kum->g = g;
     999         553 :   kum->mgi = Fl_div(m, g, ell);
    1000         553 :   get_tau(kum);
    1001         553 :   if (DEBUGLEVEL>2) err_printf("Step 4\n");
    1002         553 :   if (kum->rc)
    1003         315 :     _rnfkummer_step4(kum, d, m);
    1004             :   else
    1005         238 :   { kum->vecC = cgetg(1, t_VEC); kum->tQ = NULL; }
    1006         553 :   if (DEBUGLEVEL>2) err_printf("Step 5\n");
    1007         553 :   kum->vecW = kervirtualunit(kum, vselmer);
    1008         553 :   if (DEBUGLEVEL>2) err_printf("Step 8\n");
    1009             :   /* left inverse */
    1010         553 :   T->invexpoteta1 = QM_inv(RgXQ_matrix_pow(COMPO->p, degKz, degK, COMPO->R));
    1011         553 :   T->polnf = polnf;
    1012         553 :   T->tau = &kum->tau;
    1013         553 :   T->m = m;
    1014         553 :   T->powg = Fl_powers(g, m, ell);
    1015         553 : }
    1016             : 
    1017             : static GEN
    1018         672 : rnfkummer_ell(struct rnfkummer *kum, GEN bnr, GEN H)
    1019             : {
    1020         672 :   ulong ell = kum->ell;
    1021         672 :   GEN bnfz = kum->bnfz, nfz = bnf_get_nf(bnfz), gell = utoipos(ell);
    1022         672 :   GEN vecC = kum->vecC, vecW = kum->vecW, cycgenmod = kum->cycgenmod;
    1023         672 :   long lW = lg(vecW), rc = kum->rc, j, lSp;
    1024         672 :   toK_s *T = &kum->T;
    1025             :   GEN K, be, P, faFz, vsprk, Sp, vecAp, vecBp, matP, vecWA, vecWB, M, lambdaWB;
    1026             :   /* primes landing in H must be totally split */
    1027         672 :   GEN Lpr = get_prlist(bnr, H, gell, &faFz, kum);
    1028             : 
    1029         672 :   if (DEBUGLEVEL>2) err_printf("Step 9, 10 and 11\n");
    1030         672 :   list_Hecke(&Sp, &vsprk, nfz, faFz, gell, T->tau);
    1031             : 
    1032         672 :   if (DEBUGLEVEL>2) err_printf("Step 12\n");
    1033         672 :   lSp = lg(Sp);
    1034         672 :   vecAp = cgetg(lSp, t_VEC);
    1035         672 :   vecBp = cgetg(lSp, t_VEC);
    1036         672 :   matP  = cgetg(lSp, t_MAT);
    1037        1015 :   for (j = 1; j < lSp; j++)
    1038             :   {
    1039             :     GEN e, a;
    1040         343 :     isprincipalell(bnfz, gel(Sp,j), cycgenmod,ell,rc, &e, &a);
    1041         343 :     gel(matP,j) = e;
    1042         343 :     gel(vecBp,j) = famat_mul_shallow(famatV_zv_factorback(vecC, zv_neg(e)), a);
    1043         343 :     gel(vecAp,j) = lambdaofelt(gel(vecBp,j), T);
    1044             :   }
    1045         672 :   if (DEBUGLEVEL>2) err_printf("Step 13\n");
    1046         672 :   vecWA = shallowconcat(vecW, vecAp);
    1047         672 :   vecWB = shallowconcat(vecW, vecBp);
    1048             : 
    1049         672 :   if (DEBUGLEVEL>2) err_printf("Step 14, 15 and 17\n");
    1050         672 :   M = matlogall(nfz, vecWA, lW, kum->mgi, gell, vsprk);
    1051         672 :   if (kum->tQ)
    1052             :   {
    1053         245 :     GEN QtP = Flm_mul(kum->tQ, matP, ell);
    1054         245 :     M = vconcat(M, shallowconcat(zero_Flm(lgcols(kum->tQ)-1,lW-1), QtP));
    1055             :   }
    1056         672 :   lambdaWB = shallowconcat(lambdaofvec(vecW, T), vecAp);/*vecWB^lambda*/
    1057         672 :   M = vconcat(M, subgroup_info(bnfz, Lpr, gell, lambdaWB));
    1058         672 :   if (DEBUGLEVEL>2) err_printf("Step 16\n");
    1059         672 :   K = Flm_ker(M, ell);
    1060         672 :   if (DEBUGLEVEL>2) err_printf("Step 18\n");
    1061         672 :   be = compute_beta(gel(K,1), vecWB, gell, kum->bnfz);
    1062         672 :   P = compute_polrel(kum, be);
    1063         672 :   nfX_Z_normalize(bnr_get_nf(bnr), P);
    1064         672 :   if (DEBUGLEVEL>1) err_printf("polrel(beta) = %Ps\n", P);
    1065         672 :   return P;
    1066             : }
    1067             : 
    1068             : static void
    1069        1302 : bnrclassfield_sanitize(GEN *pbnr, GEN *pH)
    1070             : {
    1071             :   GEN T;
    1072        1302 :   bnr_subgroup_sanitize(pbnr, pH);
    1073        1267 :   T = nf_get_pol(bnr_get_nf(*pbnr));
    1074        1267 :   if (!varn(T)) pari_err_PRIORITY("bnrclassfield", T, "=", 0);
    1075        1253 : }
    1076             : 
    1077             : static GEN
    1078         861 : _rnfkummer(GEN bnr, GEN H, long prec)
    1079             : {
    1080             :   ulong ell;
    1081             :   GEN gell;
    1082             :   struct rnfkummer kum;
    1083             : 
    1084         861 :   bnrclassfield_sanitize(&bnr, &H);
    1085         854 :   gell = H? ZM_det(H): ZV_prod(bnr_get_cyc(bnr));
    1086         854 :   ell = itou(gell);
    1087         854 :   if (ell == 1) return pol_x(0);
    1088         854 :   if (!uisprime(ell)) pari_err_IMPL("rnfkummer for composite relative degree");
    1089         854 :   if (bnf_get_tuN(bnr_get_bnf(bnr)) % ell == 0)
    1090         504 :     return rnfkummersimple(bnr, H, ell);
    1091         350 :   rnfkummer_init(&kum, bnr_get_bnf(bnr), ell, prec);
    1092         350 :   return rnfkummer_ell(&kum, bnr, H);
    1093             : }
    1094             : 
    1095             : GEN
    1096         644 : rnfkummer(GEN bnr, GEN H, long prec)
    1097         644 : { pari_sp av = avma; return gerepilecopy(av, _rnfkummer(bnr, H, prec)); }
    1098             : 
    1099             : /*******************************************************************/
    1100             : /*                        bnrclassfield                            */
    1101             : /*******************************************************************/
    1102             : 
    1103             : /* TODO: could be exported */
    1104             : static void
    1105      199647 : gsetvarn(GEN x, long v)
    1106             : {
    1107             :   long i;
    1108      199647 :   switch(typ(x))
    1109             :   {
    1110        1498 :     case t_POL: case t_SER:
    1111        1498 :       setvarn(x, v); return;
    1112           0 :     case t_LIST:
    1113           0 :       x = list_data(x); if (!x) return;
    1114             :       /* fall through t_VEC */
    1115             :     case t_VEC: case t_COL: case t_MAT:
    1116      222957 :       for (i = lg(x)-1; i > 0; i--) gsetvarn(gel(x,i), v);
    1117             :   }
    1118             : }
    1119             : 
    1120             : /* emb root of pol as polmod modulo pol2, return relative polynomial */
    1121             : static GEN
    1122         196 : relative_pol(GEN pol, GEN emb, GEN pol2)
    1123             : {
    1124             :   GEN eqn, polrel;
    1125         196 :   if (degree(pol)==1) return pol2;
    1126         147 :   emb = liftpol(emb);
    1127         147 :   eqn = gsub(emb, pol_x(varn(pol)));
    1128         147 :   eqn = Q_remove_denom(eqn, NULL);
    1129         147 :   polrel = nfgcd(pol2, eqn, pol, NULL);
    1130         147 :   return RgX_Rg_div(polrel, leading_coeff(polrel));
    1131             : }
    1132             : 
    1133             : /* pol defines K/nf */
    1134             : static GEN
    1135         217 : bnrclassfield_tower(GEN bnr, GEN subgroup, GEN TB, GEN p, long finaldeg, long absolute, long prec)
    1136             : {
    1137         217 :   pari_sp av = avma;
    1138             :   GEN nf, nf2, rnf, bnf, bnf2, bnr2, q, H, dec, cyc, pk, sgpk, pol2, emb, emb2, famod, fa, Lbad, cert;
    1139             :   long i, r1, ell, sp, spk, last;
    1140             :   forprime_t iter;
    1141             : 
    1142         217 :   bnf = bnr_get_bnf(bnr);
    1143         217 :   nf = bnf_get_nf(bnf);
    1144         217 :   rnf = rnfinit0(nf, TB, 1);
    1145         217 :   nf2 = rnf_build_nfabs(rnf, prec);
    1146         217 :   gsetvarn(nf2, varn(nf_get_pol(nf)));
    1147             : 
    1148         217 :   cert = nfcertify(nf2);
    1149         217 :   if (lg(cert)>1)
    1150             :   {
    1151           0 :     rnf = rnfinit0(nf, gel(TB,1), 1);
    1152           0 :     nf2 = rnf_build_nfabs(rnf, prec);
    1153           0 :     gsetvarn(nf2, varn(nf_get_pol(nf)));
    1154             :   }
    1155             : 
    1156         217 :   r1 = nf_get_r1(nf2);
    1157         217 :   bnf2 = Buchall(nf2, nf_FORCE, prec);
    1158             : 
    1159         217 :   sp = itos(p);
    1160         217 :   spk = sp * rnf_get_degree(rnf);
    1161         217 :   pk = stoi(spk);
    1162         217 :   sgpk = hnfmodid(subgroup,pk);
    1163         217 :   last = spk==finaldeg;
    1164             : 
    1165             :   /* compute conductor */
    1166         217 :   famod = gel(bid_get_fact2(bnr_get_bid(bnr)),1);
    1167         217 :   if (lg(famod)==1)
    1168             :   {
    1169         140 :     fa = trivial_fact();
    1170         140 :     Lbad = cgetg(1, t_VECSMALL);
    1171             :   }
    1172             :   else
    1173             :   {
    1174          77 :     long j=1;
    1175          77 :     fa = cgetg(3, t_MAT);
    1176          77 :     gel(fa,1) = cgetg(lg(famod), t_VEC);
    1177          77 :     Lbad = cgetg(lg(famod), t_VEC);
    1178         189 :     for(i=1; i<lg(famod); i++)
    1179             :     {
    1180         112 :       GEN pr = gel(famod,i);
    1181         112 :       gmael(fa,1,i) = rnfidealprimedec(rnf, pr);
    1182         112 :       q = pr_get_p(pr);
    1183         112 :       if (lgefint(q) == 3) gel(Lbad,j++) = q;
    1184             :     }
    1185          77 :     setlg(Lbad,j);
    1186          77 :     Lbad = ZV_to_zv(ZV_sort_uniq(Lbad));
    1187          77 :     gel(fa,1) = shallowconcat1(gel(fa,1));
    1188          77 :     settyp(gel(fa,1), t_COL);
    1189          77 :     gel(fa,2) = cgetg(lg(gel(fa,1)), t_COL);
    1190         203 :     for (i=1; i<lg(gel(fa,1)); i++)
    1191             :     {
    1192         126 :       GEN pr = gcoeff(fa,i,1);
    1193         126 :       long e = equalii(p, pr_get_p(pr))? 1 + (pr_get_e(pr)*sp) / (sp-1): 1;
    1194         126 :       gcoeff(fa,i,2) = utoipos(e);
    1195             :     }
    1196             :   }
    1197         217 :   bnr2 = Buchraymod(bnf2, mkvec2(fa, const_vec(r1,gen_1)), nf_INIT, pk);
    1198             : 
    1199             :   /* compute subgroup */
    1200         217 :   cyc = bnr_get_cyc(bnr2);
    1201         217 :   H = Flm_image(zv_diagonal(ZV_to_Flv(cyc,sp)), sp);
    1202         217 :   u_forprime_init(&iter, 2, ULONG_MAX);
    1203       12754 :   while ((ell = u_forprime_next(&iter))) if (!zv_search(Lbad, ell))
    1204             :   {
    1205       12656 :     dec = idealprimedec_limit_f(nf, utoi(ell), 1);
    1206       24836 :     for (i=1; i<lg(dec); i++)
    1207             :     {
    1208       12180 :       GEN pr = gel(dec,i), Pr = gel(rnfidealprimedec(rnf, pr), 1);
    1209       12180 :       long f = pr_get_f(Pr) / pr_get_f(pr);
    1210       12180 :       GEN vpr = FpC_Fp_mul(bnrisprincipalmod(bnr,pr,pk,0), utoi(f), pk);
    1211       12180 :       if (gequal0(ZC_hnfrem(vpr,sgpk)))
    1212        1260 :         H = vec_append(H, ZV_to_Flv(bnrisprincipalmod(bnr2,Pr,p,0), sp));
    1213             :     }
    1214       12656 :     if (lg(H) > lg(cyc)+3)
    1215             :     {
    1216         217 :       H = Flm_image(H, sp);
    1217         217 :       if (lg(cyc)-lg(H) == 1) break;
    1218             :     }
    1219             :   }
    1220         217 :   H = hnfmodid(shallowconcat(zm_to_ZM(H), diagonal_shallow(cyc)), p);
    1221             : 
    1222             :   /* polynomial over nf2 */
    1223         217 :   pol2 = _rnfkummer(bnr2, H, prec);
    1224             :   /* absolute polynomial */
    1225         217 :   pol2 = rnfequation2(nf2, pol2);
    1226         217 :   emb2 = gel(pol2,2); /* generator of nf2 as polmod modulo pol2 */
    1227         217 :   pol2 = gel(pol2,1);
    1228             :   /* polynomial over nf */
    1229         217 :   if (!absolute || !last)
    1230             :   {
    1231         196 :     emb = rnf_get_alpha(rnf); /* generator of nf as polynomial in nf2 */
    1232         196 :     emb = poleval(emb, emb2); /* generator of nf as polmod modulo pol2 */
    1233         196 :     pol2 = relative_pol(nf_get_pol(nf), emb, pol2);
    1234             :   }
    1235         217 :   if (!last) pol2 = rnfpolredbest(nf, pol2, 0);
    1236             : 
    1237         217 :   obj_free(rnf);
    1238         217 :   pol2 = gerepilecopy(av, pol2);
    1239         217 :   if (last) return pol2;
    1240          56 :   TB = mkvec2(pol2, gel(TB,2));
    1241          56 :   return bnrclassfield_tower(bnr, subgroup, TB, p, finaldeg, absolute, prec);
    1242             : }
    1243             : 
    1244             : /* subgroups H_i of bnr s.t. bnr/H_i is cyclic and inter_i H_i = subgroup */
    1245             : static GEN
    1246         539 : cyclic_compos(GEN subgroup)
    1247             : {
    1248         539 :   pari_sp av = avma;
    1249         539 :   GEN Ui, L, pe, D = ZM_snf_group(subgroup, NULL, &Ui);
    1250         539 :   long i, l = lg(D);
    1251             : 
    1252         539 :   L = cgetg(l, t_VEC);
    1253         539 :   if (l == 1) return L;
    1254         539 :   pe = gel(D,1);
    1255        1288 :   for (i = 1; i < l; i++)
    1256         749 :     gel(L,i) = hnfmodid(shallowconcat(subgroup, vecsplice(Ui,i)),pe);
    1257         539 :   return gerepilecopy(av, L);
    1258             : }
    1259             : 
    1260             : /* p prime; set pkum=NULL if p-th root of unity in base field
    1261             :  * absolute=1 allowed if extension is cyclic with exponent>1 */
    1262             : static GEN
    1263         539 : bnrclassfield_primepower(struct rnfkummer *pkum, GEN bnr, GEN subgroup, GEN p,
    1264             :   GEN P, long absolute, long prec)
    1265             : {
    1266         539 :   GEN res, subs = cyclic_compos(subgroup);
    1267         539 :   long i, l = lg(subs);
    1268             : 
    1269         539 :   res = cgetg(l,t_VEC);
    1270        1288 :   for (i = 1; i < l; i++)
    1271             :   {
    1272         749 :     GEN H = gel(subs,i), cnd = bnrconductormod(bnr, hnfmodid(H,p), p);
    1273         749 :     GEN pol, pe, bnr2 = gel(cnd,2), Hp = gel(cnd,3);
    1274         749 :     if (pkum) pol = rnfkummer_ell(pkum, bnr2, Hp);
    1275         427 :     else      pol = rnfkummersimple(bnr2, Hp, itos(p));
    1276         749 :     pe = ZM_det_triangular(H);
    1277         749 :     if (!equalii(p,pe))
    1278         161 :       pol = bnrclassfield_tower(bnr, H, mkvec2(pol,P), p, itos(pe), absolute, prec);
    1279         749 :     gel(res,i) = pol;
    1280             :   }
    1281         539 :   return res;
    1282             : }
    1283             : 
    1284             : /* partition of v into two subsets whose products are as balanced as possible */
    1285             : /* assume v sorted */
    1286             : static GEN
    1287         105 : vecsmall_balance(GEN v)
    1288             : {
    1289             :   forvec_t T;
    1290         105 :   GEN xbounds, x, vuniq, mult, ind, prod, prodbest = gen_0, bound,
    1291         105 :       xbest = NULL, res1, res2;
    1292         105 :   long i=1, j, k1, k2;
    1293         105 :   if (lg(v) == 3) return mkvec2(mkvecsmall(1), mkvecsmall(2));
    1294          35 :   vuniq = cgetg(lg(v), t_VECSMALL);
    1295          35 :   mult = cgetg(lg(v), t_VECSMALL);
    1296          35 :   ind = cgetg(lg(v), t_VECSMALL);
    1297          35 :   vuniq[1] = v[1];
    1298          35 :   mult[1] = 1;
    1299          35 :   ind[1] = 1;
    1300         140 :   for (j=2; j<lg(v); j++)
    1301             :   {
    1302         105 :     if (v[j] == vuniq[i]) mult[i]++;
    1303             :     else
    1304             :     {
    1305          14 :       i++;
    1306          14 :       vuniq[i] = v[j];
    1307          14 :       mult[i] = 1;
    1308          14 :       ind[i] = j;
    1309             :     }
    1310             :   }
    1311          35 :   setlg(vuniq, ++i);
    1312          35 :   setlg(mult, i);
    1313          35 :   setlg(ind, i);
    1314             : 
    1315          35 :   vuniq = zv_to_ZV(vuniq);
    1316          35 :   prod = factorback2(vuniq, mult);
    1317          35 :   bound = sqrti(prod);
    1318          35 :   xbounds = cgetg(lg(mult), t_VEC);
    1319          84 :   for (i=1; i<lg(mult); i++) gel(xbounds,i) = mkvec2s(0,mult[i]);
    1320             : 
    1321          35 :   forvec_init(&T, xbounds, 0);
    1322         238 :   while ((x = forvec_next(&T)))
    1323             :   {
    1324         203 :     prod = factorback2(vuniq, x);
    1325         203 :     if (cmpii(prod,bound)<=0 && cmpii(prod,prodbest)>0)
    1326             :     {
    1327          91 :       prodbest = prod;
    1328          91 :       xbest = gcopy(x);
    1329             :     }
    1330             :   }
    1331          35 :   res1 = cgetg(lg(v), t_VECSMALL);
    1332          35 :   res2 = cgetg(lg(v), t_VECSMALL);
    1333          84 :   for (i=1,k1=1,k2=1; i<lg(xbest); i++)
    1334             :   {
    1335         105 :     for (j=0; j<itos(gel(xbest,i)); j++) res1[k1++] = ind[i]+j;
    1336         133 :     for (; j<mult[i]; j++)               res2[k2++] = ind[i]+j;
    1337             :   }
    1338          35 :   setlg(res1, k1);
    1339          35 :   setlg(res2, k2); return mkvec2(res1, res2);
    1340             : }
    1341             : 
    1342             : /* TODO nfcompositum should accept vectors of pols */
    1343             : /* assume all fields are linearly disjoint */
    1344             : /* assume the polynomials are sorted by degree */
    1345             : static GEN
    1346         371 : nfcompositumall(GEN nf, GEN L)
    1347             : {
    1348             :   GEN pol, vdeg, part;
    1349             :   long i;
    1350         371 :   if (lg(L)==2) return gel(L,1);
    1351         105 :   vdeg = cgetg(lg(L), t_VECSMALL);
    1352         385 :   for (i=1; i<lg(L); i++) vdeg[i] = degree(gel(L,i));
    1353         105 :   part = vecsmall_balance(vdeg);
    1354         105 :   pol = cgetg(3, t_VEC);
    1355         315 :   for (i = 1; i < 3; i++)
    1356             :   {
    1357         210 :     GEN L2 = vecpermute(L, gel(part,i)), T = nfcompositumall(nf, L2);
    1358         210 :     gel(pol,i) = rnfpolredbest(nf, T, 0);
    1359             :   }
    1360         105 :   return nfcompositum(nf, gel(pol,1), gel(pol,2), 2);
    1361             : }
    1362             : 
    1363             : static GEN
    1364         392 : bad_primes(GEN bnr)
    1365             : {
    1366         392 :   GEN v = bid_primes(bnr_get_bid(bnr));
    1367         392 :   GEN r = nf_get_ramified_primes(bnr_get_nf(bnr));
    1368         392 :   return ZV_sort_uniq(shallowconcat(r, v));
    1369             : }
    1370             : static struct rnfkummer **
    1371         392 : rnfkummer_initall(GEN bnr, GEN vP, long prec)
    1372             : {
    1373         392 :   long i, w, lP = lg(vP);
    1374         392 :   GEN bnf = bnr_get_bnf(bnr);
    1375             :   struct rnfkummer **vkum;
    1376             : 
    1377         392 :   w = bnf_get_tuN(bnf);
    1378         392 :   vkum = (struct rnfkummer **)stack_malloc(vP[lP-1] * sizeof(struct rnfkummer*));
    1379         868 :   for (i = 1; i < lP; i++)
    1380             :   {
    1381         476 :     long p = vP[i];
    1382         476 :     if (w % p == 0) { vkum[p] = NULL; continue; }
    1383         203 :     vkum[p] = (struct rnfkummer *)stack_malloc(sizeof(struct rnfkummer));
    1384         203 :     rnfkummer_init(vkum[p], bnf, p, prec);
    1385             :   }
    1386         392 :   return vkum;
    1387             : }
    1388             : /* fully handle a single congruence subgroup H */
    1389             : static GEN
    1390         455 : bnrclassfield_H(struct rnfkummer **vkum, GEN bnr, GEN bad, GEN H0, GEN fa, long flag,
    1391             :                 long prec)
    1392             : {
    1393         455 :   GEN PN = gel(fa,1), EN = gel(fa,2), res;
    1394         455 :   long i, lPN = lg(PN), absolute;
    1395             : 
    1396         455 :   if (lPN == 1) return pol_x(0);
    1397         455 :   absolute = flag==2 && lPN==2 && !equali1(gel(EN,1)); /* one prime, exponent > 1 */
    1398         455 :   res = cgetg(lPN, t_VEC);
    1399         994 :   for (i = 1; i < lPN; i++)
    1400             :   {
    1401         539 :     GEN p = gel(PN,i), H = hnfmodid(H0, powii(p, gel(EN,i)));
    1402         539 :     long sp = itos(p);
    1403         539 :     if (absolute) absolute = FpM_rank(H,p)==lg(H)-2; /* cyclic */
    1404         539 :     gel(res,i) = bnrclassfield_primepower(vkum[sp], bnr, H, p, bad, absolute, prec);
    1405             :   }
    1406         455 :   res = liftpol_shallow(shallowconcat1(res));
    1407         455 :   res = gen_sort_shallow(res, (void*)cmp_RgX, gen_cmp_RgX);
    1408         455 :   if (flag)
    1409             :   {
    1410         161 :     GEN nf = bnr_get_nf(bnr);
    1411         161 :     res = nfcompositumall(nf, res);
    1412         161 :     if (flag==2 && !absolute) res = rnfequation(nf, res);
    1413             :   }
    1414         455 :   return res;
    1415             : }
    1416             : /* for a vector of subgroups */
    1417             : static GEN
    1418          42 : bnrclassfieldvec(GEN bnr, GEN v, long flag, long prec)
    1419             : {
    1420          42 :   long j, lv = lg(v);
    1421          42 :   GEN vH, vfa, vP, P, w = cgetg(lv, t_VEC);
    1422          42 :   struct rnfkummer **vkum = NULL;
    1423             : 
    1424          42 :   if (lv == 1) return w;
    1425          35 :   vH = cgetg(lv, t_VEC);
    1426          35 :   vP = cgetg(lv, t_VEC);
    1427          35 :   vfa = cgetg(lv, t_VEC);
    1428         126 :   for (j = 1; j < lv; j++)
    1429             :   {
    1430          98 :     GEN N, fa, H = bnr_subgroup_check(bnr, gel(v,j), &N);
    1431          98 :     if (is_bigint(N)) pari_err_OVERFLOW("bnrclassfield [too large degree]");
    1432          91 :     if (!H) H = diagonal_shallow(bnr_get_cyc(bnr));
    1433          91 :     gel(vH,j) = H;
    1434          91 :     gel(vfa,j) = fa = Z_factor(N);
    1435          91 :     gel(vP,j) = ZV_to_zv(gel(fa, 1));
    1436             :   }
    1437          28 :   vP = shallowconcat1(vP); vecsmall_sort(vP);
    1438          28 :   vP = vecsmall_uniq_sorted(vP);
    1439          28 :   if (lg(vP) > 1) vkum = rnfkummer_initall(bnr, vP, prec);
    1440          28 :   P = bad_primes(bnr);
    1441         119 :   for (j = 1; j < lv; j++)
    1442          91 :     gel(w,j) = bnrclassfield_H(vkum, bnr, P, gel(vH,j), gel(vfa,j), flag, prec);
    1443          28 :   return w;
    1444             : }
    1445             : /* flag:
    1446             :  * 0 list of polynomials whose compositum is the extension
    1447             :  * 1 single polynomial
    1448             :  * 2 single absolute polynomial */
    1449             : GEN
    1450         497 : bnrclassfield(GEN bnr, GEN subgroup, long flag, long prec)
    1451             : {
    1452         497 :   pari_sp av = avma;
    1453             :   GEN N, fa, P;
    1454             :   struct rnfkummer **vkum;
    1455             : 
    1456         497 :   if (flag<0 || flag>2) pari_err_FLAG("bnrclassfield [must be 0,1 or 2]");
    1457         483 :   if (subgroup && typ(subgroup) == t_VEC)
    1458             :   {
    1459          49 :     if (nftyp(bnr)==typ_BNF) bnr = Buchray(bnr, gen_1, nf_INIT);
    1460          35 :     else checkbnr(bnr);
    1461          49 :     if (!char_check(bnr_get_cyc(bnr), subgroup))
    1462          42 :       return gerepilecopy(av, bnrclassfieldvec(bnr, subgroup, flag, prec));
    1463             :   }
    1464         441 :   bnrclassfield_sanitize(&bnr, &subgroup);
    1465         399 :   N = ZM_det_triangular(subgroup);
    1466         399 :   if (equali1(N)) { set_avma(av); return pol_x(0); }
    1467         378 :   if (is_bigint(N)) pari_err_OVERFLOW("bnrclassfield [too large degree]");
    1468         364 :   fa = Z_factor(N);
    1469         364 :   vkum = rnfkummer_initall(bnr, ZV_to_zv(gel(fa,1)), prec);
    1470         364 :   P = bad_primes(bnr);
    1471         364 :   return  gerepilecopy(av, bnrclassfield_H(vkum, bnr, P, subgroup, fa, flag, prec));
    1472             : }

Generated by: LCOV version 1.13