/* 

  Prime decomposition (Rosetta code) in Picat.
 
  http://rosettacode.org/wiki/Prime_decomposition
  """

  """

  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.

% Checking 2**prime-1
go =>
  foreach(P in primes(60))
     Factors = factors(2**P-1),
     println([n=2**P-1,factors=Factors])
  end,
  println(factors(1361129467683753853853498429727072845823)),
  nl.

go2 =>
   foreach(P in 2..200000)
      % println(P=factors(P))
      _ = factors(P)
   end,
   nl.

go3 =>
   println(factors(1361129467683753853853498429727072845823)),
   nl.

%
% factors of N
%
factors(N) = Factors =>
   Factors = [],
   M = N,    
   while (M mod 2 == 0) 
      Factors := Factors ++ [2],  
      M := M div 2 
   end,
   T = 3,
   while (M > 1, T < 1+(sqrt(M)))
      if M mod T == 0 then
         [Divisors, NewM] = alldivisorsM(M, T),
         Factors := Factors ++ Divisors,
         M := NewM
      end,
      T := T + 2
   end,
   if M > 1 then Factors := Factors ++ [M] end.

alldivisorsM(N,Div) = [Divisors,M] =>
   M = N,
   Divisors = [],
   while (M mod Div == 0) 
      Divisors := Divisors ++ [Div],
      M := M div Div
   end.


table
prime_cached(N) => prime(N).