# Pure random number generators

**URL:** <https://fortran-lang.discourse.group/t/pure-random-number-generators/2170>\
**Category:** Uncategorized\
**Created:** [October 27, 2021, 2:02am UTC](https://fortran-lang.discourse.group/t/pure-random-number-generators/2170 "2021-10-27T02:02:28Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![Beliavsky](https://avatars.discourse-cdn.com/v4/letter/b/ba8739/32.png) [@Beliavsky](https://fortran-lang.discourse.group/u/Beliavsky)\
**Post date:** [October 27, 2021, 2:02am UTC](https://fortran-lang.discourse.group/t/pure-random-number-generators/2170/1 "2021-10-27T02:02:28Z")

</div>

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

```auto
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](https://github.com/Beliavsky/Ziggurat/blob/main/ziggurat.f90) gives output

```auto
  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.

---

<div class="post-metadata">

**Author:** ![MarDie](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/mardie/32/54_2.png) [@MarDie](https://fortran-lang.discourse.group/u/MarDie)\
**Post date:** [October 27, 2021, 4:19am UTC](https://fortran-lang.discourse.group/t/pure-random-number-generators/2170/2 "2021-10-27T04:19:33Z")

</div>

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

```auto
# 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](https://numpy.org/doc/stable/reference/random/index.html?highlight=random#module-numpy.random).  
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)`

---

<div class="post-metadata">

**Author:** ![CRquantum](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/crquantum/32/730_2.png) [@CRquantum](https://fortran-lang.discourse.group/u/CRquantum)\
**Post date:** [October 27, 2021, 7:12am UTC](https://fortran-lang.discourse.group/t/pure-random-number-generators/2170/3 "2021-10-27T07:12:17Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![mecej4](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/mecej4/32/855_2.png) [@mecej4](https://fortran-lang.discourse.group/u/mecej4)\
**Post date:** [October 27, 2021, 8:02am UTC](https://fortran-lang.discourse.group/t/pure-random-number-generators/2170/4 "2021-10-27T08:02:05Z")

</div>

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](https://en.wikipedia.org/wiki/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.

---

<div class="post-metadata">

**Author:** ![MarDie](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/mardie/32/54_2.png) [@MarDie](https://fortran-lang.discourse.group/u/MarDie)\
**Post date:** [October 27, 2021, 2:10pm UTC](https://fortran-lang.discourse.group/t/pure-random-number-generators/2170/5 "2021-10-27T14:10:14Z")

</div>

This will make it impossible to have reproducible results.

---

<div class="post-metadata">

**Author:** ![certik](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/certik/32/4_2.png) [@certik](https://fortran-lang.discourse.group/u/certik)\
**Post date:** [October 27, 2021, 2:16pm UTC](https://fortran-lang.discourse.group/t/pure-random-number-generators/2170/6 "2021-10-27T14:16:48Z")

</div>

See here: [Pseudorandom number generator · Issue #135 · fortran-lang/stdlib · GitHub](https://github.com/fortran-lang/stdlib/issues/135)

We should add that to stdlib.

---

<div class="post-metadata">

**Author:** ![Beliavsky](https://avatars.discourse-cdn.com/v4/letter/b/ba8739/32.png) [@Beliavsky](https://fortran-lang.discourse.group/u/Beliavsky)\
**Post date:** [October 27, 2021, 5:07pm UTC](https://fortran-lang.discourse.group/t/pure-random-number-generators/2170/8 "2021-10-27T17:07:36Z")

</div>

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:

```auto
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

```

```auto
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

```auto
         min max
 -2147483306 2147482694

       min max mean mean_sq
  0.000000 1.000000 0.500141 0.333496

```

---

<div class="post-metadata">

**Author:** ![everythingfunctional](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/everythingfunctional/32/176_2.png) [@everythingfunctional](https://fortran-lang.discourse.group/u/everythingfunctional)\
**Post date:** [October 29, 2021, 12:41am UTC](https://fortran-lang.discourse.group/t/pure-random-number-generators/2170/9 "2021-10-29T00:41:30Z")

</div>

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

```auto
rand_num, next_gen = gen%next_number()

```