ziggurat.f90

Go to the documentation of this file.

! 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
!
! Modified by D. M. scott in 2016 to assist vectorisation by the compiler.
MODULE ziggurat
   IMPLICIT NONE

   PRIVATE

! Define the number of integers in a vector, N_INT, and also the number of
! pseudo random integers to be generated by each call to shr3_vec, N_GEN.
! The size of the array used to hold the pseudo random numbers must be an
! integral multiple of N_INT. N_GEN must be no greater than the size of this array.
! Also, N_GEN must be a power of two.
!
! This code is suitable for processors implementing AVX (e.g. Sandy Bridge) for
! which the vector length is 4 and arrays used in loops should be 32 byte aligned.
! See https://en.wikipedia.org/wiki/Advanced_Vector_Extensions#CPUs_with_AVX
! for details of which processors implement which instruction sets.
   INTEGER, PARAMETER  ::  n_int = 4
   INTEGER, PARAMETER  ::  multiplier = 1
   INTEGER, PARAMETER  ::  array_size = multiplier*n_int
   INTEGER, PARAMETER  ::  array_size_m1 = array_size - 1
   INTEGER, PARAMETER  ::  n_gen = array_size                ! Must <= ARRAY_SIZE
   INTEGER, PARAMETER  ::  n_gen_m1 = n_gen - 1

   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, 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.

   INTEGER, DIMENSION(0:ARRAY_SIZE_M1), SAVE :: jz_vec, jsr_vec
   INTEGER, DIMENSION(0:ARRAY_SIZE_M1), SAVE :: rand_ints
   !DIR$ ATTRIBUTES ALIGN : 32 :: jz_vec
   !DIR$ ATTRIBUTES ALIGN : 32 :: jsr_vec
   INTEGER, SAVE :: index_shr3 = 0

   PUBLIC  :: zigset, shr3, uni, rnor, rexp


CONTAINS


SUBROUTINE zigset( jsrseed )

   INTEGER, INTENT(IN)  :: jsrseed

   INTEGER  :: i

   !  Set the seeds
   !DIR$ VECTOR ALIGNED
   DO i = 0, n_gen_m1
      jsr_vec(i) = 1 + i + n_gen*jsrseed
   END DO

   !  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 N_GEN pseudo random 32-bit integers.
FUNCTION shr3_vec( ) RESULT( ival_vec )
   INTEGER, DIMENSION(0:ARRAY_SIZE_M1)  ::  ival_vec
   INTEGER :: i

   IF( .NOT. initialized ) stop 'ERROR: shr3_vec has not been initialised (call zigset).'
   !DIR$ VECTOR ALIGNED
   DO i = 0, n_gen_m1
      jz_vec(i) = jsr_vec(i)
      jsr_vec(i) = ieor( jsr_vec(i), ishft( jsr_vec(i),  13 ) )
      jsr_vec(i) = ieor( jsr_vec(i), ishft( jsr_vec(i), -17 ) )
      jsr_vec(i) = ieor( jsr_vec(i), ishft( jsr_vec(i),   5 ) )
      ival_vec(i) = jz_vec(i) + jsr_vec(i)
   END DO
   RETURN
END FUNCTION shr3_vec



!  Generate pseudo random 32-bit integers.
FUNCTION shr3( ) RESULT( ival )
   INTEGER  ::  ival

   IF (index_shr3 == 0) THEN
      rand_ints = shr3_vec()
   END IF

   ival = rand_ints( index_shr3 )

   index_shr3 = index_shr3 + 1
   index_shr3 = iand( index_shr3, n_gen_m1 ) ! Calculate modulus.

   RETURN
END FUNCTION shr3



!  Generate uniformly distributed pseudo random numbers.
FUNCTION uni( ) RESULT( fn_val )
   REAL(DP)  ::  fn_val

   fn_val = half + 0.2328306e-9_dp * shr3( )
   RETURN
END FUNCTION uni



!  Generate pseudo random numbers with a normal distribution.
FUNCTION rnor( ) RESULT( fn_val )
   REAL(DP)             ::  fn_val

   REAL(DP), PARAMETER  ::  r = 3.442620_DP
   REAL(DP)             ::  x, y

   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 pseudo random numbers with an exponential distribution.
FUNCTION rexp( ) RESULT( fn_val )
   REAL(DP)  ::  fn_val

   REAL(DP)  ::  x

   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