RANDOM - Pseudo-random number generator

#include <spral_random.h> /* or <spral.h> for all packages */

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.

Random Seed

The random number generator’s state is stored in the variable state that is common to all calls. Before its first use, state should be initialized by the user to an initial seed value. For convienience the following preprocessor macro is defined.

SPRAL_RANDOM_INITIAL_SEED 486502

Convience macro for initial random seed.

Example of use:

int state = SPRAL_RANDOM_INITIAL_SEED;

The user may change the seed at any time, for example to restore a previous value. The same sequence of calls following the (re-)initialization of seed to the same value will produce the same sequence of pseudo-random values.

Routines

double spral_random_real(int *state, bool positive)

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

Parameters:

• state – current state of RNG.

• positive – if true, sample from \((0,1)\); otherwise, sample from \((-1,1)\).

Returns:

Sampled value.

int spral_random_integer(int *state, int n)

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

Parameters:

• state – current state of the RNG.

• n – largest value in range to be sampled.

Returns:

Sampled value.

int64_t spral_random_long(int *state, int64_t n)

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

Parameters:

• state – current state of the RNG.

• n – largest value in range to be sampled.

Returns:

Sampled value.

bool spral_random_logical(int *state)

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

Parameters:

• state – current state of the RNG.

Returns:

Sampled value.

Example

The following code:

/* examples/C/random.c - Example code for SPRAL_RANDOM package */
#include "spral.h"
#include <limits.h>
#include <stdint.h>
#include <stdio.h>

#define bool2str(VAL) ( (VAL) ? "true" : "false" )

int main(void) {
   int state = SPRAL_RANDOM_INITIAL_SEED;

   // Store initial random seed so we can reuse it later
   int initial_seed = state;

   // Generate some random values
   printf("\nSome random values\n");
   printf("Sample Unif(-1,1)               = %16.12f\n",
         spral_random_real(&state, false));
   printf("Sample Unif(0,1)                = %16.12f\n",
         spral_random_real(&state, true));
   printf("Sample Unif(1, ..., 20)         = %16d\n",
         spral_random_integer(&state, 20));
   printf("Sample Unif(1, ..., 20*INT_MAX) = %16ld\n",
         spral_random_long(&state,((int64_t)20)*INT_MAX));
   printf("Sample B(1,0.5)                 = %16s\n",
         bool2str(spral_random_logical(&state)));

   // Restore initial seed
   state = initial_seed;

   // Generate the same random values
   printf("\nThe same random values again\n");
   printf("Sample Unif(-1,1)               = %16.12f\n",
         spral_random_real(&state, false));
   printf("Sample Unif(0,1)                = %16.12f\n",
         spral_random_real(&state, true));
   printf("Sample Unif(1, ..., 20)         = %16d\n",
         spral_random_integer(&state, 20));
   printf("Sample Unif(1, ..., 20*INT_MAX) = %16ld\n",
         spral_random_long(&state, ((int64_t)20)*INT_MAX));
   printf("Sample B(1,0.5)                 = %16s\n",
         bool2str(spral_random_logical(&state)));

   /* Success */
   return 0;
}

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*INT_MAX) =      33572664025
Sample B(1,0.5)                 =            false

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*INT_MAX) =      33572664025
Sample B(1,0.5)                 =            false

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 variable state stores the current value of \(X_n\).

In spral_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 spral_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 spral_random_logical()

Returns the value of the Fortran expression

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