Optimizing static calendar queues.
The calendar queue is an important implementation of a priority queue that is particularly useful in discrete event simulators. We investigate the performance of the static calendar queue that maintains N active events. The main contribution of this article is to prove that, under reasonable assumptions and with the proper parameter settings, the calendar queue data structure will have constant (independent of N) expected time per event processed. A simple formula is derived to approximate the expected time per event. The formula can be used to set the parameters of the calendar queue to achieve optimal or near optimal performance. In addition, a technique is given to calibrate a specific calendar queue implementation so that the formula can be applied in a practical setting.
INTRODUCTION
The calendar queue data structure, as described by Brown [Brown 1988], is an
important implementation of a priority queue that is useful as the event queue in
a discrete event simulator. At any time in a discrete event simulator there are
active events, where each event e has an associated event time t(e) when it is
intended to occur in simulated time. The set of events is stored in the priority
queue ordered by their associated event times. A basic simulation step consists
of nding an event e 0 which has the smallest t(e 0 ), removing the event from the
priority queue, and processing it. As a result of the processing new events may
be generated. The parameter N can vary if zero or more than one new events are
generated. Each new event e has an event time t(e) > t(e 0 ) and must be inserted
in the priority queue accordingly.
In the calendar queue the events are stored in buckets with each bucket containing
events whose times are close to each other. All the events with the smallest times
are in the same bucket so they can be accessed quickly and simulated. Any newly
generated event can be quickly put into its bucket. When the events in one bucket
are consumed, the next bucket is considered. The details of the algorithm are given
later. The calendar queue has several user controllable parameters, the bucket
width and number of buckets, that aect its performance. Brown [Brown 1988]
provided empirical evidence that the calendar queue, with its parameters properly
set, achieves expected constant time per event processed. The goal of this paper is
to prove the constant time per event of the calendar queue behavior in a reasonable
model where, for each new event e, the quantity t(e) t(e 0 ) is a nonnegative random
variable sampled from some distribution.
Generally, the number of active events may vary over time. An important case
is the static case which arises when N is a constant, such as the case of simulating
a parallel computer. In this case, each event corresponds to either the execution
of a segment of code or an idle period by one of the processors. Thus, if there are
processors, then there are exactly N active events in the priority queue. In this
paper we focus on the static calendar queue.
Even before Brown's paper [Brown 1988] the calendar queue was used in discrete
event simulators when the number of events is large. In many of these situations the
calendar queue signicantly outperforms traditional priority queue data structures
[Brown 1978; Francon et al. 1978; Knuth 1973; Sleator and Tarjan 1985; Sleator and
Tarjan 1986; Vuillemin 1978]. An interesting new development is the employment
of a calendar queue like data structure as part of the queuing mechanism of high-speed
network switches and routers [Rexford et al. 1997]. In this case the calendar
queue like data structure is implemented in hardware.
1.1 Organization
In Section 2 we dene the calendar queue data structure and the parameters that
govern its performance. In Section 3 we present the Markov chain model of calendar
queue performance. In Section 4 we present an expression that describes
the performance of the calendar queue in the innite bucket case. The bucket
width can be chosen to approximately minimize the expected time per event at a
constant. In Section 5 we describe how to choose the number of buckets without
Optimizing Static Calendar Queues  3
signicantly compromising the performance over the innite bucket calendar queue.
In Section 6 we develop a technique for calibrating a calendar queue implementation
and demonstrate the eectiveness of the technique. In Section 7 we present
our conclusions. The Appendix contains the longer technical proofs.
2. THE CALENDAR QUEUE
A calendar queue has M buckets numbered 0 to M 1, a current bucket i 0 , a bucket
width -, and a current time t 0 . We have the relationship that i
For each event e in the calendar queue, t(e)  t 0 , and event e is located in bucket
if and only if i  t(e)=- mod M < (i 1). The analogy with a calendar can be
stated by: there are M days in a year each of duration -, and today is i 0 which
started at absolute time t 0 . Each event is found on the calendar on the day it is to
occur regardless of the year.
As an example choose 30. The 8 events
have times 31, 54, 85, 98, 111, 128, 138, 251.
{ 111 128 31 { 54 { { 85 98
In this example the next event to be processed has time 31 which is in the current
bucket numbered 3. Suppose it is deleted and the new event generated has time
87. Then, the new event is placed in bucket 8 next to the event with time 85. Since
will not be processed until the current bucket has cycled around all the
buckets once. Thus, t 0 is increased by -, and the next bucket to be examined is
bucket 4 which happens to be empty. Thus, the processing of the buckets is done in
cyclic order and only the events e which are in the current cycle, t 0  t(e) <
are processed.
A calendar queue is implemented as an array of lists. The current bucket is
an index into the array, and the bucket width and current time are either integers,
xed-point or
oating-point numbers. Each bucket can be implemented in a number
of ways most typically as an unordered linked list or as an ordered linked list. In
the former case insertion into a bucket takes constant time and deletion of the
minimum from a bucket takes time proportional to the number of events in the
bucket. In the latter case insertion may take time proportional to the number of
events in the bucket, but deletion of the minimum takes constant time. The choice
of algorithm for managing the individual buckets is called the bucket discipline. In
this paper we focus on the unordered list bucket discipline.
2.1 Calendar Queue Performance
For the calendar queue, the performance measure we are most interested in is the
expected time per event, that is, the time to delete the event with minimum time and
insert the generated new event. There are two key (user controllable) parameters
in the implementation of a calendar queue that eect its performance, namely, the
bucket width - and the number of buckets M . The choice of the best - and M
depends on the number of events N and the process by which t(e) is chosen for a
newly generated event e. Assuming M is very large (innite), if - is chosen too
4  K.B. Erickson, R.E. Ladner, A. LaMarca
N number of events known parameter
mean of the jump known or estimated parameter
time per empty bucket hidden parameter determined by calibration
time per list entry hidden parameter determined by calibration
d xed time per event hidden parameter determined by calibration
- bucket width user controlled parameter
M number of buckets user controlled parameter

Table

1. Parameters of the calendar queue.
large, then the current bucket will tend to have many events which is ine-cient. On
the other hand if - is chosen too small, then there will be many empty buckets to
traverse before reaching a non-empty bucket, which again is ine-cient. Regardless
of the choice of -, if M is chosen too small, then the current bucket will again tend
to have too many events in it which are not to be processed until later visits to the
same bucket.
In order to analyze the calendar queue we make some simplifying assumptions
on the process by which t(e) is chosen for a new event e. The main assumption we
make is that the quantity t(e) t(e 0 ), called the jump, is a random variable sampled
from some distribution that has a mean , where e 0 is the event with minimum
time t(e 0 ). We will fully delineate the simplifying assumptions later. The choice of
a good - certainly depends on both  and N . As  grows so should -. As N grows
should decrease. Determining exactly how - should change as a function of  and
N to achieve optimal performance is a goal of this paper.
Assume that we have innitely many buckets. In addition to the two parameters
and N the choice of a good - also depends on three hidden implementation
parameters b, c, and d where b is the incremental time to process an empty bucket,
c is the incremental time to traverse a member of a list in search of the minimum in
the list, and d is the xed time to process an event. If m empty buckets are visited
before reaching a bucket with n events (n  1), then the time to process an event
is dened to be:
KN (-) to be the stationary expected value of bm d. Then KN (-) is
the expected time per event in the innite bucket calendar queue.
In a real implementation of a calendar queue the number of buckets M is nite.
In this case, it may happen that some events in the bucket have times that are not
within - of the current time and are not processed until much later. Dene K M
to be the expected time per event in the M bucket calendar queue. Generally,
N (-)  KN (-) because extra time may be spent traversing events in buckets
that are not processed until later. Another goal of this paper is to determine how
to choose M so that K M
N (-) is the same or only slightly larger than KN (-).

Table

1 summarizes the various parameters that aect the performance of the
calendar queue.
Optimizing Static Calendar Queues  551525
time
per
event
bucket width
Fig. 1. Graph of bucket width, -, vs. the expected time per event, KN (-), in the simulated
innite bucket calendar queue with 100 events, exponential jump with mean 1, and
time
per
event
number of buckets
Fig. 2. Graph of number of buckets, M , vs. the expected time per event, K M
N (-), for - chosen
optimally in the M bucket simulated calendar queue with 1,000 events, exponential jump with
mean 1, and

Figure

1 illustrates the existence of an optimal - for minimizing the expected time
per event. Figure 2 illustrates the eect of selection of M on the expected time
per event. The graphs in both Figures were generated by simulating the calendar
queue with an exponential jump with mean 1 and
were taken after a suitably long warm-up period and over a long enough period so
that average time per event was very stable. The simulation of Figure 1 uses an
innite number of buckets with 100 events. The simulation of Figure 2 uses the
optimal bucket width for for the innite bucket calendar queue, then
varying the number of buckets. In choosing or 3; 000 the performance
curve is almost
at approaching the performance with innitely many buckets.
6  K.B. Erickson, R.E. Ladner, A. LaMarca
3. MODELING THE CALENDAR QUEUE PERFORMANCE
To model the calendar queue performance we begin by specifying the properties of
the random variable, is the event with current minimal
time, and e is the newly generated event. We assume that  is a random variable
with density f dened on [0; 1), the nonnegative reals. We call f the jump density
and its random variable simply the jump. Successive jumps are assumed to be
mutually independent and identically distributed. Let  be the mean of the jump,
that is:
Z 1f1 F (z)g dz: (1)
where F is the distribution function of the : F
R xf(z) dz. We call F the
jump distribution.
We dene the support of the jump distribution to be
The value 1 is not excluded.
Technical Assumptions about f and F .
In order to facilitate the proofs we make several technical assumptions about f and
F that will be in force throughout, except as noted.
J1. The density f(x) > 0 for all x in the interval (0; ).
J2. The mean  is nite.
J3. There is an  0 > 0 and c 0 such that F (x)  c 0 x for all x   0 .
Assumption J2 is crucial; it guarantees the existence of a non-trivial \steady state".
Note that J3 holds if the density f is bounded in a neighborhood of 0.
3.1 The Markov Chain
We model the innite bucket calendar queue as a Markov chain b
X with state space
in [0; 1) N . For denote the state of
the chain at time t. The state (X 1 represents the positions,
relative to the beginning of the current bucket, of the N events (indexed 1 to N)
in the calendar queue at step t. A step of the calendar queue consists of examining
the current bucket, and either moving to the next bucket, if the current bucket
is empty, or removing the event with smallest time from the current bucket and
inserting a new event (with the same index) according to the jump distribution.
Accordingly, the transitions of b
are as follows: Let m be the index such that
for all i. If X m (t) < -, then for i
are independent non-negative random variables. It is assumed that
these random variables  t ; t  0, all have the same probability density f . The
parameter - is a xed non-negative real number.
We can think of X i (t) as the position of the i-th particle in an N particle system.
If no particle is in the interval [0; -), then all particles move - closer to the origin.
Otherwise, the particle closest to the origin, in the interval [0; -), jumps a random
Optimizing Static Calendar Queues  7
distance from its current position and the other particles remain stationary. Thus,
a particle in the Markov chain b
X represents an event in the innite bucket calendar
queue where the position of the particle corresponding to an event e is the quantity
. The interval [0; -) corresponds to the currently active bucket in the innite
bucket calendar queue.
It is important to note that a step of the Markov chain b
X does not correspond to
the processing of an event in the calendar queue. The processing of an event in the
calendar queue corresponds to a number of steps of the Markov chain where the
interval [0; -) is empty followed by one step where the interval [0; -) is nonempty.
to be the limiting probability, as t goes to innity, that the interval [0; -)
has exactly i particles in it. Technically, q is a function of N and -, but
we drop the N and - to simplify the notation. The quantity q 0 is the probability
that the interval [0; -) is empty. It is not obvious that q i exists for 0  i  N , so
we prove the following Lemma in Appendix A.
Lemma 3.1. If the jump density has properties J1, J2, and J3, then the limiting
probabilities exist and are independent of the initial state of b
X.
Let us also dene EN (-) to be the limiting expected number of particles in the
interval [0; -), that is,
4. EXPECTED TIME PER EVENT IN INFINITE BUCKET CASE
The expected time to process an event in the innite bucket calendar queue is
closely related to the function EN (-) as we see from the following Lemma.
Lemma 4.1. The expected time per event in the innite bucket calendar queue
is
Proof. The Markov chain b
X models the calendar queue. Thus q 0 is the portion
of buckets visited which are empty and for j > 0, q j is the portion of buckets visited
which have j events. Each empty bucket visited, which happens with probability
cost b, but does not result in nding an event to process. Each bucket
visited with j > 0 events, which happens with probability q j , has cost cj + d, and
results in nding an event to process. Thus, the expected cost per event in the
calendar queue is
d)
which yields equation (3) using equation (2).
8  K.B. Erickson, R.E. Ladner, A. LaMarca
Let us dene the following important quantity
Z -[1 F (x)] dx: (4)
The second part of equation (1) implies that 0  p  1. Note also that -[1
F (-)]=  p  -=.
In order to derive a good approximating formula for KN (-) we rst need to nd
good bounds for the quantities q i for 0  i  N . The following technical Lemma,
proved in the Appendix, Section B, provides those bounds.
Lemma 4.2. For N  2 and all - > 0 we have
and for
where B(j) is the tail of the binomial distribution for N trials with \success" parameter
p:
The simple exact formula for q 0 is interesting. It is possible to write down some
very complicated integrals which give exact expressions for the other q j , but these
are highly unwieldy and their proofs are not informative (cf. [Erickson
1999]).
It is also interesting to note that our assumption J1 requiring the probability
density f to be positive on its support can be removed, but the proofs of the
theorems become even longer. Without J1, if the density f has the property that
there is a constant c > 0 such that we have
exact expressions for q j for all j  1.
Lemma 4.2 yields the following upper and lower bounds on KN (-).
Lemma 4.3. For N  2 and all - > 0
KN (-) d  b
Proof. From (5) and (6) we have
+N-
to be the standard binomial distribution with N trials
and success parameter p. Since and b(i) has mean Np and
Optimizing Static Calendar Queues  9
second moment (Np) 2 +Np(1 p) we sum by parts to derive
which, upon substituting into equation (3) and doing a little rearranging, yields the
left side of (8). Similarly, one derives the right side of (8).
The range of - that gives good calendar queue performance is when
In this case the bounds of (8) give us a wonderfully simple, and accurate, approximating
formula for KN (-).
Theorem 4.1. If then the expected time per event in the innite
bucket calendar queue with bucket width - is
In fact, there are numbers  1 ;  2 such that for any xed
> 0 the O(N 1 ) term is
bounded by ( 2
)=N uniformly for 0 < -
=N .
The proof of this Theorem is almost an immediate consequence of (3), (5), and
but it is also postponed to the Appendix, Section C. Interestingly, the expected
time depends on the mean of the jump and not on the shape of its probability
density.
Note that one immediate consequence of Theorem 4.1 is that if the bucket width
is chosen to be =N for  in a xed interval, then the innite bucket calendar
queue has constant expected time per event performance. Indeed, a formula for the
optimal performance of the calendar queue can be derived as seen in the following
Theorem.
Theorem 4.2. The expected time per event, KN (-), achieves a global minimum
in the interval (0; 1) at - opt where
r
c
KN (- opt
2bc +O N 1
The proof of this Theorem is in the Appendix, Section D. Theorem 4.2 shows that
the optimal choice of - only depends on the ratio of b to c, the mean  of the jump,
and N .
5. CHOOSING THE NUMBER OF BUCKETS
Now that we have found how to select - so as to approximately minimize the
expected time per event in the innite bucket calendar queue our next goal is to
select M , the number of buckets, so that the M bucket calendar queue has the
same or similar performance as the innite bucket calendar queue.
For the case in which the jump distribution has nite support ( < 1), there is a
natural choice for M which guarantees that the calendar queue with M buckets has
A. LaMarca
exactly the same performance as the innite bucket calendar queue. If M  =-+1,
then it is guaranteed that in the long run all the events e in the current bucket will
have In this case, eventually each event in the current bucket will
be processed during the current visit to the bucket and not postponed until future
visits to the bucket. For the case in which the support  of the jump distribution
is either innite or is nite but =- is too large to be practical, then it will be
necessary to choose a number M which gives performance less than that of the
innite bucket calendar queue.
5.1 Expected Time per Event in the nite Bucket Case
The same Markov chain b
X can be used to analyze this case. Let L M
N (-) be the
(steady state) expected number of particles in the set
In terms of the M bucket calendar queue, if an event e has t(e) t 0 2 , then
the event is in the current bucket but is not processed. The occurrence of such an
event will cause the M bucket calendar queue to run less e-ciently than the innite
bucket calendar queue. The following Lemma quanties the dierence between the
performance of the nite and innite bucket calendar queues.
Lemma 5.1. The expected time per event in the M bucket calendar queue with
bucket width - is
Proof. In the Markov chain b
X, let q ij be the limiting probability that there are
particles in the interval [0; -) and j particles in . The probabilities q ij can be
shown to exist in the same way as we did for the probabilities q i in Lemma 3.1 by
using Corollary A.1 in the Appendix, Section A. In the M bucket calendar queue
the cost of visiting a bucket with i events whose times are in the interval [t
and j events whose times are in the set
Thus, the expected cost per event K M
But
(equation (5) of Lemma 4.2) and
(-). By using equation (3) in the proof of Lemma 4.1 we derive the equation
for K M
(-).
In the

Appendix

, Section E, we indicate how to derive the following rather horrible
looking bounds for L M
(-).
Lemma 5.2. The function L M
N (-) is bounded above by
Optimizing Static Calendar Queues  11
and bounded below by
+N-
where p and F have the same meaning as before (see equation (4)) and
Z jM-
[1 F (x)] dx;
[F (jM- y) F (jM- y)]:
It should be noted that under the hypothesis  < 1 (J2), the above series converge
and can be given bounds in terms of ; -, and M . However, using the bounds as
stated in the Lemma, we can derive a more useful asymptotic expression for L M
(-).
The Lemma is proven in the Appendix, Section F.
Lemma 5.3. If are constants, then
[1 F (jxr)]: 1
5.2 Degradation in Performance due to Finitely Many Buckets
M to be the degradation in performance in choosing M buckets instead of
innitely many buckets, that is,
If we choose - optimally, then Lemma 5.3 and Theorem 4.2 yield the following
asymptotic expression for  M .
Theorem 5.1. If M=N is constant and
c
[1 F (jM-)] (12)
The following asymptotic bound is implied by Theorem 5.1.
Theorem 5.2. If M=N is constant and
c
r c
12  K.B. Erickson, R.E. Ladner, A. LaMarca
Proof. By Theorem 5.1 it su-ces to show that
To see this let k  2, then
Z 1xf(x)dx
Z jD+D
x f(x)dx
The niteness of  implies that x[1 F (x)] ! 0 as x ! 1. Therefore, if we let
k !1, we get =D
Equation (13) shows that for a xed  > 0 (like .01) M can be chosen to be O(N)
so that  M  . In other words, one can always choose the number of buckets M
to be a multiple of N and still obtain a performance almost as good as that of the
innite bucket case.
For the interesting case of the exponential jump density
we can calculate the series in the equation (12) exactly:
e
c
Let us suppose optimally equal to
2=N allows us
to solve equation (14) for M=N when given an acceptable  M . For example, if we
choose  should be approximately 1:92 and if
M=N should be approximately 3:02. Figure 3 illustrates that asymptotic equation
provides an excellent choice of M over a wide range of N . Using our simulation
of the calendar queue we plot for a wide range of N the value of  M for each of
3:02. Again, measurements were taken after a suitably
long warm-up period and over a long enough period so that average time per event
was very stable. Both plots are relatively
at near the asymptotic values .05 and
respectively. Thus, equation (14) seems quite accurate. The bound of Theorem
5.2 is not necessarily tight because we are crudely approximating an integral. For
example, if we choose  then the formula (13) requires M=N to be at least
Optimizing Static Calendar Queues  130.030.07100 1000
degredation
Fig. 3. Graph of N vs degradation,  M , from simulations for
6. CALIBRATION OF A CALENDAR QUEUE IMPLEMENTATION
In an actual calendar queue implementation we would like to nd the best bucket
width - and number of buckets M . The preceding theory tells us how to do so
if we know the hidden implementation parameters b, c, and d. In this Section we
give a relatively simple method of estimating these parameters simply by timing
executions of the simulation for various values of - proportional to =N . The key
to the method is equation (9) for the expected time per event. We can write KN (-)
as a linear function of the unknowns b, c, and d. The general calibration method is
as follows: rst estimate M to be large enough so that the degradation in using M
buckets over innitely many is small. Second, nd K M
N (-) for a number of dierent
-'s by timing executions of the implementation, and third, use a linear least squares
approximation to nd the b, c and d that best ts the function
We illustrate this method with an example. We developed a calendar queue implementation
in C++ and ran it on a DEC alphastation 250. We chose
and an exponential jump with mean 10; 000. Just by examining the code we felt
that b, the time to process an empty bucket, was considerably larger than c, the
cost of traversing a list entry. We made an educated guess that the optimal - was
certainly greater than 5. We chose
or larger there is only a small chance that an event in the current bucket is not
processed because its time is too large. We timed the calendar queue for 20 values
of - ranging over several orders of magnitude, namely, Using
this data we used linear least squares approximation to compute
using equation (15). Figure 4 shows the curve of
equation (15) using these parameters. The Figure also shows the time per event for
200. Thus, this method accurately predicts data points that were
not used in the linear least squares approximation. It is interesting to note that
using these values of b, c, and d in equations (10) and (11), we obtain - opt  61:5
and K(- opt )  1756:38. By contrast, the best - among the 40 executions is
14  K.B. Erickson, R.E. Ladner, A. LaMarca
with execution time 1754:92.1800220026003000
time
per
event
bucket width
measured
predicted
Fig. 4. Measured and predicted expected time per event for a calibrated implementation of a
calendar queue.
Care must be taken in applying this calibration method. In the method, the
hidden parameters, b, c, and d, are measured indirectly by measuring the expected
time per event. For xed M , N , and -, the expected time per event can vary over
dierent runs because of interruptions by other processes, page faults, or other ef-
fects. However, in our experimental setting we carefully controlled the environment
so that our running times varied little for a xed parameter setting. In addition,
measurements were taken after a suitably long warm-up period and over a long
enough period so that average time per event was very stable. In a real computing
environment that cannot be controlled this calibration method might not yield such
good results.
Ideally, using a xed N , M (large enough), and  we can estimate the hidden
parameters b, c, and d which then could be used for any other N , M (large enough),
and , and -. However, because of the cache behavior of modern processors, the
values of b, c, and d are not actually constant independent of M , N , and properties
of the jump distribution other than its mean. For example, a smaller M
might achieve fewer cache misses reducing the running time and thereby eectively
lowering the values of these constants. It may be that in applying the calibration
method, - is chosen so large that the original M chosen is far larger than necessary.
In this case, it might be wise to choose a smaller M , then recalibrate the calendar
queue starting with a larger -.
In a real application of the calendar queue it is unlikely that the jumps are mutually
independent, identically distributed random variables as described in our
model. Nonetheless, the mean of the jump can be empirically estimated, the calibration
done, and equation (10) for the optimal - applied to nd a potentially good
-.
Optimizing Static Calendar Queues  15
7. CONCLUSION
We have shown that there is an expression for the expected time to process an event
in the innite bucket calendar queue and that the bucket width can be chosen
optimally. With the bucket width near the optimal bucket width the calendar
queue has expected constant time per event. The optimal bucket width depends
only on a few parameters, the incremental time to process an empty bucket (b),
the incremental time to traverse a list item (c), the mean of the jump (), and the
number of events (N ). We have shown that the number of buckets M can be chosen
to be O(N) so as to achieve minimal or almost minimal expected time per event.
Finally, we have shown that the implementation parameters can be determined by
using approximation based on the method of linear least squares.
Although the calendar queue runs very fast for certain applications it has the
disadvantage that its performance depends on the choice of parameters - and M .
An interesting problem would be to design a priority queue based on the calendar
queue that automatically determines good choices for - and M . We believe that
the calibration method described in this paper might give insight into the design
of a dynamic calendar queue where N and/or  can vary over time.




Section A of the Appendix sets up the notation and concepts that are used through-out
the Appendix.
A. INVARIANT DISTRIBUTION, POSITIVITY, AND LIMITS
Consider the Markov chain b
described in Section 3. The symbol P  shall denote
the probability measure induced on the trajectory space of the chain b
X when the
initial distribution is , and P b x shall denote trajectory space probabilities when
the chain starts at the point b x. (Note: P
R P b x d(b x), the integration being
carried out over the entire state space.) Integration (better known as expectation)
with respect to P  and P b x is denoted E  and E b
x , respectively.
stand for the set of points b such that
and let A . For an b x in the state space and a (measurable) subset A, the
one-step transition probability (T. P.) that the chain will move from b x to a point
in A is given by
(b
Z 11A (b x
the standard i th unit coordinate vector, and 1A (b x)
is the function which is 1 for b x in A and 0 otherwise. 3
3 Actually, (16) does not dene a proper transition probability at all points! Indeed, if b x is any
point with two or more coordinates that are equal and strictly less than -, then b x does not lie in
A 0 nor in B i for any i. Hence P (b x; subsets A. However, the jump distribution F
has a density so it has no atoms (discrete points of positive probability). This and the dynamical
description of the chain in Section 3 imply that two points which start at the same position
A. LaMarca
Let  be the set [0; - It follows from (16) that
(b in  and A   C .
In other words,  is an absorbing set for the chain. Moreover, the dynamical
description of the chain implies that a particle which starts outside of  will reach
in a nite (but possibly random) number of steps. (One can show, using the
method in the proof in B.3 that the number of steps required to eventually enter
has a nite expectation.) Thus  C is transient for the chain. (Of course,
measure m is an invariant measure for the chain if m is -nite and for every
measurable subset A of the state space
Z
(b x; A)mfdb xg; (17)
A Markov chain is called a Harris recurrent chain, or simply a Harris chain, if
there exists a unique, up to positive multiples, invariant measure m such that if
A is any Borel subset with m(A) > 0, then P b x ( b
x in the
state space. (The initials i.o. stand for \innitely often".) A Harris chain with an
invariant probability measure, necessarily unique, is called positive.
The state space of a Harris chain can be written as a disjoint union:
6. The sets C are
known as the recurrent cyclic classes. The integer d is nite and if d = 1, the chain
is called aperiodic.
Theorem A.1. If the jump density satises J1, J2, and J3, then the Markov
chain b
X with T.P. (16) is a positive, aperiodic, recurrent Harris chain. Its invariant
probability is concentrated on
Corollary A.1. If ' is any bounded measurable function, then for any initial
distribution , we have
lim
Z
'(b x) dm(b x); P  a.s. (18)
Note: The left-most term is a limit of averages of random quantities and the assertion
is that the limit exists with P  -probability 1 and equals the (non-random)
quantities on the right. If ' is unbounded but integrable with respect to m and
if the chain is Harris and positive, then the above limit relations remain valid at
least in the case that  is point mass at some b x or that  = m.
Proof. of Corollary A.1. The deterministic limit statements of Corollary A.1
are immediate consequences of Proposition 2.5 in Ch.6, x2, of [Revuz 1984]. The
eventually become and remain separated w. p. 1. Indeed, this occurs as soon as one of them
jumps to the right. Thus, it does no harm if we banish such points from the state space initially.
With this understanding P is indeed a transition probability on its state space.
Optimizing Static Calendar Queues  17
a.s. limit-of-averages assertion is a consequence of the ergodic Theorem for Harris
chains. See Theorem 4.3, and its companion remark, in [Revuz 1984], Ch. 4, x4.
Remark. The main signicance of aperiodicity is that it justies the existence of
the limit of E b x f( b
occurring in (18). The existence of limits of averages does
not require aperiodicity.
Proof. of Lemma 3.1. Let
exactly j components of b
x lie in [0; -)g:
the number of particles in interval [0; -) at time t
and from (18) it follows immediately that q g.
The proof of Theorem A.1 will be postponed to the very last Appendix below.
It is lengthy and somewhat tedious, but there are some interesting features.
B. PROOF OF LEMMA 4.2
B.1 The computation of q
In this Section we prove (5) of Lemma 4.2.
The sets B i , dened in the last Section, are disjoint and their union is the complement
(in ) of A 0 . Since m assigns 0 mass to [0; 1) N \  c , we have
Let be any bounded or positive function on the state space. Equation (17) has
an analogue for functions which reads:
R (b x)mfdb
R mfdb xg
R (b y)P (b x; db y).
Noting that
R
R (b y)P (b x; db
R
A0 (b x - b 1)mfdb xg, by (16), and doing
a little rearranging, we get
Z
(b x - b 1) (b x)
Z
mfdb xg
(b x)
Z
(b y)P (b x; db y)
Fix i and let (b complex number with nonnegative
real part. Then for b x in B i ,
Z
(b y)P (b x; db
R 1e z f(z) dz is the Laplace transform of F . Also, for b
x in B j with
Z
(b y)P (b x; db
All but the i th term on the right side of (20) vanishes and it becomes
[1 ()]
Z
A. LaMarca
Simplication of the left-hand side of (20) leads to:
Z
Z
Divide (21) by  and make  ! 0. The result is -mfA 0 g=  0 (0)mfB i
Equation (5) follows immediately from this and (19).
B.2 The Case
If we observe the successive positions of a single one of our N particles at only
those times at which it actually moves, we get a 1-dimensional version of the N -
dimensional chain. For
From the description of the chain in terms of the independent random variables ,
one concludes:
(i) Each sequence fX i is itself a Markov chain on the line;
(ii) These N Markov chains are mutually independent.
If one can nd an increasing sequence of times fS k g such that each S k is a common
value of every one of the u i (that is, for each k there are numbers r i (k) not
necessarily the same, such that S
mutually independent components. Here is such a sequence: let S
and for k > 0, let
The times T are the successive (random) times at which the
interval [0; -) is empty of particles (Z(T k
is obtained from
by adding the deterministic constant - to each of the components of b
it follows that the components of b
are also mutually independent. It is easy to
show that the chain induced on A 0 (or trace chain), the sequence f b
is also a Markov chain. See [Revuz 1984], Exercise 3.13, page 27.
An important point to note is that the special structure of b
X implies that for
each i the chain fX i (T k ); k  0g , coincides in law with the trace chain on [-; 1)
of an
X. It is at least intuitively clear that the trace chain
positive recurrent and has an invariant probability distribution m 0 ,
say, obtained by renormalizing the distribution m restricted to A 0 . (See [Revuz
1984], Ex. 3.13, p. 27, and Prop. 2.9, p. 93 for a formal proof.) Thus for subsets B,
But because this trace chain also has independent components, it follows that m 0
is a \product measure" built up from the invariant distributions of each of its
component chains. These component chains have identical T.P.'s, so the factors in
are the same. Let us call this common factor distribution m 10 . Once computed,
(concentrated on [-; 1)) may be used to compute the limit, as k !1, of the
Optimizing Static Calendar Queues  19
probability of nding exactly j particles in the interval [0; -) at the times S
+1. By now it should be clear that the limiting distribution of Z(S k ) is a binomial
distribution corresponding to N Bernoulli trials with parameter
(However, the limit distribution of Z(t) for t tending to innity without restriction
is not a binomial.)
The invariant distribution, let us call it m 1 rather than m, in the case
our basic chain can be calculated explicitly and then m 10 obtained from the special
case of (22). The measure m 1 turns out to be uniform on [0; -) and coincides
with the (-translate) of the stationary distribution for the renewal process with
interarrival distribution F . This stationary distribution has a density equal to the
normalized tail-sum 1 F . See [Feller 1971], XI.4. One can give queuing theory
arguments for the above description of m 1 , but, since equation (21) leads to this
result almost immediately, we use that equation to give a quick proof. In the case
equation (21), simplies to
valid for any complex number ; Re()  0. If we set is an
arbitrary integer, we nd that the left-hand side vanishes. The density assumption
implies that () 6= 1 for any  6= 0. Hence
results in the theory of Fourier series implies that we
must have dm 1 some constant C. From (5) in the case
-). For the Laplace
transform of m 1 on [-; 1), we get
Z -e x
Inverting the Laplace transforms in this equation reveals that the density g 1 of m 1
on [-; 1), for x  -, is given by g 1 From (22) it is then
clear that m 10 has the density
The upshot of the preceding is that we can now conclude that the limit distribution
of Z(S k ) is:
lim
(w.p.1), for
x
Remark. It follows from the work of the last two Sections that the measure m,
when restricted to A 0 , is a product measure because its restriction to A 0 coincides
with However, the m-measures of subsets of A 0 have no particular interest;
it is only the m-measures of the other A j (dened in the proof of Lemma 3.1) that
is required and these sets are contained in the complement of A 0 . But m restricted
to the complement of A 0 is not a product measure.
20  K.B. Erickson, R.E. Ladner, A. LaMarca
B.3 Estimates for the q j
Apart from the explicit representation of m on A 0 discussed in the last Section,
a simple expression for m on all of [0; 1) N for N > 1 is not available. This
means that, with one exception (the case F (-) = 0), we do not have simple explicit
formulae for the values of q and must resort to approximations. It
turns out, however, that the approximate formulae are quite amenable to analysis
particularly in the region of interest:
In this Section we nish the proof of the two inequalities of (6) which we henceforth
designate (LH-6), for the left side, and (RH-6) for the right.
To simplify the notation a little, the starting distribution  will be omitted, if
not forgotten, when it is not essential. For this proof we introduce the objects
I A j ( b
I A  j
where A
. The variable #(t; 0) diers from n - (t) by at most 1 because
the T k 's are the zeros of Z. Therefore, by the ergodic limit theory for b
lim
Hence
Sn
Sn
and then
Sn
Sn
a sequence of random variables fV k g (k  1) by V
shall be the total number of times t in the interval
that the counting variable Z(t) has the value j:
I A j ( b
The occurrence of the event V k > r implies that Z(S k 1 and that there are
at least r jumps, counting from the rst time in [S k 1
magnitudes smaller than -. Hence
Optimizing Static Calendar Queues  21
and therefore, assuming F (-) < 1,
(Indeed all moments, Ef(V k )  g, are nite.) This inequality, (23),
and the ergodic limit theory now yield
lim
lim
which is (RH-6).
As to (LH-6) note rst that
Hence, see (25), and (23),
lim n
which is (LH-6).
C. PROOF OF THEOREM 4.1
For the purposes of this proof let us write
The conclusion of Theorem 4.1 is equivalent to the assertion that
KN (-) d
x
c
x
uniformly on bounded x-intervals. We base the proof on (8) which states
c
x DN (-)  KN (-) d b
x
in the new notation. Fix a number
> 0 and conne x to the interval (0;
=] so
that 0 < -
=N . Let  0 and c 0 be the numbers introduced in assumption J3 in
Section 3. Keeping N > maxf2c 0
gets
22  K.B. Erickson, R.E. Ladner, A. LaMarca
Also DN (-)  1(N 2
=N.
Hence
x
c
x DN (-)
From this, and a little algebra, it is easily seen that to
nish the proof of (28), it su-ces to nd a number C 2 , depending on
, such
But, because Np  N-
x
x
x
The numbers C 1 (multiplied by c) and C 2 yield estimates for  1 and  2 mentioned
in Theorem 4.1:
D. PROOF OF THEOREM 4.2
Throughout this Section we will write a = b=c,
Moreover, there is no harm in also supposing that
Step 1. As a function of -, KN (-) is continuous on (0; 1). The reader is asked
to turn to the formula (3). To begin with the variable q 0 is =( + N-) which is
obviously continuous. The only possible discontinuous term in the formula for KN
is EN (-). However, the continuity of this function is an immediate consequence of
the following exact formula which will be discussed after the proof of the Theorem
is complete.
Lemma D.1.
Step 2. The function - 7! EN (-) is nondecreasing. One can prove this by
dierentiating the expression for E(-) of Lemma D.1 and checking that the result
is non-negative.
Optimizing Static Calendar Queues  23
Here is an outline of an alternative, but more intuitive, proof: Consider two
chains b
with the same N and jump density but with dierent bucket
. If we follow the trajectory of an individual particle in each chain,
then we nd that on average in the chain with the larger bucket size, - 2 , the particle
gets back to the interval [0; - 2 ) quicker than it would get back to the interval [0;
in the chain with the smaller bucket size, - 1 . Since a typical particle of the chain
more often in the interval [0; - 2 ) than in [0; - 1 ) of the chain b
the average
number of particles in [0; - 2 ) of b
2 is at least as large as the average number of
particles in [0; - 1 ) of b
Step 3. We next establish
lim
uniformly on [-
for each xed - 1 > 0. The explicit formula for q 0 yields that
nondecreasing function of -. By equation (3) and step 2, we nd
that the function (1 q 0 (-))KN (-) q 0 (-)b is also nondecreasing in -. Hence for
b:
For 1. From the
above inequality we have for - 1  - 2 and N  =- 1 ,
From this inequality, it follows that KN (-) goes to innity uniformly in the interval
for each xed - 1 > 0. This follows from
equation (8).
Step 4. For
2a; we have min 0<-<
Moreover, if - 1 is the minimizing - on this interval, then
where - o;N is dened above. To prove all this, let
By straightforward calculus for each N , HN (-) has a global minimum on (0; 1) at
the point - o;N . Let - opt be a value of - which gives the minimum value of KN (-) on
the given interval (0;
=N ]. By Theorem 4.1, on this interval we can nd a constant
C, depending on
, such that for all N su-ciently large,
uniformly for 0 < - <
=N . Since
is in the interval (0;
=N ]. On this
interval, KN (-) is thus sandwiched between the two convex functions, HN (-)C=N ,
both of which have a global minimum at the same point - o;N interior to the interval.
For a xed N , - opt must be between the two solutions to the equation (in -) HN (-)
. By a simple calculation one nds that the dierence
A. LaMarca
between the two solutions is O(N 3=2 ) and this yields -
any - between the two solutions one nds that KN
Step 5. The next step is to show that for any
0 , there is
such that for all N  N 1
KN
Thus, the minimum exhibited in step 4, extends to the xed interval (0; - 1 ]. This
fact and step 3 imply that for all N su-ciently large,
min
completing the proof of the Theorem.
For the moment we x - 1 such that 0 < F (- 1 ) < 1. We will choose - 1 later. By
the inequality (8) and the fact that p(-[1 F (-)]= we have for all - 1
KN (-)  cQ a=z
-. Dene LN
.
For each N , the horizontal line at height K cuts the graph of the convex
function LN at two points, the larger of which we call z
N . Thus, z
N is the larger
root of the equation
As a sequence in N , the values z
N converge to a bounded positive limit. Dene
N =QN . Hence, the sequence N-
also converges to a limit,
. Note that
tends to
2a as Q approaches 1. We choose - 1 (and hence Q) so that
<
Now, choose N 1 such that -
Since the minimum
of KN (-) is bounded above by K bounded below by the function
LN (QN=(-)), then the minimum of KN (-) in the interval (0; must already lie
in the interval (0; -
N ], and hence in the interval (0;
D.1 A discussion of the exact formula for EN (-).
The proof of this result is quite long and is based on some exact, though very
complicated, integral formulae for the q j 's. See [Erickson 1999] for the details. The
exact formula for EN (-) leads to an exact formula for KN (-), but our work has
led us to the conclusion that the excellent asymptotic formulae of Theorems 4.1
and 4.2 (and the simple inequalities of Lemma 4.2 which lead to them) are of much
greater practical use and are certainly easier to prove. For this reason we have not
included the long proof of the exact formula.
Our main use of Lemma D.1 was to shorten the proof, slightly, of the global
minimization of KN . Note that D.1 yields, immediately, the continuity of KN as a
function of -. One requires continuity in order to speak sensibly of the existence of a
minimizing -. Even without the continuity, however, the basic result of Theorem 4.2
is essentially correct; only the language used to express it needs to be changed. (One
must use the term \greatest lower bound" in place of \minimum" and one can only
assert that there are points - at which the greatest lower bound is approximately
attained.)
Optimizing Static Calendar Queues  25
We will write L M
N (-) and Z A (t) for the number of particles in set A at time
t. is a subset of [-; 1) we have
are the successive times at which
the interval [0; -) is empty of particles.) Hence,
st
Z
Z
w.p.1, where  Suppose that at time T k 1 there are
particles at positions x 1 2-). Then at there will be r
particles in [0; -) at positions x
where according as the i-th of these particles lands in or not when
it is nally removed from the interval [0; -). (This removal must occur during
Writing by the strong Markov property,
where F
is the  eld of the random variables T k
and where H t (b) is the probability that a particle starting at the origin lands in the
interval it rst jumps over t. This H satises
where U is the renewal measure. See [Feller 1971], page 369. 4
In general U(z)  Uf[0; z]g  [1 F (z)] 1 for distributions on [0; 1) so that
sup
[F (jM- x z) F (jM- x z)] Ufdzg
(Recall the denition of  2 in the statement of Theorem 5.2.) Calling the right-hand
side p  and noting that conditional on the -eld F
, the variables u i are
4 Feller denes H in terms of the open interval (-; 1) whereas we are using the closed interval.
Because F has no atoms, this dierence in denition has no consequence.
26  K.B. Erickson, R.E. Ladner, A. LaMarca
independent of  k , we have
[Z +-
[Z +-
r
[1 F (-)] 1 Z +-
, and Z fg
means Z fg (T j k. Now at the times fT j g
the particles are independent so the limiting joint distribution of Z +-
is a trinomial. (Separately, they have binomial limit distributions.) Letting k !1
we get EfZ [-;2-)
Z +-
where p is dened at (4) and  Going back to (32) with these
calculations we obtain
Z +-
Z [-;2-)
(Recall the basic property of conditional expectations EfE[  j F
Replacing p  with  2 =[1 F (-)] and q 0 with =( +N-) combining fractions and
dropping the factor 1 F (-)] which will occur in the numerator, we nally obtain
the upper bound on L M
N .
For the lower bound we have
which evaluates to the lower bound on L M
N .
F. PROOF OF LEMMA 5.3
In the following we let lim, lim sup, and lim inf stand for the limits of various
quantities as N !1 with the other variables constrained to vary as stated in the
hypothesis.
First let us note that lim F
jrx. Note that t j does not vary with N . For all su-ciently
large N , - will be so small that the intervals f(t j -; t j +-]g 1
are non-overlapping.
-], then for all N we have JN+1  JN : Hence
by continuity of F (no atoms). (The letter F stands for both the distribution
function and the induced probability measure as is customary.) From this (and the
limits
Optimizing Static Calendar Queues  27
Next
[1 F (jrx)]:
As the rst sum on the left goes to
it then follows that
[1 F (jrx)]
Using these limits in the upper bound for L M
N (-) we get
lim sup L M
[1 F (jrx)]:
Similarly, from the lower bound for L M
mu+N-
[1 F (jrx)]:
Thus the lim sup L M
so the limit of L M
exists and its value is as stated.
G. PROOF OF THEOREM A.1.
We have seen that for each k  1 and given initial positions, the N components of
are mutually independent random variables. Also, it is not hard to see that
the conditional distribution of X the distribution
of the residual waiting time at epoch - of a delayed renewal process starting at
epoch x with interarrival distribution F . See [Feller 1971], page 369, and [Erickson
1999]. (We have already seen this distribution in the proof of Lemma 5.2, Section
E, though it was was described in slightly dierent language).
Letting H s fIg denote the probability that the residual waiting time at epoch s
lies in I for a pure renewal process starting at 0 (H s f[0;
E), it follows that for any xed x > 0 and every integer k >> x=- and any Borel
I  [0; ),
Using the Markov property it then follows that for xed b
Borel sets I i  [0; ),
Y
28  K.B. Erickson, R.E. Ladner, A. LaMarca
If U denotes the renewal measure, then the assumption (a) implies that U has an
absolutely continuous part which possesses a strictly positive density on (0; 1).
But, see [Feller 1971], page 369,
Consequently, the measure I 7! H k- x fIg also has an absolutely continuous part
which is strictly positive on [0; ). The conclusion one may draw from the preceding
is that f b
X(S k )g is irreducible with respect to the measure ' N , the Lebesgue measure
in R N (restricted to (0; ) N ). See [Revuz 1984], ch. 3 x2. This implies that the
trace chain f b
X(T k )g is also ' N irreducible on its state space A 0 . Together, these
two assertions imply that the full chain b
But we can also draw additional useful conclusions from (34) and (35).
From Stone's decomposition Theorem, [Revuz 1984], ch 5, x5, we can write
is a nite measure and U 1 is absolutely continuous with a bounded
continuous density u such that lim x!1 1=. For any Borel I  [0; 1) and
Hence, by dominated convergence,
lim
Z
I
[1 F (x)] dx:
This and the product formula (34) yield that for any Borel set A  [0; 1) N
lim
0 is the product measure F 0  F 0      F 0 and F 0 is the probability
distribution on [0; ) with density f1 F (x)g=. Not only does (36) give us one of
the limit Theorems we have used earlier, but it implies that subchains the b
and b
are both Harris recurrent (with invariant probabilities m
translate by - b
implies that, with
probability 1, b
for innitely many times k whatever be the initial
position, but b
is obtained from b
adding - to each component so the
previous assertion is also correct for b
with respect to its invariant probability.
See [Revuz 1984], ch 2, x3.
Consider now the full chain b
g. If ' N (A) > 0, A  , then
the preceding makes it clear that b
will hit A with positive probability. This implies
that b
X is ' N -irreducible. According to [Revuz 1984], ch 2, Theorem 2.3, 2.5, and
Denition 2.6, either b
X is a Harris chain with a (unique up to constant multiples)
invariant measure m, or else the potential kernel is proper. The potential kernel, K,
is dened by K(b x;
Ag. If K is proper, then  can be written
as an increasing sequence of subsets Dn each of which has bounded potential. But,
eventually any such sets must have positive Lebesgue measure, and in that case
Optimizing Static Calendar Queues  29
(34) implies that P b x f b
for innitely many positions
x. But K(b x; Dn ) < 1 implies that the expected total number of hits in Dn is nite
which implies that the number of hits must be nite with probability 1. We thus
cannot have a proper potential kernel and therefore b
X must be a Harris chain.
Consider next the aperiodicity issue. Seeking a contradiction, let us suppose that
X is periodic. Let fC i g d
be the recurrent cyclic classes in the decomposition of
the state space. These subsets have positive Lebesgue measure. Without loss of
generality we may suppose that mfC 1 \AN g > 0,  to be the
smallest S k such that b
This stopping time is nite on account of
.
s be a doubly indexed sequence of independent random variables each with
distribution F . The earliest possible epoch after  at which b
can arrive in A
no more than one particle can move at any
particular step. For any integer r  1 and any Borel rectangle
A 0 we have
Y
On account of hypotheses J1, the distribution F and each of its convolutions F r
puts positive mass on every subinterval in (0; ). Hence the right-hand side of
(37) is strictly positive whenever the cylinder set A has positive Lebesgue measure.
Standard measure theory implies that this is also correct for any Borel set A
[-) of positive measure.
probability 1, at time , b
X() belongs to C 1 \ AN  C 1 . So, if the chain
is periodic, then, w.p.1, at all future epochs of the form nd the chain will
always be found in C 1 . Also, for each d 1, at times t
mod d; the chain must belong to the set C 1+k(d) which is disjoint from the other
classes including C 1 . But in (37), for any N + r  N the right hand side is
strictly positive. By choosing A in (37) to be any set in C k \ [-) N of positive
measure and letting r take on dierent ( mod d) values, we get a contradiction to
the previous assertion about belonging to disjoint sets. There is no contradiction if
X is aperiodic.
It remains to show that the invariant measures m (unique up to constant multi-
ples) of the full chain b
are nite: m() < 1.
The trace chain f b
X(TK )g is positive recurrent with invariant probability m 0 . But
m 0 is also a multiple of m restricted to A 0 . Hence, It follows from the
Renewal Theorem that the mean return time to A 0 is also nite. Let   0 be a
bounded measurable function on A 0 . Then  is m 0 -(and hence m-) summable and
lim
Z
(b x)dm 0 (b x) > 0:
denote the number of visits to A 0 by the full chain during 0  s  t.
A. LaMarca
ih
But if the chain b
were null, that is if m had innite mass, then for any bounded
m-summable function , ([Revuz 1984], Theorem 2.6, page 198),
0:
By Fatou's Lemma, one can readily see that this would contradict (38).

ACKNOWLEDGMENTS

We would like to thank the two referees and the editor for their many valuable
suggestions for improving the paper.



--R

Implementation and analysis of binomial queue algorithms.
Calendar queues: A fast o(1) priority queue implementation for the simulation event set problem.
Calendar queue expectations.
An Introduction to Probability Theory and Its Applications
Description and analysis of an efcient
The Art of Computer Programming

Markov Chains.
Scalable architectures for integrated tra-c shaping and link scheduling in high-speed atm switches


A data structure for manipulating priority queues.
--TR
Self-adjusting binary search trees
adjusting heaps
Calendar queues: a fast 0(1) priority queue implementation for the simulation event set problem
A data structure for manipulating priority queues
The Art of Computer Programming, 2nd Ed. (Addison-Wesley Series in Computer Science and Information

--CTR
Wai Teng Tang , Rick Siow Mong Goh , Ian Li-Jin Thng, Ladder queue: An O(1) priority queue structure for large-scale discrete event simulation, ACM Transactions on Modeling and Computer Simulation (TOMACS), v.15 n.3, p.175-204, July 2005
Rick Siow Mong Goh , Ian Li-Jin Thng, Twol-amalgamated priority queues, Journal of Experimental Algorithmics (JEA), v.9 n.es, 2004
Farokh Jamalyaria , Rori Rohlfs , Russell Schwartz, Queue-based method for efficient simulation of biological self-assembly systems, Journal of Computational Physics, v.204 n.1, p.100-120, 20 March 2005
priority queue;optimization;calendar queue;discrete event simulation;data structures;algorithm analysis;markov chain