diff --git a/.github/workflows/R-CMD-check.yml b/.github/workflows/R-CMD-check.yml
index bbdc2526..fd5fde60 100644
--- a/.github/workflows/R-CMD-check.yml
+++ b/.github/workflows/R-CMD-check.yml
@@ -338,7 +338,7 @@ jobs:
       !cancelled() && needs.sense-check.result == 'success' &&
       contains(fromJSON('["success", "skipped"]'), needs.codemeta.result) &&
       (github.event_name == 'pull_request' || github.event_name == 'workflow_dispatch')
-    runs-on: ubuntu-latest
+    runs-on: ubuntu-26.04 # RSPM below must name the same release
     # A hung setup or a wedged benchmark otherwise runs to the 6-hour default,
     # holding a runner the whole time.
     timeout-minutes: 45
@@ -354,7 +354,7 @@ jobs:
 
     env:
       GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
-      RSPM: "https://packagemanager.posit.co/cran/__linux__/noble/latest"
+      RSPM: "https://packagemanager.posit.co/cran/__linux__/resolute/latest"
       _R_CHECK_BUILD_VIGNETTES_: false
       _R_CHECK_CRAN_INCOMING_: false
       _R_CHECK_FORCE_SUGGESTS_: false
@@ -672,7 +672,7 @@ jobs:
     # and codemeta is skipped outside pull requests, so without this the
     # backstop never ran on push, schedule or workflow_dispatch.
     if: ${{ !cancelled() && needs.core.result == 'success' }}
-    runs-on: ubuntu-latest
+    runs-on: ubuntu-26.04
     timeout-minutes: 120  # revdepchecks are genuinely long
 
     name: revdepcheck (${{ matrix.config.deps }})
diff --git a/NEWS.md b/NEWS.md
index df06e482..718c7c07 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -35,6 +35,10 @@
 
 ## Performance
 
+- `CompareAll()` no longer takes time quadratic in the number of pairs on
+  older versions of R (e.g. 4.1), where each element access to a
+  `multiPhylo` object copied the whole list.
+
 - Trees with different leaf sets are pruned to their shared leaves in C++,
   without constructing pruned trees in R: `TreeDistance()` on 150 gene trees
   is c. 15× faster, with lower memory use (#123).
diff --git a/R/tree_distance_utilities.R b/R/tree_distance_utilities.R
index 518b7bba..4cb4976d 100644
--- a/R/tree_distance_utilities.R
+++ b/R/tree_distance_utilities.R
@@ -535,28 +535,31 @@ Ntropy <- function(...) {
 CompareAll <- function(x, Func, FUN.VALUE = Func(x[[1]], x[[1]], ...),
                         ...) {
   nTree <- length(x)
+  labels <- names(x)
   countUp <- seq_len(nTree - 1)
-  i <- x[rep(countUp, rev(countUp))]
-  j <- x[unlist(sapply(countUp, function(n) n + seq_len(nTree - n)))]
-  
+  i <- rep(countUp, rev(countUp))
+  j <- unlist(lapply(countUp, function(n) n + seq_len(nTree - n)))
+  # A plain list: on R 4.1, each `[[` on a multiPhylo copies the whole list
+  x <- lapply(seq_len(nTree), function(k) x[[k]])
+
   cluster <- getOption("TreeDist-cluster")
   ret <- if (is.null(cluster)) {
     cli_progress_bar("Comparing", total = length(i))
     vapply(seq_along(i), function(k) {
       cli_progress_update(1, .envir = parent.frame(2))
-      Func(i[[k]], j[[k]], ...)
+      Func(x[[i[k]]], x[[j[k]]], ...)
       }, FUN.VALUE)
   } else {
-    do.call("cbind", 
-            parLapply(cluster, seq_along(i), 
-                      function(k, Func) Func(i[[k]], j[[k]]),
+    do.call("cbind",
+            parLapply(cluster, seq_along(i),
+                      function(k, Func) Func(x[[i[k]]], x[[j[k]]]),
                       Func = Func))
   }
-  
+
   .WrapReturn <- function(dists) {
     structure(dists,
               Size = nTree,
-              Labels = names(x),
+              Labels = labels,
               Diag = FALSE,
               Upper = TRUE,
               class = "dist")
diff --git a/tests/testthat/_problems/test-parallel-19.R b/tests/testthat/_problems/test-parallel-19.R
new file mode 100644
index 00000000..37e11cd1
--- /dev/null
+++ b/tests/testthat/_problems/test-parallel-19.R
@@ -0,0 +1,29 @@
+# Extracted from test-parallel.R:19
+
+# setup ------------------------------------------------------------------------
+library(testthat)
+test_env <- simulate_test_env(package = "TreeDist", path = "..")
+attach(test_env, warn.conflicts = FALSE)
+
+# test -------------------------------------------------------------------------
+library("TreeTools")
+trees <- as.phylo(0:20, 20)
+suppressMessages({
+  expect_equal(options("TreeDist-cluster"), StartParallel(2))
+  cl <- getOption("TreeDist-cluster")
+  
+  parallel <- ClusteringInfoDistance(trees)
+  SetParallel(NULL)
+  expect_null(GetParallel())
+  expect_equal(ClusteringInfoDistance(trees), parallel)
+
+  expect_equal(options("TreeDist-cluster"), SetParallel(cl))
+  expect_equal(cl, getOption("TreeDist-cluster"))
+  expect_equal(cl, GetParallel())
+  parallel <- CompareAll(trees, NNIDist)
+  SetParallel(NULL)
+  expect_equal(CompareAll(trees, NNIDist), parallel)
+  expect_false(StopParallel())
+  SetParallel(cl)
+  expect_true(StopParallel())
+    })
diff --git a/tests/testthat/test-tree_distance_nni.R b/tests/testthat/test-tree_distance_nni.R
index 09ec2f7d..b33c69d2 100644
--- a/tests/testthat/test-tree_distance_nni.R
+++ b/tests/testthat/test-tree_distance_nni.R
@@ -186,13 +186,26 @@ test_that("NNIDist() lower bounds are valid", {
   testthat::expect_equal(max(exact), NNIDiameter(nTip)[, "exact"],
                          ignore_attr = TRUE)
 
-  bounds <- lapply(NNIDist(trees), as.matrix)
+  # Bounds for each pair, in the order of exact[lower.tri(exact)].
+  # Calling the C++ kernel directly avoids NNIDist()'s per-pair R overhead,
+  # which makes 7 leaves (446,040 pairs) take minutes, and an hour on R 4.1;
+  # tip labels are already in a common order.
+  pairs <- which(lower.tri(exact), arr.ind = TRUE)
+  edges <- lapply(trees, function(tr) Postorder(tr[["edge"]]))
+  bounds <- vapply(seq_len(nrow(pairs)), function(k) {
+    cpp_nni_distance(edges[[pairs[k, 1]]], edges[[pairs[k, 2]]], nTip)
+  }, integer(6))
+  if (nTip == 6) {
+    viaR <- vapply(NNIDist(trees), function(d) as.matrix(d)[pairs],
+                   double(nrow(pairs)))
+    testthat::expect_equal(t(bounds), viaR, ignore_attr = TRUE)
+    testthat::expect_equal(colnames(viaR), rownames(bounds))
+  }
   exact <- exact[lower.tri(exact)]
-  .Bound <- function(field) bounds[[field]][lower.tri(bounds[[field]])]
-  testthat::expect_true(all(.Bound("lower") <= exact))
+  testthat::expect_true(all(bounds["lower", ] <= exact))
   for (upper in c("tight_upper", "best_upper", "loose_upper", "fack_upper",
                   "li_upper")) {
-    testthat::expect_true(all(.Bound(upper) >= exact))
+    testthat::expect_true(all(bounds[upper, ] >= exact))
   }
 }
 
@@ -201,7 +214,7 @@ test_that("NNIDist() bounds bracket exact distances: 6 leaves", {
 })
 
 test_that("NNIDist() bounds bracket exact distances: 7 leaves", {
-  skip_on_cran() # c. 70 s
+  skip_on_cran() # c. 7 s
   .ExpectBoundsBracketExact(7)
 })
 
diff --git a/tests/testthat/testthat-problems.rds b/tests/testthat/testthat-problems.rds
new file mode 100644
index 00000000..47bbda63
Binary files /dev/null and b/tests/testthat/testthat-problems.rds differ