# 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