C                                                                       MAIN0010
C***************************************************************        MAIN0020
C                                                              *        MAIN0030
C                      TWO TEST EXAMPLES                       *        MAIN0040
C                      -----------------                       *        MAIN0050
C                                                              *        MAIN0060
C***************************************************************        MAIN0070
C                                                                       MAIN0080
      INTEGER P(100),AMOD(20,21),Y(32,20,1),DET(32),T,TWOT1,            MAIN0090
     1        SUM(7),GAMMA(7),DELTA(7),A(10),C(10),B                    MAIN0100
      INTEGER A1(1,5,5),B1(1,5,1),A2(3,20,20),B2(3,20,1),TEMP1,Q        MAIN0110
      REAL S(21)                                                        MAIN0120
      COMMON /PRIMEB/ P, IPRIME(100)                                    MAIN 130
      B = 10000                                                         MAIN0140
C                                                                       MAIN0150
C***************************************************************        MAIN0160
C                                                              *        MAIN0170
C                  EXAMPLE 1 - PASCAL MATRIX.                  *        MAIN0180
C                                                              *        MAIN0190
C                      A1*X = B1.                              *        MAIN0200
C                                                              *        MAIN0210
C                        T = 1                                 *        MAIN0220
C                        N = 5                                 *        MAIN0230
C                        M = 1                                 *        MAIN0240
C                                                              *        MAIN0250
C***************************************************************        MAIN0260
C                                                                       MAIN0270
        T = 1                                                           MAIN0280
        N = 5                                                           MAIN0290
        M = 1                                                           MAIN0300
        NM = N+M
C
C     **GENERATE THE MATRICES A1 AND B1.
C
        DO 1 I=1,N
           A1(1,I,1) = 1
 1         A1(1,1,I) = 1
C       OD
        DO 2 I=2,N
           DO 2 J=2,N
 2            A1(1,I,J) = A1(1,I,J-1)+A1(1,I-1,J)
C          OD
C       OD
        B1(1,1,1) = 1
        DO 3 I=2,N
 3         B1(1,I,1) = B1(1,I-1,1)+A1(1,I,N)
C       OD
C
C     **PRINT A1 AND B1 IN FIXED-RADIX FORM.
C
        WRITE(6,42)
        DO 5 K=1,T
           WRITE(6,46) K
           DO 4 I=1,N
 4            WRITE(6,43) (A1(K,I,J),J=1,N)
C          OD
 5                                                              CONTINUE
C       OD
        WRITE(6,44)
        DO 7 K=1,T
           WRITE(6,49) K
           DO 6 I=1,N
 6            WRITE(6,43) (B1(K,I,J),J=1,M)
C          OD
 7                                                              CONTINUE
C       OD
C
C     **CONVERT A1 AND B1 FROM FIXED-RADIX TO MIXED-RADIX FORM.
C
        N1 = 1
        TEMP = (FLOAT(N1)*ALOG(FLOAT(B)) + ALOG(2.0))/ALOG(FLOAT(P(1)))
        LMAX = TEMP
        IF (TEMP-FLOAT(LMAX).GT.0.01) LMAX = LMAX+1
        DO 10 I=1,N
           DO 10 J=1,N
              DO 8 K=1,N1
 8               A(K) = A1(K,I,J)
C             OD
              CALL MRADIX(A,N1,C,LMAX,L,B,IER)
              DO 9 K=1,L
 9               A1(K,I,J) = C(K)
C             OD
C          OD
C       OD
 10                                                             CONTINUE
        DO 13 I=1,N
           DO 13 J=1,M
              DO 11 K=1,N1
 11               A(K) = B1(K,I,J)
C             OD
              CALL MRADIX(A,N1,C,LMAX,L,B,IER)
              DO 12 K=1,L
 12              B1(K,I,J) = C(K)
C             OD
C          OD
C       OD
 13                                                             CONTINUE
C
C     **PRINT A1 AND B1 IN MIXED-RADIX FORM.
C
        WRITE(6,45)
        DO 15 K=1,T
           WRITE(6,46) K
           DO 14 I=1,N
 14           WRITE(6,47) (A1(K,I,J),J=1,N)
C          OD
 15                                                             CONTINUE
C       OD
        WRITE(6,48)
        DO 17 K=1,T
           WRITE(6,49) K
           DO 16 I=1,N
 16           WRITE(6,47) (B1(K,I,J),J=1,M)
C          OD
 17                                                             CONTINUE
C       OD
C
C     **COMPUTE MAXPRM.
C
        TWOT1 = 2*T+1
        CALL SUBBND(A1,B1,GAMMA,DELTA,S,SUM,T,N,M,NM,
     1              TWOT1,BOUND,MAXPRM)
        WRITE(6,50) BOUND,MAXPRM
C
C     **SOLVE THE SYSTEM OF LINEAR EQUATIONS A1*X=B1.
C
        CALL DRIVER(A1,B1,AMOD,T,N,M,NM,Y,DET,MAXPRM,B)
C
C***************************************************************
C                                                              *
C                  EXAMPLE 2 - PASCAL MATRIX.                  *
C                                                              *
C                      A2*X = B2.                              *
C                                                              *
C                        T = 3                                 *
C                        N = 20                                *
C                        M = 1                                 *
C                                                              *
C***************************************************************
C
        T = 3
        N = 20
        M = 1
        NM = N+M
C
C     **GENERATE THE MATRICES A2 AND B2.
C
        DO 18 I=1,N
           A2(3,I,1) = 1
 18        A2(3,1,I) = 1
C       OD
        DO 19 K=1,2
           DO 19 I=1,N
              A2(K,I,1) = 0
 19           A2(K,1,I) = 0
C          OD
C       OD
        DO 21 I=2,N
           DO 21 J=2,N
              Q = 0
              TEMP1 = 0
              DO 20 K=1,T
                 K1 = T-K+1
                 TEMP1 = A2(K1,I,J-1)+A2(K1,I-1,J)+Q
                 Q = TEMP1/B
 20              A2(K1,I,J) = TEMP1 - Q*B
C             OD
C          OD
C       OD
 21                                                             CONTINUE
        B2(1,1,1) = 0
        B2(2,1,1) = 0
        B2(3,1,1) = 1
        DO 23 I=2,N
           Q = 0
           TEMP1 = 0
           DO 22 K=1,T
              K1 = T-K+1
              TEMP1 = B2(K1,I-1,1)+A2(K1,I,N)+Q
              Q = TEMP1/B
 22           B2(K1,I,1) = TEMP1 - Q*B
C          OD
C       OD
 23                                                             CONTINUE
C
C     **PRINT A2 AND B2 IN FIXED-RADIX FORM.
C
        WRITE(6,42)
        DO 25 K=1,T
           WRITE(6,46) K
           K1 = T-K+1
           DO 24 I=1,N
 24           WRITE(6,43) (A2(K1,I,J),J=1,N)
C          OD
C       OD
 25                                                             CONTINUE
        WRITE(6,44)
        DO 27 K=1,T
           WRITE(6,49) K
           K1 = T-K+1
           DO 26 I=1,N
 26           WRITE(6,43) (B2(K1,I,J),J=1,M)
C          OD
C       OD
 27                                                             CONTINUE
C
C     **CONVERT A2 AND B2 FROM FIXED-RADIX TO MIXED-RADIX FORM.
C
        N1 = 3
        TEMP = (FLOAT(N1)*ALOG(FLOAT(B)) + ALOG(2.0))/ALOG(FLOAT(P(1)))
        LMAX = TEMP
        IF (TEMP-FLOAT(LMAX).GT.0.01) LMAX = LMAX+1
        DO 32 I=1,N
           DO 32 J=1,N
              DO 28 K=1,N1
 28              A(K) = A2(K,I,J)
C             OD
              CALL MRADIX(A,N1,C,LMAX,L,B,IER)
              DO 29 K=1,L
 29              A2(K,I,J) = C(K)
C             OD
              IF (L.LT.N1)                                      GOTO 30
                                                                GOTO 32
C             THEN
 30              L1 = L+1
                 DO 31 K=L1,N1
 31                 A2(K,I,J) = 0
C                OD
C          OD
C       OD
 32                                                             CONTINUE
        DO 37 I=1,N
           DO 37 J=1,M
              DO 33 K=1,N1
 33              A(K) = B2(K,I,J)
C             OD
              CALL MRADIX(A,N1,C,LMAX,L,B,IER)
              DO 34 K=1,L
 34              B2(K,I,J) = C(K)
C             OD
              IF (L.LT.N1)                                      GOTO 35
                                                                GOTO 37
C             THEN
 35              L1 = L+1
                 DO 36 K=L1,N1
 36                 B2(K,I,J) = 0
C                OD
C          OD
C       OD
 37                                                             CONTINUE
C
C     **PRINT A2 AND B2 IN MIXED-RADIX FORM.
C
        WRITE(6,45)
        DO 39 K=1,T
           WRITE(6,46) K
           DO 38 I=1,N
 38           WRITE(6,47) (A2(K,I,J),J=1,N)
C          OD
 39                                                             CONTINUE
C       OD
        WRITE(6,48)
        DO 41 K=1,T
           WRITE(6,49) K
           DO 40 I=1,N
 40           WRITE(6,47) (B2(K,I,J),J=1,M)
C          OD
 41                                                             CONTINUE
C       OD
C
C     **COMPUTE MAXPRM.
C
        TWOT1 = 2*T+1
        CALL SUBBND(A2,B2,GAMMA,DELTA,S,SUM,T,N,M,NM,
     1              TWOT1,BOUND,MAXPRM)
        WRITE(6,50) BOUND,MAXPRM
C
C     **SOLVE THE SYSTEM OF LINEAR EQUATIONS A2*X=B2.
C
        CALL DRIVER(A2,B2,AMOD,T,N,M,NM,Y,DET,MAXPRM,B)
C
 42   FORMAT(1H1,16X,5HINPUT/17X,5(1H*)/5X,
     1       38HFIXED-RADIX REPRESENTATION, BASE=10000/10X,
     2       40HA = A(1,N,N)+A(2,N,N)*BASE+...+A(T,N,N)*,
     3       11HBASE**(T-1))
 43   FORMAT(5X,20(I4,2X))
 44   FORMAT(/10X,40HB = B(1,N,M)+B(2,N,M)*BASE+...+B(T,N,M)*,
     1       11HBASE**(T-1))
 45   FORMAT(/5X,26HMIXED-RADIX REPRESENTATION/10X,
     1       40HA = A(1,N,N)+A(2,N,N)*P(1)+...+A(T,N,N)*,
     2       15HP(1)*...*P(T-1))
 46   FORMAT(/10X,2HA(,I1,5H,N,N))
 47   FORMAT(1H ,20(I6))
 48   FORMAT(/10X,40HB = B(1,N,M)+B(2,N,M)*P(1)+...+B(T,N,M)*,
     1       15HP(1)*...*P(T-1))
 49   FORMAT(/10X,2HB(,I1,5H,N,M))
 50   FORMAT(/10X,13HLOG(BOUND) = ,F8.4,2X,9HMAXPRM = ,I3)
      STOP
      END
C                                                                       DRIV0010
C***************************************************************        DRIV0020
C                                                              *        DRIV0030
C                      CALL ESOLVE.                            *        DRIV0040
C                 PRINT OUTPUTS IN MIXED-RADIX FORM.           *        DRIV0050
C                      CALL FRADIX.                            *        DRIV0060
C                 PRINT OUTPUTS IN FIXED-RADIX FORM.           *        DRIV0070
C                                                              *        DRIV0080
C***************************************************************        DRIV0090
C                                                                       DRIV0100
      SUBROUTINE DRIVER(A,B,AMOD,T,N,M,NM,Y,DET,MAXPRM,BASE)            DRIV0110
      INTEGER T,A(T,N,N),B(T,N,M),AMOD(N,NM),Y(MAXPRM,N,M),
     1        DET(MAXPRM),BASE,P(100)
      INTEGER  A1(10),C1(10)
      COMMON /PRIMEB/ P, IPRIME(100)
      LOGICAL RTEST
C
C     **SOLVE THE SYSTEM OF LINEAR EQUATIONS AX=B.
C
        RTEST = .TRUE.
        CALL ESOLVE(A,B,AMOD,T,N,M,NM,DET,Y,MAXPRM,NOCOEF,
     1              RTEST,IER)
C
C     **PRINT OUTPUT PARAMETERS.
C
        WRITE(6,6) IER,NOCOEF
C
C     **PRINT Y IN MIXED-RADIX FORM.
C
        WRITE(6,7)
        DO 2 J=1,M
           WRITE(6,8) J
           DO 1 I=1,N
 1            WRITE(6,10) (Y(K,I,J),K=1,NOCOEF)
C          OD
 2                                                              CONTINUE
C       OD
C
C     **PRINT DET IN MIXED-RADIX FORM.
C
        WRITE(6,9)
        WRITE(6,10) (DET(K),K=1,NOCOEF)
C
C     **CONVERT Y AND DET FROM MIXED-RADIX TO FIXED-RADIX FORM.
C     **PRINT Y AND DET IN FIXED-RADIX FORM.
C
        TEMP = (FLOAT(NOCOEF)*ALOG(FLOAT(P(NOCOEF))) - ALOG(2.0)) /
     1         ALOG(FLOAT(BASE))
        LMAX = TEMP
        IF (TEMP-FLOAT(LMAX).GT.0.01) LMAX = LMAX+1
        WRITE(6,11)
        DO 5 J=1,M
           WRITE(6,8) J
           DO 4 I=1,N
              DO 3 K=1,NOCOEF
 3               C1(K) = Y(K,I,J)
C             OD
              CALL FRADIX(C1,NOCOEF,A1,LMAX,L,BASE,IER)
 4            WRITE(6,10) (A1(II),II=1,L)
C          OD
 5                                                              CONTINUE
C       OD
        CALL FRADIX(DET,NOCOEF,A1,LMAX,L,BASE,IER)
        WRITE(6,12) (A1(I),I=1,L)
C
 6    FORMAT(/16X,7HOUTPUTS/16X,7(1H*)/6X,
     1       6HIER = ,I3,5X,9HNOCOEF = ,I3)
 7    FORMAT(/16X,21HY IN MIXED-RADIX FORM/10X,
     1       35HFROM LEFT TO RIGHT, THE NUMBERS ARE,
     2       25H COEF(1),...,COEF(L), AND/10X,
     3       33HY(I,J) = COEF(1)+COEF(2)*P(1)+...,
     4       23H+COEF(L)*P(1)*...P(L-1))
 8    FORMAT(/10X,6HY(L,N,,I1,1H))
 9    FORMAT(/16X,23HDET IN MIXED-RADIX FORM/10X,
     1       35HFROM LEFT TO RIGHT, THE NUMBERS ARE,
     2       24H COEF(1),...,COEF(L),AND/10X,
     3       30HDET = COEF(1)+COEF(2)*P(1)+...,
     4       23H+COEF(L)*P(1)*...P(L-1))
 10   FORMAT(16X,20(I6))
 11   FORMAT(/16X,21HY IN FIXED-RADIX FORM/10X,
     1       35HFROM LEFT TO RIGHT, THE NUMBERS ARE,
     2       24H COEF(1),...,COEF(L),AND/10X,
     3       37HY(I,J) = COEF(1)*BASE**(L-1)+COEF(2)*,
     4       23HBASE**(L-2)+...+COEF(L))
 12   FORMAT(/16X,23HDET IN FIXED-RADIX FORM/10X,
     1       35HFROM LEFT TO RIGHT, THE NUMBERS ARE,
     2       24H COEF(1),...,COEF(L),AND/10X,
     3       34HDET = COEF(1)*BASE**(L-1)+COEF(2)*,
     4       21HBASE(N-2)+...+COEF(L)/16X,10(I6))
      RETURN
      END
      SUBROUTINE ESOLVE(A,B,AMOD,T,N,M,NM,DET,Y,MAXPRM,NOCOEF,          ESOL0010
     1                  RTEST,IER)
C
C***************************************************************
C                                                              *
C THIS SUBROUTINE SOLVES EXACTLY THE SYSTEM OF LINEAR EQUATIONS*
C AX=B, WHERE A( ,N,N) AND B( ,N,M) ARE MATRICES WITH MULTIPLE-*
C PRECISION INTEGER COEFFICIENTS EXPRESSED IN MIXED-RADIX FORM.*
C THE SUBROUTINE NORMALLY RETURNS INTEGER VALUES FOR           *
C DET=DETERMINANT(A) AND Y=A(ADJOINT)*B, ALSO IN MIXED-RADIX   *
C FORM. IF THE SOLUTION X IS REQUIRED, THE USER NEED ONLY      *
C COMPUTE X=Y/DET.  (FOR CONVERSION OF INTEGERS FROM MIXED-    *
C RADIX TO FIXED-RADIX FORMS AND CONVERSELY, SUBROUTINES FRADIX*
C AND MRADIX CAN BE USED.)                                     *
C                                                              *
C  PRIME  IS A LINEAR ARRAY CONTAINING 100 DISTINCT PRIME      *
C (INPUT) INTEGERS IN ASCENDING ORDER.  THE PRIMES ARE CHOSEN  *
C         AS LARGE AS POSSIBLE SUBJECT TO THE CONDITION THAT   *
C         FOR ALL I AND J PRIME(I)*PRIME(J) DOES NOT OVERFLOW  *
C         AN INTEGER WORD.  THESE PRIMES ARE USED AS RADII FOR *
C         THE REPRESENTATION OF INTEGERS IN MIXED-RADIX FORM.  *
C IPRIME  IS A LINEAR ARRAY OF INTEGERS SUCH THAT FOR EACH K,  *
C (INPUT) IPRIME(K)* PRIME(1)*PRIME(2)*...*PRIME(K-1) = 1,     *
C                                          MODULO PRIME(K).    *
C      A  IS AN INTEGER MATRIX OF DIMENSION T BY N BY N.  THE  *
C (INPUT) FIRST DIMENSION CONTAINS THE COEFFICIENTS OF THE     *
C         MIXED-RADIX REPRESENTATION OF THE MULTIPLE-PRECISION *
C         COMPONENTS OF THE N BY N MATRIX A( ,N,N).  THAT IS,  *
C            A( ,I,J) = A(1,I,J)                               *
C                     + A(2,I,J)*PRIME(1)                      *
C                     .                                        *
C                     .                                        *
C                     .                                        *
C                     + A(T,I,J)*PRIME(1)*...*PRIME(T-1),      *
C         FOR I,J=1,2,...,N.                                   *
C      B  IS AN INTEGER MATRIX OF DIMENSION T BY N BY M WITH   *
C (INPUT) A SIMILAR NOTATIONAL CORRESPONDENCE MADE FOR A ABOVE.*
C   AMOD  IS AN N BY N+M MATRIX USED FOR TEMPORARY STORAGE OF  *
C  (TEMP) THE AUGMENTED MATRIX (A,B) MODULO THE VARIOUS PRIMES *
C         PRIME(K).  THIS MATRIX IS INCLUDED IN THE ARGUMENT   *
C         LIST ONLY TO PERMIT ITS DIMENSIONS TO BE VARIABLE.   *
C      T  IS THE NUMBER OF RADII, PRIME(1),...,PRIME(T), USED  *
C (INPUT) TO REPRESENT EACH COMPONENT OF A( ,N,N) AND B( ,N,M) *
C         IN MIXED-RADIX FORM.                                 *
C      N  IS THE NUMBER OF EQUATIONS AND THE NUMBER OF UNKNOWNS*
C (INPUT) IN THE SYSTEM.  (I.E., N IS THE SIZE OF THE SECOND   *
C         AND THIRD DIMENSIONS OF A.)                          *
C      M  IS THE NUMBER OF RIGHT-HAND VECTORS FOR WHICH THE    *
C (INPUT) SYSTEM IS TO BE SOLVED.  (I.E., M IS THE SIZE OF THE *
C         THIRD DIMENSION OF B.)                               *
C     NM  IS EQUAL TO N+M.                                     *
C (INPUT)                                                      *
C    DET  IS A VECTOR OF DIMENSION MAXPRM WHICH CONTAINS THE   *
C(OUTPUT) COEFFICIENTS OF THE MIXED-RADIX REPRESENTATION OF    *
C         DETERMINANT(A)                                       *
C            = DET(1)                                          *
C            + DET(2)*PRIME(1)                                 *
C            .                                                 *
C            .                                                 *
C            .                                                 *
C            + DET(MAXPRM)*PRIME(1)*...*PRIME(MAXPRM-1).       *
C      Y  IS A MATRIX OF DIMENSION MAXPRM BY N BY M.  THE FIRST*
C(OUTPUT) DIMENSION CONTAINS THE COEFFICIENTS OF THE MIXED-    *
C         RADIX REPRESENTATION OF Y( ,N,M) = A(ADJOINT)*B:     *
C            Y( ,I,J)                                          *
C               = Y(1,I,J)                                     *
C               + Y(2,I,J)*PRIME(1)                            *
C               .                                              *
C               .                                              *
C               .                                              *
C               + Y(MAXPRM,I,J)*PRIME(1)*...*PRIME(MAXPRM-1).  *
C MAXPRM  IS THE MAXIMUM NUMBER OF RADII, PRIME(1),...,        *
C (INPUT) PRIME(MAXPRM), REQUIRED TO REPRESENT DETERMINANT(A)  *
C         AND A(ADJOINT)*B IN MIXED-RADIX FORM.  MAXPRM SHOULD *
C         BE CHOSEN SO THAT                                    *
C            ABS(DETERMINANT(A)), ABS(Y( ,I,J))                *
C               < PRIME(1)*...*PRIME(MAXPRM)/2, I=1,2,...,N,   *
C                                               J=1,2,...,M.   *
C         CHOOSING MAXPRM = N*T + N*LOG(N)/(2*LOG(PRIME(1)))   *
C         WILL SUFFICE. (A TIGHTER BOUND CAN BE OBTAINED USING *
C         SUBROUTINE SUBBND.) MAXPRM MUST BE LESS THAN 101.    *
C NOCOEF  IS THE NUMBER OF RADII ACTUALLY USED TO REPRESENT    *
C(OUTPUT) DET AND Y( ,I,J) IN MIXED-RADIX FORM.                *
C         NOCOEF <= MAXPRM.                                    *
C  RTEST  IF RTEST = .FALSE. , MAXPRM COEFFICIENTS IN THE      *
C (INPUT)    MIXED-RADIX REPRESENTATION OF DET AND Y( ,N,M)    *
C            ARE COMPUTED.                                     *
C         IF RTEST = .TRUE. , THE ALGORITHM CONTINUES TO       *
C            ITERATE UNTIL TR  (USUALLY TR=T, BUT SOMETIMES    *
C            TR=T+1)  SUCCESSIVE ZERO COEFFICIENTS IN THE      *
C            MIXED-RADIX REPRESENTATION OF DET AND Y( ,N,M)    *
C            ARE ENCOUNTERED, THAT IS, WHENEVER FOR            *
C            K = NOCOEF+1,...,NOCOEF+TR  (NOCOEF+TR <= MAXPRM),*
C               DET(K) = 0                                     *
C               Y(K,I,J) = 0, I=1,2,...,N, J=1,2,...,M.        *
C            IF NOCOEF > 0, THE PROGRAM TERMINATES WITH THE    *
C               KNOWLEDGE THAT                                 *
C               (1) A IS SINGULAR IF DET=0, OR                 *
C               (2) Y( ,N,M)/DET IS THE SOLUTION TO AX=B IF    *
C                   DET =0. (*WARNING* IN THIS CASE, A COMMON  *
C                   FACTOR IN Y( ,N,M) AND DET MAY HAVE        *
C                   BEEN IGNORED.  THAT IS, DET MIGHT NO       *
C                   LONGER BE THE DETERMINANT OF A.)           *
C            IF NOCOEF = 0, A IS PROBABLY SINGULAR, BUT THE    *
C               PROGRAM CONTINUES TO ITERATE UNTIL MAXPRM      *
C               COEFFICIENTS OF DET AND Y( ,N,M) HAVE BEEN     *
C               COMPUTED.                                      *
C    IER  IS AN ERROR CODE WHICH IS                            *
C(OUTPUT) 0  IF Y( ,N,M)/DET IS THE SOLUTION TO AX=B.          *
C            (*WARNING*  DET MIGHT NOT BE THE DETERMINANT OF A)*
C         1  IF MAXPRM COEFFICIENTS IN THE MIXED-RADIX         *
C            REPRESENTATION OF DET AND Y( ,N,M) HAVE BEEN      *
C            COMPUTED.  SINCE THE USER CHOOSES THE VALUE OF    *
C            MAXPRM, THE PROGRAM CANNOT GUANENTEE THAT         *
C            Y( ,N,M)/DET IS THE SOLUTION TO AX=B.             *
C         2  IF A IS SINGULAR.                                 *
C         3  IF A IS SINGULAR, OR DET AND Y( ,N,M) ARE NON-ZERO*
C            MULTIPLES OF PRIME(1)*...*PRIME(MAXPRM).          *
C         4  IF THE INPUT PARAMETERS ARE INCORRECT.            *
C                                                              *
C***************************************************************
C
      INTEGER T,TR,A(T,N,N),B(T,N,M),DET(MAXPRM),Y(MAXPRM,N,M),
     1        AMOD(N,NM)
      LOGICAL RTEST
      INTEGER PT,TK,TT,TT1,DMOD,V1,V3,U1,U3,T1,T3,Q
      INTEGER PRIME(100),IPRIME(100),P,P1,P2
      LOGICAL TEST
      COMMON /PRIMEB/PRIME,IPRIME
      IER=4
      IF (MAXPRM.GT.100 .OR. N.LT.2 .OR. M.LT.1) RETURN
      IF (N+M.NE.NM) RETURN
      N1 = N+1
C
C     PROCEDURE ZERO-CRITERION *********************************
C        DETERMINED IN THIS PROCEDURE IS THE NUMBER TR OF      *
C        CONSECUTIVE ZERO COEFFICIENTS IN THE MIXED-RADIX      *
C        REPRESENTATION OF DET AND Y REQUIRED TO GUARANTEE THAT*
C        X=Y/DET IS THE SOLUTION OF AX=B. TR=T IF NORM INFINITY*
C        OF THE MATRIX A AUGMENTED BY ANY COLUMN VECTOR OF B IS*
C        BOUNDED BY 2*PRIME(1)*...*PRIME(T); OTHERWISE, TR=T+1.*
C                                                              *
         TR = T
         IF (RTEST)                                             GOTO 1
                                                                GOTO 6
C        THEN
  1         PT = 2*PRIME(T)
            I = 0
C           REPEAT ... FOR EACH ROW I
  2            I = I+1
               IF (I.LE.N .AND. TR.EQ.T)                        GOTO 3
                                                                GOTO 6
C              THEN
  3               NORM = 0
                  DO 4 J=1,N
  4                   NORM = NORM + IABS(A(T,I,J))
C                 OD
                  MAXB = 0
                  DO 5 J=1,M
  5                  IF (MAXB.LT.IABS(B(T,I,J)))
     1                   MAXB = IABS(B(T,I,J))
C                 OD
                  NORM = NORM + MAXB
                  IF (T.GT.1) NORM = NORM + N1
C                 **HERE, USED IS THE FACT THAT FOR T>1
C                 **(ABS(K(T)+1)*PRIME(1)*...*PRIME(T-1)) IS
C                 **A BOUND FOR ANY INTEGER GIVEN BY
C                 **K(1)+K(2)*PRIME(1)+...+K(T)*PRIME(1)*
C                 **...*PRIME(T-1).
                  IF (NORM.GT.PT) TR = T+1
C            CONTINUE
                                                                GOTO 2
C                                                              *
C        END ZERO-CRITERION ************************************
C
  6   NOZERO = 0
      DO 64 ITER=1,MAXPRM
C        **THE SYSTEM AX=B IS SOLVED MODULO PRIME(ITER) FOR ALL
C        **ITER = 1,2,...,MAXPRM, OR IF RTEST IS TRUE UNTIL
C        **THE NUMBER OF CONSECUTIVE ZERO COEFFICIENTS (NOZERO)
C        **IN THE MIXED-RADIX REPRESENTATION OF DET AND Y IS
C        **EQUAL TO TR.
         P = PRIME(ITER)
         P1 = P-1
         P2 = P1/2 + 1
         IP = IPRIME(ITER)
C
C        PROCEDURE MAP *****************************************
C           THIS PROCEDURE COMPUTES AMOD, THE AUGMENTED        *
C           N BY N+M MATRIX (A,B) MODULO P.  THAT IS, FOR      *
C           I=1,2,...,N, HORNER'S RULE AND MODULO P ARITHMETIC *
C           IS USED TO COMPUTE                                 *
C             AMOD(I,J) = A(1,I,J) + A(2,I,J)*PRIME(1) +       *
C                       ... + A(T,I,J)*PRIME(1)*...*PRIME(T-1),*
C                                                  J=1,2,...,N,*
C           AMOD(I,J+N) = B(1,I,J) + B(2,I,J)*PRIME(1) +       *
C                       ... + B(T,I,J)*PRIME(1)*...*PRIME(T-1),*
C                                                  J=1,2,...,M.*
C                                                              *
            TT = MIN0(T,ITER)
            TT1 = TT-1
            IF (TT.EQ.1)                                        GOTO 7
                                                                GOTO 11
C           THEN
  7            DO 10 I=1,N
                  DO 8 J=1,N
  8                  AMOD(I,J) = A(1,I,J)
C                 OD
                  DO 9 J=1,M
                     NJ = N+J
  9                  AMOD(I,NJ) = B(1,I,J)
C                 OD
C              OD
 10                                                             CONTINUE
                                                                GOTO 17
C           ELSE
 11            DO 16 I=1,N
                  DO 13 J=1,N
                     NTEMP = A(TT,I,J)
                     DO 12 K=1,TT1
                        TK = TT - K
 12                     NTEMP = MOD(PRIME(TK)*NTEMP+A(TK,I,J),P)
C                    OD
 13                  AMOD(I,J) = NTEMP
C                 OD
                  DO 15 J=1,M
                     NTEMP = B(TT,I,J)
                     DO 14 K=1,TT1
                        TK = TT - K
 14                     NTEMP = MOD(PRIME(TK)*NTEMP+B(TK,I,J),P)
C                    OD
                     NJ = N+J
 15                  AMOD(I,NJ) = NTEMP
C                 OD
C              OD
 16                                                             CONTINUE
C                                                              *
C        END MAP ***********************************************
C
C        PROCEDURE MSOLVE **************************************
C           THIS PROCEDURE SOLVES THE SYSTEM OF EQUATIONS      *
C           AX=B MODULO P BY COMPUTING                         *
C              DMOD = DETERMINANT(A) MOD P                     *
C           AND                                                *
C              A(ADJOINT)*B MOD P,                             *
C           WHICH IS TEMPORARILY STORED IN Y(ITER,I,J),        *
C           I=1,2,...,N, J=1,2,...,M.                          *
C                                                              *
C           THE GAUSSIAN ELIMINATION METHOD WITH PARTIAL       *
C           PIVOTING IS USED TO REDUCE A TO UPPER ECHELON      *
C           FORM.                                              *
C                                                              *
C           NSING ASSUMES THE VALUE                            *
C              0, IF A HAS RANK N, MOD P                       *
C              1, IF A HAS RANK N-1, MOD P                     *
C              2, IF A HAS RANK SMALLER THAN N-1, MOD P        *
C                                                              *
C           KRAM IS EQUAL TO N IF A IS NONSINGULAR MODULO P;   *
C              OTHERWISE, KRAM IS THE INDEX OF THE FIRST       *
C              COLUMN WHICH CONTAINS A ZERO PIVOT ELEMENT.     *
C                                                              *
 17         DMOD = 1
            NSING = 0
            II = 0
C           REPEAT ... FORWARD ELIMINATION WITH II AS THE PIVOT
C                  ... COLUMN
 18            II = II+1
               IF (NSING.LT.2 .AND. II.LE.N)                    GOTO 19
                                                                GOTO 34
C              THEN
 19               I = II - NSING
C                 **ELIMINATION PROCEEDS ON THE SUBMATRIX WITH
C                 **ROWS I,I+1,...,N AND COLUMNS II,II+1,...,N+M
                  I1 = I+1
                  II1 = II+1
                  K = I-1
C                 REPEAT ... SEARCH FOR NON-ZERO ELEMENT IN
C                        ... II'TH COLUMN.
 20                  K = K+1
                     IF (K.LE.N)                                GOTO 21
                                                                GOTO 22
C                       THEN
 21                        IF (AMOD(K,II).NE.0)                 GOTO 22
                                                                GOTO 20
C                          THEN
C                <------------EXIT
C                 CONTINUE
 22               IF (K.LE.N)                                   GOTO 23
                                                                GOTO 33
C                 THEN ... AMOD(K,II) IS THE PIVOT ELEMENT
 23                  IF (K.NE.I)                                GOTO 24
                                                                GOTO 26
C                    THEN ... INTERCHANGE ROWS I AND K
 24                     DO 25 J=II,NM
                           NTEMP = AMOD(I,J)
                           AMOD(I,J) = AMOD(K,J)
 25                        AMOD(K,J) = NTEMP
C                       OD
                        DMOD = -DMOD
C
C                    PROCEDURE INVERT **************************
C                       EUCLID'S EXTENDED ALGORITHM IS USED TO *
C                       FIND THE INVERSE, VI, MODULO P OF      *
C                       AMOD(I,II).                            *
 26                     V1 = 1
                        U1 = 0
                        V3 = IABS(AMOD(I,II))
                        U3 = P
C                       REPEAT
 27                        IF (V3.NE.1)                         GOTO 28
                                                                GOTO 29
C                          THEN
 28                           Q = U3/V3
                              T1 = U1 - Q*V1
                              T3 = U3 - Q*V3
                              U1 = V1
                              U3 = V3
                              V1 = T1
                              V3 = T3
                                                                GOTO 27
C                       CONTINUE
 29                     IF (AMOD(I,II).LT.0) V1 = -V1
C                                                              *
C                    END INVERT ********************************
C
                     DMOD = MOD(DMOD*AMOD(I,II),P)
                     DO 30 J=II1,NM
 30                     AMOD(I,J) = MOD(V1*AMOD(I,J),P)
C                    OD
                     IF (I.LT.N)                                GOTO 31
                                                                GOTO 34
C                    THEN
 31                     DO 32 K=I1,N
                           NTEMP = AMOD(K,II)
                           DO 32 J=II1,NM
 32                           AMOD(K,J) = MOD(AMOD(K,J) - NTEMP*
     1                                        AMOD(I,J),P)
C                          OD
C                       OD
                                                                GOTO 18
C                 ELSE ... COLUMN II IS A ZERO PIVOT COLUMN.
 33                  NSING = NSING + 1
                     IF (NSING.LT.2) KRAM = II
                                                                GOTO 18
C           CONTINUE
 34         IF (NSING.LT.2)                                     GOTO 35
                                                                GOTO 42
C           THEN ... RANK OF A MODULO P IS N OR N-1; THEREFORE
C                ... BACK SUBSTITUTE.
 35            IF (NSING.EQ.0) KRAM = N
               KRAM1 = KRAM - 1
               KRAM2 = KRAM + 1
               DO 41 J=N1,NM
                  NJ = J - N
                  IF (KRAM.NE.N)                                GOTO 36
                                                                GOTO 38
C                 THEN
 36                  DO 37 I=KRAM2,N
 37                     Y(ITER,I,NJ) = 0
C                    OD
 38               NTEMP = MOD(DMOD*AMOD(N,J),P)
                  IF (NSING.EQ.1) NTEMP = NTEMP*(-1)**(N-KRAM)
                  Y(ITER,KRAM,NJ) = NTEMP
                  DO 40 II=1,KRAM1
                     I = KRAM - II
                     I1 = I+1
                     NTEMP = 0
                     DO 39 K=I1,KRAM
 39                     NTEMP = MOD(NTEMP+AMOD(I,K)*
     1                              Y(ITER,K,NJ), P)
C                    OD
                     NTEMP = -NTEMP
                     IF (NSING.EQ.0) NTEMP = MOD(NTEMP + DMOD*
     1                                           AMOD(I,J),P)
 40                  Y(ITER,I,NJ) = NTEMP
C                 OD
C              OD
 41                                                             CONTINUE
                                                                GOTO 44
C           ELSE ... RANK OF A MODULO P IS LESS THAN N-1.
 42            DO 43 J=1,M
                  DO 43 I=1,N
 43                  Y(ITER,I,J) = 0
C                 OD
C              OD
 44         IF (NSING.NE.0) DMOD = 0
C                                                              *
C        END MSOLVE ********************************************
C
C        PROCEDURE MIXED-RADIX *********************************
C           THIS PROCEDURE COMPUTES THE ITER'TH TERMS OF THE   *
C           MIXED-RADIX REPRESENTATION OF Y( ,N,M) AND DET     *
C           USING THE CHINESE REMAINDER THEOREM.               *
C                                                              *
C           TEST IS TRUE ONLY IF THE ITER'TH TERMS OF THE      *
C           MIXED-RADIX REPRESENTATIONS OF DET AND Y ARE       *
C           ALL ZERO.                                          *
C                                                              *
C           NOTE: IN FORTRAN, THE EVALUATION B=MOD(A,P) YEILDS *
C           AN INTEGER B SUCH THAT -P<B<P. IN THIS PROCEDURE,  *
C           ANY SUCH B IS CONVERTED TO A SYMMETRIC RESIDUE C,  *
C           -P<2C<P VIA THE EVALUATION C=B-(B/P2)*P.           *
C                                                              *
            TEST = .TRUE.
            ITER1 = ITER - 1
            ITER2 = ITER - 2
            DO 52 I=1,N
               DO 52 J=1,M
                  IF (ITER.EQ.1)                                GOTO 45
                                                                GOTO 46
C                 THEN
 45                  MULT = Y(ITER,I,J)
                                                                GOTO 50
C                 ELSE
 46                  MULT = Y(ITER1,I,J)
                     IF (ITER.NE.2)                             GOTO 47
                                                                GOTO 49
C                    THEN
 47                     DO 48 LL=1,ITER2
                           L = ITER1 - LL
 48                        MULT = MOD(MULT*PRIME(L)+Y(L,I,J),P)
C                       OD
 49                  MULT = MOD(IP * MOD(Y(ITER,I,J)-MULT,P),P)
 50               Y(ITER,I,J) = MULT - (MULT/P2)*P
                  IF (TEST)                                     GOTO 51
                                                                GOTO 52
C                 THEN
 51                  IF (Y(ITER,I,J).NE.0) TEST = .FALSE.
 52                                                             CONTINUE
C              OD
C           OD
            IF (ITER.EQ.1)                                      GOTO 53
                                                                GOTO 54
C           THEN
 53            MULT = DMOD
                                                                GOTO 58
C           ELSE
 54            MULT = DET(ITER1)
               IF (ITER.NE.2)                                   GOTO 55
                                                                GOTO 57
C              THEN
 55               DO 56 LL=1,ITER2
                     L = ITER1 - LL
 56                  MULT = MOD(MULT*PRIME(L) + DET(L),P)
C                 OD
 57            MULT = MOD(IP * MOD(DMOD-MULT,P), P)
 58         DET(ITER) = MULT - (MULT/P2)*P
            IF (DET(ITER).NE.0) TEST = .FALSE.
C                                                              *
C        END MIXED-RADIX ***************************************
C
         IF (TEST)                                              GOTO 59
                                                                GOTO 60
C        THEN
 59         NOZERO = NOZERO + 1
                                                                GOTO 61
C        ELSE
 60         NOZERO = 0
 61      IF (RTEST .AND. ITER.GT.NOZERO .AND. NOZERO.GE.TR)     GOTO 62
                                                                GOTO 64
C        THEN ... NORMAL EXIT; FIRST, HOWEVER CHECK FOR ZERO
C             ... DETERMINANT.
 62         IER = 0
            NOCOEF = ITER - NOZERO
            DO 63 K=1,NOCOEF
               IF (DET(K).NE.0) RETURN
C           OD
 63                                                             CONTINUE
            IER = 2
            RETURN
 64   CONTINUE
C
C     **MAXPRM COEFFICIENTS HAVE BEEN USED, ADNORMAL EXIT.
C
      NOCOEF = MAXPRM - NOZERO
      IER = 1
      IF (NOCOEF.EQ.0) IER = 3
      RETURN
      END
      BLOCK DATA                                                        DATA0010
C                                                                       DATA0020
C***************************************************************        DATA0030
C                                                              *        DATA0040
C                  *****WARNING*****                           *        DATA0050
C                                                              *        DATA0060
C ARRAY NAMES ARE USED IN THE DATA STATEMENTS BELOW TO         *        DATA0070
C SPECIFY VALUES FOR THE ARRAYS KPRIME AND IPRIME. IN SOME     *        DATA0080
C INSTALLATION IT MAY BE NECESSARY TO EXPLICITLY LIST THE      *        DATA0090
C COMPONENTS OF THE ARRAYS, NAMELY, KPRIMEI1 , KPRIMEI2 ,...,  *        DATA0100
C KPRIMEI100  AND IPRIMEI1 , IPRIMEI2 ,...,IPRIMEI100 .        *        DATA0110
C                                                              *        DATA0120
C***************************************************************        DATA0130
C                                                                       DATA0140
      COMMON /PRIMEB/KPRIME(100),IPRIME(100)                            DATA0150
      DATA KPRIME/       45233,45247,45259,45263,45281,45289            DATA0160
     1,45293,45307,45317,45319,45329,45337,45341,45343,45361            DATA0170
     2,45377,45389,45403,45413,45427,45433,45439,45481,45491            DATA0180
     3,45497,45503,45523,45533,45541,45553,45557,45569,45587            DATA0190
     4,45589,45599,45613,45631,45641,45659,45667,45673,45677            DATA0200
     5,45691,45697,45707,45737,45751,45757,45763,45767,45779            DATA0210
     6,45817,45821,45823,45827,45833,45841,45853,45863,45869            DATA0220
     7,45887,45893,45943,45949,45953,45959,45971,45979,45989            DATA0230
     8,46021,46027,46049,46051,46061,46073,46091,46093,46099            DATA0240
     9,46103,46133,46141,46147,46153,46171,46181,46183,46187            DATA0250
     X,46199,46219,46229,46237,46261,46271,46273,46279,46301            DATA0260
     Y,46307,46309,46327,46337/                                         DATA0270
      DATA IPRIME/       00000,42015,28577,01108,29342,16641            DATA0280
     1,10405,19447,26685,39525,14116,12753,32178,01043,08857            DATA0290
     2,27911,15049,07079,33425,00804,23175,23886,44779,41942            DATA0300
     3,10171,16606,10638,17371,27195,35827,42639,01829,24658            DATA0310
     4,09023,37958,30638,06339,41270,40538,10157,11783,00457            DATA0320
     5,32947,42170,17910,33474,20017,25086,36508,37444,35543            DATA0330
     6,06993,10326,16328,26765,42083,37223,30711,09408,06635            DATA0340
     7,38421,11397,32683,17333,34245,15748,35735,23492,19302            DATA0350
     8,20076,45620,44978,09864,14832,16092,19457,24045,44950            DATA0360
     9,32872,24309,15726,43057,37766,14046,41826,19946,41363            DATA0370
     X,23967,39791,29237,18085,12952,36850,02213,30023,34871            DATA0380
     Y,42667,40410,32615,46136/                                         DATA0390
      END                                                               DATA0400
      SUBROUTINE SUBBND(A,B,GAMMA,DELTA,S,SUM,T,N,M,NM,                 SBND0010
     1                  TWOT1,BOUND,NO)
C
C***************************************************************
C                                                              *
C USING HADAMARD'S INEQUALITY, THIS SUBROUTINE COMPUTES THE    *
C MAXIMUM NUMBER OF PRIMES REQUIRED TO REPRESENT DETERMINANT(A)*
C AND A(ADJOINT)*B FOR THE SYSTEM OF LINEAR EQUATIONS AX=B,    *
C WHERE A( ,N,N) AND B( ,N,M) ARE MATRICES WITH MULTIPLE-      *
C PRECISION INTEGER COEFFICIENTS EXPRESSED IN MIXED-RADIX FORM.*
C                                                              *
C  PRIME  IS A LINEAR ARRAY CONTAINING 100 DISTINCT PRIME      *
C (INPUT) INTEGERS IN ASCENDING ORDER.  THESE PRIMES ARE THE   *
C         RADII USED IN THE REPRESENTATION OF A AND B AND OF   *
C         DETERMINANT(A) AND A(ADJOINT)*B.THE PRIMES ARE CHOSEN*
C         AS LARGE AS POSSIBLE SUBJECT TO THE CONDITION THAT   *
C         FOR ALL I AND J PRIME(I)*PRIME(J) DOES NOT OVERFLOW  *
C         AN INTEGER WORD.                                     *
C      A  IS AN INTEGER MATRIX OF DIMENSION T BY N BY N.  THE  *
C (INPUT) FIRST DIMENSION CONTAINS THE COEFFICIENTS OF THE     *
C         MIXED-RADIX REPRESENTATION OF THE MULTIPLE-PRECISION *
C         COMPONENTS OF THE N BY N MATRIX A( ,N,N).  THAT IS,  *
C            A( ,I,J) = A(1,I,J)                               *
C                     + A(2,I,J)*PRIME(1)                      *
C                     .                                        *
C                     .                                        *
C                     .                                        *
C                     + A(T,I,J)*PRIME(1)*...*PRIME(T-1),      *
C         FOR I,J=1,2,...,N.                                   *
C      B  IS AN INTEGER MATRIX OF DIMENSION T BY N BY M WITH   *
C (INPUT) A SIMILAR NOTATIONAL CORRESPONDENCE MADE FOR A ABOVE.*
C  DELTA  IS AN INTEGER ARRAY OF DIMENSION TWOT1 USED TO STORE *
C  (TEMP) THE COEFFICIENTS OF THE MULTIPLE-PRECISION INTEGERS  *
C         A( ,I,J)**2 AND B( ,I,J)**2 FOR ANY PARTICULAR       *
C         A( ,I,J) OR B( ,I,J).                                *
C  GAMMA  IS AN INTEGER ARRAY OF DIMENSION TWOT1 USED TO STORE *
C  (TEMP) THE COEFFICIENTS OF THE MULTIPLE-PRECISION PARTIAL   *
C         PRODUCTS OBTAINED WHEN A( ,I,J)**2 OR B( ,I,J)**2 ARE*
C         BEING COMPUTED.                                      *
C    SUM  IS AN INTEGER ARRAY OF DIMENSION TWOT1 USED TO STORE *
C  (TEMP) THE SUM OF (A( ,I,J)**2, I=1,...,N) OR THE SUM OF    *
C         (B( ,I,J)**2 ,I=1,...,N) FOR ANY PARTICULAR COLUMN J.*
C      S  IS A REAL ARRAY OF DIMENSION NM. FOR 1<=J<=N, S(J)   *
C  (TEMP) CONTAINS A BOUND FOR THE LOGARITHM OF SUM(A( ,I,J)**2*
C         ,I=1,...,N) AND FOR N+1<=J<=NM OF SUM(B( ,I,J)**2,   *
C         I=1,...,N).                                          *
C      T  IS THE NUMBER OF RADII, PRIME(1),...,PRIME(T), USED  *
C (INPUT) TO REPRESENT EACH COMPONENT OF A( ,N,N) AND B( ,N,M) *
C         IN MIXED-RADIX FORM.                                 *
C      N  IS THE NUMBER OF EQUATIONS AND THE NUMBER OF UNKNOWNS*
C (INPUT) IN THE SYSTEM.  (I.E., N IS THE SIZE OF THE SECOND   *
C         AND THIRD DIMENSIONS OF A.)                          *
C      M  IS THE NUMBER OF RIGHT-HAND VECTORS FOR WHICH THE    *
C (INPUT) SYSTEM IS TO BE SOLVED.  (I.E., M IS THE SIZE OF THE *
C         THIRD DIMENSION OF B.)                               *
C     NM  IS EQUAL TO N+M.                                     *
C (INPUT)                                                      *
C  TWOT1  IS EQUAL TO 2*T+1.                                   *
C (INPUT)                                                      *
C    NO   IS A BOUND FOR THE NUMBER OF PRIMES REQUIRED TO SOLVE*
C(OUTPUT) THE SYSTEM OF LINEAR EQUATIONS AX=B. THAT IS, NO IS  *
C         A BOUND FOR THE NUMBER OF PRIMES REQUIRED TO         *
C         REPRESENT DETERMINANT(A) AND A(ADJOINT)*B IN MIXED-  *
C         RADIX FORM.                                          *
C   BOUND IS EQUAL TO THE LOGARITHM OF THE PRODUCT OF THE FIRST*
C(OUTPUT) NO PRIMES.                                           *
C                                                              *
C***************************************************************
C
      INTEGER T,TWOT,TWOT1,A(T,N,N),B(T,N,M),DELTA(TWOT1),
     1        GAMMA(TWOT1),PRIME(100),SUM(TWOT1),CSUM,Q,PD
      REAL S(NM)
      INTEGER  TK,TK1
      COMMON /PRIMEB/ PRIME, IPRIME(100)
      TWOT = 2*T
C     PROCEDURE ABOUND******************************************
C        FOR EACH J, 1<=J<=N, THE SUM OF THE SQUARES OF THE    *
C        ELEMENTS IN THE  J-TH COLUMN OF A IS COMPUTED. THE    *
C        SUM, WHICH IS A MIXED-RADIX, MULTIPLE-PRECISION       *
C        INTEGER OF LENGTH AT MOST 2*T+1, IS STORED IN THE     *
C        VECTOR SUM. THE LOGARITHM OF A BOUND FOR THE SUM IS   *
C        THEN COMPUTED AND STORED IN S(J).                     *
C                                                              *
         DO 16 J=1,N
C           **COMPUTE LOG (BOUND FOR (SUM OF A( ,I,J)**2,
C           **I=1,...,N)), AND STORE THE RESULTS IN S(J).
            DO 1  K=1,TWOT1
 1             SUM(K) = 0
C           OD
            DO 9 I=1,N
C              PROCEDURE SQUARE*********************************
C                 A( ,I,J)**2 IS COMPUTED. THE COEFFICIENTS OF *
C                 THE MIXED-RADIX, MULTIPLE-PRECISION RESULT   *
C                 ARE STORED IN DELTA.                         *
C                                                              *
                  DO 2  II=1,TWOT1
                     DELTA(II) = 0
 2                   GAMMA(II) = 0
C                 OD
                  Q = 0
                  K = 0
                  L = T-1
C                 REPEAT
 3                   IF (K.LE.L)                                GOTO 4
                                                                GOTO 7
C                    THEN
C                       **COMPUTE PARTIAL PRODUCT AND STORE IN
C                       **GAMMA.
 4                      Q = 0
                        TK = T-K
                        DO 5  II=1,T
                           PD = A(TK,I,J)*A(II,I,J) + Q
                           Q = PD/PRIME(II)
 5                         GAMMA(II) = PD - Q*PRIME(II)
C                       OD
                        GAMMA(T+1) = Q
C                       **ADD PARTIAL PRODUCT TO DELTA.
                        Q = 0
                        L1 = T+K
                        TK1 = T+K+1
                        DO 6  II=1,L1
                           PD = PRIME(TK)*DELTA(II) + GAMMA(II)
     1                          + Q
                           Q = PD/PRIME(II)
 6                         DELTA(II) = PD - Q*PRIME(II)
C                       OD
                        DELTA(TK1) = Q
                        K = K+1
                                                                GOTO 3
C                 CONTINUE
C                                                              *
C              END SQUARE***************************************
C
C
C              **ACCUMULATE THE SUM OF A( ,I,J)**2.
 7             Q = 0
               DO 8 II=1,TWOT1
                  CSUM = DELTA(II)+SUM(II)+Q
                  Q = CSUM/PRIME(II)
 8                SUM(II) = CSUM - Q*PRIME(II)
C              OD
C           OD
 9                                                              CONTINUE
C
C
C           **CALCULATE ABOUND FOR THE MULTIPLE-PRECISION
C           **INTEGER IN SUM AND STORE THE LOG OF THIS
C           **BOUND IN S(J).
            K = TWOT1
C           REPEAT...FIND OUT THE LENGTH K OF THE MULTIPLE-
C                 ...PRECISION NUMBER IN SUM.
 10            IF (K.GT.1.AND.SUM(K).EQ.0)                      GOTO 11
                                                                GOTO 12
C              THEN
 11               K = K - 1
                                                                GOTO 10
C           CONTINUE
 12         S(J) = 0.0
            IF (K.GT.1)                                         GOTO 13
                                                                GOTO 15
C           THEN
 13            L = K-1
               DO 14 II=1,L
 14               S(J) = S(J) +ALOG(FLOAT(PRIME(II)))
C              OD
               S(J) = S(J)+ALOG(FLOAT(SUM(K)+1))
                                                                GOTO 16
C           ELSE
 15            IF (SUM(K).GT.0) S(J) = ALOG(FLOAT(SUM(K)))
C        OD
 16                                                             CONTINUE
C                                                              *
C     END ABOUND************************************************
C
C
C     PROCEDURE BBOUND******************************************
C        FOR EACH J, 1<=J<=M, THE SUM OF THE SQUARES OF THE    *
C        ELEMENTS IN THE  J-TH COLUMN OF B IS COMPUTED. THE    *
C        SUM, WHICH IS A MIXED-RADIX, MULTIPLE PRECISION       *
C        INTEGER OF LENGTH AT MOST 2*T+1, IS STORED IN THE     *
C        VECTOR SUM. THE LOGARITHM OF A BOUND FOR THE SUM IS   *
C        THEN COMPUTED AND STORED IN S(N+J).                   *
C                                                              *
C        NOTE: THE COMPUTATIONAL PROCEDURE IS THE SAME AS FOR  *
C              A ABOVE.                                        *
C                                                              *
         DO 32 J=1,M
            DO 17 K=1,TWOT1
 17            SUM(K) = 0
C           OD
            DO 25 I=1,N
               DO 18 II=1,TWOT1
                  DELTA(II) = 0
 18               GAMMA(II) = 0
C              OD
               Q = 0
               K = 0
               L = T-1
C              REPEAT
 19               IF (K.LE.L)                                   GOTO 20
                                                                GOTO 23
C                 THEN
 20                  Q = 0
                     TK = T-K
                     DO 21 II=1,T
                        PD = B(TK,I,J)*B(II,I,J) + Q
                        Q = PD/PRIME(II)
 21                     GAMMA(II) = PD - Q*PRIME(II)
C                    OD
                     GAMMA(T+1) = Q
                     Q = 0
                     L1 = T+K
                     TK1 = T+K+1
                     DO 22 II=1,L1
                        PD = PRIME(TK)*DELTA(II) + GAMMA(II)+Q
                        Q = PD/PRIME(II)
 22                     DELTA(II) = PD - Q*PRIME(II)
C                    OD
                     DELTA(TK1) = Q
                     K = K+1
                                                                GOTO 19
C              CONTINUE
 23            Q = 0
               DO 24 II=1,TWOT1
                  CSUM = DELTA(II)+SUM(II)+Q
                  Q = CSUM/PRIME(II)
 24               SUM(II) = CSUM - Q*PRIME(II)
C              OD
C           OD
 25                                                             CONTINUE
            K = TWOT1
C           REPEAT
 26            IF (K.GT.1.AND.SUM(K).EQ.0)                      GOTO 27
                                                                GOTO 28
C              THEN
 27               K = K - 1
                                                                GOTO 26
C           CONTINUE
 28         NJ = N+J
            S(NJ) = 0.0
            IF (K.GT.1)                                         GOTO 29
                                                                GOTO 31
C           THEN
 29            L = K-1
               DO 30 II=1,L
 30               S(NJ) = S(NJ) +ALOG(FLOAT(PRIME(II)))
C              OD
               S(NJ) = S(NJ)+ALOG(FLOAT(SUM(K)+1))
                                                                GOTO 32
C           ELSE
 31            IF (SUM(K).GT.0) S(J) = ALOG(FLOAT(SUM(K)))
C        OD
 32                                                             CONTINUE
C                                                              *
C     END BBOUND************************************************
C
C
C     **FIND MIN(S(I)), I=1,...,N.
         SMIN = S(1)
         DO 33 I=1,N
 33         IF (SMIN.GT.S(I)) SMIN = S(I)
C        OD
C     **FIND MAX(S(I)), I=N+1,...,NM.
         SMAX = S(N+1)
         L = N+1
         DO 34 I=L,NM
 34         IF (SMAX.LT.S(I)) SMAX = S(I)
C        OD
C     **IF SMIN >= SMAX, THEN BOUND = (SUM OF S(I))/2+ LOG(2.0)
C     **ELSE BOUND = ((SUM OF S(I))+SMAX-SMIN)/2+LOG(2.0),
C     **WHERE THE SUM IS TAKEN OVER I=1,...,N.
         BOUND = 0.0
         DO 35 I=1,N
 35         BOUND = BOUND+S(I)
C        OD
         IF (SMIN.GE.SMAX)                                      GOTO 36
                                                                GOTO 37
C        THEN
 36         BOUND = BOUND/2.0 + ALOG(2.0)
                                                                GOTO 38
C        ELSE
 37         BOUND = (BOUND+SMAX-SMIN)/2.0 + ALOG(2.0)
C     **CALCULATE NUMBER OF PRIMES REQUIRED.
 38      NO = 0
         SUMLOG = 0.0
C        REPEAT
 39         NO = NO+1
            IF (SUMLOG.GE.BOUND)                                GOTO 40
C           THEN
C       <------EXIT
            SUMLOG = SUMLOG + ALOG(FLOAT(PRIME(NO)))
                                                                GOTO 39
C        CONTINUE
 40   RETURN
      END
      SUBROUTINE MRADIX(A,N,C,LMAX,L,B,IER)                             MRAD0010
C
C***************************************************************
C GIVEN THE FIXED-RADIX INTEGER                                *
C                                                              *
C   A(1)*B**(N-1) + ... + A(N-1)*B + A(N),                     *
C                                                              *
C WHERE ABS(A(I)) < B, I=1,...,N,                              *
C THIS SUBROUTINE COMPUTES(FOR SOME L <= LMAX) THE COEFFICIENTS*
C C(1),C(2),...,C(L) IN ITS SYMMETRIC MIXED-RADIX              *
C REPRESENTATION:                                              *
C                                                              *
C   C(1)+C(2)*P(1) + ... + C(L)*P(1)*P(2)*...*P(L-1).          *
C                                                              *
C LMAX MUST BE GIVEN SUCH THAT                                 *
C                                                              *
C   LMAX >= (N * LOG B + LOG 2) / LOG(P(1)).                   *
C                                                              *
C IER IS AN ERROR CODE WHICH IS 1 IF THE DIMENSION PARAMETER N,*
C LMAX ARE INCORRECT, AND 0 OTHERWISE.                         *
C                                                              *
C                   ******* WARNING  *******                   *
C                                                              *
C THIS SUBROUTINE ASSUMMES THAT FOR ALL I,J                    *
C                                                              *
C   P(I)*P(J) AND P(I)*B                                       *
C                                                              *
C DO NOT OVERFLOW A COMPUTER WORD.                             *
C                                                              *
C***************************************************************
C
      INTEGER A(N),C(LMAX),P(100),Q,QTEMP,PP,B,P1
      COMMON /PRIMEB/ P, IPRIME(100)
      IER = 1
      IF (N.LT.1 .OR. LMAX.LT.1) RETURN
      L = 1
      C(1) = 0
      DO 4 I=1,N
C        **AT THIS STAGE C(1)+C(2)*P(1)+...+C(L)*P(1)*P(2)*...*
C        **P(L-1) IS THE MIXED-RADIX REPRESENTATION OF
C        **A(1)*B**(I-2)+A(2)*B**(I-3)+...+A(I-1).
         Q = A(I)
         DO 1 J=1,L
C           **COMPUTE THE FIRST L COEFFICIENTS OF THE MIXED-
C           **RADIX REPRESENTATION OF A(I)+B*(C(1)+C(2)*P(1)+
C           **...+C(L)*P(1)*P(2)*...*P(L-1)).
            PP = P(J)
            QTEMP = Q + B*C(J)
            Q = QTEMP/PP
 1          C(J) = QTEMP - Q*PP
C        OD
C        REPEAT ... CONVERT C(1)+C(2)*P(1)+...+C(L)*P(1)*
C               ... P(2)*...*P(L-1)+Q*P(1)*...*P(L)
C               ... TO MIXED-RADIX FORM.
 2          IF (Q .NE.0)                                        GOTO 3
                                                                GOTO 4
C           THEN
 3             L = L+1
               IF (L.GT.LMAX) RETURN
               PP = P(L)
               QTEMP = Q/PP
               C(L) = Q - QTEMP*PP
               Q = QTEMP
                                                                GOTO 2
C        CONTINUE
C     OD
 4                                                              CONTINUE
C     **CONVERT TO SYMMETRIC MIXED-RADIX FORM.
        Q = 0
        DO 5 I=1,L
           C(I) = C(I)+Q
           P1 = (P(I)+1)/2
           Q = C(I)/P1
 5         C(I) = C(I) - Q*P(I)
C       OD
        IF (Q.NE.0)                                             GOTO 6
                                                                GOTO 7
C       THEN
 6         L = L+1
           IF (L.GT.LMAX) RETURN
           C(L) = Q
 7      IER = 0
      RETURN
      END
      SUBROUTINE FRADIX(C,N,A,LMAX,L,B,IER)                             FRAD0010
C
C***************************************************************
C GIVEN THE MIXED-RADIX INTEGER                                *
C                                                              *
C   C(1)+C(2)*P(1)+ ... +C(N)*P(1)*...*P(N-1),                 *
C                                                              *
C WHERE ABS(C(I)) < (P(I)+1)/2, I=1,...,N,                     *
C THIS SUBROUTINE COMPUTES(FOR SOME L <= LMAX) THE COEFFICIENTS*
C A(1),A(2),...,A(L) IN ITS FIXED-RADIX REPRESENTATION:        *
C                                                              *
C   A(1)*B**(L-1)+A(2)*B*(L-2)+...+A(L),                       *
C                                                              *
C WHERE ABS(A(I)) < B.                                         *
C                                                              *
C LMAX MUST BE GIVEN SUCH THAT                                 *
C                                                              *
C   LMAX >= (N*LOG(P(N)) - LOG 2) / LOG B.                     *
C                                                              *
C IER IS AN ERROR CODE WHICH IS 1 IF THE DIMENSION PARAMETER N,*
C LMAX ARE INCORRECT, AND 0 OTHERWISE.                         *
C                                                              *
C                   ******* WARNING  *******                   *
C                                                              *
C THIS SUBROUTINE ASSUMMES THAT FOR ALL I                      *
C                                                              *
C    P(I)*B                                                    *
C                                                              *
C DOES NOT OVERFLOW A COMPUTER WORD.                           *
C                                                              *
C***************************************************************
C
      INTEGER C(N),A(LMAX),P(100),Q,QTEMP,PP,B
      COMMON /PRIMEB/ P, IPRIME(100)
      IER = 1
      IF (N.LT.1 .OR. LMAX.LT.1) RETURN
      L = 0
      DO 5 I=1,N
C        **AT THIS STAGE A(1)+A(2)*B+...+A(L)*B**(L-1) IS THE
C        **FIXED-RADIX REPRESENTATION OF C(N-I+2)+C(N-I+3)*
C        **P(N-I+2)+...+C(N)*P(N-I+2)*...*P(N-1).
         NI1 = N-I+1
         PP = P(NI1)
         Q = C(NI1)
         IF (L.GT.0)                                            GOTO 1
                                                                GOTO 3
C        THEN ... COMPUTE THE FIRST L COEFFICIENTS OF THE
C             ... FIXED-RADIX REPRESENTATION OF C(N-I+1)+
C             ... P(N-I+1)*(A(1)+A(2)*B+...+A(L)*B**(L-1)).
 1          DO 2 J=1,L
               QTEMP = A(J)*PP + Q
               Q = QTEMP/B
 2             A(J) = QTEMP - Q*B
C           OD
C        REPEAT ... CONVERT A(1)+A(2)*B+...+A(L)*B**(L-1)
C               ... +Q*B**L TO FIXED-RADIX FORM.
 3          IF (Q.NE.0)                                         GOTO 4
                                                                GOTO 5
C           THEN
 4             L = L+1
               IF (L.GT.LMAX) RETURN
               QTEMP = Q/B
               A(L) = Q - QTEMP*B
               Q = QTEMP
                                                                GOTO 3
C        CONTINUE
 5                                                              CONTINUE
C     OD
      IF (L.GE.2)                                               GOTO 6
                                                                GOTO 8
C     THEN ... REORDER THE COEFFICIENTS.
 6       L2 = L/2
         DO 7 I=1,L2
            LI1 = L-I+1
            QTEMP = A(I)
            A(I) = A(LI1)
 7          A(LI1) = QTEMP
C        OD
 8    IER = 0
      RETURN
      END
