203

What is the best approach to calculating the largest prime factor of a number?

I'm thinking the most efficient would be the following:

  1. Find lowest prime number that divides cleanly
  2. Check if result of division is prime
  3. If not, find next lowest
  4. Go to 2.

I'm basing this assumption on it being easier to calculate the small prime factors. Is this about right? What other approaches should I look into?

Edit: I've now realised that my approach is futile if there are more than 2 prime factors in play, since step 2 fails when the result is a product of two other primes, therefore a recursive algorithm is needed.

Edit again: And now I've realised that this does still work, because the last found prime number has to be the highest one, therefore any further testing of the non-prime result from step 2 would result in a smaller prime.

3
  • My approach was: (1) divide large, possible number by 2; (2) check if the large number divides evenly into it; (3) if so, check if the divided by 2 number is prime. If it is, return it. (4) Else, substract 1 from the divided by 2 number, returning to step 3. Jul 28, 2014 at 0:59
  • 1. find any number that divides clearly (for i = 2 to int(sqr(num)) ) 2. divide by that number (num = num/i) and recur until nothing is found in 1.'s interval 3. num is the largest factor Sep 11, 2015 at 21:37
  • 1
    We can Divide with small primes, and the one which is finally left, is the Largest Prime Factor (I guess)
    – user5058091
    Jul 24, 2016 at 3:05

30 Answers 30

152

Here's the best algorithm I know of (in Python)

def prime_factors(n):
    """Returns all the prime factors of a positive integer"""
    factors = []
    d = 2
    while n > 1:
        while n % d == 0:
            factors.append(d)
            n /= d
        d = d + 1

    return factors


pfs = prime_factors(1000)
largest_prime_factor = max(pfs) # The largest element in the prime factor list

The above method runs in O(n) in the worst case (when the input is a prime number).

EDIT:
Below is the O(sqrt(n)) version, as suggested in the comment. Here is the code, once more.

def prime_factors(n):
    """Returns all the prime factors of a positive integer"""
    factors = []
    d = 2
    while n > 1:
        while n % d == 0:
            factors.append(d)
            n /= d
        d = d + 1
        if d*d > n:
            if n > 1: factors.append(n)
            break
    return factors


pfs = prime_factors(1000)
largest_prime_factor = max(pfs) # The largest element in the prime factor list
17
  • 11
    Please read and/or run this code before voting it down. It works fine. Just copy and paste. As written prime_factors(1000) will return [2,2,2,5,5,5], which should be interpreted as 2^3*5^3, a.k.a. the prime factorization. Jan 7, 2009 at 17:59
  • 11
    "runs in O(sqrt(n)) in the worst case" - No, it runs in O(n) in the worst case (e.g. when n is prime.) Sep 26, 2010 at 15:23
  • 19
    Easy to make it O(sqrt(n)), you just stop the loop when d*d > n, and if n > 1 at this point then its value should be appended to the list of prime factors. Mar 19, 2012 at 4:53
  • 5
    Is there a name for this? Jul 19, 2013 at 23:44
  • 16
    since 2 is the only even prime number, so instead of adding 1 each time, you can iterate separately for d=2 and then increment it by 1 and then from d=3 onwards you can increment by 2. so it will decrease the number of iterations... :)
    – tailor_raj
    Nov 19, 2013 at 9:13
147

Actually there are several more efficient ways to find factors of big numbers (for smaller ones trial division works reasonably well).

One method which is very fast if the input number has two factors very close to its square root is known as Fermat factorisation. It makes use of the identity N = (a + b)(a - b) = a^2 - b^2 and is easy to understand and implement. Unfortunately it's not very fast in general.

The best known method for factoring numbers up to 100 digits long is the Quadratic sieve. As a bonus, part of the algorithm is easily done with parallel processing.

Yet another algorithm I've heard of is Pollard's Rho algorithm. It's not as efficient as the Quadratic Sieve in general but seems to be easier to implement.


Once you've decided on how to split a number into two factors, here is the fastest algorithm I can think of to find the largest prime factor of a number:

Create a priority queue which initially stores the number itself. Each iteration, you remove the highest number from the queue, and attempt to split it into two factors (not allowing 1 to be one of those factors, of course). If this step fails, the number is prime and you have your answer! Otherwise you add the two factors into the queue and repeat.

3
  • 3
    Pollard rho and the elliptic curve method are much better at getting rid of small prime factors of your number than the quadratic sieve. QS has about the same runtime no matter the number. Which approach is faster depends on what your number is; QS will crack hard-to-factor numbers faster while rho and ECM will crack easy-to-factor numbers faster.
    – tmyklebu
    Sep 22, 2014 at 17:02
  • Thank you for the Quadratic variation suggestion. I needed to implement this for one of my clients, the initial version I came up was something along the lines of what @mercutio suggested in his question. The quadratic solution is what is powering my client's tool now at math.tools/calculator/prime-factors .
    – dors
    Apr 3, 2017 at 14:34
  • If there is an efficient way of solving this problem, doesn't that mean that web encryption is not not secure?
    – BenKoshy
    Apr 23, 2020 at 23:55
18

My answer is based on Triptych's, but improves a lot on it. It is based on the fact that beyond 2 and 3, all the prime numbers are of the form 6n-1 or 6n+1.

var largestPrimeFactor;
if(n mod 2 == 0)
{
    largestPrimeFactor = 2;
    n = n / 2 while(n mod 2 == 0);
}
if(n mod 3 == 0)
{
    largestPrimeFactor = 3;
    n = n / 3 while(n mod 3 == 0);
}

multOfSix = 6;
while(multOfSix - 1 <= n)
{
    if(n mod (multOfSix - 1) == 0)
    {
        largestPrimeFactor = multOfSix - 1;
        n = n / largestPrimeFactor while(n mod largestPrimeFactor == 0);
    }

    if(n mod (multOfSix + 1) == 0)
    {
        largestPrimeFactor = multOfSix + 1;
        n = n / largestPrimeFactor while(n mod largestPrimeFactor == 0);
    }
    multOfSix += 6;
}

I recently wrote a blog article explaining how this algorithm works.

I would venture that a method in which there is no need for a test for primality (and no sieve construction) would run faster than one which does use those. If that is the case, this is probably the fastest algorithm here.

2
  • 12
    You can actually take this idea even further, e.g. beyond 2,3,5 all primes are of the form 30n+k (n >= 0) where k only takes those values between between 1 and 29 that are not divisible by 2,3 or 5, i.e. 7,11,13,17,19,23,29. You can even have this dynamically adapt after every few primes you found so far to 2*3*5*7*...*n+k where k must not be divisible by any of these primes (note that not all possible k need be prime, e.g. for 210n+k you have to include 121, otherwise you'd miss 331) Oct 7, 2013 at 11:51
  • 2
    I guess it should be while (multOfSix - 1 <= n) Jul 16, 2014 at 4:05
9

Similar to Triptych's answer but also different. In this example, a list or dictionary is not used. Code is written in Ruby

UPDATED

Thanks to a suggestion by Cooper Wolfe, the algorithm has been enhanced and optimized.

def largest_prime_factor(number)
  i = 2
  while number > 1
    if number % i == 0
      number /= i
    elsif i > Math.sqrt(number)
      i = number
    else
      i += 1
    end
  end
  return i
end

largest_prime_factor(600851475143)
# => 6857
4
  • Finally something readable and instantly (in js) executable at the same time. I was trying to use infinite prime list and it was already too slow on 1 million.
    – Ebuall
    Apr 7, 2018 at 16:39
  • Another else-if can significantly cut down on wasted compute: else if i > sqrt(number): i = number Jan 17, 2023 at 21:11
  • (Assigning i = number is a very indirect way to terminate the loop - I'd just return number.)
    – greybeard
    Oct 4, 2023 at 7:34
  • @greybeard Yeah, but in this case I prefer as-is Oct 4, 2023 at 11:11
8

JavaScript code:

'option strict';

function largestPrimeFactor(val, divisor = 2) { 
    let square = (val) => Math.pow(val, 2);

    while ((val % divisor) != 0 && square(divisor) <= val) {
        divisor++;
    }

    return square(divisor) <= val
        ? largestPrimeFactor(val / divisor, divisor)
        : val;
}

Usage Example:

let result = largestPrimeFactor(600851475143);

Here is an example of the code:

4

All numbers can be expressed as the product of primes, eg:

102 = 2 x 3 x 17
712 = 2 x 2 x 2 x 89

You can find these by simply starting at 2 and simply continuing to divide until the result isn't a multiple of your number:

712 / 2 = 356 .. 356 / 2 = 178 .. 178 / 2 = 89 .. 89 / 89 = 1

using this method you don't have to actually calculate any primes: they'll all be primes, based on the fact that you've already factorised the number as much as possible with all preceding numbers.

number = 712;
currNum = number;    // the value we'll actually be working with
for (currFactor in 2 .. number) {
    while (currNum % currFactor == 0) {
        // keep on dividing by this number until we can divide no more!
        currNum = currNum / currFactor     // reduce the currNum
    }
    if (currNum == 1) return currFactor;    // once it hits 1, we're done.
}
6
  • Yes, but this is horribly inefficient. Once you've divided out all the 2s, you really shouldn't try dividing by 4, or by 6, or ...; It really is much more efficient in the limit to only check primes, or use some toher algorithm.
    – wnoise
    Oct 28, 2008 at 5:30
  • 6
    +1 to offset wnoise, who I think is wrong. Trying to divide by 4 will only happen once, and will fail immediately. I don't think that's worse than removing 4 from some list of candidates, and it's certainly faster than finding all primes beforehand. Jan 7, 2009 at 16:15
  • 2
    @Beowulf. Try running this code before voting down. It returns prime factors; you just don't understand the algorithm. Jan 7, 2009 at 18:01
  • 3
    the code works ok, but is slow if the incoming number is a prime. I would also only run up to the square and increment by 2. It might be too slow for very big numbers, though.
    – blabla999
    Jan 13, 2009 at 19:11
  • 4
    blabla999 is exactly right. Example is 1234567898766700 = 2*2*5*5*12345678987667. When we've reached currFactor = 3513642, we know that 12345678987667 is prime, and should return it as the answer. Instead, this code will continue the enumeration until 12345678987667 itself. That's 3,513,642x slower than necessary.
    – Will Ness
    Apr 11, 2014 at 18:04
4
    //this method skips unnecessary trial divisions and makes 
    //trial division more feasible for finding large primes

    public static void main(String[] args) 
    {
        long n= 1000000000039L; //this is a large prime number 
        long i = 2L;
        int test = 0;

        while (n > 1)
        {
            while (n % i == 0)
            {
                n /= i;     
            }

            i++;

            if(i*i > n && n > 1) 
            {
                System.out.println(n); //prints n if it's prime
                test = 1;
                break;
            }
        }

        if (test == 0)  
            System.out.println(i-1); //prints n if it's the largest prime factor
    }
9
  • 1
    have you tried your code with 1,000,000,000,039? it should run in a blink of an eye too. Does it?
    – Will Ness
    Apr 12, 2014 at 11:57
  • 2
    You could know it in advance, without trying. 10^12 = (2*5)^12 = 2^12 * 5^12. So your while loop will go through i values of 2,2,2,2,2,2,2,2,2,2,2,2, 2,3,4,5, 2,3,4,5, 2,3,4,5, 2,3,4,5, 2,3,4,5, 2,3,4,5, 2,3,4,5, 2,3,4,5, 2,3,4,5, 2,3,4,5, 2,3,4,5, 2,3,4,5. All of 60 iterations. But for (10^12+39) there will be (10^12+38) iterations, i=2,3,4,5,6,...,10^12+39. Even if 10^10 ops take one second, 10^12 will take 100 seconds. But only 10^6 iterations are really needed, and if 10^10 ops take a second, 10^6 would take 1/10,000th of a second.
    – Will Ness
    Apr 12, 2014 at 15:51
  • 1
    Because I didn't realize (10^12+39) was a prime number which I do now. I get exactly what you're saying.
    – the_prole
    Apr 14, 2014 at 12:58
  • 1
    OK, so you can change your code so that there won't be such a big slowdown for the primes: if n = ab and a <= b, then aa <= ba = n, i.e. aa <= n. And if we've reached a+1, then n is surely a prime. (ping me if you edit your answer to incorporate this).
    – Will Ness
    Apr 14, 2014 at 13:30
  • 1
    what happens when long n = 2*1000000000039L? Does it work as fast as it should? (also, can you simplify your code by using a return; statement?). (if you want me to stop nudging you, just say so ;))
    – Will Ness
    Apr 18, 2014 at 8:37
4

The simplest solution is a pair of mutually recursive functions.

The first function generates all the prime numbers:

  1. Start with a list of all natural numbers greater than 1.
  2. Remove all numbers that are not prime. That is, numbers that have no prime factors (other than themselves). See below.

The second function returns the prime factors of a given number n in increasing order.

  1. Take a list of all the primes (see above).
  2. Remove all the numbers that are not factors of n.

The largest prime factor of n is the last number given by the second function.

This algorithm requires a lazy list or a language (or data structure) with call-by-need semantics.

For clarification, here is one (inefficient) implementation of the above in Haskell:

import Control.Monad

-- All the primes
primes = 2 : filter (ap (<=) (head . primeFactors)) [3,5..]

-- Gives the prime factors of its argument
primeFactors = factor primes
  where factor [] n = []
        factor xs@(p:ps) n =
          if p*p > n then [n]
          else let (d,r) = divMod n p in
            if r == 0 then p : factor xs d
            else factor ps n

-- Gives the largest prime factor of its argument
largestFactor = last . primeFactors

Making this faster is just a matter of being more clever about detecting which numbers are prime and/or factors of n, but the algorithm stays the same.

0
3
n = abs(number);
result = 1;
if (n mod 2 == 0) {
 result = 2;
 while (n mod 2 = 0) n /= 2;
}
for(i=3; i<sqrt(n); i+=2) {
 if (n mod i == 0) {
   result = i;
   while (n mod i = 0)  n /= i;
 }
}
return max(n,result)

There are some modulo tests that are superflous, as n can never be divided by 6 if all factors 2 and 3 have been removed. You could only allow primes for i, which is shown in several other answers here.

You could actually intertwine the sieve of Eratosthenes here:

  • First create the list of integers up to sqrt(n).
  • In the for loop mark all multiples of i up to the new sqrt(n) as not prime, and use a while loop instead.
  • set i to the next prime number in the list.

Also see this question.

0
2

I'm aware this is not a fast solution. Posting as hopefully easier to understand slow solution.

 public static long largestPrimeFactor(long n) {

        // largest composite factor must be smaller than sqrt
        long sqrt = (long)Math.ceil(Math.sqrt((double)n));

        long largest = -1;

        for(long i = 2; i <= sqrt; i++) {
            if(n % i == 0) {
                long test = largestPrimeFactor(n/i);
                if(test > largest) {
                    largest = test;
                }
            }
        }

        if(largest != -1) {
            return largest;
        }

        // number is prime
        return n;
    } 
0
2

Part 1. Pollard-Rho + Trial Division.

Inspired by your question I decided to implement my own version of factorization (and finding largest prime factor) in Python.

Probably the simplest to implement, yet quite efficient, factoring algorithm that I know is Pollard's Rho algorithm. It has a running time of O(N^(1/4)) at most which is much more faster than time of O(N^(1/2)) for trial division algorithm. Both algos have these running times only in case of composite (non-prime) number, that's why primality test should be used to filter out prime (non-factorable) numbers.

I used following algorithms in my code: Fermat Primality Test ..., Pollard's Rho Algorithm ..., Trial Division Algorithm. Fermat primality test is used before running Pollard's Rho in order to filter out prime numbers. Trial Division is used as a fallback because Pollard's Rho in very rare cases may fail to find a factor, especially for some small numbers.

Obviously after fully factorizing a number into sorted list of prime factors the largest prime factor will be the last element in this list. In general case (for any random number) I don't know of any other ways to find out largest prime factor besides fully factorizing a number.

As an example in my code I'm factoring first 190 fractional digits of Pi, code factorizes this number within 1 second, and shows largest prime factor which is 165 digits (545 bits) in size!

Try it online!

def is_fermat_probable_prime(n, *, trials = 32):
    # https://en.wikipedia.org/wiki/Fermat_primality_test
    import random
    if n <= 16:
        return n in (2, 3, 5, 7, 11, 13)
    for i in range(trials):
        if pow(random.randint(2, n - 2), n - 1, n) != 1:
            return False
    return True

def pollard_rho_factor(N, *, trials = 16):
    # https://en.wikipedia.org/wiki/Pollard%27s_rho_algorithm
    import random, math
    for j in range(trials):
        i, stage, y, x = 0, 2, 1, random.randint(1, N - 2)
        while True:
            r = math.gcd(N, x - y)
            if r != 1:
                break
            if i == stage:
                y = x
                stage <<= 1
            x = (x * x + 1) % N
            i += 1
        if r != N:
            return [r, N // r]
    return [N] # Pollard-Rho failed

def trial_division_factor(n, *, limit = None):
    # https://en.wikipedia.org/wiki/Trial_division
    fs = []
    while n & 1 == 0:
        fs.append(2)
        n >>= 1
    d = 3
    while d * d <= n and limit is None or d <= limit:
        q, r = divmod(n, d)
        if r == 0:
            fs.append(d)
            n = q
        else:
            d += 2
    if n > 1:
        fs.append(n)
    return fs

def factor(n):
    if n <= 1:
        return []
    if is_fermat_probable_prime(n):
        return [n]
    fs = trial_division_factor(n, limit = 1 << 12)
    if len(fs) >= 2:
        return sorted(fs[:-1] + factor(fs[-1]))
    fs = pollard_rho_factor(n)
    if len(fs) >= 2:
        return sorted([e1 for e0 in fs for e1 in factor(e0)])
    return trial_division_factor(n)

def demo():
    import time, math
    # http://www.math.com/tables/constants/pi.htm
    # pi = 3.
    #     1415926535 8979323846 2643383279 5028841971 6939937510 5820974944 5923078164 0628620899 8628034825 3421170679
    #     8214808651 3282306647 0938446095 5058223172 5359408128 4811174502 8410270193 8521105559 6446229489 5493038196
    #     4428810975 6659334461 2847564823 3786783165 2712019091 4564856692 3460348610 4543266482 1339360726 0249141273
    #     7245870066 0631558817 4881520920 9628292540 9171536436 7892590360 0113305305 4882046652 1384146951 9415116094
    #     3305727036 5759591953 0921861173 8193261179 3105118548 0744623799 6274956735 1885752724 8912279381 8301194912
    #     9833673362 4406566430 8602139494 6395224737 1907021798 6094370277 0539217176 2931767523 8467481846 7669405132
    #
    # n = first 190 fractional digits of Pi
    n =   1415926535_8979323846_2643383279_5028841971_6939937510_5820974944_5923078164_0628620899_8628034825_3421170679_8214808651_3282306647_0938446095_5058223172_5359408128_4811174502_8410270193_8521105559_6446229489
    print('Number:', n)
    tb = time.time()
    fs = factor(n)
    print('All Prime Factors:', fs)
    print('Largest Prime Factor:', f'({math.log2(fs[-1]):.02f} bits, {len(str(fs[-1]))} digits)', fs[-1])
    print('Time Elapsed:', round(time.time() - tb, 3), 'sec')

if __name__ == '__main__':
    demo()

Output:

Number: 1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489
All Prime Factors: [3, 71, 1063541, 153422959, 332958319, 122356390229851897378935483485536580757336676443481705501726535578690975860555141829117483263572548187951860901335596150415443615382488933330968669408906073630300473]
Largest Prime Factor: (545.09 bits, 165 digits) 122356390229851897378935483485536580757336676443481705501726535578690975860555141829117483263572548187951860901335596150415443615382488933330968669408906073630300473
Time Elapsed: 0.593 sec

Part 2. Elliptic Curve Method.

Some time later I decided to improve my post by implementing from scratch more advanced elliptic curve factorization method that is called ECM, see Wikipedia Article Lenstra Elliptic Curve Factorization.

This method is considerably faster than Pollard Rho described in Part 1 of my answer post. Time complexity of elliptic ECM method is O(exp[(Sqrt(2) + o(1)) Sqrt(ln p ln ln p)]), where p signifies smallest prime factor of a number. While time complexity of Pollard Rho is O(Sqrt(p)). So ECM is much much faster for big enough smallest P factor.

Steps of ECM method:

  1. Check if number is smaller than 2^16, then factor it through Trial Division method. Return result.

  2. Check if number is probably prime with high condifence, for that I use Fermat Test with 32 trials. To overcome Carmichael numbers you may use Miller Rabin Test instead. If number is primes return it as only factor.

  3. Generate curve parameters A, X, Y randomly and derive B from curve equation Y^2 = X^3 + AX + B (mod N). Check if curve is OK, value 4 * A ** 3 - 27 * B ** 2 should be non-zero.

  4. Generate small primes through Sieve of Eratosthenes, primes below our Bound. Each prime raise to some small power, this raised prime would be called K. I do raising into power while it is smaller than some Bound2, which is Sqrt(Bound) in my case.

  5. Compute elliptic point multiplication P = k * P, where K taken from previous step and P is (X, Y). Compute according to Wiki.

  6. Point multiplication uses Modular Inverse, which computes GCD(SomeValue, N) according to Wiki. If this GCD is not 1, then it gives non-1 factor of N, hence in this case I through an Exception and return this factor from ECM factorization algorithm.

  7. If all primes till Bound were multiplied and gave no factor then re-run ECM factorization algorithm (1.-6. above) again with another random curve and bigger Bound. In my code I take new bound by adding 256 to old bound.

Following sub-algorithms were used in my ECM code, all mentioned sub-algorithms have links to corresponding Wikipedia articles: Trial Division Factorization, Fermat Probability Test, Sieve of Eratosthenes (prime numbers generator), Euclidean Algorithm (computing Greatest Common Divisor, GCD), Extended Euclidean Algorithm (GCD with Bezu coefficients), Modular Multiplicative Inverse, Right-to-Left Binary Exponentation (for elliptic point multiplication), Elliptic Curve Arithmetics (point addition and multiplication), Lenstra Elliptic Curve Factorization.

Unlike Part 1 of my answer where I factor digits of Pi number, here in Part 2 I create a special number composed of n = Prime(24) * Prime(35) * Prime(37) * ... etc, which means number as a product of Random prime numbers of 24 bits and 35 bits and 37 and etc... This custom number is more visually impressive to show algorithm capabilities.

As in Part-1 of my answer I also use several methods to compare there speed and also to factor-out smaller factors with simpler methods. So in code below I use Trial Division + Pollard Rho + Elliptic Curve Method.

After code below see console output to figure out what nicde stuff my code outputs.

Try it online!

def is_fermat_probable_prime(n, *, trials = 32):
    # https://en.wikipedia.org/wiki/Fermat_primality_test
    import random
    if n <= 16:
        return n in (2, 3, 5, 7, 11, 13)
    for i in range(trials):
        if pow(random.randint(2, n - 2), n - 1, n) != 1:
            return False
    return True

def pollard_rho_factor(N, *, trials = 8, limit = None):
    # https://en.wikipedia.org/wiki/Pollard%27s_rho_algorithm
    import random, math
    if N <= 1:
        return []
    if is_fermat_probable_prime(N):
        return [N]
    for j in range(trials):
        i, stage, y, x = 0, 2, 1, random.randint(1, N - 2)
        while True:
            r = math.gcd(N, x - y)
            if r != 1:
                break
            if i == stage:
                y = x
                stage <<= 1
            x = (x * x + 1) % N
            i += 1
            if limit is not None and i >= limit:
                return [N] # Pollard-Rho failed
        if r != N:
            return sorted(pollard_rho_factor(r, trials = trials, limit = limit) +
                pollard_rho_factor(N // r, trials = trials, limit = limit))
    return [N] # Pollard-Rho failed

def trial_division_factor(n, *, limit = None):
    # https://en.wikipedia.org/wiki/Trial_division
    fs = []
    while n & 1 == 0:
        fs.append(2)
        n >>= 1
    d = 3
    while d * d <= n and limit is None or d <= limit:
        q, r = divmod(n, d)
        if r == 0:
            fs.append(d)
            n = q
        else:
            d += 2
    if n > 1:
        fs.append(n)
    return fs
    
def ecm_factor(N0, *, verbose = False):
    # https://en.wikipedia.org/wiki/Lenstra_elliptic-curve_factorization
    import math, random, time; gmpy2 = None
    #import gmpy2
    def GenPrimes_SieveOfEratosthenes(end):
        # https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes
        composite = [False] * (end // 2) # Allocate for odd numbers only
        for p in range(3, int(end ** 0.5 + 3), 2):
            if composite[p >> 1]:
                continue
            for i in range(p * p, end, p * 2):
                composite[i >> 1] = True
        yield 2
        for p in range(3, end, 2):
            if not composite[p >> 1]:
                yield p
    def GCD(a, b):
        # https://en.wikipedia.org/wiki/Euclidean_algorithm
        # return math.gcd(a, b)
        while b != 0:
            a, b = b, a % b
        return a
    def EGCD(a, b):
        # https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm
        if gmpy2 is None:
            ro, r, so, s = a, b, 1, 0
            while r != 0:
                ro, (q, r) = r, divmod(ro, r)
                so, s = s, so - q * s
            return ro, so, (ro - so * a) // b
        else:
            return tuple(map(int, gmpy2.gcdext(a, b)))
    def ModularInverse(a, n):
        # https://en.wikipedia.org/wiki/Modular_multiplicative_inverse
        # return pow(a, -1, n)
        g, s, t = EGCD(a, n)
        if g != 1:
            raise ValueError(a)
        return s % n
    def EllipticCurveAdd(N, A, B, X0, Y0, X1, Y1):
        # https://en.wikipedia.org/wiki/Elliptic_curve_point_multiplication
        if X0 == X1 and Y0 == Y1:
            # Double
            l = ((3 * X0 ** 2 + A) * ModularInverse(2 * Y0, N)) % N
            x = (l ** 2 - 2 * X0) % N
            y = (l * (X0 - x) - Y0) % N
        else:
            # Add
            l = ((Y1 - Y0) * ModularInverse(X1 - X0, N)) % N
            x = (l ** 2 - X0 - X1) % N
            y = (l * (X0 - x) - Y0) % N
        return x, y
    def EllipticCurveMul(N, A, B, X, Y, k):
        # https://en.wikipedia.org/wiki/Modular_exponentiation#Right-to-left_binary_method
        assert k >= 2, k
        k -= 1
        BX, BY = X, Y
        while k != 0:
            if k & 1:
                X, Y = EllipticCurveAdd(N, A, B, X, Y, BX, BY)
            BX, BY = EllipticCurveAdd(N, A, B, BX, BY, BX, BY)
            k >>= 1
        return X, Y
    
    bound_start = 1 << 9
    
    def Main(N, *, bound = bound_start, icurve = 0):
        def NextFactorECM(x):
            return Main(x, bound = bound + bound_start, icurve = icurve + 1)
        def PrimePow(p, *, bound2 = int(bound ** 0.5 + 1.01)):
            mp = p
            while True:
                mp *= p
                if mp >= bound2:
                    return mp // p
        
        if N <= 1:
            return []
            
        if N < (1 << 16):
            fs = trial_division_factor(N)
            if verbose and len(fs) >= 2:
                print('Factors from TrialDiv:', fs, flush = True)
            return fs
            
        if is_fermat_probable_prime(N):
            return [N]
        
        if verbose:
            print(f'Curve {icurve:>4},  bound 2^{math.log2(bound):>7.3f}', flush = True)
        
        while True:
            X, Y, A = [random.randrange(N) for i in range(3)]
            B = (Y ** 2 - X ** 3 - A * X) % N
            if 4 * A ** 3 - 27 * B ** 2 == 0:
                continue
            break
        
        for p in GenPrimes_SieveOfEratosthenes(bound):
            k = PrimePow(p)
            try:
                X, Y = EllipticCurveMul(N, A, B, X, Y, k)
            except ValueError as ex:
                g = GCD(ex.args[0], N)
                assert g > 1, g
                if g != N:
                    if verbose:
                        print(f'Factor from ECM: {g} ({math.log2(g):.1f} bits)', flush = True)
                    return sorted(NextFactorECM(g) + NextFactorECM(N // g))
                else:
                    return NextFactorECM(N)
        
        return NextFactorECM(N)
    
    return Main(N0)

def factor(n):
    if n <= 1:
        return []
    if is_fermat_probable_prime(n):
        return [n]
    fs1 = trial_division_factor(n, limit = 1 << 12)
    fs1, n2 = fs1[:-1], fs1[-1]
    print(len(fs1), 'factors from TrialDivision')
    fs2 = pollard_rho_factor(n2, limit = 1 << 17)
    fs2, n3 = fs2[:-1], fs2[-1]
    print(len(fs2), 'factors from PollardRho')
    fs3 = ecm_factor(n3, verbose = True)
    print(len(fs3), 'factors from ECM')
    return sorted(fs1 + fs2 + fs3)

def demo():
    import time, math, random
    def Prime(bits):
        while True:
            x = random.randrange(1 << (bits - 1), 1 << bits)
            if is_fermat_probable_prime(x):
                return x
    
    # http://www.math.com/tables/constants/pi.htm
    # pi = 3.
    #     1415926535 8979323846 2643383279 5028841971 6939937510 5820974944 5923078164 0628620899 8628034825 3421170679
    #     8214808651 3282306647 0938446095 5058223172 5359408128 4811174502 8410270193 8521105559 6446229489 5493038196
    #     4428810975 6659334461 2847564823 3786783165 2712019091 4564856692 3460348610 4543266482 1339360726 0249141273
    #     7245870066 0631558817 4881520920 9628292540 9171536436 7892590360 0113305305 4882046652 1384146951 9415116094
    #     3305727036 5759591953 0921861173 8193261179 3105118548 0744623799 6274956735 1885752724 8912279381 8301194912
    #     9833673362 4406566430 8602139494 6395224737 1907021798 6094370277 0539217176 2931767523 8467481846 7669405132
    #
    # n =   1415926535_8979323846_2643383279_5028841971_6939937510_5820974944_5923078164_0628620899_8628034825_3421170679_8214808651_3282306647_0938446095_5058223172_5359408128_4811174502_8410270193_8521105559_6446229489_5493038196_4428810975_6659334461_2847564823_3786783165_2712019091_4564856692_3460348610_4543266482_1339360726_0249141273_7245870066_0631558817_4881520920_9628292540_9171536436_7892590360_0113305305_4882046652_1384146951_9415116094_3305727036_5759591953_0921861173_8193261179_3105118548_0744623799_6274956735_1885752724_8912279381_8301194912_9833673362_4406566430_8602139494_6395224737_1907021798_6094370277_0539217176_2931767523_8467481846_7669405132
    #       141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587006606315588174881520920962829254091715364367892590360011330530548820466521384146951941511609433057270365759591953092186117381932611793105118548074462379962749567351885752724891227938183011949129833673362440656643086021394946395224737190702179860943702770539217176293176752384674818467669405132
    
    n = Prime(9) * Prime(10) * Prime(11) * Prime(20) * Prime(24) * Prime(35) * Prime(37) * Prime(38) * Prime(120)
    print('Number:', n)
    tb = time.time()
    fs = factor(n)
    print('All Prime Factors:', fs)
    print('Largest Prime Factor:', f'({math.log2(fs[-1]):.02f} bits, {len(str(fs[-1]))} digits)', fs[-1])
    if len(fs) >= 2:
        print('2nd-largest Prime Factor:', f'({math.log2(fs[-2]):.02f} bits, {len(str(fs[-2]))} digits)', fs[-2])
    print('Time Elapsed:', round(time.time() - tb, 3), 'sec')

if __name__ == '__main__':
    demo()

Console output:

Number: 3020823358956369790763854998578637168366763837218991014777892420353187988302225517459334041

3 factors from TrialDivision
2 factors from PollardRho

Curve    0,  bound 2^  9.000
Curve    1,  bound 2^ 10.000
Curve    2,  bound 2^ 10.585
Factor from ECM: 20028139561 (34.2 bits)
Curve    3,  bound 2^ 11.000
Curve    4,  bound 2^ 11.322
Curve    5,  bound 2^ 11.585
Curve    6,  bound 2^ 11.807
Curve    7,  bound 2^ 12.000
Curve    8,  bound 2^ 12.170
Factor from ECM: 96583780901 (36.5 bits)
Curve    9,  bound 2^ 12.322
Curve   10,  bound 2^ 12.459
Curve   11,  bound 2^ 12.585
Curve   12,  bound 2^ 12.700
Curve   13,  bound 2^ 12.807
Curve   14,  bound 2^ 12.907
Curve   15,  bound 2^ 13.000
Curve   16,  bound 2^ 13.087
Curve   17,  bound 2^ 13.170
Factor from ECM: 239171423261 (37.8 bits)
4 factors from ECM

All Prime Factors: [397, 1021, 1459, 754333, 16156687, 20028139561,
    96583780901, 239171423261, 905908369146483365552973334921981697]
Largest Prime Factor: (119.45 bits, 36 digits) 905908369146483365552973334921981697
2nd-largest Prime Factor: (37.80 bits, 12 digits) 239171423261

Time Elapsed: 17.156 sec
1

Python Iterative approach by removing all prime factors from the number

def primef(n):
    if n <= 3:
        return n
    if n % 2 == 0:
        return primef(n/2)
    elif n % 3 ==0:
        return primef(n/3)
    else:
        for i in range(5, int((n)**0.5) + 1, 6):
            #print i
            if n % i == 0:
                return primef(n/i)
            if n % (i + 2) == 0:
                return primef(n/(i+2))
    return n
1

I am using algorithm which continues dividing the number by it's current Prime Factor.

My Solution in python 3 :

def PrimeFactor(n):
    m = n
    while n%2==0:
        n = n//2
    if n == 1:         # check if only 2 is largest Prime Factor 
        return 2
    i = 3
    sqrt = int(m**(0.5))  # loop till square root of number
    last = 0              # to store last prime Factor i.e. Largest Prime Factor
    while i <= sqrt :
        while n%i == 0:
            n = n//i       # reduce the number by dividing it by it's Prime Factor
            last = i
        i+=2
    if n> last:            # the remaining number(n) is also Factor of number 
        return n
    else:
        return last
print(PrimeFactor(int(input()))) 

Input : 10 Output : 5

Input : 600851475143 Output : 6857

0
0

Here is my attempt in c#. The last print out is the largest prime factor of the number. I checked and it works.

namespace Problem_Prime
{
  class Program
  {
    static void Main(string[] args)
    {
      /*
       The prime factors of 13195 are 5, 7, 13 and 29.

      What is the largest prime factor of the number 600851475143 ?
       */
      long x = 600851475143;
      long y = 2;
      while (y < x)
      {
        if (x % y == 0)
        {
          // y is a factor of x, but is it prime
          if (IsPrime(y))
          {
            Console.WriteLine(y);
          }
          x /= y;
        }

        y++;

      }
      Console.WriteLine(y);
      Console.ReadLine();
    }
    static bool IsPrime(long number)
    {
      //check for evenness
      if (number % 2 == 0)
      {
        if (number == 2)
        {
          return true;
        }
        return false;
      }
      //don't need to check past the square root
      long max = (long)Math.Sqrt(number);
      for (int i = 3; i <= max; i += 2)
      {
        if ((number % i) == 0)
        {
          return false;
        }
      }
      return true;
    }

  }
}
0
#python implementation
import math
n = 600851475143
i = 2
factors=set([])
while i<math.sqrt(n):
   while n%i==0:
       n=n/i
       factors.add(i)
   i+=1
factors.add(n)
largest=max(factors)
print factors
print largest
1
  • 1
    is 25 the largest prime factor of 25?
    – Will Ness
    May 31, 2014 at 23:18
0

Calculates the largest prime factor of a number using recursion in C++. The working of the code is explained below:

int getLargestPrime(int number) {
    int factor = number; // assumes that the largest prime factor is the number itself
    for (int i = 2; (i*i) <= number; i++) { // iterates to the square root of the number till it finds the first(smallest) factor
        if (number % i == 0) { // checks if the current number(i) is a factor
            factor = max(i, number / i); // stores the larger number among the factors
            break; // breaks the loop on when a factor is found
        }
    }
    if (factor == number) // base case of recursion
        return number;
    return getLargestPrime(factor); // recursively calls itself
}
0

Here is my approach to quickly calculate the largest prime factor. It is based on fact that modified x does not contain non-prime factors. To achieve that, we divide x as soon as a factor is found. Then, the only thing left is to return the largest factor. It would be already prime.

The code (Haskell):

f max' x i | i > x = max'
           | x `rem` i == 0 = f i (x `div` i) i  -- Divide x by its factor
           | otherwise = f max' x (i + 1)  -- Check for the next possible factor

g x = f 2 x 2
1
  • but won't this try to divide with all even numbers too? Feb 17, 2016 at 2:30
0

The following C++ algorithm is not the best one, but it works for numbers under a billion and its pretty fast

#include <iostream>
using namespace std;

// ------ is_prime ------
// Determines if the integer accepted is prime or not
bool is_prime(int n){
    int i,count=0;
    if(n==1 || n==2)
      return true;
    if(n%2==0)
      return false;
    for(i=1;i<=n;i++){
    if(n%i==0)
        count++;
    }
    if(count==2)
      return true;
    else
      return false;
 }
 // ------ nextPrime -------
 // Finds and returns the next prime number
 int nextPrime(int prime){
     bool a = false;
     while (a == false){
         prime++;
         if (is_prime(prime))
            a = true;
     }
  return prime;
 }
 // ----- M A I N ------
 int main(){

      int value = 13195;
      int prime = 2;
      bool done = false;

      while (done == false){
          if (value%prime == 0){
             value = value/prime;
             if (is_prime(value)){
                 done = true;
             }
          } else {
             prime = nextPrime(prime);
          }
      }
        cout << "Largest prime factor: " << value << endl;
 }
0

Found this solution on the web by "James Wang"

public static int getLargestPrime( int number) {

    if (number <= 1) return -1;

    for (int i = number - 1; i > 1; i--) {
        if (number % i == 0) {
            number = i;
        }
    }
    return number;
}
0

Prime factor using sieve :

#include <bits/stdc++.h>
using namespace std;
#define N 10001  
typedef long long ll;
bool visit[N];
vector<int> prime;

void sieve()
{
            memset( visit , 0 , sizeof(visit));
            for( int i=2;i<N;i++ )
            {
                if( visit[i] == 0)
                {
                    prime.push_back(i);
                    for( int j=i*2; j<N; j=j+i )
                    {
                        visit[j] = 1;
                    }
                }
            }   
}
void sol(long long n, vector<int>&prime)
{
            ll ans = n;
            for(int i=0; i<prime.size() || prime[i]>n; i++)
            {
                while(n%prime[i]==0)
                {
                    n=n/prime[i];
                    ans = prime[i];
                }
            }
            ans = max(ans, n);
            cout<<ans<<endl;
}
int main() 
{
           ll tc, n;
           sieve();

           cin>>n;
           sol(n, prime);

           return 0;
}
0

Guess, there is no immediate way but performing a factorization, as examples above have done, i.e.

in a iteration you identify a "small" factor f of a number N, then continue with the reduced problem "find largest prime factor of N':=N/f with factor candidates >=f ".

From certain size of f the expected search time is less, if you do a primality test on reduced N', which in case confirms, that your N' is already the largest prime factor of initial N.

0

Here is my attempt in Clojure. Only walking the odds for prime? and the primes for prime factors ie. sieve. Using lazy sequences help producing the values just before they are needed.

(defn prime? 
  ([n]
    (let [oddNums (iterate #(+ % 2) 3)]
    (prime? n (cons 2 oddNums))))
  ([n [i & is]]
    (let [q (quot n i)
          r (mod n i)]
    (cond (< n 2)       false
          (zero? r)     false
          (> (* i i) n) true
          :else         (recur n is)))))

(def primes 
  (let [oddNums (iterate #(+ % 2) 3)]
  (lazy-seq (cons 2 (filter prime? oddNums)))))

;; Sieve of Eratosthenes
(defn sieve
  ([n] 
    (sieve primes n))
  ([[i & is :as ps] n]
    (let [q (quot n i)
          r (mod n i)]
    (cond (< n 2)       nil
          (zero? r)     (lazy-seq (cons i (sieve ps q)))
          (> (* i i) n) (when (> n 1) (lazy-seq [n]))
          :else         (recur is n)))))

(defn max-prime-factor [n]
  (last (sieve n)))
0

Recursion in C

Algorithm could be

  1. Check if n is a factor or t
  2. Check if n is prime. If so, remember n
  3. Increment n
  4. Repeat until n > sqrt(t)

Here's an example of a (tail)recursive solution to the problem in C:

#include <stdio.h>
#include <stdbool.h>

bool is_factor(long int t, long int n){
    return ( t%n == 0);
}

bool is_prime(long int n0, long int n1, bool acc){
    if ( n1 * n1 > n0 || acc < 1 )
        return acc;
    else
        return is_prime(n0, n1+2, acc && (n0%n1 != 0));
}

int gpf(long int t, long int n, long int acc){
    if (n * n > t)
        return acc;
    if (is_factor(t, n)){
        if (is_prime(n, 3, true))
            return gpf(t, n+2, n);
        else
            return gpf(t, n+2, acc);
    }
    else
        return gpf(t, n+2, acc);
}

int main(int argc, char ** argv){
    printf("%d\n", gpf(600851475143, 3, 0));
    return 0;
}

The solution is composed of three functions. One to test if the candidate is a factor, another to test if that factor is prime, and finally one to compose those two together.

Some key ideas here are:

1- Stopping the recursion at sqrt(600851475143)

2- Only test odd numbers for factorness

3- Only testing candidate factors for primeness with odd numbers

-1

It seems to me that step #2 of the algorithm given isn't going to be all that efficient an approach. You have no reasonable expectation that it is prime.

Also, the previous answer suggesting the Sieve of Eratosthenes is utterly wrong. I just wrote two programs to factor 123456789. One was based on the Sieve, one was based on the following:

1)  Test = 2 
2)  Current = Number to test 
3)  If Current Mod Test = 0 then  
3a)     Current = Current Div Test 
3b)     Largest = Test
3c)     Goto 3. 
4)  Inc(Test) 
5)  If Current < Test goto 4
6)  Return Largest

This version was 90x faster than the Sieve.

The thing is, on modern processors the type of operation matters far less than the number of operations, not to mention that the algorithm above can run in cache, the Sieve can't. The Sieve uses a lot of operations striking out all the composite numbers.

Note, also, that my dividing out factors as they are identified reduces the space that must be tested.

3
  • that's what i said, but got voted down :( I guess the problem is that if the number has a really large prime factor (such as itself), then this method must loop all the way up to that number. In a lot of cases though, this method is quite efficient.
    – nickf
    Oct 28, 2008 at 5:57
  • Reading back through yours it is the same but the first part of yours is confusing. Oct 29, 2008 at 1:56
  • Try that on this number 143816789988504044536402352738195137863656439, let me know how efficient this is...
    – Kladskull
    May 8, 2009 at 16:54
-1

Compute a list storing prime numbers first, e.g. 2 3 5 7 11 13 ...

Every time you prime factorize a number, use implementation by Triptych but iterating this list of prime numbers rather than natural integers.

0
-1

With Java:

For int values:

public static int[] primeFactors(int value) {
    int[] a = new int[31];
    int i = 0, j;
    int num = value;
    while (num % 2 == 0) {
        a[i++] = 2;
        num /= 2;
    }
    j = 3;
    while (j <= Math.sqrt(num) + 1) {
        if (num % j == 0) {
            a[i++] = j;
            num /= j;
        } else {
            j += 2;
        }
    }
    if (num > 1) {
        a[i++] = num;
    }
    int[] b = Arrays.copyOf(a, i);
    return b;
}

For long values:

static long[] getFactors(long value) {
    long[] a = new long[63];
    int i = 0;
    long num = value;
    while (num % 2 == 0) {
        a[i++] = 2;
        num /= 2;
    }
    long j = 3;
    while (j <= Math.sqrt(num) + 1) {
        if (num % j == 0) {
            a[i++] = j;
            num /= j;
        } else {
            j += 2;
        }
    }
    if (num > 1) {
        a[i++] = num;
    }
    long[] b = Arrays.copyOf(a, i);
    return b;
}
0
-3

This is probably not always faster but more optimistic about that you find a big prime divisor:

  1. N is your number
  2. If it is prime then return(N)
  3. Calculate primes up until Sqrt(N)
  4. Go through the primes in descending order (largest first)
    • If N is divisible by Prime then Return(Prime)

Edit: In step 3 you can use the Sieve of Eratosthenes or Sieve of Atkins or whatever you like, but by itself the sieve won't find you the biggest prime factor. (Thats why I wouldn't choose SQLMenace's post as an official answer...)

1
  • 1
    Don't you need to do the trial factoring to determine if it is a prime number (step 2)? Also, consider finding the largest prime factor of 15. The primes up to sqrt(15) are 2 and 3; but the largest prime factor is 5, is it not? Similarly with 20. Jul 16, 2016 at 6:52
-3

Here is the same function@Triptych provided as a generator, which has also been simplified slightly.

def primes(n):
    d = 2
    while (n > 1):
        while (n%d==0):
            yield d
            n /= d
        d += 1

the max prime can then be found using:

n= 373764623
max(primes(n))

and a list of factors found using:

list(primes(n))
0
-5

I think it would be good to store somewhere all possible primes smaller then n and just iterate through them to find the biggest divisior. You can get primes from prime-numbers.org.

Of course I assume that your number isn't too big :)

-6
#include<stdio.h>
#include<conio.h>
#include<math.h>
#include <time.h>

factor(long int n)
{
long int i,j;
while(n>=4)
 {
if(n%2==0) {  n=n/2;   i=2;   }

 else
 { i=3;
j=0;
  while(j==0)
  {
   if(n%i==0)
   {j=1;
   n=n/i;
   }
   i=i+2;
  }
 i-=2;
 }
 }
return i;
 }

 void main()
 { 
  clock_t start = clock();
  long int n,sp;
  clrscr();
  printf("enter value of n");
  scanf("%ld",&n);
  sp=factor(n);
  printf("largest prime factor is %ld",sp);

  printf("Time elapsed: %f\n", ((double)clock() - start) / CLOCKS_PER_SEC);
  getch();
 }
0

Not the answer you're looking for? Browse other questions tagged or ask your own question.