Finding the greatest prime divisor (the fastest program possible)

Question

I was checking the problems on http://projecteuler.net/

The third problem is as follows:

The prime factors of 13195 are 5, 7, 13 and 29. What is the largest prime factor of the number 600851475143 ?

My solution code is below. But it is so slow that I think it will take weeks to complete. How can I improve it? Or is Python itself too slow to solve this problem?

def IsPrime(num):
    if num < 2:
        return False
    if num == 2:
        return True
    else:
        for div in range(2,num):
            if num % div == 0:
                return False
        return True

GetInput = int (input ("Enter the number: "))


PrimeDivisors = []

for i in range(GetInput, 1, -1):
    print(i)
    if GetInput % i == 0 and IsPrime(i) is True:
        PrimeDivisors.append(i)
        break
    else:
        continue


print(PrimeDivisors)
print("The greatest prime divisor is:", max(PrimeDivisors))

For starters you only need to test up until the square root of the number that you are searching for factors from. — shuttle87, Dec 14 '13 at 09:14
[This is the fastest python implementation I have ever seen](http://code.activestate.com/recipes/117119-sieve-of-eratosthenes/) — thefourtheye, Dec 14 '13 at 09:16
@shuttle87 and you only need to test odd numbers (apart from 2). — Roger Rowland, Dec 14 '13 at 09:17
@shuttle87, Largest prime factor of 10 is 5, which is bigger than the sqrt(10) (3.16, aprox). Largest prime factor of 33 is 11, which is bigger than the sqrt(33) (5.74, aprox). I guess you're confusing this with the propriety which states that, if a number has a prime factor bigger than its sqrt, it has to have at least another one smaller than its sqrt. — Rakanoth, Dec 14 '13 at 13:34
related: [Stuck on Project Euler #3 in python](http://stackoverflow.com/q/12999706/4279). Even [brute-force approach completes in fractions of a second](http://ideone.com/fbv5V4) — jfs, Dec 15 '13 at 10:00

Blender · Accepted Answer · 2013-12-14T11:43:38.620

2

The problem with your solution is that you don't take your found prime factors into account, so you're needlessly checking factors after you've actually found the largest one. Here's my solution:

def largest_prime_factor(n):
    largest = None

    for i in range(2, n):
        while n % i == 0:
            largest = i
            n //= i

        if n == 1:
            return largest

    if n > 1:
        return n

Project Euler problems are more about mathematics than programming, so if your solution is too slow, it's probably not your language that's at fault.

Note that my solution works quickly for this specific number by chance, so it's definitely not a general solution. Faster solutions are complicated and overkill in this specific case.

edited Dec 14 '13 at 11:43

answered Dec 14 '13 at 09:51

Blender

257,973
46
399
459

The General Number Field Sieve solution linked to in this answer is definitely overkill for this problem. It works best for numbers that are at least 100 decimal digits long. – poke Dec 15 '13 at 02:33
@poke: It works for any number, especially large numbers. `msieve -v 600851475143` runs faster than my solution. – Blender Dec 15 '13 at 02:36
@Bender: My comment was framed around the fact that the quadratic sieve is better than the GNFS for numbers less than 100-130 decimal digits. – poke Dec 15 '13 at 02:55

James Mills · Answer 2 · 2013-12-14T09:41:22.577

This might not be the fastest algorithm but it's quite efficient:

def prime(x):
    if x in [0, 1]:
        return False
    for n in xrange(2, int(x ** 0.5 + 1)):
        if x % n == 0:
            return False
    return True

def primes():
    """Prime Number Generator

    Generator an infinite sequence of primes

    http://stackoverflow.com/questions/567222/simple-prime-generator-in-python
    """

    # Maps composites to primes witnessing their compositeness.
    # This is memory efficient, as the sieve is not "run forward"
    # indefinitely, but only as long as required by the current
    # number being tested.
    #
    D = {}  

    # The running integer that's checked for primeness
    q = 2  

    while True:
        if q not in D:
            # q is a new prime.
            # Yield it and mark its first multiple that isn't
            # already marked in previous iterations
            # 
            yield q        
            D[q * q] = [q]
        else:
            # q is composite. D[q] is the list of primes that
            # divide it. Since we've reached q, we no longer
            # need it in the map, but we'll mark the next 
            # multiples of its witnesses to prepare for larger
            # numbers
            # 
            for p in D[q]:
                D.setdefault(p + q, []).append(p)
            del D[q]

        q += 1

def primefactors(x):
    if x in [0, 1]:
        yield x
    elif prime(x):
        yield x
    else:
        for n in primes():
            if x % n == 0:
                yield n
                break
        for factor in primefactors(x // n):
            yield factor

Usage:

>>> list(primefactors(100))
[2, 2, 5, 5]

You didn't define `primes`. Also, it seems much simpler to use `for prime in primes` than to use that `while True`. — user2357112 supports Monica, Dec 14 '13 at 09:38
As it is, your code doesn't seem to be bugged (at least, if you consider the strange output for 0 and 1 to be okay), but it does something silly. To determine that a composite number is composite, it searches for a factor. If it finds one, it then... does another search for a factor, implemented completely differently. It also doesn't use any knowledge of previously determined factors, restarting every search from the ground up. — user2357112 supports Monica, Dec 14 '13 at 09:50
Yeah I agree the prime sieve good be improved upon. Perhaps I'll find some time to do so later :) — James Mills, Dec 14 '13 at 09:52

SzieberthAdam · Answer 3 · 2013-12-14T22:00:57.633

My code which seems enough fast to me. Using collections.defaultdict() would make the code of primes() abit cleaner but I guess the code would loose some speed due to importing it.

def primes():
    """Prime number generator."""
    n, skip = 2, {}
    while True:
        primes = skip.get(n)
        if primes:
            for p in primes:
                skip.setdefault(n + p, set()).add(p)
            del skip[n]
        else:
            yield n
            skip[n * n] = {n}
        n += 1

def un_factor(n):
    """Does an unique prime factorization on n.

    Returns an ordered tuple of (prime, prime_powers)."""
    if n == 1:
        return ()
    result = []
    for p in primes():
        (div, mod), power = divmod(n, p), 1
        while mod == 0:
            if div == 1:
                result.append((p, power))
                return tuple(result)
            n = div
            div, mod = divmod(n, p)
            if mod != 0:
                result.append((p, power))
            power += 1

Test run:

>>> un_factor(13195)
((5, 1), (7, 1), (13, 1), (29, 1))
>>> un_factor(600851475143)
((71, 1), (839, 1), (1471, 1), (6857, 1))
>>> un_factor(20)
((2, 2), (5, 1))

EDIT: Minor edits for primes() generator based on this recipe.

EDIT2: Fixed for 20.

EDIT3: Replaced greatest_prime_divisor() with un_factor().

Your code doesn't work for greatest_prime_divisor(20), for example. — Rakanoth, Dec 14 '13 at 13:35

M4rtini · Answer 4 · 2013-12-14T17:24:34.613

def getLargestFactor(n):
    maxFactor = sqrt(n)
    lastFactor = n
    while n%2 == 0:
        n  /= 2
        lastFactor = 2
    for i in xrange(3,int(maxFactor),2 ):
        if sqrt(n) < i:
             return n
        while n%i == 0 and n > 1:
            n  /= i
            lastFactor = i
    return lastFactor

This should be fairly efficient. Dividing each factor all out, this way we only find the prime factors. And using the fact that there can only be one prime factor of a number larger than sqrt(n).

Finding the greatest prime divisor (the fastest program possible)

4 Answers4

Linked