I think I'll post my answer to Problem 12 of Project Euler next:

What is the value of the first triangle number to have over five hundred divisors?

What exactly is a triangle number? It is a number given by this formula:

For this problem we actually want to generate a sequence of triangular numbers. It is most efficient to generate the n^th number by adding n to the (n-1)^th number. We can capture that in F# as follows:

let triangles =

    let rec inner m n =

        seq {

            let p = m + n

            yield p

            yield! inner p (n+1L)

        }

    seq {

        yield! inner 0L 1L

    }

The next question is how do we find the number of divisors of the number n? The brute force way would be to try all values up to n/2 and find the ones which divide evenly into the target number. I figured there had to be a better way. A little bit of research on the web showed me that if you know the prime factorization for a number then you can easily compute the number of divisors. Suppose that the a number n can be factored as follows:

where ω(n) is the number of distinct prime factors of n. Then the number of divisors of n is given by:

Fortunately, I'm no longer a math major and so I don't have to prove that this is true. For my purposes it's enough to just try it out and see if it seems reasonable. So, we can say that 12 = 2² x 3. According to the formula then the number of divisors should be (2 + 1) x (1 + 1) = 6. Listing the divisors we find the following set: {1, 2, 3, 4, 6, 12}. Good enough for me.

To solve this Project Euler problem I need to be able to factor an arbitrary number. The simplest thing is to try all the primes up to the square root of the target number and see if they divide into it evenly. As before I need a sequence of primes. Here's the code:

let primes =

    let rec sieve n table =

        let rec update n p table =

            let n = n+p

            match Map.tryFind n table with

            | None              -> Map.add n p table

            | Some(prime)       -> update n p table

        seq {

            match Map.tryFind n table with

            | None          ->

                yield n

                yield! sieve (n+1L) (table |> Map.add (n*n) n)

            | Some(prime)   ->

                yield! sieve (n+1L) (Map.remove n table |> update n prime)

        }

    sieve 2L Map.empty

I think this will be the last time I display this code. From now on just assume that when needed I can generate a list of primes with a lazy sequence.

I also need a few helper functions. Here's an integer square root:

let root (n:int64) =

    int64 (sqrt (float n))

And an integer power function (tail recursive just for fun):

let power p e  =

    let rec inner p e acc =

        if e = 0L then

            acc

        else

            inner p (e-1L) p*acc

    inner p e 1L

Finally here's a function to compute the number of times a prime divides into a number:

let exponent p n =

    let rec inner p n acc =

        if n%p = 0L then

            inner p (n/p) (acc+1L)

        else

            acc

    inner p n 0L

The next step is to compute a sequence which consists of the prime factors of a number along with each factor's exponent. Here's what I came up with:

let rec factors n =

    seq {

        if n > 1L then

            let divisor =

                primes

                    |> Seq.takeWhile (fun k -> k <= root n)

                    |> Seq.tryFind (fun k -> n%k = 0L)

            match divisor with

            | None      -> yield (n, 1L)

            | Some(k)   ->

                let e = exponent k n

                yield (k, e)

                yield! factors (n/power k e)

    }

There's some inefficiency here in that each time we recurse we start the sequence of primes over again. I'm not losing any sleep over it...

Given all this we can now compute the number of divisors of a number:

let divisors n =

    factors n

        |> Seq.map (fun k -> snd k)

        |> Seq.fold (fun acc e -> acc*(e+1L)) 1L

The final step is to find the first entry in the triangular number sequence which has more than 500 divisors.

triangles

    |> Seq.tryFind (fun n -> divisors n > 500L)

    |> printfn "%A"

It was nice to revisit this solution while putting together this post. My original solution was a bit hard to read. I think it is much clearer now.