Four Soviets Walk the Dog: Improved Bounds for Computing the Fréchet Distance

Discrete & Computational Geometry, Feb 2017

Given two polygonal curves in the plane, there are many ways to define a notion of similarity between them. One popular measure is the Fréchet distance. Since it was proposed by Alt and Godau in 1992, many variants and extensions have been studied. Nonetheless, even more than 20 years later, the original \(O(n^2 \log n)\) algorithm by Alt and Godau for computing the Fréchet distance remains the state of the art (here, n denotes the number of edges on each curve). This has led Helmut Alt to conjecture that the associated decision problem is 3SUM-hard. In recent work, Agarwal et al. show how to break the quadratic barrier for the discrete version of the Fréchet distance, where one considers sequences of points instead of polygonal curves. Building on their work, we give a randomized algorithm to compute the Fréchet distance between two polygonal curves in time \(O(n^2 \sqrt{\log n}(\log \log n)^{3/2})\) on a pointer machine and in time \(O(n^2(\log \log n)^2)\) on a word RAM. Furthermore, we show that there exists an algebraic decision tree for the decision problem of depth \(O(n^{2-\varepsilon })\), for some \(\varepsilon > 0\). We believe that this reveals an intriguing new aspect of this well-studied problem. Finally, we show how to obtain the first subquadratic algorithm for computing the weak Fréchet distance on a word RAM.

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:

https://link.springer.com/content/pdf/10.1007%2Fs00454-017-9878-7.pdf

Four Soviets Walk the Dog: Improved Bounds for Computing the Fréchet Distance

Four Soviets Walk the Dog: Improved Bounds for Computing the Fréchet Distance Kevin Buchin 0 1 2 3 Maike Buchin 0 1 2 3 Wouter Meulemans 0 1 2 3 Wolfgang Mulzer 0 1 2 3 0 Department of Mathematics and Computer Science , TU Eindhoven, PO Box 513, 5600 MB Eindhoven , The Netherlands 1 A preliminary version appeared as K. Buchin , M. Buchin, W. Meulemans, and W. Mulzer. Four Soviets Walk the Dog-with an Application to Alt's Conjecture. Proc. 25th SODA, pp. 1399-1413, 2014 2 Institut für Informatik, Freie Universität Berlin , 14195 Berlin , Germany 3 Mathematik, Ruhr Universität Bochum , 44780 Bochum , Germany Given two polygonal curves in the plane, there are many ways to define a notion of similarity between them. One popular measure is the Fréchet distance. Since it was proposed by Alt and Godau in 1992, many variants and extensions have been studied. Nonetheless, even more than 20 years later, the original O (n2 log n) algorithm by Alt and Godau for computing the Fréchet distance remains the state of the art (here, n denotes the number of edges on each curve). This has led Helmut Alt to conjecture that the associated decision problem is 3SUM-hard. In recent work, Agarwal et al. show how to break the quadratic barrier for the discrete version of the Fréchet distance, where one considers sequences of points instead of polygonal curves. Building on their work, we give a randomized algorithm to compute the Fréchet distance between two polygonal curves in time O (n2√log n(log log n)3/2) on a pointer machine and in - time O(n2(log log n)2) on a word RAM. Furthermore, we show that there exists an algebraic decision tree for the decision problem of depth O(n2−ε), for some ε > 0. We believe that this reveals an intriguing new aspect of this well-studied problem. Finally, we show how to obtain the first subquadratic algorithm for computing the weak Fréchet distance on a word RAM. 1 Introduction Shape matching is a fundamental problem in computational geometry, computer vision, and image processing. A simple version can be stated as follows: given a database D of shapes (or images) and a query shape S, find the shape in D that most resembles S. However, before we can solve this problem, we first need to address an issue: what does it mean for two shapes to be “similar”? In the mathematical literature, on can find many different notions of distance between two sets, a prominent example being the Hausdorff distance. Informally, the Hausdorff distance is defined as the maximal distance between two elements when every element of one set is mapped to the closest element in the other. It has the advantage of being simple to describe and easy to compute for discrete sets. In the context of shape matching, however, the Hausdorff distance often turns out to be unsatisfactory: it does not take the continuity of the shapes into account. There are well known examples where the distance fails to capture the similarity of shapes as perceived by human observers [6]. In order to address this issue, Alt and Godau introduced the Fréchet distance into the computational geometry literature [8,45]. They argued that the Fréchet distance is better suited as a similarity measure, and they described an O(n2 log n) time algorithm to compute it on a real RAM or pointer machine.1 Since Alt and Godau’s seminal paper, there has been a wealth of research in various directions, such as extensions to higher dimensions [7,23,26,28,33,46], approximation algorithms [9,10,37], the geodesic and the homotopic Fréchet distance [29,34,38,48], and much more [2,22,25, 35,36,51,54,55]. Most known approximation algorithms make further assumptions on the curves, and only an O(n2)-time approximation algorithm is known for arbitrary polygonal curves [24]. The Fréchet distance and its variants, such as dynamic timewarping [13], have found various applications, with recent work particularly focusing on geographic applications such as map-matching tracking data [15,63] and moving objects analysis [19,20,47]. Despite the large amount of published research, the original algorithm by Alt and Godau has not been improved, and the quadratic barrier on the running time of the associated decision problem remains unbroken. If we cannot improve on a quadratic bound for a geometric problem despite many efforts, a possible culprit may be the underlying 3SUM-hardness [44]. This situation induced Helmut Alt to make the following conjecture.2 1 For a brief overview of the different computational models in this paper, refer to Appendix. 2 Personal communication 2012, see also [6]. Conjecture 1.1 (Alt’s conjecture) Let P, Q be two polygonal curves in the plane. Then it is 3SUM-hard to decide whether the Fréchet distance between P and Q is at most 1. 1.1 Our Contribution We address the question by Agarwal et al. and show how to extend their approach to the Fréchet distance between two polygonal curves. Our algorithm requires total expected time O(n2√log n(log log n)3/2). This is the first algorithm with a running time of o(n2 log n) and constitutes the first improvement for the general case since the original paper by Alt and Godau [8]. To achieve this running time, we give the first subquadratic algorithm for the decision problem of the Fréchet distance. We emphasize that these algorithms run on a real RAM/pointer machine and do not require any bit-manipulation tricks. Therefore, our results are more in the line of Chan’s recent subcubic-time algorithms for all-pairs-shortest paths [30,31] or recent subquadratictime algorithms for min-plus convolution [16] than the subquadratic-time algorithms for 3SUM due to Baran et al. [12]. If we relax the model to allow constant time tablelookups, the running time can be improved to be almost quadratic, up to O(log log n) factors. As in Agarwal et al., our results are achieved by first giving a faster algorithm for the decision version, and then performing an appropriate search over the critical values to solve the optimization problem. In addition, we show that non-uniformly, the Fréchet distance can be computed in subquadratic time. More precisely, we prove that the decision version of the problem can be solved by an algebraic decision tree [11] of depth O(n2−ε), for some fixed ε > 0. It is, however, not clear how to implement this decision tree in subquadratic time, which hints at a discrepancy between the decision tree and the uniform complexity of the Fréchet problem. Finally, we consider the weak Fréchet distance, where we are allowed to walk backwards along the curves. In this case, our framework allows us to achieve a subquadratic algorithm on the work RAM. Refer to Table 1 for a comprehensive summary of our results. 3 It is well known that the four Russians are not actually Russian, so we refer to them as four Soviets in Table 1 Summary of results O n2√log n(log log n)3/2 We distinguish the continuous and the discrete Fréchet distance in the decision (D) and the computation (C) version. We also consider the weak continuous Fréchet distance. The computational models are the pointer machine (PM), the word RAM (WRAM) or the algebraic decision trees (DT). The old bounds are due to Alt and Godau [8] (continuous and weak) and Eiter and Mannila [39] (discrete) 1.2 Recent Developments Recently, Ben Avraham et al. [14] presented a subquadratic algorithm for the discrete Fréchet distance with shortcuts that runs in O (n4/3 log3 n) time. This running time resembles, at least superficially, our result on algebraic computation trees for the general discrete Fréchet distance. When we initially announced our results, we believed that they provided strong evidence that Alt’s conjecture is false. Indeed, for a long time it was conjectured that no subquadratic decision tree exists for 3SUM [57] and an (n2) lower bound is known in a restricted linear decision tree model [4, 40]. However, in a recent– and in our opinion quite astonishing–result, Grønlund and Pettie showed that if we allow only slightly more powerful algebraic decision trees than in the previous lower bounds, one can decide 3SUM non-uniformly in O (n2−ε) steps, for some fixed ε > 0 [52]. They also show that this leads to a general subquadratic algorithm for 3SUM, a situation very similar to the Fréchet distance as described in the present paper. Thus, despite some interesting developments, the status of Alt’s conjecture remains as open as before. However, we can now see that there exists a wide variety of efficiently solvable problems such as (in addition to 3SUM and the Fréchet distance) Sorting X + Y [41], Min- Plus- Convolution [16], or finding the Delaunay triangulation for a point set that has been sorted in two orthogonal directions [27], for which there seems to be a noticeable gap between the decision tree complexity and the uniform complexity. In our initial announcement, we also asked whether, besides 3SUM-hardness, there may be other reasons to believe that the quadratic running time for the Fréchet distance cannot be improved. Karl Bringmann provided an interesting answer to this question by showing that any algorithm for the Fréchet distance with running time O(n2−ε), for some fixed ε > 0, would violate the strong exponential time hypothesis (SETH) [17]. These results were later refined and improved to show that the lower bound holds in basically all settings (with the notable exception of the one-dimensional continuous Fréchet distance, which is still unresolved) [18]. We believe that these developments show that the Fréchet distance still holds many interesting aspects to be discovered and remains an intriguing object of further study. 2 Preliminaries and Basic Definitions Let P and Q be two polygonal curves in the plane, defined by their vertices p0, p1, . . . , pn and q0, q1, . . . , qn . Depending on the context, we interpret P and Q either as sequences of n and n edges, or as continuous functions P : [0, n] → R2 and Q : [0, n] → R2. In the latter case, we have P(i + λ) = (1 − λ) pi + λpi+1 for i = 0, . . . , n − 1 and λ ∈ [0, 1], and similarly for Q. Let be the set of all continuous and nondecreasing functions α : [0, 1] → [0, n] with α(0) = 0 and α(1) = n. The Fréchet distance between P and Q is defined as where · denotes the Euclidean distance. The classic approach to computing dF ( P, Q) uses the free-space diagram FSD( P, Q). It is defined as FSD( P, Q) := {(x , y) ∈ [0, n] × [0, n] | P(x ) − Q(y) ≤ 1}. In other words, FSD( P, Q) is the subset of the joint parameter space for P and Q where the corresponding points on the curves have distance at most 1, see Fig. 1. The structure of FSD( P, Q) is easy to describe. Let R := [0, n] × [0, n] be the ground set. We subdivide R into n2 cells C (i, j ) = [i, i + 1] × [ j, j + 1], for i, j = 0, . . . , n − 1. The cell C (i, j ) corresponds to the edge pair ei+1 and f j+1, where ei+1 is the (i + 1)th edge of P and f j+1 is the ( j + 1)th edge of Q. Then the set F (i, j ) := FSD( P, Q) ∩ C (i, j ) represents all pairs of points on ei+1 × f j+1 with distance at most 1. Elementary geometry shows that F (i, j ) is the intersection of C (i, j ) with an ellipse [8]. In particular, the set F (i, j ) is convex, and the intersection of FSD( P, Q) with the boundary of C (i, j ) consists of four (possibly empty) intervals, one on each side of ∂C (i, j ). We call these intervals the doors of C (i, j ) in FSD( P, Q). A door is said to be closed if the interval is empty, and open otherwise. A path π in FSD( P, Q) is bimonotone if it is both x - and y-monotone, i.e., every vertical and every horizontal line intersects π in at most one connected component. Alt and Godau observed that it suffices to decide whether there exists a bimonotone path from (0, 0) to (n, n) inside FSD( P, Q). We define the reachable region reach( P, Q) Fig. 1 Two polygonal curves P and Q, together with their associated free-space diagram. The reachable region reach(P, Q) is shown in blue. For example, the white area in C(2, 1), denoted F(2, 1), corresponds to all points on the third edge of P and the second edge of Q that have distance at most 1. As (5, 5) ∈ reach(P, Q), we have dF (P, Q) ≤ 1 as the set of points in FSD( P, Q) that are reachable from (0, 0) on a bimonotone path. Then, dF ( P, Q) ≤ 1 if and only if (n, n) ∈ reach( P, Q), see Fig. 1. It is not necessary to compute all of reach( P, Q): since FSD( P, Q) is convex inside each cell, we only need the intersections reach( P, Q) ∩ ∂C (i, j ). The sets defined by reach( P, Q) ∩ ∂C (i, j ) are subintervals of the doors of the free-space diagram, and they are defined by endpoints of doors in the free-space diagram in the same row or column. We call the intersection of a door with reach( P, Q) a reach-door. The reachdoors can be found in O(n2) time through a simple breadth-first-traversal of the cells [8]. In the next sections, we show how to obtain the crucial information, i.e., whether (n, n) ∈ reach( P, Q), in o(n2) time instead. 2.1 Basic Approach and Intuition In our algorithm for the decision problem, we basically want to compute reach( P, Q). But instead of propagating the reachability information cell by cell, we always group τ × τ cells (with 1 τ n) into an elementary box of cells. When processing a box, we can assume that we know which parts of the left and the bottom boundary of the box are reachable. That is, we know the reach-doors on the bottom and left boundary, and we need to compute the reach-doors on the top and right boundary of the elementary box. These reach-doors are determined by the combinatorial structure of the box. More specifically, suppose we know for every row and column the order of the door endpoints (including for the reach-doors on the left and bottom boundary). Then, we can deduce which of these door boundaries determine the reach-doors on the top and right boundary. We call the sequence of these orders, the (full) signature of the box. The total number of possible signatures is bounded by an expression in terms of τ . Thus, if we pick τ sufficiently small compared to n, we can pre-compute for all possible signatures the reach-doors on the top and right boundary, and build a data structure to query these quickly (Sect. 3). Since the reach-doors on the bottom and left boundary are required to make the signature, we initially have only incomplete signatures. In Sect. 4, we describe how to compute these efficiently. The incomplete signatures are then used to preprocess the data structure such that we can quickly find the full signature once we know the reach-doors of an elementary box. After building and preprocessing the data structure, it is possible to determine dF ( P, Q) ≤ 1 efficiently by traversing the free-space diagram elementary box by elementary box, as explained in Sect. 5. 3 Building a Lookup Table 3.1 Preprocessing an Elementary Box Before it considers the input, our algorithm builds a lookup table. As mentioned above, the purpose of this table is to speed up the computation of small parts of the free-space diagram. Let τ ∈ N be a parameter.4 The elementary box is a subdivision of [0, τ ]2 into τ columns and rows, thus τ 2 cells.5 For i, j = 0, . . . , τ − 1, we denote the cell [i, i + 1] × [ j, j + 1] with D(i, j ). We denote the left side of the boundary ∂ D(i, j ) by l(i, j ) and the bottom side by b(i, j ). Note that l(i, j ) coincides with the right side of ∂ D(i − 1, j ) and b(i, j ) with the top of ∂ D(i, j − 1). Thus, we write l(τ, j ) for the right side of D(τ − 1, j ) and b(i, τ ) for the top side of D(i, τ − 1). Figure 2 shows the elementary box. The door-order σ rj for a row j is a permutation of {s0, t0, . . . , sτ , tτ }, having 2τ + 2 elements. For i = 1, . . . , τ , the element si represents the lower endpoint of the door on l(i, j ), and ti represents the upper endpoint. The elements s0 and t0 are an exception: they describe the reach-door on the boundary l(0, j ) (i.e., its intersection with reach( P, Q)). The door-order σ rj represents the combinatorial order of these endpoints, as projected onto a vertical line, i.e., they are sorted into their vertical order. Some door-orders may encode the same combinatorial structure. In particular, when door i is closed, the exact position of si and ti in a door-order is irrelevant, as long as ti comes before si . For a closed door i (i > 0), we assign si to the upper endpoint of l(i, j ) and ti to the lower endpoint. The values of s0 and t0 are defined by the reachdoor and their relative order is thus a result of computation. We break ties between si and ti by placing si before ti , and any other ties are resolved by index. A door-order σic is defined analogously for a column i . We write x <ic y if x comes before y in σic, and x <rj y if x comes before y in σ rj . An incomplete door-order is a door-order in which s0 and t0 are omitted (i.e. the intersection of reach( P, Q) with the door is still unknown); see Fig. 3. 4 A preview for the impatient reader: we later set τ = (√log n/ log log n). 5 For now, the elementary box is a combinatorial concept. In the next section, we overlay these boxes on the free-space diagram to obtain “concrete” elementary boxes. Fig. 2 The elementary box. The cell D(2, 1) is shown white. Its boundaries—l(2, 1), l(3, 1), b(2, 1), b(2, 2)—are indicated t3 s3 Fig. 3 The door-order of a row (the vertical order of the points) encodes the combinatorial structure of the doors. The door-order for the row in the figure is s1s3s4t5t3t0s2t4s0s5t1t2. Note that s0 and t0 represent the reach-door, which is empty in this case. These are omitted in the incomplete door-order We can now define the (full) signature of the elementary box as the aggregation of the door-orders of its rows and columns. Therefore, a signature = (σ1c, . . . , στc, σ1r , . . . , στr ) consists of 2τ door-orders: one door-order σic for each column i and one doororder σ rj for each row j of the elementary box. Similarly, an incomplete signature is the aggregation of incomplete door-orders. For a given signature, we define the combinatorial reachability structure of the elementary box as follows. For each column i and for each row j , the combinatorial reachability structure indicates which door boundaries in the respective column or row define the reach-door of b(i, τ ) or l(τ, j ). Lemma 3.1 Let be a signature for the elementary box. Then we can determine the combinatorial reachability structure of the box in total time O(τ 2). Proof We use dynamic programming, very similar to the algorithm by Alt and Godau [8]. For each vertical edge l(i, j ) we define a variable l(i, j ), and for each horizontal edge b(i, j ) we define a variable b(i, j ). The l(i, j ) are pairs of the form (su , tv), representing the reach-door reach( P, Q) ∩ l(i, j ). If this reach-door is closed, then tv <rj su holds. If the reach-door is open, then it is bounded by the lower endpoint of the door on l(u, j ) and by the upper endpoint of the door on l(v, j ). (Note that in this case we have v = i .) Once again s0 and t0 are special and represent the reach-door on l(0, j ). The variables b(i, j ) are defined analogously. Fig. 4 The three cases for the recursive definition of l(i, j ). If the lower boundary is reachable, we can reach the whole right door (left). If neither the lower nor the left boundary is reachable, the right door is not reachable either (middle). Otherwise, the lower boundary is the maximum of l(i − 1, j ) and the lower boundary of the right door (right) Now we can compute l(i, j ) and b(i, j ) recursively as follows: first, we set l(0, j ) = b(i, 0) = (s0, t0), for i, j = 0, . . . , τ − 1. Next, we describe how to find l(i, j ) given l(i − 1, j ) and b(i − 1, j ), see Fig. 4. Case 1: Suppose b(i −1, j ) is open. This means that b(i −1, j ) intersects reach( P, Q), so reach( P, Q) ∩ l(i, j ) is limited only by the door on l(i, j ), and we can set l(i, j ) := (si , ti ). Case 2: If both b(i − 1, j ) and l(i − 1, j ) are closed, it is impossible to reach l(i, j ) and thus we set l(i, j ) := l(i − 1, j ). Case 3: If b(i − 1, j ) is closed and l(i − 1, j ) is open, we may be able to reach l(i, j ) via l(i − 1, j ). Let su be the lower endpoint of l(i − 1, j ). We need to pass l(i, j ) above su and si and below ti , and therefore set l(i, j ) := (max(su , si ), ti ), where the r maximum is taken according to the order < j . The recursion for the variable b(i, j ) is defined similarly. We can implement the recursion in time O (τ 2) for any given signature, for example by traversing the elementary box column by column, while processing each column from bottom to top. There are at most ((2τ + 2)!)2τ = τ O(τ 2) distinct signatures for the elementary box. We choose τ = λ√log n/ log log n for a sufficiently small constant λ > 0, so that this number becomes o(n). Thus, during the preprocessing stage we have time to enumerate all possible signatures and determine the corresponding combinatorial reachability structure inside the elementary box. This information is then stored in an appropriate data structure. 3.2 Building the Data Structure Before we describe this data structure, we first explain how the door-orders are represented. This depends on the computational model. By our choice of τ , there are o(n) distinct door-orders. On the word RAM, we represent each door-order and incomplete door-order by an integer between 1 and (2τ )!. This fits into a word of log n bits. On the pointer machine, we create a record for each door-order and incomplete door-order; we represent an order by a pointer to the corresponding record. The data structure has two stages. In the first stage, we assume we know the incomplete door-order for each row and for each column of the elementary box,6 and we wish to determine the incomplete signature. In the second stage we have obtained the reach-doors for the left and bottom sides of the elementary box, and we are looking for the full signature. The details of our method depend on the computational model. One way uses table lookup and requires the word RAM; the other way works on the pointer machine, but is a bit more involved. 3.2.1 Word RAM We organize the lookup table as a large tree T . In the first stage, each level of T corresponds to a row or column of the elementary box. Thus, there are 2τ levels. Each node has (2τ )! children, representing the possible incomplete door-orders for the next row or column. Since we represent door-orders by positive integers, each node of T may store an array for its children; we can choose the appropriate child for a given incomplete door-order in constant time. Thus, determining the incomplete signature for an elementary box requires O(τ ) steps on a word RAM. For the second stage, we again use a tree structure. Now the tree has O(τ ) layers, each with O(log τ ) levels. Again, each layer corresponds to a row or column of the elementary box. The levels inside each layer then implement a balanced binary search tree that allows us to locate the endpoints of the reach-door within the incomplete signature. Since there are 2τ endpoints, this requires O(log τ ) levels. Thus, it takes O(τ log τ ) time to find the full signature of a given elementary box. 3.2.2 Pointer machine Unlike in the word RAM model, we are not allowed to store a lookup table on every level of the tree T , and there is no way to quickly find the appropriate child for a given door-order. Instead, we must rely on batch processing to achieve a reasonable running time. Thus, suppose that during the first stage we want to find the incomplete signatures for a set B of m elementary boxes, where again for each box in B we know the incomplete door-order for each row and each column. Recall that we represent the door-order by a pointer to the corresponding record. With each such record, we store a queue of elementary boxes that is empty initially. We now simultaneously propagate the boxes in B through T , proceeding level by level. In the first level, all of B is assigned to the root of T . Then, we go through the nodes of one level of T , from left to right. Let v be the current node of T . We consider each elementary box b assigned to v. We determine the next incomplete door-order 6 In the next section, we describe how to determine the incomplete door-orders efficiently. for b, and we append b to the queue for this incomplete door-order—the queue is addressed through the corresponding record, so all elementary boxes with the same next incomplete door-order end up in the same queue. Next, we go through the nodes of the next level, again from left to right. Let v be the current node. The node v corresponds to a next incomplete door-order σ that extends the known signature of its parents. We consider the queue stored at the record for σ . By construction, the elementary boxes that should be assigned to v appear consecutively at the beginning of this queue. We remove these boxes from the queue and assign them to v . After this, all the queues are empty, and we can continue by propagating the boxes to the next level. During this procedure, we traverse each node of T a constant number of times, and in each level of the T we consider all the boxes in B. Since T has o(n) nodes, the total running time is O(n + mτ ). For the second stage, the data structure works just as in the word RAM case, because no table lookup is necessary. Again, we need O(τ log τ ) steps to process one box. After the second stage, we obtain the combinatorial reachability structure of the box in constant time since we precomputed this information for each box (Lemma 3.1). Thus, we have shown the following lemma, independently of the computational model. Lemma 3.2 For τ = λ√log n/ log log n, with a sufficiently small constant λ > 0, we can construct in o(n) time a data structure of size o(n) such that • given a set of m elementary boxes where the incomplete door-orders are known, we can find the incomplete signature of each box in total time O(n + mτ ); • given the incomplete signature and the reach-doors on the bottom and left boundary of an elementary box, we can find the full signature in O(τ log τ ) time; • given the full signature of an elementary box, we can find the combinatorial reachability structure of the box in constant time. 4 Preprocessing a Given Input Next, we perform a second preprocessing phase that considers the input curves P and Q. Our eventual goal is to compute the intersection of reach( P, Q) with the cell boundaries, taking advantage of the data structure from Sect. 3. For this, we aggregate the cells of FSD( P, Q) into (concrete) elementary boxes consisting of τ × τ cells. There are n2/τ 2 such boxes. We may avoid rounding issues by either duplicating vertices or handling a small part of FSD( P, Q) without lookup tables. The goal is to determine the signature for each elementary box S. At this point, this is not quite possible yet, since the signature depends on the intersection of reach( P, Q) with the lower and left boundary of S. Nonetheless, we can find the incomplete signature, in which the positions of s0, t0 (the reach-door) in the (incomplete) door-orders σir , σ jc are still to be determined. We aggregate the columns of FSD( P, Q) into vertical strips, each corresponding to a single column of elementary boxes (i.e., τ consecutive columns of cells in FSD( P, Q)). See Fig. 5. Let A be such a strip. It corresponds to a subcurve P of P with τ edges. The following lemma implies that we can build a data structure for A such that, given any Fig. 5 FSD(P, Q) is subdivided into n2/τ 2 elementary boxes of size τ × τ . The free-space diagram is subdivided into n2/τ 2 elementary boxes of size τ × τ . A strip is a column of elementary boxes: it corresponds to a subcurve of P with τ edges Fig. 6 By using the arrangement A defined by unit circles centered at vertices of P , we can determine the incomplete door-order of each segment s on Q. This is done by locating the dual point of s in the dual arrangement B. The dual arrangement also contains pseudolines to determine when s leaves a circle of A Lemma 4.1 Given a subcurve P with τ edges, we can compute in O (τ 6) time a data structure that requires O (τ 6) space and that allows us to determine the incomplete door-order of any line segment on Q in time O (log τ ). Proof Consider the arrangement A of unit circles whose centers are the vertices of P (see Fig. 6). The incomplete door-order of a line segment s is determined by the intersections of s with the arcs of A (and for a circle not intersecting s by whether s lies inside or outside of the circle). Let s be the line spanned by line segment s. Suppose we wiggle s . The order of intersections of s and the arcs of A changes only when s moves over a vertex of A or if s leaves or enters a circle. We use the standard duality transform that maps a line : y = ax + b to the point ∗ : (a, −b), and vice versa. Consider a unit circle C in A with center (cx , cy ). Elementary geometry shows that the set of all lines that are tangent to C from above dualizes to the curve ta∗(C ) : y = cx x − cy − √1 + x 2. Similarly, the lines that are tangent to C from below dualize to the curve tb∗(C ) : y = cx x − cy + √1 + x 2. Define C ∗ := {ta∗(C ), tb∗(C ) | C ∈ A}. Since any pair of distinct circles C1, C2 has at most four common tangents, one for each choice of above/below C1 and above/below C2, it follows that any two curves in C ∗ intersect at most once. Let V be the set of vertices in A, and let V ∗ be the lines dual to the points in V (note that |V | = O (τ 2)). Since for any vertex v ∈ V and any circle C ∈ A there are at most two tangents through v on C , each line in V ∗ intersects each curve in C ∗ at most once. Thus, the arrangement B of the curves in V ∗ ∪ C ∗ is an arrangement of pseudolines with complexity O (τ 4). Furthermore, it can be constructed in the same expected time, together with a point location structure that finds the containing cell in B of any given point in time O (log τ ) [60, Chap. 6.6.1]. Now consider a line segment s and the supporting line s . As observed in the first paragraph, the combinatorial structure of the intersection between s and A is completely determined by the cell of B that contains the dual point s∗. Thus, for every cell f (s) ∈ B, we construct a list L f (s) that represents the combinatorial structure of s ∩ A. There are O (τ 4) such lists, each having size O (τ ). We can compute L f (s) by traversing the zone of s in A. Since circles intersect at most twice and since a line intersects any circle at most twice, the zone has complexity O (τ 2α(τ )), where α(·) denotes the inverse Ackermann function [60, Thm. 5.11]. Since O (τ 2α(τ )) ⊂ O (τ 2), we can compute all lists in O (τ 6) time. Given the list L f (s), the incomplete door-order of s is determined by the position of the endpoints of s in L f (s). There are O (τ 2) possible ways for this, and we build a table T f (s) that represents them. For each entry in T f (s), we store a representative for the corresponding incomplete door-order. As described in the previous section, the representative is a positive integer in the word RAM model and a pointer to the appropriate record on a pointer machine. The total size of the data structure is O (τ 6) and it can be constructed in the same time. A query works as follows: given s, we can compute s∗ in constant time. Then we use the point location structure of B to find f (s) in O (log τ ) time. Using binary search on T f (s) (or an appropriate tree structure in the case of a point machine), we can then determine the position of the endpoints of s in the list L f (s) in O (log τ ) time. This bound holds both on the word RAM and on the pointer machine. Lemma 4.2 Given the data structure of Lemma 3.2, the incomplete signature for each elementary box can be determined in time O (nτ 5 + n2(log τ )/τ ). Proof By building and using the data structure from Lemma 4.1, we determine the incomplete door-order for each row in each vertical τ -strip in total time proportional to We repeat the procedure with the horizontal strips. Now we know for each elementary box in FSD( P, Q) the incomplete door-order for each row and each column. We use the data structure of Lemma 3.2 to combine these. As there are n2/τ 2 boxes, the number of steps is O(n2/τ + n) = O(n2/τ ). Hence, the incomplete signature for each elementary box is found in O(nτ 5 + n2(log τ )/τ ) steps. 5 Solving the Decision Problem With the data structures and preprocessing from the previous sections, we have all ingredients in place to determine whether dF ( P, Q) ≤ 1. We know for each elementary box its incomplete signature and we have a data structure to derive its full signature (and with it, the combinatorial reachability structure) when its reach-doors are known. What remains to be shown is that we can efficiently process the free-space diagram to determine whether (n, n) ∈ reach( P, Q). This is captured in the following lemma. Lemma 5.1 If the incomplete signature for each elementary box is known, we can determine whether (n, n) ∈ reach( P, Q) in time O(n2(log τ )/τ ). Proof We go through all elementary boxes of FSD( P, Q), processing them one column at a time, going from bottom to top in each column. Initially, we know the full signature for the box S in the lower left corner of FSD( P, Q). We use the signature to determine the intersections of reach( P, Q) with the upper and right boundary of S. There is a subtlety here: the signature gives us only the combinatorial reachability structure, and we need to map the resulting si , t j back to the corresponding vertices on the curves. On the word RAM, this can be done easily through table lookups. On the pointer machine, we use representative records for the si , ti elements and use O(τ ) time before processing the box to store a pointer from each representative record to the appropriate vertices on P and Q. We proceed similarly for the other boxes. By the choice of the processing order of the elementary boxes we always know the incoming reach-doors on the bottom and left boundary when processing a box. Given the incoming reach-doors, we can determine the full signature and find the structure of the outgoing reach-doors in total time O(τ log τ ), using Lemma 3.2. Again, we need O(τ ) additional time on the pointer machine to establish the mapping from the abstract si , ti elements to the concrete vertices of P and Q. In total, we spend O(τ log τ ) time per box. Thus, it takes time O(n2(log τ )/τ ) to process all boxes, as claimed. As a result, we obtain the following theorem for the pointer machine (and, by extension, for the real RAM model). For the word RAM model, we can obtain an even faster algorithm (see Sect. 6). Theorem 5.2 There is an algorithm that solves the decision version of the Fréchet problem in O(n2(log log n)3/2/√log n) time on a pointer machine. Proof Set τ = λ√log n/ log log n, for a sufficiently small constant λ > 0. The theorem follows by applying Lemmas 3.2, 4.2, and 5.1 in sequence. 6 Improved Bound on Word RAM We now explain how the running time of our algorithm can be improved if our computational model allows for constant time table-lookup. We use the same τ as above (up Fig. 7 A cluster consists of τ × τ elementary boxes, thus of τ 2 × τ 2 cells. A row R and its corresponding R for the central elementary box are indicated to a constant factor). However, we change a number of things. “Signatures” are represented differently and the data structure to obtain combinatorial reachability structures is changed accordingly. Furthermore, we aggregate elementary boxes into clusters and determine “incomplete door-orders” for multiple boxes at the same time. Finally, we walk the free-space diagram based on the clusters to decide dF ( P, Q) ≤ 1. 6.1 Clusters and Extended Signatures We introduce a second level of aggregation in the free-space diagram (see Fig. 7): a cluster is a collection of τ ×τ elementary boxes, that is, τ 2×τ 2 cells in FSD( P, Q). Let R be a row of cells in FSD( P, Q) of a certain cluster. As before, the row R corresponds to an edge e on Q and a subcurve P of P with τ 2 edges. We associate with R an ordered set Z = e0, z0, z1, z1, z2, z2, . . . , zk , zk , e1 with 2 · k + 3 elements. Here k is the number of intersections of e with the unit circles centered at the τ 2 vertices of P (all but the very first). Hence, k is bounded by 2τ 2 and |Z | is bounded by 4τ 2 + 3. The order of Z indicates the order of these intersections with e directed along Q. Elements e0 and e1 represent the endpoints of e and take a special role. In particular, these are used to represent closed doors and snap open doors to the edge e. The elements zi are placeholders for the positions of the endpoints of the reach-doors: z0 represents a possible reach-door endpoint between e0 and z1, the element z1 represents an endpoint between z1 and z2, etc. Consider a row R of an elementary box inside the row R of a cluster, corresponding to an edge e of Q. The door-index of R is an ordered set s0, t0, . . . , sτ , tτ of size 2τ + 2. Similar to a door-order, elements s0 and t0 represent the reach-door at the leftmost boundary of R ; the elements si and ti (1 ≤ i ≤ τ ) represent the door at the right boundary of the i th cell in R . However, instead of rearranging the set to indicate relative positions, the elements si and ti simply refer to elements in Z . If the door is open, they refer to the corresponding intersections with e (possibly snapped to e0 or e1). If the door is closed, si is set to e1 and ti is set to e0. The elements s0 and t0 are special, representing the reach-door, and they refer to one of the elements zi . An incomplete door-index is a door-index without s0 and t0. The advantage of a doorindex over a door-order is that the reach-door is always at the start. Hence, completing Fig. 8 (left) Every field represents the incomplete door-index of a row in an elementary box. (center) The fields are grouped into words per row in a cluster. (right) Transposition yields the desired organization, where a word represents the incomplete door-index of the rows in an elementary box an incomplete door-index to a full door-index can be done in constant time. Since a door-index has size 2τ + 2, the number of possible door-indices for R is τ O(τ ). We define the door-indices for the columns analogously. We concatenate the doorindices for the rows and the columns to obtain the indexed signature for an elementary box. Similarly, we define the incomplete indexed signature. The total number of possible indexed signatures remains τ O(τ 2). For each possible incomplete indexed signature we build a lookup table T as follows: the input is a word with 4τ fields of O(log τ ) bits each. Each field stores the positions in Z of the endpoints of the ingoing reach-doors for the elementary box: 2τ fields for the left side, 2τ fields for the lower side. The output consists of a word that represents the indices for the elements in Z that represent the outgoing reach-doors for the upper and right boundary of the box. Thus, the input of T is a word of O(τ log τ ) bits, and T has size τ O(τ ). Hence, for all incomplete indexed signatures combined, the size is τ O(τ 2) = o(n) by our choice of τ . 6.2 Preprocessing a Given Input During the preprocessing for a given input P, Q, we use superstrips consisting of τ strips. That is, a superstrip is a column of clusters and consists of τ 2 columns of the free-space diagram. Lemma 4.1 still holds, albeit with a larger constant c in place of 6. The data structure gets as input a query edge e, and it returns in O(log τ ) time a word that contains τ fields. Each field represents the incomplete door-index for e in the corresponding elementary box and thus consists of O(τ log τ ) bits. Hence, the word size is O(τ 2 log τ ) = O(log n) by our choice of τ . Thus, the total time for building a data structure for each superstrip and for processing all rows is O(n/τ 2 (τ c + n log τ )) = O(n2(log τ )/τ 2). We now have parts of the incomplete indexed signature for each elementary box packed into different words. To obtain the incomplete indexed signature, we need to rearrange the information such that the incomplete door-indices of the rows in one elementary box are in a single word. This corresponds to computing a transpose of a matrix, as is illustrated in Fig. 8. For this, we need the following lemma, which can be found—in slightly different form—in Thorup [62, Lem. 9]. Lemma 6.1 Let X be a sequence of τ words that contain τ fields each, so that X can be interpreted as a τ × τ matrix. Then we can compute in time O(τ log τ ) on a word RAM a sequence Y of τ words with τ fields each that represents the transpose of X . Proof The algorithm is recursive and solves a more general problem: let X be a sequence of a words that represents a sequence M of b different a × a matrices, such that the i th word in X contains the fields of the i th row of each matrix in M from left to right. Compute a sequence of words Y that represents the sequence M of the transposed matrices in M . The recursion works as follows: if a = 1, there is nothing to be done. Otherwise, we split X into the sequence X1 of the first a/2 words and the sequence X2 of the remaining words. X1 and X2 now represent a sequence of 2b (a/2) × (a/2) matrices, which we transpose recursively. After the recursion, we put the (a/2) × (a/2) submatrices back together in the obvious way. To finish, we need to transpose the off-diagonal submatrices. This can be done simultaneously for all matrices in time O(a), by using appropriate bit-operations (or table lookup). Hence, the running time obeys a recursion of the form T (a) = 2T (a/2) + O(a), giving T (a) = O(a log a), as desired. By applying the lemma to the words that represent τ consecutive rows in a superstrip, we obtain the incomplete door-indices of the rows for each elementary box. This takes total time proportional to We repeat this procedure for the horizontal superstrips. By using an appropriate lookup table to combine the incomplete door-indices of the rows and columns, we obtain the incomplete indexed signature for each elementary box in total time O(n2(log τ )/τ 2). 6.3 The Actual Computation We traverse the free-space diagram cluster by cluster (recall that a cluster consists of τ × τ elementary boxes). The clusters are processed column by column from left to right, and inside each column from bottom to top. Before processing a cluster, we walk along the left and lower boundary of the cluster to determine the incoming reach-doors. This is done by performing a binary search for each box on the boundary, and determining the appropriate elements zi which correspond to the incoming reachdoors. Using this information, we assemble the appropriate words that represent the incoming information for each elementary box. Since there are n2/τ 4 clusters, this step requires time O((n2/τ 4)τ 2 log τ ) = O(n2(log τ )/τ 2). We then process the elementary boxes inside the cluster, in a similar fashion. Now, however, we can process each elementary box in constant time through a single table lookup, so the total time is O(n2/τ 2). Hence, the total running time of our algorithm is O(n2(log τ )/τ 2). By our choice of τ = λ√log n/ log log n for a sufficiently small λ > 0, we obtain the following theorem. Theorem 6.2 The decision version of the Fréchet problem can be solved in O(n2(log log n)2/ log n) time on a word RAM. 7 Computing the Fréchet Distance The optimization version of the Fréchet problem, i.e., computing the Fréchet distance, can be done in O(n2 log n) time using parametric search with the decision version as a subroutine [8]. We showed that the decision problem can be solved in o(n2) time. However, this does not directly yield a faster algorithm for the optimization problem: if the running time of the decision problem is T (n), parametric search gives an O((T (n)+n2) log n) time algorithm [8]. There is an alternative randomized algorithm by Raichel and Har-Peled [49]. Their algorithm also needs O((T (n) + n2) log n) time, but below we adapt it to obtain the following lemma. Lemma 7.1 The Fréchet distance of two polygonal curves with n vertices each can be computed by a randomized algorithm in O(n22α(n) + T (n) log n) expected time, where T (n) is the time for the decision problem. Before we prove the lemma, we recall that possible values of the Fréchet distance are limited to a certain set of critical values [8]: 1. the distance between a vertex of one curve and a vertex of the other curve (vertex– vertex); 2. the distance between a vertex of one curve and an edge of the other curve (vertex– edge); and 3. for two vertices of one curve and an edge of the other curve, the distance between one of the vertices and the intersection of e with the bisector of the two vertices (if this intersection exists) (vertex–vertex–edge). If we also include vertex–vertex–edge tuples with no intersection, we can sample a critical value uniformly at random in constant time. The algorithm now works as follows (see Har-Peled and Raichel [49] for more details): first, we sample a set S of K = 4n2 critical values uniformly at random. Next, we find a , b ∈ S such that the Fréchet distance lies between a and b and such that [a , b ] contains no other value from S. In the original algorithm this is done by sorting S and performing a binary search using the decision version. Using median-finding instead, this step can be done in O(K + T (n) log K ) time. Alternatively, the running time of this step could be reduced by picking a smaller K . However, this does not improve the final bound, since it is dominated by a O(n22α(n)) term. The interval [a , b ] with high probability contains only a small number of the remaining critical values. More precisely, for K = 4n2 the probability that [a , b ] has more than 2cn ln n critical values is at most 1/nc [49, Lem. 6.2]. The remainder of the algorithm proceeds as follows: first, we find all critical values of type vertex–vertex and vertex–edge that lie inside the interval [a , b ]. This can be done in O(n2) time by checking all vertex–vertex and vertex–edge pairs. Among these values, we again use median-finding to determine the interval [a, b] ⊆ [a , b ] that contains the Fréchet distance in O(K + T (n) log K ) time. It remains to determine the critical values corresponding to vertex–vertex–edge tuples that lie in [a, b]. For this, take an edge e of P and the vertices of Q. Conceptually, we start with circles of radius a around the vertices of Q, and we increase the radii until b. During this process, we observe the evolution of the intersection points between the circle arcs and e. Because all vertex–vertex and vertex–edge events have been eliminated, each circle intersects e in either 0 or 2 points, and this does not change throughout the process. A critical value of vertex–vertex–edge type corresponds to the event that two different circles intersect e in the same point, i.e., that two intersection points meet while growing the circles. Two intersection points can meet at most once, and when they do, they exchange their order along e. This suggests the following algorithm: let Aa be the arrangement of circles with radius a around the vertices of Q, and let Ab be the concentric arrangement of circles with radius b. We determine the ordered sequence Ia of the intersection points of the circles in Aa with e, and we number them in their order along e. Next, we find the ordered sequence of intersection points Ib between e and the circles in Ab. We assign to each point in Ib the number of the corresponding intersection points in Ia . Since |Ia | = |Ib|, this gives a permutation of {1, . . . , |Ia |}, Two intersection points change their order from Ia to Ib exactly if there is a vertex–vertex–edge event in [a, b], so these events correspond to the inversions of the resulting permutation. Given that there are k such inversions, we can find them in time O(|Ia | + k) using insertion sort. Thus, the overall running time to find the critical events in [a, b], ignoring the time for computing Ia and Ib, is O(n2 + K ). It remains to show that we can quickly find Ia and Ib. We describe the algorithm for Ia . First, compute the arrangement Aa of circles with radius a around the vertices of Q. This takes O(n2) time [32]. To find the intersection order, traverse in Aa the zone of the line spanned by e. The time for the traversal is bounded by the complexity of the zone. Since the circles pairwise intersect at most twice and intersects each circle only twice, the complexity of the zone is O(n2α(n)) [60, Thm. 5.11]. Summing over all edges e, this adds a total of O(n22α(n)) to the running time. To find Ib, we proceed similarly with Ab. Thus the overall time is O(T (n) log(n) + n22α(n) + K ). The event K > 8n ln n has probability less than 1/n4, and we always have K = O(n3). Thus, this case adds o(1) to the expected running time. Given K ≤ 8n ln n, the running time is O(n log n). Lemma 7.1 follows. Theorem 7.2 now results from Lemma 7.1, Theorem 5.2, and Theorem 6.2. Theorem 7.2 The Fréchet distance of two polygonal curves with n edges each can be computed by a randomized algorithm in time O(n2√log n(log log n)3/2) on a pointer machine and in time O(n2(log log n)2) on a word RAM. 8 Discrete Fréchet Distance on the Pointer Machine As mentioned in the introduction, Agarwal et al. [1] give a subquadratic algorithm for finding the discrete Fréchet distance between two point sequences, using the word RAM. In this section, we explain how their algorithm for the decision version of the problem can be adapted to the pointer machine. This shows that, at least for the decision version, the speed-up does not come from bit-manipulation tricks but from a deeper understanding of the underlying geometric structure. Our presentation is slightly different from Agarwal et al. [1], in order to allow for a clearer comparison with our continuous algorithm. We recall the problem definition: we are given two sequences P = p1, p2, . . . , pn and Q = q1, q2, . . . , qn of n points in the plane. For δ > 0, we define a directed graph Gδ with vertex set P × Q. In Gδ, there is an edge between two vertices ( pi , q j ), ( pi , q j+1) if and only if both d( pi , q j ) ≤ δ and d( pi , q j+1) ≤ δ. The condition is similar for an edge between vertices ( pi , q j ) and ( pi+1, q j ), and vertices ( pi , q j ) and Fig. 9 The discrete Fréchet distance: (left) two point sequences P (disks) and Q (crosses) with 5 points each; (middle) the associated free-space matrix F (white = 1, gray = 0); (right) the resulting reachability matrix M. Since M55 = 1, the discrete Fréchet distance is at most 1 ( pi+1, q j+1). There are no further edges in Gδ. The discrete Fréchet distance between P and Q is the smallest δ for which Gδ has a path from ( p1, q1) to ( pn, qn). In the decision version of the problem, we are given δ > 0, and we need to decide whether there is a path from ( p1, q1) to ( pn, qn) in Gδ. We now describe a subquadratic pointer machine algorithm for the decision version. Thus, let point sequences P, Q be given, and suppose without loss of generality that δ = 1. The discrete analogue of the free-space diagram is an n × n Boolean matrix F where Fi j = 1, if d( pi , q j ) ≤ 1, and Fi j = 0, otherwise, for i, j = 1, . . . , n. We call F the free-space matrix. Similarly, the discrete analogue of the reachable region is an n × n Boolean matrix M that is defined recursively as follows: M11 = F11, and for i, j = 1, . . . , n, (i, j ) = 1, we have Mi j = 1 if and only if Fi j = 1 and at least one of Mi−1, j , Mi, j−1 or Mi−1, j−1 equals 1 (we set Mi,0 = M0, j = 0, for i, j = 1, . . . n). Then the discrete Fréchet distance between P and Q is at most 1 if and only if Mnn = 1. We call M the reachability matrix; see Fig. 9. Adapting the method of Agarwal et al. [1], we show how to use preprocessing and table lookup in order to decide whether Mnn = 1 in o(n2) steps on a pointer machine. Let τ = λ log n, for a suitable constant λ > 0. We subdivide the rows of M into k = O(n/τ ) strips, each consisting of τ consecutive rows: the first strip L1 consists of rows 1, . . . , τ , the second strip L2 consists of rows τ + 1, . . . , 2τ , and so on. Each strip Li , i = 1, . . . , k, corresponds to a contiguous subsequence Pi of τ points on P. Let Ai be the arrangement of disks obtained by drawing a unit disk around each vertex in Pi . The arrangement Ai has O(τ 2) faces. Next, let ρ = λ log n/ log log n, with λ > 0 as above. We subdivide each strip Li , i = 1, . . . , k into l = O(n/ρ) elementary boxes, each consisting of ρ consecutive columns in Li . We label the elementary boxes as Bi j , for i = 1, . . . , k and j = 1, . . . , l. As above, an elementary box Bi j has corresponding contiguous subsequences Pi of τ vertices on P and Q j of ρ vertices on Q. Now, the incomplete signature of an elementary box Bi j consists of (i) the index i of the strip that contains it; and (ii) the sequence f1, f2, . . . , fρ of faces in the disk arrangement Ai that contain the ρ vertices of Q j , in that order. The full signature of an elementary box Bi j consists of its incomplete signature plus a sequence of ρ + τ bits, that represent the entries in the reach matrix M directly above and to the left of Bi j . We call these bits the reach bits. 2 2 3 4 1 3 Fig. 10 (left) We subdivide the reachability matrix M into strips of τ rows. Each strip is subdivided into elementary boxes of ρ columns. Each elementary box corresponds to a subsequence of length τ on P and a subsequence of length ρ on Q. (right) The incomplete signature σ of an elementary box Bi j consists of the index of the containing strip and the face sequence for Q j in Ai . The full signature σ f additionally contains ρ + τ reach bits that represent the bits in M directly above and to the left of Bi j As in the continuous case, the information in the full signature suffices to determine how the reachability information propagates through the elementary box; see Fig. 10. The preprocessing phase proceeds as follows: first, we enumerate all possible incomplete signatures. For this, we need to compute all strips Li and the corresponding disk arrangements Ai , for i = 1, . . . , τ . Furthermore, we also compute a suitable point location structure for each Ai . Since there are O (n/τ ) strips, each of which consists of τ rows, this takes time O ((n/τ ) · τ 2 log τ ) = O (nτ log τ ) = O (n log n log log n). For each strip Li , since Ai has O (τ 2) faces, the number of possible face sequences f1, . . . , fρ is τ O(ρ) ≤ n1/3, by our choice of τ and ρ and for λ small enough. Thus, there are O (n4/3/ log n) incomplete signatures, and they can be enumerated in the same time. Now, for each incomplete signature σ = (i, f1, . . . , fρ ) we build a lookup-table that encodes for each possible setting of the reach bits the resulting reach bits at the bottom and the right boundary of the elementary box. There are 2ρ+τ ≤ n1/3 possible settings of the reach bits, by our choice of τ and ρ and for λ small enough. We enumerate all of them and organize them as a complete binary tree of depth ρ + τ . For each setting of the reach bits, we use the information of the incomplete signature to determine the result through a straightforward dynamic programming algorithm [1, 18, 39] in O (τ · ρ) = O (log2 n) time, and we store the result as a linked list of length ρ + τ − 1 at the leaf for the corresponding reach bits; see Fig. 11. Thus, the total time for this part of the preprocessing phase is O (n5/3 log n). Next, we determine for each elementary box Bi j its incomplete signature. For this, we use the point location structure for Ai to determine for each vertex in Q j the face of Ai that contains it. There are O (n2/τρ) elementary boxes, each Q j has ρ vertices, and one point location query takes O (log τ ) time, so the total time for this step is O ((n2/τρ) · ρ · log τ ) = O ((n2/ log n) log log n). Using this information, we can store with each elementary box a pointer to the lookup table for the corresponding incomplete signature. 1 0 1 0 1 0 Fig. 11 For each incomplete signature, we create a lookup table organized as a complete binary tree. Each leaf corresponds to a setting of the reach bits for the elementary box. In the leaves, we store a linked list of length ρ + τ − 1 that represents the contents of the reachability matrix at the bottom and at the right of the elementary box Finally, we can now use the lookup tables to propagate the reachability information through M , one elementary box at a time, as in the continuous case. The time to process one elementary box is O (ρ +τ ), because we need to traverse the corresponding lookup table to find the reach bits for the adjacent boxes. Thus, the total running time is O ((n2/τρ) · (ρ + τ )) = O (n2/ρ) = O ((n2/ log n) log log n). Thus, we get the following pointer machine version of the result by Agarwal et al. [1] Theorem 8.1 There is an algorithm that solves the decision version of the discrete Fréchet problem in O ((n2/ log n) log log n) time on a pointer machine. Remark Agarwal et al. [1] further describe how to get a faster algorithm for the decision version by aggregating the elementary boxes into larger clusters, similar to the method given in Sect. 6. This improved algorithm finally leads to a subquadratic algorithm for computing the discrete Fréchet distance. Unfortunately, as in Sect. 6, it seems that this improvement crucially relies on constant time table lookup, so it does not directly translate to the pointer machine. 9 Decision Trees Our results also have implications for the decision-tree complexity of the Fréchet problem. Since in that model we account only for comparisons between input elements, the preprocessing comes for free, and hence the size of the elementary boxes can be increased. Before we consider the continuous Fréchet problem, we first note that a similar result can be obtained easily for the discrete Fréchet problem. Proof First, we consider the decision version: we are given two sequences P = p1, . . . , pn and Q = q1, . . . , qn of n points in the plane, and we would like to decide whether the discrete Fréchet distance between P and Q is at most 1. Katz and Sharir [53] showed that we can compute a representation of the set of pairs ( pi , q j ) with pi − q j ≤ 1 in O(n4/3) steps. This information suffices to complete the reachability matrix without further comparisons. As shown by Agarwal et al. [1], one can then solve the optimization problem at the cost of another O(log n)-factor, which is absorbed into the O-notation. Given our results above, we prove an analogous statement for the continuous Fréchet distance. Theorem 9.2 There exists an algebraic decision tree for the Fréchet problem (decision version) of depth O(n2−ε), for a fixed constant ε > 0. Proof We reconsider the steps of our algorithm. The only phases that actually involve the input are the second preprocessing phase and the traversal of the elementary boxes. The reason of our choice for τ was to keep the time for the first preprocessing phase small. This is no longer a problem. By Lemmas 4.2 and 5.1, the remaining cost is bounded by O(nτ 5 + n2(log τ )/τ ). Choosing τ = n1/6, we get a decision tree of depth n · n5/6 + n2−1/6 log n. This is O(n2−(1/6) log n) = O(n2−ε), for any fixed 0 < ε < 1/6. 10 Weak Fréchet Distance The weak Fréchet distance is a variant of the Fréchet distance where we are allowed to walk backwards along the curves [8]. More precisely, let P and Q be two polygonal curves, each with n edges, and let be the set of all continuous functions α : [0, 1] → [0, n] with α(0) = 0 and α(1) = n. The weak Fréchet distance between P and Q is defined as Compared to the regular Fréchet distance, the set now also contains non-monotone functions. The weak Fréchet distance was also introduced by Alt and Godau [8], who showed how to compute it in O(n2 log n) worst-case time. We will now use our framework to obtain an algorithm that runs in o(n2) expected time on a word RAM. 10.1 A Decision Algorithm for the Pointer Machine As usual, we start with the decision version: given two polygonal curves P and Q, each with n edges, decide whether dwF( P, Q) ≤ 1. This has an easy interpretation in Fig. 12 The polygonal curves P and Q have weak Fréchet distance at most 1, but Fréchet distance larger than 1: the point (n, n) is not in reach(P, Q), but the vertices for the cells C(0, 0) and C(4, 4) are in the same connected component of G. The reachable region is shown dark blue, the edges if G are shown light blue terms of the free-space diagram. Define an undirected graph G = (V , E ) with vertex set V = {0, 1, . . . , n − 1}2. The vertex (i, j ) ∈ V corresponds to the cell C (i, j ), and there is an edge between two vertices (i, j ) and (i , j ) if and only if the two cells C (i, j ) and C (i , j ) are neighboring (i.e., if |i − i | + | j − j | = 1) and the door between them is open. Then, dwF( P, Q) ≤ 1 if and only if (i) | P(0) − Q(0)| ≤ 1; (ii) | P(n) − Q(n)| ≤ 1; and (iii) the vertices (0, 0) and (n − 1, n − 1) are in the same connected component of G, see Fig. 12. Let τ, ρ ∈ N be parameters, to be determined later. We subdivide the cells into k = O(n/τ ) vertical strips L1, . . . , Lk , each consisting of τ consecutive columns. Each strip Li corresponds to a subcurve Pi of P with τ edges. For each such subcurve Pi , we define two arrangements Ai and Bi . To obtain Ai , we take for each edge e of Pi the “stadium” ce of points with distance exactly 1 from e, and we compute the resulting arrangement. Since two distinct curves ce, ce cross in O(1) points, the complexity of Ai is O(τ 2), see Fig. 13. The arrangement Bi is the arrangement B described in the proof of Lemma 4.1, i.e., the arrangement of the curves dual to the tangent lines for the unit circles around the vertices of Pi . Next, we subdivide each strip into = O(n/ρ) elementary boxes, each consisting of ρ consecutive rows. We label the elementary boxes as Bi j , with 1 ≤ i ≤ k, 1 ≤ j ≤ . The rows of an elementary box Bi j correspond to a subcurve Q j of Q with ρ edges. The signature of Bi j consists of (i) the index i of the corresponding strip; (ii) for each vertex of Q j the face of Ai that contains it; and (ii) for each edge e of Q j the face of Bi that contains the point that is dual to the supporting line of e, plus two indices a, b ∈ {1, . . . , τ } that indicate the first and the last unit circle around a vertex of Pi that e intersects, as we walk from one endpoint to another. Given an elementary box B, the connection graph G B of B has τρ vertices, one for each cell in B, and an edge between two cells C , C of B if and only if C and C share a (horizontal or vertical) edge with an open door. The connectivity list of B is a linked list with 2τ + 2ρ − 4 entries that stores for each cell C on the boundary of B a pointer to a record that represents the connected component of the connection graph G B that contains C . Lemma 10.1 There are O(nτ 8ρ+1) different signatures. The connection graph of an elementary box Bi j depends only on its signature, and the connectivity list can be computed in O(τρ) time on a pointer machine, given the signature. Proof First, we count the signatures. There are O(n/τ ) strips. Once the strip index i is fixed, a signature consists of ρ + 1 faces of Ai , ρ faces of Bi , and ρ pairs of indices a, b ∈ {1, . . . , τ }. The arrangement Ai has O(τ 2) faces, and the arrangement Bi has O(τ 4) faces, as explained in the proof of Lemma 4.1. Finally, there are τ 2 pairs of indices. Thus, the number of possible signatures in one strip is τ 8ρ+2. In total, we get O(nτ 8ρ+1) signatures. Next, the connection graph of an elementary box Bi j is determined solely by which doors are open and which doors are closed. We explain how to deduce this information from the signature. A horizontal edge of Bi j corresponds to an edge e of Pi and a vertex q of Q j . The door is open if and only if q has distance at most 1 from e. This is determined by the face of Ai containing q. Similarly, a vertical edge of Bi j corresponds to a vertex p of Pi and an edge e of Q j . The door is open if and only if e intersects the unit circle with center p. As in Lemma 4.1, this can be inferred from the face of Bi that contains the point dual to the supporting line of e, together with the indices (a, b) of the first and last circle intersected by e. Finally, given the signature, we can build the connection graph G Bi j in O(τρ) time, assuming that the arrangements Ai and Bi provide suitable data structures. With G Bi j at hand, the connection list can be found in O(τρ) steps, using breadth first search. As usual, our strategy now is to preprocess all possible signatures and to determine the signature of each elementary box. Using this information, we can then process the elementary boxes quickly in our main algorithm. The next lemma describes the preprocessing steps. Lemma 10.2 We can determine for each elementary box Bi j a pointer to its connectivity list in total time O(nτ 8ρ+2ρ + (n2/τ ) log τ ) on a pointer machine. Proof First, we compute the arrangements Ai and Bi for each vertical strip. By Lemma 4.1, this takes O(τ 6) steps per strip, for a total of O(τ 6 · n/τ ) = O(nτ 5). Then, we enumerate all possible signatures, and we compute the connectivity list for each of them. By Lemma 10.1, this needs O(nτ 8ρ+2ρ) time. We store the signatures and their connectivity lists Next, we determine the signature for each elementary box. Fix a strip Li . For each vertex q and each edge e of Q, we determine the containing faces of Ai and Bi and the pair (a, b) that represents the first and least intersection of e with the unit circles around the vertices of Pi . This takes O(n log τ ) steps in total, using appropriate point location structures for the arrangements. Now, we use the data structure from the preprocessing to connect each elementary box to its connectivity list. A simple pointer-based structure supports one lookup in O(ρ log τ ) time. Since there are O(n/ρ) elementary boxes in one strip, we get a total running time of O(n log τ ) per strip. Since there are O(n/τ ) strips, the resulting running time is O((n2/τ ) log τ ). With the information from the preprocessing phase, we can easily solve the decision problem with a union-find data structure. Proof Let S be the set that contains the boundary cells of all elementary boxes. Then, |S| = O n2(ττρ+ρ) . We create a Union-Find data structure for S. Then, we go through all elementary boxes Bi j , and for each Bi j , we use the connectivity list to connect those subsets of boundary cells that are connected inside Bi j . This can be done with The following theorem summarizes our algorithm for the decision problem. Theorem 10.4 Let P and Q be two polygonal curves with n edges each. We can decide whether dwF ( p, q) ≤ 1 in O((n2/ log n)α(n) log log n) time on a pointer machine, where α(·) denotes the inverse Ackermann function. Proof We set τ = log n and ρ = λ log n/ log log n, for a suitable constant λ > 0. If λ is small enough, then τ 8ρ+2 = O(n4/3), and the algorithm from Lemma 10.2 runs in time O((n2/ log n) log log n). After the preprocessing is finished, we can use Lemma 10.3 to obtain the final result in time O((n2/ log n)α(n) log log n). 10.2 A Faster Decision Algorithm on the Word RAM Theorem 10.4 is too weak to obtain a subquadratic algorithm for the weak Fréchet distance. This requires the full power of the word RAM. Theorem 10.5 Let P and Q be two polygonal curves with n edges each. We can decide whether dwF ( p, q) ≤ 1 in O((n2/ log2 n)(log log n)5) time on a word RAM. Proof We set τ = λ log2 n/ log log n and ρ = λ log n/ log log n. If λ is small enough, then τ 8ρ+2 = O(n4/3), and we can perform the algorithm from Lemma 10.2 in time O((n2/ log2 n) log log n). We modify the algorithm from Lemma 10.2 slightly. Instead of a pointer-based connectivity list, we compute a packed connectivity list. It consists of O(log n) words of log n bits each. Each word stores (log n/ log log n) entries of O(log log n) bits. An entry consists of O(1) fields, each with O(log log n) bits. As in Lemma 10.2, the entries in the packed connectivity list represent the connected components of the connection graph G B for the boundary cells of a given elementary box B. This is done as follows: each entry of the connectivity list corresponds to a boundary cell C of B. In the first field, we store a unique identifier from {1, . . . , 2τ + 2ρ − 4} that identifies C . In the second field, we store the smallest identifier of any boundary cell of B that lies in the same connected component as C . The remaining fields of the entry are initialized to 0. Furthermore, we compute for each elementary box B a sequence of O((τ / log n) log log n) words that indicates for each boundary cell of B whether the corresponding door to the neighboring elementary box is open. This information can be obtained in the same time by adapting the algorithm from Lemma 10.2. Next, we group the elementary boxes into clusters. A cluster consists of log n vertically adjacent elementary boxes from a single strip. The first set of clusters come from bottommost log n elementary boxes, the second set of clusters from following log n elementary boxes, etc. The boundary and the connectivity list of a cluster are defined analogously as for an elementary box. Below, in Lemma 10.6, we show that we can compute the (pointer-based) connectivity list of a cluster in time O((τ + ρ)(log log n)4). Then, the lemma follows: there are O(n2/(τρ log n)) clusters, so the total time to find connectivity lists for all clusters is O(n2/(ρ log n)(log log n)4+ n2/(τ log n)(log log n)4) = O((n2/ log2 n)(log log n)5). After that, we can solve the decision problem in time = O as in Lemma 10.3. Lemma 10.6 Given a cluster, we can compute the connectivity list for its boundary in total time O((τ + ρ)(log log n)4) on a word RAM. Proof Our algorithm makes extensive use of the word RAM capabilities. Data is processed in packed form. That is, we usually deal with a sequence of packed words, each with O(log n/ log log n) entries. Each entry has O(1) fields of O(log log n) bits. We restrict the number of entries in each word so that the total number of bits used is at most, say, (1/3) log n. Then, we can support any reasonable binary operation on two packed words in O(1) time after o(n) preprocessing time. Indeed, we only need to build a lookup-table during the preprocessing phase. First, we observe that we can sort packed data efficiently. Claim 10.7 Suppose we are given a sequence of μ packed words, representing a sequence of O(μ(log n/ log log n)) entries. We can obtain a sequence of μ packed words that represents the same entries, sorted according to any given field, in O(μ log μ) time. Proof We adapt usual techniques for packed integer sorting [5]. We precompute a merge-operation that receives two packed words, each with their entries sorted according to a given field, and returns two packed words that represent the sorted sequence of all entries. We also precompute an operation that receives one packed word and returns a packed word with the same entries, sorted according to a given field. Then, we can perform a merge sort in O(μ log μ) time, because in each level of the recursion, the total time for merging is O(μ). Once we are down to a single word, we can sort its entries in one step. More details can be found in the literature on packed sorting [5,27]. Now we describe the strategy of the main algorithm. As explained above, a cluster consists of log n vertically adjacent elementary boxes. We number the boxes B1, . . . , Blog n, from bottom to top. For i = 1, . . . , log n, we denote by Gi the connection graph for the boxes B1, . . . , Bi . That is, Gi is obtained by taking the union ij=1 G B j of the individual connection graphs and adding edges for adjacent cells in neighboring boxes that share an open door. Our algorithm proceeds in two phases. In the first phase, we propagate the connectivity information upwards. That is, for i = 1, . . . , log n, we compute a sequence Wi of O((τ + ρ) log log n/ log n) packed words. Each entry in Wi corresponds to a cell C on the boundary of Bi . The first field stores a unique identifier for the cell, and the second entry stores an identifier of the connected component of C in Gi . In the second phase, we go in the reverse direction. For i = log n, . . . , 1, we compute a sequence Wi of O((τ + ρ) log log n/ log n) packed words that store for each cell C on the boundary of Bi an identifier of the connected component of C in Glog n. Once this information is available, the (pointerbased) connectivity list for the cluster boundary can be extracted in O(τ + ρ) time, using appropriate precomputed operations on packed words, see Fig. 14 . We begin with the first phase. The upper boundary of an elementary box are the τ cells in the topmost row of the box, the lower boundary are the τ cells in the bottommost row. From the preprocessing phase, we have a pointer to the packed connectivity lists for all elementary boxes B1, . . . , Blog n. We make local copies of Fig. 14 The algorithm for computing the connected components within a cluster. In the beginning, we know the components inside each elementary box (left). In the first phase (top right), we update the connectivity information from bottom to top, by considering the upper and lower boundaries of the vertically adjacent elementary boxes. In the second phase (bottom right), we propagate the connected components from top to bottom these packed connectivity lists, and we modify them so that the identifiers of all cells and connected components are unique. For example, we can increase all identifiers in the connectivity list for Bi by i log3 n. Using lookup-tables for these word-operations, this step can be carried out in O((τ + ρ) log log n) time for the whole cluster. The sequence W1 for B1 is exactly the packed connectivity list of G1. Now, suppose that we have computed for Bi a sequence Wi of packed words that represent the connected component in Gi for each boundary cell of Bi . We need to compute a similar sequence Wi+1 for Bi+1. For this, we need to determine how the doors between Bi and Bi+1 affect the connected components in G Bi+1 . Let H = (VH , E H ) be the graph that has one vertex for each connected component of Gi that contains a cell at the upper boundary of Bi and one vertex for each connected component of G Bi+1 that contains a cell at the lower boundary of Bi+1. The vertices of H are labeled by the identifiers of their corresponding components. Suppose that v ∈ VH represents a component Dv in Gi and w ∈ VH represents a component Dw in G Bi+1 . There is an edge between v and w in E H if and only if Dv contains a cell on the upper boundary of Bi and Dw contains a cell on the lower boundary of Bi+1 such that these two cells share an open door. The graph H has O(τ ) vertices and O(τ ) edges. Our goal is to obtain the connected components of H . From this, we will then Fig. 15 The algorithm of Hirschberg, Chandra, and Sarwate. In each step, we either connect all nodes to their largest neighbor or to their smallest neighbor. Then, we contract the resulting rooted trees and label the contracted nodes with the index of the smallest node derive the next sequence Wi+1. To obtain the components of H , we adapt the parallel connectivity algorithm of Hirschberg, Chandra, and Sarwate [50]. Proof First, we obtain a list of O((τ / log n) log log n) packed words that represent the vertices of H , as follows: we extract from Wi and from the connectivity list of Bi+1 the entries for the upper boundary of Bi and for the lower boundary of Bi+1. If all these entries are stored together in the respective packed lists, this takes O((τ / log n) log log n) time to copy the relevant data words. Then, we sort these packed sequences according to the component identifier, in time O(τ / log n) log log2 n) (Claim 10.7). Finally, we go over the sorted packed sequence to extract a sorted lists of the distinct component identifiers. This takes O((τ / log n) log log n) time, using an appropriate operation on packed words. The edges of E H can also be represented by a sequence of O((τ / log n) log log n) words. As explained above, we have from the preprocessing phase for each elementary box a sequence of O((τ / log n) log log n) words whose entries indicate the open doors along its upper boundary. From this, we can obtain in O((τ / log n) log log n) time a sequence of O((τ / log n)) log log n) words whose entries represent the edges of H , where each entry stores the component identifiers of the two endpoints of the edge. We are now ready to implement the algorithm of Hirschberg, Chandra, and Sarwate. The main steps of the algorithm are as follows, see Fig. 15. Step 1: Find for each vertex of H the neighbor with the smallest and with the largest identifier. This can be done in O((τ / log n) log log2 n) time by sorting the edge lists twice, once in lexicographic order and once in reverse lexicographic order of the identifiers of the endpoints. From these sorted lists, we can extract the desired information in the claimed time, using appropriate word operations. Step 2: Let VH be the vertices of H with at least one neighbor, and let VH ⊆ VH be the vertices v ∈ VH having a neighbor with a smaller identifier than the identifier of v. If |VH | ≥ |VH |/2, we set the successor of each v ∈ VH to the neighbor of v with the smallest identifier. Otherwise, at least half of the nodes in VH have a neighbor with a larger identifier. In this case, we let VH be the set of these nodes, and we set the successor of each v ∈ VH to the neighbor of v with the largest identifier. The successor relation defines a directed forest F on VH such that at least half of the vertices in VH are not a root in F . Given the information available from Step 1 and appropriate word operations, this step can carried out in O((τ / log n) log log n) time. Step 3: Use pointer jumping to determine for each vertex v ∈ VH the identifier of the root of the tree in F that contains v. For this, we set the successor of each v ∈ VH that does not yet have a successor to v itself. Then, for log |VH | = O(log log n) rounds, we set simultaneously for each v ∈ VH the new successor of v to the old successor of the old successor of v (pointer jumping). Each step at least halves the distance of v to its root in F , so it takes O(log |VH |) rounds until each vertex in VH has found the root of its tree in F . Each round can be implemented in O((τ / log n) log log2 n) time by sorting the vertices according to their successors. Thus, this step takes O((τ / log n) log log3 n) time in total. Step 4: Contract each tree of F into a single vertex whose identifier is the smallest identifier in the tree.. Maintain for the original vertices of H a list that gives the identifier of the contracted node that represents it. Again, this step can be carried out in O(τ / log n) log log2 n) time using sorting. After Steps 1–4, the number of non-singleton components in H has at least halved. Thus, by repeating the steps O(log |VH |) = O(log log n) times, we can identify the connected components of H . The total time of the algorithm is O((τ / log n) log log4 n), as claimed. Given the connected components of H , we can find the desired sequence Wi+1 for the boundary Bi+1. Indeed, the procedure from Claim 10.8 outputs a sequence of O((τ / log n) log log n) packed words that gives for each vertex in H an identifier of the component in H that contains it. We can use this list as a lookup table to update the identifiers of the components in the connectivity list of Bi+1. This takes O((τ + ρ)/ log n) log log2 n) time, using sorting. In summary, since we consider log n elementary boxes, the total time for the first phase if O((τ + ρ) log log4 n). The second phase is much easier. For i = logn, . . . , 2, we propagate the connectivity information from Bi+1 to Bi . For this, we need to update the identifiers of the connected components for the cells on the upper boundary of Bi using the identifiers of the connected components on the lower boundary of Bi+1, and then adjust the connectivity list Bi+1 with these new indices. Again, this takes O((τ + ρ)/ log n) log log2 n) time, using sorting. To actually compute the weak Fréchet distance, we use a simplified version of the procedure from Sect. 7. In particular, for the weak Fréchet distance, there are only critical values of the type vertex–vertex and vertex–edge, i.e., there are only O(n2) critical values. However, we aim for a subquadratic running time, so we need to perform the sampling procedure in a slightly different way. Theorem 10.9 Suppose we can answer the decision problem for the weak Fréchet distance in time T (n), for input curves P and Q with n edges each. Then, we can compute the weak Fréchet distance of P and Q in expected time O(n3/2 logc n + T (n) log n), for some fixed constant c > 0. Proof First, we sample a set S of K = 6n1/2 critical values uniformly at random. Then, we find a, b such that the weak Fréchet distance lies between a and b and such that the interval [a, b] contains no other element from S. This takes O(K + T (n) log n) time, using median finding. Similarly to Har-Peled and Raichel [49, Lem. 6.2], we see that the probability, that the interval [a, b] has more than 2γ n3/2 ln n critical values is at most 1/nγ . Indeed, there are at most n2 vertex–vertex and at most 2n2 vertex–edge events. Thus, the total number of critical values is at most L ≤ 3n2. Let U + be the next γ n3/2 ln n larger critical values after dwF( P, Q). Then, the probability that S contains no value from U + is L−|U +| K L K K −1 i=0 ≤ 1 − ≤ exp Analogously, the probability that S contains none of the next γ n3/2 ln n smaller critical values is also at most 1/2nγ , so the claim follows. For γ > 0 large enough, the contribution of this event to the expected running time is negligible. Next, we find all critical values in the interval [a, b]. For this, we must determine all vertex–vertex and vertex–edge pairs with distance in [a, b]. We report for every vertex p of P or Q the set of vertices of the other curve that lie in the annulus with radii a and b around p. Furthermore, we report for every edge e of P or Q the set of vertices of the other curve that lie in the stadium with radii a and b. This can be done efficiently with a range-searching structure for semi-algebraic sets. Agarwal, Matoušek and Sharir [3] show that we can preprocess the vertices of P and the vertices of Q into a data structure that can answer our desired range reporting queries in time O(n1/2 logc n + k), where k is the output size and c > 0 is some fixed constant. The expected preprocessing time is O(n1+ε), where ε > 0 can be made arbitrarily small. We perform O(n) queries, and the total expected output size of our queries is O(n3/2 ln n), so it takes expected time O(n3/2 logc n) to find the critical values in [a, b]. Finally, we perform a binary search on theses critical values to compute the weak Fréchet distance. This takes O(T (n) log n) time. The following theorem summarizes our results on the weak Fréchet distance. Theorem 10.10 The weak Fréchet distance of two polygonal curves, each with n edges, can be computed by a randomized algorithm in time O (n2α(n) log log n) on a pointer machine and in time O ((n2/ log n)(log log n)5) on a word RAM. We have broken the long-standing quadratic upper bound for the decision version of the Fréchet problem. Moreover, we have shown that this problem has an algebraic decision tree of depth O (n2−ε), for some ε > 0 and where n is the number of vertices of the polygonal curves. We have shown how our faster algorithm for the decision version can be used for a faster algorithm to compute the Fréchet distance. If we allow constant-time table-lookup, we obtain a running time in close reach of O (n2). This leaves us with intriguing open research questions. Can we devise a quadratic or even a slightly subquadratic algorithm for the optimization version? Can we devise such an algorithm on the word RAM, that is, with constant-time table-lookup? What can be said about approximation algorithms? Acknowledgements We would like to thank Natan Rubin for pointing out to us that Fréchet-related papers require a witty title involving a dog, Bettina Speckmann for inspiring discussions during the early stages of this research, Günter Rote for providing us with his Ipelet for the free-space diagram, and Otfried Cheong for Ipe. We would also like to thank anonymous referees for suggesting simpler proofs of Lemmas 4.1 and 7.1, for pointing out [53] for Theorem 9.1, and for suggesting a section on the weak Fréchet distance. K. Buchin and M. Buchin in part supported by European Cooperation in Science and Technology (COST) ICT Action IC0903 MOVE. M. Buchin supported in part by the Netherlands Organisation for Scientific Research (NWO) under Project No. 612.001.106. W. Meulemans supported by the Netherlands Organisation for Scientific Research (NWO) under Project No. 639.022.707. W. Mulzer supported in part by DFG Projects MU/3501/1 and MU/3501/2. Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Appendix: Computational Models The standard machine model in computational geometry is the real RAM [56]. Here, data is represented as a (countably) infinite sequence of storage cells. These cells can be of two different types: they can store real numbers or integers. The model supports standard operations on these numbers in constant time, including addition, multiplication, and elementary functions like square-root, sine or cosine. Furthermore, the integers can be used as indices to memory locations. Integers can be converted to real numbers in constant time, but we need to be careful about the reverse direction. The floor function can be used to truncate a real number to an integer, but if we were allowed to use it arbitrarily, the real RAM could solve PSPACE-complete problems The word RAM is essentially a real RAM without support for real numbers. However, on a real RAM, the integers are usually treated as atomic, whereas the word RAM allows for powerful bit-manipulation tricks [42]. More precisely, the word RAM represents the data as a sequence of w-bit words, where w = (log n). Data can be accessed arbitrarily, and standard operations, such as Boolean operations (and, xor, shl, . . .), addition, or multiplication take constant time. There are many variants of the word RAM, depending on precisely which instructions are supported in constant time [27, 42]. The general consensus seems to be that any function in AC0 is acceptable.7 However, it is always preferable to rely on a set of operations as small, and as non-exotic, as possible. Note that multiplication is not in AC0 [43], but nevertheless is often included in the word RAM instruction set [42]. The pointer machine model disallows the use of constant time table lookup, and is therefore a restriction of the (real) RAM model [59, 61]. The data structure is modeled as a directed graph G with bounded out-degree. Each node in G represents a record, with a bounded number of pointers to other records and a bounded number of (real or integer) data items. The algorithm can access data only by following pointers from the inputs (and a bounded number of global entry records); random access is not possible. The data can be manipulated through the usual real RAM operations, but without support for the floor function, for reasons mentioned above. Algebraic Computation Tree Algebraic computation trees (ACTs) [11] are the computational geometry analogue of binary decision trees, and like these they are mainly used for proving lower bounds. Let x1, . . . , xn ∈ R be the inputs. An ACT is a binary tree with two kinds of nodes: computation and branch nodes. A computation node v has one child and is labeled with an expression of the type yv = yu ⊕ yw, where ⊕ ∈ {+, −, ∗, /, √·} is a operator and yu , yw is either an input variable x1, . . . , xn or corresponds to a computation node that is an ancestor of v. A branch node has degree 2 and is labeled by yu = 0 or yu > 0, where again yu is either an input or a variable for an ancestor. A family of algebraic computation trees (Tn )n∈N solves a computational problem (like Delaunay 7 AC0 is the class of all functions f : {0, 1}∗ → {0, 1}∗ that can be computed by a family of circuits (Cn )n∈N with the following properties [11]: (i) each Cn has n inputs; (ii) there exist constants a, b, such that Cn has at most anb gates, for n ∈ N; (iii) there is a constant d such that for all n the length of the longest path from an input to an output in Cn is at most d (i.e., the circuit family has bounded depth); (iv) each gate has an arbitrary number of incoming edges (i.e., the fan-in is unbounded). triangulation or convex hull computation), if for each n ∈ N, the tree Tn accepts inputs of size n, and if for any such input x1, . . . , xn the corresponding path in Tn (where the children of the branch nodes are determined according the conditions they represent) constitutes a computation that represents the answer in the variables yv encountered during the path. Algebraic decision trees are defined as follows: we allow only branch nodes. Each branch node is labeled with a predicate of the form p(x1, . . . , xn ) = 0 or p(x1, . . . , xn ) > 0. The leaves are labeled yes or no. Fix some r ∈ {1, . . . , n}. If p is restricted to be of the form p(x1, . . . , xn ) = in=1 ai xi − b, with at most r coefficients ai = 0, we call the decision tree r -linear. Erickson [40] showed that any 3-linear decision tree for 3SUM has depth (n2). However, Grønlund and Pettie showed that there is a 4-linear decision tree of depth O (n3/2√log n) for the problem. In geometric problems, linear predicates are often much too restrictive. For example, there is no r -linear decision tree for the Fréchet problem, no matter the choice of r : with r -linear decision trees, we cannot even decide whether two given points p and q have Euclidean distance at most 1. 1. Agarwal , P.K. , Ben Avraham , R., Kaplan , H. , Sharir , M. : Computing the discrete Fréchet distance in subquadratic time . SIAM J. Comput . 43 ( 2 ), 429 - 449 ( 2014 ) 2. Agarwal , P.K. , Har-Peled , S. , Mustafa , N.H. , Wang , Y. : Near-linear time approximation algorithms for curve simplification . Algorithmica 42 ( 3-4 ), 203 - 219 ( 2005 ) 3. Agarwal , P.K. , Matoušek , J. , Sharir , M. : On range searching with semialgebraic sets . II. SIAM J. Comput . 42 ( 6 ), 2039 - 2062 ( 2013 ) 4. Ailon , N. , Chazelle , B. : Lower bounds for linear degeneracy testing . J. ACM 52 ( 2 ), 157 - 171 ( 2005 ) 5. Albers , S. , Hagerup, T.: Improved parallel integer sorting without concurrent writing . Inf. Comput . 136 ( 1 ), 25 - 51 ( 1997 ) 6. Alt , H. : The computational geometry of comparing shapes . In: Albers, S. , Alt , H. , Näher , S. (eds.) Efficient Algorithms. Lecture Notes in Computer Science , vol. 5760 , pp. 235 - 248 . Springer, Berlin ( 2009 ) 7. Alt , H. , Buchin , M. : Can we compute the similarity between surfaces? Discrete Comput . Geom. 43 ( 1 ), 78 - 99 ( 2010 ) 8. Alt , H. , Godau , M. : Computing the Fréchet distance between two polygonal curves . Int. J. Comput. Geom. Appl . 5 ( 1-2 ), 75 - 91 ( 1995 ) 9. Alt , H. , Knauer , C. , Wenk , C. : Comparison of distance measures for planar curves . Algorithmica 38 ( 1 ), 45 - 58 ( 2003 ) 10. Aronov , B. , Har-Peled , S. , Knauer , C. , Wang , Y. , Wenk , C. : Fréchet distances for curves, revisited . In: Azar, Y. , Erlebach, T. (eds.) Algorithms-ESA 2006. Lecture Notes in Computer Science , vol. 4168 , pp. 52 - 63 . Springer, Berlin ( 2006 ) 11. Arora , S. , Barak , B. : Computational Complexity . Cambridge University Press , Cambridge ( 2009 ) 12. Baran , I., Demaine , E.D. , Pa˘tras¸cu, M.: Subquadratic algorithms for 3SUM . Algorithmica 50 ( 4 ), 584 - 596 ( 2008 ) 13. Bellman , R. , Kalaba , R.: On adaptive control processes . IRE Trans. Autom. Control 4 ( 2 ), 1 - 9 ( 1959 ) 14. Ben Avraham , R., Filtser , O. , Kaplan , H. , Katz , M.J. , Sharir , M. : He discrete and semicontinuous Fréchet distance with shortcuts via approximate distance counting and selection . ACM Trans. Algorithms 11 ( 4 ), 29 ( 2015 ) 15. Brakatsoulas , S. , Pfoser , D. , Salas , R. , Wenk , C.: On map-matching vehicle tracking data . In: Böhm, K. , et al. (eds.) Proceedings of the 31st VLDB Conference , pp. 853 - 864 . ACM, New York ( 2005 ) 16. Bremner , D. , Chan, T.M., Demaine , E.D. , Erickson , J. , Hurtado , F. , Iacono , J. , Langerman , S. , Pa˘tras¸cu, M. , Taslakian , P. : Necklaces, convolutions, and X + Y . Algorithmica 69 ( 2 ), 294 - 314 ( 2012 ) 17. Bringmann , K. : Why walking the dog takes time: Fréchet distance has no strongly subquadratic algorithms unless SETH fails . In: 55th Annual IEEE Symposium on Foundations of Computer Science-FOCS 2014 , pp. 661 - 670 . IEEE Computer Society, Los Alamitos ( 2014 ) 18. Bringmann , K. , Mulzer , W.: Approximability of the discrete Fréchet distance . J. Comput. Geom . 7 ( 2 ), 46 - 76 ( 2016 ) 19. Buchin , K. , Buchin , M. , Gudmundsson , J. : Constrained free space diagrams: a tool for trajectory analysis . Int. J. GIS 24 ( 7 ), 1101 - 1125 ( 2010 ) 20. Buchin , K. , Buchin , M. , Gudmundsson , J. , Löffler , M. , Luo , J. : Detecting commuting patterns by clustering subtrajectories . Int. J. Comput. Geom. Appl . 21 ( 3 ), 253 - 282 ( 2011 ) 21. Buchin , K. , Buchin , M. , Knauer , C. , Rote , G. , Wenk , C. : How difficult is it to walk the dog ? In: Aichholzer, O. , Hackl , T. (eds.) 23rd EuroCG/FWCG , pp. 170 - 173 . Technischen Universität Graz, Graz ( 2007 ) 22. Buchin , K. , Buchin , M. , Meulemans , W. , Speckmann , B. : Locally correct Fréchet matchings . In: Epstein, L. , Ferragina , P. (eds.) Algorithms-ESA 2012. Lecture Notes in Computer Science , vol. 7501 , pp. 229 - 240 . Springer, Heidelberg ( 2012 ) 23. Buchin , K. , Buchin , M. , Schulz , A. : Fréchet distance of surfaces: some simple hard cases . In: de Berg, M. , Meyer , U. (eds.) Algorithms-ESA 2010. Part II. Lecture Notes in Computer Science , vol. 6347 , pp. 63 - 74 . Springer, Berlin ( 2010 ) 24. Buchin , K. , Buchin , M. , van Leusden , R. , Meulemans , W. , Mulzer , W.: Computing the Fréchet distance with a retractable leash . Discrete Comput. Geom . 56 ( 2 ), 315 - 336 ( 2016 ) 25. Buchin , K. , Buchin , M. , Wang , Y. : Exact algorithms for partial curve matching via the Fréchet distance . In: Proceedings of the Twentieth Annual ACM-SIAM Symposium on Discrete Algorithms , pp. 645 - 654 . SIAM, Philadelphia ( 2009 ) 26. Buchin , K. , Buchin , M. , Wenk , C. : Computing the Fréchet distance between simple polygons . Comput. Geom . 41 ( 1-2 ), 2 - 20 ( 2008 ) 27. Buchin , K. , Mulzer , W.: Delaunay triangulations in O(sort(n)) time and more . J. ACM 58 ( 2 ), 6 ( 2011 ) 28. Buchin , M. : On the Computability of the Fréchet Distance Between Triangulated Surfaces . PhD thesis , Free University Berlin, Berlin ( 2007 ). http://www.diss.fu-berlin. de/diss/receive/FUDISS_thesis_ 000000002618 29. Chambers , E.W. , Colin de Verdière , É. , Erickson , J. , Lazard , S. , Lazarus , F. , Thite , S. : Homotopic Fréchet distance between curves or, walking your dog in the woods in polynomial time . Comput. Geom . 43 ( 3 ), 295 - 311 ( 2010 ) 30. Chan, T.M.: All-pairs shortest paths with real weights in O(n3/ log n) time . Algorithmica 50 ( 2 ), 236 - 243 ( 2008 ) 31. Chan, T.M.: More algorithms for all-pairs shortest paths in weighted graphs . SIAM J. Comput . 39 ( 5 ), 2075 - 2089 ( 2010 ) 32. Chazelle , B.M. , Lee , D.T.: On a circle placement problem . Computing 36 ( 1-2 ), 1 - 16 ( 1986 ) 33. Cook IV , A.F., Driemel , A. , Har-Peled , S. , Sherette , J. , Wenk , C. : Computing the Fréchet distance between folded polygons . In: Dehne, F. , Iacono , J. , Sack , J.-R . (eds.) Algorithms and Data Structures . Lecture Notes in Computer Science , vol. 6844 , pp. 267 - 278 . Springer, Heidelberg ( 2011 ) 34. Cook IV , A.F., Wenk , C. : Geodesic Fréchet distance inside a simple polygon . ACM Trans. Algorithms 7 ( 1 ), 9 ( 2010 ) 35. de Berg, M. , Cook IV , A.F., Gudmundsson , J. : Fast Fréchet queries . Comput. Geom . 46 ( 6 ), 747 - 755 ( 2013 ) 36. Driemel , A. , Har-Peled , S. : Jaywalking your dog: computing the Fréchet distance with shortcuts . SIAM J. Comput . 42 ( 5 ), 1830 - 1866 ( 2013 ) 37. Driemel , A. , Har-Peled , S. , Wenk , C. : Approximating the Fréchet distance for realistic curves in near linear time . Discrete Comput. Geom . 48 ( 1 ), 94 - 127 ( 2012 ) 38. Efrat , A. , Guibas , L.J. , Har-Peled , S. , Mitchell, J.S.B., Murali, T.M.: New similarity measures between polylines with applications to morphing and polygon sweeping . Discrete Comput. Geom . 28 ( 4 ), 535 - 569 ( 2002 ) 39. Eiter , T. , Mannila , H. : Computing Discrete Fréchet Distance . Technical report CD-TR 94 /65, Christian Doppler Laboratory ( 1994 ) 40. Erickson , J. : Bounds for linear satisfiability problems . Chic. J. Theor. Comput. Sci . 1999 , 8 ( 1999 ) 41. Fredman , M.L. : How good is the information theory bound in sorting? Theor . Comput. Sci. 1 ( 4 ), 355 - 361 ( 1975 /76) 42. Fredman , M.L. , Willard , D.E. : Surpassing the information-theoretic bound with fusion trees . J. Comput. Syst. Sci . 47 ( 3 ), 424 - 436 ( 1993 ) 43. Fredman , M.L. , Saxe , J.B. , Sipser , M. : Parity, circuits, and the polynomial-time hierarchy . Math. Syst. Theory 17 ( 1 ), 13 - 27 ( 1984 ) 44. Gajentaan , A. , Overmars , M.H.: On a class of O(n2) problems in computational geometry . Comput. Geom . 5 ( 3 ), 165 - 185 ( 1995 ) 45. Godau , M. : A natural metric for curves-computing the distance for polygonal chains and approximation algorithms . In: Choffrut, C. , Jantzen , M. (eds.) STACS 91. Lecture Notes in Computer Science , vol. 480 , pp. 127 - 136 . Springer, Berlin ( 1991 ) 46. Godau , M. : On the Complexity of Measuring the Similarity Between Geometric Objects in Higher Dimensions . PhD thesis , Free University Berlin, Berlin ( 1998 ) 47. Gudmundsson , J. , Wolle , T.: Towards automated football analysis: algorithms and data structures . In: Proceedings of the 10th Australasian Conference on Mathematics and Computers in Sport 2010 . ANZIAM, Darwin ( 2010 ) 48. Har-Peled , S. , Nayyeri , A. , Salavatipour , M. , Sidiropoulos , A. : How to walk your dog in the mountains with no magic leash . In: Dey, T. , Whitesides , S. (eds.) Computational Geometry (SCG'12) , pp. 121 - 130 . ACM, New York ( 2012 ) 49. Har-Peled , S. , Raichel , B. : The Fréchet distance revisited and extended . ACM Trans. Algorithms 10 ( 1 ), 3 ( 2014 ) 50. Hirschberg , D.S. , Chandra , A.K. , Sarwate , D.V. : Computing connected components on parallel computers . Commun. ACM 22 ( 8 ), 461 - 464 ( 1979 ) 51. Indyk , P. : Approximate nearest neighbor algorithms for Frechet distance via product metrics . In: Computational Geometry (SCG'02) , pp. 102 - 106 . ACM, New York ( 2002 ) 52. Grønlund , A. , Pettie , S. : Threesomes, degenerates, and love triangles . In: 55th Annual IEEE Symposium on Foundations of Computer Science-FOCS 2014 , pp. 621 - 630 . IEEE Computer Society, Los Alamitos ( 2014 ) 53. Katz , M.J. , Sharir , M. : An expander-based approach to geometric optimization . SIAM J. Comput . 26 ( 5 ), 1384 - 1408 ( 1997 ) 54. Maheshwari , A. , Sack , J.-R. , Shahbaz , K. , Zarrabi-Zadeh , H. : Fréchet distance with speed limits . Comput. Geom . 44 ( 2 ), 110 - 120 ( 2011 ) 55. Maheshwari , A. , Sack , J.-R. , Shahbaz , K. , Zarrabi-Zadeh , H. : Improved algorithms for partial curve matching . In: Demetrescu, C. , Halldórsson , M.M. (eds.) Algorithms-ESA 2011. Lecture Notes in Computer Science , vol. 6942 , pp. 518 - 529 . Springer, Heidelberg ( 2011 ) 56. Preparata , F.P. , Shamos , M.I.: Computational Geometry . Texts and Monographs in Computer Science. Springer, New York ( 1985 ) 57. Pa˘tras¸cu, M.: Towards polynomial lower bounds for dynamic problems . In: STOC'10 , pp. 603 - 610 . ACM, New York ( 2010 ) 58. Schönhage , A. : On the power of random access machines . In: Maurer, H.A. (ed.) Automata, Languages and Programming. Lecture Notes in Computer Science , vol. 71 , pp. 520 - 529 . Springer, Berlin ( 1979 ) 59. Schönhage , A. : Storage modification machines . SIAM J. Comput . 9 ( 3 ), 490 - 508 ( 1980 ) 60. Sharir , M. , Agarwal , P.K. : Davenport-Schinzel Sequences and Their Geometric Applications . Cambridge University Press, Cambridge ( 1995 ) 61. Tarjan , R.E.: Efficiency of a good but not linear set union algorithm . J. Assoc. Comput. Mach . 22 ( 2 ), 215 - 225 ( 1975 ) 62. Thorup , M. : Randomized sorting in O(n log log n) time and linear space using addition, shift, and bit-wise Boolean operations . J. Algorithms 42 ( 2 ), 205 - 230 ( 2002 ) 63. Wenk , C. , Salas , R. , Pfoser , D. : Addressing the need for map-matching speed: localizing global curvematching algorithms . In: Proceedings of the 18th International Conference on Scientific and Statistical Database Management , pp. 379 - 388 . IEEE Computer Society, Los Alamitos ( 2006 )


This is a preview of a remote PDF: https://link.springer.com/content/pdf/10.1007%2Fs00454-017-9878-7.pdf

Kevin Buchin, Maike Buchin, Wouter Meulemans, Wolfgang Mulzer. Four Soviets Walk the Dog: Improved Bounds for Computing the Fréchet Distance, Discrete & Computational Geometry, 2017, 180-216, DOI: 10.1007/s00454-017-9878-7