MODULE Knuth_random

! Code converted using TO_F90 by Alan Miller
! Date: 2000-09-10  Time: 16:37:48
! Latest revision - 16 January 2003

! FORTRAN 77 version of "ran_array"
! from Seminumerical Algorithms by D E Knuth, 3rd edition (1997)
!       including the MODIFICATIONS made in the 9th printing (2002)
! ********* see the book for explanations and caveats! *********
! Author: Steve Kifowit
! http://ourworld.compuserve.com/homepages/steve_kifowit
! with modifications by Alan Miller to rnarry and rnstrt based upon
! Knuth's code.

! For Donald Knuth's Fortran 77 versions, go to:
! http://www-cs-faculty.stanford.edu/~knuth/programs
! Look for frng.f and frngdb.f

IMPLICIT NONE
INTEGER, PARAMETER  :: kk=100, ll=37, mm=2**30, tt=70, kkk=kk+kk-1
INTEGER, SAVE       :: ranx(kk)

CONTAINS


SUBROUTINE rand(u, n)

! Generate an array of n real values between 0 and 1.

REAL, INTENT(OUT)    :: u(:)
INTEGER, INTENT(IN)  :: n

! Local array
INTEGER  :: aa(n)

CALL rnarry(aa, n)
u(1:n) = SCALE( REAL(aa), -30)

RETURN
END SUBROUTINE rand



SUBROUTINE rnarry(aa, n)

! Generate an array of n integers between 0 and 2^30-1.
 
INTEGER, INTENT(OUT)  :: aa(:)
INTEGER, INTENT(IN)   :: n

! Local variables
INTEGER  :: j

aa(1:kk) = ranx(1:kk)
DO  j = kk + 1, n
  aa(j) = aa(j-kk) - aa(j-ll)
  IF (aa(j) < 0) aa(j) = aa(j) + mm
END DO
DO  j=1,ll
  ranx(j) = aa(n+j-kk) - aa(n+j-ll)
  IF (ranx(j) < 0) ranx(j) = ranx(j) + mm
END DO
DO  j=ll+1,kk
  ranx(j) = aa(n+j-kk) - ranx(j-ll)
  IF (ranx(j) < 0) ranx(j) = ranx(j) + mm
END DO

RETURN
END SUBROUTINE rnarry



SUBROUTINE rnstrt(seed)

! Initialize integer array ranx using the input seed.

INTEGER, INTENT(IN)  :: seed

! Local variables
INTEGER  :: x(kkk), j, ss, sseed, t

IF (seed < 0) THEN
  sseed = mm - 1 - MOD(-1-seed, mm)
ELSE
  sseed = MOD(seed, mm)
END IF
ss = sseed - MOD(sseed,2) + 2
DO  j=1, kk
  x(j) = ss
  ss = ISHFT(ss, 1)
  IF (ss >= mm) ss = ss - mm + 2
END DO
x(kk+1:kkk) = 0
x(2) = x(2)+1
ss = sseed
t = tt - 1
10 DO  j=kk, 2, -1
  x(j+j-1) = x(j)
END DO
DO  j = kkk, kk + 1, -1
  x(j-(kk-ll)) = x(j-(kk-ll)) - x(j)
  IF (x(j-(kk-ll)) < 0) x(j-(kk-ll)) = x(j-(kk-ll)) + mm
  x(j-kk) = x(j-kk) - x(j)
  IF (x(j-kk) < 0) x(j-kk) = x(j-kk) + mm
END DO
IF (MOD(ss,2) == 1) THEN
  DO  j=kk, 1, -1
    x(j+1) = x(j)
  END DO
  x(1) = x(kk+1)
  x(ll+1) = x(ll+1) - x(kk+1)
  IF (x(ll+1) < 0) x(ll+1) = x(ll+1) + mm
END IF
IF (ss /= 0) THEN
  ss = ISHFT(ss, -1)
ELSE
  t = t - 1
END IF
IF (t > 0) GO TO 10

DO  j=1, ll
  ranx(j+kk-ll) = x(j)
END DO
DO  j=ll+1,kk
  ranx(j-ll) = x(j)
END DO

DO  j = 1, 10
  CALL rnarry(x,kkk)
END DO

RETURN
END SUBROUTINE rnstrt

END MODULE Knuth_random



PROGRAM main
!      a rudimentary test program:
USE Knuth_random
IMPLICIT NONE

INTEGER  :: a(2009), i
REAL     :: u(1000), aver, dev, sumsq, stdev

CALL rnstrt(310952)
DO  i=1,2010
  CALL rnarry(a, 1009)
END DO
WRITE(*, '(I15)') a(1)
!                   the number should be 211101277
CALL rnstrt(310952)
DO  i=1,1010
  CALL rnarry(a, 2009)
END DO
WRITE(*, '(I15)') a(1)
!                                 again, 211101277

! Generate real variables.

CALL rand(u, 1000)
aver = 0.0
sumsq = 0.0
DO i = 1, 1000
  dev = u(i) - aver
  aver = aver + dev/i
  sumsq = sumsq + dev*(u(i) - aver)
END DO

! N.B. Dividing the sum of squares by (n-1) gives an unbiassed estimate
! of the variance; dividing by (n-1.5) gives an almost unbiassed estimate
! of the standard deviation.
stdev = SQRT(sumsq/998.5)

WRITE(*, '(a, f9.4, a, f9.4)') ' Mean = ', aver, '  Std. devn. = ', stdev
WRITE(*, '(a, f6.1, "  ", f9.4)')  &
         ' Values above should be close to ', 0.5, SQRT(1./12.)
STOP
END PROGRAM main