Replace b_complement stack arrays with std::vector for large-tree support

The two splitbit b_complement[SL_MAX_SPLITS][SL_MAX_BINS] stack arrays in
robinson_foulds_distance() and robinson_foulds_info() would overflow the
stack when compiled against TreeTools with SL_MAX_TIPS = 32768 (128 GB).

Replace with std::vector<splitbit> sized to actual dimensions (b.n_splits *
n_bins).  These are serial per-pair paths (reportMatching = TRUE), so heap
allocation cost is negligible.

Also upgrade assert() to static_assert() in tree_distances.h for the int16
width checks — these now fire at compile time rather than silently passing
in release builds.

Author: ms609

src/tree_distances.cpp

@@ -51,29 +51,32 @@ inline List robinson_foulds_distance(const RawMatrix &x, const RawMatrix &y,
 
   grf_match matching(a.n_splits, NA_INTEGER);
 
-  splitbit b_complement[SL_MAX_SPLITS][SL_MAX_BINS];
+  const int32 n_bins = a.n_bins;
+  std::vector<splitbit> b_complement(b.n_splits * n_bins);
   for (int32 i = b.n_splits; i--; ) {
+    splitbit* bc_i = &b_complement[i * n_bins];
     for (int32 bin = last_bin; bin--; ) {
-      b_complement[i][bin] = ~b.state[i][bin];
+      bc_i[bin] = ~b.state[i][bin];
     }
-    b_complement[i][last_bin] = b.state[i][last_bin] ^ unset_mask;
+    bc_i[last_bin] = b.state[i][last_bin] ^ unset_mask;
   }
 
   for (int32 ai = a.n_splits; ai--; ) {
     for (int32 bi = b.n_splits; bi--; ) {
 
       bool all_match = true;
       bool all_complement = true;
+      const splitbit* bc_bi = &b_complement[bi * n_bins];
 
-      for (int32 bin = 0; bin < a.n_bins; ++bin) {
+      for (int32 bin = 0; bin < n_bins; ++bin) {
         if ((a.state[ai][bin] != b.state[bi][bin])) {
           all_match = false;
           break;
         }
       }
       if (!all_match) {
-        for (int32 bin = 0; bin < a.n_bins; ++bin) {
-          if (a.state[ai][bin] != b_complement[bi][bin]) {
+        for (int32 bin = 0; bin < n_bins; ++bin) {
+          if (a.state[ai][bin] != bc_bi[bin]) {
             all_complement = false;
             break;
           }
@@ -105,29 +108,31 @@ inline List robinson_foulds_info(const RawMatrix &x, const RawMatrix &y,
 
   grf_match matching(a.n_splits, NA_INTEGER);
 
-  /* Dynamic allocation 20% faster for 105 tips, but VLA not permitted in C11 */
-  splitbit b_complement[SL_MAX_SPLITS][SL_MAX_BINS]; 
+  const int16 n_bins = a.n_bins;
+  std::vector<splitbit> b_complement(b.n_splits * n_bins);
   for (int16 i = 0; i < b.n_splits; i++) {
+    splitbit* bc_i = &b_complement[i * n_bins];
     for (int16 bin = 0; bin < last_bin; ++bin) {
-      b_complement[i][bin] = ~b.state[i][bin];
+      bc_i[bin] = ~b.state[i][bin];
     }
-    b_complement[i][last_bin] = b.state[i][last_bin] ^ unset_mask;
+    bc_i[last_bin] = b.state[i][last_bin] ^ unset_mask;
   }
 
   for (int16 ai = 0; ai < a.n_splits; ++ai) {
     for (int16 bi = 0; bi < b.n_splits; ++bi) {
 
       bool all_match = true, all_complement = true;
+      const splitbit* bc_bi = &b_complement[bi * n_bins];
 
-      for (int16 bin = 0; bin < a.n_bins; ++bin) {
+      for (int16 bin = 0; bin < n_bins; ++bin) {
         if ((a.state[ai][bin] != b.state[bi][bin])) {
           all_match = false;
           break;
         }
       }
       if (!all_match) {
-        for (int16 bin = 0; bin < a.n_bins; ++bin) {
-          if ((a.state[ai][bin] != b_complement[bi][bin])) {
+        for (int16 bin = 0; bin < n_bins; ++bin) {
+          if ((a.state[ai][bin] != bc_bi[bin])) {
             all_complement = false;
             break;
           }

src/tree_distances.h

@@ -39,7 +39,8 @@ namespace TreeDist {
 
   // Returns lg2_unrooted[x] - lg2_trees_matching_split(y, x - y)
   [[nodiscard]] inline double mmsi_pair_score(const int16 x, const int16 y) noexcept {
-    assert(SL_MAX_TIPS + 2 <= std::numeric_limits<int16>::max()); // verify int16 ok
+    static_assert(SL_MAX_TIPS + 2 <= std::numeric_limits<int16>::max(),
+                  "int16 too narrow for SL_MAX_TIPS");
 
     return lg2_unrooted[x] - (lg2_rooted[y] + lg2_rooted[x - y]);
   }
@@ -60,7 +61,8 @@ namespace TreeDist {
 
 
 [[nodiscard]] inline double one_overlap(const int16 a, const int16 b, const int16 n) noexcept {
-    assert(SL_MAX_TIPS + 2 <= std::numeric_limits<int16>::max()); // verify int16 ok
+    static_assert(SL_MAX_TIPS + 2 <= std::numeric_limits<int16>::max(),
+                  "int16 too narrow for SL_MAX_TIPS");
     if (a == b) {
       return lg2_rooted[a] + lg2_rooted[n - a];
     }
@@ -71,7 +73,8 @@ namespace TreeDist {
   }
 
   [[nodiscard]] inline double one_overlap_notb(const int16 a, const int16 n_minus_b, const int16 n) noexcept {
-    assert(SL_MAX_TIPS + 2 <= std::numeric_limits<int16>::max()); // verify int16 ok
+    static_assert(SL_MAX_TIPS + 2 <= std::numeric_limits<int16>::max(),
+                  "int16 too narrow for SL_MAX_TIPS");
     const int16 b = n - n_minus_b;
     if (a == b) {
       return lg2_rooted[b] + lg2_rooted[n_minus_b];
@@ -90,7 +93,8 @@ namespace TreeDist {
                        const int16 n_tips, const int16 in_a,
                        const int16 in_b, const int16 n_bins) noexcept {
 
-    assert(SL_MAX_BINS <= INT16_MAX);
+    static_assert(SL_MAX_BINS <= INT16_MAX,
+                  "int16 too narrow for SL_MAX_BINS");
 
     int16 n_ab = 0;
     for (int16 bin = 0; bin < n_bins; ++bin) {