Release calc version 2.12.0.4

cal/ellip.cal

@@ -17,8 +17,8 @@
  * 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.3 $
- * @(#) $Id: ellip.cal,v 29.3 2006/03/07 22:16:25 chongo Exp $
+ * @(#) $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
@@ -28,16 +28,17 @@
  */
 
 /*
- * Attempt to factor numbers using elliptic functions.
- *	y^2 = x^3 + a*x + b   (mod N).
+ * Attempt to factor numbers using elliptic functions:
  *
- * Many points (x,y) (mod N) are found that solve the above equation,
+ *	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 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
+ * 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
@@ -45,9 +46,9 @@
  * 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.
+ * 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.
@@ -81,15 +82,15 @@
  *
  * 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.
+ * variable ellip_N.
  */
 
 
 obj point {x, y};
-global	N;		/* number to factor */
-global	a;		/* first coefficient */
-global	b;		/* second coefficient */
-global	f;		/* found factor */
+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)
@@ -103,28 +104,28 @@ define efactor(iN, ia, B, force)
 	if (isnull(ia))
 		ia = 1;
 	obj point x;
-	a = ia;
-	b = -ia;
-	N = iN;
-	C = isqrt(N);
+	ellip_a = ia;
+	ellip_b = -ia;
+	ellip_N = iN;
+	C = isqrt(ellip_N);
 	C = 2 * C + 2 * isqrt(C) + 1;
-	f = 0;
-	while (f == 0) {
-		print "A =", a;
+	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) && (f == 0)); p += 2) {
+		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));
 		}
-		a++;
-		b--;
+		ellip_a++;
+		ellip_b--;
 	}
-	return f;
+	return ellip_f;
 }
 
 
@@ -143,16 +144,16 @@ define point_mul(p1, p2)
 	if (p1 == p2)
 		return point_square(`p1);
 	obj point r;
-	m = (minv(p2.x - p1.x, N) * (p2.y - p1.y)) % N;
+	m = (minv(p2.x - p1.x, ellip_N) * (p2.y - p1.y)) % ellip_N;
 	if (m == 0) {
-		if (f == 0)
-			f = gcd(p2.x - p1.x, N);
+		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) % N;
-	r.y = ((m * (p1.x - r.x)) - p1.y) % N;
+	r.x = (m^2 - p1.x - p2.x) % ellip_N;
+	r.y = ((m * (p1.x - r.x)) - p1.y) % ellip_N;
 	return r;
 }
 
@@ -162,16 +163,16 @@ define point_square(p)
 	local	r, m;
 
 	obj point r;
-	m = ((3 * p.x^2 + a) * minv(p.y << 1, N)) % N;
+	m = ((3 * p.x^2 + ellip_a) * minv(p.y << 1, ellip_N)) % ellip_N;
 	if (m == 0) {
-		if (f == 0)
-			f = gcd(p.y << 1, N);
+		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) % N;
-	r.y = ((m * (p.x - r.x)) - p.y) % N;
+	r.x = (m^2 - p.x - p.x) % ellip_N;
+	r.y = ((m * (p.x - r.x)) - p.y) % ellip_N;
 	return r;
 }
 
@@ -184,7 +185,7 @@ define point_pow(p, pow)
 	if (isodd(pow))
 		r = p;
 	t = p;
-	for (bit = 2; ((bit <= pow) && (f == 0)); bit <<= 1) {
+	for (bit = 2; ((bit <= pow) && (ellip_f == 0)); bit <<= 1) {
 		t = point_square(`t);
 		if (bit & pow)
 			r = point_mul(`t, `r);