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 gc_upto(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 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 30443 : set_avma(av);
93 30443 : a = random_Flx(n, vs, p);
94 30443 : } 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 1607671 : FlxqE_dbl_slope(GEN P, GEN a4, GEN T, ulong p, ulong pi, GEN *ps)
117 : {
118 : GEN x, y, Q, s;
119 1607671 : if (ell_is_inf(P) || !lgpol(gel(P,2))) return ellinf();
120 1481075 : x = gel(P,1); y = gel(P,2);
121 1481075 : if (p==3UL)
122 528465 : s = typ(a4)==t_VEC? Flxq_div_pre(Flxq_mul_pre(x, gel(a4,1), T,p,pi), y, T,p,pi)
123 528465 : : Flxq_div_pre(a4, Flx_neg(y, p), T,p,pi);
124 : else
125 : {
126 952610 : GEN sx = Flx_add(Flx_triple(Flxq_sqr_pre(x, T, p, pi), p), a4, p);
127 952610 : s = Flxq_div_pre(sx, Flx_double(y, p), T, p, pi);
128 : }
129 1481075 : Q = cgetg(3,t_VEC);
130 1481075 : gel(Q,1) = Flx_sub(Flxq_sqr_pre(s, T, p, pi), Flx_double(x, p), p);
131 1481075 : if (typ(a4)==t_VEC) gel(Q, 1) = Flx_sub(gel(Q,1), gel(a4,1), p);
132 1481075 : gel(Q,2) = Flx_sub(Flxq_mul_pre(s, Flx_sub(x, gel(Q,1), p), T, p, pi),
133 : y, p);
134 1481075 : if (ps) *ps = s;
135 1481075 : 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 732609 : 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 732609 : if (ell_is_inf(P)) return Q;
151 727998 : if (ell_is_inf(Q)) return P;
152 727739 : Px = gel(P,1); Py = gel(P,2);
153 727739 : Qx = gel(Q,1); Qy = gel(Q,2);
154 727739 : if (Flx_equal(Px, Qx))
155 : {
156 75059 : if (Flx_equal(Py, Qy))
157 4271 : return FlxqE_dbl_slope(P, a4, T, p, pi, ps);
158 : else
159 70788 : return ellinf();
160 : }
161 652680 : s = Flxq_div_pre(Flx_sub(Py, Qy, p), Flx_sub(Px, Qx, p), T, p, pi);
162 652680 : R = cgetg(3,t_VEC);
163 652680 : gel(R,1) = Flx_sub(Flx_sub(Flxq_sqr_pre(s, T, p, pi), Px, p), Qx, p);
164 652680 : if (typ(a4)==t_VEC) gel(R,1) = Flx_sub(gel(R,1), gel(a4,1), p);
165 652680 : gel(R,2) = Flx_sub(Flxq_mul_pre(s, Flx_sub(Px, gel(R,1), p), T, p, pi),
166 : Py, p);
167 652680 : if (ps) *ps = s;
168 652680 : 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 6818 : FlxqE_neg_i(GEN P, ulong p)
181 : {
182 6818 : if (ell_is_inf(P)) return P;
183 6818 : 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 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 1576929 : _FlxqE_dbl(void *E, GEN P)
210 : {
211 1576929 : struct _FlxqE *e = (struct _FlxqE *) E;
212 1576929 : return FlxqE_dbl_slope(P, e->a4, e->T, e->p, e->pi, NULL);
213 : }
214 :
215 : static GEN
216 723030 : _FlxqE_add(void *E, GEN P, GEN Q)
217 : {
218 723030 : struct _FlxqE *e = (struct _FlxqE *) E;
219 723030 : return FlxqE_add_slope(P, Q, e->a4, e->T, e->p, e->pi, NULL);
220 : }
221 :
222 : static GEN
223 6818 : _FlxqE_sub(void *E, GEN P, GEN Q)
224 : {
225 6818 : struct _FlxqE *e = (struct _FlxqE *) E;
226 6818 : 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 309540 : _FlxqE_mul(void *E, GEN P, GEN n)
231 : {
232 309540 : pari_sp av = avma;
233 309540 : struct _FlxqE *e=(struct _FlxqE *) E;
234 309540 : long s = signe(n);
235 309540 : if (!s || ell_is_inf(P)) return ellinf();
236 308859 : if (s < 0) P = FlxqE_neg(P, e->T, e->p);
237 308859 : if (is_pm1(n)) return s>0? gcopy(P): P;
238 294014 : return gc_GEN(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 76930 : random_F3xqE(GEN a2, GEN a6, GEN T)
254 : {
255 76930 : pari_sp ltop = avma;
256 : GEN x, y, rhs;
257 76930 : const ulong p = 3;
258 : do
259 : {
260 155407 : set_avma(ltop);
261 155407 : x = random_Flx(get_Flx_degree(T),get_Flx_var(T),p);
262 155407 : rhs = Flx_add(Flxq_mul(Flxq_sqr(x, T, p), Flx_add(x, a2, p), T, p), a6, p);
263 155407 : } while ((!lgpol(rhs) && !lgpol(x)) || !Flxq_issquare(rhs, T, p));
264 76930 : y = Flxq_sqrt(rhs, T, p);
265 76930 : if (!y) pari_err_PRIME("random_F3xqE", T);
266 76930 : return gc_GEN(ltop, mkvec2(x, y));
267 : }
268 :
269 : /* Finds a random nonsingular point on E */
270 : static GEN
271 151164 : random_FlxqE_pre(GEN a4, GEN a6, GEN T, ulong p, ulong pi)
272 : {
273 151164 : pari_sp ltop = avma;
274 : GEN x, x2, y, rhs;
275 151164 : if (typ(a4)==t_VEC) return random_F3xqE(gel(a4,1), a6, T);
276 : do
277 : {
278 143631 : set_avma(ltop);
279 143631 : x = random_Flx(get_Flx_degree(T),get_Flx_var(T),p);
280 143631 : x2 = Flxq_sqr_pre(x, T, p, pi); /* x^3+a4*x+a6 = x*(x^2+a4)+a6 */
281 143631 : rhs = Flx_add(Flxq_mul_pre(x, Flx_add(x2, a4, p), T, p, pi), a6, p);
282 148895 : } while ((!lgpol(rhs) && !lgpol(Flx_add(Flx_triple(x2, p), a4, p)))
283 148832 : || !Flxq_issquare(rhs, T, p));
284 74234 : y = Flxq_sqrt(rhs, T, p);
285 74234 : if (!y) pari_err_PRIME("random_FlxqE", T);
286 74234 : return gc_GEN(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 74223 : _FlxqE_rand(void *E)
294 : {
295 74223 : struct _FlxqE *e=(struct _FlxqE *) E;
296 74223 : 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 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 68922 : FlxqE_vert(GEN P, GEN Q, GEN a4, GEN T, ulong p, ulong pi)
336 : {
337 68922 : long vT = get_Flx_var(T);
338 : GEN df;
339 68922 : if (ell_is_inf(P)) return pol1_Flx(vT);
340 47085 : if (!Flx_equal(gel(Q,1), gel(P,1))) return Flx_sub(gel(Q,1), gel(P,1), p);
341 1672 : if (lgpol(gel(P,2))!=0) return pol1_Flx(vT);
342 770 : df = typ(a4)==t_VEC ? Flxq_mul_pre(gel(P,1), Flx_double(gel(a4,1), p), T,p,pi)
343 1108 : : a4;
344 1108 : 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 29232 : FlxqE_Miller_line(GEN R, GEN Q, GEN slope, GEN a4, GEN T, ulong p, ulong pi)
350 : {
351 29232 : long vT = get_Flx_var(T);
352 29232 : GEN x = gel(Q,1), y = gel(Q,2);
353 29232 : GEN tmp1 = Flx_sub(x, gel(R,1), p);
354 29232 : GEN tmp2 = Flx_add(Flxq_mul_pre(tmp1, slope, T, p, pi), gel(R,2), p);
355 29232 : if (!Flx_equal(y, tmp2)) return Flx_sub(y, tmp2, p);
356 1153 : if (lgpol(y) == 0) return pol1_Flx(vT);
357 : else
358 : {
359 599 : GEN s1, s2, a2 = typ(a4)==t_VEC ? gel(a4,1): NULL;
360 599 : GEN y2i = Flxq_inv_pre(Flx_mulu(y, 2, p), T, p, pi);
361 599 : GEN df = a2 ? Flxq_mul_pre(x, Flx_mulu(a2, 2, p), T, p, pi): a4;
362 : GEN x3, ddf;
363 599 : s1 = Flxq_mul_pre(Flx_add(Flx_triple(Flxq_sqr_pre(x, T, p, pi), p), df, p), y2i, T, p, pi);
364 599 : if (!Flx_equal(s1, slope)) return Flx_sub(s1, slope, p);
365 289 : x3 = Flx_triple(x, p);
366 289 : ddf = a2 ? Flx_add(x3, a2, p): x3;
367 289 : s2 = Flxq_mul_pre(Flx_sub(ddf, Flxq_sqr_pre(s1, T,p,pi), p), y2i, T,p,pi);
368 289 : 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 46944 : FlxqE_tangent_update(GEN R, GEN Q, GEN a4, GEN T, ulong p, ulong pi, GEN *pt_R)
376 : {
377 46944 : if (ell_is_inf(R))
378 : {
379 4080 : *pt_R = ellinf();
380 4080 : return pol1_Flx(get_Flx_var(T));
381 : }
382 42864 : else if (!lgpol(gel(R,2)))
383 : {
384 16393 : *pt_R = ellinf();
385 16393 : return FlxqE_vert(R, Q, a4, T, p, pi);
386 : } else {
387 : GEN slope;
388 26471 : *pt_R = FlxqE_dbl_slope(R, a4, T, p, pi, &slope);
389 26471 : 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 4173 : FlxqE_chord_update(GEN R, GEN P, GEN Q, GEN a4, GEN T, ulong p, ulong pi, GEN *pt_R)
397 : {
398 4173 : if (ell_is_inf(R))
399 : {
400 48 : *pt_R = gcopy(P);
401 48 : return FlxqE_vert(P, Q, a4, T, p, pi);
402 : }
403 4125 : 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 4125 : else if (Flx_equal(gel(P, 1), gel(R, 1)))
409 : {
410 1364 : 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 1364 : *pt_R = ellinf();
415 1364 : return FlxqE_vert(R, Q, a4, T, p, pi);
416 : }
417 : } else {
418 : GEN slope;
419 2761 : *pt_R = FlxqE_add_slope(P, R, a4, T, p, pi, &slope);
420 2761 : 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 46944 : FlxqE_Miller_dbl(void* E, GEN d)
432 : {
433 46944 : struct _FlxqE_miller *m = (struct _FlxqE_miller *)E;
434 46944 : ulong p = m->p, pi = m->pi;
435 46944 : GEN T = m->T, a4 = m->a4, P = m->P;
436 46944 : GEN v, line, point = gel(d,3);
437 46944 : GEN N = Flxq_sqr_pre(gel(d,1), T, p, pi);
438 46944 : GEN D = Flxq_sqr_pre(gel(d,2), T, p, pi);
439 46944 : line = FlxqE_tangent_update(point, P, a4, T, p, pi, &point);
440 46944 : N = Flxq_mul_pre(N, line, T, p, pi);
441 46944 : v = FlxqE_vert(point, P, a4, T, p, pi);
442 46944 : D = Flxq_mul_pre(D, v, T, p, pi); return mkvec3(N, D, point);
443 : }
444 :
445 : static GEN
446 4173 : FlxqE_Miller_add(void* E, GEN va, GEN vb)
447 : {
448 4173 : struct _FlxqE_miller *m = (struct _FlxqE_miller *)E;
449 4173 : ulong p = m->p, pi = m->pi;
450 4173 : GEN T = m->T, a4 = m->a4, P = m->P;
451 : GEN v, line, point;
452 4173 : GEN na = gel(va,1), da = gel(va,2), pa = gel(va,3);
453 4173 : GEN nb = gel(vb,1), db = gel(vb,2), pb = gel(vb,3);
454 4173 : GEN N = Flxq_mul_pre(na, nb, T, p, pi);
455 4173 : GEN D = Flxq_mul_pre(da, db, T, p, pi);
456 4173 : line = FlxqE_chord_update(pa, pb, P, a4, T, p, pi, &point);
457 4173 : N = Flxq_mul_pre(N, line, T, p, pi);
458 4173 : v = FlxqE_vert(point, P, a4, T, p, pi);
459 4173 : 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 17709 : FlxqE_Miller(GEN Q, GEN P, GEN m, GEN a4, GEN T, ulong p, ulong pi)
466 : {
467 17709 : pari_sp av = avma;
468 : struct _FlxqE_miller d;
469 : GEN v, N, D, g1;
470 :
471 17709 : d.a4 = a4; d.T = T; d.p = p; d.P = P; d.pi = pi;
472 17709 : g1 = pol1_Flx(get_Flx_var(T));
473 17709 : v = gen_pow_i(mkvec3(g1,g1,Q), m, (void*)&d,
474 : FlxqE_Miller_dbl, FlxqE_Miller_add);
475 17709 : N = gel(v,1); D = gel(v,2);
476 17709 : return gc_upto(av, Flxq_div_pre(N, D, T, p, pi));
477 : }
478 :
479 : GEN
480 14568 : FlxqE_weilpairing_pre(GEN P, GEN Q, GEN m, GEN a4, GEN T, ulong p, ulong pi)
481 : {
482 14568 : pari_sp av = avma;
483 : GEN N, D, w;
484 14568 : if (ell_is_inf(P) || ell_is_inf(Q)
485 11563 : || (Flx_equal(gel(P,1),gel(Q,1)) && Flx_equal(gel(P,2),gel(Q,2))))
486 5745 : return pol1_Flx(get_Flx_var(T));
487 8823 : N = FlxqE_Miller(P, Q, m, a4, T, p, pi);
488 8823 : D = FlxqE_Miller(Q, P, m, a4, T, p, pi);
489 8823 : w = Flxq_div_pre(N, D, T, p, pi); if (mpodd(m)) w = Flx_neg(w, p);
490 8823 : 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 14547 : _FlxqE_pairorder(void *E, GEN P, GEN Q, GEN m, GEN F)
505 : {
506 14547 : struct _FlxqE *e = (struct _FlxqE *) E;
507 14547 : 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 : /** 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 gc_upto(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 gc_upto(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 : s = Qp_exp(cvtop(s, p, e-1)); return modii(padic_u(s), 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 gc_upto(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 : /* These polynomials were extracted from the ECHIDNA databases
843 : * available at <http://echidna.maths.usyd.edu.au/echidna/>
844 : * and computed by David R. Kohel.
845 : * Return the matrix of the modular polynomial, set act to the parametrization,
846 : * and set dj to the opposite of the supersingular j-invariant. */
847 : static GEN
848 25032 : get_Kohel_polynomials(ulong p, GEN *act, long *dj)
849 : {
850 25032 : const long mat3[] = {-1,-36,-270};
851 25032 : const long num3[] = {1,-483,-21141,-59049};
852 25032 : const long den3[] = {1,261, 4347, -6561};
853 25032 : const long mat5[] = {-1,-30,-315,-1300,-1575};
854 25032 : const long num5[] = {-1,490,20620,158750,78125};
855 25032 : const long den5[] = {-1,-254,-4124,-12250,3125};
856 25032 : const long mat7[] = {-1,-28,-322,-1904,-5915,-8624,-4018};
857 25032 : const long num7[] = {1,-485,-24058,-343833,-2021642,-4353013,-823543};
858 25032 : const long den7[] = {1,259,5894,49119,168406,166355,-16807};
859 25032 : const long mat13[]= {-1,-26,-325,-2548,-13832,-54340,-157118,-333580,-509366,
860 : -534820,-354536,-124852,-15145};
861 25032 : const long num13[]= {1,-487,-24056,-391463,-3396483,-18047328,-61622301,
862 : -133245853,-168395656,-95422301,-4826809};
863 25032 : const long den13[]= {1,257,5896,60649,364629,1388256,3396483,5089019,4065464,
864 : 1069939,-28561};
865 25032 : switch(p)
866 : {
867 14805 : case 3:
868 14805 : *dj = 0;
869 14805 : return fill_pols(3,mat3,4,num3,den3,act);
870 10178 : case 5:
871 10178 : *dj = 0;
872 10178 : return fill_pols(5,mat5,5,num5,den5,act);
873 35 : case 7:
874 35 : *dj = 1;
875 35 : return fill_pols(7,mat7,7,num7,den7,act);
876 14 : case 13:
877 14 : *dj = 8;
878 14 : return fill_pols(13,mat13,11,num13,den13,act);
879 : }
880 : *dj=0; *act = NULL; return NULL; /* LCOV_EXCL_LINE */
881 : }
882 :
883 : long
884 25053 : zx_is_pcyc(GEN T)
885 : {
886 25053 : long i, n = degpol(T);
887 25053 : if (!uisprime(n+1)) return 0;
888 95249 : for (i = 0; i <= n; i++)
889 84147 : if (T[i+2] != 1UL) return 0;
890 11102 : return 1;
891 : }
892 :
893 : static GEN
894 25032 : Flxq_ellcard_Kohel(GEN a4, GEN a6, GEN T, ulong p)
895 : {
896 25032 : pari_sp av = avma, av2;
897 : pari_timer ti;
898 25032 : long n = get_Flx_degree(T), N = (n+4)/2, dj;
899 25032 : GEN q = powuu(p, N);
900 : GEN T2, Xm, s1, c2, t, lr, S1, sqx, Nc2, Np;
901 25032 : GEN act, phi = get_Kohel_polynomials(p, &act, &dj);
902 25032 : long ispcyc = zx_is_pcyc(get_Flx_mod(T));
903 25032 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
904 25032 : timer_start(&ti);
905 25032 : if (!ispcyc)
906 : {
907 13937 : T2 = Flx_Teichmuller(get_Flx_mod(T),p,N);
908 13937 : if (DEBUGLEVEL) timer_printf(&ti,"Teich");
909 : } else
910 11095 : T2 = Flx_to_ZX(get_Flx_mod(T));
911 :
912 25032 : T2 = FpX_get_red(T2, q); T = ZXT_to_FlxT(T2, p);
913 25032 : av2 = avma;
914 25032 : if (DEBUGLEVEL) timer_printf(&ti,"Barrett");
915 25032 : if (!ispcyc)
916 : {
917 13937 : Xm = FpXQ_powers(pol_xn(n,get_FpX_var(T2)),p-1,T2,q);
918 13937 : if (DEBUGLEVEL) timer_printf(&ti,"Xm");
919 : } else
920 11095 : Xm = cgetg(1,t_VEC);
921 25032 : s1 = Flxq_inv_pre(Flx_Fl_add(Flxq_ellj(a4,a6,T,p),dj, p),T,p,pi);
922 25032 : lr = Flxq_lroot_pre(polx_Flx(get_Flx_var(T)), T,p,pi);
923 25032 : sqx = Flxq_powers_pre(lr, p-1, T, p, pi);
924 25032 : S1 = lift_isogeny(phi, Flx_to_ZX(s1), N, Xm, T2, sqx, T,p,pi);
925 25032 : if (DEBUGLEVEL) timer_printf(&ti,"Lift isogeny");
926 25032 : c2 = getc2(act, S1, T2, q, p, N);
927 25032 : if (DEBUGLEVEL) timer_printf(&ti,"c^2");
928 25032 : if (p>3 && !ispcyc)
929 : {
930 10220 : GEN c2p = Flx_to_ZX(Flxq_inv_pre(ZX_to_Flx(c2,p),T,p,pi));
931 10220 : GEN tc2 = Teichmuller_lift(c2p,Xm, T2,sqx,T,p,pi,N);
932 10220 : if (DEBUGLEVEL) timer_printf(&ti,"Teichmuller/Fq");
933 10220 : c2 = FpX_rem(FpX_mul(tc2,c2,q),T2,q);
934 : }
935 25032 : c2 = gc_upto(av2, c2);
936 25032 : if (DEBUGLEVEL) timer_printf(&ti,"tc2");
937 25032 : Nc2 = (ispcyc? ZpXQ_sqrtnorm_pcyc: ZpXQ_sqrtnorm)(c2, T2, q, utoipos(p), N);
938 25032 : if (DEBUGLEVEL) timer_printf(&ti,"Norm");
939 25032 : Np = get_norm(a4,a6,T,p,pi,N);
940 25032 : if (p>3 && ispcyc)
941 : {
942 7 : GEN Ncpi = utoi(Fl_inv(umodiu(Nc2,p), p));
943 7 : GEN tNc2 = Zp_sqrtnlift(gen_1, subss(p,1), Ncpi, utoipos(p),N);
944 7 : if (DEBUGLEVEL) timer_printf(&ti,"Teichmuller/Fp");
945 7 : Nc2 = Fp_mul(Nc2,tNc2,q);
946 : }
947 25032 : t = Fp_center_i(Fp_mul(Nc2,Np,q),q,shifti(q,-1));
948 25032 : return gc_upto(av, subii(addiu(powuu(p,n),1),t));
949 : }
950 :
951 : /* Use Damien Robert's method */
952 : static GEN
953 21 : get_trace_Robert(GEN J, GEN phi, GEN Xm, GEN T, GEN q, ulong p, long e)
954 : {
955 21 : long n = lg(phi)-2;
956 21 : GEN K = ZpXQ_frob(J, Xm, T, q, p);
957 21 : GEN Jp = FpXQ_powers(J, n, T, q);
958 21 : GEN Kp = FpXQ_powers(K, n, T, q);
959 21 : GEN Jd = FpXC_powderiv(Jp, q);
960 21 : GEN Kd = FpXC_powderiv(Kp, q);
961 21 : GEN Dx = FpM_FpXQV_bilinear(phi, Kd, Jp, T, q);
962 21 : GEN Dy = FpM_FpXQV_bilinear(phi, Kp, Jd, T, q);
963 21 : GEN C = ZpXQ_inv(ZX_divuexact(Dy, p), T, utoi(p), e);
964 21 : return FpX_neg(FpXQ_mul(Dx, C, T, q), q);
965 : }
966 :
967 : /* in p^2, so p is tiny */
968 : static GEN
969 21 : Flxq_ellcard_Harley(GEN a4, GEN a6, GEN T, ulong p)
970 : {
971 21 : pari_sp av = avma, av2;
972 : pari_timer ti;
973 21 : long n = get_Flx_degree(T), N = (n+5)/2;
974 21 : GEN pp = utoipos(p), q = powuu(p, N);
975 : GEN T2, j, t, phi, J1, sqx, Xm, c2, tc2, c2p, Nc2, Np;
976 21 : long ispcyc = zx_is_pcyc(get_Flx_mod(T));
977 21 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p); /* = 0 here */
978 21 : timer_start(&ti);
979 21 : if (!ispcyc)
980 : {
981 14 : T2 = Flx_Teichmuller(get_Flx_mod(T),p,N);
982 14 : if (DEBUGLEVEL) timer_printf(&ti,"Teich");
983 : } else
984 7 : T2 = Flx_to_ZX(get_Flx_mod(T));
985 21 : T2 = FpX_get_red(T2, q); T = ZXT_to_FlxT(T2, p);
986 21 : av2 = avma;
987 21 : if (DEBUGLEVEL) timer_printf(&ti,"Barrett");
988 21 : if (!ispcyc)
989 : {
990 14 : Xm = FpXQ_powers(pol_xn(n,get_FpX_var(T2)),p-1,T2,q);
991 14 : if (DEBUGLEVEL) timer_printf(&ti,"Xm");
992 : } else
993 7 : Xm = cgetg(1,t_VEC);
994 21 : j = Flxq_ellj(a4,a6,T,p);
995 21 : sqx = Flxq_powers_pre(Flxq_lroot_pre(polx_Flx(T[1]), T,p,pi), p-1, T,p,pi);
996 21 : phi = polmodular_ZM(p, 0);
997 21 : if (DEBUGLEVEL) timer_printf(&ti,"phi");
998 21 : J1 = lift_isogeny(phi, Flx_to_ZX(j), N, Xm, T2,sqx,T,p,pi);
999 21 : if (DEBUGLEVEL) timer_printf(&ti,"Lift isogeny");
1000 21 : c2 = get_trace_Robert(J1, phi, Xm, T2, q, p, N);
1001 21 : q = diviuexact(q,p); N--;
1002 21 : if (DEBUGLEVEL) timer_printf(&ti,"c^2");
1003 21 : if (!ispcyc)
1004 : {
1005 14 : c2p = Flx_to_ZX(Flxq_inv_pre(ZX_to_Flx(c2,p),T,p,pi));
1006 14 : tc2 = Teichmuller_lift(c2p,Xm, T2,sqx,T,p,pi,N);
1007 14 : if (DEBUGLEVEL) timer_printf(&ti,"teichmuller");
1008 14 : c2 = FpX_rem(FpX_mul(tc2,c2,q),T2,q);
1009 : }
1010 21 : c2 = gc_upto(av2, c2);
1011 21 : q = powuu(p, N);
1012 21 : Nc2 = (ispcyc? ZpXQ_sqrtnorm_pcyc: ZpXQ_sqrtnorm)(c2, T2, q, pp, N);
1013 21 : if (DEBUGLEVEL) timer_printf(&ti,"Norm");
1014 21 : Np = get_norm(a4,a6,T,p,pi,N);
1015 21 : if (ispcyc)
1016 : {
1017 7 : GEN Ncpi = utoi(Fl_inv(umodiu(Nc2,p), p));
1018 7 : GEN tNc2 = Zp_sqrtnlift(gen_1, subss(p,1), Ncpi, pp, N);
1019 7 : if (DEBUGLEVEL) timer_printf(&ti,"Teichmuller/Fp");
1020 7 : Nc2 = Fp_mul(Nc2,tNc2,q);
1021 : }
1022 21 : t = Fp_center_i(Fp_mul(Nc2,Np,q),q,shifti(q,-1));
1023 21 : return gc_upto(av, subii(addiu(powuu(p,n),1),t));
1024 : }
1025 :
1026 : /***************************************************************************/
1027 : /* Shanks-Mestre */
1028 : /***************************************************************************/
1029 :
1030 : /* Return the lift of a (mod b), which is closest to h */
1031 : static GEN
1032 10352 : closest_lift(GEN a, GEN b, GEN h)
1033 10352 : { return addii(a, mulii(b, diviiround(subii(h,a), b))); }
1034 :
1035 : /* find multiple of order of f using Baby Step/Giant Step, f^h close to 1,
1036 : * order lies in an interval of size <= 'bound' and known mod B */
1037 : static GEN
1038 5578 : _FlxqE_order_multiple(void *E, GEN f, GEN h, GEN bound, GEN B)
1039 : {
1040 5578 : pari_sp av = avma, av1;
1041 : pari_timer Ti;
1042 5578 : long i, s = ceilsqrtdiv(bound, B) >> 1;
1043 : GEN P, F, tx, ti, fg, fh;
1044 :
1045 5578 : P = fh = _FlxqE_mul(E, f, h);
1046 5578 : if (DEBUGLEVEL >= 6) timer_start(&Ti);
1047 5578 : if (ell_is_inf(fh)) return h;
1048 5167 : F = _FlxqE_mul(E, f, B);
1049 5167 : if (s < 3)
1050 : { /* we're nearly done: naive search */
1051 1116 : GEN Q = P;
1052 1116 : for (i=1;; i++)
1053 : {
1054 3661 : P = _FlxqE_add(E, P, F); /* h.f + i.F */
1055 3661 : if (ell_is_inf(P)) return gc_upto(av, addii(h, mului(i,B)));
1056 3099 : Q = _FlxqE_sub(E, Q, F); /* h.f - i.F */
1057 3099 : if (ell_is_inf(Q)) return gc_upto(av, subii(h, mului(i,B)));
1058 : }
1059 : }
1060 4051 : tx = cgetg(s+1,t_VECSMALL); av1 = avma;
1061 43291 : for (i=1; i<=s; i++)
1062 : { /* baby steps */
1063 39572 : tx[i] = hash_GEN(gel(P, 1));
1064 39572 : P = _FlxqE_add(E, P, F); /* h.f + i.F */
1065 39572 : if (ell_is_inf(P)) return gc_upto(av, addii(h, mului(i,B)));
1066 39240 : if (gc_needed(av1,3))
1067 : {
1068 0 : if(DEBUGMEM>1) pari_warn(warnmem,"[Flxq_ellcard] baby steps, i=%ld",i);
1069 0 : P = gc_upto(av1,P);
1070 : }
1071 : }
1072 3719 : if (DEBUGLEVEL >= 6) timer_printf(&Ti,"[Flxq_ellcard] baby steps, s = %ld",s);
1073 : /* giant steps: fg = s.F */
1074 3719 : fg = gc_upto(av1, _FlxqE_sub(E, P, fh));
1075 3719 : if (ell_is_inf(fg)) return gc_upto(av, mului(s,B));
1076 3719 : ti = vecsmall_indexsort(tx); /* = permutation sorting tx */
1077 3719 : tx = perm_mul(tx,ti);
1078 3719 : if (DEBUGLEVEL >= 6) timer_printf(&Ti, "[Flxq_ellcard] sorting");
1079 3719 : av1 = avma;
1080 3719 : for (P=fg, i=1; ; i++)
1081 35493 : {
1082 39212 : long k = hash_GEN(gel(P,1)), r = zv_search(tx, k);
1083 39212 : if (r)
1084 : {
1085 7443 : while (r && tx[r] == k) r--;
1086 3719 : for (r++; r <= s && tx[r] == k; r++)
1087 : {
1088 3719 : long j = ti[r]-1;
1089 3719 : GEN Q = _FlxqE_add(E, _FlxqE_mul(E, F, stoi(j)), fh);
1090 3719 : if (DEBUGLEVEL >= 6)
1091 0 : timer_printf(&Ti, "[Flxq_ellcard] giant steps, i = %ld",i);
1092 3719 : if (Flx_equal(gel(P,1), gel(Q,1)))
1093 : {
1094 3719 : if (Flx_equal(gel(P,2), gel(Q,2))) i = -i;
1095 3719 : return gc_upto(av,addii(h, mulii(addis(mulss(s,i), j), B)));
1096 : }
1097 : }
1098 : }
1099 35493 : P = _FlxqE_add(E, P, fg);
1100 35493 : if (gc_needed(av1,3))
1101 : {
1102 0 : if(DEBUGMEM>1) pari_warn(warnmem,"[Flxq_ellcard] giants steps, i=%ld",i);
1103 0 : P = gc_upto(av1,P);
1104 : }
1105 : }
1106 : }
1107 : static GEN
1108 5578 : _FlxqE_order(void *E, GEN f, GEN h, GEN bound, GEN B)
1109 : {
1110 5578 : GEN o = _FlxqE_order_multiple(E, f, h, bound, B);
1111 5578 : return gen_order(f, o, E, &FlxqE_group);
1112 : }
1113 :
1114 : static void
1115 64365 : Flx_next(GEN t, ulong p)
1116 : {
1117 : long i;
1118 78358 : for(i=2;;i++)
1119 78358 : if (uel(t,i)==p-1) t[i]=0; else { t[i]++; break; }
1120 64365 : }
1121 :
1122 : static void
1123 64365 : Flx_renormalize_ip(GEN x, long lx)
1124 : {
1125 : long i;
1126 78358 : for (i = lx-1; i>=2; i--)
1127 71267 : if (x[i]) break;
1128 64365 : setlg(x, i+1);
1129 64365 : }
1130 :
1131 : static ulong
1132 6188 : F3xq_ellcard_naive(GEN a2, GEN a6, GEN T)
1133 : {
1134 6188 : pari_sp av = avma;
1135 6188 : long i, d = get_Flx_degree(T), lx = d+2;
1136 6188 : long q = upowuu(3, d), a;
1137 6188 : GEN x = zero_zv(lx); x[1] = get_Flx_var(T);
1138 26978 : for(a=1, i=0; i<q; i++)
1139 : {
1140 : GEN rhs;
1141 20790 : Flx_renormalize_ip(x, lx);
1142 20790 : rhs = Flx_add(Flxq_mul(Flxq_sqr(x, T, 3), Flx_add(x, a2, 3), T, 3), a6, 3);
1143 20790 : if (!lgpol(rhs)) a++; else if (Flxq_issquare(rhs, T, 3)) a+=2;
1144 20790 : Flx_next(x, 3);
1145 : }
1146 6188 : set_avma(av); return a;
1147 : }
1148 :
1149 : /* p^deg(T) is tiny */
1150 : static ulong
1151 903 : Flxq_ellcard_naive(GEN a4, GEN a6, GEN T, ulong p)
1152 : {
1153 903 : pari_sp av = avma;
1154 903 : long i, d = get_Flx_degree(T), lx = d+2;
1155 903 : long q = upowuu(p, d), a;
1156 903 : GEN x = zero_zv(lx); x[1] = get_Flx_var(T);
1157 44478 : for(a = 1, i = 0; i < q; i++)
1158 : {
1159 : GEN x2, rhs;
1160 43575 : Flx_renormalize_ip(x, lx);
1161 43575 : x2 = Flxq_sqr_pre(x, T, p, 0);
1162 43575 : rhs = Flx_add(Flxq_mul_pre(x, Flx_add(x2, a4, p), T, p, 0), a6, p);
1163 43575 : if (!lgpol(rhs)) a++; else if (Flxq_issquare(rhs,T,p)) a += 2;
1164 43575 : Flx_next(x,p);
1165 : }
1166 903 : set_avma(av); return a;
1167 : }
1168 :
1169 : static long
1170 11373 : Flxq_kronecker(GEN x, GEN T, ulong p)
1171 : {
1172 : pari_sp av;
1173 11373 : if (lgpol(x) == 0) return 0;
1174 11351 : av = avma; return gc_long(av, krouu(Flxq_norm(x, T, p), p));
1175 : }
1176 :
1177 : /* Find x such that kronecker(u = x^3+a4x+a6, p) is KRO.
1178 : * Return point [x*u,u^2] on E (KRO=1) / E^twist (KRO=-1) */
1179 : static GEN
1180 5578 : Flxq_ellpoint(long KRO, GEN a4, GEN a6, GEN T, ulong p, ulong pi)
1181 : {
1182 5578 : long v = get_Flx_var(T), n = get_Flx_degree(T);
1183 : for(;;)
1184 5795 : {
1185 11373 : GEN x = random_Flx(n, v, p), x2 = Flxq_sqr_pre(x,T,p,pi);
1186 11373 : GEN u = Flx_add(a6, Flxq_mul_pre(Flx_add(a4, x2, p), x, T,p, pi), p);
1187 11373 : if (Flxq_kronecker(u,T,p) == KRO)
1188 5578 : return mkvec2(Flxq_mul_pre(u,x, T,p,pi), Flxq_sqr_pre(u, T,p,pi));
1189 : }
1190 : }
1191 :
1192 : static GEN
1193 4774 : Flxq_ellcard_Shanks(GEN a4, GEN a6, GEN q, GEN T, ulong p)
1194 : {
1195 4774 : pari_sp av = avma;
1196 4774 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
1197 4774 : long v = get_Flx_var(T), KRO = -1;
1198 : GEN h,f, A, B;
1199 4774 : GEN q1p = addiu(q,1), q2p = shifti(q1p, 1);
1200 4774 : GEN bound = addiu(sqrti(gmul2n(q,4)), 1); /* ceil( 4sqrt(q) ) */
1201 : struct _FlxqE e;
1202 4774 : e.p = p; e.pi = pi; e.T = Flx_get_red_pre(T, p, pi);
1203 : /* once #E(Flxq) is known mod B >= bound, it is determined */
1204 4774 : switch(FlxqX_nbroots(mkpoln(4, pol1_Flx(v), pol0_Flx(v), a4, a6), T, p))
1205 : { /* how many 2-torsion points ? */
1206 2891 : case 3: A = gen_0; B = utoipos(4); break;
1207 1414 : case 1: A = gen_0; B = gen_2; break;
1208 469 : default: A = gen_1; B = gen_2; break; /* 0 */
1209 : }
1210 : for(;;)
1211 : {
1212 5578 : h = closest_lift(A, B, q1p);
1213 : /* [ux, u^2] is on E_u: y^2 = x^3 + c4 u^2 x + c6 u^3
1214 : * E_u isomorphic to E (resp. E') iff KRO = 1 (resp. -1)
1215 : * #E(F_p) = p+1 - a_p, #E'(F_p) = p+1 + a_p
1216 : *
1217 : * #E_u(Flxq) = A (mod B), h is close to #E_u(Flxq) */
1218 5578 : KRO = -KRO;
1219 5578 : f = Flxq_ellpoint(KRO, a4,a6, T,p,pi);
1220 5578 : e.a4 = Flxq_mul_pre(a4, gel(f,2), T,p,pi); /* a4 for E_u */
1221 5578 : h = _FlxqE_order((void*)&e, f, h, bound, B);
1222 : /* h | #E_u(Flxq) = A (mod B) */
1223 5578 : A = Z_chinese_all(A, gen_0, B, h, &B);
1224 5578 : if (cmpii(B, bound) >= 0) break;
1225 : /* not done, update A mod B for the _next_ curve, isomorphic to
1226 : * the quadratic twist of this one */
1227 804 : A = remii(subii(q2p,A), B); /* #E(Fq)+#E'(Fq) = 2q+2 */
1228 : }
1229 4774 : h = closest_lift(A, B, q1p);
1230 4774 : return gc_INT(av, KRO == 1? h: subii(q2p,h));
1231 : }
1232 :
1233 : static GEN
1234 20993 : F3xq_ellcard(GEN a2, GEN a6, GEN T)
1235 : {
1236 20993 : long n = get_Flx_degree(T);
1237 20993 : if (n <= 2)
1238 5887 : return utoi(F3xq_ellcard_naive(a2, a6, T));
1239 : else
1240 : {
1241 15106 : GEN q1 = addiu(powuu(3, get_Flx_degree(T)), 1), t;
1242 15106 : GEN a = Flxq_div(a6,Flxq_powu(a2,3,T,3),T,3);
1243 15106 : if (Flx_equal1(Flxq_powu(a, 8, T, 3)))
1244 : {
1245 301 : GEN P = Flxq_minpoly(a,T,3);
1246 301 : long dP = degpol(P); /* dP <= 2 */
1247 301 : ulong q = upowuu(3,dP);
1248 301 : GEN A2 = pol1_Flx(P[1]), A6 = Flx_rem(polx_Flx(P[1]), P, 3);
1249 301 : long tP = q + 1 - F3xq_ellcard_naive(A2, A6, P);
1250 301 : t = elltrace_extension(stoi(tP), n/dP, utoi(q));
1251 301 : if (umodiu(t, 3)!=1) t = negi(t);
1252 301 : return Flx_equal1(a2) || Flxq_issquare(a2,T,3) ? subii(q1,t): addii(q1,t);
1253 : }
1254 14805 : else return Flxq_ellcard_Kohel(mkvec(a2), a6, T, 3);
1255 : }
1256 : }
1257 :
1258 : static GEN
1259 11151 : Flxq_ellcard_Satoh(GEN a4, GEN a6, GEN j, GEN T, ulong p)
1260 : {
1261 11151 : long n = get_Flx_degree(T);
1262 11151 : if (n <= 2)
1263 903 : return utoi(Flxq_ellcard_naive(a4, a6, T, p));
1264 : else
1265 : {
1266 10248 : GEN jp = Flxq_powu(j, p, T, p);
1267 10248 : GEN s = Flx_add(j, jp, p);
1268 10248 : if (degpol(s) <= 0)
1269 : { /* it is assumed j not in F_p */
1270 0 : GEN m = Flxq_mul(j, jp, T, p);
1271 0 : if (degpol(m) <= 0)
1272 : {
1273 0 : GEN q = sqru(p);
1274 0 : GEN q1 = addiu(powuu(p, get_Flx_degree(T)), 1);
1275 0 : GEN sk = Flx_Fl_add(Flx_neg(j, p), 1728%p, p);
1276 0 : GEN sA4 = Flx_triple(Flxq_mul(sk, j, T, p), p);
1277 0 : GEN u = Flxq_div(a4, sA4, T, p);
1278 0 : ulong ns = lgpol(s) ? Fl_neg(s[2], p): 0UL;
1279 0 : GEN P = mkvecsmall4(T[1], m[2], ns, 1L);
1280 : GEN A4, A6, t, tP;
1281 0 : Flxq_ellj_to_a4a6(polx_Flx(T[1]), P, p, &A4, &A6);
1282 0 : tP = addis(q, 1 - Flxq_ellcard_naive(A4, A6, P, p));
1283 0 : t = elltrace_extension(tP, n>>1, q);
1284 0 : return Flxq_is2npower(u, 2, T, p) ? subii(q1,t): addii(q1,t);
1285 : }
1286 : }
1287 10248 : if (p<=7 || p==13) return Flxq_ellcard_Kohel(a4, a6, T, p);
1288 21 : else return Flxq_ellcard_Harley(a4, a6, T, p);
1289 : }
1290 : }
1291 :
1292 : static GEN
1293 0 : Flxq_ellcard_Kedlaya(GEN a4, GEN a6, GEN T, ulong p)
1294 : {
1295 0 : pari_sp av = avma;
1296 0 : GEN H = mkpoln(4, gen_1, gen_0, Flx_to_ZX(a4), Flx_to_ZX(a6));
1297 0 : GEN Tp = Flx_to_ZX(get_Flx_mod(T));
1298 0 : long n = degpol(Tp), e = ((p < 16 ? n+1: n)>>1)+1;
1299 0 : GEN M = ZlXQX_hyperellpadicfrobenius(H, Tp, p, e);
1300 0 : GEN N = ZpXQM_prodFrobenius(M, Tp, utoipos(p), e);
1301 0 : GEN q = powuu(p, e);
1302 0 : GEN tp = Fq_add(gcoeff(N,1,1), gcoeff(N,2,2), Tp, q);
1303 0 : GEN t = Fp_center_i(typ(tp)==t_INT ? tp: leading_coeff(tp), q, shifti(q,-1));
1304 0 : return gc_upto(av, subii(addiu(powuu(p, n), 1), t));
1305 : }
1306 :
1307 : GEN
1308 57313 : Flxq_ellj(GEN a4, GEN a6, GEN T, ulong p)
1309 : {
1310 57313 : pari_sp av=avma;
1311 57313 : if (p==3)
1312 : {
1313 : GEN J;
1314 14805 : if (typ(a4)!=t_VEC) return pol0_Flx(get_Flx_var(T));
1315 14805 : J = Flxq_div(Flxq_powu(gel(a4,1),3, T, p),Flx_neg(a6,p), T, p);
1316 14805 : return gc_leaf(av, J);
1317 : }
1318 : else
1319 : {
1320 42508 : pari_sp av=avma;
1321 42508 : GEN a43 = Flxq_mul(a4,Flxq_sqr(a4,T,p),T,p);
1322 42508 : GEN a62 = Flxq_sqr(a6,T,p);
1323 42508 : GEN num = Flx_mulu(a43,6912,p);
1324 42508 : GEN den = Flx_add(Flx_mulu(a43,4,p),Flx_mulu(a62,27,p),p);
1325 42508 : return gc_leaf(av, Flxq_div(num, den, T, p));
1326 : }
1327 : }
1328 :
1329 : void
1330 0 : Flxq_ellj_to_a4a6(GEN j, GEN T, ulong p, GEN *pt_a4, GEN *pt_a6)
1331 : {
1332 0 : ulong zagier = 1728 % p;
1333 0 : if (lgpol(j)==0)
1334 0 : { *pt_a4 = pol0_Flx(T[1]); *pt_a6 =pol1_Flx(T[1]); }
1335 0 : else if (lgpol(j)==1 && uel(j,2) == zagier)
1336 0 : { *pt_a4 = pol1_Flx(T[1]); *pt_a6 =pol0_Flx(T[1]); }
1337 : else
1338 : {
1339 0 : GEN k = Flx_Fl_add(Flx_neg(j, p), zagier, p);
1340 0 : GEN kj = Flxq_mul(k, j, T, p);
1341 0 : GEN k2j = Flxq_mul(kj, k, T, p);
1342 0 : *pt_a4 = Flx_triple(kj, p);
1343 0 : *pt_a6 = Flx_double(k2j, p);
1344 : }
1345 0 : }
1346 :
1347 : static GEN
1348 9212 : F3xq_ellcardj(GEN a4, GEN a6, GEN T, GEN q, long n)
1349 : {
1350 9212 : const ulong p = 3;
1351 : ulong t;
1352 9212 : GEN q1 = addiu(q,1);
1353 9212 : GEN na4 = Flx_neg(a4,p), ra4;
1354 9212 : if (!Flxq_issquare(na4,T,p))
1355 4536 : return q1;
1356 4676 : ra4 = Flxq_sqrt(na4,T,p);
1357 4676 : t = Flxq_trace(Flxq_div(a6,Flxq_mul(na4,ra4,T,p),T,p),T,p);
1358 4676 : if (n%2==1)
1359 : {
1360 : GEN q3;
1361 2499 : if (t==0) return q1;
1362 812 : q3 = powuu(p,(n+1)>>1);
1363 812 : return (t==1)^(n%4==1) ? subii(q1,q3): addii(q1,q3);
1364 : }
1365 : else
1366 : {
1367 2177 : GEN q22, q2 = powuu(p,n>>1);
1368 2177 : GEN W = Flxq_pow(a4,shifti(q,-2),T,p);
1369 2177 : long s = (W[2]==1)^(n%4==2);
1370 2177 : if (t!=0) return s ? addii(q1,q2): subii(q1, q2);
1371 2177 : q22 = shifti(q2,1);
1372 2177 : return s ? subii(q1,q22): addii(q1, q22);
1373 : }
1374 : }
1375 :
1376 : static GEN
1377 15827 : Flxq_ellcardj(GEN a4, GEN a6, ulong j, GEN T, GEN q, ulong p, long n)
1378 : {
1379 15827 : GEN q1 = addiu(q,1);
1380 15827 : if (j==0)
1381 : {
1382 : ulong w;
1383 : GEN W, t, N;
1384 5677 : if (umodiu(q,6)!=1) return q1;
1385 4277 : N = Fp_ffellcard(gen_0,gen_1,q,n,utoipos(p));
1386 4277 : t = subii(q1, N);
1387 4277 : W = Flxq_pow(a6,diviuexact(shifti(q,-1), 3),T,p);
1388 4277 : if (degpol(W)>0) /*p=5 mod 6*/
1389 1302 : return Flx_equal1(Flxq_powu(W,3,T,p)) ? addii(q1,shifti(t,-1)):
1390 434 : subii(q1,shifti(t,-1));
1391 3409 : w = W[2];
1392 3409 : if (w==1) return N;
1393 2688 : if (w==p-1) return addii(q1,t);
1394 : else /*p=1 mod 6*/
1395 : {
1396 1918 : GEN u = shifti(t,-1), v = sqrtint(diviuexact(subii(q,sqri(u)),3));
1397 1918 : GEN a = addii(u,v), b = shifti(v,1);
1398 1918 : if (Fl_powu(w,3,p)==1)
1399 : {
1400 959 : if (Fl_add(umodiu(a,p),Fl_mul(w,umodiu(b,p),p),p)==0)
1401 434 : return subii(q1,subii(shifti(b,1),a));
1402 : else
1403 525 : return addii(q1,addii(a,b));
1404 : }
1405 : else
1406 : {
1407 959 : if (Fl_sub(umodiu(a,p),Fl_mul(w,umodiu(b,p),p),p)==0)
1408 434 : return subii(q1,subii(a,shifti(b,1)));
1409 : else
1410 525 : return subii(q1,addii(a,b));
1411 : }
1412 : }
1413 10150 : } else if (j==1728%p)
1414 : {
1415 : ulong w;
1416 : GEN W, N, t;
1417 5663 : if (mod4(q)==3) return q1;
1418 4263 : W = Flxq_pow(a4,shifti(q,-2),T,p);
1419 4263 : if (degpol(W)>0) return q1; /*p=3 mod 4*/
1420 3633 : w = W[2];
1421 3633 : N = Fp_ffellcard(gen_1,gen_0,q,n,utoipos(p));
1422 3633 : if(w==1) return N;
1423 2541 : t = subii(q1, N);
1424 2541 : if(w==p-1) return addii(q1, t);
1425 : else /*p=1 mod 4*/
1426 : {
1427 1386 : GEN u = shifti(t,-1), v = sqrtint(subii(q,sqri(u)));
1428 1386 : if (Fl_add(umodiu(u,p),Fl_mul(w,umodiu(v,p),p),p)==0)
1429 693 : return subii(q1,shifti(v,1));
1430 : else
1431 693 : return addii(q1,shifti(v,1));
1432 : }
1433 : } else
1434 : {
1435 4487 : ulong g = Fl_div(j, Fl_sub(1728%p, j, p), p);
1436 4487 : GEN N = Fp_ffellcard(utoi(Fl_triple(g,p)),utoi(Fl_double(g,p)),q,n,utoipos(p));
1437 4487 : GEN l = Flxq_mul(Flx_triple(a6,p),Flx_double(a4,p),T,p);
1438 4487 : if (Flxq_issquare(l,T,p)) return N;
1439 2940 : return subii(shifti(q1,1),N);
1440 : }
1441 : }
1442 :
1443 : static GEN
1444 454 : Flxq_ffellcard(GEN a4, GEN a6, GEN M, GEN q, GEN T, ulong p, long n)
1445 : {
1446 454 : long m = degpol(M);
1447 454 : GEN j = polx_Flx(M[1]);
1448 454 : GEN g = Flxq_div(j, mkvecsmall3(M[1],1728%p,p-1), M, p);
1449 454 : GEN N = Flxq_ellcard(Flx_triple(g, p), Flx_double(g, p), M, p);
1450 454 : GEN qm = powuu(p, m), q1 = addiu(q, 1), qm1 = addiu(qm, 1);
1451 454 : GEN l = Flxq_mul(Flx_triple(a6,p), Flx_double(a4,p), T, p);
1452 454 : GEN te = elltrace_extension(subii(qm1, N), n/m, qm);
1453 454 : return Flxq_issquare(l,T,p) ? subii(q1, te): addii(q1, te);
1454 : }
1455 :
1456 : static GEN
1457 62675 : Flxq_ellcard_i(GEN a4, GEN a6, GEN T, ulong p)
1458 : {
1459 62675 : long n = get_Flx_degree(T);
1460 62675 : GEN J, M, q = powuu(p, n);
1461 62675 : if (typ(a4)==t_VEC)
1462 20993 : return F3xq_ellcard(gel(a4,1), a6, T);
1463 41682 : if (p==3)
1464 9212 : return F3xq_ellcardj(a4, a6, T, q, n);
1465 32470 : if (degpol(a4)<=0 && degpol(a6)<=0)
1466 210 : return Fp_ffellcard(utoi(Flx_eval(a4,0,p)),utoi(Flx_eval(a6,0,p)),q,n,utoipos(p));
1467 32260 : J = Flxq_ellj(a4,a6,T,p);
1468 32260 : if (degpol(J)<=0)
1469 15827 : return Flxq_ellcardj(a4,a6,lgpol(J)?J[2]:0,T,q,p,n);
1470 16433 : M = Flxq_minpoly(J, T, p);
1471 16433 : if (degpol(M) < n)
1472 454 : return Flxq_ffellcard(a4, a6, M, q, T, p, n);
1473 15979 : if (p <= 7)
1474 10997 : return Flxq_ellcard_Satoh(a4, a6, J, T, p);
1475 4982 : if (cmpis(q,100)<0)
1476 0 : return utoi(Flxq_ellcard_naive(a4, a6, T, p));
1477 4982 : if (p == 13 || (7*p <= (ulong)10*n && (BITS_IN_LONG==64 || p <= 103)))
1478 154 : return Flxq_ellcard_Satoh(a4, a6, J, T, p);
1479 4828 : if (p <= (ulong)2*n)
1480 0 : return Flxq_ellcard_Kedlaya(a4, a6, T, p);
1481 4828 : if (expi(q)<=62)
1482 4774 : return Flxq_ellcard_Shanks(a4, a6, q, T, p);
1483 : else
1484 54 : return Fq_ellcard_SEA(Flx_to_ZX(a4),Flx_to_ZX(a6),q,Flx_to_ZX(T),utoipos(p),0);
1485 : }
1486 :
1487 : GEN
1488 62675 : Flxq_ellcard(GEN a4, GEN a6, GEN T, ulong p)
1489 : {
1490 62675 : pari_sp av = avma;
1491 62675 : return gc_INT(av, Flxq_ellcard_i(a4, a6, T, p));
1492 : }
1493 :
1494 : static long
1495 350 : Fl_ellj_trace(ulong j, ulong p)
1496 : {
1497 : ulong a4, a6;
1498 350 : Fl_ellj_to_a4a6(j, p, &a4, &a6);
1499 350 : return Fl_elltrace(a4, a6, p);
1500 : }
1501 :
1502 : /* Given ordinary E/Fq, a prime ell, and the height of the ell-volcano
1503 : * containing j(E) (= v_ell(conductor of Z[pi_E]) returns the height of j(E)
1504 : * on its ell-volcano (= v_ell(conductor of the order End(E)). */
1505 : static long
1506 147 : Fl_ellheightabovefloor(ulong j, long ell, long e, ulong p)
1507 : {
1508 147 : pari_sp av = avma;
1509 : GEN Xp, G, phi, phix, j0, j1;
1510 : long h, i, nj1;
1511 147 : if (e==0) return 0;
1512 147 : if (j==0 || j==1728%p) return e;
1513 119 : phi = ZXX_to_FlxX(polmodular_ZXX(ell, 0, 0, 1), p, 1);
1514 119 : phix = FlxY_evalx(phi, j, p);
1515 119 : Xp = Flx_Frobenius(phix, p);
1516 119 : G = Flx_gcd(Flx_sub(Xp, polx_Flx(0), p), phix, p);
1517 119 : nj1 = degpol(G);
1518 119 : if (nj1 < ell) return 0;
1519 112 : if (e==1 || nj1 != ell+1) return e;
1520 21 : j1 = Flx_roots(G, p);
1521 21 : nj1 = lg(j1)-1;
1522 21 : if (nj1 < 3) return 0;
1523 21 : j0 = mkvecsmall3(j,j,j);
1524 42 : for (h = 1; ; h++)
1525 133 : for(i = 1; i <= 3; i++)
1526 : {
1527 112 : GEN P = Flx_div_by_X_x(FlxY_evalx(phi, uel(j1,i), p), uel(j0,i), p, NULL);
1528 112 : GEN r = Flx_roots(P, p);
1529 112 : if (lg(r) == 1) return gc_long(av, h);
1530 91 : j0[i] = j1[i];
1531 91 : j1[i] = r[1];
1532 : }
1533 : }
1534 :
1535 : /* Given an ordinary elliptic curve E/Fp and an integer h, returns
1536 : * D = disc(End(E)) assuming h(D) = h, using the approach sketched in
1537 : * Remark 13. If the algorithm returns 0 it has proved that h(D) != h, but it
1538 : * is under no obligation to do so and is allowed to return any value when the
1539 : * assumption h(d) = h is false. */
1540 : static long
1541 350 : Fl_end13(ulong j, ulong h, ulong p)
1542 : {
1543 : ulong D0, v, h0;
1544 : long i, lL, lc, lD, nc;
1545 : GEN D, DF, cs, L, vP, vE;
1546 350 : ulong t = Fl_ellj_trace(j, p);
1547 :
1548 350 : D0 = coredisc2u_fact(factoru(4*p-t*t), -1, &vP, &vE);
1549 350 : h0 = itou(classno(stoi(-D0)));
1550 350 : if (h % h0 != 0) return 0;
1551 350 : h /= h0;
1552 350 : D = divisorsu_fact_factored(mkmat2(vP,vE));
1553 350 : DF = gel(D,2); D = gel(D,1);
1554 350 : lD = lg(D); v = D[lD-1];
1555 350 : cs = cgetg(lD,t_VECSMALL); nc = 0;
1556 1848 : for (i = 1; i < lD; i++)
1557 : {
1558 1498 : GEN F = gel(DF,i);
1559 1498 : ulong w = uquadclassnoF_fact(D0, -1, gel(F,1), gel(F,2));
1560 1498 : if (w == h) uel(cs,++nc) = v / uel(D,i);
1561 : }
1562 350 : if (nc==0) return 0;
1563 350 : if (nc==1) { v /= uel(cs,1); return -D0*v*v; }
1564 147 : L = cgetg(nc+1, t_VEC);
1565 448 : for (i = 1; i <= nc; i++) gel(L,i) = gel(factoru(uel(cs,i)), 1);
1566 147 : L = vecsmall_uniq(shallowconcat1(L));
1567 147 : lL = lg(L); lc = nc+1;
1568 147 : for (i = 1; i < lL; i++)
1569 : {
1570 147 : ulong ell = L[i];
1571 147 : long k, e = Fl_ellheightabovefloor(j, ell, z_lval(v,ell), p);
1572 448 : for (k = 1; k < lc; k++)
1573 301 : if(cs[k] && z_lval(cs[k], ell) != e) { cs[k] = 0; nc--; }
1574 147 : if (nc==0) return 0;
1575 147 : if (nc==1)
1576 : {
1577 161 : for (k = 1; k < lc; k++)
1578 161 : if (cs[k]) { v /= uel(cs,k); return -D0*v*v; }
1579 : }
1580 : }
1581 0 : return 0;
1582 : }
1583 :
1584 : INLINE int
1585 350 : RgX_is_monic_ZX(GEN pol)
1586 350 : { return RgX_is_ZX(pol) && ZX_is_monic(pol); }
1587 :
1588 : long
1589 357 : polisclass(GEN H)
1590 : {
1591 357 : pari_sp av = avma, btop;
1592 357 : long h = degpol(H), hl, i, pmin, vH = varn(H), vh;
1593 : double lmin;
1594 : ulong p;
1595 : GEN h2list;
1596 : forprime_t T;
1597 :
1598 357 : if (typ(H)!= t_POL) pari_err_TYPE("polsisclass",H);
1599 357 : if (h <= 0 || !RgX_is_monic_ZX(H)) return 0;
1600 350 : vh = vals(h);
1601 350 : h2list = cgetg(vh+2, t_VECSMALL); hl = 1;
1602 1071 : for (i = 0; i <= vh; i++)
1603 : {
1604 721 : ulong d = 1UL<<i;
1605 721 : if (((d-h)&1)==0) h2list[hl++] = d;
1606 : }
1607 350 : setlg(h2list, hl);
1608 350 : lmin = h * (log(log(h+2))+2);
1609 350 : pmin = 33 * ceil(lmin*lmin);
1610 350 : u_forprime_init(&T, pmin, ULONG_MAX);
1611 350 : btop = avma;
1612 3339 : while((p = u_forprime_next(&T)))
1613 : {
1614 : ulong r;
1615 : long D, nroots;
1616 3339 : GEN Xp, G, Hp = ZX_to_Flx(H,p);
1617 3339 : if (!Flx_is_squarefree(Hp, p)) { set_avma(btop); continue; }
1618 3339 : Xp = Flx_Frobenius(Hp, p);
1619 3339 : G = Flx_gcd(Flx_sub(Xp, polx_Flx(evalvarn(vH)), p), Hp, p);
1620 3339 : nroots = degpol(G);
1621 3339 : if (nroots==0) { set_avma(btop); continue; }
1622 1274 : if (nroots < h && !zv_search(h2list,nroots)) return gc_long(av, 0);
1623 1274 : r = Flx_oneroot(G, p);
1624 1274 : if (Fp_elljissupersingular(utoi(r), utoi(p))) { set_avma(btop); continue; }
1625 350 : D = Fl_end13(r, h, p);
1626 350 : if (D && gequal(H, polclass(stoi(D), 0, vH))) return gc_long(av, D);
1627 0 : return gc_long(av, 0);
1628 : }
1629 0 : pari_err_BUG("polisclass");
1630 : return 0; /* LCOV_EXCL_LINE */
1631 : }
|
