10
\$\begingroup\$

Previous question:

Checking for the nth prime number between 1 and 1,000,000, must compute 1,000,000th in under 3 minutes

My mission is to find the nth prime number, based on user input 1-1,000,000. It also must be able to compute the 1,000,000th prime in under 3 minutes. My code is relatively fast, however I need to make it calculate the 1,000,000th prime in 180 seconds or less. As of now it is taking 430 seconds, although before heavy updating last night it was taking about 20 minutes.

I am trying to accomplish this without importing any libraries (time's just for record keeping) and without making use of the sieve of Eratosthenes. This is an assignment that I have been given, and I have been instructed to not use the help of either of those options.

import time
goal = int(input("Enter a number n, 1-1,000,000 to find the nth prime number."))
start = time.time()
if goal > 4:
    lst = [2, 3, 5, 7]
    chk = 11
    primes = 4
    while primes < goal:
        j = chk**.5
        for n in lst:
            if n > j:
                lst.append(chk)
                primes += 1
                break
            if chk % n == 0:
                break
        chk += 2

else:
    if goal == 1:
        print("2")
    if goal == 2:
        print("3")
    if goal == 3:
        print("5")
    if goal == 4:
        print("7")
print("Prime", str(goal), "is", str(chk - 2))
elapsed = (time.time() - start)
print(elapsed, "\n")
\$\endgroup\$
  • \$\begingroup\$ Time is of course relative. On my box your code takes 112 seconds to compute p[1000000]=15485863 :) \$\endgroup\$ – Hagen von Eitzen Oct 14 '14 at 11:36
  • 2
    \$\begingroup\$ Using pypy instead of cpython gives a 10x performance increase for me on this code. Although I think that's not what you wanted. Exchanging the if-clauses in the for loop should give a tiny improvement, since the chk % n == 0 clause will be true more often, saving the time for the other comparison. \$\endgroup\$ – Mathias Oct 14 '14 at 12:36
  • 3
    \$\begingroup\$ I'm not sure what your goal is with not importing libraries, but if you want a fast way to get the nth prime of the first million primes, just use a lookup table with 1,000,000 entries. \$\endgroup\$ – Neil Oct 14 '14 at 17:25
  • 1
    \$\begingroup\$ I suppose the tongue-in-cheek solution of "use C" is out of the question? The point is, running in 3 minutes is quite an arbitrary requirement. \$\endgroup\$ – David Z Oct 14 '14 at 23:21
  • 2
    \$\begingroup\$ The goal of finding a result in 3 minutes or less makes no sense to me. You could achieve this with the existing code by just running it on a faster computer. Conversely, running the same program on a very slow computer might make the 3 minutes physically unachievable. \$\endgroup\$ – barbecue Oct 15 '14 at 1:53
17
\$\begingroup\$

Trial division is entirely the wrong algorithm to use when you want to find many primes. The Sieve of Eratosthenes is much more efficient. It seems that you know that too, but have chosen a slower algorithm anyway, and now you are unhappy that it exceeds your expected running time.

Division is slower than multiplication and addition. Square root are even more work to compute. The Sieve involves neither of those expensive operations.

Nevertheless, I'll play along. If you want to do trial division, you could at least write the code better:

  • Eliminate special cases. You have hard-coded 2, 3, 5, 7 twice. The output format for small primes is inconsistent as well.
  • Eliminate redundant bookkeeping. Since primes is just the length of lst, why not use len(lst)? That frees up the variable name primes, which would be a more appropriate name for the list itself.
  • Name your variables better. Use candidate instead of chk. Since n is always prime, call it p instead.
import time
goal = int(input("Enter a number n, 1-1,000,000 to find the nth prime number."))
start = time.time()

primes = [2, 3, 5, 7]
candidate = 11
while goal >= len(primes):
    limit = candidate ** 0.5
    for p in primes:
        if p > limit:
            primes.append(candidate)
            break
        if candidate % p == 0:
            break
    candidate += 2

print("Prime", str(goal), "is", str(primes[goal - 1]))
elapsed = (time.time() - start)
print(elapsed, "\n")

This isn't going to change the performance, but at least you would have more readable code, and less of it to manage.

\$\endgroup\$
11
\$\begingroup\$

I don't understand the rationale for not following Eratosthenes approach (which you do anyway); in any case incrementing chk by 6 (like Eratosthenes did) instead of 2 will give you a desired boost.

Code review relevant notes:

Naming. chk and lst are not very meaningful names. Consider renaming.

Code structure. If the code written cannot be reused, it is a wasted effort. Don't be shy of functions. Just give your algorithm a meaningful name.

\$\endgroup\$
  • \$\begingroup\$ Could you expand a little on incrementing chk by 6, and how that would work? It causes many numbers to not be marked as prime and reports the 1000th prime as 12105. \$\endgroup\$ – DrakkorNoir Oct 14 '14 at 6:26
  • \$\begingroup\$ @DrakkorNoir If you increment by two you have to check for numbers that are 0 or 1 modulo two. The 0 case is trivial and can be omitted. If increment by 6 you have to check for numbers that are 0,1,2,3,4, or 5 modulo 6. The 0,2,3 and 4 cases are trivial and can be omitted. You are left with 1 and 5. So when you start with chk=11 you need to check chk and chk+2 on each pass. \$\endgroup\$ – Taemyr Oct 14 '14 at 8:02
  • 1
    \$\begingroup\$ @DrakkoNoir In both cases the benefit is that at each step you need to do less trial divisions. So when you increment by two you should start your trial division with 3 - No number that you hit will be divisible by 2, and when you increment by 6 you should start your trial division on 5 - No number that you check should be divisible by 2 or 3. \$\endgroup\$ – Taemyr Oct 14 '14 at 11:55
  • \$\begingroup\$ Could I see a working example of this? I'm trying to follow the logic, and I know it should be simple. \$\endgroup\$ – DrakkorNoir Oct 14 '14 at 17:49
  • \$\begingroup\$ chk and lst are reserved on Linux console I think, not sure if that is relevant to the code though, if they are input parameters or something like that it might be worth noting... \$\endgroup\$ – Malachi Oct 14 '14 at 19:02
6
\$\begingroup\$

In addition to the many excellent points by others, avoid this kind of code:

if goal == 1:
    print("2")
if goal == 2:
    print("3")

The two conditions are mutually exclusive: if goal is 1, it cannot be 2, so in that case it shouldn't be checked for 2. Use an elif in such situations, to make sure that if the first condition was true, the 2nd one isn't evaluated unnecessarily:

if goal == 1:
    print("2")
elif goal == 2:
    print("3")

Note that this is a different situation here:

for n in lst:
    if n > j:
        # ...
        break
    if chk % n == 0:
        # ...

Here it's fine to have two ifs, because when the first one is true, thanks to the break statement in it the second one won't be executed.

\$\endgroup\$
6
\$\begingroup\$

Various other answers have suggested alternative solutions (the Sieve method is by far the most efficient, but that is not the question here). I would like to clarify @vnp's solution as it is the most relevant.

All prime numbers must be of the form \$6n+1\$ or \$6n-1\$. This is because in any list of 6 consecutive integers, only 6n-1 and 6n+1 are not divisible by 2 or 3:

\$6n-2 % 2 = 0\$, \$6n-1\$, \$6n % 6 = 0\$, \$6n+1\$, \$6n+2 % 2 = 0\$, \$6n+3 % 3 = 0\$

This is equivalent to applying the first two levels of the Sieve method, but it is conceptually simple enough to apply without having to understand the Sieve algorithm.

The resulting code is then:

def is_prime(chk, lst):
  # lst must not be an empty array
  # lst[-1] must be larger than chk**.5
  j = chk**.5
  for n in lst:
    if n > j:
      return True
    if chk % n == 0:
      return False

if goal == 1:
  result = 2
elif goal == 2:
  result = 3
elif goal == 3:
  result = 5
elif goal == 4:
  result = 7
else:             # goal > 4
    lst = [5, 7]
    # we use chk=6n, because ALL primes must be 6n-1 or 6n+1:
    chk = 6       # the number we have just checked
    primes = 4    # the last number in lst, lst[-1]=7, is the 4th prime
    while primes < goal:
      chk += 6
      chk1 = chk - 1    # check 6n-1
      if is_prime(chk1, lst):
        lst.append(chk1)
        primes += 1
        if primes==goal:
          result = chk1
          break
      chk2 = chk + 1    # check 6n+1
      if is_prime(chk2, lst):
        lst.append(chk2)
        primes += 1
        if primes==goal:
          result = chk2
          break

The timing, input and output has been omitted.

This results in a speed-up by a factor of ~1.5 from your original code.

  • Goal=100,000: 5.7s vs 3.6s
  • Goal=200,000: 16s vs 10s
  • Goal=500,000: 61s vs 40s
  • Goal=1,000,000: 154s vs 96s
\$\endgroup\$
  • \$\begingroup\$ You can't post code in a comment anyway, so an answer is better. As an FYI, you will be able to comment on your own posts, even without reputation. \$\endgroup\$ – Brythan Oct 15 '14 at 3:41
3
\$\begingroup\$

Here you have a complete rewrite of your code with my functional style:

  • I created a function called called is_prime and a function called handle_input, you can reuse them if you wish.

  • I put the code that "displays stuff" under a __name__ == "__main__" condition so that if you want you can import this as a module and run only the parts you like...

  • I put a documentation string at the start of the functions so that if you forget what they do you can just read that string and remember.

Why I use documentation

Documentation helps you write maintainable code. Here I explain: very simple and basic code can be written in any way:

isprime = lambda x: not [i for i in range(2,x) if x%i == 0] # Very slow but works

and if you look at it closely you will understand how this one-liner works. This has some use, when I want to write in the interactive shell a prime-checking function in seconds I write this. Anyway in anything that is going to remain in my computer I would use the is_prime function I defined below because it is much simpler and straitghforward, and I would add a line of documentation on it (in such a simple case documentation may be redundant but I just like it).

  • Now imagine half a year is passed what would you rather see 1 or 2?:

1)

isprime = lambda x: not [i for i in range(2,x) if x%i == 0]
# The programme continues ...

2)

def is_prime(n):
    """
    Computes if n is prime using optimazed trial division.
    """
    if n % 2 == 0 and not n == 2:
        return False
    else:
        for i in range(3,int(n**0.5+1),2):
            if n % i == 0:
                return False
    return True
# The program continues

Why I use functions

You can use a piece of your old program in a new program

I can define a shift function:

def shift(letter,value):
    """
    Shifts a letter in the alphabet by the value,
    if the alphabet ends go back to the start.
    >>> shift('a',5)
    f
    >>> "".join([shift(i,20) for i in "hello"])
    'byffi'
    """
    current_letter_value = ALPHABET.find(letter)
    end_value = current_letter_value + value
    return cycle_get(ALPHABET,end_value)

OR

I can inline the code

#some code
    current_letter_value = ALPHABET.find(letter)
    end_value = current_letter_value + value
    end_value = cycle_get(ALPHABET,end_value)
#some code

I prefer using the function because this function is useful both in a Caesar cipher and in a polyalphabetic cipher, thus allowing me to reuse code, not bore myself and go straitgh to what is new. Also if I forget about what it does (see: Why I use documentation) I can just read the documentation string and remember

Anyway

If you prefer imperative style programming (i.e. no functions just statments) go ahead and you will be fine as long as the program is not huge.

It follows a re-write of the program following my style, you may or may not like it, see the above to understand why I like coding this way.

def is_prime(n):
    """
    Computes if n is prime using optimazed trial division.
    """
    if n % 2 == 0 and not n == 2:
        return False
    else:
        for i in range(3,int(n**0.5+1),2):
            if n % i == 0:
                return False
    return True

def handle_input(primes):
    """
    Provides a user interface to get the n-th prime.
    -1 is because non programmers users will use 1 indexing
    while computers use 0 indexing: -1 converts from 1 indexing to
    0 indexing.
    """
    index = int(input("What prime number do you want? "))
    print(primes[index-1])
    print() # Empty line for better interface

if __name__ == "__main__":
    print("Wait some seconds for initialization...")
    primes = [i for i in range(2,16*10**6) if is_prime(i)]
    while 1:
        handle_input(primes)

Weirdly this programme takes 8 minutes to initialize on my machine. I expected it to be faster.

\$\endgroup\$
  • \$\begingroup\$ You're supposed to look not for primes below one million, but the one million first primes. Therefore we need primes to be initialized to [i for i in range(2, 16*10**6) if is_prime(i)]. Shouldn't take you from 9 seconds to 3 minutes, though. \$\endgroup\$ – Arthur Oct 14 '14 at 9:03
3
\$\begingroup\$
  • Since you use chk += 2 to skip over even numbers, you don't need to test division by 2. You can take 2 off lst.
  • You could do the same with 3 if you increment chk by 4 every other time to skip the odd multiples of 3.
  • Convert j to int before the inner loop. Otherwise the n > j comparison has to convert n to float every time.
  • Better yet, avoid the n > j comparison in the inner loop by maintaining the list of possible divisors (odd primes up to the square root of current candidate) separately from the list of primes found so far.
  • Put the code in a function because local variables are faster in Python.
\$\endgroup\$
1
\$\begingroup\$

Your code is a sieve analogous to the sieve of Eratosthenes, only much more poorly designed. For the record, the following standard sieve finds the nth prime in less than 2 seconds on my machine:

public static long primeNth( int n ){
    boolean[] sieve = new boolean[100000000];
    int prime = 2;
    int ctPrimes = 0;
    while( true ){
        ctPrimes++;
        if( ctPrimes == n ) return prime;
        int multiple = 1;
        while( true ){
            int composite = multiple * prime;
            if( composite >= 100000000 ) break;
            sieve[composite] = true;
            multiple++;
        }
        prime++;
        while( sieve[prime] ) prime++; 
    }
}

I would recommend using more explanatory variable names.

  • Appending to an array is expensive; it is better to pre-allocate the array to the size needed.

  • Square roots are expensive.

Since there is no way to predict a prime number, you have to find every single prime to N. Therefore, there is really no option other than to use a sieve, except simply trying every single number which is obviously a lot less efficient than a sieve.

\$\endgroup\$
  • \$\begingroup\$ If 100 million isn't enough, you can make better bit effenciency (for a small, but perhaps measurable cost) by using a bit set instead of a boolean array. However, if 100m is enough, you might as well use the boolean array since it will be faster. \$\endgroup\$ – corsiKa Oct 14 '14 at 23:39
1
\$\begingroup\$

The best way to improve the performance by far:

change the limit number to sqrt(candidate) .

Assuming that the number is not prime, there being another number that divides it means that a / b = c. If there is no integer sqrt then we can conclude that one of them (b or c) must be a root.

Don't forget to convert the sqrt(candidate) to int.

\$\endgroup\$
  • 7
    \$\begingroup\$ He already does that: j = chk**.5 \$\endgroup\$ – Mathias Oct 14 '14 at 12:32
0
\$\begingroup\$

Even if you are not allowed to use the Sieve of Eratosthenes, you might be able to get away with sneaking in another kind of sieve.

For example, the linear-time algorithm for computing a table of Euler's totient function indirectly finds all primes. In its outer loop a prime will find its table slot empty whereas a composite will already have a value. It is several orders of magnitude more efficient than trial division and it is most definitely not based on the Sieve of Eratosthenes (or some half-assed inferior reinvention of it, like the Sieve of Sundaram). See Linear time Euler's Totient Function calculation.

The topic Calculate the number of primes up to n (a speed challenge on the Programming Puzzles & Code Golf site) contains a wealth of ideas, information and working code that can help you find the nth prime with record-breaking speed, if you are so inclined.

\$\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.