/*
        ecm.c
        Kimball Martin
        CMSC 443
        5/17/97

        An interactive C implementation of H. Lenstra's elliptic curve factoring 
	method using A. Lenstra's LIP package.  The interaction allows the user
	to choose the curve (restricted to Weierstrass normal form) with which 
	to attempt factoring.

        To compile:
		make lip.o
		cc  ecm.c lip.o -o ecm -lm
        Usage:
                ecm 
        and the program will prompt the user for the data.
*/

#include "lip.h"

int ecm(verylong, verylong, verylong, verylong, verylong);
int modmult(verylong, verylong, verylong, verylong, verylong, verylong);
int padd(verylong, verylong, verylong*, verylong*, verylong, verylong);
int pdouble(verylong*, verylong*, verylong, verylong);

/* ecm(n, B, X, Y, a) does one iteration of the elliptic curve method.
   Parameters: n, the number to factor,
               B, the upper bound for the factor base,
               X, the x-value of the starting point,
               Y, the y-value of the starting point,
               a, the coefficient on the x term of the curve in Weierstrass normal
                   form:        y^3 = x^3 + ax + b
   Notes:  k = B!  Some speed-up could be obtained if we change k to
		lcm[1, 2, ..., B] for large n and large B.  

   modmult(X, Y, k, a, b, n) adds the point (X, Y) k times to itself on the
      elliptic curve defined by a and b (using Weierstrass form above) over Z_n.
   Parameters: X, the x-value of the starting point,
               Y, the y-value of the starting point,
               k, the scalar multiple for (X, Y),
               a, the coefficient on the x term of the curve in Weierstrass form,
               b, the constant term of the curve in Weierstrass form,
               n, the modulus for the curve.

   padd(X2, Y2, &X3, &Y3, a, n) adds the point (X2, Y2) to the point (X3, Y3) on
      the curve with linear coefficient a over Z_n, so (X3, Y3) = (X2, Y2) + (X3, Y3)
   Parameters: X2, the x-value of the first point
	       Y2, the y-value of the first point
	       X3, the x-value of the second point
	       Y3, the y-value of the second point
	       a, the coefficient on the x term of the curve in Weierstrass form,
	       n, the modulus for the curve.
   Notes:  X2, Y2, X3 and Y3 are all expected to be different integers, i.e. they
		may have the same value, but not point to the same memory location.
		To add a point to itself, use pdouble (see below).

   pdouble(&X, &Y, a, n) sets (X, Y) = (X, Y) + (X, Y) over the curve in Z_N with 
      linear coefficient a.
   Parameters: X, the x-value of the point
	       Y, the y-value of the point
	       a, the coefficient on the x term of the curve in Weierstrass form,
	       n, the modulus for the curve.

*/


void main() {
  char c;
  verylong n = 0, B = 0, X = 0, Y = 0, a = 0;

  printf("Enter integer n to factor:  ");	/* read in data */
  zread(&n);
  do {
    fflush(stdin);
    printf("Enter factoring bound B:  ");
    zread(&B);
    fflush(stdin);
    printf("Enter the x-value of the starting point : ");
    zread(&X);
    fflush(stdin);
    printf("Enter the y-value of the starting point : ");
    zread(&Y);
    fflush(stdin);
    printf("Enter the linear coefficient for the curve: ");
    zread(&a);

    if (ecm(n, B, X, Y, a) == 0) printf("ECM failed.\n");
    fflush(stdin);				/* try to factor */
    printf("Continue [y/n]?");
    scanf("%c", &c);
  } while (c != 'n');
}

int ecm(verylong n, verylong B, verylong X,verylong Y,verylong a) {
  verylong b = 0, g = 0, k = 0, tmp = 0, tmp2 = 0;

  zsqmod(Y, n, &b);				/* Let b = Y^2 - X^3 - aX (mod n) */
  zsexpmod(X, 3, n, &tmp);
  zsubmod(b, tmp, n, &tmp2);
  zmulmod(a, X, n, &tmp);
  zsubmod(tmp2, tmp, n, &b);

  zsexpmod(a, 3, n, &g);			/* Let g = gcd(4a^3 + 27b^2, n) */
  zsmulmod(g, 4, n, &g);
  zsqmod(b, n, &tmp);
  zsmulmod(tmp, 27, n, &tmp);
  zaddmod(g, tmp, n, &g);
  zgcd(g, n, &g);

  if (zscompare(g, 1) != 0) {			/* If g <> 1 then our curve is not */
    printf("Not an elliptic curve.\n");		/* an elliptic curve. */
    if (zcompare(g, n) != 0) {			/* However if g < n then we have a */
      printf("n factors into:\n");		/* non-trivial factor of n. */
      zwriteln(g);
      zdiv(n, g, &g, &tmp);
      zwriteln(g);
      return 1;
    }
    return 0;
  }

						/* Let k = B! */
  for(zone(&k), zone(&tmp); zcompare(tmp, B) != 1; zsadd(tmp, 1, &tmp))  {
    zmul(k, tmp, &tmp2);
    zcopy(tmp2, &k);
  }
  return modmult(X, Y, k, a, b, n);

}

int modmult(verylong X, verylong Y, verylong k, verylong a, verylong b, verylong n) {
  int f = 0;
  verylong e = 0, X2 = 0, X3 = 0, Y2 = 0, Y3 = 0, tmp = 0;

  zcopy(X, &X3);
  zcopy(Y, &Y3);

  while(zscompare(k, 1) > 0) {			/* compute (X3, Y3) = k(X,Y) mod n */
    zintoz(2, &e);				/* by the repeated doubling method */
    zcopy(X, &X2);
    zcopy(Y, &Y2);
    for(;;) {
      if (pdouble(&X2, &Y2, a, n) == 0) return 1;
      z2mul(e, &tmp);
      if (zcompare(tmp, k) < 0) zcopy(tmp, &e);
      else break;
    }
    if (f == 0) {
      zcopy(X2, &X3);
      zcopy(Y2, &Y3);
      f = 1;
    }
    else
      if (padd(X2, Y2, &X3, &Y3, a, n) == 0) return 1;
    zsub(k, e, &k);
  }

  return 0;

}

int padd(verylong X2, verylong Y2, verylong *X3, verylong *Y3, verylong a, verylong n) {
  verylong m = 0, tmp = 0, tmp2 = 0;

  if (zcompare(X2, *X3) == 0) {			/* If X2 = X3, */
    zaddmod(Y2, *Y3, n, &tmp);			/* Let m = (3 * X2^2 + a)/(Y2 + Y3) */ 

    if (ziszero(tmp)) {				/* If tmp = 0 */
      printf("Addition failed with no non-trivial divisors.\n");
      return 0;					/* exit with failure */
    }
    zinv(tmp, n, &tmp2);
    zmulmod(tmp, tmp2, n, &tmp);
    if (zscompare(tmp, 1) != 0) {		/* If inverse failed, we found a */
      printf("n factors into:\n");		/* non-trivial factor of n. */
      zwriteln(tmp2);
      zdiv(n, tmp2, &tmp, &tmp2);
      zwriteln(tmp);
      return 0;					/* and exit with failure */
    }

    zsqmod(X2, n, &m);				/* If we can add successfully */
    zsmulmod(m, 3, n, &m);			/* Finish computing m */
    zaddmod(m, a, n, &m);
    zmulmod(m, tmp2, n, &m);

  }

  else {					/* If X2 <> X3, */
    zsubmod(*X3, X2, n, &tmp);			/* Let m = (Y3 - Y2)/(X3 - X2) mod n */

    if (ziszero(tmp)) {				/* If tmp = 0 */
      printf("Addition failed with no non-trivial divisors.\n");
      return 0;					/* exit with failure */
    }
    zinv(tmp, n, &tmp2);
    zmulmod(tmp, tmp2, n, &tmp);
    if (zscompare(tmp, 1) != 0) {		/* If inverse failed, we found a */
      printf("n factors into:\n");		/* non-trivial factor of n. */
      zwriteln(tmp2);
      zdiv(n, tmp2, &tmp, &tmp2);
      zwriteln(tmp);
      return 0;					/* and exit with failure */
    }

    zsubmod(*Y3, Y2, n, &m);			/* If we can add successfully */
    zmulmod(m, tmp2, n, &m);			/* Finish computing m */
    
  }
  
  zaddmod(X2, *X3, n, &tmp);			/* Let X3' = m^2 - X2 - X3 */
  zcopy(*X3, &tmp2);
  zsqmod(m, n, X3);
  zsubmod(*X3, tmp, n, X3);
  
  zsubmod(X2, *X3, n, Y3);			/* Let Y3' = m(X2 - X3') - Y2 */
  zmulmod(m, *Y3, n, Y3);
  zsubmod(*Y3, Y2, n, Y3);

  return 1;

}

int pdouble(verylong *X, verylong *Y, verylong a, verylong n) {
  verylong X0 = 0, Y0 = 0;

  zcopy(*X, &X0);
  zcopy(*Y, &Y0);
  return padd(X0, Y0, X, Y, a, n);
}