1 /* parametrizations.c - functions to compute coefficients of the curve from
2 parameter and to choose random parameter.
3
4 Copyright 2012 Cyril Bouvier.
5
6 This file is part of the ECM Library.
7
8 The ECM Library is free software; you can redistribute it and/or modify
9 it under the terms of the GNU Lesser General Public License as published by
10 the Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12
13 The ECM Library is distributed in the hope that it will be useful, but
14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16 License for more details.
17
18 You should have received a copy of the GNU Lesser General Public License
19 along with the ECM Library; see the file COPYING.LIB. If not, see
20 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
22
23 #include "ecm-gmp.h"
24 #include "ecm-impl.h"
25
26 #if 0
27 /* this function is useful in debug mode to print residues */
28 static void
29 mpres_print (mpres_t x, char* name, mpmod_t n)
30 {
31 mp_size_t m, xn;
32 mpres_t t;
33 mpres_init(t, n);
34 mpz_set_ui(t, 1);
35 mpres_mul (t, x, t, n);
36
37 xn = SIZ(t);
38 m = ABSIZ(t);
39 MPN_NORMALIZE(PTR(t), m);
40 SIZ(t) = xn >= 0 ? m : -m;
41 gmp_printf ("%s=%Zd\n", name, t);
42 SIZ(t) = xn;
43 mpres_clear (t, n);
44 }
45 #endif
46
47 static void
dbl_param(mpres_t x,mpres_t y,mpres_t z,mpres_t t,mpres_t u,mpres_t v,mpmod_t n)48 dbl_param (mpres_t x, mpres_t y, mpres_t z, mpres_t t, mpres_t u, mpres_t v,
49 mpmod_t n)
50 {
51 mpres_mul (z, y, z, n); /* Y1*Z1 */
52 mpres_mul_ui (z, z, 2, n); /* Z3 = 2*Y1*Z1 */
53
54 mpres_sqr (u, x, n); /* A = X1*X1 */
55 mpres_sqr (t, y, n); /* B = Y1*Y1 */
56 mpres_sqr (y, t, n); /* C = B^2 */
57 mpres_add (v, x, t, n); /* X1+B */
58 mpres_sqr (v, v, n); /* (X1+B)^2 */
59 mpres_sub (v, v, u, n); /* (X1+B)^2-A */
60 mpres_sub (v, v, y, n); /* (X1+B)^2-A-C */
61 mpres_mul_ui (v, v, 2, n); /* D = 2*((X1+B)^2-A-C) */
62 mpres_mul_ui (u, u, 3, n); /* E = 3*A */
63 mpres_sqr (t, u, n); /* F = E^2 */
64
65 mpres_mul_ui (x, v, 2, n); /* 2*D */
66 mpres_sub (x, t, x, n); /* X3 = F-2*D */
67
68 mpres_sub (v, v, x, n); /* D-X3 */
69 mpres_mul_ui (y, y, 8, n); /* 8*C */
70 mpres_mul (t, u, v, n); /* E*(D-X3) */
71 mpres_sub (y, t, y, n); /* Y3 = E*(D-X3)-8*C */
72 }
73
74 /*Add sgn*P=(-3:sgn*3:1) to Q=(x:y:z) */
75 static void
add_param(mpres_t x,mpres_t y,mpres_t z,int sgn,mpres_t t,mpres_t u,mpres_t v,mpres_t w,mpmod_t n)76 add_param (mpres_t x, mpres_t y, mpres_t z, int sgn, mpres_t t, mpres_t u,
77 mpres_t v, mpres_t w, mpmod_t n)
78 {
79 mpres_sqr (t, z, n); /* Z1Z1 = Z1^2 */
80 mpres_mul_ui (u, t, 3, n);
81 mpres_neg (u, u, n); /* U2 = X2*Z1Z1 with X2=-3 */
82 mpres_mul (v, z, t, n); /* Z1*Z1Z1 */
83 mpres_mul_ui (v, v, 3, n); /* S2 = Y2*Z1*Z1Z1 with Y2=sgn*3 */
84 if (sgn == -1)
85 mpres_neg (v, v, n); /* S2 = Y2*Z1*Z1Z1 with Y2=sgn*3 */
86 mpres_sub (u, u, x, n); /* H = U2-X1 */
87 mpres_sqr (w, u, n); /* HH = H^2 */
88
89 mpres_add (z, z, u, n); /* Z1+H */
90 mpres_sqr (z, z, n); /* (Z1+H)^2 */
91 mpres_sub (z, z, t, n); /* (Z1+H)^2-Z1Z1 */
92 mpres_sub (z, z, w, n); /* Z3 = (Z1+H)^2-Z1Z1-HH */
93
94 mpres_mul_ui (t, w, 4, n); /* I = 4*HH */
95 mpres_mul (u, u, t, n); /* J = H*I */
96 mpres_sub (v, v, y, n); /* S2-Y1 */
97 mpres_mul_ui (v, v, 2, n); /* r = 2*(S2-Y1) */
98 mpres_mul (t, x, t, n); /* V = X1*I */
99 mpres_sqr (x, v, n); /* r^2 */
100 mpres_mul_ui (w, t, 2, n); /* 2*V */
101 mpres_sub (x, x, u, n); /* r^2-J */
102 mpres_sub (x, x, w, n); /* X3 = r^2-J-2*V */
103
104 mpres_sub (w, t, x, n); /* V-X3 */
105 mpres_mul (y, y, u, n); /* Y1*J */
106 mpres_mul_ui (y, y, 2, n); /* 2*Y1*J */
107 mpres_mul (w, v, w, n); /* r*(V-X3) */
108 mpres_sub (y, w, y, n); /* Y3=r*(V-X3)-2*Y1*J */
109 }
110
111 /* computes s*(x:y:z) mod n, where t, u, v, w are temporary variables */
112 static void
addchain_param(mpres_t x,mpres_t y,mpres_t z,mpz_t s,mpres_t t,mpres_t u,mpres_t v,mpres_t w,mpmod_t n)113 addchain_param (mpres_t x, mpres_t y, mpres_t z, mpz_t s, mpres_t t,
114 mpres_t u, mpres_t v, mpres_t w, mpmod_t n)
115 {
116 if (mpz_cmp_ui (s, 1) == 0)
117 {
118 mpres_set_si (x, -3, n);
119 mpres_set_ui (y, 3, n);
120 mpres_set_ui (z, 1, n);
121 }
122 else if (mpz_cmp_ui (s, 3) == 0)
123 {
124 mpz_sub_ui (s, s, 1);
125 addchain_param (x, y, z, s, t, u, v, w, n);
126 add_param (x, y, z, +1, t, u, v, w, n);
127 }
128 else if (mpz_divisible_2exp_p (s, 1))
129 {
130 mpz_tdiv_q_2exp (s, s, 1);
131 addchain_param (x, y, z, s, t, u, v, w, n);
132 dbl_param (x, y, z, t, u, v, n);
133 }
134 else if (mpz_congruent_ui_p (s, 1, 4))
135 {
136 mpz_sub_ui (s, s, 1);
137 addchain_param (x, y, z, s, t, u, v, w, n);
138 add_param (x, y, z, +1, t, u, v, w, n);
139 }
140 else /* (s % 4 == 3) and s != 3 */
141 {
142 mpz_add_ui (s, s, 1);
143 addchain_param (x, y, z, s, t, u, v, w, n);
144 add_param (x, y, z, -1, t, u, v, w, n);
145 }
146 }
147
148 /*
149 get_curve_from_param* functions compute curve coefficient A and a starting
150 point (x::1) from a given sigma value
151
152 If a factor of n was found during the process, returns
153 ECM_FACTOR_FOUND_STEP1 (and factor in f).
154 If a invalid value of sigma is given, returns ECM_ERROR,
155 Returns ECM_NO_FACTOR_FOUND otherwise.
156 */
157
158
159
160 /* Parametrization ECM_PARAM_SUYAMA */
161 /* (sigma mod N) should be different from 0, 1, 3, 5, 5/3, -1, -3, -5, -5/3 */
162 int
get_curve_from_param0(mpz_t f,mpres_t A,mpres_t x,mpz_t sigma,mpmod_t n)163 get_curve_from_param0 (mpz_t f, mpres_t A, mpres_t x, mpz_t sigma, mpmod_t n)
164 {
165 mpres_t t, u, v, b, z;
166 mpz_t tmp;
167
168 mpres_init (t, n);
169 mpres_init (u, n);
170 mpres_init (v, n);
171 mpres_init (b, n);
172 mpres_init (z, n);
173 mpz_init (tmp);
174
175 mpz_mod (tmp, sigma, n->orig_modulus);
176 /* TODO add -5 -3 -1 and +/- 5/3 */
177 if (mpz_cmp_ui (tmp, 5) == 0 || mpz_cmp_ui (tmp, 3) == 0 ||
178 mpz_cmp_ui (tmp, 1) == 0 || mpz_sgn (tmp) == 0)
179 return ECM_ERROR;
180
181 mpres_set_z (u, sigma, n);
182 mpres_mul_ui (v, u, 4, n); /* v = (4*sigma) mod n */
183 mpres_sqr (t, u, n);
184 mpres_sub_ui (u, t, 5, n); /* u = (sigma^2-5) mod n */
185 mpres_sqr (t, u, n);
186 mpres_mul (x, t, u, n); /* x = (u^3) mod n */
187 mpres_sqr (t, v, n);
188 mpres_mul (z, t, v, n); /* z = (v^3) mod n */
189 mpres_mul (t, x, v, n);
190 mpres_mul_ui (b, t, 4, n); /* b = (4*x*v) mod n */
191 mpres_mul_ui (t, u, 3, n);
192 mpres_sub (u, v, u, n); /* u' = v-u */
193 mpres_add (v, t, v, n); /* v' = (3*u+v) mod n */
194 mpres_sqr (t, u, n);
195 mpres_mul (u, t, u, n); /* u'' = ((v-u)^3) mod n */
196 mpres_mul (A, u, v, n); /* a = (u'' * v') mod n =
197 ((v-u)^3 * (3*u+v)) mod n */
198
199 /* Normalize b and z to 1 */
200 mpres_mul (v, b, z, n);
201 if (!mpres_invert (u, v, n)) /* u = (b*z)^(-1) (mod n) */
202 {
203 mpres_gcd (f, v, n);
204 mpres_clear (t, n);
205 mpres_clear (u, n);
206 mpres_clear (v, n);
207 mpres_clear (b, n);
208 mpres_clear (z, n);
209 mpz_clear (tmp);
210 if (mpz_cmp (f, n->orig_modulus) == 0)
211 return ECM_ERROR;
212 else
213 return ECM_FACTOR_FOUND_STEP1;
214 }
215
216 mpres_mul (v, u, b, n); /* v = z^(-1) (mod n) */
217 mpres_mul (x, x, v, n); /* x = x * z^(-1) */
218
219 mpres_mul (v, u, z, n); /* v = b^(-1) (mod n) */
220 mpres_mul (t, A, v, n);
221 mpres_sub_ui (A, t, 2, n);
222
223 mpres_clear (t, n);
224 mpres_clear (u, n);
225 mpres_clear (v, n);
226 mpres_clear (b, n);
227 mpres_clear (z, n);
228 mpz_clear (tmp);
229
230 return ECM_NO_FACTOR_FOUND;
231 }
232
233 /* Parametrization ECM_PARAM_BATCH_SQUARE */
234 /* Only work for 64-bit machines */
235 /* d = (sigma^2/2^64 mod N) should be different from 0, 1, -1/8 */
236 int
get_curve_from_param1(mpres_t A,mpres_t x0,mpz_t sigma,mpmod_t n)237 get_curve_from_param1 (mpres_t A, mpres_t x0, mpz_t sigma, mpmod_t n)
238 {
239 int i;
240 mpz_t tmp;
241 mpz_init (tmp);
242
243 ASSERT (GMP_NUMB_BITS == 64);
244
245 mpz_mul (tmp, sigma, sigma); /* tmp = sigma^2*/
246
247 /* A=4*d-2 with d = sigma^2/2^GMP_NUMB_BITS*/
248 /* Compute d = sigma^2/2^GMP_NUMB_BITS */
249 for (i = 0; i < GMP_NUMB_BITS; i++)
250 {
251 if (mpz_tstbit (tmp, 0) == 1)
252 mpz_add (tmp, tmp, n->orig_modulus);
253 mpz_div_2exp (tmp, tmp, 1);
254 }
255
256 mpz_mod (tmp, tmp, n->orig_modulus);
257 /* TODO add d!=-1/8*/
258 if (mpz_sgn (tmp) == 0 || mpz_cmp_ui (tmp, 1) == 0)
259 return ECM_ERROR;
260
261 mpz_mul_2exp (tmp, tmp, 2); /* 4d */
262 mpz_sub_ui (tmp, tmp, 2); /* 4d-2 */
263
264 mpres_set_z (A, tmp, n);
265 mpres_set_ui (x0, 2, n);
266
267 mpz_clear(tmp);
268 return ECM_NO_FACTOR_FOUND;
269 }
270
271 /* Parametrization ECM_PARAM_BATCH_2 */
272 /* 2 <= sigma */
273 /* Compute (x:y:z) = sigma*(-3:3:1) on the elliptic curve y^2 = x^3 + 36
274 using Jacobian coordinates; formulae were found at
275 https://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-0.html
276 Then we let x3 = (3*x+y+6)/(2*(y-3)), A = -(3*x3^4+6*x3^2-1)/(4*x3^3) and
277 x0 = 2. The Sage code below gives the factored group order:
278
279 def FindGroupOrderParam2(p,sigma):
280 K = GF(p)
281 E = EllipticCurve(K,[0,36])
282 P = sigma*E(-3,3)
283 x,y = P.xy()
284 x3 = (3*x+y+6)/(2*(y-3))
285 A = -(3*x3^4+6*x3^2-1)/(4*x3^3)
286 d = K((A+2)/4)
287 a = K(4*d-2)
288 b = K(16*d+2)
289 E = EllipticCurve(K,[0,a/b,0,1/b^2,0])
290 return factor(E.cardinality())
291 */
292 int
get_curve_from_param2(mpz_t f,mpres_t A,mpres_t x0,mpz_t sigma,mpmod_t n)293 get_curve_from_param2 (mpz_t f, mpres_t A, mpres_t x0, mpz_t sigma, mpmod_t n)
294 {
295 mpres_t t, u, v, w, x, y, z;
296 mpz_t k;
297 int ret = ECM_NO_FACTOR_FOUND;
298
299 mpres_init (t, n);
300 mpres_init (u, n);
301 mpres_init (v, n);
302 mpres_init (w, n);
303 mpres_init (x, n);
304 mpres_init (y, n);
305 mpres_init (z, n);
306 mpz_init (k);
307
308 mpz_set (k, sigma);
309
310 if (mpz_cmp_ui (k, 2) < 0)
311 {
312 ret = ECM_ERROR;
313 goto clear_and_exit;
314 }
315
316 addchain_param (x, y, z, k, t, u, v, w, n);
317
318 /* Now (x:y:z) = k*P */
319
320 if (!mpres_invert (u, z, n))
321 {
322 mpres_gcd (f, z, n);
323 ret = ECM_FACTOR_FOUND_STEP1;
324 goto clear_and_exit;
325 }
326
327 mpres_sqr (v, u, n);
328 mpres_mul (u, v, u, n);
329 mpres_mul (x, x, v, n);
330 mpres_mul (y, y, u, n);
331
332 mpres_sub_ui (t, y, 3, n);
333 mpres_mul_ui (t, t, 2, n);
334
335 if (!mpres_invert (u, t, n))
336 {
337 mpres_gcd (f, t, n);
338 ret = ECM_FACTOR_FOUND_STEP1;
339 goto clear_and_exit;
340 }
341
342 mpres_mul_ui (w, x, 3, n);
343 mpres_add (w, w, y, n);
344 mpres_add_ui (w, w, 6, n);
345 mpres_mul (x, w, u, n); /* Now x contains x_3 */
346
347 /* A=-(3*x3^4+6*x3^2-1)/(4*x3^3) */
348 mpres_sqr (u, x, n);
349 mpres_mul (v, u, x, n);
350 mpres_sqr (w, u, n);
351
352 mpres_mul_ui (u, u, 6, n);
353 mpres_neg (u, u, n);
354 mpres_mul_ui (v, v, 4, n);
355 mpres_mul_ui (w, w, 3, n);
356 mpres_neg (w, w, n);
357
358 if (!mpres_invert (t, v, n))
359 {
360 mpres_gcd (f, v, n);
361 ret = ECM_FACTOR_FOUND_STEP1;
362 goto clear_and_exit;
363 }
364
365 mpres_add (w, w, u, n);
366 mpres_add_ui (w, w, 1, n);
367 mpres_mul (A, w, t, n);
368 mpz_mod (A, A, n->orig_modulus);
369
370 mpres_set_ui (x0, 2, n);
371
372 clear_and_exit:
373 mpres_clear (t, n);
374 mpres_clear (u, n);
375 mpres_clear (v, n);
376 mpres_clear (w, n);
377 mpres_clear (x, n);
378 mpres_clear (y, n);
379 mpres_clear (z, n);
380 mpz_clear (k);
381
382 return ret;
383 }
384
385 /* Parametrization ECM_PARAM_BATCH_32BITS_D */
386 /* d = (sigma/2^32 mod N) should be different from 0, 1, -1/8 */
387 int
get_curve_from_param3(mpres_t A,mpres_t x0,mpz_t sigma,mpmod_t n)388 get_curve_from_param3 (mpres_t A, mpres_t x0, mpz_t sigma, mpmod_t n)
389 {
390 int i;
391 mpz_t tmp;
392 mpz_t two32;
393 mpz_init (two32);
394 mpz_ui_pow_ui (two32, 2, 32);
395 mpz_init (tmp);
396
397 /* sigma < 2^32 (it was generated for 32-bit machines) */
398 /* To use it on a 64-bits machines one should multiplied it by 2^32 */
399 if (GMP_NUMB_BITS == 64)
400 mpz_mul (tmp, sigma, two32);
401 else
402 mpz_set (tmp, sigma);
403
404 /* A=4*d-2 with d = sigma/2^GMP_NUMB_BITS*/
405 /* Compute d = sigma/2^GMP_NUMB_BITS */
406 for (i = 0; i < GMP_NUMB_BITS; i++)
407 {
408 if (mpz_tstbit (tmp, 0) == 1)
409 mpz_add (tmp, tmp, n->orig_modulus);
410 mpz_div_2exp (tmp, tmp, 1);
411 }
412
413 mpz_mod (tmp, tmp, n->orig_modulus);
414 /* TODO add d!=-1/8*/
415 if (mpz_sgn (tmp) == 0 || mpz_cmp_ui (tmp, 1) == 0)
416 return ECM_ERROR;
417
418 mpz_mul_2exp (tmp, tmp, 2); /* 4d */
419 mpz_sub_ui (tmp, tmp, 2); /* 4d-2 */
420
421 mpres_set_z (A, tmp, n);
422 mpres_set_ui (x0, 2, n);
423
424 mpz_clear(tmp);
425 mpz_clear (two32);
426 return ECM_NO_FACTOR_FOUND;
427 }
428
429 int
get_curve_from_random_parameter(mpz_t f,mpres_t A,mpres_t x,mpz_t sigma,int param,mpmod_t modulus,gmp_randstate_t rng)430 get_curve_from_random_parameter (mpz_t f, mpres_t A, mpres_t x, mpz_t sigma,
431 int param, mpmod_t modulus, gmp_randstate_t rng)
432 {
433 int ret;
434
435 /* initialize the random number generator if not already done */
436 init_randstate (rng);
437 do
438 {
439 if (param == ECM_PARAM_SUYAMA)
440 {
441 mpz_urandomb (sigma, rng, 64);
442 ret = get_curve_from_param0 (f, A, x, sigma, modulus);
443 }
444 else if (param == ECM_PARAM_BATCH_SQUARE)
445 {
446 mpz_urandomb (sigma, rng, 32);
447 ret = get_curve_from_param1 (A, x, sigma, modulus);
448 }
449 else if (param == ECM_PARAM_BATCH_2)
450 {
451 mpz_urandomb (sigma, rng, 64);
452 ret = get_curve_from_param2 (f, A, x, sigma, modulus);
453 }
454 else if (param == ECM_PARAM_BATCH_32BITS_D)
455 {
456 mpz_urandomb (sigma, rng, 32);
457 ret = get_curve_from_param3 (A, x, sigma, modulus);
458 }
459 else
460 return ECM_ERROR;
461 } while (ret == ECM_ERROR);
462
463 return ret;
464 }
465
466 int
get_default_param(int sigma_is_A,double B1done,int repr)467 get_default_param (int sigma_is_A, double B1done, int repr)
468 {
469 /* if B1done is not the default value, use ECM_PARAM_SUYAMA, since
470 ECM_PARAM_BATCH* requires B1done is the default */
471 if (!ECM_IS_DEFAULT_B1_DONE(B1done))
472 return ECM_PARAM_SUYAMA;
473
474 if (sigma_is_A == 1 || sigma_is_A == -1)
475 {
476 /* For now we keep the default values in order not to compute the
477 expected number of curves. But it will do stage 1 as ECM_PARAM_SUYAMA */
478 return ECM_PARAM_DEFAULT;
479 }
480
481 /* ECM_PARAM_BATCH* requires ECM_MOD_MODMULN */
482 if (repr != ECM_MOD_MODMULN)
483 return ECM_PARAM_SUYAMA;
484
485 if (GMP_NUMB_BITS == 64)
486 return ECM_PARAM_BATCH_SQUARE;
487 else
488 return ECM_PARAM_BATCH_32BITS_D;
489 }
490