Creative Scratchpad

It gets more and more interesting with every new problem, let's have a look at problem No. 3:

The prime factors of 13195 are 5, 7, 13 and 29.

What is the largest prime factor of the number 600851475143 ?

Next on are prime numbers combined with a search for the factors of a number. The task is to find the largest prime factor of a number, which makes it suitable for some optimizations. The smallest prime factor of an even number is 2, and even though the given number is not even let's take it into account for the sake of reusability. Let's use the function from the previous task for this.

  fn isEven nr = bit.and 1L nr == 0

Checking for primarility is quite easy. We've already skipped number 2 with the previous step so the iteration can start at 3 and continue until we reach a square root of the number. To catch multples of two the function is checking if the number is even and returns true only if it is the number two itself (as we know, this won't happen here, yet it doesn't slow the function down to the point where we would consider leaving it out).

The reason for stopping at the square root is that if there would be any factor bigger than the square root, we would have already found a smaller factor before and returned false, as the bigger number would need a smaller number to be multiplied by to result in the original number. If there's no such smaller number we don't need to check for bigger ones and we can safely return true. As the input value is bigger than 2147483647 I'm using integer64, and as stated last time in that case nr^0.5 returns integer64 as well and doesn't need any rounding down. Since we can leave out all even numbers, the step for the loop is 2. If the number is divisible by any of the numbers in the loop without remainder, it isn't prime.

fn isPrime nr =
(
    if isEven nr do
        return nr == 2

    local rootNr = nr^0.5
    for i in 3L to rootNr by 2L
        where mod nr i == 0 do
            return false
    true
)

Now it's time for some little recursion. Let's take the number and search for the smallest prime factor we can find using the exact same method as before, trial division with the help of modulo function. Then we divide the number by it and try to check for primarility of this result, if we don't succeed, we try to find the smallest prime factor of the result and so on. I've made smallestPrimeFactor a separate function to increase readability:

fn smallestPrimeFactor nr =
(
    if isEven nr do
        return 2L
    for i = 3L to nr by 2L
        where mod nr i == 0 AND isPrime i do
            exit with i as integer64
)

fn largestPrimeFactor nr =
(
    local result
    if NOT isPrime nr then
    (
        local limit = nr/(smallestPrimeFactor nr)
        if NOT isPrime limit then
            result = largestPrimeFactor limit
        else result = limit
    )
    else result = nr
    result
)

The only thing left is to put everything together and ask for a result:

fn isEven nr = bit.and 1L nr == 0

fn isPrime nr =
(
    if isEven nr do
        return nr == 2

    local rootNr = nr^0.5
    for i in 3L to rootNr by 2L
        where mod nr i == 0 do
            return false
    true
)

fn smallestPrimeFactor nr =
(
    if isEven nr do
        return 2L
    for i = 3L to nr by 2L
        where mod nr i == 0 AND isPrime i do
            exit with i as integer64
)

fn largestPrimeFactor nr =
(
    local result
    if NOT isPrime nr then
    (
        local limit = nr/(smallestPrimeFactor nr)
        if NOT isPrime limit then
            result = largestPrimeFactor limit
        else result = limit
    )
    else result = nr
    result
)

largestPrimeFactor 600851475143L

You might have noticed one thing, there are many repetitions in the code, particularly the for loop part. Could we somehow merge it to one function instead of doing the same calculation several times? What do we actually do? In the largestPrimeFactor function we are dividing the number by prime factors and repeat the process for the result. To find the smallest prime factor we use trial division and we use the exact same trial division in the isPrime function with the only difference in the upper bound.

Well then, as we end up using trial division for all the partial results anyway, we can save ourselves some work and do it in the top-level function. We can use the vary same loop as before, checking only for odd numbers and catching even numbers at the start of the function, dividing them by two and recursing. If the number can be divided by any of the numbers in the loop, let's divide it and recurse, if not, return it:

fn isEven nr = bit.and 1L nr == 0

fn largestPrimeFactor nr =
(
    if isEven nr AND nr > 2 do
        return largestPrimeFactor (nr/2L)

    local rootNr = nr^0.5

    for i = 3L to rootNr by 2L
        where mod nr i == 0 do
            return largestPrimeFactor (nr/i)
    nr
)

largestPrimeFactor 600851475143L

DISCLAIMER: All scripts and snippets are provided as is under Creative Commons Zero (public domain, no restrictions) license. The author and this blog cannot be held liable for any loss caused as a result of inaccuracy or error within these web pages. Use at your own risk.