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