/*
  Author:  Pate Williams (c) 1997

  Morrison and Brillhart's continued fraction
  factoring method. See "Prime Numbers and
  Computer Methods for Factorization" by Hans
  Riesel second edition pages 193-200.
*/

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

#define BITS_PER_LONG   32
#define BITS_PER_LONG_1 31
#define FB_MAX 256
#define FB_MAX1 257
#define N_MAX 50000
#define SIEVE_SIZE 8192

struct factor {long expon; verylong prime;};

long get_bit(long i, long *sieve)
{
  long b = i % BITS_PER_LONG;
  long c = i / BITS_PER_LONG;

  return (sieve[c] >> (BITS_PER_LONG_1 - b)) & 1;
}

void set_bit(long i, long v, long *sieve)
{
  long b = i % BITS_PER_LONG;
  long c = i / BITS_PER_LONG;
  long mask = 1 << (BITS_PER_LONG_1 - b);

  if (v == 1)
    sieve[c] |= mask;
  else
    sieve[c] &= ~mask;
}

void Sieve(long n, long *sieve)
{
  long c, i, inc;

  set_bit(0l, 0l, sieve);
  set_bit(1l, 0l, sieve);
  set_bit(2l, 1l, sieve);
  for (i = 3; i <= n; i++)
    set_bit(i, i & 1, sieve);
  c = 3;
  do {
    i = c * c, inc = c + c;
    while (i <= n) {
      set_bit(i, 0l, sieve);
      i += inc;
    }
    c += 2;
    while (!get_bit(c, sieve)) c++;
  } while (c * c <= n);
}

long get_k(verylong zN, long f_max, long k_max,
           long *f_sieve, long *k_sieve,
           long *fb, long *fb_size, verylong *zkN)
{
  long K, c_max = 0, k = 1;
  long count, f_count = 0, p;
  verylong zp = 0;

  p = 2;
  while (p <= f_max) {
    while (!get_bit(p, f_sieve)) p++;
    f_count++;
    p++;
  }
  while (k <= k_max) {
    count = 0;
    zsmul(zN, k, zkN);
    p = 1;
    do {
      while (!get_bit(p, f_sieve)) p++;
      zintoz(p, &zp);
      if (zjacobi(*zkN, zp) != - 1)
        fb[count++] = p;
      p++;
    } while (p <= f_max);
    if (count > c_max) {
      K = k;
      c_max = count;
      if (count == f_count) break;
    }
    k++;
    if (k <= k_max)
      while (!get_bit(k, k_sieve)) k++;
  }
  count = 0;
  zsmul(zN, K, zkN);
  p = 1;
  do {
    while (!get_bit(p, f_sieve)) p++;
    zintoz(p, &zp);
    if (zjacobi(*zkN, zp) != - 1)
      fb[count++] = p;
    p++;
  } while (count < f_count && p < f_max);
  *fb_size = count;
  zfree(&zp);
  return K;
}

int factor(verylong zn, long fb_size,
           long *fb, long *expon)
{
  long c, i, p;
  verylong za = 0, zb = 0;

  for (i = 0; i < fb_size; i++) expon[i] = 0;
  zcopy(zn, &za);
  if (zscompare(za, 0) < 0) {
    expon[0] = 1;
    znegate(&za);
  }
  for (i = 1; i < fb_size; i++) {
    p = fb[i];
    if (zsdiv(za, p, &zb) == 0) {
      c = 1;
      zcopy(zb, &za);
      while (zsdiv(za, p, &zb) == 0) {
        c++;
        zcopy(zb, &za);
      }
      expon[i] = c;
      if (zscompare(za, 1) == 0) {
        zfree(&za);
        zfree(&zb);
        return 1;
      }
    }
  }
  zfree(&za);
  zfree(&zb);
  return 0;
}

int continued_fraction(long fb_size, verylong zN)
{
  int factor_found = 0, found;
  long D, count = 0, i, j, k, n, s;
  long f_count = 0, f_max = 2, k_max = 10000;
  long f_sieve[SIEVE_SIZE], k_sieve[SIEVE_SIZE];
  long e, fb_size1 = fb_size + 1;
  long c[FB_MAX1], expon[FB_MAX1], fb[FB_MAX1], X[FB_MAX1];
  long e_matrix[FB_MAX1][FB_MAX], matrix[FB_MAX1][FB_MAX];
  struct factor f[32] = {{0, 0}}, t = {0, 0};
  verylong za = 0, zb = 0, zc = 0, zd = 0, zq = 0;
  verylong zx = 0, zy = 0;
  verylong zA0 = 0, zA1 = 0, zB0 = 0, zB1 = 0;
  verylong zP0 = 0, zP1 = 0, zQ0 = 0, zQ1 = 0;
  verylong zb0 = 0, zb1 = 0;
  verylong zkN = 0, zAm1 = 0, zBm1 = 0, zsqrtkN = 0;
  verylong *zA, *zQ;

  zA = calloc(fb_size1, sizeof(verylong));
  zQ = calloc(fb_size1, sizeof(verylong));
  if (!zA || !zQ) {
    fprintf(stderr, "fatal error - memory allocation\n");
    fprintf(stderr, "from continued_fraction\n");
    exit(1);
  }
  zpstart2();
  for (i = 0; i < fb_size; i++) f_max = zpnext();
  f_max++;
  Sieve(f_max, f_sieve);
  Sieve(k_max, k_sieve);
  get_k(zN, f_max, k_max, f_sieve, k_sieve,
        fb, &fb_size, &zkN);
  zsqrt(zkN, &zsqrtkN, &za);
  zcopy(zsqrtkN, &zA0);
  zone(&zAm1);
  zone(&zB0);
  zzero(&zBm1);
  zone(&zQ0);
  zzero(&zP0);
  for (n = 1; !factor_found && n <= N_MAX; n++) {
    zadd(zsqrtkN, zP0, &za);
    zdiv(za, zQ0, &zb0, &zb);
    zmul(zb0, zQ0, &za);
    zsub(za, zP0, &zP1);
    zsq(zP1, &za);
    zsub(zkN, za, &zb);
    zdiv(zb, zQ0, &zQ1, &za);
    zcopy(zP1, &zP0);
    zcopy(zQ1, &zQ0);
    zadd(zsqrtkN, zP0, &za);
    zdiv(za, zQ0, &zb1, &zb);
    zmul(zb1, zA0, &za);
    zaddmod(za, zAm1, zN, &zA1);
    zmul(zb1, zB0, &za);
    zaddmod(za, zBm1, zN, &zB1);
    zcopy(zA0, &zA[count]);
    zcopy(zA0, &zAm1);
    zcopy(zA1, &zA0);
    zcopy(zB0, &zBm1);
    zcopy(zB1, &zB0);
    zcopy(zQ1, &zq);
    if (n & 1) znegate(&zq);
    if (factor(zq, fb_size, fb, expon)) {
      printf("\b\b\b%ld", count + 1);
      zcopy(zq, &zQ[count]);
      for (i = 0; i < fb_size; i++) {
        e_matrix[count][i] = expon[i];
        matrix[count][i] = expon[i] & 1;
      }
      count++;
    }
    if (count == fb_size1) {
      for (i = 0; i < fb_size1; i++) c[i] = - 1;
      for (k = 0; !factor_found && k < fb_size1; k++) {
        found = 0;
        j = 0;
        while (!found && j < fb_size) {
          found = matrix[k][j] != 0 && c[j] < 0;
          if (!found) j++;
        }
        if (found) {
          matrix[k][j] = 1;
          for (i = 0; i < fb_size; i++) {
            if (i != j ) {
              D = matrix[k][i];
              matrix[k][i] = 0;
              for (s = k + 1; s < fb_size1; s++)
                matrix[s][i] = (matrix[s][i] + D * matrix[s][j]) & 1;
            }
          }
          c[j] = k;
        }
        else {
          for (j = 0; j < fb_size1; j++) X[j] = 0;
          for (j = 0; j < fb_size1; j++) {
            if (j == k) X[j] = 1;
            else
              for (s = 0; s < fb_size; s++)
                if (c[s] == j) X[j] = matrix[k][s];
          }
          zone(&za);
          zone(&zb);
          for (j = 0; j < fb_size1; j++) {
            if (X[j] != 0) {
              zmul(za, zA[j], &zc);
              zcopy(zc, &za);
              zmul(zb, zQ[j], &zc);
              zcopy(zc, &zb);
            }
          }
          zsqrt(zb, &zy, &zd);
          zmod(za, zN, &zx);
          zsub(zx, zy, &zc);
          zgcd(zc, zN, &zd);
          if (zscompare(zd, 1) != 0 && zcompare(zd, zN) != 0) {
            e = 0;
            do {
              e++;
              zdiv(zN, zd, &za, &zb);
              zcopy(za, &zN);
              zdiv(zN, zd, &za, &zb);
            } while (zscompare(zb, 0) == 0);
            f[f_count].expon = e;
            zcopy(zd, &f[f_count].prime);
            f_count++;
            factor_found = zscompare(zN, 1) == 0;
            if (!factor_found && zprobprime(zN, 5)) {
              factor_found = 1;
              f[f_count].expon = 1;
              zcopy(zN, &f[f_count++].prime);
            }
          }
          count = 0;
        }
      }
    }
  }
  for (i = 0; i < f_count - 1; i++) {
    for (j = i + 1; j < f_count; j++) {
      if (zcompare(f[i].prime, f[j].prime) > 0) {
        t.expon = f[i].expon;
        zcopy(f[i].prime, &t.prime);
        f[i].expon = f[j].expon;
        zcopy(f[j].prime, &f[i].prime);
        f[j].expon = t.expon;
        zcopy(t.prime, &f[j].prime);
      }
    }
  }
  printf("\nfactors:\n");
  for (i = 0; i < f_count; i++) {
    zwrite(f[i].prime);
    if (f[i].expon != 1) printf(" ^ %ld\n", f[i].expon);
    else printf("\n");
  }
  for (i = 0; i < fb_size1; i++) {
    zfree(&zA[i]);
    zfree(&zQ[i]);
  }
  free(zA);
  free(zQ);
  for (i = 0; i < f_count; i++) zfree(&f[i].prime);
  zfree(&t.prime);
  zfree(&za);
  zfree(&zb);
  zfree(&zc);
  zfree(&zd);
  zfree(&zq);
  zfree(&zx);
  zfree(&zy);
  zfree(&zA0);
  zfree(&zA1);
  zfree(&zB0);
  zfree(&zB1);
  zfree(&zP0);
  zfree(&zP1);
  zfree(&zQ0);
  zfree(&zQ1);
  zfree(&zb0);
  zfree(&zb1);
  zfree(&zkN);
  zfree(&zAm1);
  zfree(&zBm1);
  zfree(&zsqrtkN);
  return 0;
}

int main(void)
{
  long fb_size;
  verylong zN = 0;

  for (;;) {
    do {
      printf("number of factor bases (10 - %ld) = ", FB_MAX);
      scanf("%ld", &fb_size);
    } while (fb_size < 10 || fb_size > FB_MAX);
    printf("enter the number to be factored below or 0 to quit:\n");
    zread(&zN);
    if (zscompare(zN, 0) == 0) break;
    zwrite(zN);
    if (!zprobprime(zN, 5)) {
      printf(" is composite\n");
      continued_fraction(fb_size, zN);
    }
    else
      printf(" is prime\n");
  }
  zfree(&zN);
  return 0;
}