Within the context of this document, a base random number generator (BRNG) is a
+mathematical algorithm that, given an initial state, produces a sequence (or stream)
+of variates (or values) uniformly distributed over the semi-open interval (0,1].
+Note that this definition means that the value 1.0 may be returned, but the value
+0.0 will not. The period of the BRNG is defined as the maximum number of values
+that can be generated before the sequence starts to repeat. The initial state of
+a BRNG is often called the seed.
+A pseudorandom number generator (PRNG) is a BRNG that produces a stream of variates
+that are independent and statistically indistinguishable from a random sequence.
+A PRNG has several advantages over a true random number generator in that the
+generated sequence is repeatable, has known mathematical properties and is usually
+much quicker to generate.
+A quasi-random number generator (QRNG) is like a PRNG but the variates generated
+are not statistically independent, being designed to give a more even distribution
+in multidimensional space. Many books on statistics and computer science have good
+introductions to PRNGs and QRNGs, see for example Knuth [1] or Banks [2].
+All the BRNGs supplied in the AOCL-RNG library are PRNGs.
+In addition to standard PRNGs, some applications require cryptographically secure generators.
+A PRNG is said to be cryptographically secure if there is no polynomial-time algorithm which,
+on input of the first l bits of the output sequence, can predict the (l + 1)st bit of the
+sequence with probability significantly greater than 0.5. This is equivalent to saying there
+exists no polynomial-time algorithm that can correctly distinguish between an output sequence
+from the PRNG and a truly random sequence of the same length with probability significantly
+greater than 0.5 [3].
+A distribution generator is a routine that takes variates generated from a BRNG and
+transforms them into variates from a specified distribution, for example the Gaussian
+(Normal) distribution.
+The AOCL-RNG library contains six base generators, and twenty-three distribution generators.
+In addition, users can supply a custom-built generator as the base generator for all the
+distribution generators.
+The base generators were tested using the Big Crush, Small Crush and Pseudo Diehard test
+suites from the TestU01 software library <8th doc>.
+
+Base Generators
+The six base generators (BRNG) supplied with the AOCL-RNG library are;
+the NAG Basic Generator, a series of Wichmann-Hill Generator,
+the Mersenne Twister Generator, L’Ecuyer’s Combined Recursive Generator (MRG32k3a),
+the Blum-Blum-Shub Generator and the SIMD oriented Fast Mersenne Twister Generator (SFMT) generator.
+Some of the generators have been slightly modified from their usual form
+to make them consistent between themselves. For instance, the Wichmann-Hill
+generators in standard form may return exactly 0.0 but not exactly 1.0.
+In this library, we return 1.0 \(\mathcal x\) to convert the value x
+into the semi-open interval (0, 1] without affecting any other randomness
+properties. The original Mersenne Twister algorithm returns an exact zero
+about one time in a few billion; the AOCL-RNG implementation returns a tiny
+non-zero number as surrogate for zero. Same goes for SFMT as well.
+If a single stream of variates is required it is recommended that the
+Mersenne Twister base generator is used. This generator combines speed
+with good statistical properties and an extremely long period. SFMT
+provides all features of Mersenne Twister along with better speed.
+The NAG basic generator is another quick generator suitable for generating
+a single stream. However it has a shorter period than the Mersenne Twister
+and being a linear congruential generator, its statistical properties are
+not as good.
+If 273 or fewer multiple streams, with a period of up to 280 are required
+then it is recommended that the Wichmann-Hill generators are used. For more
+streams or multiple streams with a longer period it is recommended that the
+L’Ecuyer combined recursive generator is used in combination with the skip
+ahead routine. Generating multiple streams of variates by skipping ahead is
+generally quicker than generating the streams using the leap frog method.
+The Blum-Blum-Shub generator should only be used if a cryptographically
+secure generator is required. This generator is extremely slow and has
+poor statistical properties when used as a base generator for any of the
+distributional generators.
+We recommend to use AOCL-SecureRNG library wherever
+cryptographically secure random numbers generator is required.
+Refer Secure RNG whitepaper
+for more information on AMD’s secure random number generator.
+
+Initialization of the Base Generators
+A random number generator must be initialized before use.
+Three routines are supplied within the library for this purpose:
+Initialization Functiontions, User’s own Base Generator Hooks, and Blum-Blum-Shub Generator Functions.
+Both double and single precision versions of all RNG routines are supplied.
+Double precision names are prefixed by DRAND, and single precision by SRAND.
+Note that if a generator has been initialized using the relevant double precision
+routine, then the double precision versions of the distribution generators must also
+be used, and vice versa. This even applies to generators with no double or single
+precision parameters; for example, a call of Discrete Univariate Distribution Functions must be preceded
+by a call to one of the double precision initializers (typically DRANDINITIALIZE).
+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.
+The statistical properties of a sequence of random numbers are only guaranteed
+within the sequence, and not between sequences provided by the same generator.
+Therefore it is likely that repeated initialization will render the numbers
+obtained less, rather than more, independent. In most cases there should only
+be a single call to one of the initialization routines, per application, and
+this call must be made before any variates are generated. One example of where
+multiple initialization may be required is briefly touched upon in
+Multiple Stream Generators.
+In order to initialize the Blum-Blum-Shub generator a number of additional
+parameters, as well as an initial state (seed), are required. Although this
+generator can be initialized through the DRANDINITIALIZE routine it is
+recommended to use Blum-Blum-Shub routines instead.
+DRANDINITIALIZE / SRANDINITIALIZE
+Initialize one of the six supplied base generators:
+NAG basic generator,
+Wichmann-Hill generator,
+Mersenne Twister Generator,
+L’Ecuyer’s Combined Recursive Generator (MRG32k3a),
+Blum-Blum-Shub Generator,
+SIMD-oriented Fast Mersenne Twister (SFMT) Generator.
+(Note that SRANDINITIALIZE is the single precision version of DRANDINITIALIZE.
+The argument lists of both routines are identical except that any double precision arguments
+of DRANDINITIALIZE are replaced in SRANDINITIALIZE by single precision arguments,
+i.e. type REAL in Fortran or type float in C).
+For API information refer Initialization Functiontions
+C Generate 100 values from the Beta distribution
+ INTEGER LSTATE,N
+ PARAMETER (LSTATE=16,N=100)
+ INTEGER I,INFO,SEED(1),STATE(LSTATE)
+ DOUBLE PRECISION A,B
+ DOUBLE PRECISION X(N)
+
+C Set the seed
+ SEED(1) = 1234
+
+C Read in the distributional parameters
+ READ(5,*) A,B
+
+C Initialize the STATE vector
+ CALL DRANDINITIALIZE(1,1,SEED,1,STATE,LSTATE,INFO)
+
+C Generate N variates from the Beta distribution CALL
+ DRANDBETA(N,A,B,STATE,X,INFO)
+
+C Print the results
+ WRITE(6,*) (X(I),I=1,N)
+
+
Blum-Blum-Shub Generator Functions initialization routine for the Blum-Blum-Shub generator.
+Unlike the other base generators supplied with the library, the Blum-Blum-Shub generator
+requires two additional parameters, p and q, as well as an initial state, s.
+The parameters p, q and s can be of an arbitrary size.
+In order to avoid overflow, these values are expressed as a polynomial in B,
+where B = 224. For example, p can be factored into a polynomial of order lp, with
+p = p1 + p2B + p3B2 + · · · + pl Blp−1. Similarly, q = q1 + q2B + q3B2 + · · · + ql Blq −1 and
+s = s1 + s2B + s3B2 + · · · + sl Bls−1.
+(Note that SRANDINITIALIZEBBS is the single precision version of DRANDINITIALIZEBBS.
+The argument lists of both routines are identical except that any double precision arguments
+of DRANDINITIALIZEBBS are replaced in SRANDINITIALIZEBBS by single precision arguments,
+i.e. type REAL in Fortran or type float in C).
+For API information refer Blum-Blum-Shub Generator Functions
+
+
+Calling the Base Generators
+With the exception of the Blum-Blum-Shub generator, there are no interfaces for
+direct access to the base generators. All of the base generators return variates
+uniformly distributed over the semi-open interval (0, 1]. This functionality can
+be accessed using the uniform distributional generator DRANDUNIFORM, with parameter
+A = 0.0 and parameter B = 1.0. The base generator used is, as usual, selected
+during the initialization process Initialization of the Base Generators.
+To directly access the Blum-Blum-Shub generator, use the routine Blum-Blum-Shub Generator Functions.
+DRANDBLUMBLUMSHUB / SRANDBLUMBLUMSHUB
+Allows direct access to the bit stream generated by the Blum-Blum-Shub generator.
+(Note that SRANDBLUMBLUMSHUB is the single precision version of DRANDBLUMBLUMSHUB.
+The argument lists of both routines are identical except that any double precision arguments
+of DRANDBLUMBLUMSHUB are replaced in SRANDBLUMBLUMSHUB by single precision arguments,
+i.e. type REAL in Fortran or type float in C).
+For API information refer Blum-Blum-Shub Generator Functions
+
+
+NAG Basic Generator
+The NAG basic generator is a linear congruential generator (LCG) and,
+like all LCGs, has the form:
+
+\[\begin{split}\begin{align}
+ x_i& = a_1 x_{i−1} \ mod \ m_1, \\
+ u_i& = \frac{ x_i}{m_1},
+\end{align}\end{split}\]
+where the \(u_i, i = 1, 2, · · ·\) form the required sequence.
+The NAG basic generator takes \(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 numerical 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
+LCGs, and, even though the period is relatively short compared
+to many of the newer generators, it is sufficiently large for
+many practical problems.
+
+
+Wichmann-Hill Generator
+The Wichmann-Hill base generator uses a combination of four
+linear congruential generators (LCGs) and has the form:
+
+\[\begin{split} \begin{align}
+ w_i& = a_1 w_{i−1}\ mod \ m_1 \\
+ x_i& = a_2 x_{i−1}\ mod\ m_2 \\
+ y_i& = a_3 y_{i−1}\ mod\ m_3z_i \\
+ & = a_4 z_{i−1}\ mod\ m_4 \\
+ u_i& = \mathbf[\frac{w_i}{m_1} \mathbf+ \frac{x_i}{m_2} \mathbf+
+ \frac{y_i}{m_3} \mathbf+ \frac{z_i}{m_4}]\ mod \ 1,
+\end{align}\end{split}\]
+where the \(u_i, i = 1, 2, · · ·\) form the required sequence.
+There are 273 sets of parameters, {\(a_i, m_i\) : i = 1, 2, 3, 4}, to choose from.
+These values have been selected so that the resulting generators are independent and
+have a period of approximately \(2^{80}\).
+
+
+Mersenne Twister Generator
+The Mersenne Twister is a twisted generalized feedback shift register generator.
+The algorithm is as follows:
+
+Set some arbitrary initial values \(x_1, x_2, · · · , x_r\), each consisting
+of w bits.
Letting
+
+\[\begin{split}A = \begin{pmatrix}
+ 0 & I_{w−1} \\
+ a_w & a_{w−1}...a_1 \\
+ \end{pmatrix}\end{split}\]
+where \(\mathbf I_{w−1}\) is the \((w-1)\times(w-1)\) 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
+
+\[\mathcal X_{i+r} = \mathbf(\mathcal X_{i+s} \oplus \mathbf(
+ \mathcal X_{i}^{(w:(l+1))} | \mathcal X_{i+1}^{l-1})A )\]
+where \(\mathcal X_{i}^{(w:(l+1))} | \mathcal X_{i+1}^{(l:1)}\) indicates
+the concatenation of the most significant (upper) \(w − l\) bits of
+\(\mathcal X_i\) and the least significant (lower) \(\mathcal l\)
+bits of \(\mathcal X_{i+1}\).
+
+
+\[\begin{split}\begin{align}
+ \mathcal Z& = \mathcal X_{i+r} ⊕ (\mathcal X_{i+r} \gg t1) \\
+ \mathcal Z& = \mathcal Z ⊕ ((\mathcal Z \ll t2)\ AND\ m_1) \\
+ \mathcal Z& = \mathcal Z ⊕ ((\mathcal Z \ll t3)\ AND\ m_2) \\
+ \mathcal Z& = \mathcal Z ⊕ (\mathcal Z \gg t_4) \\
+ \mathcal u_{i+r}& = \frac{\mathcal Z}{(2^w − 1)}, \\
+\end{align}\end{split}\]
+where \(t_1, t_2, t_3\) and \(t_4\) are integers and \(m_1\) and
+\(m_2\) are bit-masks and >>t and <<t represent a t bit
+shift right and left respectively, is bit-wise exclusively or (xor) operation
+and AND is a bit-wise and operation.
+The \(u_{i+r} : i = 1, 2,…\) then form a pseudo-random sequence,
+with \(u_i \in (0, 1)\), for all i.
+This implementation of the Mersenne Twister uses the following values
+for the algorithmic constants:
+
+\[\begin{split}\begin{align}
+ w& = 32 \\
+ a& = 0x9908b0df \\
+ l& = 31 \\
+ r& = 624 \\
+ s& = 397 \\
+ t_1& = 11 \\
+ t_2& = 7 \\
+ t_3& = 15 \\
+ t_4& = 18 \\
+ m_1& = 0x9d2c5680 \\
+ m_2& = 0xefc60000
+\end{align}\end{split}\]
+where the notation \(0 \times DD \dotsc`\) indicates the bit pattern of
+the integer whose hexadecimal representation is \(DD \dotsc\).
+This algorithm has a period length of approximately \(2^{19,937} − 1\) and
+has been shown to be uniformly distributed in 623 dimensions.
+
+
+SIMD oriented Fast Mersenne Twister Generator
+SIMD-oriented Fast Mersenne Twister (SFMT) is a new variant of Mersenne Twister.
+SFMT is a Linear Feedbacked Shift Register generator that generates 128-bit
+pseudorandom integer recursively.
+The algorithm is as follows:
+
+Set some arbitrary initial values \(W_0, W_1, · · · , W_{N-1}\), each
+consisting of 128 bits.
Perform recursive operation:
+
+\[g(W_0,…,W_{N-1}) = W_0A ⊕ W_MB ⊕W_{N-2}C ⊕W_{N-1}D\]
+Where \(W_0, W_M,…\) are 128-bit integers, and A,B,C,D are sparse
+\(128 \times 128\) matrices over (0,1) for which
+\(W_A,W_B,W_C,W_D\) can be computed.
+The degree of recursion N is [19937/128] = 156, and the linear
+transformations A,B,C,D are as follows.
+
+\[\begin{split}\begin{align}
+ wA& = (w^{128} \ll 8) ⊕ w; \verb| w is considered as a single 128-bit integer.| \\
+ wB& = (w^{32} \gg 11)\ \& \ (BFFFFFF6 BFFAFFFF DDFECB7F DFFFFFEF); \\
+ & \verb|w is considered as a quadruple of 32-bit integers for right-shift operation| \\
+ wC& = (w^{128} \gg 8); \verb| w is considered as a single 128-bit integer.| \\
+ wD& = (w^{32} \ll 18); \verb| w is considered as a quadruple of 32-bit integer.| \\
+\end{align}\end{split}\]
+This algorithm has a period length of approximately \(2^{19,937} − 1\) and
+has better equidistribution property than Mersenne Twister.
+
+
+L’Ecuyer’s Combined Recursive Generator
+The base generator referred to as L’Ecuyer’s combined recursive generator is
+referred to as MRG32k3a in and combines two multiple recursive generators:
+
+\[\begin{split}\begin{align}
+ x_i& = a_{11}x_{i−1} + a_{12}x_{i−2} + a_{13}x_{i−3}\ mod \ m_1 \\
+ y_i& = a_{21}y_{i−1} + a_{22}y_{i−2} + a_{23}y_{i−3}\ mod \ m_2 \\
+ z_i& = x_i − y_i \ mod \ m_1 \\
+ u_i& =\frac{z_i}{m_1} ,
+\end{align}\end{split}\]
+where the \(u_i, i = 1, 2, \dotsc\) form the required sequence and
+\(a_{11} = 0, a_{12} = 1403580,\)
+\(a_{13} =−810728, m_1 = 2^{32} − 209, a_{21} = 527612, a_{22} = 0,
+a_{23} = −1370589 \ and \ m_2 = 2^{32} − 22853\).
+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}\).
+
+
+Blum-Blum-Shub Generator
+The Blum-Blum-Shub pseudo random number generator is cryptographically secure
+under the assumption that the quadratic residuosity problem is intractable.
+The algorithm consists of the following:
+
+Generate two large and distinct primes, p and q, each congruent to 3 mod 4.
+Define m = pq.
Select a seed s taking a value between 1 and m-1, such that the greatest
+common divisor between s and m is 1. Let \(x_0 = s^2\ mod \ m\).
+For \(I = 1, 2,\dotsc`\) generate:
+
+\[\begin{split}\begin{align}
+ x_i& = x_2\ mod\ m \\
+ z_i& = v \verb| least significant bits of | x_i
+\end{align}\end{split}\]
+where \(v \ge 1\).
+
+
+
+User Supplied Generator
+All of the distributional generators, require a base generator which returns a
+uniformly distributed value in the semi-open interval (0, 1] and AOCL-RNG library
+includes several such generators. However, for greater flexibility, the library
+routines allow the user to register their own base generator function. This
+user-supplied generator then becomes the base generator for all of the
+distribution generators.
+A user supplied generator comes in the form of two routines,
+one to initialize the generator and one to generate a set of
+uniformly distributed values in the semi-open interval (0, 1].
+These two routines can be named anything, but are referred to as
+UINI for the initialization routine and UGEN for the generation
+routine in the following documentation.
+User’s own Base Generator Hooks can be used to register User-supplied generator.
+Once registered the generator can be accessed and used in the same manner as
+the library supplied base generators.
+Refer code snippet provided below to use
+User-supplied base generator.
+C Generate 100 values from the Uniform distribution using
+C a user supplied base generator
+INTEGER LSTATE,N
+PARAMETER (LSTATE=16,N=100)
+INTEGER I,INFO,NSKIP,SEED(1),STATE(LSTATE) INTEGER X(N)
+DOUBLE PRECISION A,B
+
+C Set the seed SEED(1) = 1234
+C Set the distributional parameters
+A = 0.0D0
+B = 1.0D0
+
+C Initialize the base generator. Here RNGNB0GND is a user
+C supplied generator and RNGNB0INI is its initializer
+CALL DRANDINITIALIZEUSER(RNGNB0INI,RNGNB0GND,1,0,SEED,
+ * LSEED,STATE,LSTATE,INFO)
+
+C Generate N variates from the Univariate distribution CALL
+DRANDUNIFORM(N,A,B,STATE,X,LDX,INFO)
+
+C Print the results
+WRITE(6,*) (X(I),I=1,N)
+
+
UINI (GENID,SUBID,SEED,LSEED,STATE,LSTATE,INFO)
+INTEGER GENID can be any the ID associated with the generator.
+INTEGER SUBID can be the sub-ID associated with the generator.
+INTEGER SEED(LSEED) is an array containing the initial seed for your generator.
+INTEGER LSEED is either the size of the SEED array, or a value < 1.
+On output if LSEED < 1 on entry, LSEED must be set to the required size of the SEED array.
+This allows a caller of UINI to query the required size.
+INTEGER STATE(LSTATE) if LSTATE < 1 on entry, STATE should be unchanged.
+Otherwise, STATE is a state vector holding internal details required by your generator.
+On exit from UINI, the array STATE must hold the following information:
+STATE(1) = ESTATE, where ESTATE is your minimum allowed size of array STATE.
+STATE(2) = MAGIC, where MAGIC is a magic number of your own choice.
+This can be used by your routine UGEN as a check that UINI has previously been called.
+STATE(3) = GENID STATE(4) = SUBID
+STATE(5) … STATE(ESTATE-1) = internal state values required by your generator routine UGEN;
+for example, the current value of your seed.
+STATE(ESTATE) = MAGIC, i.e. the same value as STATE(2).
+INTEGER LSTATE [Input/Output] is either the size of the STATE array, or a value < 1.
+On output if LSTATE < 1 on entry, LSTATE should be set to the required size of the STATE array,
+i.e. the value ESTATE as described above. This allows the caller of UINI to query the required size.
+Constraint: either LSTATE < 1 or LSTATE≥ ESTATE.
+INTEGER INFO [Output] returns an error code, to be used in whatever way you wish;
+for example, to flag an incorrect argument to UINI.If no error is encountered, UINI must set INFO to 0.
+UGEN (N,STATE,X,INFO).
+INTEGER N [Input] is the number of random numbers to be generated.
+INTEGER STATE(*) [Input/Output] is the internal state of your generator.
+DOUBLE PRECISION X(N) [Output] is the array of N uniform distributed random numbers,
+each in the semi-open interval (0.0, 1.0], i.e. 1.0 is a legitimate return value, but 0.0 is not.
+INTEGER INFO [Output] is a flag which you can use to signal an error in the call of UGEN;
+for example, if UGEN is called without being initialized by UINI.
+
+
+
+Multiple Stream Generators
+Multiple stream generation method
+It is often advantageous to be able to generate variates from multiple,
+independent, streams. For example when running a simulation in parallel
+on several processors. There are four ways of generating multiple streams
+using the routines available in the AOCL-RNG library:
+
+
+
(a) Using different seeds
+
(b) Using different sequences
+
(c) Block-Splitting or Skipping-Ahead
+
(d) Leap Frogging
+
+
+The four methods are detailed in the following sections. Of the four,
+(a) should be avoided in most cases, (b) is only really of any practical
+use when using the Wichmann-Hill generator, and is then still limited to
+273 streams. Both block-splitting and leap-frogging work using 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 block-splitting requires the user to know (approximately) the maximum
+number of variates required from each stream. Block-splitting requires no a-priori
+information on the number of streams required. In contrast leap-frogging requires
+the user to know the maximum number of streams required, prior to generating the
+first value.
+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 Initialization of the Base Generators.
+In addition, for more complicated generators like a L’Ecuyer’s multiple recursive
+generator leap-frogging can increase the time required to generate each variate
+compared to block-splitting. The additional time required by block-splitting occurs
+at the initialization stage, and not at the variate generation stage.
+Therefore in most instances block-splitting would be the preferred method for
+generating multiple sequences.
+
+Using Different Seeds
+A different sequence of variates can be generated from the same base generator
+by initializing the generator using a different set of seeds. Of the four methods
+for creating multiple streams described here, this is the least satisfactory.
+As mentioned in Initialization of the Base Generators, the statistical properties
+of the base generators are only guaranteed within sequences, not between sequences.
+For example, sequences generated from different starting points may overlap if the
+initial values are not far enough apart. The potential for overlapping sequences is
+reduced if the period of the generator being used is large. Although there is no
+guarantee of the independence of the sequences, due to its extremely large period,
+using the Mersenne Twister with random starting values is unlikely to lead to problems,
+especially if the number of sequences required is small. This is the only way in which
+multiple sequences can be generated with the AOCL-RNG library using the Mersenne Twister
+as the base generator.
+If the statistical properties of different sequences must be provable then one of
+the other methods should be adopted.
+
+
+Using Different Generators
+Independent sequences of variates can be generated using different base generators
+for each sequence. For example, sequence 1 can be generated using the NAG basic generator,
+sequence 2 using the L’Ecuyer’s Combined Recursive generator, sequence 3 using the
+Mersenne Twister. The Wichmann-Hill generator implemented in the library is in fact
+a series of 273 independent generators. The particular sub-generator being used can
+be selected using the SUBID variable Initialization Functiontions. Therefore, in total,
+277 independent streams can be generated with each using an independent generator
+(273 Wichmann-Hill generators, and 4 additional base generators).
+
+
+Skip Ahead
+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 larger than
+the maximum number of variates required from any of the sequences.For example:
+
+\[\frac{x_1,x_2, \dotsi,x_n}{block 1}
+\frac{x_{n+1},x_{n+2}, \dotsi,x_{2n}}{block 2}
+\frac{x_{2n+1},x_{2n+2}, \dotsi,x_{3n}}{block 3},\ etc\]
+where \(x_1, x_2\), is the sequence produced by the generator of interest.
+Each of the k blocks provide an independent sequence.
+The block splitting algorithm therefore requires the sequence to be advanced a
+large number of places. Due to their form this can be done efficiently for
+linear congruential generators and multiple congruential generators.
+The AOCL-RNG library provides block-splitting for the NAG Basic generator,
+the Wichmann-Hill generators and L’Ecuyer’s Combined Recursive generator.
+C Generate 3 * 100 values from the Uniform distribution
+C Multiple streams generated using the Skip Ahead method
+ INTEGER LSTATE,N
+ PARAMETER (LSTATE=16,N=100)
+ INTEGER I,INFO,NSKIP
+ INTEGER SEED(1),STATE1(LSTATE),STATE2(LSTATE),STATE3(LSTATE)
+ DOUBLE PRECISION X1(N),X2(N),X3(N)
+ DOUBLE PRECISION A,B
+
+C Set the seed
+ SEED(1) = 1234
+
+C Set the distributional parameters
+ A = 0.0D0
+ B = 1.0D0
+
+C Initialize the STATE1 vector
+ CALL DRANDINITIALIZE(1,1,SEED,1,STATE1,LSTATE,INFO)
+
+C Copy the STATE1 vector into other state vectors
+ DO 20 I = 1,LSTATE
+ STATE2(I) = STATE1(I)
+ STATE3(I) = STATE1(I)
+20 CONTINUE
+
+C Calculate how many places we want to skip, this
+C should be >> than the number of variates we
+C wish to generate from each stream
+ NSKIP = N * N
+
+C Advance each stream, first does not need changing
+ CALL DRANDSKIPAHEAD(NSKIP,STATE2,INFO)
+ CALL DRANDSKIPAHEAD(2*NSKIP,STATE3,INFO)
+
+C Generate 3 sets of N variates from the Univariate distribution
+ CALL DRANDUNIFORM(N,A,B,STATE1,X1,LDX,INFO)
+ CALL DRANDUNIFORM(N,A,B,STATE2,X2,LDX,INFO)
+ CALL DRANDUNIFORM(N,A,B,STATE3,X3,LDX,INFO)
+
+C Print the results
+ DO 40 I = 1,N
+ WRITE(6,*) X1(I),X2(I),X3(I)
+40 CONTINUE
+
+
Leap Frog
+Independent sequences of variates can 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:
+
+\[\begin{split}\begin{align}
+ Subsequence 1& : x_1, x_{k+1}, x_{2k+1}, \dotsc \\
+ Subsequence 2& : x_2, x_{k+2}, x_{2k+2}, \dotsc \\
+ \ddots \\
+ Subsequence k& : x_k, x2k, x3k,\dotsc
+\end{align}\end{split}\]
+each subsequence is then provides an independent stream.
+The leap-frog algorithm therefore requires the generation of every kth
+variate of a sequence. Due to their form this can be done efficiently for
+linear congruential generators and multiple congruential generators.
+The library provides leap-frogging for the NAG Basic generator,
+the Wichmann-Hill generators and L’Ecuyer’s Combined Recursive generator.
+As an illustrative example, a brief description of the algebra behind the implementation
+of the leap-frog algorithm (and block-splitting algorithm) for a linear congruential
+generator (LCG) will be given. A linear congruential generator has the form
+\(x_{i+1}\ =\ a_1 x_i\ mod\ m1\). The recursive nature of a LCG means that
+
+\[\begin{split}\begin{align}
+ x_{i+v}& = a_1 x_{i+v−1}\ mod\ m_1 \\
+ & = a_1(a_1 x_{i+v−2}\ mod\ m1)\ mod\ m_1 \\
+ & = a_2 x_{i+v−2}\ mod\ m_1 \\
+ & = a_1^v x_i\ mod\ m_1 \\
+\end{align}\end{split}\]
+The sequence can be quickly advanced v places by multiplying the current
+state \((x_i)\) by \(a_1^v\ mod\ m_1\), hence allowing block-splitting.
+Leap-frogging is implemented by using \(a_k\), where k is the number of
+streams required, in place of \(a_1\) in the standard LCG recursive formula.
+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_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 26mphasized by
+the lattice structure of LCGs.
+Note that, due to rounding, a sequence generated using leap-frogging and a sequence
+constructed by taking every kth value from a set of variates generated
+without leap-frogging may differ slightly. These differences should only affect
+the least significant digit.
+C Generate 3 * 100 values from the Uniform distribution
+C Multiple streams generated using the Leap Frog method
+ INTEGER LSTATE,N
+ PARAMETER (LSTATE=16,N=100)
+ INTEGER I,INFO
+ INTEGER SEED(1),STATE1(LSTATE),STATE2(LSTATE),STATE3(LSTATE)
+ DOUBLE PRECISION X1(N),X2(N),X3(N)
+ DOUBLE PRECISION A,B
+
+C Set the seed
+ SEED(1) = 1234
+
+C Set the distributional parameters
+ A = 0.0D0
+ B = 1.0D0
+
+C Initialize the STATE1 vector
+ CALL DRANDINITIALIZE(1,1,SEED,1,STATE1,LSTATE,INFO)
+
+C Copy the STATE1 vector into other state vectors
+ DO 20 I = 1,LSTATE
+ STATE2(I) = STATE1(I)
+ STATE3(I) = STATE1(I)
+20 CONTINUE
+
+C Update each stream so they generate every 3rd value
+ CALL DRANDLEAPFROG(3,1,STATE1,INFO)
+ CALL DRANDLEAPFROG(3,2,STATE2,INFO)
+ CALL DRANDLEAPFROG(3,3,STATE3,INFO)
+
+C Generate 3 sets of N variates from the Univariate distribution
+ CALL DRANDUNIFORM(N,A,B,STATE1,X1,LDX,INFO)
+ CALL DRANDUNIFORM(N,A,B,STATE2,X2,LDX,INFO)
+ CALL DRANDUNIFORM(N,A,B,STATE3,X3,LDX,INFO)
+
+C Print the results
+ DO 40 I = 1,N
+ WRITE(6,*) X1(I),X2(I),X3(I)
+40 CONTINUE
+
+
+