Line data Source code
1 : /* torsions.c - ECM with special torsion curves
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 : #include "torsions.h"
19 :
20 : #define DEBUG_TORSION 0
21 :
22 : /* We use three variants of Weierstrass parametrization:
23 : CW (complete): y^2+a1*x*y+a3*y=x^3+a2*x^2+a4*x+a6
24 : MW (medium) : y^2=x^3+a2*x^2+a4*x+a6
25 : SW (short) : y^2=x^3+a4*x+a6
26 :
27 : A Kubert curve is the special case Y^2+(1-c)*X*Y-b*Y = X^3-b*X^2
28 :
29 : Generally, we build a curve under the SW form, with affine law, meaning
30 : that constructed points will be [x, y, 1].
31 : */
32 :
33 : /********** utilities **********/
34 :
35 : void
36 330 : mod_div_2(mpz_t x, mpz_t n)
37 : {
38 330 : if(mpz_tstbit(x, 0)){
39 : /* x is odd, x/2 = (x+N)/2 */
40 220 : mpz_add(x, x, n);
41 220 : mpz_div_2exp(x, x, 1);
42 : }
43 : else
44 : /* x is even, x/2 is easy */
45 110 : mpz_div_2exp(x, x, 1);
46 330 : }
47 :
48 : /* r <- q mod N.
49 : Return value: 1 if den invertible, 0 if factor found; in this case
50 : gcd(den(q), N) is put in r.
51 : */
52 : int
53 2070 : mod_from_rat2(mpz_t r, mpz_t num, mpz_t den, mpz_t N)
54 : {
55 2070 : int ret = 1;
56 :
57 2070 : if(mpz_invert(r, den, N) == 0){
58 100 : mpz_gcd(r, den, N);
59 100 : ret = 0;
60 : }
61 : else{
62 1970 : mpz_mul(r, r, num);
63 1970 : mpz_mod(r, r, N);
64 : }
65 2070 : return ret;
66 : }
67 :
68 : int
69 1440 : mod_from_rat_str(mpz_t r, char *str, mpz_t N)
70 : {
71 : mpq_t q;
72 : int ret;
73 :
74 1440 : mpq_init(q);
75 1440 : mpq_set_str(q, str, 10);
76 1440 : ret = mod_from_rat2(r, mpq_numref(q), mpq_denref (q), N);
77 1440 : mpq_clear(q);
78 1440 : return ret;
79 : }
80 :
81 : /* From a curve in Kubert form Y^2+(1-c)*X*Y-b*Y = X^3-b*X^2
82 : to a Weierstrass form y^2 = X^3 + a2 * X^2 + a4 * X + a6
83 : where y = Y+((1-c)*X-b)/2
84 : WE:=[0,(1/4*c^2+1/4-1/2*c-b),0,(1/2*c*b-1/2*b),1/4*b^2]);
85 : We compute:
86 : a2 = 1/4*c^2+1/4-1/2*c-b = ((c-1)/2)^2-b
87 : a4 = 1/2*c*b-1/2*b = b*(c-1)/2
88 : a6 = (b/2)^2
89 : TODO: rewrite this with MediumW, etc.
90 : */
91 : void
92 90 : KW2W246(mpz_t a2, mpz_t a4, mpz_t a6, mpz_t b, mpz_t c, mpz_t n, int compute_a6)
93 : {
94 : /** a4 <- (c-1)/2 **/
95 90 : mpz_sub_si(a4, c, 1);
96 90 : mod_div_2(a4, n);
97 : /** a2 <- a4^2-b **/
98 90 : mpz_mul(a2, a4, a4);
99 90 : mpz_sub(a2, a2, b);
100 90 : mpz_mod(a2, a2, n);
101 : /** a4 <- a4*b **/
102 90 : mpz_mul(a4, a4, b);
103 90 : mpz_mod(a4, a4, n);
104 90 : if(compute_a6 != 0){
105 70 : mpz_set(a6, b);
106 70 : mod_div_2(a6, n);
107 70 : mpz_mul(a6, a6, a6);
108 70 : mpz_mod(a6, a6, n);
109 : }
110 : #if DEBUG_TORSION >= 2
111 : gmp_printf("N:=%Zd;\n", n);
112 : gmp_printf("b:=%Zd;\n", b);
113 : gmp_printf("c:=%Zd;\n", c);
114 : gmp_printf("a2:=%Zd;\n", a2);
115 : gmp_printf("a4:=%Zd;\n", a4);
116 : printf("a6:=RatMod(b^2/4, N);\n");
117 : if(compute_a6 != 0)
118 : gmp_printf("a6:=%Zd;\n", a6);
119 : #endif
120 90 : }
121 :
122 : static int
123 70 : check_weierstrass(mpz_t A, mpz_t B, mpz_t X, mpz_t Y, mpz_t tmp1, mpz_t tmp2,
124 : mpz_t n)
125 : {
126 70 : mpz_mul(tmp1, Y, Y);
127 70 : mpz_mul(tmp2, X, X);
128 70 : mpz_add(tmp2, tmp2, A);
129 70 : mpz_mul(tmp2, tmp2, X);
130 70 : mpz_add(tmp2, tmp2, B);
131 70 : mpz_sub(tmp1, tmp1, tmp2);
132 70 : mpz_mod(tmp1, tmp1, n);
133 70 : return mpz_sgn(tmp1) == 0;
134 : }
135 :
136 : /* Weierstrass (a2, a4, a6) to (A, B)
137 : A = (a4-1/3*a2^2)
138 : B = -1/3*a4*a2 + 2/27*a2^3 + a6
139 : = -1/3*a2*(a4-2/9*a2^2) + a6
140 : X = x+a2/3
141 : Y = y
142 : INPUT: if x0 == NULL, we have no point to translate
143 : if B == NULL, we do not need and we do not compute B
144 : REM: we assume gcd(n, 3) = 1.
145 : */
146 : void
147 70 : MediumWeierstrassToShortWeierstrass(mpz_t A, mpz_t B, mpz_t X, mpz_t Y,
148 : mpz_t a2, mpz_t a4, mpz_t a6,
149 : mpz_t x0, mpz_t y0, mpz_t n)
150 : {
151 : mpz_t tmp1, tmp2, tmp3;
152 :
153 70 : mpz_init(tmp1);
154 70 : mpz_init(tmp2);
155 : /* tmp2 <- a2/3 */
156 70 : mpz_init_set_si(tmp3, 3);
157 70 : mod_from_rat2(tmp2, a2, tmp3, n);
158 70 : if(X != NULL && x0 != NULL){
159 : /* wx0 = x0 + a2/3 */
160 70 : mpz_add(X, tmp2, x0);
161 70 : mpz_mod(X, X, n);
162 : }
163 70 : if(Y != NULL && y0 != NULL){
164 70 : mpz_set(Y, y0);
165 70 : mpz_mod(Y, Y, n);
166 : }
167 : /* A = a4-1/3*a2^2 = a4 - a2 * (a2/3) */
168 : /** tmp1 <- tmp2*a2 = a2^2/3 */
169 70 : mpz_mul(tmp1, a2, tmp2);
170 70 : mpz_mod(tmp1, tmp1, n);
171 70 : mpz_sub(A, a4, tmp1);
172 70 : mpz_mod(A, A, n);
173 70 : if(B != NULL){
174 : /* B = -1/3*a2*(a4-2/9*a2^2) + a6 */
175 : /** B <- 2/9*a2^2 = 2 * (a2^2/3) / 3 **/
176 70 : mod_from_rat2(B, tmp1, tmp3, n);
177 70 : mpz_mul_si(B, B, 2);
178 70 : mpz_sub(B, a4, B);
179 70 : mpz_mul(B, B, tmp2);
180 70 : mpz_sub(B, a6, B);
181 70 : mpz_mod(B, B, n);
182 : }
183 : #if DEBUG_TORSION >= 2
184 : gmp_printf("N:=%Zd;\n", n);
185 : gmp_printf("a2:=%Zd; a4:=%Zd; a6:=%Zd;\n", a2, a4, a6);
186 : gmp_printf("A:=%Zd; B:=%Zd;\n", A, B);
187 : if(X != NULL && x0 != NULL){
188 : gmp_printf("x:=%Zd;\n", x0);
189 : gmp_printf("X:=%Zd;\n", X);
190 : }
191 : if(Y != NULL && y0 != NULL){
192 : gmp_printf("y:=%Zd;\n", Y);
193 : printf("(y^2-x^3-a2*x^2-a4*x-a6) mod N;\n");
194 : printf("(y^2-X^3-A*X-B) mod N;\n");
195 : }
196 : #endif
197 70 : mpz_clear(tmp1);
198 70 : mpz_clear(tmp2);
199 70 : mpz_clear(tmp3);
200 70 : }
201 :
202 : /*
203 : The original Kubert curve E(b, c) is y^2+(1-c)*x*y-b*y = x^3-b*x^2
204 : The medium Weierstrass form is y^2=x^3+a2*x^2+a4*x+a6 with point (x0, y0);
205 : we convert this to short Weierstrass form:
206 : E: Y^2 = X^3 + A * X + B
207 : and point P=(X, Y) on E.
208 : */
209 : void
210 70 : kubert_to_weierstrass(mpz_t A, mpz_t B, mpz_t X, mpz_t Y,
211 : mpz_t b, mpz_t c, mpz_t x0, mpz_t y0, mpz_t n)
212 : {
213 : mpz_t a2, a4, a6;
214 :
215 70 : mpz_init(a2);
216 70 : mpz_init(a4);
217 70 : mpz_init(a6);
218 70 : KW2W246(a2, a4, a6, b, c, n, 1);
219 : /* second conversion */
220 70 : MediumWeierstrassToShortWeierstrass(A, B, X, Y, a2, a4, a6, x0, y0, n);
221 : #if DEBUG_TORSION >= 2
222 : gmp_printf("a2:=%Zd; a4:=%Zd; a6:=%Zd; A:=%Zd; B:=%Zd;\n", a2, a4, a6,A,B);
223 : gmp_printf("X:=%Zd; Y:=%Zd;\n", X, Y);
224 : #endif
225 70 : mpz_clear(a2);
226 70 : mpz_clear(a4);
227 70 : mpz_clear(a6);
228 70 : }
229 :
230 : static int
231 120 : forbidden(char *torsion, int u){
232 120 : if(strcmp(torsion, "Z10") == 0)
233 60 : return u == 1 || u == 2;
234 60 : else if(strcmp(torsion, "Z3xZ3") == 0)
235 60 : return u == 2;
236 0 : return 0;
237 : }
238 :
239 : /* Kubert: put b = c.
240 : SIDE EFFECT: tE[0..nE[ and tP[0..nE[ receive a curve of torsion Z5
241 : and a point on it using parameters [smin..smax[.
242 : OUTPUT: ECM_NO_FACTOR_FOUND or ECM_FACTOR_FOUND_STEP1 if a factor is found.
243 : */
244 : int
245 20 : build_curves_with_torsion_Z5(mpz_t f, mpmod_t n,
246 : ell_curve_t *tE, ell_point_t *tP,
247 : int smin, int smax, int nE)
248 : {
249 : mpz_t A, B, X, Y;
250 20 : int s, ret = ECM_NO_FACTOR_FOUND, nc = 0;
251 : mpz_t x0, y0, c, tmp;
252 :
253 20 : mpz_init(A);
254 20 : mpz_init(B);
255 20 : mpz_init(X);
256 20 : mpz_init(Y);
257 20 : mpz_init(x0);
258 20 : mpz_init(y0);
259 20 : mpz_init(c);
260 20 : mpz_init(tmp);
261 20 : for(s = smin; s < smax; s++){
262 20 : mpz_set_si(x0, s);
263 : /* c:=1/2*x0*(4*x0+1)/(3*x0+1); */
264 : /* y0 <- 2*(3*x0+1) */
265 20 : mpz_mul_si(y0, x0, 3);
266 20 : mpz_add_si(y0, y0, 1);
267 : /* tmp <- x0*(4*x0+1) */
268 20 : mpz_add(tmp, y0, x0);
269 20 : mpz_mul(tmp, tmp, x0);
270 20 : mpz_add(y0, y0, y0);
271 20 : if(mod_from_rat2(c, tmp, y0, n->orig_modulus) == 0){
272 10 : printf("factor found during Z5_init\n");
273 10 : mpz_gcd(f, c, n->orig_modulus);
274 10 : ret = ECM_FACTOR_FOUND_STEP1;
275 10 : break;
276 : }
277 : /* y0:=x0*(x0+1)*(4*x0+1)/4/(3*x0+1) = (x0+1)*c/2 */
278 10 : mpz_add_si(y0, x0, 1);
279 10 : mpz_mul(y0, y0, c);
280 10 : mpz_mod(y0, y0, n->orig_modulus);
281 10 : mod_div_2(y0, n->orig_modulus);
282 : #if DEBUG_TORSION >= 2
283 : gmp_printf("x0:=%Zd;\nc:=%Zd;\ny0:=%Zd;\n", x0, c, y0);
284 : printf("cr:=1/2*x0*(4*x0+1)/(3*x0+1);\n");
285 : #endif
286 : /* P:=WE![x0, y0, 1]; */
287 10 : kubert_to_weierstrass(A, B, X, Y, c, c, x0, y0, n->orig_modulus);
288 10 : if(check_weierstrass(A, B, X, Y, tmp, x0, n->orig_modulus) == 0){
289 0 : printf("#!# check_weierstrass false\n");
290 0 : ret = ECM_ERROR;
291 0 : break;
292 : }
293 10 : ell_curve_init(tE[nc], ECM_EC_TYPE_WEIERSTRASS, ECM_LAW_AFFINE,n);
294 10 : mpz_set(tE[nc]->a4, A);
295 10 : mpz_set(tE[nc]->a6, B);
296 10 : ell_point_init(tP[nc], tE[nc], n);
297 10 : mpz_set(tP[nc]->x, X);
298 10 : mpz_set(tP[nc]->y, Y);
299 10 : nc++;
300 10 : if(nc >= nE)
301 10 : break;
302 : }
303 20 : mpz_clear(A);
304 20 : mpz_clear(B);
305 20 : mpz_clear(X);
306 20 : mpz_clear(Y);
307 20 : mpz_clear(x0);
308 20 : mpz_clear(y0);
309 20 : mpz_clear(c);
310 20 : mpz_clear(tmp);
311 20 : return ret;
312 : }
313 :
314 : /*
315 : E_aux: T^2 = S^3 + A * S + B
316 : => quartic QC: Y^2 = X^4 - 6 * A2 * X^2 + 4 * A1 * X + A0, with
317 : X = (T-A1/2)/(S-A2), Y = -X^2 + 2 * S + A2.
318 : => quartic y^2 = f(x) = a4*x^4+...+a0, where
319 : x = x0+y0/(X-cte), where cte = f'(x0)/4/y0
320 : y = Y/y0*(x-x0)^2 = Y*y0/(X-cte)^2
321 : INPUT: (s, t) is a point on E_aux; (x0, y0) a point on QC.
322 : SIDE EFFECT: x, y contain a point on the elliptic curve.
323 : OUTPUT: 1 if no pb occurred,
324 : 0 if a factor was found and put in f
325 : */
326 : int
327 100 : cubic_to_quartic(mpz_t f, mpz_t n, mpz_t x, mpz_t y,
328 : mpz_t s, mpz_t t, mpz_t A2, mpz_t A1div2,
329 : mpz_t x0, mpz_t y0, mpz_t cte)
330 : {
331 : mpz_t X, Y;
332 100 : int ret = 1;
333 :
334 100 : mpz_init(X);
335 100 : mpz_init(Y);
336 : /* X <- (t-A1/2)/(s-A2) */
337 100 : mpz_sub(x, t, A1div2);
338 100 : mpz_sub(y, s, A2);
339 100 : if(mod_from_rat2(X, x, y, n) == 0){
340 20 : mpz_set(f, X);
341 20 : ret = 0;
342 : }
343 : else{
344 : /* Y <- -X^2 + 2 * s + A2 */
345 80 : mpz_mul(Y, X, X);
346 80 : mpz_sub(Y, A2, Y);
347 80 : mpz_add(Y, Y, s);
348 80 : mpz_add(Y, Y, s);
349 80 : mpz_mod(Y, Y, n);
350 : /* x <- x0+y0/(X-cte) */
351 80 : mpz_sub(X, X, cte);
352 80 : mpz_mod(X, X, n);
353 80 : if(mpz_invert(f, X, n) == 0){
354 10 : mpz_gcd(f, X, n);
355 10 : ret = 0;
356 : }
357 : else{
358 : /* x <- y0/(X-cte) */
359 70 : mpz_mul(x, f, y0);
360 70 : mpz_mod(x, x, n);
361 : /* y <- x/(X-cte) = y0/(X-cte)^2 */
362 70 : mpz_mul(y, x, f);
363 70 : mpz_mod(y, y, n);
364 70 : mpz_mul(y, y, Y);
365 70 : mpz_mod(y, y, n);
366 70 : mpz_add(x, x, x0);
367 70 : mpz_mod(x, x, n);
368 : }
369 : }
370 100 : mpz_clear(X);
371 100 : mpz_clear(Y);
372 100 : return ret;
373 : }
374 :
375 : int
376 160 : build_curves_with_torsion_aux(ell_curve_t Eaux, ell_point_t Paux,
377 : mpz_t A2, mpz_t A1div2, mpz_t x0, mpz_t y0,
378 : mpz_t cte,
379 : char *sa4, char *sa6, char *sPx, char *sPy,
380 : char *sA2, char *sA1div2, char *sx0, char *sy0,
381 : char *scte, mpmod_t n, mpres_t tmp)
382 : {
383 : mpz_t f;
384 :
385 160 : mpz_init(f);
386 160 : mod_from_rat_str(f, sa4, n->orig_modulus);
387 160 : mpres_set_z(tmp, f, n);
388 160 : ell_curve_init_set(Eaux, ECM_EC_TYPE_WEIERSTRASS, ECM_LAW_AFFINE, tmp, n);
389 160 : mod_from_rat_str(f, sa6, n->orig_modulus);
390 160 : mpres_set_z(Eaux->a6, f, n);
391 160 : ell_point_init(Paux, Eaux, n);
392 160 : mod_from_rat_str(f, sPx, n->orig_modulus);
393 160 : mpres_set_z(Paux->x, f, n);
394 160 : mod_from_rat_str(f, sPy, n->orig_modulus);
395 160 : mpres_set_z(Paux->y, f, n);
396 : #if DEBUG_TORSION >= 2
397 : printf("Paux:=");
398 : pt_print(Eaux, Paux, n);
399 : printf(";\n");
400 : #endif
401 160 : mod_from_rat_str(A2, sA2, n->orig_modulus);
402 160 : mod_from_rat_str(A1div2, sA1div2, n->orig_modulus);
403 160 : mod_from_rat_str(x0, sx0, n->orig_modulus);
404 160 : mod_from_rat_str(y0, sy0, n->orig_modulus);
405 160 : mod_from_rat_str(cte, scte, n->orig_modulus);
406 160 : mpz_clear(f);
407 160 : return 1;
408 : }
409 :
410 : /*
411 : SIDE EFFECT: tE[0..nE[ and tP[0..nE[ receive a curve of torsion Z7
412 : and a point on it using parameters [umin..umax[.
413 : OUTPUT: ECM_NO_FACTOR_FOUND or ECM_FACTOR_FOUND_STEP1 if a factor is found.
414 : tE[i], tP[i] are built in raw modular form, not Montgomery form.
415 : REM: we assume gcd(n, 6).
416 : */
417 : int
418 70 : build_curves_with_torsion_Z7(mpz_t fac, mpmod_t n,
419 : ell_curve_t *tE, ell_point_t *tP,
420 : int umin, int umax, int nE)
421 : {
422 70 : int u, ret = ECM_NO_FACTOR_FOUND, nc = 0;
423 : mpz_t A2, A1div2, x0, y0, cte, d, c, b, kx0, ky0, A, B, X, Y;
424 : mpres_t tmp;
425 : ell_curve_t E;
426 : ell_point_t P, Q;
427 :
428 70 : mpz_init(A2);
429 70 : mpz_init(A1div2);
430 70 : mpz_init(cte);
431 70 : mpz_init(x0);
432 70 : mpz_init(y0);
433 70 : mpz_init(A);
434 70 : mpz_init(B);
435 70 : mpz_init(X);
436 70 : mpz_init(Y);
437 : /* Eaux = "1295/48", "-1079/864" */
438 : /* Paux = "2185/12", "-2458" */
439 : /* Y^2 = X^4-1/2*X^2-8*X-1727/16 */
440 70 : mpres_init(tmp, n);
441 70 : build_curves_with_torsion_aux(E, P, A2, A1div2, x0, y0, cte,
442 : "1295/48", "-1079/864",
443 : "2185/12", "-2458",
444 : "1/12", "-1",
445 : "-1", "8", "-7/2",
446 : n, tmp);
447 70 : mpz_init(d);
448 70 : mpz_init(c);
449 70 : mpz_init(b);
450 70 : mpz_init(kx0);
451 70 : mpz_init(ky0);
452 70 : ell_point_init(Q, E, n);
453 70 : mpz_set_si(d, umin-1);
454 70 : if(ell_point_mul(fac, Q, d, P, E, n) == 0){
455 10 : printf("found factor during init of Q in Z7\n");
456 10 : ret = ECM_FACTOR_FOUND_STEP1;
457 : }
458 70 : for(u = umin; (ret != ECM_FACTOR_FOUND_STEP1) && u < umax; u++){
459 : /* update Q */
460 60 : if(ell_point_add(fac, Q, P, Q, E, n) == 0){
461 20 : printf("found factor during update of Q in Z7\n");
462 20 : ret = ECM_FACTOR_FOUND_STEP1;
463 20 : break;
464 : }
465 40 : if(ell_point_is_on_curve(Q, E, n) == 0){
466 0 : printf("#!# Q=[%d]P is not on E\n", u);
467 : // ell_point_print(Q, E, n); printf("\n");
468 0 : ret = ECM_ERROR;
469 0 : break;
470 : }
471 : /* come back to plain (not Montgomery) residues */
472 40 : mpres_get_z(b, Q->x, n);
473 40 : mpres_get_z(c, Q->y, n);
474 : #if DEBUG_TORSION >= 2
475 : gmp_printf("b:=%Zd; c:=%Zd;\n", b, c);
476 : printf("(s, t)[%d]:=", u);
477 : pt_print(E, Q, n);
478 : printf(";\n");
479 : #endif
480 40 : if(cubic_to_quartic(fac, n->orig_modulus, d, ky0, b, c,
481 : A2, A1div2, x0, y0, cte) == 0){
482 10 : printf("found factor in Z7 (cubic_to_quartic)\n");
483 10 : ret = ECM_FACTOR_FOUND_STEP1;
484 10 : break;
485 : }
486 : /* (d, ky0) is a point on y^2 = x^4-18*x^3+13*x^2-28*x+4 */
487 : /* d:=x; */
488 : /* x0:=-2*d; */
489 30 : mpz_mul_si(kx0, d, -2);
490 30 : mpz_mod(kx0, kx0, n->orig_modulus);
491 : /* y0:=d*y/2; */
492 30 : mpz_mul(ky0, ky0, d);
493 30 : mpz_mod(ky0, ky0, n->orig_modulus);
494 30 : mod_div_2(ky0, n->orig_modulus);
495 : /* c:=d^2-d; */
496 30 : mpz_mul(c, d, d);
497 30 : mpz_sub(c, c, d);
498 30 : mpz_mod(c, c, n->orig_modulus);
499 : /* b:=c*d; */
500 30 : mpz_mul(b, c, d);
501 30 : mpz_mod(b, b, n->orig_modulus);
502 : /* to short Weierstrass form */
503 30 : kubert_to_weierstrass(A, B, X, Y, b, c, kx0, ky0, n->orig_modulus);
504 30 : if(check_weierstrass(A, B, X, Y, tmp, x0, n->orig_modulus) == 0){
505 0 : ret = ECM_ERROR;
506 0 : break;
507 : }
508 30 : ell_curve_init(tE[nc], ECM_EC_TYPE_WEIERSTRASS, ECM_LAW_AFFINE,n);
509 30 : mpz_set(tE[nc]->a4, A);
510 30 : mpz_set(tE[nc]->a6, B);
511 30 : ell_point_init(tP[nc], tE[nc], n);
512 30 : mpz_set(tP[nc]->x, X);
513 30 : mpz_set(tP[nc]->y, Y);
514 : #if DEBUG_TORSION >= 2
515 : gmp_printf("E[%d]:=[%Zd];\n", nc, tE[nc]->a4);
516 : gmp_printf("P[%d]:=[%Zd, %Zd, %Zd];\n",
517 : nc, tP[nc]->x, tP[nc]->y, tP[nc]->z);
518 : #endif
519 30 : nc++;
520 30 : if(nc >= nE)
521 30 : break;
522 : }
523 70 : mpz_clear(A2);
524 70 : mpz_clear(A1div2);
525 70 : mpz_clear(x0);
526 70 : mpz_clear(y0);
527 70 : mpz_clear(cte);
528 70 : mpz_clear(A);
529 70 : mpz_clear(B);
530 70 : mpz_clear(X);
531 70 : mpz_clear(Y);
532 70 : ell_point_clear(P, E, n);
533 70 : ell_point_clear(Q, E, n);
534 70 : ell_curve_clear(E, n);
535 70 : mpz_clear(d);
536 70 : mpz_clear(c);
537 70 : mpz_clear(b);
538 70 : mpz_clear(kx0);
539 70 : mpz_clear(ky0);
540 70 : mpres_clear(tmp, n);
541 70 : return ret;
542 : }
543 :
544 : /*
545 : SIDE EFFECT: tE[0..nE[ and tP[0..nE[ receive a curve of torsion Z9
546 : and a point on it using parameters [umin..umax[.
547 : OUTPUT: ECM_NO_FACTOR_FOUND or ECM_FACTOR_FOUND_STEP1 if a factor is found.
548 : tE[i], tP[i] are built in raw modular form, not Montgomery form.
549 : REM: we assume gcd(n, 6).
550 : */
551 : int
552 40 : build_curves_with_torsion_Z9(mpz_t fac, mpmod_t n, ell_curve_t *tE,
553 : ell_point_t *tP, int umin, int umax, int nE)
554 : {
555 40 : int u, ret = ECM_NO_FACTOR_FOUND, nc = 0;
556 : mpz_t A2, A1div2, x0, y0, cte, d, c, b, kx0, ky0, A, B, X, Y, f;
557 : mpres_t tmp;
558 : ell_curve_t E;
559 : ell_point_t P, Q;
560 :
561 40 : mpz_init(A2);
562 40 : mpz_init(A1div2);
563 40 : mpz_init(cte);
564 40 : mpz_init(x0);
565 40 : mpz_init(y0);
566 40 : mpz_init(A);
567 40 : mpz_init(B);
568 40 : mpz_init(X);
569 40 : mpz_init(Y);
570 : /* Eaux = [-9, 9] */
571 : /* Paux = [1, 1, 1] */
572 : /* Y^2 = X^4-24*X-36 */
573 40 : mpres_init(tmp, n);
574 40 : build_curves_with_torsion_aux(E, P, A2, A1div2, x0, y0, cte,
575 : "-9", "9", "1", "1",
576 : "0", "3",
577 : "2", "3", "0",
578 : n, tmp);
579 40 : mpz_init(f);
580 40 : mpz_init(d);
581 40 : mpz_init(c);
582 40 : mpz_init(b);
583 40 : mpz_init(kx0);
584 40 : mpz_init(ky0);
585 40 : ell_point_init(Q, E, n);
586 40 : mpz_set_si(d, umin-1);
587 40 : if(ell_point_mul(fac, Q, d, P, E, n) == 0){
588 10 : printf("found factor during init of Q in Z9\n");
589 10 : ret = ECM_FACTOR_FOUND_STEP1;
590 : }
591 40 : for(u = umin; (ret != ECM_FACTOR_FOUND_STEP1) && u < umax; u++){
592 : /* update Q */
593 30 : if(ell_point_add(fac, Q, P, Q, E, n) == 0){
594 10 : printf("found factor during update of Q in Z9\n");
595 10 : ret = ECM_FACTOR_FOUND_STEP1;
596 10 : break;
597 : }
598 : #if DEBUG_TORSION >= 2
599 : printf("(s, t)[%d]:=", u);
600 : pt_print(E, Q, n);
601 : printf(";\n");
602 : #endif
603 20 : if(ell_point_is_on_curve(Q, E, n) == 0){
604 0 : printf("#!# Q=[%d]P is not on E\n", u);
605 0 : ret = ECM_ERROR;
606 0 : break;
607 : }
608 20 : mpres_get_z(b, Q->x, n);
609 20 : mpres_get_z(c, Q->y, n);
610 20 : if(cubic_to_quartic(fac, n->orig_modulus, f, ky0, b, c,
611 : A2, A1div2, x0, y0, cte) == 0){
612 10 : printf("found factor in Z9 (cubic_2_quartic)\n");
613 10 : ret = ECM_FACTOR_FOUND_STEP1;
614 10 : break;
615 : }
616 : /* f:=x; */
617 : /* d:=f*(f-1)+1; */
618 10 : mpz_sub_si(d, f, 1);
619 10 : mpz_mul(d, d, f);
620 10 : mpz_add_si(d, d, 1);
621 10 : mpz_mod(d, d, n->orig_modulus);
622 : /* c:=f*(d-1); */
623 10 : mpz_sub_si(c, d, 1);
624 10 : mpz_mul(c, c, f);
625 10 : mpz_mod(c, c, n->orig_modulus);
626 : /* kx0:=(2*f-1)*f^2; */
627 : /** b <- f^2 **/
628 10 : mpz_mul(b, f, f);
629 10 : mpz_mod(b, b, n->orig_modulus);
630 10 : mpz_mul_si(kx0, f, 2);
631 10 : mpz_sub_si(kx0, kx0, 1);
632 10 : mpz_mul(kx0, kx0, b);
633 10 : mpz_mod(kx0, kx0, n->orig_modulus);
634 : /* ky0:=y*f^4/2; */
635 : /** b <- b^2 = f^4 **/
636 10 : mpz_mul(b, b, b);
637 10 : mpz_mul(ky0, ky0, b);
638 10 : mpz_mod(ky0, ky0, n->orig_modulus);
639 10 : mod_div_2(ky0, n->orig_modulus);
640 : /* b:=c*d; */
641 10 : mpz_mul(b, c, d);
642 10 : mpz_mod(b, b, n->orig_modulus);
643 : #if DEBUG_TORSION >= 2
644 : gmp_printf("f=%Zd d=%Zd c=%Zd b=%Zd\n", f, d, c, b);
645 : gmp_printf("kx0=%Zd ky0=%Zd\n", kx0, ky0);
646 : #endif
647 : /* to short Weierstrass form */
648 10 : kubert_to_weierstrass(A, B, X, Y, b, c, kx0, ky0, n->orig_modulus);
649 10 : if(check_weierstrass(A, B, X, Y, tmp, x0, n->orig_modulus) == 0){
650 0 : ret = ECM_ERROR;
651 0 : break;
652 : }
653 10 : ell_curve_init(tE[nc], ECM_EC_TYPE_WEIERSTRASS, ECM_LAW_AFFINE, n);
654 10 : mpz_set(tE[nc]->a4, A);
655 10 : mpz_set(tE[nc]->a6, B);
656 10 : ell_point_init(tP[nc], tE[nc], n);
657 10 : mpz_set(tP[nc]->x, X);
658 10 : mpz_set(tP[nc]->y, Y);
659 10 : nc++;
660 10 : if(nc >= nE)
661 10 : break;
662 : }
663 40 : mpz_clear(A);
664 40 : mpz_clear(B);
665 40 : mpz_clear(X);
666 40 : mpz_clear(Y);
667 40 : mpz_clear(A2);
668 40 : mpz_clear(A1div2);
669 40 : mpz_clear(x0);
670 40 : mpz_clear(y0);
671 40 : mpz_clear(cte);
672 40 : ell_point_clear(P, E, n);
673 40 : ell_point_clear(Q, E, n);
674 40 : mpz_clear(f);
675 40 : mpz_clear(d);
676 40 : mpz_clear(c);
677 40 : mpz_clear(b);
678 40 : mpz_clear(kx0);
679 40 : mpz_clear(ky0);
680 40 : mpres_clear(tmp, n);
681 40 : return ret;
682 : }
683 :
684 : int
685 50 : build_curves_with_torsion_Z10(mpz_t fac, mpmod_t n, ell_curve_t *tE,
686 : ell_point_t *tP, int umin, int umax, int nE)
687 : {
688 50 : int u, ret = ECM_NO_FACTOR_FOUND, nc = 0;
689 : mpz_t A2, A1div2, x0, y0, cte, d, c, b, kx0, ky0, A, B, X, Y;
690 : mpz_t f;
691 : mpres_t tmp;
692 : ell_curve_t E;
693 : ell_point_t P, Q;
694 :
695 50 : mpz_init(A2);
696 50 : mpz_init(A1div2);
697 50 : mpz_init(cte);
698 50 : mpz_init(x0);
699 50 : mpz_init(y0);
700 50 : mpz_init(A);
701 50 : mpz_init(B);
702 50 : mpz_init(X);
703 50 : mpz_init(Y);
704 : /* Eaux = [2/3, -53/108] */
705 : /* Paux = [2/3, 1/2, 1] */
706 : /* Y^2 = X^4-4*X^2-4*X-4 */
707 50 : mpres_init(tmp, n);
708 50 : build_curves_with_torsion_aux(E, P, A2, A1div2, x0, y0, cte,
709 : "2/3", "-53/108", "2/3", "1/2",
710 : "2/3", "-1/2",
711 : "0", "1", "-2",
712 : n, tmp);
713 50 : mpz_init(f);
714 50 : mpz_init(d);
715 50 : mpz_init(c);
716 50 : mpz_init(b);
717 50 : mpz_init(kx0);
718 50 : mpz_init(ky0);
719 50 : ell_point_init(Q, E, n);
720 60 : for(u = umin; u < umax; u++){
721 60 : if(forbidden("Z10", u))
722 10 : continue;
723 : /* update Qaux */
724 50 : mpz_set_si(d, u);
725 50 : if(ell_point_mul(fac, Q, d, P, E, n) == 0){
726 10 : printf("found factor during update of Q in Z10\n");
727 10 : ret = ECM_FACTOR_FOUND_STEP1;
728 10 : break;
729 : }
730 : #if DEBUG_TORSION >= 2
731 : printf("(s, t)[%d]:=", u);
732 : pt_print(E, Q, n);
733 : printf(";\n");
734 : #endif
735 40 : if(ell_point_is_on_curve(Q, E, n) == 0){
736 0 : printf("#!# Q=[%d]P is not on E\n", u);
737 0 : ret = ECM_ERROR;
738 0 : break;
739 : }
740 40 : mpres_get_z(b, Q->x, n);
741 40 : mpres_get_z(c, Q->y, n);
742 40 : if(cubic_to_quartic(fac, n->orig_modulus, f, ky0, b, c,
743 : A2, A1div2, x0, y0, cte) == 0){
744 10 : printf("found factor in Z10 (cubic_2_quartic)\n");
745 10 : ret = ECM_FACTOR_FOUND_STEP1;
746 10 : break;
747 : }
748 : /* f:=x; */
749 : /* d:=f^2/(f-(f-1)^2) = -f^2/(f^2-3*f+1) */
750 : /** b <- -f^2 **/
751 30 : mpz_mul(b, f, f);
752 30 : mpz_neg(b, b);
753 30 : mpz_mod(b, b, n->orig_modulus);
754 : /* c = f^2-3*f+1 = f*(f-3)+1 */
755 30 : mpz_sub_si(c, f, 3);
756 30 : mpz_mul(c, c, f);
757 30 : mpz_add_si(c, c, 1);
758 30 : mpz_mod(c, c, n->orig_modulus);
759 30 : if(mod_from_rat2(d, b, c, n->orig_modulus) == 0){
760 10 : printf("inverse found in Z10 (d)\n");
761 10 : mpz_set(fac, d);
762 10 : ret = ECM_FACTOR_FOUND_STEP1;
763 10 : break;
764 : }
765 : /* c:=f*(d-1); */
766 20 : mpz_sub_si(c, d, 1);
767 20 : mpz_mul(c, c, f);
768 20 : mpz_mod(c, c, n->orig_modulus);
769 : /* ky0:=y*f^4/(f^2-3*f+1)^2/2; = num/den */
770 : /* it seems that ky0 = y*d^2/2 */
771 20 : mpz_mul(b, ky0, d);
772 20 : mpz_mul(b, b, d);
773 20 : mpz_mod(b, b, n->orig_modulus);
774 20 : mpz_set_si(fac, 2);
775 20 : mod_from_rat2(ky0, b, fac, n->orig_modulus);
776 : /* kx0:=-f*d; */
777 20 : mpz_mul(kx0, f, d);
778 20 : mpz_neg(kx0, kx0);
779 20 : mpz_mod(kx0, kx0, n->orig_modulus);
780 : /* b:=c*d; */
781 20 : mpz_mul(b, c, d);
782 20 : mpz_mod(b, b, n->orig_modulus);
783 : #if DEBUG_TORSION >= 2
784 : gmp_printf("f:=%Zd; d:=%Zd; c:=%Zd; b:=%Zd;\n", f, d, c, b);
785 : gmp_printf("kx0:=%Zd; ky0:=%Zd;\n", kx0, ky0);
786 : #endif
787 : /* to short Weierstrass form */
788 20 : kubert_to_weierstrass(A, B, X, Y, b, c, kx0, ky0, n->orig_modulus);
789 20 : if(check_weierstrass(A, B, X, Y, tmp, x0, n->orig_modulus) == 0){
790 0 : ret = ECM_ERROR;
791 0 : break;
792 : }
793 20 : ell_curve_init(tE[nc], ECM_EC_TYPE_WEIERSTRASS, ECM_LAW_AFFINE, n);
794 20 : mpz_set(tE[nc]->a4, A);
795 20 : mpz_set(tE[nc]->a6, B);
796 20 : ell_point_init(tP[nc], tE[nc], n);
797 20 : mpz_set(tP[nc]->x, X);
798 20 : mpz_set(tP[nc]->y, Y);
799 20 : nc++;
800 20 : if(nc >= nE)
801 20 : break;
802 : }
803 : #if DEBUG_TORSION >= 2
804 : if(ret != ECM_ERROR && nc > 0){
805 : printf("Curves built\n");
806 : pt_many_print(tE, tP, nE, n);
807 : }
808 : #endif
809 50 : mpz_clear(A);
810 50 : mpz_clear(B);
811 50 : mpz_clear(X);
812 50 : mpz_clear(Y);
813 50 : mpz_clear(A2);
814 50 : mpz_clear(A1div2);
815 50 : mpz_clear(x0);
816 50 : mpz_clear(y0);
817 50 : mpz_clear(cte);
818 50 : ell_point_clear(P, E, n);
819 50 : ell_point_clear(Q, E, n);
820 50 : ell_curve_clear(E, n);
821 50 : mpres_clear(tmp, n);
822 50 : mpz_clear(d);
823 50 : mpz_clear(c);
824 50 : mpz_clear(b);
825 50 : mpz_clear(kx0);
826 50 : mpz_clear(ky0);
827 50 : mpz_clear(f);
828 50 : return ret;
829 : }
830 :
831 : /* Warning: b and a have the Montgomery meaning in this function.
832 : All tE[i] will be in Montgomery form: B*Y^2 = X^3 + A * X^2 + X.
833 : */
834 : int
835 80 : build_curves_with_torsion_Z2xZ8(mpz_t fac, mpmod_t n,
836 : ell_curve_t *tE, ell_point_t *tP,
837 : int umin, int umax, int nE)
838 : {
839 80 : int u, nc = 0, ret = ECM_NO_FACTOR_FOUND;
840 : mpz_t tmp, a, b, alpha, beta, c, d, kx0, ky0, wx0, mb;
841 : mpres_t tmp2;
842 : ell_curve_t E;
843 : ell_point_t P, Q;
844 :
845 80 : mpz_init(alpha);
846 80 : mpz_init(beta);
847 80 : mpz_init(tmp);
848 80 : mpz_init(a);
849 80 : mpz_init(b);
850 80 : mpz_init(c);
851 80 : mpz_init(d);
852 80 : mpz_init(kx0);
853 80 : mpz_init(ky0);
854 80 : mpz_init(wx0);
855 80 : mpz_init(mb);
856 :
857 : /* Eaux = [-8, -32] */
858 : /* Paux = [12, 40, 1] */
859 80 : mpres_init(tmp2, n);
860 80 : mpz_set_str(fac, "-8", 10);
861 80 : mpres_set_z(tmp2, fac, n);
862 80 : ell_curve_init_set(E, ECM_EC_TYPE_WEIERSTRASS, ECM_LAW_AFFINE, tmp2, n);
863 80 : ell_point_init(P, E, n);
864 80 : mpz_set_str(fac, "12", 10);
865 80 : mpres_set_z(P->x, fac, n);
866 80 : mpz_set_str(fac, "40", 10);
867 80 : mpres_set_z(P->y, fac, n);
868 80 : mpz_set_ui(P->z, 1);
869 :
870 80 : ell_point_init(Q, E, n);
871 80 : mpz_set_si(d, umin-1);
872 80 : if(ell_point_mul(fac, Q, d, P, E, n) == 0){
873 10 : printf("found factor during init of Q in Z2xZ8\n");
874 10 : ret = ECM_FACTOR_FOUND_STEP1;
875 : }
876 80 : for(u = umin; (ret != ECM_FACTOR_FOUND_STEP1) && u < umax; u++){
877 : /* update Q */
878 70 : if(ell_point_add(fac, Q, P, Q, E, n) == 0){
879 10 : printf("found factor during update of Q in Z2xZ8\n");
880 10 : ret = ECM_FACTOR_FOUND_STEP1;
881 10 : break;
882 : }
883 : #if DEBUG_TORSION >= 2
884 : printf("(s, t)[%d]:=", u);
885 : pt_print(E, Q, n);
886 : printf(";\n");
887 : #endif
888 : /* beta <- (y+25)/(x-9) */
889 60 : mpres_get_z(a, Q->x, n);
890 60 : mpres_get_z(b, Q->y, n);
891 60 : mpz_mod(wx0, a, n->orig_modulus);
892 60 : mpz_sub_si(a, a, 9);
893 60 : mpz_mod(a, a, n->orig_modulus);
894 60 : mpz_add_si(b, b, 25);
895 60 : mpz_mod(b, b, n->orig_modulus);
896 60 : if(mod_from_rat2(beta, b, a, n->orig_modulus) == 0){
897 10 : printf("found factor in Z2xZ8 (beta)\n");
898 10 : mpz_set(fac, beta);
899 10 : ret = ECM_FACTOR_FOUND_STEP1;
900 10 : break;
901 : }
902 50 : mpz_add_si(tmp, beta, 1);
903 50 : mpz_mod(tmp, tmp, n->orig_modulus);
904 : /* alpha <- 1/(beta+1) */
905 50 : if(mpz_invert(alpha, tmp, n->orig_modulus) == 0){
906 10 : printf("found factor in Z2xZ8 (alpha)\n");
907 10 : mpz_gcd(fac, tmp, n->orig_modulus);
908 10 : ret = ECM_FACTOR_FOUND_STEP1;
909 10 : break;
910 : }
911 : /** d <- 8*alpha^2-1;
912 : d = -(beta^2+2*beta-7)/(beta+1)^2
913 : **/
914 40 : mpz_mul(d, alpha, alpha);
915 40 : mpz_mul_si(d, d, 8);
916 40 : mpz_sub_si(d, d, 1);
917 40 : mpz_mod(d, d, n->orig_modulus);
918 : /* d:=2*alpha*(4*alpha+1)/d; */
919 40 : mpz_mul_si(c, alpha, 4);
920 40 : mpz_add_si(c, c, 1);
921 40 : mpz_mul(c, c, alpha);
922 40 : mpz_mul_si(c, c, 2);
923 40 : mpz_mod(c, c, n->orig_modulus);
924 40 : if(mod_from_rat2(fac, c, d, n->orig_modulus) == 0){
925 : // the only possibility is d = 0 mod p or 8*alpha^2-1 = 0 mod p
926 10 : printf("found factor in Z2xZ8 (d)\n");
927 10 : ret = ECM_FACTOR_FOUND_STEP1;
928 10 : break;
929 : }
930 30 : mpz_set(d, fac);
931 : /* c:=(2*d-1)*(d-1)/d;*/
932 30 : mpz_sub_si(fac, d, 1);
933 : /** kx0 <- 2*d-1 **/
934 30 : mpz_mul_si(kx0, d, 2);
935 30 : mpz_sub_si(kx0, kx0, 1);
936 30 : mpz_mul(fac, fac, kx0);
937 30 : mpz_mod(fac, fac, n->orig_modulus);
938 30 : if(mod_from_rat2(c, fac, d, n->orig_modulus) == 0){
939 : // this is possible only if d = 0 mod p or
940 : // 2*alpha*(4*alpha+1) = 0 mod p
941 10 : printf("found factor in Z2xZ8 (d2)\n");
942 10 : mpz_set(fac, c);
943 10 : ret = ECM_FACTOR_FOUND_STEP1;
944 10 : break;
945 : }
946 : /* b = c*d */
947 20 : mpz_mul(b, c, d);
948 20 : mpz_mod(b, b, n->orig_modulus);
949 : /* kx0:=-(2*d-1)/4;*/
950 20 : mod_div_2(kx0, n->orig_modulus);
951 20 : mod_div_2(kx0, n->orig_modulus);
952 20 : mpz_mul_si(kx0, kx0, -1);
953 20 : mpz_mod(kx0, kx0, n->orig_modulus);
954 : /* ky0:=(c/8)*(-beta^2+2*uP[1]+9); */
955 20 : mpz_mul(fac, beta, beta);
956 20 : mpz_set(a, wx0);
957 20 : mpz_sub(fac, a, fac);
958 20 : mpz_add(fac, fac, a);
959 20 : mpz_add_si(fac, fac, 9);
960 20 : mpz_mul(fac, fac, c);
961 20 : mpz_mod(fac, fac, n->orig_modulus);
962 20 : mod_div_2(fac, n->orig_modulus);
963 20 : mod_div_2(fac, n->orig_modulus);
964 20 : mod_div_2(fac, n->orig_modulus);
965 : /* ky0:=ky0/(beta^2+2*beta-7); */
966 20 : mpz_add_si(tmp, beta, 2);
967 20 : mpz_mul(tmp, tmp, beta);
968 20 : mpz_sub_si(tmp, tmp, 7);
969 20 : mpz_mod(tmp, tmp, n->orig_modulus);
970 : /* as proven above, we cannot have tmp non invertible at that point */
971 20 : mod_from_rat2(ky0, fac, tmp, n->orig_modulus);
972 20 : KW2W246(fac, a, NULL, b, c, n->orig_modulus, 0);
973 : #if DEBUG_TORSION >= 2
974 : gmp_printf("kwx0:=%Zd;\n", kx0);
975 : gmp_printf("kwy0:=%Zd;\n", ky0);
976 : printf("(kwy0^2-(kwx0^3+a2*kwx0^2+a4*kwx0+a6)) mod N;\n");
977 : #endif
978 : /* wx0:=kx0+a2/3; */
979 20 : mpz_set_si(tmp, 3);
980 20 : mod_from_rat2(wx0, fac, tmp, n->orig_modulus);
981 20 : mpz_add(wx0, wx0, kx0);
982 20 : mpz_mod(wx0, wx0, n->orig_modulus);
983 : /* mb:=-1/(d-1)^2; */
984 20 : mpz_sub_si(tmp, d, 1);
985 20 : mpz_mul(tmp, tmp, tmp);
986 20 : mpz_mod(tmp, tmp, n->orig_modulus);
987 20 : mpz_neg(tmp, tmp);
988 20 : if(mpz_invert(mb, tmp, n->orig_modulus) == 0){
989 10 : printf("found factor in Z2xZ8 (mb)\n");
990 10 : mpz_gcd(fac, tmp, n->orig_modulus);
991 10 : ret = ECM_FACTOR_FOUND_STEP1;
992 10 : break;
993 : }
994 : /* ma:=-1/4*(8*d^4-16*d^3+16*d^2-8*d+1)/(d-1)^2/d^2;
995 : :=mb*(8*d^4-16*d^3+16*d^2-8*d+1)/(4*d^2)
996 : */
997 10 : mpz_set_si(fac, 8); /* num */
998 10 : mpz_mul(fac, fac, d); mpz_add_si(fac, fac, -16);
999 10 : mpz_mul(fac, fac, d); mpz_add_si(fac, fac, 16);
1000 10 : mpz_mul(fac, fac, d); mpz_add_si(fac, fac, -8);
1001 10 : mpz_mul(fac, fac, d); mpz_add_si(fac, fac, 1);
1002 : #if 0
1003 : mpz_sub_si(tmp, d, 1); /* den */
1004 : mpz_mul(tmp, tmp, d);
1005 : mpz_mul(tmp, tmp, tmp);
1006 : mpz_mul_si(tmp, tmp, -4);
1007 : mpz_mod(tmp, tmp, n->orig_modulus);
1008 : #else
1009 10 : mpz_mul(fac, fac, mb);
1010 : /* one day, we could save 1/d computation again */
1011 10 : mpz_mul(tmp, d, d);
1012 10 : mpz_mul_si(tmp, tmp, 4);
1013 : #endif
1014 : /* to Montgomery form */
1015 10 : ell_curve_init(tE[nc], ECM_EC_TYPE_MONTGOMERY, ECM_LAW_HOMOGENEOUS,n);
1016 10 : ell_point_init(tP[nc], tE[nc], n);
1017 : /* this cannot yield a factor, since d is invertible at that point */
1018 10 : mod_from_rat2(tE[nc]->a2, fac, tmp, n->orig_modulus);
1019 : /* not really needed, but useful for debug */
1020 10 : mpz_set_ui(tE[nc]->a4, 1);
1021 10 : mpz_set_ui(tE[nc]->a6, 0);
1022 : /* mx:=mb*wx0-ma/3; */
1023 10 : mpz_mul(fac, mb, wx0);
1024 10 : mpz_set_si(tmp, 3);
1025 10 : mod_from_rat2(tP[nc]->x, tE[nc]->a2, tmp, n->orig_modulus);
1026 10 : mpz_sub(tP[nc]->x, fac, tP[nc]->x);
1027 10 : mpz_mod(tP[nc]->x, tP[nc]->x, n->orig_modulus);
1028 : /* my:=mb*ky0; */
1029 : #if DEBUG_TORSION >= 2
1030 : gmp_printf("N:=%Zd;\n", n->orig_modulus);
1031 : gmp_printf("ma:=%Zd;\n", tE[nc]->a2);
1032 : gmp_printf("mb:=%Zd;\n", mb);
1033 : gmp_printf("kx0:=%Zd;\n", kx0);
1034 : gmp_printf("ky0:=%Zd;\n", ky0);
1035 : gmp_printf("mx0:=%Zd;\n", tP[nc]->x);
1036 : mpz_mul(tmp, mb, ky0);
1037 : mpz_mod(tmp, tmp, n->orig_modulus);
1038 : gmp_printf("my0:=%Zd;\n", tmp);
1039 : printf("chk:=(mb*my0^2-mx0^3-ma*mx0^2-mx0) mod N;\n");
1040 : #endif
1041 10 : nc++;
1042 10 : if(nc >= nE)
1043 10 : break;
1044 : }
1045 : #if DEBUG_TORSION >= 2
1046 : printf("Curves built\n");
1047 : pt_many_print(tE, tP, nE, n);
1048 : #endif
1049 80 : ell_point_clear(P, E, n);
1050 80 : ell_point_clear(Q, E, n);
1051 80 : ell_curve_clear(E, n);
1052 80 : mpz_clear(mb);
1053 80 : mpz_clear(tmp);
1054 80 : mpz_clear(a);
1055 80 : mpz_clear(b);
1056 80 : mpz_clear(c);
1057 80 : mpz_clear(d);
1058 80 : mpz_clear(alpha);
1059 80 : mpz_clear(beta);
1060 80 : mpz_init(kx0);
1061 80 : mpz_init(ky0);
1062 80 : mpz_init(wx0);
1063 80 : mpres_clear(tmp2, n);
1064 80 : return ret;
1065 : }
1066 :
1067 : /* Z3xZ3 over Q(sqrt(-3)). Interesting if we know that p | N is s.t.
1068 : p = 1 mod 3.
1069 : Source: Dujella and Najman, arxiv:1201.0266v1
1070 : A more simpler and more efficient stuff, using Hessian form. */
1071 : int
1072 50 : build_curves_with_torsion_Z3xZ3(mpz_t f, mpmod_t n,
1073 : ell_curve_t *tE, ell_point_t *tP,
1074 : int umin, int umax, int nE)
1075 : {
1076 50 : int u, nc = 0, ret = ECM_NO_FACTOR_FOUND;
1077 : mpz_t u0, v0, D, num, den;
1078 :
1079 50 : mpz_init(u0);
1080 50 : mpz_init(num);
1081 50 : mpz_init(den);
1082 50 : mpz_init(D);
1083 50 : mpz_init_set_si(v0, umin-1); /* to prevent u0 = v0 */
1084 60 : for(u = umin; u < umax; u++){
1085 60 : if(forbidden("Z3xZ3", u))
1086 10 : continue;
1087 50 : mpz_set_si(u0, u);
1088 : /* D:=RatMod((u0^3+v0^3+1)/(3*u0*v0), N); */
1089 50 : mpz_mul(num, u0, u0);
1090 50 : mpz_mul(num, num, u0);
1091 50 : mpz_mul(den, v0, v0);
1092 50 : mpz_mul(den, den, v0);
1093 50 : mpz_add(num, num, den);
1094 50 : mpz_add_si(num, num, 1);
1095 50 : mpz_mod(num, num, n->orig_modulus);
1096 50 : if(mpz_sgn(num) == 0)
1097 0 : continue;
1098 :
1099 50 : mpz_mul(den, u0, v0);
1100 50 : mpz_mul_si(den, den, 3);
1101 50 : mpz_mod(den, den, n->orig_modulus);
1102 :
1103 50 : if(mod_from_rat2(D, num, den, n->orig_modulus) == 0){
1104 10 : printf("found factor in Z3xZ3 (D)\n");
1105 10 : mpz_set(f, D);
1106 10 : ret = ECM_FACTOR_FOUND_STEP1;
1107 10 : break;
1108 : }
1109 40 : mpz_mul(num, D, D);
1110 40 : mpz_mul(num, num, D);
1111 40 : mpz_mod(num, num, n->orig_modulus);
1112 40 : if(mpz_cmp_ui(num, 1) == 0){
1113 10 : printf("D^3=1 => singluar curve\n");
1114 10 : ret = ECM_ERROR;
1115 10 : break;
1116 : }
1117 30 : ell_curve_init_set(tE[nc],ECM_EC_TYPE_HESSIAN,ECM_LAW_HOMOGENEOUS,D,n);
1118 30 : ell_point_init(tP[nc], tE[nc], n);
1119 30 : mpz_set(tP[nc]->x, u0);
1120 30 : mpz_set(tP[nc]->y, v0);
1121 30 : mpz_set_ui(tP[nc]->z, 1);
1122 30 : nc++;
1123 30 : if(nc >= nE)
1124 30 : break;
1125 : }
1126 50 : mpz_clear(u0);
1127 50 : mpz_clear(v0);
1128 50 : mpz_clear(D);
1129 50 : mpz_clear(num);
1130 50 : mpz_clear(den);
1131 50 : return ret;
1132 : }
1133 :
1134 : /* For a small price, add a 2-torsion point, also over Q(sqrt(-3)). */
1135 : int
1136 30 : build_curves_with_torsion_Z3xZ6(mpz_t f, mpmod_t n,
1137 : ell_curve_t *tE, ell_point_t *tP,
1138 : int umin, int umax, int nE)
1139 : {
1140 30 : int u, nc = 0, ret = ECM_NO_FACTOR_FOUND;
1141 : ell_curve_t E;
1142 : ell_point_t P, Q;
1143 : mpres_t tmp, num, den, tk, sk;
1144 : mpz_t t;
1145 :
1146 30 : mpz_init(t);
1147 30 : mpz_init(num);
1148 30 : mpz_init(den);
1149 30 : mpz_init(tk);
1150 30 : mpz_init(sk);
1151 : /* Eaux:=EllipticCurve([0, -4]); */
1152 : /* Paux:=Eaux![2, 2, 1]; */
1153 30 : mpres_init(tmp, n);
1154 30 : mpres_set_ui(tmp, 0, n);
1155 30 : ell_curve_init_set(E, ECM_EC_TYPE_WEIERSTRASS, ECM_LAW_AFFINE, tmp, n);
1156 30 : ell_point_init(P, E, n);
1157 30 : mpz_set_str(f, "2", 10);
1158 30 : mpres_set_z(P->x, f, n);
1159 30 : mpz_set_str(f, "2", 10);
1160 30 : mpres_set_z(P->y, f, n);
1161 30 : mpz_set_ui(P->z, 1);
1162 :
1163 30 : ell_point_init(Q, E, n);
1164 30 : for(u = umin; u < umax; u++){
1165 : /* update Qaux */
1166 30 : mpz_set_si(f, u);
1167 30 : if(ell_point_mul(f, Q, f, P, E, n) == 0){
1168 10 : printf("found factor in Z3xZ6 (update of Q)\n");
1169 10 : ret = ECM_FACTOR_FOUND_STEP1;
1170 10 : break;
1171 : }
1172 : #if DEBUG_TORSION >= 2
1173 : printf("(s, t)[%d]:=", u);
1174 : pt_print(E, Q, n);
1175 : printf(";\n");
1176 : #endif
1177 20 : mpres_get_z(tk, Q->x, n);
1178 20 : mpres_get_z(sk, Q->y, n);
1179 : #if 0 /* useless in affine form? */
1180 : mpres_get_z(t, Q->z, n);
1181 : if(mpz_invert(f, t, n->orig_modulus) == 0){
1182 : printf("found factor in Z3xZ6 (normalization)\n");
1183 : mpz_gcd(f, t, n->orig_modulus);
1184 : break;
1185 : }
1186 : mpz_mul(tk, tk, f);
1187 : mpz_mod(tk, tk, n->orig_modulus);
1188 : mpz_mul(sk, sk, f);
1189 : mpz_mod(sk, sk, n->orig_modulus);
1190 : #endif
1191 : /* t:=RatMod(-tk/2, N); */
1192 20 : mpz_mul_si(t, tk, -1);
1193 20 : mod_div_2(t, n->orig_modulus);
1194 : /* D:=RatMod((2*t^3+1)/3/t^2, N); */
1195 20 : mpz_mul(den, t, t);
1196 20 : mpz_mod(den, den, n->orig_modulus);
1197 20 : mpz_mul(num, den, t);
1198 20 : mpz_mul_si(num, num, 2);
1199 20 : mpz_add_si(num, num, 1);
1200 20 : mpz_mod(num, num, n->orig_modulus);
1201 20 : mpz_mul_si(den, den, 3);
1202 20 : mpz_mod(den, den, n->orig_modulus);
1203 20 : ell_curve_init(tE[nc], ECM_EC_TYPE_HESSIAN, ECM_LAW_HOMOGENEOUS, n);
1204 20 : ell_point_init(tP[nc], tE[nc], n);
1205 20 : if(mod_from_rat2(tE[nc]->a4, num, den, n->orig_modulus) == 0){
1206 : /* only if t = 0, which seems hard */
1207 10 : printf("found factor in Z3xZ6 (D)\n");
1208 10 : mpz_set(f, tE[nc]->a4);
1209 10 : ret = ECM_FACTOR_FOUND_STEP1;
1210 10 : break;
1211 : }
1212 : #if DEBUG_TORSION >= 1
1213 : gmp_printf("D%d:=%Zd;\n", nc, tE[nc]->a4);
1214 : #endif
1215 : /* u0:=RatMod(sk/tk, N);
1216 : if tk was not invertible, it would have been caught before
1217 : */
1218 10 : mod_from_rat2(tP[nc]->x, sk, tk, n->orig_modulus);
1219 : /* v0:=-1; */
1220 10 : mpz_sub_si(tP[nc]->y, n->orig_modulus, 1);
1221 10 : mpz_set_ui(tP[nc]->z, 1);
1222 10 : nc++;
1223 10 : if(nc >= nE)
1224 10 : break;
1225 : }
1226 30 : mpz_clear(t);
1227 30 : mpz_clear(num);
1228 30 : mpz_clear(den);
1229 30 : mpz_clear(sk);
1230 30 : mpz_clear(tk);
1231 30 : mpres_clear(tmp, n);
1232 30 : return ret;
1233 : }
1234 :
1235 : /* JKL: K = Q(sqrt(-3), sqrt(8*t^3+1), t in Q, t != 0, 1, -1/2;
1236 : mu = (2*t^3+1)/(3*t^2) => parameter for Hessian form.
1237 : Tors(E) = Z/6xZ/6.
1238 : A "specified" point is (0:-1:1), but does it have infinite order?
1239 : Also: twisted Hessian is a*X^3+Y^3+Z^3=d*X*Y*Z, d/a=3*mu.
1240 : See JKL-ECM in ANTS-XII.
1241 : */
1242 :
1243 : /* Original source is Brier + Clavier.
1244 : We can build curves in Montgomery form directly...
1245 : Useful if one knows that all p | n are 1 mod 4 (Cunningham, etc.).
1246 : */
1247 : int
1248 30 : build_curves_with_torsion_Z4xZ4(mpz_t f, mpmod_t n, ell_curve_t *tE,
1249 : ell_point_t *tP,
1250 : int smin, int smax, int nE)
1251 : {
1252 : mpz_t tau, lambda, nu2, tmp, b, x0;
1253 30 : int nu, nc = 0, ret = ECM_NO_FACTOR_FOUND;
1254 :
1255 30 : mpz_init(tau);
1256 30 : mpz_init(lambda);
1257 30 : mpz_init(nu2);
1258 30 : mpz_init(tmp);
1259 30 : mpz_init(b);
1260 30 : mpz_init(x0);
1261 30 : for(nu = smin; nu < smax; nu++){
1262 30 : mpz_set_si(nu2, nu*nu);
1263 : /* tau:=(nu^2+3)/2/nu; */
1264 30 : mpz_add_si(lambda, nu2, 3);
1265 30 : mpz_set_si(tmp, 2*nu);
1266 30 : if(mod_from_rat2(tau, lambda, tmp, n->orig_modulus) == 0){
1267 10 : printf("Factor found during init of Z4xZ4 (tau)\n");
1268 10 : mpz_set(f, tau);
1269 10 : ret = ECM_FACTOR_FOUND_STEP1;
1270 10 : break;
1271 : }
1272 : /* lambda:=8*nu^3; */
1273 20 : mpz_mul_si(lambda, nu2, 8*nu);
1274 20 : mpz_mod(lambda, lambda, n->orig_modulus);
1275 : /* A:=-27*lambda^4*(tau^8+14*tau^4+1); */
1276 : /* B:=54*lambda^6*(tau^12-33*tau^8-33*tau^4+1); */
1277 : /* x0:=3*(3*nu^12+34*nu^10+117*nu^8+316*nu^6+1053*nu^4+2754*nu^2+2187); */
1278 : /* y0:=27*(nu^2-3)*(nu^2+1)*(nu^2+9)*(nu^6+5*nu^4+15*nu^2+27)^2; */
1279 : /* P = (x0, y0) is a point on Y^2 = X^3+A*X+B */
1280 :
1281 : /* Montgomery form: there are several possible mb */
1282 : /* mb:=1/(9*lambda^2*(tau^4-1));
1283 : lambda is invertible iff nu is;
1284 : tau^4-1 = (tau-1)(tau+1)(tau^2+1)
1285 : */
1286 20 : mpz_powm_ui(x0, tau, 4, n->orig_modulus);
1287 20 : mpz_sub_si(x0, x0, 1);
1288 20 : mpz_mod(x0, x0, n->orig_modulus);
1289 20 : mpz_mul(tmp, x0, lambda);
1290 20 : mpz_mul(tmp, tmp, lambda);
1291 20 : mpz_mul_si(tmp, tmp, 9);
1292 20 : if(mpz_invert(b, tmp, n->orig_modulus) == 0){
1293 10 : printf("Factor found during init of Z4xZ4 (mb)\n");
1294 10 : mpz_gcd(f, tmp, n->orig_modulus);
1295 10 : ret = ECM_FACTOR_FOUND_STEP1;
1296 10 : break;
1297 : }
1298 : /* ma:=-2*(tau^4+1)/(tau^4-1); at this point: invertible! */
1299 10 : mpz_add_si(tmp, x0, 2);
1300 10 : mpz_mul_si(tmp, tmp, -2);
1301 10 : mpz_mod(tmp, tmp, n->orig_modulus);
1302 : /* to Montgomery form */
1303 10 : ell_curve_init(tE[nc], ECM_EC_TYPE_MONTGOMERY, ECM_LAW_HOMOGENEOUS, n);
1304 10 : ell_point_init(tP[nc], tE[nc], n);
1305 10 : mod_from_rat2(tE[nc]->a4, tmp, x0, n->orig_modulus);
1306 : /* now compute real x0 */
1307 : /* x0:=3*(3*nu^12+34*nu^10+117*nu^8+316*nu^6+1053*nu^4+2754*nu^2+2187); */
1308 10 : mpz_set_si(x0, 3);
1309 10 : mpz_mul(x0, x0, nu2); mpz_add_si(x0, x0, 34);
1310 10 : mpz_mul(x0, x0, nu2); mpz_add_si(x0, x0, 117);
1311 10 : mpz_mul(x0, x0, nu2); mpz_add_si(x0, x0, 316);
1312 10 : mpz_mul(x0, x0, nu2); mpz_add_si(x0, x0, 1053);
1313 10 : mpz_mul(x0, x0, nu2); mpz_add_si(x0, x0, 2754);
1314 10 : mpz_mul(x0, x0, nu2); mpz_add_si(x0, x0, 2187);
1315 10 : mpz_mul_si(x0, x0, 3);
1316 10 : mpz_mod(x0, x0, n->orig_modulus);
1317 : #if DEBUG_TORSION >= 2
1318 : gmp_printf("N:=%Zd;\n", n);
1319 : printf("nu:=%d;\n", nu);
1320 : gmp_printf("tau:=%Zd;\n", tau);
1321 : gmp_printf("lambda:=%Zd;\n", lambda);
1322 : gmp_printf("a:=%Zd;\n", tE[nc]->a4);
1323 : gmp_printf("x0:=%Zd;\n", x0);
1324 : #endif
1325 : /* x:=b*x0-a/3; not needed: y:=b*y0 */
1326 10 : mpz_set_si(tmp, 3);
1327 10 : mod_from_rat2(tP[nc]->x, tE[nc]->a4, tmp, n->orig_modulus);
1328 10 : mpz_mul(b, b, x0);
1329 10 : mpz_mod(b, b, n->orig_modulus);
1330 10 : mpz_sub(tP[nc]->x, b, tP[nc]->x);
1331 10 : mpz_mod(tP[nc]->x, tP[nc]->x, n->orig_modulus);
1332 10 : nc++;
1333 10 : if(nc >= nE)
1334 10 : break;
1335 : }
1336 30 : mpz_clear(tau);
1337 30 : mpz_clear(lambda);
1338 30 : mpz_clear(nu2);
1339 30 : mpz_clear(tmp);
1340 30 : mpz_clear(b);
1341 30 : mpz_clear(x0);
1342 30 : if(ret != ECM_FACTOR_FOUND_STEP1 && nc < nE){
1343 0 : printf("Not enough curves generated\n");
1344 0 : return ECM_ERROR;
1345 : }
1346 30 : return ret;
1347 : }
1348 :
1349 : /* Assuming we can generate curves with given torsion using parameter s
1350 : in interval [smin..smax[.
1351 : */
1352 : int
1353 380 : build_curves_with_torsion(mpz_t f, mpmod_t n, ell_curve_t *tE, ell_point_t *tP,
1354 : char *torsion, int smin, int smax, int nE)
1355 : {
1356 380 : int ret = 0;
1357 :
1358 : /* over Q: see Atkin-Morain, Math. Comp., 1993 */
1359 380 : if(strcmp(torsion, "Z5") == 0)
1360 20 : return build_curves_with_torsion_Z5(f, n, tE, tP, smin, smax, nE);
1361 360 : else if(strcmp(torsion, "Z7") == 0)
1362 70 : return build_curves_with_torsion_Z7(f, n, tE, tP, smin, smax, nE);
1363 290 : else if(strcmp(torsion, "Z9") == 0)
1364 40 : return build_curves_with_torsion_Z9(f, n, tE, tP, smin, smax, nE);
1365 250 : else if(strcmp(torsion, "Z10") == 0)
1366 50 : return build_curves_with_torsion_Z10(f, n, tE, tP, smin, smax, nE);
1367 200 : else if(strcmp(torsion, "Z2xZ8") == 0)
1368 80 : return build_curves_with_torsion_Z2xZ8(f, n, tE, tP, smin, smax, nE);
1369 : /* no longer over Q */
1370 : /** interesting when p = 1 mod 3 **/
1371 120 : else if(strcmp(torsion, "Z3xZ3") == 0) /* over Q(sqrt(-3)) */
1372 50 : return build_curves_with_torsion_Z3xZ3(f, n, tE, tP, smin, smax, nE);
1373 70 : else if(strcmp(torsion, "Z3xZ6") == 0) /* over Q(sqrt(-3)) */
1374 30 : return build_curves_with_torsion_Z3xZ6(f, n, tE, tP, smin, smax, nE);
1375 : /** interesting when p = 1 mod 4 **/
1376 40 : else if(strcmp(torsion, "Z4xZ4") == 0) /* over Q(sqrt(-1)) */
1377 30 : return build_curves_with_torsion_Z4xZ4(f, n, tE, tP, smin, smax, nE);
1378 : else{
1379 10 : printf("Unknown torsion group: %s\n", torsion);
1380 10 : ret = ECM_ERROR;
1381 : }
1382 10 : return ret;
1383 : }
1384 :
1385 : /* E is a curve with given torsion and (x, y) a point on E mod n.
1386 : OUTPUT: ECM_NO_FACTOR_FOUND if everything went ok
1387 : ECM_FACTOR_FOUND_STEP1 in case a factor was found when building E.
1388 : REM: E is defined over Z, not in mpres_t.
1389 : */
1390 : int
1391 380 : build_curves_with_torsion2(mpz_t f, mpz_t n, ell_curve_t E,
1392 : mpz_t x, mpz_t y, char *torsion,
1393 : mpz_t sigma)
1394 : {
1395 : ell_curve_t tE[1];
1396 : ell_point_t tP[1];
1397 : mpmod_t modulus;
1398 : int ret, smin, smax;
1399 :
1400 380 : smin = (int)mpz_get_si(sigma);
1401 380 : smax = smin+10;
1402 380 : mpmod_init(modulus, n, ECM_MOD_DEFAULT);
1403 380 : ret = build_curves_with_torsion(f, modulus, tE, tP, torsion, smin, smax,1);
1404 380 : if(ret == ECM_NO_FACTOR_FOUND){
1405 130 : E->type = tE[0]->type;
1406 130 : E->law = tE[0]->law;
1407 130 : mpz_set(E->a2, tE[0]->a2);
1408 130 : mpz_set(E->a4, tE[0]->a4);
1409 130 : mpz_set(E->a6, tE[0]->a6);
1410 130 : mpz_set(x, tP[0]->x);
1411 130 : mpz_set(y, tP[0]->y);
1412 130 : ell_point_clear(tP[0], tE[0], modulus);
1413 130 : ell_curve_clear(tE[0], modulus);
1414 : }
1415 380 : mpmod_clear(modulus);
1416 380 : return ret;
1417 : }
|
