last minor loop optimization

2026-05-28 18:55:25 +02:00
parent 37ce96041d
commit 7f235ebb20
2 changed files with 595 additions and 0 deletions
@@ -0,0 +1,118 @@
+# 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.