2
\$\begingroup\$

Problem description:

By listing the first six prime numbers: 2, 3, 5, 7, 11, and 13, we can see that the 6th prime is 13.

What is the 10001st prime number?

Prime number:

A prime number is a whole number greater than 1 whose only factors are 1 and itself. A factor is a whole number that can be divided evenly into another number.

Example:

The first 25 prime numbers (all the prime numbers less than 100) are:

2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97

The best solution is with Sieve Of Eratosthenes

Sieve Of Eratosthenes:

In mathematics, the sieve of Eratosthenes is an ancient algorithm for finding all prime numbers up to any given limit.

It does so by iteratively marking as composite (i.e., not prime) the multiples of each prime, starting with the first prime number, 2.

The multiples of a given prime are generated as a sequence of numbers starting from that prime, with constant difference between them that is equal to that prime.

This is the sieve's key distinction from using trial division to sequentially test each candidate number for divisibility by each prime.

Once all the multiples of each discovered prime have been marked as composites, the remaining unmarked numbers are primes.

The sieve of Eratosthenes is one of the most efficient ways to find all primes smaller than n when n is smaller than 10 million or so.

Sieve Of Eratosthenes(Step by Step):

1- Create a list of consecutive integers from 2 through n: (2, 3, 4, ..., n).

2- Initially, let p equal 2, the smallest prime number.

3- Enumerate the multiples of p by counting in increments of p from 2p to n, and mark them in the list (these will be 2p, 3p, 4p, ...; the p itself should not be marked).

4- Find the smallest number in the list greater than p that is not marked. If there was no such number, stop. Otherwise, let p now equal this new number (which is the next prime), and repeat from step 3.

5- When the algorithm terminates, the numbers remaining not marked in the list are all the primes below n.

My Solution

This is my solution for problem 7 of Project Euler using Python:

 def sieve_of_eratosthenes(n): 
  '''method for finding all primes up to 
     a given natural n.
  '''
    is_prime = [True]*n
    is_prime[0] = False
    is_prime[1] = False
    
    for i in range(2,int(math.sqrt(n)+1)):  
        index = i*2
        while index < n:
            is_prime[index] = False
            index = index+i
    prime = []
    for i in range(n):
        if is_prime[i] == True:
            prime.append(i)
    return prime

How could my code be improved?

\$\endgroup\$
7
  • 5
    \$\begingroup\$ Where does your program print the 10001st prime number? \$\endgroup\$
    – Martin R
    May 29 at 18:53
  • \$\begingroup\$ @MartinR - This is the sieve_of_sratosthenes function. print the 10001st prime number out of function by: nth_prime = 10001 Max_bound = nth_primemath.log(nth_prime) + nth_primemath.log(math.log(nth_prime)) print ("the 10001st prime number is:", Sieve_of_Eratosthenes(int(Max_bound))[10000]) \$\endgroup\$ May 30 at 5:49
  • 1
    \$\begingroup\$ How do you choose n ? It is not 10,001. There are gaps in the list of primes -- you need to terminate when length(prime) reaches 10,001, not when n does. \$\endgroup\$ May 30 at 14:50
  • \$\begingroup\$ @Paul_Pedant - After completing the function: nth_prime = 10001 \$\endgroup\$ May 31 at 11:04
  • \$\begingroup\$ But 10001 is NOT the 10001st prime! It is not even PRIME - it is 73 * 137. 7 is not the seventh prime. 2, 3, 5, 7. Seven is the fourth prime. The 10001st prime is 104743. \$\endgroup\$ May 31 at 14:39

3 Answers 3

7
\$\begingroup\$

Initial improvements

Don't use int(math.sqrt(n)). Python's integers can be very large, and math.sqrt(n) needs to first convert the integer to a floating point value before the square-root is calculated, which can lead to inaccurate results when n exceeds \$2^{53}\$. While you are likely nowhere near that limit, the function math.isqrt(n) (introduced in Python 3.8) avoids the conversion to floating point, and is actually shorter to type.

Your while loop ...

        index = i*2
        while index < n:
            is_prime[index] = False
            index = index+i

... can be replaced with a for loop using Python's range object:

        for index in range(i * 2, n, i):
            is_prime[index] = False

This should be faster than your while loop, due to the looping condition and increment all handled internally in the range() object.

The final for loop used to collect primes could be replaced with list comprehension using the enumerate function:

    primes = [i for i, flag in enumerate(is_prime) if flag]

Functional improvements

2 is the only even prime number. It is worth while handling it separately, and only testing odd numbers in the sieve.

Since only odd primes are being considered in the sieve, you can double your "crossing out" increment to 2p.

Also, crossing out multiples 2p, 3p, 4p, ... is inefficient, since all odd multiples below \$p^2\$ will have been crossed out during prior steps.

As pointed out by vnp, you're missing "step 4". You don't need to cross out multiples of composite numbers, since they've already been crossed out.

    is_prime = [False, True] * (n // 2)
    is_prime[1] = False
    is_prime[2] = True

    for i in range(3, math.isqrt(n) + 1, 2):
        if is_prime[i]:
            for index in range(i * i, n, 2 * i):
                is_prime[index] = False

Tool changes

https://pypi.org/project/bitarray/

There is a module called bitarray which can be installed. Instead of using a list to store the is_prime flags, the bitarray can store the flags packed into individual bits of one long buffer. This is a phenomenal memory optimization.

Moreover, using bitarray further simplifies the crossing-out of multiples to a slice assignment operation:

import bitarray

...

    is_prime = bitarray.bitarray(n)
    is_prime.setall(False)
    is_prime[2] = True
    is_prime[3::2] = True

    ...

    for i in range(3, math.isqrt(n) + 1, 2):
        if is_prime[i]:
            is_prime[i*i::2*i] = False

    ...

\$\endgroup\$
3
  • \$\begingroup\$ Is it not even better, from a performance perspective, to swap out for i in range(3, math.isqrt(n) + 1, 2): with a manually initialised and incremented i, and letting while i**2 <= n: do the looping? \$\endgroup\$
    – Arthur
    May 30 at 10:13
  • \$\begingroup\$ I think more descriptive identifiers would increase the comprehension readability: primes = [number for number, number_is_prime in enumerate(is_prime) if number_is_prime] \$\endgroup\$ May 30 at 13:16
  • 1
    \$\begingroup\$ @Arthur Based on my timing, your while i*i <=n: i += 2 loop is 13 times slower than my for i in range(3, math.isqrt(n) + 1, 2): pass loop. One square-root is way faster than thousands of \$i^2\$ operations, and since Python is an interpreted language, the overhead of rolling-your-own loop with while and an increment statement will never be as efficient as Python's dedicated range() object. \$\endgroup\$
    – AJNeufeld
    May 31 at 20:32
3
\$\begingroup\$

The code does not implement Step 4. Before going into elimination loop, you have to test is_prime[index]. It it is False, don't bother.

Step 3 is suboptimal. Every multiple of \$p\$ less than \$p^2\$ has a factor less than \$p\$, and therefore was already marked as non-prime. It is safe to start elimination directly from \$p^2\$.

This observation also allows to remove a call to math.sqrt, which looks very out of place in this problem.

\$\endgroup\$
1
  • \$\begingroup\$ How can my code be improved with your idea? \$\endgroup\$ May 30 at 6:02
2
\$\begingroup\$

Your code is clear in intent and structure. Some suggestions on performance:

  • only keep an index of odd numbers. This halfs memory footprint and looping.
  • start applying the sieve (is_prime[index] = False) from index * index. All other multiples of index are already False (for example 7 * index is already sieved by 7.)
  • the loop at the end producing the primes seems a bit clunky. 4 lines for returning a filter on an index/range. Could be something like: primes = [index for index in range(n) if is_prime[index]]

Finally: you could use type hinting… That seems to be the current standard in Python.

\$\endgroup\$
1
  • \$\begingroup\$ How can my code be improved with use type hinting? \$\endgroup\$ May 30 at 6:03

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.