SUBROUTINE RNARRY(AA,N)
C       FORTRAN 77 version of "ran_array"
C       I'm not sure that this is the debugged version, but it's pretty close
      IMPLICIT DOUBLE PRECISION (A,R,X,Y)
      DIMENSION AA(*)
      PARAMETER (KK=100)
      PARAMETER (LL=37)
      COMMON /RSTATE/ RANX(KK)
      Z=3.14
      ZZ=IDINT(Z)
      DO 1 J=1,KK
 1    AA(J)=RANX(J)
      DO 2 J=KK+1,N
         Y=AA(J-KK)+AA(J-LL)
         AA(J)=Y-IDINT(Y)
 2    CONTINUE
      DO 3 J=1,LL
         Y=AA(N+J-KK)+AA(N+J-LL)
         RANX(J)=Y-IDINT(Y)
 3    CONTINUE
      DO 4 J=LL+1,KK
         Y=AA(N+J-KK)+RANX(J-LL)
         RANX(J)=Y-IDINT(Y)
 4    CONTINUE
      END

      SUBROUTINE RNSTRT(SEED)
      IMPLICIT DOUBLE PRECISION (A,R,U-Z)
      IMPLICIT INTEGER (T)
      PARAMETER (KK=100)
      PARAMETER (LL=37)
      PARAMETER (ULP=1d0/(2d0**52))
      PARAMETER (TT=70)
      PARAMETER (KKK=KK+KK-1)
      INTEGER SEED,S
      DOUBLE PRECISION SS
      DIMENSION X(KKK), XL(KKK)
      COMMON /RSTATE/ RANX(KK)
      SS=2d0*ULP*(SEED+2)
      DO 1 J=1,KK
         X(J)=SS
         XL(J)=0
         SS=SS+SS
         IF (SS .GE. 1d0) SS=SS-1d0+2*ULP
 1    CONTINUE
      X(2)=X(2)+ULP
      XL(2)=ULP
      DO 2 J=KK+1,KKK
 2    X(J)=0
      S=SEED
      T=TT-1
 10   DO 12 J=KK,2,-1
         XL(J+J-1)=XL(J)
 12      X(J+J-1)=X(J)
      DO 13 J=KKK,KK-LL+1,-2
         XL(KKK-J+2)=0
 13      X(KKK-J+2)=X(J)-XL(J)
      DO 14 J=KKK,KK+1,-1
         IF (XL(J) .NE. 0) THEN
            XL(J-(KK-LL))=ULP-XL(J-(KK-LL))
            Y=X(J-(KK-LL))+X(J)
            X(J-(KK-LL))=Y-IDINT(Y)
            XL(J-KK)=ULP-XL(J-KK)
            Y=X(J-KK)+X(J)
            X(J-KK)=Y-IDINT(Y)
         END IF
 14   CONTINUE
C     this part is redundant but useful while I'm debugging
      DO 15 J=1,KKK
      IF (DMOD(X(J),2*ULP) .NE. XL(J)) THEN
         PRINT '(I15)', J
      END IF
 15   CONTINUE
      IF (MOD(S,2) .EQ. 1) THEN
         DO 16 J=KK,1,-1
            XL(J+1)=XL(J)
 16         X(J+1)=X(J)
         X(1)=X(KK+1)
         XL(1)=XL(KK+1)
         IF (XL(KK+1) .NE. 0) THEN
            XL(LL+1)=ULP-XL(LL+1)
            Y=X(LL+1)+X(KK+1)
            X(LL+1)=Y-IDINT(Y)
         END IF
      END IF
      IF (S .NE. 0) THEN
         S=S/2
      ELSE
         T=T-1
      END IF
      IF (T .GT. 0) GO TO 10
      DO 20 J=1,LL
 20   RANX(J+KK-LL)=X(J)
      DO 21 J=LL+1,KK
 21   RANX(J-LL)=X(J)
      END



       PROGRAM MAIN
       IMPLICIT DOUBLE PRECISION (A)
       DIMENSION A(1000)
       DO 1 I=0,6
          CALL RNSTRT(I)
          CALL RNARRY(A,1000)
          PRINT '(F20.18)',A(1)
          CALL RNARRY(A,1000)
          PRINT '(F20.18)',A(1)
 1     CONTINUE 
       END