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

Generated by: LCOV version 1.16