The NAG basic generator is a linear congruential generator (LCG) and, like all linear congruential generators, has the form: where the $u_{\mathit{i}}$, for $\mathit{i}=1,2,\dots$, form the required sequence.

The NAG basic generator uses $a_1={13}^{13}$ and $m_1=2^{59}$, which gives a period of approximately $2^{57}$.

This generator has been part of the NAG Library since Mark 6 and as such has been widely used. It suffers from no known problems, other than those due to the lattice structure inherent in all linear congruential generators, and, even though the period is relatively short compared to many of the newer generators, it is sufficiently large for many practical problems.

The performance of the NAG basic generator has been analysed by the Spectral Test, see Section 3.3.4 of Knuth (1981), yielding the following results in the notation of Knuth (1981).

The right-hand column gives an upper bound for the values of ${\nu }_n$ attainable by any multiplicative congruential generator working modulo $2^{59}$.

An informal interpretation of the quantities ${\nu }_n$ is that consecutive $n$-tuples are statistically uncorrelated to an accuracy of $1/{\nu }_n$. This is a theoretical result; in practice the degree of randomness is usually much greater than the above figures might support. More details are given in Knuth (1981), and in the references cited therein.

Note that the achievable accuracy drops rapidly as the number of dimensions increases. This is a property of all multiplicative congruential generators and is the reason why very long periods are needed even for samples of only a few random numbers.

This series of Wichmann–Hill base generators (see Maclaren (1989)) use a combination of four linear congruential generators and has the form: where the $u_{\mathit{i}}$, for $\mathit{i}=1,2,\dots$, form the required sequence. The NAG Library implementation includes 273 sets of parameters, $a_{\mathit{j}},m_{\mathit{j}}$, for $\mathit{j}=1,2,3,4$, to choose from.

The constants $a_i$ are in the range 112 to 127 and the constants $m_j$ are prime numbers in the range $16718909$ to $16776971$, which are close to $2^{24}=16777216$. These constants have been chosen so that each of the resulting 273 generators are essentially independent, all calculations can be carried out in 32-bit integer arithmetic and the generators give good results with the spectral test, see Knuth (1981) and Maclaren (1989). The period of each of these generators would be at least $2^{92}$ if it were not for common factors between $\left(m_1-1\right)$, $\left(m_2-1\right)$, $\left(m_3-1\right)$ and $\left(m_4-1\right)$. However, each generator should still have a period of at least $2^{80}$. Further discussion of the properties of these generators is given in Maclaren (1989).

This Wichmann–Hill base generator (see Wichmann and Hill (2006)) is of the same form as that described in [Wichmann–Hill I Generator], i.e., a combination of four linear congruential generators. In this case $a_1=11600$, $m_1=2147483579$, $a_2=47003$, $m_2=2147483543$, $a_3=23000$, $m_3=2147483423$, $a_4=33000$, $m_4=2147483123$.

Unlike in the original Wichmann–Hill generator, these values are too large to carry out the calculations detailed in (1) using 32-bit integer arithmetic, however, if then setting gives and $W_i$ can be calculated in 32-bit integer arithmetic. Similar expressions exist for $x_i$, $y_i$ and $z_i$. The period of this generator is approximately $2^{121}$.

Further details of implementing this algorithm and its properties are given in Wichmann and Hill (2006). This paper also gives some useful guidelines on testing PRNGs.

The Mersenne Twister (see Matsumoto and Nishimura (1998)) is a twisted generalized feedback shift register generator. The algorithm underlying the Mersenne Twister is as follows:

(i)	Set some arbitrary initial values $x_1,x_2,\dots ,x_r$, each consisting of $w$ bits.
(ii)	Letting $$A=\left(\begin{array}{cc}0& I_{w-1}\\ a_w& a_{w-1}\cdots a_1\end{array}\right)\text{,}$$ where $I_{w-1}$ is the $\left(w-1\right)\times \left(w-1\right)$ identity matrix and each of the $a_i,i=1$ to $w$ take a value of either $0$ or $1$ (i.e., they can be represented as bits). Define $$x_{i+r}=\left(x_{i+s}\oplus \left(x_i^{\left(\omega :\left(l+1\right)\right)}|x_{i+1}^{\left(\mathrm{l}:1\right)}\right)A\right)\text{,}$$ where $x_i^{\left(\omega :\left(l+1\right)\right)}|x_{i+1}^{\left(\mathrm{l}:1\right)}$ indicates the concatenation of the most significant (upper) $w-l$ bits of $x_i$ and the least significant (lower) $l$ bits of $x_{i+1}$.
(iii)	Perform the following operations sequentially: $$\begin{array}{ccc}z& =& x_{i+r}\oplus \left(x_{i+r}\gg t_1\right)\\ z& =& z\oplus \left(\left(z\ll t_2\right)\text{ AND }{m}_1\right)\\ z& =& z\oplus \left(\left(z\ll t_3\right)\text{ AND }{m}_2\right)\\ z& =& z\oplus \left(z\gg t_4\right)\\ u_{i+r}& =& z/\left(2^w-1\right)\text{,}\end{array}$$ where $t_1$, $t_2$, $t_3$ and $t_4$ are integers and $m_1$ and $m_2$ are bit-masks and ‘$\gg t$’ and ‘$\ll t$’ represent a $t$ bit shift right and left respectively, $\oplus$ is bit-wise exclusively or (xor) operation and ‘AND’ is a bit-wise and operation.

The $u_{\mathit{i}+r}$, for $\mathit{i}=1,2,\dots$, form the required sequence. The supplied implementation of the Mersenne Twister uses the following values for the algorithmic constants: where the notation 0xDD $\dots$ indicates the bit pattern of the integer whose hexadecimal representation is DD $\dots$.

This algorithm has a period length of approximately $2^{\mathrm{19,937}}-1$ and has been shown to be uniformly distributed in 623 dimensions (see Matsumoto and Nishimura (1998)).

The ACORN generator is a special case of a multiple recursive generator (see Wikramaratna (1989) and Wikramaratna (2007)). The algorithm underlying ACORN is as follows:

(i)	Choose an integer value $k\ge 1$.
(ii)	Choose an integer value $M$, and an integer seed $Y_0^{\left(0\right)}$, such that $0<Y_0^{\left(0\right)}<M$ and $Y_0^{\left(0\right)}$ and $M$ are relatively prime.
(iii)	Choose an arbitrary set of $k$ initial integer values, $Y_0^{\left(1\right)},Y_0^{\left(2\right)},\dots ,Y_0^{\left(k\right)}$, such that $0\le Y_0^{\left(m\right)}<M$, for all $m=1,2,\dots ,k$.
(iv)	Perform the following sequentially: $$\begin{array}{c}{Y}_i^{\left(m\right)}=\left(Y_i^{\left(m-1\right)}+Y_{i-1}^{\left(m\right)}\right)\mathrm{ mod }M\end{array}$$ for $m=1,2,\dots ,k$.
(v)	Set $u_i=Y_i^{\left(k\right)}/M$.

The $u_{\mathit{i}}$, for $\mathit{i}=1,2,\dots$, then form a pseudorandom sequence, with $u_i\in \left[0,1\right)$, for all $i$.

Although you can choose any value for $k$, $M$, $Y_0^{\left(0\right)}$ and the $Y_0^{\left(m\right)}$, within the constraints mentioned in (i) to (iii) above, it is recommended that $k\ge 10$, $M$ is chosen to be a large power of two with $M\ge 2^{60}$ and $Y_0^{\left(0\right)}$ is chosen to be odd.

The period of the ACORN generator, with the modulus $M$ equal to a power of two, and an odd value for $Y_0^{\left(0\right)}$ has been shown to be an integer multiple of $M$ (see Wikramaratna (1992)). Therefore, increasing $M$ will give a series with a longer period.

The base generator L'Ecuyer MRG32k3a (see L'Ecuyer and Simard (2002)) combines two multiple recursive generators: where $a_{11}=0$, $a_{12}=1403580$, $a_{13}=-810728$, $m_1=2^{32}-209$, $a_{21}=527612$, $a_{22}=0$, $a_{23}=-1370589$, $m_2=2^{32}-22853$, and $u_i,i=1,2,\dots$ form the required sequence. If $d=m_1$ then $u_i\in \left(0,1\right]$ else if $d=m_1+1$ then $u_i\in \left(0,1\right)$. Combining the two multiple recursive generators (MRG) results in sequences with better statistical properties in high dimensions and longer periods compared with those generated from a single MRG. The combined generator described above has a period length of approximately $2^{191}$.

Fix two times $t_0<T$ and let $W={\left(W_t\right)}_{0\le t\le T-t_0}$ be a standard $d$-dimensional Wiener process on the interval $\left[0,T-t_0\right]$. Recall that the terms Wiener process and Brownian motion are often used interchangeably.

A standard $d$-dimensional Brownian bridge $B={\left(B_t\right)}_{t_0\le t\le T}$ on $\left[t_0,T\right]$ is defined (see Revuz and Yor (1999)) as

The process is continuous, starts at zero at time $t_0$ and ends at zero at time $T$. It is Gaussian, has zero mean and has a covariance structure given by for any $s\le t$ in $\left[t_0,T\right]$ where $I_d$ is the $d$-dimensional identity matrix. The Brownian bridge is often called a non-free or ‘pinned’ Wiener process since it is forced to be $0$ at time $T$, but is otherwise very similar to a standard Wiener process.

We can generalize this construction as follows. Fix points $x,w\in {ℝ}^d$, let $\Sigma$ be a $d\times d$ covariance matrix and choose any $d\times d$ matrix $C$ such that $C{C}^{\mathrm{T}}=\Sigma$. The generalized $d$-dimensional Brownian bridge $X={\left(X_t\right)}_{t_0\le t\le T}$ is defined by setting for all $t\in \left[t_0,T\right]$. The process $X$ is continuous, starts at $x$ at time $t_0$ and ends at $w$ at time $T$. It has mean $\left(\left(t-t_0\right)w+\left(T-t\right)x\right)/\left(T-t_0\right)$ and covariance structure for all $s\le t$ in $\left[t_0,T\right]$. This is a non-free Wiener process since it is forced to be equal to $w$ at time $T$. However if we set $w=x+C{W}_{T-t_0}$, then $X$ simplifies to for all $t\in \left[t_0,T\right]$ which is nothing other than a $d$-dimensional Wiener process with covariance given by $\Sigma$.

Figure 1 shows two sample paths for a two-dimensional free Wiener process $X={\left(X_t^1,X_t^2\right)}_{0\le t\le 2}$. The correlation coefficient between the one-dimensional processes $X^1$ and $X^2$ at any time is $\rho =0.80$. Note that the red and green paths in each figure are uncorrelated, however it is fairly evident that the two red paths are correlated, and that the two green paths are correlated (when one path increases so does the other, and vice versa).

Figure 2 shows two sample paths for a two-dimensional non-free Wiener process. The process starts at $\left(0,0\right)$ and ends at $\left(1,-1\right)$. The correlation coefficient between the one-dimensional processes is again $\rho =0.80$. The red and green paths in each figure are uncorrelated, while the two red paths tend to increase and decrease together, as do the two green paths. Both Figure 1 and Figure 2 were constructed using (G05XBF not in this release).

The ideas above can also be used to construct sample paths of a free or non-free Wiener process (recall that a non-free Wiener process is the Brownian bridge process outlined above). Fix two times $t_0<T$ and let ${\left(t_i\right)}_{1\le i\le N}$ be any set of time points satisfying $t_0<t_1<t_2<\cdots <t_N<T$. Let ${\left(X_{t_i}\right)}_{1\le i\le N}$ denote a $d$-dimensional (free or non-free) Wiener sample path at these times. These values can be generated by the so-called Brownian bridge algorithm (see Glasserman (2004)) which works as follows. From any two known points $X_{t_i}$ at time $t_{\mathrm{i}}$ and $X_{t_k}$ at time $t_k$ with $t_i<t_k$, a new point $X_{t_j}$ can be interpolated at any time $t_j\in \left(t_i,t_k\right)$ by setting where $Z$ is a $d$-dimensional standard Normal random variable and $C$ is any $d\times d$ matrix such that $C{C}^{\mathrm{T}}$ is the desired covariance structure for the (free or non-free) Wiener process $X$. Clearly this algorithm is iterative in nature. All that is needed to complete the specification is to fix the start point $X_{t_0}$ and end point $X_T$, and to specify how successive interpolation times $t_j$ are chosen. For $X$ to behave like a usual (free) Wiener process we should set $X_{t_0}$ equal to some value $x\in {ℝ}^d$ and then set $X_T=x+C\sqrt{T-t_0}Z$ where $Z$ is any $d$-dimensional standard Normal random variable. However when it comes to deciding how the successive interpolation times $t_j$ should be chosen, there is virtually no restriction. Any method of choosing which $t_j\in \left(t_i,t_k\right)$ to interpolate next is equally valid, provided $t_i$ is the nearest known point to the left of $t_j$ and $t_k$ is the nearest known point to the right of $t_j$. In other words, the interpolation interval $\left(t_i,t_k\right)$ must not contain any other known points, otherwise the covariance structure of the process will be incorrect.

The order in which the successive interpolation times $t_j$ are chosen is called the bridge construction order. Since all construction orders will produce a correct process, the question arises whether one construction order should be preferred over another. When the $Z$ values are drawn from a pseudorandom generator, the answer is typically no. However the bridge algorithm is frequently used with quasi-random numbers, and in this case the bridge construction order can be important.

Consider the one-dimensional case of a free Wiener process where $d=C=1$. The Brownian bridge is frequently combined with low-discrepancy (quasi-random) sequences to perform quasi-Monte–Carlo integration. Quasi-random points $Z^1,Z^2,Z^3,\dots$ are generated from the standard Normal distribution, where each quasi-random point $Z^i=\left(Z_1^i,Z_2^i,\cdots ,Z_D^i\right)$ consists of $D$ one-dimensional values. The process $X$ starts at $X_{t_0}=x$ which is known. There remain $N+1$ time points at which the bridge is to be computed, namely ${\left(X_{t_i}\right)}_{1\le i\le N}$ and $X_T$ (recall we are considering a free Wiener process). In this case $D$ is set equal to $N+1$, so that $N+1$ dimensional quasi-random points are generated. A single quasi-random point is used to construct one Wiener sample path.

The question is how to use the dimension values of each $N+1$ dimensional quasi-random point. Often the ‘lower’ dimension values ($Z_1^i,Z_2^i$, etc.) display better uniformity properties than the ‘higher’ dimension values ($Z_{N+1}^i,Z_N^i$, etc.) so that the ‘lower’ dimension values should be used to construct the most important sections of the sample path. For example, consider a model which is particularly sensitive to the behaviour of the underlying process at time $3$. When constructing the sample paths, one would therefore ensure that time $3$ was one of the interpolation points of the bridge, and that a ‘lower’ dimension value was used in (2) to construct the corresponding bridge point $X_3$. Indeed, one would most likely also ensure that time $X_3$ was one of the first bridge points that was constructed: ‘lower’ dimension values would be used to construct both the left and right bridge points used in (2) to interpolate $X_3$, so that the distribution of $X_3$ benefits as much as possible from the uniformity properties of the quasi-random sequence. For further discussions in this regard we refer to Glasserman (2004). These remarks extend readily to the case of a non-free Wiener process.

The Brownian bridge algorithm, especially when combined with quasi-random variates, is frequently used to obtain numerical solutions to stochastic differential equations (SDEs) driven by (free or non-free) Wiener processes. The quasi-random variates produce a family of Wiener sample paths which cover the space of all Wiener sample paths fairly evenly. This is analogous to the way in which a two-dimensional quasi-random sequence covers the unit square ${\left[0,1\right]}^2$ evenly. When solving SDEs one is typically interested in the increments of the driving Wiener process between two time points, rather than the value of the process at a particular time point. [Brownian Bridge] contains details on which methods can be used to obtain such Wiener increments.

A different sequence of variates can be generated from the same base generator by initializing the generator using a different set of seeds. The statistical properties of the base generators are only guaranteed within, not between sequences. For example, two sequences generated from two different starting points may overlap if these initial values are not far enough apart. The potential for overlapping sequences is reduced if the period of the generator being used is large. In general, of the four methods for creating multiple streams described here, this is the least satisfactory.

The one exception to this is the Wichmann–Hill II generator. The Wichmann and Hill (2006) paper describes a method of generating blocks of variates, with lengths up to $2^{90}$, by fixing the first three seed values of the generator ($w_0$, $x_0$ and $y_0$), and setting $z_0$ to a different value for each stream required. This is similar to the skip-ahead method described in [Multiple Streams via Skip-ahead], in that the full sequence of the Wichmann–Hill II generator is split into a number of different blocks, in this case with a fixed length of $2^{90}$. But without the computationally intensive initialization usually required for the skip-ahead method.

Independent sequences of variates can be generated using a different base generator for each sequence. For example, sequence 1 can be generated using the NAG basic generator, sequence 2 using Mersenne Twister, sequence 3 the ACORN generator and sequence 4 using L'Ecuyer generator. The Wichmann–Hill I generator implemented in this chapter is, in fact, a series of 273 independent generators. The particular sub-generator to use is selected using the subid variable. Therefore, in total, 278 independent streams can be generated with each using a different generator (273 Wichmann–Hill I generators, and 5 additional base generators).

Independent sequences of variates can be generated from a single base generator through the use of block-splitting, or skipping-ahead. This method consists of splitting the sequence into $k$ non-overlapping blocks, each of length $n$, where $n$ is no smaller than the maximum number of variates required from any of the sequences. For example, where $x_1,x_2,\dots$ is the sequence produced by the generator of interest. Each of the $k$ blocks provide an independent sequence.

The skip-ahead algorithm therefore requires the sequence to be advanced a large number of places, as to generate values from say, block $b$, you must skip over the $\left(b-1\right)n$ values in the first $b-1$ blocks. Due to their form this can be done efficiently for linear congruential generators and multiple congruential generators. A skip-ahead algorithm is also provided for the Mersenne Twister generator.

Although skip-ahead requires some additional computation at the initialization stage (to ‘fast forward’ the sequence) no additional computation is required at the generation stage.

This method of producing multiple streams can also be used for the Sobol and Niederreiter quasi-random number generator via the parameter iskip in g05yl.

Independent sequences of variates can also be generated from a single base generator through the use of leap-frogging. This method involves splitting the sequence from a single generator into $k$ disjoint subsequences. For example: where $x_1,x_2,\dots$ is the sequence produced by the generator of interest. Each of the $k$ subsequences then provides an independent stream of variates.

The leap-frog algorithm therefore requires the generation of every $k$th variate from the base generator. Due to their form this can be done efficiently for linear congruential generators and multiple congruential generators. A leap-frog algorithm is provided for the NAG Basic generator, both the Wichmann–Hill I and Wichmann–Hill II generators and L'Ecuyer generator.

It is known that, dependent on the number of streams required, leap-frogging can lead to sequences with poor statistical properties, especially when applied to linear congruential generators. In addition, leap-frogging can increase the time required to generate each variate. Therefore leap-frogging should be avoided unless absolutely necessary.

As an illustrative example, a brief description of the algebra behind the implementation of the leap-frog and skip-ahead algorithms for a linear congruential generator is given. A linear congruential generator has the form $x_{i+1}=a_1{x}_i\mathrm{ mod }{m}_1$. The recursive nature of a linear congruential generator means that

The sequence can therefore be quickly advanced $v$ places by multiplying the current state ($x_i$) by $a_1^v\mathrm{ mod }{m}_1$, hence skipping the sequence ahead. Leap-frogging can be implemented by using $a_1^k$, where $k$ is the number of streams required, in place of $a_1$ in the standard linear congruential generator recursive formula, in order to advance $k$ places, rather than one, at each iteration.

In a linear congruential generator the multiplier $a_1$ is constructed so that the generator has good statistical properties in, for example, the spectral test. When using leap-frogging to construct multiple streams this multiplier is replaced with $a_1^k$, and there is no guarantee that this new multiplier will have suitable properties especially as the value of $k$ depends on the number of streams required and so is likely to change depending on the application. This problem can be emphasized by the lattice structure of linear congruential generators. Similiarly, the value of $a_1$ is often chosen such that the computation $a_1{x}_i\mathrm{ mod }{m}_1$ can be performed efficiently. When $a_1$ is replaced by $a_1^k$, this is often no longer the case.

Note that, due to rounding, when using a distributional generator, a sequence generated using leap-frogging and a sequence constructed by taking every $k$ value from a set of variates generated without leap-frogging may differ slightly. These differences should only affect the least significant digit.

Skipping ahead with the Mersenne Twister generator is based on the definition of a $k\times k$ (where $k=19937$) transition matrix, $A$, over the finite field ${𝔽}_2$ (with elements 0 and 1). Multiplying $A$ by the current state $x_n$, represented as a vector of bits, produces the next state vector $x_{n+1}$:

Thus, skipping ahead $v$ places in a sequence is equivalent to multiplying by $A^v$:

Since calculating $A^v$ by a standard square and multiply algorithm is $\mathit{O}\left(k^3\log \left(v\right)\right)$ and requires over 47MB of memory (see Haramoto et al. (2008)), an indirect calculation is performed which relies on a property of the characteristic polynomial $p\left(z\right)$ of $A$, namely that $p\left(A\right)=0$. We then define and observe that for a polynomial $q\left(z\right)$. Since $p\left(A\right)=0$, we have that $g\left(A\right)=A^v$ and

This polynomial evaluation can be performed using Horner's method: which reduces the problem to advancing the generator $k-1$ places from state $x_n$ and adding (where addition is as defined over ${𝔽}_2$) the intermediate states for which $a_i$ is nonzero.

There are therefore two stages to skipping the Mersenne Twister ahead $v$ places:

(i)	Calculate the coefficients of the polynomial $g\left(z\right)=z^v\mathrm{ mod }p\left(z\right)$;
(ii)	advance the sequence $k-1$ places from the starting state and add the intermediate states that correspond to nonzero coefficients in the polynomial calculated in the first step.

The resulting state is that for position $v$ in the sequence.

The cost of calculating the polynomial is $\mathit{O}\left(k^2\log \left(v\right)\right)$ and the cost of applying it to state is constant. Skip ahead functionality is typically used in order to generate $n$ independent pseudorandom number streams (e.g., for separate threads of computation). There are two options for generating the $n$ states:

(i)	On the master thread calculate the polynomial for a skip ahead distance of $v$ and apply this polynomial to state $n$ times, after each iteration $j$ saving the current state for later usage by thread $j$.
(ii)	Have each thread $j$ independently and in parallel with other threads calculate the polynomial for a distance of $\left(j+1\right)v$ and apply to the original state.

Since ${\displaystyle \underset{v\to \infty }{\lim }}\log \left(v\right)=\log nv$, then for large $v$ the cost of generating the polynomial for a skip ahead distance of $nv$ (i.e., the calculation performed by thread $n-1$ in option (ii) above) is approximately the same as generating that for a distance of $v$ (i.e., the calculation performed by thread $0$). However, only one application to state need be made per thread, and if $n$ is sufficiently large the cost of applying the polynomial to state becomes the dominant cost in option (i), in which case it is desirable to use option (ii). Tests have shown that as a guideline it becomes worthwhile to switch from option (i) to option (ii) for approximately $n>30$.

Leap frog calculations with the Mersenne Twister are performed by computing the sequence fully up to the required size and discarding the redundant numbers for a given stream.

Prior to generating any variates the base generator must be initialized. Two utility methods are provided for this, (G05KFF not in this release) and (G05KGF not in this release), both of which allow any of the base generators to be chosen.

(G05KFF not in this release) selects and initializes a base generator to a repeatable (when executed serially) state: two calls of (G05KFF not in this release) with the same argument-values will result in the same subsequent sequences of random numbers (when both generated serially).

(G05KGF not in this release) selects and initializes a base generator to a non-repeatable state in such a way that different calls of (G05KGF not in this release), either in the same run or different runs of the program, will almost certainly result in different subsequent sequences of random numbers.

No utilities for saving, retrieving or copying the current state of a generator have been provided. All of the information on the current state of a generator (or stream, if multiple streams are being used) is stored in the integer array state and as such this array can be treated as any other integer array, allowing for easy copying, restoring, etc.

As mentioned in [Multiple Streams via Different Initial Values (Seeds)], it is important to note that the statistical properties of pseudorandom numbers are only guaranteed within sequences and not between sequences produced by the same generator. Repeated initialization will thus render the numbers obtained less rather than more independent. In a simple case there should be only one call to (G05KFF not in this release) (G05KGF not in this release) and this call should be before any call to an actual generation method.

If a single sequence is required then it is recommended that the Mersenne Twister is used as the base generator ($\mathbf{genid}=3$). This generator is fast, has an extremely long period and has been shown to perform well on various test suites, see Matsumoto and Nishimura (1998), L'Ecuyer and Simard (2002) and Wichmann and Hill (2006) for example.

When choosing a base generator, the period of the chosen generator should be borne in mind. A good rule of thumb is never to use more numbers than the square root of the period in any one experiment as the statistical properties are impaired. For closely related reasons, breaking numbers down into their bit patterns and using individual bits may also cause trouble.

If the Wichmann–Hill II base generator is being used, and a period of $2^{90}$ is sufficient, then the method described in [Multiple Streams via Different Initial Values (Seeds)] can be used. If a different generator is used, or a longer period length is required then generating multiple streams by altering the initial values should be avoided.

Using a different generator works well if less than 277 streams are required.

Of the remaining two methods, both skip-ahead and leap-frogging use the sequence from a single generator, both guarantee that the different sequences will not overlap and both can be scaled to an arbitrary number of streams. Leap-frogging requires no a-priori knowledge about the number of variates being generated, whereas skip-ahead requires you to know (approximately) the maximum number of variates required from each stream. Skip-ahead requires no a-priori information on the number of streams required. In contrast leap-frogging requires you to know the maximum number of streams required, prior to generating the first value. Of these two, if possible, skip-ahead should be used in preference to leap-frogging. Both methods required additional computation compared with generating a single sequence, but for skip-ahead this computation occurs only at initialization. For leap-frogging additional computation is required both at initialization and during the generation of the variates. In addition, as mentioned in [Multiple Streams via Leap-frog], using leap-frogging can, in some instances, change the statistical properties of the sequences being generated.

Leap-frogging is performed by calling g05kh after the initialization method ( (G05KFF not in this release) (G05KGF not in this release)). For skip-ahead, either g05kj or g05kk can be called. Of these, g05kk restricts the amount being skipped to a power of $2$, but allows for a large ‘skip’ to be performed.

After calling one of the copula methods the inverse cumulative distribution function (CDF) can be applied to convert the uniform marginal distribution into the required form. Scalar and vector routines for evaluating the CDF, for a range of distributions, are supplied in G01 class. It should be noted that these routines are often described as computing the ‘deviates’ of the distribution.

When using the inverse CDF methods from G01 class it should be noted that some are limited in the number of significant figures they return. This may affect the statistical properties of the resulting sequence of variates. Section 7 of the individual method documentation will give a discussion of the accuracy of the particular algorithm being used and any available alternatives.