Wednesday, April 15, 2009

Adventures in F#: Euler Problem 12

So the question is "What is the first triangle number to have over 500 divisors?"

My first attempt at solving this question included a brute force method of creating a sequence containing all of the potential divisors for a particular number, and then, using the modulus operator, filter this list down to only the items where the result of the modulus operation is zero. I tried to perform this on every item within the sequence of triangle numbers, but this took longer than I was willing to wait, so I started researching some math formulas.

It turns out that you can calculate the number of divisors for a particular number by using a rather simple formula. Here's an example.

Let's say that I have the number 36, and I want to find the number of divisors for this number. First, I would attempt to divide 36 by 2 evenly, which would result in 18. Next I would try to divide 18 by 2 evenly, which would result in 9. Next, I would find that 9 cannot be divided into two evenly, so I would then increase my attempt by 1 and try again, so I would try to divide 9 by 3, and I would find that it divides cleanly, resulting in 3. Finally, I would divide 3 by 3, and would return 1, at which point I would stop. All in all, this would look like this:

36 / 2 = 18
18 / 2 = 9
9 / 3 = 3
3 / 3 = 1

So, if you look at all of the divisors in this process, you would see that some of them would be repeated. Thus, I can use exponent notation to represent the number 36, as follows:

22 * 32 = 36

Now, you'll notice my exponents, 2 and 2. If I take these exponent values, add 1 and then multiply the two together, we'll find the number of possible divisors for 36:

(2 + 1) * (2 + 1) = 9

What are those 9 divisors?

1, 2, 3, 4, 6, 9, 12, 18, 36

So this is the formula for finding the number of divisors for any particular number. In F#, I wrote this as a Seq.unfold on a recursive function that would return the prime factorials for a particular number. I then skipped the first element (because the first element was always 2), and executed the Seq.count_by on the resulting function, thus returning a list of tuples where the first item in the tuple was the factorial, and the second item was the number of times that factorial appeared within the sequence. Next, I mapped a new sequence, retrieving the second value of the tuple and adding one to it, and lastly, I used Seq.reduce to find the product of all the items in the resulting sequence. The output was then the number of divisors for a particular number, x. This function was named divisorcount.

Next, I used Seq.unfold to create a lazy sequence of triangle numbers.

Lastly, I retrieved my result by calling Seq.first and testing each of the elements of the triangle number sequence to determine if the return value was over 500.

So finally, without further ado, here's my code. It returns in a little over 2 seconds.

module Problem12

let divisorcount x =
if x = 1 then 0
let rec div (x:int * int) =
if (fst x = 0) then None
else if (snd x = 1) then Some(0, snd x)
else match (snd x % fst x) with
| 0 -> Some(fst x, snd x / fst x)
| _ -> div (fst x + 1, snd x)
Seq.unfold (fun x ->
match (div x) with
| None -> None
| Some(o) -> Some(fst x, o)) (2, x)
|> Seq.skip 1
|> Seq.count_by (fun x -> x)
|> (fun x -> snd x + 1)
|> Seq.reduce (fun x prod -> x * prod)

let triangles =
Seq.unfold (fun x -> Some(snd x + fst x, (snd x + fst x, snd x + 1))) (0, 1)

let result =
|> Seq.first (fun x -> if divisorcount x > 500 then Some(x) else None)
|> Option.get

let print = divisorcount 9 |> printfn "%A"