Help with Sieve of Atkin algorithm

[munkyeetr]

I am trying to code the Sieve of Atkin algorithm for finding prime numbers. I haven't been able to find any existing VB code, so I have based my algorithm on this C++ implementation.

My function is not correctly marking all the primes. For example, when I try to find the prime numbers to 100 the function only returns:

2 3 5 7 13 19 29 53 67 85

To my eye the code looks correct (in relation to the C++ it is based on), and I have tried to cross check it against the pseudocode in the Wikipedia Article (above), and again to me it looks correct. My math isn't very strong, so I am not sure if I am misinterpreting some of the symbols.

Some new sets of eyes and any advice would be appreciated.

Private Function AtkinSieve(ByVal limit As Integer) As Integer()
    Dim root As Integer = Math.Ceiling(Math.Sqrt(limit))
    Dim sieve(limit) As Boolean
    Dim primes(1) As Integer
    Dim insert As Integer = 2
    primes(0) = 2
    primes(1) = 3

    For z = 0 To limit
        sieve(z) = False
    Next

    ' Main part of Sieve of Atkin
    Dim x As Integer = 1
    Dim y As Integer = 1
    While x <= root
        While y <= root
            Dim n As Integer = (4 * (x * x)) + (y * y)
            If (n <= limit And (n Mod 12 = 1 Or n Mod 12 = 5)) Then sieve(n) ^= True
            n = (3 * (x * x)) + (y * y)
            If (n <= limit And (n Mod 12 = 7)) Then sieve(n) ^= True
            n = (3 * (x * x)) - (y * y)
            If (x > y And n <= limit And n Mod 12 = 11) Then sieve(n) ^= True
            y += 1
        End While
        x += 1
    End While

    ' Mark all multiples of squares as non-prime
    Dim r As Integer = 5
    Dim i As Integer = 0
    While r <= root
        If sieve(r) Then
            i = (r ^ 2)
            While i < limit
                sieve(i) = False
                i += (r ^ 2)
            End While
        End If
        r += 1
    End While

    ' Add into the prime array
    Dim a As Integer = 5
    While a < limit
        If sieve(a) Then
            Redim Preserve primes(insert)
            primes(insert) = a
            insert += 1
        End If
        a += 1
    End While

    Return primes
End Function

 

[munkyeetr]

So, I've made a little progress debugging this and found that in main part of the algorithm (the x,y loops) I had declared y outside the x loop and it wasn't being reinstantiated to 1 with each iteration of x. It should read (and I am making this change in the original post as well):

' Main part of Sieve of Atkin
    Dim x As Integer = 1
    While x <= root
        [B]Dim y As Integer = 1[/B]
        While y <= root
            Dim n As Integer = (4 * (x * x)) + (y * y)
            If (n <= limit And (n Mod 12 = 1 Or n Mod 12 = 5)) Then sieve(n) ^= True
            n = (3 * (x * x)) + (y * y)
            If (n <= limit And (n Mod 12 = 7)) Then sieve(n) ^= True
            n = (3 * (x * x)) - (y * y)
            If (x > y And n <= limit And n Mod 12 = 11) Then sieve(n) ^= True
            y += 1
        End While
        x += 1
    End While

I thought I had fixed it and it seemed to return all the primes when I re-ran it, but it actually returns too many numbers as prime.

Output calculating primes to 100:

2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 65 67 71 73 79 83 85 89 91 97

65, 85 and 91 shouldn't be marked as prime.

Still looking for fresh eyes on this. Thanks.

 

[munkyeetr]

Anyway, after a little bit more sleuthing I found that the issue was in the conversion from C++ to VB, specifically this piece of code:

While y <= root
            Dim n As Integer = (4 * (x * x)) + (y * y)
            If (n <= limit And (n Mod 12 = 1 Or n Mod 12 = 5)) Then sieve(n) ^= True
            n = (3 * (x * x)) + (y * y)
            If (n <= limit And (n Mod 12 = 7)) Then sieve(n) ^= True
            n = (3 * (x * x)) - (y * y)
            If (x > y And n <= limit And n Mod 12 = 11) Then sieve(n) ^= True
            y += 1
        End While

I had copied the '^= True' pieces from the original C++ and it didn't flag as a syntax error. I played around quite a bit (no pun intended) with those statements. I tried leaving off the ^ and tried to substitute Xor, etc in many different ways, but it made no difference in my code output so I had pretty much excluded those parts as the problem.

Until I fixed the problem with the x,y loop (see above post).

Now by flipping the bit with

sieve(n) = not sieve(n)

the algorithm works.

Marking as SOLVED.

 

[mscir]

I would declare all variables, x, y, and n outside of the loops.

I found this related discussion

From MSDN:

Even if the scope of a variable is limited to a block, its lifetime is still that of the entire procedure. If you enter the block more than once during the procedure, each block variable retains its previous value. To avoid unexpected results in such a case, it is wise to initialize block variables at the beginning of the block.
This article might help you (variable scope and lifetime)
Scope in Visual Basic

.net - variable declaration inside a loop - Stack Overflow

 

[munkyeetr]

Thanks mscir, I'll take a look at that.