/* Elliptic Curve Method: toplevel and stage 1 routines. Copyright 2001, 2002, 2003, 2004, 2005 Paul Zimmermann and Alexander Kruppa. This file is part of the ECM Library. The ECM Library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation; either version 2.1 of the License, or (at your option) any later version. The ECM Library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. You should have received a copy of the GNU Lesser General Public License along with the ECM Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ #include #include #include "ecm-impl.h" #include #if HAVE_LIMITS_H # include #else # define ULONG_MAX __GMP_ULONG_MAX #endif /****************************************************************************** * * * Elliptic Curve Method * * * ******************************************************************************/ #define mpz_mulmod5(r,s1,s2,m,t) { mpz_mul(t,s1,s2); mpz_mod(r, t, m); } /* Computes curve parameter A and a starting point (x:1) from a given sigma value. If a factor of n was found during the process, returns 1 (and factor in f), 0 otherwise. */ static int get_curve_from_sigma (mpz_t f, mpres_t A, mpres_t x, mpz_t sigma, mpmod_t n) { mpres_t t, u, v, b, z; MEMORY_TAG; mpres_init (t, n); MEMORY_TAG; mpres_init (u, n); MEMORY_TAG; mpres_init (v, n); MEMORY_TAG; mpres_init (b, n); MEMORY_TAG; mpres_init (z, n); MEMORY_UNTAG; mpres_set_z (u, sigma, n); mpres_mul_ui (v, u, 4, n); /* v = (4*sigma) mod n */ mpres_mul (t, u, u, n); mpres_sub_ui (u, t, 5, n); /* u = (sigma^2-5) mod n */ mpres_mul (t, u, u, n); mpres_mul (x, t, u, n); /* x = (u^3) mod n */ mpres_mul (t, v, v, n); mpres_mul (z, t, v, n); /* z = (v^3) mod n */ mpres_mul (t, x, v, n); mpres_mul_ui (b, t, 4, n); /* b = (4*x*v) mod n */ mpres_mul_ui (t, u, 3, n); mpres_sub (u, v, u, n); /* u' = v-u */ mpres_add (v, t, v, n); /* v' = (3*u+v) mod n */ mpres_mul (t, u, u, n); mpres_mul (u, t, u, n); /* u'' = ((v-u)^3) mod n */ mpres_mul (A, u, v, n); /* a = (u'' * v') mod n = ((v-u)^3 * (3*u+v)) mod n */ /* Normalize b and z to 1 */ mpres_mul (v, b, z, n); if (!mpres_invert (u, v, n)) /* u = (b*z)^(-1) (mod n) */ { mpres_gcd (f, v, n); mpres_clear (t, n); mpres_clear (u, n); mpres_clear (v, n); mpres_clear (b, n); mpres_clear (z, n); return 1; } mpres_mul (v, u, b, n); /* v = z^(-1) (mod n) */ mpres_mul (x, x, v, n); /* x = x * z^(-1) */ mpres_mul (v, u, z, n); /* v = b^(-1) (mod n) */ mpres_mul (t, A, v, n); mpres_sub_ui (A, t, 2, n); mpres_clear (t, n); mpres_clear (u, n); mpres_clear (v, n); mpres_clear (b, n); mpres_clear (z, n); return 0; } /* switch from Montgomery's form g*y^2 = x^3 + a*x^2 + x to Weierstrass' form Y^2 = X^3 + A*X + B by change of variables x -> g*X-a/3, y -> g*Y. We have A = (3-a^2)/(3g^2), X = (3x+a)/(3g), Y = y/g. */ static int montgomery_to_weierstrass (mpz_t f, mpres_t x, mpres_t y, mpres_t A, mpmod_t n) { mpres_t g; MEMORY_TAG; mpres_init (g, n); MEMORY_UNTAG; mpres_add (g, x, A, n); mpres_mul (g, g, x, n); mpres_add_ui (g, g, 1, n); mpres_mul (g, g, x, n); /* g = x^3+a*x^2+x (y=1) */ mpres_mul_ui (y, g, 3, n); mpres_mul (y, y, g, n); /* y = 3g^2 */ if (!mpres_invert (y, y, n)) /* y = 1/(3g^2) temporarily */ { mpres_gcd (f, y, n); mpres_clear (g, n); return 1; } /* update x */ mpres_mul_ui (x, x, 3, n); /* 3x */ mpres_add (x, x, A, n); /* 3x+a */ mpres_mul (x, x, g, n); /* (3x+a)*g */ mpres_mul (x, x, y, n); /* (3x+a)/(3g) */ /* update A */ mpres_mul (A, A, A, n); /* a^2 */ mpres_sub_ui (A, A, 3, n); mpres_neg (A, A, n); /* 3-a^2 */ mpres_mul (A, A, y, n); /* (3-a^2)/(3g^2) */ /* update y */ mpres_mul_ui (g, g, 3, n); /* 3g */ mpres_mul (y, y, g, n); /* (3g)/(3g^2) = 1/g */ mpres_clear (g, n); return 0; } /* adds Q=(x2:z2) and R=(x1:z1) and puts the result in (x3:z3), using 6 muls (4 muls and 2 squares), and 6 add/sub. One assumes that Q-R=P or R-Q=P where P=(x:z). - n : number to factor - u, v, w : auxiliary variables Modifies: x3, z3, u, v, w. (x3,z3) may be identical to (x2,z2) and to (x,z) */ void add3 (mpres_t x3, mpres_t z3, mpres_t x2, mpres_t z2, mpres_t x1, mpres_t z1, mpres_t x, mpres_t z, mpmod_t n, mpres_t u, mpres_t v, mpres_t w) { mpres_sub (u, x2, z2, n); mpres_add (v, x1, z1, n); /* u = x2-z2, v = x1+z1 */ mpres_mul (u, u, v, n); /* u = (x2-z2)*(x1+z1) */ mpres_add (w, x2, z2, n); mpres_sub (v, x1, z1, n); /* w = x2+z2, v = x1-z1 */ mpres_mul (v, w, v, n); /* v = (x2+z2)*(x1-z1) */ mpres_add (w, u, v, n); /* w = 2*(x1*x2-z1*z2) */ mpres_sub (v, u, v, n); /* v = 2*(x2*z1-x1*z2) */ mpres_mul (w, w, w, n); /* w = 4*(x1*x2-z1*z2)^2 */ mpres_mul (v, v, v, n); /* v = 4*(x2*z1-x1*z2)^2 */ if (x == x3) /* same variable: in-place variant */ { /* x3 <- w * z mod n z3 <- x * v mod n */ mpres_mul (z3, w, z, n); mpres_mul (x3, x, v, n); mpres_swap (x3, z3, n); } else { mpres_mul (x3, w, z, n); /* x3 = 4*z*(x1*x2-z1*z2)^2 mod n */ mpres_mul (z3, x, v, n); /* z3 = 4*x*(x2*z1-x1*z2)^2 mod n */ } /* mul += 6; */ } /* computes 2P=(x2:z2) from P=(x1:z1), with 5 muls (3 muls and 2 squares) and 4 add/sub. - n : number to factor - b : (a+2)/4 mod n - t, u, v, w : auxiliary variables */ void duplicate (mpres_t x2, mpres_t z2, mpres_t x1, mpres_t z1, mpmod_t n, mpres_t b, mpres_t u, mpres_t v, mpres_t w) { mpres_add (u, x1, z1, n); mpres_mul (u, u, u, n); /* u = (x1+z1)^2 mod n */ mpres_sub (v, x1, z1, n); mpres_mul (v, v, v, n); /* v = (x1-z1)^2 mod n */ mpres_mul (x2, u, v, n); /* x2 = u*v = (x1^2 - z1^2)^2 mod n */ mpres_sub (w, u, v, n); /* w = u-v = 4*x1*z1 */ mpres_mul (u, w, b, n); /* u = w*b = ((A+2)/4*(4*x1*z1)) mod n */ mpres_add (u, u, v, n); /* u = (x1-z1)^2+(A+2)/4*(4*x1*z1) */ mpres_mul (z2, w, u, n); /* z2 = ((4*x1*z1)*((x1-z1)^2+(A+2)/4*(4*x1*z1))) mod n */ } /* multiply P=(x:z) by e and puts the result in (x:z). */ void ecm_mul (mpres_t x, mpres_t z, mpz_t e, mpmod_t n, mpres_t b) { size_t l; int negated = 0; mpres_t x0, z0, x1, z1, u, v, w; /* In Montgomery coordinates, the point at infinity is (0::0) */ if (mpz_sgn (e) == 0) { mpz_set_ui (x, 0); mpz_set_ui (z, 0); return; } /* The negative of a point (x:y:z) is (x:-y:u). Since we do not compute y, e*(x::z) == (-e)*(x::z). */ if (mpz_sgn (e) < 0) { negated = 1; mpz_neg (e, e); } if (mpz_cmp_ui (e, 1) == 0) goto ecm_mul_end; MEMORY_TAG; mpres_init (x0, n); MEMORY_TAG; mpres_init (z0, n); MEMORY_TAG; mpres_init (x1, n); MEMORY_TAG; mpres_init (z1, n); MEMORY_TAG; mpres_init (u, n); MEMORY_TAG; mpres_init (v, n); MEMORY_TAG; mpres_init (w, n); MEMORY_UNTAG; l = mpz_sizeinbase (e, 2) - 1; /* l >= 1 */ mpres_set (x0, x, n); mpres_set (z0, z, n); duplicate (x1, z1, x0, z0, n, b, u, v, w); /* invariant: (P1,P0) = ((k+1)P, kP) where k = floor(e/2^l) */ while (l-- > 0) { if (mpz_tstbit (e, l)) /* k, k+1 -> 2k+1, 2k+2 */ { add3 (x0, z0, x0, z0, x1, z1, x, z, n, u, v, w); /* 2k+1 */ duplicate (x1, z1, x1, z1, n, b, u, v, w); /* 2k+2 */ } else /* k, k+1 -> 2k, 2k+1 */ { add3 (x1, z1, x1, z1, x0, z0, x, z, n, u, v, w); /* 2k+1 */ duplicate (x0, z0, x0, z0, n, b, u, v, w); /* 2k */ } } mpres_set (x, x0, n); mpres_set (z, z0, n); mpres_clear (x0, n); mpres_clear (z0, n); mpres_clear (x1, n); mpres_clear (z1, n); mpres_clear (u, n); mpres_clear (v, n); mpres_clear (w, n); ecm_mul_end: /* Undo negation to avoid changing the caller's e value */ if (negated) mpz_neg (e, e); } #define ADD 6.0 /* number of multiplications in an addition */ #define DUP 5.0 /* number of multiplications in a duplicate */ /* returns the number of modular multiplications for computing V_n from V_r * V_{n-r} - V_{n-2r}. ADD is the cost of an addition DUP is the cost of a duplicate */ static double lucas_cost (unsigned long n, double v) { unsigned long d, e, r; double c; /* cost */ d = n; r = (unsigned long) ((double) d / v + 0.5); if (r >= n) return (ADD * (double) n); d = n - r; e = 2 * r - n; c = DUP + ADD; /* initial duplicate and final addition */ while (d != e) { if (d < e) { r = d; d = e; e = r; } if (4 * d <= 5 * e && ((d + e) % 3) == 0) { /* condition 1 */ d = (2 * d - e) / 3; e = (e - d) / 2; c += 3.0 * ADD; /* 3 additions */ } else if (4 * d <= 5 * e && (d - e) % 6 == 0) { /* condition 2 */ d = (d - e) / 2; c += ADD + DUP; /* one addition, one duplicate */ } else if (d <= 4 * e) { /* condition 3 */ d -= e; c += ADD; /* one addition */ } else if ((d + e) % 2 == 0) { /* condition 4 */ d = (d - e) / 2; c += ADD + DUP; /* one addition, one duplicate */ } /* now d+e is odd */ else if (d % 2 == 0) { /* condition 5 */ d /= 2; c += ADD + DUP; /* one addition, one duplicate */ } /* now d is odd and e is even */ else if (d % 3 == 0) { /* condition 6 */ d = d / 3 - e; c += 3.0 * ADD + DUP; /* three additions, one duplicate */ } else if ((d + e) % 3 == 0) { /* condition 7 */ d = (d - 2 * e) / 3; c += 3.0 * ADD + DUP; /* three additions, one duplicate */ } else if ((d - e) % 3 == 0) { /* condition 8 */ d = (d - e) / 3; c += 3.0 * ADD + DUP; /* three additions, one duplicate */ } else /* necessarily e is even: catches all cases */ { /* condition 9 */ e /= 2; c += ADD + DUP; /* one addition, one duplicate */ } } return c; } /* computes kP from P=(xA:zA) and puts the result in (xA:zA). Assumes k>2. WARNING! The calls to add3() assume that the two input points are distinct, which is not neccessarily satisfied. The result can be that in rare cases the point at infinity (z==0) results when it shouldn't. A test case is echo 33554520197234177 | ./ecm -sigma 2046841451 373 1 which finds the prime even though it shouldn't (23^2=529 divides order). This is not a problem for ECM since at worst we'll find a factor we shouldn't have found. For other purposes (i.e. primality proving) this would have to be fixed first. */ static void prac (mpres_t xA, mpres_t zA, unsigned long k, mpmod_t n, mpres_t b, mpres_t u, mpres_t v, mpres_t w, mpres_t xB, mpres_t zB, mpres_t xC, mpres_t zC, mpres_t xT, mpres_t zT, mpres_t xT2, mpres_t zT2) { unsigned long d, e, r, i = 0; double c, cmin; __mpz_struct *tmp; #define NV 10 static double val[NV] = { 1.6180339887498948, 1.7236067977499790, 1.6183471196562281, 1.6179144065288179, 1.6124299495094950, 1.6328398060887063, 1.6201819808074158, 1.5801787282954641, 1.6172146165344039, 1.3819660112501052 }; /* chooses the best value of v */ for (d = 0, cmin = ADD * (double) k; d < NV; d++) { c = lucas_cost (k, val[d]); if (c < cmin) { cmin = c; i = d; } } d = k; r = (unsigned long) ((double) d / val[i] + 0.5); /* first iteration always begins by Condition 3, then a swap */ d = k - r; e = 2 * r - k; mpres_set (xB, xA, n); mpres_set (zB, zA, n); /* B=A */ mpres_set (xC, xA, n); mpres_set (zC, zA, n); /* C=A */ duplicate (xA, zA, xA, zA, n, b, u, v, w); /* A = 2*A */ while (d != e) { if (d < e) { r = d; d = e; e = r; mpres_swap (xA, xB, n); mpres_swap (zA, zB, n); } /* do the first line of Table 4 whose condition qualifies */ if (4 * d <= 5 * e && ((d + e) % 3) == 0) { /* condition 1 */ d = (2 * d - e) / 3; e = (e - d) / 2; add3 (xT, zT, xA, zA, xB, zB, xC, zC, n, u, v, w); /* T = f(A,B,C) */ add3 (xT2, zT2, xT, zT, xA, zA, xB, zB, n, u, v, w); /* T2 = f(T,A,B) */ add3 (xB, zB, xB, zB, xT, zT, xA, zA, n, u, v, w); /* B = f(B,T,A) */ mpres_swap (xA, xT2, n); mpres_swap (zA, zT2, n); /* swap A and T2 */ } else if (4 * d <= 5 * e && (d - e) % 6 == 0) { /* condition 2 */ d = (d - e) / 2; add3 (xB, zB, xA, zA, xB, zB, xC, zC, n, u, v, w); /* B = f(A,B,C) */ duplicate (xA, zA, xA, zA, n, b, u, v, w); /* A = 2*A */ } else if (d <= (4 * e)) { /* condition 3 */ d -= e; add3 (xT, zT, xB, zB, xA, zA, xC, zC, n, u, v, w); /* T = f(B,A,C) */ /* circular permutation (B,T,C) */ tmp = xB; xB = xT; xT = xC; xC = tmp; tmp = zB; zB = zT; zT = zC; zC = tmp; } else if ((d + e) % 2 == 0) { /* condition 4 */ d = (d - e) / 2; add3 (xB, zB, xB, zB, xA, zA, xC, zC, n, u, v, w); /* B = f(B,A,C) */ duplicate (xA, zA, xA, zA, n, b, u, v, w); /* A = 2*A */ } /* now d+e is odd */ else if (d % 2 == 0) { /* condition 5 */ d /= 2; add3 (xC, zC, xC, zC, xA, zA, xB, zB, n, u, v, w); /* C = f(C,A,B) */ duplicate (xA, zA, xA, zA, n, b, u, v, w); /* A = 2*A */ } /* now d is odd, e is even */ else if (d % 3 == 0) { /* condition 6 */ d = d / 3 - e; duplicate (xT, zT, xA, zA, n, b, u, v, w); /* T = 2*A */ add3 (xT2, zT2, xA, zA, xB, zB, xC, zC, n, u, v, w); /* T2 = f(A,B,C) */ add3 (xA, zA, xT, zT, xA, zA, xA, zA, n, u, v, w); /* A = f(T,A,A) */ add3 (xT, zT, xT, zT, xT2, zT2, xC, zC, n, u, v, w); /* T = f(T,T2,C) */ /* circular permutation (C,B,T) */ tmp = xC; xC = xB; xB = xT; xT = tmp; tmp = zC; zC = zB; zB = zT; zT = tmp; } else if ((d + e) % 3 == 0) { /* condition 7 */ d = (d - 2 * e) / 3; add3 (xT, zT, xA, zA, xB, zB, xC, zC, n, u, v, w); /* T = f(A,B,C) */ add3 (xB, zB, xT, zT, xA, zA, xB, zB, n, u, v, w); /* B = f(T,A,B) */ duplicate (xT, zT, xA, zA, n, b, u, v, w); add3 (xA, zA, xA, zA, xT, zT, xA, zA, n, u, v, w); /* A = 3*A */ } else if ((d - e) % 3 == 0) { /* condition 8 */ d = (d - e) / 3; add3 (xT, zT, xA, zA, xB, zB, xC, zC, n, u, v, w); /* T = f(A,B,C) */ add3 (xC, zC, xC, zC, xA, zA, xB, zB, n, u, v, w); /* C = f(A,C,B) */ mpres_swap (xB, xT, n); mpres_swap (zB, zT, n); /* swap B and T */ duplicate (xT, zT, xA, zA, n, b, u, v, w); add3 (xA, zA, xA, zA, xT, zT, xA, zA, n, u, v, w); /* A = 3*A */ } else /* necessarily e is even here */ { /* condition 9 */ e /= 2; add3 (xC, zC, xC, zC, xB, zB, xA, zA, n, u, v, w); /* C = f(C,B,A) */ duplicate (xB, zB, xB, zB, n, b, u, v, w); /* B = 2*B */ } } add3 (xA, zA, xA, zA, xB, zB, xC, zC, n, u, v, w); ASSERT(d == 1); } /* Input: x is initial point A is curve parameter in Montgomery's form: g*y^2*z = x^3 + a*x^2*z + x*z^2 n is the number to factor B1 is the stage 1 bound Output: If a factor is found, it is returned in x. Otherwise, x contains the x-coordinate of the point computed in stage 1 (with z coordinate normalized to 1). Return value: non-zero iff a factor is found */ static int ecm_stage1 (mpz_t f, mpres_t x, mpres_t A, mpmod_t n, double B1, double *B1done, mpz_t go, int (*stop_asap)(void)) { mpres_t b, z, u, v, w, xB, zB, xC, zC, xT, zT, xT2, zT2; double q, r; int ret = 0; MEMORY_TAG; mpres_init (b, n); MEMORY_TAG; mpres_init (z, n); MEMORY_TAG; mpres_init (u, n); MEMORY_TAG; mpres_init (v, n); MEMORY_TAG; mpres_init (w, n); MEMORY_TAG; mpres_init (xB, n); MEMORY_TAG; mpres_init (zB, n); MEMORY_TAG; mpres_init (xC, n); MEMORY_TAG; mpres_init (zC, n); MEMORY_TAG; mpres_init (xT, n); MEMORY_TAG; mpres_init (zT, n); MEMORY_TAG; mpres_init (xT2, n); MEMORY_TAG; mpres_init (zT2, n); MEMORY_UNTAG; mpres_set_ui (z, 1, n); mpres_add_ui (b, A, 2, n); mpres_div_2exp (b, b, 2, n); /* b == (A0+2)*B/4, where B=2^(k*GMP_NUMB_LIMB) for MODMULN or REDC, B=1 otherwise */ #ifndef FULL_REDUCTION mpres_semi_normalize (b, mpz_size (n->orig_modulus)); #endif /* preload group order */ if (go != NULL) ecm_mul (x, z, go, n, b); /* prac() wants multiplicands > 2 */ for (r = 2.0; r <= B1; r *= 2.0) if (r > *B1done) duplicate (x, z, x, z, n, b, u, v, w); /* We'll do 3 manually, too (that's what ecm4 did..) */ for (r = 3.0; r <= B1; r *= 3.0) if (r > *B1done) { duplicate (xB, zB, x, z, n, b, u, v, w); add3 (x, z, x, z, xB, zB, x, z, n, u, v, w); } q = getprime (2.0); /* Puts 3.0 into q. Next call gives 5.0 */ for (q = getprime (q); q <= B1; q = getprime (q)) { for (r = q; r <= B1; r *= q) if (r > *B1done) prac (x, z, (unsigned long) q, n, b, u, v, w, xB, zB, xC, zC, xT, zT, xT2, zT2); if (mpres_is_zero (z, n)) { outputf (OUTPUT_VERBOSE, "Reached point at infinity, %.0f divides " "group order\n", q); break; } if (stop_asap != NULL && (*stop_asap) ()) { outputf (OUTPUT_NORMAL, "Interrupted at prime %.0f\n", q); break; } } getprime (FREE_PRIME_TABLE); /* free the prime tables, and reinitialize */ *B1done = (q < B1) ? q : B1; /* Normalize z to 1 */ #ifndef FULL_REDUCTION mpres_normalize (z); /* needed for gcd */ #endif if (!mpres_invert (u, z, n)) /* Factor found? */ { mpres_gcd (f, z, n); ret = 1; } mpres_mul (x, x, u, n); mpres_clear (zT2, n); mpres_clear (xT2, n); mpres_clear (zT, n); mpres_clear (xT, n); mpres_clear (zC, n); mpres_clear (xC, n); mpres_clear (zB, n); mpres_clear (xB, n); mpres_clear (w, n); mpres_clear (v, n); mpres_clear (u, n); mpres_clear (z, n); mpres_clear (b, n); return ret; } /* choose "optimal" S according to step 2 range B2 */ int choose_S (mpz_t B2len) { if (mpz_cmp_d (B2len, 1e7) < 0) return 1; /* x^1 */ else if (mpz_cmp_d (B2len, 1e8) < 0) return 2; /* x^2 */ else if (mpz_cmp_d (B2len, 1e9) < 0) return -3; /* Dickson(3) */ else if (mpz_cmp_d (B2len, 1e10) < 0) return -6; /* Dickson(6) */ else if (mpz_cmp_d (B2len, 3e11) < 0) return -12; /* Dickson(12) */ else return -30; /* Dickson(30) */ } static void print_expcurves (mpz_t B2min, mpz_t effB2, unsigned long dF, unsigned long k, int S, int clear) { double prob; int i; char sep; if (!test_verbose (OUTPUT_VERBOSE)) return; rhoinit (256, 10); outputf (OUTPUT_VERBOSE, "Expected number of curves to find a factor " "of n digits:\n20\t25\t30\t35\t40\t45\t50\t55\t60\t65\n"); for (i = 20; i <= 65; i += 5) { sep = (i < 65) ? '\t' : '\n'; prob = ecmprob (mpz_get_d (B2min), mpz_get_d (effB2), pow (10., i - .5), (double) dF * dF * k, S); if (prob > 1. / 10000000) outputf (OUTPUT_VERBOSE, "%.0f%c", floor (1. / prob + .5), sep); else if (prob > 0.) outputf (OUTPUT_VERBOSE, "%.2g%c", floor (1. / prob + .5), sep); else outputf (OUTPUT_VERBOSE, "Inf%c", sep); } if (clear) rhoinit (1, 0); /* Free memory of rhotable */ } static void print_exptime (mpz_t B2min, mpz_t effB2, unsigned long dF, unsigned long k, int S, double tottime, int clear) { double prob, exptime; int i; char sep; if (!test_verbose (OUTPUT_VERBOSE)) return; rhoinit (256, 10); outputf (OUTPUT_VERBOSE, "Expected time to find a factor of n digits:\n" "20\t25\t30\t35\t40\t45\t50\t55\t60\t65\n"); for (i = 20; i <= 65; i += 5) { sep = (i < 65) ? '\t' : '\n'; prob = ecmprob (mpz_get_d (B2min), mpz_get_d (effB2), pow (10., i - .5), (double) dF * dF * k, S); exptime = (prob > 0.) ? tottime / prob : HUGE_VAL; outputf (OUTPUT_TRACE, "Digits: %d, Total time: %.0f, probability: " "%g, expected time: %.0f\n", i, tottime, prob, exptime); if (exptime < 1000.) outputf (OUTPUT_VERBOSE, "%.0fms%c", exptime, sep); else if (exptime < 60000.) /* One minute */ outputf (OUTPUT_VERBOSE, "%.2fs%c", exptime / 1000., sep); else if (exptime < 3600000.) /* One hour */ outputf (OUTPUT_VERBOSE, "%.2fm%c", exptime / 60000., sep); else if (exptime < 86400000.) /* One day */ outputf (OUTPUT_VERBOSE, "%.2fh%c", exptime / 3600000., sep); else if (exptime < 31536000000.) /* One year */ outputf (OUTPUT_VERBOSE, "%.2fd%c", exptime / 86400000., sep); else if (exptime < 31536000000000.) /* One thousand years */ outputf (OUTPUT_VERBOSE, "%.2fy%c", exptime / 31536000000., sep); else if (exptime < 31536000000000000.) /* One million years */ outputf (OUTPUT_VERBOSE, "%.0fy%c", exptime / 31536000000., sep); else if (prob > 0.) outputf (OUTPUT_VERBOSE, "%.1gy%c", exptime / 31536000000., sep); else outputf (OUTPUT_VERBOSE, "Inf%c", sep); } if (clear) rhoinit (1, 0); /* Free memory of rhotable */ } /* Input: x is starting point or zero sigma is sigma value (if x is set to zero) or A parameter (if x is non-zero) of curve n is the number to factor go is the initial group order to preload B1, B2 are the stage 1/stage 2 bounds, respectively B2min the lower bound for stage 2 B2scale is the stage 2 scale factor k is the number of blocks to do in stage 2 S is the degree of the Suyama-Brent extension for stage 2 verbose is verbosity level: 0 no output, 1 normal output, 2 diagnostic output. sigma_is_a: If true, the sigma parameter contains the curve's A value Output: f is the factor found. Return value: ECM_FACTOR_FOUND_STEPn if a factor was found, ECM_NO_FACTOR_FOUND if no factor was found, ECM_ERROR in case of error. */ int ecm (mpz_t f, mpz_t x, mpz_t sigma, mpz_t n, mpz_t go, double *B1done, double B1, mpz_t B2min_parm, mpz_t B2_parm, double B2scale, unsigned long k, const int S, int verbose, int repr, int use_ntt, int sigma_is_A, FILE *os, FILE* es, char *TreeFilename, double maxmem, double stage1time, gmp_randstate_t rng, int (*stop_asap)(void)) { int youpi = ECM_NO_FACTOR_FOUND; int base2 = 0; /* If n is of form 2^n[+-]1, set base to [+-]n */ int Fermat = 0; /* If base2 > 0 is a power of 2, set Fermat to base2 */ int po2 = 0; /* Whether we should use power-of-2 poly degree */ long st; mpmod_t modulus; curve P; mpz_t B2min, B2; /* Local B2, B2min to avoid changing caller's values */ unsigned long dF; root_params_t root_params; set_verbose (verbose); ECM_STDOUT = (os == NULL) ? stdout : os; ECM_STDERR = (es == NULL) ? stdout : es; /* if n is even, return 2 */ if (mpz_divisible_2exp_p (n, 1)) { mpz_set_ui (f, 2); return ECM_FACTOR_FOUND_STEP1; } /* now n is odd */ /* check that B1 is not too large */ if (B1 > (double) ULONG_MAX) { outputf (OUTPUT_ERROR, "Error, maximal step 1 bound for ECM is %lu.\n", ULONG_MAX); return ECM_ERROR; } st = cputime (); /* See what kind of number we have as that may influence optimal parameter selection. Test for base 2 number */ if (repr != ECM_MOD_NOBASE2) base2 = (abs (repr) >= 16) ? repr : isbase2 (n, BASE2_THRESHOLD); /* For for a Fermat number (base2 a positive power of 2) */ for (Fermat = base2; Fermat > 0 && (Fermat & 1) == 0; Fermat >>= 1); if (Fermat == 1) { Fermat = base2; po2 = 1; } else Fermat = 0; if (base2) { if (mpmod_init_BASE2 (modulus, base2, n) == ECM_ERROR) return ECM_ERROR; } else if (repr == ECM_MOD_MPZ) mpmod_init_MPZ (modulus, n); else if (repr == ECM_MOD_MODMULN) mpmod_init_MODMULN (modulus, n); else if (repr == ECM_MOD_REDC) mpmod_init_REDC (modulus, n); else /* automatic choice of general reduction method */ mpmod_init (modulus, n, -1); MEMORY_TAG; mpres_init (P.x, modulus); MEMORY_TAG; mpres_init (P.y, modulus); MEMORY_TAG; mpres_init (P.A, modulus); mpres_set_z (P.x, x, modulus); mpres_set_ui (P.y, 1, modulus); MEMORY_TAG; mpz_init_set (B2min, B2min_parm); MEMORY_TAG; mpz_init_set (B2, B2_parm); MEMORY_TAG; mpz_init (root_params.i0); MEMORY_UNTAG; /* set second stage bound B2: when using polynomial multiplication of complexity n^alpha, stage 2 has complexity about B2^(alpha/2), and we want stage 2 to take about half of stage 1, thus we choose B2 = (c*B1)^(2/alpha). Experimentally, c=1/4 seems to work well. For Toom-Cook 3, this gives alpha=log(5)/log(3), and B2 ~ (c*B1)^1.365. For Toom-Cook 4, this gives alpha=log(7)/log(4), and B2 ~ (c*B1)^1.424. */ /* We take the cost of P+1 stage 1 to be about twice that of P-1. Since nai"ve P+1 and ECM cost respectively 2 and 11 multiplies per addition and duplicate, and both are optimized with PRAC, we can assume the ratio remains about 11/2. */ /* Also scale B2 by what the user said (or by the default scaling of 1.0) */ if (ECM_IS_DEFAULT_B2(B2)) mpz_set_d (B2, B2scale * pow (ECM_COST * B1, DEFAULT_B2_EXPONENT)); /* set B2min */ if (mpz_sgn (B2min) < 0) mpz_set_d (B2min, B1); /* Let bestD determine parameters for root generation and the effective B2 */ if (use_ntt) po2 = 1; if (bestD (&root_params, &k, &dF, B2min, B2, po2, use_ntt, maxmem, (TreeFilename != NULL), modulus) == ECM_ERROR) { youpi = ECM_ERROR; goto end_of_ecm; } /* Set default degree for Brent-Suyama extension */ /* We try to keep the time used by the Brent-Suyama extension at about 10% of the stage 2 time */ /* Degree S Dickson polys and x^S are equally fast for ECM, so we go for the better Dickson polys whenever possible. For S == 1, 2, they behave identically. */ root_params.S = S; if (root_params.S == ECM_DEFAULT_S) { if (Fermat > 0) { /* For Fermat numbers, default is 1 (no Brent-Suyama) */ root_params.S = 1; } else { mpz_t t; MEMORY_TAG; mpz_init (t); MEMORY_UNTAG; mpz_sub (t, B2, B2min); root_params.S = choose_S (t); mpz_clear (t); } } if (sigma_is_A == 0) { /* if sigma=0, generate it at random */ if (mpz_sgn (sigma) == 0) { mpz_urandomb (sigma, rng, 32); mpz_add_ui (sigma, sigma, 6); } /* sigma contains sigma value, A and x values must be computed */ /* If x is specified by user, keep that value */ youpi = get_curve_from_sigma (f, P.A, P.x, sigma, modulus); if (youpi != ECM_NO_FACTOR_FOUND) goto end_of_ecm; } else { /* sigma contains the A value */ mpres_set_z (P.A, sigma, modulus); /* TODO: make a valid, random starting point in case none was given */ /* Problem: this may be as hard as factoring as we'd need to determine whether x^3 + a*x^2 + x is a quadratic residue or not */ /* For now, we'll just chicken out. */ if (mpz_sgn (x) == 0) { outputf (OUTPUT_ERROR, "Error, -A requires a starting point (-x0 x).\n"); youpi = ECM_ERROR; goto end_of_ecm; } } /* If a nonzero value is given in x, then we use it as the starting point */ if (mpz_sgn (x) != 0) mpres_set_z (P.x, x, modulus); /* Now that the parameters are decided, print info */ outputf (OUTPUT_NORMAL, "Using B1=%1.0f, B2=", B1); if (mpz_cmp_d (B2min, B1) == 0) outputf (OUTPUT_NORMAL, "%Zd", B2); else outputf (OUTPUT_NORMAL, "%Zd-%Zd", B2min, B2); outputf (OUTPUT_NORMAL, ", polynomial "); if (root_params.S > 0) outputf (OUTPUT_NORMAL, "x^%u", root_params.S); else outputf (OUTPUT_NORMAL, "Dickson(%u)", -root_params.S); outputf (OUTPUT_NORMAL, ", sigma=%Zd\n", sigma); #if 0 outputf (OUTPUT_VERBOSE, "b2=%1.0f, dF=%lu, k=%lu, d=%lu, d2=%lu, i0=%Zd\n", b2, dF, k, root_params.d1, root_params.d2, root_params.i0); #else outputf (OUTPUT_VERBOSE, "dF=%lu, k=%lu, d=%lu, d2=%lu, i0=%Zd\n", dF, k, root_params.d1, root_params.d2, root_params.i0); #endif if (test_verbose (OUTPUT_RESVERBOSE)) { mpz_t t; MEMORY_TAG; mpz_init (t); MEMORY_UNTAG; mpres_get_z (t, P.A, modulus); outputf (OUTPUT_RESVERBOSE, "a=%Zd\n", t); mpres_get_z (t, P.x, modulus); outputf (OUTPUT_RESVERBOSE, "starting point: x=%Zd\n", t); mpz_clear (t); } if (go != NULL && mpz_cmp_ui (go, 1) > 0) outputf (OUTPUT_VERBOSE, "initial group order: %Zd\n", go); print_expcurves (B2min, B2, dF, k, root_params.S, 0); #ifdef HAVE_GWNUM /* Right now, we only do base 2 numbers with GWNUM */ if (base2 != 0 && B1 >= *B1done) youpi = gw_ecm_stage1 (f, &P, modulus, B1, B1done, go); /* At this point B1 == *B1done unless interrupted, or no GWNUM ecm_stage1 is available */ if (youpi != ECM_NO_FACTOR_FOUND) goto end_of_ecm; #endif if (B1 > *B1done) youpi = ecm_stage1 (f, P.x, P.A, modulus, B1, B1done, go, stop_asap); if (stage1time > 0.) { const long st2 = elltime (st, cputime ()); const long s1t = (long) (stage1time * 1000.); outputf (OUTPUT_NORMAL, "Step 1 took %ldms (%ld in this run, %ld from previous runs)\n", st2 + s1t, st2, s1t); } else outputf (OUTPUT_NORMAL, "Step 1 took %ldms\n", elltime (st, cputime ())); /* Store end-of-stage-1 residue in x in case we write it to a save file, before P.x is converted to Weierstrass form */ mpres_get_z (x, P.x, modulus); if (youpi != ECM_NO_FACTOR_FOUND) goto end_of_ecm; if (test_verbose (OUTPUT_RESVERBOSE)) { mpz_t t; MEMORY_TAG; mpz_init (t); MEMORY_UNTAG; mpres_get_z (t, P.x, modulus); outputf (OUTPUT_RESVERBOSE, "x=%Zd\n", t); mpz_clear (t); } /* In case of a signal, we'll exit after the residue is printed. If no save file is specified, the user may still resume from the residue */ if (stop_asap != NULL && (*stop_asap) ()) goto end_of_ecm; #ifdef MONT_ROOTS /* If we want to use Montgomery form for generating the roots for stage 2, convert to Weierstrass form only if S != 1 */ if (root_params.S != 1) { mpz_t t; outputf (OUTPUT_DEVVERBOSE, "ecm: Converting to Weierstrass form\n"); youpi = montgomery_to_weierstrass (f, P.x, P.y, P.A, modulus); if (test_verbose (OUTPUT_TRACE)) { MEMORY_TAG; mpz_init (t); MEMORY_UNTAG; mpres_get_z (t, P.x, modulus); outputf (OUTPUT_TRACE, "P = (%Zd, ", t); mpres_get_z (t, P.y, modulus); outputf (OUTPUT_TRACE, "%Zd)\n", t); mpz_clear (t); } } #else youpi = montgomery_to_weierstrass (f, P.x, P.y, P.A, modulus); #endif if (youpi == ECM_NO_FACTOR_FOUND && mpz_cmp (B2, B2min) >= 0) youpi = stage2 (f, &P, modulus, dF, k, &root_params, ECM_ECM, use_ntt, TreeFilename, stop_asap); if (youpi == ECM_NO_FACTOR_FOUND && (stop_asap == NULL || !(*stop_asap)())) print_exptime (B2min, B2, dF, k, root_params.S, (long) (stage1time * 1000.) + elltime (st, cputime ()), 1); end_of_ecm: mpres_clear (P.A, modulus); mpres_clear (P.y, modulus); mpres_clear (P.x, modulus); mpmod_clear (modulus); mpz_clear (root_params.i0); mpz_clear (B2); mpz_clear (B2min); return youpi; }