SUBROUTINE RNFARR(AA,N)
C       FORTRAN 77 version of "ranf_array"
C       from Seminumerical Algorithms by D E Knuth, 3rd edition (1997)
C       ********* see the book for explanations and caveats! *********
C     Copyright (c) 2000 by the author;
C     copy it if you like, but don't change it!
      IMPLICIT DOUBLE PRECISION (A,R,X,Y)
      DIMENSION AA(*)
      PARAMETER (KK=100)
      PARAMETER (LL=37)
      COMMON /RSTATE/ RANX(KK)
      SAVE /RSTATE/
      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 RNFSTR(SEED)
      IMPLICIT DOUBLE PRECISION (A,R,U-Z)
      IMPLICIT INTEGER (T)
      PARAMETER (KK=100)
      PARAMETER (LL=37)
      PARAMETER (MM=2**30)
      PARAMETER (ULP=1D0/(2D0**52))
      PARAMETER (TT=70)
      PARAMETER (KKK=KK+KK-1)
      INTEGER SEED,S,SSEED
      DOUBLE PRECISION SS
      DIMENSION X(KKK), XL(KKK)
      COMMON /RSTATE/ RANX(KK)
      SAVE /RSTATE/
      IF (SEED .LT. 0) THEN
         SSEED=MM-1-MOD(-1-SEED,MM)
      ELSE
         SSEED=MOD(SEED,MM)
      END IF
      SS=2D0*ULP*(SSEED+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
      X(J)=0
 2    XL(J)=0
      S=SSEED
      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
C      a rudimentary test program:
       IMPLICIT DOUBLE PRECISION (A)
       DIMENSION A(2009)
       EXTERNAL RNFARR, RNFSTR
       CALL RNFSTR(310952)
       DO 1 I=1,2010
          CALL RNFARR(A,1009)
 1     CONTINUE 
       PRINT '(F22.20)',A(1)
C                   the number should be 0.27452626307394156768
       CALL RNFSTR(310952)
       DO 2 I=1,1010
          CALL RNFARR(A,2009)
 2     CONTINUE 
       PRINT '(F22.20)',A(1)
C                                 again, 0.27452626307394156768
       END