/* 

  Prime numbers and divisors (constraint model) in Picat.

  - is_prime(N,IsPrime): IsPrime #= 1 if N is a prime.

  - divisors(N, NumDivisors, Divisor): calculate the divisors of N
    This does not scale well.

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

*/

import cp.

main => go.


% Just prime checking. This scales much better.
go =>
  M = 1_000_000,
  N :: 2..M,

  IsPrime :: 0..1,
  is_prime(N,IsPrime),
  % IsPrime #= 1,

  % N #= 123_456,

  solve([N,IsPrime]),

  println(N=IsPrime),
  PrimeQ = cond(prime(N),1,0),
  if IsPrime != PrimeQ then
    println([wrong,N,is_prime=IsPrime,prime_q=PrimeQ])
  end,
  fail,
  nl.

%
% This does not scales well.
%
go2 =>
  M = 1200,
  NN :: 1..M,

  NumDivisors :: 2..1+(M div 2),
  % Divisor[I] = 1: I is a divisor of NN, 0 if not
  Divisor = new_list(M),
  Divisor :: 0..1,

  % is a prime
  IsPrime :: 0..1,

  % NN = 100,
  divisors(NN, NumDivisors, Divisor),
   
  IsPrime #= 1 #<=> NumDivisors #= 2,
  % IsPrime #= 0 #<=> NumDivisors #> 2,

  Vars = Divisor ++ [NN,NumDivisors,IsPrime],
  solve($[degree,updown],Vars),

  println(nn=NN),
  println(num_divisors=NumDivisors),
  % println(divisor=Divisor),  
  println(divisor=[I : I in 1..M, Divisor[I] == 1]),
  println(is_prime=NN=IsPrime),
  nl,
  fail,

  nl.


%
% Calculate the number of divisors of n
%
divisors(N, NumDivisors, D) =>
  Len = D.len,
  UbN = fd_max(N),
  T = (UbN div 2)+1,

  % is a divisor?
  foreach(I in 2..Len)
    D[I] #= 1 #<=> N mod I #= 0
  end,
  foreach(I in T+1..D.len-1)
   I #!= N #=> D[I] #= 0
  end,
  D[1] #= 1,
  % D[N] #= 1,
  element(N,D,1),
  NumDivisors #= sum(D).


is_prime(N, IsPrime) =>
  UbN = fd_max(N),
  T = ceiling(sqrt(UbN+1)),
  % standard brute force approach
  ( (N #!= 2 #/\ N mod 2 #= 0) #\/
    sum([ I #!= N #/\ N mod I #= 0 : I in 3..2..T]) #> 0)  #<=> IsPrime #= 0.