Something went wrong. Try again.
Distances on Directed Graphs in R
Something went wrong. Try again.
123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335133613371338133913401341134213431344134513461347134813491350135113521353135413551356135713581359136013611362136313641365136613671368136913701371137213731374137513761377137813791380138113821383138413851386138713881389139013911392139313941395139613971398139914001401140214031404140514061407140814091410141114121413141414151416141714181419142014211422142314241425142614271428142914301431143214331434143514361437143814391440144114421443144414451446144714481449145014511452145314541455145614571458145914601461146214631464146514661467146814691470147114721473147414751476147714781479148014811482/*** Trinomial Heap Implementation ***//* * Mark Padgham, adapted from code by Shane Saunders */
/* This version is implemented using the node-pair pointer structure; that is, * nodes have a partner pointer, so that nodes can be paired. In this * implementation, a nodes child pointer points to its highest dimension child * node. */
#include "triheap_ext.h"#include <cmath>//#if SHOW_trih_ext#include <cstdio>#include <stdexcept>//#endif
#define TRUE 1#define FALSE 0
/*--- TriHeapExt (public methods) -------------------------------------------*/
/* --- Constructor --- * creates and returns a pointer to a trinomial heap. Argument * maxNodes specifies the maximum number of nodes the heap can contain. */TriHeapExt::TriHeapExt(size_t n){#if SHOW_trih_ext Rcpp::Rcout << "init, "; Rcpp::Rcout.flush ();#endif
/* The maximum number of nodes and the maximum number of trees allowed. */ maxNodes = n; maxTrees = 1 + static_cast <size_t> (log(static_cast <double> (n)) / log (3.0));
/* The tolerance trigger of the heap. That is, t+1; the number used for * detecting when we are over the tolerance and need to cleanup active * nodes. */ activeLimit = maxTrees;
/* Allocate space for an array of pointers to trees, and nodes in the heap. * Initialise all array entries to zero, that is, NULL pointers. */ trees = new TriHeapExtNode *[maxTrees]; for (size_t i = 0; i < maxTrees; i++) trees[i] = nullptr;
nodes = new TriHeapExtNode *[n]; for (size_t i = 0; i < n; i++) nodes[i] = nullptr;
/* Allocate space for: * - an unordered array of pointers to active nodes. * - an array of pointers to queues of active nodes. * - an array of pointers to entries in the candidates queue. */ activeNodes = new TriHeapExtNode *[activeLimit]; for (size_t i = 0; i < activeLimit; i++) activeNodes[i] = nullptr;
activeQueues = new ActiveItem *[activeLimit-1]; for (size_t i = 0; i < activeLimit-1; i++) activeQueues[i] = nullptr;
candidateItems = new CandidateItem *[activeLimit-1]; for (size_t i = 0; i < activeLimit-1; i++) candidateItems[i] = nullptr;
/* We begin with no nodes in the heap. */ itemCount = 0; activeCount = 0; candQueueHead = nullptr;
/* The value of the heap helps to keep track of the maximum dimension while * nodes are inserted or deleted. */ treeSum = 0;
/* For experimental purposes, we keep a count of the number of key * comparisons. */ compCount = 0;}
/* --- Destructor ---*/TriHeapExt::~TriHeapExt(){#if SHOW_trih_ext Rcpp::Rcout << "free, "; Rcpp::Rcout.flush ();#endif
for (size_t i = 0; i < maxNodes; i++) { delete nodes[i]; }
delete [] nodes; delete [] trees; delete [] activeNodes; delete [] activeQueues; delete [] candidateItems;}
/* --- insert() --- * inserts item $item$ with associated key $k$ into the heap. */void TriHeapExt::insert(size_t item, double k){ TriHeapExtNode *newNode;
#if SHOW_trih_ext Rcpp::Rcout << "insert, "; Rcpp::Rcout.flush ();#endif
/* Create an initialise the new node. The parent pointer will be set to * NULL by meld(). */ newNode = new TriHeapExtNode; newNode->child = nullptr; newNode->extra = FALSE; newNode->left = newNode->right = nullptr; newNode->partner = nullptr;
newNode->activeEntry = nullptr;
newNode->dim = 0; newNode->item = item; newNode->key = k;
/* Maintain a pointer to item's new node in the heap. */ nodes[item] = newNode;
/* Meld the new node into the heap. */ meld(newNode);
/* Update the heap's node count. */ itemCount++;}
/* --- deleteMin() --- * Deletes and returns the minimum item from the heap. */size_t TriHeapExt::deleteMin(){ static TriHeapExtNode meldListHeader;
TriHeapExtNode *minNode, *child, *next, *partner; TriHeapExtNode *tail, *breakNode, *firstChild; TriHeapExtNode *l, *parent, *childZero, *childHigher; TriHeapExtNode *nextPartner, *nextParent, *nextFirstChild; TriHeapExtNode *node, *nodePartner; TriHeapExtNode *nextChildZero = nullptr, *nextChildHigher = nullptr; double k, k2; size_t d, nextDim, v, item; size_t wasExtra;
#if SHOW_trih_ext dump(); Rcpp::Rcout << "deleteMin, "; Rcpp::Rcout.flush ();#endif
/* First we determine the maximum dimension tree in the heap. */ v = treeSum; d = 0; while(v) { v = v >> 1; d++; } d--;
/* Now locate the root node with the smallest key, scanning from the * maximum dimension root position, down to dimension 0 root position. * At the same time, scan active nodes. Note that the maximum dimension of * a active node is one less than the maximum dimension main trunk since we * never get active nodes on a main trunk. */ minNode = trees[d]; k = minNode->key; while(d > 0) { d--; next = trees[d]; if(next) { compCount++; if((k2 = next->key) < k) { k = k2; minNode = next; } }
} size_t i = activeCount; while(i > 0) { i--; next = activeNodes[i]; compCount++; if((k2 = next->key) < k) { k = k2; minNode = next; } }
#if SHOW_trih_ext Rcpp::Rcout << "on vertex no " << minNode->item; Rcpp::Rcout.flush ();#endif
/* * The node-pair tree containing the minimum node must be broken up into * a linked list of node-pair sub-trees to be melded back into the heap. * This linked list first includes any children trees broken of the minimum * node. * */ tail = &meldListHeader;
/* A nodes child pointer always points to its highest dimension child, * so child->right is the smallest dimension. For melding the linked list * starting at child->right, we terminate the circular link with a NULL * pointer. * * For the linked list, only the right pointer is used, so the value of the * left pointer will be left undefined. */ child = minNode->child; if(child) { tail->right = child->right; node = tail = child;
/* Nodes in this break up may go from active to inactive. */ do { node = node->right; nodePartner = node->partner;
if(node->activeEntry) { deactivate(node); if(nodePartner->activeEntry) deactivate(nodePartner); } /* Note that we do not need to check a second child for activeness * if the first child was not active. */
} while(node != tail); }
/* * The linked list also includes higher trees that become broken up in * cases where a non-root (i.e. active) minimum node was removed. * */
/* During higher break-up breakNode points to the node that `appears' to * have been removed. Initially this is minNode. */ d = minNode->dim; breakNode = minNode; partner = breakNode->partner; firstChild = breakNode->extra ? partner : breakNode; parent = firstChild->parent; /* parent pointer of the broken node. */
/* If the minimum node was active then deactivate it, and propagate * break-up through the tree until the main trunk is reached. * Nodes on the main trunk are treated like first child and second child * with no parent. */ if(parent) {
deactivate(minNode);
childZero = parent->child->right; childHigher = firstChild->right;
/* If `partner' is active then deactivate it and determine the order of * `partner' and `parent' in the trunk resulting from break-up. Otherwise * parent is the first node on the trunk. */ if(partner->activeEntry) { deactivate(partner); compCount++; if(partner->key < parent->key) { tail->right = partner; tail = partner; } else { tail->right = parent; tail = parent; } } else { tail->right = parent; tail = parent; }
while(1) { /* At the start of the loop links are not yet updated to account for * the node that appears to have been removed. * parent - points to the parent of the last broken node. * partner - points to the partner of the last broken node. * firstChild - is the first child on the trunk of the broken node. * The smaller of parent and partner is pointed to by the linked list, * Since the linked list is updated in advance. */
/* Make the partner of breakNode `parent's partner. Remember the * old value of the parents partner pointer for the next dimension. * Also remember the value of parent->dim. The code has been written * this way because certain pointer updates overwrite pointers that * are needed later on. * Note, if we must deactivate parent, this must be done before * changing its dimension. */ if(parent->activeEntry) deactivate(parent); nextDim = parent->dim; parent->dim = d; nextPartner = parent->partner; parent->partner = partner; partner->partner = parent;
/* For the last node pair added to the linked list * (i.e. (parent,partner) or (partner,parent)), we ensure that the * correct node is labelled as extra. */ wasExtra = parent->extra; tail->extra = FALSE; tail->partner->extra = TRUE;
/* Obtain future values of pointer variables now. This is done because * constructing the linked list overwrites pointers that are followed. */ if(wasExtra) { nextFirstChild = nextPartner; } else { nextFirstChild = parent; } nextParent = nextFirstChild->parent; if(nextParent) { nextChildZero = nextParent->child->right; nextChildHigher = nextFirstChild->right; }
/* Add the child pairs of `parent' that have a greater dimension than * firstChild to the list of broken trunks. Keep a pointer to * parent->child->right for use below. */ if(parent->child != firstChild) { node = tail; tail->right = childHigher; tail = parent->child;
/* Nodes in this break up may go from active to inactive. */ do { node = node->right; nodePartner = node->partner; if(node->activeEntry) { deactivate(node); if(nodePartner->activeEntry) deactivate(nodePartner); } } while(node != tail); }
/* Update the list of children of `parent' to only include those of * lower dimension than firstChild. */ if(d) { /* Lower dimension children exist. Note that tail currently points * to the highest dimension child. */ l = firstChild->left; l->right = childZero; childZero->left = l; parent->child = l; } else { /* No lower dimension children. */ parent->child = nullptr; }
/* Now continue break up at 1 dimension higher up by treating `parent' * as the node that has been broken. * Note that on ending, nextParent will be NULL, and we only require * `partner' and d to be updated before exiting. */ partner = nextPartner; d = nextDim; if(!nextParent) break; breakNode = parent; parent = nextParent; firstChild = nextFirstChild;
/* Pointers to the dimension zero child of parent, and the next highest * dimension child from firstChild. */ childZero = nextChildZero; childHigher = nextChildHigher;
/* Update the linked list in advance to point to the next breakNode, * usually `parent', but possibly `partner' if `partner' is active. */
/* `partner' may be active, so we may need to swap the order of * `partner' and `parent' in the trunk resulting from break-up. */ if(partner->activeEntry) { deactivate(partner); compCount++; if(partner->key < parent->key) { /* We make the linked list point to `partner' instead of * `parent', and make parent an extra node. */ tail->right = partner; tail = partner; continue; /* Back to start of loop. */ } }
tail->right = parent; tail = parent;
}} /* if-while */
/* Terminate the list of trees to be melded. */ tail->right = nullptr;
/* Break up always propagates up to the main trunk level. After break up * the length the main trunk decreases by one. The current tree position * will become empty unless breakNode has a partner node. */ if(partner) { partner->partner = nullptr;
if(partner->extra) { partner->extra = FALSE; partner->parent = nullptr; partner->left = partner->right = partner; trees[d] = partner; } } else { trees[d] = nullptr; treeSum -= (1 << d); } itemCount--;
/* Meld the linked list of trunks resulting from break-up into the main * trunk level of the heap. */ if(meldListHeader.right) meld(meldListHeader.right);
/* Record the vertex no to return. */ item = minNode->item;
/* Delete the old minimum node. */ nodes[item] = nullptr; delete minNode;
return item;}
/* --- decreaseKey() --- * Decrease the key of item $item$ in the heap to the value $newValue$. It is * up to the user to ensure that newValue is in-fact less than or equal to the * current value. */void TriHeapExt::decreaseKey(size_t item, double newValue){ TriHeapExtNode *v, *v2, *w, *w2, *p, *above, *partner, *activeNode; TriHeapExtNode *l, *r, *lowChild, *highChild, *node; ActiveItem *activeEntry; size_t d;
#if SHOW_trih_ext dump(); Rcpp::Rcout << "decreaseKey on vn = " << item << " (" << newValue << ")"; Rcpp::Rcout.flush ();#endif
/* Pointer v points to the decreased node. */ v = nodes[item]; v->key = newValue; d = v->dim; /* dimension */
v2 = v->partner; /* partner */ if(v->extra) { p = v2->parent; /* parent */ above = v2; } else { above = p = v->parent; /* parent */ }
/* Determine if rearrangement is necessary */
/* If v is a root node, then rearrangement is not necessary. */ if(!above) return;
/* If v is already active then rearrangement is not necessary. */ if(v->activeEntry) { /* If v is a second child and its key has became less than its partner, * which will also be active, then swap to maintain the correct * ordering. */ if(v->extra && v->key < above->key) { v->extra = FALSE; v2->extra = TRUE; replaceChild(v2, v); } return; }
if(p) { /* Non-main trunk. */
if(v->extra) { /* v is a second child. If necessary, we swap with the first * child to maintain the correct ordering. */
compCount++; if(v->key < above->key) { /* swap */ v->extra = FALSE; v2->extra = TRUE; replaceChild(v2, v); } else { /* If v remains as a second child, then only activate it if * the first child is active. */ if(v2->activeEntry) activate(v);
if(activeCount == activeLimit) goto rearrange; /* see below */
return; } }
/* At this point, v is a first child. */
activate(v);
/* If the number of active nodes is at tolerance level we must perform * some rearrangement. */ if(activeCount == activeLimit) goto rearrange; /* see below */
/* Otherwise it is okay to exit. */ return;
}
/* Otherwise, v is the second node on a main trunk. */
compCount++; if(v->key < above->key) { /* If v is smaller, we swap it with the first child (i.e. the root * node) to maintain the correct ordering. */ v->extra = FALSE; v2->extra = TRUE; v->parent = nullptr; v->left = v->right = v; trees[d] = v; return; }
/* Otherwise, no rearrangement is required since heap order is still * satisfied. */ return;
/*------------------------------*/rearrange:
/* Get a candidate for rearrangement. */ d = candQueueHead->dim; activeEntry = activeQueues[d];
activeNode = activeEntry->node; if(activeNode->extra) { v = activeNode->partner; v2 = activeNode; } else { v = activeNode; v2 = activeNode->partner; } p = v->parent;
/* If we have two active nodes on the same trunk. */ if(v2->activeEntry) { deactivate(v2);
/* If the 2nd child of these is less than the parent, then * do promotion. */ if(v2->key < v->parent->key) goto promote;
return; }
/* Try the second trunk. */ activeNode = activeEntry->next->node; if(activeNode->extra) { w = activeNode->partner; w2 = activeNode; } else { w = activeNode; w2 = activeNode->partner; }
/* If we have two active nodes on the same trunk. */ if(w2->activeEntry) { deactivate(w2);
/* If the 2nd child of these is less than the parent, then * do promotion. */ if(w2->key < w->parent->key) { v = w; v2 = w2; p = w->parent; goto promote; }
return; }
/* The final rearrangement always pairs v with w and v2 with w2. */ v->partner = w; w->partner = v; v2->partner = w2; w2->partner = v2;
/* Determine the ordering in the rearrangement. */ compCount++; if(v2->key < w2->key) { /* Make the (1st child, 2nd child) pair (v2,w2), replacing (v,v2). */ v2->extra = FALSE; replaceChild(v, v2);
compCount++; if(v->key < w->key) { /* Make the pair (v,w), replacing (w,w2). */ w->extra = TRUE; replaceChild(w,v); deactivate(w);
compCount++; if(w->key < v->parent->key) { /* Both v and w are inconsistent, continue with promotion. * Update v2, and p; */ v2 = w; p = v->parent; goto promote; /* see below */ }
/* Otherwise, although w was active, it is not inconsistent. In * this case we do not do promotion. Instead, only v remains an * active node, although v may not be inconsistent. * (Avoids a key comparison to check if v is consistent) */ return; } else { /* Make the pair (w,v), replacing (w,w2). */ v->extra = TRUE; deactivate(v);
compCount++; if(v->key < w->parent->key) { /* Both v and w are inconsistent, so continue with promotion. * Update v, v2, and p. */ v2 = v; v = w; p = w->parent; goto promote; /* see below */ }
/* Only w is possibly inconsistent, so it remains an active node. */ return; } } else { /* Make the pair (w2,v2), replacing (w,w2). */ w2->extra = FALSE; replaceChild(w, w2);
compCount++; if(v->key < w->key) { /* Make the pair (v,w), replacing (v,v2). */ w->extra = TRUE; deactivate(w);
compCount++; if(w->key < v->parent->key) { /* Both v and w are inconsistent, so continue with promotion. * Update v2. */ v2 = w; goto promote; /* see below */ }
/* Only v is possibly inconsistent, so it remains an active node. */ return; } else { /* Make the pair (w,v), replacing (v,v2). */ v->extra = TRUE; replaceChild(v,w); deactivate(v);
compCount++; if(v->key < w->parent->key) { /* Both v and w are inconsistent, so continue with promotion. * Update v, v2. */ v2 = v; v = w; goto promote; }
/* Only w is possibly inconsistent, so it remains an active node. */ return; } }
/*------------------------------*/promote: /* Promotion Code. Reached if the loop has not exited. Uses variables * p, v, v2, and d. Must ensure that on the trunk (p,v,v2), both v and v2 * are inconsistent nodes, so that (v,v2,p) will give heap ordering. * Node v should still be active, and node v2 should have been deactivated. * Node v will later become active at a higher dimension. */ deactivate(v);
/* First we make v2 a child node of v. */ v2->extra = FALSE; addChild(v, v2);
/* Then v replaces p. Any child nodes of p that have a higher dimension * than v will become child nodes of v. Only child nodes of lower * dimension than v will be left under p. */
v->dim = p->dim; partner = v->partner = p->partner; highChild = p->child;
if(d) { /* v has lower dimension siblings. */ l = p->child = v->left; r = highChild->right; l->right = r; r->left = l; } else { /* v was an only child. */ p->child = nullptr; }
if(highChild != v) { /* v has higher dimension siblings. Add them to the list of v's * children, and update their parent pointer. Note that v currently * has at least one child, since v2 was made a child of v. */ lowChild = v->right; l = v->child; r = v->child->right; l->right = lowChild; lowChild->left = l; r->left = highChild; highChild->right = r;
v->child = highChild;
node = v; do { node = node->right; node->parent = v; } while(node != highChild); }
/* partner may be NULL if p is a root node, so don't update * partner->partner yet. See update below. */
/* If p was an extra node no further pointer updates are required. */
/* Further pointer updates are only needed if p is a first child or a root * node. */ if(!p->extra) { if(p->parent) { /* p is non-root node and a first child. */ partner->partner = v; replaceChild(p, v); } else { /* p is a root node, so update the tree pointer. */ if(partner) partner->partner = v; trees[p->dim] = v; v->left = v->right = v; v->parent = nullptr; }
/* p will become an extra node, see below. */ p->extra = TRUE; } else { /* If p was an extra node then v becomes an extra node. */ partner->partner = v; v->extra = TRUE; }
/* Finally, make p the partner of node v2 (i.e. the 2nd child of node v). * Note we cant update the dimension of p yet until we deactivate it. */ v2->partner = p; p->partner = v2;
/* If v is a second child we may need to swap it with the first child to * maintain correct ordering. If v remains as a second child, we do not * make it active, unless the first child is also active. */ if(v->extra) { v2 = v->partner;
compCount++; if(v->key < v2->key) { /* swap */ v->extra = FALSE; v2->extra = TRUE;
/* Note there is a special case, where a main trunk is reached. */ if(v2->parent) { /* Non-main trunk */ replaceChild(v2, v); } else { /* Main trunk: v replaces v2 as the root node. */ v->extra = FALSE; v2->extra = TRUE; v->parent = nullptr; v->left = v->right = v; trees[v->dim] = v;
/* Don't make v active. * If necessary, deactivate p (which v replaced). */ if(p->activeEntry) deactivate(p); p->dim = d; return; } } else if(!v2->activeEntry) { /* Don't make v active. This is always the case for main trunks. * If necessary, deactivate p (which v replaced). */ if(p->activeEntry) deactivate(p); p->dim = d; return; } } else { /* Don't swap, and if at main trunk level, don't make v active either. */ if(!v->parent) { /* If necessary, deactivate p (which v replaced). */ if(p->activeEntry) deactivate(p); p->dim = d; return; } }
/* The result of the above code never puts active nodes on a main trunk. * If at main trunk level, the function will have exited above. */
/* If p was an active node, it no longer is, and v takes its place as an * active node. */ if(p->activeEntry) { replaceActive(p, v); } else { /* Otherwise make v active. */ activate(v); } p->dim = d;}
/*--- TriHeapExt (private methods) ------------------------------------------*/
/* --- meld() --- * Melds the linked list of trees pointed to by $treeList$ into the heap * pointed to by $h$. This function uses the 'right' sibling pointer of nodes * to traverse the linked list from lower dimension nodes to higher * dimension nodes. It expects the last nodes `right' pointer to be NULL. */void TriHeapExt::meld(TriHeapExtNode *treeList){ TriHeapExtNode *next, *addTree; TriHeapExtNode *carryTree; size_t d;
#if SHOW_trih_ext Rcpp::Rcout << "meld - "; Rcpp::Rcout.flush ();#endif
/* addTree points to the tree to be merged. */ addTree = treeList;
carryTree = nullptr;
do { /* addTree() gets merged into the heap, and also carryTree if one * exists from a previous merge. */
/* Keep a pointer to the next tree and remove sibling and parent links * from the current tree. The dimension of the next tree is always * one greater than the dimension of the previous tree, so this merging * is like an addition of two ternary numbers. * * Note that if addTree is NULL and the loop has not exited, then * there is only a carryTree to be merged, so treat it like addTree. */ next = nullptr; if(addTree) { next = addTree->right; addTree->right = addTree->left = addTree; addTree->parent = nullptr; } else { addTree = carryTree; carryTree = nullptr; }
#if SHOW_trih_ext Rcpp::Rcout << addTree->item << ", "; Rcpp::Rcout.flush ();#endif
/* First we merge addTree with carryTree, if there is one. Note that * carryTree contains only one node in its main trunk, and addTree * has at most two, so the result is at most one 3-node trunk, which is * treated as a 1-node main trunk one dimension higher up. */ if(carryTree) { compCount += merge(&addTree, &carryTree); }
/* After the merge, if addTree is NULL, then the resulting tree * pointed to by carryTree carries to higher entry, so we do not need * to merge anything into the existing main trunk. * If addTree is not NULL we add it to the existing main trunk. */ if(addTree) { d = addTree->dim; if(trees[d]) { /* Nodes already in this main trunk position, so merge. */ compCount += merge(&trees[d], &addTree); if(!trees[d]) treeSum -= (1 << d); carryTree = addTree; } else { /* No nodes in this main trunk position, so use addTree. */ trees[d] = addTree; treeSum += (1 << d); } }
/* Obtain a pointer to the next tree to add. */ addTree = next;
/* We continue if there is still a node in the list to be merged, or * a carry tree remains to be merged. */ } while(addTree || carryTree);}
/* --- activate() --- * Add an inactive node to the queue of active nodes. It is up to the user of * this function to ensure that the node was not already active. */void TriHeapExt::activate(TriHeapExtNode *node){ size_t d; ActiveItem *activeEntry, *first, *last; CandidateItem *candidate, *lastC;
/* Note that we maintain doubly linked lists. This allows active items * and candidates to be removed in O(1) time. The circular linking allows * us to access the last element in the list in O(1) time. */
/* Add n to the array of active nodes. */ size_t i = activeCount++; activeNodes[i] = node;
/* Create an entry for n in the list of pointers to active nodes. */ activeEntry = new ActiveItem; activeEntry->node = node; activeEntry->position = i; node->activeEntry = activeEntry;
/* Insertion depends on whether the list is empty or not. */ d = node->dim;
if((first = activeQueues[d])) { /* At least one item already. */ last = first->prev; last->next = activeEntry; activeEntry->prev = last; activeEntry->next = first; first->prev = activeEntry;
/* If there was originally one item, but now two, then insert a new * entry into the candidates queue. */ if(first == last) { candidate = new CandidateItem; candidate->dim = d; candidateItems[d] = candidate;
if(candQueueHead) { /* At least one candidate is already in the queue. */ lastC = candQueueHead->prev; lastC->next = candidate; candidate->prev = lastC; candidate->next = candQueueHead; candQueueHead->prev = candidate; } else { /* Inserting into empty queue. */ candQueueHead = candidate; candidate->next = candidate->prev = candidate; } } } else { /* empty */ activeQueues[d] = activeEntry; activeEntry->next = activeEntry->prev = activeEntry; }}
/* --- deactivate() --- * Remove an active node from the set of active nodes, and make it inactive. */void TriHeapExt::deactivate(TriHeapExtNode *node){ ActiveItem *activeEntry, *next, *prev, *first, *second; TriHeapExtNode *topActive; CandidateItem *candidate, *nextCandidate, *prevCandidate; size_t d;
/* Obtain pointers to the corresponding entry in the active node structure * and remove the node from the array of active nodes. */
activeEntry = node->activeEntry;
size_t i = --activeCount; topActive = activeNodes[i]; activeNodes[activeEntry->position] = topActive; topActive->activeEntry->position = activeEntry->position; activeNodes[i] = nullptr;
node->activeEntry = nullptr; d = node->dim; first = activeQueues[d]; second = first->next;
/* Update the list according to the amount and position of existing list * items. This is a circular doubly linked list. */ if(second != first) { /* There are at least two items in the list. */ prev = activeEntry->prev; next = activeEntry->next;
/* May need to change pointer to first item. */ if(activeEntry == first) { activeQueues[d] = second; }
/* If there were only two nodes, we need to remove the candidate entry * from the candidates queue. */ if(second->next == first) { /* remove candidate. */ candidate = candidateItems[d]; candidateItems[d] = nullptr; nextCandidate = candidate->next; prevCandidate = candidate->prev;
/* May need to change pointer to the first item. */ if(nextCandidate == candidate) { candQueueHead = nullptr; } else { if(candQueueHead == candidate) { candQueueHead = nextCandidate; } prevCandidate->next = nextCandidate; nextCandidate->prev = prevCandidate; }
delete candidate; }
prev->next = next; next->prev = prev; } else { /* There is only one item in the list. */ activeQueues[d] = nullptr; }
delete activeEntry;}
/* --- replaceActive() --- * Replaces one active node with another by simply updating the entry to point * to the other node. This is much quicker than performing deactivate() and * activate(). */void TriHeapExt::replaceActive( TriHeapExtNode *oldNode, TriHeapExtNode *newNode){ ActiveItem *activeEntry;
activeEntry = oldNode->activeEntry; newNode->activeEntry = activeEntry; oldNode->activeEntry = nullptr; activeNodes[activeEntry->position] = newNode;
activeEntry->node = newNode;}
/*--- TriHeapExt (node manipulation methods) --------------------------------*/
/* --- merge() --- * Merges the two trunks pointed to by *a and *b, returning the sum trunk * through `a' and any carry tree through `b'. When this function is used, * both parameters `a' and `b' refer to either a 1-node or 2-node trunk. * * Returns the number of key comparisons used. */size_t TriHeapExt::merge(TriHeapExtNode **a, TriHeapExtNode **b){ TriHeapExtNode *tree, *nextTree, *other, *nextOther; size_t c;
/* Number of comparisons. */ c = 0;
/* `tree' always points to the node with the lowest key. * To begin with, `tree' points to the smaller head node, and `other' * points to the head node of the other trunk. */ if((*a)->key <= (*b)->key) { tree = (*a); other = (*b); } else { tree = (*b); other = (*a); } c++;
/* nextTree points to the next node on the trunk that `tree' is the head * of (if there is another node). * nextOther points to the next node on the trunk that `other' is the head * of (if there is another node). */ nextTree = tree->partner; nextOther = other->partner;
/* The merging depends on the existence of nodes and the values of keys. */ if(!nextTree) { /* nextTree does not exist, so we simply make `other' the child of * `tree'. If nextOther exist the resulting 3-node trunk is a carry * tree. */
if(nextOther) { addChild(tree, other); tree->dim++; *a = nullptr; *b = tree; } else { tree->partner = other; other->partner = tree; other->extra = TRUE;
*a = tree; *b = nullptr; } } else if(!nextOther) { /* nextTree exists but nextOther does not, so the linked order of * nextTree and `other' in the resulting 3-node trunk depends on the * values of keys. The resulting 3-node trunk becomes a carry tree. */
tree->partner = nullptr; other->partner = nextTree; nextTree->partner = other;
if(other->key < nextTree->key) { addChild(tree, other); } else { nextTree->extra = FALSE; other->extra = TRUE; addChild(tree, nextTree); }
tree->dim++;
c++; *a = nullptr; *b = tree; } else { /* Otherwise, both nextTree and nextOther exist. The result consists * of a 1 node trunk plus the 3-node trunk which becomes a carry tree. * We two trunks are made up as (tree, other, nextOther) * and (nextTree). This uses no key comparisons. */
tree->partner = nullptr; nextTree->partner = nullptr; nextTree->extra = FALSE; nextTree->left = nextTree->right = nextTree; nextTree->parent = nullptr;
addChild(tree, other);
tree->dim++;
*a = nextTree; *b = tree; }
return c;}
/* --- addChild() --- * Adds a new child node pointed to by $c$, to the parent node pointed to by * $p$. The user should ensure that the correct dimension child is being * added. */void TriHeapExt::addChild(TriHeapExtNode *p, TriHeapExtNode *c){ TriHeapExtNode *l, *r;
/* If p already has child nodes we must update the sibling pointers. * Otherwise only initialise the left and right pointers of the added * child. */ if((l = p->child)) { r = l->right; c->left = l; c->right = r; r->left = c; l->right = c; } else { c->left = c->right = c; }
p->child = c; c->parent = p;}
/* --- replaceChild() --- * Replaces child node pointed to by $old$ and its sub-tree with the child node * pointed to by $new$ and its sub-tree. */void TriHeapExt::replaceChild( TriHeapExtNode *oldNode, TriHeapExtNode *newNode){ TriHeapExtNode *parent, *l, *r;
r = oldNode->right;
/* If `oldNode' is an only child we only need to initialise the sibling * pointers of the new node. Otherwise we update sibling pointers of other * child nodes. */ if(r == oldNode) { newNode->right = newNode->left = newNode; } else { l = oldNode->left; l->right = newNode; r->left = newNode; newNode->left = l; newNode->right = r; }
/* Update parent pointer of the new node and possibly the child pointer * of the parent node. */ parent = oldNode->parent; newNode->parent = parent; if(parent->child == oldNode) parent->child = newNode;}
/*--- TriHeapExt (debugging) ------------------------------------------------*/
/* --- dumpNodes() ---- * Recursively print the nodes of a trinomial heap. */void TriHeapExt::dumpNodes(TriHeapExtNode *node, size_t level){ //#if SHOW_trih_ext TriHeapExtNode *childNode, *partner; size_t childCount;
/* Print leading whitespace for this level. */#if SHOW_trih_ext for (size_t i = 0; i < level; i++) Rcpp::Rcout << " ";
Rcpp::Rcout << node->item << "(" << node->key << ")"; Rcpp::Rcout.flush (); if(node->activeEntry) Rcpp::Rcout << '*'; Rcpp::Rcout << std::endl;#endif
if((childNode = node->child)) { childNode = node->child->right;
childCount = 0;
do { dumpNodes(childNode, level+1); if(childNode->dim != childCount) { throw std::runtime_error ("error(dim)"); } if(childNode->parent != node) { throw std::runtime_error ("error(parent)"); } if(childNode->activeEntry == nullptr && childNode->key < node->key) { throw std::runtime_error ("error(key)"); } childNode = childNode->right; childCount++; } while(childNode != node->child->right);
if(childCount != node->dim) { throw std::runtime_error ("error(childCount)"); } } else { if(node->dim != 0) { throw std::runtime_error ("error(dim)"); } }
if((partner=node->partner)) { if(node->extra==partner->extra) { throw std::runtime_error ("error(extra?)"); } if(partner->extra) { if(partner->dim != node->dim) { throw std::runtime_error ("error(dim)"); } if(partner->activeEntry && ! node->activeEntry) { throw std::runtime_error ("error(active)"); }
dumpNodes(partner, level); if(partner->key < node->key) { throw std::runtime_error ("error(key)"); } } } else if(node->parent) { throw std::runtime_error ("error(no partner)"); }
if(node->activeEntry) { if(node->activeEntry->node != node) { throw std::runtime_error ("error(active entry wrong)"); } } //#endif}
/* --- dump() --- * Print out a trinomial heap. */void TriHeapExt::dump() const{#if SHOW_trih_ext TriHeapExtNode *node; TriHeapExtNode *firstChild; size_t c;
Rcpp::Rcout << std::endl << "value = " << treeSum << std::endl << "array entries 0..maxTrees ="; for (size_t i=0; i<maxTrees; i++) { bool is1or0 = trees [i] ? 1 : 0; Rcpp::Rcout << " " << is1or0; } Rcpp::Rcout << std::endl << "active nodes ="; for (size_t i=0; i<activeCount; i++) { if((node = activeNodes[i])) { Rcpp::Rcout << " " << node->item; firstChild = node->extra ? node->partner : node; if(!firstChild->parent) { throw std::runtime_error ("error(main trunk)"); } if(!node->activeEntry) { throw std::runtime_error ("error(inactive)"); } if(node->activeEntry->node != node) { throw std::runtime_error ("error(active entry wrong)"); } if(activeNodes[node->activeEntry->position] != node) { throw std::runtime_error ("error(active entry index wrong)"); }
c = 0; for(size_t j = i+1; j < activeCount; j++) { if(activeNodes[j] == node) c++; } if(c) { throw std::runtime_error ("error(repeated)"); } } else { throw std::runtime_error ("error(missing active mode)"); } } Rcpp::Rcout << std::endl << std::endl; for (size_t i=0; i<maxTrees; i++) { if((node = trees[i])) { Rcpp::Rcout << "tree " << i << std::endl << std::endl; dumpNodes(node, 0); Rcpp::Rcout << std::endl; } }#endif}
/*---------------------------------------------------------------------------*/