First triangle number to have over five hundred divisors

This was another interesting problem at Project Euler (Problem 12). Interesting because the naΓ―ve solution to this was all too trivial but slow, which forced me to seek out a better approach and I finally ended up learning something new πŸ™‚

The nth triangle number is defined as the sum of all natural numbers till n. Well, that’s definitely trivial to calculate. It’s basically the sum of first n natural numbers and can be calculated using the well known formula:Β Β 

So, all that remains is to calculate the divisors and we all know how to do that right? Just count the numbers from 2 to half (or square root, if you prefer) the triangle number that divide the triangle number. So, here’s the code I started with:

// Sum of first n natural numbers
let triangle_number (n:bigint) =
    (n * (n + 1I)) / 2I ;;
let divisors (n:bigint) =
    let rec divisor_loop (i:bigint) (a_list:bigint list) =
        if (i = 1I) then
            (i :: a_list)
        else if ((n % i) = 0I) then
            (divisor_loop (i - 1I) (i :: a_list))
            (divisor_loop (i - 1I) a_list)
    (divisor_loop (n / 2I) [n;]) ;;

As it turned out, using this method, it takes a really looooong time to even find the divisors of a large number, let alone use it for searching the triangle number with over 500 divisors.

So, I started looking for alternate methods, and finally found one here: (you’ll have to scroll to the end of the page to get to the part where they explain finding the number of divisors for a number)

The idea is simple actually. Take 5050 for example (the 100th triangle number). It’s prime factors are: 2, 5, 5 and 101. Now reduce the list to a list of unique numbers with their repetition represented in their exponents. So, the list 2, 5, 5, 101 becomes: 21, 52,1011. Now, add 1 to every exponent and multiply the resulting numbers. So, that would give us: (1+1) * (2 + 1) * (1 + 1) = 12. That’s the number of divisors that 5050 has. And, of course, that can be verified using the divisors method:

(divisors 5050I);; -> [1I; 2I; 5I; 10I; 25I; 50I; 101I; 202I; 505I; 1010I; 2525I; 5050I]

Since I’d already written a function to find the prime factors to solve Problem 3, solving this one was just a matter of writing code to search through a list of triangle numbers to find the first one with more than 500 divisors. Here’s the final code I ended up with:

// Sum of first n natural numbers
let triangle_number (n:bigint) =
    (n * (n + 1I)) / 2I ;;
let divides (n:bigint) (i:bigint) =
    n % i = 0I ;;
let is_prime (n:bigint) =
    let rec prime_check_loop (i:bigint) =
        ( i > (n / 2I) ) || ((not (i |> divides n)) && prime_check_loop (i + 1I))
    prime_check_loop 2I ;;
let rec prime_factors (n:bigint) (factors:bigint list) =
    if (is_prime n) then
        n :: factors
        let mutable i:bigint = 2I;
        while (n % i <> 0I) && (i <= n/2I) do
    for factor in factors do
        if (factor_map.ContainsKey(factor)) then
            factor_map.[factor] prime_factor_map
    let exponents = factor_map.Values
    let mutable divisor_count = 1I
    for exponent in exponents do

This took less than a minute to run and running it gives:

(triangle_number_search 500I);; -> 76576500I

So, the first triangle number to have over 500 divisors is: 76576500
(btw, it has 576 divisors)

(and as you can see, I’m still not very proficient with F# πŸ˜€ )

Sum of all the primes below two million

I’ve written about Project Euler in my previous post (which I guess was the day before yesterday). I was going pretty well with solving the problems until I reached this one – to calculate the sum of all primes below 2 million. That’s the 10th problem in Project Euler and the one I’ve had to spend the most time on till now. It’s a fairly trivial problem, if you exclude the scale that is. If you’re planning to calculate 2 million primes with your regular algorithm which checks each number for primality, then you better run it on a supercomputer πŸ˜€

I’m solving these problems in F#, and this’s what I started with:

let divides (n:bigint) (x:bigint) =
    n % x = 0I;;
let is_prime (n:bigint) =
    let rec prime_check_loop (i:bigint) =
        ( i > (n / 2I) ) || ((not (i |> divides n)) && prime_check_loop (i + 1I))
    prime_check_loop 2I;;
let rec next_prime (from:bigint) =
    if from < 2I then 2I
    else if is_prime (from + 1I) then (from + 1I)
    else next_prime (from + 1I) ;;

As a quick test, I tried calculating the next prime after 2 million and got the result fairly immediately. Happy with that, I started writing the summation code.

let sum_of_primes_below (limit:bigint) =
    let rec add_prime (n:bigint) (acc:bigint) =
        let next = next_prime n
        if next < limit then
            add_prime next (acc + next)
    add_prime 0I 0I;;

Did a quick check for the sum of primes below 10, 20, 200, 2000 and 20000. The last one took a bit of time, but that didn’t really concern me much. After all, it’s my computer doing the job right? πŸ˜€

So, ask it to calculate the sum of all primes below 2 million and left it for a minute. I switch back to the F# interactive window and see that it’s still happily calculating. I check the fsi.exe process, and it’s maxing out one of the cores completely and taking around 50M. So, I leave it for a little longer. Now, I start getting worried about the range of BigInt and whether the result would overflow. So, I quickly google around and see that BigInt is for arbitrarily large integers. Then found that the SQL Sever BigInt is actually a 64 bit int, so guess even the .net type is the same, and that maxes out at 264 – 1 = 9,223,372,036,854,775,807. hmm.. Is the sum larger than that? It’s already been around 4 minutes and I’m impatient. So I google around for the solution and found a Yahoo! Answers page which says that the sum of all primes below 2 million is 142,913,828,922. That’s good. coz, now at least I’m confident that it won’t overflow and keep looping forever.

I knew that the Sieve of Eratosthenes can be used to calculate primes, but always took it lightly. I mean, how often do you go about calculating an entire series of primes? Turns out, this is one of those times πŸ˜€

I realized that Sieve of Eratosthenes is probably the best solution for this problem, but since my code has already been running for over 10 min, I was hoping to get an answer anytime soon, so didn’t bother implementing the sieve. When I didn’t get the answer for another 20 mnutes, I had given up all hopes on my solution and I knew it was time I implemented the sieve. But just for curiosity, I tried estimating how long it’d take to get to the solution with my current method.

Assuming 1 sec per primality check, that’s about 2 million seconds = 23.1481481 days. Not good. Not good at all πŸ™ I thought that must be ridiculously wrong, since 2 million still didn’t seem like that large a number to handle. So, I tried running my code for 2000 (two thousand) and that took about a second. Then I ran it for 20000 (twenty thousand) and that took about 30 seconds. hmm.. That’s about a 30 times increase for a 10 times increase in the input size. So, for 200000 (two hundred thousand), it should take approximately 900 sec which’s 15 minutes. But well, I must’ve been estimating it completely wrong since even after 30 minutes, it was still calculating! πŸ™ But, even considering just 15 minutes for the two hundred thousand and extrapolating it, it’d have taken at least 7.5 hours! (which I did test and that estimate’s wrong too πŸ™ I left it overnight and woke up to see it still calculating.. )

So, finally I decided to generate primes with the sieve method and sum them up for the result. So, I ended up with this:

let prime_series_sieve (limit:bigint) =
    let series = List.to_array [0I..limit]
    series.SetValue(0I, 1)
    let rec eliminate_multiples (n:bigint) (i:bigint) =
        let index = (i * n) // Just to reduce as many calculations as possible
        if index < bigint.Parse(series.Length.ToString()) then
            series.SetValue(0I, (bigint.ToInt64(index)))
            eliminate_multiples n (i + 1I)
    for n in [2I..(limit/2I)] do
        eliminate_multiples n 2I
let sum_of_primes_under (limit:bigint) =
    Array.sum (prime_series_sieve limit);;

After a few quick tests for the some smaller numbers, I let it calculate the for 2 million and hurray! It gave me the right answer! It took about 2.5 minutes and had a memory footprint of around 250 M, but I was happy just to get the answer πŸ˜€

Of course it can still be improved. The wikipedia page suggests a few optimizations like Wheel Factorization. But I’m not gonna worry about that for now. I’ve got F# to learn and a whole lotta problems to solve πŸ˜€