NumPy Project Euler Problem 3

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))[0]

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))[0]

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))[0]
 
#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))[0]
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.