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 162652 : to_principal_unit(GEN nf, GEN x, GEN pr, GEN sprk)
158 : {
159 162652 : 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 162652 : 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 115870 : FpX_issquare(GEN q, GEN p)
453 : {
454 115870 : GEN F = FpX_factor_squarefree(q,p);
455 115870 : long i, l = lg(F);
456 154762 : for (i = 1; i < l; i+=2)
457 115583 : if (degpol(gel(F,i)) > 0) return 0;
458 39179 : return 1;
459 : }
460 :
461 : static GEN
462 115422 : hyperell_red(GEN q, GEN p)
463 : {
464 : GEN Q;
465 115422 : long v = ZX_pvalrem(q, p, &Q);
466 115422 : if (v == 1) return q;
467 115422 : return odd(v)? ZX_Z_mul(Q, p): Q;
468 : }
469 :
470 : static GEN
471 117249 : hyperell_reg_point(GEN q, GEN p)
472 : {
473 : GEN Q, F;
474 117249 : long i, l, v = ZX_pvalrem(q, p, &Q);
475 117249 : if (v != 1) q = odd(v)? ZX_Z_mul(Q, p): Q;
476 117249 : if (!odd(v))
477 : {
478 115870 : GEN qr = FpX_red(q, p);
479 115870 : if (!FpX_issquare(qr,p) || Fp_issquare(leading_coeff(qr), p))
480 115422 : return mkvec2s(0,1);
481 : }
482 1827 : F = FpX_roots(Q, p); l = lg(F);
483 1834 : for (i = 1; i < l; i++)
484 : {
485 1827 : GEN r = gel(F,i), s = hyperell_reg_point(ZX_affine(q, p, r), p);
486 1827 : if (s)
487 1820 : retmkvec2(addii(r,mulii(gel(s,1),p)), mulii(gel(s,2),p));
488 : }
489 7 : return NULL;
490 : }
491 :
492 : static GEN
493 1331 : hyperell_lift(GEN q, GEN x, GEN p)
494 : {
495 : long e;
496 1331 : GEN qy2 = ZX_Z_sub(q, sqri(p));
497 1331 : for (e = 2; ; e<<=1)
498 785 : {
499 2116 : pari_sp av = avma;
500 2116 : GEN z = ZpX_liftroot(qy2, x, p, e);
501 2116 : if (signe(z) == 0) z = powiu(p, e);
502 2116 : if (Zp_issquare(poleval(q, z), p)) return z;
503 785 : set_avma(av);
504 : }
505 : }
506 :
507 : static GEN
508 57711 : affine_apply(GEN r, GEN x)
509 : {
510 57711 : return addii(mulii(gel(r,2),x), gel(r,1));
511 : }
512 :
513 : static GEN
514 57711 : Qp_hyperell_solve_odd(GEN q, GEN p)
515 : {
516 57711 : GEN qi = RgX_recip_shallow(q);
517 57711 : GEN r = hyperell_reg_point(q, p), qr = NULL, qrp = NULL;
518 57711 : GEN s = hyperell_reg_point(qi, p), qs = NULL, qsp = NULL;
519 57711 : if (!r && !s) return NULL;
520 57711 : if (r)
521 : {
522 57711 : qr = hyperell_red(ZX_affine(q, gel(r,2), gel(r,1)), p);
523 57711 : qrp = FpX_deriv(qr, p);
524 : }
525 57711 : if (s)
526 : {
527 57711 : qs = hyperell_red(ZX_affine(qi, gel(s,2), gel(s,1)), p);
528 57711 : qsp = FpX_deriv(qs, p);
529 : }
530 : while(1)
531 15011 : {
532 72722 : GEN x = randomi(p);
533 72722 : if (r)
534 : {
535 72722 : GEN y2 = FpX_eval(qr, x, p);
536 72722 : if (Fp_issquare(y2,p))
537 : {
538 47605 : if (signe(y2))
539 43325 : return affine_apply(r,x);
540 4280 : if (signe(FpX_eval(qrp, x, p)))
541 : {
542 1016 : x = hyperell_lift(qr, x, p);
543 1016 : return affine_apply(r,x);
544 : }
545 : }
546 : }
547 28381 : if (s)
548 : {
549 28381 : GEN y2 = FpX_eval(qs, x, p);
550 28381 : if (Fp_issquare(y2,p))
551 : {
552 15187 : if (signe(x)==0) x = p;
553 15187 : if (signe(y2))
554 13055 : return ginv(affine_apply(s,x));
555 2132 : if (signe(FpX_eval(qsp, x, p)))
556 : {
557 315 : GEN xl = hyperell_lift(qs, x, p);
558 315 : return ginv(affine_apply(s,xl));
559 : }
560 : }
561 : }
562 : }
563 : }
564 :
565 : static GEN
566 469 : Q2_hyperell_lift(GEN p, GEN q, long x, long y)
567 : {
568 : GEN T, h;
569 : long e;
570 469 : if (y==0) y = 2;
571 469 : T = ZX_sub(p, ZX_Z_add(ZX_mulu(q, y), sqru(y)));
572 469 : h = ZX_add(ZX_sqr(q), ZX_shifti(p, 2));
573 469 : for (e = 3; ; e++)
574 840 : {
575 1309 : pari_sp av = avma;
576 1309 : GEN r = ZpX_liftroot(T, utoi(x), gen_2, e);
577 1309 : if (signe(r) == 0) r = int2n(e);
578 1309 : if (Zp_issquare(poleval(h, r), gen_2)) return r;
579 840 : set_avma(av);
580 : }
581 : return NULL;/*LCOV_EXCL_LINE*/
582 : }
583 :
584 : static GEN
585 2827 : Q2_hyperell_regpoint(GEN P, GEN Q)
586 : {
587 : long x;
588 2827 : GEN p = ZX_to_F2x(P), dp = F2x_deriv(p);
589 2827 : GEN q = ZX_to_F2x(Q), dq = F2x_deriv(q);
590 :
591 4659 : for (x = 0; x <= 1; x++)
592 : {
593 4128 : long px = F2x_eval(p,x), qx = F2x_eval(q,x);
594 : long dpx, dqx;
595 4128 : if (qx == 1)
596 : {
597 1981 : if(px == 0) return x==0 ? gen_2: gen_1;
598 154 : continue;
599 : }
600 2147 : dpx = F2x_eval(dp,x);
601 2147 : dqx = F2x_eval(dq,x);
602 2147 : if (odd(dqx * px + dpx))
603 469 : return Q2_hyperell_lift(P, Q, x, px);
604 : }
605 531 : return NULL;
606 : }
607 :
608 : static GEN
609 2827 : Q2_hyperell_solve_affine(GEN p, GEN q)
610 : {
611 2827 : pari_sp av = avma;
612 : GEN R, p4, q4;
613 : long x;
614 : while(1)
615 2317 : {
616 : GEN pp, p0, p1;
617 5144 : long vp = ZX_lval(p, 2);
618 5144 : long vq = ZX_lval(q, 2);
619 5144 : long w = minss(vp>>1, vq);
620 5144 : p = ZX_shifti(p, -2*w);
621 5144 : q = ZX_shifti(q, -w);
622 5144 : if (ZX_lval(q,2)==0) break;
623 2590 : pp = RgX_splitting(p,2); p0 = gel(pp,1); p1 = gel(pp,2);
624 2590 : if (ZX_lval(p1,2)==0 || ZX_lval(p0,2)>=1) break;
625 2317 : p = ZX_sub(p, ZX_mul(p0, ZX_add(q, p0)));
626 2317 : q = ZX_add(q, ZX_shifti(p0, 1));
627 : }
628 2827 : R = Q2_hyperell_regpoint(p, q);
629 2827 : if (R) return gerepileuptoint(av, R);
630 531 : p4 = ZX_to_Flx(p,4);
631 531 : q4 = ZX_to_Flx(q,4);
632 678 : for (x = 0; x <= 1; x++)
633 : {
634 671 : ulong px = Flx_eval(p4, x, 4);
635 671 : ulong qx = Flx_eval(q4, x, 4);
636 671 : if (px == 0 || (1+qx+3*px)%4==0)
637 : {
638 524 : GEN p2 = ZX_affine(p, gen_2, utoi(x));
639 524 : GEN q2 = ZX_affine(q, gen_2, utoi(x));
640 524 : GEN S = Q2_hyperell_solve_affine(p2, q2);
641 524 : 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 60007 : hyperell_local_solve(GEN q, GEN p)
658 : {
659 60007 : if (equaliu(p,2))
660 2296 : return Q2_hyperell_solve(q);
661 57711 : return Qp_hyperell_solve_odd(q, p);
662 : }
663 :
664 : /*******************************************************************/
665 : /* BINARY QUARTIC */
666 : /*******************************************************************/
667 : static int
668 104825 : Qp_issquare(GEN a, GEN p)
669 : {
670 104825 : GEN b = typ(a) == t_INT? a: mulii(gel(a,1), gel(a,2));
671 104825 : return Zp_issquare(b, p);
672 : }
673 :
674 : static GEN
675 18403 : quartic_IJ(GEN g)
676 : {
677 18403 : GEN a = gel(g, 6), b = gel(g, 5), c = gel(g, 4), d = gel(g, 3), e = gel(g, 2);
678 18403 : GEN ae = gmul(a,e), bd = gmul(b,d), c2 = gsqr(c);
679 : /* 12ae - 3bd + c^2 */
680 18403 : GEN I = gadd(gsub(gmulsg(12, ae), gmulsg(3, bd)), c2);
681 : /* c(72ae + 9bd - 2c^2) - 27ad^2 - 27eb^2*/
682 18403 : 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 18403 : 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 12761 : quartic_cubic(GEN g, long v)
699 : {
700 12761 : GEN a = gel(g, 6), b = gel(g, 5), c = gel(g, 4);
701 12761 : GEN a3 = gdivgu(a,3);
702 12761 : return deg1pol(gmul2n(a3,2), gsub(gsqr(b),gmul2n(gmul(a3,c),3)), v);
703 : }
704 :
705 : static GEN
706 6804 : quarticinv_pol(GEN IJ)
707 : {
708 6804 : GEN I = gel(IJ,1), J = gel(IJ,2);
709 6804 : 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 4508 : quartic_disc(GEN q)
724 : {
725 4508 : GEN IJ = quartic_IJ(q), I = gel(IJ,1), J = gel(IJ,2);
726 4508 : return gsub(gmul2n(gpowgs(I, 3), 2), gsqr(J));
727 : }
728 :
729 : static GEN
730 7091 : quartic_minim_all(GEN F, GEN discF)
731 : {
732 7091 : GEN IJ = quartic_IJ(F), I = gel(IJ,1), J = gel(IJ,2);
733 7091 : GEN g = Z_ppo(gcdii(I,J), gel(discF,1));
734 7091 : GEN plist = ZV_sort_uniq_shallow(shallowconcat(gel(absZ_factor(g),1),gel(discF,2)));
735 7091 : GEN W, C, PQ = hyperellminimalmodel(F, &C, plist);
736 7091 : GEN P = gel(PQ,1), Q = gel(PQ,2);
737 7091 : if (signe(Q)==0)
738 812 : W = mkvec2(P, C);
739 : else
740 6279 : W = mkvec2(ZX_add(ZX_shifti(P,2),ZX_sqr(Q)),mkvec2(shifti(gel(C,1),-1),gel(C,2)));
741 7091 : return W;
742 : }
743 :
744 : /*******************************************************************/
745 : /* Cassels' pairing */
746 : /*******************************************************************/
747 :
748 : static GEN
749 6097 : nfsqrt(GEN nf, GEN z)
750 : {
751 6097 : long v = fetch_var_higher();
752 6097 : GEN R = nfroots(nf, deg2pol_shallow(gen_m1, gen_0, z, v));
753 6097 : delete_var();
754 6097 : return lg(R)==1 ? NULL: gel(R,1);
755 : }
756 :
757 : static GEN
758 6097 : nfsqrt_safe(GEN nf, GEN z)
759 : {
760 6097 : GEN r = nfsqrt(nf, z);
761 6097 : if (!r) pari_err_BUG("ellrank");
762 6097 : return r;
763 : }
764 :
765 : static GEN
766 2296 : vecnfsqrtmod(GEN x, GEN P)
767 8393 : { 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 525 : cassels_oo_solve_i(GEN q, GEN g)
780 : {
781 525 : long dg = degpol(g);
782 : GEN D, a, b, c;
783 :
784 525 : 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 525 : cassels_oo_solve(GEN q, GEN g)
802 525 : { pari_sp av = avma; return gc_int(av, cassels_oo_solve_i(q, g)); }
803 :
804 : static GEN
805 60007 : cassels_Qp_solve(GEN q, GEN gam, GEN p)
806 : {
807 60007 : pari_sp av = avma;
808 60007 : GEN a = hyperell_local_solve(q, p);
809 60007 : GEN c = poleval(gam,a);
810 : long e;
811 60007 : if (!gequal0(c)) return c;
812 70 : for (e = 2; ; e++)
813 0 : {
814 70 : GEN b = gadd(a, powiu(p,e));
815 70 : if (Qp_issquare(poleval(q, b), p))
816 : {
817 70 : c = poleval(gam, b);
818 70 : if (!gequal0(c)) return gerepileupto(av,c);
819 : }
820 : }
821 : }
822 :
823 : static GEN
824 4592 : to_ZX(GEN a, long v) { return typ(a)==t_INT? scalarpol_shallow(a,v): a; }
825 :
826 : static GEN
827 4508 : quartic_findunit(GEN D, GEN q)
828 : {
829 4508 : GEN T = quarticinv_pol(quartic_IJ(q));
830 : while(1)
831 1365 : {
832 5873 : pari_sp av = avma;
833 5873 : GEN z = quartic_cubic(q,0);
834 5873 : if (signe(QXQ_norm(z,T)))
835 4508 : 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 525 : if (signe(leading_coeff(gam)) < 0) gam = ZX_neg(gam);
863 525 : if (cassels_oo_solve(q1, gam) < 0) e = 1;
864 : }
865 62303 : for (i = 1; i < lF; i++)
866 : {
867 60007 : GEN p = gel(F, i);
868 60007 : GEN c = cassels_Qp_solve(q1, gam, p);
869 60007 : if (hilbert(c, a, p) < 0) e = !e;
870 : }
871 2296 : return gc_long(av,e);
872 : }
873 :
874 : static GEN
875 301 : matcassels(GEN F, GEN M)
876 : {
877 301 : long i, j, n = lg(M)-1;
878 301 : GEN C = zero_F2m_copy(n,n);
879 301 : pari_sp av = avma;
880 1911 : for (i = 1; i <= n; i++)
881 : {
882 1610 : GEN Mii = gcoeff(M,i,i);
883 1610 : if (isintzero(Mii)) continue;
884 4508 : for (j = 1; j < i; j++)
885 : {
886 3500 : GEN Mjj = gcoeff(M,j,j);
887 3500 : if (!isintzero(Mjj) && casselspairing(F, Mii, Mjj, gcoeff(M,i,j)))
888 336 : { F2m_set(C,i,j); F2m_set(C,j,i); }
889 : }
890 : }
891 301 : 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 1701 : ell2pol(GEN ell)
907 1701 : { 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 4353 : projratpointxz(GEN pol, long lim, GEN *py)
912 : {
913 : pari_timer ti;
914 : GEN P;
915 4353 : if (issquareall(leading_coeff(pol), py)) return mkcol2(gen_1, gen_0);
916 4304 : if (DEBUGLEVEL) timer_start(&ti);
917 4304 : P = hyperellratpoints(pol, stoi(lim), 1);
918 4304 : if (DEBUGLEVEL) timer_printf(&ti,"hyperellratpoints(%ld)",lim);
919 4304 : if (lg(P) == 1) return NULL;
920 1735 : P = gel(P,1); *py = gel(P,2); return mkcol2(gel(P,1), gen_1);
921 : }
922 :
923 : /* P a list of integers (actually primes) one of which divides x; return
924 : * the first one */
925 : static GEN
926 623 : first_divisor(GEN x, GEN P)
927 : {
928 623 : long i, n = lg(P);
929 623 : for (i = 1; i < n; i++)
930 623 : if (dvdii(x, gel(P,i))) return gel(P,i);
931 0 : return gel(P,i);
932 : }
933 :
934 : /* find point (x:y:z) on y^2 = pol, return [x,z]~ and set *py = y */
935 : static GEN
936 1694 : projratpointxz2(GEN pol, long lim, GEN *py)
937 : {
938 1694 : pari_sp av = avma;
939 1694 : GEN list = mkvec(mkvec4(pol, matid(2), gen_1, gen_1));
940 : long i, j, c;
941 :
942 4501 : for (i = 1, c = 1; i < lg(list); i++,c++)
943 : {
944 2842 : GEN K, k, ff, co, p, M, C, r, pol, L = gel(list, i);
945 : long lr;
946 :
947 2842 : list = vecsplice(list, i); i--;
948 2842 : pol = Q_primitive_part(gel(L,1), &K);
949 2842 : M = gel(L,2);
950 2842 : K = K? mulii(gel(L,3), K): gel(L,3);
951 2842 : C = gel(L,4);
952 2842 : if (Z_issquareall(K, &k))
953 : {
954 : GEN xz, y, aux, U;
955 3094 : if (c==1) continue;
956 910 : pol = ZX_hyperellred(pol, &U);
957 910 : if (DEBUGLEVEL>1) err_printf(" reduced quartic(%ld): Y^2 = %Ps\n", i, pol);
958 910 : xz = projratpointxz(pol, lim, &y); if (!xz) continue;
959 35 : *py = gmul(y, mulii(C, k));
960 35 : aux = RgM_RgC_mul(ZM2_mul(M, U), xz);
961 70 : if (gequal0(gel(aux, 2))) return mkcol2(gel(aux,1), gen_0);
962 35 : *py = gdiv(*py, gpowgs(gel(aux, 2), degpol(pol)>>1));
963 35 : return mkcol2(gdiv(gel(aux, 1), gel(aux, 2)), gen_1);
964 : }
965 623 : ff = Z_factor(K); co = core2(mkvec2(K, ff)); K = gel(co,1); /* > 1 */
966 623 : p = first_divisor(K, gel(ff,1));
967 623 : K = diviiexact(K, p);
968 623 : C = mulii(mulii(C, gel(co,2)), p);
969 : /* root at infinity */
970 623 : if (dvdii(leading_coeff(pol), p))
971 : {
972 329 : GEN U = mkmat2(gel(M,1), ZC_Z_mul(gel(M,2), p));
973 329 : if (equali1(content(U)))
974 : {
975 147 : GEN t = ZX_rescale(pol, p);
976 147 : list = vec_append(list, mkvec4(ZX_Z_divexact(t,p), U, K, C));
977 : }
978 : }
979 623 : r = FpC_center(FpX_roots(pol, p), p, shifti(p,-1)); lr = lg(r);
980 1701 : for (j = 1; j < lr; j++)
981 : {
982 1078 : GEN U = ZM2_mul(M, mkmat22(p, gel(r, j), gen_0, gen_1));
983 1078 : if (equali1(content(U)))
984 : {
985 1050 : GEN t = ZX_unscale_div(ZX_translate(pol, gel(r,j)), p);
986 1050 : list = vec_append(list, mkvec4(t, U, K, C));
987 : }
988 : }
989 623 : if (gc_needed(av, 1)) gerepileall(av, 2, &pol, &list);
990 : }
991 1659 : return NULL;
992 : }
993 :
994 : static GEN
995 8316 : polrootsmodpn(GEN pol, GEN p)
996 : {
997 8316 : pari_sp av = avma, av2;
998 8316 : long j, l, i = 1, vd = Z_pval(ZX_disc(pol), p);
999 : GEN v, r, P;
1000 :
1001 8316 : if (!vd) { set_avma(av); retmkvec(zerovec(2)); }
1002 6090 : pol = Q_primpart(pol);
1003 6090 : P = gpowers0(p, vd-1, p); av2 = avma;
1004 6090 : v = FpX_roots(pol, p); l = lg(v);
1005 17752 : for (j = 1; j < l; j++) gel(v,j) = mkvec2(gel(v,j), gen_1);
1006 69923 : while (i < lg(v))
1007 : {
1008 63833 : GEN pol2, pe, roe = gel(v, i), ro = gel(roe,1);
1009 63833 : long e = itou(gel(roe,2));
1010 :
1011 64820 : if (e >= vd) { i++; continue; }
1012 49070 : pe = gel(P, e);
1013 49070 : (void)ZX_pvalrem(ZX_affine(pol, pe, ro), p, &pol2);
1014 49070 : r = FpX_roots(pol2, p); l = lg(r);
1015 49070 : if (l == 1) { i++; continue; }
1016 100254 : for (j = 1; j < l; j++)
1017 52171 : 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 48083 : if (l > 2) v = shallowconcat(v, vecslice(r, 2, l-1));
1020 48083 : gel(v, i) = gel(r, 1);
1021 48083 : if (gc_needed(av2, 1)) gerepileall(av2, 1, &v);
1022 : }
1023 6090 : if (lg(v) == 1) { set_avma(av); retmkvec(zerovec(2)); }
1024 6090 : 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 1596945 : kpmodsquares1(GEN nf, GEN z, GEN sprk)
1035 : {
1036 1596945 : GEN modpr = (typ(sprk) == t_VEC)? gmael(sprk, 4, 1): sprk;
1037 1596945 : GEN pr = modpr_get_pr(modpr), p = pr_get_p(pr);
1038 1596945 : long v = nfvalrem(nf, z, pr, &z);
1039 1596945 : if (equaliu(p, 2))
1040 : {
1041 : GEN c;
1042 156744 : z = to_principal_unit(nf, z, pr, sprk); /* = 1 mod pr */
1043 156744 : c = ZV_to_Flv(sprk_log_prk1(nf, z, sprk), 2);
1044 156744 : return vecsmall_prepend(c, odd(v));
1045 : }
1046 : else
1047 : {
1048 1440201 : GEN T = modpr_get_T(modpr);
1049 1440201 : long c = !Fq_issquare(nf_to_Fq(nf, z, modpr), T, p);
1050 1440201 : return mkvecsmall2(odd(v), c);
1051 : }
1052 : }
1053 :
1054 : static GEN
1055 782719 : kpmodsquares(GEN vnf, GEN z, GEN PP)
1056 : {
1057 782719 : pari_sp av = avma;
1058 782719 : long i, j, l = lg(vnf);
1059 782719 : GEN dz, vchar = cgetg(l, t_VEC);
1060 782719 : z = Q_remove_denom(z, &dz); if (dz) z = ZX_Z_mul(z, dz);
1061 1859564 : for (i = 1; i < l; i++)
1062 : {
1063 1076845 : GEN nf = gel(vnf, i), pp = gel(PP, i);
1064 1076845 : GEN kp, delta = ZX_rem(z, nf_get_pol(nf));
1065 1076845 : long lp = lg(pp);
1066 1076845 : kp = cgetg(lp, t_VEC);
1067 2673790 : for (j = 1; j < lp; j++) gel(kp, j) = kpmodsquares1(nf, delta, gel(pp,j));
1068 1076845 : gel(vchar, i) = shallowconcat1(kp);
1069 : }
1070 782719 : return gerepileuptoleaf(av, shallowconcat1(vchar));
1071 : }
1072 :
1073 : static GEN
1074 16632 : veckpmodsquares(GEN x, GEN vnf, GEN PP)
1075 777707 : { 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 3898 : ellabs(GEN P)
1083 3898 : { return ell_is_inf(P) ? P: mkvec2(gel(P,1), Q_abs_shallow(gel(P,2))); }
1084 : static GEN
1085 4654 : 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 812 : elltwistequation(GEN ell, GEN K)
1090 : {
1091 : GEN K2, a2, a4, a6;
1092 812 : 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 1617 : elltwistpoints(GEN x, GEN K)
1110 : {
1111 : GEN K2;
1112 1617 : if (!K || gequal1(K)) return x;
1113 126 : K2 = gsqr(K);
1114 350 : pari_APPLY_same(elltwistpoint(gel(x,i), K, K2))
1115 : }
1116 :
1117 : /* \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */
1118 : /* FUNCTIONS FOR NUMBER FIELDS \\ */
1119 : /* \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */
1120 :
1121 : /* return a set S2 of prime ideals disjoint from S such that
1122 : * Cl_{S \cup S2}(K) has no p-torsion */
1123 : static GEN
1124 434 : bestS(GEN bnf,GEN S, ulong p)
1125 : {
1126 434 : GEN v, S2, h = bnf_get_no(bnf), cyc = bnf_get_cyc(bnf);
1127 : long i, lS2;
1128 : ulong l, vD;
1129 : forprime_t P;
1130 :
1131 434 : if (!dvdiu(h, p)) return cgetg(1,t_VEC);
1132 315 : if (!S)
1133 : {
1134 0 : v = diagonal_shallow(cyc);
1135 0 : vD = Z_lval(h, p);
1136 : }
1137 : else
1138 : {
1139 315 : long lS = lg(S);
1140 315 : v = cgetg(lS,t_MAT);
1141 994 : for (i = 1; i < lS; i++) gel(v,i) = isprincipal(bnf, gel(S,i));
1142 315 : v = ZM_hnfmodid(v, cyc);
1143 315 : vD = Z_lval(ZM_det(v), p); if (!vD) return cgetg(1, t_VEC);
1144 : }
1145 301 : S2 = cgetg(vD+2, t_VEC); lS2 = 1;
1146 301 : u_forprime_init(&P,2,ULONG_MAX);
1147 2569 : while ((l = u_forprime_next(&P)))
1148 : {
1149 2569 : pari_sp av = avma;
1150 2569 : GEN w, Sl = idealprimedec(bnf, utoi(l));
1151 2569 : long nSl = lg(Sl)-1;
1152 : ulong vDl;
1153 4074 : for (i = 1; i < nSl; i++) /* remove one prime ideal */
1154 : {
1155 1806 : w = ZM_hnf(shallowconcat(v, isprincipal(bnf, gel(Sl,i))));
1156 1806 : vDl = Z_lval(ZM_det(w), p);
1157 1806 : if (vDl < vD)
1158 : {
1159 1239 : gel(S2,lS2++) = gel(Sl,i);
1160 1239 : vD = vDl; v = w; av = avma;
1161 1239 : if (!vD) { setlg(S2, lS2); return S2; }
1162 : }
1163 : }
1164 2268 : set_avma(av);
1165 : }
1166 : return NULL;/*LCOV_EXCL_LINE*/
1167 : }
1168 :
1169 : static GEN
1170 868 : nfC_prV_val(GEN nf, GEN G, GEN P)
1171 : {
1172 868 : long i, j, lG = lg(G), lP = lg(P);
1173 868 : GEN M = cgetg(lG, t_MAT);
1174 5362 : for (i = 1; i < lG; i++)
1175 : {
1176 4494 : GEN V = cgetg(lP, t_COL);
1177 19285 : for (j = 1; j < lP; j++)
1178 14791 : gel(V,j) = gpnfvalrem(nf, gel(G,i), gel(P,j), NULL);
1179 4494 : gel(M,i) = V;
1180 : }
1181 868 : return M;
1182 : }
1183 :
1184 : static GEN
1185 5558 : _factorbackmod(GEN nf, GEN g, GEN e, ulong p)
1186 : {
1187 5558 : GEN y = nffactorback(nf, g, e), den;
1188 5558 : GEN z = nfmul(nf, y, nfsqr(nf, idealredmodpower(nf, y, p, 0)));
1189 5558 : z = Q_remove_denom(z, &den);
1190 5558 : if (den)
1191 : {
1192 28 : if (p != 2) den = powiu(den, p-1);
1193 28 : z = gmul(z, den);
1194 : }
1195 5558 : return z;
1196 : }
1197 : static GEN
1198 434 : nfV_factorbackmod(GEN nf, GEN x, ulong p)
1199 : {
1200 434 : long i, l = lg(x);
1201 434 : GEN v = cgetg(l, t_VEC);
1202 3752 : for (i = 1; i < l; i++)
1203 : {
1204 3318 : GEN y = gel(x,i), g = gel(y,1), e = gel(y,2);
1205 3318 : gel(v,i) = _factorbackmod(nf, g, ZV_to_Flv(e,p), p);
1206 : }
1207 434 : return v;
1208 : }
1209 : static GEN
1210 434 : nfV_zm_factorback(GEN nf, GEN G, GEN x, long p)
1211 2674 : { pari_APPLY_type(t_VEC, _factorbackmod(nf, G, gel(x,i), p)) }
1212 :
1213 : static GEN
1214 434 : prV_ZM_factorback(GEN nf, GEN S, GEN x)
1215 2674 : { pari_APPLY_type(t_VEC,idealfactorback(nf, S, gel(x,i), 0)) }
1216 :
1217 : /* shortcut for bnf = Q and p = 2 */
1218 : static GEN
1219 791 : bnfselmerQ(GEN S)
1220 : {
1221 791 : GEN g = vec_prepend(prV_primes(S), gen_m1), e;
1222 791 : long n = lg(S)-1;
1223 791 : e = n? shallowconcat(zerocol(n), matid(n)): zeromat(0, 1);
1224 791 : return mkvec3(g, e, const_vec(n + 1, gen_1));
1225 : }
1226 :
1227 : static GEN
1228 434 : bnfselmer(GEN bnf, GEN S)
1229 : {
1230 434 : const long p = 2;
1231 434 : pari_sp av = avma;
1232 434 : GEN nf = bnf_get_nf(bnf), S2, S3, e, f, e2, kerval, LS2gen, LS2fu, LS2all;
1233 434 : long n = lg(S)-1, n3, n2all, r;
1234 :
1235 434 : S2 = bestS(bnf, S, p);
1236 434 : S3 = shallowconcat(S, S2);
1237 434 : LS2all = nfV_factorbackmod(nf, gel(bnfunits(bnf, S3), 1), p);
1238 434 : n3 = lg(S3)-1; n2all = lg(LS2all)-1;
1239 434 : LS2gen = vecslice(LS2all,1,n3);
1240 434 : LS2fu = vecslice(LS2all,n3+1, n2all-1);
1241 434 : e2 = nfC_prV_val(nf, LS2gen, S2);
1242 434 : kerval = Flm_ker(ZM_to_Flm(e2, p), p);
1243 434 : LS2gen = nfV_zm_factorback(nf, LS2gen, kerval, p);
1244 434 : e = nfC_prV_val(nf, LS2gen, S);
1245 434 : e2 = ZM_divexactu(ZM_zm_mul(e2, kerval), p);
1246 434 : f = prV_ZM_factorback(nf, S2, e2);
1247 434 : LS2gen = shallowconcat(LS2fu, LS2gen);
1248 434 : LS2gen = nfV_to_scalar_or_alg(nf, vec_prepend(LS2gen, bnf_get_tuU(bnf)));
1249 434 : r = n2all - n3;
1250 434 : e = shallowconcat(zeromat(n, r), e);
1251 434 : f = shallowconcat(const_vec(r, gen_1), f);
1252 434 : return gerepilecopy(av, mkvec3(LS2gen,e,f));
1253 : }
1254 :
1255 : static GEN
1256 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, norms, 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 :
1277 245 : QS2gen = vec_prepend(SlistQ, gen_m1);
1278 245 : bad = ZV_sort_uniq_shallow(shallowconcat(factD, SlistQ));
1279 245 : Hlist = ZV_search(bad, gen_2)? bad: vec_prepend(bad, gen_2);
1280 245 : l = lg(QS2gen);
1281 245 : lHlist = lg(Hlist);
1282 245 : H = cgetg(l, t_MAT);
1283 2254 : for (j = 1; j < l; j++)
1284 : {
1285 2009 : GEN v = cgetg(lHlist, t_VECSMALL);
1286 29939 : for (i = 1; i < lHlist; i++)
1287 27930 : v[i] = hilbertii(D, gel(QS2gen, j), gel(Hlist, i)) < 0;
1288 2009 : gel(H, j) = Flv_to_F2v(v);
1289 : }
1290 245 : KerH = F2m_ker(H); l = lg(KerH);
1291 245 : norms = cgetg(l, t_VEC);
1292 1358 : for (i = 1; i < l; i++)
1293 1113 : gel(norms, i) = factorback2(QS2gen, F2c_to_ZC(gel(KerH, i)));
1294 245 : LS2gen = cgetg(l, t_VEC);
1295 245 : chpol = QXQ_charpoly(gel(nf_get_zk(nf), 2), nf_get_pol(nf), 0);
1296 245 : b = gdivgu(negi(gel(chpol, 3)), 2);
1297 245 : c = gel(chpol, 2);
1298 245 : Q = mkmat3(mkcol3(gen_1, b, gen_0),
1299 : mkcol3(b, c, gen_0),
1300 : mkcol3(gen_0, gen_0, gen_0));
1301 1358 : for (k = 1; k < l; k++)
1302 : {
1303 : GEN sol;
1304 1113 : gcoeff(Q, 3, 3) = negi(gel(norms, k));
1305 1113 : sol = qfsolve(Q); /* must be solvable */
1306 1113 : sol = Q_primpart(mkcol2(gel(sol,1), gel(sol,2)));
1307 1113 : gel(LS2gen, k) = basistoalg(nf, sol);
1308 : }
1309 245 : if (equalis(D, -4)) G = bad;
1310 : else
1311 : {
1312 245 : G = vecsplice(bad, ZV_search(bad, veclast(factD)));
1313 245 : G = vec_prepend(G, gen_m1);
1314 : }
1315 245 : LS2gen = shallowconcat(G, LS2gen);
1316 245 : l = lg(SlistQ); S2 = cgetg(l, t_VEC);
1317 245 : if (l > 1)
1318 : {
1319 2002 : for (i = 1; i < l; i++) gel(S2, i) = idealprimedec(nf, gel(SlistQ, i));
1320 238 : S2 = setminus(shallowconcat1(S2), S);
1321 : }
1322 245 : kerval = get_kerval(nf, S2, LS2gen); l = lg(kerval);
1323 245 : gen = cgetg(l, t_VEC);
1324 245 : e = cgetg(l, t_MAT);
1325 245 : f = cgetg(l, t_VEC);
1326 2639 : for (i = 1; i < l; i++)
1327 : {
1328 2394 : GEN id, ei, x = nffactorback(nf, LS2gen, F2v_to_Flv(gel(kerval, i)));
1329 2394 : gel(e,i) = ei = cgetg(lS, t_COL);
1330 33040 : for (j = 1; j < lS; j++) gel(ei,j) = stoi(idealval(nf, x, gel(S,j)));
1331 2394 : id = idealdiv(nf, x, idealfactorback(nf, S, gel(e,i), 0));
1332 2394 : if (!idealispower(nf, id, 2, &gel(f,i))) pari_err_BUG("nf2selmer_quad");
1333 2394 : gel(gen, i) = nf_to_scalar_or_alg(nf, x);
1334 : }
1335 245 : return gerepilecopy(ltop, mkvec3(gen, e, f));
1336 : }
1337 :
1338 : static GEN
1339 840 : makevbnf(GEN ell, long prec)
1340 : {
1341 840 : GEN vbnf, P = gel(ZX_factor(ell2pol(ell)), 1);
1342 840 : long k, l = lg(P);
1343 840 : vbnf = cgetg(l, t_VEC);
1344 2303 : for (k = 1; k < l; k++)
1345 : {
1346 1463 : GEN t = gel(P,k);
1347 1463 : gel(vbnf,k) = degpol(t) <= 2? nfinit(t, prec): Buchall(t, nf_FORCE, prec);
1348 : }
1349 840 : return vbnf;
1350 : }
1351 :
1352 : static GEN
1353 861 : kernorm(GEN G, GEN S, GEN pol, GEN signs)
1354 : {
1355 861 : long i, j, lS = lg(S), lG = lg(G), lv = signs? lS+2: lS+1;
1356 861 : GEN M = cgetg(lG, t_MAT);
1357 13720 : for (j = 1; j < lG; j++)
1358 : {
1359 12859 : GEN v, N = QXQ_norm(gel(G,j), pol);
1360 12859 : gel(M, j) = v = cgetg(lv, t_VECSMALL);
1361 12859 : v[1] = gsigne(N) < 0;
1362 181601 : for (i = 1; i < lS; i++) v[i+1] = odd(Q_pvalrem(N, gel(S,i), &N));
1363 12859 : if (signs) v[i+1] = signs[j];
1364 : }
1365 861 : return Flm_ker(M, 2);
1366 : }
1367 :
1368 : /* \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */
1369 : /* FUNCTIONS FOR 2-DESCENT \\ */
1370 : /* \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */
1371 : /* vector of t_VEC; return total number of entries */
1372 : static long
1373 8316 : RgVV_nb(GEN v)
1374 : {
1375 8316 : long i, l = lg(v), n = 0;
1376 24115 : for (i = 1; i < l; i++) n += lg(gel(v,i)) - 1;
1377 8316 : return n;
1378 : }
1379 :
1380 : /* return an Fp basis */
1381 : static GEN
1382 8316 : elllocalimage(GEN pol, GEN K, GEN vnf, GEN p, GEN pp, GEN pts)
1383 : {
1384 8316 : long n, p2 = RgVV_nb(pp), prank = equaliu(p, 2)? p2: p2 - 1;
1385 8316 : GEN R = polrootsmodpn(pol, p), bound = addiu(p, 6);
1386 :
1387 8316 : for(n = 1;; n++)
1388 21644 : {
1389 : pari_sp btop;
1390 : GEN x, y2, d;
1391 29960 : pts = Flm_image(pts, 2); if (lg(pts)-1 == prank) break;
1392 21644 : if ((n & 0xf) == 0) bound = mulii(bound, p);
1393 21644 : btop = avma;
1394 : do
1395 : {
1396 104923 : GEN r = gel(R, random_Fl(lg(R)-1)+1);
1397 104923 : long pprec = random_Fl(itou(gel(r, 2)) + 3) - 2; /* >= -2 */
1398 104923 : set_avma(btop);
1399 104923 : x = gadd(gel(r, 1), gmul(powis(p, pprec), randomi(bound)));
1400 104923 : y2 = gmul(K, poleval(pol, x));
1401 104923 : } while (gequal0(y2) || !Qp_issquare(y2, p));
1402 21644 : d = deg1pol_shallow(negi(K), gmul(K, x), 0);
1403 21644 : pts = vec_append(pts, kpmodsquares(vnf, d, pp));
1404 : }
1405 8316 : return pts;
1406 : }
1407 :
1408 : static GEN
1409 861 : ellLS2image(GEN pol, GEN vP, GEN K, GEN vpol, GEN vcrt)
1410 : {
1411 861 : long i, l = lg(vP);
1412 : GEN v;
1413 :
1414 861 : if (l == 1) return cgetg(1, t_VEC);
1415 728 : v = cgetg(l, t_VEC);
1416 59899 : for (i = 1; i < l; i++)
1417 : {
1418 59171 : GEN P = gel(vP, i), x, t;
1419 59171 : if (ell_is_inf(P)) { gel(v, i) = gen_1; continue; }
1420 59171 : x = gel(P,1);
1421 59171 : t = deg1pol_shallow(negi(K), gmul(K, x), 0);
1422 59171 : if (gequal0(gel(P,2)))
1423 : { /* 2-torsion, x now integer and a root of exactly one linear vpol[k]=T */
1424 567 : long k, lp = lg(vpol);
1425 : GEN a;
1426 1001 : for (k = 1; k < lp; k++)
1427 : {
1428 1001 : GEN T = gel(vpol,k), z = gel(T,2);
1429 1001 : if (absequalii(x, z) && signe(x) == -signe(z)) break; /* T = X-x */
1430 : }
1431 567 : a = ZX_Z_eval(ZX_deriv(pol), x);
1432 567 : t = gadd(a, gmul(gel(vcrt,k), gsub(t, a))); /* a mod T, t mod pol/T*/
1433 : }
1434 59171 : gel(v, i) = t;
1435 : }
1436 728 : return v;
1437 : }
1438 :
1439 : /* find small points on ell; 2-torsion points must be returned first */
1440 : static GEN
1441 861 : ellsearchtrivialpoints(GEN ell, GEN lim, GEN help)
1442 : {
1443 861 : pari_sp av = avma;
1444 861 : GEN tors2 = gel(elltors_psylow(ell,2),3);
1445 861 : GEN triv = lim ? ellratpoints(ell, lim, 0): cgetg(1,t_VEC);
1446 861 : if (help) triv = shallowconcat(triv, help);
1447 861 : return gerepilecopy(av, shallowconcat(tors2, triv));
1448 : }
1449 :
1450 : GEN
1451 63 : ellrankinit(GEN ell, long prec)
1452 : {
1453 63 : pari_sp av = avma;
1454 : GEN urst;
1455 63 : checkell_Q(ell); ell = ellminimalbmodel(ell, &urst);
1456 63 : return gerepilecopy(av, mkvec3(ell, urst, makevbnf(ell, prec)));
1457 : }
1458 :
1459 : INLINE GEN
1460 252252 : ZV_isneg(GEN x) { pari_APPLY_long(signe(gel(x,i)) < 0) }
1461 :
1462 : static GEN
1463 72030 : selmersign(GEN x, GEN vpol, GEN L)
1464 161056 : { pari_APPLY_same(ZV_isneg(nfeltsign(gel(x, i), RgX_rem(L, gel(vpol, i)), NULL))) }
1465 :
1466 : static GEN
1467 1722 : matselmersign(GEN vnf, GEN vpol, GEN x)
1468 73752 : { pari_APPLY_type(t_MAT, shallowconcat1(selmersign(vnf, vpol, gel(x, i)))) }
1469 :
1470 : static GEN
1471 117342 : _trace(GEN z, GEN T)
1472 : {
1473 117342 : long n = degpol(T)-1;
1474 117342 : if (degpol(z) < n) return gen_0;
1475 68058 : return gdiv(gel(z, 2+n), gel(T, 3+n));
1476 : }
1477 : static GEN
1478 19557 : tracematrix(GEN zc, GEN b, GEN T)
1479 : {
1480 : long i, j;
1481 19557 : GEN q = cgetg(4, t_MAT);
1482 78228 : for (j = 1; j <= 3; j++) gel(q,j) = cgetg(4, t_COL);
1483 78228 : for (j = 1; j <= 3; j++)
1484 : {
1485 117342 : for (i = 1; i < j; i++) gcoeff(q,i,j) = gcoeff(q,j,i) =
1486 58671 : _trace(QXQ_mul(zc, QXQ_mul(gel(b,i), gel(b,j), T), T), T);
1487 58671 : gcoeff(q,i,i) = _trace(QXQ_mul(zc, QXQ_sqr(gel(b,i), T), T), T);
1488 : }
1489 19557 : return q;
1490 : }
1491 :
1492 : static GEN
1493 21273 : RgXV_cxeval(GEN x, GEN r, GEN ri)
1494 85092 : { pari_APPLY_same(RgX_cxeval(gel(x,i), r, ri)) }
1495 :
1496 : static GEN
1497 7091 : redquadric(GEN base, GEN q2, GEN pol, GEN zc)
1498 : {
1499 7091 : long i, l, prec = nbits2prec(2*gexpo(q2)) + 1;
1500 7091 : GEN s = NULL, R = roots(pol, prec);
1501 7091 : l = lg(R);
1502 28364 : for (i = 1; i < l; ++i)
1503 : {
1504 21273 : GEN r = gel(R,i), ri = gexpo(r) > 1? ginv(r): NULL;
1505 21273 : GEN b = RgXV_cxeval(base, r, ri), z = RgX_cxeval(zc, r, ri);
1506 21273 : GEN M = RgC_RgV_mulrealsym(RgV_Rg_mul(b, gabs(z, prec)), gconj(b));
1507 21273 : s = s? RgM_add(s, M): M;
1508 : }
1509 7091 : return lllgram(s);
1510 : }
1511 :
1512 : static GEN
1513 19716 : RgX_homogenous_evaldeg(GEN P, GEN A, GEN B)
1514 : {
1515 19716 : long i, d = degpol(P), e = lg(B)-1;
1516 19716 : GEN s = gmul(gel(P, d+2), gel(B,e-d));
1517 61794 : for (i = d-1; i >= 0; i--)
1518 42078 : s = gadd(gmul(s, A), gmul(gel(B,e-i), gel(P,i+2)));
1519 19716 : return s;
1520 : }
1521 :
1522 : static GEN
1523 5375 : RgXV_homogenous_evaldeg(GEN x, GEN a, GEN b)
1524 21500 : { pari_APPLY_same(RgX_homogenous_evaldeg(gel(x,i), a, b)) }
1525 :
1526 : static void
1527 42 : check_oncurve(GEN ell, GEN v)
1528 : {
1529 42 : long i, l = lg(v);
1530 49 : for (i = 1; i < l; i++)
1531 : {
1532 7 : GEN P = gel(v, i);
1533 7 : if (lg(P) != 3 || !oncurve(ell,P)) pari_err_TYPE("ellrank",P);
1534 : }
1535 42 : }
1536 :
1537 : static GEN
1538 18998 : basis(GEN nf, GEN x, GEN crt, GEN pol)
1539 : {
1540 18998 : long i, l = lg(x);
1541 18998 : GEN b = cgetg(l, t_COL);
1542 42371 : for (i = 1; i < l; i++)
1543 : {
1544 23373 : GEN z = nf_to_scalar_or_alg(nf, gel(x, i));
1545 23373 : gel(b, i) = grem(gsub(z, gmul(crt, z)), pol); /* z mod T, 0 mod (pol/T) */
1546 : }
1547 18998 : return b;
1548 : }
1549 :
1550 : static GEN
1551 861 : vecsmallbasis(GEN x, GEN vcrt, GEN pol)
1552 2331 : { pari_APPLY_same(basis(gel(x,i), nf_get_zk(gel(x,i)), gel(vcrt,i), pol)) }
1553 :
1554 : static GEN
1555 266406 : ZC_shifti(GEN x, long n) { pari_APPLY_type(t_COL, shifti(gel(x,i), n)) }
1556 :
1557 : /* true nf */
1558 : static GEN
1559 17528 : selmerbasis(GEN nf, GEN ek, GEN sqrtLS2, GEN factLS2,
1560 : GEN badprimes, GEN crt, GEN pol)
1561 : {
1562 17528 : GEN sqrtzc = idealfactorback(nf, sqrtLS2, zv_neg(ek), 0);
1563 17528 : GEN E = ZC_shifti(ZM_zc_mul(factLS2, ek), -1);
1564 :
1565 17528 : if (ZV_equal0(E))
1566 15744 : sqrtzc = idealhnf_shallow(nf, sqrtzc);
1567 : else
1568 1784 : sqrtzc = idealmul(nf, sqrtzc, idealfactorback(nf, badprimes, ZC_neg(E), 0));
1569 17528 : return basis(nf, sqrtzc, crt, pol);
1570 : }
1571 :
1572 483 : static long randu(void) { return random_Fl(127) - 63; }
1573 : static GEN
1574 161 : randS(GEN b)
1575 : {
1576 483 : return gadd(gmulgs(gel(b,1), randu()),
1577 322 : gadd(gmulgs(gel(b,2), randu()), gmulgs(gel(b,3), randu())));
1578 : }
1579 :
1580 : static GEN
1581 6930 : liftselmerinit(GEN expo, GEN vnf, GEN sqrtLS2, GEN factLS2,
1582 : GEN badprimes, GEN vcrt, GEN pol)
1583 : {
1584 6930 : long n = lg(vnf)-1, k, t;
1585 6930 : GEN b = cgetg(n+1, t_VEC);
1586 24458 : for (k = t = 1; k <= n; k++)
1587 : {
1588 17528 : GEN fak = gel(factLS2,k), ek;
1589 17528 : long m = lg(fak)-1;
1590 17528 : ek = vecslice(expo, t, t + m-1); t += m;
1591 17528 : gel(b,k) = selmerbasis(gel(vnf,k), ek, gel(sqrtLS2,k),
1592 17528 : fak, gel(badprimes,k), gel(vcrt,k), pol);
1593 : }
1594 6930 : return shallowconcat1(b);
1595 : }
1596 :
1597 : static GEN
1598 3591 : liftselmer_cover(GEN b, GEN expo, GEN LS2, GEN pol, GEN discF, GEN K)
1599 : {
1600 : GEN P, Q, Q4, R, den, q0, q1, q2, xz, x, y, y2m, U, param, newb;
1601 : GEN ttheta, tttheta, zc, polprime;
1602 : GEN QM, zden;
1603 3591 : zc = RgXQV_factorback(LS2, expo, pol);
1604 3591 : if (typ(zc)==t_INT) zc = scalarpol(zc, varn(pol));
1605 3591 : ttheta = RgX_shift_shallow(pol,-2);
1606 3591 : tttheta = RgX_shift_shallow(pol, -1);
1607 3591 : polprime = ZX_deriv(pol);
1608 3591 : q2 = Q_primpart(tracematrix(zc, b, pol));
1609 3591 : U = redquadric(b, q2, pol, QXQ_div(zc, polprime, pol));
1610 3591 : q2 = qf_apply_RgM(q2, U);
1611 3591 : newb = RgV_RgM_mul(b, U);
1612 3591 : param = Q_primpart(qfparam(q2, qfsolve(q2), 1));
1613 3591 : param = RgM_to_RgXV_reverse(shallowtrans(param), 0);
1614 3591 : q1 = RgM_neg(tracematrix(QXQ_mul(zc, ttheta, pol), newb, pol));
1615 3591 : q1 = Q_remove_denom(qfeval(q1, param), &den);
1616 3591 : if (den) q1 = ZX_Z_mul(q1, den);
1617 3591 : if (!equali1(K)) q1 = RgX_Rg_mul(q1, K);
1618 3591 : QM = quartic_minim_all(q1, discF);
1619 3591 : q1 = gel(QM,1);
1620 3591 : zden = gmael(QM,2,1);
1621 3591 : Q = ZX_hyperellred(q1, &R);
1622 3591 : R = gmul(gmael(QM,2,2), R);
1623 3591 : if (DEBUGLEVEL>1) err_printf(" reduced quartic: Y^2 = %Ps\n", Q);
1624 3591 : xz = mkcol2(pol_x(0),gen_1);
1625 3591 : P = RgM_RgC_mul(R, xz); x = gel(P,1); y = gel(P,2);
1626 3591 : param = RgXV_homogenous_evaldeg(param, x, gpowers(y, 2));
1627 3591 : param = gmul(param, gdiv(den? mulii(K, den): K, zden));
1628 3591 : q0 = tracematrix(QXQ_mul(zc, tttheta, pol), newb, pol);
1629 3591 : x = gdiv(qfeval(q0, param), K);
1630 3591 : Q4 = gpowers(Q,4);
1631 3591 : y2m = gmul(RgX_homogenous_evaldeg(pol, x, Q4), Q);
1632 3591 : if (!issquareall(gdiv(y2m, K), &y))
1633 0 : pari_err_BUG("liftselmer_cover");
1634 3591 : y = gdiv(y, gel(Q4,2));
1635 3591 : P = mkvec2(gdiv(gmul(K,x),pol_xn(2,1)),gdiv(gmul(gsqr(K),y),pol_xn(3,1)));
1636 3591 : return mkvec2(Q,P);
1637 : }
1638 :
1639 : static GEN
1640 3339 : liftselmer(GEN b, GEN expo, GEN sbase, GEN LS2, GEN pol, GEN discF, GEN K, long ntry, GEN *pt_Q)
1641 : {
1642 3339 : pari_sp av = avma, av2;
1643 3339 : long t, lim = ntry * LIM1;
1644 : GEN ttheta, tttheta, z, polprime;
1645 : hashtable h;
1646 3339 : hash_init_GEN(&h, ntry, ZX_equal, 1);
1647 3339 : z = RgXQV_factorback(LS2, expo, pol);
1648 3339 : ttheta = RgX_shift_shallow(pol,-2);
1649 3339 : tttheta = RgX_shift_shallow(pol, -1);
1650 3339 : polprime = ZX_deriv(pol); av2 = avma;
1651 4005 : for (t = 1; t <= ntry; t++, set_avma(av2))
1652 : {
1653 : GEN P, Q, Qk, R, den, q0, q1, q2, xz, x, y, zz, zc, U, param, newb, zden, QM;
1654 : long idx;
1655 3500 : if (t == 1) zc = z;
1656 : else
1657 : {
1658 : GEN r;
1659 161 : do r = randS(sbase); while (degpol(QX_gcd(r, pol)));
1660 161 : zc = QXQ_mul(z, QXQ_sqr(r, pol), pol);
1661 : }
1662 3500 : q2 = Q_primpart(tracematrix(zc, b, pol));
1663 3500 : U = redquadric(b, q2, pol, QXQ_div(zc, polprime, pol));
1664 4166 : if (lg(U) < 4) continue;
1665 3500 : q2 = qf_apply_RgM(q2, U);
1666 3500 : newb = RgV_RgM_mul(b, U);
1667 :
1668 3500 : param = Q_primpart(qfparam(q2, qfsolve(q2), 1));
1669 3500 : param = RgM_to_RgXV_reverse(shallowtrans(param), 0);
1670 3500 : q1 = RgM_neg(tracematrix(QXQ_mul(zc, ttheta, pol), newb, pol));
1671 3500 : q1 = Q_remove_denom(qfeval(q1, param), &den);
1672 3500 : if (den) q1 = ZX_Z_mul(q1, den);
1673 3500 : if (!equali1(K)) q1 = RgX_Rg_mul(q1, K);
1674 3500 : QM = quartic_minim_all(q1, discF);
1675 3500 : q1 = gel(QM,1);
1676 3500 : zden = gmael(QM,2,1);
1677 3500 : Q = ZX_hyperellred(q1, &R);
1678 3500 : R = gmul(gmael(QM,2,2), R);
1679 3500 : if (pt_Q) *pt_Q = Q;
1680 3500 : Qk = shallowcopy(Q);
1681 3500 : (void) ZX_canon_neg(Qk);
1682 3500 : if (hash_haskey_long(&h, (void*)Qk, &idx)) continue;
1683 3443 : hash_insert_long(&h, (void*)Qk, 1); av2 = avma;
1684 3443 : if (DEBUGLEVEL>1) err_printf(" reduced quartic: Y^2 = %Ps\n", Q);
1685 :
1686 3443 : xz = projratpointxz(Q, lim, &zz);
1687 3443 : if (!xz)
1688 : {
1689 1694 : xz = projratpointxz2(Q, lim, &zz);
1690 1694 : if (!xz)
1691 : {
1692 3443 : if (pt_Q) return NULL; else continue;
1693 : }
1694 : }
1695 1784 : P = RgM_RgC_mul(R, xz); x = gel(P,1); y = gel(P,2);
1696 1784 : param = RgXV_homogenous_evaldeg(param, x, gpowers(y, 2));
1697 1784 : param = gmul(param, gdiv(den? mulii(K, den): K, gmul(zz, zden)));
1698 1784 : q0 = tracematrix(QXQ_mul(zc, tttheta, pol), newb, pol);
1699 1784 : x = gdiv(qfeval(q0, param), K);
1700 1784 : if (!issquareall(gdiv(poleval(pol, x), K), &y)) /* K y^2 = pol(x) */
1701 0 : pari_err_BUG("ellrank");
1702 1784 : P = mkvec2(x, y);
1703 1784 : if (DEBUGLEVEL) err_printf("Found point: %Ps\n", P);
1704 1784 : if (pt_Q) *pt_Q = gen_0;
1705 1784 : return gerepilecopy(av, P);
1706 : }
1707 505 : return NULL;
1708 : }
1709 :
1710 : static void
1711 1617 : gtoset_inplace(GEN x)
1712 1617 : { gen_sort_inplace(x, (void*)&cmp_universal, cmp_nodata, NULL); }
1713 :
1714 : /* FIXME: export */
1715 : static void
1716 8316 : setlgall(GEN x, long L)
1717 : {
1718 8316 : long i, l = lg(x);
1719 180201 : for(i = 1; i < l; i++) setlg(gel(x,i), L);
1720 8316 : }
1721 :
1722 : static long
1723 8316 : dim_selmer(GEN p, GEN pol, GEN K, GEN vnf, GEN LS2, GEN helpLS2,
1724 : GEN *selmer, GEN *LS2chars, GEN *helpchars)
1725 : {
1726 : pari_sp av;
1727 8316 : long dim, k, lvnf = lg(vnf);
1728 8316 : GEN X, L, LS2image, helpimage, pp = cgetg(lvnf, t_VEC);
1729 8316 : int pis2 = equaliu(p, 2);
1730 :
1731 24115 : for (k = 1; k < lvnf; k++)
1732 : {
1733 15799 : GEN v, nf = gel(vnf,k), PR = idealprimedec(nf, p);
1734 15799 : long j, l = lg(PR);
1735 15799 : gel(pp, k) = v = cgetg(l, t_VEC);
1736 35329 : for (j = 1; j < l; j++)
1737 : {
1738 19530 : GEN pr = gel(PR,j);
1739 19530 : gel(v,j) = pis2? log_prk_init(nf, pr, 1 + 2 * pr_get_e(pr), NULL)
1740 19530 : : zkmodprinit(nf, pr);
1741 : }
1742 : }
1743 8316 : LS2image = veckpmodsquares(LS2, vnf, pp);
1744 8316 : *LS2chars = vconcat(*LS2chars, LS2image);
1745 8316 : helpimage = veckpmodsquares(helpLS2, vnf, pp);
1746 8316 : *helpchars = vconcat(*helpchars, helpimage);
1747 8316 : av = avma;
1748 8316 : L = elllocalimage(pol, K, vnf, p, pp, helpimage);
1749 8316 : X = Flm_ker(shallowconcat(LS2image, L), 2); setlgall(X, lg(LS2image));
1750 : /* intersect(LS2image, locim) = LS2image.X */
1751 8316 : *selmer = Flm_intersect_i(*selmer, shallowconcat(Flm_ker(LS2image,2), X), 2);
1752 8316 : *selmer = gerepileupto(av, Flm_image(*selmer, 2));
1753 8316 : dim = lg(*selmer)-1; return (dim == Flm_rank(helpimage,2))? dim: -1;
1754 : }
1755 :
1756 : /* Assume there are 3 real roots, if K>0 return the smallest, otherwise the largest */
1757 : static long
1758 483 : get_row(GEN vnf, GEN K)
1759 : {
1760 483 : long k, sK = signe(K), n = lg(vnf)-1;
1761 : GEN R;
1762 483 : if (n == 1) return sK > 0? 1: 3;
1763 287 : if (n == 2)
1764 : {
1765 105 : GEN P = nf_get_pol(gel(vnf,2));
1766 105 : GEN z = negi(constant_coeff(nf_get_pol(gel(vnf,1))));
1767 105 : GEN y = poleval(P,z);
1768 105 : GEN b = gel(P,3), a = gel(P,4);
1769 105 : if (signe(y) != signe(a))
1770 : /* 1 is between 2 and 3 */
1771 7 : return sK > 0? 2: 3;
1772 98 : else if (cmpii(mulii(z,mulis(a,-2)), b) == signe(a))
1773 21 : return sK > 0? 1: 3;
1774 : else
1775 77 : return sK > 0? 2: 1;
1776 : }
1777 182 : R = cgetg(4, t_VEC);
1778 728 : for (k = 1; k <= 3; k++) gel(R, k) = gel(nf_get_roots(gel(vnf,k)), 1);
1779 182 : return sK > 0? vecindexmin(R): vecindexmax(R);
1780 : }
1781 :
1782 : static GEN
1783 861 : ell2selmer(GEN ell, GEN ell_K, GEN help, GEN K, GEN vbnf,
1784 : long effort, long flag, long prec)
1785 : {
1786 : GEN KP, pol, vnf, vpol, vcrt, sbase, LS2, factLS2, sqrtLS2, signs;
1787 : GEN selmer, helpLS2, LS2chars, helpchars, newselmer, factdisc, badprimes;
1788 : GEN helplist, listpoints, etors2, p, covers, disc, discF;
1789 861 : long i, k, n, tors2, mwrank, dim, nbpoints, lfactdisc, t, u, sha2 = 0;
1790 : forprime_t T;
1791 :
1792 861 : pol = ell2pol(ell);
1793 861 : help = ellsearchtrivialpoints(ell_K, flag ? NULL:muluu(LIMTRIV,effort+1), help);
1794 861 : help = elltwistpoints(help, ginv(K)); /* points on K Y^2 = pol(X) */
1795 861 : n = lg(vbnf) - 1; tors2 = n - 1;
1796 861 : etors2 = vecslice(help,1, tors2);
1797 861 : gtoset_inplace(etors2);
1798 861 : KP = gel(absZ_factor(K), 1);
1799 861 : disc = ZX_disc(pol);
1800 861 : factdisc = mkvec3(KP, mkcol(gen_2), gel(absZ_factor(disc), 1));
1801 861 : factdisc = ZV_sort_uniq_shallow(shallowconcat1(factdisc));
1802 861 : discF = mkvec2(gmul(K,disc), factdisc);
1803 861 : badprimes = cgetg(n+1, t_VEC);
1804 861 : vnf = cgetg(n+1, t_VEC);
1805 861 : vpol = cgetg(n+1, t_VEC);
1806 861 : vcrt = cgetg(n+1, t_VEC);
1807 861 : LS2 = cgetg(n+1, t_VEC);
1808 861 : factLS2 = cgetg(n+1, t_VEC);
1809 861 : sqrtLS2 = cgetg(n+1, t_VEC);
1810 2331 : for (k = 1; k <= n; k++)
1811 : {
1812 1470 : GEN T, Tinv, id, f, sel, L, S, nf, NF = gel(vbnf,k);
1813 : long i, lk;
1814 1470 : nf = (n == 1)? bnf_get_nf(NF): NF;
1815 1470 : gel(vnf, k) = nf;
1816 1470 : gel(vpol, k) = T = nf_get_pol(nf);
1817 1470 : Tinv = RgX_div(pol, gel(vpol, k));
1818 1470 : gel(vcrt, k) = QX_mul(QXQ_inv(T, Tinv), T); /* 0 mod T, 1 mod pol/T */
1819 :
1820 1470 : id = idealadd(nf, nf_get_index(nf), ZX_deriv(T));
1821 1470 : f = nf_pV_to_prV(nf, KP); settyp(f, t_COL);
1822 1470 : f = mkvec3(gel(idealfactor(nf, Tinv), 1),
1823 1470 : gel(idealfactor(nf, id), 1), f);
1824 1470 : gel(badprimes, k) = S = gtoset(shallowconcat1(f));
1825 1470 : if (n == 1)
1826 : {
1827 434 : sel = bnfselmer(NF, S);
1828 434 : obj_free(NF); /* units */
1829 : }
1830 1036 : else if (degpol(T) == 1)
1831 791 : sel = bnfselmerQ(S);
1832 : else /* degree 2 */
1833 245 : sel = nf2selmer_quad(NF, S);
1834 1470 : gel(LS2, k) = L = gel(sel, 1); lk = lg(L);
1835 1470 : gel(factLS2, k) = gel(sel, 2);
1836 1470 : gel(sqrtLS2, k) = gel(sel, 3);
1837 14329 : for (i = 1; i < lk; i++)
1838 : {
1839 12859 : GEN z = gel(L,i); /* z mod T, 1 mod (pol/T) */
1840 12859 : gel(L,i) = grem(gadd(z, gmul(gsubsg(1,z), gel(vcrt,k))), pol);
1841 : }
1842 : }
1843 861 : sbase = shallowconcat1(vecsmallbasis(vnf, vcrt, pol));
1844 861 : if (DEBUGLEVEL>2) err_printf(" local badprimes = %Ps\n", badprimes);
1845 861 : LS2 = shallowconcat1(LS2);
1846 861 : helpLS2 = ellLS2image(pol, help, K, vpol, vcrt);
1847 861 : LS2chars = matselmersign(vnf, vpol, LS2);
1848 861 : helpchars = matselmersign(vnf, vpol, helpLS2);
1849 861 : signs = NULL;
1850 861 : if (signe(ell_get_disc(ell)) > 0) signs = Flm_row(LS2chars, get_row(vnf,K));
1851 861 : selmer = kernorm(LS2, factdisc, pol, signs);
1852 861 : forprime_init(&T, gen_2, NULL); lfactdisc = lg(factdisc); dim = -1;
1853 7126 : for (i = 1; dim < 0 && i < lfactdisc; i++)
1854 6265 : dim = dim_selmer(gel(factdisc,i), pol, K, vnf, LS2, helpLS2,
1855 : &selmer,&LS2chars,&helpchars);
1856 2912 : while (dim < 0 && Flm_rank(Flm_mul(LS2chars, selmer, 2), 2) < lg(selmer)-1)
1857 : {
1858 2513 : while ((p = forprime_next(&T)) && ZV_search(factdisc, p));
1859 2051 : dim = dim_selmer(p, pol, K, vnf, LS2, helpLS2,
1860 : &selmer,&LS2chars,&helpchars);
1861 : }
1862 861 : helplist = gel(Flm_indexrank(helpchars,2), 2);
1863 861 : help = shallowextract(help, helplist);
1864 861 : helpchars = shallowextract(helpchars, helplist);
1865 861 : helpLS2 = shallowextract(helpLS2, helplist);
1866 861 : dim = lg(selmer)-1;
1867 861 : if (DEBUGLEVEL) err_printf("Selmer rank: %ld\n", dim);
1868 861 : newselmer = Flm_invimage(Flm_mul(LS2chars, selmer, 2), helpchars, 2);
1869 861 : nbpoints = lg(help) - 1;
1870 861 : if (flag==1)
1871 : {
1872 49 : GEN u = nbpoints? Flm_mul(selmer,Flm_suppl(newselmer,2), 2): selmer;
1873 49 : long l = lg(u);
1874 49 : GEN z = cgetg(l, t_VEC);
1875 154 : for (i = 1; i < l; i++) gel(z,i) = RgXQV_factorback(LS2, gel(u,i), pol);
1876 49 : return mkvec2(mkvec3(vnf,sbase,pol), z);
1877 : }
1878 812 : else if (flag==2)
1879 : {
1880 56 : GEN u = nbpoints ? Flm_mul(selmer,Flm_suppl(newselmer,2), 2): selmer;
1881 56 : long l = lg(u);
1882 56 : GEN P = cgetg(l, t_VEC), b;
1883 147 : for (i = 1; i < l; i++)
1884 : {
1885 91 : b = liftselmerinit(gel(u,i), vnf, sqrtLS2, factLS2, badprimes, vcrt, pol);
1886 91 : gel(P,i) = liftselmer_cover(b, gel(u,i), LS2, pol, discF, K);
1887 : }
1888 56 : return P;
1889 : }
1890 756 : listpoints = vec_lengthen(help, dim); /* points on ell */
1891 756 : covers = zerovec(dim);
1892 5677 : for (i=1; i <= dim; i++)
1893 : {
1894 4921 : GEN b, P, expo, vec = vecsmall_ei(dim, i);
1895 4921 : if (Flm_Flc_invimage(newselmer, vec, 2)) continue;
1896 2499 : expo = Flm_Flc_mul(selmer, vec, 2);
1897 2499 : b = liftselmerinit(expo, vnf, sqrtLS2, factLS2, badprimes, vcrt, pol);
1898 2499 : P = liftselmer(b, expo, sbase, LS2, pol, discF, K, 1, &gel(covers,i));
1899 2499 : if (P)
1900 : {
1901 1449 : gel(listpoints, ++nbpoints) = P; /* new point on ell */
1902 1449 : gel(newselmer, nbpoints) = vec;
1903 1449 : setlg(newselmer, nbpoints+1);
1904 : }
1905 : }
1906 756 : if (nbpoints < dim)
1907 : {
1908 : long i, j;
1909 301 : GEN M = cgetg(dim+1, t_MAT), selker;
1910 301 : GEN D = mulii(muliu(absi(disc), 27*4096), powiu(K,6));
1911 301 : GEN FD = ZV_sort_uniq_shallow(shallowconcat1(mkvec2(mkcol3s(3,5,7), factdisc)));
1912 :
1913 1911 : for (i = 1; i <= dim; i++) gel(M,i) = cgetg(dim+1, t_COL);
1914 1911 : for (i = 1; i <= dim; i++)
1915 9240 : for (j = 1; j <= i; j++)
1916 : {
1917 : GEN Q;
1918 7630 : if (isintzero(gel(covers,i)))
1919 3122 : Q = gen_0;
1920 4508 : else if (i==j)
1921 1008 : Q = quartic_findunit(D, gel(covers,i));
1922 : else
1923 : {
1924 3500 : GEN e = Flv_add(gel(selmer,i), gel(selmer,j), 2);
1925 3500 : GEN b = liftselmerinit(e, vnf, sqrtLS2, factLS2, badprimes, vcrt, pol);
1926 3500 : Q = quartic_findunit(D, gel(liftselmer_cover(b, e, LS2, pol, discF, K),1));
1927 : }
1928 7630 : gmael(M,j,i) = gmael(M,i,j) = Q;
1929 : }
1930 301 : selker = F2m_to_Flm(F2m_ker(matcassels(FD, M)));
1931 301 : sha2 = dim - (lg(selker)-1);
1932 301 : dim = lg(selker)-1;
1933 1141 : for (t=1, u=1; nbpoints < dim && effort > 0; t++)
1934 : {
1935 840 : pari_sp btop = avma;
1936 : GEN expo, b, P, vec;
1937 1119 : do vec = Flm_Flc_mul(selker,random_Flv(dim, 2), 2);
1938 1119 : while (zv_equal0(vec) || Flm_Flc_invimage(newselmer, vec, 2));
1939 840 : expo = Flm_Flc_mul(selmer, vec, 2);
1940 840 : b = liftselmerinit(expo, vnf, sqrtLS2, factLS2, badprimes, vcrt, pol);
1941 840 : P = liftselmer(b, expo, sbase, LS2, pol, discF, K, u, NULL);
1942 840 : if (P)
1943 : {
1944 335 : gel(listpoints, ++nbpoints) = P;
1945 335 : gel(newselmer, nbpoints) = vec;
1946 335 : setlg(newselmer, nbpoints+1);
1947 505 : } else set_avma(btop);
1948 840 : if (t == dim) { t = 0; u++; effort--; }
1949 : }
1950 : }
1951 756 : setlg(listpoints, nbpoints+1);
1952 756 : mwrank = nbpoints - tors2;
1953 756 : if (odd(dim - nbpoints)) mwrank++;
1954 756 : gtoset_inplace(listpoints);
1955 756 : listpoints = setminus(listpoints, etors2);
1956 756 : listpoints = elltwistpoints(listpoints, K);
1957 756 : listpoints = vecellabs(ellQ_genreduce(ell_K, listpoints, NULL, prec));
1958 756 : return mkvec4(utoi(mwrank), utoi(dim-tors2), utoi(sha2), listpoints);
1959 : }
1960 :
1961 : GEN
1962 49 : ell2selmer_basis(GEN ell, GEN *cb, long prec)
1963 : {
1964 49 : GEN E = ellminimalbmodel(ell, cb);
1965 49 : GEN S = ell2selmer(E, E, NULL, gen_1, makevbnf(E, prec), 0, 1, prec);
1966 49 : obj_free(E); return S;
1967 : }
1968 :
1969 : static void
1970 98 : ellrank_get_nudur(GEN E, GEN F, GEN *nu, GEN *du, GEN *r)
1971 : {
1972 98 : GEN ea2 = ell_get_a2(E), ea2t = ell_get_a2(F);
1973 98 : GEN ec4 = ell_get_c4(E), ec4t = ell_get_c4(F);
1974 98 : GEN ec6 = ell_get_c6(E), ec6t = ell_get_c6(F);
1975 : GEN N, D, d;
1976 98 : if (signe(ec4)==0)
1977 : {
1978 21 : if (!Z_ispowerall(mulii(ec6,sqri(ec6t)),3,&N))
1979 7 : pari_err_TYPE("ellrank",F);
1980 14 : D = ec6t;
1981 : }
1982 77 : else if (signe(ec6)==0)
1983 : {
1984 21 : if (!Z_issquareall(mulii(ec4,ec4t),&N))
1985 7 : pari_err_TYPE("ellrank",F);
1986 14 : D = ec4t;
1987 : }
1988 : else
1989 : {
1990 56 : GEN d46 = mulii(gcdii(ec4,ec4t),gcdii(ec6,ec6t));
1991 56 : N = diviiexact(mulii(ec6,ec4t),d46);
1992 56 : D = diviiexact(mulii(ec6t,ec4),d46);
1993 : }
1994 84 : d = gcdii(N, D);
1995 84 : *nu = diviiexact(N, d); *du = diviiexact(D, d);
1996 84 : *r = diviuexact(subii(mulii(*nu,ea2t),mulii(*du,ea2)),3);
1997 84 : }
1998 :
1999 : static GEN
2000 854 : ellrank_flag(GEN e, long effort, GEN help, long flag, long prec)
2001 : {
2002 854 : pari_sp ltop = avma;
2003 : GEN urst, v, vbnf, eK;
2004 854 : GEN et = NULL, K = gen_1, nu = NULL, du = NULL, r = NULL, urstK = NULL;
2005 854 : long newell = 0;
2006 :
2007 854 : if (lg(e)==3 && typ(e)==t_VEC) { et = gel(e,2); e = gel(e,1); }
2008 854 : if (lg(e)==4 && typ(e)==t_VEC)
2009 : {
2010 126 : vbnf = gel(e,3); urst = gel(e,2);
2011 126 : e = gel(e,1); checkell_Q(e);
2012 126 : if (!ell_is_integral(e)) pari_err_TYPE("ellrank [nonintegral model]",e);
2013 119 : if (signe(ell_get_a1(e))) pari_err_TYPE("ellrank [a1 != 0]", e);
2014 112 : if (signe(ell_get_a3(e))) pari_err_TYPE("ell2rank [a3 != 0]", e);
2015 : }
2016 : else
2017 : {
2018 728 : checkell_Q(e);
2019 728 : e = ellminimalbmodel(e, &urst);
2020 728 : newell = 1;
2021 728 : vbnf = makevbnf(e, prec);
2022 : }
2023 833 : if (et)
2024 : {
2025 105 : checkell_Q(et);
2026 105 : if (!gequal(ell_get_j(et),ell_get_j(e))) pari_err_TYPE("ellrank",et);
2027 98 : et = ellminimalbmodel(et, &urst);
2028 98 : ellrank_get_nudur(e, et, &nu, &du, &r);
2029 84 : K = mulii(nu, du);
2030 84 : urstK = mkvec4(nu, mulii(nu,r), gen_0, gen_0);
2031 : }
2032 812 : if (help)
2033 : {
2034 42 : if (urst) help = ellchangepoint(help, urst);
2035 42 : if (et) help = ellchangepointinv(help, urstK);
2036 : }
2037 812 : eK = elltwistequation(e, K);
2038 : /* help is a vector of rational points [x,y] satisfying K y^2 = pol(x) */
2039 : /* [Kx, K^2y] is on eK: Y^2 = K^3 pol(X/K) */
2040 812 : if (help) check_oncurve(eK, help);
2041 812 : v = ell2selmer(e, eK, help, K, vbnf, effort, flag, prec);
2042 812 : if (flag==0)
2043 : {
2044 756 : if (et) gel(v,4) = ellchangepoint(gel(v,4), urstK);
2045 756 : if (urst) gel(v,4) = ellchangepointinv(gel(v,4), urst);
2046 : }
2047 : else
2048 : {
2049 56 : long i, l = lg(v);
2050 147 : for (i = 1; i < l; i++)
2051 : {
2052 91 : if (et) gmael(v,i,2) = ellchangepoint(gmael(v,i,2), urstK);
2053 91 : if (urst) gmael(v,i,2) = ellchangepointinv(gmael(v,i,2), urst);
2054 : }
2055 : }
2056 812 : if (newell) obj_free(e);
2057 812 : if (eK != e) obj_free(eK);
2058 812 : return gerepilecopy(ltop, v);
2059 : }
2060 :
2061 : GEN
2062 798 : ellrank(GEN e, long effort, GEN help, long prec)
2063 : {
2064 798 : return ellrank_flag(e, effort, help, 0, prec);
2065 : }
2066 :
2067 : GEN
2068 56 : ell2cover(GEN ell, long prec)
2069 : {
2070 56 : return ellrank_flag(ell, 0, NULL, 2, prec);
2071 : }
|
