17 Pseudo-Random Numbers

 17.1 Setting the seed (replacements for G05CBF and G05CCF)
 17.2 Data type of the random numbers
 17.3 Replacements for G05CAF and G05DAF
 17.4 Replacement for G05DBF
 17.5 Replacements for G05DDF, G05DRF, and G05FFF

The routines for creating pseudo-random numbers in this library all have a period of 226 and 6–7 digits accuracy. They are based upon code by Ahrens, Dieter, & Grube. They use a multiplicative congruential generator which is certainly not the state of the art and may not be suitable for critical or sophisticated use.

The routines are:

17.1 Setting the seed (replacements for G05CBF and G05CCF)

Before any random numbers can be selected, a seed must be set using PDA_RNSED. The integer seed should satisfy the relationship

seed = 4 k + 1

where k is a non-negative integer. A fixed seed gives rise to a reproducible sequence of pseudo-random numbers.

For a non-repeatable sequence, there is no equivalent to NAG routine G05CCF because the system clock used to create the seed is not accessible portably in Fortran, and PDA is independent of other libraries. However, the following code has the desired effect.

        SEED = TICKS + PID
        SEED = MOD( SEED, VAL__MAXI / 4 ) * 4 + 1
        SEED = MOD( SEED, 2**28 )

PSX_TIME returns the time in units of clock ticks since some arbitrary time. See SUN/121 for more details and linking instructions. The above code also permits storage of the chosen seed.

17.2 Data type of the random numbers

There is a major difference between the PDA random-number routines and those provided in the standard NAG library: in general the former are single-precision functions, whereas the latter are double precision. However, PDA_RNPOI and the corresponding G05DRF are both integer functions.

17.3 Replacements for G05CAF and G05DAF

Like G05CAF, PDA_RAND has a dummy argument demanded by the Fortran standard. It is convenient to set it to zero. Here is an example where two random numbers are drawn from a uniform distribution between 0 and 1. In this example a fixed seed is used, but you could use the computer’s clock to create a random seed (see Section 17.1).

        REAL PDA_RAND, VALUES( 2 )
  *  Use a fixed seed of 1.
        SEED = 1
  *  Obtain two random numbers from a uniform distribution between 0
  *  and 1.
        VALUE( 1 ) = PDA_RAND( 0.0 )
        VALUE( 2 ) = PDA_RAND( 0.0 )

The EXTERNAL statement is recommended, although in many cases it will be unnecessary. To obtain in the range [a,b] as provided by G05DAF, merely apply the following relationship.

random value = (ba)  PDA_RAND(0.0) + a

17.4 Replacement for G05DBF

PDA_RNEXP is only a partial replacement for G05DBF in that it computes pseudo-random numbers from ex, whereas G05DBF uses the function 1 aex/a. Thus its argument is also a dummy mandated by the Fortran standard.

17.5 Replacements for G05DDF, G05DRF, and G05FFF

The following code shows the remaining three routines in action.

  *  Use a fixed seed of 1001.
        SEED = 1001
  *  Obtain a random number from a Normal distribution of mean 4.2 and
  *  standard deviation 0.15
        VALUE( 1 ) = PDA_RNNOR( 4.2, 0.15 )
  *  Obtain a random number from a Poisson distribution of mean 3.4.
        VALUE( 2 ) = PDA_RNPOI( 3.4 )
  *  Obtain a random number from a Gamma-function distribution of mean
  *  1.2.
        VALUE( 3 ) = PDA_RNGAM( 1.2 )

Apart from the change of data type, calls to G05DDF can be replaced with PDA_RNNOR using the same arguments. PDA_RNPOI is in effect a renamed G05DRF.

PDA_RNGAM only has one argument—the mean—of the Gamma function, whereas G05FFF has a second scaling parameter similar in role to the a argument of G05DBF. G05FFF also generates a vector of pseudo-random numbers.