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.

Answer 1

Start at the large values and count down. We are looking for largest not smallest. We don't need to find all of them The largest prime factor is a bit less than sqrt(n)+1 So we can start there. I wanted to skip numbers that are not 6k±1. So we need to mod the sqrt by 6 and adjust it to the nearest 6k+1 number. Then we can test 6k+1 and 6k-1 all the way down to 5. (31,29,25,23,19,17,13,11,7,5) So the for loop subtracts 6 every time.

def is_prime(n):  
    if (n <= 3): return n > 1      
    if not (n%6 == 1 or n%6 == 5): return False      
    for p,q in ((i,i+2) for i in range(5,int(n**.5)+1,6)): 
        if (n%p == 0 or n%q == 0):  return False          
    return True  

def gpf(n):  # greatest Prime Factor
    a = int(n**.5)+1
    start = a
    if a%6 != 1: # 6k+1? 
        start = (a-5)+(a%6-((a%2)*2))+4 # make it 6k+1
    for p,q in ((i,i-2) for i in range(start,5,-6)): 
        if n % p == 0 and is_prime(p): return p # 6k+1
        if n % q == 0 and is_prime(q): return q # 6k-1
    if n%3 == 0: return 3  # last resort 3 and then 2
    if n%2 == 0: return 2
    return n # n is prime return it

If you want to use sympy libraries you can:

from sympy import prevprime

def find_largest_prime_factor(n):
   p = (n**.5)+1
   while p > 2:
      p = prevprime(p)
      if n%p == 0: return p
   return n # n is prime    

Answer 2

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

Answer 3

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

Answer 4

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

Answer 5

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.

Answer 6

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

Answer 7

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.

Answer 8

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.

Answer 9

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

Answer 10

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

Answer 11

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

Answer 12

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.

Answer 13

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

Answer 14

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:

Answer 15

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

Answer 16

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

Answer 17

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
}

Answer 18

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

Answer 19

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.

Answer 20

#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

Answer 21

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

Answer 22

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

Answer 23

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

  }
}

Answer 24

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

Answer 25

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

Answer 26

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

Answer 27

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.

Answer 28

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

Answer 29

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