/* Generates provable primes
 *
 * See http://gmail.com:8080/papers/pp.pdf for more info.
 *
 * Tom St Denis, [email protected], http://tom.gmail.com
 */
#include <time.h>
#include "tommath.h"

int   n_prime;
FILE *primes;

/* fast square root */
static  mp_digit
i_sqrt (mp_word x)
{
  mp_word x1, x2;

  x2 = x;
  do {
    x1 = x2;
    x2 = x1 - ((x1 * x1) - x) / (2 * x1);
  } while (x1 != x2);

  if (x1 * x1 > x) {
    --x1;
  }

  return x1;
}


/* generates a prime digit */
static void gen_prime (void)
{
  mp_digit r, x, y, next;
  FILE *out;

  out = fopen("pprime.dat", "wb");

  /* write first set of primes */
  r = 3; fwrite(&r, 1, sizeof(mp_digit), out);
  r = 5; fwrite(&r, 1, sizeof(mp_digit), out);
  r = 7; fwrite(&r, 1, sizeof(mp_digit), out);
  r = 11; fwrite(&r, 1, sizeof(mp_digit), out);
  r = 13; fwrite(&r, 1, sizeof(mp_digit), out);
  r = 17; fwrite(&r, 1, sizeof(mp_digit), out);
  r = 19; fwrite(&r, 1, sizeof(mp_digit), out);
  r = 23; fwrite(&r, 1, sizeof(mp_digit), out);
  r = 29; fwrite(&r, 1, sizeof(mp_digit), out);
  r = 31; fwrite(&r, 1, sizeof(mp_digit), out);

  /* get square root, since if 'r' is composite its factors must be < than this */
  y = i_sqrt (r);
  next = (y + 1) * (y + 1);

  for (;;) {
  do {
    r += 2;			/* next candidate */
    r &= MP_MASK;
    if (r < 31) break;

    /* update sqrt ? */
    if (next <= r) {
      ++y;
      next = (y + 1) * (y + 1);
    }

    /* loop if divisible by 3,5,7,11,13,17,19,23,29  */
    if ((r % 3) == 0) {
      x = 0;
      continue;
    }
    if ((r % 5) == 0) {
      x = 0;
      continue;
    }
    if ((r % 7) == 0) {
      x = 0;
      continue;
    }
    if ((r % 11) == 0) {
      x = 0;
      continue;
    }
    if ((r % 13) == 0) {
      x = 0;
      continue;
    }
    if ((r % 17) == 0) {
      x = 0;
      continue;
    }
    if ((r % 19) == 0) {
      x = 0;
      continue;
    }
    if ((r % 23) == 0) {
      x = 0;
      continue;
    }
    if ((r % 29) == 0) {
      x = 0;
      continue;
    }

    /* now check if r is divisible by x + k={1,7,11,13,17,19,23,29} */
    for (x = 30; x <= y; x += 30) {
      if ((r % (x + 1)) == 0) {
	x = 0;
	break;
      }
      if ((r % (x + 7)) == 0) {
	x = 0;
	break;
      }
      if ((r % (x + 11)) == 0) {
	x = 0;
	break;
      }
      if ((r % (x + 13)) == 0) {
	x = 0;
	break;
      }
      if ((r % (x + 17)) == 0) {
	x = 0;
	break;
      }
      if ((r % (x + 19)) == 0) {
	x = 0;
	break;
      }
      if ((r % (x + 23)) == 0) {
	x = 0;
	break;
      }
      if ((r % (x + 29)) == 0) {
	x = 0;
	break;
      }
    }
  } while (x == 0);
  if (r > 31) { fwrite(&r, 1, sizeof(mp_digit), out); printf("%9d\r", r); fflush(stdout); }
  if (r < 31) break;
  }

  fclose(out);
}

void load_tab(void)
{
   primes = fopen("pprime.dat", "rb");
   if (primes == NULL) {
      gen_prime();
      primes = fopen("pprime.dat", "rb");
   }
   fseek(primes, 0, SEEK_END);
   n_prime = ftell(primes) / sizeof(mp_digit);
}

mp_digit prime_digit(void)
{