# Optimization: ghost-self trick + r_sum fma fuse

- **Author**: Linus
- **Predecessor**: variant_float_inlining_SoA_AVX2_polyexp_rsqrt.cpp
- **File(s)**: variants/variant_float_inlining_SoA_AVX2_polyexp_rsqrt_loopopt.cpp
- **Status**: positive

## What this changes

Kill the per-block branchless mask in the j-loop. Before the loop we save
`pos[i]`, park it at `1e9`, then broadcast i's real position from the saved
copy. The j==i lane now reads like a padding lane (huge distance ->
~1e-37 force), so it self-cancels with no mask. That removes the
`m_neq_i`/`m_in_range` cmps, the active-AND, and the two safety-ANDs on
`f_x` / `f_y` before the accumulator fmadd. Restore `pos[i]` right
after the j-loop so the rest of the simulation sees the real positions.

Second change in the same variant: fold the two sqrt-recovery muls and the
add for `r_sum = |r_ab| + |r_ab2|` into one mul + one fmadd. The two norms
are only used in this one place, so we don't need them as named values.

## Code (relevant snippet)

```cpp
// ghost-self trick: park pedestrian i at 1e9 for its own j-loop, broadcasting i's
// real position from a saved copy. then the j==i lane reads 1e9 like a padding lane
// -> ~1e-37 force, self-cancels with no mask. restored right after the loop.
const float saved_pix = pos_x[i];
const float saved_piy = pos_y[i];
const __m256 pix = _mm256_set1_ps(saved_pix);
const __m256 piy = _mm256_set1_ps(saved_piy);
const __m256 eix = _mm256_set1_ps(e_x[i]);
const __m256 eiy = _mm256_set1_ps(e_y[i]);
pos_x[i] = 1e9f;
pos_y[i] = 1e9f;

__m256 acc_x = _mm256_setzero_ps();
__m256 acc_y = _mm256_setzero_ps();

for (int j = 0; j < n_p_pad; j += 8) {
  // no more j_idx / m_neq_i / m_in_range / active

  __m256 pjx = _mm256_loadu_ps(&pos_x[j]);
  __m256 pjy = _mm256_loadu_ps(&pos_y[j]);
  __m256 sbe_jx = _mm256_loadu_ps(&sb_e_x[j]);
  __m256 sbe_jy = _mm256_loadu_ps(&sb_e_y[j]);
  __m256 sb_jsq = _mm256_loadu_ps(&sb_sq[j]);

  __m256 r_ab_x = _mm256_sub_ps(pix, pjx);
  __m256 r_ab_y = _mm256_sub_ps(piy, pjy);
  __m256 r_ab_n2 = _mm256_fmadd_ps(r_ab_x, r_ab_x, _mm256_mul_ps(r_ab_y, r_ab_y));
  r_ab_n2 = _mm256_max_ps(r_ab_n2, v_eps_sq);
  __m256 inv_r_ab_n = _mm256_rsqrt_ps(r_ab_n2);

  __m256 r_ab2_x = _mm256_sub_ps(r_ab_x, sbe_jx);
  __m256 r_ab2_y = _mm256_sub_ps(r_ab_y, sbe_jy);
  __m256 r_ab2_n2 = _mm256_fmadd_ps(r_ab2_x, r_ab2_x, _mm256_mul_ps(r_ab2_y, r_ab2_y));
  r_ab2_n2 = _mm256_max_ps(r_ab2_n2, v_eps_sq);
  __m256 inv_r_ab2_n = _mm256_rsqrt_ps(r_ab2_n2);

  // r_sum fuse: two recovery muls + add -> one mul + one fmadd
  __m256 r_sum = _mm256_fmadd_ps(r_ab2_n2, inv_r_ab2_n,
                                 _mm256_mul_ps(r_ab_n2, inv_r_ab_n));

  // ... force-scaling unchanged ...

  // no f_x_safe / f_y_safe AND before the accumulator
  acc_x = _mm256_fmadd_ps(f_x, w_val, acc_x);
  acc_y = _mm256_fmadd_ps(f_y, w_val, acc_y);
}

// restore the real positions
pos_x[i] = saved_pix;
pos_y[i] = saved_piy;
```

## Benchmark

Conditions: N_ITER=1000, -O3 -march=native -flto.

### x86

**Machine**: TBD

**Frequency**: TBD

| N    | W(n) [Flops] | T(n) [Cycles] | Performance [Flops/Cycle] | Diff vs Predecessor |
| ---- | ------------ | ------------- | ------------------------- | ------------------- |
| 50   |              |               |                           |                     |
| 100  |              |               |                           |                     |
| 200  |              |               |                           |                     |
| 500  |              |               |                           |                     |
| 1000 |              |               |                           |                     |

### ARM (Apple)

no avx2 here

## Analysis

ILP tricks (j-unroll, i-tiling) didn't help earlier, so the win had to come
from fewer FP ops, not from running more in parallel. Ghost-self + r_sum
fuse just drop ops. Asm diff (gcc 13.3, -O3 -march=native):

```
                vcmpps  vandps  vfmadd  vrsqrtps
predecessor:       4      4      33       3
variant:           2      1      33       3
```

The mask is gone, FMA/rsqrt counts unchanged. The AND that mattered was
the one between `f_x` and the accumulator fmadd, since that's on the
loop-carried chain. 

## Takeaway

Less work, same precision. The mask was there to handle j==i but the
ghost-self lane handles it for free.