Finding Primes
October 9th, 2006
I’ve been getting into Project Euler a little bit more recently. I found this puzzle really interesting:
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?
My Original Solver
To my knowledge, there isn’t anyway of telling if a number is prime except from dividing it by all the numbers between 1 and the number and seeing if they are factors. So I implemented the following in Python:
# Python Script to find 1001th prime number# Original Method import datetimestart = datetime.datetime.utcnow() print ‘\nOriginal Method’ i = 0n = 1 while i < 1001: n = n + 1 fail = 0 for j in range(2,n): if n % j == 0: fail = 1 if fail == 0: i = i + 1 #print str(i) + ‘. ‘ + str(n) if i == 1001: print ‘The 1001st prime number is ‘ + str(n) print ‘Time Taken: ‘+str(datetime.datetime.utcnow() - start)
Note the script above calculates the 1001st prime rather than the 10001st prime as the challenge states. When I used a similar script to find the 10,001st prime it took quite a long time on my computer so using the 1,001st to benchmark makes it a lot easier.
Unscientific testing gives the following output after 3 runs (around 21 seconds to solve):
Original MethodThe 1001st prime number is 7927Time Taken: 0:00:18.344000 Original MethodThe 1001st prime number is 7927Time Taken: 0:00:21.047000 Original MethodThe 1001st prime number is 7927Time Taken: 0:00:25.078000
Stopping
I noticed in the above code that even when the script finds a factor of a number, it’d still continue searching every number above that to see if they were factors. I put in a check:
for j in range(2,n): if fail == 0: if n % j == 0: fail = 1
This gave the following output after 3 runs:
Original Method + Factor StopThe 1001st prime number is 7927Time Taken: 0:00:09.547000 Original Method + Factor StopThe 1001st prime number is 7927Time Taken: 0:00:11.813000 Original Method + Factor StopThe 1001st prime number is 7927Time Taken: 0:00:13.750000
This halved the time it took to find the 1001st prime number.
Square Roots
With the remaining prime numbers, the script was still searching every single number between 1 and that number in order to determine whether it was a factor. For example with 13, there is no point in checking whether 8, 9, 10, 11 or 12 are factors because they would have to be multiplied by a number between 1 and 2 to get 13 and obviously there are no integers between 1 and 2.
After thinking about it, I realised that there was no point in searching any of the numbers greater than the square root of the number. With the number 12, this means you only have to search up to 3.46 (ceil to 4). Sure, 6 is a factor but you need to multiply it by 2 to make 12. Any factor above the square root of the number will have to be multiplied by a factor smaller than the square root to make the number.
I changed the script so it’d only search numbers up to the square root:
n = n + 1 fail = 0 for j in range(2,(n**0.5)+1): if fail == 0:
This had a drastic effect on the time required:
Stop at square root methodThe 1001st prime number is 7927Time Taken: 0:00:00.313000 Stop at square root methodThe 1001st prime number is 7927Time Taken: 0:00:00.516000 Stop at square root methodThe 1001st prime number is 7927Time Taken: 0:00:00.453000
From 21 seconds to 0.4 seconds - a speed up of 5000%.
The square root method managed to calculate the 10,001st prime in 11 seconds - a far cry from the 20 minutes or so required using the original script.
- Maths , Programming
- Comments(7)
Couldn’t you just find it for 1001 and times it by 10?
hahaha I crack myself up!
If you want to be a bit smarter, only use primes-up-to-sqrt(M) as denominators, too.
Note that this leads to two mutually recursive functions, but runs hellish fast because of the way primes thin-out over increasing N.
http://pig.sty.nu/temp/isprime-q.txt
(Oh yeah. I’m using a hash in there of previously-computed primes; and I print out true or false for 10 numbers either side of N. Benchmark it! Let there be gnuplot graphs!)
Interesting Tim
Only checking existing primes sounds good; I’m gonna try implementing that. I also got told in IRC that all primes above 2 and 3 are either 1 greater or smaller than a multiple of 6.
You can floor the square root, rather than ceil it, which will save (a very tiny amount of) time.
You can see the 6n±1 by doing a Sieve of Erastothenes six wide. So that would winnow the search pretty well.
The other speedup is that, as Tim said, every natural number greater than one has a unique prime factorisation, so you can do the modulus of the primes you’ve already found, 2, 3, 5, … up to √n. You don’t actually have to divide by either 2 or three if you’ve already selected only numbers matching 6n±1. This is part of the fundamental theorem of arithmetic.
I’m not sure if already mentioned but wouldn’t save it quite some time only checking every 2nd number (by changing n = n + 1 to n = n + 2)? Except of 2, every single prime is odd - of course, otherwise it’s dividable by 2 hence it can’t be a prime.
I don’t have phyton installed, would be nice to know if it really saves that much time since it fails anyway after dividing by 2. Might only make a difference if you go for really big numbers…
Greetings
xeen
Well, yes they are: in order for a number not to be divisible by either 2 or 3, it must not be divisible by 2*3; if you consider n%6, then the n for which n%6 ==2 or ==4 are obviously divisible by 2, and n%6 ==3 => it’s divisible by 3, so you’re left with the n for which n%6 is 1 or 5.
This is tantamount to saying that 2 & 3 are the two most-frequent denominators that can "catch out" a composite the quickest when you start counting from 2. Not exactly surprising and so I don’t expect that peeling-off the "is n%6 == 1 or 5?" test is going to save you anything over and above simply starting from 2 and 3 and 5… :)
However, you can scale this further. By the time you’ve examined the number 2, you know every second number is out of the running (duh). By the time you’ve examined the number 3, you know every 3rd number is out of the running. By the time you’ve examined the number 5, …. You’re now looking at http://en.wikipedia.org/wiki/Sieve_of_eratosthenes .
Now, if you were to implement the Sieve inefficiently, you’d set aside a huge array to store [1..N] in, and mark boxes prime, composite, or unknown as you count up from 2. This is where the hash (dict) in my script comes in: as it goes along computing isprime 4, isprime 5, isprime 6, … etc, it memoizes the results for later. So you can loop your n from 2..N happily knowing that, while each primality test will involve a whole bunch more tests-for-primes-in-the-range-2..sqrt(n), as time goes on those will become simple instantaneous lookups rather than calculations.
If you want it to look pretty, plot http://en.wikipedia.org/wiki/Ulam’s_spiral as you go along :)
Hi,
The following statement caught my attention: "To my knowledge, there isn’t anyway of telling if a number is prime except from dividing it by all the numbers between 1 and the number and seeing if they are factors."
There is actually a way, described at site http://www.primestructure.com. This does not work for all integers, but it does work for infinitely many ranges of integers.
Marc