/*
  Author:  Pate Williams (c) 1997

  Special number field sieve factorization of
  r ^ 3 + 2. See "Prime Numbers and Computer
  Methods for Factorization" by Hans Riesel
  second edition pages 214-216.
*/

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

/* default: -200, +200, +100, 10, 10, 20, 21 */
/* crashes: -500, +500, +250, 15, 15, 30, 31 */

#define A_MIN -200
#define A_MAX +200
#define B_MAX +100
#define FB1_SIZE 10  /* number of rational primes */
#define FB2_SIZE 10  /* number of algebraic primes */
#define FB_SIZE  20  /* FB1_SIZE + FB2_SIZE */
#define FB_SIZE1 21  /* FB_SIZE + 1 */

struct element {long a, b, c;};
struct zelement {verylong a, b, c;};
struct factor {long expon; verylong prime;};

void norm(verylong zN, struct zelement e, verylong *nm)
{
  verylong a = 0, b = 0, c = 0;
  verylong za = 0, zb = 0, zc = 0;

  zcopy(e.a, &a);
  zcopy(e.b, &b);
  zcopy(e.c, &c);
  zsexpmod(a, 3, zN, &za);
  zsexpmod(b, 3, zN, &zb);
  zsmulmod(zb, 2, zN, &zc);
  zsubmod(za, zc, zN, &zb);
  zsexpmod(c, 3, zN, &za);
  zsmulmod(za, 4, zN, &zc);
  zaddmod(zb, zc, zN, &za);
  zmulmod(a, b, zN, &zb);
  zmulmod(zb, c, zN, &zc);
  zsmulmod(zc, 6, zN, &zb);
  zaddmod(za, zb, zN, nm);
  zfree(&a);
  zfree(&b);
  zfree(&c);
  zfree(&za);
  zfree(&zb);
  zfree(&zc);
}

void multiply(verylong zN, struct zelement e1,
              struct zelement e2,
              struct zelement *product)
{
  verylong za = 0, zb = 0, zc = 0;

  zmulmod(e1.a, e2.a, zN, &za);
  zmulmod(e1.b, e2.c, zN, &zb);
  zsmulmod(zb, 2, zN, &zc);
  zsubmod(za, zc, zN, &zb);
  zmulmod(e1.c, e2.b, zN, &za);
  zsmulmod(za, 2, zN, &zc);
  zsubmod(zb, zc, zN, &product->a);
  zmulmod(e1.a, e2.b, zN, &za);
  zmulmod(e1.b, e2.a, zN, &zb);
  zaddmod(za, zb, zN, &zc);
  zmulmod(e1.c, e2.c, zN, &za);
  zsmulmod(za, 2, zN, &zb);
  zsubmod(zc, zb, zN, &product->b);
  zmulmod(e1.a, e2.c, zN, &za);
  zmulmod(e1.b, e2.b, zN, &zb);
  zaddmod(za, zb, zN, &zc);
  zmulmod(e1.c, e2.a, zN, &za);
  zaddmod(zc, za, zN, &product->c);
  zfree(&za);
  zfree(&zb);
  zfree(&zc);
}

void power(verylong zN, struct zelement x,
           long b, struct zelement *a)
/* returns a = x ^ b */
{
  struct zelement s = {0, 0, 0}, t = {0, 0, 0};

  zone(&a->a);
  zzero(&a->b);
  zzero(&a->c);
  zcopy(x.a, &s.a);
  zcopy(x.b, &s.b);
  zcopy(x.c, &s.c);
  while (b != 0) {
    if (b & 1) {
      multiply(zN, *a, s, &t);
      zcopy(t.a, &a->a);
      zcopy(t.b, &a->b);
      zcopy(t.c, &a->c);
    }
    b >>= 1;
    if (b != 0) {
      multiply(zN, s, s, &t);
      zcopy(t.a, &s.a);
      zcopy(t.b, &s.b);
      zcopy(t.c, &s.c);
    }
  }
  zfree(&s.a);
  zfree(&s.b);
  zfree(&s.c);
  zfree(&t.a);
  zfree(&t.b);
  zfree(&t.c);
}

int rational_factor(long a, long b,
                    verylong zr, long *FB1,
                    long *expon)
{
  long c, i, p;
  verylong za = 0, zn = 0;

  zsmul(zr, b, &za);
  zsadd(za, a, &zn);
  for (i = 0; i < FB1_SIZE; i++) expon[i] = 0;
  if (zscompare(zn, 0) < 0) {
    expon[0] = 1;
    znegate(&zn);
  }
  if (zscompare(zn, 1) <= 0) return 0;
  for (i = 1; i < FB1_SIZE; i++) {
    p = FB1[i];
    if (zsdiv(zn, p, &za) == 0) {
      c = 1;
      zcopy(za, &zn);
      while (zsdiv(zn, p, &za) == 0) {
        c++;
        zcopy(za, &zn);
      }
      expon[i] = c;
      if (zscompare(zn, 1) == 0) {
        zfree(&za);
        zfree(&zn);
        return 1;
      }
    }
  }
  zfree(&za);
  zfree(&zn);
  return 0;
}

int divide(verylong zN,
           struct zelement alpha,
           struct zelement beta,
           struct zelement *gamma)
{
  int flag = 0;
  verylong a = 0, b = 0, c = 0, d = 0, e = 0;
  verylong f = 0, g = 0, h = 0, i = 0, x = 0;
  verylong aa = 0, ab = 0, ac = 0;
  verylong ba = 0, bb = 0, bc = 0;
  verylong dd = 0, ee = 0, ff = 0;
  verylong za = 0, zb = 0, zc = 0, zd = 0, ze = 0;
  verylong zf = 0;
  verylong zN1 = 0;

  zcopy(alpha.a, &aa);
  zcopy(alpha.b, &ab);
  zcopy(alpha.c, &ac);
  zcopy(beta.a, &ba);
  zcopy(beta.b, &bb);
  zcopy(beta.c, &bc);
  zcopy(aa, &g);
  zcopy(ab, &h);
  zcopy(ac, &i);
  zcopy(ba, &d);
  zcopy(bb, &e);
  zcopy(bc, &f);
  if (zcompare(aa, ba) == 0 && zcompare(ab, bb) == 0 &&
      zcompare(ac, bc) == 0) {
    zfree(&d);
    zfree(&e);
    zfree(&f);
    zfree(&g);
    zfree(&h);
    zfree(&i);
    zfree(&aa);
    zfree(&ab);
    zfree(&ac);
    zfree(&ba);
    zfree(&bb);
    zfree(&bc);
    zone(&gamma->a);
    zzero(&gamma->b);
    zzero(&gamma->c);
    return 1;
  }
  znegate(&aa);
  znegate(&ab);
  znegate(&ac);
  if (zcompare(aa, ba) == 0 && zcompare(ab, bb) == 0 &&
      zcompare(ac, bc) == 0) {
    zfree(&d);
    zfree(&e);
    zfree(&f);
    zfree(&g);
    zfree(&h);
    zfree(&i);
    zfree(&aa);
    zfree(&ab);
    zfree(&ac);
    zfree(&ba);
    zfree(&bb);
    zfree(&bc);
    zone(&gamma->a);
    znegate(&gamma->a);
    zzero(&gamma->b);
    zzero(&gamma->c);
    return 1;
  }
  zsexpmod(d, 3, zN, &za);
  zsexpmod(e, 3, zN, &zb);
  zsexpmod(f, 3, zN, &zc);
  zsmulmod(zb, 2, zN, &zd);
  zsmulmod(zc, 4, zN, &ze);
  zsubmod(za, zd, zN, &zf);
  zaddmod(zf, ze, zN, &za);
  zmulmod(d, e, zN, &zb);
  zmulmod(zb, f, zN, &zc);
  zsmulmod(zc, 8, zN, &zd);
  zaddmod(za, zd, zN, &x);
  zmulmod(d, d, zN, &dd);
  zmulmod(e, e, zN, &ee);
  zmulmod(f, f, zN, &ff);
  zmulmod(dd, g, zN, &za);
  zmulmod(ee, h, zN, &zb);
  zsmulmod(zb, 2, zN, &zc);
  zsubmod(za, zc, zN, &zb);
  zmulmod(ff, i, zN, &za);
  zsmulmod(za, 4, zN, &zc);
  zaddmod(zb, zc, zN, &za);
  zmulmod(d, f, zN, &zb);
  zmulmod(zb, h, zN, &zc);
  zsmulmod(zc, 2, zN, &zd);
  zaddmod(za, zd, zN, &zb);
  zmulmod(e, f, zN, &za);
  zmulmod(za, g, zN, &zc);
  zsmulmod(zc, 2, zN, &zd);
  zaddmod(zb, zd, zN, &za);
  zmulmod(d, e, zN, &zb);
  zmulmod(zb, i, zN, &zc);
  zsmulmod(zc, 2, zN, &zd);
  zaddmod(za, zd, zN, &a);
  zmulmod(dd, h, zN, &za);
  zmulmod(ee, i, zN, &zb);
  zsmulmod(zb, 2, zN, &zc);
  zsubmod(za, zc, zN, &zb);
  zmulmod(ff, g, zN, &za);
  zsmulmod(za, 2, zN, &zc);
  zsubmod(zb, zc, zN, &za);
  zmulmod(d, e, zN, &zb);
  zmulmod(zb, g, zN, &zc);
  zsubmod(za, zc, zN, &zb);
  zmulmod(d, f, zN, &za);
  zmulmod(za, i, zN, &zc);
  zsmulmod(zc, 2, zN, &zd);
  zaddmod(zb, zd, zN, &za);
  zmulmod(e, f, zN, &zb);
  zmulmod(zb, h, zN, &zc);
  zsmulmod(zc, 2, zN, &zd);
  zaddmod(za, zd, zN, &b);
  zmulmod(dd, i, zN, &za);
  zmulmod(ee, g, zN, &zb);
  zaddmod(za, zb, zN, &zc);
  zmulmod(ff, h, zN, &za);
  zsmulmod(za, 2, zN, &zb);
  zsubmod(zc, zb, zN, &za);
  zmulmod(e, f, zN, &zb);
  zmulmod(zb, i, zN, &zc);
  zsmulmod(zc, 2, zN, &zb);
  zaddmod(za, zb, zN, &zc);
  zmulmod(d, e, zN, &za);
  zmulmod(za, h, zN, &zb);
  zsubmod(zc, zb, zN, &za);
  zmulmod(d, f, zN, &zb);
  zmulmod(zb, g, zN, &zc);
  zsubmod(za, zc, zN, &c);
  zsadd(zN, - 1, &zN1);
  if (zscompare(x, 0)  != 0 && zscompare(x, 1) != 0 &&
      zcompare(x, zN1) != 0) {
    zrshift(zN, 1, &zf);
    if (zcompare(x, zf) > 0) {
      zsub(x, zN, &za);
      zcopy(za, &x);
    }
    if (zcompare(a, zf) > 0) {
      zsub(a, zN, &za);
      zcopy(za, &zb);
      zabs(&zb);
      if (zcompare(zb, zN) == 0) zzero(&za);
      zcopy(za, &a);
    }
    if (zcompare(b, zf) > 0) {
      zsub(b, zN, &za);
      zcopy(za, &zb);
      zabs(&zb);
      if (zcompare(zb, zN) == 0) zzero(&za);
      zcopy(za, &b);
    }
    if (zcompare(c, zf) > 0) {
      zsub(c, zN, &za);
      zcopy(za, &zb);
      zabs(&zb);
      if (zcompare(zb, zN) == 0) zzero(&za);
      zcopy(za, &c);
    }
    zdiv(a, x, &za, &zb);
    if (zscompare(zb, 0) == 0) {
      zcopy(za, &a);
      zdiv(b, x, &za, &zb);
      if (zscompare(zb, 0) == 0) {
        zcopy(za, &b);
        zdiv(c, x, &za, &zb);
        if (zscompare(zb, 0) == 0) {
          zcopy(za, &c);
          if (zscompare(a, 0) < 0) {
            zadd(a, zN, &za);
            zmod(za,  zN, &zb);
            zcopy(zb, &a);
          }
          if (zscompare(b, 0) < 0) {
            zadd(b, zN, &za);
            zmod(za, zN, &zb);
            zcopy(zb, &b);
          }
          if (zscompare(c, 0) < 0) {
            zadd(c, zN, &za);
            zmod(za, zN, &zb);
            zcopy(zb, &c);
          }
          flag = 1;
        }
      }
    }
  }
  else {
    zzero(&a);
    zzero(&b);
    zzero(&c);
  }
  zcopy(a, &gamma->a);
  zcopy(b, &gamma->b);
  zcopy(c, &gamma->c);
  zfree(&a);
  zfree(&b);
  zfree(&c);
  zfree(&d);
  zfree(&e);
  zfree(&f);
  zfree(&g);
  zfree(&h);
  zfree(&i);
  zfree(&x);
  zfree(&aa);
  zfree(&ab);
  zfree(&ac);
  zfree(&ba);
  zfree(&bb);
  zfree(&bc);
  zfree(&dd);
  zfree(&ee);
  zfree(&ff);
  zfree(&za);
  zfree(&zb);
  zfree(&zc);
  zfree(&zd);
  zfree(&ze);
  zfree(&zf);
  zfree(&zN1);
  return flag;
}

int algebraic_factor(verylong zN, struct zelement *FB2,
                     long *expon, struct zelement e)
{
  int flag = 0;
  long count, i;
  static struct element u[8] = {{- 1, 1, - 1},
                                {5, - 4, 3},
                                {- 19, 15, - 12},
                                {73, - 58, 46},
                                {- 281, 223, - 177},
                                {1081, - 858, 681},
                                {- 4159, 3301, - 2620},
                                {16001, - 12700, 10080}};
  struct zelement n = {0, 0, 0}, p = {0, 0, 0};
  struct zelement q = {0, 0, 0}, unit = {0, 0, 0};
  verylong nm = 0, za = 0, zb = 0, zc = 0;
  verylong zN1 = 0;

  for (i = 0; i < FB2_SIZE; i++) expon[i] = 0;
  zsadd(zN, - 1, &zN1);
  norm(zN, e, &nm);
  if (zscompare(nm, 0) == 0) {
    zfree(&nm);
    return 0;
  }
  zcopy(e.a, &n.a);
  zcopy(e.b, &n.b);
  zcopy(e.c, &n.c);
  for (i = 2; i < FB2_SIZE; i++) {
    zcopy(FB2[i].a, &p.a);
    zcopy(FB2[i].b, &p.b);
    zcopy(FB2[i].c, &p.c);
    if (divide(zN, n, p, &q)) {
      count = 1;
      zcopy(q.a, &n.a);
      zcopy(q.b, &n.b);
      zcopy(q.c, &n.c);
      while (divide(zN, n, p, &q)) {
        count++;
        zcopy(q.a, &n.a);
        zcopy(q.b, &n.b);
        zcopy(q.c, &n.c);
      }
      expon[i] = count;
      if (zscompare(n.a, 1) == 0 &&
          zscompare(n.b, 0) == 0 &&
          zscompare(n.c, 0) == 0) {
        expon[0] = 0;
        zfree(&nm);
        zfree(&n.a);
        zfree(&n.b);
        zfree(&n.c);
        zfree(&p.a);
        zfree(&p.b);
        zfree(&p.c);
        zfree(&q.a);
        zfree(&q.c);
        return 1;
      }
      if (zcompare(n.a, zN1) == 0 &&
          zscompare(n.b, 0)  == 0 &&
          zscompare(n.c, 0)  == 0) {
        expon[0] = 1;
        zfree(&nm);
        zfree(&n.a);
        zfree(&n.b);
        zfree(&n.c);
        zfree(&p.a);
        zfree(&p.b);
        zfree(&p.c);
        zfree(&q.a);
        zfree(&q.c);
        return 1;
      }
    }
  }
  norm(zN, n, &nm);
  if (zscompare(nm, 1) == 0 || zcompare(nm, zN1) == 0) {
    for (i = 0; i < 8; i++) {
      zintoz(u[i].a, &za);
      zintoz(u[i].b, &zb);
      zintoz(u[i].c, &zc);
      if (zcompare(n.a, unit.a) == 0 &&
          zcompare(n.b, unit.b) == 0 &&
          zcompare(n.c, unit.c) == 0) {
        zfree(&nm);
        zfree(&za);
        zfree(&zb);
        zfree(&zc);
        zfree(&n.a);
        zfree(&n.b);
        zfree(&n.c);
        zfree(&p.a);
        zfree(&p.b);
        zfree(&p.c);
        zfree(&q.a);
        zfree(&q.c);
        zfree(&unit.a);
        zfree(&unit.b);
        zfree(&unit.c);
        return 0;
      }
      znegate(&n.a);
      znegate(&n.b);
      znegate(&n.c);
      if (zcompare(n.a, unit.a) == 0 &&
          zcompare(n.b, unit.b) == 0 &&
          zcompare(n.c, unit.c) == 0) {
        zfree(&nm);
        zfree(&za);
        zfree(&zb);
        zfree(&zc);
        zfree(&n.a);
        zfree(&n.b);
        zfree(&n.c);
        zfree(&p.a);
        zfree(&p.b);
        zfree(&p.c);
        zfree(&q.a);
        zfree(&q.c);
        zfree(&unit.a);
        zfree(&unit.b);
        zfree(&unit.c);
        return 0;
      }
      znegate(&n.a);
      znegate(&n.b);
      znegate(&n.c);
    }
    zcopy(FB2[1].a, &p.a);
    zcopy(FB2[1].b, &p.b);
    zcopy(FB2[1].c, &p.c);
    if (divide(zN, n, p, &q)) {
      count = 1;
      zcopy(q.a, &n.a);
      zcopy(q.b, &n.b);
      zcopy(q.c, &n.c);
      if ((zscompare(n.a, 1) != 0  ||
          zcompare(n.a, zN1) != 0) ||
          zscompare(n.b, 0)  != 0  ||
          zscompare(n.c, 0)  != 0) {
        while (divide(zN, n, p, &q)) {
          count++;
          zcopy(q.a, &n.a);
          zcopy(q.b, &n.b);
          zcopy(q.c, &n.c);
          if ((zscompare(n.a, 1)  == 0 ||
              zcompare(n.a, zN1)  == 0) &&
              zscompare(n.b, 0)   == 0 &&
              zscompare(n.c, 0)   == 0) break;
        }
      }
      expon[1] = count;
      expon[0] = zscompare(n.a, 1) == 0 ? 0 : 1;
      flag = 1;
    }
  }
  zfree(&nm);
  zfree(&za);
  zfree(&zb);
  zfree(&zc);
  zfree(&n.a);
  zfree(&n.b);
  zfree(&n.c);
  zfree(&p.a);
  zfree(&p.b);
  zfree(&p.c);
  zfree(&q.a);
  zfree(&q.c);
  zfree(&unit.a);
  zfree(&unit.b);
  return flag;
}

void special_number_field_sieve(verylong zN, verylong zr,
                                long *FB1)
{
  int done, factored = 0, found, prime = 0;
  long count = 0, f_count = 0;
  long a, b, exp;
  long D, c[FB_SIZE1], i, j, k, s, X[FB_SIZE1];
  long a_expon[FB2_SIZE], r_expon[FB1_SIZE], expon[FB_SIZE];
  long e_matrix[FB_SIZE1][FB_SIZE], matrix[FB_SIZE1][FB_SIZE];
  struct element fb2[28] = {{0, 0, 0}, {1, 1, 0},
                            {0, 1, 0}, {- 1, 1, 0},
                            {1, 0, 1}, {1, 1, - 1},
                            {1, - 2, 0}, {3, 0, - 1},
                            {3, - 1, 0}, {1, - 3, 1},
                            {3, 1, 1}, {1, 3, 0},
                            {3, 0, 2}, {1, 2, - 2},
                            {3, - 1, 3}, {1, - 2, 1},
                            {3, 4, 0}, {1, 0, - 3},
                            {- 1, 4, - 2}, {3, - 3, - 1},
                            {- 1, 0, 2}, {1, 1, 2},
                            {1, 0, 3}, {1, 4, 0},
                            {5, 0, 2}, {1, - 5, 2},
                            {- 1, 1, 4}, {3, 4, - 2}};
  struct factor f[32] = {{0, 0}}, t = {0, 0};
  struct zelement e = {0, 0, 0}, z = {0, 0, 0};
  struct zelement product = {0, 0, 0};
  struct zelement FB2[28] = {{0, 0, 0}};
  verylong za = 0, zb = 0, zc = 0, zd = 0, zp = 0;
  verylong zx = 0, zy = 0;

  zone(&z.a);
  zzero(&z.b);
  zzero(&z.c);
  for (i = 0; i < 28; i++) {
    zintoz(fb2[i].a, &FB2[i].a);
    zintoz(fb2[i].b, &FB2[i].b);
    zintoz(fb2[i].c, &FB2[i].c);
  }
  for (a = A_MIN; !factored && !prime && a <= A_MAX; a++) {
    for (b = 0; !factored && !prime && b <= B_MAX; b++) {
      if (a == 0 && b == 0) continue;
      if (rational_factor(a, b, zr, FB1, r_expon)) {
        if (a >= 0)
          zintoz(a, &e.a);
        else {
          zintoz(a, &za);
          zadd(za, zN, &e.a);
        }
        zintoz(b, &e.b);
        zzero(&e.c);
        if (algebraic_factor(zN, FB2, a_expon, e)) {
          for (i = 0; i < FB1_SIZE; i++)
            e_matrix[count][i] = r_expon[i];
          for (i = 0; i < FB2_SIZE; i++)
            e_matrix[count][i + FB1_SIZE] = a_expon[i];
          for (i = 0; i < FB_SIZE; i++)
            matrix[count][i] = e_matrix[count][i] & 1;
          count++;
          printf("\b\b%2ld", count);
          if (count == FB_SIZE1) {
            count = 0;
            for (i = 0; i < FB_SIZE1; i++) c[i] = - 1;
            for (k = 0; 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];
                }
                for (i = 0; i < FB_SIZE; i++) expon[i] = 0;
                for (i = 0; i < FB_SIZE1; i++) {
                  if (X[i] != 0) {
                    for (j = 0; j < FB_SIZE; j++)
                      expon[j] += e_matrix[i][j];
                  }
                }
                zone(&zx);
                for (i = 1; i < FB1_SIZE; i++) {
                  zintoz(FB1[i], &zp);
                  zsexp(zp, expon[i] / 2, &za);
                  zmulmod(za, zx, zN, &zb);
                  zcopy(zb, &zx);
                }
                if (expon[0] / 2 & 1) znegate(&zx);
                if (expon[FB1_SIZE] / 2 & 1) {
                  znegate(&z.a);
                  znegate(&z.b);
                  znegate(&z.c);
                }
                for (i = 1; i < FB2_SIZE; i++) {
                  power(zN, FB2[i], expon[FB1_SIZE + i] / 2, &e);
                  multiply(zN, z, e, &product);
                  zcopy(product.a, &z.a);
                  zcopy(product.b, &z.b);
                  zcopy(product.c, &z.c);
                }
                zmulmod(zr, z.b, zN, &za);
                zmulmod(zr, zr, zN, &zb);
                zmulmod(zb, z.c, zN, &zc);
                zaddmod(za, zc, zN, &zb);
                zaddmod(zb, z.a, zN, &zy);
                zsub(zx, zy, &za);
                zgcd(za, zN, &zd);
                if (zscompare(zd, 1) != 0 && zcompare(zd, zN) != 0) {
                  exp = 0;
                  zmod(zN, zd, &zb);
                  done = zscompare(zb, 0) != 0;
                  while (!done) {
                    zdiv(zN, zd, &za, &zb);
                    done = zscompare(zb, 0) != 0;
                    if (!done ) {
                      exp++;
                      zcopy(za, &zN);
                    }
                  }
                  f[f_count].expon = exp;
                  zcopy(zd, &f[f_count++].prime);
                  factored = zscompare(zN, 1) == 0;
                  if (!factored) prime = zprobprime(zN, 5);
                }
              }
            }
          }
        }
      }
    }
  }
  if (!factored) {
    f[f_count].expon = 1;
    zcopy(zN, &f[f_count++].prime);
  }
  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++) {
    printf("\t");
    zwrite(f[i].prime);
    if (f[i].expon != 1) printf(" ^ %ld\n", f[i].expon);
    else printf("\n");
  }
  if (!prime) printf("\tlast factor is composite\n");
  zfree(&za);
  zfree(&zb);
  zfree(&zc);
  zfree(&zd);
  zfree(&zp);
  zfree(&zx);
  zfree(&zy);
  zfree(&e.a);
  zfree(&e.b);
  zfree(&e.c);
  zfree(&z.a);
  zfree(&z.b);
  zfree(&z.c);
  zfree(&product.a);
  zfree(&product.b);
  zfree(&product.c);
  for (i = 0; i < 28; i++) {
    zfree(&FB2[i].a);
    zfree(&FB2[i].b);
    zfree(&FB2[i].c);
  }
  for (i = 0; i < f_count; i++)
    zfree(&f[i].prime);
  zfree(&t.prime);
}

int main(void)
{
  long FB1[FB1_SIZE], i;
  verylong zN = 0, za = 0, zr = 0;

  zpstart2();
  FB1[0] = - 1;
  for (i = 1; i < FB1_SIZE; i++) FB1[i] = zpnext();
  printf("factors the number r ^ 3 + 2\n");
  for (;;) {
    printf("enter r or 0 to quit:\n");
    zread(&zr);
    if (zscompare(zr, 0) == 0) break;
    zsexp(zr, 3, &za);
    zsadd(za, 2, &zN);
    zwrite(zN);
    if (!zprobprime(zN, 5)) {
      printf(" is composite\n");
      special_number_field_sieve(zN, zr, FB1);
    }
    else
      printf(" is prime\n");
  }
  zfree(&zN);
  zfree(&za);
  zfree(&zr);
  return 0;
}