Find next prime given all prior

Question

I'm writing a recursive infinite prime number generator, and I'm almost sure I can optimize it better.

Right now, aside from a lookup table of the first dozen primes, each call to the recursive function receives a list of all previous primes.

Since it's a lazy generator, right now I'm just filtering out any number that is modulo 0 for any of the previous primes, and taking the first unfiltered result. (The check I'm using short-circuits, so the first time a previous prime divides the current number evenly it aborts with that information.)

Right now, my performance degrades around when searching for the 400th prime (37,813). I'm looking for ways to use the unique fact that I have a list of all prior primes, and am only searching for the next, to improve my filtering algorithm. (Most information I can find offers non-lazy sieves to find primes under a limit, or ways to find the p_nth prime given p_n-1, not optimizations to find p_n given 2...p_n-1 primes.)

For example, I know that the p_nth prime must reside in the range (p_n-1 + 1)...(p_n-1+p_n-2). Right now I start my filtering of integers at p_n-1 + 2 (since p_n-1 + 1 can only be prime for p_n-1 = 2, which is precomputed). But since this is a lazy generator, knowing the terminal bounds of the range (p_n-1+p_n-2) doesn't help me filter anything.

What can I do to filter more effectively given all previous primes?

Code Sample

  @doc """
  Creates an infinite stream of prime numbers.

    iex> Enum.take(primes, 5)
    [2, 3, 5, 7, 11]

    iex> Enum.take_while(primes, fn(n) -> n < 25 end)
    [2, 3, 5, 7, 11, 13, 17, 19, 23]

  """
  @spec primes :: Stream.t
  def primes do
    Stream.unfold( [], fn primes ->
      next = next_prime(primes)
      { next, [next | primes] }
    end )
  end

  defp next_prime([]),      do: 2
  defp next_prime([2 | _]), do: 3
  defp next_prime([3 | _]), do: 5
  defp next_prime([5 | _]), do: 7
  # ... etc

  defp next_prime(primes) do
    start = Enum.first(primes) + 2
    Enum.first(
      Stream.drop_while(
        Integer.stream(from: start, step: 2),
        fn number ->
          Enum.any?(primes, fn prime ->
            rem(number, prime) == 0
          end )
        end
      )
    )
  end

The primes function starts with an empty array, gets the next prime for it (2 initially), and then 1) emits it from the Stream and 2) Adds it to the top the primes stack used in the next call. (I'm sure this stack is the source of some slowdown.)

The next_primes function takes in that stack. Starting from the last known prime+2, it creates an infinite stream of integers, and drops each integer that divides evenly by any known prime for the list, and then returns the first occurrence.

This is, I suppose, something similar to a lazy incremental Eratosthenes's sieve.

You can see some basic attempts at optimization: I start checking at p_n-1+2, and I step over even numbers.

I tried a more verbatim Eratosthenes's sieve by just passing the Integer.stream through each calculation, and after finding a prime, wrapping the Integer.stream in a new Stream.drop_while that filtered just multiples of that prime out. But since Streams are implemented as anonymous functions, that mutilated the call stack.

It's worth noting that I'm not assuming you need all prior primes to generate the next one. I just happen to have them around, thanks to my implementation.

Answer 1

For any number k you only need to try division with primes up to and including √k. This is because any prime larger than √k would need to be multiplied with a prime smaller than √k.

Proof:

√k * √k = k so (a+√k) * √k > k (for all 0<a<(k-√k)). From this follows that (a+√k) divides k iff there is another divisor smaller than √k.

This is commonly used to speed up finding primes tremendously.

Answer 2

You don't need all prior primes, just those below the square root of your current production point are enough, when generating composites from primes by the sieve of Eratosthenes algorithm.

This greatly reduces the memory requirements. The primes are then simply those odd numbers which are not among the composites.

Each prime p produces a chain of its multiples, starting from its square, enumerated with the step of 2p (because we work only with odd numbers). These multiples, each with its step value, are stored in a dictionary, thus forming a priority queue. Only the primes up to the square root of the current candidate are present in this priority queue (the same memory requirement as that of a segmented sieve of E.).

Symbolically, the sieve of Eratosthenes is

     P = {3,5,7,9, ...} \ &bigcup; {{p², p²+2p, p²+4p, p²+6p, ...} | p in P}

Each odd prime generates a stream of its multiples by repeated addition; all these streams merged together give us all the odd composites; and primes are all the odd numbers without the composites (and the one even prime number, 2).

In Python (can be read as an executable pseudocode, hopefully),

def postponed_sieve():      # postponed sieve, by Will Ness,   
    yield 2; yield 3;       #   https://stackoverflow.com/a/10733621/849891
    yield 5; yield 7;       # original code David Eppstein / Alex Martelli 
    D = {}                  #   2002, http://code.activestate.com/recipes/117119
    ps = (p for p in postponed_sieve())  # a separate Primes Supply:
    p = ps.next() and ps.next()          #   (3) a Prime to add to dict
    q = p*p                              #   (9) when its sQuare is 
    c = 9                                #   the next Candidate
    while True:
        if c not in D:                # not a multiple of any prime seen so far:
            if c < q: yield c         #   a prime, or
            else:   # (c==q):         #   the next prime's square:
                add(D,c + 2*p,2*p)    #     (9+6,6 : 15,21,27,33,...)
                p=ps.next()           #     (5)
                q=p*p                 #     (25)
        else:                         # 'c' is a composite:
            s = D.pop(c)              #   step of increment
            add(D,c + s,s)            #   next multiple, same step
        c += 2                        # next odd candidate

def add(D,x,s):                          # make no multiple keys in Dict
    while x in D: x += s                 # increment by the given step
    D[x] = s  

Once a prime is produced, it can be forgotten. A separate prime supply is taken from a separate invocation of the same generator, recursively, to maintain the dictionary. And the prime supply for that one is taken from another, recursively as well. Each needs to be supplied only up to the square root of its production point, so very few generators are needed overall (on the order of log log N generators), and their sizes are asymptotically insignificant (sqrt(N), sqrt( sqrt(N) ), etc).