# Sieve of Sundaram

Every integer can be written in the form $2^i*v$ where $v$ is odd. This algorithm generates a list of numbers that are only odd. We want to ensure that all the composite odd numbers are marked out.

In the final list, all numbers will be of the form $2^0 v$. So, we only need to worry about odd factors. An odd number will be composite if it is of the form $(2i + 1)(2j + 1)$, where $i$ and $j$ are strictly positive integers. So, a composite odd number $q = 2k + 1 = 4ij + 2i + 2j + 1$. This means $k = i + j + 2ij$. Whenever $k$ is of the form, $i + j + 2ij, 2k+ 1$ is composite and whenever an odd integer $q$ is composite, the condition holds for $k, q = 2k + 1$.

All we seek to do is mark all numbers that can be written as $i +j + 2ij$. And for any unmarked integer $k$, we know that $2k + 1$ will be prime. Note - We need to hardcode 2 since this algorithm only generates the odd primes.

#include <stdio.h>
#include <stdlib.h>
#define target 50000

//Function prototypes
void initial_marking(long int[]);
void marking(long int[]);
void make_prime_list(long int[],long int[]);
void print_prime(long int[]);

int main()
{
//Half of the numbers from 1 to target will be even so the number of primes will never exceed target/2.
//The algorithm only worries about the first half of the range. That's why list also goes till target/2.
long list[target/2], primes[target/2];

//Sieve of Sundaram
initial_marking(list);
//Crossing out all numbers of the form i + j + 2ij
marking(list);

//Making the list of primes
make_prime_list(primes, list);

//Displaying the list
print_primes(primes);

return 0;
}

//Initially, everything is marked 1
void initial_marking(long int list[])
{
int i;

for(i = 1; i <= target/2; i++)
{
list[i] = 1;
}
}

//Crossing out numbers of the form i + j + 2ij
void marking(long int list[])
{
int i, j, crossed_out_num;

//All numbers of the form i + j + 2ij need to be marked out
for(i = 1; i <= target/2; i++)
{
for(j = 1; j <= i; j++)
{
crossed_out_num = i + j + 2*i*j;

if(crossed_out_num <= target/2)//If this condition isn't there, we might me marking numbers not in list.
{
list[crossed_out_num] = 0;
}
}
}
}

//Making the list of primes
void make_prime_list(long int primes[],long int list[])
{
int i = 0, primeCount = 0;

//We need to put 2 in ourselves because the algorithm 'only' generates all the odd primes
primes[primeCount++] = 2;
for(i = 1; i <= target/2; i++)
{
//Checking if the number is not crossed out and adding 2i + 1 to the list if it is unmarked.
if(list[i] == 1)
{
primes[primeCount++] = (2*i + 1);
}
}

//Putting 0 to mark the end of the array
primes[primeCount] = 0;
}

void print_primes(long int primes[])
{
int i;

//0 indicates the end of the array
for(i = 0; primes[i] != 0; i++)
{
printf("%ld\t",primes[i]);
}
}


Let $f(i,j) = i + j + 2ij$

Then $f(i,j+1) = i + j + 1 + 2i(j+1) = i+j+2ij+2i+1 = f(i,j)+2i+1$

Using this, the marking() function (which is the highest complexity function) can be rewritten as such:

void marking(long int list[])
{
long long int i, crossed_out_num, i_limit, cross_limit, increment;
i_limit = target/6; // i > i_limit will never give f(i,j) in range

//All numbers of the form i + j + 2ij need to be marked out
for(i = 1; i <= i_limit; i++)
{
increment = 2*i + 1;

// f(i,j) <= target/2 and j <= i
cross_limit = target/2;
if(cross_limit > 2*i*(i+1)) cross_limit = 2*i*(i+1);

// Find crossed_out_num directly by starting with j=1 and using the above result
for(crossed_out_num = 3*i + 1; crossed_out_num <= cross_limit; crossed_out_num += increment)
{
list[crossed_out_num] = 0;
}
}
}


This should speed up the code considerably and avoid overflow when i becomes large.

• Can you explain how you arrived at target/6 ? I also have another question. My program seems to work till limit is 60 thousand but doesn't work when I set limit greater than that. Is it because of time complexity or because of something else ? Like array being too large or something Commented Feb 3, 2017 at 15:40
• For a given value of i, the smallest value of f(i,j) is 3i+1. We want it to be less than target/2, hence the limit of imax=target/6. By "doesn't work" do you mean wrong answer or crash? If wrong answer that might be because of i*j overflow. Crash would be because you are allocating array memory in main() so it goes to stack. Either move the arrays global or use heap-based allocation using malloc(). Commented Feb 3, 2017 at 15:45
• I'll try moving it to main. Thanks. Are you a mathematician ? Commented Feb 3, 2017 at 16:06
• Can you explain why declaring it in main would cause a problem Commented Feb 3, 2017 at 16:08
• Putting it in global doesn't work. Advise me about the heap based allocation. Commented Feb 3, 2017 at 16:09

    long list[target/2], primes[target/2];


Elsewhere it's long int, here it's long. Be consistent.

list only ever contains values 0 and 1. Ideally it would use a 1-bit integer type (and you could consider writing a basic bitset implementation if it's not in the library). It's certainly wasteful to use long int.

Integer types in C are a mess, which is why <stdint.h> was created. This might be controversial, but my opinion is that code written after the year 2000 should use int, long int, long long int, etc. only when required to do so for compatibility with older libraries.

The primes will never be negative, so it makes more sense to use unsigned types than signed.

## Don't reinvent the wheel

initial_marking is unnecessary. With C-99 you could just initialise list as {} to set them all to 0 and reverse the meaning of the mark; with legacy versions of C you could use memset to set them all to 1.

            if(crossed_out_num <= target/2)//If this condition isn't there, we might me marking numbers not in list.
{
list[crossed_out_num] = 0;
}


list was initialised as long list[target/2]. C uses 0-indexing. So that condition should be < target/2, not <= target/2: as it stands, it can access memory beyond the end of the array.

In this particular case you probably get away with it, because the following chunk of memory probably belongs to primes, but that's not something it's wise to rely on.

There is a similar bug in make_prime_list.

## Use descriptive names

The name list tells me nothing that I wouldn't already know from the type (i.e. that it's an array). What do the elements of the array mean? In this case it seems to be is_prime (or is_half_prime if you want to be pedantic).

marking isn't as bad, but is still fairly undescriptive. sieve would be an improvement; sundaram_sieve would be better still.

print_primes, on the other hand, goes too far. There's nothing in the implementation which limits it to printing primes. It's really print_integers.

## Do you need to statically allocate so much memory?

The number of primes will indeed never exceed target/2, but you could go further and observe that they will never exceed 2 + target/3 (by considerations of values modulo 6), etc. However, the sieve tells you how many primes there are before you try to extract them. Consider this simple modification:

size_t sundaram_sieve(uint8_t is_prime[])
{
size_t primeCount = 1; // Allow for the special case 2
uint32_t i, delta, composite;

//All numbers of the form i + j + 2ij need to be marked out
for (i = 1; i < target/2; i++)
{
primeCount += is_prime[i];

delta = 2*i + 1;
for (composite = 2*i*(i+1); composite < target/2; composite += delta)
{
is_prime[composite] = 0;
}
}

return primeCount;
}


Now primes can be calloced to exactly the size required.

• Thanks. I didn't know about the new data types and that you can initialize the entire array at once. I don't recognize the data types you use at the end. What are they ? Commented Feb 13, 2017 at 23:16
• Do you recommend using the byte data type for list ? Commented Feb 14, 2017 at 0:12
• @user230452, I really recommend writing a bitset implementation as a good exercise, but in the example refactor I've used uint8_t, which is a byte datatype on most modern architectures. The types I've used are (mainly) those from <stdint.h>; size_t is the exception, and that's because it's intended as a parameter to calloc. Commented Feb 14, 2017 at 9:11
• Are you a programmer ? Can you explain what a bit set is ? Is it a binary array ? Commented Feb 14, 2017 at 14:09
• I am. Yes, a bitset is basically a binary array. The typical implementation would be to have an array of some integer type (e.g. uint8_t[] or uint32_t[]) and to set/clear/test individual bits. Commented Feb 14, 2017 at 14:23