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 38273 : ellinf_FpJ(void)
35 38273 : { 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 1711518 : 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 1711518 : if (signe(gel(P,3)) == 0) return ellinf_FpJ();
46 :
47 1704060 : X1 = gel(P,1); Y1 = gel(P,2); Z1 = gel(P,3);
48 :
49 1704060 : XX = Fp_sqr(X1, p);
50 1704060 : YY = Fp_sqr(Y1, p);
51 1704060 : YYYY = Fp_sqr(YY, p);
52 1704060 : ZZ = Fp_sqr(Z1, p);
53 1704060 : S = Fp_double(Fp_sub(Fp_sqr(Fp_add(X1,YY,p), p), Fp_add(XX,YYYY,p), p), p);
54 1704060 : M = Fp_addmul(Fp_mulu(XX, 3, p), a4, Fp_sqr(ZZ, p), p);
55 1704060 : T = Fp_sub(Fp_sqr(M, p), Fp_double(S, p), p);
56 1704060 : Q = cgetg(4, t_VEC);
57 1704060 : gel(Q,1) = T;
58 1704060 : gel(Q,2) = Fp_sub(Fp_mul(M, Fp_sub(S, T, p), p), Fp_mulu(YYYY, 8, p), p);
59 1704060 : gel(Q,3) = Fp_sub(Fp_sqr(Fp_add(Y1, Z1, p), p), Fp_add(YY, ZZ, p), p);
60 1704060 : 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 510907 : 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 510907 : if (signe(gel(Q,3)) == 0) return gcopy(P);
72 510907 : if (signe(gel(P,3)) == 0) return gcopy(Q);
73 :
74 493657 : X1 = gel(P,1); Y1 = gel(P,2); Z1 = gel(P,3);
75 493657 : X2 = gel(Q,1); Y2 = gel(Q,2); Z2 = gel(Q,3);
76 :
77 493657 : Z1Z1 = Fp_sqr(Z1, p);
78 493657 : Z2Z2 = Fp_sqr(Z2, p);
79 493657 : U1 = Fp_mul(X1, Z2Z2, p);
80 493657 : U2 = Fp_mul(X2, Z1Z1, p);
81 493657 : S1 = mulii(Y1, Fp_mul(Z2, Z2Z2, p));
82 493657 : S2 = mulii(Y2, Fp_mul(Z1, Z1Z1, p));
83 493657 : H = Fp_sub(U2, U1, p);
84 493657 : r = Fp_double(Fp_sub(S2, S1, p), p);
85 :
86 : /* If points are equal we must double. */
87 493657 : if (signe(H)== 0) {
88 31709 : if (signe(r) == 0)
89 : /* Points are equal so double. */
90 894 : return FpJ_dbl(P, a4, p);
91 : else
92 30815 : return ellinf_FpJ();
93 : }
94 461948 : I = Fp_sqr(Fp_double(H, p), p);
95 461948 : J = Fp_mul(H, I, p);
96 461948 : V = Fp_mul(U1, I, p);
97 461948 : W = Fp_sub(Fp_sqr(r, p), Fp_add(J, Fp_double(V, p), p), p);
98 461948 : R = cgetg(4, t_VEC);
99 461948 : gel(R,1) = W;
100 461948 : gel(R,2) = Fp_sub(mulii(r, subii(V, W)),
101 : shifti(mulii(S1, J), 1), p);
102 461948 : 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 461948 : 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 169828 : FpE_to_FpJ(GEN P)
115 : {
116 169828 : return ell_is_inf(P) ? ellinf_FpJ()
117 169828 : : mkvec3(icopy(gel(P,1)),icopy(gel(P,2)), gen_1);
118 : }
119 :
120 : GEN
121 169396 : FpJ_to_FpE(GEN P, GEN p)
122 : {
123 169396 : if (signe(gel(P,3)) == 0) return ellinf();
124 : else
125 : {
126 138932 : GEN Z = Fp_inv(gel(P,3), p);
127 138920 : GEN Z2 = Fp_sqr(Z, p), Z3 = Fp_mul(Z, Z2, p);
128 138920 : 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 1710624 : _FpJ_dbl(void *E, GEN P)
135 : {
136 1710624 : struct _FpE *ell = (struct _FpE *) E;
137 1710624 : return FpJ_dbl(P, ell->a4, ell->p);
138 : }
139 : static GEN
140 510907 : _FpJ_add(void *E, GEN P, GEN Q)
141 : {
142 510907 : struct _FpE *ell=(struct _FpE *) E;
143 510907 : return FpJ_add(P, Q, ell->a4, ell->p);
144 : }
145 : static GEN
146 4908 : _FpJ_mul(void *E, GEN P, GEN n)
147 : {
148 4908 : pari_sp av = avma;
149 4908 : struct _FpE *e=(struct _FpE *) E;
150 4908 : long s = signe(n);
151 4908 : if (!s || signe(gel(P,3))==0) return ellinf_FpJ();
152 4908 : if (s < 0) P = FpJ_neg(P, e->p);
153 4908 : if (is_pm1(n)) return s > 0 ? gcopy(P): P;
154 4908 : return gc_GEN(av, gen_pow_i(P, n, e, &_FpJ_dbl, &_FpJ_add));
155 : }
156 :
157 : GEN
158 4908 : FpJ_mul(GEN P, GEN n, GEN a4, GEN p)
159 : {
160 : struct _FpE E;
161 4908 : E.a4= a4; E.p = p;
162 4908 : 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 1530 : RgE_to_FpE(GEN x, GEN p)
177 : {
178 1530 : if (ell_is_inf(x)) return x;
179 1530 : retmkvec2(Rg_to_Fp(gel(x,1),p),Rg_to_Fp(gel(x,2),p));
180 : }
181 :
182 : GEN
183 126 : FpE_to_mod(GEN x, GEN p)
184 : {
185 126 : if (ell_is_inf(x)) return x;
186 72 : retmkvec2(Fp_to_mod(gel(x,1),p),Fp_to_mod(gel(x,2),p));
187 : }
188 :
189 : GEN
190 702 : FpE_changepoint(GEN P, GEN ch, GEN p)
191 : {
192 702 : pari_sp av = avma;
193 : GEN c, z, u, r, s, t, v, v2, v3;
194 702 : if (ell_is_inf(P)) return P;
195 648 : if (lgefint(p) == 3)
196 : {
197 611 : ulong pp = p[2];
198 611 : z = Fle_changepoint(ZV_to_Flv(P, pp), ZV_to_Flv(ch, pp), pp);
199 611 : return gc_upto(av, Flv_to_ZV(z));
200 : }
201 37 : u = gel(ch,1); r = gel(ch,2); s = gel(ch,3); t = gel(ch,4);
202 37 : v = Fp_inv(u, p); v2 = Fp_sqr(v,p); v3 = Fp_mul(v,v2,p);
203 37 : c = Fp_sub(gel(P,1),r,p);
204 37 : z = cgetg(3,t_VEC);
205 37 : gel(z,1) = Fp_mul(v2, c, p);
206 37 : gel(z,2) = Fp_mul(v3, Fp_sub(gel(P,2), Fp_add(Fp_mul(s,c, p),t, p),p),p);
207 37 : return gc_upto(av, z);
208 : }
209 :
210 : GEN
211 1530 : FpE_changepointinv(GEN P, GEN ch, GEN p)
212 : {
213 : GEN u, r, s, t, u2, u3, c, z;
214 1530 : if (ell_is_inf(P)) return P;
215 1530 : if (lgefint(p) == 3)
216 : {
217 1485 : ulong pp = p[2];
218 1485 : z = Fle_changepointinv(ZV_to_Flv(P, pp), ZV_to_Flv(ch, pp), pp);
219 1485 : return Flv_to_ZV(z);
220 : }
221 45 : u = gel(ch,1); r = gel(ch,2); s = gel(ch,3); t = gel(ch,4);
222 45 : u2 = Fp_sqr(u, p); u3 = Fp_mul(u,u2,p);
223 45 : c = Fp_mul(u2, gel(P,1), p);
224 45 : z = cgetg(3, t_VEC);
225 45 : gel(z,1) = Fp_add(c,r,p);
226 45 : gel(z,2) = Fp_add(Fp_mul(u3,gel(P,2),p), Fp_add(Fp_mul(s,c,p), t, p), p);
227 45 : return z;
228 : }
229 :
230 : static GEN
231 360 : random_nonsquare_Fp(GEN p)
232 : {
233 360 : pari_sp av = avma;
234 : GEN a;
235 360 : switch(mod8(p))
236 : { /* easy special cases */
237 360 : 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 222330 : FpE_dbl_slope(GEN P, GEN a4, GEN p, GEN *slope)
258 : {
259 : GEN x, y, Q;
260 222330 : if (ell_is_inf(P) || !signe(gel(P,2))) return ellinf();
261 113934 : x = gel(P,1); y = gel(P,2);
262 113934 : *slope = Fp_div(Fp_add(Fp_mulu(Fp_sqr(x,p), 3, p), a4, p),
263 : Fp_mulu(y, 2, p), p);
264 113934 : Q = cgetg(3,t_VEC);
265 113934 : gel(Q, 1) = Fp_sub(Fp_sqr(*slope, p), Fp_mulu(x, 2, p), p);
266 113934 : gel(Q, 2) = Fp_sub(Fp_mul(*slope, Fp_sub(x, gel(Q, 1), p), p), y, p);
267 113934 : return Q;
268 : }
269 :
270 : GEN
271 221804 : FpE_dbl(GEN P, GEN a4, GEN p)
272 : {
273 221804 : pari_sp av = avma;
274 : GEN slope;
275 221804 : return gc_upto(av, FpE_dbl_slope(P,a4,p,&slope));
276 : }
277 :
278 : static GEN
279 788427 : FpE_add_slope(GEN P, GEN Q, GEN a4, GEN p, GEN *slope)
280 : {
281 : GEN Px, Py, Qx, Qy, R;
282 788427 : if (ell_is_inf(P)) return Q;
283 788007 : if (ell_is_inf(Q)) return P;
284 788007 : Px = gel(P,1); Py = gel(P,2);
285 788007 : Qx = gel(Q,1); Qy = gel(Q,2);
286 788007 : if (equalii(Px, Qx))
287 : {
288 503 : if (equalii(Py, Qy))
289 485 : return FpE_dbl_slope(P, a4, p, slope);
290 : else
291 18 : return ellinf();
292 : }
293 787504 : *slope = Fp_div(Fp_sub(Py, Qy, p), Fp_sub(Px, Qx, p), p);
294 787504 : R = cgetg(3,t_VEC);
295 787504 : gel(R, 1) = Fp_sub(Fp_sub(Fp_sqr(*slope, p), Px, p), Qx, p);
296 787504 : gel(R, 2) = Fp_sub(Fp_mul(*slope, Fp_sub(Px, gel(R, 1), p), p), Py, p);
297 787504 : return R;
298 : }
299 :
300 : GEN
301 788423 : FpE_add(GEN P, GEN Q, GEN a4, GEN p)
302 : {
303 788423 : pari_sp av = avma;
304 : GEN slope;
305 788423 : 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 310717 : FpE_neg(GEN P, GEN p)
317 : {
318 310717 : if (ell_is_inf(P)) return ellinf();
319 310717 : 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 221804 : _FpE_dbl(void *E, GEN P)
332 : {
333 221804 : struct _FpE *ell = (struct _FpE *) E;
334 221804 : return FpE_dbl(P, ell->a4, ell->p);
335 : }
336 :
337 : static GEN
338 769152 : _FpE_add(void *E, GEN P, GEN Q)
339 : {
340 769152 : struct _FpE *ell=(struct _FpE *) E;
341 769152 : return FpE_add(P, Q, ell->a4, ell->p);
342 : }
343 :
344 : static GEN
345 763038 : _FpE_mul(void *E, GEN P, GEN n)
346 : {
347 763038 : pari_sp av = avma;
348 763038 : struct _FpE *e=(struct _FpE *) E;
349 763038 : long s = signe(n);
350 : GEN Q;
351 763038 : if (!s || ell_is_inf(P)) return ellinf();
352 763032 : if (s<0) P = FpE_neg(P, e->p);
353 763032 : if (is_pm1(n)) return s>0? gcopy(P): P;
354 391200 : if (equalis(n,2)) return _FpE_dbl(E, P);
355 169396 : Q = gen_pow_i(FpE_to_FpJ(P), n, e, &_FpJ_dbl, &_FpJ_add);
356 169396 : return gc_upto(av, FpJ_to_FpE(Q, e->p));
357 : }
358 :
359 : GEN
360 349 : FpE_mul(GEN P, GEN n, GEN a4, GEN p)
361 : {
362 : struct _FpE E;
363 349 : E.a4 = a4; E.p = p;
364 349 : return _FpE_mul(&E, P, n);
365 : }
366 :
367 : /* Finds a random nonsingular point on E */
368 :
369 : GEN
370 161711 : random_FpE(GEN a4, GEN a6, GEN p)
371 : {
372 161711 : pari_sp ltop = avma;
373 : GEN x, x2, y, rhs;
374 : do
375 : {
376 282010 : set_avma(ltop);
377 282010 : x = randomi(p); /* x^3+a4*x+a6 = x*(x^2+a4)+a6 */
378 282010 : x2 = Fp_sqr(x, p);
379 282010 : rhs = Fp_add(Fp_mul(x, Fp_add(x2, a4, p), p), a6, p);
380 30126 : } while ((!signe(rhs) && !signe(Fp_add(Fp_mulu(x2,3,p),a4,p)))
381 312136 : || kronecker(rhs, p) < 0);
382 161711 : y = Fp_sqrt(rhs, p);
383 161711 : if (!y) pari_err_PRIME("random_FpE", p);
384 161711 : return gc_GEN(ltop, mkvec2(x, y));
385 : }
386 :
387 : static GEN
388 159665 : _FpE_rand(void *E)
389 : {
390 159665 : struct _FpE *e=(struct _FpE *) E;
391 159665 : 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 774 : get_FpE_group(void ** pt_E, GEN a4, GEN a6, GEN p)
398 : {
399 774 : struct _FpE *e = (struct _FpE *) stack_malloc(sizeof(struct _FpE));
400 774 : e->a4 = a4; e->a6 = a6; e->p = p;
401 774 : *pt_E = (void *) e;
402 774 : return &FpE_group;
403 : }
404 :
405 : GEN
406 643 : FpE_order(GEN z, GEN o, GEN a4, GEN p)
407 : {
408 643 : pari_sp av = avma;
409 : struct _FpE e;
410 : GEN r;
411 643 : if (lgefint(p) == 3)
412 : {
413 537 : ulong pp = p[2];
414 537 : 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 643 : return gc_INT(av, r);
423 : }
424 :
425 : GEN
426 42 : FpE_log(GEN a, GEN b, GEN o, GEN a4, GEN p)
427 : {
428 42 : pari_sp av = avma;
429 : struct _FpE e;
430 : GEN r;
431 42 : if (lgefint(p) == 3)
432 : {
433 42 : ulong pp = p[2];
434 42 : 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 42 : 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 62544 : FpE_weilpairing(GEN P, GEN Q, GEN m, GEN a4, GEN p)
590 : {
591 62544 : pari_sp av = avma;
592 : GEN N, D, w;
593 62544 : if (ell_is_inf(P) || ell_is_inf(Q) || ZV_equal(P,Q)) return gen_1;
594 41418 : if (lgefint(p)==3 && lgefint(m)==3)
595 : {
596 41396 : ulong pp = p[2];
597 41396 : GEN Pp = ZV_to_Flv(P, pp), Qp = ZV_to_Flv(Q, pp);
598 41396 : ulong w = Fle_weilpairing(Pp, Qp, itou(m), umodiu(a4, pp), pp);
599 41396 : 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 174 : FpE_tatepairing(GEN P, GEN Q, GEN m, GEN a4, GEN p)
610 : {
611 174 : if (ell_is_inf(P) || ell_is_inf(Q)) return gen_1;
612 174 : if (lgefint(p)==3 && lgefint(m)==3)
613 : {
614 174 : pari_sp av = avma;
615 174 : ulong pp = p[2];
616 174 : GEN Pp = ZV_to_Flv(P, pp), Qp = ZV_to_Flv(Q, pp);
617 174 : ulong w = Fle_tatepairing(Pp, Qp, itou(m), umodiu(a4, pp), pp);
618 174 : 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 564509 : is_CMj(long J, GEN jn, GEN jd, GEN p)
632 564509 : { 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 14579 : u2_is_CMj(ulong a, ulong b, GEN jn, GEN jd, GEN p)
637 : {
638 14579 : GEN mJ = uu32toi(a,b);
639 14579 : return dvdii(addii(mulii(jd,mJ), jn), p);
640 : }
641 : #endif
642 :
643 : static long
644 45513 : 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 45513 : CHECK(-3, 0);
648 45421 : CHECK(-4, 1728);
649 45321 : CHECK(-7, -3375);
650 45114 : CHECK(-8, 8000);
651 44931 : CHECK(-11, -32768);
652 44724 : CHECK(-12, 54000);
653 44526 : CHECK(-16, 287496);
654 44370 : CHECK(-19, -884736);
655 44184 : CHECK(-27, -12288000);
656 43992 : CHECK(-28, 16581375);
657 43818 : CHECK(-43, -884736000);
658 : #ifdef LONG_IS_64BIT
659 36350 : CHECK(-67, -147197952000L);
660 36245 : CHECK(-163, -262537412640768000L);
661 : #else
662 7300 : if (u2_is_CMj(0x00000022UL,0x45ae8000UL,jn,jd,p)) return -67;
663 7279 : if (u2_is_CMj(0x03a4b862UL,0xc4b40000UL,jn,jd,p)) return -163;
664 : #endif
665 : #undef CHECK
666 43363 : 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 59212 : FqX_quad_root(GEN x, GEN T, GEN p)
678 : {
679 59212 : GEN b = gel(x,3), c = gel(x,2);
680 59212 : GEN D = Fq_sub(Fq_sqr(b, T, p), Fq_mulu(c,4, T, p), T, p);
681 59212 : GEN s = Fq_sqrt(D,T, p);
682 59212 : if (!s) return NULL;
683 56502 : return Fq_halve(Fq_sub(s, b, T, p), T, p);
684 : }
685 :
686 : static GEN
687 954 : FpX_quad_root(GEN x, GEN p)
688 : {
689 954 : GEN s, b = gel(x,3), c = gel(x,2);
690 954 : GEN D = Fp_sub(Fp_sqr(b, p), shifti(c,2), p);
691 954 : if (kronecker(D,p) == -1) return NULL;
692 601 : s = Fp_sqrt(D,p);
693 601 : 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 4037 : Fq_path_extends_to_floor(GEN j_prev, GEN j, GEN T, GEN p, GEN Phi2, long max_len)
701 : {
702 4037 : pari_sp ltop = avma;
703 4037 : 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 25559 : for (d = 1; d <= max_len; ++d)
709 : {
710 65050 : for (i = 1; i < l; i++)
711 : {
712 43528 : GEN Phi2_j = FqX_div_by_X_x(FqXY_evalx(Phi2, gel(j,i), T, p), gel(j_prev,i), T, p, NULL);
713 43528 : GEN j_next = FqX_quad_root(Phi2_j, T, p);
714 43528 : if (!j_next)
715 2710 : return gc_long(ltop, 1);
716 40818 : gel(j_prev,i) = gel(j, i); gel(j,i) = j_next;
717 : }
718 21522 : if (gc_needed(ltop, 2))
719 0 : (void)gc_all(ltop, 2, &j, &j_prev);
720 : }
721 1327 : return gc_long(ltop, 0);
722 : }
723 :
724 : static long
725 353 : 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 353 : pari_sp ltop = avma;
728 353 : 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 482 : for (d = 1; d <= max_len; ++d)
734 : {
735 1083 : for (i = 1; i < l; i++)
736 : {
737 954 : GEN Phi2_j = FpX_div_by_X_x(FpXY_evalx(Phi2, gel(j,i), p), gel(j_prev,i), p, NULL);
738 954 : GEN j_next = FpX_quad_root(Phi2_j, p);
739 954 : if (!j_next)
740 : {
741 353 : *pt_j = gel(j,i);
742 353 : *pt_j_prev = gel(j_prev,i);
743 353 : return 1;
744 : }
745 601 : gel(j_prev,i) = gel(j, i); gel(j,i) = j_next;
746 : }
747 129 : 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 2245 : Fp_jissupersingular(GEN j, GEN p)
756 : {
757 2245 : long max_path_len = expi(p)+1;
758 2245 : GEN Phi2 = FpXX_red(polmodular_ZXX(2,0,0,1), p);
759 2245 : GEN Phi2_j = FpXY_evalx(Phi2, j, p);
760 2245 : GEN roots = FpX_roots(Phi2_j, p);
761 2245 : long nbroots = lg(roots)-1;
762 2245 : 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 2245 : if (nbroots==0)
767 566 : return 0;
768 1679 : S = deg2pol_shallow(gen_1, gen_0, Fp_neg(Fp_2gener(p),p),1);
769 1679 : if (nbroots==1 && FpX_is_squarefree(Phi2_j, p))
770 1326 : { j_prev = j; j = FqX_quad_root(FpX_div_by_X_x(Phi2_j, gel(roots,1), p, NULL), S, p); }
771 : else
772 353 : if (!Fp_path_extends_to_floor(const_vec(nbroots,j), roots, p, Phi2, max_path_len, &j, &j_prev))
773 0 : return 1;
774 1679 : return !Fq_path_extends_to_floor(mkvec(j_prev), mkvec(j), S, p, Phi2, max_path_len);
775 : }
776 :
777 : static int
778 12000 : jissupersingular(GEN j, GEN S, GEN p)
779 : {
780 12000 : long max_path_len = expi(p)+1;
781 12000 : GEN Phi2 = FpXX_red(polmodular_ZXX(2,0,0,1), p);
782 12000 : GEN Phi2_j = FqXY_evalx(Phi2, j, S, p);
783 12000 : GEN roots = FpXQX_roots(Phi2_j, S, p);
784 12000 : 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 12000 : if (nbroots==0 || (nbroots==1 && FqX_is_squarefree(Phi2_j, S, p)))
789 9642 : return 0;
790 : else
791 2358 : return !Fq_path_extends_to_floor(const_vec(nbroots,j), roots, S, p, Phi2, max_path_len);
792 : }
793 :
794 : int
795 3081 : Fp_elljissupersingular(GEN j, GEN p)
796 : {
797 : long CM;
798 3081 : if (abscmpiu(p, 5) <= 0) return signe(j) == 0; /* valid if p <= 5 */
799 2961 : CM = Fp_ellj_get_CM(j, gen_1, p);
800 2961 : if (CM < 0) return krosi(CM, p) < 0; /* valid if p > 3 */
801 : else
802 2245 : 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 656776 : Fl_elltrace_naive(ulong a4, ulong a6, ulong p)
814 : {
815 : ulong i, j;
816 656776 : long a = 0;
817 : long d0, d1, d2, d3;
818 656776 : GEN k = const_vecsmall(p, -1);
819 656776 : k[1] = 0;
820 119129170 : for (i=1, j=1; i < p; i += 2, j = Fl_add(j, i, p))
821 118472394 : k[j+1] = 1;
822 656776 : d0 = 6%p; d1 = d0; d2 = Fl_add(a4, 1, p); d3 = a6;
823 656776 : for(i=0;; i++)
824 : {
825 237601564 : a -= k[1+d3];
826 237601564 : if (i==p-1) break;
827 236944788 : d3 = Fl_add(d3, d2, p);
828 236944788 : d2 = Fl_add(d2, d1, p);
829 236944788 : d1 = Fl_add(d1, d0, p);
830 : }
831 656776 : 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 218941 : closest_lift(GEN a, GEN b, GEN c)
865 : {
866 218941 : 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 6228 : Fl_ellpoint(long KRO, ulong *px, ulong c4, ulong c6, ulong p)
898 : {
899 6228 : ulong t, u, x = *px;
900 : for(;;)
901 : {
902 12241 : if (++x >= p) pari_err_PRIME("ellap",utoi(p));
903 12241 : t = Fl_add(c4, Fl_sqr(x,p), p);
904 12241 : u = Fl_add(c6, Fl_mul(x, t, p), p);
905 12241 : if (krouu(u,p) == KRO) break;
906 : }
907 6228 : *px = x;
908 6228 : 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 12916581 : compare_multiples(const void *A, const void *B)
1144 : {
1145 12916581 : multiple *a = (multiple*)A, *b = (multiple*)B;
1146 12916581 : 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 225011 : uclosest_lift(ulong a, ulong b, ulong c)
1154 : {
1155 225011 : ulong x = (b + ((c-a) << 1)) / (b << 1);
1156 225011 : return a + b * x;
1157 : }
1158 :
1159 : static long
1160 194716 : Fle_dbl_inplace(GEN P, ulong a4, ulong p)
1161 : {
1162 : ulong x, y, slope;
1163 194716 : if (!P[2]) return 1;
1164 194692 : x = P[1]; y = P[2];
1165 194692 : slope = Fl_div(Fl_add(Fl_triple(Fl_sqr(x,p), p), a4, p),
1166 : Fl_double(y, p), p);
1167 194692 : P[1] = Fl_sub(Fl_sqr(slope, p), Fl_double(x, p), p);
1168 194692 : P[2] = Fl_sub(Fl_mul(slope, Fl_sub(x, P[1], p), p), y, p);
1169 194692 : return 0;
1170 : }
1171 :
1172 : static long
1173 4921807 : Fle_add_inplace(GEN P, GEN Q, ulong a4, ulong p)
1174 : {
1175 : ulong Px, Py, Qx, Qy, slope;
1176 4921807 : if (ell_is_inf(Q)) return 0;
1177 4921807 : Px = P[1]; Py = P[2];
1178 4921807 : Qx = Q[1]; Qy = Q[2];
1179 4921807 : if (Px==Qx)
1180 204539 : return Py==Qy ? Fle_dbl_inplace(P, a4, p): 1;
1181 4717268 : slope = Fl_div(Fl_sub(Py, Qy, p), Fl_sub(Px, Qx, p), p);
1182 4717268 : P[1] = Fl_sub(Fl_sub(Fl_sqr(slope, p), Px, p), Qx, p);
1183 4717268 : P[2] = Fl_sub(Fl_mul(slope, Fl_sub(Px, P[1], p), p), Py, p);
1184 4717268 : 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 218795 : 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 218795 : pari_sp av = avma;
1196 : multiple *table;
1197 :
1198 218795 : if (!c6) {
1199 12 : GEN ap = ap_j1728(utoi(c4), utoipos(p));
1200 12 : return gc_long(av, p+1 - itos(ap));
1201 : }
1202 :
1203 218783 : pordmin = (ulong)(1 + 4*sqrt((double)p));
1204 218783 : p1p = p+1;
1205 218783 : p2p = p1p << 1;
1206 218783 : x = 0; KRO = 0;
1207 218783 : switch(Flx_nbroots(mkvecsmall5(0L, c6,c4,0L,1L), p))
1208 : {
1209 44434 : case 3: A = 0; B = 4; break;
1210 106965 : case 1: A = 0; B = 2; break;
1211 67384 : default: A = 1; B = 2; break; /* 0 */
1212 : }
1213 : for(;;)
1214 : { /* see comments in Fp_ellcard_Shanks */
1215 225011 : h = uclosest_lift(A, B, p1p);
1216 225011 : if (!KRO) /* first time, initialize */
1217 : {
1218 218783 : KRO = krouu(c6,p); /* != 0 */
1219 218783 : f = mkvecsmall2(0, Fl_sqr(c6,p));
1220 : }
1221 : else
1222 : {
1223 6228 : KRO = -KRO;
1224 6228 : f = Fl_ellpoint(KRO, &x, c4,c6,p);
1225 : }
1226 225011 : cp4 = Fl_mul(c4, f[2], p);
1227 225011 : fh = Fle_mulu(f, h, cp4, p);
1228 225011 : if (ell_is_inf(fh)) goto FOUND;
1229 :
1230 219219 : s = (ulong) (sqrt(((double)pordmin)/B) / 2);
1231 219219 : if (!s) s = 1;
1232 219219 : table = (multiple *) stack_malloc((s+1) * sizeof(multiple));
1233 219219 : F = Fle_mulu(f, B, cp4, p);
1234 2844224 : for (i=0; i < s; i++)
1235 : {
1236 2634852 : table[i].x = fh[1];
1237 2634852 : table[i].y = fh[2];
1238 2634852 : table[i].i = i;
1239 2634852 : if (Fle_add_inplace(fh, F, cp4, p)) { h += B*(i+1); goto FOUND; }
1240 : }
1241 209372 : qsort(table,s,sizeof(multiple),compare_multiples);
1242 209372 : fg = Fle_mulu(F, s, cp4, p); ftest = zv_copy(fg);
1243 209372 : 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 2496327 : for (i=1; ; i++)
1248 : {
1249 2496327 : l=0; r=s;
1250 17407724 : while (l<r)
1251 : {
1252 14911397 : ulong m = (l+r) >> 1;
1253 14911397 : if (table[m].x < uel(ftest,1)) l=m+1; else r=m;
1254 : }
1255 2496327 : if (r < s && table[r].x == uel(ftest,1)) break;
1256 2286955 : if (Fle_add_inplace(ftest, fg, cp4, p)) pari_err_PRIME("ellap",utoi(p));
1257 : }
1258 209372 : h += table[r].i * B;
1259 209372 : if (table[r].y == uel(ftest,2))
1260 108734 : h -= s * i * B;
1261 : else
1262 100638 : h += s * i * B;
1263 225011 : FOUND:
1264 225011 : h = itou(Fle_order(f, utoipos(h), cp4, p));
1265 : /* h | #E_u(Fp) = A (mod B) */
1266 : {
1267 : GEN C;
1268 225011 : A = itou( Z_chinese_all(gen_0, utoi(A), utoipos(h), utoipos(B), &C) );
1269 225011 : if (abscmpiu(C, pordmin) >= 0) { /* uclosest_lift could overflow */
1270 218783 : h = itou( closest_lift(utoi(A), C, utoipos(p1p)) );
1271 218783 : break;
1272 : }
1273 6228 : B = itou(C);
1274 : }
1275 6228 : A = (p2p - A) % B; set_avma(av);
1276 : }
1277 218783 : 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 85878 : ap_j0(GEN a6,GEN p)
1284 : {
1285 : GEN a, b, e, d;
1286 85878 : if (umodiu(p,3) != 1) return gen_0;
1287 36402 : (void)cornacchia2(utoipos(27),p, &a,&b);
1288 36402 : if (umodiu(a, 3) == 1) a = negi(a);
1289 36402 : d = mulis(a6,-108);
1290 36402 : e = diviuexact(shifti(p,-1), 3); /* (p-1) / 6 */
1291 36402 : return centermod(mulii(a, Fp_pow(d, e, p)), p);
1292 : }
1293 : static GEN
1294 2465256 : ap_j1728(GEN a4,GEN p)
1295 : {
1296 : GEN a, b, e;
1297 2465256 : if (mod4(p) != 1) return gen_0;
1298 1131618 : (void)cornacchia2(utoipos(4),p, &a,&b);
1299 1131618 : if (Mod4(a)==0) a = b;
1300 1131618 : if (Mod2(a)==1) a = shifti(a,1);
1301 1131618 : if (Mod8(a)==6) a = negi(a);
1302 1131618 : e = shifti(p,-2); /* (p-1) / 4 */
1303 1131618 : return centermod(mulii(a, Fp_pow(a4, e, p)), p);
1304 : }
1305 : static GEN
1306 114 : ap_j8000(GEN a6, GEN p)
1307 : {
1308 : GEN a, b;
1309 114 : long r = mod8(p), s = 1;
1310 114 : if (r != 1 && r != 3) return gen_0;
1311 42 : (void)cornacchia2(utoipos(8),p, &a,&b);
1312 42 : switch(Mod16(a)) {
1313 12 : case 2: case 6: if (Mod4(b)) s = -s;
1314 12 : break;
1315 30 : case 10: case 14: if (!Mod4(b)) s = -s;
1316 30 : break;
1317 : }
1318 42 : if (kronecker(mulis(a6, 42), p) < 0) s = -s;
1319 42 : return s > 0? a: negi(a);
1320 : }
1321 : static GEN
1322 120 : ap_j287496(GEN a6, GEN p)
1323 : {
1324 : GEN a, b;
1325 120 : long s = 1;
1326 120 : if (mod4(p) != 1) return gen_0;
1327 60 : (void)cornacchia2(utoipos(4),p, &a,&b);
1328 60 : if (Mod4(a)==0) a = b;
1329 60 : if (Mod2(a)==1) a = shifti(a,1);
1330 60 : if (Mod8(a)==6) s = -s;
1331 60 : if (krosi(2,p) < 0) s = -s;
1332 60 : if (kronecker(mulis(a6, -14), p) < 0) s = -s;
1333 60 : return s > 0? a: negi(a);
1334 : }
1335 : static GEN
1336 1200 : ap_cm(int CM, long A6B, GEN a6, GEN p)
1337 : {
1338 : GEN a, b;
1339 1200 : long s = 1;
1340 1200 : if (krosi(CM,p) < 0) return gen_0;
1341 552 : (void)cornacchia2(utoipos(-CM),p, &a, &b);
1342 552 : if ((CM&3) == 0) CM >>= 2;
1343 552 : if ((krois(a, -CM) > 0) ^ (CM == -7)) s = -s;
1344 552 : if (kronecker(mulis(a6,A6B), p) < 0) s = -s;
1345 552 : return s > 0? a: negi(a);
1346 : }
1347 : static GEN
1348 426468 : ec_ap_cm(int CM, GEN a4, GEN a6, GEN p)
1349 : {
1350 426468 : switch(CM)
1351 : {
1352 24954 : case -3: return ap_j0(a6, p);
1353 400080 : case -4: return ap_j1728(a4, p);
1354 114 : case -8: return ap_j8000(a6, p);
1355 120 : case -16: return ap_j287496(a6, p);
1356 138 : case -7: return ap_cm(CM, -2, a6, p);
1357 138 : case -11: return ap_cm(CM, 21, a6, p);
1358 150 : case -12: return ap_cm(CM, 22, a6, p);
1359 126 : case -19: return ap_cm(CM, 1, a6, p);
1360 132 : case -27: return ap_cm(CM, 253, a6, p);
1361 126 : case -28: return ap_cm(-7, -114, a6, p); /* yes, -7 ! */
1362 138 : case -43: return ap_cm(CM, 21, a6, p);
1363 126 : case -67: return ap_cm(CM, 217, a6, p);
1364 126 : case -163:return ap_cm(CM, 185801, a6, p);
1365 0 : default: return NULL;
1366 : }
1367 : }
1368 :
1369 : static GEN
1370 42636 : Fp_ellj_nodiv(GEN a4, GEN a6, GEN p)
1371 : {
1372 42636 : GEN a43 = Fp_mulu(Fp_powu(a4, 3, p), 4, p);
1373 42636 : GEN a62 = Fp_mulu(Fp_sqr(a6, p), 27, p);
1374 42636 : return mkvec2(Fp_mulu(a43, 1728, p), Fp_add(a43, a62, p));
1375 : }
1376 :
1377 : GEN
1378 84 : Fp_ellj(GEN a4, GEN a6, GEN p)
1379 : {
1380 84 : pari_sp av = avma;
1381 : GEN z;
1382 84 : 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 84 : z = Fp_ellj_nodiv(a4, a6, p);
1388 84 : return gc_INT(av,Fp_div(gel(z,1),gel(z,2),p));
1389 : }
1390 :
1391 : void
1392 936 : Fp_ellj_to_a4a6(GEN j, GEN p, GEN *pt_a4, GEN *pt_a6)
1393 : {
1394 936 : j = modii(j, p);
1395 936 : if (signe(j) == 0) { *pt_a4 = gen_0; *pt_a6 = gen_1; }
1396 642 : else if (equaliu(j,umodui(1728,p))) { *pt_a4 = gen_1; *pt_a6 = gen_0; }
1397 : else
1398 : {
1399 504 : GEN k = Fp_sub(utoi(1728), j, p);
1400 504 : GEN kj = Fp_mul(k, j, p);
1401 504 : GEN k2j = Fp_mul(kj, k, p);
1402 504 : *pt_a4 = Fp_mulu(kj, 3, p);
1403 504 : *pt_a6 = Fp_double(k2j, p);
1404 : }
1405 936 : }
1406 :
1407 : static GEN /* Only compute a mod p, so assume p>=17 */
1408 2168640 : Fp_ellcard_CM(GEN a4, GEN a6, GEN p)
1409 : {
1410 2168640 : pari_sp av = avma;
1411 : GEN a;
1412 2168640 : if (!signe(a4)) a = ap_j0(a6,p);
1413 2107716 : else if (!signe(a6)) a = ap_j1728(a4,p);
1414 : else
1415 : {
1416 42552 : GEN j = Fp_ellj_nodiv(a4, a6, p);
1417 42552 : long CM = Fp_ellj_get_CM(gel(j,1), gel(j,2), p);
1418 42552 : if (!CM) return gc_NULL(av);
1419 1434 : a = ec_ap_cm(CM,a4,a6,p);
1420 : }
1421 2127522 : return gc_INT(av, subii(addiu(p,1),a));
1422 : }
1423 :
1424 : GEN
1425 2429244 : Fp_ellcard(GEN a4, GEN a6, GEN p)
1426 : {
1427 2429244 : long lp = expi(p);
1428 2429244 : ulong pp = p[2];
1429 2429244 : if (lp < 11)
1430 260604 : return utoi(pp+1 - Fl_elltrace_naive(umodiu(a4,pp), umodiu(a6,pp), pp));
1431 2168640 : { GEN a = Fp_ellcard_CM(a4,a6,p); if (a) return a; }
1432 41118 : if (lp >= 56)
1433 744 : return Fp_ellcard_SEA(a4, a6, p, 0);
1434 40374 : if (lp <= BITS_IN_LONG-2)
1435 40295 : 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 533326 : Fl_elltrace(ulong a4, ulong a6, ulong p)
1441 : {
1442 : pari_sp av;
1443 : long lp;
1444 : GEN a;
1445 533326 : if (p < (1<<11)) return Fl_elltrace_naive(a4, a6, p);
1446 178500 : lp = expu(p);
1447 178500 : 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 999048 : Fl_elltrace_CM(long CM, ulong a4, ulong a6, ulong p)
1453 : {
1454 : pari_sp av;
1455 : GEN a;
1456 999048 : if (!CM) return Fl_elltrace(a4,a6,p);
1457 466380 : if (p < (1<<11)) return Fl_elltrace_naive(a4, a6, p);
1458 425034 : av = avma; a = ec_ap_cm(CM, utoi(a4), utoi(a6), utoipos(p));
1459 425034 : return gc_long(av, itos(a));
1460 : }
1461 :
1462 : static GEN
1463 62334 : _FpE_pairorder(void *E, GEN P, GEN Q, GEN m, GEN F)
1464 : {
1465 62334 : struct _FpE *e = (struct _FpE *) E;
1466 62334 : return Fp_order(FpE_weilpairing(P,Q,m,e->a4,e->p), F, e->p);
1467 : }
1468 :
1469 : GEN
1470 103470 : Fp_ellgroup(GEN a4, GEN a6, GEN N, GEN p, GEN *pt_m)
1471 : {
1472 : struct _FpE e;
1473 103470 : e.a4=a4; e.a6=a6; e.p=p;
1474 103470 : return gen_ellgroup(N, subiu(p,1), pt_m, (void*)&e, &FpE_group, _FpE_pairorder);
1475 : }
1476 :
1477 : GEN
1478 492 : Fp_ellgens(GEN a4, GEN a6, GEN ch, GEN D, GEN m, GEN p)
1479 : {
1480 : GEN P;
1481 492 : pari_sp av = avma;
1482 : struct _FpE e;
1483 492 : e.a4=a4; e.a6=a6; e.p=p;
1484 492 : switch(lg(D)-1)
1485 : {
1486 408 : case 1:
1487 408 : P = gen_gener(gel(D,1), (void*)&e, &FpE_group);
1488 408 : P = mkvec(FpE_changepoint(P, ch, p));
1489 408 : break;
1490 84 : default:
1491 84 : P = gen_ellgens(gel(D,1), gel(D,2), m, (void*)&e, &FpE_group, _FpE_pairorder);
1492 84 : gel(P,1) = FpE_changepoint(gel(P,1), ch, p);
1493 84 : gel(P,2) = FpE_changepoint(gel(P,2), ch, p);
1494 84 : break;
1495 : }
1496 492 : 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 838 : RgE_to_FpXQE(GEN x, GEN T, GEN p)
1513 : {
1514 838 : if (ell_is_inf(x)) return x;
1515 838 : retmkvec2(Rg_to_FpXQ(gel(x,1),T,p),Rg_to_FpXQ(gel(x,2),T,p));
1516 : }
1517 :
1518 : GEN
1519 1610 : FpXQE_changepoint(GEN x, GEN ch, GEN T, GEN p)
1520 : {
1521 1610 : pari_sp av = avma;
1522 : GEN p1,z,u,r,s,t,v,v2,v3;
1523 1610 : if (ell_is_inf(x)) return x;
1524 808 : u = gel(ch,1); r = gel(ch,2);
1525 808 : s = gel(ch,3); t = gel(ch,4);
1526 808 : v = FpXQ_inv(u, T, p); v2 = FpXQ_sqr(v, T, p); v3 = FpXQ_mul(v,v2, T, p);
1527 808 : p1 = FpX_sub(gel(x,1),r, p);
1528 808 : z = cgetg(3,t_VEC);
1529 808 : gel(z,1) = FpXQ_mul(v2, p1, T, p);
1530 808 : 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 808 : return gc_upto(av, z);
1532 : }
1533 :
1534 : GEN
1535 838 : FpXQE_changepointinv(GEN x, GEN ch, GEN T, GEN p)
1536 : {
1537 : GEN u, r, s, t, X, Y, u2, u3, u2X, z;
1538 838 : if (ell_is_inf(x)) return x;
1539 838 : X = gel(x,1); Y = gel(x,2);
1540 838 : u = gel(ch,1); r = gel(ch,2);
1541 838 : s = gel(ch,3); t = gel(ch,4);
1542 838 : u2 = FpXQ_sqr(u, T, p); u3 = FpXQ_mul(u,u2, T, p);
1543 838 : u2X = FpXQ_mul(u2,X, T, p);
1544 838 : z = cgetg(3, t_VEC);
1545 838 : gel(z,1) = FpX_add(u2X,r, p);
1546 838 : gel(z,2) = FpX_add(FpXQ_mul(u3,Y, T, p), FpX_add(FpXQ_mul(s,u2X, T, p), t, p), p);
1547 838 : return z;
1548 : }
1549 :
1550 : static GEN
1551 720 : random_nonsquare_FpXQ(GEN T, GEN p)
1552 : {
1553 720 : pari_sp av = avma;
1554 720 : long n = degpol(T), v = varn(T);
1555 : GEN a;
1556 720 : if (odd(n))
1557 : {
1558 360 : GEN z = cgetg(3, t_POL);
1559 360 : z[1] = evalsigne(1) | evalvarn(v);
1560 360 : gel(z,2) = random_nonsquare_Fp(p); return z;
1561 : }
1562 : do
1563 : {
1564 678 : set_avma(av);
1565 678 : a = random_FpX(n, v, p);
1566 678 : } while (FpXQ_issquare(a, T, p));
1567 360 : return a;
1568 : }
1569 :
1570 : void
1571 720 : FpXQ_elltwist(GEN a4, GEN a6, GEN T, GEN p, GEN *pt_a4, GEN *pt_a6)
1572 : {
1573 720 : GEN d = random_nonsquare_FpXQ(T, p);
1574 720 : GEN d2 = FpXQ_sqr(d, T, p), d3 = FpXQ_mul(d2, d, T, p);
1575 720 : *pt_a4 = FpXQ_mul(a4, d2, T, p);
1576 720 : *pt_a6 = FpXQ_mul(a6, d3, T, p);
1577 720 : }
1578 :
1579 : static GEN
1580 230033 : FpXQE_dbl_slope(GEN P, GEN a4, GEN T, GEN p, GEN *slope)
1581 : {
1582 : GEN x, y, Q;
1583 230033 : if (ell_is_inf(P) || !signe(gel(P,2))) return ellinf();
1584 228630 : x = gel(P,1); y = gel(P,2);
1585 228630 : *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 228630 : Q = cgetg(3,t_VEC);
1588 228630 : gel(Q, 1) = FpX_sub(FpXQ_sqr(*slope, T, p), FpX_mulu(x, 2, p), p);
1589 228630 : gel(Q, 2) = FpX_sub(FpXQ_mul(*slope, FpX_sub(x, gel(Q, 1), p), T, p), y, p);
1590 228630 : return Q;
1591 : }
1592 :
1593 : GEN
1594 226517 : FpXQE_dbl(GEN P, GEN a4, GEN T, GEN p)
1595 : {
1596 226517 : pari_sp av = avma;
1597 : GEN slope;
1598 226517 : return gc_upto(av, FpXQE_dbl_slope(P,a4,T,p,&slope));
1599 : }
1600 :
1601 : static GEN
1602 215680 : FpXQE_add_slope(GEN P, GEN Q, GEN a4, GEN T, GEN p, GEN *slope)
1603 : {
1604 : GEN Px, Py, Qx, Qy, R;
1605 215680 : if (ell_is_inf(P)) return Q;
1606 215668 : if (ell_is_inf(Q)) return P;
1607 215668 : Px = gel(P,1); Py = gel(P,2);
1608 215668 : Qx = gel(Q,1); Qy = gel(Q,2);
1609 215668 : if (ZX_equal(Px, Qx))
1610 : {
1611 518 : if (ZX_equal(Py, Qy))
1612 6 : return FpXQE_dbl_slope(P, a4, T, p, slope);
1613 : else
1614 512 : return ellinf();
1615 : }
1616 215150 : *slope = FpXQ_div(FpX_sub(Py, Qy, p), FpX_sub(Px, Qx, p), T, p);
1617 215150 : R = cgetg(3,t_VEC);
1618 215150 : gel(R, 1) = FpX_sub(FpX_sub(FpXQ_sqr(*slope, T, p), Px, p), Qx, p);
1619 215150 : gel(R, 2) = FpX_sub(FpXQ_mul(*slope, FpX_sub(Px, gel(R, 1), p), T, p), Py, p);
1620 215150 : return R;
1621 : }
1622 :
1623 : GEN
1624 215200 : FpXQE_add(GEN P, GEN Q, GEN a4, GEN T, GEN p)
1625 : {
1626 215200 : pari_sp av = avma;
1627 : GEN slope;
1628 215200 : 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 62855 : FpXQE_neg(GEN P, GEN T, GEN p)
1640 : {
1641 : (void) T;
1642 62855 : if (ell_is_inf(P)) return ellinf();
1643 62855 : 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 226517 : _FpXQE_dbl(void *E, GEN P)
1657 : {
1658 226517 : struct _FpXQE *ell = (struct _FpXQE *) E;
1659 226517 : return FpXQE_dbl(P, ell->a4, ell->T, ell->p);
1660 : }
1661 : static GEN
1662 215200 : _FpXQE_add(void *E, GEN P, GEN Q)
1663 : {
1664 215200 : struct _FpXQE *ell=(struct _FpXQE *) E;
1665 215200 : return FpXQE_add(P, Q, ell->a4, ell->T, ell->p);
1666 : }
1667 : static GEN
1668 70102 : _FpXQE_mul(void *E, GEN P, GEN n)
1669 : {
1670 70102 : pari_sp av = avma;
1671 70102 : struct _FpXQE *e=(struct _FpXQE *) E;
1672 70102 : long s = signe(n);
1673 70102 : if (!s || ell_is_inf(P)) return ellinf();
1674 70102 : if (s<0) P = FpXQE_neg(P, e->T, e->p);
1675 70102 : if (is_pm1(n)) return s>0? gcopy(P): P;
1676 7168 : return gc_GEN(av, gen_pow_i(P, n, e, &_FpXQE_dbl, &_FpXQE_add));
1677 : }
1678 :
1679 : GEN
1680 802 : FpXQE_mul(GEN P, GEN n, GEN a4, GEN T, GEN p)
1681 : {
1682 : struct _FpXQE E;
1683 802 : E.a4= a4; E.T = T; E.p = p;
1684 802 : return _FpXQE_mul(&E, P, n);
1685 : }
1686 :
1687 : /* Finds a random nonsingular point on E */
1688 :
1689 : GEN
1690 930 : random_FpXQE(GEN a4, GEN a6, GEN T, GEN p)
1691 : {
1692 930 : pari_sp ltop = avma;
1693 : GEN x, x2, y, rhs;
1694 930 : long v = get_FpX_var(T), d = get_FpX_degree(T);
1695 : do
1696 : {
1697 1890 : set_avma(ltop);
1698 1890 : x = random_FpX(d,v,p); /* x^3+a4*x+a6 = x*(x^2+a4)+a6 */
1699 1890 : x2 = FpXQ_sqr(x, T, p);
1700 1890 : 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 1890 : || !FpXQ_issquare(rhs, T, p));
1703 930 : y = FpXQ_sqrt(rhs, T, p);
1704 930 : if (!y) pari_err_PRIME("random_FpE", p);
1705 930 : return gc_GEN(ltop, mkvec2(x, y));
1706 : }
1707 :
1708 : static GEN
1709 129 : _FpXQE_rand(void *E)
1710 : {
1711 129 : struct _FpXQE *e=(struct _FpXQE *) E;
1712 129 : 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 14 : get_FpXQE_group(void ** pt_E, GEN a4, GEN a6, GEN T, GEN p)
1719 : {
1720 14 : struct _FpXQE *e = (struct _FpXQE *) stack_malloc(sizeof(struct _FpXQE));
1721 14 : e->a4 = a4; e->a6 = a6; e->T = T; e->p = p;
1722 14 : *pt_E = (void *) e;
1723 14 : return &FpXQE_group;
1724 : }
1725 :
1726 : GEN
1727 12 : FpXQE_order(GEN z, GEN o, GEN a4, GEN T, GEN p)
1728 : {
1729 12 : pari_sp av = avma;
1730 : struct _FpXQE e;
1731 12 : e.a4=a4; e.T=T; e.p=p;
1732 12 : 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 4104 : FpXQE_vert(GEN P, GEN Q, GEN a4, GEN T, GEN p)
1754 : {
1755 4104 : long vT = get_FpX_var(T);
1756 4104 : if (ell_is_inf(P))
1757 60 : return pol_1(get_FpX_var(T));
1758 4044 : if (!ZX_equal(gel(Q, 1), gel(P, 1)))
1759 4044 : 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 3990 : FpXQE_Miller_line(GEN R, GEN Q, GEN slope, GEN a4, GEN T, GEN p)
1767 : {
1768 3990 : long vT = get_FpX_var(T);
1769 3990 : GEN x = gel(Q, 1), y = gel(Q, 2);
1770 3990 : GEN tmp1 = FpX_sub(x, gel(R, 1), p);
1771 3990 : GEN tmp2 = FpX_add(FpXQ_mul(tmp1, slope, T, p), gel(R, 2), p);
1772 3990 : if (!ZX_equal(y, tmp2))
1773 3990 : 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 3564 : FpXQE_tangent_update(GEN R, GEN Q, GEN a4, GEN T, GEN p, GEN *pt_R)
1794 : {
1795 3564 : if (ell_is_inf(R))
1796 : {
1797 6 : *pt_R = ellinf();
1798 6 : return pol_1(get_FpX_var(T));
1799 : }
1800 3558 : else if (!signe(gel(R,2)))
1801 : {
1802 48 : *pt_R = ellinf();
1803 48 : return FpXQE_vert(R, Q, a4, T, p);
1804 : } else {
1805 : GEN slope;
1806 3510 : *pt_R = FpXQE_dbl_slope(R, a4, T, p, &slope);
1807 3510 : 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 486 : FpXQE_chord_update(GEN R, GEN P, GEN Q, GEN a4, GEN T, GEN p, GEN *pt_R)
1817 : {
1818 486 : if (ell_is_inf(R))
1819 : {
1820 0 : *pt_R = gcopy(P);
1821 0 : return FpXQE_vert(P, Q, a4, T, p);
1822 : }
1823 486 : 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 486 : else if (ZX_equal(gel(P, 1), gel(R, 1)))
1829 : {
1830 6 : 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 6 : *pt_R = ellinf();
1835 6 : return FpXQE_vert(R, Q, a4, T, p);
1836 : }
1837 : } else {
1838 : GEN slope;
1839 480 : *pt_R = FpXQE_add_slope(P, R, a4, T, p, &slope);
1840 480 : 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 3564 : FpXQE_Miller_dbl(void* E, GEN d)
1847 : {
1848 3564 : struct _FpXQE_miller *m = (struct _FpXQE_miller *)E;
1849 3564 : GEN p = m->p;
1850 3564 : GEN T = m->T, a4 = m->a4, P = m->P;
1851 : GEN v, line;
1852 3564 : GEN N = FpXQ_sqr(gel(d,1), T, p);
1853 3564 : GEN D = FpXQ_sqr(gel(d,2), T, p);
1854 3564 : GEN point = gel(d,3);
1855 3564 : line = FpXQE_tangent_update(point, P, a4, T, p, &point);
1856 3564 : N = FpXQ_mul(N, line, T, p);
1857 3564 : v = FpXQE_vert(point, P, a4, T, p);
1858 3564 : D = FpXQ_mul(D, v, T, p); return mkvec3(N, D, point);
1859 : }
1860 :
1861 : static GEN
1862 486 : FpXQE_Miller_add(void* E, GEN va, GEN vb)
1863 : {
1864 486 : struct _FpXQE_miller *m = (struct _FpXQE_miller *)E;
1865 486 : GEN p = m->p;
1866 486 : GEN T = m->T, a4 = m->a4, P = m->P;
1867 : GEN v, line, point;
1868 486 : GEN na = gel(va,1), da = gel(va,2), pa = gel(va,3);
1869 486 : GEN nb = gel(vb,1), db = gel(vb,2), pb = gel(vb,3);
1870 486 : GEN N = FpXQ_mul(na, nb, T, p);
1871 486 : GEN D = FpXQ_mul(da, db, T, p);
1872 486 : line = FpXQE_chord_update(pa, pb, P, a4, T, p, &point);
1873 486 : N = FpXQ_mul(N, line, T, p);
1874 486 : v = FpXQE_vert(point, P, a4, T, p);
1875 486 : 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 54 : FpXQE_Miller(GEN Q, GEN P, GEN m, GEN a4, GEN T, GEN p)
1882 : {
1883 54 : pari_sp av = avma;
1884 : struct _FpXQE_miller d;
1885 : GEN v, N, D, g1;
1886 :
1887 54 : d.a4 = a4; d.T = T; d.p = p; d.P = P;
1888 54 : g1 = pol_1(get_FpX_var(T));
1889 54 : v = gen_pow_i(mkvec3(g1,g1,Q), m, (void*)&d,
1890 : FpXQE_Miller_dbl, FpXQE_Miller_add);
1891 54 : N = gel(v,1); D = gel(v,2);
1892 54 : return gc_upto(av, FpXQ_div(N, D, T, p));
1893 : }
1894 :
1895 : GEN
1896 24 : FpXQE_weilpairing(GEN P, GEN Q, GEN m, GEN a4, GEN T, GEN p)
1897 : {
1898 24 : pari_sp av = avma;
1899 : GEN N, D, w;
1900 24 : if (ell_is_inf(P) || ell_is_inf(Q) || ZXV_equal(P,Q))
1901 0 : return pol_1(get_FpX_var(T));
1902 24 : N = FpXQE_Miller(P, Q, m, a4, T, p);
1903 24 : D = FpXQE_Miller(Q, P, m, a4, T, p);
1904 24 : w = FpXQ_div(N, D, T, p);
1905 24 : if (mpodd(m)) w = FpX_neg(w, p);
1906 24 : return gc_upto(av, w);
1907 : }
1908 :
1909 : GEN
1910 6 : FpXQE_tatepairing(GEN P, GEN Q, GEN m, GEN a4, GEN T, GEN p)
1911 : {
1912 6 : if (ell_is_inf(P) || ell_is_inf(Q)) return pol_1(get_FpX_var(T));
1913 6 : return FpXQE_Miller(P, Q, m, a4, T, p);
1914 : }
1915 :
1916 : /***********************************************************************/
1917 : /** **/
1918 : /** issupersingular **/
1919 : /** **/
1920 : /***********************************************************************/
1921 :
1922 : GEN
1923 1473 : FpXQ_ellj(GEN a4, GEN a6, GEN T, GEN p)
1924 : {
1925 1473 : if (absequaliu(p,3)) return pol_0(get_FpX_var(T));
1926 : else
1927 : {
1928 1473 : pari_sp av=avma;
1929 1473 : GEN a43 = FpXQ_mul(a4,FpXQ_sqr(a4,T,p),T,p);
1930 1473 : GEN a62 = FpXQ_sqr(a6,T,p);
1931 1473 : GEN num = FpX_mulu(a43,6912,p);
1932 1473 : GEN den = FpX_add(FpX_mulu(a43,4,p),FpX_mulu(a62,27,p),p);
1933 1473 : return gc_leaf(av, FpXQ_div(num, den, T, p));
1934 : }
1935 : }
1936 :
1937 : static GEN
1938 28740 : FpXQ_is_quad(GEN x, GEN T, GEN p)
1939 : {
1940 28740 : pari_sp av = avma;
1941 : GEN K;
1942 28740 : long d = degpol(T);
1943 28740 : x = FpXQ_red(x,T,p);
1944 28740 : if (lgpol(x)<=1) return NULL;
1945 28740 : if (d==2) return FpXQ_minpoly(x, T, p);
1946 28740 : if (odd(degpol(T))) return NULL;
1947 28740 : K = FpM_ker(FpXQ_matrix_pow(x, d, 3, T, p), p);
1948 28740 : if (lg(K)!=2) return gc_NULL(av);
1949 504 : return RgV_to_RgX(gel(K,1), get_FpX_var(T));
1950 : }
1951 :
1952 : int
1953 141870 : FpXQ_elljissupersingular(GEN j, GEN T, GEN p)
1954 : {
1955 141870 : 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 141870 : ulong d = get_FpX_degree(T);
1962 : GEN S;
1963 141870 : if (degpol(j) <= 0) return Fp_elljissupersingular(constant_coeff(j), p);
1964 141090 : j = FpXQ_red(j, T, p);
1965 141090 : if (degpol(j) <= 0) return gc_bool(ltop, Fp_elljissupersingular(constant_coeff(j), p));
1966 : /* Now j is not in F_p */
1967 141090 : if (abscmpiu(p, 5) <= 0) return gc_bool(ltop,0); /* j != 0*/
1968 141084 : if (odd(d)) return 0;
1969 : /* Set S so that FF_p[T]/(S) is isomorphic to FF_{p^2}: */
1970 40236 : if (d == 2)
1971 11496 : S = T;
1972 : else /* d > 2 */
1973 : {
1974 28740 : S = FpXQ_is_quad(j, T, p);
1975 28740 : if (!S) return gc_bool(ltop,0);
1976 504 : j = pol_x(varn(S));
1977 : }
1978 12000 : return gc_bool(ltop, jissupersingular(j,S,p));
1979 : }
1980 :
1981 : int
1982 900 : Fq_elljissupersingular(GEN j, GEN T, GEN p)
1983 822 : { return typ(j)==t_INT? Fp_elljissupersingular(j, p)
1984 1722 : : FpXQ_elljissupersingular(j, T, p); }
1985 :
1986 : /* p > 5 prime; return d such that (-d/p) = -1 */
1987 : static ulong
1988 1014 : find_inert_disc(GEN p)
1989 : {
1990 1014 : long s = mod4(p) == 1? -1: 1; /* - (-1/p) */
1991 1014 : ulong d = 3;
1992 : while(1)
1993 : {
1994 1020 : if (kroui(d,p) == s) return d; /* = 3 mod (16) */
1995 510 : d++;
1996 510 : if (kroui(d>>2,p) == s) return d; /* = 4 mod (16) */
1997 228 : d += 3;
1998 228 : if (kroui(d,p) == s) return d; /* = 7 mod (16) */
1999 90 : d++;
2000 90 : if (kroui(d>>2,p) == s) return d; /* = 8 mod (16) */
2001 30 : d += 3;
2002 30 : if (kroui(d,p) == s) return d; /* = 11 mod (16) */
2003 6 : d += 4;
2004 6 : if (kroui(d,p) == s) return d; /* = 15 mod (16) */
2005 6 : d += 4;
2006 : }
2007 : }
2008 :
2009 : /* p > 5 */
2010 : static GEN
2011 1014 : ellsupersingularj_easy_FpXQ(GEN T, GEN p)
2012 : {
2013 1014 : long d = find_inert_disc(p);
2014 1014 : GEN R = FpXQX_roots(polclass(stoi(-d), 0, 0), T, p);
2015 1014 : return gel(R,1);
2016 : }
2017 :
2018 : GEN
2019 1032 : ellsupersingularj_FpXQ(GEN T, GEN p)
2020 : {
2021 : GEN j, j2, R, Phi2;
2022 : long i, ep, lp;
2023 1032 : if (cmpiu(p, 5) <= 0) return pol_0(get_FpX_var(T));
2024 1014 : j2 = ellsupersingularj_easy_FpXQ(T, p);
2025 1014 : Phi2 = polmodular_ZXX(2,0,0,1);
2026 1014 : R = FpXQX_roots(FqXY_evalx(Phi2, j2, T, p), T, p);
2027 1014 : j = gel(R,1+random_Fl(lg(R)-1));
2028 1014 : ep = expi(p); lp = ep + random_Fl(ep);
2029 15372 : for (i = 1; i <= lp; i++)
2030 : {
2031 14358 : GEN Phi2_j = FqX_div_by_X_x(FqXY_evalx(Phi2, j, T, p), j2, T, p, NULL);
2032 14358 : R = FqX_quad_root(Phi2_j, T, p);
2033 14358 : if (!R) pari_err_PRIME("ellsupersingularj",p);
2034 14358 : j2 = j; j = random_bits(1) ? R: Fq_neg(Fq_add(gel(Phi2_j,3), R, T, p), T, p);
2035 : }
2036 1014 : return j;
2037 : }
2038 :
2039 : /***********************************************************************/
2040 : /** **/
2041 : /** Point counting **/
2042 : /** **/
2043 : /***********************************************************************/
2044 :
2045 : GEN
2046 13260 : elltrace_extension(GEN t, long n, GEN q)
2047 : {
2048 13260 : pari_sp av = avma;
2049 13260 : GEN v = RgX_to_RgC(RgXQ_powu(pol_x(0), n, mkpoln(3,gen_1,negi(t),q)),2);
2050 13260 : GEN te = addii(shifti(gel(v,1),1), mulii(t,gel(v,2)));
2051 13260 : return gc_INT(av, te);
2052 : }
2053 :
2054 : GEN
2055 12606 : Fp_ffellcard(GEN a4, GEN a6, GEN q, long n, GEN p)
2056 : {
2057 12606 : pari_sp av = avma;
2058 12606 : GEN ap = subii(addiu(p, 1), Fp_ellcard(a4, a6, p));
2059 12606 : GEN te = elltrace_extension(ap, n, p);
2060 12606 : return gc_INT(av, subii(addiu(q, 1), te));
2061 : }
2062 :
2063 : static GEN
2064 1446 : FpXQ_ellcardj(GEN a4, GEN a6, GEN j, GEN T, GEN q, GEN p, long n)
2065 : {
2066 1446 : GEN q1 = addiu(q,1);
2067 1446 : if (signe(j)==0)
2068 : {
2069 : GEN W, w, t, N;
2070 480 : if (umodiu(q,6)!=1) return q1;
2071 360 : N = Fp_ffellcard(gen_0,gen_1,q,n,p);
2072 360 : t = subii(q1, N);
2073 360 : W = FpXQ_pow(a6,diviuexact(shifti(q,-1), 3),T,p);
2074 360 : if (degpol(W)>0) /*p=5 mod 6*/
2075 90 : return ZX_equal1(FpXQ_powu(W,3,T,p)) ? addii(q1,shifti(t,-1)):
2076 30 : subii(q1,shifti(t,-1));
2077 300 : w = modii(gel(W,2),p);
2078 300 : if (equali1(w)) return N;
2079 222 : if (equalii(w,subiu(p,1))) return addii(q1,t);
2080 : else /*p=1 mod 6*/
2081 : {
2082 144 : GEN u = shifti(t,-1), v = sqrtint(diviuexact(subii(q,sqri(u)),3));
2083 144 : GEN a = addii(u,v), b = shifti(v,1);
2084 144 : if (equali1(Fp_powu(w,3,p)))
2085 : {
2086 72 : if (dvdii(addmulii(a, w, b), p))
2087 48 : return subii(q1,subii(shifti(b,1),a));
2088 : else
2089 24 : return addii(q1,addii(a,b));
2090 : }
2091 : else
2092 : {
2093 72 : if (dvdii(submulii(a, w, b), p))
2094 48 : return subii(q1,subii(a,shifti(b,1)));
2095 : else
2096 24 : return subii(q1,addii(a,b));
2097 : }
2098 : }
2099 966 : } else if (equalii(j,modsi(1728,p)))
2100 : {
2101 : GEN w, W, N, t;
2102 486 : if (mod4(q)==3) return q1;
2103 366 : W = FpXQ_pow(a4,shifti(q,-2),T,p);
2104 366 : if (degpol(W)>0) return q1; /*p=3 mod 4*/
2105 306 : w = modii(gel(W,2),p);
2106 306 : N = Fp_ffellcard(gen_1,gen_0,q,n,p);
2107 306 : if (equali1(w)) return N;
2108 204 : t = subii(q1, N);
2109 204 : if (equalii(w,subiu(p,1))) return addii(q1,t);
2110 : else /*p=1 mod 4*/
2111 : {
2112 96 : GEN u = shifti(t,-1), v = sqrtint(subii(q,sqri(u)));
2113 96 : if (dvdii(addmulii(u, w, v), p))
2114 48 : return subii(q1,shifti(v,1));
2115 : else
2116 48 : return addii(q1,shifti(v,1));
2117 : }
2118 : } else
2119 : {
2120 480 : GEN g = Fp_div(j, Fp_sub(utoi(1728), j, p), p);
2121 480 : GEN l = FpXQ_div(FpX_mulu(a6,3,p),FpX_mulu(a4,2,p),T,p);
2122 480 : GEN N = Fp_ffellcard(Fp_mulu(g,3,p),Fp_double(g,p),q,n,p);
2123 480 : if (FpXQ_issquare(l,T,p)) return N;
2124 240 : return subii(shifti(q1,1),N);
2125 : }
2126 : }
2127 :
2128 : static GEN
2129 7 : FpXQ_ffellcard(GEN a4, GEN a6, GEN M, GEN q, GEN T, GEN p, long n)
2130 : {
2131 7 : long m = degpol(M);
2132 7 : GEN j = pol_x(get_FpX_var(T));
2133 7 : GEN g = FpXQ_div(j, Fp_FpX_sub(utoi(1728), j, p), M, p);
2134 7 : GEN N = FpXQ_ellcard(FpX_mulu(g,3,p),FpX_mulu(g,2,p),M,p);
2135 7 : GEN qm = powiu(p, m), q1 = addiu(q, 1), qm1 = addiu(qm, 1);
2136 7 : GEN l = FpXQ_mul(FpX_mulu(a6,3,p),FpX_mulu(a4,2,p),T,p);
2137 7 : GEN te = elltrace_extension(subii(qm1, N), n/m, qm);
2138 7 : return FpXQ_issquare(l,T,p) ? subii(q1, te): addii(q1, te);
2139 : }
2140 :
2141 : static int
2142 6 : FpXQ_is4power(GEN x, GEN T, GEN p)
2143 : {
2144 6 : long d = get_FpX_degree(T);
2145 6 : if (lg(x) == 2 || absequalui(2, p)) return 1;
2146 6 : if (Mod4(p)==1)
2147 6 : 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 6 : FpXQ_ellcard_supersingular(GEN a4, GEN a6, GEN T, GEN p)
2157 : {
2158 6 : pari_sp av = avma;
2159 6 : long d = get_FpX_degree(T);
2160 : GEN r;
2161 6 : 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 6 : else if (signe(a4)==0)
2164 0 : r = FpXQ_ellcardj(a4, a6, gen_0, T, powiu(p, d), p, d);
2165 6 : 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 6 : long qm4 = (odd(d>>1) && Mod4(p)==3);
2171 6 : if (odd(d)) return gen_0;
2172 6 : q2 = powiu(p, d>>1); q = sqri(q2);
2173 6 : t = shifti(q2, 1);
2174 6 : 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 12 : r = qm4 ^ FpXQ_is4power(D, T, p) ? subii(addiu(q, 1), t)
2177 6 : : addii(addiu(q, 1), t);
2178 : }
2179 6 : return gc_INT(av, r);
2180 : }
2181 :
2182 : GEN
2183 18 : Fq_ellcard_supersingular(GEN a4, GEN a6, GEN T, GEN p)
2184 18 : { return T ? FpXQ_ellcard_supersingular(a4, a6, T, p) : addiu(p, 1); }
2185 :
2186 : static GEN
2187 7293 : FpXQ_ellcard_i(GEN a4, GEN a6, GEN T, GEN p)
2188 : {
2189 7293 : long n = get_FpX_degree(T);
2190 7293 : GEN q = powiu(p, n);
2191 7293 : if (degpol(a4)<=0 && degpol(a6)<=0)
2192 654 : return Fp_ffellcard(constant_coeff(a4),constant_coeff(a6),q,n,p);
2193 6639 : if (lgefint(p)==3)
2194 : {
2195 5166 : ulong pp = p[2];
2196 5166 : return Flxq_ellcard(ZX_to_Flx(a4,pp),ZX_to_Flx(a6,pp),ZX_to_Flx(T,pp),pp);
2197 : }
2198 : else
2199 : {
2200 1473 : GEN J = FpXQ_ellj(a4,a6,T,p), M;
2201 1473 : if (degpol(J) <= 0)
2202 1446 : return FpXQ_ellcardj(a4,a6,constant_coeff(J),T,q,p,n);
2203 27 : M = FpXQ_minpoly(J,T,p);
2204 27 : if (degpol(M) < degpol(T))
2205 7 : return FpXQ_ffellcard(a4, a6, M, q, T, p, n);
2206 20 : return Fq_ellcard_SEA(a4, a6, q, T, p, 0);
2207 : }
2208 : }
2209 :
2210 : GEN
2211 7293 : FpXQ_ellcard(GEN a4, GEN a6, GEN T, GEN p)
2212 : {
2213 7293 : pari_sp av = avma;
2214 7293 : return gc_INT(av, FpXQ_ellcard_i(a4, a6, T, p));
2215 : }
2216 :
2217 : static GEN
2218 18 : _FpXQE_pairorder(void *E, GEN P, GEN Q, GEN m, GEN F)
2219 : {
2220 18 : struct _FpXQE *e = (struct _FpXQE *) E;
2221 18 : return FpXQ_order(FpXQE_weilpairing(P,Q,m,e->a4,e->T,e->p), F, e->T, e->p);
2222 : }
2223 :
2224 : GEN
2225 13 : FpXQ_ellgroup(GEN a4, GEN a6, GEN N, GEN T, GEN p, GEN *pt_m)
2226 : {
2227 : struct _FpXQE e;
2228 13 : GEN q = powiu(p, get_FpX_degree(T));
2229 13 : e.a4=a4; e.a6=a6; e.T=T; e.p=p;
2230 13 : return gen_ellgroup(N, subiu(q,1), pt_m, (void*)&e, &FpXQE_group, _FpXQE_pairorder);
2231 : }
2232 :
2233 : GEN
2234 7 : FpXQ_ellgens(GEN a4, GEN a6, GEN ch, GEN D, GEN m, GEN T, GEN p)
2235 : {
2236 : GEN P;
2237 7 : pari_sp av = avma;
2238 : struct _FpXQE e;
2239 7 : e.a4=a4; e.a6=a6; e.T=T; e.p=p;
2240 7 : switch(lg(D)-1)
2241 : {
2242 7 : case 1:
2243 7 : P = gen_gener(gel(D,1), (void*)&e, &FpXQE_group);
2244 7 : P = mkvec(FpXQE_changepoint(P, ch, T, p));
2245 7 : 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 7 : return gc_GEN(av, P);
2253 : }
|
