/* 

  GCD and LCM in Picat.

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

*/

import cp.


main => go.

go ?=>
  N = 100,

  A :: 1..N,
  B :: 1..N,
  GCD :: 1..N,
  LCM :: 1..N*N,

  % GCD #= 1,
  LCM #= 60,

  gcd(A,B,GCD),
  lcm(A,B,LCM),  

  A #< B,

  Vars = [A,B,GCD,LCM],
  % solve($[ff,split,max(LCM)],Vars), % quite slow
  % solve($[ff,split,max(GCD)],Vars), % quite slow
  solve($[ff,split],Vars),  
  
  println([A,B,gcd=GCD,lcm=LCM]),
  fail,
   
  nl.

go => true.


go2 ?=>
  N = 1000,

  A :: 1..N,
  B :: 1..N,
  GCD :: 1..N,
  
  % gcd(A,B,GCD),
  % gcd2(A,B,GCD),
  euclid(A,B,GCD),
  
  GCD #= 8,
  % A #!= B,

  solve($[ff],[GCD,A,B]),
  % gcd(A,B,GCD),
  println([a=A,b=B,gcd=GCD]),
  fail,

  nl.
go2 => true.

%
% Use Picat's non-CP gcd/2 _after_ solve.
%
go3 ?=>
  % First solve the "main" problem
  N = 100000,
  [A,B] :: 1..N,

  A #< B,
  A * 100 + B * 10 #> 1000,
  (A * B) mod 8 #= 1,
  
  solve([A,B]),
  
  % Then use Picat's non-CP gcd/2 _after_ solve.
  GCD = gcd(A,B),
  println([a=A,b=B,gcd=GCD]),
  fail,
  nl.

go3 => true.


% From
% https://stackoverflow.com/questions/22450582/program-for-finding-gcd-in-prolog
% ported to CP
gcd2(X, 0, GCD):- X #= GCD.
gcd2(0, X, GCD):- X #= GCD.
gcd2(X, Y, GCD):- X #=< Y, Z #= Y - X, gcd2(X, Z, GCD).
gcd2(X, Y, GCD):- gcd2(Y, X, GCD).

% From
% https://stackoverflow.com/questions/66016948/inverting-a-function-via-clpfd-in-pure-prolog
euclid(A,A,A).
euclid(A,B,R) :- A #< B, C #= B-A, euclid(A,C,R).
euclid(A,B,R) :- A #> B, C #= A-B, euclid(C,B,R).

euclid(A,A,1,A).
euclid(A,B,P,R) :- A #< B, C #= B-A, P #= 2*Q, Q #> 0, euclid(A,C,Q,R).
euclid(A,B,P,R) :- A #> B, C #= A-B, P #= 2*Q+1, Q #> 0, euclid(C,B,Q,R).

%
% gcd(X,Y,G)
% G is the gcd of X and Y
% hakank: This is my decomposition.
% Not very efficient for larger domains.
%
gcd(X, Y, GCD) =>
  fd_min_max(X,_MinX,MaxX),
  fd_min_max(Y,_MinY,MaxY),
  PMax = min(MaxX,MaxY),
  P :: 1..PMax,

  foreach (I in 1..PMax)
    I #= P #=> (X mod I #= 0 #/\ Y mod I #= 0 #/\ GCD #= I),
    foreach(J in I+1..PMax) 
      I#>= P #=> (X mod J #!= 0 #\/ Y mod J #!= 0)
    end
  end.


%
% lcm(X,Y,LCM)
% LCM is the lcm of x and y
%
% Not very efficient for larger domains.
%
lcm(X, Y, LCM) =>
  fd_min_max(X,_MinX,MaxX),
  fd_min_max(Y,_MinY,MaxY),
  PMax = min(MaxX,MaxY),
  GCD :: 1..PMax,
  gcd(X,Y,GCD),
  LCM #= X*Y div GCD.