Lenstra’s Elliptic Curve Method of Factoring Large Integers Using the Free Large Integer Package (FreeLIP) by James Pate Williams, Jr.

/*
	Author:  Pate Williams (c) 1997 - 2009
	
	Algorithm 10.3.3 (Lenstra's ECM). See "A Course
	in Computational Algebraic Number Theory" by
	Henri Cohen page 488.
*/

#include <stdlib.h>
#include <string.h>
#include <time.h>
#include "lip.h"

#ifndef CLK_TCK
#define CLK_TCK CLOCKS_PER_SEC
#endif

/*
	if you creating a Windows only program
	the following function can be edited
	out or completely removed since strrev
	is defined in Windows Visual Studio C
*/

char * strrev(char *s)
{
	int i, length = strlen(s), j = length - 1;
	char *t = malloc((length + 1) * sizeof(char));
	
	if (t == NULL)
		return NULL;

	for (i = 0; i < length; i++, j--)
		t[j] = s[i];
		
	t[length] = '\0';
	
	strcpy(s, t);
	free(t);
	
	return s;
}

#define CURVES 1024l

struct point {verylong zx, zy, zz;};
struct factor {int expon; verylong prime;};

static char cexp[32][32], cfac[32][512];
static int cnt, primeCount;

int partial_addition(verylong za,
					 verylong zn,
					 struct point P,
					 struct point Q,
					 struct point *R,
					 verylong *zd)

/* returns 0 if sum is found or 1 if divisor is found */

{
	int value = 0;
	verylong zb = 0, zc = 0, zl = 0, zs = 0, zt = 0;
	verylong zx = 0, zy = 0, zy2 = 0;
	
	if (zcompare(P.zx, Q.zx) == 0 &&
		zscompare(P.zy, 0l) == 0 &&
		zscompare(Q.zy, 0l) == 0) {
		zzero(&R->zx);
		zone(&R->zy);
		return 0;
	}

	zsub(zn, Q.zy, &zb);
	
	if (zcompare(P.zx, Q.zx) == 0 &&
		zcompare(P.zy, zb) == 0) {
		zzero(&R->zx);
		zone(&R->zy);
		zfree(&zb);
		return 0;
	}
	
	if (zscompare(P.zx, 0l) == 0 &&
		zscompare(P.zy, 1l) == 0 &&
		zscompare(P.zz, 0l) == 0) {
		
		/* O + Q = Q */
		
		zcopy(Q.zx, &R->zx);
		zcopy(Q.zy, &R->zy);
		zcopy(Q.zz, &R->zz);
		value = 0;
	}

	else if (zscompare(Q.zx, 0l) == 0 &&
		zscompare(Q.zy, 1l) == 0 &&
		zscompare(Q.zz, 0l) == 0) {

		/* P + O = P */
		
		zcopy(P.zx, &R->zx);
		zcopy(P.zy, &R->zy);
		zcopy(P.zz, &R->zz);
		value = 0;
	}

	else {
		
		/* P != O and Q != O */
		
		zcopy(Q.zy, &zy2);
		znegate(&zy2);
		if (zcompare(P.zx, Q.zx) == 0  &&
			zcompare(P.zy, zy2)  == 0) {
			zzero(&R->zx);
			zone(&R->zy);
			zzero(&R->zz);
		}

		else {
			if (zcompare(P.zx, Q.zx) != 0) {
				zsubmod(P.zx, Q.zx, zn, &zx);
				zexteucl(zx, &zs, zn, &zt, zd);
				if (zscompare(*zd, 1l) != 0) goto L1;
				zsubmod(P.zy, Q.zy, zn, &zy);
				zmulmod(zs, zy, zn, &zl);
			}

			else {
				zaddmod(P.zy, Q.zy, zn, &zy);
				zexteucl(zy, &zs, zn, &zt, zd);
				if (zscompare(*zd, 1l) != 0) goto L1;
				zmulmod(P.zx, P.zx, zn, &zb);
				zsmulmod(zb, 3l, zn, &zc);
				zaddmod(zc, za, zn, &zb);
				zmulmod(zs, zb, zn, &zl);
			}

			zmulmod(zl, zl, zn, &zb);
			zaddmod(P.zx, Q.zx, zn, &zc);
			zsubmod(zb, zc, zn, &zx);
			zcopy(zx, &R->zx);
			zsubmod(zx, P.zx, zn, &zb);
			zmulmod(zl, zb, zn, &zc);
			zaddmod(zc, P.zy, zn, &zy);
			znegate(&zy);
			zcopy(zy, &R->zy);
			zone(&R->zz);
			goto L2;
L1:
			value = 1;
L2:
			;
		}
	}
	zfree(&zb);
	zfree(&zc);
	zfree(&zl);
	zfree(&zs);
	zfree(&zt);
	zfree(&zx);
	zfree(&zy);
	zfree(&zy2);
	return value;
}

int multiply(long k,
			 verylong za,
			 verylong zn,
			 struct point P,
			 struct point *R,
			 verylong *zd)
{
	int value = 0;
	struct point A, S, T;
	
	A.zx = A.zy = A.zz = S.
	  zx = S.zy = S.zz = 0;
	T.zx = T.zy = T.zz = 0;

	zzero(&R->zx);
	zone(&R->zy);
	zzero(&R->zz);
	zcopy(P.zx, &S.zx);
	zcopy(P.zy, &S.zy);
	zcopy(P.zz, &S.zz);
	
	while (!value && k != 0) {
		if (k & 1) {
			value = partial_addition(za, zn, *R, S, &A, zd);
			zcopy(A.zx, &R->zx);
			zcopy(A.zy, &R->zy);
			zcopy(A.zz, &R->zz);
		}

		k >>= 1;
		
		if (!value && k != 0) {
			value = partial_addition(za, zn, S, S, &T, zd);
			zcopy(T.zx, &S.zx);
			zcopy(T.zy, &S.zy);
			zcopy(T.zz, &S.zz);
		}
	}

	if (zscompare(R->zy, 0l) < 0) {
		zadd(R->zy, zn, &A.zy);
		zcopy(A.zy, &R->zy);
	}

	zfree(&A.zx);
	zfree(&A.zy);
	zfree(&A.zz);
	zfree(&S.zx);
	zfree(&S.zy);
	zfree(&S.zz);
	zfree(&T.zx);
	zfree(&T.zy);
	zfree(&T.zz);
	return value;
}

/* the following definition limits the L3 loop */

int LenstrasECM(verylong *zN, verylong *zg)
{
	int expon = 0, found;
	long B = 2000000l, i, j, k = 0, l, *p, q, q1;
	long giveUp = z2log(*zN);
	struct point x, y;
	verylong za[CURVES], zb = 0, zd = 0;
	
	for (i = 0; i < CURVES; i++)
		za[i] = 0;
	
	x.zx = x.zy = x.zz = y.zx = y.zy = y.zz = 0;
	
	zpstart2();

	do {
		q = zpnext();
		k++;
	} while (q < B);
	
	p = calloc(k, sizeof(long));
	
	if (!p) {
		expon = - 1;
		goto L4;
	}
	
	zpstart2();
	
	for (i = 0; i < k; i++)
		p[i] = zpnext();
	
	for (i = 0; i < CURVES; i++)
		zrandomb(*zN, &za[i]);

L2:

	zone(&x.zx);
	zone(&x.zy);
	zzero(&x.zz);
	i = - 1;

L3:
	
	i++;
	
	if (i == giveUp) {
		expon = 0;
		goto L4;
	}
	
	if (i == k) {
	
		for (i = 0; i < CURVES; i++)
			zrandomb(*zN, &za[i]);

		goto L2;
	}

	q = p[i];
	q1 = q;
	l = B / q;
	
	while (q1 <= l)
		q1 *= q;
	
	found = 0;
	
	for (j = 0; !found && j < CURVES; j++)
		found = multiply(q1, za[j], *zN, x, &y, &zd);
	
	if (!found)
		goto L3;
	
	zcopy(y.zx, &x.zx);
	zcopy(y.zy, &x.zy);
	zcopy(y.zz, &x.zz);
	zgcd(zd, *zN, zg);

	if (zcompare(*zg, *zN) == 0) {
		for (j = 0; j < CURVES; j++)
			zrandomb(*zN, &za[j]);
	
		goto L2;
	}

	if (!zprobprime(*zg, 5l))
		goto L4;

	expon = 0;

	do {
		zdiv(*zN, *zg, &zb, &zd);
		zcopy(zb, zN);
		zmod(*zN, *zg, &zd);
		expon++;
	} while (zscompare(zd, 0l) == 0);

L4:
	
	for (i = 0; i < CURVES; i++)
		zfree(&za[i]);
	
	zfree(&zb);
	zfree(&zd);
	zfree(&x.zx);
	zfree(&x.zy);
	zfree(&x.zz);
	zfree(&y.zx);
	zfree(&y.zy);
	zfree(&y.zz);
	return expon;
}

int ECMVLConvert(verylong zx, char *str, int length)
{
	verylong zy = 0;
	int a, result = 1, i = 0;
	
	while (zscompare(zx, 0) > 0)
	{
		zsdiv(zx, 10, &zy);
		a = zsmod(zx, 10);
		zcopy(zy, &zx);
		
		if (i < length)
			str[i++] = a + '0';
		
		else
		{
			result = 0;
			break;
		}
	}
	
	str[i] = '\0';
	strrev(str);
	
	zfree(&zy);
	return result;
}

int ECMLast(char exp[32][32], char fac[32][512], struct factor *f, int *cnt)
{
	int i;
	
	for (i = 0; i < *cnt; i++)
	{
		if (ECMVLConvert(f[i].prime, fac[i], 512) == 1)
			sprintf(exp[i], "%d", f[i].expon);
		
		else
			return 0;
	}
	
	return 1;
}

int DoECM(char exp[32][32], char fac[32][512], char *numStr,
		  int *cnt, int base, int iExpon, int addend)
{
	struct factor f[32];
	int i, j, result;
	int cant, expon, one, pri;
	verylong zn = 0, zN = 0, zd = 0, zg = 0;
	
	if (strlen(numStr) != 0)
		zstrtoz(numStr, &zN);
	
	else
	{
		static verylong zb = 0, zp = 0;
		
		zintoz(base, &zb);
		
		zsexp(zb, iExpon, &zp);
		zsadd(zp, addend, &zN);
	}
	
	zcopy(zN, &zn);

	for (i = 0; i < 32; i++)
	{
		f[i].expon = 0;
		f[i].prime = 0;
	}
	
	*cnt = 0;
	
	f[*cnt].expon = 1;
	zcopy(zN, &f[*cnt].prime);
	*cnt = *cnt + 1;
	
	if (z2log(zn) > 512)
	{
		result = 3;
		goto L1;
	}
	
	if (zprobprime(zN, 5l))
	{
		result = 2;
		goto L1;
	}
	
	cant = pri = 0;
	
	do {

		expon = LenstrasECM(&zN, &zg);

		if (expon == - 1)
		{
			result = - 1;
			goto L1;
		}

		if (zprobprime(zg, 5l))
		{
			f[*cnt].expon = expon;
			zcopy(zg, &f[*cnt].prime);
			*cnt = *cnt + 1;
		}

		one = zscompare(zN, 1l) == 0;

		if (zprobprime(zN, 5l))
		{
			f[*cnt].expon = 1;
			zcopy(zN, &f[*cnt].prime);
			*cnt = *cnt + 1;
			break;
		}

	} while (!one);
	
	/* selection sort the factors */
	
	for (i = 1; i < *cnt - 1; i++)
	{
		for (j = i + 1; j < *cnt; j++)
		{
			if (zcompare(f[i].prime, f[j].prime) > 0)
			{
				int temp;
				verylong zt = 0;
					
				zcopy(f[i].prime, &zt);
				zcopy(f[j].prime, &f[i].prime);
				zcopy(zt, &f[j].prime);
				
				temp = f[i].expon;
				f[i].expon = f[j].expon;
				f[j].expon = temp;
			}
		}
	}

	result = ECMLast(exp, fac, f, cnt);

L1:

	if (result == 2)
	{
		int temp = ECMLast(exp, fac, f, cnt);
		
		if (temp == 0)
			result = 3;
	}
	
	zfree(&zn);
	zfree(&zN);
	zfree(&zd);
	zfree(&zg);

	return result;
}

int PrimeCheck(char fac[32][512], int cnt)
{
	int i, result = 0;
	verylong zp = 0;
	
	for (i = 0; i < cnt; i++)
	{	
		zstrtoz(fac[i], &zp);

		if (zprobprime(zp, 5l))
			result++;
	}

	return result;
}

int ECMMethod(char *numStr, int base, int expon, int addend)
{
	int i, result;
	
	cnt = 0;
	
	for (i = 0; i < 32; i++)
	{
		cexp[i][0] = '\0';
		cfac[i][0] = '\0';
	}
	
	result = DoECM(cexp, cfac, numStr, &cnt, base, expon, addend);
	
	primeCount = PrimeCheck(cfac, cnt);

	return result;
}

int main(int argc, char *argv[])
{
	while (1) {
		clock_t time0;
		char numStr[1024] = {0}, *bPtr, *ePtr;
		double time;
		int c, base, expon, addend, i = 0, positive, result;

		printf("enter the number below (0 to quit):\n");
		
		c = getchar();
		
		while (c != '\n' && i < 1023)
		{
			numStr[i++] = c;
			c = getchar();
		}
		
		numStr[i] = '\0';
		
		bPtr = strchr(numStr, '^');

		if (bPtr != NULL)
		{
			*bPtr = '\0';
			base = atoi(numStr);
			ePtr = strchr(bPtr + 1, '+');
			
			if (ePtr != NULL)
				positive = 1;

			else {
				ePtr = strchr(bPtr + 1, '-');
				
				if (ePtr != NULL)
					positive = 0;
			}

			if (ePtr != NULL)
				*ePtr = '\0';
			
			expon = atoi(bPtr + 1);
			
			if (ePtr != NULL) {
				addend = atoi(ePtr + 1);

				if (positive == 0)
					addend = - addend;
			}
			
			else
				addend = 0;

			numStr[0] = '\0';
		}
		
		else
		{
			addend = atoi(numStr);

			if (addend == 0)
				break;
		}

		printf("\n");

		time0 = clock();

		result = ECMMethod(numStr, base, expon, addend);

		if (result == 0)
			printf("can't completely factor this number\n");
		
		else if (result == 1)
		{
			printf("%s\nnumber is composite factors:\n", cfac[0]);
			
			for (i = 1; i < cnt; i++)
				if (strcmp(cexp[i], "1") == 0)
					printf("\t%s\n", cfac[i]);
				else
					printf("\t%s ^ %s\n", cfac[i], cexp[i]);
		}
		
		else if (result == 2)
			printf("%s\nnumber is prime:\n", cfac[0]);

		else if (result == 3)
			printf("number or prime factor is too large\n");

		time = (clock() - time0) / (double) CLK_TCK;
		printf("total time required: %f seconds\n", time);
		printf("\n");
	}

	return 0;
}