MODULE Ecuyer_random
! L'Ecuyer's 1996 random number generator.
! Fortran version by Alan.Miller @ vic.cmis.csiro.au
! N.B. This version is compatible with Lahey's ELF90
! http://www.ozemail.com.au/~milleraj
! Latest revision - 30 March 1999

IMPLICIT NONE
INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(14, 60)

! These are unsigned integers in the C version
INTEGER, SAVE :: s1 = 1234, s2 = -4567, s3 = 7890

CONTAINS

SUBROUTINE init_seeds(i1, i2, i3)
IMPLICIT NONE

INTEGER, INTENT(IN) :: i1, i2, i3

s1 = i1
s2 = i2
s3 = i3
IF (IAND(s1,-2) == 0) s1 = i1 - 1023
IF (IAND(s2,-8) == 0) s2 = i2 - 1023
IF (IAND(s3,-16) == 0) s3 = i3 - 1023

RETURN
END SUBROUTINE init_seeds



FUNCTION taus88() RESULT(random_numb)
! Generates a random number between 0 and 1.  Translated from C function in:
! Reference:
! L'Ecuyer, P. (1996) `Maximally equidistributed combined Tausworthe
! generators', Math. of Comput., 65, 203-213.

! The cycle length is claimed to be about 2^(88) or about 3E+26.
! Actually - (2^31 - 1).(2^29 - 1).(2^28 - 1).

IMPLICIT NONE
REAL (dp) :: random_numb

INTEGER   :: b

! N.B. ISHFT(i,j) is a bitwise (non-circular) shift operation;
!      to the left if j > 0, otherwise to the right.

b  = ISHFT( IEOR( ISHFT(s1,13), s1), -19)
s1 = IEOR( ISHFT( IAND(s1,-2), 12), b)
b  = ISHFT( IEOR( ISHFT(s2,2), s2), -25)
s2 = IEOR( ISHFT( IAND(s2,-8), 4), b)
b  = ISHFT( IEOR( ISHFT(s3,3), s3), -11)
s3 = IEOR( ISHFT( IAND(s3,-16), 17), b)
random_numb = IEOR( IEOR(s1,s2), s3) * 2.3283064365E-10_dp + 0.5_dp

RETURN
END FUNCTION taus88

END MODULE Ecuyer_random



PROGRAM t_taus88
USE Ecuyer_random
IMPLICIT NONE
INTEGER              :: i1, i2, i3, i, j, k, l, m, freq(0:15,0:15,0:15,0:15)
INTEGER, ALLOCATABLE :: seed(:)
REAL (dp)            :: x, chi_sq, expctd, deg_freedom, stdev, lower, upper

WRITE(*, *)'Enter 3 integers as seeds: '
READ(*, *) i1, i2, i3
CALL init_seeds(i1, i2, i3)

x = taus88()
i = 16*x
x = taus88()
j = 16*x
x = taus88()
k = 16*x
freq = 0
DO m= 1, 10000000
  x = taus88()
  l = 16*x
  freq(i,j,k,l) = freq(i,j,k,l) + 1
  i = j
  j = k
  k = l
END DO

expctd = REAL( SUM(freq) ) / (16. * 16. * 16. * 16.)

chi_sq = 0.0_dp
DO i = 0, 15
  DO j = 0, 15
    DO k = 0, 15
      chi_sq = chi_sq + SUM( (REAL(freq(i,j,k,:)) - expctd)**2 ) / expctd
    END DO
  END DO
END DO

deg_freedom = (16. * 16. * 16. * 16.) - 1.
WRITE(*, '(a, f10.1, a, f6.0, a)') ' Chi-squared = ', chi_sq, ' with ',  &
                                  deg_freedom, ' deg. of freedom'

! Now repeat the exercise with the compiler's random number generator.

CALL RANDOM_SEED(size=k)
ALLOCATE( seed(k) )
IF (i1 /= 0) THEN
  seed(1) = i1
ELSE
  seed(1) = 1234567
END IF
IF (k >= 2) seed(2) = i2
IF (k >= 3) seed(3) = i3
DO i = 4, k
  seed(i) = seed(i-3)
END DO
CALL RANDOM_SEED(put=seed)

CALL RANDOM_NUMBER(x)
i = 16*x
CALL RANDOM_NUMBER(x)
j = 16*x
CALL RANDOM_NUMBER(x)
k = 16*x
freq = 0
DO m= 1, 10000000
  CALL RANDOM_NUMBER(x)
  l = 16*x
  freq(i,j,k,l) = freq(i,j,k,l) + 1
  i = j
  j = k
  k = l
END DO

chi_sq = 0.0_dp
DO i = 0, 15
  DO j = 0, 15
    DO k = 0, 15
      chi_sq = chi_sq + SUM( (REAL(freq(i,j,k,:)) - expctd)**2 ) / expctd
    END DO
  END DO
END DO

deg_freedom = (16. * 16. * 16. * 16.) - 1.
WRITE(*, '(a, f10.1, a, f6.0, a)') ' Chi-squared = ', chi_sq, ' with ',  &
                                  deg_freedom, ' deg. of freedom'

! Calculate rough limits for chi-squared based upon a normal approximation.
stdev = SQRT(2. * deg_freedom)
upper = deg_freedom + 2.*stdev
lower = deg_freedom - 2.*stdev
WRITE(*, '(a, f8.1, a, f8.1)') ' Approx. 95% limits for chi-squared: ',  &
            lower, ' - ', upper

STOP
END PROGRAM t_taus88