Line data Source code
1 : /* Copyright (C) 2014 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 : /* Not so fast arithmetic with points over elliptic curves over Fl */
19 :
20 : /***********************************************************************/
21 : /** Flj **/
22 : /***********************************************************************/
23 : /* Jacobian coordinates: we represent a projective point (x:y:z) on E by
24 : * [z*x, z^2*y, z]. Not the fastest representation available for the problem,
25 : * but easy to implement and up to 60% faster than the school-book method. */
26 :
27 : /* Cost: 1M + 8S + 1*a + 10add + 1*8 + 2*2 + 1*3.
28 : * http://www.hyperelliptic.org/EFD/g1p/auto-shortw-jacobian.html#doubling-dbl-2007-bl */
29 : INLINE void
30 26404747 : Flj_dbl_indir_pre(GEN P, GEN Q, ulong a4, ulong p, ulong pi)
31 : {
32 : ulong X1, Y1, Z1;
33 : ulong XX, YY, YYYY, ZZ, S, M, T;
34 :
35 26404747 : X1 = P[1]; Y1 = P[2]; Z1 = P[3];
36 :
37 26404747 : if (Z1 == 0) { Q[1] = X1; Q[2] = Y1; Q[3] = Z1; return; }
38 :
39 26235965 : XX = Fl_sqr_pre(X1, p, pi);
40 26235664 : YY = Fl_sqr_pre(Y1, p, pi);
41 26217238 : YYYY = Fl_sqr_pre(YY, p, pi);
42 26214350 : ZZ = Fl_sqr_pre(Z1, p, pi);
43 26210967 : S = Fl_double(Fl_sub(Fl_sqr_pre(Fl_add(X1, YY, p), p, pi),
44 : Fl_add(XX, YYYY, p), p), p);
45 26215100 : M = Fl_addmul_pre(Fl_triple(XX, p), a4, Fl_sqr_pre(ZZ, p, pi), p, pi);
46 26241584 : T = Fl_sub(Fl_sqr_pre(M, p, pi), Fl_double(S, p), p);
47 26229785 : Q[1] = T;
48 26229785 : Q[2] = Fl_sub(Fl_mul_pre(M, Fl_sub(S, T, p), p, pi),
49 : Fl_double(Fl_double(Fl_double(YYYY, p), p), p), p);
50 26221163 : Q[3] = Fl_sub(Fl_sqr_pre(Fl_add(Y1, Z1, p), p, pi),
51 : Fl_add(YY, ZZ, p), p);
52 : }
53 :
54 : INLINE void
55 21759525 : Flj_dbl_pre_inplace(GEN P, ulong a4, ulong p, ulong pi)
56 : {
57 21759525 : Flj_dbl_indir_pre(P, P, a4, p, pi);
58 21755772 : }
59 :
60 : GEN
61 4647804 : Flj_dbl_pre(GEN P, ulong a4, ulong p, ulong pi)
62 : {
63 4647804 : GEN Q = cgetg(4, t_VECSMALL);
64 4647409 : Flj_dbl_indir_pre(P, Q, a4, p, pi);
65 4645937 : return Q;
66 : }
67 :
68 : /* Cost: 11M + 5S + 9add + 4*2.
69 : * http://www.hyperelliptic.org/EFD/g1p/auto-shortw-jacobian.html#addition-add-2007-bl */
70 : INLINE void
71 8256894 : Flj_add_indir_pre(GEN P, GEN Q, GEN R, ulong a4, ulong p, ulong pi)
72 : {
73 : ulong X1, Y1, Z1, X2, Y2, Z2;
74 : ulong Z1Z1, Z2Z2, U1, U2, S1, S2, H, I, J, r, V, W;
75 8256894 : X1 = P[1]; Y1 = P[2]; Z1 = P[3];
76 8256894 : X2 = Q[1]; Y2 = Q[2]; Z2 = Q[3];
77 :
78 8256894 : if (Z2 == 0) { R[1] = X1; R[2] = Y1; R[3] = Z1; return; }
79 8256582 : if (Z1 == 0) { R[1] = X2; R[2] = Y2; R[3] = Z2; return; }
80 :
81 8239913 : Z1Z1 = Fl_sqr_pre(Z1, p, pi);
82 8240067 : Z2Z2 = Fl_sqr_pre(Z2, p, pi);
83 8237988 : U1 = Fl_mul_pre(X1, Z2Z2, p, pi);
84 8238674 : U2 = Fl_mul_pre(X2, Z1Z1, p, pi);
85 8239160 : S1 = Fl_mul_pre(Y1, Fl_mul_pre(Z2, Z2Z2, p, pi), p, pi);
86 8238865 : S2 = Fl_mul_pre(Y2, Fl_mul_pre(Z1, Z1Z1, p, pi), p, pi);
87 8238689 : H = Fl_sub(U2, U1, p);
88 8239188 : r = Fl_double(Fl_sub(S2, S1, p), p);
89 :
90 8241212 : if (H == 0) {
91 1132431 : if (r == 0) Flj_dbl_indir_pre(P, R, a4, p, pi); /* P = Q */
92 1125621 : else { R[1] = R[2] = 1; R[3] = 0; } /* P = -Q */
93 1132430 : return;
94 : }
95 7108781 : I = Fl_sqr_pre(Fl_double(H, p), p, pi);
96 7108803 : J = Fl_mul_pre(H, I, p, pi);
97 7108602 : V = Fl_mul_pre(U1, I, p, pi);
98 7108293 : W = Fl_sub(Fl_sqr_pre(r, p, pi), Fl_add(J, Fl_double(V, p), p), p);
99 7108313 : R[1] = W;
100 7108313 : R[2] = Fl_sub(Fl_mul_pre(r, Fl_sub(V, W, p), p, pi),
101 : Fl_double(Fl_mul_pre(S1, J, p, pi), p), p);
102 7108419 : R[3] = Fl_mul_pre(Fl_sub(Fl_sqr_pre(Fl_add(Z1, Z2, p), p, pi),
103 : Fl_add(Z1Z1, Z2Z2, p), p), H, p, pi);
104 : }
105 :
106 : INLINE void
107 8164526 : Flj_add_pre_inplace(GEN P, GEN Q, ulong a4, ulong p, ulong pi)
108 8164526 : { Flj_add_indir_pre(P, Q, P, a4, p, pi); }
109 :
110 : GEN
111 91285 : Flj_add_pre(GEN P, GEN Q, ulong a4, ulong p, ulong pi)
112 : {
113 91285 : GEN R = cgetg(4, t_VECSMALL);
114 91285 : Flj_add_indir_pre(P, Q, R, a4, p, pi);
115 91285 : return R;
116 : }
117 :
118 : GEN
119 2882344 : Flj_neg(GEN Q, ulong p)
120 2882344 : { return mkvecsmall3(Q[1], Fl_neg(Q[2], p), Q[3]); }
121 :
122 : typedef struct {
123 : ulong pbits, nbits; /* Positive bits and negative bits */
124 : ulong lnzb; /* Leading nonzero bit */
125 : } naf_t;
126 :
127 : /* Return the signed binary representation (i.e. the Non-Adjacent Form
128 : * in base 2) of a; a = x.pbits - x.nbits (+ 2^BILif < 0; this
129 : * exceptional case can happen if a > 2^(BIL-1)) */
130 : static void
131 5177797 : naf_repr(naf_t *x, ulong a)
132 : {
133 5177797 : ulong pbits = 0, nbits = 0, c0 = 0, c1, a0;
134 : long t, i;
135 :
136 47117405 : for (i = 0; a; a >>= 1, ++i) {
137 41939608 : a0 = a & 1;
138 41939608 : c1 = (c0 + a0 + ((a & 2) >> 1)) >> 1;
139 41939608 : t = c0 + a0 - (c1 << 1);
140 41939608 : if (t < 0)
141 6703956 : nbits |= (1UL << i);
142 35235652 : else if (t > 0)
143 7639000 : pbits |= (1UL << i);
144 41939608 : c0 = c1;
145 : }
146 5177797 : c1 = c0 >> 1;
147 5177797 : t = c0 - (c1 << 1);
148 : /* since a >= 0, we have t >= 0; if i = BIL, pbits (virtually) overflows;
149 : * that leading overflowed bit is implied and not recorded in pbits */
150 5177797 : if (t > 0 && i != BITS_IN_LONG) pbits |= (1UL << i);
151 5177797 : x->pbits = pbits;
152 5177797 : x->nbits = nbits;
153 5177797 : x->lnzb = t? i-2: i-3;
154 5177797 : }
155 :
156 : /* Standard left-to-right signed double-and-add to compute [n]P. */
157 : static GEN
158 4709381 : Flj_mulu_pre_naf(GEN P, ulong n, ulong a4, ulong p, ulong pi, const naf_t *x)
159 : {
160 : GEN R, Pi;
161 : ulong pbits, nbits;
162 : ulong m;
163 :
164 4709381 : if (n == 0) return mkvecsmall3(1, 1, 0);
165 4709381 : if (n == 1) return Flv_copy(P);
166 :
167 4647844 : R = Flj_dbl_pre(P, a4, p, pi);
168 4645895 : if (n == 2) return R;
169 :
170 3797824 : pbits = x->pbits;
171 3797824 : nbits = x->nbits;
172 3797824 : Pi = nbits? Flj_neg(P, p): NULL;
173 3798779 : m = (1UL << x->lnzb);
174 25557821 : for ( ; m; m >>= 1) {
175 21758301 : Flj_dbl_pre_inplace(R, a4, p, pi);
176 21755701 : if (m & pbits)
177 3508830 : Flj_add_pre_inplace(R, P, a4, p, pi);
178 18246871 : else if (m & nbits)
179 4659407 : Flj_add_pre_inplace(R, Pi, a4, p, pi);
180 : }
181 3799520 : return gc_const((pari_sp)R, R);
182 : }
183 :
184 : GEN
185 3462042 : Flj_mulu_pre(GEN P, ulong n, ulong a4, ulong p, ulong pi)
186 : {
187 3462042 : naf_t x; naf_repr(&x, n);
188 3462692 : return Flj_mulu_pre_naf(P, n, a4, p, pi, &x);
189 : }
190 :
191 : struct _Flj { ulong a4, p, pi; };
192 :
193 : static GEN
194 91286 : _Flj_add(void *E, GEN P, GEN Q)
195 : {
196 91286 : struct _Flj *ell=(struct _Flj *) E;
197 91286 : return Flj_add_pre(P, Q, ell->a4, ell->p, ell->pi);
198 : }
199 :
200 : static GEN
201 106607 : _Flj_mul(void *E, GEN P, GEN n)
202 : {
203 106607 : struct _Flj *ell = (struct _Flj *) E;
204 106607 : long s = signe(n);
205 : GEN Q;
206 106607 : if (s==0) return mkvecsmall3(1, 1, 0);
207 106607 : Q = Flj_mulu_pre(P, itou(n), ell->a4, ell->p, ell->pi);
208 106609 : return s>0 ? Q : Flj_neg(Q, ell->p);
209 : }
210 : static GEN
211 0 : _Flj_one(void *E)
212 0 : { (void) E; return mkvecsmall3(1, 1, 0); }
213 :
214 : GEN
215 15322 : FljV_factorback_pre(GEN P, GEN L, ulong a4, ulong p, ulong pi)
216 : {
217 : struct _Flj E;
218 15322 : E.a4 = a4; E.p = p; E.pi = pi;
219 15322 : return gen_factorback(P, L, (void*)&E, &_Flj_add, &_Flj_mul, &_Flj_one);
220 : }
221 :
222 : ulong
223 294603 : Flj_order_ufact(GEN P, ulong n, GEN fa, ulong a4, ulong p, ulong pi)
224 : {
225 294603 : pari_sp av = avma;
226 294603 : GEN T = gel(fa,1), E = gel(fa,2);
227 294603 : long i, l = lg(T);
228 294603 : ulong res = 1;
229 :
230 831160 : for (i = 1; i < l; i++, set_avma(av))
231 : {
232 664351 : ulong j, t = T[i], e = E[i];
233 664351 : GEN b = P;
234 664351 : naf_t x; naf_repr(&x, t);
235 664446 : if (l != 2) b = Flj_mulu_pre(b, n / upowuu(t,e), a4, p, pi);
236 1911140 : for (j = 0; j < e && b[3]; j++) b = Flj_mulu_pre_naf(b, t, a4, p, pi, &x);
237 664274 : if (b[3]) return 0;
238 536472 : res *= upowuu(t, j);
239 : }
240 166809 : return res;
241 : }
242 :
243 : GEN
244 2270219 : Fle_to_Flj(GEN P)
245 4540177 : { return ell_is_inf(P) ? mkvecsmall3(1UL, 1UL, 0UL):
246 2270154 : mkvecsmall3(P[1], P[2], 1UL); }
247 :
248 : GEN
249 68257 : Flj_to_Fle(GEN P, ulong p)
250 : {
251 68257 : if (P[3] == 0) return ellinf();
252 : else
253 : {
254 67557 : ulong Z = Fl_inv(P[3], p);
255 67557 : ulong Z2 = Fl_sqr(Z, p);
256 67557 : ulong X3 = Fl_mul(P[1], Z2, p);
257 67557 : ulong Y3 = Fl_mul(P[2], Fl_mul(Z, Z2, p), p);
258 67557 : return mkvecsmall2(X3, Y3);
259 : }
260 : }
261 :
262 : GEN
263 2285210 : Flj_to_Fle_pre(GEN P, ulong p, ulong pi)
264 : {
265 2285210 : if (P[3] == 0) return ellinf();
266 : else
267 : {
268 1582939 : ulong Z = Fl_inv(P[3], p);
269 1583636 : ulong Z2 = Fl_sqr_pre(Z, p, pi);
270 1583394 : ulong X3 = Fl_mul_pre(P[1], Z2, p, pi);
271 1583298 : ulong Y3 = Fl_mul_pre(P[2], Fl_mul_pre(Z, Z2, p, pi), p, pi);
272 1583206 : return mkvecsmall2(X3, Y3);
273 : }
274 : }
275 :
276 : INLINE void
277 8284416 : random_Fle_pre_indir(ulong a4, ulong a6, ulong p, ulong pi,
278 : ulong *pt_x, ulong *pt_y)
279 : {
280 : ulong x, x2, y, rhs;
281 : do
282 : {
283 8284416 : x = random_Fl(p); /* x^3+a4*x+a6 = x*(x^2+a4)+a6 */
284 8292794 : x2 = Fl_sqr_pre(x, p, pi);
285 8271257 : rhs = Fl_addmul_pre(a6, x, Fl_add(x2, a4, p), p, pi);
286 8273151 : } while ((!rhs && !Fl_add(Fl_triple(x2,p),a4,p)) || krouu(rhs, p) < 0);
287 4150627 : y = Fl_sqrt_pre(rhs, p, pi);
288 4151435 : *pt_x = x; *pt_y = y;
289 4151435 : }
290 :
291 : GEN
292 500226 : random_Flj_pre(ulong a4, ulong a6, ulong p, ulong pi)
293 : {
294 : ulong x, y;
295 500226 : random_Fle_pre_indir(a4, a6, p, pi, &x, &y);
296 500213 : return mkvecsmall3(x, y, 1);
297 : }
298 :
299 : GEN
300 153237 : Flj_changepointinv_pre(GEN P, GEN ch, ulong p, ulong pi)
301 : {
302 : ulong c, u, r, s, t, u2, u3, z2;
303 153237 : ulong x = uel(P,1), y = uel(P,2), z = uel(P,3);
304 : GEN w;
305 153237 : if (z == 0) return Flv_copy(P);
306 153209 : u = ch[1]; r = ch[2];
307 153209 : s = ch[3]; t = ch[4];
308 153209 : u2 = Fl_sqr_pre(u, p, pi); u3 = Fl_mul_pre(u, u2, p, pi);
309 153211 : c = Fl_mul_pre(u2, x, p, pi);
310 153207 : z2 = Fl_sqr_pre(z, p, pi);
311 153208 : w = cgetg(4, t_VECSMALL);
312 153210 : uel(w,1) = Fl_add(c, Fl_mul_pre(r, z2, p, pi), p);
313 153210 : uel(w,2) = Fl_add(Fl_mul_pre(u3 ,y, p, pi),
314 : Fl_mul_pre(z, Fl_add(Fl_mul_pre(s,c,p,pi),
315 : Fl_mul_pre(z2,t,p,pi), p), p, pi), p);
316 153208 : uel(w,3) = z;
317 153208 : return w;
318 : }
319 :
320 : /***********************************************************************/
321 : /** Fle **/
322 : /***********************************************************************/
323 : GEN
324 16041 : Fle_changepoint(GEN P, GEN ch, ulong p)
325 : {
326 : ulong c, u, r, s, t, v, v2, v3;
327 : GEN z;
328 16041 : if (ell_is_inf(P)) return ellinf();
329 16041 : u = ch[1]; r = ch[2];
330 16041 : s = ch[3]; t = ch[4];
331 16041 : v = Fl_inv(u, p); v2 = Fl_sqr(v,p); v3 = Fl_mul(v,v2,p);
332 16041 : c = Fl_sub(uel(P,1),r,p);
333 16041 : z = cgetg(3,t_VECSMALL);
334 16041 : z[1] = Fl_mul(v2, c, p);
335 16041 : z[2] = Fl_mul(v3, Fl_sub(uel(P,2), Fl_add(Fl_mul(s,c, p),t, p),p),p);
336 16041 : return z;
337 : }
338 :
339 : GEN
340 134231 : Fle_changepointinv(GEN P, GEN ch, ulong p)
341 : {
342 : ulong c, u, r, s, t, u2, u3;
343 : GEN z;
344 134231 : if (ell_is_inf(P)) return ellinf();
345 133531 : u = ch[1]; r = ch[2];
346 133531 : s = ch[3]; t = ch[4];
347 133531 : u2 = Fl_sqr(u, p); u3 = Fl_mul(u,u2,p);
348 133531 : c = Fl_mul(u2,uel(P,1), p);
349 133531 : z = cgetg(3, t_VECSMALL);
350 133531 : z[1] = Fl_add(c,r,p);
351 133531 : z[2] = Fl_add(Fl_mul(u3,uel(P,2),p), Fl_add(Fl_mul(s,c,p), t, p), p);
352 133531 : return z;
353 : }
354 : static GEN
355 826673 : Fle_dbl_slope(GEN P, ulong a4, ulong p, ulong *slope)
356 : {
357 : ulong x, y, Qx, Qy;
358 826673 : if (ell_is_inf(P) || !P[2]) return ellinf();
359 686181 : x = P[1]; y = P[2];
360 686181 : *slope = Fl_div(Fl_add(Fl_triple(Fl_sqr(x,p), p), a4, p),
361 : Fl_double(y, p), p);
362 686192 : Qx = Fl_sub(Fl_sqr(*slope, p), Fl_double(x, p), p);
363 686186 : Qy = Fl_sub(Fl_mul(*slope, Fl_sub(x, Qx, p), p), y, p);
364 686177 : return mkvecsmall2(Qx, Qy);
365 : }
366 :
367 : GEN
368 580410 : Fle_dbl(GEN P, ulong a4, ulong p)
369 : {
370 : ulong slope;
371 580410 : return Fle_dbl_slope(P,a4,p,&slope);
372 : }
373 :
374 : static GEN
375 1518834 : Fle_add_slope(GEN P, GEN Q, ulong a4, ulong p, ulong *slope)
376 : {
377 : ulong Px, Py, Qx, Qy, Rx, Ry;
378 1518834 : if (ell_is_inf(P)) return Q;
379 1518831 : if (ell_is_inf(Q)) return P;
380 1518818 : Px = P[1]; Py = P[2];
381 1518818 : Qx = Q[1]; Qy = Q[2];
382 1518818 : if (Px==Qx) return Py==Qy ? Fle_dbl_slope(P, a4, p, slope): ellinf();
383 1384969 : *slope = Fl_div(Fl_sub(Py, Qy, p), Fl_sub(Px, Qx, p), p);
384 1385131 : Rx = Fl_sub(Fl_sub(Fl_sqr(*slope, p), Px, p), Qx, p);
385 1385079 : Ry = Fl_sub(Fl_mul(*slope, Fl_sub(Px, Rx, p), p), Py, p);
386 1385074 : return mkvecsmall2(Rx, Ry);
387 : }
388 :
389 : GEN
390 1506344 : Fle_add(GEN P, GEN Q, ulong a4, ulong p)
391 : {
392 : ulong slope;
393 1506344 : return Fle_add_slope(P,Q,a4,p,&slope);
394 : }
395 :
396 : static GEN
397 131208 : Fle_neg(GEN P, ulong p)
398 : {
399 131208 : if (ell_is_inf(P)) return P;
400 131208 : return mkvecsmall2(P[1], Fl_neg(P[2], p));
401 : }
402 :
403 : GEN
404 0 : Fle_sub(GEN P, GEN Q, ulong a4, ulong p)
405 : {
406 0 : pari_sp av = avma;
407 : ulong slope;
408 0 : return gerepileupto(av, Fle_add_slope(P, Fle_neg(Q, p), a4, p, &slope));
409 : }
410 :
411 : struct _Fle { ulong a4, a6, p; };
412 :
413 : static GEN
414 0 : _Fle_dbl(void *E, GEN P)
415 : {
416 0 : struct _Fle *ell = (struct _Fle *) E;
417 0 : return Fle_dbl(P, ell->a4, ell->p);
418 : }
419 :
420 : static GEN
421 365533 : _Fle_add(void *E, GEN P, GEN Q)
422 : {
423 365533 : struct _Fle *ell=(struct _Fle *) E;
424 365533 : return Fle_add(P, Q, ell->a4, ell->p);
425 : }
426 :
427 : GEN
428 2827920 : Fle_mulu(GEN P, ulong n, ulong a4, ulong p)
429 : {
430 : ulong pi;
431 2827920 : if (!n || ell_is_inf(P)) return ellinf();
432 2827838 : if (n==1) return zv_copy(P);
433 2815112 : if (n==2) return Fle_dbl(P, a4, p);
434 2269982 : pi = get_Fl_red(p);
435 2270227 : return Flj_to_Fle_pre(Flj_mulu_pre(Fle_to_Flj(P), n, a4, p, pi), p, pi);
436 : }
437 :
438 : static GEN
439 2345713 : _Fle_mul(void *E, GEN P, GEN n)
440 : {
441 2345713 : pari_sp av = avma;
442 2345713 : struct _Fle *e=(struct _Fle *) E;
443 2345713 : long s = signe(n);
444 : GEN Q;
445 2345713 : if (!s || ell_is_inf(P)) return ellinf();
446 2326865 : if (s < 0) P = Fle_neg(P, e->p);
447 2326865 : if (is_pm1(n)) return s > 0? zv_copy(P): P;
448 1973569 : Q = (lgefint(n)==3) ? Fle_mulu(P, uel(n,2), e->a4, e->p):
449 0 : gen_pow(P, n, (void*)e, &_Fle_dbl, &_Fle_add);
450 1973405 : return s > 0? Q: gerepileuptoleaf(av, Q);
451 : }
452 :
453 : GEN
454 28959 : Fle_mul(GEN P, GEN n, ulong a4, ulong p)
455 : {
456 : struct _Fle E;
457 28959 : E.a4 = a4; E.p = p;
458 28959 : return _Fle_mul(&E, P, n);
459 : }
460 :
461 : /* Finds a random nonsingular point on E */
462 : GEN
463 3649637 : random_Fle_pre(ulong a4, ulong a6, ulong p, ulong pi)
464 : {
465 : ulong x, y;
466 3649637 : random_Fle_pre_indir(a4, a6, p, pi, &x, &y);
467 3651588 : return mkvecsmall2(x, y);
468 : }
469 :
470 : GEN
471 26040 : random_Fle(ulong a4, ulong a6, ulong p)
472 26040 : { return random_Fle_pre(a4, a6, p, get_Fl_red(p)); }
473 :
474 : static GEN
475 0 : _Fle_rand(void *E)
476 : {
477 0 : struct _Fle *e=(struct _Fle *) E;
478 0 : return random_Fle(e->a4, e->a6, e->p);
479 : }
480 :
481 : static const struct bb_group Fle_group={_Fle_add,_Fle_mul,_Fle_rand,hash_GEN,zv_equal,ell_is_inf,NULL};
482 :
483 : GEN
484 288788 : Fle_order(GEN z, GEN o, ulong a4, ulong p)
485 : {
486 288788 : pari_sp av = avma;
487 : struct _Fle e;
488 288788 : e.a4=a4;
489 288788 : e.p=p;
490 288788 : return gerepileuptoint(av, gen_order(z, o, (void*)&e, &Fle_group));
491 : }
492 :
493 : GEN
494 54327 : Fle_log(GEN a, GEN b, GEN o, ulong a4, ulong p)
495 : {
496 54327 : pari_sp av = avma;
497 : struct _Fle e;
498 54327 : e.a4=a4;
499 54327 : e.p=p;
500 54327 : return gerepileuptoint(av, gen_PH_log(a, b, o, (void*)&e, &Fle_group));
501 : }
502 :
503 : ulong
504 0 : Fl_ellj(ulong a4, ulong a6, ulong p)
505 : {
506 0 : if (SMALL_ULONG(p))
507 : { /* a43 = 4 a4^3 */
508 0 : ulong a43 = Fl_double(Fl_double(Fl_mul(a4, Fl_sqr(a4, p), p), p), p);
509 : /* a62 = 27 a6^2 */
510 0 : ulong a62 = Fl_mul(Fl_sqr(a6, p), 27 % p, p);
511 0 : ulong z1 = Fl_mul(a43, 1728 % p, p);
512 0 : ulong z2 = Fl_add(a43, a62, p);
513 0 : return Fl_div(z1, z2, p);
514 : }
515 0 : return Fl_ellj_pre(a4, a6, p, get_Fl_red(p));
516 : }
517 :
518 : void
519 176118 : Fl_ellj_to_a4a6(ulong j, ulong p, ulong *pt_a4, ulong *pt_a6)
520 : {
521 176118 : ulong zagier = 1728 % p;
522 176118 : if (j == 0) { *pt_a4 = 0; *pt_a6 =1; }
523 176104 : else if (j == zagier) { *pt_a4 = 1; *pt_a6 =0; }
524 : else
525 : {
526 176090 : ulong k = Fl_sub(zagier, j, p);
527 176089 : ulong kj = Fl_mul(k, j, p);
528 176090 : ulong k2j = Fl_mul(kj, k, p);
529 176077 : *pt_a4 = Fl_triple(kj, p);
530 176075 : *pt_a6 = Fl_double(k2j, p);
531 : }
532 176098 : }
533 :
534 : ulong
535 2118669 : Fl_elldisc_pre(ulong a4, ulong a6, ulong p, ulong pi)
536 : { /* D = -(4A^3 + 27B^2) */
537 : ulong t1, t2;
538 2118669 : t1 = Fl_mul_pre(a4, Fl_sqr_pre(a4, p, pi), p, pi);
539 2116534 : t1 = Fl_double(Fl_double(t1, p), p);
540 2115988 : t2 = Fl_mul_pre(27 % p, Fl_sqr_pre(a6, p, pi), p, pi);
541 2116602 : return Fl_neg(Fl_add(t1, t2, p), p);
542 : }
543 :
544 : ulong
545 0 : Fl_elldisc(ulong a4, ulong a6, ulong p)
546 : {
547 0 : if (SMALL_ULONG(p))
548 : { /* D = -(4A^3 + 27B^2) */
549 : ulong t1, t2;
550 0 : t1 = Fl_mul(a4, Fl_sqr(a4, p), p);
551 0 : t1 = Fl_double(Fl_double(t1, p), p);
552 0 : t2 = Fl_mul(27 % p, Fl_sqr(a6, p), p);
553 0 : return Fl_neg(Fl_add(t1, t2, p), p);
554 : }
555 0 : return Fl_elldisc_pre(a4, a6, p, get_Fl_red(p));
556 : }
557 :
558 : void
559 107856 : Fl_elltwist_disc(ulong a4, ulong a6, ulong D, ulong p, ulong *pa4, ulong *pa6)
560 : {
561 107856 : ulong D2 = Fl_sqr(D, p);
562 107855 : *pa4 = Fl_mul(a4, D2, p);
563 107855 : *pa6 = Fl_mul(a6, Fl_mul(D, D2, p), p);
564 107856 : }
565 :
566 : void
567 0 : Fl_elltwist(ulong a4, ulong a6, ulong p, ulong *pt_a4, ulong *pt_a6)
568 0 : { Fl_elltwist_disc(a4, a6, nonsquare_Fl(p), p, pt_a4, pt_a6); }
569 :
570 : static void
571 52505059 : Fle_dbl_sinv_pre_inplace(GEN P, ulong a4, ulong sinv, ulong p, ulong pi)
572 : {
573 : ulong x, y, slope;
574 52505059 : if (uel(P,1)==p) return;
575 52259939 : if (!P[2]) { P[1] = p; return; }
576 52110534 : x = P[1]; y = P[2];
577 52110534 : slope = Fl_mul_pre(Fl_add(Fl_triple(Fl_sqr_pre(x, p, pi), p), a4, p),
578 : sinv, p, pi);
579 52083572 : P[1] = Fl_sub(Fl_sqr_pre(slope, p, pi), Fl_double(x, p), p);
580 51986940 : P[2] = Fl_sub(Fl_mul_pre(slope, Fl_sub(x, P[1], p), p, pi), y, p);
581 : }
582 :
583 : static void
584 7208536 : Fle_add_sinv_pre_inplace(GEN P, GEN Q, ulong a4, ulong sinv, ulong p, ulong pi)
585 : {
586 : ulong Px, Py, Qx, Qy, slope;
587 7208536 : if (uel(P,1)==p) { P[1] = Q[1]; P[2] = Q[2]; }
588 7208536 : if (ell_is_inf(Q)) return;
589 7210086 : Px = P[1]; Py = P[2];
590 7210086 : Qx = Q[1]; Qy = Q[2];
591 7210086 : if (Px==Qx)
592 : {
593 29371 : if (Py==Qy) Fle_dbl_sinv_pre_inplace(P, a4, sinv, p, pi);
594 14068 : else P[1] = p;
595 29372 : return;
596 : }
597 7180715 : slope = Fl_mul_pre(Fl_sub(Py, Qy, p), sinv, p, pi);
598 7171356 : P[1] = Fl_sub(Fl_sub(Fl_sqr_pre(slope, p, pi), Px, p), Qx, p);
599 7160506 : P[2] = Fl_sub(Fl_mul_pre(slope, Fl_sub(Px, P[1], p), p, pi), Py, p);
600 : }
601 :
602 : static void
603 8061755 : Fle_sub_sinv_pre_inplace(GEN P, GEN Q, ulong a4, ulong sinv, ulong p, ulong pi)
604 : {
605 : ulong Px, Py, Qx, Qy, slope;
606 8061755 : if (uel(P,1)==p) { P[1] = Q[1]; P[2] = Fl_neg(Q[2], p); }
607 8061755 : if (ell_is_inf(Q)) return;
608 8064239 : Px = P[1]; Py = P[2];
609 8064239 : Qx = Q[1]; Qy = Q[2];
610 8064239 : if (Px==Qx)
611 : {
612 35091 : if (Py==Qy) P[1] = p;
613 : else
614 13719 : Fle_dbl_sinv_pre_inplace(P, a4, sinv, p, pi);
615 35091 : return;
616 : }
617 8029148 : slope = Fl_mul_pre(Fl_add(Py, Qy, p), sinv, p, pi);
618 8015375 : P[1] = Fl_sub(Fl_sub(Fl_sqr_pre(slope, p, pi), Px, p), Qx, p);
619 8004217 : P[2] = Fl_sub(Fl_mul_pre(slope, Fl_sub(Px, P[1], p), p, pi), Py, p);
620 : }
621 :
622 : static long
623 68051327 : skipzero(long n) { return n ? n:1; }
624 :
625 : void
626 1788321 : FleV_add_pre_inplace(GEN P, GEN Q, GEN a4, ulong p, ulong pi)
627 : {
628 1788321 : long i, l=lg(a4);
629 1788321 : GEN sinv = cgetg(l, t_VECSMALL);
630 9018384 : for(i=1; i<l; i++)
631 7228987 : uel(sinv,i) = umael(P,i,1)==p? 1: skipzero(Fl_sub(mael(P,i,1), mael(Q,i,1), p));
632 1789397 : Flv_inv_pre_inplace(sinv, p, pi);
633 8993326 : for (i=1; i<l; i++)
634 7207654 : Fle_add_sinv_pre_inplace(gel(P,i), gel(Q,i), uel(a4,i), uel(sinv,i), p, pi);
635 1785672 : }
636 :
637 : void
638 2067201 : FleV_sub_pre_inplace(GEN P, GEN Q, GEN a4, ulong p, ulong pi)
639 : {
640 2067201 : long i, l=lg(a4);
641 2067201 : GEN sinv = cgetg(l, t_VECSMALL);
642 10150319 : for(i=1; i<l; i++)
643 8083494 : uel(sinv,i) = umael(P,i,1)==p? 1: skipzero(Fl_sub(mael(P,i,1), mael(Q,i,1), p));
644 2066825 : Flv_inv_pre_inplace(sinv, p, pi);
645 10125384 : for (i=1; i<l; i++)
646 8060940 : Fle_sub_sinv_pre_inplace(gel(P,i), gel(Q,i), uel(a4,i), uel(sinv,i), p, pi);
647 2064444 : }
648 :
649 : void
650 13557813 : FleV_dbl_pre_inplace(GEN P, GEN a4, ulong p, ulong pi)
651 : {
652 13557813 : long i, l=lg(a4);
653 13557813 : GEN sinv = cgetg(l, t_VECSMALL);
654 66982649 : for(i=1; i<l; i++)
655 53439006 : uel(sinv,i) = umael(P,i,1)==p? 1: skipzero(Fl_double(umael(P,i,2), p));
656 13543643 : Flv_inv_pre_inplace(sinv, p, pi);
657 66038528 : for(i=1; i<l; i++)
658 52501799 : Fle_dbl_sinv_pre_inplace(gel(P,i), uel(a4,i), uel(sinv,i), p, pi);
659 13536729 : }
660 :
661 : static void
662 1052359 : FleV_mulu_pre_naf_inplace(GEN P, ulong n, GEN a4, ulong p, ulong pi, const naf_t *x)
663 : {
664 1052359 : pari_sp av = avma;
665 : ulong pbits, nbits, m;
666 : GEN R;
667 1052359 : if (n == 1) return;
668 :
669 1052359 : R = P; P = gcopy(P);
670 1052153 : FleV_dbl_pre_inplace(R, a4, p, pi);
671 1050851 : if (n == 2) return;
672 :
673 1050757 : pbits = x->pbits;
674 1050757 : nbits = x->nbits;
675 1050757 : m = (1UL << x->lnzb);
676 13568354 : for ( ; m; m >>= 1) {
677 12517430 : FleV_dbl_pre_inplace(R, a4, p, pi);
678 12509946 : if (m & pbits)
679 1788320 : FleV_add_pre_inplace(R, P, a4, p, pi);
680 10721626 : else if (m & nbits)
681 2067207 : FleV_sub_pre_inplace(R, P, a4, p, pi);
682 : }
683 1050924 : set_avma(av);
684 : }
685 :
686 : void
687 1051922 : FleV_mulu_pre_inplace(GEN P, ulong n, GEN a4, ulong p, ulong pi)
688 : {
689 1051922 : naf_t x; naf_repr(&x, n);
690 1052367 : FleV_mulu_pre_naf_inplace(P, n, a4, p, pi, &x);
691 1050921 : }
692 :
693 : /***********************************************************************/
694 : /** Pairings **/
695 : /** Derived from APIP by Jerome Milan, 2012 **/
696 : /***********************************************************************/
697 : static ulong
698 403230 : Fle_vert(GEN P, GEN Q, ulong a4, ulong p)
699 : {
700 403230 : if (ell_is_inf(P))
701 151503 : return 1;
702 251727 : if (uel(Q, 1) != uel(P, 1))
703 232512 : return Fl_sub(uel(Q, 1), uel(P, 1), p);
704 19215 : if (uel(P,2) != 0 ) return 1;
705 13923 : return Fl_inv(Fl_add(Fl_triple(Fl_sqr(uel(P,1),p), p), a4, p), p);
706 : }
707 :
708 : static ulong
709 130810 : Fle_Miller_line(GEN R, GEN Q, ulong slope, ulong a4, ulong p)
710 : {
711 130810 : ulong x = uel(Q, 1), y = uel(Q, 2);
712 130810 : ulong tmp1 = Fl_sub(x, uel(R, 1), p);
713 130810 : ulong tmp2 = Fl_add(Fl_mul(tmp1, slope, p), uel(R,2), p);
714 130810 : if (y != tmp2)
715 120254 : return Fl_sub(y, tmp2, p);
716 10556 : if (y == 0)
717 6678 : return 1;
718 : else
719 : {
720 : ulong s1, s2;
721 3878 : ulong y2i = Fl_inv(Fl_double(y, p), p);
722 3878 : s1 = Fl_mul(Fl_add(Fl_triple(Fl_sqr(x, p), p), a4, p), y2i, p);
723 3878 : if (s1 != slope)
724 2108 : return Fl_sub(s1, slope, p);
725 1770 : s2 = Fl_mul(Fl_sub(Fl_triple(x, p), Fl_sqr(s1, p), p), y2i, p);
726 1770 : return s2 ? s2: y2i;
727 : }
728 : }
729 :
730 : /* Computes the equation of the line tangent to R and returns its
731 : * evaluation at the point Q. Also doubles the point R. */
732 : static ulong
733 250299 : Fle_tangent_update(GEN R, GEN Q, ulong a4, ulong p, GEN *pt_R)
734 : {
735 250299 : if (ell_is_inf(R)) { *pt_R = ellinf(); return 1; }
736 217949 : else if (uel(R,2) == 0) { *pt_R = ellinf(); return Fle_vert(R, Q, a4, p); }
737 : else
738 : {
739 : ulong slope;
740 118319 : *pt_R = Fle_dbl_slope(R, a4, p, &slope);
741 118319 : return Fle_Miller_line(R, Q, slope, a4, p);
742 : }
743 : }
744 :
745 : /* Computes the equation of the line through R and P, and returns its
746 : * evaluation at the point Q. Also adds P to the point R. */
747 : static ulong
748 32896 : Fle_chord_update(GEN R, GEN P, GEN Q, ulong a4, ulong p, GEN *pt_R)
749 : {
750 32896 : if (ell_is_inf(R)) { *pt_R = P; return Fle_vert(P, Q, a4, p); }
751 32014 : else if (ell_is_inf(P)) { *pt_R = R; return Fle_vert(R, Q, a4, p); }
752 32014 : else if (uel(P, 1) == uel(R, 1))
753 : {
754 19523 : if (uel(P, 2) == uel(R, 2)) return Fle_tangent_update(R, Q, a4, p, pt_R);
755 19523 : else { *pt_R = ellinf(); return Fle_vert(R, Q, a4, p); }
756 : }
757 : else
758 : {
759 : ulong slope;
760 12491 : *pt_R = Fle_add_slope(P, R, a4, p, &slope);
761 12491 : return Fle_Miller_line(R, Q, slope, a4, p);
762 : }
763 : }
764 :
765 : struct _Fle_miller { ulong p, a4; GEN P; };
766 : static GEN
767 250299 : Fle_Miller_dbl(void* E, GEN d)
768 : {
769 250299 : struct _Fle_miller *m = (struct _Fle_miller *)E;
770 250299 : ulong p = m->p, a4 = m->a4;
771 250299 : GEN P = m->P;
772 : ulong v, line;
773 250299 : ulong N = Fl_sqr(umael(d,1,1), p);
774 250299 : ulong D = Fl_sqr(umael(d,1,2), p);
775 250299 : GEN point = gel(d,2);
776 250299 : line = Fle_tangent_update(point, P, a4, p, &point);
777 250299 : N = Fl_mul(N, line, p);
778 250299 : v = Fle_vert(point, P, a4, p);
779 250299 : D = Fl_mul(D, v, p); return mkvec2(mkvecsmall2(N, D), point);
780 : }
781 : static GEN
782 32896 : Fle_Miller_add(void* E, GEN va, GEN vb)
783 : {
784 32896 : struct _Fle_miller *m = (struct _Fle_miller *)E;
785 32896 : ulong p = m->p, a4= m->a4;
786 32896 : GEN P = m->P, point;
787 : ulong v, line;
788 32896 : ulong na = umael(va,1,1), da = umael(va,1,2);
789 32896 : ulong nb = umael(vb,1,1), db = umael(vb,1,2);
790 32896 : GEN pa = gel(va,2), pb = gel(vb,2);
791 32896 : ulong N = Fl_mul(na, nb, p);
792 32896 : ulong D = Fl_mul(da, db, p);
793 32896 : line = Fle_chord_update(pa, pb, P, a4, p, &point);
794 32896 : N = Fl_mul(N, line, p);
795 32896 : v = Fle_vert(point, P, a4, p);
796 32896 : D = Fl_mul(D, v, p); return mkvec2(mkvecsmall2(N, D), point);
797 : }
798 :
799 : /* Returns the Miller function f_{m, Q} evaluated at the point P using
800 : * the standard Miller algorithm. */
801 : static ulong
802 118271 : Fle_Miller(GEN Q, GEN P, ulong m, ulong a4, ulong p)
803 : {
804 118271 : pari_sp av = avma;
805 : struct _Fle_miller d;
806 : GEN v;
807 : ulong N, D;
808 :
809 118271 : d.a4 = a4; d.p = p; d.P = P;
810 118271 : v = gen_powu_i(mkvec2(mkvecsmall2(1,1), Q), m, (void*)&d,
811 : Fle_Miller_dbl, Fle_Miller_add);
812 118271 : N = umael(v,1,1); D = umael(v,1,2);
813 118271 : return gc_ulong(av, Fl_div(N, D, p));
814 : }
815 :
816 : ulong
817 51824 : Fle_weilpairing(GEN P, GEN Q, ulong m, ulong a4, ulong p)
818 : {
819 51824 : pari_sp ltop = avma;
820 : ulong N, D, w;
821 51824 : if (ell_is_inf(P) || ell_is_inf(Q) || zv_equal(P,Q)) return 1;
822 50697 : N = Fle_Miller(P, Q, m, a4, p);
823 50697 : D = Fle_Miller(Q, P, m, a4, p);
824 50697 : w = Fl_div(N, D, p);
825 50697 : if (odd(m)) w = Fl_neg(w, p);
826 50697 : return gc_ulong(ltop, w);
827 : }
828 :
829 : ulong
830 16877 : Fle_tatepairing(GEN P, GEN Q, ulong m, ulong a4, ulong p)
831 : {
832 16877 : if (ell_is_inf(P) || ell_is_inf(Q)) return 1;
833 16877 : return Fle_Miller(P, Q, m, a4, p);
834 : }
835 :
836 : GEN
837 11438 : Fl_ellptors(ulong l, ulong N, ulong a4, ulong a6, ulong p)
838 : {
839 11438 : long v = z_lval(N, l);
840 : ulong N0, N1;
841 : GEN F;
842 11438 : if (v==0) return cgetg(1,t_VEC);
843 11438 : N0 = upowuu(l, v); N1 = N/N0;
844 11438 : F = mkmat2(mkcols(l), mkcols(v));
845 : while(1)
846 1582 : {
847 13020 : pari_sp av2 = avma;
848 : GEN P, Q;
849 : ulong d, s, t;
850 :
851 13020 : P = Fle_mulu(random_Fle(a4, a6, p), N1, a4, p);
852 13020 : Q = Fle_mulu(random_Fle(a4, a6, p), N1, a4, p);
853 13020 : s = itou(Fle_order(P, F, a4, p));
854 13020 : t = itou(Fle_order(Q, F, a4, p));
855 13020 : if (s < t) { swap(P,Q); lswap(s,t); }
856 13020 : if (s == N0) retmkvec(Fle_mulu(P, s/l, a4, p));
857 2394 : d = Fl_order(Fle_weilpairing(P, Q, s, a4, p), s, p);
858 2394 : if (s*d == N0)
859 812 : retmkvec2(Fle_mulu(P, s/l, a4, p), Fle_mulu(Q, t/l, a4, p));
860 1582 : set_avma(av2);
861 : }
862 : }
|
