Finding prime factors

Prime factors (http://en.wikipedia.org/wiki/Prime_factor) are prime numbers that divide an integer exactly without a remainder. Finding prime factors seems almost impossible to crack. However, using the right algorithm—Fermat's factorization method (http://en.wikipedia.org/wiki/Fermat%27s_factorization_method) and NumPy—it becomes very easy. The idea is to factor a number N into two numbers c and d, according to the following equation:

Finding prime factors

We can apply the factorization recursively, until we get the required prime factors.

How to do it...

The algorithm requires us to try a number of trial values for a.

  1. Create an array of trial values.

    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 in terms of memory requirements. 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

    We used the ceil function to return the ceiling of the input, element-wise.

  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 zero fractions, where the fractional part is 0:

    indices = numpy.where(fractions == 0)
  4. Find the first occurrence of a zero fraction.

    Actually, we only need the first occurrence of a zero fraction. First, call the NumPy take function with the indices array from the previous step to get the values of zero fractions. Now we need to "flatten" this array with the NumPy ravel function:

    a = numpy.ravel(numpy.take(a, indices))[0]

The following is the entire code needed to solve the problem of finding the largest prime factor of the number 600851475143:

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 output for this code is the following:

1234169 486847
1471 839
6857 71

How it works...

We applied the Fermat factorization recursively using the NumPy functions ceil, modf, where, ravel, and take. The description of these functions is as follows:

Function

Description

ceil

Calculates the ceiling of array elements.

modf

Returns the fractional and integral part of floating point numbers.

where

Returns array indices based on condition.

ravel

Returns a flattened array.

take

Takes element from an array.

..................Content has been hidden....................

You can't read the all page of ebook, please click here login for view all page.
Reset