Pure random number generators

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

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)

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.

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.

This will make it impossible to have reproducible results.

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

We should add that to stdlib.

I doubt you can do it in pure Fortran (pun intended). The constraints on pure
will cause you problems (e.g., C1589 and C1594 may come into play). Now,
you can likely fake it with ISO bind C and real random number generator.

% cat bar.c
#include <stdint.h>
#include <stdio.h>
float bar(void)
  FILE *fp;
  union {
    uint32_t i;
    float x;
  } u;
  fp = fopen("/dev/random","r");
  fread(&u.i, 1, sizeof(int), fp);
  u.i = ((u.i << 9) >> 9) | 0x3f800000;  /* Little endian */
  return (u.x - 1);

module pngm
   implicit none
   public png
      pure function png() result(r)
            pure function bar() bind(c,name='bar')
               real bar
            end function bar
         end interface
         real r
         r = bar()
      end function png
end module pngm

program foo
   use pngm
   implicit none
   integer i
   do i = 1, 10
      print *, png()
   end do
end program foo

This, of course, will be slow due to opening and closing /dev/random
on each call. You might be able to have a C png_init() routine to open
/dev/random once and leave the FILE *fp open. This also assumes
/dev/random is backed by a real random number generator.

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
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)
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

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