C     ALGORITHM 406 COLLECTED ALGORITHMS FROM ACM.
C     ALGORITHM APPEARED IN COMM. ACM, VOL. 14, NO. 03,
C     P. 180.
      SUBROUTINE EXACT(A,N,IN,B,M,IM,IMPIN,IMIN1,NDIGIT,
     1KPRIME,NOPRIM,NO2,X,DET,IER,MULTL,LCOUNT,ATEMP,MM,
     2RY,W,V)
      DIMENSION A(IN,IN),B(IN,IM),X(IN,IM),ATEMP(IN,IMPIN),
     1MULTL(IMIN1,NOPRIM),MM(NOPRIM),RY(IMIN1),
     2KPRIME(NOPRIM),W(NO2),V(NO2)
      INTEGER A,B,ATEMP,PP,W,V
      COMMON/MLEN/IB,PP,NZ,IS,IFLAG,IQUIT,NORES
C
C         THIS SUBROUTINE SOLVES THE MATRIX EQUATION AX=B
C         FOR X AND FOR THE EXACT SOLUTION, Y=A(ADJ)*B
C         AND DET A. RESIDUE ARITHMETIC IS USED TO OBTAIN
C         THE SOLUTION.
C
C         A  IS THE N BY N COEFFICIENT MATRIX AND MUST BE
C            OF TYPE INTEGER.
C         N  IS THE ORDER OF THE MATRIX A (N GREATER THAN 1).
C        IN  IS A DIMENSION PARAMETER WHICH DEFINES THE
C            DIMENSION OF A.  IT MUST BE EQUAL TO OR GREATER
C            THAN N.
C         B  IS THE N BY M MATRIX OF THE RIGHT-HAND SIDE AND
C            MUST BE OF TYPE INTEGER.
C         M  IS THE NUMBER OF COLUMNS OF B AND X (M GREATER
C            THAN 0).
C        IM  IS A DIMENSION PARAMETER WHICH DEFINES THE
C            SECOND DIMENSION OF THE 2-DIMENSIONAL ARRAYS
C            B AND X.  IT MUST BE EQUAL TO OR GREATER THAN M.
C     IMPIN  IS A DIMENSION PARAMETER WHICH IS IM + IN.
C     IMIN1  IS A DIMENSION PARAMETER WHICH IS IM * IN + 1.
C    NDIGIT  IS THE NUMBER OF DIGITS STORED IN EACH WORD
C            DURING MULTILENGTH ARITHMETIC OPERATIONS.  IT IS
C            MACHINE DEPENDENT AND MUST BE CHOSEN SO THAT
C            10 ** (2 * NDIGIT) IS LESS THAN OR EQUAL TO THE
C            LARGEST REPRESENTABLE INTEGER FOR THE COMPUTER
C            BEING USED.
C    KPRIME  IS THE LINEAR ARRAY OF NOPRIM MODULI.  THE
C            MODULI MUST BE PRIMES, CHOSED AS LARGE AS
C            POSSIBLE AND SO THAT KPRIME(I) * KPRIME(J) DOES
C            NOT OVERFLOW AN INTEGER WORD, FOR ALL I AND J.
C    NOPRIM  IS A DIMENSION PARAMETER WHICH DENOTES THE
C            NUMBER OF PRIMES (MODULI) STORED IN KPRIME.
C       NO2  IS A DIMENSION PARAMETER WHICH IS 2*NOPRIM.
C         X  IS THE N BY M FLOATING-POINT MATRIX WITH THE
C            M SOLUTION VECTORS AS COLUMNS.  IT IS THE
C            ROUNDED QUOTIENT OF THE RATIONAL COMPONENTS
C            OF X.
C       DET  IS THE FLOATING POINT DETERMINANT OF A.
C       IER  IS AN ERROR CODE WHICH IS
C              0 IF THE SYSTEM IS SOLVED SATISFACTORILY,
C              1 IF THERE ARE NOT ENOUGH MODULI AVAILABLE
C                TO SOLVE THE SYSTEM,
C              2 IF THE COEFFICIENT MATRIX IS SINGULAR
C                MOLULO EACH OF THE NOPRIM MODULI (IN WHICH
C                CASE X AND DET ARE NOT COMPUTED),
C              3 IF ONE OR MORE OF THE INPUT INTEGER
C                ARGUMENTS IS INCORRECT (I.E. N,M,IN,IM,
C                IMPIN,IMIN1,NO2).
C     MULTL  IS THE MATRIX IN WHICH THE MULTILENGTH DIGITS
C            OF Y(I,J) AND DET A ARE STORED.  THE ELEMENTS
C            OF Y ARE STORED BY COLUMNS IN THE FIRST M * N
C            ROWS OF MULTL, AND DET A IS STORED IN THE
C            (M * N + 1)TH ROW.  LOW ORDER DIGITS ARE IN
C            COLUMN ONE OF MULTL, AND HIGHEST ORDER DIGITS
C            ARE IN COLUMN LCOUNT.  IT SHOULD BE DIMENSIONED
C            IMIN1 BY NOPRIM.
C    LCOUNT  IS THE COLUMN NUMBER IN MULTL WHICH CONTAINS
C            THE HIGHEST ORDER MULTILENGTH DIGITS.
C     ATEMP  IS THE IN BY IMPIN MATRIX OF TYPE INTEGER USED
C            BY EXACT TO HOLD THE AUGMENTED MATRIX (A,B)
C            IN RESIDUE FORM.
C        MM  IS THE LINEAR ARRAY USED BY SOLVE TO HOLD THE
C            MODULI WHICH WERE USED TO SOLVE THE SYSTEM OF
C            EQUATIONS.  IT SHOULD BE DIMENSIONED THE SAME
C            AS KPRIME.
C        RY  IS THE LINEAR ARRAY USED BY EXACT TO STORE THE
C            FLOATING-POINT ELEMENTS OF Y AND THE FLOATING-
C            POINT DETERMINANT OF A.  ITS DIMENSION SHOULD
C            BE IMIN1.
C         W  IS THE LINEAR ARRAY OF TYPE INTEGER USED RY
C            MLTLTH TO HOLD A MULTILENGTH NUMBER WHILE
C            PERFORMING MULTILENGTH ARITHMETIC OPERATIONS
C            ON IT.  IT SHOULD BE DIMENSIONED NO2.
C         V  IS THE LINEAR ARRAY OF TYPE INTEGER USED BY
C            CHECK FOR COMPARING THE VALUES OF TWO MULTILENGTH
C            NUMBERS.  IT SHOULD BE DIMENSIONED THE SAME AS W.
C         SAME AS W
C
C      CHECK INPUT PARAMETERS FOR CONSISTENCY
       IF(N .LE.1 .OR. N .GT. IN)  GO TO 80
       IF(M .LE.0 .OR. M .GT. IM)    GO TO 80
        IF(IMPIN .NE. (IM+IN))  GO TO 80
        IF(IMIN1 .NE. (IM*IN+1))  GO TO 80
        IF(NO2 .NE. (2*NOPRIM))  GO TO 80
       NORES =M*N+1
      IB=10**NDIGIT
      SUMLOG=0.
C NZ IS THE NUMBER OF PRIMES FOR WHICH
C THE RESIDUE SYSTEM IS SINGULAR
      NZ=0
C IF IQUIT IS NOT EQUAL TO 0, THEN A IS SINGULAR
C MODULO EACH OF THE STORED PRIMES
      IQUIT=0
C IS WILL COUNT THE NUMBER OF PRIMES
C USED SUCCESSFULLY
      IS=1
C ICOUNT WILL COUNT THE NUMBER OF PRIMES TRIED
      ICOUNT=1
C COMPUTE A LOWER BOUND ON THE NUMBER
C OF REQUIRED MODULI
      CALL LOGBND(A,N,IN,B,M,IM,BOUND)
C COMPUTE RESIDUE OF A AND B
C AND STORE BOTH IN ATEMP
   10 PP=KPRIME(ICOUNT)
      P=PP
      DO 20 I=1,N
   20     ATEMP(I,J)=MOD(A(I,J),PP)
C    THERE WAS A UNNEEDED CARD HERE
      DO 30 I=1,N
        DO 30 J=1,M
          JJ=N+J
   30     ATEMP(I,JJ)=MOD(B(I,J),PP)
      IFLAG=0
C SOLVE THE RESIDUE SYSTEM AX=B (MOD PP)
C FOR Y=A(ADJ)*B (MOD PP) AND DET (MOD PP)
C AND STORE RESULTS IN MULTL
      CALL SOLVE(ATEMP,MULTL,N,IN,MM,M,IMPIN,IMIN1,NOPRIM)
C IF IQUIT IS NOT EQUAL TO 0, THEN THE SYSTEM IS
C SINGULAR MODULO EACH OF THE STORED PRIMES,
C AND HENCE, CANNOT BE SOLVED BY THIS PROGRAM.
C RETURN AN ERROR CODE OF 2.
      IF(IQUIT .EQ. 0) GO TO 40
      IER=2
      RETURN
C IF IFLAG IS NOT EQUAL TO 0, THEN A IS SINGULAR
C MODULO KPRIME(ICOUNT).  CHOOSE ANOTHER PRIME,
C I.E. KPRIME(ICOUNT+1), AND TRY TO SOLVE
C THE SYSTEM AGAIN.
   40 IF(IFLAG .NE. 0) GO TO 50
      SUMLOG=SUMLOG+ALOG(P)
C TEST TO SEE IF THE REQUIRED NUMBER
C OF PRIMES HAVE BEEN USED
      IF(SUMLOG .GE. BOUND) GO TO 60
      IS=IS+1
   50 ICOUNT=ICOUNT+1
C IF ALL PRIMES HAVE BEEN TRIED
C AND STILL ANOTHER IS REQUIRED,
C COMBINE RESULTS AND CHECK SOLUTION
      IF(ICOUNT .LE. NOPRIM) GO TO 10
      IS=IS-1
C COMBINE RESULTS BY CONVERSION
C TO SYMMETRIC MIXED-RADIX
   60 CALL MXRADX(MULTL,MM,RY,LCOUNT,NDIGIT,IMIN1,NOPRIM,
     1NO2,W)
C CHECK SOLUTION BY COMPUTING AY AND DB
      CALL CHECK(A,N,IN,B,M,IM,IER,MULTL,IMIN1,NOPRIM,
     1NO2,W,V)
      IF(IER .EQ. 1) RETURN
C COMPUTE THE SOLUTION X = (1/DET)*Y
      DET=RY(NORES)
      INDEX=0
      DO 70 J=1,M
        DO 70 I=1,N
          INDEX=INDEX+1
   70     X(I,J)=RY(INDEX)/DET
      RETURN
C RETURN ERROR CODE OF 3 FOR INCONSISTENT
C INPUT PARAMETERS
   80 IER=3
      RETURN
      END
      SUBROUTINE SOLVE(ATEMP,MULTL,N,IN,MM,M,IMPIN,IMIN1,
     1NOPRIM)
      DIMENSION MM(NOPRIM),MULTL(IMIN1,NOPRIM),ATEMP(IN,IMPIN)
      INTEGER ATEMP,PP,RESIDU
      COMMON/MLEN/IB,PP,NZ,IS,IFLAG,IQUIT,NORES
C THIS SUBROUTINE SOLVES THE RESIDUE SYSTEM
C AX=B (MOD PP) FOR Y (MOD PP) AND DET (MOD PP).
      IDET=1
C FIND A PIVOTAL ELEMENT RELATIVELY PRIME TO PP
      MPN=M+N
      DO 80 J=1,N
        DO 10 I=J,N
          IF(MOD(ATEMP(I,J),PP) .NE. 0) GO TO 20
          IF(I .EQ. N) GO TO 100
   10     CONTINUE
C PERMUTE ROWS I AND J
   20   IF(I .EQ. J) GO TO 40
        IDET=-IDET
        DO 30 JJ=J,MPN
          ITEMP=ATEMP(J,JJ)
          ATEMP(J,JJ)=ATEMP(I,JJ)
   30     ATEMP(I,JJ)=ITEMP
C ACCUMULATE DETERMINANT
   40   IDET=IDET*ATEMP(J,J)
        IDET=MOD(IDET,PP)
C FIND INVERSE OF PIVOTAL ELEMENT
        IX=INVERS(ATEMP(J,J),PP)
C MULTIPLY ROW J BY INVERSE OF PIVOTAL ELEMENT
        DO 50 JJ=J,MPN
          ITEMP=ATEMP(J,JJ)*IX
   50     ATEMP(J,JJ)=MOD(ITEMP,PP)
C REPLACE LTH ROW BY LTH ROW-JTH ROW, (L NOT EQUAL J)
        DO 70 L=1,N
          IF(L .EQ. J) GO TO 70
          IK=ATEMP(L,J)
          DO 60 JJ=J,MPN
            ITEMP=ATEMP(J,JJ)*IK
            ITEMP=MOD(ITEMP,PP)
            ITEMP=ATEMP(L,JJ)-ITEMP
            ATEMP(L,JJ)=MOD(ITEMP,PP)
   60       CONTINUE
   70     CONTINUE
   80   CONTINUE
C STORE SYMMETRIC RESIDUE DIGITS IN MULTL,
C AND MODULUS IN MM
      N1=N+1
      INDEX=0
      DO 90 J=N1,MPN
        DO 90 I=1,N
          INDEX=INDEX+1
          ITEMP=ATEMP(I,J)*IDET
   90     MULTL(INDEX,IS)=RESIDU(ITEMP,PP)
      MULTL(NORES,IS)=RESIDU(IDET,PP)
      MM(IS)=PP
      RETURN
  100 NZ=NZ+1
      IFLAG=1
C TEST TO SEE IF ALL PRIMES HAVE FAILED
      IF(NZ .GT. NOPRIM-1) IQUIT=1
      RETURN
      END
      SUBROUTINE MXRADX(MULTL,MM,RY,LCOUNT,NDIGIT,IMIN1,
     1NOPRIM,NO2,W)
      DIMENSION MM(NOPRIM),MULTL(IMIN1,NOPRIM),RY(IMIN1),
     1W(NO2)
      INTEGER RESIDU,W,PP
      DOUBLE PRECISION ACC,ACC1,ACC2,TEX
      COMMON/MLEN/IB,PP,NZ,IS,IFLAG,IQUIT,NORES
C SUBROUTINE MXRADX COMPUTES THE SYMMETRIC
C MIXED-RADIX DIGITS OF Y AND DET A FROM
C THEIR SYMMETRIC RESIDUE DIGITS.
      IF(IS .EQ. 1) GO TO 150
C COMPUTE SYMMETRIC MIXED-RADIX DIGITS
C AND STORE THEM IN MULTL
      DO 10 I=2,IS
        KK=I-1
        DO 10 J=I,IS
          IX=INVERS(MM(KK),MM(J))
          DO 10 K=1,NORES
            ITEMP=MULTL(K,J)-MULTL(K,I-1)
            ITEMP=ITEMP*IX
   10       MULTL(K,J)=RESIDU(ITEMP,MM(J))
C COMPUTE Y AND D FROM
C THEIR SYMMETRIC MIXED-RADIX DIGITS
C USING MULTILENGTH ARITHMETIC
      LCOUNT=0
   20 DO 140 I=1,NORES
        W(1)=1
        DO 30 K=2,IS
   30     W(K)=0
C COMPUTE Y(I)=(...(MULTL(I,IS)*MM(IS-1)+
C MULTL(I,IS-1))*MM(IS-2)+...+MULTL(I,2))
C MM(1)+MULTL(I,1)
        W(1)=W(1)*MULTL(I,IS)*MM(IS-1)
        CALL MLTLTH(NO2,W)
        J=IS
   40   J=J-1
        IF(J .LE. 1) GO TO 60
        W(1)=W(1)+MULTL(I,J)
        CALL MLTLTH(NO2,W)
        DO 50 K=1,IS
   50     W(K)=W(K)*MM(J-1)
        CALL MLTLTH(NO2,W)
        GO TO 40
   60   W(1)=W(1)+MULTL(I,J)
        CALL MLTLTH(NO2,W)
C STORE MULTILENGTH DIGITS OF Y(I)
C IN MULTL(I,J),J=1,IS
C STORE MULTILENGTH DIGITS OF DET A
C IN MULTL(NORES,J),J=1,IS
        DO 70 J=1,IS
   70     MULTL(I,J)=W(J)
C COMPUTE Y(I) IN FLOATING-PT. FROM MULTILENGTH DIGITS
        K=IS
   80   IF(W(K) .NE. 0) GO TO 90
        IF(K .EQ. 1) GO TO 100
        K=K-1
        GO TO 80
   90   IF(K .LE. 1) GO TO 100
        ACC=W(K)*IB+W(K-1)
        TEX=NDIGIT*(K-2)
        GO TO 110
  100   RY(I)=W(1)
        GO TO 130
  110   IF(K .LE. 2) GO TO 120
        ACC1=W(K-2)
        ACC2=10.D0**NDIGIT
        ACC=ACC+ACC1/ACC2
  120   ACC1=1.0D+1**TEX
        RY(I)=ACC*ACC1
  130   IF(K .LE. LCOUNT) GO TO 140
        LCOUNT=K
  140   CONTINUE
      RETURN
  150 DO 160 I=1,NORES
  160   RY(I)=MULTL(I,1)
      RETURN
      END
      SUBROUTINE MLTLTH(NO2,W)
      DIMENSION W(NO2)
      INTEGER W,PP
      COMMON/MLEN/IB,PP,NZ,IS,IFLAG,IQUIT,NORES
      IF(IS .EQ. 1) RETURN
      L=IS-1
C DISTRIBUTE THE DIGITS IN W SO THAT
C EACH ELEMENT OF W CONTAINS NDIGIT DIGITS
      DO 10 K=1,L
        W(K+1)=W(K)/IB+W(K+1)
   10   W(K)=-W(K)/IB*IB+W(K)
      K=IS
C ALL THE ELEMENTS OF W SHOULD HAVE THE SAME SIGN.
   20 IF(W(K))60,30,40
   30 IF(K .EQ. 1) RETURN
      K=K-1
      GO TO 20
   40 DO 50 K=1,L
        IF(W(K) .GE. 0) GO TO 50
        W(K)=W(K)+IB
        W(K+1)=W(K+1)-1
   50   CONTINUE
      RETURN
   60 DO 70 K=1,L
        IF(W(K) .LE. 0) GO TO 70
        W(K)=W(K)-IB
        W(K+1)=W(K+1)+1
   70   CONTINUE
      RETURN
      END
      SUBROUTINE CHECK(A,N,IN,B,M,IM,IER,MULTL,IMIN1,
     1NOPRIM,NO2,W,V)
      DIMENSION V(NO2),MULTL(IMIN1,NOPRIM),A(IN,IN),
     1B(IN,IM),W(NO2)
      INTEGER W,V,A,B,PP
      COMMON/MLEN/IB,PP,NZ,IS,IFLAG,IQUIT,NORES
C SUBROUTINE CHECK CHECKS THE SOLUTION BY COMPUTING
C A*Y AND (DET A)*B AND COMPARING THE RESULTS.
C Y IS STORED BY COLUMNS IN MULTL
C DET IS STORED IN MULTL(NORES,I),I=1,IS
      LL=IS
      KK=IS+1
      DO 90 I=1,N
        INDEX=0
        DO 90 L=1,M
C MULTIPLY ROW I OF A BY COLUMN L OF Y
          DO 10 K=1,NO2
   10       W(K)=0
          IS=KK
          DO 40 J=1,N
            INDEX=INDEX+1
            JJ=A(I,J)/IB
            II=-JJ*IB+A(I,J)
            IF(LL .EQ. 1) GO TO 30
            DO 20 K=2,LL
              W(K)=W(K)+MULTL(INDEX,K)*II+MULTL(INDEX,K-1)*JJ
              CALL MLTLTH(NO2,W)
   20         CONTINUE
   30       W(1)=II*MULTL(INDEX,1)+W(1)
            W(KK)=JJ*MULTL(INDEX,LL)+W(KK)
            CALL MLTLTH(NO2,W)
            IF(IABS(W(IS)) .LT. IB) GO TO 40
            IF(IS .GE. NO2) GO TO 40
            IS=IS+1
            W(IS)=W(IS-1)/IB
            W(IS-1)=-W(IS)*IB+W(IS-1)
   40       CONTINUE
C STORE THE PRODUCT IN V
          DO 50 K=1,IS
   50       V(K)=W(K)
C MULTIPLY B(I,L) BY DET AND STORE IN W
          JJ=B(I,L)/IB
          II=-JJ*IB+B(I,L)
          IF(LL .EQ. 1) GO TO 70
          DO 60 K=2,LL
   60       W(K)=MULTL(NORES,K)*II+MULTL(NORES,K-1)*JJ
   70     W(1)=II*MULTL(NORES,1)
          W(KK)=JJ*MULTL(NORES,LL)+W(KK)
          CALL MLTLTH(NO2,W)
C TEST EQUALITY OF W AND V
          DO 80 J=1,IS
            IF(W(J) .NE. V(J)) GO TO 100
   80       CONTINUE
   90     CONTINUE
C IF SOLUTION CHECKS, RETURN IER=0
C ELSE, RETURN IER=1
      IER=0
      IS=LL
      RETURN
  100 IER=1
      RETURN
      END
      FUNCTION INVERS(K,M)
C INVERS COMPUTES AN INVERSE OF K (MOD M)
C BY THE EUCLIDEAN ALGORITHN
      I=K
      L=M
      J=1
      INVERS=0
   10 KK = I/L
      NN=MOD(I,L)
      IF(NN .EQ. 0) GO TO 20
      I=L
      L=NN
      NN=-KK*INVERS+J
      J=INVERS
      INVERS=NN
      GO TO 10
   20 IF(L .GE. 0) GO TO 30
      INVERS=-INVERS
C RETURN A POSITIVE VALUE
   30 IF(INVERS .GE. 0) RETURN
      INVERS=M+INVERS
      RETURN
      END
      INTEGER FUNCTION RESIDU(K,M)
      RESIDU=MOD(K,M)
C THE FUNCTION RESIDU COMPUTES THE SYMMETRIC
C RESIDUE OF K (MOD M)
C I.E.  -M/2 LESS THAN RESIDU LESS THAN M/2
      IF(RESIDU)10,20,30
   10 IF(2*RESIDU+M .GE. 0) RETURN
      RESIDU=RESIDU+M
   20 RETURN
   30 IF(-2*RESIDU+M .GE. 0) RETURN
      RESIDU=RESIDU-M
      RETURN
      END
      SUBROUTINE LOGBND(A,N,IN,B,M,IM,BOUND)
      DIMENSION A(IN,IN),B(IN,IM)
      INTEGER A,B
C BOUND IS A LOWER BOUND FOR THE
C LOG OF THE PRODUCT OF THE MODULI
      BOUND=0.
      DO 20 I=1,N
        ALPHA=0.
        DO 10 J=1,N
          TEMP=A(I,J)
          TEMP=TEMP*TEMP
   10     ALPHA=ALPHA+TEMP
   20   BOUND=BOUND+ALOG(ALPHA)
      BOUND=BOUND/2.
      DO 30 J=1,M
        DO 30 I=1,N
          ALPHA=IABS(B(I,J))
          IF(ALPHA .EQ. 0.) GO TO 30
          BOUND=BOUND+ALOG(ALPHA)
   30     CONTINUE
   40 BOUND=BOUND+ALOG(2.)
      RETURN
      END