`spral_random` - Pseudo-random number generator

Purpose

This package generates pseudo-random numbers using a linear congruential generator. It should generate the same random numbers using any standards compliant Fortran compiler on any architecture so long as the default integer and real kinds are the same.

The seed can optionally be observed or specified by the user. Otherwise a default seed of 486502 is used.

Version history

2016-09-08 Version 1.1.0: Add support for long integers
2014-04-07 Version 1.0.0: Initial release

Routines

Random Number Generation

function random_real(state[, positive])

Return a real uniformly at random from the interval \((-1,1)\) (positive =.true.) or \((0,1)\) (positive =.false.).

Parameters:: state [random_state ,inout] :: current state of RNG.
Options:: positive [logical ,in,default=.false.] :: if .true., sample from \((0,1)\); otherwise, sample from \((-1,1)\).
Return:: random_real [real] :: Sampled value.

function random_integer(state, n)

Return an integer uniformly at random from the interval \([1,n]\).

Parameters:

state [random_state ,inout] :: current state of the RNG.
n [integer(kind),in] :: largest value in range to be sampled. kind may be either default or long integer (return type will match).

Return:

random_integer [integer(kind)] :: Sampled value.

function random_logical(state)

Return a logical with equal probability of being .true. or .false..

Parameters:: state [random_state ,inout] :: current state of the RNG.
Return:: random_logical [logical] :: Sampled value.

Get/Set Random Seed

function random_get_seed(state)

Return the current random seed stored in state.

The stream of random numbers generated after this call can be reproduced through the same sequence of calls after seed has been passed to random_set_seed().

Parameters:: state [random_state ,in] :: state variable to extract state from
Return:: random_get_seed :: current random seed

subroutine random_set_seed(state, seed)

Set the random seed stored in the state variable

Parameters:

state [random_state ,inout] :: state variable to set seed for.
seed [integer ,in] :: new seed.

Data Types

type random_state: State of the random number generator. Compontents are not available to the user, but may be examined and altered through calls to random_get_seed() and random_set_seed() respectively.

Example

The following code:

! examples/Fortran/random.f90 - Example code for SPRAL_RANDOM package
program random_example
   use spral_random
   implicit none

   integer, parameter :: long = selected_int_kind(18)

   type(random_state) :: state
   integer :: seed

   ! Store initial random seed so we can reuse it later
   seed = random_get_seed(state)

   ! Generate some random values
   write(*,"(a)") "Some random values"
   write(*,"(a,f16.12)") "Sample Unif(-1,1)               = ", &
      random_real(state)
   write(*,"(a,f16.12)") "Sample Unif(0,1)                = ", &
      random_real(state, positive=.true.)
   write(*,"(a,i16)") "Sample Unif(1, ..., 20)         = ", &
      random_integer(state, 20)
   write(*,"(a,i16)") "Sample Unif(1, ..., 20*huge(0)) = ", &
      random_integer(state, 20_long*huge(0))
   write(*,"(a,l16)") "Sample B(1,0.5)                 = ", &
      random_logical(state)

   ! Restore initial seed
   call random_set_seed(state, seed)

   ! Generate the same random values
   write(*,"(/a)") "The same random values again"
   write(*,"(a,f16.12)") "Sample Unif(-1,1)               = ", &
      random_real(state)
   write(*,"(a,f16.12)") "Sample Unif(0,1)                = ", &
      random_real(state, positive=.true.)
   write(*,"(a,i16)") "Sample Unif(1, ..., 20)         = ", &
      random_integer(state, 20)
   write(*,"(a,i16)") "Sample Unif(1, ..., 20*huge(0)) = ", &
      random_integer(state, 20_long*huge(0))
   write(*,"(a,l16)") "Sample B(1,0.5)                 = ", &
      random_logical(state)

end program random_example

Produces the following output:

Some random values
Sample Unif(-1,1)               =   0.951878630556
Sample Unif(0,1)                =   0.395779648796
Sample Unif(1, ..., 20)         =                3
Sample Unif(1, ..., 20*huge(0)) =      33572664025
Sample B(1,0.5)                 =                F

The same random values again
Sample Unif(-1,1)               =   0.951878630556
Sample Unif(0,1)                =   0.395779648796
Sample Unif(1, ..., 20)         =                3
Sample Unif(1, ..., 20*huge(0)) =      33572664025
Sample B(1,0.5)                 =                F

Method

We use a linear congruential generator of the following form:

\[X_{n+1} = (aX_n + c)\quad \mathrm{mod}\; m\]

with the following constants

\[\begin{split}a = 1103515245, \\ c = 12345, \\ m = 2^{31}.\end{split}\]

According to Wikipedia, this is the same as used in glibc.

The LCG is evolved before each sample is taken, and the sample is based on the new value.

The routines random_get_seed() and random_set_seed() allow the user to get and set the current value of \(X_n\). The default seed is \(X_0 = 486502\).

In `random_real()`

Samples from \(\mathrm{Unif}(0,1)\) are generated as

\[\frac{\mathrm{real}(X_n)}{\mathrm{real}(m)},\]

and samples from \(\mathrm{Unif}(-1,1)\) are generated as

\[1−\frac{\mathrm{real}(2X_n)}{\mathrm{real}(m)}.\]

In `random_integer()`

Samples from \(\mathrm{Unif}(1,\ldots,n)\) are generated as

\[\mathrm{int}\left(X_n\frac{\mathrm{real}(n)}{\mathrm{real}(m)}\right) + 1\]

In `random_logical()`

Returns the value of the Fortran expression

(1 .eq. random_integer(state,2))

spral_random - Pseudo-random number generator