/*****************************************************************************/
/* */
/* ACM ALGORITHM 288. SOLUTION OF SIMULTANEOUS LINEAR DIOPHANTINE EQUATIONS */
/* */
/* W. A. Blankinship, National Security Agency, Ft. Geo. G. Meade, Md. */
/* Received 19 May 1965 and 17 September 1965. */
/* */
/*****************************************************************************/
/* Translated to PL/I by R. A. Vowels. 4 December 2006. */
/* A is (1:m,1:n) and contains the coefficients of the equations. */
/* b is (1:m) and contains the right-hand sides of the equations. */
/* x is (1:n) and contains the solutions. */
/* Y is (1:n,1:n) and contains the linearly independent auxiliary solutions. */
/* k is the number of linearly independent solutions stored in Y. */
/* m is the number of equations. */
/* n is the number of unknowns. */
/* A Diophantine solution is obtained only if SOLVE_INTEGER returns '1'B and */
/* the value of d returned is 1. */
SOLVE_INTEGER:
PROCEDURE (A, x, b, d, m, n, k, Y) RETURNS (BIT(1)) OPTIONS (REORDER);
declare (A(*,*), x(*), b(*)) fixed decimal (prec),
(d, m, n, k) fixed binary, Y(*,*) fixed decimal (prec);
declare mat(n+1, m+n+1) fixed decimal (prec);
declare (i, j, rank, s) fixed binary;
declare (np1, mp1) fixed binary;
/* Make -b the first row of mat. */
/* The second and subsequent rows of mat are the transpose of A. */
do j = 1 to m;
mat(1,j) = -b(j);
do i = 1 to n; mat(i+1,j) = A(j,i); end;
end;
/* Thus, each equation is stored in one column, with the RHS */
/* at the top of the column. */
/* You read the equations from the bottom of the column to the */
/* top. */
/* Generate a unit matrix in the columns to the right of */
/* (b and the tranpose of A). */
do j = 1 to n+1;
do i = 1 to n+1;
if i = j then mat(i,j+m) = 1; else mat(i,j+m) = 0;
end;
end;
np1 = n + 1; mp1 = m + 1; /* for DR PL/I. */
rank = INTRANK (mat, np1, mp1, n);
d = mat(rank, m+binary(1));
if d = 0 then return ('0'B);
do i = rank to m;
if mat(rank, i) ^= 0 then return ('0'B);
end;
if d < 0 then s = -1; else s = 1;
d = s * d;
/* above line can be replaced by d = abs(d); */
k = n - rank + 1;
do i = 1 to n;
x(i) = mat (rank, m+i+1) * s;
do j = 1 to k;
Y(j,i) = mat(rank+j, m+i+1);
end;
end;
return ('1'B);
END SOLVE_INTEGER;