Line data Source code
1 : /* Copyright (C) 2012 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 : #include "pari.h"
16 : #include "paripriv.h"
17 :
18 : #define DEBUGLEVEL DEBUGLEVEL_ellcard
19 :
20 : /* Not so fast arithmetic with points over elliptic curves over Fq,
21 : small characteristic. */
22 :
23 : /***********************************************************************/
24 : /** **/
25 : /** FlxqE **/
26 : /** **/
27 : /***********************************************************************/
28 : /* These functions deal with point over elliptic curves over Fq defined
29 : * by an equation of the form y^2=x^3+a4*x+a6. Most of the time a6 is omitted
30 : * since it can be recovered from any point on the curve. */
31 :
32 : GEN
33 65804 : RgE_to_FlxqE(GEN x, GEN T, ulong p)
34 : {
35 65804 : if (ell_is_inf(x)) return x;
36 65804 : retmkvec2(Rg_to_Flxq(gel(x,1),T,p), Rg_to_Flxq(gel(x,2),T,p));
37 : }
38 :
39 : GEN
40 157822 : FlxqE_changepoint(GEN x, GEN ch, GEN T, ulong p)
41 : {
42 157822 : pari_sp av = avma;
43 : GEN p1, p2, z, u, r, s, t, v, v2, v3;
44 : ulong pi;
45 157822 : if (ell_is_inf(x)) return x;
46 93810 : pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
47 93810 : u = gel(ch,1); r = gel(ch,2);
48 93810 : s = gel(ch,3); t = gel(ch,4);
49 93810 : v = Flxq_inv_pre(u, T, p, pi);
50 93810 : v2 = Flxq_sqr_pre(v, T, p, pi);
51 93810 : v3 = Flxq_mul_pre(v,v2, T, p, pi);
52 93810 : p1 = Flx_sub(gel(x,1), r, p);
53 93810 : p2 = Flx_sub(gel(x,2), Flx_add(Flxq_mul_pre(s, p1, T, p, pi),t, p), p);
54 93810 : z = cgetg(3,t_VEC);
55 93810 : gel(z,1) = Flxq_mul_pre(v2, p1, T, p, pi);
56 93810 : gel(z,2) = Flxq_mul_pre(v3, p2, T, p, pi);
57 93810 : return gerepileupto(av, z);
58 : }
59 :
60 : GEN
61 65804 : FlxqE_changepointinv(GEN x, GEN ch, GEN T, ulong p)
62 : {
63 65804 : pari_sp av = avma;
64 : GEN p1, p2, u, r, s, t, X, Y, u2, u3, u2X, z;
65 : ulong pi;
66 65804 : if (ell_is_inf(x)) return x;
67 65804 : pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
68 65804 : X = gel(x,1); Y = gel(x,2);
69 65804 : u = gel(ch,1); r = gel(ch,2);
70 65804 : s = gel(ch,3); t = gel(ch,4);
71 65804 : u2 = Flxq_sqr_pre(u, T, p, pi);
72 65804 : u3 = Flxq_mul_pre(u,u2, T, p, pi);
73 65804 : u2X = Flxq_mul_pre(u2,X, T, p, pi);
74 65804 : p1 = Flxq_mul_pre(u3,Y, T, p, pi);
75 65804 : p2 = Flx_add(Flxq_mul_pre(s,u2X, T, p, pi), t, p);
76 65804 : z = cgetg(3, t_VEC);
77 65804 : gel(z,1) = Flx_add(u2X, r, p);
78 65804 : gel(z,2) = Flx_add(p1, p2, p);
79 65804 : return gerepileupto(av, z);
80 : }
81 :
82 : static GEN
83 22834 : nonsquare_Flxq(GEN T, ulong p)
84 : {
85 22834 : pari_sp av = avma;
86 22834 : long n = degpol(T), vs = T[1];
87 : GEN a;
88 22834 : if (odd(n))
89 7686 : return mkvecsmall2(vs, nonsquare_Fl(p));
90 : do
91 : {
92 30485 : set_avma(av);
93 30485 : a = random_Flx(n, vs, p);
94 30485 : } while (Flxq_issquare(a, T, p));
95 15148 : return a;
96 : }
97 :
98 : void
99 22834 : Flxq_elltwist(GEN a, GEN a6, GEN T, ulong p, GEN *pt_a, GEN *pt_a6)
100 : {
101 22834 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
102 22834 : GEN d = nonsquare_Flxq(T, p);
103 22834 : GEN d2 = Flxq_sqr_pre(d, T, p, pi), d3 = Flxq_mul_pre(d2, d, T, p, pi);
104 22834 : if (typ(a)==t_VECSMALL)
105 : {
106 15232 : *pt_a = Flxq_mul_pre(a, d2, T, p, pi);
107 15232 : *pt_a6 = Flxq_mul_pre(a6, d3, T, p, pi);
108 : } else
109 : {
110 7602 : *pt_a = mkvec(Flxq_mul_pre(gel(a,1), d, T, p, pi));
111 7602 : *pt_a6 = Flxq_mul_pre(a6, d3, T, p, pi);
112 : }
113 22834 : }
114 :
115 : static GEN
116 1594453 : FlxqE_dbl_slope(GEN P, GEN a4, GEN T, ulong p, ulong pi, GEN *ps)
117 : {
118 : GEN x, y, Q, s;
119 1594453 : if (ell_is_inf(P) || !lgpol(gel(P,2))) return ellinf();
120 1477510 : x = gel(P,1); y = gel(P,2);
121 1477510 : if (p==3UL)
122 531265 : s = typ(a4)==t_VEC? Flxq_div_pre(Flxq_mul_pre(x, gel(a4,1), T,p,pi), y, T,p,pi)
123 531265 : : Flxq_div_pre(a4, Flx_neg(y, p), T,p,pi);
124 : else
125 : {
126 946245 : GEN sx = Flx_add(Flx_triple(Flxq_sqr_pre(x, T, p, pi), p), a4, p);
127 946245 : s = Flxq_div_pre(sx, Flx_double(y, p), T, p, pi);
128 : }
129 1477510 : Q = cgetg(3,t_VEC);
130 1477510 : gel(Q,1) = Flx_sub(Flxq_sqr_pre(s, T, p, pi), Flx_double(x, p), p);
131 1477510 : if (typ(a4)==t_VEC) gel(Q, 1) = Flx_sub(gel(Q,1), gel(a4,1), p);
132 1477510 : gel(Q,2) = Flx_sub(Flxq_mul_pre(s, Flx_sub(x, gel(Q,1), p), T, p, pi),
133 : y, p);
134 1477510 : if (ps) *ps = s;
135 1477510 : return Q;
136 : }
137 :
138 : GEN
139 0 : FlxqE_dbl(GEN P, GEN a4, GEN T, ulong p)
140 : {
141 0 : pari_sp av = avma;
142 0 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
143 0 : return gerepileupto(av, FlxqE_dbl_slope(P,a4, T, p, pi, NULL));
144 : }
145 :
146 : static GEN
147 713511 : FlxqE_add_slope(GEN P, GEN Q, GEN a4, GEN T, ulong p, ulong pi, GEN *ps)
148 : {
149 : GEN Px, Py, Qx, Qy, R, s;
150 713511 : if (ell_is_inf(P)) return Q;
151 708906 : if (ell_is_inf(Q)) return P;
152 708720 : Px = gel(P,1); Py = gel(P,2);
153 708720 : Qx = gel(Q,1); Qy = gel(Q,2);
154 708720 : if (Flx_equal(Px, Qx))
155 : {
156 53679 : if (Flx_equal(Py, Qy))
157 4310 : return FlxqE_dbl_slope(P, a4, T, p, pi, ps);
158 : else
159 49369 : return ellinf();
160 : }
161 655041 : s = Flxq_div_pre(Flx_sub(Py, Qy, p), Flx_sub(Px, Qx, p), T, p, pi);
162 655041 : R = cgetg(3,t_VEC);
163 655041 : gel(R,1) = Flx_sub(Flx_sub(Flxq_sqr_pre(s, T, p, pi), Px, p), Qx, p);
164 655041 : if (typ(a4)==t_VEC) gel(R,1) = Flx_sub(gel(R,1), gel(a4,1), p);
165 655041 : gel(R,2) = Flx_sub(Flxq_mul_pre(s, Flx_sub(Px, gel(R,1), p), T, p, pi),
166 : Py, p);
167 655041 : if (ps) *ps = s;
168 655041 : return R;
169 : }
170 :
171 : GEN
172 0 : FlxqE_add(GEN P, GEN Q, GEN a4, GEN T, ulong p)
173 : {
174 0 : pari_sp av = avma;
175 0 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
176 0 : return gerepileupto(av, FlxqE_add_slope(P,Q,a4, T,p,pi, NULL));
177 : }
178 :
179 : static GEN
180 6802 : FlxqE_neg_i(GEN P, ulong p)
181 : {
182 6802 : if (ell_is_inf(P)) return P;
183 6802 : return mkvec2(gel(P,1), Flx_neg(gel(P,2), p));
184 : }
185 :
186 : GEN
187 1193 : FlxqE_neg(GEN P, GEN T, ulong p)
188 : {
189 : (void) T;
190 1193 : if (ell_is_inf(P)) return ellinf();
191 1193 : return mkvec2(gcopy(gel(P,1)), Flx_neg(gel(P,2), p));
192 : }
193 :
194 : GEN
195 0 : FlxqE_sub(GEN P, GEN Q, GEN a4, GEN T, ulong p)
196 : {
197 0 : pari_sp av = avma;
198 0 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
199 0 : return gerepileupto(av, FlxqE_add_slope(P, FlxqE_neg_i(Q, p), a4, T,p,pi, NULL));
200 : }
201 :
202 : struct _FlxqE
203 : {
204 : GEN a4, a6, T;
205 : ulong p, pi;
206 : };
207 :
208 : static GEN
209 1566276 : _FlxqE_dbl(void *E, GEN P)
210 : {
211 1566276 : struct _FlxqE *e = (struct _FlxqE *) E;
212 1566276 : return FlxqE_dbl_slope(P, e->a4, e->T, e->p, e->pi, NULL);
213 : }
214 :
215 : static GEN
216 703860 : _FlxqE_add(void *E, GEN P, GEN Q)
217 : {
218 703860 : struct _FlxqE *e = (struct _FlxqE *) E;
219 703860 : return FlxqE_add_slope(P, Q, e->a4, e->T, e->p, e->pi, NULL);
220 : }
221 :
222 : static GEN
223 6802 : _FlxqE_sub(void *E, GEN P, GEN Q)
224 : {
225 6802 : struct _FlxqE *e = (struct _FlxqE *) E;
226 6802 : return FlxqE_add_slope(P, FlxqE_neg_i(Q,e->p), e->a4, e->T,e->p,e->pi, NULL);
227 : }
228 :
229 : static GEN
230 265633 : _FlxqE_mul(void *E, GEN P, GEN n)
231 : {
232 265633 : pari_sp av = avma;
233 265633 : struct _FlxqE *e=(struct _FlxqE *) E;
234 265633 : long s = signe(n);
235 265633 : if (!s || ell_is_inf(P)) return ellinf();
236 264928 : if (s < 0) P = FlxqE_neg(P, e->T, e->p);
237 264928 : if (is_pm1(n)) return s>0? gcopy(P): P;
238 251685 : return gerepilecopy(av, gen_pow_i(P, n, e, &_FlxqE_dbl, &_FlxqE_add));
239 : }
240 :
241 : GEN
242 64068 : FlxqE_mul(GEN P, GEN n, GEN a4, GEN T, ulong p)
243 : {
244 : struct _FlxqE E;
245 64068 : E.a4= a4; E.T = T; E.p = p; E.pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
246 64068 : return _FlxqE_mul(&E, P, n);
247 : }
248 :
249 : /* 3*x^2+2*a2*x = -a2*x, and a2!=0 */
250 :
251 : /* Finds a random nonsingular point on E */
252 : static GEN
253 76839 : random_F3xqE(GEN a2, GEN a6, GEN T)
254 : {
255 76839 : pari_sp ltop = avma;
256 : GEN x, y, rhs;
257 76839 : const ulong p = 3;
258 : do
259 : {
260 151221 : set_avma(ltop);
261 151221 : x = random_Flx(get_Flx_degree(T),get_Flx_var(T),p);
262 151221 : rhs = Flx_add(Flxq_mul(Flxq_sqr(x, T, p), Flx_add(x, a2, p), T, p), a6, p);
263 151221 : } while ((!lgpol(rhs) && !lgpol(x)) || !Flxq_issquare(rhs, T, p));
264 76839 : y = Flxq_sqrt(rhs, T, p);
265 76839 : if (!y) pari_err_PRIME("random_F3xqE", T);
266 76839 : return gerepilecopy(ltop, mkvec2(x, y));
267 : }
268 :
269 : /* Finds a random nonsingular point on E */
270 : static GEN
271 151101 : random_FlxqE_pre(GEN a4, GEN a6, GEN T, ulong p, ulong pi)
272 : {
273 151101 : pari_sp ltop = avma;
274 : GEN x, x2, y, rhs;
275 151101 : if (typ(a4)==t_VEC) return random_F3xqE(gel(a4,1), a6, T);
276 : do
277 : {
278 143553 : set_avma(ltop);
279 143553 : x = random_Flx(get_Flx_degree(T),get_Flx_var(T),p);
280 143553 : x2 = Flxq_sqr_pre(x, T, p, pi); /* x^3+a4*x+a6 = x*(x^2+a4)+a6 */
281 143553 : rhs = Flx_add(Flxq_mul_pre(x, Flx_add(x2, a4, p), T, p, pi), a6, p);
282 149300 : } while ((!lgpol(rhs) && !lgpol(Flx_add(Flx_triple(x2, p), a4, p)))
283 149258 : || !Flxq_issquare(rhs, T, p));
284 74262 : y = Flxq_sqrt(rhs, T, p);
285 74262 : if (!y) pari_err_PRIME("random_FlxqE", T);
286 74262 : return gerepilecopy(ltop, mkvec2(x, y));
287 : }
288 : GEN
289 76941 : random_FlxqE(GEN a4, GEN a6, GEN T, ulong p)
290 76941 : { return random_FlxqE_pre(a4, a6, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
291 :
292 : static GEN
293 74160 : _FlxqE_rand(void *E)
294 : {
295 74160 : struct _FlxqE *e=(struct _FlxqE *) E;
296 74160 : return random_FlxqE_pre(e->a4, e->a6, e->T, e->p, e->pi);
297 : }
298 :
299 : static const struct bb_group FlxqE_group={_FlxqE_add,_FlxqE_mul,_FlxqE_rand,hash_GEN,zvV_equal,ell_is_inf, NULL};
300 :
301 : const struct bb_group *
302 61 : get_FlxqE_group(void ** pt_E, GEN a4, GEN a6, GEN T, ulong p)
303 : {
304 61 : struct _FlxqE *e = (struct _FlxqE *) stack_malloc(sizeof(struct _FlxqE));
305 61 : e->a4 = a4; e->a6 = a6;
306 61 : e->pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
307 61 : e->p = p;
308 61 : e->T = Flx_get_red_pre(T, p, e->pi);
309 61 : *pt_E = (void *)e; return &FlxqE_group;
310 : }
311 :
312 : GEN
313 1470 : FlxqE_order(GEN z, GEN o, GEN a4, GEN T, ulong p)
314 : {
315 1470 : pari_sp av = avma;
316 : struct _FlxqE e;
317 1470 : e.a4 = a4; e.T = T; e.p = p; e.pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
318 1470 : return gerepileuptoint(av, gen_order(z, o, (void*)&e, &FlxqE_group));
319 : }
320 :
321 : GEN
322 49 : FlxqE_log(GEN a, GEN b, GEN o, GEN a4, GEN T, ulong p)
323 : {
324 49 : pari_sp av = avma;
325 : struct _FlxqE e;
326 49 : e.a4 = a4; e.T = T; e.p = p; e.pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
327 49 : return gerepileuptoint(av, gen_PH_log(a, b, o, (void*)&e, &FlxqE_group));
328 : }
329 :
330 : /***********************************************************************/
331 : /** Pairings **/
332 : /***********************************************************************/
333 : /* Derived from APIP by Jerome Milan, 2012 */
334 : static GEN
335 65236 : FlxqE_vert(GEN P, GEN Q, GEN a4, GEN T, ulong p, ulong pi)
336 : {
337 65236 : long vT = get_Flx_var(T);
338 : GEN df;
339 65236 : if (ell_is_inf(P)) return pol1_Flx(vT);
340 44082 : if (!Flx_equal(gel(Q,1), gel(P,1))) return Flx_sub(gel(Q,1), gel(P,1), p);
341 1846 : if (lgpol(gel(P,2))!=0) return pol1_Flx(vT);
342 798 : df = typ(a4)==t_VEC ? Flxq_mul_pre(gel(P,1), Flx_double(gel(a4,1), p), T,p,pi)
343 1246 : : a4;
344 1246 : return Flxq_inv_pre(Flx_add(Flx_triple(Flxq_sqr_pre(gel(P,1), T,p, pi), p),
345 : df, p), T, p, pi);
346 : }
347 :
348 : static GEN
349 26716 : FlxqE_Miller_line(GEN R, GEN Q, GEN slope, GEN a4, GEN T, ulong p, ulong pi)
350 : {
351 26716 : long vT = get_Flx_var(T);
352 26716 : GEN x = gel(Q,1), y = gel(Q,2);
353 26716 : GEN tmp1 = Flx_sub(x, gel(R,1), p);
354 26716 : GEN tmp2 = Flx_add(Flxq_mul_pre(tmp1, slope, T, p, pi), gel(R,2), p);
355 26716 : if (!Flx_equal(y, tmp2)) return Flx_sub(y, tmp2, p);
356 1230 : if (lgpol(y) == 0) return pol1_Flx(vT);
357 : else
358 : {
359 607 : GEN s1, s2, a2 = typ(a4)==t_VEC ? gel(a4,1): NULL;
360 607 : GEN y2i = Flxq_inv_pre(Flx_mulu(y, 2, p), T, p, pi);
361 607 : GEN df = a2 ? Flxq_mul_pre(x, Flx_mulu(a2, 2, p), T, p, pi): a4;
362 : GEN x3, ddf;
363 607 : s1 = Flxq_mul_pre(Flx_add(Flx_triple(Flxq_sqr_pre(x, T, p, pi), p), df, p), y2i, T, p, pi);
364 607 : if (!Flx_equal(s1, slope)) return Flx_sub(s1, slope, p);
365 244 : x3 = Flx_triple(x, p);
366 244 : ddf = a2 ? Flx_add(x3, a2, p): x3;
367 244 : s2 = Flxq_mul_pre(Flx_sub(ddf, Flxq_sqr_pre(s1, T,p,pi), p), y2i, T,p,pi);
368 244 : return lgpol(s2)!=0 ? s2: y2i;
369 : }
370 : }
371 :
372 : /* Computes the equation of the line tangent to R and returns its
373 : * evaluation at the point Q. Also doubles the point R. */
374 : static GEN
375 43768 : FlxqE_tangent_update(GEN R, GEN Q, GEN a4, GEN T, ulong p, ulong pi, GEN *pt_R)
376 : {
377 43768 : if (ell_is_inf(R))
378 : {
379 3942 : *pt_R = ellinf();
380 3942 : return pol1_Flx(get_Flx_var(T));
381 : }
382 39826 : else if (!lgpol(gel(R,2)))
383 : {
384 15959 : *pt_R = ellinf();
385 15959 : return FlxqE_vert(R, Q, a4, T, p, pi);
386 : } else {
387 : GEN slope;
388 23867 : *pt_R = FlxqE_dbl_slope(R, a4, T, p, pi, &slope);
389 23867 : return FlxqE_Miller_line(R, Q, slope, a4, T, p, pi);
390 : }
391 : }
392 :
393 : /* Computes the equation of the line through R and P, and returns its
394 : * evaluation at the point Q. Also adds P to the point R. */
395 : static GEN
396 4179 : FlxqE_chord_update(GEN R, GEN P, GEN Q, GEN a4, GEN T, ulong p, ulong pi, GEN *pt_R)
397 : {
398 4179 : if (ell_is_inf(R))
399 : {
400 77 : *pt_R = gcopy(P);
401 77 : return FlxqE_vert(P, Q, a4, T, p, pi);
402 : }
403 4102 : else if (ell_is_inf(P))
404 : {
405 0 : *pt_R = gcopy(R);
406 0 : return FlxqE_vert(R, Q, a4, T, p, pi);
407 : }
408 4102 : else if (Flx_equal(gel(P, 1), gel(R, 1)))
409 : {
410 1253 : if (Flx_equal(gel(P, 2), gel(R, 2)))
411 0 : return FlxqE_tangent_update(R, Q, a4, T, p, pi, pt_R);
412 : else
413 : {
414 1253 : *pt_R = ellinf();
415 1253 : return FlxqE_vert(R, Q, a4, T, p, pi);
416 : }
417 : } else {
418 : GEN slope;
419 2849 : *pt_R = FlxqE_add_slope(P, R, a4, T, p, pi, &slope);
420 2849 : return FlxqE_Miller_line(R, Q, slope, a4, T, p, pi);
421 : }
422 : }
423 :
424 : struct _FlxqE_miller
425 : {
426 : ulong p, pi;
427 : GEN T, a4, P;
428 : };
429 :
430 : static GEN
431 43768 : FlxqE_Miller_dbl(void* E, GEN d)
432 : {
433 43768 : struct _FlxqE_miller *m = (struct _FlxqE_miller *)E;
434 43768 : ulong p = m->p, pi = m->pi;
435 43768 : GEN T = m->T, a4 = m->a4, P = m->P;
436 43768 : GEN v, line, point = gel(d,3);
437 43768 : GEN N = Flxq_sqr_pre(gel(d,1), T, p, pi);
438 43768 : GEN D = Flxq_sqr_pre(gel(d,2), T, p, pi);
439 43768 : line = FlxqE_tangent_update(point, P, a4, T, p, pi, &point);
440 43768 : N = Flxq_mul_pre(N, line, T, p, pi);
441 43768 : v = FlxqE_vert(point, P, a4, T, p, pi);
442 43768 : D = Flxq_mul_pre(D, v, T, p, pi); return mkvec3(N, D, point);
443 : }
444 :
445 : static GEN
446 4179 : FlxqE_Miller_add(void* E, GEN va, GEN vb)
447 : {
448 4179 : struct _FlxqE_miller *m = (struct _FlxqE_miller *)E;
449 4179 : ulong p = m->p, pi = m->pi;
450 4179 : GEN T = m->T, a4 = m->a4, P = m->P;
451 : GEN v, line, point;
452 4179 : GEN na = gel(va,1), da = gel(va,2), pa = gel(va,3);
453 4179 : GEN nb = gel(vb,1), db = gel(vb,2), pb = gel(vb,3);
454 4179 : GEN N = Flxq_mul_pre(na, nb, T, p, pi);
455 4179 : GEN D = Flxq_mul_pre(da, db, T, p, pi);
456 4179 : line = FlxqE_chord_update(pa, pb, P, a4, T, p, pi, &point);
457 4179 : N = Flxq_mul_pre(N, line, T, p, pi);
458 4179 : v = FlxqE_vert(point, P, a4, T, p, pi);
459 4179 : D = Flxq_mul_pre(D, v, T, p, pi); return mkvec3(N, D, point);
460 : }
461 :
462 : /* Returns the Miller function f_{m, Q} evaluated at the point P using
463 : * the standard Miller algorithm. */
464 : static GEN
465 17135 : FlxqE_Miller(GEN Q, GEN P, GEN m, GEN a4, GEN T, ulong p, ulong pi)
466 : {
467 17135 : pari_sp av = avma;
468 : struct _FlxqE_miller d;
469 : GEN v, N, D, g1;
470 :
471 17135 : d.a4 = a4; d.T = T; d.p = p; d.P = P; d.pi = pi;
472 17135 : g1 = pol1_Flx(get_Flx_var(T));
473 17135 : v = gen_pow_i(mkvec3(g1,g1,Q), m, (void*)&d,
474 : FlxqE_Miller_dbl, FlxqE_Miller_add);
475 17135 : N = gel(v,1); D = gel(v,2);
476 17135 : return gerepileupto(av, Flxq_div_pre(N, D, T, p, pi));
477 : }
478 :
479 : GEN
480 14019 : FlxqE_weilpairing_pre(GEN P, GEN Q, GEN m, GEN a4, GEN T, ulong p, ulong pi)
481 : {
482 14019 : pari_sp av = avma;
483 : GEN N, D, w;
484 14019 : if (ell_is_inf(P) || ell_is_inf(Q)
485 11245 : || (Flx_equal(gel(P,1),gel(Q,1)) && Flx_equal(gel(P,2),gel(Q,2))))
486 5483 : return pol1_Flx(get_Flx_var(T));
487 8536 : N = FlxqE_Miller(P, Q, m, a4, T, p, pi);
488 8536 : D = FlxqE_Miller(Q, P, m, a4, T, p, pi);
489 8536 : w = Flxq_div_pre(N, D, T, p, pi); if (mpodd(m)) w = Flx_neg(w, p);
490 8536 : return gerepileupto(av, w);
491 : }
492 : GEN
493 21 : FlxqE_weilpairing(GEN P, GEN Q, GEN m, GEN a4, GEN T, ulong p)
494 21 : { return FlxqE_weilpairing_pre(P,Q,m,a4,T,p, SMALL_ULONG(p)?0:get_Fl_red(p)); }
495 :
496 : GEN
497 63 : FlxqE_tatepairing(GEN P, GEN Q, GEN m, GEN a4, GEN T, ulong p)
498 : {
499 63 : if (ell_is_inf(P) || ell_is_inf(Q)) return pol1_Flx(get_Flx_var(T));
500 63 : return FlxqE_Miller(P, Q, m, a4, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p));
501 : }
502 :
503 : static GEN
504 13998 : _FlxqE_pairorder(void *E, GEN P, GEN Q, GEN m, GEN F)
505 : {
506 13998 : struct _FlxqE *e = (struct _FlxqE *) E;
507 13998 : return Flxq_order(FlxqE_weilpairing_pre(P,Q,m,e->a4,e->T,e->p,e->pi), F, e->T, e->p);
508 : }
509 :
510 : GEN
511 22294 : Flxq_ellgroup(GEN a4, GEN a6, GEN N, GEN T, ulong p, GEN *pt_m)
512 : {
513 : struct _FlxqE e;
514 22294 : GEN q = powuu(p, get_Flx_degree(T));
515 22294 : e.a4=a4; e.a6=a6; e.T=T; e.p=p; e.pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
516 22294 : return gen_ellgroup(N, subiu(q,1), pt_m, (void*)&e, &FlxqE_group, _FlxqE_pairorder);
517 : }
518 :
519 : GEN
520 14370 : Flxq_ellgens(GEN a4, GEN a6, GEN ch, GEN D, GEN m, GEN T, ulong p)
521 : {
522 : GEN P;
523 14370 : pari_sp av = avma;
524 : struct _FlxqE e;
525 14370 : e.a4=a4; e.a6=a6; e.T=T; e.p=p; e.pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
526 14370 : switch(lg(D)-1)
527 : {
528 63 : case 0:
529 63 : return cgetg(1,t_VEC);
530 11801 : case 1:
531 11801 : P = gen_gener(gel(D,1), (void*)&e, &FlxqE_group);
532 11801 : P = mkvec(FlxqE_changepoint(P, ch, T, p));
533 11801 : break;
534 2506 : default:
535 2506 : P = gen_ellgens(gel(D,1), gel(D,2), m, (void*)&e, &FlxqE_group, _FlxqE_pairorder);
536 2506 : gel(P,1) = FlxqE_changepoint(gel(P,1), ch, T, p);
537 2506 : gel(P,2) = FlxqE_changepoint(gel(P,2), ch, T, p);
538 2506 : break;
539 : }
540 14307 : return gerepilecopy(av, P);
541 : }
542 : /***********************************************************************/
543 : /** Point counting **/
544 : /***********************************************************************/
545 :
546 : /* assume a and n are coprime */
547 : static GEN
548 76251 : RgX_circular_shallow(GEN P, long a, long n)
549 : {
550 76251 : long i, l = lgpol(P);
551 76251 : GEN Q = cgetg(2+n,t_POL);
552 76251 : Q[1] = P[1];
553 512330 : for(i=0; i<l; i++)
554 436079 : gel(Q,2+(i*a)%n) = gel(P,2+i);
555 168693 : for( ; i<n; i++)
556 92442 : gel(Q,2+(i*a)%n) = gen_0;
557 76251 : return normalizepol_lg(Q,2+n);
558 : }
559 :
560 : static GEN
561 76251 : ZpXQ_frob_cyc(GEN x, GEN T, GEN q, ulong p)
562 : {
563 76251 : long n = get_FpX_degree(T);
564 76251 : return FpX_rem(RgX_circular_shallow(x,p,n+1), T, q);
565 : }
566 :
567 : static GEN
568 113610 : ZpXQ_frob(GEN x, GEN Xm, GEN T, GEN q, ulong p)
569 : {
570 113610 : if (lg(Xm)==1)
571 43428 : return ZpXQ_frob_cyc(x, T, q, p);
572 : else
573 : {
574 70182 : long n = get_FpX_degree(T);
575 70182 : GEN V = RgX_blocks(RgX_inflate(x, p), n, p);
576 70182 : GEN W = ZXV_dotproduct(V, Xm);
577 70182 : return FpX_rem(W, T, q);
578 : }
579 : }
580 :
581 : struct _lift_lin
582 : {
583 : ulong p, pi;
584 : GEN sqx, Tp, ai, Xm;
585 : };
586 :
587 : static GEN
588 84077 : _lift_invl(void *E, GEN x)
589 : {
590 84077 : struct _lift_lin *d = (struct _lift_lin *) E;
591 84077 : GEN T = d->Tp;
592 84077 : ulong p = d->p, pi = d->pi;
593 84077 : GEN xai = Flxq_mul_pre(ZX_to_Flx(x, p), d->ai, T, p, pi);
594 84077 : return Flx_to_ZX(Flxq_lroot_fast_pre(xai, d->sqx, T, p, pi));
595 : }
596 : static GEN
597 23744 : _lift_lin(void *E, GEN F, GEN x2, GEN q)
598 : {
599 23744 : struct _lift_lin *d = (struct _lift_lin *) E;
600 23744 : pari_sp av = avma;
601 23744 : GEN T = gel(F,3), Xm = gel(F,4);
602 23744 : GEN y2 = ZpXQ_frob(x2, Xm, T, q, d->p);
603 23744 : GEN lin = FpX_add(ZX_mul(gel(F,1), y2), ZX_mul(gel(F,2), x2), q);
604 23744 : return gerepileupto(av, FpX_rem(lin, T, q));
605 : }
606 :
607 : static GEN
608 181041 : FpM_FpXV_bilinear(GEN P, GEN X, GEN Y, GEN p)
609 : {
610 181041 : pari_sp av = avma;
611 181041 : GEN s = ZX_mul(FpXV_FpC_mul(X,gel(P,1),p),gel(Y,1));
612 181041 : long i, l = lg(P);
613 851487 : for(i=2; i<l; i++)
614 670446 : s = ZX_add(s, ZX_mul(FpXV_FpC_mul(X,gel(P,i),p),gel(Y,i)));
615 181041 : return gerepileupto(av, FpX_red(s, p));
616 : }
617 :
618 : static GEN
619 181041 : FpM_FpXQV_bilinear(GEN P, GEN X, GEN Y, GEN T, GEN p)
620 181041 : { return FpX_rem(FpM_FpXV_bilinear(P,X,Y,p),T,p); }
621 :
622 : static GEN
623 120708 : FpXC_powderiv(GEN M, GEN p)
624 : {
625 : long i, l;
626 120708 : long v = varn(gel(M,2));
627 120708 : GEN m = cgetg_copy(M, &l);
628 120708 : gel(m,1) = pol_0(v);
629 120708 : gel(m,2) = pol_1(v);
630 447132 : for(i=2; i<l-1; i++)
631 326424 : gel(m,i+1) = FpX_Fp_mul(gel(M,i),utoi(i), p);
632 120708 : return m;
633 : }
634 :
635 : struct _lift_iso
636 : {
637 : GEN phi, Xm, T, sqx, Tp;
638 : ulong p, pi;
639 : };
640 :
641 : static GEN
642 60333 : _lift_iter(void *E, GEN x2, GEN q)
643 : {
644 60333 : struct _lift_iso *d = (struct _lift_iso *) E;
645 60333 : ulong p = d->p;
646 60333 : long n = lg(d->phi)-2;
647 60333 : GEN TN = FpXT_red(d->T, q), XN = FpXV_red(d->Xm, q);
648 60333 : GEN y2 = ZpXQ_frob(x2, XN, TN, q, p);
649 60333 : GEN xp = FpXQ_powers(x2, n, TN, q);
650 60333 : GEN yp = FpXQ_powers(y2, n, TN, q);
651 60333 : GEN V = FpM_FpXQV_bilinear(d->phi,xp,yp,TN,q);
652 60333 : return mkvec3(V,xp,yp);
653 : }
654 :
655 : static GEN
656 60333 : _lift_invd(void *E, GEN V, GEN v, GEN qM, long M)
657 : {
658 60333 : struct _lift_iso *d = (struct _lift_iso *) E;
659 : struct _lift_lin e;
660 60333 : ulong p = d->p, pi = d->pi;
661 60333 : GEN TM = FpXT_red(d->T, qM), XM = FpXV_red(d->Xm, qM);
662 60333 : GEN xp = FpXV_red(gel(v,2), qM);
663 60333 : GEN yp = FpXV_red(gel(v,3), qM);
664 60333 : GEN Dx = FpM_FpXQV_bilinear(d->phi, FpXC_powderiv(xp, qM), yp, TM, qM);
665 60333 : GEN Dy = FpM_FpXQV_bilinear(d->phi, xp, FpXC_powderiv(yp, qM), TM, qM);
666 60333 : GEN F = mkvec4(Dy, Dx, TM, XM);
667 60333 : e.ai = Flxq_inv_pre(ZX_to_Flx(Dy,p),d->Tp, p, pi);
668 60333 : e.sqx = d->sqx; e.Tp = d->Tp; e.p=p; e.pi=pi; e.Xm = XM;
669 60333 : return gen_ZpX_Dixon(F,V,qM,utoipos(p),M,(void*) &e, _lift_lin, _lift_invl);
670 : }
671 :
672 : static GEN
673 25053 : lift_isogeny(GEN phi, GEN x0, long n, GEN Xm, GEN T, GEN sqx, GEN Tp,
674 : ulong p, ulong pi)
675 : {
676 : struct _lift_iso d;
677 25053 : d.phi = phi; d.Xm = Xm; d.T = T;
678 25053 : d.sqx = sqx; d.Tp = Tp; d.p = p; d.pi = pi;
679 25053 : return gen_ZpX_Newton(x0, utoipos(p), n,(void*)&d, _lift_iter, _lift_invd);
680 : }
681 :
682 : static GEN
683 25032 : getc2(GEN act, GEN X, GEN T, GEN q, ulong p, long N)
684 : {
685 25032 : GEN A1 = RgV_to_RgX(gel(act,1),0), A2 = RgV_to_RgX(gel(act,2),0);
686 25032 : long n = brent_kung_optpow(maxss(degpol(A1),degpol(A2)),2,1);
687 25032 : GEN xp = FpXQ_powers(X,n,T,q);
688 25032 : GEN P = FpX_FpXQV_eval(A1, xp, T, q);
689 25032 : GEN Q = FpX_FpXQV_eval(A2, xp, T, q);
690 25032 : return ZpXQ_div(P, Q, T, q, utoipos(p), N);
691 : }
692 :
693 : struct _ZpXQ_norm
694 : {
695 : long n;
696 : GEN T, p;
697 : };
698 :
699 : static GEN
700 32823 : ZpXQ_norm_mul(void *E, GEN x, GEN y)
701 : {
702 32823 : struct _ZpXQ_norm *D = (struct _ZpXQ_norm*)E;
703 32823 : GEN P = gel(x,1), Q = gel(y,1);
704 32823 : long a = mael(x,2,1), b = mael(y,2,1);
705 32823 : retmkvec2(FpXQ_mul(P,ZpXQ_frob_cyc(Q, D->T, D->p, a), D->T, D->p),
706 : mkvecsmall((a*b)%D->n));
707 : }
708 : static GEN
709 22715 : ZpXQ_norm_sqr(void *E, GEN x) { return ZpXQ_norm_mul(E, x, x); }
710 :
711 : /* Assume T = Phi_(n) and n prime */
712 : GEN
713 11340 : ZpXQ_norm_pcyc(GEN x, GEN T, GEN q, GEN p)
714 : {
715 : GEN z;
716 : struct _ZpXQ_norm D;
717 11340 : long d = get_FpX_degree(T);
718 11340 : D.T = T; D.p = q; D.n = d+1;
719 11340 : if (d==1) return ZX_copy(x);
720 11340 : z = mkvec2(x,mkvecsmall(p[2]));
721 11340 : z = gen_powu_i(z,d,(void*)&D,ZpXQ_norm_sqr,ZpXQ_norm_mul);
722 11340 : return gmael(z,1,2);
723 : }
724 :
725 : /* Assume T = Phi_(n) and n prime */
726 : static GEN
727 11102 : ZpXQ_sqrtnorm_pcyc(GEN x, GEN T, GEN q, GEN p, long e)
728 : {
729 11102 : GEN z = ZpXQ_norm_pcyc(x, T, q, p);
730 11102 : return Zp_sqrtlift(z,Fp_sqrt(z,p),p,e);
731 : }
732 :
733 : /* Assume a = 1 [p], return the square root of the norm */
734 : static GEN
735 13951 : ZpXQ_sqrtnorm(GEN a, GEN T, GEN q, GEN p, long e)
736 : {
737 13951 : GEN s = Fp_div(FpXQ_trace(ZpXQ_log(a, T, p, e), T, q), gen_2, q);
738 13951 : return modii(gel(Qp_exp(cvtop(s, p, e-1)),4), q);
739 : }
740 :
741 : struct _teich_lin
742 : {
743 : ulong p, pi;
744 : GEN sqx, Tp;
745 : long m;
746 : };
747 :
748 : static GEN
749 29512 : _teich_invl(void *E, GEN x)
750 : {
751 29512 : struct _teich_lin *d = (struct _teich_lin *) E;
752 29512 : ulong p = d->p, pi = d->pi;
753 29512 : return Flx_to_ZX(Flxq_lroot_fast_pre(ZX_to_Flx(x,p), d->sqx, d->Tp, p, pi));
754 : }
755 :
756 : static GEN
757 8953 : _teich_lin(void *E, GEN F, GEN x2, GEN q)
758 : {
759 8953 : struct _teich_lin *d = (struct _teich_lin *) E;
760 8953 : pari_sp av = avma;
761 8953 : GEN T = gel(F,2), Xm = gel(F,3);
762 8953 : GEN y2 = ZpXQ_frob(x2, Xm, T, q, d->p);
763 8953 : GEN lin = FpX_sub(y2, ZX_mulu(ZX_mul(gel(F,1), x2), d->p), q);
764 8953 : return gerepileupto(av, FpX_rem(lin, T, q));
765 : }
766 :
767 : struct _teich_iso
768 : {
769 : GEN Xm, T, sqx, Tp;
770 : ulong p, pi;
771 : };
772 :
773 : static GEN
774 20559 : _teich_iter(void *E, GEN x2, GEN q)
775 : {
776 20559 : struct _teich_iso *d = (struct _teich_iso *) E;
777 20559 : ulong p = d->p;
778 20559 : GEN TN = FpXT_red(d->T, q), XN = FpXV_red(d->Xm, q);
779 20559 : GEN y2 = ZpXQ_frob(x2, XN, TN, q, d->p);
780 20559 : GEN x1 = FpXQ_powu(x2, p-1, TN, q);
781 20559 : GEN xp = FpXQ_mul(x2, x1, TN, q);
782 20559 : GEN V = FpX_sub(y2,xp,q);
783 20559 : return mkvec2(V,x1);
784 : }
785 :
786 : static GEN
787 20559 : _teich_invd(void *E, GEN V, GEN v, GEN qM, long M)
788 : {
789 20559 : struct _teich_iso *d = (struct _teich_iso *) E;
790 : struct _teich_lin e;
791 20559 : ulong p = d->p;
792 20559 : GEN TM = FpXT_red(d->T, qM), XM = FpXV_red(d->Xm, qM);
793 20559 : GEN x1 = FpX_red(gel(v,2), qM);
794 20559 : GEN F = mkvec3(x1, TM, XM);
795 20559 : e.sqx = d->sqx; e.Tp = d->Tp; e.p = p; e.pi = d->pi;
796 20559 : return gen_ZpX_Dixon(F,V,qM,utoipos(p),M,(void*) &e, _teich_lin, _teich_invl);
797 : }
798 :
799 : static GEN
800 10234 : Teichmuller_lift(GEN x, GEN Xm, GEN T, GEN sqx, GEN Tp, ulong p, ulong pi,
801 : long N)
802 : {
803 : struct _teich_iso d;
804 10234 : d.Xm = Xm; d.T = T; d.sqx = sqx; d.Tp = Tp; d.p = p; d.pi = pi;
805 10234 : return gen_ZpX_Newton(x,utoipos(p), N,(void*)&d, _teich_iter, _teich_invd);
806 : }
807 :
808 : static GEN
809 25053 : get_norm(GEN a4, GEN a6, GEN T, ulong p, ulong pi, long N)
810 : {
811 25053 : long sv=T[1];
812 : GEN a;
813 25053 : if (p==3) a = gel(a4,1);
814 : else
815 : {
816 10248 : GEN P = mkpoln(4, pol1_Flx(sv), pol0_Flx(sv), a4, a6);
817 10248 : a = gel(FlxqX_powu_pre(P, p>>1, T,p,pi), 2+p-1);
818 : }
819 25053 : return Zp_sqrtnlift(gen_1,subss(p,1),utoi(Flxq_norm(a,T,p)),utoipos(p), N);
820 : }
821 :
822 : static GEN
823 25032 : fill_pols(long n, const long *v, long m, const long *vn,
824 : const long *vd, GEN *act)
825 : {
826 : long i, j;
827 25032 : long d = upowuu(n,12/(n-1));
828 25032 : GEN N, D, M = zeromatcopy(n+1,n+1);
829 25032 : gmael(M,1,n+1) = gen_1;
830 120764 : for (i = 2; i <= n+1; i++)
831 339486 : for (j = i-1; j <= n; j++) gmael(M,i,j) = mulis(powuu(d,i-2), v[j-i+1]);
832 25032 : N = cgetg(m+1,t_COL);
833 25032 : D = cgetg(m+1,t_COL);
834 135541 : for(i = 1; i <= m; i++)
835 : {
836 110509 : gel(N,i) = stoi(*vn++);
837 110509 : gel(D,i) = stoi(*vd++);
838 : }
839 25032 : *act = mkmat2(N,D); return M;
840 : }
841 :
842 : /*
843 : These polynomials were extracted from the ECHIDNA databases
844 : available at <http://echidna.maths.usyd.edu.au/echidna/>
845 : and computed by David R. Kohel.
846 : Return the matrix of the modular polynomial, set act to the parametrization,
847 : and set dj to the opposite of the supersingular j-invariant.
848 : */
849 : static GEN
850 25032 : get_Kohel_polynomials(ulong p, GEN *act, long *dj)
851 : {
852 25032 : const long mat3[] = {-1,-36,-270};
853 25032 : const long num3[] = {1,-483,-21141,-59049};
854 25032 : const long den3[] = {1,261, 4347, -6561};
855 25032 : const long mat5[] = {-1,-30,-315,-1300,-1575};
856 25032 : const long num5[] = {-1,490,20620,158750,78125};
857 25032 : const long den5[] = {-1,-254,-4124,-12250,3125};
858 25032 : const long mat7[] = {-1,-28,-322,-1904,-5915,-8624,-4018};
859 25032 : const long num7[] = {1,-485,-24058,-343833,-2021642,-4353013,-823543};
860 25032 : const long den7[] = {1,259,5894,49119,168406,166355,-16807};
861 25032 : const long mat13[]= {-1,-26,-325,-2548,-13832,-54340,-157118,-333580,-509366,
862 : -534820,-354536,-124852,-15145};
863 25032 : const long num13[]= {1,-487,-24056,-391463,-3396483,-18047328,-61622301,
864 : -133245853,-168395656,-95422301,-4826809};
865 25032 : const long den13[]= {1,257,5896,60649,364629,1388256,3396483,5089019,4065464,
866 : 1069939,-28561};
867 25032 : switch(p)
868 : {
869 14805 : case 3:
870 14805 : *dj = 0;
871 14805 : return fill_pols(3,mat3,4,num3,den3,act);
872 10178 : case 5:
873 10178 : *dj = 0;
874 10178 : return fill_pols(5,mat5,5,num5,den5,act);
875 35 : case 7:
876 35 : *dj = 1;
877 35 : return fill_pols(7,mat7,7,num7,den7,act);
878 14 : case 13:
879 14 : *dj = 8;
880 14 : return fill_pols(13,mat13,11,num13,den13,act);
881 : }
882 : *dj=0; *act = NULL; return NULL; /* LCOV_EXCL_LINE */
883 : }
884 :
885 : long
886 32494 : zx_is_pcyc(GEN T)
887 : {
888 32494 : long i, n = degpol(T);
889 32494 : if (!uisprime(n+1)) return 0;
890 99148 : for (i = 0; i <= n; i++)
891 87808 : if (T[i+2] != 1UL) return 0;
892 11340 : return 1;
893 : }
894 :
895 : static GEN
896 25032 : Flxq_ellcard_Kohel(GEN a4, GEN a6, GEN T, ulong p)
897 : {
898 25032 : pari_sp av = avma, av2;
899 : pari_timer ti;
900 25032 : long n = get_Flx_degree(T), N = (n+4)/2, dj;
901 25032 : GEN q = powuu(p, N);
902 : GEN T2, Xm, s1, c2, t, lr, S1, sqx, Nc2, Np;
903 25032 : GEN act, phi = get_Kohel_polynomials(p, &act, &dj);
904 25032 : long ispcyc = zx_is_pcyc(get_Flx_mod(T));
905 25032 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
906 25032 : timer_start(&ti);
907 25032 : if (!ispcyc)
908 : {
909 13937 : T2 = Flx_Teichmuller(get_Flx_mod(T),p,N);
910 13937 : if (DEBUGLEVEL) timer_printf(&ti,"Teich");
911 : } else
912 11095 : T2 = Flx_to_ZX(get_Flx_mod(T));
913 :
914 25032 : T2 = FpX_get_red(T2, q); T = ZXT_to_FlxT(T2, p);
915 25032 : av2 = avma;
916 25032 : if (DEBUGLEVEL) timer_printf(&ti,"Barrett");
917 25032 : if (!ispcyc)
918 : {
919 13937 : Xm = FpXQ_powers(pol_xn(n,get_FpX_var(T2)),p-1,T2,q);
920 13937 : if (DEBUGLEVEL) timer_printf(&ti,"Xm");
921 : } else
922 11095 : Xm = cgetg(1,t_VEC);
923 25032 : s1 = Flxq_inv_pre(Flx_Fl_add(Flxq_ellj(a4,a6,T,p),dj, p),T,p,pi);
924 25032 : lr = Flxq_lroot_pre(polx_Flx(get_Flx_var(T)), T,p,pi);
925 25032 : sqx = Flxq_powers_pre(lr, p-1, T, p, pi);
926 25032 : S1 = lift_isogeny(phi, Flx_to_ZX(s1), N, Xm, T2, sqx, T,p,pi);
927 25032 : if (DEBUGLEVEL) timer_printf(&ti,"Lift isogeny");
928 25032 : c2 = getc2(act, S1, T2, q, p, N);
929 25032 : if (DEBUGLEVEL) timer_printf(&ti,"c^2");
930 25032 : if (p>3 && !ispcyc)
931 : {
932 10220 : GEN c2p = Flx_to_ZX(Flxq_inv_pre(ZX_to_Flx(c2,p),T,p,pi));
933 10220 : GEN tc2 = Teichmuller_lift(c2p,Xm, T2,sqx,T,p,pi,N);
934 10220 : if (DEBUGLEVEL) timer_printf(&ti,"Teichmuller/Fq");
935 10220 : c2 = FpX_rem(FpX_mul(tc2,c2,q),T2,q);
936 : }
937 25032 : c2 = gerepileupto(av2, c2);
938 25032 : if (DEBUGLEVEL) timer_printf(&ti,"tc2");
939 25032 : Nc2 = (ispcyc? ZpXQ_sqrtnorm_pcyc: ZpXQ_sqrtnorm)(c2, T2, q, utoipos(p), N);
940 25032 : if (DEBUGLEVEL) timer_printf(&ti,"Norm");
941 25032 : Np = get_norm(a4,a6,T,p,pi,N);
942 25032 : if (p>3 && ispcyc)
943 : {
944 7 : GEN Ncpi = utoi(Fl_inv(umodiu(Nc2,p), p));
945 7 : GEN tNc2 = Zp_sqrtnlift(gen_1, subss(p,1), Ncpi, utoipos(p),N);
946 7 : if (DEBUGLEVEL) timer_printf(&ti,"Teichmuller/Fp");
947 7 : Nc2 = Fp_mul(Nc2,tNc2,q);
948 : }
949 25032 : t = Fp_center_i(Fp_mul(Nc2,Np,q),q,shifti(q,-1));
950 25032 : return gerepileupto(av, subii(addiu(powuu(p,n),1),t));
951 : }
952 :
953 : /* Use Damien Robert's method */
954 : static GEN
955 21 : get_trace_Robert(GEN J, GEN phi, GEN Xm, GEN T, GEN q, ulong p, long e)
956 : {
957 21 : long n = lg(phi)-2;
958 21 : GEN K = ZpXQ_frob(J, Xm, T, q, p);
959 21 : GEN Jp = FpXQ_powers(J, n, T, q);
960 21 : GEN Kp = FpXQ_powers(K, n, T, q);
961 21 : GEN Jd = FpXC_powderiv(Jp, q);
962 21 : GEN Kd = FpXC_powderiv(Kp, q);
963 21 : GEN Dx = FpM_FpXQV_bilinear(phi, Kd, Jp, T, q);
964 21 : GEN Dy = FpM_FpXQV_bilinear(phi, Kp, Jd, T, q);
965 21 : GEN C = ZpXQ_inv(ZX_divuexact(Dy, p), T, utoi(p), e);
966 21 : return FpX_neg(FpXQ_mul(Dx, C, T, q), q);
967 : }
968 :
969 : /* in p^2, so p is tiny */
970 : static GEN
971 21 : Flxq_ellcard_Harley(GEN a4, GEN a6, GEN T, ulong p)
972 : {
973 21 : pari_sp av = avma, av2;
974 : pari_timer ti;
975 21 : long n = get_Flx_degree(T), N = (n+5)/2;
976 21 : GEN pp = utoipos(p), q = powuu(p, N);
977 : GEN T2, j, t, phi, J1, sqx, Xm, c2, tc2, c2p, Nc2, Np;
978 21 : long ispcyc = zx_is_pcyc(get_Flx_mod(T));
979 21 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p); /* = 0 here */
980 21 : timer_start(&ti);
981 21 : if (!ispcyc)
982 : {
983 14 : T2 = Flx_Teichmuller(get_Flx_mod(T),p,N);
984 14 : if (DEBUGLEVEL) timer_printf(&ti,"Teich");
985 : } else
986 7 : T2 = Flx_to_ZX(get_Flx_mod(T));
987 21 : T2 = FpX_get_red(T2, q); T = ZXT_to_FlxT(T2, p);
988 21 : av2 = avma;
989 21 : if (DEBUGLEVEL) timer_printf(&ti,"Barrett");
990 21 : if (!ispcyc)
991 : {
992 14 : Xm = FpXQ_powers(pol_xn(n,get_FpX_var(T2)),p-1,T2,q);
993 14 : if (DEBUGLEVEL) timer_printf(&ti,"Xm");
994 : } else
995 7 : Xm = cgetg(1,t_VEC);
996 21 : j = Flxq_ellj(a4,a6,T,p);
997 21 : sqx = Flxq_powers_pre(Flxq_lroot_pre(polx_Flx(T[1]), T,p,pi), p-1, T,p,pi);
998 21 : phi = polmodular_ZM(p, 0);
999 21 : if (DEBUGLEVEL) timer_printf(&ti,"phi");
1000 21 : J1 = lift_isogeny(phi, Flx_to_ZX(j), N, Xm, T2,sqx,T,p,pi);
1001 21 : if (DEBUGLEVEL) timer_printf(&ti,"Lift isogeny");
1002 21 : c2 = get_trace_Robert(J1, phi, Xm, T2, q, p, N);
1003 21 : q = diviuexact(q,p); N--;
1004 21 : if (DEBUGLEVEL) timer_printf(&ti,"c^2");
1005 21 : if (!ispcyc)
1006 : {
1007 14 : c2p = Flx_to_ZX(Flxq_inv_pre(ZX_to_Flx(c2,p),T,p,pi));
1008 14 : tc2 = Teichmuller_lift(c2p,Xm, T2,sqx,T,p,pi,N);
1009 14 : if (DEBUGLEVEL) timer_printf(&ti,"teichmuller");
1010 14 : c2 = FpX_rem(FpX_mul(tc2,c2,q),T2,q);
1011 : }
1012 21 : c2 = gerepileupto(av2, c2);
1013 21 : q = powuu(p, N);
1014 21 : Nc2 = (ispcyc? ZpXQ_sqrtnorm_pcyc: ZpXQ_sqrtnorm)(c2, T2, q, pp, N);
1015 21 : if (DEBUGLEVEL) timer_printf(&ti,"Norm");
1016 21 : Np = get_norm(a4,a6,T,p,pi,N);
1017 21 : if (ispcyc)
1018 : {
1019 7 : GEN Ncpi = utoi(Fl_inv(umodiu(Nc2,p), p));
1020 7 : GEN tNc2 = Zp_sqrtnlift(gen_1, subss(p,1), Ncpi, pp, N);
1021 7 : if (DEBUGLEVEL) timer_printf(&ti,"Teichmuller/Fp");
1022 7 : Nc2 = Fp_mul(Nc2,tNc2,q);
1023 : }
1024 21 : t = Fp_center_i(Fp_mul(Nc2,Np,q),q,shifti(q,-1));
1025 21 : return gerepileupto(av, subii(addiu(powuu(p,n),1),t));
1026 : }
1027 :
1028 : /***************************************************************************/
1029 : /* Shanks-Mestre */
1030 : /***************************************************************************/
1031 :
1032 : /* Return the lift of a (mod b), which is closest to h */
1033 : static GEN
1034 10346 : closest_lift(GEN a, GEN b, GEN h)
1035 10346 : { return addii(a, mulii(b, diviiround(subii(h,a), b))); }
1036 :
1037 : /* find multiple of order of f using Baby Step/Giant Step, f^h close to 1,
1038 : * order lies in an interval of size <= 'bound' and known mod B */
1039 : static GEN
1040 5572 : _FlxqE_order_multiple(void *E, GEN f, GEN h, GEN bound, GEN B)
1041 : {
1042 5572 : pari_sp av = avma, av1;
1043 : pari_timer Ti;
1044 5572 : long i, s = ceilsqrtdiv(bound, B) >> 1;
1045 : GEN P, F, tx, ti, fg, fh;
1046 :
1047 5572 : P = fh = _FlxqE_mul(E, f, h);
1048 5572 : if (DEBUGLEVEL >= 6) timer_start(&Ti);
1049 5572 : if (ell_is_inf(fh)) return h;
1050 5146 : F = _FlxqE_mul(E, f, B);
1051 5146 : if (s < 3)
1052 : { /* we're nearly done: naive search */
1053 1093 : GEN Q = P;
1054 1093 : for (i=1;; i++)
1055 : {
1056 3633 : P = _FlxqE_add(E, P, F); /* h.f + i.F */
1057 3633 : if (ell_is_inf(P)) return gerepileupto(av, addii(h, mului(i,B)));
1058 3080 : Q = _FlxqE_sub(E, Q, F); /* h.f - i.F */
1059 3080 : if (ell_is_inf(Q)) return gerepileupto(av, subii(h, mului(i,B)));
1060 : }
1061 : }
1062 4053 : tx = cgetg(s+1,t_VECSMALL); av1 = avma;
1063 43295 : for (i=1; i<=s; i++)
1064 : { /* baby steps */
1065 39573 : tx[i] = hash_GEN(gel(P, 1));
1066 39573 : P = _FlxqE_add(E, P, F); /* h.f + i.F */
1067 39573 : if (ell_is_inf(P)) return gerepileupto(av, addii(h, mului(i,B)));
1068 39242 : if (gc_needed(av1,3))
1069 : {
1070 0 : if(DEBUGMEM>1) pari_warn(warnmem,"[Flxq_ellcard] baby steps, i=%ld",i);
1071 0 : P = gerepileupto(av1,P);
1072 : }
1073 : }
1074 3722 : if (DEBUGLEVEL >= 6) timer_printf(&Ti,"[Flxq_ellcard] baby steps, s = %ld",s);
1075 : /* giant steps: fg = s.F */
1076 3722 : fg = gerepileupto(av1, _FlxqE_sub(E, P, fh));
1077 3722 : if (ell_is_inf(fg)) return gerepileupto(av, mului(s,B));
1078 3722 : ti = vecsmall_indexsort(tx); /* = permutation sorting tx */
1079 3722 : tx = perm_mul(tx,ti);
1080 3722 : if (DEBUGLEVEL >= 6) timer_printf(&Ti, "[Flxq_ellcard] sorting");
1081 3722 : av1 = avma;
1082 3722 : for (P=fg, i=1; ; i++)
1083 35528 : {
1084 39250 : long k = hash_GEN(gel(P,1)), r = zv_search(tx, k);
1085 39250 : if (r)
1086 : {
1087 7448 : while (r && tx[r] == k) r--;
1088 3722 : for (r++; r <= s && tx[r] == k; r++)
1089 : {
1090 3722 : long j = ti[r]-1;
1091 3722 : GEN Q = _FlxqE_add(E, _FlxqE_mul(E, F, stoi(j)), fh);
1092 3722 : if (DEBUGLEVEL >= 6)
1093 0 : timer_printf(&Ti, "[Flxq_ellcard] giant steps, i = %ld",i);
1094 3722 : if (Flx_equal(gel(P,1), gel(Q,1)))
1095 : {
1096 3722 : if (Flx_equal(gel(P,2), gel(Q,2))) i = -i;
1097 3722 : return gerepileupto(av,addii(h, mulii(addis(mulss(s,i), j), B)));
1098 : }
1099 : }
1100 : }
1101 35528 : P = _FlxqE_add(E, P, fg);
1102 35528 : if (gc_needed(av1,3))
1103 : {
1104 0 : if(DEBUGMEM>1) pari_warn(warnmem,"[Flxq_ellcard] giants steps, i=%ld",i);
1105 0 : P = gerepileupto(av1,P);
1106 : }
1107 : }
1108 : }
1109 : static GEN
1110 5572 : _FlxqE_order(void *E, GEN f, GEN h, GEN bound, GEN B)
1111 : {
1112 5572 : GEN o = _FlxqE_order_multiple(E, f, h, bound, B);
1113 5572 : return gen_order(f, o, E, &FlxqE_group);
1114 : }
1115 :
1116 : static void
1117 64190 : Flx_next(GEN t, ulong p)
1118 : {
1119 : long i;
1120 78141 : for(i=2;;i++)
1121 78141 : if (uel(t,i)==p-1) t[i]=0; else { t[i]++; break; }
1122 64190 : }
1123 :
1124 : static void
1125 64190 : Flx_renormalize_ip(GEN x, long lx)
1126 : {
1127 : long i;
1128 78141 : for (i = lx-1; i>=2; i--)
1129 71057 : if (x[i]) break;
1130 64190 : setlg(x, i+1);
1131 64190 : }
1132 :
1133 : static ulong
1134 6188 : F3xq_ellcard_naive(GEN a2, GEN a6, GEN T)
1135 : {
1136 6188 : pari_sp av = avma;
1137 6188 : long i, d = get_Flx_degree(T), lx = d+2;
1138 6188 : long q = upowuu(3, d), a;
1139 6188 : GEN x = zero_zv(lx); x[1] = get_Flx_var(T);
1140 26978 : for(a=1, i=0; i<q; i++)
1141 : {
1142 : GEN rhs;
1143 20790 : Flx_renormalize_ip(x, lx);
1144 20790 : rhs = Flx_add(Flxq_mul(Flxq_sqr(x, T, 3), Flx_add(x, a2, 3), T, 3), a6, 3);
1145 20790 : if (!lgpol(rhs)) a++; else if (Flxq_issquare(rhs, T, 3)) a+=2;
1146 20790 : Flx_next(x, 3);
1147 : }
1148 6188 : set_avma(av); return a;
1149 : }
1150 :
1151 : /* p^deg(T) is tiny */
1152 : static ulong
1153 896 : Flxq_ellcard_naive(GEN a4, GEN a6, GEN T, ulong p)
1154 : {
1155 896 : pari_sp av = avma;
1156 896 : long i, d = get_Flx_degree(T), lx = d+2;
1157 896 : long q = upowuu(p, d), a;
1158 896 : GEN x = zero_zv(lx); x[1] = get_Flx_var(T);
1159 44296 : for(a = 1, i = 0; i < q; i++)
1160 : {
1161 : GEN x2, rhs;
1162 43400 : Flx_renormalize_ip(x, lx);
1163 43400 : x2 = Flxq_sqr_pre(x, T, p, 0);
1164 43400 : rhs = Flx_add(Flxq_mul_pre(x, Flx_add(x2, a4, p), T, p, 0), a6, p);
1165 43400 : if (!lgpol(rhs)) a++; else if (Flxq_issquare(rhs,T,p)) a += 2;
1166 43400 : Flx_next(x,p);
1167 : }
1168 896 : set_avma(av); return a;
1169 : }
1170 :
1171 : static long
1172 11553 : Flxq_kronecker(GEN x, GEN T, ulong p)
1173 : {
1174 : pari_sp av;
1175 11553 : if (lgpol(x) == 0) return 0;
1176 11533 : av = avma; return gc_long(av, krouu(Flxq_norm(x, T, p), p));
1177 : }
1178 :
1179 : /* Find x such that kronecker(u = x^3+a4x+a6, p) is KRO.
1180 : * Return point [x*u,u^2] on E (KRO=1) / E^twist (KRO=-1) */
1181 : static GEN
1182 5572 : Flxq_ellpoint(long KRO, GEN a4, GEN a6, GEN T, ulong p, ulong pi)
1183 : {
1184 5572 : long v = get_Flx_var(T), n = get_Flx_degree(T);
1185 : for(;;)
1186 5981 : {
1187 11553 : GEN x = random_Flx(n, v, p), x2 = Flxq_sqr_pre(x,T,p,pi);
1188 11553 : GEN u = Flx_add(a6, Flxq_mul_pre(Flx_add(a4, x2, p), x, T,p, pi), p);
1189 11553 : if (Flxq_kronecker(u,T,p) == KRO)
1190 5572 : return mkvec2(Flxq_mul_pre(u,x, T,p,pi), Flxq_sqr_pre(u, T,p,pi));
1191 : }
1192 : }
1193 :
1194 : static GEN
1195 4774 : Flxq_ellcard_Shanks(GEN a4, GEN a6, GEN q, GEN T, ulong p)
1196 : {
1197 4774 : pari_sp av = avma;
1198 4774 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
1199 4774 : long v = get_Flx_var(T), KRO = -1;
1200 : GEN h,f, A, B;
1201 4774 : GEN q1p = addiu(q,1), q2p = shifti(q1p, 1);
1202 4774 : GEN bound = addiu(sqrti(gmul2n(q,4)), 1); /* ceil( 4sqrt(q) ) */
1203 : struct _FlxqE e;
1204 4774 : e.p = p; e.pi = pi; e.T = Flx_get_red_pre(T, p, pi);
1205 : /* once #E(Flxq) is known mod B >= bound, it is determined */
1206 4774 : switch(FlxqX_nbroots(mkpoln(4, pol1_Flx(v), pol0_Flx(v), a4, a6), T, p))
1207 : { /* how many 2-torsion points ? */
1208 2891 : case 3: A = gen_0; B = utoipos(4); break;
1209 1414 : case 1: A = gen_0; B = gen_2; break;
1210 469 : default: A = gen_1; B = gen_2; break; /* 0 */
1211 : }
1212 : for(;;)
1213 : {
1214 5572 : h = closest_lift(A, B, q1p);
1215 : /* [ux, u^2] is on E_u: y^2 = x^3 + c4 u^2 x + c6 u^3
1216 : * E_u isomorphic to E (resp. E') iff KRO = 1 (resp. -1)
1217 : * #E(F_p) = p+1 - a_p, #E'(F_p) = p+1 + a_p
1218 : *
1219 : * #E_u(Flxq) = A (mod B), h is close to #E_u(Flxq) */
1220 5572 : KRO = -KRO;
1221 5572 : f = Flxq_ellpoint(KRO, a4,a6, T,p,pi);
1222 5572 : e.a4 = Flxq_mul_pre(a4, gel(f,2), T,p,pi); /* a4 for E_u */
1223 5572 : h = _FlxqE_order((void*)&e, f, h, bound, B);
1224 : /* h | #E_u(Flxq) = A (mod B) */
1225 5572 : A = Z_chinese_all(A, gen_0, B, h, &B);
1226 5572 : if (cmpii(B, bound) >= 0) break;
1227 : /* not done, update A mod B for the _next_ curve, isomorphic to
1228 : * the quadratic twist of this one */
1229 798 : A = remii(subii(q2p,A), B); /* #E(Fq)+#E'(Fq) = 2q+2 */
1230 : }
1231 4774 : h = closest_lift(A, B, q1p);
1232 4774 : return gerepileuptoint(av, KRO == 1? h: subii(q2p,h));
1233 : }
1234 :
1235 : static GEN
1236 20993 : F3xq_ellcard(GEN a2, GEN a6, GEN T)
1237 : {
1238 20993 : long n = get_Flx_degree(T);
1239 20993 : if (n <= 2)
1240 5887 : return utoi(F3xq_ellcard_naive(a2, a6, T));
1241 : else
1242 : {
1243 15106 : GEN q1 = addiu(powuu(3, get_Flx_degree(T)), 1), t;
1244 15106 : GEN a = Flxq_div(a6,Flxq_powu(a2,3,T,3),T,3);
1245 15106 : if (Flx_equal1(Flxq_powu(a, 8, T, 3)))
1246 : {
1247 301 : GEN P = Flxq_minpoly(a,T,3);
1248 301 : long dP = degpol(P); /* dP <= 2 */
1249 301 : ulong q = upowuu(3,dP);
1250 301 : GEN A2 = pol1_Flx(P[1]), A6 = Flx_rem(polx_Flx(P[1]), P, 3);
1251 301 : long tP = q + 1 - F3xq_ellcard_naive(A2, A6, P);
1252 301 : t = elltrace_extension(stoi(tP), n/dP, utoi(q));
1253 301 : if (umodiu(t, 3)!=1) t = negi(t);
1254 301 : return Flx_equal1(a2) || Flxq_issquare(a2,T,3) ? subii(q1,t): addii(q1,t);
1255 : }
1256 14805 : else return Flxq_ellcard_Kohel(mkvec(a2), a6, T, 3);
1257 : }
1258 : }
1259 :
1260 : static GEN
1261 11144 : Flxq_ellcard_Satoh(GEN a4, GEN a6, GEN j, GEN T, ulong p)
1262 : {
1263 11144 : long n = get_Flx_degree(T);
1264 11144 : if (n <= 2)
1265 896 : return utoi(Flxq_ellcard_naive(a4, a6, T, p));
1266 : else
1267 : {
1268 10248 : GEN jp = Flxq_powu(j, p, T, p);
1269 10248 : GEN s = Flx_add(j, jp, p);
1270 10248 : if (degpol(s) <= 0)
1271 : { /* it is assumed j not in F_p */
1272 0 : GEN m = Flxq_mul(j, jp, T, p);
1273 0 : if (degpol(m) <= 0)
1274 : {
1275 0 : GEN q = sqru(p);
1276 0 : GEN q1 = addiu(powuu(p, get_Flx_degree(T)), 1);
1277 0 : GEN sk = Flx_Fl_add(Flx_neg(j, p), 1728%p, p);
1278 0 : GEN sA4 = Flx_triple(Flxq_mul(sk, j, T, p), p);
1279 0 : GEN u = Flxq_div(a4, sA4, T, p);
1280 0 : ulong ns = lgpol(s) ? Fl_neg(s[2], p): 0UL;
1281 0 : GEN P = mkvecsmall4(T[1], m[2], ns, 1L);
1282 : GEN A4, A6, t, tP;
1283 0 : Flxq_ellj_to_a4a6(polx_Flx(T[1]), P, p, &A4, &A6);
1284 0 : tP = addis(q, 1 - Flxq_ellcard_naive(A4, A6, P, p));
1285 0 : t = elltrace_extension(tP, n>>1, q);
1286 0 : return Flxq_is2npower(u, 2, T, p) ? subii(q1,t): addii(q1,t);
1287 : }
1288 : }
1289 10248 : if (p<=7 || p==13) return Flxq_ellcard_Kohel(a4, a6, T, p);
1290 21 : else return Flxq_ellcard_Harley(a4, a6, T, p);
1291 : }
1292 : }
1293 :
1294 : static GEN
1295 0 : Flxq_ellcard_Kedlaya(GEN a4, GEN a6, GEN T, ulong p)
1296 : {
1297 0 : pari_sp av = avma;
1298 0 : GEN H = mkpoln(4, gen_1, gen_0, Flx_to_ZX(a4), Flx_to_ZX(a6));
1299 0 : GEN Tp = Flx_to_ZX(get_Flx_mod(T));
1300 0 : long n = degpol(Tp), e = ((p < 16 ? n+1: n)>>1)+1;
1301 0 : GEN M = ZlXQX_hyperellpadicfrobenius(H, Tp, p, e);
1302 0 : GEN N = ZpXQM_prodFrobenius(M, Tp, utoipos(p), e);
1303 0 : GEN q = powuu(p, e);
1304 0 : GEN tp = Fq_add(gcoeff(N,1,1), gcoeff(N,2,2), Tp, q);
1305 0 : GEN t = Fp_center_i(typ(tp)==t_INT ? tp: leading_coeff(tp), q, shifti(q,-1));
1306 0 : return gerepileupto(av, subii(addiu(powuu(p, n), 1), t));
1307 : }
1308 :
1309 : GEN
1310 57299 : Flxq_ellj(GEN a4, GEN a6, GEN T, ulong p)
1311 : {
1312 57299 : pari_sp av=avma;
1313 57299 : if (p==3)
1314 : {
1315 : GEN J;
1316 14805 : if (typ(a4)!=t_VEC) return pol0_Flx(get_Flx_var(T));
1317 14805 : J = Flxq_div(Flxq_powu(gel(a4,1),3, T, p),Flx_neg(a6,p), T, p);
1318 14805 : return gerepileuptoleaf(av, J);
1319 : }
1320 : else
1321 : {
1322 42494 : pari_sp av=avma;
1323 42494 : GEN a43 = Flxq_mul(a4,Flxq_sqr(a4,T,p),T,p);
1324 42494 : GEN a62 = Flxq_sqr(a6,T,p);
1325 42494 : GEN num = Flx_mulu(a43,6912,p);
1326 42494 : GEN den = Flx_add(Flx_mulu(a43,4,p),Flx_mulu(a62,27,p),p);
1327 42494 : return gerepileuptoleaf(av, Flxq_div(num, den, T, p));
1328 : }
1329 : }
1330 :
1331 : void
1332 0 : Flxq_ellj_to_a4a6(GEN j, GEN T, ulong p, GEN *pt_a4, GEN *pt_a6)
1333 : {
1334 0 : ulong zagier = 1728 % p;
1335 0 : if (lgpol(j)==0)
1336 0 : { *pt_a4 = pol0_Flx(T[1]); *pt_a6 =pol1_Flx(T[1]); }
1337 0 : else if (lgpol(j)==1 && uel(j,2) == zagier)
1338 0 : { *pt_a4 = pol1_Flx(T[1]); *pt_a6 =pol0_Flx(T[1]); }
1339 : else
1340 : {
1341 0 : GEN k = Flx_Fl_add(Flx_neg(j, p), zagier, p);
1342 0 : GEN kj = Flxq_mul(k, j, T, p);
1343 0 : GEN k2j = Flxq_mul(kj, k, T, p);
1344 0 : *pt_a4 = Flx_triple(kj, p);
1345 0 : *pt_a6 = Flx_double(k2j, p);
1346 : }
1347 0 : }
1348 :
1349 : static GEN
1350 9205 : F3xq_ellcardj(GEN a4, GEN a6, GEN T, GEN q, long n)
1351 : {
1352 9205 : const ulong p = 3;
1353 : ulong t;
1354 9205 : GEN q1 = addiu(q,1);
1355 9205 : GEN na4 = Flx_neg(a4,p), ra4;
1356 9205 : if (!Flxq_issquare(na4,T,p))
1357 4823 : return q1;
1358 4382 : ra4 = Flxq_sqrt(na4,T,p);
1359 4382 : t = Flxq_trace(Flxq_div(a6,Flxq_mul(na4,ra4,T,p),T,p),T,p);
1360 4382 : if (n%2==1)
1361 : {
1362 : GEN q3;
1363 2261 : if (t==0) return q1;
1364 812 : q3 = powuu(p,(n+1)>>1);
1365 812 : return (t==1)^(n%4==1) ? subii(q1,q3): addii(q1,q3);
1366 : }
1367 : else
1368 : {
1369 2121 : GEN q22, q2 = powuu(p,n>>1);
1370 2121 : GEN W = Flxq_pow(a4,shifti(q,-2),T,p);
1371 2121 : long s = (W[2]==1)^(n%4==2);
1372 2121 : if (t!=0) return s ? addii(q1,q2): subii(q1, q2);
1373 2121 : q22 = shifti(q2,1);
1374 2121 : return s ? subii(q1,q22): addii(q1, q22);
1375 : }
1376 : }
1377 :
1378 : static GEN
1379 15820 : Flxq_ellcardj(GEN a4, GEN a6, ulong j, GEN T, GEN q, ulong p, long n)
1380 : {
1381 15820 : GEN q1 = addiu(q,1);
1382 15820 : if (j==0)
1383 : {
1384 : ulong w;
1385 : GEN W, t, N;
1386 5677 : if (umodiu(q,6)!=1) return q1;
1387 4277 : N = Fp_ffellcard(gen_0,gen_1,q,n,utoipos(p));
1388 4277 : t = subii(q1, N);
1389 4277 : W = Flxq_pow(a6,diviuexact(shifti(q,-1), 3),T,p);
1390 4277 : if (degpol(W)>0) /*p=5 mod 6*/
1391 1323 : return Flx_equal1(Flxq_powu(W,3,T,p)) ? addii(q1,shifti(t,-1)):
1392 441 : subii(q1,shifti(t,-1));
1393 3395 : w = W[2];
1394 3395 : if (w==1) return N;
1395 2653 : if (w==p-1) return addii(q1,t);
1396 : else /*p=1 mod 6*/
1397 : {
1398 1862 : GEN u = shifti(t,-1), v = sqrtint(diviuexact(subii(q,sqri(u)),3));
1399 1862 : GEN a = addii(u,v), b = shifti(v,1);
1400 1862 : if (Fl_powu(w,3,p)==1)
1401 : {
1402 931 : if (Fl_add(umodiu(a,p),Fl_mul(w,umodiu(b,p),p),p)==0)
1403 469 : return subii(q1,subii(shifti(b,1),a));
1404 : else
1405 462 : return addii(q1,addii(a,b));
1406 : }
1407 : else
1408 : {
1409 931 : if (Fl_sub(umodiu(a,p),Fl_mul(w,umodiu(b,p),p),p)==0)
1410 469 : return subii(q1,subii(a,shifti(b,1)));
1411 : else
1412 462 : return subii(q1,addii(a,b));
1413 : }
1414 : }
1415 10143 : } else if (j==1728%p)
1416 : {
1417 : ulong w;
1418 : GEN W, N, t;
1419 5656 : if (mod4(q)==3) return q1;
1420 4256 : W = Flxq_pow(a4,shifti(q,-2),T,p);
1421 4256 : if (degpol(W)>0) return q1; /*p=3 mod 4*/
1422 3612 : w = W[2];
1423 3612 : N = Fp_ffellcard(gen_1,gen_0,q,n,utoipos(p));
1424 3612 : if(w==1) return N;
1425 2534 : t = subii(q1, N);
1426 2534 : if(w==p-1) return addii(q1, t);
1427 : else /*p=1 mod 4*/
1428 : {
1429 1400 : GEN u = shifti(t,-1), v = sqrtint(subii(q,sqri(u)));
1430 1400 : if (Fl_add(umodiu(u,p),Fl_mul(w,umodiu(v,p),p),p)==0)
1431 700 : return subii(q1,shifti(v,1));
1432 : else
1433 700 : return addii(q1,shifti(v,1));
1434 : }
1435 : } else
1436 : {
1437 4487 : ulong g = Fl_div(j, Fl_sub(1728%p, j, p), p);
1438 4487 : GEN N = Fp_ffellcard(utoi(Fl_triple(g,p)),utoi(Fl_double(g,p)),q,n,utoipos(p));
1439 4487 : GEN l = Flxq_mul(Flx_triple(a6,p),Flx_double(a4,p),T,p);
1440 4487 : if (Flxq_issquare(l,T,p)) return N;
1441 2940 : return subii(shifti(q1,1),N);
1442 : }
1443 : }
1444 :
1445 : static GEN
1446 454 : Flxq_ffellcard(GEN a4, GEN a6, GEN M, GEN q, GEN T, ulong p, long n)
1447 : {
1448 454 : long m = degpol(M);
1449 454 : GEN j = polx_Flx(M[1]);
1450 454 : GEN g = Flxq_div(j, mkvecsmall3(M[1],1728%p,p-1), M, p);
1451 454 : GEN N = Flxq_ellcard(Flx_triple(g, p), Flx_double(g, p), M, p);
1452 454 : GEN qm = powuu(p, m), q1 = addiu(q, 1), qm1 = addiu(qm, 1);
1453 454 : GEN l = Flxq_mul(Flx_triple(a6,p), Flx_double(a4,p), T, p);
1454 454 : GEN te = elltrace_extension(subii(qm1, N), n/m, qm);
1455 454 : return Flxq_issquare(l,T,p) ? subii(q1, te): addii(q1, te);
1456 : }
1457 :
1458 : static GEN
1459 62661 : Flxq_ellcard_i(GEN a4, GEN a6, GEN T, ulong p)
1460 : {
1461 62661 : long n = get_Flx_degree(T);
1462 62661 : GEN J, M, q = powuu(p, n);
1463 62661 : if (typ(a4)==t_VEC)
1464 20993 : return F3xq_ellcard(gel(a4,1), a6, T);
1465 41668 : if (p==3)
1466 9205 : return F3xq_ellcardj(a4, a6, T, q, n);
1467 32463 : if (degpol(a4)<=0 && degpol(a6)<=0)
1468 217 : return Fp_ffellcard(utoi(Flx_eval(a4,0,p)),utoi(Flx_eval(a6,0,p)),q,n,utoipos(p));
1469 32246 : J = Flxq_ellj(a4,a6,T,p);
1470 32246 : if (degpol(J)<=0)
1471 15820 : return Flxq_ellcardj(a4,a6,lgpol(J)?J[2]:0,T,q,p,n);
1472 16426 : M = Flxq_minpoly(J, T, p);
1473 16426 : if (degpol(M) < n)
1474 454 : return Flxq_ffellcard(a4, a6, M, q, T, p, n);
1475 15972 : if (p <= 7)
1476 10990 : return Flxq_ellcard_Satoh(a4, a6, J, T, p);
1477 4982 : if (cmpis(q,100)<0)
1478 0 : return utoi(Flxq_ellcard_naive(a4, a6, T, p));
1479 4982 : if (p == 13 || (7*p <= (ulong)10*n && (BITS_IN_LONG==64 || p <= 103)))
1480 154 : return Flxq_ellcard_Satoh(a4, a6, J, T, p);
1481 4828 : if (p <= (ulong)2*n)
1482 0 : return Flxq_ellcard_Kedlaya(a4, a6, T, p);
1483 4828 : if (expi(q)<=62)
1484 4774 : return Flxq_ellcard_Shanks(a4, a6, q, T, p);
1485 : else
1486 54 : return Fq_ellcard_SEA(Flx_to_ZX(a4),Flx_to_ZX(a6),q,Flx_to_ZX(T),utoipos(p),0);
1487 : }
1488 :
1489 : GEN
1490 62661 : Flxq_ellcard(GEN a4, GEN a6, GEN T, ulong p)
1491 : {
1492 62661 : pari_sp av = avma;
1493 62661 : return gerepileuptoint(av, Flxq_ellcard_i(a4, a6, T, p));
1494 : }
1495 :
1496 : static long
1497 350 : Fl_ellj_trace(ulong j, ulong p)
1498 : {
1499 : ulong a4, a6;
1500 350 : Fl_ellj_to_a4a6(j, p, &a4, &a6);
1501 350 : return Fl_elltrace(a4, a6, p);
1502 : }
1503 :
1504 : /* Given ordinary E/Fq, a prime ell, and the height of the ell-volcano
1505 : * containing j(E) (= v_ell(conductor of Z[pi_E]) returns the height of j(E)
1506 : * on its ell-volcano (= v_ell(conductor of the order End(E)). */
1507 : static long
1508 147 : Fl_ellheightabovefloor(ulong j, long ell, long e, ulong p)
1509 : {
1510 147 : pari_sp av = avma;
1511 : GEN Xp, G, phi, phix, j0, j1;
1512 : long h, i, nj1;
1513 147 : if (e==0) return 0;
1514 147 : if (j==0 || j==1728%p) return e;
1515 119 : phi = ZXX_to_FlxX(polmodular_ZXX(ell, 0, 0, 1), p, 1);
1516 119 : phix = FlxY_evalx(phi, j, p);
1517 119 : Xp = Flx_Frobenius(phix, p);
1518 119 : G = Flx_gcd(Flx_sub(Xp, polx_Flx(0), p), phix, p);
1519 119 : nj1 = degpol(G);
1520 119 : if (nj1 < ell) return 0;
1521 112 : if (e==1 || nj1 != ell+1) return e;
1522 21 : j1 = Flx_roots(G, p);
1523 21 : nj1 = lg(j1)-1;
1524 21 : if (nj1 < 3) return 0;
1525 21 : j0 = mkvecsmall3(j,j,j);
1526 42 : for (h = 1; ; h++)
1527 133 : for(i = 1; i <= 3; i++)
1528 : {
1529 112 : GEN P = Flx_div_by_X_x(FlxY_evalx(phi, uel(j1,i), p), uel(j0,i), p, NULL);
1530 112 : GEN r = Flx_roots(P, p);
1531 112 : if (lg(r) == 1) return gc_long(av, h);
1532 91 : j0[i] = j1[i];
1533 91 : j1[i] = r[1];
1534 : }
1535 : }
1536 :
1537 : /* Given an ordinary elliptic curve E/Fp and an integer h, returns
1538 : * D = disc(End(E)) assuming h(D) = h, using the approach sketched in
1539 : * Remark 13. If the algorithm returns 0 it has proved that h(D) != h, but it
1540 : * is under no obligation to do so and is allowed to return any value when the
1541 : * assumption h(d) = h is false. */
1542 : static long
1543 350 : Fl_end13(ulong j, ulong h, ulong p)
1544 : {
1545 : ulong D0, v, h0;
1546 : long i, lL, lc, lD, nc;
1547 : GEN D, DF, cs, L, vP, vE;
1548 350 : ulong t = Fl_ellj_trace(j, p);
1549 :
1550 350 : D0 = coredisc2u_fact(factoru(4*p-t*t), -1, &vP, &vE);
1551 350 : h0 = itou(classno(stoi(-D0)));
1552 350 : if (h % h0 != 0) return 0;
1553 350 : h /= h0;
1554 350 : D = divisorsu_fact_factored(mkmat2(vP,vE));
1555 350 : DF = gel(D,2); D = gel(D,1);
1556 350 : lD = lg(D); v = D[lD-1];
1557 350 : cs = cgetg(lD,t_VECSMALL); nc = 0;
1558 1848 : for (i = 1; i < lD; i++)
1559 : {
1560 1498 : GEN F = gel(DF,i);
1561 1498 : ulong w = uquadclassnoF_fact(D0, -1, gel(F,1), gel(F,2));
1562 1498 : if (w == h) uel(cs,++nc) = v / uel(D,i);
1563 : }
1564 350 : if (nc==0) return 0;
1565 350 : if (nc==1) { v /= uel(cs,1); return -D0*v*v; }
1566 147 : L = cgetg(nc+1, t_VEC);
1567 448 : for (i = 1; i <= nc; i++) gel(L,i) = gel(factoru(uel(cs,i)), 1);
1568 147 : L = vecsmall_uniq(shallowconcat1(L));
1569 147 : lL = lg(L); lc = nc+1;
1570 147 : for (i = 1; i < lL; i++)
1571 : {
1572 147 : ulong ell = L[i];
1573 147 : long k, e = Fl_ellheightabovefloor(j, ell, z_lval(v,ell), p);
1574 448 : for (k = 1; k < lc; k++)
1575 301 : if(cs[k] && z_lval(cs[k], ell) != e) { cs[k] = 0; nc--; }
1576 147 : if (nc==0) return 0;
1577 147 : if (nc==1)
1578 : {
1579 161 : for (k = 1; k < lc; k++)
1580 161 : if (cs[k]) { v /= uel(cs,k); return -D0*v*v; }
1581 : }
1582 : }
1583 0 : return 0;
1584 : }
1585 :
1586 : INLINE int
1587 350 : RgX_is_monic_ZX(GEN pol)
1588 350 : { return RgX_is_ZX(pol) && ZX_is_monic(pol); }
1589 :
1590 : long
1591 357 : polisclass(GEN H)
1592 : {
1593 357 : pari_sp av = avma, btop;
1594 357 : long h = degpol(H), hl, i, pmin, vH = varn(H), vh;
1595 : double lmin;
1596 : ulong p;
1597 : GEN h2list;
1598 : forprime_t T;
1599 :
1600 357 : if (typ(H)!= t_POL) pari_err_TYPE("polsisclass",H);
1601 357 : if (h <= 0 || !RgX_is_monic_ZX(H)) return 0;
1602 350 : vh = vals(h);
1603 350 : h2list = cgetg(vh+2, t_VECSMALL); hl = 1;
1604 1071 : for (i = 0; i <= vh; i++)
1605 : {
1606 721 : ulong d = 1UL<<i;
1607 721 : if (((d-h)&1)==0) h2list[hl++] = d;
1608 : }
1609 350 : setlg(h2list, hl);
1610 350 : lmin = h * (log(log(h+2))+2);
1611 350 : pmin = 33 * ceil(lmin*lmin);
1612 350 : u_forprime_init(&T, pmin, ULONG_MAX);
1613 350 : btop = avma;
1614 3339 : while((p = u_forprime_next(&T)))
1615 : {
1616 : ulong r;
1617 : long D, nroots;
1618 3339 : GEN Xp, G, Hp = ZX_to_Flx(H,p);
1619 3339 : if (!Flx_is_squarefree(Hp, p)) { set_avma(btop); continue; }
1620 3339 : Xp = Flx_Frobenius(Hp, p);
1621 3339 : G = Flx_gcd(Flx_sub(Xp, polx_Flx(evalvarn(vH)), p), Hp, p);
1622 3339 : nroots = degpol(G);
1623 3339 : if (nroots==0) { set_avma(btop); continue; }
1624 1274 : if (nroots < h && !zv_search(h2list,nroots)) return gc_long(av, 0);
1625 1274 : r = Flx_oneroot(G, p);
1626 1274 : if (Fp_elljissupersingular(utoi(r), utoi(p))) { set_avma(btop); continue; }
1627 350 : D = Fl_end13(r, h, p);
1628 350 : if (D && gequal(H, polclass(stoi(D), 0, vH))) return gc_long(av, D);
1629 0 : return gc_long(av, 0);
1630 : }
1631 0 : pari_err_BUG("polisclass");
1632 : return 0; /* LCOV_EXCL_LINE */
1633 : }
|
