diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md
index c5bf88fe..026525a3 100644
--- a/docs/CHANGELOG.md
+++ b/docs/CHANGELOG.md
@@ -18,6 +18,7 @@ This is an UNSTABLE release. It includes many experimental features, as well as
 * HDF5 and image stream: `mirror_y` records whether the assembled image is mirrored in Y relative to the detector's raw readout.
 * rugnux: an image integrated in pyFAI through the `.poni` file written by `--mode calibration` now comes out with the correct azimuth. Radial integration is unchanged.
 * rugnux: the `.poni` file declares pyFAI's `orientation`, which needs pyFAI 2024.01 or newer.
+* rugnux: scaling and merging are faster, with identical output.
 * The per-image resolution estimate now predicts the resolution the merged data reach, rather than reporting the highest-resolution spot found; rugnux reports the run's value as `SPOT_RESOLUTION_ESTIMATE` in its report.
 * rugnux: fixing the space group with `-S` no longer prevents the lattice from being found; the group is applied to scaling and merging rather than to the indexing search.
 * rugnux: the detector geometry is also logged in XDS's convention (`ORGX`/`ORGY`, detector axis vectors, rotation axis), so it can be compared with an XDS refinement.
diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp
index 53c7f5a8..97df6bbb 100644
--- a/image_analysis/scale_merge/RotationScaleMerge.cpp
+++ b/image_analysis/scale_merge/RotationScaleMerge.cpp
@@ -927,24 +927,30 @@ int RotationScaleMerge::ComputeAsuGroups(const HKLKeyGenerator &keygen) {
             eligible[r] = 1;
         }
     });
-    std::vector<int32_t> idx;
-    idx.reserve(n_run);
+    // The key travels with the run it belongs to, so the sort and the run-detection below read the
+    // key they are standing on instead of chasing it back through `key` - an indirect compare is a
+    // cache miss per comparison, and there are a few tens of millions of them. Which of two runs
+    // with the SAME key ends up first is left to the sort, and cannot matter: the packed key holds
+    // h, k, l and the hand exactly, so every run in a tie reduces to the same ASU reflection.
+    struct RunKey { uint64_t key; int32_t run; };
+    std::vector<RunKey> sorted;
+    sorted.reserve(n_run);
     for (int r = 0; r < n_run; ++r)
-        if (eligible[r]) idx.push_back(r);
-    std::sort(idx.begin(), idx.end(), [&](int32_t a, int32_t b) { return key[a] < key[b]; });
+        if (eligible[r]) sorted.push_back({key[r], r});
+    std::sort(sorted.begin(), sorted.end(), [](const RunKey &a, const RunKey &b) { return a.key < b.key; });
 
     group_h.clear(); group_k.clear(); group_l.clear();
     int n_groups = 0;
-    for (size_t j = 0; j < idx.size(); ++j) {
-        if (j == 0 || key[idx[j]] != key[idx[j - 1]]) {
-            const int r = idx[j];
+    for (size_t j = 0; j < sorted.size(); ++j) {
+        if (j == 0 || sorted[j].key != sorted[j - 1].key) {
+            const int r = sorted[j].run;
             const auto hkl = keygen(rawrun_h[r], rawrun_k[r], rawrun_l[r]);
             group_h.push_back(hkl.plus ? hkl.h : -hkl.h);
             group_k.push_back(hkl.plus ? hkl.k : -hkl.k);
             group_l.push_back(hkl.plus ? hkl.l : -hkl.l);
             ++n_groups;
         }
-        rawrun_group[idx[j]] = n_groups - 1;
+        rawrun_group[sorted[j].run] = n_groups - 1;
     }
 
     // Stamp the ASU-group id per obs from its raw hkl + the precomputed finiteness. For the GPU we build a
@@ -1001,15 +1007,27 @@ int RotationScaleMerge::ComputeAsuGroups(const HKLKeyGenerator &keygen) {
         for (auto &f : f1) f.get();
 
         // CSR starts + convert hist[t][g] into chunk t's write base for group g (exclusive prefix over t).
+        // Only the prefix ACROSS groups is sequential; the sum over the chunks inside one group, and
+        // handing that group's chunks their bases, are per-group and go on all threads. The single loop
+        // this replaces walked nt separate histograms one column at a time on one thread, which is nt
+        // cache misses per group. Same sums in the same order, so the same numbers.
         std::vector<int32_t> gstart(n_groups), gcount(n_groups);
+        const size_t hist_work = static_cast<size_t>(n_groups) * static_cast<size_t>(nt);
+        ParallelChunks(n_groups, ThreadsForWork(hist_work, nthreads), [&](int glo, int ghi) {
+            for (int g = glo; g < ghi; ++g) {
+                int c = 0;
+                for (int t = 0; t < nt; ++t) c += hist[t][g];
+                gcount[g] = c;
+            }
+        });
         int acc = 0;
-        for (int g = 0; g < n_groups; ++g) {
-            int base = acc;
-            gstart[g] = acc;
-            for (int t = 0; t < nt; ++t) { const int c = hist[t][g]; hist[t][g] = base; base += c; }
-            gcount[g] = base - acc;
-            acc = base;
-        }
+        for (int g = 0; g < n_groups; ++g) { gstart[g] = acc; acc += gcount[g]; }
+        ParallelChunks(n_groups, ThreadsForWork(hist_work, nthreads), [&](int glo, int ghi) {
+            for (int g = glo; g < ghi; ++g) {
+                int base = gstart[g];
+                for (int t = 0; t < nt; ++t) { const int c = hist[t][g]; hist[t][g] = base; base += c; }
+            }
+        });
 
         // Pass 2 (parallel): each chunk fills its obs into gperm at its per-group base (stable).
         std::vector<int32_t> gperm(acc);
@@ -2594,16 +2612,28 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
                 if (a.dev2 != b.dev2) return a.dev2 < b.dev2;
                 return a.d < b.d;
             };
-            const std::function<void(size_t, size_t, int, int)> split =
-                    [&](size_t lo, size_t hi, int b0, int b1) {
+            // The two halves a split leaves behind are disjoint ranges whose contents the nth_element
+            // above them has already fixed, so a child sees the same sequence whether it runs before,
+            // after or beside its sibling and the bins come out the same either way. Run them side by
+            // side while the range is still worth a thread; below that they stay on this one.
+            constexpr size_t SPLIT_MIN_PARALLEL = 1u << 17;
+            const std::function<void(size_t, size_t, int, int, int)> split =
+                    [&](size_t lo, size_t hi, int b0, int b1, int depth) {
                         if (b0 >= b1) return;
                         const int mid = (b0 + b1) / 2;
                         const size_t k = static_cast<size_t>(mid) * per;
                         std::nth_element(smp.begin() + lo, smp.begin() + k, smp.begin() + hi, by_I2);
-                        split(lo, k, b0, mid);
-                        split(k, hi, mid + 1, b1);
+                        if (depth > 0 && hi - lo > SPLIT_MIN_PARALLEL) {
+                            auto left = std::async(std::launch::async,
+                                                   [&, lo, k, b0, mid, depth] { split(lo, k, b0, mid, depth - 1); });
+                            split(k, hi, mid + 1, b1, depth - 1);
+                            left.get();
+                        } else {
+                            split(lo, k, b0, mid, 0);
+                            split(k, hi, mid + 1, b1, 0);
+                        }
                     };
-            split(0, smp.size(), 1, n_bins);
+            split(0, smp.size(), 1, n_bins, nthreads > 1 ? 3 : 0);
             std::vector<double> bs2(n_bins), bI2(n_bins), bd2(n_bins);
             ParallelChunks(n_bins, ThreadsForWork(smp.size(), nthreads, 8 * 32768), [&](int blo, int bhi) {
                 for (int bin = blo; bin < bhi; ++bin) {