/*
  Author:  Pate Williams (c) 1997

  Algorithm 10.2.2 (Schnorr-Lenstra). See "A Course
  in Computational Algebraic Number Theory" by Henri
  Cohen pages 482-483. Factoring method for numbers
  less than equal to 10 ^ 60.
*/

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

struct zform { verylong za, zb, zc; };

int compare_zform(const void *f1, const void *f2)

{
  struct zform *zform1 = (struct zform *) f1;
  struct zform *zform2 = (struct zform *) f2;
  int i = zcompare(zform1->za, zform2->za), j;

  if (i != 0) return i;
  j = zcompare(zform1->zb, zform2->zb);
  if (j != 0) return j;
  return zcompare(zform1->zc, zform2->zc);
}

void copy_zform(struct zform zform1, struct zform *zform2)
{
  zcopy(zform1.za, &zform2->za);
  zcopy(zform1.zb, &zform2->zb);
  zcopy(zform1.zc, &zform2->zc);
}

void free_zform(struct zform *zform1)
{
  zfree(&zform1->za);
  zfree(&zform1->zb);
  zfree(&zform1->zc);
}

void get_zform(verylong *zD, struct zform *zform1)
{
  verylong zd = 0, ze = 0, zf = 0;

  printf("a = "); zread(&zform1->za);
  printf("b = "); zread(&zform1->zb);
  printf("c = "); zread(&zform1->zc);
  zmul(zform1->za, zform1->zc, &zd);
  zlshift(zd, 2l, &ze);
  zsq(zform1->zb, &zf);
  zsub(zf, ze, zD);
  printf("D = "); zwriteln(*zD);
  zfree(&zd);
  zfree(&ze);
  zfree(&zf);
}

void init_zform(struct zform *zform1)
{
  zform1->za = zform1->zb = zform1->zc = 0;
}

void put_zform(struct zform zform1)
{
  verylong zD = 0, zd = 0, ze = 0, zf = 0;

  zmul(zform1.za, zform1.zc, &zd);
  zlshift(zd, 2l, &ze);
  zsq(zform1.zb, &zf);
  zsub(zf, ze, &zD);
  printf("a = "); zwriteln(zform1.za);
  printf("b = "); zwriteln(zform1.zb);
  printf("c = "); zwriteln(zform1.zc);
  printf("D = "); zwriteln(zD);
  zfree(&zD);
  zfree(&zd);
  zfree(&ze);
  zfree(&zf);
}

struct zform *search_zform(long n, struct zform key,
                           struct zform *array)
{
  long i;
  struct zform *address;

  for (i = 0; i < n; i++) {
    address = &array[i];
    if (compare_zform(&key, address) == 0) return address;
  }
  return 0;
}

void reduce(verylong *za, verylong *zb, verylong *zc)
{
  verylong zd = 0, ze = 0, zf = 0, zq = 0, zr = 0;

  zcopy(*za, &zd);
  znegate(&zd);
  if (zcompare(zd, *zb) < 0 && zcompare(*zb, *za) <= 0)
    goto L3;
  L2:
  zlshift(*za, 1l, &zd);
  zdiv(*zb, zd, &zq, &zr);
  if (zscompare(zr, 0l) < 0) {
    if (zscompare(zd, 0l) < 0) znegate(&zd);
    zadd(zr, zd, &ze);
    znegate(&zd);
    zcopy(ze, &zr);
  }
  if (zcompare(zr, *za) > 0) {
    zsub(zr, zd, &ze);
    zcopy(ze, &zr);
    zsadd(zq, 1l, &ze);
    zcopy(ze, &zq);
  }
  zadd(*zb, zr, &zd);
  zmul(zd, zq, &ze);
  zrshift(ze, 1l, &zf);
  zsub(*zc, zf, &zd);
  zcopy(zd, zc);
  zcopy(zr, zb);
  L3:
  if (zcompare(*za, *zc) > 0) {
    znegate(zb);
    zcopy(*za, &zd);
    zcopy(*zc, za);
    zcopy(zd, zc);
    goto L2;
  }
  if (zcompare(*za, *zc) == 0 && zscompare(*zb, 0l) < 0)
    znegate(zb);
  zfree(&zd);
  zfree(&ze);
  zfree(&zf);
  zfree(&zq);
  zfree(&zr);
}

void PARTEUCL(verylong zL, verylong za,
              verylong zb, verylong *zv2,
              verylong *zv3, verylong *zd,
              verylong *zz)
{
  verylong zq = 0, zr = 0, zs = 0, zt = 0, zv = 0;
  verylong zt2 = 0, zt3 = 0;

  zzero(&zv);
  zcopy(za, zd);
  zone(zv2);
  zcopy(zb, zv3);
  zzero(zz);
  L2:
  zcopy(*zv3, &zr);
  zabs(&zr);
  if (zcompare(zr, zL) > 0) goto L3;
  if (zodd(*zz)) {
    znegate(zv2);
    znegate(zv3);
  }
  zfree(&zq);
  zfree(&zr);
  zfree(&zs);
  zfree(&zt);
  zfree(&zv);
  zfree(&zt2);
  zfree(&zt3);
  return;
  L3:
  zdiv(*zd, *zv3, &zq, &zt3);
  if (zscompare(zt3, 0l) < 0) {
    zadd(zt3, zr, &zt);
    zcopy(zt, &zt3);
  }
  zmul(zq, *zv2, &zs);
  zsub(zv, zs, &zt2);
  zcopy(*zv2, &zv);
  zcopy(*zv3, zd);
  zcopy(zt2, zv2);
  zcopy(zt3, zv3);
  zsadd(*zz, 1l, &zt);
  zcopy(zt, zz);
  goto L2;
}

void NUDUPL(verylong zL, struct zform zform1,
            struct zform *zform2)
{
  verylong zA = 0, zB = 0, zC = 0, za = 0, zb = 0;
  verylong zc = 0, zd = 0, ze = 0, zg = 0, zr = 0;
  verylong zs = 0, zt = 0, zu = 0, zv = 0, zz = 0;
  verylong zC1 = 0, za2 = 0, zb2 = 0, zc2 = 0;
  verylong zd1 = 0, zv2 = 0, zv3 = 0;

  zcopy(zform1.za, &za);
  zcopy(zform1.zb, &zb);
  zcopy(zform1.zc, &zc);
  zexteucl(zb, &zu, za, &zv, &zd1);
  zdiv(za, zd1, &zA, &zr);
  zdiv(zb, zd1, &zB, &zr);
  zcopy(zc, &zr);
  znegate(&zr);
  zmulmod(zr, zu, zA, &zC);
  zsub(zA, zC, &zC1);
  if (zcompare(zC1, zC) < 0) {
    zcopy(zC1, &zC);
    znegate(&zC);
  }
  PARTEUCL(zL, zA, zC, &zv2, &zv3, &zd, &zz);
  if (zscompare(zz, 0l) == 0) {
    zmul(zB, zv3, &zr);
    zadd(zr, zc, &zs);
    zdiv(zs, zd, &zg, &zr);
    zsq(zd, &za2);
    zsq(zv3, &zc2);
    zadd(zd, zv3, &zr);
    zsq(zr, &zs);
    zadd(zb, zs, &zt);
    zsub(zt, za2, &zr);
    zsub(zr, zc2, &zb2);
    zmul(zg, zd1, &zr);
    zadd(zc2, zr, &zs);
    zcopy(zs, &zc2);
  }
  else {
    zmul(zc, zv, &zr);
    zmul(zB, zd, &zs);
    zadd(zr, zs, &zt);
    zdiv(zt, zA, &ze, &zr);
    zmul(ze, zv2, &zr);
    zsub(zr, zB, &zs);
    zdiv(zs, zv, &zg, &zt);
    zmul(zv, zg, &zs);
    zadd(zr, zs, &zb2);
    if (zscompare(zd1, 1l) > 0) {
      zmul(zd1, zb2, &zr);
      zcopy(zr, &zb2);
      zmul(zd1, zv, &zr);
      zcopy(zr, &zv);
      zmul(zd1, zv2, &zr);
      zcopy(zr, &zv2);
    }
    zsq(zd, &za2);
    zsq(zv3, &zc2);
    zadd(zd, zv3, &zr);
    zsq(zr, &zs);
    zadd(zb2, zs, &zr);
    zsub(zr, za2, &zs);
    zsub(zs, zc2, &zb2);
    zmul(ze, zv, &zr);
    zadd(za2, zr, &zs);
    zcopy(zs, &za2);
    zmul(zg, zv2, &zr);
    zadd(zc2, zr, &zs);
    zcopy(zs, &zc2);
  }
  reduce(&za2, &zb2, &zc2);
  zcopy(za2, &zform2->za);
  zcopy(zb2, &zform2->zb);
  zcopy(zc2, &zform2->zc);
  zfree(&zA);
  zfree(&zB);
  zfree(&zC);
  zfree(&za);
  zfree(&zb);
  zfree(&zc);
  zfree(&zd);
  zfree(&ze);
  zfree(&zg);
  zfree(&zr);
  zfree(&zs);
  zfree(&zt);
  zfree(&zu);
  zfree(&zv);
  zfree(&zz);
  zfree(&zC1);
  zfree(&za2);
  zfree(&zb2);
  zfree(&zc2);
  zfree(&zd1);
  zfree(&zv2);
  zfree(&zv3);
}

void composition(struct zform zform1, struct zform zform2,
                 struct zform *zform3)
{
  verylong za1 = 0, za2 = 0, zb1 = 0, zb2 = 0;
  verylong zc1 = 0, zc2 = 0, zd1 = 0, zv1 = 0;
  verylong zv2 = 0, zx1 = 0, zx2 = 0, zy1 = 0;
  verylong zy2 = 0;
  verylong za = 0, zb = 0, zc = 0, zd = 0, zn = 0;
  verylong zr = 0, zs = 0, zu = 0, zv = 0;

  zcopy(zform1.za, &za1);
  zcopy(zform1.zb, &zb1);
  zcopy(zform1.zc, &zc1);
  zcopy(zform2.za, &za2);
  zcopy(zform2.zb, &zb2);
  zcopy(zform2.zc, &zc2);
  if (zcompare(za1, za2) > 0) {
    /* excange forms 1 and 2 */
    zcopy(zform2.za, &za1);
    zcopy(zform2.zb, &zb1);
    zcopy(zform2.zc, &zc1);
    zcopy(zform1.za, &za2);
    zcopy(zform1.zb, &zb2);
    zcopy(zform1.zc, &zc2);
  }
  zadd(zb1, zb2, &za);
  zrshift(za, 1l, &zs);
  zsub(zb2, zs, &zn);
  zmod(za2, za1, &zr);
  if (zscompare(zr, 0l) == 0) {
    zzero(&zy1);
    zcopy(za1, &zd);
  }
  else {
    zexteucl(za2, &zu, za1, &zv, &zd);
    zcopy(zu, &zy1);
  }
  zmod(zs, zd, &zr);
  if (zscompare(zr, 0l) == 0) {
    zone(&zy2);
    znegate(&zy2);
    zzero(&zx2);
    zcopy(zd, &zd1);
  }
  else {
    zexteucl(zs, &zu, zd, &zv, &zd1);
    zcopy(zu, &zx2);
    zcopy(zv, &zy2);
    znegate(&zy2);
  }
  zdiv(za1, zd1, &zv1, &zr);
  zdiv(za2, zd1, &zv2, &zr);
  zmul(zy1, zy2, &za);
  zmul(za, zn, &zb);
  zmul(zx2, zc2, &zc);
  zsubmod(zb, zc, zv1, &zr);
  zmul(zv2, zr, &za);
  zlshift(za, 1l, &zb);
  zmul(zv1, zv2, &zform3->za);
  zadd(zb2, zb, &zform3->zb);
  zmul(zc2, zd1, &zb);
  zadd(zb2, za, &zc);
  zmul(zc, zr, &zd);
  zadd(zb, zd, &za);
  zdiv(za, zv1, &zform3->zc, &zr);
  reduce(&zform3->za, &zform3->zb, &zform3->zc);
  zfree(&za1);
  zfree(&za2);
  zfree(&zb1);
  zfree(&zb2);
  zfree(&zc1);
  zfree(&zc2);
  zfree(&zd1);
  zfree(&zv1);
  zfree(&zv2);
  zfree(&zx1);
  zfree(&zx2);
  zfree(&zy1);
  zfree(&zy2);
  zfree(&za);
  zfree(&zb);
  zfree(&zc);
  zfree(&zd);
  zfree(&zn);
  zfree(&zr);
  zfree(&zs);
  zfree(&zu);
  zfree(&zv);
}

void zform_pow(long e, struct zform x,
               struct zform *y,
               verylong zL)
{
  struct zform s, t, u;

  init_zform(&s);
  init_zform(&t);
  init_zform(&u);
  if (e == 1) copy_zform(x, y);
  else {
    copy_zform(x, &t);
    zintoz(1l, &s.za);
    zcopy(x.zb, &s.zb);
    znegate(&s.zb);
    zmul(x.za, x.zc, &s.zc);
    reduce(&s.za, &s.zb, &s.zc);
    while (e > 0) {
      if (e & 1) {
        composition(s, t, &u);
        copy_zform(u, &s);
      }
      e >>= 1;
      if (e > 0) {
        NUDUPL(zL, t, &u);
        copy_zform(u, &t);
      }
    }
    copy_zform(s, y);
  }
  free_zform(&s);
  free_zform(&t);
  free_zform(&u);
}

void zsquare_root(verylong za, verylong zp, verylong *zr)
{
  long i, s = 0, t;
  verylong zb = 0, zc = 0, zd = 0, ze = 0, zi = 0;
  verylong zq = 0, zs = 0, zt = 0, zx = 0, zy = 0;

  if (zcompare(za, zp) < 0) zcopy(za, &ze);
  else
    zmod(za, zp, &ze);
  if (zscompare(ze, 0l) == 0)
    zzero(zr);
  else if (zjacobi(ze, zp) == - 1)
    zzero(zr);
  else {
    do
      zrandomb(zp, &zb);
    while (zscompare(zb, 0l) == 0 || zjacobi(zb, zp) != - 1);
    zsadd(zp, - 1l, &zq);
    zcopy(zq, &zy);
    do {
      zrshift(zy, 1l, &zt);
      s++;
      t = zodd(zt);
      zcopy(zt, &zy);
    } while (!t);
    zinvmod(ze, zp, &zi);
    zexpmod(zb, zt, zp, &zc);
    zsadd(zt, 1l, &zs);
    zrshift(zs, 1l, &zx);
    zexpmod(ze, zx, zp, zr);
    for (i = 1; i < s; i++) {
      zsq(*zr, &zs);
      zmul(zs, zi, &zx);
      zsexpmod(zx, pow(2, s - i - 1), zp, &zd);
      if (zcompare(zd, zq) == 0) {
        zmulmod(*zr, zc, zp, &zx);
        zcopy(zx, zr);
      }
      zmulmod(zc, zc, zp, &zx);
      zcopy(zx, &zc);
    }
  }
  zfree(&zb);
  zfree(&zc);
  zfree(&zd);
  zfree(&ze);
  zfree(&zi);
  zfree(&zq);
  zfree(&zs);
  zfree(&zt);
  zfree(&zx);
}

void zsquare_roots(verylong za, verylong zq,
                   verylong *zx1, verylong *zx2,
                   verylong *zy1, verylong *zy2)
{
  long i, n, r;
  verylong zb = 0, zc = 0, zd = 0, ze = 0, zg = 0;
  verylong zn = 0, zp = 0, zr = 0, zs = 0, zx = 0;
  verylong zy = 0, zz = 0;

  zintoz(4l, &zp);
  zcopy(za, &zn);
  n = zsmod(zn, 4l);
  for (i = 1; i < 4; i++)
    if ((i * i) % 4 == n) r = i;
  zintoz(r, &zr);
  if (zcompare(za, zq) < 0)
    zcopy(za, &zn);
  else
    zmod(za, zq, &zn);
  zsquare_root(zn, zq, &zs);
  zexteucl(zp, &zc, zq, &zd, &zg);
  if (zscompare(zg, 1l) == 0) {
    zmul(zp, zq, &zn);
    zmul(zr, zq, &zb);
    zmul(zs, zc, &ze);
    zmul(ze, zp, &zz);
    zaddmod(zb, zz, zn, &zx);
    zmul(zr, zq, &zb);
    zmul(zs, zc, &ze);
    zmul(ze, zp, &zz);
    zsubmod(zb, zz, zn, &zy);
    zmod(zx, zn, zx1);
    znegate(&zx);
    zmod(zx, zn, zx2);
    zmod(zy, zn, zy1);
    znegate(&zy);
    zmod(zy, zn, zy2);
  }
  else {
    zzero(zx1);
    zzero(zx2);
    zzero(zy1);
    zzero(zy2);
  }
  zfree(&zb);
  zfree(&zc);
  zfree(&zd);
  zfree(&ze);
  zfree(&zg);
  zfree(&zn);
  zfree(&zp);
  zfree(&zr);
  zfree(&zs);
  zfree(&zx);
  zfree(&zy);
  zfree(&zz);
}

int Schnorr_Lenstra(verylong *zN, verylong *zf)
{
  double L, logz, z;
  int expon = 0, found;
  long B, K = 1, c, e, e1, i, k = 0, l, *p, q, q1;
  struct zform unit, x, y;
  verylong zD = 0, zE = 0, zL = 0, zKN = 0;
  verylong za = 0, zb = 0, zc = 0, zp = 0, zq = 0;
  verylong zx2 = 0, zy1 = 0, zy2 = 0;

  init_zform(&unit);
  init_zform(&x);
  init_zform(&y);
  z = zdoub(*zN);
  logz = log(z);
  L = exp(sqrt(logz * log(logz)));
  B = sqrt(L);
  e = log(sqrt(L)) / log(2.0);
  zpstart2();
  do {
    q = zpnext();
    k++;
  } while (q <= B);
  p = calloc(k, sizeof(long));
  if (p == 0) {
    fprintf(stderr, "fatal error\ninsufficient memory\n");
    exit(1);
  }
  zpstart2();
  for (i = 0; i < k; i++) p[i] = zpnext();
  zpstart();
  L2:
  zsmul(*zN, K, &zKN);
  if (zsmod(zKN, 4l) == 3) zcopy(zKN, &zD);
  else zlshift(zKN, 2l, &zD);
  zintoz(pow(zdoub(zD), 0.25), &zL);
  znegate(&zD);
  L3:
  do zintoz(zpnext(), &zp); while (zjacobi(zD, zp) != 1);
  zlshift(zp, 2l, &zq);
  zmod(zD, zq, &zE);
  zsquare_roots(zE, zp, &zb, &zx2, &zy1, &zy2);
  zcopy(zp, &x.za);
  zcopy(zb, &x.zb);
  zsq(zb, &za);
  zsub(za, zD, &zc);
  zdiv(zc, zq, &x.zc, &za);
  zone(&unit.za);
  zcopy(x.zb, &unit.zb);
  znegate(&unit.zb);
  zmul(x.za, x.zc, &unit.zc);
  reduce(&unit.za, &unit.zb, &unit.zc);
  c = 0, i = 0;
  L4:
  i++;
  if (i == k) {
    K++;
    goto L2;
  }
  q = p[i];
  q1 = q;
  l = B / q;
  while (q1 <= l) q1 *= q;
  zform_pow(q1, x, &y, zL);
  copy_zform(y, &x);
  if (i < k - 1 && ++c < 20) goto L4;
  e1 = 0;
  found = 0;
  while (!found && e1 < e) {
    composition(x, x, &y);
    found = compare_zform(&y, &unit) == 0;
    if (!found) {
      copy_zform(y, &x);
      e1++;
    }
  }
  if (!found) {
    c = 0;
    goto L4;
  }
  if (zscompare(x.za, 1l) == 0 ||
      zcompare(x.za, *zN) == 0) goto L3;
  zmod(*zN, x.za, &za);
  if (zscompare(za, 0l) != 0) goto L3;
  zcopy(x.za, zf);
  do {
    expon++;
    zdiv(*zN, *zf, &za, &zb);
    zcopy(za, zN);
    zmod(*zN, *zf, &zb);
  } while (zscompare(zb, 0l) == 0);
  free_zform(&unit);
  free_zform(&x);
  free_zform(&y);
  zfree(&zD);
  zfree(&zE);
  zfree(&zL);
  zfree(&zKN);
  zfree(&za);
  zfree(&zb);
  zfree(&zc);
  zfree(&zp);
  zfree(&zq);
  zfree(&zx2);
  zfree(&zy1);
  zfree(&zy2);
  return expon;
}

int main(void)
{
  char answer[256];
  int e, one;
  verylong zN = 0, zf = 0;

  do {
    printf("enter the number to be factored below:\n");
    zread(&zN);
    printf("factors:\n");
    do {
      e = Schnorr_Lenstra(&zN, &zf);
      zwrite(zf);
      if (e > 1) printf(" ^ %d\n", e);
      else printf("\n");
      one = zscompare(zN, 1l) == 0;
    } while (!one && !zprobprime(zN, 5l));
    if (!one) zwriteln(zN);
    printf("another number (n or y)? ");
    scanf("%s", answer);
  } while (tolower(answer[0]) == 'y');
  zfree(&zN);
  zfree(&zf);
  return 0;
}