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= 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