/* * Copyright (c) 1993 David I. Bell * Permission is granted to use, distribute, or modify this source, * provided that this copyright notice remains intact. * * Attempt to factor numbers using elliptic functions. * y^2 = x^3 + a*x + b (mod N). * * Many points (x,y) (mod N) are found that solve the above equation, * starting from a trivial solution and 'multiplying' that point together * to generate high powers of the point, looking for such a point whose * order contains a common factor with N. The order of the group of points * varies almost randomly within a certain interval for each choice of a * and b, and thus each choice provides an independent opportunity to * factor N. To generate a trivial solution, a is chosen and then b is * selected so that (1,1) is a solution. The multiplication is done using * the basic fact that the equation is a cubic, and so if a line hits the * curve in two rational points, then the third intersection point must * also be rational. Thus by drawing lines between known rational points * the number of rational solutions can be made very large. When modular * arithmetic is used, solving for the third point requires the taking of a * modular inverse (instead of division), and if this fails, then the GCD * of the failing value and N provides a factor of N. This description is * only an approximation, read "A Course in Number Theory and Cryptography" * by Neal Koblitz for a good explanation. * * factor(iN, ia, B, force) * iN is the number to be factored. * ia is the initial value of a in the equation, and each successive * value of a is an independent attempt at factoring (default 1). * B is the limit of the primes that make up the high power that the * point is raised to for each factoring attempt (default 100). * force is a flag to attempt to factor numbers even if they are * thought to already be prime (default FALSE). * * Making B larger makes the power the point being raised to contain more * prime factors, thus increasing the chance that the order of the point * will be made up of those factors. The higher B is then, the greater * the chance that any individual attempt will find a factor. However, * a higher B also slows down the number of independent functions being * examined. The order of the point for any particular function might * contain a large prime and so won't succeed even for a really large B, * whereas the next function might have an order which is quickly found. * So you want to trade off the depth of a particular search with the * number of searches made. For example, for factoring 30 digits, I make * B be about 1000 (probably still too small). * * If you have lots of machines available, then you can run parallel * factoring attempts for the same number by giving different starting * values of ia for each machine (e.g. 1000, 2000, 3000). * * The output as the function is running is (occasionally) the value of a * when a new function is started, the prime that is being included in the * high power being calculated, and the current point which is the result * of the powers so far. * * If a factor is found, it is returned and is also saved in the global * variable f. The number being factored is also saved in the global * variable N. */ obj point {x, y}; global N; /* number to factor */ global a; /* first coefficient */ global b; /* second coefficient */ global f; /* found factor */ define factor(iN, ia, B, force) { local C, x, p; if (!force && ptest(iN, 50)) return 1; if (isnull(B)) B = 100; if (isnull(ia)) ia = 1; obj point x; a = ia; b = -ia; N = iN; C = isqrt(N); C = 2 * C + 2 * isqrt(C) + 1; f = 0; while (f == 0) { print "A =", a; x.x = 1; x.y = 1; print 2, x; x = x ^ (2 ^ (highbit(C) + 1)); for (p = 3; ((p < B) && (f == 0)); p += 2) { if (!ptest(p, 1)) continue; print p, x; x = x ^ (p ^ ((highbit(C) // highbit(p)) + 1)); } a++; b--; } return f; } define point_print(p) { print "(" : p.x : "," : p.y : ")" :; } define point_mul(p1, p2) { local r, m; if (p2 == 1) return p1; if (p1 == p2) return point_square(&p1); obj point r; m = (minv(p2.x - p1.x, N) * (p2.y - p1.y)) % N; if (m == 0) { if (f == 0) f = gcd(p2.x - p1.x, N); r.x = 1; r.y = 1; return r; } r.x = (m^2 - p1.x - p2.x) % N; r.y = ((m * (p1.x - r.x)) - p1.y) % N; return r; } define point_square(p) { local r, m; obj point r; m = ((3 * p.x^2 + a) * minv(p.y << 1, N)) % N; if (m == 0) { if (f == 0) f = gcd(p.y << 1, N); r.x = 1; r.y = 1; return r; } r.x = (m^2 - p.x - p.x) % N; r.y = ((m * (p.x - r.x)) - p.y) % N; return r; } define point_pow(p, pow) { local bit, r, t; r = 1; if (isodd(pow)) r = p; t = p; for (bit = 2; ((bit <= pow) && (f == 0)); bit <<= 1) { t = point_square(&t); if (bit & pow) r = point_mul(&t, &r); } return r; } global lib_debug; if (!isnum(lib_debug) || lib_debug>0) print "factor(N, I, B, force) defined";