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