/* 

  Abundant, deficient and perfect number classifications (Rosetta code) in Picat.

  https://rosettacode.org/wiki/Abundant,_deficient_and_perfect_number_classifications
  """
  These define three classifications of positive integers based on their proper 
  divisors.

  Let P(n) be the sum of the proper divisors of n, where the proper divisors of n are 
  all positive divisors of   n   other than   n   itself.

    - if P(n) < n  then n is classed as deficient   (OEIS A005100).
    - if P(n) == n then n is classed as perfect   (OEIS A000396).
    - if P(n) > n  then n is classed as abundant   (OEIS A005101).

  Example

  6   has proper divisors of   1,   2,   and   3.

  1 + 2 + 3 = 6,   so   6   is classed as a perfect number.


  Task

  Calculate how many of the integers 1 to 20,000 (inclusive) are in each of the three classes.

  Show the results here. 
  """

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

*/


main => go.

go =>
  
  Classes = new_map([deficient=0,perfect=0,abundant=0]),
  foreach(N in 1..20000)
    C = classify(N),
    Classes.put(C,Classes.get(C)+1)
  end,
  % println(Classes.to_list().sort(1)),
  println(Classes),
  nl.

go2 =>
  
  Classes = new_map([deficient=0,perfect=0,abundant=0]),
  foreach(N in 1..20000)
    C = classify2(N,sum_divisors(N)),
    Classes.put(C,Classes.get(C)+1)
  end,
  % println(Classes.to_list().sort(1)),
  println(Classes),
  nl.


classify(N) = Class =>
 S = sum_divisors(N),
 if S < N then 
   Class1 = deficient
 elseif S = N then 
   Class1 = perfect
 elseif S > N then
   Class1 = abundant
 end,
 Class = Class1.

classify2(N,S) = C, S <  N => C = deficient.
classify2(N,S) = C, S == N => C = perfect.
classify2(N,S) = C, S >  N => C = abundant.

sum_divisors(N) = Sum =>
  sum_divisors(2,N,cond(N>1,1,0),Sum).

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

% Part 1: 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).

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


%
% some alternative - and slower - variant of sum_divisors
%

sum_divisors2(N) = Sum =>
  D = floor(sqrt(N)),
  Sum = cond(N>1,1,0),
  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 = cond(N>1,1,0),
  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) = cond(N>1,1,0)+sum([I + cond(I != N div I,N div I, 0) 
                          : I in 2..floor(sqrt(N)), N mod I == 0]).