Pure random number generators

8

As @mecej4 observed, the code can be simplified if one is only generating uniform variates. I have implemented pure subroutines (not functions) to generate single values or sequences of random 32-bit integers or uniform variates, with an integer :: intent(in out) argument representing the state:

module ziggurat_pure_mod
! Marsaglia, G. & Tsang, W.W. (2000) `The ziggurat method for generating random variables', J. Statist. Software, v5(8).
implicit none
private
public             :: dp, random_int_32, random_uni
interface random_int_32
   module procedure random_int_32_scalar, random_int_32_vec
end interface random_int_32
interface random_uni
   module procedure random_uni_scalar,random_uni_vec
end interface random_uni
integer, parameter :: dp = selected_real_kind(12, 60)
contains
!
pure elemental subroutine random_int_32_scalar(jsr,iran)
! generate random 32-bit integers
integer, intent(in out) :: jsr  ! state of RNG
integer, intent(out)    :: iran ! random integer
integer                 :: jz
jz   = jsr
jsr  = ieor(jsr, ishft(jsr,  13))
jsr  = ieor(jsr, ishft(jsr, -17))
jsr  = ieor(jsr, ishft(jsr,   5))
iran = jz + jsr
end subroutine random_int_32_scalar
!
pure subroutine random_int_32_vec(jsr,iran)
! generate random 32-bit integers
integer, intent(in out) :: jsr     ! state of RNG
integer, intent(out)    :: iran(:) ! random integers
integer                 :: i
do i=1,size(iran)
   call random_int_32_scalar(jsr,iran(i))
end do
end subroutine random_int_32_vec
!
pure elemental subroutine random_uni_scalar(jsr,xran)
integer      , intent(in out) :: jsr  ! state of RNG
real(kind=dp), intent(out)    :: xran ! random uniform variate
integer                       :: iran
call random_int_32(jsr,iran)
xran = 0.2328306e-9_dp*iran + 0.5_dp
end subroutine random_uni_scalar
!
pure subroutine random_uni_vec(jsr,xran)
integer      , intent(in out) :: jsr     ! state of RNG
real(kind=dp), intent(out)    :: xran(:) ! random uniform variates
integer                       :: i
do i=1,size(xran)
   call random_uni_scalar(jsr,xran(i))
end do
end subroutine random_uni_vec
end module ziggurat_pure_mod

program xzig_pure
! 10/27/2021 12:54 AM driver for random_int_32 and random_uni
use ziggurat_pure_mod, only: dp, random_int_32, random_uni
implicit none
integer, parameter :: seed = 123, nran = 10000000
integer            :: iran(nran),state
real(kind=dp)      :: xran(nran)
state = seed
call random_int_32(state,iran)
call random_uni(state,xran)
write (*,"(2a12)") "min","max"
write (*,"(2i12)") minval(iran),maxval(iran)
write (*,"(/,4a10)") "min","max","mean","mean_sq"
write (*,"(4f10.6)") minval(xran),maxval(xran),sum(xran)/nran,sum(xran**2)/nran
end program xzig_pure

Compiling with gfortran ziggurat_pure.f90 xzig_pure.f90 or the analogous, gfortran, Intel Fortran, flang, and g95 all give output

         min         max
 -2147483306  2147482694

       min       max      mean   mean_sq
  0.000000  1.000000  0.500141  0.333496
1 Like