/*
  Author:  Pate Williams (c) 1997

  Algorithm 8.4.1 (Lehman). See "A Course in
  Computational Algebraic Number Theory" by
  Henri Cohen page 425.
*/

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

struct factor {int expon; verylong prime;};

void add_factor(int e, verylong zp,
                struct factor *f, int *count)
{
  int done = 0, found, i, j;

  for (i = 0; !done && i < *count; i++) {
    found = zcompare(zp, f[i].prime);
    if (found <= 0) done = 1, j = i;
  }
  if (found < 0) {
    for (i = *count - 1; i >= j; i--) {
      f[i + 1].expon = f[i].expon;
      zcopy(f[i].prime, &f[i + 1].prime);
    }
    *count = *count + 1;
    f[j].expon = e;
    zcopy(zp, &f[j].prime);
  }
  else if (found == 0) f[j].expon++;
  else {
    f[*count].expon = e;
    zcopy(zp, &f[*count].prime);
    *count = *count + 1;
  }
}

int square_test(verylong zn, verylong *zq)
{
  int square = 1;
  long q11[11], q63[63], q64[64], q65[65];
  long k, r, t;
  verylong zd = 0;

  for (k = 0; k < 11; k++) q11[k] = 0;
  for (k = 0; k < 6; k++) q11[(k * k) % 11] = 1;
  for (k = 0; k < 63; k++) q63[k] = 0;
  for (k = 0; k < 32; k++) q63[(k * k) % 63] = 1;
  for (k = 0; k < 64; k++) q64[k] = 0;
  for (k = 0; k < 32; k++) q64[(k * k) % 64] = 1;
  for (k = 0; k < 65; k++) q65[k] = 0;
  for (k = 0; k < 33; k++) q65[(k * k) % 65] = 1;
  t = zsmod(zn, 64l);
  r = zsmod(zn, 45045);
  if (q64[t] == 0) square = 0;
  if (square && q63[r % 63] == 0) square = 0;
  if (square && q65[r % 65] == 0) square = 0;
  if (square && q11[r % 11] == 0) square = 0;
  if (square) {
    zsqrt(zn, zq, &zd);
    if (zscompare(zd, 0l) != 0) square = 0;
  }
  zfree(&zd);
  return square;
}

int trial_division(long B, verylong *zN, verylong *zf)
{
  int e = 0;
  long p;
  verylong za = 0, zb = 0;

  zcopy(*zN, &za);
  zpstart2();
  do {
    p = zpnext();
    if (zsmod(za, p) == 0) {
      do {
        zsdiv(za, p, &zb);
        zcopy(zb, &za);
        e++;
      } while (zsmod(za, p) == 0);
      zintoz(p, zf);
      zcopy(za, zN);
    }
  } while (!e && p <= B);
  zfree(&za);
  zfree(&zb);
  return e;
}

int Lehman(verylong *zN, verylong *zf)
{
  long B = pow(zdoub(*zN), 1.0 / 3.0), a, k = 0, m, r;
  int e = trial_division(B, zN, zf);
  verylong za = 0, zb = 0, zc = 0, zh = 0, zl = 0;
  verylong zr = 0, zs = 0;

  if (!e) {
    L2:
    k++;
    if (k > B) goto L4;
    if (!(k & 1)) {
      zone(&zr);
      m = 2;
    }
    else {
      zsadd(*zN, k, &zr);
      m = 4;
    }
    zsmul(*zN, k, &za);
    zlshift(za, 2l, &zl);
    zintoz(B, &zh);
    zsq(zh, &zb);
    zadd(zl, zb, &zh);
    zsqrt(zl, &za, &zb);
    zsq(za, &zs);
    while (zcompare(zs, zl) < 0) {
      zsadd(za, 1l, &zb);
      zcopy(zb, &za);
      zsq(za, &zs);
    }
    r = zsmod(zr, m);
    while (!e && zcompare(zs, zh) <= 0) {
      a = zsmod(za, m);
      if (a == r) {
        zsub(zs, zl, &zc);
        if (zscompare(zc, 0l) == 0 || square_test(zc, &zb)) {
          zadd(za, zb, &zc);
          zgcd(zc, *zN, zf);
          do {
            zdiv(*zN, *zf, &za, &zb);
            zcopy(za, zN);
            zmod(za, *zf, &zb);
            e++;
          } while (zscompare(zb, 0l) == 0);
        }
      }
      if (!e) {
        zsadd(za, 1l, &zb);
        zcopy(zb, &za);
        zsq(za, &zs);
      }
    }
    if (!e) goto L2;
    L4:
  }
  zfree(&za);
  zfree(&zb);
  zfree(&zc);
  zfree(&zh);
  zfree(&zl);
  zfree(&zr);
  zfree(&zs);
  return e;
}

int main(void)
{
  char answer[256];
  int c, count, e, i;
  verylong zN = 0, zf = 0;
  struct factor f[32];

  for (i = 0; i < 32; i++) {
    f[i].expon = 0;
    f[i].prime = 0;
  }
  do {
    count = 0;
    printf("enter the number to be factored below:\n");
    zread(&zN);
    if (!zprobprime(zN, 5l)) {
      do {
        e = Lehman(&zN, &zf);
        c = zscompare(zN, 1l) == 0;
        if (e) add_factor(e, zf, f, &count);
      } while (!c && e);
      printf("factors:\n");
      if (!c) add_factor(1, zN, f, &count);
      for (i = 0; i < count; i++) {
        zwrite(f[i].prime);
        if (f[i].expon != 1) printf(" ^ %d\n", f[i].expon);
        else printf("\n");
      }
    }
    else printf("number is prime\n");
    printf("another number (n or y)? ");
    scanf("%s", answer);
  } while (tolower(answer[0]) == 'y');
  zfree(&zN);
  zfree(&zf);
  for (i = 0; i < 32; i++) zfree(&f[i].prime);
  return 0;
}