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 1743 : RgE_to_FlxqE(GEN x, GEN T, ulong p)
34 : {
35 1743 : if (ell_is_inf(x)) return x;
36 1743 : retmkvec2(Rg_to_Flxq(gel(x,1),T,p), Rg_to_Flxq(gel(x,2),T,p));
37 : }
38 :
39 : GEN
40 93775 : FlxqE_changepoint(GEN x, GEN ch, GEN T, ulong p)
41 : {
42 93775 : pari_sp av = avma;
43 : GEN p1, p2, z, u, r, s, t, v, v2, v3;
44 : ulong pi;
45 93775 : if (ell_is_inf(x)) return x;
46 93768 : pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
47 93768 : u = gel(ch,1); r = gel(ch,2);
48 93768 : s = gel(ch,3); t = gel(ch,4);
49 93768 : v = Flxq_inv_pre(u, T, p, pi);
50 93768 : v2 = Flxq_sqr_pre(v, T, p, pi);
51 93768 : v3 = Flxq_mul_pre(v,v2, T, p, pi);
52 93768 : p1 = Flx_sub(gel(x,1), r, p);
53 93768 : p2 = Flx_sub(gel(x,2), Flx_add(Flxq_mul_pre(s, p1, T, p, pi),t, p), p);
54 93768 : z = cgetg(3,t_VEC);
55 93768 : gel(z,1) = Flxq_mul_pre(v2, p1, T, p, pi);
56 93768 : gel(z,2) = Flxq_mul_pre(v3, p2, T, p, pi);
57 93768 : return gc_upto(av, z);
58 : }
59 :
60 : GEN
61 1743 : FlxqE_changepointinv(GEN x, GEN ch, GEN T, ulong p)
62 : {
63 1743 : pari_sp av = avma;
64 : GEN p1, p2, u, r, s, t, X, Y, u2, u3, u2X, z;
65 : ulong pi;
66 1743 : if (ell_is_inf(x)) return x;
67 1743 : pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
68 1743 : X = gel(x,1); Y = gel(x,2);
69 1743 : u = gel(ch,1); r = gel(ch,2);
70 1743 : s = gel(ch,3); t = gel(ch,4);
71 1743 : u2 = Flxq_sqr_pre(u, T, p, pi);
72 1743 : u3 = Flxq_mul_pre(u,u2, T, p, pi);
73 1743 : u2X = Flxq_mul_pre(u2,X, T, p, pi);
74 1743 : p1 = Flxq_mul_pre(u3,Y, T, p, pi);
75 1743 : p2 = Flx_add(Flxq_mul_pre(s,u2X, T, p, pi), t, p);
76 1743 : z = cgetg(3, t_VEC);
77 1743 : gel(z,1) = Flx_add(u2X, r, p);
78 1743 : gel(z,2) = Flx_add(p1, p2, p);
79 1743 : return gc_upto(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 30118 : set_avma(av);
93 30118 : a = random_Flx(n, vs, p);
94 30118 : } 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 752214 : FlxqE_dbl_slope(GEN P, GEN a4, GEN T, ulong p, ulong pi, GEN *ps)
117 : {
118 : GEN x, y, Q, s;
119 752214 : if (ell_is_inf(P) || !lgpol(gel(P,2))) return ellinf();
120 717556 : x = gel(P,1); y = gel(P,2);
121 717556 : if (p==3UL)
122 263190 : s = typ(a4)==t_VEC? Flxq_div_pre(Flxq_mul_pre(x, gel(a4,1), T,p,pi), y, T,p,pi)
123 263190 : : Flxq_div_pre(a4, Flx_neg(y, p), T,p,pi);
124 : else
125 : {
126 454366 : GEN sx = Flx_add(Flx_triple(Flxq_sqr_pre(x, T, p, pi), p), a4, p);
127 454366 : s = Flxq_div_pre(sx, Flx_double(y, p), T, p, pi);
128 : }
129 717556 : Q = cgetg(3,t_VEC);
130 717556 : gel(Q,1) = Flx_sub(Flxq_sqr_pre(s, T, p, pi), Flx_double(x, p), p);
131 717556 : if (typ(a4)==t_VEC) gel(Q, 1) = Flx_sub(gel(Q,1), gel(a4,1), p);
132 717556 : gel(Q,2) = Flx_sub(Flxq_mul_pre(s, Flx_sub(x, gel(Q,1), p), T, p, pi),
133 : y, p);
134 717556 : if (ps) *ps = s;
135 717556 : 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 gc_upto(av, FlxqE_dbl_slope(P,a4, T, p, pi, NULL));
144 : }
145 :
146 : static GEN
147 414410 : 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 414410 : if (ell_is_inf(P)) return Q;
151 411510 : if (ell_is_inf(Q)) return P;
152 411496 : Px = gel(P,1); Py = gel(P,2);
153 411496 : Qx = gel(Q,1); Qy = gel(Q,2);
154 411496 : if (Flx_equal(Px, Qx))
155 : {
156 22157 : if (Flx_equal(Py, Qy))
157 4121 : return FlxqE_dbl_slope(P, a4, T, p, pi, ps);
158 : else
159 18036 : return ellinf();
160 : }
161 389339 : s = Flxq_div_pre(Flx_sub(Py, Qy, p), Flx_sub(Px, Qx, p), T, p, pi);
162 389339 : R = cgetg(3,t_VEC);
163 389339 : gel(R,1) = Flx_sub(Flx_sub(Flxq_sqr_pre(s, T, p, pi), Px, p), Qx, p);
164 389339 : if (typ(a4)==t_VEC) gel(R,1) = Flx_sub(gel(R,1), gel(a4,1), p);
165 389339 : gel(R,2) = Flx_sub(Flxq_mul_pre(s, Flx_sub(Px, gel(R,1), p), T, p, pi),
166 : Py, p);
167 389339 : if (ps) *ps = s;
168 389339 : 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 gc_upto(av, FlxqE_add_slope(P,Q,a4, T,p,pi, NULL));
177 : }
178 :
179 : static GEN
180 6742 : FlxqE_neg_i(GEN P, ulong p)
181 : {
182 6742 : if (ell_is_inf(P)) return P;
183 6742 : return mkvec2(gel(P,1), Flx_neg(gel(P,2), p));
184 : }
185 :
186 : GEN
187 1060 : FlxqE_neg(GEN P, GEN T, ulong p)
188 : {
189 : (void) T;
190 1060 : if (ell_is_inf(P)) return ellinf();
191 1060 : 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 gc_upto(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 722088 : _FlxqE_dbl(void *E, GEN P)
210 : {
211 722088 : struct _FlxqE *e = (struct _FlxqE *) E;
212 722088 : return FlxqE_dbl_slope(P, e->a4, e->T, e->p, e->pi, NULL);
213 : }
214 :
215 : static GEN
216 404701 : _FlxqE_add(void *E, GEN P, GEN Q)
217 : {
218 404701 : struct _FlxqE *e = (struct _FlxqE *) E;
219 404701 : return FlxqE_add_slope(P, Q, e->a4, e->T, e->p, e->pi, NULL);
220 : }
221 :
222 : static GEN
223 6742 : _FlxqE_sub(void *E, GEN P, GEN Q)
224 : {
225 6742 : struct _FlxqE *e = (struct _FlxqE *) E;
226 6742 : 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 215795 : _FlxqE_mul(void *E, GEN P, GEN n)
231 : {
232 215795 : pari_sp av = avma;
233 215795 : struct _FlxqE *e=(struct _FlxqE *) E;
234 215795 : long s = signe(n);
235 215795 : if (!s || ell_is_inf(P)) return ellinf();
236 215115 : if (s < 0) P = FlxqE_neg(P, e->T, e->p);
237 215115 : if (is_pm1(n)) return s>0? gcopy(P): P;
238 199980 : return gc_GEN(av, gen_pow_i(P, n, e, &_FlxqE_dbl, &_FlxqE_add));
239 : }
240 :
241 : GEN
242 7 : FlxqE_mul(GEN P, GEN n, GEN a4, GEN T, ulong p)
243 : {
244 : struct _FlxqE E;
245 7 : E.a4= a4; E.T = T; E.p = p; E.pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
246 7 : 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 77165 : random_F3xqE(GEN a2, GEN a6, GEN T)
254 : {
255 77165 : pari_sp ltop = avma;
256 : GEN x, y, rhs;
257 77165 : const ulong p = 3;
258 : do
259 : {
260 152676 : set_avma(ltop);
261 152676 : x = random_Flx(get_Flx_degree(T),get_Flx_var(T),p);
262 152676 : rhs = Flx_add(Flxq_mul(Flxq_sqr(x, T, p), Flx_add(x, a2, p), T, p), a6, p);
263 152676 : } while ((!lgpol(rhs) && !lgpol(x)) || !Flxq_issquare(rhs, T, p));
264 77165 : y = Flxq_sqrt(rhs, T, p);
265 77165 : if (!y) pari_err_PRIME("random_F3xqE", T);
266 77165 : return gc_GEN(ltop, mkvec2(x, y));
267 : }
268 :
269 : /* Finds a random nonsingular point on E */
270 : static GEN
271 151548 : random_FlxqE_pre(GEN a4, GEN a6, GEN T, ulong p, ulong pi)
272 : {
273 151548 : pari_sp ltop = avma;
274 : GEN x, x2, y, rhs;
275 151548 : if (typ(a4)==t_VEC) return random_F3xqE(gel(a4,1), a6, T);
276 : do
277 : {
278 144091 : set_avma(ltop);
279 144091 : x = random_Flx(get_Flx_degree(T),get_Flx_var(T),p);
280 144091 : x2 = Flxq_sqr_pre(x, T, p, pi); /* x^3+a4*x+a6 = x*(x^2+a4)+a6 */
281 144091 : rhs = Flx_add(Flxq_mul_pre(x, Flx_add(x2, a4, p), T, p, pi), a6, p);
282 149414 : } while ((!lgpol(rhs) && !lgpol(Flx_add(Flx_triple(x2, p), a4, p)))
283 149351 : || !Flxq_issquare(rhs, T, p));
284 74383 : y = Flxq_sqrt(rhs, T, p);
285 74383 : if (!y) pari_err_PRIME("random_FlxqE", T);
286 74383 : return gc_GEN(ltop, mkvec2(x, y));
287 : }
288 : GEN
289 76955 : random_FlxqE(GEN a4, GEN a6, GEN T, ulong p)
290 76955 : { return random_FlxqE_pre(a4, a6, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
291 :
292 : static GEN
293 74593 : _FlxqE_rand(void *E)
294 : {
295 74593 : struct _FlxqE *e=(struct _FlxqE *) E;
296 74593 : 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 54 : get_FlxqE_group(void ** pt_E, GEN a4, GEN a6, GEN T, ulong p)
303 : {
304 54 : struct _FlxqE *e = (struct _FlxqE *) stack_malloc(sizeof(struct _FlxqE));
305 54 : e->a4 = a4; e->a6 = a6;
306 54 : e->pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
307 54 : e->p = p;
308 54 : e->T = Flx_get_red_pre(T, p, e->pi);
309 54 : *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 gc_INT(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 gc_INT(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 68328 : FlxqE_vert(GEN P, GEN Q, GEN a4, GEN T, ulong p, ulong pi)
336 : {
337 68328 : long vT = get_Flx_var(T);
338 : GEN df;
339 68328 : if (ell_is_inf(P)) return pol1_Flx(vT);
340 46748 : if (!Flx_equal(gel(Q,1), gel(P,1))) return Flx_sub(gel(Q,1), gel(P,1), p);
341 1776 : if (lgpol(gel(P,2))!=0) return pol1_Flx(vT);
342 734 : df = typ(a4)==t_VEC ? Flxq_mul_pre(gel(P,1), Flx_double(gel(a4,1), p), T,p,pi)
343 1166 : : a4;
344 1166 : 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 28972 : FlxqE_Miller_line(GEN R, GEN Q, GEN slope, GEN a4, GEN T, ulong p, ulong pi)
350 : {
351 28972 : long vT = get_Flx_var(T);
352 28972 : GEN x = gel(Q,1), y = gel(Q,2);
353 28972 : GEN tmp1 = Flx_sub(x, gel(R,1), p);
354 28972 : GEN tmp2 = Flx_add(Flxq_mul_pre(tmp1, slope, T, p, pi), gel(R,2), p);
355 28972 : if (!Flx_equal(y, tmp2)) return Flx_sub(y, tmp2, p);
356 1108 : if (lgpol(y) == 0) return pol1_Flx(vT);
357 : else
358 : {
359 525 : GEN s1, s2, a2 = typ(a4)==t_VEC ? gel(a4,1): NULL;
360 525 : GEN y2i = Flxq_inv_pre(Flx_mulu(y, 2, p), T, p, pi);
361 525 : GEN df = a2 ? Flxq_mul_pre(x, Flx_mulu(a2, 2, p), T, p, pi): a4;
362 : GEN x3, ddf;
363 525 : s1 = Flxq_mul_pre(Flx_add(Flx_triple(Flxq_sqr_pre(x, T, p, pi), p), df, p), y2i, T, p, pi);
364 525 : if (!Flx_equal(s1, slope)) return Flx_sub(s1, slope, p);
365 206 : x3 = Flx_triple(x, p);
366 206 : ddf = a2 ? Flx_add(x3, a2, p): x3;
367 206 : s2 = Flxq_mul_pre(Flx_sub(ddf, Flxq_sqr_pre(s1, T,p,pi), p), y2i, T,p,pi);
368 206 : 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 46211 : FlxqE_tangent_update(GEN R, GEN Q, GEN a4, GEN T, ulong p, ulong pi, GEN *pt_R)
376 : {
377 46211 : if (ell_is_inf(R))
378 : {
379 3938 : *pt_R = ellinf();
380 3938 : return pol1_Flx(get_Flx_var(T));
381 : }
382 42273 : else if (!lgpol(gel(R,2)))
383 : {
384 16268 : *pt_R = ellinf();
385 16268 : return FlxqE_vert(R, Q, a4, T, p, pi);
386 : } else {
387 : GEN slope;
388 26005 : *pt_R = FlxqE_dbl_slope(R, a4, T, p, pi, &slope);
389 26005 : 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 4409 : FlxqE_chord_update(GEN R, GEN P, GEN Q, GEN a4, GEN T, ulong p, ulong pi, GEN *pt_R)
397 : {
398 4409 : if (ell_is_inf(R))
399 : {
400 67 : *pt_R = gcopy(P);
401 67 : return FlxqE_vert(P, Q, a4, T, p, pi);
402 : }
403 4342 : 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 4342 : else if (Flx_equal(gel(P, 1), gel(R, 1)))
409 : {
410 1375 : if (Flx_equal(gel(P, 2), gel(R, 2)))
411 1 : return FlxqE_tangent_update(R, Q, a4, T, p, pi, pt_R);
412 : else
413 : {
414 1374 : *pt_R = ellinf();
415 1374 : return FlxqE_vert(R, Q, a4, T, p, pi);
416 : }
417 : } else {
418 : GEN slope;
419 2967 : *pt_R = FlxqE_add_slope(P, R, a4, T, p, pi, &slope);
420 2967 : 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 46210 : FlxqE_Miller_dbl(void* E, GEN d)
432 : {
433 46210 : struct _FlxqE_miller *m = (struct _FlxqE_miller *)E;
434 46210 : ulong p = m->p, pi = m->pi;
435 46210 : GEN T = m->T, a4 = m->a4, P = m->P;
436 46210 : GEN v, line, point = gel(d,3);
437 46210 : GEN N = Flxq_sqr_pre(gel(d,1), T, p, pi);
438 46210 : GEN D = Flxq_sqr_pre(gel(d,2), T, p, pi);
439 46210 : line = FlxqE_tangent_update(point, P, a4, T, p, pi, &point);
440 46210 : N = Flxq_mul_pre(N, line, T, p, pi);
441 46210 : v = FlxqE_vert(point, P, a4, T, p, pi);
442 46210 : D = Flxq_mul_pre(D, v, T, p, pi); return mkvec3(N, D, point);
443 : }
444 :
445 : static GEN
446 4409 : FlxqE_Miller_add(void* E, GEN va, GEN vb)
447 : {
448 4409 : struct _FlxqE_miller *m = (struct _FlxqE_miller *)E;
449 4409 : ulong p = m->p, pi = m->pi;
450 4409 : GEN T = m->T, a4 = m->a4, P = m->P;
451 : GEN v, line, point;
452 4409 : GEN na = gel(va,1), da = gel(va,2), pa = gel(va,3);
453 4409 : GEN nb = gel(vb,1), db = gel(vb,2), pb = gel(vb,3);
454 4409 : GEN N = Flxq_mul_pre(na, nb, T, p, pi);
455 4409 : GEN D = Flxq_mul_pre(da, db, T, p, pi);
456 4409 : line = FlxqE_chord_update(pa, pb, P, a4, T, p, pi, &point);
457 4409 : N = Flxq_mul_pre(N, line, T, p, pi);
458 4409 : v = FlxqE_vert(point, P, a4, T, p, pi);
459 4409 : 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 17575 : FlxqE_Miller(GEN Q, GEN P, GEN m, GEN a4, GEN T, ulong p, ulong pi)
466 : {
467 17575 : pari_sp av = avma;
468 : struct _FlxqE_miller d;
469 : GEN v, N, D, g1;
470 :
471 17575 : d.a4 = a4; d.T = T; d.p = p; d.P = P; d.pi = pi;
472 17575 : g1 = pol1_Flx(get_Flx_var(T));
473 17575 : v = gen_pow_i(mkvec3(g1,g1,Q), m, (void*)&d,
474 : FlxqE_Miller_dbl, FlxqE_Miller_add);
475 17575 : N = gel(v,1); D = gel(v,2);
476 17575 : return gc_upto(av, Flxq_div_pre(N, D, T, p, pi));
477 : }
478 :
479 : GEN
480 14292 : FlxqE_weilpairing_pre(GEN P, GEN Q, GEN m, GEN a4, GEN T, ulong p, ulong pi)
481 : {
482 14292 : pari_sp av = avma;
483 : GEN N, D, w;
484 14292 : if (ell_is_inf(P) || ell_is_inf(Q)
485 11475 : || (Flx_equal(gel(P,1),gel(Q,1)) && Flx_equal(gel(P,2),gel(Q,2))))
486 5536 : return pol1_Flx(get_Flx_var(T));
487 8756 : N = FlxqE_Miller(P, Q, m, a4, T, p, pi);
488 8756 : D = FlxqE_Miller(Q, P, m, a4, T, p, pi);
489 8756 : w = Flxq_div_pre(N, D, T, p, pi); if (mpodd(m)) w = Flx_neg(w, p);
490 8756 : return gc_upto(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 14271 : _FlxqE_pairorder(void *E, GEN P, GEN Q, GEN m, GEN F)
505 : {
506 14271 : struct _FlxqE *e = (struct _FlxqE *) E;
507 14271 : 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 gc_GEN(av, P);
541 : }
542 :
543 : static GEN
544 31549 : ell_to_a4a6_F3xq(GEN E, GEN T)
545 : {
546 : GEN a1, a3, b2, b4, b6;
547 31549 : a1 = Rg_to_Flxq(ell_get_a1(E),T,3);
548 31549 : a3 = Rg_to_Flxq(ell_get_a3(E),T,3);
549 31549 : b2 = Rg_to_Flxq(ell_get_b2(E),T,3);
550 31549 : b4 = Rg_to_Flxq(ell_get_b4(E),T,3);
551 31549 : b6 = Rg_to_Flxq(ell_get_b6(E),T,3);
552 31549 : if(lgpol(b2)) /* ordinary case */
553 : {
554 21896 : GEN b4b2 = Flxq_div(b4,b2,T,3);
555 21896 : GEN a6 = Flx_sub(b6,Flxq_mul(b4b2,Flx_add(b4,Flxq_sqr(b4b2,T,3),3),T,3),3);
556 21896 : retmkvec3(mkvec(b2), a6,
557 : mkvec4(Fl_to_Flx(1,T[1]),b4b2,Flx_neg(a1,3),Flx_neg(a3,3)));
558 : }
559 : else /* super-singular case */
560 9653 : retmkvec3(Flx_neg(b4, 3), b6,
561 : mkvec4(Fl_to_Flx(1,T[1]),zero_Flx(T[1]), Flx_neg(a1,3), Flx_neg(a3,3)));
562 : }
563 :
564 : GEN
565 73874 : ell_to_a4a6_Flxq(GEN E, GEN T, ulong p)
566 : {
567 : GEN a1, a3, b2, c4, c6;
568 73874 : if(p==3) return ell_to_a4a6_F3xq(E, T);
569 42325 : a1 = Rg_to_Flxq(ell_get_a1(E),T,p);
570 42325 : a3 = Rg_to_Flxq(ell_get_a3(E),T,p);
571 42325 : b2 = Rg_to_Flxq(ell_get_b2(E),T,p);
572 42325 : c4 = Rg_to_Flxq(ell_get_c4(E),T,p);
573 42325 : c6 = Rg_to_Flxq(ell_get_c6(E),T,p);
574 42325 : retmkvec3(Flx_neg(Flx_mulu(c4, 27, p), p), Flx_neg(Flx_mulu(c6, 54, p), p),
575 : mkvec4(Fl_to_Flx(6%p,T[1]), Flx_triple(b2,p), Flx_triple(a1,p),
576 : Flx_mulu(a3,108,p)));
577 : }
578 : /***********************************************************************/
579 : /** Point counting **/
580 : /***********************************************************************/
581 :
582 : /* assume a and n are coprime */
583 : static GEN
584 76251 : RgX_circular_shallow(GEN P, long a, long n)
585 : {
586 76251 : long i, l = lgpol(P);
587 76251 : GEN Q = cgetg(2+n,t_POL);
588 76251 : Q[1] = P[1];
589 512330 : for(i=0; i<l; i++)
590 436079 : gel(Q,2+(i*a)%n) = gel(P,2+i);
591 168693 : for( ; i<n; i++)
592 92442 : gel(Q,2+(i*a)%n) = gen_0;
593 76251 : return normalizepol_lg(Q,2+n);
594 : }
595 :
596 : static GEN
597 76251 : ZpXQ_frob_cyc(GEN x, GEN T, GEN q, ulong p)
598 : {
599 76251 : long n = get_FpX_degree(T);
600 76251 : return FpX_rem(RgX_circular_shallow(x,p,n+1), T, q);
601 : }
602 :
603 : static GEN
604 114016 : ZpXQ_frob(GEN x, GEN Xm, GEN T, GEN q, ulong p)
605 : {
606 114016 : if (lg(Xm)==1)
607 43428 : return ZpXQ_frob_cyc(x, T, q, p);
608 : else
609 : {
610 70588 : long n = get_FpX_degree(T);
611 70588 : GEN V = RgX_blocks(RgX_inflate(x, p), n, p);
612 70588 : GEN W = ZXV_dotproduct(V, Xm);
613 70588 : return FpX_rem(W, T, q);
614 : }
615 : }
616 :
617 : struct _lift_lin
618 : {
619 : ulong p, pi;
620 : GEN sqx, Tp, ai, Xm;
621 : };
622 :
623 : static GEN
624 84280 : _lift_invl(void *E, GEN x)
625 : {
626 84280 : struct _lift_lin *d = (struct _lift_lin *) E;
627 84280 : GEN T = d->Tp;
628 84280 : ulong p = d->p, pi = d->pi;
629 84280 : GEN xai = Flxq_mul_pre(ZX_to_Flx(x, p), d->ai, T, p, pi);
630 84280 : return Flx_to_ZX(Flxq_lroot_fast_pre(xai, d->sqx, T, p, pi));
631 : }
632 : static GEN
633 23884 : _lift_lin(void *E, GEN F, GEN x2, GEN q)
634 : {
635 23884 : struct _lift_lin *d = (struct _lift_lin *) E;
636 23884 : pari_sp av = avma;
637 23884 : GEN T = gel(F,3), Xm = gel(F,4);
638 23884 : GEN y2 = ZpXQ_frob(x2, Xm, T, q, d->p);
639 23884 : GEN lin = FpX_add(ZX_mul(gel(F,1), y2), ZX_mul(gel(F,2), x2), q);
640 23884 : return gc_upto(av, FpX_rem(lin, T, q));
641 : }
642 :
643 : static GEN
644 181258 : FpM_FpXV_bilinear(GEN P, GEN X, GEN Y, GEN p)
645 : {
646 181258 : pari_sp av = avma;
647 181258 : GEN s = ZX_mul(FpXV_FpC_mul(X,gel(P,1),p),gel(Y,1));
648 181258 : long i, l = lg(P);
649 861854 : for(i=2; i<l; i++)
650 680596 : s = ZX_add(s, ZX_mul(FpXV_FpC_mul(X,gel(P,i),p),gel(Y,i)));
651 181258 : return gc_upto(av, FpX_red(s, p));
652 : }
653 :
654 : static GEN
655 181258 : FpM_FpXQV_bilinear(GEN P, GEN X, GEN Y, GEN T, GEN p)
656 181258 : { return FpX_rem(FpM_FpXV_bilinear(P,X,Y,p),T,p); }
657 :
658 : static GEN
659 120862 : FpXC_powderiv(GEN M, GEN p)
660 : {
661 : long i, l;
662 120862 : long v = varn(gel(M,2));
663 120862 : GEN m = cgetg_copy(M, &l);
664 120862 : gel(m,1) = pol_0(v);
665 120862 : gel(m,2) = pol_1(v);
666 454328 : for(i=2; i<l-1; i++)
667 333466 : gel(m,i+1) = FpX_Fp_mul(gel(M,i),utoi(i), p);
668 120862 : return m;
669 : }
670 :
671 : struct _lift_iso
672 : {
673 : GEN phi, Xm, T, sqx, Tp;
674 : ulong p, pi;
675 : };
676 :
677 : static GEN
678 60396 : _lift_iter(void *E, GEN x2, GEN q)
679 : {
680 60396 : struct _lift_iso *d = (struct _lift_iso *) E;
681 60396 : ulong p = d->p;
682 60396 : long n = lg(d->phi)-2;
683 60396 : GEN TN = FpXT_red(d->T, q), XN = FpXV_red(d->Xm, q);
684 60396 : GEN y2 = ZpXQ_frob(x2, XN, TN, q, p);
685 60396 : GEN xp = FpXQ_powers(x2, n, TN, q);
686 60396 : GEN yp = FpXQ_powers(y2, n, TN, q);
687 60396 : GEN V = FpM_FpXQV_bilinear(d->phi,xp,yp,TN,q);
688 60396 : return mkvec3(V,xp,yp);
689 : }
690 :
691 : static GEN
692 60396 : _lift_invd(void *E, GEN V, GEN v, GEN qM, long M)
693 : {
694 60396 : struct _lift_iso *d = (struct _lift_iso *) E;
695 : struct _lift_lin e;
696 60396 : ulong p = d->p, pi = d->pi;
697 60396 : GEN TM = FpXT_red(d->T, qM), XM = FpXV_red(d->Xm, qM);
698 60396 : GEN xp = FpXV_red(gel(v,2), qM);
699 60396 : GEN yp = FpXV_red(gel(v,3), qM);
700 60396 : GEN Dx = FpM_FpXQV_bilinear(d->phi, FpXC_powderiv(xp, qM), yp, TM, qM);
701 60396 : GEN Dy = FpM_FpXQV_bilinear(d->phi, xp, FpXC_powderiv(yp, qM), TM, qM);
702 60396 : GEN F = mkvec4(Dy, Dx, TM, XM);
703 60396 : e.ai = Flxq_inv_pre(ZX_to_Flx(Dy,p),d->Tp, p, pi);
704 60396 : e.sqx = d->sqx; e.Tp = d->Tp; e.p=p; e.pi=pi; e.Xm = XM;
705 60396 : return gen_ZpX_Dixon(F,V,qM,utoipos(p),M,(void*) &e, _lift_lin, _lift_invl);
706 : }
707 :
708 : static GEN
709 25067 : lift_isogeny(GEN phi, GEN x0, long n, GEN Xm, GEN T, GEN sqx, GEN Tp,
710 : ulong p, ulong pi)
711 : {
712 : struct _lift_iso d;
713 25067 : d.phi = phi; d.Xm = Xm; d.T = T;
714 25067 : d.sqx = sqx; d.Tp = Tp; d.p = p; d.pi = pi;
715 25067 : return gen_ZpX_Newton(x0, utoipos(p), n,(void*)&d, _lift_iter, _lift_invd);
716 : }
717 :
718 : static GEN
719 25032 : getc2(GEN act, GEN X, GEN T, GEN q, ulong p, long N)
720 : {
721 25032 : GEN A1 = RgV_to_RgX(gel(act,1),0), A2 = RgV_to_RgX(gel(act,2),0);
722 25032 : long n = brent_kung_optpow(maxss(degpol(A1),degpol(A2)),2,1);
723 25032 : GEN xp = FpXQ_powers(X,n,T,q);
724 25032 : GEN P = FpX_FpXQV_eval(A1, xp, T, q);
725 25032 : GEN Q = FpX_FpXQV_eval(A2, xp, T, q);
726 25032 : return ZpXQ_div(P, Q, T, q, utoipos(p), N);
727 : }
728 :
729 : struct _ZpXQ_norm
730 : {
731 : long n;
732 : GEN T, p;
733 : };
734 :
735 : static GEN
736 32823 : ZpXQ_norm_mul(void *E, GEN x, GEN y)
737 : {
738 32823 : struct _ZpXQ_norm *D = (struct _ZpXQ_norm*)E;
739 32823 : GEN P = gel(x,1), Q = gel(y,1);
740 32823 : long a = mael(x,2,1), b = mael(y,2,1);
741 32823 : retmkvec2(FpXQ_mul(P,ZpXQ_frob_cyc(Q, D->T, D->p, a), D->T, D->p),
742 : mkvecsmall((a*b)%D->n));
743 : }
744 : static GEN
745 22715 : ZpXQ_norm_sqr(void *E, GEN x) { return ZpXQ_norm_mul(E, x, x); }
746 :
747 : /* Assume T = Phi_(n) and n prime */
748 : GEN
749 11340 : ZpXQ_norm_pcyc(GEN x, GEN T, GEN q, GEN p)
750 : {
751 : GEN z;
752 : struct _ZpXQ_norm D;
753 11340 : long d = get_FpX_degree(T);
754 11340 : D.T = T; D.p = q; D.n = d+1;
755 11340 : if (d==1) return ZX_copy(x);
756 11340 : z = mkvec2(x,mkvecsmall(p[2]));
757 11340 : z = gen_powu_i(z,d,(void*)&D,ZpXQ_norm_sqr,ZpXQ_norm_mul);
758 11340 : return gmael(z,1,2);
759 : }
760 :
761 : /* Assume T = Phi_(n) and n prime */
762 : static GEN
763 11102 : ZpXQ_sqrtnorm_pcyc(GEN x, GEN T, GEN q, GEN p, long e)
764 : {
765 11102 : GEN z = ZpXQ_norm_pcyc(x, T, q, p);
766 11102 : return Zp_sqrtlift(z,Fp_sqrt(z,p),p,e);
767 : }
768 :
769 : /* Assume a = 1 [p], return the square root of the norm */
770 : static GEN
771 13965 : ZpXQ_sqrtnorm(GEN a, GEN T, GEN q, GEN p, long e)
772 : {
773 13965 : GEN s = Fp_div(FpXQ_trace(ZpXQ_log(a, T, p, e), T, q), gen_2, q);
774 13965 : s = Qp_exp(cvtop(s, p, e-1)); return modii(padic_u(s), q);
775 : }
776 :
777 : struct _teich_lin
778 : {
779 : ulong p, pi;
780 : GEN sqx, Tp;
781 : long m;
782 : };
783 :
784 : static GEN
785 29701 : _teich_invl(void *E, GEN x)
786 : {
787 29701 : struct _teich_lin *d = (struct _teich_lin *) E;
788 29701 : ulong p = d->p, pi = d->pi;
789 29701 : return Flx_to_ZX(Flxq_lroot_fast_pre(ZX_to_Flx(x,p), d->sqx, d->Tp, p, pi));
790 : }
791 :
792 : static GEN
793 9079 : _teich_lin(void *E, GEN F, GEN x2, GEN q)
794 : {
795 9079 : struct _teich_lin *d = (struct _teich_lin *) E;
796 9079 : pari_sp av = avma;
797 9079 : GEN T = gel(F,2), Xm = gel(F,3);
798 9079 : GEN y2 = ZpXQ_frob(x2, Xm, T, q, d->p);
799 9079 : GEN lin = FpX_sub(y2, ZX_mulu(ZX_mul(gel(F,1), x2), d->p), q);
800 9079 : return gc_upto(av, FpX_rem(lin, T, q));
801 : }
802 :
803 : struct _teich_iso
804 : {
805 : GEN Xm, T, sqx, Tp;
806 : ulong p, pi;
807 : };
808 :
809 : static GEN
810 20622 : _teich_iter(void *E, GEN x2, GEN q)
811 : {
812 20622 : struct _teich_iso *d = (struct _teich_iso *) E;
813 20622 : ulong p = d->p;
814 20622 : GEN TN = FpXT_red(d->T, q), XN = FpXV_red(d->Xm, q);
815 20622 : GEN y2 = ZpXQ_frob(x2, XN, TN, q, d->p);
816 20622 : GEN x1 = FpXQ_powu(x2, p-1, TN, q);
817 20622 : GEN xp = FpXQ_mul(x2, x1, TN, q);
818 20622 : GEN V = FpX_sub(y2,xp,q);
819 20622 : return mkvec2(V,x1);
820 : }
821 :
822 : static GEN
823 20622 : _teich_invd(void *E, GEN V, GEN v, GEN qM, long M)
824 : {
825 20622 : struct _teich_iso *d = (struct _teich_iso *) E;
826 : struct _teich_lin e;
827 20622 : ulong p = d->p;
828 20622 : GEN TM = FpXT_red(d->T, qM), XM = FpXV_red(d->Xm, qM);
829 20622 : GEN x1 = FpX_red(gel(v,2), qM);
830 20622 : GEN F = mkvec3(x1, TM, XM);
831 20622 : e.sqx = d->sqx; e.Tp = d->Tp; e.p = p; e.pi = d->pi;
832 20622 : return gen_ZpX_Dixon(F,V,qM,utoipos(p),M,(void*) &e, _teich_lin, _teich_invl);
833 : }
834 :
835 : static GEN
836 10248 : Teichmuller_lift(GEN x, GEN Xm, GEN T, GEN sqx, GEN Tp, ulong p, ulong pi,
837 : long N)
838 : {
839 : struct _teich_iso d;
840 10248 : d.Xm = Xm; d.T = T; d.sqx = sqx; d.Tp = Tp; d.p = p; d.pi = pi;
841 10248 : return gen_ZpX_Newton(x,utoipos(p), N,(void*)&d, _teich_iter, _teich_invd);
842 : }
843 :
844 : static GEN
845 25067 : get_norm(GEN a4, GEN a6, GEN T, ulong p, ulong pi, long N)
846 : {
847 25067 : long sv=T[1];
848 : GEN a;
849 25067 : if (p==3) a = gel(a4,1);
850 : else
851 : {
852 10262 : GEN P = mkpoln(4, pol1_Flx(sv), pol0_Flx(sv), a4, a6);
853 10262 : a = gel(FlxqX_powu_pre(P, p>>1, T,p,pi), 2+p-1);
854 : }
855 25067 : return Zp_sqrtnlift(gen_1,subss(p,1),utoi(Flxq_norm(a,T,p)),utoipos(p), N);
856 : }
857 :
858 : static GEN
859 25032 : fill_pols(long n, const long *v, long m, const long *vn,
860 : const long *vd, GEN *act)
861 : {
862 : long i, j;
863 25032 : long d = upowuu(n,12/(n-1));
864 25032 : GEN N, D, M = zeromatcopy(n+1,n+1);
865 25032 : gmael(M,1,n+1) = gen_1;
866 120764 : for (i = 2; i <= n+1; i++)
867 339486 : for (j = i-1; j <= n; j++) gmael(M,i,j) = mulis(powuu(d,i-2), v[j-i+1]);
868 25032 : N = cgetg(m+1,t_COL);
869 25032 : D = cgetg(m+1,t_COL);
870 135541 : for(i = 1; i <= m; i++)
871 : {
872 110509 : gel(N,i) = stoi(*vn++);
873 110509 : gel(D,i) = stoi(*vd++);
874 : }
875 25032 : *act = mkmat2(N,D); return M;
876 : }
877 :
878 : /* These polynomials were extracted from the ECHIDNA databases
879 : * available at <http://echidna.maths.usyd.edu.au/echidna/>
880 : * and computed by David R. Kohel.
881 : * Return the matrix of the modular polynomial, set act to the parametrization,
882 : * and set dj to the opposite of the supersingular j-invariant. */
883 : static GEN
884 25032 : get_Kohel_polynomials(ulong p, GEN *act, long *dj)
885 : {
886 25032 : const long mat3[] = {-1,-36,-270};
887 25032 : const long num3[] = {1,-483,-21141,-59049};
888 25032 : const long den3[] = {1,261, 4347, -6561};
889 25032 : const long mat5[] = {-1,-30,-315,-1300,-1575};
890 25032 : const long num5[] = {-1,490,20620,158750,78125};
891 25032 : const long den5[] = {-1,-254,-4124,-12250,3125};
892 25032 : const long mat7[] = {-1,-28,-322,-1904,-5915,-8624,-4018};
893 25032 : const long num7[] = {1,-485,-24058,-343833,-2021642,-4353013,-823543};
894 25032 : const long den7[] = {1,259,5894,49119,168406,166355,-16807};
895 25032 : const long mat13[]= {-1,-26,-325,-2548,-13832,-54340,-157118,-333580,-509366,
896 : -534820,-354536,-124852,-15145};
897 25032 : const long num13[]= {1,-487,-24056,-391463,-3396483,-18047328,-61622301,
898 : -133245853,-168395656,-95422301,-4826809};
899 25032 : const long den13[]= {1,257,5896,60649,364629,1388256,3396483,5089019,4065464,
900 : 1069939,-28561};
901 25032 : switch(p)
902 : {
903 14805 : case 3:
904 14805 : *dj = 0;
905 14805 : return fill_pols(3,mat3,4,num3,den3,act);
906 10178 : case 5:
907 10178 : *dj = 0;
908 10178 : return fill_pols(5,mat5,5,num5,den5,act);
909 35 : case 7:
910 35 : *dj = 1;
911 35 : return fill_pols(7,mat7,7,num7,den7,act);
912 14 : case 13:
913 14 : *dj = 8;
914 14 : return fill_pols(13,mat13,11,num13,den13,act);
915 : }
916 : *dj=0; *act = NULL; return NULL; /* LCOV_EXCL_LINE */
917 : }
918 :
919 : long
920 25067 : zx_is_pcyc(GEN T)
921 : {
922 25067 : long i, n = degpol(T);
923 25067 : if (!uisprime(n+1)) return 0;
924 95263 : for (i = 0; i <= n; i++)
925 84161 : if (T[i+2] != 1UL) return 0;
926 11102 : return 1;
927 : }
928 :
929 : static GEN
930 25032 : Flxq_ellcard_Kohel(GEN a4, GEN a6, GEN T, ulong p)
931 : {
932 25032 : pari_sp av = avma, av2;
933 : pari_timer ti;
934 25032 : long n = get_Flx_degree(T), N = (n+4)/2, dj;
935 25032 : GEN q = powuu(p, N);
936 : GEN T2, Xm, s1, c2, t, lr, S1, sqx, Nc2, Np;
937 25032 : GEN act, phi = get_Kohel_polynomials(p, &act, &dj);
938 25032 : long ispcyc = zx_is_pcyc(get_Flx_mod(T));
939 25032 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
940 25032 : timer_start(&ti);
941 25032 : if (!ispcyc)
942 : {
943 13937 : T2 = Flx_Teichmuller(get_Flx_mod(T),p,N);
944 13937 : if (DEBUGLEVEL) timer_printf(&ti,"Teich");
945 : } else
946 11095 : T2 = Flx_to_ZX(get_Flx_mod(T));
947 :
948 25032 : T2 = FpX_get_red(T2, q); T = ZXT_to_FlxT(T2, p);
949 25032 : av2 = avma;
950 25032 : if (DEBUGLEVEL) timer_printf(&ti,"Barrett");
951 25032 : if (!ispcyc)
952 : {
953 13937 : Xm = FpXQ_powers(pol_xn(n,get_FpX_var(T2)),p-1,T2,q);
954 13937 : if (DEBUGLEVEL) timer_printf(&ti,"Xm");
955 : } else
956 11095 : Xm = cgetg(1,t_VEC);
957 25032 : s1 = Flxq_inv_pre(Flx_Fl_add(Flxq_ellj(a4,a6,T,p),dj, p),T,p,pi);
958 25032 : lr = Flxq_lroot_pre(polx_Flx(get_Flx_var(T)), T,p,pi);
959 25032 : sqx = Flxq_powers_pre(lr, p-1, T, p, pi);
960 25032 : S1 = lift_isogeny(phi, Flx_to_ZX(s1), N, Xm, T2, sqx, T,p,pi);
961 25032 : if (DEBUGLEVEL) timer_printf(&ti,"Lift isogeny");
962 25032 : c2 = getc2(act, S1, T2, q, p, N);
963 25032 : if (DEBUGLEVEL) timer_printf(&ti,"c^2");
964 25032 : if (p>3 && !ispcyc)
965 : {
966 10220 : GEN c2p = Flx_to_ZX(Flxq_inv_pre(ZX_to_Flx(c2,p),T,p,pi));
967 10220 : GEN tc2 = Teichmuller_lift(c2p,Xm, T2,sqx,T,p,pi,N);
968 10220 : if (DEBUGLEVEL) timer_printf(&ti,"Teichmuller/Fq");
969 10220 : c2 = FpX_rem(FpX_mul(tc2,c2,q),T2,q);
970 : }
971 25032 : c2 = gc_upto(av2, c2);
972 25032 : if (DEBUGLEVEL) timer_printf(&ti,"tc2");
973 25032 : Nc2 = (ispcyc? ZpXQ_sqrtnorm_pcyc: ZpXQ_sqrtnorm)(c2, T2, q, utoipos(p), N);
974 25032 : if (DEBUGLEVEL) timer_printf(&ti,"Norm");
975 25032 : Np = get_norm(a4,a6,T,p,pi,N);
976 25032 : if (p>3 && ispcyc)
977 : {
978 7 : GEN Ncpi = utoi(Fl_inv(umodiu(Nc2,p), p));
979 7 : GEN tNc2 = Zp_sqrtnlift(gen_1, subss(p,1), Ncpi, utoipos(p),N);
980 7 : if (DEBUGLEVEL) timer_printf(&ti,"Teichmuller/Fp");
981 7 : Nc2 = Fp_mul(Nc2,tNc2,q);
982 : }
983 25032 : t = Fp_center_i(Fp_mul(Nc2,Np,q),q,shifti(q,-1));
984 25032 : return gc_upto(av, subii(addiu(powuu(p,n),1),t));
985 : }
986 :
987 : /* Use Damien Robert's method */
988 : static GEN
989 35 : get_trace_Robert(GEN J, GEN phi, GEN Xm, GEN T, GEN q, ulong p, long e)
990 : {
991 35 : long n = lg(phi)-2;
992 35 : GEN K = ZpXQ_frob(J, Xm, T, q, p);
993 35 : GEN Jp = FpXQ_powers(J, n, T, q);
994 35 : GEN Kp = FpXQ_powers(K, n, T, q);
995 35 : GEN Jd = FpXC_powderiv(Jp, q);
996 35 : GEN Kd = FpXC_powderiv(Kp, q);
997 35 : GEN Dx = FpM_FpXQV_bilinear(phi, Kd, Jp, T, q);
998 35 : GEN Dy = FpM_FpXQV_bilinear(phi, Kp, Jd, T, q);
999 35 : GEN C = ZpXQ_inv(ZX_divuexact(Dy, p), T, utoi(p), e);
1000 35 : return FpX_neg(FpXQ_mul(Dx, C, T, q), q);
1001 : }
1002 :
1003 : /* in p^2, so p is tiny */
1004 : static GEN
1005 35 : Flxq_ellcard_Harley(GEN a4, GEN a6, GEN T, ulong p)
1006 : {
1007 35 : pari_sp av = avma, av2;
1008 : pari_timer ti;
1009 35 : long n = get_Flx_degree(T), N = (n+5)/2;
1010 35 : GEN pp = utoipos(p), q = powuu(p, N);
1011 : GEN T2, j, t, phi, J1, sqx, Xm, c2, tc2, c2p, Nc2, Np;
1012 35 : long ispcyc = zx_is_pcyc(get_Flx_mod(T));
1013 35 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p); /* = 0 here */
1014 35 : timer_start(&ti);
1015 35 : if (!ispcyc)
1016 : {
1017 28 : T2 = Flx_Teichmuller(get_Flx_mod(T),p,N);
1018 28 : if (DEBUGLEVEL) timer_printf(&ti,"Teich");
1019 : } else
1020 7 : T2 = Flx_to_ZX(get_Flx_mod(T));
1021 35 : T2 = FpX_get_red(T2, q); T = ZXT_to_FlxT(T2, p);
1022 35 : av2 = avma;
1023 35 : if (DEBUGLEVEL) timer_printf(&ti,"Barrett");
1024 35 : if (!ispcyc)
1025 : {
1026 28 : Xm = FpXQ_powers(pol_xn(n,get_FpX_var(T2)),p-1,T2,q);
1027 28 : if (DEBUGLEVEL) timer_printf(&ti,"Xm");
1028 : } else
1029 7 : Xm = cgetg(1,t_VEC);
1030 35 : j = Flxq_ellj(a4,a6,T,p);
1031 35 : sqx = Flxq_powers_pre(Flxq_lroot_pre(polx_Flx(T[1]), T,p,pi), p-1, T,p,pi);
1032 35 : phi = polmodular_ZM(p, 0);
1033 35 : if (DEBUGLEVEL) timer_printf(&ti,"phi");
1034 35 : J1 = lift_isogeny(phi, Flx_to_ZX(j), N, Xm, T2,sqx,T,p,pi);
1035 35 : if (DEBUGLEVEL) timer_printf(&ti,"Lift isogeny");
1036 35 : c2 = get_trace_Robert(J1, phi, Xm, T2, q, p, N);
1037 35 : q = diviuexact(q,p); N--;
1038 35 : if (DEBUGLEVEL) timer_printf(&ti,"c^2");
1039 35 : if (!ispcyc)
1040 : {
1041 28 : c2p = Flx_to_ZX(Flxq_inv_pre(ZX_to_Flx(c2,p),T,p,pi));
1042 28 : tc2 = Teichmuller_lift(c2p,Xm, T2,sqx,T,p,pi,N);
1043 28 : if (DEBUGLEVEL) timer_printf(&ti,"teichmuller");
1044 28 : c2 = FpX_rem(FpX_mul(tc2,c2,q),T2,q);
1045 : }
1046 35 : c2 = gc_upto(av2, c2);
1047 35 : q = powuu(p, N);
1048 35 : Nc2 = (ispcyc? ZpXQ_sqrtnorm_pcyc: ZpXQ_sqrtnorm)(c2, T2, q, pp, N);
1049 35 : if (DEBUGLEVEL) timer_printf(&ti,"Norm");
1050 35 : Np = get_norm(a4,a6,T,p,pi,N);
1051 35 : if (ispcyc)
1052 : {
1053 7 : GEN Ncpi = utoi(Fl_inv(umodiu(Nc2,p), p));
1054 7 : GEN tNc2 = Zp_sqrtnlift(gen_1, subss(p,1), Ncpi, pp, N);
1055 7 : if (DEBUGLEVEL) timer_printf(&ti,"Teichmuller/Fp");
1056 7 : Nc2 = Fp_mul(Nc2,tNc2,q);
1057 : }
1058 35 : t = Fp_center_i(Fp_mul(Nc2,Np,q),q,shifti(q,-1));
1059 35 : return gc_upto(av, subii(addiu(powuu(p,n),1),t));
1060 : }
1061 :
1062 : /***************************************************************************/
1063 : /* Shanks-Mestre */
1064 : /***************************************************************************/
1065 :
1066 : /* Return the lift of a (mod b), which is closest to h */
1067 : static GEN
1068 10317 : closest_lift(GEN a, GEN b, GEN h)
1069 10317 : { return addii(a, mulii(b, diviiround(subii(h,a), b))); }
1070 :
1071 : /* find multiple of order of f using Baby Step/Giant Step, f^h close to 1,
1072 : * order lies in an interval of size <= 'bound' and known mod B */
1073 : static GEN
1074 5543 : _FlxqE_order_multiple(void *E, GEN f, GEN h, GEN bound, GEN B)
1075 : {
1076 5543 : pari_sp av = avma, av1;
1077 : pari_timer Ti;
1078 5543 : long i, s = ceilsqrtdiv(bound, B) >> 1;
1079 : GEN P, F, tx, ti, fg, fh;
1080 :
1081 5543 : P = fh = _FlxqE_mul(E, f, h);
1082 5543 : if (DEBUGLEVEL >= 6) timer_start(&Ti);
1083 5543 : if (ell_is_inf(fh)) return h;
1084 5146 : F = _FlxqE_mul(E, f, B);
1085 5146 : if (s < 3)
1086 : { /* we're nearly done: naive search */
1087 1085 : GEN Q = P;
1088 1085 : for (i=1;; i++)
1089 : {
1090 3574 : P = _FlxqE_add(E, P, F); /* h.f + i.F */
1091 3574 : if (ell_is_inf(P)) return gc_upto(av, addii(h, mului(i,B)));
1092 3010 : Q = _FlxqE_sub(E, Q, F); /* h.f - i.F */
1093 3010 : if (ell_is_inf(Q)) return gc_upto(av, subii(h, mului(i,B)));
1094 : }
1095 : }
1096 4061 : tx = cgetg(s+1,t_VECSMALL); av1 = avma;
1097 43335 : for (i=1; i<=s; i++)
1098 : { /* baby steps */
1099 39603 : tx[i] = hash_GEN(gel(P, 1));
1100 39603 : P = _FlxqE_add(E, P, F); /* h.f + i.F */
1101 39603 : if (ell_is_inf(P)) return gc_upto(av, addii(h, mului(i,B)));
1102 39274 : if (gc_needed(av1,3))
1103 : {
1104 0 : if(DEBUGMEM>1) pari_warn(warnmem,"[Flxq_ellcard] baby steps, i=%ld",i);
1105 0 : P = gc_upto(av1,P);
1106 : }
1107 : }
1108 3732 : if (DEBUGLEVEL >= 6) timer_printf(&Ti,"[Flxq_ellcard] baby steps, s = %ld",s);
1109 : /* giant steps: fg = s.F */
1110 3732 : fg = gc_upto(av1, _FlxqE_sub(E, P, fh));
1111 3732 : if (ell_is_inf(fg)) return gc_upto(av, mului(s,B));
1112 3732 : ti = vecsmall_indexsort(tx); /* = permutation sorting tx */
1113 3732 : tx = perm_mul(tx,ti);
1114 3732 : if (DEBUGLEVEL >= 6) timer_printf(&Ti, "[Flxq_ellcard] sorting");
1115 3732 : av1 = avma;
1116 3732 : for (P=fg, i=1; ; i++)
1117 35556 : {
1118 39288 : long k = hash_GEN(gel(P,1)), r = zv_search(tx, k);
1119 39288 : if (r)
1120 : {
1121 7467 : while (r && tx[r] == k) r--;
1122 3732 : for (r++; r <= s && tx[r] == k; r++)
1123 : {
1124 3732 : long j = ti[r]-1;
1125 3732 : GEN Q = _FlxqE_add(E, _FlxqE_mul(E, F, stoi(j)), fh);
1126 3732 : if (DEBUGLEVEL >= 6)
1127 0 : timer_printf(&Ti, "[Flxq_ellcard] giant steps, i = %ld",i);
1128 3732 : if (Flx_equal(gel(P,1), gel(Q,1)))
1129 : {
1130 3732 : if (Flx_equal(gel(P,2), gel(Q,2))) i = -i;
1131 3732 : return gc_upto(av,addii(h, mulii(addis(mulss(s,i), j), B)));
1132 : }
1133 : }
1134 : }
1135 35556 : P = _FlxqE_add(E, P, fg);
1136 35556 : if (gc_needed(av1,3))
1137 : {
1138 0 : if(DEBUGMEM>1) pari_warn(warnmem,"[Flxq_ellcard] giants steps, i=%ld",i);
1139 0 : P = gc_upto(av1,P);
1140 : }
1141 : }
1142 : }
1143 : static GEN
1144 5543 : _FlxqE_order(void *E, GEN f, GEN h, GEN bound, GEN B)
1145 : {
1146 5543 : GEN o = _FlxqE_order_multiple(E, f, h, bound, B);
1147 5543 : return gen_order(f, o, E, &FlxqE_group);
1148 : }
1149 :
1150 : static void
1151 64512 : Flx_next(GEN t, ulong p)
1152 : {
1153 : long i;
1154 78554 : for(i=2;;i++)
1155 78554 : if (uel(t,i)==p-1) t[i]=0; else { t[i]++; break; }
1156 64512 : }
1157 :
1158 : static void
1159 64512 : Flx_renormalize_ip(GEN x, long lx)
1160 : {
1161 : long i;
1162 78554 : for (i = lx-1; i>=2; i--)
1163 71414 : if (x[i]) break;
1164 64512 : setlg(x, i+1);
1165 64512 : }
1166 :
1167 : static ulong
1168 6237 : F3xq_ellcard_naive(GEN a2, GEN a6, GEN T)
1169 : {
1170 6237 : pari_sp av = avma;
1171 6237 : long i, d = get_Flx_degree(T), lx = d+2;
1172 6237 : long q = upowuu(3, d), a;
1173 6237 : GEN x = zero_zv(lx); x[1] = get_Flx_var(T);
1174 27174 : for(a=1, i=0; i<q; i++)
1175 : {
1176 : GEN rhs;
1177 20937 : Flx_renormalize_ip(x, lx);
1178 20937 : rhs = Flx_add(Flxq_mul(Flxq_sqr(x, T, 3), Flx_add(x, a2, 3), T, 3), a6, 3);
1179 20937 : if (!lgpol(rhs)) a++; else if (Flxq_issquare(rhs, T, 3)) a+=2;
1180 20937 : Flx_next(x, 3);
1181 : }
1182 6237 : set_avma(av); return a;
1183 : }
1184 :
1185 : /* p^deg(T) is tiny */
1186 : static ulong
1187 903 : Flxq_ellcard_naive(GEN a4, GEN a6, GEN T, ulong p)
1188 : {
1189 903 : pari_sp av = avma;
1190 903 : long i, d = get_Flx_degree(T), lx = d+2;
1191 903 : long q = upowuu(p, d), a;
1192 903 : GEN x = zero_zv(lx); x[1] = get_Flx_var(T);
1193 44478 : for(a = 1, i = 0; i < q; i++)
1194 : {
1195 : GEN x2, rhs;
1196 43575 : Flx_renormalize_ip(x, lx);
1197 43575 : x2 = Flxq_sqr_pre(x, T, p, 0);
1198 43575 : rhs = Flx_add(Flxq_mul_pre(x, Flx_add(x2, a4, p), T, p, 0), a6, p);
1199 43575 : if (!lgpol(rhs)) a++; else if (Flxq_issquare(rhs,T,p)) a += 2;
1200 43575 : Flx_next(x,p);
1201 : }
1202 903 : set_avma(av); return a;
1203 : }
1204 :
1205 : static long
1206 11240 : Flxq_kronecker(GEN x, GEN T, ulong p)
1207 : {
1208 : pari_sp av;
1209 11240 : if (lgpol(x) == 0) return 0;
1210 11213 : av = avma; return gc_long(av, krouu(Flxq_norm(x, T, p), p));
1211 : }
1212 :
1213 : /* Find x such that kronecker(u = x^3+a4x+a6, p) is KRO.
1214 : * Return point [x*u,u^2] on E (KRO=1) / E^twist (KRO=-1) */
1215 : static GEN
1216 5543 : Flxq_ellpoint(long KRO, GEN a4, GEN a6, GEN T, ulong p, ulong pi)
1217 : {
1218 5543 : long v = get_Flx_var(T), n = get_Flx_degree(T);
1219 : for(;;)
1220 5697 : {
1221 11240 : GEN x = random_Flx(n, v, p), x2 = Flxq_sqr_pre(x,T,p,pi);
1222 11240 : GEN u = Flx_add(a6, Flxq_mul_pre(Flx_add(a4, x2, p), x, T,p, pi), p);
1223 11240 : if (Flxq_kronecker(u,T,p) == KRO)
1224 5543 : return mkvec2(Flxq_mul_pre(u,x, T,p,pi), Flxq_sqr_pre(u, T,p,pi));
1225 : }
1226 : }
1227 :
1228 : static GEN
1229 4774 : Flxq_ellcard_Shanks(GEN a4, GEN a6, GEN q, GEN T, ulong p)
1230 : {
1231 4774 : pari_sp av = avma;
1232 4774 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
1233 4774 : long v = get_Flx_var(T), KRO = -1;
1234 : GEN h,f, A, B;
1235 4774 : GEN q1p = addiu(q,1), q2p = shifti(q1p, 1);
1236 4774 : GEN bound = addiu(sqrti(gmul2n(q,4)), 1); /* ceil( 4sqrt(q) ) */
1237 : struct _FlxqE e;
1238 4774 : e.p = p; e.pi = pi; e.T = Flx_get_red_pre(T, p, pi);
1239 : /* once #E(Flxq) is known mod B >= bound, it is determined */
1240 4774 : switch(FlxqX_nbroots(mkpoln(4, pol1_Flx(v), pol0_Flx(v), a4, a6), T, p))
1241 : { /* how many 2-torsion points ? */
1242 2891 : case 3: A = gen_0; B = utoipos(4); break;
1243 1414 : case 1: A = gen_0; B = gen_2; break;
1244 469 : default: A = gen_1; B = gen_2; break; /* 0 */
1245 : }
1246 : for(;;)
1247 : {
1248 5543 : h = closest_lift(A, B, q1p);
1249 : /* [ux, u^2] is on E_u: y^2 = x^3 + c4 u^2 x + c6 u^3
1250 : * E_u isomorphic to E (resp. E') iff KRO = 1 (resp. -1)
1251 : * #E(F_p) = p+1 - a_p, #E'(F_p) = p+1 + a_p
1252 : *
1253 : * #E_u(Flxq) = A (mod B), h is close to #E_u(Flxq) */
1254 5543 : KRO = -KRO;
1255 5543 : f = Flxq_ellpoint(KRO, a4,a6, T,p,pi);
1256 5543 : e.a4 = Flxq_mul_pre(a4, gel(f,2), T,p,pi); /* a4 for E_u */
1257 5543 : h = _FlxqE_order((void*)&e, f, h, bound, B);
1258 : /* h | #E_u(Flxq) = A (mod B) */
1259 5543 : A = Z_chinese_all(A, gen_0, B, h, &B);
1260 5543 : if (cmpii(B, bound) >= 0) break;
1261 : /* not done, update A mod B for the _next_ curve, isomorphic to
1262 : * the quadratic twist of this one */
1263 769 : A = remii(subii(q2p,A), B); /* #E(Fq)+#E'(Fq) = 2q+2 */
1264 : }
1265 4774 : h = closest_lift(A, B, q1p);
1266 4774 : return gc_INT(av, KRO == 1? h: subii(q2p,h));
1267 : }
1268 :
1269 : static GEN
1270 21042 : F3xq_ellcard(GEN a2, GEN a6, GEN T)
1271 : {
1272 21042 : long n = get_Flx_degree(T);
1273 21042 : if (n <= 2)
1274 5936 : return utoi(F3xq_ellcard_naive(a2, a6, T));
1275 : else
1276 : {
1277 15106 : GEN q1 = addiu(powuu(3, get_Flx_degree(T)), 1), t;
1278 15106 : GEN a = Flxq_div(a6,Flxq_powu(a2,3,T,3),T,3);
1279 15106 : if (Flx_equal1(Flxq_powu(a, 8, T, 3)))
1280 : {
1281 301 : GEN P = Flxq_minpoly(a,T,3);
1282 301 : long dP = degpol(P); /* dP <= 2 */
1283 301 : ulong q = upowuu(3,dP);
1284 301 : GEN A2 = pol1_Flx(P[1]), A6 = Flx_rem(polx_Flx(P[1]), P, 3);
1285 301 : long tP = q + 1 - F3xq_ellcard_naive(A2, A6, P);
1286 301 : t = elltrace_extension(stoi(tP), n/dP, utoi(q));
1287 301 : if (umodiu(t, 3)!=1) t = negi(t);
1288 301 : return Flx_equal1(a2) || Flxq_issquare(a2,T,3) ? subii(q1,t): addii(q1,t);
1289 : }
1290 14805 : else return Flxq_ellcard_Kohel(mkvec(a2), a6, T, 3);
1291 : }
1292 : }
1293 :
1294 : static GEN
1295 11165 : Flxq_ellcard_Satoh(GEN a4, GEN a6, GEN j, GEN T, ulong p)
1296 : {
1297 11165 : long n = get_Flx_degree(T);
1298 11165 : if (n <= 2)
1299 903 : return utoi(Flxq_ellcard_naive(a4, a6, T, p));
1300 : else
1301 : {
1302 10262 : GEN jp = Flxq_powu(j, p, T, p);
1303 10262 : GEN s = Flx_add(j, jp, p);
1304 10262 : if (degpol(s) <= 0)
1305 : { /* it is assumed j not in F_p */
1306 0 : GEN m = Flxq_mul(j, jp, T, p);
1307 0 : if (degpol(m) <= 0)
1308 : {
1309 0 : GEN q = sqru(p);
1310 0 : GEN q1 = addiu(powuu(p, get_Flx_degree(T)), 1);
1311 0 : GEN sk = Flx_Fl_add(Flx_neg(j, p), 1728%p, p);
1312 0 : GEN sA4 = Flx_triple(Flxq_mul(sk, j, T, p), p);
1313 0 : GEN u = Flxq_div(a4, sA4, T, p);
1314 0 : ulong ns = lgpol(s) ? Fl_neg(s[2], p): 0UL;
1315 0 : GEN P = mkvecsmall4(T[1], m[2], ns, 1L);
1316 : GEN A4, A6, t, tP;
1317 0 : Flxq_ellj_to_a4a6(polx_Flx(T[1]), P, p, &A4, &A6);
1318 0 : tP = addis(q, 1 - Flxq_ellcard_naive(A4, A6, P, p));
1319 0 : t = elltrace_extension(tP, n>>1, q);
1320 0 : return Flxq_is2npower(u, 2, T, p) ? subii(q1,t): addii(q1,t);
1321 : }
1322 : }
1323 10262 : if (p<=7 || p==13) return Flxq_ellcard_Kohel(a4, a6, T, p);
1324 35 : else return Flxq_ellcard_Harley(a4, a6, T, p);
1325 : }
1326 : }
1327 :
1328 : #if 0
1329 : /* Not used, slower than alternative, can be simulated by hyperellcharpoly */
1330 : static GEN
1331 : Flxq_ellcard_Kedlaya(GEN a4, GEN a6, GEN T, ulong p)
1332 : {
1333 : pari_sp av = avma;
1334 : GEN H = mkpoln(4, gen_1, gen_0, Flx_to_ZX(a4), Flx_to_ZX(a6));
1335 : GEN Tp = Flx_to_ZX(get_Flx_mod(T));
1336 : long n = degpol(Tp), e = ((p < 16 ? n+1: n)>>1)+1;
1337 : GEN M = ZlXQX_hyperellpadicfrobenius(H, Tp, p, e);
1338 : GEN N = ZpXQM_prodFrobenius(M, Tp, utoipos(p), e);
1339 : GEN q = powuu(p, e);
1340 : GEN tp = Fq_add(gcoeff(N,1,1), gcoeff(N,2,2), Tp, q);
1341 : GEN t = Fp_center_i(typ(tp)==t_INT ? tp: leading_coeff(tp), q, shifti(q,-1));
1342 : return gc_upto(av, subii(addiu(powuu(p, n), 1), t));
1343 : }
1344 : #endif
1345 :
1346 : GEN
1347 57333 : Flxq_ellj(GEN a4, GEN a6, GEN T, ulong p)
1348 : {
1349 57333 : pari_sp av=avma;
1350 57333 : if (p==3)
1351 : {
1352 : GEN J;
1353 14805 : if (typ(a4)!=t_VEC) return pol0_Flx(get_Flx_var(T));
1354 14805 : J = Flxq_div(Flxq_powu(gel(a4,1),3, T, p),Flx_neg(a6,p), T, p);
1355 14805 : return gc_leaf(av, J);
1356 : }
1357 : else
1358 : {
1359 42528 : pari_sp av=avma;
1360 42528 : GEN a43 = Flxq_mul(a4,Flxq_sqr(a4,T,p),T,p);
1361 42528 : GEN a62 = Flxq_sqr(a6,T,p);
1362 42528 : GEN num = Flx_mulu(a43,6912,p);
1363 42528 : GEN den = Flx_add(Flx_mulu(a43,4,p),Flx_mulu(a62,27,p),p);
1364 42528 : return gc_leaf(av, Flxq_div(num, den, T, p));
1365 : }
1366 : }
1367 :
1368 : void
1369 0 : Flxq_ellj_to_a4a6(GEN j, GEN T, ulong p, GEN *pt_a4, GEN *pt_a6)
1370 : {
1371 0 : ulong zagier = 1728 % p;
1372 0 : if (lgpol(j)==0)
1373 0 : { *pt_a4 = pol0_Flx(T[1]); *pt_a6 =pol1_Flx(T[1]); }
1374 0 : else if (lgpol(j)==1 && uel(j,2) == zagier)
1375 0 : { *pt_a4 = pol1_Flx(T[1]); *pt_a6 =pol0_Flx(T[1]); }
1376 : else
1377 : {
1378 0 : GEN k = Flx_Fl_add(Flx_neg(j, p), zagier, p);
1379 0 : GEN kj = Flxq_mul(k, j, T, p);
1380 0 : GEN k2j = Flxq_mul(kj, k, T, p);
1381 0 : *pt_a4 = Flx_triple(kj, p);
1382 0 : *pt_a6 = Flx_double(k2j, p);
1383 : }
1384 0 : }
1385 :
1386 : static GEN
1387 9240 : F3xq_ellcardj(GEN a4, GEN a6, GEN T, GEN q, long n)
1388 : {
1389 9240 : const ulong p = 3;
1390 : ulong t;
1391 9240 : GEN q1 = addiu(q,1);
1392 9240 : GEN na4 = Flx_neg(a4,p), ra4;
1393 9240 : if (!Flxq_issquare(na4,T,p))
1394 4826 : return q1;
1395 4414 : ra4 = Flxq_sqrt(na4,T,p);
1396 4414 : t = Flxq_trace(Flxq_div(a6,Flxq_mul(na4,ra4,T,p),T,p),T,p);
1397 4414 : if (n%2==1)
1398 : {
1399 : GEN q3;
1400 2381 : if (t==0) return q1;
1401 840 : q3 = powuu(p,(n+1)>>1);
1402 840 : return (t==1)^(n%4==1) ? subii(q1,q3): addii(q1,q3);
1403 : }
1404 : else
1405 : {
1406 2033 : GEN q22, q2 = powuu(p,n>>1);
1407 2033 : GEN W = Flxq_pow(a4,shifti(q,-2),T,p);
1408 2033 : long s = (W[2]==1)^(n%4==2);
1409 2033 : if (t!=0) return s ? addii(q1,q2): subii(q1, q2);
1410 2033 : q22 = shifti(q2,1);
1411 2033 : return s ? subii(q1,q22): addii(q1, q22);
1412 : }
1413 : }
1414 :
1415 : static GEN
1416 15826 : Flxq_ellcardj(GEN a4, GEN a6, ulong j, GEN T, GEN q, ulong p, long n)
1417 : {
1418 15826 : GEN q1 = addiu(q,1);
1419 15826 : if (j==0)
1420 : {
1421 : ulong w;
1422 : GEN W, t, N;
1423 5677 : if (umodiu(q,6)!=1) return q1;
1424 4277 : N = Fp_ffellcard(gen_0,gen_1,q,n,utoipos(p));
1425 4277 : t = subii(q1, N);
1426 4277 : W = Flxq_pow(a6,diviuexact(shifti(q,-1), 3),T,p);
1427 4277 : if (degpol(W)>0) /*p=5 mod 6*/
1428 1542 : return Flx_equal1(Flxq_powu(W,3,T,p)) ? addii(q1,shifti(t,-1)):
1429 514 : subii(q1,shifti(t,-1));
1430 3249 : w = W[2];
1431 3249 : if (w==1) return N;
1432 2556 : if (w==p-1) return addii(q1,t);
1433 : else /*p=1 mod 6*/
1434 : {
1435 1814 : GEN u = shifti(t,-1), v = sqrtint(diviuexact(subii(q,sqri(u)),3));
1436 1814 : GEN a = addii(u,v), b = shifti(v,1);
1437 1814 : if (Fl_powu(w,3,p)==1)
1438 : {
1439 907 : if (Fl_add(umodiu(a,p),Fl_mul(w,umodiu(b,p),p),p)==0)
1440 452 : return subii(q1,subii(shifti(b,1),a));
1441 : else
1442 455 : return addii(q1,addii(a,b));
1443 : }
1444 : else
1445 : {
1446 907 : if (Fl_sub(umodiu(a,p),Fl_mul(w,umodiu(b,p),p),p)==0)
1447 452 : return subii(q1,subii(a,shifti(b,1)));
1448 : else
1449 455 : return subii(q1,addii(a,b));
1450 : }
1451 : }
1452 10149 : } else if (j==1728%p)
1453 : {
1454 : ulong w;
1455 : GEN W, N, t;
1456 5662 : if (mod4(q)==3) return q1;
1457 4262 : W = Flxq_pow(a4,shifti(q,-2),T,p);
1458 4262 : if (degpol(W)>0) return q1; /*p=3 mod 4*/
1459 3566 : w = W[2];
1460 3566 : N = Fp_ffellcard(gen_1,gen_0,q,n,utoipos(p));
1461 3566 : if(w==1) return N;
1462 2539 : t = subii(q1, N);
1463 2539 : if(w==p-1) return addii(q1, t);
1464 : else /*p=1 mod 4*/
1465 : {
1466 1450 : GEN u = shifti(t,-1), v = sqrtint(subii(q,sqri(u)));
1467 1450 : if (Fl_add(umodiu(u,p),Fl_mul(w,umodiu(v,p),p),p)==0)
1468 725 : return subii(q1,shifti(v,1));
1469 : else
1470 725 : return addii(q1,shifti(v,1));
1471 : }
1472 : } else
1473 : {
1474 4487 : ulong g = Fl_div(j, Fl_sub(1728%p, j, p), p);
1475 4487 : GEN N = Fp_ffellcard(utoi(Fl_triple(g,p)),utoi(Fl_double(g,p)),q,n,utoipos(p));
1476 4487 : GEN l = Flxq_mul(Flx_triple(a6,p),Flx_double(a4,p),T,p);
1477 4487 : if (Flxq_issquare(l,T,p)) return N;
1478 2940 : return subii(shifti(q1,1),N);
1479 : }
1480 : }
1481 :
1482 : static GEN
1483 454 : Flxq_ffellcard(GEN a4, GEN a6, GEN M, GEN q, GEN T, ulong p, long n)
1484 : {
1485 454 : long m = degpol(M);
1486 454 : GEN j = polx_Flx(M[1]);
1487 454 : GEN g = Flxq_div(j, mkvecsmall3(M[1],1728%p,p-1), M, p);
1488 454 : GEN N = Flxq_ellcard(Flx_triple(g, p), Flx_double(g, p), M, p);
1489 454 : GEN qm = powuu(p, m), q1 = addiu(q, 1), qm1 = addiu(qm, 1);
1490 454 : GEN l = Flxq_mul(Flx_triple(a6,p), Flx_double(a4,p), T, p);
1491 454 : GEN te = elltrace_extension(subii(qm1, N), n/m, qm);
1492 454 : return Flxq_issquare(l,T,p) ? subii(q1, te): addii(q1, te);
1493 : }
1494 :
1495 : static GEN
1496 62759 : Flxq_ellcard_i(GEN a4, GEN a6, GEN T, ulong p)
1497 : {
1498 62759 : long n = get_Flx_degree(T);
1499 62759 : GEN J, M, q = powuu(p, n);
1500 62759 : if (typ(a4)==t_VEC)
1501 21042 : return F3xq_ellcard(gel(a4,1), a6, T);
1502 41717 : if (p==3)
1503 9240 : return F3xq_ellcardj(a4, a6, T, q, n);
1504 32477 : if (degpol(a4)<=0 && degpol(a6)<=0)
1505 211 : return Fp_ffellcard(utoi(Flx_eval(a4,0,p)),utoi(Flx_eval(a6,0,p)),q,n,utoipos(p));
1506 32266 : J = Flxq_ellj(a4,a6,T,p);
1507 32266 : if (degpol(J)<=0)
1508 15826 : return Flxq_ellcardj(a4,a6,lgpol(J)?J[2]:0,T,q,p,n);
1509 16440 : M = Flxq_minpoly(J, T, p);
1510 16440 : if (degpol(M) < n)
1511 454 : return Flxq_ffellcard(a4, a6, M, q, T, p, n);
1512 15986 : if (p <= 7)
1513 10997 : return Flxq_ellcard_Satoh(a4, a6, J, T, p);
1514 4989 : if (cmpis(q,100)<0)
1515 0 : return utoi(Flxq_ellcard_naive(a4, a6, T, p));
1516 4989 : if (p == 13 || ((double) 6*p <= n*log(p)*log(p) && (BITS_IN_LONG==64 || p <= 103)))
1517 168 : return Flxq_ellcard_Satoh(a4, a6, J, T, p);
1518 : #if 0
1519 : /* No good case ? */
1520 : if (...)
1521 : return Flxq_ellcard_Kedlaya(a4, a6, T, p);
1522 : #endif
1523 4821 : if (expi(q)<=62)
1524 4774 : return Flxq_ellcard_Shanks(a4, a6, q, T, p);
1525 : else
1526 47 : return Fq_ellcard_SEA(Flx_to_ZX(a4),Flx_to_ZX(a6),q,Flx_to_ZX(T),utoipos(p),0);
1527 : }
1528 :
1529 : GEN
1530 62759 : Flxq_ellcard(GEN a4, GEN a6, GEN T, ulong p)
1531 : {
1532 62759 : pari_sp av = avma;
1533 62759 : return gc_INT(av, Flxq_ellcard_i(a4, a6, T, p));
1534 : }
1535 :
1536 : static long
1537 371 : Fl_ellj_trace(ulong j, ulong p)
1538 : {
1539 : ulong a4, a6;
1540 371 : Fl_ellj_to_a4a6(j, p, &a4, &a6);
1541 371 : return Fl_elltrace(a4, a6, p);
1542 : }
1543 :
1544 : /* Given ordinary E/Fq, a prime ell, and the height of the ell-volcano
1545 : * containing j(E) (= v_ell(conductor of Z[pi_E]) returns the height of j(E)
1546 : * on its ell-volcano (= v_ell(conductor of the order End(E)). */
1547 : static long
1548 154 : Fl_ellheightabovefloor(ulong j, long ell, long e, ulong p)
1549 : {
1550 154 : pari_sp av = avma;
1551 : GEN Xp, G, phi, phix, j0, j1;
1552 : long h, i, nj1;
1553 154 : if (e==0) return 0;
1554 154 : if (j==0 || j==1728%p) return e;
1555 126 : phi = ZXX_to_FlxX(polmodular_ZXX(ell, 0, 0, 1), p, 1);
1556 126 : phix = FlxY_evalx(phi, j, p);
1557 126 : Xp = Flx_Frobenius(phix, p);
1558 126 : G = Flx_gcd(Flx_sub(Xp, polx_Flx(0), p), phix, p);
1559 126 : nj1 = degpol(G);
1560 126 : if (nj1 < ell) return 0;
1561 112 : if (e==1 || nj1 != ell+1) return e;
1562 21 : j1 = Flx_roots(G, p);
1563 21 : nj1 = lg(j1)-1;
1564 21 : if (nj1 < 3) return 0;
1565 21 : j0 = mkvecsmall3(j,j,j);
1566 42 : for (h = 1; ; h++)
1567 133 : for(i = 1; i <= 3; i++)
1568 : {
1569 112 : GEN P = Flx_div_by_X_x(FlxY_evalx(phi, uel(j1,i), p), uel(j0,i), p, NULL);
1570 112 : GEN r = Flx_roots(P, p);
1571 112 : if (lg(r) == 1) return gc_long(av, h);
1572 91 : j0[i] = j1[i];
1573 91 : j1[i] = r[1];
1574 : }
1575 : }
1576 :
1577 : /* Given an ordinary elliptic curve E/Fp and an integer h, returns
1578 : * D = disc(End(E)) assuming h(D) = h, using the approach sketched in
1579 : * Remark 13. If the algorithm returns 0 it has proved that h(D) != h, but it
1580 : * is under no obligation to do so and is allowed to return any value when the
1581 : * assumption h(d) = h is false. */
1582 : static long
1583 371 : Fl_end13(ulong j, ulong h, ulong p)
1584 : {
1585 : ulong D0, v, h0;
1586 : long i, lL, lc, lD, nc;
1587 : GEN D, DF, cs, L, vP, vE;
1588 371 : ulong t = Fl_ellj_trace(j, p);
1589 :
1590 371 : D0 = coredisc2u_fact(factoru(4*p-t*t), -1, &vP, &vE);
1591 371 : h0 = itou(classno(stoi(-D0)));
1592 371 : if (h % h0 != 0) return 0;
1593 357 : h /= h0;
1594 357 : D = divisorsu_fact_factored(mkmat2(vP,vE));
1595 357 : DF = gel(D,2); D = gel(D,1);
1596 357 : lD = lg(D); v = D[lD-1];
1597 357 : cs = cgetg(lD,t_VECSMALL); nc = 0;
1598 1897 : for (i = 1; i < lD; i++)
1599 : {
1600 1540 : GEN F = gel(DF,i);
1601 1540 : ulong w = uquadclassnoF_fact(D0, -1, gel(F,1), gel(F,2));
1602 1540 : if (w == h) uel(cs,++nc) = v / uel(D,i);
1603 : }
1604 357 : if (nc==0) return 0;
1605 357 : if (nc==1) { v /= uel(cs,1); return -D0*v*v; }
1606 154 : L = cgetg(nc+1, t_VEC);
1607 469 : for (i = 1; i <= nc; i++) gel(L,i) = gel(factoru(uel(cs,i)), 1);
1608 154 : L = vecsmall_uniq(shallowconcat1(L));
1609 154 : lL = lg(L); lc = nc+1;
1610 154 : for (i = 1; i < lL; i++)
1611 : {
1612 154 : ulong ell = L[i];
1613 154 : long k, e = Fl_ellheightabovefloor(j, ell, z_lval(v,ell), p);
1614 469 : for (k = 1; k < lc; k++)
1615 315 : if(cs[k] && z_lval(cs[k], ell) != e) { cs[k] = 0; nc--; }
1616 154 : if (nc==0) return 0;
1617 154 : if (nc==1)
1618 : {
1619 175 : for (k = 1; k < lc; k++)
1620 175 : if (cs[k]) { v /= uel(cs,k); return -D0*v*v; }
1621 : }
1622 : }
1623 0 : return 0;
1624 : }
1625 :
1626 : INLINE int
1627 392 : RgX_is_monic_ZX(GEN pol)
1628 392 : { return RgX_is_ZX(pol) && ZX_is_monic(pol); }
1629 :
1630 : long
1631 399 : polisclass(GEN H)
1632 : {
1633 399 : pari_sp av = avma, btop;
1634 399 : long h = degpol(H), hl, i, pmin, vH = varn(H), vh;
1635 : double lmin;
1636 : ulong p;
1637 : GEN h2list;
1638 : forprime_t T;
1639 :
1640 399 : if (typ(H)!= t_POL) pari_err_TYPE("polsisclass",H);
1641 399 : if (h <= 0 || !RgX_is_monic_ZX(H)) return 0;
1642 371 : vh = vals(h);
1643 371 : h2list = cgetg(vh+2, t_VECSMALL); hl = 1;
1644 1155 : for (i = 0; i <= vh; i++)
1645 : {
1646 784 : ulong d = 1UL<<i;
1647 784 : if (((d-h)&1)==0) h2list[hl++] = d;
1648 : }
1649 371 : setlg(h2list, hl);
1650 371 : lmin = h * (log(log(h+2))+2);
1651 371 : pmin = 33 * ceil(lmin*lmin);
1652 371 : u_forprime_init(&T, pmin, ULONG_MAX);
1653 371 : btop = avma;
1654 3465 : while((p = u_forprime_next(&T)))
1655 : {
1656 : ulong r;
1657 : long D, nroots;
1658 3465 : GEN Xp, G, Hp = ZX_to_Flx(H,p);
1659 3465 : if (!Flx_is_squarefree(Hp, p)) { set_avma(btop); continue; }
1660 3465 : Xp = Flx_Frobenius(Hp, p);
1661 3465 : G = Flx_gcd(Flx_sub(Xp, polx_Flx(evalvarn(vH)), p), Hp, p);
1662 3465 : nroots = degpol(G);
1663 3465 : if (nroots==0) { set_avma(btop); continue; }
1664 1295 : if (nroots < h && !zv_search(h2list,nroots)) return gc_long(av, 0);
1665 1295 : r = Flx_oneroot(G, p);
1666 1295 : if (Fp_elljissupersingular(utoi(r), utoi(p))) { set_avma(btop); continue; }
1667 371 : D = Fl_end13(r, h, p);
1668 371 : if (D && gequal(H, polclass(stoi(D), 0, vH))) return gc_long(av, D);
1669 21 : return gc_long(av, 0);
1670 : }
1671 0 : pari_err_BUG("polisclass");
1672 : return 0; /* LCOV_EXCL_LINE */
1673 : }
|
