Random Number Generation Utility Routines

All of the random number generators in this chapter depend on the generation of uniform (0, 1) numbers, which is done by the routine RNUN, or by its function analog RNUNF. These basic generators use either a multiplicative congruential method or a generalized feedback shift register (GFSR) method, or the Mersenne Twister method to yield a subsequence of a fixed cyclic sequence. The beginning of the subsequence is determined by the seed.

The utility routines for the random number generators allow the user to select the type of the generator (or to determine the type of the generator being used) and to set or retrieve the seed.

Selection of the Type of the Generator

The uniform pseudorandom number generators use a multiplicative congruential method, with or without shuffling or a GFSR method, or the Mersenne Twister method. Routine RNOPT determines which method is used; and in the case of a multiplicative congruential method, it determines the value of the multiplier and whether or not to use shuffling. The description of RNUN may provide some guidance in the choice of the form of the generator. If no selection is made explicitly, the generators use the multiplier 16807 without shuffling. This form of the generator has been in use for some time (see Lewis, Goodman, and Miller, 1969). This is the generator formerly known as GGUBS in the IMSL Library. It is the “minimal standard generator” discussed by Park and Miller (1988).

Both of the Mersenne Twister generators have a period of 219937 ‑ 1 and a 624‑dimensional equidistribution property. See Matsumoto et al. 1998 for details.

The IMSL Mersenne Twister generators are derived from code copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, All rights reserved. It is subject to the following notice:

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

The IMSL 32‑bit Mersenne Twister generator is based on the Matsumoto and Nishimura code ‘mt19937ar’ and the 64‑bit code is based on ‘mt19937‑64’.

The selection of the type of generator is made by calling the routine RNOPT, choosing one of nine different options.

RNOPT

CALL RNOPT (IOPT)

The argument is:

IOPT — The indicator of the generator. (Input)

The random number generator is either a multiplicative congruential generator with modulus 231 ‑ 1 or a GFSR generator or Mersenne Twister. IOPT is used to choose the multiplier and whether or not shuffling is done, or is used to choose the GFSR method, or is used to choose the Mersenne Twister generator.

If no selection is made explicitly, a multiplicative generator (with multiplier 16807) is used (equivalent to IOPT = 1).

The type of generator being used can be determined by calling the routine RNOPG.

RNOPG

CALL RNOPG (IOPT)

The argument is:

IOPT, which is an output variable in RNOPG.

Setting the Seed

Before using any of the random number generators, the generator must be initialized by selecting a seed, or starting value. The user does not have to do this, but it can done by calling the routine RNSET. If the user does not select a seed, one is generated using the system clock. A seed needs to be selected only once in a program unless there is some desire to maintain two separate streams of random numbers.

RNSET

CALL RNSET (ISEED)

The argument is:

ISEED — The seed of the random number generator. (Input)

ISEED must be in the range (0, 2147483646). If ISEED is zero (or if RNSET is not called before the generation of random numbers begins), a value is computed using the system clock; and, hence, the results of programs using the IMSL random number generators will be different at different times.

Stopping and Restarting Simulations and Controlling More Than One Stream of Random Numbers

For most purposes, even if several simulations are being run in the same program or if the simulation is being conducted in blocks, it is best to use the sequence of uniform random deviates just as produced by RNUN or RNUNF without concern for from where in the underlying cyclic sequence the numbers are coming.

If, however, the simulations are being conducted incrementally or if simulations are being run in parallel, it may be necessary to exercise more control over the sequence. The routines that are used in stopping and restarting simulations come in pairs, one to get the current value and one to set the value. The argument for each pair is the same within the pair; it is output in one case and input in the other. (RNSET is an exception since it is often used at the beginning of a simulation before any seed is ever set.) If a nonshuffled form of the multiplicative congruential generators is used (that is IOPT in RNOPT is 1, 3, or 5), the only thing that must be controlled is the seed, so in this case the only routines needed are

RNSET

Initializes the seed used in the generators

RNGET

Retrieves the current value of the seed

The usages are:

CALL RNSET (ISEED) (ISEED is input)

CALL RNGET (ISEED) (ISEED is output)

ISEED is an integer in the range 1 to 2147483646 (except, as noted above, it can be input to RNSET as 0 to indicate that the system clock is used to generate a seed).

If a shuffled generator, the GFSR generator, or a Mersenne Twister generator is used, in addition to controlling the seed as described above, another array must be maintained if the user wishes to stop and restart the simulation. It is a floating‑point array for the shuffled generators and an integer array for the GFSR generator and Mersenne Twister generator. The routines are:

RNSES

Initializes the table used in the shuffled generators.

RNGES

Retrieves the current table used in the shuffled generators.

RNSEF

Initializes the table used in the GFSR generator.

RNGEF

Retrieves the current table used in the GFSR generator.

RNIN32

Initializes the table used in the 32‑bit Mersenne Twister generator using an array.

RNSE32

Sets the current table used in the 32‑bit Mersenne Twister generator.

RNGE32

Retrieves the current table used in the 32‑bit Mersenne Twister generator.

RNIN64

Initializes the table used in the 64‑bit Mersenne Twister generator using an array.

RNSE64

Sets the current table used in the 64‑bit Mersenne Twister generator.

RNGE64

Retrieves the current table used in the 64‑bit Mersenne Twister generator.

There are different tables used in the single and double precision versions of the shuffled generators, so RNGES and RNSES can also be used in double precision.

The usages are:

CALL RNGES (TABLE) (TABLE is output.)

CALL RNSES (TABLE) (TABLE is input.)

CALL RNGEF (IARRAY) (IARRAY is output.)

CALL RNSEF (IARRAY) (IARRAY is input.)

CALL RNGE32 (MTABLE32) (MTABLE is output.)

CALL RNSE32 (MTABLE32) (MTABLE is input.)

CALL RNGE64 (MTABLE64) (MTABLE is output.)

CALL RNSE64 (MTABLE64) (MTABLE is input.)

The arguments are:

TABLE — Array of length 128 used in the shuffled generators.

IARRAY — Array of length 1565 used in the GFSR generators.

MTABLE32 — Array of length 625 used in the 32‑bit Mersenne Twister generators.

MTABLE64 — Array of length 313 used in the 64‑bit Mersenne Twister generators.

The values in both TABLE and IARRAY are initialized by the IMSL random number generators. The values are all positive in both arrays except if the user wishes to reinitialize the array, in which case the first element of the array is input as a nonpositive value. (Usually, one should avoid reinitializing these arrays, but it might be necessary sometimes in restarting a simulation.) If the first element of TABLE or IARRAY is set to a nonpositive value on the call to RNSES or RNSEF, on the next invocation of a routine to generate random numbers using shuffling (if RNSES) or a GFSR method (if RNSEF), the appropriate array will be reinitialized.

In addition to controlling separate streams of random numbers, sometimes it is desirable to insure from the beginning that two streams do not overlap. This can be done with the congruential generators that do not do shuffling by using RNISD to get a seed that will generate random numbers beginning 100,000 numbers farther along.

RNISD

CALL RNISD (ISEED1, ISEED2)

The arguments are:

ISEED1 — The seed that yields the first stream. (Input)

ISEED2 — The seed that yields a stream beginning 100,000 numbers beyond the stream that begins with ISEED1. (Output)

Given a seed, ISEED1, RNISD determines another seed, ISEED2, such that if one of the IMSL multiplicative congruential generators, using no shuffling, went through 100,000 generations starting with ISEED1, the next number in that sequence would be the first number in the sequence that begins with the seed ISEED2. This can be described more simply by stating that RN1 and RN2 in the following sequence of FORTRAN are assigned the same values.

       CALL RNISD(ISEED1, ISEED2)
       CALL RNSET(ISEED1)
       DO 10 I = 1, 100000
          RN1 = RNUNF()
    10 CONTINUE
       RN1 = RNUNF()
       CALL RNSET(ISEED2)
       RN2 = RNUNF()

To obtain seeds that generate sequences with beginning values separated by 200,000 numbers, call RNISD twice:

CALL RNISD(ISEED1, ISEED2)
CALL RNISD(ISEED2, ISEED2)

Note that RNISD works only when a multiplicative congruential generator without shuffling is used. This means that either the routine RNOPT has not been called at all or that it has been last called with IOPT taking a value of 1, 3, or 5.

For many of the IMSL generators for nonuniform distributions that do not use the inverse CDF method, the distance between the sequences generated starting with ISEED1 and starting with ISEED2 may be less than 100,000. This is because the nonuniform generators that use other techniques may require more than one uniform deviate for each output deviate.

The reason that one may want two seeds that generate sequences a known distance apart is for blocking Monte Carlo experiments or for running parallel streams.

Examples

Example 1: Selecting the Type of Generator and Stopping and Restarting the Simulations

In this example, three separate simulation streams are used, each with a different form of the generator. Each stream is stopped and restarted. (Although this example is obviously an artificial one, there may be reasons for maintaining separate streams and stopping and restarting them because of the nature of the usage of the random numbers coming from the separate streams.)

      USE IMSL_LIBRARIES
 
      IMPLICIT   NONE
      INTEGER    I, IARRAY(1565), ISEED1, ISEED2, ISEED7, NOUT, NR
      REAL       R(5), TABLE(128)
!
      CALL UMACH (2, NOUT)
      NR     = 5
      ISEED1 = 123457
      ISEED2 = 123457
      ISEED7 = 123457
!                                 Begin first stream, IOPT = 1 (by
!                                 default)
      CALL RNSET (ISEED1)
      CALL RNUN (R)
      CALL RNGET (ISEED1)
      WRITE (NOUT,99997) (R(I),I=1,NR), ISEED1
!                                 Begin second stream, IOPT = 2
      CALL RNOPT (2)
      CALL RNSET (ISEED2)
      CALL RNUN (R)
      CALL RNGET (ISEED2)
      CALL RNGES (TABLE)
      WRITE (NOUT,99998) (R(I),I=1,NR), ISEED2
!                                 Begin third stream, IOPT = 7
      CALL RNOPT (7)
      CALL RNSET (ISEED7)
      CALL RNUN (R)
      CALL RNGET (ISEED7)
      CALL RNGEF (IARRAY)
      WRITE (NOUT,99999) (R(I),I=1,NR), ISEED7
!                                 Reinitialize seed
!                                 Resume first stream
      CALL RNOPT (1)
      CALL RNSET (ISEED1)
      CALL RNUN (R)
      CALL RNGET (ISEED1)
      WRITE (NOUT,99997) (R(I),I=1,NR), ISEED1
!                                 Reinitialize seed and table for
!                                 shuffling
!                                 Resume second stream
      CALL RNOPT (2)
      CALL RNSET (ISEED2)
      CALL RNSES (TABLE)
      CALL RNUN (R)
      CALL RNGET (ISEED2)
      WRITE (NOUT,99998) (R(I),I=1,NR), ISEED2
!                                 Reinitialize seed and table for GFSR
!                                 Resume third stream
      CALL RNOPT (7)
      CALL RNSET (ISEED7)
      CALL RNSEF (IARRAY)
      CALL RNUN (R)
      CALL RNGET (ISEED7)
      WRITE (NOUT,99999) (R(I),I=1,NR), ISEED7
!
99997 FORMAT (/, '  First stream  ', 5F8.4, /, '  Output seed = ', &
            I11)
99998 FORMAT (/, '  Second stream ', 5F8.4, /, '  Output seed = ', &
            I11)
99999 FORMAT (/, '  Third stream  ', 5F8.4, /, '  Output seed = ', &
            I11)
!
      END

Output

First stream    0.9662  0.2607  0.7663  0.5693  0.8448
Output seed =   1814256879
Second stream   0.7095  0.1861  0.4794  0.6038  0.3790
Output seed =   1965912801
Third stream    0.3914  0.0263  0.7622  0.0281  0.8997
Output seed =   1932158269
First stream    0.0443  0.9872  0.6014  0.8964  0.3809
Output seed =   817878095
Second stream   0.2557  0.4788  0.2258  0.3455  0.5811
Output seed =   2108806573
Third stream    0.7519  0.5084  0.9070  0.0910  0.6917
Output seed =   1485334679

Example 2: Determining Seeds for Separate Streams

In this example, RNISD is used to determine seeds for 4 separate streams, each 200,000 numbers apart, for a multiplicative congruential generator without shuffling. (Since RNOPT is not invoked to select a generator, the multiplier is 16807.) To get each seed requires two invocations of RNISD. All of the streams are non‑overlapping, since the period of the underlying generator is 2,147,483,646.

      USE UMACH_INT
      USE RNISD_INT
 
      IMPLICIT   NONE
      INTEGER    ISEED1, ISEED2, ISEED3, ISEED4, NOUT
!
      CALL UMACH (2, NOUT)
      ISEED1 = 123457
      CALL RNISD (ISEED1, ISEED2)
      CALL RNISD (ISEED2, ISEED2)
      CALL RNISD (ISEED2, ISEED3)
      CALL RNISD (ISEED3, ISEED3)
      CALL RNISD (ISEED3, ISEED4)
      CALL RNISD (ISEED4, ISEED4)
      WRITE (NOUT,99999) ISEED1, ISEED2, ISEED3, ISEED4
!
99999 FORMAT ('  Seeds for four separate streams: ', /, '  ', 4I11)
!
      END

Output

Seeds for four separate streams:
     123457 2016130173   85016329  979156171