/*
  Author:  Pate Williams (c) 1997

  Algorithm 8.5.2 (Pollard rho). See "A Course in
  Computational Algebraic Number Theory" by Henri
  Cohen page 429.
*/

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

#define BOUND 1000000l

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 trial_division(verylong *zN, struct factor *f, int *count)
{
  int e, one = 0, pr = 0;
  long B, p;
  verylong za = 0, zb = 0, zp = 0;

  zsqrt(*zN, &za, &zb);
  if (zscompare(za, BOUND) > 0) B = BOUND;
  else B = ztoint(za);
  zcopy(*zN, &za);
  zpstart2();
  do {
    p = zpnext();
    if (zsmod(za, p) == 0) {
      e = 0;
      do {
        zsdiv(za, p, &zb);
        zcopy(zb, &za);
        e++;
      } while (zsmod(za, p) == 0);
      zintoz(p, &zp);
      add_factor(e, zp, f, count);
      zcopy(za, zN);
      one = zscompare(*zN, 1l) == 0;
      if (!one) pr = zprobprime(*zN, 5l);
    }
  } while (!one && !pr && p <= B);
  if (pr) add_factor(1, *zN, f, count);
  zfree(&za);
  zfree(&zb);
  zfree(&zp);
  return p <= B;
}

void Brent(verylong zN, struct factor *f, int *count)
{
  int e, one, pr;
  long c, i, k, l;
  verylong zP = 0, za = 0, zb = 0, zg = 0, zn = 0;
  verylong zx = 0, zx1 = 0, zy = 0;

  zcopy(zN, &zn);
  do {
    c = 0, k = l = 1;
    zone(&zP);
    zintoz(2l, &zy);
    zintoz(2l, &zx);
    zintoz(2l, &zx1);
    L2:
    zsq(zx, &za);
    zsadd(za, 1l, &zb);
    zmod(zb, zn, &zx);
    zsub(zx1, zx, &za);
    zmulmod(za, zP, zn, &zb);
    zcopy(zb, &zP);
    if (++c == 20) {
      zgcd(zP, zn, &zg);
      if (zscompare(zg, 1l) > 0) goto L4;
      zcopy(zx, &zy);
      c = 0;
    }
    if (--k != 0) goto L2;
    zgcd(zP, zn, &zg);
    if (zscompare(zg, 1l) > 0) goto L4;
    zcopy(zx, &zx1);
    k = l, l *= 2;
    for (i = 0; i < k; i++) {
      zsq(zx, &za);
      zsadd(za, 1l, &zb);
      zmod(zb, zn, &zx);
    }
    zcopy(zx, &zy);
    c = 0;
    goto L2;
    L4:
    do {
      zsq(zy, &za);
      zsadd(za, 1l, &zb);
      zmod(zb, zn, &zy);
      zsub(zx1, zy, &za);
      zgcd(za, zn, &zg);
    } while (zscompare(zg, 1l) == 0);
    if (zcompare(zg, zn) == 0) {
      fprintf(stderr, "fatal error\nBrent's method failed\n");
      exit(1);
    }
    if (!zprobprime(zg, 5l)) {
      zcopy(zg, &za);
      if (!trial_division(&zg, f, count)) {
        fprintf(stderr, "fatal error\ncould not trial divide\n");
        exit(1);
      }
      zcopy(za, &zg);
      zdiv(zn, zg, &za, &zb);
      zcopy(za, &zn);
    }
    else {
      e = 0;
      do {
        zdiv(zn, zg, &za, &zb);
        zcopy(za, &zn);
        zmod(zn, zg, &za);
        e++;
      } while (zscompare(za, 0l) == 0);
      add_factor(e, zg, f, count);
    }
    one = zscompare(zn, 1l) == 0;
    if (!one) pr = zprobprime(zn, 5l);
  } while (!one && !pr);
  if (!one) add_factor(1, zn, f, count);
  zfree(&zP);
  zfree(&za);
  zfree(&zb);
  zfree(&zg);
  zfree(&zn);
  zfree(&zx);
  zfree(&zx1);
  zfree(&zy);
}

int main(void)
{
  char answer[256];
  int count, 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)) {
      printf("factors:\n");
      Brent(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;
}