Implementation of a portable and reproducible parallel pseudorandom number generator

Implementation of a Portable and Reproducible Parallel Pseudorandom Number Generator Daniel V. Pryor, Steven A . Cuccaro, Michael Mascagni, M. L. Robinson Supercomputing Research Center, 1.D.A.

Abstract We descrabe an detaal the parallel amplementataon of a famaly of addatave lagged-Fabonacca pseudorandom number generators The theoretzcal structure of these generators as explozted t o preserve thew well-known randomness properties and to provade a parallel system of dastanct cycles. The algorzthm presented here solves the reproduczbalaty problem for a f a r larger class of parallel Monte Carlo applacataons than has been prevzously possable. I n partacular, Monte Carlo applacataons that undergo “splattang” can be coded i o be reproducable, zndependent both of the number of processors and the executaon order of the parallel processes. A labrary of portable C routanes (avaalable f r o m the authors) that amplements these adeas as also descrabed.

[7, 31. Given these simple conditions, the maximum possible period is obtained if and only if at least one of the initial e residues modulo 2” is odd. In other words, if the least significant bits of the seed are all zero, the maximum possible period is not obtained. For this generator there are (2‘ - 1)2‘(m-1) seeds that satisfy the condition for giving the maximum possible period. Since each of these seeds is in a cycle of maximum possible period, there must be

different maximal or full-period cycles. Each of these full-period cycles will be called an equivalence class (EC). In a previous paper, [8], the authors described a way to enumerate all E ECs based on the reduction of a seed from a given EC into a single unique seed representative of the entire EC. In addition, an explicit construction of all of the EC representatives was shown to have the form:

1. Introduction In Knuth’s well known exposition on pseudorandom number generation [5], several methods of generation are considered. Among these is the additive laggedFibonacci pseudorandom number generator:

m.s.b bE2

which we shall denote as F ( t , IC, m ) . In Marsaglia’s empirical study of pseudorandom number generators [6], this generator was among many considered. Overall, it did well on all of Marsaglia’s “stringent” tests, save for the “non-overlapping birthday spacing test.” However, Marsaglia noted that by choosing a generator with a large length, e, improvements are seen in this test. It is important to note that while the lagged-Fibonacci generators did have one shortcoming, in general they performed better than linear congruential generators on the “stringent” tests. This generator has a maximum possible period of (2‘ - 1)2“+l under certain easy to check conditions, Typeset by

(3)

...

bl 0 0

I I ‘S“b:

bot-1

bot-:!

2t-1 ~ e - 2

Here the b b i t long vector of least significant bits, bo, can be computed in advance for F ( l , k , m ) . We call this representation of F(!, 6 , m ) ,i.e, the bottom row of zeros and the precomputed column of least significant bits, its “canonical form.” This canonical form tableau leaves exactly (e - l ) ( m - 1) bits left to be specified ECs. to select one of the 2(‘-’)(’’+’) The simplest way to parallelize this generator is to associate the K t h parallel process with EC number IC.

d,+g-w

311

1063-9535/94$4.000 1994 IEEE

U

:::

The numbering is done lexicographically in the rectangle of boxes shown above. (Whether numbering is in row-column order or column-row order does not matter.) This “naive” approach has certain advantages. If we associate the EC indices with the nodes on a binary tree, any dynamic need for a new stream of pseudorandom numbers may be filled using a successor node on the tree. This is a local computation and can be implemented to be independent of the number of parallel processors and the order of parallel execution, [8]. This approach is quite prudent from the quality perspective, as the relationship between the ECs close on the binary tree minimizes the exponential sum crosscorrelation among the full-period sequences] [8]. However, starting all of the pseudorandom sequences from seeds that are in canonical form leads to so-called “flat spots.” This is because the seeding values used for small EC numbers are numerically small themselves. Thus the initial segments will start and remain numerically small for an unacceptably long stretch. Worse than that, since all of the ECs start from their canonical forms in the “naive” implementation, all low-numbered ECs will suffer from flat spots that are lined up with respect to their cycles. So not only is a low-numbered EC initially distorted, all of those with similar EC number will be similarly distorted. To illustrate the phenomenon of flat spots, Figure 1 shows the normalized initial sequences of two different ECs started from the EC canonical form.

om 0.60

...

A=

0 0 0 1

.. .. . .

... ... ... ... . ..

0 1 0 ... 0 0 0 ... 0 0 0 ... 0 0 0 ... .. .

.. .

.

..

.

.

.

.

.

.

.

.

0 0 0 0 0 0 0 0 0

. . . 0 0 0 0 0 0

. .. .. . . .. .. .. ... 0 0 0 ... 1 0 0 ... 0 0 0 ... 0 1 0

.. .. ..

.

I “‘P

I

0.00

a50

1.ca

1.50

2.00

0.00

am

I .00

150

2.00

0.00

am

1.00

1.33

200

0.20

0.00

I .a,

0.50 0.00 4.50

FIGURES l ~1~~ , AND IC. The initial sequence from two different ECs started at their EC canonical form. Upper plot (1A) is with EC number K = 0. Middle plot (1B) is with K = 1. Lower plot (1C) is the difference between 1A and 1B. The units for the abscissa is e x m and the random numbers are scaled to be between 0 and 1.

Thus by run up we mean that if we are in EC number Ir(n), , n ) eliminates the flat, spot problem and the initial cross correlation of neighboring sequences. Furthermore, this form of initialization can be individualized on a per-run basis by selecting a “global seed,” g, that maps each n to a different unique value n, in the above expression. Thus for any particular computation, g is set at the beginning and xo is initialized by first setting n = ( n @ g) 1 (where “@” denotes the bitwise exclusive OR operator) and . . . r2(ii), r(ii),n).’ then setting s, = ( I ’ I 5 ( i i ) , rI4(ii), Note that the rightmost entry of s is n rather than n. There is no reason to prefer one or the other in the present context, but we use n to conform to our usage in subsequent discussions, where it does make a difference. This definition of n requires that g E (0, 1, . . . , 231 - 3}, so that n (and hence all rk(ii) values) cannot be 0. This gives us an initialization for the Oth EC that is not all zeros, even if g is chosen to be 0. In addressing problem (2), note that the above seeding method could provide initialization for up to N = 231 - 2 unique full-period cycles. If the generators can be assigned statically, then very few applications would require this many independent generators. However, there are a number of applications which require that processes split, spawning one or more child

processes. These spawnings occur in tree structures that, while sparse, may be deep and very unbalanced, thus precluding the possibility of pre-allocating all possible child sequence numbers that may be needed. So we need a way for generators to spawn new generators with EC indices much higher than 231 - 2. Clearly such indices are available, there being 2496of them for F ( 17,5,32). We are now ready to describe the full implementation of this system of generators. Each EC is specified by the 16-long integer array e = (e15, e 1 4 , . . . , e l , eo), where, as noted above, e i E (0, 1 , 2 , .. ., 231 - 1). This representation of an EC number is meant to act like a 496-bit unsigned integer. In fact, we number the elements of e in descending order to remind us that it does represent a single number. Define n similarly. Now let G, denote a generator which is thought of as occupying node number n of the binary tree in Figure 2. 0

+

‘Other mappings of

n

to h are certainly possible.

314

4k

4k+l 4k+2

4k+3

FIGURE2. Enumeration of a binary tree.

G, is composed of two arrays: x, the current contents of the Fibonacci tableau (initialized with the canonical form seed as described above) and e,, the node number of the next child generator to be spawned by G,. Thus G, = (x;e,). What we are trying to accomplish is to construct, for each of the 2496nodes of the tree, a unique initial seed in canonical form. Additionally, we are providing each existing generator a means, based only on its own local information, of deciding on which node of the tree to spawn its next child generator.

Suppose a program requires N initial generators which will run independently of one another. For this algorithm the global seed g is chosen as any integer expressible in 30 bits, a slight modification of its previous range. (Remember, g is used only as an easy - single integer - way to introduce variation into a Monte Carlo program on a per-run basis. It’s default could very well be set to 0, without loss of “random looking” behavior.) Since we are initially setting up N generators, we will be placing them at nodes n = 0 , 1 , 2 , . . . , nmaz= N - 1 of the binary tree. a) set n = {(no & bitwise AND;

p30-

l))@g}+

1 where (‘&” denotes

. . , r(n)@nl,n o ) b) set s = (r15(ii)$n15, r14(h)@n14,. c) set e, = 2j(2n 1) for the smallest j such that 2j(2n 1) > n m a r .

+

+

+

Here the multiplication and addition 2n 1 is understood t o be done in the 496-bit representation of n, stored as explained above. We must now convince ourselves that this construction of seeds for ECs does indeed provide 2496distinct seeds. We are given n = (n15,. . .,no) and g. We 1. The constant form ii = {(no & (230 - 1)) @ g} g is used to change n from run to run. Thus without loss of generality we can set g = 0. Thus we have ii = (no & (230 - 1)) 1. If two ECs have different no’s, they have different s vectors. If two ECs have the same no’s, they will have the same ii’s. However, one of the raj’s, j E (1,. . . ,15} must be different if the two canonical forms are different, so their s vectors must be different. Keep in mind that the use of the LCG (the rk(ii) terms in step b) is to map, in a simple fashion, the initial seeds of neighboring tree nodes into seeds which are quite unlike each other. If we dispensed with this step by setting them equal to 0, then we would still achieve the goal of initializing distinct Fibonacci cycles. There are obviously an enormous number of permutation mappings on the set of 496-bit integers that could be used to “randomize” the initial seeds of neighboring generators. This particular LCG rule was c h e sen as our mapping because it is specifically known to have good “bitwise” randomness properties, particularly as related to the initial seed of Fibonacci generators, [l].Also, its implementation, as described in [9], is simple and extremely portable. We know of no 32bit environment that does not execute this algorithm as intended. Our only objection to this method of initialization is that it requires m - 1 to be a multiple of 31. For practical applications, we are comfortable with the constraint this seeding procedure entails.

+

+

Incidentally, an alternative strategy for accomplishing the same thing (without the use of a mapping) would be to advance each generator a fixed number of steps. Our computations indicate that, to provide good separation of neighboring generators, the number of steps required would be about e x m = 544 for the family F(17,5,32). For generators F ( l , 6,m) where k is close to e, the required number should be on the order of 2 x t x m. For example, if we go back t o Figure 1 we see that the initial flat spot is about half a unit of e x m in length.2 To produce a random number, a generator’s x array is updated according to equation (l), and the result is copied to a temporary location where it is right-shifted by one bit. This right shift eliminates the non-random least significant bit and scales the returned random number to a value in (0, 1 , 2 , .. ., 231 - 1). This sequence suffers no loss of period, as the returned values are the high 31 bits of a 32-bit generator. In other words, the period remains (217 - l)23*. Now suppose it is determined at a later time that the process accessing generator G, needs to spawn M new processes, each with its own generator. Each new generator must be different from the other new ones and different from any t h a t some other existing generator may spawn. The new tree nodes assigned in this operation are: e n , 2en, (Zen

+ I), de,,

(de,

+ I), (4en + 2), . . . ,n m a z

until M new tree nodes (EC indices) have been determined. In words, we are spawning generators on the nodes of the numerically smallest M-node subtree below, and including, node en. Having been assigned node numbers on the tree, the M new generators are now initialized by the same procedure used above to initialize the first N generators. Finally, the e, value for the original generator and each of the new generators is replaced by successive doublings, until the new values are greater than nmao;in this way any new child generators spawned will be different from all previously created generators. A small example will help to illustrate the spawning process just described. Suppose that 5 generators are needed at the start of a job. These would be placed at nodes 0,1,2,3, and 4 of the tree. Applying rule (c) we would have: e o = 8, el = 6,

e2

= 5, e3 = 7,

e4

=9

Now assume that, sometime later, generator 0 spawns 4 more children and that generator 3 spawns 6 more.

e

2The fact that the flat spot dissipates in about x m steps was the motivation for scaling the abscissa in Figure 1 by x m.

315

e

The new children spawned by generator 0 will have tree node numbers 8, 16, 17, and 32, and we would have

ea = 34, e16 = 33, e17 = 35,

e32

Table 1 Selected Canonical Forms

= 65

e o would be updated to 64. The new children spawned

by generator 3 will have tree node numbers 7, 14, 15, 28, 29, and 30, with e values e7

= 60,

e14

= 58,

e29

e15

= 31,

e28

= 57,

= 59, e g o = 61.

The value of e3 would be updated to 56. We have computationally determined canonical forms for a large number of F ( t , k, .) families. Indeed, we have these forms for Fibonacci generators corresponding to all primitive trinomials of degree less than 256, as well as several trinomials of greater degree. In Table 1 we list selected canonical configurations and note that for most (e, k ) pairs there are several least significant bit patterns that allow the ECs to be enumerated within the same rectangle of bits. Note also that L and S are independent of m [8]. Our practice is to denote as THE canonical form the numerically smallest such bit pattern which can be written with the least number of contiguous ones in a word of zeros. Thus for F ( 1 7 , 5 , 3 2 ) the canonical least significant bits are all zero, with a single one in position 10, there being no single ones in positions 0 through 9 which satisfy the definition of canonical form. In cases where no single one bit will work (only about 1 per cent of those computed) we have found canonical forms involving two consecutive one bits. For example, the canonical least significant bits for F ( 5 , 3 , .) are 00-1-1-0. Thus it appears, based on our computational results to date, that we can designate Fibonacci generators in canonical form completely as a quintuple of integers F ( ! , k , m , L , S ) , where L is the number of contiguous ones and S is the position number of the first one bit. We write the previous two examples as F ( 1 7 , 5 , 3 2 , 1 , 1 0 )and F ( 5 , 3 , ., 2, l ) , respectively.

3. The C implementation The library which implements the algorithms described above was written in the C language to take advantage of its intrinsic bit-shifting and memory management capabilities. (A set of wrapper routines is provided for FORTRAN users.) The library contains routines designed to deal with both the shared and distributed memory paradigms, and a pair of utility

e

k

3 5 10 17 35 55 71 93 127 158 521

2 3 7 5 2 24 65 91 97 128 168

m

L

S

.

1 2 1 1 1 1 1 2 1 1 1

0 1 7 10 0 11 1 1 21 63 83

routines to simplify movement of a generator from one processor to another using message passing. Let us assume that a particular program begins by initializing N generators, where N E { 1 , . . . , 231}. There are two methods for creating the initial generator or set of generators and assigning them to nodes. When this initial set is created using i n i t x n g s , a vector of generators assigned nodes 0 through N - 1 is returned. This routine is intended for shared memory systems; each element in the vector of pointers to the generators may be associated with a specific process from a set of tasks and the processes run in parallel. (Of course, more than one generator may be given to a particular process in cases where this would be useful.) For distributed memory systems, we have the routine i n i t m g d . To initialize N generators, the routine is called N times, each time provided with a different value of ni E (0, . . . , N } , the node index to initialize. Each call to the routine returns a pointer to the generator associated with that node index. The normal use of this routine is to have each virtual processor on a distributed system initialize its own generator or subset of generators. Once the initial set of generators has been created, random numbers may be obtained from the initialized nodes by calling the routine g e t m , which takes as argument a generator pointer. If in the course of the program a need for more generators develops, the routine spawning is used. The algorithm for this routine was described in 52. It is passed the pointer to the generator of the node which is spawning, and the next child of that node is spawned. If multiple nodes are being initialized, the next two nodes are obtained by spawning one from the original node and one from its child; the next four are obtained by spawning one from each

316

of the four nodes from the last step. In this way a full subtree whose root is the original node is generated, and the tree structure prevents repetition of generators. It also ensures that the random number generation will not be responsible for non-reproducibility of the results of the program, since the same generators may be assigned to the same processes independent of the order in which the processes are executed. For message passing systems, we provide a pair of routines, pacblng and unpacbmg. The routine pacblng creates an integer array containing all of the information needed to construct the generator in its current state. The array of integers may then be passed to another processor using whatever native message passing routines are provided to the user. The unpacking routine, unpackmg, can then reverse the packing process, creating a generator and a pointer to its location. (It is left to the user to free the array.) To prevent inadvertent duplication and to conserve memory, a third provided routine, free-mg, removes the generator and frees its pointer.

4. Discussion and conclusions We have presented the details of the parallel implementation of an additive lagged-Fibonacci pseudorandom number generator. This generator has good empirical randomness properties, and its parallel implementation is based on theoretical results, [8, 31, and empirical observations of its pseudorandom behavior within and between independent streams. We believe this family of generators solves, to a much greater extent than previously possible, the reproducibility problem in Monte Carlo applications, especially those applications in which process splitting occurs. We have produced a machine independent library of C subroutines that should run on any 32-bit architecture without modification. Considerable care has gone into this implementation to avoid obvious pitfalls. However, extra refinements of the current implementation are possible. Those things we have chosen not t o do fall into two categories: (1) those that are best done in an architecturally dependent way and (2) those generalizations that seemed unnecessary for the current F(1, k,32) implementation. Optimization of code for specific architectures is of this first type. While it is very easy to vectorize this generator for the production of pseudorandom numbers within a single EC, [2], this implementation makes no attempt to do so. It is hoped that the vendors of parallel processors will find our current library sufficiently useful that they will undertake machine specific

implementations where efficiency is built upon the current implementation's correctness. Refinements of the second type include allowing for variations in m, the number of bits in the pseudorandom integers. Our current implementation uses a modulo 231 - 1 LCG for reordering the ECs via I'(.). If one were to use an (m- 1)-bit maximal-period shift-register sequence to instantiate the function I?( .) in a general F(1, k,m) implementation, then all the properties of the current implementation would be retained.3 Namely, the (m - 1)-bit register will cycle over all 2m-1 - 1 non-zero register contents. Given an implementation that allows for arbitrary m, within some bounds, a native single and double precision floating-point implementation would then be possible. While we do provide floating-point pseudorandom numbers in the current implementation, they are computed by first producing a pseudorandom integer and then dividing by the modulus. More efficient is to use the floating-point recursion wt = W t - k wt-t (mod 1). In this way floating-point numbers can be directly produced without division. In addition, optimization that utilizes floating-point hardware can also be undertaken. However, a certain care must be taken to insure that the numerical roundoff does not cause a sequence to wander through many different ECs. Elsewhere, [3,8], it was shown that with an f-bit mantissa, generators may be correctly implemented in floatingpoint by setting m = f - 1 and initializing the least significant bit in the mantissa to zero.

+

3The required ( m - 1)-bit maximal-periodshift-register sequence necessitatesfindinga primitivepolynomialmodulo 2 and of degree m - 1. While primitive trinomialsmodulo 2 of all degrees do not exist, primitivepolynomialsmodulo2 of all degrees do exist, [4, 51.

317

Appendix: Subroutines and usage Initialization: i n t * * i n i t m g s ( i n t ngen, i n t length, i n t seed) i n t * i n i t m g - d ( i n t gennum, i n t length, i n t seed, int totalgen)

These routines initialize the system of generators and return either a single pointer or a vector of pointers to the generators initialized. In normal use, either the routine i n i t m g - s O would be called once, or the routine i n i t m g - d o would be called t o t a l g e n times with a different value of gennum each time. Although i n i t m g - d o is intended for initializing on multiple processors in a distributed-memory system and i n i t m g - s 0 for shared-memory systems, the user may implement whichever is more convenient. For example, i n i t m g - d o is normally called from each processor on a distributed system; however, i n i t m g - s () may be called from one processor to initialize a set of generators, after which the generators can be passed to the other nodes using the message passing routines provided. Both routines return the default case (one generator of length 17) if all parameters have the value 0. If there is insufficient memory to allocate the generators, the program exits.

i n t **spawnmg(int *genptr , i n t nspawned) This routine returns an nspawned-long vector of

pointers to generators. The particular generators that are initialized are determined using a branching algorithm which guarantees that generators will not be used twice. If there is insufficient memory to allocate the generators, the program exits. Inputs: - a pointer to a generator that can be accessed by the calling process (else indeterminate results) nspawned- the number of generators initialized and the number of pointers to be returned. Range: positive integers (else NULL pointer returned)

genptr

Return of random number: i n t g e t m ( i n t *genptr)

This routine returns the 31 most significant bits of the next number in the sequence produced by a generator. Inputs: - a pointer to a generator that can be accessed by the calling process (else indeterminate results)

genptr

Deletion of a generator:

Inputs: One of a set of allowed values for the l variable. (A k value that gives maximal period is selected.) If l is zero and the given value of length is not allowed, l is set to the default value. If f2 is non-zero and not equal to length, a zero input defaults to the current l value and any other input causes a fatal error. Range: {2,3,5,7,10,11,15,17,29,31,35, 55,63,127} (default = 17) seed - Permutes the positions of the generators on the binary enumeration tree, allowing the user to vary the generators used in successive runs. If e is non-zero and g is not equal to seed, a 0 input defaults to the current gseed value and any other input causes a fatal error. Note: g is described in §2. Range: any integer (converted to value modulo 230 - 1) ngen ( i n i t m g - s ) and t o t a l g e n ( i n i t s n g - d ) - The number of generators to be initialized. Range: positive integers (default = 1) length

-

void f r e e m g ( i n t *genptr)

This routine, provided for memory management purposes, deallocates the memory for a generator which will no longer be used. Inputs: genptr - a pointer to a generator that can be accessed by the calling process (else indeterminate results, possibly a fatal error)

Message Passing utilities: i n t * p a c k m g ( i n t *genptr)

This routine deallocates the generator pointed to by 3)-long vector of integers which contains the information needed by the u n p a c k i n g 0 routine to recreate the generator. The first three elements of this vector are the values of e, IC and g s e e d , respectively. If genptr does not point to a generator, a fatal error may result from deallocation. genptr - a pointer to a generator that can be accessed by the calling process (else indeterminate results, possibly a fatal error) genptr and allocates and returns a (2!

+

gennum ( i n i t i n g - d ) - The index of the generator

i n t *unpackmg(int *packed)

returned by this call. Range: [O,total_gen) (fatal error otherwise)

This routine initializes a generator that is described by the input vector packed, previously created by the

318

packlngo routine, and returns a pointer to the new

Inputs:

generator. If the process has a non-zero t value before the call to unpackxng0 and any of the t, k and gseed values passed do not agree with their current values, a fatal error occurs. packed - a pointer to a vector of integers, created by the p a c k m g o routine and presumably passed to the current process from another process. Other Utilities (possibly useful in debugging):

ngen

- the total number of generators initialized and the number of pointers to generators returned. Range: positive integers l l a g - the value of t for the generators. If the current value of t is non-zero and different from Uag, a fatal error occurs. Range: integers > 1 klag - the value of k for the generators. If the current value of t is non-zero and the current value of k is different from klag, a fatal error occurs. Range: [l,Ilag) f i l l - pointer to an array of ngen*llag integers. The first l l a g integers will become the initial values for the first.generator, and so on. Range: any integer

int g e t l l a g m g ( v o i d ) int getJtlagmg(void) int getseedxng(void)

These routines return the values of e, k and gseed for the calling process. If a generator has not yet been initialized in this process, each of these routines will return the value 0. int int int int

g e t h p t r m g ( i n t *genptr) * g e t f i l l i n g ( i n t *genptr) *getmode-indexmg (int *genptr) *getnext-indexmg (int *genptr)

References N.S. Altman, Bit-wise behavior of random number gen-

These routines allow the user t o examine a particular generator accessible t o the calling process. The routine g e t f i l 1 m g O returns an !-long vector of integers (the current contents of the register holding past values of the sequence), while g e t h p t r m g O returns an integer value which indexes the next element of this vector to be updated. getnode-indexmgo and getnextAndexmg0 return (e - l)-long vectors of integers, which specify the positions on the branching tree of the current generator or the next generator it will spawn, respectively.

erators, SIAM J. Sci. Stat. Comput. 9 (1988), 941-949. R. P. Brent, Uniform Random Number Generators for Supercomputers, Proceedings Fifth Australian Supercomputer Conference, SASC Organizing Committee, 1992, pp. 95-104. R. P. Brent, On the periods of generalized Fibonacci recurrences, in the Press, Math. Comput. (1994). S. W. Golomb, Shift Register Sequences, Revised Edition, Aegean Park Press, Laguna Hills, California, 1982. D. E. Knuth, The Art of Computer Progmmming, Vol. 2: Seminumerical Algorithms, Second edition, Addison-Wesley, Reading, Massachusetts, 1981. G . Marsagha, A current view of random number generators, Computing Science and Statistics: Proceedings of the XVIth Symposium on the Interface, 1985, pp. 310. G. Marsagha and L.-H. Tsay, Matrices and the structure of mndom number sequences, Linear Alg. and A p plic. 67 (1985), 147-156. M. Mascagni, S. A. Cuccaro, D. V. Pryor and M. L. Robinson, A fast, high-quality, and reproducible laggedFibonacci pseudorandom number generator, submitted, SRC-TR-94-115 (1994). S. K. Park and K . W. Miller, Random number generators: good ones are hard to find, Comm of the ACM 31 (1988), 1192-1201.

Inputs: genptr - a pointer to a generator that can be ac-

cessed by the calling process (else indeterminate results)

For experts only: int * * i n i t s p e c i f y m g ( i n t ngen, int l l a g , int klag, int * f i l l )

This routine allows the user to initialize laggedFibonacci generators with complete generality and to use the g e t m 0 , packlngo and unpackmg0 routines in the library t o handle these generators. The value returned is a pointer to a vector of generators (or NULL if there is insufficient memory). The routine spamumgO cannot be used t o initialize new generators when i n i t s p e c i f y m g 0 has been used, and therefore no guarantees as to the independence of these generators can be made. Any input value not in the allowed range results in a fatal error.

Supercomputing Research Center; InstitUte for Defense Analyses; 17100 Science Drive; Bowie, Maryland 20715-4300 USA E-mail address: pryofisuper.org, cuccaroQsuper.org, rascagniQsuper.org,robinsonQsuper.org

319