Distances on Directed Graphs in R
mpadge.tngl.sh
1/* File triheap.h - Trinomial Heap
2 * ----------------------------------------------------------------------------
3 * Mark Padgham, adapted from code by Shane Saunders
4 */
5
6/* This version is implemented using the node-pair pointer structure; that is,
7 * nodes have a partner pointer, so that nodes can be paired. In this
8 * implementation, the child pointer points to the highest dimension child
9 * node.
10 */
11
12#include "triheap.h"
13#include <cmath>
14#if SHOW_trih
15#include <cstdio>
16#include <Rcpp.h>
17#endif
18
19
20#define TRUE 1
21#define FALSE 0
22
23/*--- TriHeap ---------------------------------------------------------------*/
24
25/* --- Constructor ---
26 * Allocates a trinomial heap capable of holding $n$ items.
27 */
28TriHeap::TriHeap(size_t n)
29{
30#if SHOW_trih
31 Rcpp::Rcout << "init, ";
32 Rcpp::Rcout.flush ();
33#endif
34
35 /* The maximum number of nodes and the maximum number of trees allowed. */
36 maxNodes = n;
37 maxTrees = 1 + static_cast <size_t> (log(static_cast <double> (n)) /
38 log (3.0));
39
40 /* Allocate space for an array of pointers to trees, and nodes in the heap.
41 * Initialise all array entries to zero, that is, NULL pointers.
42 */
43 trees = new TriHeapNode *[maxTrees];
44 for (size_t i = 0; i < maxTrees; i++)
45 trees[i] = std::nullptr_t ();
46
47 active = new TriHeapNode *[maxTrees];
48 for (size_t i = 0; i < maxTrees; i++)
49 active[i] = std::nullptr_t ();
50
51 nodes = new TriHeapNode *[n];
52 for (size_t i = 0; i < n; i++)
53 nodes[i] = std::nullptr_t ();
54
55 /* We begin with no nodes in the heap. */
56 itemCount = 0;
57
58 /* The value of the heap helps to keep track of the maximum dimension while
59 * nodes are inserted or deleted.
60 */
61 treeSum = 0;
62
63 /* For experimental purposes, we keep a count of the number of key
64 * comparisons.
65 */
66 compCount = 0;
67}
68
69/* --- Destructor ---
70*/
71TriHeap::~TriHeap()
72{
73#if SHOW_trih
74 Rcpp::Rcout << "free, ";
75 Rcpp::Rcout.flush ();
76#endif
77
78 for (size_t i = 0; i < maxNodes; i++) {
79 delete nodes[i];
80 }
81
82 delete [] nodes;
83 delete [] trees;
84 delete [] active;
85
86#if SHOW_trih
87 Rcpp::Rcout << "free-exited, ";
88 Rcpp::Rcout.flush ();
89#endif
90}
91
92/* --- insert() ---
93 * Inserts a new item, $item$, with assoicated key, $k$, into the heap.
94 */
95void TriHeap::insert(size_t item, double k)
96{
97 TriHeapNode *newNode;
98
99#if SHOW_trih
100 Rcpp::Rcout << "insert, ";
101 Rcpp::Rcout.flush ();
102#endif
103
104 /* Create an initialise the new node. The parent pointer will be set to
105 * NULL by meld().
106 */
107 newNode = new TriHeapNode;
108 newNode->child = std::nullptr_t ();
109 newNode->extra = FALSE;
110 newNode->left = newNode->right = std::nullptr_t ();
111 newNode->partner = std::nullptr_t ();
112
113 newNode->dim = 0;
114 newNode->item = item;
115 newNode->key = k;
116
117 /* Maintain a pointer to item's new node in the heap. */
118 nodes[item] = newNode;
119
120 /* Meld the new node into the heap. */
121 meld(newNode);
122
123 /* Update the heap's node count. */
124 itemCount++;
125
126#if SHOW_trih
127 Rcpp::Rcout << "insert-exited, ";
128 Rcpp::Rcout.flush ();
129#endif
130}
131
132/* --- deleteMin() ---
133 * Deletes and returns the minimum node from the heap.
134 */
135size_t TriHeap::deleteMin()
136{
137 TriHeapNode *minNode, *child, *next, *partner;
138 TriHeapNode *head, *tail, *breakNode, *firstChild;
139 TriHeapNode *l, *parent, *childZero, *childHigher;
140 TriHeapNode *ptr, *nextPartner, *nextParent, *nextFirstChild;
141 TriHeapNode *nextChildZero, *nextChildHigher;
142 double k, k2;
143 size_t d, nextDim, v, item;
144 size_t wasExtra;
145
146#if SHOW_trih
147 Rcpp::Rcout << "deleteMin, ";
148 Rcpp::Rcout.flush ();
149#endif
150
151 /* First we determine the maximum dimension tree in the heap. */
152 v = treeSum;
153 d = 0;
154 while(v) {
155 v = v >> 1;
156 d++;
157 };
158 d--;
159
160 /* Now locate the root node with the smallest key, scanning from the
161 * maximum dimension root position, down to dimension 0 root position.
162 * At the same time, scan active nodes. Note that the maximum dimension of
163 * a active node is one less than the maximum dimension main trunk since we
164 * never get active nodes on a main trunk
165 */
166 minNode = trees[d];
167 k = minNode->key;
168 while(d > 0) {
169 d--;
170 next = trees[d];
171 if(next) {
172 if((k2 = next->key) < k) {
173 k = k2;
174 minNode = next;
175 }
176 compCount++;
177 }
178 next = active[d];
179 if(next) {
180 if((k2 = next->key) < k) {
181 k = k2;
182 minNode = next;
183 }
184 compCount++;
185 }
186 }
187
188#if SHOW_trih
189 Rcpp::Rcout << "on vertex no " << minNode->item << ", ";
190 Rcpp::Rcout.flush ();
191#endif
192
193 /* The resulting break-up depends on whether or not minNode is a root node
194 * or an active node.
195 */
196
197 /* An active node may have been destroyed. */
198 if(minNode->parent) {
199 active[minNode->dim] = std::nullptr_t ();
200 }
201
202 /* During break-up breakNode points to the node that `appears' to have
203 * been removed. Initially this is minNode.
204 */
205 breakNode = minNode;
206 d = breakNode->dim;
207 partner = breakNode->partner;
208 firstChild = breakNode->extra ? partner : breakNode;
209 parent = firstChild->parent; /* parent pointer of the broken node. */
210
211 /* A nodes child pointer always points to its highest dimension child,
212 * so child->right is the smallest dimension. For melding the linked list
213 * starting at child->right, we terminate the circular link with a NULL
214 * pointer.
215 *
216 * For the linked list, only the right pointer is used, so the value of the
217 * left pointer will be left undefined.
218 */
219 child = minNode->child;
220 if(child) {
221 head = child->right;
222 ptr = tail = child;
223
224 /* Nodes in this break up may go from active to inactive. */
225 do {
226 ptr = ptr->right;
227 if(active[ptr->dim] == ptr)
228 active[ptr->dim] = std::nullptr_t ();
229 } while(ptr != tail);
230
231 if(parent) {
232 tail->right = parent;
233 tail = parent;
234 }
235 }
236 else {
237 head = tail = parent;
238 }
239
240
241 /* Propagate break-up through the tree until the main trunk is reached.
242 * Nodes on the main trunk are treated like first child and second child
243 * with no parent. Active nodes never occur as a second child.
244 */
245 if(parent) {
246 childZero = parent->child->right;
247 childHigher = firstChild->right;
248 while(1) {
249 /* At the start of the loop links are not yet updated to account for
250 * the node that appears to have been removed.
251 * parent - points to the parent of the last broken node.
252 * partner - points to the partner of the last broken node.
253 * firstChild - is the first child on the trunk of the broken node.
254 * The smaller of parent and partner is pointed to by the linked list,
255 * Since the linked list is updated in advance.
256 */
257
258
259 /* Make the partner of breakNode `parent's partner. Remember the
260 * old value of the parents partner pointer for the next dimension.
261 * Also remember the value of parent->dim.
262 */
263 nextDim = parent->dim;
264 parent->dim = d;
265 nextPartner = parent->partner;
266 parent->partner = partner;
267 partner->partner = parent;
268
269 /* For the last node pair added to the linked list
270 * (i.e. (parent,partner) or (partner,parent)), we ensure that the
271 * correct node is labelled as extra.
272 */
273 wasExtra = parent->extra;
274 tail->extra = FALSE;
275 tail->partner->extra = TRUE;
276
277 /* Obtain future values of pointer variables now. This is done because
278 * constructing the linked list overwrites pointers that are followed.
279 */
280 if(wasExtra) {
281 nextFirstChild = nextPartner;
282 }
283 else {
284 nextFirstChild = parent;
285 }
286 nextParent = nextFirstChild->parent;
287 nextChildZero = std::nullptr_t ();
288 nextChildHigher = std::nullptr_t ();
289 if(nextParent) {
290 nextChildZero = nextParent->child->right;
291 nextChildHigher = nextFirstChild->right;
292 }
293
294
295 /* Add the child pairs of `parent' that have a greater dimension than
296 * firstChild to the list of broken trunks. Keep a pointer to
297 * parent->child->right for use below.
298 */
299 if(parent->child != firstChild) {
300 ptr = tail;
301 tail->right = childHigher;
302 tail = parent->child;
303
304 /* Nodes in this break up may go from active to inactive.
305 */
306 do {
307 ptr = ptr->right;
308 if(active[ptr->dim] == ptr)
309 active[ptr->dim] = std::nullptr_t ();
310 } while(ptr != tail);
311 }
312
313
314 /* Update the list of children of `parent' to only include those of
315 * lower dimension than firstChild.
316 */
317 if(d) {
318 /* Lower dimension children exist. Note that tail currently points
319 * to the highest dimension child.
320 */
321 l = firstChild->left;
322 l->right = childZero;
323 childZero->left = l;
324 parent->child = l;
325 }
326 else {
327 /* No lower dimension children. */
328 parent->child = std::nullptr_t ();
329 }
330
331
332 /* Now continue break up at 1 dimension higher up by treating `parent'
333 * as the node that has been broken.
334 * Note that on ending, nextParent will be NULL, and we only require
335 * `partner' and d to be updated before exiting.
336 */
337 partner = nextPartner;
338 d = nextDim;
339 if(!nextParent) break;
340 breakNode = parent;
341 parent = nextParent;
342 firstChild = nextFirstChild;
343
344
345 /* Pointers to the dimension zero child of parent, and the next highest
346 * dimension child from firstChild.
347 */
348 childZero = nextChildZero;
349 childHigher = nextChildHigher;
350
351
352 /* Update the linked list in advance to point to the next breakNode,
353 * usually `parent', but possibly `partner' if `partner' is active.
354 */
355 if(wasExtra) {
356
357 /* Since breakNode is not active, `partner' may be, so we may
358 * need to swap the order of `partner' and `parent' in the trunk
359 * resulting from break-up.
360 */
361 if(active[d] == partner) {
362 active[d] = std::nullptr_t ();
363 if(partner->key < parent->key) {
364 /* We make the linked list point to `partner' instead of
365 * `parent', and make parent an extra node.
366 */
367 tail->right = partner;
368 tail = partner;
369 continue; /* Back to start of loop. */
370 }
371 compCount++;
372 }
373 }
374 else {
375 if(active[d] == breakNode)
376 active[d] = std::nullptr_t ();
377 }
378 tail->right = parent;
379 tail = parent;
380
381 }} /* if-while */
382
383
384 /* Break up always propagates up to the main trunk level. After break up
385 * the length the main trunk decreases by one. The current tree position
386 * will become empty unless breakNode has a partner node.
387 */
388 if(partner) {
389 partner->partner = std::nullptr_t ();
390
391 if(partner->extra) {
392 partner->extra = FALSE;
393 partner->parent = std::nullptr_t ();
394 partner->left = partner->right = partner;
395 trees[d] = partner;
396 }
397 }
398 else {
399 trees[d] = std::nullptr_t ();
400 treeSum -= (1 << d);
401 }
402 itemCount--;
403
404 /* Meld the linked list of trunks resulting from break-up into the main
405 * trunk level of the heap.
406 */
407 if(head) {
408 tail->right = std::nullptr_t ();
409 meld(head);
410 }
411
412 /* Record the vertex no to return. */
413 item = minNode->item;
414
415 /* Delete the old minimum node. */
416 nodes[item] = std::nullptr_t ();
417 delete minNode;
418
419#if SHOW_trih
420 Rcpp::Rcout << "deleteMin-exited, ";
421 Rcpp::Rcout.flush ();
422#endif
423
424 return item;
425}
426
427/* --- decreaseKey() ---
428 * This mthod decreases the key of the node corresponding to $item$ to
429 * $newValue$. It is up to the user to ensure that $newValue$ is in-fact less
430 * than or equal to the current value.
431 */
432void TriHeap::decreaseKey(size_t item, double newValue)
433{
434 TriHeapNode *v, *v2, *w, *w2, *p, *above, *partner, *activeNode;
435 TriHeapNode *l, *r, *lowChild, *highChild, *ptr;
436 size_t d;
437
438#if SHOW_trih
439 Rcpp::Rcout << "decreaseKey on vn = " << item << "(" << newValue << "), ";
440 Rcpp::Rcout.flush ();
441#endif
442
443
444 /* Pointer v points to the decreased node. */
445 v = nodes[item];
446 v->key = newValue;
447 d = v->dim; /* dimension */
448
449
450 /* This loop allows rearrangement to propagate to higher dimensions. */
451 while(1) {
452
453
454 v2 = v->partner; /* partner */
455 if(v->extra) {
456 p = v2->parent; /* parent */
457 above = v2;
458 }
459 else {
460 above = p = v->parent; /* parent */
461 }
462
463
464 /* Determine if rearrangement is necessary */
465
466 /* If v is a root node, then rearrangement is not necessary. */
467 if(!above) return;
468
469 if(p) {
470 /* Non-main trunk. */
471
472 activeNode = active[d];
473
474 if(v->extra) {
475 /* v is a second child. If necessary, we swap with the first
476 * child to maintain the correct ordering.
477 */
478
479 compCount++;
480 if(v->key < above->key) {
481 /* swap */
482 v->extra = FALSE;
483 v2->extra = TRUE;
484 replaceChild(v2, v);
485
486 /* If v2 is inconsistent try promotion. By checking if
487 * v2 was an active node we can avoid key comparison, since
488 * v2 can only be inconsistent if it is active.
489 */
490 if(activeNode == v2) {
491
492 compCount++;
493 if(v2->key < p->key) goto promote; /* see below */
494
495 /* At this point, v2 is active but is consistent, so v
496 * will replace it as the active node, regardless of
497 * whether v is inconsistent or not.
498 */
499 active[d] = v;
500 return;
501 }
502
503 /* If there is some other active node for this dimension.
504 * try rearrangement with it.
505 */
506 if(activeNode) goto rearrange; /* see below */
507
508 /* Otherwise we can make v an active node, regardless of
509 * whether v is inconsistent or not.
510 */
511 active[d] = v;
512 return;
513 }
514
515 /* Otherwise don't swap */
516
517 /* We may need to promote, but only if both v and v2 are
518 * inconsistent. By checking if v2 is already active, we
519 * can avoid key comparison since v2 won't be inconsistent
520 * if it is not active.
521 */
522 if(activeNode == v2) {
523 /* If v is inconsistent, then v2 is also, so promote. */
524 compCount++;
525 if(v->key < p->key) {
526 /* Swap v and v2. */
527 v2 = v;
528 v = above;
529 goto promote;
530 }
531
532 /* Otherwise promotion is not necessary. */
533 return;
534 }
535
536 /* At this point, promotion is not necessary. */
537 return;
538 }
539
540 /* Otherwise, v is a first child. */
541
542 if(activeNode) {
543 /* If v is already the active node leave it. */
544 if(activeNode == v) return;
545
546 /* Otherwise some other node is active, so try
547 * rearrangement.
548 */
549 goto rearrange; /* see below */
550 }
551
552 /* At this point, there is no active node, so make v the active
553 * node, as it is possibly inconsistent.
554 */
555 active[d] = v;
556 return;
557 }
558
559 /* Otherwise, v is the second node on a main trunk. */
560
561 compCount++;
562 if(v->key < above->key) {
563 /* If v is smaller, we swap it with the first child
564 * (i.e. the root node) to maintain the correct ordering.
565 */
566 v->extra = FALSE;
567 v2->extra = TRUE;
568 v->parent = std::nullptr_t ();
569 v->left = v->right = v;
570 trees[d] = v;
571 return;
572 }
573
574 /* Otherwise, no rearrangement is required since heap order is still
575 * satisfied.
576 */
577 return;
578
579 /*------------------------------*/
580rearrange:
581
582 /* At this stage v is a first child and possibly inconsistent, but not
583 * an active node. Node v2 is consistent. Rearrangement occurs with
584 * the currently active node for this dimension, which at this point
585 * must exist.
586 */
587
588
589 /* The current active node and its partner are w and w2 respectively.
590 */
591 w = activeNode;
592 w2 = w->partner;
593
594 /* The final rearrangement always pairs v with w and v2 with w2. */
595 v->partner = w; w->partner = v;
596 v2->partner = w2; w2->partner = v2;
597
598 /* Determine the ordering in the rearrangement. */
599 compCount++;
600 if(v2->key < w2->key) {
601 /* Make the (1st child, 2nd child) pair (v2,w2), replacing (v,v2).
602 */
603 v2->extra = FALSE;
604 replaceChild(v, v2);
605
606 compCount++;
607 if(v->key < w->key) {
608 /* Make the pair (v,w), replacing (w,w2). */
609 w->extra = TRUE;
610 replaceChild(w,v);
611
612 compCount++;
613 if(w->key < v->parent->key) {
614 /* Both v and w are inconsistent, continue with
615 * promotion. Update v2, and p;
616 */
617 v2 = w;
618 p = v->parent;
619 goto promote; /* see below */
620 }
621
622 /* Otherwise, although w is active, it is not inconsistent.
623 * In this case we do not do promotion. Instead, v replaces w
624 * as the active node, although v may not be inconsistent.
625 * (Avoids a key comparison to check if v is consistent)
626 */
627 active[d] = v;
628 return;
629 }
630 else {
631 /* Make the pair (w,v), replacing (w,w2). */
632 v->extra = TRUE;
633
634 compCount++;
635 if(v->key < w->parent->key) {
636 /* Both v and w are inconsistent, so continue with
637 * promotion. Update v, v2, and p.
638 */
639 v2 = v;
640 v = w;
641 p = w->parent;
642 goto promote; /* see below */
643 }
644
645 /* Only w is possibly inconsistent, so it remains the active
646 * node.
647 */
648 return;
649 }
650 }
651 else {
652 /* Make the pair (w2,v2), replacing (w,w2). */
653 w2->extra = FALSE;
654 replaceChild(w, w2);
655
656 compCount++;
657 if(v->key < w->key) {
658 /* Make the pair (v,w), replacing (v,v2). */
659 w->extra = TRUE;
660
661 compCount++;
662 if(w->key < v->parent->key) {
663 /* Both v and w are inconsistent, so continue with
664 * promotion. Update v2.
665 */
666 v2 = w;
667 goto promote; /* see below */
668 }
669
670 /* Only v is possibly inconsistent, so it becomes the active
671 * node.
672 */
673 active[d] = v;
674 return;
675 }
676 else {
677 /* Make the pair (w,v), replacing (v,v2). */
678 v->extra = TRUE;
679 replaceChild(v,w);
680
681 compCount++;
682 if(v->key < w->parent->key) {
683 /* Both v and w are inconsistent, so continue with
684 * promotion. Update v, v2.
685 */
686 v2 = v;
687 v = w;
688 goto promote;
689 }
690
691 /* Only w is possibly inconsistent, so it remains the
692 * active node.
693 */
694 return;
695 }
696 }
697
698 /*------------------------------*/
699promote:
700 /* Promotion Code. Reached if the loop has not exited. Uses variables
701 * p, v, v2, and d. Must ensure that on the trunk (p,v,v2), both v and
702 * v2 are inconsistent nodes, so that (v,v2,p) will give heap
703 * ordering.
704 */
705
706 /* First we make v2 a child node of v. */
707 v2->extra = FALSE;
708 addChild(v, v2);
709
710 /* Then v replaces p. Any child nodes of p that have a higher
711 * dimension than v will become child nodes of v. Only child nodes of
712 * lower dimension than v will be left under p.
713 */
714
715 v->dim = p->dim;
716 partner = v->partner = p->partner;
717 highChild = p->child;
718
719 if(d) {
720 /* v has lower dimension siblings. */
721 l = p->child = v->left;
722 r = highChild->right;
723 l->right = r;
724 r->left = l;
725 }
726 else {
727 /* v was an only child. */
728 p->child = std::nullptr_t ();
729 }
730
731 if(highChild != v) {
732 /* v has higher dimension siblings. Add them to the list of v's
733 * children, and update their parent pointer. Note that v
734 * currently has at least one child, since v2 was made a child
735 * of v.
736 */
737 lowChild = v->right;
738 l = v->child;
739 r = v->child->right;
740 l->right = lowChild;
741 lowChild->left = l;
742 r->left = highChild;
743 highChild->right = r;
744
745 v->child = highChild;
746
747 ptr = v;
748 do {
749 ptr = ptr->right;
750 ptr->parent = v;
751 } while(ptr != highChild);
752 }
753
754 /* partner may be NULL if p is a root node, so don't update
755 * partner->partner yet. See update below.
756 */
757
758 /* If p was an extra node no further pointer updates are required. */
759
760 /* Further pointer updates are only needed if p is a first child or a
761 * root node.
762 */
763 if(!p->extra) {
764 if(p->parent) {
765 /* p is non-root node and a first child. */
766 partner->partner = v;
767 replaceChild(p, v);
768 }
769 else {
770 /* p is a root node, so update the tree pointer. */
771 if(partner) partner->partner = v;
772 trees[p->dim] = v;
773 v->left = v->right = v;
774 v->parent = std::nullptr_t ();
775 }
776
777 /* p will become an extra node, see below. */
778 p->extra = TRUE;
779 }
780 else {
781 /* If p was an extra node then v becomes an extra node. */
782 partner->partner = v;
783 v->extra = TRUE;
784 }
785
786 /* Finally, make p the partner of node v2 (i.e. the 2nd child of
787 * node v).
788 */
789 p->dim = d;
790 v2->partner = p;
791 p->partner = v2;
792
793 /* The result of promotion is to release the active node. */
794 active[d] = std::nullptr_t ();
795
796 d = v->dim; /* next dimension, moving up */
797
798 /* If p was an active node, it no longer is, and v takes its place
799 * as an active node.
800 */
801 if(active[d] == p) {
802 active[d] = v;
803 return;
804 }
805
806 /* Rearrangement continues at a higher dimension, since the position of
807 * v has moved up, meaning v may be an inconsistent node.
808 */
809 }
810
811#if SHOW_trih
812 Rcpp::Rcout << "decrease_key-exited, ";
813 Rcpp::Rcout.flush ();
814#endif
815
816}
817
818/*--- TriHeap (private methods) ---------------------------------------------*/
819
820/* --- meld() ---
821 * Melds the linked list of trees pointed to by *treeList into the heap.
822 * This function uses the `right' sibling pointer of nodes to traverse the
823 * linked list from lower dimension nodes to higher dimension nodes.
824 * It expects the last nodes `right' pointer to be NULL.
825 */
826void TriHeap::meld(TriHeapNode *treeList)
827{
828 TriHeapNode *next, *addTree;
829 TriHeapNode *carryTree;
830 size_t d;
831
832#if SHOW_trih
833 Rcpp::Rcout << "meld - ";
834 Rcpp::Rcout.flush ();
835#endif
836
837 /* addTree points to the tree to be merged. */
838 addTree = treeList;
839
840 carryTree = std::nullptr_t ();
841
842 do {
843 /* addTree() gets merged into the heap, and also carryTree if one
844 * exists from a previous merge.
845 */
846
847 /* Keep a pointer to the next tree and remove sibling and parent links
848 * from the current tree. The dimension of the next tree is always
849 * one greater than the dimension of the previous tree, so this merging
850 * is like an addition of two ternary numbers.
851 *
852 * Note that if addTree is NULL and the loop has not exited, then
853 * there is only a carryTree to be merged, so treat it like addTree.
854 */
855 next = std::nullptr_t ();
856 if(addTree) {
857 next = addTree->right;
858 addTree->right = addTree->left = addTree;
859 addTree->parent = std::nullptr_t ();
860 }
861 else {
862 addTree = carryTree;
863 carryTree = std::nullptr_t ();
864 }
865
866#if SHOW_trih
867 Rcpp::Rcout << addTree->item << ", ";
868 Rcpp::Rcout.flush ();
869#endif
870
871 /* First we merge addTree with carryTree, if there is one. Note that
872 * carryTree contains only one node in its main trunk, and addTree
873 * has at most two, so the result is at most one 3-node trunk, which is
874 * treated as a 1-node main trunk one dimension higher up.
875 */
876 if(carryTree) {
877 compCount += merge(&addTree, &carryTree);
878 }
879
880 /* After the merge, if addTree is NULL, then the resulting tree
881 * pointed to by carryTree carries to higher entry, so we do not need
882 * to merge anything into the existing main trunk.
883 * If addTree is not NULL we add it to the existing main trunk.
884 */
885 if(addTree) {
886 d = addTree->dim;
887 if(trees[d]) {
888 /* Nodes already in this main trunk position, so merge. */
889 compCount += merge(&trees[d], &addTree);
890 if(!trees[d]) treeSum -= (1 << d);
891 carryTree = addTree;
892 }
893 else {
894 /* No nodes in this main trunk position, so use addTree. */
895 trees[d] = addTree;
896 treeSum += (1 << d);
897 }
898 }
899
900 /* Obtain a pointer to the next tree to add. */
901 addTree = next;
902
903 /* We continue if there is still a node in the list to be merged, or
904 * a carry tree remains to be merged.
905 */
906 } while(addTree || carryTree);
907
908
909#if SHOW_trih
910 Rcpp::Rcout << "meld-exited, ";
911 Rcpp::Rcout.flush ();
912#endif
913
914}
915
916/*--- TriHeap (node manipulation methods) -----------------------------------*/
917
918/* --- merge() ---
919 * Merges the two trunks pointed to by *a and *b, returning the
920 * sum trunk through `a' and any carry tree through `b'.
921 * When this function is used, both parameters `a' and `b' refer to either
922 * a 1-node or 2-node trunk.
923 *
924 * Returns the number of key comparisons used.
925 */
926size_t TriHeap::merge(TriHeapNode **a, TriHeapNode **b)
927{
928 TriHeapNode *tree, *nextTree, *other, *nextOther;
929 size_t c;
930
931 /* Number of comparisons. */
932 c = 0;
933
934 /* `tree' always points to the node with the lowest key.
935 * To begin with, `tree' points to the smaller head node, and `other'
936 * points to the head node of the other trunk.
937 */
938 if((*a)->key <= (*b)->key) {
939 tree = (*a);
940 other = (*b);
941 }
942 else {
943 tree = (*b);
944 other = (*a);
945 }
946 c++;
947
948 /* nextTree points to the next node on the trunk that `tree' is the head
949 * of (if there is another node).
950 * nextOther points to the next node on the trunk that `other' is the head
951 * of (if there is another node).
952 */
953 nextTree = tree->partner;
954 nextOther = other->partner;
955
956 /* The merging depends on the existence of nodes and the values of keys. */
957 if(!nextTree) {
958 /* nextTree does not exist, so we simply make `other' the child of
959 * `tree'. If nextOther exist the resulting 3-node trunk is a carry
960 * tree.
961 */
962
963 if(nextOther) {
964 addChild(tree, other);
965 tree->dim++;
966 *a = std::nullptr_t ();
967 *b = tree;
968 }
969 else {
970 tree->partner = other;
971 other->partner = tree;
972 other->extra = TRUE;
973
974 *a = tree;
975 *b = std::nullptr_t ();
976 }
977 }
978 else if(!nextOther) {
979 /* nextTree exists but nextOther does not, so the linked order of
980 * nextTree and `other' in the resulting 3-node trunk depends on the
981 * values of keys. The resulting 3-node trunk becomes a carry tree.
982 */
983
984 tree->partner = std::nullptr_t ();
985 other->partner = nextTree;
986 nextTree->partner = other;
987
988 if(other->key < nextTree->key) {
989 addChild(tree, other);
990 }
991 else {
992 nextTree->extra = FALSE;
993 other->extra = TRUE;
994 addChild(tree, nextTree);
995 }
996
997 tree->dim++;
998
999 c++;
1000 *a = std::nullptr_t ();
1001 *b = tree;
1002 }
1003 else {
1004 /* Otherwise, both nextTree and nextOther exist. The result consists
1005 * of a 1 node trunk plus the 3-node trunk which becomes a carry tree.
1006 * We two trunks are made up as (tree, other, nextOther)
1007 * and (nextTree). This uses no key comparisons.
1008 */
1009
1010 tree->partner = std::nullptr_t ();
1011 nextTree->partner = std::nullptr_t ();
1012 nextTree->extra = FALSE;
1013 nextTree->left = nextTree->right = nextTree;
1014 nextTree->parent = std::nullptr_t ();
1015
1016 addChild(tree, other);
1017
1018 tree->dim++;
1019
1020 *a = nextTree; *b = tree;
1021 }
1022
1023 return c;
1024}
1025
1026/* --- addChild() ---
1027 * Adds the child node pointeed to by $c$ to the parent node pointed to by $p$.
1028 * The user should ensure that the correct dimension child is being added.
1029 */
1030void TriHeap::addChild(TriHeapNode *p, TriHeapNode *c)
1031{
1032 TriHeapNode *l, *r;
1033
1034 /* If $p$ already has child nodes we must update the sibling pointers.
1035 * Otherwise only initialise the left and right pointers of the added
1036 * child.
1037 */
1038 if((l = p->child)) {
1039 r = l->right;
1040 c->left = l;
1041 c->right = r;
1042 r->left = c;
1043 l->right = c;
1044 }
1045 else {
1046 c->left = c->right = c;
1047 }
1048
1049 p->child = c;
1050 c->parent = p;
1051}
1052
1053
1054/* --- replaceChild() ---
1055 * Replaces child node pointed to by $oldNode$ and its associated sub-tree with the
1056 * child node pointed to by $newNode$ and its associated sub-tree.
1057 */
1058void TriHeap::replaceChild(TriHeapNode *oldNode, TriHeapNode *newNode)
1059{
1060 TriHeapNode *parent, *l, *r;
1061
1062 r = oldNode->right;
1063
1064 /* If $oldNode$ is an only child we only need to initialise the sibling
1065 * pointers of the new node. Otherwise we update sibling pointers of other
1066 * child nodes.
1067 */
1068 if(r == oldNode) {
1069 newNode->right = newNode->left = newNode;
1070 }
1071 else {
1072 l = oldNode->left;
1073 l->right = newNode;
1074 r->left = newNode;
1075 newNode->left = l;
1076 newNode->right = r;
1077 }
1078
1079 /* Update parent pointer of the new node and possibly the child pointer
1080 * of the parent node.
1081 */
1082 parent = oldNode->parent;
1083 newNode->parent = parent;
1084 if(parent->child == oldNode) parent->child = newNode;
1085}
1086
1087/*--- TriHeap (debugging) ---------------------------------------------------*/
1088
1089/* --- dumpNodes() ---
1090 * Recursively print the nodes of a trinomial heap.
1091 */
1092//TriHeapNode **Active;
1093void TriHeap::dumpNodes(TriHeapNode *node, size_t level)
1094{
1095#if SHOW_trih
1096 TriHeapNode *childNode, *partnerNode;
1097 size_t childCount;
1098
1099 /* Print leading whitespace for this level. */
1100 for (size_t i = 0; i < level; i++)
1101 Rcpp::Rcout << " ";
1102
1103 Rcpp::Rcout << node->item << "(" << node->key << ")" << std::endl;
1104
1105 if((childNode = node->child)) {
1106 childNode = node->child->right;
1107
1108 childCount = 0;
1109
1110 do {
1111 dumpNodes(childNode, level+1);
1112 if(childNode->dim != childCount) {
1113 for (size_t i = 0; i < level+1; i++)
1114 Rcpp::Rcout << " ";
1115 throw std::runtime_error ("error(dim)");
1116 }
1117 if(childNode->parent != node) {
1118 for (size_t i = 0; i < level+1; i++)
1119 Rcpp::Rcout << " ";
1120 throw std::runtime_error ("error(parent)");
1121 }
1122 if(Active[childNode->dim] != childNode &&
1123 childNode->key < node->key) {
1124 for (size_t i = 0; i < level; i++)
1125 Rcpp::Rcout << " ";
1126 throw std::runtime_error ("error(key)");
1127 }
1128 childNode = childNode->right;
1129 childCount++;
1130 } while(childNode != node->child->right);
1131
1132 if(childCount != node->dim) {
1133 for (size_t i = 0; i < level; i++)
1134 Rcpp::Rcout << " ";
1135 throw std::runtime_error ("error(childCount)");
1136 }
1137 }
1138 else {
1139 if(node->dim != 0) {
1140 for (size_t i = 0; i < level; i++)
1141 Rcpp::Rcout << " ";
1142 throw std::runtime_error ("error(dim)");
1143 }
1144 }
1145
1146 if((partner=node->partner)) {
1147 if(node->extra==partner->extra) {
1148 for (size_t i = 0; i < level; i++)
1149 Rcpp::Rcout << " ";
1150 Rcpp::Rcout << partner->item;
1151 throw std::runtime_error (" - error(extra?)");
1152 }
1153 if(partner->extra) {
1154 if(partner->dim != node->dim) {
1155 for (size_t i = 0; i < level; i++)
1156 Rcpp::Rcout << " ";
1157 Rcpp::Rcout << partner->item;
1158 throw std::runtime_error (" - error(dim)");
1159 }
1160
1161 dumpNodes(partner, level);
1162 if(partner->key < node->key) {
1163 for (size_t i = 0; i < level; i++)
1164 Rcpp::Rcout << " ";
1165 throw std::runtime_error ("error(key)");
1166 }
1167 }
1168 }
1169 else if(node->parent) {
1170 for (size_t i = 0; i < level; i++)
1171 Rcpp::Rcout << " ";
1172 throw std::runtime_error ("error(no partner)");
1173 }
1174#endif
1175}
1176
1177/* --- dump() ---
1178 * Print out a trinomial heap.
1179 */
1180void TriHeap::dump() const
1181{
1182#if SHOW_trih
1183 TriHeapNode *node;
1184
1185 Active = active;
1186
1187 Rcpp::Rcout << std::endl << "value = " << treeSum << std::endl;
1188 Rcpp::Rcout << "array entries 0..maxTrees =";
1189 for (size_t i=0; i<maxTrees; i++) {
1190 bool tempval = trees [i] ? 1 : 0
1191 Rcpp::Rcout << " " << tempval;
1192 }
1193 Rcpp::Rcout << std::endl << "active nodes =";
1194 for (size_t i=0; i<maxTrees-1; i++) {
1195 if(active[i]) {
1196 Rcpp::Rcout << " " << active [i]->item;
1197 if(active[i]->dim != i) {
1198 throw std::runtime_error ("-error(dim)");
1199 }
1200 if(active[i]->extra) {
1201 throw std::runtime_error ("-error(extra)");
1202 }
1203 if(!active[i]->parent) {
1204 throw std::runtime_error ("-error(parent)");
1205 }
1206 }
1207 else {
1208 Rcpp::Rcout << " _";
1209 }
1210 }
1211 Rcpp::Rcout << std::endl << std::endl;
1212 for (size_t i=0; i<maxTrees; i++) {
1213 if((node = trees[i])) {
1214 Rcpp::Rcout << "tree " << i << std::endl << std::endl;
1215 dumpNodes(node, 0);
1216 Rcpp::Rcout << std::endl;
1217 }
1218 }
1219 Rcpp::Rcout.flush ();
1220#endif
1221}
1222/*---------------------------------------------------------------------------*/