std::array<float32x4_t, 4> v;
for (size_t j = 0; j < 4; ++j) v[j] = vld1q_f32(ptr + i + 4 * j);
for (size_t j = 0; j < 4; ++j) vst1q_f32(ptr + i + 4 * j, f(v[j]));
}

for (size_t i = 0; i < n; i += 16) {
std::array<uint32_t, 4> sk;
for (size_t j = 0; j < 4; ++j) sk[j] = s[i / 4 + j];
std::array<size_t, 4> off{};
for (size_t j = 1; j < 4; ++j) off[j] = off[j - 1] + (sk[j - 1] >> 4);
std::array<float32x4_t, 4> v;
for (size_t j = 0; j < 4; ++j) v[j] = vld1q_f32(ptr + off[j]);
std::array<float32x4_t, 4> a;
for (size_t j = 0; j < 4; ++j) a[j] = vld1q_f32(dst + i + 4 * j);
std::array<uint8x16_t, 4> index;
for (size_t j = 0; j < 4; ++j) index[j] = vld1q_u8(exp_tbl[sk[j] & 15].data());
std::array<uint8x16x2_t, 4> tbl;
for (size_t j = 0; j < 4; ++j) tbl[j] = {{vreinterpretq_u8_f32(v[j]), vreinterpretq_u8_f32(a[j])}};
for (size_t j = 0; j < 4; ++j) vst1q_f32(dst + i + 4 * j, vreinterpretq_f32_u8(vqtbl2q_u8(tbl[j], index[j])));
ptr += off[3] + (sk[3] >> 4);
}
return size;
}

void tiled_detour(float* dst, const size_t n, auto f, auto mask) {
const auto w = vld1q_u32(weights.data());
for (size_t i = 0; i < n; i += tile)
detour<false>(dst + i, tile, w, f, mask);
}

compress is the same as
in [my previous post](https://www.reddit.com/r/cpp/comments/1uqscnr/stream_compaction_on_neon_vectorizing_copy_if_by/).

expand_table: for true lanes it selects the next element from the compressed register (bytes from `[0, 15]`),
and for false lanes, selects the same bytes + 16.
Then tbl2 on {processed, original}, the same trick as in compress basically

Also:

- expand is skipped for free on empty tiles
- `s` is saved for free to avoid recomputing addv
- cache-sized tiling.
- Empirically `tile=4096` is optimal.

BSL speed doesn't depend on density: `sqrt 32.1`, `sin35 - 9.7`, `pow10 1.02`.

| thd | 0 | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 | 0.6 | 0.7 | 0.8 | 0.9 | 1 |
|:-------------|--------:|--------:|--------:|--------:|--------:|--------:|-------:|-------:|-------:|-------:|------:|
| detour sqrt | 38.03\* | 18.81 | 17.78 | 16.82 | 15.98 | 15.21 | 14.51 | 13.85 | 13.24 | 12.71 | 12.30 |
| detour sin35 | 38.16\* | 16.70\* | 14.46\* | 12.71\* | 11.36\* | 10.25\* | 9.35 | 8.59 | 7.95 | 7.39 | 6.95 |
| detour pow10 | 38.09\* | 6.83\* | 4.13\* | 2.96\* | 2.31\* | 1.89\* | 1.60\* | 1.39\* | 1.22\* | 1.10\* | 0.99 |

For cheap sqrt BSL is always better (except thd = 0, obviously).
And for very expensive pow, detour is better (except thd = 1, of course).

## When detour wins

Define:

- `B` = BSL(`vld + vbsl + vst`) overhead per register.
- `D` = detour(compress + expand) overhead per register.
- `T` = cost of `f` per register (we assume `cost of f >> cost of mask`, affects only calibration accuracy).
- `d` = fraction of selected elements

BSL applies `f` to every register. detour applies it only to `d` of them, so it saves `T * (1 - d)`.
Detour wins when the saving outweighs `D - B`.

`T * (1 - d) > D - B`

So detour pays off when `d < d_max = 1 - (D - B) / T`.
`B` and `D` depend only on hardware, so let's premeasure them (I have `B = 0.3`, `D = 0.845` ns/register)

`T` we measure over the first few tiles, timing BSL. And from it, we also find `d_max`:

template<size_t tile>
double bsl_calibrate(float* dst, const size_t len, auto f, auto mask) {
double ns = 1e18;
for (size_t i = 0; i < len; i += tile) {
const auto st = std::chrono::high_resolution_clock::now();
bsl<false>(dst, tile, f, mask);
const auto ed = std::chrono::high_resolution_clock::now();
ns = std::min(ns, std::chrono::duration<double, std::nano>(ed - st).count());
dst += tile;
}
const double t = std::max(1e-9, ns / (tile / 4.0) - B);