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, " "))', r
end program main

Compile and run the example with:

$ gfortran10 -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:

$ gfortran10 -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

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

    ! Print some random numbers between 0 and 10,000.
    do i = 1, 10
        call random_number(r)
        print '(i0)', int(r * 10000)
    end do
contains
    function urandom_seed(n, stat)
        !! Returns a seed array filled with random values from `/dev/urandom`.
        integer, intent(in)            :: n
        integer, intent(out), optional :: stat
        integer                        :: urandom_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) urandom_seed

        close (fu)
    end function urandom_seed
end program main

The example outputs ten random integers:

$ gfortran10 -o urandom urandom.f90
$ ./urandom
985
5948
1279
556
3512
1125
8403
1776
6224
6034

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
            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
    implicit none
    use :: openssl
    integer          :: i, rc
    character(len=1) :: bytes(10)

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

    do i = 1, size(bytes)
        print *, ichar(bytes(i))
    end do
end program main

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

$ gfortran10 -o crypto crypto.f90 -lcrypto

The example application just prints the random bytes [0 … 255] to console:

$ ./crypto
          78
         231
         118
          40
          28
         239
         100
          21
           9
         226

References

< Sleep [Index] Unicode >