Line data Source code
1 : /* Copyright (C) 2009 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 Fp */
21 :
22 : /***********************************************************************/
23 : /** **/
24 : /** FpJ **/
25 : /** **/
26 : /***********************************************************************/
27 : /* Arithmetic is implemented using Jacobian coordinates, representing
28 : * a projective point (x : y : z) on E by [z*x , z^2*y , z]. This is
29 : * probably not the fastest representation available for the given
30 : * problem, but they're easy to implement and up to 60% faster than
31 : * the school-book method used in FpE_mulu(). */
32 :
33 : static GEN
34 51085 : ellinf_FpJ(void)
35 51085 : { return mkvec3(gen_1, gen_1, gen_0); }
36 :
37 : /* Cost: 1M + 8S + 1*a + 10add + 1*8 + 2*2 + 1*3.
38 : * Source: http://www.hyperelliptic.org/EFD/g1p/auto-shortw-jacobian.html#doubling-dbl-2007-bl */
39 : GEN
40 6354539 : FpJ_dbl(GEN P, GEN a4, GEN p)
41 : {
42 : GEN X1, Y1, Z1;
43 : GEN XX, YY, YYYY, ZZ, S, M, T, Q;
44 :
45 6354539 : if (signe(gel(P,3)) == 0) return ellinf_FpJ();
46 :
47 6345302 : X1 = gel(P,1); Y1 = gel(P,2); Z1 = gel(P,3);
48 :
49 6345302 : XX = Fp_sqr(X1, p);
50 6335360 : YY = Fp_sqr(Y1, p);
51 6335616 : YYYY = Fp_sqr(YY, p);
52 6336018 : ZZ = Fp_sqr(Z1, p);
53 6336902 : S = Fp_double(Fp_sub(Fp_sqr(Fp_add(X1,YY,p), p), Fp_add(XX,YYYY,p), p), p);
54 6256941 : M = Fp_addmul(Fp_mulu(XX, 3, p), a4, Fp_sqr(ZZ, p), p);
55 6323820 : T = Fp_sub(Fp_sqr(M, p), Fp_double(S, p), p);
56 6284178 : Q = cgetg(4, t_VEC);
57 6322584 : gel(Q,1) = T;
58 6322584 : gel(Q,2) = Fp_sub(Fp_mul(M, Fp_sub(S, T, p), p), Fp_mulu(YYYY, 8, p), p);
59 6285875 : gel(Q,3) = Fp_sub(Fp_sqr(Fp_add(Y1, Z1, p), p), Fp_add(YY, ZZ, p), p);
60 6283287 : return Q;
61 : }
62 :
63 : /* Cost: 11M + 5S + 9add + 4*2.
64 : * Source: http://www.hyperelliptic.org/EFD/g1p/auto-shortw-jacobian.html#addition-add-2007-bl */
65 : GEN
66 1149042 : FpJ_add(GEN P, GEN Q, GEN a4, GEN p)
67 : {
68 : GEN X1, Y1, Z1, X2, Y2, Z2;
69 : GEN Z1Z1, Z2Z2, U1, U2, S1, S2, H, I, J, r, V, W, R;
70 :
71 1149042 : if (signe(gel(Q,3)) == 0) return gcopy(P);
72 1149042 : if (signe(gel(P,3)) == 0) return gcopy(Q);
73 :
74 1129032 : X1 = gel(P,1); Y1 = gel(P,2); Z1 = gel(P,3);
75 1129032 : X2 = gel(Q,1); Y2 = gel(Q,2); Z2 = gel(Q,3);
76 :
77 1129032 : Z1Z1 = Fp_sqr(Z1, p);
78 1129074 : Z2Z2 = Fp_sqr(Z2, p);
79 1128924 : U1 = Fp_mul(X1, Z2Z2, p);
80 1129075 : U2 = Fp_mul(X2, Z1Z1, p);
81 1129089 : S1 = mulii(Y1, Fp_mul(Z2, Z2Z2, p));
82 1127600 : S2 = mulii(Y2, Fp_mul(Z1, Z1Z1, p));
83 1127601 : H = Fp_sub(U2, U1, p);
84 1128282 : r = Fp_double(Fp_sub(S2, S1, p), p);
85 :
86 : /* If points are equal we must double. */
87 1127642 : if (signe(H)== 0) {
88 42982 : if (signe(r) == 0)
89 : /* Points are equal so double. */
90 1134 : return FpJ_dbl(P, a4, p);
91 : else
92 41848 : return ellinf_FpJ();
93 : }
94 1084660 : I = Fp_sqr(Fp_double(H, p), p);
95 1086077 : J = Fp_mul(H, I, p);
96 1086149 : V = Fp_mul(U1, I, p);
97 1086142 : W = Fp_sub(Fp_sqr(r, p), Fp_add(J, Fp_double(V, p), p), p);
98 1085101 : R = cgetg(4, t_VEC);
99 1085827 : gel(R,1) = W;
100 1085827 : gel(R,2) = Fp_sub(mulii(r, subii(V, W)),
101 : shifti(mulii(S1, J), 1), p);
102 1086068 : gel(R,3) = Fp_mul(Fp_sub(Fp_sqr(Fp_add(Z1, Z2, p), p),
103 : Fp_add(Z1Z1, Z2Z2, p), p), H, p);
104 1086089 : return R;
105 : }
106 :
107 : GEN
108 0 : FpJ_neg(GEN Q, GEN p)
109 : {
110 0 : return mkvec3(icopy(gel(Q,1)), Fp_neg(gel(Q,2), p), icopy(gel(Q,3)));
111 : }
112 :
113 : GEN
114 206433 : FpE_to_FpJ(GEN P)
115 : {
116 206433 : return ell_is_inf(P) ? ellinf_FpJ()
117 206433 : : mkvec3(icopy(gel(P,1)),icopy(gel(P,2)), gen_1);
118 : }
119 :
120 : GEN
121 205900 : FpJ_to_FpE(GEN P, GEN p)
122 : {
123 205900 : if (signe(gel(P,3)) == 0) return ellinf();
124 : else
125 : {
126 164479 : GEN Z = Fp_inv(gel(P,3), p);
127 164453 : GEN Z2 = Fp_sqr(Z, p), Z3 = Fp_mul(Z, Z2, p);
128 164453 : retmkvec2(Fp_mul(gel(P,1), Z2, p), Fp_mul(gel(P,2), Z3, p));
129 : }
130 : }
131 :
132 : struct _FpE { GEN p,a4,a6; };
133 : static GEN
134 6334473 : _FpJ_dbl(void *E, GEN P)
135 : {
136 6334473 : struct _FpE *ell = (struct _FpE *) E;
137 6334473 : return FpJ_dbl(P, ell->a4, ell->p);
138 : }
139 : static GEN
140 1148560 : _FpJ_add(void *E, GEN P, GEN Q)
141 : {
142 1148560 : struct _FpE *ell=(struct _FpE *) E;
143 1148560 : return FpJ_add(P, Q, ell->a4, ell->p);
144 : }
145 : static GEN
146 5865 : _FpJ_mul(void *E, GEN P, GEN n)
147 : {
148 5865 : pari_sp av = avma;
149 5865 : struct _FpE *e=(struct _FpE *) E;
150 5865 : long s = signe(n);
151 5865 : if (!s || signe(gel(P,3))==0) return ellinf_FpJ();
152 5866 : if (s < 0) P = FpJ_neg(P, e->p);
153 5866 : if (is_pm1(n)) return s > 0 ? gcopy(P): P;
154 5866 : return gerepilecopy(av, gen_pow_i(P, n, e, &_FpJ_dbl, &_FpJ_add));
155 : }
156 :
157 : GEN
158 5865 : FpJ_mul(GEN P, GEN n, GEN a4, GEN p)
159 : {
160 : struct _FpE E;
161 5865 : E.a4= a4; E.p = p;
162 5865 : return _FpJ_mul(&E, P, n);
163 : }
164 :
165 : /***********************************************************************/
166 : /** **/
167 : /** FpE **/
168 : /** **/
169 : /***********************************************************************/
170 : /* These functions deal with point over elliptic curves over Fp defined
171 : * by an equation of the form y^2=x^3+a4*x+a6.
172 : * Most of the time a6 is omitted since it can be recovered from any point
173 : * on the curve. */
174 :
175 : GEN
176 2738 : RgE_to_FpE(GEN x, GEN p)
177 : {
178 2738 : if (ell_is_inf(x)) return x;
179 2738 : retmkvec2(Rg_to_Fp(gel(x,1),p),Rg_to_Fp(gel(x,2),p));
180 : }
181 :
182 : GEN
183 1058 : FpE_to_mod(GEN x, GEN p)
184 : {
185 1058 : if (ell_is_inf(x)) return x;
186 995 : retmkvec2(Fp_to_mod(gel(x,1),p),Fp_to_mod(gel(x,2),p));
187 : }
188 :
189 : GEN
190 1730 : FpE_changepoint(GEN P, GEN ch, GEN p)
191 : {
192 1730 : pari_sp av = avma;
193 : GEN c, z, u, r, s, t, v, v2, v3;
194 1730 : if (ell_is_inf(P)) return P;
195 1667 : if (lgefint(p) == 3)
196 : {
197 719 : ulong pp = p[2];
198 719 : z = Fle_changepoint(ZV_to_Flv(P, pp), ZV_to_Flv(ch, pp), pp);
199 719 : return gerepileupto(av, Flv_to_ZV(z));
200 : }
201 948 : u = gel(ch,1); r = gel(ch,2); s = gel(ch,3); t = gel(ch,4);
202 948 : v = Fp_inv(u, p); v2 = Fp_sqr(v,p); v3 = Fp_mul(v,v2,p);
203 948 : c = Fp_sub(gel(P,1),r,p);
204 947 : z = cgetg(3,t_VEC);
205 948 : gel(z,1) = Fp_mul(v2, c, p);
206 948 : gel(z,2) = Fp_mul(v3, Fp_sub(gel(P,2), Fp_add(Fp_mul(s,c, p),t, p),p),p);
207 947 : return gerepileupto(av, z);
208 : }
209 :
210 : GEN
211 2738 : FpE_changepointinv(GEN P, GEN ch, GEN p)
212 : {
213 : GEN u, r, s, t, u2, u3, c, z;
214 2738 : if (ell_is_inf(P)) return P;
215 2738 : if (lgefint(p) == 3)
216 : {
217 1738 : ulong pp = p[2];
218 1738 : z = Fle_changepointinv(ZV_to_Flv(P, pp), ZV_to_Flv(ch, pp), pp);
219 1738 : return Flv_to_ZV(z);
220 : }
221 1000 : u = gel(ch,1); r = gel(ch,2); s = gel(ch,3); t = gel(ch,4);
222 1000 : u2 = Fp_sqr(u, p); u3 = Fp_mul(u,u2,p);
223 1000 : c = Fp_mul(u2, gel(P,1), p);
224 1000 : z = cgetg(3, t_VEC);
225 1000 : gel(z,1) = Fp_add(c,r,p);
226 1000 : gel(z,2) = Fp_add(Fp_mul(u3,gel(P,2),p), Fp_add(Fp_mul(s,c,p), t, p), p);
227 1000 : return z;
228 : }
229 :
230 : static GEN
231 420 : random_nonsquare_Fp(GEN p)
232 : {
233 420 : pari_sp av = avma;
234 : GEN a;
235 420 : switch(mod8(p))
236 : { /* easy special cases */
237 420 : case 3: case 5: return gen_2;
238 0 : case 7: return subiu(p, 1);
239 : }
240 : do
241 : {
242 0 : set_avma(av);
243 0 : a = randomi(p);
244 0 : } while (kronecker(a, p) >= 0);
245 0 : return a;
246 : }
247 :
248 : void
249 0 : Fp_elltwist(GEN a4, GEN a6, GEN p, GEN *pt_a4, GEN *pt_a6)
250 : {
251 0 : GEN d = random_nonsquare_Fp(p), d2 = Fp_sqr(d, p), d3 = Fp_mul(d2, d, p);
252 0 : *pt_a4 = Fp_mul(a4, d2, p);
253 0 : *pt_a6 = Fp_mul(a6, d3, p);
254 0 : }
255 :
256 : static GEN
257 292458 : FpE_dbl_slope(GEN P, GEN a4, GEN p, GEN *slope)
258 : {
259 : GEN x, y, Q;
260 292458 : if (ell_is_inf(P) || !signe(gel(P,2))) return ellinf();
261 134711 : x = gel(P,1); y = gel(P,2);
262 134711 : *slope = Fp_div(Fp_add(Fp_mulu(Fp_sqr(x,p), 3, p), a4, p),
263 : Fp_mulu(y, 2, p), p);
264 134711 : Q = cgetg(3,t_VEC);
265 134711 : gel(Q, 1) = Fp_sub(Fp_sqr(*slope, p), Fp_mulu(x, 2, p), p);
266 134711 : gel(Q, 2) = Fp_sub(Fp_mul(*slope, Fp_sub(x, gel(Q, 1), p), p), y, p);
267 134711 : return Q;
268 : }
269 :
270 : GEN
271 291864 : FpE_dbl(GEN P, GEN a4, GEN p)
272 : {
273 291864 : pari_sp av = avma;
274 : GEN slope;
275 291864 : return gerepileupto(av, FpE_dbl_slope(P,a4,p,&slope));
276 : }
277 :
278 : static GEN
279 916539 : FpE_add_slope(GEN P, GEN Q, GEN a4, GEN p, GEN *slope)
280 : {
281 : GEN Px, Py, Qx, Qy, R;
282 916539 : if (ell_is_inf(P)) return Q;
283 916049 : if (ell_is_inf(Q)) return P;
284 916049 : Px = gel(P,1); Py = gel(P,2);
285 916049 : Qx = gel(Q,1); Qy = gel(Q,2);
286 916049 : if (equalii(Px, Qx))
287 : {
288 574 : if (equalii(Py, Qy))
289 553 : return FpE_dbl_slope(P, a4, p, slope);
290 : else
291 21 : return ellinf();
292 : }
293 915475 : *slope = Fp_div(Fp_sub(Py, Qy, p), Fp_sub(Px, Qx, p), p);
294 915475 : R = cgetg(3,t_VEC);
295 915475 : gel(R, 1) = Fp_sub(Fp_sub(Fp_sqr(*slope, p), Px, p), Qx, p);
296 915475 : gel(R, 2) = Fp_sub(Fp_mul(*slope, Fp_sub(Px, gel(R, 1), p), p), Py, p);
297 915475 : return R;
298 : }
299 :
300 : GEN
301 916535 : FpE_add(GEN P, GEN Q, GEN a4, GEN p)
302 : {
303 916535 : pari_sp av = avma;
304 : GEN slope;
305 916535 : return gerepileupto(av, FpE_add_slope(P,Q,a4,p,&slope));
306 : }
307 :
308 : static GEN
309 0 : FpE_neg_i(GEN P, GEN p)
310 : {
311 0 : if (ell_is_inf(P)) return P;
312 0 : return mkvec2(gel(P,1), Fp_neg(gel(P,2), p));
313 : }
314 :
315 : GEN
316 362489 : FpE_neg(GEN P, GEN p)
317 : {
318 362489 : if (ell_is_inf(P)) return ellinf();
319 362489 : return mkvec2(gcopy(gel(P,1)), Fp_neg(gel(P,2), p));
320 : }
321 :
322 : GEN
323 0 : FpE_sub(GEN P, GEN Q, GEN a4, GEN p)
324 : {
325 0 : pari_sp av = avma;
326 : GEN slope;
327 0 : return gerepileupto(av, FpE_add_slope(P, FpE_neg_i(Q, p), a4, p, &slope));
328 : }
329 :
330 : static GEN
331 291864 : _FpE_dbl(void *E, GEN P)
332 : {
333 291864 : struct _FpE *ell = (struct _FpE *) E;
334 291864 : return FpE_dbl(P, ell->a4, ell->p);
335 : }
336 :
337 : static GEN
338 897344 : _FpE_add(void *E, GEN P, GEN Q)
339 : {
340 897344 : struct _FpE *ell=(struct _FpE *) E;
341 897344 : return FpE_add(P, Q, ell->a4, ell->p);
342 : }
343 :
344 : static GEN
345 933510 : _FpE_mul(void *E, GEN P, GEN n)
346 : {
347 933510 : pari_sp av = avma;
348 933510 : struct _FpE *e=(struct _FpE *) E;
349 933510 : long s = signe(n);
350 : GEN Q;
351 933510 : if (!s || ell_is_inf(P)) return ellinf();
352 933503 : if (s<0) P = FpE_neg(P, e->p);
353 933503 : if (is_pm1(n)) return s>0? gcopy(P): P;
354 497794 : if (equalis(n,2)) return _FpE_dbl(E, P);
355 205930 : Q = gen_pow_i(FpE_to_FpJ(P), n, e, &_FpJ_dbl, &_FpJ_add);
356 205900 : return gerepileupto(av, FpJ_to_FpE(Q, e->p));
357 : }
358 :
359 : GEN
360 1319 : FpE_mul(GEN P, GEN n, GEN a4, GEN p)
361 : {
362 : struct _FpE E;
363 1319 : E.a4 = a4; E.p = p;
364 1319 : return _FpE_mul(&E, P, n);
365 : }
366 :
367 : /* Finds a random nonsingular point on E */
368 :
369 : GEN
370 190660 : random_FpE(GEN a4, GEN a6, GEN p)
371 : {
372 190660 : pari_sp ltop = avma;
373 : GEN x, x2, y, rhs;
374 : do
375 : {
376 333429 : set_avma(ltop);
377 333429 : x = randomi(p); /* x^3+a4*x+a6 = x*(x^2+a4)+a6 */
378 333429 : x2 = Fp_sqr(x, p);
379 333429 : rhs = Fp_add(Fp_mul(x, Fp_add(x2, a4, p), p), a6, p);
380 35654 : } while ((!signe(rhs) && !signe(Fp_add(Fp_mulu(x2,3,p),a4,p)))
381 369082 : || kronecker(rhs, p) < 0);
382 190660 : y = Fp_sqrt(rhs, p);
383 190660 : if (!y) pari_err_PRIME("random_FpE", p);
384 190660 : return gerepilecopy(ltop, mkvec2(x, y));
385 : }
386 :
387 : static GEN
388 188203 : _FpE_rand(void *E)
389 : {
390 188203 : struct _FpE *e=(struct _FpE *) E;
391 188203 : return random_FpE(e->a4, e->a6, e->p);
392 : }
393 :
394 : static const struct bb_group FpE_group={_FpE_add,_FpE_mul,_FpE_rand,hash_GEN,ZV_equal,ell_is_inf,NULL};
395 :
396 : const struct bb_group *
397 903 : get_FpE_group(void ** pt_E, GEN a4, GEN a6, GEN p)
398 : {
399 903 : struct _FpE *e = (struct _FpE *) stack_malloc(sizeof(struct _FpE));
400 903 : e->a4 = a4; e->a6 = a6; e->p = p;
401 903 : *pt_E = (void *) e;
402 903 : return &FpE_group;
403 : }
404 :
405 : GEN
406 736 : FpE_order(GEN z, GEN o, GEN a4, GEN p)
407 : {
408 736 : pari_sp av = avma;
409 : struct _FpE e;
410 : GEN r;
411 736 : if (lgefint(p) == 3)
412 : {
413 630 : ulong pp = p[2];
414 630 : r = Fle_order(ZV_to_Flv(z, pp), o, umodiu(a4,pp), pp);
415 : }
416 : else
417 : {
418 106 : e.a4 = a4;
419 106 : e.p = p;
420 106 : r = gen_order(z, o, (void*)&e, &FpE_group);
421 : }
422 736 : return gerepileuptoint(av, r);
423 : }
424 :
425 : GEN
426 49 : FpE_log(GEN a, GEN b, GEN o, GEN a4, GEN p)
427 : {
428 49 : pari_sp av = avma;
429 : struct _FpE e;
430 : GEN r;
431 49 : if (lgefint(p) == 3)
432 : {
433 49 : ulong pp = p[2];
434 49 : r = Fle_log(ZV_to_Flv(a,pp), ZV_to_Flv(b,pp), o, umodiu(a4,pp), pp);
435 : }
436 : else
437 : {
438 0 : e.a4 = a4;
439 0 : e.p = p;
440 0 : r = gen_PH_log(a, b, o, (void*)&e, &FpE_group);
441 : }
442 49 : return gerepileuptoint(av, r);
443 : }
444 :
445 : /***********************************************************************/
446 : /** **/
447 : /** Pairings **/
448 : /** **/
449 : /***********************************************************************/
450 :
451 : /* Derived from APIP from and by Jerome Milan, 2012 */
452 :
453 : static GEN
454 140 : FpE_vert(GEN P, GEN Q, GEN a4, GEN p)
455 : {
456 140 : if (ell_is_inf(P))
457 51 : return gen_1;
458 89 : if (!equalii(gel(Q, 1), gel(P, 1)))
459 87 : return Fp_sub(gel(Q, 1), gel(P, 1), p);
460 2 : if (signe(gel(P,2))!=0) return gen_1;
461 2 : return Fp_inv(Fp_add(Fp_mulu(Fp_sqr(gel(P,1),p), 3, p), a4, p), p);
462 : }
463 :
464 : static GEN
465 45 : FpE_Miller_line(GEN R, GEN Q, GEN slope, GEN a4, GEN p)
466 : {
467 45 : GEN x = gel(Q, 1), y = gel(Q, 2);
468 45 : GEN tmp1 = Fp_sub(x, gel(R, 1), p);
469 45 : GEN tmp2 = Fp_add(Fp_mul(tmp1, slope, p), gel(R,2), p);
470 45 : if (!equalii(y, tmp2))
471 44 : return Fp_sub(y, tmp2, p);
472 1 : if (signe(y) == 0)
473 1 : return gen_1;
474 : else
475 : {
476 : GEN s1, s2;
477 0 : GEN y2i = Fp_inv(Fp_mulu(y, 2, p), p);
478 0 : s1 = Fp_mul(Fp_add(Fp_mulu(Fp_sqr(x, p), 3, p), a4, p), y2i, p);
479 0 : if (!equalii(s1, slope))
480 0 : return Fp_sub(s1, slope, p);
481 0 : s2 = Fp_mul(Fp_sub(Fp_mulu(x, 3, p), Fp_sqr(s1, p), p), y2i, p);
482 0 : return signe(s2)!=0 ? s2: y2i;
483 : }
484 : }
485 :
486 : /* Computes the equation of the line tangent to R and returns its
487 : evaluation at the point Q. Also doubles the point R.
488 : */
489 :
490 : static GEN
491 92 : FpE_tangent_update(GEN R, GEN Q, GEN a4, GEN p, GEN *pt_R)
492 : {
493 92 : if (ell_is_inf(R))
494 : {
495 7 : *pt_R = ellinf();
496 7 : return gen_1;
497 : }
498 85 : else if (signe(gel(R,2)) == 0)
499 : {
500 44 : *pt_R = ellinf();
501 44 : return FpE_vert(R, Q, a4, p);
502 : } else {
503 : GEN slope;
504 41 : *pt_R = FpE_dbl_slope(R, a4, p, &slope);
505 41 : return FpE_Miller_line(R, Q, slope, a4, p);
506 : }
507 : }
508 :
509 : /* Computes the equation of the line through R and P, and returns its
510 : evaluation at the point Q. Also adds P to the point R.
511 : */
512 :
513 : static GEN
514 4 : FpE_chord_update(GEN R, GEN P, GEN Q, GEN a4, GEN p, GEN *pt_R)
515 : {
516 4 : if (ell_is_inf(R))
517 : {
518 0 : *pt_R = gcopy(P);
519 0 : return FpE_vert(P, Q, a4, p);
520 : }
521 4 : else if (ell_is_inf(P))
522 : {
523 0 : *pt_R = gcopy(R);
524 0 : return FpE_vert(R, Q, a4, p);
525 : }
526 4 : else if (equalii(gel(P, 1), gel(R, 1)))
527 : {
528 0 : if (equalii(gel(P, 2), gel(R, 2)))
529 0 : return FpE_tangent_update(R, Q, a4, p, pt_R);
530 : else {
531 0 : *pt_R = ellinf();
532 0 : return FpE_vert(R, Q, a4, p);
533 : }
534 : } else {
535 : GEN slope;
536 4 : *pt_R = FpE_add_slope(P, R, a4, p, &slope);
537 4 : return FpE_Miller_line(R, Q, slope, a4, p);
538 : }
539 : }
540 :
541 : struct _FpE_miller { GEN p, a4, P; };
542 : static GEN
543 92 : FpE_Miller_dbl(void* E, GEN d)
544 : {
545 92 : struct _FpE_miller *m = (struct _FpE_miller *)E;
546 92 : GEN p = m->p, a4 = m->a4, P = m->P;
547 : GEN v, line;
548 92 : GEN N = Fp_sqr(gel(d,1), p);
549 92 : GEN D = Fp_sqr(gel(d,2), p);
550 92 : GEN point = gel(d,3);
551 92 : line = FpE_tangent_update(point, P, a4, p, &point);
552 92 : N = Fp_mul(N, line, p);
553 92 : v = FpE_vert(point, P, a4, p);
554 92 : D = Fp_mul(D, v, p); return mkvec3(N, D, point);
555 : }
556 : static GEN
557 4 : FpE_Miller_add(void* E, GEN va, GEN vb)
558 : {
559 4 : struct _FpE_miller *m = (struct _FpE_miller *)E;
560 4 : GEN p = m->p, a4= m->a4, P = m->P;
561 : GEN v, line, point;
562 4 : GEN na = gel(va,1), da = gel(va,2), pa = gel(va,3);
563 4 : GEN nb = gel(vb,1), db = gel(vb,2), pb = gel(vb,3);
564 4 : GEN N = Fp_mul(na, nb, p);
565 4 : GEN D = Fp_mul(da, db, p);
566 4 : line = FpE_chord_update(pa, pb, P, a4, p, &point);
567 4 : N = Fp_mul(N, line, p);
568 4 : v = FpE_vert(point, P, a4, p);
569 4 : D = Fp_mul(D, v, p); return mkvec3(N, D, point);
570 : }
571 :
572 : /* Returns the Miller function f_{m, Q} evaluated at the point P using
573 : * the standard Miller algorithm. */
574 : static GEN
575 44 : FpE_Miller(GEN Q, GEN P, GEN m, GEN a4, GEN p)
576 : {
577 44 : pari_sp av = avma;
578 : struct _FpE_miller d;
579 : GEN v, N, D;
580 :
581 44 : d.a4 = a4; d.p = p; d.P = P;
582 44 : v = gen_pow_i(mkvec3(gen_1,gen_1,Q), m, (void*)&d,
583 : FpE_Miller_dbl, FpE_Miller_add);
584 44 : N = gel(v,1); D = gel(v,2);
585 44 : return gerepileuptoint(av, Fp_div(N, D, p));
586 : }
587 :
588 : GEN
589 73914 : FpE_weilpairing(GEN P, GEN Q, GEN m, GEN a4, GEN p)
590 : {
591 73914 : pari_sp av = avma;
592 : GEN N, D, w;
593 73914 : if (ell_is_inf(P) || ell_is_inf(Q) || ZV_equal(P,Q)) return gen_1;
594 49404 : if (lgefint(p)==3 && lgefint(m)==3)
595 : {
596 49382 : ulong pp = p[2];
597 49382 : GEN Pp = ZV_to_Flv(P, pp), Qp = ZV_to_Flv(Q, pp);
598 49382 : ulong w = Fle_weilpairing(Pp, Qp, itou(m), umodiu(a4, pp), pp);
599 49382 : return gc_utoi(av, w);
600 : }
601 22 : N = FpE_Miller(P, Q, m, a4, p);
602 22 : D = FpE_Miller(Q, P, m, a4, p);
603 22 : w = Fp_div(N, D, p);
604 22 : if (mpodd(m)) w = Fp_neg(w, p);
605 22 : return gerepileuptoint(av, w);
606 : }
607 :
608 : GEN
609 203 : FpE_tatepairing(GEN P, GEN Q, GEN m, GEN a4, GEN p)
610 : {
611 203 : if (ell_is_inf(P) || ell_is_inf(Q)) return gen_1;
612 203 : if (lgefint(p)==3 && lgefint(m)==3)
613 : {
614 203 : pari_sp av = avma;
615 203 : ulong pp = p[2];
616 203 : GEN Pp = ZV_to_Flv(P, pp), Qp = ZV_to_Flv(Q, pp);
617 203 : ulong w = Fle_tatepairing(Pp, Qp, itou(m), umodiu(a4, pp), pp);
618 203 : return gc_utoi(av,w);
619 : }
620 0 : return FpE_Miller(P, Q, m, a4, p);
621 : }
622 :
623 : /***********************************************************************/
624 : /** **/
625 : /** CM by principal order **/
626 : /** **/
627 : /***********************************************************************/
628 :
629 : /* is jn/jd = J (mod p) */
630 : static int
631 651310 : is_CMj(long J, GEN jn, GEN jd, GEN p)
632 651310 : { return dvdii(subii(mulis(jd,J), jn), p); }
633 : #ifndef LONG_IS_64BIT
634 : /* is jn/jd = -(2^32 a + b) (mod p) */
635 : static int
636 14407 : u2_is_CMj(ulong a, ulong b, GEN jn, GEN jd, GEN p)
637 : {
638 14407 : GEN mJ = uu32toi(a,b);
639 14407 : return dvdii(addii(mulii(jd,mJ), jn), p);
640 : }
641 : #endif
642 :
643 : static long
644 52620 : Fp_ellj_get_CM(GEN jn, GEN jd, GEN p)
645 : {
646 : #define CHECK(CM,J) if (is_CMj(J,jn,jd,p)) return CM;
647 52620 : CHECK(-3, 0);
648 52495 : CHECK(-4, 1728);
649 52386 : CHECK(-7, -3375);
650 52128 : CHECK(-8, 8000);
651 51904 : CHECK(-11, -32768);
652 51661 : CHECK(-12, 54000);
653 51424 : CHECK(-16, 287496);
654 51250 : CHECK(-19, -884736);
655 51027 : CHECK(-27, -12288000);
656 50805 : CHECK(-28, 16581375);
657 50630 : CHECK(-43, -884736000);
658 : #ifdef LONG_IS_64BIT
659 43253 : CHECK(-67, -147197952000L);
660 43124 : CHECK(-163, -262537412640768000L);
661 : #else
662 7214 : if (u2_is_CMj(0x00000022UL,0x45ae8000UL,jn,jd,p)) return -67;
663 7193 : if (u2_is_CMj(0x03a4b862UL,0xc4b40000UL,jn,jd,p)) return -163;
664 : #endif
665 : #undef CHECK
666 50155 : return 0;
667 : }
668 :
669 : /***********************************************************************/
670 : /** **/
671 : /** issupersingular **/
672 : /** **/
673 : /***********************************************************************/
674 :
675 : /* assume x reduced mod p, monic. Return one root, or NULL if irreducible */
676 : static GEN
677 73563 : FqX_quad_root(GEN x, GEN T, GEN p)
678 : {
679 73563 : GEN b = gel(x,3), c = gel(x,2);
680 73563 : GEN D = Fq_sub(Fq_sqr(b, T, p), Fq_mulu(c,4, T, p), T, p);
681 73563 : GEN s = Fq_sqrt(D,T, p);
682 73563 : if (!s) return NULL;
683 70385 : return Fq_halve(Fq_sub(s, b, T, p), T, p);
684 : }
685 :
686 : static GEN
687 1229 : FpX_quad_root(GEN x, GEN p)
688 : {
689 1229 : GEN s, b = gel(x,3), c = gel(x,2);
690 1229 : GEN D = Fp_sub(Fp_sqr(b, p), shifti(c,2), p);
691 1229 : if (kronecker(D,p) == -1) return NULL;
692 781 : s = Fp_sqrt(D,p);
693 781 : return Fp_halve(Fp_sub(s, b, p), p);
694 : }
695 :
696 : /* pol is the modular polynomial of level 2 modulo p.
697 : *
698 : * (T, p) defines the field FF_{p^2} in which j_prev and j live. */
699 : static long
700 4878 : Fq_path_extends_to_floor(GEN j_prev, GEN j, GEN T, GEN p, GEN Phi2, long max_len)
701 : {
702 4878 : pari_sp ltop = avma;
703 4878 : long d, i, l = lg(j);
704 :
705 : /* A path made its way to the floor if (i) its length was cut off
706 : * before reaching max_path_len, or (ii) it reached max_path_len but
707 : * only has one neighbour. */
708 32215 : for (d = 1; d <= max_len; ++d)
709 : {
710 81746 : for (i = 1; i < l; i++)
711 : {
712 54409 : GEN Phi2_j = FqX_div_by_X_x(FqXY_evalx(Phi2, gel(j,i), T, p), gel(j_prev,i), T, p, NULL);
713 54409 : GEN j_next = FqX_quad_root(Phi2_j, T, p);
714 54409 : if (!j_next)
715 3178 : return gc_long(ltop, 1);
716 51231 : gel(j_prev,i) = gel(j, i); gel(j,i) = j_next;
717 : }
718 27337 : if (gc_needed(ltop, 2))
719 0 : gerepileall(ltop, 2, &j, &j_prev);
720 : }
721 1700 : return gc_long(ltop, 0);
722 : }
723 :
724 : static long
725 448 : Fp_path_extends_to_floor(GEN j_prev, GEN j, GEN p, GEN Phi2, long max_len, GEN *pt_j, GEN *pt_j_prev)
726 : {
727 448 : pari_sp ltop = avma;
728 448 : long d, i, l = lg(j);
729 :
730 : /* A path made its way to the floor if (i) its length was cut off
731 : * before reaching max_path_len, or (ii) it reached max_path_len but
732 : * only has one neighbour. */
733 603 : for (d = 1; d <= max_len; ++d)
734 : {
735 1384 : for (i = 1; i < l; i++)
736 : {
737 1229 : GEN Phi2_j = FpX_div_by_X_x(FpXY_evalx(Phi2, gel(j,i), p), gel(j_prev,i), p, NULL);
738 1229 : GEN j_next = FpX_quad_root(Phi2_j, p);
739 1229 : if (!j_next)
740 : {
741 448 : *pt_j = gel(j,i);
742 448 : *pt_j_prev = gel(j_prev,i);
743 448 : return 1;
744 : }
745 781 : gel(j_prev,i) = gel(j, i); gel(j,i) = j_next;
746 : }
747 155 : if (gc_needed(ltop, 2))
748 0 : gerepileall(ltop, 2, &j, &j_prev);
749 : }
750 0 : return gc_long(ltop, 0);
751 : }
752 :
753 :
754 : static int
755 2737 : Fp_jissupersingular(GEN j, GEN p)
756 : {
757 2737 : long max_path_len = expi(p)+1;
758 2737 : GEN Phi2 = FpXX_red(polmodular_ZXX(2,0,0,1), p);
759 2737 : GEN Phi2_j = FpXY_evalx(Phi2, j, p);
760 2737 : GEN roots = FpX_roots(Phi2_j, p);
761 2737 : long nbroots = lg(roots)-1;
762 2737 : GEN S, j_prev = NULL;
763 :
764 : /* Every node in a supersingular L-volcano has L + 1 neighbours. */
765 : /* Note: a multiple root only occur when j has CM by sqrt(-15). */
766 2737 : if (nbroots==0)
767 665 : return 0;
768 2072 : S = deg2pol_shallow(gen_1, gen_0, Fp_neg(Fp_2gener(p),p),1);
769 2072 : if (nbroots==1 && FpX_is_squarefree(Phi2_j, p))
770 1624 : { j_prev = j; j = FqX_quad_root(FpX_div_by_X_x(Phi2_j, gel(roots,1), p, NULL), S, p); }
771 : else
772 448 : if (!Fp_path_extends_to_floor(const_vec(nbroots,j), roots, p, Phi2, max_path_len, &j, &j_prev))
773 0 : return 1;
774 2072 : return !Fq_path_extends_to_floor(mkvec(j_prev), mkvec(j), S, p, Phi2, max_path_len);
775 : }
776 :
777 : static int
778 14055 : jissupersingular(GEN j, GEN S, GEN p)
779 : {
780 14055 : long max_path_len = expi(p)+1;
781 14055 : GEN Phi2 = FpXX_red(polmodular_ZXX(2,0,0,1), p);
782 14055 : GEN Phi2_j = FqXY_evalx(Phi2, j, S, p);
783 14055 : GEN roots = FpXQX_roots(Phi2_j, S, p);
784 14055 : long nbroots = lg(roots)-1;
785 :
786 : /* Every node in a supersingular L-volcano has L + 1 neighbours. */
787 : /* Note: a multiple root only occur when j has CM by sqrt(-15). */
788 14055 : if (nbroots==0 || (nbroots==1 && FqX_is_squarefree(Phi2_j, S, p)))
789 11249 : return 0;
790 : else
791 2806 : return !Fq_path_extends_to_floor(const_vec(nbroots,j), roots, S, p, Phi2, max_path_len);
792 : }
793 :
794 : int
795 3711 : Fp_elljissupersingular(GEN j, GEN p)
796 : {
797 : long CM;
798 3711 : if (abscmpiu(p, 5) <= 0) return signe(j) == 0; /* valid if p <= 5 */
799 3571 : CM = Fp_ellj_get_CM(j, gen_1, p);
800 3571 : if (CM < 0) return krosi(CM, p) < 0; /* valid if p > 3 */
801 : else
802 2737 : return Fp_jissupersingular(j, p);
803 : }
804 :
805 : /***********************************************************************/
806 : /** **/
807 : /** Cardinal **/
808 : /** **/
809 : /***********************************************************************/
810 :
811 : /*assume a4,a6 reduced mod p odd */
812 : static ulong
813 723360 : Fl_elltrace_naive(ulong a4, ulong a6, ulong p)
814 : {
815 : ulong i, j;
816 723360 : long a = 0;
817 : long d0, d1, d2, d3;
818 723360 : GEN k = const_vecsmall(p, -1);
819 723384 : k[1] = 0;
820 129234085 : for (i=1, j=1; i < p; i += 2, j = Fl_add(j, i, p))
821 128510793 : k[j+1] = 1;
822 723292 : d0 = 6%p; d1 = d0; d2 = Fl_add(a4, 1, p); d3 = a6;
823 723290 : for(i=0;; i++)
824 : {
825 253454045 : a -= k[1+d3];
826 253454045 : if (i==p-1) break;
827 252730750 : d3 = Fl_add(d3, d2, p);
828 252732093 : d2 = Fl_add(d2, d1, p);
829 252743713 : d1 = Fl_add(d1, d0, p);
830 : }
831 723295 : return a;
832 : }
833 :
834 : /* z1 <-- z1 + z2, with precomputed inverse */
835 : static void
836 305694 : FpE_add_ip(GEN z1, GEN z2, GEN a4, GEN p, GEN p2inv)
837 : {
838 : GEN p1,x,x1,x2,y,y1,y2;
839 :
840 305694 : x1 = gel(z1,1); y1 = gel(z1,2);
841 305694 : x2 = gel(z2,1); y2 = gel(z2,2);
842 305694 : if (x1 == x2)
843 67 : p1 = Fp_add(a4, mulii(x1,mului(3,x1)), p);
844 : else
845 305627 : p1 = Fp_sub(y2,y1, p);
846 :
847 305694 : p1 = Fp_mul(p1, p2inv, p);
848 305694 : x = Fp_sub(sqri(p1), addii(x1,x2), p);
849 305694 : y = Fp_sub(mulii(p1,subii(x1,x)), y1, p);
850 305694 : affii(x, x1);
851 305694 : affii(y, y1);
852 305694 : }
853 :
854 : /* make sure *x has lgefint >= k */
855 : static void
856 19038 : _fix(GEN x, long k)
857 : {
858 19038 : GEN y = (GEN)*x;
859 19038 : if (lgefint(y) < k) { GEN p1 = cgeti(k); affii(y,p1); *x = (long)p1; }
860 19038 : }
861 :
862 : /* Return the lift of a (mod b), which is closest to c */
863 : static GEN
864 254894 : closest_lift(GEN a, GEN b, GEN c)
865 : {
866 254894 : return addii(a, mulii(b, diviiround(subii(c,a), b)));
867 : }
868 :
869 : static long
870 78 : get_table_size(GEN pordmin, GEN B)
871 : {
872 78 : pari_sp av = avma;
873 78 : GEN t = ceilr( sqrtr( divri(itor(pordmin, DEFAULTPREC), B) ) );
874 78 : if (is_bigint(t))
875 0 : pari_err_OVERFLOW("ellap [large prime: install the 'seadata' package]");
876 78 : set_avma(av);
877 78 : return itos(t) >> 1;
878 : }
879 :
880 : /* Find x such that kronecker(u = x^3+c4x+c6, p) is KRO.
881 : * Return point [x*u,u^2] on E (KRO=1) / E^twist (KRO=-1) */
882 : static GEN
883 0 : Fp_ellpoint(long KRO, ulong *px, GEN c4, GEN c6, GEN p)
884 : {
885 0 : ulong x = *px;
886 : GEN u;
887 : for(;;)
888 : {
889 0 : x++; /* u = x^3 + c4 x + c6 */
890 0 : u = modii(addii(c6, mului(x, addii(c4, sqru(x)))), p);
891 0 : if (kronecker(u,p) == KRO) break;
892 : }
893 0 : *px = x;
894 0 : return mkvec2(modii(mului(x,u),p), Fp_sqr(u,p));
895 : }
896 : static GEN
897 7227 : Fl_ellpoint(long KRO, ulong *px, ulong c4, ulong c6, ulong p)
898 : {
899 7227 : ulong t, u, x = *px;
900 : for(;;)
901 : {
902 14220 : if (++x >= p) pari_err_PRIME("ellap",utoi(p));
903 14220 : t = Fl_add(c4, Fl_sqr(x,p), p);
904 14219 : u = Fl_add(c6, Fl_mul(x, t, p), p);
905 14219 : if (krouu(u,p) == KRO) break;
906 : }
907 7227 : *px = x;
908 7227 : return mkvecsmall2(Fl_mul(x,u,p), Fl_sqr(u,p));
909 : }
910 :
911 : static GEN ap_j1728(GEN a4,GEN p);
912 : /* compute a_p using Shanks/Mestre + Montgomery's trick. Assume p > 457 */
913 : static GEN
914 78 : Fp_ellcard_Shanks(GEN c4, GEN c6, GEN p)
915 : {
916 : pari_timer T;
917 : long *tx, *ty, *ti, pfinal, i, j, s, KRO, nb;
918 : ulong x;
919 78 : pari_sp av = avma, av2;
920 : GEN p1, P, mfh, h, F,f, fh,fg, pordmin, u, v, p1p, p2p, A, B, a4, pts;
921 78 : tx = NULL;
922 78 : ty = ti = NULL; /* gcc -Wall */
923 :
924 78 : if (!signe(c6)) {
925 0 : GEN ap = ap_j1728(c4, p);
926 0 : return gerepileuptoint(av, subii(addiu(p,1), ap));
927 : }
928 :
929 78 : if (DEBUGLEVEL >= 6) timer_start(&T);
930 : /* once #E(Fp) is know mod B >= pordmin, it is completely determined */
931 78 : pordmin = addiu(sqrti(gmul2n(p,4)), 1); /* ceil( 4sqrt(p) ) */
932 78 : p1p = addiu(p, 1);
933 78 : p2p = shifti(p1p, 1);
934 78 : x = 0; KRO = 0;
935 : /* how many 2-torsion points ? */
936 78 : switch(FpX_nbroots(mkpoln(4, gen_1, gen_0, c4, c6), p))
937 : {
938 9 : case 3: A = gen_0; B = utoipos(4); break;
939 31 : case 1: A = gen_0; B = gen_2; break;
940 38 : default: A = gen_1; B = gen_2; break; /* 0 */
941 : }
942 : for(;;)
943 : {
944 78 : h = closest_lift(A, B, p1p);
945 78 : if (!KRO) /* first time, initialize */
946 : {
947 78 : KRO = kronecker(c6,p);
948 78 : f = mkvec2(gen_0, Fp_sqr(c6,p));
949 : }
950 : else
951 : {
952 0 : KRO = -KRO;
953 0 : f = Fp_ellpoint(KRO, &x, c4,c6,p);
954 : }
955 : /* [ux, u^2] is on E_u: y^2 = x^3 + c4 u^2 x + c6 u^3
956 : * E_u isomorphic to E (resp. E') iff KRO = 1 (resp. -1)
957 : * #E(F_p) = p+1 - a_p, #E'(F_p) = p+1 + a_p
958 : *
959 : * #E_u(Fp) = A (mod B), h is close to #E_u(Fp) */
960 78 : a4 = modii(mulii(c4, gel(f,2)), p); /* c4 for E_u */
961 78 : fh = FpE_mul(f, h, a4, p);
962 78 : if (ell_is_inf(fh)) goto FOUND;
963 :
964 78 : s = get_table_size(pordmin, B);
965 : /* look for h s.t f^h = 0 */
966 78 : if (!tx)
967 : { /* first time: initialize */
968 78 : tx = newblock(3*(s+1));
969 78 : ty = tx + (s+1);
970 78 : ti = ty + (s+1);
971 : }
972 78 : F = FpE_mul(f,B,a4,p);
973 78 : *tx = evaltyp(t_VECSMALL) | evallg(s+1);
974 :
975 : /* F = B.f */
976 78 : P = gcopy(fh);
977 78 : if (s < 3)
978 : { /* we're nearly done: naive search */
979 0 : GEN q1 = P, mF = FpE_neg(F, p); /* -F */
980 0 : for (i=1;; i++)
981 : {
982 0 : P = FpE_add(P,F,a4,p); /* h.f + i.F */
983 0 : if (ell_is_inf(P)) { h = addii(h, mului(i,B)); goto FOUND; }
984 0 : q1 = FpE_add(q1,mF,a4,p); /* h.f - i.F */
985 0 : if (ell_is_inf(q1)) { h = subii(h, mului(i,B)); goto FOUND; }
986 : }
987 : }
988 : /* Baby Step/Giant Step */
989 78 : nb = minss(128, s >> 1); /* > 0. Will do nb pts at a time: faster inverse */
990 78 : pts = cgetg(nb+1, t_VEC);
991 78 : j = lgefint(p);
992 9597 : for (i=1; i<=nb; i++)
993 : { /* baby steps */
994 9519 : gel(pts,i) = P; /* h.f + (i-1).F */
995 9519 : _fix(P+1, j); tx[i] = mod2BIL(gel(P,1));
996 9519 : _fix(P+2, j); ty[i] = mod2BIL(gel(P,2));
997 9519 : P = FpE_add(P,F,a4,p); /* h.f + i.F */
998 9519 : if (ell_is_inf(P)) { h = addii(h, mului(i,B)); goto FOUND; }
999 : }
1000 78 : mfh = FpE_neg(fh, p);
1001 78 : fg = FpE_add(P,mfh,a4,p); /* h.f + nb.F - h.f = nb.F */
1002 78 : if (ell_is_inf(fg)) { h = mului(nb,B); goto FOUND; }
1003 78 : u = cgetg(nb+1, t_VEC);
1004 78 : av2 = avma; /* more baby steps, nb points at a time */
1005 1356 : while (i <= s)
1006 : {
1007 : long maxj;
1008 164235 : for (j=1; j<=nb; j++) /* adding nb.F (part 1) */
1009 : {
1010 162957 : P = gel(pts,j); /* h.f + (i-nb-1+j-1).F */
1011 162957 : gel(u,j) = subii(gel(fg,1), gel(P,1));
1012 162957 : if (!signe(gel(u,j))) /* sum = 0 or doubling */
1013 : {
1014 1 : long k = i+j-2;
1015 1 : if (equalii(gel(P,2),gel(fg,2))) k -= 2*nb; /* fg == P */
1016 1 : h = addii(h, mulsi(k,B)); goto FOUND;
1017 : }
1018 : }
1019 1278 : v = FpV_inv(u, p);
1020 1278 : maxj = (i-1 + nb <= s)? nb: s % nb;
1021 160545 : for (j=1; j<=maxj; j++,i++) /* adding nb.F (part 2) */
1022 : {
1023 159267 : P = gel(pts,j);
1024 159267 : FpE_add_ip(P,fg, a4,p, gel(v,j));
1025 159267 : tx[i] = mod2BIL(gel(P,1));
1026 159267 : ty[i] = mod2BIL(gel(P,2));
1027 : }
1028 1278 : set_avma(av2);
1029 : }
1030 77 : P = FpE_add(gel(pts,j-1),mfh,a4,p); /* = (s-1).F */
1031 77 : if (ell_is_inf(P)) { h = mului(s-1,B); goto FOUND; }
1032 77 : if (DEBUGLEVEL >= 6)
1033 0 : timer_printf(&T, "[Fp_ellcard_Shanks] baby steps, s = %ld",s);
1034 :
1035 : /* giant steps: fg = s.F */
1036 77 : fg = FpE_add(P,F,a4,p);
1037 77 : if (ell_is_inf(fg)) { h = mului(s,B); goto FOUND; }
1038 77 : pfinal = mod2BIL(p); av2 = avma;
1039 : /* Goal of the following: sort points by increasing x-coordinate hash.
1040 : * Done in a complicated way to avoid allocating a large temp vector */
1041 77 : p1 = vecsmall_indexsort(tx); /* = permutation sorting tx */
1042 168784 : for (i=1; i<=s; i++) ti[i] = tx[p1[i]];
1043 : /* ti = tx sorted */
1044 168784 : for (i=1; i<=s; i++) { tx[i] = ti[i]; ti[i] = ty[p1[i]]; }
1045 : /* tx is sorted. ti = ty sorted */
1046 168784 : for (i=1; i<=s; i++) { ty[i] = ti[i]; ti[i] = p1[i]; }
1047 : /* ty is sorted. ti = permutation sorting tx */
1048 77 : if (DEBUGLEVEL >= 6) timer_printf(&T, "[Fp_ellcard_Shanks] sorting");
1049 77 : set_avma(av2);
1050 :
1051 77 : gaffect(fg, gel(pts,1));
1052 9440 : for (j=2; j<=nb; j++) /* pts[j] = j.fg = (s*j).F */
1053 : {
1054 9363 : P = FpE_add(gel(pts,j-1),fg,a4,p);
1055 9363 : if (ell_is_inf(P)) { h = mulii(mulss(s,j), B); goto FOUND; }
1056 9363 : gaffect(P, gel(pts,j));
1057 : }
1058 : /* replace fg by nb.fg since we do nb points at a time */
1059 77 : set_avma(av2);
1060 77 : fg = gcopy(gel(pts,nb)); /* copy: we modify (temporarily) pts[nb] below */
1061 77 : av2 = avma;
1062 :
1063 77 : for (i=1,j=1; ; i++)
1064 152075 : {
1065 152152 : GEN ftest = gel(pts,j);
1066 152152 : long m, l = 1, r = s+1;
1067 : long k, k2, j2;
1068 :
1069 152152 : set_avma(av2);
1070 152152 : k = mod2BIL(gel(ftest,1));
1071 1930966 : while (l < r)
1072 : {
1073 1778814 : m = (l+r) >> 1;
1074 1778814 : if (tx[m] < k) l = m+1; else r = m;
1075 : }
1076 152152 : if (r <= s && tx[r] == k)
1077 : {
1078 154 : while (r && tx[r] == k) r--;
1079 77 : k2 = mod2BIL(gel(ftest,2));
1080 77 : for (r++; r <= s && tx[r] == k; r++)
1081 77 : if (ty[r] == k2 || ty[r] == pfinal - k2)
1082 : { /* [h+j2] f == +/- ftest (= [i.s] f)? */
1083 77 : j2 = ti[r] - 1;
1084 77 : if (DEBUGLEVEL >=6)
1085 0 : timer_printf(&T, "[Fp_ellcard_Shanks] giant steps, i = %ld",i);
1086 77 : P = FpE_add(FpE_mul(F,stoi(j2),a4,p),fh,a4,p);
1087 77 : if (equalii(gel(P,1), gel(ftest,1)))
1088 : {
1089 77 : if (equalii(gel(P,2), gel(ftest,2))) i = -i;
1090 77 : h = addii(h, mulii(addis(mulss(s,i), j2), B));
1091 77 : goto FOUND;
1092 : }
1093 : }
1094 : }
1095 152075 : if (++j > nb)
1096 : { /* compute next nb points */
1097 1149 : long save = 0; /* gcc -Wall */;
1098 147576 : for (j=1; j<=nb; j++)
1099 : {
1100 146427 : P = gel(pts,j);
1101 146427 : gel(u,j) = subii(gel(fg,1), gel(P,1));
1102 146427 : if (gel(u,j) == gen_0) /* occurs once: i = j = nb, P == fg */
1103 : {
1104 67 : gel(u,j) = shifti(gel(P,2),1);
1105 67 : save = fg[1]; fg[1] = P[1];
1106 : }
1107 : }
1108 1149 : v = FpV_inv(u, p);
1109 147576 : for (j=1; j<=nb; j++)
1110 146427 : FpE_add_ip(gel(pts,j),fg,a4,p, gel(v,j));
1111 1149 : if (i == nb) { fg[1] = save; }
1112 1149 : j = 1;
1113 : }
1114 : }
1115 78 : FOUND: /* found a point of exponent h on E_u */
1116 78 : h = FpE_order(f, h, a4, p);
1117 : /* h | #E_u(Fp) = A (mod B) */
1118 78 : A = Z_chinese_all(A, gen_0, B, h, &B);
1119 78 : if (cmpii(B, pordmin) >= 0) break;
1120 : /* not done: update A mod B for the _next_ curve, isomorphic to
1121 : * the quadratic twist of this one */
1122 0 : A = remii(subii(p2p,A), B); /* #E(Fp)+#E'(Fp) = 2p+2 */
1123 : }
1124 78 : if (tx) killblock(tx);
1125 78 : h = closest_lift(A, B, p1p);
1126 78 : return gerepileuptoint(av, KRO==1? h: subii(p2p,h));
1127 : }
1128 :
1129 : typedef struct
1130 : {
1131 : ulong x,y,i;
1132 : } multiple;
1133 :
1134 : static int
1135 15372502 : compare_multiples(multiple *a, multiple *b) { return a->x > b->x? 1:a->x<b->x?-1:0; }
1136 :
1137 : /* find x such that h := a + b x is closest to c and return h:
1138 : * x = round((c-a) / b) = floor( (2(c-a) + b) / 2b )
1139 : * Assume 0 <= a < b < c and b + 2c < 2^BIL */
1140 : static ulong
1141 261971 : uclosest_lift(ulong a, ulong b, ulong c)
1142 : {
1143 261971 : ulong x = (b + ((c-a) << 1)) / (b << 1);
1144 261971 : return a + b * x;
1145 : }
1146 :
1147 : static long
1148 227178 : Fle_dbl_inplace(GEN P, ulong a4, ulong p)
1149 : {
1150 : ulong x, y, slope;
1151 227178 : if (!P[2]) return 1;
1152 227150 : x = P[1]; y = P[2];
1153 227150 : slope = Fl_div(Fl_add(Fl_triple(Fl_sqr(x,p), p), a4, p),
1154 : Fl_double(y, p), p);
1155 227153 : P[1] = Fl_sub(Fl_sqr(slope, p), Fl_double(x, p), p);
1156 227150 : P[2] = Fl_sub(Fl_mul(slope, Fl_sub(x, P[1], p), p), y, p);
1157 227149 : return 0;
1158 : }
1159 :
1160 : static long
1161 5795603 : Fle_add_inplace(GEN P, GEN Q, ulong a4, ulong p)
1162 : {
1163 : ulong Px, Py, Qx, Qy, slope;
1164 5795603 : if (ell_is_inf(Q)) return 0;
1165 5795626 : Px = P[1]; Py = P[2];
1166 5795626 : Qx = Q[1]; Qy = Q[2];
1167 5795626 : if (Px==Qx)
1168 238639 : return Py==Qy ? Fle_dbl_inplace(P, a4, p): 1;
1169 5556987 : slope = Fl_div(Fl_sub(Py, Qy, p), Fl_sub(Px, Qx, p), p);
1170 5558004 : P[1] = Fl_sub(Fl_sub(Fl_sqr(slope, p), Px, p), Qx, p);
1171 5557100 : P[2] = Fl_sub(Fl_mul(slope, Fl_sub(Px, P[1], p), p), Py, p);
1172 5556643 : return 0;
1173 : }
1174 :
1175 : /* assume 99 < p < 2^(BIL-1) - 2^((BIL+1)/2) and e has good reduction at p.
1176 : * Should use Barett reduction + multi-inverse. See Fp_ellcard_Shanks() */
1177 : static long
1178 254746 : Fl_ellcard_Shanks(ulong c4, ulong c6, ulong p)
1179 : {
1180 : GEN f, fh, fg, ftest, F;
1181 : ulong i, l, r, s, h, x, cp4, p1p, p2p, pordmin,A,B;
1182 : long KRO;
1183 254746 : pari_sp av = avma;
1184 : multiple *table;
1185 :
1186 254746 : if (!c6) {
1187 14 : GEN ap = ap_j1728(utoi(c4), utoipos(p));
1188 14 : return gc_long(av, p+1 - itos(ap));
1189 : }
1190 :
1191 254732 : pordmin = (ulong)(1 + 4*sqrt((double)p));
1192 254732 : p1p = p+1;
1193 254732 : p2p = p1p << 1;
1194 254732 : x = 0; KRO = 0;
1195 254732 : switch(Flx_nbroots(mkvecsmall5(0L, c6,c4,0L,1L), p))
1196 : {
1197 51715 : case 3: A = 0; B = 4; break;
1198 124410 : case 1: A = 0; B = 2; break;
1199 78620 : default: A = 1; B = 2; break; /* 0 */
1200 : }
1201 : for(;;)
1202 : { /* see comments in Fp_ellcard_Shanks */
1203 261972 : h = uclosest_lift(A, B, p1p);
1204 261971 : if (!KRO) /* first time, initialize */
1205 : {
1206 254744 : KRO = krouu(c6,p); /* != 0 */
1207 254742 : f = mkvecsmall2(0, Fl_sqr(c6,p));
1208 : }
1209 : else
1210 : {
1211 7227 : KRO = -KRO;
1212 7227 : f = Fl_ellpoint(KRO, &x, c4,c6,p);
1213 : }
1214 261970 : cp4 = Fl_mul(c4, f[2], p);
1215 261970 : fh = Fle_mulu(f, h, cp4, p);
1216 261955 : if (ell_is_inf(fh)) goto FOUND;
1217 :
1218 255749 : s = (ulong) (sqrt(((double)pordmin)/B) / 2);
1219 255749 : if (!s) s = 1;
1220 255749 : table = (multiple *) stack_malloc((s+1) * sizeof(multiple));
1221 255752 : F = Fle_mulu(f, B, cp4, p);
1222 3346203 : for (i=0; i < s; i++)
1223 : {
1224 3101934 : table[i].x = fh[1];
1225 3101934 : table[i].y = fh[2];
1226 3101934 : table[i].i = i;
1227 3101934 : if (Fle_add_inplace(fh, F, cp4, p)) { h += B*(i+1); goto FOUND; }
1228 : }
1229 244269 : qsort(table,s,sizeof(multiple),(QSCOMP)compare_multiples);
1230 244278 : fg = Fle_mulu(F, s, cp4, p); ftest = zv_copy(fg);
1231 244270 : if (ell_is_inf(ftest)) {
1232 0 : if (!uisprime(p)) pari_err_PRIME("ellap",utoi(p));
1233 0 : pari_err_BUG("ellap (f^(i*s) = 1)");
1234 : }
1235 2938614 : for (i=1; ; i++)
1236 : {
1237 2938614 : l=0; r=s;
1238 20639107 : while (l<r)
1239 : {
1240 17700493 : ulong m = (l+r) >> 1;
1241 17700493 : if (table[m].x < uel(ftest,1)) l=m+1; else r=m;
1242 : }
1243 2938614 : if (r < s && table[r].x == uel(ftest,1)) break;
1244 2694344 : if (Fle_add_inplace(ftest, fg, cp4, p)) pari_err_PRIME("ellap",utoi(p));
1245 : }
1246 244270 : h += table[r].i * B;
1247 244270 : if (table[r].y == uel(ftest,2))
1248 126869 : h -= s * i * B;
1249 : else
1250 117401 : h += s * i * B;
1251 261966 : FOUND:
1252 261966 : h = itou(Fle_order(f, utoipos(h), cp4, p));
1253 : /* h | #E_u(Fp) = A (mod B) */
1254 : {
1255 : GEN C;
1256 261964 : A = itou( Z_chinese_all(gen_0, utoi(A), utoipos(h), utoipos(B), &C) );
1257 261963 : if (abscmpiu(C, pordmin) >= 0) { /* uclosest_lift could overflow */
1258 254737 : h = itou( closest_lift(utoi(A), C, utoipos(p1p)) );
1259 254741 : break;
1260 : }
1261 7227 : B = itou(C);
1262 : }
1263 7227 : A = (p2p - A) % B; set_avma(av);
1264 : }
1265 254741 : return gc_long(av, KRO==1? h: p2p-h);
1266 : }
1267 :
1268 : /** ellap from CM (original code contributed by Mark Watkins) **/
1269 :
1270 : static GEN
1271 85147 : ap_j0(GEN a6,GEN p)
1272 : {
1273 : GEN a, b, e, d;
1274 85147 : if (umodiu(p,3) != 1) return gen_0;
1275 42281 : (void)cornacchia2(utoipos(27),p, &a,&b);
1276 42435 : if (umodiu(a, 3) == 1) a = negi(a);
1277 42435 : d = mulis(a6,-108);
1278 42372 : e = diviuexact(shifti(p,-1), 3); /* (p-1) / 6 */
1279 42303 : return centermod(mulii(a, Fp_pow(d, e, p)), p);
1280 : }
1281 : static GEN
1282 2642444 : ap_j1728(GEN a4,GEN p)
1283 : {
1284 : GEN a, b, e;
1285 2642444 : if (mod4(p) != 1) return gen_0;
1286 1320221 : (void)cornacchia2(utoipos(4),p, &a,&b);
1287 1320221 : if (Mod4(a)==0) a = b;
1288 1320221 : if (Mod2(a)==1) a = shifti(a,1);
1289 1320221 : if (Mod8(a)==6) a = negi(a);
1290 1320221 : e = shifti(p,-2); /* (p-1) / 4 */
1291 1320221 : return centermod(mulii(a, Fp_pow(a4, e, p)), p);
1292 : }
1293 : static GEN
1294 126 : ap_j8000(GEN a6, GEN p)
1295 : {
1296 : GEN a, b;
1297 126 : long r = mod8(p), s = 1;
1298 126 : if (r != 1 && r != 3) return gen_0;
1299 49 : (void)cornacchia2(utoipos(8),p, &a,&b);
1300 49 : switch(Mod16(a)) {
1301 14 : case 2: case 6: if (Mod4(b)) s = -s;
1302 14 : break;
1303 35 : case 10: case 14: if (!Mod4(b)) s = -s;
1304 35 : break;
1305 : }
1306 49 : if (kronecker(mulis(a6, 42), p) < 0) s = -s;
1307 49 : return s > 0? a: negi(a);
1308 : }
1309 : static GEN
1310 140 : ap_j287496(GEN a6, GEN p)
1311 : {
1312 : GEN a, b;
1313 140 : long s = 1;
1314 140 : if (mod4(p) != 1) return gen_0;
1315 70 : (void)cornacchia2(utoipos(4),p, &a,&b);
1316 70 : if (Mod4(a)==0) a = b;
1317 70 : if (Mod2(a)==1) a = shifti(a,1);
1318 70 : if (Mod8(a)==6) s = -s;
1319 70 : if (krosi(2,p) < 0) s = -s;
1320 70 : if (kronecker(mulis(a6, -14), p) < 0) s = -s;
1321 70 : return s > 0? a: negi(a);
1322 : }
1323 : static GEN
1324 1344 : ap_cm(int CM, long A6B, GEN a6, GEN p)
1325 : {
1326 : GEN a, b;
1327 1344 : long s = 1;
1328 1344 : if (krosi(CM,p) < 0) return gen_0;
1329 644 : (void)cornacchia2(utoipos(-CM),p, &a, &b);
1330 644 : if ((CM&3) == 0) CM >>= 2;
1331 644 : if ((krois(a, -CM) > 0) ^ (CM == -7)) s = -s;
1332 644 : if (kronecker(mulis(a6,A6B), p) < 0) s = -s;
1333 644 : return s > 0? a: negi(a);
1334 : }
1335 : static GEN
1336 497483 : ec_ap_cm(int CM, GEN a4, GEN a6, GEN p)
1337 : {
1338 497483 : switch(CM)
1339 : {
1340 29113 : case -3: return ap_j0(a6, p);
1341 466760 : case -4: return ap_j1728(a4, p);
1342 126 : case -8: return ap_j8000(a6, p);
1343 140 : case -16: return ap_j287496(a6, p);
1344 154 : case -7: return ap_cm(CM, -2, a6, p);
1345 147 : case -11: return ap_cm(CM, 21, a6, p);
1346 168 : case -12: return ap_cm(CM, 22, a6, p);
1347 147 : case -19: return ap_cm(CM, 1, a6, p);
1348 154 : case -27: return ap_cm(CM, 253, a6, p);
1349 140 : case -28: return ap_cm(-7, -114, a6, p); /* yes, -7 ! */
1350 147 : case -43: return ap_cm(CM, 21, a6, p);
1351 147 : case -67: return ap_cm(CM, 217, a6, p);
1352 140 : case -163:return ap_cm(CM, 185801, a6, p);
1353 0 : default: return NULL;
1354 : }
1355 : }
1356 :
1357 : static GEN
1358 49146 : Fp_ellj_nodiv(GEN a4, GEN a6, GEN p)
1359 : {
1360 49146 : GEN a43 = Fp_mulu(Fp_powu(a4, 3, p), 4, p);
1361 49152 : GEN a62 = Fp_mulu(Fp_sqr(a6, p), 27, p);
1362 49154 : return mkvec2(Fp_mulu(a43, 1728, p), Fp_add(a43, a62, p));
1363 : }
1364 :
1365 : GEN
1366 98 : Fp_ellj(GEN a4, GEN a6, GEN p)
1367 : {
1368 98 : pari_sp av = avma;
1369 : GEN z;
1370 98 : if (lgefint(p) == 3)
1371 : {
1372 0 : ulong pp = p[2];
1373 0 : return utoi(Fl_ellj(umodiu(a4,pp), umodiu(a6,pp), pp));
1374 : }
1375 98 : z = Fp_ellj_nodiv(a4, a6, p);
1376 98 : return gerepileuptoint(av,Fp_div(gel(z,1),gel(z,2),p));
1377 : }
1378 :
1379 : void
1380 1106 : Fp_ellj_to_a4a6(GEN j, GEN p, GEN *pt_a4, GEN *pt_a6)
1381 : {
1382 1106 : j = modii(j, p);
1383 1106 : if (signe(j) == 0) { *pt_a4 = gen_0; *pt_a6 = gen_1; }
1384 791 : else if (equaliu(j,umodui(1728,p))) { *pt_a4 = gen_1; *pt_a6 = gen_0; }
1385 : else
1386 : {
1387 616 : GEN k = Fp_sub(utoi(1728), j, p);
1388 616 : GEN kj = Fp_mul(k, j, p);
1389 616 : GEN k2j = Fp_mul(kj, k, p);
1390 616 : *pt_a4 = Fp_mulu(kj, 3, p);
1391 616 : *pt_a6 = Fp_double(k2j, p);
1392 : }
1393 1106 : }
1394 :
1395 : static GEN /* Only compute a mod p, so assume p>=17 */
1396 2280696 : Fp_ellcard_CM(GEN a4, GEN a6, GEN p)
1397 : {
1398 2280696 : pari_sp av = avma;
1399 : GEN a;
1400 2280696 : if (!signe(a4)) a = ap_j0(a6,p);
1401 2224681 : else if (!signe(a6)) a = ap_j1728(a4,p);
1402 : else
1403 : {
1404 49011 : GEN j = Fp_ellj_nodiv(a4, a6, p);
1405 49049 : long CM = Fp_ellj_get_CM(gel(j,1), gel(j,2), p);
1406 49027 : if (!CM) return gc_NULL(av);
1407 1610 : a = ec_ap_cm(CM,a4,a6,p);
1408 : }
1409 2233408 : return gerepileuptoint(av, subii(addiu(p,1),a));
1410 : }
1411 :
1412 : GEN
1413 2542433 : Fp_ellcard(GEN a4, GEN a6, GEN p)
1414 : {
1415 2542433 : long lp = expi(p);
1416 2542405 : ulong pp = p[2];
1417 2542405 : if (lp < 11)
1418 261776 : return utoi(pp+1 - Fl_elltrace_naive(umodiu(a4,pp), umodiu(a6,pp), pp));
1419 2280629 : { GEN a = Fp_ellcard_CM(a4,a6,p); if (a) return a; }
1420 47418 : if (lp >= 56)
1421 868 : return Fp_ellcard_SEA(a4, a6, p, 0);
1422 46550 : if (lp <= BITS_IN_LONG-2)
1423 46471 : return utoi(Fl_ellcard_Shanks(umodiu(a4,pp), umodiu(a6,pp), pp));
1424 79 : return Fp_ellcard_Shanks(a4, a6, p);
1425 : }
1426 :
1427 : long
1428 621591 : Fl_elltrace(ulong a4, ulong a6, ulong p)
1429 : {
1430 : pari_sp av;
1431 : long lp;
1432 : GEN a;
1433 621591 : if (p < (1<<11)) return Fl_elltrace_naive(a4, a6, p);
1434 208254 : lp = expu(p);
1435 208254 : if (lp <= minss(56, BITS_IN_LONG-2)) return p+1-Fl_ellcard_Shanks(a4, a6, p);
1436 0 : av = avma; a = subui(p+1, Fp_ellcard(utoi(a4), utoi(a6), utoipos(p)));
1437 0 : return gc_long(av, itos(a));
1438 : }
1439 : long
1440 1164937 : Fl_elltrace_CM(long CM, ulong a4, ulong a6, ulong p)
1441 : {
1442 : pari_sp av;
1443 : GEN a;
1444 1164937 : if (!CM) return Fl_elltrace(a4,a6,p);
1445 544110 : if (p < (1<<11)) return Fl_elltrace_naive(a4, a6, p);
1446 495873 : av = avma; a = ec_ap_cm(CM, utoi(a4), utoi(a6), utoipos(p));
1447 495873 : return gc_long(av, itos(a));
1448 : }
1449 :
1450 : static GEN
1451 73669 : _FpE_pairorder(void *E, GEN P, GEN Q, GEN m, GEN F)
1452 : {
1453 73669 : struct _FpE *e = (struct _FpE *) E;
1454 73669 : return Fp_order(FpE_weilpairing(P,Q,m,e->a4,e->p), F, e->p);
1455 : }
1456 :
1457 : GEN
1458 120715 : Fp_ellgroup(GEN a4, GEN a6, GEN N, GEN p, GEN *pt_m)
1459 : {
1460 : struct _FpE e;
1461 120715 : e.a4=a4; e.a6=a6; e.p=p;
1462 120715 : return gen_ellgroup(N, subiu(p,1), pt_m, (void*)&e, &FpE_group, _FpE_pairorder);
1463 : }
1464 :
1465 : GEN
1466 574 : Fp_ellgens(GEN a4, GEN a6, GEN ch, GEN D, GEN m, GEN p)
1467 : {
1468 : GEN P;
1469 574 : pari_sp av = avma;
1470 : struct _FpE e;
1471 574 : e.a4=a4; e.a6=a6; e.p=p;
1472 574 : switch(lg(D)-1)
1473 : {
1474 476 : case 1:
1475 476 : P = gen_gener(gel(D,1), (void*)&e, &FpE_group);
1476 476 : P = mkvec(FpE_changepoint(P, ch, p));
1477 476 : break;
1478 98 : default:
1479 98 : P = gen_ellgens(gel(D,1), gel(D,2), m, (void*)&e, &FpE_group, _FpE_pairorder);
1480 98 : gel(P,1) = FpE_changepoint(gel(P,1), ch, p);
1481 98 : gel(P,2) = FpE_changepoint(gel(P,2), ch, p);
1482 98 : break;
1483 : }
1484 574 : return gerepilecopy(av, P);
1485 : }
1486 :
1487 : /* Not so fast arithmetic with points over elliptic curves over FpXQ */
1488 :
1489 : /***********************************************************************/
1490 : /** **/
1491 : /** FpXQE **/
1492 : /** **/
1493 : /***********************************************************************/
1494 : /* These functions deal with point over elliptic curves over FpXQ defined
1495 : * by an equation of the form y^2=x^3+a4*x+a6.
1496 : * Most of the time a6 is omitted since it can be recovered from any point
1497 : * on the curve. */
1498 :
1499 : GEN
1500 976 : RgE_to_FpXQE(GEN x, GEN T, GEN p)
1501 : {
1502 976 : if (ell_is_inf(x)) return x;
1503 976 : retmkvec2(Rg_to_FpXQ(gel(x,1),T,p),Rg_to_FpXQ(gel(x,2),T,p));
1504 : }
1505 :
1506 : GEN
1507 1876 : FpXQE_changepoint(GEN x, GEN ch, GEN T, GEN p)
1508 : {
1509 1876 : pari_sp av = avma;
1510 : GEN p1,z,u,r,s,t,v,v2,v3;
1511 1876 : if (ell_is_inf(x)) return x;
1512 942 : u = gel(ch,1); r = gel(ch,2);
1513 942 : s = gel(ch,3); t = gel(ch,4);
1514 942 : v = FpXQ_inv(u, T, p); v2 = FpXQ_sqr(v, T, p); v3 = FpXQ_mul(v,v2, T, p);
1515 942 : p1 = FpX_sub(gel(x,1),r, p);
1516 942 : z = cgetg(3,t_VEC);
1517 942 : gel(z,1) = FpXQ_mul(v2, p1, T, p);
1518 942 : gel(z,2) = FpXQ_mul(v3, FpX_sub(gel(x,2), FpX_add(FpXQ_mul(s,p1, T, p),t, p), p), T, p);
1519 942 : return gerepileupto(av, z);
1520 : }
1521 :
1522 : GEN
1523 976 : FpXQE_changepointinv(GEN x, GEN ch, GEN T, GEN p)
1524 : {
1525 : GEN u, r, s, t, X, Y, u2, u3, u2X, z;
1526 976 : if (ell_is_inf(x)) return x;
1527 976 : X = gel(x,1); Y = gel(x,2);
1528 976 : u = gel(ch,1); r = gel(ch,2);
1529 976 : s = gel(ch,3); t = gel(ch,4);
1530 976 : u2 = FpXQ_sqr(u, T, p); u3 = FpXQ_mul(u,u2, T, p);
1531 976 : u2X = FpXQ_mul(u2,X, T, p);
1532 976 : z = cgetg(3, t_VEC);
1533 976 : gel(z,1) = FpX_add(u2X,r, p);
1534 976 : gel(z,2) = FpX_add(FpXQ_mul(u3,Y, T, p), FpX_add(FpXQ_mul(s,u2X, T, p), t, p), p);
1535 976 : return z;
1536 : }
1537 :
1538 : static GEN
1539 840 : random_nonsquare_FpXQ(GEN T, GEN p)
1540 : {
1541 840 : pari_sp av = avma;
1542 840 : long n = degpol(T), v = varn(T);
1543 : GEN a;
1544 840 : if (odd(n))
1545 : {
1546 420 : GEN z = cgetg(3, t_POL);
1547 420 : z[1] = evalsigne(1) | evalvarn(v);
1548 420 : gel(z,2) = random_nonsquare_Fp(p); return z;
1549 : }
1550 : do
1551 : {
1552 791 : set_avma(av);
1553 791 : a = random_FpX(n, v, p);
1554 791 : } while (FpXQ_issquare(a, T, p));
1555 420 : return a;
1556 : }
1557 :
1558 : void
1559 840 : FpXQ_elltwist(GEN a4, GEN a6, GEN T, GEN p, GEN *pt_a4, GEN *pt_a6)
1560 : {
1561 840 : GEN d = random_nonsquare_FpXQ(T, p);
1562 840 : GEN d2 = FpXQ_sqr(d, T, p), d3 = FpXQ_mul(d2, d, T, p);
1563 840 : *pt_a4 = FpXQ_mul(a4, d2, T, p);
1564 840 : *pt_a6 = FpXQ_mul(a6, d3, T, p);
1565 840 : }
1566 :
1567 : static GEN
1568 267975 : FpXQE_dbl_slope(GEN P, GEN a4, GEN T, GEN p, GEN *slope)
1569 : {
1570 : GEN x, y, Q;
1571 267975 : if (ell_is_inf(P) || !signe(gel(P,2))) return ellinf();
1572 266398 : x = gel(P,1); y = gel(P,2);
1573 266398 : *slope = FpXQ_div(FpX_add(FpX_mulu(FpXQ_sqr(x, T, p), 3, p), a4, p),
1574 : FpX_mulu(y, 2, p), T, p);
1575 266398 : Q = cgetg(3,t_VEC);
1576 266398 : gel(Q, 1) = FpX_sub(FpXQ_sqr(*slope, T, p), FpX_mulu(x, 2, p), p);
1577 266398 : gel(Q, 2) = FpX_sub(FpXQ_mul(*slope, FpX_sub(x, gel(Q, 1), p), T, p), y, p);
1578 266398 : return Q;
1579 : }
1580 :
1581 : GEN
1582 263873 : FpXQE_dbl(GEN P, GEN a4, GEN T, GEN p)
1583 : {
1584 263873 : pari_sp av = avma;
1585 : GEN slope;
1586 263873 : return gerepileupto(av, FpXQE_dbl_slope(P,a4,T,p,&slope));
1587 : }
1588 :
1589 : static GEN
1590 251880 : FpXQE_add_slope(GEN P, GEN Q, GEN a4, GEN T, GEN p, GEN *slope)
1591 : {
1592 : GEN Px, Py, Qx, Qy, R;
1593 251880 : if (ell_is_inf(P)) return Q;
1594 251866 : if (ell_is_inf(Q)) return P;
1595 251866 : Px = gel(P,1); Py = gel(P,2);
1596 251866 : Qx = gel(Q,1); Qy = gel(Q,2);
1597 251866 : if (ZX_equal(Px, Qx))
1598 : {
1599 674 : if (ZX_equal(Py, Qy))
1600 7 : return FpXQE_dbl_slope(P, a4, T, p, slope);
1601 : else
1602 667 : return ellinf();
1603 : }
1604 251192 : *slope = FpXQ_div(FpX_sub(Py, Qy, p), FpX_sub(Px, Qx, p), T, p);
1605 251192 : R = cgetg(3,t_VEC);
1606 251192 : gel(R, 1) = FpX_sub(FpX_sub(FpXQ_sqr(*slope, T, p), Px, p), Qx, p);
1607 251192 : gel(R, 2) = FpX_sub(FpXQ_mul(*slope, FpX_sub(Px, gel(R, 1), p), T, p), Py, p);
1608 251192 : return R;
1609 : }
1610 :
1611 : GEN
1612 251320 : FpXQE_add(GEN P, GEN Q, GEN a4, GEN T, GEN p)
1613 : {
1614 251320 : pari_sp av = avma;
1615 : GEN slope;
1616 251320 : return gerepileupto(av, FpXQE_add_slope(P,Q,a4,T,p,&slope));
1617 : }
1618 :
1619 : static GEN
1620 0 : FpXQE_neg_i(GEN P, GEN p)
1621 : {
1622 0 : if (ell_is_inf(P)) return P;
1623 0 : return mkvec2(gel(P,1), FpX_neg(gel(P,2), p));
1624 : }
1625 :
1626 : GEN
1627 73329 : FpXQE_neg(GEN P, GEN T, GEN p)
1628 : {
1629 : (void) T;
1630 73329 : if (ell_is_inf(P)) return ellinf();
1631 73329 : return mkvec2(gcopy(gel(P,1)), FpX_neg(gel(P,2), p));
1632 : }
1633 :
1634 : GEN
1635 0 : FpXQE_sub(GEN P, GEN Q, GEN a4, GEN T, GEN p)
1636 : {
1637 0 : pari_sp av = avma;
1638 : GEN slope;
1639 0 : return gerepileupto(av, FpXQE_add_slope(P, FpXQE_neg_i(Q, p), a4, T, p, &slope));
1640 : }
1641 :
1642 : struct _FpXQE { GEN a4,a6,T,p; };
1643 : static GEN
1644 263873 : _FpXQE_dbl(void *E, GEN P)
1645 : {
1646 263873 : struct _FpXQE *ell = (struct _FpXQE *) E;
1647 263873 : return FpXQE_dbl(P, ell->a4, ell->T, ell->p);
1648 : }
1649 : static GEN
1650 251320 : _FpXQE_add(void *E, GEN P, GEN Q)
1651 : {
1652 251320 : struct _FpXQE *ell=(struct _FpXQE *) E;
1653 251320 : return FpXQE_add(P, Q, ell->a4, ell->T, ell->p);
1654 : }
1655 : static GEN
1656 81862 : _FpXQE_mul(void *E, GEN P, GEN n)
1657 : {
1658 81862 : pari_sp av = avma;
1659 81862 : struct _FpXQE *e=(struct _FpXQE *) E;
1660 81862 : long s = signe(n);
1661 81862 : if (!s || ell_is_inf(P)) return ellinf();
1662 81862 : if (s<0) P = FpXQE_neg(P, e->T, e->p);
1663 81862 : if (is_pm1(n)) return s>0? gcopy(P): P;
1664 8441 : return gerepilecopy(av, gen_pow_i(P, n, e, &_FpXQE_dbl, &_FpXQE_add));
1665 : }
1666 :
1667 : GEN
1668 934 : FpXQE_mul(GEN P, GEN n, GEN a4, GEN T, GEN p)
1669 : {
1670 : struct _FpXQE E;
1671 934 : E.a4= a4; E.T = T; E.p = p;
1672 934 : return _FpXQE_mul(&E, P, n);
1673 : }
1674 :
1675 : /* Finds a random nonsingular point on E */
1676 :
1677 : GEN
1678 1072 : random_FpXQE(GEN a4, GEN a6, GEN T, GEN p)
1679 : {
1680 1072 : pari_sp ltop = avma;
1681 : GEN x, x2, y, rhs;
1682 1072 : long v = get_FpX_var(T), d = get_FpX_degree(T);
1683 : do
1684 : {
1685 2241 : set_avma(ltop);
1686 2241 : x = random_FpX(d,v,p); /* x^3+a4*x+a6 = x*(x^2+a4)+a6 */
1687 2241 : x2 = FpXQ_sqr(x, T, p);
1688 2241 : rhs = FpX_add(FpXQ_mul(x, FpX_add(x2, a4, p), T, p), a6, p);
1689 0 : } while ((!signe(rhs) && !signe(FpX_add(FpX_mulu(x2,3,p), a4, p)))
1690 2241 : || !FpXQ_issquare(rhs, T, p));
1691 1072 : y = FpXQ_sqrt(rhs, T, p);
1692 1072 : if (!y) pari_err_PRIME("random_FpE", p);
1693 1072 : return gerepilecopy(ltop, mkvec2(x, y));
1694 : }
1695 :
1696 : static GEN
1697 138 : _FpXQE_rand(void *E)
1698 : {
1699 138 : struct _FpXQE *e=(struct _FpXQE *) E;
1700 138 : return random_FpXQE(e->a4, e->a6, e->T, e->p);
1701 : }
1702 :
1703 : static const struct bb_group FpXQE_group={_FpXQE_add,_FpXQE_mul,_FpXQE_rand,hash_GEN,ZXV_equal,ell_is_inf};
1704 :
1705 : const struct bb_group *
1706 16 : get_FpXQE_group(void ** pt_E, GEN a4, GEN a6, GEN T, GEN p)
1707 : {
1708 16 : struct _FpXQE *e = (struct _FpXQE *) stack_malloc(sizeof(struct _FpXQE));
1709 16 : e->a4 = a4; e->a6 = a6; e->T = T; e->p = p;
1710 16 : *pt_E = (void *) e;
1711 16 : return &FpXQE_group;
1712 : }
1713 :
1714 : GEN
1715 14 : FpXQE_order(GEN z, GEN o, GEN a4, GEN T, GEN p)
1716 : {
1717 14 : pari_sp av = avma;
1718 : struct _FpXQE e;
1719 14 : e.a4=a4; e.T=T; e.p=p;
1720 14 : return gerepileuptoint(av, gen_order(z, o, (void*)&e, &FpXQE_group));
1721 : }
1722 :
1723 : GEN
1724 0 : FpXQE_log(GEN a, GEN b, GEN o, GEN a4, GEN T, GEN p)
1725 : {
1726 0 : pari_sp av = avma;
1727 : struct _FpXQE e;
1728 0 : e.a4=a4; e.T=T; e.p=p;
1729 0 : return gerepileuptoint(av, gen_PH_log(a, b, o, (void*)&e, &FpXQE_group));
1730 : }
1731 :
1732 : /***********************************************************************/
1733 : /** **/
1734 : /** Pairings **/
1735 : /** **/
1736 : /***********************************************************************/
1737 :
1738 : /* Derived from APIP from and by Jerome Milan, 2012 */
1739 :
1740 : static GEN
1741 4788 : FpXQE_vert(GEN P, GEN Q, GEN a4, GEN T, GEN p)
1742 : {
1743 4788 : long vT = get_FpX_var(T);
1744 4788 : if (ell_is_inf(P))
1745 70 : return pol_1(get_FpX_var(T));
1746 4718 : if (!ZX_equal(gel(Q, 1), gel(P, 1)))
1747 4718 : return FpX_sub(gel(Q, 1), gel(P, 1), p);
1748 0 : if (signe(gel(P,2))!=0) return pol_1(vT);
1749 0 : return FpXQ_inv(FpX_add(FpX_mulu(FpXQ_sqr(gel(P,1), T, p), 3, p),
1750 : a4, p), T, p);
1751 : }
1752 :
1753 : static GEN
1754 4655 : FpXQE_Miller_line(GEN R, GEN Q, GEN slope, GEN a4, GEN T, GEN p)
1755 : {
1756 4655 : long vT = get_FpX_var(T);
1757 4655 : GEN x = gel(Q, 1), y = gel(Q, 2);
1758 4655 : GEN tmp1 = FpX_sub(x, gel(R, 1), p);
1759 4655 : GEN tmp2 = FpX_add(FpXQ_mul(tmp1, slope, T, p), gel(R, 2), p);
1760 4655 : if (!ZX_equal(y, tmp2))
1761 4655 : return FpX_sub(y, tmp2, p);
1762 0 : if (signe(y) == 0)
1763 0 : return pol_1(vT);
1764 : else
1765 : {
1766 : GEN s1, s2;
1767 0 : GEN y2i = FpXQ_inv(FpX_mulu(y, 2, p), T, p);
1768 0 : s1 = FpXQ_mul(FpX_add(FpX_mulu(FpXQ_sqr(x, T, p), 3, p), a4, p), y2i, T, p);
1769 0 : if (!ZX_equal(s1, slope))
1770 0 : return FpX_sub(s1, slope, p);
1771 0 : s2 = FpXQ_mul(FpX_sub(FpX_mulu(x, 3, p), FpXQ_sqr(s1, T, p), p), y2i, T, p);
1772 0 : return signe(s2)!=0 ? s2: y2i;
1773 : }
1774 : }
1775 :
1776 : /* Computes the equation of the line tangent to R and returns its
1777 : evaluation at the point Q. Also doubles the point R.
1778 : */
1779 :
1780 : static GEN
1781 4158 : FpXQE_tangent_update(GEN R, GEN Q, GEN a4, GEN T, GEN p, GEN *pt_R)
1782 : {
1783 4158 : if (ell_is_inf(R))
1784 : {
1785 7 : *pt_R = ellinf();
1786 7 : return pol_1(get_FpX_var(T));
1787 : }
1788 4151 : else if (!signe(gel(R,2)))
1789 : {
1790 56 : *pt_R = ellinf();
1791 56 : return FpXQE_vert(R, Q, a4, T, p);
1792 : } else {
1793 : GEN slope;
1794 4095 : *pt_R = FpXQE_dbl_slope(R, a4, T, p, &slope);
1795 4095 : return FpXQE_Miller_line(R, Q, slope, a4, T, p);
1796 : }
1797 : }
1798 :
1799 : /* Computes the equation of the line through R and P, and returns its
1800 : evaluation at the point Q. Also adds P to the point R.
1801 : */
1802 :
1803 : static GEN
1804 567 : FpXQE_chord_update(GEN R, GEN P, GEN Q, GEN a4, GEN T, GEN p, GEN *pt_R)
1805 : {
1806 567 : if (ell_is_inf(R))
1807 : {
1808 0 : *pt_R = gcopy(P);
1809 0 : return FpXQE_vert(P, Q, a4, T, p);
1810 : }
1811 567 : else if (ell_is_inf(P))
1812 : {
1813 0 : *pt_R = gcopy(R);
1814 0 : return FpXQE_vert(R, Q, a4, T, p);
1815 : }
1816 567 : else if (ZX_equal(gel(P, 1), gel(R, 1)))
1817 : {
1818 7 : if (ZX_equal(gel(P, 2), gel(R, 2)))
1819 0 : return FpXQE_tangent_update(R, Q, a4, T, p, pt_R);
1820 : else
1821 : {
1822 7 : *pt_R = ellinf();
1823 7 : return FpXQE_vert(R, Q, a4, T, p);
1824 : }
1825 : } else {
1826 : GEN slope;
1827 560 : *pt_R = FpXQE_add_slope(P, R, a4, T, p, &slope);
1828 560 : return FpXQE_Miller_line(R, Q, slope, a4, T, p);
1829 : }
1830 : }
1831 :
1832 : struct _FpXQE_miller { GEN p, T, a4, P; };
1833 : static GEN
1834 4158 : FpXQE_Miller_dbl(void* E, GEN d)
1835 : {
1836 4158 : struct _FpXQE_miller *m = (struct _FpXQE_miller *)E;
1837 4158 : GEN p = m->p;
1838 4158 : GEN T = m->T, a4 = m->a4, P = m->P;
1839 : GEN v, line;
1840 4158 : GEN N = FpXQ_sqr(gel(d,1), T, p);
1841 4158 : GEN D = FpXQ_sqr(gel(d,2), T, p);
1842 4158 : GEN point = gel(d,3);
1843 4158 : line = FpXQE_tangent_update(point, P, a4, T, p, &point);
1844 4158 : N = FpXQ_mul(N, line, T, p);
1845 4158 : v = FpXQE_vert(point, P, a4, T, p);
1846 4158 : D = FpXQ_mul(D, v, T, p); return mkvec3(N, D, point);
1847 : }
1848 :
1849 : static GEN
1850 567 : FpXQE_Miller_add(void* E, GEN va, GEN vb)
1851 : {
1852 567 : struct _FpXQE_miller *m = (struct _FpXQE_miller *)E;
1853 567 : GEN p = m->p;
1854 567 : GEN T = m->T, a4 = m->a4, P = m->P;
1855 : GEN v, line, point;
1856 567 : GEN na = gel(va,1), da = gel(va,2), pa = gel(va,3);
1857 567 : GEN nb = gel(vb,1), db = gel(vb,2), pb = gel(vb,3);
1858 567 : GEN N = FpXQ_mul(na, nb, T, p);
1859 567 : GEN D = FpXQ_mul(da, db, T, p);
1860 567 : line = FpXQE_chord_update(pa, pb, P, a4, T, p, &point);
1861 567 : N = FpXQ_mul(N, line, T, p);
1862 567 : v = FpXQE_vert(point, P, a4, T, p);
1863 567 : D = FpXQ_mul(D, v, T, p); return mkvec3(N, D, point);
1864 : }
1865 :
1866 : /* Returns the Miller function f_{m, Q} evaluated at the point P using
1867 : * the standard Miller algorithm. */
1868 : static GEN
1869 63 : FpXQE_Miller(GEN Q, GEN P, GEN m, GEN a4, GEN T, GEN p)
1870 : {
1871 63 : pari_sp av = avma;
1872 : struct _FpXQE_miller d;
1873 : GEN v, N, D, g1;
1874 :
1875 63 : d.a4 = a4; d.T = T; d.p = p; d.P = P;
1876 63 : g1 = pol_1(get_FpX_var(T));
1877 63 : v = gen_pow_i(mkvec3(g1,g1,Q), m, (void*)&d,
1878 : FpXQE_Miller_dbl, FpXQE_Miller_add);
1879 63 : N = gel(v,1); D = gel(v,2);
1880 63 : return gerepileupto(av, FpXQ_div(N, D, T, p));
1881 : }
1882 :
1883 : GEN
1884 31 : FpXQE_weilpairing(GEN P, GEN Q, GEN m, GEN a4, GEN T, GEN p)
1885 : {
1886 31 : pari_sp av = avma;
1887 : GEN N, D, w;
1888 31 : if (ell_is_inf(P) || ell_is_inf(Q) || ZXV_equal(P,Q))
1889 3 : return pol_1(get_FpX_var(T));
1890 28 : N = FpXQE_Miller(P, Q, m, a4, T, p);
1891 28 : D = FpXQE_Miller(Q, P, m, a4, T, p);
1892 28 : w = FpXQ_div(N, D, T, p);
1893 28 : if (mpodd(m)) w = FpX_neg(w, p);
1894 28 : return gerepileupto(av, w);
1895 : }
1896 :
1897 : GEN
1898 7 : FpXQE_tatepairing(GEN P, GEN Q, GEN m, GEN a4, GEN T, GEN p)
1899 : {
1900 7 : if (ell_is_inf(P) || ell_is_inf(Q)) return pol_1(get_FpX_var(T));
1901 7 : return FpXQE_Miller(P, Q, m, a4, T, p);
1902 : }
1903 :
1904 : /***********************************************************************/
1905 : /** **/
1906 : /** issupersingular **/
1907 : /** **/
1908 : /***********************************************************************/
1909 :
1910 : GEN
1911 1718 : FpXQ_ellj(GEN a4, GEN a6, GEN T, GEN p)
1912 : {
1913 1718 : if (absequaliu(p,3)) return pol_0(get_FpX_var(T));
1914 : else
1915 : {
1916 1718 : pari_sp av=avma;
1917 1718 : GEN a43 = FpXQ_mul(a4,FpXQ_sqr(a4,T,p),T,p);
1918 1718 : GEN a62 = FpXQ_sqr(a6,T,p);
1919 1718 : GEN num = FpX_mulu(a43,6912,p);
1920 1718 : GEN den = FpX_add(FpX_mulu(a43,4,p),FpX_mulu(a62,27,p),p);
1921 1718 : return gerepileuptoleaf(av, FpXQ_div(num, den, T, p));
1922 : }
1923 : }
1924 :
1925 : static GEN
1926 33530 : FpXQ_is_quad(GEN x, GEN T, GEN p)
1927 : {
1928 33530 : pari_sp av = avma;
1929 : GEN K;
1930 33530 : long d = degpol(T);
1931 33530 : x = FpXQ_red(x,T,p);
1932 33530 : if (lgpol(x)<=1) return NULL;
1933 33530 : if (d==2) return FpXQ_minpoly(x, T, p);
1934 33530 : if (odd(degpol(T))) return NULL;
1935 33530 : K = FpM_ker(FpXQ_matrix_pow(x, d, 3, T, p), p);
1936 33530 : if (lg(K)!=2) return gc_NULL(av);
1937 588 : return RgV_to_RgX(gel(K,1), get_FpX_var(T));
1938 : }
1939 :
1940 : int
1941 165515 : FpXQ_elljissupersingular(GEN j, GEN T, GEN p)
1942 : {
1943 165515 : pari_sp ltop = avma;
1944 :
1945 : /* All supersingular j-invariants are in FF_{p^2}, so we first check
1946 : * whether j is in FF_{p^2}. If d is odd, then FF_{p^2} is not a
1947 : * subfield of FF_{p^d} so the j-invariants are all in FF_p. Hence
1948 : * the j-invariants are in FF_{p^{2 - e}}. */
1949 165515 : ulong d = get_FpX_degree(T);
1950 : GEN S;
1951 165515 : if (degpol(j) <= 0) return Fp_elljissupersingular(constant_coeff(j), p);
1952 164660 : j = FpXQ_red(j, T, p);
1953 164660 : if (degpol(j) <= 0) return gc_bool(ltop, Fp_elljissupersingular(constant_coeff(j), p));
1954 : /* Now j is not in F_p */
1955 164660 : if (abscmpiu(p, 5) <= 0) return gc_bool(ltop,0); /* j != 0*/
1956 164653 : if (odd(d)) return 0;
1957 : /* Set S so that FF_p[T]/(S) is isomorphic to FF_{p^2}: */
1958 46997 : if (d == 2)
1959 13467 : S = T;
1960 : else /* d > 2 */
1961 : {
1962 33530 : S = FpXQ_is_quad(j, T, p);
1963 33530 : if (!S) return gc_bool(ltop,0);
1964 588 : j = pol_x(varn(S));
1965 : }
1966 14055 : return gc_bool(ltop, jissupersingular(j,S,p));
1967 : }
1968 :
1969 : int
1970 1050 : Fq_elljissupersingular(GEN j, GEN T, GEN p)
1971 959 : { return typ(j)==t_INT? Fp_elljissupersingular(j, p)
1972 2009 : : FpXQ_elljissupersingular(j, T, p); }
1973 :
1974 : /* p > 5 prime; return d such that (-d/p) = -1 */
1975 : static ulong
1976 1183 : find_inert_disc(GEN p)
1977 : {
1978 1183 : long s = mod4(p) == 1? -1: 1; /* - (-1/p) */
1979 1183 : ulong d = 3;
1980 : while(1)
1981 : {
1982 1190 : if (kroui(d,p) == s) return d; /* = 3 mod (16) */
1983 595 : d++;
1984 595 : if (kroui(d>>2,p) == s) return d; /* = 4 mod (16) */
1985 266 : d += 3;
1986 266 : if (kroui(d,p) == s) return d; /* = 7 mod (16) */
1987 105 : d++;
1988 105 : if (kroui(d>>2,p) == s) return d; /* = 8 mod (16) */
1989 35 : d += 3;
1990 35 : if (kroui(d,p) == s) return d; /* = 11 mod (16) */
1991 7 : d += 4;
1992 7 : if (kroui(d,p) == s) return d; /* = 15 mod (16) */
1993 7 : d += 4;
1994 : }
1995 : }
1996 :
1997 : /* p > 5 */
1998 : static GEN
1999 1183 : ellsupersingularj_easy_FpXQ(GEN T, GEN p)
2000 : {
2001 1183 : long d = find_inert_disc(p);
2002 1183 : GEN R = FpXQX_roots(polclass(stoi(-d), 0, 0), T, p);
2003 1183 : return gel(R,1);
2004 : }
2005 :
2006 : GEN
2007 1204 : ellsupersingularj_FpXQ(GEN T, GEN p)
2008 : {
2009 : GEN j, j2, R, Phi2;
2010 : long i, ep, lp;
2011 1204 : if (cmpiu(p, 5) <= 0) return pol_0(get_FpX_var(T));
2012 1183 : j2 = ellsupersingularj_easy_FpXQ(T, p);
2013 1183 : Phi2 = polmodular_ZXX(2,0,0,1);
2014 1183 : R = FpXQX_roots(FqXY_evalx(Phi2, j2, T, p), T, p);
2015 1183 : j = gel(R,1+random_Fl(lg(R)-1));
2016 1183 : ep = expi(p); lp = ep + random_Fl(ep);
2017 18713 : for (i = 1; i <= lp; i++)
2018 : {
2019 17530 : GEN Phi2_j = FqX_div_by_X_x(FqXY_evalx(Phi2, j, T, p), j2, T, p, NULL);
2020 17530 : R = FqX_quad_root(Phi2_j, T, p);
2021 17530 : if (!R) pari_err_PRIME("ellsupersingularj",p);
2022 17530 : j2 = j; j = random_bits(1) ? R: Fq_neg(Fq_add(gel(Phi2_j,3), R, T, p), T, p);
2023 : }
2024 1183 : return j;
2025 : }
2026 :
2027 : /***********************************************************************/
2028 : /** **/
2029 : /** Point counting **/
2030 : /** **/
2031 : /***********************************************************************/
2032 :
2033 : GEN
2034 15540 : elltrace_extension(GEN t, long n, GEN q)
2035 : {
2036 15540 : pari_sp av = avma;
2037 15540 : GEN v = RgX_to_RgC(RgXQ_powu(pol_x(0), n, mkpoln(3,gen_1,negi(t),q)),2);
2038 15540 : GEN te = addii(shifti(gel(v,1),1), mulii(t,gel(v,2)));
2039 15540 : return gerepileuptoint(av, te);
2040 : }
2041 :
2042 : GEN
2043 14777 : Fp_ffellcard(GEN a4, GEN a6, GEN q, long n, GEN p)
2044 : {
2045 14777 : pari_sp av = avma;
2046 14777 : GEN ap = subii(addiu(p, 1), Fp_ellcard(a4, a6, p));
2047 14777 : GEN te = elltrace_extension(ap, n, p);
2048 14777 : return gerepileuptoint(av, subii(addiu(q, 1), te));
2049 : }
2050 :
2051 : static GEN
2052 1687 : FpXQ_ellcardj(GEN a4, GEN a6, GEN j, GEN T, GEN q, GEN p, long n)
2053 : {
2054 1687 : GEN q1 = addiu(q,1);
2055 1687 : if (signe(j)==0)
2056 : {
2057 : GEN W, w, t, N;
2058 560 : if (umodiu(q,6)!=1) return q1;
2059 420 : N = Fp_ffellcard(gen_0,gen_1,q,n,p);
2060 420 : t = subii(q1, N);
2061 420 : W = FpXQ_pow(a6,diviuexact(shifti(q,-1), 3),T,p);
2062 420 : if (degpol(W)>0) /*p=5 mod 6*/
2063 105 : return ZX_equal1(FpXQ_powu(W,3,T,p)) ? addii(q1,shifti(t,-1)):
2064 35 : subii(q1,shifti(t,-1));
2065 350 : w = modii(gel(W,2),p);
2066 350 : if (equali1(w)) return N;
2067 259 : if (equalii(w,subiu(p,1))) return addii(q1,t);
2068 : else /*p=1 mod 6*/
2069 : {
2070 168 : GEN u = shifti(t,-1), v = sqrtint(diviuexact(subii(q,sqri(u)),3));
2071 168 : GEN a = addii(u,v), b = shifti(v,1);
2072 168 : if (equali1(Fp_powu(w,3,p)))
2073 : {
2074 84 : if (dvdii(addmulii(a, w, b), p))
2075 56 : return subii(q1,subii(shifti(b,1),a));
2076 : else
2077 28 : return addii(q1,addii(a,b));
2078 : }
2079 : else
2080 : {
2081 84 : if (dvdii(submulii(a, w, b), p))
2082 56 : return subii(q1,subii(a,shifti(b,1)));
2083 : else
2084 28 : return subii(q1,addii(a,b));
2085 : }
2086 : }
2087 1127 : } else if (equalii(j,modsi(1728,p)))
2088 : {
2089 : GEN w, W, N, t;
2090 567 : if (mod4(q)==3) return q1;
2091 427 : W = FpXQ_pow(a4,shifti(q,-2),T,p);
2092 427 : if (degpol(W)>0) return q1; /*p=3 mod 4*/
2093 357 : w = modii(gel(W,2),p);
2094 357 : N = Fp_ffellcard(gen_1,gen_0,q,n,p);
2095 357 : if (equali1(w)) return N;
2096 238 : t = subii(q1, N);
2097 238 : if (equalii(w,subiu(p,1))) return addii(q1,t);
2098 : else /*p=1 mod 4*/
2099 : {
2100 112 : GEN u = shifti(t,-1), v = sqrtint(subii(q,sqri(u)));
2101 112 : if (dvdii(addmulii(u, w, v), p))
2102 56 : return subii(q1,shifti(v,1));
2103 : else
2104 56 : return addii(q1,shifti(v,1));
2105 : }
2106 : } else
2107 : {
2108 560 : GEN g = Fp_div(j, Fp_sub(utoi(1728), j, p), p);
2109 560 : GEN l = FpXQ_div(FpX_mulu(a6,3,p),FpX_mulu(a4,2,p),T,p);
2110 560 : GEN N = Fp_ffellcard(Fp_mulu(g,3,p),Fp_double(g,p),q,n,p);
2111 560 : if (FpXQ_issquare(l,T,p)) return N;
2112 280 : return subii(shifti(q1,1),N);
2113 : }
2114 : }
2115 :
2116 : static GEN
2117 8 : FpXQ_ffellcard(GEN a4, GEN a6, GEN M, GEN q, GEN T, GEN p, long n)
2118 : {
2119 8 : long m = degpol(M);
2120 8 : GEN j = pol_x(get_FpX_var(T));
2121 8 : GEN g = FpXQ_div(j, Fp_FpX_sub(utoi(1728), j, p), M, p);
2122 8 : GEN N = FpXQ_ellcard(FpX_mulu(g,3,p),FpX_mulu(g,2,p),M,p);
2123 8 : GEN qm = powiu(p, m), q1 = addiu(q, 1), qm1 = addiu(qm, 1);
2124 8 : GEN l = FpXQ_mul(FpX_mulu(a6,3,p),FpX_mulu(a4,2,p),T,p);
2125 8 : GEN te = elltrace_extension(subii(qm1, N), n/m, qm);
2126 8 : return FpXQ_issquare(l,T,p) ? subii(q1, te): addii(q1, te);
2127 : }
2128 :
2129 : static int
2130 7 : FpXQ_is4power(GEN x, GEN T, GEN p)
2131 : {
2132 7 : long d = get_FpX_degree(T);
2133 7 : if (lg(x) == 2 || absequalui(2, p)) return 1;
2134 7 : if (Mod4(p)==1)
2135 7 : return equali1(Fp_pow(FpXQ_norm(x,T,p),shifti(p,-2), p));
2136 0 : if (odd(d))
2137 0 : return FpXQ_issquare(x, T, p);
2138 0 : return ZX_equal1(FpXQ_pow(x, shifti(powiu(p, d),-2), T, p));
2139 : }
2140 :
2141 : /* http://www.numdam.org/article/ASENS_1969_4_2_4_521_0.pdf */
2142 :
2143 : GEN
2144 7 : FpXQ_ellcard_supersingular(GEN a4, GEN a6, GEN T, GEN p)
2145 : {
2146 7 : pari_sp av = avma;
2147 7 : long d = get_FpX_degree(T);
2148 : GEN r;
2149 7 : if (equaliu(p,3))
2150 0 : r = Flxq_ellcard(ZX_to_Flx(a4,3), ZX_to_Flx(a6,3), ZXT_to_FlxT(T,3), 3);
2151 7 : else if (signe(a4)==0)
2152 0 : r = FpXQ_ellcardj(a4, a6, gen_0, T, powiu(p, d), p, d);
2153 7 : else if (signe(a6)==0)
2154 0 : r = FpXQ_ellcardj(a4, a6, modsi(1728,p), T, powiu(p, d), p, d);
2155 : else
2156 : {
2157 : GEN q, q2, t, D;
2158 7 : long qm4 = (odd(d>>1) && Mod4(p)==3);
2159 7 : if (odd(d)) return gen_0;
2160 7 : q2 = powiu(p, d>>1); q = sqri(q2);
2161 7 : t = shifti(q2, 1);
2162 7 : D = FpX_sub(FpX_Fp_mul(FpXQ_powu(a4,3,T,p), stoi(-4), p),
2163 : FpX_mulu(FpXQ_sqr(a6,T,p), 27, p), p);
2164 14 : r = qm4 ^ FpXQ_is4power(D, T, p) ? subii(addiu(q, 1), t)
2165 7 : : addii(addiu(q, 1), t);
2166 : }
2167 7 : return gerepileuptoint(av, r);
2168 : }
2169 :
2170 : GEN
2171 21 : Fq_ellcard_supersingular(GEN a4, GEN a6, GEN T, GEN p)
2172 21 : { return T ? FpXQ_ellcard_supersingular(a4, a6, T, p) : addiu(p, 1); }
2173 :
2174 : static GEN
2175 8571 : FpXQ_ellcard_i(GEN a4, GEN a6, GEN T, GEN p)
2176 : {
2177 8571 : long n = get_FpX_degree(T);
2178 8571 : GEN q = powiu(p, n);
2179 8571 : if (degpol(a4)<=0 && degpol(a6)<=0)
2180 833 : return Fp_ffellcard(constant_coeff(a4),constant_coeff(a6),q,n,p);
2181 7738 : if (lgefint(p)==3)
2182 : {
2183 6020 : ulong pp = p[2];
2184 6020 : return Flxq_ellcard(ZX_to_Flx(a4,pp),ZX_to_Flx(a6,pp),ZX_to_Flx(T,pp),pp);
2185 : }
2186 : else
2187 : {
2188 1718 : GEN J = FpXQ_ellj(a4,a6,T,p), M;
2189 1718 : if (degpol(J) <= 0)
2190 1687 : return FpXQ_ellcardj(a4,a6,constant_coeff(J),T,q,p,n);
2191 31 : M = FpXQ_minpoly(J,T,p);
2192 31 : if (degpol(M) < degpol(T))
2193 8 : return FpXQ_ffellcard(a4, a6, M, q, T, p, n);
2194 23 : return Fq_ellcard_SEA(a4, a6, q, T, p, 0);
2195 : }
2196 : }
2197 :
2198 : GEN
2199 8571 : FpXQ_ellcard(GEN a4, GEN a6, GEN T, GEN p)
2200 : {
2201 8571 : pari_sp av = avma;
2202 8571 : return gerepileuptoint(av, FpXQ_ellcard_i(a4, a6, T, p));
2203 : }
2204 :
2205 : static GEN
2206 24 : _FpXQE_pairorder(void *E, GEN P, GEN Q, GEN m, GEN F)
2207 : {
2208 24 : struct _FpXQE *e = (struct _FpXQE *) E;
2209 24 : return FpXQ_order(FpXQE_weilpairing(P,Q,m,e->a4,e->T,e->p), F, e->T, e->p);
2210 : }
2211 :
2212 : GEN
2213 15 : FpXQ_ellgroup(GEN a4, GEN a6, GEN N, GEN T, GEN p, GEN *pt_m)
2214 : {
2215 : struct _FpXQE e;
2216 15 : GEN q = powiu(p, get_FpX_degree(T));
2217 15 : e.a4=a4; e.a6=a6; e.T=T; e.p=p;
2218 15 : return gen_ellgroup(N, subiu(q,1), pt_m, (void*)&e, &FpXQE_group, _FpXQE_pairorder);
2219 : }
2220 :
2221 : GEN
2222 8 : FpXQ_ellgens(GEN a4, GEN a6, GEN ch, GEN D, GEN m, GEN T, GEN p)
2223 : {
2224 : GEN P;
2225 8 : pari_sp av = avma;
2226 : struct _FpXQE e;
2227 8 : e.a4=a4; e.a6=a6; e.T=T; e.p=p;
2228 8 : switch(lg(D)-1)
2229 : {
2230 8 : case 1:
2231 8 : P = gen_gener(gel(D,1), (void*)&e, &FpXQE_group);
2232 8 : P = mkvec(FpXQE_changepoint(P, ch, T, p));
2233 8 : break;
2234 0 : default:
2235 0 : P = gen_ellgens(gel(D,1), gel(D,2), m, (void*)&e, &FpXQE_group, _FpXQE_pairorder);
2236 0 : gel(P,1) = FpXQE_changepoint(gel(P,1), ch, T, p);
2237 0 : gel(P,2) = FpXQE_changepoint(gel(P,2), ch, T, p);
2238 0 : break;
2239 : }
2240 8 : return gerepilecopy(av, P);
2241 : }
|
