Solution of a Diophantine Equation

% 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 is abs(A), 
  AbsB is 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):-!,
  fwrite(s, 0, 0, `Indeterminate`), nl, fail.
solution_exists(0, _):-!,
  fwrite(s, 0, 0, `No solution`), nl, fail.
solution_exists(GCD, C):-0 =\= C mod GCD, !,
  fwrite(s, 0, 0, `No solution`), nl, fail.
solution_exists(_, _).

LPA Index     Home Page

Valid HTML 4.01 Strict