/*
  Author:  Pate Williams (c) 1997

  Aldleman-Pomerance-Rumely-Cohen-Lenstra primality
  test. See "Implementation of a New Primality Test",
  Mathematics of Computation, Volume 48, Number 177,
  January 1987, pages 103 - 121 by H. Cohen and
  A. K. Lenstra. Table generation.
*/

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

#define t 25200l

struct power {long x, p;};
struct long_factor {long expon, prime;};
struct very_factor {int expon; verylong prime;};

void long_trial_division(long bound, verylong zn,
                         struct long_factor *f,
                         long *count)
{
  int one = 0, pri = 0, value;
  long e, p;
  verylong za = 0, zb = 0;

  *count = 0;
  zcopy(zn, &za);
  zpstart2();
  do {
    p = zpnext();
    if (zsmod(za, p) == 0) {
      e = 0;
      do {
        e++;
        zsdiv(za, p, &zb);
        zcopy(zb, &za);
      } while (zsmod(za, p) == 0);
      f[*count].expon = e;
      f[*count].prime = p;
      *count = *count + 1;
      one = zscompare(za, 1l) == 0;
      if (!one) pri = zprobprime(za, 5l);
    }
  } while (!one && !pri && p <= bound);
  if (pri) {
    f[*count].expon = 1;
    f[*count].prime = ztoint(za);
    *count = *count + 1;
    value = zscompare(za, LONG_MAX) > 0;
  }
  if (value || p > bound) {
    fprintf(stderr, "fatal error\ncan't factor number");
    exit(1);
  }
  zfree(&za);
  zfree(&zb);
}

void compute_e_of_t(verylong zt, long *Q,
                    verylong *ze, long *Q_count)
{
  int found;
  long e, i, lf_count, p, q;
  struct long_factor lf[16];
  verylong zq = 0, zr = 0;

  *Q_count = 0;
  long_trial_division(1000000l, zt, lf, &lf_count);
  printf("t = "); zwriteln(zt);
  printf("prime factorization of t\n");
  for (i = 0; i < lf_count; i++) {
    printf("%2ld", lf[i].prime);
    if (lf[i].expon == 1) printf("\n");
    else printf(" ^ %ld\n", lf[i].expon);
  }
  zintoz(2l, ze);
  zpstart2();
  do {
    p = zpnext();
    q = p - 1;
    if (zsmod(zt, q) == 0) {
      Q[*Q_count] = p;
      *Q_count = *Q_count + 1;
      e = 1;
      found = 0;
      for (i = 0; !found && i < lf_count; i++) {
        if (p == lf[i].prime) {
          found = 1;
          e += lf[i].expon;
        }
      }
      zintoz(p, &zq);
      zsexp(zq, e, &zr);
      zmul(*ze, zr, &zq);
      zcopy(zq, ze);
    }
  } while (zscompare(zt, q) >= 0);
  zfree(&zq);
  zfree(&zr);
}

void very_trial_division(long bound, verylong zn,
                         struct very_factor *f,
                         long *count)
{
  int e, one = 0, pri = 0;
  long p;
  verylong za = 0, zb = 0;

  *count = 0;
  zcopy(zn, &za);
  zpstart2();
  do {
    p = zpnext();
    if (zsmod(za, p) == 0) {
      e = 0;
      do {
        e++;
        zsdiv(za, p, &zb);
        zcopy(zb, &za);
      } while (zsmod(za, p) == 0);
      f[*count].expon = e;
      zintoz(p, &f[*count].prime);
      *count = *count + 1;
      one = zscompare(za, 1l) == 0;
      if (!one) pri = zprobprime(za, 5l);
    }
  } while (!one && !pri && p <= bound);
  if (pri) {
    f[*count].expon = 1;
    zcopy(za, &f[*count].prime);
    *count = *count + 1;
  }
  if (p > bound) {
    fprintf(stderr, "fatal error\ncan't factor number");
    exit(1);
  }
  zfree(&za);
  zfree(&zb);
}

void primitive_root(verylong zp, verylong *za)
{
  long i, k;
  struct very_factor f[16];
  verylong zb = 0, ze = 0, zq = 0, zr = 0, zp1 = 0;

  for (i = 0; i < 16; i++) f[i].prime = 0;
  zsadd(zp, - 1l, &zp1);
  very_trial_division(1000000l, zp1, f, &k);
  zone(za);
  L2:
  i = 0;
  zsadd(*za, 1l, &zb);
  zcopy(zb, za);
  L3:
  zdiv(zp1, f[i].prime, &zq, &zr);
  zexpmod(*za, zq, zp, &ze);
  if (zscompare(ze, 1l) == 0) goto L2;
  if (++i < k) goto L3;
  zfree(&zb);
  zfree(&ze);
  zfree(&zq);
  zfree(&zr);
  zfree(&zp1);
  for (i = 0; i < 16; i++) zfree(&f[i].prime);
}

int fcmp(const void *g1, const void *g2)
{
  struct power *p1 = (struct power *) g1;
  struct power *p2 = (struct power *) g2;

  if (p1->p < p2->p) return - 1;
  if (p1->p > p2->p) return + 1;
  return 0;
}

void build_f_table(long q, verylong zg, verylong zq, FILE *file)
{
  long f, q1 = q - 1, q2 = q - 2, x;
  long *h = calloc(q, sizeof(long));
  struct power *g, key, *ptr;
  verylong zr = 0, zs = 0, zt = 0;

  g = calloc(q, sizeof(struct power));
  if (!g || !h) {
    fprintf(stderr, "fatal error\ninsufficient memory\n");
    exit(1);
  }
  zone(&zr);
  g[0].x = 0;
  g[0].p = 1;
  for (x = 1; x <= q2; x++) {
    zsexpmod(zg, x, zq, &zs);
    zsubmod(zr, zs, zq, &zt);
    g[x].x = x;
    g[x].p = ztoint(zs);
    h[x] = ztoint(zt);
  }
  qsort(g, q1, sizeof(struct power), fcmp);
  fwrite(&q, sizeof(q), 1, file);
  for (x = 1; x <= q2; x++) {
    key.x = x;
    key.p = h[x];
    ptr = bsearch(&key, g, q1, sizeof(struct power), fcmp);
    if (!ptr) {
      fprintf(stderr, "fatal error\ncan't find logarith\n");
      exit(1);
    }
    f = ptr->x;
    fwrite(&f, sizeof(f), 1, file);
  }
  free(g);
  free(h);
  zfree(&zr);
  zfree(&zs);
  zfree(&zt);
}

void generate_f_table(long q_count, long *q)
{
  long i, i0, p;
  FILE *file = fopen("f_table.dat", "w+b");
  verylong zg = 0, zp = 0;

  if (!file) {
    fprintf(stderr, "fatal error\ncan't open file\n");
    exit(1);
  }
  if (q[0] == 2) i0 = 1; else i0 = 0;
  for (i = i0; i < q_count; i++) {
    printf("\b\b%ld", i);
    p = q[i];
    zintoz(p, &zp);
    primitive_root(zp, &zg);
    build_f_table(p, zg, zp, file);
  }
  zfree(&zg);
  zfree(&zp);
  fclose(file);
}

void generate_j_table(long q_count, long q0)
{
  long A, B, *a, *f, i, k, l, m, p, pk;
  long pk1, pk3, q, q1, r, r0, x;
  FILE *f_file = fopen("f_table.dat", "r+b");
  FILE *j_file = fopen("j_table.dat", "w+b");

  if (!f_file || !j_file) {
    fprintf(stderr, "fatal error\ncan't open file(s)\n");
    exit(1);
  }
  if (q0 == 2) r0 = 1; else r0 = 0;
  for (r = r0; r < q_count; r++) {
    printf("\b\b%ld", r);
    fread(&q, sizeof(q), 1, f_file);
    f = calloc(q, sizeof(long));
    if (!f) {
      fprintf(stderr, "fatal error\ninsufficient memory\n");
      exit(1);
    }
    for (i = 1; i <= q - 2; i++)
      fread(&f[i], sizeof(f[i]), 1, f_file);
    q1 = q - 1;
    zpstart2();
    do {
      p = zpnext();
      if (q1 % p == 0) {
        k = 0;
        do {
          k++;
          q1 /= p;
        } while (q1 % p == 0);
      }
      pk = pow(p, k);
      pk1 = pow(p, k - 1);
      m = (p - 1) * pk1;
      a = calloc(m, sizeof(long));
      if (!a) {
        fprintf(stderr, "fatal error\ninsufficient memory\n");
        exit(1);
      }
      if (pk != 2) {
        for (i = 0; i < m; i++) a[i] = 0;
        A = B = 1;
        for (x = 1; x <= q - 2; x++) {
          l = (A * x + B * f[x]) % pk;
          if (l < m) a[l]++;
          else
            for (i = 1; i < p; i++)
              a[((l - i) * pk1) % l]--;
        }
        fwrite(&m, sizeof(m), 1, j_file);
        for (i = 0; i < m; i++)
          fwrite(&a[i], sizeof(a[i]), 1, j_file);
      }
      else if (p == 2 && k >= 3) {
        for (i = 0; i < m; i++) a[i] = 0;
        A = 2, B = 1;
        for (x = 1; x <= q - 2; x++) {
          l = (A * x + B * f[x]) % pk;
          if (l < m) a[l]++;
          else
            for (i = 1; i < p; i++)
              a[((l - i) * pk1) % l]--;
        }
        fwrite(&m, sizeof(m), 1, j_file);
        for (i = 0; i < m; i++)
          fwrite(&a[i], sizeof(a[i]), 1, j_file);
        for (i = 0; i < m; i++) a[i] = 0;
        pk3 = pow(2, k - 3);
        A = 3 * pk3, B = pk3;
        for (x = 1; x <= q - 2; x++) {
          l = (A * x + B * f[x]) % pk;
          if (l < m) a[l]++;
          else
            for (i = 1; i < p; i++)
              a[((l - i) * pk1) % l]--;
        }
        fwrite(&m, sizeof(m), 1, j_file);
        for (i = 0; i < m; i++)
          fwrite(&a[i], sizeof(a[i]), 1, j_file);
      }
      free(a);
    } while (q1 != 1);
  }
  free(f);
  fclose(f_file);
  fclose(j_file);
}

int main(void)
{
  long q_count, q[64];
  verylong ze = 0, zt = 0;

  zintoz(t, &zt);
  compute_e_of_t(zt, q, &ze, &q_count);
  printf("e(t) = "); zwriteln(ze);
  printf("log10(e(t)) = %lf\n", zln(ze) / log(10.0));
  printf("number of q primes = %ld\n", q_count);
  generate_f_table(q_count, q);
  printf("\ngenerated f_table\n");
  generate_j_table(q_count, q[0]);
  printf("\ngenerated j_table\n");
  zfree(&ze);
  zfree(&zt);
  return 0;
}