Line data Source code
1 : /* addlaws.c - various addition laws for ECM
2 : Author: F. Morain
3 : */
4 :
5 : #include <stdio.h>
6 : #include <stdlib.h>
7 : #include <string.h>
8 : #include <math.h>
9 :
10 : #include <gmp.h> /* GMP header file */
11 :
12 : #include "ecm.h" /* ecm header file */
13 : #include "ecm-impl.h"
14 : #include "ecm-ecm.h"
15 : #include "mpmod.h"
16 :
17 : #include "addlaws.h"
18 :
19 : #if DEBUG_ADD_LAWS >= 1
20 : void
21 : print_mpz_from_mpres(mpres_t x, mpmod_t n)
22 : {
23 : mpz_t tmp;
24 :
25 : mpz_init(tmp);
26 : mpres_get_z(tmp, x, n);
27 : gmp_printf("%Zd", tmp);
28 : mpz_clear(tmp);
29 : }
30 : #endif
31 :
32 : /******************** Weierstrass section ********************/
33 :
34 : void
35 50 : pt_w_set_to_zero(ell_point_t P, mpmod_t n)
36 : {
37 50 : mpres_set_ui(P->x, 0, n);
38 50 : mpres_set_ui(P->y, 1, n);
39 50 : mpres_set_ui(P->z, 0, n);
40 50 : }
41 :
42 : int
43 7436860 : pt_w_is_zero(mpres_t z, mpmod_t n)
44 : {
45 7436860 : return mpres_is_zero(z, n);
46 : }
47 :
48 : void
49 470 : pt_w_set(mpres_t x0, mpres_t y0, mpres_t z0,
50 : mpres_t x, mpres_t y, mpres_t z,
51 : ATTRIBUTE_UNUSED mpmod_t n)
52 : {
53 470 : mpres_set(x0, x, n);
54 470 : mpres_set(y0, y, n);
55 470 : mpres_set(z0, z, n);
56 470 : }
57 :
58 : #if DEBUG_ADD_LAWS >= 1
59 : void
60 : pt_w_print(mpres_t x, mpres_t y, mpres_t z, ell_curve_t E, mpmod_t n)
61 : {
62 : printf("[");
63 : print_mpz_from_mpres(x, n);
64 : printf(", ");
65 : print_mpz_from_mpres(y, n);
66 : printf(", ");
67 : if(E->type == ECM_EC_TYPE_WEIERSTRASS && E->law == ECM_LAW_AFFINE)
68 : gmp_printf("%Zd", z);
69 : else
70 : print_mpz_from_mpres(z, n);
71 : printf("]");
72 : }
73 : #endif
74 :
75 : /* [x0, y0, z0] <- [x1, y1, z1] + [x2, y2, z2] using lambda=num/den
76 : with buffer inv.
77 :
78 : (lambda*x+mu)^2+a1*x*(lambda*x+mu)+a3*(lambda*x+mu)=x^3+a2*x^2+...
79 : x^3+(a2-lambda^2-a1*lambda)*x^2+... = 0
80 : x1+x2+x3 = lambda^2+a1*lambda-a2.
81 : y3 = lambda*(x1-x3)-y1-a1*x3-a3
82 : */
83 : static int
84 2449090 : pt_w_common_aff(mpz_t f, mpres_t x0, mpres_t y0, mpres_t z0,
85 : mpres_t x1, mpres_t y1,
86 : mpres_t x2, mpres_t a1, mpres_t a3, mpres_t a2,
87 : mpmod_t n, mpres_t num, mpres_t den, mpres_t lambda)
88 : {
89 2449090 : if(mpres_invert(lambda, den, n) == 0){
90 130 : mpres_gcd(f, den, n);
91 130 : return 0;
92 : }
93 : /** lambda = num/den **/
94 2448960 : mpres_mul(lambda, lambda, num, n);
95 : /** num <- (lambda+a1)*lambda **/
96 2448960 : mpres_add(num, lambda, a1, n);
97 2448960 : mpres_mul(num, num, lambda, n);
98 2448960 : mpres_sub(num, num, a2, n);
99 : /** x0 = den <- num-x1-x2 **/
100 2448960 : mpres_sub(den, num, x1, n);
101 2448960 : mpres_sub(den, den, x2, n);
102 : /** y0 = num <- lambda*(x1-x0)-(y1+a1*x0+a3) **/
103 2448960 : mpres_sub(num, x1, den, n);
104 2448960 : mpres_mul(num, num, lambda, n);
105 2448960 : mpres_sub(y0, num, y1, n);
106 2448960 : mpres_sub(y0, y0, a3, n);
107 2448960 : mpres_mul(x0, a1, den, n);
108 2448960 : mpres_sub(y0, y0, x0, n);
109 : /** finish **/
110 2448960 : mpres_set(x0, den, n);
111 2448960 : mpz_set_ui(z0, 1); /* just in case */
112 2448960 : return 1;
113 : }
114 :
115 : /* [x3, y3, z3] <- [2] * [x1, y1, z1] */
116 : int
117 3390100 : pt_w_duplicate(mpz_t f, mpres_t x3, mpres_t y3, mpres_t z3,
118 : mpres_t x1, mpres_t y1, mpres_t z1,
119 : mpmod_t n, ell_curve_t E)
120 : {
121 3390100 : if(pt_w_is_zero(z1, n) == 1){
122 450 : pt_w_set(x3, y3, z3, x1, y1, z1, n);
123 450 : return 1;
124 : }
125 3389650 : if(E->type == ECM_EC_TYPE_WEIERSTRASS && E->law == ECM_LAW_AFFINE){
126 : /* buf[1] <- 2*y1+a1*x1+a3 */
127 1608890 : mpres_mul(E->buf[1], E->a1, x1, n);
128 1608890 : mpres_add(E->buf[1], E->buf[1], E->a3, n);
129 1608890 : mpres_add(E->buf[1], E->buf[1], y1, n);
130 1608890 : mpres_add(E->buf[1], E->buf[1], y1, n);
131 1608890 : if(mpres_is_zero(E->buf[1], n)){
132 : /* buf1 = 0 <=> P is a [2]-torsion point */
133 0 : mpres_set_ui(x3, 0, n);
134 0 : mpres_set_ui(y3, 1, n);
135 0 : mpres_set_ui(z3, 0, n);
136 0 : return 1;
137 : }
138 : /* buf[0] <- 3*x^2+2*a2*x+a4-a1*y = (3*x+2*a2)*x+a4-a1*y */
139 1608890 : mpres_mul_ui(E->buf[0], x1, 3, n);
140 1608890 : mpres_add(E->buf[0], E->buf[0], E->a2, n);
141 1608890 : mpres_add(E->buf[0], E->buf[0], E->a2, n);
142 1608890 : mpres_mul(E->buf[0], E->buf[0], x1, n);
143 1608890 : mpres_add(E->buf[0], E->buf[0], E->a4, n);
144 1608890 : mpres_mul(E->buf[2], E->a1, y1, n);
145 1608890 : mpres_sub(E->buf[0], E->buf[0], E->buf[2], n);
146 1608890 : return pt_w_common_aff(f, x3, y3, z3, x1, y1, x1,
147 1608890 : E->a1, E->a3, E->a2, n,
148 1608890 : E->buf[0], E->buf[1], E->buf[2]);
149 : }
150 1780760 : else if(E->type == ECM_EC_TYPE_WEIERSTRASS
151 1780760 : && E->law == ECM_LAW_HOMOGENEOUS){
152 : /* source is dbl-2007-bl: 5M + 6S + 1*a + 7add + 3*2 + 1*3 */
153 : /* mapping: h = buf[0], w = buf[1], s = buf[2], RR = buf[3], B = buf[4];*/
154 : /* h:=X1^2 mod p; # S*/
155 1780760 : mpres_sqr(E->buf[0], x1, n);
156 : /* w:=Z1^2 mod p;*/
157 1780760 : mpres_sqr(E->buf[1], z1, n);
158 : /* w:=a*w mod p;*/
159 1780760 : mpres_mul(E->buf[1], E->buf[1], E->a4, n);
160 : /* s:=3*h mod p; # *3*/
161 1780760 : mpres_mul_ui(E->buf[2], E->buf[0], 3, n);
162 : /* w:=w+s mod p;*/
163 1780760 : mpres_add(E->buf[1], E->buf[1], E->buf[2], n);
164 : /* s:=Y1*Z1 mod p;*/
165 1780760 : mpres_mul(E->buf[2], y1, z1, n);
166 : /* s:=2*s mod p;*/
167 1780760 : mpres_mul_ui(E->buf[2], E->buf[2], 2, n);
168 : /* Z3:=s^2 mod p;*/
169 1780760 : mpres_sqr(z3, E->buf[2], n);
170 : /* Z3:=s*Z3 mod p;*/
171 1780760 : mpres_mul(z3, z3, E->buf[2], n);
172 : /* RR:=Y1*s mod p; # M*/
173 1780760 : mpres_mul(E->buf[3], y1, E->buf[2], n);
174 : /* B:=X1+RR mod p; # add*/
175 1780760 : mpres_add(E->buf[4], x1, E->buf[3], n);
176 : /* B:=B^2 mod p;*/
177 1780760 : mpres_sqr(E->buf[4], E->buf[4], n);
178 : /* RR:=RR^2 mod p; # S*/
179 1780760 : mpres_sqr(E->buf[3], E->buf[3], n);
180 : /* B:=B-h mod p;*/
181 1780760 : mpres_sub(E->buf[4], E->buf[4], E->buf[0], n);
182 : /* B:=B-RR mod p;*/
183 1780760 : mpres_sub(E->buf[4], E->buf[4], E->buf[3], n);
184 : /* h:=w^2 mod p;*/
185 1780760 : mpres_sqr(E->buf[0], E->buf[1], n);
186 : /* X3:=2*B mod p;*/
187 1780760 : mpres_mul_ui(x3, E->buf[4], 2, n);
188 : /* h:=h-X3 mod p;*/
189 1780760 : mpres_sub(E->buf[0], E->buf[0], x3, n);
190 : /* X3:=h*s mod p; # M*/
191 1780760 : mpres_mul(x3, E->buf[0], E->buf[2], n);
192 : /* s:=B-h mod p;*/
193 1780760 : mpres_sub(E->buf[2], E->buf[4], E->buf[0], n);
194 : /* s:=w*s mod p;*/
195 1780760 : mpres_mul(E->buf[2], E->buf[2], E->buf[1], n);
196 : /* Y3:=2*RR mod p;*/
197 1780760 : mpres_mul_ui(y3, E->buf[3], 2, n);
198 : /* Y3:=s-Y3 mod p;*/
199 1780760 : mpres_sub(y3, E->buf[2], y3, n);
200 1780760 : return 1;
201 : }
202 0 : return 0;
203 : }
204 :
205 : /* [x3, y3, z3] <- [x1, y1, z1] + [x2, y2, z2]; P3 can be either P1 or P2. */
206 : int
207 1772920 : pt_w_add(mpz_t f, mpres_t x3, mpres_t y3, mpres_t z3,
208 : mpres_t x1, mpres_t y1, mpres_t z1,
209 : mpres_t x2, mpres_t y2, mpres_t z2,
210 : mpmod_t n, ell_curve_t E)
211 : {
212 1772920 : if(pt_w_is_zero(z1, n)){
213 0 : pt_w_set(x3, y3, z3, x2, y2, z2, n);
214 0 : return 1;
215 : }
216 1772920 : else if(pt_w_is_zero(z2, n)){
217 20 : pt_w_set(x3, y3, z3, x1, y1, z1, n);
218 20 : return 1;
219 : }
220 1772900 : if(E->type == ECM_EC_TYPE_WEIERSTRASS && E->law == ECM_LAW_AFFINE)
221 840220 : if(mpres_equal(x1, x2, n) && mpres_equal(y1, y2, n))
222 20 : return pt_w_duplicate(f, x3, y3, z3, x1, y1, z1, n, E);
223 : else{
224 840200 : mpres_sub(E->buf[0], y1, y2, n);
225 840200 : mpres_sub(E->buf[1], x1, x2, n);
226 840200 : return pt_w_common_aff(f, x3, y3, z3, x1, y1, x2,
227 840200 : E->a1, E->a3, E->a2,
228 840200 : n, E->buf[0], E->buf[1], E->buf[2]);
229 : }
230 932680 : else if(E->type == ECM_EC_TYPE_WEIERSTRASS
231 932680 : && E->law == ECM_LAW_HOMOGENEOUS){
232 : /* Cohen-Miyaji-Ono: 12M+2S+6add+1*2 */
233 : /* mapping: y1z2 = buf, AA = buf+1, u = buf+2, v = buf+3, R = buf+4, */
234 : /* vvv = buf+5; */
235 : #if DEBUG_ADD_LAWS >= 2
236 : printf("y1="); print_mpz_from_mpres(y1, n); printf("\n");
237 : printf("y2="); print_mpz_from_mpres(y2, n); printf("\n");
238 : printf("z1="); print_mpz_from_mpres(z1, n); printf("\n");
239 : printf("z2="); print_mpz_from_mpres(z2, n); printf("\n");
240 : #endif
241 : /* Y1Z2:=Y1*Z2 mod p; # M*/
242 932680 : mpres_mul(E->buf[0], y1, z2, n);
243 : /* A:=X1*Z2 mod p; # M*/
244 932680 : mpres_mul(E->buf[1], x1, z2, n);
245 : /* u:=Y2*Z1 mod p;*/
246 932680 : mpres_mul(E->buf[2], y2, z1, n);
247 : /* u:=u-Y1Z2 mod p;*/
248 932680 : mpres_sub(E->buf[2], E->buf[2], E->buf[0], n);
249 : /* v:=X2*Z1 mod p;*/
250 932680 : mpres_mul(E->buf[3], x2, z1, n);
251 : /* v:=v-A mod p;*/
252 932680 : mpres_sub(E->buf[3], E->buf[3], E->buf[1], n);
253 932680 : if(mpz_sgn(E->buf[2]) == 0 && mpz_sgn(E->buf[3]) == 0){
254 : /* u = 0 <=> Y2*Z1 = Y1*Z2 <=> Y2/Z2 = Y1/Z1*/
255 : /* v = 0 <=> X2*Z1 = X1*Z2 <=> X2/Z2 = X1/Z1*/
256 0 : return pt_w_duplicate(f, x3, y3, z3, x1, y1, z1, n, E);
257 : }
258 : /* Z3:=Z1*Z2 mod p; # M*/
259 932680 : mpres_mul(z3, z1, z2, n);
260 : /* X3:=u^2 mod p;*/
261 932680 : mpres_sqr(x3, E->buf[2], n);
262 : /* X3:=X3*Z3 mod p;*/
263 932680 : mpres_mul(x3, x3, z3, n);
264 : /* R:=v^2 mod p;*/
265 932680 : mpres_sqr(E->buf[4], E->buf[3], n);
266 : /* vvv:=v*R mod p;*/
267 932680 : mpres_mul(E->buf[5], E->buf[3], E->buf[4], n);
268 : /* R:=R*A mod p;*/
269 932680 : mpres_mul(E->buf[4], E->buf[4], E->buf[1], n);
270 : /* Y3:=2*R mod p; # *2*/
271 932680 : mpres_mul_ui(y3, E->buf[4], 2, n);
272 : /* A:=X3-vvv mod p;*/
273 932680 : mpres_sub(E->buf[1], x3, E->buf[5], n);
274 : /* A:=A-Y3 mod p;*/
275 932680 : mpres_sub(E->buf[1], E->buf[1], y3, n);
276 : /* X3:=v*A mod p; # M*/
277 932680 : mpres_mul(x3, E->buf[3], E->buf[1], n);
278 : /* Y3:=R-A mod p;*/
279 932680 : mpres_sub(y3, E->buf[4], E->buf[1], n);
280 : /* Y3:=u*Y3 mod p;*/
281 932680 : mpres_mul(y3, y3, E->buf[2], n);
282 : /* A:=vvv*Y1Z2 mod p;*/
283 932680 : mpres_mul(E->buf[1], E->buf[5], E->buf[0], n);
284 : /* Y3:=Y3-A mod p;*/
285 932680 : mpres_sub(y3, y3, E->buf[1], n);
286 : /* Z3:=vvv*Z3 mod p; # M*/
287 932680 : mpres_mul(z3, z3, E->buf[5], n);
288 932680 : return 1;
289 : }
290 0 : return 0;
291 : }
292 :
293 : #if USE_ADD_SUB_CHAINS > 0
294 : /* [x3, y3, z3] <- [x1, y1, z1] - [x2, y2, z2]; P3 != P1, P3 != P2.
295 : -P2 ~ -(x2/z2, y2/z2, 1) = (x2/z2, -y2/z2-a1*x/z2-a3, 1)
296 : ~ (x2, -y2-a1*x2-a3*z2, z2).
297 : */
298 : int
299 : pt_w_sub(mpz_t f, mpres_t x3, mpres_t y3, mpres_t z3,
300 : mpres_t x1, mpres_t y1, mpres_t z1,
301 : mpres_t x2, mpres_t y2, mpres_t z2,
302 : mpmod_t n, ell_curve_t E)
303 : {
304 : int res = 1;
305 :
306 : if(E->law == ECM_LAW_HOMOGENEOUS){
307 : /* FIXME: does not work for complete equation! */
308 : mpres_neg(y2, y2, n);
309 : res = pt_w_add(f, x3, y3, z3, x1, y1, z1, x2, y2, z2, n, E);
310 : mpres_neg(y2, y2, n);
311 : }
312 : else if(E->law == ECM_LAW_AFFINE){
313 : /* buf[3] not used in law, so use it */
314 : mpres_mul(E->buf[3], E->a1, x2, n);
315 : mpres_add(E->buf[3], E->buf[3], E->a3, n);
316 : mpres_add(E->buf[3], E->buf[3], y2, n);
317 : mpres_neg(E->buf[3], E->buf[3], n);
318 : res = pt_w_add(f, x3, y3, z3, x1, y1, z1, x2, E->buf[3], z2, n, E);
319 : }
320 : return res;
321 : }
322 : #endif
323 :
324 : /******************** projective Hessian form ********************/
325 :
326 : /* U^3+V^3+W^3 = 3*D*U*V*W, D^3 <> 1.
327 : O_H = [1:-1:0]
328 : -[u:v:w] = [v:u:w]
329 : Warning: there can exist two other points at infinity, namely
330 : [1:-omega:0] and [1:-omega^2:0] where omega^3 = 1.
331 : */
332 : int
333 557970 : hessian_is_zero(ell_point_t P, ATTRIBUTE_UNUSED ell_curve_t E, mpmod_t n)
334 : {
335 : mpres_t tmp;
336 : int ret;
337 :
338 557970 : if(mpz_sgn(P->z) != 0)
339 557950 : return 0;
340 20 : mpres_init(tmp, n);
341 20 : mpres_add(tmp, P->x, P->y, n);
342 20 : ret = mpz_sgn(tmp) == 0;
343 : #if 0
344 : if(ret)
345 : gmp_printf("found a third root of unity? %Zd/%Zd\n", P->x, P->y);
346 : #endif
347 20 : mpres_clear(tmp, n);
348 20 : return ret;
349 : }
350 :
351 : void
352 10 : hessian_set_to_zero(ell_point_t P, ATTRIBUTE_UNUSED ell_curve_t E, mpmod_t n)
353 : {
354 10 : mpres_set_si(P->x, 1, n);
355 10 : mpres_set_si(P->y, -1, n);
356 10 : mpres_set_si(P->z, 0, n);
357 10 : }
358 :
359 : #if DEBUG_ADD_LAWS >= 1
360 : void
361 : hessian_print(ell_point_t P, ell_curve_t E, mpmod_t n)
362 : {
363 : pt_w_print(P->x, P->y, P->z, E, n);
364 : }
365 : #endif
366 :
367 : #if USE_ADD_SUB_CHAINS > 0
368 : /* -[u:v:w] = [v:u:w] */
369 : void
370 : hessian_negate(ell_point_t P, ATTRIBUTE_UNUSED ell_curve_t E, ATTRIBUTE_UNUSED mpmod_t n)
371 : {
372 : mpz_swap(P->x, P->y); /* humf */
373 : }
374 : #endif
375 :
376 : /* TODO: decrease the number of buffers? */
377 : int
378 446640 : hessian_duplicate(ell_point_t R, ell_point_t P,
379 : ATTRIBUTE_UNUSED ell_curve_t E, mpmod_t n)
380 : {
381 : /* A = buf[0], ..., G = buf[6], H = buf[7], J = buf[8] */
382 : /* A:=P[1]^2 mod N; */
383 446640 : mpres_mul(E->buf[0], P->x, P->x, n);
384 : /* B:=P[2]^2 mod N; */
385 446640 : mpres_mul(E->buf[1], P->y, P->y, n);
386 : /* C:=P[3]^2 mod N; */
387 446640 : mpres_mul(E->buf[2], P->z, P->z, n);
388 : /* D:=(A+B) mod N; */
389 446640 : mpres_add(E->buf[3], E->buf[0], E->buf[1], n);
390 : /* E:=(A+C) mod N; */
391 446640 : mpres_add(E->buf[4], E->buf[0], E->buf[2], n);
392 : /* F:=(B+C) mod N; */
393 446640 : mpres_add(E->buf[5], E->buf[1], E->buf[2], n);
394 : /* G:=((P[1]+P[2])^2-D) mod N; */
395 446640 : mpres_add(E->buf[6], P->x, P->y, n);
396 446640 : mpres_mul(E->buf[6], E->buf[6], E->buf[6], n);
397 446640 : mpres_sub(E->buf[6], E->buf[6], E->buf[3], n);
398 : /* H:=((P[1]+P[3])^2-E) mod N; */
399 446640 : mpres_add(E->buf[7], P->x, P->z, n);
400 446640 : mpres_mul(E->buf[7], E->buf[7], E->buf[7], n);
401 446640 : mpres_sub(E->buf[7], E->buf[7], E->buf[4], n);
402 : /* J:=((P[2]+P[3])^2-F) mod N; */
403 446640 : mpres_add(E->buf[8], P->y, P->z, n);
404 446640 : mpres_mul(E->buf[8], E->buf[8], E->buf[8], n);
405 446640 : mpres_sub(E->buf[8], E->buf[8], E->buf[5], n);
406 : /* R->x = ((J-G)*(H+2*E)) mod N */
407 446640 : mpres_sub(E->buf[0], E->buf[8], E->buf[6], n);
408 446640 : mpres_add(E->buf[1], E->buf[7], E->buf[4], n);
409 446640 : mpres_add(E->buf[1], E->buf[1], E->buf[4], n);
410 446640 : mpres_mul(R->x, E->buf[0], E->buf[1], n);
411 : /* R->y = ((G-H)*(J+2*F)) mod N */
412 446640 : mpres_sub(E->buf[0], E->buf[6], E->buf[7], n);
413 446640 : mpres_add(E->buf[1], E->buf[8], E->buf[5], n);
414 446640 : mpres_add(E->buf[1], E->buf[1], E->buf[5], n);
415 446640 : mpres_mul(R->y, E->buf[0], E->buf[1], n);
416 : /* R->z = ((H-J)*(G+2*D)) mod N */
417 446640 : mpres_sub(E->buf[0], E->buf[7], E->buf[8], n);
418 446640 : mpres_add(E->buf[1], E->buf[6], E->buf[3], n);
419 446640 : mpres_add(E->buf[1], E->buf[1], E->buf[3], n);
420 446640 : mpres_mul(R->z, E->buf[0], E->buf[1], n);
421 446640 : return 1;
422 : }
423 :
424 : /* TODO: reduce the number of buffers? */
425 : int
426 237720 : hessian_plus(ell_point_t R, ell_point_t P, ell_point_t Q,
427 : ATTRIBUTE_UNUSED ell_curve_t E, mpmod_t n)
428 : {
429 : /* P = [T1,T2,T3], Q = [T4,T5,T6] */
430 : /* P = Q <=> T1/T3=T4/T6 and T2/T3=T5/T6
431 : <=> T1*T6=T3*T4 and T2*T6=T3*T5
432 : */
433 : /* T1 = buf[0], ..., T7 = buf[6] */
434 : /* T7:=(T1*T6) mod N; */
435 237720 : mpres_mul(E->buf[6], P->x, Q->z, n);
436 : /* T1:=(T1*T5) mod N; */
437 237720 : mpres_mul(E->buf[0], P->x, Q->y, n);
438 : /* T5:=(T3*T5) mod N; */
439 237720 : mpres_mul(E->buf[4], P->z, Q->y, n);
440 : /* T3:=(T3*T4) mod N; */
441 237720 : mpres_mul(E->buf[2], P->z, Q->x, n);
442 : /* T4:=(T2*T4) mod N; */
443 237720 : mpres_mul(E->buf[3], P->y, Q->x, n);
444 : /* T2:=(T2*T6) mod N; */
445 237720 : mpres_mul(E->buf[1], P->y, Q->z, n);
446 :
447 237720 : if(mpres_equal(E->buf[6], E->buf[2], n)
448 0 : && mpres_equal(E->buf[4], E->buf[1], n))
449 : /* as a matter of that, P = Q and we need duplicate */
450 0 : return hessian_duplicate(R, P, E, n);
451 :
452 : /* T6:=(T2*T7) mod N; */
453 237720 : mpres_mul(E->buf[5], E->buf[1], E->buf[6], n);
454 : /* T2:=(T2*T4) mod N; */
455 237720 : mpres_mul(E->buf[1], E->buf[1], E->buf[3], n);
456 : /* T4:=(T3*T4) mod N; */
457 237720 : mpres_mul(E->buf[3], E->buf[2], E->buf[3], n);
458 : /* T3:=(T3*T5) mod N; */
459 237720 : mpres_mul(E->buf[2], E->buf[2], E->buf[4], n);
460 : /* T5:=(T1*T5) mod N; */
461 237720 : mpres_mul(E->buf[4], E->buf[0], E->buf[4], n);
462 : /* T1:=(T1*T7) mod N; */
463 237720 : mpres_mul(E->buf[0], E->buf[0], E->buf[6], n);
464 : /* T1:=(T1-T4) mod N; */
465 237720 : mpres_sub(R->y, E->buf[0], E->buf[3], n);
466 : /* T2:=(T2-T5) mod N; */
467 237720 : mpres_sub(R->x, E->buf[1], E->buf[4], n);
468 : /* T3:=(T3-T6) mod N; */
469 237720 : mpres_sub(R->z, E->buf[2], E->buf[5], n);
470 : /* return [T2, T1, T3]; */
471 237720 : return 1;
472 : }
473 :
474 : int
475 237720 : hessian_add(ell_point_t R, ell_point_t P, ell_point_t Q, ell_curve_t E, mpmod_t n)
476 : {
477 237720 : if(hessian_is_zero(P, E, n)){
478 0 : ell_point_set(R, Q, E, n);
479 0 : return 1;
480 : }
481 237720 : else if(hessian_is_zero(Q, E, n)){
482 0 : ell_point_set(R, P, E, n);
483 0 : return 1;
484 : }
485 : else
486 237720 : return hessian_plus(R, P, Q, E, n);
487 : }
488 :
489 : #if USE_ADD_SUB_CHAINS > 0
490 : int
491 : hessian_sub(ell_point_t R, ell_point_t P, ell_point_t Q, ell_curve_t E, mpmod_t n)
492 : {
493 : int ret;
494 :
495 : hessian_negate(Q, E, n);
496 : ret = hessian_add(R, P, Q, E, n);
497 : hessian_negate(Q, E, n);
498 : return ret;
499 : }
500 : #endif
501 :
502 : /* switch from X^3+Y^3+1=3*D*X*Y to Y^2=X^3+A*X+B
503 : A:=-27*D*(D^3+8);
504 : B:=54*(D^6-20*D^3-8);
505 : xi:=12*(D^3-1)/(D*u+v+1);
506 : x:=-9*D^2+xi*u;
507 : y:=3*xi*(v-1);
508 : OUTPUT: If a factor is found during the inversion, it is put in f and
509 : ECM_FACTOR_FOUND_STEP1 is returned. Otherwise, ECM_NO_FACTOR_FOUND is
510 : returned.
511 : SIDE-EFFECT: (x, y, D) <- (x_on_W, y_on_W, A_of_W)
512 : */
513 : int
514 30 : hessian_to_weierstrass(mpz_t f, mpres_t x, mpres_t y, mpres_t D, mpmod_t n)
515 : {
516 : mpres_t D3, A, xi, tmp1, tmp2;
517 30 : int ret = ECM_NO_FACTOR_FOUND;
518 :
519 : #if DEBUG_ADD_LAWS >= 1
520 : printf("P:=[");
521 : print_mpz_from_mpres(x, n);
522 : printf(", ");
523 : print_mpz_from_mpres(y, n);
524 : printf(", 1];\n");
525 : printf("D:=");
526 : print_mpz_from_mpres(D, n);
527 : printf(";\n");
528 : #endif
529 : /* D3 <- D^3 */
530 30 : mpres_init(D3, n);
531 30 : mpres_mul(D3, D, D, n);
532 30 : mpres_mul(D3, D3, D, n);
533 : /* finish A */
534 30 : mpres_init(A, n);
535 30 : mpres_add_ui(A, D3, 8, n);
536 30 : mpres_mul(A, A, D, n);
537 30 : mpres_mul_ui(A, A, 27, n);
538 30 : mpres_neg(A, A, n);
539 : /* compute xi */
540 30 : mpres_init(xi, n);
541 30 : mpres_init(tmp1, n);
542 30 : mpres_mul(tmp1, D, x, n);
543 30 : mpres_add(tmp1, tmp1, y, n);
544 30 : mpres_add_ui(tmp1, tmp1, 1, n);
545 30 : mpres_init(tmp2, n);
546 30 : mpres_sub_ui(tmp2, D3, 1, n);
547 30 : mpres_mul_ui(tmp2, tmp2, 12, n);
548 30 : if(mpres_invert(xi, tmp1, n) == 0){
549 0 : mpres_gcd(f, tmp1, n);
550 0 : ret = ECM_FACTOR_FOUND_STEP1;
551 : }
552 : else{
553 30 : mpres_mul(xi, xi, tmp2, n);
554 : /* compute x */
555 30 : mpres_mul(tmp1, D, D, n);
556 30 : mpres_mul_ui(tmp1, tmp1, 9, n);
557 30 : mpres_mul(tmp2, xi, x, n);
558 30 : mpres_sub(x, tmp2, tmp1, n);
559 : /* compute y */
560 30 : mpres_sub_ui(tmp1, y, 1, n);
561 30 : mpres_mul(tmp1, tmp1, xi, n);
562 30 : mpres_mul_ui(y, tmp1, 3, n);
563 30 : mpres_set(D, A, n);
564 : #if DEBUG_ADD_LAWS >= 1
565 : printf("WP:=[");
566 : print_mpz_from_mpres(x, n);
567 : printf(", ");
568 : print_mpz_from_mpres(y, n);
569 : printf(", 1];\n");
570 : printf("WA:=");
571 : print_mpz_from_mpres(D, n);
572 : printf(";\nWB:=(WP[2]^2-WP[1]^3-WA*WP[1]) mod N;WE:=[WA, WB];\n");
573 : #endif
574 : }
575 30 : mpres_clear(A, n);
576 30 : mpres_clear(D3, n);
577 30 : mpres_clear(xi, n);
578 30 : mpres_clear(tmp1, n);
579 30 : mpres_clear(tmp2, n);
580 30 : return ret;
581 : }
582 :
583 : int
584 30 : mult_by_3(mpz_t f, mpres_t x, mpres_t y, mpres_t A, mpmod_t n)
585 : {
586 : ell_curve_t E;
587 : ell_point_t P, Q;
588 30 : int ret = ECM_NO_FACTOR_FOUND;
589 : mpz_t e;
590 :
591 30 : ell_curve_init_set(E, ECM_EC_TYPE_WEIERSTRASS, ECM_LAW_AFFINE, A, n);
592 30 : ell_point_init(P, E, n);
593 30 : mpres_set(P->x, x, n);
594 30 : mpres_set(P->y, y, n);
595 30 : mpres_set_ui(P->z, 1, n);
596 30 : ell_point_init(Q, E, n);
597 30 : mpz_init_set_ui(e, 3);
598 30 : if(ell_point_mul(f, Q, e, P, E, n) != 0){
599 30 : mpres_set(x, Q->x, n);
600 30 : mpres_set(y, Q->y, n);
601 : }
602 30 : mpz_clear(e);
603 30 : ell_point_clear(Q, E, n);
604 30 : ell_point_clear(P, E, n);
605 30 : ell_curve_clear(E, n);
606 30 : return ret;
607 : }
608 :
609 : /******************** projective twisted Hessian form ********************/
610 :
611 : /* a*U^3+V^3+W^3 = d*U*V*W
612 : O_E = [0:-1:1]
613 : -[U:V:W]=[U:W:V]
614 : */
615 : int
616 1256190 : twisted_hessian_is_zero(ell_point_t P, ATTRIBUTE_UNUSED ell_curve_t E, mpmod_t n)
617 : {
618 : mpres_t tmp;
619 : int ret;
620 :
621 1256190 : if(mpz_sgn(P->x) != 0)
622 1256190 : return 0;
623 0 : mpres_init(tmp, n);
624 0 : mpres_add(tmp, P->y, P->z, n);
625 0 : ret = mpz_sgn(tmp) == 0;
626 : #if 0
627 : if(ret)
628 : gmp_printf("found a third root of unity? %Zd/%Zd\n", P->x, P->y);
629 : #endif
630 0 : mpres_clear(tmp, n);
631 0 : return ret;
632 : }
633 :
634 : void
635 0 : twisted_hessian_set_to_zero(ell_point_t P, ATTRIBUTE_UNUSED ell_curve_t E, mpmod_t n)
636 : {
637 0 : mpres_set_si(P->x, 0, n);
638 0 : mpres_set_si(P->y, -1, n);
639 0 : mpres_set_si(P->z, 1, n);
640 0 : }
641 :
642 : #if DEBUG_ADD_LAWS >= 1
643 : void
644 : twisted_hessian_print(ell_point_t P, ell_curve_t E, mpmod_t n)
645 : {
646 : pt_w_print(P->x, P->y, P->z, E, n);
647 : }
648 : #endif
649 :
650 : #if USE_ADD_SUB_CHAINS > 0
651 : /* -[u:v:w] = [u:w:v] */
652 : void
653 : twisted_hessian_negate(ell_point_t P, ATTRIBUTE_UNUSED ell_curve_t E, ATTRIBUTE_UNUSED mpmod_t n)
654 : {
655 : mpz_swap(P->y, P->z); /* humf */
656 : }
657 : #endif
658 :
659 : /* TODO: decrease the number of buffers? */
660 : /* 6M+2S+1M_d: better when d is small */
661 : int
662 1036340 : twisted_hessian_duplicate(ell_point_t R, ell_point_t P,
663 : ATTRIBUTE_UNUSED ell_curve_t E, mpmod_t n)
664 : {
665 : /* R = buf[0], ..., W = buf[5], C = buf[6], D = buf[7], E = buf[8] */
666 : /* R:=Y1+Z1;*/
667 1036340 : mpres_add(E->buf[0], P->y, P->z, n);
668 : /* S:=Y1-Z1;*/
669 1036340 : mpres_sub(E->buf[1], P->y, P->z, n);
670 : /* T:=R^2 mod N;*/
671 1036340 : mpres_sqr(E->buf[2], E->buf[0], n);
672 : /* U:=S^2 mod N;*/
673 1036340 : mpres_sqr(E->buf[3], E->buf[1], n);
674 : /* V:=T+3*U;*/
675 1036340 : mpres_add(E->buf[4], E->buf[2], E->buf[3], n);
676 1036340 : mpres_add(E->buf[4], E->buf[4], E->buf[3], n);
677 1036340 : mpres_add(E->buf[4], E->buf[4], E->buf[3], n);
678 : /* W:=3*T+U;*/
679 1036340 : mpres_add(E->buf[5], E->buf[3], E->buf[2], n);
680 1036340 : mpres_add(E->buf[5], E->buf[5], E->buf[2], n);
681 1036340 : mpres_add(E->buf[5], E->buf[5], E->buf[2], n);
682 : /* C:=(R*V) mod N;*/
683 1036340 : mpres_mul(E->buf[6], E->buf[0], E->buf[4], n);
684 : /* D:=(S*W) mod N;*/
685 1036340 : mpres_mul(E->buf[7], E->buf[1], E->buf[5], n);
686 : /* E:=(3*C-E0[2]*X1*(W-V)) mod N;*/
687 1036340 : mpres_sub(E->buf[8], E->buf[5], E->buf[4], n);
688 1036340 : mpres_mul(E->buf[8], E->buf[8], P->x, n);
689 1036340 : mpres_mul(E->buf[8], E->buf[8], E->a6, n);
690 1036340 : mpres_sub(E->buf[8], E->buf[6], E->buf[8], n);
691 1036340 : mpres_add(E->buf[8], E->buf[8], E->buf[6], n);
692 1036340 : mpres_add(E->buf[8], E->buf[8], E->buf[6], n);
693 : /* X3:=(-2*X1*D) mod N;*/
694 1036340 : mpres_mul(R->x, P->x, E->buf[7], n);
695 1036340 : mpres_add(R->x, R->x, R->x, n);
696 1036340 : mpres_neg(R->x, R->x, n);
697 : /* Y3:=((D+E)*Z1) mod N;*/
698 1036340 : mpres_add(E->buf[0], E->buf[7], E->buf[8], n);
699 1036340 : mpres_mul(E->buf[1], E->buf[0], P->z, n);
700 : /* Z3:=((D-E)*Y1) mod N;*/
701 1036340 : mpres_sub(E->buf[0], E->buf[7], E->buf[8], n);
702 1036340 : mpres_mul(R->z, E->buf[0], P->y, n);
703 1036340 : mpres_set(R->y, E->buf[1], n);
704 1036340 : return 1;
705 : }
706 :
707 : /* TODO: reduce the number of buffers? */
708 : int
709 549780 : twisted_hessian_plus(ell_point_t R, ell_point_t P, ell_point_t Q,
710 : ATTRIBUTE_UNUSED ell_curve_t E, mpmod_t n)
711 : {
712 : /* A = buf[0], ... F = buf[5], G = [6], H = [7], J = [8] */
713 : // A:=X1*Z2 mod N;
714 549780 : mpres_mul(E->buf[0], P->x, Q->z, n);
715 : // B:=Z1*Z2 mod N;
716 549780 : mpres_mul(E->buf[1], P->z, Q->z, n);
717 : // C:=Y1*X2 mod N;
718 549780 : mpres_mul(E->buf[2], P->y, Q->x, n);
719 : // D:=Y1*Y2 mod N;
720 549780 : mpres_mul(E->buf[3], P->y, Q->y, n);
721 : // E:=Z1*Y2 mod N;
722 549780 : mpres_mul(E->buf[4], P->z, Q->y, n);
723 : // F:=E0[1]*X1*X2 mod N;
724 549780 : mpres_mul(E->buf[5], P->x, Q->x, n);
725 549780 : mpres_mul(E->buf[5], E->buf[5], E->a4, n);
726 : // Hisil
727 : // G := (D+B)*(A-C) mod N;
728 549780 : mpres_add(E->buf[9], E->buf[3], E->buf[1], n);
729 549780 : mpres_sub(E->buf[6], E->buf[0], E->buf[2], n);
730 549780 : mpres_mul(E->buf[6], E->buf[6], E->buf[9], n);
731 : // H := (D-B)*(A+C) mod N;
732 549780 : mpres_sub(E->buf[9], E->buf[3], E->buf[1], n);
733 549780 : mpres_add(E->buf[7], E->buf[0], E->buf[2], n);
734 549780 : mpres_mul(E->buf[7], E->buf[7], E->buf[9], n);
735 : // J := (D+F)*(A-E) mod N;
736 549780 : mpres_add(E->buf[9], E->buf[3], E->buf[5], n);
737 549780 : mpres_sub(E->buf[8], E->buf[0], E->buf[4], n);
738 549780 : mpres_mul(E->buf[8], E->buf[8], E->buf[9], n);
739 : // K := (D-F)*(A+E) mod N;
740 : // this is the last use of A, so that K -> buf[0]
741 549780 : mpres_sub(E->buf[9], E->buf[3], E->buf[5], n);
742 549780 : mpres_add(E->buf[0], E->buf[0], E->buf[4], n);
743 549780 : mpres_mul(E->buf[0], E->buf[0], E->buf[9], n);
744 : // X3 := G-H
745 549780 : mpres_sub(R->x, E->buf[6], E->buf[7], n);
746 : // Y3 := K-J
747 549780 : mpres_sub(R->y, E->buf[0], E->buf[8], n);
748 : // Z3 := (J+K-G-H-2*(B-F)*(C+E)) mod N;
749 549780 : mpres_sub(E->buf[9], E->buf[1], E->buf[5], n);
750 549780 : mpres_add(R->z, E->buf[2], E->buf[4], n);
751 549780 : mpres_mul(R->z, R->z, E->buf[9], n);
752 549780 : mpres_add(R->z, R->z, R->z, n);
753 549780 : mpres_add(R->z, R->z, E->buf[7], n);
754 549780 : mpres_add(R->z, R->z, E->buf[6], n);
755 549780 : mpres_sub(R->z, E->buf[0], R->z, n);
756 549780 : mpres_add(R->z, R->z, E->buf[8], n);
757 549780 : if(mpz_sgn(R->x) == 0 && mpz_sgn(R->y) == 0 && mpz_sgn(R->z) == 0){
758 : // iff (X2:Y2:Z2)=(Z1:gamma^2*X1:gamma*Y1), gamma^3 = a
759 0 : fprintf(stderr, "GASP: X3, Y3 and Z3 are 0\n");
760 0 : exit(-1);
761 : #if 0
762 : // TODO: rewrite with above quantities!
763 : X3:=(X1^2*Y2*Z2-X2^2*Y1*Z1) mod N;
764 : // A*X1*Y2-C*X2*Z1 = A*U-C*V
765 : Y3:=(Z1^2*X2*Y2-Z2^2*X1*Y1) mod N;
766 : // E*Z1*X2-A*Z2*Y1 = E*V-A*W
767 : Z3:=(Y1^2*X2*Z2-Y2^2*X1*Z1) mod N;
768 : // C*Y1*Z2-E*Y2*X1 = C*W-E*U
769 :
770 : // X3 = Y1*(a*X1^3-Z1^3)
771 : // Y3 = g^2*X1*(Z1^3-Y1^3)
772 : // Z3 = g*Z1*(Y1^3-Z1^3)
773 :
774 : // can be made faster with a = aa^3, since then g = aa and we
775 : // can share many things
776 :
777 : #endif
778 : }
779 549780 : return 1;
780 : }
781 :
782 : int
783 549780 : twisted_hessian_add(ell_point_t R, ell_point_t P, ell_point_t Q, ell_curve_t E, mpmod_t n)
784 : {
785 549780 : if(twisted_hessian_is_zero(P, E, n)){
786 0 : ell_point_set(R, Q, E, n);
787 0 : return 1;
788 : }
789 549780 : else if(twisted_hessian_is_zero(Q, E, n)){
790 0 : ell_point_set(R, P, E, n);
791 0 : return 1;
792 : }
793 : else
794 549780 : return twisted_hessian_plus(R, P, Q, E, n);
795 : }
796 :
797 : #if USE_ADD_SUB_CHAINS > 0
798 : int
799 : twisted_hessian_sub(ell_point_t R, ell_point_t P, ell_point_t Q, ell_curve_t E, mpmod_t n)
800 : {
801 : int ret;
802 :
803 : twisted_hessian_negate(Q, E, n);
804 : ret = twisted_hessian_add(R, P, Q, E, n);
805 : twisted_hessian_negate(Q, E, n);
806 : return ret;
807 : }
808 : #endif
809 :
810 : /* INPUT: a*x^3+y^3+1 = d*x*y
811 : OUTPUT: Y^2 = X^3+A*X+B
812 : If a=c^3, then curve isom to Hessian (c*x)^3+y^3+1=3*(d/(3*c))*(c*x)*y
813 : SIDE EFFECT: (x, y, c) <- (x_on_W, y_on_W, A_of_W)
814 : */
815 : int
816 10 : twisted_hessian_to_weierstrass(mpz_t f, mpres_t x, mpres_t y, mpres_t c, mpres_t d, mpmod_t n)
817 : {
818 10 : int ret = ECM_NO_FACTOR_FOUND;
819 : mpres_t tmp;
820 :
821 : #if DEBUG_ADD_LAWS >= 2
822 : printf("x_tH="); print_mpz_from_mpres(x, n); printf("\n");
823 : printf("y_tH="); print_mpz_from_mpres(y, n); printf("\n");
824 : printf("c_tH="); print_mpz_from_mpres(c, n); printf("\n");
825 : printf("d_tH="); print_mpz_from_mpres(d, n); printf("\n");
826 : #endif
827 10 : mpres_init(tmp, n);
828 10 : mpres_mul_ui(tmp, c, 3, n);
829 10 : if(mpres_invert(tmp, tmp, n) == 0){
830 0 : mpres_gcd(f, tmp, n);
831 0 : ret = ECM_FACTOR_FOUND_STEP1;
832 : }
833 : else{
834 10 : mpres_mul(x, x, c, n);
835 10 : mpres_mul(c, tmp, d, n);
836 : /* from x^3+y^3+1=3*c*x*y to Weierstrass stuff */
837 10 : ret = hessian_to_weierstrass(f, x, y, c, n);
838 : #if DEBUG_ADD_LAWS >= 2
839 : printf("A_W="); print_mpz_from_mpres(c, n); printf("\n");
840 : printf("x_W="); print_mpz_from_mpres(x, n); printf("\n");
841 : printf("y_W="); print_mpz_from_mpres(y, n); printf("\n");
842 : #endif
843 : }
844 10 : mpres_clear(tmp, n);
845 10 : return ret;
846 : }
847 :
848 : /******************** generic ec's ********************/
849 :
850 : void
851 376550 : ell_point_init(ell_point_t P, ell_curve_t E, mpmod_t n)
852 : {
853 376550 : mpres_init(P->x, n);
854 376550 : mpres_init(P->y, n);
855 376550 : mpres_init(P->z, n);
856 376550 : if(E->type == ECM_EC_TYPE_WEIERSTRASS){
857 254300 : if(E->law == ECM_LAW_AFFINE)
858 118780 : mpz_set_ui(P->z, 1); /* humf */
859 135520 : else if(E->law == ECM_LAW_HOMOGENEOUS)
860 135520 : mpres_set_ui(P->z, 1, n);
861 : }
862 122250 : else if(E->type == ECM_EC_TYPE_HESSIAN
863 79510 : || E->type == ECM_EC_TYPE_TWISTED_HESSIAN)
864 122230 : mpres_set_ui(P->z, 1, n);
865 376550 : }
866 :
867 : /* TODO: change this according to E->type */
868 : void
869 376200 : ell_point_clear(ell_point_t P, ATTRIBUTE_UNUSED ell_curve_t E, mpmod_t n)
870 : {
871 376200 : mpres_clear(P->x, n);
872 376200 : mpres_clear(P->y, n);
873 376200 : mpres_clear(P->z, n);
874 376200 : }
875 :
876 : #if DEBUG_ADD_LAWS >= 1
877 : void
878 : ell_point_print(ell_point_t P, ell_curve_t E, mpmod_t n)
879 : {
880 : if(E->type == ECM_EC_TYPE_WEIERSTRASS)
881 : pt_w_print(P->x, P->y, P->z, E, n);
882 : else if(E->type == ECM_EC_TYPE_HESSIAN)
883 : hessian_print(P, E, n);
884 : else if(E->type == ECM_EC_TYPE_TWISTED_HESSIAN)
885 : twisted_hessian_print(P, E, n);
886 : }
887 : #endif
888 :
889 : /* TODO: should depend on E->type... */
890 : void
891 1129820 : ell_point_set(ell_point_t Q, ell_point_t P,
892 : ATTRIBUTE_UNUSED ell_curve_t E, ATTRIBUTE_UNUSED mpmod_t n)
893 : {
894 1129820 : mpres_set(Q->x, P->x, n);
895 1129820 : mpres_set(Q->y, P->y, n);
896 1129820 : mpres_set(Q->z, P->z, n);
897 1129820 : }
898 :
899 : void
900 2006 : ell_curve_init(ell_curve_t E, int etype, int law, mpmod_t n)
901 : {
902 : int i;
903 :
904 2006 : E->type = etype;
905 2006 : E->law = law;
906 2006 : mpres_init(E->a1, n);
907 2006 : mpres_init(E->a3, n);
908 2006 : mpres_init(E->a2, n);
909 2006 : mpres_init(E->a4, n);
910 2006 : mpres_init(E->a6, n);
911 2006 : mpres_set_ui(E->a1, 0, n);
912 2006 : mpres_set_ui(E->a3, 0, n);
913 2006 : mpres_set_ui(E->a2, 0, n);
914 2006 : mpres_set_ui(E->a4, 0, n);
915 2006 : mpres_set_ui(E->a6, 0, n);
916 22066 : for(i = 0; i < EC_W_NBUFS; i++)
917 20060 : mpres_init (E->buf[i], n);
918 2006 : }
919 :
920 : void
921 330 : ell_curve_init_set(ell_curve_t E, int etype, int law, mpres_t A, mpmod_t n)
922 : {
923 330 : ell_curve_init(E, etype, law, n);
924 330 : mpres_set(E->a4, A, n);
925 330 : }
926 :
927 : void
928 1566 : ell_curve_set_z(ell_curve_t E, ell_curve_t zE, mpmod_t n)
929 : {
930 1566 : ell_curve_init(E, zE->type, zE->law, n);
931 1566 : mpres_set_z(E->a1, zE->a1, n);
932 1566 : mpres_set_z(E->a3, zE->a3, n);
933 1566 : mpres_set_z(E->a2, zE->a2, n);
934 1566 : mpres_set_z(E->a4, zE->a4, n);
935 1566 : mpres_set_z(E->a6, zE->a6, n);
936 : #if 0
937 : E->disc = zE->disc;
938 : if(E->disc != 0){
939 : mpres_init(E->sq[0], n);
940 : mpres_set_z(E->sq[0], zE->sq[0], n);
941 : }
942 : #endif
943 1566 : }
944 :
945 : void
946 1926 : ell_curve_clear(ell_curve_t E, mpmod_t n)
947 : {
948 : int i;
949 :
950 1926 : mpres_clear(E->a4, n);
951 21186 : for(i = 0; i < EC_W_NBUFS; i++)
952 19260 : mpres_clear (E->buf[i], n);
953 : /* TODO: case of sq */
954 1926 : }
955 :
956 : #if DEBUG_ADD_LAWS >= 1
957 : void
958 : ell_curve_print(ell_curve_t E, mpmod_t n)
959 : {
960 : if(E->type == ECM_EC_TYPE_WEIERSTRASS){
961 : printf("["); print_mpz_from_mpres(E->a1, n);
962 : printf(", "); print_mpz_from_mpres(E->a3, n);
963 : printf(", "); print_mpz_from_mpres(E->a2, n);
964 : printf(", "); print_mpz_from_mpres(E->a4, n);
965 : printf(", "); print_mpz_from_mpres(E->a6, n); printf("];\n");
966 : }
967 : else if(E->type == ECM_EC_TYPE_HESSIAN){
968 : printf("D:="); print_mpz_from_mpres(E->a4, n); printf(";\n");
969 : printf("E:=[D];\n");
970 : }
971 : else if(E->type == ECM_EC_TYPE_TWISTED_HESSIAN){
972 : printf("a:="); print_mpz_from_mpres(E->a4, n); printf(";\n");
973 : printf("d:="); print_mpz_from_mpres(E->a6, n); printf(";\n");
974 : printf("E:=[a, d];\n");
975 : }
976 : }
977 : #endif
978 :
979 : /* OUTPUT: 1 if P = O_E, 0 otherwise. */
980 : int
981 740080 : ell_point_is_zero(ell_point_t P, ell_curve_t E, mpmod_t n)
982 : {
983 740080 : if(E->type == ECM_EC_TYPE_WEIERSTRASS)
984 500920 : return pt_w_is_zero(P->z, n);
985 239160 : else if(E->type == ECM_EC_TYPE_HESSIAN)
986 82530 : return hessian_is_zero(P, E, n);
987 156630 : else if(E->type == ECM_EC_TYPE_TWISTED_HESSIAN)
988 156630 : return twisted_hessian_is_zero(P, E, n);
989 0 : return ECM_ERROR;
990 : }
991 :
992 : void
993 60 : ell_point_set_to_zero(ell_point_t P, ell_curve_t E, mpmod_t n)
994 : {
995 60 : if(E->type == ECM_EC_TYPE_WEIERSTRASS)
996 50 : pt_w_set_to_zero(P, n);
997 10 : else if(E->type == ECM_EC_TYPE_HESSIAN)
998 10 : hessian_set_to_zero(P, E, n);
999 0 : else if(E->type == ECM_EC_TYPE_TWISTED_HESSIAN)
1000 0 : twisted_hessian_set_to_zero(P, E, n);
1001 60 : }
1002 :
1003 : int
1004 100 : ell_point_is_on_curve(ell_point_t P, ell_curve_t E, mpmod_t n)
1005 : {
1006 100 : int ok = 1;
1007 :
1008 100 : if(ell_point_is_zero(P, E, n))
1009 0 : return 1;
1010 100 : if(E->type == ECM_EC_TYPE_WEIERSTRASS){
1011 : mpres_t tmp1, tmp2;
1012 :
1013 100 : mpres_init(tmp1, n);
1014 100 : mpres_init(tmp2, n);
1015 100 : if(E->law == ECM_LAW_AFFINE){
1016 : /* y^2+a1*x*y+a3*y = x^3+a2*x^2+a4*x+a6? */
1017 100 : mpres_mul(tmp1, E->a1, P->x, n);
1018 100 : mpres_add(tmp1, tmp1, P->y, n);
1019 100 : mpres_add(tmp1, tmp1, E->a3, n);
1020 100 : mpres_mul(tmp1, tmp1, P->y, n);
1021 :
1022 100 : mpres_add(tmp2, E->a2, P->x, n);
1023 100 : mpres_mul(tmp2, tmp2, P->x, n);
1024 100 : mpres_add(tmp2, tmp2, E->a4, n);
1025 100 : mpres_mul(tmp2, tmp2, P->x, n);
1026 100 : mpres_add(tmp2, tmp2, E->a6, n);
1027 : }
1028 : #if 0 // useless for the time being
1029 : else{
1030 : /* y^2*z+a1*x*y*z+a3*y*z^2 = x^3+a2*x^2*z+a4*x*z^2+a6*z^3? */
1031 : /* y*z*(y+a1*x+a3*z) = ((x+a2*z)*x+a4*z^2)*x+a6*z^3? */
1032 : mpres_t tmp3;
1033 :
1034 : mpres_mul(tmp1, E->a1, P->x, n); /* a1*x */
1035 : mpres_add(tmp1, tmp1, P->y, n); /* a1*x+y */
1036 : mpres_mul(tmp2, E->a3, P->z, n); /* a3*z */
1037 : mpres_add(tmp1, tmp1, tmp2, n); /* y+a1*x+a3*z */
1038 : mpres_mul(tmp1, tmp1, P->y, n); /* y*(...) */
1039 : mpres_mul(tmp1, tmp1, P->z, n); /* lhs */
1040 :
1041 : mpres_init(tmp3, n);
1042 : mpres_mul(tmp2, E->a2, P->z, n); /* a2*z */
1043 : mpres_add(tmp2, tmp2, P->x, n); /* x+a2*z */
1044 : mpres_mul(tmp2, tmp2, P->x, n); /* (x+a2*z)*x */
1045 : mpres_mul(tmp3, E->a4, P->z, n); /* a4*z */
1046 : mpres_mul(tmp3, tmp3, P->z, n); /* a4*z^2 */
1047 : mpres_add(tmp2, tmp2, tmp3, n); /* (x+a2*z)*x+a4*z^2 */
1048 : mpres_mul(tmp2, tmp2, P->x, n); /* (...)*x */
1049 : mpres_mul(tmp3, P->z, P->z, n); /* z^2 */
1050 : mpres_mul(tmp3, tmp3, P->z, n); /* z^3 */
1051 : mpres_mul(tmp3, tmp3, E->a6, n); /* a6*z^3 */
1052 : mpres_add(tmp2, tmp2, tmp3, n); /* rhs */
1053 : mpres_clear(tmp3, n);
1054 : }
1055 : #endif
1056 100 : ok = mpres_equal(tmp1, tmp2, n);
1057 :
1058 100 : mpres_clear(tmp1, n);
1059 100 : mpres_clear(tmp2, n);
1060 : }
1061 0 : else if(E->type == ECM_EC_TYPE_HESSIAN){
1062 : /* TODO */
1063 : }
1064 0 : else if(E->type == ECM_EC_TYPE_TWISTED_HESSIAN){
1065 : /* TODO */
1066 : }
1067 100 : return ok;
1068 : }
1069 :
1070 : #if DEBUG_ADD_LAWS >= 1
1071 : static void
1072 : ell_point_check(ell_point_t P, ell_curve_t E, mpmod_t n)
1073 : {
1074 : if(ell_point_is_on_curve(P, E, n) == 0){
1075 : printf("Point not on curve\n");
1076 : printf("E:=");
1077 : ell_curve_print(E, n);
1078 : printf("P:=");
1079 : pt_print(E, P, n);
1080 : printf("\n");
1081 : exit(-1);
1082 : }
1083 : }
1084 : #endif
1085 :
1086 : #if DEBUG_ADD_LAWS >= 1
1087 : int
1088 : ell_point_equal(ell_point_t P, ell_point_t Q, ell_curve_t E, mpmod_t n)
1089 : {
1090 : int ret = 1;
1091 :
1092 : if(E->type == ECM_EC_TYPE_WEIERSTRASS){
1093 : if(E->law == ECM_LAW_AFFINE)
1094 : return mpres_equal(P->x, Q->x, n)
1095 : && mpres_equal(P->y, Q->y, n)
1096 : && mpres_equal(P->z, Q->z, n);
1097 : else if(E->law == ECM_LAW_HOMOGENEOUS){
1098 : mpres_t tmp1, tmp2;
1099 :
1100 : mpres_init(tmp1, n);
1101 : mpres_init(tmp2, n);
1102 : mpres_mul(tmp1, P->x, Q->z, n);
1103 : mpres_mul(tmp2, P->z, Q->x, n);
1104 : if(mpres_equal(tmp1, tmp2, n) == 0){
1105 : printf("Px/Pz != Qx/Qz\n");
1106 : ret = 0;
1107 : exit(-1);
1108 : }
1109 : else{
1110 : mpres_mul(tmp1, P->y, Q->z, n);
1111 : mpres_mul(tmp2, P->z, Q->y, n);
1112 : if(mpres_equal(tmp1, tmp2, n) == 0){
1113 : printf("Py/Pz != Qy/Qz\n");
1114 : ret = 0;
1115 : exit(-1);
1116 : }
1117 : }
1118 : mpres_clear(tmp1, n);
1119 : mpres_clear(tmp2, n);
1120 : }
1121 : }
1122 : return ret;
1123 : }
1124 : #endif
1125 :
1126 : /* OUTPUT: 1 if everything ok, 0 otherwise */
1127 : int
1128 2560420 : ell_point_add(mpz_t f, ell_point_t R, ell_point_t P, ell_point_t Q, ell_curve_t E, mpmod_t n)
1129 : {
1130 2560420 : if(E->type == ECM_EC_TYPE_WEIERSTRASS)
1131 1772920 : return pt_w_add(f, R->x, R->y, R->z, P->x, P->y, P->z,
1132 1772920 : Q->x, Q->y, Q->z, n, E);
1133 787500 : else if(E->type == ECM_EC_TYPE_HESSIAN)
1134 237720 : return hessian_add(R, P, Q, E, n);
1135 549780 : else if(E->type == ECM_EC_TYPE_TWISTED_HESSIAN)
1136 549780 : return twisted_hessian_add(R, P, Q, E, n);
1137 : else
1138 0 : return ECM_ERROR;
1139 : }
1140 :
1141 : #if USE_ADD_SUB_CHAINS > 0
1142 : /* R <- P-Q */
1143 : int
1144 : ell_point_sub(mpz_t f, ell_point_t R, ell_point_t P, ell_point_t Q, ell_curve_t E, mpmod_t n)
1145 : {
1146 : if(E->type == ECM_EC_TYPE_WEIERSTRASS)
1147 : return pt_w_sub(f, R->x, R->y, R->z, P->x, P->y, P->z,
1148 : Q->x, Q->y, Q->z, n, E);
1149 : else if(E->type == ECM_EC_TYPE_HESSIAN)
1150 : return hessian_sub(R, P, Q, E, n);
1151 : else if(E->type == ECM_EC_TYPE_TWISTED_HESSIAN)
1152 : return twisted_hessian_sub(R, P, Q, E, n);
1153 : else
1154 : return ECM_ERROR;
1155 : }
1156 : #endif
1157 :
1158 : int
1159 4873060 : ell_point_duplicate(mpz_t f, ell_point_t R, ell_point_t P, ell_curve_t E, mpmod_t n)
1160 : {
1161 : #if DEBUG_ADD_LAWS >= 2
1162 : printf("E:=");
1163 : ell_curve_print(E, n);
1164 : #endif
1165 4873060 : if(E->type == ECM_EC_TYPE_WEIERSTRASS)
1166 3390080 : return pt_w_duplicate(f, R->x, R->y, R->z, P->x, P->y, P->z, n, E);
1167 1482980 : else if(E->type == ECM_EC_TYPE_HESSIAN)
1168 446640 : return hessian_duplicate(R, P, E, n);
1169 1036340 : else if(E->type == ECM_EC_TYPE_TWISTED_HESSIAN)
1170 1036340 : return twisted_hessian_duplicate(R, P, E, n);
1171 : else
1172 0 : return ECM_ERROR;
1173 : }
1174 :
1175 : void
1176 20 : ell_point_negate(ell_point_t P, ell_curve_t E, mpmod_t n)
1177 : {
1178 : #if DEBUG_ADD_LAWS >= 2
1179 : printf("P:="); ell_point_print(P, E, n); printf(";\n");
1180 : #endif
1181 20 : if(ell_point_is_zero(P, E, n) == 0){
1182 20 : if(E->type == ECM_EC_TYPE_WEIERSTRASS){
1183 20 : if(E->law == ECM_LAW_HOMOGENEOUS){
1184 : /* FIXME: does not work for complete equation! */
1185 0 : mpres_neg(P->y, P->y, n);
1186 : }
1187 20 : else if(E->law == ECM_LAW_AFFINE){
1188 : /* (-P).y = -P.y-a1*P.x-a3 */
1189 20 : if(mpz_sgn(E->a1) != 0
1190 20 : || mpz_sgn(E->a3) != 0
1191 20 : || mpz_sgn(E->a2) != 0){ /* FIXME */
1192 0 : printf("GROUMF\n");
1193 0 : exit(-1);
1194 : }
1195 20 : mpres_neg(P->y, P->y, n);
1196 : }
1197 : }
1198 : #if USE_ADD_SUB_CHAINS > 0
1199 : else if(E->type == ECM_EC_TYPE_HESSIAN)
1200 : hessian_negate(P, E, n);
1201 : else if(E->type == ECM_EC_TYPE_TWISTED_HESSIAN)
1202 : twisted_hessian_negate(P, E, n);
1203 : #endif
1204 : }
1205 : #if DEBUG_ADD_LAWS >= 2
1206 : printf("neg(P):="); ell_point_print(P, E, n); printf(";\n");
1207 : #endif
1208 20 : }
1209 :
1210 : /* Q <- [e]*P
1211 : Return value: 0 if a factor is found, and the factor is in Q->x,
1212 : 1 otherwise.
1213 : */
1214 : int
1215 375860 : ell_point_mul_plain (mpz_t f, ell_point_t Q, mpz_t e, ell_point_t P, ell_curve_t E, mpmod_t n)
1216 : {
1217 : size_t l;
1218 375860 : int negated = 0, status = 1;
1219 : ell_point_t P0;
1220 :
1221 375860 : if(ell_point_is_zero(P, E, n) != 0){
1222 300 : ell_point_set(Q, P, E, n);
1223 300 : return 1;
1224 : }
1225 :
1226 375560 : if (mpz_sgn (e) == 0)
1227 : {
1228 20 : ell_point_set_to_zero(Q, E, n);
1229 20 : return 1;
1230 : }
1231 :
1232 375540 : if (mpz_sgn (e) < 0)
1233 : {
1234 10 : negated = 1;
1235 10 : mpz_neg (e, e);
1236 10 : ell_point_negate(P, E, n); /* since the point is non-zero */
1237 : }
1238 :
1239 375540 : if (mpz_cmp_ui (e, 1) == 0){
1240 290 : ell_point_set(Q, P, E, n);
1241 290 : goto ell_point_mul_plain_end;
1242 : }
1243 :
1244 375250 : l = mpz_sizeinbase (e, 2) - 1; /* l >= 1 */
1245 :
1246 375250 : ell_point_init(P0, E, n);
1247 375250 : ell_point_set(P0, P, E, n);
1248 :
1249 : #if DEBUG_ADD_LAWS >= 2
1250 : printf("P:="); ell_point_print(P, E, n); printf(";\n");
1251 : #endif
1252 5244900 : while (l-- > 0)
1253 : {
1254 : #if DEBUG_ADD_LAWS >= 2
1255 : printf("P0:="); ell_point_print(P0, E, n); printf(";\n");
1256 : #endif
1257 4869740 : if(ell_point_duplicate (f, P0, P0, E, n) == 0)
1258 : {
1259 20 : status = 0;
1260 20 : break;
1261 : }
1262 : #if DEBUG_ADD_LAWS >= 2
1263 : printf("Rdup:="); ell_point_print(P0, E, n); printf(";\n");
1264 : printf("dup:=ProjEcmDouble(P0, E, N); ProjEcmEqual(dup, Rdup, N);\n");
1265 : #endif
1266 4869720 : if (mpz_tstbit (e, l))
1267 : {
1268 2560260 : if(ell_point_add (f, P0, P0, P, E, n) == 0)
1269 : {
1270 70 : status = 0;
1271 70 : break;
1272 : }
1273 : #if DEBUG_ADD_LAWS >= 2
1274 : printf("Radd:="); ell_point_print(P0, E, n); printf(";\n");
1275 : printf("Padd:=ProjEcmAdd(P, Rdup, E, N); ProjEcmEqual(Padd, Radd, N);\n");
1276 : #endif
1277 : }
1278 : }
1279 :
1280 375250 : ell_point_set(Q, P0, E, n);
1281 375250 : ell_point_clear(P0, E, n);
1282 375540 : ell_point_mul_plain_end:
1283 :
1284 : /* Undo negation to avoid changing the caller's e value */
1285 375540 : if (negated){
1286 10 : mpz_neg (e, e);
1287 10 : ell_point_negate(P, E, n);
1288 : }
1289 375540 : return status;
1290 : }
1291 :
1292 : int
1293 375860 : ell_point_mul(mpz_t f, ell_point_t Q, mpz_t e, ell_point_t P, ell_curve_t E, mpmod_t n)
1294 : {
1295 : #if 1 /* keeping it simple */
1296 375860 : return ell_point_mul_plain(f, Q, e, P, E, n);
1297 : #else
1298 : return ell_point_mul_add_sub(f, Q, e, P, E, n);
1299 : #endif
1300 : }
1301 :
|
