/*
  Author:  Pate Williams (c) 1997

  Algorithm 10.3.3 (Lenstra's ECM). See "A Course
  in Computational Algebraic Number Theory" by
  Henri Cohen page 488.
*/

#include <ctype.h>
#include <stdio.h>
#include <stdlib.h>
#include "lip.h"

#define CURVES 50l

struct point {verylong zx, zy, zz;};

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);
  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;
}

int LenstrasECM(verylong *zN, verylong *zg)
{
  int expon = 0, found;
  long B = 2000000l, i, j, k = 0, l, *p, q, q1;
  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) {
    fprintf(stderr, "fatal error\ninsufficient memory\n");
    exit(1);
  }
  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 == 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;
  }
  do {
    zdiv(*zN, *zg, &zb, &zd);
    zcopy(zb, zN);
    zmod(*zN, *zg, &zd);
    expon++;
  } while (zscompare(zd, 0l) == 0);
  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 main(void)
{
  int cant, expon, one, pri;
  char answer[256];
  verylong zN = 0, zd = 0, zg = 0;

  do {
    printf("enter the number to be factored below:\n");
    zread(&zN);
    if (zprobprime(zN, 5l))
      printf("number is prime\n");
    else {
      printf("number is composite\n");
      zintoz(6l, &zg);
      zgcd(zN, zg, &zd);
      if (zscompare(zd, 1l) != 0)
        printf("number can't be factored using Lenstra's ECM\n");
      else {
        printf("factors:\n");
        cant = 0, pri = 0;
        do {
          expon = LenstrasECM(&zN, &zg);
          zwrite(zg);
          if (expon == 1) printf("\n");
          else printf(" ^ %d\n", expon);
          one = zscompare(zN, 1l) == 0;
          if (!one) pri = zprobprime(zN, 5l);
          zintoz(6l, &zg);
          zgcd(zN, zg, &zd);
          if (zscompare(zd, 1l) != 0) {
            printf("number can't be factored using Lenstra's ECM\n");
            cant = 1;
          }
        } while (!cant && !one && !pri);
        if (pri) zwriteln(zN);
        if (cant) {
          zwriteln(zN);
          printf("last factor is composite\n");
        }
      }
    }
    printf("another number (n or y)? ");
    scanf("%s", answer);
  } while (tolower(answer[0]) == 'y');
  zfree(&zN);
  zfree(&zd);
  zfree(&zg);
  return 0;
}