Distances on Directed Graphs in R
mpadge.tngl.sh
1/* File heap23.c - 2-3 Heap
2 * ----------------------------------------------------------------------------
3 * Mark Padgham, adapted from code by Shane Saunders
4 */
5
6/* This version is implemented using the same pointer structure as a Fibonacci
7 * heap; that is, nodes have a parent pointer, and a child pointer which points
8 * to a linked list of children constructed from the left and right pointers of
9 * nodes. In this implementation, the child pointer points to the highest
10 * dimension child node.
11 */
12#include <cmath>
13#include <cstdio>
14#if defined(TTHEAP_DUMP) && TTHEAP_DUMP > 0
15#include <Rcpp.h>
16#endif
17#include "heap23.h"
18
19/*--- Heap23 (public methods) -----------------------------------------------*/
20
21/* --- Constructor ---
22 * Allocates a 2-3 heap capable of holding $n$ items.
23 */
24Heap23::Heap23(size_t n)
25{
26#if defined(TTHEAP_DUMP) && TTHEAP_DUMP > 0
27 Rcpp::Rcout << "init, ";
28 Rcpp::Rcout.flush ();
29#endif
30
31 /* The maximum number of nodes and the maximum number of trees allowed.
32 */
33 maxNodes = n;
34 maxTrees = static_cast <size_t>(0.5 +
35 log(static_cast <double> (n) + 1.0) / log (2.0));
36
37 /* Allocate space for an array of pointers to trees, and nodes in the heap.
38 * Initialise all array entries to zero, that is, NULL pointers.
39 */
40 trees = new Heap23Node *[maxTrees];
41 for(size_t i = 0; i < maxTrees; i++) trees[i] = std::nullptr_t ();
42
43 nodes = new Heap23Node *[n];
44 for(size_t i = 0; i < n; i++) nodes[i] = std::nullptr_t ();
45
46 /* We begin with no nodes in the heap. */
47 itemCount = 0;
48
49 /* The value of the heap helps to keep track of the maximum rank as nodes
50 * are inserted and deleted.
51 */
52 treeSum = 0;
53
54 /* For experimental purposes, we keep a count of the number of key
55 * comparisons.
56 */
57 compCount = 0;
58}
59
60/* --- Destructor ---
61*/
62Heap23::~Heap23()
63{
64#if defined(TTHEAP_DUMP) && TTHEAP_DUMP > 0
65 Rcpp::Rcout << "free, ";
66 Rcpp::Rcout.flush ();
67#endif
68
69 for(size_t i = 0; i < maxNodes; i++) delete nodes[i];
70
71 delete [] nodes;
72 delete [] trees;
73#if defined(TTHEAP_DUMP) && TTHEAP_DUMP > 0
74 Rcpp::Rcout << "free-exited, ";
75 Rcpp::Rcout.flush ();
76#endif
77}
78
79/* --- insert() ---
80 * Inserts item $item$, with associated key $k$, into the heap.
81 */
82void Heap23::insert(size_t item, double k)
83{
84 Heap23Node *newNode;
85
86#if defined(TTHEAP_DUMP) && TTHEAP_DUMP > 0
87 Rcpp::Rcout << "insert, ";
88 dump();
89 Rcpp::Rcout.flush ();
90#endif
91
92 /* Create an initialise the new node. The parent pointer will be set to
93 * NULL by meld().
94 */
95 newNode = new Heap23Node;
96 newNode->child = std::nullptr_t ();
97 newNode->left = newNode->right = std::nullptr_t ();
98 newNode->dim = 0;
99 newNode->item = item;
100 newNode->key = k;
101
102 /* Maintain a pointer to item's new node in the heap. */
103 nodes[item] = newNode;
104
105 /* Meld the new node into the heap. */
106 meld(newNode);
107
108 /* Update the heap's node count. */
109 itemCount++;
110
111#if defined(TTHEAP_DUMP) && TTHEAP_DUMP > 0
112 Rcpp::Rcout << "insert-exited, ";
113 dump();
114 Rcpp::Rcout.flush ();
115#endif
116}
117
118/* --- deleteMin() ---
119 * Delete and return the minimum item from the heap.
120 */
121size_t Heap23::deleteMin()
122{
123 Heap23Node *minNode, *child, *next;
124 double k, k2;
125 size_t r, v, item;
126
127#if defined(TTHEAP_DUMP) && TTHEAP_DUMP > 0
128 Rcpp::Rcout << "deleteMin, ";
129 Rcpp::Rcout.flush ();
130#endif
131
132 /* First we determine the maximum rank tree in the heap. */
133 v = treeSum;
134 // MP Note: r was an int formerly intialised at -1, but it's used as a
135 // direct array index below, so really should be unsigned! The r-- at end
136 // fixes that.
137 r = 0;
138 while(v) {
139 v = v >> 1;
140 r++;
141 };
142 r--;
143
144 /* Now locate the root node with the smallest key, scanning from the
145 * maximum rank root position, down to rank 0 root position.
146 */
147 minNode = trees[r];
148 k = minNode->key;
149 while(r > 0) {
150 r--;
151 next = trees[r];
152 if(next) {
153 if((k2 = next->key) < k) {
154 k = k2;
155 minNode = next;
156 }
157 compCount++;
158 }
159 }
160
161 /* We remove the minimum node from the heap but keep a pointer to it. */
162 r = minNode->dim;
163 trees[r] = std::nullptr_t ();
164 treeSum -= (1 << r);
165 itemCount--;
166
167 /* A nodes child pointer always points to the child with the highest rank,
168 * so child->right is the smallest rank. For melding the linked list
169 * starting at child->right we terminate the circular link with a NULL
170 * pointer.
171 */
172 child = minNode->child;
173 if(child) {
174 next = child->right;
175 next->left = child->right = std::nullptr_t ();
176 meld(next);
177 }
178
179 /* Record the vertex no to return. */
180 item = minNode->item;
181
182 /* Delete the old minimum node. */
183 nodes[item] = std::nullptr_t ();
184 delete minNode;
185
186#if defined(TTHEAP_DUMP) && TTHEAP_DUMP > 0
187 Rcpp::Rcout << "deleteMin-exited, ";
188 Rcpp::Rcout.flush ();
189#endif
190
191 return item;
192}
193
194/* --- decreaseKey() ---
195 * Decrease the key of item $item$ in the heap to the value $newValue$. It
196 * is the users reponsibility to ensure that newValue is in-fact less than or
197 * equal to the current value.
198 */
199void Heap23::decreaseKey(size_t item, double newValue)
200{
201 Heap23Node *cutNode, *parent;
202
203#if defined(TTHEAP_DUMP) && TTHEAP_DUMP > 0
204 Rcpp::Rcout << "decreaseKey on vn = " << item;
205 Rcpp::Rcout.flush ();
206#endif
207
208 /* Obtain a pointer to the decreased node and its parent and child.*/
209 cutNode = nodes[item];
210 parent = cutNode->parent;
211 cutNode->key = newValue;
212
213 /* No reinsertion occurs if the node changed was a root. */
214 if(!parent) {
215#if defined(TTHEAP_DUMP) && TTHEAP_DUMP > 0
216 Rcpp::Rcout << "decreaseKey-exited, ";
217 Rcpp::Rcout.flush ();
218#endif
219 return;
220 }
221
222 /* Now remove the node and its tree and reinsert it. */
223 removeNode(cutNode);
224 cutNode->right = cutNode->left = std::nullptr_t ();
225 meld(cutNode);
226
227#if defined(TTHEAP_DUMP) && TTHEAP_DUMP > 0
228 Rcpp::Rcout << "decreaseKey-exited, ";
229 Rcpp::Rcout.flush ();
230#endif
231
232}
233
234/*--- Heap23 (private methods) ----------------------------------------------*/
235
236/* --- meld() ---
237 * Melds the linked list of trees pointed to by $treeList$ into the heap
238 * pointed to by $h$. This function uses the $right$ sibling pointer
239 * of nodes to traverse the linked list from lower dimension nodes to higher
240 * dimension nodes. It expects the last nodes $right$ pointer to be NULL.
241 */
242void Heap23::meld(Heap23Node *treeList)
243{
244 Heap23Node *next, *addTree;
245 Heap23Node *carryTree;
246 size_t d;
247
248#if defined(TTHEAP_DUMP) && TTHEAP_DUMP > 0
249 Rcpp::Rcout << "meld - ";
250 Rcpp::Rcout.flush ();
251#endif
252
253 /* addTree points to the current tree to be merged. */
254 addTree = treeList;
255
256 carryTree = std::nullptr_t ();
257
258 do {
259 /* addTree() gets merged into the heap, and also carryTree if one
260 * exists from a previous merge.
261 */
262
263 /* Keep a pointer to the next tree and remove sibling and parent links
264 * from the current tree. The dimension of the next tree is always
265 * one greater than the dimension of the previous tree, so this merging
266 * is like an addition of two ternary numbers.
267 *
268 * Note that if addTree is NULL and the loop has not exited, then
269 * there is only a carryTree to be merged, so treat it like addTree.
270 */
271 next = std::nullptr_t ();
272 if(addTree) {
273 next = addTree->right;
274 addTree->right = addTree->left = addTree;
275 addTree->parent = std::nullptr_t ();
276 }
277 else {
278 addTree = carryTree;
279 carryTree = std::nullptr_t ();
280 }
281
282#if defined(TTHEAP_DUMP) && TTHEAP_DUMP > 0
283 Rcpp::Rcout << addTree->item << " ";
284 Rcpp::Rcout.flush ();
285#endif
286
287 /* First we merge addTree with carryTree, if there is one. Note that
288 * carryTree contains only one node in its main trunk, and addTree
289 * has at most two, so the result is at most one 3-node trunk, which is
290 * treated as a 1-node main trunk one dimension higher up.
291 */
292 if(carryTree) {
293 compCount += merge(&addTree, &carryTree);
294 }
295
296 /* After the merge, if addTree is NULL, then the resulting tree
297 * pointed to by carryTree carries to higher entry, so we do not need
298 * to merge anything into the existing main trunk.
299 * If addTree is not NULL we add it to the existing main trunk.
300 */
301 if(addTree) {
302 d = addTree->dim;
303 if(trees[d]) {
304 /* Nodes already in this main trunk position, so merge. */
305 compCount += merge(&trees[d], &addTree);
306 if(!trees[d]) treeSum -= (1 << d);
307 carryTree = addTree;
308 }
309 else {
310 /* No nodes in this main trunk position, so use addTree. */
311 trees[d] = addTree;
312 treeSum += (1 << d);
313 }
314 }
315
316 /* Obtain a pointer to the next tree to add. */
317 addTree = next;
318
319 /* We continue if there is still a node in the list to be merged, or
320 * a carry tree remains to be merged.
321 */
322 } while(addTree || carryTree);
323
324
325#if defined(TTHEAP_DUMP) && TTHEAP_DUMP > 0
326 Rcpp::Rcout << "meld-exited, ";
327 Rcpp::Rcout.flush ();
328#endif
329
330}
331
332/* --- removeNode() ---
333 * Removes the node pointed to by $rNode$, and its corresponding sub-tree, from
334 * the heap. If necessary, this causes rearrangement of $rNode$'s work space.
335 */
336void Heap23::removeNode(Heap23Node *rNode)
337{
338 Heap23Node *parent, *child, *ax, *bx, *ap, *bp, *b1, *c, *p;
339 size_t d, d1;
340
341 parent = rNode->parent;
342 child = rNode->child;
343 d = rNode->dim;
344
345 /* If this node is an extra node we simply cut the link between it and its
346 * parent and update its sibling pointers.
347 */
348 if(d == parent->dim) {
349 trimExtraNode(rNode);
350 }
351 /* Else if its child is an extra node then use its child to replace it. */
352 else if(child && child->dim == d) {
353
354 /* First we remove the child. */
355 trimExtraNode(child);
356
357 /* Now we put the child in rNodes position. */
358 replaceNode(rNode, child);
359
360 }
361 /* Otherwise we need some rearrangement of the workspace. */
362 else {
363
364 /* Look at up to two similar nodes in the work space and determine if
365 * they have an extra node under them. Nodes relative to the node
366 * being removed are pointed to by the pointers ax, ap, bx, and bp.
367 * If a similar trunk lies immediately below cutNode's trunk in the
368 * work space, then either ax or ap will be set to point to the node on
369 * the end of that trunk. The same applies for bx and bp, but with a
370 * similar trunk immediately above in the work space. the 'x' pointers
371 * are set if there is a 3rd (i.e. extra) node on the trunk. Otherwise
372 * the 'p' pointer is set to point to the 2nd node. Pointers will be
373 * set to null if a trunk does not exist or they are not used.
374 */
375
376 /* Check for nodes on a similar trunk above in the work space. */
377 p = rNode->parent->left;
378 if (p->dim == d) {
379 c = p->child;
380 if(c && c->dim == d) {
381 ax = c; ap = std::nullptr_t ();
382 }
383 else {
384 ap = p; ax = std::nullptr_t ();
385 }
386 }
387 else {
388 ax = ap = std::nullptr_t ();
389 }
390
391 /* Check for nodes on a similar trunk below in the work space. */
392 d1 = d + 1;
393 p = rNode->right;
394 if (p->dim == d1) {
395 p = p->child;
396 if(p->dim == d1) p = p->left;
397
398 c = p->child;
399 if(c && c->dim == d) {
400 bx = c; bp = std::nullptr_t ();
401 }
402 else {
403 bp = p; bx = std::nullptr_t ();
404 }
405 }
406 else {
407 bx = bp = std::nullptr_t ();
408 }
409
410
411 if(bx) {
412
413 /* First break `bx's parent link and sibling links. */
414 trimExtraNode(bx);
415
416 /* Then we insert bx in rNodes place. */
417 replaceNode(rNode, bx);
418 }
419 else if(bp) {
420
421 b1 = bp->parent;
422
423 /* Recursively remove b1. */
424 removeNode(b1);
425 b1->dim = d;
426
427 replaceNode(rNode, b1);
428
429 /* It may improve speed by using trimExtraNode() when recursion can be
430 * avoided.
431 */
432 }
433 else if(ax) {
434
435 /* Bend the tree to modify its shape then remove rNode. */
436 swapTrunks(ax->parent, parent);
437 trimExtraNode(rNode);
438 }
439 else if(ap) {
440
441 /* Bend the tree, so that the node to be relocated, parent, has the
442 * larger key value.
443 */
444 if(parent->key < ap->key) {
445 swapTrunks(ap, parent);
446 p = parent;
447 parent = ap;
448 ap = p;
449 }
450 compCount++;
451
452 trimExtraNode(rNode);
453 removeNode(parent);
454 parent->dim = d;
455
456 /* Make parent the child of ap. */
457 addChild(ap, parent);
458 }
459 else {
460 /* The work space only has rNode node and parent. This only
461 * occurs when parent is a root node, so after removing rNode we
462 * demote parent to a lower dimension main trunk.
463 */
464
465 /* Note that parent is a root node and has dimension d + 1. */
466 trees[d+1] = std::nullptr_t ();
467 treeSum -= (1 << (d+1));
468
469 parent->dim = d;
470 trimExtraNode(rNode);
471 parent->left = parent->right = std::nullptr_t ();
472
473 meld(parent);
474 }
475 }
476}
477
478/*--- Heap23 (node manipulation methods) ------------------------------------*/
479
480/* --- merge() ---
481 * Merges the two trunks pointed to by $*a$ and $*b$, returning the sum trunk
482 * through $a$ and any carry tree through $b$. When this function is used,
483 * both parameters $a$ and $b$ refer to either a 1-node or 2-node trunk.
484 *
485 * Returns the number of key comparisons used.
486 */
487size_t Heap23::merge(Heap23Node **a, Heap23Node **b)
488{
489 Heap23Node *tree, *nextTree, *other, *nextOther;
490 size_t c;
491
492 /* Number of comparisons. */
493 c = 0;
494
495 /* `tree' always points to the node with the lowest key.
496 * To begin with, `tree' points to the smaller head node, and `other'
497 * points to the head node of the other trunk.
498 */
499 if((*a)->key <= (*b)->key) {
500 tree = (*a);
501 other = (*b);
502 }
503 else {
504 tree = (*b);
505 other = (*a);
506 }
507 c++;
508
509 /* nextTree points to the next node on the trunk that `tree' is the head
510 * of (if there is another node).
511 * nextOther points to the next node on the trunk that `other' is the head
512 * of (if there is another node).
513 */
514 nextTree = tree->child;
515 if(nextTree && nextTree->dim != other->dim)
516 nextTree = std::nullptr_t ();
517 nextOther = other->child;
518 if(nextOther && nextOther->dim != other->dim)
519 nextOther = std::nullptr_t ();
520
521 /* The merging depends on the existence of nodes and the values of keys. */
522 if(!nextTree) {
523 /* nextTree does not exist, so we simply make `other' the child of
524 * `tree'. If nextOther exist the resulting 3-node trunk is a carry
525 * tree.
526 */
527
528 addChild(tree, other);
529 if(nextOther) {
530 tree->dim++;
531 *a = std::nullptr_t ();
532 *b = tree;
533 }
534 else {
535 *a = tree;
536 *b = std::nullptr_t ();
537 }
538 }
539 else if(!nextOther) {
540 /* nextTree exists but nextOther does not, so the linked order of
541 * nextTree and `other' in the resulting 3-node trunk depends on the
542 * values of keys. The resulting 3-node trunk becomes a carry tree.
543 */
544
545 if(nextTree->key <= other->key) {
546 addChild(nextTree, other);
547 }
548 else {
549 replaceNode(nextTree, other);
550 addChild(other, nextTree);
551 }
552 c++;
553 tree->dim++;
554 *a = std::nullptr_t ();
555 *b = tree;
556 }
557 else {
558 /* Otherwise, both nextTree and nextOther exist. The result consists
559 * of a 1 node trunk plus the 3-node trunk which becomes a carry tree.
560 * We two trunks are made up as (tree, other, nextOther)
561 * and (nextTree). This uses no key comparisons.
562 */
563
564 replaceNode(nextTree, other);
565 nextTree->left = nextTree->right = nextTree;
566 nextTree->parent = std::nullptr_t ();
567 tree->dim++;
568 *a = nextTree; *b = tree;
569 }
570
571 return c;
572}
573
574/* --- trimExtraNode() ---
575 * Trims the extra node pointed to by $x$ from the trunk to
576 * which it belongs.
577 */
578void Heap23::trimExtraNode(Heap23Node *x)
579{
580 Heap23Node *l, *r;
581
582 if(x->dim == 0) {
583 /* A dimension 0 node is an only child, so cutting it leaves no
584 * children.
585 */
586
587 x->parent->child = std::nullptr_t ();
588 }
589 else {
590 /* Otherwise, sibling pointers of other child nodes must be updated. */
591
592 l = x->left;
593 r = x->right;
594 l->right = r;
595 r->left = l;
596
597 x->parent->child = l;
598 }
599}
600
601/* --- swapTrunks() ---
602 * Where a node in an (i)th trunk, $lowNode$, and a node in an (i+1)th trunk,
603 * $highNode$, share the same parent, this function is used for swapping them.
604 */
605void Heap23::swapTrunks(Heap23Node *lowNode, Heap23Node *highNode)
606{
607 size_t d;
608 Heap23Node *parent, *l, *r;
609
610 /* The dimensions of the two nodes are exchanged. */
611 d = lowNode->dim;
612 lowNode->dim = highNode->dim;
613 highNode->dim = d;
614
615 /* Obtain a pointer to the parent of both nodes. */
616 parent = highNode->parent;
617
618 /* If the left sibling of lowNode is not highNode, we need to update sibling
619 * pointers. Otherwise, the child pointer of the common parent now
620 * points to lowNode.
621 */
622 if((l = lowNode->left) != highNode) {
623
624 /* Update sibling pointers. */
625 r = highNode->right;
626 highNode->left = l;
627 lowNode->right = r;
628 highNode->right = lowNode;
629 lowNode->left = highNode;
630 l->right = highNode;
631 r->left = lowNode;
632
633 /* Determine if the child pointer of the common parent will need to be
634 * updated.
635 */
636 if(parent->child == highNode) {
637 parent->child = lowNode;
638 }
639 }
640 else {
641 parent->child = lowNode;
642 }
643}
644
645/* --- addChild() ---
646 * Makes node $c$ and its tree a child of node $p$.
647 */
648void Heap23::addChild(Heap23Node *p, Heap23Node *c)
649{
650 Heap23Node *l, *r;
651
652 /* If p already has child nodes we must update the sibling pointers.
653 * Otherwise only initialise the left and right pointers of the added
654 * child.
655 */
656 if((l = p->child)) {
657 r = l->right;
658 c->left = l;
659 c->right = r;
660 r->left = c;
661 l->right = c;
662 }
663 else {
664 c->left = c->right = c;
665 }
666
667 p->child = c;
668 c->parent = p;
669}
670
671/* --- replaceNode() ---
672 * Replaces node $old$ and its sub-tree with node $new$ and its sub-tree.
673 */
674void Heap23::replaceNode(Heap23Node *oldNode, Heap23Node *newNode)
675{
676 Heap23Node *parent, *l, *r;
677
678 l = oldNode->left;
679 r = oldNode->right;
680
681 /* If `oldNode' is an only child we only need to initialise the sibling
682 * pointers of the new node. Otherwise we update sibling pointers of other
683 * child nodes.
684 */
685 if(r == oldNode) {
686 newNode->right = newNode->left = newNode;
687 }
688 else {
689 l->right = newNode;
690 r->left = newNode;
691 newNode->left = l;
692 newNode->right = r;
693 }
694
695 /* Update parent pointer of the new node and possibly the child pointer
696 * of the parent node.
697 */
698 parent = oldNode->parent;
699 newNode->parent = parent;
700 if(parent->child == oldNode) parent->child = newNode;
701}
702
703/*--- Heap23 (debugging) ----------------------------------------------------*/
704
705/* --- dumpNopdes() ---
706 * Recursively print the nodes of a 2-3 heap.
707 */
708void Heap23::dumpNodes(Heap23Node *node, size_t level)
709{
710#if defined(TTHEAP_DUMP) && TTHEAP_DUMP > 0
711 Heap23Node *childNode, *partner;
712 size_t childCount;
713
714 /* Print leading whitespace for this level. */
715 for(size_t i = 0; i < level; i++)
716 Rcpp::Rcout << " ";
717
718 Rcpp::Rcout << node->item << "(" << node->key << ")" << std::endl;
719
720 if((childNode = node->child)) {
721 childNode = node->child->right;
722
723 childCount = 0;
724
725 do {
726 dumpNodes(childNode, level+1);
727 if(childNode->dim != childCount) {
728 for(size_t i = 0; i < level+1; i++)
729 Rcpp::Rcout << " ";
730 throw std::error ("error(dim)");
731 }
732 if(childNode->parent != node) {
733 for(size_t i = 0; i < level+1; i++)
734 Rcpp::Rcout << " ";
735 throw std::error ("error(parent)");
736 }
737 if(childNode->key < node->key) {
738 for(size_t i = 0; i < level+1; i++)
739 Rcpp::Rcout << " ";
740 throw std::error ("error(key)");
741 }
742 childNode = childNode->right;
743 childCount++;
744 } while(childNode != node->child->right);
745
746 if(childCount != node->dim && childCount != node->dim + 1) {
747 for(size_t i = 0; i < level; i++)
748 Rcpp::Rcout << " ";
749 throw std::error ("error(childCount)");
750 }
751 }
752 else {
753 if(node->dim != 0) {
754 for(size_t i = 0; i < level; i++)
755 Rcpp::Rcout << " ";
756 throw std::error ("error(dim)");
757 }
758 }
759#endif
760}
761
762/* --- dump() ---
763 * Print a text representation of the heap to the standard output.
764 */
765void Heap23::dump() const
766{
767#if defined(TTHEAP_DUMP) && TTHEAP_DUMP > 0
768 Heap23Node *node;
769
770 Rcpp::Rcout << std::endl << "value = " << treeSum << std::endl;
771 Rcpp::Rcout << "array entries 0..maxTrees =";
772 for(size_t i=0; i<maxTrees; i++) {
773 bool tempval = trees [i] ? 1 : 0;
774 Rcpp::Rcout << " " << tempval;
775 }
776 Rcpp::Rcout << std::endl << std::endl;
777 for(size_t i=0; i<maxTrees; i++) {
778 if((node = trees[i])) {
779 Rcpp::Rcout << "tree " << i << std::endl << std::endl;
780 dumpNodes(node, 0);
781 Rcpp::Rcout << std::endl;
782 }
783 }
784 Rcpp::Rcout.flush ();
785#endif
786}
787
788/*---------------------------------------------------------------------------*/