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 : /** **/
22 : /** Flj **/
23 : /** **/
24 : /***********************************************************************/
25 :
26 : /* Arithmetic is implemented using Jacobian coordinates, representing
27 : * a projective point (x : y : z) on E by [z*x , z^2*y , z]. This is
28 : * probably not the fastest representation available for the given
29 : * problem, but they're easy to implement and up to 60% faster than
30 : * the school-book method. */
31 :
32 : /* Cost: 1M + 8S + 1*a + 10add + 1*8 + 2*2 + 1*3.
33 : * http://www.hyperelliptic.org/EFD/g1p/auto-shortw-jacobian.html#doubling-dbl-2007-bl */
34 : INLINE void
35 23678220 : Flj_dbl_indir_pre(GEN P, GEN Q, ulong a4, ulong p, ulong pi)
36 : {
37 : ulong X1, Y1, Z1;
38 : ulong XX, YY, YYYY, ZZ, S, M, T;
39 :
40 23678220 : X1 = P[1]; Y1 = P[2]; Z1 = P[3];
41 :
42 23678220 : if (Z1 == 0) { Q[1] = X1; Q[2] = Y1; Q[3] = Z1; return; }
43 :
44 23517642 : XX = Fl_sqr_pre(X1, p, pi);
45 23519681 : YY = Fl_sqr_pre(Y1, p, pi);
46 23516032 : YYYY = Fl_sqr_pre(YY, p, pi);
47 23515119 : ZZ = Fl_sqr_pre(Z1, p, pi);
48 23514069 : S = Fl_double(Fl_sub(Fl_sqr_pre(Fl_add(X1, YY, p), p, pi),
49 : Fl_add(XX, YYYY, p), p), p);
50 23509309 : M = Fl_addmul_pre(Fl_triple(XX, p), a4, Fl_sqr_pre(ZZ, p, pi), p, pi);
51 23519101 : T = Fl_sub(Fl_sqr_pre(M, p, pi), Fl_double(S, p), p);
52 23514853 : Q[1] = T;
53 23514853 : Q[2] = Fl_sub(Fl_mul_pre(M, Fl_sub(S, T, p), p, pi),
54 : Fl_double(Fl_double(Fl_double(YYYY, p), p), p), p);
55 23510805 : Q[3] = Fl_sub(Fl_sqr_pre(Fl_add(Y1, Z1, p), p, pi),
56 : Fl_add(YY, ZZ, p), p);
57 : }
58 :
59 : INLINE void
60 20044061 : Flj_dbl_pre_inplace(GEN P, ulong a4, ulong p, ulong pi)
61 : {
62 20044061 : Flj_dbl_indir_pre(P, P, a4, p, pi);
63 20044303 : }
64 :
65 : GEN
66 3630341 : Flj_dbl_pre(GEN P, ulong a4, ulong p, ulong pi)
67 : {
68 3630341 : GEN Q = cgetg(4, t_VECSMALL);
69 3630008 : Flj_dbl_indir_pre(P, Q, a4, p, pi);
70 3629905 : return Q;
71 : }
72 :
73 : /* Cost: 11M + 5S + 9add + 4*2.
74 : * http://www.hyperelliptic.org/EFD/g1p/auto-shortw-jacobian.html#addition-add-2007-bl */
75 : INLINE void
76 7190557 : Flj_add_indir_pre(GEN P, GEN Q, GEN R, ulong a4, ulong p, ulong pi)
77 : {
78 : ulong X1, Y1, Z1, X2, Y2, Z2;
79 : ulong Z1Z1, Z2Z2, U1, U2, S1, S2, H, I, J, r, V, W;
80 7190557 : X1 = P[1]; Y1 = P[2]; Z1 = P[3];
81 7190557 : X2 = Q[1]; Y2 = Q[2]; Z2 = Q[3];
82 :
83 7190557 : if (Z2 == 0) { R[1] = X1; R[2] = Y1; R[3] = Z1; return; }
84 7190382 : if (Z1 == 0) { R[1] = X2; R[2] = Y2; R[3] = Z2; return; }
85 :
86 7177847 : Z1Z1 = Fl_sqr_pre(Z1, p, pi);
87 7177984 : Z2Z2 = Fl_sqr_pre(Z2, p, pi);
88 7177561 : U1 = Fl_mul_pre(X1, Z2Z2, p, pi);
89 7177599 : U2 = Fl_mul_pre(X2, Z1Z1, p, pi);
90 7177705 : S1 = Fl_mul_pre(Y1, Fl_mul_pre(Z2, Z2Z2, p, pi), p, pi);
91 7177749 : S2 = Fl_mul_pre(Y2, Fl_mul_pre(Z1, Z1Z1, p, pi), p, pi);
92 7177912 : H = Fl_sub(U2, U1, p);
93 7177546 : r = Fl_double(Fl_sub(S2, S1, p), p);
94 :
95 7177454 : if (H == 0) {
96 545923 : if (r == 0) Flj_dbl_indir_pre(P, R, a4, p, pi); /* P = Q */
97 541261 : else { R[1] = R[2] = 1; R[3] = 0; } /* P = -Q */
98 545923 : return;
99 : }
100 6631531 : I = Fl_sqr_pre(Fl_double(H, p), p, pi);
101 6632112 : J = Fl_mul_pre(H, I, p, pi);
102 6632020 : V = Fl_mul_pre(U1, I, p, pi);
103 6632115 : W = Fl_sub(Fl_sqr_pre(r, p, pi), Fl_add(J, Fl_double(V, p), p), p);
104 6631847 : R[1] = W;
105 6631847 : R[2] = Fl_sub(Fl_mul_pre(r, Fl_sub(V, W, p), p, pi),
106 : Fl_double(Fl_mul_pre(S1, J, p, pi), p), p);
107 6632073 : R[3] = Fl_mul_pre(Fl_sub(Fl_sqr_pre(Fl_add(Z1, Z2, p), p, pi),
108 : Fl_add(Z1Z1, Z2Z2, p), p), H, p, pi);
109 : }
110 :
111 : INLINE void
112 7100255 : Flj_add_pre_inplace(GEN P, GEN Q, ulong a4, ulong p, ulong pi)
113 7100255 : { Flj_add_indir_pre(P, Q, P, a4, p, pi); }
114 :
115 : GEN
116 89508 : Flj_add_pre(GEN P, GEN Q, ulong a4, ulong p, ulong pi)
117 : {
118 89508 : GEN R = cgetg(4, t_VECSMALL);
119 89500 : Flj_add_indir_pre(P, Q, R, a4, p, pi);
120 89504 : return R;
121 : }
122 :
123 : GEN
124 2297678 : Flj_neg(GEN Q, ulong p)
125 2297678 : { return mkvecsmall3(Q[1], Fl_neg(Q[2], p), Q[3]); }
126 :
127 : typedef struct {
128 : ulong pbits, nbits; /* Positive bits and negative bits */
129 : ulong lnzb; /* Leading nonzero bit */
130 : } naf_t;
131 :
132 : /* Return the signed binary representation (i.e. the Non-Adjacent Form
133 : * in base 2) of a; a = x.pbits - x.nbits (+ 2^BILif < 0; this
134 : * exceptional case can happen if a > 2^(BIL-1)) */
135 : static void
136 3952815 : naf_repr(naf_t *x, ulong a)
137 : {
138 3952815 : ulong pbits = 0, nbits = 0, c0 = 0, c1, a0;
139 : long t, i;
140 :
141 38542469 : for (i = 0; a; a >>= 1, ++i) {
142 34589654 : a0 = a & 1;
143 34589654 : c1 = (c0 + a0 + ((a & 2) >> 1)) >> 1;
144 34589654 : t = c0 + a0 - (c1 << 1);
145 34589654 : if (t < 0)
146 5359880 : nbits |= (1UL << i);
147 29229774 : else if (t > 0)
148 6189075 : pbits |= (1UL << i);
149 34589654 : c0 = c1;
150 : }
151 3952815 : c1 = c0 >> 1;
152 3952815 : t = c0 - (c1 << 1);
153 : /* since a >= 0, we have t >= 0; if i = BIL, pbits (virtually) overflows;
154 : * that leading overflowed bit is implied and not recorded in pbits */
155 3952815 : if (t > 0 && i != BITS_IN_LONG) pbits |= (1UL << i);
156 3952815 : x->pbits = pbits;
157 3952815 : x->nbits = nbits;
158 3952815 : x->lnzb = t? i-2: i-3;
159 3952815 : }
160 :
161 : /* Standard left-to-right signed double-and-add to compute [n]P. */
162 : static GEN
163 3697429 : Flj_mulu_pre_naf(GEN P, ulong n, ulong a4, ulong p, ulong pi, const naf_t *x)
164 : {
165 : GEN R, Pi;
166 : ulong pbits, nbits;
167 : ulong m;
168 :
169 3697429 : if (n == 0) return mkvecsmall3(1, 1, 0);
170 3697372 : if (n == 1) return Flv_copy(P);
171 :
172 3630379 : R = Flj_dbl_pre(P, a4, p, pi);
173 3629860 : if (n == 2) return R;
174 :
175 2888393 : pbits = x->pbits;
176 2888393 : nbits = x->nbits;
177 2888393 : Pi = nbits? Flj_neg(P, p): NULL;
178 2888298 : m = (1UL << x->lnzb);
179 22933498 : for ( ; m; m >>= 1) {
180 20044992 : Flj_dbl_pre_inplace(R, a4, p, pi);
181 20043675 : if (m & pbits)
182 3126374 : Flj_add_pre_inplace(R, P, a4, p, pi);
183 16917301 : else if (m & nbits)
184 3976003 : Flj_add_pre_inplace(R, Pi, a4, p, pi);
185 : }
186 2888506 : return gc_const((pari_sp)R, R);
187 : }
188 :
189 : GEN
190 2604862 : Flj_mulu_pre(GEN P, ulong n, ulong a4, ulong p, ulong pi)
191 : {
192 2604862 : naf_t x; naf_repr(&x, n);
193 2605293 : return Flj_mulu_pre_naf(P, n, a4, p, pi, &x);
194 : }
195 :
196 : struct _Flj { ulong a4, p, pi; };
197 :
198 : static GEN
199 89509 : _Flj_add(void *E, GEN P, GEN Q)
200 : {
201 89509 : struct _Flj *ell=(struct _Flj *) E;
202 89509 : return Flj_add_pre(P, Q, ell->a4, ell->p, ell->pi);
203 : }
204 :
205 : static GEN
206 104619 : _Flj_mul(void *E, GEN P, GEN n)
207 : {
208 104619 : struct _Flj *ell = (struct _Flj *) E;
209 104619 : long s = signe(n);
210 : GEN Q;
211 104619 : if (s==0) return mkvecsmall3(1, 1, 0);
212 104619 : Q = Flj_mulu_pre(P, itou(n), ell->a4, ell->p, ell->pi);
213 104633 : return s>0 ? Q : Flj_neg(Q, ell->p);
214 : }
215 : static GEN
216 0 : _Flj_one(void *E)
217 0 : { (void) E; return mkvecsmall3(1, 1, 0); }
218 :
219 : GEN
220 15087 : FljV_factorback_pre(GEN P, GEN L, ulong a4, ulong p, ulong pi)
221 : {
222 : struct _Flj E;
223 15087 : E.a4 = a4; E.p = p; E.pi = pi;
224 15087 : return gen_factorback(P, L, (void*)&E, &_Flj_add, &_Flj_mul, &_Flj_one);
225 : }
226 :
227 : ulong
228 253556 : Flj_order_ufact(GEN P, ulong n, GEN F, ulong a4, ulong p, ulong pi)
229 : {
230 253556 : pari_sp av = avma;
231 253556 : ulong res = 1;
232 : long i, nfactors;
233 : GEN primes, exps;
234 :
235 253556 : primes = gel(F, 1);
236 253556 : nfactors = lg(primes);
237 253556 : exps = gel(F, 2);
238 :
239 721484 : for (i = 1; i < nfactors; ++i) {
240 577384 : ulong q, pp = primes[i];
241 577384 : long j, k, ei = exps[i];
242 : naf_t x;
243 : GEN b;
244 :
245 1369951 : for (q = pp, j = 1; j < ei; ++j) q *= pp;
246 577384 : b = Flj_mulu_pre(P, n / q, a4, p, pi);
247 :
248 577384 : naf_repr(&x, pp);
249 1669517 : for (j = 0; j < ei && b[3]; ++j)
250 1092133 : b = Flj_mulu_pre_naf(b, pp, a4, p, pi, &x);
251 577384 : if (b[3]) return 0;
252 1158076 : for (k = 0; k < j; ++k) res *= pp;
253 467928 : set_avma(av);
254 : }
255 144100 : return res;
256 : }
257 :
258 : GEN
259 1549125 : Fle_to_Flj(GEN P)
260 3098046 : { return ell_is_inf(P) ? mkvecsmall3(1UL, 1UL, 0UL):
261 1549062 : mkvecsmall3(P[1], P[2], 1UL);
262 : }
263 :
264 : GEN
265 61600 : Flj_to_Fle(GEN P, ulong p)
266 : {
267 61600 : if (P[3] == 0) return ellinf();
268 : else
269 : {
270 60942 : ulong Z = Fl_inv(P[3], p);
271 60942 : ulong Z2 = Fl_sqr(Z, p);
272 60942 : ulong X3 = Fl_mul(P[1], Z2, p);
273 60942 : ulong Y3 = Fl_mul(P[2], Fl_mul(Z, Z2, p), p);
274 60942 : return mkvecsmall2(X3, Y3);
275 : }
276 : }
277 :
278 : GEN
279 1563902 : Flj_to_Fle_pre(GEN P, ulong p, ulong pi)
280 : {
281 1563902 : if (P[3] == 0) return ellinf();
282 : else
283 : {
284 1382563 : ulong Z = Fl_inv(P[3], p);
285 1383034 : ulong Z2 = Fl_sqr_pre(Z, p, pi);
286 1382786 : ulong X3 = Fl_mul_pre(P[1], Z2, p, pi);
287 1382784 : ulong Y3 = Fl_mul_pre(P[2], Fl_mul_pre(Z, Z2, p, pi), p, pi);
288 1382824 : return mkvecsmall2(X3, Y3);
289 : }
290 : }
291 :
292 : INLINE void
293 6655815 : random_Fle_pre_indir(ulong a4, ulong a6, ulong p, ulong pi,
294 : ulong *pt_x, ulong *pt_y)
295 : {
296 : ulong x, x2, y, rhs;
297 : do
298 : {
299 6655815 : x = random_Fl(p); /* x^3+a4*x+a6 = x*(x^2+a4)+a6 */
300 6655898 : x2 = Fl_sqr_pre(x, p, pi);
301 6655877 : rhs = Fl_addmul_pre(a6, x, Fl_add(x2, a4, p), p, pi);
302 6655889 : } while ((!rhs && !Fl_add(Fl_triple(x2,p),a4,p)) || krouu(rhs, p) < 0);
303 3338093 : y = Fl_sqrt_pre(rhs, p, pi);
304 3338090 : *pt_x = x; *pt_y = y;
305 3338090 : }
306 :
307 : GEN
308 432888 : random_Flj_pre(ulong a4, ulong a6, ulong p, ulong pi)
309 : {
310 : ulong x, y;
311 432888 : random_Fle_pre_indir(a4, a6, p, pi, &x, &y);
312 432869 : return mkvecsmall3(x, y, 1);
313 : }
314 :
315 : GEN
316 156218 : Flj_changepointinv_pre(GEN P, GEN ch, ulong p, ulong pi)
317 : {
318 : ulong c, u, r, s, t, u2, u3, z2;
319 156218 : ulong x = uel(P,1), y = uel(P,2), z = uel(P,3);
320 : GEN w;
321 156218 : if (z == 0) return Flv_copy(P);
322 156218 : u = ch[1]; r = ch[2];
323 156218 : s = ch[3]; t = ch[4];
324 156218 : u2 = Fl_sqr_pre(u, p, pi); u3 = Fl_mul_pre(u, u2, p, pi);
325 156222 : c = Fl_mul_pre(u2, x, p, pi);
326 156209 : z2 = Fl_sqr_pre(z, p, pi);
327 156210 : w = cgetg(4, t_VECSMALL);
328 156200 : uel(w,1) = Fl_add(c, Fl_mul_pre(r, z2, p, pi), p);
329 156206 : uel(w,2) = Fl_add(Fl_mul_pre(u3 ,y, p, pi),
330 : Fl_mul_pre(z, Fl_add(Fl_mul_pre(s,c,p,pi),
331 : Fl_mul_pre(z2,t,p,pi), p), p, pi), p);
332 156214 : uel(w,3) = z;
333 156214 : return w;
334 : }
335 :
336 : /***********************************************************************/
337 : /** **/
338 : /** Fle **/
339 : /** **/
340 : /***********************************************************************/
341 : GEN
342 15807 : Fle_changepoint(GEN P, GEN ch, ulong p)
343 : {
344 : ulong c, u, r, s, t, v, v2, v3;
345 : GEN z;
346 15807 : if (ell_is_inf(P)) return ellinf();
347 15807 : u = ch[1]; r = ch[2];
348 15807 : s = ch[3]; t = ch[4];
349 15807 : v = Fl_inv(u, p); v2 = Fl_sqr(v,p); v3 = Fl_mul(v,v2,p);
350 15807 : c = Fl_sub(uel(P,1),r,p);
351 15807 : z = cgetg(3,t_VECSMALL);
352 15807 : z[1] = Fl_mul(v2, c, p);
353 15807 : z[2] = Fl_mul(v3, Fl_sub(uel(P,2), Fl_add(Fl_mul(s,c, p),t, p),p),p);
354 15807 : return z;
355 : }
356 :
357 : GEN
358 70681 : Fle_changepointinv(GEN P, GEN ch, ulong p)
359 : {
360 : ulong c, u, r, s, t, u2, u3;
361 : GEN z;
362 70681 : if (ell_is_inf(P)) return ellinf();
363 70023 : u = ch[1]; r = ch[2];
364 70023 : s = ch[3]; t = ch[4];
365 70023 : u2 = Fl_sqr(u, p); u3 = Fl_mul(u,u2,p);
366 70023 : c = Fl_mul(u2,uel(P,1), p);
367 70023 : z = cgetg(3, t_VECSMALL);
368 70023 : z[1] = Fl_add(c,r,p);
369 70023 : z[2] = Fl_add(Fl_mul(u3,uel(P,2),p), Fl_add(Fl_mul(s,c,p), t, p), p);
370 70023 : return z;
371 : }
372 : static GEN
373 572784 : Fle_dbl_slope(GEN P, ulong a4, ulong p, ulong *slope)
374 : {
375 : ulong x, y, Qx, Qy;
376 572784 : if (ell_is_inf(P) || !P[2]) return ellinf();
377 503317 : x = P[1]; y = P[2];
378 503317 : *slope = Fl_div(Fl_add(Fl_triple(Fl_sqr(x,p), p), a4, p),
379 : Fl_double(y, p), p);
380 503323 : Qx = Fl_sub(Fl_sqr(*slope, p), Fl_double(x, p), p);
381 503310 : Qy = Fl_sub(Fl_mul(*slope, Fl_sub(x, Qx, p), p), y, p);
382 503299 : return mkvecsmall2(Qx, Qy);
383 : }
384 :
385 : GEN
386 428142 : Fle_dbl(GEN P, ulong a4, ulong p)
387 : {
388 : ulong slope;
389 428142 : return Fle_dbl_slope(P,a4,p,&slope);
390 : }
391 :
392 : static GEN
393 1154967 : Fle_add_slope(GEN P, GEN Q, ulong a4, ulong p, ulong *slope)
394 : {
395 : ulong Px, Py, Qx, Qy, Rx, Ry;
396 1154967 : if (ell_is_inf(P)) return Q;
397 1154962 : if (ell_is_inf(Q)) return P;
398 1154952 : Px = P[1]; Py = P[2];
399 1154952 : Qx = Q[1]; Qy = Q[2];
400 1154952 : if (Px==Qx) return Py==Qy ? Fle_dbl_slope(P, a4, p, slope): ellinf();
401 1033836 : *slope = Fl_div(Fl_sub(Py, Qy, p), Fl_sub(Px, Qx, p), p);
402 1034027 : Rx = Fl_sub(Fl_sub(Fl_sqr(*slope, p), Px, p), Qx, p);
403 1033975 : Ry = Fl_sub(Fl_mul(*slope, Fl_sub(Px, Rx, p), p), Py, p);
404 1033963 : return mkvecsmall2(Rx, Ry);
405 : }
406 :
407 : GEN
408 1150016 : Fle_add(GEN P, GEN Q, ulong a4, ulong p)
409 : {
410 : ulong slope;
411 1150016 : return Fle_add_slope(P,Q,a4,p,&slope);
412 : }
413 :
414 : static GEN
415 125846 : Fle_neg(GEN P, ulong p)
416 : {
417 125846 : if (ell_is_inf(P)) return P;
418 125846 : return mkvecsmall2(P[1], Fl_neg(P[2], p));
419 : }
420 :
421 : GEN
422 0 : Fle_sub(GEN P, GEN Q, ulong a4, ulong p)
423 : {
424 0 : pari_sp av = avma;
425 : ulong slope;
426 0 : return gerepileupto(av, Fle_add_slope(P, Fle_neg(Q, p), a4, p, &slope));
427 : }
428 :
429 : struct _Fle { ulong a4, a6, p; };
430 :
431 : static GEN
432 0 : _Fle_dbl(void *E, GEN P)
433 : {
434 0 : struct _Fle *ell = (struct _Fle *) E;
435 0 : return Fle_dbl(P, ell->a4, ell->p);
436 : }
437 :
438 : static GEN
439 364490 : _Fle_add(void *E, GEN P, GEN Q)
440 : {
441 364490 : struct _Fle *ell=(struct _Fle *) E;
442 364490 : return Fle_add(P, Q, ell->a4, ell->p);
443 : }
444 :
445 : GEN
446 1981595 : Fle_mulu(GEN P, ulong n, ulong a4, ulong p)
447 : {
448 : ulong pi;
449 1981595 : if (!n || ell_is_inf(P)) return ellinf();
450 1981518 : if (n==1) return zv_copy(P);
451 1969912 : if (n==2) return Fle_dbl(P, a4, p);
452 1548878 : pi = get_Fl_red(p);
453 1549082 : return Flj_to_Fle_pre(Flj_mulu_pre(Fle_to_Flj(P), n, a4, p, pi), p, pi);
454 : }
455 :
456 : static GEN
457 1454622 : _Fle_mul(void *E, GEN P, GEN n)
458 : {
459 1454622 : pari_sp av = avma;
460 1454622 : struct _Fle *e=(struct _Fle *) E;
461 1454622 : long s = signe(n);
462 : GEN Q;
463 1454622 : if (!s || ell_is_inf(P)) return ellinf();
464 1442854 : if (s < 0) P = Fle_neg(P, e->p);
465 1442854 : if (is_pm1(n)) return s > 0? zv_copy(P): P;
466 1138573 : Q = (lgefint(n)==3) ? Fle_mulu(P, uel(n,2), e->a4, e->p):
467 6 : gen_pow(P, n, (void*)e, &_Fle_dbl, &_Fle_add);
468 1138508 : return s > 0? Q: gerepileuptoleaf(av, Q);
469 : }
470 :
471 : GEN
472 238 : Fle_mul(GEN P, GEN n, ulong a4, ulong p)
473 : {
474 : struct _Fle E;
475 238 : E.a4 = a4; E.p = p;
476 238 : return _Fle_mul(&E, P, n);
477 : }
478 :
479 : /* Finds a random nonsingular point on E */
480 : GEN
481 2905221 : random_Fle_pre(ulong a4, ulong a6, ulong p, ulong pi)
482 : {
483 : ulong x, y;
484 2905221 : random_Fle_pre_indir(a4, a6, p, pi, &x, &y);
485 2905221 : return mkvecsmall2(x, y);
486 : }
487 :
488 : GEN
489 22456 : random_Fle(ulong a4, ulong a6, ulong p)
490 22456 : { return random_Fle_pre(a4, a6, p, get_Fl_red(p)); }
491 :
492 : static GEN
493 0 : _Fle_rand(void *E)
494 : {
495 0 : struct _Fle *e=(struct _Fle *) E;
496 0 : return random_Fle(e->a4, e->a6, e->p);
497 : }
498 :
499 : static const struct bb_group Fle_group={_Fle_add,_Fle_mul,_Fle_rand,hash_GEN,zv_equal,ell_is_inf,NULL};
500 :
501 : GEN
502 284380 : Fle_order(GEN z, GEN o, ulong a4, ulong p)
503 : {
504 284380 : pari_sp av = avma;
505 : struct _Fle e;
506 284380 : e.a4=a4;
507 284380 : e.p=p;
508 284380 : return gerepileuptoint(av, gen_order(z, o, (void*)&e, &Fle_group));
509 : }
510 :
511 : GEN
512 50477 : Fle_log(GEN a, GEN b, GEN o, ulong a4, ulong p)
513 : {
514 50477 : pari_sp av = avma;
515 : struct _Fle e;
516 50477 : e.a4=a4;
517 50477 : e.p=p;
518 50477 : return gerepileuptoint(av, gen_PH_log(a, b, o, (void*)&e, &Fle_group));
519 : }
520 :
521 : ulong
522 0 : Fl_ellj(ulong a4, ulong a6, ulong p)
523 : {
524 0 : if (SMALL_ULONG(p))
525 : { /* a43 = 4 a4^3 */
526 0 : ulong a43 = Fl_double(Fl_double(Fl_mul(a4, Fl_sqr(a4, p), p), p), p);
527 : /* a62 = 27 a6^2 */
528 0 : ulong a62 = Fl_mul(Fl_sqr(a6, p), 27 % p, p);
529 0 : ulong z1 = Fl_mul(a43, 1728 % p, p);
530 0 : ulong z2 = Fl_add(a43, a62, p);
531 0 : return Fl_div(z1, z2, p);
532 : }
533 0 : return Fl_ellj_pre(a4, a6, p, get_Fl_red(p));
534 : }
535 :
536 : void
537 132840 : Fl_ellj_to_a4a6(ulong j, ulong p, ulong *pt_a4, ulong *pt_a6)
538 : {
539 132840 : ulong zagier = 1728 % p;
540 132840 : if (j == 0) { *pt_a4 = 0; *pt_a6 =1; }
541 132840 : else if (j == zagier) { *pt_a4 = 1; *pt_a6 =0; }
542 : else
543 : {
544 132840 : ulong k = Fl_sub(zagier, j, p);
545 132840 : ulong kj = Fl_mul(k, j, p);
546 132841 : ulong k2j = Fl_mul(kj, k, p);
547 132842 : *pt_a4 = Fl_triple(kj, p);
548 132840 : *pt_a6 = Fl_double(k2j, p);
549 : }
550 132841 : }
551 :
552 : ulong
553 1706957 : Fl_elldisc_pre(ulong a4, ulong a6, ulong p, ulong pi)
554 : { /* D = -(4A^3 + 27B^2) */
555 : ulong t1, t2;
556 1706957 : t1 = Fl_mul_pre(a4, Fl_sqr_pre(a4, p, pi), p, pi);
557 1706957 : t1 = Fl_double(Fl_double(t1, p), p);
558 1706957 : t2 = Fl_mul_pre(27 % p, Fl_sqr_pre(a6, p, pi), p, pi);
559 1706957 : return Fl_neg(Fl_add(t1, t2, p), p);
560 : }
561 :
562 : ulong
563 0 : Fl_elldisc(ulong a4, ulong a6, ulong p)
564 : {
565 0 : if (SMALL_ULONG(p))
566 : { /* D = -(4A^3 + 27B^2) */
567 : ulong t1, t2;
568 0 : t1 = Fl_mul(a4, Fl_sqr(a4, p), p);
569 0 : t1 = Fl_double(Fl_double(t1, p), p);
570 0 : t2 = Fl_mul(27 % p, Fl_sqr(a6, p), p);
571 0 : return Fl_neg(Fl_add(t1, t2, p), p);
572 : }
573 0 : return Fl_elldisc_pre(a4, a6, p, get_Fl_red(p));
574 : }
575 :
576 : void
577 93547 : Fl_elltwist_disc(ulong a4, ulong a6, ulong D, ulong p, ulong *pa4, ulong *pa6)
578 : {
579 93547 : ulong D2 = Fl_sqr(D, p);
580 93547 : *pa4 = Fl_mul(a4, D2, p);
581 93547 : *pa6 = Fl_mul(a6, Fl_mul(D, D2, p), p);
582 93548 : }
583 :
584 : void
585 0 : Fl_elltwist(ulong a4, ulong a6, ulong p, ulong *pt_a4, ulong *pt_a6)
586 0 : { Fl_elltwist_disc(a4, a6, nonsquare_Fl(p), p, pt_a4, pt_a6); }
587 :
588 : static void
589 42450151 : Fle_dbl_sinv_pre_inplace(GEN P, ulong a4, ulong sinv, ulong p, ulong pi)
590 : {
591 : ulong x, y, slope;
592 42450151 : if (uel(P,1)==p) return;
593 42237169 : if (!P[2]) { P[1] = p; return; }
594 42108740 : x = P[1]; y = P[2];
595 42108740 : slope = Fl_mul_pre(Fl_add(Fl_triple(Fl_sqr_pre(x, p, pi), p), a4, p),
596 : sinv, p, pi);
597 42108740 : P[1] = Fl_sub(Fl_sqr_pre(slope, p, pi), Fl_double(x, p), p);
598 42108740 : P[2] = Fl_sub(Fl_mul_pre(slope, Fl_sub(x, P[1], p), p, pi), y, p);
599 : }
600 :
601 : static void
602 5623762 : Fle_add_sinv_pre_inplace(GEN P, GEN Q, ulong a4, ulong sinv, ulong p, ulong pi)
603 : {
604 : ulong Px, Py, Qx, Qy, slope;
605 5623762 : if (uel(P,1)==p) { P[1] = Q[1]; P[2] = Q[2]; }
606 5623762 : if (ell_is_inf(Q)) return;
607 5623762 : Px = P[1]; Py = P[2];
608 5623762 : Qx = Q[1]; Qy = Q[2];
609 5623762 : if (Px==Qx)
610 : {
611 24515 : if (Py==Qy) Fle_dbl_sinv_pre_inplace(P, a4, sinv, p, pi);
612 11338 : else P[1] = p;
613 24515 : return;
614 : }
615 5599247 : slope = Fl_mul_pre(Fl_sub(Py, Qy, p), sinv, p, pi);
616 5599247 : P[1] = Fl_sub(Fl_sub(Fl_sqr_pre(slope, p, pi), Px, p), Qx, p);
617 5599247 : P[2] = Fl_sub(Fl_mul_pre(slope, Fl_sub(Px, P[1], p), p, pi), Py, p);
618 : }
619 :
620 : static void
621 6278391 : Fle_sub_sinv_pre_inplace(GEN P, GEN Q, ulong a4, ulong sinv, ulong p, ulong pi)
622 : {
623 : ulong Px, Py, Qx, Qy, slope;
624 6278391 : if (uel(P,1)==p) { P[1] = Q[1]; P[2] = Fl_neg(Q[2], p); }
625 6278391 : if (ell_is_inf(Q)) return;
626 6278391 : Px = P[1]; Py = P[2];
627 6278391 : Qx = Q[1]; Qy = Q[2];
628 6278391 : if (Px==Qx)
629 : {
630 28823 : if (Py==Qy) P[1] = p;
631 : else
632 10991 : Fle_dbl_sinv_pre_inplace(P, a4, sinv, p, pi);
633 28823 : return;
634 : }
635 6249568 : slope = Fl_mul_pre(Fl_add(Py, Qy, p), sinv, p, pi);
636 6249568 : P[1] = Fl_sub(Fl_sub(Fl_sqr_pre(slope, p, pi), Px, p), Qx, p);
637 6249568 : P[2] = Fl_sub(Fl_mul_pre(slope, Fl_sub(Px, P[1], p), p, pi), Py, p);
638 : }
639 :
640 : static long
641 54091689 : skipzero(long n) { return n ? n:1; }
642 :
643 : void
644 1244832 : FleV_add_pre_inplace(GEN P, GEN Q, GEN a4, ulong p, ulong pi)
645 : {
646 1244832 : long i, l=lg(a4);
647 1244832 : GEN sinv = cgetg(l, t_VECSMALL);
648 6868594 : for(i=1; i<l; i++)
649 5623762 : uel(sinv,i) = umael(P,i,1)==p? 1: skipzero(Fl_sub(mael(P,i,1), mael(Q,i,1), p));
650 1244832 : Flv_inv_pre_inplace(sinv, p, pi);
651 6868594 : for (i=1; i<l; i++)
652 5623762 : Fle_add_sinv_pre_inplace(gel(P,i), gel(Q,i), uel(a4,i), uel(sinv,i), p, pi);
653 1244832 : }
654 :
655 : void
656 1415250 : FleV_sub_pre_inplace(GEN P, GEN Q, GEN a4, ulong p, ulong pi)
657 : {
658 1415250 : long i, l=lg(a4);
659 1415250 : GEN sinv = cgetg(l, t_VECSMALL);
660 7693641 : for(i=1; i<l; i++)
661 6278391 : uel(sinv,i) = umael(P,i,1)==p? 1: skipzero(Fl_sub(mael(P,i,1), mael(Q,i,1), p));
662 1415250 : Flv_inv_pre_inplace(sinv, p, pi);
663 7693641 : for (i=1; i<l; i++)
664 6278391 : Fle_sub_sinv_pre_inplace(gel(P,i), gel(Q,i), uel(a4,i), uel(sinv,i), p, pi);
665 1415250 : }
666 :
667 : void
668 9612189 : FleV_dbl_pre_inplace(GEN P, GEN a4, ulong p, ulong pi)
669 : {
670 9612189 : long i, l=lg(a4);
671 9612189 : GEN sinv = cgetg(l, t_VECSMALL);
672 52038172 : for(i=1; i<l; i++)
673 42425983 : uel(sinv,i) = umael(P,i,1)==p? 1: skipzero(Fl_double(umael(P,i,2), p));
674 9612189 : Flv_inv_pre_inplace(sinv, p, pi);
675 52038172 : for(i=1; i<l; i++)
676 42425983 : Fle_dbl_sinv_pre_inplace(gel(P,i), uel(a4,i), uel(sinv,i), p, pi);
677 9612189 : }
678 :
679 : static void
680 770561 : FleV_mulu_pre_naf_inplace(GEN P, ulong n, GEN a4, ulong p, ulong pi, const naf_t *x)
681 : {
682 770561 : pari_sp av = avma;
683 : ulong pbits, nbits, m;
684 : GEN R;
685 770561 : if (n == 1) return;
686 :
687 770561 : R = P; P = gcopy(P);
688 770561 : FleV_dbl_pre_inplace(R, a4, p, pi);
689 770561 : if (n == 2) return;
690 :
691 770505 : pbits = x->pbits;
692 770505 : nbits = x->nbits;
693 770505 : m = (1UL << x->lnzb);
694 9612133 : for ( ; m; m >>= 1) {
695 8841628 : FleV_dbl_pre_inplace(R, a4, p, pi);
696 8841628 : if (m & pbits)
697 1244832 : FleV_add_pre_inplace(R, P, a4, p, pi);
698 7596796 : else if (m & nbits)
699 1415250 : FleV_sub_pre_inplace(R, P, a4, p, pi);
700 : }
701 770505 : set_avma(av);
702 : }
703 :
704 : void
705 770561 : FleV_mulu_pre_inplace(GEN P, ulong n, GEN a4, ulong p, ulong pi)
706 : {
707 770561 : naf_t x; naf_repr(&x, n);
708 770561 : FleV_mulu_pre_naf_inplace(P, n, a4, p, pi, &x);
709 770561 : }
710 :
711 : /***********************************************************************/
712 : /** **/
713 : /** Pairings **/
714 : /** **/
715 : /***********************************************************************/
716 :
717 : /* Derived from APIP from and by Jerome Milan, 2012 */
718 :
719 : static ulong
720 98099 : Fle_vert(GEN P, GEN Q, ulong a4, ulong p)
721 : {
722 98099 : if (ell_is_inf(P))
723 33768 : return 1;
724 64331 : if (uel(Q, 1) != uel(P, 1))
725 60880 : return Fl_sub(uel(Q, 1), uel(P, 1), p);
726 3451 : if (uel(P,2) != 0 ) return 1;
727 2382 : return Fl_inv(Fl_add(Fl_triple(Fl_sqr(uel(P,1),p), p), a4, p), p);
728 : }
729 :
730 : static ulong
731 34232 : Fle_Miller_line(GEN R, GEN Q, ulong slope, ulong a4, ulong p)
732 : {
733 34232 : ulong x = uel(Q, 1), y = uel(Q, 2);
734 34232 : ulong tmp1 = Fl_sub(x, uel(R, 1), p);
735 34232 : ulong tmp2 = Fl_add(Fl_mul(tmp1, slope, p), uel(R,2), p);
736 34232 : if (y != tmp2)
737 32672 : return Fl_sub(y, tmp2, p);
738 1560 : if (y == 0)
739 1009 : return 1;
740 : else
741 : {
742 : ulong s1, s2;
743 551 : ulong y2i = Fl_inv(Fl_double(y, p), p);
744 551 : s1 = Fl_mul(Fl_add(Fl_triple(Fl_sqr(x, p), p), a4, p), y2i, p);
745 551 : if (s1 != slope)
746 243 : return Fl_sub(s1, slope, p);
747 308 : s2 = Fl_mul(Fl_sub(Fl_triple(x, p), Fl_sqr(s1, p), p), y2i, p);
748 308 : return s2 ? s2: y2i;
749 : }
750 : }
751 :
752 : /* Computes the equation of the line tangent to R and returns its
753 : evaluation at the point Q. Also doubles the point R.
754 : */
755 :
756 : static ulong
757 53834 : Fle_tangent_update(GEN R, GEN Q, ulong a4, ulong p, GEN *pt_R)
758 : {
759 53834 : if (ell_is_inf(R))
760 : {
761 4145 : *pt_R = ellinf();
762 4145 : return 1;
763 : }
764 49689 : else if (uel(R,2) == 0)
765 : {
766 20397 : *pt_R = ellinf();
767 20397 : return Fle_vert(R, Q, a4, p);
768 : } else {
769 : ulong slope;
770 29292 : *pt_R = Fle_dbl_slope(R, a4, p, &slope);
771 29292 : return Fle_Miller_line(R, Q, slope, a4, p);
772 : }
773 : }
774 :
775 : /* Computes the equation of the line through R and P, and returns its
776 : evaluation at the point Q. Also adds P to the point R.
777 : */
778 :
779 : static ulong
780 14404 : Fle_chord_update(GEN R, GEN P, GEN Q, ulong a4, ulong p, GEN *pt_R)
781 : {
782 14404 : if (ell_is_inf(R))
783 : {
784 238 : *pt_R = gcopy(P);
785 238 : return Fle_vert(P, Q, a4, p);
786 : }
787 14166 : else if (ell_is_inf(P))
788 : {
789 0 : *pt_R = gcopy(R);
790 0 : return Fle_vert(R, Q, a4, p);
791 : }
792 14166 : else if (uel(P, 1) == uel(R, 1))
793 : {
794 9226 : if (uel(P, 2) == uel(R, 2))
795 0 : return Fle_tangent_update(R, Q, a4, p, pt_R);
796 : else {
797 9226 : *pt_R = ellinf();
798 9226 : return Fle_vert(R, Q, a4, p);
799 : }
800 : } else {
801 : ulong slope;
802 4940 : *pt_R = Fle_add_slope(P, R, a4, p, &slope);
803 4940 : return Fle_Miller_line(R, Q, slope, a4, p);
804 : }
805 : }
806 :
807 : struct _Fle_miller { ulong p, a4; GEN P; };
808 : static GEN
809 53834 : Fle_Miller_dbl(void* E, GEN d)
810 : {
811 53834 : struct _Fle_miller *m = (struct _Fle_miller *)E;
812 53834 : ulong p = m->p, a4 = m->a4;
813 53834 : GEN P = m->P;
814 : ulong v, line;
815 53834 : ulong N = Fl_sqr(umael(d,1,1), p);
816 53834 : ulong D = Fl_sqr(umael(d,1,2), p);
817 53834 : GEN point = gel(d,2);
818 53834 : line = Fle_tangent_update(point, P, a4, p, &point);
819 53834 : N = Fl_mul(N, line, p);
820 53834 : v = Fle_vert(point, P, a4, p);
821 53834 : D = Fl_mul(D, v, p); return mkvec2(mkvecsmall2(N, D), point);
822 : }
823 : static GEN
824 14404 : Fle_Miller_add(void* E, GEN va, GEN vb)
825 : {
826 14404 : struct _Fle_miller *m = (struct _Fle_miller *)E;
827 14404 : ulong p = m->p, a4= m->a4;
828 14404 : GEN P = m->P, point;
829 : ulong v, line;
830 14404 : ulong na = umael(va,1,1), da = umael(va,1,2);
831 14404 : ulong nb = umael(vb,1,1), db = umael(vb,1,2);
832 14404 : GEN pa = gel(va,2), pb = gel(vb,2);
833 14404 : ulong N = Fl_mul(na, nb, p);
834 14404 : ulong D = Fl_mul(da, db, p);
835 14404 : line = Fle_chord_update(pa, pb, P, a4, p, &point);
836 14404 : N = Fl_mul(N, line, p);
837 14404 : v = Fle_vert(point, P, a4, p);
838 14404 : D = Fl_mul(D, v, p); return mkvec2(mkvecsmall2(N, D), point);
839 : }
840 :
841 : /* Returns the Miller function f_{m, Q} evaluated at the point P using
842 : * the standard Miller algorithm. */
843 : static ulong
844 29385 : Fle_Miller(GEN Q, GEN P, ulong m, ulong a4, ulong p)
845 : {
846 29385 : pari_sp av = avma;
847 : struct _Fle_miller d;
848 : GEN v;
849 : ulong N, D;
850 :
851 29385 : d.a4 = a4; d.p = p; d.P = P;
852 29385 : v = gen_powu_i(mkvec2(mkvecsmall2(1,1), Q), m, (void*)&d,
853 : Fle_Miller_dbl, Fle_Miller_add);
854 29385 : N = umael(v,1,1); D = umael(v,1,2);
855 29385 : return gc_ulong(av, Fl_div(N, D, p));
856 : }
857 :
858 : ulong
859 8865 : Fle_weilpairing(GEN P, GEN Q, ulong m, ulong a4, ulong p)
860 : {
861 8865 : pari_sp ltop = avma;
862 : ulong N, D, w;
863 8865 : if (ell_is_inf(P) || ell_is_inf(Q) || zv_equal(P,Q)) return 1;
864 8158 : N = Fle_Miller(P, Q, m, a4, p);
865 8158 : D = Fle_Miller(Q, P, m, a4, p);
866 8158 : w = Fl_div(N, D, p);
867 8158 : if (odd(m)) w = Fl_neg(w, p);
868 8158 : return gc_ulong(ltop, w);
869 : }
870 :
871 : ulong
872 13069 : Fle_tatepairing(GEN P, GEN Q, ulong m, ulong a4, ulong p)
873 : {
874 13069 : if (ell_is_inf(P) || ell_is_inf(Q)) return 1;
875 13069 : return Fle_Miller(P, Q, m, a4, p);
876 : }
877 :
878 : GEN
879 10178 : Fl_ellptors(ulong l, ulong N, ulong a4, ulong a6, ulong p)
880 : {
881 10178 : long v = z_lval(N, l);
882 : ulong N0, N1;
883 : GEN F;
884 10178 : if (v==0) return cgetg(1,t_VEC);
885 10178 : N0 = upowuu(l, v); N1 = N/N0;
886 10178 : F = mkmat2(mkcols(l), mkcols(v));
887 : while(1)
888 1050 : {
889 11228 : pari_sp av2 = avma;
890 : GEN P, Q;
891 : ulong d, s, t;
892 :
893 11228 : P = Fle_mulu(random_Fle(a4, a6, p), N1, a4, p);
894 11228 : Q = Fle_mulu(random_Fle(a4, a6, p), N1, a4, p);
895 11228 : s = itou(Fle_order(P, F, a4, p));
896 11228 : t = itou(Fle_order(Q, F, a4, p));
897 11228 : if (s < t) { swap(P,Q); lswap(s,t); }
898 11228 : if (s == N0) retmkvec(Fle_mulu(P, s/l, a4, p));
899 1645 : d = Fl_order(Fle_weilpairing(P, Q, s, a4, p), s, p);
900 1645 : if (s*d == N0)
901 595 : retmkvec2(Fle_mulu(P, s/l, a4, p), Fle_mulu(Q, t/l, a4, p));
902 1050 : set_avma(av2);
903 : }
904 : }
|
