1
2 #include <stdio.h>
3 #include <stdlib.h>
4 #include <string.h>
5 #include <math.h>
6 #include <gmp.h>
7
8 #include "ptypes.h"
9 #include "ecm.h"
10 #include "utility.h"
11 #include "prime_iterator.h"
12
13 #define USE_PRAC
14
15 #define TEST_FOR_2357(n, f) \
16 { \
17 if (mpz_divisible_ui_p(n, 2)) { mpz_set_ui(f, 2); return 1; } \
18 if (mpz_divisible_ui_p(n, 3)) { mpz_set_ui(f, 3); return 1; } \
19 if (mpz_divisible_ui_p(n, 5)) { mpz_set_ui(f, 5); return 1; } \
20 if (mpz_divisible_ui_p(n, 7)) { mpz_set_ui(f, 7); return 1; } \
21 if (mpz_cmp_ui(n, 121) < 0) { return 0; } \
22 }
23
24 /* P3 = P1 + P2 */
_ec_add_AB(mpz_t n,struct ec_affine_point P1,struct ec_affine_point P2,struct ec_affine_point * P3,mpz_t m,mpz_t t1,mpz_t t2)25 static void _ec_add_AB(mpz_t n,
26 struct ec_affine_point P1,
27 struct ec_affine_point P2,
28 struct ec_affine_point *P3,
29 mpz_t m,
30 mpz_t t1,
31 mpz_t t2)
32 {
33 if (!mpz_cmp(P1.x, P2.x)) {
34 mpz_add(t2, P1.y, P2.y);
35 mpz_mod(t1, t2, n);
36 if (!mpz_cmp_ui(t1, 0) ) {
37 mpz_set_ui(P3->x, 0);
38 mpz_set_ui(P3->y, 1);
39 return;
40 }
41 }
42
43 mpz_sub(t1, P2.x, P1.x);
44 mpz_mod(t2, t1, n);
45
46 /* t1 = 1/deltay mod n */
47 if (!mpz_invert(t1, t2, n)) {
48 /* We've found a factor! In multiply, gcd(mult,n) will be a factor. */
49 mpz_set_ui(P3->x, 0);
50 mpz_set_ui(P3->y, 1);
51 return;
52 }
53
54 mpz_sub(m, P2.y, P1.y);
55 mpz_mod(t2, m, n); /* t2 = deltay mod n */
56 mpz_mul(m, t1, t2);
57 mpz_mod(m, m, n); /* m = deltay / deltax mod n */
58
59 /* x3 = m^2 - x1 - x2 mod n */
60 mpz_mul(t1, m, m);
61 mpz_sub(t2, t1, P1.x);
62 mpz_sub(t1, t2, P2.x);
63 mpz_mod(P3->x, t1, n);
64 /* y3 = m(x1 - x3) - y1 mod n */
65 mpz_sub(t1, P1.x, P3->x);
66 mpz_mul(t2, m, t1);
67 mpz_sub(t1, t2, P1.y);
68 mpz_mod(P3->y, t1, n);
69 }
70
71 /* P3 = 2*P1 */
_ec_add_2A(mpz_t a,mpz_t n,struct ec_affine_point P1,struct ec_affine_point * P3,mpz_t m,mpz_t t1,mpz_t t2)72 static void _ec_add_2A(mpz_t a,
73 mpz_t n,
74 struct ec_affine_point P1,
75 struct ec_affine_point *P3,
76 mpz_t m,
77 mpz_t t1,
78 mpz_t t2)
79 {
80 /* m = (3x1^2 + a) * (2y1)^-1 mod n */
81 mpz_mul_ui(t1, P1.y, 2);
82 if (!mpz_invert(m, t1, n)) {
83 mpz_set_ui(P3->x, 0);
84 mpz_set_ui(P3->y, 1);
85 return;
86 }
87 mpz_mul_ui(t1, P1.x, 3);
88 mpz_mul(t2, t1, P1.x);
89 mpz_add(t1, t2, a);
90 mpz_mul(t2, m, t1);
91 mpz_tdiv_r(m, t2, n);
92
93 /* x3 = m^2 - 2x1 mod n */
94 mpz_mul(t1, m, m);
95 mpz_mul_ui(t2, P1.x, 2);
96 mpz_sub(t1, t1, t2);
97 mpz_tdiv_r(P3->x, t1, n);
98
99 /* y3 = m(x1 - x3) - y1 mod n */
100 mpz_sub(t1, P1.x, P3->x);
101 mpz_mul(t2, t1, m);
102 mpz_sub(t1, t2, P1.y);
103 mpz_tdiv_r(P3->y, t1, n);
104 }
105
ec_affine_multiply(mpz_t a,mpz_t k,mpz_t n,struct ec_affine_point P,struct ec_affine_point * R,mpz_t d)106 int ec_affine_multiply(mpz_t a, mpz_t k, mpz_t n, struct ec_affine_point P, struct ec_affine_point *R, mpz_t d)
107 {
108 int found = 0;
109 struct ec_affine_point A, B, C;
110 mpz_t t, t2, t3, mult;
111
112 mpz_init(A.x); mpz_init(A.y);
113 mpz_init(B.x); mpz_init(B.y);
114 mpz_init(C.x); mpz_init(C.y);
115 mpz_init(t); mpz_init(t2); mpz_init(t3);
116 mpz_init_set_ui(mult, 1); /* holds intermediates, gcd at end */
117
118 mpz_set(A.x, P.x); mpz_set(A.y, P.y);
119 mpz_set_ui(B.x, 0); mpz_set_ui(B.y, 1);
120
121 /* Binary ladder multiply. */
122 while (mpz_cmp_ui(k, 0) > 0) {
123 if (mpz_odd_p(k)) {
124 mpz_sub(t, B.x, A.x);
125 mpz_mul(t2, mult, t);
126 mpz_mod(mult, t2, n);
127
128 if ( !mpz_cmp_ui(A.x, 0) && !mpz_cmp_ui(A.y, 1) ) {
129 /* nothing */
130 } else if ( !mpz_cmp_ui(B.x, 0) && !mpz_cmp_ui(B.y, 1) ) {
131 mpz_set(B.x, A.x); mpz_set(B.y, A.y);
132 } else {
133 _ec_add_AB(n, A, B, &C, t, t2, t3);
134 /* If the add failed to invert, then we have a factor. */
135 mpz_set(B.x, C.x); mpz_set(B.y, C.y);
136 }
137 mpz_sub_ui(k, k, 1);
138 } else {
139 mpz_mul_ui(t, A.y, 2);
140 mpz_mul(t2, mult, t);
141 mpz_mod(mult, t2, n);
142
143 _ec_add_2A(a, n, A, &C, t, t2, t3);
144 mpz_set(A.x, C.x); mpz_set(A.y, C.y);
145 mpz_tdiv_q_2exp(k, k, 1);
146 }
147 }
148 mpz_gcd(d, mult, n);
149 found = (mpz_cmp_ui(d, 1) && mpz_cmp(d, n));
150
151 mpz_tdiv_r(R->x, B.x, n);
152 mpz_tdiv_r(R->y, B.y, n);
153
154 mpz_clear(mult);
155 mpz_clear(t); mpz_clear(t2); mpz_clear(t3);
156 mpz_clear(A.x); mpz_clear(A.y);
157 mpz_clear(B.x); mpz_clear(B.y);
158 mpz_clear(C.x); mpz_clear(C.y);
159
160 return found;
161 }
162
_GMP_ecm_factor_affine(mpz_t n,mpz_t f,UV B1,UV ncurves)163 int _GMP_ecm_factor_affine(mpz_t n, mpz_t f, UV B1, UV ncurves)
164 {
165 mpz_t a, mk;
166 struct ec_affine_point X, Y;
167 UV B, curve, q;
168
169 TEST_FOR_2357(n, f);
170
171 mpz_init(a); mpz_init(mk);
172 mpz_init(X.x); mpz_init(X.y);
173 mpz_init(Y.x); mpz_init(Y.y);
174
175 for (B = 100; B < B1*5; B *= 5) {
176 if (B*5 > 2*B1) B = B1;
177 for (curve = 0; curve < ncurves; curve++) {
178 PRIME_ITERATOR(iter);
179 mpz_isaac_urandomm(a, n);
180 mpz_set_ui(X.x, 0); mpz_set_ui(X.y, 1);
181 for (q = 2; q < B; q = prime_iterator_next(&iter)) {
182 UV k = q;
183 UV kmin = B / q;
184
185 while (k <= kmin)
186 k *= q;
187
188 mpz_set_ui(mk, k);
189 if (ec_affine_multiply(a, mk, n, X, &Y, f)) {
190 prime_iterator_destroy(&iter);
191 mpz_clear(a);
192 mpz_clear(X.x); mpz_clear(X.y);
193 mpz_clear(Y.x); mpz_clear(Y.y);
194 return 1;
195 }
196 mpz_set(X.x, Y.x); mpz_set(X.y, Y.y);
197 /* Check that we're not starting over */
198 if ( !mpz_cmp_ui(X.x, 0) && !mpz_cmp_ui(X.y, 1) )
199 break;
200 }
201 prime_iterator_destroy(&iter);
202 }
203 }
204
205 mpz_clear(a); mpz_clear(mk);
206 mpz_clear(X.x); mpz_clear(X.y);
207 mpz_clear(Y.x); mpz_clear(Y.y);
208
209 return 0;
210 }
211
212
213 /*******************************************************************/
214
215 /* A better ECM, with a stage 2.
216 * Heavily inspired by GMP-ECM, written by Paul Zimmermann (1998),
217 * especially the stage 2 method. Also see "The elliptic curve
218 * integer factorization method" by Bosma and Lenstra as well as many
219 * other articles.
220 */
221
222 static mpz_t b, ecn; /* used throughout ec mult */
223 static mpz_t u, v, w; /* temporaries */
224 static mpz_t x1, z1, x2, z2; /* used by ec_mult and stage2 */
225
226 #define mpz_mulmod(r, a, b, n, t) \
227 do { mpz_mul(t, a, b); mpz_mod(r, t, n); } while (0)
228
229 /* (x2:z2) = (x1:z1) + (x2:z2) */
ec_add(mpz_t x2,mpz_t z2,mpz_t x1,mpz_t z1,mpz_t xinit)230 static void ec_add(mpz_t x2, mpz_t z2, mpz_t x1, mpz_t z1, mpz_t xinit)
231 {
232 mpz_sub(u, x2, z2);
233 mpz_add(v, x1, z1);
234 mpz_mulmod(u, u, v, ecn, w); /* u = (x2 - z2) * (x1 + z1) % n */
235
236 mpz_add(v, x2, z2);
237 mpz_sub(w, x1, z1);
238 mpz_mulmod(v, v, w, ecn, x2); /* v = (x2 + z2) * (x1 - z1) % n */
239
240 mpz_add(w, u, v);
241 mpz_mulmod(x2, w, w, ecn, z2); /* x2 = (u+v)^2 % n */
242
243 mpz_sub(w, u, v);
244 mpz_mulmod(z2, w, w, ecn, v); /* z2 = (u-v)^2 % n */
245
246 mpz_mulmod(z2, xinit, z2, ecn, v); /* z2 *= X1. */
247 /* Per Montgomery 1987, we set Z1 to 1, so no need for x2 *= Z1 */
248 /* 5 mulmods, 6 adds */
249 }
250
251 /* This version assumes no normalization, so uses an extra mulmod. */
252 /* (xout:zout) = (x1:z1) + (x2:z2) */
ec_add3(mpz_t xout,mpz_t zout,mpz_t x1,mpz_t z1,mpz_t x2,mpz_t z2,mpz_t xin,mpz_t zin)253 static void ec_add3(mpz_t xout, mpz_t zout,
254 mpz_t x1, mpz_t z1,
255 mpz_t x2, mpz_t z2,
256 mpz_t xin, mpz_t zin)
257 {
258 mpz_sub(u, x2, z2);
259 mpz_add(v, x1, z1);
260 mpz_mulmod(u, u, v, ecn, w); /* u = (x2 - z2) * (x1 + z1) % n */
261
262 mpz_add(v, x2, z2);
263 mpz_sub(w, x1, z1);
264 mpz_mulmod(v, v, w, ecn, v); /* v = (x2 + z2) * (x1 - z1) % n */
265
266 mpz_add(w, u, v); /* w = u+v */
267 mpz_sub(v, u, v); /* v = u-v */
268
269 mpz_mulmod(w, w, w, ecn, u); /* w = (u+v)^2 % n */
270 mpz_mulmod(v, v, v, ecn, u); /* v = (u-v)^2 % n */
271
272 mpz_set(u, xin);
273 mpz_mulmod(xout, w, zin, ecn, w);
274 mpz_mulmod(zout, v, u, ecn, w);
275 /* 6 mulmods, 6 adds */
276 }
277
278 /* (x2:z2) = 2(x1:z1) */
ec_double(mpz_t x2,mpz_t z2,mpz_t x1,mpz_t z1)279 static void ec_double(mpz_t x2, mpz_t z2, mpz_t x1, mpz_t z1)
280 {
281 mpz_add(u, x1, z1);
282 mpz_mulmod(u, u, u, ecn, w); /* u = (x1+z1)^2 % n */
283
284 mpz_sub(v, x1, z1);
285 mpz_mulmod(v, v, v, ecn, w); /* v = (x1-z1)^2 % n */
286
287 mpz_mulmod(x2, u, v, ecn, w); /* x2 = uv % n */
288
289 mpz_sub(w, u, v); /* w = u-v = 4(x1 * z1) */
290 mpz_mulmod(u, b, w, ecn, z2);
291 mpz_add(u, u, v); /* u = (v+b*w) mod n */
292 mpz_mulmod(z2, w, u, ecn, v); /* z2 = (w*u) mod n */
293 /* 5 mulmods, 4 adds */
294 }
295
296 /* See http://alexandria.tue.nl/extra2/200311829.pdf for lots of discussion
297 * and algorithms for various addition chains.
298 */
299
300 #ifndef USE_PRAC
301
ec_mult(UV k,mpz_t x,mpz_t z)302 static void ec_mult(UV k, mpz_t x, mpz_t z)
303 {
304 int l, r;
305
306 r = --k; l = -1; while (r != 1) { r >>= 1; l++; }
307 if (k & ( UVCONST(1)<<l)) {
308 ec_double(x2, z2, x, z);
309 ec_add3(x1, z1, x2, z2, x, z, x, z);
310 ec_double(x2, z2, x2, z2);
311 } else {
312 ec_double(x1, z1, x, z);
313 ec_add3(x2, z2, x, z, x1, z1, x, z);
314 }
315 l--;
316 while (l >= 1) {
317 if (k & ( UVCONST(1)<<l)) {
318 ec_add3(x1, z1, x1, z1, x2, z2, x, z);
319 ec_double(x2, z2, x2, z2);
320 } else {
321 ec_add3(x2, z2, x2, z2, x1, z1, x, z);
322 ec_double(x1, z1, x1, z1);
323 }
324 l--;
325 }
326 if (k & 1) {
327 ec_double(x, z, x2, z2);
328 } else {
329 ec_add3(x, z, x2, z2, x1, z1, x, z);
330 }
331 }
332
333 #else
334
335 /* PRAC
336 *
337 * "Evaluating recurrences of form X_{m+n} = f(X_m, X_n, X_{m-n}) via
338 * Lucas chains, Peter L. Montgomery, Dec 1983 (revised Jan 1992).
339 * https://cr.yp.to/bib/1992/montgomery-lucas.pdf
340 * http://research.microsoft.com/en-us/um/people/petmon/lucas.pdf
341 *
342 * "20 years of ECM" by Paul Zimmermann, 2006.
343 *
344 * Code derived from GMP-ECM (Zimmermann & Kruppa)
345 */
346
347 /* PRAC, details from GMP-ECM, algorithm from Montgomery */
348 /* See "20 years of ECM" by Paul Zimmermann for more info */
349 static mpz_t x3, z3, x4, z4; /* used by prac */
350 #define ADD 6 /* number of multiplications in an addition */
351 #define DUP 5 /* number of multiplications in a double */
352
353 /* Returns the number of mulmods */
lucas_cost(UV n,double val)354 static UV lucas_cost(UV n, double val)
355 {
356 UV c, d, e, r;
357
358 d = n;
359 r = (UV) ( ((double)d / val) + 0.5 );
360 if (r >=n )
361 return(ADD*n);
362 d = n - r;
363 e = 2 * r - n;
364 c = DUP + ADD; /* initial double and final add */
365 while (d != e) {
366 if (d < e) { r = d; d = e; e = r; }
367 if (4 * d <= 5 * e && ((d + e) % 3) == 0) {
368 /*C1*/ d = (2 * d - e) / 3; e = (e - d) / 2; c += 3 * ADD; /* 3 adds */
369 } else if (4 * d <= 5 * e && (d - e) % 6 == 0) {
370 /*C2*/ d = (d - e) / 2; c += ADD + DUP; /* one add, one double */
371 } else if (d <= 4 * e) {
372 /*C3*/ d -= e; c += ADD; /* one add */
373 } else if ((d + e) % 2 == 0) {
374 /*C4*/ d = (d - e) / 2; c += ADD + DUP; /* one add, one double */
375 } else if (d % 2 == 0) {
376 /*C5*/ d /= 2; c += ADD + DUP; /* one add, one double */
377 } else if (d % 3 == 0) {
378 /*C6*/ d = d / 3 - e; c += 3 * ADD + DUP; /* three adds, one double */
379 } else if ((d + e) % 3 == 0) {
380 /*C7*/ d = (d - 2 * e) / 3; c += 3 * ADD + DUP; /* three adds, one double */
381 } else if ((d - e) % 3 == 0) {
382 /*C8*/ d = (d - e) / 3; c += 3 * ADD + DUP; /* three adds, one double */
383 } else {
384 /*C9*/ e /= 2; c += ADD + DUP; /* one add, one double */
385 }
386 }
387 return(c);
388 }
389
390 #define SWAP(a, b) \
391 t = x##a; x##a = x##b; x##b = t; t = z##a; z##a = z##b; z##b = t;
392
393 /* PRAC: computes kP from P=(x:z) and puts the result in (x:z). Assumes k>2.*/
ec_mult(UV k,mpz_t x,mpz_t z)394 static void ec_mult(UV k, mpz_t x, mpz_t z)
395 {
396 unsigned int d, e, r, i;
397 __mpz_struct *xA, *zA, *xB, *zB, *xC, *zC, *xT, *zT, *xT2, *zT2, *t;
398
399 static double const val[] =
400 {1.61803398875, 1.72360679775, 1.618347119656, 1.617914406529,
401 1.58017872826};
402
403 /* chooses the best value of val */
404 r = ADD * k;
405 i = 0;
406 for (d = 0; d < 5; d++) {
407 e = lucas_cost(k, val[d]);
408 if (e < r) { r = e; i = d; }
409 }
410 r = (unsigned int)((double)k / val[i] + 0.5);
411 /* A=(x:z) B=(x1:z1) C=(x2:z2) T=T1=(x3:z3) T2=(x4:z4) */
412 xA=x; zA=z; xB=x1; zB=z1; xC=x2; zC=z2; xT=x3; zT=z3; xT2=x4; zT2=z4;
413 /* first iteration always begins by Condition 3, then a swap */
414 d = k - r;
415 e = 2 * r - k;
416 mpz_set(xB,xA); mpz_set(zB,zA); /* B=A */
417 mpz_set(xC,xA); mpz_set(zC,zA); /* C=A */
418 ec_double(xA,zA,xA,zA); /* A=2*A */
419 while (d != e) {
420 if (d < e) {
421 r = d; d = e; e = r;
422 SWAP(A,B);
423 }
424 /* do the first line of Table 4 whose condition qualifies */
425 if (4 * d <= 5 * e && ((d + e) % 3) == 0) { /* condition 1 */
426 d = (2 * d - e) / 3;
427 e = (e - d) / 2;
428 ec_add3(xT,zT,xA,zA,xB,zB,xC,zC); /* T = f(A,B,C) */
429 ec_add3(xT2,zT2,xT,zT,xA,zA,xB,zB); /* T2= f(T,A,B) */
430 ec_add3(xB,zB,xB,zB,xT,zT,xA,zA); /* B = f(B,T,A) */
431 SWAP(A,T2);
432 } else if (4 * d <= 5 * e && (d - e) % 6 == 0) { /* condition 2 */
433 d = (d - e) / 2;
434 ec_add3(xB,zB,xA,zA,xB,zB,xC,zC); /* B = f(A,B,C) */
435 ec_double(xA,zA,xA,zA); /* A = 2*A */
436 } else if (d <= (4 * e)) { /* condition 3 */
437 d -= e;
438 ec_add3(xC,zC,xB,zB,xA,zA,xC,zC); /* C = f(B,A,C) */
439 SWAP(B,C);
440 } else if ((d + e) % 2 == 0) { /* condition 4 */
441 d = (d - e) / 2;
442 ec_add3(xB,zB,xB,zB,xA,zA,xC,zC); /* B = f(B,A,C) */
443 ec_double(xA,zA,xA,zA); /* A = 2*A */
444 } else if (d % 2 == 0) { /* condition 5 */
445 d /= 2;
446 ec_add3(xC,zC,xC,zC,xA,zA,xB,zB); /* C = f(C,A,B) */
447 ec_double(xA,zA,xA,zA); /* A = 2*A */
448 } else if (d % 3 == 0) { /* condition 6 */
449 d = d / 3 - e;
450 ec_double(xT,zT,xA,zA); /* T = 2*A */
451 ec_add3(xT2,zT2,xA,zA,xB,zB,xC,zC); /* T2= f(A,B,C) */
452 ec_add3(xA,zA,xT,zT,xA,zA,xA,zA); /* A = f(T,A,A) */
453 ec_add3(xC,zC,xT,zT,xT2,zT2,xC,zC); /* C = f(T,T2,C) */
454 SWAP(B,C);
455 } else if ((d + e) % 3 == 0) { /* condition 7 */
456 d = (d - 2 * e) / 3;
457 ec_add3(xT,zT,xA,zA,xB,zB,xC,zC); /* T = f(A,B,C) */
458 ec_add3(xB,zB,xT,zT,xA,zA,xB,zB); /* B = f(T1,A,B) */
459 ec_double(xT,zT,xA,zA);
460 ec_add3(xA,zA,xA,zA,xT,zT,xA,zA); /* A = 3*A */
461 } else if ((d - e) % 3 == 0) { /* condition 8 */
462 d = (d - e) / 3;
463 ec_add3(xT,zT,xA,zA,xB,zB,xC,zC); /* T = f(A,B,C) */
464 ec_add3(xC,zC,xC,zC,xA,zA,xB,zB); /* C = f(A,C,B) */
465 SWAP(B,T);
466 ec_double(xT,zT,xA,zA);
467 ec_add3(xA,zA,xA,zA,xT,zT,xA,zA); /* A = 3*A */
468 } else { /* condition 9 */
469 e /= 2;
470 ec_add3(xC,zC,xC,zC,xB,zB,xA,zA); /* C = f(C,B,A) */
471 ec_double(xB,zB,xB,zB); /* B = 2*B */
472 }
473 }
474 ec_add3(xA,zA,xA,zA,xB,zB,xC,zC);
475 if (x!=xA) { mpz_set(x,xA); mpz_set(z,zA); }
476 }
477
478 #endif /* PRAC */
479
480 #define NORMALIZE(f, u, v, x, z, n) \
481 mpz_gcdext(f, u, NULL, z, n); \
482 found = mpz_cmp_ui(f, 1); \
483 if (found) break; \
484 mpz_mulmod(x, x, u, n, v); \
485 mpz_set_ui(z, 1);
486
ec_stage2(UV B1,UV B2,mpz_t x,mpz_t z,mpz_t f)487 static int ec_stage2(UV B1, UV B2, mpz_t x, mpz_t z, mpz_t f)
488 {
489 UV D, i, m;
490 mpz_t* nqx = 0;
491 mpz_t g, one;
492 int found;
493 PRIME_ITERATOR(iter);
494
495 do {
496 NORMALIZE(f, u, v, x, z, ecn);
497
498 D = sqrt( (double)B2 / 2.0 );
499 if (D%2) D++;
500
501 /* We really only need half of these. Only even values used. */
502 Newz(0, nqx, 2*D+1, mpz_t);
503 mpz_init_set(nqx[1], x);
504 mpz_init_set_ui(g, 1);
505 mpz_init_set_ui(one, 1);
506
507 for (i = 2; i <= 2*D; i++) {
508 if (i % 2) {
509 mpz_set(x2, nqx[(i+1)/2]); mpz_set_ui(z2, 1);
510 ec_add(x2, z2, nqx[(i-1)/2], one, x);
511 } else {
512 ec_double(x2, z2, nqx[i/2], one);
513 }
514 mpz_init_set(nqx[i], x2);
515 NORMALIZE(f, u, v, nqx[i], z2, ecn);
516 }
517 if (found) break;
518
519 mpz_set(x1, x);
520 mpz_set(z1, z);
521 mpz_set(x, nqx[2*D-1]);
522 mpz_set_ui(z, 1);
523
524 /* See Zimmermann, "20 Years of ECM" slides, 2006, page 11-12 */
525 for (m = 1; m < B2+D; m += 2*D) {
526 if (m != 1) {
527 mpz_set(x2, x1);
528 mpz_set(z2, z1);
529 ec_add(x1, z1, nqx[2*D], one, x);
530 NORMALIZE(f, u, v, x1, z1, ecn);
531 mpz_set(x, x2); mpz_set(z, z2);
532 }
533 if (m+D > B1 && m >= D) {
534 prime_iterator_setprime(&iter, m-D-1);
535 for (i = prime_iterator_next(&iter); i < m; i = prime_iterator_next(&iter)) {
536 /* if (m+D-i<1 || m+D-i>2*D) croak("index %lu range\n",i-(m-D)); */
537 mpz_sub(w, x1, nqx[m+D-i]);
538 mpz_mulmod(g, g, w, ecn, u);
539 }
540 for ( ; i <= m+D; i = prime_iterator_next(&iter)) {
541 if (i > m && !prime_iterator_isprime(&iter, m+m-i)) {
542 /* if (i-m<1 || i-m>2*D) croak("index %lu range\n",i-(m-D)); */
543 mpz_sub(w, x1, nqx[i-m]);
544 mpz_mulmod(g, g, w, ecn, u);
545 }
546 }
547 mpz_gcd(f, g, ecn);
548 found = mpz_cmp_ui(f, 1);
549 if (found) break;
550 }
551 }
552 } while (0);
553 prime_iterator_destroy(&iter);
554
555 if (nqx != 0) {
556 for (i = 1; i <= 2*D; i++) {
557 if (nqx[i] != 0)
558 mpz_clear(nqx[i]);
559 }
560 Safefree(nqx);
561 mpz_clear(g);
562 mpz_clear(one);
563 }
564 if (found && !mpz_cmp(f, ecn)) found = 0;
565 return (found) ? 2 : 0;
566 }
567
_GMP_ecm_factor_projective(mpz_t n,mpz_t f,UV B1,UV B2,UV ncurves)568 int _GMP_ecm_factor_projective(mpz_t n, mpz_t f, UV B1, UV B2, UV ncurves)
569 {
570 mpz_t sigma, a, x, z;
571 UV i, curve, q, k;
572 int found = 0;
573 int _verbose = get_verbose_level();
574
575 TEST_FOR_2357(n, f);
576
577 if (B2 < B1) B2 = 100*B1; /* time(S1) == time(S2) ~ 125 */
578
579 mpz_init_set(ecn, n);
580 mpz_init(x); mpz_init(z); mpz_init(a); mpz_init(b);
581 mpz_init(u); mpz_init(v); mpz_init(w); mpz_init(sigma);
582 mpz_init(x1); mpz_init(z1); mpz_init(x2); mpz_init(z2);
583 #ifdef USE_PRAC
584 mpz_init(x3); mpz_init(z3); mpz_init(x4); mpz_init(z4);
585 #endif
586
587 if (_verbose>2) gmp_printf("# ecm trying %Zd (B1=%lu B2=%lu ncurves=%lu)\n", n, (unsigned long)B1, (unsigned long)B2, (unsigned long)ncurves);
588
589 for (curve = 0; curve < ncurves; curve++) {
590 PRIME_ITERATOR(iter);
591 do {
592 mpz_isaac_urandomm(sigma, n);
593 } while (mpz_cmp_ui(sigma, 5) <= 0);
594 mpz_mul_ui(w, sigma, 4);
595 mpz_mod(v, w, n); /* v = 4σ */
596
597 mpz_mul(x, sigma, sigma);
598 mpz_sub_ui(w, x, 5);
599 mpz_mod(u, w, n); /* u = σ^2-5 */
600
601 mpz_mul(x, u, u);
602 mpz_mulmod(x, x, u, n, w); /* x = u^3 */
603
604 mpz_mul(z, v, v);
605 mpz_mulmod(z, z, v, n, w); /* z = v^3 */
606
607 mpz_mul(b, x, v);
608 mpz_mul_ui(w, b, 4);
609 mpz_mod(b, w, n); /* b = 4 u^3 v */
610
611 mpz_sub(a, v, u);
612 mpz_mul(w, a, a);
613 mpz_mulmod(w, w, a, n, w);
614
615 mpz_mul_ui(a, u, 3);
616 mpz_add(a, a, v);
617 mpz_mul(w, w, a);
618 mpz_mod(a, w, n); /* a = ((v-u)^3 * (3*u + v)) % n */
619
620 mpz_gcdext(f, u, NULL, b, n);
621 found = mpz_cmp_ui(f, 1);
622 if (found) { if (!mpz_cmp(f, n)) { found = 0; continue; } break; }
623 mpz_mul(a, a, u);
624
625 mpz_sub_ui(a, a, 2);
626 mpz_mod(a, a, n);
627
628 mpz_add_ui(b, a, 2);
629 if (mpz_mod_ui(w, b, 2)) mpz_add(b, b, n);
630 mpz_tdiv_q_2exp(b, b, 1);
631 if (mpz_mod_ui(w, b, 2)) mpz_add(b, b, n);
632 mpz_tdiv_q_2exp(b, b, 1);
633
634 /* Use sigma to collect possible factors */
635 mpz_set_ui(sigma, 1);
636
637 /* Stage 1 */
638 for (q = 2; q < B1; q *= 2)
639 ec_double(x, z, x, z);
640 mpz_mulmod(sigma, sigma, x, ecn, w);
641 i = 15;
642 for (q = prime_iterator_next(&iter); q < B1; q = prime_iterator_next(&iter)) {
643 /* PRAC is a little faster with:
644 * for (k = 1; k <= B1/q; k *= q)
645 * ec_mult(q, x, z);
646 * but binary multiplication is much slower that way. */
647 for (k = q; k <= B1/q; k *= q) ;
648 ec_mult(k, x, z);
649 mpz_mulmod(sigma, sigma, x, ecn, w);
650 if (i++ % 32 == 0) {
651 mpz_gcd(f, sigma, ecn);
652 if (mpz_cmp_ui(f, 1)) break;
653 }
654 }
655 prime_iterator_destroy(&iter);
656
657 /* Find factor in S1 */
658 do { NORMALIZE(f, u, v, x, z, n); } while (0);
659 if (!found) {
660 mpz_gcd(f, sigma, ecn);
661 found = mpz_cmp_ui(f, 1);
662 }
663 if (found) { if (!mpz_cmp(f, n)) { found = 0; continue; } break; }
664
665 /* Stage 2 */
666 if (!found && B2 > B1)
667 found = ec_stage2(B1, B2, x, z, f);
668
669 if (found) { if (!mpz_cmp(f, n)) { found = 0; continue; } break; }
670 }
671 if (_verbose>2) {
672 if (found) gmp_printf("# ecm: %Zd in stage %d\n", f, found);
673 else gmp_printf("# ecm: no factor\n");
674 }
675
676 mpz_clear(ecn);
677 mpz_clear(x); mpz_clear(z); mpz_clear(a); mpz_clear(b);
678 mpz_clear(u); mpz_clear(v); mpz_clear(w); mpz_clear(sigma);
679 mpz_clear(x1); mpz_clear(z1); mpz_clear(x2); mpz_clear(z2);
680 #ifdef USE_PRAC
681 mpz_clear(x3); mpz_clear(z3); mpz_clear(x4); mpz_clear(z4);
682 #endif
683
684 return found;
685 }
686