A simple implementation of the Sieve of Eratosthenes, a simple, ancient algorithm for finding all prime numbers up to a specified integer. More information can be found on wikipedia, here.
We use a bit array to efficiently store the state of each number. If a number is 0, its a prime. If a number is 1, it isn't. Simple, no? As an example, after the sieve runs, prime[0] will look like this:
| Bit Index | 7 | 6 | 5 | 4 | 3 | 2 | 1 | 0 |
|---|---|---|---|---|---|---|---|---|
| Value | 0 | 1 | 0 | 1 | 0 | 0 | 1 | 1 |
| Decimal | 7 | 6 | 5 | 4 | 3 | 2 | 1 | 0 |
| Meaning | Prime | Non-prime | Prime | Non-prime | Prime | Prime | Non-prime | Non-prime |
The next byte, prime[1], would look like this:
| Bit Index | 7 | 6 | 5 | 4 | 3 | 2 | 1 | 0 |
|---|---|---|---|---|---|---|---|---|
| Value | 0 | 0 | 0 | 1 | 0 | 1 | 1 | 1 |
| Decimal | 15 | 14 | 13 | 12 | 11 | 10 | 9 | 8 |
| Meaning | Non-prime | Non-prime | Prime | Non-prime | Prime | Non-prime | Non-prime | Non-prime |
Benchmarks are done with -O3, average is over 400 iterations, and no output. Printing takes another ~8s at least.
| Language | CPU Speed | Average for 10 000 000 |
|---|---|---|
| C | 1 Ghz | 0.175s |
This sieve will do up to 1 000 000 000 in under 34 seconds on a decent machine.
#include <stdlib.h> #include <stdio.h> #include <math.h> #define NumberToFind 10000000 // The maximum number to find #define LocalAlloc 1 // Allocate storage on the stack #define SetBit(d,c) ((d) |= (1L << (c))) // Set a bit to 1 #define TestBit(d,c) ((d) & (1L << (c))) // Check if a bit is set #if !LocalAlloc unsigned char *prime = NULL; #endif int main(void) { #if LocalAlloc unsigned char prime[(NumberToFind / 8)+1] = {0}; #else prime = (unsigned char*)malloc((NumberToFind / 8)+1); #endif int x,y,t = NumberToFind; prime[0] |= 3; // Set the first two bits to 1 (00000011) - 0 and 1 are both non-prime. for(x = 2; x < (int)(sqrt(NumberToFind)+1); x++) { if(TestBit(prime[x/8],((x % 8))) == 0) { for(y = x*x; y < NumberToFind; y+=x) SetBit(prime[y/8],(y % 8)); } } // Print out the results for(x = 0; x < NumberToFind+1; x++) { if(TestBit(prime[x/8],(x % 8)) == 0) printf("%d\n",x); } // Wait for the enter key to exit fgetc(stdin); #if !LocalAlloc free(prime); #endif return 0; }