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