Sunday, January 18, 2009

Erlang solution for Euler problem 12

In project Euler problem 12 we are asked to find the value of the first triangle number to have over five hundred divisors.

I implemented a brute force approach which generates all triangle numbers. For each triangle number T, the program computes the number of divisors. When it finds the first triangle number to have more than N = 500 divisors, the program terminates.

My first implementation (the code is not show below) computed the number of divisors of a given triangle number T by trying to divide it by all numbers from 1 tru T. This allowed me to solve the problem, but the run-time was much more than the allowed one minute.

My second implementation (see "approach 1" below) only tried all numbers from 1 to sqrt(T). This is based on the observation that if T = x*y then one of the two factors (x or y) must be below sqrt(T) and and the other must be above sqrt(T). Thus, if we find a factor x below sqrt(T), we know there are two factors because there must be a corresponding factor y above sqrt(T). The code must be careful to count sqrt(T) only once as a factor if T happens to be a square number. This approach solved the problem in 2.8 seconds for N = 500.

Once I had submitted my solution to the project Euler website and got access to the problem notes, I discovered the existence of a more efficient method to compute the number of factors of a number which is called the Tau function. My third implementation (see "approach 2" below) used this approach. This approach was almost twice as fast when the required number of divisors N is 500. For greater number of divisors, the difference between approach 1 and 2 was greater as shown in the following graph.



The code for approach 1:


%% Project Euler, problem 12
%%
%% What is the value of the first triangle number to have over five hundred divisors?

-module(euler12b).
-author('Cayle Spandon').

-export([solve/0, solve/1, count_divisors/1]).

% Try only the divisors between 1 and sqrt(Number).
%
count_divisors(Number) ->
Try = 1,
MaxTry = trunc(math:sqrt(Number)),
Count = 0,
count_divisors(Number, Try, MaxTry, Count).

count_divisors(Number, Try, MaxTry, Count) ->
if
Try > MaxTry ->
Count;
Number rem Try == 0 ->
if
Try == MaxTry ->
Found = 1;
true ->
Found = 2 % every divisor below sqrt(Number) has a counter-part above it
end,
count_divisors(Number, Try + 1, MaxTry, Count + Found);
true ->
count_divisors(Number, Try + 1, MaxTry, Count)
end.

solve() ->
solve(500).

solve(WantedFactors) ->
solve(3, 3, 0, WantedFactors).

solve(Try, Step, BestSoFar, WantedFactors) ->
Count = count_divisors(Try),
if
Count > BestSoFar ->
NewBestSoFar = Count;
true ->
NewBestSoFar = BestSoFar
end,
case Count > WantedFactors of
true ->
Try;
false ->
solve(Try + Step, Step + 1, NewBestSoFar, WantedFactors)
end.


The code for approach 2:


%% Project Euler, problem 12
%%
%% What is the value of the first triangle number to have over five hundred divisors?
%%
%% Implementation uses the Tau function.
%% See http://mathschallenge.net/index.php?section=faq&ref=number/number_of_divisors for details.

-module(euler12c).
-author('Cayle Spandon').

-export([solve/0, solve/1, prime_factors/1, count_divisors/1, make_count_list/1]).

prime_factors(Number) ->
prime_factors(Number, 2).

prime_factors(Number, Try) ->
if
Try > Number ->
[];
Number rem Try == 0 ->
[Try | prime_factors(Number div Try, Try)];
true ->
prime_factors(Number, Try + 1)
end.

make_count_list(List) ->
make_count_list(List, []).

make_count_list([], Results) ->
lists:reverse(Results);

make_count_list(List = [Head | _Tail], Results) ->
{RemainingList, Count} = count_nr(List, Head),
make_count_list(RemainingList, [Count | Results]).

count_nr(List = [Head | Tail], Nr) ->
if
Head == Nr ->
{RemainingList, Count} = count_nr(Tail, Nr),
{RemainingList, Count + 1};
true ->
{List, 0}
end;

count_nr([], _Nr) ->
{[], 0}.

count_divisors(Number) ->
CountList = make_count_list(prime_factors(Number)),
lists:foldl(fun(Nr, Acc) -> Acc * (Nr + 1) end, 1, CountList).

solve() ->
solve(500).

solve(WantedFactors) ->
solve(1, 2, WantedFactors).

solve(Try, Step, WantedFactors) ->
case count_divisors(Try) > WantedFactors of
true ->
Try;
false ->
solve(Try + Step, Step + 1, WantedFactors)
end.

No comments: