# Finds all the prime numbers up to n using a strange pattern

Here's an odd pattern I found for sieving out composite numbers to find primes. I wrote it myself, originally thinking I made an indexing error, but it turned out to work after removing squares.

import bitarray
def sieve(n):
"""bitarray sieve"""
b=bitarray.bitarray(*(n+1))
x=int(n**.5)+1
for k in range(2,x):
if b[k]:
f=k+1 #this line is indexing error
e=n//k +1
for g in range(f,e):
if b[g]:
b[g*k]=0
return sorted(list(to_numbs(b)-squares_set(b)))

def squares_set(ba):
"""this returns a set of all primes**(2n).
note 0**2==0, 1**2==1"""
return {i**2 for i in range(int(len(ba)**.5)+1) if ba[i]}

def to_numbs(b_ar):
"""converts the bit array to a set of numbers"""
return {i for i in range(len(b_ar)) if b_ar[i]}
if __name__=="__main__":
print(sieve(100))


I'm not entirely sure what causes this pattern to work so I figured it would be fun to post here.

The Sieve of Eratosthenes works by going through and marking all numbers that can be divided into, essentially.

Take your input "n", for which you want to find all the prime numbers up to it. The sieve starts at 2, and multiplies 2 by increasing amounts until the multiplication is about to exceed the value "n" - in your program, this is accomplished through the integer division and then iteration up to the calculated result. So, as it multiplies the values together, it marks the result as non-prime. Meaning that 2, 4, 6, 8, 10... etc are all set in the bitarray as 0 (meaning non-prime in this context). It moves on to the next value, so it starts at 3. The process repeats as it marks all multiples as non-prime. Time to move on to the next number - 4. As 4 has been marked as non-prime, it skips over 4 and moves on to 5. This continues until it finally reaches the number which is the value of its square root, so when you are using the sieve to get primes below 100, that means the algorithm checks the multiples of values from 2 to sqrt(100), so 2 to 10. It doesn't need to go past this value as all prime values above that share multiples which values that have already been checked. So, for example, if we moved past 10 to check multiples of 11, we'll get 22 (which is divisible by 2, so has already been marked non-prime) and 33 (divisible by 3). The values which remain with a value of 1 are prime.

tl;dr The algorithm goes through all prime values below or equal to sqrt(n) and marks their multiples up to n as non-prime. The remaining numbers must be prime.

EDIT: The code above adds 1 to the div operation, which I didn't mention in my explanation.

EDIT 2: This specific modification of the sieve is essentially the same as the original sieve, however rather than removing square values during the normal removal process it removes them separately at the end.

• If you put a print(gk) after b[gk]=0, you will see every number that it marks composite in the order it marks them. You will also see that this order is not what you've described. - edit: new user- ignore the italics – mikeLundquist Jul 24 '16 at 12:43
• Throwing some prints into the code, it shows on my end near exactly how I described - review link to see each calculation that causes a value to be set to 0. It starts at 2 and multiplies in increasing amounts until it reaches n//2, which is 50 (actually, the code adds 1 after the div operation, however I didn't mention that in my original post). In the output, you'll see that each result isn't always increasing by 2 - this is because of the check b[g], which checks for primes that have already been found. – VortixDev Jul 24 '16 at 13:23
• If b[g*k] would check if composites were already found. There are many examples of composites that are deleted after one of their factors is looped (for example 28). – mikeLundquist Jul 24 '16 at 15:08
• additionally, a sieve of eratosthenes wouldn't loop over 4. – mikeLundquist Jul 24 '16 at 15:11
• Not sure about other examples, but 28 was only marked as non-prime when it was calculated using g*k. If you check the file, you'll see that once 28 was calculated it was never seen again in the file, because it b[g] didn't evaluate to true. As for the sieve looping over 4, that is because of f=k+1 starting the multiplication at 2*3 = 6, meaning 2*2 = 4 isn't calculated – VortixDev Jul 24 '16 at 15:35
def sieve(n):
"""bitarray sieve"""
b=bitarray.bitarray(*(n+1))


That's not a helpful comment. It's obviously a bitarray, but what does it mean that bit i is clear or set?

            f=k+1 #this line is indexing error
e=n//k +1
for g in range(f,e):
if b[g]:
b[g*k]=0


The test is unnecessary unless clearing a bit is significantly more expensive than testing. And once the test is removed, the bookkeeping can be done by range, giving

            for g in range(k*(k+1), n+1, k):
b[g]=0


(at which point it should be obvious that there's an out-by-one error, so you could fix it).

    return sorted(list(to_numbs(b)-squares_set(b)))


Why should it be necessary to sort? If it is necessary then I would definitely want to change the way the sequence is generated to make it unnecessary. And at that point it might not be necessary to coerce to list:

    squares = set({i*i for i in range(x)})
return {i for i in range(n+1) if b[i] and not (i in squares)}


Although obviously it's preferable to fix the out-by-one error and ditch squares entirely.

• good point on the sorted. – mikeLundquist Nov 11 '16 at 0:58
• testing is less expensive than multiplication, not clearing (the multiplication also explains why this isn't an off by one error. I'm multiplying, not adding.). Honestly I really need to make a YouTube video stepping through the flow of this, it's not the same as the sieve of Eratosthenes. ps this was an edit but it timed out – mikeLundquist Nov 11 '16 at 1:04

The 'indexing error' causes the crossing-off sequence for a prime p to start at (p + 1) * p instead of at p * p, which is why the code needs a wasteful separate pass for removing the squares of primes. VortixDev has already explained this several times. Moreover, since squares of primes aren't eliminated when they should, the code does superfluous crossing-off passes for the squares of primes.

P.S.: didn't see that the question is several months old - I picked it off the top of the active listing.

• @Peter: You are right, thank you for pointing out this silly mistake (too much language-hopping, I guess). I've removed the offending part of the text. – DarthGizka Nov 10 '16 at 15:42
• no, seriously, if you change it to p instead of p+1 it doesn't work. – mikeLundquist Nov 11 '16 at 1:05
• @user: That's because the superfluous test if b[g] messes things up. If it were if b[g*k] then it would be harmless, a standard pattern employed when packed bit writing is much more expensive than packed bit testing. But it isn't. Instead it interferes with the crossing-off sequence of the sieve. Eliminate the test or change it to if b[g*k] and your sieve should start working normally. You really should heed Peter Taylor's remark on that score though (step additively from k*k to n with stride k instead of multiplying, i.e. range(k * k, n + 1, k)) – DarthGizka Nov 11 '16 at 6:22
• The test keeps you from eliminating multiples g * k of the current prime k that have been crossed off because g is a multiple of a smaller prime (crossed off earlier). Many of the skipped crossings-off are redundant, but not all. I guess you'll find if you simplify your code by removing the test and changing the crossing-off to additive stepping then it'll actually get faster - the speed supremacy of the Sieve of Eratosthenes derives from its simplicity, so don't sabotage that by replacing b[j] = 0 with if b[g*k]: b[g*k] = 0. Especially in an interpreted language like Python. – DarthGizka Nov 11 '16 at 6:34
• The point of this isn't really speed, the point is more to show this weird pattern in the non-prime numbers. – mikeLundquist Nov 16 '16 at 17:42