Pure random number generators

1

Are there pure random number generators in Fortran? They could be used within pure procedures. You would pass the function a seed and it would return a vector of deviates that depends on the seed. The main program

program xzig
! 10/26/2021 01:38 PM demonstrate random number generator that is effectively pure
use ziggurat_mod, only: uni
implicit none
integer, parameter :: n = 5, nseeds = 2, seeds(nseeds) = [654321,123456]
integer :: iter, iseed
do iter=1,3
   do iseed=1,nseeds
      write (*,"(i8,1000f8.4)") seeds(iseed),uni(n,seeds(iseed))
   end do
end do
end program xzig

compiled with ziggurat.f90 gives output

  654321  0.3187  0.5605  0.5417  0.9146  0.5146
  123456  0.2089  0.2956  0.7174  0.3179  0.6881
  654321  0.3187  0.5605  0.5417  0.9146  0.5146
  123456  0.2089  0.2956  0.7174  0.3179  0.6881
  654321  0.3187  0.5605  0.5417  0.9146  0.5146
  123456  0.2089  0.2956  0.7174  0.3179  0.6881

so function uni(n,seed) is effectively pure in the sense of returning the same output based on its input. It cannot be declared pure, however.

1 Like
2

I see no reason why a self-implemented pseudo random number generator should not be pure. The build-in random number generator (random_number) is also a pseudo random number generator. However, it is based on the the value of a seed aka ‘state of the world’ and therefore has side-effects. random_seed sets the state of the world depending on the given seed and is therefore also non-pure.

numpy.random's solution is to have an object that has its own state

# Do this (new version)
from numpy.random import default_rng
rng = default_rng()
vals = rng.standard_normal(10)
more_vals = rng.standard_normal(10)

# instead of this (legacy version)
from numpy import random
vals = random.standard_normal(10)
more_vals = random.standard_normal(10)

Example taken from Random sampling (numpy.random) — NumPy v1.21 Manual.
In legacy mode, the state is global, in the new version the state is attached to the object called rng.

edit: it is important that the seed needs to be an argument to the call of the random number generator, otherwise it cannot be pure. Something like rnd(1:10) = random_numbers(seed=42,N=10)

3

I have the same question.
If there are pure RNGs, some of the loops can immediately be done by do concurrent, which can be great.

4

The following is an observation that is not related to the issue that @Beliavsky raised.

It seems to me that if ziggurat.f90 is used only for repeated invocations of uni(n,seed) the call to subroutine zigset in uni_vec_given_seed could be replaced by the simple assignment jsr = jsrseed, which is the first executable statement in zigset(). The rest of the code in zigset() is superfluous and wastefully expensive when only uniformly distributed random numbers are going to be generated. In other words, why expend the effort of setting up the Ziggurat when it is not going to be used?

[My comments are based on visual inspection of ziggurat.f90, and I could be mistaken.]

[Added 3 hours later] I confirmed my conjecture by inserting RETURN after the first executable statement in subroutine ZIGSET of ziggurat.f90 before building and running the program. The results were the same as what @Beliavsky reported above.

5

This will make it impossible to have reproducible results.

6

See here: Pseudorandom number generator · Issue #135 · fortran-lang/stdlib · GitHub

We should add that to stdlib.

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
9

In purely functional languages, random number generators are implemented where generating a random number returns a random number and a new generator (i.e. the state after updating the seed). Asking for a second number from the original generator will give you the same number as the previous time, because it doesn’t depend on any outside state. It’s an interesting technique, and is absolutely something that could be implemented in Fortran, it just wouldn’t seem very Fortranic to most people. If we had destructuring syntax a call might look something like

rand_num, next_gen = gen%next_number()
1 Like