Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - kummer.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.12.1 lcov report (development 25406-bf255ab81b) Lines: 895 906 98.8 %
Date: 2020-06-04 05:59:24 Functions: 65 65 100.0 %
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.13