Line data Source code
1 : /* Elliptic Curve Method: toplevel and stage 1 routines.
2 :
3 : Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011,
4 : 2012, 2016 Paul Zimmermann, Alexander Kruppa, Cyril Bouvier, David Cleaver.
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 <stdio.h>
24 : #include <stdlib.h>
25 : #include <string.h>
26 : #include "ecm-impl.h"
27 : #include "getprime_r.h"
28 : #include <math.h>
29 :
30 : #ifdef HAVE_ADDLAWS
31 : #include "addlaws.h"
32 : #endif
33 :
34 : #ifdef HAVE_LIMITS_H
35 : # include <limits.h>
36 : #else
37 : # define ULONG_MAX __GMP_ULONG_MAX
38 : #endif
39 :
40 : #ifdef TIMING_CRT
41 : extern int mpzspv_from_mpzv_slow_time, mpzspv_to_mpzv_time,
42 : mpzspv_normalise_time;
43 : #endif
44 :
45 : /* the following factor takes into account the smaller expected smoothness
46 : for Montgomery's curves (batch mode) with respect to Suyama's curves */
47 : /* For param 1 we use A=4d-2 with d a square (see main.c). In that
48 : case, Cyril Bouvier and Razvan Barbulescu have shown that the average
49 : expected torsion is that of a generic Suyama curve multiplied by the
50 : constant 2^(1/3)/(3*3^(1/128)) */
51 : #define EXTRA_SMOOTHNESS_SQUARE 0.416384512396064
52 : /* For A=4d-2 (param 3) for d a random integer, the average expected torsion
53 : is that of a generic Suyama curve multiplied by the constant
54 : 1/(3*3^(1/128)) */
55 : #define EXTRA_SMOOTHNESS_32BITS_D 0.330484606500389
56 : /******************************************************************************
57 : * *
58 : * Elliptic Curve Method *
59 : * *
60 : ******************************************************************************/
61 :
62 : void duplicate (mpres_t, mpres_t, mpres_t, mpres_t, mpmod_t, mpres_t,
63 : mpres_t, mpres_t, mpres_t) ATTRIBUTE_HOT;
64 : void add3 (mpres_t, mpres_t, mpres_t, mpres_t, mpres_t, mpres_t, mpres_t,
65 : mpres_t, mpmod_t, mpres_t, mpres_t, mpres_t) ATTRIBUTE_HOT;
66 :
67 : #define mpz_mulmod5(r,s1,s2,m,t) { mpz_mul(t,s1,s2); mpz_mod(r, t, m); }
68 :
69 : /* switch from Montgomery's form g*y^2 = x^3 + a*x^2 + x
70 : to Weierstrass' form Y^2 = X^3 + A*X + B
71 : by change of variables x -> g*X-a/3, y -> g*Y.
72 : We have A = (3-a^2)/(3g^2), X = (3x+a)/(3g), Y = y/g.
73 : If a factor is found during the modular inverse, returns
74 : ECM_FACTOR_FOUND_STEP1 and the factor in f, otherwise returns
75 : ECM_NO_FACTOR_FOUND.
76 :
77 : The input value of y is the y0 value in the initial equation:
78 : g*y0^2 = x0^3 + a*x0^2 + x0.
79 : */
80 : int
81 978 : montgomery_to_weierstrass (mpz_t f, mpres_t x, mpres_t y, mpres_t A, mpmod_t n)
82 : {
83 : mpres_t g;
84 :
85 978 : mpres_init (g, n);
86 978 : mpres_add (g, x, A, n);
87 978 : mpres_mul (g, g, x, n);
88 978 : mpres_add_ui (g, g, 1, n);
89 978 : mpres_mul (g, g, x, n); /* g = x^3+a*x^2+x (y=1) */
90 978 : mpres_mul_ui (y, g, 3, n);
91 978 : mpres_mul (y, y, g, n); /* y = 3g^2 */
92 978 : if (!mpres_invert (y, y, n)) /* y = 1/(3g^2) temporarily */
93 : {
94 0 : mpres_gcd (f, y, n);
95 0 : mpres_clear (g, n);
96 0 : return ECM_FACTOR_FOUND_STEP1;
97 : }
98 :
99 : /* update x */
100 978 : mpres_mul_ui (x, x, 3, n); /* 3x */
101 978 : mpres_add (x, x, A, n); /* 3x+a */
102 978 : mpres_mul (x, x, g, n); /* (3x+a)*g */
103 978 : mpres_mul (x, x, y, n); /* (3x+a)/(3g) */
104 :
105 : /* update A */
106 978 : mpres_sqr (A, A, n); /* a^2 */
107 978 : mpres_sub_ui (A, A, 3, n);
108 978 : mpres_neg (A, A, n); /* 3-a^2 */
109 978 : mpres_mul (A, A, y, n); /* (3-a^2)/(3g^2) */
110 :
111 : /* update y */
112 978 : mpres_mul_ui (g, g, 3, n); /* 3g */
113 978 : mpres_mul (y, y, g, n); /* (3g)/(3g^2) = 1/g */
114 :
115 978 : mpres_clear (g, n);
116 978 : return ECM_NO_FACTOR_FOUND;
117 : }
118 :
119 : /* adds Q=(x2:z2) and R=(x1:z1) and puts the result in (x3:z3),
120 : using 6 muls (4 muls and 2 squares), and 6 add/sub.
121 : One assumes that Q-R=P or R-Q=P where P=(x:z).
122 : - n : number to factor
123 : - u, v, w : auxiliary variables
124 : Modifies: x3, z3, u, v, w.
125 : (x3,z3) may be identical to (x2,z2) and to (x,z)
126 : */
127 : void
128 474034543 : add3 (mpres_t x3, mpres_t z3, mpres_t x2, mpres_t z2, mpres_t x1, mpres_t z1,
129 : mpres_t x, mpres_t z, mpmod_t n, mpres_t u, mpres_t v, mpres_t w)
130 : {
131 474034543 : mpres_sub (u, x2, z2, n);
132 474034543 : mpres_add (v, x1, z1, n); /* u = x2-z2, v = x1+z1 */
133 :
134 474034543 : mpres_mul (u, u, v, n); /* u = (x2-z2)*(x1+z1) */
135 :
136 474034543 : mpres_add (w, x2, z2, n);
137 474034543 : mpres_sub (v, x1, z1, n); /* w = x2+z2, v = x1-z1 */
138 :
139 474034543 : mpres_mul (v, w, v, n); /* v = (x2+z2)*(x1-z1) */
140 :
141 474034543 : mpres_add (w, u, v, n); /* w = 2*(x1*x2-z1*z2) */
142 474034543 : mpres_sub (v, u, v, n); /* v = 2*(x2*z1-x1*z2) */
143 :
144 474034543 : mpres_sqr (w, w, n); /* w = 4*(x1*x2-z1*z2)^2 */
145 474034543 : mpres_sqr (v, v, n); /* v = 4*(x2*z1-x1*z2)^2 */
146 :
147 474034543 : if (x == x3) /* same variable: in-place variant */
148 : {
149 : /* x3 <- w * z mod n
150 : z3 <- x * v mod n */
151 1494556 : mpres_mul (z3, w, z, n);
152 1494556 : mpres_mul (x3, x, v, n);
153 1494556 : mpres_swap (x3, z3, n);
154 : }
155 : else
156 : {
157 472539987 : mpres_mul (x3, w, z, n); /* x3 = 4*z*(x1*x2-z1*z2)^2 mod n */
158 472539987 : mpres_mul (z3, x, v, n); /* z3 = 4*x*(x2*z1-x1*z2)^2 mod n */
159 : }
160 : /* mul += 6; */
161 474034543 : }
162 :
163 : /* computes 2P=(x2:z2) from P=(x1:z1), with 5 muls (3 muls and 2 squares)
164 : and 4 add/sub.
165 : - n : number to factor
166 : - b : (a+2)/4 mod n
167 : - t, u, v, w : auxiliary variables
168 : */
169 : void
170 45066559 : duplicate (mpres_t x2, mpres_t z2, mpres_t x1, mpres_t z1, mpmod_t n,
171 : mpres_t b, mpres_t u, mpres_t v, mpres_t w)
172 : {
173 45066559 : mpres_add (u, x1, z1, n);
174 45066559 : mpres_sqr (u, u, n); /* u = (x1+z1)^2 mod n */
175 45066559 : mpres_sub (v, x1, z1, n);
176 45066559 : mpres_sqr (v, v, n); /* v = (x1-z1)^2 mod n */
177 45066559 : mpres_mul (x2, u, v, n); /* x2 = u*v = (x1^2 - z1^2)^2 mod n */
178 45066559 : mpres_sub (w, u, v, n); /* w = u-v = 4*x1*z1 */
179 45066559 : mpres_mul (u, w, b, n); /* u = w*b = ((A+2)/4*(4*x1*z1)) mod n */
180 45066559 : mpres_add (u, u, v, n); /* u = (x1-z1)^2+(A+2)/4*(4*x1*z1) */
181 45066559 : mpres_mul (z2, w, u, n); /* z2 = ((4*x1*z1)*((x1-z1)^2+(A+2)/4*(4*x1*z1))) mod n */
182 45066559 : }
183 :
184 : /* multiply P=(x:z) by e and puts the result in (x:z). */
185 : void
186 975 : ecm_mul (mpres_t x, mpres_t z, mpz_t e, mpmod_t n, mpres_t b)
187 : {
188 : size_t l;
189 975 : int negated = 0;
190 : mpres_t x0, z0, x1, z1, u, v, w;
191 :
192 : /* In Montgomery coordinates, the point at infinity is (0::0) */
193 975 : if (mpz_sgn (e) == 0)
194 : {
195 0 : mpz_set_ui (x, 0);
196 0 : mpz_set_ui (z, 0);
197 0 : return;
198 : }
199 :
200 : /* The negative of a point (x:y:z) is (x:-y:u). Since we do not compute
201 : y, e*(x::z) == (-e)*(x::z). */
202 975 : if (mpz_sgn (e) < 0)
203 : {
204 0 : negated = 1;
205 0 : mpz_neg (e, e);
206 : }
207 :
208 975 : if (mpz_cmp_ui (e, 1) == 0)
209 850 : goto ecm_mul_end;
210 :
211 125 : mpres_init (x0, n);
212 125 : mpres_init (z0, n);
213 125 : mpres_init (x1, n);
214 125 : mpres_init (z1, n);
215 125 : mpres_init (u, n);
216 125 : mpres_init (v, n);
217 125 : mpres_init (w, n);
218 :
219 125 : l = mpz_sizeinbase (e, 2) - 1; /* l >= 1 */
220 :
221 125 : mpres_set (x0, x, n);
222 125 : mpres_set (z0, z, n);
223 125 : duplicate (x1, z1, x0, z0, n, b, u, v, w);
224 :
225 : /* invariant: (P1,P0) = ((k+1)P, kP) where k = floor(e/2^l) */
226 :
227 4117 : while (l-- > 0)
228 : {
229 3992 : if (ecm_tstbit (e, l)) /* k, k+1 -> 2k+1, 2k+2 */
230 : {
231 1957 : add3 (x0, z0, x0, z0, x1, z1, x, z, n, u, v, w); /* 2k+1 */
232 1957 : duplicate (x1, z1, x1, z1, n, b, u, v, w); /* 2k+2 */
233 : }
234 : else /* k, k+1 -> 2k, 2k+1 */
235 : {
236 2035 : add3 (x1, z1, x1, z1, x0, z0, x, z, n, u, v, w); /* 2k+1 */
237 2035 : duplicate (x0, z0, x0, z0, n, b, u, v, w); /* 2k */
238 : }
239 : }
240 :
241 125 : mpres_set (x, x0, n);
242 125 : mpres_set (z, z0, n);
243 :
244 125 : mpres_clear (x0, n);
245 125 : mpres_clear (z0, n);
246 125 : mpres_clear (x1, n);
247 125 : mpres_clear (z1, n);
248 125 : mpres_clear (u, n);
249 125 : mpres_clear (v, n);
250 125 : mpres_clear (w, n);
251 :
252 975 : ecm_mul_end:
253 :
254 : /* Undo negation to avoid changing the caller's e value */
255 975 : if (negated)
256 0 : mpz_neg (e, e);
257 : }
258 :
259 : #define ADD 6.0 /* number of multiplications in an addition */
260 : #define DUP 5.0 /* number of multiplications in a duplicate */
261 :
262 : /* returns the number of modular multiplications for computing
263 : V_n from V_r * V_{n-r} - V_{n-2r}.
264 : ADD is the cost of an addition
265 : DUP is the cost of a duplicate
266 : */
267 : static double
268 174398870 : lucas_cost (ecm_uint n, double v)
269 : {
270 : ecm_uint d, e, r;
271 : double c; /* cost */
272 :
273 174398870 : d = n;
274 174398870 : r = (ecm_uint) ((double) d * v + 0.5);
275 174398870 : if (r >= n)
276 0 : return (ADD * (double) n);
277 174398870 : d = n - r;
278 174398870 : e = 2 * r - n;
279 174398870 : c = DUP + ADD; /* initial duplicate and final addition */
280 4501529164 : while (d != e)
281 : {
282 4327130294 : if (d < e)
283 : {
284 2913599944 : r = d;
285 2913599944 : d = e;
286 2913599944 : e = r;
287 : }
288 4327130294 : if (d - e <= e / 4 && ((d + e) % 3) == 0)
289 : { /* condition 1 */
290 98320044 : d = (2 * d - e) / 3;
291 98320044 : e = (e - d) / 2;
292 98320044 : c += 3.0 * ADD; /* 3 additions */
293 : }
294 4228810250 : else if (d - e <= e / 4 && (d - e) % 6 == 0)
295 : { /* condition 2 */
296 16810291 : d = (d - e) / 2;
297 16810291 : c += ADD + DUP; /* one addition, one duplicate */
298 : }
299 4211999959 : else if ((d + 3) / 4 <= e)
300 : { /* condition 3 */
301 3801244637 : d -= e;
302 3801244637 : c += ADD; /* one addition */
303 : }
304 410755322 : else if ((d + e) % 2 == 0)
305 : { /* condition 4 */
306 180017362 : d = (d - e) / 2;
307 180017362 : c += ADD + DUP; /* one addition, one duplicate */
308 : }
309 : /* now d+e is odd */
310 230737960 : else if (d % 2 == 0)
311 : { /* condition 5 */
312 152315821 : d /= 2;
313 152315821 : c += ADD + DUP; /* one addition, one duplicate */
314 : }
315 : /* now d is odd and e is even */
316 78422139 : else if (d % 3 == 0)
317 : { /* condition 6 */
318 32058854 : d = d / 3 - e;
319 32058854 : c += 3.0 * ADD + DUP; /* three additions, one duplicate */
320 : }
321 46363285 : else if ((d + e) % 3 == 0)
322 : { /* condition 7 */
323 29927311 : d = (d - 2 * e) / 3;
324 29927311 : c += 3.0 * ADD + DUP; /* three additions, one duplicate */
325 : }
326 16435974 : else if ((d - e) % 3 == 0)
327 : { /* condition 8 */
328 5057068 : d = (d - e) / 3;
329 5057068 : c += 3.0 * ADD + DUP; /* three additions, one duplicate */
330 : }
331 : else /* necessarily e is even: catches all cases */
332 : { /* condition 9 */
333 11378906 : e /= 2;
334 11378906 : c += ADD + DUP; /* one addition, one duplicate */
335 : }
336 : }
337 :
338 174398870 : return c;
339 : }
340 :
341 :
342 : /* computes kP from P=(xA:zA) and puts the result in (xA:zA). Assumes k>2.
343 : WARNING! The calls to add3() assume that the two input points are distinct,
344 : which is not neccessarily satisfied. The result can be that in rare cases
345 : the point at infinity (z==0) results when it shouldn't. A test case is
346 : echo 33554520197234177 | ./ecm -sigma 2046841451 373 1
347 : which finds the prime even though it shouldn't (23^2=529 divides order).
348 : This is not a problem for ECM since at worst we'll find a factor we
349 : shouldn't have found. For other purposes (i.e. primality proving) this
350 : would have to be fixed first.
351 : */
352 :
353 : static void
354 17439887 : prac (mpres_t xA, mpres_t zA, ecm_uint k, mpmod_t n, mpres_t b,
355 : mpres_t u, mpres_t v, mpres_t w, mpres_t xB, mpres_t zB, mpres_t xC,
356 : mpres_t zC, mpres_t xT, mpres_t zT, mpres_t xT2, mpres_t zT2)
357 : {
358 17439887 : ecm_uint d, e, r, i = 0, nv;
359 : double c, cmin;
360 : __mpz_struct *tmp;
361 : #define NV 10
362 : /* 1/val[0] = the golden ratio (1+sqrt(5))/2, and 1/val[i] for i>0
363 : is the real number whose continued fraction expansion is all 1s
364 : except for a 2 in i+1-st place */
365 : static double val[NV] =
366 : { 0.61803398874989485, 0.72360679774997897, 0.58017872829546410,
367 : 0.63283980608870629, 0.61242994950949500, 0.62018198080741576,
368 : 0.61721461653440386, 0.61834711965622806, 0.61791440652881789,
369 : 0.61807966846989581};
370 :
371 : /* for small n, it makes no sense to try 10 different Lucas chains */
372 17439887 : nv = mpz_size ((mpz_ptr) n);
373 17439887 : if (nv > NV)
374 17439887 : nv = NV;
375 :
376 17439887 : if (nv > 1)
377 : {
378 : /* chooses the best value of v */
379 191838757 : for (d = 0, cmin = ADD * (double) k; d < nv; d++)
380 : {
381 174398870 : c = lucas_cost (k, val[d]);
382 174398870 : if (c < cmin)
383 : {
384 42406179 : cmin = c;
385 42406179 : i = d;
386 : }
387 : }
388 : }
389 :
390 17439887 : d = k;
391 17439887 : r = (ecm_uint) ((double) d * val[i] + 0.5);
392 :
393 : /* first iteration always begins by Condition 3, then a swap */
394 17439887 : d = k - r;
395 17439887 : e = 2 * r - k;
396 17439887 : mpres_set (xB, xA, n);
397 17439887 : mpres_set (zB, zA, n); /* B=A */
398 17439887 : mpres_set (xC, xA, n);
399 17439887 : mpres_set (zC, zA, n); /* C=A */
400 17439887 : duplicate (xA, zA, xA, zA, n, b, u, v, w); /* A = 2*A */
401 463036741 : while (d != e)
402 : {
403 445596854 : if (d < e)
404 : {
405 327365332 : r = d;
406 327365332 : d = e;
407 327365332 : e = r;
408 327365332 : mpres_swap (xA, xB, n);
409 327365332 : mpres_swap (zA, zB, n);
410 : }
411 : /* do the first line of Table 4 whose condition qualifies */
412 445596854 : if (d - e <= e / 4 && ((d + e) % 3) == 0)
413 : { /* condition 1 */
414 4251087 : d = (2 * d - e) / 3;
415 4251087 : e = (e - d) / 2;
416 4251087 : add3 (xT, zT, xA, zA, xB, zB, xC, zC, n, u, v, w); /* T = f(A,B,C) */
417 4251087 : add3 (xT2, zT2, xT, zT, xA, zA, xB, zB, n, u, v, w); /* T2 = f(T,A,B) */
418 4251087 : add3 (xB, zB, xB, zB, xT, zT, xA, zA, n, u, v, w); /* B = f(B,T,A) */
419 4251087 : mpres_swap (xA, xT2, n);
420 4251087 : mpres_swap (zA, zT2, n); /* swap A and T2 */
421 : }
422 441345767 : else if (d - e <= e / 4 && (d - e) % 6 == 0)
423 : { /* condition 2 */
424 121314 : d = (d - e) / 2;
425 121314 : add3 (xB, zB, xA, zA, xB, zB, xC, zC, n, u, v, w); /* B = f(A,B,C) */
426 121314 : duplicate (xA, zA, xA, zA, n, b, u, v, w); /* A = 2*A */
427 : }
428 441224453 : else if ((d + 3) / 4 <= e)
429 : { /* condition 3 */
430 413741080 : d -= e;
431 413741080 : add3 (xT, zT, xB, zB, xA, zA, xC, zC, n, u, v, w); /* T = f(B,A,C) */
432 : /* circular permutation (B,T,C) */
433 413741080 : tmp = xB;
434 413741080 : xB = xT;
435 413741080 : xT = xC;
436 413741080 : xC = tmp;
437 413741080 : tmp = zB;
438 413741080 : zB = zT;
439 413741080 : zT = zC;
440 413741080 : zC = tmp;
441 : }
442 27483373 : else if ((d + e) % 2 == 0)
443 : { /* condition 4 */
444 19842261 : d = (d - e) / 2;
445 19842261 : add3 (xB, zB, xB, zB, xA, zA, xC, zC, n, u, v, w); /* B = f(B,A,C) */
446 19842261 : duplicate (xA, zA, xA, zA, n, b, u, v, w); /* A = 2*A */
447 : }
448 : /* now d+e is odd */
449 7641112 : else if (d % 2 == 0)
450 : { /* condition 5 */
451 6398731 : d /= 2;
452 6398731 : add3 (xC, zC, xC, zC, xA, zA, xB, zB, n, u, v, w); /* C = f(C,A,B) */
453 6398731 : duplicate (xA, zA, xA, zA, n, b, u, v, w); /* A = 2*A */
454 : }
455 : /* now d is odd, e is even */
456 1242381 : else if (d % 3 == 0)
457 : { /* condition 6 */
458 997070 : d = d / 3 - e;
459 997070 : duplicate (xT, zT, xA, zA, n, b, u, v, w); /* T = 2*A */
460 997070 : add3 (xT2, zT2, xA, zA, xB, zB, xC, zC, n, u, v, w); /* T2 = f(A,B,C) */
461 997070 : add3 (xA, zA, xT, zT, xA, zA, xA, zA, n, u, v, w); /* A = f(T,A,A) */
462 997070 : add3 (xT, zT, xT, zT, xT2, zT2, xC, zC, n, u, v, w); /* T = f(T,T2,C) */
463 : /* circular permutation (C,B,T) */
464 997070 : tmp = xC;
465 997070 : xC = xB;
466 997070 : xB = xT;
467 997070 : xT = tmp;
468 997070 : tmp = zC;
469 997070 : zC = zB;
470 997070 : zB = zT;
471 997070 : zT = tmp;
472 : }
473 245311 : else if ((d + e) % 3 == 0)
474 : { /* condition 7 */
475 245301 : d = (d - 2 * e) / 3;
476 245301 : add3 (xT, zT, xA, zA, xB, zB, xC, zC, n, u, v, w); /* T = f(A,B,C) */
477 245301 : add3 (xB, zB, xT, zT, xA, zA, xB, zB, n, u, v, w); /* B = f(T,A,B) */
478 245301 : duplicate (xT, zT, xA, zA, n, b, u, v, w);
479 245301 : add3 (xA, zA, xA, zA, xT, zT, xA, zA, n, u, v, w); /* A = 3*A */
480 : }
481 10 : else if ((d - e) % 3 == 0)
482 : { /* condition 8 */
483 10 : d = (d - e) / 3;
484 10 : add3 (xT, zT, xA, zA, xB, zB, xC, zC, n, u, v, w); /* T = f(A,B,C) */
485 10 : add3 (xC, zC, xC, zC, xA, zA, xB, zB, n, u, v, w); /* C = f(A,C,B) */
486 10 : mpres_swap (xB, xT, n);
487 10 : mpres_swap (zB, zT, n); /* swap B and T */
488 10 : duplicate (xT, zT, xA, zA, n, b, u, v, w);
489 10 : add3 (xA, zA, xA, zA, xT, zT, xA, zA, n, u, v, w); /* A = 3*A */
490 : }
491 : else /* necessarily e is even here */
492 : { /* condition 9 */
493 0 : e /= 2;
494 0 : add3 (xC, zC, xC, zC, xB, zB, xA, zA, n, u, v, w); /* C = f(C,B,A) */
495 0 : duplicate (xB, zB, xB, zB, n, b, u, v, w); /* B = 2*B */
496 : }
497 : }
498 :
499 17439887 : add3 (xA, zA, xA, zA, xB, zB, xC, zC, n, u, v, w);
500 :
501 : ASSERT(d == 1);
502 17439887 : }
503 :
504 : /* Input: x is initial point
505 : A is curve parameter in Montgomery's form:
506 : g*y^2*z = x^3 + a*x^2*z + x*z^2
507 : n is the number to factor
508 : B1 is the stage 1 bound
509 : Output: If a factor is found, it is returned in f.
510 : Otherwise, x contains the x-coordinate of the point computed
511 : in stage 1 (with z coordinate normalized to 1).
512 : B1done is set to B1 if stage 1 completed normally,
513 : or to the largest prime processed if interrupted, but never
514 : to a smaller value than B1done was upon function entry.
515 : Return value: ECM_FACTOR_FOUND_STEP1 if a factor, otherwise
516 : ECM_NO_FACTOR_FOUND
517 : */
518 : static int
519 975 : ecm_stage1 (mpz_t f, mpres_t x, mpres_t A, mpmod_t n, double B1,
520 : double *B1done, mpz_t go, int (*stop_asap)(void),
521 : char *chkfilename)
522 : {
523 : mpres_t b, z, u, v, w, xB, zB, xC, zC, xT, zT, xT2, zT2;
524 : uint64_t p, r, last_chkpnt_p;
525 975 : int ret = ECM_NO_FACTOR_FOUND;
526 : long last_chkpnt_time;
527 : prime_info_t prime_info;
528 :
529 975 : prime_info_init (prime_info);
530 :
531 975 : mpres_init (b, n);
532 975 : mpres_init (z, n);
533 975 : mpres_init (u, n);
534 975 : mpres_init (v, n);
535 975 : mpres_init (w, n);
536 975 : mpres_init (xB, n);
537 975 : mpres_init (zB, n);
538 975 : mpres_init (xC, n);
539 975 : mpres_init (zC, n);
540 975 : mpres_init (xT, n);
541 975 : mpres_init (zT, n);
542 975 : mpres_init (xT2, n);
543 975 : mpres_init (zT2, n);
544 :
545 975 : last_chkpnt_time = cputime ();
546 :
547 975 : mpres_set_ui (z, 1, n);
548 :
549 975 : mpres_add_ui (b, A, 2, n);
550 975 : mpres_div_2exp (b, b, 2, n); /* b == (A0+2)*B/4 */
551 :
552 : /* preload group order */
553 975 : if (go != NULL)
554 975 : ecm_mul (x, z, go, n, b);
555 :
556 : /* prac() wants multiplicands > 2 */
557 11969 : for (r = 2; r <= B1; r *= 2)
558 10994 : if (r > *B1done)
559 10994 : duplicate (x, z, x, z, n, b, u, v, w);
560 :
561 : /* We'll do 3 manually, too (that's what ecm4 did..) */
562 7849 : for (r = 3; r <= B1; r *= 3)
563 6874 : if (r > *B1done)
564 : {
565 6874 : duplicate (xB, zB, x, z, n, b, u, v, w);
566 6874 : add3 (x, z, x, z, xB, zB, x, z, n, u, v, w);
567 : }
568 :
569 975 : last_chkpnt_p = 3;
570 975 : p = getprime_mt (prime_info); /* Puts 3 into p. Next call gives 5 */
571 17391227 : for (p = getprime_mt (prime_info); p <= B1; p = getprime_mt (prime_info))
572 : {
573 34830199 : for (r = p; r <= B1; r *= p)
574 17439887 : if (r > *B1done)
575 17439887 : prac (x, z, (ecm_uint) p, n, b, u, v, w, xB, zB, xC, zC, xT,
576 : zT, xT2, zT2);
577 :
578 17390312 : if (mpres_is_zero (z, n))
579 : {
580 60 : outputf (OUTPUT_VERBOSE, "Reached point at infinity, %.0f divides "
581 : "group orders\n", p);
582 60 : break;
583 : }
584 :
585 17390252 : if (stop_asap != NULL && (*stop_asap) ())
586 : {
587 0 : outputf (OUTPUT_NORMAL, "Interrupted at prime %.0f\n", p);
588 0 : break;
589 : }
590 :
591 17390252 : if (chkfilename != NULL && p > last_chkpnt_p + 10000 &&
592 0 : elltime (last_chkpnt_time, cputime ()) > CHKPNT_PERIOD)
593 : {
594 0 : writechkfile (chkfilename, ECM_ECM, MAX(p, *B1done), n, A, x, NULL, z);
595 0 : last_chkpnt_p = p;
596 0 : last_chkpnt_time = cputime ();
597 : }
598 : }
599 :
600 : /* If stage 1 finished normally, p is the smallest prime >B1 here.
601 : In that case, set to B1 */
602 975 : if (p > B1)
603 915 : p = B1;
604 :
605 975 : if (p > *B1done)
606 975 : *B1done = p;
607 :
608 975 : if (chkfilename != NULL)
609 10 : writechkfile (chkfilename, ECM_ECM, *B1done, n, A, x, NULL, z);
610 :
611 975 : prime_info_clear (prime_info);
612 :
613 975 : if (!mpres_invert (u, z, n)) /* Factor found? */
614 : {
615 222 : mpres_gcd (f, z, n);
616 222 : ret = ECM_FACTOR_FOUND_STEP1;
617 : }
618 975 : mpres_mul (x, x, u, n);
619 :
620 975 : mpres_clear (zT2, n);
621 975 : mpres_clear (xT2, n);
622 975 : mpres_clear (zT, n);
623 975 : mpres_clear (xT, n);
624 975 : mpres_clear (zC, n);
625 975 : mpres_clear (xC, n);
626 975 : mpres_clear (zB, n);
627 975 : mpres_clear (xB, n);
628 975 : mpres_clear (w, n);
629 975 : mpres_clear (v, n);
630 975 : mpres_clear (u, n);
631 975 : mpres_clear (z, n);
632 975 : mpres_clear (b, n);
633 :
634 975 : return ret;
635 : }
636 :
637 : #define DEBUG_EC_W 0
638 :
639 : #ifdef HAVE_ADDLAWS
640 : /* Input: when Etype == ECM_EC_TYPE_WEIERSTRASS*:
641 : (x, y) is initial point
642 : A is curve parameter in Weierstrass's form:
643 : Y^2 = X^3 + A*X + B, where B = y^2-(x^3+A*x) is implicit
644 : when Etype == ECM_EC_TYPE_TWISTED_HESSIAN:
645 : (x, y) is initial point
646 : A=a/d is curve parameter in Hessian form: a*X^3+Y^3+Z^3=d*X*Y*Z
647 : n is the number to factor
648 : B1 is the stage 1 bound
649 : batch_s = prod(p^e <= B1) if != 1
650 : Output: If a factor is found, it is returned in f.
651 : Otherwise, (x, y) contains the point computed in stage 1.
652 : B1done is set to B1 if stage 1 completed normally,
653 : or to the largest prime processed if interrupted, but never
654 : to a smaller value than B1done was upon function entry.
655 : Return value: ECM_FACTOR_FOUND_STEP1 if a factor is found, otherwise
656 : ECM_NO_FACTOR_FOUND
657 : */
658 : static int
659 280 : ecm_stage1_W (mpz_t f, ell_curve_t E, ell_point_t P, mpmod_t n,
660 : double B1, double *B1done, mpz_t batch_s, mpz_t go,
661 : int (*stop_asap)(void), char *chkfilename)
662 : {
663 : mpres_t xB;
664 : ell_point_t Q;
665 280 : uint64_t p = 0, r, last_chkpnt_p;
666 280 : int ret = ECM_NO_FACTOR_FOUND, status;
667 : long last_chkpnt_time;
668 : prime_info_t prime_info;
669 :
670 280 : prime_info_init (prime_info);
671 :
672 280 : mpres_init (xB, n);
673 :
674 280 : ell_point_init(Q, E, n);
675 :
676 280 : last_chkpnt_time = cputime ();
677 :
678 : #if DEBUG_EC_W >= 2
679 : gmp_printf("N:=%Zd;\n", n->orig_modulus);
680 : printf("E:="); ell_curve_print(E, n);
681 : printf("E:=[E[4], E[5]];\n");
682 : printf("P:="); ell_point_print(P, E, n); printf("; Q:=P;\n");
683 : #endif
684 : /* preload group order */
685 280 : if (go != NULL){
686 280 : if (ell_point_mul (f, Q, go, P, E, n) == 0){
687 0 : ret = ECM_FACTOR_FOUND_STEP1;
688 0 : goto end_of_stage1_w;
689 : }
690 280 : ell_point_set(P, Q, E, n);
691 : }
692 : #if DEBUG_EC_W >= 1
693 : gmp_printf("go:=%Zd;\n", go);
694 : printf("goP:="); ell_point_print(P, E, n); printf(";\n");
695 : #endif
696 280 : if(mpz_cmp_ui(batch_s, 1) == 0){
697 280 : outputf (OUTPUT_VERBOSE, "Using traditional approach to Step 1\n");
698 3600 : for (r = 2; r <= B1; r *= 2)
699 3320 : if (r > *B1done){
700 3320 : if(ell_point_duplicate (f, Q, P, E, n) == 0){
701 0 : ret = ECM_FACTOR_FOUND_STEP1;
702 0 : goto end_of_stage1_w;
703 : }
704 3320 : ell_point_set(P, Q, E, n);
705 : #if DEBUG_EC_W >= 2
706 : printf("P%ld:=", (long)r); ell_point_print(P, E, n); printf(";\n");
707 : printf("Q:=EcmMult(2, Q, E, N);\n");
708 : printf("(Q[1]*P%ld[3]-Q[3]*P%ld[1]) mod N;\n",(long)r,(long)r);
709 : #endif
710 : }
711 :
712 280 : last_chkpnt_p = 3;
713 364210 : for (p = getprime_mt (prime_info); p <= B1; p = getprime_mt (prime_info)){
714 739250 : for (r = p; r <= B1; r *= p){
715 375280 : if (r > *B1done){
716 375280 : mpz_set_ui(f, (ecm_uint)p);
717 375280 : status = ell_point_mul (f, Q, f, P, E, n);
718 375280 : if(status == 0){
719 : }
720 375240 : else if(E->law == ECM_LAW_HOMOGENEOUS){
721 257570 : if(E->type == ECM_EC_TYPE_TWISTED_HESSIAN)
722 79430 : mpres_gcd(f, Q->x, n);
723 : else
724 178140 : mpres_gcd(f, Q->z, n);
725 : // gmp_printf("gcd=%Zd\n", f);
726 257570 : if(mpz_cmp(f, n->orig_modulus) < 0
727 257260 : && mpz_cmp_ui(f, 1) > 0)
728 110 : status = 0;
729 : }
730 375280 : if(status == 0){
731 150 : ret = ECM_FACTOR_FOUND_STEP1;
732 150 : goto end_of_stage1_w;
733 : }
734 : #if DEBUG_EC_W >= 2
735 : printf("R%ld:=", (long)r); ell_point_print(Q, E, n);
736 : printf(";\nQ:=EcmMult(%ld, Q, E, N);\n", (long)p);
737 : printf("(Q[1]*R%ld[3]-Q[3]*R%ld[1]) mod N;\n",(long)r,(long)r);
738 : #endif
739 375130 : ell_point_set(P, Q, E, n);
740 : }
741 : }
742 363970 : if (ell_point_is_zero (P, E, n)){
743 40 : outputf (OUTPUT_VERBOSE, "Reached point at infinity, "
744 : "%.0f divides group orders\n", p);
745 40 : break;
746 : }
747 :
748 363930 : if (stop_asap != NULL && (*stop_asap) ()){
749 0 : outputf (OUTPUT_NORMAL, "Interrupted at prime %.0f\n", p);
750 0 : break;
751 : }
752 :
753 363930 : if (chkfilename != NULL && p > last_chkpnt_p + 10000 &&
754 0 : elltime (last_chkpnt_time, cputime ()) > CHKPNT_PERIOD){
755 0 : writechkfile (chkfilename, ECM_ECM, MAX(p, *B1done),
756 0 : n, E->a4, P->x, P->y, P->z);
757 0 : last_chkpnt_p = p;
758 0 : last_chkpnt_time = cputime ();
759 : }
760 : }
761 : }
762 : else{
763 : #if USE_ADD_SUB_CHAINS == 0 /* keeping it simple */
764 0 : if (ell_point_mul (f, Q, batch_s, P, E, n) == 0){
765 0 : ret = ECM_FACTOR_FOUND_STEP1;
766 0 : goto end_of_stage1_w;
767 : }
768 : #else
769 : /* batch mode and special coding... */
770 : short *S = NULL;
771 : size_t iS;
772 : int w;
773 : add_sub_unpack(&w, &S, &iS, batch_s);
774 : if (ell_point_mul_add_sub_with_S(f, Q, P, E, n, w, S, iS) == 0){
775 : ret = ECM_FACTOR_FOUND_STEP1;
776 : }
777 : #endif
778 0 : ell_point_set(P, Q, E, n);
779 0 : p = B1;
780 : }
781 280 : end_of_stage1_w:
782 : /* If stage 1 finished normally, p is the smallest prime > B1 here.
783 : In that case, set to B1 */
784 280 : if (p > B1)
785 90 : p = B1;
786 :
787 280 : if (p > *B1done)
788 280 : *B1done = p;
789 :
790 280 : if (chkfilename != NULL)
791 0 : writechkfile (chkfilename, ECM_ECM, *B1done, n, E->a4, P->x, P->y,P->z);
792 280 : prime_info_clear (prime_info);
793 : #if DEBUG_EC_W >= 2
794 : printf("lastP="); ell_point_print(P, E, n); printf("\n");
795 : #endif
796 280 : if(ret != ECM_FACTOR_FOUND_STEP1){
797 130 : if(ell_point_is_zero(P, E, n) == 1){
798 : /* too bad */
799 40 : ell_point_set_to_zero(P, E, n);
800 40 : mpz_set(f, n->orig_modulus);
801 40 : ret = ECM_FACTOR_FOUND_STEP1;
802 : }
803 : else{
804 : /* for affine case, z = 1 anyway */
805 90 : if(E->law == ECM_LAW_HOMOGENEOUS){
806 60 : if (!mpres_invert (xB, P->z, n)){ /* Factor found? */
807 0 : mpres_gcd (f, P->z, n);
808 0 : gmp_printf("# factor found during normalization: %Zd\n", f);
809 0 : ret = ECM_FACTOR_FOUND_STEP1;
810 : }
811 : else{
812 : /* normalize to get (x:y:1) valid in W or H form... */
813 : #if DEBUG_EC_W >= 2
814 : mpres_get_z(f, xB, n); gmp_printf("1/z=%Zd\n", f);
815 : #endif
816 60 : mpres_mul (P->x, P->x, xB, n);
817 60 : mpres_mul (P->y, P->y, xB, n);
818 : #if DEBUG_EC_W >= 2
819 : mpres_get_z(f, P->x, n); gmp_printf("x/z=%Zd\n", f);
820 : mpres_get_z(f, P->y, n); gmp_printf("y/z=%Zd\n", f);
821 : #endif
822 60 : mpres_set_ui (P->z, 1, n);
823 : }
824 : }
825 : }
826 : }
827 :
828 280 : mpres_clear (xB, n);
829 280 : ell_point_clear(Q, E, n);
830 :
831 280 : return ret;
832 : }
833 : #endif
834 :
835 : /* choose "optimal" S according to step 2 range B2 */
836 : int
837 1524 : choose_S (mpz_t B2len)
838 : {
839 1524 : if (mpz_cmp_d (B2len, 1e7) < 0)
840 1133 : return 1; /* x^1 */
841 391 : else if (mpz_cmp_d (B2len, 1e8) < 0)
842 281 : return 2; /* x^2 */
843 110 : else if (mpz_cmp_d (B2len, 1e9) < 0)
844 53 : return -3; /* Dickson(3) */
845 57 : else if (mpz_cmp_d (B2len, 1e10) < 0)
846 46 : return -6; /* Dickson(6) */
847 11 : else if (mpz_cmp_d (B2len, 3e11) < 0)
848 10 : return -12; /* Dickson(12) */
849 : else
850 1 : return -30; /* Dickson(30) */
851 : }
852 :
853 : #define DIGITS_START 35
854 : #define DIGITS_INCR 5
855 : #define DIGITS_END 80
856 :
857 : void
858 57 : print_expcurves (double B1, const mpz_t B2, unsigned long dF, unsigned long k,
859 : int S, int param)
860 : {
861 : double prob;
862 : int i, j;
863 : char sep, outs[128], flt[16];
864 : double smoothness_correction;
865 :
866 57 : if (param == ECM_PARAM_SUYAMA || param == ECM_PARAM_BATCH_2)
867 40 : smoothness_correction = 1.0;
868 17 : else if (param == ECM_PARAM_BATCH_SQUARE)
869 17 : smoothness_correction = EXTRA_SMOOTHNESS_SQUARE;
870 0 : else if (param == ECM_PARAM_BATCH_32BITS_D)
871 0 : smoothness_correction = EXTRA_SMOOTHNESS_32BITS_D;
872 : else /* This case should never happen */
873 0 : smoothness_correction = 0.0;
874 :
875 627 : for (i = DIGITS_START, j = 0; i <= DIGITS_END; i += DIGITS_INCR)
876 570 : j += sprintf (outs + j, "%u%c", i, (i < DIGITS_END) ? '\t' : '\n');
877 57 : outs[j] = '\0';
878 57 : outputf (OUTPUT_VERBOSE, "Expected number of curves to find a factor "
879 : "of n digits (assuming one exists):\n%s", outs);
880 627 : for (i = DIGITS_START; i <= DIGITS_END; i += DIGITS_INCR)
881 : {
882 570 : sep = (i < DIGITS_END) ? '\t' : '\n';
883 570 : prob = ecmprob (B1, mpz_get_d (B2),
884 : /* smoothness depends on the parametrization */
885 570 : pow (10., i - .5) / smoothness_correction,
886 570 : (double) dF * dF * k, S);
887 570 : if (prob > 1. / 10000000)
888 40 : outputf (OUTPUT_VERBOSE, "%.0f%c", floor (1. / prob + .5), sep);
889 530 : else if (prob > 0.)
890 : {
891 : /* on Windows: 2.6e+009 or 2.6e+025 (8 characters in string) */
892 : /* on Linux : 2.6e+09 or 2.6e+25 (7 characters in string) */
893 : /* This will normalize the output so that the Windows ouptut will match the Linux output */
894 64 : if (sprintf (flt, "%.2g", floor (1. / prob + .5)) == 8)
895 0 : memmove (&flt[5], &flt[6], strlen(flt) - 5);
896 64 : outputf (OUTPUT_VERBOSE, "%s%c", flt, sep);
897 : }
898 : else
899 466 : outputf (OUTPUT_VERBOSE, "Inf%c", sep);
900 : }
901 57 : }
902 :
903 : void
904 9 : print_exptime (double B1, const mpz_t B2, unsigned long dF, unsigned long k,
905 : int S, double tottime, int param)
906 : {
907 : double prob, exptime;
908 : int i, j;
909 : char sep, outs[128];
910 : double smoothness_correction;
911 :
912 9 : if (param == ECM_PARAM_SUYAMA || param == ECM_PARAM_BATCH_2)
913 0 : smoothness_correction = 1.0;
914 9 : else if (param == ECM_PARAM_BATCH_SQUARE)
915 9 : smoothness_correction = EXTRA_SMOOTHNESS_SQUARE;
916 0 : else if (param == ECM_PARAM_BATCH_32BITS_D)
917 0 : smoothness_correction = EXTRA_SMOOTHNESS_32BITS_D;
918 : else /* This case should never happen */
919 0 : smoothness_correction = 0.0;
920 :
921 99 : for (i = DIGITS_START, j = 0; i <= DIGITS_END; i += DIGITS_INCR)
922 90 : j += sprintf (outs + j, "%u%c", i, (i < DIGITS_END) ? '\t' : '\n');
923 9 : outs[j] = '\0';
924 9 : outputf (OUTPUT_VERBOSE, "Expected time to find a factor of n digits:\n%s",
925 : outs);
926 99 : for (i = DIGITS_START; i <= DIGITS_END; i += DIGITS_INCR)
927 : {
928 90 : sep = (i < DIGITS_END) ? '\t' : '\n';
929 90 : prob = ecmprob (B1, mpz_get_d (B2),
930 : /* in batch mode, the extra smoothness is smaller */
931 90 : pow (10., i - .5) / smoothness_correction,
932 90 : (double) dF * dF * k, S);
933 90 : exptime = (prob > 0.) ? tottime / prob : HUGE_VAL;
934 90 : outputf (OUTPUT_TRACE, "Digits: %d, Total time: %.0f, probability: "
935 : "%g, expected time: %.0f\n", i, tottime, prob, exptime);
936 90 : if (exptime < 1000.)
937 0 : outputf (OUTPUT_VERBOSE, "%.0fms%c", exptime, sep);
938 90 : else if (exptime < 60000.) /* One minute */
939 0 : outputf (OUTPUT_VERBOSE, "%.2fs%c", exptime / 1000., sep);
940 90 : else if (exptime < 3600000.) /* One hour */
941 0 : outputf (OUTPUT_VERBOSE, "%.2fm%c", exptime / 60000., sep);
942 90 : else if (exptime < 86400000.) /* One day */
943 9 : outputf (OUTPUT_VERBOSE, "%.2fh%c", exptime / 3600000., sep);
944 81 : else if (exptime < 31536000000.) /* One year */
945 23 : outputf (OUTPUT_VERBOSE, "%.2fd%c", exptime / 86400000., sep);
946 58 : else if (exptime < 31536000000000.) /* One thousand years */
947 16 : outputf (OUTPUT_VERBOSE, "%.2fy%c", exptime / 31536000000., sep);
948 42 : else if (exptime < 31536000000000000.) /* One million years */
949 16 : outputf (OUTPUT_VERBOSE, "%.0fy%c", exptime / 31536000000., sep);
950 26 : else if (prob > 0.)
951 16 : outputf (OUTPUT_VERBOSE, "%.1gy%c", exptime / 31536000000., sep);
952 : else
953 10 : outputf (OUTPUT_VERBOSE, "Inf%c", sep);
954 : }
955 9 : }
956 :
957 : /* y should be NULL for P+1, and P-1, it contains the y coordinate for the
958 : Weierstrass form for ECM (when sigma_is_A = -1). */
959 : void
960 2021 : print_B1_B2_poly (int verbosity, int method, double B1, double B1done,
961 : mpz_t B2min_param, mpz_t B2min, mpz_t B2, int S, mpz_t sigma,
962 : int sigma_is_A, int Etype,
963 : mpz_t y, int param, unsigned int nb_curves)
964 : {
965 : ASSERT ((method == ECM_ECM) || (y == NULL));
966 : ASSERT ((-1 <= sigma_is_A) && (sigma_is_A <= 1));
967 : ASSERT (param != ECM_PARAM_DEFAULT || sigma_is_A == 1 || sigma_is_A == -1);
968 :
969 2021 : if (test_verbose (verbosity))
970 : {
971 2001 : outputf (verbosity, "Using ");
972 2001 : if (ECM_IS_DEFAULT_B1_DONE(B1done))
973 1906 : outputf (verbosity, "B1=%1.0f, ", B1);
974 : else
975 95 : outputf (verbosity, "B1=%1.0f-%1.0f, ", B1done, B1);
976 2001 : if (mpz_sgn (B2min_param) < 0)
977 1825 : outputf (verbosity, "B2=%Zd", B2);
978 : else
979 176 : outputf (verbosity, "B2=%Zd-%Zd", B2min, B2);
980 :
981 2001 : if (S > 0)
982 1891 : outputf (verbosity, ", polynomial x^%u", S);
983 110 : else if (S < 0)
984 110 : outputf (verbosity, ", polynomial Dickson(%u)", -S);
985 :
986 : /* don't print in resume case, since x0 is saved in resume file */
987 2001 : if (method == ECM_ECM)
988 : {
989 1470 : if (sigma_is_A == 1)
990 86 : outputf (verbosity, ", A=%Zd", sigma);
991 1384 : else if (sigma_is_A == 0)
992 : {
993 1074 : if (nb_curves > 1)
994 : {
995 0 : outputf (verbosity, ", sigma=%d:%Zd", param, sigma);
996 0 : mpz_add_ui (sigma, sigma, nb_curves-1);
997 0 : outputf (verbosity, "-%d:%Zd", param, sigma);
998 0 : mpz_sub_ui (sigma, sigma, nb_curves-1);
999 0 : outputf (verbosity, " (%u curves)", nb_curves);
1000 : }
1001 : else
1002 1074 : outputf (verbosity, ", sigma=%d:%Zd", param, sigma);
1003 : }
1004 : else{
1005 310 : if (Etype == ECM_EC_TYPE_WEIERSTRASS)
1006 190 : outputf (verbosity, ", Weierstrass(A=%Zd,y=%Zd)", sigma, y);
1007 120 : else if (Etype == ECM_EC_TYPE_HESSIAN)
1008 70 : outputf (verbosity, ", Hessian(D=%Zd,y=%Zd)", sigma, y);
1009 50 : else if (Etype == ECM_EC_TYPE_TWISTED_HESSIAN)
1010 30 : outputf (verbosity, ", twisted Hessian(y=%Zd)", y);
1011 : }
1012 : }
1013 531 : else if (ECM_IS_DEFAULT_B1_DONE(B1done))
1014 : /* in case of P-1 or P+1, we store the initial point in sigma */
1015 486 : outputf (verbosity, ", x0=%Zd", sigma);
1016 :
1017 2001 : outputf (verbosity, "\n");
1018 : }
1019 2021 : }
1020 :
1021 : /* Compute parameters for stage 2 */
1022 : int
1023 1566 : set_stage_2_params (mpz_t B2, mpz_t B2_parm, mpz_t B2min, mpz_t B2min_parm,
1024 : root_params_t *root_params, double B1,
1025 : unsigned long *k, const int S, int use_ntt, int *po2,
1026 : unsigned long *dF, char *TreeFilename, double maxmem,
1027 : int Fermat, mpmod_t modulus)
1028 : {
1029 1566 : mpz_set (B2min, B2min_parm);
1030 1566 : mpz_set (B2, B2_parm);
1031 :
1032 1566 : mpz_init (root_params->i0);
1033 :
1034 : /* set second stage bound B2: when using polynomial multiplication of
1035 : complexity n^alpha, stage 2 has complexity about B2^(alpha/2), and
1036 : we want stage 2 to take about half of stage 1, thus we choose
1037 : B2 = (c*B1)^(2/alpha). Experimentally, c=1/4 seems to work well.
1038 : For Toom-Cook 3, this gives alpha=log(5)/log(3), and B2 ~ (c*B1)^1.365.
1039 : For Toom-Cook 4, this gives alpha=log(7)/log(4), and B2 ~ (c*B1)^1.424. */
1040 :
1041 : /* We take the cost of P+1 stage 1 to be about twice that of P-1.
1042 : Since nai"ve P+1 and ECM cost respectively 2 and 11 multiplies per
1043 : addition and duplicate, and both are optimized with PRAC, we can
1044 : assume the ratio remains about 11/2. */
1045 :
1046 : /* Also scale B2 by what the user said (or by the default scaling of 1.0) */
1047 :
1048 1566 : if (ECM_IS_DEFAULT_B2(B2))
1049 939 : mpz_set_d (B2, pow (ECM_COST * B1, DEFAULT_B2_EXPONENT));
1050 :
1051 : /* set B2min */
1052 1566 : if (mpz_sgn (B2min) < 0)
1053 1430 : mpz_set_d (B2min, B1);
1054 :
1055 : /* Let bestD determine parameters for root generation and the
1056 : effective B2 */
1057 :
1058 1566 : if (use_ntt)
1059 1393 : *po2 = 1;
1060 :
1061 1566 : root_params->d2 = 0; /* Enable automatic choice of d2 */
1062 1566 : if (bestD (root_params, k, dF, B2min, B2, *po2, use_ntt, maxmem,
1063 : (TreeFilename != NULL), modulus) == ECM_ERROR)
1064 20 : return ECM_ERROR;
1065 :
1066 : /* Set default degree for Brent-Suyama extension */
1067 : /* We try to keep the time used by the Brent-Suyama extension
1068 : at about 10% of the stage 2 time */
1069 : /* Degree S Dickson polys and x^S are equally fast for ECM, so we go for
1070 : the better Dickson polys whenever possible. For S == 1, 2, they behave
1071 : identically. */
1072 :
1073 1546 : root_params->S = S;
1074 1546 : if (root_params->S == ECM_DEFAULT_S)
1075 : {
1076 1536 : if (Fermat > 0)
1077 : {
1078 : /* For Fermat numbers, default is 1 (no Brent-Suyama) */
1079 12 : root_params->S = 1;
1080 : }
1081 : else
1082 : {
1083 : mpz_t t;
1084 1524 : mpz_init (t);
1085 1524 : mpz_sub (t, B2, B2min);
1086 1524 : root_params->S = choose_S (t);
1087 1524 : mpz_clear (t);
1088 : }
1089 : }
1090 1546 : return ECM_NO_FACTOR_FOUND;
1091 : }
1092 :
1093 : /* Input: x is starting point or zero
1094 : y is used for Weierstrass curves (even in Step1)
1095 : sigma is sigma value (if x is set to zero) or
1096 : A parameter (if x is non-zero) of curve
1097 : n is the number to factor
1098 : go is the initial group order to preload
1099 : B1, B2 are the stage 1/stage 2 bounds, respectively
1100 : B2min the lower bound for stage 2
1101 : k is the number of blocks to do in stage 2
1102 : S is the degree of the Suyama-Brent extension for stage 2
1103 : verbose is verbosity level: 0 no output, 1 normal output,
1104 : 2 diagnostic output.
1105 : sigma_is_A: If true, the sigma parameter contains the curve's A value
1106 : Etype
1107 : zE is a curve that is used when a special torsion group was used; in
1108 : that case, (x, y) must be a point on E.
1109 : Output: f is the factor found.
1110 : Return value: ECM_FACTOR_FOUND_STEPn if a factor was found,
1111 : ECM_NO_FACTOR_FOUND if no factor was found,
1112 : ECM_ERROR in case of error.
1113 : (x, y) contains the new point at the end of Stage 1.
1114 : */
1115 : int
1116 1588 : ecm (mpz_t f, mpz_t x, mpz_t y, int param, mpz_t sigma, mpz_t n, mpz_t go,
1117 : double *B1done, double B1, mpz_t B2min_parm, mpz_t B2_parm,
1118 : unsigned long k, const int S, int verbose, int repr, int nobase2step2,
1119 : int use_ntt, int sigma_is_A, ell_curve_t zE,
1120 : FILE *os, FILE* es, char *chkfilename, char
1121 : *TreeFilename, double maxmem, double stage1time, gmp_randstate_t rng, int
1122 : (*stop_asap)(void), mpz_t batch_s, double *batch_last_B1_used,
1123 : ATTRIBUTE_UNUSED double gw_k, ATTRIBUTE_UNUSED unsigned long gw_b,
1124 : ATTRIBUTE_UNUSED unsigned long gw_n, ATTRIBUTE_UNUSED signed long gw_c)
1125 : {
1126 1588 : int youpi = ECM_NO_FACTOR_FOUND;
1127 1588 : int base2 = 0; /* If n is of form 2^n[+-]1, set base to [+-]n */
1128 1588 : int Fermat = 0; /* If base2 > 0 is a power of 2, set Fermat to base2 */
1129 1588 : int po2 = 0; /* Whether we should use power-of-2 poly degree */
1130 : long st;
1131 : mpmod_t modulus;
1132 : curve P;
1133 : ell_curve_t E;
1134 : #ifdef HAVE_ADDLAWS
1135 : ell_point_t PE;
1136 : #endif
1137 : mpz_t B2min, B2; /* Local B2, B2min to avoid changing caller's values */
1138 : unsigned long dF;
1139 : root_params_t root_params;
1140 :
1141 : /* 1: sigma contains A from Montgomery form By^2 = x^3 + Ax^2 + x
1142 : 0: sigma contains sigma
1143 : -1: sigma contains A from Weierstrass form y^2 = x^3 + Ax + B,
1144 : and y contains y */
1145 : ASSERT((-1 <= sigma_is_A) && (sigma_is_A <= 1));
1146 : ASSERT((GMP_NUMB_BITS == 32) || (GMP_NUMB_BITS == 64));
1147 :
1148 1588 : set_verbose (verbose);
1149 1588 : ECM_STDOUT = (os == NULL) ? stdout : os;
1150 1588 : ECM_STDERR = (es == NULL) ? stdout : es;
1151 :
1152 : #ifdef MPRESN_NO_ADJUSTMENT
1153 : /* When no adjustment is made in mpresn_ functions, N should be smaller
1154 : than B^n/16 */
1155 : if (mpz_sizeinbase (n, 2) > mpz_size (n) * GMP_NUMB_BITS - 4)
1156 : {
1157 : outputf (OUTPUT_ERROR, "Error, N should be smaller than B^n/16\n");
1158 : return ECM_ERROR;
1159 : }
1160 : #endif
1161 :
1162 : /* if a batch mode is requested by the user, this implies ECM_MOD_MODMULN */
1163 1588 : if (repr == ECM_MOD_DEFAULT && IS_BATCH_MODE(param))
1164 154 : repr = ECM_MOD_MODMULN;
1165 :
1166 : /* choose the arithmetic used before the parametrization, since for divisors
1167 : of 2^n+/-1 the default choice param=1 might not be optimal */
1168 1588 : if (mpmod_init (modulus, n, repr) != 0)
1169 10 : return ECM_ERROR;
1170 :
1171 1578 : repr = modulus->repr;
1172 :
1173 : /* If the parametrization is not given, choose it. */
1174 1578 : if (param == ECM_PARAM_DEFAULT)
1175 41 : param = get_default_param (sigma_is_A, *B1done, repr);
1176 : /* when dealing with several input numbers, if we had already computed
1177 : batch_s, but the new number uses the base-2 representation, then we
1178 : are forced to use ECM_PARAM_SUYAMA, and we reset batch_s to 1 to avoid
1179 : the error "-bsaves/-bloads makes sense in batch mode only" below */
1180 1578 : if (param == ECM_PARAM_SUYAMA)
1181 1065 : mpz_set_ui (batch_s, 1);
1182 :
1183 : /* In batch mode,
1184 : we force repr=MODMULN,
1185 : B1done should be either the default value or greater than B1
1186 : x should be either 0 (undetermined) or 2 */
1187 1578 : if (IS_BATCH_MODE(param))
1188 : {
1189 183 : if (repr != ECM_MOD_MODMULN)
1190 : {
1191 4 : outputf (OUTPUT_ERROR, "Error, with param %d, repr should be "
1192 : "ECM_MOD_MODMULN.\n", param);
1193 4 : return ECM_ERROR;
1194 : }
1195 :
1196 179 : if (!ECM_IS_DEFAULT_B1_DONE(*B1done) && *B1done < B1)
1197 : {
1198 0 : outputf (OUTPUT_ERROR, "Error, cannot resume with param %d, except "
1199 : "for doing only stage 2\n", param);
1200 0 : return ECM_ERROR;
1201 : }
1202 :
1203 179 : if (sigma_is_A >= 0 && mpz_sgn (x) != 0 && mpz_cmp_ui (x, 2) != 0)
1204 : {
1205 8 : if (ECM_IS_DEFAULT_B1_DONE(*B1done))
1206 : {
1207 8 : outputf (OUTPUT_ERROR, "Error, x0 should be equal to 2 with this "
1208 : "parametrization\n");
1209 8 : return ECM_ERROR;
1210 : }
1211 : }
1212 : }
1213 :
1214 : /* check that if ECM_PARAM_BATCH_SQUARE is used, GMP_NUMB_BITS == 64 */
1215 : if (param == ECM_PARAM_BATCH_SQUARE && GMP_NUMB_BITS == 32)
1216 : {
1217 : outputf (OUTPUT_ERROR, "Error, parametrization ECM_PARAM_BATCH_SQUARE "
1218 : "works only with GMP_NUMB_BITS=64\n");
1219 : return ECM_ERROR;
1220 : }
1221 :
1222 : /* check that B1 is not too large */
1223 1566 : if (B1 > (double) ECM_UINT_MAX)
1224 : {
1225 0 : outputf (OUTPUT_ERROR, "Error, maximal step 1 bound for ECM is %lu.\n",
1226 : ECM_UINT_MAX);
1227 0 : return ECM_ERROR;
1228 : }
1229 :
1230 : /* See what kind of number we have as that may influence optimal parameter
1231 : selection. Test for base 2 number. Note: this was already done by
1232 : mpmod_init. */
1233 :
1234 1566 : if (modulus->repr == ECM_MOD_BASE2)
1235 161 : base2 = modulus->bits;
1236 :
1237 : /* For a Fermat number (base2 a positive power of 2) */
1238 1663 : for (Fermat = base2; Fermat > 0 && (Fermat & 1) == 0; Fermat >>= 1);
1239 1566 : if (Fermat == 1)
1240 : {
1241 12 : Fermat = base2;
1242 12 : po2 = 1;
1243 : }
1244 : else
1245 1554 : Fermat = 0;
1246 :
1247 1566 : mpres_init (P.x, modulus);
1248 1566 : mpres_init (P.y, modulus);
1249 1566 : mpres_init (P.A, modulus);
1250 :
1251 : #ifdef HAVE_ADDLAWS
1252 1566 : ell_curve_set_z (E, zE, modulus);
1253 : #else
1254 : E->type = ECM_EC_TYPE_MONTGOMERY;
1255 : #endif
1256 :
1257 1566 : mpz_init (B2);
1258 1566 : mpz_init (B2min);
1259 1566 : youpi = set_stage_2_params (B2, B2_parm, B2min, B2min_parm,
1260 : &root_params, B1, &k, S, use_ntt,
1261 : &po2, &dF, TreeFilename, maxmem, Fermat,modulus);
1262 :
1263 : /* if the user gave B2, print that B2 on the Using B1=..., B2=... line */
1264 1566 : if(!ECM_IS_DEFAULT_B2(B2_parm))
1265 627 : mpz_set (B2, B2_parm);
1266 :
1267 1566 : if (youpi == ECM_ERROR)
1268 20 : goto end_of_ecm;
1269 :
1270 1546 : if (sigma_is_A == 0)
1271 : {
1272 1140 : if (mpz_sgn (sigma) == 0)
1273 : {
1274 29 : youpi = get_curve_from_random_parameter (f, P.A, P.x, sigma, param,
1275 : modulus, rng);
1276 :
1277 29 : if (youpi == ECM_ERROR)
1278 : {
1279 10 : outputf (OUTPUT_ERROR, "Error, invalid parametrization.\n");
1280 10 : goto end_of_ecm;
1281 : }
1282 : }
1283 : else /* Compute A and x0 from given sigma values */
1284 : {
1285 1111 : if (param == ECM_PARAM_SUYAMA)
1286 1007 : youpi = get_curve_from_param0 (f, P.A, P.x, sigma, modulus);
1287 104 : else if (param == ECM_PARAM_BATCH_SQUARE)
1288 16 : youpi = get_curve_from_param1 (P.A, P.x, sigma, modulus);
1289 88 : else if (param == ECM_PARAM_BATCH_2)
1290 56 : youpi = get_curve_from_param2 (f, P.A, P.x, sigma, modulus);
1291 32 : else if (param == ECM_PARAM_BATCH_32BITS_D)
1292 32 : youpi = get_curve_from_param3 (P.A, P.x, sigma, modulus);
1293 : else
1294 : {
1295 0 : outputf (OUTPUT_ERROR, "Error, invalid parametrization.\n");
1296 0 : youpi = ECM_ERROR;
1297 0 : goto end_of_ecm;
1298 : }
1299 :
1300 1111 : if (youpi != ECM_NO_FACTOR_FOUND)
1301 : {
1302 36 : if (youpi == ECM_ERROR)
1303 0 : outputf (OUTPUT_ERROR, "Error, invalid value of sigma.\n");
1304 :
1305 36 : goto end_of_ecm;
1306 : }
1307 : }
1308 : }
1309 406 : else if (sigma_is_A == 1)
1310 : {
1311 : /* sigma contains the A value */
1312 96 : mpres_set_z (P.A, sigma, modulus);
1313 : /* TODO: make a valid, random starting point in case none was given */
1314 : /* Problem: this may be as hard as factoring as we'd need to determine
1315 : whether x^3 + a*x^2 + x is a quadratic residue or not */
1316 : /* For now, we'll just chicken out. */
1317 : /* Except for batch mode where we know that x0=2 */
1318 96 : if (mpz_sgn (x) == 0)
1319 : {
1320 66 : if (IS_BATCH_MODE(param))
1321 56 : mpres_set_ui (P.x, 2, modulus);
1322 : else
1323 : {
1324 10 : outputf (OUTPUT_ERROR,
1325 : "Error, -A requires a starting point (-x0 x).\n");
1326 10 : youpi = ECM_ERROR;
1327 10 : goto end_of_ecm;
1328 : }
1329 : }
1330 : else
1331 30 : mpres_set_z (P.x, x, modulus);
1332 : }
1333 :
1334 : /* Print B1, B2, polynomial and sigma */
1335 1490 : print_B1_B2_poly (OUTPUT_NORMAL, ECM_ECM, B1, *B1done, B2min_parm, B2min,
1336 : B2, root_params.S, sigma, sigma_is_A, E->type,
1337 : y, param, 0);
1338 :
1339 : #if 0
1340 : outputf (OUTPUT_VERBOSE, "b2=%1.0f, dF=%lu, k=%lu, d=%lu, d2=%lu, i0=%Zd\n",
1341 : b2, dF, k, root_params.d1, root_params.d2, root_params.i0);
1342 : #else
1343 1490 : outputf (OUTPUT_VERBOSE, "dF=%lu, k=%lu, d=%lu, d2=%lu, i0=%Zd\n",
1344 : dF, k, root_params.d1, root_params.d2, root_params.i0);
1345 : #endif
1346 :
1347 1490 : if (sigma_is_A == -1) /* Weierstrass or Hessian form.
1348 : It is known that all curves in Weierstrass form do
1349 : not admit a Montgomery form. However, we could
1350 : be interested in performing some plain Step 1
1351 : on some special curves.
1352 : */
1353 : {
1354 310 : if (param != ECM_PARAM_TORSION)
1355 : {
1356 180 : mpres_set_z (P.A, sigma, modulus); /* sigma contains A */
1357 180 : mpres_set_z (P.x, x, modulus);
1358 180 : mpres_set_z (P.y, y, modulus);
1359 : }
1360 : else
1361 : {
1362 : #ifdef HAVE_ADDLAWS
1363 130 : if(E->type == ECM_EC_TYPE_WEIERSTRASS)
1364 : #endif
1365 70 : mpres_set_z(P.A, zE->a4, modulus);
1366 : #ifdef HAVE_ADDLAWS
1367 60 : else if(E->type == ECM_EC_TYPE_MONTGOMERY)
1368 20 : mpres_set_z(P.A, zE->a2, modulus);
1369 : #endif
1370 130 : mpres_set_z (P.x, x, modulus);
1371 130 : mpres_set_z (P.y, y, modulus);
1372 : }
1373 : }
1374 :
1375 1490 : if (test_verbose (OUTPUT_RESVERBOSE))
1376 : {
1377 : mpz_t t;
1378 :
1379 30 : mpz_init (t);
1380 30 : mpres_get_z (t, P.A, modulus);
1381 30 : outputf (OUTPUT_RESVERBOSE, "A=%Zd\n", t);
1382 30 : mpres_get_z (t, P.x, modulus);
1383 30 : outputf (OUTPUT_RESVERBOSE, "starting point: x0=%Zd\n", t);
1384 : #ifdef HAVE_ADDLAWS
1385 30 : if (E->type == ECM_EC_TYPE_WEIERSTRASS)
1386 : {
1387 0 : mpres_get_z (t, P.y, modulus);
1388 0 : outputf (OUTPUT_RESVERBOSE, " y0=%Zd\n", t);
1389 : }
1390 : #endif
1391 30 : mpz_clear (t);
1392 : }
1393 :
1394 1490 : if (go != NULL && mpz_cmp_ui (go, 1) > 0)
1395 135 : outputf (OUTPUT_VERBOSE, "initial group order: %Zd\n", go);
1396 :
1397 1490 : if (test_verbose (OUTPUT_VERBOSE))
1398 : {
1399 57 : if (mpz_cmp_d (B2min, B1) != 0)
1400 : {
1401 0 : outputf (OUTPUT_VERBOSE,
1402 : "Can't compute success probabilities for B1 <> B2min\n");
1403 : }
1404 57 : else if (param == ECM_PARAM_DEFAULT)
1405 : {
1406 0 : outputf (OUTPUT_VERBOSE, "Can't compute success probabilities "
1407 : "for this parametrization.\n");
1408 : }
1409 : else
1410 : {
1411 57 : rhoinit (256, 10);
1412 57 : print_expcurves (B1, B2, dF, k, root_params.S, param);
1413 : }
1414 : }
1415 :
1416 : /* Compute s for the batch mode */
1417 1490 : if (IS_BATCH_MODE(param) && ECM_IS_DEFAULT_B1_DONE(*B1done) &&
1418 153 : (B1 != *batch_last_B1_used || mpz_cmp_ui (batch_s, 1) <= 0))
1419 : {
1420 137 : *batch_last_B1_used = B1;
1421 :
1422 137 : st = cputime ();
1423 : /* construct the batch exponent */
1424 137 : compute_s (batch_s, B1, NULL);
1425 137 : outputf (OUTPUT_VERBOSE, "Computing batch product (of %" PRIu64
1426 : " bits) of primes up to B1=%1.0f took %ldms\n",
1427 : mpz_sizeinbase (batch_s, 2), B1,
1428 : elltime (st, cputime ()));
1429 : }
1430 :
1431 1490 : st = cputime ();
1432 :
1433 : #ifdef HAVE_GWNUM
1434 : /* We will only use GWNUM for numbers of the form k*b^n+c */
1435 :
1436 : if (gw_b != 0 && B1 >= *B1done && param == ECM_PARAM_SUYAMA)
1437 : youpi = gw_ecm_stage1 (f, &P, modulus, B1, B1done, go, gw_k, gw_b, gw_n, gw_c);
1438 :
1439 : /* At this point B1 == *B1done unless interrupted, or no GWNUM ecm_stage1
1440 : is available */
1441 :
1442 : if (youpi != ECM_NO_FACTOR_FOUND)
1443 : goto end_of_ecm_rhotable;
1444 : #endif /* HAVE_GWNUM */
1445 :
1446 1490 : if (B1 > *B1done || mpz_cmp_ui (go, 1) > 0)
1447 : {
1448 1407 : if (IS_BATCH_MODE(param))
1449 : /* FIXME: go, stop_asap and chkfilename are ignored in batch mode */
1450 152 : youpi = ecm_stage1_batch (f, P.x, P.A, modulus, B1, B1done,
1451 : param, batch_s);
1452 : else{
1453 : #ifdef HAVE_ADDLAWS
1454 1255 : if(E->type == ECM_EC_TYPE_MONTGOMERY)
1455 : #endif
1456 975 : youpi = ecm_stage1 (f, P.x, P.A, modulus, B1, B1done, go,
1457 : stop_asap, chkfilename);
1458 : #ifdef HAVE_ADDLAWS
1459 : else{
1460 280 : ell_point_init(PE, E, modulus);
1461 280 : mpres_set(PE->x, P.x, modulus);
1462 280 : mpres_set(PE->y, P.y, modulus);
1463 280 : youpi = ecm_stage1_W (f, E, PE, modulus, B1, B1done, batch_s,
1464 : go, stop_asap, chkfilename);
1465 280 : mpres_set(P.x, PE->x, modulus);
1466 280 : mpres_set(P.y, PE->y, modulus);
1467 : }
1468 : #endif
1469 : }
1470 : }
1471 83 : else if (mpz_sgn (x) != 0)
1472 : {
1473 : /* when x <> 0, we initialize P to (x:y) */
1474 62 : mpres_set_z (P.x, x, modulus);
1475 62 : mpres_set_z (P.y, y, modulus);
1476 : }
1477 :
1478 1490 : if (stage1time > 0.)
1479 : {
1480 10 : const long st2 = elltime (st, cputime ());
1481 10 : const long s1t = (long) (stage1time * 1000.);
1482 10 : outputf (OUTPUT_NORMAL,
1483 : "Step 1 took %ldms (%ld in this run, %ld from previous runs)\n",
1484 : st2 + s1t, st2, s1t);
1485 : }
1486 : else
1487 1480 : outputf (OUTPUT_NORMAL, "Step 1 took %ldms\n", elltime (st, cputime ()));
1488 :
1489 : /* Store end-of-stage-1 residue in x in case we write it to a save file,
1490 : before P.x is (perhaps) converted to Weierstrass form */
1491 :
1492 1490 : mpres_get_z (x, P.x, modulus);
1493 : #ifdef HAVE_ADDLAWS
1494 1490 : if (E->type == ECM_EC_TYPE_WEIERSTRASS
1495 1300 : || E->type == ECM_EC_TYPE_HESSIAN
1496 1230 : || E->type == ECM_EC_TYPE_TWISTED_HESSIAN)
1497 290 : mpres_get_z (y, P.y, modulus);
1498 : #endif
1499 :
1500 1490 : if (youpi != ECM_NO_FACTOR_FOUND)
1501 412 : goto end_of_ecm_rhotable;
1502 :
1503 1078 : if (test_verbose (OUTPUT_RESVERBOSE))
1504 : {
1505 : mpz_t t;
1506 :
1507 30 : mpz_init (t);
1508 30 : mpres_get_z (t, P.x, modulus);
1509 30 : outputf (OUTPUT_RESVERBOSE, "x=%Zd\n", t);
1510 : #ifdef HAVE_ADDLAWS
1511 30 : if (E->type == ECM_EC_TYPE_WEIERSTRASS)
1512 : {
1513 0 : mpres_get_z (t, P.y, modulus);
1514 0 : outputf (OUTPUT_RESVERBOSE, "y=%Zd\n", t);
1515 : }
1516 : #endif
1517 30 : mpz_clear (t);
1518 : }
1519 :
1520 : /* In case of a signal, we'll exit after the residue is printed. If no save
1521 : file is specified, the user may still resume from the residue */
1522 1078 : if (stop_asap != NULL && (*stop_asap) ())
1523 0 : goto end_of_ecm_rhotable;
1524 :
1525 : /* If using 2^k +/-1 modulus and 'nobase2step2' flag is set,
1526 : set default (-nobase2) modular method and remap P.x, P.y, and P.A */
1527 1078 : if (modulus->repr == ECM_MOD_BASE2 && nobase2step2)
1528 : {
1529 : mpz_t x_t, y_t, A_t;
1530 :
1531 7 : mpz_init (x_t);
1532 7 : mpz_init (y_t);
1533 7 : mpz_init (A_t);
1534 :
1535 7 : mpz_mod (x_t, P.x, modulus->orig_modulus);
1536 7 : mpz_mod (y_t, P.y, modulus->orig_modulus);
1537 7 : mpz_mod (A_t, P.A, modulus->orig_modulus);
1538 :
1539 7 : mpmod_clear (modulus);
1540 :
1541 7 : repr = ECM_MOD_NOBASE2;
1542 7 : if (mpmod_init (modulus, n, repr) != 0) /* reset modulus for nobase2 */
1543 0 : return ECM_ERROR;
1544 :
1545 : /* remap x, y, and A for new modular method */
1546 7 : mpres_set_z (P.x, x_t, modulus);
1547 7 : mpres_set_z (P.y, y_t, modulus);
1548 7 : mpres_set_z (P.A, A_t, modulus);
1549 :
1550 7 : mpz_clear (x_t);
1551 7 : mpz_clear (y_t);
1552 7 : mpz_clear (A_t);
1553 : }
1554 :
1555 : #ifdef HAVE_ADDLAWS
1556 1078 : if (E->type == ECM_EC_TYPE_MONTGOMERY)
1557 : #endif
1558 978 : youpi = montgomery_to_weierstrass (f, P.x, P.y, P.A, modulus);
1559 : #ifdef HAVE_ADDLAWS
1560 100 : else if (E->type == ECM_EC_TYPE_HESSIAN)
1561 : {
1562 20 : youpi = hessian_to_weierstrass (f, P.x, P.y, P.A, modulus);
1563 20 : if(youpi == ECM_NO_FACTOR_FOUND)
1564 : /* due to that non-trivial kernel(?) */
1565 20 : youpi = mult_by_3(f, P.x, P.y, P.A, modulus);
1566 : }
1567 80 : else if (E->type == ECM_EC_TYPE_TWISTED_HESSIAN)
1568 : {
1569 : mpz_t c, rm;
1570 10 : mpz_init(c);
1571 10 : mpz_init(rm);
1572 10 : mpres_get_z(rm, E->a4, modulus);
1573 10 : mpz_rootrem(c, rm, rm, 3);
1574 10 : if(mpz_sgn(rm) != 0){
1575 0 : printf("ECM_EC_TYPE_TWISTED_HESSIAN: not a cube!\n");
1576 0 : exit(-1);
1577 : }
1578 10 : mpres_set_z(P.A, c, modulus);
1579 10 : mpz_clear(c);
1580 10 : mpz_clear(rm);
1581 10 : youpi = twisted_hessian_to_weierstrass (f, P.x, P.y, P.A, E->a6, modulus);
1582 10 : if(youpi == ECM_NO_FACTOR_FOUND){
1583 : /* due to that non-trivial kernel(?) */
1584 10 : youpi = mult_by_3(f, P.x, P.y, P.A, modulus);
1585 : }
1586 : }
1587 : #endif
1588 :
1589 1078 : if (test_verbose (OUTPUT_RESVERBOSE) && youpi == ECM_NO_FACTOR_FOUND &&
1590 30 : mpz_cmp (B2, B2min) >= 0)
1591 : {
1592 : mpz_t t;
1593 :
1594 30 : mpz_init (t);
1595 30 : mpres_get_z (t, P.x, modulus);
1596 30 : outputf (OUTPUT_RESVERBOSE, "After switch to Weierstrass form, "
1597 : "P=(%Zd", t);
1598 30 : mpres_get_z (t, P.y, modulus);
1599 30 : outputf (OUTPUT_RESVERBOSE, ", %Zd)\n", t);
1600 30 : mpres_get_z (t, P.A, modulus);
1601 30 : outputf (OUTPUT_RESVERBOSE, "on curve Y^2 = X^3 + %Zd * X + b\n", t);
1602 30 : mpz_clear (t);
1603 : }
1604 :
1605 1078 : P.disc = 0; /* FIXME: should disappear one day */
1606 :
1607 1078 : if (youpi == ECM_NO_FACTOR_FOUND && mpz_cmp (B2, B2min) >= 0)
1608 1048 : youpi = stage2 (f, &P, modulus, dF, k, &root_params, use_ntt,
1609 : TreeFilename, 0, stop_asap);
1610 : #ifdef TIMING_CRT
1611 : printf ("mpzspv_from_mpzv_slow: %dms\n", mpzspv_from_mpzv_slow_time);
1612 : printf ("mpzspv_to_mpzv: %dms\n", mpzspv_to_mpzv_time);
1613 : printf ("mpzspv_normalise: %dms\n", mpzspv_normalise_time);
1614 : #endif
1615 :
1616 30 : end_of_ecm_rhotable:
1617 1490 : if (test_verbose (OUTPUT_VERBOSE))
1618 : {
1619 57 : if (mpz_cmp_d (B2min, B1) == 0 && param != ECM_PARAM_DEFAULT)
1620 : {
1621 57 : if (youpi == ECM_NO_FACTOR_FOUND &&
1622 0 : (stop_asap == NULL || !(*stop_asap)()))
1623 9 : print_exptime (B1, B2, dF, k, root_params.S,
1624 9 : (long) (stage1time * 1000.) +
1625 9 : elltime (st, cputime ()), param);
1626 57 : rhoinit (1, 0); /* Free memory of rhotable */
1627 : }
1628 : }
1629 :
1630 1433 : end_of_ecm:
1631 : #ifdef HAVE_ADDLAWS
1632 1566 : ell_curve_clear(E, modulus);
1633 : #endif
1634 1566 : mpres_clear (P.A, modulus);
1635 1566 : mpres_clear (P.y, modulus);
1636 1566 : mpres_clear (P.x, modulus);
1637 1566 : mpz_clear (root_params.i0);
1638 1566 : mpz_clear (B2);
1639 1566 : mpz_clear (B2min);
1640 1566 : mpmod_clear (modulus);
1641 1566 : return youpi;
1642 : }
|
