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 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 - modules - kummer.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.8.0 lcov report (development 19628-9774e23) Lines: 720 837 86.0 %
Date: 2016-10-01 05:54:29 Functions: 51 56 91.1 %
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         861 : prank(GEN cyc, long ell)
      44             : {
      45             :   long i;
      46        2485 :   for (i=1; i<lg(cyc); i++)
      47        1883 :     if (smodis(gel(cyc,i),ell)) break;
      48         861 :   return i-1;
      49             : }
      50             : 
      51             : /* increment y, which runs through [0,d-1]^(k-1). Return 0 when done. */
      52             : static int
      53         112 : increment(GEN y, long k, long d)
      54             : {
      55         112 :   long i = k, j;
      56             :   do
      57             :   {
      58         154 :     if (--i == 0) return 0;
      59         126 :     y[i]++;
      60         126 :   } while (y[i] >= d);
      61          84 :   for (j = i+1; j < k; j++) y[j] = 0;
      62          84 :   return 1;
      63             : }
      64             : 
      65             : static int
      66         399 : ok_congruence(GEN X, ulong ell, long lW, GEN vecMsup)
      67             : {
      68             :   long i, l;
      69         399 :   if (zv_equal0(X)) return 0;
      70         399 :   l = lg(X);
      71         721 :   for (i=lW; i<l; i++)
      72         343 :     if (X[i] == 0) return 0;
      73         378 :   l = lg(vecMsup);
      74         602 :   for (i=1; i<l; i++)
      75         245 :     if (zv_equal0(Flm_Flc_mul(gel(vecMsup,i),X, ell))) return 0;
      76         357 :   return 1;
      77             : }
      78             : 
      79             : static int
      80         105 : ok_sign(GEN X, GEN msign, GEN arch)
      81             : {
      82         105 :   return zv_equal(Flm_Flc_mul(msign, X, 2), arch);
      83             : }
      84             : 
      85             : /* REDUCTION MOD ell-TH POWERS */
      86             : 
      87             : #if 0
      88             : static GEN
      89             : logarch2arch(GEN x, long r1, long prec)
      90             : {
      91             :   long i, lx;
      92             :   GEN y = cgetg_copy(x, &lx);
      93             :   if (typ(x) == t_MAT)
      94             :   {
      95             :     for (i=1; i<lx; i++) gel(y,i) = logarch2arch(gel(x,i), r1, prec);
      96             :   }
      97             :   else
      98             :   {
      99             :     for (i=1; i<=r1;i++) gel(y,i) = gexp(gel(x,i),prec);
     100             :     for (   ; i<lx; i++) gel(y,i) = gexp(gmul2n(gel(x,i),-1),prec);
     101             :   }
     102             :   return y;
     103             : }
     104             : #endif
     105             : 
     106             : /* make be integral by multiplying by t in (Q^*)^ell */
     107             : static GEN
     108         322 : reduce_mod_Qell(GEN bnfz, GEN be, GEN gell)
     109             : {
     110             :   GEN c;
     111         322 :   be = nf_to_scalar_or_basis(bnfz, be);
     112         322 :   be = Q_primitive_part(be, &c);
     113         322 :   if (c)
     114             :   {
     115         169 :     GEN d, fa = factor(c);
     116         169 :     gel(fa,2) = FpC_red(gel(fa,2), gell);
     117         169 :     d = factorback(fa);
     118         169 :     be = typ(be) == t_INT? mulii(be,d): ZC_Z_mul(be, d);
     119             :   }
     120         322 :   return be;
     121             : }
     122             : 
     123             : /* return q, q^n r = x, v_pr(r) < n for all pr. Insist q is a genuine n-th
     124             :  * root (i.e r = 1) if strict != 0. */
     125             : static GEN
     126        1148 : idealsqrtn(GEN nf, GEN x, GEN gn, int strict)
     127             : {
     128        1148 :   long i, l, n = itos(gn);
     129             :   GEN fa, q, Ex, Pr;
     130             : 
     131        1148 :   fa = idealfactor(nf, x);
     132        1148 :   Pr = gel(fa,1); l = lg(Pr);
     133        1148 :   Ex = gel(fa,2); q = NULL;
     134        3035 :   for (i=1; i<l; i++)
     135             :   {
     136        1887 :     long ex = itos(gel(Ex,i));
     137        1887 :     GEN e = stoi(ex / n);
     138        1887 :     if (strict && ex % n) pari_err_SQRTN("idealsqrtn", fa);
     139        1887 :     if (q) q = idealmulpowprime(nf, q, gel(Pr,i), e);
     140         386 :     else   q = idealpow(nf, gel(Pr,i), e);
     141             :   }
     142        1148 :   return q? q: gen_1;
     143             : }
     144             : 
     145             : static GEN
     146         322 : reducebeta(GEN bnfz, GEN b, GEN ell)
     147             : {
     148         322 :   long prec = nf_get_prec(bnfz);
     149         322 :   GEN y, elllogfu, nf = bnf_get_nf(bnfz), fu = bnf_get_fu_nocheck(bnfz);
     150             : 
     151         322 :   if (DEBUGLEVEL>1) err_printf("reducing beta = %Ps\n",b);
     152         322 :   b = reduce_mod_Qell(nf, b, ell);
     153             :   /* reduce l-th root */
     154         322 :   y = idealsqrtn(nf, b, ell, 0); /* (b) = y^ell I, I integral */
     155         322 :   if (typ(y) == t_MAT && !is_pm1(gcoeff(y,1,1)))
     156             :   {
     157         148 :     GEN T = idealred(nf, mkvec2(y, gen_1)), t = gel(T,2);
     158             :     /* (t)*T[1] = y, T[1] integral and small */
     159         148 :     if (gcmp(idealnorm(nf,t), gen_1) > 0)
     160         134 :       b = nfmul(nf, b, nfpow(nf, t, negi(ell)));
     161             :   }
     162         322 :   if (DEBUGLEVEL>1) err_printf("beta reduced via ell-th root = %Ps\n",b);
     163             :   /* log. embeddings of fu^ell */
     164         322 :   elllogfu = RgM_Rg_mul(real_i(bnf_get_logfu(bnfz)), ell);
     165             :   for (;;)
     166             :   {
     167         343 :     GEN emb, z = get_arch_real(nf, b, &emb, prec);
     168         343 :     if (z)
     169             :     {
     170         322 :       GEN ex = RgM_Babai(elllogfu, z);
     171         322 :       if (ex)
     172             :       {
     173         322 :         b = nfdiv(nf, b, nffactorback(nf, fu, RgC_Rg_mul(ex,ell)));
     174         322 :         break;
     175             :       }
     176             :     }
     177          21 :     prec = precdbl(prec);
     178          21 :     if (DEBUGLEVEL) pari_warn(warnprec,"reducebeta",prec);
     179          21 :     nf = nfnewprec_shallow(nf,prec);
     180          21 :   }
     181         322 :   if (DEBUGLEVEL>1) err_printf("beta LLL-reduced mod U^l = %Ps\n",b);
     182         322 :   return b;
     183             : }
     184             : 
     185             : static GEN
     186        1666 : tauofalg(GEN x, tau_s *tau) {
     187        1666 :   long tx = typ(x);
     188        1666 :   if (tx == t_POLMOD) { x = gel(x,2); tx = typ(x); }
     189        1666 :   if (tx == t_POL) x = RgX_RgXQ_eval(x, tau->x, tau->R);
     190        1666 :   return mkpolmod(x, tau->R);
     191             : }
     192             : 
     193             : static tau_s *
     194         231 : get_tau(tau_s *tau, GEN nf, compo_s *C, long g)
     195             : {
     196         231 :   GEN bas = nf_get_zk(nf), U, Uzk;
     197         231 :   long i, l = lg(bas);
     198             : 
     199             :   /* compute action of tau: q^g + kp */
     200         231 :   U = RgX_add(RgXQ_powu(C->q, g, C->R), RgX_muls(C->p, C->k));
     201         231 :   U = RgX_RgXQ_eval(C->rev, U, C->R);
     202             : 
     203         231 :   tau->x  = U;
     204         231 :   tau->R  = C->R;
     205         231 :   Uzk = cgetg(l, t_MAT);
     206        1603 :   for (i=1; i<l; i++)
     207        1372 :     gel(Uzk,i) = algtobasis(nf, tauofalg(gel(bas,i), tau));
     208         231 :   tau->zk = Uzk; return tau;
     209             : }
     210             : 
     211             : static GEN tauoffamat(GEN x, tau_s *tau);
     212             : 
     213             : static GEN
     214        9716 : tauofelt(GEN x, tau_s *tau)
     215             : {
     216        9716 :   switch(typ(x))
     217             :   {
     218        8344 :     case t_COL: return RgM_RgC_mul(tau->zk, x);
     219        1078 :     case t_MAT: return tauoffamat(x, tau);
     220         294 :     default: return tauofalg(x, tau);
     221             :   }
     222             : }
     223             : static GEN
     224        1246 : tauofvec(GEN x, tau_s *tau)
     225             : {
     226             :   long i, l;
     227        1246 :   GEN y = cgetg_copy(x, &l);
     228        1246 :   for (i=1; i<l; i++) gel(y,i) = tauofelt(gel(x,i), tau);
     229        1246 :   return y;
     230             : }
     231             : /* [x, tau(x), ..., tau^m(x)] */
     232             : static GEN
     233         630 : powtau(GEN x, long m, tau_s *tau)
     234             : {
     235         630 :   GEN y = cgetg(m+1, t_VEC);
     236             :   long i;
     237         630 :   gel(y,1) = x;
     238         630 :   for (i=2; i<=m; i++) gel(y,i) = tauofelt(gel(y,i-1), tau);
     239         630 :   return y;
     240             : }
     241             : /* x^lambda */
     242             : static GEN
     243         539 : lambdaofelt(GEN x, toK_s *T)
     244             : {
     245         539 :   tau_s *tau = T->tau;
     246         539 :   long i, m = T->m;
     247         539 :   GEN y = cgetg(1, t_MAT), powg = T->powg; /* powg[i] = g^i */
     248        1372 :   for (i=1; i<m; i++)
     249             :   {
     250         833 :     y = famat_mul(y, famat_pow(x, gel(powg,m-i)));
     251         833 :     x = tauofelt(x, tau);
     252             :   }
     253         539 :   return famat_mul(y, x);
     254             : }
     255             : static GEN
     256         448 : lambdaofvec(GEN x, toK_s *T)
     257             : {
     258             :   long i, l;
     259         448 :   GEN y = cgetg_copy(x, &l);
     260         448 :   for (i=1; i<l; i++) gel(y,i) = lambdaofelt(gel(x,i), T);
     261         448 :   return y;
     262             : }
     263             : 
     264             : static GEN
     265        1078 : tauoffamat(GEN x, tau_s *tau)
     266             : {
     267        1078 :   return mkmat2(tauofvec(gel(x,1), tau), gel(x,2));
     268             : }
     269             : 
     270             : static GEN
     271         140 : tauofideal(GEN id, tau_s *tau)
     272             : {
     273         140 :   return ZM_hnfmodid(RgM_mul(tau->zk, id), gcoeff(id, 1,1));
     274             : }
     275             : 
     276             : static int
     277         644 : isprimeidealconj(GEN nfz, GEN P, GEN Q, tau_s *tau)
     278             : {
     279         644 :   GEN p = pr_get_p(P);
     280         644 :   GEN x = pr_get_gen(P);
     281         644 :   if (!equalii(p, pr_get_p(Q))
     282         490 :    || pr_get_e(P) != pr_get_e(Q)
     283         490 :    || pr_get_f(P) != pr_get_f(Q)) return 0;
     284         483 :   if (ZV_equal(x, pr_get_gen(Q))) return 1;
     285             :   for(;;)
     286             :   {
     287        1295 :     if (ZC_prdvd(nfz,x,Q)) return 1;
     288         980 :     x = FpC_red(tauofelt(x, tau), p);
     289         980 :     if (ZC_prdvd(nfz,x,P)) return 0;
     290         812 :   }
     291             : }
     292             : 
     293             : static int
     294         980 : isconjinprimelist(GEN nfz, GEN S, GEN pr, tau_s *tau)
     295             : {
     296             :   long i, l;
     297             : 
     298         980 :   if (!tau) return 0;
     299         756 :   l = lg(S);
     300        1085 :   for (i=1; i<l; i++)
     301         644 :     if (isprimeidealconj(nfz, gel(S,i),pr,tau)) return 1;
     302         441 :   return 0;
     303             : }
     304             : 
     305             : /* assume x in basistoalg form */
     306             : static GEN
     307        1036 : downtoK(toK_s *T, GEN x)
     308             : {
     309        1036 :   long degKz = lg(T->invexpoteta1) - 1;
     310        1036 :   GEN y = gmul(T->invexpoteta1, Rg_to_RgC(lift_shallow(x), degKz));
     311        1036 :   return gmodulo(gtopolyrev(y,varn(T->polnf)), T->polnf);
     312             : }
     313             : 
     314             : static GEN
     315           0 : no_sol(long all, long i)
     316             : {
     317           0 :   if (!all) pari_err_BUG(stack_sprintf("kummer [bug%ld]", i));
     318           0 :   return cgetg(1,t_VEC);
     319             : }
     320             : 
     321             : static GEN
     322         308 : get_gell(GEN bnr, GEN subgp, long all)
     323             : {
     324             :   GEN gell;
     325         308 :   if (all && all != -1) return utoipos(labs(all));
     326         287 :   if (!subgp) return ZV_prod(bnr_get_cyc(bnr));
     327         287 :   gell = det(subgp);
     328         287 :   if (typ(gell) != t_INT) pari_err_TYPE("rnfkummer",gell);
     329         287 :   return gell;
     330             : }
     331             : 
     332             : typedef struct {
     333             :   GEN Sm, Sml1, Sml2, Sl, ESml2;
     334             : } primlist;
     335             : 
     336             : static int
     337         301 : build_list_Hecke(primlist *L, GEN nfz, GEN fa, GEN gothf, GEN gell, tau_s *tau)
     338             : {
     339             :   GEN listpr, listex, pr, factell;
     340         301 :   long vp, i, l, ell = itos(gell), degKz = nf_get_degree(nfz);
     341             : 
     342         301 :   if (!fa) fa = idealfactor(nfz, gothf);
     343         301 :   listpr = gel(fa,1);
     344         301 :   listex = gel(fa,2); l = lg(listpr);
     345         301 :   L->Sm  = vectrunc_init(l);
     346         301 :   L->Sml1= vectrunc_init(l);
     347         301 :   L->Sml2= vectrunc_init(l);
     348         301 :   L->Sl  = vectrunc_init(l+degKz);
     349         301 :   L->ESml2=vecsmalltrunc_init(l);
     350        1050 :   for (i=1; i<l; i++)
     351             :   {
     352         749 :     pr = gel(listpr,i);
     353         749 :     vp = itos(gel(listex,i));
     354         749 :     if (!equalii(pr_get_p(pr), gell))
     355             :     {
     356         574 :       if (vp != 1) return 1;
     357         574 :       if (!isconjinprimelist(nfz, L->Sm,pr,tau)) vectrunc_append(L->Sm,pr);
     358             :     }
     359             :     else
     360             :     {
     361         175 :       long e = pr_get_e(pr), vd = (vp-1)*(ell-1)-ell*e;
     362         175 :       if (vd > 0) return 4;
     363         175 :       if (vd==0)
     364             :       {
     365           7 :         if (!isconjinprimelist(nfz, L->Sml1,pr,tau)) vectrunc_append(L->Sml1, pr);
     366             :       }
     367             :       else
     368             :       {
     369         168 :         if (vp==1) return 2;
     370         168 :         if (!isconjinprimelist(nfz, L->Sml2,pr,tau))
     371             :         {
     372         168 :           vectrunc_append(L->Sml2, pr);
     373         168 :           vecsmalltrunc_append(L->ESml2, vp);
     374             :         }
     375             :       }
     376             :     }
     377             :   }
     378         301 :   factell = idealprimedec(nfz,gell); l = lg(factell);
     379         707 :   for (i=1; i<l; i++)
     380             :   {
     381         406 :     pr = gel(factell,i);
     382         406 :     if (!idealval(nfz,gothf,pr))
     383         231 :       if (!isconjinprimelist(nfz, L->Sl,pr,tau)) vectrunc_append(L->Sl, pr);
     384             :   }
     385         301 :   return 0; /* OK */
     386             : }
     387             : 
     388             : /* Return a Flm */
     389             : static GEN
     390         560 : logall(GEN nf, GEN vec, long lW, long mginv, long ell, GEN pr, long ex)
     391             : {
     392         560 :   GEN m, M, bid = Idealstarprk(nf, pr, ex, nf_INIT);
     393         560 :   long ellrank, i, l = lg(vec);
     394             : 
     395         560 :   ellrank = prank(bid_get_cyc(bid), ell);
     396         560 :   M = cgetg(l,t_MAT);
     397        2289 :   for (i=1; i<l; i++)
     398             :   {
     399        1729 :     m = ideallog(nf, gel(vec,i), bid);
     400        1729 :     setlg(m, ellrank+1);
     401        1729 :     if (i < lW) m = gmulsg(mginv, m);
     402        1729 :     gel(M,i) = ZV_to_Flv(m, ell);
     403             :   }
     404         560 :   return M;
     405             : }
     406             : 
     407             : /* compute the u_j (see remark 5.2.15.) */
     408             : static GEN
     409         301 : get_u(GEN cyc, long rc, GEN gell)
     410             : {
     411         301 :   long i, l = lg(cyc);
     412         301 :   GEN u = cgetg(l,t_VEC);
     413         301 :   for (i=1; i<=rc; i++) gel(u,i) = gen_0;
     414         301 :   for (   ; i<  l; i++) gel(u,i) = Fp_inv(gel(cyc,i), gell);
     415         301 :   return u;
     416             : }
     417             : 
     418             : /* alg. 5.2.15. with remark */
     419             : static GEN
     420         413 : isprincipalell(GEN bnfz, GEN id, GEN cycgen, GEN u, GEN gell, long rc)
     421             : {
     422         413 :   long i, l = lg(cycgen);
     423         413 :   GEN logdisc, b, y = bnfisprincipal0(bnfz, id, nf_FORCE|nf_GENMAT);
     424             : 
     425         413 :   logdisc = FpC_red(gel(y,1), gell);
     426         413 :   b = gel(y,2);
     427         609 :   for (i=rc+1; i<l; i++)
     428             :   {
     429         196 :     GEN e = modii(mulii(gel(logdisc,i),gel(u,i)), gell);
     430         196 :     if (signe(e)) b = famat_mul(b, famat_pow(gel(cycgen,i), e));
     431             :   }
     432         413 :   setlg(logdisc,rc+1); return mkvec2(logdisc, b);
     433             : }
     434             : 
     435             : static GEN
     436        1218 : famat_factorback(GEN v, GEN e)
     437             : {
     438        1218 :   long i, l = lg(e);
     439        1218 :   GEN V = cgetg(1, t_MAT);
     440        4305 :   for (i=1; i<l; i++)
     441        3087 :     if (signe(gel(e,i))) V = famat_mul(V, famat_pow(gel(v,i), gel(e,i)));
     442        1218 :   return V;
     443             : }
     444             : 
     445             : static GEN
     446         322 : compute_beta(GEN X, GEN vecWB, GEN ell, GEN bnfz)
     447             : {
     448             :   GEN BE, be;
     449         322 :   BE = famat_reduce(famat_factorback(vecWB, zv_to_ZV(X)));
     450         322 :   gel(BE,2) = centermod(gel(BE,2), ell);
     451         322 :   be = nffactorback(bnfz, BE, NULL);
     452         322 :   be = reducebeta(bnfz, be, ell);
     453         322 :   if (DEBUGLEVEL>1) err_printf("beta reduced = %Ps\n",be);
     454         322 :   return be;
     455             : }
     456             : 
     457             : static GEN
     458         301 : get_Selmer(GEN bnf, GEN cycgen, long rc)
     459             : {
     460         301 :   GEN U = bnf_build_units(bnf), tu = gel(U,1), fu = vecslice(U, 2, lg(U)-1);
     461         301 :   return shallowconcat(shallowconcat(fu,mkvec(tu)), vecslice(cycgen,1,rc));
     462             : }
     463             : 
     464             : GEN
     465       19747 : lift_if_rational(GEN x)
     466             : {
     467             :   long lx, i;
     468             :   GEN y;
     469             : 
     470       19747 :   switch(typ(x))
     471             :   {
     472        2541 :     default: break;
     473             : 
     474             :     case t_POLMOD:
     475       11697 :       y = gel(x,2);
     476       11697 :       if (typ(y) == t_POL)
     477             :       {
     478        2890 :         long d = degpol(y);
     479        2890 :         if (d > 0) return x;
     480         644 :         return (d < 0)? gen_0: gel(y,2);
     481             :       }
     482        8807 :       return y;
     483             : 
     484        2233 :     case t_POL: lx = lg(x);
     485        2233 :       for (i=2; i<lx; i++) gel(x,i) = lift_if_rational(gel(x,i));
     486        2233 :       break;
     487        3276 :     case t_VEC: case t_COL: case t_MAT: lx = lg(x);
     488        3276 :       for (i=1; i<lx; i++) gel(x,i) = lift_if_rational(gel(x,i));
     489             :   }
     490        8050 :   return x;
     491             : }
     492             : 
     493             : /* A column vector representing a subgroup of prime index */
     494             : static GEN
     495           0 : grptocol(GEN H)
     496             : {
     497           0 :   long i, j, l = lg(H);
     498           0 :   GEN col = cgetg(l, t_VECSMALL);
     499           0 :   for (i = 1; i < l; i++)
     500             :   {
     501           0 :     ulong ell = itou( gcoeff(H,i,i) );
     502           0 :     if (ell == 1) col[i] = 0; else { col[i] = ell-1; break; }
     503             :   }
     504           0 :   for (j=i; ++j < l; ) col[j] = itou( gcoeff(H,i,j) );
     505           0 :   return col;
     506             : }
     507             : 
     508             : /* Reorganize kernel basis so that the tests of ok_congruence can be ok
     509             :  * for y[ncyc]=1 and y[1..ncyc]=1 */
     510             : static GEN
     511           0 : fix_kernel(GEN K, GEN M, GEN vecMsup, long lW, long ell)
     512             : {
     513           0 :   pari_sp av = avma;
     514           0 :   long i, j, idx, ffree, dK = lg(K)-1;
     515           0 :   GEN Ki, Kidx = cgetg(dK+1, t_VECSMALL);
     516             : 
     517             :   /* First step: Gauss elimination on vectors lW...lg(M) */
     518           0 :   for (idx = lg(K), i=lg(M); --i >= lW; )
     519             :   {
     520           0 :     for (j=dK; j > 0; j--) if (coeff(K, i, j)) break;
     521           0 :     if (!j)
     522             :     { /* Do our best to ensure that K[dK,i] != 0 */
     523           0 :       if (coeff(K, i, dK)) continue;
     524           0 :       for (j = idx; j < dK; j++)
     525           0 :         if (coeff(K, i, j) && coeff(K, Kidx[j], dK) != ell - 1)
     526           0 :           Flv_add_inplace(gel(K,dK), gel(K,j), ell);
     527             :     }
     528           0 :     if (j != --idx) swap(gel(K, j), gel(K, idx));
     529           0 :     Kidx[idx] = i;
     530           0 :     if (coeff(K,i,idx) != 1)
     531           0 :       Flv_Fl_div_inplace(gel(K,idx), coeff(K,i,idx), ell);
     532           0 :     Ki = gel(K,idx);
     533           0 :     if (coeff(K,i,dK) != 1)
     534             :     {
     535           0 :       ulong t = Fl_sub(coeff(K,i,dK), 1, ell);
     536           0 :       Flv_sub_inplace(gel(K,dK), Flv_Fl_mul(Ki, t, ell), ell);
     537             :     }
     538           0 :     for (j = dK; --j > 0; )
     539             :     {
     540           0 :       if (j == idx) continue;
     541           0 :       if (coeff(K,i,j))
     542           0 :         Flv_sub_inplace(gel(K,j), Flv_Fl_mul(Ki, coeff(K,i,j), ell), ell);
     543             :     }
     544             :   }
     545             :   /* ffree = first vector that is not "free" for the scalar products */
     546           0 :   ffree = idx;
     547             :   /* Second step: for each hyperplane equation in vecMsup, do the same
     548             :    * thing as before. */
     549           0 :   for (i=1; i < lg(vecMsup); i++)
     550             :   {
     551           0 :     GEN Msup = gel(vecMsup,i);
     552             :     ulong dotprod;
     553           0 :     if (lgcols(Msup) != 2) continue;
     554           0 :     Msup = zm_row(Msup, 1);
     555           0 :     for (j=ffree; --j > 0; )
     556             :     {
     557           0 :       dotprod = Flv_dotproduct(Msup, gel(K,j), ell);
     558           0 :       if (dotprod)
     559             :       {
     560           0 :         if (j != --ffree) swap(gel(K, j), gel(K, ffree));
     561           0 :         if (dotprod != 1) Flv_Fl_div_inplace(gel(K, ffree), dotprod, ell);
     562           0 :         break;
     563             :       }
     564             :     }
     565           0 :     if (!j)
     566             :     { /* Do our best to ensure that vecMsup.K[dK] != 0 */
     567           0 :       if (Flv_dotproduct(Msup, gel(K,dK), ell) == 0)
     568             :       {
     569           0 :         for (j = ffree-1; j <= dK; j++)
     570           0 :           if (Flv_dotproduct(Msup, gel(K,j), ell)
     571           0 :               && coeff(K,Kidx[j],dK) != ell-1)
     572           0 :             Flv_add_inplace(gel(K,dK), gel(K,j), ell);
     573             :       }
     574           0 :       continue;
     575             :     }
     576           0 :     Ki = gel(K,ffree);
     577           0 :     dotprod = Flv_dotproduct(Msup, gel(K,dK), ell);
     578           0 :     if (dotprod != 1)
     579             :     {
     580           0 :       ulong t = Fl_sub(dotprod,1,ell);
     581           0 :       Flv_sub_inplace(gel(K,dK), Flv_Fl_mul(Ki,t,ell), ell);
     582             :     }
     583           0 :     for (j = dK; --j > 0; )
     584             :     {
     585           0 :       if (j == ffree) continue;
     586           0 :       dotprod = Flv_dotproduct(Msup, gel(K,j), ell);
     587           0 :       if (dotprod) Flv_sub_inplace(gel(K,j), Flv_Fl_mul(Ki,dotprod,ell), ell);
     588             :     }
     589             :   }
     590           0 :   if (ell == 2)
     591             :   {
     592           0 :     for (i = ffree, j = ffree-1; i <= dK && j; i++, j--)
     593           0 :     { swap(gel(K,i), gel(K,j)); }
     594             :   }
     595             :   /* Try to ensure that y = vec_ei(n, i) gives a good candidate */
     596           0 :   for (i = 1; i < dK; i++) Flv_add_inplace(gel(K,i), gel(K,dK), ell);
     597           0 :   return gerepilecopy(av, K);
     598             : }
     599             : 
     600             : static GEN
     601           0 : Flm_init(long m, long n)
     602             : {
     603           0 :   GEN M = cgetg(n+1, t_MAT);
     604           0 :   long i; for (i = 1; i <= n; i++) gel(M,i) = cgetg(m+1, t_VECSMALL);
     605           0 :   return M;
     606             : }
     607             : static void
     608           0 : Flv_fill(GEN v, GEN y)
     609             : {
     610           0 :   long i, l = lg(y);
     611           0 :   for (i = 1; i < l; i++) v[i] = y[i];
     612           0 : }
     613             : 
     614             : static GEN
     615         504 : get_badbnf(GEN bnf)
     616             : {
     617             :   long i, l;
     618         504 :   GEN bad = gen_1, gen = bnf_get_gen(bnf);
     619         504 :   l = lg(gen);
     620         945 :   for (i = 1; i < l; i++)
     621             :   {
     622         441 :     GEN g = gel(gen,i);
     623         441 :     bad = lcmii(bad, gcoeff(g,1,1));
     624             :   }
     625         504 :   return bad;
     626             : }
     627             : /* Let K base field, L/K described by bnr (conductor f) + H. Return a list of
     628             :  * primes coprime to f*ell of degree 1 in K whose images in Cl_f(K) generate H:
     629             :  * thus they all split in Lz/Kz; t in Kz is such that
     630             :  * t^(1/p) generates Lz => t is an ell-th power in k(pr) for all such primes.
     631             :  * Restrict to primes not dividing
     632             :  * - the index fz of the polynomial defining Kz, or
     633             :  * - the modulus, or
     634             :  * - ell, or
     635             :  * - a generator in bnf.gen or bnfz.gen */
     636             : static GEN
     637         287 : get_prlist(GEN bnr, GEN H, ulong ell, GEN bnfz)
     638             : {
     639         287 :   pari_sp av0 = avma;
     640             :   forprime_t T;
     641             :   ulong p;
     642             :   GEN L, nf, cyc, bad, cond, condZ, Hsofar;
     643         287 :   L = cgetg(1, t_VEC);
     644         287 :   cyc = bnr_get_cyc(bnr);
     645         287 :   nf = bnr_get_nf(bnr);
     646             : 
     647         287 :   cond = gel(bnr_get_mod(bnr), 1);
     648         287 :   condZ = gcoeff(cond,1,1);
     649         287 :   bad = get_badbnf(bnr_get_bnf(bnr));
     650         287 :   if (bnfz)
     651             :   {
     652         217 :     GEN badz = lcmii(get_badbnf(bnfz), nf_get_index(bnf_get_nf(bnfz)));
     653         217 :     bad = mulii(bad,badz);
     654             :   }
     655         287 :   bad = lcmii(muliu(condZ, ell), bad);
     656             :   /* restrict to primes not dividing bad */
     657             : 
     658         287 :   u_forprime_init(&T, 2, ULONG_MAX);
     659         287 :   Hsofar = cgetg(1, t_MAT);
     660        6480 :   while ((p = u_forprime_next(&T)))
     661             :   {
     662             :     GEN LP;
     663             :     long i, l;
     664        6193 :     if (p == ell || !umodiu(bad, p)) continue;
     665        5240 :     LP = idealprimedec_limit_f(nf, utoipos(p), 1);
     666        5240 :     l = lg(LP);
     667        8779 :     for (i = 1; i < l; i++)
     668             :     {
     669        3826 :       pari_sp av = avma;
     670        3826 :       GEN M, P = gel(LP,i), v = bnrisprincipal(bnr, P, 0);
     671        3826 :       if (!hnf_invimage(H, v)) { avma = av; continue; }
     672         855 :       M = shallowconcat(Hsofar, v);
     673         855 :       M = ZM_hnfmodid(M, cyc);
     674         855 :       if (ZM_equal(M, Hsofar)) continue;
     675         511 :       L = shallowconcat(L, mkvec(P));
     676         511 :       Hsofar = M;
     677             :       /* the primes in L generate H */
     678         511 :       if (ZM_equal(M, H)) return gerepilecopy(av0, L);
     679             :     }
     680             :   }
     681           0 :   pari_err_BUG("rnfkummer [get_prlist]");
     682           0 :   return NULL;
     683             : }
     684             : /*Lprz list of prime ideals in Kz that must split completely in Lz/Kz, vecWA
     685             :  * generators for the S-units used to build the Kummer generators. Return
     686             :  * matsmall M such that \prod WA[j]^x[j] ell-th power mod pr[i] iff
     687             :  * \sum M[i,j] x[j] = 0 (mod ell) */
     688             : static GEN
     689         287 : subgroup_info(GEN bnfz, GEN Lprz, long ell, GEN vecWA)
     690             : {
     691         287 :   GEN nfz = bnf_get_nf(bnfz), M, gell = utoipos(ell), Lell = mkvec(gell);
     692         287 :   long i, j, l = lg(vecWA), lz = lg(Lprz);
     693         287 :   M = cgetg(l, t_MAT);
     694         287 :   for (j=1; j<l; j++) gel(M,j) = cgetg(lz, t_VECSMALL);
     695         798 :   for (i=1; i < lz; i++)
     696             :   {
     697         511 :     GEN pr = gel(Lprz,i), EX = subiu(pr_norm(pr), 1);
     698         511 :     GEN N, g,T,p, prM = idealhnf(nfz, pr);
     699         511 :     GEN modpr = zk_to_Fq_init(nfz, &pr,&T,&p);
     700         511 :     long v = Z_lvalrem(divis(EX,ell), ell, &N) + 1; /* Norm(pr)-1 = N * ell^v */
     701         511 :     GEN ellv = powuu(ell, v);
     702         511 :     g = gener_Fq_local(T,p, Lell);
     703         511 :     g = Fq_pow(g,N, T,p); /* order ell^v */
     704        2387 :     for (j=1; j < l; j++)
     705             :     {
     706        1876 :       GEN logc, c = gel(vecWA,j);
     707        1876 :       if (typ(c) == t_MAT) /* famat */
     708        1078 :         c = famat_makecoprime(nfz, gel(c,1), gel(c,2), pr, prM, EX);
     709        1876 :       c = nf_to_Fq(nfz, c, modpr);
     710        1876 :       c = Fq_pow(c, N, T,p);
     711        1876 :       logc = Fq_log(c, g, ellv, T,p);
     712        1876 :       ucoeff(M, i,j) = umodiu(logc, ell);
     713             :     }
     714             :   }
     715         287 :   return M;
     716             : }
     717             : 
     718             : /* if all!=0, give all equations of degree 'all'. Assume bnr modulus is the
     719             :  * conductor */
     720             : static GEN
     721          70 : rnfkummersimple(GEN bnr, GEN subgroup, GEN gell, long all)
     722             : {
     723             :   long ell, i, j, degK, dK;
     724             :   long lSml2, lSl2, lSp, rc, lW;
     725             :   long prec;
     726          70 :   long rk=0, ncyc=0;
     727          70 :   GEN mat=NULL, matgrp=NULL, xell, be1 = NULL;
     728          70 :   long firstpass = all<0;
     729             : 
     730             :   GEN bnf,nf,bid,ideal,arch,cycgen;
     731             :   GEN cyc;
     732             :   GEN Sp,listprSp,matP;
     733          70 :   GEN res=NULL,u,M,K,y,vecMsup,vecW,vecWB,vecBp,msign;
     734             :   primlist L;
     735             : 
     736          70 :   bnf = bnr_get_bnf(bnr); (void)bnf_build_units(bnf);
     737          70 :   nf  = bnf_get_nf(bnf);
     738          70 :   degK = nf_get_degree(nf);
     739             : 
     740          70 :   bid = bnr_get_bid(bnr);
     741          70 :   ideal= bid_get_ideal(bid);
     742          70 :   arch = bid_get_arch(bid); /* this is the conductor */
     743          70 :   ell = itos(gell);
     744          70 :   i = build_list_Hecke(&L, nf, gel(bid,3), ideal, gell, NULL);
     745          70 :   if (i) return no_sol(all,i);
     746             : 
     747          70 :   lSml2 = lg(L.Sml2)-1;
     748          70 :   Sp = shallowconcat(L.Sm, L.Sml1); lSp = lg(Sp)-1;
     749          70 :   listprSp = shallowconcat(L.Sml2, L.Sl); lSl2 = lg(listprSp)-1;
     750             : 
     751          70 :   cycgen = bnf_build_cycgen(bnf);
     752          70 :   cyc = bnf_get_cyc(bnf); rc = prank(cyc, ell);
     753             : 
     754          70 :   vecW = get_Selmer(bnf, cycgen, rc);
     755          70 :   u = get_u(cyc, rc, gell);
     756             : 
     757          70 :   vecBp = cgetg(lSp+1, t_VEC);
     758          70 :   matP  = cgetg(lSp+1, t_MAT);
     759         189 :   for (j=1; j<=lSp; j++)
     760             :   {
     761         119 :     GEN L = isprincipalell(bnf,gel(Sp,j), cycgen,u,gell,rc);
     762         119 :     gel( matP,j) = gel(L,1);
     763         119 :     gel(vecBp,j) = gel(L,2);
     764             :   }
     765          70 :   vecWB = shallowconcat(vecW, vecBp);
     766             : 
     767          70 :   prec = DEFAULTPREC +
     768          70 :       nbits2extraprec(((degK-1) * (gexpo(vecWB) + gexpo(nf_get_M(nf)))));
     769          70 :   if (nf_get_prec(nf) < prec) nf = nfnewprec_shallow(nf, prec);
     770          70 :   msign = nfsign(nf, vecWB);
     771          70 :   arch = ZV_to_zv(arch);
     772             : 
     773          70 :   vecMsup = cgetg(lSml2+1,t_VEC);
     774          70 :   M = NULL;
     775         175 :   for (i=1; i<=lSl2; i++)
     776             :   {
     777         105 :     GEN pr = gel(listprSp,i);
     778         105 :     long e = pr_get_e(pr), z = ell * (e / (ell-1));
     779             : 
     780         105 :     if (i <= lSml2)
     781             :     {
     782          35 :       z += 1 - L.ESml2[i];
     783          35 :       gel(vecMsup,i) = logall(nf, vecWB, 0,0, ell, pr,z+1);
     784             :     }
     785         105 :     M = vconcat(M, logall(nf, vecWB, 0,0, ell, pr,z));
     786             :   }
     787          70 :   lW = lg(vecW);
     788          70 :   M = vconcat(M, shallowconcat(zero_Flm(rc,lW-1), ZM_to_Flm(matP, ell)));
     789          70 :   if (!all)
     790             :   { /* primes landing in subgroup must be totally split */
     791          70 :     GEN Lpr = get_prlist(bnr, subgroup, ell, NULL);
     792          70 :     GEN M2 = subgroup_info(bnf, Lpr, ell, vecWB);
     793          70 :     M = vconcat(M, M2);
     794             :   }
     795          70 :   K = Flm_ker(M, ell);
     796          70 :   dK = lg(K)-1;
     797             : 
     798          70 :   if (all < 0)
     799           0 :     K = fix_kernel(K, M, vecMsup, lW, ell);
     800             : 
     801          70 :   y = cgetg(dK+1,t_VECSMALL);
     802          70 :   if (all) res = cgetg(1,t_VEC); /* in case all = 1 */
     803          70 :   if (all < 0)
     804             :   {
     805           0 :     ncyc = dK;
     806           0 :     mat = Flm_init(dK, ncyc);
     807           0 :     if (all == -1) matgrp = Flm_init(lg(bnr_get_cyc(bnr)), ncyc+1);
     808           0 :     rk = 0;
     809             :   }
     810          70 :   xell = pol_xn(ell, 0);
     811             :   do {
     812          70 :     dK = lg(K)-1;
     813         147 :     while (dK)
     814             :     {
     815          77 :       for (i=1; i<dK; i++) y[i] = 0;
     816          77 :       y[i] = 1; /* y = [0,...,0,1,0,...,0], 1 at dK'th position */
     817             :       do
     818             :       {
     819         147 :         pari_sp av = avma;
     820         147 :         GEN be, P=NULL, X;
     821         147 :         if (all < 0)
     822             :         {
     823           0 :           Flv_fill(gel(mat, rk+1), y);
     824           0 :           setlg(mat, rk+2);
     825           0 :           if (Flm_rank(mat, ell) <= rk) continue;
     826             :         }
     827         147 : FOUND:  X = Flm_Flc_mul(K, y, ell);
     828         147 :         if (ok_congruence(X, ell, lW, vecMsup) && ok_sign(X, msign, arch))
     829             :         {/* be satisfies all congruences, x^ell - be is irreducible, signature
     830             :           * and relative discriminant are correct */
     831          70 :           if (all < 0) rk++;
     832          70 :           be = compute_beta(X, vecWB, gell, bnf);
     833          70 :           be = nf_to_scalar_or_alg(nf, be);
     834          70 :           if (typ(be) == t_POL) be = mkpolmod(be, nf_get_pol(nf));
     835          70 :           if (all == -1)
     836             :           {
     837           0 :             pari_sp av2 = avma;
     838           0 :             GEN Kgrp, colgrp = grptocol(rnfnormgroup(bnr, gsub(xell, be)));
     839           0 :             if (ell != 2)
     840             :             {
     841           0 :               if (rk == 1) be1 = be;
     842             :               else
     843             :               { /* Compute the pesky scalar */
     844           0 :                 GEN K2, C = cgetg(4, t_MAT);
     845           0 :                 gel(C,1) = gel(matgrp,1);
     846           0 :                 gel(C,2) = colgrp;
     847           0 :                 gel(C,3) = grptocol(rnfnormgroup(bnr, gsub(xell, gmul(be1,be))));
     848           0 :                 K2 = Flm_ker(C, ell);
     849           0 :                 if (lg(K2) != 2) pari_err_BUG("linear algebra");
     850           0 :                 K2 = gel(K2,1);
     851           0 :                 if (K2[1] != K2[2])
     852           0 :                   Flv_Fl_mul_inplace(colgrp, Fl_div(K2[2],K2[1],ell), ell);
     853             :               }
     854             :             }
     855           0 :             Flv_fill(gel(matgrp,rk), colgrp);
     856           0 :             setlg(matgrp, rk+1);
     857           0 :             Kgrp = Flm_ker(matgrp, ell);
     858           0 :             if (lg(Kgrp) == 2)
     859             :             {
     860           0 :               setlg(gel(Kgrp,1), rk+1);
     861           0 :               y = Flm_Flc_mul(mat, gel(Kgrp,1), ell);
     862           0 :               all = 0; goto FOUND;
     863             :             }
     864           0 :             avma = av2;
     865             :           }
     866             :           else
     867             :           {
     868          70 :             P = gsub(xell, be);
     869          70 :             if (all)
     870           0 :               res = shallowconcat(res, gerepileupto(av, P));
     871             :             else
     872             :             {
     873          70 :               if (ZM_equal(rnfnormgroup(bnr,P),subgroup)) return P; /*DONE*/
     874           0 :               avma = av; continue;
     875             :             }
     876             :           }
     877           0 :           if (all < 0 && rk == ncyc) return res;
     878           0 :           if (firstpass) break;
     879             :         }
     880          77 :         else avma = av;
     881          77 :       } while (increment(y, dK, ell));
     882           7 :       y[dK--] = 0;
     883             :     }
     884           0 :   } while (firstpass--);
     885           0 :   return all? res: gen_0;
     886             : }
     887             : 
     888             : /* alg. 5.3.11 (return only discrete log mod ell) */
     889             : static GEN
     890         826 : isvirtualunit(GEN bnf, GEN v, GEN cycgen, GEN cyc, GEN gell, long rc)
     891             : {
     892         826 :   GEN L, b, eps, y, q, nf = bnf_get_nf(bnf);
     893         826 :   GEN w = idealsqrtn(nf, v, gell, 1);
     894         826 :   long i, l = lg(cycgen);
     895             : 
     896         826 :   L = bnfisprincipal0(bnf, w, nf_GENMAT|nf_FORCE);
     897         826 :   q = gel(L,1);
     898         826 :   if (gequal0(q)) { eps = v; y = q; }
     899             :   else
     900             :   {
     901         140 :     b = gel(L,2);
     902         140 :     y = cgetg(l,t_COL);
     903         420 :     for (i=1; i<l; i++)
     904         280 :       gel(y,i) = diviiexact(mulii(gell,gel(q,i)), gel(cyc,i));
     905         140 :     eps = famat_mul(famat_factorback(cycgen, y), famat_pow(b, gell));
     906         140 :     eps = famat_mul(famat_inv(eps), v);
     907             :   }
     908         826 :   setlg(y, rc+1);
     909         826 :   b = bnfisunit(bnf,eps);
     910         826 :   if (lg(b) == 1) pari_err_BUG("isvirtualunit");
     911         826 :   return shallowconcat(lift_shallow(b), y);
     912             : }
     913             : 
     914             : /* J a vector of elements in nfz = relative extension of nf by polrel,
     915             :  * return the Steinitz element attached to the module generated by J */
     916             : static GEN
     917         651 : Stelt(GEN nf, GEN J, GEN polrel)
     918             : {
     919         651 :   long i, l = lg(J), vx = varn(polrel);
     920         651 :   GEN A = cgetg(l, t_VEC), I = cgetg(l, t_VEC);
     921        4487 :   for (i = 1; i < l; i++)
     922             :   {
     923        3836 :     GEN v = gel(J,i);
     924        3836 :     if (typ(v) == t_POL) { v = RgX_rem(v, polrel); setvarn(v,vx); }
     925        3836 :     gel(A,i) = v;
     926        3836 :     gel(I,i) = gen_1;
     927             :   }
     928         651 :   A = RgV_to_RgM(A, degpol(polrel));
     929         651 :   return idealprod(nf, gel(nfhnf(nf, mkvec2(A,I)),2));
     930             : }
     931             : 
     932             : static GEN
     933         126 : polrelKzK(toK_s *T, GEN x)
     934             : {
     935         126 :   GEN P = roots_to_pol(powtau(x, T->m, T->tau), 0);
     936         126 :   long i, l = lg(P);
     937         126 :   for (i=2; i<l; i++) gel(P,i) = downtoK(T, gel(P,i));
     938         126 :   return P;
     939             : }
     940             : 
     941             : /* N: Cl_m(Kz) --> Cl_m(K), lift subgroup from bnr to bnrz using Algo 4.1.11 */
     942             : static GEN
     943         126 : invimsubgroup(GEN bnrz, GEN bnr, GEN subgroup, toK_s *T)
     944             : {
     945             :   long l, j;
     946             :   GEN P, cyc, gen, U, polrel, StZk;
     947         126 :   GEN nf = bnr_get_nf(bnr), nfz = bnr_get_nf(bnrz);
     948         126 :   GEN polz = nf_get_pol(nfz), zkz = nf_get_zk(nfz);
     949             : 
     950         126 :   polrel = polrelKzK(T, pol_x(varn(polz)));
     951         126 :   StZk = Stelt(nf, zkz, polrel);
     952         126 :   cyc = bnr_get_cyc(bnrz); l = lg(cyc);
     953         126 :   gen = bnr_get_gen(bnrz);
     954         126 :   P = cgetg(l,t_MAT);
     955         651 :   for (j=1; j<l; j++)
     956             :   {
     957         525 :     GEN g, id = idealhnf_shallow(nfz, gel(gen,j));
     958         525 :     g = Stelt(nf, RgV_RgM_mul(zkz, id), polrel);
     959         525 :     g = idealdiv(nf, g, StZk); /* N_{Kz/K}(gen[j]) */
     960         525 :     gel(P,j) = isprincipalray(bnr, g);
     961             :   }
     962         126 :   (void)ZM_hnfall_i(shallowconcat(P, subgroup), &U, 1);
     963         126 :   setlg(U, l); for (j=1; j<l; j++) setlg(U[j], l);
     964         126 :   return ZM_hnfmodid(U, cyc);
     965             : }
     966             : 
     967             : static GEN
     968         252 : pol_from_Newton(GEN S)
     969             : {
     970         252 :   long i, k, l = lg(S);
     971         252 :   GEN C = cgetg(l+1, t_VEC), c = C + 1;
     972         252 :   gel(c,0) = gen_1;
     973         252 :   gel(c,1) = gel(S,1); /* gen_0 in our case */
     974         882 :   for (k = 2; k < l; k++)
     975             :   {
     976         630 :     GEN s = gel(S,k);
     977         630 :     for (i = 2; i < k-1; i++) s = gadd(s, gmul(gel(S,i), gel(c,k-i)));
     978         630 :     gel(c,k) = gdivgs(s, -k);
     979             :   }
     980         252 :   return gtopoly(C, 0);
     981             : }
     982             : 
     983             : /* - mu_b = sum_{0 <= i < m} floor(r_b r_{d-1-i} / ell) tau^i */
     984             : static GEN
     985         602 : get_mmu(long b, GEN r, long ell)
     986             : {
     987         602 :   long i, m = lg(r)-1;
     988         602 :   GEN M = cgetg(m+1, t_VEC);
     989         602 :   for (i = 0; i < m; i++) gel(M,i+1) = stoi((r[b + 1] * r[m - i]) / ell);
     990         602 :   return M;
     991             : }
     992             : 
     993             : /* coeffs(x, a..b) in variable v >= varn(x) */
     994             : static GEN
     995        5964 : split_pol(GEN x, long v, long a, long b)
     996             : {
     997        5964 :   long i, l = degpol(x);
     998        5964 :   GEN y = x + a, z;
     999             : 
    1000        5964 :   if (l < b) b = l;
    1001        5964 :   if (a > b || varn(x) != v) return pol_0(v);
    1002        5320 :   l = b-a + 3;
    1003        5320 :   z = cgetg(l, t_POL); z[1] = x[1];
    1004        5320 :   for (i = 2; i < l; i++) gel(z,i) = gel(y,i);
    1005        5320 :   return normalizepol_lg(z, l);
    1006             : }
    1007             : 
    1008             : /* return (den_a * z) mod (v^ell - num_a/den_a), assuming deg(z) < 2*ell
    1009             :  * allow either num/den to be NULL (= 1) */
    1010             : static GEN
    1011        2982 : mod_Xell_a(GEN z, long v, long ell, GEN num_a, GEN den_a)
    1012             : {
    1013        2982 :   GEN z1 = split_pol(z, v, ell, degpol(z));
    1014        2982 :   GEN z0 = split_pol(z, v, 0,   ell-1); /* z = v^ell z1 + z0*/
    1015        2982 :   if (den_a) z0 = gmul(den_a, z0);
    1016        2982 :   if (num_a) z1 = gmul(num_a, z1);
    1017        2982 :   return gadd(z0, z1);
    1018             : }
    1019             : static GEN
    1020         854 : to_alg(GEN nfz, GEN c, long v)
    1021             : {
    1022             :   GEN z;
    1023         854 :   if (typ(c) != t_COL) return c;
    1024         602 :   z = gmul(nf_get_zk(nfz), c);
    1025         602 :   if (typ(z) == t_POL) setvarn(z, v);
    1026         602 :   return z;
    1027             : }
    1028             : 
    1029             : /* th. 5.3.5. and prop. 5.3.9. */
    1030             : static GEN
    1031         252 : compute_polrel(GEN nfz, toK_s *T, GEN be, long g, long ell)
    1032             : {
    1033         252 :   long i, k, m = T->m, vT = fetch_var(), vz = fetch_var();
    1034             :   GEN r, powtaubet, S, p1, root, num_t, den_t, nfzpol, powtau_prim_invbe;
    1035             :   GEN prim_Rk, C_Rk, prim_root, C_root, prim_invbe, C_invbe;
    1036             :   pari_timer ti;
    1037             : 
    1038         252 :   r = cgetg(m+1,t_VECSMALL); /* r[i+1] = g^i mod ell */
    1039         252 :   r[1] = 1;
    1040         252 :   for (i=2; i<=m; i++) r[i] = (r[i-1] * g) % ell;
    1041         252 :   powtaubet = powtau(be, m, T->tau);
    1042         252 :   if (DEBUGLEVEL>1) { err_printf("Computing Newton sums: "); timer_start(&ti); }
    1043         252 :   prim_invbe = Q_primitive_part(nfinv(nfz, be), &C_invbe);
    1044         252 :   powtau_prim_invbe = powtau(prim_invbe, m, T->tau);
    1045             : 
    1046         252 :   root = cgetg(ell + 2, t_POL);
    1047         252 :   root[1] = evalsigne(1) | evalvarn(0);
    1048         252 :   for (i = 0; i < ell; i++) gel(root,2+i) = gen_0;
    1049         854 :   for (i = 0; i < m; i++)
    1050             :   { /* compute (1/be) ^ (-mu) instead of be^mu [mu << 0].
    1051             :      * 1/be = C_invbe * prim_invbe */
    1052         602 :     GEN mmu = get_mmu(i, r, ell);
    1053             :     /* p1 = prim_invbe ^ -mu */
    1054         602 :     p1 = to_alg(nfz, nffactorback(nfz, powtau_prim_invbe, mmu), vz);
    1055         602 :     if (C_invbe) p1 = gmul(p1, powgi(C_invbe, RgV_sumpart(mmu, m)));
    1056             :     /* root += zeta_ell^{r_i} T^{r_i} be^mu_i */
    1057         602 :     gel(root, 2 + r[i+1]) = monomial(p1, r[i+1], vT);
    1058             :   }
    1059             :   /* Other roots are as above with z_ell --> z_ell^j.
    1060             :    * Treat all contents (C_*) and principal parts (prim_*) separately */
    1061         252 :   prim_Rk = prim_root = Q_primitive_part(root, &C_root);
    1062         252 :   C_Rk = C_root;
    1063             : 
    1064         252 :   r = vecsmall_reverse(r); /* theta^ell = be^( sum tau^a r_{d-1-a} ) */
    1065             :   /* Compute modulo X^ell - 1, T^ell - t, nfzpol(vz) */
    1066         252 :   p1 = to_alg(nfz, nffactorback(nfz, powtaubet, r), vz);
    1067         252 :   num_t = Q_remove_denom(p1, &den_t);
    1068             : 
    1069         252 :   nfzpol = leafcopy(nf_get_pol(nfz));
    1070         252 :   setvarn(nfzpol, vz);
    1071         252 :   S = cgetg(ell+1, t_VEC); /* Newton sums */
    1072         252 :   gel(S,1) = gen_0;
    1073         882 :   for (k = 2; k <= ell; k++)
    1074             :   { /* compute the k-th Newton sum */
    1075         630 :     pari_sp av = avma;
    1076         630 :     GEN z, D, Rk = gmul(prim_Rk, prim_root);
    1077         630 :     C_Rk = mul_content(C_Rk, C_root);
    1078         630 :     Rk = mod_Xell_a(Rk, 0, ell, NULL, NULL); /* mod X^ell - 1 */
    1079        3010 :     for (i = 2; i < lg(Rk); i++)
    1080             :     {
    1081        2380 :       if (typ(gel(Rk,i)) != t_POL) continue;
    1082        2352 :       z = mod_Xell_a(gel(Rk,i), vT, ell, num_t,den_t); /* mod T^ell - t */
    1083        2352 :       gel(Rk,i) = RgXQX_red(z, nfzpol); /* mod nfz.pol */
    1084             :     }
    1085         630 :     if (den_t) C_Rk = mul_content(C_Rk, ginv(den_t));
    1086         630 :     prim_Rk = Q_primitive_part(Rk, &D);
    1087         630 :     C_Rk = mul_content(C_Rk, D); /* root^k = prim_Rk * C_Rk */
    1088             : 
    1089             :     /* Newton sum is ell * constant coeff (in X), which has degree 0 in T */
    1090         630 :     z = polcoeff_i(prim_Rk, 0, 0);
    1091         630 :     z = polcoeff_i(z      , 0,vT);
    1092         630 :     z = downtoK(T, gmulgs(z, ell));
    1093         630 :     if (C_Rk) z = gmul(z, C_Rk);
    1094         630 :     gerepileall(av, C_Rk? 3: 2, &z, &prim_Rk, &C_Rk);
    1095         630 :     if (DEBUGLEVEL>1) { err_printf("%ld(%ld) ", k, timer_delay(&ti)); err_flush(); }
    1096         630 :     gel(S,k) = z;
    1097             :   }
    1098         252 :   if (DEBUGLEVEL>1) err_printf("\n");
    1099         252 :   (void)delete_var();
    1100         252 :   (void)delete_var(); return pol_from_Newton(S);
    1101             : }
    1102             : 
    1103             : /* lift elt t in nf to nfz, algebraic form */
    1104             : static GEN
    1105         350 : lifttoKz(GEN nf, GEN t, compo_s *C)
    1106             : {
    1107         350 :   GEN x = nf_to_scalar_or_alg(nf, t);
    1108         350 :   if (typ(x) != t_POL) return x;
    1109         350 :   return RgX_RgXQ_eval(x, C->p, C->R);
    1110             : }
    1111             : /* lift ideal id in nf to nfz */
    1112             : static GEN
    1113         231 : ideallifttoKz(GEN nfz, GEN nf, GEN id, compo_s *C)
    1114             : {
    1115         231 :   GEN I = idealtwoelt(nf,id);
    1116         231 :   GEN x = nf_to_scalar_or_alg(nf, gel(I,2));
    1117         231 :   if (typ(x) != t_POL) return gel(I,1);
    1118         147 :   gel(I,2) = algtobasis(nfz, RgX_RgXQ_eval(x, C->p, C->R));
    1119         147 :   return idealhnf_two(nfz,I);
    1120             : }
    1121             : /* lift ideal pr in nf to ONE prime in nfz (the others are conjugate under tau
    1122             :  * and bring no further information on e_1 W). Assume pr coprime to
    1123             :  * index of both nf and nfz, and unramified in Kz/K (minor simplification) */
    1124             : static GEN
    1125         378 : prlifttoKz(GEN nfz, GEN nf, GEN pr, compo_s *C)
    1126             : {
    1127         378 :   GEN F, p = pr_get_p(pr), t = pr_get_gen(pr), T = nf_get_pol(nfz);
    1128         378 :   if (nf_get_degree(nf) != 1)
    1129             :   { /* restrict to primes above pr */
    1130         350 :     t = Q_primpart( lifttoKz(nf,t,C) );
    1131         350 :     T = FpX_gcd(FpX_red(T,p), FpX_red(t,p), p);
    1132         350 :     T = FpX_normalize(T, p);
    1133             :   }
    1134         378 :   F = FpX_factor(T, p);
    1135         378 :   return idealprimedec_kummer(nfz,gcoeff(F,1,1), pr_get_e(pr), p);
    1136             : }
    1137             : static GEN
    1138         217 : get_przlist(GEN L, GEN nfz, GEN nf, compo_s *C)
    1139             : {
    1140             :   long i, l;
    1141         217 :   GEN M = cgetg_copy(L, &l);
    1142         217 :   for (i = 1; i < l; i++) gel(M,i) = prlifttoKz(nfz, nf, gel(L,i), C);
    1143         217 :   return M;
    1144             : }
    1145             : 
    1146             : static void
    1147         231 : compositum_red(compo_s *C, GEN P, GEN Q)
    1148             : {
    1149         231 :   GEN p, q, a, z = gel(compositum2(P, Q),1);
    1150         231 :   a = gel(z,1);
    1151         231 :   p = gel(gel(z,2), 2);
    1152         231 :   q = gel(gel(z,3), 2);
    1153         231 :   C->k = itos( gel(z,4) );
    1154             :   /* reduce R. FIXME: should be polredbest(a, 1), but breaks rnfkummer bench */
    1155         231 :   z = polredabs0(a, nf_ORIG|nf_PARTIALFACT);
    1156         231 :   C->R = gel(z,1);
    1157         231 :   a = gel(gel(z,2), 2);
    1158         231 :   C->p = RgX_RgXQ_eval(p, a, C->R);
    1159         231 :   C->q = RgX_RgXQ_eval(q, a, C->R);
    1160         231 :   C->rev = QXQ_reverse(a, C->R);
    1161         231 :   if (DEBUGLEVEL>1) err_printf("polred(compositum) = %Ps\n",C->R);
    1162         231 : }
    1163             : 
    1164             : /* replace P->C^(-deg P) P(xC) for the largest integer C such that coefficients
    1165             :  * remain algebraic integers. Lift *rational* coefficients */
    1166             : static void
    1167         252 : nfX_Z_normalize(GEN nf, GEN P)
    1168             : {
    1169             :   long i, l;
    1170         252 :   GEN C, Cj, PZ = cgetg_copy(P, &l);
    1171         252 :   PZ[1] = P[1];
    1172        1386 :   for (i = 2; i < l; i++) /* minor variation on RgX_to_nfX (create PZ) */
    1173             :   {
    1174        1134 :     GEN z = nf_to_scalar_or_basis(nf, gel(P,i));
    1175        1134 :     if (typ(z) == t_INT)
    1176         775 :       gel(PZ,i) = gel(P,i) = z;
    1177             :     else
    1178         359 :       gel(PZ,i) = ZV_content(z);
    1179             :   }
    1180         252 :   (void)ZX_Z_normalize(PZ, &C);
    1181             : 
    1182         504 :   if (C == gen_1) return;
    1183          57 :   Cj = C;
    1184         242 :   for (i = l-2; i > 1; i--)
    1185             :   {
    1186         185 :     if (i != l-2) Cj = mulii(Cj, C);
    1187         185 :     gel(P,i) = gdiv(gel(P,i), Cj);
    1188             :   }
    1189             : }
    1190             : 
    1191             : /* v[i] = g^i mod ell, 1 <= i < m */
    1192             : static GEN
    1193         231 : Fl_powers_FpV(ulong g, long m, ulong ell)
    1194             : {
    1195         231 :   GEN v = cgetg(m, t_VEC);
    1196         231 :   ulong gi = g % ell;
    1197             :   long i;
    1198         231 :   for (i=1; i<m; i++) { gel(v,i) = utoi(gi); gi = (gi*g) % ell; }
    1199         231 :   return v;
    1200             : }
    1201             : 
    1202             : static GEN
    1203         308 : _rnfkummer(GEN bnr, GEN subgroup, long all, long prec)
    1204             : {
    1205             :   long ell, i, j, m, d, dK, dc, rc, ru, rv, g, mginv, degK, degKz, vnf;
    1206             :   long lSp, lSml2, lSl2, lW;
    1207             :   GEN polnf,bnf,nf,bnfz,nfz,bid,ideal,cycgen,gell,p1,p2,vselmer;
    1208             :   GEN cyc,gen;
    1209             :   GEN Q,idealz,gothf;
    1210         308 :   GEN res=NULL,u,M,K,y,vecMsup,vecW,vecWA,vecWB,vecB,vecC,vecAp,vecBp;
    1211             :   GEN matP, Sp, listprSp, Tc, Tv, P;
    1212             :   primlist L;
    1213             :   toK_s T;
    1214             :   tau_s _tau, *tau;
    1215             :   compo_s COMPO;
    1216             :   pari_timer t;
    1217         308 :   long rk=0, ncyc=0;
    1218         308 :   GEN mat = NULL;
    1219         308 :   long firstpass = all<0;
    1220             : 
    1221         308 :   if (DEBUGLEVEL) timer_start(&t);
    1222         308 :   checkbnr(bnr);
    1223         308 :   bnf = bnr_get_bnf(bnr);
    1224         308 :   nf  = bnf_get_nf(bnf);
    1225         308 :   polnf = nf_get_pol(nf); vnf = varn(polnf);
    1226         308 :   if (!vnf) pari_err_PRIORITY("rnfkummer", polnf, "=", 0);
    1227             :   /* step 7 */
    1228         308 :   p1 = bnrconductor_i(bnr, subgroup, 2);
    1229         308 :   if (DEBUGLEVEL) timer_printf(&t, "[rnfkummer] conductor");
    1230         308 :   bnr      = gel(p1,2);
    1231         308 :   subgroup = gel(p1,3);
    1232         308 :   gell = get_gell(bnr,subgroup,all);
    1233         308 :   ell = itos(gell);
    1234         308 :   if (ell == 1) return pol_x(0);
    1235         308 :   if (!uisprime(ell)) pari_err_IMPL("kummer for composite relative degree");
    1236         308 :   if (all && all != -1 && umodiu(bnr_get_no(bnr), ell))
    1237           7 :     return cgetg(1, t_VEC);
    1238         301 :   if (bnf_get_tuN(bnf) % ell == 0)
    1239          70 :     return rnfkummersimple(bnr, subgroup, gell, all);
    1240             : 
    1241         231 :   if (all == -1) all = 0;
    1242         231 :   bid = bnr_get_bid(bnr);
    1243         231 :   ideal = bid_get_ideal(bid);
    1244             :   /* step 1 of alg 5.3.5. */
    1245         231 :   if (DEBUGLEVEL>2) err_printf("Step 1\n");
    1246         231 :   compositum_red(&COMPO, polnf, polcyclo(ell,vnf));
    1247             :   /* step 2 */
    1248         231 :   if (DEBUGLEVEL>2) err_printf("Step 2\n");
    1249         231 :   if (DEBUGLEVEL) timer_printf(&t, "[rnfkummer] compositum");
    1250         231 :   degK  = degpol(polnf);
    1251         231 :   degKz = degpol(COMPO.R);
    1252         231 :   m = degKz / degK;
    1253         231 :   d = (ell-1) / m;
    1254         231 :   g = (long)Fl_powu(pgener_Fl(ell), d, ell);
    1255         231 :   if (Fl_powu((ulong)g, m, ell*ell) == 1) g += ell;
    1256             :   /* ord(g) = m in all (Z/ell^k)^* */
    1257             :   /* step 3 */
    1258         231 :   if (DEBUGLEVEL>2) err_printf("Step 3\n");
    1259             :   /* could factor disc(R) using th. 2.1.6. */
    1260         231 :   bnfz = Buchall(COMPO.R, nf_FORCE, maxss(prec,BIGDEFAULTPREC));
    1261         231 :   if (DEBUGLEVEL) timer_printf(&t, "[rnfkummer] bnfinit(Kz)");
    1262         231 :   cycgen = bnf_build_cycgen(bnfz);
    1263         231 :   nfz = bnf_get_nf(bnfz);
    1264         231 :   cyc = bnf_get_cyc(bnfz); rc = prank(cyc,ell);
    1265         231 :   gen = bnf_get_gen(bnfz);
    1266         231 :   u = get_u(cyc, rc, gell);
    1267             : 
    1268         231 :   vselmer = get_Selmer(bnfz, cycgen, rc);
    1269         231 :   if (DEBUGLEVEL) timer_printf(&t, "[rnfkummer] Selmer group");
    1270         231 :   ru = (degKz>>1)-1;
    1271         231 :   rv = rc+ru+1;
    1272         231 :   tau = get_tau(&_tau, nfz, &COMPO, g);
    1273             : 
    1274             :   /* step 4 */
    1275         231 :   if (DEBUGLEVEL>2) err_printf("Step 4\n");
    1276         231 :   vecB=cgetg(rc+1,t_VEC);
    1277         231 :   Tc=cgetg(rc+1,t_MAT);
    1278         371 :   for (j=1; j<=rc; j++)
    1279             :   {
    1280         140 :     p1 = tauofideal(gel(gen,j), tau);
    1281         140 :     p1 = isprincipalell(bnfz, p1, cycgen,u,gell,rc);
    1282         140 :     gel(Tc,j)  = gel(p1,1);
    1283         140 :     gel(vecB,j)= gel(p1,2);
    1284             :   }
    1285             : 
    1286         231 :   vecC = cgetg(rc+1,t_VEC);
    1287         231 :   if (rc)
    1288             :   {
    1289         126 :     for (j=1; j<=rc; j++) gel(vecC,j) = cgetg(1, t_MAT);
    1290         126 :     p1 = cgetg(m,t_VEC);
    1291         126 :     gel(p1,1) = matid(rc);
    1292         126 :     for (j=2; j<=m-1; j++) gel(p1,j) = gmul(gel(p1,j-1),Tc);
    1293         126 :     p2 = vecB;
    1294         294 :     for (j=1; j<=m-1; j++)
    1295             :     {
    1296         168 :       GEN z = FpM_red(gmulsg((j*d)%ell,gel(p1,m-j)), gell);
    1297         168 :       p2 = tauofvec(p2, tau);
    1298         364 :       for (i=1; i<=rc; i++)
    1299         196 :         gel(vecC,i) = famat_mul(gel(vecC,i), famat_factorback(p2,gel(z,i)));
    1300             :     }
    1301         126 :     for (i=1; i<=rc; i++) gel(vecC,i) = famat_reduce(gel(vecC,i));
    1302             :   }
    1303             :   /* step 5 */
    1304         231 :   if (DEBUGLEVEL>2) err_printf("Step 5\n");
    1305         231 :   Tv = cgetg(rv+1,t_MAT);
    1306        1057 :   for (j=1; j<=rv; j++)
    1307             :   {
    1308         826 :     p1 = tauofelt(gel(vselmer,j), tau);
    1309         826 :     if (typ(p1) == t_MAT) /* famat */
    1310         140 :       p1 = nffactorback(nfz, gel(p1,1), FpC_red(gel(p1,2),gell));
    1311         826 :     gel(Tv,j) = isvirtualunit(bnfz, p1, cycgen,cyc,gell,rc);
    1312             :   }
    1313         231 :   P = FpM_ker(RgM_Rg_add_shallow(Tv, stoi(-g)), gell);
    1314         231 :   lW = lg(P); vecW = cgetg(lW,t_VEC);
    1315         231 :   for (j=1; j<lW; j++) gel(vecW,j) = famat_factorback(vselmer, gel(P,j));
    1316             :   /* step 6 */
    1317         231 :   if (DEBUGLEVEL>2) err_printf("Step 6\n");
    1318         231 :   Q = FpM_ker(RgM_Rg_add_shallow(shallowtrans(Tc), stoi(-g)), gell);
    1319             :   /* step 8 */
    1320         231 :   if (DEBUGLEVEL>2) err_printf("Step 8\n");
    1321         231 :   p1 = RgXQ_matrix_pow(COMPO.p, degKz, degK, COMPO.R);
    1322         231 :   T.invexpoteta1 = RgM_inv(p1); /* left inverse */
    1323         231 :   T.polnf = polnf;
    1324         231 :   T.tau = tau;
    1325         231 :   T.m = m;
    1326         231 :   T.powg = Fl_powers_FpV(g, m, ell);
    1327             : 
    1328         231 :   idealz = ideallifttoKz(nfz, nf, ideal, &COMPO);
    1329         231 :   if (umodiu(gcoeff(ideal,1,1), ell)) gothf = idealz;
    1330             :   else
    1331             :   { /* ell | N(ideal) */
    1332         126 :     GEN bnrz = Buchray(bnfz, idealz, nf_INIT|nf_GEN);
    1333         126 :     GEN subgroupz = invimsubgroup(bnrz, bnr, subgroup, &T);
    1334         126 :     gothf = bnrconductor_i(bnrz,subgroupz,0);
    1335             :   }
    1336             :   /* step 9, 10, 11 */
    1337         231 :   if (DEBUGLEVEL>2) err_printf("Step 9, 10 and 11\n");
    1338         231 :   i = build_list_Hecke(&L, nfz, NULL, gothf, gell, tau);
    1339         231 :   if (i) return no_sol(all,i);
    1340             : 
    1341         231 :   lSml2 = lg(L.Sml2)-1;
    1342         231 :   Sp = shallowconcat(L.Sm, L.Sml1); lSp = lg(Sp)-1;
    1343         231 :   listprSp = shallowconcat(L.Sml2, L.Sl); lSl2 = lg(listprSp)-1;
    1344             : 
    1345             :   /* step 12 */
    1346         231 :   if (DEBUGLEVEL>2) err_printf("Step 12\n");
    1347         231 :   vecAp = cgetg(lSp+1, t_VEC);
    1348         231 :   vecBp = cgetg(lSp+1, t_VEC);
    1349         231 :   matP  = cgetg(lSp+1, t_MAT);
    1350             : 
    1351         385 :   for (j=1; j<=lSp; j++)
    1352             :   {
    1353             :     GEN e, a;
    1354         154 :     p1 = isprincipalell(bnfz, gel(Sp,j), cycgen,u,gell,rc);
    1355         154 :     e = gel(p1,1); gel(matP,j) = e;
    1356         154 :     a = gel(p1,2);
    1357         154 :     gel(vecBp,j) = famat_mul(famat_factorback(vecC, gneg(e)), a);
    1358             :   }
    1359         231 :   vecAp = lambdaofvec(vecBp, &T);
    1360             :   /* step 13 */
    1361         231 :   if (DEBUGLEVEL>2) err_printf("Step 13\n");
    1362         231 :   vecWA = shallowconcat(vecW, vecAp);
    1363         231 :   vecWB = shallowconcat(vecW, vecBp);
    1364             : 
    1365             :   /* step 14, 15, and 17 */
    1366         231 :   if (DEBUGLEVEL>2) err_printf("Step 14, 15 and 17\n");
    1367         231 :   mginv = (m * Fl_inv(g,ell)) % ell;
    1368         231 :   vecMsup = cgetg(lSml2+1,t_VEC);
    1369         231 :   M = NULL;
    1370         518 :   for (i=1; i<=lSl2; i++)
    1371             :   {
    1372         287 :     GEN pr = gel(listprSp,i);
    1373         287 :     long e = pr_get_e(pr), z = ell * (e / (ell-1));
    1374             : 
    1375         287 :     if (i <= lSml2)
    1376             :     {
    1377         133 :       z += 1 - L.ESml2[i];
    1378         133 :       gel(vecMsup,i) = logall(nfz, vecWA,lW,mginv,ell, pr,z+1);
    1379             :     }
    1380         287 :     M = vconcat(M, logall(nfz, vecWA,lW,mginv,ell, pr,z));
    1381             :   }
    1382         231 :   dc = lg(Q)-1;
    1383         231 :   if (dc)
    1384             :   {
    1385         105 :     GEN QtP = gmul(shallowtrans(Q), matP);
    1386         105 :     M = vconcat(M, shallowconcat(zero_Flm(dc,lW-1), ZM_to_Flm(QtP,ell)));
    1387             :   }
    1388         231 :   if (!M) M = zero_Flm(1, lSp + lW - 1);
    1389             : 
    1390         231 :   if (!all)
    1391             :   { /* primes landing in subgroup must be totally split */
    1392         217 :     GEN lambdaWB = shallowconcat(lambdaofvec(vecW, &T), vecAp);/*vecWB^lambda*/
    1393         217 :     GEN Lpr = get_prlist(bnr, subgroup, ell, bnfz);
    1394         217 :     GEN Lprz= get_przlist(Lpr, nfz, nf, &COMPO);
    1395         217 :     GEN M2 = subgroup_info(bnfz, Lprz, ell, lambdaWB);
    1396         217 :     M = vconcat(M, M2);
    1397             :   }
    1398             :   /* step 16 */
    1399         231 :   if (DEBUGLEVEL>2) err_printf("Step 16\n");
    1400         231 :   K = Flm_ker(M, ell);
    1401         231 :   if (all < 0)
    1402           0 :     K = fix_kernel(K, M, vecMsup, lW, ell);
    1403             :   /* step 18 & ff */
    1404         231 :   if (DEBUGLEVEL>2) err_printf("Step 18\n");
    1405         231 :   dK = lg(K)-1;
    1406         231 :   y = cgetg(dK+1,t_VECSMALL);
    1407         231 :   if (all) res = cgetg(1, t_VEC);
    1408         231 :   if (DEBUGLEVEL) timer_printf(&t, "[rnfkummer] candidate list");
    1409         231 :   if (all < 0) { ncyc = dK; rk = 0; mat = zero_Flm(lg(M)-1, ncyc); }
    1410             : 
    1411             :   do {
    1412         231 :     dK = lg(K)-1;
    1413         483 :     while (dK)
    1414             :     {
    1415         238 :       for (i=1; i<dK; i++) y[i] = 0;
    1416         238 :       y[i] = 1; /* y = [0,...,0,1,0,...,0], 1 at dK'th position */
    1417             :       do
    1418             :       { /* cf. algo 5.3.18 */
    1419         252 :         GEN H, be, P, X = Flm_Flc_mul(K, y, ell);
    1420         252 :         if (ok_congruence(X, ell, lW, vecMsup))
    1421             :         {
    1422         252 :           pari_sp av = avma;
    1423         252 :           if (all < 0)
    1424             :           {
    1425           0 :             gel(mat, rk+1) = X;
    1426           0 :             if (Flm_rank(mat,ell) <= rk) continue;
    1427           0 :             rk++;
    1428             :           }
    1429         252 :           be = compute_beta(X, vecWB, gell, bnfz);
    1430         252 :           P = compute_polrel(nfz, &T, be, g, ell);
    1431         252 :           nfX_Z_normalize(nf, P);
    1432         252 :           if (DEBUGLEVEL>1) err_printf("polrel(beta) = %Ps\n", P);
    1433         252 :           if (!all) {
    1434         217 :             H = rnfnormgroup(bnr, P);
    1435         217 :             if (ZM_equal(subgroup, H)) return P; /* DONE */
    1436           0 :             avma = av; continue;
    1437             :           } else {
    1438          35 :             GEN P0 = Q_primpart(lift_shallow(P));
    1439          35 :             GEN g = nfgcd(P0, RgX_deriv(P0), polnf, nf_get_index(nf));
    1440          35 :             if (degpol(g)) continue;
    1441          35 :             H = rnfnormgroup(bnr, P);
    1442          35 :             if (!ZM_equal(subgroup,H) && !bnrisconductor(bnr,H)) continue;
    1443             :           }
    1444          35 :           P = gerepilecopy(av, P);
    1445          35 :           res = shallowconcat(res, P);
    1446          35 :           if (all < 0 && rk == ncyc) return res;
    1447          35 :           if (firstpass) break;
    1448             :         }
    1449          35 :       } while (increment(y, dK, ell));
    1450          21 :       y[dK--] = 0;
    1451             :     }
    1452          14 :   } while (firstpass--);
    1453          14 :   if (all) return res;
    1454           0 :   return gen_0; /* FAIL */
    1455             : }
    1456             : 
    1457             : GEN
    1458         308 : rnfkummer(GEN bnr, GEN subgroup, long all, long prec)
    1459             : {
    1460         308 :   pari_sp av = avma;
    1461         308 :   return gerepilecopy(av, _rnfkummer(bnr, subgroup, all, prec));
    1462             : }

Generated by: LCOV version 1.11