[erlang-questions] LAST RESUME of Pytagoriam Numbers!!!

Ivan Carmenates García

unread,

Jan 15, 2006, 4:11:42 AM1/15/06

to erlang-q...@erlang.org

Compiled using c(py) normal in the shell.

Feel free to make a new Resume on a DualCore CPU!!!

Now we are re-talking about!!!

Tested on Intel(R) Celeron CPU 2.13 GHz, 192 of Ram

*** Pytagoriam's Numbers REPORT! ***

Using N = 300

*** THE ORIGINAL IMPLEMENTATIONS (ME!)****
-py1: 14860000 mcs = 14.86 s
-py2: 14843999 mcs = 14.843999 s
-py3: 12749999 mcs = 12.749999 s

*** Tony's implementations ****
-pythag1: 11702999 mcs = 11.702999 s
-pythag2: 1983999 mcs = 1.983999 s

*** Hynek's implementation ****
Simpler and about 5% faster version
-pythag3: 1952999 mcs = 1.952999 s

*** Edmond's implementation, using parallelism****
-py2E: 14843999 mcs = 14.843999 s

*** Willem's implementation****
-wpythag2: 1749999 mcs = 1.749999 s

*** Hynek's new implementation****
-pythag4: 1140999 mcs = 1.140999 s

*** Willem's new implementation in parallel by Edmond!!!****
-wpythag2P: 1795999 mcs = 1.795999 s

*** Morten's implementation****
-pythag5: 46999 mcs = 0.046999 s

*** Richard's improvement****
-py3R: 46999 mcs = 0.046999 s *** NEW ***

*** Joe's improvement****
-py3a: 1577999 mcs = 1.577999 s *** NEW ***

Comparisons in results agains 'pythag1' Tony's
function
py1 returns the same than 'pythag1'
py2 returns the same than 'pythag1'
py3 returns the same than 'pythag1'
pythag1 returns the same than 'pythag1'
pythag2 returns the same than 'pythag1'
pythag3 returns the same than 'pythag1'
py2E returns the same than 'pythag1'
wpythag2 returns the same than 'pythag1'
pythag4 returns the same than 'pythag1'
wpythag2P returns the same than 'pythag1'
pythag5 returns the same than 'pythag1'
py3R returns the same than 'pythag1' <- now are we talking about!!! (is
the best!!!)
py3a returns the same than 'pythag1'

NOTE: Time took using timer:tc/3 function.

%******************* START CODE *********************

-module(py).
-compile([export_all]).

start(N)->
%% My implementations.
{T1,R1} = timer:tc(py,py1, [N]),
{T2,R2} = timer:tc(py,py2,[N]),
{T3,R3} = timer:tc(py,py3,[N]),

%% Tony's improvement of the original form 3.
{T4,R4} = timer:tc(py,pythag1,[N]),
%% Tony's implementation.
{T5,R5} = timer:tc(py,pythag2,[N]),

%% Hynek's implementation.
%% Simpler and about 5% faster version:
{T6,R6} = timer:tc(py,pythag3,[N]),

%% Edmond's implementation using parallelism.
{T7,R7} = timer:tc(py,py2E,[N]),

%% Willem's implementation.
{T8,R8} = timer:tc(py,wpythag2,[N]),

%% Hynek's new version
{T9,R9} = timer:tc(py,pythag4,[N]),

%% Willem's implementation in parallel by Edmond.
{T10,R10} = timer:tc(py,wpythag2P,[N]),

%% Morten's implementation.
{T11,R11} = timer:tc(py,pythag5,[N]),

%% Richard's improvement.
{T12,R12} = timer:tc(py,py3R,[N]),

%% Joe's improvement.
{T13,R13} = timer:tc(py,py3a,[N]),

io:format("~n *** Pytagoriam's Numbers REPORT! ***~n~n"),
io:format("Using N = ~p~n", [N]),

io:format("~n*** THE ORIGINAL IMPLEMENTATIONS (ME!)****~n"),
io:format("-py1: ~p mcs = ~p s~n", [T1,T1/1000000]),
io:format("-py2: ~p mcs = ~p s~n", [T2,T2/1000000]),
io:format("-py3: ~p mcs = ~p s~n", [T3,T3/1000000]),

io:format("~n*** Tony's implementations ****~n"),
io:format("-pythag1: ~p mcs = ~p s~n", [T4,T4/1000000]),
io:format("-pythag2: ~p mcs = ~p s~n", [T5,T5/1000000]),

io:format("~n*** Hynek's implementation ****~n"),
io:format("Simpler and about 5% faster version~n"),
io:format("-pythag3: ~p mcs = ~p s~n", [T6,T6/1000000]),

io:format("~n*** Edmond's implementation, using parallelism****~n"),
io:format("-py2E: ~p mcs = ~p s~n", [T7,T7/1000000]),

io:format("~n*** Willem's implementation****~n"),
io:format("-wpythag2: ~p mcs = ~p s~n", [T8,T8/1000000]),

io:format("~n*** Hynek's new implementation****~n"),
io:format("-pythag4: ~p mcs = ~p s~n", [T9,T9/1000000]),

io:format("~n*** Willem's new implementation in parallel by
Edmond!!!****~n"),
io:format("-wpythag2P: ~p mcs = ~p s~n", [T10,T10/1000000]),

io:format("~n*** Morten's implementation****~n"),
io:format("-pythag5: ~p mcs = ~p s~n", [T11,T11/1000000]),

io:format("~n*** Richard's improvement****~n"),
io:format("-py3R: ~p mcs = ~p s~n", [T12,T12/1000000]),

io:format("~n*** Joe's improvement****~n"),
io:format("-py3a: ~p mcs = ~p s~n", [T13,T13/1000000]),

io:format("~nComparisons in results agains 'pythag1' Tony's
function~n"),
Rs = [{py1,R1}, {py2,R2}, {py3,R3}, {pythag1,R4}, {pythag2,R5},
{pythag3,R6},
{py2E,R7}, {wpythag2,R8}, {pythag4,R9}, {wpythag2P,R10},
{pythag5,R11}, {py3R,R12},
{py3a,R13}],
lists:foreach(fun({Name, R})->
if
(R=:=R4)->
io:format("~p returns the same than 'pythag1'~n", [Name]);
true->
io:format("~p NO returns the same than 'pythag1'~n",
[Name])
end
end, Rs),
io:format("~nNOTE: Time took using timer:tc/3 function.~n").

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% The original form 1.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
py1(Max)->
L = lists:seq(1, Max),
lists:foldr(
fun(A, Acc3)->
lists:foldr(
fun(B, Acc2)->
lists:foldr(
fun(C, Acc)->
case ((A*A + B*B =:= C*C) andalso (A+B+C =<
Max)) of
true->
[{A,B,C}|Acc];
false->
Acc
end
end
, Acc2, L)
end
, Acc3, L)
end
, [], L).

%% The original form 2.
py2(Max)->
lists:reverse(fora(1, [], Max)).

fora(A, Acc, Max)->
Acc1 = forb(A,1, Acc, Max),
case A< Max of
true->
fora(A+1, Acc1, Max);
false->
Acc1
end.

forb(A,B, Acc, Max)->
Acc1 = forc(A,B,1, Acc, Max),
case B< Max of
true->
forb(A,B+1, Acc1, Max);
false->
Acc1
end.

forc(A,B,C, Acc, Max)->
Acc1 = case (A*A + B*B =:= C*C) andalso (A+B+C =< Max) of
true->
[{A,B,C}|Acc];
_->
Acc
end,
case C< Max of
true->
forc(A,B,C+1, Acc1, Max);
false->
Acc1
end.

%% The original form 3.
py3(Max)->
[{A,B,C} ||
A<-lists:seq(1, Max),
B<-lists:seq(1, Max),
C<-lists:seq(1, Max),
A*A + B*B =:= C*C,
A+B+C =< Max].

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Tony's improvement of the original form 3.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
pythag1(N) ->
L = lists:seq(1,N),
[ {A,B,C} ||
A<- L,
B<- L,
C<- L,
A+B+C =< N,
A*A+B*B =:= C*C].

%% Tony's implementation.
pythag2(N) ->
lists:reverse(pythan2_A(1, N, [])).

pythan2_A(A, N, Acc) when A> N -> Acc;
pythan2_A(A, N, Acc) -> pythan2_A(A+1,N,pythan2_B(A, 1, N, Acc)).

pythan2_B(A, B, N, Acc) when A+B> N -> Acc;
pythan2_B(A, B, N, Acc) -> pythan2_B(A,B+1,N,pythan2_C(A, B, 1, N, Acc)).

pythan2_C(A, B, C, N, Acc) when A+B+C> N -> Acc;
pythan2_C(A, B, C, N, Acc) ->
if A*A+B*B =:= C*C ->
pythan2_C(A, B, C+1, N, [{A,B,C}|Acc]);
true ->
pythan2_C(A, B, C+1, N, Acc)
end.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Hynek's implementation.
%% Simpler and about 5% faster version:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
pythag3(N) when is_integer(N) -> pythag3(N,1).

pythag3(N, A) when A+2> N -> [];
pythag3(N, A) -> pythag3(N, A, 1).

pythag3(N, A, B) when A+B+1> N -> pythag3(N, A+1);
pythag3(N, A, B) -> pythag3(N, A, B, 1).

pythag3(N, A, B, C) when A+B+C> N -> pythag3(N, A, B+1);
pythag3(N, A, B, C) when A*A + B*B =:= C*C -> [{A, B, C}|pythag3(N, A,
B, C+1)];
pythag3(N, A, B, C) -> pythag3(N, A, B, C+1).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Edmond's implementation using parallelism.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%---- START CODE ----

py2E(Max)->
lists:flatten(lpmap(fun(A) ->
forbE(A, 1, [], Max)
end, lists:seq(1, Max), ordered)).

forbE(A, B, Acc, Max) ->
Acc1 = forcE(A, B, 1, Acc, Max),
case B< Max of
true -> forbE(A, B+1, Acc1, Max);
false -> Acc1
end.

forcE(A, B, C, Acc, Max) ->
Acc1 = case (A*A + B*B =:= C*C) andalso (A+B+C =< Max) of
true -> [{A,B,C}|Acc];
_ -> Acc
end,
case C< Max of
true-> forcE(A, B, C+1, Acc1, Max);
false-> Acc1
end.

pythag2E(N)->
lists:flatten(lpmap(fun(A) ->
pythan2_BE(A, 1, N, [])
end, lists:seq(1, N), ordered)).

pythan2_AE(A, N, Acc) when A> N -> Acc;
pythan2_AE(A, N, Acc) -> pythan2_AE(A+1,N,pythan2_BE(A, 1, N, Acc)).

pythan2_BE(A, B, N, Acc) when A+B> N -> Acc;
pythan2_BE(A, B, N, Acc) -> pythan2_BE(A,B+1,N,pythan2_CE(A, B, 1, N,
Acc)).

pythan2_CE(A, B, C, N, Acc) when A+B+C> N -> Acc;
pythan2_CE(A, B, C, N, Acc) ->
if A*A+B*B =:= C*C ->
pythan2_CE(A, B, C+1, N, [{A,B,C}|Acc]);
true ->
pythan2_CE(A, B, C+1, N, Acc)
end.

%% @spec lpmap(fun(), list(), (atom() = ordered|unordered)) -> list()
%% @doc Spawns a process for each element in list L, performs specified
%% function F against each in parallel and then returns results
either
%% same order as L (ordered) or in any order (unordered).
%% NB: See also lpmap/4.

lpmap(F, L, ordered) ->
Ref = erlang:make_ref(),
Pids = [lpmap_spawn_link(self(), Ref, F, I) || I<- L],
lpmap_gather_ordered(Pids, Ref, [], 0, void);
lpmap(F, L, unordered) ->
Ref = erlang:make_ref(),
lists:foreach(fun(I) ->
lpmap_spawn_link(self(), Ref, F, I)
end, L),
lpmap_gather_unordered(length(L), Ref, [], 0, void).

%% @spec lpmap(fun(), integer(), list(), (atom() =
ordered|unordered)) ->
list()
%% @doc Same as lpmap/3 except ensures only a maximum of MaxPs parallel
%% processes execute function F at any one time (i.e. first takes
MaxPs
%% items from list, executes F in parallel against each, then as
each
%% process returns, spawns another process on next item in L as
long as
%% active processes are less than MaxPs).
%% NB: See also lpmap/4.

lpmap(F, L, MaxPs, ordered) when MaxPs>0 ->
Ref = erlang:make_ref(),
{HPids, TPids} = if
length(L)> MaxPs -> lists:split(MaxPs, L);
true -> {L, []}
end,
Pids = [lpmap_spawn_link(self(), Ref, F, I) || I<- HPids],
lpmap_gather_ordered(Pids, Ref, TPids, MaxPs, F);
lpmap(F, L, MaxPs, unordered) when MaxPs>0 ->
Ref = erlang:make_ref(),
{HPids, TPids} = if
length(L)> MaxPs -> lists:split(MaxPs, L);
true -> {L, []}
end,
lists:foreach(fun(I) ->
lpmap_spawn_link(self(), Ref, F, I)
end, HPids),
lpmap_gather_unordered(length(HPids), Ref, TPids, MaxPs, F).

%% lpmap internal functions

lpmap_spawn_link(Parent, Ref, F, I) ->
spawn_link(fun() ->
Parent ! {self(), Ref, F(I)}
end).

lpmap_gather_ordered([], _Ref, [], _MaxPs, _F) ->
[];
lpmap_gather_ordered([HPid|TPids], Ref, L, MaxPs, F) ->
receive
{HPid, Ref, Ret} when length(TPids)<MaxPs, L=/=[] ->
[H | T] = L,
[Ret | lpmap_gather_ordered(
lists:append(TPids, [lpmap_spawn_link(self(), Ref, F, H)]),
Ref, T, MaxPs, F)];
{HPid, Ref, Ret} ->
[Ret | lpmap_gather_ordered(TPids, Ref, L, MaxPs, F)]
end.

lpmap_gather_unordered(0, _Ref, [], _MaxPs, _F) ->
[];
lpmap_gather_unordered(NPs, Ref, L, MaxPs, F) ->
receive
{_Pid, Ref, Ret} when NPs-1<MaxPs, L=/=[] ->
[H | T] = L,
lpmap_spawn_link(self(), Ref, F, H),
[Ret | lpmap_gather_unordered(NPs, Ref, T, MaxPs, F)];
{_Pid, Ref, Ret} ->
[Ret | lpmap_gather_unordered(NPs-1, Ref, L, MaxPs, F)]
end.

%%---- END CODE -----

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Willem's implementation.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

wpythag2(N) ->
L = [{A, A*A} || A<- lists:seq(1,N)],
lists:flatten([forAllBs(A, A2, L, N) || {A, A2}<- L]).

forAllBs(A, A2, L, N) ->
[forAllCs(A, B, A + B, A2 + B2, L, N) || {B, B2}<- L, A + B< N].

forAllCs(A, B, AB, A2B2, L, N) ->
[{A, B, C} || {C, C2}<- L, A2B2 =:= C2, AB + C =< N].

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Hynek's new version
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

pythag4(N) when is_integer(N) -> pythag4(N,1).

pythag4(N, A) when A+2> N -> [];
pythag4(N, A) -> pythag4(N, A, A*A, 1).

pythag4(N, A, _A2, B) when A+B+1> N -> pythag4(N, A+1);
pythag4(N, A, A2, B) -> pythag4(N, A, A2, B, B*B, 1).

pythag4(N, A, A2, B, _B2, C) when A+B+C> N -> pythag4(N, A, A2, B+1);
pythag4(N, A, A2, B, B2, C) when A2 + B2 =:= C*C ->
[{A, B, C}|pythag4(N, A, A2, B, B2, C+1)];
pythag4(N, A, A2, B, B2, C) -> pythag4(N, A, A2, B, B2, C+1).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Willem's implementation in parallel by Hynek
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

wpythag2P(N) ->
L = [{A, A*A} || A<- lists:seq(1,N)], % For all A's
lists:flatten(lpmap(fun({A, A2}) -> % For all B's in parallel
[forAllCsWH(A, B, A + B, A2 + B2, L, N)
|| {B, B2}<- L, A + B< N]
end, L, 2000, ordered)).

forAllCsWH(A, B, AB, A2B2, L, N) ->
[{A, B, C} || {C, C2}<- L, A2B2 =:= C2, AB + C =< N].

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Morten's implementation.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

pythag5(N) when is_integer(N) ->
Primes = sieve(N div 2),
M1M2s = incorporate_primes([{1,1}], N, Primes),
lists:usort(lists:flatten([ [{A,B,C}, {B,A,C}] || {M1, M2}<- M1M2s,
M1> M2, A<- [M1-M2], B<- [2*round(math:sqrt(M1*M2))], C<- [M1+M2],
A+B+C =< N])).

sieve(N) when is_integer(N) ->
erase(),
sieve(N,2).

sieve(N, K) when K>= N ->
[X || X<- lists:seq(2, N), erase(X) == undefined];
sieve(N, K) ->
cross_off(K, K, N div K - 1),
sieve(N, find_next_in_sieve(K + 1)).

cross_off(_K, _Current, 0) ->
ok;
cross_off(K, Current, Left) ->
Next = Current + K,
put(Next, out),
cross_off(K, Next, Left - 1).

find_next_in_sieve(K) ->
case get(K) of
undefined ->
K;
_ ->
find_next_in_sieve(K+1)
end.

incorporate_prime(M1M2s, N, P) ->
lists:flatten([incorporate_prime_single({M1,M2}, N, P)|| {M1, M2}<-
M1M2s]).

incorporate_prime_single({M1,M2}, N, P) ->
Evens = [{X, Y} || X<- incorporate_prime_even(M1, N, P), Y<-
incorporate_prime_even(M2, N, P)],
Odds = [{X, Y} || X<- incorporate_prime_odd(M1, N, P), Y<-
incorporate_prime_odd(M2, N, P)],
Evens ++ Odds.

incorporate_prime_even(M, N, P) ->
incorporate_prime(M, N, P, []).

incorporate_prime_odd(M, N, P) ->
incorporate_prime(M * P, N, P, []).

incorporate_prime(M, N, _P, Acc) when M> N/2 ->
Acc;
incorporate_prime(M, N, P, Acc) ->
incorporate_prime(M * P * P, N, P, [M|Acc]).

incorporate_primes(M1M2s, _N, []) ->
M1M2s;
incorporate_primes(M1M2s, N, [P|Rest]) ->
M1M2s_new = incorporate_prime(M1M2s, N, P),
incorporate_primes(M1M2s_new, N, Rest).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Richard's improvement.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

py3R(N)->
[{A,B,C} ||
A <- lists:seq(1, N div 2),
B <- lists:seq(1, N - A),
C <- [trunc(math:sqrt(A * A + B * B))],
A + B + C =< N,
A*A + B*B =:= C*C].

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Joe's improvement.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
py3a(Max) ->
N = Max div 2,
[{A,B,C} ||
A <- lists:seq(1,N+1),
B <- lists:seq(1,Max-A),
C <- lists:seq(1,Max-A-B),
A*A + B*B =:= C*C].

%************************* END CODE*************************

Cheers,

Ivan.

=======================================================================
Este mensaje ha sido enviado mediante el servicio de correo electronico
que ofrece la Federacion de Radioaficionados de Cuba a sus miembros para
respaldar el cumplimiento de los objetivos de la organizacion y su
politica informativa. La persona que envia este correo asume el
compromiso de usar el servicio a tales fines y cumplir con las
regulaciones establecidas.

Edmond Begumisa

unread,

Nov 15, 2010, 8:53:30 PM11/15/10

to erlang-q...@erlang.org, Ivan Carmenates García, Richard O'Keefe, Håkan Huss \" \"

Below is a parallelised version of Richard/Håkan's blazing and yet
easy-to-read algo using a max of 2000 processes (find lpmap function on
previous thread.)

On my dual-core it's about 20% faster for when N is small (300) and gets
up to 40% faster when N is large (3,000). Evens out from there. I'm
getting about 4.5s for N = 5,000!

I must say, I'm floored. I need to learn much, much, more about Erlang. I
just wasn't convinced it would be any good at this! And their algo retains
Garcia's original beauty too! It reads much like the high-school math
problems.

py3R_p(N)->
lists:flatten(lpmap(fun(A) ->
[{A, B, C} ||

B <- lists:seq(1, N - A),
C <- [trunc(math:sqrt(A * A + B * B))],
A + B + C =< N,
A*A + B*B =:= C*C]

end, lists:seq(1, N div 2), 2000, ordered)).

- Edmond -

Håkan Huss" <hus...@gmail.com>

--
Using Opera's revolutionary e-mail client: http://www.opera.com/mail/

________________________________________________________________
erlang-questions (at) erlang.org mailing list.
See http://www.erlang.org/faq.html
To unsubscribe; mailto:erlang-questio...@erlang.org

Hynek Vychodil

unread,

Nov 20, 2010, 1:27:48 PM11/20/10

to Edmond Begumisa, erlang-q...@erlang.org

(I'm just retuning from Stockholm and sitting at the Vienna airport
and hacking for the fun, so ...)
This version should be little bit faster than py3R (~30% in bytecode
and ~20% in native for N=1000):

pythag4a(N) when is_integer(N) -> pythag4a(N,1).

pythag4a(N, A) when A+2 > N -> [];
pythag4a(N, A) when A rem 2 =:= 0 -> pythag4a(N, A, A*A, 1, 1);
pythag4a(N, A) -> pythag4a(N, A, A*A, 2, 2).

pythag4a(N, A, _A2, B, _S) when A+B+1 > N -> pythag4a(N, A+1);
pythag4a(N, A, A2, B, S) ->
B2 = B*B,
C = round(math:sqrt(A2+B2)),
if
A+B+C =< N, A2 + B2 =:= C*C ->
[{A, B, C}|pythag4a(N, A, A2, B+S, S)];
true ->
pythag4a(N, A, A2, B+S, S)
end.

--
--Hynek (Pichi) Vychodil

Analyze your data in minutes. Share your insights instantly. Thrill
your boss. Be a data hero!
Try GoodData now for free: www.gooddata.com

Reply all

Reply to author

Forward