#! /bin/sh
# This is a shell archive.  Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file".  To overwrite existing
# files, type "sh file -c".  You can also feed this as standard input via
# unshar, or by typing "sh <file", e.g..  If this archive is complete, you
# will see the following message at the end:
#		"End of shell archive."
# Contents:  README discrete.f gfssrt.f lattest.f mul.f normgam.f
#   proot.f
# Wrapped by ripley@puffin.stats on Sun Aug  9 14:06:02 1992
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f 'README' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'README'\"
else
echo shar: Extracting \"'README'\" \(115 characters\)
sed "s/^X//" >'README' <<'END_OF_FILE'
XThis directory contains the Fortran programs from Appendix B of 
X
XB.D. Ripley (1987) `Stochastic Simulation' Wiley
END_OF_FILE
if test 115 -ne `wc -c <'README'`; then
    echo shar: \"'README'\" unpacked with wrong size!
fi
# end of 'README'
fi
if test -f 'discrete.f' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'discrete.f'\"
else
echo shar: Extracting \"'discrete.f'\" \(1285 characters\)
sed "s/^X//" >'discrete.f' <<'END_OF_FILE'
X      SUBROUTINE SETIND (P,M,IND,MI)
X      REAL P(M),P0
X      INTEGER IND(MI),M,MI,I,K
X      I = 1
X      DO 20 K = 1,MI
X        P0 = REAL(K-1)/REAL(MI)
X   10   IF (P(I) .GE. P0) GO TO 20
X        I = I+1
X        GO TO 10
X   20   IND(K) = I
X      RETURN
X      END
X
X      FUNCTION DISRV(P,M,IND,MI)
X      REAL P(M),U
X      INTEGER DISRV, IND(MI),M,MI,I
X      U = RND()
X      I = IND(INT(MI*U)+1)
X   10 IF (P(I) .GE. U) GO TO 20
X      I = I+1
X      GO TO 10
X   20 DISRV = I
X      RETURN
X      END
X
X      SUBROUTINE SETAL (P,M,A,Q,W)
X      REAL P(M),Q(M)
X      INTEGER M,A(M),W(M),I,J,NN,NP,S
X      NN = 0
X      NP = M+1
X      DO 10 I = 1,M
X        Q(I) = M*P(I)
X        IF (Q(I) .LT. 1.0) THEN
X          NN = NN+1
X          W(NN) = I
X        ELSE
X          NP = NP-1
X          W(NP) = I
X        ENDIF
X   10   CONTINUE
X      DO 20 S = 1,M-1
X        I = W(S)
X        J = W(NP)
X        A(I) = J
X        Q(J) = Q(J)+Q(I)-1.0
X        IF (Q(J) .LT. 1.0) NP = NP+1
X   20   CONTINUE
X       A(W(M)) = W(M)
X      DO 30 I = 1,M
X   30   Q(I) = Q(I)+I-1
X      RETURN
X      END 
X
X      FUNCTION ALRV (A,Q,M)
X      INTEGER ALRV,A(M),M,I
X      REAL Q(M),U
X      U = M*RND()
X      I = 1+INT(U)
X      IF ( U .LE. Q(I)) THEN
X        ALRV = I
X      ELSE
X        ALRV = A(I)
X      ENDIF
X      RETURN
X      END
END_OF_FILE
if test 1285 -ne `wc -c <'discrete.f'`; then
    echo shar: \"'discrete.f'\" unpacked with wrong size!
fi
# end of 'discrete.f'
fi
if test -f 'gfssrt.f' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'gfssrt.f'\"
else
echo shar: Extracting \"'gfssrt.f'\" \(1360 characters\)
sed "s/^X//" >'gfssrt.f' <<'END_OF_FILE'
X      SUBROUTINE GFSRT (A,P,N)
XC     set pmax as required
X      INTEGER P,N,A(pmax,pmax),EOR,I,J,L,T
X      DO 60 I = 1,N
X        J = I
X        IF (A(I,I) .EQ. 0) THEN
X   10     J = J+1
X          IF (A(J,I) .GT. 0) GO TO 20
X          IF (J .LT. P) GO TO 10
X          PRINT *,'RANK DEFICIENT'
X          RETURN
X   20     DO 30 L = I,N
X            T = A(J,L)
X            A(J,L) = A(I,L)
X   30       A(I,L) = T
X        ENDIF
X          DO 50 J = I+1,P
X            IF (A(J,I) .GT. 0) THEN
X              DO 40 L = I,N
X   40           A(J,L) = EOR(A(J,L),A(I,L))
X            ENDIF
X   50     CONTINUE
X   60   CONTINUE
X      PRINT *,'FULL RANK'
X      RETURN
X      END
X
X      FUNCTION EOR(I,J)
X      INTEGER EOR,I,J
X      EOR = MOD(I+J,2)
X      RETURN
X      END
X
X      SUBROUTINE GFSRS (Y,L,P,K)
XC     set pmax as required
X      INTEGER L,P,K,Y(P+K-1),A(pmax,pmax),D
X      N = K*L
X      IF (N .GT. P) THEN
X        PRINT *,'MUST FAIL'
X        RETURN
X      ENDIF
X      M = 2**(L-1)
X      DO 30 I = 1,P
X        I1 = 0
X        DO 20 J = 1,K
X          D = Y(I+J-1)
X          M1 = M
X          DO 10 J1 = 1,L
X            I1 = I1+1
X            IF (D .GE. M1) THEN
X              A(I,I1) = 1
X              D = D-M1
X            ELSE
X              A(I,I1) = 0
X            ENDIF
X   10       M1 = M1/2
X   20     CONTINUE
X   30   CONTINUE
X      CALL GFSRT (A,P,N)
X      RETURN
X      END
END_OF_FILE
if test 1360 -ne `wc -c <'gfssrt.f'`; then
    echo shar: \"'gfssrt.f'\" unpacked with wrong size!
fi
# end of 'gfssrt.f'
fi
if test -f 'lattest.f' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'lattest.f'\"
else
echo shar: Extracting \"'lattest.f'\" \(6743 characters\)
sed "s/^X//" >'lattest.f' <<'END_OF_FILE'
X      SUBROUTINE LATT (R,MULT,M)
X      LOGICAL TEST2,TEST1
X      INTEGER R,K
X      DOUBLE PRECISION X(8,8),LEN(8),MULT,M
X      DO 20 K = 2,R
X        CALL INIT (X,LEN,K,MULT,M)
X   10   IF (TEST1(X,LEN,K)) GO TO 10
X        IF (TEST2(X,LEN,K)) GO TO 10
X        CALL RES (X,LEN,K,M)
X   20   CONTINUE
X      RETURN
X      END
XC
X      SUBROUTINE INIT (X,LEN,R,MULT,M)
X      INTEGER R,I,K
X      DOUBLE PRECISION X(8,8),LEN(8),L,MULT,M,A,MUL
X      IF (R .EQ. 2) THEN
X      L = 1.0
X      DO 10 K = 1,8
X        X(1,K) = L
X   10   L = MUL(L,MULT,M)
X      DO 20 I = 2,8
X        DO 20 K = 1,8
X          X(I,K) = 0.0
X          IF (I .EQ. K) X(I,K) = M
X   20     CONTINUE
X      ELSE
X        A = X(1,R)
X        DO 25 I = 1,R-1
X   25     X(I,R) = MUL(A,X(I,1),M)
X      ENDIF
X      DO 40 I = 1,R
X        A = 0.0
X        DO 30 K = 1,R
X   30     A = A+X(I,K)*X(I,K)
X   40   LEN(I) = A
X      RETURN
X      END
XC
X      FUNCTION TEST1 (X,LEN,R)
X      LOGICAL TEST1
X      INTEGER R,I,J,K,I1,I2,L,NCHNGS
X      DOUBLE PRECISION X(8,8),LEN(8),XY,A,B,MX
X      DATA MX/maxint/
X      NCHNGS = 0
XC     try each pair I,J in turn    
X      DO 40 I = 2,R
X        DO 40 J = 1,I-1
X   10     XY = 0.0
X          DO 20 K = 1,R
X   20       XY = XY+X(I,K)*X(J,K)
X          IF (LEN(I) .LE. LEN(J)) THEN
X            I1 = I
X            I2 = J
X          ELSE
X            I1 = J
X            I2 = I
X          ENDIF
XC     round halves towards zero          
X          A = XY/LEN(I1)
X          L = INT(ABS(A)+0.499999999)
X          IF (A .LT. 0.0) L = -L
X          IF (L .EQ. 0) GO TO 40
X          NCHNGS = NCHNGS+1
X          A = 0.0
X          DO 30 K = 1,R
X            B = L*X(I1,K)
X            IF (ABS(B) .GT. MX) PRINT *,'Accuracy loss'
X            B = X(I2,K)-B
X            X(I2,K) = B
X   30       A = A+B*B
X          LEN(I2) = A
X          IF (LEN(I2) .LT. LEN(I1)) GO TO 10
X   40     CONTINUE
X      TEST1 = NCHNGS .GT. 0
X      RETURN
X      END
XC
X      FUNCTION TEST2 (X,LEN,R)
X      LOGICAL TEST2
X      INTEGER I,I1,J,K,L,L1,R,S,S1,T1,CON(3),IN(8),PTR(8)
X      DOUBLE PRECISION X(8,8),LEN(8),T(8),A,AL,B,MX
X      DATA CON/0,-1,+1/
X      DATA MX/maxint/
X      TEST2 = .FALSE.
X      IF (R .EQ. 2) RETURN 
X      DO 10 I = 1,R
X   10   PTR(I) = I
XC     PTR is pointer to vectors in length order
X      DO 30 I = 1,R-1
X        L = I
X        AL = LEN(PTR(L))
X        DO 20 J = I+1,R
X        L1 = PTR(J)
X          IF (LEN(L1) .LT. AL) THEN
X            L = J
X            AL = LEN(L1)
X          ENDIF
X   20     CONTINUE
X        IF (L .NE. I) THEN
X          L1 = PTR(L)
X          PTR(L) = PTR(I)
X          PTR(I) = L1
X        ENDIF
X   30   CONTINUE
XC     try Minkowski's test
X      DO 100 I = 3,R
X        I1 = PTR(I)
X        AL = LEN(I1)-0.5
X        IN(I) = 1
X        DO 90 S = 1,3**(I-1)-1
X          S1 = S
X          DO 50 L = 1,I-1
X            T1= S1/3
X            IN(L) = CON(S1-3*T1+1)
X   50       S1 = T1
X          A = 0.0
X          DO 70 J = 1,R
X            B = 0.0
X            DO 60 K = 1,I
X              B = B+X(PTR(K),J)*IN(K)
X              IF (ABS(B) .GT. MX) PRINT *,'Possible accuracy loss'
X   60         CONTINUE
X            T(J) = B
X   70       A = A+B*B
X          IF (A .GT. AL) GO TO 90
X          LEN(PTR(I)) = A
X          DO 80 J = 1,R
X   80       X(I1,J) = T(J)
X          TEST2 = .TRUE.
X          PRINT *,'TEST2 SUCCESS'
X          RETURN
X   90     CONTINUE
X  100   CONTINUE
X      RETURN
X      END
XC
X      SUBROUTINE RES (X,LEN,R,M)
X      DOUBLE PRECISION X(8,8),LEN(8),A(8,8),M,LL,LU
X      REAL NU,RATIO,AA,B,UL
X      INTEGER R,Z(8),I,J
X      LL = LEN(1)
X      LU = LL
X      DO 10 I = 2,R
X        LL = MIN (LL,LEN(I))
X   10   LU = MAX (LU,LEN(I))
X      RATIO = SQRT(LU/LL)
X      NU = SQRT(LL)
X      IF (R .EQ. 2) GO TO 50
X      CALL INV (X,A,R,M)
X      B = 1.0E20
X      DO 30 I = 1,R
X        AA = 0.0
X        DO 20 J = 1,R
X   20     AA = AA+A(J,I)*A(J,I)
X   30   B = MIN(AA,B)
X      NU = SQRT(B)
X      DO 40 I = 1,R
X   40   Z(I) = INT(NU*SQRT(LEN(I))/M)
X      CALL SEARCH(A,R,NU,Z,M)
X   50 UL = SQRT(LU)/M
X      PRINT 1000,R,RATIO,UL,NU
X 1000 FORMAT (' DIM',I2,' RATIO ',1PG10.3,' LMAX ',1PE10.2,
X     & ' NU ',1PE10.2)
X      RETURN
X      END
XC
X      SUBROUTINE INV (X,Y,R,M)
XC     invert X/M to Y by Gaussian elimination 
XC     with partial pivoting
X      INTEGER R,H,I,J,K,N
X      DOUBLE PRECISION X(8,8),Y(8,8),M,W(8,16),S
X      N = R+R
X      DO 20 I = 1,R
X        DO 10 J = 1,R
X   10     W(I,J) = X(I,J)/M
X        DO 20 J = 1,R
X          IF (I .EQ. J) THEN
X            W(I,J+R) = 1.0
X          ELSE
X            W(I,J+R) = 0.0
X          ENDIF
X   20     CONTINUE
X      DO 60 J = 1,R-1
X        S = ABS(W(J,J))
X        K = J
X        DO 30 H = J+1,R
X          IF (ABS(W(H,J)) .GT. S) THEN
X            S = ABS(W(H,J))
X            K = H
X          ENDIF
X   30     CONTINUE
X        IF (K .NE. J) THEN
X          DO 40 I = J,N
X            S = W(K,I)
X            W(K,I) = W(J,I)
X   40       W(J,I) = S
X        ENDIF
X        DO 50 K = J+1,R
X          W(K,J) = W(K,J)/W(J,J)
X          DO 50 I = J+1,N
X   50       W(K,I) = W(K,I)-W(K,J)*W(J,I)
X   60   CONTINUE
X      DO 90 I = R+1,N
X        W(R,I) = W(R,I)/W(R,R)
X        DO 80 J = R-1,1,-1
X          S = W(J,I)
X          DO 70 K = J+1,R
X   70       S = S-W(J,K)*W(K,I)
X   80     W(J,I) = S/W(J,J)
X   90   CONTINUE
X      DO 100 I = 1,R
X        DO 100 J = 1,R
X  100     Y(I,J) = W(I,J+R)
X      RETURN
X      END
XC
X      SUBROUTINE SEARCH (A,R,NU,Z,M)
X      INTEGER R,Z(8),T(8),I,J,K
X      REAL NU,AC
X      DOUBLE PRECISION M, A(8,8), Y(8), AA
X      DO 10 I = 1,R
X        T(I) = 0
X   10   Y(I) = 0.0
X      K = R
X   20 IF (T(K) .EQ. Z(K)) GO TO 80
X      T(K) = T(K)+1
X      DO 30 J = 1,R
X   30   Y(J) = Y(J)+A(J,K)
X   40 K = K+1
X      IF (K .GT. R) GO TO 60
X      T(K) = -Z(K)
X      IF ( Z(K) .NE. 0) THEN
X        DO 50 J = 1,R
X   50     Y(J) = Y(J)-2*Z(K)*A(J,K)
X      ENDIF
X      GO TO 40
X   60 AA = 0.0
X      DO 70 J = 1,R
X   70   AA = AA+Y(J)*Y(J)
X      AC = SQRT(AA)
X      IF (AC .LT. 0.9999*NU) THEN
X        PRINT *,'SEARCH SUCCESS'
X        NU = AC
X      ENDIF
X   80 K = K-1
X      IF (K.GE.1) GO TO 20
X      RETURN
X      END
XC
X      FUNCTION MUL(A,B,C)
XC     forms A*B MOD C for A,B < C <= 2^50
X      DOUBLE PRECISION A,B,C,MUL,A1,A2,B1,B2,D,M,MM
XC     set M so 2*M*M <= maxint   
X      M = 2.0D0**25
X      MM = M*M
X      D = (A*B)/C
X      CALL SPLIT (A,B,A1,A2)
X      CALL SPLIT (C,D,B1,B2)
X      MUL = A1-B1 + (A2-B2)*MM
X      RETURN
X      END
XC
X      SUBROUTINE SPLIT (A,B,C,D)
X      DOUBLE PRECISION A,B,C,D,A1,A2,B1,B2,M,MM,AC,AD,C1,C2
X      M = 2**25
X      MM = M*M
X      A2 = INT(A/M)
X      A1 = NINT(A-A2*M)
X      B2 = INT(B/M)
X      B1 = NINT(B-B2*M)
X      AC = A2*B1 + A1*B2
X      C2 = INT(AC/M)
X      C1 = AC - M*C2
X      AC = C1*M + A1*B1
X      AD = INT(AC/MM)
X      C = AC - AD*MM
X      D = AD + C2 + A2*B2
X      RETURN
X      END
END_OF_FILE
if test 6743 -ne `wc -c <'lattest.f'`; then
    echo shar: \"'lattest.f'\" unpacked with wrong size!
fi
# end of 'lattest.f'
fi
if test -f 'mul.f' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'mul.f'\"
else
echo shar: Extracting \"'mul.f'\" \(811 characters\)
sed "s/^X//" >'mul.f' <<'END_OF_FILE'
X      FUNCTION MUL(A,B,C)
XC     forms A*B MOD C for A,B < C <= 2^50
X      DOUBLE PRECISION A,B,C,MUL,A1,A2,B1,B2,D,M,MM
XC     set M so that 2*M*M <= maxint
X      M = 2.0D0**25
X      MM = M*M
X      D = (A*B)/C
X      CALL SPLIT (A,B,A1,A2)
X      CALL SPLIT (C,D,B1,B2)
X      MUL = A1-B1 + (A2-B2)*MM
X      RETURN
X      END
XC
X      SUBROUTINE SPLIT (A,B,C,D)
XC     forms A*B = C + D*MM for 0<= C,D < MM
X      DOUBLE PRECISION A,B,C,D,A1,A2,B1,B2,M,MM,AC,AD,C1,C2
XC     set M to the same value as in MUL
X      M = 2.0D0**25
X      MM = M*M
X      A2 = dINT(A/M)
X      A1 = dINT(A-A2*M)
X      B2 = dINT(B/M)
X      B1 = dINT(B-B2*M)
X      AC = A2*B1 + A1*B2
X      C2 = dINT(AC/M)
X      C1 = AC - M*C2
X      AC = C1*M + A1*B1
X      AD = dINT(AC/MM)
X      C = AC - AD*MM
X      D = AD + C2 + A2*B2
X      RETURN
X      END
X
END_OF_FILE
if test 811 -ne `wc -c <'mul.f'`; then
    echo shar: \"'mul.f'\" unpacked with wrong size!
fi
# end of 'mul.f'
fi
if test -f 'normgam.f' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'normgam.f'\"
else
echo shar: Extracting \"'normgam.f'\" \(2881 characters\)
sed "s/^X//" >'normgam.f' <<'END_OF_FILE'
X      FUNCTION POLAR()
X      REAL AN,E,POLAR,V1,V2,W
X      INTEGER IR
X      SAVE IR,AN
X      DATA IR/0/
X      IF (IR.EQ.0) THEN
X   10   V1 = 2.*RND()-1.0
X        V2 = 2.*RND()-1.0
X        W = V1*V1 +V2*V2
X        IF (W.GT.1.) GO TO 10
X        E = SQRT((-2.0*LOG(W))/W)
X        AN = V1*E
X        IR = 1
X        POLAR = V2*E
X      ELSE
X        IR = 0
X        POLAR = AN
X      ENDIF
X      RETURN
X      END
X
X      FUNCTION NRATIO()
X      REAL NRATIO,U,V,X,Z
X   10 U = RND()
X      V = 0.8578*(2.*RND()-1.)
X      X = V/U
X      Z = 0.25*X*X
X      IF (Z .LT. 1.-U) GO TO 20
X      IF (Z .GT. (0.259/U+0.35)) GO TO 10
X      IF (Z .GT. -LOG(U)) GO TO 10
X   20 NRATIO = X
X      RETURN
X      END
X
X      FUNCTION NMB()
X      REAL NMB,AV,G,U,U1,U2,V,W
X      U = RND()
X      IF (U.LE.0.8638) THEN
X        NMB = 2.0*(RND()+RND()+RND()-1.5)
X        RETURN
X      ENDIF
X      IF (U.LE.0.9745) THEN
X        NMB = 1.5*(RND()+RND()-1.0)
X        RETURN
X      ENDIF
X      IF (U.LE.0.9973002039) THEN
X   10   V = 6.0*RND()-3.0
X        AV = ABS(V)
X        G = 17.49731196*EXP(-0.5*V*V)
X        IF (AV .LT. 1.0) THEN
X          G = G-4.73570326*(3.0-V*V)
X        ELSE
X          G = G-2.36785163*(3.0-AV)*(3.0-AV)
X        ENDIF
X        IF (AV .LT. 1.5) G = G-2.157875*(1.5-AV)
X        IF (0.358*RND() .GT. G) GO TO 10
X        NMB = V
X        RETURN
X      ENDIF
X   20 U1 = 2.0*RND()-1.0
X      U2 = 2.0*RND()-1.0
X      W = U1*U1+U2*U2
X      IF (W .GE. 1.0) GO TO 20
X      W = SQRT((9.0-2.0*LOG(W))/W)
X      IF (ABS(U1*W) .GT. 3.0) THEN
X        NMB = U1*W
X        RETURN
X      ENDIF
X      IF (ABS(U2*W) .LE. 3.0) GO TO 20
X      NMB = U2*W
X      RETURN
X      END
X
X      FUNCTION EXPRV()
X      REAL EXPRV,A,U,U0,USTAR
X      A = 0.0
X   10 U = RND()
X      U0 = U
X   20 USTAR = RND()
X      IF (U .LT. USTAR) GO TO 30
X      U = RND()
X      IF (U .LT. USTAR) GO TO 20
X      A = A+1.0
X      GO TO 10
X   30 EXPRV = A+U0
X      RETURN
X      END
X
X      FUNCTION GS (ALPHA)
X      REAL GS,ALPHA,B,P,X
X      DATA E/2.71828182/
X      B = (ALPHA+E)/E
X   10 P = B*RND()
X      IF (P.GT.1.0) GO TO 20
X      X = P**(1./ALPHA)
X      IF (X .GT. -LOG(RND())) GO TO 10
X      GS = X
X      RETURN
X   20 X = -LOG((B-P)/ALPHA)
X      IF (X**(ALPHA-1.0) .LT. RND()) GO TO 10
X      GS = X
X      RETURN
X      END
X
X      FUNCTION GCF (ALPHA)
X      REAL GCF,ALPHA,AA,APREV,C1,C2,C3,C4,C5,U1,U2,W,X
X      SAVE APREV
X      DATA APREV/0.0/
X      IF (ALPHA .EQ. APREV) GO TO 10
X      C1 = ALPHA-1.0
X      AA = 1.0/C1
X      C2 = AA*(ALPHA-1.0/(6.0*ALPHA))
X      C3 = 2.0*AA
X      C4 = C3+2.0
X      IF (ALPHA .GT. 2.5) C5 = 1.0/SQRT(ALPHA)
X   10 U1 = RND()
X      U2 = RND()
X      IF (ALPHA .LE. 2.5) GO TO 20
X      U1 = U2+C5*(1.0-1.86*U1)
X      IF (U1.LE.0.0 .OR. U1.GE.1.0) GO TO 10
X   20 W = C2*U2/U1
X      IF (C3*U1+W+1.0/W .LT. C4) GO TO 30
X      IF (C3*LOG(U1)-LOG(W)+W .GE. 1.0) GO TO 10
X   30 GCF = C1*W
X      APREV = ALPHA
X      RETURN
X      END
END_OF_FILE
if test 2881 -ne `wc -c <'normgam.f'`; then
    echo shar: \"'normgam.f'\" unpacked with wrong size!
fi
# end of 'normgam.f'
fi
if test -f 'proot.f' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'proot.f'\"
else
echo shar: Extracting \"'proot.f'\" \(1518 characters\)
sed "s/^X//" >'proot.f' <<'END_OF_FILE'
X      PROGRAM PROOT
X      DOUBLE PRECISION A(50),AA,FAC,M,MUL,MULT,PI,PWR
X      INTEGER I,J,R
X      PRINT *,'Multiplier,Modulus,No of factors of M-1'
X      READ *,MULT,M,R
X      AA = MULT
X      A(1) = MULT
X      DO 10 I = 2,50
X        AA = MUL (AA,AA,M)
X   10   A(I) = AA
X      DO 30 I = 1,R
X        PRINT *,'Factor',I
X        READ *,FAC
X        PWR = (M-1)/FAC
X        AA = 1.0D0
X        DO 20 J = 50,1,-1
X          PI = 2.0D0**(J-1)
X          IF (PWR .GE. PI) THEN
X            PWR = PWR-PI
X            AA = MUL (AA,A(J),M)
X          ENDIF
X   20     CONTINUE
X        IF (AA .EQ. 1.0D0) THEN
X          PRINT *,'FAILED'
X        ELSE
X          PRINT *,'OK for this factor'
X        ENDIF
X   30   CONTINUE
X      END
X
X      FUNCTION MUL(A,B,C)
XC     forms A*B MOD C for A,B < C <= 2^50
X      DOUBLE PRECISION A,B,C,MUL,A1,A2,B1,B2,D,M,MM
XC     set M so that 2*M*M <= maxint
X      M = 2.0D0**25
X      MM = M*M
X      D = (A*B)/C
X      CALL SPLIT (A,B,A1,A2)
X      CALL SPLIT (C,D,B1,B2)
X      MUL = A1-B1 + (A2-B2)*MM
X      RETURN
X      END
XC
X      SUBROUTINE SPLIT (A,B,C,D)
XC     forms A*B = C + D*MM for 0<= C,D < MM
X      DOUBLE PRECISION A,B,C,D,A1,A2,B1,B2,M,MM,AC,AD,C1,C2
XC     set M to the same value as in MUL
X      M = 2.0D0**25
X      MM = M*M
X      A2 = INT(A/M)
X      A1 = INT(A-A2*M)
X      B2 = INT(B/M)
X      B1 = INT(B-B2*M)
X      AC = A2*B1 + A1*B2
X      C2 = INT(AC/M)
X      C1 = AC - M*C2
X      AC = C1*M + A1*B1
X      AD = INT(AC/MM)
X      C = AC - AD*MM
X      D = AD + C2 + A2*B2
X      RETURN
X      END
X
END_OF_FILE
if test 1518 -ne `wc -c <'proot.f'`; then
    echo shar: \"'proot.f'\" unpacked with wrong size!
fi
# end of 'proot.f'
fi
echo shar: End of shell archive.
exit 0