/*
  Author:  Pate Williams (c) 1997

  "Exercise 7.1. Sieving for a prime in
  (10 ^ 60, 10 ^ 60 + 10000)." -Hans Riesel-
  See "Prime Numbers and Computer Methods
  for Factorization" by Hans Riesel second
  edition pages 234-235.
*/

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

#define BITS_PER_LONG   32
#define BITS_PER_LONG_1 31
#define SIEVE_LIMIT 10000
#define FACTOR_LIMIT 100000l
#define SIEVE_SIZE 8192

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(0, 0, sieve);
  set_bit(1, 0, sieve);
  set_bit(2, 1, 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, 0, sieve);
      i += inc;
    }
    c += 2;
    while (!get_bit(c, sieve)) c++;
  } while (c * c <= n);
}

void Sieve_Interval(verylong zm, verylong zn, long *sieve)
{
  long b, base[SIEVE_SIZE], c, i, inc, l, s;
  verylong za = 0, zd = 0, zi = 0, zs = 0;

  zsqrt(zn, &zs, &zd);
  if (zscompare(zs, SIEVE_LIMIT) <= 0)
    s = ztoint(zs);
  else s = 10000;
  Sieve(s, base);
  zsub(zn, zm, &zd);
  zsadd(zd, 1, &zi);
  l = ztoint(zi);
  zcopy(zm, &zi);
  for (i = 0; i < l; i++) {
    b = zodd(zi) ? 1 : 0;
    set_bit(i, b, sieve);
    zsadd(zi, 1, &zd);
    zcopy(zd, &zi);
  }
  c = 3;
  do {
    if (zsdiv(zm, c, &za) == 0)
      zsmul(za, c, &zi);
    else {
      zsmul(za, c, &zd);
      zsadd(zd, c, &zi);
    }
    if (zscompare(zi, c) == 0) {
      zsadd(zi, c, &zd);
      zcopy(zd, &zi);
    }
    inc = c;
    while (zcompare(zi, zn) <= 0) {
      zsub(zi, zm, &zd);
      set_bit(ztoint(zd), 0, sieve);
      zsadd(zi, inc, &zd);
      zcopy(zd, &zi);
    }
    c += 2;
    while (!get_bit(c, base)) c++;
  } while (c * c <= SIEVE_LIMIT);
  zfree(&za);
  zfree(&zd);
  zfree(&zi);
  zfree(&zs);
}

int factor(verylong zn)
{
  long p;
  verylong za = 0, zb = 0;

  zcopy(zn, &za);
  zpstart2();
  p = zpnext();
  while (p <= FACTOR_LIMIT) {
    if (zsdiv(za, p, &zb) == 0) {
      zcopy(zb, &za);
      while (zsdiv(za, p, &zb) == 0)
        zcopy(zb, &za);
      if (zscompare(za, 1) == 0) {
        zfree(&za);
        zfree(&zb);
        return 1;
      }
    }
    p = zpnext();
  }
  zfree(&za);
  zfree(&zb);
  return 0;
}

int main(void)
{
  long b = 10, count = 0, d = 10000, e = 60, i;
  long sieve[SIEVE_SIZE];
  verylong za = 0, zb = 0, zm = 0, zn = 0, zp = 0;

  zintoz(b, &zb);
  zsexp(zb, e, &za);
  zsadd(za, 1, &zm);
  zsadd(zm, d, &zn);
  Sieve_Interval(zm, zn, sieve);
  for (i = 0; i < d + 1; i++)
    if (get_bit(i, sieve)) count++;
  printf("number of pseudoprimes = %ld\n", count);
  count = 0;
  for (i = 0; count < 15 && i < d + 1; i++) {
    if (get_bit(i, sieve)) {
      zsadd(zm, i, &zp);
      if (zprobprime(zp, 10)) {
        zsadd(zp, - 1, &za);
        zsadd(zp, + 1, &zb);
        if (!factor(za) && !factor(zb)) {
          printf("10 ^ 60 + %4ld\n", i + 1);
          count++;
        }
      }
    }
  }
  zfree(&za);
  zfree(&zb);
  zfree(&zm);
  zfree(&zn);
  zfree(&zp);
  return 0;
}