// 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 #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*/