# 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.