/*
  Author:  Pate Williams (c) 1997

  "Exercise 1.3. Computing pi(x).
  Incorporate your function phi(x, a) from Exercise
  1.2 above and the function Lehmer into a computer
  program which reads x and prints pi(x). Check your
  program by comparing its results with some entries
  from Table 3." -Hans Riesel-
  See "Prime Numbers and Computer Method for
  Factorization second edition page 23.
*/

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

#define A 5
#define G 100000l
#define MA 2310
#define MA2 1155
#define LENGTH 1155
#define PHIMA 480
#define BITS_PER_LONG   32
#define BITS_PER_LONG_1 31
#define SIZE 8192

long prime[A] = {2, 3, 5, 7, 11};
long sieve[SIZE] = {0}, table[SIZE] = {0};

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 init_table(void)
{
  long i, j, p, s = 0;

  for (i = 2; i <= LENGTH; i++)
    set_bit(i, i & 1, sieve);
  for (i = 0; i < A; i++) {
    p = prime[i];
    j = p;
    while (j <= LENGTH) {
      set_bit(j, 0, sieve);
      j += p;
    }
  }
  set_bit(1, 1, sieve);
  for (i = 1; i <= LENGTH; i++) {
    if (get_bit(i, sieve)) table[i] = ++s;
    else table[i] = s;
  }
  Sieve(G, sieve);
}

long phi5(long x)
{
  long absx = labs(x), r = absx % MA, z;

  z = (absx / MA) * PHIMA;
  if (r < MA2)
    z += table[r];
  else
    z += PHIMA - table[MA - r];
  return z * (absx / x);
}

long phi6(long x)
{
  return phi5(x) - phi5(x / 13);
}

long phi7(long x)
{
  return phi6(x) - phi6(x / 17);
}

long phi8(long x)
{
  return phi7(x) - phi7(x / 19);
}

long phi9(long x)
{
  return phi8(x) - phi8(x / 23);
}

long phi10(long x)
{
  return phi9(x) - phi9(x / 29);
}

long phi(long x, long a)
{
  if (a == 5) return phi5(x);
  if (a == 6) return phi6(x);
  if (a == 7) return phi7(x);
  if (a == 8) return phi8(x);
  if (a == 9) return phi9(x);
  return phi10(x);
}

long pi(long x)
{
  long i, s = 0;

  for (i = 2; i <= x; i++)
    if (get_bit(i, sieve)) s++;
  return s;
}

long i_prime(long i)
{
  long j = 1, p = 2;

  for (;;) {
    if (j == i) return p;
    j++, p++;
    while (!get_bit(p, sieve)) p++;
  }
}

long Lehmer(long x)
{
  long a, b, bi, c, i, j, sum, w, z;

  z = sqrt(x) + 0.5;
  b = pi(z);
  c = pi(exp(log(x) / 3.0) + 0.5);
  a = pi(sqrt(z) + 0.5);
  sum = phi(x, a) + (b + a - 2) * (b - a + 1) / 2;
  for (i = a + 1; i <= b; i++) {
    w = x / i_prime(i);
    sum -= pi(w);
    if (i <= c) {
      bi = pi(sqrt(w) + 0.5);
      for (j = i; j <= bi; j++)
        sum = sum - pi(w / i_prime(j)) + j - 1;
    }
  }
  return sum;
}

int main(void)
{
  long x;

  init_table();
  printf("x = ");
  scanf("%ld", &x);
  printf("pi(%ld) = %ld\n", x, Lehmer(x));
  return 0;
}