Project Euler Problem 3 seems almost impossible to crack. However, using the right algorithm – Fermat’s factorization method and NumPy, it becomes very easy. The idea is to factor a number N into two numbers c and d.
N = c * d = (a + b) (a - b) = a ** 2 - b ** 2
We can apply the factorization recursively, until we get the required prime factors.
1. Create array of trial values
The algorithm requires us to try a number of trial values for a. It makes sense to create a NumPy array and eliminate the need for loops. However, you should be careful to not create an array that is too big. On my system an array of a million elements seems to be just the right size.
a = numpy.ceil(numpy.sqrt(n)) lim = min(n, LIM) a = numpy.arange(a, a + lim) b2 = a ** 2 - n
2. Get the fractional part of the b array
We are now supposed to check whether b is a square. Use the NumPy modf function to get the fractional part of the b array.
fractions = numpy.modf(numpy.sqrt(b2))
3. Find 0 fractions
Call the NumPy where function to find the indices of 0 fractions.
indices = numpy.where(fractions == 0)
4. Find the first occurrence of a 0 fraction
Actually we only need the first occurrence of a 0 fraction. First, call the NumPy take function with the indices array from the previous step to get the a values of 0 fractions. Now we need to “flatten” this array with the NumPy ravel function.
a = numpy.ravel(numpy.take(a, indices))
Below is the entire code needed to solve the problem.
import numpy #The prime factors of 13195 are 5, 7, 13 and 29. #What is the largest prime factor of the number 600851475143 ? N = 600851475143 LIM = 10 ** 6 def factor(n): #1. Create array of trial values a = numpy.ceil(numpy.sqrt(n)) lim = min(n, LIM) a = numpy.arange(a, a + lim) b2 = a ** 2 - n #2. Check whether b is a square fractions = numpy.modf(numpy.sqrt(b2)) #3. Find 0 fractions indices = numpy.where(fractions == 0) #4. Find the first occurrence of a 0 fraction a = numpy.ravel(numpy.take(a, indices)) a = int(a) b = numpy.sqrt(a ** 2 - n) b = int(b) c = a + b d = a - b if c == 1 or d == 1: return print c, d factor(c) factor(d) factor(N)
The code prints the following output:
1234169 486847 1471 839 6857 71
Have a go
I would suggest writing an unit test to check the outcome. Use for example trial division. To get you started here is some of the code of a possible test function. Please fill in the dots.
def is_prime(n): a = numpy.arange(2, n) return len(a[...]) == 0
Also maybe you should try finetuning the size of the trial values array, but be careful.
If you liked this post and are interested in NumPy check out NumPy Beginner’s Guide by yours truly.