/* 

  Miller-Rabin primality test (Rosetta code) in Picat.

  http://rosettacode.org/wiki/Miller-Rabin_primality_test
  """
  The Miller–Rabin primality test or Rabin–Miller primality test is a primality test: an algorithm 
  which determines whether a given number is prime or not.

  The algorithm, as modified by Michael O. Rabin to avoid the generalized Riemann hypothesis, 
  is a probabilistic algorithm.

  The pseudocode, from Wikipedia is:

  Input: n > 2, an odd integer to be tested for primality;
         k, a parameter that determines the accuracy of the test
  Output: composite if n is composite, otherwise probably prime
  write n − 1 as 2s·d with d odd by factoring powers of 2 from n − 1
  LOOP: repeat k times:
     pick a randomly in the range [2, n − 1]
     x ← ad mod n
     if x = 1 or x = n − 1 then do next LOOP
     for r = 1 .. s − 1
       x ← x2 mod n
       if x = 1 then return composite
       if x = n − 1 then do next LOOP
     return composite
  return probably prime

    The nature of the test involves big numbers, so the use of "big numbers" libraries (or similar 
    features of the language of your choice) are suggested, but not mandatory.
    Deterministic variants of the test exist and can be implemented as extra (not mandatory 
    to complete the task) 
  """

  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.


go =>
 
  println([N : N in 3..100, miller_rabin(N,5)]),
  nl.

go2 =>
  _ = random2(),
  % garbage_collect(200_000_000),
  N = 123456791234567891234567,
  println(n=N),
  if miller_rabin(N,20) then
    println(ok1)
  else 
    println(not_ok1)
  end,
  nl.


go3 =>
  _ = random2(),
  % garbage_collect(200_000_000),
  member(Prime,primes(100_000_000)),
  Prime > 2,
  println(prime=Prime),
  if miller_rabin(Prime,20) then
    println(ok1)
  else 
    println(not_ok1),
    fail
  end,
  fail,
  nl.

% From the C++ solution
miller_rabin(N,K) =>
  if N <= 2 ; N mod 2 = 0 then
    false
  end,
  S = N - 1,
  while (S mod 2 == 0)  
    S := S >> 1
    % S := S // 2
  end,
  Check = true,
  foreach(_I in 1..K, break(Check==false))
    A = random(1,N-1),
    Temp = S,
    Mod = pow_mod(A,Temp,N),
    while (Temp != N-1, Mod != 1, Mod != N-1)
      Mod := (Mod*Mod) mod N,
      Temp := Temp * 2
    end,
    if Mod != N-1, Temp mod 2 == 0 then
      Check := false
    end
  end,
  Check=true.