5
\$\begingroup\$

I wrote a simple prime sieve in Python. But it ended up looking... ugly. I included a few of the most basic optimizations. When posting this question SE suggested a bunch of questions on the same topic, which lead me to some improvements. But few were in generator form.

It's the optimizations' fault.
Unfortunately I can't continue an outer loop from an inner loop, like I can in JavaScript. I very much disagree with Guido here, as it makes the inner loop clunky.

from itertools import count

def sieve():
    primes = [2]
    for candidate in count(start=3, step=2):
        cont = False
        n = 0
        while primes[n]**2 <= candidate:  # You only need to check up to the square root of a number.
            if candidate % primes[n] == 0:
                cont = True  # outer
                break
            n = n + 1
        if cont:
            cont = False
            continue
        yield primes[-1]
        primes.append(candidate)

Can we make this more concise without changing the basic logic? For example this is extremely concise, but it doesn't have some of the optimizations my code does.

For fun, I wrote the same logic in Javascript. This looks cleaner due to being able to continue an outer loop. But the lack of negative indices is a step backwards.

function* sieve() {
    let primes = [2]
    counter: for (let candidate = 3;; candidate+=2) {
        for (let n = 1; primes[n]**2 <= candidate; n++)
            if (candidate % primes[n] == 0)
                continue counter
        yield primes[primes.length - 1]
        primes.push(candidate)
    }
}
\$\endgroup\$
  • \$\begingroup\$ Please pick one language per question. Do you want the Python or the JavaScript reviewed? \$\endgroup\$ – Peilonrayz May 12 at 18:30
  • \$\begingroup\$ Hm, if I have to choose, I guess I'll go with python. \$\endgroup\$ – martixy May 12 at 18:31
  • 3
    \$\begingroup\$ Note that what you are implementing here is not the Sieve of Eratosthenes, but rather the oft-mistaken and much slower Trial Division algorithm. The Sieve has no divisions nor modulo operations in it and accomplishes everything through additions or, at worst, some multiplications. \$\endgroup\$ – RBarryYoung May 13 at 13:17
  • \$\begingroup\$ You can compare Trial Division (en.wikipedia.org/wiki/Trial_division) to the Sieve of Eratosthenes here: en.wikipedia.org/wiki/Sieve_of_Eratosthenes. \$\endgroup\$ – RBarryYoung May 13 at 13:28
6
\$\begingroup\$
  1. We can move the while loop into another function.

    def is_prime(candidate, primes):
        n = 0
        while primes[n]**2 <= candidate:
            if candidate % primes[n] == 0:
                return False
            n += 1
        return True
    
    def sieve():
        primes = [2]
        for candidate in count(start=3, step=2):
            if not is_prime(candidate, primes):
                continue
            yield primes[-1]
            primes.append(candidate)
    
  2. We can use itertools.takewhile to expose the while loop as an iterable.

    def is_prime(candidate, primes):
        for prime in itertools.takewhile(lambda p: p**2 <= candidate, primes):
            if candidate % prime == 0:
                return False
        return True
    
  3. We can use any to make is_prime easier to read.

    If for any of the values, candidate % prime == 0, are true the result is true. If none of them are then it is false. Since we want it to be the other way around we can just use not.

    def is_prime(candidate, primes):
        return not any(
            candidate % prime == 0
            for prime in itertools.takewhile(lambda p: p**2 <= candidate, primes)
        )
    
  4. We can move is_prime back into the first function.

    def sieve():
        primes = [2]
        for candidate in count(start=3, step=2):
            if not not any(
                candidate % prime == 0
                for prime in itertools.takewhile(lambda p: p**2 <= candidate, primes)
            ):
                continue
            yield primes[-1]
            primes.append(candidate)
    
  5. We can swap the if to cancel the continue and a not.

  6. We can swap the any and == 0, with not all.
  7. We can cancel the double not.
def sieve():
    primes = [2]
    for candidate in count(start=3, step=2):
        if all(
            candidate % prime
            for prime in itertools.takewhile(lambda p: p**2 <= candidate, primes)
        ):
            yield primes[-1]
            primes.append(candidate)

At the expense of readability you can get the following trade offs.

  1. There is a chance that using int(candidate ** 0.5).__ge__ is faster than the lambda.
| improve this answer | |
\$\endgroup\$
  • \$\begingroup\$ While this is a good explanation of the approach, I am sad to report that it is also wrong - the logic is broken. It thinks 9 is a prime. More generally, it considers all odd squares prime. \$\endgroup\$ – martixy May 13 at 4:35
  • \$\begingroup\$ Though the fix is simple - the lambda condition must be <=, not <. I do like how mathematically / logically this is more expressive than the other answer. It is also slower(15.1s vs 12.9s over n=200K). TBH I can't decide which answer is better... \$\endgroup\$ – martixy May 13 at 4:40
  • \$\begingroup\$ @martixy Ah yes I dropped the equals there, but somehow kept it in the 'extra' part - __ge__. Thank you. You didn't highlight that you wanted this super optimized, however if you do then another one of my answers likely performs better than either of these answers. \$\endgroup\$ – Peilonrayz May 13 at 4:48
  • \$\begingroup\$ Not specifically "super optimized", I was just listing pros and cons. \$\endgroup\$ – martixy May 13 at 4:55
  • 1
    \$\begingroup\$ I like everything in this answer except the implicit conversion of a remainder to a truthy value: to me, candidate % prime != 0 is vastly more readable, despite being a minimal change. \$\endgroup\$ – Konrad Rudolph May 13 at 9:09
5
\$\begingroup\$

This outputs the same as your function:

import array
from itertools import count


def sieve():
    primes = array.array("L", [2])  # L: unsigned long int
    for candidate in count(start=3, step=2):
        n = 0
        while primes[n] ** 2 <= candidate:
            if candidate % primes[n] == 0:
                break
            n += 1
        else:  # nobreak
            yield primes[-1]
            primes.append(candidate)


n = 0
primes = sieve()

while True:
    prime = next(primes)
    print(prime)
    n += 1

I did not touch the algorithm itself because I am unfamiliar, but Python (tested on 3.8.2, but also runs on 2.7.18, as you seem to need given your tag ) has the while/else construct that can help you with your control flow here.

A break break, break <label> or similar thinkable constructs have been proposed, as you linked as well, but rejected. A possible remedy is to extract sub-routines into functions and to use their return statements for control flow/multi-level breaking. In your case, that is not needed.

The while/else construct is unfamiliar to many. In fact, it is so alien that Guido would not implement it again nowadays. For now, it is best to think of the else as nobreak: the while loop finished normally (its condition evaluated to False) and exited. Since no break occurred, the else block is executed.

In the above case, if break is hit in the if block, the else is skipped: no yield occurs, and since after the else block, there is no code left, a continue for the outer for loop is implied and not explicitly necessary, since there is nothing else to do anyway.


A frequent example is in the form of for/else (which works like the while/else) when looking for a hit, like in your case:

for file in files:
    if file == file_looked_for:
        file.do_something()
        break
else:  # nobreak
    # code to handle file not found

So while, like you, I have been tripped up by the lack of advanced break functionalities, I agree with Guido. Apart from the unfortunate naming of else in the while construct (more discussion here), it can do the job just fine. Trying to break through multiple levels is an occasion to rethink the implementation.


Other observations:

  • primes[n] has to call __getitem__, which happens twice. This is constant-time, but the following will probably provide a speed-up (but requires the "walrus" operator from Python 3.8):

      while (
          (prime := primes[n]) ** 2 <= candidate
      ):  # You only need to check up to the square root of a number.
          if candidate % prime == 0:
    

    primes no longer has to be looked up, the simple, local prime suffices for the second call.

  • n = n + 1 can be n += 1. However, this is possibly slower.
  • Finally, the built-in array.array can be much faster than a list. This is because lists can hold arbitrary objects, whereas an array has to be initiliazed for a specific one (here: L for long int, which will last a while). As a positive consequence, the array can then be optimized accordingly.

    Note how array.array, in this case, is a 1:1 drop-in for the previous list. As such, the code did not have to change. As such, it is also easy for you to revert the change if unwanted.

| improve this answer | |
\$\endgroup\$
  • 1
    \$\begingroup\$ You can change the calling code in your first code block to just for prime in sieve(): print(prime). \$\endgroup\$ – Peilonrayz May 12 at 19:46
  • \$\begingroup\$ Yes you are right. There was a guard clause to break at n == 20 that felt more natural in a while True, but I removed it before posting. \$\endgroup\$ – Alex Povel May 12 at 19:50
  • 2
    \$\begingroup\$ In the future you may wish to do for i, prime in enumerate(sieve()): if i == 20: break it's more 'Pythonic' ;) \$\endgroup\$ – Peilonrayz May 12 at 19:51
  • 1
    \$\begingroup\$ alien > You can say that again. I initially thought you'd messed up your indentation before reading the explanation. While the construct looks odd (and else is a questionable choice of words), it does make it pretty neat - and faster. For the other observations: Walrus operator - this is just a straight up a holy crap, wtf python? moment - why is the same thing that works for a million other languages not enough for python? Why does there have to be a special operator? Also provides at best a marginal improvement (12.9s > 12.6s over n=200K). Continued... \$\endgroup\$ – martixy May 13 at 3:54
  • 1
    \$\begingroup\$ And I am sad to report that array.array actually results in a slowdown. For n += 1 - I know python doesn't have ++ and I was under the mistaken impression shorthand assignment is also not in there. Regardless, if n += 1 is slower, it is not observable over n=200K. Thank you for the review. It's a good explanation of a language construct I've never encountered before. The other observations also fit in the spirit of the question, even if they end up being hit or miss. Python is just like that special child(walrus - seriously? I can't even...) \$\endgroup\$ – martixy May 13 at 4:03
2
\$\begingroup\$

Three main things vex me about your code:

  • the just generated prime isn't returned,
  • the boolean test a the end of the while loop
  • repeated prime "square and test" for candidacy.

Not returning generated prime

Your code generates 3, and yields 2, then it generates 5 and yields 3, then it generates 7 and yields 5, then generates 11 and yields 7, and so on.

This happens because you half treat 2 as a special case. You initialize the primes array with it. But to return it, you use yield primes[-1] just like every other prime.

If you treated it completely as a special case, and yield it right off the hop, you could yield candidate at the end of the loop, thus returning the prime you just computed.

def sieve():
    primes = [2]
    yield 2

    for candidate in count(start=3, step=2):
        ...
        yield candidate
        primes.append(candidate)

Unnecessary boolean test at end of while

A while loop is often used for searching. If the value is found, the while loop is escaped via a break statement. If the while loop completes without ever breaking, the condition searched for was never found, and something different needs to happen. In Python, this is the while ... else statement:

def sieve():
    primes = [2]
    yield 2
    for candidate in count(start=3, step=2):
        n = 0
        while primes[n]**2 <= candidate:  # Only check up to the square root of number.
            if candidate % primes[n] == 0:
                break
            n = n + 1
        else:
            yield candidate
            primes.append(candidate)

Repeated prime "square and test" for candidacy.

How often is the primes[n]**2 <= candidate done?

If candidate is just over 10,000, and is prime, then we will be squaring all primes less than 100, and testing that they are less than candidate. Then, we do the same thing for candidate + 2, and the results will be the same. No prime number less than 100, squared, will ever be greater than candidate once candidate exceeds 10,000 ... so this is all busy work, repeating the same test over and over.

What you need is to partition your primes list into two parts: primes less or equal to the square-root of candidate, and primes greater the square-root of candidate.

You can do this in several ways. The smallest change would be to keep track of a count of "small" primes. As candidate gets larger by 2, you would only need add at most one more prime into the "small" primes bucket:

def sieve():
    primes = [2]
    yield 2

    small_primes = 0
    for candidate in count(start=3, step=2):

        if primes[small_primes] ** 2 <= candidate:
            small_primes += 1

        for n in range(small_primes):
             if candidate % primes[n] == 0:
                break
        else:
            yield candidate
            primes.append(candidate)

Now how often is primes[small_primes] ** 2 <= candidate being done? Once per candidate! This has got to be an improvement. Also, all n = 0 and n = n + 1 code has been absorbed into for n in range(small_primes), and having Python do this work is faster than coding it ourselves.

Other improvements

Odd numbers

Why are we test-dividing all of our candidates by primes[0] == 2? By design, they are all odd, and can never be evenly divided by 2.

        for n in range(1, small_primes):    # Skip divide-by-2 tests

All

As mentioned by Peilonrayz, Python has an any() function, though I think all() is more appropriate here.

def sieve():
    primes = [2]
    yield 2

    small_primes = 0
    for candidate in count(start=3, step=2):

        if primes[small_primes] ** 2 <= candidate:
            small_primes += 1

        if all(candidate % primes[n] != 0 for n in range(1, small_primes)):
            yield candidate
            primes.append(candidate)

Maintain separate lists

Instead of small_primes being a count of the number of primes less than the square-root of the candidate, what if it actually was a list of the small prime numbers? And instead of adding prime candidates to that list, we add to a large_primes list? Then we could move primes from the large_primes to the small_primes as the square-root of the candidate increases.

Optimizations:

  • deque for large_primes
  • Omit 2 from the small_primes list,
  • Cache the large_prime[0] ** 2 value, to avoid repeatedly squaring the same quantity.

Resulting code:

from itertools import count
from collections import deque

def sieve():
    yield 2
    yield 3

    small_primes = []
    large_primes = deque((3,))
    next_prime_squared = large_primes[0] ** 2

    for candidate in count(start=5, step=2):

        if candidate >= next_prime_squared:
            small_primes.append(large_primes.popleft())
            next_prime_squared = large_primes[0] ** 2

        if all(candidate % prime != 0 for prime in small_primes):
            yield candidate
            large_primes.append(candidate)

Time Comparisons

Time (in seconds) for generating 100 to 100,000 primes: enter image description here

Timing code:

import array
from timeit import timeit
from itertools import count, takewhile
from collections import deque
import matplotlib.pyplot as plt

def martixy():
    primes = [2]
    for candidate in count(start=3, step=2):
        cont = False
        n = 0
        while primes[n]**2 <= candidate:  # You only need to check up to the square root of a number.
            if candidate % primes[n] == 0:
                cont = True  # outer
                break
            n = n + 1
        if cont:
            cont = False
            continue
        yield primes[-1]
        primes.append(candidate)

def alex_povel():
    primes = array.array("L", [2])  # L: unsigned long int
    for candidate in count(start=3, step=2):
        n = 0
        while primes[n] ** 2 <= candidate:
            if candidate % primes[n] == 0:
                break
            n += 1
        else:  # nobreak
            yield primes[-1]
            primes.append(candidate)

def peilonrayz():
    primes = [2]
    for candidate in count(start=3, step=2):
        if all(
            candidate % prime
            for prime in takewhile(lambda p: p**2 <= candidate, primes)
        ):
            yield primes[-1]
            primes.append(candidate)

def ajneufeld():
    yield 2
    yield 3

    small_primes = []
    large_primes = deque((3,))
    next_prime_squared = large_primes[0] ** 2

    for candidate in count(start=5, step=2):

        if candidate >= next_prime_squared:
            small_primes.append(large_primes.popleft())
            next_prime_squared = large_primes[0] ** 2

        if all(candidate % prime != 0 for prime in small_primes):
            yield candidate
            large_primes.append(candidate)

def test(candidate, limit):
    sieve = candidate()
    for _ in range(limit):
        next(sieve)

if __name__ == '__main__':
    candidates = (martixy, alex_povel, peilonrayz, ajneufeld)
    limits = [int(10 ** (power * 0.25)) for power in range(8, 21)]

    fig, ax = plt.subplots()

    for candidate in candidates:
        print("Testing", candidate.__name__)
        times = [ timeit(lambda: test(candidate, limit), number=1) for limit in limits ]
        ax.plot(limits, times, '-+', label=candidate.__name__)

    ax.legend()
    plt.show()
| improve this answer | |
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.