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 = gc_upto(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 165718 : to_principal_unit(GEN nf, GEN x, GEN pr, GEN sprk)
158 : {
159 165718 : if (pr_get_f(pr) != 1)
160 : {
161 18200 : GEN prk = gel(sprk,3);
162 18200 : x = nfpowmodideal(nf, x, gmael(sprk,5,1), prk);
163 : }
164 165718 : 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 119336 : FpX_issquare(GEN q, GEN p)
453 : {
454 119336 : GEN F = FpX_factor_squarefree(q,p);
455 119336 : long i, l = lg(F);
456 161658 : for (i = 1; i < l; i+=2)
457 120582 : if (degpol(gel(F,i)) > 0) return 0;
458 41076 : return 1;
459 : }
460 :
461 : static GEN
462 118580 : hyperell_red(GEN q, GEN p)
463 : {
464 : GEN Q;
465 118580 : long v = ZX_pvalrem(q, p, &Q);
466 118580 : if (v == 1) return q;
467 118580 : return odd(v)? ZX_Z_mul(Q, p): Q;
468 : }
469 :
470 : static GEN
471 121659 : hyperell_reg_point(GEN q, GEN p)
472 : {
473 : GEN Q, F;
474 121659 : long i, l, v = ZX_pvalrem(q, p, &Q);
475 121659 : if (v != 1) q = odd(v)? ZX_Z_mul(Q, p): Q;
476 121659 : if (!odd(v))
477 : {
478 119336 : GEN qr = FpX_red(q, p);
479 119336 : if (!FpX_issquare(qr,p) || Fp_issquare(leading_coeff(qr), p))
480 118580 : return mkvec2s(0,1);
481 : }
482 3079 : F = FpX_roots(Q, p); l = lg(F);
483 3086 : for (i = 1; i < l; i++)
484 : {
485 3079 : GEN r = gel(F,i), s = hyperell_reg_point(ZX_affine(q, p, r), p);
486 3079 : if (s)
487 3072 : retmkvec2(addii(r,mulii(gel(s,1),p)), mulii(gel(s,2),p));
488 : }
489 7 : return NULL;
490 : }
491 :
492 : static GEN
493 1432 : hyperell_lift(GEN q, GEN x, GEN p)
494 : {
495 : long e;
496 1432 : GEN qy2 = ZX_Z_sub(q, sqri(p));
497 1432 : for (e = 2; ; e<<=1)
498 702 : {
499 2134 : pari_sp av = avma;
500 2134 : GEN z = ZpX_liftroot(qy2, x, p, e);
501 2134 : if (signe(z) == 0) z = powiu(p, e);
502 2134 : if (Zp_issquare(poleval(q, z), p)) return z;
503 702 : set_avma(av);
504 : }
505 : }
506 :
507 : static GEN
508 59290 : affine_apply(GEN r, GEN x)
509 : {
510 59290 : return addii(mulii(gel(r,2),x), gel(r,1));
511 : }
512 :
513 : static GEN
514 59290 : Qp_hyperell_solve_odd(GEN q, GEN p)
515 : {
516 59290 : GEN qi = RgX_recip_shallow(q);
517 59290 : GEN r = hyperell_reg_point(q, p), qr = NULL, qrp = NULL;
518 59290 : GEN s = hyperell_reg_point(qi, p), qs = NULL, qsp = NULL;
519 59290 : if (!r && !s) return NULL;
520 59290 : if (r)
521 : {
522 59290 : qr = hyperell_red(ZX_affine(q, gel(r,2), gel(r,1)), p);
523 59290 : qrp = FpX_deriv(qr, p);
524 : }
525 59290 : if (s)
526 : {
527 59290 : qs = hyperell_red(ZX_affine(qi, gel(s,2), gel(s,1)), p);
528 59290 : qsp = FpX_deriv(qs, p);
529 : }
530 : while(1)
531 14882 : {
532 74172 : GEN x = randomi(p);
533 74172 : if (r)
534 : {
535 74172 : GEN y2 = FpX_eval(qr, x, p);
536 74172 : if (Fp_issquare(y2,p))
537 : {
538 49015 : if (signe(y2))
539 44589 : return affine_apply(r,x);
540 4426 : if (signe(FpX_eval(qrp, x, p)))
541 : {
542 1105 : x = hyperell_lift(qr, x, p);
543 1105 : return affine_apply(r,x);
544 : }
545 : }
546 : }
547 28478 : if (s)
548 : {
549 28478 : GEN y2 = FpX_eval(qs, x, p);
550 28478 : if (Fp_issquare(y2,p))
551 : {
552 15617 : if (signe(x)==0) x = p;
553 15617 : if (signe(y2))
554 13269 : return ginv(affine_apply(s,x));
555 2348 : if (signe(FpX_eval(qsp, x, p)))
556 : {
557 327 : GEN xl = hyperell_lift(qs, x, p);
558 327 : return ginv(affine_apply(s,xl));
559 : }
560 : }
561 : }
562 : }
563 : }
564 :
565 : static GEN
566 490 : Q2_hyperell_lift(GEN p, GEN q, long x, long y)
567 : {
568 : GEN T, h;
569 : long e;
570 490 : if (y==0) y = 2;
571 490 : T = ZX_sub(p, ZX_Z_add(ZX_mulu(q, y), sqru(y)));
572 490 : h = ZX_add(ZX_sqr(q), ZX_shifti(p, 2));
573 490 : for (e = 3; ; e++)
574 952 : {
575 1442 : pari_sp av = avma;
576 1442 : GEN r = ZpX_liftroot(T, utoi(x), gen_2, e);
577 1442 : if (signe(r) == 0) r = int2n(e);
578 1442 : if (Zp_issquare(poleval(h, r), gen_2)) return r;
579 952 : set_avma(av);
580 : }
581 : return NULL;/*LCOV_EXCL_LINE*/
582 : }
583 :
584 : static GEN
585 2947 : Q2_hyperell_regpoint(GEN P, GEN Q)
586 : {
587 : long x;
588 2947 : GEN p = ZX_to_F2x(P), dp = F2x_deriv(p);
589 2947 : GEN q = ZX_to_F2x(Q), dq = F2x_deriv(q);
590 :
591 4907 : for (x = 0; x <= 1; x++)
592 : {
593 4326 : long px = F2x_eval(p,x), qx = F2x_eval(q,x);
594 : long dpx, dqx;
595 4326 : if (qx == 1)
596 : {
597 2100 : if(px == 0) return x==0 ? gen_2: gen_1;
598 224 : continue;
599 : }
600 2226 : dpx = F2x_eval(dp,x);
601 2226 : dqx = F2x_eval(dq,x);
602 2226 : if (odd(dqx * px + dpx))
603 490 : return Q2_hyperell_lift(P, Q, x, px);
604 : }
605 581 : return NULL;
606 : }
607 :
608 : static GEN
609 2947 : Q2_hyperell_solve_affine(GEN p, GEN q)
610 : {
611 2947 : pari_sp av = avma;
612 : GEN R, p4, q4;
613 : long x;
614 : while(1)
615 2394 : {
616 : GEN pp, p0, p1;
617 5341 : long vp = ZX_lval(p, 2);
618 5341 : long vq = ZX_lval(q, 2);
619 5341 : long w = minss(vp>>1, vq);
620 5341 : p = ZX_shifti(p, -2*w);
621 5341 : q = ZX_shifti(q, -w);
622 5341 : if (ZX_lval(q,2)==0) break;
623 2667 : pp = RgX_splitting(p,2); p0 = gel(pp,1); p1 = gel(pp,2);
624 2667 : if (ZX_lval(p1,2)==0 || ZX_lval(p0,2)>=1) break;
625 2394 : p = ZX_sub(p, ZX_mul(p0, ZX_add(q, p0)));
626 2394 : q = ZX_add(q, ZX_shifti(p0, 1));
627 : }
628 2947 : R = Q2_hyperell_regpoint(p, q);
629 2947 : if (R) return gc_INT(av, R);
630 581 : p4 = ZX_to_Flx(p,4);
631 581 : q4 = ZX_to_Flx(q,4);
632 763 : for (x = 0; x <= 1; x++)
633 : {
634 742 : ulong px = Flx_eval(p4, x, 4);
635 742 : ulong qx = Flx_eval(q4, x, 4);
636 742 : if (px == 0 || (1+qx+3*px)%4==0)
637 : {
638 560 : GEN p2 = ZX_affine(p, gen_2, utoi(x));
639 560 : GEN q2 = ZX_affine(q, gen_2, utoi(x));
640 560 : GEN S = Q2_hyperell_solve_affine(p2, q2);
641 560 : if (S) return gc_INT(av, addiu(shifti(S,1),x));
642 : }
643 : }
644 21 : return gc_NULL(av);
645 : }
646 :
647 : static GEN
648 2366 : Q2_hyperell_solve(GEN P)
649 : {
650 2366 : long v = varn(P);
651 2366 : GEN S = Q2_hyperell_solve_affine(P, pol_0(v));
652 2366 : if (!S) S = ginv(Q2_hyperell_solve_affine(RgX_recip(P), pol_0(v)));
653 2366 : return S;
654 : }
655 :
656 : static GEN
657 61656 : hyperell_local_solve(GEN q, GEN p)
658 : {
659 61656 : if (equaliu(p,2))
660 2366 : return Q2_hyperell_solve(q);
661 59290 : return Qp_hyperell_solve_odd(q, p);
662 : }
663 :
664 : /*******************************************************************/
665 : /* BINARY QUARTIC */
666 : /*******************************************************************/
667 : static int
668 107604 : Qp_issquare(GEN a, GEN p)
669 : {
670 107604 : GEN b = typ(a) == t_INT? a: mulii(gel(a,1), gel(a,2));
671 107604 : return Zp_issquare(b, p);
672 : }
673 :
674 : static GEN
675 18780 : quartic_IJ(GEN g)
676 : {
677 18780 : GEN a = gel(g, 6), b = gel(g, 5), c = gel(g, 4), d = gel(g, 3), e = gel(g, 2);
678 18780 : GEN ae = gmul(a,e), bd = gmul(b,d), c2 = gsqr(c);
679 : /* 12ae - 3bd + c^2 */
680 18780 : GEN I = gadd(gsub(gmulsg(12, ae), gmulsg(3, bd)), c2);
681 : /* c(72ae + 9bd - 2c^2) - 27ad^2 - 27eb^2*/
682 18780 : 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 18780 : return mkvec2(I, J);
685 : }
686 :
687 : static GEN
688 2366 : quartic_hessiandd(GEN g)
689 : {
690 2366 : GEN a = gel(g, 6), b = gel(g, 5), c = gel(g, 4), d = gel(g, 3), e = gel(g, 2);
691 2366 : GEN a8 = gmul2n(a, 3), p4 = gsub(gmulsg(3, gsqr(b)), gmul(a8, c));
692 2366 : GEN p3 = gsub(gmul(b, c), gmul(gmulsg(6, a), d));
693 2366 : GEN p2 = gsub(gmulsg(8, gsqr(c)), gmulsg(12, gadd(gmul(b, d), gmul(a8, e))));
694 2366 : return deg2pol_shallow(gmulgu(p4,12), gmulgu(p3,24), p2, 1);
695 : }
696 :
697 : static GEN
698 13090 : quartic_cubic(GEN g, long v)
699 : {
700 13090 : GEN a = gel(g, 6), b = gel(g, 5), c = gel(g, 4);
701 13090 : GEN a3 = gdivgu(a,3);
702 13090 : return deg1pol(gmul2n(a3,2), gsub(gsqr(b),gmul2n(gmul(a3,c),3)), v);
703 : }
704 :
705 : static GEN
706 6972 : quarticinv_pol(GEN IJ)
707 : {
708 6972 : GEN I = gel(IJ,1), J = gel(IJ,2);
709 6972 : return mkpoln(4, gen_1, gen_0, gmulgs(I,-3), J);
710 : }
711 : static GEN
712 2366 : quartic_H(GEN g, GEN *pT)
713 : {
714 2366 : GEN a = gel(g, 6), b = gel(g, 5), c = gel(g, 4);
715 2366 : GEN IJ = quartic_IJ(g), I = gel(IJ, 1);
716 2366 : GEN ddh = quartic_hessiandd(g);
717 2366 : GEN ddg = deg2pol_shallow(gmulgu(a,12), gmulgu(b,6), gmulgu(c,2), 1);
718 2366 : *pT = quarticinv_pol(IJ);
719 2366 : return deg2pol_shallow(stoi(-8), gmul2n(ddg,2), gadd(ddh,gmul2n(I,3)),0);
720 : }
721 :
722 : static GEN
723 4606 : quartic_disc(GEN q)
724 : {
725 4606 : GEN IJ = quartic_IJ(q), I = gel(IJ,1), J = gel(IJ,2);
726 4606 : return gsub(gmul2n(gpowgs(I, 3), 2), gsqr(J));
727 : }
728 :
729 : static GEN
730 7202 : quartic_minim_all(GEN F, GEN discF)
731 : {
732 7202 : GEN IJ = quartic_IJ(F), I = gel(IJ,1), J = gel(IJ,2);
733 7202 : GEN g = Z_ppo(gcdii(I,J), gel(discF,1));
734 7202 : GEN plist = ZV_sort_uniq_shallow(shallowconcat(gel(absZ_factor(g),1),gel(discF,2)));
735 7202 : GEN W, C, PQ = hyperellminimalmodel(F, &C, plist);
736 7202 : GEN P = gel(PQ,1), Q = gel(PQ,2);
737 7202 : if (signe(Q)==0)
738 798 : W = mkvec2(P, C);
739 : else
740 6404 : W = mkvec2(ZX_add(ZX_shifti(P,2),ZX_sqr(Q)),mkvec2(shifti(gel(C,1),-1),gel(C,2)));
741 7202 : return W;
742 : }
743 :
744 : /*******************************************************************/
745 : /* Cassels' pairing */
746 : /*******************************************************************/
747 :
748 : static GEN
749 6230 : nfsqrt(GEN nf, GEN z)
750 : {
751 6230 : long v = fetch_var_higher();
752 6230 : GEN R = nfroots(nf, deg2pol_shallow(gen_m1, gen_0, z, v));
753 6230 : delete_var();
754 6230 : return lg(R)==1 ? NULL: gel(R,1);
755 : }
756 :
757 : static GEN
758 6230 : nfsqrt_safe(GEN nf, GEN z)
759 : {
760 6230 : GEN r = nfsqrt(nf, z);
761 6230 : if (!r) pari_err_BUG("ellrank");
762 6230 : return r;
763 : }
764 :
765 : static GEN
766 2366 : vecnfsqrtmod(GEN x, GEN P)
767 8596 : { pari_APPLY_same(gmodulo(nfsqrt_safe(gel(x,i), P), gel(x,i))) }
768 :
769 : static GEN
770 2366 : enfsqrt(GEN T, GEN P)
771 : {
772 2366 : GEN F = gel(ZX_factor(T),1);
773 2366 : 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 602 : cassels_oo_solve_i(GEN q, GEN g)
780 : {
781 602 : long dg = degpol(g);
782 : GEN D, a, b, c;
783 :
784 602 : if (dg == 0 || signe(leading_coeff(q)) > 0) return 1;
785 329 : if (signe(gel(q,2)) > 0) return signe(gel(g,2));
786 301 : 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 301 : if (dg == 1) return ZX_sturmpart(q, mkvec2(mkmoo(), gdiv(negi(c), b)))? -1: 1;
790 287 : a = gel(g,4); D = subii(sqri(b), shifti(mulii(a,c), 2)); /* g = ax^2+bx+c */
791 287 : 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 266 : 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 252 : return (signe(gel(q,2)) > 0 ||
798 518 : ZX_sturmpart(ZX_graeffe(q), mkvec2(gen_0, D)))? -1: 1;
799 : }
800 : static int
801 602 : cassels_oo_solve(GEN q, GEN g)
802 602 : { pari_sp av = avma; return gc_int(av, cassels_oo_solve_i(q, g)); }
803 :
804 : static GEN
805 61656 : cassels_Qp_solve(GEN q, GEN gam, GEN p)
806 : {
807 61656 : pari_sp av = avma;
808 61656 : GEN a = hyperell_local_solve(q, p);
809 61656 : GEN c = poleval(gam,a);
810 : long e;
811 61656 : if (!gequal0(c)) return c;
812 42 : for (e = 2; ; e++)
813 0 : {
814 42 : GEN b = gadd(a, powiu(p,e));
815 42 : if (Qp_issquare(poleval(q, b), p))
816 : {
817 42 : c = poleval(gam, b);
818 42 : if (!gequal0(c)) return gc_upto(av,c);
819 : }
820 : }
821 : }
822 :
823 : static GEN
824 4732 : to_ZX(GEN a, long v) { return typ(a)==t_INT? scalarpol_shallow(a,v): a; }
825 :
826 : static GEN
827 4606 : quartic_findunit(GEN D, GEN q)
828 : {
829 4606 : GEN T = quarticinv_pol(quartic_IJ(q));
830 : while(1)
831 1386 : {
832 5992 : pari_sp av = avma;
833 5992 : GEN z = quartic_cubic(q,0);
834 5992 : if (signe(QXQ_norm(z,T)))
835 4606 : return absequalii(quartic_disc(q), D)? q: ZX_shifti(q, 2);
836 1386 : set_avma(av);
837 1386 : 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 2366 : casselspairing(GEN FD, GEN q1, GEN q2, GEN q3)
848 : {
849 2366 : pari_sp av = avma;
850 2366 : GEN T, H = quartic_H(q1, &T);
851 2366 : GEN z1 = quartic_cubic(q1, 0);
852 2366 : GEN z2 = quartic_cubic(q2, 0);
853 2366 : GEN z3 = quartic_cubic(q3, 0);
854 2366 : GEN m = to_ZX(enfsqrt(T, QXQ_mul(QX_mul(z1,z2),z3,T)), 0);
855 2366 : GEN Hm = RgXQ_mul(QXQ_div(m, z1, T), H, T); /* deg(Hm) >= 2 */
856 2366 : GEN gam = to_ZX(Q_primpart(gel(Hm,4)),1);
857 2366 : GEN a = leading_coeff(q2), Fa = gel(absZ_factor(a),1);
858 2366 : GEN F = ZV_sort_uniq_shallow(shallowconcat1(mkvec2(Fa, FD)));
859 2366 : long i, e = 0, lF = lg(F);
860 2366 : if (signe(a) <= 0)
861 : {
862 602 : if (signe(leading_coeff(gam)) < 0) gam = ZX_neg(gam);
863 602 : if (cassels_oo_solve(q1, gam) < 0) e = 1;
864 : }
865 64022 : for (i = 1; i < lF; i++)
866 : {
867 61656 : GEN p = gel(F, i);
868 61656 : GEN c = cassels_Qp_solve(q1, gam, p);
869 61656 : if (hilbert(c, a, p) < 0) e = !e;
870 : }
871 2366 : return gc_long(av,e);
872 : }
873 :
874 : static GEN
875 315 : matcassels(GEN F, GEN M)
876 : {
877 315 : long i, j, n = lg(M)-1;
878 315 : GEN C = zero_F2m_copy(n,n);
879 315 : pari_sp av = avma;
880 1974 : for (i = 1; i <= n; i++)
881 : {
882 1659 : GEN Mii = gcoeff(M,i,i);
883 1659 : if (isintzero(Mii)) continue;
884 4606 : for (j = 1; j < i; j++)
885 : {
886 3549 : GEN Mjj = gcoeff(M,j,j);
887 3549 : if (!isintzero(Mjj) && casselspairing(F, Mii, Mjj, gcoeff(M,i,j)))
888 399 : { F2m_set(C,i,j); F2m_set(C,j,i); }
889 : }
890 : }
891 315 : if (DEBUGLEVEL>=2) err_printf("Cassels Pairing: %Ps\n", F2m_to_ZM(C));
892 315 : return gc_const(av, C);
893 : }
894 :
895 : /*******************************************************************/
896 : /* ELLRANK */
897 : /*******************************************************************/
898 : /* This section is a port by Bill Allombert of ellQ.gp by Denis Simon
899 : * Copyright (C) 2019 Denis Simon
900 : * Distributed under the terms of the GNU General Public License (GPL) */
901 :
902 : /* \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */
903 : /* FUNCTIONS FOR POLYNOMIALS \\ */
904 : /* \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */
905 :
906 : static GEN
907 1757 : ell2pol(GEN ell)
908 1757 : { return mkpoln(4, gen_1, ell_get_a2(ell), ell_get_a4(ell), ell_get_a6(ell)); }
909 :
910 : /* find point (x:y:z) on y^2 = pol, return [x,z]~ and set *py = y */
911 : static GEN
912 4504 : projratpointxz(GEN pol, long lim, GEN *py)
913 : {
914 : pari_timer ti;
915 : GEN P;
916 4504 : if (issquareall(leading_coeff(pol), py)) return mkcol2(gen_1, gen_0);
917 4455 : if (DEBUGLEVEL) timer_start(&ti);
918 4455 : P = hyperellratpoints(pol, stoi(lim), 1);
919 4455 : if (DEBUGLEVEL) timer_printf(&ti,"hyperellratpoints(%ld)",lim);
920 4455 : if (lg(P) == 1) return NULL;
921 1736 : P = gel(P,1); *py = gel(P,2); return mkcol2(gel(P,1), gen_1);
922 : }
923 :
924 : /* N = K f^2 where K = core(N) is assumed > 1; let p = P^-(K)
925 : * Return p, K/p and f*p */
926 : static GEN
927 972 : first_core_prime(GEN N, GEN *Kp, GEN *fp)
928 : {
929 972 : GEN v, p = NULL, fa = Z_factor(N), P = gel(fa,1), E = gel(fa,2);
930 972 : long i, j, l = lg(P);
931 972 : for (i = 1; i < l; i++)
932 972 : if (mpodd(gel(E,i))) { p = gel(P,i); break; }
933 972 : v = cgetg(l-1, t_VEC);
934 3751 : for (j = 1, i++; i < l; i++)
935 2779 : if (mpodd(gel(E,i))) gel(v,j++) = gel(P,i);
936 972 : setlg(v, j); *Kp = ZV_prod(v);
937 972 : *fp = mulii(p, sqrtint(diviiexact(N, mulii(*Kp,p))));
938 972 : return p;
939 : }
940 :
941 : /* find point (x:y:z) on y^2 = T, return [x,z]~ and set *py = y */
942 : static GEN
943 1784 : projratpointxz2(GEN T, long lim, long effort, GEN *py)
944 : {
945 1784 : pari_sp av = avma;
946 1784 : GEN list = mkvec(mkvec4(T, matid(2), gen_1, gen_1));
947 1784 : long j, c, ntry = effort * 10;
948 :
949 5035 : for (c = 1; lg(list) > 1 && c <= ntry; c++)
950 : {
951 : GEN K, k, f, p, M, C, r, pol, L;
952 : pari_sp av2;
953 : long lr;
954 :
955 3293 : if (gc_needed(av, 1))
956 : {
957 0 : if (DEBUGMEM > 1)
958 0 : pari_warn(warnmem, "projratpointxz2: #list = %ld",lg(list)-1);
959 0 : list = gc_GEN(av, list);
960 : }
961 3293 : L = gel(list, 1);
962 3293 : list = vecsplice(list, 1);
963 3293 : av2 = avma;
964 3293 : pol = Q_primitive_part(gel(L,1), &K);
965 3293 : M = gel(L,2); /* content = 1 */
966 3293 : K = K? mulii(gel(L,3), K): gel(L,3);
967 3293 : C = gel(L,4);
968 3293 : if (Z_issquareall(K, &k))
969 : {
970 : GEN xz, y, aux, U;
971 3256 : if (c==1) { set_avma(av2); continue; }
972 977 : pol = ZX_hyperellred(pol, &U);
973 977 : if (DEBUGLEVEL>1) err_printf(" reduced quartic: Y^2 = %Ps\n", pol);
974 977 : xz = projratpointxz(pol, lim, &y); if (!xz) { set_avma(av2); continue; }
975 42 : *py = gmul(y, mulii(C, k));
976 42 : aux = RgM_RgC_mul(ZM2_mul(M, U), xz);
977 84 : if (gequal0(gel(aux, 2))) return mkcol2(gel(aux,1), gen_0);
978 42 : *py = gdiv(*py, gpowgs(gel(aux, 2), degpol(pol)>>1));
979 42 : return mkcol2(gdiv(gel(aux, 1), gel(aux, 2)), gen_1);
980 : }
981 972 : p = first_core_prime(K, &K, &f);
982 972 : C = mulii(C, f);
983 : /* root at infinity */
984 1300 : if (dvdii(leading_coeff(pol), p) &&
985 510 : (!dvdii(gcoeff(M,1,1), p) || !dvdii(gcoeff(M,2,1), p)))
986 : {
987 146 : GEN U = mkmat2(gel(M,1), ZC_Z_mul(gel(M,2), p));
988 146 : GEN t = ZX_Z_divexact(ZX_rescale(pol, p), p);
989 146 : list = vec_append(list, mkvec4(t, U, K, C));
990 : }
991 972 : r = FpC_center(FpX_roots(pol, p), p, shifti(p,-1)); lr = lg(r);
992 3028 : for (j = 1; j < lr; j++)
993 : {
994 2056 : pari_sp av3 = avma;
995 2056 : GEN t, U = ZM2_mul(M, mkmat22(p, gel(r,j), gen_0, gen_1));
996 2056 : if (dvdii(gcoeff(U,1,2), p) && dvdii(gcoeff(U,2,2), p))
997 28 : { set_avma(av3); continue; } /* content(U) != 1 (in fact p) */
998 2028 : t = ZX_unscale_div(ZX_translate(pol, gel(r,j)), p);
999 2028 : list = vec_append(list, mkvec4(t, U, K, C));
1000 : }
1001 : }
1002 1742 : return NULL;
1003 : }
1004 :
1005 : static GEN
1006 8540 : polrootsmodpn(GEN pol, GEN p)
1007 : {
1008 8540 : pari_sp av = avma, av2;
1009 8540 : long j, l, i = 1, vd = Z_pval(ZX_disc(pol), p);
1010 : GEN v, r, P;
1011 :
1012 8540 : if (!vd) { set_avma(av); retmkvec(zerovec(2)); }
1013 6314 : pol = Q_primpart(pol);
1014 6314 : P = gpowers0(p, vd-1, p); av2 = avma;
1015 6314 : v = FpX_roots(pol, p); l = lg(v);
1016 18270 : for (j = 1; j < l; j++) gel(v,j) = mkvec2(gel(v,j), gen_1);
1017 72849 : while (i < lg(v))
1018 : {
1019 66535 : GEN pol2, pe, roe = gel(v, i), ro = gel(roe,1);
1020 66535 : long e = itou(gel(roe,2));
1021 :
1022 67529 : if (e >= vd) { i++; continue; }
1023 51282 : pe = gel(P, e);
1024 51282 : (void)ZX_pvalrem(ZX_affine(pol, pe, ro), p, &pol2);
1025 51282 : r = FpX_roots(pol2, p); l = lg(r);
1026 51282 : if (l == 1) { i++; continue; }
1027 104867 : for (j = 1; j < l; j++)
1028 54579 : gel(r, j) = mkvec2(addii(ro, mulii(pe, gel(r, j))), utoi(e + 1));
1029 : /* roots with higher precision = ro + r*p^(e+1) */
1030 50288 : if (l > 2) v = shallowconcat(v, vecslice(r, 2, l-1));
1031 50288 : gel(v, i) = gel(r, 1);
1032 50288 : if (gc_needed(av2, 1)) (void)gc_all(av2, 1, &v);
1033 : }
1034 6314 : if (lg(v) == 1) { set_avma(av); retmkvec(zerovec(2)); }
1035 6314 : return gc_GEN(av, v);
1036 : }
1037 :
1038 : /* \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */
1039 : /* FUNCTIONS FOR LOCAL COMPUTATIONS \\ */
1040 : /* \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */
1041 :
1042 : /* z is integral; sprk true (pr | 2) [t_VEC] or modpr structure (pr | p odd)
1043 : * [t_COL] */
1044 : static GEN
1045 1626751 : kpmodsquares1(GEN nf, GEN z, GEN sprk)
1046 : {
1047 1626751 : GEN modpr = (typ(sprk) == t_VEC)? gmael(sprk, 4, 1): sprk;
1048 1626751 : GEN pr = modpr_get_pr(modpr), p = pr_get_p(pr);
1049 1626751 : long v = nfvalrem(nf, z, pr, &z);
1050 1626751 : if (equaliu(p, 2))
1051 : {
1052 : GEN c;
1053 159810 : z = to_principal_unit(nf, z, pr, sprk); /* = 1 mod pr */
1054 159810 : c = ZV_to_Flv(sprk_log_prk1(nf, z, sprk), 2);
1055 159810 : return vecsmall_prepend(c, odd(v));
1056 : }
1057 : else
1058 : {
1059 1466941 : GEN T = modpr_get_T(modpr);
1060 1466941 : long c = !Fq_issquare(nf_to_Fq(nf, z, modpr), T, p);
1061 1466941 : return mkvecsmall2(odd(v), c);
1062 : }
1063 : }
1064 :
1065 : static GEN
1066 794248 : kpmodsquares(GEN vnf, GEN z, GEN PP)
1067 : {
1068 794248 : pari_sp av = avma;
1069 794248 : long i, j, l = lg(vnf);
1070 794248 : GEN dz, vchar = cgetg(l, t_VEC);
1071 794248 : z = Q_remove_denom(z, &dz); if (dz) z = ZX_Z_mul(z, dz);
1072 1895712 : for (i = 1; i < l; i++)
1073 : {
1074 1101464 : GEN nf = gel(vnf, i), pp = gel(PP, i);
1075 1101464 : GEN kp, delta = ZX_rem(z, nf_get_pol(nf));
1076 1101464 : long lp = lg(pp);
1077 1101464 : kp = cgetg(lp, t_VEC);
1078 2728215 : for (j = 1; j < lp; j++) gel(kp, j) = kpmodsquares1(nf, delta, gel(pp,j));
1079 1101464 : gel(vchar, i) = shallowconcat1(kp);
1080 : }
1081 794248 : return gc_uptoleaf(av, shallowconcat1(vchar));
1082 : }
1083 :
1084 : static GEN
1085 17080 : veckpmodsquares(GEN x, GEN vnf, GEN PP)
1086 789180 : { pari_APPLY_type(t_MAT, kpmodsquares(vnf, gel(x, i), PP)) }
1087 :
1088 : /* \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */
1089 : /* GENERIC FUNCTIONS FOR ELLIPTIC CURVES \\ */
1090 : /* \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */
1091 :
1092 : static GEN
1093 3913 : ellabs(GEN P)
1094 3913 : { return ell_is_inf(P) ? P: mkvec2(gel(P,1), Q_abs_shallow(gel(P,2))); }
1095 : static GEN
1096 4697 : vecellabs(GEN x) { pari_APPLY_same(ellabs(gel(x,i))) }
1097 :
1098 : /* y^2 = x^3 + K a2 x + K^2 a4 x + K^3 a6 */
1099 : static GEN
1100 840 : elltwistequation(GEN ell, GEN K)
1101 : {
1102 : GEN K2, a2, a4, a6;
1103 840 : if (!K || equali1(K)) return ell;
1104 84 : K2 = sqri(K);
1105 84 : a2 = mulii(ell_get_a2(ell), K);
1106 84 : a4 = mulii(ell_get_a4(ell), K2);
1107 84 : a6 = mulii(ell_get_a6(ell), mulii(K, K2));
1108 84 : return ellinit(mkvec5(gen_0, a2, gen_0, a4, a6), NULL, DEFAULTPREC);
1109 : }
1110 :
1111 : /* P=[x,y] a point on Ky^2 = pol(x); [Kx, K^2y] point on Y^2 = K^3pol(X/K) */
1112 : static GEN
1113 224 : elltwistpoint(GEN P, GEN K, GEN K2)
1114 : {
1115 224 : if (ell_is_inf(P)) return ellinf();
1116 224 : return mkvec2(gmul(gel(P,1), K), gmul(gel(P,2), K2));
1117 : }
1118 :
1119 : static GEN
1120 1673 : elltwistpoints(GEN x, GEN K)
1121 : {
1122 : GEN K2;
1123 1673 : if (!K || gequal1(K)) return x;
1124 126 : K2 = gsqr(K);
1125 350 : pari_APPLY_same(elltwistpoint(gel(x,i), K, K2))
1126 : }
1127 :
1128 : /* \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */
1129 : /* FUNCTIONS FOR NUMBER FIELDS \\ */
1130 : /* \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */
1131 :
1132 : /* return a set S2 of prime ideals disjoint from S such that
1133 : * Cl_{S \cup S2}(K) has no p-torsion */
1134 : static GEN
1135 441 : bestS(GEN bnf,GEN S, ulong p)
1136 : {
1137 441 : GEN v, S2, h = bnf_get_no(bnf), cyc = bnf_get_cyc(bnf);
1138 : long i, lS2;
1139 : ulong l, vD;
1140 : forprime_t P;
1141 :
1142 441 : if (!dvdiu(h, p)) return cgetg(1,t_VEC);
1143 322 : if (!S)
1144 : {
1145 0 : v = diagonal_shallow(cyc);
1146 0 : vD = Z_lval(h, p);
1147 : }
1148 : else
1149 : {
1150 322 : long lS = lg(S);
1151 322 : v = cgetg(lS,t_MAT);
1152 1022 : for (i = 1; i < lS; i++) gel(v,i) = isprincipal(bnf, gel(S,i));
1153 322 : v = ZM_hnfmodid(v, cyc);
1154 322 : vD = Z_lval(ZM_det(v), p); if (!vD) return cgetg(1, t_VEC);
1155 : }
1156 301 : S2 = cgetg(vD+2, t_VEC); lS2 = 1;
1157 301 : u_forprime_init(&P,2,ULONG_MAX);
1158 2569 : while ((l = u_forprime_next(&P)))
1159 : {
1160 2569 : pari_sp av = avma;
1161 2569 : GEN w, Sl = idealprimedec(bnf, utoi(l));
1162 2569 : long nSl = lg(Sl)-1;
1163 : ulong vDl;
1164 4074 : for (i = 1; i < nSl; i++) /* remove one prime ideal */
1165 : {
1166 1806 : w = ZM_hnf(shallowconcat(v, isprincipal(bnf, gel(Sl,i))));
1167 1806 : vDl = Z_lval(ZM_det(w), p);
1168 1806 : if (vDl < vD)
1169 : {
1170 1239 : gel(S2,lS2++) = gel(Sl,i);
1171 1239 : vD = vDl; v = w; av = avma;
1172 1239 : if (!vD) { setlg(S2, lS2); return S2; }
1173 : }
1174 : }
1175 2268 : set_avma(av);
1176 : }
1177 : return NULL;/*LCOV_EXCL_LINE*/
1178 : }
1179 :
1180 : static GEN
1181 882 : nfC_prV_val(GEN nf, GEN G, GEN P)
1182 : {
1183 882 : long i, j, lG = lg(G), lP = lg(P);
1184 882 : GEN M = cgetg(lG, t_MAT);
1185 5418 : for (i = 1; i < lG; i++)
1186 : {
1187 4536 : GEN V = cgetg(lP, t_COL);
1188 19390 : for (j = 1; j < lP; j++)
1189 14854 : gel(V,j) = gpnfvalrem(nf, gel(G,i), gel(P,j), NULL);
1190 4536 : gel(M,i) = V;
1191 : }
1192 882 : return M;
1193 : }
1194 :
1195 : static GEN
1196 5614 : _factorbackmod(GEN nf, GEN g, GEN e, ulong p)
1197 : {
1198 5614 : GEN y = nffactorback(nf, g, e), den;
1199 5614 : GEN z = nfmul(nf, y, nfsqr(nf, idealredmodpower(nf, y, p, 0)));
1200 5614 : z = Q_remove_denom(z, &den);
1201 5614 : if (den)
1202 : {
1203 0 : if (p != 2) den = powiu(den, p-1);
1204 0 : z = gmul(z, den);
1205 : }
1206 5614 : return z;
1207 : }
1208 : static GEN
1209 441 : nfV_factorbackmod(GEN nf, GEN x, ulong p)
1210 : {
1211 441 : long i, l = lg(x);
1212 441 : GEN v = cgetg(l, t_VEC);
1213 3794 : for (i = 1; i < l; i++)
1214 : {
1215 3353 : GEN y = gel(x,i), g = gel(y,1), e = gel(y,2);
1216 3353 : gel(v,i) = _factorbackmod(nf, g, ZV_to_Flv(e,p), p);
1217 : }
1218 441 : return v;
1219 : }
1220 : static GEN
1221 441 : nfV_zm_factorback(GEN nf, GEN G, GEN x, long p)
1222 2702 : { pari_APPLY_type(t_VEC, _factorbackmod(nf, G, gel(x,i), p)) }
1223 :
1224 : static GEN
1225 441 : prV_ZM_factorback(GEN nf, GEN S, GEN x)
1226 2702 : { pari_APPLY_type(t_VEC,idealfactorback(nf, S, gel(x,i), 0)) }
1227 :
1228 : /* shortcut for bnf = Q and p = 2 */
1229 : static GEN
1230 826 : bnfselmerQ(GEN S)
1231 : {
1232 826 : GEN g = vec_prepend(prV_primes(S), gen_m1), e;
1233 826 : long n = lg(S)-1;
1234 826 : e = n? shallowconcat(zerocol(n), matid(n)): zeromat(0, 1);
1235 826 : return mkvec3(g, e, const_vec(n + 1, gen_1));
1236 : }
1237 :
1238 : static GEN
1239 441 : bnfselmer(GEN bnf, GEN S)
1240 : {
1241 441 : const long p = 2;
1242 441 : pari_sp av = avma;
1243 441 : GEN nf = bnf_get_nf(bnf), S2, S3, e, f, e2, kerval, LS2gen, LS2fu, LS2all;
1244 441 : long n = lg(S)-1, n3, n2all, r;
1245 :
1246 441 : S2 = bestS(bnf, S, p);
1247 441 : S3 = shallowconcat(S, S2);
1248 441 : LS2all = nfV_factorbackmod(nf, gel(bnfunits(bnf, S3), 1), p);
1249 441 : n3 = lg(S3)-1; n2all = lg(LS2all)-1;
1250 441 : LS2gen = vecslice(LS2all,1,n3);
1251 441 : LS2fu = vecslice(LS2all,n3+1, n2all-1);
1252 441 : e2 = nfC_prV_val(nf, LS2gen, S2);
1253 441 : kerval = Flm_ker(ZM_to_Flm(e2, p), p);
1254 441 : LS2gen = nfV_zm_factorback(nf, LS2gen, kerval, p);
1255 441 : e = nfC_prV_val(nf, LS2gen, S);
1256 441 : e2 = ZM_divexactu(ZM_zm_mul(e2, kerval), p);
1257 441 : f = prV_ZM_factorback(nf, S2, e2);
1258 441 : LS2gen = shallowconcat(LS2fu, LS2gen);
1259 441 : LS2gen = nfV_to_scalar_or_alg(nf, vec_prepend(LS2gen, bnf_get_tuU(bnf)));
1260 441 : r = n2all - n3;
1261 441 : e = shallowconcat(zeromat(n, r), e);
1262 441 : f = shallowconcat(const_vec(r, gen_1), f);
1263 441 : return gc_GEN(av, mkvec3(LS2gen,e,f));
1264 : }
1265 :
1266 : static GEN
1267 259 : get_kerval(GEN nf, GEN S2, GEN LS2gen)
1268 : {
1269 259 : long i, j, lS2 = lg(S2), l = lg(LS2gen);
1270 259 : GEN e = cgetg(l, t_MAT);
1271 3864 : for (i = 1; i < l; i++)
1272 : {
1273 3605 : GEN v = cgetg(lS2, t_VECSMALL);
1274 25739 : for (j=1; j < lS2; j++) v[j] = odd(idealval(nf, gel(LS2gen, i), gel(S2,j)));
1275 3605 : gel(e, i) = Flv_to_F2v(v);
1276 : }
1277 259 : return F2m_ker(e);
1278 : }
1279 : static GEN
1280 259 : nf2selmer_quad(GEN nf, GEN S)
1281 : {
1282 259 : pari_sp ltop = avma;
1283 259 : GEN D = nf_get_disc(nf), factD = nf_get_ramified_primes(nf);
1284 259 : GEN SlistQ = prV_primes(S), QS2gen, gen, Hlist, H, KerH, LS2gen;
1285 : GEN chpol, Q, kerval, S2, G, e, f, b, c, bad;
1286 259 : long lS = lg(S), l, lHlist, i, j, k;
1287 259 : GEN nfa = nf_get_ramified_primes(nf);
1288 :
1289 259 : QS2gen = vec_prepend(SlistQ, gen_m1);
1290 259 : bad = ZV_sort_uniq_shallow(shallowconcat(factD, SlistQ));
1291 259 : Hlist = ZV_search(bad, gen_2)? bad: vec_prepend(bad, gen_2);
1292 259 : l = lg(QS2gen);
1293 259 : lHlist = lg(Hlist);
1294 259 : H = cgetg(l, t_MAT);
1295 2450 : for (j = 1; j < l; j++)
1296 : {
1297 2191 : GEN v = cgetg(lHlist, t_VECSMALL);
1298 33705 : for (i = 1; i < lHlist; i++)
1299 31514 : v[i] = hilbertii(D, gel(QS2gen, j), gel(Hlist, i)) < 0;
1300 2191 : gel(H, j) = Flv_to_F2v(v);
1301 : }
1302 259 : KerH = F2m_ker(H); l = lg(KerH);
1303 259 : LS2gen = cgetg(l, t_VEC);
1304 259 : chpol = QXQ_charpoly(gel(nf_get_zk(nf), 2), nf_get_pol(nf), 0);
1305 259 : b = negi(gel(chpol, 3));
1306 259 : c = shifti(gel(chpol, 2),1);
1307 259 : Q = mkmat3(mkcol3(gen_2, b, gen_0),
1308 : mkcol3(b, c, gen_0),
1309 : mkcol3(gen_0, gen_0, gen_0));
1310 1470 : for (k = 1; k < l; k++)
1311 : {
1312 1211 : GEN P = RgV_F2v_extract_shallow(QS2gen, gel(KerH, k));
1313 1211 : GEN F = ZV_sort_uniq_shallow(shallowconcat(nfa, P)), sol;
1314 1211 : gcoeff(Q, 3, 3) = mulis(ZV_prod(P), -2);
1315 1211 : sol = qfsolve(mkvec2(Q, F)); /* must be solvable */
1316 1211 : sol = Q_primpart(mkcol2(gel(sol,1), gel(sol,2)));
1317 1211 : gel(LS2gen, k) = basistoalg(nf, sol);
1318 : }
1319 259 : if (equalis(D, -4)) G = bad;
1320 : else
1321 : {
1322 259 : G = vecsplice(bad, ZV_search(bad, veclast(factD)));
1323 259 : G = vec_prepend(G, gen_m1);
1324 : }
1325 259 : LS2gen = shallowconcat(G, LS2gen);
1326 259 : l = lg(SlistQ); S2 = cgetg(l, t_VEC);
1327 259 : if (l > 1)
1328 : {
1329 2184 : for (i = 1; i < l; i++) gel(S2, i) = idealprimedec(nf, gel(SlistQ, i));
1330 252 : S2 = setminus(shallowconcat1(S2), S);
1331 : }
1332 259 : kerval = get_kerval(nf, S2, LS2gen); l = lg(kerval);
1333 259 : gen = cgetg(l, t_VEC);
1334 259 : e = cgetg(l, t_MAT);
1335 259 : f = cgetg(l, t_VEC);
1336 2905 : for (i = 1; i < l; i++)
1337 : {
1338 2646 : GEN id, ei, x = nffactorback(nf, LS2gen, F2v_to_Flv(gel(kerval, i)));
1339 2646 : gel(e,i) = ei = cgetg(lS, t_COL);
1340 40726 : for (j = 1; j < lS; j++) gel(ei,j) = stoi(idealval(nf, x, gel(S,j)));
1341 2646 : id = idealdiv(nf, x, idealfactorback(nf, S, gel(e,i), 0));
1342 2646 : if (!idealispower(nf, id, 2, &gel(f,i))) pari_err_BUG("nf2selmer_quad");
1343 2646 : gel(gen, i) = nf_to_scalar_or_alg(nf, x);
1344 : }
1345 259 : return gc_GEN(ltop, mkvec3(gen, e, f));
1346 : }
1347 :
1348 : static GEN
1349 868 : makevbnf(GEN ell, long prec)
1350 : {
1351 868 : GEN vbnf, P = gel(ZX_factor(ell2pol(ell)), 1);
1352 868 : long k, l = lg(P);
1353 868 : vbnf = cgetg(l, t_VEC);
1354 2387 : for (k = 1; k < l; k++)
1355 : {
1356 1519 : GEN t = gel(P,k);
1357 1519 : gel(vbnf,k) = degpol(t) <= 2? nfinit(t, prec): Buchall(t, nf_FORCE, prec);
1358 : }
1359 868 : return vbnf;
1360 : }
1361 :
1362 : static GEN
1363 889 : kernorm(GEN G, GEN S, GEN pol, GEN signs)
1364 : {
1365 889 : long i, j, lS = lg(S), lG = lg(G), lv = signs? lS+2: lS+1;
1366 889 : GEN M = cgetg(lG, t_MAT);
1367 14315 : for (j = 1; j < lG; j++)
1368 : {
1369 13426 : GEN v, N = QXQ_norm(gel(G,j), pol);
1370 13426 : gel(M, j) = v = cgetg(lv, t_VECSMALL);
1371 13426 : v[1] = gsigne(N) < 0;
1372 191506 : for (i = 1; i < lS; i++) v[i+1] = odd(Q_pvalrem(N, gel(S,i), &N));
1373 13426 : if (signs) v[i+1] = signs[j];
1374 : }
1375 889 : return Flm_ker(M, 2);
1376 : }
1377 :
1378 : /* \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */
1379 : /* FUNCTIONS FOR 2-DESCENT \\ */
1380 : /* \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */
1381 : /* vector of t_VEC; return total number of entries */
1382 : static long
1383 8540 : RgVV_nb(GEN v)
1384 : {
1385 8540 : long i, l = lg(v), n = 0;
1386 24794 : for (i = 1; i < l; i++) n += lg(gel(v,i)) - 1;
1387 8540 : return n;
1388 : }
1389 :
1390 : /* return an Fp basis */
1391 : static GEN
1392 8540 : elllocalimage(GEN pol, GEN K, GEN vnf, GEN p, GEN pp, GEN pts)
1393 : {
1394 8540 : long n, p2 = RgVV_nb(pp), prank = equaliu(p, 2)? p2: p2 - 1;
1395 8540 : GEN R = polrootsmodpn(pol, p), bound = addiu(p, 6);
1396 :
1397 8540 : for(n = 1;; n++)
1398 22148 : {
1399 : pari_sp btop;
1400 : GEN x, y2, d;
1401 30688 : pts = Flm_image(pts, 2); if (lg(pts)-1 == prank) break;
1402 22148 : if ((n & 0xf) == 0) bound = mulii(bound, p);
1403 22148 : btop = avma;
1404 : do
1405 : {
1406 107737 : GEN r = gel(R, random_Fl(lg(R)-1)+1);
1407 107737 : long pprec = random_Fl(itou(gel(r, 2)) + 3) - 2; /* >= -2 */
1408 107737 : set_avma(btop);
1409 107737 : x = gadd(gel(r, 1), gmul(powis(p, pprec), randomi(bound)));
1410 107737 : y2 = gmul(K, poleval(pol, x));
1411 107737 : } while (gequal0(y2) || !Qp_issquare(y2, p));
1412 22148 : d = deg1pol_shallow(negi(K), gmul(K, x), 0);
1413 22148 : pts = vec_append(pts, kpmodsquares(vnf, d, pp));
1414 : }
1415 8540 : return pts;
1416 : }
1417 :
1418 : static GEN
1419 889 : ellLS2image(GEN pol, GEN vP, GEN K, GEN vpol, GEN vcrt)
1420 : {
1421 889 : long i, l = lg(vP);
1422 : GEN v;
1423 :
1424 889 : if (l == 1) return cgetg(1, t_VEC);
1425 756 : v = cgetg(l, t_VEC);
1426 60263 : for (i = 1; i < l; i++)
1427 : {
1428 59507 : GEN P = gel(vP, i), x, t;
1429 59507 : if (ell_is_inf(P)) { gel(v, i) = gen_1; continue; }
1430 59507 : x = gel(P,1);
1431 59507 : t = deg1pol_shallow(negi(K), gmul(K, x), 0);
1432 59507 : if (gequal0(gel(P,2)))
1433 : { /* 2-torsion, x now integer and a root of exactly one linear vpol[k]=T */
1434 623 : long k, lp = lg(vpol);
1435 : GEN a;
1436 1099 : for (k = 1; k < lp; k++)
1437 : {
1438 1099 : GEN T = gel(vpol,k), z = gel(T,2);
1439 1099 : if (absequalii(x, z) && signe(x) == -signe(z)) break; /* T = X-x */
1440 : }
1441 623 : a = ZX_Z_eval(ZX_deriv(pol), x);
1442 623 : t = gadd(a, gmul(gel(vcrt,k), gsub(t, a))); /* a mod T, t mod pol/T*/
1443 : }
1444 59507 : gel(v, i) = t;
1445 : }
1446 756 : return v;
1447 : }
1448 :
1449 : /* find small points on ell; 2-torsion points must be returned first */
1450 : static GEN
1451 889 : ellsearchtrivialpoints(GEN ell, GEN lim, GEN help)
1452 : {
1453 889 : pari_sp av = avma;
1454 889 : GEN tors2 = gel(elltors_psylow(ell,2),3);
1455 889 : GEN triv = lim ? ellratpoints(ell, lim, 0): cgetg(1,t_VEC);
1456 889 : if (help) triv = shallowconcat(triv, help);
1457 889 : return gc_GEN(av, shallowconcat(tors2, triv));
1458 : }
1459 :
1460 : GEN
1461 63 : ellrankinit(GEN ell, long prec)
1462 : {
1463 63 : pari_sp av = avma;
1464 : GEN urst;
1465 63 : checkell_Q(ell); ell = ellminimalbmodel(ell, &urst);
1466 63 : return gc_GEN(av, mkvec3(ell, urst, makevbnf(ell, prec)));
1467 : }
1468 :
1469 : INLINE GEN
1470 256109 : ZV_isneg(GEN x) { pari_APPLY_long(signe(gel(x,i)) < 0) }
1471 :
1472 : static GEN
1473 72933 : selmersign(GEN x, GEN vpol, GEN L)
1474 164115 : { pari_APPLY_same(ZV_isneg(nfeltsign(gel(x, i), RgX_rem(L, gel(vpol, i)), NULL))) }
1475 :
1476 : static GEN
1477 1778 : matselmersign(GEN vnf, GEN vpol, GEN x)
1478 74711 : { pari_APPLY_type(t_MAT, shallowconcat1(selmersign(vnf, vpol, gel(x, i)))) }
1479 :
1480 : static GEN
1481 118974 : _trace(GEN z, GEN T)
1482 : {
1483 118974 : long n = degpol(T)-1;
1484 118974 : if (degpol(z) < n) return gen_0;
1485 69236 : return gdiv(gel(z, 2+n), gel(T, 3+n));
1486 : }
1487 : static GEN
1488 19829 : tracematrix(GEN zc, GEN b, GEN T)
1489 : {
1490 : long i, j;
1491 19829 : GEN q = cgetg(4, t_MAT);
1492 79316 : for (j = 1; j <= 3; j++) gel(q,j) = cgetg(4, t_COL);
1493 79316 : for (j = 1; j <= 3; j++)
1494 : {
1495 118974 : for (i = 1; i < j; i++) gcoeff(q,i,j) = gcoeff(q,j,i) =
1496 59487 : _trace(QXQ_mul(zc, QXQ_mul(gel(b,i), gel(b,j), T), T), T);
1497 59487 : gcoeff(q,i,i) = _trace(QXQ_mul(zc, QXQ_sqr(gel(b,i), T), T), T);
1498 : }
1499 19829 : return q;
1500 : }
1501 :
1502 : static GEN
1503 21630 : RgXV_cxeval(GEN x, GEN r, GEN ri)
1504 86520 : { pari_APPLY_same(RgX_cxeval(gel(x,i), r, ri)) }
1505 :
1506 : static GEN
1507 7202 : redquadric(GEN base, GEN pol, GEN zc)
1508 : {
1509 7202 : pari_sp av = avma;
1510 7202 : long prec = nbits2prec(gexpo(pol)+gexpo(zc)) + EXTRAPRECWORD;
1511 : for (;;)
1512 8 : {
1513 7210 : GEN R = roots(pol, prec), s = NULL;
1514 7210 : long i, l = lg(R);
1515 28840 : for (i = 1; i < l; ++i)
1516 : {
1517 21630 : GEN r = gel(R,i), ri = gexpo(r) > 1? ginv(r): NULL;
1518 21630 : GEN b = RgXV_cxeval(base, r, ri), z = RgX_cxeval(zc, r, ri);
1519 21630 : GEN M = RgC_RgV_mulrealsym(RgV_Rg_mul(b, gabs(z, prec)), gconj(b));
1520 21630 : s = s? RgM_add(s, M): M;
1521 : }
1522 7210 : s = RgM_Cholesky(s, prec);
1523 7210 : if (s) return gc_upto(av, lll(s));
1524 8 : prec = precdbl(prec); set_avma(av);
1525 : }
1526 : }
1527 :
1528 : static GEN
1529 19915 : RgX_homogenous_evaldeg(GEN P, GEN A, GEN B)
1530 : {
1531 19915 : long i, d = degpol(P), e = lg(B)-1;
1532 19915 : GEN s = gmul(gel(P, d+2), gel(B,e-d));
1533 62377 : for (i = d-1; i >= 0; i--)
1534 42462 : s = gadd(gmul(s, A), gmul(gel(B,e-i), gel(P,i+2)));
1535 19915 : return s;
1536 : }
1537 :
1538 : static GEN
1539 5425 : RgXV_homogenous_evaldeg(GEN x, GEN a, GEN b)
1540 21700 : { pari_APPLY_same(RgX_homogenous_evaldeg(gel(x,i), a, b)) }
1541 :
1542 : static void
1543 42 : check_oncurve(GEN ell, GEN v)
1544 : {
1545 42 : long i, l = lg(v);
1546 49 : for (i = 1; i < l; i++)
1547 : {
1548 7 : GEN P = gel(v, i);
1549 7 : if (!checkellpt_i(P) || ell_is_inf(P)
1550 7 : || !RgV_is_QV(P)
1551 7 : || !ellisoncurve_i(ell,P)) pari_err_TYPE("ellrank",P);
1552 : }
1553 42 : }
1554 :
1555 : static GEN
1556 19178 : basis(GEN nf, GEN x, GEN crt, GEN pol)
1557 : {
1558 19178 : long i, l = lg(x);
1559 19178 : GEN b = cgetg(l, t_COL);
1560 42884 : for (i = 1; i < l; i++)
1561 : {
1562 23706 : GEN z = nf_to_scalar_or_alg(nf, gel(x, i));
1563 23706 : gel(b, i) = grem(gsub(z, gmul(crt, z)), pol); /* z mod T, 0 mod (pol/T) */
1564 : }
1565 19178 : return b;
1566 : }
1567 :
1568 : static GEN
1569 889 : vecsmallbasis(GEN x, GEN vcrt, GEN pol)
1570 2415 : { pari_APPLY_same(basis(gel(x,i), nf_get_zk(gel(x,i)), gel(vcrt,i), pol)) }
1571 :
1572 : static GEN
1573 270372 : ZC_shifti(GEN x, long n) { pari_APPLY_type(t_COL, shifti(gel(x,i), n)) }
1574 :
1575 : /* true nf */
1576 : static GEN
1577 17652 : selmerbasis(GEN nf, GEN ek, GEN sqrtLS2, GEN factLS2,
1578 : GEN badprimes, GEN crt, GEN pol)
1579 : {
1580 17652 : GEN sqrtzc = idealfactorback(nf, sqrtLS2, zv_neg(ek), 0);
1581 17652 : GEN E = ZC_shifti(ZM_zc_mul(factLS2, ek), -1);
1582 :
1583 17652 : if (ZV_equal0(E))
1584 15812 : sqrtzc = idealhnf_shallow(nf, sqrtzc);
1585 : else
1586 1840 : sqrtzc = idealmul(nf, sqrtzc, idealfactorback(nf, badprimes, ZC_neg(E), 0));
1587 17652 : return basis(nf, sqrtzc, crt, pol);
1588 : }
1589 :
1590 567 : static long randu(void) { return random_Fl(127) - 63; }
1591 : static GEN
1592 189 : randS(GEN b)
1593 : {
1594 567 : return gadd(gmulgs(gel(b,1), randu()),
1595 378 : gadd(gmulgs(gel(b,2), randu()), gmulgs(gel(b,3), randu())));
1596 : }
1597 :
1598 : static GEN
1599 7013 : liftselmerinit(GEN expo, GEN vnf, GEN sqrtLS2, GEN factLS2,
1600 : GEN badprimes, GEN vcrt, GEN pol)
1601 : {
1602 7013 : long n = lg(vnf)-1, k, t;
1603 7013 : GEN b = cgetg(n+1, t_VEC);
1604 24665 : for (k = t = 1; k <= n; k++)
1605 : {
1606 17652 : GEN fak = gel(factLS2,k), ek;
1607 17652 : long m = lg(fak)-1;
1608 17652 : ek = vecslice(expo, t, t + m-1); t += m;
1609 17652 : gel(b,k) = selmerbasis(gel(vnf,k), ek, gel(sqrtLS2,k),
1610 17652 : fak, gel(badprimes,k), gel(vcrt,k), pol);
1611 : }
1612 7013 : return shallowconcat1(b);
1613 : }
1614 :
1615 : static GEN
1616 7202 : qf_disc_fact(GEN q, GEN L)
1617 : {
1618 7202 : GEN D = ZM_det(q), P, E;
1619 7202 : GEN H = Z_smoothen(D, L, &P, &E);
1620 7202 : if (H)
1621 455 : P = shallowconcat(P, gel(Z_factor(H),1));
1622 7202 : return mkvec2(q, P);
1623 : }
1624 :
1625 : static GEN
1626 3640 : liftselmer_cover(GEN b, GEN expo, GEN LS2, GEN pol, GEN discF, GEN K)
1627 : {
1628 : GEN P, Q, Q4, R, den, q0, q1, q2, xz, x, y, y2m, U, param, newb;
1629 : GEN ttheta, tttheta, zc, polprime;
1630 : GEN QM, zden;
1631 3640 : zc = RgXQV_factorback(LS2, expo, pol);
1632 3640 : if (typ(zc)==t_INT) zc = scalarpol(zc, varn(pol));
1633 3640 : ttheta = RgX_shift_shallow(pol,-2);
1634 3640 : tttheta = RgX_shift_shallow(pol, -1);
1635 3640 : polprime = ZX_deriv(pol);
1636 3640 : q2 = Q_primpart(tracematrix(zc, b, pol));
1637 3640 : U = redquadric(b, pol, QXQ_div(zc, polprime, pol));
1638 3640 : q2 = qf_RgM_apply(q2, U);
1639 3640 : newb = RgV_RgM_mul(b, U);
1640 3640 : param = Q_primpart(qfparam(q2, qfsolve(qf_disc_fact(q2,gel(discF,2))), 1));
1641 3640 : param = RgM_to_RgXV_reverse(shallowtrans(param), 0);
1642 3640 : q1 = RgM_neg(tracematrix(QXQ_mul(zc, ttheta, pol), newb, pol));
1643 3640 : q1 = Q_remove_denom(qfeval(q1, param), &den);
1644 3640 : if (den) q1 = ZX_Z_mul(q1, den);
1645 3640 : if (!equali1(K)) q1 = RgX_Rg_mul(q1, K);
1646 3640 : QM = quartic_minim_all(q1, discF);
1647 3640 : q1 = gel(QM,1);
1648 3640 : zden = gmael(QM,2,1);
1649 3640 : Q = ZX_hyperellred(q1, &R);
1650 3640 : R = gmul(gmael(QM,2,2), R);
1651 3640 : if (DEBUGLEVEL>1) err_printf(" reduced quartic: Y^2 = %Ps\n", Q);
1652 3640 : xz = mkcol2(pol_x(0),gen_1);
1653 3640 : P = RgM_RgC_mul(R, xz); x = gel(P,1); y = gel(P,2);
1654 3640 : param = RgXV_homogenous_evaldeg(param, x, gpowers(y, 2));
1655 3640 : param = gmul(param, gdiv(den? mulii(K, den): K, zden));
1656 3640 : q0 = tracematrix(QXQ_mul(zc, tttheta, pol), newb, pol);
1657 3640 : x = gdiv(qfeval(q0, param), K);
1658 3640 : Q4 = gpowers(Q,4);
1659 3640 : y2m = gmul(RgX_homogenous_evaldeg(pol, x, Q4), Q);
1660 3640 : if (!issquareall(gdiv(y2m, K), &y))
1661 0 : pari_err_BUG("liftselmer_cover");
1662 3640 : y = gdiv(y, gel(Q4,2));
1663 3640 : P = mkvec2(gdiv(gmul(K,x),pol_xn(2,1)),gdiv(gmul(gsqr(K),y),pol_xn(3,1)));
1664 3640 : return mkvec2(Q,P);
1665 : }
1666 :
1667 : static GEN
1668 3373 : liftselmer(GEN b, GEN expo, GEN sbase, GEN LS2, GEN pol, GEN discF, GEN K, long ntry, GEN *pt_Q)
1669 : {
1670 3373 : pari_sp av = avma, av2;
1671 3373 : long t, lim = ntry * LIM1;
1672 : GEN ttheta, tttheta, z, polprime;
1673 : hashtable h;
1674 3373 : hash_init_GEN(&h, ntry, ZX_equal, 1);
1675 3373 : z = RgXQV_factorback(LS2, expo, pol);
1676 3373 : ttheta = RgX_shift_shallow(pol,-2);
1677 3373 : tttheta = RgX_shift_shallow(pol, -1);
1678 3373 : polprime = ZX_deriv(pol); av2 = avma;
1679 4051 : for (t = 1; t <= ntry; t++, set_avma(av2))
1680 : {
1681 : GEN P, Q, Qk, R, den, q0, q1, q2, xz, x, y, zz, zc, U, param, newb, zden, QM;
1682 : long idx;
1683 3562 : if (t == 1) zc = z;
1684 : else
1685 : {
1686 : GEN r;
1687 189 : do r = randS(sbase); while (degpol(QX_gcd(r, pol)));
1688 189 : zc = QXQ_mul(z, QXQ_sqr(r, pol), pol);
1689 : }
1690 3562 : q2 = Q_primpart(tracematrix(zc, b, pol));
1691 3562 : U = redquadric(b, pol, QXQ_div(zc, polprime, pol));
1692 4240 : if (lg(U) < 4) continue;
1693 3562 : q2 = qf_RgM_apply(q2, U);
1694 3562 : newb = RgV_RgM_mul(b, U);
1695 3562 : param = Q_primpart(qfparam(q2, qfsolve(qf_disc_fact(q2,gel(discF,2))), 1));
1696 3562 : param = RgM_to_RgXV_reverse(shallowtrans(param), 0);
1697 3562 : q1 = RgM_neg(tracematrix(QXQ_mul(zc, ttheta, pol), newb, pol));
1698 3562 : q1 = Q_remove_denom(qfeval(q1, param), &den);
1699 3562 : if (den) q1 = ZX_Z_mul(q1, den);
1700 3562 : if (!equali1(K)) q1 = RgX_Rg_mul(q1, K);
1701 3562 : QM = quartic_minim_all(q1, discF);
1702 3562 : q1 = gel(QM,1);
1703 3562 : zden = gmael(QM,2,1);
1704 3562 : Q = ZX_hyperellred(q1, &R);
1705 3562 : R = gmul(gmael(QM,2,2), R);
1706 3562 : if (pt_Q) *pt_Q = Q;
1707 3562 : Qk = shallowcopy(Q);
1708 3562 : (void) ZX_canon_neg(Qk);
1709 3562 : if (hash_haskey_long(&h, (void*)Qk, &idx)) continue;
1710 3527 : hash_insert_long(&h, (void*)Qk, 1); av2 = avma;
1711 3527 : if (DEBUGLEVEL>1) err_printf(" reduced quartic: Y^2 = %Ps\n", Q);
1712 :
1713 3527 : xz = projratpointxz(Q, lim, &zz);
1714 3527 : if (!xz)
1715 : {
1716 1784 : xz = projratpointxz2(Q, lim, ntry, &zz);
1717 1784 : if (!xz)
1718 : {
1719 3527 : if (pt_Q) return NULL; else continue;
1720 : }
1721 : }
1722 1785 : P = RgM_RgC_mul(R, xz); x = gel(P,1); y = gel(P,2);
1723 1785 : param = RgXV_homogenous_evaldeg(param, x, gpowers(y, 2));
1724 1785 : param = gmul(param, gdiv(den? mulii(K, den): K, gmul(zz, zden)));
1725 1785 : q0 = tracematrix(QXQ_mul(zc, tttheta, pol), newb, pol);
1726 1785 : x = gdiv(qfeval(q0, param), K);
1727 1785 : if (!issquareall(gdiv(poleval(pol, x), K), &y)) /* K y^2 = pol(x) */
1728 0 : pari_err_BUG("ellrank");
1729 1785 : P = mkvec2(x, y);
1730 1785 : if (DEBUGLEVEL) err_printf("Found point: %Ps\n", P);
1731 1785 : if (pt_Q) *pt_Q = gen_0;
1732 1785 : return gc_GEN(av, P);
1733 : }
1734 489 : return NULL;
1735 : }
1736 :
1737 : static void
1738 784 : gtoset_inplace(GEN x)
1739 784 : { gen_sort_inplace(x, (void*)&cmp_universal, cmp_nodata, NULL); }
1740 :
1741 : /* FIXME: export */
1742 : static void
1743 8540 : setlgall(GEN x, long L)
1744 : {
1745 8540 : long i, l = lg(x);
1746 188972 : for(i = 1; i < l; i++) setlg(gel(x,i), L);
1747 8540 : }
1748 :
1749 : static long
1750 8540 : dim_selmer(GEN p, GEN pol, GEN K, GEN vnf, GEN LS2, GEN helpLS2,
1751 : GEN *selmer, GEN *LS2chars, GEN *helpchars)
1752 : {
1753 : pari_sp av;
1754 8540 : long dim, k, lvnf = lg(vnf);
1755 8540 : GEN X, L, LS2image, helpimage, pp = cgetg(lvnf, t_VEC);
1756 8540 : int pis2 = equaliu(p, 2);
1757 :
1758 24794 : for (k = 1; k < lvnf; k++)
1759 : {
1760 16254 : GEN v, nf = gel(vnf,k), PR = idealprimedec(nf, p);
1761 16254 : long j, l = lg(PR);
1762 16254 : gel(pp, k) = v = cgetg(l, t_VEC);
1763 36358 : for (j = 1; j < l; j++)
1764 : {
1765 20104 : GEN pr = gel(PR,j);
1766 20104 : gel(v,j) = pis2? log_prk_init(nf, pr, 1 + 2 * pr_get_e(pr), NULL)
1767 20104 : : zkmodprinit(nf, pr);
1768 : }
1769 : }
1770 8540 : LS2image = veckpmodsquares(LS2, vnf, pp);
1771 8540 : *LS2chars = vconcat(*LS2chars, LS2image);
1772 8540 : helpimage = veckpmodsquares(helpLS2, vnf, pp);
1773 8540 : *helpchars = vconcat(*helpchars, helpimage);
1774 8540 : av = avma;
1775 8540 : L = elllocalimage(pol, K, vnf, p, pp, helpimage);
1776 8540 : X = Flm_ker(shallowconcat(LS2image, L), 2); setlgall(X, lg(LS2image));
1777 : /* intersect(LS2image, locim) = LS2image.X */
1778 8540 : *selmer = Flm_intersect_i(*selmer, shallowconcat(Flm_ker(LS2image,2), X), 2);
1779 8540 : *selmer = gc_upto(av, Flm_image(*selmer, 2));
1780 8540 : dim = lg(*selmer)-1; return (dim == Flm_rank(helpimage,2))? dim: -1;
1781 : }
1782 :
1783 : /* Assume there are 3 real roots, if K>0 return the smallest, otherwise the largest */
1784 : static long
1785 490 : get_row(GEN vnf, GEN K)
1786 : {
1787 490 : long k, sK = signe(K), n = lg(vnf)-1;
1788 : GEN R;
1789 490 : if (n == 1) return sK > 0? 1: 3;
1790 294 : if (n == 2)
1791 : {
1792 105 : GEN P = nf_get_pol(gel(vnf,2));
1793 105 : GEN z = negi(constant_coeff(nf_get_pol(gel(vnf,1))));
1794 105 : GEN y = poleval(P,z);
1795 105 : GEN b = gel(P,3), a = gel(P,4);
1796 105 : if (signe(y) != signe(a))
1797 : /* 1 is between 2 and 3 */
1798 7 : return sK > 0? 2: 3;
1799 98 : else if (cmpii(mulii(z,mulis(a,-2)), b) == signe(a))
1800 21 : return sK > 0? 1: 3;
1801 : else
1802 77 : return sK > 0? 2: 1;
1803 : }
1804 189 : R = cgetg(4, t_VEC);
1805 756 : for (k = 1; k <= 3; k++) gel(R, k) = gel(nf_get_roots(gel(vnf,k)), 1);
1806 189 : return sK > 0? vecindexmin(R): vecindexmax(R);
1807 : }
1808 :
1809 : static GEN
1810 1337 : Z_factor_addprimes(GEN N, GEN P)
1811 : {
1812 : GEN Pnew, Enew;
1813 1337 : GEN H = Z_smoothen(N, P, &Pnew, &Enew);
1814 1337 : return H ? shallowconcat(P, gel(Z_factor(H),1)): P;
1815 : }
1816 :
1817 : static GEN
1818 889 : vbnf_discfactors(GEN vbnf)
1819 : {
1820 889 : long n = lg(vbnf)-1;
1821 889 : switch (n)
1822 : {
1823 441 : case 1:
1824 : {
1825 441 : GEN nf = bnf_get_nf(gel(vbnf,1));
1826 441 : GEN L = shallowtrans(nf_get_ramified_primes(nf));
1827 441 : return Z_factor_addprimes(nf_get_index(nf), L);
1828 : }
1829 259 : case 2:
1830 : {
1831 259 : GEN nf = gel(vbnf,2), R;
1832 259 : GEN L = shallowtrans(nf_get_ramified_primes(nf));
1833 259 : L = Z_factor_addprimes(nf_get_index(nf), L);
1834 259 : R = absi(ZX_resultant(nf_get_pol(gel(vbnf,1)), nf_get_pol(nf)));
1835 259 : return Z_factor_addprimes(R, L);
1836 : }
1837 189 : case 3:
1838 : {
1839 189 : GEN P1 = nf_get_pol(gel(vbnf,1));
1840 189 : GEN P2 = nf_get_pol(gel(vbnf,2));
1841 189 : GEN P3 = nf_get_pol(gel(vbnf,3));
1842 189 : GEN L = gel(Z_factor(absi(ZX_resultant(P1,P2))),1);
1843 189 : L = Z_factor_addprimes(absi(ZX_resultant(P2,P3)),L);
1844 189 : return Z_factor_addprimes(absi(ZX_resultant(P3,P1)),L);
1845 : }
1846 0 : default: pari_err_BUG("vbnf_discfactors");
1847 : return NULL;/*LCOV_EXCL_LINE*/
1848 : }
1849 : }
1850 :
1851 : static GEN
1852 889 : ell2selmer(GEN ell, GEN ell_K, GEN help, GEN K, GEN vbnf,
1853 : long effort, long flag, long prec)
1854 : {
1855 : GEN KP, pol, vnf, vpol, vcrt, sbase, LS2, factLS2, sqrtLS2, signs;
1856 : GEN selmer, helpLS2, LS2chars, helpchars, newselmer, factdisc, badprimes;
1857 : GEN helplist, listpoints, p, covers, disc, discF;
1858 889 : long i, k, n, tors2, mwrank, dim, nbpoints, lfactdisc, t, u, sha2 = 0;
1859 : forprime_t T;
1860 :
1861 889 : pol = ell2pol(ell);
1862 889 : help = ellsearchtrivialpoints(ell_K, flag ? NULL:muluu(LIMTRIV,effort+1), help);
1863 889 : help = elltwistpoints(help, ginv(K)); /* points on K Y^2 = pol(X) */
1864 889 : n = lg(vbnf) - 1; tors2 = n - 1;
1865 889 : KP = gel(absZ_factor(K), 1);
1866 889 : disc = ZX_disc(pol);
1867 889 : factdisc = mkvec3(KP, mkcol(gen_2), vbnf_discfactors(vbnf));
1868 889 : factdisc = ZV_sort_uniq_shallow(shallowconcat1(factdisc));
1869 889 : discF = mkvec2(gmul(K,disc), factdisc);
1870 889 : badprimes = cgetg(n+1, t_VEC);
1871 889 : vnf = cgetg(n+1, t_VEC);
1872 889 : vpol = cgetg(n+1, t_VEC);
1873 889 : vcrt = cgetg(n+1, t_VEC);
1874 889 : LS2 = cgetg(n+1, t_VEC);
1875 889 : factLS2 = cgetg(n+1, t_VEC);
1876 889 : sqrtLS2 = cgetg(n+1, t_VEC);
1877 2415 : for (k = 1; k <= n; k++)
1878 : {
1879 1526 : GEN T, Tinv, id, f, sel, L, S, nf, NF = gel(vbnf,k);
1880 : long i, lk;
1881 1526 : nf = (n == 1)? bnf_get_nf(NF): NF;
1882 1526 : gel(vnf, k) = nf;
1883 1526 : gel(vpol, k) = T = nf_get_pol(nf);
1884 1526 : Tinv = RgX_div(pol, gel(vpol, k));
1885 1526 : gel(vcrt, k) = QX_mul(QXQ_inv(T, Tinv), T); /* 0 mod T, 1 mod pol/T */
1886 :
1887 1526 : id = idealadd(nf, nf_get_index(nf), ZX_deriv(T));
1888 1526 : f = nf_pV_to_prV(nf, KP); settyp(f, t_COL);
1889 1526 : f = mkvec3(gel(idealfactor_partial(nf, Tinv, factdisc), 1),
1890 1526 : gel(idealfactor(nf, id), 1), f);
1891 1526 : gel(badprimes, k) = S = gtoset(shallowconcat1(f));
1892 1526 : if (n == 1)
1893 : {
1894 441 : sel = bnfselmer(NF, S);
1895 441 : obj_free(NF); /* units */
1896 : }
1897 1085 : else if (degpol(T) == 1)
1898 826 : sel = bnfselmerQ(S);
1899 : else /* degree 2 */
1900 259 : sel = nf2selmer_quad(NF, S);
1901 1526 : gel(LS2, k) = L = gel(sel, 1); lk = lg(L);
1902 1526 : gel(factLS2, k) = gel(sel, 2);
1903 1526 : gel(sqrtLS2, k) = gel(sel, 3);
1904 14952 : for (i = 1; i < lk; i++)
1905 : {
1906 13426 : GEN z = gel(L,i); /* z mod T, 1 mod (pol/T) */
1907 13426 : gel(L,i) = grem(gadd(z, gmul(gsubsg(1,z), gel(vcrt,k))), pol);
1908 : }
1909 : }
1910 889 : sbase = shallowconcat1(vecsmallbasis(vnf, vcrt, pol));
1911 889 : if (DEBUGLEVEL>2) err_printf(" local badprimes = %Ps\n", badprimes);
1912 889 : LS2 = shallowconcat1(LS2);
1913 889 : helpLS2 = ellLS2image(pol, help, K, vpol, vcrt);
1914 889 : LS2chars = matselmersign(vnf, vpol, LS2);
1915 889 : helpchars = matselmersign(vnf, vpol, helpLS2);
1916 889 : signs = NULL;
1917 889 : if (signe(ell_get_disc(ell)) > 0) signs = Flm_row(LS2chars, get_row(vnf,K));
1918 889 : selmer = kernorm(LS2, factdisc, pol, signs);
1919 889 : forprime_init(&T, gen_2, NULL); lfactdisc = lg(factdisc); dim = -1;
1920 7378 : for (i = 1; dim < 0 && i < lfactdisc; i++)
1921 6489 : dim = dim_selmer(gel(factdisc,i), pol, K, vnf, LS2, helpLS2,
1922 : &selmer,&LS2chars,&helpchars);
1923 2940 : while (dim < 0 && Flm_rank(Flm_mul(LS2chars, selmer, 2), 2) < lg(selmer)-1)
1924 : {
1925 2513 : while ((p = forprime_next(&T)) && ZV_search(factdisc, p));
1926 2051 : dim = dim_selmer(p, pol, K, vnf, LS2, helpLS2,
1927 : &selmer,&LS2chars,&helpchars);
1928 : }
1929 : /* transpose the matrix to get the solution that contains 1..tors2*/
1930 889 : helplist = gel(Flm_indexrank(Flm_transpose(helpchars),2), 1);
1931 889 : help = shallowextract(help, helplist);
1932 889 : helpchars = shallowextract(helpchars, helplist);
1933 889 : dim = lg(selmer)-1;
1934 889 : if (DEBUGLEVEL) err_printf("Selmer rank: %ld\n", dim);
1935 889 : newselmer = Flm_invimage(Flm_mul(LS2chars, selmer, 2), helpchars, 2);
1936 889 : if (!newselmer) pari_err_BUG("ellrank, wrong nfeltsign");
1937 889 : nbpoints = lg(help) - 1;
1938 889 : if (flag==1)
1939 : {
1940 49 : GEN u = nbpoints? Flm_mul(selmer,Flm_suppl(newselmer,2), 2): selmer;
1941 49 : long l = lg(u);
1942 49 : GEN z = cgetg(l, t_VEC);
1943 154 : for (i = 1; i < l; i++) gel(z,i) = RgXQV_factorback(LS2, gel(u,i), pol);
1944 49 : return mkvec2(mkvec3(vnf,sbase,pol), z);
1945 : }
1946 840 : else if (flag==2)
1947 : {
1948 56 : GEN u = nbpoints ? Flm_mul(selmer,Flm_suppl(newselmer,2), 2): selmer;
1949 56 : long l = lg(u);
1950 56 : GEN P = cgetg(l, t_VEC), b;
1951 147 : for (i = 1; i < l; i++)
1952 : {
1953 91 : b = liftselmerinit(gel(u,i), vnf, sqrtLS2, factLS2, badprimes, vcrt, pol);
1954 91 : gel(P,i) = liftselmer_cover(b, gel(u,i), LS2, pol, discF, K);
1955 : }
1956 56 : return P;
1957 : }
1958 784 : listpoints = vec_lengthen(help, dim); /* points on ell */
1959 784 : covers = zerovec(dim);
1960 5789 : for (i=1; i <= dim; i++)
1961 : {
1962 5005 : GEN b, P, expo, vec = vecsmall_ei(dim, i);
1963 5005 : if (Flm_Flc_invimage(newselmer, vec, 2)) continue;
1964 2555 : expo = Flm_Flc_mul(selmer, vec, 2);
1965 2555 : b = liftselmerinit(expo, vnf, sqrtLS2, factLS2, badprimes, vcrt, pol);
1966 2555 : P = liftselmer(b, expo, sbase, LS2, pol, discF, K, 1, &gel(covers,i));
1967 2555 : if (P)
1968 : {
1969 1456 : gel(listpoints, ++nbpoints) = P; /* new point on ell */
1970 1456 : gel(newselmer, nbpoints) = vec;
1971 1456 : setlg(newselmer, nbpoints+1);
1972 : }
1973 : }
1974 784 : if (nbpoints < dim)
1975 : {
1976 : long i, j;
1977 315 : GEN M = cgetg(dim+1, t_MAT), selker;
1978 315 : GEN D = mulii(muliu(absi(disc), 27*4096), powiu(K,6));
1979 315 : GEN FD = ZV_sort_uniq_shallow(shallowconcat1(mkvec2(mkcol3s(3,5,7), factdisc)));
1980 :
1981 1974 : for (i = 1; i <= dim; i++) gel(M,i) = cgetg(dim+1, t_COL);
1982 1974 : for (i = 1; i <= dim; i++)
1983 9415 : for (j = 1; j <= i; j++)
1984 : {
1985 : GEN Q;
1986 7756 : if (isintzero(gel(covers,i)))
1987 3150 : Q = gen_0;
1988 4606 : else if (i==j)
1989 1057 : Q = quartic_findunit(D, gel(covers,i));
1990 : else
1991 : {
1992 3549 : GEN e = Flv_add(gel(selmer,i), gel(selmer,j), 2);
1993 3549 : GEN b = liftselmerinit(e, vnf, sqrtLS2, factLS2, badprimes, vcrt, pol);
1994 3549 : Q = quartic_findunit(D, gel(liftselmer_cover(b, e, LS2, pol, discF, K),1));
1995 : }
1996 7756 : gmael(M,j,i) = gmael(M,i,j) = Q;
1997 : }
1998 315 : selker = F2m_to_Flm(F2m_ker(matcassels(FD, M)));
1999 315 : sha2 = dim - (lg(selker)-1);
2000 315 : dim = lg(selker)-1;
2001 1133 : for (t=1, u=1; nbpoints < dim && effort > 0; t++)
2002 : {
2003 818 : pari_sp btop = avma;
2004 : GEN expo, b, P, vec;
2005 1091 : do vec = Flm_Flc_mul(selker,random_Flv(dim, 2), 2);
2006 1091 : while (zv_equal0(vec) || Flm_Flc_invimage(newselmer, vec, 2));
2007 818 : expo = Flm_Flc_mul(selmer, vec, 2);
2008 818 : b = liftselmerinit(expo, vnf, sqrtLS2, factLS2, badprimes, vcrt, pol);
2009 818 : P = liftselmer(b, expo, sbase, LS2, pol, discF, K, u, NULL);
2010 818 : if (P)
2011 : {
2012 329 : gel(listpoints, ++nbpoints) = P;
2013 329 : gel(newselmer, nbpoints) = vec;
2014 329 : setlg(newselmer, nbpoints+1);
2015 489 : } else set_avma(btop);
2016 818 : if (t == dim) { t = 0; u++; effort--; }
2017 : }
2018 : }
2019 784 : setlg(listpoints, nbpoints+1);
2020 784 : mwrank = nbpoints - tors2;
2021 784 : if (odd(dim - nbpoints)) mwrank++;
2022 784 : listpoints = vecslice(listpoints,tors2+1, nbpoints);
2023 784 : listpoints = elltwistpoints(listpoints, K);
2024 784 : listpoints = vecellabs(ellQ_genreduce(ell_K, listpoints, NULL, prec));
2025 784 : gtoset_inplace(listpoints);
2026 784 : return mkvec4(utoi(mwrank), utoi(dim-tors2), utoi(sha2), listpoints);
2027 : }
2028 :
2029 : GEN
2030 49 : ell2selmer_basis(GEN ell, GEN *cb, long prec)
2031 : {
2032 49 : GEN E = ellminimalbmodel(ell, cb);
2033 49 : GEN S = ell2selmer(E, E, NULL, gen_1, makevbnf(E, prec), 0, 1, prec);
2034 49 : obj_free(E); return S;
2035 : }
2036 :
2037 : static void
2038 98 : ellrank_get_nudur(GEN E, GEN F, GEN *nu, GEN *du, GEN *r)
2039 : {
2040 98 : GEN ea2 = ell_get_a2(E), ea2t = ell_get_a2(F);
2041 98 : GEN ec4 = ell_get_c4(E), ec4t = ell_get_c4(F);
2042 98 : GEN ec6 = ell_get_c6(E), ec6t = ell_get_c6(F);
2043 : GEN N, D, d;
2044 98 : if (signe(ec4)==0)
2045 : {
2046 21 : if (!Z_ispowerall(mulii(ec6,sqri(ec6t)),3,&N))
2047 7 : pari_err_TYPE("ellrank",F);
2048 14 : D = ec6t;
2049 : }
2050 77 : else if (signe(ec6)==0)
2051 : {
2052 21 : if (!Z_issquareall(mulii(ec4,ec4t),&N))
2053 7 : pari_err_TYPE("ellrank",F);
2054 14 : D = ec4t;
2055 : }
2056 : else
2057 : {
2058 56 : GEN d46 = mulii(gcdii(ec4,ec4t),gcdii(ec6,ec6t));
2059 56 : N = diviiexact(mulii(ec6,ec4t),d46);
2060 56 : D = diviiexact(mulii(ec6t,ec4),d46);
2061 : }
2062 84 : d = gcdii(N, D);
2063 84 : *nu = diviiexact(N, d); *du = diviiexact(D, d);
2064 84 : *r = diviuexact(subii(mulii(*nu,ea2t),mulii(*du,ea2)),3);
2065 84 : }
2066 :
2067 : static GEN
2068 882 : ellrank_flag(GEN e, long effort, GEN help, long flag, long prec)
2069 : {
2070 882 : pari_sp ltop = avma;
2071 : GEN urst, v, vbnf, eK;
2072 882 : GEN et = NULL, K = gen_1, nu = NULL, du = NULL, r = NULL, urstK = NULL;
2073 882 : long newell = 0;
2074 :
2075 882 : if (lg(e)==3 && typ(e)==t_VEC) { et = gel(e,2); e = gel(e,1); }
2076 882 : if (lg(e)==4 && typ(e)==t_VEC)
2077 : {
2078 126 : vbnf = gel(e,3); urst = gel(e,2);
2079 126 : e = gel(e,1); checkell_Q(e);
2080 126 : if (!ell_is_integral(e)) pari_err_TYPE("ellrank [nonintegral model]",e);
2081 119 : if (signe(ell_get_a1(e))) pari_err_TYPE("ellrank [a1 != 0]", e);
2082 112 : if (signe(ell_get_a3(e))) pari_err_TYPE("ell2rank [a3 != 0]", e);
2083 : }
2084 : else
2085 : {
2086 756 : checkell_Q(e);
2087 756 : e = ellminimalbmodel(e, &urst);
2088 756 : newell = 1;
2089 756 : vbnf = makevbnf(e, prec);
2090 : }
2091 861 : if (et)
2092 : {
2093 105 : checkell_Q(et);
2094 105 : if (!gequal(ell_get_j(et),ell_get_j(e))) pari_err_TYPE("ellrank",et);
2095 98 : et = ellminimalbmodel(et, &urst);
2096 98 : ellrank_get_nudur(e, et, &nu, &du, &r);
2097 84 : K = mulii(nu, du);
2098 84 : urstK = mkvec4(nu, mulii(nu,r), gen_0, gen_0);
2099 : }
2100 840 : if (help)
2101 : {
2102 42 : if (urst) help = ellchangepoint(help, urst);
2103 42 : if (et) help = ellchangepointinv(help, urstK);
2104 : }
2105 840 : eK = elltwistequation(e, K);
2106 : /* help is a vector of rational points [x,y] satisfying K y^2 = pol(x) */
2107 : /* [Kx, K^2y] is on eK: Y^2 = K^3 pol(X/K) */
2108 840 : if (help) check_oncurve(eK, help);
2109 840 : v = ell2selmer(e, eK, help, K, vbnf, effort, flag, prec);
2110 840 : if (flag==0)
2111 : {
2112 784 : if (et) gel(v,4) = ellchangepoint(gel(v,4), urstK);
2113 784 : if (urst) gel(v,4) = ellchangepointinv(gel(v,4), urst);
2114 : }
2115 : else
2116 : {
2117 56 : long i, l = lg(v);
2118 147 : for (i = 1; i < l; i++)
2119 : {
2120 91 : if (et) gmael(v,i,2) = ellchangepoint(gmael(v,i,2), urstK);
2121 91 : if (urst) gmael(v,i,2) = ellchangepointinv(gmael(v,i,2), urst);
2122 : }
2123 : }
2124 840 : if (newell) obj_free(e);
2125 840 : if (eK != e) obj_free(eK);
2126 840 : return gc_GEN(ltop, v);
2127 : }
2128 :
2129 : GEN
2130 826 : ellrank(GEN e, long effort, GEN help, long prec)
2131 : {
2132 826 : return ellrank_flag(e, effort, help, 0, prec);
2133 : }
2134 :
2135 : GEN
2136 56 : ell2cover(GEN ell, long prec)
2137 : {
2138 56 : return ellrank_flag(ell, 0, NULL, 2, prec);
2139 : }
|
