[erlang-questions] Erlang shows its slow face!
Morten Krogh
mk@REDACTED
Sun Nov 14 17:26:32 CET 2010
Hi
Here is another version of the same program using a different approach.
if A^2 + B^2 = C^2, then A and B cannot both be odd, so at least one of the must be even, say, B.
Write B^2 = (C+A)*(C-A)
Now A and C are both even or both odd, so we can write, (B/2)^2 = M1 * M2, where M1 = (C+A)/2, M2 = (C-A)/2
where M1 and M2 are integers. Let p be a prime number and let alpha1 be the highest number such that p^alpha1 divides M1, and similarly for alpha2,
then alpha1 and alpha2 must both be even or both be odd, since the left hand side (B/2)^2 only has even powers of any prime factor.
On the other hand any M1, M2 built in this way will lead to a pythagorean triple, as can be seen by just plugging in.
If A+B+C <= N, then no prime number greater than N/2 can occur.
This lead to the following algorithm.
1. Generate all prime numbers up to N/2. I used Eratosthenes sieve for that.
2. Generate all M1 and M2, by distributing even or odd powers as explained above. Only keep numbers that do not already break A+B+C <= N, and M1 > M2.
3. Put them all in a list, permute A and B, e.g. {3,4,5} and {4,3,5}.
4. unique sort it. Unique because of the A and B permutation, and sort to make the outpout the same as all the other implementations.
The code is pasted below.
The runtime for N=300 is 5 ms on a computer where pythag3 is 80 ms.
But it is for larger values of N that this algorithm will really pull away.
For N=3000, it uses 680 ms compared to 76 s for pythag3.
The earlier programs were all N^3, this one seems, by testing, to be N^2. I haven't thought theoretically about the complexity of this algotithm.
I haven't optimised this program.
Cheers,
Morten.
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).
More information about the erlang-questions
mailing list