/* 

  Amicable pairs (Rosetta code) in Picat.

  http://rosettacode.org/wiki/Amicable_pairs
  """
  Two integers N and M are said to be amicable pairs if N != M and the sum of the proper 
  divisors of N (sum(propDivs(N))) = M as well as sum(propDivs(M)) = N.

  For example 1184 and 1210 are an amicable pair (with proper divisors 
  1, 2, 4, 8, 16, 32, 37, 74, 148, 296, 592 and 
  1, 2, 5, 10, 11, 22, 55, 110, 121, 242, 605 respectively).

  Task

  Calculate and show here the Amicable pairs below 20,000; (there are eight). 
  """

  This Picat model was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Picat page: http://www.hakank.org/picat/

*/


% import util.
% import cp.


main => go.


go =>
  N = 20000-1,
  garbage_collect(200_000_000),

  println(amicable1),
  time(amicable1(N)),

  % initialize_table is needed to clear the table cache
  initialize_table,
  println(amicable2),
  time(amicable2(N)),

  initialize_table,
  println(amicable3),
  time(amicable3(N)),

  initialize_table,
  println(amicable4),
  time(amicable4(N)),

  initialize_table,
  println(amicable5),
  time(amicable5(N)),


  nl.


%
% cf http://hakank.org/picat/euler21.pi
% 0.236s
%
amicable1(N) => 
  Pairs = new_map(),
  foreach(A in 1..N) 
     B = sum_divisors(A),
     C = sum_divisors(B),
     if A != B, A == C then
        Pairs.put([A,B].sort(),1)
     end
  end,
  println(Pairs.keys().sort()),

  nl.


%
% as list comprehension
% 0.224s
%
amicable2(N) => 
  S = [[A,B].sort() : A in 1..N, 
         B = sum_divisors(A), 
         C = sum_divisors(B),
         A != B, A == C].remove_dups(),
  println(S).

%
% while loop
% 0.224s
%
amicable3(N) => 
  A = 1,
  while(A <= N)
     B = sum_divisors(A),
     if A < B, A == sum_divisors(B) then
       println([A,B])
     end,
     A := A + 1
  end,
  nl.

%
% foreach loop, everything in the condition
% 0.224s
%
amicable4(N) => 
  foreach(A in 1..N, B = sum_divisors(A), A < B, A == sum_divisors(B))
     println([A,B])
  end,

  nl.

%
% another list comprehension variant
% 0.22s
%
amicable5(N) => 
  println([[A,B] : A in 1..N, B = sum_divisors(A), A < B, A == sum_divisors(B)]).




%
% This recursive version is slightly faster than sum_divisors2/1.
%
table
sum_divisors(N) = Sum =>
  sum_divisors(2,N,1,Sum).

sum_divisors(I,N,Sum0,Sum), I > floor(sqrt(N)) =>
  Sum = Sum0.

% I is a divisor of N
sum_divisors(I,N,Sum0,Sum), N mod I == 0 =>
  Sum1 = Sum0 + I,
  (I != N div I -> 
    Sum2 = Sum1 + N div I 
    ; 
    Sum2 = Sum1
  ),
  sum_divisors(I+1,N,Sum2,Sum).

% I is no divisor of N.
sum_divisors(I,N,Sum0,Sum) =>
  sum_divisors(I+1,N,Sum0,Sum).



table
sum_divisors2(N) = Sum =>
    D = floor(sqrt(N)),
    Sum1 = 1,
    foreach(I in 2..D, N mod I == 0) 
        Sum1 := Sum1+I,
        if I != N div I then
            Sum1 := Sum1 + N div I
        end
    end,
    Sum = Sum1.