Line data Source code
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
48 1224 : 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 1224 : mpres_mul (z, y, z, n); /* Y1*Z1 */
52 1224 : mpres_mul_ui (z, z, 2, n); /* Z3 = 2*Y1*Z1 */
53 :
54 1224 : mpres_sqr (u, x, n); /* A = X1*X1 */
55 1224 : mpres_sqr (t, y, n); /* B = Y1*Y1 */
56 1224 : mpres_sqr (y, t, n); /* C = B^2 */
57 1224 : mpres_add (v, x, t, n); /* X1+B */
58 1224 : mpres_sqr (v, v, n); /* (X1+B)^2 */
59 1224 : mpres_sub (v, v, u, n); /* (X1+B)^2-A */
60 1224 : mpres_sub (v, v, y, n); /* (X1+B)^2-A-C */
61 1224 : mpres_mul_ui (v, v, 2, n); /* D = 2*((X1+B)^2-A-C) */
62 1224 : mpres_mul_ui (u, u, 3, n); /* E = 3*A */
63 1224 : mpres_sqr (t, u, n); /* F = E^2 */
64 :
65 1224 : mpres_mul_ui (x, v, 2, n); /* 2*D */
66 1224 : mpres_sub (x, t, x, n); /* X3 = F-2*D */
67 :
68 1224 : mpres_sub (v, v, x, n); /* D-X3 */
69 1224 : mpres_mul_ui (y, y, 8, n); /* 8*C */
70 1224 : mpres_mul (t, u, v, n); /* E*(D-X3) */
71 1224 : mpres_sub (y, t, y, n); /* Y3 = E*(D-X3)-8*C */
72 1224 : }
73 :
74 : /*Add sgn*P=(-3:sgn*3:1) to Q=(x:y:z) */
75 : static void
76 296 : 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 296 : mpres_sqr (t, z, n); /* Z1Z1 = Z1^2 */
80 296 : mpres_mul_ui (u, t, 3, n);
81 296 : mpres_neg (u, u, n); /* U2 = X2*Z1Z1 with X2=-3 */
82 296 : mpres_mul (v, z, t, n); /* Z1*Z1Z1 */
83 296 : mpres_mul_ui (v, v, 3, n); /* S2 = Y2*Z1*Z1Z1 with Y2=sgn*3 */
84 296 : if (sgn == -1)
85 144 : mpres_neg (v, v, n); /* S2 = Y2*Z1*Z1Z1 with Y2=sgn*3 */
86 296 : mpres_sub (u, u, x, n); /* H = U2-X1 */
87 296 : mpres_sqr (w, u, n); /* HH = H^2 */
88 :
89 296 : mpres_add (z, z, u, n); /* Z1+H */
90 296 : mpres_sqr (z, z, n); /* (Z1+H)^2 */
91 296 : mpres_sub (z, z, t, n); /* (Z1+H)^2-Z1Z1 */
92 296 : mpres_sub (z, z, w, n); /* Z3 = (Z1+H)^2-Z1Z1-HH */
93 :
94 296 : mpres_mul_ui (t, w, 4, n); /* I = 4*HH */
95 296 : mpres_mul (u, u, t, n); /* J = H*I */
96 296 : mpres_sub (v, v, y, n); /* S2-Y1 */
97 296 : mpres_mul_ui (v, v, 2, n); /* r = 2*(S2-Y1) */
98 296 : mpres_mul (t, x, t, n); /* V = X1*I */
99 296 : mpres_sqr (x, v, n); /* r^2 */
100 296 : mpres_mul_ui (w, t, 2, n); /* 2*V */
101 296 : mpres_sub (x, x, u, n); /* r^2-J */
102 296 : mpres_sub (x, x, w, n); /* X3 = r^2-J-2*V */
103 :
104 296 : mpres_sub (w, t, x, n); /* V-X3 */
105 296 : mpres_mul (y, y, u, n); /* Y1*J */
106 296 : mpres_mul_ui (y, y, 2, n); /* 2*Y1*J */
107 296 : mpres_mul (w, v, w, n); /* r*(V-X3) */
108 296 : mpres_sub (y, w, y, n); /* Y3=r*(V-X3)-2*Y1*J */
109 296 : }
110 :
111 : /* computes s*(x:y:z) mod n, where t, u, v, w are temporary variables */
112 : static void
113 1576 : 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 1576 : if (mpz_cmp_ui (s, 1) == 0)
117 : {
118 56 : mpres_set_si (x, -3, n);
119 56 : mpres_set_ui (y, 3, n);
120 56 : mpres_set_ui (z, 1, n);
121 : }
122 1520 : else if (mpz_cmp_ui (s, 3) == 0)
123 : {
124 8 : mpz_sub_ui (s, s, 1);
125 8 : addchain_param (x, y, z, s, t, u, v, w, n);
126 8 : add_param (x, y, z, +1, t, u, v, w, n);
127 : }
128 1512 : else if (mpz_divisible_2exp_p (s, 1))
129 : {
130 1224 : mpz_tdiv_q_2exp (s, s, 1);
131 1224 : addchain_param (x, y, z, s, t, u, v, w, n);
132 1224 : dbl_param (x, y, z, t, u, v, n);
133 : }
134 288 : else if (mpz_congruent_ui_p (s, 1, 4))
135 : {
136 144 : mpz_sub_ui (s, s, 1);
137 144 : addchain_param (x, y, z, s, t, u, v, w, n);
138 144 : add_param (x, y, z, +1, t, u, v, w, n);
139 : }
140 : else /* (s % 4 == 3) and s != 3 */
141 : {
142 144 : mpz_add_ui (s, s, 1);
143 144 : addchain_param (x, y, z, s, t, u, v, w, n);
144 144 : add_param (x, y, z, -1, t, u, v, w, n);
145 : }
146 1576 : }
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
163 1017 : 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 1017 : mpres_init (t, n);
169 1017 : mpres_init (u, n);
170 1017 : mpres_init (v, n);
171 1017 : mpres_init (b, n);
172 1017 : mpres_init (z, n);
173 1017 : mpz_init (tmp);
174 :
175 1017 : mpz_mod (tmp, sigma, n->orig_modulus);
176 : /* TODO add -5 -3 -1 and +/- 5/3 */
177 1017 : if (mpz_cmp_ui (tmp, 5) == 0 || mpz_cmp_ui (tmp, 3) == 0 ||
178 1017 : mpz_cmp_ui (tmp, 1) == 0 || mpz_sgn (tmp) == 0)
179 0 : return ECM_ERROR;
180 :
181 1017 : mpres_set_z (u, sigma, n);
182 1017 : mpres_mul_ui (v, u, 4, n); /* v = (4*sigma) mod n */
183 1017 : mpres_sqr (t, u, n);
184 1017 : mpres_sub_ui (u, t, 5, n); /* u = (sigma^2-5) mod n */
185 1017 : mpres_sqr (t, u, n);
186 1017 : mpres_mul (x, t, u, n); /* x = (u^3) mod n */
187 1017 : mpres_sqr (t, v, n);
188 1017 : mpres_mul (z, t, v, n); /* z = (v^3) mod n */
189 1017 : mpres_mul (t, x, v, n);
190 1017 : mpres_mul_ui (b, t, 4, n); /* b = (4*x*v) mod n */
191 1017 : mpres_mul_ui (t, u, 3, n);
192 1017 : mpres_sub (u, v, u, n); /* u' = v-u */
193 1017 : mpres_add (v, t, v, n); /* v' = (3*u+v) mod n */
194 1017 : mpres_sqr (t, u, n);
195 1017 : mpres_mul (u, t, u, n); /* u'' = ((v-u)^3) mod n */
196 1017 : 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 1017 : mpres_mul (v, b, z, n);
201 1017 : if (!mpres_invert (u, v, n)) /* u = (b*z)^(-1) (mod n) */
202 : {
203 20 : mpres_gcd (f, v, n);
204 20 : mpres_clear (t, n);
205 20 : mpres_clear (u, n);
206 20 : mpres_clear (v, n);
207 20 : mpres_clear (b, n);
208 20 : mpres_clear (z, n);
209 20 : mpz_clear (tmp);
210 20 : if (mpz_cmp (f, n->orig_modulus) == 0)
211 0 : return ECM_ERROR;
212 : else
213 20 : return ECM_FACTOR_FOUND_STEP1;
214 : }
215 :
216 997 : mpres_mul (v, u, b, n); /* v = z^(-1) (mod n) */
217 997 : mpres_mul (x, x, v, n); /* x = x * z^(-1) */
218 :
219 997 : mpres_mul (v, u, z, n); /* v = b^(-1) (mod n) */
220 997 : mpres_mul (t, A, v, n);
221 997 : mpres_sub_ui (A, t, 2, n);
222 :
223 997 : mpres_clear (t, n);
224 997 : mpres_clear (u, n);
225 997 : mpres_clear (v, n);
226 997 : mpres_clear (b, n);
227 997 : mpres_clear (z, n);
228 997 : mpz_clear (tmp);
229 :
230 997 : 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
237 25 : get_curve_from_param1 (mpres_t A, mpres_t x0, mpz_t sigma, mpmod_t n)
238 : {
239 : mpz_t tmp;
240 25 : mpz_init (tmp);
241 :
242 : ASSERT (GMP_NUMB_BITS == 64);
243 :
244 : /* A=4*d-2 with d = sigma^2/2^64 */
245 : /* Compute d = sigma^2/2^64 */
246 25 : mpz_ui_pow_ui(tmp, 2, 64);
247 25 : mpz_invert(tmp, tmp, n->orig_modulus);
248 :
249 : /* tmp = sigma^2/2^64 */
250 25 : mpz_mul (tmp, tmp, sigma);
251 25 : mpz_mul (tmp, tmp, sigma);
252 :
253 25 : mpz_mod (tmp, tmp, n->orig_modulus);
254 : /* TODO add d!=-1/8*/
255 25 : if (mpz_sgn (tmp) == 0 || mpz_cmp_ui (tmp, 1) == 0)
256 0 : return ECM_ERROR;
257 :
258 25 : mpz_mul_2exp (tmp, tmp, 2); /* 4d */
259 25 : mpz_sub_ui (tmp, tmp, 2); /* 4d-2 */
260 :
261 25 : mpres_set_z (A, tmp, n);
262 25 : mpres_set_ui (x0, 2, n);
263 :
264 25 : mpz_clear(tmp);
265 25 : return ECM_NO_FACTOR_FOUND;
266 : }
267 :
268 : /* Parametrization ECM_PARAM_BATCH_2 */
269 : /* 2 <= sigma */
270 : /* Compute (x:y:z) = sigma*(-3:3:1) on the elliptic curve y^2 = x^3 + 36
271 : using Jacobian coordinates; formulae were found at
272 : https://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-0.html
273 : Then we let x3 = (3*x+y+6)/(2*(y-3)), A = -(3*x3^4+6*x3^2-1)/(4*x3^3) and
274 : x0 = 2. The Sage code below gives the factored group order:
275 :
276 : def FindGroupOrderParam2(p,sigma):
277 : K = GF(p)
278 : E = EllipticCurve(K,[0,36])
279 : P = sigma*E(-3,3)
280 : x,y = P.xy()
281 : x3 = (3*x+y+6)/(2*(y-3))
282 : A = -(3*x3^4+6*x3^2-1)/(4*x3^3)
283 : d = K((A+2)/4)
284 : a = K(4*d-2)
285 : b = K(16*d+2)
286 : E = EllipticCurve(K,[0,a/b,0,1/b^2,0])
287 : return factor(E.cardinality())
288 : */
289 : int
290 56 : get_curve_from_param2 (mpz_t f, mpres_t A, mpres_t x0, mpz_t sigma, mpmod_t n)
291 : {
292 : mpres_t t, u, v, w, x, y, z;
293 : mpz_t k;
294 56 : int ret = ECM_NO_FACTOR_FOUND;
295 :
296 56 : mpres_init (t, n);
297 56 : mpres_init (u, n);
298 56 : mpres_init (v, n);
299 56 : mpres_init (w, n);
300 56 : mpres_init (x, n);
301 56 : mpres_init (y, n);
302 56 : mpres_init (z, n);
303 56 : mpz_init (k);
304 :
305 56 : mpz_set (k, sigma);
306 :
307 56 : if (mpz_cmp_ui (k, 2) < 0)
308 : {
309 0 : ret = ECM_ERROR;
310 0 : goto clear_and_exit;
311 : }
312 :
313 56 : addchain_param (x, y, z, k, t, u, v, w, n);
314 :
315 : /* Now (x:y:z) = k*P */
316 :
317 56 : if (!mpres_invert (u, z, n))
318 : {
319 8 : mpres_gcd (f, z, n);
320 8 : ret = ECM_FACTOR_FOUND_STEP1;
321 8 : goto clear_and_exit;
322 : }
323 :
324 48 : mpres_sqr (v, u, n);
325 48 : mpres_mul (u, v, u, n);
326 48 : mpres_mul (x, x, v, n);
327 48 : mpres_mul (y, y, u, n);
328 :
329 48 : mpres_sub_ui (t, y, 3, n);
330 48 : mpres_mul_ui (t, t, 2, n);
331 :
332 48 : if (!mpres_invert (u, t, n))
333 : {
334 8 : mpres_gcd (f, t, n);
335 8 : ret = ECM_FACTOR_FOUND_STEP1;
336 8 : goto clear_and_exit;
337 : }
338 :
339 40 : mpres_mul_ui (w, x, 3, n);
340 40 : mpres_add (w, w, y, n);
341 40 : mpres_add_ui (w, w, 6, n);
342 40 : mpres_mul (x, w, u, n); /* Now x contains x_3 */
343 :
344 : /* A=-(3*x3^4+6*x3^2-1)/(4*x3^3) */
345 40 : mpres_sqr (u, x, n);
346 40 : mpres_mul (v, u, x, n);
347 40 : mpres_sqr (w, u, n);
348 :
349 40 : mpres_mul_ui (u, u, 6, n);
350 40 : mpres_neg (u, u, n);
351 40 : mpres_mul_ui (v, v, 4, n);
352 40 : mpres_mul_ui (w, w, 3, n);
353 40 : mpres_neg (w, w, n);
354 :
355 40 : if (!mpres_invert (t, v, n))
356 : {
357 0 : mpres_gcd (f, v, n);
358 0 : ret = ECM_FACTOR_FOUND_STEP1;
359 0 : goto clear_and_exit;
360 : }
361 :
362 40 : mpres_add (w, w, u, n);
363 40 : mpres_add_ui (w, w, 1, n);
364 40 : mpres_mul (A, w, t, n);
365 40 : mpz_mod (A, A, n->orig_modulus);
366 :
367 40 : mpres_set_ui (x0, 2, n);
368 :
369 56 : clear_and_exit:
370 56 : mpres_clear (t, n);
371 56 : mpres_clear (u, n);
372 56 : mpres_clear (v, n);
373 56 : mpres_clear (w, n);
374 56 : mpres_clear (x, n);
375 56 : mpres_clear (y, n);
376 56 : mpres_clear (z, n);
377 56 : mpz_clear (k);
378 :
379 56 : return ret;
380 : }
381 :
382 : /* Parametrization ECM_PARAM_BATCH_32BITS_D */
383 : /* d = (sigma/2^32 mod N) should be different from 0, 1, -1/8 */
384 : int
385 32 : get_curve_from_param3 (mpres_t A, mpres_t x0, mpz_t sigma, mpmod_t n)
386 : {
387 : mpz_t tmp;
388 32 : mpz_init (tmp);
389 :
390 : /* A=4*d-2 with d = sigma/2^32*/
391 : /* Compute d = sigma/2^32 */
392 32 : mpz_ui_pow_ui (tmp, 2, 32);
393 32 : mpz_invert (tmp, tmp, n->orig_modulus);
394 32 : mpz_mul (tmp, sigma, tmp);
395 32 : mpz_mod (tmp, tmp, n->orig_modulus);
396 :
397 : /* TODO add d!=-1/8*/
398 32 : if (mpz_sgn (tmp) == 0 || mpz_cmp_ui (tmp, 1) == 0)
399 0 : return ECM_ERROR;
400 :
401 32 : mpz_mul_2exp (tmp, tmp, 2); /* 4d */
402 32 : mpz_sub_ui (tmp, tmp, 2); /* 4d-2 */
403 :
404 32 : mpres_set_z (A, tmp, n);
405 32 : mpres_set_ui (x0, 2, n);
406 :
407 32 : mpz_clear(tmp);
408 32 : return ECM_NO_FACTOR_FOUND;
409 : }
410 :
411 : int
412 29 : get_curve_from_random_parameter (mpz_t f, mpres_t A, mpres_t x, mpz_t sigma,
413 : int param, mpmod_t modulus, gmp_randstate_t rng)
414 : {
415 : int ret;
416 :
417 : /* initialize the random number generator if not already done */
418 29 : init_randstate (rng);
419 : do
420 : {
421 29 : if (param == ECM_PARAM_SUYAMA)
422 : {
423 10 : mpz_urandomb (sigma, rng, 64);
424 10 : ret = get_curve_from_param0 (f, A, x, sigma, modulus);
425 : }
426 19 : else if (param == ECM_PARAM_BATCH_SQUARE)
427 : {
428 9 : mpz_urandomb (sigma, rng, 32);
429 9 : ret = get_curve_from_param1 (A, x, sigma, modulus);
430 : }
431 10 : else if (param == ECM_PARAM_BATCH_2)
432 : {
433 0 : mpz_urandomb (sigma, rng, 64);
434 0 : ret = get_curve_from_param2 (f, A, x, sigma, modulus);
435 : }
436 10 : else if (param == ECM_PARAM_BATCH_32BITS_D)
437 : {
438 0 : mpz_urandomb (sigma, rng, 32);
439 0 : ret = get_curve_from_param3 (A, x, sigma, modulus);
440 : }
441 : else
442 10 : return ECM_ERROR;
443 19 : } while (ret == ECM_ERROR);
444 :
445 19 : return ret;
446 : }
447 :
448 : int
449 41 : get_default_param (int sigma_is_A, double B1done, int repr)
450 : {
451 : /* if B1done is not the default value, use ECM_PARAM_SUYAMA, since
452 : ECM_PARAM_BATCH* requires B1done is the default */
453 41 : if (!ECM_IS_DEFAULT_B1_DONE(B1done))
454 0 : return ECM_PARAM_SUYAMA;
455 :
456 41 : if (sigma_is_A == 1 || sigma_is_A == -1)
457 : {
458 : /* For now we keep the default values in order not to compute the
459 : expected number of curves. But it will do stage 1 as ECM_PARAM_SUYAMA */
460 20 : return ECM_PARAM_DEFAULT;
461 : }
462 :
463 : /* ECM_PARAM_BATCH* requires ECM_MOD_MODMULN */
464 21 : if (repr != ECM_MOD_MODMULN)
465 18 : return ECM_PARAM_SUYAMA;
466 :
467 : if (GMP_NUMB_BITS == 64)
468 3 : return ECM_PARAM_BATCH_SQUARE;
469 : else
470 : return ECM_PARAM_BATCH_32BITS_D;
471 : }
|
