My favourite prime quantity generator

The sieve of Eratosthenes

To grasp what this code does, we should always first begin with the essential Sieve
of Eratosthenes; when you’re aware of it, be happy to skip this part.

The Sieve of Eratosthenes is a well known
algorithm from historic Greek instances for locating all of the primes under a sure
quantity fairly effectively utilizing a tabular illustration. This animation
from Wikipedia explains it fairly nicely:

Beginning with the primary prime (2) it marks all its multiples till the requested
restrict. It then takes the subsequent unmarked quantity, assumes it is a prime (as a result of it
is just not a a number of of a smaller prime), and marks its multiples, and so forth
till all of the multiples under the restrict are marked. The remaining
unmarked numbers are primes.

This is a well-commented, primary Python implementation:

def gen_primes_upto(n):
    """Generates a sequence of primes < n.

    Makes use of the total sieve of Eratosthenes with O(n) reminiscence.
    """
    if n == 2:
        return

    # Initialize desk; True means "prime", initially assuming all numbers
    # are prime.
    desk = [True] * n
    sqrtn = int(math.ceil(math.sqrt(n)))

    # Beginning with 2, for every True (prime) quantity I within the desk, mark all
    # its multiples as composite (beginning with I*I, since earlier multiples
    # ought to have already been marked as multiples of smaller primes).
    # On the finish of this course of, the remaining True objects within the desk are
    # primes, and the False objects are composites.
    for i in vary(2, sqrtn):
        if desk[i]:
            for j in vary(i * i, n, i):
                desk[j] = False

    # Yield all of the primes within the desk.
    yield 2
    for i in vary(3, n, 2):
        if desk[i]:
            yield i

After we need a listing of all of the primes under some recognized restrict,
gen_primes_upto is nice, and performs pretty nicely. There are two points
with it, although:

1. We now have to know what the restrict is forward of time; this is not all the time attainable
   or handy.
2. Its reminiscence utilization is excessive – O(n); this may be considerably optimized,
   nonetheless; see the bonus part on the finish of the submit for particulars.

The infinite prime generator

Again to the infinite prime generator that is within the focus of this submit. Right here is
its code once more, now with some feedback:

def gen_primes():
    """Generate an infinite sequence of prime numbers."""
    # Maps composites to primes witnessing their compositeness.
    D = {}

    # The operating integer that is checked for primeness
    q = 2

    whereas True:
        if q not in D:
            # q is a brand new prime.
            # Yield it and mark its first a number of that is not
            # already marked in earlier iterations
            D[q * q] = [q]
            yield q
        else:
            # q is composite. D[q] holds among the primes that
            # divide it. Since we have reached q, we not
            # want it within the map, however we'll mark the subsequent
            # multiples of its witnesses to arrange for bigger
            # numbers
            for p in D[q]:
                D.setdefault(p + q, []).append(p)
            del D[q]

        q += 1

The important thing to the algorithm is the map D. It holds all of the primes encountered
to this point, however not as keys! Relatively, they’re saved as values, with the keys being
the subsequent composite quantity they divide. This lets this system keep away from having to
divide every quantity it encounters by all of the primes recognized to this point – it could merely
look within the map. A quantity that is not within the map is a brand new prime, and the best way
the map updates is just not not like the sieve of Eratosthenes – when a composite is
eliminated, we add the subsequent composite a number of of the identical prime(s). That is
assured to cowl all of the composite numbers, whereas prime numbers ought to by no means
be keys in D.

I extremely suggest instrumenting this operate with some printouts and operating
via a pattern invocation – it makes it simple to grasp how the algorithm
makes progress.

In comparison with the total sieve gen_primes_upto, this operate does not require
us to know the restrict forward of time – it should preserve producing prime numbers advert
infinitum (however will run out of reminiscence finally). As for reminiscence utilization, the
D map has all of the primes in it someplace, however every one seems solely as soon as.
So its measurement is , the place is the
Prime-counting function,
the variety of primes smaller or equal to n. This may be
approximated by .

I do not bear in mind the place I first noticed this method talked about, however all of the
breadcrumbs result in this ActiveState Recipe by David Eppstein from
manner again in 2002.

Optimizing the generator

I actually like gen_primes; it is quick, simple to grasp and offers me as
many primes as I would like with out forcing me to know what restrict to make use of, and its
reminiscence utilization is way more cheap than the full-blown sieve of Eratosthenes.
It’s, nonetheless, additionally fairly sluggish, over 5x slower than gen_primes_upto.

The aforementioned ActiveState Recipe thread has a number of optimization concepts;
here is a model that comes with concepts from Alex Martelli, Tim Hochberg and
Wolfgang Beneicke:

def gen_primes_opt():
    yield 2
    D = {}
    for q in itertools.rely(3, step=2):
        p = D.pop(q, None)
        if not p:
            D[q * q] = q
            yield q
        else:
            x = q + p + p  # get odd multiples
            whereas x in D:
                x += p + p
            D[x] = p

The optimizations are:

See Also

Vulkan Foliage rendering utilizing GPU Instancing

1. As an alternative of holding an inventory as the worth of D, simply have a single quantity.
   In circumstances the place we’d like a couple of witness to a composite, discover the subsequent
   a number of of the witness and assign that as a substitute (that is the whereas x in D
   internal loop within the else clause). This can be a bit like utilizing linear probing
   in a hash desk as a substitute of getting an inventory per bucket.
2. Skip even numbers by beginning with 2 after which continuing from 3 in steps
   of two.
3. The loop assigning the subsequent a number of of witnesses might land on even numbers
   (when p and p are each odd). So as a substitute leap to q + p + p
   immediately, which is assured to be odd.

With these in place, the operate is greater than 3x quicker than earlier than, and is
now solely inside 40% or so of gen_primes_upto, whereas remaining quick and
fairly clear.

There are even fancier algorithms that use fascinating mathematical tips to do
much less work. This is an approach by Will Ness and Tim Peters (sure, that Tim Peters) that is
reportedly quicker. It makes use of the wheels concept from this paper by Sorenson. Some extra
particulars on this method can be found here. This algorithm is each quicker and
consumes much less reminiscence; then again, it is not quick and easy.

To be trustworthy, it all the time feels a bit odd to me to painfully optimize Python code,
when switching languages supplies vastly larger advantages. For instance, I threw
collectively the identical algorithms using Go
and its experimental iterator support; it is 3x quicker than the
Python model, with little or no effort (although the brand new Go iterators and
yield features are nonetheless within the proposal stage and are not optimized). I
cannot attempt to rewrite it in C++ or Rust for now, as a result of lack of generator
assist; the yield assertion is what makes this code so good and chic,
and various idioms are a lot much less handy.

Bonus: segmented sieve of Eratosthenes

The Wikipedia article on the sieve of Eratosthenes mentions a segmented
approach, which
can also be described within the Sorenson paper in part 5.

The primary perception is that we solely want the primes as much as to
have the ability to sieve a desk all the best way to N. This ends in a sieve that makes use of
solely reminiscence. This is a commented Python implementation:

def gen_primes_upto_segmented(n):
    """Generates a sequence of primes < n.

    Makes use of the segmented sieve or Eratosthenes algorithm with O(√n) reminiscence.
    """
    # Simplify boundary circumstances by hard-coding some small primes.
    if n < 11:
        for p in [2, 3, 5, 7]:
            if p < n:
                yield p
        return

    # We break the vary [0..n) into segments of size √n
    segsize = int(math.ceil(math.sqrt(n)))

    # Find the primes in the first segment by calling the basic sieve on that
    # segment (its memory usage will be O(√n)). We'll use these primes to
    # sieve all subsequent segments.
    baseprimes = list(gen_primes_upto(segsize))
    for bp in baseprimes:
        yield bp

    for segstart in range(segsize, n, segsize):
        # Create a new table of size √n for each segment; the old table
        # is thrown away, so the total memory use here is √n
        # seg[i] represents the quantity segstart+i
        seg = [True] * segsize

        for bp in baseprimes:
            # The primary a number of of bp on this section may be calculated utilizing
            # modulo.
            first_multiple = (
                segstart if segstart % bp == 0 else segstart + bp - segstart % bp
            )
            # Mark all multiples of bp within the section as composite.
            for q in vary(first_multiple, segstart + segsize, bp):
                seg[q % len(seg)] = False

        # Sieving is finished; yield all composites within the section (iterating solely
        # over the odd ones).
        begin = 1 if segstart % 2 == 0 else 0
        for i in vary(begin, len(seg), 2):
            if seg[i]:
                if segstart + i >= n:
                    break
                yield segstart + i

Code

The total code for this submit – together with checks and benchmarks – is offered
on GitHub.