/* 

  Diophantic equations in Picat v3.

  The Prolog code is from Colin Barker's collection of Prolog programs.
  http://colin.barker.pagesperso-orange.fr/lpa/dioph.htm

  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 ?=>
  L1 = findall([X1,Y1],dioph(98, -199, 5, X1, Y1)),
  println(l1=L1),
  nl,

  % From https://mathworld.wolfram.com/DiophantineEquation.html
  L2 = findall([X2,Y2],dioph(1027, 712, 1, X2, Y2)),
  println(l2=L2),
  nl,

  % https://brilliant.org/wiki/linear-diophantine-equations-one-equation/
  L3 = findall([X3,Y3],dioph(141, 34, 30, X3, Y3)),
  println(l3=L3),

  % This has no solution
  L4 = findall([X4,Y4],dioph(2, 2, 3, X4, Y4)),
  println(l4=L4),

  
  nl.

go => true.


% From Colin Barker's collection of Prolog programs.
% http://colin.barker.pagesperso-orange.fr/lpa/dioph.htm

% """
% dioph(+A, +B, +C, ?X0, ?Y0) is true if x = X0 + B * n and
%   y = Y0 - A * n (for any integer n) is the general solution
%   of the Diophantine equation Ax + By = C.
% e.g. dioph(98, -199, 5, -136, -67).
% """
dioph(A, B, C, X, Y):-
  AbsA = abs(A),
  AbsB = abs(B),
  gcd_extended(AbsA, AbsB, X0, Y0, GCD),
  solution_exists(GCD, C),
  A0 is A // GCD,
  B0 is B // GCD,
  C0 is C // GCD,
  X1 is C0 * X0 * sign(A),
  Y1 is C0 * Y0 * sign(B),
  C1 is X1 // B0,
  X is X1 - C1 * B0,
  Y is Y1 + C1 * A0.

% """
% gcd_extended(+X, +Y, ?U1, ?U2, ?Z) is true if X*U1 + Y*U2 = Z = gcd(X, Y).
%   e.g. gcd_extended(24140, 16762, -234, 337, 34).
% """
gcd_extended(X, Y, U1, U2, Z):-
  gcd_extended_1(Y, 1, 0,  X, 0, 1,  Z, U2, U1).

gcd_extended_1(0, _, _,  U3, U2, U1,  U3, U2, U1):-!.
gcd_extended_1(V3, V2, V1,  U3, U2, U1,  U3P, U2P, U1P):-
  Q  is U3 // V3,
  W3 is U3 mod V3,
  W2 is U2 - V2 * Q,
  W1 is U1 - V1 * Q,
  gcd_extended_1(W3, W2, W1,  V3, V2, V1,  U3P, U2P, U1P).

solution_exists(0, 0):-!,
  write('Indeterminate'), nl, fail.
solution_exists(0, _):-!,
  write('No solution'), nl, fail.
solution_exists(GCD, C):-0 =\= C mod GCD, !,
  write('No solution'), nl, fail.
solution_exists(_, _).