/* 

  Primality by trial division (Rosetta code) in Picat.

  http://rosettacode.org/wiki/Primality_by_trial_division
  """
  Write a boolean function that tells whether a given integer is prime. Remember that 
  1 and all non-positive numbers are not prime.

  Use trial division. Even numbers over two may be eliminated right away. A loop from 3 
  to sqrt(n) will suffice, but other loops are allowed. 
  """
  

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

*/

import util.

main => go.


% imperative: 0.983s
go =>
  println([I : I in 1..100, is_prime1(I)]),
  nl,
  foreach(P in 1..6) 
     Primes = [ I : I in 1..10**P, is_prime1(I)],
     println([10**P,Primes.len])
  end,
  nl.

% recursive: is_prime2/1: 3.22s
go2 =>
  foreach(P in 1..6) 
     Primes = [ I : I in 1..10**P, is_prime2(I)],
     println([10**P,Primes.len])
  end,
  nl.

% functional: is_prime3/1: 0.978s 
go3 =>
  foreach(P in 1..6) 
     Primes = [ I : I in 1..10**P, is_prime3(I)],
     println([10**P,Primes.len])
  end,
  nl.

% Using prim2/1: 2.115s
go4 =>
  foreach(P in 1..6) 
     Primes = [ I : I in 1..10**P, prime2(I)],
     println([10**P,Primes.len])
  end,
  nl.


% iterative
is_prime1(N) => 
  if N == 2 then
    true
  elseif N <= 1 ; N mod 2 == 0 then
    false
  else
    foreach(I in 3..2..ceiling(sqrt(N)))
       N mod I > 0
    end
  end.


% recursive
is_prime2(N) => 
  (N == 2 ; is_prime2b(N,3)).

is_prime2b(N,_Div), N <= 1 => false.
is_prime2b(2,_Div) => true.
is_prime2b(N,_Div), N mod 2 == 0 => false.
is_prime2b(N,Div), Div > ceiling(sqrt(N)) => true.
is_prime2b(N,Div), Div > 2 => 
 (N mod Div == 0 -> 
    false 
  ; 
    is_prime2b(N,Div+2)
  ).

% Functional
is_prime3(2) => true.
is_prime3(3) => true.
is_prime3(P) => P > 3, P mod 2 =\= 0, not has_factor3(P,3).  

has_factor3(N,L), N mod L == 0 => true.
has_factor3(N,L) => L * L < N, L2 = L + 2, has_factor3(N,L2).


% {{trans|Prolog}}
% prime2(N) can be used to generate primes until memory is exhausted
% Difference from Prolog implementation: Picat does not support between/3 with
% inf as upper bound, so we have to place a high number
prime2(2).
prime2(N) :-
  between(3, 2**156+1, N),
  N mod 2 = 1,              % odd
  M is floor(sqrt(N+1)),    % round-off paranoia 
  Max is (M-1) // 2,        % integer division
  foreach(I in 1..Max) N mod (2*I+1) > 0 end.