! Marsaglia & Tsang generator for random normals & random exponentials. ! Translated from C by Alan Miller (amiller@bigpond.net.au) ! Marsaglia, G. & Tsang, W.W. (2000) `The ziggurat method for generating ! random variables', J. Statist. Software, v5(8). ! This is an electronic journal which can be downloaded from: ! http://www.jstatsoft.org/v05/i08 ! N.B. It is assumed that all integers are 32-bit. ! N.B. The value of M2 has been halved to compensate for the lack of ! unsigned integers in Fortran. ! Latest version - 1 January 2001 MODULE Ziggurat IMPLICIT NONE PRIVATE INTEGER, PARAMETER :: DP=SELECTED_REAL_KIND( 12, 60 ) REAL(DP), PARAMETER :: m1=2147483648.0_DP, m2=2147483648.0_DP, & half=0.5_DP REAL(DP) :: dn=3.442619855899_DP, tn=3.442619855899_DP, & vn=0.00991256303526217_DP, & q, de=7.697117470131487_DP, & te=7.697117470131487_DP, & ve=0.003949659822581572_DP INTEGER, SAVE :: iz, jz, jsr=123456789, kn(0:127), & ke(0:255), hz REAL(DP), SAVE :: wn(0:127), fn(0:127), we(0:255), fe(0:255) LOGICAL, SAVE :: initialized=.FALSE. PUBLIC :: zigset, shr3, uni, rnor, rexp CONTAINS SUBROUTINE zigset( jsrseed ) INTEGER, INTENT(IN) :: jsrseed INTEGER :: i ! Set the seed jsr = jsrseed ! Tables for RNOR q = vn*EXP(half*dn*dn) kn(0) = (dn/q)*m1 kn(1) = 0 wn(0) = q/m1 wn(127) = dn/m1 fn(0) = 1.0_DP fn(127) = EXP( -half*dn*dn ) DO i = 126, 1, -1 dn = SQRT( -2.0_DP * LOG( vn/dn + EXP( -half*dn*dn ) ) ) kn(i+1) = (dn/tn)*m1 tn = dn fn(i) = EXP(-half*dn*dn) wn(i) = dn/m1 END DO ! Tables for REXP q = ve*EXP( de ) ke(0) = (de/q)*m2 ke(1) = 0 we(0) = q/m2 we(255) = de/m2 fe(0) = 1.0_DP fe(255) = EXP( -de ) DO i = 254, 1, -1 de = -LOG( ve/de + EXP( -de ) ) ke(i+1) = m2 * (de/te) te = de fe(i) = EXP( -de ) we(i) = de/m2 END DO initialized = .TRUE. RETURN END SUBROUTINE zigset ! Generate random 32-bit integers FUNCTION shr3( ) RESULT( ival ) INTEGER :: ival jz = jsr jsr = IEOR( jsr, ISHFT( jsr, 13 ) ) jsr = IEOR( jsr, ISHFT( jsr, -17 ) ) jsr = IEOR( jsr, ISHFT( jsr, 5 ) ) ival = jz + jsr RETURN END FUNCTION shr3 ! Generate uniformly distributed random numbers FUNCTION uni( ) RESULT( fn_val ) REAL(DP) :: fn_val fn_val = half + 0.2328306e-9_DP * shr3( ) RETURN END FUNCTION uni ! Generate random normals FUNCTION rnor( ) RESULT( fn_val ) REAL(DP) :: fn_val REAL(DP), PARAMETER :: r = 3.442620_DP REAL(DP) :: x, y IF( .NOT. initialized ) CALL zigset( jsr ) hz = shr3( ) iz = IAND( hz, 127 ) IF( ABS( hz ) < kn(iz) ) THEN fn_val = hz * wn(iz) ELSE DO IF( iz == 0 ) THEN DO x = -0.2904764_DP* LOG( uni( ) ) y = -LOG( uni( ) ) IF( y+y >= x*x ) EXIT END DO fn_val = r+x IF( hz <= 0 ) fn_val = -fn_val RETURN END IF x = hz * wn(iz) IF( fn(iz) + uni( )*(fn(iz-1)-fn(iz)) < EXP(-half*x*x) ) THEN fn_val = x RETURN END IF hz = shr3( ) iz = IAND( hz, 127 ) IF( ABS( hz ) < kn(iz) ) THEN fn_val = hz * wn(iz) RETURN END IF END DO END IF RETURN END FUNCTION rnor ! Generate random exponentials FUNCTION rexp( ) RESULT( fn_val ) REAL(DP) :: fn_val REAL(DP) :: x IF( .NOT. initialized ) CALL Zigset( jsr ) jz = shr3( ) iz = IAND( jz, 255 ) IF( ABS( jz ) < ke(iz) ) THEN fn_val = ABS(jz) * we(iz) RETURN END IF DO IF( iz == 0 ) THEN fn_val = 7.69711 - LOG( uni( ) ) RETURN END IF x = ABS( jz ) * we(iz) IF( fe(iz) + uni( )*(fe(iz-1) - fe(iz)) < EXP( -x ) ) THEN fn_val = x RETURN END IF jz = shr3( ) iz = IAND( jz, 255 ) IF( ABS( jz ) < ke(iz) ) THEN fn_val = ABS( jz ) * we(iz) RETURN END IF END DO RETURN END FUNCTION rexp END MODULE ziggurat PROGRAM test_ziggurat USE Ziggurat IMPLICIT NONE INTEGER, PARAMETER :: DP = SELECTED_REAL_KIND(12, 60) REAL :: t1, t2 REAL(DP) :: x, chisq, expctd INTEGER :: i, ten_million=10000000, freq(0:999), j ! Set the random number seed WRITE(*, '(a)', ADVANCE='NO') ' Enter a non-zero integer: ' READ(*, *) i CALL zigset( i ) ! First time the generation of uniform, normal & exponential generators CALL CPU_TIME( t1 ) DO i = 1, ten_million x = uni( ) END DO CALL CPU_TIME( t2 ) WRITE(*, '(a, i10, a, f9.2, a)') ' Time to generate ', ten_million, & ' uniform random nos. = ', t2-t1, 'sec.' CALL CPU_TIME( t1 ) DO i = 1, ten_million x = rnor( ) END DO CALL CPU_TIME( t2 ) WRITE(*, '(a, i10, a, f9.2, a)') ' Time to generate ', ten_million, & ' normal random nos. = ', t2-t1, 'sec.' CALL CPU_TIME( t1 ) DO i = 1, ten_million x = rexp( ) END DO CALL CPU_TIME( t2 ) WRITE(*, '(a, i10, a, f9.2, a)') ' Time to generate ', ten_million, & ' exponential random nos. = ', t2-t1, 'sec.' WRITE(*, *) ! Now look at the distributions generated WRITE(*, *) 'Std. devn. of chi-squared with 999 d.of.f. = 44.70' WRITE(*, *) 'Values of chi-squared below should be within about 90. of 999.' freq = 0 DO i = 1, ten_million j = 1000 * uni( ) freq(j) = freq(j) +1 END DO chisq = 0.0_DP expctd = ten_million / 1000 DO j = 0, 999 chisq = chisq + (freq(j) - expctd)**2 / expctd END DO WRITE(*, '(a, f9.2, a)') ' Uniform distribution, chi-squared = ', chisq, & ' with 999 deg. of freedom' freq = 0 DO i = 1, ten_million j = 1000 * alnorm( rnor( ), .FALSE. ) freq(j) = freq(j) +1 END DO chisq = 0.0_DP expctd = ten_million / 1000 DO j = 0, 999 chisq = chisq + (freq(j) - expctd)**2 / expctd END DO WRITE(*, '(a, f9.2, a)') ' Normal distribution, chi-squared = ', chisq, & ' with 999 deg. of freedom' freq = 0 DO i = 1, ten_million j = 1000 * EXP( -rexp( ) ) IF( j > 999 .OR. j < 0 ) THEN WRITE(*, '(a, 2i10)') ' i, j = ', i, j ELSE freq(j) = freq(j) +1 END IF END DO chisq = 0.0_DP expctd = ten_million / 1000 DO j = 0, 999 chisq = chisq + (freq(j) - expctd)**2 / expctd END DO WRITE(*, '(a, f9.2, a)') ' Exponential distribution, chi-squared = ', & chisq, ' with 999 deg. of freedom' STOP CONTAINS ! Algorithm AS66 Applied Statistics (1973) vol.22, no.3 ! Evaluates the tail area of the standardised normal curve ! from x to infinity if upper is .true. or ! from minus infinity to x if upper is .false. ! ELF90-compatible version by Alan Miller ! Latest revision - 29 November 2001 FUNCTION alnorm( x, upper ) RESULT( fn_val ) REAL(DP), INTENT(IN) :: x LOGICAL, INTENT(IN) :: upper REAL(DP) :: fn_val ! Local variables REAL(DP), PARAMETER :: zero=0.0_DP, one=1.0_DP, half=0.5_DP, con=1.28_DP REAL(DP) :: z, y LOGICAL :: up ! Machine dependent constants REAL(DP), PARAMETER :: ltone = 7.0_DP, utzero = 18.66_DP REAL(DP), PARAMETER :: p = 0.398942280444_DP, q = 0.39990348504_DP, & r = 0.398942280385_DP, a1 = 5.75885480458_DP, & a2 = 2.62433121679_DP, a3 = 5.92885724438_DP, & b1 = -29.8213557807_DP, b2 = 48.6959930692_DP, & c1 = -3.8052E-8_DP, c2 = 3.98064794E-4_DP, & c3 = -0.151679116635_DP, c4 = 4.8385912808_DP, & c5 = 0.742380924027_DP, c6 = 3.99019417011_DP, & d1 = 1.00000615302_DP, d2 = 1.98615381364_DP, & d3 = 5.29330324926_DP, d4 = -15.1508972451_DP, & d5 = 30.789933034_DP up = upper z = x IF( z < zero ) THEN up = .NOT. up z = -z END IF IF( z <= ltone .OR. (up .AND. z <= utzero) ) THEN y = half*z*z IF( z > con ) THEN fn_val = r*EXP( -y )/(z+c1+d1/(z+c2+d2/(z+c3+d3/(z+c4+d4/(z+c5+d5/(z+c6)))))) ELSE fn_val = half - z*(p-q*y/(y+a1+b1/(y+a2+b2/(y+a3)))) END IF ELSE fn_val = zero END IF IF( .NOT. up ) fn_val = one - fn_val RETURN END FUNCTION alnorm END PROGRAM test_ziggurat