Algorithm to find Largest prime factor of a number

Question

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.

Culliton

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. – Bucharest

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 – Bithynia

We can Divide with small primes, and the one which is finally left, is the Largest Prime Factor (I guess) – Bilabial

Answer 1

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.

Mockery

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. – Vestavestal

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 @Culliton suggested in his question. The quadratic solution is what is powering my client's tool now at math.tools/calculator/prime-factors . – Wieland

If there is an efficient way of solving this problem, doesn't that mean that web encryption is not not secure? – Secure

Answer 2

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

Handstand

I do believe this does not work, nor does it retrieve prime-factors, but factors eventually. Where in this code are prime numbers used? d is only natural numbers upcounted, nothing prime there. – Jerrelljerri

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. – Handstand

I'm sorry, did an error in converting the code to C#, put one line more into the second while loop. Undone the downvote. – Jerrelljerri

"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.) – Porthole

@Triptych Sheldon is correct, this is O(n) for prime values of n. – Tibia

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. – Unspeakable

Is there a name for this? – Monet

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... :) – Matchmaker

@Monet I think you would call this approach "Trial Division" or "Brute Force". – Firebrat

This algorithm is horribly inefficient. And the claim that it's O(sqrt(n)) is misleading: for problems such as integer factorization it's customary to measure complexity in the number of input digits, not the size of the input. Given than the input n has m=log(n) digits the complexity of this algorithm is O(exp(m/2)). – Chromomere

@avi The trick is to realize that because this loop starts at the first prime (d=2) and removes all factors, d will always be a prime if n % d == 0. Say n = 120. 120's prime factors are 2,2,2,3,5. So when d = 2, the (while n % d == 0) strips all 2's from the number, and you're left with some other prime factors. The key idea is those prime factors will be larger than 2. The other cool part is the loop terminates at sqrt(n), because there cannot be any prime factors larger than sqrt(n). – Berkeley

@avi or anyone else who is still confused: another way to think of it is first: The smallest divisor will always be prime. Take 12 for example. 2 is the smallest divisor. Then note that the list of all prime numbers times each other will equal the total. So e.g. we know 2 is a prime. So we can divide 12 by 2 = 6. We know the rest of the primes multiplied by each other must equal 6. Lucas explained the rest well. – Nadinenadir

What is the trick that causes d = d + 1 to be ignored for numbers that are not a prime? – Guimond

The 'trick' is the the line n /= d. It removes the prime factor d from n for all subsequent calculations in the loop. – Hefty

I think this method will only generate a set of prime factors only, the result set can't have any composite since we start with the first prime 2 in increasing order always, which is neat – Joung

Inefficient, but it works well enough for my purposes. But it can return a real number (199.0) instead of an integer (199). Simple fix: replace n/=d with n//=d. – Dodds

@Redsandro, VinceOSullivan: there is no 'trick'. d = d + 1 causes us to brute-force all candidate numbers d <= sqrt(n). (If d is composite, then n cannot still be divisible by d since we've already removed all factors < d. But this algorithm still attempts the trial division while n % d == 0 anyway, even if d is composite. That's why this is very inefficient. The higher d gets, the more likely it's composite and hence doesn't need to be tested.) – Kraus

Answer 3

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.

Cristiano

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) – Benzoic

I guess it should be while (multOfSix - 1 <= n) – Mortie

Answer 4

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

Heterogeneity

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. – Liftoff

Another else-if can significantly cut down on wasted compute: else if i > sqrt(number): i = number – Cassette

(Assigning i = number is a very indirect way to terminate the loop - I'd just return number.) – Tricho

@Tricho Yeah, but in this case I prefer as-is – Anorthite

Answer 5

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:

Moonlit

Answer 6

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.

Messene

Answer 7

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.
}

Germain

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. – Anastomosis

+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. – Handstand

Same as with Triptych - nothing with prime-numbers here, same code, other language. – Jerrelljerri

@Beowulf. Try running this code before voting down. It returns prime factors; you just don't understand the algorithm. – Handstand

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. – Roping

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. – Roam

Answer 8

    //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
    }

Squad

have you tried your code with 1,000,000,000,039? it should run in a blink of an eye too. Does it? – Roam

I don't know Wily. I will have to try. – Squad

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. – Roam

Sorry Will, I didn't understand if you were being sarcastic. The speed of my code increases linearly with the size of the number. I don't see why my program would run any slower if I used one trillion and thirty nine. I wrote the code for general purpose factoring, not for specific numbers. I just used 600 billion as an example. So while your statement is technically correct, your method won't work for every number. Or will it? – Squad

Because I didn't realize (10^12+39) was a prime number which I do now. I get exactly what you're saying. – Squad

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). – Roam

I am writing my thesis and don't very much time right now. As soon as I implement the Fermat/trial-divsion combination I will ping you... How DO I ping you? – Squad

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 ;)) – Roam

@WillNess So isn't the second edit in Triptych's post incorrect? Shouldn't he be using a counter as I did in my code? – Squad

Answer 9

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.

Treharne

Answer 10

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;
    } 

Crystlecs

Answer 11

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

Clementia

Answer 12

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

Nickeliferous

Answer 13

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

Airfield

Answer 14

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;
    }

  }
}

Brashear

Answer 15

#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

Rouse

is 25 the largest prime factor of 25? – Roam

Answer 16

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
}

Laurettalaurette

Answer 17

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

Doxology

but won't this try to divide with all even numbers too? – Toolmaker

Answer 18

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;
 }

Approver

Answer 19

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;
}

Prorogue

Answer 20

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;
}

Aeromarine

Answer 21

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)))

Bedbug

Answer 22

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.

Stavros

Answer 23

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

Tarsier

Answer 24

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.

Torbernite

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. – Germain

Reading back through yours it is the same but the first part of yours is confusing. – Torbernite

Try that on this number 143816789988504044536402352738195137863656439, let me know how efficient this is... – Benge

Answer 25

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.

Aloise

Answer 26

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;
}

Gayle

Answer 27

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...)

Rose

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. – Adenocarcinoma

Answer 28

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))

Phosphatase

Answer 29

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 :)

Methodius

Answer 30

#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();
 }

Railroader