Linus/variant lut

docs/optimization_lut_exp.md

@@ -0,0 +1,138 @@
+# Optimization: LUT for exp()
+
+- **Author**: Linus
+- **Predecessor**: variant_baseline.cpp
+- **File(s)**: variants/variant_lut_exp.cpp
+- **Status**: positive at small N (10-15%), tied at large N
+
+## What this changes
+
+Replaced the exp() calls in pedestrian_repulsion_dt and
+obstacle_repulsion_dt with a precomputed lookup table over the argument
+range [-40, 0]. The table stores exp at evenly spaced sample points and we
+interpolate with a cubic Hermite spline between them. Below -40 the value
+is below FP noise relative to the force accumulator, so we just return 0.
+
+The thing that makes this work cleanly: for exp(x), the derivative IS the
+function value (d/dx exp(x) = exp(x)). So at every sample point we know
+the exact slope without having to estimate it with a centered difference.
+That gives O(h^4) accuracy with just 2 LUT loads per call.
+
+## Code (relevant snippet)
+
+```cpp
+constexpr int EXP_LUT_SIZE = 69000;
+constexpr double EXP_LUT_MIN = -40.0;
+constexpr double EXP_LUT_MAX = 0.0;
+constexpr double EXP_LUT_INV_STEP = (EXP_LUT_SIZE - 1) / (EXP_LUT_MAX - EXP_LUT_MIN);
+
+static double exp_lut[EXP_LUT_SIZE];
+
+// Populated before main() via a static initializer.
+struct ExpLutInit {
+  ExpLutInit() {
+    const double step = (EXP_LUT_MAX - EXP_LUT_MIN) / (EXP_LUT_SIZE - 1);
+    for (int i = 0; i < EXP_LUT_SIZE; i++) {
+      exp_lut[i] = std::exp(EXP_LUT_MIN + i * step);
+    }
+  }
+};
+static ExpLutInit exp_lut_init_helper;
+
+inline double fast_exp(double x) {
+  if (x <= EXP_LUT_MIN) return 0.0;
+  if (x >= EXP_LUT_MAX) return 1.0;
+  double pos = (x - EXP_LUT_MIN) * EXP_LUT_INV_STEP;
+  int idx = static_cast<int>(pos);
+  if (idx >= EXP_LUT_SIZE - 1) return 1.0;
+  double t = pos - idx;
+  double y0 = exp_lut[idx];
+  double y1 = exp_lut[idx + 1];
+  // For exp, slope = value. In normalized parameter space, slope = h * y.
+  constexpr double h_step = (EXP_LUT_MAX - EXP_LUT_MIN) / (EXP_LUT_SIZE - 1);
+  double t2 = t * t;
+  double t3 = t2 * t;
+  double h00 =  2.0*t3 - 3.0*t2 + 1.0;
+  double h10 =  t3 - 2.0*t2 + t;
+  double h01 = -2.0*t3 + 3.0*t2;
+  double h11 =  t3 - t2;
+  return y0 * (h00 + h_step * h10) + y1 * (h01 + h_step * h11);
+}
+```
+
+Then in pedestrian_repulsion_dt and obstacle_repulsion_dt, the libm
+`exp(...)` is replaced with `fast_exp(...)`.
+
+## Benchmark
+
+Conditions: min of 3 runs, N_ITER=1000, -O3 -march=native -flto.
+
+On my laptop (Intel ULTRA 9 185H @ 2.2 GHz, WSL Ubuntu): about
+15% faster at N=50, 12% at N=100, 8% at N=200, then roughly tied at N=500
+and N=1000. The shrinking speedup matches a cache story (see Analysis).
+
+
+### x86
+
+**Machine**:
+
+```
+<paste machine info here>
+```
+
+**Frequency**: <e.g., 2 GHz>
+
+| N    | W(n) [Flops] | T(n) [Cycles] | Performance [Flops/Cycle] | Diff vs Baseline |
+| ---- | ------------ | ------------- | ------------------------- | ---------------- |
+| 50   |              |               |                           |                  |
+| 100  |              |               |                           |                  |
+| 200  |              |               |                           |                  |
+| 500  |              |               |                           |                  |
+| 1000 |              |               |                           |                  |
+
+### ARM (Apple)
+
+**Machine**:
+
+```
+<paste machine info here>
+```
+
+**Frequency**: <e.g., 3.2 GHz>
+
+| N    | W(n) [Flops] | T(n) [Cycles] | Performance [Flops/Cycle] | Diff vs Baseline |
+| ---- | ------------ | ------------- | ------------------------- | ---------------- |
+| 50   |              |               |                           |                  |
+| 100  |              |               |                           |                  |
+| 200  |              |               |                           |                  |
+| 500  |              |               |                           |                  |
+| 1000 |              |               |                           |                  |
+
+## Analysis
+
+The LUT replaces a libm exp (about 30 cycles) with a table lookup plus a
+small cubic polynomial (about 18 cycles). At small N the table fits in L2
+and the speedup is real, around 10-15%. As N grows the pedestrian array
+starts competing with the LUT for L2 cache, and the LUT lookups end
+up in L3 instead, which is about as expensive as the libm exp we were
+trying to replace. 
+
+About verify: the LUT introduces a small per-call error, and chaos in the
+simulation amplifies it over 4000 timesteps, but still passes our strict verification
+
+If we are more tolerant with the "error" its possible to use a much smaller LUT. Maybe we could try that out
+
+## Takeaway
+
+LUT for exp works at small N where the table fits in cache. At large N the
+cache pressure kills the gain. A smaller LUT (or range reduction into a
+tiny table) would probably do better across the whole N range without the
+cache penalty.
+
+## References
+
+- Tang, P. T. P. (1989). "Table-driven implementation of the exponential
+  function in IEEE floating-point arithmetic." ACM TOMS 15(2), 144-157.
+  https://dl.acm.org/doi/10.1145/63522.214389
+- Cubic Hermite spline. Wikipedia.
+  https://en.wikipedia.org/wiki/Cubic_Hermite_spline