Random Numbers

Fortran 90 introduced a pseudo-random number generator (PRNG) to the language standard. The routine random_number() returns a single number or an array of numbers from the uniform distribution over the range 0 ≤ x < 1.

GNU Fortran also provides the older extensions srand(), rand() and irand() for backwards compatibility with g77. These routines access an independent PRNG. Additionally, several implementations of the Mersenne Twister (MT19937) have been written in Fortran that return random 32-bit or 64-bit values.

Fortran 90 also added the random_seed() routine that initialises the PRNG. Without seeding, random_number() will always return the same sequence of numbers.

call random_seed([size][, put][, get])
Argument Type Description
size integer The returned minimum size of the array used with put and get (optional).
put integer, dimension(*) The array with seed values to use (optional).
get integer, dimension(*) The used array with seed values (optional).

The GNU Fortran implementation uses the Xoshiro256** pseudo-random number generator, with a period of 2256 – 1. Threads defined by OpenMP directives have their own random number state. The routine random_number() returns either a random real number or an array of random numbers.

call random_number(harvest)
Argument Type Description
harvest real[, dimension(*)] Returned scalar or array with random number(s).

The example initialises the PRNG and fills and array with random numbers:

! random.f90
program main
    implicit none
    real :: r(5)

    call random_seed()
    call random_number(r)

    print '(5(f8.6, 1x))', r
end program main

Compile and run the example with:

$ gfortran13 -o random random.f90
$ ./random
0.325213 0.037541 0.241239 0.077683 0.846473

Fortran 2018

In Fortran 2018, the routine random_init() has been added to initialise the state of the PRNG with respect to Co-Arrays.

call random_init(repeatable, image_distinct)
Argument Type Description
repeatable logical If .true., the seed is set to a processor-dependent value that is the same each time random_init() is called from the same instance of program execution. If .false., the seed is set to a processor-dependent value.
image_distinct logical If .true., the seed is set to a processor-dependent value that is distinct from the seed set by a call to random_init() in another process. If .false., the seed depends on which image called random_init().

The routine random_init() allows us to generate the very same random numbers on each initialisation:

! random2018.f90
program main
    implicit none
    real :: r(5)

    call random_init(repeatable=.true., image_distinct=.true.)
    call random_number(r)
    print '(5(f8.6))', r

    call random_init(repeatable=.true., image_distinct=.true.)
    call random_number(r)
    print '(5(f8.6))', r
end program main

The array is filled with the same random numbers despite re-initialisation:

$ gfortran13 -o random2018 random2018.f90
$ ./random2018
0.083351 0.977660 0.863393 0.669568 0.659231
0.083351 0.977660 0.863393 0.669568 0.659231

Seeding

For higher entropy, we can feed a custom seed array to random_seed(), with data obtained from /dev/urandom:

! urandom.f90
program main
    implicit none
    integer :: i, n
    real    :: r(64)

    ! Initialise random number generator.
    call random_seed(size=n)               ! Get size of seed array.
    call random_seed(put=urandom_seed(n))  ! Put seed array into PRNG.

    ! Generate and output values.
    call random_number(r)
    print '(8(f8.6, 1x))', [ (r(i), i = 1, size(r)) ]
contains
    function urandom_seed(n, stat) result(seed)
        !! Returns a seed array filled with random values from `/dev/urandom`.
        integer, intent(in)            :: n
        integer, intent(out), optional :: stat
        integer                        :: seed(n)

        integer :: fu, rc

        open (access='stream', action='read', file='/dev/urandom', &
              form='unformatted', iostat=rc, newunit=fu)
        if (present(stat)) stat = rc
        if (rc == 0) read (fu) seed
        close (fu)
    end function urandom_seed
end program main

The example program outputs 64 random values:

$ gfortran13 -o urandom urandom.f90
$ ./urandom
0.474257 0.743894 0.944154 0.483623 0.600250 0.315834 0.778809 0.910940
0.538556 0.133123 0.613712 0.619409 0.426239 0.300121 0.989344 0.204518
0.744600 0.119839 0.033373 0.665467 0.795407 0.529277 0.868635 0.488096
0.940779 0.474789 0.680964 0.208645 0.681481 0.886648 0.027110 0.860485
0.826086 0.242641 0.759755 0.102684 0.808770 0.957726 0.359169 0.488995
0.626006 0.464330 0.721632 0.069854 0.248335 0.448519 0.376246 0.805929
0.559363 0.198283 0.137510 0.299755 0.283269 0.067363 0.222422 0.551090
0.661703 0.671904 0.254210 0.674276 0.599526 0.958398 0.852826 0.975271

OpenSSL

The RAND_bytes() function of OpenSSL returns cryptographically strong pseudo-random bytes. We first need to implement an appropriate ISO C binding interface to the C function. A character array can then be passed to the function interface to be filled with random values:

! crypto.f90
module openssl
    use, intrinsic :: iso_c_binding, only: c_char, c_int
    implicit none

    interface
        ! int RAND_bytes(unsigned char *buf, int num)
        function c_rand_bytes(buf, num) bind(c, name='RAND_bytes')
            import :: c_char, c_int
            implicit none
            character(kind=c_char), intent(inout)     :: buf(*)
            integer(kind=c_int),    intent(in), value :: num
            integer(kind=c_int)                       :: c_rand_bytes
        end function c_rand_bytes
    end interface
end module openssl

program main
    use :: openssl
    implicit none
    character :: bytes(128)
    integer   :: i, rc

    rc = c_rand_bytes(bytes, size(bytes))
    if (rc /= 1) stop 'Error: RAND_bytes() failed'

    print '(16(i4))', [ (ichar(bytes(i)), i = 1, size(bytes)) ]
end program main

The example has to be linked against OpenSSL’s crypto library, using flag -lcrypto:

$ gfortran13 -o crypto crypto.f90 -lcrypto

The program then just prints 128 random bytes [0 … 255] to console:

$ ./crypto
  98  85  79  19 129 221  19 228 182 145 102 113 161  23 139 216
 215  35 144 184  83 234 167 235  85  67 242 145 128  23 227 247
 113  13 115 166 170  35  44  59  53 242 231 191 139 157  84 138
 140 150  48 182 100  28 241  33 234 111  45 106 246 112 163 211
  51 108 237 223 178  65 206  23 126 248  41 152   8 221  54  63
   9 242 160 145 217 156 234 214 175  17 197 171 131  10  10  41
 215 179 186   3  39 246 139 126 231 219 183 177  35  94 188 184
  52 134 189 145 213 186 235  74 203  37 253  69 212 126   2  46

Further Reading

< Sleep [Index] Unicode >