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 :
29 : /* Theses functions deal with point over elliptic curves over Fq defined
30 : * by an equation of the form y^2=x^3+a4*x+a6.
31 : * Most of the time a6 is omitted since it can be recovered from any point
32 : * on the curve.
33 : */
34 :
35 : GEN
36 63966 : RgE_to_FlxqE(GEN x, GEN T, ulong p)
37 : {
38 63966 : if (ell_is_inf(x)) return x;
39 63966 : retmkvec2(Rg_to_Flxq(gel(x,1),T,p),Rg_to_Flxq(gel(x,2),T,p));
40 : }
41 :
42 : GEN
43 154174 : FlxqE_changepoint(GEN x, GEN ch, GEN T, ulong p)
44 : {
45 154174 : pari_sp av = avma;
46 : GEN p1,z,u,r,s,t,v,v2,v3;
47 154174 : if (ell_is_inf(x)) return x;
48 91986 : u = gel(ch,1); r = gel(ch,2);
49 91986 : s = gel(ch,3); t = gel(ch,4);
50 91986 : v = Flxq_inv(u, T, p); v2 = Flxq_sqr(v, T, p); v3 = Flxq_mul(v,v2, T, p);
51 91986 : p1 = Flx_sub(gel(x,1),r, p);
52 91986 : z = cgetg(3,t_VEC);
53 91986 : gel(z,1) = Flxq_mul(v2, p1, T, p);
54 91986 : gel(z,2) = Flxq_mul(v3, Flx_sub(gel(x,2), Flx_add(Flxq_mul(s, p1, T, p),t, p), p), T, p);
55 91986 : return gerepileupto(av, z);
56 : }
57 :
58 : GEN
59 63966 : FlxqE_changepointinv(GEN x, GEN ch, GEN T, ulong p)
60 : {
61 : GEN u, r, s, t, X, Y, u2, u3, u2X, z;
62 63966 : if (ell_is_inf(x)) return x;
63 63966 : X = gel(x,1); Y = gel(x,2);
64 63966 : u = gel(ch,1); r = gel(ch,2);
65 63966 : s = gel(ch,3); t = gel(ch,4);
66 63966 : u2 = Flxq_sqr(u, T, p); u3 = Flxq_mul(u,u2, T, p);
67 63966 : u2X = Flxq_mul(u2,X, T, p);
68 63966 : z = cgetg(3, t_VEC);
69 63966 : gel(z,1) = Flx_add(u2X,r, p);
70 63966 : gel(z,2) = Flx_add(Flxq_mul(u3,Y, T, p), Flx_add(Flxq_mul(s,u2X, T, p), t, p), p);
71 63966 : return z;
72 : }
73 :
74 : static GEN
75 22834 : nonsquare_Flxq(GEN T, ulong p)
76 : {
77 22834 : pari_sp av = avma;
78 22834 : long n = degpol(T), vs = T[1];
79 : GEN a;
80 22834 : if (odd(n))
81 7686 : return mkvecsmall2(vs, nonsquare_Fl(p));
82 : do
83 : {
84 30821 : set_avma(av);
85 30821 : a = random_Flx(n, vs, p);
86 30821 : } while (Flxq_issquare(a, T, p));
87 15148 : return a;
88 : }
89 :
90 : void
91 22834 : Flxq_elltwist(GEN a, GEN a6, GEN T, ulong p, GEN *pt_a, GEN *pt_a6)
92 : {
93 22834 : GEN d = nonsquare_Flxq(T, p);
94 22834 : GEN d2 = Flxq_sqr(d, T, p), d3 = Flxq_mul(d2, d, T, p);
95 22834 : if (typ(a)==t_VECSMALL)
96 : {
97 15232 : *pt_a = Flxq_mul(a, d2, T, p);
98 15232 : *pt_a6 = Flxq_mul(a6, d3, T, p);
99 : } else
100 : {
101 7602 : *pt_a = mkvec(Flxq_mul(gel(a,1), d, T, p));
102 7602 : *pt_a6 = Flxq_mul(a6, d3, T, p);
103 : }
104 22834 : }
105 :
106 : static GEN
107 1288706 : FlxqE_dbl_slope(GEN P, GEN a4, GEN T, ulong p, GEN *slope)
108 : {
109 : GEN x, y, Q;
110 1288706 : if (ell_is_inf(P) || !lgpol(gel(P,2))) return ellinf();
111 1186009 : x = gel(P,1); y = gel(P,2);
112 1186009 : if (p==3UL)
113 534415 : *slope = typ(a4)==t_VEC ? Flxq_div(Flxq_mul(x, gel(a4, 1), T, p), y, T, p)
114 534415 : : Flxq_div(a4, Flx_neg(y, p), T, p);
115 : else
116 : {
117 651594 : GEN sx = Flx_add(Flx_triple(Flxq_sqr(x, T, p), p), a4, p);
118 651594 : *slope = Flxq_div(sx, Flx_double(y, p), T, p);
119 : }
120 1186009 : Q = cgetg(3,t_VEC);
121 1186009 : gel(Q, 1) = Flx_sub(Flxq_sqr(*slope, T, p), Flx_double(x, p), p);
122 1186009 : if (typ(a4)==t_VEC) gel(Q, 1) = Flx_sub(gel(Q, 1), gel(a4, 1), p);
123 1186009 : gel(Q, 2) = Flx_sub(Flxq_mul(*slope, Flx_sub(x, gel(Q, 1), p), T, p), y, p);
124 1186009 : return Q;
125 : }
126 :
127 : GEN
128 1259733 : FlxqE_dbl(GEN P, GEN a4, GEN T, ulong p)
129 : {
130 1259733 : pari_sp av = avma;
131 : GEN slope;
132 1259733 : return gerepileupto(av, FlxqE_dbl_slope(P,a4, T, p,&slope));
133 : }
134 :
135 : static GEN
136 538237 : FlxqE_add_slope(GEN P, GEN Q, GEN a4, GEN T, ulong p, GEN *slope)
137 : {
138 : GEN Px, Py, Qx, Qy, R;
139 538237 : if (ell_is_inf(P)) return Q;
140 534575 : if (ell_is_inf(Q)) return P;
141 534281 : Px = gel(P,1); Py = gel(P,2);
142 534281 : Qx = gel(Q,1); Qy = gel(Q,2);
143 534281 : if (Flx_equal(Px, Qx))
144 : {
145 47986 : if (Flx_equal(Py, Qy))
146 1413 : return FlxqE_dbl_slope(P, a4, T, p, slope);
147 : else
148 46573 : return ellinf();
149 : }
150 486295 : *slope = Flxq_div(Flx_sub(Py, Qy, p), Flx_sub(Px, Qx, p), T, p);
151 486295 : R = cgetg(3,t_VEC);
152 486295 : gel(R, 1) = Flx_sub(Flx_sub(Flxq_sqr(*slope, T, p), Px, p), Qx, p);
153 486295 : if (typ(a4)==t_VEC) gel(R, 1) = Flx_sub(gel(R, 1),gel(a4, 1), p);
154 486295 : gel(R, 2) = Flx_sub(Flxq_mul(*slope, Flx_sub(Px, gel(R, 1), p), T, p), Py, p);
155 486295 : return R;
156 : }
157 :
158 : GEN
159 533819 : FlxqE_add(GEN P, GEN Q, GEN a4, GEN T, ulong p)
160 : {
161 533819 : pari_sp av = avma;
162 : GEN slope;
163 533819 : return gerepileupto(av, FlxqE_add_slope(P,Q,a4, T, p,&slope));
164 : }
165 :
166 : static GEN
167 1391 : FlxqE_neg_i(GEN P, ulong p)
168 : {
169 1391 : if (ell_is_inf(P)) return P;
170 1391 : return mkvec2(gel(P,1), Flx_neg(gel(P,2), p));
171 : }
172 :
173 : GEN
174 490 : FlxqE_neg(GEN P, GEN T, ulong p)
175 : {
176 : (void) T;
177 490 : if (ell_is_inf(P)) return ellinf();
178 490 : return mkvec2(gcopy(gel(P,1)), Flx_neg(gel(P,2), p));
179 : }
180 :
181 : GEN
182 1391 : FlxqE_sub(GEN P, GEN Q, GEN a4, GEN T, ulong p)
183 : {
184 1391 : pari_sp av = avma;
185 : GEN slope;
186 1391 : return gerepileupto(av, FlxqE_add_slope(P, FlxqE_neg_i(Q, p), a4, T, p, &slope));
187 : }
188 :
189 : struct _FlxqE
190 : {
191 : GEN a4, a6;
192 : GEN T;
193 : ulong p;
194 : };
195 :
196 : static GEN
197 1259733 : _FlxqE_dbl(void *E, GEN P)
198 : {
199 1259733 : struct _FlxqE *ell = (struct _FlxqE *) E;
200 1259733 : return FlxqE_dbl(P, ell->a4, ell->T, ell->p);
201 : }
202 :
203 : static GEN
204 524988 : _FlxqE_add(void *E, GEN P, GEN Q)
205 : {
206 524988 : struct _FlxqE *ell=(struct _FlxqE *) E;
207 524988 : return FlxqE_add(P, Q, ell->a4, ell->T, ell->p);
208 : }
209 :
210 : static GEN
211 228296 : _FlxqE_mul(void *E, GEN P, GEN n)
212 : {
213 228296 : pari_sp av = avma;
214 228296 : struct _FlxqE *e=(struct _FlxqE *) E;
215 228296 : long s = signe(n);
216 228296 : if (!s || ell_is_inf(P)) return ellinf();
217 228078 : if (s<0) P = FlxqE_neg(P, e->T, e->p);
218 228078 : if (is_pm1(n)) return s>0? gcopy(P): P;
219 221464 : return gerepilecopy(av, gen_pow_i(P, n, e, &_FlxqE_dbl, &_FlxqE_add));
220 : }
221 :
222 : GEN
223 65437 : FlxqE_mul(GEN P, GEN n, GEN a4, GEN T, ulong p)
224 : {
225 : struct _FlxqE E;
226 65437 : E.a4= a4; E.T = T; E.p = p;
227 65437 : return _FlxqE_mul(&E, P, n);
228 : }
229 :
230 : /* 3*x^2+2*a2*x = -a2*x, and a2!=0 */
231 :
232 : /* Finds a random nonsingular point on E */
233 : static GEN
234 77546 : random_F3xqE(GEN a2, GEN a6, GEN T)
235 : {
236 77546 : pari_sp ltop = avma;
237 : GEN x, y, rhs;
238 77546 : const ulong p = 3;
239 : do
240 : {
241 154728 : set_avma(ltop);
242 154728 : x = random_Flx(get_Flx_degree(T),get_Flx_var(T),p);
243 154728 : rhs = Flx_add(Flxq_mul(Flxq_sqr(x, T, p), Flx_add(x, a2, p), T, p), a6, p);
244 154728 : } while ((!lgpol(rhs) && !lgpol(x)) || !Flxq_issquare(rhs, T, p));
245 77546 : y = Flxq_sqrt(rhs, T, p);
246 77546 : if (!y) pari_err_PRIME("random_F3xqE", T);
247 77546 : return gerepilecopy(ltop, mkvec2(x, y));
248 : }
249 :
250 : /* Finds a random nonsingular point on E */
251 : GEN
252 145336 : random_FlxqE(GEN a4, GEN a6, GEN T, ulong p)
253 : {
254 145336 : pari_sp ltop = avma;
255 : GEN x, x2, y, rhs;
256 145336 : if (typ(a4)==t_VEC)
257 77546 : return random_F3xqE(gel(a4,1), a6, T);
258 : do
259 : {
260 132842 : set_avma(ltop);
261 132842 : x = random_Flx(get_Flx_degree(T),get_Flx_var(T),p);
262 132842 : x2 = Flxq_sqr(x, T, p); /* x^3+a4*x+a6 = x*(x^2+a4)+a6 */
263 132842 : rhs = Flx_add(Flxq_mul(x, Flx_add(x2, a4, p), T, p), a6, p);
264 134809 : } while ((!lgpol(rhs) && !lgpol(Flx_add(Flx_triple(x2, p), a4, p)))
265 134767 : || !Flxq_issquare(rhs, T, p));
266 67790 : y = Flxq_sqrt(rhs, T, p);
267 67790 : if (!y) pari_err_PRIME("random_FlxqE", T);
268 67790 : return gerepilecopy(ltop, mkvec2(x, y));
269 : }
270 :
271 : static GEN
272 70212 : _FlxqE_rand(void *E)
273 : {
274 70212 : struct _FlxqE *ell=(struct _FlxqE *) E;
275 70212 : return random_FlxqE(ell->a4, ell->a6, ell->T, ell->p);
276 : }
277 :
278 : static const struct bb_group FlxqE_group={_FlxqE_add,_FlxqE_mul,_FlxqE_rand,hash_GEN,zvV_equal,ell_is_inf, NULL};
279 :
280 : const struct bb_group *
281 34 : get_FlxqE_group(void ** pt_E, GEN a4, GEN a6, GEN T, ulong p)
282 : {
283 34 : struct _FlxqE *e = (struct _FlxqE *) stack_malloc(sizeof(struct _FlxqE));
284 34 : e->a4 = a4; e->a6 = a6; e->T = Flx_get_red(T, p); e->p = p;
285 34 : *pt_E = (void *) e;
286 34 : return &FlxqE_group;
287 : }
288 :
289 : GEN
290 2729 : FlxqE_order(GEN z, GEN o, GEN a4, GEN T, ulong p)
291 : {
292 2729 : pari_sp av = avma;
293 : struct _FlxqE e;
294 2729 : e.a4=a4; e.T=T; e.p=p;
295 2729 : return gerepileuptoint(av, gen_order(z, o, (void*)&e, &FlxqE_group));
296 : }
297 :
298 : GEN
299 49 : FlxqE_log(GEN a, GEN b, GEN o, GEN a4, GEN T, ulong p)
300 : {
301 49 : pari_sp av = avma;
302 : struct _FlxqE e;
303 49 : e.a4=a4; e.T=T; e.p=p;
304 49 : return gerepileuptoint(av, gen_PH_log(a, b, o, (void*)&e, &FlxqE_group));
305 : }
306 :
307 : /***********************************************************************/
308 : /** **/
309 : /** Pairings **/
310 : /** **/
311 : /***********************************************************************/
312 :
313 : /* Derived from APIP from and by Jerome Milan, 2012 */
314 :
315 : static GEN
316 67918 : FlxqE_vert(GEN P, GEN Q, GEN a4, GEN T, ulong p)
317 : {
318 67918 : long vT = get_Flx_var(T);
319 : GEN df;
320 67918 : if (ell_is_inf(P))
321 20627 : return pol1_Flx(vT);
322 47291 : if (!Flx_equal(gel(Q, 1), gel(P, 1)))
323 45480 : return Flx_sub(gel(Q, 1), gel(P, 1), p);
324 1811 : if (lgpol(gel(P,2))!=0) return pol1_Flx(vT);
325 672 : df = typ(a4)==t_VEC ? Flxq_mul(gel(P,1), Flx_mulu(gel(a4, 1), 2, p), T, p)
326 1142 : : a4;
327 1142 : return Flxq_inv(Flx_add(Flx_mulu(Flxq_sqr(gel(P,1), T, p), 3, p),
328 : df, p), T, p);
329 : }
330 :
331 : static GEN
332 30587 : FlxqE_Miller_line(GEN R, GEN Q, GEN slope, GEN a4, GEN T, ulong p)
333 : {
334 30587 : long vT = get_Flx_var(T);
335 30587 : GEN x = gel(Q, 1), y = gel(Q, 2);
336 30587 : GEN tmp1 = Flx_sub(x, gel(R, 1), p);
337 30587 : GEN tmp2 = Flx_add(Flxq_mul(tmp1, slope, T, p), gel(R, 2), p);
338 30587 : if (!Flx_equal(y, tmp2))
339 29452 : return Flx_sub(y, tmp2, p);
340 1135 : if (lgpol(y) == 0)
341 571 : return pol1_Flx(vT);
342 : else
343 : {
344 564 : GEN s1, s2, a2 = typ(a4)==t_VEC ? gel(a4,1): NULL;
345 564 : GEN y2i = Flxq_inv(Flx_mulu(y, 2, p), T, p);
346 564 : GEN df = a2 ? Flxq_mul(x, Flx_mulu(a2, 2, p), T, p): a4;
347 : GEN x3, ddf;
348 564 : s1 = Flxq_mul(Flx_add(Flx_mulu(Flxq_sqr(x, T, p), 3, p), df, p), y2i, T, p);
349 564 : if (!Flx_equal(s1, slope))
350 330 : return Flx_sub(s1, slope, p);
351 234 : x3 = Flx_mulu(x, 3, p);
352 234 : ddf = a2 ? Flx_add(x3, a2, p): x3;
353 234 : s2 = Flxq_mul(Flx_sub(ddf, Flxq_sqr(s1, T, p), p), y2i, T, p);
354 234 : return lgpol(s2)!=0 ? s2: y2i;
355 : }
356 : }
357 :
358 : /* Computes the equation of the line tangent to R and returns its
359 : * evaluation at the point Q. Also doubles the point R. */
360 : static GEN
361 46863 : FlxqE_tangent_update(GEN R, GEN Q, GEN a4, GEN T, ulong p, GEN *pt_R)
362 : {
363 46863 : if (ell_is_inf(R))
364 : {
365 4077 : *pt_R = ellinf();
366 4077 : return pol1_Flx(get_Flx_var(T));
367 : }
368 42786 : else if (!lgpol(gel(R,2)))
369 : {
370 15226 : *pt_R = ellinf();
371 15226 : return FlxqE_vert(R, Q, a4, T, p);
372 : } else {
373 : GEN slope;
374 27560 : *pt_R = FlxqE_dbl_slope(R, a4, T, p, &slope);
375 27560 : return FlxqE_Miller_line(R, Q, slope, a4, T, p);
376 : }
377 : }
378 :
379 : /* Computes the equation of the line through R and P, and returns its
380 : * evaluation at the point Q. Also adds P to the point R. */
381 : static GEN
382 4435 : FlxqE_chord_update(GEN R, GEN P, GEN Q, GEN a4, GEN T, ulong p, GEN *pt_R)
383 : {
384 4435 : if (ell_is_inf(R))
385 : {
386 77 : *pt_R = gcopy(P);
387 77 : return FlxqE_vert(P, Q, a4, T, p);
388 : }
389 4358 : else if (ell_is_inf(P))
390 : {
391 0 : *pt_R = gcopy(R);
392 0 : return FlxqE_vert(R, Q, a4, T, p);
393 : }
394 4358 : else if (Flx_equal(gel(P, 1), gel(R, 1)))
395 : {
396 1331 : if (Flx_equal(gel(P, 2), gel(R, 2)))
397 7 : return FlxqE_tangent_update(R, Q, a4, T, p, pt_R);
398 : else
399 : {
400 1324 : *pt_R = ellinf();
401 1324 : return FlxqE_vert(R, Q, a4, T, p);
402 : }
403 : } else {
404 : GEN slope;
405 3027 : *pt_R = FlxqE_add_slope(P, R, a4, T, p, &slope);
406 3027 : return FlxqE_Miller_line(R, Q, slope, a4, T, p);
407 : }
408 : }
409 :
410 : struct _FlxqE_miller
411 : {
412 : ulong p;
413 : GEN T, a4, P;
414 : };
415 :
416 : static GEN
417 46856 : FlxqE_Miller_dbl(void* E, GEN d)
418 : {
419 46856 : struct _FlxqE_miller *m = (struct _FlxqE_miller *)E;
420 46856 : ulong p = m->p;
421 46856 : GEN T = m->T, a4 = m->a4, P = m->P;
422 46856 : GEN v, line, point = gel(d,3);
423 46856 : GEN N = Flxq_sqr(gel(d,1), T, p);
424 46856 : GEN D = Flxq_sqr(gel(d,2), T, p);
425 46856 : line = FlxqE_tangent_update(point, P, a4, T, p, &point);
426 46856 : N = Flxq_mul(N, line, T, p);
427 46856 : v = FlxqE_vert(point, P, a4, T, p);
428 46856 : D = Flxq_mul(D, v, T, p); return mkvec3(N, D, point);
429 : }
430 :
431 : static GEN
432 4435 : FlxqE_Miller_add(void* E, GEN va, GEN vb)
433 : {
434 4435 : struct _FlxqE_miller *m = (struct _FlxqE_miller *)E;
435 4435 : ulong p = m->p;
436 4435 : GEN T = m->T, a4 = m->a4, P = m->P;
437 : GEN v, line, point;
438 4435 : GEN na = gel(va,1), da = gel(va,2), pa = gel(va,3);
439 4435 : GEN nb = gel(vb,1), db = gel(vb,2), pb = gel(vb,3);
440 4435 : GEN N = Flxq_mul(na, nb, T, p);
441 4435 : GEN D = Flxq_mul(da, db, T, p);
442 4435 : line = FlxqE_chord_update(pa, pb, P, a4, T, p, &point);
443 4435 : N = Flxq_mul(N, line, T, p);
444 4435 : v = FlxqE_vert(point, P, a4, T, p);
445 4435 : D = Flxq_mul(D, v, T, p); return mkvec3(N, D, point);
446 : }
447 :
448 : /* Returns the Miller function f_{m, Q} evaluated at the point P using
449 : * the standard Miller algorithm. */
450 : static GEN
451 16473 : FlxqE_Miller(GEN Q, GEN P, GEN m, GEN a4, GEN T, ulong p)
452 : {
453 16473 : pari_sp av = avma;
454 : struct _FlxqE_miller d;
455 : GEN v, N, D, g1;
456 :
457 16473 : d.a4 = a4; d.T = T; d.p = p; d.P = P;
458 16473 : g1 = pol1_Flx(get_Flx_var(T));
459 16473 : v = gen_pow_i(mkvec3(g1,g1,Q), m, (void*)&d,
460 : FlxqE_Miller_dbl, FlxqE_Miller_add);
461 16473 : N = gel(v,1); D = gel(v,2);
462 16473 : return gerepileupto(av, Flxq_div(N, D, T, p));
463 : }
464 :
465 : GEN
466 12792 : FlxqE_weilpairing(GEN P, GEN Q, GEN m, GEN a4, GEN T, ulong p)
467 : {
468 12792 : pari_sp av = avma;
469 : GEN N, D, result;
470 12792 : if (ell_is_inf(P) || ell_is_inf(Q)
471 10085 : || (Flx_equal(gel(P,1),gel(Q,1)) && Flx_equal(gel(P,2),gel(Q,2))))
472 4587 : return pol1_Flx(get_Flx_var(T));
473 8205 : N = FlxqE_Miller(P, Q, m, a4, T, p);
474 8205 : D = FlxqE_Miller(Q, P, m, a4, T, p);
475 8205 : result = Flxq_div(N, D, T, p);
476 8205 : if (mpodd(m)) result = Flx_neg(result, p);
477 8205 : return gerepileupto(av, result);
478 : }
479 :
480 : GEN
481 63 : FlxqE_tatepairing(GEN P, GEN Q, GEN m, GEN a4, GEN T, ulong p)
482 : {
483 63 : if (ell_is_inf(P) || ell_is_inf(Q)) return pol1_Flx(get_Flx_var(T));
484 63 : return FlxqE_Miller(P, Q, m, a4, T, p);
485 : }
486 :
487 : static GEN
488 12771 : _FlxqE_pairorder(void *E, GEN P, GEN Q, GEN m, GEN F)
489 : {
490 12771 : struct _FlxqE *e = (struct _FlxqE *) E;
491 12771 : return Flxq_order(FlxqE_weilpairing(P,Q,m,e->a4,e->T,e->p), F, e->T, e->p);
492 : }
493 :
494 : GEN
495 15609 : Flxq_ellgroup(GEN a4, GEN a6, GEN N, GEN T, ulong p, GEN *pt_m)
496 : {
497 : struct _FlxqE e;
498 15609 : GEN q = powuu(p, get_Flx_degree(T));
499 15609 : e.a4=a4; e.a6=a6; e.T=T; e.p=p;
500 15609 : return gen_ellgroup(N, subiu(q,1), pt_m, (void*)&e, &FlxqE_group, _FlxqE_pairorder);
501 : }
502 :
503 : GEN
504 14363 : Flxq_ellgens(GEN a4, GEN a6, GEN ch, GEN D, GEN m, GEN T, ulong p)
505 : {
506 : GEN P;
507 14363 : pari_sp av = avma;
508 : struct _FlxqE e;
509 14363 : e.a4=a4; e.a6=a6; e.T=T; e.p=p;
510 14363 : switch(lg(D)-1)
511 : {
512 63 : case 0:
513 63 : return cgetg(1,t_VEC);
514 11794 : case 1:
515 11794 : P = gen_gener(gel(D,1), (void*)&e, &FlxqE_group);
516 11794 : P = mkvec(FlxqE_changepoint(P, ch, T, p));
517 11794 : break;
518 2506 : default:
519 2506 : P = gen_ellgens(gel(D,1), gel(D,2), m, (void*)&e, &FlxqE_group, _FlxqE_pairorder);
520 2506 : gel(P,1) = FlxqE_changepoint(gel(P,1), ch, T, p);
521 2506 : gel(P,2) = FlxqE_changepoint(gel(P,2), ch, T, p);
522 2506 : break;
523 : }
524 14300 : return gerepilecopy(av, P);
525 : }
526 : /***********************************************************************/
527 : /** **/
528 : /** Point counting **/
529 : /** **/
530 : /***********************************************************************/
531 :
532 : /* assume a and n are coprime */
533 : static GEN
534 76251 : RgX_circular_shallow(GEN P, long a, long n)
535 : {
536 76251 : long i, l = lgpol(P);
537 76251 : GEN Q = cgetg(2+n,t_POL);
538 76251 : Q[1] = P[1];
539 512330 : for(i=0; i<l; i++)
540 436079 : gel(Q,2+(i*a)%n) = gel(P,2+i);
541 168693 : for( ; i<n; i++)
542 92442 : gel(Q,2+(i*a)%n) = gen_0;
543 76251 : return normalizepol_lg(Q,2+n);
544 : }
545 :
546 : static GEN
547 76251 : ZpXQ_frob_cyc(GEN x, GEN T, GEN q, ulong p)
548 : {
549 76251 : long n = get_FpX_degree(T);
550 76251 : return FpX_rem(RgX_circular_shallow(x,p,n+1), T, q);
551 : }
552 :
553 : static GEN
554 113526 : ZpXQ_frob(GEN x, GEN Xm, GEN T, GEN q, ulong p)
555 : {
556 113526 : if (lg(Xm)==1)
557 43428 : return ZpXQ_frob_cyc(x, T, q, p);
558 : else
559 : {
560 70098 : long n = get_FpX_degree(T);
561 70098 : GEN V = RgX_blocks(RgX_inflate(x, p), n, p);
562 70098 : GEN W = ZXV_dotproduct(V, Xm);
563 70098 : return FpX_rem(W, T, q);
564 : }
565 : }
566 :
567 : struct _lift_lin
568 : {
569 : ulong p;
570 : GEN sqx, Tp;
571 : GEN ai, Xm;
572 : };
573 :
574 84035 : static GEN _lift_invl(void *E, GEN x)
575 : {
576 84035 : struct _lift_lin *d = (struct _lift_lin *) E;
577 84035 : GEN T = d->Tp;
578 84035 : ulong p = d->p;
579 84035 : GEN xai = Flxq_mul(ZX_to_Flx(x, p), d->ai, T, p);
580 84035 : return Flx_to_ZX(Flxq_lroot_fast(xai, d->sqx, T, p));
581 : }
582 :
583 23744 : static GEN _lift_lin(void *E, GEN F, GEN x2, GEN q)
584 : {
585 23744 : struct _lift_lin *d = (struct _lift_lin *) E;
586 23744 : pari_sp av = avma;
587 23744 : GEN T = gel(F,3), Xm = gel(F,4);
588 23744 : GEN y2 = ZpXQ_frob(x2, Xm, T, q, d->p);
589 23744 : GEN lin = FpX_add(ZX_mul(gel(F,1), y2), ZX_mul(gel(F,2), x2), q);
590 23744 : return gerepileupto(av, FpX_rem(lin, T, q));
591 : }
592 :
593 : static GEN
594 180915 : FpM_FpXV_bilinear(GEN P, GEN X, GEN Y, GEN p)
595 : {
596 180915 : pari_sp av = avma;
597 180915 : GEN s = ZX_mul(FpXV_FpC_mul(X,gel(P,1),p),gel(Y,1));
598 180915 : long i, l = lg(P);
599 850311 : for(i=2; i<l; i++)
600 669396 : s = ZX_add(s, ZX_mul(FpXV_FpC_mul(X,gel(P,i),p),gel(Y,i)));
601 180915 : return gerepileupto(av, FpX_red(s, p));
602 : }
603 :
604 : static GEN
605 180915 : FpM_FpXQV_bilinear(GEN P, GEN X, GEN Y, GEN T, GEN p)
606 : {
607 180915 : return FpX_rem(FpM_FpXV_bilinear(P,X,Y,p),T,p);
608 : }
609 :
610 : static GEN
611 120624 : FpXC_powderiv(GEN M, GEN p)
612 : {
613 : long i, l;
614 120624 : long v = varn(gel(M,2));
615 120624 : GEN m = cgetg_copy(M, &l);
616 120624 : gel(m,1) = pol_0(v);
617 120624 : gel(m,2) = pol_1(v);
618 446432 : for(i=2; i<l-1; i++)
619 325808 : gel(m,i+1) = FpX_Fp_mul(gel(M,i),utoi(i), p);
620 120624 : return m;
621 : }
622 :
623 : struct _lift_iso
624 : {
625 : GEN phi;
626 : GEN Xm,T;
627 : GEN sqx, Tp;
628 : ulong p;
629 : };
630 :
631 : static GEN
632 60291 : _lift_iter(void *E, GEN x2, GEN q)
633 : {
634 60291 : struct _lift_iso *d = (struct _lift_iso *) E;
635 60291 : ulong p = d->p;
636 60291 : long n = lg(d->phi)-2;
637 60291 : GEN TN = FpXT_red(d->T, q), XN = FpXV_red(d->Xm, q);
638 60291 : GEN y2 = ZpXQ_frob(x2, XN, TN, q, p);
639 60291 : GEN xp = FpXQ_powers(x2, n, TN, q);
640 60291 : GEN yp = FpXQ_powers(y2, n, TN, q);
641 60291 : GEN V = FpM_FpXQV_bilinear(d->phi,xp,yp,TN,q);
642 60291 : return mkvec3(V,xp,yp);
643 : }
644 :
645 : static GEN
646 60291 : _lift_invd(void *E, GEN V, GEN v, GEN qM, long M)
647 : {
648 60291 : struct _lift_iso *d = (struct _lift_iso *) E;
649 : struct _lift_lin e;
650 60291 : ulong p = d->p;
651 60291 : GEN TM = FpXT_red(d->T, qM), XM = FpXV_red(d->Xm, qM);
652 60291 : GEN xp = FpXV_red(gel(v,2), qM);
653 60291 : GEN yp = FpXV_red(gel(v,3), qM);
654 60291 : GEN Dx = FpM_FpXQV_bilinear(d->phi, FpXC_powderiv(xp, qM), yp, TM, qM);
655 60291 : GEN Dy = FpM_FpXQV_bilinear(d->phi, xp, FpXC_powderiv(yp, qM), TM, qM);
656 60291 : GEN F = mkvec4(Dy, Dx, TM, XM);
657 60291 : e.ai = Flxq_inv(ZX_to_Flx(Dy,p),d->Tp,p);
658 60291 : e.sqx = d->sqx; e.Tp = d->Tp; e.p=p; e.Xm = XM;
659 60291 : return gen_ZpX_Dixon(F,V,qM,utoipos(p),M,(void*) &e, _lift_lin, _lift_invl);
660 : }
661 :
662 : static GEN
663 25032 : lift_isogeny(GEN phi, GEN x0, long n, GEN Xm, GEN T, GEN sqx, GEN Tp, ulong p)
664 : {
665 : struct _lift_iso d;
666 25032 : d.phi=phi;
667 25032 : d.Xm=Xm; d.T=T;
668 25032 : d.sqx=sqx; d.Tp=Tp; d.p=p;
669 25032 : return gen_ZpX_Newton(x0, utoipos(p), n,(void*)&d, _lift_iter, _lift_invd);
670 : }
671 :
672 : static GEN
673 25011 : getc2(GEN act, GEN X, GEN T, GEN q, ulong p, long N)
674 : {
675 25011 : GEN A1 = RgV_to_RgX(gel(act,1),0), A2 = RgV_to_RgX(gel(act,2),0);
676 25011 : long n = brent_kung_optpow(maxss(degpol(A1),degpol(A2)),2,1);
677 25011 : GEN xp = FpXQ_powers(X,n,T,q);
678 25011 : GEN P = FpX_FpXQV_eval(A1, xp, T, q);
679 25011 : GEN Q = FpX_FpXQV_eval(A2, xp, T, q);
680 25011 : return ZpXQ_div(P, Q, T, q, utoipos(p), N);
681 : }
682 :
683 : struct _ZpXQ_norm
684 : {
685 : long n;
686 : GEN T, p;
687 : };
688 :
689 : static GEN
690 32823 : ZpXQ_norm_mul(void *E, GEN x, GEN y)
691 : {
692 32823 : struct _ZpXQ_norm *D = (struct _ZpXQ_norm*)E;
693 32823 : GEN P = gel(x,1), Q = gel(y,1);
694 32823 : long a = mael(x,2,1), b = mael(y,2,1);
695 32823 : retmkvec2(FpXQ_mul(P,ZpXQ_frob_cyc(Q, D->T, D->p, a), D->T, D->p),
696 : mkvecsmall((a*b)%D->n));
697 : }
698 :
699 : static GEN
700 22715 : ZpXQ_norm_sqr(void *E, GEN x)
701 : {
702 22715 : return ZpXQ_norm_mul(E, x, x);
703 : }
704 :
705 : /* Assume T = Phi_(n) and n prime */
706 : GEN
707 11340 : ZpXQ_norm_pcyc(GEN x, GEN T, GEN q, GEN p)
708 : {
709 : GEN z;
710 : struct _ZpXQ_norm D;
711 11340 : long d = get_FpX_degree(T);
712 11340 : D.T = T; D.p = q; D.n = d+1;
713 11340 : if (d==1) return ZX_copy(x);
714 11340 : z = mkvec2(x,mkvecsmall(p[2]));
715 11340 : z = gen_powu_i(z,d,(void*)&D,ZpXQ_norm_sqr,ZpXQ_norm_mul);
716 11340 : return gmael(z,1,2);
717 : }
718 :
719 : /* Assume T = Phi_(n) and n prime */
720 : static GEN
721 11102 : ZpXQ_sqrtnorm_pcyc(GEN x, GEN T, GEN q, GEN p, long e)
722 : {
723 11102 : GEN z = ZpXQ_norm_pcyc(x, T, q, p);
724 11102 : return Zp_sqrtlift(z,Fp_sqrt(z,p),p,e);
725 : }
726 :
727 : /* Assume a = 1 [p], return the square root of the norm */
728 : static GEN
729 13930 : ZpXQ_sqrtnorm(GEN a, GEN T, GEN q, GEN p, long e)
730 : {
731 13930 : GEN s = Fp_div(FpXQ_trace(ZpXQ_log(a, T, p, e), T, q), gen_2, q);
732 13930 : return modii(gel(Qp_exp(cvtop(s, p, e-1)),4), q);
733 : }
734 :
735 : struct _teich_lin
736 : {
737 : ulong p;
738 : GEN sqx, Tp;
739 : long m;
740 : };
741 :
742 : static GEN
743 29470 : _teich_invl(void *E, GEN x)
744 : {
745 29470 : struct _teich_lin *d = (struct _teich_lin *) E;
746 29470 : ulong p = d->p;
747 29470 : GEN T = d->Tp;
748 29470 : return Flx_to_ZX(Flxq_lroot_fast(ZX_to_Flx(x, p), d->sqx, T, p));
749 : }
750 :
751 : static GEN
752 8953 : _teich_lin(void *E, GEN F, GEN x2, GEN q)
753 : {
754 8953 : struct _teich_lin *d = (struct _teich_lin *) E;
755 8953 : pari_sp av = avma;
756 8953 : GEN T = gel(F,2), Xm = gel(F,3);
757 8953 : GEN y2 = ZpXQ_frob(x2, Xm, T, q, d->p);
758 8953 : GEN lin = FpX_sub(y2, ZX_mulu(ZX_mul(gel(F,1), x2), d->p), q);
759 8953 : return gerepileupto(av, FpX_rem(lin, T, q));
760 : }
761 :
762 : struct _teich_iso
763 : {
764 : GEN Xm, T;
765 : GEN sqx, Tp;
766 : ulong p;
767 : };
768 :
769 : static GEN
770 20517 : _teich_iter(void *E, GEN x2, GEN q)
771 : {
772 20517 : struct _teich_iso *d = (struct _teich_iso *) E;
773 20517 : ulong p = d->p;
774 20517 : GEN TN = FpXT_red(d->T, q), XN = FpXV_red(d->Xm, q);
775 20517 : GEN y2 = ZpXQ_frob(x2, XN, TN, q, d->p);
776 20517 : GEN x1 = FpXQ_powu(x2, p-1, TN, q);
777 20517 : GEN xp = FpXQ_mul(x2, x1, TN, q);
778 20517 : GEN V = FpX_sub(y2,xp,q);
779 20517 : return mkvec2(V,x1);
780 : }
781 :
782 : static GEN
783 20517 : _teich_invd(void *E, GEN V, GEN v, GEN qM, long M)
784 : {
785 20517 : struct _teich_iso *d = (struct _teich_iso *) E;
786 : struct _teich_lin e;
787 20517 : ulong p = d->p;
788 20517 : GEN TM = FpXT_red(d->T, qM), XM = FpXV_red(d->Xm, qM);
789 20517 : GEN x1 = FpX_red(gel(v,2), qM);
790 20517 : GEN F = mkvec3(x1, TM, XM);
791 20517 : e.sqx = d->sqx; e.Tp = d->Tp; e.p=p;
792 20517 : return gen_ZpX_Dixon(F,V,qM,utoipos(p),M,(void*) &e, _teich_lin, _teich_invl);
793 : }
794 :
795 : static GEN
796 10213 : Teichmuller_lift(GEN x, GEN Xm, GEN T, GEN sqx, GEN Tp, ulong p, long N)
797 : {
798 : struct _teich_iso d;
799 10213 : d.Xm = Xm; d.T = T; d.sqx = sqx; d.Tp = Tp; d.p = p;
800 10213 : return gen_ZpX_Newton(x,utoipos(p), N,(void*)&d, _teich_iter, _teich_invd);
801 : }
802 :
803 : static GEN
804 25032 : get_norm(GEN a4, GEN a6, GEN T, ulong p, long N)
805 : {
806 25032 : long sv=T[1];
807 : GEN a;
808 25032 : if (p==3) a = gel(a4,1);
809 : else
810 : {
811 10227 : GEN P = mkpoln(4, pol1_Flx(sv), pol0_Flx(sv), a4, a6);
812 10227 : a = gel(FlxqX_powu(P,p>>1,T,p),2+p-1);
813 : }
814 25032 : return Zp_sqrtnlift(gen_1,subss(p,1),utoi(Flxq_norm(a,T,p)),utoipos(p), N);
815 : }
816 :
817 : static GEN
818 25011 : fill_pols(long n, const long *v, long m, const long *vn,
819 : const long *vd, GEN *act)
820 : {
821 : long i, j;
822 25011 : long d = upowuu(n,12/(n-1));
823 25011 : GEN N, D, M = zeromatcopy(n+1,n+1);
824 25011 : gmael(M,1,n+1) = gen_1;
825 120568 : for(i=2;i<=n+1;i++)
826 338373 : for(j=i-1;j<=n;j++)
827 242816 : gmael(M,i,j) = mulis(powuu(d,i-2),v[j-i+1]);
828 25011 : N = cgetg(m+1,t_COL);
829 25011 : D = cgetg(m+1,t_COL);
830 135359 : for(i=1;i<=m;i++)
831 : {
832 110348 : gel(N,i) = stoi(*vn++);
833 110348 : gel(D,i) = stoi(*vd++);
834 : }
835 25011 : *act = mkmat2(N,D);
836 25011 : return M;
837 : }
838 :
839 : /*
840 : These polynomials were extracted from the ECHIDNA databases
841 : available at <http://echidna.maths.usyd.edu.au/echidna/>
842 : and computed by David R. Kohel.
843 : Return the matrix of the modular polynomial, set act to the parametrization,
844 : and set dj to the opposite of the supersingular j-invariant.
845 : */
846 : static GEN
847 25011 : get_Kohel_polynomials(ulong p, GEN *act, long *dj)
848 : {
849 25011 : const long mat3[] = {-1,-36,-270};
850 25011 : const long num3[] = {1,-483,-21141,-59049};
851 25011 : const long den3[] = {1,261, 4347, -6561};
852 25011 : const long mat5[] = {-1,-30,-315,-1300,-1575};
853 25011 : const long num5[] = {-1,490,20620,158750,78125};
854 25011 : const long den5[] = {-1,-254,-4124,-12250,3125};
855 25011 : const long mat7[] = {-1,-28,-322,-1904,-5915,-8624,-4018};
856 25011 : const long num7[] = {1,-485,-24058,-343833,-2021642,-4353013,-823543};
857 25011 : const long den7[] = {1,259,5894,49119,168406,166355,-16807};
858 25011 : const long mat13[]= {-1,-26,-325,-2548,-13832,-54340,-157118,-333580,-509366,
859 : -534820,-354536,-124852,-15145};
860 25011 : const long num13[]= {1,-487,-24056,-391463,-3396483,-18047328,-61622301,
861 : -133245853,-168395656,-95422301,-4826809};
862 25011 : const long den13[]= {1,257,5896,60649,364629,1388256,3396483,5089019,4065464,
863 : 1069939,-28561};
864 25011 : switch(p)
865 : {
866 14805 : case 3:
867 14805 : *dj = 0;
868 14805 : return fill_pols(3,mat3,4,num3,den3,act);
869 10171 : case 5:
870 10171 : *dj = 0;
871 10171 : return fill_pols(5,mat5,5,num5,den5,act);
872 28 : case 7:
873 28 : *dj = 1;
874 28 : return fill_pols(7,mat7,7,num7,den7,act);
875 7 : case 13:
876 7 : *dj = 8;
877 7 : return fill_pols(13,mat13,11,num13,den13,act);
878 : }
879 : *dj=0; *act = NULL; return NULL; /* LCOV_EXCL_LINE */
880 : }
881 :
882 : long
883 32476 : zx_is_pcyc(GEN T)
884 : {
885 32476 : long i, n = degpol(T);
886 32476 : if (!uisprime(n+1))
887 11840 : return 0;
888 99148 : for (i=0; i<=n; i++)
889 87808 : if (T[i+2]!=1UL)
890 9296 : return 0;
891 11340 : return 1;
892 : }
893 :
894 : static GEN
895 25011 : Flxq_ellcard_Kohel(GEN a4, GEN a6, GEN T, ulong p)
896 : {
897 25011 : pari_sp av = avma, av2;
898 : pari_timer ti;
899 25011 : long n = get_Flx_degree(T), N = (n+4)/2, dj;
900 25011 : GEN q = powuu(p, N);
901 : GEN T2, Xm, s1, c2, t, lr;
902 : GEN S1, sqx;
903 : GEN Nc2, Np;
904 25011 : GEN act, phi = get_Kohel_polynomials(p, &act, &dj);
905 25011 : long ispcyc = zx_is_pcyc(get_Flx_mod(T));
906 25011 : timer_start(&ti);
907 25011 : if (!ispcyc)
908 : {
909 13916 : T2 = Flx_Teichmuller(get_Flx_mod(T),p,N);
910 13916 : if (DEBUGLEVEL) timer_printf(&ti,"Teich");
911 : } else
912 11095 : T2 = Flx_to_ZX(get_Flx_mod(T));
913 25011 : T2 = FpX_get_red(T2, q); T = ZXT_to_FlxT(T2, p);
914 25011 : av2 = avma;
915 25011 : if (DEBUGLEVEL) timer_printf(&ti,"Barrett");
916 25011 : if (!ispcyc)
917 : {
918 13916 : Xm = FpXQ_powers(pol_xn(n,get_FpX_var(T2)),p-1,T2,q);
919 13916 : if (DEBUGLEVEL) timer_printf(&ti,"Xm");
920 : } else
921 11095 : Xm = cgetg(1,t_VEC);
922 25011 : s1 = Flxq_inv(Flx_Fl_add(Flxq_ellj(a4,a6,T,p),dj, p),T,p);
923 25011 : lr = Flxq_lroot(polx_Flx(get_Flx_var(T)), T, p);
924 25011 : sqx = Flxq_powers(lr, p-1, T, p);
925 25011 : S1 = lift_isogeny(phi, Flx_to_ZX(s1), N, Xm, T2, sqx, T ,p);
926 25011 : if (DEBUGLEVEL) timer_printf(&ti,"Lift isogeny");
927 25011 : c2 = getc2(act, S1, T2, q, p, N);
928 25011 : if (DEBUGLEVEL) timer_printf(&ti,"c^2");
929 25011 : if (p>3 && !ispcyc)
930 : {
931 10199 : GEN c2p = Flx_to_ZX(Flxq_inv(ZX_to_Flx(c2,p),T,p));
932 10199 : GEN tc2 = Teichmuller_lift(c2p,Xm, T2,sqx,T,p,N);
933 10199 : if (DEBUGLEVEL) timer_printf(&ti,"Teichmuller/Fq");
934 10199 : c2 = FpX_rem(FpX_mul(tc2,c2,q),T2,q);
935 : }
936 25011 : c2 = gerepileupto(av2, c2);
937 25011 : if (DEBUGLEVEL) timer_printf(&ti,"tc2");
938 25011 : Nc2 = (ispcyc? ZpXQ_sqrtnorm_pcyc: ZpXQ_sqrtnorm)(c2, T2, q, utoipos(p), N);
939 25011 : if (DEBUGLEVEL) timer_printf(&ti,"Norm");
940 25011 : Np = get_norm(a4,a6,T,p,N);
941 25011 : if (p>3 && ispcyc)
942 : {
943 7 : GEN Ncpi = utoi(Fl_inv(umodiu(Nc2,p), p));
944 7 : GEN tNc2 = Zp_sqrtnlift(gen_1, subss(p,1), Ncpi, utoipos(p),N);
945 7 : if (DEBUGLEVEL) timer_printf(&ti,"Teichmuller/Fp");
946 7 : Nc2 = Fp_mul(Nc2,tNc2,q);
947 : }
948 25011 : t = Fp_center_i(Fp_mul(Nc2,Np,q),q,shifti(q,-1));
949 25011 : return gerepileupto(av, subii(addiu(powuu(p,n),1),t));
950 : }
951 :
952 : /* Use Damien Robert method */
953 :
954 : static GEN
955 21 : get_trace_Robert(GEN J, GEN phi, GEN Xm, GEN T, GEN q, ulong p, long e)
956 : {
957 21 : long n = lg(phi)-2;
958 21 : GEN K = ZpXQ_frob(J, Xm, T, q, p);
959 21 : GEN Jp = FpXQ_powers(J, n, T, q);
960 21 : GEN Kp = FpXQ_powers(K, n, T, q);
961 21 : GEN Jd = FpXC_powderiv(Jp, q);
962 21 : GEN Kd = FpXC_powderiv(Kp, q);
963 21 : GEN Dx = FpM_FpXQV_bilinear(phi, Kd, Jp, T, q);
964 21 : GEN Dy = FpM_FpXQV_bilinear(phi, Kp, Jd, T, q);
965 21 : GEN C = ZpXQ_inv(ZX_divuexact(Dy, p), T, utoi(p), e);
966 21 : return FpX_neg(FpXQ_mul(Dx, C, T, q), q);
967 : }
968 :
969 : static GEN
970 21 : Flxq_ellcard_Harley(GEN a4, GEN a6, GEN T, ulong p)
971 : {
972 21 : pari_sp av = avma, av2;
973 : pari_timer ti;
974 21 : long n = get_Flx_degree(T), N = (n+5)/2;
975 21 : GEN pp = utoipos(p), q = powuu(p, N);
976 : GEN T2, j, t, phi;
977 : GEN J1,sqx,Xm;
978 : GEN c2, tc2, c2p, Nc2, Np;
979 21 : long ispcyc = zx_is_pcyc(get_Flx_mod(T));
980 21 : timer_start(&ti);
981 21 : if (!ispcyc)
982 : {
983 14 : T2 = Flx_Teichmuller(get_Flx_mod(T),p,N);
984 14 : if (DEBUGLEVEL) timer_printf(&ti,"Teich");
985 : } else
986 7 : T2 = Flx_to_ZX(get_Flx_mod(T));
987 21 : T2 = FpX_get_red(T2, q); T = ZXT_to_FlxT(T2, p);
988 21 : av2 = avma;
989 21 : if (DEBUGLEVEL) timer_printf(&ti,"Barrett");
990 21 : if (!ispcyc)
991 : {
992 14 : Xm = FpXQ_powers(pol_xn(n,get_FpX_var(T2)),p-1,T2,q);
993 14 : if (DEBUGLEVEL) timer_printf(&ti,"Xm");
994 : } else
995 7 : Xm = cgetg(1,t_VEC);
996 21 : j = Flxq_ellj(a4,a6,T,p);
997 21 : sqx = Flxq_powers(Flxq_lroot(polx_Flx(T[1]), T, p), p-1, T, p);
998 21 : phi = polmodular_ZM(p, 0);
999 21 : if (DEBUGLEVEL) timer_printf(&ti,"phi");
1000 21 : J1 = lift_isogeny(phi, Flx_to_ZX(j), N, Xm, T2,sqx,T,p);
1001 21 : if (DEBUGLEVEL) timer_printf(&ti,"Lift isogeny");
1002 21 : c2 = get_trace_Robert(J1, phi, Xm, T2, q, p, N);
1003 21 : q = diviuexact(q,p); N--;
1004 21 : if (DEBUGLEVEL) timer_printf(&ti,"c^2");
1005 21 : if (!ispcyc)
1006 : {
1007 14 : c2p = Flx_to_ZX(Flxq_inv(ZX_to_Flx(c2,p),T,p));
1008 14 : tc2 = Teichmuller_lift(c2p,Xm, T2,sqx,T,p,N);
1009 14 : if (DEBUGLEVEL) timer_printf(&ti,"teichmuller");
1010 14 : c2 = FpX_rem(FpX_mul(tc2,c2,q),T2,q);
1011 : }
1012 21 : c2 = gerepileupto(av2, c2);
1013 21 : q = powuu(p, N);
1014 21 : Nc2 = (ispcyc? ZpXQ_sqrtnorm_pcyc: ZpXQ_sqrtnorm)(c2, T2, q, pp, N);
1015 21 : if (DEBUGLEVEL) timer_printf(&ti,"Norm");
1016 21 : Np = get_norm(a4,a6,T,p,N);
1017 21 : if (ispcyc)
1018 : {
1019 7 : GEN Ncpi = utoi(Fl_inv(umodiu(Nc2,p), p));
1020 7 : GEN tNc2 = Zp_sqrtnlift(gen_1, subss(p,1), Ncpi, pp, N);
1021 7 : if (DEBUGLEVEL) timer_printf(&ti,"Teichmuller/Fp");
1022 7 : Nc2 = Fp_mul(Nc2,tNc2,q);
1023 : }
1024 21 : t = Fp_center_i(Fp_mul(Nc2,Np,q),q,shifti(q,-1));
1025 21 : return gerepileupto(av, subii(addiu(powuu(p,n),1),t));
1026 : }
1027 :
1028 : /***************************************************************************/
1029 : /* */
1030 : /* Shanks Mestre */
1031 : /* */
1032 : /***************************************************************************/
1033 :
1034 : /* Return the lift of a (mod b), which is closest to h */
1035 : static GEN
1036 2309 : closest_lift(GEN a, GEN b, GEN h)
1037 : {
1038 2309 : return addii(a, mulii(b, diviiround(subii(h,a), b)));
1039 : }
1040 :
1041 : static GEN
1042 1273 : FlxqE_find_order(GEN f, GEN h, GEN bound, GEN B, GEN a4, GEN T, ulong p)
1043 : {
1044 1273 : pari_sp av = avma, av1;
1045 : pari_timer Ti;
1046 1273 : long s = itos( gceil(gsqrt(gdiv(bound,B),DEFAULTPREC)) ) >> 1;
1047 : GEN tx, ti;
1048 1273 : GEN fh = FlxqE_mul(f, h, a4, T, p);
1049 1273 : GEN F, P = fh, fg;
1050 : long i;
1051 1273 : if (DEBUGLEVEL >= 6) timer_start(&Ti);
1052 1273 : if (ell_is_inf(fh)) return h;
1053 1149 : F = FlxqE_mul(f, B, a4, T, p);
1054 1149 : if (s < 3)
1055 : { /* we're nearly done: naive search */
1056 246 : GEN Q = P;
1057 246 : for (i=1;; i++)
1058 : {
1059 690 : P = FlxqE_add(P, F, a4, T, p); /* h.f + i.F */
1060 690 : if (ell_is_inf(P)) return gerepileupto(av, addii(h, mului(i,B)));
1061 620 : Q = FlxqE_sub(Q, F, a4, T, p); /* h.f - i.F */
1062 620 : if (ell_is_inf(Q)) return gerepileupto(av, subii(h, mului(i,B)));
1063 : }
1064 : }
1065 903 : tx = cgetg(s+1,t_VECSMALL);
1066 : /* Baby Step/Giant Step */
1067 903 : av1 = avma;
1068 5224 : for (i=1; i<=s; i++)
1069 : { /* baby steps */
1070 4453 : tx[i] = hash_GEN(gel(P, 1));
1071 4453 : P = FlxqE_add(P, F, a4, T, p); /* h.f + i.F */
1072 4453 : if (ell_is_inf(P)) return gerepileupto(av, addii(h, mului(i,B)));
1073 4321 : if (gc_needed(av1,3))
1074 : {
1075 0 : if(DEBUGMEM>1) pari_warn(warnmem,"[Flxq_ellcard] baby steps, i=%ld",i);
1076 0 : P = gerepileupto(av1,P);
1077 : }
1078 : }
1079 771 : if (DEBUGLEVEL >= 6) timer_printf(&Ti, "[Flxq_ellcard] baby steps, s = %ld",s);
1080 : /* giant steps: fg = s.F */
1081 771 : fg = gerepileupto(av1, FlxqE_sub(P, fh, a4, T, p));
1082 771 : if (ell_is_inf(fg)) return gerepileupto(av,mului(s,B));
1083 771 : ti = vecsmall_indexsort(tx); /* = permutation sorting tx */
1084 771 : tx = perm_mul(tx,ti);
1085 771 : if (DEBUGLEVEL >= 6) timer_printf(&Ti, "[Flxq_ellcard] sorting");
1086 771 : av1 = avma;
1087 771 : for (P=fg, i=1; ; i++)
1088 2917 : {
1089 3688 : long k = hash_GEN(gel(P,1));
1090 3688 : long r = zv_search(tx, k);
1091 3688 : if (r)
1092 : {
1093 1543 : while (r && tx[r] == k) r--;
1094 771 : for (r++; r <= s && tx[r] == k; r++)
1095 : {
1096 771 : long j = ti[r]-1;
1097 771 : GEN Q = FlxqE_add(FlxqE_mul(F, stoi(j), a4, T, p), fh, a4, T, p);
1098 771 : if (DEBUGLEVEL >= 6)
1099 0 : timer_printf(&Ti, "[Flxq_ellcard] giant steps, i = %ld",i);
1100 771 : if (Flx_equal(gel(P,1), gel(Q,1)))
1101 : {
1102 771 : if (Flx_equal(gel(P,2), gel(Q,2))) i = -i;
1103 771 : return gerepileupto(av,addii(h, mulii(addis(mulss(s,i), j), B)));
1104 : }
1105 : }
1106 : }
1107 2917 : P = FlxqE_add(P,fg,a4,T,p);
1108 2917 : if (gc_needed(av1,3))
1109 : {
1110 0 : if(DEBUGMEM>1) pari_warn(warnmem,"[Flxq_ellcard] giants steps, i=%ld",i);
1111 0 : P = gerepileupto(av1,P);
1112 : }
1113 : }
1114 : }
1115 :
1116 : static void
1117 33614 : Flx_next(GEN t, ulong p)
1118 : {
1119 : long i;
1120 33614 : for(i=2;;i++)
1121 41573 : if (uel(t,i)==p-1)
1122 7959 : t[i]=0;
1123 : else
1124 : {
1125 33614 : t[i]++;
1126 33614 : break;
1127 : }
1128 33614 : }
1129 :
1130 : static void
1131 33614 : Flx_renormalize_ip(GEN x, long lx)
1132 : {
1133 : long i;
1134 41573 : for (i = lx-1; i>=2; i--)
1135 38605 : if (x[i]) break;
1136 33614 : setlg(x, i+1);
1137 33614 : }
1138 :
1139 : static ulong
1140 2282 : F3xq_ellcard_naive(GEN a2, GEN a6, GEN T)
1141 : {
1142 2282 : pari_sp av = avma;
1143 2282 : long i, d = get_Flx_degree(T), lx = d+2;
1144 2282 : long q = upowuu(3, d), a;
1145 2282 : GEN x = zero_zv(lx); x[1] = get_Flx_var(T);
1146 11354 : for(a=1, i=0; i<q; i++)
1147 : {
1148 : GEN rhs;
1149 9072 : Flx_renormalize_ip(x, lx);
1150 9072 : rhs = Flx_add(Flxq_mul(Flxq_sqr(x, T, 3), Flx_add(x, a2, 3), T, 3), a6, 3);
1151 9072 : if (!lgpol(rhs)) a++; else if (Flxq_issquare(rhs, T, 3)) a+=2;
1152 9072 : Flx_next(x, 3);
1153 : }
1154 2282 : set_avma(av);
1155 2282 : return a;
1156 : }
1157 :
1158 : static ulong
1159 686 : Flxq_ellcard_naive(GEN a4, GEN a6, GEN T, ulong p)
1160 : {
1161 686 : pari_sp av = avma;
1162 686 : long i, d = get_Flx_degree(T), lx = d+2;
1163 686 : long q = upowuu(p, d), a;
1164 686 : GEN x = zero_zv(lx); x[1] = get_Flx_var(T);
1165 25228 : for(a=1, i=0; i<q; i++)
1166 : {
1167 : GEN x2, rhs;
1168 24542 : Flx_renormalize_ip(x, lx);
1169 24542 : x2 = Flxq_sqr(x, T, p);
1170 24542 : rhs = Flx_add(Flxq_mul(x, Flx_add(x2, a4, p), T, p), a6, p);
1171 24542 : if (!lgpol(rhs)) a++; else if (Flxq_issquare(rhs,T,p)) a+=2;
1172 24542 : Flx_next(x,p);
1173 : }
1174 686 : set_avma(av);
1175 686 : return a;
1176 : }
1177 :
1178 : /* assume T irreducible mod p, m = (q-1)/(p-1) */
1179 : static long
1180 2611 : Flxq_kronecker(GEN x, GEN m, GEN T, ulong p)
1181 : {
1182 : pari_sp av;
1183 : ulong z;
1184 2611 : if (lgpol(x) == 0) return 0;
1185 2598 : av = avma; z = Flxq_pow(x, m, T, p)[2];
1186 2598 : return gc_long(av, krouu(z, p));
1187 : }
1188 :
1189 : /* Find x such that kronecker(u = x^3+a4x+a6, p) is KRO.
1190 : * Return point [x*u,u^2] on E (KRO=1) / E^twist (KRO=-1) */
1191 : static GEN
1192 2611 : Flxq_ellpoint(long KRO, GEN a4, GEN a6, GEN m, long n, long vn, GEN T, ulong p)
1193 : {
1194 : for(;;)
1195 1338 : {
1196 2611 : GEN x = random_Flx(n,vn,p);
1197 2611 : GEN u = Flx_add(a6, Flxq_mul(Flx_add(a4, Flxq_sqr(x,T,p), p), x, T,p), p);
1198 2611 : if (Flxq_kronecker(u, m,T,p) == KRO)
1199 1273 : return mkvec2(Flxq_mul(u,x, T,p), Flxq_sqr(u, T,p));
1200 : }
1201 : }
1202 :
1203 : static GEN
1204 1036 : Flxq_ellcard_Shanks(GEN a4, GEN a6, GEN q, GEN T, ulong p)
1205 : {
1206 1036 : pari_sp av = avma;
1207 1036 : long vn = get_Flx_var(T), n = get_Flx_degree(T), KRO = -1;
1208 : GEN h,f, ta4, A, B, m;
1209 1036 : GEN q1p = addiu(q,1), q2p = shifti(q1p, 1);
1210 1036 : GEN bound = addiu(sqrti(gmul2n(q,4)), 1); /* ceil( 4sqrt(q) ) */
1211 : /* once #E(Flxq) is know mod B >= bound, it is completely determined */
1212 : /* how many 2-torsion points ? */
1213 1036 : switch(FlxqX_nbroots(mkpoln(4, pol1_Flx(vn), pol0_Flx(vn), a4, a6), T, p))
1214 : {
1215 385 : case 3: A = gen_0; B = utoipos(4); break;
1216 245 : case 1: A = gen_0; B = gen_2; break;
1217 406 : default: A = gen_1; B = gen_2; break; /* 0 */
1218 : }
1219 1036 : m = diviuexact(subiu(powuu(p,n), 1), p-1);
1220 : for(;;)
1221 : {
1222 1273 : h = closest_lift(A, B, q1p);
1223 : /* [ux, u^2] is on E_u: y^2 = x^3 + c4 u^2 x + c6 u^3
1224 : * E_u isomorphic to E (resp. E') iff KRO = 1 (resp. -1)
1225 : * #E(F_p) = p+1 - a_p, #E'(F_p) = p+1 + a_p
1226 : *
1227 : * #E_u(Flxq) = A (mod B), h is close to #E_u(Flxq) */
1228 1273 : KRO = -KRO;
1229 1273 : f = Flxq_ellpoint(KRO, a4,a6, m,n,vn, T,p);
1230 :
1231 1273 : ta4 = Flxq_mul(a4, gel(f,2), T, p); /* a4 for E_u */
1232 1273 : h = FlxqE_find_order(f, h, bound, B, ta4,T,p);
1233 1273 : h = FlxqE_order(f, h, ta4, T, p);
1234 : /* h | #E_u(Flxq) = A (mod B) */
1235 1273 : A = Z_chinese_all(A, gen_0, B, h, &B);
1236 1273 : if (cmpii(B, bound) >= 0) break;
1237 : /* not done, update A mod B for the _next_ curve, isomorphic to
1238 : * the quadratic twist of this one */
1239 237 : A = remii(subii(q2p,A), B); /* #E(Fq)+#E'(Fq) = 2q+2 */
1240 : }
1241 1036 : h = closest_lift(A, B, q1p);
1242 1036 : return gerepileuptoint(av, KRO == 1? h: subii(q2p,h));
1243 : }
1244 :
1245 : static GEN
1246 17087 : F3xq_ellcard(GEN a2, GEN a6, GEN T)
1247 : {
1248 17087 : long n = get_Flx_degree(T);
1249 17087 : if (n <= 2)
1250 1981 : return utoi(F3xq_ellcard_naive(a2, a6, T));
1251 : else
1252 : {
1253 15106 : GEN q1 = addiu(powuu(3, get_Flx_degree(T)), 1), t;
1254 15106 : GEN a = Flxq_div(a6,Flxq_powu(a2,3,T,3),T,3);
1255 15106 : if (Flx_equal1(Flxq_powu(a, 8, T, 3)))
1256 : {
1257 301 : GEN P = Flxq_minpoly(a,T,3);
1258 301 : long dP = degpol(P); /* dP <= 2 */
1259 301 : ulong q = upowuu(3,dP);
1260 301 : GEN A2 = pol1_Flx(P[1]), A6 = Flx_rem(polx_Flx(P[1]), P, 3);
1261 301 : long tP = q + 1 - F3xq_ellcard_naive(A2, A6, P);
1262 301 : t = elltrace_extension(stoi(tP), n/dP, utoi(q));
1263 301 : if (umodiu(t, 3)!=1) t = negi(t);
1264 301 : return Flx_equal1(a2) || Flxq_issquare(a2,T,3) ? subii(q1,t): addii(q1,t);
1265 : }
1266 14805 : else return Flxq_ellcard_Kohel(mkvec(a2), a6, T, 3);
1267 : }
1268 : }
1269 :
1270 : static GEN
1271 10913 : Flxq_ellcard_Satoh(GEN a4, GEN a6, GEN j, GEN T, ulong p)
1272 : {
1273 10913 : long n = get_Flx_degree(T);
1274 10913 : if (n <= 2)
1275 406 : return utoi(Flxq_ellcard_naive(a4, a6, T, p));
1276 : else
1277 : {
1278 10507 : GEN jp = Flxq_powu(j, p, T, p);
1279 10507 : GEN s = Flx_add(j, jp, p);
1280 10507 : if (degpol(s) <= 0)
1281 : { /* it is assumed j not in F_p */
1282 280 : GEN m = Flxq_mul(j, jp, T, p);
1283 280 : if (degpol(m) <= 0)
1284 : {
1285 280 : GEN q = sqru(p);
1286 280 : GEN q1 = addiu(powuu(p, get_Flx_degree(T)), 1);
1287 280 : GEN sk = Flx_Fl_add(Flx_neg(j, p), 1728%p, p);
1288 280 : GEN sA4 = Flx_triple(Flxq_mul(sk, j, T, p), p);
1289 280 : GEN u = Flxq_div(a4, sA4, T, p);
1290 280 : ulong ns = lgpol(s) ? Fl_neg(s[2], p): 0UL;
1291 280 : GEN P = mkvecsmall4(T[1], m[2], ns, 1L);
1292 : GEN A4, A6, t, tP;
1293 280 : Flxq_ellj_to_a4a6(polx_Flx(T[1]), P, p, &A4, &A6);
1294 280 : tP = addis(q, 1 - Flxq_ellcard_naive(A4, A6, P, p));
1295 280 : t = elltrace_extension(tP, n>>1, q);
1296 280 : return Flxq_is2npower(u, 2, T, p) ? subii(q1,t): addii(q1,t);
1297 : }
1298 : }
1299 10227 : if (p<=7 || p==13 ) return Flxq_ellcard_Kohel(a4, a6, T, p);
1300 21 : else return Flxq_ellcard_Harley(a4, a6, T, p);
1301 : }
1302 : }
1303 :
1304 : static GEN
1305 0 : Flxq_ellcard_Kedlaya(GEN a4, GEN a6, GEN T, ulong p)
1306 : {
1307 0 : pari_sp av = avma;
1308 0 : GEN H = mkpoln(4, gen_1, gen_0, Flx_to_ZX(a4), Flx_to_ZX(a6));
1309 0 : GEN Tp = Flx_to_ZX(get_Flx_mod(T));
1310 0 : long n = degpol(Tp), e = ((p < 16 ? n+1: n)>>1)+1;
1311 0 : GEN M = ZlXQX_hyperellpadicfrobenius(H, Tp, p, e);
1312 0 : GEN N = ZpXQM_prodFrobenius(M, Tp, utoipos(p), e);
1313 0 : GEN q = powuu(p, e);
1314 0 : GEN tp = Fq_add(gcoeff(N,1,1), gcoeff(N,2,2), Tp, q);
1315 0 : GEN t = Fp_center_i(typ(tp)==t_INT ? tp: leading_coeff(tp), q, shifti(q,-1));
1316 0 : return gerepileupto(av, subii(addiu(powuu(p, n), 1), t));
1317 : }
1318 :
1319 : GEN
1320 51820 : Flxq_ellj(GEN a4, GEN a6, GEN T, ulong p)
1321 : {
1322 51820 : pari_sp av=avma;
1323 51820 : if (p==3)
1324 : {
1325 : GEN J;
1326 14805 : if (typ(a4)!=t_VEC) return pol0_Flx(get_Flx_var(T));
1327 14805 : J = Flxq_div(Flxq_powu(gel(a4,1),3, T, p),Flx_neg(a6,p), T, p);
1328 14805 : return gerepileuptoleaf(av, J);
1329 : }
1330 : else
1331 : {
1332 37015 : pari_sp av=avma;
1333 37015 : GEN a43 = Flxq_mul(a4,Flxq_sqr(a4,T,p),T,p);
1334 37015 : GEN a62 = Flxq_sqr(a6,T,p);
1335 37015 : GEN num = Flx_mulu(a43,6912,p);
1336 37015 : GEN den = Flx_add(Flx_mulu(a43,4,p),Flx_mulu(a62,27,p),p);
1337 37015 : return gerepileuptoleaf(av, Flxq_div(num, den, T, p));
1338 : }
1339 : }
1340 :
1341 : void
1342 280 : Flxq_ellj_to_a4a6(GEN j, GEN T, ulong p, GEN *pt_a4, GEN *pt_a6)
1343 : {
1344 280 : ulong zagier = 1728 % p;
1345 280 : if (lgpol(j)==0)
1346 0 : { *pt_a4 = pol0_Flx(T[1]); *pt_a6 =pol1_Flx(T[1]); }
1347 280 : else if (lgpol(j)==1 && uel(j,2) == zagier)
1348 0 : { *pt_a4 = pol1_Flx(T[1]); *pt_a6 =pol0_Flx(T[1]); }
1349 : else
1350 : {
1351 280 : GEN k = Flx_Fl_add(Flx_neg(j, p), zagier, p);
1352 280 : GEN kj = Flxq_mul(k, j, T, p);
1353 280 : GEN k2j = Flxq_mul(kj, k, T, p);
1354 280 : *pt_a4 = Flx_triple(kj, p);
1355 280 : *pt_a6 = Flx_double(k2j, p);
1356 : }
1357 280 : }
1358 :
1359 : static GEN
1360 6426 : F3xq_ellcardj(GEN a4, GEN a6, GEN T, GEN q, long n)
1361 : {
1362 6426 : const ulong p = 3;
1363 : ulong t;
1364 6426 : GEN q1 = addiu(q,1);
1365 6426 : GEN na4 = Flx_neg(a4,p), ra4;
1366 6426 : if (!Flxq_issquare(na4,T,p))
1367 3234 : return q1;
1368 3192 : ra4 = Flxq_sqrt(na4,T,p);
1369 3192 : t = Flxq_trace(Flxq_div(a6,Flxq_mul(na4,ra4,T,p),T,p),T,p);
1370 3192 : if (n%2==1)
1371 : {
1372 : GEN q3;
1373 1078 : if (t==0) return q1;
1374 301 : q3 = powuu(p,(n+1)>>1);
1375 301 : return (t==1)^(n%4==1) ? subii(q1,q3): addii(q1,q3);
1376 : }
1377 : else
1378 : {
1379 2114 : GEN q22, q2 = powuu(p,n>>1);
1380 2114 : GEN W = Flxq_pow(a4,shifti(q,-2),T,p);
1381 2114 : long s = (W[2]==1)^(n%4==2);
1382 2114 : if (t!=0) return s ? addii(q1,q2): subii(q1, q2);
1383 2114 : q22 = shifti(q2,1);
1384 2114 : return s ? subii(q1,q22): addii(q1, q22);
1385 : }
1386 : }
1387 :
1388 : static GEN
1389 14812 : Flxq_ellcardj(GEN a4, GEN a6, ulong j, GEN T, GEN q, ulong p, long n)
1390 : {
1391 14812 : GEN q1 = addiu(q,1);
1392 14812 : if (j==0)
1393 : {
1394 : ulong w;
1395 : GEN W, t, N;
1396 5600 : if (umodiu(q,6)!=1) return q1;
1397 4200 : N = Fp_ffellcard(gen_0,gen_1,q,n,utoipos(p));
1398 4200 : t = subii(q1, N);
1399 4200 : W = Flxq_pow(a6,diviuexact(shifti(q,-1), 3),T,p);
1400 4200 : if (degpol(W)>0) /*p=5 mod 6*/
1401 1428 : return Flx_equal1(Flxq_powu(W,3,T,p)) ? addii(q1,shifti(t,-1)):
1402 476 : subii(q1,shifti(t,-1));
1403 3248 : w = W[2];
1404 3248 : if (w==1) return N;
1405 2576 : if (w==p-1) return addii(q1,t);
1406 : else /*p=1 mod 6*/
1407 : {
1408 1904 : GEN u = shifti(t,-1), v = sqrtint(diviuexact(subii(q,sqri(u)),3));
1409 1904 : GEN a = addii(u,v), b = shifti(v,1);
1410 1904 : if (Fl_powu(w,3,p)==1)
1411 : {
1412 952 : if (Fl_add(umodiu(a,p),Fl_mul(w,umodiu(b,p),p),p)==0)
1413 511 : return subii(q1,subii(shifti(b,1),a));
1414 : else
1415 441 : return addii(q1,addii(a,b));
1416 : }
1417 : else
1418 : {
1419 952 : if (Fl_sub(umodiu(a,p),Fl_mul(w,umodiu(b,p),p),p)==0)
1420 511 : return subii(q1,subii(a,shifti(b,1)));
1421 : else
1422 441 : return subii(q1,addii(a,b));
1423 : }
1424 : }
1425 9212 : } else if (j==1728%p)
1426 : {
1427 : ulong w;
1428 : GEN W, N, t;
1429 5607 : if (mod4(q)==3) return q1;
1430 4207 : W = Flxq_pow(a4,shifti(q,-2),T,p);
1431 4207 : if (degpol(W)>0) return q1; /*p=3 mod 4*/
1432 3521 : w = W[2];
1433 3521 : N = Fp_ffellcard(gen_1,gen_0,q,n,utoipos(p));
1434 3521 : if(w==1) return N;
1435 2464 : t = subii(q1, N);
1436 2464 : if(w==p-1) return addii(q1, t);
1437 : else /*p=1 mod 4*/
1438 : {
1439 1400 : GEN u = shifti(t,-1), v = sqrtint(subii(q,sqri(u)));
1440 1400 : if (Fl_add(umodiu(u,p),Fl_mul(w,umodiu(v,p),p),p)==0)
1441 700 : return subii(q1,shifti(v,1));
1442 : else
1443 700 : return addii(q1,shifti(v,1));
1444 : }
1445 : } else
1446 : {
1447 3605 : ulong g = Fl_div(j, Fl_sub(1728%p, j, p), p);
1448 3605 : GEN l = Flxq_div(Flx_triple(a6,p),Flx_double(a4,p),T,p);
1449 3605 : GEN N = Fp_ffellcard(utoi(Fl_triple(g,p)),utoi(Fl_double(g,p)),q,n,utoipos(p));
1450 3605 : if (Flxq_issquare(l,T,p)) return N;
1451 2184 : return subii(shifti(q1,1),N);
1452 : }
1453 : }
1454 :
1455 : GEN
1456 50518 : Flxq_ellcard(GEN a4, GEN a6, GEN T, ulong p)
1457 : {
1458 50518 : pari_sp av = avma;
1459 50518 : long n = get_Flx_degree(T);
1460 50518 : GEN J, r, q = powuu(p, n);
1461 50518 : if (typ(a4)==t_VEC)
1462 17087 : r = F3xq_ellcard(gel(a4,1), a6, T);
1463 33431 : else if (p==3)
1464 6426 : r = F3xq_ellcardj(a4, a6, T, q, n);
1465 27005 : else if (degpol(a4)<=0 && degpol(a6)<=0)
1466 217 : r = Fp_ffellcard(utoi(Flx_eval(a4,0,p)),utoi(Flx_eval(a6,0,p)),q,n,utoipos(p));
1467 26788 : else if (degpol(J=Flxq_ellj(a4,a6,T,p))<=0)
1468 14812 : r = Flxq_ellcardj(a4,a6,lgpol(J)?J[2]:0,T,q,p,n);
1469 11976 : else if (p <= 7)
1470 10843 : r = Flxq_ellcard_Satoh(a4, a6, J, T, p);
1471 1133 : else if (cmpis(q,100)<0)
1472 0 : r = utoi(Flxq_ellcard_naive(a4, a6, T, p));
1473 1133 : else if (p == 13 || (7*p <= (ulong)10*n && (BITS_IN_LONG==64 || p <= 103)))
1474 70 : r = Flxq_ellcard_Satoh(a4, a6, J, T, p);
1475 1063 : else if (p <= (ulong)2*n)
1476 0 : r = Flxq_ellcard_Kedlaya(a4, a6, T, p);
1477 1063 : else if (expi(q)<=62)
1478 1036 : r = Flxq_ellcard_Shanks(a4, a6, q, T, p);
1479 : else
1480 27 : r = Fq_ellcard_SEA(Flx_to_ZX(a4),Flx_to_ZX(a6),q,Flx_to_ZX(T),utoipos(p),0);
1481 50518 : return gerepileuptoint(av, r);
1482 : }
|
