SoA implementation for float and AVX2 includind the polyexp fix

docs/optimization_float_inlining_SoA_AVX2.md

@@ -0,0 +1,98 @@
+# Optimization: AVX2 vectorization of the ped-ped loop (float, 8-wide)
+
+- **Author**: Linus
+- **Predecessor**: variant_float_inlining_SoA.cpp
+- **File(s)**: variants/variant_float_inlining_SoA_AVX2.cpp
+- **Status**: TBD
+
+## What this changes
+
+Float port of the double `variant_inlining_7_SoA_AVX2`. Inner ped-ped j
+loop runs on `__m256` (8 floats) instead of the double version's `__m256d`
+(4 doubles). 8 ped-ped interactions per j step. SoA layout from the
+predecessor makes the j-side loads contiguous.
+
+Same masking trick as the double variant: build an `active` mask from
+`j_idx != i` AND `j_idx < n_p`, blend a safe `1.0` into the denominator
+where `active` is zero so the i==j and padding lanes don't NaN out, then
+`vandps` the per-lane `f_x`/`f_y` against `active` before accumulating
+into `acc_x`/`acc_y`.
+
+n_p is padded to a multiple of 8 (`n_p_pad = (n_p + 7) & ~7`) so the j
+loop has no scalar tail. Padding positions are set to `1e9f` and `sb_e`,
+`sb_sq` to `0` so any unmasked contribution from those lanes would still
+be ~zero through `expf(very negative) ≈ 0`.
+
+EPS clamps from the float scalar predecessor are kept and lifted to SIMD
+form (`_mm256_max_ps(r_ab_n, v_eps)` etc.). Same correctness reason as
+the scalar version: float precision around close encounters can collapse
+norms to zero.
+
+The `expf` calls are still scalar inside the j loop: 8 of them per j
+step, manually unrolled into a scratch `tmp_arg[8]`/`tmp_val[8]` array.
+gcc refused to unroll the loop on its own (it serialised all 8 calls on
+`xmm0`), the manual unroll at least lets them issue independently. Kill
+these next with the polyexp variant.
+
+## Code (relevant snippet)
+
+```cpp
+for (int j = 0; j < n_p_pad; j += 8) {
+  __m256 j_idx = _mm256_add_ps(_mm256_set1_ps((float)j), lane_off);
+  __m256 active = _mm256_and_ps(
+      _mm256_cmp_ps(j_idx, v_i,  _CMP_NEQ_OQ),
+      _mm256_cmp_ps(j_idx, v_np, _CMP_LT_OQ));
+
+  __m256 pjx    = _mm256_loadu_ps(&pos_x[j]);
+  __m256 pjy    = _mm256_loadu_ps(&pos_y[j]);
+  // ... sb_e_x, sb_e_y, sb_sq
+
+  // body: r_ab, r_ab2, two_b, scalar expf x 8, smc, sm_r_ab, sm_r_ab2,
+  // f_x, f_y, w_val branch, AND with active, fmadd into acc_x/acc_y.
+}
+```
+
+See `variants/variant_float_inlining_SoA_AVX2.cpp` for the full inner
+loop.
+
+## Benchmark
+
+Conditions: min of N runs (TBD), N_ITER=1000, -O3 -march=native -flto.
+Diff is vs `variant_float_inlining_SoA`.
+
+### x86
+
+**Machine**: TBD
+
+**Frequency**: TBD
+
+| N    | W(n) [Flops] | T(n) [Cycles] | Performance [Flops/Cycle] | Diff vs Predecessor |
+| ---- | ------------ | ------------- | ------------------------- | ------------------- |
+| 50   |              |               |                           |                     |
+| 100  |              |               |                           |                     |
+| 200  |              |               |                           |                     |
+| 500  |              |               |                           |                     |
+
+### ARM (Apple)
+
+Not applicable, this variant uses AVX2 intrinsics.
+
+## Analysis
+
+Same shape of win as the double AVX2 variant, just at 8-wide. Scalar
+`expf` is still the dominant cost, polyexp kills it next.
+
+One catch: this variant only hits its real performance if FTZ is
+enabled in `main-float.cpp` (that gets added on the polyexp MR). Without
+it the active-lane chain produces subnormal intermediates (the `smc =
+numer/denom` step is the one that lands there for typical mid-range pair
+distances) and every following FMA stalls 100+ cycles in microcode. With
+FTZ on it runs as expected, without it benchmarks look 5-10x slow.
+
+Verify passes against `baseline_float_n60_i500.csv` at 1cm tolerance.
+Local max drift ~3mm.
+
+## Takeaway
+
+ Polyexp is the next step. Note:
+benchmarks are only meaningful with FTZ/DAZ on in main-float.cpp.