[erlang-questions] Erlang shows its slow face!

[Morten Krogh]

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).