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
.