Bug
At the end of your program, the b array is all full of 1s, which means that you didn't find any primes. The problem is here:
num_t j = i;
Because you start your j loop at i, you will mark i (which is prime) as non-prime. You should start your j loop at i*i instead (see below).
Don't get tricky
Why write this:
for(num_t i = 0; i ^ N; ++i)
when you can write this:
for(num_t i = 0; i < N; ++i)
I had to stare at your code for a long time just to figure out whether it was correct, for a simple loop!
Also, you did it here:
for(num_t i = 3; i ^ 1 ^ N; i += 2) {
which should be:
for(num_t i = 3; i < N; i += 2) {
You can afford to waste 2 bytes
Here, you use index-2 to save 2 bytes:
if(!b[i - 2]) {
// printf("%llu\n", i);
const num_t increment = i << 1;
num_t j = i;
while(j < N) {
b[j - 2] = 1;
j += increment;
}
}
There's no need to save 2 bytes. Just use b[index] instead of b[index-2]. You already allocated the full array anyways.
Speed improvements
There are three things I can see that could speed up your program:
- In your
j loop, you can start j at i*i instead of just i. All non-primes less than i*i will already be marked.
- You can stop the
i loop at sqrt(N) instead of at N, because once you reach sqrt(N), the second loop won't mark any more entries.
Since you are finding primes up to 10^9, you are using 1GB of memory for the b array. You should get better performance if you use less memory, because the processor cache will be utilized more effectively if your array is smaller.
a) You can only track odd numbers, which will reduce your memory usage to 1/2 of the original (512MB).
b) You can use 1 bit per number instead of 1 byte per number. This will reduce your memory to 1/8 of the original. Combined with (a), it will be 1/16 the original, or 64MB. This is a much better than 1GB.
Example using 1 bit per odd number
#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#define N 1000000000
int main(void)
{
int arraySize = ((N + 63)/64 + 1) * sizeof(uint32_t);
uint32_t *primes = malloc(arraySize);
// Each bit in primes is an odd number.
// Bit 0 = 1, bit 1 = 3, bit 2 = 5, etc.
memset(primes, 0xff, arraySize);
// 1 is not a prime.
primes[0] &= ~0x1;
int sqrt_N = sqrt(N);
for(int i = 3; i <= sqrt_N; i += 2) {
int iIndex = i >> 6;
int iBit = (1 << ((i >> 1) & 31));
if ((primes[iIndex] & iBit) != 0) {
int increment = i+i;
for (int j = i * i; j < N; j += increment) {
int jIndex = j >> 6;
int jBit = (1 << ((j >> 1) & 31));
primes[jIndex] &= ~jBit;
}
}
}
// Count the number of primes in order to verify that the above worked.
// Start count at 1 to include 2 as the first prime, since we are only
// going to count odd primes.
int count = 1;
for (int i = 3; i < N; i += 2) {
int iIndex = i >> 6;
int iBit = (1 << ((i >> 1) & 31));
if (primes[iIndex] & iBit)
count++;
}
printf("%d\n", count);
free(primes);
return 0;
}
Edit: Even less memory usage
Pete Kirkham suggested using even less memory by using 2 out of every 6 numbers. In other words, using 1 bit per every 3 numbers instead of 1 bit per every 2 numbers. At first I was skeptical because this required using a division in the inner loop. However, after coding it up, it turned out to be faster. The code is quite a bit trickier however, because the inner loops need to avoid any multiples of 3, because all multiples of 3 are no longer in the primes array:
#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#define N 1000000000
int main(void)
{
int arraySize = (N/24 + 1);
uint32_t *primes = malloc(arraySize);
// The bits in primes follow this pattern:
//
// Bit 0 = 5, bit 1 = 7, bit 2 = 11, bit 3 = 13, bit 4 = 17, etc.
//
// For even bits, bit n represents 5 + 6*n
// For odd bits, bit n represents 1 + 6*n
memset(primes , 0xff, arraySize);
int sqrt_N = sqrt(N);
for(int i = 5; i <= sqrt_N; i += 4) {
int iBitNumber = i / 3 - 1;
int iIndex = iBitNumber >> 5;
int iBit = 1 << (iBitNumber & 31);
if ((primes[iIndex] & iBit) != 0) {
int increment = i+i;
for (int j = i * i; j < N; j += increment) {
int jBitNumber = j / 3 - 1;
int jIndex = jBitNumber >> 5;
int jBit = 1 << (jBitNumber & 31);
primes[jIndex] &= ~jBit;
j += increment;
if (j >= N)
break;
jBitNumber = j / 3 - 1;
jIndex = jBitNumber >> 5;
jBit = 1 << (jBitNumber & 31);
primes[jIndex] &= ~jBit;
// Skip multiple of 3.
j += increment;
}
}
i += 2;
iBit <<= 1;
if ((primes[iIndex] & iBit) != 0) {
int increment = i+i;
for (int j = i * i; j < N; j += increment) {
int jBitNumber = j / 3 - 1;
int jIndex = jBitNumber >> 5;
int jBit = 1 << (jBitNumber & 31);
primes[jIndex] &= ~jBit;
// Skip multiple of 3.
j += increment;
j += increment;
if (j >= N)
break;
jBitNumber = j / 3 - 1;
jIndex = jBitNumber >> 5;
jBit = 1 << (jBitNumber & 31);
primes[jIndex] &= ~jBit;
}
}
}
// Initial count includes 2, 3.
int count=2;
for (int i=5;i<N;i+=6) {
int iBitNumber = i / 3 - 1;
int iIndex = iBitNumber >> 5;
int iBit = 1 << (iBitNumber & 31);
if (primes[iIndex] & iBit) {
count++;
}
iBit <<= 1;
if (primes[iIndex] & iBit) {
count++;
}
}
printf("%d\n", count);
free(primes);
return 0;
}
I'm sure the variant that works on multiples of 30 would be even faster than this, although the code would be even more complicated.
Timings
Here are the speeds of the various programs I ran on my computer:
Original program : 13.33 seconds
OP fixed (i*i, sqrt(n)) : 7.75 seconds
1 bit per odd number : 3.90 seconds
1 bit per 3 numbers : 2.77 seconds
Note that I did not include the time it took to count the primes to verify the result. I modified each program to count the number of primes just to make sure that they all worked, which is how I discovered that the original program had a bug in it.
