calc/cal/ellip.cal

/*
 * ellip - attempt to factor numbers using elliptic functions
 *
 * Copyright (C) 1999  David I. Bell
 *
 * Calc is open software; you can redistribute it and/or modify it under
 * the terms of the version 2.1 of the GNU Lesser General Public License
 * as published by the Free Software Foundation.
 *
 * Calc 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.
 *
 * A copy of version 2.1 of the GNU Lesser General Public License is
 * distributed with calc under the filename COPYING-LGPL.  You should have
 * received a copy with calc; if not, write to Free Software Foundation, Inc.
 * 59 Temple Place, Suite 330, Boston, MA  02111-1307, USA.
 *
 * @(#) $Revision: 29.4 $
 * @(#) $Id: ellip.cal,v 29.4 2006/06/20 09:29:16 chongo Exp $
 * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/ellip.cal,v $
 *
 * Under source code control:	1990/02/15 01:50:33
 * File existed as early as:	before 1990
 *
 * Share and enjoy!  :-)	http://www.isthe.com/chongo/tech/comp/calc/
 */

/*
 * Attempt to factor numbers using elliptic functions:
 *
 *	y^2 = x^3 + a*x + b   (mod ellip_N).
 *
 * Many points (x,y) (mod ellip_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 ellip_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 ellip_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 ellip_N provides a factor of ellip_N.
 * This description is only an approximation, read "A Course in Number
 * Theory and Cryptography" by Neal Koblitz for a good explanation.
 *
 * efactor(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 ellip_N.
 */


obj point {x, y};
global	ellip_N;		/* number to factor */
global	ellip_a;		/* first coefficient */
global	ellip_b;		/* second coefficient */
global	ellip_f;		/* found factor */


define efactor(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;
	ellip_a = ia;
	ellip_b = -ia;
	ellip_N = iN;
	C = isqrt(ellip_N);
	C = 2 * C + 2 * isqrt(C) + 1;
	ellip_f = 0;
	while (ellip_f == 0) {
		print "A =", ellip_a;
		x.x = 1;
		x.y = 1;
		print 2, x;
		x = x ^ (2 ^ (highbit(C) + 1));
		for (p = 3; ((p < B) && (ellip_f == 0)); p += 2) {
			if (!ptest(p, 1))
				continue;
			print p, x;
			x = x ^ (p ^ ((highbit(C) // highbit(p)) + 1));
		}
		ellip_a++;
		ellip_b--;
	}
	return ellip_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, ellip_N) * (p2.y - p1.y)) % ellip_N;
	if (m == 0) {
		if (ellip_f == 0)
			ellip_f = gcd(p2.x - p1.x, ellip_N);
		r.x = 1;
		r.y = 1;
		return r;
	}
	r.x = (m^2 - p1.x - p2.x) % ellip_N;
	r.y = ((m * (p1.x - r.x)) - p1.y) % ellip_N;
	return r;
}


define point_square(p)
{
	local	r, m;

	obj point r;
	m = ((3 * p.x^2 + ellip_a) * minv(p.y << 1, ellip_N)) % ellip_N;
	if (m == 0) {
		if (ellip_f == 0)
			ellip_f = gcd(p.y << 1, ellip_N);
		r.x = 1;
		r.y = 1;
		return r;
	}
	r.x = (m^2 - p.x - p.x) % ellip_N;
	r.y = ((m * (p.x - r.x)) - p.y) % ellip_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) && (ellip_f == 0)); bit <<= 1) {
		t = point_square(`t);
		if (bit & pow)
			r = point_mul(`t, `r);
	}
	return r;
}