Will Clausen's Website

Categories
WCW

Solution: Project Euler Problem 7

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!

Leave a Reply

Your email address will not be published.