Prime Number Generation

If I asked you right now: what is the 10001st prime number?  It would be very hard to come up with an answer on the spot. (this is question #7 on Project Euler). However, if one did want to find the answer, one method to use would be simple “trial division”.


The idea behind trial division is simple: since a prime must not have any factors other than itself and 1, just divide all the numbers from 2 to N-1 from the number (call it N) to be tested. If the division carries out and has no remainders, then the divided number must be a factor of the number we are checking and so the number we are checking is not prime.

Simple right? However, think of how many operations that would take to test a single number: one would need N-2 operations to loop from 2 to N-1, N-2 operations to carry out the division, N-2 operations to check if the division was carried out, all of which contribute to a massive computation time.

Luckily, there are methods to optimize our trial division prime checker. For instance, a commonly known fact about primes is that one only needs to check up to the integer before $\sqrt{n}$ instead of N-1 (a simple proof is offered in the next post). Furthermore, if we check right off the bat if the number is divisible by 2 or 3, the we will not need to check every second and third number from 4 to $\sqrt{n}$. This is really wheel factorization with two spokes. Such a type of trial division also leads one to the realization that all primes greater than 3 can be expressed in the form $6k +/- 1$, an interesting fact proven by the great mathmatician Leonhard Euler.

We are now set to create our prime checking algorithm!

The tricks we will employ are:

1. Checking only up to $\sqrt{n}$.

2. Checking first if the number is divisible by 2 or 3

3. Starting from 5 and checking every number of the form $6k +/- 1$.

Excellent. Now we are ready to write our algorithm. I will be writing in pseudocode.

Pseudocode:
function prime?[n]

     if n mod 2 or n mod 3 == 0

           return false // if the number is divisible by 2 or 3, it is not prime (only works for numbers above 3)

           else

                limit = sqrt(n) // bounds our testing limit to under $\sqrt{n}$

               for (i = 5, i <= limit, i+= 6) // creates a FOR loop that starts at 5, stops when exceeding the limit, and increments the counter variable (i) by 6 each iteration

                       if n mod i or n mod (i+2) = 0

return false // if the function is divisible by numbers of the form $6k+/-1$, the number is not prime.

                       end if

                end for

     end if

     return true // the number has passed all the tests we have given it; it is a prime

end function prime?[n]

And there it is! Just run this function with a counter until it reaches 10001 if you want to find the 10001st prime! Of course, there is a faster way to create a list of primes up to 10001st, namely using the sieve of Eratosthenes. This very effective and elegant algorithm may be covered in future blog entries. Keep readin’ til then.

********

References:

1. Project Euler (the pseudocode for the algorithm was heavily borrowed from PE).

2. The Prime Number Glossary (for interesting and useful facts about primes).