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 24964-97b24cd8ac) Lines: 999 1086 92.0 %
Date: 2020-01-22 05:56:41 Functions: 69 71 97.2 %
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        3857 : prank(GEN cyc, long ell)
      44             : {
      45        3857 :   long i, l = lg(cyc);
      46       10654 :   for (i=1; i < l; i++)
      47        7238 :     if (umodiu(gel(cyc,i),ell)) break;
      48        3857 :   return i-1;
      49             : }
      50             : 
      51             : /* increment y, which runs through [0,d-1]^(k-1). Return 0 when done. */
      52             : static int
      53         476 : increment(GEN y, long k, long d)
      54             : {
      55         476 :   long i = k, j;
      56             :   do
      57             :   {
      58         714 :     if (--i == 0) return 0;
      59         630 :     y[i]++;
      60         630 :   } while (y[i] >= d);
      61         392 :   for (j = i+1; j < k; j++) y[j] = 0;
      62         392 :   return 1;
      63             : }
      64             : 
      65             : static int
      66        1890 : ok_congruence(GEN X, ulong ell, long lW, GEN vecMsup)
      67             : {
      68             :   long i, l;
      69        1890 :   l = lg(X);
      70        3213 :   for (i=lW; i<l; i++)
      71        1512 :     if (X[i] == 0) return 0;
      72        1701 :   if (lW >= l && zv_equal0(X)) return 0;
      73        1701 :   l = lg(vecMsup);
      74        2030 :   for (i=1; i<l; i++)
      75         329 :     if (zv_equal0(Flm_Flc_mul(gel(vecMsup,i),X, ell))) return 0;
      76        1701 :   return 1;
      77             : }
      78             : 
      79             : static int
      80        1155 : ok_sign(GEN X, GEN msign, GEN arch)
      81             : {
      82        1155 :   return zv_equal(Flm_Flc_mul(msign, X, 2), arch);
      83             : }
      84             : 
      85             : /* REDUCTION MOD ell-TH POWERS */
      86             : 
      87             : /* make be integral by multiplying by t in (Q^*)^ell */
      88             : static GEN
      89        1554 : reduce_mod_Qell(GEN bnfz, GEN be, GEN gell)
      90             : {
      91             :   GEN c;
      92        1554 :   be = nf_to_scalar_or_basis(bnfz, be);
      93        1554 :   be = Q_primitive_part(be, &c);
      94        1554 :   if (c)
      95             :   {
      96         924 :     GEN d, fa = factor(c);
      97         924 :     gel(fa,2) = FpC_red(gel(fa,2), gell);
      98         924 :     d = factorback(fa);
      99         924 :     be = typ(be) == t_INT? mulii(be,d): ZC_Z_mul(be, d);
     100             :   }
     101        1554 :   return be;
     102             : }
     103             : 
     104             : /* return q, q^n r = x, v_pr(r) < n for all pr */
     105             : static GEN
     106        1554 : idealsqrtn(GEN nf, GEN x, GEN gn)
     107             : {
     108        1554 :   long i, l, n = itos(gn);
     109             :   GEN fa, q, Ex, Pr;
     110             : 
     111        1554 :   fa = idealfactor(nf, x);
     112        1554 :   Pr = gel(fa,1); l = lg(Pr);
     113        1554 :   Ex = gel(fa,2); q = NULL;
     114       10173 :   for (i=1; i<l; i++)
     115             :   {
     116        8619 :     long ex = itos(gel(Ex,i));
     117        8619 :     GEN e = stoi(ex / n);
     118        8619 :     if (q) q = idealmulpowprime(nf, q, gel(Pr,i), e);
     119        1386 :     else   q = idealpow(nf, gel(Pr,i), e);
     120             :   }
     121        1554 :   return q? q: gen_1;
     122             : }
     123             : 
     124             : static GEN
     125        1554 : reducebeta(GEN bnfz, GEN b, GEN ell)
     126             : {
     127        1554 :   GEN y, cb, fu, nf = bnf_get_nf(bnfz);
     128             : 
     129        1554 :   if (DEBUGLEVEL>1) err_printf("reducing beta = %Ps\n",b);
     130        1554 :   b = reduce_mod_Qell(nf, b, ell);
     131             :   /* reduce l-th root */
     132        1554 :   y = idealsqrtn(nf, b, ell); /* (b) = y^ell I, I integral */
     133        1554 :   if (typ(y) == t_MAT && !is_pm1(gcoeff(y,1,1)))
     134             :   {
     135         847 :     GEN T = idealred(nf, mkvec2(y, gen_1)), t = gel(T,2);
     136             :     /* (t)*T[1] = y, T[1] integral and small */
     137         847 :     if (gcmp(idealnorm(nf,t), gen_1) > 0)
     138         623 :       b = nfmul(nf, b, nfpow(nf, t, negi(ell)));
     139             :   }
     140        1554 :   if (DEBUGLEVEL>1) err_printf("beta reduced via ell-th root = %Ps\n",b);
     141        1554 :   b = Q_primitive_part(b, &cb);
     142        1554 :   if (cb)
     143             :   {
     144         525 :     y = nfroots(nf, gsub(monomial(gen_1, itou(ell), fetch_var_higher()),
     145             :                          basistoalg(nf,b)));
     146         525 :     delete_var();
     147             :   }
     148        1554 :   if (cb && lg(y) != 1) b = cb;
     149        1400 :   else if ((fu = bnf_build_cheapfu(bnfz)))
     150             :   { /* log. embeddings of fu^ell */
     151        1400 :     GEN logfu = bnf_get_logfu(bnfz);
     152        1400 :     GEN elllogfu = RgM_Rg_mul(real_i(logfu), ell);
     153        1400 :     long prec = nf_get_prec(nf);
     154             :     for (;;)
     155          14 :     {
     156        1414 :       GEN ex, z = nflogembed(nf, b, NULL, prec);
     157        1414 :       if (z)
     158             :       {
     159        1400 :         ex = RgM_Babai(elllogfu, z);
     160        1400 :         if (ex)
     161             :         {
     162        1400 :           if (ZV_equal0(ex)) break;
     163         393 :           y = nffactorback(nf, fu, RgC_Rg_mul(ex,ell));
     164         393 :           b = nfdiv(nf, b, y); break;
     165             :         }
     166             :       }
     167          14 :       prec = precdbl(prec);
     168          14 :       if (DEBUGLEVEL) pari_warn(warnprec,"reducebeta",prec);
     169          14 :       nf = nfnewprec_shallow(nf,prec);
     170             :     }
     171        1400 :     if (cb) b = gmul(b, cb);
     172             :   }
     173        1554 :   if (DEBUGLEVEL>1) err_printf("beta LLL-reduced mod U^l = %Ps\n",b);
     174        1554 :   return b;
     175             : }
     176             : 
     177             : /* FIXME: remove */
     178             : static GEN
     179        4167 : tauofalg(GEN x, tau_s *tau) {
     180        4167 :   long tx = typ(x);
     181        4167 :   if (tx == t_POLMOD) { x = gel(x,2); tx = typ(x); }
     182        4167 :   if (tx == t_POL) x = RgX_RgXQ_eval(x, tau->x, tau->R);
     183        4167 :   return mkpolmod(x, tau->R);
     184             : }
     185             : 
     186             : struct rnfkummer
     187             : {
     188             :   GEN bnfz, u, vecC, Q, vecW;
     189             :   ulong g, ell;
     190             :   long rc;
     191             :   compo_s COMPO;
     192             :   tau_s tau;
     193             :   toK_s T;
     194             : };
     195             : 
     196             : /* set kum->tau; compute Gal(K(\zeta_l)/K) */
     197             : static void
     198         497 : get_tau(struct rnfkummer *kum)
     199             : { /* compute action of tau: q^g + kp */
     200         497 :   compo_s *C = &kum->COMPO;
     201         497 :   GEN U = RgX_add(RgXQ_powu(C->q, kum->g, C->R), RgX_muls(C->p, C->k));
     202         497 :   kum->tau.x  = RgX_RgXQ_eval(C->rev, U, C->R);
     203         497 :   kum->tau.R  = C->R;
     204         497 :   kum->tau.zk = nfgaloismatrix(bnf_get_nf(kum->bnfz), kum->tau.x);
     205         497 : }
     206             : 
     207             : static GEN tauoffamat(GEN x, tau_s *tau);
     208             : 
     209             : static GEN
     210       52669 : tauofelt(GEN x, tau_s *tau)
     211             : {
     212       52669 :   switch(typ(x))
     213             :   {
     214       46444 :     case t_COL: return RgM_RgC_mul(tau->zk, x);
     215        2058 :     case t_MAT: return tauoffamat(x, tau);
     216        4167 :     default: return tauofalg(x, tau);
     217             :   }
     218             : }
     219             : static GEN
     220        2408 : tauofvec(GEN x, tau_s *tau)
     221             : {
     222             :   long i, l;
     223        2408 :   GEN y = cgetg_copy(x, &l);
     224        2408 :   for (i=1; i<l; i++) gel(y,i) = tauofelt(gel(x,i), tau);
     225        2408 :   return y;
     226             : }
     227             : /* [x, tau(x), ..., tau^(m-1)(x)] */
     228             : static GEN
     229        1225 : powtau(GEN x, long m, tau_s *tau)
     230             : {
     231        1225 :   GEN y = cgetg(m+1, t_VEC);
     232             :   long i;
     233        1225 :   gel(y,1) = x;
     234        1225 :   for (i=2; i<=m; i++) gel(y,i) = tauofelt(gel(y,i-1), tau);
     235        1225 :   return y;
     236             : }
     237             : /* x^lambda */
     238             : static GEN
     239        1302 : lambdaofelt(GEN x, toK_s *T)
     240             : {
     241        1302 :   tau_s *tau = T->tau;
     242        1302 :   long i, m = T->m;
     243        1302 :   GEN y = trivial_fact(), powg = T->powg; /* powg[i] = g^i */
     244        2954 :   for (i=1; i<m; i++)
     245             :   {
     246        1652 :     y = famat_mulpows_shallow(y, x, uel(powg,m-i+1));
     247        1652 :     x = tauofelt(x, tau);
     248             :   }
     249        1302 :   return famat_mul_shallow(y, x);
     250             : }
     251             : static GEN
     252         511 : lambdaofvec(GEN x, toK_s *T)
     253             : {
     254             :   long i, l;
     255         511 :   GEN y = cgetg_copy(x, &l);
     256         511 :   for (i=1; i<l; i++) gel(y,i) = lambdaofelt(gel(x,i), T);
     257         511 :   return y;
     258             : }
     259             : 
     260             : static GEN
     261        2058 : tauoffamat(GEN x, tau_s *tau)
     262             : {
     263        2058 :   return mkmat2(tauofvec(gel(x,1), tau), gel(x,2));
     264             : }
     265             : 
     266             : static GEN
     267         364 : tauofideal(GEN id, tau_s *tau)
     268             : {
     269         364 :   return ZM_hnfmodid(RgM_mul(tau->zk, id), gcoeff(id, 1,1));
     270             : }
     271             : 
     272             : static int
     273         938 : isprimeidealconj(GEN P, GEN Q, tau_s *tau)
     274             : {
     275         938 :   GEN p = pr_get_p(P);
     276         938 :   GEN x = pr_get_gen(P);
     277         938 :   if (!equalii(p, pr_get_p(Q))
     278         784 :    || pr_get_e(P) != pr_get_e(Q)
     279         784 :    || pr_get_f(P) != pr_get_f(Q)) return 0;
     280         756 :   if (ZV_equal(x, pr_get_gen(Q))) return 1;
     281             :   for(;;)
     282             :   {
     283        2674 :     if (ZC_prdvd(x,Q)) return 1;
     284        1267 :     x = FpC_red(tauofelt(x, tau), p);
     285        1267 :     if (ZC_prdvd(x,P)) return 0;
     286             :   }
     287             : }
     288             : 
     289             : static int
     290        3668 : isconjinprimelist(GEN S, GEN pr, tau_s *tau)
     291             : {
     292             :   long i, l;
     293             : 
     294        3668 :   if (!tau) return 0;
     295        1442 :   l = lg(S);
     296        1932 :   for (i=1; i<l; i++)
     297         938 :     if (isprimeidealconj(gel(S,i),pr,tau)) return 1;
     298         994 :   return 0;
     299             : }
     300             : 
     301             : /* assume x in basistoalg form */
     302             : static GEN
     303        1673 : downtoK(toK_s *T, GEN x)
     304             : {
     305        1673 :   long degKz = lg(T->invexpoteta1) - 1;
     306        1673 :   GEN y = gmul(T->invexpoteta1, Rg_to_RgC(lift_shallow(x), degKz));
     307        1673 :   return gmodulo(gtopolyrev(y,varn(T->polnf)), T->polnf);
     308             : }
     309             : 
     310             : static GEN
     311           0 : no_sol(long all, long i)
     312             : {
     313           0 :   if (!all) pari_err_BUG(stack_sprintf("kummer [bug%ld]", i));
     314           0 :   return cgetg(1,t_VEC);
     315             : }
     316             : 
     317             : static GEN
     318         987 : get_gell(GEN bnr, GEN subgp, long all)
     319             : {
     320         987 :   if (all && all != -1) return utoipos(labs(all));
     321         952 :   if (!subgp) return ZV_prod(bnr_get_cyc(bnr));
     322         952 :   return det(subgp);
     323             : }
     324             : 
     325             : typedef struct {
     326             :   GEN Sm, Sml1, Sml2, Sl, ESml2;
     327             : } primlist;
     328             : 
     329             : static int
     330        1519 : build_list_Hecke(primlist *L, GEN nfz, GEN fa, GEN gothf, long ell, tau_s *tau)
     331             : {
     332             :   GEN listpr, listex, pr, factell;
     333        1519 :   long vp, i, l, degKz = nf_get_degree(nfz);
     334             : 
     335        1519 :   if (!fa) fa = idealfactor(nfz, gothf);
     336        1519 :   listpr = gel(fa,1);
     337        1519 :   listex = gel(fa,2); l = lg(listpr);
     338        1519 :   L->Sm  = vectrunc_init(l);
     339        1519 :   L->Sml1= vectrunc_init(l);
     340        1519 :   L->Sml2= vectrunc_init(l);
     341        1519 :   L->Sl  = vectrunc_init(l+degKz);
     342        1519 :   L->ESml2=vecsmalltrunc_init(l);
     343        3374 :   for (i=1; i<l; i++)
     344             :   {
     345        1855 :     pr = gel(listpr,i);
     346        1855 :     vp = itos(gel(listex,i));
     347        1855 :     if (!equaliu(pr_get_p(pr), ell))
     348             :     {
     349        1484 :       if (vp != 1) return 1;
     350        1484 :       if (!isconjinprimelist(L->Sm,pr,tau)) vectrunc_append(L->Sm,pr);
     351             :     }
     352             :     else
     353             :     {
     354         371 :       long e = pr_get_e(pr), vd = (vp-1)*(ell-1)-ell*e;
     355         371 :       if (vd > 0) return 4;
     356         371 :       if (vd==0)
     357             :       {
     358          91 :         if (!isconjinprimelist(L->Sml1,pr,tau)) vectrunc_append(L->Sml1, pr);
     359             :       }
     360             :       else
     361             :       {
     362         280 :         if (vp==1) return 2;
     363         280 :         if (!isconjinprimelist(L->Sml2,pr,tau))
     364             :         {
     365         280 :           vectrunc_append(L->Sml2, pr);
     366         280 :           vecsmalltrunc_append(L->ESml2, vp);
     367             :         }
     368             :       }
     369             :     }
     370             :   }
     371        1519 :   factell = idealprimedec(nfz,utoipos(ell)); l = lg(factell);
     372        3703 :   for (i=1; i<l; i++)
     373             :   {
     374        2184 :     pr = gel(factell,i);
     375        2184 :     if (!idealval(nfz,gothf,pr) && !isconjinprimelist(L->Sl,pr,tau))
     376        1806 :       vectrunc_append(L->Sl, pr);
     377             :   }
     378        1519 :   return 0; /* OK */
     379             : }
     380             : 
     381             : /* Return a Flm */
     382             : static GEN
     383        2366 : logall(GEN nf, GEN vec, long lW, long mginv, long ell, GEN pr, long ex)
     384             : {
     385        2366 :   pari_sp av = avma;
     386        2366 :   GEN m, M, sprk = log_prk_init(nf, pr, ex);
     387        2366 :   long ellrank, i, l = lg(vec);
     388             : 
     389        2366 :   ellrank = prank(gel(sprk,1), ell);
     390        2366 :   M = cgetg(l,t_MAT);
     391       10689 :   for (i=1; i<l; i++)
     392             :   {
     393        8323 :     m = log_prk(nf, gel(vec,i), sprk);
     394        8323 :     setlg(m, ellrank+1);
     395        8323 :     if (i < lW) m = gmulsg(mginv, m);
     396        8323 :     gel(M,i) = ZV_to_Flv(m, ell);
     397             :   }
     398        2366 :   return gerepilecopy(av, M);
     399             : }
     400             : 
     401             : /* compute the u_j (see remark 5.2.15.) */
     402             : static GEN
     403        1491 : get_u(GEN cyc, long rc, ulong ell)
     404             : {
     405        1491 :   long i, l = lg(cyc);
     406        1491 :   GEN u = cgetg(l,t_VECSMALL);
     407        1491 :   for (i=1; i<=rc; i++) uel(u,i) = 0;
     408        1491 :   for (   ; i < l; i++) uel(u,i) = Fl_inv(umodiu(gel(cyc,i), ell), ell);
     409        1491 :   return u;
     410             : }
     411             : 
     412             : /* alg. 5.2.15. with remark */
     413             : static void
     414        1498 : isprincipalell(GEN bnfz, GEN id, GEN cycgen, GEN u, ulong ell, long rc,
     415             :                GEN *pv, GEN *pb)
     416             : {
     417        1498 :   long i, l = lg(cycgen);
     418        1498 :   GEN v, b, db, y = bnfisprincipal0(bnfz, id, nf_FORCE|nf_GENMAT);
     419             : 
     420        1498 :   v = ZV_to_Flv(gel(y,1), ell);
     421        1498 :   b = gel(y,2);
     422        1498 :   if (typ(b) == t_COL)
     423             :   {
     424        1448 :     b = Q_remove_denom(gel(y,2), &db);
     425        1448 :     if (db) b = famat_mulpows_shallow(b, db, -1);
     426             :   }
     427        1743 :   for (i = rc+1; i < l; i++)
     428             :   {
     429         245 :     ulong e = Fl_mul( uel(v,i), uel(u,i), ell);
     430         245 :     b = famat_mulpows_shallow(b, gel(cycgen,i), e);
     431             :   }
     432        1498 :   setlg(v,rc+1); *pv = v; *pb = b;
     433        1498 : }
     434             : 
     435             : static GEN
     436        1554 : compute_beta(GEN X, GEN vecWB, GEN ell, GEN bnfz)
     437             : {
     438             :   GEN BE, be;
     439        1554 :   BE = famat_reduce(famatV_zv_factorback(vecWB, X));
     440        1554 :   gel(BE,2) = centermod(gel(BE,2), ell);
     441        1554 :   be = nffactorback(bnfz, BE, NULL);
     442        1554 :   be = reducebeta(bnfz, be, ell);
     443        1554 :   if (DEBUGLEVEL>1) err_printf("beta reduced = %Ps\n",be);
     444        1554 :   return be;
     445             : }
     446             : 
     447             : static GEN
     448        1491 : futu(GEN bnf)
     449             : {
     450        1491 :   GEN fu, tu, SUnits = bnf_get_sunits(bnf);
     451        1491 :   if (SUnits)
     452             :   {
     453         833 :     tu = nf_to_scalar_or_basis(bnf_get_nf(bnf), bnf_get_tuU(bnf));
     454         833 :     fu = bnf_compactfu(bnf);
     455             :   }
     456             :   else
     457             :   {
     458         658 :     GEN U = bnf_build_units(bnf);
     459         658 :     tu = gel(U,1); fu = vecslice(U, 2, lg(U)-1);
     460             :   }
     461        1491 :   return vec_append(fu, tu);
     462             : }
     463             : static GEN
     464        1491 : get_Selmer(GEN bnf, GEN cycgen, long rc)
     465        1491 : { return shallowconcat(futu(bnf), vecslice(cycgen,1,rc)); }
     466             : 
     467             : GEN
     468       56371 : lift_if_rational(GEN x)
     469             : {
     470             :   long lx, i;
     471             :   GEN y;
     472             : 
     473       56371 :   switch(typ(x))
     474             :   {
     475        8764 :     default: break;
     476             : 
     477             :     case t_POLMOD:
     478       32466 :       y = gel(x,2);
     479       32466 :       if (typ(y) == t_POL)
     480             :       {
     481       12019 :         long d = degpol(y);
     482       12019 :         if (d > 0) return x;
     483        2219 :         return (d < 0)? gen_0: gel(y,2);
     484             :       }
     485       20447 :       return y;
     486             : 
     487        6895 :     case t_POL: lx = lg(x);
     488        6895 :       for (i=2; i<lx; i++) gel(x,i) = lift_if_rational(gel(x,i));
     489        6895 :       break;
     490        8246 :     case t_VEC: case t_COL: case t_MAT: lx = lg(x);
     491        8246 :       for (i=1; i<lx; i++) gel(x,i) = lift_if_rational(gel(x,i));
     492             :   }
     493       23905 :   return x;
     494             : }
     495             : 
     496             : /* A column vector representing a subgroup of prime index */
     497             : static GEN
     498           0 : grptocol(GEN H)
     499             : {
     500           0 :   long i, j, l = lg(H);
     501           0 :   GEN col = cgetg(l, t_VECSMALL);
     502           0 :   for (i = 1; i < l; i++)
     503             :   {
     504           0 :     ulong ell = itou( gcoeff(H,i,i) );
     505           0 :     if (ell == 1) col[i] = 0; else { col[i] = ell-1; break; }
     506             :   }
     507           0 :   for (j=i; ++j < l; ) col[j] = itou( gcoeff(H,i,j) );
     508           0 :   return col;
     509             : }
     510             : 
     511             : /* Reorganize kernel basis so that the tests of ok_congruence can be ok
     512             :  * for y[ncyc]=1 and y[1..ncyc]=1 */
     513             : static GEN
     514          14 : fix_kernel(GEN K, GEN M, GEN vecMsup, long lW, long ell)
     515             : {
     516          14 :   pari_sp av = avma;
     517          14 :   long i, j, idx, ffree, dK = lg(K)-1;
     518          14 :   GEN Ki, Kidx = cgetg(dK+1, t_VECSMALL);
     519             : 
     520             :   /* First step: Gauss elimination on vectors lW...lg(M)-1 */
     521          28 :   for (idx = lg(K), i = lg(M)-1; i >= lW; i--)
     522             :   {
     523          14 :     for (j = dK; j > 0; j--) if (coeff(K, i, j)) break;
     524          14 :     if (!j || j == dK) continue;
     525             :     /* ensure that K[i,dK] != 0 */
     526           0 :     for (j = idx; j < dK; j++)
     527           0 :       if (coeff(K, i, j) && coeff(K, Kidx[j], dK) != ell - 1)
     528           0 :         Flv_add_inplace(gel(K,dK), gel(K,j), ell);
     529           0 :     idx--;
     530           0 :     if (j != idx) swap(gel(K, j), gel(K, idx));
     531           0 :     Kidx[idx] = i;
     532           0 :     if (coeff(K,i,idx) != 1)
     533           0 :       Flv_Fl_div_inplace(gel(K,idx), coeff(K,i,idx), ell);
     534           0 :     Ki = gel(K,idx);
     535           0 :     if (coeff(K,i,dK) != 1)
     536             :     {
     537           0 :       ulong t = Fl_sub(coeff(K,i,dK), 1, ell);
     538           0 :       Flv_sub_inplace(gel(K,dK), Flv_Fl_mul(Ki, t, ell), ell);
     539             :     }
     540           0 :     for (j = dK; --j > 0; )
     541             :     {
     542           0 :       if (j == idx) continue;
     543           0 :       if (coeff(K,i,j))
     544           0 :         Flv_sub_inplace(gel(K,j), Flv_Fl_mul(Ki, coeff(K,i,j), ell), ell);
     545             :     }
     546             :   }
     547             :   /* ffree = first vector that is not "free" for the scalar products */
     548          14 :   ffree = idx;
     549             :   /* Second step: for each hyperplane equation in vecMsup, do the same
     550             :    * thing as before. */
     551          14 :   for (i=1; i < lg(vecMsup); i++)
     552             :   {
     553           0 :     GEN Msup = gel(vecMsup,i);
     554             :     ulong dotprod;
     555           0 :     if (lgcols(Msup) != 2) continue;
     556           0 :     Msup = zm_row(Msup, 1);
     557           0 :     for (j=ffree; j > 0; j--)
     558             :     {
     559           0 :       dotprod = Flv_dotproduct(Msup, gel(K,j), ell);
     560           0 :       if (dotprod)
     561             :       {
     562           0 :         if (j != --ffree) swap(gel(K, j), gel(K, ffree));
     563           0 :         if (dotprod != 1) Flv_Fl_div_inplace(gel(K, ffree), dotprod, ell);
     564           0 :         break;
     565             :       }
     566             :     }
     567           0 :     if (!j)
     568             :     { /* Do our best to ensure that vecMsup.K[dK] != 0 */
     569           0 :       if (Flv_dotproduct(Msup, gel(K,dK), ell) == 0)
     570             :       {
     571           0 :         for (j = ffree-1; j <= dK; j++)
     572           0 :           if (Flv_dotproduct(Msup, gel(K,j), ell)
     573           0 :               && coeff(K,Kidx[j],dK) != ell-1)
     574           0 :             Flv_add_inplace(gel(K,dK), gel(K,j), ell);
     575             :       }
     576           0 :       continue;
     577             :     }
     578           0 :     Ki = gel(K,ffree);
     579           0 :     dotprod = Flv_dotproduct(Msup, gel(K,dK), ell);
     580           0 :     if (dotprod != 1)
     581             :     {
     582           0 :       ulong t = Fl_sub(dotprod,1,ell);
     583           0 :       Flv_sub_inplace(gel(K,dK), Flv_Fl_mul(Ki,t,ell), ell);
     584             :     }
     585           0 :     for (j = dK; j > 0; j--)
     586             :     {
     587           0 :       if (j == ffree) continue;
     588           0 :       dotprod = Flv_dotproduct(Msup, gel(K,j), ell);
     589           0 :       if (dotprod) Flv_sub_inplace(gel(K,j), Flv_Fl_mul(Ki,dotprod,ell), ell);
     590             :     }
     591             :   }
     592          14 :   if (ell == 2)
     593             :   {
     594          14 :     for (i = ffree, j = ffree-1; i <= dK && j; i++, j--)
     595           0 :     { swap(gel(K,i), gel(K,j)); }
     596             :   }
     597             :   /* Try to ensure that y = vec_ei(n, i) gives a good candidate */
     598          14 :   for (i = 1; i < dK; i++) Flv_add_inplace(gel(K,i), gel(K,dK), ell);
     599          14 :   return gerepilecopy(av, K);
     600             : }
     601             : 
     602             : static GEN
     603          14 : Flm_init(long m, long n)
     604             : {
     605          14 :   GEN M = cgetg(n+1, t_MAT);
     606          14 :   long i; for (i = 1; i <= n; i++) gel(M,i) = cgetg(m+1, t_VECSMALL);
     607          14 :   return M;
     608             : }
     609             : static void
     610         308 : Flv_fill(GEN v, GEN y)
     611             : {
     612         308 :   long i, l = lg(y);
     613         308 :   for (i = 1; i < l; i++) v[i] = y[i];
     614         308 : }
     615             : 
     616             : static GEN
     617        2002 : get_badbnf(GEN bnf)
     618             : {
     619             :   long i, l;
     620        2002 :   GEN bad = gen_1, gen = bnf_get_gen(bnf);
     621        2002 :   l = lg(gen);
     622        3626 :   for (i = 1; i < l; i++)
     623             :   {
     624        1624 :     GEN g = gel(gen,i);
     625        1624 :     bad = lcmii(bad, gcoeff(g,1,1));
     626             :   }
     627        2002 :   return bad;
     628             : }
     629             : /* Let K base field, L/K described by bnr (conductor f) + H. Return a list of
     630             :  * primes coprime to f*ell of degree 1 in K whose images in Cl_f(K) generate H:
     631             :  * thus they all split in Lz/Kz; t in Kz is such that
     632             :  * t^(1/p) generates Lz => t is an ell-th power in k(pr) for all such primes.
     633             :  * Restrict to primes not dividing
     634             :  * - the index fz of the polynomial defining Kz, or
     635             :  * - the modulus, or
     636             :  * - ell, or
     637             :  * - a generator in bnf.gen or bnfz.gen */
     638             : static GEN
     639        1491 : get_prlist(GEN bnr, GEN H, ulong ell, GEN bnfz)
     640             : {
     641        1491 :   pari_sp av0 = avma;
     642             :   forprime_t T;
     643             :   ulong p;
     644             :   GEN L, nf, cyc, bad, cond, condZ, Hsofar;
     645        1491 :   L = cgetg(1, t_VEC);
     646        1491 :   cyc = bnr_get_cyc(bnr);
     647        1491 :   nf = bnr_get_nf(bnr);
     648             : 
     649        1491 :   cond = gel(bnr_get_mod(bnr), 1);
     650        1491 :   condZ = gcoeff(cond,1,1);
     651        1491 :   bad = get_badbnf(bnr_get_bnf(bnr));
     652        1491 :   if (bnfz)
     653             :   {
     654         511 :     GEN badz = lcmii(get_badbnf(bnfz), nf_get_index(bnf_get_nf(bnfz)));
     655         511 :     bad = mulii(bad,badz);
     656             :   }
     657        1491 :   bad = lcmii(muliu(condZ, ell), bad);
     658             :   /* restrict to primes not dividing bad */
     659             : 
     660        1491 :   u_forprime_init(&T, 2, ULONG_MAX);
     661        1491 :   Hsofar = cgetg(1, t_MAT);
     662       21098 :   while ((p = u_forprime_next(&T)))
     663             :   {
     664             :     GEN LP;
     665             :     long i, l;
     666       19607 :     if (p == ell || !umodiu(bad, p)) continue;
     667       15931 :     LP = idealprimedec_limit_f(nf, utoipos(p), 1);
     668       15931 :     l = lg(LP);
     669       24703 :     for (i = 1; i < l; i++)
     670             :     {
     671       10263 :       pari_sp av = avma;
     672       10263 :       GEN M, P = gel(LP,i), v = bnrisprincipal(bnr, P, 0);
     673       10263 :       if (!hnf_invimage(H, v)) { set_avma(av); continue; }
     674        3031 :       M = shallowconcat(Hsofar, v);
     675        3031 :       M = ZM_hnfmodid(M, cyc);
     676        3031 :       if (ZM_equal(M, Hsofar)) continue;
     677        2345 :       L = vec_append(L, P);
     678        2345 :       Hsofar = M;
     679             :       /* the primes in L generate H */
     680        2345 :       if (ZM_equal(M, H)) return gerepilecopy(av0, L);
     681             :     }
     682             :   }
     683           0 :   pari_err_BUG("rnfkummer [get_prlist]");
     684           0 :   return NULL;
     685             : }
     686             : /*Lprz list of prime ideals in Kz that must split completely in Lz/Kz, vecWA
     687             :  * generators for the S-units used to build the Kummer generators. Return
     688             :  * matsmall M such that \prod WA[j]^x[j] ell-th power mod pr[i] iff
     689             :  * \sum M[i,j] x[j] = 0 (mod ell) */
     690             : static GEN
     691        1491 : subgroup_info(GEN bnfz, GEN Lprz, long ell, GEN vecWA)
     692             : {
     693        1491 :   GEN nfz = bnf_get_nf(bnfz), M, gell = utoipos(ell), Lell = mkvec(gell);
     694        1491 :   long i, j, l = lg(vecWA), lz = lg(Lprz);
     695        1491 :   M = cgetg(l, t_MAT);
     696        1491 :   for (j=1; j<l; j++) gel(M,j) = cgetg(lz, t_VECSMALL);
     697        3836 :   for (i=1; i < lz; i++)
     698             :   {
     699        2345 :     GEN pr = gel(Lprz,i), EX = subiu(pr_norm(pr), 1);
     700        2345 :     GEN N, g,T,p, prM = idealhnf(nfz, pr);
     701        2345 :     GEN modpr = zk_to_Fq_init(nfz, &pr,&T,&p);
     702        2345 :     long v = Z_lvalrem(divis(EX,ell), ell, &N) + 1; /* Norm(pr)-1 = N * ell^v */
     703        2345 :     GEN ellv = powuu(ell, v);
     704        2345 :     g = gener_Fq_local(T,p, Lell);
     705        2345 :     g = Fq_pow(g,N, T,p); /* order ell^v */
     706       11186 :     for (j=1; j < l; j++)
     707             :     {
     708        8841 :       GEN logc, c = gel(vecWA,j);
     709        8841 :       if (typ(c) == t_MAT) /* famat */
     710        5103 :         c = famat_makecoprime(nfz, gel(c,1), gel(c,2), pr, prM, EX);
     711        8841 :       c = nf_to_Fq(nfz, c, modpr);
     712        8841 :       c = Fq_pow(c, N, T,p);
     713        8841 :       logc = Fq_log(c, g, ellv, T,p);
     714        8841 :       ucoeff(M, i,j) = umodiu(logc, ell);
     715             :     }
     716             :   }
     717        1491 :   return M;
     718             : }
     719             : 
     720             : /* if all>0, give all equations of degree 'all'. Assume bnr modulus is the
     721             :  * conductor */
     722             : static GEN
     723         994 : rnfkummersimple(GEN bnr, GEN subgroup, long ell, long all)
     724             : {
     725         994 :   long i, j, degK, dK, lSml2, lSl2, lSp, rc, lW, prec, rk = 0, ncyc = 0;
     726         994 :   long firstpass = all<0;
     727             :   GEN bnf, nf,bid, ideal, arch, cycgen, cyc, Sp, prSp, matP;
     728             :   GEN gell, xell, u, M, K, y, vecMsup, vecW, vecWB, vecBp, msign;
     729         994 :   GEN mat = NULL, matgrp = NULL, be1 = NULL, res = NULL;
     730             :   primlist L;
     731             : 
     732         994 :   bnf = bnr_get_bnf(bnr); if (!bnf_get_sunits(bnf)) bnf_build_units(bnf);
     733         994 :   nf  = bnf_get_nf(bnf);
     734         994 :   degK = nf_get_degree(nf);
     735             : 
     736         994 :   bid = bnr_get_bid(bnr);
     737         994 :   ideal= bid_get_ideal(bid);
     738         994 :   arch = bid_get_arch(bid); /* this is the conductor */
     739         994 :   i = build_list_Hecke(&L, nf, bid_get_fact2(bid), ideal, ell, NULL);
     740         994 :   if (i) return no_sol(all,i);
     741             : 
     742         994 :   lSml2 = lg(L.Sml2);
     743         994 :   Sp = shallowconcat(L.Sm, L.Sml1); lSp = lg(Sp);
     744         994 :   prSp = shallowconcat(L.Sml2, L.Sl); lSl2 = lg(prSp);
     745             : 
     746         994 :   cycgen = bnf_build_cycgen(bnf);
     747         994 :   cyc = bnf_get_cyc(bnf); rc = prank(cyc, ell);
     748             : 
     749         994 :   vecW = get_Selmer(bnf, cycgen, rc);
     750         994 :   u = get_u(cyc, rc, ell);
     751             : 
     752         994 :   vecBp = cgetg(lSp, t_VEC);
     753         994 :   matP  = cgetg(lSp, t_MAT);
     754        1841 :   for (j = 1; j < lSp; j++)
     755         847 :     isprincipalell(bnf,gel(Sp,j), cycgen,u,ell,rc, &gel(matP,j), &gel(vecBp,j));
     756         994 :   vecWB = shallowconcat(vecW, vecBp);
     757             : 
     758         994 :   prec = DEFAULTPREC +
     759         994 :       nbits2extraprec(((degK-1) * (gexpo(vecWB) + gexpo(nf_get_M(nf)))));
     760         994 :   if (nf_get_prec(nf) < prec) nf = nfnewprec_shallow(nf, prec);
     761         994 :   msign = nfsign(nf, vecWB);
     762         994 :   arch = ZV_to_zv(arch);
     763             : 
     764         994 :   vecMsup = cgetg(lSml2,t_VEC);
     765         994 :   M = NULL;
     766        2373 :   for (i = 1; i < lSl2; i++)
     767             :   {
     768        1379 :     GEN pr = gel(prSp,i);
     769        1379 :     long e = pr_get_e(pr), z = ell * (e / (ell-1));
     770             : 
     771        1379 :     if (i < lSml2)
     772             :     {
     773         140 :       z += 1 - L.ESml2[i];
     774         140 :       gel(vecMsup,i) = logall(nf, vecWB, 0,0, ell, pr,z+1);
     775             :     }
     776        1379 :     M = vconcat(M, logall(nf, vecWB, 0,0, ell, pr,z));
     777             :   }
     778         994 :   lW = lg(vecW);
     779         994 :   M = vconcat(M, shallowconcat(zero_Flm(rc,lW-1), matP));
     780         994 :   if (!all)
     781             :   { /* primes landing in subgroup must be totally split */
     782         980 :     GEN Lpr = get_prlist(bnr, subgroup, ell, NULL);
     783         980 :     GEN M2 = subgroup_info(bnf, Lpr, ell, vecWB);
     784         980 :     M = vconcat(M, M2);
     785             :   }
     786         994 :   K = Flm_ker(M, ell);
     787         994 :   if (all < 0) K = fix_kernel(K, M, vecMsup, lW, ell);
     788         994 :   dK = lg(K)-1;
     789         994 :   y = cgetg(dK+1,t_VECSMALL);
     790         994 :   if (all) res = cgetg(1,t_VEC); /* in case all = 1 */
     791         994 :   if (all < 0)
     792             :   {
     793          14 :     ncyc = dK; rk = 0; mat = Flm_init(dK, ncyc);
     794          14 :     if (all == -1) matgrp = Flm_init(lg(bnr_get_cyc(bnr)), ncyc+1);
     795             :   }
     796         994 :   xell = pol_xn(ell, 0);
     797         994 :   gell = utoipos(ell);
     798             :   do {
     799        1008 :     dK = lg(K)-1;
     800        2107 :     while (dK)
     801             :     {
     802        1071 :       for (i=1; i<dK; i++) y[i] = 0;
     803        1071 :       y[i] = 1; /* y = [0,...,0,1,0,...,0], 1 at i'th position */
     804             :       do
     805             :       {
     806        1449 :         pari_sp av = avma;
     807        1449 :         GEN be, P=NULL, X;
     808        1449 :         if (all < 0)
     809             :         {
     810         308 :           Flv_fill(gel(mat, rk+1), y);
     811         308 :           setlg(mat, rk+2);
     812         308 :           if (Flm_rank(mat, ell) <= rk) continue;
     813             :         }
     814        1344 : FOUND:  X = Flm_Flc_mul(K, y, ell);
     815        1344 :         if (ok_congruence(X, ell, lW, vecMsup) && ok_sign(X, msign, arch))
     816             :         {/* be satisfies all congruences, x^ell - be is irreducible, signature
     817             :           * and relative discriminant are correct */
     818        1008 :           if (all < 0) rk++;
     819        1008 :           be = compute_beta(X, vecWB, gell, bnf);
     820        1008 :           be = nf_to_scalar_or_alg(nf, be);
     821        1008 :           if (typ(be) == t_POL) be = mkpolmod(be, nf_get_pol(nf));
     822        1008 :           if (all == -1)
     823             :           {
     824           0 :             pari_sp av2 = avma;
     825           0 :             GEN Kgrp, colgrp = grptocol(rnfnormgroup(bnr, gsub(xell, be)));
     826           0 :             if (ell != 2)
     827             :             {
     828           0 :               if (rk == 1) be1 = be;
     829             :               else
     830             :               { /* Compute the pesky scalar */
     831           0 :                 GEN K2, C = cgetg(4, t_MAT);
     832           0 :                 gel(C,1) = gel(matgrp,1);
     833           0 :                 gel(C,2) = colgrp;
     834           0 :                 gel(C,3) = grptocol(rnfnormgroup(bnr, gsub(xell, gmul(be1,be))));
     835           0 :                 K2 = Flm_ker(C, ell);
     836           0 :                 if (lg(K2) != 2) pari_err_BUG("linear algebra");
     837           0 :                 K2 = gel(K2,1);
     838           0 :                 if (K2[1] != K2[2])
     839           0 :                   Flv_Fl_mul_inplace(colgrp, Fl_div(K2[2],K2[1],ell), ell);
     840             :               }
     841             :             }
     842           0 :             Flv_fill(gel(matgrp,rk), colgrp);
     843           0 :             setlg(matgrp, rk+1);
     844           0 :             Kgrp = Flm_ker(matgrp, ell);
     845           0 :             if (lg(Kgrp) == 2)
     846             :             {
     847           0 :               setlg(gel(Kgrp,1), rk+1);
     848           0 :               y = Flm_Flc_mul(mat, gel(Kgrp,1), ell);
     849           0 :               all = 0; goto FOUND;
     850             :             }
     851           0 :             set_avma(av2);
     852             :           }
     853             :           else
     854             :           {
     855        1008 :             P = gsub(xell, be);
     856        1008 :             if (all)
     857          28 :               res = shallowconcat(res, gerepileupto(av, P));
     858             :             else
     859             :             {
     860         980 :               if (ZM_equal(rnfnormgroup(bnr,P),subgroup)) return P; /*DONE*/
     861           0 :               set_avma(av); continue;
     862             :             }
     863             :           }
     864          28 :           if (all < 0 && rk == ncyc) return res;
     865          28 :           if (firstpass) break;
     866             :         }
     867         336 :         else set_avma(av);
     868         441 :       } while (increment(y, dK, ell));
     869          91 :       y[dK--] = 0;
     870             :     }
     871          28 :   } while (firstpass--);
     872          14 :   return all? res: gen_0;
     873             : }
     874             : 
     875             : static ulong
     876       15960 : nf_to_logFl(GEN nf, GEN x, GEN modpr, ulong g, ulong q, ulong ell, ulong p)
     877             : {
     878       15960 :   x = nf_to_Fp_coprime(nf, x, modpr);
     879       15960 :   return Fl_log(Fl_powu(umodiu(x, p), q, p), g, ell, p);
     880             : }
     881             : static GEN
     882        4984 : nfV_to_logFlv(GEN nf, GEN x, GEN modpr, ulong g, ulong q, ulong ell, ulong p)
     883        4984 : { pari_APPLY_long(nf_to_logFl(nf, gel(x,i), modpr, g, q, ell, p)); }
     884             : 
     885             : /* Compute e_1 Cl(K)/Cl(K)^ell. If u = w^ell a virtual unit, compute
     886             :  * discrete log mod ell on units.gen + bnf.gen (efficient variant of algo
     887             :  * 5.3.11) by finding ru degree 1 primes Pj coprime to everything, and gj
     888             :  * in k(Pj)^* of order ell such that
     889             :  *      log_gj(u_i^((Pj.p - 1) / ell) mod Pj), j = 1..ru
     890             :  * has maximal F_ell rank ru then solve linear system */
     891             : static GEN
     892         497 : kervirtualunit(struct rnfkummer *kum, GEN vselmer)
     893             : {
     894         497 :   GEN bnf = kum->bnfz, cyc = bnf_get_cyc(bnf), nf = bnf_get_nf(bnf);
     895         497 :   GEN B, vy, vz, M, U1, U2, vtau, vell, SUnits = bnf_get_sunits(bnf);
     896         497 :   long i, j, r, l = lg(vselmer), rc = kum->rc, ru = l-1 - rc, ell = kum->ell;
     897         497 :   long LIMC = SUnits? itou(gel(SUnits,4)): 1;
     898             :   ulong p;
     899             :   forprime_t T;
     900             : 
     901         497 :   vtau = cgetg(l, t_VEC);
     902         497 :   vell = cgetg(l, t_VEC);
     903        2212 :   for (j = 1; j < l; j++)
     904             :   {
     905        1715 :     GEN t = gel(vselmer,j);
     906        1715 :     if (typ(t) == t_MAT)
     907             :     {
     908             :       GEN ct;
     909        1218 :       t = nffactorback(bnf, gel(t,1), ZV_to_Flv(gel(t,2), ell));
     910        1218 :       t = Q_primitive_part(t, &ct);
     911        1218 :       if (ct)
     912             :       {
     913         665 :         GEN F = Q_factor(ct);
     914         665 :         ct = factorback2(gel(F,1), ZV_to_Flv(gel(F,2), ell));
     915         665 :         t = (typ(t) == t_INT)? ct: ZC_Z_mul(t, ct);
     916             :       }
     917             :     }
     918        1715 :     gel(vell,j) = t; /* integral, not to far from primitive */
     919        1715 :     gel(vtau,j) = tauofelt(t, &kum->tau);
     920             :   }
     921         497 :   U1 = vecslice(vell, 1, ru); /* units */
     922         497 :   U2 = vecslice(vell, ru+1, ru+rc); /* cycgen (mod ell-th powers) */
     923         497 :   B = nf_get_index(nf); /* bad primes; from 1 to ru are LIMC-units */
     924         497 :   for (i = 1; i <= rc; i++) B = mulii(B, nfnorm(nf, gel(U2,i)));
     925         497 :   if (LIMC > 1)
     926             :   {
     927         497 :     GEN F = absZ_factor_limit(B, LIMC), P = gel(F,1);
     928         497 :     long lP = lg(P);
     929         497 :     B = (lP > 1)? gel(P,lP-1): gen_1;
     930             :   }
     931         497 :   vy = cgetg(l, t_MAT);
     932         497 :   for (j = 1; j <= ru; j++) gel(vy,j) = zero_Flv(rc); /* units */
     933         861 :   for (     ; j < l; j++)
     934             :   {
     935         364 :     GEN y, w, u = gel(vtau, j); /* virtual unit */
     936         364 :     if (!idealispower(nf, u, ell, &w)) pari_err_BUG("kervirtualunit");
     937         364 :     y = isprincipal(bnf, w); setlg(y, rc+1);
     938         364 :     if (!ZV_equal0(y))
     939         980 :       for (i = 1; i <= rc; i++)
     940         616 :         gel(y,i) = diviiexact(mului(ell,gel(y,i)), gel(cyc,i));
     941         364 :     gel(vy,j) = ZV_to_Flv(y, ell);
     942             :   }
     943         497 :   u_forprime_arith_init(&T, LIMC+1, ULONG_MAX, 1, ell);
     944         497 :   M = cgetg(ru+1, t_MAT); r = 1; setlg(M,2);
     945         497 :   vz = cgetg(ru+1, t_MAT);
     946        2541 :   while ((p = u_forprime_next(&T))) if (umodiu(B,p))
     947             :   {
     948        2009 :     GEN P = idealprimedec_limit_f(nf, utoipos(p), 1);
     949        2009 :     long nP = lg(P)-1;
     950        2009 :     ulong g = rootsof1_Fl(ell, p), q = p / ell; /* (p-1) / ell */
     951        3794 :     for (i = 1; i <= nP; i++)
     952             :     {
     953        2282 :       GEN modpr = zkmodprinit(nf, gel(P,i));
     954             :       GEN z, v2;
     955        2282 :       gel(M, r) = nfV_to_logFlv(nf, U1, modpr, g, q, ell, p); /* log futu */
     956        2282 :       if (Flm_rank(M, ell) < r) continue; /* discard */
     957             : 
     958        1351 :       v2 = nfV_to_logFlv(nf, U2, modpr, g, q, ell, p); /* log alpha[1..rc] */
     959        1351 :       gel(vz, r) = z = nfV_to_logFlv(nf, vtau, modpr, g, q, ell, p);
     960        2443 :       for (j = ru+1; j < l; j++)
     961        1092 :         uel(z,j) = Fl_sub(uel(z,j), Flv_dotproduct(v2, gel(vy,j), ell), ell);
     962        1351 :       if (r == ru) break;
     963         854 :       r++; setlg(M, r+1);
     964             :     }
     965        2009 :     if (i < nP) break;
     966             :   }
     967         497 :   if (r != ru) pari_err_BUG("kervirtualunit");
     968             :   /* Solve prod_k U[k]^x[j,k] = vtau[j] / prod_i alpha[i]^vy[j,i] mod (K^*)^ell
     969             :    * for 1 <= j <= #vtau. I.e. for a fixed j: M x[j] = vz[j] (mod ell) */
     970         497 :   M = Flm_inv(Flm_transpose(M), ell);
     971         497 :   vz = Flm_transpose(vz); /* now ru x #vtau */
     972        2212 :   for (j = 1; j < l; j++)
     973        1715 :     gel(vy,j) = shallowconcat(Flm_Flc_mul(M, gel(vz,j), ell), gel(vy,j));
     974         497 :   return Flm_ker(Flm_Fl_sub(vy, kum->g, ell), ell);
     975             : }
     976             : 
     977             : /* J a vector of elements in nfz = relative extension of nf by polrel,
     978             :  * return the Steinitz element attached to the module generated by J */
     979             : static GEN
     980         672 : Stelt(GEN nf, GEN J, GEN polrel)
     981             : {
     982         672 :   long i, l = lg(J), vx = varn(polrel);
     983         672 :   GEN A = cgetg(l, t_VEC), I = cgetg(l, t_VEC);
     984        4550 :   for (i = 1; i < l; i++)
     985             :   {
     986        3878 :     GEN v = gel(J,i);
     987        3878 :     if (typ(v) == t_POL) { v = RgX_rem(v, polrel); setvarn(v,vx); }
     988        3878 :     gel(A,i) = v;
     989        3878 :     gel(I,i) = gen_1;
     990             :   }
     991         672 :   A = RgV_to_RgM(A, degpol(polrel));
     992         672 :   return idealprod(nf, gel(nfhnf(nf, mkvec2(A,I)),2));
     993             : }
     994             : 
     995             : static GEN
     996         133 : polrelKzK(toK_s *T, GEN x)
     997             : {
     998         133 :   GEN P = roots_to_pol(powtau(x, T->m, T->tau), 0);
     999         133 :   long i, l = lg(P);
    1000         133 :   for (i=2; i<l; i++) gel(P,i) = downtoK(T, gel(P,i));
    1001         133 :   return P;
    1002             : }
    1003             : 
    1004             : /* N: Cl_m(Kz) --> Cl_m(K), lift subgroup from bnr to bnrz using Algo 4.1.11 */
    1005             : static GEN
    1006         133 : invimsubgroup(GEN bnrz, GEN bnr, GEN subgroup, toK_s *T)
    1007             : {
    1008             :   long l, j;
    1009             :   GEN P, cyc, gen, U, polrel, StZk;
    1010         133 :   GEN nf = bnr_get_nf(bnr), nfz = bnr_get_nf(bnrz);
    1011         133 :   GEN polz = nf_get_pol(nfz), zkzD = nf_get_zkprimpart(nfz);
    1012             : 
    1013         133 :   polrel = polrelKzK(T, pol_x(varn(polz)));
    1014         133 :   StZk = Stelt(nf, zkzD, polrel);
    1015         133 :   cyc = bnr_get_cyc(bnrz); l = lg(cyc);
    1016         133 :   gen = bnr_get_gen(bnrz);
    1017         133 :   P = cgetg(l,t_MAT);
    1018         672 :   for (j=1; j<l; j++)
    1019             :   {
    1020         539 :     GEN g, id = idealhnf_shallow(nfz, gel(gen,j));
    1021         539 :     g = Stelt(nf, RgV_RgM_mul(zkzD, id), polrel);
    1022         539 :     g = idealdiv(nf, g, StZk); /* N_{Kz/K}(gen[j]) */
    1023         539 :     gel(P,j) = isprincipalray(bnr, g);
    1024             :   }
    1025         133 :   (void)ZM_hnfall_i(shallowconcat(P, subgroup), &U, 1);
    1026         133 :   setlg(U, l); for (j=1; j<l; j++) setlg(U[j], l);
    1027         133 :   return ZM_hnfmodid(U, cyc);
    1028             : }
    1029             : 
    1030             : static GEN
    1031         546 : pol_from_Newton(GEN S)
    1032             : {
    1033         546 :   long i, k, l = lg(S);
    1034         546 :   GEN C = cgetg(l+1, t_VEC), c = C + 1;
    1035         546 :   gel(c,0) = gen_1;
    1036         546 :   gel(c,1) = gel(S,1); /* gen_0 in our case */
    1037        1792 :   for (k = 2; k < l; k++)
    1038             :   {
    1039        1246 :     GEN s = gel(S,k);
    1040        1246 :     for (i = 2; i < k-1; i++) s = gadd(s, gmul(gel(S,i), gel(c,k-i)));
    1041        1246 :     gel(c,k) = gdivgs(s, -k);
    1042             :   }
    1043         546 :   return gtopoly(C, 0);
    1044             : }
    1045             : 
    1046             : /* - mu_b = sum_{0 <= i < m} floor(r_b r_{d-1-i} / ell) tau^i */
    1047             : static GEN
    1048        1218 : get_mmu(long b, GEN r, long ell)
    1049             : {
    1050        1218 :   long i, m = lg(r)-1;
    1051        1218 :   GEN M = cgetg(m+1, t_VEC);
    1052        1218 :   for (i = 0; i < m; i++) gel(M,i+1) = stoi((r[b + 1] * r[m - i]) / ell);
    1053        1218 :   return M;
    1054             : }
    1055             : 
    1056             : /* coeffs(x, a..b) in variable v >= varn(x) */
    1057             : static GEN
    1058       11116 : split_pol(GEN x, long v, long a, long b)
    1059             : {
    1060       11116 :   long i, l = degpol(x);
    1061       11116 :   GEN y = x + a, z;
    1062             : 
    1063       11116 :   if (l < b) b = l;
    1064       11116 :   if (a > b || varn(x) != v) return pol_0(v);
    1065        9856 :   l = b-a + 3;
    1066        9856 :   z = cgetg(l, t_POL); z[1] = x[1];
    1067        9856 :   for (i = 2; i < l; i++) gel(z,i) = gel(y,i);
    1068        9856 :   return normalizepol_lg(z, l);
    1069             : }
    1070             : 
    1071             : /* return (den_a * z) mod (v^ell - num_a/den_a), assuming deg(z) < 2*ell
    1072             :  * allow either num/den to be NULL (= 1) */
    1073             : static GEN
    1074        5558 : mod_Xell_a(GEN z, long v, long ell, GEN num_a, GEN den_a)
    1075             : {
    1076        5558 :   GEN z1 = split_pol(z, v, ell, degpol(z));
    1077        5558 :   GEN z0 = split_pol(z, v, 0,   ell-1); /* z = v^ell z1 + z0*/
    1078        5558 :   if (den_a) z0 = gmul(den_a, z0);
    1079        5558 :   if (num_a) z1 = gmul(num_a, z1);
    1080        5558 :   return gadd(z0, z1);
    1081             : }
    1082             : static GEN
    1083        1764 : to_alg(GEN nfz, GEN c, long v)
    1084             : {
    1085             :   GEN z, D;
    1086        1764 :   if (typ(c) != t_COL) return c;
    1087        1218 :   z = gmul(nf_get_zkprimpart(nfz), c);
    1088        1218 :   if (typ(z) == t_POL) setvarn(z, v);
    1089        1218 :   D = nf_get_zkden(nfz);
    1090        1218 :   if (!equali1(D)) z = RgX_Rg_div(z, D);
    1091        1218 :   return z;
    1092             : }
    1093             : 
    1094             : /* th. 5.3.5. and prop. 5.3.9. */
    1095             : static GEN
    1096         546 : compute_polrel(struct rnfkummer *kum, GEN be)
    1097             : {
    1098         546 :   toK_s *T = &kum->T;
    1099         546 :   long i, k, ell = kum->ell, m = T->m, vT = fetch_var(), vz = fetch_var();
    1100             :   GEN r, powtaubet, S, p1, root, num_t, den_t, nfzpol, powtau_prim_invbe;
    1101             :   GEN prim_Rk, C_Rk, prim_root, C_root, prim_invbe, C_invbe;
    1102         546 :   GEN nfz = bnf_get_nf(kum->bnfz);
    1103             :   pari_timer ti;
    1104             : 
    1105         546 :   r = cgetg(m+1,t_VECSMALL); /* r[i+1] = g^i mod ell */
    1106         546 :   r[1] = 1;
    1107         546 :   for (i=2; i<=m; i++) r[i] = (r[i-1] * kum->g) % ell;
    1108         546 :   powtaubet = powtau(be, m, T->tau);
    1109         546 :   if (DEBUGLEVEL>1) { err_printf("Computing Newton sums: "); timer_start(&ti); }
    1110         546 :   prim_invbe = Q_primitive_part(nfinv(nfz, be), &C_invbe);
    1111         546 :   powtau_prim_invbe = powtau(prim_invbe, m, T->tau);
    1112             : 
    1113         546 :   root = cgetg(ell + 2, t_POL);
    1114         546 :   root[1] = evalsigne(1) | evalvarn(0);
    1115         546 :   for (i = 0; i < ell; i++) gel(root,2+i) = gen_0;
    1116        1764 :   for (i = 0; i < m; i++)
    1117             :   { /* compute (1/be) ^ (-mu) instead of be^mu [mu << 0].
    1118             :      * 1/be = C_invbe * prim_invbe */
    1119        1218 :     GEN mmu = get_mmu(i, r, ell);
    1120             :     /* p1 = prim_invbe ^ -mu */
    1121        1218 :     p1 = to_alg(nfz, nffactorback(nfz, powtau_prim_invbe, mmu), vz);
    1122        1218 :     if (C_invbe) p1 = gmul(p1, powgi(C_invbe, RgV_sumpart(mmu, m)));
    1123             :     /* root += zeta_ell^{r_i} T^{r_i} be^mu_i */
    1124        1218 :     gel(root, 2 + r[i+1]) = monomial(p1, r[i+1], vT);
    1125             :   }
    1126             :   /* Other roots are as above with z_ell --> z_ell^j.
    1127             :    * Treat all contents (C_*) and principal parts (prim_*) separately */
    1128         546 :   prim_Rk = prim_root = Q_primitive_part(root, &C_root);
    1129         546 :   C_Rk = C_root;
    1130             : 
    1131         546 :   r = vecsmall_reverse(r); /* theta^ell = be^( sum tau^a r_{d-1-a} ) */
    1132             :   /* Compute modulo X^ell - 1, T^ell - t, nfzpol(vz) */
    1133         546 :   p1 = to_alg(nfz, nffactorback(nfz, powtaubet, r), vz);
    1134         546 :   num_t = Q_remove_denom(p1, &den_t);
    1135             : 
    1136         546 :   nfzpol = leafcopy(nf_get_pol(nfz));
    1137         546 :   setvarn(nfzpol, vz);
    1138         546 :   S = cgetg(ell+1, t_VEC); /* Newton sums */
    1139         546 :   gel(S,1) = gen_0;
    1140        1792 :   for (k = 2; k <= ell; k++)
    1141             :   { /* compute the k-th Newton sum */
    1142        1246 :     pari_sp av = avma;
    1143        1246 :     GEN z, D, Rk = gmul(prim_Rk, prim_root);
    1144        1246 :     C_Rk = mul_content(C_Rk, C_root);
    1145        1246 :     Rk = mod_Xell_a(Rk, 0, ell, NULL, NULL); /* mod X^ell - 1 */
    1146        5586 :     for (i = 2; i < lg(Rk); i++)
    1147             :     {
    1148        4340 :       if (typ(gel(Rk,i)) != t_POL) continue;
    1149        4312 :       z = mod_Xell_a(gel(Rk,i), vT, ell, num_t,den_t); /* mod T^ell - t */
    1150        4312 :       gel(Rk,i) = RgXQX_red(z, nfzpol); /* mod nfz.pol */
    1151             :     }
    1152        1246 :     if (den_t) C_Rk = mul_content(C_Rk, ginv(den_t));
    1153        1246 :     prim_Rk = Q_primitive_part(Rk, &D);
    1154        1246 :     C_Rk = mul_content(C_Rk, D); /* root^k = prim_Rk * C_Rk */
    1155             : 
    1156             :     /* Newton sum is ell * constant coeff (in X), which has degree 0 in T */
    1157        1246 :     z = polcoef_i(prim_Rk, 0, 0);
    1158        1246 :     z = polcoef_i(z      , 0,vT);
    1159        1246 :     z = downtoK(T, gmulgs(z, ell));
    1160        1246 :     if (C_Rk) z = gmul(z, C_Rk);
    1161        1246 :     gerepileall(av, C_Rk? 3: 2, &z, &prim_Rk, &C_Rk);
    1162        1246 :     if (DEBUGLEVEL>1) err_printf("%ld(%ld) ", k, timer_delay(&ti));
    1163        1246 :     gel(S,k) = z;
    1164             :   }
    1165         546 :   if (DEBUGLEVEL>1) err_printf("\n");
    1166         546 :   (void)delete_var();
    1167         546 :   (void)delete_var(); return pol_from_Newton(S);
    1168             : }
    1169             : 
    1170             : /* lift elt t in nf to nfz, algebraic form */
    1171             : static GEN
    1172         651 : lifttoKz(GEN nf, GEN t, compo_s *C)
    1173             : {
    1174         651 :   GEN x = nf_to_scalar_or_alg(nf, t);
    1175         651 :   if (typ(x) != t_POL) return x;
    1176         651 :   return RgX_RgXQ_eval(x, C->p, C->R);
    1177             : }
    1178             : /* lift ideal id in nf to nfz */
    1179             : static GEN
    1180         525 : ideallifttoKz(GEN nfz, GEN nf, GEN id, compo_s *C)
    1181             : {
    1182         525 :   GEN I = idealtwoelt(nf,id);
    1183         525 :   GEN x = nf_to_scalar_or_alg(nf, gel(I,2));
    1184         525 :   if (typ(x) != t_POL) return gel(I,1);
    1185         343 :   gel(I,2) = algtobasis(nfz, RgX_RgXQ_eval(x, C->p, C->R));
    1186         343 :   return idealhnf_two(nfz,I);
    1187             : }
    1188             : /* lift ideal pr in nf to ONE prime in nfz (the others are conjugate under tau
    1189             :  * and bring no further information on e_1 W). Assume pr coprime to
    1190             :  * index of both nf and nfz, and unramified in Kz/K (minor simplification) */
    1191             : static GEN
    1192         756 : prlifttoKz(GEN nfz, GEN nf, GEN pr, compo_s *C)
    1193             : {
    1194         756 :   GEN F, p = pr_get_p(pr), t = pr_get_gen(pr), T = nf_get_pol(nfz);
    1195         756 :   if (nf_get_degree(nf) != 1)
    1196             :   { /* restrict to primes above pr */
    1197         651 :     t = Q_primpart( lifttoKz(nf,t,C) );
    1198         651 :     T = FpX_gcd(FpX_red(T,p), FpX_red(t,p), p);
    1199         651 :     T = FpX_normalize(T, p);
    1200             :   }
    1201         756 :   F = FpX_factor(T, p);
    1202         756 :   return idealprimedec_kummer(nfz,gcoeff(F,1,1), pr_get_e(pr), p);
    1203             : }
    1204             : static GEN
    1205         511 : get_przlist(GEN L, GEN nfz, GEN nf, compo_s *C)
    1206             : {
    1207             :   long i, l;
    1208         511 :   GEN M = cgetg_copy(L, &l);
    1209         511 :   for (i = 1; i < l; i++) gel(M,i) = prlifttoKz(nfz, nf, gel(L,i), C);
    1210         511 :   return M;
    1211             : }
    1212             : 
    1213             : static void
    1214         497 : compositum_red(compo_s *C, GEN P, GEN Q)
    1215             : {
    1216         497 :   GEN p, q, a, z = gel(compositum2(P, Q),1);
    1217         497 :   a = gel(z,1);
    1218         497 :   p = gel(gel(z,2), 2);
    1219         497 :   q = gel(gel(z,3), 2);
    1220         497 :   C->k = itos( gel(z,4) );
    1221             :   /* reduce R. FIXME: should be polredbest(a, 1), but breaks rnfkummer bench */
    1222         497 :   z = polredabs0(a, nf_ORIG|nf_PARTIALFACT);
    1223         497 :   C->R = gel(z,1);
    1224         497 :   a = gel(gel(z,2), 2);
    1225         497 :   C->p = RgX_RgXQ_eval(p, a, C->R);
    1226         497 :   C->q = RgX_RgXQ_eval(q, a, C->R);
    1227         497 :   C->rev = QXQ_reverse(a, C->R);
    1228         497 : }
    1229             : 
    1230             : /* replace P->C^(-deg P) P(xC) for the largest integer C such that coefficients
    1231             :  * remain algebraic integers. Lift *rational* coefficients */
    1232             : static void
    1233         546 : nfX_Z_normalize(GEN nf, GEN P)
    1234             : {
    1235             :   long i, l;
    1236         546 :   GEN C, Cj, PZ = cgetg_copy(P, &l);
    1237         546 :   PZ[1] = P[1];
    1238        2884 :   for (i = 2; i < l; i++) /* minor variation on RgX_to_nfX (create PZ) */
    1239             :   {
    1240        2338 :     GEN z = nf_to_scalar_or_basis(nf, gel(P,i));
    1241        2338 :     if (typ(z) == t_INT)
    1242        1666 :       gel(PZ,i) = gel(P,i) = z;
    1243             :     else
    1244         672 :       gel(PZ,i) = ZV_content(z);
    1245             :   }
    1246         546 :   (void)ZX_Z_normalize(PZ, &C);
    1247             : 
    1248         546 :   if (C == gen_1) return;
    1249         119 :   Cj = C;
    1250         504 :   for (i = l-2; i > 1; i--)
    1251             :   {
    1252         385 :     if (i != l-2) Cj = mulii(Cj, C);
    1253         385 :     gel(P,i) = gdiv(gel(P,i), Cj);
    1254             :   }
    1255             : }
    1256             : 
    1257             : /* set kum->vecC, kum->Q */
    1258             : static void
    1259         497 : _rnfkummer_step4(struct rnfkummer *kum, GEN cycgen, long d, long m)
    1260             : {
    1261         497 :   long i, j, rc = kum->rc;
    1262         497 :   GEN vecC, vecB = cgetg(rc+1,t_VEC), Tc = cgetg(rc+1,t_MAT);
    1263         497 :   GEN gen = bnf_get_gen(kum->bnfz), u = kum->u;
    1264         497 :   ulong ell = kum->ell;
    1265         861 :   for (j=1; j<=rc; j++)
    1266             :   {
    1267         364 :     GEN p1 = tauofideal(gel(gen,j), &kum->tau);
    1268         364 :     isprincipalell(kum->bnfz, p1, cycgen,u,ell,rc, &gel(Tc,j), &gel(vecB,j));
    1269             :   }
    1270             : 
    1271         497 :   kum->vecC = vecC = const_vec(rc, trivial_fact());
    1272         497 :   if (rc)
    1273             :   {
    1274         280 :     GEN p1 = Flm_powers(Tc, m-2, ell), p2 = vecB;
    1275         630 :     for (j = 1; j < m; j++)
    1276             :     {
    1277         350 :       GEN z = Flm_Fl_mul(gel(p1,m-j), Fl_mul(j,d,ell), ell);
    1278         350 :       p2 = tauofvec(p2, &kum->tau);
    1279         798 :       for (i = 1; i <= rc; i++)
    1280         896 :         gel(vecC,i) = famat_mul_shallow(gel(vecC,i),
    1281         448 :                                         famatV_zv_factorback(p2, gel(z,i)));
    1282             :     }
    1283         280 :     for (i = 1; i <= rc; i++) gel(vecC,i) = famat_reduce(gel(vecC,i));
    1284             :   }
    1285         497 :   kum->Q = Flm_ker(Flm_Fl_sub(Flm_transpose(Tc), kum->g, ell), ell);
    1286         497 : }
    1287             : 
    1288             : static GEN
    1289         497 : _rnfkummer_step5(struct rnfkummer *kum, GEN vselmer)
    1290             : {
    1291         497 :   GEN W = kervirtualunit(kum, vselmer);
    1292         497 :   long j, l = lg(W);
    1293        1463 :   for (j = 1; j < l; j++)
    1294         966 :     gel(W,j) = famat_reduce(famatV_zv_factorback(vselmer, gel(W,j)));
    1295         497 :   settyp(W, t_VEC); return W;
    1296             : }
    1297             : 
    1298             : static GEN
    1299         525 : _rnfkummer_step18(struct rnfkummer *kum, GEN bnr, GEN subgroup, GEN M,
    1300             :      GEN vecWB, GEN vecMsup, long all)
    1301             : {
    1302         525 :   ulong ell = kum->ell;
    1303         525 :   long i, dK, ncyc = 0, firstpass = all < 0, rk = 0, lW = lg(kum->vecW);
    1304         525 :   GEN K, y, nf = bnr_get_nf(bnr), gell = utoipos(ell), res = NULL, mat = NULL;
    1305             : 
    1306         525 :   K = Flm_ker(M, ell);
    1307         525 :   if (all < 0) K = fix_kernel(K, M, vecMsup, lW, ell);
    1308         525 :   if (DEBUGLEVEL>2) err_printf("Step 18\n");
    1309         525 :   dK = lg(K)-1;
    1310         525 :   y = cgetg(dK+1,t_VECSMALL);
    1311         525 :   if (all) res = cgetg(1, t_VEC);
    1312         525 :   if (all < 0) { ncyc = dK; rk = 0; mat = zero_Flm(lg(M)-1, ncyc); }
    1313             : 
    1314             :   do {
    1315         525 :     dK = lg(K)-1;
    1316        1071 :     while (dK)
    1317             :     {
    1318         532 :       for (i=1; i<dK; i++) y[i] = 0;
    1319         532 :       y[i] = 1; /* y = [0,...,0,1,0,...,0], 1 at dK'th position */
    1320             :       do
    1321             :       { /* cf. algo 5.3.18 */
    1322         546 :         GEN H, be, P, X = Flm_Flc_mul(K, y, ell);
    1323         546 :         if (ok_congruence(X, ell, lW, vecMsup))
    1324             :         {
    1325         546 :           pari_sp av = avma;
    1326         546 :           if (all < 0)
    1327             :           {
    1328           0 :             gel(mat, rk+1) = X;
    1329           0 :             if (Flm_rank(mat,ell) <= rk) continue;
    1330           0 :             rk++;
    1331             :           }
    1332         546 :           be = compute_beta(X, vecWB, gell, kum->bnfz);
    1333         546 :           P = compute_polrel(kum, be);
    1334         546 :           nfX_Z_normalize(nf, P);
    1335         546 :           if (DEBUGLEVEL>1) err_printf("polrel(beta) = %Ps\n", P);
    1336         546 :           if (!all) {
    1337         511 :             H = rnfnormgroup(bnr, P);
    1338         511 :             if (ZM_equal(subgroup, H)) return P; /* DONE */
    1339           0 :             set_avma(av); continue;
    1340             :           } else {
    1341          35 :             GEN P0 = Q_primpart(lift_shallow(P));
    1342          35 :             GEN g = nfgcd(P0, RgX_deriv(P0), nf_get_pol(nf), nf_get_index(nf));
    1343          35 :             if (degpol(g)) continue;
    1344          35 :             H = rnfnormgroup(bnr, P);
    1345          35 :             if (!ZM_equal(subgroup,H) && !bnrisconductor(bnr,H)) continue;
    1346             :           }
    1347          35 :           P = gerepilecopy(av, P);
    1348          35 :           res = shallowconcat(res, P);
    1349          35 :           if (all < 0 && rk == ncyc) return res;
    1350          35 :           if (firstpass) break;
    1351             :         }
    1352          35 :       } while (increment(y, dK, ell));
    1353          21 :       y[dK--] = 0;
    1354             :     }
    1355          14 :   } while (firstpass--);
    1356          14 :   if (!res) pari_err_BUG("kummer [no solution]");
    1357          14 :   return res;
    1358             : }
    1359             : 
    1360             : /* alg 5.3.5 */
    1361             : static void
    1362         497 : rnfkummer_init(struct rnfkummer *kum, GEN bnf, ulong ell, long prec)
    1363             : {
    1364         497 :   compo_s *COMPO = &kum->COMPO;
    1365         497 :   toK_s *T = &kum->T;
    1366         497 :   GEN nf  = bnf_get_nf(bnf), polnf = nf_get_pol(nf), vselmer, bnfz, cyc, cycgen;
    1367             :   long degK, degKz, m, d;
    1368             :   ulong g;
    1369             :   pari_timer ti;
    1370         497 :   if (DEBUGLEVEL>2) err_printf("Step 1\n");
    1371         497 :   compositum_red(COMPO, polnf, polcyclo(ell, varn(polnf)));
    1372         497 :   if (DEBUGLEVEL)
    1373             :   {
    1374           0 :     timer_printf(&ti, "[rnfkummer] compositum");
    1375           0 :     if (DEBUGLEVEL>1) err_printf("polred(compositum) = %Ps\n",COMPO->R);
    1376             :   }
    1377         497 :   if (DEBUGLEVEL>2) err_printf("Step 2\n");
    1378         497 :   degK  = degpol(polnf);
    1379         497 :   degKz = degpol(COMPO->R);
    1380         497 :   m = degKz / degK;
    1381         497 :   d = (ell-1) / m;
    1382         497 :   g = Fl_powu(pgener_Fl(ell), d, ell);
    1383         497 :   if (Fl_powu(g, m, ell*ell) == 1) g += ell;
    1384             :   /* ord(g) = m in all (Z/ell^k)^* */
    1385         497 :   if (DEBUGLEVEL>2) err_printf("Step 3\n");
    1386             :   /* could factor disc(R) using th. 2.1.6. */
    1387         497 :   kum->bnfz = bnfz = Buchall(COMPO->R, nf_FORCE, maxss(prec,BIGDEFAULTPREC));
    1388         497 :   if (DEBUGLEVEL) timer_printf(&ti, "[rnfkummer] bnfinit(Kz)");
    1389         497 :   cycgen = bnf_build_cycgen(bnfz);
    1390         497 :   cyc = bnf_get_cyc(bnfz);
    1391         497 :   kum->ell = ell;
    1392         497 :   kum->rc = prank(cyc, ell);
    1393         497 :   kum->u = get_u(cyc, kum->rc, ell);
    1394         497 :   kum->g = g;
    1395             : 
    1396         497 :   vselmer = get_Selmer(bnfz, cycgen, kum->rc);
    1397         497 :   get_tau(kum);
    1398         497 :   if (DEBUGLEVEL>2) err_printf("Step 4\n");
    1399         497 :   _rnfkummer_step4(kum, cycgen, d, m);
    1400         497 :   if (DEBUGLEVEL>2) err_printf("Step 5\n");
    1401         497 :   kum->vecW = _rnfkummer_step5(kum, vselmer);
    1402         497 :   if (DEBUGLEVEL>2) err_printf("Step 8\n");
    1403             :   /* left inverse */
    1404         497 :   T->invexpoteta1 = RgM_inv(RgXQ_matrix_pow(COMPO->p, degKz, degK, COMPO->R));
    1405         497 :   T->polnf = polnf;
    1406         497 :   T->tau = &kum->tau;
    1407         497 :   T->m = m;
    1408         497 :   T->powg = Fl_powers(g, m, ell);
    1409         497 : }
    1410             : 
    1411             : static GEN
    1412         525 : rnfkummer_ell(struct rnfkummer *kum, GEN bnr, GEN subgroup, long all)
    1413             : {
    1414         525 :   ulong ell = kum->ell, mginv;
    1415         525 :   GEN bnfz = kum->bnfz, nfz = bnf_get_nf(bnfz), cycgen = bnf_build_cycgen(bnfz);
    1416         525 :   GEN nf = bnr_get_nf(bnr), bid = bnr_get_bid(bnr), ideal = bid_get_ideal(bid);
    1417         525 :   GEN vecC = kum->vecC, vecW = kum->vecW, u = kum->u, Q = kum->Q;
    1418         525 :   long lW = lg(vecW), rc = kum->rc, i, j, lSml2, lSp, lSl2, dc;
    1419         525 :   toK_s *T = &kum->T;
    1420             :   primlist L;
    1421             :   GEN gothf, idealz, Sp, prSp, vecAp, vecBp, matP, vecWA, vecWB, vecMsup, M;
    1422             : 
    1423         525 :   idealz = ideallifttoKz(nfz, nf, ideal, &kum->COMPO);
    1424         525 :   if (umodiu(gcoeff(ideal,1,1), ell)) gothf = idealz;
    1425             :   else
    1426             :   { /* ell | N(ideal) */
    1427         133 :     GEN bnrz = Buchray(bnfz, idealz, nf_INIT|nf_GEN);
    1428         133 :     GEN subgroupz = invimsubgroup(bnrz, bnr, subgroup, T);
    1429         133 :     gothf = bnrconductor_i(bnrz,subgroupz,0);
    1430             :   }
    1431         525 :   if (DEBUGLEVEL>2) err_printf("Step 9, 10 and 11\n");
    1432         525 :   i = build_list_Hecke(&L, nfz, NULL, gothf, ell, T->tau);
    1433         525 :   if (i) return no_sol(all,i);
    1434             : 
    1435         525 :   lSml2 = lg(L.Sml2);
    1436         525 :   Sp = shallowconcat(L.Sm, L.Sml1); lSp = lg(Sp);
    1437         525 :   prSp = shallowconcat(L.Sml2, L.Sl); lSl2 = lg(prSp);
    1438             : 
    1439         525 :   if (DEBUGLEVEL>2) err_printf("Step 12\n");
    1440         525 :   vecAp = cgetg(lSp, t_VEC);
    1441         525 :   vecBp = cgetg(lSp, t_VEC);
    1442         525 :   matP  = cgetg(lSp, t_MAT);
    1443         812 :   for (j = 1; j < lSp; j++)
    1444             :   {
    1445             :     GEN e, a;
    1446         287 :     isprincipalell(bnfz, gel(Sp,j), cycgen,u,ell,rc, &e, &a);
    1447         287 :     gel(matP,j) = e;
    1448         287 :     gel(vecBp,j) = famat_mul_shallow(famatV_zv_factorback(vecC, zv_neg(e)), a);
    1449         287 :     gel(vecAp,j) = lambdaofelt(gel(vecBp,j), T);
    1450             :   }
    1451         525 :   if (DEBUGLEVEL>2) err_printf("Step 13\n");
    1452         525 :   vecWA = shallowconcat(vecW, vecAp);
    1453         525 :   vecWB = shallowconcat(vecW, vecBp);
    1454             : 
    1455         525 :   if (DEBUGLEVEL>2) err_printf("Step 14, 15 and 17\n");
    1456         525 :   mginv = Fl_div(T->m, kum->g, ell);
    1457         525 :   vecMsup = cgetg(lSml2,t_VEC);
    1458         525 :   M = NULL;
    1459        1232 :   for (i = 1; i < lSl2; i++)
    1460             :   {
    1461         707 :     GEN pr = gel(prSp,i);
    1462         707 :     long e = pr_get_e(pr), z = ell * (e / (ell-1));
    1463             : 
    1464         707 :     if (i < lSml2)
    1465             :     {
    1466         140 :       z += 1 - L.ESml2[i];
    1467         140 :       gel(vecMsup,i) = logall(nfz, vecWA,lW,mginv,ell, pr,z+1);
    1468             :     }
    1469         707 :     M = vconcat(M, logall(nfz, vecWA,lW,mginv,ell, pr,z));
    1470             :   }
    1471         525 :   dc = lg(Q)-1;
    1472         525 :   if (dc)
    1473             :   {
    1474         168 :     GEN QtP = Flm_mul(Flm_transpose(Q), matP, ell);
    1475         168 :     M = vconcat(M, shallowconcat(zero_Flm(dc,lW-1), QtP));
    1476             :   }
    1477         525 :   if (!M) M = zero_Flm(1, lSp-1 + lW-1);
    1478             : 
    1479         525 :   if (!all)
    1480             :   { /* primes landing in subgroup must be totally split */
    1481         511 :     GEN lambdaWB = shallowconcat(lambdaofvec(vecW, T), vecAp);/*vecWB^lambda*/
    1482         511 :     GEN Lpr = get_prlist(bnr, subgroup, ell, bnfz);
    1483         511 :     GEN Lprz= get_przlist(Lpr, nfz, nf, &kum->COMPO);
    1484         511 :     M = vconcat(M, subgroup_info(bnfz, Lprz, ell, lambdaWB));
    1485             :   }
    1486         525 :   if (DEBUGLEVEL>2) err_printf("Step 16\n");
    1487         525 :   return _rnfkummer_step18(kum,bnr,subgroup, M, vecWB, vecMsup, all);
    1488             : }
    1489             : 
    1490             : static void
    1491        1330 : bnrclassfield_sanitize(GEN *pbnr, GEN *pH)
    1492             : {
    1493        1330 :   GEN cyc, cnd, bnr = *pbnr, H = *pH, T;
    1494             : 
    1495        1330 :   if (nftyp(bnr)==typ_BNF) bnr = Buchray(bnr, gen_1, nf_INIT);
    1496        1281 :   else checkbnr(bnr);
    1497        1316 :   T = nf_get_pol(bnr_get_nf(bnr));
    1498        1316 :   if (!varn(T)) pari_err_PRIORITY("bnrclassfield", T, "=", 0);
    1499        1302 :   cyc = bnr_get_cyc(bnr);
    1500        1302 :   if (!H) H = gen_0;
    1501        1302 :   switch(typ(H))
    1502             :   {
    1503         301 :     case t_INT: H = scalarmat_shallow(H, lg(cyc)-1);
    1504        1295 :     case t_MAT: H = hnfmodid(H, cyc); break;
    1505           7 :     default: pari_err_TYPE("bnrclassfield [subgroup]", H);
    1506             :   }
    1507        1281 :   cnd = bnrconductor_i(bnr, H, 2);
    1508        1281 :   *pbnr = gel(cnd,2);
    1509        1281 :   *pH = gel(cnd,3);
    1510        1281 : }
    1511             : 
    1512             : static GEN
    1513         994 : _rnfkummer(GEN bnr, GEN subgroup, long all, long prec)
    1514             : {
    1515             :   ulong ell;
    1516             :   GEN gell;
    1517             :   struct rnfkummer kum;
    1518             : 
    1519         994 :   bnrclassfield_sanitize(&bnr, &subgroup);
    1520         987 :   gell = get_gell(bnr,subgroup,all);
    1521         987 :   if (typ(gell) != t_INT) pari_err_TYPE("rnfkummer",gell);
    1522         987 :   ell = itou(gell);
    1523         987 :   if (ell == 1) return pol_x(0);
    1524         987 :   if (!uisprime(ell)) pari_err_IMPL("rnfkummer for composite relative degree");
    1525         987 :   if (all && all != -1 && umodiu(bnr_get_no(bnr), ell))
    1526           7 :     return cgetg(1, t_VEC);
    1527         980 :   if (bnf_get_tuN(bnr_get_bnf(bnr)) % ell == 0)
    1528         616 :     return rnfkummersimple(bnr, subgroup, ell, all);
    1529         364 :   rnfkummer_init(&kum, bnr_get_bnf(bnr), ell, prec);
    1530         364 :   return rnfkummer_ell(&kum, bnr, subgroup, all == -1? 0: all);
    1531             : }
    1532             : 
    1533             : GEN
    1534         784 : rnfkummer(GEN bnr, GEN subgroup, long all, long prec)
    1535             : {
    1536         784 :   pari_sp av = avma;
    1537         784 :   return gerepilecopy(av, _rnfkummer(bnr, subgroup, all, prec));
    1538             : }
    1539             : 
    1540             : /*******************************************************************/
    1541             : /*                        bnrclassfield                            */
    1542             : /*******************************************************************/
    1543             : 
    1544             : /* TODO: could be exported */
    1545             : static void
    1546      199038 : gsetvarn(GEN x, long v)
    1547             : {
    1548             :   long i;
    1549      199038 :   switch(typ(x))
    1550             :   {
    1551             :     case t_POL: case t_SER:
    1552        1477 :       setvarn(x, v); return;
    1553             :     case t_LIST:
    1554           0 :       x = list_data(x); if (!x) return;
    1555             :       /* fall through t_VEC */
    1556             :     case t_VEC: case t_COL: case t_MAT:
    1557       23275 :       for (i = lg(x)-1; i > 0; i--) gsetvarn(gel(x,i), v);
    1558             :   }
    1559             : }
    1560             : 
    1561             : /* emb root of pol as polmod modulo pol2, return relative polynomial */
    1562             : static GEN
    1563         189 : relative_pol(GEN pol, GEN emb, GEN pol2)
    1564             : {
    1565             :   GEN eqn, polrel;
    1566         189 :   if (degree(pol)==1) return pol2;
    1567         147 :   emb = liftpol(emb);
    1568         147 :   eqn = gsub(emb, pol_x(varn(pol)));
    1569         147 :   eqn = Q_remove_denom(eqn, NULL);
    1570         147 :   polrel = nfgcd(pol2, eqn, pol, NULL);
    1571         147 :   return RgX_Rg_div(polrel, leading_coeff(polrel));
    1572             : }
    1573             : 
    1574             : /* pol defines K/nf */
    1575             : static GEN
    1576         210 : bnrclassfield_tower(GEN bnr, GEN subgroup, GEN TB, GEN p, long finaldeg, long absolute, long prec)
    1577             : {
    1578         210 :   pari_sp av = avma;
    1579             :   GEN nf, nf2, rnf, bnf, bnf2, bnr2, q, H, dec, cyc, pk, sgpk, pol2, emb, emb2, famod, fa, Lbad, cert;
    1580             :   long i, r1, ell, sp, spk, last;
    1581             :   forprime_t iter;
    1582             : 
    1583         210 :   bnf = bnr_get_bnf(bnr);
    1584         210 :   nf = bnf_get_nf(bnf);
    1585         210 :   rnf = rnfinit0(nf, TB, 1);
    1586         210 :   nf2 = rnf_build_nfabs(rnf, prec);
    1587         210 :   gsetvarn(nf2, varn(nf_get_pol(nf)));
    1588             : 
    1589         210 :   cert = nfcertify(nf2);
    1590         210 :   if (lg(cert)>1)
    1591             :   {
    1592           0 :     rnf = rnfinit0(nf, gel(TB,1), 1);
    1593           0 :     nf2 = rnf_build_nfabs(rnf, prec);
    1594           0 :     gsetvarn(nf2, varn(nf_get_pol(nf)));
    1595             :   }
    1596             : 
    1597         210 :   r1 = nf_get_r1(nf2);
    1598         210 :   bnf2 = Buchall(nf2, nf_FORCE, prec);
    1599             : 
    1600         210 :   sp = itos(p);
    1601         210 :   spk = sp * rnf_get_degree(rnf);
    1602         210 :   pk = stoi(spk);
    1603         210 :   sgpk = hnfmodid(subgroup,pk);
    1604         210 :   last = spk==finaldeg;
    1605             : 
    1606             :   /* compute conductor */
    1607         210 :   famod = gel(bid_get_fact2(bnr_get_bid(bnr)),1);
    1608         210 :   if (lg(famod)==1)
    1609             :   {
    1610         140 :     fa = trivial_fact();
    1611         140 :     Lbad = cgetg(1, t_VECSMALL);
    1612             :   }
    1613             :   else
    1614             :   {
    1615          70 :     long j=1;
    1616          70 :     fa = cgetg(3, t_MAT);
    1617          70 :     gel(fa,1) = cgetg(lg(famod), t_VEC);
    1618          70 :     Lbad = cgetg(lg(famod), t_VEC);
    1619         175 :     for(i=1; i<lg(famod); i++)
    1620             :     {
    1621         105 :       GEN pr = gel(famod,i);
    1622         105 :       gmael(fa,1,i) = rnfidealprimedec(rnf, pr);
    1623         105 :       q = pr_get_p(pr);
    1624         105 :       if (lgefint(q) == 3) gel(Lbad,j++) = q;
    1625             :     }
    1626          70 :     setlg(Lbad,j);
    1627          70 :     Lbad = ZV_to_zv(ZV_sort_uniq(Lbad));
    1628          70 :     gel(fa,1) = shallowconcat1(gel(fa,1));
    1629          70 :     settyp(gel(fa,1), t_COL);
    1630          70 :     gel(fa,2) = cgetg(lg(gel(fa,1)), t_COL);
    1631         189 :     for (i=1; i<lg(gel(fa,1)); i++)
    1632             :     {
    1633         119 :       GEN pr = gcoeff(fa,i,1);
    1634         119 :       long e = equalii(p, pr_get_p(pr))? 1 + (pr_get_e(pr)*sp) / (sp-1): 1;
    1635         119 :       gcoeff(fa,i,2) = utoipos(e);
    1636             :     }
    1637             :   }
    1638         210 :   bnr2 = Buchray(bnf2, mkvec2(fa, const_vec(r1,gen_1)), nf_INIT);
    1639             : 
    1640             :   /* compute subgroup */
    1641         210 :   cyc = bnr_get_cyc(bnr2);
    1642         210 :   H = Flm_image(zv_diagonal(ZV_to_Flv(cyc,sp)), sp);
    1643         210 :   u_forprime_init(&iter, 2, ULONG_MAX);
    1644       12768 :   while ((ell = u_forprime_next(&iter))) if (!zv_search(Lbad, ell))
    1645             :   {
    1646       12460 :     dec = idealprimedec_limit_f(nf, utoi(ell), 1);
    1647       24444 :     for (i=1; i<lg(dec); i++)
    1648             :     {
    1649       11984 :       GEN pr = gel(dec,i), Pr = gel(rnfidealprimedec(rnf, pr), 1);
    1650       11984 :       long f = pr_get_f(Pr) / pr_get_f(pr);
    1651       11984 :       GEN vpr = FpC_Fp_mul(isprincipalray(bnr, pr), utoi(f), pk);
    1652       11984 :       if (gequal0(ZC_hnfrem(vpr,sgpk)))
    1653        1225 :         H = vec_append(H, ZV_to_Flv(isprincipalray(bnr2, Pr), sp));
    1654             :     }
    1655       12460 :     if (lg(H) > lg(cyc)+3)
    1656             :     {
    1657         210 :       H = Flm_image(H, sp);
    1658         210 :       if (lg(cyc)-lg(H) == 1) break;
    1659             :     }
    1660             :   }
    1661         210 :   H = hnfmodid(shallowconcat(zm_to_ZM(H), diagonal_shallow(cyc)), p);
    1662             : 
    1663             :   /* polynomial over nf2 */
    1664         210 :   pol2 = _rnfkummer(bnr2, H, 0, prec);
    1665             :   /* absolute polynomial */
    1666         210 :   pol2 = rnfequation2(nf2, pol2);
    1667         210 :   emb2 = gel(pol2,2); /* generator of nf2 as polmod modulo pol2 */
    1668         210 :   pol2 = gel(pol2,1);
    1669             :   /* polynomial over nf */
    1670         210 :   if (!absolute || !last)
    1671             :   {
    1672         189 :     emb = rnf_get_alpha(rnf); /* generator of nf as polynomial in nf2 */
    1673         189 :     emb = poleval(emb, emb2); /* generator of nf as polmod modulo pol2 */
    1674         189 :     pol2 = relative_pol(nf_get_pol(nf), emb, pol2);
    1675             :   }
    1676         210 :   if (!last) pol2 = rnfpolredbest(nf, pol2, 0);
    1677             : 
    1678         210 :   obj_free(rnf);
    1679         210 :   pol2 = gerepilecopy(av, pol2);
    1680         210 :   if (last) return pol2;
    1681          56 :   TB = mkvec2(pol2, gel(TB,2));
    1682          56 :   return bnrclassfield_tower(bnr, subgroup, TB, p, finaldeg, absolute, prec);
    1683             : }
    1684             : 
    1685             : /* subgroups H_i of bnr s.t. bnr/H_i is cyclic and inter_i H_i = subgroup */
    1686             : static GEN
    1687         364 : cyclic_compos(GEN subgroup)
    1688             : {
    1689         364 :   pari_sp av = avma;
    1690         364 :   GEN U, L, pe, D = smithclean( ZM_snfall_i(subgroup, &U, NULL, 1) );
    1691         364 :   long i, l = lg(D);
    1692             : 
    1693         364 :   L = cgetg(l, t_VEC);
    1694         364 :   if (l == 1) return L;
    1695         364 :   pe = gel(D,1); U = matinvmod(U, pe);
    1696         903 :   for (i = 1; i < l; i++)
    1697         539 :     gel(L,i) = hnfmodid(shallowconcat(subgroup, vecsplice(U,i)),pe);
    1698         364 :   return gerepilecopy(av, L);
    1699             : }
    1700             : 
    1701             : /* set kum=NULL if roots of unity already in base field */
    1702             : /* absolute=1 allowed if extension is cyclic with exponent>1 */
    1703             : static GEN
    1704         364 : bnrclassfield_primepower(struct rnfkummer *ptkum, GEN bnr, GEN subgroup, GEN p,
    1705             :   GEN P, long absolute, long prec)
    1706             : {
    1707         364 :   GEN res, subs = cyclic_compos(subgroup);
    1708         364 :   long i, l = lg(subs);
    1709             : 
    1710         364 :   res = cgetg(l,t_VEC);
    1711         903 :   for (i = 1; i < l; i++)
    1712             :   {
    1713         539 :     GEN H = gel(subs,i), cnd = bnrconductor_i(bnr, hnfmodid(H,p), 2);
    1714         539 :     GEN pol, pe, bnr2 = gel(cnd,2), Hp = gel(cnd,3);
    1715         539 :     if (ptkum)  pol = rnfkummer_ell(ptkum, bnr2, Hp, 0);
    1716         378 :     else        pol = rnfkummersimple(bnr2, Hp, itos(p), 0);
    1717         539 :     pe = ZM_det_triangular(H);
    1718         539 :     if (!equalii(p,pe))
    1719         154 :       pol = bnrclassfield_tower(bnr, H, mkvec2(pol,P), p, itos(pe), absolute, prec);
    1720         539 :     gel(res,i) = pol;
    1721             :   }
    1722         364 :   return res;
    1723             : }
    1724             : 
    1725             : /* partition of v into two subsets whose products are as balanced as possible */
    1726             : /* assume v sorted */
    1727             : static GEN
    1728          98 : vecsmall_balance(GEN v)
    1729             : {
    1730             :   forvec_t T;
    1731          98 :   GEN xbounds, x, vuniq, mult, ind, prod, prodbest = gen_0, bound,
    1732          98 :       xbest = NULL, res1, res2;
    1733          98 :   long i=1, j, k1, k2;
    1734          98 :   if (lg(v) == 3) return mkvec2(mkvecsmall(1), mkvecsmall(2));
    1735          35 :   vuniq = cgetg(lg(v), t_VECSMALL);
    1736          35 :   mult = cgetg(lg(v), t_VECSMALL);
    1737          35 :   ind = cgetg(lg(v), t_VECSMALL);
    1738          35 :   vuniq[1] = v[1];
    1739          35 :   mult[1] = 1;
    1740          35 :   ind[1] = 1;
    1741         140 :   for (j=2; j<lg(v); j++)
    1742             :   {
    1743         105 :     if (v[j] == vuniq[i]) mult[i]++;
    1744             :     else
    1745             :     {
    1746          14 :       i++;
    1747          14 :       vuniq[i] = v[j];
    1748          14 :       mult[i] = 1;
    1749          14 :       ind[i] = j;
    1750             :     }
    1751             :   }
    1752          35 :   setlg(vuniq, ++i);
    1753          35 :   setlg(mult, i);
    1754          35 :   setlg(ind, i);
    1755             : 
    1756          35 :   vuniq = zv_to_ZV(vuniq);
    1757          35 :   prod = factorback2(vuniq, mult);
    1758          35 :   bound = sqrti(prod);
    1759          35 :   xbounds = cgetg(lg(mult), t_VEC);
    1760          35 :   for (i=1; i<lg(mult); i++) gel(xbounds,i) = mkvec2s(0,mult[i]);
    1761             : 
    1762          35 :   forvec_init(&T, xbounds, 0);
    1763         273 :   while ((x = forvec_next(&T)))
    1764             :   {
    1765         203 :     prod = factorback2(vuniq, x);
    1766         203 :     if (cmpii(prod,bound)<=0 && cmpii(prod,prodbest)>0)
    1767             :     {
    1768          91 :       prodbest = prod;
    1769          91 :       xbest = gcopy(x);
    1770             :     }
    1771             :   }
    1772          35 :   res1 = cgetg(lg(v), t_VECSMALL);
    1773          35 :   res2 = cgetg(lg(v), t_VECSMALL);
    1774          84 :   for (i=1,k1=1,k2=1; i<lg(xbest); i++)
    1775             :   {
    1776          49 :     for (j=0; j<itos(gel(xbest,i)); j++) res1[k1++] = ind[i]+j;
    1777          49 :     for (; j<mult[i]; j++)               res2[k2++] = ind[i]+j;
    1778             :   }
    1779          35 :   setlg(res1, k1);
    1780          35 :   setlg(res2, k2); return mkvec2(res1, res2);
    1781             : }
    1782             : 
    1783             : /* TODO nfcompositum should accept vectors of pols */
    1784             : /* assume all fields are linearly disjoint */
    1785             : /* assume the polynomials are sorted by degree */
    1786             : static GEN
    1787         280 : nfcompositumall(GEN nf, GEN L)
    1788             : {
    1789             :   GEN pol, vdeg, part;
    1790             :   long i;
    1791         280 :   if (lg(L)==2) return gel(L,1);
    1792          98 :   vdeg = cgetg(lg(L), t_VECSMALL);
    1793          98 :   for (i=1; i<lg(L); i++) vdeg[i] = degree(gel(L,i));
    1794          98 :   part = vecsmall_balance(vdeg);
    1795          98 :   pol = cgetg(3, t_VEC);
    1796         294 :   for (i = 1; i < 3; i++)
    1797             :   {
    1798         196 :     GEN L2 = vecpermute(L, gel(part,i)), T = nfcompositumall(nf, L2);
    1799         196 :     gel(pol,i) = rnfpolredbest(nf, T, 0);
    1800             :   }
    1801          98 :   return nfcompositum(nf, gel(pol,1), gel(pol,2), 2);
    1802             : }
    1803             : 
    1804             : /* flag:
    1805             :  * 0 list of polynomials whose compositum is the extension
    1806             :  * 1 single polynomial
    1807             :  * 2 single absolute polynomial */
    1808             : GEN
    1809         350 : bnrclassfield(GEN bnr, GEN subgroup, long flag, long prec)
    1810             : {
    1811         350 :   pari_sp av = avma;
    1812             :   GEN N, fa, res, bnf, nf, P, PN, Pmod, EN;
    1813             :   long i, absolute, lPN;
    1814             :   struct rnfkummer kum;
    1815         350 :   if (flag<0 || flag>2) pari_err_FLAG("bnrclassfield [must be 0,1 or 2]");
    1816         336 :   bnrclassfield_sanitize(&bnr, &subgroup);
    1817             : 
    1818         294 :   N = ZM_det_triangular(subgroup);
    1819         294 :   if (equali1(N)) { set_avma(av); return pol_x(0); }
    1820         287 :   fa = Z_factor(N);
    1821         287 :   PN = gel(fa,1); lPN = lg(PN);
    1822         287 :   EN = gel(fa,2);
    1823         287 :   if (lgefint(gel(PN,lPN-1)) > 3)
    1824           7 :     pari_err_OVERFLOW("bnrclassfield [extension of too large degree]");
    1825         280 :   bnf = bnr_get_bnf(bnr);
    1826         280 :   nf = bnf_get_nf(bnf);
    1827             : 
    1828             :   /* one prime, exponent > 1 */
    1829         280 :   absolute = flag==2 && lPN==2 && !equali1(gel(EN,1));
    1830             : 
    1831         280 :   Pmod = leafcopy(gel(bid_get_fact(bnr_get_bid(bnr)),1));
    1832         280 :   for (i=1; i<lg(Pmod); i++) gel(Pmod,i) = pr_get_p(gel(Pmod,i));
    1833         280 :   settyp(Pmod, t_VEC);
    1834         280 :   P = ZV_sort_uniq(shallowconcat(nf_get_ramified_primes(nf), Pmod));
    1835             : 
    1836         280 :   res = cgetg(lPN, t_VEC);
    1837         644 :   for (i = 1; i < lPN; i++)
    1838             :   {
    1839         364 :     struct rnfkummer *pkum = NULL;
    1840         364 :     GEN p = gel(PN,i), H = hnfmodid(subgroup, powii(p, gel(EN,i)));
    1841         364 :     long sp = itos(p);
    1842         364 :     absolute &= (lg(H)==2 || equali1(gcoeff(H,2,2))); /* cyclic */
    1843         364 :     if (bnf_get_tuN(bnf) % sp)
    1844             :     {
    1845         133 :       pkum = &kum;
    1846         133 :       rnfkummer_init(pkum, bnf, sp, prec);
    1847             :     }
    1848         364 :     gel(res,i) = bnrclassfield_primepower(pkum, bnr, H, p, P, absolute, prec);
    1849             :   }
    1850         280 :   res = liftpol_shallow(shallowconcat1(res));
    1851         280 :   res = gen_sort(res, (void*)cmp_RgX, gen_cmp_RgX);
    1852         280 :   if (flag)
    1853             :   {
    1854          84 :     res = nfcompositumall(nf, res);
    1855          84 :     if (flag==2 && !absolute) res = rnfequation(nf, res);
    1856             :   }
    1857         280 :   return gerepileupto(av,res);
    1858             : }

Generated by: LCOV version 1.13