CPU profile fit: the Gaussian's cells in a loop that vectorises, output unchanged

Building each reflection's Gaussian profile was a third of the CPU pass B's samples (cytochrome C),
one scalar BraggExp per cell. BraggExp branched to its two ends before the arithmetic, and the
branches - with GCC's -ftrapping-math, which forbids running float arithmetic on a path that does not
need it - kept every loop over it scalar. It now tells the ends apart on the bits of x, an integer
comparison, runs the same arithmetic on x (or on 0 out of range) and masks the result: the same bits for
every float but NaN, checked over all 2^32 inputs against the old function, scalar and vectorised.
The cells are filled a row at a time and summed afterwards in the same order, so the loop vectorises
(8 lanes, AVX2) and the profile is unchanged. BraggScale2 chooses its operands instead of a branch;
scaling by 1 is exact. The Metal kernel shares the header and gets the same arithmetic.

cytochrome C --no-gpu: p.hkl and MTZ data identical; user cycles -7% against the previous commit on
the same (loaded, load 16-32) machine.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

image_analysis/bragg_integration/BraggExp.h

@@ -18,24 +18,38 @@
 
 #ifdef __METAL_VERSION__
 inline float BraggBitsToFloat(int i) { return as_type<float>(i); }
+inline int BraggFloatToBits(float f) { return as_type<int>(f); }
 #else
 #include <bit>
 #include <cstdint>
 inline float BraggBitsToFloat(int32_t i) { return std::bit_cast<float>(i); }
+inline int32_t BraggFloatToBits(float f) { return std::bit_cast<int32_t>(f); }
 #endif
 
-// p * 2^k, rounded once where the result is subnormal.
+// p * 2^k, rounded once where the result is subnormal: below k = -125 the scaling is split in two. The
+// two cases differ only in their operands - the normal one scales by 1 the second time, which is exact -
+// so a loop over it has no branch and vectorises.
 inline float BraggScale2(float p, int k) {
-    if (k < -125)
-        return (p * BraggBitsToFloat((k + 100 + 127) << 23)) * BraggBitsToFloat((-100 + 127) << 23);
-    return p * BraggBitsToFloat((k + 127) << 23);
+    const bool subnormal = k < -125;
+    const float first = BraggBitsToFloat(((subnormal ? k + 100 : k) + 127) << 23);
+    const float second = subnormal ? BraggBitsToFloat((-100 + 127) << 23) : 1.0f;
+    return (p * first) * second;
 }
 
 // Following Cody & Waite (1980) Software Manual for the Elementary Functions; Dekker (1971) Numer. Math.
 // 18, 224-242; Knuth, The Art of Computer Programming Vol. 2, 4.2.2 (2Sum).
-inline float BraggExp(float x) {
-    if (x < -104.0f) return 0.0f;                                  // below half the smallest subnormal
-    if (x > 88.72283f) return BraggBitsToFloat(0x7f800000);        // overflows
+// The two ends are told apart on the bits of x and applied to the bits of the result, with integer masks
+// rather than branches or selects: the compiler moves arithmetic behind a branch or a select that only
+// one side needs, and may not run float arithmetic on a path that does not need it, which either way
+// keeps a loop over BraggExp from vectorising. Out of range the arithmetic runs on 0 and its result is
+// masked off; in range it is the same arithmetic on the same x. (A negative float's bits are a negative
+// integer that grows as the float falls; a positive one's grow with it.)
+inline float BraggExp(float x_in) {
+    const int bits = BraggFloatToBits(x_in);
+    const bool under = (bits < 0) & (bits > BraggFloatToBits(-104.0f));   // below half the smallest subnormal
+    const bool over = bits > BraggFloatToBits(88.72283f);                  // overflows
+    const int in_range = -static_cast<int>(!(under | over));              // every bit set, or none
+    const float x = BraggBitsToFloat(bits & in_range);
     const float shifter = 12582912.0f;                             // 1.5 * 2^23: adding it rounds to an integer
     const float kf = (x * 1.44269504f + shifter) - shifter;
     const int k = static_cast<int>(kf);
@@ -59,5 +73,6 @@ inline float BraggExp(float x) {
     const float e1 = (rr - (s1 - bv1)) + (hz - bv1);
     const float t1 = 1.0f + s1;                                    // 1 + s1 = t1 + e3 (Fast2Sum, |s1| < 1)
     const float e3 = (1.0f - t1) + s1;
-    return BraggScale2(t1 + (e3 + (e1 + c)), k);
+    const float e = BraggScale2(t1 + (e3 + (e1 + c)), k);
+    return BraggBitsToFloat((BraggFloatToBits(e) & in_range) | (static_cast<int>(over) * 0x7f800000));
 }

image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp

@@ -486,13 +486,14 @@ BraggFitOut BraggIntegrationEngineCPU::FitProfile(const Sampler &img, const Owne
             prefetch_row(img.address(static_cast<size_t>(y) * W + gx0),
                          img.address(static_cast<size_t>(y) * W + gx1));
         cells.resize(static_cast<size_t>(Gf) * Gf);
+        // The cells a row at a time, a loop that vectorises; then their sum, in the same order.
+        for (int dy = -Rf; dy <= Rf; ++dy) {
+            float *row = cells.data() + static_cast<size_t>(dy + Rf) * Gf;
+            for (int dx = -Rf; dx <= Rf; ++dx)
+                row[dx + Rf] = BraggGaussian(dx, dy, s.fx, s.fy, s.ux, s.uy, s.inv_2s2r, s.inv_2s2t);
+        }
         float gs = 0.0f;
-        for (int dy = -Rf; dy <= Rf; ++dy)
-            for (int dx = -Rf; dx <= Rf; ++dx) {
-                const float g = BraggGaussian(dx, dy, s.fx, s.fy, s.ux, s.uy, s.inv_2s2r, s.inv_2s2t);
-                cells[(dy + Rf) * Gf + (dx + Rf)] = g;
-                gs += g;
-            }
+        for (const float g : cells) gs += g;
         const float inv_gs = 1.0f / gs;
         for (float &p : cells) p *= inv_gs;
         P = cells.data();