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