diff --git a/inst/include/TreeTools/ClusterTable.h b/inst/include/TreeTools/ClusterTable.h
index 2bb770d22..ccdae4dc8 100644
--- a/inst/include/TreeTools/ClusterTable.h
+++ b/inst/include/TreeTools/ClusterTable.h
@@ -4,14 +4,20 @@
 #include <array> /* for array */
 #include <bitset> /* for bitset */
 #include <vector> /* for vector */
+#include <sstream> /* for ostringstream */
 #include <Rcpp/Lightest>
 #include "assert.h" /* for ASSERT */
-#include "types.h" /* for int16 */
+#include "types.h" /* for int32 */
 #include "root_tree.h" /* for root_on_node */
 
 #define UNINIT -999
 #define INF TreeTools::INTX_MAX
 
+#define CT_ASSERT_CAN_PUSH()                                   \
+  ASSERT(static_cast<size_t>(Spos + CT_STACK_SIZE) <= S.size())
+
+#define CT_ASSERT_CAN_POP() ASSERT(Spos >= CT_STACK_SIZE)
+
 #define CT_PUSH(a, b, c, d)                                      \
   S[Spos++] = (a);                                               \
   S[Spos++] = (b);                                               \
@@ -36,8 +42,17 @@ const int_fast32_t CT_MAX_LEAVES = 16383;
 
 namespace TreeTools {
 
-  constexpr int_fast32_t ct_max_leaves = 16383;
-  constexpr int_fast32_t ct_stack_size = 4;
+  // Use stack allocation for trees up to this size for optimal performance
+  // Using 'const' instead of 'constexpr' to ensure external linkage for packages
+  // that depend on TreeTools (e.g., TreeDist). constexpr can cause ODR issues.
+  const int_fast32_t ct_stack_threshold = 8192;
+  // Old hard limit, kept for backward compatibility with TreeDist 2.9.2
+  // NOTE: This constant is deprecated - new code should use ct_max_leaves_heap
+  // External packages may still reference this constant for compatibility
+  const int_fast32_t ct_max_leaves = 16383;
+  // New increased limit with heap allocation
+  const int_fast32_t ct_max_leaves_heap = 100000;
+  const int_fast32_t ct_stack_size = 4;
   
   template <typename T>
   inline void resize_uninitialized(std::vector<T>& v, std::size_t n) {
@@ -50,33 +65,77 @@ namespace TreeTools {
     }
   }
   
+  // Dynamic bitset that uses stack allocation for small sizes, heap for large
+  class DynamicBitset {
+  private:
+    std::bitset<ct_stack_threshold + 1> stack_bits;
+    std::vector<unsigned char> heap_bits;
+    bool use_heap;
+    std::size_t size_;
+    
+  public:
+    DynamicBitset() : use_heap(false), size_(0) {}
+    
+    void resize(std::size_t n) {
+      size_ = n;
+      if (n > ct_stack_threshold) {
+        use_heap = true;
+        heap_bits.assign(n, 0);
+      } else {
+        use_heap = false;
+        stack_bits.reset();
+      }
+    }
+    
+    bool operator[](std::size_t idx) const {
+      return use_heap ? heap_bits[idx] != 0 : stack_bits[idx];
+    }
+    
+    void set(std::size_t idx, bool value = true) {
+      if (use_heap) {
+        heap_bits[idx] = value ? 1 : 0;
+      } else {
+        stack_bits[idx] = value;
+      }
+    }
+    
+    void reset() {
+      if (use_heap) {
+        std::fill(heap_bits.begin(), heap_bits.end(), 0);
+      } else {
+        stack_bits.reset();
+      }
+    }
+    
+    std::size_t size() const { return size_; }
+  };
+  
   class ClusterTable {
     
     struct ClusterRow {
-      int16 L;
-      int16 R;
+      int32 L;
+      int32 R;
     };
     
-    int16 n_edge;
-    int16 n_internal;
-    int16 n_leaves;
-    int16 n_shared = 0;
-    int16 enumeration = 0;
-    int16 v_j;
-    int16 Tlen;
-    int16 Tlen_short;
-    int16 X_ROWS;
-    std::vector<int16> internal_label;
-    int16 *internal_label_ptr = nullptr;
-    std::vector<int16> leftmost_leaf;
-    std::vector<int16> T;
-    int16* T_ptr = nullptr;
-    std::vector<int16> visited_nth;
+    int32 n_edge;
+    int32 n_internal;
+    int32 n_leaves;
+    int32 n_shared = 0;
+    int32 enumeration = 0;
+    int32 v_j;
+    int32 Tlen;
+    int32 Tlen_short;
+    int32 X_ROWS;
+    std::vector<int32> internal_label;
+    int32 *internal_label_ptr = nullptr;
+    std::vector<int32> leftmost_leaf;
+    std::vector<int32> T;
+    int32* T_ptr = nullptr;
+    std::vector<int32> visited_nth;
     std::vector<ClusterRow> x_rows;
-    // Using bitset; can obtain a ~1% speedup using vector of ULLs
-    // Retaining slower code as easier to read.
-    // See branch ct-xswitch for implementation
-    std::bitset<ct_max_leaves + 1> Xswitch;
+    // Dynamic bitset that uses stack allocation for small trees,
+    // heap allocation for large trees
+    DynamicBitset Xswitch;
     // Track number of set switches (excluding index 0)
     std::size_t xswitch_set_count = 0;
     
@@ -84,36 +143,36 @@ namespace TreeTools {
   public:
     ClusterTable(Rcpp::List); // i.e. PREPARE(T)
 
-    [[nodiscard]] inline bool is_leaf(const int16 v) noexcept {
+    [[nodiscard]] inline bool is_leaf(const int32 v) noexcept {
       return v <= n_leaves;
     }
     
     // Required by TreeDist 2.9.2
-    // TODO Remove in later version, to prefer is_leaf(int16 v)
-    [[nodiscard]] inline bool is_leaf(const int16 *v) noexcept {
+    // TODO Remove in later version, to prefer is_leaf(int32 v)
+    [[nodiscard]] inline bool is_leaf(const int32 *v) noexcept {
       return *v <= n_leaves;
     }
 
-    [[nodiscard]] inline const int16 edges() noexcept {
+    [[nodiscard]] inline const int32 edges() noexcept {
       return n_edge;
     }
 
-    [[nodiscard]] inline const int16 leaves() noexcept {
+    [[nodiscard]] inline const int32 leaves() noexcept {
       return n_leaves;
     }
 
-    inline void ENTER(int16 v, int16 w) noexcept {
+    inline void ENTER(int32 v, int32 w) noexcept {
       *T_ptr = v;
       ++T_ptr;
       *T_ptr = w;
       ++T_ptr;
     }
 
-    [[nodiscard]] inline int16 N() noexcept {
+    [[nodiscard]] inline int32 N() noexcept {
       return n_leaves;
     }
 
-    [[nodiscard]] inline int16 M() noexcept {
+    [[nodiscard]] inline int32 M() noexcept {
       return n_internal;
     }
 
@@ -123,12 +182,12 @@ namespace TreeTools {
       T_ptr = T.data();
     }
 
-    inline void READT(int16 *v, int16 *w) {
+    inline void READT(int32 *v, int32 *w) {
       *v = *T_ptr++;
       *w = *T_ptr++;
     }
 
-    inline void NVERTEX(int16 *v, int16 *w) noexcept {
+    inline void NVERTEX(int32 *v, int32 *w) noexcept {
       if (T_ptr != T.data() + Tlen) {
         READT(v, w);
         v_j = *v;
@@ -138,7 +197,7 @@ namespace TreeTools {
       }
     }
 
-    inline void NVERTEX_short(int16 *v, int16 *w) noexcept {
+    inline void NVERTEX_short(int32 *v, int32 *w) noexcept {
       // Don't count all-tips or all-ingroup: vertices 0, ROOT, Ingp.
       if (T_ptr != T.data() + Tlen_short) {
         READT(v, w);
@@ -149,69 +208,69 @@ namespace TreeTools {
       }
     }
 
-    inline int16 LEFTLEAF() noexcept {
+    inline int32 LEFTLEAF() noexcept {
       // If NVERTEX has returned entry <vj, wj> in T, the leftmost leaf in the
       // subtree rooted at vj has entry <vk, wk> where k = j - wj.
       // This function procedure returns Vk as its value.
       return leftmost_leaf[v_j - 1];
     }
 
-    inline void SET_LEFTMOST(int16 index, int16 val) noexcept {
+    inline void SET_LEFTMOST(int32 index, int32 val) noexcept {
       leftmost_leaf[index - 1] = val;
     }
 
-    [[nodiscard]] inline int16 GET_LEFTMOST(int16 index) noexcept {
+    [[nodiscard]] inline int32 GET_LEFTMOST(int32 index) noexcept {
       return leftmost_leaf[index - 1];
     }
 
     // Procedures to manipulate cluster tables, per Table 4 of Day 1985.
 
-    inline int16 ENCODE(const int v) noexcept {
+    inline int32 ENCODE(const int v) noexcept {
       // This function procedure returns as its value the internal label
       // assigned to leaf v
       // MS note: input = v; output = X[v, 3]
       return internal_label_ptr[v];
     }
 
-    inline int16 DECODE(const int16 internal_relabeling) noexcept {
+    inline int32 DECODE(const int32 internal_relabeling) noexcept {
       // MS: input = X[v, 3], output = v
       return visited_nth[internal_relabeling - 1];
     }
 
-    inline void VISIT_LEAF (const int16* leaf, int16* n_visited) {
+    inline void VISIT_LEAF (const int32* leaf, int32* n_visited) {
       visited_nth[(*n_visited)++] = *leaf;
       internal_label[*leaf] = *n_visited;
     }
 
     Rcpp::IntegerVector X_decode() {
       Rcpp::IntegerVector ret(N());
-      for (int16 i = n_leaves; i--; ) {
+      for (int32 i = n_leaves; i--; ) {
         ret(i) = DECODE(i + 1);
       }
       return ret;
     }
     
-    [[nodiscard]] inline int16 X_left(int16 row) noexcept {
+    [[nodiscard]] inline int32 X_left(int32 row) noexcept {
       ASSERT(row > 0);
       ASSERT(row <= X_ROWS);
-      ASSERT(x_rows[row - 1].L < std::numeric_limits<int16>::max());
+      ASSERT(x_rows[row - 1].L < std::numeric_limits<int32>::max());
       return x_rows[row - 1].L;
     }
     
-    [[nodiscard]] inline int16 X_right(int16 row) noexcept {
+    [[nodiscard]] inline int32 X_right(int32 row) noexcept {
       ASSERT(row > 0);
       ASSERT(row <= X_ROWS);
-      ASSERT(x_rows[row - 1].R < std::numeric_limits<int16>::max());
+      ASSERT(x_rows[row - 1].R < std::numeric_limits<int32>::max());
       return x_rows[row - 1].R;
     }
     
-    inline void setX_left(int16 row, int16 value) noexcept {
+    inline void setX_left(int32 row, int32 value) noexcept {
       ASSERT(row > 0);
       ASSERT(row <= X_ROWS);
       x_rows[row - 1].L = value;
     }
     
-    inline void setX_right(int16 row, int16 value) noexcept {
+    inline void setX_right(int32 row, int32 value) noexcept {
       ASSERT(row > 0);
       ASSERT(row <= X_ROWS);
       x_rows[row - 1].R = value;
@@ -219,7 +278,7 @@ namespace TreeTools {
     
     Rcpp::IntegerMatrix X_contents() noexcept {
       Rcpp::IntegerMatrix ret(X_ROWS, 2);
-      for (int16 i = X_ROWS; i--; ) {
+      for (int32 i = X_ROWS; i--; ) {
         ret(i, 0) = x_rows[i].L;
         ret(i, 1) = x_rows[i].R;
       }
@@ -228,38 +287,38 @@ namespace TreeTools {
 
     
     // Required by TreeDist 2.9.2
-    // TODO Remove in later version, to prefer CLUSTONL(int16 L, R)
-    [[nodiscard]] inline bool CLUSTONL(int16* L, int16* R) noexcept {
+    // TODO Remove in later version, to prefer CLUSTONL(int32 L, R)
+    [[nodiscard]] inline bool CLUSTONL(int32* L, int32* R) noexcept {
       return X_left(*L) == *L && X_right(*L) == *R;
     }
     
     // Required by TreeDist 2.9.2
-    // TODO Remove in later version, to prefer CLUSTONR(int16 L, R)
-    [[nodiscard]] inline bool CLUSTONR(int16* L, int16* R) noexcept {
+    // TODO Remove in later version, to prefer CLUSTONR(int32 L, R)
+    [[nodiscard]] inline bool CLUSTONR(int32* L, int32* R) noexcept {
       return X_left(*R) == *L && X_right(*R) == *R;
     }
     
     // Required by TreeDist 2.9.2
-    // TODO Remove in later version, to prefer ISCLUST(int16 L, R)
-    [[nodiscard]] inline bool ISCLUST(int16* L, int16* R) noexcept {
+    // TODO Remove in later version, to prefer ISCLUST(int32 L, R)
+    [[nodiscard]] inline bool ISCLUST(int32* L, int32* R) noexcept {
       // This function procedure returns value true if cluster <L,R> is in X;
       // otherwise it returns value false
       return CLUSTONL(*L, *R) || CLUSTONR(*L, *R);
     }
 
-    [[nodiscard]] inline bool CLUSTONL(int16 L, int16 R) noexcept {
+    [[nodiscard]] inline bool CLUSTONL(int32 L, int32 R) noexcept {
       ASSERT(L > 0 && L <= X_ROWS);
       const ClusterRow &r = x_rows[L - 1];
       return (r.L == L) & (r.R == R);
     }
     
-    [[nodiscard]] inline bool CLUSTONR(int16 L, int16 R) noexcept {
+    [[nodiscard]] inline bool CLUSTONR(int32 L, int32 R) noexcept {
       ASSERT(R > 0 && R <= X_ROWS);
       const ClusterRow &r = x_rows[R - 1];
       return (r.L == L) & (r.R == R);
     }
     
-    [[nodiscard]] inline bool ISCLUST(int16 L, int16 R) noexcept {
+    [[nodiscard]] inline bool ISCLUST(int32 L, int32 R) noexcept {
       // This function procedure returns value true if cluster <L,R> is in X;
       // otherwise it returns value false
       ASSERT(L > 0 && L <= X_ROWS);
@@ -281,25 +340,28 @@ namespace TreeTools {
     }
     
     // Required by TreeDist 2.9.2
-    // TODO Remove in later version, to prefer SETSWX(int16 row)
-    inline void SETSWX(int16* row) noexcept {
+    // TODO Remove in later version, to prefer SETSWX(int32 row)
+    inline void SETSWX(int32* row) noexcept {
       // Only increment our counter on a 0 -> 1 transition
       const auto idx = static_cast<std::size_t>(*row);
+      ASSERT(idx > 0 && idx <= static_cast<size_t>(X_ROWS));
       if (!Xswitch[idx]) {
-        Xswitch[idx] = true;
+        Xswitch.set(idx, true);
         ++xswitch_set_count;
       }
     }
     
     inline void SETSWX(std::size_t row) noexcept {
       // Only increment our counter on a 0 -> 1 transition
+      ASSERT(row > 0 && row <= static_cast<size_t>(X_ROWS));
       if (!Xswitch[row]) {
-        Xswitch[row] = true;
+        Xswitch.set(row, true);
         ++xswitch_set_count;
       }
     }
 
-    [[nodiscard]] inline bool GETSWX(int16* row) noexcept {
+    [[nodiscard]] inline bool GETSWX(int32* row) noexcept {
+      ASSERT(*row > 0 && *row <= X_ROWS);
       return Xswitch[*row];
     }
 
@@ -308,10 +370,10 @@ namespace TreeTools {
     }
     
     // Required by TreeDist 2.9.2
-    // TODO Remove in later version, to prefer SETSW(int16 L, R)
-    inline void SETSW(int16 *L, int16 *R) noexcept {
-      const int16 l = *L;
-      const int16 r = *R;
+    // TODO Remove in later version, to prefer SETSW(int32 L, R)
+    inline void SETSW(int32 *L, int32 *R) noexcept {
+      const int32 l = *L;
+      const int32 r = *R;
       // If <L,R> is a cluster in X, 
       // this procedure sets the cluster switch for <L,R>.
       if (CLUSTONL(l, r)) {
@@ -323,7 +385,7 @@ namespace TreeTools {
       }
     }
     
-    inline void SETSW(int16 L, int16 R) noexcept {
+    inline void SETSW(int32 L, int32 R) noexcept {
       // If <L,R> is a cluster in X, 
       // this procedure sets the cluster switch for <L,R>.
       if (CLUSTONL(L, R)) {
@@ -339,7 +401,7 @@ namespace TreeTools {
       // This procedure inspects every cluster switch in X.
       // If the switch for cluster <L,R> is cleared, UPDATE deletes <L,R>
       // from X; thereafter ISCLUST(X,L,R) will return the value false.
-      for (int16 i = X_ROWS; i--; ) {
+      for (int32 i = X_ROWS; i--; ) {
         if (!(Xswitch[i])) {
           x_rows[i].L = 0;
           x_rows[i].R = 0;
@@ -347,7 +409,7 @@ namespace TreeTools {
       }
     }
 
-    [[nodiscard]] inline int16 SHARED() noexcept {
+    [[nodiscard]] inline int32 SHARED() noexcept {
       // Used by COMCLUST in TreeDist::Day_1985.cpp
       return n_shared;
     }
@@ -361,14 +423,14 @@ namespace TreeTools {
       enumeration = 0;
     }
 
-    inline void NCLUS(int16* L, int16* R) noexcept {
+    inline void NCLUS(int32* L, int32* R) noexcept {
       // This procedure returns the next cluster <L,R> in the current
       // enumeration of clusters in X.
       // If m clusters are in X, they are returned by the first m invocations
       // of NCLUS after initialization by XRESET; thereafter NCLUS returns the
       // invalid cluster <0,0>.
       ASSERT(enumeration < X_ROWS);
-      const int16 row = static_cast<int16>(enumeration + 1); // rows are 1-based
+      const int32 row = static_cast<int32>(enumeration + 1); // rows are 1-based
       *L = X_left(row);
       *R = X_right(row);
       ++enumeration;
@@ -384,21 +446,26 @@ namespace TreeTools {
     // BEGIN
     n_internal = rooted["Nnode"]; // = M
     Rcpp::CharacterVector leaf_labels = rooted["tip.label"];
-    if (leaf_labels.length() > int(ct_max_leaves)) {
-      Rcpp::stop("Tree has too many leaves. "
-                 "Contact the 'TreeTools' maintainer.");
-    }
-    ASSERT(ct_max_leaves <= std::numeric_limits<int16>::max());
-    n_leaves = int16(leaf_labels.length()); // = N
-    if (double(edge.nrow()) > double(std::numeric_limits<int16>::max())) {
-      Rcpp::stop("Tree has too many edges. "
-                 "Contact the 'TreeTools' maintainer.");
-    }
-    n_edge = int16(edge.nrow());
-    const int16 n_vertex = M() + N();
+    if (leaf_labels.length() > int(ct_max_leaves_heap)) {
+      std::ostringstream msg;
+      msg << "Tree has too many leaves (>" << ct_max_leaves_heap << "). "
+          << "Contact the 'TreeTools' maintainer.";
+      Rcpp::stop(msg.str());
+    }
+    ASSERT(ct_max_leaves_heap <= std::numeric_limits<int32>::max());
+    n_leaves = int32(leaf_labels.length()); // = N
+    if (double(edge.nrow()) > double(std::numeric_limits<int32>::max())) {
+      std::ostringstream msg;
+      msg << "Tree has too many edges (" << edge.nrow() << " > " << 
+        std::numeric_limits<int32>::max() << "). " << 
+          "Contact the 'TreeTools' maintainer.";
+      Rcpp::stop(msg.str());
+    }
+    n_edge = int32(edge.nrow());
+    const int32 n_vertex = M() + N();
     Tlen = 2 * n_vertex;
     Tlen_short = Tlen - (2 * 3);
-    T = std::vector<int16>(Tlen);
+    T = std::vector<int32>(Tlen);
     T_ptr = T.data();
 
     resize_uninitialized(leftmost_leaf, n_vertex);
@@ -408,21 +475,21 @@ namespace TreeTools {
     ASSERT(visited_nth.size()   >= static_cast<size_t>(n_leaves));
     ASSERT(internal_label.size()>= static_cast<size_t>(1 + n_leaves));
     internal_label_ptr = internal_label.data();
-    int16 n_visited = 0;
-    std::vector<int16> weights(1 + n_vertex);
+    int32 n_visited = 0;
+    std::vector<int32> weights(1 + n_vertex);
 
-    for (int16 i = 1; i != n_leaves + 1; ++i) {
+    for (int32 i = 1; i != n_leaves + 1; ++i) {
       SET_LEFTMOST(i, i); // Line 402
       weights[i] = 0;
     }
-    for (int16 i = 1 + n_leaves; i != 1 + n_vertex; ++i) {
+    for (int32 i = 1 + n_leaves; i != 1 + n_vertex; ++i) {
       SET_LEFTMOST(i, 0);
       weights[i] = 0;
     }
-    for (int16 i = n_edge; i--; ) {
-      const int16
-        parent_i = int16(edge(i, 0)),
-        child_i = int16(edge(i, 1))
+    for (int32 i = n_edge; i--; ) {
+      const int32
+        parent_i = int32(edge(i, 0)),
+        child_i = int32(edge(i, 1))
       ;
       if (!GET_LEFTMOST(parent_i)) {
         SET_LEFTMOST(parent_i, GET_LEFTMOST(child_i));
@@ -436,11 +503,15 @@ namespace TreeTools {
         ENTER(child_i, weights[child_i]);
       }
     }
-    ENTER(int16(edge(0, 0)), weights[edge(0, 0)]);
+    ENTER(int32(edge(0, 0)), weights[edge(0, 0)]);
 
     // BUILD Cluster table
     X_ROWS = n_leaves;
     x_rows = std::vector<ClusterRow>(X_ROWS);
+    
+    // Initialize Xswitch with appropriate size
+    // Uses stack allocation for small trees, heap for large ones
+    Xswitch.resize(n_leaves + 1);
 
     // This procedure constructs in X descriptions of the clusters in a
     // rooted tree described by the postorder sequence T with weights,
@@ -449,11 +520,11 @@ namespace TreeTools {
     // simply described by a pair <L,R> of internal labels.
 
     TRESET();
-    for (int16 i = 1; i != N(); ++i) {
+    for (int32 i = 1; i != N(); ++i) {
       setX_left(i, 0);
       setX_right(i, 0);
     }
-    int16 leafcode = 0, v, w, L, R = UNINIT, loc;
+    int32 leafcode = 0, v, w, L, R = UNINIT, loc;
 
     NVERTEX(&v, &w);
     while (v) {
diff --git a/src/consensus.cpp b/src/consensus.cpp
index 02b02c47c..03cfec261 100644
--- a/src/consensus.cpp
+++ b/src/consensus.cpp
@@ -9,19 +9,21 @@ using namespace Rcpp;
 #include <vector> /* for vector */
 
 using TreeTools::ct_stack_size;
-using TreeTools::ct_max_leaves;
+using TreeTools::ct_stack_threshold;
+using TreeTools::ct_max_leaves_heap;
 
-// trees is a list of objects of class phylo, all with the same tip labels
-// (try RenumberTips(trees, trees[[1]]))
-// Per #168, unexpected behaviour if root position differs in non-preorder trees
-// Further investigation could be beneficial; for now, suggest applying
-// the function to preorder trees only.
-// [[Rcpp::export]]
-RawMatrix consensus_tree(const List trees, const NumericVector p) {
-  int16 v = 0;
-  int16 w = 0;
-  int16 L, R, N, W;
-  int16 L_j, R_j, N_j, W_j;
+// Helper template function to perform consensus computation
+// Uses StackContainer for the S array (either std::array or std::vector)
+template<typename StackContainer>
+RawMatrix consensus_tree_impl(
+  const List& trees,
+  const NumericVector& p,
+  StackContainer& S
+) {
+  int32 v = 0;
+  int32 w = 0;
+  int32 L, R, N, W;
+  int32 L_j, R_j, N_j, W_j;
   
   const int32 n_trees = trees.length();
   const int32 frac_thresh = int32(n_trees * p[0]) + 1;
@@ -37,7 +39,6 @@ RawMatrix consensus_tree(const List trees, const NumericVector p) {
   const int32 ntip_3  = n_tip - 3;
   const int32 nbin    = (n_tip + 7) / 8;  // bytes per row in packed output
   
-  std::array<int32, ct_stack_size * ct_max_leaves> S;
   std::vector<int32> split_count(n_tip, 1);
   
   // Packed output: each row has nbin bytes
@@ -54,25 +55,29 @@ RawMatrix consensus_tree(const List trees, const NumericVector p) {
     std::fill(split_count.begin(), split_count.end(), 1);
     
     for (int32 j = i + 1; j < n_trees; ++j) {
+      ASSERT(tables[i].N() == tables[j].N());
       
       tables[i].CLEAR();
       
       tables[j].TRESET();
       tables[j].READT(&v, &w);
       
-      int16 j_pos = 0;
+      int32 j_pos = 0;
       int32 Spos = 0; // Empty the stack S. Used in CT_PUSH / CT_POP macros.
       
       do {
         if (CT_IS_LEAF(v)) {
+          CT_ASSERT_CAN_PUSH();
           CT_PUSH(tables[i].ENCODE(v), tables[i].ENCODE(v), 1, 1);
         } else {
+          CT_ASSERT_CAN_POP();
           CT_POP(L, R, N, W_j);
           
           W = 1 + W_j;
           w = w - W_j;
           
           while (w) {
+            CT_ASSERT_CAN_POP();
             CT_POP(L_j, R_j, N_j, W_j);
             if (L_j < L) L = L_j;
             if (R_j > R) R = R_j;
@@ -81,6 +86,7 @@ RawMatrix consensus_tree(const List trees, const NumericVector p) {
             w = w - W_j;
           }
           
+          CT_ASSERT_CAN_PUSH();
           CT_PUSH(L, R, N, W);
           
           ++j_pos;
@@ -136,3 +142,33 @@ RawMatrix consensus_tree(const List trees, const NumericVector p) {
     return ret;
   }
 }
+
+// trees is a list of objects of class phylo, all with the same tip labels
+// (try RenumberTips(trees, trees[[1]]))
+// Per #168, unexpected behaviour if root position differs in non-preorder trees
+// Further investigation could be beneficial; for now, suggest applying
+// the function to preorder trees only.
+// [[Rcpp::export]]
+RawMatrix consensus_tree(const List trees, const NumericVector p) {
+  // First, peek at the tree size to determine allocation strategy
+  // We'll create a temporary ClusterTable just to check the size
+  try {
+    TreeTools::ClusterTable temp_table(Rcpp::List(trees(0)));
+    const int32 n_tip = temp_table.N();
+    
+    if (n_tip <= ct_stack_threshold) {
+      // Small tree: use stack-allocated array
+      std::array<int32, ct_stack_size * ct_stack_threshold> S;
+      return consensus_tree_impl(trees, p, S);
+    } else {
+      // Large tree: use heap-allocated vector
+      std::vector<int32> S(ct_stack_size * n_tip);
+      return consensus_tree_impl(trees, p, S);
+    }
+  } catch(const std::exception& e) {
+    Rcpp::stop(e.what());
+  }
+  
+  ASSERT(false && "Unreachable code in consensus_tree");
+  return RawMatrix(0, 0);
+}
diff --git a/src/splits_to_tree.cpp b/src/splits_to_tree.cpp
index cbd3e5aba..11ed9eab4 100644
--- a/src/splits_to_tree.cpp
+++ b/src/splits_to_tree.cpp
@@ -6,9 +6,9 @@ using namespace Rcpp;
 #include "../inst/include/TreeTools/renumber_tree.h"
 using namespace TreeTools;
 
-inline void insert_ancestor(const int16 tip, const int16 *next_node,
-                            int16* parent,
-                            int16* patriarch) {
+inline void insert_ancestor(const int32 tip, const int32 *next_node,
+                            int32* parent,
+                            int32* patriarch) {
   if (patriarch[tip]) {
     parent[patriarch[tip]] = *next_node;
   } else {
@@ -19,7 +19,7 @@ inline void insert_ancestor(const int16 tip, const int16 *next_node,
 
 // [[Rcpp::export]]
 IntegerMatrix splits_to_edge(const RawMatrix splits, const IntegerVector nTip) {
-  const int16 n_tip = int16(nTip[0]);
+  const int32 n_tip = int32(nTip[0]);
   if (splits.nrow() == 0) {
     IntegerMatrix ret(n_tip, 2);
     for (int i = 0; i < n_tip; ++i) {
@@ -34,16 +34,16 @@ IntegerMatrix splits_to_edge(const RawMatrix splits, const IntegerVector nTip) {
   const bool use_heap = (n_tip > SL_MAX_TIPS) || (x.n_splits > SL_MAX_SPLITS);
   
   // Stack allocation for small trees (fast path)
-  alignas(64) std::array<int16, SL_MAX_TIPS + SL_MAX_SPLITS> stack_parent{};
-  alignas(64) std::array<int16, SL_MAX_TIPS> stack_patriarch{};
+  alignas(64) std::array<int32, SL_MAX_TIPS + SL_MAX_SPLITS> stack_parent{};
+  alignas(64) std::array<int32, SL_MAX_TIPS> stack_patriarch{};
   
   // Heap allocation for large trees
-  std::vector<int16> heap_parent;
-  std::vector<int16> heap_patriarch;
+  std::vector<int32> heap_parent;
+  std::vector<int32> heap_patriarch;
   
   // Pointers to active storage
-  int16* parent;
-  int16* patriarch;
+  int32* parent;
+  int32* patriarch;
   
   if (use_heap) {
     const size_t parent_size = static_cast<size_t>(n_tip) + 
@@ -58,39 +58,39 @@ IntegerMatrix splits_to_edge(const RawMatrix splits, const IntegerVector nTip) {
   }
   
   // Allocate split_order appropriately
-  std::vector<int16> split_order(x.n_splits);
+  std::vector<int32> split_order(x.n_splits);
   std::iota(split_order.begin(), split_order.end(), 0);
   std::sort(split_order.begin(), split_order.end(),
-            [&in_split = x.in_split](int16 a, int16 b) {
+            [&in_split = x.in_split](int32 a, int32 b) {
               return in_split[a] > in_split[b];
             });
 
-  int16 next_node = n_tip;
-  for (int16 split = x.n_splits; split--; ) { // Work back through splits
+  int32 next_node = n_tip;
+  for (int32 split = x.n_splits; split--; ) { // Work back through splits
     if (split > 0) {
       __builtin_prefetch(&x.state[split_order[split - 1]][0], 0, 3);
     }
     
-    for (int16 bin = 0; bin < x.n_bins; ++bin) {
+    for (int32 bin = 0; bin < x.n_bins; ++bin) {
       splitbit chunk = x.state[split_order[split]][bin];
       if (!chunk) continue;
-      const int16 base_tip = bin * SL_BIN_SIZE;
+      const int32 base_tip = bin * SL_BIN_SIZE;
       while (chunk) {
-        const int16 bin_tip = __builtin_ctzll(chunk); // count trailing zeros
-        const int16 tip = base_tip + bin_tip;
+        const int32 bin_tip = __builtin_ctzll(chunk); // count trailing zeros
+        const int32 tip = base_tip + bin_tip;
         insert_ancestor(tip, &next_node, parent, patriarch);
         chunk &= chunk - 1; // clear lowest set bit
       }
     }
     ++next_node;
   }
-  for (int16 tip = 0; tip < n_tip; ++tip) {
+  for (int32 tip = 0; tip < n_tip; ++tip) {
     insert_ancestor(tip, &next_node, parent, patriarch);
   }
 
-  const int16 n_edge = n_tip + x.n_splits;
+  const int32 n_edge = n_tip + x.n_splits;
   IntegerVector edge1(n_edge), edge2(n_edge);
-  for (int16 i = 0; i < n_edge; ++i) {
+  for (int32 i = 0; i < n_edge; ++i) {
     edge1[i] = parent[i] + 1;
     edge2[i] = i + 1;
   }
diff --git a/tests/testthat/test-ClusterTable.R b/tests/testthat/test-ClusterTable.R
index 8ffac908b..5b2a8644b 100644
--- a/tests/testthat/test-ClusterTable.R
+++ b/tests/testthat/test-ClusterTable.R
@@ -1,9 +1,13 @@
 test_that("ClusterTable fails gracefully", {
-  bigTree <- PectinateTree(2^14 + 1)
+  bigTree <- PectinateTree(100001)
   expect_error(
     as.ClusterTable(bigTree),
-    "Tree has too many leaves. Contact the .TreeTools. maintainer."
+    "Tree has too many leaves.*100000"
   )
+  
+  mediumTree <- PectinateTree(20000)
+  ct <- as.ClusterTable(mediumTree)
+  expect_equal(attr(ct, "nTip"), 20000)
 })
 
 test_that("ClusterTable class behaves", {
diff --git a/tests/testthat/test-consensus.R b/tests/testthat/test-consensus.R
index 074405987..0e42dd62c 100644
--- a/tests/testthat/test-consensus.R
+++ b/tests/testthat/test-consensus.R
@@ -17,10 +17,16 @@ test_that("Consensus() errors", {
   expect_equal(Consensus(c(halfTree)), halfTree)
   expect_equal(Consensus(list(halfTree)), halfTree)
   
+  # Test that trees larger than heap limit fail gracefully
   expect_error(
-    Consensus(c(BalancedTree(33333), PectinateTree(33333))),
-    "too many leaves"
+    Consensus(c(BalancedTree(100001), PectinateTree(100001))),
+    "too many leaves.*100000"
   )
+  
+  skip_on_cran() # Slow!
+  largeTree <- BalancedTree(33333)
+  consensus_large <- Consensus(c(largeTree, largeTree))
+  expect_equal(NTip(consensus_large), 33333)
 })
 
 test_that("Consensus()", {