Line data Source code
1 : /* Elliptic Curve Method implementation: stage 2 routines.
2 :
3 : Copyright 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2012 Paul Zimmermann,
4 : Alexander Kruppa, Pierrick Gaudry, Dave Newman.
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 <stdlib.h>
24 : #include <math.h>
25 : #include "ecm-impl.h"
26 :
27 : /* R_i <- q_i * S, 0 <= i < n, where q_i are large integers, S is a point on
28 : an elliptic curve. Uses max(bits in q_i) modular inversions (one less if
29 : max(q_i) is a power of 2). Needs up to n+2 cells in T.
30 : Returns whether factor was found or not found, factor goes into p.
31 : No error can occur.
32 : We should have n > 0.
33 : */
34 :
35 : static int
36 2086 : multiplyW2n (mpz_t p, point *R, curve *S, mpz_t *q, const unsigned int n,
37 : mpmod_t modulus, mpres_t u, mpres_t v, mpres_t *T,
38 : unsigned long *tot_muls, unsigned long *tot_gcds)
39 : {
40 : ecm_uint i, maxbit, k; /* k is the number of values to batch invert */
41 2086 : ecm_uint l, t, muls = 0, gcds = 0;
42 : #ifdef WANT_EXPCOST
43 : unsigned int hamweight = 0;
44 : #endif
45 2086 : int youpi = ECM_NO_FACTOR_FOUND;
46 : mpz_t flag; /* Used as bit field, keeps track of which R[i] contain partial results */
47 : point s; /* 2^t * S */
48 : mpz_t signs; /* Used as bit field, i-th bit is set iff q[i]<0 */
49 : #ifdef WANT_ASSERT
50 : mpz_t __dummy; /* used for local computations */
51 : #endif
52 :
53 : ASSERT(n > 0);
54 :
55 : /* We assume the code below works also when S is the neutral element. */
56 :
57 2086 : mpz_init2 (flag, n);
58 2086 : mpz_init2 (signs, n);
59 2086 : mpres_init (s.x, modulus);
60 2086 : mpres_init (s.y, modulus);
61 2086 : mpres_set (s.x, S->x, modulus);
62 2086 : mpres_set (s.y, S->y, modulus);
63 :
64 : /* Set maxbit to index of highest set bit among all the q[i] */
65 : /* Index of highest bit of q is sizeinbase(q, 2) - 1 */
66 2086 : maxbit = 0;
67 57324 : for (i = 0; i < n; i++)
68 : {
69 : /* We'll first compute positive multiples and change signs later */
70 55238 : if (mpz_sgn (q[i]) < 0)
71 : {
72 5292 : mpz_setbit (signs, i);;
73 5292 : mpz_neg (q[i], q[i]);
74 : }
75 :
76 : /* Multiplier == 0? Then set result to neutral element */
77 55238 : if (mpz_sgn (q[i]) == 0)
78 : {
79 20 : mpres_set_ui (R[i].x, 0, modulus);
80 20 : mpres_set_ui (R[i].y, 0, modulus);
81 : }
82 : #ifdef WANT_EXPCOST
83 : else
84 : hamweight += mpz_popcount (q[i]) - 1;
85 : #endif
86 55238 : if ((t = mpz_sizeinbase (q[i], 2) - 1) > maxbit)
87 4361 : maxbit = t;
88 : }
89 :
90 : #ifdef WANT_EXPCOST
91 : outputf (OUTPUT_ALWAYS, "Expecting %d multiplications and %d extgcds\n",
92 : 4 * (maxbit) + 6 * hamweight - 3, maxbit + 1); /* maxbit is floor(log_2(max(q_i))) */
93 : #endif
94 :
95 40533 : for (t = 0; t <= maxbit && !youpi; t++) /* Examine t-th bit of the q[i] */
96 : {
97 : /* See which values need inverting and put them into T[]. Keep number
98 : of those values in k */
99 38457 : k = 0;
100 :
101 : /* Will we have to double s at the end of this pass? If yes,
102 : schedule 2*s.y for inverting */
103 38457 : if (t < maxbit)
104 36381 : mpres_add (T[k++], s.y, s.y, modulus);
105 :
106 2561943 : for (i = 0; i < n && !youpi; i++)
107 2523486 : if (ecm_tstbit (q[i], t)) /* If q[i] & (1<<t), we'll add s to R[i] */
108 920110 : if (ecm_tstbit (flag, i)) /* Does R[i] contain a partial result yet ? */
109 : { /* If Yes: need actual point addition so */
110 864892 : mpres_sub (T[k], s.x, R[i].x, modulus); /* schedule (s.x-R[i].x) for inverting */
111 864892 : if (k > 0)
112 862816 : mpres_mul (T[k], T[k], T[k - 1], modulus);
113 864892 : k++;
114 : } /* If No: we'll simply set R[i] to s later on, nothing tbd here */
115 :
116 : /* So there are k values in need of inverting, call them v[m], 0 <= m < k. */
117 : /* Here T[m], 0 <= m < k, contains v[0]*...*v[m] */
118 :
119 : /* Put inverse of the product of all scheduled values in T[k]*/
120 38457 : if (k > 0)
121 : {
122 38457 : muls += 3 * (k - 1);
123 38457 : gcds++;
124 38457 : if (!mpres_invert (T[k], T[k - 1], modulus))
125 : {
126 : /* If a factor was found, put factor in p,
127 : flag success and bail out of loop */
128 10 : if (p != NULL)
129 10 : mpres_gcd (p, T[k - 1], modulus);
130 10 : youpi = ECM_FACTOR_FOUND_STEP2;
131 10 : break;
132 : }
133 : }
134 :
135 : /* T[k] now contains 1/(v[0]*...*v[k - 1]),
136 : T[m], 0 <= m < k, still contain v[0]*...*v[m] */
137 :
138 38447 : l = k - 1;
139 :
140 2561853 : for (i = n; i-- > 0; ) /* Go through the R[i] again, backwards */
141 2523406 : if (ecm_tstbit (q[i], t))
142 : {
143 920079 : if (ecm_tstbit (flag, i))
144 : {
145 : /* T[k] contains 1/(v[0]*...*v[l]) */
146 864861 : if (l > 0) /* need to separate the values */
147 : {
148 : /* T[l - 1] has v[0]*...*v[l-1] */
149 862785 : mpres_mul (T[l], T[l - 1], T[k], modulus); /* So T[l] now has 1/v[l] == 1/(s.x - R[i].x) */
150 862785 : mpres_sub (u, s.x, R[i].x, modulus);
151 862785 : mpres_mul (T[k], T[k], u, modulus); /* T[k] now has 1/(v[0]*...*v[l - 1]) */
152 : }
153 : else
154 : {
155 : /* T[k] contains 1/v[0] */
156 2076 : mpres_set (T[0], T[k], modulus);
157 : }
158 :
159 : /* 1/(s.x - R[i].x) is in T[l] */
160 : #ifdef WANT_ASSERT
161 : mpres_sub (u, s.x, R[i].x, modulus);
162 : mpres_mul (u, u, T[l], modulus);
163 : mpz_init(__dummy);
164 : mpres_get_z (__dummy, u, modulus);
165 : mpz_mod (__dummy, __dummy, modulus->orig_modulus);
166 : if (mpz_cmp_ui (__dummy, 1) != 0)
167 : outputf (OUTPUT_ERROR, "Error, (s.x - R[%d].x) * T[%d] == "
168 : "%Zd\n", i, l, __dummy);
169 : mpz_clear(__dummy);
170 : #endif
171 :
172 864861 : mpres_sub (u, s.y, R[i].y, modulus); /* U = y2 - y1 */
173 864861 : mpres_mul (T[l], T[l], u, modulus); /* T[l] = (y2-y1)/(x2-x1) = lambda */
174 864861 : mpres_sqr (u, T[l], modulus); /* U = lambda^2 */
175 864861 : mpres_sub (u, u, R[i].x, modulus); /* U = lambda^2 - x1 */
176 864861 : mpres_sub (R[i].x, u, s.x, modulus); /* x3 = lambda^2 - x1 - x2 */
177 864861 : mpres_sub (u, s.x, R[i].x, modulus); /* U = x2 - x3 */
178 864861 : mpres_mul (u, u, T[l], modulus); /* U = lambda*(x2 - x3) */
179 864861 : mpres_sub (R[i].y, u, s.y, modulus); /* y3 = lambda*(x2 - x3) - y2 */
180 864861 : muls += 3;
181 864861 : l--;
182 : }
183 : else /* R[i] does not contain a partial result. */
184 : {
185 55218 : mpres_set (R[i].x, s.x, modulus); /* Just set R[i] to s */
186 55218 : mpres_set (R[i].y, s.y, modulus);
187 55218 : mpz_setbit (flag, i); /* and flag it as used */
188 : }
189 : }
190 :
191 38447 : if (t < maxbit) /* Double s */
192 : {
193 : ASSERT(l==0);
194 : #ifdef WANT_ASSERT
195 : mpres_add (u, s.y, s.y, modulus);
196 : mpres_mul (u, u, T[k], modulus);
197 : mpz_init(__dummy);
198 : mpres_get_z (__dummy, u, modulus);
199 : mpz_mod (__dummy, __dummy, modulus->orig_modulus);
200 : if (mpz_cmp_ui (__dummy, 1) != 0)
201 : outputf (OUTPUT_ERROR, "Error, at t==%d, 2*s.y / (2*s.y) == %Zd\n",
202 : t, __dummy);
203 : mpz_clear(__dummy);
204 : #endif
205 :
206 : /* 1/(2*s.y) is in T[k] */
207 36371 : mpres_sqr (u, s.x, modulus); /* U = X^2 */
208 36371 : mpres_mul_ui (u, u, 3, modulus); /* U = 3*X^2 */
209 36371 : mpres_add (u, u, S->A, modulus); /* U = 3*X^2 + A */
210 36371 : mpres_mul (T[k], T[k], u, modulus); /* T = (3*X^2 + A) / (2*Y) = lambda */
211 36371 : mpres_sqr (u, T[k], modulus); /* U = lambda^2 */
212 36371 : mpres_sub (u, u, s.x, modulus); /* U = lambda^2 - X */
213 36371 : mpres_sub (u, u, s.x, modulus); /* U = lambda^2 - 2*X = s.x' */
214 36371 : mpres_sub (v, s.x, u, modulus); /* V = s.x - s.x' */
215 36371 : mpres_mul (v, v, T[k], modulus); /* V = lambda*(s.x - s.x') */
216 36371 : mpres_sub (s.y, v, s.y, modulus); /* s.y' = lambda*(s.x - s.x') - s.y */
217 36371 : mpres_set (s.x, u, modulus);
218 36371 : muls += 4;
219 : }
220 : }
221 :
222 2086 : mpres_clear (s.y, modulus);
223 2086 : mpres_clear (s.x, modulus);
224 2086 : mpz_clear (flag);
225 :
226 2086 : if (tot_muls != NULL)
227 2086 : *tot_muls += muls;
228 2086 : if (tot_gcds != NULL)
229 2086 : *tot_gcds += gcds;
230 :
231 : /* Now take inverse points (negative y-coordinate) where q[i] was < 0 */
232 57324 : for (i = 0; i < n; i++)
233 55238 : if (ecm_tstbit (signs, i))
234 : {
235 5292 : mpz_neg (R[i].y, R[i].y);
236 5292 : mpz_neg (q[i], q[i]);
237 : }
238 :
239 2086 : mpz_clear (signs);
240 2086 : return youpi;
241 : }
242 :
243 :
244 : /* Input: Points X[0]..X[(n+1)*m-1]
245 : T is used for temporary values and needs to have (n-1)*m+1 entries.
246 :
247 : Performs the following loop with only one gcdext, using Montgomery's trick:
248 : for (i=0;i<m;i++)
249 : for (j=0;j<n;j++)
250 : (x[j+(n+1)*i] : y[j+(n+1)*i]) += (x[j+1+(n+1)*i] : y[j+1+(n+1)*i])
251 :
252 : Uses one inversion and 6*n*m-3 multiplications for n*m > 0.
253 : Processes neutral (zero), identical and negative points correctly.
254 :
255 : Return factor found or not (no error can occur here).
256 : */
257 :
258 : static int
259 287573 : addWnm (mpz_t p, point *X, curve *S, mpmod_t modulus, unsigned int m,
260 : unsigned int n, mpres_t *T, unsigned long *tot_muls,
261 : unsigned long *tot_gcds)
262 : {
263 : unsigned int k, l;
264 : int i, j;
265 :
266 287573 : if (n == 0 || m == 0)
267 0 : return ECM_NO_FACTOR_FOUND;
268 :
269 287573 : k = 0;
270 4457640 : for (i = m - 1; i >= 0; i--) /* Go through the m different lists */
271 48235726 : for (j = n - 1; j >= 0; j--) /* Go through each list backwards */
272 : { /* And prepare the values to be inverted */
273 : point *X1, *X2;
274 44065659 : X1 = X + i * (n + 1) + j;
275 44065659 : X2 = X + i * (n + 1) + j + 1;
276 :
277 : /* If either element is the neutral element, nothing tbd here */
278 88131289 : if ((mpres_is_zero (X1->x, modulus) && mpres_is_zero (X1->y, modulus)) ||
279 44065639 : (mpres_is_zero (X2->x, modulus) && mpres_is_zero (X2->y, modulus)))
280 38 : continue;
281 :
282 44065621 : mpres_sub (T[k], X2->x, X1->x, modulus); /* Schedule X2.x - X1.x */
283 :
284 44065621 : if (mpres_is_zero (T[k], modulus)) /* If both x-cordinates are identical */
285 : {
286 : /* Are the points identical? Compare y coordinates: */
287 34 : mpres_sub (T[k], X2->y, X1->y, modulus);
288 34 : if (mpres_is_zero (T[k], modulus))
289 : {
290 : /* Yes, we need to double. Schedule 2*X[...].y */
291 25 : mpres_add (T[k], X1->y, X1->y, modulus);
292 : }
293 : else /* No, they are inverses. Nothing tbd here */
294 : {
295 : #ifdef WANT_ASSERT
296 : /* Check that the y coordinates are mutual negatives */
297 : mpres_add (T[k], X2->y, X1->y, modulus);
298 : ASSERT (mpres_is_zero (T[k], modulus));
299 : #endif
300 9 : continue;
301 : }
302 : }
303 :
304 44065612 : if (k > 0)
305 43778059 : mpres_mul (T[k], T[k], T[k - 1], modulus);
306 44065612 : k++;
307 : }
308 :
309 : /* v_m = X[i * (n + 1) + j] - X[i * (n + 1) + j + 1], 0 <= j < n,
310 : and m = i * n + j */
311 : /* Here T[m] = v_0 * ... * v_m, 0 <= m < k */
312 :
313 287573 : if (k > 0 && !mpres_invert (T[k], T[k - 1], modulus))
314 : {
315 11 : if (p != NULL)
316 11 : mpres_gcd (p, T[k - 1], modulus);
317 11 : if (tot_muls != NULL)
318 11 : (*tot_muls) += m * n - 1;
319 11 : if (tot_gcds != NULL)
320 11 : (*tot_gcds) ++;
321 11 : return ECM_FACTOR_FOUND_STEP2;
322 : }
323 :
324 : /* T[k] = 1/(v_0 * ... * v_m), 0 <= m < k */
325 :
326 287562 : l = k - 1;
327 :
328 4457526 : for (i = 0; (unsigned) i < m; i++)
329 48235316 : for (j = 0; (unsigned) j < n; j++)
330 : {
331 : point *X1, *X2;
332 44065352 : X1 = X + i * (n + 1) + j;
333 44065352 : X2 = X + i * (n + 1) + j + 1;
334 :
335 : /* Is X1 the neutral element? */
336 44065352 : if (mpres_is_zero (X1->x, modulus) && mpres_is_zero (X1->y, modulus))
337 : {
338 : /* Yes, set X1 to X2 */
339 29 : mpres_set (X1->x, X2->x, modulus);
340 29 : mpres_set (X1->y, X2->y, modulus);
341 29 : continue;
342 : }
343 :
344 : /* Is X2 the neutral element? If so, X1 stays the same */
345 44065323 : if (mpres_is_zero (X2->x, modulus) && mpres_is_zero (X2->y, modulus))
346 9 : continue;
347 :
348 : /* Are the x-coordinates identical? */
349 44065314 : mpres_sub (T[k + 1], X2->x, X1->x, modulus);
350 44065314 : if (mpres_is_zero (T[k + 1], modulus))
351 : {
352 : /* Are the points inverses of each other? */
353 34 : mpres_sub (T[k + 1], X2->y, X1->y, modulus);
354 34 : if (!mpres_is_zero (T[k + 1], modulus))
355 : {
356 : /* Yes. Set X1 to neutral element */
357 9 : mpres_set_ui (X1->x, 0, modulus);
358 9 : mpres_set_ui (X1->y, 0, modulus);
359 9 : continue;
360 : }
361 : /* No, we need to double. Restore T[k+1] */
362 25 : mpres_sub (T[k + 1], X2->x, X1->x, modulus);
363 : }
364 :
365 44065305 : if (l == 0)
366 287542 : mpz_set (T[0], T[k]);
367 : else
368 43777763 : mpres_mul (T[l], T[k], T[l - 1], modulus);
369 : /* T_l = 1/(v_0 * ... * v_l) * (v_0 * ... * v_{l-1}) = 1/v_l */
370 :
371 :
372 44065305 : if (mpres_is_zero (T[k + 1], modulus)) /* Identical points, so double X1 */
373 : {
374 25 : if (l > 0)
375 : {
376 6 : mpres_add (T[k + 1], X1->y, X1->y, modulus); /* T[k+1] = v_{l} */
377 6 : mpres_mul (T[k], T[k], T[k + 1], modulus);
378 : /* T_k = 1/(v_0 * ... * v_l) * v_l = 1/(v_0 * ... * v_{l-1}) */
379 : }
380 :
381 25 : mpres_sqr (T[k + 1], X1->x, modulus);
382 25 : mpres_mul_ui (T[k + 1], T[k + 1], 3, modulus);
383 25 : mpres_add (T[k + 1], T[k + 1], S->A, modulus);
384 25 : mpres_mul (T[l], T[k + 1], T[l], modulus); /* T[l] = lambda */
385 25 : mpres_sqr (T[k + 1], T[l], modulus); /* T1 = lambda^2 */
386 25 : mpres_sub (T[k + 1], T[k + 1], X1->x, modulus); /* T1 = lambda^2 - x1 */
387 25 : mpres_sub (X1->x, T[k + 1], X2->x, modulus); /* X1.x = lambda^2 - x1 - x2 = x3 */
388 25 : mpres_sub (T[k + 1], X2->x, X1->x, modulus); /* T1 = x2 - x3 */
389 25 : mpres_mul (T[k + 1], T[k + 1], T[l], modulus); /* T1 = lambda*(x2 - x3) */
390 25 : mpres_sub (X1->y, T[k + 1], X2->y, modulus); /* Y1 = lambda*(x2 - x3) - y2 = y3 */
391 : }
392 : else
393 : {
394 44065280 : if (l > 0)
395 : {
396 43777757 : mpres_mul (T[k], T[k], T[k + 1], modulus);
397 : /* T_k = 1/(v_0 * ... * v_l) * v_l = 1/(v_0 * ... * v_{l-1}) */
398 : }
399 :
400 44065280 : mpres_sub (T[k + 1], X2->y, X1->y, modulus); /* T1 = y2 - y1 */
401 44065280 : mpres_mul (T[l], T[l], T[k + 1], modulus); /* Tl = (y2 - y1) / (x2 - x1) = lambda */
402 44065280 : mpres_sqr (T[k + 1], T[l], modulus); /* T1 = lambda^2 */
403 44065280 : mpres_sub (T[k + 1], T[k + 1], X1->x, modulus); /* T1 = lambda^2 - x1 */
404 44065280 : mpres_sub (X1->x, T[k + 1], X2->x, modulus); /* X1.x = lambda^2 - x1 - x2 = x3 */
405 44065280 : mpres_sub (T[k + 1], X2->x, X1->x, modulus); /* T1 = x2 - x3 */
406 44065280 : mpres_mul (T[k + 1], T[k + 1], T[l], modulus); /* T1 = lambda*(x2 - x3) */
407 44065280 : mpres_sub (X1->y, T[k + 1], X2->y, modulus); /* Y1 = lambda*(x2 - x3) - y2 = y3 */
408 : }
409 :
410 44065305 : l--;
411 : }
412 :
413 287562 : if (tot_muls != NULL)
414 287562 : (*tot_muls) += 6 * m * n - 3;
415 287562 : if (tot_gcds != NULL)
416 287562 : (*tot_gcds) ++;
417 :
418 287562 : return ECM_NO_FACTOR_FOUND;
419 : }
420 :
421 : /* puts in F[0..dF-1] the successive values of
422 :
423 : Dickson_{S, a} (j * d2) * s where s is a point on the elliptic curve
424 :
425 : for j == 1 mod 6, j and d1 coprime.
426 : Returns non-zero iff a factor was found (then stored in f)
427 : or an error occurred.
428 : */
429 :
430 : int
431 1048 : ecm_rootsF (mpz_t f, listz_t F, root_params_t *root_params,
432 : unsigned long dF, curve *s, mpmod_t modulus)
433 : {
434 : unsigned long i;
435 1048 : unsigned long muls = 0, gcds = 0;
436 : long st;
437 1048 : int youpi = ECM_NO_FACTOR_FOUND;
438 : listz_t coeffs;
439 : ecm_roots_state_t state;
440 1048 : progression_params_t *params = &state.params; /* for less typing */
441 : mpz_t t;
442 :
443 1048 : if (dF == 0)
444 0 : return ECM_NO_FACTOR_FOUND;
445 :
446 1048 : st = cputime ();
447 :
448 : /* Relative cost of point add during init and computing roots assumed =1 */
449 1048 : init_roots_params (params, root_params->S, root_params->d1, root_params->d2,
450 : 1.0);
451 :
452 1048 : outputf (OUTPUT_DEVVERBOSE, "ecm_rootsF: state: nr = %d, dsieve = %d, "
453 : "size_fd = %d, S = %d, dickson_a = %d\n",
454 : params->nr, params->dsieve, params->size_fd, params->S,
455 : params->dickson_a);
456 :
457 : /* Init finite differences tables */
458 1048 : mpz_init (t); /* t = 0 */
459 1048 : coeffs = init_progression_coeffs (t, params->dsieve, root_params->d2,
460 : 1, 6, params->S, params->dickson_a);
461 1048 : mpz_clear (t);
462 :
463 1048 : if (coeffs == NULL) /* error */
464 : {
465 0 : youpi = ECM_ERROR;
466 0 : goto clear;
467 : }
468 :
469 : /* The highest coefficient is the same for all progressions, so set them
470 : to one for all but the first progression, later we copy the point.
471 : FIXME: can we avoid the multiplication of those points in multiplyW2n()
472 : below?
473 : */
474 9646 : for (i = params->S + 1; i < params->size_fd; i += params->S + 1)
475 8598 : mpz_set_ui (coeffs[i + params->S], 1);
476 :
477 : /* Allocate memory for fd[] and T[] */
478 :
479 1048 : state.fd = (point *) malloc (params->size_fd * sizeof (point));
480 1048 : if (state.fd == NULL)
481 : {
482 0 : youpi = ECM_ERROR;
483 0 : goto exit_ecm_rootsF;
484 : }
485 29550 : for (i = 0; i < params->size_fd; i++)
486 : {
487 28502 : outputf (OUTPUT_TRACE, "ecm_rootsF: coeffs[%d] = %Zd\n", i, coeffs[i]);
488 28502 : mpres_init (state.fd[i].x, modulus);
489 28502 : mpres_init (state.fd[i].y, modulus);
490 : }
491 :
492 1048 : state.T = (mpres_t *) malloc ((params->size_fd + 4) * sizeof (mpres_t));
493 1048 : if (state.T == NULL)
494 : {
495 0 : youpi = ECM_ERROR;
496 0 : goto ecm_rootsF_clearfdi;
497 : }
498 33742 : for (i = 0 ; i < params->size_fd + 4; i++)
499 32694 : mpres_init (state.T[i], modulus);
500 :
501 : /* Multiply fd[] = s * coeffs[] */
502 :
503 1048 : youpi = multiplyW2n (f, state.fd, s, coeffs, params->size_fd, modulus,
504 1048 : state.T[0], state.T[1], state.T + 2, &muls, &gcds);
505 1048 : if (youpi == ECM_FACTOR_FOUND_STEP2)
506 10 : outputf (OUTPUT_VERBOSE, "Found factor while computing coeff[] * X\n");
507 :
508 1048 : if (youpi == ECM_ERROR)
509 0 : goto clear;
510 :
511 : /* Copy the point corresponding to the highest coefficient of the first
512 : progression to the other progressions */
513 9646 : for (i = params->S + 1; i < params->size_fd; i += params->S + 1)
514 : {
515 8598 : mpres_set (state.fd[i + params->S].x, state.fd[params->S].x, modulus);
516 8598 : mpres_set (state.fd[i + params->S].y, state.fd[params->S].y, modulus);
517 : }
518 :
519 1048 : clear_list (coeffs, params->size_fd);
520 1048 : coeffs = NULL;
521 :
522 1048 : if (test_verbose (OUTPUT_VERBOSE))
523 : {
524 57 : unsigned int st1 = cputime ();
525 57 : outputf (OUTPUT_VERBOSE,
526 : "Initializing tables of differences for F took %ldms",
527 : elltime (st, st1));
528 57 : outputf (OUTPUT_DEVVERBOSE, ", %lu muls and %lu extgcds", muls, gcds);
529 57 : outputf (OUTPUT_VERBOSE, "\n");
530 57 : st = st1;
531 57 : muls = 0;
532 57 : gcds = 0;
533 : }
534 :
535 : /* Now for the actual calculation of the roots. */
536 :
537 1950491 : for (i = 0; i < dF && !youpi;)
538 : {
539 : /* Is this a rsieve value where we computed Dickson(j * d2) * X? */
540 1949443 : if (gcd ((unsigned long) params->rsieve,
541 1949443 : (unsigned long) params->dsieve) == 1UL)
542 : {
543 : /* Did we use every progression since the last update? */
544 1357551 : if (params->next == params->nr)
545 : {
546 : /* Yes, time to update again */
547 88697 : youpi = addWnm (f, state.fd, s, modulus, params->nr, params->S,
548 : state.T, &muls, &gcds);
549 : ASSERT(youpi != ECM_ERROR); /* no error can occur in addWnm */
550 88697 : params->next = 0;
551 88697 : if (youpi == ECM_FACTOR_FOUND_STEP2)
552 0 : outputf (OUTPUT_VERBOSE,
553 : "Found factor while computing roots of F\n");
554 : }
555 :
556 : /* Is this a j value where we want Dickson(j * d2) * X as a root? */
557 1357551 : if (gcd ((unsigned long) params->rsieve, root_params->d1)
558 : == 1UL)
559 1135590 : mpres_get_z (F[i++],
560 1135590 : state.fd[params->next * (params->S + 1)].x, modulus);
561 :
562 1357551 : params->next ++;
563 : }
564 1949443 : params->rsieve += 6;
565 : }
566 :
567 1048 : clear:
568 33742 : for (i = 0 ; i < params->size_fd + 4; i++)
569 32694 : mpres_clear (state.T[i], modulus);
570 1048 : free (state.T);
571 :
572 1048 : ecm_rootsF_clearfdi:
573 29550 : for (i = 0; i < params->size_fd; i++)
574 : {
575 28502 : mpres_clear (state.fd[i].x, modulus);
576 28502 : mpres_clear (state.fd[i].y, modulus);
577 : }
578 1048 : free (state.fd);
579 :
580 1048 : exit_ecm_rootsF:
581 1048 : if (youpi)
582 10 : return youpi; /* error or factor found */
583 :
584 1038 : outputf (OUTPUT_VERBOSE, "Computing roots of F took %ldms",
585 : elltime (st, cputime ()));
586 1038 : outputf (OUTPUT_DEVVERBOSE, ", %ld muls and %ld extgcds", muls, gcds);
587 1038 : outputf (OUTPUT_VERBOSE, "\n");
588 :
589 1038 : return ECM_NO_FACTOR_FOUND;
590 : }
591 :
592 : /* Perform the necessary initialization to allow computation of
593 :
594 : Dickson_{S, a}(s+n*d) * P , where P is a point on the elliptic curve
595 :
596 : for successive n, where Dickson_{S, a} is the degree S Dickson
597 : polynomial with parameter a. For a == 0, Dickson_{S, a} (x) = x^S.
598 :
599 : If a factor is found during the initialisation, NULL is returned and the
600 : factor in f. If an error occurred, NULL is returned and f is -1.
601 : */
602 :
603 : ecm_roots_state_t *
604 1038 : ecm_rootsG_init (mpz_t f, curve *X, root_params_t *root_params,
605 : unsigned long dF, unsigned long blocks, mpmod_t modulus)
606 : {
607 : unsigned int k, phid2;
608 1038 : unsigned long muls = 0, gcds = 0;
609 : listz_t coeffs;
610 : ecm_roots_state_t *state;
611 : progression_params_t *params; /* for less typing */
612 1038 : int youpi = 0;
613 : unsigned int T_inv;
614 : double bestnr;
615 1038 : long st = 0;
616 :
617 : ASSERT (gcd (root_params->d1, root_params->d2) == 1UL);
618 :
619 1038 : if (test_verbose (OUTPUT_VERBOSE))
620 57 : st = cputime ();
621 :
622 1038 : state = (ecm_roots_state_t *) malloc (sizeof (ecm_roots_state_t));
623 1038 : if (state == NULL)
624 : {
625 0 : mpz_set_si (f, -1);
626 0 : return NULL;
627 : }
628 1038 : params = &(state->params);
629 :
630 : /* If S < 0, use degree |S| Dickson poly, otherwise use x^S */
631 1038 : params->dickson_a = (root_params->S < 0) ? -1 : 0;
632 1038 : params->S = abs (root_params->S);
633 :
634 : /* Estimate the cost of a modular inversion (in unit of time per
635 : modular multiplication) */
636 1038 : if (modulus->repr == ECM_MOD_BASE2)
637 79 : T_inv = 18;
638 : else
639 959 : T_inv = 6;
640 :
641 : /* Guesstimate a value for the number of disjoint progressions to use */
642 1038 : bestnr = -(4. + T_inv) + sqrt(12. * (double) dF * (double) blocks *
643 1038 : (T_inv - 3.) * log (2. * root_params->d1) / log (2.) - (4. + T_inv) *
644 1038 : (4. + T_inv));
645 1038 : bestnr /= 6. * (double) (params->S) * log (2. * root_params->d1) / log (2.0);
646 :
647 1038 : outputf (OUTPUT_TRACE, "ecm_rootsG_init: bestnr = %f\n", bestnr);
648 :
649 1038 : if (bestnr < 1.)
650 50 : params->nr = 1;
651 : else
652 988 : params->nr = (unsigned int) (bestnr + .5);
653 :
654 1038 : phid2 = eulerphi (root_params->d2);
655 :
656 : /* Round up params->nr to multiple of eulerphi(d2) */
657 1038 : if (phid2 > 1)
658 963 : params->nr = ((params->nr + (phid2 - 1)) / phid2) * phid2;
659 :
660 1038 : params->size_fd = params->nr * (params->S + 1);
661 :
662 1038 : outputf (OUTPUT_DEVVERBOSE, "ecm_rootsG_init: i0=%Zd, d1=%lu, d2=%lu, "
663 1038 : "dF=%lu, blocks=%lu, S=%u, T_inv = %d, nr=%d\n", root_params->i0,
664 : root_params->d1, root_params->d2, dF, blocks, params->S,
665 : T_inv, params->nr);
666 :
667 1038 : state->X = X;
668 1038 : params->next = 0;
669 1038 : params->dsieve = 1; /* We only init progressions coprime to d2,
670 : so nothing to be skipped */
671 1038 : params->rsieve = 0;
672 :
673 1038 : coeffs = init_progression_coeffs (root_params->i0, root_params->d2,
674 1038 : root_params->d1, params->nr / phid2,
675 : 1, params->S, params->dickson_a);
676 :
677 1038 : if (coeffs == NULL) /* error */
678 : {
679 0 : free (state);
680 0 : mpz_set_si (f, -1);
681 0 : return NULL;
682 : }
683 :
684 1038 : state->fd = (point *) malloc (params->size_fd * sizeof (point));
685 1038 : if (state->fd == NULL)
686 : {
687 0 : clear_list (coeffs, params->size_fd);
688 0 : free (state);
689 0 : mpz_set_si (f, -1);
690 0 : return NULL;
691 : }
692 27774 : for (k = 0; k < params->size_fd; k++)
693 : {
694 26736 : mpres_init (state->fd[k].x, modulus);
695 26736 : mpres_init (state->fd[k].y, modulus);
696 : }
697 :
698 1038 : state->size_T = params->size_fd + 4;
699 1038 : state->T = (mpres_t *) malloc (state->size_T * sizeof (mpres_t));
700 1038 : if (state->T == NULL)
701 : {
702 0 : for (k = 0; k < params->size_fd; k++)
703 : {
704 0 : mpres_clear (state->fd[k].x, modulus);
705 0 : mpres_clear (state->fd[k].y, modulus);
706 : }
707 0 : clear_list (coeffs, params->size_fd);
708 0 : free (state);
709 0 : mpz_set_si (f, -1);
710 0 : return NULL;
711 : }
712 31926 : for (k = 0; k < state->size_T; k++)
713 30888 : mpres_init (state->T[k], modulus);
714 :
715 10918 : for (k = params->S + 1; k < params->size_fd; k += params->S + 1)
716 9880 : mpz_set_ui (coeffs[k + params->S], 1);
717 :
718 1038 : if (test_verbose (OUTPUT_TRACE))
719 246 : for (k = 0; k < params->size_fd; k++)
720 236 : outputf (OUTPUT_TRACE, "ecm_rootsG_init: coeffs[%d] == %Zd\n",
721 236 : k, coeffs[k]);
722 :
723 1038 : youpi = multiplyW2n (f, state->fd, X, coeffs, params->size_fd, modulus,
724 1038 : state->T[0], state->T[1], state->T + 2, &muls, &gcds);
725 1038 : if (youpi == ECM_ERROR)
726 0 : mpz_set_si (f, -1); /* fall through */
727 :
728 10918 : for (k = params->S + 1; k < params->size_fd; k += params->S + 1)
729 : {
730 9880 : mpres_set (state->fd[k + params->S].x, state->fd[params->S].x, modulus);
731 9880 : mpres_set (state->fd[k + params->S].y, state->fd[params->S].y, modulus);
732 : }
733 :
734 1038 : clear_list (coeffs, params->size_fd);
735 1038 : coeffs = NULL;
736 :
737 1038 : if (youpi != ECM_NO_FACTOR_FOUND) /* factor found or error */
738 : {
739 0 : if (youpi == ECM_FACTOR_FOUND_STEP2)
740 0 : outputf (OUTPUT_VERBOSE, "Found factor while computing fd[]\n");
741 :
742 0 : ecm_rootsG_clear (state, modulus);
743 :
744 : /* Signal that a factor was found, or an error occurred (f=-1) */
745 0 : state = NULL;
746 : }
747 : else
748 : {
749 1038 : if (test_verbose (OUTPUT_VERBOSE))
750 : {
751 57 : st = elltime (st, cputime ());
752 57 : outputf (OUTPUT_VERBOSE,
753 : "Initializing table of differences for G took %ldms", st);
754 57 : outputf (OUTPUT_DEVVERBOSE, ", %lu muls and %lu extgcds",
755 : muls, gcds);
756 57 : outputf (OUTPUT_VERBOSE, "\n");
757 : }
758 : }
759 :
760 1038 : return state;
761 : }
762 :
763 : void
764 1038 : ecm_rootsG_clear (ecm_roots_state_t *state, ATTRIBUTE_UNUSED mpmod_t modulus)
765 : {
766 : unsigned int k;
767 :
768 27774 : for (k = 0; k < state->params.size_fd; k++)
769 : {
770 26736 : mpres_clear (state->fd[k].x, modulus);
771 26736 : mpres_clear (state->fd[k].y, modulus);
772 : }
773 1038 : free (state->fd);
774 :
775 1038 : if(state->size_T != 0){
776 31926 : for (k = 0; k < state->size_T; k++)
777 30888 : mpres_clear (state->T[k], modulus);
778 1038 : free (state->T);
779 : }
780 :
781 1038 : free (state);
782 1038 : }
783 :
784 : /* Puts in G the successive values of
785 :
786 : Dickson_{S, a}(s+j*k) P
787 :
788 : where P is a point on the elliptic curve,
789 : 0<= j <= dF-1, k is the 'd' value from ecm_rootsG_init()
790 : and s is the 's' value of ecm_rootsG_init() or where a previous
791 : call to ecm_rootsG has left off.
792 :
793 : Returns non-zero iff a factor was found (then stored in f).
794 : Cannot return an error.
795 : */
796 :
797 : int
798 3362 : ecm_rootsG (mpz_t f, listz_t G, unsigned long dF, ecm_roots_state_t *state,
799 : mpmod_t modulus)
800 : {
801 : unsigned long i;
802 3362 : unsigned long muls = 0, gcds = 0;
803 3362 : int youpi = ECM_NO_FACTOR_FOUND;
804 : long st;
805 3362 : point *fd = state->fd; /* to save typing */
806 3362 : progression_params_t *params = &(state->params); /* for less typing */
807 :
808 3362 : st = cputime ();
809 :
810 3362 : outputf (OUTPUT_TRACE, "ecm_rootsG: dF = %lu, state: nr = %u, next = %u, "
811 : "S = %u, dsieve = %u, rsieve = %u,\n\tdickson_a = %d\n",
812 : dF, params->nr, params->next, params->S, params->dsieve,
813 : params->rsieve, params->dickson_a);
814 :
815 2830091 : for (i = 0; i < dF;)
816 : {
817 : /* Did we use every progression since the last update? */
818 2826740 : if (params->next == params->nr)
819 : {
820 : /* Yes, time to update again */
821 198876 : youpi = addWnm (f, fd, state->X, modulus, params->nr, params->S,
822 : state->T, &muls, &gcds);
823 : ASSERT(youpi != ECM_ERROR); /* no error can occur in addWnm */
824 198876 : params->next = 0;
825 :
826 198876 : if (youpi == ECM_FACTOR_FOUND_STEP2)
827 : {
828 11 : outputf (OUTPUT_VERBOSE, "Found factor while computing G[]\n");
829 11 : break;
830 : }
831 : }
832 :
833 : /* Is this a root we should skip? (Take only if gcd == 1) */
834 2826729 : if (gcd ((unsigned long) params->rsieve,
835 2826729 : (unsigned long) params->dsieve) == 1UL)
836 : {
837 2826729 : mpres_get_z (G[i++], (fd + params->next * (params->S + 1))->x,
838 : modulus);
839 2826729 : outputf (OUTPUT_TRACE,
840 : "ecm_rootsG: storing d1*%u*X = %Zd in G[%lu]\n",
841 2826729 : params->rsieve, G[i - 1], i);
842 : }
843 :
844 2826729 : params->next ++;
845 2826729 : params->rsieve ++;
846 : }
847 :
848 3362 : outputf (OUTPUT_VERBOSE, "Computing roots of G took %ldms",
849 : elltime (st, cputime ()));
850 3362 : outputf (OUTPUT_DEVVERBOSE, ", %lu muls and %lu extgcds", muls, gcds);
851 3362 : outputf (OUTPUT_VERBOSE, "\n");
852 :
853 3362 : return youpi;
854 : }
|
