/*

  Perfect numbers in Picat.

  See http://rosettacode.org/wiki/Perfect_numbers
  """
  A number is perfect if the sum of its factors is equal to 
  twice the number. An equivalent condition is that n is perfect 
  if the sum of n's factors that are less than n is equal to n. 
  """

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

*/

main => go.

% Simple version:
go =>
   println("perfect1:/1"),
   println(perfect1=[I : I in 1..10_000, perfect1(I)]),
   nl.
   

%
% Using the formula: 2^(p-1)*(2^p-1) for prime p
% 
% 6 (prime 2)
% 28 (prime 3)
% 496 (prime 5)
% 8128 (prime 7)
% 33550336 (prime 13)
% 8589869056 (prime 17)
% 137438691328 (prime 19)
% 2305843008139952128 (prime 31)
% 
% CPU time 118.039 seconds. Backtracks: 0
%
% Here we use the faster perfect2/1. 
go2 =>
   % Check perfect number candidates generated by the
   % formula 2^(p-1)*(2^p-1)
   println("Using the formula: 2^(p-1)*(2^p-1) for prime p"),
   foreach(P in primes(32))
      PF=perfectf(P),
      % Check that it is really a perfect number
      if perfect2(PF) then
        printf("%w (prime %w)\n",PF,P)
      end
   end,
   nl.

% Using the list of the first 48 primes that generates perfect numbers
% according to the formula 2^(p-1)*(2^p-1)
% from https://en.wikipedia.org/wiki/Perfect_number
%
% Print the length of the number if >= 80 digits
%
% prime 2: 6
% prime 3: 28
% prime 5: 496
% prime 7: 8128
% prime 13: 33550336
% prime 17: 8589869056
% prime 19: 137438691328
% prime 31: 2305843008139952128
% prime 61: 2658455991569831744654692615953842176
% prime 89: 191561942608236107294793378084303638130997321548169216
% prime 107: 13164036458569648337239753460458722910223472318386943117783728128
% prime 127: 14474011154664524427946373126085988481573677491474835889066354349131199152128
% prime 521: len = 314
% prime 607: len = 366
% prime 1279: len = 770
% prime 2203: len = 1327
% prime 2281: len = 1373
% prime 3217: len = 1937
% prime 4253: len = 2561
% prime 4423: len = 2663
% prime 9689: len = 5834
% prime 9941: len = 5985
% prime 11213: len = 6751
% prime 19937: len = 12003
% prime 21701: len = 13066
% prime 23209: len = 13973
% prime 44497: len = 26790
% prime 86243: len = 51924
% prime 110503: len = 66530
% prime 132049: len = 79502
% prime 216091: len = 130100

% CPU time 64.956 seconds. Backtracks: 0

go3 =>
  ValidP = [2, 3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127, 521, 607,
            1279, 2203, 2281, 3217, 4253, 4423, 9689, 9941, 11213,
            19937, 21701, 23209, 44497, 86243, 110503, 132049, 216091,
            
            756839, 859433, 1257787, 1398269, 2976221, 3021377, 6972593,
            13466917, 20996011, 24036583, 25964951, 30402457, 32582657,
            37156667, 42643801, 43112609, 57885161],
  foreach(P in ValidP, P <= 216091 )
      printf("prime %w: ", P),
      PF = perfectf(P),
      Len = PF.to_string.len,
      % Len = int_len(PF),
      if Len < 80 then
        println(PF)
      else
         println(len=Len)
      end,
      if Len >= 100_000 then
        fail
      end
  end,
  nl.
  

% Formula for perfect number candidates:
%   2^(p-1)*(2^p-1) where p is a prime
% 
perfectf(P) = (2**(P-1))*((2**P)-1).

% Definitions of perfect

% Very slow
perfect1(N) => sum(divisors(N)) == N.

% Faster
perfect2(N) => sum_divisors(N) == N.

% Very slow
divisors(N) = [J: J in 1..1+N div 2, N mod J == 0].

% Sum of divisors: Faster
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 not a divisor of N.
sum_divisors(I,N,Sum0,Sum) =>
  sum_divisors(I+1,N,Sum0,Sum).


int_len1(V) = Len =>
  Len1=1,
  while (V > 9)
     Len1 := Len1 + 1,
     V := V div 10
  end,
  Len = Len1.

% log10/1 Doesn't work with large integers
int_len2(N) = floor(log10(N))+1.

int_len(N) = Digits => 
    if N <= 999999999999997 then
       Digits = int_len1(N)
    else
        Digits1 = 15,
        while (N >= 10**Digits1)
            Digits1 := Digits1 + 1
        end,
        Digits = Digits1
    end.