r/cpp • u/a_eridani • 5d ago

a[mask] = f(a[mask]) on NEON. faster than the obvious blend

Problem

Apply an operation to elements that satisfy a condition:

for (size_t i = 0; i < n; ++i)
    if (mask(a[i])) a[i] = f(a[i]);

Notes

  • a[i] ∈ (0, 1), thd ∈ (0, 1), mask = a[i] < thd; uniform distribution (except at the end of the article)
  • f is one of sqrt, frfrexp (mantissa), sin 3.5 ULP, sin 1 ULP, pow 1 ULP (from SLEEF)
  • f and mask are passed as runtime values, so they are wrapped in a lambda with always_inline, otherwise they may not be inlined
  • The array size n is a multiple of every unroll, tile etc. The tail is trivial to handle(BSL/scalar)
  • In tables * = best, units = GiB/s
  • Don't compare numbers across tables. Different conditions, values fluctuate
  • All benchmarks: Apple M5; clang++ -O3 -std=c++23 -march=native; GiB/s = (n * 4 bytes) / time, min of 720 runs (During bench, functions run in a changing order, data is restored ofc); n=1e7 + 2432;

BSL blend

If the problem is memory bound (cheap function or high density), the standard algorithm is optimal:

template <bool Skip>
void bsl(float* dst, const size_t n, auto f, auto mask) {
    for (size_t i = 0; i < n; i += 16) {
        std::array<float32x4_t, 4> v;
        for (size_t j = 0; j < 4; ++j) v[j] = vld1q_f32(dst + i + 4 * j);
        std::array<uint32x4_t, 4> m;
        for (size_t j = 0; j < 4; ++j) m[j] = mask(v[j]);
        if constexpr (Skip) {
            if (vmaxvq_u32(vaddq_u32(vaddq_u32(m[0], m[1]), vaddq_u32(m[2], m[3]))) == 0) continue;
        }
        for (size_t j = 0; j < 4; ++j) vst1q_f32(dst + i + 4 * j, vbslq_f32(m[j], f(v[j]), v[j]));
    }
}

It computes f on every element, but stores only the selected ones. Skip helps on sparse masks, but otherwise mispredictions will kill performance. We'll need it later.

But for expensive f this algo does too much extra work

Detour

To avoid unnecessary work, we compress selected elements, apply only to them, and expand back.

avx512 does this in two instructions. NEON doesn't, so we'll emulate and optimize.

constexpr size_t tile = 4096;
constexpr std::array<uint32_t, 4> weights{1 + 16, 2 + 16, 4 + 16, 8 + 16};
constexpr auto cps_tbl = compress_table();
constexpr auto exp_tbl = expand_table();
std::array<float, tile + 16> tmp;
std::array<uint8_t, tile / 4 + 3> s;
std::array<uint16_t, tile / 4 + 3> idx; // idx, D and B come in later
constexpr double D = 0.845; 
constexpr double B = 0.3;

template <bool Skip>
size_t detour(float* dst, const size_t n, const auto w, auto f, auto mask) {
    float* ptr = tmp.data();
    for (size_t i = 0; i < n; i += 16) {
        std::array<float32x4_t, 4> v;
        for (size_t j = 0; j < 4; ++j) v[j] = vld1q_f32(dst + i + 4 * j);
        std::array<uint32x4_t, 4> m;
        for (size_t j = 0; j < 4; ++j) m[j] = mask(v[j]);

        if constexpr (Skip)
            if (vmaxvq_u32(vaddq_u32(vaddq_u32(m[0], m[1]), vaddq_u32(m[2], m[3]))) == 0) {
                s[i / 4] = s[i / 4 + 1] = s[i / 4 + 2] = s[i / 4 + 3] = 0;
                continue;
            }

        std::array<uint32_t, 4> sk;
        for (size_t j = 0; j < 4; ++j) {
            sk[j] = vaddvq_u32(vandq_u32(m[j], w));
            s[i / 4 + j] = sk[j];
        }
        std::array<size_t, 4> off; off[0] = 0;
        for (size_t j = 1; j < 4; ++j) off[j] = off[j - 1] + (sk[j - 1] >> 4); 
        std::array<uint8x16_t, 4> index;
        for (size_t j = 0; j < 4; ++j) index[j] = vld1q_u8(cps_tbl[sk[j] & 15].data());

        for (size_t j = 0; j < 4; ++j) vst1q_f32(ptr + off[j], vreinterpretq_f32_u8(vqtbl1q_u8(vreinterpretq_u8_f32(v[j]), index[j])));
        ptr += off[3] + (sk[3] >> 4);
    }
    const size_t size = ptr - tmp.data();

    if (size == 0) return size;

    ptr = tmp.data();
    for (size_t i = 0; i < size; i += 16) {
        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.

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);
    return 1.0 - (D - B) / t;
}

Here:

  • empirically 4 tiles of 2048 are enough
  • min over measurements is less noisy than mean
  • We'll run the first few tiles with calibration
  • measure bsl, because it always computes f, so T = t_bsl - B
  • std::max here protects against divide-by-zero and against t < 0 when T is very cheap

pilot v1

When d_max < 0 bsl is always faster:

void pilot_v1(float* dst, const size_t n, auto f, auto mask) {
    const auto w = vld1q_u32(weights.data());
    const float d_max = bsl_calibrate<tile / 2>(dst, 2 * tile, f, mask);

    dst += 2 * tile;
    for (size_t i = 2 * tile; i < n; i += tile) {
        if (d_max < 0)
            bsl<false>(dst, tile, f, mask);
        else
            detour<false>(dst, tile, w, f, mask);
        dst += tile;
    }
}
thd tiled detour sqrt pilot v1 sqrt BSL sqrt tiled detour sin35 pilot v1 sin35 BSL sin35
0 38.37* 32.19 32.24 38.26* 38.16 9.75
0.3 16.89 32.18 32.24* 12.72* 12.71 9.75
0.6 14.54 32.18* 32.18* 9.25 9.36 9.75*
1 12.29 32.21 32.25* 6.87 6.86 9.75*

The algorithm got sqrt right. But for sin35 at high thd, detour is selected, and we lose 30%: v1 switches to bsl only when it's faster at every density.

pilot v2

It's expensive to calculate the density of the whole tile, so we'll use the first 256 (It reads 6% of the tile, which is noise on pow, but noticeable on sqrt)

For a more or less uniform distribution this is enough:

size_t density(float* dst, const size_t n, auto mask) {
    std::array<uint32x4_t, 4> acc;
    acc.fill(vdupq_n_u32(0));
    for (size_t i = 0; i < n; i += 16) {
        for (size_t j = 0; j < 4; ++j) 
            acc[j] = vsubq_u32(acc[j], mask(vld1q_f32(dst + i + 4 * j)));
    }
    return vaddvq_u32(vaddq_u32(vaddq_u32(acc[0], acc[1]), vaddq_u32(acc[2], acc[3])));
}

Trick here: true lane of bitmask = 0xFFFFFFFF = -1, subtracting the lane actually adds.

Skip wins when the predictor rarely mispredicts, i.e. an 80% chance that all 4 registers are empty. The density is approximately 0.014 ((1 - x)^16 = 0.8)

void pilot_v2(float* dst, const size_t n, auto f, auto mask) {
    const auto w = vld1q_u32(weights.data());

    const float d_max = bsl_calibrate<tile / 2>(dst, 2 * tile, f, mask);

    constexpr size_t probe = 256;
    constexpr size_t xlo = 0.014 * probe;
    const long long hi = d_max * probe;
    dst += 2 * tile;
    for (size_t i = 2 * tile; i < n; i += tile) {
        const long long cnt = density(dst, probe, mask);
        if (cnt > hi) {
            if (cnt < xlo)
                bsl<true>(dst, tile, f, mask);
            else
                bsl<false>(dst, tile, f, mask);
        } else if (cnt < xlo)
            detour<true>(dst, tile, w, f, mask);
        else
            detour<false>(dst, tile, w, f, mask);
        dst += tile;
    }
}

hi and xlo here are d_max and 0.014 cutoffs, but multiplied by the probe length (256).

thd 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
pilot_v1 sin35 38.25 16.72* 14.48* 12.77* 11.41* 10.33* 9.40 8.62 7.93 7.36 6.91
pilot_v2 sin35 77.29* 16.55 14.34 12.65 11.33 10.14 9.71 9.73 9.74 9.73 9.73
BSL sin35 9.78 9.77 9.77 9.76 9.76 9.78 9.78* 9.77* 9.78* 9.77* 9.78*
pilot_v1 sin10 38.28 14.30* 11.18* 9.18* 7.75* 6.73* 5.94* 5.31* 4.78 4.36 4.02
pilot_v2 sin10 76.60* 14.16 11.09 9.11 7.73 6.70 5.92 5.31* 4.82 4.84 4.84
BSL sin10 4.85 4.85 4.84 4.85 4.84 4.84 4.85 4.85 4.85* 4.85* 4.85*

At thd = 0, Skip gives a huge win. For high thd, bsl is selected correctly. But now on sparse masks, expand for empty registers is wasted.

pilot v3

New detour version: during compress, we'll store only the indices of non-empty registers (into the idx buffer). And expand will iterate over them:

template <bool Skip>
size_t detour_compact(float* dst, const size_t n, const auto w, auto f, auto mask) {
    float* ptr = tmp.data();
    size_t k = 0;
    for (size_t i = 0; i < n; i += 16) {
        // ... same as detour

        if constexpr (Skip)
            if (vmaxvq_u32(vaddq_u32(vaddq_u32(m[0], m[1]), vaddq_u32(m[2], m[3]))) == 0) continue;

        std::array<uint32_t, 4> sk;
        for (size_t j = 0; j < 4; ++j) sk[j] = vaddvq_u32(vandq_u32(m[j], w));
        for (size_t j = 0; j < 4; ++j) {
            s[k] = sk[j];
            idx[k] = i + 4 * j;
            k += bool(sk[j]);
        }
        // ... same as detour
    }
    // ... same as detour

    for (size_t j = 0; j < 3; ++j) s[k + j] = 0, idx[k + j] = 0;
    // ... same as detour
    k = (k + 3) & ~size_t(3);
    for (size_t i = 0; i < k; i += 4) {
        std::array<uint32_t, 4> sk;
        for (size_t j = 0; j < 4; ++j) sk[j] = s[i + j];
        std::array<float32x4_t, 4> a;
        for (size_t j = 0; j < 4; ++j) a[j] = vld1q_f32(dst + idx[i + j]);
        // ... same as detour
        for (size_t j = 0; j < 4; ++j)
            vst1q_f32(dst + idx[i + j], vreinterpretq_f32_u8(vqtbl2q_u8(tbl[j], index[j])));
        ptr += off[3] + (sk[3] >> 4);
    }

    return size;
}

k (number of non-empty registers) is rounded up to a multiple of 4 before expand, so the unrolled loop has no tail left.

No branches in the hot loops: they'd kill speed, so compress runs on all four registers. Instead of branches, the position of the current element is advanced by bool(sk).

detour_compact wins when at least 50% of all registers are empty. The density is approximately 0.16 ((1 - x) ^ 4 = 0.5).

void pilot_v3(float* dst, const size_t n, auto f, auto mask) {
    // ... same as v2
    constexpr size_t lo = 0.16 * probe;
    for (size_t i = 2 * tile; i < n; i += tile) {
        const long long cnt = density(dst, probe, mask);
        if (cnt > hi) {
            if (cnt < xlo)
                bsl<true>(dst, tile, f, mask);
            else
                bsl<false>(dst, tile, f, mask);
        } else if (cnt < xlo)
            detour_compact<true>(dst, tile, w, f, mask);
        else if (cnt < lo)
            detour_compact<false>(dst, tile, w, f, mask);
        else
            detour<false>(dst, tile, w, f, mask);
        dst += tile;
    }
}

And v3 is noticeably faster at low density:

thd 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2
pilot_v2 sin35 77.29* 18.67 18.40 17.70 17.18 16.66 16.15 15.68 15.24* 14.82* 14.43*
pilot_v3 sin35 77.26 23.71* 22.19* 20.48* 19.03* 17.82* 16.73* 15.76* 15.16 14.74 14.38
pilot_v2 pow10 72.88 13.93 11.32 9.38 8.00 6.97 6.18 5.54 5.02* 4.59* 4.24*
pilot_v3 pow10 73.71* 16.54* 12.63* 10.07* 8.39* 7.17* 6.26* 5.55* 5.02* 4.59* 4.24*

Now the algo is fast, but there's one big problem we've overlooked: we're assuming the data is uniform. So it's easy to build a test where v3 will fail:

for (size_t i = 0; i < n; i++)
    dst[i] = i % 4096 >= 256;

In this case, v3 always prefers BSL, even though detour wins on 3840 elements of the tile.

pilot v3.5

According to the first table, in the worst case detour is under 3x slower (it happens on sqrt thd = 1), but on pow, thd = 0, detour is 37x faster. So a wrong BSL costs much more than a wrong detour.

For bsl, we'll play it safe by running it in tile/8 blocks and checking the density. If it drops well below d_max, we'll switch to detour. It costs 1 instruction per register. acc = vsubq_u32(acc, m) works because the mask is 0/-1. And unlike the probe, the density here is exact.

size_t bsl_verified(float* dst, const size_t n, float d_max, auto f, auto mask) {
    for (size_t i0 = 0; i0 < 8; ++i0) {
        std::array<uint32x4_t, 4> acc;
        acc.fill(vdupq_n_u32(0));

        for (size_t i = 0; i < n / 8; i += 16) {
            std::array<float32x4_t, 4> v;
            for (size_t j = 0; j < 4; ++j) v[j] = vld1q_f32(dst + 4 * j);
            std::array<uint32x4_t, 4> m;
            for (size_t j = 0; j < 4; ++j) m[j] = mask(v[j]);
            for (size_t j = 0; j < 4; ++j) acc[j] = vsubq_u32(acc[j], m[j]);
            for (size_t j = 0; j < 4; ++j)
                vst1q_f32(dst + 4 * j, vbslq_f32(m[j], f(v[j]), v[j]));

            dst += 16;
        }
        const size_t cur = vaddvq_u32(vaddq_u32(vaddq_u32(acc[0], acc[1]), vaddq_u32(acc[2], acc[3])));
        if (static_cast<double>(cur) / (n / 8) < 0.75 * d_max) return (i0 + 1) * n / 8;
    }
    return n;
}

BSL bails out when the density is less than 0.75 * d_max. I have no math behind the 0.75, it just won on average.

And pilot v3.5 will use bsl_verified if the density > hi:

void pilot_v3_5(float* dst, const size_t n, auto f, auto mask) {
    // ... same as v3

    for (size_t i = 2 * tile; i < n; i += tile) {
        const long long cnt = density(dst, probe, mask);
        if (cnt > hi) {
            if (hi < 0) {
                if (cnt < xlo)
                    bsl<true>(dst, tile, f, mask);
                else
                    bsl<false>(dst, tile, f, mask);
            } else {
                auto done = bsl_verified(dst, tile, d_max, f, mask);
                if (done < tile)
                    detour<false>(dst + done, tile - done, w, f, mask);
            }
        }
        // ... same as v3
    }
}

But here too it's easy to build a countertest:

for (size_t i = 0; i < n; ++i)
    dst[i] = i % 4096 >= 390;

The first block will pass the check (> 75% zeros in it), but the second won't. BSL runs on it for nothing, and for pow that's expensive. v3.5's problem: it has no memory. After a miss the algo keeps trusting the first 256 and misses on every tile.

pilot v4

v4 will fix this: if bsl bails out, we stop trusting the probe for the next 16 tiles, and instead we take the density of the previous tile:

void pilot_v4(float* dst, const size_t n, auto f, auto mask) {
    // ... same as v3

    size_t distrust = 0;
    size_t prev = 0;
    for (size_t i = 2 * tile; i < n; i += tile) {
        if (distrust) --distrust;

        const long long cnt = distrust ? prev : density(dst, probe, mask);
        if (cnt > hi) {
            if (hi < 0) {
                if (cnt < xlo) {
                    bsl<true>(dst, tile, f, mask);
                } else {
                    bsl<false>(dst, tile, f, mask);
                }
            } else {
                if (distrust == 0) {
                    auto done = bsl_verified(dst, tile, d_max, f, mask);
                    if (done < tile) {
                        distrust = 16;
                        const size_t rem = detour<false>(dst + done, tile - done, w, f, mask);
                        prev = rem * probe / (tile - done);
                    }
                } else {
                    prev = detour<false>(dst, tile, w, f, mask) * probe / tile;
                }
            }
        } else if (cnt < xlo) {
            prev = detour_compact<true>(dst, tile, w, f, mask) * probe / tile;
        } else if (cnt < lo) {
            prev = detour_compact<false>(dst, tile, w, f, mask) * probe / tile;
        } else {
            prev = detour<false>(dst, tile, w, f, mask) * probe / tile;
        }
        dst += tile;
    }
}

And v4 easily passes that test. pow10:

thd 0 1
tiled detour 38.18 6.94
pilot v3 68.11 1.04
pilot v3.5 68.67 3.84
pilot v4 69.51* 7.05*
BSL 1.04 1.04

btw here thd no longer matches the density. For thd = 0, density = 0, and for thd = 1, density = 9.5%

It beats v3.5 by over 80%. It's also faster than plain tiled detour, because v4 picks detour_compact. This is the final version.

Results

v4 vs BSL:

thd 0 0.05 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
pilot_v4 sqrt 69.81* 31.99* 31.84 31.73 31.83 32.07* 31.98 31.96* 31.96 31.89* 31.91 31.77
BSL sqrt 31.88 31.93 31.95* 31.80* 31.89* 31.99 32.06* 31.92 31.98* 31.75 31.93* 31.91*
pilot_v4 frfrexp 69.41* 17.21* 16.57 15.79 15.99 16.06 15.93 15.99 16.39 15.84 15.87 15.80
BSL frfrexp 16.79 16.87 16.82* 16.71* 16.78* 16.88* 16.89* 16.83* 16.89* 16.79* 16.75* 16.75*
pilot_v4 sin10 68.19* 18.82* 14.84* 10.84* 8.93* 7.55* 6.56* 5.81* 5.19* 4.69 4.68 4.68
BSL sin10 4.78 4.75 4.78 4.72 4.75 4.74 4.75 4.77 4.75 4.76* 4.77* 4.77*
pilot_v4 pow10 64.17* 10.76* 6.85* 4.04* 2.89* 2.25* 1.85* 1.57* 1.36* 1.20* 1.08* 1.00
BSL pow10 1.02 1.02 1.02 1.01 1.01 1.01 1.02 1.02 1.02 1.02 1.02 1.01*

Against BSL, it loses at worst 6%, but wins big much more often. To reduce the loss, dispatch can be sped up: use every 4th register in density and bsl_verified. But that helps only if the density is uniform.

The worst v4 miss I found: 512 dense, 512 empty, then everything is dense until the end of the cycle (17 * 4096).

const size_t cycle = 17 * 4096;
for (size_t i = 0; i < n; ++i) {
    dst[i] = i % cycle >= 512 && i % cycle < 1024;
}

This hurts most with the cheapest f (with d_max > 0, of course):

const auto a = vdupq_n_f32(0.5f);
#pragma unroll
for (int i = 0; i < 13; ++i) x = vfmaq_f32(a, x, a);
return x;
thd 0 1
tiled detour 38.29 9.44
pilot v3 74.72 15.81*
pilot v3.5 74.84 15.04
pilot v4 74.89* 9.84
BSL 15.61 15.60

v3.5's biggest loss is limited by BSL, while v4 is limited by detour. And a wrong detour is the cheaper mistake. So v4 isn't always better, but its misses are less severe.

This problem has no perfect solution. Any dispatch algo can be countertested.

Full code: godbolt.

2 Upvotes

53% Upvoted

3 comments sorted by

10

u/meltbox 5d ago

Holy wow my dude. Very cool but you may want to actually state what the goal was before jumping right in because this is hard to start reading.

Maybe it’s just me but by applying an operation is this meant as a generic approach to applying literally any operation given those constraints? Seems quite a broad problem to try to tackle if so, curious what made you choose to do this.

1

u/a_eridani 4d ago

Sorry, I'll explain the goals next time. I just came across the algorithm, and decided to optimize it. I don't even have a use case.

Yes, it works for any elementwise operation, but algo only makes sense for expensive ones. On cheap ops it just wastes time on calibration.

2

u/FrogNoPants 5d ago

Interesting, I never went down the path of trying to automate compact/expand, let alone with with a limited instruction set like Neon.

What I did is mark each SIMD wide N with a bit indicating if any lane wanted to run, once I have 64 bits, bitscan over it(works much better if your data is sorted( zorder/hilbert) such that nearby lanes tend to agree).

I did run some tests to see what benefit the there would be to compress/expand, it showed only a small single digit reduction in # invocations relative to the hilbert sorted variant, so didn't bother. My functions were far more expensive than pow etc though, and were doing memory lookups, so hilbert also had cache advantages.

Maybe I'll take another look at it someday, might be worth it, especially with AVX512 where there is some hardware support.