# Reproducible use of random\_number with a single integer seed

**URL:** <https://fortran-lang.discourse.group/t/reproducible-use-of-random-number-with-a-single-integer-seed/3593>\
**Category:** Uncategorized\
**Created:** [May 29, 2022, 11:38pm UTC](https://fortran-lang.discourse.group/t/reproducible-use-of-random-number-with-a-single-integer-seed/3593 "2022-05-29T23:38:35Z")\
**Posts on this page:** 4\
**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:** [May 29, 2022, 11:38pm UTC](https://fortran-lang.discourse.group/t/reproducible-use-of-random-number-with-a-single-integer-seed/3593/1 "2022-05-29T23:38:35Z")

</div>

It is inconvenient that Fortran compilers require different, and possibly large, numbers of random seeds for use in

`call random_seed(put=seeds(:))`

Fortran 2018 introduced [random\_init()](https://gcc.gnu.org/onlinedocs/gfortran/RANDOM_005fINIT.html), but I think this went too far in the other direction, not allowing the user to specify a single integer seed. I suggest that a convenient but flexible way to use `random_number` to generate reproduce random numbers is to

(1) get an integer seed from the user  
(2) use that seed to generate the random integers needed for a call to `random_seed`  
(3) call `random_seed` with those integers

This is done in the program below, which gives reproducible results for a given compiler:

```auto
module random_init_mod
implicit none
private
public :: random_init_zig
contains
pure elemental subroutine random_int_32_scalar(jsr,iran)
! Marsaglia, G. & Tsang, W.W. (2000) `The ziggurat method for generating random variables', J. Statist. Software, v5(8).
! 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
!
subroutine random_init_zig(iseed)
! given an integer seed, generate the remaining seeds needed for
! call random_seed(put=seeds)
integer, intent(in) :: iseed
integer :: i, jseed, nseeds
integer, allocatable :: seeds(:)
logical, parameter :: print_seeds = .true. ! .false. in production
jseed = iseed
call random_seed(size=nseeds)
allocate (seeds(nseeds))
do i=1,nseeds
   call random_int_32_scalar(jseed,seeds(i))
end do
if (print_seeds) print "('seeds =',*(1x,i0))", seeds
call random_seed(put=seeds)
end subroutine random_init_zig
end module random_init_mod

program xrandom_init_zig
use random_init_mod, only: random_init_zig
implicit none
integer, parameter :: dp = kind(1.0d0), nran = 10
real(kind=dp) :: xran(nran)
integer :: iseed
do iseed=1,10
   call random_init_zig(iseed)
   call random_number(xran)
   print "('uniform deviates =',*(f6.3))", xran
end do
end program xrandom_init_zig

```

> **reproducible gfortran output on Windows and WSL2:**
>
> ```auto
> seeds = 270370 67905058 -1579897146 -1339932140 -1588678368 -1150782559 1377930986 1068191692
> uniform deviates = 0.615 0.439 0.702 0.149 0.711 0.285 0.860 0.386 0.472 0.783
> seeds = 540740 134794308 832136324 -1926131320 1195958560 343972669 419821180 1002759439
> uniform deviates = 0.060 0.210 0.750 0.667 0.851 0.127 0.199 0.409 0.939 0.599
> seeds = 811110 202697318 -1067648438 639733904 -552181740 -1979160386 -364455626 402533703
> uniform deviates = 0.872 0.649 0.299 0.443 0.127 0.549 0.648 0.378 0.274 0.845
> seeds = 1081480 269588616 1664272648 442704656 -1903050176 688961179 614142560 672690764
> uniform deviates = 0.238 0.911 0.423 0.729 0.526 0.847 0.042 0.218 0.419 0.130
> seeds = 1351850 337493674 -486131770 -1537080060 519975720 -1753709156 892851866 1680700356
> uniform deviates = 0.915 0.799 0.367 0.501 0.902 0.759 0.017 0.644 0.357 0.264
> seeds = 1622220 404378828 -1834287732 525546640 -914734280 1982893426 2084519396 1635897477
> uniform deviates = 0.797 0.514 0.772 0.018 0.824 0.929 0.422 0.133 0.265 0.931
> seeds = 1892590 472281838 57627466 -1637739640 1629864516 658030227 -1938003378 1053349649
> uniform deviates = 0.621 0.894 0.096 0.138 0.845 0.288 0.167 0.383 0.757 0.879
> seeds = 2162960 539177232 -966421967 893788835 -1707914490 1702002840 1393213960 -866571821
> uniform deviates = 0.310 0.699 0.107 0.696 0.534 0.777 0.387 0.930 0.953 0.909
> seeds = 2433330 607082290 1606025461 -664329183 576779922 -67613895 1403627202 -1478215093
> uniform deviates = 0.052 0.045 0.871 0.692 0.889 0.790 0.768 0.046 0.367 0.331
> seeds = 2703700 673971540 -1208060235 2104898599 -698672110 1406281669 2077150676 -1847819002
> uniform deviates = 0.229 0.400 0.781 0.325 0.695 0.913 0.109 0.533 0.477 0.944
> 
> ```

Should the `random_init()` intrinsic be expanded to allow a single integer seed to be given? Compiler writers would generate the seeds needed for a call to `random_seed()` in a deterministic way when the `REPEATABLE` argument is `.true.`.

---

<div class="post-metadata">

**Author:** ![Abhilash](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/abhilash/32/1635_2.png) [@Abhilash](https://fortran-lang.discourse.group/u/Abhilash)\
**Post date:** [May 30, 2022, 12:32am UTC](https://fortran-lang.discourse.group/t/reproducible-use-of-random-number-with-a-single-integer-seed/3593/2 "2022-05-30T00:32:09Z")

</div>

Do the seed values have to have some property to be “good”?

How about setting all seeds to the same value like below

```fortran
program main
    integer, allocatable:: seed(:)
    integer nseed, stat
    integer, parameter :: dp = kind(1.0d0), nran = 10000
    real(kind=dp) :: xran(nran)

    call random_seed(size = nseed)
    allocate(seed(nseed))

    seed=333331
    call random_seed(put=seed)
    deallocate(seed)
    call random_number(xran)

    open(unit=1,file='rand.bin', iostat=stat,status='unknown')
    if(stat==0) close(1,status="delete")

    open(unit=1,file='rand.bin', access='stream')
    write(1)xran
    close(1)
end program main

```

It produced the following distribution:

 ![rand](https://global.discourse-cdn.com/free1/uploads/fortran_lang/original/2X/3/3537ae36de3d41da818e9a70393003986f575467.png)

---

<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:** [May 30, 2022, 12:40am UTC](https://fortran-lang.discourse.group/t/reproducible-use-of-random-number-with-a-single-integer-seed/3593/3 "2022-05-30T00:40:22Z")

</div>

> [@Abhilash](#):
>
> Do the seed values have to have some property to be “good”?

The seeds should not have low entropy, [according to NAG](https://www.nag.com/nagware/np/r62_doc/manual/compiler_2_12.html). I think I read something similar about the seeds for the Intel Fortran RNG. The aim of my `random_init_zig` subroutine is to let the user specify high-entropy seeds by supplying a single integer.

> The random number generator supplied as the intrinsic subroutine RANDOM\_NUMBER is the “Mersenne Twister”.
> 
> Note that this generator has a large state (630 32-bit integers) and an extremely long period (approx 106000), and therefore it is **strongly recommended** that the RANDOM\_SEED routine only be used with a PUT argument that is the value returned by a previous call with GET; i.e., only to repeat a previous sequence. This is because if a user-specified seed has low entropy (likely since there are 630 values to be supplied), it is highly likely to set the generator to an apparently-low-entropy part of the sequence.
> 
> If you do want to provide your own seed (and thus entropy), you should store your values in the initial elements of the seed array and set all the remaining elements to zero — trailing zero elements will be ignored and not used to initialise the generator. Note that the seed is a random bitstream, and is therefore expected to have approximately half of its bits nonzero (thus providing many small integer values will likely result in a low-entropy part of the Mersenne Twister sequence being reached).

---

<div class="post-metadata">

**Author:** ![Abhilash](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/abhilash/32/1635_2.png) [@Abhilash](https://fortran-lang.discourse.group/u/Abhilash)\
**Post date:** [May 30, 2022, 1:34am UTC](https://fortran-lang.discourse.group/t/reproducible-use-of-random-number-with-a-single-integer-seed/3593/4 "2022-05-30T01:34:04Z")

</div>

Thank you! I think your subroutine will be helpful indeed. I have also found specifying the seed to the RNG inconvenient. Having a robust procedure is great.

As I am writing this, I have another bright idea to boot strap the `random_number` using itself as shown below to get different integer seeds:

```fortran
program main
    integer, parameter :: dp = kind(1.0d0)
    integer nseed
    integer, allocatable:: seed(:)
    real(kind=dp), allocatable:: tmp(:)

    call random_seed(size = nseed)
    allocate(seed(nseed))
    allocate(tmp(nseed))

    seed=333331
    call random_seed(put=seed)
    call random_number(tmp) !maybe put this statement in a short do loop
    seed=floor(huge(seed)*tmp)
    call random_seed(put=seed)

    do i=1,5
        call random_number(tmp)
        call random_seed(get=seed)
        print "(100(I12,x))", seed
    enddo

end program main

```

Hopefully I will get around to doing a small quantitative analysis of effect of seed choice some day 🙂