In the following discussions, the phrases “random numbers,” “random deviates,” “deviates,” and “variates” are used interchangeably. The phrase “pseudorandom” is sometimes used to emphasize that the numbers generated are not really “random” since they result from a deterministic process. The usefulness of pseudorandom numbers derives from the similarity, in a statistical sense, of samples of the pseudorandom numbers to samples of observations from the specified distributions. In short, while the pseudorandom numbers are completely deterministic and repeatable, they simulate the realizations of independent and identically distributed random variables.

The random number generators in this chapter use either a multiplicative congruential method, or a generalized feedback shift register (GFSR) method, or a Mersenne Twister method. The selection of the type of generator is made by calling the routine RNOPT. If no selection is made explicitly, a multiplicative generator (with multiplier 16807) is used. Whatever distribution is being simulated, uniform (0, 1) numbers are first generated and then transformed if necessary. The generation of the uniform (0, 1) numbers is done by the routine RNUN, or by its function analog RNUNF. These routines are portable in the sense that, given the same seed and for a given type of generator, they produce the same sequence in all computer/compiler environments. There are many other issues that must be considered in developing programs for the methods described below (see Gentle 1981 and 1990).

The form of the multiplicative congruential generators is

Each xi is then scaled into the unit interval (0, 1). If the multiplier, c, is a primitive root modulo 231 ‑ 1 (which is a prime), then the generator will have maximal period of 231 ‑ 2. There are several other considerations, however. The lattice structure induced by congruential generators (see Marsaglia 1968) can be assessed by the lattice test of Marsaglia (1972) or the spectral test of Coveyou and MacPherson (1967) (see also Knuth 1981, pages 89‑113). Also, empirical studies, such as by Fishman and Moore (1982 and 1986), indicate that different values of multipliers, all of which perform well under the lattice test and the spectral test, may yield quite different performances where the criterion is similarity of samples generated to samples from a true uniform distribution.

There are three possible choices for c in the IMSL generators: 16807 (which is 75), 397204094 (which is 2 ⋅ 72 ⋅ 4053103), and 950706376 (which is 23 ⋅ 118838297). The selection is made by the routine RNOPT. The choice of 16807 will result in the fastest execution time (see Gentle 1981), but Fishman and Moore’s studies would seem to indicate that the performance of 950706376 is best among these three choices. If no selection is made explicitly, the routines use the multiplier 16807, which has been in use for some time (Lewis, Goodman, and Miller 1969). It is the “minimal standard generator” discussed by Park and Miller (1988).

The user can also select a shuffled version of the multiplicative congruential generators using RNOPT. The shuffled generators use a scheme due to Learmonth and Lewis (1973a). In this scheme, a table is filled with the first 128 uniform (0, 1) numbers resulting from the simple multiplicative congruential generator. Then, for each xi from the simple generator, the low‑order bits of xi are used to select a random integer, j, from 1 to 128. The j‑th entry in the table is then delivered as the random number; and xi, after being scaled into the unit interval, is inserted into the j‑th position in the table.

The GFSR generator uses the recursion Xt = Xt−1563 ⊕ Xt−96. This generator, which is different from earlier GFSR generators, was proposed by Fushimi (1990), who discusses the theory behind the generator and reports on several empirical tests of it. Background discussions on this type of generator can be found in Kennedy and Gentle (1980), pages 150‑162.

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 seed of the generator can be set in RNSET and can be retrieved by RNGET. Prior to invoking any generator in this chapter , the user can call RNSET to initialize the seed, which is an integer variable taking a value between 1 and 2147483646. If it is not initialized by RNSET, a random seed is obtained from the system clock. Once it is initialized, the seed need not be set again. The seed is updated and passed from one routine to another by means of a named COMMON block, R2NCOM.

If the user wishes to restart a simulation, RNGET can be used to obtain the final seed value of one run to be used as the starting value in a subsequent run. Also, if two random number streams are desired in one run, RNSET and RNGET can be used before and after the invocations of the generators in each stream. If a shuffled generator or the GFSR generator is used, in addition to resetting the seed, the user must also reset some values in a table. For the shuffled generators, this is done using the routines RNGES and RNSES and for the GFSR generator, the table is retrieved and set by the routines RNGEF and RNSEF. The tables for the shuffled generators are separate for single and double precision; so, if precisions are mixed in a program, it is necessary to manage each precision separately for the shuffled generators.

The generation of the uniform (0,1) numbers is done by the routine RNUN or by its function analog RNUNF. The particular generator selected in RNOPT, that is, the value of the multiplier and whether shuffling is done or whether the GFSR generator is used, affects the speed of RNUN and RNUNF. The smaller multiplier (16807, selected by IOPT = 1) is faster than the other multipliers. The multiplicative congruential generators that do not shuffle are faster than the ones that do. The GFSR generator is roughly as fast as the fastest multiplicative congruential generator, but the initialization for it (required only on the first invocation) takes longer than the generation of thousands of uniform random numbers. Precise statements of relative speeds depend on the computing system.

Whether RNUN or RNUNF is used also has an effect on the speed due to the overhead in invoking an external routine, or due to the program’s inability to optimize computations by holding some operands in registers. This effect, of course, may be different in different environments. On an array processor or other computers with pipelined instructions, RNUN is likely to be considerably faster than RNUNF when several random numbers are to be generated at one time. In the case of array processors, the multiplicative congruential generators in RNUN are coded to generate subsequences in larger blocks (see Gentle 1990).

The basic uniform (0, 1) generators RNUN or RNUNF are used by all other routines in this chapter. If, for some reason, the user would prefer a different basic uniform generator, routines named “RNUN” and “RNUNF” can be written so that they include the named COMMON, through which the seed is passed, and that calls the user’s custom generator. The named COMMON is

The user’s “RNUN” and “RNUNF” can pass the seed through any of the variables, but the routines RNSET and RNGET expect the seed to be in ICEED. (The user should not expect to use any utility routines other than RNSET and RNGET if customized versions of RNUN or RNUNF are used.) The double precision versions of the nonuniform generators, such as DRNBET and DRNGAM (RNGAM), use the double precision versions of the uniform generators, DRNUN (RNUN) and DRNUNF (RNUNF), so to use the double precision nonuniform generators with customized uniform generators, the user would supply routines to replace DRNUN and DRNUNF.

The nonuniform generators use a variety of transformation procedures. All of the transformations used are exact (mathematically). The most straightforward transformation is the inverse CDF technique, but it is often less efficient than others involving acceptance/rejection and mixtures. See Kennedy and Gentle (1980) for discussion of these and other techniques.

Many of the nonuniform generators in this chapter use different algorithms depending on the values of the parameters of the distributions. This is particularly true of the generators for discrete distributions. Schmeiser (1983) gives an overview of techniques for generating deviates from discrete distributions.

Although, as noted above, the uniform generators yield the same sequences on different computers, because of rounding, the nonuniform generators that use acceptance/rejection may occasionally produce different sequences on different computer/compiler environments.

Although the generators for nonuniform distributions use fast algorithms, if a very large number of deviates from a fixed distribution are to be generated, it might be worthwhile to consider a table sampling method, as implemented in the routines RNGDA, RNGDS, RNGDT, RNGCS, and RNGCT. After an initialization stage, which may take some time, the actual generation may proceed very fast.

For those generators, such as RNCHY and RNNOR, that use the inverse CDF technique, it is possible to generate any set of order statistics directly by use of a customized uniform generator, as discussed above, by generating order statistics in a custom “RNUN” or “RNUNF”. In some routines that employ an inverse CDF technique, such as RNEXP and RNWIB, instead of directly using the uniform (0, 1) deviate u from RNUN, the uniform (0, 1) deviate 1 ‑ u is used. In such routines the i-th order statistic from the uniform will yield the (n + 1 ‑ i)-th order statistic from the nonuniform distribution.

A similar technique can be used to get antithetic variates. For each uniform deviate u, a second deviate 1 ‑ u could be produced by a custom “RNUN” or “RNUNF”. As with order statistics, this technique would only be reasonable for routines that use the inverse CDF technique.

Extensive empirical tests of some of the uniform random number generators available in RNUN and RNUNF are reported by Fishman and Moore (1982 and 1986). Results of tests on the generator using the multiplier 16807 with and without shuffling are reported by Learmonth and Lewis (1973b). If the user wishes to perform additional tests, the routines in Chapter 7, “Tests of Goodness of Fit and Randomness” may be of use. The user may also wish to compute some basic statistics or to make some plots of the output of the random number generator being used. The routines in Chapter 1, “Basic Statistics” and Chapter 16, “Line Printer Graphics” may be used for this purpose. Often in Monte Carlo applications, it is appropriate to construct an ad hoc test that is sensitive to departures that are important in the given application. For example, in using Monte Carlo methods to evaluate a one‑dimensional integral, autocorrelations of order one may not be harmful, but they may be disastrous in evaluating a two‑dimensional integral. Although generally the routines in this chapter for generating random deviates from nonuniform distributions use exact methods, and, hence, their quality depends almost solely on the quality of the underlying uniform generator, it is often advisable to employ an ad hoc test of goodness of fit for the transformations that are to be applied to the deviates from the nonuniform generator.

With release 7.0, three new subroutines associated with copulas have been added to the Fortran Numerical Library . A copula is a multivariate cumulative probability distribution (CDF) whose arguments are random variables uniformly distributed on the interval [0, 1] corresponding to the probabilities (variates) associated with arbitrarily distributed marginal deviates. The copula structure allows the multivariate CDF to be partitioned into the copula, which has associated with it information characterizing the dependence among the marginal variables, and the set of separate marginal deviates, each of which has its own distribution structure.

Two subroutines, RNMVGC and RNMVTC, allow the user to specify a correlation structure (in the form of a Cholesky matrix) which can be used to imprint correlation information on a sequence of multivariate random vectors. Each call to one of these methods returns a random vector whose elements (variates) are each uniformly distributed on the interval [0, 1] and correlated according to a user‑specified Cholesky matrix. These variate vector sequences may then be inverted to marginal deviate sequences whose distributions and imprinted correlations are user‑specified.

Method RNMVGC generates a random Gaussian copula vector by inverting a vector of uniform [0, 1] random numbers to a N(0, 1) deviate vector, imprinting the N(0,1) vector with the correlation information by multiplying it with the Cholesky matrix, and then using the N(0,1) CDF to map the Cholesky‑imprinted deviate vector back to a vector of imprinted uniform [0, 1] variates.

The third copula subroutine, CANCOR, extracts a “canonical correlation” matrix from a sequence of multivariate deviate vectors whose component marginals are arbitrarily distributed. This is accomplished by first extracting the empirical CDF from each of the marginal deviate sequences and then using this empirical CDF to map the deviates to uniform [0, 1] variates which are then inverted to N(0, 1) deviates. Each element Ci j of the canonical correlation matrix can then be extracted by averaging the products zi t zj t of N(0, 1) deviates i and j over the t‑indexed sequence. The utility of subroutine CANCOR is that because the canonical correlation matrix is derived from N(0, 1) deviates, the correlation is unbiased, i.e. undistorted by the arbitrary marginal distribution structures of the original deviate vector sequences. This is important in such financial applications as portfolio optimization, where correlation is used to estimate and minimize risk.

The use of subroutines RNMVGC, RNMVTC, and CANCOR is illustrated in the examples following subroutines RNMVGC and RNMVTC. The example following RNMVGC first uses method RNMVGC to create a correlation imprinted sequence of random deviate vectors and then uses method CANCOR to extract the correlation matrix from the imprinted sequence of vectors. Similarly, The example following RNMVTC first uses method RNMVTC to create a correlation imprinted sequence of random deviate vectors and then uses method CANCOR to extract the correlation matrix from the imprinted sequence of vectors.

The generators for continuous distributions are available in both single and double precision versions. This is merely for the convenience of the user; the double precision versions should not be considered more “accurate,” except possibly for the multivariate distributions.

The names of all of the routines for random number generation begin with “RN” for single precision and “DRN” for double precision. In most routines, the first argument, NR, is the number of variates to generate; and the last variable, either R or IR, is the vector of random variates.

Error handling and workspace allocation in the routines for random number generation are done somewhat differently than in most other IMSL routines. In general, there is less error checking than in other routines since there is more emphasis on speed in the random number generation routines. Simple checks for gross errors are made in all routines; and the routines for setup do complete checking since it is assumed that they would not be called frequently. Some routines, such as those that construct tables or interpolate from tables, require that the user explicitly provide some work arrays.