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)

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

5.0 KiB Raw Blame History