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

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.

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.

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.

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

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.

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

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

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

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

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.

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

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.

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.

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

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