RANDOM: Pseudo-random Number Generator

Major version history

2014-04-07 Version 1.0.0

Initial release

2.1 Installation

Please see the SPRAL install documentation.

2.2 Usage overview

2.2.1 Calling sequences

Access to the package requires a USE statement

   use spral_random

The following procedures are available to the user:

• random_real() generates a real uniformally at random from the interval $\left(-1,1\right)$ or $\left(0,1\right)$.
• random_integer() generates an integer unformally at random from the interval $\left[1,\dots ,n\right]$.
• random_logical() generates a random boolean value.
• random_get_seed() returns the current random seed to the user.
• random_set_seed() sets the current random seed.

2.2.2 Derived types

The user must employ the derived type random_state defined by the package to store the state. The following pseudo-code illustrates how to declare a scalar of this type.

The components of random_state are not available to the user, but may be examined and altered through calls to random_get_seed() and random_set_seed() respectively.

2.2.3 Optional arguments

We use square brackets [ ] to indicate optional arguments. In each call, optional arguments appear last in argument list. Since we reserve the right to add additional optional arguments in future releases of the code, we strongly recommend that all optional arguments be called by keyword, not by position.

2.2.4 Integer, real and package types

INTEGER denotes default INTEGER and INTEGER(long) denotes INTEGER(kind=selected_int_kind(18)).

REAL denotes double precision real. We also use the term package type to mean the same.

2.3 Random Generation Subroutines

2.3.1 random_real()

To generate a real uniformly at random from the interval $\left(-1,1\right)$ or $\left(0,1\right)$,
sample = random_real(state[, positive])

The function returns a scalar value sample of type REAL that is a sample from $Unif\left(0,1\right)$ if positive is present with the value .true. or from $Unif\left(-1,1\right)$ otherwise.

state

is an INTENT(INOUT) scalar of type random_state. It contains the current state of the random number generator, and need not be initialized by the user.

positive

is an optional INTENT(IN) scalar of type LOGICAL. If present with the value .true., the sample will be returned from the interval $\left(0,1\right)$. Otherwise, the sample will be returned from the interval $\left(-1,1\right)$.

2.3.2 random_integer()

To generate an integer uniformly at random from the interval $\left[1,n\right]$,
sample = random_integer(state, n)

The function returns a scalar value sample of type INTEGER that is a sample from $Unif\left(1,\dots ,n\right)$.

state

is an INTENT(INOUT) scalar of type random_state. It contains the current state of the random number generator, and need not be initialized by the user.

n

is an optional INTENT(IN) scalar of type INTEGER. It specifies the maximum value that sample may take.

2.3.3 random_logical()

To generate a random logical with equal probability of being .true. or .false.,
sample = random_logical(state)

The function returns a scalar value sample of type LOGICAL that is a sample from $B\left(1,0.5\right)$, i.e. has an equal probability of being .true. or .false..

state

is an INTENT(INOUT) scalar of type random_state. It contains the current state of the random number generator, and need not be initialized by the user.

2.4 Seed control subroutines

2.4.1 random_get_seed()

To obtain the current random seed,
seed = random_get_seed(state)

The function returns a scalar value seed of type INTEGER that is the current 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().

state

is an INTENT(IN) scalar of type random_state. It contains the current state of the random number generator, and need not be initialized by the user.

2.4.2 random_set_seed()

To set the current random seed,
call random_get_seed(state, seed)

state

is an INTENT(INOUT) scalar of type random_state to be reinitialized with the specified seed.

seed

is an INTENT(IN) scalar of type INTEGER that specifies the new seed.

2.5 Method

2.5.1 Pseudo-random number generation

We use a linear congruential generator of the following form:

$$X_{n+1}=\left(a{X}_n+c\right)modm$$

with the following constants

$$\begin{array}{rcll}a& =& 1103515245,& \\ c& =& 12345,& \\ m& =& 2^{31}.& \end{array}$$

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, $X_0=486502$.

2.5.2 random_real()

If positive is present with value .true., a sample from $Unif\left(0,1\right)$ is generated as

$$\frac{\mathtt{\text{real}}\left(X_n\right)}{\mathtt{\text{real}}\left(m\right)},$$

otherwise, a sample from $Unif\left(-1,1\right)$ is generated as

$$1-\frac{\mathtt{\text{real}}\left(2X_n\right)}{\mathtt{\text{real}}\left(m\right)}.$$

2.5.3 random_int()

A random sample from the discrete distribution $Unif\left(1,\dots ,n\right)$ is generated as

$$\mathtt{\text{int}}\left(X_n\frac{\mathtt{\text{real}}\left(n\right)}{\mathtt{\text{real}}\left(m\right)}\right)+1.$$

2.5.4 random_logical()

A random logical value is generated by evaluating the expression

$$\left(1.eq.random\text{\_}integer\left(state,2\right)\right).$$

2.6 Example

The following example code:

Produces the following output: