The Problem:
What is the 10 001st prime number?
My Solution (in Python):
# Author: Will Clausen
#
# Date: Jan 7, 2013
#
# This program will solve Problem 7 from Project Euler
import math
import time
# This function solves problem 7 from Project Euler.
# This function uses a primality tester implemented with
# trial division. This solution is the second fastest.
# See the other implementations below.
def problem7(limit):
primesSeen = 1
check = 3
while primesSeen < limit:
if isPrime(check, 30):
primesSeen += 1
lastPrime = check
check += 2
return lastPrime
# This function solves also solves problem 7, but
# uses the Miller-Rabin primality test to check
# if a number is prime instead of trial division.
# This solution is not as fast as either of the
# other solutions.
def problem7_2(limit):
primesSeen = 1
check = 3
while primesSeen < limit:
if probIsPrime(check, 30):
primesSeen += 1
lastPrime = check
check += 2
return lastPrime
# This solution to problem 7 uses the Prime Number Theorem
# (hidden in a couple helper functions below) and the sieve
# of Eratosthenes to quickly solve the problem. This is the
# fastest solution.
def problem7_3(limit):
LoP = firstPrimes(limit)
return LoP[-1]
# Helper function to generate primes less than an integer
# input. The sieve of Eratosthenes is the method used to
# generate the primes.
def primesLessThan(num):
# Initialize some variables
stop = (num-1)/2
LoP = [2]
oddNums = [True]*stop
i = 0
p = 3
j = num
#Begin sieving
while math.pow(p, 2)< num:
# Check for prime
if oddNums[i] == True:
LoP += [p]
j = int((math.pow(p, 2)-3))/2
#sift on p
while j < stop:
oddNums[j] = False
j += p
i += 1
p += 2
# Get the remaining primes
while i < stop:
if oddNums[i] == True:
LoP += [p]
i += 1
p += 2
return LoP
# Function that returns the first numPrimes primes in a list
def firstPrimes(numPrimes):
# Use the Prime Number Theorem here to help figure out about where the
# last prime is
# First (slightly over)estimate where the numPrimes prime is using the Prime Number Theorem.
estimate = int(numPrimes*(math.log(numPrimes) + math.log(math.log(numPrimes))))
# Get a list of the primes less than the estimate
bigList = primesLessThan(estimate)
# Cut the list at the number needed and return that
result = bigList[:numPrimes]
return result
# This function determines is an integer input is prime.
# The method used is trial division. The runtime of this
# function is bounded by (Num^(1/2)). It beats the
# Miller-Rabin primality test for numbers around less than
# 1,000,000 (in terms of runtime).
def isPrime(Num):
if (Num == 2):
return True
if (Num % 2 == 0):
return False
sqrtNum = int(math.sqrt(Num)) + 1
for x in range(3, sqrtNum, 2):
if Num % x == 0:
return False
return True
# This is a python implementation of the Miller-Rabin
# primality test. Trial division beats the Miller-Rabin
# method for numbers under 1,000,000, with numchecks = 10
# (actually a relatively small value)
def probIsPrime(num, numChecks):
# Some optimizations
if num == 2:
return True
if num % 2 == 0:
return False
# Here I use a list of the first 100 primes to use
# as bases when checking primality. This will allow for
# accurate calculations for primes of any reasonable size
# (or even an unreasonable size, like 10^100)
LoP = [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, 101, 103, 107, 109, 113,
127, 131, 137, 139, 149, 151, 157, 163, 167, 173,
179, 181, 191, 193, 197, 199, 211, 223, 227, 229,
233, 239, 241, 251, 257, 263, 269, 271, 277, 281,
283, 293, 307, 311, 313, 317, 331, 337, 347, 349,
353, 359, 367, 373, 379, 383, 389, 397, 401, 409,
419, 421, 431, 433, 439, 443, 449, 457, 461, 463,
467, 479, 487, 491, 499, 503, 509, 521, 523, 541]
i = 0
# Unfortunately, using numbers greater than num for checking
# causes errors, so some checking is necessary here to ensure
# accuracy
if (num < LoP[numChecks]):
LoP = range(2, num)
# The meat of the computation
while numChecks > 0:
check = LoP[i]
i += 1
# Be careful not to index out of the list of primes...
if (i == len(LoP)):
i = 0
if not strongPseudoprimeTest(num, check):
break
numChecks -= 1
return numChecks == 0
# Function to check primality. It effectively employs Fermat's
# Little Theorem in a creative way to deterimine if a number,
# (the check) demonstrates the compositeness of a number (the num).
def strongPseudoprimeTest(num, check):
# Initialize some helpful variables, nice names are unnecessary.
d = num - 1
s = 0
numLessOne = num - 1
# If d is even divide it by two and increment s
while d % 2 == 0:
d /= 2
s += 1
t = pow(check, d, num)
if ((t == 1) or (t == numLessOne)):
return True
s -= 1
while s > 0:
s -= 1
t = pow(t, 2, num)
if (t == numLessOne):
return True
return False
# Solution: 104743
Feedback is always appreciated!