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

Generated by: LCOV version 1.13