/*
  Author:  Pate Williams (c) 1997

  Exercise VI.4.3
  "For the following values of p and B, find (using
  a computer if necessary) the fraction of the
  integers between p + 1 - 2 sqrt(p) and
  p + 1 + 2 sqrt(p) which have no prime divisors
  greater than B: (a) p = 109, B = 3;
  (b) p = 109, B = 19; (c) p = 1009, B = 19;
  (d) p = 1009, B = 97; (e) p = 9973, B = 97."
  -Neal Koblitz-
  See "A Course in Number Theory and Cryptography"
  by Neal Koblitz second edition page 199.
*/

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

#define BITS_PER_LONG   32
#define BITS_PER_LONG_1 31
#define SIZE 32

struct factor {long expon, 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(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 factor(long n, long *sieve, struct factor *f,
            long *count)
/* factor using trial division */
{
  int one = 0;
  long e, p = 2;

  *count = 0;
  while (!one) {
    while (!get_bit(p, sieve)) p++;
    if (n % p == 0) {
      e = 0;
      do {
        e++;
        n = n / p;
      } while (n % p == 0);
      f[*count].expon = e;
      f[*count].prime = p;
      *count = *count + 1;
      one = n == 1;
    }
    p++;
  }
}

int main(void)
{
  long B[5] = {3, 19, 19, 97, 97};
  long p[5] = {109, 109, 1009, 1009, 9973};
  long Bi, a, b, c, count, i, m, n, pi, s;
  long sieve[1024] = {0};
  struct factor f[32];

  Sieve(15000, sieve);
  for (i = 0; i < 5; i++) {
    Bi = B[i];
    pi = p[i];
    s = 2 * sqrt(pi);
    a = pi + 1 - s;
    b = pi + 1 + s;
    n = b - a + 1;
    c = 0;
    for (m = a; m <= b; m++) {
      if (!get_bit(m, sieve)) {
        factor(m, sieve, f, &count);
        if (count == 1 && f[0].prime <= Bi) c++;
        else if (f[count - 1].prime <= Bi) c++;
      }
    }
    printf("%3ld out of %3ld\n", c, n);
  }
  return 0;
}