// Calculate number of primes pi(2^n) for n=2..32 (see A007053)
// Calculate max prime <= 2^n for n=2..32 (see A014234)

// C++ source. Compiles happily under g++. Assumes int >= 32 bits.

#include <stdio.h>

#define sieve_log2_bits_per_word 5
#define sieve_root_bits 16
#define sieve_bits (sieve_root_bits * 2)
#define sieve_square_root (1 << sieve_root_bits)
#define sieve_word_index(odd_x) ((odd_x) >> (sieve_log2_bits_per_word + 1))
#define sieve_bit_mask(odd_x) (1 << (((odd_x) >> 1) & ((1 << sieve_log2_bits_per_word) - 1)))
#define sieve_test(odd_x) (sieve_words[sieve_word_index(odd_x)] & sieve_bit_mask(odd_x))
#define sieve_clear(odd_x) (sieve_words[sieve_word_index(odd_x)] &= ~sieve_bit_mask(odd_x))

unsigned int sieve_words[1 << (sieve_bits - sieve_log2_bits_per_word - 1)];

int main(int argc, char *argv[]) {
  for (int i = sizeof(sieve_words) / sizeof(sieve_words[0]); --i >= 0; ) {
    sieve_words[i] = (unsigned int) -1;
  }

  sieve_clear(1);
  for (int i = 3; i < sieve_square_root; i += 2) {
    if (((i + 1) & i) == 0) {
      int log2i = 0;
      for (int j = i; j > 0; j>>=1) { log2i++; }
      printf("Initializing sieve -- 2^%d\n", log2i);
    }
    if (sieve_test(i)) {
      unsigned int delta = i * 2;
      for (unsigned int j = i + delta; j > delta; j += delta) {
        sieve_clear(j);
      }
    }
  }
  int prime_pi = 1; // 2
  unsigned int i;
  unsigned int last_p;
  for (i = 3; i > 1; i += 2) {
    if (sieve_test(i)) {
      prime_pi++;
      last_p = i;
    }
    if (((i + 1) & i) == 0) {
      int log2i = 0;
      for (unsigned int j = i; j > 0; j>>=1) { log2i++; }
      printf("%d primes < 2^%d -- last prime %u\n", prime_pi, log2i, last_p);
    }
  }
  return 0;
}
/* Fred Curtis (fred(AT)f2 org), Dec 08 2009*/