MODULE Lin_Feedback_Shift_Reg ! L'Ecuyer's 1999 random number generator. ! Fortran version by Alan.Miller @ vic.cmis.csiro.au ! This version is for 64-bit integers and assumes that KIND=8 identifies them. ! It has only been tested using Lahey's LF95 compiler ! http://users.bigpond.net.au/amiller/ ! http://www.ozemail.com.au/~milleraj ! Latest revision - 12 January 2001 IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(14, 60) ! These are unsigned integers in the C version INTEGER (KIND=8), SAVE :: s1 = 153587801, s2 = -759022222, s3 = 1288503317, & s4 = -1718083407, s5 = -123456789 CONTAINS SUBROUTINE init_seeds(i1, i2, i3, i4, i5) IMPLICIT NONE INTEGER (KIND=8), INTENT(IN) :: i1, i2, i3, i4, i5 s1 = i1 s2 = i2 s3 = i3 s4 = i4 s5 = i5 IF (IAND(s1, -2) == 0) s1 = i1 - 8388607 IF (IAND(s2, -512) == 0) s2 = i2 - 8388607 IF (IAND(s3, -4096) == 0) s3 = i3 - 8388607 IF (IAND(s4, -131072) == 0) s4 = i4 - 8388607 IF (IAND(s5,-8388608) == 0) s5 = i5 - 8388607 RETURN END SUBROUTINE init_seeds FUNCTION lfsr258() 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^(258) or about 4.6 x 10^77. ! Actually - (2^63 - 1).(2^55 - 1).(2^52 - 1).(2^47 - 1).(2^41 - 1) IMPLICIT NONE REAL (dp) :: random_numb INTEGER (KIND=8) :: 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,1), s1), -53) s1 = IEOR( ISHFT( IAND(s1,-2), 10), b) b = ISHFT( IEOR( ISHFT(s2,24), s2), -50) s2 = IEOR( ISHFT( IAND(s2,-512), 5), b) b = ISHFT( IEOR( ISHFT(s3,3), s3), -23) s3 = IEOR( ISHFT( IAND(s3,-4096), 29), b) b = ISHFT( IEOR( ISHFT(s4,5), s4), -24) s4 = IEOR( ISHFT( IAND(s4,-131072), 23), b) b = ISHFT( IEOR( ISHFT(s5,3), s5), -33) s5 = IEOR( ISHFT( IAND(s5,-8388608), 8), b) ! The constant below is the reciprocal of (2^64 - 1) random_numb = IEOR( IEOR( IEOR( IEOR(s1,s2), s3), s4), s5) & * 5.4210108624275221E-20_dp + 0.5_dp RETURN END FUNCTION lfsr258 END MODULE Lin_Feedback_Shift_Reg PROGRAM t_lfsr258 USE Lin_Feedback_Shift_Reg IMPLICIT NONE INTEGER (KIND=8) :: i1, i2, i3, i4, i5 INTEGER (KIND=4) :: i, j, k, l, m, freq(0:15,0:15,0:15,0:15) INTEGER (KIND=8), ALLOCATABLE :: seed(:) REAL (dp) :: x, chi_sq, expctd, deg_freedom, stdev, & lower, upper WRITE(*, *)'Enter 5 integers as seeds: ' READ(*, *) i1, i2, i3, i4, i5 CALL init_seeds(i1, i2, i3, i4, i5) x = lfsr258() i = 16*x x = lfsr258() j = 16*x x = lfsr258() k = 16*x freq = 0 DO m= 1, 10000000 x = lfsr258() 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_lfsr258