Faffing round the edges Its easier than getting stuck in

Euler 7: No sieve ‘cos I can’t see how.

Project Euler's Problem 7 is back to the meaty business of generating primes. Now I know that lots of languages can do this for you but where is the fun of that? Likewise I could probably look up an algorithm on Wikipedia but only after I have an answer.

Answer three gave me a prime checker that I just reused here. I need to find a faster way to do this later but for now it seems adequate.

I tidied my prime checker to make it easier (for me) to read:

%%
%test for primality
%%
is_prime2(N) when N =:= 2 -> true;
is_prime2(N) when N rem 2 =:= 0 -> false;
is_prime2(N) ->
	Limit = round(math:sqrt(N) + 1),
	is_prime_tester(N, Limit, 3).

is_prime_tester(N, Limit, CandidateFactor) when CandidateFactor < Limit ->
	case N rem CandidateFactor =:= 0 of
		true -> false;
		false ->
			is_prime_tester(N, Limit, CandidateFactor+2)
	end;
is_prime_tester(_N, _Limit, _CandidateFactor) ->
	true.

With that to reuse the rest of problem 7 is just:

prob7(N) ->
	gen_primes(N, [2], 3).

gen_primes(N, Primes, _Candidate) when N =:= length(Primes) -> Primes;
gen_primes(N, Primes, Candidate) ->
	case is_prime2(Candidate) of
		true ->
			gen_primes(N, [Candidate|Primes], Candidate+2);
		false ->
			gen_primes(N, Primes, Candidate+2)
		end.

A colleague at work has a faster Java version (and he is my new hero) but I love hwo quick (and easy) it is for me to knock up an Erlang one. It is such a readable, expressive language.

Euler 6: At last I find it easy!

Project Euler problem 6:

Find the difference between the sum of the squares of the first one hundred natural numbers and the square of the sum.

Ok...this makes sense to me. No special, magic formulae in a foreign script, no head twisting numbers that remain inscrutable even with 1000s of digits (primes, WTF?) Just something that lends itself to learning the functional programming paradigm of programming (I said paradigm, in a blog, check me).

Enough. My solution (for what it is worth):

prob6(N) ->
	L = lists:seq(1, N),
	square_of_sum(L) - sum_of_squares(L).

square_of_sum(L) ->
	Sum = lists:sum(L),
	Sum * Sum.

sum_of_squares(L) ->
	lists:foldl(fun(X, Acc) -> Acc + (X*X) end, 0, L).

I am sure there is a better way and if I spot it I will post it but I'm pretty pleased with that for a beginner in Erlang and maths.

Euler 5. Cop out

Probelm 5 from Project Euler didn't float my boat.

2520 is the smallest number that can be divided by each of the numbers from 1 to 10 without any remainder.

What is the smallest number that is evenly divisible by all of the numbers from 1 to 20?

I just fiddled with my calculator for a bit and this is what I ended up with:

prob5() ->
	N = (2 * 3 * 4 * 5 * 6 * 7 * 11 * 13 * 17 * 19),
	Factors = euler_too:factorise( N ),
	io:format("~p has factors ~p~n", [N, Factors]).

So...I learned nothing...Or I learned about the Lowest Common Multiple and Prime Factors but I have yet had the time to go back and create a function that can calculate the lcm for any set of numbers. The Euler solution pdf that one gains access to with a correct answer is going to help me with that I am sure.

Like all software, one you have a working solution it is very hard to go back and make it better since you move ever forward so maybe I should take the time to get it right next time. I'm doing it for fun after all.

Euler 4: The long and the short

I have been quiet as I have been very busy at work and my garage got broken into and my bikes stolen but I have also done the Euler's up to 10 now. Here is two versions of problem 4:

A palindromic number reads the same both ways. The largest palindrome made from the product of two 2-digit numbers is 9009 = 91 99.

Find the largest palindrome made from the product of two 3-digit numbers.

Sounds simple? Well my first attempt I made a proper meal of it. Here is the long listing:

prob4() ->
	Candidates = lists:seq(998001, 1000, -1),
	C = lists:filter(fun(X) -> is_palindrome(X) end, Candidates),
	PalinFacs = lists:foldl(fun(X, Acc) -> [{X,factorise(X)}|Acc] end, [], C),
	highest_palindrome(lists:reverse(PalinFacs)).

highest_palindrome([{Palindrome, Factors}|T]) ->
	ThreeDigitFactors = lists:filter(fun(X) -> is_three_digit(X) end, Factors),
	TwoThreeDigitFactors = lists:filter(fun(X) -> is_three_digit(Palindrome div X) end, ThreeDigitFactors),
	case length(TwoThreeDigitFactors) =/= 0 of
		true ->
			{Palindrome, lists:zipwith(fun(X, Y) -> {Palindrome div X, Palindrome div Y} end, TwoThreeDigitFactors, lists:reverse(TwoThreeDigitFactors))};
		false ->
			highest_palindrome(T)
	end.

is_three_digit(N) when N > 99, N < 1000 -> true;
is_three_digit(_N) -> false.

is_palindrome(N) ->
	A = integer_to_list(N),
	B = lists:reverse(A),
	(list_to_integer(A) =:= list_to_integer(B)).

I mean what was I thinking? What is all that "TwoThreeDigitFactors" stuff? I must have been thinking about something 'cos there is code there. I even use the factorising function from Euler 3 which is handy. So it got the answer so it "works". And fast enough too (well under second). I was trying to avoid the brute force of just looping around 100-999 twice and then looking for palindromes.

The right answer gave me access to the forum where a Haskell solution led me to this c'n'p erlang version:

prob4_alt() ->
	lists:max([A*B || A <- lists:seq(100,999), B <- lists:seq(100,999), is_palindrome(A*B)]).

Which is quicker too. Proves the Knuth maxim:

Premature optimization is the root of all evil

Ah. You live and learn. And you learn best through experience, innit?!

Euler Problem 3

Project Euler problem 3 was a real challenge for me. I don't have much maths education (though I am getting more all the time) so the factoring of a number to find it's prime factors meant learning a lot (from wikipedia mainly).

One thing I noticed in my code compared to a lot of the solution on Project Euler is that it is longer. It is also less specific. I have a Prime test and a factorization function that I use (and have re-used on euler 4 and 5).

Anyway the code.

problem_3(N) ->
	lists:max(lists:filter(fun(X) -> is_prime2(X) end, factorise(N))).

This gets the largest prime number from the list of factors for N (in the case of the problem 600851475143)

This function gets the factors for N:-

factorise(N) ->
	Limit = round(math:sqrt(N)+1),
	factorise(N, Limit, 1, []).

factorise(N, Limit, Limit, Factors) ->
	case length(Factors) =:= 0 of
		true -> [1,N];
		false ->
			lists:sort(lists:merge([N div A || A <- Factors], Factors))
	end;
factorise(N, Limit, Candidate, Factors) ->
	case N rem Candidate =:= 0 of
		true -> factorise(N, Limit, Candidate +1, [Candidate|Factors]);
		false -> factorise(N, Limit, Candidate +1, Factors)
	end.

I'm sure the mathematicians know a more efficient factoring algorithm and I hope to learn one from the solution pdf on euler when I have the time (after this weekends Bug-o-rama maybe).

This code is my prime test:-

is_prime2(N) ->
	case N =:= 2 of
		true -> true;
		false ->
			case N rem 2 =:= 0 of
				true -> false;
				false->
					Limit = round(math:sqrt(N) + 1),
					is_prime_tester(N, Limit, 3)
			end
		end.

is_prime_tester(N, Limit, CandidateFactor) when CandidateFactor < Limit ->
	case N rem CandidateFactor =:= 0 of
		true -> false;
		false ->
			is_prime_tester(N, Limit, CandidateFactor+2)
	end;
is_prime_tester(_N, _Limit, _CandidateFactor) ->
	true.

This is like the factoring algorithm really. It is recursive and tests if a number can be divided. It misses out all the even numbers as there is no point in testing against them after testing against 2. Although, to be fair I could just factor the number and if it only has two factors (itself and 1) then it is prime but this was a tad faster. I think I need to look at refactoring my code in the light of how short some of my colleagues code is (see Chris's python solution on his blog) although this erlang effort still performs well for numbers 3 orders of magnitude larger.

I really enjoyed this problem and i am sure there are more optimal ways to test for primality and to factor. This example on the interweb is blazingly fast for even biggish numbers.

Project Euler problem 1 parallel processing

Ok, so I have managed problem 3 finally but before I write about that I thought I'd write about my parallel version of problem 1.

This was just an exercise for me to learn Elrang so I am sure it is sub-optimal in many ways. The bottleneck in my solution is the generation of the terms in the line:

lists:sum(lists:seq(0, N, 3)) +
    lists:sum(lists:seq(0,N, 5)) -
    lists:sum(lists:seq(0, N, 15)).

So I thought the best thing to do was parallelize the term generation.

gen_term() ->
	receive
		{From, {gen, Limit, Incr, Action}} ->
			From ! {self(), lists:sum(lists:seq(0, Limit, Incr)), Action}
	end.

The process dies as soon as it has done its job. All this process does is calculated a sequence from 0 to Limit in steps of Incr (EG the sequence from 0, 999 in steps of 3). It passes through the value of Action so that another process I use for collecting and acting on the terms can process the result. The Pid of the collecting process is passed as From.

This process:

prob_1_accumulater(Sum) ->
	receive
		{From, Amount, Action} ->
			case Action of
		 		add -> prob_1_accumulater(Sum + Amount);
				sub -> prob_1_accumulater(Sum - Amount)
			end
	end.

Is the process that performs the calculation on the terms provided by the gen_term process. When the process is started Sum is set to zero. Then depending on the value of Action the process restarts itself with the result of Sum [+/-] Amount. The atoms add and sub are the keys to the action.

The function that controls all this is:

prob_1_para(N) ->
	AccumPid = spawn(?MODULE, prob_1_accumulater,[0]),
	Pid3 = spawn(fun gen_term/0),
	Pid5 = spawn(fun gen_term/0),
	Pid15 = spawn(fun gen_term/0),
	Pid3 ! {AccumPid, {gen, N, 3, add}},
	Pid5 ! {AccumPid, {gen, N, 5, add}},
	Pid15 ! {AccumPid, {gen, N, 15, sub}}.

All that happens here is that a process is created for calculating the terms and a process is created for each term we wish to generate. The prob_1_accumulater process prints the final total when all the terms have been gathered.

I am sure that there is a more elegant solution to this but the code above is a revelation to a poor old Java developer. To create what amounts to a 3 server instances for generating terms and one server instance to collect them in so few lines of code is amazing. The code above does run faster when I start Erlang's beam in SMP mode:

erl -smp

Running the sequential solution for 99999999 takes about 27-30 seconds. Running the parallel function takes 18-20 seconds. Not a massive overall improvement but still good. It is great to see the monitors for both cpu cores spring to life. With 3 cpus I guess it would be faster still. Still very much enjoying Erlang.

Erlang for Project Euler problem 2.

Ok. So now I am looking at problem 2 from Project Euler.

Each new term in the Fibonacci sequence is generated by adding the previous two terms. By starting with 1 and 2, the first 10 terms will be:

1, 2, 3, 5, 8, 13, 21, 34, 55, 89, ...

Find the sum of all the even-valued terms in the sequence which do not exceed four million

Fibonacci series? I have heard of it...Quick look at wikipedia and I at least know what it is. My rather brute force approach is (thanks mainly to the rather excellent LiteratePrograms.org. There was actually a (purportedly) more optimised version on the site but since I did not understand it there was nothing I could do to appropriate it):-

prob_2() ->
	L = lists:filter(fun(X) -> (X rem 2 =:= 0) end, fib(33)),
	lists:sum(L).

%used by the problem 2
fib_acc( 0, _Result, _Next, List) -> List;
fib_acc( Iter, Result, Next, List) when Iter > 0 ->
        fib_acc( Iter -1, Next, Result + Next,[Next|List]).

fib(N) -> fib_acc( N, 0, 1,[]).

I guess this is rather sub optimal since I have to first run the Fib sequence function a few times and figure out which Fib number is the closest to 4 million. Then there is the hack where I have hardcoded that number (33). But it gets me the answer so I get to learn loads more by reading the forum.

I am really let down by my lack of maths at school. Sure some of that is my fault for being a cocky, arrogant jerk at 13 but then again the teachers should have been used to that.

Ok, so after reading a bit I have these two more optimal approaches:-

prob_2_alt(Limit) ->
	summer(Limit, 0, 1, 1).

summer(Limit, Sum, _A, B) when B >= Limit -> Sum;
summer(Limit, Sum, A, B) when B rem 2 =:= 0->
	summer(Limit, Sum+B, B, A+B );
summer(Limit, Sum, A, B) ->
	summer(Limit, Sum, B, A+B).

And this (which is just a novice's erlanging of the Euler pdf on the problem):-

prob_2_alt2(Limit) ->
	summer2(Limit, 0,  1, 2).

summer2(Limit, Sum,  _B, C) when C >= Limit -> Sum;
summer2(Limit, Sum,  B, C) ->
	summer2(Limit, Sum+C,  C+B+C, C+B+C+B+C).

Not sure how much I have learned here except that I have had my admiration of Erlang further stoked. I understand the more elegant Euler solutions and I'm pleased I was able to brute force my way in to get an answer and learn a lot more about the problem. Thankfully I understand the answers completely so feel satisfied if a tad dim.

I am not sure how pretty my Java Impl of the above would be. Might try later...then again, I might not.

Erlang for Project Euler problem 1.

I saw this post and it made me have a crack at Project Euler problem 1.

If we list all the natural numbers below 10 that are multiples of 3 or 5, we get 3, 5, 6 and 9. The sum of these multiples is 23.

Find the sum of all the multiples of 3 or 5 below 1000.

I  thought I'd better have a go before reading the post. My first crack was a bit long:-

-module (sum_3_5).

-compile(export_all).

sum_3_5() ->
	L = lists:seq(1, 999),
	lists:sum(lists:filter(fun worker/1, L)).

worker(Elem) ->
	by_3(Elem) or by_5(Elem).		

by_3(Elem) ->
	case Elem rem 3 of
		0 -> true;
		_ -> false
	end.

by_5(Elem) ->
	case Elem rem 5 of
		0 -> true;
		_ -> false
	end.

But it worked. Second attempt (2 minutes later) is much more terse:-

prob_1(N) ->
	lists:sum([A || A <- lists:seq(1,N), (A rem 3 =:=0) or (A rem 5 =:= 0)]).

I had one more go at it to try and speed it up for large values of N:-

prob_1_alt(N) ->
	lists:sum(lists:seq(0, N, 3)) +
        lists:sum(lists:seq(0,N, 5)) -
        lists:sum(lists:seq(0, N, 15)).

Which was a lot faster for large numbers. Is it worth creating a process to calculate each term in the above equation? Will that speed it up more for very large values of N? I might try later. But I'll probably try Project Euler 1 in Squeak, Clojure and Scheme first. And problem 2 in Erlang.