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-bordeaux1.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 16375-9f41ae0) Lines: 774 889 87.1 %
Date: 2014-04-19 Functions: 55 60 91.7 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 385 573 67.2 %

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

Generated by: LCOV version 1.9