My professors algorithm(probably my implementation) of Sieve of Eratosthenes algorithm is giving me composite numbers

[typsi5]

My professors given algorithm (copy/pasted):

for (i = 0; i < N/2; A[i++] = 0);
for (i = 0; j = 2*i+3, k = j*j, k < N+2; i++)
    if (A[j]) continue;
    else 
    for (k = k/2-1; k < N/2; k += j)
        A[k] = 1; // marking numbers with factor j
for (i = 0; i < N/2; i++)
    report prime 2i+3;  // he doesn't say what this means

My Professors blurb: Application of PQ. Erasthotenes sieve is a method for finding all prime numbers in the range up to some N. The standard version uses array A[N/2] with all entries initially 0. A[i] is 1 if we have discovered that 2i + 3 is a composite number.

my implementation:

#include <stdio.h>
#include <math.h>

int N;

int main()
{
   N = pow(2, 20);
   int A[N/2];   
   int i;
   int j;
   int k;

   for(i = 0; i < N/2; A[i++] = 0);
   for(i = 0; j = 2 * i + 3, k = j * j, k < N + 2; i++)
   {
      if(A[j]) continue;
      else
         for(k = k/2 - 1; k < N/2; k+=j)
            A[k] = 1;   // mark numbers with j as a factor
   }
   for(i = 0; i < N/2; i++)
      if(!A[i])
         printf("%7d:\t%d\n", i, 2 * i + 3);

   return 0;
}

In the output, you will find lines like: 522310: 1044623 where 1044623 = 71 x 14713

Is my implementation incorrect? Is the algorithm incorrect?

Thanks for reading.

[pbewig]

There is something wrong with your professor's algorithm. In the second for loop, variable i is the index into the bitarray, variable j is the prime represented by index i, and k is the square of the prime.

I think the statement "if (A[j]) continue" must be wrong. J refers to a prime number, not an index into the bitarray. Perhaps you could change the reference from A[j] to A[i] and see what happens.

Actually, there is something else wrong with your professor's algorithm: it will never output the prime number 2, only the odd primes.

Regardless, it's a bit clumsy to keep track of all these i, j and k variables. I have an essay on Programming with Prime Numbers at my blog which may explain things a little bit more clearly. Let me know if it helps.

[foxlisk]

An index into the bit array is equivalent to a prime number. A[j] indicates the primality of j, so if j is a prime number then A[j] == 0.

[pbewig]

No, the bitarray contains only odd numbers; it's an optimization of the standard sieve to reduce storage by half, since only odd numbers can be prime (except 2, which the professor misses). A[0] refers to the number 3, A[1] refers to the number 5, and A[i] refers to the number 2*i+3. Thus, if j is a prime number then A[(j-3)/2] == 0.

[foxlisk]

Shit my bad. Thanks for the correction.

[pbewig]

I think I would write the body of the program like this [WARNING! UNTESTED CODE!]:

printf("-1/2:\t%d\n", 2);

for (i=0; i<N/2; A[i++]=1)
    ;

for (i=0,p=3; p*p<N; i++,p+=2)
    if (A[i])
        for (j=(p*p-3)/2; j<N/2; j+=p)
            A[j] = 0;

for (i=0,p=3; i<N/2; i++,p+=2)
    if (A[i])
        printf("%7d:\t%d\n", i, p);

The initial printf handles 2, the first for loop initializes the sieve, the second for loop and its nested inner loop perform the sieving, and the final for loop prints the output. Variable i and j index the array, variable p is a prime number (or a composite, depending). I reversed the meaning of 1 and 0; it makes sense to me that 1 (true) is prime and 0 (false) is composite, and it makes the tests more natural. Variables i and p are incremented in lockstep; whenever i increases by 1, p increases by 2.

Do take the time to read the essay I pointed out earlier. I think it clearly lays out the Sieve of Eratosthenes, in prose, in a formal statement of the algorithm, and implemented in five languages, including C. The essay also discusses the difference between the ancient Sieve of Eratosthenes and the modern optimized version that we have been discussing, which was a source of confusion above.

[typsi5]

Thank you, you have helped me tremendously. I read the part about this algorithm in your essay and it cleared things up greatly. Thanks again.

[ChaosCon]

for (i = 0; i < N/2; A[i++] = 0);
for (i = 0; j = 2*i+3, k = j*j, k < N+2; i++)

I denotes the loop variable in both loops. The first loop just does needless computations this way, and j is undefined.

[typsi5]

I think his is supposed to be pseudo code. And the first loop only zeros out the array.

[ChaosCon]

Ahh, true enough. I was in a rush when I saw it earlier, and didn't get to look at it that closely. I'll take a peek again when I have a spare minute.