MODULE Lin_Feedback_Shift_Reg
! L'Ecuyer's 1999 random number generator.
! Fortran version by Alan.Miller @ vic.cmis.csiro.au
! N.B. This version is compatible with Lahey's ELF90
! This version requires that the default integer type is of 32-bits
! http://www.ozemail.com.au/~milleraj
! http://users.bigpond.net.au/amiller/
! Latest revision - 12 January 2001

IMPLICIT NONE
INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(14, 60)
! These are unsigned integers in the C version
INTEGER, SAVE      :: s1 = 153587801, s2 = -759022222, s3 = 1288503317, &
                      s4 = -1718083407

CONTAINS

SUBROUTINE init_seeds(i1, i2, i3, i4)
IMPLICIT NONE

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

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

RETURN
END SUBROUTINE init_seeds



FUNCTION lfsr113() RESULT(random_numb)
! Generates a random number between 0 and 1.  Translated from C function in:
! Reference:
! L'Ecuyer, P. (1999) `Tables of maximally equidistributed combined LFSR
! generators', Math. of Comput., 68, 261-269.

! The cycle length is claimed to be about 2^(113) or about 10^(34).
! Actually - (2^31 - 1).(2^29 - 1).(2^28 - 1).(2^25 - 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,6), s1), -13)
s1 = IEOR( ISHFT( IAND(s1,-2), 18), b)
b  = ISHFT( IEOR( ISHFT(s2,2), s2), -27)
s2 = IEOR( ISHFT( IAND(s2,-8), 2), b)
b  = ISHFT( IEOR( ISHFT(s3,13), s3), -21)
s3 = IEOR( ISHFT( IAND(s3,-16), 7), b)
b  = ISHFT( IEOR( ISHFT(s4,3), s4), -12)
s4 = IEOR( ISHFT( IAND(s4,-128), 13), b)

! The constant below is the reciprocal of (2^32 - 1)
random_numb = IEOR( IEOR( IEOR(s1,s2), s3), s4) * 2.3283064365E-10_dp + 0.5_dp

RETURN
END FUNCTION lfsr113

END MODULE Lin_Feedback_Shift_Reg



PROGRAM t_lfsr113
USE Lin_Feedback_Shift_Reg
IMPLICIT NONE
INTEGER              :: i1, i2, i3, i4, 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 4 integers as seeds: '
READ(*, *) i1, i2, i3, i4
CALL init_seeds(i1, i2, i3, i4)

x = lfsr113()
i = 16*x
x = lfsr113()
j = 16*x
x = lfsr113()
k = 16*x
freq = 0
DO m= 1, 10000000
  x = lfsr113()
  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

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

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_lfsr113