/************************************************************************/ /* */ /* ACM ALGORITHM 290 - LINEAR EQUATIONS, EXACT SOLUTIONS */ /* */ /* J. Boothroyd, University of Tasmania, Hobart, Tasmania, Australia. */ /* */ /* 7 September 1965, revised 21 March 1966. */ /* */ /************************************************************************/ /* Translated to PL/I by R. A. Vowels, 24 November 2006. */ /* Solves the matrix equations Ax = b for A = a(1:n, 1:n) */ /* and x, b(1:n) where the elements of A, b are small integers */ /* and the results are required as ratios of integers. The solution */ /* vector overwrites b and has values given by det A * x where det A */ /* is the determinant of A and x is the true solution vector. */ /* [That is, x = b / det A] */ /* The user is warned that this procedure, of limited though useful */ /* application, is not a substitiue for other well-established */ /* methods of solving general sets of linear equations owing to the */ /* inherent danger of integer overflow. This may occur in the */ /* reduction if the elements of the matrix are large or in the back[ward]*/ /* substitution if the determinant and/or the elements of the right- */ /* hand side are large and may even occur with small elements and */ /* determinant if the order of the matrix and the nature of the */ /* equations combine to produce large solution values. Four */ /* devices intended to avoid integer overflow are incorporated. */ /* These are, (1) choice of column pivots having the smallest */ /* non-zero absolute value, (2) division by previous pivots (both after */ /* Fox, L., "An Introduction to Numerical Linear Algebra", Oxford */ /* U. Press, New York, 1965, p. 82), and (3) the local procedures */ /* CROSSMPY and ABDIVC which respectively evaluate integer expressions */ /* of the form (a*b-c*d)/e and a*b/c by performing */ /* the divisions before the multiplications. The output parameter */ /* det yields the determinant of A. If A is singular det = 0. */ /* Translator's note: */ /* Boothroyd's algorithm used binary integers for the data and all */ /* intermediate working. */ /* This version uses decimal arithmetic to PREC digits, and consequently*/ /* the risk of overflow is somewhat reduced compared to Boothroyd's */ /* algorithm. Use the preprocessor to set the number of digits to */ /* any value between 15 and 31. The risk of overflow should be almost */ /* non-existant. */ (fixedoverflow, size): EXACTLE: PROCEDURE(a, b, det) OPTIONS (REORDER); declare (a(*,*), b(*)) fixed decimal (prec), det fixed decimal (prec); declare (piv, pivot, sum, arii, aki) fixed decimal (prec), (i, j, k, m, n, pivi, ri, rk) fixed binary; declare r(hbound(a,1)) fixed binary, zpiv bit (1) aligned; declare false bit (1) static initial ('0'b), true bit (1) static initial ('1'b); n = hbound(a,1); ON FIXEDOVERFLOW SNAP BEGIN; det = 0; go to out; END; ON SIZE SNAP BEGIN; det = 0; go to out; END; crossmpy: procedure (a, b, c, d, e) returns (fixed decimal (prec)); declare (a, b, c, d, e) fixed decimal (prec); declare (qab, qcd, r, res) fixed decimal (prec); if abs(a) > abs(b) then do; qab = divide(a, e, prec); r = a - qab * e; qab = qab * b; res = r * b; end; else do; qab = divide (b, e, prec); r = b - qab * e; qab = qab * a; res = r * a; end; if abs(c) > abs(d) then do; qcd = divide(c, e, prec); r = c - qcd * e; qcd = qcd * d; res = res - r * d; end; else do; qcd = divide (d, e, prec); r = d - qcd * e; qcd = qcd * c; res = res - r * c; end; return ( qab - qcd + divide (res, e, prec)); end crossmpy; abdivc: procedure (a, b, c, sum) returns (fixed decimal (prec)); declare (a, b, c, sum) fixed decimal (prec); declare (q, r, temp) fixed decimal (prec); if abs(a) > abs(b) then do; q = divide (a, c, prec); temp = q * b; r = a - c * q; q = divide(b, c, prec, 0); sum = sum + (b - q*c) * r; end; else do; q = divide (b, c, prec); temp = q * a; r = b - c * q; q = divide (a, c, prec); sum = sum + (a - q*c) * r; end; return (temp + q * r); end abdivc; permb: procedure (b, r, n); declare b(*) fixed decimal (prec), (r(*), n) fixed binary; declare (i, k) fixed binary, w fixed decimal (prec); do i = n to 2 by -1; k = r(i); L: if k ^= i then do; if k > i then do; k = r(k); go to L; end; w = b(i); b(i) = b(k); b(k) = w; end; end; end permb; m = 1; do i = 1 to n; r(i) = i; end; do i = 1 to n; pivot = 0; zpiv = true; do k = i to n; aki = abs(a(r(k), i)); if zpiv & aki > 0 | aki ^= 0 & aki < abs(pivot) then do; zpiv = false; pivi = k; pivot = a(r(k), i); end; end; if pivot = 0 then do; det = 0; RETURN; end; ri = r(pivi); r(pivi) = r(i); r(i) = ri; if pivi ^= i then m = -m; do k = i + 1 to n; rk = r(k); aki = a(rk, i); do j = i+1 to n; if i = 1 then a(rk, j) = a(rk, j) * pivot - aki * a(ri, j); else a(rk, j) = crossmpy (a(rk, j), pivot, aki, a(ri, j), piv); end; if i = 1 then b(rk) = b(rk) * pivot - aki * b(ri); else b(rk) = crossmpy (b(rk), pivot, aki, b(ri), piv); end; piv = pivot; end; ri = r(n); if m ^= 1 then do; det = -a(ri, n); aki = det; b(ri) = -b(ri); end; else do; det = a(ri, n); aki = det; end; do i = n-1 to 1 by -1; ri = r(i); arii = a(ri,i); sum = 0; piv = abdivc (b(ri), aki, arii, sum); sum = -sum; do j = i+1 to n; piv = piv - abdivc (b(r(j)), a(ri,j), arii, sum); end; b(ri) = piv - divide (sum, arii, prec); end; call permb (b, r, n); out: END EXACTLE;