diff options
Diffstat (limited to 'blossom5-v2.05.src/LCA.h')
| -rw-r--r-- | blossom5-v2.05.src/LCA.h | 297 |
1 files changed, 297 insertions, 0 deletions
diff --git a/blossom5-v2.05.src/LCA.h b/blossom5-v2.05.src/LCA.h new file mode 100644 index 0000000..3426f7e --- /dev/null +++ b/blossom5-v2.05.src/LCA.h @@ -0,0 +1,297 @@ +/* + +LCA.h - Imlementation of LCA data structure described in + +O. Berkman, U. Vishkin: Recursive Star-Tree Parallel Data Structure. SIAM J. Comput. 22(2): 221-242 (1993) + +It takes O(n \log n) space (and O(n \log n) time for construction). The memory requirement can be +reduced to O(n) by defining LCA_BLOCKS. The implementation then follows the approach described in + +J. Fischer, V. Heun: Theoretical and Practical Improvements on the RMQ-Problem, with Applications to LCA and LCE. + Proceedings of the 17th Annual Symposium on Combinatorial Pattern Matching + (CPM'06), Lecture Notes in Computer Science 4009, 36-48, Springer-Verlag, 2006. + +(although in-block queries are computed naively in O(K) worst-case time, where K is the size of the block, +rather than as described by Fischer and Heun). + +Written by Vladimir Kolmogorov, vnk@adastral.ucl.ac.uk. + +// example for +// 4 +// / \ +// 2 3 +// / \ +// 0 1 +// (Note: the ordering of nodes is in the tree preorder!!!) + + LCATree* lca = new LCATree(5); + char A0, A1, A2, A3, A4; + + lca->Add(&A0, &A2); + lca->Add(&A1, &A2); + lca->Add(&A2, &A4); + lca->Add(&A3, &A4); + lca->AddRoot(&A4); + + int result = lca->GetLCA(1, 3); // should be 4 + delete lca; + + +*/ + +#ifndef GNAKDLATHJSTHAJSRNAKSJDA +#define GNAKDLATHJSTHAJSRNAKSJDA + +#include <assert.h> +#include <stdio.h> +#include <stdlib.h> +#include <string.h> + +//#define LCA_BLOCKS + +class LCATree +{ +public: + typedef void* NodeId; // can be any type, e.g. int. (The code checks NodeId's only for equalities.) + typedef int PreorderId; + + LCATree(int node_num_max); + ~LCATree(); + + // construct tree. Nodes must be added in the tree preorder!!! + // First call returns 0, second returns 1, and so on. + PreorderId Add(NodeId i, NodeId i_parent); + PreorderId AddRoot(NodeId i); // completes tree construction + + PreorderId GetLCA(PreorderId i, PreorderId j); + // Let i0=i, j0=j be the input nodes, and let r = LCA(i0,j0). + // This function sets i and j to be the immediate children of r + // such that i is a descendant of i0 and j is a descendant of j0. + // There must hold i0!=r, j0!=r. + void GetPenultimateNodes(PreorderId& i, PreorderId& j); + + +////////////////////////////////////////////////////////////////////////// +private: + int n, n_max, K, k_max; + int** array; + + NodeId* buf0; + int* buf1; + NodeId* parent_current; + int* child_current; + + int* parents; + + int _GetLCA(int i, int j); // same as GetLCA, but assumes that i<j + int GetLCADirect(int i, int j); +}; + +inline LCATree::LCATree(int _n_max) : n_max(_n_max), array(NULL) +{ +#ifdef LCA_BLOCKS + K = -2; + n = n_max; + while (n > 0) { K ++; n /= 2; } + if (K < 1) K = 1; +#else + K = 1; + n = 0; +#endif + parents = new int[n_max]; + buf0 = new NodeId[n_max]; + buf1 = new int[n_max]; + parent_current = buf0; + child_current = buf1; +} + +inline LCATree::~LCATree() +{ + int k; + delete [] parents; + if (buf0) delete [] buf0; + if (buf1) delete [] buf1; + if (array) + { + for (k=1; k<=k_max; k++) delete [] array[k]; + delete [] array; + } +} + +inline LCATree::PreorderId LCATree::Add(NodeId i, NodeId i_parent) +{ + assert(n < n_max); + + if (n == 0) + { + *parent_current = i; + *(++ parent_current) = i_parent; + parents[0] = -1; + } + else + { + if (i == *parent_current) + { + int c = *child_current --; + while ( 1 ) + { + int c_next = parents[c]; + parents[c] = n; + if (c_next < 0) break; + c = c_next; + } + parent_current --; + } + if (i_parent == *parent_current) parents[n] = *child_current; + else + { + *(++ parent_current) = i_parent; + parents[n] = -1; + child_current ++; + } + } + *child_current = n; + return n ++; +} + + +inline LCATree::PreorderId LCATree::AddRoot(NodeId i) +{ + assert(n < n_max); + + if (n > 0) + { + if (i != *parent_current || parent_current != buf0+1) + { + printf("Error in LCATree construction: wrong sequence of calls!\n"); + exit(1); + } + int c = *child_current --; + while ( 1 ) + { + int c_next = parents[c]; + parents[c] = n; + if (c_next < 0) break; + c = c_next; + } + child_current ++; + } + parents[n++] = -1; + + delete [] buf0; + buf0 = NULL; + delete [] buf1; + buf1 = NULL; + + // initialize array + int b, k = 1, block_num = (n-1)/K+1; + if (block_num < 3) return n-1; + int d = (block_num-1)/4; + while (d) { k ++; d >>= 1; } + k_max = k; + + array = new int*[k_max+1]; + array[0] = parents; + for (k=1, d=2; k<=k_max; k++, d*=2) + { + array[k] = new int[block_num-d]; + if (k == 1) + { + for (b=0; b<block_num-d; b++) array[1][b] = GetLCADirect((b+1)*K-1, (b+d)*K); + } + else + { + for (b=0; b<block_num-d; b++) + { + int i = array[k-1][b]; + int j = array[k-1][b+d/2]; + if (i < j) i = j; + j = array[1][b+d/2-1]; + array[k][b] = (i > j) ? i : j; + } + } + } + + return n-1; +} + +/////////////////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////////////////// + +inline int LCATree::GetLCADirect(int i, int j) +{ + while (i < j) i = parents[i]; + return i; +} + + +inline int LCATree::_GetLCA(int i, int j) +{ +#ifdef LCA_BLOCKS + + int bi = i/K, bj = j/K; + if (bi == bj) return GetLCADirect(i, j); + int i_last = (bi+1)*K-1, j_first = bj*K; + i = GetLCADirect(i, i_last); + j = GetLCADirect(j_first, j); + if (i < j) i = j; + // set j = LCA(i_last, j_first) + if (j_first - i_last == 1) j = parents[i_last]; + else + { + int k = 1, d = (bj-bi)/4; + while (d) { k ++; d >>= 1; } + int diff = 1<<k; + //assert(bi+diff <= bj && bi+diff>bj-diff); + j = (array[k][bi] > array[k][bj-diff]) ? array[k][bi] : array[k][bj-diff]; + } + return (i > j) ? i : j; + +#else + + if (j == i) return i; + + int k = 0, d = (j-i)/2; + while (d) { k ++; d >>= 1; } + int diff = 1<<k; + //assert(i+diff <= j && i+diff>j-diff); + return (array[k][i] > array[k][j-diff]) ? array[k][i] : array[k][j-diff]; + +#endif +} + +/////////////////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////////////////// + +inline LCATree::PreorderId LCATree::GetLCA(PreorderId i, PreorderId j) +{ + if (i > j) { PreorderId k = i; i = j; j = k; } + return _GetLCA(i, j); +} + +inline void LCATree::GetPenultimateNodes(PreorderId& _i, PreorderId& _j) +{ + int i, j, d, swap; + if (_i < _j) { i = _i; j = _j; swap = 0; } + else { i = _j; j = _i; swap = 1; } + int r = _GetLCA(i, j); + assert(i!=r && j!=r); + while (parents[i] != r) + { + int i0 = parents[i]; + d = (j - i0)/2; + while ( (i=_GetLCA(i0, i0+d)) == r ) d /= 2; + } + while (parents[j] != r) + { + int j0 = parents[j]; + d = (r - j0)/2; + while ( (j=_GetLCA(j0, j0+d)) == r ) d /= 2; + } + if (swap == 0) { _i = i; _j = j; } + else { _j = i; _i = j; } +} + +#endif |