/* 

  Coprimes in Picat.

  https://www.mathsisfun.com/definitions/relatively-prime.html
  """
  Relatively primes

  When two numbers have no common factors other than 1.

  In other words there is no value that you could divide them both by exactly (without any remainder).

  21 and 22 are relatively prime:
  * The factors of 21 are 1, 3, 7 and 21
  * The factors of 22 are 1, 2, 11 and 22
  (the only common factor is 1)

  But 21 and 24 are NOT relatively prime:
  * The factors of 21 are 1, 3, 7 and 21
  * The factors of 24 are 1, 2, 3, 4, 6, 8, 12 and 24
  (the common factors are 1 AND 3)

  Relatively Prime is also called "coprime" or "mutually prime".
  """

  Probability Fact on Twitter
  https://twitter.com/ProbFact/status/1198354294942175233
  """
  The probability that two random integers are relatively prime is 6/π^2.
  """
  
  Mathematica:
  In[1]:= 1/Zeta[2]                                                                                         

            6
  Out[1]= ---
             2
           Pi


  This is simulated in go2/0.


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

*/

import ordset.

main => go.


go =>
  check_coprime(21,22),
  check_coprime(21,24),
  check_coprime(2,10),  
  nl.


%
% The probability that two random integers are coprime:
% 
% n = 10000
% numCoprimes = 6090
% numCoprimes = 0.609
% theoretical = 0.607927101854027
%
go2 =>
  % _ = random2(),
  N = 20000,
  println(n=N),
  NumCoprimes = 0,
  foreach(_ in 1..N)
    A = 2+random() mod N,
    B = 2+random() mod N,
    if coprime2(A,B) then
       NumCoprimes := NumCoprimes + 1
    end
  end,
  println(numCoprimes=NumCoprimes),
  println(numCoprimes=(NumCoprimes/N)),  
  println(theoretical=(6/math.pi**2)),
  nl.


%
% Check the three different implementations of coprime.
% Seems fine.
% 
go3 =>
  N = 1000,
  foreach(I in 1..N, J in 1..N)
    C1 = cond(coprime(I,J), true, false),
    C2 = cond(coprime2(I,J), true, false),
    C3 = cond(coprime3(I,J), true, false),
    if C1 != C2 ; C1 != C3 ; C2 != C3 then
      println([I,J, C1, C2, C3])
    end
  end,
  nl.

%
% Check the two divisors functions:
%
% divisors/1: CPU time 2.078 seconds.
% divisors/2: CPU time 1.085 seconds.
%
go4 =>
   N = 10000,
   time(test_divisors1(N)),
   time(test_divisors2(N)),
   nl.
     


table
divisors(N) = [ I : I in 1..(N div 2)+1, N mod I == 0].

% Faster version
table
divisors2(N) = [I : I in 2..2, N mod 2 == 0] ++ [ I : I in 3..2..(N div 2)+1, N mod I == 0].


test_divisors1(N) =>
  foreach(I in 1..N)
    _ = divisors(I)
  end.

test_divisors2(N) =>
  foreach(I in 1..N)
    _ = divisors2(I)
  end.


check_coprime(A,B) => 
  if coprime3(A,B) then
    printf("%d and %d are coprime\n", A, B)
  else
    printf("%d and %d are not coprime\n", A, B)  
  end.


% Slower
coprime(A,B) =>
  DivisorsA = divisors(A),
  DivisorsB = divisors(B),
  intersection(DivisorsA,DivisorsB) == [1].
 
% Faster than coprime/2
coprime2(A,B) =>
  DivisorsA = divisors2(A),
  DivisorsB = divisors2(B),
  intersection(DivisorsA,DivisorsB) == [].


% Slightly faster than coprime2/2
coprime3(A,B) =>
  N = 2,
  if A mod N == 0, B mod N == 0 then
    false
  end,
  N := 3,
  Check = true,
  while (N <= min(A,B) div 2, Check==true)
    if A mod N == 0, B mod N == 0 then
       Check := false
    end,
    N := N+2
  end,
  Check == true.