Skip to content

Instantly share code, notes, and snippets.

@hjr3
Created March 12, 2012 06:26
Show Gist options
  • Save hjr3/2020292 to your computer and use it in GitHub Desktop.
Save hjr3/2020292 to your computer and use it in GitHub Desktop.
Find all prime numbers in a given range using Sieve of Eratosthenes
// gcc -std=c99 -Wall -Werror -g -o primes primes.c -lm
#include <stdio.h>
#include <string.h>
#include <limits.h>
#include <math.h>
void primes(int length)
{
unsigned int bits_per_int = (sizeof(unsigned int) * CHAR_BIT);
unsigned int mask_size = 0;
for (unsigned int x = bits_per_int; x > 1; x >>= 1) {
mask_size++;
}
bits_per_int--;
unsigned int arr_length;
unsigned int max;
arr_length = (length / bits_per_int) + 1;
unsigned int primes[arr_length];
unsigned int index;
unsigned int bit;
unsigned int n;
memset(primes, 0, arr_length);
// only need to go sqrt of the array length
// since all divisors larger than the sqrt will
// have already been processed
max = sqrt(length);
for (int i = 2; i <= max; i++) {
n = i * i;
index = i >> mask_size;
bit = 1 << (i & bits_per_int);
if ((primes[index] & bit) == 0) {
while (n < length) {
n += i;
index = n >> mask_size;
bit = 1 << (n & bits_per_int);
primes[index] |= bit;
}
}
}
for (int i = 0; i < length; i++) {
index = i >> mask_size;
bit = 1 << (i & bits_per_int);
if ((primes[index] & bit) == 0) {
printf("%d\n", i);
}
}
}
int main(int argc, char *argv[])
{
primes(120);
return 0;
}
@charlesgarvin
Copy link

Yeah, that's the sieve. A few issues I noticed and places to improve:

  1. Bit shifting is not well defined on signed types, always use unsigned.
  2. You should avoid magic numbers like 5 and 31. You can replace 5 with a small loop (see below). You can calculate how many bits in an integer with bits_per_int = sizeof(unsigned int) * CHAR_BIT; CHAR_BIT is in limits.h and tells you how many bits per byte. This makes your code work with 32 and 64 bit machines.
  3. Your array is much bigger than it needs to be. sizeof gives you number of bytes, not bits, so you're getting (assuming 32-bit) length/4 instead of length/bits_per_int.
  4. Use memset to set primes to all-zero-bits, it's probably quite a bit quicker than your loop.
  5. You don't need an if/else for setting the index. index >> 5 works fine if n < 31 (which should probably be < 32 anyway). The bits will just shift off the end and you'll get 0.
  6. 1 is not a prime number.
  7. No need for argc and argv if you're not using them, and it's good practice to actually return an int at the end of main.

// something like this should do it
unsigned int mask_size;
unsigned int x;
for (x = bits_per_int, mask_size = 0; x > 1; x >>= 1) {
mask_size++;
}

mask_size will correctly be 5 for a 32-bit system, 6 for a 64-bit system, etc.

@hjr3
Copy link
Author

hjr3 commented Mar 12, 2012 via email

@charlesgarvin
Copy link

Could be. Don't know what your weird results were and what exactly you were trying to do with (n < 31) >> 5. It's common to do "sign extension" or "sign preservation" when shifting negative numbers, so the effect of the shift is equivalent to a multiply/divide by 2 (as it is for positive numbers).

A couple of optimizations:

  1. You can start n = i * i; since all lower multiples of i will have been crossed off by previous values of i. E.g., no need to start at 6 when crossing off multiples of 3, since 2 already took care of it.
  2. You can stop your outer loop at sqrt(length) instead of half (but make sure the loop actually goes up to and includes square_root), since all divisors larger than square_root will have been crossed off by the smaller divisor pair. E.g. if length is 1000, the largest number you need to cross off multiples of is 31.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment