Time-Space Trade-offs for Triangulations and Voronoi Diagrams
Matias Korman1 , Wolfgang Mulzer∗,2 , André van Renssen3,4 , Marcel Roeloffzen3,4 , Paul
Seiferth∗,2 , and Yannik Stein∗,2

arXiv:1507.03403v2 [cs.CG] 8 Nov 2016

1

Tohoku University, Sendai, Japan., mati@dais.is.tohoku.ac.jp
2
Institut für Informatik, Freie Universität Berlin, Germany.
{mulzer,pseiferth,yannikstein}@inf.fu-berlin.de
3
National Institute of Informatics (NII), Tokyo, Japan., {andre,marcel}@nii.ac.jp
4
JST, ERATO, Kawarabayashi Large Graph Project.

Abstract
Let S be a planar n-point set. A triangulation for S is a maximal plane straight-line graph
with vertex set S. The Voronoi diagram for S is the subdivision of the plane into cells such
that all points in a cell have the same nearest neighbor in S. Classically, both structures can
be computed in O(n log n) time and O(n) space. We study the situation when the available
workspace is limited: given a parameter s ∈ {1, . . . , n}, an s-workspace algorithm has readonly access to an input array with the points from S in arbitrary order, and it may use only
O(s) additional words of Θ(log n) bits for reading and writing intermediate data. The output
should then be written to a write-only structure. We describe a deterministic s-workspace
algorithm for computing an arbitrary triangulation of S in time O(n2 /s + n log n log s) and
a randomized s-workspace algorithm for finding the Voronoi diagram of S in expected time
O((n2 /s) log s + n log s log∗ s).

1

Introduction

Since the early days of computer science, a major concern has been to cope with strong memory
constraints. This started in the ’70s [22] when memory was expensive. Nowadays, a major motivation comes from a proliferation of small embedded devices where large memory is neither feasible
nor desirable (e.g., due to constraints on budget, power, size, or simply to make the device less
attractive to thieves).
Even when memory size is not an issue, we might want to limit the number of write operations: one can read flash memory quickly, but writing (or even reordering) data is slow and may
reduce the lifetime of the storage system; write-access to removable memory may be limited for
technical or security reasons (e.g., when using read-only media such as DVDs or to prevent leaking
information about the algorithm). Similar problems occur when concurrent algorithms access data
simultaneously. A natural way to address this is to consider algorithms that do not modify the
input.
∗

WS and PS were supported in part by DFG Grants MU 3501/1 and MU 3501/2. YS was supported by the DFG
within the research training group “Methods for Discrete Structures” (GRK 1408).

1

The exact setting may vary, but there is a common theme: the input resides in read-only memory,
the output must be written to a write-only structure, and we can use O(s) additional variables to
find the solution (for a parameter s). The goal is to design algorithms whose running time decreases
as s grows, giving a time-space trade-off [23]. One of the first problems considered in this model
is sorting [19, 20]. Here, the time-space product is known to be Ω(n2 ) [8], and matching upper
bounds for the case b ∈ Ω(log n) ∩ O(n/ log n) were obtained by Pagter and Rauhe [21] (b denotes
the available workspace in bits).
Our current notion of memory constrained algorithms was introduced to computational geometry by Asano et al. [4], who showed how to compute many classic geometric structures with
O(1) workspace (related models were studied before [9]). Later, time-space trade-offs were given
for problems on simple polygons, e.g., shortest paths [1], visibility [6], or the convex hull of the
vertices [5].
We consider a model in which the set S of n points is in an array such that random access to
each input point is possible, but we may not change or even reorder the input. Additionally, we have
O(s) variables (for a parameter s ∈ {1, . . . , n}). We assume that each variable or pointer contains a
data word of Θ(log n) bits. Other than this, the model allows the usual word RAM operations. In
this setting we study two problems: computing an arbitrary triangulation for S and computing the
Voronoi diagram VD(S) for S. Since the output cannot be stored explicitly, the goal is to report
the edges of the triangulation or the vertices of VD(S) successively, in no particular order. Dually,
the latter goal may be phrased in terms of Delaunay triangulations. We focus on Voronoi diagrams,
as they lead to a more natural presentation.
Both problems can be solved in O(n2 ) time with O(1) workspace [4] or in O(n log n) time with
O(n) workspace [7]. However, to the best of our knowledge, no trade-offs were known before. Our
triangulation algorithm achieves a running time of O(n2 /s + n log n log s) using O(s) variables.
A key ingredient is the recent time-space trade-off by Asano and Kirkpatrick for triangulating a
special type of simple polygons [3]. This also lets us obtain significantly better running times for
the case that the input is sorted in x-order; see Section 2. For Voronoi diagrams, we use random
sampling to find the result in expected time O((n2 log s)/s + n log s log∗ s)); see Section 3. Together
with recent work of Har-Peled [16], this appears to be one of the first uses of random sampling
to obtain space-time trade-offs for geometric algorithms. The sorting lower bounds also apply to
triangulations and Voronoi diagrams (since we can reduce the former to the latter). This implies
that our second algorithm is almost optimal.

2

A Time-Space Trade-Off for General Triangulations

In this section we describe an algorithm that outputs the edges of a triangulation for a given point
set S in arbitrary order. For ease in the presentation we first assume that S is presented in sorted
order. In this case, a time-space trade-off follows quite readily from known results. We then show
how to generalize this for arbitrary inputs, which requires a careful adaptation of the existing data
structures.

2.1

Sorted Input

Suppose the input points S = {q1 , . . . , qn } are stored in increasing order of x-coordinate and that
all x-coordinates are distinct, i.e., xi < xi+1 for 1 ≤ i < n, where xi denotes the x-coordinate of qi .
2

A crucial ingredient in our algorithm is a recent result by Asano and Kirkpatrick for triangulating
monotone mountains 1 (or mountains for short). A mountain is a simple polygon with vertex
sequence v1 , v2 , . . . , vk such that the x-coordinates of the vertices increase monotonically. The edge
v1 vk is called the base. Mountains can be triangulated very efficiently with bounded workspace.
Theorem 2.1 (Lemma 3 in [3], rephrased). Let H be a mountain with n vertices, stored in sorted
x-order in read-only memory. Let s ∈ {2, . . . , n}. We can report the edges of a triangulation of H
in O(n logs n) time and using O(s) words of space.
Since S is given in x-order, the edges qi qi+1 , for 1 ≤ i < n, form a monotone simple polygonal
chain. Let Part(S) be the subdivision obtained by the union of this chain with the edges of the
convex hull of S (denoted by conv(S)). A convex hull edge is long if the difference between its
indices is at least two (i.e., the endpoints are not consecutive). The following lemma (illustrated in
Fig. 1) lets us decompose the problem into smaller pieces.

Figure 1: Any face of Part(S) is a mountain that is uniquely associated with a long convex hull
edge.
Lemma 2.2. Any bounded face of Part(S) is a mountain whose base is a long convex hull edge.
Moreover, no point of S lies in more than four faces of Part(S).
Proof. Any point qi ∈ S has at most four neighbors in Part(S): qi−1 , qi+1 , its predecessor and its
successor along the convex hull (if qi lies on conv(S)). Thus, no point of S belongs to more than
four faces of Part(S).
Next we show that every face F of Part(S) is a mountain with a long convex-hull edge as its
base. The boundary of F contains at least one long convex-hull edge e = (qi , qj ) (i < j), as other
edges connect only consecutive vertices. Since the monotone path qi , . . . , qj forms a cycle with the
edge e and since the boundary of F is a simple polygon, we conclude that e is the only long convexhull edge bounding F . Recall that e is a convex hull edge, and thus all points qi+1 , . . . , qj−1 lie on
one side of e and form a monotone chain (and in particular F is a mountain with base e).
The algorithm for sorted input is now very simple. We compute the edges of the convex hull
(starting from the leftmost point and proceeding in clockwise order). Whenever a long edge would
be reported, we pause the convex hull algorithm, and we triangulate the corresponding mountain.
Once the mountain has been triangulated, we resume with the convex hull algorithm until all convex
1

Also known as unimonotone polygons [15].

3

hull edges have been computed. The trade-off now follows from already existing trade-offs in the
various subroutines.
Theorem 2.3. Let S be a set of n points, sorted in x-order. We can report the edges of a triangulation of S in O(n2 ) time using O(1) variables, in O(n2 log n/2s ) time using O(s) variables (for any
s ∈ Ω(log log n) ∩ o(log n)), and in O(n logp n) time using O(p logp n) variables (for any 2 ≤ p ≤ n).
Proof. Correctness follows from Lemma 2.2, so we focus on the performance analysis. The main
steps are: (i) computing the convex hull of a point set given in x-order; and (ii) triangulating a
mountain.
By Theorem 2.1, we can triangulate a mountain Fi with ni vertices in time O(ni logs ni ) with
O(s) variables. We do not need to store Fi explicitly, since its vertices constitute a consecutive
subsequence of S and can be specified by the two endpoints of the base. No vertex appears in
more
than four mountains by Lemma 2.2, so the total time for triangulating the mountains is
P
O(n
i logs ni ) = O(n logs n). By reusing space, we can ensure that the total space requirement is
i
O(s).
Now we bound the time for computing conv(S). This algorithm is paused to triangulate mountains, but overall it is executed only once. There are several convex hull algorithms for sorted point
sets under memory constraints. If s ∈ Θ(1), we can use gift-wrapping (Jarvis march [17]), which
runs in O(n2 ) time. Barba et al. [5] provided a different algorithm that runs in O(n2 log n/2s ) time
using O(s) variables (for any s ∈ o(log n)).2 This approach is desirable for s ∈ Ω(log log n)∩o(log n).
As soon as s = Ω(log n), we can use the approach of Chan and Chen [10]. This algorithm runs in
O(n logp n) time and uses O(p logp n) variables, for any 2 ≤ p ≤ n. Regardless of the size of the
workspace, the time for computing the convex hull dominates the time needed for triangulating all
mountains.
A similar approach is unlikely to work for the Delaunay triangulation, since knowing the x-order
of the input does not help in computing it [14].

2.2

General Input

The algorithm from Section 2.1 uses the sorted order in two ways. Firstly, the convex-hull algorithms
of Barba et al. [5] and of Chan and Chen [10] work only for simple polygons (e.g., for sorted input).
Instead, we use the algorithm by Darwish and Elmasry [13] that gives the upper (or lower) convex
hull of any sequence of n points in O(n2 /(s log n) + n log n) time with O(s) variables3 , matching
known lower bounds. Secondly, and more importantly, the Asano-Kirkpatrick (AK) algorithm for
triangulating a mountain requires the input to be sorted. To address this issue, we simulate sorted
input using multiple heap structures. This requires a close examination of how the AK-algorithm
accesses its input.
Let F be a mountain with n vertices. Let F ↑ and F ↓ denote the vertices of F in ascending and
in descending x-order. The AK-algorithm has two phases, one focused on F ↑ and the other one on
2

In fact, Barba et al. show how to compute the convex hull of a simple polygon, but also show that both problems
are equivalent. The monotone chain can be completed to a polygon by adding a vertex with a very high or low
y-coordinate.
3
Darwish and Elmasry [13] state a running time of O(n2 /s + n log n), but they measure workspace in bits, while
we use words.

4

F ↓ .4 Each pass computes a portion of the triangulation edges, uses O(s) variables, and scans the
input Θ(logs n) times. We focus on the approach for F ↑ .
As mentioned, the algorithm uses Θ(logs n) rounds. In round i, it partitions F into blocks of
O(|F |/si ) consecutive points that are processed from left to right. Each block is further subdivided
into O(s) sub-blocks b1 , . . . , bk of size O(|F |/si+1 ). The algorithm does two scans over the subblocks. The first scan processes the elements from left to right. Whenever the first scan finishes
reading a sub-block bi , the algorithm makes bi active and creates a pointer li to the rightmost
element of bi . The second scan goes from right to left and is concurrent to the first scan. In each
step, it reads the element at li in the rightmost active sub-block bi , and it decreases li by one. If li
leaves bi , then bi becomes inactive. As the first scan creates new active sub-blocks as it proceeds,
the second scan may jump between sub-blocks. The interested reader may find a more detailed
description in A.
To provide the input for the AK-algorithm, we need the heap by Asano et al. [2]. For completeness, we briefly restate its properties here.
Lemma 2.4 ([2]). Let S be a set of n points. There is a heap that supports insert and extract-min
(resp. extract-max) in O (n/(s log n) + log s)D(n)) time using O(s) variables, where D(n) is the
time to decide whether a given element currently resides in the heap (is alive).5
Proof. We first describe the data structure. Then we discuss how to perform insertions and extractmin operations.
We partition the input into s log n consecutive buckets of equal size, and we build a complete
binary tree T over the buckets. Let v be a node of T with height h. Then, there are 2h buckets
below v in T . We store 2h information bits in v to specify the minimum alive element below v.
The first h bits identify the bucket containing the minimum. We further divide this bucket into
2h consecutive parts of equal size, called quantiles. The second h bits in v specify the quantile
containing the minimum. If 2h > log n, we use log n bits to specify the minimum directly. Hence,
the total number of bits is bounded by
log(s log n)

X
h=0

s log n
min{2h, log n} = O(s log n).
2h

Therefore we need O(s) variables in total.
Let v be a node with height h. To find the minimum alive element in T below v, we use the 2h
information bits stored in v. First, we identify the bucket containing
the minimum and the correct

quantile within this bucket. This quantile contains O 2h s nlog n elements. For each element in the
quantile, we decide inD(n) time whether it is alive, and we return the minimum such element. This
takes O 2h s nlog n D(n) time in total.
insert: Assume we want to insert an element x that is at position i in the input array. Let v be the
parent of the leaf of T corresponding to the bucket that contains x. We update the information
bits at each node u on the root path starting at v. To do so, we use the information bits in u to
4

AK reduce triangulation to the next smaller right neighbor (NSR) and the next smaller left neighbor (NSL)
problem and present an algorithm for NSR if the input is in x-order. This implies an NSL-algorithm by reading the
input in reverse.
5
The bounds in [2] do not include the factor D(n) since the authors studied a setting similar to Lemma 2.5 where
it takes O(1) time to decide whether an element is alive.

5

find the minimum element in the buckets covered by u, as described above. Then we compare
it with x. If x is larger, we are done and we stop the insertion. Otherwise, we update the
information bits at u to the bucket and quantile that contain x. If we reach and update the
root node, we also update the pointer that points to the minimum element in the heap. The

work per node is dominated by the costs for finding the minimum, which is O 2h s nlog n D(n) .
Thus, the total cost for insertion is bounded by
log(s log n)

X
h=0


 n
n
D(n)
.
D(n)
=
O
s log n
2h s log n

extract-min: First we use the pointer to the minimum alive element to determine the element x
to return. Then we use a similar update strategy as for insertions. Let v be the leaf node
corresponding to the bucket of x. We first update the information bits of v by scanning through
the whole bucket of v and determining the smallest alive element. Since a bucket contains
O(n/s log n) elements, this needs time O(n/(s log n)D(n)). Then we update the information
bits of each node u on the path for v as follows: let v1 and v2 be the two children of u. We
determine the minimum alive element in the buckets covered by v1 and v2 , take the smaller
one, and use it to update the information bits at u. Once we reach the root, we also update
the pointer to the minimum element of the heapto the new minimum element of the root.
n
The total time again is bounded by O s log
n D(n) .

Lemma 2.5 ([2]). Let S be a set of n points. We can build a heap with all elements in S in O(n)
time that supports extract-min in O n/(s log n) + log n) time using O(s) variables.
Proof. The construction time is given in [2]. To decide in O(1) time if some x ∈ S is alive, we store
the last extracted minimum m and test whether x > m.
We now present the complete algorithm. We show how to subdivide S into mountains Fi and
how to run the AK-algorithm on each Fi↑ . By reversing the order, the same discussion applies to Fi↓ .
Sorted input is emulated by two heaps H1 , H2 for S according to x-order. By Lemma 2.5, each heap
uses O(s) space, can be constructed in O(n) time, and supports extract-min in O(n/(s log n)+log n)
worst-case time. We will use H1 to determine the size of the next mountain Fi and H2 to process
the points of Fi .
We execute the convex hull algorithm with Θ(s) space until it reports the next convex hull edge
pq. Throughout the execution of the algorithm, heaps H1 and H2 contain exactly the points to the
right of p. We repeatedly extract the minimum of H1 until q becomes the minimum element. Let
k be the number of removed points.
If k = 1, then pq is short. We extract the minimum of H2 , and we continue with the convex
hull algorithm. If k ≥ 2, then Lemma 2.2 shows that pq is the base of a mountain F that consists
of all points between p and q. These are exactly the k + 1 smallest elements in H2 (including p
and q). If k ≤ s, we extract them from H2 , and we triangulate F in memory. If k > s, we execute
the AK-algorithm on F using O(s) variables. At the beginning of the ith round, we create a copy
H(i) of H2 , i.e., we duplicate the O(s) variables that determine the state of H2 . Further, we create
an empty max-heap H(ii) using O(s) variables to provide input for the second scan. To be able
6

0 of H . Whenever the AK-algorithm requests
to reread a sub-block, we create a further copy H(i)
2
the next point in the first scan, we simply extract the minimum of H(i) . When a sub-block is fully
0 to reread the elements and insert them into H
read, we use H(i)
(ii) . Now, the rightmost element of
all active sub-blocks corresponds exactly to the maximum of H(ii) . One step in the second scan is
equivalent to an extract-max on H(ii) .
0 , and H
At the end of a round, we delete H(i) , H(i)
(ii) , so that the space can be reused in the
next round. Once the AK-algorithm finishes, we repeatedly extract the minimum of H2 until we
reach q.

Theorem 2.6. We can report the edges of a triangulation of a set S of n points in time O(n2 /s +
n log n log s) using O(s) additional variables.
Proof. Similarly as before, correctness directly follows from Lemma 2.2 and the correctness of the
AK-algorithm. The bound on the space usage is immediate.
Computing the convex hull now needs O(n2 /(s log n) + n log n) time [13]. By Lemma 2.5, the
heaps H1 and H2 can be constructed in O(n) time. During execution, we perform n extract-min
operations on each heap, requiring O(n2 /(s log n) + n log n) time in total.
Let Fj be a mountain with nj vertices that is discovered by the convex hull algorithm. If
nj ≤ s, then Fj is triangulated in memory in O(nj ) time, and the total time for such mountains
is O(n). If nj > s, then the AK-algorithm runs in O(nj logs nj ) time. We must also account
for providing the input for the algorithm. For this, consider some round i ≥ 1. We copy H2 to
H(i) in O(s) time. This time can be charged to the first scan, since nj > s. Furthermore, we
perform nj extract-min operations on H(i) . Hence the total time to provide input for the first scan
is O(nj n/(s log n) + nj log n).
0 of H . Again, the time for this can be charged to
For the second scan, we create another copy H(i)
2
0 which takes O(n n/(s log n)+n log n)
the scan. Also, we perform nj extract-min operations on H(i)
j
j
time. Additionally, we insert each fully-read block into H(ii) . The main problem is to determine if an
element in H(ii) is alive: there are at most O(s) active sub-blocks. For each active sub-block bi , we
know the first element yi and the element zi that li points to. An element is alive if and only if it is
in the interval [yi , zi ] for some active bi . This can be checked in O(log s) time. Thus, by Lemma 2.4,
each insert and extract-max on H(ii) takes O (n/(s log n) + log s) log s) time. Since each element is
inserted once, the total time to provide input to the second scan is O(nj log(s)(n/(s log n) + log s)).
This dominates the time for the first scan. There areO(logs nj ) rounds, so we can triangulate
Fj in time O nj logs nj + nj log(nj ) n/(s log n) + log s . Summing over all Fj , the total time is
O(n2 /s + n log n log s).

3

Voronoi Diagrams

Given a planar n-point set S, we would like to find the vertices of VD(S). Let K = {p1 , p2 , p3 } be a
triangle with S ∩ K = ∅, S ⊆ conv(K), and so that all vertices of VD(S) are vertices of VD(S ∪ K).
For example, we can set K = {(−κ, −κ), (−κ, κ), (0, κ)} for some large κ > 0. Since the desired
properties hold for all large enough κ, we do not need to find an explicit value for it. Instead,
whenever we want to evaluate a predicate involving points from K, we can take the result obtained
for κ → ∞.
Our algorithm relies on random sampling. First, we show how to take a random sample from
S with small workspace. One of many possible approaches is the following one that ensures a
7

worst-case guarantee:
Lemma 3.1. We can sample a uniform random subset R ⊆ S of size s in time O(n + s log s) and
space O(s).
Proof. The sampling algorithm consists of two phases. In the first phase, we sample a random
sequence I of s distinct numbers from [n].6 The phase proceeds in s rounds. At the beginning of
round k, for k = 1, . . . , s, we have already sampled a sequence I of k − 1 numbers from [n], and we
would like to pick an element from [n] \ I uniformly at random. We store I in a binary search tree
T . We maintain the invariant that T stores with each element x ∈ [n − k + 1] ∩ I a replacement
ρx ∈ {n − k + 2, . . . , n} \ I such that [n] \ I = ([n − k + 1] \ I) ∪ {ρx | x ∈ [n − k + 1] ∩ I}, see
Figure 2. In round k, we sample a random number x from [n − k + 1], and we check in T whether
n−k+1

1

n

Figure 2: Sampling a random sequence I from [n]. At the beginning of round k, we have already
sampled k − 1 elements (shown in gray). Each element x ∈ [n − k + 1] ∩ I has a replacement
ρx ∈ {n − k + 2, . . . , n} \ I (indicated by the arrows). In round k, we pick a random number
x ∈ [n − k + 1]. If x is already contained in I, we add ρx to I. Otherwise, we add x.
x ∈ I. If not, we add x to I (and T ), otherwise, we add ρx to I (and T ). By the invariant, we add
a uniform random element from [n] \ I to I.
It remains to update the replacements, see Figure 3. If x = n − k + 1, we do not need a
replacement for x. Now suppose x < n − k + 1. If n − k + 1 6∈ I, we set ρx = n − k + 1. Otherwise,
we set ρx = ρn−k+1 . This ensures that the invariant holds at the beginning of round k + 1, and it
takes O(log s) time and O(s) space. We continue for s rounds. At the end of the first phase, any
sequence of s distinct numbers in [n] is sampled with equal probability. Furthermore, the phase
takes O(s log s) time and O(s) space.
In the second phase, we scan through S to obtain the elements whose positions correspond to
the numbers in I. This requires O(n) time and O(s) space.
We use Lemma 3.1 to find a random sample R ⊆ S of size s. We compute VD(R∪K), triangulate
the bounded cells and construct a planar point location structure for the triangulation. This takes
O(s log s) time and O(s) space [18]. By our choice of K, all Voronoi cells for points in R are
bounded, and every point in S lies in a bounded Voronoi cell. Given a vertex v ∈ VD(R ∪ K), the
conflict circle of v is the largest circle with center v and no point from R ∪ K in its interior. The
conflict set Bv of v contains all points from S that lie in the conflict circle of v, and the conflict
size bv of v is |Bv |. We scan through S to find the conflict size bv for each vertex v ∈ VD(R ∪ K):
every Voronoi vertex has a counter that is initially 0. For each p ∈ S \ (R ∪ K), we use the point
location structure to find the triangle ∆ of VD(R ∪ K) that contains it. At least one vertex v of
∆ is in conflict with p. Starting from v, we walk along the edges of VD(R ∪ K) to find all Voronoi
vertices in conflict with p (recall that these vertices induce a connected component in VD(R ∪ K)).
We increment the counters of all these vertices. This may take a long time in the worst case, so we
6

We write [n] for the set {1, . . . , n}.

8

1

n−k+1

n

1

n−k+1

n

1

n−k

n

1

n−k

n

x

x

1

n−k+1

n

1

n−k

n

x

Figure 3: Finding a replacement for x. If x = n − k + 1, we do not need a replacement for x in the
next round (top left). If n − k + 1 is not sampled yet, we can make it the replacement for x (top
right). Otherwise, we make the old replacement for n − k + 1 the new replacement for x (bottom).
impose an upper bound on the total work. For this, we choose a threshold M . When the sum of the
conflict counters exceeds M , we start over with a new sample R. The total time for one attempt
is O(n log s + M ), and below we prove that for M = Θ(n), the success probability is at least 3/4.
Next, we pick another threshold T , and we compute for each vertex v of VD(R ∪ K) the excess
tv = bv s/n.
P The excess measures how far the vertex deviates from the desired conflict size n/s. We
check if v∈VD(R∪K) tv log tv ≤ T . If not, we start over with a new sample. Below, we prove that
for T = Θ(s), the success probability is at least 3/4. The total success probability is 1/2, and the
expected number
in expected time O(n log s + s log s), we can find a sample
P of attempts is 2. Thus, P
R ⊆ S with v∈VD(R∪K) bv = O(n) and v∈VD(R∪K) tv log tv = O(s).
We now analyze the success probabilities, using the classic Clarkson-Shor method [12]. We begin
with a variant of the Chazelle-Friedman bound [11].
Lemma 3.2. Let X be a planar point set of size m, and let Y ⊂ R2 with |Y | ≤ 3 and X ∩ Y = ∅.
For fixed p ∈ (0, 1], let R ⊆ X be a random subset of size pm and let R0 ⊆ X be a random subset
of size p0 m, for p0 = p/2. Suppose that p0 m ≥ 4. Fix u ⊂ X ∪ Y with |u| = 3, and let vu be the
Voronoi vertex defined by u. Let bu be the number of points from X ∪ Y in the interior of the circle
with center vu and with the points from u on the boundary. Then,
Pr[vu ∈ VD(R ∪ Y )] ≤ 64e−pbu /2 Pr[vu ∈ VD(R0 ∪ Y )].
Proof. Let σ = Pr[vu ∈ VD(R ∪ Y )] and σ 0 = Pr[vu ∈ VD(R0 ∪ Y )]. The vertex vu is in VD(R ∪ Y )
precisely if u ⊆ R ∪ Y and Bu ∩ (R ∪ Y ) = ∅, where Bu are the points from X ∪ Y inside the circle
with center vu and with the points from u on the boundary. If Bu ∩ Y 6= ∅, then σ = σ 0 = 0, and
the lemma holds. Thus,
assume that Bu ⊆ X. Let du = |u ∩ X|, the number of points in u from X.

m−bu −du
There are pm−du ways to choose a pm-subset from X that avoids all points in Bu and contains
all points of u ∩ X, so
9

  
m
m − bu − du
σ=
pm − du
pm
, Qpm−1
Qpm−du −1
(m − bu − du − j)
j=0
j=0 (m − j)
= Qpm−du −1
Qpm−1
(pm − du − j)
j=0 (pm − j)
j=0


=

dY
u −1

pm−d
Yu −1

j=0

pm − j
·
m−j

du

pm−d
Yu −1 

≤p

j=0

m − bu − du − j
m − du − j

bu
1−
m − du − j

j=0


.

Similarly, we get
0

σ =

dY
u −1
i=0

p0 m − i
m−i

p0 m−d
Yu −1 
j=0

bu
1−
m − du − j


,

and since p0 m ≥ 4 and i ≤ 2, it follows that
 0 du
p
σ ≥
2
0

p0 m−d
Yu −1 
j=0

1−

bu
m − du − j


.

Therefore, since p0 = p/2,
σ
≤
σ0



2p
p0

du

pm−d
Yu −1
j=p0 m−du



bu
1−
m − du − j





bu pm/2
≤ 64 1 −
≤ 64 e−pbu /2 .
m

We can now bound the total expected conflict size.
hP
i
Lemma 3.3. We have E
b
= O(n).
v
v∈VD(R∪K)
Proof. By expanding the expectation, we get


X
X
E
bv  =
Pr[vu ∈ VD(R ∪ K)]bu ,
v∈VD(R∪K)

u⊂S∪K,|u|=3

with vu being the Voronoi vertex of u and bu its conflict size. By Lemma 3.2 with X = S, Y = K
and p = s/n, this is
≤

X
u⊂S∪K,|u|=3

64e−pbu /2 Pr[vu ∈ VD(R0 ∪ K)]bu ,

10

where R0 ⊆ S is a sample of size s/2. We bound this as
≤

≤

∞
X

X

i=0 u⊂S∪K,|u|=3
)
bu ∈[ pi , i+1
p

1
p

X
u⊂S∪K,|u|=3

64e−i/2 (i + 1)
Pr[vu ∈ VD(R0 ∪ K)]
p

Pr[vu ∈ VD(R0 ∪ K)]

∞
X

64e−i/2 (i + 1)

i=0

= O(s/p) = O(n),
since
O(1).

P

u⊂S∪K,|u|=3 Pr[vu

∈ VD(R0 ∪ K)] = O(s) is the size of VD(R0 ∪ K) and

P∞

i=0 e

−i/2 (i + 1)

=

By P
Lemma 3.3 and Markov’s inequality, we can conclude that there is an M = Θ(n) such that
Pr[ v∈VD(R∪K) bv > M ] ≤ 1/4.
hP
i
Lemma 3.4. E
v∈VD(R∪K) tv log tv = O(s).
Proof. By Lemma 3.2 with X = S, Y = K, and p = s/n,


X
X
E
tv log tv  =
Pr[vu ∈ VD(R ∪ K)] tu log tu
v∈VD(R∪K)

u⊂S∪K,|u|=3

≤
≤

≤

X
u⊂S∪K,|u|=3
∞
X

64e−pbu /2 Pr[vu ∈ VD(R0 ∪ K)]tu log tu

X

i=0 u⊂S∪K,|u|=3
)
bu ∈[ pi , i+1
p
∞
X

i

64e− 2 (i + 1)2

i=0

i

64e− 2 (i + 1)2 Pr[vu ∈ VD(R0 ∪ K)]

X
u⊂S∪K,|u|=3

Pr[vu ∈ VD(R0 ∪ K)]

= O(s).

By Markov’s inequality and Lemma 3.4, we can conclude that there is a T = Θ(s) with
P
Pr[ v∈VD(R∪K) tv log tv ≥ T ] ≤ 1/4. This finishes the first sampling phase.
The next goal is to sample for each vertex v with tv ≥ 2 a random subset Rv ⊆ Bv of size
min{αtv log tv , bv } for large enough α > 0 (recall that Bv is the conflict set of v and that bv = |Bv |).
Lemma 3.5. In total time O(n log s), we can sample for each vertex v ∈ VD(R ∪ K) with tv ≥ 2 a
random subset Rv ⊆ Bv of size min{αtv log tv , bv }.
Proof. First, we sample for each vertex v with tv ≥ 2 a sequence Iv of αtv log tv distinct numbers
from {1, . . . , bv }. For this, we use the first phase of the algorithm from the proof of Lemma 3.1 for
11

each such vertex, but without reusing the space. As explained in the proof of Lemma 3.1, this takes
total time
!
!
X
X
O
(tv log tv ) log(tv log tv ) = O
(tv log tv ) log s = O(s log s),
v

v

P

since v tv log tv = O(s), and in particular tv log tv P
= O(s) for each vertex v (note that the constant
in the O-notation is independent of v). Also, since v tv log tv = O(s), the total space requirement
is O(s).
After that, we scan through S. For each vertex v, we have a counter cv , initialized to 0. For
each p ∈ S, we find the conflict vertices of p, and for each conflict vertex v, we increment cv . If cv
appears in the corresponding set Iv , we add p to Rv . The total running time is O(n log s), as we do
one point location for each input point and the total conflict size is O(n).
We next show that for a fixed vertex v ∈ VD(R ∪ K), with constant probability, all vertices in
VD(Rv ) have conflict size n/s with respect to Bv .
Lemma 3.6. Let v ∈ VD(R ∪ K) with tv ≥ 2, and let Rv ⊆ Bv be the sample for v. The expected
number of vertices v 0 in VD(Rv ) with at least n/s points from Bv in their conflict circle is at most
1/4.
Proof. If Rv = Bv , the lemma holds, so assume that αtv log tv < bv . Recall that tv = bv s/n. We
have
"
#
X
X
E
1 =
Pr[vu0 ∈ VD(Rv )],
v 0 ∈VD(Rv )
b0v0 ≥n/s

u⊂Bv ,|u|=3
b0u ≥n/s

where b0u denotes the number of points from Bv inside the circle with center vu0 and with the points
from u on the boundary. Using Lemma 3.2 with X = Bv , Y = ∅, and p = (αtv log tv )/bv =
α(s/n) log tv , this is
≤

X
u⊂Bv ,|u|=3
b0u ≥n/s

0

64e−pbu /2 Pr[vu0 ∈ VD(Rv0 )]

≤ 64e−(α/2) log tv
=

O(t−α/2
tv
v

X
u⊂Bv ,|u|=3

Pr[vu0 ∈ VD(Rv0 )]

log tv ) ≤ 1/4,

for α large enough (remember that tv ≥ 2).
By Lemma 3.6 and Markov’s inequality, the probability that all vertices from VD(Rv ) have at
most n/s points from Bv in their conflict circles is at least 3/4. If so, we call v good, see Figure 4.
Scanning through S, we can identify the good vertices in time O(n log s) and space O(s). Let s0 be
the size of VD(R ∪ K). If we have less than s0 /2 good vertices, we repeat the process. Since the
expected number of good vertices is 3s0 /4, the probability that there are at least s0 /2 good vertices
is at least 1/2, by Markov’s inequality. Thus, in expectation, we need to perform the sampling
12

Figure 4: A Voronoi Diagram of the sampled set R (left). The two red square vertices of VD(R ∪ K)
are not good and we need to resample within their conflict list (the blue crosses) and compute the
new Voronoi Diagram (right).
twice. For the remaining vertices, we repeat the process, but now we take two samples per vertex,
decreasing the failure probability to 1/4. We repeat the process, taking in each round the maximum
number of samples that fit into the work space. In general, if we have s0 /ai active vertices in round
i, we can take ai samples per vertex, resulting in a failure probability of 2−ai . Thus, the expected
number of active vertices in round i + 1 is s0 /ai+1 = s0 /(ai 2ai ). After O(log∗ s) rounds, all vertices
are good. To summarize:
Lemma 3.7. In total expected time O(n log s log∗ s) and space O(s),
Pwe can find sets R ⊆ S and
Rv ⊂ Bv for each vertex v ∈ VD(R ∪ K) such that (i) |R| = s: (ii) v∈VD(R∪K) |Rv | = O(s); and
(iii) for every Rv , all vertices of VD(Rv ) have at most n/s points from Bv in their conflict circle.
S
We set R2 = R ∪ v∈VD(R∪K) Rv . By Lemma 3.7, |R2 | = O(s). We compute VD(R2 ∪ K)
and triangulate its bounded cells. For a triangle ∆ of the triangulation, let
S3r ∈ R2 ∪ K be the site
whose cell contains ∆, and v1 , v2 , v3 the vertices of ∆. We set B∆ = {r} ∪ i=1 Bvi . Using the next
lemma, we show that |B∆ | = O(n/s).
Lemma 3.8. Let S ⊂ R2 and ∆ = {v1 , v2 , v3 } a triangle in the triangulation of VD(S). Let x ∈ ∆.
Then any circle C with center x that contains no points from S is covered by the conflict circles of
v1 , v2 and v3 .
Proof. Let p ∈ C and let r ∈ S be the site whose cell contains ∆. We show that p is contained in
the conflict circle of v1 , v2 , or v3 . Consider the bisector B of p and r. Since C contains p but not
r, we have d(x, p) < d(x, r), so x lies on the same side of B as p. Since x ∈ ∆, at least one of v1 ,
v2 , v3 , is on the same side of B as p; say v1 . This means that d(v1 , p) < d(v1 , r), so p lies inside the
circle around v1 with r on the boundary. This is precisely the conflict circle of v1 .
Lemma 3.9. Any triangle ∆ in the triangulation of VD(R2 ∪ K) has |B∆ | = O(n/s).
Proof. Let v be a vertex of ∆. We show that bv = O(n/s). Let ∆R = {v1 , v2 , v3 } be
S the triangle
in the triangulation of VD(R) that contains v. By Lemma 3.8, we have Bv ⊆ 3i=1 Bvi . We
consider the intersections Bv ∩ Bvi , for i = 1, 2, 3. If tvi < 2, then bvi = O(n/s) and |Bv ∩ Bvi | =
O(n/s). Otherwise, we have sampled a set Rvi for vi . Let ∆i = {w1 , w2 , w3 } beS
the triangle in the
triangulation of VD(Rvi ) that contains v. Again, by Lemma 3.8, we have Bv ⊆ 3j=1 Bwj and thus
S
also Bv ∩ Bvi ⊆ 3j=1 Bwj ∩ Bvi . However, by construction of Rvi , |Bwj ∩ Bvi | is at most n/s for
j = 1, 2, 3. Hence, |Bv ∩ Bvi | = O(n/s) and bv = O(n/s).
13

The following lemma enables us to compute the Voronoi diagram of R2 ∪ K locally for each
triangle ∆ in the triangulation of VD(R2 ∪ K) by only considering sites in B∆ . It is a direct
consequence of Lemma 3.8.
Lemma 3.10. For every triangle ∆ in the triangulation of VD(R2 ∪ K), we have VD(S ∪ K) ∩ ∆ =
VD(B∆ ) ∩ ∆.
Theorem 3.11. Let S be a planar n-point set. In expected time O((n2 /s) log s + n log s log∗ s) and
space O(s), we can compute all Voronoi vertices of S.
Proof. We compute a set R2 as above. This takes O(n log s log∗ s) time and space O(s). We
triangulate the bounded cells of VD(R2 ∪ K) and compute a point location structure for the result.
Since there are O(s) triangles, we can store the resulting triangulation in the workspace. Now, the
goal is to compute simultaneously for all triangles ∆ the Voronoi diagram VD(B∆ ) and to output
all Voronoi vertices that lie in ∆ and are defined by points from S. By Lemma 3.10, this gives all
Voronoi vertices of VD(S).
Given a planar m-point set X, the algorithm by Asano et al. finds all vertices of VD(X) in O(m)
scans over the input, with constant workspace [4]. We can perform a simultaneous scan for all sets
B∆ by determining for each point in S all sets B∆ that contain it. This takes total time O(n log s),
since we need one point location for each p ∈ S and since the total size of the B∆ ’s is O(n). We
need O(max∆ |B∆ |) = O(n/s) such scans, so the second part of the algorithm needs O((n2 /s) log s)
time.
As mentioned in the introduction, Theorem 3.11 also lets us report all edges of the Delaunay
triangulation of S in the same time bound: by duality, the three sites that define a vertex of VD(S)
also define a triangle for the Delaunay triangulation. Thus, whenever we discover a vertex of VD(S),
we can instead output the corresponding Delaunay edges, while using a consistent tie-breaking rule
to make sure that every edge is reported only once.
Acknowledgments This work began while W. Mulzer, P. Seiferth, and Y. Stein visited the
Tokuyama Laboratory at Tohoku University. We would like to thank Takeshi Tokuyama and all
members of the lab for their hospitality and for creating a conducive and stimulating research
environment.

References
[1] T. Asano, K. Buchin, M. Buchin, M. Korman, W. Mulzer, G. Rote, and A. Schulz. Memoryconstrained algorithms for simple polygons. Comput. Geom., 46(8):959–969, 2013.
[2] T. Asano, A. Elmasry, and J. Katajainen. Priority queues and sorting for read-only data. In
Theory and applications of models of computation (TAMC), pages 32–41, 2013.
[3] T. Asano and D. Kirkpatrick. Time-space tradeoffs for all-nearest-larger-neighbors problems.
In Proc. 13th Int. Symp. Algorithms and Data Structures (WADS), pages 61–72, 2013.
[4] T. Asano, W. Mulzer, G. Rote, and Y. Wang. Constant-work-space algorithms for geometric
problems. J. of Comput. Geom., 2(1):46–68, 2011.

14

[5] L. Barba, M. Korman, S. Langerman, K. Sadakane, and R. I. Silveira. Space–time trade-offs
for stack-based algorithms. Algorithmica, 72(4):1097–1129, 2015.
[6] L. Barba, M. Korman, S. Langerman, and R. I. Silveira. Computing the visibility polygon
using few variables. Comput. Geom., 47(9):918–926, 2013.
[7] M. de Berg, O. Cheong, M. van Kreveld, and M. Overmars. Computational geometry. Algorithms and applications. Springer-Verlag, third edition, 2008.
[8] A. Borodin and S. Cook. A time-space tradeoff for sorting on a general sequential model of
computation. SIAM J. Comput., 11:287–297, 1982.
[9] H. Brönnimann, T. M. Chan, and E. Y. Chen. Towards in-place geometric algorithms and data
structures. In Proc. 20th Annu. Sympos. Comput. Geom. (SoCG), pages 239–246, 2004.
[10] T. M. Chan and E. Y. Chen. Multi-pass geometric algorithms. Discrete Comput. Geom.,
37(1):79–102, 2007.
[11] B. Chazelle and J. Friedman. A deterministic view of random sampling and its use in geometry.
Combinatorica, 10(3):229–249, 1990.
[12] K. L. Clarkson and P. W. Shor. Applications of random sampling in computational geometry,
II. Discrete Comput. Geom., 4:387–421, 1989.
[13] O. Darwish and A. Elmasry. Optimal time-space tradeoff for the 2D convex-hull problem. In
Proc. 22nd Annu. European Sympos. Algorithms (ESA), pages 284–295, 2014.
[14] H. Djidjev and A. Lingas. On computing Voronoi diagrams for sorted point sets. Internat. J.
Comput. Geom. Appl., 5(3):327–337, 1995.
[15] A. Fournier and D. Y. Montuno. Triangulating simple polygons and equivalent problems. ACM
Transactions on Graphics, 3:153–174, 1984.
[16] S. Har-Peled. Shortest path in a polygon using sublinear space. J. of Comput. Geom., 7(2):19–
45, 2016.
[17] R. Jarvis. On the identification of the convex hull of a finite set of points in the plane. Inform.
Process. Lett., 2(1):18–21, 1973.
[18] D. Kirkpatrick. Optimal search in planar subdivisions. SIAM J. Comput., 12(1):28–35, 1983.
[19] J. I. Munro and M. Paterson. Selection and sorting with limited storage. Theoret. Comput.
Sci., 12:315–323, 1980.
[20] J. I. Munro and V. Raman. Selection from read-only memory and sorting with minimum data
movement. Theoret. Comput. Sci., 165(2):311–323, 1996.
[21] J. Pagter and T. Rauhe. Optimal time-space trade-offs for sorting. In Proc. 39th Annu. IEEE
Sympos. Found. Comput. Sci. (FOCS), pages 264–268, 1998.
[22] I. Pohl. A minimum storage algorithm for computing the median. Technical Report RC2701,
IBM, 1969.
15

[23] J. E. Savage. Models of computation—exploring the power of computing. Addison-Wesley, 1998.

16

A

The Asano-Kirkpatrick Algorithm

We give more details on the algorithm of Asano and Kirkpatrick [3]. Let F be a mountain with
vertices q1 , . . . , qn sorted in x-order and base q1 qn . We define the height h(qi ) of qi , i = 1, . . . , n, as
the distance from qi to the line through the base. Let A = (q1 , . . . , qn ) be the input array. A vertex
qr is the nearest-smaller-right-neighbor (NSR) of a vertex ql if (i) l < r; (ii) h(ql ) > h(qr ); and (iii)
h(ql ) ≤ h(qk ) for l < k < r. We call (ql , qr ) a NSR-pair, with left endpoint ql and right endpoint
qr . Nearest-smaller-left-neighbors (NSL) and NSL-pairs are defined similarly. Let R be the set of
all NSR-pairs and L be the set of all NSL pairs. Asano and Kirkpatrick show that the edges R ∪ L
triangulate F . We describe the algorithm for computing R. The algorithm for L is the same, but
it reads the input in reverse.
Let s denote the space parameter. The algorithm runs in logs n rounds. In round i, i =
0, . . . , logs n − 1, we partition A into si consecutive blocks of size n/si . Each block B is further
partitioned into s consecutive sub-blocks b1 , . . . , bs of size n/si+1 . In each round, we compute only
NSR-pairs with endpoints in different sub-blocks of the same block. We handle each block B
individually as follows. The sub-blocks of B are visited from left to right. When we visit a subblock bj , we compute all NSR-pairs with a right endpoint in bj and a left endpoint in the sub-blocks
b1 , . . . , bj−1 . Initially, we visit the first sub-block b1 and we push a pointer to the rightmost element
in b1 onto a stack S. We call a sub-block with a pointer in S active. Assume now that we have
already visited sub-blocks b1 , . . . , bj−1 . Let l be the topmost pointer in S, referring to an element ql
in bj 0 , j 0 < j. Furthermore, let r be a pointer to the leftmost element qr in bj . If h(ql ) > h(qr ), we
output (ql , qr ) and we decrement l until we find the first element whose height is smaller than the
current h(ql ). If l leaves bj 0 , this sub-block becomes inactive and we remove l from S. We continue
with the new topmost pointer as our new l. On the other hand, if h(ql ) ≤ h(qr ), we increment r
by one. We continue until either r leaves bj or S becomes empty. Then we push a pointer to the
rightmost element in bj onto S and proceed to the next sub-block.
In each round, the algorithm reads the complete input once in x-order. In addition, the algorithm
reads at most once each active sub-blocks in reverse order. Note that a sub-block becomes active
only once.

17

