Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - ellrank.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.14.0 lcov report (development 27775-aca467eab2) Lines: 1107 1185 93.4 %
Date: 2022-07-03 07:33:15 Functions: 105 112 93.8 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2020  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             : #include "pari.h"
      15             : #include "paripriv.h"
      16             : 
      17             : #define DEBUGLEVEL DEBUGLEVEL_ellrank
      18             : 
      19             : static long LIM1 = 10000;
      20             : static long LIMTRIV = 10000;
      21             : 
      22             : /*******************************************************************/
      23             : /*               NFHILBERT and LOCAL SOLUBILITY                    */
      24             : /*       adapted from Denis Simon's original C implementation      */
      25             : /*******************************************************************/
      26             : /* p > 2, T ZX, p prime, x t_INT */
      27             : static long
      28           0 : lemma6(GEN T, GEN p, long nu, GEN x)
      29             : {
      30             :   long la, mu;
      31           0 :   GEN y = ZX_Z_eval(T, x);
      32             : 
      33           0 :   if (Zp_issquare(y, p)) return 1;
      34           0 :   la = Z_pval(y, p); y = ZX_Z_eval(ZX_deriv(T), x);
      35           0 :   if (!signe(y)) return la >= (nu<<1)? 0: -1;
      36           0 :   mu = Z_pval(y,p); if (la > (mu<<1)) return 1;
      37           0 :   return (la >= (nu<<1) && mu >= nu)? 0: -1;
      38             : }
      39             : static long
      40           0 : lemma7_aux(long nu, long la, long r)
      41             : {
      42           0 :   long nu2 = nu << 1;
      43           0 :   return (la >= nu2 || (la == nu2 - 2 && r == 1))? 0: -1;
      44             : }
      45             : /* p = 2, T ZX, x t_INT: return 1 = yes, -1 = no, 0 = inconclusive */
      46             : static long
      47           0 : lemma7(GEN T, long nu, GEN x)
      48             : {
      49             :   long r, la, mu;
      50           0 :   GEN y = ZX_Z_eval(T, x);
      51             : 
      52           0 :   if (!signe(y)) return 1;
      53           0 :   la = Z_lvalrem(y, 2, &y);
      54           0 :   r = Mod8(y); if (!odd(la) && r == 1) return 1;
      55           0 :   r &= 3; /* T(x) / 2^oo mod 4 */
      56           0 :   y = ZX_Z_eval(ZX_deriv(T), x);
      57           0 :   if (!signe(y)) return lemma7_aux(nu, la, r);
      58           0 :   mu = vali(y); if (la > mu<<1) return 1;
      59           0 :   if (nu <= mu) return lemma7_aux(nu, la, r);
      60             :   /* la <= 2mu, mu < nu */
      61           0 :   if (!odd(la) && mu + nu - la <= (r == 1? 2: 1)) return 1;
      62           0 :   return -1;
      63             : }
      64             : 
      65             : /* T a ZX, p a prime, pnu = p^nu, x0 t_INT */
      66             : static long
      67           0 : zpsol(GEN T, GEN p, long nu, GEN pnu, GEN x0)
      68             : {
      69             :   long i, res;
      70           0 :   pari_sp av = avma, btop;
      71             :   GEN x, pnup;
      72             : 
      73           0 :   res = absequaliu(p,2)? lemma7(T,nu,x0): lemma6(T,p,nu,x0);
      74           0 :   set_avma(av);
      75           0 :   if (res== 1) return 1;
      76           0 :   if (res==-1) return 0;
      77           0 :   x = x0; pnup = mulii(pnu,p);
      78           0 :   btop = avma;
      79           0 :   for (i=0; i < itos(p); i++)
      80             :   {
      81           0 :     x = addii(x,pnu);
      82           0 :     if (zpsol(T,p,nu+1,pnup,x)) return gc_long(av,1);
      83           0 :     if (gc_needed(btop, 2))
      84             :     {
      85           0 :       x = gerepileupto(btop, x);
      86           0 :       if (DEBUGMEM > 1)
      87           0 :         pari_warn(warnmem, "hyperell_locally_soluble: %ld/%Ps",i,p);
      88             :     }
      89             :   }
      90           0 :   return gc_long(av,0);
      91             : }
      92             : 
      93             : /* return 1 if equation y^2=T(x) has a rational p-adic solution (possibly
      94             :  * infinite), 0 otherwise. */
      95             : long
      96           0 : hyperell_locally_soluble(GEN T,GEN p)
      97             : {
      98           0 :   pari_sp av = avma;
      99             :   long res;
     100           0 :   if (typ(T)!=t_POL) pari_err_TYPE("hyperell_locally_soluble",T);
     101           0 :   if (typ(p)!=t_INT) pari_err_TYPE("hyperell_locally_soluble",p);
     102           0 :   RgX_check_ZX(T, "hyperell_locally_soluble");
     103           0 :   res = zpsol(T,p,0,gen_1,gen_0) || zpsol(RgX_recip_i(T), p, 1, p, gen_0);
     104           0 :   return gc_long(av, res);
     105             : }
     106             : 
     107             : /* is t a square in (O_K/pr) ? Assume v_pr(t) = 0 */
     108             : static long
     109         140 : quad_char(GEN nf, GEN t, GEN pr)
     110             : {
     111         140 :   GEN T, p, modpr = zk_to_Fq_init(nf, &pr,&T,&p);
     112         140 :   return Fq_issquare(nf_to_Fq(nf,t,modpr), T, p)? 1: -1;
     113             : }
     114             : /* quad_char(x), x in Z, nonzero mod p */
     115             : static long
     116         161 : Z_quad_char(GEN x, GEN pr)
     117             : {
     118         161 :   long f = pr_get_f(pr);
     119         161 :   if (!odd(f)) return 1;
     120         154 :   return kronecker(x, pr_get_p(pr));
     121             : }
     122             : 
     123             : /* (pr,2) = 1. return 1 if x in Z_K is a square in Z_{K_pr}, 0 otherwise.
     124             :  * modpr = zkmodprinit(nf,pr) */
     125             : static long
     126           0 : psquarenf(GEN nf,GEN x,GEN pr,GEN modpr)
     127             : {
     128           0 :   pari_sp av = avma;
     129           0 :   GEN p = pr_get_p(pr);
     130             :   long v;
     131             : 
     132           0 :   x = nf_to_scalar_or_basis(nf, x);
     133           0 :   if (typ(x) == t_INT) {
     134           0 :     if (!signe(x)) return 1;
     135           0 :     v = Z_pvalrem(x, p, &x) * pr_get_e(pr);
     136           0 :     if (v&1) return 0;
     137           0 :     v = (Z_quad_char(x, pr) == 1);
     138             :   } else {
     139           0 :     v = ZC_nfvalrem(x, pr, &x);
     140           0 :     if (v&1) return 0;
     141           0 :     v = (quad_char(nf, x, modpr) == 1);
     142             :   }
     143           0 :   return gc_long(av,v);
     144             : }
     145             : 
     146             : static long
     147        5908 : ZV_iseven(GEN zlog)
     148             : {
     149        5908 :   long i, l = lg(zlog);
     150       21546 :   for (i = 1; i < l; i++)
     151       21413 :     if (mpodd(gel(zlog,i))) return 0;
     152         133 :   return 1;
     153             : }
     154             : 
     155             : /* pr | 2, project to principal units (trivializes later discrete log) */
     156             : static GEN
     157      162435 : to_principal_unit(GEN nf, GEN x, GEN pr, GEN sprk)
     158             : {
     159      162435 :   if (pr_get_f(pr) != 1)
     160             :   {
     161       18172 :     GEN prk = gel(sprk,3);
     162       18172 :     x = nfpowmodideal(nf, x, gmael(sprk,5,1), prk);
     163             :   }
     164      162435 :   return x;
     165             : }
     166             : /* pr | 2. Return 1 if x in Z_K is square in Z_{K_pr}, 0 otherwise */
     167             : static int
     168         511 : psquare2nf(GEN nf, GEN x, GEN pr, GEN sprk)
     169             : {
     170         511 :   long v = nfvalrem(nf, x, pr, &x);
     171         511 :   if (v == LONG_MAX) return 1; /* x = 0 */
     172             :   /* (x,pr) = 1 */
     173         511 :   if (odd(v)) return 0;
     174         490 :   x = to_principal_unit(nf, x, pr, sprk); /* = 1 mod pr */
     175         490 :   return ZV_iseven(sprk_log_prk1(nf, x, sprk));
     176             : }
     177             : 
     178             : /* pr above an odd prime */
     179             : static long
     180           0 : lemma6nf(GEN nf, GEN T, GEN pr, long nu, GEN x, GEN modpr)
     181             : {
     182             :   long la, mu;
     183           0 :   GEN y = nfpoleval(nf, T, x);
     184             : 
     185           0 :   if (psquarenf(nf,y,pr,modpr)) return 1;
     186           0 :   la = nfval(nf, y, pr); y = nfpoleval(nf, RgX_deriv(T), x);
     187           0 :   if (gequal0(y)) return la >= (nu<<1)? 0: -1;
     188           0 :   mu = nfval(nf, y, pr); if (la > (mu<<1)) return 1;
     189           0 :   return (la >= (nu<<1) && mu >= nu)? 0: -1;
     190             : }
     191             : /* pr above 2 */
     192             : static long
     193        5642 : lemma7nf(GEN nf, GEN T, GEN pr, long nu, GEN x, GEN sprk)
     194             : {
     195             :   long i, res, la, mu, q, e, v;
     196        5642 :   GEN M, y, gpx, loggx = NULL, gx = nfpoleval(nf, T, x);
     197             : 
     198        5642 :   la = nfvalrem(nf, gx, pr, &gx); /* gx /= pi^la, pi a pr-uniformizer */
     199        5642 :   if (la == LONG_MAX) return 1;
     200        5635 :   if (!odd(la))
     201             :   {
     202        5418 :     gx = to_principal_unit(nf, gx, pr, sprk); /* now 1 mod pr */
     203        5418 :     loggx = sprk_log_prk1(nf, gx, sprk); /* cheap */
     204        5418 :     if (ZV_iseven(loggx)) return 1;
     205             :   }
     206        5509 :   gpx = nfpoleval(nf, RgX_deriv(T), x);
     207        5509 :   mu = gequal0(gpx)? la+nu+1 /* oo */: nfval(nf,gpx,pr);
     208             : 
     209        5509 :   if (la > (mu << 1)) return 1;
     210        5509 :   if (nu > mu)
     211             :   {
     212          35 :     if (odd(la)) return -1;
     213          35 :     q = mu+nu-la; res = 1;
     214             :   }
     215             :   else
     216             :   {
     217        5474 :     q = (nu << 1) - la;
     218        5474 :     if (q <= 0) return 0;
     219        5173 :     if (odd(la)) return -1;
     220        4977 :     res = 0;
     221             :   }
     222             :   /* la even */
     223        5012 :   e = pr_get_e(pr);
     224        5012 :   if (q > e << 1)  return -1;
     225        4935 :   if (q == 1) return res;
     226             : 
     227             :   /* gx = 1 mod pr; square mod pi^q ? */
     228        4935 :   v = nfvalrem(nf, nfadd(nf, gx, gen_m1), pr, &y);
     229        4935 :   if (v >= q) return res;
     230             :   /* 1 + pi^v y = (1 + pi^vz z)^2 mod pr^q ? v < q <= 2e => vz < e => vz = vy/2
     231             :    * => y = x^2 mod pr^(min(q-v, e+v/2)), (y,pr) = 1 */
     232        4690 :   if (odd(v)) return -1;
     233             :   /* e > 1 */
     234        2037 :   M = cgetg(2*e+1 - q + 1, t_VEC);
     235        4074 :   for (i = q+1; i <= 2*e+1; i++) gel(M, i-q) = sprk_log_gen_pr(nf, sprk, i);
     236        2037 :   M = ZM_hnfmodid(shallowconcat1(M), gen_2);
     237        2037 :   return hnf_solve(M,loggx)? res: -1;
     238             : }
     239             : /* zinit either a sprk (pr | 2) or a modpr structure (pr | p odd).
     240             :    pnu = pi^nu, pi a uniformizer */
     241             : static long
     242        5642 : zpsolnf(GEN nf,GEN T,GEN pr,long nu,GEN pnu,GEN x0,GEN repr,GEN zinit)
     243             : {
     244             :   long i, res;
     245        5642 :   pari_sp av = avma;
     246             :   GEN pnup;
     247             : 
     248        5642 :   res = typ(zinit) == t_VEC? lemma7nf(nf,T,pr,nu,x0,zinit)
     249        5642 :                            : lemma6nf(nf,T,pr,nu,x0,zinit);
     250        5642 :   set_avma(av);
     251        5642 :   if (res== 1) return 1;
     252        5509 :   if (res==-1) return 0;
     253         574 :   pnup = nfmul(nf, pnu, pr_get_gen(pr));
     254         574 :   nu++;
     255        5558 :   for (i=1; i<lg(repr); i++)
     256             :   {
     257        5250 :     GEN x = nfadd(nf, x0, nfmul(nf,pnu,gel(repr,i)));
     258        5250 :     if (zpsolnf(nf,T,pr,nu,pnup,x,repr,zinit)) return gc_long(av,1);
     259             :   }
     260         308 :   return gc_long(av,0);
     261             : }
     262             : 
     263             : /* Let y = copy(x); y[k] := j; return y */
     264             : static GEN
     265        3206 : ZC_add_coeff(GEN x, long k, long j)
     266        3206 : { GEN y = shallowcopy(x); gel(y, k) = utoipos(j); return y; }
     267             : 
     268             : /* system of representatives for Zk/pr */
     269             : static GEN
     270         252 : repres(GEN nf, GEN pr)
     271             : {
     272         252 :   long f = pr_get_f(pr), N = nf_get_degree(nf), p = itos(pr_get_p(pr));
     273         252 :   long i, j, k, pi, pf = upowuu(p, f);
     274         252 :   GEN perm = pr_basis_perm(nf, pr), rep = cgetg(pf+1,t_VEC);
     275             : 
     276         252 :   gel(rep,1) = zerocol(N);
     277        1141 :   for (pi=i=1; i<=f; i++,pi*=p)
     278             :   {
     279         889 :     long t = perm[i];
     280        1778 :     for (j=1; j<p; j++)
     281        4095 :       for (k=1; k<=pi; k++) gel(rep, j*pi+k) = ZC_add_coeff(gel(rep,k), t, j);
     282             :   }
     283         252 :   return rep;
     284             : }
     285             : 
     286             : /* = 1 if equation y^2 = z^deg(T) * T(x/z) has a pr-adic rational solution
     287             :  * (possibly (1,y,0) = oo), 0 otherwise.
     288             :  * coeffs of T are algebraic integers in nf */
     289             : static long
     290         259 : locally_soluble(GEN nf,GEN T,GEN pr)
     291             : {
     292             :   GEN repr, zinit;
     293             : 
     294         259 :   if (typ(T)!=t_POL) pari_err_TYPE("nf_hyperell_locally_soluble",T);
     295         259 :   if (gequal0(T)) return 1;
     296         259 :   checkprid(pr);
     297         259 :   if (absequaliu(pr_get_p(pr), 2))
     298             :   { /* tough case */
     299         259 :     zinit = log_prk_init(nf, pr, 1+2*pr_get_e(pr), NULL);
     300         259 :     if (psquare2nf(nf,constant_coeff(T),pr,zinit)) return 1;
     301         252 :     if (psquare2nf(nf, leading_coeff(T),pr,zinit)) return 1;
     302             :   }
     303             :   else
     304             :   {
     305           0 :     zinit = zkmodprinit(nf, pr);
     306           0 :     if (psquarenf(nf,constant_coeff(T),pr,zinit)) return 1;
     307           0 :     if (psquarenf(nf, leading_coeff(T),pr,zinit)) return 1;
     308             :   }
     309         252 :   repr = repres(nf,pr); /* FIXME: inefficient if Npr is large */
     310         392 :   return zpsolnf(nf, T, pr, 0, gen_1, gen_0, repr, zinit) ||
     311         140 :          zpsolnf(nf, RgX_recip_i(T), pr, 1, pr_get_gen(pr),
     312             :                  gen_0, repr, zinit);
     313             : }
     314             : long
     315         259 : nf_hyperell_locally_soluble(GEN nf,GEN T,GEN pr)
     316             : {
     317         259 :   pari_sp av = avma;
     318         259 :   return gc_long(av, locally_soluble(nf, T, pr));
     319             : }
     320             : 
     321             : /* return a * denom(a)^2, as an 'liftalg' */
     322             : static GEN
     323         518 : den_remove(GEN nf, GEN a)
     324             : {
     325             :   GEN da;
     326         518 :   a = nf_to_scalar_or_basis(nf, a);
     327         518 :   switch(typ(a))
     328             :   {
     329          49 :     case t_INT: return a;
     330           0 :     case t_FRAC: return mulii(gel(a,1), gel(a,2));
     331         469 :     case t_COL:
     332         469 :       a = Q_remove_denom(a, &da);
     333         469 :       if (da) a = ZC_Z_mul(a, da);
     334         469 :       return nf_to_scalar_or_alg(nf, a);
     335           0 :     default: pari_err_TYPE("nfhilbert",a);
     336             :       return NULL;/*LCOV_EXCL_LINE*/
     337             :   }
     338             : }
     339             : 
     340             : static long
     341         259 : hilb2nf(GEN nf,GEN a,GEN b,GEN p)
     342             : {
     343         259 :   pari_sp av = avma;
     344             :   GEN pol;
     345         259 :   a = den_remove(nf, a);
     346         259 :   b = den_remove(nf, b);
     347         259 :   pol = mkpoln(3, a, gen_0, b);
     348             :   /* varn(nf.pol) = 0, pol is not a valid GEN  [as in Pol([x,x], x)].
     349             :    * But it is only used as a placeholder, hence it is not a problem */
     350         259 :   return gc_long(av, nf_hyperell_locally_soluble(nf,pol,p)? 1: -1);
     351             : }
     352             : 
     353             : /* local quadratic Hilbert symbol (a,b)_pr, for a,b (nonzero) in nf */
     354             : static long
     355         567 : nfhilbertp(GEN nf, GEN a, GEN b, GEN pr)
     356             : {
     357             :   GEN t;
     358             :   long va, vb;
     359         567 :   pari_sp av = avma;
     360             : 
     361         567 :   if (absequaliu(pr_get_p(pr), 2)) return hilb2nf(nf,a,b,pr);
     362             : 
     363             :   /* pr not above 2, compute t = tame symbol */
     364         308 :   va = nfval(nf,a,pr);
     365         308 :   vb = nfval(nf,b,pr);
     366         308 :   if (!odd(va) && !odd(vb)) return gc_long(av,1);
     367             :   /* Trick: pretend the exponent is 2, result is OK up to squares ! */
     368         301 :   t = famat_makecoprime(nf, mkvec2(a,b), mkvec2s(vb, -va),
     369             :                         pr, pr_hnf(nf, pr), gen_2);
     370             :   /* quad. symbol is image of t = (-1)^(v(a)v(b)) a^v(b) b^(-v(a))
     371             :    * by the quadratic character  */
     372         301 :   switch(typ(t))
     373             :   {
     374         147 :     default: /* t_COL */
     375         147 :       if (!ZV_isscalar(t)) break;
     376           7 :       t = gel(t,1); /* fall through */
     377         161 :     case t_INT:
     378         161 :       if (odd(va) && odd(vb)) t = negi(t);
     379         161 :       return gc_long(av,  Z_quad_char(t, pr));
     380             :   }
     381         140 :   if (odd(va) && odd(vb)) t = ZC_neg(t);
     382         140 :   return gc_long(av, quad_char(nf, t, pr));
     383             : }
     384             : 
     385             : /* Global quadratic Hilbert symbol (a,b):
     386             :  *  =  1 if X^2 - aY^2 - bZ^2 has a point in projective plane
     387             :  *  = -1 otherwise
     388             :  * a, b should be nonzero */
     389             : long
     390          21 : nfhilbert(GEN nf, GEN a, GEN b)
     391             : {
     392          21 :   pari_sp av = avma;
     393             :   long i, l;
     394             :   GEN S, S2, Sa, Sb, sa, sb;
     395             : 
     396          21 :   a = nf_to_scalar_or_basis(nf, a);
     397          21 :   b = nf_to_scalar_or_basis(nf, b);
     398             :   /* local solutions in real completions ? [ error in nfsign if arg is 0 ]*/
     399          21 :   sa = nfsign(nf, a);
     400          21 :   sb = nfsign(nf, b); l = lg(sa);
     401          35 :   for (i=1; i<l; i++)
     402          21 :     if (sa[i] && sb[i])
     403             :     {
     404           7 :       if (DEBUGLEVEL>3)
     405           0 :         err_printf("nfhilbert not soluble at real place %ld\n",i);
     406           7 :       return gc_long(av,-1);
     407             :     }
     408             : 
     409             :   /* local solutions in finite completions ? (pr | 2ab)
     410             :    * primes above 2 are toughest. Try the others first */
     411          14 :   Sa = idealfactor(nf, a);
     412          14 :   Sb = idealfactor(nf, b);
     413          14 :   S2 = idealfactor(nf, gen_2);
     414          14 :   S = merge_factor(Sa, Sb, (void*)&cmp_prime_ideal, &cmp_nodata);
     415          14 :   S = merge_factor(S,  S2, (void*)&cmp_prime_ideal, &cmp_nodata);
     416          14 :   S = gel(S,1);
     417             :   /* product of all hilbertp is 1 ==> remove one prime (above 2!) */
     418          28 :   for (i=lg(S)-1; i>1; i--)
     419          21 :     if (nfhilbertp(nf,a,b,gel(S,i)) < 0)
     420             :     {
     421           7 :       if (DEBUGLEVEL>3)
     422           0 :         err_printf("nfhilbert not soluble at finite place %Ps\n",S[i]);
     423           7 :       return gc_long(av,-1);
     424             :     }
     425           7 :   return gc_long(av,1);
     426             : }
     427             : 
     428             : long
     429         581 : nfhilbert0(GEN nf,GEN a,GEN b,GEN p)
     430             : {
     431         581 :   nf = checknf(nf);
     432         581 :   if (p) {
     433         560 :     checkprid(p);
     434         560 :     if (gequal0(a)) pari_err_DOMAIN("nfhilbert", "a", "=", gen_0, a);
     435         553 :     if (gequal0(b)) pari_err_DOMAIN("nfhilbert", "b", "=", gen_0, b);
     436         546 :     return nfhilbertp(nf,a,b,p);
     437             :   }
     438          21 :   return nfhilbert(nf,a,b);
     439             : }
     440             : 
     441             : /*******************************************************************/
     442             : /*                      HYPERELL_LOCAL_SOLVE                       */
     443             : /*******************************************************************/
     444             : 
     445             : /* Based on
     446             : T.A. Fisher and G.F. Sills
     447             : Local solubility and height bounds for coverings of elliptic curves
     448             : https://www.dpmms.cam.ac.uk/~taf1000/papers/htbounds.pdf
     449             : */
     450             : 
     451             : static int
     452      114394 : FpX_issquare(GEN q, GEN p)
     453             : {
     454      114394 :   GEN F = FpX_factor_squarefree(q,p);
     455      114394 :   long i, l = lg(F);
     456      153230 :   for (i = 1; i < l; i+=2)
     457      114107 :     if (degpol(gel(F,i)) > 0) return 0;
     458       39123 :   return 1;
     459             : }
     460             : 
     461             : static GEN
     462      113946 : hyperell_red(GEN q, GEN p)
     463             : {
     464             :   GEN Q;
     465      113946 :   long v = ZX_pvalrem(q, p, &Q);
     466      113946 :   if (v == 1) return q;
     467      113946 :   return odd(v)? ZX_Z_mul(Q, p): Q;
     468             : }
     469             : 
     470             : static GEN
     471      115773 : hyperell_reg_point(GEN q, GEN p)
     472             : {
     473             :   GEN Q, F;
     474      115773 :   long i, l, v = ZX_pvalrem(q, p, &Q);
     475      115773 :   if (v != 1) q = odd(v)? ZX_Z_mul(Q, p): Q;
     476      115773 :   if (!odd(v))
     477             :   {
     478      114394 :     GEN qr = FpX_red(q, p);
     479      114394 :     if (!FpX_issquare(qr,p) || Fp_issquare(leading_coeff(qr), p))
     480      113946 :       return mkvec2s(0,1);
     481             :   }
     482        1827 :   F = FpX_roots(Q, p); l = lg(F);
     483        1834 :   for (i = 1; i < l; i++)
     484             :   {
     485        1827 :     GEN r = gel(F,i), s = hyperell_reg_point(ZX_affine(q, p, r), p);
     486        1827 :     if (s)
     487        1820 :       retmkvec2(addii(r,mulii(gel(s,1),p)), mulii(gel(s,2),p));
     488             :   }
     489           7 :   return NULL;
     490             : }
     491             : 
     492             : static GEN
     493        1225 : hyperell_lift(GEN q, GEN x, GEN p)
     494             : {
     495             :   long e;
     496        1225 :   GEN qy2 = ZX_Z_sub(q, sqri(p));
     497        1225 :   for (e = 2; ; e<<=1)
     498         742 :   {
     499        1967 :     pari_sp av = avma;
     500        1967 :     GEN z = ZpX_liftroot(qy2, x, p, e);
     501        1967 :     if (signe(z) == 0) z = powiu(p, e);
     502        1967 :     if (Zp_issquare(poleval(q, z), p)) return z;
     503         742 :     set_avma(av);
     504             :   }
     505             : }
     506             : 
     507             : static GEN
     508       56973 : affine_apply(GEN r, GEN x)
     509             : {
     510       56973 :   return addii(mulii(gel(r,2),x), gel(r,1));
     511             : }
     512             : 
     513             : static GEN
     514       56973 : Qp_hyperell_solve_odd(GEN q, GEN p)
     515             : {
     516       56973 :   GEN qi = RgX_recip_shallow(q);
     517       56973 :   GEN r = hyperell_reg_point(q,  p), qr = NULL, qrp;
     518       56973 :   GEN s = hyperell_reg_point(qi, p), qs = NULL, qsp;
     519       56973 :   if (!r && !s) return NULL;
     520       56973 :   if (r)
     521             :   {
     522       56973 :     qr = hyperell_red(ZX_affine(q,  gel(r,2), gel(r,1)), p);
     523       56973 :     qrp = FpX_deriv(qr, p);
     524             :   }
     525       56973 :   if (s)
     526             :   {
     527       56973 :     qs = hyperell_red(ZX_affine(qi, gel(s,2), gel(s,1)), p);
     528       56973 :     qsp = FpX_deriv(qs, p);
     529             :   }
     530             :   while(1)
     531       14707 :   {
     532       71680 :     GEN x = randomi(p);
     533       71680 :     if (r)
     534             :     {
     535       71680 :       GEN y2 = FpX_eval(qr, x, p);
     536       71680 :       if (Fp_issquare(y2,p))
     537             :       {
     538       47005 :          if (signe(y2))
     539       42938 :            return affine_apply(r,x);
     540        4067 :          if (signe(FpX_eval(qrp, x, p)))
     541             :          {
     542         924 :            x = hyperell_lift(qr, x, p);
     543         924 :            return affine_apply(r,x);
     544             :          }
     545             :       }
     546             :     }
     547       27818 :     if (s)
     548             :     {
     549       27818 :       GEN y2 = FpX_eval(qs, x, p);
     550       27818 :       if (Fp_issquare(y2,p))
     551             :       {
     552       14833 :          if (signe(x)==0) x = p;
     553       14833 :          if (signe(y2))
     554       12810 :            return ginv(affine_apply(s,x));
     555        2023 :          if (signe(FpX_eval(qsp, x, p)))
     556             :          {
     557         301 :            GEN xl = hyperell_lift(qs, x, p);
     558         301 :            return ginv(affine_apply(s,xl));
     559             :          }
     560             :       }
     561             :     }
     562             :   }
     563             : }
     564             : 
     565             : static GEN
     566         399 : Q2_hyperell_lift(GEN p, GEN q, long x, long y)
     567             : {
     568             :   GEN T, h;
     569             :   long e;
     570         399 :   if (y==0) y = 2;
     571         399 :   T = ZX_sub(p, ZX_Z_add(ZX_mulu(q, y), sqru(y)));
     572         399 :   h = ZX_add(ZX_sqr(q), ZX_shifti(p, 2));
     573         399 :   for (e = 3; ; e++)
     574         861 :   {
     575        1260 :     pari_sp av = avma;
     576        1260 :     GEN r = ZpX_liftroot(T, utoi(x), gen_2, e);
     577        1260 :     if (signe(r) == 0) r = int2n(e);
     578        1260 :     if (Zp_issquare(poleval(h, r), gen_2)) return r;
     579         861 :     set_avma(av);
     580             :   }
     581             :   return NULL;/*LCOV_EXCL_LINE*/
     582             : }
     583             : 
     584             : static GEN
     585        2751 : Q2_hyperell_regpoint(GEN P, GEN Q)
     586             : {
     587             :   long x;
     588        2751 :   GEN p = ZX_to_F2x(P), dp = F2x_deriv(p);
     589        2751 :   GEN q = ZX_to_F2x(Q), dq = F2x_deriv(q);
     590             : 
     591        4557 :   for (x = 0; x <= 1; x++)
     592             :   {
     593        4032 :     long px = F2x_eval(p,x), qx = F2x_eval(q,x);
     594             :     long dpx, dqx;
     595        4032 :     if (qx == 1)
     596             :     {
     597        1981 :       if(px == 0) return x==0 ? gen_2: gen_1;
     598         154 :       continue;
     599             :     }
     600        2051 :     dpx = F2x_eval(dp,x);
     601        2051 :     dqx = F2x_eval(dq,x);
     602        2051 :     if (odd(dqx * px + dpx))
     603         399 :       return Q2_hyperell_lift(P, Q, x, px);
     604             :   }
     605         525 :   return NULL;
     606             : }
     607             : 
     608             : static GEN
     609        2751 : Q2_hyperell_solve_affine(GEN p, GEN q)
     610             : {
     611        2751 :   pari_sp av = avma;
     612             :   GEN R, p4, q4;
     613             :   long x;
     614             :   while(1)
     615        2247 :   {
     616             :     GEN pp, p0, p1;
     617        4998 :     long vp = ZX_lval(p, 2);
     618        4998 :     long vq = ZX_lval(q, 2);
     619        4998 :     long w = minss(vp>>1, vq);
     620        4998 :     p = ZX_shifti(p, -2*w);
     621        4998 :     q = ZX_shifti(q, -w);
     622        4998 :     if (ZX_lval(q,2)==0) break;
     623        2520 :     pp = RgX_splitting(p,2); p0 = gel(pp,1); p1 = gel(pp,2);
     624        2520 :     if (ZX_lval(p1,2)==0 || ZX_lval(p0,2)>=1) break;
     625        2247 :     p = ZX_sub(p, ZX_mul(p0, ZX_add(q, p0)));
     626        2247 :     q = ZX_add(q, ZX_shifti(p0, 1));
     627             :   }
     628        2751 :   R = Q2_hyperell_regpoint(p, q);
     629        2751 :   if (R) return gerepileuptoint(av, R);
     630         525 :   p4 = ZX_to_Flx(p,4);
     631         525 :   q4 = ZX_to_Flx(q,4);
     632         672 :   for (x = 0; x <= 1; x++)
     633             :   {
     634         665 :     ulong px = Flx_eval(p4, x, 4);
     635         665 :     ulong qx = Flx_eval(q4, x, 4);
     636         665 :     if (px == 0 || (1+qx+3*px)%4==0)
     637             :     {
     638         518 :       GEN p2 = ZX_affine(p, gen_2, utoi(x));
     639         518 :       GEN q2 = ZX_affine(q, gen_2, utoi(x));
     640         518 :       GEN S = Q2_hyperell_solve_affine(p2, q2);
     641         518 :       if (S) return gerepileuptoint(av, addiu(shifti(S,1),x));
     642             :     }
     643             :   }
     644           7 :   return gc_NULL(av);
     645             : }
     646             : 
     647             : static GEN
     648        2226 : Q2_hyperell_solve(GEN P)
     649             : {
     650        2226 :   long v = varn(P);
     651        2226 :   GEN S = Q2_hyperell_solve_affine(P, pol_0(v));
     652        2226 :   if (!S) S = ginv(Q2_hyperell_solve_affine(RgX_recip(P), pol_0(v)));
     653        2226 :   return S;
     654             : }
     655             : 
     656             : static GEN
     657       59199 : hyperell_local_solve(GEN q, GEN p)
     658             : {
     659       59199 :   if (equaliu(p,2))
     660        2226 :     return Q2_hyperell_solve(q);
     661       56973 :   return Qp_hyperell_solve_odd(q, p);
     662             : }
     663             : 
     664             : /*******************************************************************/
     665             : /*                         BINARY QUARTIC                          */
     666             : /*******************************************************************/
     667             : static int
     668      105101 : Qp_issquare(GEN a, GEN p)
     669             : {
     670      105101 :   GEN b = typ(a) == t_INT? a: mulii(gel(a,1), gel(a,2));
     671      105101 :   return Zp_issquare(b, p);
     672             : }
     673             : 
     674             : static GEN
     675       17913 : quartic_IJ(GEN g)
     676             : {
     677       17913 :   GEN a = gel(g, 6), b = gel(g, 5), c = gel(g, 4), d = gel(g, 3), e = gel(g, 2);
     678       17913 :   GEN ae = gmul(a,e), bd = gmul(b,d), c2 = gsqr(c);
     679             :   /* 12ae - 3bd + c^2 */
     680       17913 :   GEN I = gadd(gsub(gmulsg(12, ae), gmulsg(3, bd)), c2);
     681             :   /* c(72ae + 9bd - 2c^2) - 27ad^2 - 27eb^2*/
     682       17913 :   GEN J = gsub(gmul(c, gsub(gadd(gmulsg(72,ae), gmulsg(9,bd)), gmul2n(c2,1))),
     683             :                gmulsg(27, gadd(gmul(a, gsqr(d)), gmul(gsqr(b), e))));
     684       17913 :   return mkvec2(I, J);
     685             : }
     686             : 
     687             : static GEN
     688        2226 : quartic_hessiandd(GEN g)
     689             : {
     690        2226 :   GEN a = gel(g, 6), b = gel(g, 5), c = gel(g, 4), d = gel(g, 3), e = gel(g, 2);
     691        2226 :   GEN a8 = gmul2n(a, 3), p4 = gsub(gmulsg(3, gsqr(b)), gmul(a8, c));
     692        2226 :   GEN p3 = gsub(gmul(b, c), gmul(gmulsg(6, a), d));
     693        2226 :   GEN p2 = gsub(gmulsg(8, gsqr(c)), gmulsg(12, gadd(gmul(b, d), gmul(a8, e))));
     694        2226 :   return deg2pol_shallow(gmulgu(p4,12), gmulgu(p3,24), p2, 1);
     695             : }
     696             : 
     697             : static GEN
     698       12131 : quartic_cubic(GEN g, long v)
     699             : {
     700       12131 :   GEN a = gel(g, 6), b = gel(g, 5), c = gel(g, 4);
     701       12131 :   GEN a3 = gdivgu(a,3);
     702       12131 :   return deg1pol(gmul2n(a3,2), gsub(gsqr(b),gmul2n(gmul(a3,c),3)), v);
     703             : }
     704             : 
     705             : static GEN
     706        6594 : quarticinv_pol(GEN IJ)
     707             : {
     708        6594 :   GEN I = gel(IJ,1), J = gel(IJ,2);
     709        6594 :   return mkpoln(4, gen_1, gen_0, gmulgs(I,-3), J);
     710             : }
     711             : static GEN
     712        2226 : quartic_H(GEN g, GEN *pT)
     713             : {
     714        2226 :   GEN a = gel(g, 6), b = gel(g, 5), c = gel(g, 4);
     715        2226 :   GEN IJ = quartic_IJ(g), I = gel(IJ, 1);
     716        2226 :   GEN ddh = quartic_hessiandd(g);
     717        2226 :   GEN ddg = deg2pol_shallow(gmulgu(a,12), gmulgu(b,6), gmulgu(c,2), 1);
     718        2226 :   *pT = quarticinv_pol(IJ);
     719        2226 :   return deg2pol_shallow(stoi(-8), gmul2n(ddg,2), gadd(ddh,gmul2n(I,3)),0);
     720             : }
     721             : 
     722             : static GEN
     723        4368 : quartic_disc(GEN q)
     724             : {
     725        4368 :   GEN IJ = quartic_IJ(q), I = gel(IJ,1), J = gel(IJ,2);
     726        4368 :   return gsub(gmul2n(gpowgs(I, 3), 2), gsqr(J));
     727             : }
     728             : 
     729             : static GEN
     730        6951 : quartic_minim_all(GEN F, GEN discF)
     731             : {
     732        6951 :   GEN IJ = quartic_IJ(F), I = gel(IJ,1), J = gel(IJ,2);
     733        6951 :   GEN g = Z_ppo(gcdii(I,J), gel(discF,1));
     734        6951 :   GEN plist = ZV_sort_uniq(shallowconcat(gel(absZ_factor(g),1), gel(discF,2)));
     735        6951 :   GEN W, C, PQ = hyperellminimalmodel(F, &C, plist);
     736        6951 :   GEN P = gel(PQ,1), Q = gel(PQ,2);
     737        6951 :   if (signe(Q)==0)
     738         812 :     W = mkvec2(P, C);
     739             :   else
     740        6139 :     W = mkvec2(ZX_add(ZX_shifti(P,2),ZX_sqr(Q)),mkvec2(shifti(gel(C,1),-1),gel(C,2)));
     741        6951 :   return W;
     742             : }
     743             : 
     744             : /*******************************************************************/
     745             : /*                         Cassels' pairing                        */
     746             : /*******************************************************************/
     747             : 
     748             : static GEN
     749        5957 : nfsqrt(GEN nf, GEN z)
     750             : {
     751        5957 :   long v = fetch_var_higher();
     752        5957 :   GEN R = nfroots(nf, deg2pol_shallow(gen_m1, gen_0, z, v));
     753        5957 :   delete_var();
     754        5957 :   return lg(R)==1 ? NULL: gel(R,1);
     755             : }
     756             : 
     757             : static GEN
     758        5957 : nfsqrt_safe(GEN nf, GEN z)
     759             : {
     760        5957 :   GEN r = nfsqrt(nf, z);
     761        5957 :   if (!r) pari_err_BUG("ellrank");
     762        5957 :   return r;
     763             : }
     764             : 
     765             : static GEN
     766        2226 : vecnfsqrtmod(GEN x, GEN P)
     767        8183 : { pari_APPLY_same(gmodulo(nfsqrt_safe(gel(x,i), P), gel(x,i))) }
     768             : 
     769             : static GEN
     770        2226 : enfsqrt(GEN T, GEN P)
     771             : {
     772        2226 :   GEN F = gel(ZX_factor(T),1);
     773        2226 :   return liftpol(chinese1(vecnfsqrtmod(F,P)));
     774             : }
     775             : 
     776             : /* Quartic q, at most quadratic g s.t. lc(g) > 0. There exist a real r s.t.
     777             :  * q(r) > 0 and g(r) != 0. Return sign(g(r)) */
     778             : static int
     779         476 : cassels_oo_solve_i(GEN q, GEN g)
     780             : {
     781         476 :   long dg = degpol(g);
     782             :   GEN D, a, b, c;
     783             : 
     784         476 :   if (dg == 0 || signe(leading_coeff(q)) > 0) return 1;
     785         217 :   if (signe(gel(q,2)) > 0) return signe(gel(g,2));
     786         189 :   c = gel(g,2); b = gel(g,3);
     787             :   /* g = bx+c, b>0, is negative on I=]-oo,-c/b[: if q has a root there,
     788             :    * then g(r) < 0. Else it has the sign of q(oo) < 0 on I */
     789         189 :   if (dg == 1) return ZX_sturmpart(q, mkvec2(mkmoo(), gdiv(negi(c), b)))? -1: 1;
     790         175 :   a = gel(g,4); D = subii(sqri(b), shifti(mulii(a,c), 2)); /* g = ax^2+bx+c */
     791         175 :   if (signe(D) <= 0) return 1; /* sign(g) = 1 is constant */
     792             :   /* Rescale q and g: x->(x - b)/2a; roots of new g are \pm sqrt(D) */
     793         161 :   q = ZX_translate(ZX_rescale(q, shifti(a,1)), negi(b));
     794             :   /* Now g has sign -1 in I=[-sqrt(D),sqrt(D)] and 1 elsewhere.
     795             :    * Check if q vanishes in I <=> Graeffe(q) vanishes on [0,D].
     796             :    * If so or if q(0) > 0 we take r in there; else r is outside of I */
     797         154 :   return (signe(gel(q,2)) > 0 ||
     798         315 :           ZX_sturmpart(ZX_graeffe(q), mkvec2(gen_0, D)))? -1: 1;
     799             : }
     800             : static int
     801         476 : cassels_oo_solve(GEN q, GEN g)
     802         476 : { pari_sp av = avma; return gc_int(av, cassels_oo_solve_i(q, g)); }
     803             : 
     804             : static GEN
     805       59199 : cassels_Qp_solve(GEN q, GEN gam, GEN p)
     806             : {
     807       59199 :   pari_sp av = avma;
     808       59199 :   GEN a = hyperell_local_solve(q, p);
     809       59199 :   GEN c = poleval(gam,a);
     810             :   long e;
     811       59199 :   if (!gequal0(c)) return c;
     812          70 :   for (e = 2; ; e++)
     813           0 :   {
     814          70 :     GEN b = gadd(a, powiu(p,e));
     815          70 :     if (Qp_issquare(poleval(q, b), p))
     816             :     {
     817          70 :       c = poleval(gam, b);
     818          70 :       if (!gequal0(c)) return gerepileupto(av,c);
     819             :     }
     820             :   }
     821             : }
     822             : 
     823             : static GEN
     824        4452 : to_ZX(GEN a, long v) { return typ(a)==t_INT? scalarpol_shallow(a,v): a; }
     825             : 
     826             : static GEN
     827        4368 : quartic_findunit(GEN D, GEN q)
     828             : {
     829        4368 :   GEN T = quarticinv_pol(quartic_IJ(q));
     830             :   while(1)
     831        1085 :   {
     832        5453 :     pari_sp av = avma;
     833        5453 :     GEN z = quartic_cubic(q,0);
     834        5453 :     if (signe(QXQ_norm(z,T)))
     835        4368 :       return absequalii(quartic_disc(q), D)? q: ZX_shifti(q, 2);
     836        1085 :     set_avma(av);
     837        1085 :     q = ZX_translate(RgX_recip(q), gen_1);
     838             :   }
     839             : }
     840             : 
     841             : /* Crude implementation of an algorithm by Tom Fisher
     842             :  * On binary quartics and the Cassels-Tate pairing
     843             :  * https://www.dpmms.cam.ac.uk/~taf1000/papers/bq-ctp.pdf */
     844             : 
     845             : /* FD = primes | 2*3*5*7*D, q1,q2,q3 have discriminant D */
     846             : static long
     847        2226 : casselspairing(GEN FD, GEN q1, GEN q2, GEN q3)
     848             : {
     849        2226 :   pari_sp av = avma;
     850        2226 :   GEN T, H = quartic_H(q1, &T);
     851        2226 :   GEN z1 = quartic_cubic(q1, 0);
     852        2226 :   GEN z2 = quartic_cubic(q2, 0);
     853        2226 :   GEN z3 = quartic_cubic(q3, 0);
     854        2226 :   GEN m = to_ZX(enfsqrt(T, QXQ_mul(QX_mul(z1,z2),z3,T)), 0);
     855        2226 :   GEN Hm = RgXQ_mul(QXQ_div(m, z1, T), H, T); /* deg(Hm) >= 2 */
     856        2226 :   GEN gam = to_ZX(Q_primpart(gel(Hm,4)),1);
     857        2226 :   GEN a = leading_coeff(q2), Fa = gel(absZ_factor(a),1);
     858        2226 :   GEN F = ZV_sort_uniq(shallowconcat1(mkvec2(Fa, FD)));
     859        2226 :   long i, e = 0, lF = lg(F);
     860        2226 :   if (signe(a) <= 0)
     861             :   {
     862         476 :     if (signe(leading_coeff(gam)) < 0) gam = ZX_neg(gam);
     863         476 :     if (cassels_oo_solve(q1, gam) < 0) e = 1;
     864             :   }
     865       61425 :   for (i = 1; i < lF; i++)
     866             :   {
     867       59199 :     GEN p = gel(F, i);
     868       59199 :     GEN c = cassels_Qp_solve(q1, gam, p);
     869       59199 :     if (hilbert(c, a, p) < 0) e = !e;
     870             :   }
     871        2226 :   return gc_long(av,e);
     872             : }
     873             : 
     874             : static GEN
     875         294 : matcassels(GEN F, GEN M)
     876             : {
     877         294 :   long i, j, n = lg(M)-1;
     878         294 :   GEN C = zero_F2m_copy(n,n);
     879         294 :   pari_sp av = avma;
     880        1862 :   for (i = 1; i <= n; i++)
     881             :   {
     882        1568 :     GEN Mii = gcoeff(M,i,i);
     883        1568 :     if (isintzero(Mii)) continue;
     884        4368 :     for (j = 1; j < i; j++)
     885             :     {
     886        3395 :       GEN Mjj = gcoeff(M,j,j);
     887        3395 :       if (!isintzero(Mjj) && casselspairing(F, Mii, Mjj, gcoeff(M,i,j)))
     888         336 :       { F2m_set(C,i,j); F2m_set(C,j,i); }
     889             :     }
     890             :   }
     891         294 :   return gc_const(av, C);
     892             : }
     893             : 
     894             : /*******************************************************************/
     895             : /*                         ELLRANK                                 */
     896             : /*******************************************************************/
     897             : /* This section is a port by Bill Allombert of ellQ.gp by Denis Simon
     898             :  * Copyright (C) 2019 Denis Simon
     899             :  * Distributed under the terms of the GNU General Public License (GPL) */
     900             : 
     901             : /* \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */
     902             : /*    FUNCTIONS FOR POLYNOMIALS                \\ */
     903             : /* \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */
     904             : 
     905             : static GEN
     906        1687 : ell2pol(GEN ell)
     907        1687 : { return mkpoln(4, gen_1, ell_get_a2(ell), ell_get_a4(ell), ell_get_a6(ell)); }
     908             : 
     909             : /* find point (x:y:z) on y^2 = pol, return [x,z]~ and set *py = y */
     910             : static GEN
     911        4319 : projratpointxz(GEN pol, long lim, GEN *py)
     912             : {
     913             :   pari_timer ti;
     914             :   GEN P;
     915        4319 :   if (issquareall(leading_coeff(pol), py)) return mkcol2(gen_1, gen_0);
     916        4270 :   if (DEBUGLEVEL) timer_start(&ti);
     917        4270 :   P = hyperellratpoints(pol, stoi(lim), 1);
     918        4270 :   if (DEBUGLEVEL) timer_printf(&ti,"hyperellratpoints(%ld)",lim);
     919        4270 :   if (lg(P) == 1) return NULL;
     920        1736 :   P = gel(P,1); *py = gel(P,2); return mkcol2(gel(P,1), gen_1);
     921             : }
     922             : 
     923             : /* P a list of integers (actually primes) one of which divides x; return
     924             :  * the first one */
     925             : static GEN
     926         623 : first_divisor(GEN x, GEN P)
     927             : {
     928         623 :   long i, n = lg(P);
     929         623 :   for (i = 1; i < n; i++)
     930         623 :     if (dvdii(x, gel(P,i))) return gel(P,i);
     931           0 :   return gel(P,i);
     932             : }
     933             : 
     934             : /* find point (x:y:z) on y^2 = pol, return [x,z]~ and set *py = y */
     935             : static GEN
     936        1659 : projratpointxz2(GEN pol, long lim, GEN *py)
     937             : {
     938        1659 :   pari_sp av = avma;
     939        1659 :   GEN list = mkvec(mkvec4(pol, matid(2), gen_1, gen_1));
     940             :   long i, j, c;
     941             : 
     942        4431 :   for (i = 1, c = 1; i < lg(list); i++,c++)
     943             :   {
     944        2807 :     GEN K, k, ff, co, p, M, C, r, pol, L = gel(list, i);
     945             :     long lr;
     946             : 
     947        2807 :     list = vecsplice(list, i); i--;
     948        2807 :     pol = Q_primitive_part(gel(L,1), &K);
     949        2807 :     M = gel(L,2);
     950        2807 :     K = K? mulii(gel(L,3), K): gel(L,3);
     951        2807 :     C = gel(L,4);
     952        2807 :     if (Z_issquareall(K, &k))
     953             :     {
     954             :       GEN xz, y, aux, U;
     955        3059 :       if (c==1) continue;
     956         910 :       pol = ZX_hyperellred(pol, &U);
     957         910 :       if (DEBUGLEVEL>1) err_printf("  reduced quartic(%ld): Y^2 = %Ps\n", i, pol);
     958         910 :       xz = projratpointxz(pol, lim, &y); if (!xz) continue;
     959          35 :       *py = gmul(y, mulii(C, k));
     960          35 :       aux = RgM_RgC_mul(ZM2_mul(M, U), xz);
     961          70 :       if (gequal0(gel(aux, 2))) return mkcol2(gel(aux,1), gen_0);
     962          35 :       *py = gdiv(*py, gpowgs(gel(aux, 2), degpol(pol)>>1));
     963          35 :       return mkcol2(gdiv(gel(aux, 1), gel(aux, 2)), gen_1);
     964             :     }
     965         623 :     ff = Z_factor(K); co = core2(mkvec2(K, ff)); K = gel(co,1); /* > 1 */
     966         623 :     p = first_divisor(K, gel(ff,1));
     967         623 :     K = diviiexact(K, p);
     968         623 :     C = mulii(mulii(C, gel(co,2)), p);
     969             :     /* root at infinity */
     970         623 :     if (dvdii(leading_coeff(pol), p))
     971             :     {
     972         329 :       GEN U = mkmat2(gel(M,1), ZC_Z_mul(gel(M,2), p));
     973         329 :       if (equali1(content(U)))
     974             :       {
     975         147 :         GEN t = ZX_rescale(pol, p);
     976         147 :         list = vec_append(list, mkvec4(ZX_Z_divexact(t,p), U, K, C));
     977             :       }
     978             :     }
     979         623 :     r = FpC_center(FpX_roots(pol, p), p, shifti(p,-1)); lr = lg(r);
     980        1701 :     for (j = 1; j < lr; j++)
     981             :     {
     982        1078 :       GEN U = ZM2_mul(M, mkmat22(p, gel(r, j), gen_0, gen_1));
     983        1078 :       if (equali1(content(U)))
     984             :       {
     985        1050 :         GEN t = ZX_unscale_div(ZX_translate(pol, gel(r,j)), p);
     986        1050 :         list = vec_append(list, mkvec4(t, U, K, C));
     987             :       }
     988             :     }
     989         623 :     if (gc_needed(av, 1)) gerepileall(av, 2, &pol, &list);
     990             :   }
     991        1624 :   return NULL;
     992             : }
     993             : 
     994             : static GEN
     995        8267 : polrootsmodpn(GEN pol, GEN p)
     996             : {
     997        8267 :   pari_sp av = avma, av2;
     998        8267 :   long j, l, i = 1, vd = Z_pval(ZX_disc(pol), p);
     999             :   GEN v, r, P;
    1000             : 
    1001        8267 :   if (!vd) { set_avma(av); retmkvec(zerovec(2)); }
    1002        6041 :   pol = Q_primpart(pol);
    1003        6041 :   P = gpowers0(p, vd-1, p); av2 = avma;
    1004        6041 :   v = FpX_roots(pol, p); l = lg(v);
    1005       17605 :   for (j = 1; j < l; j++) gel(v,j) = mkvec2(gel(v,j), gen_1);
    1006       69531 :   while (i < lg(v))
    1007             :   {
    1008       63490 :     GEN pol2, pe, roe = gel(v, i), ro = gel(roe,1);
    1009       63490 :     long e = itou(gel(roe,2));
    1010             : 
    1011       64470 :     if (e >= vd) { i++; continue; }
    1012       48832 :     pe = gel(P, e);
    1013       48832 :     (void)ZX_pvalrem(ZX_affine(pol, pe, ro), p, &pol2);
    1014       48832 :     r = FpX_roots(pol2, p); l = lg(r);
    1015       48832 :     if (l == 1) { i++; continue; }
    1016       99778 :     for (j = 1; j < l; j++)
    1017       51926 :       gel(r, j) = mkvec2(addii(ro, mulii(pe, gel(r, j))), utoi(e + 1));
    1018             :     /* roots with higher precision = ro + r*p^(e+1) */
    1019       47852 :     if (l > 2) v = shallowconcat(v, vecslice(r, 2, l-1));
    1020       47852 :     gel(v, i) = gel(r, 1);
    1021       47852 :     if (gc_needed(av2, 1)) gerepileall(av2, 1, &v);
    1022             :   }
    1023        6041 :   if (lg(v) == 1) { set_avma(av); retmkvec(zerovec(2)); }
    1024        6041 :   return gerepilecopy(av, v);
    1025             : }
    1026             : 
    1027             : /* \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */
    1028             : /*    FUNCTIONS FOR LOCAL COMPUTATIONS         \\ */
    1029             : /* \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */
    1030             : 
    1031             : /* z is integral; sprk true (pr | 2) [t_VEC] or modpr structure (pr | p odd)
    1032             :  * [t_COL] */
    1033             : static GEN
    1034     1595855 : kpmodsquares1(GEN nf, GEN z, GEN sprk)
    1035             : {
    1036     1595855 :   GEN modpr = (typ(sprk) == t_VEC)? gmael(sprk, 4, 1): sprk;
    1037     1595855 :   GEN pr = modpr_get_pr(modpr), p = pr_get_p(pr);
    1038     1595855 :   long v = nfvalrem(nf, z, pr, &z);
    1039     1595855 :   if (equaliu(p, 2))
    1040             :   {
    1041             :     GEN c;
    1042      156527 :     z = to_principal_unit(nf, z, pr, sprk); /* = 1 mod pr */
    1043      156527 :     c = ZV_to_Flv(sprk_log_prk1(nf, z, sprk), 2);
    1044      156527 :     return vecsmall_prepend(c, odd(v));
    1045             :   }
    1046             :   else
    1047             :   {
    1048     1439328 :     GEN T = modpr_get_T(modpr);
    1049     1439328 :     long c = !Fq_issquare(nf_to_Fq(nf, z, modpr), T, p);
    1050     1439328 :     return mkvecsmall2(odd(v), c);
    1051             :   }
    1052             : }
    1053             : 
    1054             : static GEN
    1055      782272 : kpmodsquares(GEN vnf, GEN z, GEN PP)
    1056             : {
    1057      782272 :   pari_sp av = avma;
    1058      782272 :   long i, j, l = lg(vnf);
    1059      782272 :   GEN dz, vchar = cgetg(l, t_VEC);
    1060      782272 :   z = Q_remove_denom(z, &dz); if (dz) z = ZX_Z_mul(z, dz);
    1061     1858047 :   for (i = 1; i < l; i++)
    1062             :   {
    1063     1075775 :     GEN nf = gel(vnf, i), pp = gel(PP, i);
    1064     1075775 :     GEN kp, delta = ZX_rem(z, nf_get_pol(nf));
    1065     1075775 :     long lp = lg(pp);
    1066     1075775 :     kp = cgetg(lp, t_VEC);
    1067     2671630 :     for (j = 1; j < lp; j++) gel(kp, j) = kpmodsquares1(nf, delta, gel(pp,j));
    1068     1075775 :     gel(vchar, i) = shallowconcat1(kp);
    1069             :   }
    1070      782272 :   return gerepileuptoleaf(av, shallowconcat1(vchar));
    1071             : }
    1072             : 
    1073             : static GEN
    1074       16534 : veckpmodsquares(GEN x, GEN vnf, GEN PP)
    1075      777119 : { pari_APPLY_type(t_MAT, kpmodsquares(vnf, gel(x, i), PP)) }
    1076             : 
    1077             : /* \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */
    1078             : /*    GENERIC FUNCTIONS FOR ELLIPTIC CURVES    \\ */
    1079             : /* \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */
    1080             : 
    1081             : static GEN
    1082        3899 : ellabs(GEN P)
    1083        3899 : { return ell_is_inf(P) ? P: mkvec2(gel(P,1), Q_abs_shallow(gel(P,2))); }
    1084             : static GEN
    1085        4648 : vecellabs(GEN x) { pari_APPLY_same(ellabs(gel(x,i))) }
    1086             : 
    1087             : /* y^2 = x^3 + K a2 x + K^2 a4 x + K^3 a6 */
    1088             : static GEN
    1089         805 : elltwistequation(GEN ell, GEN K)
    1090             : {
    1091             :   GEN K2, a2, a4, a6;
    1092         805 :   if (!K || equali1(K)) return ell;
    1093          84 :   K2 = sqri(K);
    1094          84 :   a2 = mulii(ell_get_a2(ell), K);
    1095          84 :   a4 = mulii(ell_get_a4(ell), K2);
    1096          84 :   a6 = mulii(ell_get_a6(ell), mulii(K, K2));
    1097          84 :   return ellinit(mkvec5(gen_0, a2, gen_0, a4, a6), NULL, DEFAULTPREC);
    1098             : }
    1099             : 
    1100             : /* P=[x,y] a point on Ky^2 =  pol(x); [Kx, K^2y] point on Y^2 = K^3pol(X/K) */
    1101             : static GEN
    1102         224 : elltwistpoint(GEN P, GEN K, GEN K2)
    1103             : {
    1104         224 :   if (ell_is_inf(P)) return ellinf();
    1105         224 :   return mkvec2(gmul(gel(P,1), K), gmul(gel(P,2), K2));
    1106             : }
    1107             : 
    1108             : static GEN
    1109        1603 : elltwistpoints(GEN x, GEN K)
    1110             : {
    1111             :   GEN K2;
    1112        1603 :   if (!K || gequal1(K)) return x;
    1113         126 :   K2 = gsqr(K);
    1114         350 :   pari_APPLY_same(elltwistpoint(gel(x,i), K, K2))
    1115             : }
    1116             : 
    1117             : /* \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */
    1118             : /*    FUNCTIONS FOR NUMBER FIELDS              \\ */
    1119             : /* \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */
    1120             : 
    1121             : /* return a set S2 of prime ideals disjoint from S such that
    1122             :  * Cl_{S \cup S2}(K) has no p-torsion */
    1123             : static GEN
    1124         434 : bestS(GEN bnf,GEN S, ulong p)
    1125             : {
    1126         434 :   GEN v, S2, h = bnf_get_no(bnf), cyc = bnf_get_cyc(bnf);
    1127             :   long i, lS2;
    1128             :   ulong l, vD;
    1129             :   forprime_t P;
    1130             : 
    1131         434 :   if (!dvdiu(h, p)) return cgetg(1,t_VEC);
    1132         315 :   if (!S)
    1133             :   {
    1134           0 :     v = diagonal_shallow(cyc);
    1135           0 :     vD = Z_lval(h, p);
    1136             :   }
    1137             :   else
    1138             :   {
    1139         315 :     long lS = lg(S);
    1140         315 :     v = cgetg(lS,t_MAT);
    1141         994 :     for (i = 1; i < lS; i++) gel(v,i) = isprincipal(bnf, gel(S,i));
    1142         315 :     v = ZM_hnfmodid(v, cyc);
    1143         315 :     vD = Z_lval(ZM_det(v), p); if (!vD) return cgetg(1, t_VEC);
    1144             :   }
    1145         301 :   S2 = cgetg(vD+2, t_VEC); lS2 = 1;
    1146         301 :   u_forprime_init(&P,2,ULONG_MAX);
    1147        2569 :   while ((l = u_forprime_next(&P)))
    1148             :   {
    1149        2569 :     pari_sp av = avma;
    1150        2569 :     GEN w, Sl = idealprimedec(bnf, utoi(l));
    1151        2569 :     long nSl = lg(Sl)-1;
    1152             :     ulong vDl;
    1153        4074 :     for (i = 1; i < nSl; i++) /* remove one prime ideal */
    1154             :     {
    1155        1806 :       w = ZM_hnf(shallowconcat(v, isprincipal(bnf, gel(Sl,i))));
    1156        1806 :       vDl = Z_lval(ZM_det(w), p);
    1157        1806 :       if (vDl < vD)
    1158             :       {
    1159        1239 :         gel(S2,lS2++) = gel(Sl,i);
    1160        1239 :         vD = vDl; v = w; av = avma;
    1161        1239 :         if (!vD) { setlg(S2, lS2); return S2; }
    1162             :       }
    1163             :     }
    1164        2268 :     set_avma(av);
    1165             :   }
    1166             :   return NULL;/*LCOV_EXCL_LINE*/
    1167             : }
    1168             : 
    1169             : static GEN
    1170         868 : nfC_prV_val(GEN nf, GEN G, GEN P)
    1171             : {
    1172         868 :   long i, j, lG = lg(G), lP = lg(P);
    1173         868 :   GEN M = cgetg(lG, t_MAT);
    1174        5362 :   for (i = 1; i < lG; i++)
    1175             :   {
    1176        4494 :     GEN V = cgetg(lP, t_COL);
    1177       19285 :     for (j = 1; j < lP; j++)
    1178       14791 :       gel(V,j) = gpnfvalrem(nf, gel(G,i), gel(P,j), NULL);
    1179        4494 :     gel(M,i) = V;
    1180             :   }
    1181         868 :   return M;
    1182             : }
    1183             : 
    1184             : static GEN
    1185        5558 : _factorbackmod(GEN nf, GEN g, GEN e, ulong p)
    1186             : {
    1187        5558 :   GEN y = nffactorback(nf, g, e), den;
    1188        5558 :   GEN z = nfmul(nf, y, nfsqr(nf, idealredmodpower(nf, y, p, 0)));
    1189        5558 :   z = Q_remove_denom(z, &den);
    1190        5558 :   if (den)
    1191             :   {
    1192          35 :     if (p != 2) den = powiu(den, p-1);
    1193          35 :     z = gmul(z, den);
    1194             :   }
    1195        5558 :   return z;
    1196             : }
    1197             : static GEN
    1198         434 : nfV_factorbackmod(GEN nf, GEN x, ulong p)
    1199             : {
    1200         434 :   long i, l = lg(x);
    1201         434 :   GEN v = cgetg(l, t_VEC);
    1202        3752 :   for (i = 1; i < l; i++)
    1203             :   {
    1204        3318 :     GEN y = gel(x,i), g = gel(y,1), e = gel(y,2);
    1205        3318 :     gel(v,i) = _factorbackmod(nf, g, ZV_to_Flv(e,p), p);
    1206             :   }
    1207         434 :   return v;
    1208             : }
    1209             : static GEN
    1210         434 : nfV_zm_factorback(GEN nf, GEN G, GEN x, long p)
    1211        2674 : { pari_APPLY_type(t_VEC, _factorbackmod(nf, G, gel(x,i), p)) }
    1212             : 
    1213             : static GEN
    1214         434 : prV_ZM_factorback(GEN nf, GEN S, GEN x)
    1215        2674 : { pari_APPLY_type(t_VEC,idealfactorback(nf, S, gel(x,i), 0)) }
    1216             : 
    1217             : /* shortcut for bnf = Q and p = 2 */
    1218             : static GEN
    1219         784 : bnfselmerQ(GEN S)
    1220             : {
    1221         784 :   GEN g = vec_prepend(prV_primes(S), gen_m1), e;
    1222         784 :   long n = lg(S)-1;
    1223         784 :   e = n? shallowconcat(zerocol(n), matid(n)): zeromat(0, 1);
    1224         784 :   return mkvec3(g, e, const_vec(n + 1, gen_1));
    1225             : }
    1226             : 
    1227             : static GEN
    1228         434 : bnfselmer(GEN bnf, GEN S)
    1229             : {
    1230         434 :   const long p = 2;
    1231         434 :   pari_sp av = avma;
    1232         434 :   GEN nf = bnf_get_nf(bnf), S2, S3, e, f, e2, kerval, LS2gen, LS2fu, LS2all;
    1233         434 :   long n = lg(S)-1, n3, n2all, r;
    1234             : 
    1235         434 :   S2 = bestS(bnf, S, p);
    1236         434 :   S3 = shallowconcat(S, S2);
    1237         434 :   LS2all = nfV_factorbackmod(nf, gel(bnfunits(bnf, S3), 1), p);
    1238         434 :   n3 = lg(S3)-1; n2all = lg(LS2all)-1;
    1239         434 :   LS2gen = vecslice(LS2all,1,n3);
    1240         434 :   LS2fu  = vecslice(LS2all,n3+1, n2all-1);
    1241         434 :   e2 = nfC_prV_val(nf, LS2gen, S2);
    1242         434 :   kerval = Flm_ker(ZM_to_Flm(e2, p), p);
    1243         434 :   LS2gen = nfV_zm_factorback(nf, LS2gen, kerval, p);
    1244         434 :   e = nfC_prV_val(nf, LS2gen, S);
    1245         434 :   e2 = ZM_divexactu(ZM_zm_mul(e2, kerval), p);
    1246         434 :   f = prV_ZM_factorback(nf, S2, e2);
    1247         434 :   LS2gen = shallowconcat(LS2fu, LS2gen);
    1248         434 :   LS2gen = nfV_to_scalar_or_alg(nf, vec_prepend(LS2gen, bnf_get_tuU(bnf)));
    1249         434 :   r = n2all - n3;
    1250         434 :   e = shallowconcat(zeromat(n, r), e);
    1251         434 :   f = shallowconcat(const_vec(r, gen_1), f);
    1252         434 :   return gerepilecopy(av, mkvec3(LS2gen,e,f));
    1253             : }
    1254             : 
    1255             : static GEN
    1256         238 : get_kerval(GEN nf, GEN S2, GEN LS2gen)
    1257             : {
    1258         238 :   long i, j, lS2 = lg(S2), l = lg(LS2gen);
    1259         238 :   GEN e = cgetg(l, t_MAT);
    1260        3521 :   for (i = 1; i < l; i++)
    1261             :   {
    1262        3283 :     GEN v = cgetg(lS2, t_VECSMALL);
    1263       25151 :     for (j=1; j < lS2; j++) v[j] = odd(idealval(nf, gel(LS2gen, i), gel(S2,j)));
    1264        3283 :     gel(e, i) = Flv_to_F2v(v);
    1265             :   }
    1266         238 :   return F2m_ker(e);
    1267             : }
    1268             : static GEN
    1269         238 : nf2selmer_quad(GEN nf, GEN S)
    1270             : {
    1271         238 :   pari_sp ltop = avma;
    1272         238 :   GEN D = nf_get_disc(nf), factD = nf_get_ramified_primes(nf);
    1273         238 :   GEN SlistQ = prV_primes(S), QS2gen, gen, Hlist, H, KerH, norms, LS2gen;
    1274             :   GEN chpol, Q, kerval, S2, G, e, f, b, c, bad;
    1275         238 :   long lS = lg(S), l, lHlist, i, j, k;
    1276             : 
    1277         238 :   QS2gen = vec_prepend(SlistQ, gen_m1);
    1278         238 :   bad = ZV_sort_uniq(shallowconcat(factD, SlistQ));
    1279         238 :   Hlist = ZV_search(bad, gen_2)? bad: vec_prepend(bad, gen_2);
    1280         238 :   l = lg(QS2gen);
    1281         238 :   lHlist = lg(Hlist);
    1282         238 :   H = cgetg(l, t_MAT);
    1283        2219 :   for (j = 1; j < l; j++)
    1284             :   {
    1285        1981 :     GEN v = cgetg(lHlist, t_VECSMALL);
    1286       29715 :     for (i = 1; i < lHlist; i++)
    1287       27734 :       v[i] = hilbertii(D, gel(QS2gen, j), gel(Hlist, i)) < 0;
    1288        1981 :     gel(H, j) = Flv_to_F2v(v);
    1289             :   }
    1290         238 :   KerH = F2m_ker(H); l = lg(KerH);
    1291         238 :   norms = cgetg(l, t_VEC);
    1292        1344 :   for (i = 1; i < l; i++)
    1293        1106 :     gel(norms, i) = factorback2(QS2gen, F2c_to_ZC(gel(KerH, i)));
    1294         238 :   LS2gen = cgetg(l, t_VEC);
    1295         238 :   chpol = QXQ_charpoly(gel(nf_get_zk(nf), 2), nf_get_pol(nf), 0);
    1296         238 :   b = gdivgu(negi(gel(chpol, 3)), 2);
    1297         238 :   c = gel(chpol, 2);
    1298         238 :   Q = mkmat3(mkcol3(gen_1, b, gen_0),
    1299             :              mkcol3(b, c, gen_0),
    1300             :              mkcol3(gen_0, gen_0, gen_0));
    1301        1344 :   for (k = 1; k < l; k++)
    1302             :   {
    1303             :     GEN sol;
    1304        1106 :     gcoeff(Q, 3, 3) = negi(gel(norms, k));
    1305        1106 :     sol = qfsolve(Q); /* must be solvable */
    1306        1106 :     sol = Q_primpart(mkcol2(gel(sol,1), gel(sol,2)));
    1307        1106 :     gel(LS2gen, k) = basistoalg(nf, sol);
    1308             :   }
    1309         238 :   if (equalis(D, -4)) G = bad;
    1310             :   else
    1311             :   {
    1312         238 :     G = vecsplice(bad, ZV_search(bad, gel(factD, lg(factD)-1)));
    1313         238 :     G = vec_prepend(G, gen_m1);
    1314             :   }
    1315         238 :   LS2gen = shallowconcat(G, LS2gen);
    1316         238 :   l = lg(SlistQ); S2 = cgetg(l, t_VEC);
    1317         238 :   if (l > 1)
    1318             :   {
    1319        1974 :     for (i = 1; i < l; i++) gel(S2, i) = idealprimedec(nf, gel(SlistQ, i));
    1320         231 :     S2 = setminus(shallowconcat1(S2), S);
    1321             :   }
    1322         238 :   kerval = get_kerval(nf, S2, LS2gen); l = lg(kerval);
    1323         238 :   gen = cgetg(l, t_VEC);
    1324         238 :   e = cgetg(l, t_MAT);
    1325         238 :   f = cgetg(l, t_VEC);
    1326        2576 :   for (i = 1; i < l; i++)
    1327             :   {
    1328        2338 :     GEN id, ei, x = nffactorback(nf, LS2gen, F2v_to_Flv(gel(kerval, i)));
    1329        2338 :     gel(e,i) = ei = cgetg(lS, t_COL);
    1330       32704 :     for (j = 1; j < lS; j++) gel(ei,j) = stoi(idealval(nf, x, gel(S,j)));
    1331        2338 :     id = idealdiv(nf, x, idealfactorback(nf, S, gel(e,i), 0));
    1332        2338 :     if (!idealispower(nf, id, 2, &gel(f,i))) pari_err_BUG("nf2selmer_quad");
    1333        2338 :     gel(gen, i) = nf_to_scalar_or_alg(nf, x);
    1334             :   }
    1335         238 :   return gerepilecopy(ltop, mkvec3(gen, e, f));
    1336             : }
    1337             : 
    1338             : static GEN
    1339         833 : makevbnf(GEN ell, long prec)
    1340             : {
    1341         833 :   GEN vbnf, P = gel(ZX_factor(ell2pol(ell)), 1);
    1342         833 :   long k, l = lg(P);
    1343         833 :   vbnf = cgetg(l, t_VEC);
    1344        2282 :   for (k = 1; k < l; k++)
    1345             :   {
    1346        1449 :     GEN t = gel(P,k);
    1347        1449 :     gel(vbnf,k) = degpol(t) <= 2? nfinit(t, prec): Buchall(t, nf_FORCE, prec);
    1348             :   }
    1349         833 :   return vbnf;
    1350             : }
    1351             : 
    1352             : static GEN
    1353         854 : kernorm(GEN G, GEN S, GEN pol, GEN signs)
    1354             : {
    1355         854 :   long i, j, lS = lg(S), lG = lg(G), lv = signs? lS+2: lS+1;
    1356         854 :   GEN M = cgetg(lG, t_MAT);
    1357       13650 :   for (j = 1; j < lG; j++)
    1358             :   {
    1359       12796 :     GEN v, N = QXQ_norm(gel(G,j), pol);
    1360       12796 :     gel(M, j) = v = cgetg(lv, t_VECSMALL);
    1361       12796 :     v[1] = gsigne(N) < 0;
    1362      181097 :     for (i = 1; i < lS; i++) v[i+1] = odd(Q_pvalrem(N, gel(S,i), &N));
    1363       12796 :     if (signs) v[i+1] = signs[j];
    1364             :   }
    1365         854 :   return Flm_ker(M, 2);
    1366             : }
    1367             : 
    1368             : /* \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */
    1369             : /*    FUNCTIONS FOR 2-DESCENT                  \\ */
    1370             : /* \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */
    1371             : /* vector of t_VEC; return total number of entries */
    1372             : static long
    1373        8267 : RgVV_nb(GEN v)
    1374             : {
    1375        8267 :   long i, l = lg(v), n = 0;
    1376       23968 :   for (i = 1; i < l; i++) n += lg(gel(v,i)) - 1;
    1377        8267 :   return n;
    1378             : }
    1379             : 
    1380             : /* return an Fp basis */
    1381             : static GEN
    1382        8267 : elllocalimage(GEN pol, GEN K, GEN vnf, GEN p, GEN pp, GEN pts)
    1383             : {
    1384        8267 :   long n, p2 = RgVV_nb(pp), prank = equaliu(p, 2)? p2: p2 - 1;
    1385        8267 :   GEN R = polrootsmodpn(pol, p), bound = addiu(p, 6);
    1386             : 
    1387        8267 :   for(n = 1;; n++)
    1388       21687 :   {
    1389             :     pari_sp btop;
    1390             :     GEN x, y2, d;
    1391       29954 :     pts = Flm_image(pts, 2); if (lg(pts)-1 == prank) break;
    1392       21687 :     if ((n & 0xf) == 0) bound = mulii(bound, p);
    1393       21687 :     btop = avma;
    1394             :     do
    1395             :     {
    1396      105199 :       GEN r = gel(R, random_Fl(lg(R)-1)+1);
    1397      105199 :       long pprec = random_Fl(itou(gel(r, 2)) + 3) - 2; /* >= -2 */
    1398      105199 :       set_avma(btop);
    1399      105199 :       x = gadd(gel(r, 1), gmul(powis(p, pprec), randomi(bound)));
    1400      105199 :       y2 = gmul(K, poleval(pol, x));
    1401      105199 :     } while (gequal0(y2) || !Qp_issquare(y2, p));
    1402       21687 :     d = deg1pol_shallow(negi(K), gmul(K, x), 0);
    1403       21687 :     pts = vec_append(pts, kpmodsquares(vnf, d, pp));
    1404             :   }
    1405        8267 :   return pts;
    1406             : }
    1407             : 
    1408             : static GEN
    1409         854 : ellLS2image(GEN pol, GEN vP, GEN K, GEN vpol, GEN vcrt)
    1410             : {
    1411         854 :   long i, l = lg(vP);
    1412             :   GEN v;
    1413             : 
    1414         854 :   if (l == 1) return cgetg(1, t_VEC);
    1415         721 :   v = cgetg(l, t_VEC);
    1416       59885 :   for (i = 1; i < l; i++)
    1417             :   {
    1418       59164 :     GEN P = gel(vP, i), x, t;
    1419       59164 :     if (ell_is_inf(P)) { gel(v, i) = gen_1; continue; }
    1420       59164 :     x = gel(P,1);
    1421       59164 :     t = deg1pol_shallow(negi(K), gmul(K, x), 0);
    1422       59164 :     if (gequal0(gel(P,2)))
    1423             :     { /* 2-torsion, x now integer and a root of exactly one linear vpol[k]=T */
    1424         560 :       long k, lp = lg(vpol);
    1425             :       GEN a;
    1426         994 :       for (k = 1; k < lp; k++)
    1427             :       {
    1428         994 :         GEN T = gel(vpol,k), z = gel(T,2);
    1429         994 :         if (absequalii(x, z) && signe(x) == -signe(z)) break; /* T = X-x */
    1430             :       }
    1431         560 :       a = ZX_Z_eval(ZX_deriv(pol), x);
    1432         560 :       t = gadd(a, gmul(gel(vcrt,k), gsub(t, a))); /* a mod T, t mod pol/T*/
    1433             :     }
    1434       59164 :     gel(v, i) = t;
    1435             :   }
    1436         721 :   return v;
    1437             : }
    1438             : 
    1439             : /* find small points on ell; 2-torsion points must be returned first */
    1440             : static GEN
    1441         854 : ellsearchtrivialpoints(GEN ell, GEN lim, GEN help)
    1442             : {
    1443         854 :   pari_sp av = avma;
    1444         854 :   GEN tors2 = gel(elltors_psylow(ell,2),3);
    1445         854 :   GEN triv = lim ? ellratpoints(ell, lim, 0): cgetg(1,t_VEC);
    1446         854 :   if (help) triv = shallowconcat(triv, help);
    1447         854 :   return gerepilecopy(av, shallowconcat(tors2, triv));
    1448             : }
    1449             : 
    1450             : GEN
    1451          63 : ellrankinit(GEN ell, long prec)
    1452             : {
    1453          63 :   pari_sp av = avma;
    1454             :   GEN urst;
    1455          63 :   checkell_Q(ell); ell = ellminimalbmodel(ell, &urst);
    1456          63 :   return gerepilecopy(av, mkvec3(ell, urst, makevbnf(ell, prec)));
    1457             : }
    1458             : 
    1459             : INLINE GEN
    1460      251902 : ZV_isneg(GEN x) { pari_APPLY_long(signe(gel(x,i)) < 0) }
    1461             : 
    1462             : static GEN
    1463       71960 : selmersign(GEN x, GEN vpol, GEN L)
    1464      160846 : { pari_APPLY_same(ZV_isneg(nfeltsign(gel(x, i), RgX_rem(L, gel(vpol, i)), NULL))) }
    1465             : 
    1466             : static GEN
    1467        1708 : matselmersign(GEN vnf, GEN vpol, GEN x)
    1468       73668 : { pari_APPLY_type(t_MAT, shallowconcat1(selmersign(vnf, vpol, gel(x, i)))) }
    1469             : 
    1470             : static GEN
    1471      115038 : _trace(GEN z, GEN T)
    1472             : {
    1473      115038 :   long n = degpol(T)-1;
    1474      115038 :   if (degpol(z) < n) return gen_0;
    1475       66507 :   return gdiv(gel(z, 2+n), gel(T, 3+n));
    1476             : }
    1477             : static GEN
    1478       19173 : tracematrix(GEN zc, GEN b, GEN T)
    1479             : {
    1480             :   long i, j;
    1481       19173 :   GEN q = cgetg(4, t_MAT);
    1482       76692 :   for (j = 1; j <= 3; j++) gel(q,j) = cgetg(4, t_COL);
    1483       76692 :   for (j = 1; j <= 3; j++)
    1484             :   {
    1485      115038 :     for (i = 1; i < j; i++) gcoeff(q,i,j) = gcoeff(q,j,i) =
    1486       57519 :       _trace(QXQ_mul(zc, QXQ_mul(gel(b,i), gel(b,j), T), T), T);
    1487       57519 :     gcoeff(q,i,i) = _trace(QXQ_mul(zc, QXQ_sqr(gel(b,i), T), T), T);
    1488             :   }
    1489       19173 :   return q;
    1490             : }
    1491             : 
    1492             : static GEN
    1493       20853 : RgXV_cxeval(GEN x, GEN r, GEN ri)
    1494       83412 : { pari_APPLY_same(RgX_cxeval(gel(x,i), r, ri)) }
    1495             : 
    1496             : static GEN
    1497        6951 : redquadric(GEN base, GEN q2, GEN pol, GEN zc)
    1498             : {
    1499        6951 :   long i, l, prec = nbits2prec(2*gexpo(q2)) + 1;
    1500        6951 :   GEN s = NULL, R = roots(pol, prec);
    1501        6951 :   l = lg(R);
    1502       27804 :   for (i = 1; i < l; ++i)
    1503             :   {
    1504       20853 :     GEN r = gel(R,i), ri = gexpo(r) > 1? ginv(r): NULL;
    1505       20853 :     GEN b = RgXV_cxeval(base, r, ri), z = RgX_cxeval(zc, r, ri);
    1506       20853 :     GEN M = RgC_RgV_mulrealsym(RgV_Rg_mul(b, gabs(z, prec)), gconj(b));
    1507       20853 :     s = s? RgM_add(s, M): M;
    1508             :   }
    1509        6951 :   return lllgram(s);
    1510             : }
    1511             : 
    1512             : static GEN
    1513       19299 : RgX_homogenous_evaldeg(GEN P, GEN A, GEN B)
    1514             : {
    1515       19299 :   long i, d = degpol(P), e = lg(B)-1;
    1516       19299 :   GEN s = gmul(gel(P, d+2), gel(B,e-d));
    1517       60550 :   for (i = d-1; i >= 0; i--)
    1518       41251 :     s = gadd(gmul(s, A), gmul(gel(B,e-i), gel(P,i+2)));
    1519       19299 :   return s;
    1520             : }
    1521             : 
    1522             : static GEN
    1523        5271 : RgXV_homogenous_evaldeg(GEN x, GEN a, GEN b)
    1524       21084 : { pari_APPLY_same(RgX_homogenous_evaldeg(gel(x,i), a, b)) }
    1525             : 
    1526             : static void
    1527          42 : check_oncurve(GEN ell, GEN v)
    1528             : {
    1529          42 :   long i, l = lg(v);
    1530          49 :   for (i = 1; i < l; i++)
    1531             :   {
    1532           7 :     GEN P = gel(v, i);
    1533           7 :     if (lg(P) != 3 || !oncurve(ell,P)) pari_err_TYPE("ellrank",P);
    1534             :   }
    1535          42 : }
    1536             : 
    1537             : static GEN
    1538       18704 : basis(GEN nf, GEN x, GEN crt, GEN pol)
    1539             : {
    1540       18704 :   long i, l = lg(x);
    1541       18704 :   GEN b = cgetg(l, t_COL);
    1542       41636 :   for (i = 1; i < l; i++)
    1543             :   {
    1544       22932 :     GEN z = nf_to_scalar_or_alg(nf, gel(x, i));
    1545       22932 :     gel(b, i) = grem(gsub(z, gmul(crt, z)), pol); /* z mod T, 0 mod (pol/T) */
    1546             :   }
    1547       18704 :   return b;
    1548             : }
    1549             : 
    1550             : static GEN
    1551         854 : vecsmallbasis(GEN x, GEN vcrt, GEN pol)
    1552        2310 : { pari_APPLY_same(basis(gel(x,i), nf_get_zk(gel(x,i)), gel(vcrt,i), pol)) }
    1553             : 
    1554             : static GEN
    1555      265426 : ZC_shifti(GEN x, long n) { pari_APPLY_type(t_COL, shifti(gel(x,i), n)) }
    1556             : 
    1557             : /* true nf */
    1558             : static GEN
    1559       17248 : selmerbasis(GEN nf, GEN ek, GEN sqrtLS2, GEN factLS2,
    1560             :             GEN badprimes, GEN crt, GEN pol)
    1561             : {
    1562       17248 :   GEN sqrtzc = idealfactorback(nf, sqrtLS2, zv_neg(ek), 0);
    1563       17248 :   GEN E = ZC_shifti(ZM_zc_mul(factLS2, ek), -1);
    1564             : 
    1565       17248 :   if (ZV_equal0(E))
    1566       15471 :     sqrtzc = idealhnf_shallow(nf, sqrtzc);
    1567             :   else
    1568        1777 :     sqrtzc = idealmul(nf, sqrtzc, idealfactorback(nf, badprimes, ZC_neg(E), 0));
    1569       17248 :   return basis(nf, sqrtzc, crt, pol);
    1570             : }
    1571             : 
    1572         483 : static long randu(void) { return random_Fl(127) - 63; }
    1573             : static GEN
    1574         161 : randS(GEN b)
    1575             : {
    1576         483 :   return gadd(gmulgs(gel(b,1), randu()),
    1577         322 :               gadd(gmulgs(gel(b,2), randu()), gmulgs(gel(b,3), randu())));
    1578             : }
    1579             : 
    1580             : static GEN
    1581        6790 : liftselmerinit(GEN expo, GEN vnf, GEN sqrtLS2, GEN factLS2,
    1582             :                GEN badprimes, GEN vcrt, GEN pol)
    1583             : {
    1584        6790 :   long n = lg(vnf)-1, k, t;
    1585        6790 :   GEN b = cgetg(n+1, t_VEC);
    1586       24038 :   for (k = t = 1; k <= n; k++)
    1587             :   {
    1588       17248 :     GEN fak = gel(factLS2,k), ek;
    1589       17248 :     long m = lg(fak)-1;
    1590       17248 :     ek = vecslice(expo, t, t + m-1); t += m;
    1591       17248 :     gel(b,k) = selmerbasis(gel(vnf,k), ek, gel(sqrtLS2,k),
    1592       17248 :                            fak, gel(badprimes,k), gel(vcrt,k), pol);
    1593             :   }
    1594        6790 :   return shallowconcat1(b);
    1595             : }
    1596             : 
    1597             : static GEN
    1598        3486 : liftselmer_cover(GEN b, GEN expo, GEN LS2, GEN pol, GEN discF, GEN K)
    1599             : {
    1600             :   GEN P, Q, Q4, R, den, q0, q1, q2, xz, x, y, y2m, U, param, newb;
    1601             :   GEN ttheta, tttheta, zc, polprime;
    1602             :   GEN QM, zden;
    1603        3486 :   zc = RgXQV_factorback(LS2, expo, pol);
    1604        3486 :   if (typ(zc)==t_INT) zc = scalarpol(zc, varn(pol));
    1605        3486 :   ttheta = RgX_shift_shallow(pol,-2);
    1606        3486 :   tttheta = RgX_shift_shallow(pol, -1);
    1607        3486 :   polprime = ZX_deriv(pol);
    1608        3486 :   q2 = Q_primpart(tracematrix(zc, b, pol));
    1609        3486 :   U = redquadric(b, q2, pol, QXQ_div(zc, polprime, pol));
    1610        3486 :   q2 = qf_apply_RgM(q2, U);
    1611        3486 :   newb = RgV_RgM_mul(b, U);
    1612        3486 :   param = Q_primpart(qfparam(q2, qfsolve(q2), 1));
    1613        3486 :   param = RgM_to_RgXV_reverse(shallowtrans(param), 0);
    1614        3486 :   q1 = RgM_neg(tracematrix(QXQ_mul(zc, ttheta, pol), newb, pol));
    1615        3486 :   q1 = Q_remove_denom(qfeval(q1, param), &den);
    1616        3486 :   if (den) q1 = ZX_Z_mul(q1, den);
    1617        3486 :   if (!equali1(K)) q1 = RgX_Rg_mul(q1, K);
    1618        3486 :   QM = quartic_minim_all(q1, discF);
    1619        3486 :   q1 = gel(QM,1);
    1620        3486 :   zden = gmael(QM,2,1);
    1621        3486 :   Q = ZX_hyperellred(q1, &R);
    1622        3486 :   R = gmul(gmael(QM,2,2), R);
    1623        3486 :   if (DEBUGLEVEL>1) err_printf("  reduced quartic: Y^2 = %Ps\n", Q);
    1624        3486 :   xz = mkcol2(pol_x(0),gen_1);
    1625        3486 :   P = RgM_RgC_mul(R, xz); x = gel(P,1); y = gel(P,2);
    1626        3486 :   param = RgXV_homogenous_evaldeg(param, x, gpowers(y, 2));
    1627        3486 :   param = gmul(param, gdiv(den? mulii(K, den): K, zden));
    1628        3486 :   q0 = tracematrix(QXQ_mul(zc, tttheta, pol), newb, pol);
    1629        3486 :   x = gdiv(qfeval(q0, param), K);
    1630        3486 :   Q4 = gpowers(Q,4);
    1631        3486 :   y2m = gmul(RgX_homogenous_evaldeg(pol, x, Q4), Q);
    1632        3486 :   if (!issquareall(gdiv(y2m, K), &y))
    1633           0 :     pari_err_BUG("liftselmer_cover");
    1634        3486 :   y = gdiv(y, gel(Q4,2));
    1635        3486 :   P = mkvec2(gdiv(gmul(K,x),pol_xn(2,1)),gdiv(gmul(gsqr(K),y),pol_xn(3,1)));
    1636        3486 :   return mkvec2(Q,P);
    1637             : }
    1638             : 
    1639             : static GEN
    1640        3304 : liftselmer(GEN b, GEN expo, GEN sbase, GEN LS2, GEN pol, GEN discF, GEN K, long ntry, GEN *pt_Q)
    1641             : {
    1642        3304 :   pari_sp av = avma, av2;
    1643        3304 :   long t, lim = ntry * LIM1;
    1644             :   GEN ttheta, tttheta, z, polprime;
    1645             :   hashtable h;
    1646        3304 :   hash_init_GEN(&h, ntry, ZX_equal, 1);
    1647        3304 :   z = RgXQV_factorback(LS2, expo, pol);
    1648        3304 :   ttheta = RgX_shift_shallow(pol,-2);
    1649        3304 :   tttheta = RgX_shift_shallow(pol, -1);
    1650        3304 :   polprime = ZX_deriv(pol); av2 = avma;
    1651        3969 :   for (t = 1; t <= ntry; t++, set_avma(av2))
    1652             :   {
    1653             :     GEN P, Q, Qk, R, den, q0, q1, q2, xz, x, y, zz, zc, U, param, newb, zden, QM;
    1654             :     long idx;
    1655        3465 :     if (t == 1) zc = z;
    1656             :     else
    1657             :     {
    1658             :       GEN r;
    1659         161 :       do r = randS(sbase); while (degpol(QX_gcd(r, pol)));
    1660         161 :       zc = QXQ_mul(z, QXQ_sqr(r, pol), pol);
    1661             :     }
    1662        3465 :     q2 = Q_primpart(tracematrix(zc, b, pol));
    1663        3465 :     U = redquadric(b, q2, pol, QXQ_div(zc, polprime, pol));
    1664        4130 :     if (lg(U) < 4) continue;
    1665        3465 :     q2 = qf_apply_RgM(q2, U);
    1666        3465 :     newb = RgV_RgM_mul(b, U);
    1667             : 
    1668        3465 :     param = Q_primpart(qfparam(q2, qfsolve(q2), 1));
    1669        3465 :     param = RgM_to_RgXV_reverse(shallowtrans(param), 0);
    1670        3465 :     q1 = RgM_neg(tracematrix(QXQ_mul(zc, ttheta, pol), newb, pol));
    1671        3465 :     q1 = Q_remove_denom(qfeval(q1, param), &den);
    1672        3465 :     if (den) q1 = ZX_Z_mul(q1, den);
    1673        3465 :     if (!equali1(K)) q1 = RgX_Rg_mul(q1, K);
    1674        3465 :     QM = quartic_minim_all(q1, discF);
    1675        3465 :     q1 = gel(QM,1);
    1676        3465 :     zden = gmael(QM,2,1);
    1677        3465 :     Q = ZX_hyperellred(q1, &R);
    1678        3465 :     R = gmul(gmael(QM,2,2), R);
    1679        3465 :     if (pt_Q) *pt_Q = Q;
    1680        3465 :     Qk = shallowcopy(Q);
    1681        3465 :     (void) ZX_canon_neg(Qk);
    1682        3465 :     if (hash_haskey_long(&h, (void*)Qk, &idx)) continue;
    1683        3409 :     hash_insert_long(&h, (void*)Qk, 1); av2 = avma;
    1684        3409 :     if (DEBUGLEVEL>1) err_printf("  reduced quartic: Y^2 = %Ps\n", Q);
    1685             : 
    1686        3409 :     xz = projratpointxz(Q, lim, &zz);
    1687        3409 :     if (!xz)
    1688             :     {
    1689        1659 :       xz = projratpointxz2(Q, lim, &zz);
    1690        1659 :       if (!xz)
    1691             :       {
    1692        3409 :         if (pt_Q) return NULL; else continue;
    1693             :       }
    1694             :     }
    1695        1785 :     P = RgM_RgC_mul(R, xz); x = gel(P,1); y = gel(P,2);
    1696        1785 :     param = RgXV_homogenous_evaldeg(param, x, gpowers(y, 2));
    1697        1785 :     param = gmul(param, gdiv(den? mulii(K, den): K, gmul(zz, zden)));
    1698        1785 :     q0 = tracematrix(QXQ_mul(zc, tttheta, pol), newb, pol);
    1699        1785 :     x = gdiv(qfeval(q0, param), K);
    1700        1785 :     if (!issquareall(gdiv(poleval(pol, x), K), &y)) /* K y^2 = pol(x) */
    1701           0 :       pari_err_BUG("ellrank");
    1702        1785 :     P = mkvec2(x, y);
    1703        1785 :     if (DEBUGLEVEL) err_printf("Found point: %Ps\n", P);
    1704        1785 :     if (pt_Q) *pt_Q = gen_0;
    1705        1785 :     return gerepilecopy(av, P);
    1706             :   }
    1707         504 :   return NULL;
    1708             : }
    1709             : 
    1710             : static void
    1711        1603 : gtoset_inplace(GEN x)
    1712        1603 : { gen_sort_inplace(x, (void*)&cmp_universal, cmp_nodata, NULL); }
    1713             : 
    1714             : /* FIXME: export */
    1715             : static void
    1716        8267 : setlgall(GEN x, long L)
    1717             : {
    1718        8267 :   long i, l = lg(x);
    1719      179760 :   for(i = 1; i < l; i++) setlg(gel(x,i), L);
    1720        8267 : }
    1721             : 
    1722             : static long
    1723        8267 : dim_selmer(GEN p, GEN pol, GEN K, GEN vnf, GEN LS2, GEN helpLS2,
    1724             :            GEN *selmer, GEN *LS2chars, GEN *helpchars)
    1725             : {
    1726             :   pari_sp av;
    1727        8267 :   long dim, k, lvnf = lg(vnf);
    1728        8267 :   GEN X, L, LS2image, helpimage, pp = cgetg(lvnf, t_VEC);
    1729        8267 :   int pis2 = equaliu(p, 2);
    1730             : 
    1731       23968 :   for (k = 1; k < lvnf; k++)
    1732             :   {
    1733       15701 :     GEN v, nf = gel(vnf,k), PR = idealprimedec(nf, p);
    1734       15701 :     long j, l = lg(PR);
    1735       15701 :     gel(pp, k) = v = cgetg(l, t_VEC);
    1736       35119 :     for (j = 1; j < l; j++)
    1737             :     {
    1738       19418 :       GEN pr = gel(PR,j);
    1739       19418 :       gel(v,j) = pis2? log_prk_init(nf, pr, 1 + 2 * pr_get_e(pr), NULL)
    1740       19418 :                      : zkmodprinit(nf, pr);
    1741             :     }
    1742             :   }
    1743        8267 :   LS2image = veckpmodsquares(LS2, vnf, pp);
    1744        8267 :   *LS2chars = vconcat(*LS2chars, LS2image);
    1745        8267 :   helpimage = veckpmodsquares(helpLS2, vnf, pp);
    1746        8267 :   *helpchars = vconcat(*helpchars, helpimage);
    1747        8267 :   av = avma;
    1748        8267 :   L = elllocalimage(pol, K, vnf, p, pp, helpimage);
    1749        8267 :   X = Flm_ker(shallowconcat(LS2image, L), 2); setlgall(X, lg(LS2image));
    1750             :   /* intersect(LS2image, locim) = LS2image.X */
    1751        8267 :   *selmer = Flm_intersect_i(*selmer, shallowconcat(Flm_ker(LS2image,2), X), 2);
    1752        8267 :   *selmer = gerepileupto(av, Flm_image(*selmer, 2));
    1753        8267 :   dim = lg(*selmer)-1; return (dim == Flm_rank(helpimage,2))? dim: -1;
    1754             : }
    1755             : static long
    1756         476 : get_row(GEN vnf, GEN K)
    1757             : {
    1758         476 :   long k, sK = signe(K), n = lg(vnf)-1;
    1759             :   GEN R;
    1760         476 :   if (n == 1) return sK > 0? 1: 3;
    1761         280 :   if (n == 2)
    1762             :   {
    1763          98 :     GEN a = gel(nf_get_roots(gel(vnf,1)), 1);
    1764          98 :     GEN b = gel(nf_get_roots(gel(vnf,2)), sK > 0? 1: 2);
    1765          98 :     if (cmprr(a, b) < 0) return sK > 0? 1: 3;
    1766          84 :     return sK > 0? 2: 1;
    1767             :   }
    1768         182 :   R = cgetg(4, t_VEC);
    1769         728 :   for (k = 1; k <= 3; k++) gel(R, k) = gel(nf_get_roots(gel(vnf,k)), 1);
    1770         182 :   return sK > 0? vecindexmin(R): vecindexmax(R);
    1771             : }
    1772             : 
    1773             : static GEN
    1774         854 : ell2selmer(GEN ell, GEN ell_K, GEN help, GEN K, GEN vbnf,
    1775             :            long effort, long flag, long prec)
    1776             : {
    1777             :   GEN KP, pol, vnf, vpol, vcrt, sbase, LS2, factLS2, sqrtLS2, signs;
    1778             :   GEN selmer, helpLS2, LS2chars, helpchars, newselmer, factdisc, badprimes;
    1779             :   GEN helplist, listpoints, etors2, p, covers, disc, discF;
    1780         854 :   long i, k, n, tors2, mwrank, dim, nbpoints, lfactdisc, t, u, sha2 = 0;
    1781             :   forprime_t T;
    1782             : 
    1783         854 :   pol = ell2pol(ell);
    1784         854 :   help = ellsearchtrivialpoints(ell_K, flag ? NULL:muluu(LIMTRIV,effort+1), help);
    1785         854 :   help = elltwistpoints(help, ginv(K)); /* points on K Y^2 = pol(X) */
    1786         854 :   n = lg(vbnf) - 1; tors2 = n - 1;
    1787         854 :   etors2 = vecslice(help,1, tors2);
    1788         854 :   gtoset_inplace(etors2);
    1789         854 :   KP = gel(absZ_factor(K), 1);
    1790         854 :   disc = ZX_disc(pol);
    1791         854 :   factdisc = mkvec3(KP, mkcol(gen_2), gel(absZ_factor(disc), 1));
    1792         854 :   factdisc = ZV_sort_uniq(shallowconcat1(factdisc));
    1793         854 :   discF = mkvec2(gmul(K,disc), factdisc);
    1794         854 :   badprimes = cgetg(n+1, t_VEC);
    1795         854 :   vnf = cgetg(n+1, t_VEC);
    1796         854 :   vpol = cgetg(n+1, t_VEC);
    1797         854 :   vcrt = cgetg(n+1, t_VEC);
    1798         854 :   LS2 = cgetg(n+1, t_VEC);
    1799         854 :   factLS2 = cgetg(n+1, t_VEC);
    1800         854 :   sqrtLS2 = cgetg(n+1, t_VEC);
    1801        2310 :   for (k = 1; k <= n; k++)
    1802             :   {
    1803        1456 :     GEN T, Tinv, id, f, sel, L, S, nf, NF = gel(vbnf,k);
    1804             :     long i, lk;
    1805        1456 :     nf = (n == 1)? bnf_get_nf(NF): NF;
    1806        1456 :     gel(vnf, k) = nf;
    1807        1456 :     gel(vpol, k) = T = nf_get_pol(nf);
    1808        1456 :     Tinv = RgX_div(pol, gel(vpol, k));
    1809        1456 :     gel(vcrt, k) = QX_mul(QXQ_inv(T, Tinv), T); /* 0 mod T, 1 mod pol/T */
    1810             : 
    1811        1456 :     id = idealadd(nf, nf_get_index(nf), ZX_deriv(T));
    1812        1456 :     f = nf_pV_to_prV(nf, KP); settyp(f, t_COL);
    1813        1456 :     f = mkvec3(gel(idealfactor(nf, Tinv), 1),
    1814        1456 :                gel(idealfactor(nf, id), 1), f);
    1815        1456 :     gel(badprimes, k) = S = gtoset(shallowconcat1(f));
    1816        1456 :     if (n == 1)
    1817             :     {
    1818         434 :       sel = bnfselmer(NF, S);
    1819         434 :       obj_free(NF); /* units */
    1820             :     }
    1821        1022 :     else if (degpol(T) == 1)
    1822         784 :       sel = bnfselmerQ(S);
    1823             :     else /* degree 2 */
    1824         238 :       sel = nf2selmer_quad(NF, S);
    1825        1456 :     gel(LS2, k) = L = gel(sel, 1); lk = lg(L);
    1826        1456 :     gel(factLS2, k) = gel(sel, 2);
    1827        1456 :     gel(sqrtLS2, k) = gel(sel, 3);
    1828       14252 :     for (i = 1; i < lk; i++)
    1829             :     {
    1830       12796 :       GEN z = gel(L,i); /* z mod T, 1 mod (pol/T) */
    1831       12796 :       gel(L,i) = grem(gadd(z, gmul(gsubsg(1,z), gel(vcrt,k))), pol);
    1832             :     }
    1833             :   }
    1834         854 :   sbase = shallowconcat1(vecsmallbasis(vnf, vcrt, pol));
    1835         854 :   if (DEBUGLEVEL>2) err_printf("   local badprimes = %Ps\n", badprimes);
    1836         854 :   LS2 = shallowconcat1(LS2);
    1837         854 :   helpLS2 = ellLS2image(pol, help, K, vpol, vcrt);
    1838         854 :   LS2chars = matselmersign(vnf, vpol, LS2);
    1839         854 :   helpchars = matselmersign(vnf, vpol, helpLS2);
    1840         854 :   signs = NULL;
    1841         854 :   if (signe(ell_get_disc(ell)) > 0) signs = Flm_row(LS2chars, get_row(vnf,K));
    1842         854 :   selmer = kernorm(LS2, factdisc, pol, signs);
    1843         854 :   forprime_init(&T, gen_2, NULL); lfactdisc = lg(factdisc); dim = -1;
    1844        7070 :   for (i = 1; dim < 0 && i < lfactdisc; i++)
    1845        6216 :     dim = dim_selmer(gel(factdisc,i), pol, K, vnf, LS2, helpLS2,
    1846             :                      &selmer,&LS2chars,&helpchars);
    1847        2905 :   while (dim < 0 && Flm_rank(Flm_mul(LS2chars, selmer, 2), 2) < lg(selmer)-1)
    1848             :   {
    1849        2513 :     while ((p = forprime_next(&T)) && ZV_search(factdisc, p));
    1850        2051 :     dim = dim_selmer(p, pol, K, vnf, LS2, helpLS2,
    1851             :                      &selmer,&LS2chars,&helpchars);
    1852             :   }
    1853         854 :   helplist = gel(Flm_indexrank(helpchars,2), 2);
    1854         854 :   help = shallowextract(help, helplist);
    1855         854 :   helpchars = shallowextract(helpchars, helplist);
    1856         854 :   helpLS2 = shallowextract(helpLS2, helplist);
    1857         854 :   dim = lg(selmer)-1;
    1858         854 :   if (DEBUGLEVEL) err_printf("Selmer rank: %ld\n", dim);
    1859         854 :   newselmer = Flm_invimage(Flm_mul(LS2chars, selmer, 2), helpchars, 2);
    1860         854 :   nbpoints = lg(help) - 1;
    1861         854 :   if (flag==1)
    1862             :   {
    1863          49 :     GEN u = nbpoints? Flm_mul(selmer,Flm_suppl(newselmer,2), 2): selmer;
    1864          49 :     long l = lg(u);
    1865          49 :     GEN z = cgetg(l, t_VEC);
    1866         154 :     for (i = 1; i < l; i++) gel(z,i) = RgXQV_factorback(LS2, gel(u,i), pol);
    1867          49 :     return mkvec2(mkvec3(vnf,sbase,pol), z);
    1868             :   }
    1869         805 :   else if (flag==2)
    1870             :   {
    1871          56 :     GEN u = nbpoints ? Flm_mul(selmer,Flm_suppl(newselmer,2), 2): selmer;
    1872          56 :     long l = lg(u);
    1873          56 :     GEN P = cgetg(l, t_VEC), b;
    1874         147 :     for (i = 1; i < l; i++)
    1875             :     {
    1876          91 :       b = liftselmerinit(gel(u,i), vnf, sqrtLS2, factLS2, badprimes, vcrt, pol);
    1877          91 :       gel(P,i) = liftselmer_cover(b, gel(u,i), LS2, pol, discF, K);
    1878             :     }
    1879          56 :     return P;
    1880             :   }
    1881         749 :   listpoints = vec_lengthen(help, dim); /* points on ell */
    1882         749 :   covers = zerovec(dim);
    1883        5628 :   for (i=1; i <= dim; i++)
    1884             :   {
    1885        4879 :     GEN b, P, expo, vec = vecsmall_ei(dim, i);
    1886        4879 :     if (Flm_Flc_invimage(newselmer, vec, 2)) continue;
    1887        2464 :     expo = Flm_Flc_mul(selmer, vec, 2);
    1888        2464 :     b = liftselmerinit(expo, vnf, sqrtLS2, factLS2, badprimes, vcrt, pol);
    1889        2464 :     P = liftselmer(b, expo, sbase, LS2, pol, discF, K, 1, &gel(covers,i));
    1890        2464 :     if (P)
    1891             :     {
    1892        1449 :       gel(listpoints, ++nbpoints) = P; /* new point on ell */
    1893        1449 :       gel(newselmer, nbpoints) = vec;
    1894        1449 :       setlg(newselmer, nbpoints+1);
    1895             :     }
    1896             :   }
    1897         749 :   if (nbpoints < dim)
    1898             :   {
    1899             :     long i, j;
    1900         294 :     GEN M = cgetg(dim+1, t_MAT), selker;
    1901         294 :     GEN D = mulii(muliu(absi(disc), 27*4096), powiu(K,6));
    1902         294 :     GEN FD = ZV_sort_uniq(shallowconcat1(mkvec2(mkcol3s(3,5,7), factdisc)));
    1903             : 
    1904        1862 :     for (i = 1; i <= dim; i++) gel(M,i) = cgetg(dim+1, t_COL);
    1905        1862 :     for (i = 1; i <= dim; i++)
    1906        9051 :       for (j = 1; j <= i; j++)
    1907             :       {
    1908             :         GEN Q;
    1909        7483 :         if (isintzero(gel(covers,i)))
    1910        3115 :           Q = gen_0;
    1911        4368 :         else if (i==j)
    1912         973 :           Q = quartic_findunit(D, gel(covers,i));
    1913             :         else
    1914             :         {
    1915        3395 :           GEN e = Flv_add(gel(selmer,i), gel(selmer,j), 2);
    1916        3395 :           GEN b = liftselmerinit(e, vnf, sqrtLS2, factLS2, badprimes, vcrt, pol);
    1917        3395 :           Q = quartic_findunit(D, gel(liftselmer_cover(b, e, LS2, pol, discF, K),1));
    1918             :         }
    1919        7483 :         gmael(M,j,i) = gmael(M,i,j) = Q;
    1920             :       }
    1921         294 :     selker = F2m_to_Flm(F2m_ker(matcassels(FD, M)));
    1922         294 :     sha2 = dim - (lg(selker)-1);
    1923         294 :     dim = lg(selker)-1;
    1924        1134 :     for (t=1, u=1; nbpoints < dim && effort > 0; t++)
    1925             :     {
    1926         840 :       pari_sp btop = avma;
    1927             :       GEN expo, b, P, vec;
    1928        1120 :       do vec = Flm_Flc_mul(selker,random_Flv(dim, 2), 2);
    1929        1120 :       while (zv_equal0(vec) || Flm_Flc_invimage(newselmer, vec, 2));
    1930         840 :       expo = Flm_Flc_mul(selmer, vec, 2);
    1931         840 :       b = liftselmerinit(expo, vnf, sqrtLS2, factLS2, badprimes, vcrt, pol);
    1932         840 :       P = liftselmer(b, expo, sbase, LS2, pol, discF, K, u, NULL);
    1933         840 :       if (P)
    1934             :       {
    1935         336 :         gel(listpoints, ++nbpoints) = P;
    1936         336 :         gel(newselmer, nbpoints) = vec;
    1937         336 :         setlg(newselmer, nbpoints+1);
    1938         504 :       } else set_avma(btop);
    1939         840 :       if (t == dim) { t = 0; u++; effort--; }
    1940             :     }
    1941             :   }
    1942         749 :   setlg(listpoints, nbpoints+1);
    1943         749 :   mwrank = nbpoints - tors2;
    1944         749 :   if (odd(dim - nbpoints)) mwrank++;
    1945         749 :   gtoset_inplace(listpoints);
    1946         749 :   listpoints = setminus(listpoints, etors2);
    1947         749 :   listpoints = elltwistpoints(listpoints, K);
    1948         749 :   listpoints = vecellabs(ellQ_genreduce(ell_K, listpoints, NULL, prec));
    1949         749 :   return mkvec4(utoi(mwrank), utoi(dim-tors2), utoi(sha2), listpoints);
    1950             : }
    1951             : 
    1952             : GEN
    1953          49 : ell2selmer_basis(GEN ell, GEN *cb, long prec)
    1954             : {
    1955          49 :   GEN E = ellminimalbmodel(ell, cb);
    1956          49 :   GEN S = ell2selmer(E, E, NULL, gen_1, makevbnf(E, prec), 0, 1, prec);
    1957          49 :   obj_free(E); return S;
    1958             : }
    1959             : 
    1960             : static void
    1961          98 : ellrank_get_nudur(GEN E, GEN F, GEN *nu, GEN *du, GEN *r)
    1962             : {
    1963          98 :   GEN ea2 = ell_get_a2(E), ea2t = ell_get_a2(F);
    1964          98 :   GEN ec4 = ell_get_c4(E), ec4t = ell_get_c4(F);
    1965          98 :   GEN ec6 = ell_get_c6(E), ec6t = ell_get_c6(F);
    1966             :   GEN N, D, d;
    1967          98 :   if (signe(ec4)==0)
    1968             :   {
    1969          21 :     if (!Z_ispowerall(mulii(ec6,sqri(ec6t)),3,&N))
    1970           7 :       pari_err_TYPE("ellrank",F);
    1971          14 :     D = ec6t;
    1972             :   }
    1973          77 :   else if (signe(ec6)==0)
    1974             :   {
    1975          21 :     if (!Z_issquareall(mulii(ec4,ec4t),&N))
    1976           7 :       pari_err_TYPE("ellrank",F);
    1977          14 :     D = ec4t;
    1978             :   }
    1979             :   else
    1980             :   {
    1981          56 :     GEN d46 = mulii(gcdii(ec4,ec4t),gcdii(ec6,ec6t));
    1982          56 :     N = diviiexact(mulii(ec6,ec4t),d46);
    1983          56 :     D = diviiexact(mulii(ec6t,ec4),d46);
    1984             :   }
    1985          84 :   d = gcdii(N, D);
    1986          84 :   *nu = diviiexact(N, d); *du = diviiexact(D, d);
    1987          84 :   *r  = diviuexact(subii(mulii(*nu,ea2t),mulii(*du,ea2)),3);
    1988          84 : }
    1989             : 
    1990             : static GEN
    1991         847 : ellrank_flag(GEN e, long effort, GEN help, long flag, long prec)
    1992             : {
    1993         847 :   pari_sp ltop = avma;
    1994             :   GEN urst, v, vbnf, eK;
    1995         847 :   GEN et = NULL, K = gen_1, nu = NULL, du = NULL, r = NULL, urstK = NULL;
    1996         847 :   long newell = 0;
    1997             : 
    1998         847 :   if (lg(e)==3 && typ(e)==t_VEC) { et = gel(e,2); e = gel(e,1); }
    1999         847 :   if (lg(e)==4 && typ(e)==t_VEC)
    2000             :   {
    2001         126 :     vbnf = gel(e,3); urst = gel(e,2);
    2002         126 :     e = gel(e,1); checkell_Q(e);
    2003         126 :     if (!ell_is_integral(e)) pari_err_TYPE("ellrank [nonintegral model]",e);
    2004         119 :     if (signe(ell_get_a1(e))) pari_err_TYPE("ellrank [a1 != 0]", e);
    2005         112 :     if (signe(ell_get_a3(e))) pari_err_TYPE("ell2rank [a3 != 0]", e);
    2006             :   }
    2007             :   else
    2008             :   {
    2009         721 :     checkell_Q(e);
    2010         721 :     e = ellminimalbmodel(e, &urst);
    2011         721 :     newell = 1;
    2012         721 :     vbnf = makevbnf(e, prec);
    2013             :   }
    2014         826 :   if (et)
    2015             :   {
    2016         105 :     checkell_Q(et);
    2017         105 :     if (!gequal(ell_get_j(et),ell_get_j(e))) pari_err_TYPE("ellrank",et);
    2018          98 :     et = ellminimalbmodel(et, &urst);
    2019          98 :     ellrank_get_nudur(e, et, &nu, &du, &r);
    2020          84 :     K = mulii(nu, du);
    2021          84 :     urstK = mkvec4(nu, mulii(nu,r), gen_0, gen_0);
    2022             :   }
    2023         805 :   if (help)
    2024             :   {
    2025          42 :     if (urst) help = ellchangepoint(help, urst);
    2026          42 :     if (et) help = ellchangepointinv(help, urstK);
    2027             :   }
    2028         805 :   eK = elltwistequation(e, K);
    2029             :   /* help is a vector of rational points [x,y] satisfying K y^2 = pol(x) */
    2030             :   /* [Kx, K^2y] is on eK: Y^2 = K^3 pol(X/K)  */
    2031         805 :   if (help) check_oncurve(eK, help);
    2032         805 :   v = ell2selmer(e, eK, help, K, vbnf, effort, flag, prec);
    2033         805 :   if (flag==0)
    2034             :   {
    2035         749 :     if (et)   gel(v,4) = ellchangepoint(gel(v,4), urstK);
    2036         749 :     if (urst) gel(v,4) = ellchangepointinv(gel(v,4), urst);
    2037             :   }
    2038             :   else
    2039             :   {
    2040          56 :     long i, l = lg(v);
    2041         147 :     for (i = 1; i < l; i++)
    2042             :     {
    2043          91 :       if (et)   gmael(v,i,2) = ellchangepoint(gmael(v,i,2), urstK);
    2044          91 :       if (urst) gmael(v,i,2) = ellchangepointinv(gmael(v,i,2), urst);
    2045             :     }
    2046             :   }
    2047         805 :   if (newell) obj_free(e);
    2048         805 :   if (eK != e) obj_free(eK);
    2049         805 :   return gerepilecopy(ltop, v);
    2050             : }
    2051             : 
    2052             : GEN
    2053         791 : ellrank(GEN e, long effort, GEN help, long prec)
    2054             : {
    2055         791 :   return ellrank_flag(e, effort, help, 0, prec);
    2056             : }
    2057             : 
    2058             : GEN
    2059          56 : ell2cover(GEN ell, long prec)
    2060             : {
    2061          56 :   return ellrank_flag(ell, 0, NULL, 2, prec);
    2062             : }

Generated by: LCOV version 1.13