/* 

  Prime decomposition in Picat v3.

  http://rosettacode.org/wiki/Prime_decomposition

  This is a port of the Prolog solution.


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

*/

main => go.


go ?=>
  time(prime_decomp(9007199254740991, L1)),
  println(L1),

  time(prime_decomp(576460752303423487, L2)),
  println(L2),

  time(prime_decomp(1361129467683753853853498429727072845823, L3)),
  
  println(L3),
  nl.
  
go => true.


go2 =>
  foreach(P in primes(60))
     PP = 2**P-1,
     println([p=P,pp=PP]),
     time(prime_decomp(PP,L)),
     println(L),
     nl
  end,  
  nl.

go2 => true.

prime_decomp(N, L) :-
	SN is sqrt(N),
	prime_decomp_1(N, SN, 2, [], L).
 
 
prime_decomp_1(1, _, _, L, L) :- !.
 
% Special case for 2, increment 1
prime_decomp_1(N, SN, D, L, LF) :-
	(   0 is N mod D ->
	    Q is N // D,
	    SQ is ceiling(sqrt(Q)),
	    prime_decomp_1(Q, SQ, D, [D |L], LF)
	;
	    D1 is D+1,
	    (	D1 > SN ->
	        LF = [N |L]
	    ;
	        prime_decomp_2(N, SN, D1, L, LF)
	    )
	).
 
% General case, increment 2
prime_decomp_2(1, _, _, L, L) :- !.
 
prime_decomp_2(N, SN, D, L, LF) :-
	(   0 is N mod D ->
	    Q is N // D,
	    SQ is ceiling(sqrt(Q)),
	    prime_decomp_2(Q, SQ, D, [D |L], LF);
	    D1 is D+2,
	    (	D1 > SN ->
	        LF = [N |L]
	    ;
	        prime_decomp_2(N, SN, D1, L, LF)
	    )
	).