Problem 7 – The nth prime

by

I decided to try my hand at problem 7, which is to determine the n-th prime, for some chosen n. I’m sure there’s more math out there on this subject than I could ever read and understand, so I’ll take a pretty simple view on things.

It doesn’t get much simpler, in terms of primes algorithms, than the Sieve of Eratosthenes. I’ll be entirely naive about things, and not even toss out evens from my list, and start my list at 0, even though I know the first prime is 2. This keeps the array indexing trivial (let’s focus on one problem at a time)

def sieveto(n):
“”” Determine all the primes up to n

The basic sieve of Eratosthenes, with some memory of the
implementation from numericalrecipes.wordpress.com
“””
isprime = [True for n in xrange(0,n+1)]
isprime[0],isprime[1] = False, False
idx = 2
while(idx <= int(math.sqrt(n))): for fi in xrange(idx*idx,n+1,idx): isprime[fi] = False idx += 1 while(not isprime[idx]): idx += 1 return [i for i in xrange(2,n+1) if isprime[i]] [/sourcecode] The difficulty in this problem is that we don't know the bound. In fact, you could take the problem to be "find the bound that you need to run the sieve to." So, I decided to build some initial string of primes, up to some chosen bound, and then keep extending my list until I had enough primes. If I already have a list of all the primes up to some prime $latex p$, then it's easy to get a list of all the primes up to $latex p^2$. When I'd run the sieve up to $latex p^2$, the primes, whose multiples I'd cross out, are all those up to $latex p$. So I just need to go through the primes I already know, cross out all their multiples, and I'll get up to $latex p^2$. [sourcecode language="python"] def extend(smpr): """ Take a list of primes, find primes up to the square of the largest Assumes the input list is in increasing order, and contains all of the primes up to the largest prime it contains (smpr[-1]). Return a list of the new primes """ origmax = smpr[-1] # lets make an isprime array, so that in the end # isprime[n] iff n+origmax+1 is prime # so basically, isprime indexes 0..origmax^2 with 0..origmax removed rngbound = origmax * (origmax - 1) isprime = [True for n in xrange(0,rngbound)] for p in smpr: st = 0 # starting index to start cross things out at if(not ((origmax + 1) % p == 0)): st = p*(((origmax+1)//p)+1) - (origmax+1) for idx in xrange(st,rngbound, p): isprime[idx] = False # remember to shift things to their appropriate value, before returning return map(lambda x:x+origmax+1, [n for n in xrange(0,rngbound) if isprime[n]]) [/sourcecode] So, I throw these two functions together, by first running "sieveto" on some chosen bound, and then continually "extend"ing the result until I've got a list that's big enough to tell me the $latex n$-th prime. [sourcecode language="python"] goal,smallb = 100, 50 primes = sieveto(smallb) while(len(primes) < goal): primes += extend(primes) print primes[goal-1] [/sourcecode] The real trick with this idea is picking the correct value for "smallb". The best value would be something just a little bit bigger than the prime we're looking for. Alternatively, something around the square root of the prime, or the fourth root, etc. The difference can be drastic. If our goal is to find the 100th prime, starting with a "smallb" of 28 takes about 350 times as long as starting with 29 (admittedly, we're talking less than half a second on my computer, even in the worst case). The reason is that starting with 28, we run "extend" twice, ending up with a list of the first 23925 primes (when we only needed the first 100). Running, instead, with "smallb" as 29, we only extend once, generating a list of 146 primes. Of course, starting with "smallb" at 541, we don't need to extend at all, since 541 is the 100th prime.

Tags: , ,

3 Responses to “Problem 7 – The nth prime”

  1. Jaime Says:

    And you can always introduce some mathematical savvy into your guess for the upper bound, making the extension unnecessary…

  2. sumidiot Says:

    Yeah, I had thought about using the x/ln(x) approximation for pi(x). Of course, this problem is asking for something like pi^{-1}(x), and x/ln(x) is not an easy function to invert. So my initial reaction was not to try. But now that you’ve pushed me, I realize that you just pick a guess, see if guess/ln(guess) is big enough, and if not, guess bigger (say… double your guess?). Pretty easy to code that up.

    I guess what I should take from this, and my experience in problem 8, is to try the things I initially discard (x/ln(x) here, division in problem 8).

    Thanks for the push.

  3. David Says:

    I just stumbled across this post while looking for a way to approximate the nth prime, where the approximation will always be greater to or equal to the actual nth prime. That is, an upper bound to the nth prime. (I also want it for the same reason as you – to plug into a sieve to generate the first n primes.)

    I found this short article in The College Mathematics Journal:
    http://mathdl.maa.org/images/cms_upload/jaroma03200545640.pdf

    It gives the result that for n >= 7022:

    p_n <= n ln n + n (ln ln n – 0.9385)

    For my own code, I used the following upper bound for 6 <= n < 7022

    p_n <= n ln n + n ln ln n

    (From http://en.wikipedia.org/wiki/Prime_number_theorem#Approximations_for_the_nth_prime_number )

    For the first 5 primes I just used an array.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s


%d bloggers like this: