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

Generated by: LCOV version 1.13