Sieve of Sundaram

Question

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.

Please provide feedback.

#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]);
    }
}

Answer 1

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.

Answer 2

Check your types

    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.

Check your edge conditions

            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.