/* 

  Euler #23 in Picat.

  """
  A perfect number is a number for which the sum of its proper divisors 
  is exactly equal to the number. For example, the sum of the proper divisors 
  of 28 would be 1 + 2 + 4 + 7 + 14 = 28, which means that 28 is a perfect number.

  A number n is called deficient if the sum of its proper divisors is less than 
  n and it is called abundant if this sum exceeds n.

  As 12 is the smallest abundant number, 1 + 2 + 3 + 4 + 6 = 16, the smallest number 
  that can be written as the sum of two abundant numbers is 24. By mathematical 
  analysis, it can be shown that all integers greater than 28123 can be written 
  as the sum of two abundant numbers. However, this upper limit cannot be reduced 
  any further by analysis even though it is known that the greatest number that 
  cannot be expressed as the sum of two abundant numbers is less than this limit.

  Find the sum of all the positive integers which cannot be written as the sum of 
  two abundant numbers.
  """ 

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

*/

main => time(go).

go => euler23.

% inspired by a C++ solution
% 1.29s
euler23 =>
    Limit = 20161,
    Arr = new_array(Limit), bind_vars(Arr,1), 
    foreach(I in 2..Limit,J in I*2..I..Limit)
       Arr[J] := Arr[J] + I 
    end,

    Abundant = [I: I in 12..Limit, Arr[I] > I],

    /*
    foreach(A in Abundant, B in Abundant,A <= B, A+B < Limit)
       Arr[A + B] := 0
    end,
    */

    % Faster
    foreach(A in Abundant)
       Found = 0,
       foreach(B in Abundant, A <= B, Found == 0)
         if A+B < Limit then
            Arr[A + B] := 0
         else 
            Found := 1
         end
       end
    end,

    println(sum([I : I in 1..Limit, Arr[I] != 0])).


% 1.374s
euler23a1 =>
  % N = 28123, 
  % From http://mathworld.wolfram.com/AbundantNumber.html: 
  %  "Every number greater than 20161 can be expressed as a sum of two abundant numbers."
  N = 20161, 
  Abundant = [A : A in 1.. N, sum_divisors5(A) > A],
  % Vec = [0 : _I in 1..N].to_array(),
  Vec = new_array(N), bind_vars(Vec,0),
  foreach(A in Abundant, B in Abundant, A >= B, A+B <= N) 
    Vec[A+B] := 1
  end,
  println(sum([A : A in 1..N, Vec[A] == 0])).


% 1.38s
euler23a =>
  N = 20161, 
  Abundant = [A : A in 1.. N, sum_divisors2(A) > A],
  Vec = new_array(N), bind_vars(Vec,0),
  % foreach(I in 1..N) Vec[I] := 0 end,
  foreach(A in Abundant, B in Abundant) 
    if A >= B, A+B <= N then
      Vec[A+B] := 1
    end
  end,
  println(sum([A : A in 1..N, Vec[A] == 0])).


% 1.38s
euler23b =>
  N = 20161, 
  Abundant = [A : A in 1.. N, sum_divisors4(A) > A],
  Vec = new_array(N), bind_vars(Vec,0),
  % foreach(I in 1..N) Vec[I] := 0 end,
  foreach(A in Abundant, B in Abundant, A >= B, A+B <= N) 
    Vec[A+B] := 1
  end,
  println(sum([A : A in 1..N, Vec[A] == 0])).


% using map: 2.3s
euler23c =>
  N = 20161, 
  Abundant = [A : A in 1.. N, sum_divisors5(A) > A],
  Vec = new_map(),
  foreach(A in Abundant, B in Abundant, A >= B, A+B <= N) 
    Vec.put(A+B,1)
  end,
  println(sum([A : A in 1..N, not Vec.has_key(A)])).

% testing Vec var() instead: 1.396s
euler23d =>
  N = 20161, 
  Abundant = [A : A in 1.. N, sum_divisors2(A) > A],
  Vec = new_array(N), % we don't initialize this with 0
  foreach(A in Abundant, B in Abundant) 
    if A >= B, A+B <= N then
      Vec[A+B] := 1
    end
  end,
  println(sum([A : A in 1..N, var(Vec[A])])).


% 1.44s, slightly different approach (a litle slower)
euler23e =>
  N = 20161, 
  Abundant = [A : A in 1.. N, sum_divisors5(A) > A],
  Vec = (1..N).to_array(),
  foreach(A in Abundant, B in Abundant, A+B <= N, A >= B) 
    Vec[A+B] := 0
  end,
  println(sum([Vec[I] : I in 1..N])).

% using only list is very slow
euler23f =>
  N = 20161, 
  Abundant = [A : A in 1.. N, sum_divisors5(A) > A],
  Vec = (1..N),
  foreach(A in Abundant, B in Abundant, A+B <= N, A >= B) 
    Vec[A+B] := 0
  end,
  println(sum(Vec)).

% out of memory
euler23g =>
  N = 20161,
  Abundant = [A : A in 1.. N, sum_divisors5(A) > A],
  Vec = (1..N),
  foreach(A in Abundant, B in Abundant, A+B <= N, A >= B) 
    Vec := delete(Vec,A+B)
  end,
  println(sum(Vec)).



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


sum_divisors3(N) = Sum =>
  Sum = 1,
  foreach(I in 2..floor(sqrt(N)), N mod I == 0) 
     Sum := Sum + I + cond(I != N div I,N div I, 0)
  end.

% as a sum
sum_divisors4(N) = 1+sum([I + cond(I != N div I,N div I, 0) 
                          : I in 2..floor(sqrt(N)), N mod I == 0]). 


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

% Part 0: base case
sum_divisors5(I,N,Sum0,Sum), I > floor(sqrt(N)) =>
  Sum = Sum0.

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

% Part 2: I is no divisor of N.
sum_divisors5(I,N,Sum0,Sum) =>
  sum_divisors5(I+1,N,Sum0,Sum).


% very slow
is_abundant(J) =>
  Divisors = 0,
  Loop = 1,
  foreach(N in 2..J, N*N <= J, Loop == 1)
    if J mod N == 0 then
       Divisors := Divisors + N,
       if N < J div N then
          Divisors := Divisors + J div N
       end,
       if J <= Divisors then
         Loop := 0
       end
    end
  end,
  J <= Divisors.