From 79a80f086cdcf445e956d49b08696b94d65ffb4f Mon Sep 17 00:00:00 2001
From: Filip Leonarski <filip.leonarski@psi.ch>
Date: Wed, 4 Nov 2026 14:42:29 +0100
Subject: [PATCH] Twin-composite hypothesis: count the spots per class of hkl
 mod n, not per pair of forms

Whether a reflection lies on an index-n sub-lattice depends only on its indices
mod n, so TwinCompositeAt histograms the validation spots over the n^3 classes
once and counts each congruent pair of forms over the classes on the first form's
sub-lattice. The counts are the same integers as before, so the decision and the
log are unchanged; the index-2..10 loop drops from ~0.45 s to ~10 ms a cell, and
it runs for the standing cell and every multiple in every first pass
(on a median open set ~3 s of the first passes, single-threaded).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01UqXYdL4qBGYRiSXSGXfCLu
---
 rugnux/RugnuxFirstPass.cpp | 25 +++++++++++++++++++++----
 1 file changed, 21 insertions(+), 4 deletions(-)

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;
         }