16 Jun 2011

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.
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)
                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.
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).


  1. Up until very recently, your polynomial-time algorithm was pretty much the deterministic primality tester we had. However in 2002, Agrawal, Kayal and Saxena published this


    out of the blue, which does the job in /polylogarithmic/-time.

    The fascinating thing is that it requires no maths more recent than Euler, yet no-one found it for 200 years.

  2. @Andy Jones

    Ah yes :) I hear about the AKS class of algorithms every now and then. It's a great algorithm, but I think for most people it would be a little less intuitive to implement since it requires a Euler's totient function.

    On another note, most prime checkers in software like Mathematica are made up of many probabilistic prime checking algorithms (just a tidbit of interest :) ).