diff --git a/rugnux/RugnuxFirstPass.cpp b/rugnux/RugnuxFirstPass.cpp
index a7b552acc..5605c45d0 100644
--- a/rugnux/RugnuxFirstPass.cpp
+++ b/rugnux/RugnuxFirstPass.cpp
@@ -352,6 +352,21 @@ std::optional<TwinCompositeEvidence> TwinCompositeAt(const CrystalLattice &prim,
                 return true;
         return false;
     };
+    // Whether a reflection is on a sub-lattice depends only on its indices mod n, so the spots are
+    // counted once per class of (h, k, l) mod n, and the forms are asked of the classes.
+    const int n_classes = n * n * n;
+    const auto class_hkl = [n](int c) { return std::array<int, 3>{c / (n * n), c / n % n, c % n}; };
+    std::vector<int> per_class(n_classes, 0);
+    for (const auto &m : hkl)
+        per_class[(mod(m[0]) * n + mod(m[1])) * n + mod(m[2])]++;
+    std::vector<std::vector<int>> classes_on(forms.size());   // the classes on each form's sub-lattice
+    std::vector<int> spots_on(forms.size(), 0);
+    for (size_t i = 0; i < forms.size(); i++)
+        for (int c = 0; c < n_classes; c++)
+            if (on(forms[i], class_hkl(c))) {
+                classes_on[i].push_back(c);
+                spots_on[i] += per_class[c];
+            }
     // The congruent pair whose union holds the most spots, the stronger domain first.
     std::optional<TwinCompositeEvidence> best;
     for (size_t i = 0; i < forms.size(); i++)
@@ -359,10 +374,12 @@ std::optional<TwinCompositeEvidence> TwinCompositeAt(const CrystalLattice &prim,
             if (!congruent(cells[i], cells[j]))
                 continue;
             TwinCompositeEvidence t{.f1 = forms[i], .f2 = forms[j]};
-            for (const auto &m : hkl) {
-                const bool o1 = on(forms[i], m), o2 = on(forms[j], m);
-                (o1 && o2 ? t.both : o1 ? t.only1 : o2 ? t.only2 : t.neither)++;
-            }
+            for (const int c : classes_on[i])
+                if (on(forms[j], class_hkl(c)))
+                    t.both += per_class[c];
+            t.only1 = spots_on[i] - t.both;
+            t.only2 = spots_on[j] - t.both;
+            t.neither = static_cast<int>(hkl.size()) - spots_on[i] - spots_on[j] + t.both;
             if (!best || t.neither < best->neither)
                 best = t;
         }