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)

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

4.3 KiB Raw Blame History