Line data Source code
1 : /* Copyright (C) 2000 The PARI group.
2 :
3 : This file is part of the PARI/GP package.
4 :
5 : PARI/GP is free software; you can redistribute it and/or modify it under the
6 : terms of the GNU General Public License as published by the Free Software
7 : Foundation; either version 2 of the License, or (at your option) any later
8 : version. It is distributed in the hope that it will be useful, but WITHOUT
9 : ANY WARRANTY WHATSOEVER.
10 :
11 : Check the License for details. You should have received a copy of it, along
12 : with the package; see the file 'COPYING'. If not, write to the Free Software
13 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14 :
15 : /********************************************************************/
16 : /** **/
17 : /** TORSION OF ELLIPTIC CURVES over NUMBER FIELDS **/
18 : /** **/
19 : /********************************************************************/
20 : #include "pari.h"
21 : #include "paripriv.h"
22 : static int
23 623 : smaller_x(GEN p, GEN q)
24 : {
25 623 : int s = abscmpii(denom_i(p), denom_i(q));
26 623 : return (s<0 || (s==0 && abscmpii(numer_i(p),numer_i(q)) < 0));
27 : }
28 :
29 : /* best generator in cycle of length k */
30 : static GEN
31 686 : best_in_cycle(GEN e, GEN p, long k)
32 : {
33 686 : GEN p0 = p,q = p;
34 : long i;
35 :
36 1057 : for (i=2; i+i<k; i++)
37 : {
38 371 : q = elladd(e,q,p0);
39 371 : if (ugcd(i,k)==1 && smaller_x(gel(q,1), gel(p,1))) p = q;
40 : }
41 686 : return (gsigne(ec_dmFdy_evalQ(e,p)) < 0)? ellneg(e,p): p;
42 : }
43 :
44 : /* <p,q> = E_tors, possibly NULL (= oo), p,q independent unless NULL
45 : * order p = k, order q = 2 unless NULL */
46 : static GEN
47 2226 : tors(GEN e, long k, GEN p, GEN q, GEN v)
48 : {
49 : GEN r;
50 2226 : if (q)
51 : {
52 280 : long n = k>>1;
53 280 : GEN p1, best = q, np = ellmul(e,p,utoipos(n));
54 280 : if (n % 2 && smaller_x(gel(np,1), gel(best,1))) best = np;
55 280 : p1 = elladd(e,q,np);
56 280 : if (smaller_x(gel(p1,1), gel(best,1))) q = p1;
57 259 : else if (best == np) { p = elladd(e,p,q); q = np; }
58 280 : p = best_in_cycle(e,p,k);
59 280 : if (v)
60 : {
61 0 : p = ellchangepointinv(p,v);
62 0 : q = ellchangepointinv(q,v);
63 : }
64 280 : r = cgetg(4,t_VEC);
65 280 : gel(r,1) = utoipos(2*k);
66 280 : gel(r,2) = mkvec2(utoipos(k), gen_2);
67 280 : gel(r,3) = mkvec2copy(p, q);
68 : }
69 : else
70 : {
71 1946 : if (p)
72 : {
73 406 : p = best_in_cycle(e,p,k);
74 406 : if (v) p = ellchangepointinv(p,v);
75 406 : r = cgetg(4,t_VEC);
76 406 : gel(r,1) = utoipos(k);
77 406 : gel(r,2) = mkvec( gel(r,1) );
78 406 : gel(r,3) = mkvec( gcopy(p) );
79 : }
80 : else
81 : {
82 1540 : r = cgetg(4,t_VEC);
83 1540 : gel(r,1) = gen_1;
84 1540 : gel(r,2) = cgetg(1,t_VEC);
85 1540 : gel(r,3) = cgetg(1,t_VEC);
86 : }
87 : }
88 2226 : return r;
89 : }
90 :
91 : /* Finds a multiplicative upper bound for #E_tor (p-Sylow if p != 0);
92 : * assume integral model */
93 : static long
94 2296 : torsbound(GEN e, ulong p)
95 : {
96 2296 : GEN D = ell_get_disc(e);
97 2296 : pari_sp av = avma, av2;
98 2296 : long m, b, bold, nb, CM = ellQ_get_CM(e);
99 : forprime_t S;
100 2296 : nb = expi(D) >> 3; /* number of primes to try ~ 1 prime every 8 bits in D */
101 2296 : switch (p)
102 : {
103 343 : case 0: b = 5040; break;
104 938 : case 2: b = 16; break;
105 49 : case 3: b = 9; break;
106 84 : case 5: case 7: b = p; break;
107 882 : default: return 1;
108 : }
109 1414 : bold = b;
110 1414 : m = 0;
111 : /* p > 2 has good reduction => E(Q) injects in E(Fp) */
112 1414 : (void)u_forprime_init(&S, 3, ULONG_MAX);
113 1414 : av2 = avma;
114 23730 : while (m < nb || (b > 12 && b != 16))
115 : {
116 22897 : ulong p = u_forprime_next(&S);
117 22897 : if (!p) pari_err_BUG("torsbound [ran out of primes]");
118 22897 : if (!umodiu(D, p)) continue;
119 :
120 18242 : b = ugcd(b, p+1 - ellap_CM_fast(e,p,CM));
121 18242 : set_avma(av2);
122 18242 : if (b == 1) break;
123 17661 : if (b == bold) m++; else { bold = b; m = 0; }
124 : }
125 1414 : return gc_long(av,b);
126 : }
127 :
128 : /* return a rational point of order pk = p^k on E, or NULL if E(Q)[k] = O.
129 : * *fk is either NULL (pk = 4 or prime) or elldivpol(p^(k-1)).
130 : * Set *fk to elldivpol(p^k) */
131 : static GEN
132 630 : tpoint(GEN E, long pk, GEN *fk)
133 : {
134 630 : GEN f = elldivpol(E,pk,0), g = *fk, v;
135 : long i, l;
136 630 : *fk = f;
137 630 : if (g) f = RgX_div(f, g);
138 630 : v = nfrootsQ(f); l = lg(v);
139 966 : for (i = 1; i < l; i++)
140 : {
141 791 : GEN x = gel(v,i);
142 791 : GEN y = ellordinate(E,x,0);
143 791 : if (lg(y) != 1) return mkvec2(x,gel(y,1));
144 : }
145 175 : return NULL;
146 : }
147 : /* return E(Q)[2] */
148 : static GEN
149 651 : t2points(GEN E, GEN *f2)
150 : {
151 : long i, l;
152 : GEN v;
153 651 : *f2 = ec_bmodel(E,0);
154 651 : v = nfrootsQ(*f2); l = lg(v);
155 1848 : for (i = 1; i < l; i++)
156 : {
157 1197 : GEN x = gel(v,i);
158 1197 : GEN y = ellordinate(E,x,0);
159 1197 : if (lg(y) != 1) gel(v,i) = mkvec2(x,gel(y,1));
160 : }
161 651 : return v;
162 : }
163 :
164 : /* psylow = 0 or prime (return p-Sylow subgroup) */
165 : static GEN
166 2226 : ellQtors(GEN E, long psylow)
167 : {
168 2226 : GEN T2 = NULL, p, P, Q, v;
169 : long v2, r2, B;
170 :
171 2226 : E = ellintegralmodel_i(E, &v);
172 2226 : B = torsbound(E, psylow); /* #E_tor | B */
173 2226 : if (B == 1) return tors(E,1,NULL,NULL, v);
174 763 : v2 = vals(B); /* bound for v_2(point order) */
175 763 : B >>= v2;
176 763 : p = const_vec(9, NULL);
177 763 : r2 = 0;
178 763 : if (v2) {
179 : GEN f;
180 651 : T2 = t2points(E, &f);
181 651 : switch(lg(T2)-1)
182 : {
183 14 : case 0: v2 = 0; break;
184 357 : case 1: r2 = 1; if (v2 == 4) v2 = 3; break;
185 280 : default: r2 = 2; v2--; break; /* 3 */
186 : }
187 651 : if (v2) gel(p,2) = gel(T2,1);
188 : /* f = f2 */
189 651 : if (v2 > 1) { gel(p,4) = tpoint(E,4, &f); if (!gel(p,4)) v2 = 1; }
190 : /* if (v2>1) now f = f4 */
191 651 : if (v2 > 2) { gel(p,8) = tpoint(E,8, &f); if (!gel(p,8)) v2 = 2; }
192 : }
193 763 : B <<= v2;
194 763 : if (B % 3 == 0) {
195 91 : GEN f3 = NULL;
196 91 : gel(p,3) = tpoint(E,3,&f3);
197 91 : if (!gel(p,3)) B /= (B%9)? 3: 9;
198 91 : if (gel(p,3) && B % 9 == 0)
199 : {
200 7 : gel(p,9) = tpoint(E,9,&f3);
201 7 : if (!gel(p,9)) B /= 3;
202 : }
203 : }
204 763 : if (B % 5 == 0) {
205 63 : GEN junk = NULL;
206 63 : gel(p,5) = tpoint(E,5,&junk);
207 63 : if (!gel(p,5)) B /= 5;
208 : }
209 763 : if (B % 7 == 0) {
210 35 : GEN junk = NULL;
211 35 : gel(p,7) = tpoint(E,7,&junk);
212 35 : if (!gel(p,7)) B /= 7;
213 : }
214 : /* B is the exponent of E_tors(Q), r2 is the rank of its 2-Sylow,
215 : * for i > 1, p[i] is a point of order i if one exists and i is a prime power
216 : * and NULL otherwise */
217 763 : if (r2 == 2) /* 2 cyclic factors */
218 : { /* C2 x C2 */
219 280 : if (B == 2) return tors(E,2, gel(T2,1), gel(T2,2), v);
220 119 : else if (B == 6)
221 : { /* C2 x C6 */
222 14 : P = elladd(E, gel(p,3), gel(T2,1));
223 14 : Q = gel(T2,2);
224 : }
225 : else
226 : { /* C2 x C4 or C2 x C8 */
227 105 : P = gel(p, B);
228 105 : Q = gel(T2,2);
229 105 : if (gequal(Q, ellmul(E, P, utoipos(B>>1)))) Q = gel(T2,1);
230 : }
231 : }
232 : else /* cyclic */
233 : {
234 483 : Q = NULL;
235 483 : if (v2)
236 : {
237 357 : if (B>>v2 == 1)
238 308 : P = gel(p, B);
239 : else
240 49 : P = elladd(E, gel(p, B>>v2), gel(p,1<<v2));
241 : }
242 126 : else P = gel(p, B);
243 : }
244 602 : return tors(E,B, P, Q, v);
245 : }
246 :
247 : /* either return one prime of degree 1 above p or NULL (none or expensive) */
248 : static GEN
249 227213 : primedec_deg1(GEN K, GEN p)
250 : {
251 227213 : GEN r, T, f = nf_get_index(K);
252 227213 : if (dvdii(f,p)) return NULL;
253 227136 : T = nf_get_pol(K);
254 227136 : r = FpX_oneroot(T, p); if (!r) return NULL;
255 135387 : r = deg1pol_shallow(gen_1, Fp_neg(r,p), varn(T));
256 135387 : return idealprimedec_kummer(K, r, 1, p);
257 : }
258 :
259 : /* Bound for the elementary divisors of the torsion group of elliptic curve
260 : * (p-Sylow if psylow = p is not 0)
261 : * Reduce the curve modulo some small good primes */
262 : static GEN
263 29148 : nftorsbound(GEN E, ulong psylow)
264 : {
265 : pari_sp av;
266 29148 : long k = 0, g;
267 29148 : GEN B1 = gen_0, B2 = gen_0, K = ellnf_get_nf(E);
268 29148 : GEN D = ell_get_disc(E), ND = idealnorm(K,D);
269 : forprime_t S;
270 29148 : if (typ(ND) == t_FRAC) ND = gel(ND,1);
271 29148 : ND = mulii(ND, Q_denom(vecslice(E,1,5)));
272 29148 : g = maxss(5, expi(ND) >> 3);
273 29148 : if (g > 20) g = 20;
274 : /* P | p such that e(P/p) < p-1 => E(K) injects in E(k(P)) [otherwise
275 : * we may lose some p-torsion]*/
276 29148 : (void)u_forprime_init(&S, 3, ULONG_MAX);
277 29148 : av = avma;
278 229096 : while (k < g) /* k = number of good primes already used */
279 : {
280 225043 : ulong p = u_forprime_next(&S);
281 : GEN P, gp;
282 : long j, l;
283 225043 : if (!umodiu(ND,p)) continue;
284 171745 : gp = utoipos(p);
285 : /* primes of degree 1 are easier and give smaller bounds */
286 171745 : if (typ(D) != t_POLMOD) /* E/Q */
287 : {
288 169162 : P = primedec_deg1(K, gp); /* single P|p has all the information */
289 169162 : if (!P) continue;
290 106428 : P = mkvec(P);
291 : }
292 : else
293 2583 : P = idealprimedec_limit_f(K, utoipos(p), 1);
294 109011 : l = lg(P);
295 192164 : for (j = 1; j < l; j++,k++)
296 : {
297 108248 : GEN Q = gel(P,j), EQ, cyc;
298 : long n;
299 108248 : if ((ulong)pr_get_e(Q) >= p-1) continue;
300 108248 : EQ = ellinit(E,zkmodprinit(K,Q),0);
301 108248 : cyc = ellgroup(EQ, NULL);
302 108248 : n = lg(cyc)-1;
303 108248 : if (n == 0) return mkvec2(gen_1,gen_1);
304 108241 : B1 = gcdii(B1,gel(cyc,1));
305 108241 : B2 = (n == 1)? gen_1: gcdii(B2,gel(cyc,2));
306 108241 : obj_free(EQ);
307 : /* division by 2 is cheap when it fails, no need to have a sharp bound */
308 108241 : if (psylow==0 && Z_ispow2(B1)) return mkvec2(B1,B2);
309 : }
310 83916 : if ((g & 15) == 0) gerepileall(av, 2, &B1, &B2);
311 : }
312 4053 : if (abscmpiu(B2, 2) > 0)
313 : { /* if E(K) has full n-torsion then K contains the n-th roots of 1 */
314 91 : GEN n = gel(nfrootsof1(K), 1);
315 91 : B2 = gcdii(B2,n);
316 : }
317 4053 : if (psylow)
318 : {
319 0 : B1 = powuu(psylow, Z_lval(B1, psylow));
320 0 : B2 = powuu(psylow, Z_lval(B2, psylow));
321 : }
322 4053 : return mkvec2(B1,B2);
323 : }
324 :
325 : /* Checks whether the point P != oo is divisible by n in E(K), where xn is
326 : * [phi_n, psi_n^2]. If so, returns a point Q such that nQ = P or -P.
327 : * Else, returns NULL */
328 : static GEN
329 364 : ellnfis_divisible_by(GEN E, GEN K, GEN P, GEN xn)
330 : {
331 364 : GEN R = nfroots(K, RgX_sub(RgX_Rg_mul(gel(xn,2), gel(P,1)), gel(xn,1)));
332 364 : long i, l = lg(R);
333 364 : for (i = 1; i < l; i++)
334 : {
335 252 : GEN x = gel(R,i), y = ellordinate(E,x,0);
336 252 : if (lg(y) != 1) return mkvec2(x, gel(y,1));
337 : }
338 112 : return NULL;
339 : }
340 :
341 : /* Q is not the point at infinity; do we have Q = [n]Q' on curve C/K ?
342 : * return P such that [n] P = \pm Q or NULL if none exist */
343 : static GEN
344 49 : ellnfis_divisible_by_i(GEN C, GEN K, GEN Q, GEN n, long v)
345 : {
346 : GEN xn;
347 49 : if (!isprimepower(absi_shallow(n), NULL))
348 : {
349 14 : GEN f = absZ_factor(n), P = gel(f,1), E = gel(f,2);
350 14 : long i, l = lg(P);
351 14 : for (i = 1; i < l; i++)
352 : {
353 0 : GEN q = powii(gel(P,i), gel(E,i));
354 0 : xn = ellxn(C, itou(q), v);
355 0 : if (!(Q = ellnfis_divisible_by(C, K, Q, xn))) return 0;
356 : }
357 14 : return Q;
358 : }
359 35 : xn = ellxn(C, itou(n), v);
360 35 : return ellnfis_divisible_by(C, K, Q, xn);
361 : }
362 :
363 : static GEN ellorder_nf(GEN E, GEN P, GEN B);
364 : long
365 273 : ellisdivisible(GEN E, GEN Q, GEN n, GEN *pQ)
366 : {
367 273 : pari_sp av = avma;
368 273 : GEN P = NULL, N = NULL, K = NULL;
369 : long v;
370 :
371 273 : checkell(E);
372 273 : if (!checkellpt_i(Q)) pari_err_TYPE("ellisdivisible",Q);
373 273 : switch(ell_get_type(E))
374 : {
375 252 : case t_ELL_Q: break;
376 21 : case t_ELL_NF: K = ellnf_get_nf(E); break;
377 0 : default: pari_err_TYPE("ellisdivisible",E);
378 : }
379 273 : if (ell_is_inf(Q))
380 : {
381 7 : if (pQ) *pQ = ellinf();
382 7 : return 1;
383 : }
384 266 : switch(typ(n))
385 : {
386 126 : case t_INT:
387 126 : if (!signe(n)) return 0;
388 119 : if (is_bigint(n))
389 : { /* only possible if P is torsion */
390 28 : GEN x, o, ex = cyc_get_expo(abgrp_get_cyc(elltors(E)));
391 28 : if (equali1(ex)) return gc_long(av, 0);
392 28 : o = K? ellorder_nf(E, Q, ex): ellorder(E, Q, NULL);
393 28 : if (isintzero(o)) return gc_long(av, 0);
394 21 : x = Z_ppo(n, o);
395 21 : Q = ellmul(E, Q, Fp_inv(x, o));
396 21 : n = diviiexact(n, x); /* all primes dividing n divide o = ord(Q) | ex */
397 21 : if (!dvdii(diviiexact(ex,o), n)) return gc_long(av, 0);
398 7 : if (is_bigint(n)) pari_err_IMPL("ellisdivisible for huge torsion");
399 7 : if (!ellisdivisible(E, Q, n, pQ)) return gc_long(av, 0);
400 7 : if (!pQ) return gc_long(av, 1);
401 0 : *pQ = gerepilecopy(av, *pQ); return 1;
402 : }
403 91 : if (!K)
404 : { /* could use elltors instead of a multiple ? */
405 70 : ulong nn = itou(n), n2 = u_ppo(nn, torsbound(E, 0));
406 70 : if (n2 > 1)
407 : { /* n2 coprime to torsion */
408 42 : if (!(Q = ellQ_isdivisible(E, Q, n2))) return 0;
409 28 : if (signe(n) < 0) Q = ellneg(E, Q);
410 28 : if (nn == n2)
411 : {
412 28 : if (pQ) *pQ = Q;
413 28 : return 1;
414 : }
415 : /* we may have changed n into |n/n2| (and Q in -Q if n < 0) */
416 0 : n = utoipos(nn/n2);
417 : }
418 : }
419 49 : v = fetch_var_higher();
420 49 : P = ellnfis_divisible_by_i(E, K, Q, n, v);
421 49 : delete_var(); N = n; break;
422 140 : case t_VEC:
423 140 : if (lg(n) == 3 && typ(gel(n,1)) == t_POL && typ(gel(n,2)) == t_POL)
424 : { /* ellxn */
425 140 : long d, d2 = degpol(gel(n,1));
426 140 : if (d2 < 0) return gc_long(av, 0);
427 140 : if (!uissquareall(d2,(ulong*)&d)) pari_err_TYPE("ellisdivisible",n);
428 140 : P = ellnfis_divisible_by(E, K, Q, n);
429 140 : N = utoi(d); break;
430 : } /* fall through */
431 : default:
432 0 : pari_err_TYPE("ellisdivisible",n);
433 0 : break;
434 : }
435 189 : if (!P) return gc_long(av, 0);
436 154 : if (!pQ) return gc_long(av, 1);
437 133 : if (gequal(Q, ellmul(E, P, N)))
438 70 : P = gerepilecopy(av, P);
439 : else
440 63 : P = gerepileupto(av, ellneg(E, P));
441 133 : *pQ = P; return 1;
442 : }
443 :
444 : /* 2-torsion point of abscissa x */
445 : static GEN
446 126 : tor2(GEN E, GEN x) { return mkvec2(x, gmul2n(gneg(ec_h_evalx(E,x)), -1)); }
447 :
448 : static GEN
449 35 : ptor0(void)
450 35 : { return mkvec2(mkvec(gen_1),cgetg(1,t_VEC)); }
451 : static GEN
452 126 : ptor1(long p, long n, GEN P)
453 126 : { return mkvec2(mkvec(powuu(p,n)), mkvec(P)); }
454 : static GEN
455 98 : ptor2(long p, long n1, long n2, GEN P1, GEN P2)
456 98 : { return mkvec2(mkvec2(powuu(p,n1), powuu(p,n2)), mkvec2(P1,P2)); }
457 :
458 : /* Computes the p-primary torsion in E(K). Assume that p is small, should use
459 : * Weil pairing otherwise.
460 : * N1, N2 = upper bounds on the integers n1 >= n2 such that
461 : * E(K)[p^oo] = Z/p^n1 x Z/p^n2
462 : * Returns [cyc,gen], where E(K)[p^oo] = sum Z/cyc[i] gen[i] */
463 : static GEN
464 259 : ellnftorsprimary(GEN E, long p, long N1, long N2, long v)
465 : {
466 259 : GEN X, P1, P2, Q1, Q2, xp, K = ellnf_get_nf(E);
467 : long n1, n2;
468 :
469 : /* compute E[p] = < P1 > or < P1, P2 > */
470 259 : X = nfroots(K, elldivpol(E,p,v));
471 259 : if (lg(X) == 1) return ptor0();
472 231 : P2 = ellinf();
473 231 : if (p==2)
474 : {
475 77 : P1 = tor2(E, gel(X,1));
476 77 : if (lg(X) > 2) P2 = tor2(E, gel(X,2)); /* E[2] = (Z/2Z)^2 */
477 : }
478 : else
479 : {
480 154 : long j, l = lg(X), nT, a;
481 154 : GEN T = vectrunc_init(l);
482 539 : for(j=1; j < l; j++)
483 : {
484 385 : GEN a = gel(X,j), Y = ellordinate(E,a,0);
485 385 : if (lg(Y) != 1) vectrunc_append(T, mkvec2(a,gel(Y,1)));
486 : }
487 154 : nT = lg(T)-1;
488 154 : if (!nT) return ptor0();
489 147 : P1 = gel(T,1);
490 147 : a = (p-1)/2;
491 147 : if (nT != a)
492 : { /* E[p] = (Z/pZ)^2 */
493 49 : GEN Z = cgetg(a+1,t_VEC), Q1 = P1;
494 : long k;
495 49 : gel(Z,1) = Q1;
496 49 : for (k=2; k <= a; k++) gel(Z,k) = elladd(E,Q1,P1);
497 49 : gen_sort_inplace(Z, (void*)&cmp_universal, &cmp_nodata, NULL);
498 49 : while (tablesearch(Z, gel(T,k), &cmp_universal)) k++;
499 49 : P2 = gel(T,k);
500 : }
501 : }
502 224 : xp = ellxn(E, p, v);
503 :
504 224 : if (ell_is_inf(P2))
505 : { /* E[p^oo] is cyclic, start from P1 and divide by p while possible */
506 140 : for (n1 = 1; n1 < N1; n1++)
507 : {
508 14 : GEN Q = ellnfis_divisible_by(E,K,P1,xp);
509 14 : if (!Q) break;
510 14 : P1 = Q;
511 : }
512 126 : return ptor1(p, n1, P1);
513 : }
514 :
515 : /* E[p] = (Z/pZ)^2, compute n2 and E[p^n2] */
516 98 : Q1 = NULL;
517 119 : for (n2 = 1; n2 < N2; n2++)
518 : {
519 21 : Q1 = ellnfis_divisible_by(E,K,P1,xp);
520 21 : Q2 = ellnfis_divisible_by(E,K,P2,xp);
521 21 : if (!Q1 || !Q2) break;
522 21 : P1 = Q1;
523 21 : P2 = Q2;
524 : }
525 :
526 : /* compute E[p^oo] = < P1, P2 > */
527 98 : n1 = n2;
528 98 : if (n2 == N2)
529 : {
530 98 : if (N1 == N2) return ptor2(p, n2,n2, P1,P2);
531 42 : Q1 = ellnfis_divisible_by(E,K,P1,xp);
532 : }
533 42 : if (Q1) { P1 = Q1; n1++; }
534 : else
535 : {
536 28 : Q2 = ellnfis_divisible_by(E,K,P2,xp);
537 28 : if (Q2) { P2 = P1; P1 = Q2; n1++; }
538 : else
539 : {
540 : long k;
541 35 : for (k = 1; k < p; k++)
542 : {
543 28 : P1 = elladd(E,P1,P2);
544 28 : Q1 = ellnfis_divisible_by(E,K,P1,xp);
545 28 : if (Q1) { P1 = Q1; n1++; break; }
546 : }
547 28 : if (k == p) return ptor2(p, n2,n2, P1,P2);
548 : }
549 : }
550 : /* P1,P2 of order p^n1,p^n2 with n1=n2+1.
551 : * Keep trying to divide P1 + k P2 with 0 <= k < p by p */
552 56 : while (n1 < N1)
553 : {
554 21 : Q1 = ellnfis_divisible_by(E,K,P1,xp);
555 21 : if (Q1) { P1 = Q1; n1++; }
556 : else
557 : {
558 : long k;
559 14 : for (k = 1; k < p; k++)
560 : {
561 14 : P1 = elladd(E,P1,P2);
562 14 : Q1 = ellnfis_divisible_by(E,K,P1,xp);
563 14 : if (Q1) { P1 = Q1; n1++; break; }
564 : }
565 14 : if (k == p) break;
566 : }
567 : }
568 35 : return ptor2(p, n1,n2, P1,P2);
569 : }
570 :
571 : /* P affine point */
572 : static GEN
573 259 : nfpt(GEN e, GEN P)
574 : {
575 259 : GEN T = nf_get_pol(ellnf_get_nf(e));
576 259 : GEN x = gel(P,1), y = gel(P,2);
577 259 : long tx = typ(x), ty = typ(y);
578 259 : if (tx == ty) return P;
579 98 : if (tx != t_POLMOD) x = mkpolmod(x,T); else y = mkpolmod(y,T);
580 98 : return mkvec2(x,y);
581 : }
582 : /* Computes the torsion subgroup of E(K), as [order, cyc, gen] */
583 : static GEN
584 189 : ellnftors(GEN e, ulong psylow)
585 : {
586 189 : GEN B = nftorsbound(e, psylow), B1 = gel(B,1), B2 = gel(B,2), d1,d2, P1,P2;
587 189 : GEN f = Z_factor(B1), P = gel(f,1), E = gel(f,2);
588 189 : long i, l = lg(P), v = fetch_var_higher();
589 :
590 189 : d1 = d2 = gen_1; P1 = P2 = ellinf();
591 448 : for (i=1; i<l; i++)
592 : {
593 259 : long p = itos(gel(P,i)); /* Compute p-primary torsion */
594 259 : long N1 = itos(gel(E,i)); /* >= n1 */
595 259 : long N2 = Z_lval(B2,p); /* >= n2 */
596 259 : GEN T = ellnftorsprimary(e, p, N1, N2, v), cyc = gel(T,1), gen = gel(T,2);
597 259 : if (is_pm1(gel(cyc,1))) continue;
598 : /* update generators P1,P2 and their respective orders d1,d2 */
599 224 : P1 = elladd(e, P1, gel(gen,1)); d1 = mulii(d1, gel(cyc,1));
600 224 : if (lg(cyc) > 2)
601 98 : { P2 = elladd(e, P2, gel(gen,2)); d2 = mulii(d2, gel(cyc,2)); }
602 : }
603 189 : (void)delete_var();
604 189 : if (is_pm1(d1)) return mkvec3(gen_1,cgetg(1,t_VEC),cgetg(1,t_VEC));
605 168 : if (is_pm1(d2)) return mkvec3(d1, mkvec(d1), mkvec(nfpt(e,P1)));
606 91 : return mkvec3(mulii(d1,d2), mkvec2(d1,d2), mkvec2(nfpt(e,P1),nfpt(e,P2)));
607 : }
608 :
609 : GEN
610 462 : elltors(GEN e)
611 : {
612 462 : pari_sp av = avma;
613 462 : GEN t = NULL;
614 462 : checkell(e);
615 462 : switch(ell_get_type(e))
616 : {
617 273 : case t_ELL_Q: t = ellQtors(e, 0); break;
618 189 : case t_ELL_NF: t = ellnftors(e, 0); break;
619 0 : case t_ELL_Fp:
620 0 : case t_ELL_Fq: return ellgroup0(e,NULL,1);
621 0 : default: pari_err_TYPE("elltors",e);
622 : }
623 462 : return gerepilecopy(av, t);
624 : }
625 :
626 : GEN
627 0 : elltors0(GEN e, long flag) { (void)flag; return elltors(e); }
628 :
629 : GEN
630 1953 : elltors_psylow(GEN e, ulong p)
631 : {
632 1953 : pari_sp av = avma;
633 1953 : GEN t = NULL;
634 1953 : checkell(e);
635 1953 : switch(ell_get_type(e))
636 : {
637 1953 : case t_ELL_Q: t = ellQtors(e, p); break;
638 0 : case t_ELL_NF: t = ellnftors(e, p); break;
639 0 : default: pari_err_TYPE("elltors_psylow",e);
640 : }
641 1953 : return gerepilecopy(av, t);
642 : }
643 :
644 : /********************************************************************/
645 : /** **/
646 : /** ORDER OF POINTS over NUMBER FIELDS **/
647 : /** **/
648 : /********************************************************************/
649 : /* E a t_ELL_Q (use Mazur's theorem) */
650 : long
651 52648 : ellorder_Q(GEN E, GEN P)
652 : {
653 52648 : pari_sp av = avma;
654 : GEN dx, dy, d4, d6, D, Pp, Q;
655 : forprime_t S;
656 : ulong a4, p;
657 : long k;
658 52648 : if (ell_is_inf(P)) return 1;
659 52648 : if (gequal(P, ellneg(E,P))) return 2;
660 :
661 52627 : dx = Q_denom(gel(P,1));
662 52627 : dy = Q_denom(gel(P,2));
663 52627 : if (ell_is_integral(E)) /* integral model, try Nagell Lutz */
664 52564 : if (abscmpiu(dx, 4) > 0 || abscmpiu(dy, 8) > 0) return 0;
665 :
666 35277 : d4 = Q_denom(ell_get_c4(E));
667 35277 : d6 = Q_denom(ell_get_c6(E));
668 35277 : D = ell_get_disc (E);
669 : /* choose not too small prime p dividing neither a coefficient of the
670 : short Weierstrass form nor of P and leading to good reduction */
671 35277 : u_forprime_init(&S, 100003, ULONG_MAX);
672 35277 : while ( (p = u_forprime_next(&S)) )
673 35277 : if (umodiu(d4, p) && umodiu(d6, p) && Rg_to_Fl(D, p)
674 35277 : && umodiu(dx, p) && umodiu(dy, p)) break;
675 :
676 : /* transform E into short Weierstrass form Ep modulo p and P to Pp on Ep,
677 : * check whether the order of Pp on Ep is <= 12 */
678 35277 : Pp = point_to_a4a6_Fl(E, P, p, &a4);
679 35277 : for (Q = Fle_dbl(Pp, a4, p), k = 2;
680 422708 : !ell_is_inf(Q) && k <= 12;
681 387431 : Q = Fle_add(Q, Pp, a4, p), k++) /* empty */;
682 :
683 35277 : if (k == 13) k = 0;
684 : else
685 : { /* check whether [k]P = O over Q. Save potentially costly last elladd */
686 : GEN R;
687 84 : Q = ellmul(E, P, utoipos(k>>1));
688 84 : R = odd(k)? elladd(E, P,Q): Q;
689 84 : if (!gequal(Q, ellneg(E,R))) k = 0;
690 : }
691 35277 : return gc_long(av,k);
692 : }
693 : /* E a t_ELL_NF, B a multiplicative bound for exponent(E_tor) or NULL */
694 : static GEN
695 28994 : ellorder_nf(GEN E, GEN P, GEN B)
696 : {
697 28994 : GEN K = ellnf_get_nf(E);
698 28994 : pari_sp av = avma;
699 : GEN dx, dy, d4, d6, D, ND, Ep, Pp, Q, gp, modpr, pr, T, k;
700 : forprime_t S;
701 : ulong a4, p;
702 28994 : if (ell_is_inf(P)) return gen_1;
703 28994 : if (gequal(P, ellneg(E,P))) return gen_2;
704 :
705 28959 : if (!B) B = gel(nftorsbound(E, 0), 1);
706 28959 : dx = Q_denom(gel(P,1));
707 28959 : dy = Q_denom(gel(P,2));
708 28959 : d4 = Q_denom(ell_get_c4(E));
709 28959 : d6 = Q_denom(ell_get_c6(E));
710 28959 : D = ell_get_disc(E);
711 28959 : ND = idealnorm(K,D);
712 28959 : if (typ(ND) == t_FRAC) ND = gel(ND,1);
713 :
714 : /* choose not too small prime p of degree 1 dividing neither a coefficient of
715 : * the short Weierstrass form nor of P and leading to good reduction */
716 28959 : u_forprime_init(&S, 100003, ULONG_MAX);
717 58051 : while ( (p = u_forprime_next(&S)) )
718 : {
719 58051 : if (!umodiu(d4, p) || !umodiu(d6, p) || !umodiu(ND, p)
720 58051 : || !umodiu(dx, p) || !umodiu(dy, p)) continue;
721 58051 : gp = utoipos(p);
722 58051 : pr = primedec_deg1(K, gp);
723 58051 : if (pr) break;
724 : }
725 :
726 28959 : modpr = nf_to_Fq_init(K, &pr,&T,&gp);
727 28959 : Ep = ellinit(E, pr, 0);
728 28959 : Pp = nfV_to_FqV(P, K, modpr);
729 :
730 : /* transform E into short Weierstrass form Ep modulo p and P to Pp on Ep,
731 : * check whether the order of Pp on Ep divides B */
732 28959 : Pp = point_to_a4a6_Fl(Ep, Pp, p, &a4);
733 28959 : if (!ell_is_inf(Fle_mul(Pp, B, a4, p))) { set_avma(av); return gen_0; }
734 147 : k = Fle_order(Pp, B, a4, p);
735 : { /* check whether [k]P = O over K. Save potentially costly last elladd */
736 : GEN R;
737 147 : Q = ellmul(E, P, shifti(k,-1));
738 147 : R = mod2(k)? elladd(E, P,Q): Q;
739 147 : if (!gequal(Q, ellneg(E,R))) k = gen_0;
740 : }
741 147 : return gerepileuptoint(av, k);
742 : }
743 :
744 : GEN
745 31556 : ellorder(GEN E, GEN P, GEN o)
746 : {
747 31556 : pari_sp av = avma;
748 31556 : GEN fg, r, E0 = E;
749 31556 : checkell(E);
750 31556 : if (!checkellpt_i(P)) pari_err_TYPE("ellorder",P);
751 31556 : if (ell_is_inf(P)) return gen_1;
752 31542 : if (ell_get_type(E)==t_ELL_Q)
753 : {
754 1897 : long tx = typ(gel(P,1)), ty = typ(gel(P,2));
755 1897 : GEN p = NULL;
756 1897 : if (is_rational_t(tx) && is_rational_t(ty)) return utoi(ellorder_Q(E, P));
757 1778 : if (tx == t_INTMOD || tx == t_FFELT) p = gel(P,1);
758 1778 : if (!p && (ty == t_INTMOD || ty == t_FFELT)) p = gel(P,2);
759 1778 : if (p)
760 : {
761 1778 : E = ellinit(E,p,0);
762 1778 : if (lg(E)==1) pari_err_IMPL("ellorder for curve with singular reduction");
763 : }
764 : }
765 31416 : if (ell_get_type(E)==t_ELL_NF) return ellorder_nf(E, P, NULL);
766 2422 : checkell_Fq(E);
767 2422 : fg = ellff_get_field(E);
768 2422 : if (!o) o = ellff_get_o(E);
769 2422 : if (typ(fg)==t_FFELT)
770 1764 : r = FF_ellorder(E, P, o);
771 : else
772 : {
773 658 : GEN p = fg, e = ellff_get_a4a6(E);
774 658 : GEN Pp = FpE_changepointinv(RgE_to_FpE(P,p), gel(e,3), p);
775 658 : r = FpE_order(Pp, o, gel(e,1), p);
776 : }
777 2422 : if (E != E0) obj_free(E);
778 2422 : return gerepileuptoint(av, r);
779 : }
780 :
781 : GEN
782 0 : orderell(GEN e, GEN z) { return ellorder(e,z,NULL); }
|
