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

Question

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([1]*(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.

Answer 1

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.

Answer 2

def sieve(n):
    """bitarray sieve"""
    b=bitarray.bitarray([1]*(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.

Answer 3

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.