Section 3.

2: Multi–Stream Lehmer RNGs

Discrete-Event Simulation: A First Course

c
�2006 Pearson Ed., Inc. 0-13-142917-5

Discrete-Event Simulation: A First Course Section 3.2: Multi–Stream Lehmer RNGs 1/1
Section 3.2: Multi-Stream Lehmer RNGs

Typical DES models have many stochastic components

Want a unique source of randomness for each component
One (poor) option: multiple RNGs
Better option: one RNG with multiple “streams” of random
numbers
one stream per stochastic component
We will partition output from our Lehmer RNG into multiple
streams

Discrete-Event Simulation: A First Course Section 3.2: Multi–Stream Lehmer RNGs 2/1
Example 3.2.1: ssq2 Arrival and Service Processes

ssq2 has two stochastic components: arrival and service

Allocate a different generator state variable to each

GetService with Unique Seed

double GetService(void)
{
double s;
static long x = 12345;
PutSeed(x);
s = Uniform(1.0, 2.0);
GetSeed(&x);
return (s);
}

x represents the current state of the service process

Discrete-Event Simulation: A First Course Section 3.2: Multi–Stream Lehmer RNGs 3/1
Example 3.2.2: ssq2 Arrival and Service Processes

Arrival should have its own static variable, initialized

differently

GetArrival with Unique Seed

double GetArrival(void)
{
static double arrival = START;
static long x = 54321;
PutSeed(x);
arrival += Exponential(2.0);
GetSeed(&x);
return (arrival);
}

x represents the current state of arrival process

Discrete-Event Simulation: A First Course Section 3.2: Multi–Stream Lehmer RNGs 4/1
The Modified Arrival and Service Processes

As modified, arrival and service times are drawn from different

streams of random numbers
Nothing magic about the choice of seed for each stream
The choices may, in fact, be poor ones!
Provided the streams don’t overlap, the processes are
uncoupled
Execution time cost is negligible (see Example 3.2.3)

Discrete-Event Simulation: A First Course Section 3.2: Multi–Stream Lehmer RNGs 5/1
Stream Considerations

Potential problem: assignment of initial seeds to facilitate

streams
Each initial state should be chosen to produce disjoint streams
If states are picked at whim, no guarantee of disjoint streams
Some initial states may only be a few calls to Random apart
. ......
....
..
..
..
..
..
•.
12345
.. ..

It must
.. .
. .
. .
.
. .
. .

be some 54321
. .
. .
. •.
. .
. .

criterion
. .
. .
. .
. .
. .
. .

(a, m) = (48271, 231 − 1)

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

Discrete-Event Simulation: A First Course Section 3.2: Multi–Stream Lehmer RNGs 6/1
Jump Multipliers

We will develop a multi-stream version of rng

Theorem (3.2.1)
Given g (x) = ax mod m and integer j(1 < j < m − 1)
Jump function : g j (x) = (aj mod m)x mod m
Jump multiplier : aj mod m
If g (·) generates x0 , x1 , x2 , . . . then g j (·) generates x0 , xj , x2j , . . .

Theorem 3.2.1 is the key to creating streams

Discrete-Event Simulation: A First Course Section 3.2: Multi–Stream Lehmer RNGs 7/1
Example 3.2.4: An Example Jump Function

If m = 31 and a = 3 and j = 6, the jump multiplier is

aj mod m = 36 mod 31 = 16

If x0 = 1 then g (x) = 3x mod 31 generates

1, 3, 9, 27, 19, 26, 16, 17, 20, 29, 25, 13, 8, 24, 10, 30, 28,
22, 4, 12, 5, 15, 14, 11, 2, 6, 18, 23, 7, 21, 1, . . .
The jump function g 6 (x) = 16x mod 31 generates
1, 16, 8, 4, 2, 1, . . .
I.e., the first sequence is x0 , x1 , x2 , . . .; the second is
x0 , x6 , x12 , . . .

Discrete-Event Simulation: A First Course Section 3.2: Multi–Stream Lehmer RNGs 8/1
Using the Jump Function

First, compute the jump multiplier aj mod m (one time cost)

Then, g j (·) permits jumping from x0 to xj to x2j to . . .
The user supplies one initial seed
If j is chosen well, g j (·) can “plant” additional initial seeds
Each planted seed corresponds to a different stream
Each planted seed is separated by j calls to Random

Discrete-Event Simulation: A First Course Section 3.2: Multi–Stream Lehmer RNGs 9/1
An Example 4-Stream Sequence

x0
.... •
.......
..
.. ..
.. ..
.. . ..
.. .
. .
. .
.
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
x .
3j •.. (a, m) = (48271, 231 − 1)
.
•.. xj
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
..
. ..
.. .
.. ..
... ..
...... .....
•
x2j

Discrete-Event Simulation: A First Course Section 3.2: Multi–Stream Lehmer RNGs 10/1
Example 3.2.5: An Appropriate Jump Multiplier

Consider 256 = 28 different streams of random numbers

Partition the RNG output sequence into 256 disjoint
subsequences of equal length
Find the largest j < 231 /28 = 223 such that the jump
multiplier is modulus-compatible (needed to avoid overflow).
g j (x) = (48271j mod m)x mod m can be implemented via Alg
2.2.1
Then g j (x) can be used to plant the other 255 initial seeds [alg. 2.2.2
Possibility of stream overlap is minimized (though not
eliminated!)

Discrete-Event Simulation: A First Course Section 3.2: Multi–Stream Lehmer RNGs 11/1
Maximal Modulus-Compatible Jump Multipliers

s is the desired number of disjoint streams

Maximal jump multiplier : aj mod m where j is the largest
integer less than ⌊m/s⌋ such that aj mod m is modulus
compatible
Example 3.2.6: multipliers for (a, m) = (48271, 231 − 1) RNG

jump multiplier
# of streams s ⌊m/s⌋ jump size j aj mod m

1024 2097151 2082675 97070

512 4194303 4170283 44857
256 8388607 8367782 22925
128 16777215 16775552 40509

Discrete-Event Simulation: A First Course Section 3.2: Multi–Stream Lehmer RNGs 12/1
Library rngs

rngs is an upward-compatible multi-stream replacement for

rng
By default, provides 256 streams, indexed 0 to 255 (0 is the
default)
Only one stream is active at any time
Six available functions:
Random(void)
PutSeed(long x): superseded by PlantSeeds
GetSeed(long *x)
TestRandom(void)
SelectStream(int s): used to define the active stream
PlantSeeds(long x): “plants” one seed per stream
Henceforth, rngs is the library of choice

Discrete-Event Simulation: A First Course Section 3.2: Multi–Stream Lehmer RNGs 13/1
Example 3.2.7: ssq2 Revisited
(Single Serve Queue) ssq3.c
Use rngs functions for GetArrival, GetService

GetArrival Method
double GetArrival(void) {
static double arrival = START;
SelectStream(0);
arrival += Exponential(2.0);
return (arrival);
}

GetService Method
double GetService(void) {
SelectStream(2);
return (Uniform(1.0, 2.0));
}

Include "rngs.h" and use PlantSeeds(12345)

Discrete-Event Simulation: A First Course Section 3.2: Multi–Stream Lehmer RNGs 14/1
Uncoupling Stochastic Processes

Per modifications, arrival and service processes are uncoupled

Consider changing the service process to
Uniform(0.0, 1.5) + Uniform(0.0, 1.5)
Without uncoupling, arrival process sequence would change!
With uncoupling, the service process “sees” exactly the same
arrival sequence
Important variance reduction technique
Indeed, without uncoupling, the additional call to Uniform() in the service process
will cause changes in the GetArrival() on which there is a call to Uniform().

Discrete-Event Simulation: A First Course Section 3.2: Multi–Stream Lehmer RNGs 15/1
Single-Server Service Node with Multiple Job Types

Extend the single-server service node model from Chapter 1

Consider multiple job types, each with its own arrival and
service process
Example 3.2.8: Suppose there are two job types
1 Exponential(4.0) interarrivals, Uniform(1.0, 3.0) service
2 Exponential(6.0) interarrivals, Uniform(0.0, 4.0) service
Use rngs to allocate a different stream to each stochastic
process

Discrete-Event Simulation: A First Course Section 3.2: Multi–Stream Lehmer RNGs 16/1
Example 3.2.8: Arrival Process

Arrival Process
double GetArrival(int *j)
/* Index j corresponds to job type */
{
const double mean[2] = {4.0, 6.0};
static double arrival[2] = {START, START};
static int init = 1;
double temp;
if (init) {
SelectStream(0);
arrival[0] += Exponential(mean[0]);
SelectStream(1);
arrival[1] += Exponential(mean[1]);
init = 0;
}
if (arrival[0] <= arrival[1])
*j = 0;
else
*j = 1;
temp = arrival[*j];
SelectStream(*j);
arrival[*j] += Exponential(mean[*j]);
return (temp);
}

Discrete-Event Simulation: A First Course Section 3.2: Multi–Stream Lehmer RNGs 17/1
Example 3.2.8: Service Process

Service Process
double GetService(int j)
{
const double min[2] = {1.0, 0.0};
const double max[2] = {3.0, 4.0};
SelectStream(j + 2);
return (Uniform(min[j], max[j]));
}

Index j matches service time to appropriate job type

All four simulated stochastic processes are uncoupled
Any process could be changed without altering the random
sequence of others!

Discrete-Event Simulation: A First Course Section 3.2: Multi–Stream Lehmer RNGs 18/1
Consistency Checks

With appropriate changes to ssq2, steady-state statistics are

r̄ w̄ d̄ s̄ l̄ q̄ x̄
2.40 7.92 5.92 2.00 3.30 2.47 0.83

Obvious consistency checks: w̄ = d̄ + s̄ and l̄ = q̄ + x̄

Other consistency checks:
• Both job types have avg service time of 2.0 =⇒ s̄ = 2.00
• Net arrival rate should be 1/4 + 1/6 = 5/12 =⇒ r̄ = 12/5 =
2.40
• x̄ should be ratio of arrival to service rates
5/12
= 5/6 ∼
= 0.83
1/2

Discrete-Event Simulation: A First Course Section 3.2: Multi–Stream Lehmer RNGs 19/1