Reproducible use of random_number with a single integer seed

1

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(), 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:

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: 
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..

2 Likes
2

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

How about setting all seeds to the same value like below

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
rand1200Ă—1200 11.2 KB

3

The seeds should not have low entropy, according to NAG. 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).

4

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:

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 :slight_smile: