DECLUSTER OUTPUT 1

Copyright \begin{matrix}c \\ ◯ \end{matrix}1994, 1995, 1996, 1997, 1998 David Neto This library is free software; you can redistribute it and/or modify it under the terms of the GNU Library General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public License for more details. You should have received a copy of the GNU Library General Public License along with this library; if not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. You may contact David Neto via email at netod@acm.org, or with greater latency at

Department of Computer Science University of Toronto 10 King’s College Rd. Toronto, Ontario M5S 3G4 Canada

Log : decluster.w, v Revision 1.62 1998/10/16 20:41:41 neto Export the edge comparison function too. Revision 1.61 1998/10/10 21:27:40 neto Updated the test code to use newer name for NONE edge type. Revision 1.60 1998/09/18 17:32:01 neto Added references to Bentley et al work in geometric MSTs. Revision 1.59 1998/07/16 21:58:55 neto Added the LGPL notice in each file. Revision 1.58 1998/06/19 16:23:39 neto Fixed call to pq create (to pq create size) Revision 1.57 1998/06/19 15:01:50 neto Used pq instead of dict for priority queues. Revision 1.56 1998/05/23 16:37:26 neto Factor out the plain MST algorithm (MST) so it can be used with custom cost functions, e.g. by Held-Karp lower bounding process in module ascend. Revision 1.55 1997/12/17 20:01:11 neto Removed const qualifier on tree returned by decluster topology tree. Revision 1.54 1997/12/17 17:07:13 neto Changed compile-time DECLUSTER CONSERVE MEMORY to runtime variable decluster discard topology tree. Fixed a date from 1998 (duh!) to 1988. Revision 1.53 1997/12/13 21:30:07 neto Fixed a copyo. Revision 1.52 1997/12/13 21:28:15 neto decluster mst now takes the tsp instance argument instead of relying on a global variable. Revision 1.51 1997/12/06 20:06:05 neto Be a little defensive when freeing the tree. Revision 1.50 1997/10/18 14:49:49 neto Now use E2 supports() to know when we can use a kd-tree. Revision 1.49 1997/10/17 21:47:55 neto Changed coord2d structure. Revision 1.48 1997/09/27 18:05:55 neto Fixed RCS log behaviour. Revision 1.47 1997/09/05 20:11:31 neto Change augmenting chain to augmenting cycle. improve the comments in the introduction a little bit. Revision 1.46 1997/08/15 20:18:25 neto Added Index major section. Revision 1.45 1997/06/26 19:26:14 neto Better documentation about declustering motivation. Revision 1.44 1997/06/17 13:33:19 neto Make DECLUSTER DEBUG a zero/non-zero test instead of a defined/undefined test. Revision 1.43 1997/06/16 20:08:24 neto Fixed a TeX bug. Revision 1.42 1997/06/13 17:59:32 neto Added support for DECLUSTER CONSERVE MEMORY. If defined, this frees T prime after preprocessing. Revision 1.41 1997/06/13 15:44:49 neto Comment on the bit trick I just made more safe. Revision 1.40 1997/06/13 15:16:27 neto (t + 1) >> 1 is unsafe when t = 2^{31} − 1 because of possible sign extension on the right shift. Replace it with t ⊕ (t >> 1). (Here t is a low bits mask.) Also, this might be faster because it replaces an AC^{0} operation (t+1) with an NC^{0} operation (exclusive-or). Revision 1.39 1997/06/12 19:55:19 neto Made declustertest a little quieter, and nicer for usual verbose==100 output. Revision 1.38 1997/06/12 19:45:04 neto Fixed bug in lca: last step was comparing inlabel numbers of xpp, ypp instead of level numbers. A transcripting error, pretty much. :) The LCA passes the tests on both the simple matrix and Euclidean instances. I’ll encapsulate the debugging output now. Revision 1.37 1997/06/12 18:50:21 neto Testing on one example shows that inlabel, ascendant, level, parent of head, head, parent entries are computed correctly. Revision 1.36 1997/06/11 20:04:49 neto Fixed a depth-first search bookkeeping bug. Revision 1.35 1997/06/11 19:37:37 neto Fixed visibility and ordering of prototypes for print tree. Revision 1.34 1997/06/11 19:31:44 neto Refined when and how print tree is called. Revision 1.33 1997/06/11 19:19:06 neto Must preprocess the MST! In preprocessing: component must be valid for n+n-1 entries. : must make component[w]=NO PARENT. Added function decluster print tree.

Revision 1.32 1997/06/10 22:06:56 neto MST stuff now works, including test. Added code to test LCA. LCA does not work. :( Revision 1.31 1997/06/10 20:51:47 neto Fixed 2 bugs in plain prim. Fixed 2 bugs in testing code. plain prim: mark cities as members of the component. plain prim: in updating arrays, I set infinity wrong. testing code: setup needs number of cities, not number of edges in MST. Revision 1.30 1997/06/10 20:30:24 neto Finished coding self-test. Plain prim fails. It adds self-edges. Revision 1.29 1997/06/04 22:27:08 neto Fixed a macro placement. Revision 1.28 1997/06/04 22:26:04 neto More code for selftest. Revision 1.27 1997/06/04 21:37:33 neto Begin a small test program. Revision 1.26 1997/06/04 21:06:55 neto Now it compiles. Revision 1.25 1997/06/03 22:40:38 neto Removed a bad trailing comma in pow 2. Revision 1.24 1997/06/03 22:35:46 neto Added a comment about whether we should actually store the cost in the topology tree. Revision 1.23 1997/06/03 22:32:39 neto Moved the bit twiddling forward so that it doesn’t split the digestion code. Completed the cluster distance code. Revision 1.22 1997/06/03 22:13:02 neto We don’t need index rightmost 1 Revision 1.21 1997/06/03 22:01:17 neto We don’t need halfmasks. Revision 1.20 1997/06/03 21:58:40 neto Completed coding xpp and ypp code. Created a macro for copying high 1 bit down. Reduced sectional clutter at the end of the web file. Revision 1.19 1997/06/03 19:00:17 neto Finish coding parent of head of inlabel chain. Revision 1.18 1997/06/03 18:54:04 neto Fixed accesses to edges of T prime. Use ‘parent of head of inlabel chain’ instead of ‘head of inlabel chain’. Revision 1.17 1997/06/03 18:50:14 neto Added code for finding inlabel of z, the LCA of x and y. Revision 1.16 1997/06/02 15:50:47 neto Fixed two bugs. First, head levels weren’t being checked in the initialization phase. Second, I was missing the lowest high bit in LCA code. Revision 1.15 1997/05/28 17:31:24 neto Make make happy so I can make a distribution. Revision 1.14 1997/05/28 17:28:25 neto Started on doing LCA queries.

§1  DECLUSTER  DECLUSTERING 5

Revision 1.13 1997/05/27 23:57:02 neto Finished code to create T prime. Created code to compute level, inlabel, and ascendant numbers. Revision 1.12 1997/05/23 22:46:52 neto Missed an at sign... Revision 1.11 1997/05/23 22:46:14 neto Fill out null sections so I can pack this up and back it up to school. Revision 1.10 1997/05/23 22:37:46 neto Fixed two typos. Revision 1.9 1997/05/23 22:34:55 neto Middle of coding up the transformation from MST T to topology tree T’. I need to think more before I code. I may have to go to a full-blown union-find algorithm, but then again, maybe not. Revision 1.8 1997/05/23 20:35:20 neto Bridges drawn from the priority queue might no longer be bridges. Fixed. Still no compile. Revision 1.7 1997/05/23 20:12:00 neto Added accounting for computing the MST length. Revision 1.6 1997/05/23 19:13:54 neto New major section: Interface and code preamble. Include dict.h and kdtree.h Revision 1.5 1997/05/23 19:08:07 neto Implemented Prim’s algorithm in the geometric case. Still not compiled. Revision 1.4 1997/05/23 18:05:06 neto Edit for readability. Skip over cities that are in the component (in Prim’s algorithm). Revision 1.3 1997/05/22 19:22:05 neto Fixed memory accounting. Fixed fencepost error. Revision 1.2 1997/05/22 19:13:52 neto Completed description of declustering idea. Created the interface Now writing the MST routine. Finished coding the general Prim algorithm. (Perhaps reference Knuth’s MILES SPAN). Revision 1.1 1997/05/21 16:57:50 neto Initial revision 1. Declustering. The Lin-Kernighan algorithm can perform quite badly on clustered instances. This becomes evident when one compares its performance on uniform Euclidean instances with its performance on clustered instances such as instance dsj1000. I described the problem and suggested a remedy for it in my research proposal, “Topology and the traveling salesman problem”, March 17, 1997.

This module implements the function decluster d that computes the cluster distance between any pair of cities. It is the key ingredient in the “efficient declustering” fix that I propose for Lin-Kernighan.

2. To understand the problem and its solution, we need to review the Lin-Kernighan heuristic, whose salient points are as follows. The heuristic searches for a cycle in the underlying graph consisting of 2k edges e_{1}, …, e_{2k} with the following properties:

(i) Edge e_{j} is the pair of cities (t_{j}, t_{j+1}), except that edge e_{2k} is the pair of cities (t_{2k}, t_{1}). (ii) Odd-numbered edges are in the current tour; even-numbered edges are not. Removing the odd

edges from the current tour and adding the even edges results in another tour. The 2k edges therefore form an augmenting cycle in the graph. (CHECK NOMENCLATURE)

(iii) The total weight of the odd edges exceeds the total weight of the even edges. That is,

∑_{j odd} cost(e_{j}) − ∑_{j even} cost(e_{j}) > 0.

Although many ideas go into making the search both extensive and efficient, one of the principal ideas is the cumulative gain criterion. The reasoning goes as follows. Taking any cycle satisfying the above three criteria, we may rotate it by twos (taking e′_{j} = e_{j+2}, and so on) so that we eventually find a cycle in which all the partial sums of weights are positive:

cost(e_{1}) > 0
cost(e_{1}) − cost(e_{2}) > 0
cost(e_{1}) − cost(e_{2}) + cost(e_{3}) > 0
cost(e_{1}) − cost(e_{2}) + cost(e_{3}) − cost(e_{4}) > 0

and so on. In general, for any 1 ≤ s ≤ 2k we have

∑_{j odd, j≤s} cost(e_{j}) − ∑_{j even, j≤s} cost(e_{j}) > 0.

3. The proof is easy and is obvious once one draws a picture. Let b be the least partial sum in the original list, and say it occurs after subtracting cost of edge e_{2j}. Let a be the gain made over the last 2k − 2j edges, i.e., the total sum minus b. Rotate the cycle so that e_{2j} appears last in the list of edges, i.e., e′_{2k} = e_{2j} and so on. All the partial sums will now be positive, because the least partial sum appears at the end of the list (to where e_{2j} has migrated) and is equal to the original total sum, a + b, which we know is positive.

6  DECLUSTERING  DECLUSTER §4

4. Improving cycles t are built incrementally from t_{1}, by first removing edge e_{1} = (t_{1}, t_{2}) to form a Hamiltonian path. From then on, city t_{1} anchors the head of the path while we repeatedly flip the tail of the path. Each flip of the tail adds two cities to the sequence of cities t (an odd city and an even city), and hence adds two edges to the sequence of edges e (an even edge and an odd edge). The cumulative gain criterion is applied only after an odd edge is added.

Actually, instead of checking that the cumulative sum is positive, one checks that the cumulative sum is greater than the gain given by the best tour improvement seen so far in this sequence, denoted by best gain. In symbols, the actual cumulative gain criterion is:

∑_{j odd, j<2i} cost(e_{j}) − ∑_{j even, j≤2i} cost(e_{j}) > best gain.

5. The cumulative gain criterion in the usual Lin-Kernighan algorithm deliberately ignores the cost of immediately closing up the Hamiltonian path to form a tour. That is, cost(t_{2i}, t_{1}) is not removed from the cumulative gain because there might be a way to extend the sequence a bit farther and achieve a closing-up cost much lower than cost(t_{2i}, t_{1}).

But there’s a problem with this. If the tail t_{2i} somehow crosses into a cluster of the graph distant from the cluster containing t_{1}, then we end up with an overly-optimistic cumulative gain. This makes the search go on for much longer than it should. We may never gain enough potential to cross the great divide back into the territory surrounding the anchor city t_{1}.

My declustering idea is to discount the cumulative gain by the cluster distance between t_{2i} and t_{1}. Given a path P from u to v, the bottleneck cost of P is the cost of the most expensive edge in P. Then the cluster distance between u and v is the minimum bottleneck cost over all paths from u to v. We denote this value by d(u, v).

The cumulative gain criterion applied when considering adding t_{2i−1} and t_{2i} becomes:

∑_{j odd, j<2i} cost(e_{j}) − ∑_{j even, j≤2i} cost(e_{j}) − d(t_{2i}, t_{1}) > best gain.

6. The rotation trick shows that using the original cumulative gain criterion doesn’t eliminate improvements that satisfy the three Lin-Kernighan criteria given above. That is, applying the cumulative gain criterion can’t hurt.

However, the declustering technique —discounting the cumulative gain function by the cluster distance between t_{2i} and t_{1}— definitely can hurt. That is, it might eliminate valid improvement sequences. (INSERT MORE DETAILS. Also, see research proposal.)

7. If the underlying cost function is postitive-definite and symmetric, then the cluster distance is positivedefinite, symmetric, and also satisfies the triangle inequality. In short, under such conditions d is a metric. (A function f is positive-definite if f(x, y) ≥ 0 with equality if and only if x = y; function f satisfies the triangle inequality if for any x, y, and z, we have f(x, y) ≤ f(x, z) + f(z, y).)

We’re quite lucky because d may be computed very efficiently. In fact with about O(n log n) preprocessing of the graph, online queries of d may be processed in constant time! This remarkable feat is the job of this module.

8. Several ingredients are required to make the declustering code run fast.

First, the cluster distance between vertices u and v in graph G = (V, E) is also equal to the bottleneck cost along the unique simple path from u to v in a minimum spanning tree T = (V, E_{T}) of G. (The proof is left as an exercise.)

Second, we construct a rooted binary tree T′ = (V_{T′}, E_{T′}) from T by making the vertices of T be the leaves of T′, and the edges of T be the internal nodes of T′: V_{T′} = V ∪ E_{T}. Run Kruskal’s minimum spanning tree algorithm on T. Each subtree of T′ corresponds to a connected component of T formed during the execution of Kruskal’s algorithm on T: adding edge (i, j) during reconstruction of T corresponds to adding internal node (i, j) to T′ with its two children being the current components containing i and j. Because Kruskal’s algorithm adds edges in order of non-decreasing cost, T′ is built bottom-up, with more expensive edges toward the root: the cost of the edge labeling an internal node (i, j) of T′ is as great as the cost of any edge labeling a node in the subtree rooted at (i, j). Finally, the bottleneck cost of the unique path from u to v in T is just the weight of the edge labeling the least common ancestor of u and v in T′.

For any rooted tree R, with O(n log n) time preprocessing and O(n) extra space, online least common ancestor (LCA) queries in R may be answered in constant time (INSERT REFERENCE: Schieber1993). This technique is applied to tree T′ from the previous paragraph. The LCA algorithm I give here runs in constant time assuming a fixed word length. If machine word size is unrestricted and given as parameter w, the algorithm given here runs in O(log w) time. (The main barrier to asymptotically faster code is the problem of writing standard C code to find the most significant 1 bit in a word. See http://www.cs.utoronto.ca/~neto/code/

first1bit.html.)

8  INTERFACE AND CODE PREAMBLE  DECLUSTER §9

9. Interface and code preamble. The outline of this module is as follows: #include #include

⟨ System headers 12 ⟩

⟨ Early module headers 15 ⟩

⟨ Module headers 14 ⟩

⟨ Global variables 36 ⟩

⟨ Module definitions 41 ⟩

⟨ Module types 48 ⟩

⟨ Module variables 17 ⟩

⟨ Module subroutines 44 ⟩

⟨ Subroutines 16 ⟩

const char ∗decluster rcs id = ␣ ␣ ␣ ␣ ␣ ␣ ␣ ␣ ;

10. This module provides the following interface.

Procedures decluster setup and decluster cleanup are the usual intialization and shutdown routines. They allocate the persistent space for the rest of the routines. Procedure decluster setup returns the space into which a minimum spanning tree can be written. Procedure decluster cleanup tree frees that space, given its pointer. It is a separate call because the tree is no longer needed after preprocessing, so it can be freed early. In fact, this module throws away its reference to that storage, so the caller must free it with a call to decluster cleanup tree.

Procedure decluster mst builds a minimum spanning tree for the current graph. It is passed a pointer to a block of space into which it should write the tree. It returns the weight of the tree.

We also export the comparison function decluster edge cmp so that others can build priority queues with the same results as here.

Sometimes we want to compute many minimum spanning trees with our own modified cost function, and without taking advantage of k-d trees. In that case, use decluster mst custom, which takes an n-node tree node array into which it should write the answer, two utility arrays of length n, and the custom cost function. It returns the length of the resulting tree.

Procedure decluster preprocess digests a given tree so that least common ancestor queries can be answered quickly by decluster d . It overwrites the tree it is given. Digestion indeed!

Separating the minimum spanning tree building out into a separate call allows the user to substitute their own minimum spanning tree. This is useful in those cases where special knowledge of the input allows them to build the tree very quickly.

Function decluster d computes cluster distance d efficiently. It uses the data structures built by procedure decluster preprocess .

§11  DECLUSTER  INTERFACE AND CODE PREAMBLE 9

11. Some parts of the interface supports uses outside just computing cluster distances.

Procedure decluster topology tree returns a pointer to the topology tree that is preprocessed to answer decluster distance queries. It contains much digested information about the structure of the minimum spanning tree of the instance.

Ordinarily, the topology tree is discarded immediately after preprocessing. Set variable decluster discard topology tree to zero to preserve it for further use. In particular, if decluster discard topology tree is non-zero, decluster topology tree will return Λ after decluster preprocess is called.

Procedure decluster print tree outputs the structure of an arbitrary valid decluster tree t. ⟨ Exported subroutines 11 ⟩ ≡

decluster tree t ∗decluster setup(int n);

void decluster cleanup tree(decluster tree t ∗T);

void decluster cleanup(void);

length t decluster mst (tsp instance t ∗tsp instance, decluster tree t ∗T);

length t decluster mst custom(decluster tree t ∗T , int ∗from, length t ∗dist , length t(∗cost )(int,

int)); int decluster edge cmp(const void ∗a, const void ∗b); void decluster preprocess (decluster tree t ∗T); length t decluster d (int u, int v); decluster tree t ∗decluster topology tree(void); void decluster print tree(FILE ∗out , decluster tree t const ∗t, const char ∗name);

This code is used in section 13.

12. We will be using many routines from external libraries. The interfaces to those routines are described in the following headers. ⟨ System headers 12 ⟩ ≡ #include #include #include

This code is used in section 9.

13. The exported interface is contained in the decluster.h header file, which has the following form. ⟨ 13 ⟩ ≡

extern const char ∗decluster rcs id ;

⟨ Exported types 19 ⟩

⟨ Exported variables 37 ⟩

⟨ Exported subroutines 11 ⟩

14. To ensure consistency between the interface and the implementation, we include our own header. ⟨ Module headers 14 ⟩ ≡ #include

See also sections 27 and 39. This code is used in section 9.

15. Up front, we know we’ll need interfaces to the error checking and memory allocation modules, the length t type (from length), and the cost function and internal access to the internals of the graph (from

read).

⟨ Early module headers 15 ⟩ ≡ #include #include #include #include

This code is used in section 9.

10  SETUP AND CLEANUP  DECLUSTER §16

16. Setup and cleanup. Procedure decluster setup just allocates the resources required by this module, and sets up some module-level convenience variables. It returns space into which a minimum spanning tree should be written. ⟨ Subroutines 16 ⟩ ≡

decluster tree t ∗decluster setup(int the n)

{ decluster tree t ∗mst ; n = the n; ⟨ Allocate space for a spanning tree mst 20 ⟩ ⟨ Other setup code 31 ⟩ return mst ;

}

See also sections 18, 21, 22, 23, 29, 35, 68, 70, and 71. This code is used in section 9.

17. We should declare n. ⟨ Module variables 17 ⟩ ≡

static int n; /∗ The number of cities. ∗/

See also sections 34, 42, 43, 49, 61, and 91. This code is used in section 9.

18. Deallocation is simple (especially with named sections!). ⟨ Subroutines 16 ⟩ +≡

void decluster cleanup(void)

{ ⟨ Other cleanup code 33 ⟩ n = 0;

§19  DECLUSTER  MINIMUM SPANNING TREES 11

} 19. Minimum spanning trees. We’ll be manipulating spanning trees throughout this module. We begin by defining a spanning tree type. A spanning tree on n vertices has n − 1 edges. Each edge is a pair of vertices. The tree is an array of edges, each with their associated cost. ⟨ Exported types 19 ⟩ ≡

typedef struct { int city [2]; length t cost ;

} decluster edge t; typedef struct {

int n;

decluster edge t ∗edge; /∗ An array of n edges ∗/ } decluster tree t;

This code is used in section 13.

20. Since a tree is just an array of edges, we can allocate it in one swoop. ⟨ Allocate space for a spanning tree mst 20 ⟩ ≡

mst = new of (decluster tree t); mst n = n − 1; mst edge = new arr of (decluster edge t, n − 1);

This code is used in section 16.

21. Freeing the tree is just as easy. ⟨ Subroutines 16 ⟩ +≡

void decluster cleanup tree(decluster tree t ∗T) {

if (T) { /∗ Be a little defensive. ∗/ size t r = (T n) ∗ sizeof (decluster edge t) + sizeof (decluster tree t); T n = 0; free mem(T edge); free mem(T); mem deduct (r);

12  MINIMUM SPANNING TREES  DECLUSTER §22

} } 22. We can choose among different minimum spanning tree algorithms according to our knowledge of the underlying structure of the graph. In particular, if the graph is embedded in Euclidean space then we can use Bentley’s k-d trees for semidynamic point sets, implemented in module kdtree, to speed things up considerably. Otherwise we say the graph is unstructured and we use a general minimum spanning tree algorithm.

Our graphs are complete, i.e., there is an edge between each pair of vertices. Prim’s MST algorithm is efficient for unstructured complete graphs because it can be programmed to examine each edge only once.

Prim’s algorithm begins with an empty tree T_{0} = (V, {}), and designates a seed vertex v_{0}. It builds up a single component from V_{0} = {v_{0}} to the entire vertex set V. At step i it forms the new version of the component V_{i+1} from V_{i} by adding the shortest edge between a vertex v ∈ V_{i} in the current version of the component and a vertex v_{i+1} ∈ V \ V_{i} not in the current version of the component.

Oddly enough, Prim is also good for the Euclidean case because it builds up a single component from nothing to the entire tree. We can filter out many of the edges from consideration by using the k-d tree’s ability to hide cities. We hide cities in the current version of the component, i.e., at step i all vertices in V_{i} are hidden. That way nearest neighbour searches will avoid edges that would create cycles. ⟨ Subroutines 16 ⟩ +≡

length t decluster mst (tsp instance t ∗tsp instance, decluster tree t ∗T)

{ extern int verbose; length t total len = 0; errorif (T n ≠ n − 1, ␣ ␣ ␣ ␣ ␣ ␣ ␣

␣ ␣ ␣ ␣ ␣ , T n, n − 1);

if (E2 supports (tsp instance)) {

⟨ Build an MST using k-d trees 26 ⟩ } else { /∗ Use a plain MST algorithm. ∗/

int ∗from = new arr of (int, n);

length t ∗dist = new arr of (length t, n);

total len = decluster mst custom(T , from, dist , cost ); free mem(from); mem deduct (n ∗ sizeof (int)); free mem(dist ); mem deduct (n ∗ sizeof (length t));

}

return total len; } 23. Let’s implement a plain MST algorithm, one which uses only the given cost function and the knowledge that all possible edges are in the graph. We use a generic version of Prim’s algorithm, outlined in the previous section.

Value from is a pointer to an array of n integers, and dist is a pointer to an array of n length values. If vertex i is not in the component, then from[i] is the vertex in the component that is closest to i, and

dist [i] ≡ cost (from[i], i). Otherwise, from[i] = −1.

Variable short len is the minimum distance between a vertex in the component and a vertex outside the component. Variable short to is the outside vertex closest to the component.

We use city 0 as the seed city. ⟨ Subroutines 16 ⟩ +≡

length t decluster mst custom(decluster tree t ∗T , int ∗from, length t ∗dist , length t(∗cost )(int,

int))

{ int i, next edge, short to, n; /∗ Override file-scope n. ∗/ length t short len, total len = 0; from[0] = −1; errorif (T ≡ Λ ∨ T n < 0, ); n = T n + 1; for (i = 1, short len = , short to = −1; i < n; i^{++}) {

from[i] = 0;

dist [i] = cost (0, i);

if (short len > dist [i]) { short len = dist [i]; short to = i;

} } for (next edge = 0; next edge < n − 1; next edge^{++}) {

⟨ Add the edge to short to 24 ⟩

⟨ Update distance arrays and short to 25 ⟩ } return total len;

}

24. ⟨ Add the edge to short to 24 ⟩ ≡

if (verbose ≥ 1000) printf ( ␣ ␣ ␣ ␣ length t spec ,

from[short to], short to, length t pcast (short len)); T edge[next edge].city [0] = short to; T edge[next edge].city [1] = from[short to]; T edge[next edge].cost = short len; total len += short len; from[short to] = −1;

This code is used in section 23.

25. Now that short to has been added to the component, update arrays from and dist . That is, see if short to is closer to an outside vertex than any of the other inside vertices.

We’ll also be updating short to, so first we copy the value of short to to new inside city . Scanning through all the cities wastes effort because we touch cities that are in the component. We could

save time by collecting these cities in one spot. But that would entail still other bookkeeping. Hopefully this simpler version is faster because of its simplicity.

This section completes the general version of Prim’s algorithm. ⟨ Update distance arrays and short to 25 ⟩ ≡

{ const int new inside city = short to; short len = ; for (i = 1; i < n; i^{++}) {

length t d; if (from[i] ≡ −1) continue; /∗ City i is already in the component. ∗/ d = cost (new inside city , i); if (d < dist [i]) {

dist [i] = d;

from[i] = new inside city ; } if (dist [i] < short len) {

short len = dist [i];

short to = i; }

} }

This code is used in section 23.

26. Now we’re ready to implement Prim’s algorithm for geometric inputs. Specifically, we’ll use Bentley’s k-d trees for semidynamic point sets. Throughout the execution of the algorithm, cities in the component will be hidden, and cities outside the component will be visible. This allows us to eliminate the vast majority of edges from consideration.

(It might be nice to measure this. I suspect that this scheme will work very well for uniform instances, because the shortest bridge drawn from the priority queue will probably be consistently short. The scheme might not work so well for clustered instances.)

(Hmmm. Check J. L. Bentley and J. H. Friedman, Fast Algorithms for Constructing Minimal Spanning Trees in Coordinate Spaces, IEEE Transactions on Computers, C-27:2, pp. 97–105, 1978. Also, J. L. Bentley and J. B. Saxe An Analysis of Two Heuristics for theEuclidean Travelling Salesman Problem, pp. 41–49, in 18th Annual Allerton Conference on Communication, Control, and Computing, October 1980.)

We’ll keep a priority queue bridge pq of edges that have one endpoint in the component (recorded in field city [0]) and the other endpoint outside of the component (recorded in field city [1]). We’ll eventually have an entry for each of n − 1 vertices, so we preallocate an array bridge of n of them. (We don’t know in advance which vertex will be the last to enter the component.)

Initially the priority queue contains just the edge from vertex 0 to its nearest neighbour. I’m assuming that the k-d tree is initialized to have no cities hidden. We clean up after ourselves by unhiding all of them at the end. ⟨ Build an MST using k-d trees 26 ⟩ ≡

{ decluster edge t ∗bridge = new arr of (decluster edge t, n); pq t ∗bridge pq = pq create size(decluster edge cmp, n); int next edge; char ∗is in component = new arr of zero(char, n); errorif (bridge pq ≡ Λ, ␣ ␣ ␣ ␣ ); E2 hide(0); bridge[0].city [0] = 0; bridge[0].city [1] = E2 nn(0); bridge[0].cost = cost (0, bridge[0].city [1]); is in component [0] = 1; pq insert (bridge pq , &bridge[0]); for (next edge = 0; next edge < n − 1; next edge^{++}) {

int in, out ;

decluster edge t ∗short bridge;

while (1) { /∗ Get a valid bridge. ∗/ short bridge = pq delete min(bridge pq ); in = short bridge city [0]; /∗ The city already in the component, and therefore hidden. ∗/ out = short bridge city [1]; /∗ The city possibly outside the component. ∗/ if (¬is in component [out ]) break; bridge[in].city [1] = E2 nn(in); bridge[in].cost = cost (in, bridge[in].city [1]); pq insert (bridge pq , bridge + in);

}

T edge[next edge] = ∗short bridge; total len += short bridge cost ; E2 hide(out ); is in component [out ] = 1; bridge[in].city [1] = E2 nn(in); bridge[in].cost = cost (in, bridge[in].city [1]); pq insert (bridge pq , bridge + in); bridge[out ].city [0] = out ; bridge[out ].city [1] = E2 nn(out ); bridge[out ].cost = cost (out , bridge[out ].city [1]); pq insert (bridge pq , bridge + out );

}

pq destroy (bridge pq );

free mem(is in component ); mem deduct (n ∗ sizeof (char)); free mem(bridge); mem deduct (n ∗ sizeof (decluster edge t)); E2 unhide all ( );

}

This code is used in section 22.

27. We need to import the interfaces to priority queues and to k-d trees. ⟨ Module headers 14 ⟩ +≡ #include #include

28. One last thing we need is the comparison function decluster edge cmp for edges so the priority queue can give edges to us in the right order. Function decluster edge cmp has the same interface as required by the C library function qsort .

We force determinacy by comparing pointer values in case the lengths are the same. (I don’t think this is necessary in the use in Prim’s algorithm above, but it can’t hurt.) ⟨ Edge comparison function 28 ⟩ ≡

int decluster edge cmp(const void ∗a, const void ∗b)

{ length t len diff = ((const decluster edge t ∗) a) cost − ((const decluster edge t ∗) b) cost ; return len diff < 0 ? −1 : (len diff > 0 ? 1 : ((int) (((const decluster edge t ∗) a) − ((const

decluster edge t ∗) b)))); }

This code is used in section 29.

29. I also use this function in the test routines, but it is local to both modules. ⟨ Subroutines 16 ⟩ +≡

§30  DECLUSTER  TOPOLOGY TREES 17

⟨ Edge comparison function 28 ⟩ 30. Topology trees. Now, given a minimum spanning tree T such as the one constructed in the previous sections, we’d like to create the related topology tree T′ as described in the opening sections. When T has n vertices, T′ has 2n − 1 vertices. Each is an entry in an array of decluster edge t nodes. The city array in each node is aliased to the name child , and stores child indices. (T′ is binary.) A value of −1 indicates no child.

The first n entries of T′ correspond to the vertices of the original graph. They have zero cost and no children.

The LCA queries don’t actually need T prime, so we can throw it away if need be. We defend against a memory leak in case we reallocate for a different size.

#define child city

#define (−1) ⟨ Allocate and pre-initialize T prime 30 ⟩ ≡

if (T prime ≡ Λ ∨ T prime n ≠ n + n − 1) { int i; ⟨ Clean up T prime 32 ⟩ /∗ Defend against a memory leak. ∗/ T prime = new of (decluster tree t); T prime n = n + n − 1; T prime edge = new arr of (decluster edge t, n + n − 1); for (i = 0; i < n; i^{++}) {

T prime edge[i].child [0] = ; T prime edge[i].child [1] = ; T prime edge[i].cost = 0;

} }

This code is used in sections 31 and 35.

31. ⟨ Other setup code 31 ⟩ ≡

⟨ Allocate and pre-initialize T prime 30 ⟩

See also sections 50 and 62. This code is used in section 16.

32. ⟨ Clean up T prime 32 ⟩ ≡

if (T prime) { const int n = T prime n; free mem(T prime edge); mem deduct (n ∗ sizeof (decluster edge t)); free mem(T prime); mem deduct (sizeof (decluster tree t)); T prime = Λ;

}

This code is used in sections 30, 33, and 35.

33. ⟨ Other cleanup code 33 ⟩ ≡

⟨ Clean up T prime 32 ⟩

See also sections 51 and 63. This code is used in section 18.

18  TOPOLOGY TREES  DECLUSTER §34

34. We must define T prime. ⟨ Module variables 17 ⟩ +≡

static decluster tree t ∗T prime = Λ;

35. The first n entries of T prime are filled out at memory allocation time. We then copy T into the last n − 1 positions of T prime in increasing order while making the appropriate links.

Tree T prime is digested into an array of four integers: level , inlabel , ascendant , and head.

In the normal case we throw away T prime after preprocessing because it isn’t needed for LCA queries. In this case we may have to reallocate T prime here. ⟨ Subroutines 16 ⟩ +≡

void decluster preprocess (decluster tree t ∗T)

{ errorif (T n ≠ n − 1, ␣ ␣ ␣ ␣ ␣ ␣ , T n, n − 1); if (T prime ≡ Λ) ⟨ Allocate and pre-initialize T prime 30 ⟩ ⟨ Fill the last n − 1 entries of T prime with data from T 38 ⟩ print tree(T prime, ); ⟨ Copy cost fields from T prime 52 ⟩ ⟨ Compute level numbers 53 ⟩ ⟨ Compute inlabel numbers 54 ⟩ ⟨ Compute ascendant numbers 58 ⟩ ⟨ Compute parent of head numbers 60 ⟩ ⟨ Verbose: print the preprocessing data 88 ⟩ if (decluster discard topology tree) {

⟨ Clean up T prime 32 ⟩ }

}

36. First off, in most cases we want to save space by throwing away topology tree T prime once we’re through processing it. But sometimes we want to keep it around. We make it optional through the globally visible variable decluster discard topology tree. ⟨ Global variables 36 ⟩ ≡

int decluster discard topology tree = 1;

This code is used in section 9.

37. ⟨ Exported variables 37 ⟩ ≡

extern int decluster discard topology tree;

This code is used in section 13.

38. We add each edge of T in ascending order. Essentially we simulate Kruskal’s algorithm, but with the knowledge that every edge ends up in the output. Tree T prime in some sense records the history of the Kruskal algorithm computation.

Array component stores the evolving component structure of the topology tree. Entry component [i] points to the root of the component containing vertex i. Entries may go out of date by the addition of new edges, but they can always be updated by following links up the tree. A value of indicates that a vertex is a root. #define (−2) ⟨ Fill the last n − 1 entries of T prime with data from T 38 ⟩ ≡

{ int r, w, i, ∗component = new arr of (int, n + n); /∗ r is a read cursor, w is a write cursor ∗/ for (i = 0; i < n; i^{++}) component [i] = ; sort (T edge, (size t) (n − 1), sizeof (decluster edge t), decluster edge cmp); print tree(T , ); for (r = 0, w = n; r < n − 1; r^{++}, w^{++}) {

T prime edge[w] = T edge[r];

⟨ Link T prime edge[w] with its two children and update component 40 ⟩ } free mem(component ); mem deduct ((n + n) ∗ sizeof (int));

}

This code is used in section 35.

39. We need to import the sorting routine from module lk. The user has the option of changing it. ⟨ Module headers 14 ⟩ +≡ #include

40. Does an amortized argument prove that the following runs in linear time over the entire sequence? It is similar to the classic path compression solution to the merge/find problem that is a standard component in implementations of Kruskal’s algorithm. That has worst-case time that is almost linear, O(nα(n)). But I’m confused as to whether it does less or more work. I’ll have to think about it more. (Some thoughts: Cf. combining tree in asyncrhonous PRAM (Nishimura); linear time heapify.)

We save some buffer space and some time because we already know where the nodes will all be pointing to. They’ll point to the vertex labeled with the edge we’re adding, namely T prime edge[w]. ⟨ Link T prime edge[w] with its two children and update component 40 ⟩ ≡

{ int i, here, parent ; component [w] = ; for (i = 0; i < 2; i^{++}) {

here = T prime edge[w].city [i];

while ((parent = component [here]) ≠ ) { component [here] = w; here = parent ;

}

component [here] = w;

T prime edge[w].child [i] = here; /∗ child [i] is aliased to city [i] ∗/ }

}

This code is used in section 38.

20  BIT TWIDDLING  DECLUSTER §41

41. Bit twiddling. The least common ancestor code below requires some fast bit manipulation. In several places we need to copy the top 1 bit of a word down through all the lower bits. We can do this

very quickly without any branch instructions. ⟨ Module definitions 41 ⟩ ≡ #if ≡ 8 #define copy 1 down(X)

(((X) |= (X) ≫ 1), ((X) |= (X) ≫ 2), ((X) |= (X) ≫ 4), ((X) |= (X) ≫ 8), ((X) |= (X) ≫ 16), ((X) |= (X) ≫ 32))

#else /∗ ¬ ≡ 8. Then hopefully integers are at most 4 bytes wide. ∗/ #define copy 1 down(X)

(((X) |= (X) ≫ 1), ((X) |= (X) ≫ 2), ((X) |= (X) ≫ 4), ((X) |= (X) ≫ 8), ((X) |= (X) ≫ 16))

#endif

See also section 75. This code is used in section 9.

42. Array pow 2 stores powers of two. ⟨ Module variables 17 ⟩ +≡

static const int pow 2 [ ] = { #1,^{#}2,^{#}4,^{#}8,

#10,^{#}20,^{#}40,^{#}80,
#100,^{#}200,^{#}400,^{#}800,
#1000,^{#}2000,^{#}4000,^{#}8000,
#10000,^{#}20000,^{#}40000,^{#}80000,
#100000,^{#}200000,^{#}400000,^{#}800000,
#1000000,^{#}2000000,^{#}4000000,^{#}8000000,
#10000000,^{#}20000000,^{#}40000000,^{#}80000000

#if ≡ 8

,^{#}100000000,^{#}200000000,^{#}400000000,^{#}800000000,
#1000000000,^{#}2000000000,^{#}4000000000,^{#}8000000000,
#10000000000,^{#}20000000000,^{#}40000000000,^{#}80000000000,
#100000000000,^{#}200000000000,^{#}400000000000,^{#}800000000000,
#1000000000000,^{#}2000000000000,^{#}4000000000000,^{#}8000000000000,
#10000000000000,^{#}20000000000000,^{#}40000000000000,^{#}80000000000000,
#100000000000000,^{#}200000000000000,^{#}400000000000000,^{#}800000000000000,
#1000000000000000,^{#}2000000000000000,^{#}4000000000000000,^{#}8000000000000000

#endif };

§43  DECLUSTER  BIT TWIDDLING 21

43. It will be useful to find logarithms of small numbers quickly. ⟨ Module variables 17 ⟩ +≡

static const int floor log 2 small [16] = { /∗ dummy ∗/ 0,

0,

1, 1,

2, 2, 2, 2,

3, 3, 3, 3, 3, 3, 3, 3};

44. The first bit manipulation routine I’ll define is floor log 2 (x), which computes ⌊log_{2}(x)⌋, as you might imagine. Here I assume that the argument x is a positive integer.

We use a binary search on the bits to find the highest 1 bit. Let the grey width be the number of bits of x that might be set to 1, counting from the least significant bit. Initially, the grey width of x is the word size. Each iteration cuts grey width by half. When x is known to be small, we just use table lookup.

The loop is unrolled as the body of a switch statement; we depend on each case falling through to the next. Unrolling loops is as old as the sun (Sun?), but unrolling a binary search in this way is pretty neat. I first saw this “unrolling a binary search” trick in Knuth’s book Literate Programming. It is also reminscent of “Duff’s device”; see the C Language FAQ, found on comp.lang.c.

The table lookup at the end saves a couple of possible branches at the cost of a possible cache miss. Is it worth it?

As for all the shifting, the shifting in the if conditions are constant expressions and so will be precomputed by the compiler. The shifts of variable x are by powers of two, and so benefit from fast barrel shifters in modern processors. ⟨ Module subroutines 44 ⟩ ≡ #if > 8

#error  ␣   ␣         ␣        ␣        ␣  ␣   ␣    ␣  ␣

#endif static inline int floor log 2 (unsigned int x) {

int ans = 0; #if ≡ 8

if (x & (^{#}ffffffff ≪ 32)) ans += 32, x ≫= 32; #endif

if (x & (^{#}ffff ≪ 16)) ans += 16, x ≫= 16; if (x & (^{#}ff ≪ 8)) ans += 8, x ≫= 8; if (x & (^{#}f ≪ 4)) ans += 4, x ≫= 4; ans += floor log 2 small [x]; return ans ;

}

See also section 64. This code is used in section 9.

45. We will be masking out portions of words.

Macro lo mask lets us keep only the bottom bits of a word: w & lo mask (i) is the bottom i bits of w. Macro hi mask lets us discard the bottom bits of a word: w & hi mask (i) is all but the bottom i bits of w.

#define lo mask (X) (pow 2 [(X)] − 1) #define hi mask (X) (∼lo mask (i))

22  LEAST COMMON ANCESTORS  DECLUSTER §46

46. Least common ancestors. Now we’re ready to get to the meat of the module: least common ancestors. In a rooted tree, vertex u is an ancestor of vertex v if u appears on the path from u to the root; we also say that v is a descendant of u. Equivalently, u is an ancestor of v if v appears in the subtree rooted at u. Note that each node is its own ancestor and its own descendant.

Any two vertices in a rooted tree have a non-empty set of common ancestors; in fact, the root is an ancestor of every node. The least common ancestor of u and v, denoted by LCA(u, v), is that ancestor of both u and v which is a descendant of all the ancestors of common to both u and v. In other words, LCA(u, v) is the root of the smallest subtree in which both u and v appear.

47. Now that we’ve fixed our nomenclature, let’s think about algorithms to compute least common ancestors. A simple algorithm to find LCA(u, v) would be to write down the list of vertices on the paths from both u and v to the root and to find their convergence point. In the worst case this would take time proportional to the depth of the tree, which could be very bad.

One might envision a scheme where we could digest the tree into a corresponding balanced tree and run the simple algorithm on the balanced tree. But this would still run in logarithmic time.

It turns out that the predigestion idea can be worked into a scheme so that the LCA queries can be answered in constant time. (“Constant time” assumes a fixed word size. Some people might complain about this. For machine words with w bits, the LCA query algorithm presented here runs in O(log w) time.) The first known algorithm to achieve constant-time queries is presented in D. Harel and R. E. Tarjan, “Fast algorithms for finding nearest common ancestors. SIAM Journal on Computing, 13(2):338–355, May 1984. Here I use a simpler algorithm, a serial version of a parallel algorithm: see B. Schieber and U. Vishkin “On finding lowest common ancestors: simplification and parallelization. SIAM Journal on Computing, 17(6):1253–1262, December 1988, and also Chapter 6, “Parallel Lowest Common Ancestor Computation” in the book The Synthesis of Parallel Algorithms, (INSERT PUBLISHER), 1993, edited by J. Reif. I’ve used this later paper as the basis for this implementation.

48. The tree is digested into three integer arrays indexed by vertex number: level , inlabel , and ascendant . We attempt to reduce cache misses when processing LCA queries by grouping these values for each vertex in a single structure, digest t .

I also include a cost field in the digest structure because I throw T prime away when not debugging. See also the head numbers below.

⟨ Module types 48 ⟩ ≡

typedef struct { int level , inlabel , ascendant ; length t cost ;

} digest t;

This code is used in section 9.

49. ⟨ Module variables 17 ⟩ +≡

§50  DECLUSTER  LEAST COMMON ANCESTORS 23

static digest t ∗digest ; 50. The first n correspond to vertices of the original graph, and therefore have zero cost. We initialize those fields at allocation time instead of processing time because they are a function of n only. ⟨ Other setup code 31 ⟩ +≡

digest = new arr of (digest t, n + n − 1); {

int i; for (i = 0; i < n; i^{++}) {

digest [i].cost = 0; }

}

51. ⟨ Other cleanup code 33 ⟩ +≡

free mem(digest ); mem deduct ((n + n − 1) ∗ sizeof (digest t));

52. ⟨ Copy cost fields from T prime 52 ⟩ ≡

{ int i; for (i = n; i < n + n − 1; i^{++}) digest [i].cost = T prime edge[i].cost ;

}

This code is used in section 35.

53. Entry digest [i].level is just the depth of vertex i. We compute all of level numbers in linear time with a breadth first search. The search begins at the root of the tree, T prime edge[n + n − 2].

Array queue serves as a last-in first-out queue. It is indexed by a read cursor r and a write cursor w. ⟨ Compute level numbers 53 ⟩ ≡

{ int ∗queue = new arr of (int, n + n − 1), r, w, i, ch; digest [n + n − 2].level = 0; queue[0] = n + n − 2; for (r = 0, w = 1; r < w; r^{++}) {

const int here = queue[r], cur level = digest [here].level + 1; for (i = 0; i < 2; i^{++})

if ((ch = T prime edge[here].child [i]) ≠ ) { digest [ch].level = cur level ; queue[w^{++}] = ch;

} } free mem(queue); mem deduct ((n + n − 1) ∗ sizeof (int));

}

This code is used in section 35.

54. The inlabel number of a vertex is a function of its preorder number and the size of subtree rooted there. These are stored in arrays preorder and size, naturally. Value preorder [i] is one more than the number of vertices visited before vertex i in a preorder traversal of the tree. Value size[i] is the number of vertices in the subtree rooted at i. ⟨ Compute inlabel numbers 54 ⟩ ≡

{ int ∗preorder = new arr of (int, n + n − 1), ∗size = new arr of (int, n + n − 1); int preorder number = 0;

#define

⟨ Do a depth-first search of T prime 55 ⟩ #undef

free mem(preorder ); free mem(size); mem deduct (2 ∗ (n + n − 1) ∗ sizeof (int));

}

This code is used in section 35.

55. Both preorder and size statistics are computed via a depth-first search of the tree. I’ll separate the bookkeeping from the real action by using subsections. This way I can use the same CWEB bookkeeping code twice.

Array in stores the roots of the subtrees we’re currently exploring. It behaves as a stack. Entry cur child [i] tells us which child (0 or 1) of vertex in[i] is currently being explored.

We know the graph is acyclic, so we don’t need to mark nodes. Also, the graph is a binary tree. We begin at the root, vertex n + n − 2.

⟨ Do a depth-first search of T prime 55 ⟩ ≡

{ int ∗in = new arr of (int, n + n − 1), ∗cur child = new arr of (int, n + n − 1); int top, here, next child ; top = 0; in[top] = n + n − 2; cur child [top] = −1; while (top ≥ 0) {

here = in[top];

switch (cur child [top]) {

case −1: /∗ Visit myself, then fall through. ∗/

⟨ Visit here the first time 56 ⟩ case 0: /∗ Visit the children. ∗/

next child = T prime edge[here].child [^{++}cur child [top]];

if (next child ≠ ) { /∗ Push that child on the stack. ∗/ in[^{++}top] = next child ; cur child [top] = −1;

}

break; default: ⟨ Visit here the last time 57 ⟩ /∗ … and pop the stack. ∗/

top^{−}; }

} free mem(in); free mem(cur child ); mem deduct (2 ∗ (n + n − 1) ∗ sizeof (int));

}

This code is used in sections 54 and 58.

56. For inlabel numbers, we need to update the preorder number as we enter a node. ⟨ Visit here the first time 56 ⟩ ≡ #if defined ( )

preorder [here] = ^{++}preorder number ; #endif

See also section 59. This code is used in section 55.

57. For inlabel numbers, we need to update the size array as we leave a node.

I’ll define inlabel [here] but not explain its use. The preorder numbers of all the nodes in the subtree of T prime rooted at here range from preorder [here] through preorder [here] + size[here] − 1. In fact, together they form precisely that interval, no more and no less. Value inlabel [here] is the value in that interval having most trailing 0 bits when expressed in binary.

Collectively the inlabel fields form an injective mapping from the N: = n + n − 1 vertices of T prime into a complete binary tree B of depth ⌊log_{2} N⌋ + 1.

See below for how inlabel s are used.

With care I can elide the taking of logarithms, as I’ve done below in the LCA computation. Maybe later. Or maybe never, since this is just an initialization phase, so timing is not so critical. ⟨ Visit here the last time 57 ⟩ ≡ #if defined ( )

{ const int child0 = T prime edge[here].child [0], child1 = T prime edge[here].child [1]; size[here] = 1 + (child0 ≡ ? 0 : size[child0 ]) + (child1 ≡ ? 0 : size[child1 ]);

}

{ const unsigned int last = preorder [here] + size[here] − 1, i = floor log 2 ((preorder [here] − 1) ⊕ last ); digest [here].inlabel = hi mask (i) & last ; /∗ All but the bottom i bits of last . ∗/

} #endif

This code is used in section 55.

58. The ascendant numbers may also be computed during a depth-first search. Again, see Schieber. The inlabel numbers must be computed first. ⟨ Compute ascendant numbers 58 ⟩ ≡ #define

⟨ Do a depth-first search of T prime 55 ⟩ #undef

This code is used in section 35.

59. Assuming x is positive, value (x ⊕ (x − 1)) & x is the rightmost 1 bit of x, but isolated. ⟨ Visit here the first time 56 ⟩ +≡ #if defined ( )

if (top ≡ 0) { /∗ The root. ∗/

digest [here].ascendant = digest [here].inlabel ; } else {

const int p = in[top − 1], ap = digest [p].ascendant , ip = digest [p].inlabel ; /∗ Parent data. ∗/ if (digest [here].inlabel ≡ ip) digest [here].ascendant = ap; else {

const int ih = digest [here].inlabel ;

digest [here].ascendant = ap + ((ih ⊕ (ih − 1)) & ih); }

} #endif 60. Schieber also prescribes that we compute a head table telling us the head of each inlabel chain. That is head [i] is the shallowest vertex in T prime with inlabel i.

We only ever want to know the parent of a head, so instead I compute a parent of head with the following defining equation:

parent of head [inlabel w ] ≡ parent (head [inlabel w ]).

This saves an indirection.

All the parent of head fields can be defined in four easy passes over the array. But we do need the head array during initialization. It is indexed by inlabel number, so it ranges from 1

through n + n. We also compute a parent array, indexed by vertex (in T prime). ⟨ Compute parent of head numbers 60 ⟩ ≡

{ int i, j, ∗head = new arr of (int, n + n), ∗parent = new arr of (int, n + n − 1); for (i = 0; i < n + n; i^{++}) head [i] = −1; /∗ Empty each bin. ∗/ for (i = 0; i < n + n − 1; i^{++}) { /∗ Compute heads of inlabel chains ∗/

const int ii = digest [i].inlabel , hii = head [ii ];

if (hii ≡ −1 ∨ digest [i].level < digest [hii ].level ) head [ii ] = i; } for (i = 0; i < n + n − 1; i^{++}) parent [i] = , parent of head [i] = ;

/∗ Make parents null. ∗/ for (i = 0; i < n + n − 1; i^{++}) /∗ Compute parents. ∗/

for (j = 0; j < 2; j^{++}) { const int child = T prime edge[i].child [j]; if (child ≠ ) parent [child ] = i;

} for (i = 0; i < n + n − 1; i^{++}) /∗ Record parents of heads of chains. ∗/

parent of head [digest [i].inlabel ] = parent [head [digest [i].inlabel ]]; ⟨ Verbose: print parent 89 ⟩ ⟨ Verbose: print head 90 ⟩ free mem(head ); mem deduct ((n + n) ∗ sizeof (int)); free mem(parent ); mem deduct ((n + n − 1) ∗ sizeof (int));

}

This code is used in section 35.

61. We need to define parent of head . ⟨ Module variables 17 ⟩ +≡

int ∗parent of head ;

62. It must be allocated… ⟨ Other setup code 31 ⟩ +≡

parent of head = new arr of (int, n + n);

63…. and freed. ⟨ Other cleanup code 33 ⟩ +≡

free mem(parent of head ); mem deduct ((n + n) ∗ sizeof (int));

28  LEAST COMMON ANCESTOR QUERIES  DECLUSTER §64

64. Least common ancestor queries. We’re now ready to process least common ancestor queries. ⟨ Module subroutines 44 ⟩ +≡

static inline int decluster lca(int x, int y)

{ const digest t xd = digest [x], yd = digest [y]; if (xd .inlabel ≡ yd .inlabel ) return (xd .level ≤ yd .level ) ? x : y; else { /∗ The LCA is neither x nor y. ∗/ const int xil = xd .inlabel , yil = yd .inlabel ; int b, inlabel z , common, jmask , xpp, ypp; ⟨ Set b to the LCA of C_{x} and C_{y} in the complete tree B 65 ⟩ ⟨ Find inlabel z , common, and jmask 66 ⟩ ⟨ Find xpp and ypp 67 ⟩ return (digest [xpp].level ≤ digest [ypp].level ) ? xpp : ypp;

} }

65. The digestion of T prime decomposes it into (contiguous) chains of vertices. Chain C_{u} is the chain containing vertex u of T prime. The inlabel numbers form an injective mapping from these chains of T prime into a complete binary tree B. Tree B is not explicitly stored, it’s existence is implicit in the inlabel numbers. Futhermore, inlabel [u] is the ordinal position of the vertex C_{u} of B in an in-order traversal of B, counting from 1.

Finding least common ancestors in binary trees is easy, given in-order numbers.

For clarity, Schieber describes things in terms of logarithms, but in the interest of speed I implemented the computation without logarithms.

The expression of a 1 followed by many zeroes, namely, t⊕(t ≫ 1), used to be (t+1) ≫ 1. The old code is unsafe when t is the largest positive int, e.g., 2^{31} − 1, because adding 1 to this value makes a large negative number. Then the right shift might extend the sign down and preserve the negative sign bit.

The new code may even be faster because it replaces an AC^{0} operation (addition) and an NC^{0} operation (a single bit shift) with NC^{0} operations (exclusive-or, logical and, and a single bit shift). On today’s machines, who knows? #define min(A, B) ((A) < (B) ? (A) : (B)) #define max (A, B) ((A) > (B) ? (A) : (B)) ⟨ Set b to the LCA of C_{x} and C_{y} in the complete tree B 65 ⟩ ≡

{ const unsigned int xfuzz = xil ⊕ (xil − 1), yfuzz = yil ⊕ (yil − 1); const unsigned int lomask = max (xfuzz , yfuzz ) ≫ 1; int hi diff = (xil ⊕ yil ) & (∼lomask ); if (hi diff ) { /∗ Neither ancestor nor descendant. ∗/

int t = hi diff ;

copy 1 down(t); /∗ Sort of “sign extend” the top 1 bit of t down through the lower bits. ∗/ b = (∼t) & xil ; /∗ Keep only the top few bits common to both xil and yil . ∗/ b |= t ⊕ (t ≫ 1); /∗ But append a 1 followed by many zeros. ∗/

}

else { /∗ One is an ancestor of the other. ∗/

b = (xfuzz ≥ yfuzz ) ? xil : yil ; } ⟨ Verbose: print b 92 ⟩

}

This code is used in section 64.

§66  DECLUSTER  LEAST COMMON ANCESTOR QUERIES 29

66. Let z = LCA(x, y). This section computes inlabel z :=inlabel [z] without directly computing z itself. We use the ascendant numbers and some bit masking. Again Schieber describes things in terms of logarithms, but I have elided the logarithms.

In declarative form, inlabel [z]:=inlabel [x] & ∼lo mask (j) where j = index rightmost 1 (u), u = common & ∼lo mask (i), common = ascendant [x] & ascendant [y], and i = index rightmost 1 (b).

Value jmask is j low-order 1 bits. ⟨ Find inlabel z , common, and jmask 66 ⟩ ≡

{ const int imask = (b ⊕ (b − 1)) ≫ 1; /∗ Mask of bits corresponding to low 0’s in b. ∗/ const int u = (common = xd .ascendant & yd .ascendant ) & ∼imask ; jmask = (u ⊕ (u − 1)) ≫ 1; /∗ Mask of bits corresponding to low 0’s in u. ∗/ inlabel z = (xd .inlabel & ∼jmask ) | (jmask + 1); /∗ yd .inlabel would also do. ∗/ ⟨ Verbose: print inlabel z , common, and jmask 93 ⟩

}

This code is used in section 64.

67. Vertex xpp is either x itself or is the parent of the head of the inlabel chain at inlabel w x . Value inlabel w x is all but the bottom k + 1 bits of inlabel [x], followed by a 1 and k zeroes.

Vertex ypp is similar. ⟨ Find xpp and ypp 67 ⟩ ≡

if (xd .inlabel ≡ inlabel z ) xpp = x; else {

int inlabel w x ; int kmask = xd .ascendant & jmask ; copy 1 down(kmask ); inlabel w x = (∼kmask ) & xil ; /∗ Keep only the top few bits of inlabel [x]. ∗/ inlabel w x |= kmask ⊕ (kmask ≫ 1); /∗ But append a 1 followed by many zeros. ∗/ xpp = parent of head [inlabel w x ]; ⟨ Verbose: print inlabel w x , xpp 94 ⟩

} if (yd .inlabel ≡ inlabel z ) ypp = y; else {

int inlabel w y ; int kmask = yd .ascendant & jmask ; copy 1 down(kmask ); inlabel w y = (∼kmask ) & yil ; /∗ Keep only the top few bits of inlabel [y]. ∗/ inlabel w y |= kmask ⊕ (kmask ≫ 1); /∗ But append a 1 followed by many zeros. ∗/ ypp = parent of head [inlabel w y ]; ⟨ Verbose: print inlabel w y , ypp 95 ⟩

}

This code is used in section 64.

68. Given the least common ancestor function, computing the cluster distance d is very easy. We might consider reducing memory requirements by computing the cost on the fly instead of storing it

in the topology tree T prime. ⟨ Subroutines 16 ⟩ +≡

length t decluster d (int u, int v) {

return digest [decluster lca(u, v)].cost ; }

30  EXTRAS  DECLUSTER §69

69. Extras. Here are the extra support routines. That is, these routines are not required for supporting decluster distance queries.

70. We may need to see what T prime looks like. The caller can safely change its contents if want to, because we’ve already preprocessed it. ⟨ Subroutines 16 ⟩ +≡

decluster tree t ∗decluster topology tree(void) {

return T prime; }

71. It is be useful to be able to output a tree. ⟨ Subroutines 16 ⟩ +≡

void decluster print tree(FILE ∗out , decluster tree t const ∗t, const char ∗name) {

if (t) { int n = t n, i; const char ∗print name = name ? name : ;

errorif (t ≡ Λ,  ␣     ␣ ␣    ␣  );
errorif (n < 0,  ␣    ␣  ␣    ␣  ␣  ␣  , print name, n);

fprintf (out , , print name, t n); for (i = 0; i < n; i^{++}) {

fprintf (out , ␣ ␣ ␣ length t spec , i, t edge[i].city [0], t edge[i].city [1],

length t pcast (t edge[i].cost )); }

} else {

fprintf (out , ␣ ␣ ␣ , name); fprintf (out ,

    ␣    ␣     ␣     ␣    ␣        ␣                               ␣  ␣  );

§72  DECLUSTER  TESTING 31

} } 72. Testing. Since I don’t claim to understand the least common ancestor code completely, let’s test it. ⟨ 72 ⟩ ≡ #include #include #include #include #include #include #include #include #include #include #include #include

tsp instance t ∗tsp instance;

void(∗sort )(void ∗a, size t n, size t es , int(∗cmp)(const void ∗, const void ∗)); int verbose = 100; int noround = 0; ⟨ Test module definitions 76 ⟩ ⟨ Test module variables 74 ⟩ ⟨ Test module subroutines 73 ⟩ int main(int argc, char ∗∗argv ) {

length t len;

decluster tree t ∗mst ;

tsp instance = new of (tsp instance t); sort = dsort ; mem usage reset ( ); resource setup(100); ⟨ Test on a simple cost matrix 79 ⟩ ⟨ Test case clean up 87 ⟩ ⟨ Test on a Euclidean instance 83 ⟩ ⟨ Test case clean up 87 ⟩ free mem(tsp instance); return 0; /∗ All ok. ∗/

32  TESTING  DECLUSTER §73

} 73. Here’s a simple cost matrix. It is equivalent to a Euclidean instance with six points on a line. ⟨ Test module subroutines 73 ⟩ ≡

static length t simple mat cost (const int i, const int j);

static length t simple mat cost (const int i, const int j)

{ extern int verbose; static length t cost matrix [6][6] = {

{0, 10, 60, 71, 91, 121}, /∗ distance from 0 ∗/ {10, 0, 50, 61, 81, 111}, /∗ distance from 1 ∗/ {60, 50, 0, 11, 31, 61}, /∗ distance from 2 ∗/ {71, 61, 11, 0, 20, 50}, /∗ distance from 3 ∗/ {91, 81, 31, 20, 0, 30}, /∗ distance from 4 ∗/ {121, 111, 61, 50, 30, 0} /∗ distance from 5 ∗/ };

errorif (i < 0 ∨ i > 5,  ␣      ␣   ␣  ␣     ␣  , i);
errorif (j < 0 ∨ j > 5,  ␣      ␣   ␣  ␣     ␣  , j);
if (verbose ≥ 1000) printf ( ␣␣␣  length t spec  , i, j, cost matrix [i][j]);

return cost matrix [i][j]; }

See also sections 77 and 78. This code is used in section 72.

74. Here is the minimum spanning tree we expect. It is unique up to edge reordering. ⟨ Test module variables 74 ⟩ ≡

static decluster edge t simple mat MST edges [ ] = {{{0, 1}, 10}, {{2, 3}, 11}, {{3, 4}, 20}, {{4, 5}, 30},

{{1, 2}, 50}}; static int simple mat n = 1 + sizeof (simple mat MST edges )/sizeof (decluster edge t);

See also sections 80, 84, and 85. This code is used in section 72.

75. ⟨ Module definitions 41 ⟩ +≡ #if #define print tree(A, B)decluster print tree (stdout , A, B) #else #define print tree (A, B) #endif

76. ⟨ Test module definitions 76 ⟩ ≡ #if #define print tree decluster print tree #else #define print tree (A, B) #endif

This code is used in section 72.

77. Since the minimum spanning trees I use are unique only up to reordering, it is convenient to have a canonicalizing function. ⟨ Test module subroutines 73 ⟩ +≡ #define swap(type, A, B)

{ typetmp; tmp = (A); (A) = (B); (B) = tmp;

} /∗ Messes up if A or B is the name tmp, or if either causes side effects. ∗/ static void canonicalize(decluster tree t ∗t) {

int n = t n, i; extern int verbose;

errorif (t ≡ Λ,  ␣      ␣ ␣    ␣  );
errorif (n < 0,  ␣    ␣     ␣  ␣ ␣  , n);

if (verbose ≥ 1000) print tree(t, ); dsort (t edge, (size t) n, sizeof (decluster edge t), decluster edge cmp); for (i = 0; i < n; i^{++}) {

if (t edge[i].city [0] > t edge[i].city [1]) swap(int, t edge[i].city [0], t edge[i].city [1]); } if (verbose ≥ 1000) print tree(t, );

} 78. ⟨ Test module subroutines 73 ⟩ +≡ #define max (A, B) ((A) > (B) ? (A) : (B))

static void compare trees (decluster tree t ∗a, decluster tree t ∗b) {

int i, n, ok = 1; if (verbose ≥ 200) printf ( , a n); if (verbose ≥ 200) printf ( , b n); if (a n ≠ b n) ok = 0, printf ( ␣ ␣ ␣ ␣ ␣ ␣ ␣ , a n, b n); n = a n; if (ok )

for (i = 0; i < n; i^{++}) { decluster edge t ea = a edge[i], eb = b edge[i]; if (ea.city [0] ≠ eb.city [0]) ok = 0,

printf ( ␣ ␣ ␣ ,

i, ea.city [0], i, eb.city [0]);

if (ea.city [1] ≠ eb.city [1]) ok = 0, printf ( ␣ ␣ ␣ ,

i, ea.city [1], i, eb.city [1]);

if (ea.cost ≠ eb.cost )

ok = 0, printf ( ␣ length t spec ␣ ␣

length t spec , i, length t pcast (ea.cost ), i, length t pcast (eb.cost ));

} if (¬ok ) {

int n = max (a n, b n); printf ( , a n, b n); for (i = 0; i < n; i^{++}) {

decluster edge t ea = a edge[i], eb = b edge[i]; printf ( ␣ , i); if (i < a n) printf ( ␣ ␣ length t spec, ea.city [0], ea.city [1], ea.cost ); else printf ( ); if (i < b n) printf ( ␣ ␣ length t spec, eb.city [0], eb.city [1], eb.cost ); printf ( );

} errorif (1, ␣ ␣ );

} } 79. Variable cost is owned by the READ module. ⟨ Test on a simple cost matrix 79 ⟩ ≡

{ decluster tree t simple mat MST ; simple mat MST .n = simple mat n − 1; simple mat MST .edge = simple mat MST edges ; printf ( ␣ ␣ ); tsp instance name = dup string ( ); tsp instance comment = dup string ( ␣ ␣ ␣ ); tsp instance n = simple mat MST .n + 1; tsp instance edge weight type = ; tsp instance edge weight format = ; tsp instance edge weights = Λ; /∗ We supply our own cost function. ∗/ tsp instance short edge weights = Λ; tsp instance coord = Λ; tsp instance xmin = 0; tsp instance xmax = 0; tsp instance ymin = 0; tsp instance ymax = 0; mst = decluster setup(simple mat n); printf ( ); errorif (mst ≡ Λ, ␣ ␣ ␣ ␣ ␣ ); cost = simple mat cost ; len = decluster mst (tsp instance, mst ); printf ( ); canonicalize(&simple mat MST ); printf ( ); canonicalize(mst ); printf ( ); compare trees (&simple mat MST , mst ); printf ( ); decluster preprocess (mst ); printf ( ); decluster cleanup tree(mst ); printf ( ); printf ( ␣ );

}

See also section 81. This code is used in section 72.

80. Now let’s test the LCA’s exhaustively. ⟨ Test module variables 74 ⟩ +≡

static length t simple mat lca len[6][6] = {  /∗ Row 0 ∗/ {0, 10, 50, 50, 50, 50},  /∗ Row 1 ∗/
{10, 0, 50, 50, 50, 50},  /∗ Row 2 ∗/ {50, 50, 0, 11, 20, 30},  /∗ Row 3 ∗/ {50, 50, 11, 0, 20, 30},

/∗ Row 4 ∗/ {50, 50, 20, 20, 0, 30}, /∗ Row 5 ∗/ {50, 50, 30, 30, 30, 0}}; 81. ⟨ Test on a simple cost matrix 79 ⟩ +≡

{ int n = simple mat n; print tree(decluster topology tree( ), ␣ ␣ );

#define lca table simple mat lca len

⟨ Test lca table 82 ⟩ #undef lca table

printf ( ␣ ); }

82. ⟨ Test lca table 82 ⟩ ≡

{ int i, j, k; for (k = 0; k < 2; k^{++}) { /∗ Twice, for good measure. ∗/

printf ( );

for (i = 0; i < n; i^{++}) { printf ( ); for (j = 0; j < n; j^{++}) {

const length t d = decluster d (i, j);

printf ( );

if (verbose ≥ 500) printf ( ␣ ␣ length t spec, i, j, length t pcast (d)); if (d ≠ lca table[i][j]) {

printf ( ␣ ␣ ␣ ␣ length t spec ␣ ␣ length t spec, i, j, length t pcast (lca table[i][j]), length t pcast (d)); errorif (1, ␣ );

}

} } printf ( );

} }

This code is used in sections 81 and 86.

83. ⟨ Test on a Euclidean instance 83 ⟩ ≡

{ decluster tree t euc MST ; euc MST .n = euc n − 1; euc MST .edge = euc MST edges ; printf ( ␣ ␣ ); tsp instance name = dup string ( ); tsp instance comment = dup string ( ␣ ␣ ␣ ); tsp instance n = euc MST .n + 1; tsp instance edge weight type = ; tsp instance edge weight format = ; tsp instance edge weights = Λ; tsp instance short edge weights = Λ; tsp instance coord = euc coord ; tsp instance xmin = 0; tsp instance xmax = 0; tsp instance ymin = 0; tsp instance ymax = 121; E2 create(tsp instance); printf ( ); mst = decluster setup(euc n); printf ( ); errorif (mst ≡ Λ, ␣ ␣ ␣ ␣ ␣ ); switch to(tsp instance); len = decluster mst (tsp instance, mst ); printf ( ); canonicalize(&euc MST ); printf ( ); canonicalize(mst ); printf ( ); compare trees (&euc MST , mst ); printf ( ); decluster preprocess (mst ); printf ( ); decluster cleanup tree(mst ); printf ( ); E2 destroy ( ); printf ( ); printf ( ␣ );

}

See also section 86. This code is used in section 72.

84. Here is the minimum spanning tree we expect. It is unique up to edge reordering. ⟨ Test module variables 74 ⟩ +≡

static coord 2d euc coord [ ] = {{{0, 0}}, {{0, 10}}, {{0, 85}}, {{0, 96}}, {{0, 116}}, {{0, 146}}}; static decluster edge t euc MST edges [ ] = {{{0, 1}, 10}, {{2, 3}, 11}, {{3, 4}, 20}, {{4, 5}, 30}, {{1, 2},

75}}; static int euc n = 1 + sizeof (euc MST edges )/sizeof (decluster edge t); 85. Now test LCA for the Euclidean instance. ⟨ Test module variables 74 ⟩ +≡

static length t euc lca len[6][6] = { /∗ Row 0 ∗/ {0, 10, 75, 75, 75, 75}, /∗ Row 1 ∗/ {10, 0, 75, 75, 75, 75}, /∗ Row 2 ∗/ {75, 75, 0, 11, 20, 30}, /∗ Row 3 ∗/ {75, 75, 11, 0, 20, 30}, /∗ Row 4 ∗/ {75, 75, 20, 20, 0, 30}, /∗ Row 5 ∗/ {75, 75, 30, 30, 30, 0}};

86. ⟨ Test on a Euclidean instance 83 ⟩ +≡

{ int n = euc n; print tree(decluster topology tree( ), ␣ );

#define lca table euc lca len

⟨ Test lca table 82 ⟩ #undef lca table

printf ( ␣ ); }

87. ⟨ Test case clean up 87 ⟩ ≡

decluster cleanup( );

printf ( );

This code is used in section 72.

§88  DECLUSTER  VERBOSE OUTPUT 39

88. Verbose output. ⟨ Verbose: print the preprocessing data 88 ⟩ ≡ #if

{ int i; printf ( ␣␣␣ ␣ ␣ ␣ , , , , ); for (i = 0; i < n + n − 1; i^{++}) {

printf ( ␣ ␣ ␣ ␣ length t spec , i, digest [i].level , digest [i].inlabel ,

digest [i].ascendant , length t pcast (digest [i].cost )); } printf ( ␣ ␣ ␣ ␣ ); for (i = 1; i < n + n; i^{++}) {

printf ( ␣ ␣␣␣ , i, parent of head [i]); }

} #endif

This code is used in section 35.

89. ⟨ Verbose: print parent 89 ⟩ ≡ #if

{ int i; printf ( ␣ ); for (i = 0; i < n + n − 1; i^{++}) printf ( ␣ ␣␣␣ , i, parent [i]);

} #endif

This code is used in section 60.

90. ⟨ Verbose: print head 90 ⟩ ≡ #if

{ int i; printf ( ␣ ); for (i = 1; i < n + n; i^{++}) printf ( ␣ ␣␣␣ , i, head [i]);

} #endif

This code is used in section 60.

91. ⟨ Module variables 17 ⟩ +≡ #if

static char pp[ ] = ␣␣␣␣␣␣␣␣␣␣ ; /∗ Print prefix. ∗/ #endif

40  VERBOSE OUTPUT  DECLUSTER §92

92. ⟨ Verbose: print b 92 ⟩ ≡ #if

if (verbose ≥ 2000) { printf ( ␣ ␣ ␣ ␣ , pp, x, xd .level , xd .inlabel , xd .ascendant ); printf ( ␣ ␣ ␣ ␣ , pp, y, yd .level , yd .inlabel , yd .ascendant ); printf ( ␣ , pp, b);

} #endif

This code is used in section 65.

93. ⟨ Verbose: print inlabel z , common, and jmask 93 ⟩ ≡ #if

if (verbose ≥ 2000) {

printf ( ␣ ␣ ␣ , pp, inlabel z , common, jmask ); }

#endif

This code is used in section 66.

94. ⟨ Verbose: print inlabel w x , xpp 94 ⟩ ≡ #if

if (verbose ≥ 2000) {

printf ( ␣ ␣ , pp, inlabel w x , xpp); }

#endif

This code is used in section 67.

95. ⟨ Verbose: print inlabel w y , ypp 95 ⟩ ≡ #if

if (verbose ≥ 2000) {

printf ( ␣ ␣ , pp, inlabel w y , ypp); }

#endif

This code is used in section 67.

§96  DECLUSTER  INDEX 41

96. Index.

a: \underline{11}, \underline{28}, \underline{72}, \underline{78}. ancestor: 46. anchor city: 5. ans : \underline{44}. ap: \underline{59}. argc: \underline{72}. argv : \underline{72}. ascendant : 35, \underline{48}, 58, 59, 66, 67, 88, 92. ascendant numbers: 48. b: \underline{11}, \underline{28}, \underline{64}, \underline{78}. barrel shifters: 44. best gain: 4, 5. binary search: 44. bottleneck cost: 5. bridge: \underline{26}. bridge pq : \underline{26}. canonicalize: \underline{77}, 79, 83. ch: \underline{53}. child : \underline{30}, 40, 53, 55, 57, \underline{60}. child0 : \underline{57}. child1 : \underline{57}. city : \underline{19}, 24, 26, 30, 40, 71, 77, 78. cluster distance: 1, 5, 8. cmp: 72. comment : 79, 83. common: \underline{64}, 66, 93. compare trees : \underline{78}, 79, 83. component : \underline{38}, 40. coord : 79, 83. copy 1 down: \underline{41}, 65, 67. cost : 11, 15, \underline{19}, 22, 23, 24, 25, 26, 28, 30, \underline{48},

50, 52, 68, 71, 78, 79, 88. cost matrix : \underline{73}. cumulative gain criterion: 2, 5. cur child : \underline{55}. cur level : \underline{53}. d: \underline{25}, \underline{82}. decluster cleanup: 10, \underline{11}, \underline{18}, 87. decluster cleanup tree: 10, \underline{11}, \underline{21}, 79, 83. decluster d : 1, 10, \underline{11}, \underline{68}, 82. : 75, 76, 88, 89, 90, 91, 92,

93, 94, 95. decluster discard topology tree: 11, 35, \underline{36}, \underline{37}. decluster edge cmp: 10, \underline{11}, 26, \underline{28}, 38, 77. decluster lca: \underline{64}, 68. decluster mst : 10, \underline{11}, \underline{22}, 79, 83. decluster mst custom: 10, \underline{11}, 22, \underline{23}. decluster preprocess : 10, \underline{11}, \underline{35}, 79, 83. decluster print tree: \underline{11}, \underline{71}, 75, 76. decluster rcs id : \underline{9}, \underline{13}. decluster setup: 10, \underline{11}, \underline{16}, 79, 83. decluster topology tree: \underline{11}, \underline{70}, 81, 86. descendant: 46. : \underline{58}, 59. : \underline{54}, 56, 57. digest : \underline{49}, 50, 51, 52, 53, 57, 59, 60, 64, 68, 88. digest t: \underline{48}, 49, 50, 51, 64. dist : \underline{11}, \underline{22}, \underline{23}, 25. dsort : 72, 77. dup string : 79, 83. ea: \underline{78}. eb: \underline{78}. edge: \underline{19}, 20, 21, 24, 26, 30, 32, 38, 40, 52, 53,

55, 57, 60, 71, 77, 78, 79, 83. edge weight format : 79, 83. edge weight type: 79, 83. edge weights : 79, 83. errorif : 22, 23, 26, 35, 71, 73, 77, 78, 79, 82, 83. es : \underline{72}. euc coord : 83, \underline{84}. euc lca len: \underline{85}, 86. euc MST : \underline{83}. euc MST edges : 83, \underline{84}. euc n: 83, \underline{84}, 86. : 83. : 79. E2 create: 83. E2 destroy : 83. E2 hide: 26. E2 nn: 26. E2 supports : 22. E2 unhide all : 26. floor log 2 : \underline{44}, 57. floor log 2 small : \underline{43}, 44. fprintf : 71. free mem: 21, 22, 26, 32, 38, 51, 53, 54, 55,

60, 63, 72. from: \underline{11}, \underline{22}, \underline{23}, 24, 25. : 79. grey width: 44. head : 48, \underline{60}, 90. head numbers: 48. here: \underline{40}, \underline{53}, \underline{55}, 56, 57, 59. hi diff : \underline{65}. hi mask : \underline{45}, 57. hii : \underline{60}. i: \underline{23}, \underline{30}, \underline{38}, \underline{40}, \underline{50}, \underline{52}, \underline{53}, \underline{57}, \underline{60}, \underline{71}, \underline{73}, \underline{77},

\underline{78}, \underline{82}, \underline{88}, \underline{89}, \underline{90}.

ih: \underline{59}. ii : \underline{60}. imask : \underline{66}. in: \underline{26}, \underline{55}, 59.

42  INDEX  DECLUSTER §96

in-order traversal: 65. index rightmost 1 : 66. : 23, 25. inlabel : 35, \underline{48}, 56, 57, 58, 59, 60, 64, 65, 66,

67, 88, 92. inlabel numbers: 48, 65. inlabel w : 60. inlabel w x : \underline{67}, 94. inlabel w y : \underline{67}, 95. inlabel z : \underline{64}, 66, 67, 93. ip: \underline{59}. is in component : \underline{26}. j: \underline{60}, \underline{73}, \underline{82}. jmask : \underline{64}, 66, 67, 93. k: \underline{82}. kmask : \underline{67}. Kruskal’s algorithm: 8. last : \underline{57}. LCA: 8. lca table: \underline{81}, 82, \underline{86}. least common ancestor: 8, 46. len: \underline{72}, 79, 83. len diff : \underline{28}. length t pcast : 24, 71, 78, 82, 88. length t spec: 24, 71, 73, 78, 82, 88. level : 35, \underline{48}, 53, 60, 64, 88, 92. level numbers: 48. lo mask : \underline{45}, 66. lomask : \underline{65}. main: \underline{72}. max : \underline{65}, \underline{78}. mem deduct : 21, 22, 26, 32, 38, 51, 53, 54,

55, 60, 63. mem usage reset : 72. memory leak: 30. metric: 7. min: \underline{65}. minimum spanning tree: 8. mst : \underline{16}, 20, \underline{72}, 79, 83. MST: 8. N: \underline{57}. n: \underline{11}, \underline{17}, \underline{19}, \underline{23}, \underline{32}, \underline{71}, \underline{72}, \underline{77}, \underline{78}, \underline{81}, \underline{86}. name: \underline{11}, \underline{71}, 79, 83. new arr of : 20, 22, 26, 30, 38, 50, 53, 54, 55,

60, 62. new arr of zero: 26. new inside city : \underline{25}. new of : 20, 30, 72. next child : \underline{55}. next edge: \underline{23}, 24, \underline{26}. : \underline{30}, 53, 55, 57, 60. : 83.

: \underline{38}, 40, 60. noround : \underline{72}. ok : \underline{78}. out : \underline{11}, \underline{26}, \underline{71}. p: \underline{59}. parent : \underline{40}, \underline{60}, 89. parent of head : 60, \underline{61}, 62, 63, 67, 88. positive-definite: 7. pow 2 : \underline{42}, 45. pp: \underline{91}, 92, 93, 94, 95. pq create size: 26. pq delete min: 26. pq destroy : 26. pq insert : 26. preorder : \underline{54}, 55, 56, 57. preorder number : \underline{54}, 56. print name: \underline{71}. print tree: 35, 38, \underline{75}, \underline{76}, 77, 81, 86. printf : 24, 73, 78, 79, 81, 82, 83, 86, 87, 88,

89, 90, 92, 93, 94, 95. qsort : 28. queue: \underline{53}. r: \underline{21}, \underline{38}, \underline{53}. resource setup: 72. short bridge: \underline{26}. short edge weights : 79, 83. short len: \underline{23}, 24, 25. short to: \underline{23}, 24, 25. simple mat cost : \underline{73}, 79. simple mat lca len: \underline{80}, 81. simple mat MST : \underline{79}. simple mat MST edges : \underline{74}, 79. simple mat n: \underline{74}, 79, 81. size: \underline{54}, 55, 57. : 41, 42, 44. : 44. sort : 38, 72. stdout : 75. swap: \underline{77}. switch to: 83. T: \underline{11}, \underline{21}, \underline{22}, \underline{23}, \underline{35}. t: \underline{11}, \underline{65}, \underline{71}, \underline{77}. T prime: 30, 32, \underline{34}, 35, 36, 38, 40, 48, 52, 53,

55, 57, 60, 65, 68, 70. the n: \underline{16}. tmp: 77. top: \underline{55}, 59. total len: \underline{22}, \underline{23}, 24, 26. tree B: 65. triangle inequality: 7. tsp instance: \underline{11}, \underline{22}, \underline{72}, 79, 83. type: 77. u: \underline{11}, \underline{66}, \underline{68}. v: \underline{11}, \underline{68}. verbose: \underline{22}, 24, \underline{72}, \underline{73}, \underline{77}, 78, 82, 92, 93, 94, 95. w: \underline{38}, \underline{53}. word width: 45. x: \underline{44}, \underline{64}. xd : \underline{64}, 66, 67, 92. xfuzz : \underline{65}. xil : \underline{64}, 65, 67. xmax : 79, 83. xmin: 79, 83. xpp: \underline{64}, 67, 94. y: \underline{64}. yd : \underline{64}, 66, 67, 92. yfuzz : \underline{65}. yil : \underline{64}, 65, 67. ymax : 79, 83. ymin: 79, 83. ypp: \underline{64}, 67, 95.

⟨ Add the edge to short to 24 ⟩ Used in section 23.

⟨ Allocate and pre-initialize T prime 30 ⟩ Used in sections 31 and 35. ⟨ Allocate space for a spanning tree mst 20 ⟩ Used in section 16. ⟨ Build an MST using k-d trees 26 ⟩ Used in section 22.

⟨ Clean up T prime 32 ⟩ Used in sections 30, 33, and 35. ⟨ Compute ascendant numbers 58 ⟩ Used in section 35. ⟨ Compute inlabel numbers 54 ⟩ Used in section 35. ⟨ Compute level numbers 53 ⟩ Used in section 35. ⟨ Compute parent of head numbers 60 ⟩ Used in section 35. ⟨ Copy cost fields from T prime 52 ⟩ Used in section 35. ⟨ Do a depth-first search of T prime 55 ⟩ Used in sections 54 and 58.

⟨ Early module headers 15 ⟩ Used in section 9. ⟨ Edge comparison function 28 ⟩ Used in section 29. ⟨ Exported subroutines 11 ⟩ Used in section 13.

⟨ Exported types 19 ⟩ Used in section 13.

⟨ Exported variables 37 ⟩ Used in section 13. ⟨ Fill the last n − 1 entries of T prime with data from T 38 ⟩ Used in section 35.

⟨ Find inlabel z , common, and jmask 66 ⟩ Used in section 64. ⟨ Find xpp and ypp 67 ⟩ Used in section 64.

⟨ Global variables 36 ⟩ Used in section 9. ⟨ Link T prime edge[w] with its two children and update component 40 ⟩ Used in section 38.

⟨ Module definitions 41, 75 ⟩ Used in section 9. ⟨ Module headers 14, 27, 39 ⟩ Used in section 9. ⟨ Module subroutines 44, 64 ⟩ Used in section 9. ⟨ Module types 48 ⟩ Used in section 9. ⟨ Module variables 17, 34, 42, 43, 49, 61, 91 ⟩ Used in section 9. ⟨ Other cleanup code 33, 51, 63 ⟩ Used in section 18. ⟨ Other setup code 31, 50, 62 ⟩ Used in section 16.

⟨ Set b to the LCA of C_{x} and C_{y} in the complete tree B 65 ⟩ Used in section 64.

⟨ Subroutines 16, 18, 21, 22, 23, 29, 35, 68, 70, 71 ⟩ Used in section 9. ⟨ System headers 12 ⟩ Used in section 9. ⟨ Test case clean up 87 ⟩ Used in section 72.

⟨ Test module definitions 76 ⟩ Used in section 72.

⟨ Test module subroutines 73, 77, 78 ⟩ Used in section 72. ⟨ Test module variables 74, 80, 84, 85 ⟩ Used in section 72.

⟨ Test on a Euclidean instance 83, 86 ⟩ Used in section 72. ⟨ Test on a simple cost matrix 79, 81 ⟩ Used in section 72.

⟨ Test lca table 82 ⟩ Used in sections 81 and 86.

⟨ Update distance arrays and short to 25 ⟩ Used in section 23. ⟨ Verbose: print the preprocessing data 88 ⟩ Used in section 35.

⟨ Verbose: print b 92 ⟩ Used in section 65. ⟨ Verbose: print head 90 ⟩ Used in section 60. ⟨ Verbose: print inlabel w x , xpp 94 ⟩ Used in section 67. ⟨ Verbose: print inlabel w y , ypp 95 ⟩ Used in section 67.

⟨ Verbose: print inlabel z , common, and jmask 93 ⟩ Used in section 66.

⟨ Verbose: print parent 89 ⟩ Used in section 60. ⟨ Visit here the first time 56, 59 ⟩ Used in section 55.

⟨ Visit here the last time 57 ⟩ Used in section 55. ⟨ 13 ⟩ ⟨ 72 ⟩

DECLUSTER

Section Page

Declustering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 5 Interface and code preamble . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 8 Setup and cleanup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 10 Minimum spanning trees . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 11 Topology trees . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 17 Bit twiddling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 20 Least common ancestors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 22 Least common ancestor queries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 28 Extras . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 30 Testing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 31 Verbose output . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 39 Index . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 41