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