/***************************************************************************/
/* */
/* ACM ALGORITHM 406: Exact solution of Linear Equations */
/* using Residue Arithmetic */
/* */
/* Jo Ann Howell, The University of Texas at Austin, */
/* Center for Numerical Analysis, Austin, Tx 78712, U.S.A. */
/* 23 March 1970, 2 July 1970, amended 6 October 1971. */
/* */
/***************************************************************************/
/* Translated to Fortran 90 by R. A. Vowels, 21 November 2007. */
/* Translated then to PL/I by R. A. Vowels, 23 November 2007. */
/* Produces the same results as the description for the algorithm, */
/* for the 5 x 5 matrix with 5 right-hand sides: */
/* 5
5
5 300 -2100 4200 -2520 1 0 0 0 0
-60 -2400 18900 -40320 25200 0 1 0 0 0
210 6300 -52920 117600 -75600 0 0 1 0 0
-280 -6720 58800 -134400 88200 0 0 0 1 0
126 2520 -22680 52920 -35280 0 0 0 0 1
Also gives the correct results with the 10 x 10 example (see below).
Note that the values for B given are incorrect.
The correct values are B = [0 0 0 0 0 0 0 0 0 1].
The above obtained with FIXED BINARY (31) for TYPE (see below).
10
1
10 9 8 7 6 5 4 3 2 10000000 0
9 9 8 7 6 5 4 3 2 1 0
8 8 8 7 6 5 4 3 2 1 0
7 7 7 7 6 5 4 3 2 1 0
6 6 6 6 6 5 4 3 2 1 0
5 5 5 5 5 5 4 3 2 1 0
4 4 4 4 4 4 4 3 2 1 0
3 3 3 3 3 3 3 3 2 1 0
2 2 2 2 2 2 2 2 2 1 0
1 1 1 1 1 1 1 1 1 1 1
DETERMINANT= 1.00000000000000000E+0000
============= Multilength digits for DETERMINANT OF A =========
0 0001
=============== X ===============
-1.99999980000000000E+0007
1.99999980000000000E+0007
0.00000000000000000E+0000
0.00000000000000000E+0000
0.00000000000000000E+0000
0.00000000000000000E+0000
0.00000000000000000E+0000
0.00000000000000000E+0000
-1.00000000000000000E+0000
2.00000000000000000E+0000
=============== Multilength digits of Y ===============
-1999 -9998
1999 9998
0 0000
0 0000
0 0000
0 0000
0 0000
0 0000
0 -0001
0 0002
*/
/* Notes: */
/* Arguments redundant in PL/I have been deleted from the list of */
/* parameters in procedure EXACT. */
/* All float variables changed to double precision. */
/* When FIXED DECIMAL (31) is used in the next statement but one, */
/* it is necessary to change the ten prime numbers in KPRIME, */
/* and to change NDIGIT to 15. */
/* The following values of KPRIME with NDIGIT = 4 may be used if TYPE is */
/* FIXED BINARY (31): */
/* 46237, 46261, 46271, 46273, 46279, 46301, 46307, 46309, 46327, 46337 */
/* The following values of KPRIME with NDIGIT = 7 may be used if TYPE is */
/* FIXED DECIMAL (15): */
/* 10000019, 10000079, 10000103, 10000121, 10000139, */
/* 10000141, 10000169, 10000189, 10000223, 10000229 */
/* The following values of KPRIME with NDIGIT = 15 may be used if TYPE */
/* is FIXED DECIMAL (31): */
/* 55555399, 55555429, 55555459, 55555481, 55555517, */
/* 55555523, 55555531, 55555543, 55555547, 55555553 */
/* The following values of KPRIME with NDIGIT = 15 may be used if TYPE */
/* is FIXED DECIMAL (31): */
/* 3162277660000373, 3162277660000391, 3162277660000429, */
/* 3162277660000469, 3162277660000583, 3162277660000607, */
/* 3162277660000649, 3162277660000717, 3162277660000777, */
/* 3162277660000813 */
% DECLARE TYPE CHARACTER;
% TYPE = 'FIXED DECIMAL (31)';
/* This subroutine solves the matrix equation AX=B for X and for the exact */
/* solution, Y=A(adj)*B and DET A. Residue arithmetic is used to obtain */
/* the solution. */
EXACT: PROCEDURE (A, B, KPRIME, X, DET, IER,
MULTL, LCOUNT, ATEMP, RY, W, V) OPTIONS (REORDER);
DECLARE ( N, IN, M, IM, IMPIN, IMIN1, NDIGIT, NOPRIM, IER, LCOUNT )
FIXED BINARY (31);
DECLARE ( A(*,*), B(*,*), ATEMP(HBOUND(A,1), HBOUND(A,1)+HBOUND(B,2)),
MULTL(HBOUND(A,1)*HBOUND(B,2)+1,HBOUND(KPRIME,1)), MM(HBOUND(KPRIME,1)),
KPRIME(*), W(*), V(*) ) TYPE;
DECLARE ( X(HBOUND(A,1),HBOUND(B,2)), DET,
RY(HBOUND(A,1)*HBOUND(B,2)+1) ) FLOAT (18);
/* A is the N by N coefficient matrix and must be of type integer. */
/* N is the order of the matrix A (N greater than 1). */
/* IN is a dimension parameter which defines the dimensions of A. */
/* It must be equal to or greater than N. */
/* B is the N by M matrix of the right-hand side and must be */
/* of type integer. */
/* M is the number of columns of B and X (M greater than 0). */
/* IM is a dimension parameter which defines the second dimension */
/* of the 2-dimensional arrays B and X. It must be equal to or */
/* greater then M. */
/* IMPIN is a dimension parameter which is IM + IN. */
/* IMIN1 is a dimension parameter which is IM * IN + 1. */
/* NDIGIT is the number of digits stored in each word during */
/* multilength arithmetic operations. It is machine dependent */
/* and must be chosen so that 10 ** (2 * NDIGIT) is less than */
/* or equal to the largest representable integer for the */
/* computer being used. (This is now computed automatically.) */
/* KPRIME is the linear array of NOPRIM moduli. The moduli must be */
/* primes, chosen as large as possible and so that the square of */
/* the largest prime does not overflow an integer word. */
/* NOPRIM is a dimension parameter which denotes the number of */
/* primes (moduli) stored in KPRIME. */
/* NO2 is a dimension parameter which is 2*NOPRIM. */
/* (It is no longer required to be passed as an argument.) */
/* X is the M by N floating-point matrix with the M solution */
/* vectors as columns. It is the rounded quotient of the */
/* rational components of X. */
/* DET is the floating-point determinant of A. */
/* IER is an error code which is */
/* 0 if the system is solved satisfactorily, */
/* 1 if there are not enough moduli available to solve the */
/* system, */
/* 2 if the coefficient matrix is singular modulo each of the */
/* NOPRIM moduli (in which case X and DET are not computed), */
/* 3 if one or more of the input integer arguments is */
/* incorrect (i.e., N, M, IN, IM, IMPIN, IMIN1, NO2). */
/* MULTL is the matrix in which the multilength digits of Y(I,J) */
/* and DET A are stored. The elements of Y are stored by */
/* columns in the first M * N rows of MULTL, and DET A is */
/* stored in the (M * N + 1)th row. Low order digits are in */
/* column 1 of MULTL, and the highest order digits are in */
/* column LCOUNT. It should be dimensioned IMIN1 by NOPRIM. */
/* LCOUNT is the column number in MULTL which contains the highest */
/* order multilength digits. */
/* ATEMP is the IN by IMPIN of type integer used by EXACT to hold */
/* the augmented matrix (A,B) in residue form. */
/* MM is the linear array used by SOLVE to hold the moduli that */
/* were used to solve the system of equations. It should be */
/* dimensioned the same as KPRIME. */
/* RY is the linear array used by EXACT to store the floating- */
/* point elements of Y and the floating-point determinant of */
/* A. Its dimension should be IMIN1. */
/* W is the linear array of type integer used by MLTLTH to hold */
/* a multilength number while performing multilength */
/* ariithmetic operations on it. It shoudl be dimensioned NO2. */
/* V is the linear array of type integer used by CHECK for */
/* comparing the values of two multilength numbers. It should */
/* be dimensioned the same as W. */
/* Local variables. */
DECLARE ( SUMLOG, BOUND, P ) FLOAT (18);
DECLARE ( I, J, ICOUNT, INDEX, NO2 ) FIXED BINARY (31);
DECLARE ( NZ, IS, IFLAG, IQUIT, NORES ) EXTERNAL FIXED BINARY (31),
( IB, PP ) TYPE EXTERNAL;
ON SUBSCRIPTRANGE SNAP BEGIN;
PUT SKIP DATA (I, J);
END;
NDIGIT = NDIG();
PUT SKIP DATA (NDIGIT);
NOPRIM = HBOUND(KPRIME,1);
NO2 = 2*HBOUND(KPRIME,1);
IN, N = HBOUND(A,1); IM, M = HBOUND(B,2);
IMPIN = HBOUND(A,2) + HBOUND(B,2);
IMIN1 = M*N+1;
/* Check input parameters for consistency. */
IF N ^= HBOUND(A,2) THEN GO TO L80;
IF N ^= HBOUND(B,1) THEN GO TO L80;
IF N <= 1 THEN GO TO L80;
IF M <= 0 THEN GO TO L80;
IF IMPIN ^= IM+IN THEN GO TO L80;
IF IMIN1 ^= IM*IN+1 THEN GO TO L80;
NORES = M*N+1;
IB = 10**NDIGIT;
SUMLOG = 0;
NZ = 0; /* NZ is the number of primes for which */
/* the residue system is singular. */
IQUIT = 0; /* If IQUIT is not equal to 0, then A is */
/* singular modulo each of the stored primes. */
IS = 1; /* IS will count the number of primes used */
/* successfully. */
ICOUNT = 1; /* ICOUNT will count the number of primes tried. */
/* Compute a lower bound on the number of required moduli. */
CALL LOGBND (A, N, IN, B, M, IM, BOUND);
/* Compute the residue of A and B and store both in ATEMP. */
L10:
PP = KPRIME(ICOUNT);
P = PP;
DO I = 1 TO N;
DO J = 1 TO N;
ATEMP(I,J) = REM(A(I,J), PP);
END;
END;
DO I = 1 TO N;
DO J = 1 TO M;
ATEMP(I, N+J) = REM(B(I,J), PP);
END;
END;
IFLAG = 0;
/* Solve the residue system AX = B (mod PP) for */
/* Y = A(adj)*B (mod PP) and DET (mod PP) and store results in MULTL. */
CALL SOLVE (ATEMP, MULTL, N, IN, MM, M, IMPIN, IMIN1, NOPRIM);
/* If IQUIT is not equal to 0, then the system is singular modulo */
/* each of the stored primes, and hence, cannot be solved by this */
/* program. */
IF IQUIT ^= 0 THEN
DO;
IER = 2;
RETURN;
END;
/* If IFLAG is not equal to 0, then A is singular modulo */
/* KPRIME(ICOUNT). Choose another prime, i.e. KPRIME(ICOUNT+1), */
/* and try to solve the system again. */
IF IFLAG ^= 0 THEN GO TO L50;
SUMLOG = SUMLOG + LOG(P);
/* Test to see whether the required number of primes has been used. */
IF SUMLOG >= BOUND THEN GO TO L60;
IS = IS + 1;
L50:
ICOUNT = ICOUNT + 1;
/* If all primes have been tried and still another is required, */
/* combine results and check solution. */
IF ICOUNT <= NOPRIM THEN GO TO L10;
IS = IS - 1;
/* Combime results by conversion to symmetric mixed-radix. */
L60:
CALL MIXED_RADIX (MULTL, MM, RY, LCOUNT, NDIGIT, IMIN1, NOPRIM, NO2, W);
/* Check the solution by computing AY and DB. */
CALL CHECK (A, N, IN, B, M, IM, IER, MULTL, IMIN1, NOPRIM, NO2, W, V);
IF IER = 1 THEN
RETURN;
/* Compute the solution X = (1/DET) * Y. */
DET = RY(NORES);
IF DET = 0 THEN
DO;
IER = 2; RETURN;
END;
INDEX = 0;
DO J = 1 TO M;
DO I = 1 TO N;
INDEX = INDEX + 1;
X(I,J) = RY(INDEX)/DET;
END;
END;
RETURN;
/* Return error code of 3 for inconsistent input parameters. */
L80:
IER = 3;
RETURN;
/* This function returns the largest value NDIGIT such that */
/* 10**(2*NDIGIT) does not overflow an integer. */
(FIXEDOVERFLOW, SIZE):
NDIG: PROCEDURE RETURNS (FIXED BINARY);
DECLARE NDIGIT FIXED BINARY;
DECLARE K TYPE;
ON FIXEDOVERFLOW GO TO FOUND;
ON SIZE GO TO FOUND;
DO NDIGIT = 1 TO 100;
K = 10.00000000E0**(2*NDIGIT);
END;
FOUND:
RETURN (NDIGIT-1);
END NDIG;
END EXACT;
/* This subroutine solves the residue system */
/* AX = B (mod PP) for Y (mod PP) and DET (mod PP). */
SOLVE: PROCEDURE (ATEMP, MULTL, N, IN, MM, M, IMPIN, IMIN1, NOPRIM) OPTIONS (REORDER);
DECLARE ( ATEMP(*,*), MM(*), MULTL(*,*) ) TYPE;
DECLARE ( N, IN, M, IMPIN, IMIN1, NOPRIM ) FIXED BINARY (31);
DECLARE ( NZ, IS, IFLAG, IQUIT, NORES ) EXTERNAL FIXED BINARY (31),
( IB, PP ) TYPE EXTERNAL;
DECLARE ( MPN, L, I, J, JJ, INDEX ) FIXED BINARY (31);
DECLARE ( IDET, ITEMP, IK, IX) TYPE;
ON SUBSCRIPTRANGE SNAP BEGIN;
PUT SKIP DATA (I, J);
END;
IDET = 1;
/* Find a pivotal element relatively prime to PP. */
MPN = M + N;
DO J = 1 TO N;
DO I = J TO N;
IF REM(ATEMP(I,J), PP) ^= 0 THEN LEAVE;
IF I = N THEN GO TO L100;
END;
/* Permute rows I and J. */
IF I = J THEN GO TO L40;
IDET = -IDET;
DO JJ = J TO MPN;
ITEMP = ATEMP(J,JJ);
ATEMP(J,JJ) = ATEMP(I,JJ);
ATEMP(I,JJ) = ITEMP;
END;
/* Accumulate determinant. */
L40:
IDET = IDET * ATEMP(J,J);
IDET = REM(IDET, PP);
/* Find inverse of pivotal element. */
IX = INVERSE (ATEMP(J,J), PP);
/* Multiply row J by the inverse of pivotal element. */
DO JJ = J TO MPN;
ITEMP = ATEMP(J,JJ)*IX;
ATEMP(J,JJ) = REM(ITEMP, PP);
END;
/* Replace Lth row by Lth row - Jth row, (L not equal J). */
DO L = 1 TO N;
IF L = J THEN
ITERATE;
IK = ATEMP(L,J);
DO JJ = J TO MPN;
ITEMP = ATEMP(J,JJ)*IK;
ITEMP = REM(ITEMP, PP);
ITEMP = ATEMP(L,JJ) - ITEMP;
ATEMP(L,JJ) = REM(ITEMP, PP);
END;
END;
END;
/* Store symmetric residue digits in MULTL, and modulus in MM. */
INDEX = 0;
DO J = N+1 TO MPN;
DO I = 1 TO N;
INDEX = INDEX + 1;
ITEMP = ATEMP(I,J) * IDET;
MULTL(INDEX,IS) = RESIDUE(ITEMP, PP);
END;
END;
MULTL(NORES,IS) = RESIDUE(IDET, PP);
MM(IS) = PP;
RETURN;
L100:
NZ = NZ + 1;
IFLAG = 1;
/* Test to see if all primes have failed. */
IF NZ > NOPRIM-1 THEN
IQUIT = 1;
END SOLVE;
/* Subroutine MIXED_RADIX computes the symmetric mixed-radix of Y and */
/* DET A from their symmetric residue digits. */
MIXED_RADIX: PROCEDURE (MULTL, MM, RY, LCOUNT, NDIGIT, IMIN1, NOPRIM, NO2, W) OPTIONS (REORDER);
DECLARE ( MULTL(*,*), MM(*), W(*) ) TYPE,
( LCOUNT, NDIGIT, IMIN1, NOPRIM, NO2 ) FIXED BINARY (31);
DECLARE ( RY(*) ) FLOAT (18);
DECLARE ( ACC, ACC1, ACC2, TEX ) FLOAT (18);
DECLARE ( NZ, IS, IFLAG, IQUIT, NORES ) EXTERNAL FIXED BINARY (31),
( IB, PP ) TYPE EXTERNAL;
DECLARE ( I, J, K, KK ) FIXED BINARY (31);
DECLARE (ITEMP, IX) TYPE;
ON SUBSCRIPTRANGE SNAP BEGIN;
PUT SKIP DATA (I, J, K, KK);
END;
ON ERROR PUT SKIP DATA (I, J, K, KK);
IF IS = 1 THEN
DO;
DO I = 1 TO NORES;
RY(I) = MULTL(I,I);
END;
LCOUNT = 1;
RETURN;
END;
/* Compute symmetric mixed-radix digits & store them in MULTL. */
DO I = 2 TO IS;
KK = I - 1;
DO J = I TO IS;
IX = INVERSE(MM(KK), MM(J));
DO K = 1 TO NORES;
ITEMP = MULTL(K,J) - MULTL(K,I-1);
ITEMP = ITEMP * IX;
MULTL(K,J) = RESIDUE(ITEMP,MM(J));
END;
END;
END;
/* Compute Y and D from their symmetric mixed-radix digits */
/* using multilength arithmetic. */
LCOUNT = 0;
DO I = 1 TO NORES;
W(1) = 1 /* ! IS THIS A 1 OR I ? Looks like 1, see 5 lines down. */;
DO K = 2 TO IS;
W(K) = 0;
END;
/* Compute Y(I) = (...(MULTL(I,IS)*MM(IS-1) + MULTL(I,IS-1))*MM(IS-2) + */
/* ... + MULTL(I,2))*MM(1) + MULTL(I,1). */
W(1) = W(1)*MULTL(I,IS) * MM(IS-1); /* ! Try I in place of 1 in W(I). */
CALL MLTLTH (W);
J = IS;
DO FOREVER;
J = J - 1;
IF J <= 1 THEN
LEAVE;
W(1) = W(1) + MULTL(I,J); /* <<<<<<<<< I OR 1 ? */
CALL MLTLTH (W);
DO K = 1 TO IS;
W(K) = W(K) * MM(J-1);
END;
CALL MLTLTH (W);
END;
W(1) = W(1) + MULTL(I,J) /* Is this I or 1 ? */;
CALL MLTLTH (W);
/* Store multilength digits of Y(I) in MULTL(I,J), J = 1,IS. */
/* Store multilength digits of DET A in MULTL(NORES,J), J=1,IS. */
DO J = 1 TO IS;
MULTL(I,J) = W(J);
END;
/* Compute Y(I) in floating-point from multilength digits. */
K = IS;
L80:
IF W(K) ^= 0 THEN GO TO L90;
IF K = 1 THEN GO TO L100;
K = K - 1;
GO TO L80;
L100:
RY(I) = W(1); /* With W(1), program finishes and gives correct */
/* results for 10 x 10 example. */
GO TO L130;
L90:
IF K <= 1 THEN GO TO L100;
/* ACC = FLOAT(W(K)) * IB + W(K-1); /* This was necessary during */
/* testing, but not now. */
ACC = W(K) * IB + W(K-1);
TEX = NDIGIT * (K-2);
IF K <= 2 THEN GO TO L120;
ACC1 = W(K-2);
ACC2 = 10.0000000000000000E0**NDIGIT;
IF ACC2 = 0 THEN
DO;
PUT SKIP LIST ( '*** ERROR, DIVISION BY ZERO');
ACC = 1;
END;
ELSE
ACC = ACC + ACC1/ACC2;
L120:
ACC1 = 1.00000000000000000E+1**TEX;
RY(I) = ACC*ACC1;
L130:
IF K > LCOUNT THEN LCOUNT = K;
END;
RETURN;
END MIXED_RADIX;
/* MLTLTH is called after arithmetic is performed on any of the digits */
/* in W. This procedure adjusts the content of W so that each element */
/* contains NDIGIT digits. */
MLTLTH: PROCEDURE (W) OPTIONS (REORDER);
DECLARE W(*) TYPE;
DECLARE ( K, L ) FIXED BINARY (31), IW TYPE;
DECLARE ( NZ, IS, IFLAG, IQUIT, NORES ) EXTERNAL FIXED BINARY (31),
( IB, PP) TYPE EXTERNAL;
ON SUBSCRIPTRANGE SNAP BEGIN;
PUT SKIP DATA (K);
END;
IF IS = 1 THEN RETURN;
L = IS - 1;
/* Distribute the digits in W so that each element of W contains */
/* NDIGIT digits. */
DO K = 1 TO IS-1;
IW = W(K)/IB;
W(K+1) = IW + W(K+1);
W(K) = -IW * IB + W(K);
END;
K = IS;
/* All the elements of W should have the same sign. */
L20:
IF W(K) < 0 THEN GO TO L60;
IF W(K) > 0 THEN GO TO L40;
IF K = 1 THEN
RETURN;
K = K - 1;
GO TO L20;
L40:
DO K = 1 TO L;
IF W(K) >= 0 THEN
ITERATE;
W(K) = W(K) + IB;
W(K+1) = W(K+1) - 1;
END;
RETURN;
L60:
DO K = 1 TO L;
IF W(K) <= 0 THEN
ITERATE;
W(K) = W(K) - IB;
W(K+1) = W(K+1) + 1;
END;
END MLTLTH;
/* Subroutine CHECK checks the solution by computing A*Y and (DET A)*B */
/* and comparing the results. */
/* Y is stored by columns in MULTL. */
/* DET is stored in MULTL(NORES,I), I = 1, IS). */
CHECK: PROCEDURE (A, N, IN, B, M, IM, IER, MULTL, IMIN1, NOPRIM, NO2, W, V) OPTIONS (REORDER);
DECLARE ( A(IN, IN), B(IN,IM), MULTL(IMIN1,NOPRIM), W(NO2), V(NO2) ) TYPE;
DECLARE ( N, IN, M, IM, IER, IMIN1, NOPRIM, NO2 ) FIXED BINARY (31);
DECLARE ( I, J, K, L, KK, LL, INDEX ) FIXED BINARY (31);
DECLARE ( II, JJ) TYPE;
DECLARE ( NZ, IS, IFLAG, IQUIT, NORES ) EXTERNAL FIXED BINARY (31),
( IB, PP) TYPE EXTERNAL;
ON SUBSCRIPTRANGE SNAP BEGIN;
PUT SKIP DATA (I, J);
END;
LL = IS;
KK = IS + 1;
DO I = 1 TO N;
INDEX = 0;
DO L = 1 TO M;
/* Mutiply row I of A by column L of Y. */
W = 0;
IS = KK;
DO J = 1 TO N;
INDEX = INDEX + 1;
JJ = A(I,J) / IB;
II = -JJ*IB + A(I,J);
DO K = 2 TO LL;
W(K) = W(K) + MULTL(INDEX,K)*II + MULTL(INDEX,K-1)*JJ;
CALL MLTLTH (W);
END;
W(1) = II*MULTL(INDEX,1) + W(1);
W(KK) = JJ*MULTL(INDEX,LL) + W(KK);
CALL MLTLTH (W);
IF ABS(W(IS)) < IB THEN
ITERATE;
IF IS >= NO2 THEN
ITERATE;
IS = IS + 1;
W(IS) = W(IS-1) / IB;
W(IS-1) = -W(IS)*IB + W(IS-1);
END;
/* Store the product in V. */
DO K = 1 TO IS;
V(K) = W(K);
END;
/* Multiply B(I,L) by DET and store in W. */
JJ = B(I,L) / IB;
II = -JJ*IB + B(I,L);
DO K = 2 TO LL;
W(K) = MULTL(NORES,K)*II + MULTL(NORES,K-1)*JJ;
END;
W(1) = II*MULTL(NORES,1);
W(KK) = JJ*MULTL(NORES,LL) + W(KK);
CALL MLTLTH (W);
/* Test equality of W and V. */
DO J = 1 TO IS;
IF W(J) ^= V(J) THEN GO TO L100;
END;
END;
END;
/* If solution checks, return IER = 0, else return IER = 1. */
IER = 0;
IS = LL;
RETURN;
L100:
IER = 1;
END CHECK;
/* INVERSE computes an inverse of K (mod M) by the Euclidian algorithm. */
INVERSE: PROCEDURE (K, M) RETURNS ( TYPE );
DECLARE ( K, M ) TYPE;
DECLARE ( I, J, KK, L, NN, Result ) TYPE;
I = K; L = M;
J = 1;
Result = 0;
L10:
KK = I/L;
NN = REM(I, L);
IF NN = 0 THEN GO TO L20;
I = L;
L = NN;
NN = -KK*Result + J;
J = Result;
Result = NN;
GO TO L10;
L20:
IF L >= 0 THEN GO TO L30;
Result = -Result;
/* Return a positive value. */
L30:
IF Result >= 0 THEN
RETURN (Result);
RETURN ( M + Result );
END INVERSE;
/* The function RESIDUE computes the symmetric residue of K (mod M) */
/* i.e. -M/2 less than RESIDUE less than M/2. */
RESIDUE: PROCEDURE (K, M) RETURNS ( TYPE );
DECLARE ( K, M ) TYPE;
DECLARE Result TYPE;
Result = REM(K, M);
IF Result < 0 THEN
DO;
IF 2*Result + M >= 0 THEN
RETURN (Result);
Result = Result + M;
END;
ELSE IF Result = 0 THEN
RETURN (Result);
ELSE
DO;
IF -2*Result+M >= 0 THEN
RETURN (Result);
Result = Result - M;
END;
RETURN (Result);
END RESIDUE;
/* BOUND is a lower bound for the log of the product of the moduli. */
LOGBND: PROCEDURE (A, N, IN, B, M, IM, BOUND) OPTIONS (REORDER);
DECLARE ( A(*,*), B(*,*) ) TYPE,
( N, IN, M, IM ) FIXED BINARY (31);
DECLARE ( BOUND ) FLOAT (18);
DECLARE ( ALPHA, TEMP ) FLOAT (18);
DECLARE ( I, J ) FIXED BINARY (31);
ON SUBSCRIPTRANGE SNAP BEGIN;
PUT SKIP DATA (I, J);
END;
/* BOUND is a lower bound for the log of the product of the moduli. */
BOUND = 0;
DO I = 1 TO N;
ALPHA = 0;
DO J = 1 TO N;
TEMP = A(I,J);
TEMP = TEMP*TEMP;
ALPHA = ALPHA + TEMP;
END;
BOUND = BOUND + LOG(ALPHA);
END;
BOUND = BOUND/2;
DO J = 1 TO M;
DO I = 1 TO N;
ALPHA = ABS(B(I,J));
IF ALPHA = 0 THEN
ITERATE;
BOUND = BOUND + LOG(ALPHA);
END;
END;
BOUND = BOUND + LOG(2.00000000000000000E0);
END LOGBND;