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 49719 : ellinf_FpJ(void)
35 49719 : { 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 6384496 : 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 6384496 : if (signe(gel(P,3)) == 0) return ellinf_FpJ();
46 :
47 6375794 : X1 = gel(P,1); Y1 = gel(P,2); Z1 = gel(P,3);
48 :
49 6375794 : XX = Fp_sqr(X1, p);
50 6400202 : YY = Fp_sqr(Y1, p);
51 6400242 : YYYY = Fp_sqr(YY, p);
52 6396446 : ZZ = Fp_sqr(Z1, p);
53 6398107 : S = Fp_double(Fp_sub(Fp_sqr(Fp_add(X1,YY,p), p), Fp_add(XX,YYYY,p), p), p);
54 6346828 : M = Fp_addmul(Fp_mulu(XX, 3, p), a4, Fp_sqr(ZZ, p), p);
55 6384308 : T = Fp_sub(Fp_sqr(M, p), Fp_double(S, p), p);
56 6363598 : Q = cgetg(4, t_VEC);
57 6367748 : gel(Q,1) = T;
58 6367748 : gel(Q,2) = Fp_sub(Fp_mul(M, Fp_sub(S, T, p), p), Fp_mulu(YYYY, 8, p), p);
59 6364606 : gel(Q,3) = Fp_sub(Fp_sqr(Fp_add(Y1, Z1, p), p), Fp_add(YY, ZZ, p), p);
60 6362344 : 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 1150159 : 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 1150159 : if (signe(gel(Q,3)) == 0) return gcopy(P);
72 1150159 : if (signe(gel(P,3)) == 0) return gcopy(Q);
73 :
74 1130033 : X1 = gel(P,1); Y1 = gel(P,2); Z1 = gel(P,3);
75 1130033 : X2 = gel(Q,1); Y2 = gel(Q,2); Z2 = gel(Q,3);
76 :
77 1130033 : Z1Z1 = Fp_sqr(Z1, p);
78 1130452 : Z2Z2 = Fp_sqr(Z2, p);
79 1130287 : U1 = Fp_mul(X1, Z2Z2, p);
80 1130376 : U2 = Fp_mul(X2, Z1Z1, p);
81 1130382 : S1 = mulii(Y1, Fp_mul(Z2, Z2Z2, p));
82 1129647 : S2 = mulii(Y2, Fp_mul(Z1, Z1Z1, p));
83 1129703 : H = Fp_sub(U2, U1, p);
84 1129767 : r = Fp_double(Fp_sub(S2, S1, p), p);
85 :
86 : /* If points are equal we must double. */
87 1129485 : if (signe(H)== 0) {
88 42060 : if (signe(r) == 0)
89 : /* Points are equal so double. */
90 1043 : return FpJ_dbl(P, a4, p);
91 : else
92 41017 : return ellinf_FpJ();
93 : }
94 1087425 : I = Fp_sqr(Fp_double(H, p), p);
95 1088420 : J = Fp_mul(H, I, p);
96 1088338 : V = Fp_mul(U1, I, p);
97 1088284 : W = Fp_sub(Fp_sqr(r, p), Fp_add(J, Fp_double(V, p), p), p);
98 1087680 : R = cgetg(4, t_VEC);
99 1087772 : gel(R,1) = W;
100 1087772 : gel(R,2) = Fp_sub(mulii(r, subii(V, W)),
101 : shifti(mulii(S1, J), 1), p);
102 1087884 : 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 1088168 : 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 204042 : FpE_to_FpJ(GEN P)
115 : {
116 204042 : return ell_is_inf(P) ? ellinf_FpJ()
117 204042 : : mkvec3(icopy(gel(P,1)),icopy(gel(P,2)), gen_1);
118 : }
119 :
120 : GEN
121 203539 : FpJ_to_FpE(GEN P, GEN p)
122 : {
123 203539 : if (signe(gel(P,3)) == 0) return ellinf();
124 : else
125 : {
126 162886 : GEN Z = Fp_inv(gel(P,3), p);
127 162860 : GEN Z2 = Fp_sqr(Z, p), Z3 = Fp_mul(Z, Z2, p);
128 162860 : 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 6385766 : _FpJ_dbl(void *E, GEN P)
135 : {
136 6385766 : struct _FpE *ell = (struct _FpE *) E;
137 6385766 : return FpJ_dbl(P, ell->a4, ell->p);
138 : }
139 : static GEN
140 1150009 : _FpJ_add(void *E, GEN P, GEN Q)
141 : {
142 1150009 : struct _FpE *ell=(struct _FpE *) E;
143 1150009 : return FpJ_add(P, Q, ell->a4, ell->p);
144 : }
145 : static GEN
146 5712 : _FpJ_mul(void *E, GEN P, GEN n)
147 : {
148 5712 : pari_sp av = avma;
149 5712 : struct _FpE *e=(struct _FpE *) E;
150 5712 : long s = signe(n);
151 5712 : if (!s || signe(gel(P,3))==0) return ellinf_FpJ();
152 5712 : if (s < 0) P = FpJ_neg(P, e->p);
153 5712 : if (is_pm1(n)) return s > 0 ? gcopy(P): P;
154 5711 : return gc_GEN(av, gen_pow_i(P, n, e, &_FpJ_dbl, &_FpJ_add));
155 : }
156 :
157 : GEN
158 5712 : FpJ_mul(GEN P, GEN n, GEN a4, GEN p)
159 : {
160 : struct _FpE E;
161 5712 : E.a4= a4; E.p = p;
162 5712 : 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 gc_upto(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 948 : 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 948 : return gc_upto(av, z);
208 : }
209 :
210 : GEN
211 2736 : FpE_changepointinv(GEN P, GEN ch, GEN p)
212 : {
213 : GEN u, r, s, t, u2, u3, c, z;
214 2736 : if (ell_is_inf(P)) return P;
215 2736 : 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 998 : u = gel(ch,1); r = gel(ch,2); s = gel(ch,3); t = gel(ch,4);
222 998 : u2 = Fp_sqr(u, p); u3 = Fp_mul(u,u2,p);
223 998 : c = Fp_mul(u2, gel(P,1), p);
224 998 : z = cgetg(3, t_VEC);
225 998 : gel(z,1) = Fp_add(c,r,p);
226 998 : gel(z,2) = Fp_add(Fp_mul(u3,gel(P,2),p), Fp_add(Fp_mul(s,c,p), t, p), p);
227 997 : 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 288746 : FpE_dbl_slope(GEN P, GEN a4, GEN p, GEN *slope)
258 : {
259 : GEN x, y, Q;
260 288746 : if (ell_is_inf(P) || !signe(gel(P,2))) return ellinf();
261 132881 : x = gel(P,1); y = gel(P,2);
262 132881 : *slope = Fp_div(Fp_add(Fp_mulu(Fp_sqr(x,p), 3, p), a4, p),
263 : Fp_mulu(y, 2, p), p);
264 132881 : Q = cgetg(3,t_VEC);
265 132881 : gel(Q, 1) = Fp_sub(Fp_sqr(*slope, p), Fp_mulu(x, 2, p), p);
266 132881 : gel(Q, 2) = Fp_sub(Fp_mul(*slope, Fp_sub(x, gel(Q, 1), p), p), y, p);
267 132881 : return Q;
268 : }
269 :
270 : GEN
271 288152 : FpE_dbl(GEN P, GEN a4, GEN p)
272 : {
273 288152 : pari_sp av = avma;
274 : GEN slope;
275 288152 : return gc_upto(av, FpE_dbl_slope(P,a4,p,&slope));
276 : }
277 :
278 : static GEN
279 916619 : FpE_add_slope(GEN P, GEN Q, GEN a4, GEN p, GEN *slope)
280 : {
281 : GEN Px, Py, Qx, Qy, R;
282 916619 : if (ell_is_inf(P)) return Q;
283 916129 : if (ell_is_inf(Q)) return P;
284 916129 : Px = gel(P,1); Py = gel(P,2);
285 916129 : Qx = gel(Q,1); Qy = gel(Q,2);
286 916129 : 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 915555 : *slope = Fp_div(Fp_sub(Py, Qy, p), Fp_sub(Px, Qx, p), p);
294 915555 : R = cgetg(3,t_VEC);
295 915555 : gel(R, 1) = Fp_sub(Fp_sub(Fp_sqr(*slope, p), Px, p), Qx, p);
296 915555 : gel(R, 2) = Fp_sub(Fp_mul(*slope, Fp_sub(Px, gel(R, 1), p), p), Py, p);
297 915555 : return R;
298 : }
299 :
300 : GEN
301 916615 : FpE_add(GEN P, GEN Q, GEN a4, GEN p)
302 : {
303 916615 : pari_sp av = avma;
304 : GEN slope;
305 916615 : return gc_upto(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 362490 : FpE_neg(GEN P, GEN p)
317 : {
318 362490 : if (ell_is_inf(P)) return ellinf();
319 362490 : 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 gc_upto(av, FpE_add_slope(P, FpE_neg_i(Q, p), a4, p, &slope));
328 : }
329 :
330 : static GEN
331 288152 : _FpE_dbl(void *E, GEN P)
332 : {
333 288152 : struct _FpE *ell = (struct _FpE *) E;
334 288152 : 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 925527 : _FpE_mul(void *E, GEN P, GEN n)
346 : {
347 925527 : pari_sp av = avma;
348 925527 : struct _FpE *e=(struct _FpE *) E;
349 925527 : long s = signe(n);
350 : GEN Q;
351 925527 : if (!s || ell_is_inf(P)) return ellinf();
352 925520 : if (s<0) P = FpE_neg(P, e->p);
353 925520 : if (is_pm1(n)) return s>0? gcopy(P): P;
354 491718 : if (equalis(n,2)) return _FpE_dbl(E, P);
355 203566 : Q = gen_pow_i(FpE_to_FpJ(P), n, e, &_FpJ_dbl, &_FpJ_add);
356 203539 : return gc_upto(av, FpJ_to_FpE(Q, e->p));
357 : }
358 :
359 : GEN
360 1318 : FpE_mul(GEN P, GEN n, GEN a4, GEN p)
361 : {
362 : struct _FpE E;
363 1318 : E.a4 = a4; E.p = p;
364 1318 : return _FpE_mul(&E, P, n);
365 : }
366 :
367 : /* Finds a random nonsingular point on E */
368 :
369 : GEN
370 188685 : random_FpE(GEN a4, GEN a6, GEN p)
371 : {
372 188685 : pari_sp ltop = avma;
373 : GEN x, x2, y, rhs;
374 : do
375 : {
376 329033 : set_avma(ltop);
377 329033 : x = randomi(p); /* x^3+a4*x+a6 = x*(x^2+a4)+a6 */
378 329033 : x2 = Fp_sqr(x, p);
379 329033 : rhs = Fp_add(Fp_mul(x, Fp_add(x2, a4, p), p), a6, p);
380 35146 : } while ((!signe(rhs) && !signe(Fp_add(Fp_mulu(x2,3,p),a4,p)))
381 364179 : || kronecker(rhs, p) < 0);
382 188685 : y = Fp_sqrt(rhs, p);
383 188685 : if (!y) pari_err_PRIME("random_FpE", p);
384 188685 : return gc_GEN(ltop, mkvec2(x, y));
385 : }
386 :
387 : static GEN
388 186277 : _FpE_rand(void *E)
389 : {
390 186277 : struct _FpE *e=(struct _FpE *) E;
391 186277 : 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 737 : FpE_order(GEN z, GEN o, GEN a4, GEN p)
407 : {
408 737 : pari_sp av = avma;
409 : struct _FpE e;
410 : GEN r;
411 737 : if (lgefint(p) == 3)
412 : {
413 631 : ulong pp = p[2];
414 631 : 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 737 : return gc_INT(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 gc_INT(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 gc_INT(av, Fp_div(N, D, p));
586 : }
587 :
588 : GEN
589 72970 : FpE_weilpairing(GEN P, GEN Q, GEN m, GEN a4, GEN p)
590 : {
591 72970 : pari_sp av = avma;
592 : GEN N, D, w;
593 72970 : if (ell_is_inf(P) || ell_is_inf(Q) || ZV_equal(P,Q)) return gen_1;
594 48322 : if (lgefint(p)==3 && lgefint(m)==3)
595 : {
596 48300 : ulong pp = p[2];
597 48300 : GEN Pp = ZV_to_Flv(P, pp), Qp = ZV_to_Flv(Q, pp);
598 48300 : ulong w = Fle_weilpairing(Pp, Qp, itou(m), umodiu(a4, pp), pp);
599 48300 : 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 gc_INT(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 654049 : is_CMj(long J, GEN jn, GEN jd, GEN p)
632 654049 : { 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 14425 : u2_is_CMj(ulong a, ulong b, GEN jn, GEN jd, GEN p)
637 : {
638 14425 : GEN mJ = uu32toi(a,b);
639 14425 : return dvdii(addii(mulii(jd,mJ), jn), p);
640 : }
641 : #endif
642 :
643 : static long
644 52672 : 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 52672 : CHECK(-3, 0);
648 52564 : CHECK(-4, 1728);
649 52445 : CHECK(-7, -3375);
650 52199 : CHECK(-8, 8000);
651 51991 : CHECK(-11, -32768);
652 51763 : CHECK(-12, 54000);
653 51534 : CHECK(-16, 287496);
654 51357 : CHECK(-19, -884736);
655 51135 : CHECK(-27, -12288000);
656 50911 : CHECK(-28, 16581375);
657 50728 : CHECK(-43, -884736000);
658 : #ifdef LONG_IS_64BIT
659 43320 : CHECK(-67, -147197952000L);
660 43192 : CHECK(-163, -262537412640768000L);
661 : #else
662 7223 : if (u2_is_CMj(0x00000022UL,0x45ae8000UL,jn,jd,p)) return -67;
663 7202 : if (u2_is_CMj(0x03a4b862UL,0xc4b40000UL,jn,jd,p)) return -163;
664 : #endif
665 : #undef CHECK
666 50215 : 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 71283 : FqX_quad_root(GEN x, GEN T, GEN p)
678 : {
679 71283 : GEN b = gel(x,3), c = gel(x,2);
680 71283 : GEN D = Fq_sub(Fq_sqr(b, T, p), Fq_mulu(c,4, T, p), T, p);
681 71283 : GEN s = Fq_sqrt(D,T, p);
682 71283 : if (!s) return NULL;
683 68105 : return Fq_halve(Fq_sub(s, b, T, p), T, p);
684 : }
685 :
686 : static GEN
687 1230 : FpX_quad_root(GEN x, GEN p)
688 : {
689 1230 : GEN s, b = gel(x,3), c = gel(x,2);
690 1230 : GEN D = Fp_sub(Fp_sqr(b, p), shifti(c,2), p);
691 1230 : if (kronecker(D,p) == -1) return NULL;
692 782 : s = Fp_sqrt(D,p);
693 782 : 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 4860 : Fq_path_extends_to_floor(GEN j_prev, GEN j, GEN T, GEN p, GEN Phi2, long max_len)
701 : {
702 4860 : pari_sp ltop = avma;
703 4860 : 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 31987 : for (d = 1; d <= max_len; ++d)
709 : {
710 80090 : for (i = 1; i < l; i++)
711 : {
712 52963 : GEN Phi2_j = FqX_div_by_X_x(FqXY_evalx(Phi2, gel(j,i), T, p), gel(j_prev,i), T, p, NULL);
713 52963 : GEN j_next = FqX_quad_root(Phi2_j, T, p);
714 52963 : if (!j_next)
715 3178 : return gc_long(ltop, 1);
716 49785 : gel(j_prev,i) = gel(j, i); gel(j,i) = j_next;
717 : }
718 27127 : if (gc_needed(ltop, 2))
719 0 : (void)gc_all(ltop, 2, &j, &j_prev);
720 : }
721 1682 : 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 615 : for (d = 1; d <= max_len; ++d)
734 : {
735 1397 : for (i = 1; i < l; i++)
736 : {
737 1230 : GEN Phi2_j = FpX_div_by_X_x(FpXY_evalx(Phi2, gel(j,i), p), gel(j_prev,i), p, NULL);
738 1230 : GEN j_next = FpX_quad_root(Phi2_j, p);
739 1230 : 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 782 : gel(j_prev,i) = gel(j, i); gel(j,i) = j_next;
746 : }
747 167 : if (gc_needed(ltop, 2))
748 0 : (void)gc_all(ltop, 2, &j, &j_prev);
749 : }
750 0 : return gc_long(ltop, 0);
751 : }
752 :
753 :
754 : static int
755 2767 : Fp_jissupersingular(GEN j, GEN p)
756 : {
757 2767 : long max_path_len = expi(p)+1;
758 2767 : GEN Phi2 = FpXX_red(polmodular_ZXX(2,0,0,1), p);
759 2767 : GEN Phi2_j = FpXY_evalx(Phi2, j, p);
760 2767 : GEN roots = FpX_roots(Phi2_j, p);
761 2767 : long nbroots = lg(roots)-1;
762 2767 : 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 2767 : if (nbroots==0)
767 665 : return 0;
768 2102 : S = deg2pol_shallow(gen_1, gen_0, Fp_neg(Fp_2gener(p),p),1);
769 2102 : if (nbroots==1 && FpX_is_squarefree(Phi2_j, p))
770 1654 : { 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 2102 : return !Fq_path_extends_to_floor(mkvec(j_prev), mkvec(j), S, p, Phi2, max_path_len);
775 : }
776 :
777 : static int
778 14007 : jissupersingular(GEN j, GEN S, GEN p)
779 : {
780 14007 : long max_path_len = expi(p)+1;
781 14007 : GEN Phi2 = FpXX_red(polmodular_ZXX(2,0,0,1), p);
782 14007 : GEN Phi2_j = FqXY_evalx(Phi2, j, S, p);
783 14007 : GEN roots = FpXQX_roots(Phi2_j, S, p);
784 14007 : 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 14007 : if (nbroots==0 || (nbroots==1 && FqX_is_squarefree(Phi2_j, S, p)))
789 11249 : return 0;
790 : else
791 2758 : return !Fq_path_extends_to_floor(const_vec(nbroots,j), roots, S, p, Phi2, max_path_len);
792 : }
793 :
794 : int
795 3759 : Fp_elljissupersingular(GEN j, GEN p)
796 : {
797 : long CM;
798 3759 : if (abscmpiu(p, 5) <= 0) return signe(j) == 0; /* valid if p <= 5 */
799 3619 : CM = Fp_ellj_get_CM(j, gen_1, p);
800 3619 : if (CM < 0) return krosi(CM, p) < 0; /* valid if p > 3 */
801 : else
802 2767 : 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 723442 : Fl_elltrace_naive(ulong a4, ulong a6, ulong p)
814 : {
815 : ulong i, j;
816 723442 : long a = 0;
817 : long d0, d1, d2, d3;
818 723442 : GEN k = const_vecsmall(p, -1);
819 723469 : k[1] = 0;
820 130148107 : for (i=1, j=1; i < p; i += 2, j = Fl_add(j, i, p))
821 129424673 : k[j+1] = 1;
822 723434 : d0 = 6%p; d1 = d0; d2 = Fl_add(a4, 1, p); d3 = a6;
823 723434 : for(i=0;; i++)
824 : {
825 255764440 : a -= k[1+d3];
826 255764440 : if (i==p-1) break;
827 255041026 : d3 = Fl_add(d3, d2, p);
828 254927227 : d2 = Fl_add(d2, d1, p);
829 254886349 : d1 = Fl_add(d1, d0, p);
830 : }
831 723414 : 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 19196 : _fix(GEN x, long k)
857 : {
858 19196 : GEN y = (GEN)*x;
859 19196 : if (lgefint(y) < k) { GEN p1 = cgeti(k); affii(y,p1); *x = (long)p1; }
860 19196 : }
861 :
862 : /* Return the lift of a (mod b), which is closest to c */
863 : static GEN
864 254895 : closest_lift(GEN a, GEN b, GEN c)
865 : {
866 254895 : return addii(a, mulii(b, diviiround(subii(c,a), b)));
867 : }
868 :
869 : static long
870 79 : get_table_size(GEN pordmin, GEN B)
871 : {
872 79 : pari_sp av = avma;
873 79 : GEN t = ceilr( sqrtr( divri(itor(pordmin, DEFAULTPREC), B) ) );
874 79 : if (is_bigint(t))
875 0 : pari_err_OVERFLOW("ellap [large prime: install the 'seadata' package]");
876 79 : set_avma(av);
877 79 : 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 14220 : u = Fl_add(c6, Fl_mul(x, t, p), p);
905 14220 : 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 : /* y <- x, both are pairs of t_INT */
912 : static void
913 9440 : affii2(GEN x, GEN y)
914 : {
915 9440 : affii(gel(x,1), gel(y,1));
916 9440 : affii(gel(x,2), gel(y,2));
917 9440 : }
918 :
919 : static GEN ap_j1728(GEN a4,GEN p);
920 : /* compute a_p using Shanks/Mestre + Montgomery's trick. Assume p > 457 */
921 : static GEN
922 79 : Fp_ellcard_Shanks(GEN c4, GEN c6, GEN p)
923 : {
924 : pari_timer T;
925 : long *tx, *ty, *ti, pfinal, i, j, s, KRO, nb;
926 : ulong x;
927 79 : pari_sp av = avma, av2;
928 : GEN p1, P, mfh, h, F,f, fh,fg, pordmin, u, v, p1p, p2p, A, B, a4, pts;
929 79 : tx = NULL;
930 79 : ty = ti = NULL; /* gcc -Wall */
931 :
932 79 : if (!signe(c6)) {
933 0 : GEN ap = ap_j1728(c4, p);
934 0 : return gc_INT(av, subii(addiu(p,1), ap));
935 : }
936 :
937 79 : if (DEBUGLEVEL >= 6) timer_start(&T);
938 : /* once #E(Fp) is know mod B >= pordmin, it is completely determined */
939 79 : pordmin = addiu(sqrti(gmul2n(p,4)), 1); /* ceil( 4sqrt(p) ) */
940 79 : p1p = addiu(p, 1);
941 79 : p2p = shifti(p1p, 1);
942 79 : x = 0; KRO = 0;
943 : /* how many 2-torsion points ? */
944 79 : switch(FpX_nbroots(mkpoln(4, gen_1, gen_0, c4, c6), p))
945 : {
946 9 : case 3: A = gen_0; B = utoipos(4); break;
947 32 : case 1: A = gen_0; B = gen_2; break;
948 38 : default: A = gen_1; B = gen_2; break; /* 0 */
949 : }
950 : for(;;)
951 : {
952 79 : h = closest_lift(A, B, p1p);
953 79 : if (!KRO) /* first time, initialize */
954 : {
955 79 : KRO = kronecker(c6,p);
956 79 : f = mkvec2(gen_0, Fp_sqr(c6,p));
957 : }
958 : else
959 : {
960 0 : KRO = -KRO;
961 0 : f = Fp_ellpoint(KRO, &x, c4,c6,p);
962 : }
963 : /* [ux, u^2] is on E_u: y^2 = x^3 + c4 u^2 x + c6 u^3
964 : * E_u isomorphic to E (resp. E') iff KRO = 1 (resp. -1)
965 : * #E(F_p) = p+1 - a_p, #E'(F_p) = p+1 + a_p
966 : *
967 : * #E_u(Fp) = A (mod B), h is close to #E_u(Fp) */
968 79 : a4 = modii(mulii(c4, gel(f,2)), p); /* c4 for E_u */
969 79 : fh = FpE_mul(f, h, a4, p);
970 79 : if (ell_is_inf(fh)) goto FOUND;
971 :
972 79 : s = get_table_size(pordmin, B);
973 : /* look for h s.t f^h = 0 */
974 79 : if (!tx)
975 : { /* first time: initialize */
976 79 : tx = newblock(3*(s+1));
977 79 : ty = tx + (s+1);
978 79 : ti = ty + (s+1);
979 : }
980 79 : F = FpE_mul(f,B,a4,p);
981 79 : *tx = evaltyp(t_VECSMALL) | evallg(s+1);
982 :
983 : /* F = B.f */
984 79 : P = gcopy(fh);
985 79 : if (s < 3)
986 : { /* we're nearly done: naive search */
987 0 : GEN q1 = P, mF = FpE_neg(F, p); /* -F */
988 0 : for (i=1;; i++)
989 : {
990 0 : P = FpE_add(P,F,a4,p); /* h.f + i.F */
991 0 : if (ell_is_inf(P)) { h = addii(h, mului(i,B)); goto FOUND; }
992 0 : q1 = FpE_add(q1,mF,a4,p); /* h.f - i.F */
993 0 : if (ell_is_inf(q1)) { h = subii(h, mului(i,B)); goto FOUND; }
994 : }
995 : }
996 : /* Baby Step/Giant Step */
997 79 : nb = minss(128, s >> 1); /* > 0. Will do nb pts at a time: faster inverse */
998 79 : pts = cgetg(nb+1, t_VEC);
999 79 : j = lgefint(p);
1000 9677 : for (i=1; i<=nb; i++)
1001 : { /* baby steps */
1002 9598 : gel(pts,i) = P; /* h.f + (i-1).F */
1003 9598 : _fix(P+1, j); tx[i] = mod2BIL(gel(P,1));
1004 9598 : _fix(P+2, j); ty[i] = mod2BIL(gel(P,2));
1005 9598 : P = FpE_add(P,F,a4,p); /* h.f + i.F */
1006 9598 : if (ell_is_inf(P)) { h = addii(h, mului(i,B)); goto FOUND; }
1007 : }
1008 79 : mfh = FpE_neg(fh, p);
1009 79 : fg = FpE_add(P,mfh,a4,p); /* h.f + nb.F - h.f = nb.F */
1010 79 : if (ell_is_inf(fg)) { h = mului(nb,B); goto FOUND; }
1011 79 : u = cgetg(nb+1, t_VEC);
1012 79 : av2 = avma; /* more baby steps, nb points at a time */
1013 1357 : while (i <= s)
1014 : {
1015 : long maxj;
1016 164239 : for (j=1; j<=nb; j++) /* adding nb.F (part 1) */
1017 : {
1018 162961 : P = gel(pts,j); /* h.f + (i-nb-1+j-1).F */
1019 162961 : gel(u,j) = subii(gel(fg,1), gel(P,1));
1020 162961 : if (!signe(gel(u,j))) /* sum = 0 or doubling */
1021 : {
1022 2 : long k = i+j-2;
1023 2 : if (equalii(gel(P,2),gel(fg,2))) k -= 2*nb; /* fg == P */
1024 2 : h = addii(h, mulsi(k,B)); goto FOUND;
1025 : }
1026 : }
1027 1278 : v = FpV_inv(u, p);
1028 1278 : maxj = (i-1 + nb <= s)? nb: s % nb;
1029 160545 : for (j=1; j<=maxj; j++,i++) /* adding nb.F (part 2) */
1030 : {
1031 159267 : P = gel(pts,j);
1032 159267 : FpE_add_ip(P,fg, a4,p, gel(v,j));
1033 159267 : tx[i] = mod2BIL(gel(P,1));
1034 159267 : ty[i] = mod2BIL(gel(P,2));
1035 : }
1036 1278 : set_avma(av2);
1037 : }
1038 77 : P = FpE_add(gel(pts,j-1),mfh,a4,p); /* = (s-1).F */
1039 77 : if (ell_is_inf(P)) { h = mului(s-1,B); goto FOUND; }
1040 77 : if (DEBUGLEVEL >= 6)
1041 0 : timer_printf(&T, "[Fp_ellcard_Shanks] baby steps, s = %ld",s);
1042 :
1043 : /* giant steps: fg = s.F */
1044 77 : fg = FpE_add(P,F,a4,p);
1045 77 : if (ell_is_inf(fg)) { h = mului(s,B); goto FOUND; }
1046 77 : pfinal = mod2BIL(p); av2 = avma;
1047 : /* Goal of the following: sort points by increasing x-coordinate hash.
1048 : * Done in a complicated way to avoid allocating a large temp vector */
1049 77 : p1 = vecsmall_indexsort(tx); /* = permutation sorting tx */
1050 168784 : for (i=1; i<=s; i++) ti[i] = tx[p1[i]];
1051 : /* ti = tx sorted */
1052 168784 : for (i=1; i<=s; i++) { tx[i] = ti[i]; ti[i] = ty[p1[i]]; }
1053 : /* tx is sorted. ti = ty sorted */
1054 168784 : for (i=1; i<=s; i++) { ty[i] = ti[i]; ti[i] = p1[i]; }
1055 : /* ty is sorted. ti = permutation sorting tx */
1056 77 : if (DEBUGLEVEL >= 6) timer_printf(&T, "[Fp_ellcard_Shanks] sorting");
1057 77 : set_avma(av2);
1058 :
1059 77 : affii2(fg, gel(pts,1));
1060 9440 : for (j=2; j<=nb; j++) /* pts[j] = j.fg = (s*j).F */
1061 : {
1062 9363 : P = FpE_add(gel(pts,j-1),fg,a4,p);
1063 9363 : if (ell_is_inf(P)) { h = mulii(mulss(s,j), B); goto FOUND; }
1064 9363 : affii2(P, gel(pts,j));
1065 : }
1066 : /* replace fg by nb.fg since we do nb points at a time */
1067 77 : set_avma(av2);
1068 77 : fg = gcopy(gel(pts,nb)); /* copy: we modify (temporarily) pts[nb] below */
1069 77 : av2 = avma;
1070 :
1071 77 : for (i=1,j=1; ; i++)
1072 152075 : {
1073 152152 : GEN ftest = gel(pts,j);
1074 152152 : long m, l = 1, r = s+1;
1075 : long k, k2, j2;
1076 :
1077 152152 : set_avma(av2);
1078 152152 : k = mod2BIL(gel(ftest,1));
1079 1930966 : while (l < r)
1080 : {
1081 1778814 : m = (l+r) >> 1;
1082 1778814 : if (tx[m] < k) l = m+1; else r = m;
1083 : }
1084 152152 : if (r <= s && tx[r] == k)
1085 : {
1086 154 : while (r && tx[r] == k) r--;
1087 77 : k2 = mod2BIL(gel(ftest,2));
1088 77 : for (r++; r <= s && tx[r] == k; r++)
1089 77 : if (ty[r] == k2 || ty[r] == pfinal - k2)
1090 : { /* [h+j2] f == +/- ftest (= [i.s] f)? */
1091 77 : j2 = ti[r] - 1;
1092 77 : if (DEBUGLEVEL >=6)
1093 0 : timer_printf(&T, "[Fp_ellcard_Shanks] giant steps, i = %ld",i);
1094 77 : P = FpE_add(FpE_mul(F,stoi(j2),a4,p),fh,a4,p);
1095 77 : if (equalii(gel(P,1), gel(ftest,1)))
1096 : {
1097 77 : if (equalii(gel(P,2), gel(ftest,2))) i = -i;
1098 77 : h = addii(h, mulii(addis(mulss(s,i), j2), B));
1099 77 : goto FOUND;
1100 : }
1101 : }
1102 : }
1103 152075 : if (++j > nb)
1104 : { /* compute next nb points */
1105 1149 : long save = 0; /* gcc -Wall */;
1106 147576 : for (j=1; j<=nb; j++)
1107 : {
1108 146427 : P = gel(pts,j);
1109 146427 : gel(u,j) = subii(gel(fg,1), gel(P,1));
1110 146427 : if (gel(u,j) == gen_0) /* occurs once: i = j = nb, P == fg */
1111 : {
1112 67 : gel(u,j) = shifti(gel(P,2),1);
1113 67 : save = fg[1]; fg[1] = P[1];
1114 : }
1115 : }
1116 1149 : v = FpV_inv(u, p);
1117 147576 : for (j=1; j<=nb; j++)
1118 146427 : FpE_add_ip(gel(pts,j),fg,a4,p, gel(v,j));
1119 1149 : if (i == nb) { fg[1] = save; }
1120 1149 : j = 1;
1121 : }
1122 : }
1123 79 : FOUND: /* found a point of exponent h on E_u */
1124 79 : h = FpE_order(f, h, a4, p);
1125 : /* h | #E_u(Fp) = A (mod B) */
1126 79 : A = Z_chinese_all(A, gen_0, B, h, &B);
1127 79 : if (cmpii(B, pordmin) >= 0) break;
1128 : /* not done: update A mod B for the _next_ curve, isomorphic to
1129 : * the quadratic twist of this one */
1130 0 : A = remii(subii(p2p,A), B); /* #E(Fp)+#E'(Fp) = 2p+2 */
1131 : }
1132 79 : if (tx) killblock(tx);
1133 79 : h = closest_lift(A, B, p1p);
1134 79 : return gc_INT(av, KRO==1? h: subii(p2p,h));
1135 : }
1136 :
1137 : typedef struct
1138 : {
1139 : ulong x,y,i;
1140 : } multiple;
1141 :
1142 : static int
1143 15379662 : compare_multiples(const void *A, const void *B)
1144 : {
1145 15379662 : multiple *a = (multiple*)A, *b = (multiple*)B;
1146 15379662 : return a->x > b->x? 1: a->x < b->x? -1: 0;
1147 : }
1148 :
1149 : /* find x such that h := a + b x is closest to c and return h:
1150 : * x = round((c-a) / b) = floor( (2(c-a) + b) / 2b )
1151 : * Assume 0 <= a < b < c and b + 2c < 2^BIL */
1152 : static ulong
1153 261973 : uclosest_lift(ulong a, ulong b, ulong c)
1154 : {
1155 261973 : ulong x = (b + ((c-a) << 1)) / (b << 1);
1156 261973 : return a + b * x;
1157 : }
1158 :
1159 : static long
1160 227176 : Fle_dbl_inplace(GEN P, ulong a4, ulong p)
1161 : {
1162 : ulong x, y, slope;
1163 227176 : if (!P[2]) return 1;
1164 227148 : x = P[1]; y = P[2];
1165 227148 : slope = Fl_div(Fl_add(Fl_triple(Fl_sqr(x,p), p), a4, p),
1166 : Fl_double(y, p), p);
1167 227152 : P[1] = Fl_sub(Fl_sqr(slope, p), Fl_double(x, p), p);
1168 227147 : P[2] = Fl_sub(Fl_mul(slope, Fl_sub(x, P[1], p), p), y, p);
1169 227145 : return 0;
1170 : }
1171 :
1172 : static long
1173 5797654 : Fle_add_inplace(GEN P, GEN Q, ulong a4, ulong p)
1174 : {
1175 : ulong Px, Py, Qx, Qy, slope;
1176 5797654 : if (ell_is_inf(Q)) return 0;
1177 5797603 : Px = P[1]; Py = P[2];
1178 5797603 : Qx = Q[1]; Qy = Q[2];
1179 5797603 : if (Px==Qx)
1180 238637 : return Py==Qy ? Fle_dbl_inplace(P, a4, p): 1;
1181 5558966 : slope = Fl_div(Fl_sub(Py, Qy, p), Fl_sub(Px, Qx, p), p);
1182 5559703 : P[1] = Fl_sub(Fl_sub(Fl_sqr(slope, p), Px, p), Qx, p);
1183 5559132 : P[2] = Fl_sub(Fl_mul(slope, Fl_sub(Px, P[1], p), p), Py, p);
1184 5558862 : return 0;
1185 : }
1186 :
1187 : /* assume 99 < p < 2^(BIL-1) - 2^((BIL+1)/2) and e has good reduction at p.
1188 : * Should use Barett reduction + multi-inverse. See Fp_ellcard_Shanks() */
1189 : static long
1190 254758 : Fl_ellcard_Shanks(ulong c4, ulong c6, ulong p)
1191 : {
1192 : GEN f, fh, fg, ftest, F;
1193 : ulong i, l, r, s, h, x, cp4, p1p, p2p, pordmin,A,B;
1194 : long KRO;
1195 254758 : pari_sp av = avma;
1196 : multiple *table;
1197 :
1198 254758 : if (!c6) {
1199 14 : GEN ap = ap_j1728(utoi(c4), utoipos(p));
1200 14 : return gc_long(av, p+1 - itos(ap));
1201 : }
1202 :
1203 254744 : pordmin = (ulong)(1 + 4*sqrt((double)p));
1204 254744 : p1p = p+1;
1205 254744 : p2p = p1p << 1;
1206 254744 : x = 0; KRO = 0;
1207 254744 : switch(Flx_nbroots(mkvecsmall5(0L, c6,c4,0L,1L), p))
1208 : {
1209 51716 : case 3: A = 0; B = 4; break;
1210 124413 : case 1: A = 0; B = 2; break;
1211 78619 : default: A = 1; B = 2; break; /* 0 */
1212 : }
1213 : for(;;)
1214 : { /* see comments in Fp_ellcard_Shanks */
1215 261975 : h = uclosest_lift(A, B, p1p);
1216 261973 : if (!KRO) /* first time, initialize */
1217 : {
1218 254746 : KRO = krouu(c6,p); /* != 0 */
1219 254751 : f = mkvecsmall2(0, Fl_sqr(c6,p));
1220 : }
1221 : else
1222 : {
1223 7227 : KRO = -KRO;
1224 7227 : f = Fl_ellpoint(KRO, &x, c4,c6,p);
1225 : }
1226 261977 : cp4 = Fl_mul(c4, f[2], p);
1227 261976 : fh = Fle_mulu(f, h, cp4, p);
1228 261970 : if (ell_is_inf(fh)) goto FOUND;
1229 :
1230 255763 : s = (ulong) (sqrt(((double)pordmin)/B) / 2);
1231 255763 : if (!s) s = 1;
1232 255763 : table = (multiple *) stack_malloc((s+1) * sizeof(multiple));
1233 255762 : F = Fle_mulu(f, B, cp4, p);
1234 3347911 : for (i=0; i < s; i++)
1235 : {
1236 3103633 : table[i].x = fh[1];
1237 3103633 : table[i].y = fh[2];
1238 3103633 : table[i].i = i;
1239 3103633 : if (Fle_add_inplace(fh, F, cp4, p)) { h += B*(i+1); goto FOUND; }
1240 : }
1241 244278 : qsort(table,s,sizeof(multiple),compare_multiples);
1242 244283 : fg = Fle_mulu(F, s, cp4, p); ftest = zv_copy(fg);
1243 244270 : if (ell_is_inf(ftest)) {
1244 0 : if (!uisprime(p)) pari_err_PRIME("ellap",utoi(p));
1245 0 : pari_err_BUG("ellap (f^(i*s) = 1)");
1246 : }
1247 2938961 : for (i=1; ; i++)
1248 : {
1249 2938961 : l=0; r=s;
1250 20641797 : while (l<r)
1251 : {
1252 17702836 : ulong m = (l+r) >> 1;
1253 17702836 : if (table[m].x < uel(ftest,1)) l=m+1; else r=m;
1254 : }
1255 2938961 : if (r < s && table[r].x == uel(ftest,1)) break;
1256 2694682 : if (Fle_add_inplace(ftest, fg, cp4, p)) pari_err_PRIME("ellap",utoi(p));
1257 : }
1258 244279 : h += table[r].i * B;
1259 244279 : if (table[r].y == uel(ftest,2))
1260 126878 : h -= s * i * B;
1261 : else
1262 117401 : h += s * i * B;
1263 261975 : FOUND:
1264 261975 : h = itou(Fle_order(f, utoipos(h), cp4, p));
1265 : /* h | #E_u(Fp) = A (mod B) */
1266 : {
1267 : GEN C;
1268 261967 : A = itou( Z_chinese_all(gen_0, utoi(A), utoipos(h), utoipos(B), &C) );
1269 261962 : if (abscmpiu(C, pordmin) >= 0) { /* uclosest_lift could overflow */
1270 254735 : h = itou( closest_lift(utoi(A), C, utoipos(p1p)) );
1271 254750 : break;
1272 : }
1273 7227 : B = itou(C);
1274 : }
1275 7227 : A = (p2p - A) % B; set_avma(av);
1276 : }
1277 254750 : return gc_long(av, KRO==1? h: p2p-h);
1278 : }
1279 :
1280 : /** ellap from CM (original code contributed by Mark Watkins) **/
1281 :
1282 : static GEN
1283 85164 : ap_j0(GEN a6,GEN p)
1284 : {
1285 : GEN a, b, e, d;
1286 85164 : if (umodiu(p,3) != 1) return gen_0;
1287 42270 : (void)cornacchia2(utoipos(27),p, &a,&b);
1288 42443 : if (umodiu(a, 3) == 1) a = negi(a);
1289 42444 : d = mulis(a6,-108);
1290 42413 : e = diviuexact(shifti(p,-1), 3); /* (p-1) / 6 */
1291 42388 : return centermod(mulii(a, Fp_pow(d, e, p)), p);
1292 : }
1293 : static GEN
1294 2642444 : ap_j1728(GEN a4,GEN p)
1295 : {
1296 : GEN a, b, e;
1297 2642444 : if (mod4(p) != 1) return gen_0;
1298 1320221 : (void)cornacchia2(utoipos(4),p, &a,&b);
1299 1320221 : if (Mod4(a)==0) a = b;
1300 1320221 : if (Mod2(a)==1) a = shifti(a,1);
1301 1320221 : if (Mod8(a)==6) a = negi(a);
1302 1320221 : e = shifti(p,-2); /* (p-1) / 4 */
1303 1320221 : return centermod(mulii(a, Fp_pow(a4, e, p)), p);
1304 : }
1305 : static GEN
1306 126 : ap_j8000(GEN a6, GEN p)
1307 : {
1308 : GEN a, b;
1309 126 : long r = mod8(p), s = 1;
1310 126 : if (r != 1 && r != 3) return gen_0;
1311 49 : (void)cornacchia2(utoipos(8),p, &a,&b);
1312 49 : switch(Mod16(a)) {
1313 14 : case 2: case 6: if (Mod4(b)) s = -s;
1314 14 : break;
1315 35 : case 10: case 14: if (!Mod4(b)) s = -s;
1316 35 : break;
1317 : }
1318 49 : if (kronecker(mulis(a6, 42), p) < 0) s = -s;
1319 49 : return s > 0? a: negi(a);
1320 : }
1321 : static GEN
1322 140 : ap_j287496(GEN a6, GEN p)
1323 : {
1324 : GEN a, b;
1325 140 : long s = 1;
1326 140 : if (mod4(p) != 1) return gen_0;
1327 70 : (void)cornacchia2(utoipos(4),p, &a,&b);
1328 70 : if (Mod4(a)==0) a = b;
1329 70 : if (Mod2(a)==1) a = shifti(a,1);
1330 70 : if (Mod8(a)==6) s = -s;
1331 70 : if (krosi(2,p) < 0) s = -s;
1332 70 : if (kronecker(mulis(a6, -14), p) < 0) s = -s;
1333 70 : return s > 0? a: negi(a);
1334 : }
1335 : static GEN
1336 1344 : ap_cm(int CM, long A6B, GEN a6, GEN p)
1337 : {
1338 : GEN a, b;
1339 1344 : long s = 1;
1340 1344 : if (krosi(CM,p) < 0) return gen_0;
1341 644 : (void)cornacchia2(utoipos(-CM),p, &a, &b);
1342 644 : if ((CM&3) == 0) CM >>= 2;
1343 644 : if ((krois(a, -CM) > 0) ^ (CM == -7)) s = -s;
1344 644 : if (kronecker(mulis(a6,A6B), p) < 0) s = -s;
1345 644 : return s > 0? a: negi(a);
1346 : }
1347 : static GEN
1348 497483 : ec_ap_cm(int CM, GEN a4, GEN a6, GEN p)
1349 : {
1350 497483 : switch(CM)
1351 : {
1352 29113 : case -3: return ap_j0(a6, p);
1353 466760 : case -4: return ap_j1728(a4, p);
1354 126 : case -8: return ap_j8000(a6, p);
1355 140 : case -16: return ap_j287496(a6, p);
1356 154 : case -7: return ap_cm(CM, -2, a6, p);
1357 147 : case -11: return ap_cm(CM, 21, a6, p);
1358 168 : case -12: return ap_cm(CM, 22, a6, p);
1359 147 : case -19: return ap_cm(CM, 1, a6, p);
1360 154 : case -27: return ap_cm(CM, 253, a6, p);
1361 140 : case -28: return ap_cm(-7, -114, a6, p); /* yes, -7 ! */
1362 147 : case -43: return ap_cm(CM, 21, a6, p);
1363 147 : case -67: return ap_cm(CM, 217, a6, p);
1364 140 : case -163:return ap_cm(CM, 185801, a6, p);
1365 0 : default: return NULL;
1366 : }
1367 : }
1368 :
1369 : static GEN
1370 49117 : Fp_ellj_nodiv(GEN a4, GEN a6, GEN p)
1371 : {
1372 49117 : GEN a43 = Fp_mulu(Fp_powu(a4, 3, p), 4, p);
1373 49116 : GEN a62 = Fp_mulu(Fp_sqr(a6, p), 27, p);
1374 49116 : return mkvec2(Fp_mulu(a43, 1728, p), Fp_add(a43, a62, p));
1375 : }
1376 :
1377 : GEN
1378 56 : Fp_ellj(GEN a4, GEN a6, GEN p)
1379 : {
1380 56 : pari_sp av = avma;
1381 : GEN z;
1382 56 : if (lgefint(p) == 3)
1383 : {
1384 0 : ulong pp = p[2];
1385 0 : return utoi(Fl_ellj(umodiu(a4,pp), umodiu(a6,pp), pp));
1386 : }
1387 56 : z = Fp_ellj_nodiv(a4, a6, p);
1388 56 : return gc_INT(av,Fp_div(gel(z,1),gel(z,2),p));
1389 : }
1390 :
1391 : void
1392 1050 : Fp_ellj_to_a4a6(GEN j, GEN p, GEN *pt_a4, GEN *pt_a6)
1393 : {
1394 1050 : j = modii(j, p);
1395 1050 : if (signe(j) == 0) { *pt_a4 = gen_0; *pt_a6 = gen_1; }
1396 686 : else if (equaliu(j,umodui(1728,p))) { *pt_a4 = gen_1; *pt_a6 = gen_0; }
1397 : else
1398 : {
1399 546 : GEN k = Fp_sub(utoi(1728), j, p);
1400 546 : GEN kj = Fp_mul(k, j, p);
1401 546 : GEN k2j = Fp_mul(kj, k, p);
1402 546 : *pt_a4 = Fp_mulu(kj, 3, p);
1403 546 : *pt_a6 = Fp_double(k2j, p);
1404 : }
1405 1050 : }
1406 :
1407 : static GEN /* Only compute a mod p, so assume p>=17 */
1408 2280764 : Fp_ellcard_CM(GEN a4, GEN a6, GEN p)
1409 : {
1410 2280764 : pari_sp av = avma;
1411 : GEN a;
1412 2280764 : if (!signe(a4)) a = ap_j0(a6,p);
1413 2224722 : else if (!signe(a6)) a = ap_j1728(a4,p);
1414 : else
1415 : {
1416 49052 : GEN j = Fp_ellj_nodiv(a4, a6, p);
1417 49053 : long CM = Fp_ellj_get_CM(gel(j,1), gel(j,2), p);
1418 49058 : if (!CM) return gc_NULL(av);
1419 1610 : a = ec_ap_cm(CM,a4,a6,p);
1420 : }
1421 2233480 : return gc_INT(av, subii(addiu(p,1),a));
1422 : }
1423 :
1424 : GEN
1425 2542635 : Fp_ellcard(GEN a4, GEN a6, GEN p)
1426 : {
1427 2542635 : long lp = expi(p);
1428 2542627 : ulong pp = p[2];
1429 2542627 : if (lp < 11)
1430 261900 : return utoi(pp+1 - Fl_elltrace_naive(umodiu(a4,pp), umodiu(a6,pp), pp));
1431 2280727 : { GEN a = Fp_ellcard_CM(a4,a6,p); if (a) return a; }
1432 47447 : if (lp >= 56)
1433 868 : return Fp_ellcard_SEA(a4, a6, p, 0);
1434 46579 : if (lp <= BITS_IN_LONG-2)
1435 46501 : return utoi(Fl_ellcard_Shanks(umodiu(a4,pp), umodiu(a6,pp), pp));
1436 79 : return Fp_ellcard_Shanks(a4, a6, p);
1437 : }
1438 :
1439 : long
1440 621556 : Fl_elltrace(ulong a4, ulong a6, ulong p)
1441 : {
1442 : pari_sp av;
1443 : long lp;
1444 : GEN a;
1445 621556 : if (p < (1<<11)) return Fl_elltrace_naive(a4, a6, p);
1446 208254 : lp = expu(p);
1447 208254 : if (lp <= minss(56, BITS_IN_LONG-2)) return p+1-Fl_ellcard_Shanks(a4, a6, p);
1448 0 : av = avma; a = subui(p+1, Fp_ellcard(utoi(a4), utoi(a6), utoipos(p)));
1449 0 : return gc_long(av, itos(a));
1450 : }
1451 : long
1452 1164937 : Fl_elltrace_CM(long CM, ulong a4, ulong a6, ulong p)
1453 : {
1454 : pari_sp av;
1455 : GEN a;
1456 1164937 : if (!CM) return Fl_elltrace(a4,a6,p);
1457 544110 : if (p < (1<<11)) return Fl_elltrace_naive(a4, a6, p);
1458 495873 : av = avma; a = ec_ap_cm(CM, utoi(a4), utoi(a6), utoipos(p));
1459 495873 : return gc_long(av, itos(a));
1460 : }
1461 :
1462 : static GEN
1463 72725 : _FpE_pairorder(void *E, GEN P, GEN Q, GEN m, GEN F)
1464 : {
1465 72725 : struct _FpE *e = (struct _FpE *) E;
1466 72725 : return Fp_order(FpE_weilpairing(P,Q,m,e->a4,e->p), F, e->p);
1467 : }
1468 :
1469 : GEN
1470 120715 : Fp_ellgroup(GEN a4, GEN a6, GEN N, GEN p, GEN *pt_m)
1471 : {
1472 : struct _FpE e;
1473 120715 : e.a4=a4; e.a6=a6; e.p=p;
1474 120715 : return gen_ellgroup(N, subiu(p,1), pt_m, (void*)&e, &FpE_group, _FpE_pairorder);
1475 : }
1476 :
1477 : GEN
1478 574 : Fp_ellgens(GEN a4, GEN a6, GEN ch, GEN D, GEN m, GEN p)
1479 : {
1480 : GEN P;
1481 574 : pari_sp av = avma;
1482 : struct _FpE e;
1483 574 : e.a4=a4; e.a6=a6; e.p=p;
1484 574 : switch(lg(D)-1)
1485 : {
1486 476 : case 1:
1487 476 : P = gen_gener(gel(D,1), (void*)&e, &FpE_group);
1488 476 : P = mkvec(FpE_changepoint(P, ch, p));
1489 476 : break;
1490 98 : default:
1491 98 : P = gen_ellgens(gel(D,1), gel(D,2), m, (void*)&e, &FpE_group, _FpE_pairorder);
1492 98 : gel(P,1) = FpE_changepoint(gel(P,1), ch, p);
1493 98 : gel(P,2) = FpE_changepoint(gel(P,2), ch, p);
1494 98 : break;
1495 : }
1496 574 : return gc_GEN(av, P);
1497 : }
1498 :
1499 : /* Not so fast arithmetic with points over elliptic curves over FpXQ */
1500 :
1501 : /***********************************************************************/
1502 : /** **/
1503 : /** FpXQE **/
1504 : /** **/
1505 : /***********************************************************************/
1506 : /* These functions deal with point over elliptic curves over FpXQ defined
1507 : * by an equation of the form y^2=x^3+a4*x+a6.
1508 : * Most of the time a6 is omitted since it can be recovered from any point
1509 : * on the curve. */
1510 :
1511 : GEN
1512 976 : RgE_to_FpXQE(GEN x, GEN T, GEN p)
1513 : {
1514 976 : if (ell_is_inf(x)) return x;
1515 976 : retmkvec2(Rg_to_FpXQ(gel(x,1),T,p),Rg_to_FpXQ(gel(x,2),T,p));
1516 : }
1517 :
1518 : GEN
1519 1876 : FpXQE_changepoint(GEN x, GEN ch, GEN T, GEN p)
1520 : {
1521 1876 : pari_sp av = avma;
1522 : GEN p1,z,u,r,s,t,v,v2,v3;
1523 1876 : if (ell_is_inf(x)) return x;
1524 942 : u = gel(ch,1); r = gel(ch,2);
1525 942 : s = gel(ch,3); t = gel(ch,4);
1526 942 : v = FpXQ_inv(u, T, p); v2 = FpXQ_sqr(v, T, p); v3 = FpXQ_mul(v,v2, T, p);
1527 942 : p1 = FpX_sub(gel(x,1),r, p);
1528 942 : z = cgetg(3,t_VEC);
1529 942 : gel(z,1) = FpXQ_mul(v2, p1, T, p);
1530 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);
1531 942 : return gc_upto(av, z);
1532 : }
1533 :
1534 : GEN
1535 976 : FpXQE_changepointinv(GEN x, GEN ch, GEN T, GEN p)
1536 : {
1537 : GEN u, r, s, t, X, Y, u2, u3, u2X, z;
1538 976 : if (ell_is_inf(x)) return x;
1539 976 : X = gel(x,1); Y = gel(x,2);
1540 976 : u = gel(ch,1); r = gel(ch,2);
1541 976 : s = gel(ch,3); t = gel(ch,4);
1542 976 : u2 = FpXQ_sqr(u, T, p); u3 = FpXQ_mul(u,u2, T, p);
1543 976 : u2X = FpXQ_mul(u2,X, T, p);
1544 976 : z = cgetg(3, t_VEC);
1545 976 : gel(z,1) = FpX_add(u2X,r, p);
1546 976 : gel(z,2) = FpX_add(FpXQ_mul(u3,Y, T, p), FpX_add(FpXQ_mul(s,u2X, T, p), t, p), p);
1547 976 : return z;
1548 : }
1549 :
1550 : static GEN
1551 840 : random_nonsquare_FpXQ(GEN T, GEN p)
1552 : {
1553 840 : pari_sp av = avma;
1554 840 : long n = degpol(T), v = varn(T);
1555 : GEN a;
1556 840 : if (odd(n))
1557 : {
1558 420 : GEN z = cgetg(3, t_POL);
1559 420 : z[1] = evalsigne(1) | evalvarn(v);
1560 420 : gel(z,2) = random_nonsquare_Fp(p); return z;
1561 : }
1562 : do
1563 : {
1564 791 : set_avma(av);
1565 791 : a = random_FpX(n, v, p);
1566 791 : } while (FpXQ_issquare(a, T, p));
1567 420 : return a;
1568 : }
1569 :
1570 : void
1571 840 : FpXQ_elltwist(GEN a4, GEN a6, GEN T, GEN p, GEN *pt_a4, GEN *pt_a6)
1572 : {
1573 840 : GEN d = random_nonsquare_FpXQ(T, p);
1574 840 : GEN d2 = FpXQ_sqr(d, T, p), d3 = FpXQ_mul(d2, d, T, p);
1575 840 : *pt_a4 = FpXQ_mul(a4, d2, T, p);
1576 840 : *pt_a6 = FpXQ_mul(a6, d3, T, p);
1577 840 : }
1578 :
1579 : static GEN
1580 269139 : FpXQE_dbl_slope(GEN P, GEN a4, GEN T, GEN p, GEN *slope)
1581 : {
1582 : GEN x, y, Q;
1583 269139 : if (ell_is_inf(P) || !signe(gel(P,2))) return ellinf();
1584 267496 : x = gel(P,1); y = gel(P,2);
1585 267496 : *slope = FpXQ_div(FpX_add(FpX_mulu(FpXQ_sqr(x, T, p), 3, p), a4, p),
1586 : FpX_mulu(y, 2, p), T, p);
1587 267496 : Q = cgetg(3,t_VEC);
1588 267496 : gel(Q, 1) = FpX_sub(FpXQ_sqr(*slope, T, p), FpX_mulu(x, 2, p), p);
1589 267496 : gel(Q, 2) = FpX_sub(FpXQ_mul(*slope, FpX_sub(x, gel(Q, 1), p), T, p), y, p);
1590 267496 : return Q;
1591 : }
1592 :
1593 : GEN
1594 265037 : FpXQE_dbl(GEN P, GEN a4, GEN T, GEN p)
1595 : {
1596 265037 : pari_sp av = avma;
1597 : GEN slope;
1598 265037 : return gc_upto(av, FpXQE_dbl_slope(P,a4,T,p,&slope));
1599 : }
1600 :
1601 : static GEN
1602 252060 : FpXQE_add_slope(GEN P, GEN Q, GEN a4, GEN T, GEN p, GEN *slope)
1603 : {
1604 : GEN Px, Py, Qx, Qy, R;
1605 252060 : if (ell_is_inf(P)) return Q;
1606 252046 : if (ell_is_inf(Q)) return P;
1607 252046 : Px = gel(P,1); Py = gel(P,2);
1608 252046 : Qx = gel(Q,1); Qy = gel(Q,2);
1609 252046 : if (ZX_equal(Px, Qx))
1610 : {
1611 687 : if (ZX_equal(Py, Qy))
1612 7 : return FpXQE_dbl_slope(P, a4, T, p, slope);
1613 : else
1614 680 : return ellinf();
1615 : }
1616 251359 : *slope = FpXQ_div(FpX_sub(Py, Qy, p), FpX_sub(Px, Qx, p), T, p);
1617 251359 : R = cgetg(3,t_VEC);
1618 251359 : gel(R, 1) = FpX_sub(FpX_sub(FpXQ_sqr(*slope, T, p), Px, p), Qx, p);
1619 251359 : gel(R, 2) = FpX_sub(FpXQ_mul(*slope, FpX_sub(Px, gel(R, 1), p), T, p), Py, p);
1620 251359 : return R;
1621 : }
1622 :
1623 : GEN
1624 251500 : FpXQE_add(GEN P, GEN Q, GEN a4, GEN T, GEN p)
1625 : {
1626 251500 : pari_sp av = avma;
1627 : GEN slope;
1628 251500 : return gc_upto(av, FpXQE_add_slope(P,Q,a4,T,p,&slope));
1629 : }
1630 :
1631 : static GEN
1632 0 : FpXQE_neg_i(GEN P, GEN p)
1633 : {
1634 0 : if (ell_is_inf(P)) return P;
1635 0 : return mkvec2(gel(P,1), FpX_neg(gel(P,2), p));
1636 : }
1637 :
1638 : GEN
1639 73329 : FpXQE_neg(GEN P, GEN T, GEN p)
1640 : {
1641 : (void) T;
1642 73329 : if (ell_is_inf(P)) return ellinf();
1643 73329 : return mkvec2(gcopy(gel(P,1)), FpX_neg(gel(P,2), p));
1644 : }
1645 :
1646 : GEN
1647 0 : FpXQE_sub(GEN P, GEN Q, GEN a4, GEN T, GEN p)
1648 : {
1649 0 : pari_sp av = avma;
1650 : GEN slope;
1651 0 : return gc_upto(av, FpXQE_add_slope(P, FpXQE_neg_i(Q, p), a4, T, p, &slope));
1652 : }
1653 :
1654 : struct _FpXQE { GEN a4,a6,T,p; };
1655 : static GEN
1656 265037 : _FpXQE_dbl(void *E, GEN P)
1657 : {
1658 265037 : struct _FpXQE *ell = (struct _FpXQE *) E;
1659 265037 : return FpXQE_dbl(P, ell->a4, ell->T, ell->p);
1660 : }
1661 : static GEN
1662 251500 : _FpXQE_add(void *E, GEN P, GEN Q)
1663 : {
1664 251500 : struct _FpXQE *ell=(struct _FpXQE *) E;
1665 251500 : return FpXQE_add(P, Q, ell->a4, ell->T, ell->p);
1666 : }
1667 : static GEN
1668 81874 : _FpXQE_mul(void *E, GEN P, GEN n)
1669 : {
1670 81874 : pari_sp av = avma;
1671 81874 : struct _FpXQE *e=(struct _FpXQE *) E;
1672 81874 : long s = signe(n);
1673 81874 : if (!s || ell_is_inf(P)) return ellinf();
1674 81874 : if (s<0) P = FpXQE_neg(P, e->T, e->p);
1675 81874 : if (is_pm1(n)) return s>0? gcopy(P): P;
1676 8453 : return gc_GEN(av, gen_pow_i(P, n, e, &_FpXQE_dbl, &_FpXQE_add));
1677 : }
1678 :
1679 : GEN
1680 934 : FpXQE_mul(GEN P, GEN n, GEN a4, GEN T, GEN p)
1681 : {
1682 : struct _FpXQE E;
1683 934 : E.a4= a4; E.T = T; E.p = p;
1684 934 : return _FpXQE_mul(&E, P, n);
1685 : }
1686 :
1687 : /* Finds a random nonsingular point on E */
1688 :
1689 : GEN
1690 1081 : random_FpXQE(GEN a4, GEN a6, GEN T, GEN p)
1691 : {
1692 1081 : pari_sp ltop = avma;
1693 : GEN x, x2, y, rhs;
1694 1081 : long v = get_FpX_var(T), d = get_FpX_degree(T);
1695 : do
1696 : {
1697 2208 : set_avma(ltop);
1698 2208 : x = random_FpX(d,v,p); /* x^3+a4*x+a6 = x*(x^2+a4)+a6 */
1699 2208 : x2 = FpXQ_sqr(x, T, p);
1700 2208 : rhs = FpX_add(FpXQ_mul(x, FpX_add(x2, a4, p), T, p), a6, p);
1701 0 : } while ((!signe(rhs) && !signe(FpX_add(FpX_mulu(x2,3,p), a4, p)))
1702 2208 : || !FpXQ_issquare(rhs, T, p));
1703 1081 : y = FpXQ_sqrt(rhs, T, p);
1704 1081 : if (!y) pari_err_PRIME("random_FpE", p);
1705 1081 : return gc_GEN(ltop, mkvec2(x, y));
1706 : }
1707 :
1708 : static GEN
1709 147 : _FpXQE_rand(void *E)
1710 : {
1711 147 : struct _FpXQE *e=(struct _FpXQE *) E;
1712 147 : return random_FpXQE(e->a4, e->a6, e->T, e->p);
1713 : }
1714 :
1715 : static const struct bb_group FpXQE_group={_FpXQE_add,_FpXQE_mul,_FpXQE_rand,hash_GEN,ZXV_equal,ell_is_inf};
1716 :
1717 : const struct bb_group *
1718 16 : get_FpXQE_group(void ** pt_E, GEN a4, GEN a6, GEN T, GEN p)
1719 : {
1720 16 : struct _FpXQE *e = (struct _FpXQE *) stack_malloc(sizeof(struct _FpXQE));
1721 16 : e->a4 = a4; e->a6 = a6; e->T = T; e->p = p;
1722 16 : *pt_E = (void *) e;
1723 16 : return &FpXQE_group;
1724 : }
1725 :
1726 : GEN
1727 14 : FpXQE_order(GEN z, GEN o, GEN a4, GEN T, GEN p)
1728 : {
1729 14 : pari_sp av = avma;
1730 : struct _FpXQE e;
1731 14 : e.a4=a4; e.T=T; e.p=p;
1732 14 : return gc_INT(av, gen_order(z, o, (void*)&e, &FpXQE_group));
1733 : }
1734 :
1735 : GEN
1736 0 : FpXQE_log(GEN a, GEN b, GEN o, GEN a4, GEN T, GEN p)
1737 : {
1738 0 : pari_sp av = avma;
1739 : struct _FpXQE e;
1740 0 : e.a4=a4; e.T=T; e.p=p;
1741 0 : return gc_INT(av, gen_PH_log(a, b, o, (void*)&e, &FpXQE_group));
1742 : }
1743 :
1744 : /***********************************************************************/
1745 : /** **/
1746 : /** Pairings **/
1747 : /** **/
1748 : /***********************************************************************/
1749 :
1750 : /* Derived from APIP from and by Jerome Milan, 2012 */
1751 :
1752 : static GEN
1753 4788 : FpXQE_vert(GEN P, GEN Q, GEN a4, GEN T, GEN p)
1754 : {
1755 4788 : long vT = get_FpX_var(T);
1756 4788 : if (ell_is_inf(P))
1757 70 : return pol_1(get_FpX_var(T));
1758 4718 : if (!ZX_equal(gel(Q, 1), gel(P, 1)))
1759 4718 : return FpX_sub(gel(Q, 1), gel(P, 1), p);
1760 0 : if (signe(gel(P,2))!=0) return pol_1(vT);
1761 0 : return FpXQ_inv(FpX_add(FpX_mulu(FpXQ_sqr(gel(P,1), T, p), 3, p),
1762 : a4, p), T, p);
1763 : }
1764 :
1765 : static GEN
1766 4655 : FpXQE_Miller_line(GEN R, GEN Q, GEN slope, GEN a4, GEN T, GEN p)
1767 : {
1768 4655 : long vT = get_FpX_var(T);
1769 4655 : GEN x = gel(Q, 1), y = gel(Q, 2);
1770 4655 : GEN tmp1 = FpX_sub(x, gel(R, 1), p);
1771 4655 : GEN tmp2 = FpX_add(FpXQ_mul(tmp1, slope, T, p), gel(R, 2), p);
1772 4655 : if (!ZX_equal(y, tmp2))
1773 4655 : return FpX_sub(y, tmp2, p);
1774 0 : if (signe(y) == 0)
1775 0 : return pol_1(vT);
1776 : else
1777 : {
1778 : GEN s1, s2;
1779 0 : GEN y2i = FpXQ_inv(FpX_mulu(y, 2, p), T, p);
1780 0 : s1 = FpXQ_mul(FpX_add(FpX_mulu(FpXQ_sqr(x, T, p), 3, p), a4, p), y2i, T, p);
1781 0 : if (!ZX_equal(s1, slope))
1782 0 : return FpX_sub(s1, slope, p);
1783 0 : s2 = FpXQ_mul(FpX_sub(FpX_mulu(x, 3, p), FpXQ_sqr(s1, T, p), p), y2i, T, p);
1784 0 : return signe(s2)!=0 ? s2: y2i;
1785 : }
1786 : }
1787 :
1788 : /* Computes the equation of the line tangent to R and returns its
1789 : evaluation at the point Q. Also doubles the point R.
1790 : */
1791 :
1792 : static GEN
1793 4158 : FpXQE_tangent_update(GEN R, GEN Q, GEN a4, GEN T, GEN p, GEN *pt_R)
1794 : {
1795 4158 : if (ell_is_inf(R))
1796 : {
1797 7 : *pt_R = ellinf();
1798 7 : return pol_1(get_FpX_var(T));
1799 : }
1800 4151 : else if (!signe(gel(R,2)))
1801 : {
1802 56 : *pt_R = ellinf();
1803 56 : return FpXQE_vert(R, Q, a4, T, p);
1804 : } else {
1805 : GEN slope;
1806 4095 : *pt_R = FpXQE_dbl_slope(R, a4, T, p, &slope);
1807 4095 : return FpXQE_Miller_line(R, Q, slope, a4, T, p);
1808 : }
1809 : }
1810 :
1811 : /* Computes the equation of the line through R and P, and returns its
1812 : evaluation at the point Q. Also adds P to the point R.
1813 : */
1814 :
1815 : static GEN
1816 567 : FpXQE_chord_update(GEN R, GEN P, GEN Q, GEN a4, GEN T, GEN p, GEN *pt_R)
1817 : {
1818 567 : if (ell_is_inf(R))
1819 : {
1820 0 : *pt_R = gcopy(P);
1821 0 : return FpXQE_vert(P, Q, a4, T, p);
1822 : }
1823 567 : else if (ell_is_inf(P))
1824 : {
1825 0 : *pt_R = gcopy(R);
1826 0 : return FpXQE_vert(R, Q, a4, T, p);
1827 : }
1828 567 : else if (ZX_equal(gel(P, 1), gel(R, 1)))
1829 : {
1830 7 : if (ZX_equal(gel(P, 2), gel(R, 2)))
1831 0 : return FpXQE_tangent_update(R, Q, a4, T, p, pt_R);
1832 : else
1833 : {
1834 7 : *pt_R = ellinf();
1835 7 : return FpXQE_vert(R, Q, a4, T, p);
1836 : }
1837 : } else {
1838 : GEN slope;
1839 560 : *pt_R = FpXQE_add_slope(P, R, a4, T, p, &slope);
1840 560 : return FpXQE_Miller_line(R, Q, slope, a4, T, p);
1841 : }
1842 : }
1843 :
1844 : struct _FpXQE_miller { GEN p, T, a4, P; };
1845 : static GEN
1846 4158 : FpXQE_Miller_dbl(void* E, GEN d)
1847 : {
1848 4158 : struct _FpXQE_miller *m = (struct _FpXQE_miller *)E;
1849 4158 : GEN p = m->p;
1850 4158 : GEN T = m->T, a4 = m->a4, P = m->P;
1851 : GEN v, line;
1852 4158 : GEN N = FpXQ_sqr(gel(d,1), T, p);
1853 4158 : GEN D = FpXQ_sqr(gel(d,2), T, p);
1854 4158 : GEN point = gel(d,3);
1855 4158 : line = FpXQE_tangent_update(point, P, a4, T, p, &point);
1856 4158 : N = FpXQ_mul(N, line, T, p);
1857 4158 : v = FpXQE_vert(point, P, a4, T, p);
1858 4158 : D = FpXQ_mul(D, v, T, p); return mkvec3(N, D, point);
1859 : }
1860 :
1861 : static GEN
1862 567 : FpXQE_Miller_add(void* E, GEN va, GEN vb)
1863 : {
1864 567 : struct _FpXQE_miller *m = (struct _FpXQE_miller *)E;
1865 567 : GEN p = m->p;
1866 567 : GEN T = m->T, a4 = m->a4, P = m->P;
1867 : GEN v, line, point;
1868 567 : GEN na = gel(va,1), da = gel(va,2), pa = gel(va,3);
1869 567 : GEN nb = gel(vb,1), db = gel(vb,2), pb = gel(vb,3);
1870 567 : GEN N = FpXQ_mul(na, nb, T, p);
1871 567 : GEN D = FpXQ_mul(da, db, T, p);
1872 567 : line = FpXQE_chord_update(pa, pb, P, a4, T, p, &point);
1873 567 : N = FpXQ_mul(N, line, T, p);
1874 567 : v = FpXQE_vert(point, P, a4, T, p);
1875 567 : D = FpXQ_mul(D, v, T, p); return mkvec3(N, D, point);
1876 : }
1877 :
1878 : /* Returns the Miller function f_{m, Q} evaluated at the point P using
1879 : * the standard Miller algorithm. */
1880 : static GEN
1881 63 : FpXQE_Miller(GEN Q, GEN P, GEN m, GEN a4, GEN T, GEN p)
1882 : {
1883 63 : pari_sp av = avma;
1884 : struct _FpXQE_miller d;
1885 : GEN v, N, D, g1;
1886 :
1887 63 : d.a4 = a4; d.T = T; d.p = p; d.P = P;
1888 63 : g1 = pol_1(get_FpX_var(T));
1889 63 : v = gen_pow_i(mkvec3(g1,g1,Q), m, (void*)&d,
1890 : FpXQE_Miller_dbl, FpXQE_Miller_add);
1891 63 : N = gel(v,1); D = gel(v,2);
1892 63 : return gc_upto(av, FpXQ_div(N, D, T, p));
1893 : }
1894 :
1895 : GEN
1896 28 : FpXQE_weilpairing(GEN P, GEN Q, GEN m, GEN a4, GEN T, GEN p)
1897 : {
1898 28 : pari_sp av = avma;
1899 : GEN N, D, w;
1900 28 : if (ell_is_inf(P) || ell_is_inf(Q) || ZXV_equal(P,Q))
1901 0 : return pol_1(get_FpX_var(T));
1902 28 : N = FpXQE_Miller(P, Q, m, a4, T, p);
1903 28 : D = FpXQE_Miller(Q, P, m, a4, T, p);
1904 28 : w = FpXQ_div(N, D, T, p);
1905 28 : if (mpodd(m)) w = FpX_neg(w, p);
1906 28 : return gc_upto(av, w);
1907 : }
1908 :
1909 : GEN
1910 7 : FpXQE_tatepairing(GEN P, GEN Q, GEN m, GEN a4, GEN T, GEN p)
1911 : {
1912 7 : if (ell_is_inf(P) || ell_is_inf(Q)) return pol_1(get_FpX_var(T));
1913 7 : return FpXQE_Miller(P, Q, m, a4, T, p);
1914 : }
1915 :
1916 : /***********************************************************************/
1917 : /** **/
1918 : /** issupersingular **/
1919 : /** **/
1920 : /***********************************************************************/
1921 :
1922 : GEN
1923 1718 : FpXQ_ellj(GEN a4, GEN a6, GEN T, GEN p)
1924 : {
1925 1718 : if (absequaliu(p,3)) return pol_0(get_FpX_var(T));
1926 : else
1927 : {
1928 1718 : pari_sp av=avma;
1929 1718 : GEN a43 = FpXQ_mul(a4,FpXQ_sqr(a4,T,p),T,p);
1930 1718 : GEN a62 = FpXQ_sqr(a6,T,p);
1931 1718 : GEN num = FpX_mulu(a43,6912,p);
1932 1718 : GEN den = FpX_add(FpX_mulu(a43,4,p),FpX_mulu(a62,27,p),p);
1933 1718 : return gc_uptoleaf(av, FpXQ_div(num, den, T, p));
1934 : }
1935 : }
1936 :
1937 : static GEN
1938 33530 : FpXQ_is_quad(GEN x, GEN T, GEN p)
1939 : {
1940 33530 : pari_sp av = avma;
1941 : GEN K;
1942 33530 : long d = degpol(T);
1943 33530 : x = FpXQ_red(x,T,p);
1944 33530 : if (lgpol(x)<=1) return NULL;
1945 33530 : if (d==2) return FpXQ_minpoly(x, T, p);
1946 33530 : if (odd(degpol(T))) return NULL;
1947 33530 : K = FpM_ker(FpXQ_matrix_pow(x, d, 3, T, p), p);
1948 33530 : if (lg(K)!=2) return gc_NULL(av);
1949 588 : return RgV_to_RgX(gel(K,1), get_FpX_var(T));
1950 : }
1951 :
1952 : int
1953 165515 : FpXQ_elljissupersingular(GEN j, GEN T, GEN p)
1954 : {
1955 165515 : pari_sp ltop = avma;
1956 :
1957 : /* All supersingular j-invariants are in FF_{p^2}, so we first check
1958 : * whether j is in FF_{p^2}. If d is odd, then FF_{p^2} is not a
1959 : * subfield of FF_{p^d} so the j-invariants are all in FF_p. Hence
1960 : * the j-invariants are in FF_{p^{2 - e}}. */
1961 165515 : ulong d = get_FpX_degree(T);
1962 : GEN S;
1963 165515 : if (degpol(j) <= 0) return Fp_elljissupersingular(constant_coeff(j), p);
1964 164612 : j = FpXQ_red(j, T, p);
1965 164612 : if (degpol(j) <= 0) return gc_bool(ltop, Fp_elljissupersingular(constant_coeff(j), p));
1966 : /* Now j is not in F_p */
1967 164612 : if (abscmpiu(p, 5) <= 0) return gc_bool(ltop,0); /* j != 0*/
1968 164605 : if (odd(d)) return 0;
1969 : /* Set S so that FF_p[T]/(S) is isomorphic to FF_{p^2}: */
1970 46949 : if (d == 2)
1971 13419 : S = T;
1972 : else /* d > 2 */
1973 : {
1974 33530 : S = FpXQ_is_quad(j, T, p);
1975 33530 : if (!S) return gc_bool(ltop,0);
1976 588 : j = pol_x(varn(S));
1977 : }
1978 14007 : return gc_bool(ltop, jissupersingular(j,S,p));
1979 : }
1980 :
1981 : int
1982 1050 : Fq_elljissupersingular(GEN j, GEN T, GEN p)
1983 959 : { return typ(j)==t_INT? Fp_elljissupersingular(j, p)
1984 2009 : : FpXQ_elljissupersingular(j, T, p); }
1985 :
1986 : /* p > 5 prime; return d such that (-d/p) = -1 */
1987 : static ulong
1988 1183 : find_inert_disc(GEN p)
1989 : {
1990 1183 : long s = mod4(p) == 1? -1: 1; /* - (-1/p) */
1991 1183 : ulong d = 3;
1992 : while(1)
1993 : {
1994 1190 : if (kroui(d,p) == s) return d; /* = 3 mod (16) */
1995 595 : d++;
1996 595 : if (kroui(d>>2,p) == s) return d; /* = 4 mod (16) */
1997 266 : d += 3;
1998 266 : if (kroui(d,p) == s) return d; /* = 7 mod (16) */
1999 105 : d++;
2000 105 : if (kroui(d>>2,p) == s) return d; /* = 8 mod (16) */
2001 35 : d += 3;
2002 35 : if (kroui(d,p) == s) return d; /* = 11 mod (16) */
2003 7 : d += 4;
2004 7 : if (kroui(d,p) == s) return d; /* = 15 mod (16) */
2005 7 : d += 4;
2006 : }
2007 : }
2008 :
2009 : /* p > 5 */
2010 : static GEN
2011 1183 : ellsupersingularj_easy_FpXQ(GEN T, GEN p)
2012 : {
2013 1183 : long d = find_inert_disc(p);
2014 1183 : GEN R = FpXQX_roots(polclass(stoi(-d), 0, 0), T, p);
2015 1183 : return gel(R,1);
2016 : }
2017 :
2018 : GEN
2019 1204 : ellsupersingularj_FpXQ(GEN T, GEN p)
2020 : {
2021 : GEN j, j2, R, Phi2;
2022 : long i, ep, lp;
2023 1204 : if (cmpiu(p, 5) <= 0) return pol_0(get_FpX_var(T));
2024 1183 : j2 = ellsupersingularj_easy_FpXQ(T, p);
2025 1183 : Phi2 = polmodular_ZXX(2,0,0,1);
2026 1183 : R = FpXQX_roots(FqXY_evalx(Phi2, j2, T, p), T, p);
2027 1183 : j = gel(R,1+random_Fl(lg(R)-1));
2028 1183 : ep = expi(p); lp = ep + random_Fl(ep);
2029 17849 : for (i = 1; i <= lp; i++)
2030 : {
2031 16666 : GEN Phi2_j = FqX_div_by_X_x(FqXY_evalx(Phi2, j, T, p), j2, T, p, NULL);
2032 16666 : R = FqX_quad_root(Phi2_j, T, p);
2033 16666 : if (!R) pari_err_PRIME("ellsupersingularj",p);
2034 16666 : j2 = j; j = random_bits(1) ? R: Fq_neg(Fq_add(gel(Phi2_j,3), R, T, p), T, p);
2035 : }
2036 1183 : return j;
2037 : }
2038 :
2039 : /***********************************************************************/
2040 : /** **/
2041 : /** Point counting **/
2042 : /** **/
2043 : /***********************************************************************/
2044 :
2045 : GEN
2046 15540 : elltrace_extension(GEN t, long n, GEN q)
2047 : {
2048 15540 : pari_sp av = avma;
2049 15540 : GEN v = RgX_to_RgC(RgXQ_powu(pol_x(0), n, mkpoln(3,gen_1,negi(t),q)),2);
2050 15540 : GEN te = addii(shifti(gel(v,1),1), mulii(t,gel(v,2)));
2051 15540 : return gc_INT(av, te);
2052 : }
2053 :
2054 : GEN
2055 14777 : Fp_ffellcard(GEN a4, GEN a6, GEN q, long n, GEN p)
2056 : {
2057 14777 : pari_sp av = avma;
2058 14777 : GEN ap = subii(addiu(p, 1), Fp_ellcard(a4, a6, p));
2059 14777 : GEN te = elltrace_extension(ap, n, p);
2060 14777 : return gc_INT(av, subii(addiu(q, 1), te));
2061 : }
2062 :
2063 : static GEN
2064 1687 : FpXQ_ellcardj(GEN a4, GEN a6, GEN j, GEN T, GEN q, GEN p, long n)
2065 : {
2066 1687 : GEN q1 = addiu(q,1);
2067 1687 : if (signe(j)==0)
2068 : {
2069 : GEN W, w, t, N;
2070 560 : if (umodiu(q,6)!=1) return q1;
2071 420 : N = Fp_ffellcard(gen_0,gen_1,q,n,p);
2072 420 : t = subii(q1, N);
2073 420 : W = FpXQ_pow(a6,diviuexact(shifti(q,-1), 3),T,p);
2074 420 : if (degpol(W)>0) /*p=5 mod 6*/
2075 105 : return ZX_equal1(FpXQ_powu(W,3,T,p)) ? addii(q1,shifti(t,-1)):
2076 35 : subii(q1,shifti(t,-1));
2077 350 : w = modii(gel(W,2),p);
2078 350 : if (equali1(w)) return N;
2079 259 : if (equalii(w,subiu(p,1))) return addii(q1,t);
2080 : else /*p=1 mod 6*/
2081 : {
2082 168 : GEN u = shifti(t,-1), v = sqrtint(diviuexact(subii(q,sqri(u)),3));
2083 168 : GEN a = addii(u,v), b = shifti(v,1);
2084 168 : if (equali1(Fp_powu(w,3,p)))
2085 : {
2086 84 : if (dvdii(addmulii(a, w, b), p))
2087 56 : return subii(q1,subii(shifti(b,1),a));
2088 : else
2089 28 : return addii(q1,addii(a,b));
2090 : }
2091 : else
2092 : {
2093 84 : if (dvdii(submulii(a, w, b), p))
2094 56 : return subii(q1,subii(a,shifti(b,1)));
2095 : else
2096 28 : return subii(q1,addii(a,b));
2097 : }
2098 : }
2099 1127 : } else if (equalii(j,modsi(1728,p)))
2100 : {
2101 : GEN w, W, N, t;
2102 567 : if (mod4(q)==3) return q1;
2103 427 : W = FpXQ_pow(a4,shifti(q,-2),T,p);
2104 427 : if (degpol(W)>0) return q1; /*p=3 mod 4*/
2105 357 : w = modii(gel(W,2),p);
2106 357 : N = Fp_ffellcard(gen_1,gen_0,q,n,p);
2107 357 : if (equali1(w)) return N;
2108 238 : t = subii(q1, N);
2109 238 : if (equalii(w,subiu(p,1))) return addii(q1,t);
2110 : else /*p=1 mod 4*/
2111 : {
2112 112 : GEN u = shifti(t,-1), v = sqrtint(subii(q,sqri(u)));
2113 112 : if (dvdii(addmulii(u, w, v), p))
2114 56 : return subii(q1,shifti(v,1));
2115 : else
2116 56 : return addii(q1,shifti(v,1));
2117 : }
2118 : } else
2119 : {
2120 560 : GEN g = Fp_div(j, Fp_sub(utoi(1728), j, p), p);
2121 560 : GEN l = FpXQ_div(FpX_mulu(a6,3,p),FpX_mulu(a4,2,p),T,p);
2122 560 : GEN N = Fp_ffellcard(Fp_mulu(g,3,p),Fp_double(g,p),q,n,p);
2123 560 : if (FpXQ_issquare(l,T,p)) return N;
2124 280 : return subii(shifti(q1,1),N);
2125 : }
2126 : }
2127 :
2128 : static GEN
2129 8 : FpXQ_ffellcard(GEN a4, GEN a6, GEN M, GEN q, GEN T, GEN p, long n)
2130 : {
2131 8 : long m = degpol(M);
2132 8 : GEN j = pol_x(get_FpX_var(T));
2133 8 : GEN g = FpXQ_div(j, Fp_FpX_sub(utoi(1728), j, p), M, p);
2134 8 : GEN N = FpXQ_ellcard(FpX_mulu(g,3,p),FpX_mulu(g,2,p),M,p);
2135 8 : GEN qm = powiu(p, m), q1 = addiu(q, 1), qm1 = addiu(qm, 1);
2136 8 : GEN l = FpXQ_mul(FpX_mulu(a6,3,p),FpX_mulu(a4,2,p),T,p);
2137 8 : GEN te = elltrace_extension(subii(qm1, N), n/m, qm);
2138 8 : return FpXQ_issquare(l,T,p) ? subii(q1, te): addii(q1, te);
2139 : }
2140 :
2141 : static int
2142 7 : FpXQ_is4power(GEN x, GEN T, GEN p)
2143 : {
2144 7 : long d = get_FpX_degree(T);
2145 7 : if (lg(x) == 2 || absequalui(2, p)) return 1;
2146 7 : if (Mod4(p)==1)
2147 7 : return equali1(Fp_pow(FpXQ_norm(x,T,p),shifti(p,-2), p));
2148 0 : if (odd(d))
2149 0 : return FpXQ_issquare(x, T, p);
2150 0 : return ZX_equal1(FpXQ_pow(x, shifti(powiu(p, d),-2), T, p));
2151 : }
2152 :
2153 : /* http://www.numdam.org/article/ASENS_1969_4_2_4_521_0.pdf */
2154 :
2155 : GEN
2156 7 : FpXQ_ellcard_supersingular(GEN a4, GEN a6, GEN T, GEN p)
2157 : {
2158 7 : pari_sp av = avma;
2159 7 : long d = get_FpX_degree(T);
2160 : GEN r;
2161 7 : if (equaliu(p,3))
2162 0 : r = Flxq_ellcard(ZX_to_Flx(a4,3), ZX_to_Flx(a6,3), ZXT_to_FlxT(T,3), 3);
2163 7 : else if (signe(a4)==0)
2164 0 : r = FpXQ_ellcardj(a4, a6, gen_0, T, powiu(p, d), p, d);
2165 7 : else if (signe(a6)==0)
2166 0 : r = FpXQ_ellcardj(a4, a6, modsi(1728,p), T, powiu(p, d), p, d);
2167 : else
2168 : {
2169 : GEN q, q2, t, D;
2170 7 : long qm4 = (odd(d>>1) && Mod4(p)==3);
2171 7 : if (odd(d)) return gen_0;
2172 7 : q2 = powiu(p, d>>1); q = sqri(q2);
2173 7 : t = shifti(q2, 1);
2174 7 : D = FpX_sub(FpX_Fp_mul(FpXQ_powu(a4,3,T,p), stoi(-4), p),
2175 : FpX_mulu(FpXQ_sqr(a6,T,p), 27, p), p);
2176 14 : r = qm4 ^ FpXQ_is4power(D, T, p) ? subii(addiu(q, 1), t)
2177 7 : : addii(addiu(q, 1), t);
2178 : }
2179 7 : return gc_INT(av, r);
2180 : }
2181 :
2182 : GEN
2183 21 : Fq_ellcard_supersingular(GEN a4, GEN a6, GEN T, GEN p)
2184 21 : { return T ? FpXQ_ellcard_supersingular(a4, a6, T, p) : addiu(p, 1); }
2185 :
2186 : static GEN
2187 8578 : FpXQ_ellcard_i(GEN a4, GEN a6, GEN T, GEN p)
2188 : {
2189 8578 : long n = get_FpX_degree(T);
2190 8578 : GEN q = powiu(p, n);
2191 8578 : if (degpol(a4)<=0 && degpol(a6)<=0)
2192 833 : return Fp_ffellcard(constant_coeff(a4),constant_coeff(a6),q,n,p);
2193 7745 : if (lgefint(p)==3)
2194 : {
2195 6027 : ulong pp = p[2];
2196 6027 : return Flxq_ellcard(ZX_to_Flx(a4,pp),ZX_to_Flx(a6,pp),ZX_to_Flx(T,pp),pp);
2197 : }
2198 : else
2199 : {
2200 1718 : GEN J = FpXQ_ellj(a4,a6,T,p), M;
2201 1718 : if (degpol(J) <= 0)
2202 1687 : return FpXQ_ellcardj(a4,a6,constant_coeff(J),T,q,p,n);
2203 31 : M = FpXQ_minpoly(J,T,p);
2204 31 : if (degpol(M) < degpol(T))
2205 8 : return FpXQ_ffellcard(a4, a6, M, q, T, p, n);
2206 23 : return Fq_ellcard_SEA(a4, a6, q, T, p, 0);
2207 : }
2208 : }
2209 :
2210 : GEN
2211 8578 : FpXQ_ellcard(GEN a4, GEN a6, GEN T, GEN p)
2212 : {
2213 8578 : pari_sp av = avma;
2214 8578 : return gc_INT(av, FpXQ_ellcard_i(a4, a6, T, p));
2215 : }
2216 :
2217 : static GEN
2218 21 : _FpXQE_pairorder(void *E, GEN P, GEN Q, GEN m, GEN F)
2219 : {
2220 21 : struct _FpXQE *e = (struct _FpXQE *) E;
2221 21 : return FpXQ_order(FpXQE_weilpairing(P,Q,m,e->a4,e->T,e->p), F, e->T, e->p);
2222 : }
2223 :
2224 : GEN
2225 15 : FpXQ_ellgroup(GEN a4, GEN a6, GEN N, GEN T, GEN p, GEN *pt_m)
2226 : {
2227 : struct _FpXQE e;
2228 15 : GEN q = powiu(p, get_FpX_degree(T));
2229 15 : e.a4=a4; e.a6=a6; e.T=T; e.p=p;
2230 15 : return gen_ellgroup(N, subiu(q,1), pt_m, (void*)&e, &FpXQE_group, _FpXQE_pairorder);
2231 : }
2232 :
2233 : GEN
2234 8 : FpXQ_ellgens(GEN a4, GEN a6, GEN ch, GEN D, GEN m, GEN T, GEN p)
2235 : {
2236 : GEN P;
2237 8 : pari_sp av = avma;
2238 : struct _FpXQE e;
2239 8 : e.a4=a4; e.a6=a6; e.T=T; e.p=p;
2240 8 : switch(lg(D)-1)
2241 : {
2242 8 : case 1:
2243 8 : P = gen_gener(gel(D,1), (void*)&e, &FpXQE_group);
2244 8 : P = mkvec(FpXQE_changepoint(P, ch, T, p));
2245 8 : break;
2246 0 : default:
2247 0 : P = gen_ellgens(gel(D,1), gel(D,2), m, (void*)&e, &FpXQE_group, _FpXQE_pairorder);
2248 0 : gel(P,1) = FpXQE_changepoint(gel(P,1), ch, T, p);
2249 0 : gel(P,2) = FpXQE_changepoint(gel(P,2), ch, T, p);
2250 0 : break;
2251 : }
2252 8 : return gc_GEN(av, P);
2253 : }
|
