Asking a bit-permutation what it does (and cloning it)

Given some black-box permutation that operates on an array of integers, it's easy to find out what the permutation does: feed it an array containing the number 0 .. n-1, apply the permutation to it, then the result tells you exactly what the permutation does. And now that you have that resulting array that describes the permutation, you no long need the black-box permutation - you can apply the permutation by using that array. But what about bit-permutations? You can't exactly feed them an array containing 0 .. n-1 .. or can you?

What you can do is feed in an integer constructed out of all the least significant bits of that array, permute it, and deconstruct it into an array. Repeat for every relevant bit, performing log₂(n) applications of the black-box bit-permutation. Both taking bit-slices of an array like that and deconstructing the results back into a "normal" array correspond to some form of bit-matrix transpose, so we won't be messing around with individual bits one at a time. For example, for a 64-bit bit-permutation, for which I derived the relevant transposes in Permuting bits with GF2P8AFFINEQB (or you could have them automatically generated using this little tool):

__m512i Transpose8x64(__m512i x) {
    __m512i s1 = _mm512_set_epi8(
        7, 15, 23, 31, 39, 47, 55, 63,
        6, 14, 22, 30, 38, 46, 54, 62,
        5, 13, 21, 29, 37, 45, 53, 61,
        4, 12, 20, 28, 36, 44, 52, 60,
        3, 11, 19, 27, 35, 43, 51, 59,
        2, 10, 18, 26, 34, 42, 50, 58,
        1, 9, 17, 25, 33, 41, 49, 57,
        0, 8, 16, 24, 32, 40, 48, 56);
    __m512i m2 = _mm512_set1_epi64(0x8040201008040201);
    x = _mm512_permutexvar_epi8(s1, x);
    x = _mm512_gf2p8affine_epi64_epi8(m2, x, 0);
    return x;
}

std::array<uint8_t, 64> analyze_permutation(std::function<__m512i(__m512i)> perm)
{
    __m512i input = _mm512_setr_epi64(
        0xAAAAAAAA'AAAAAAAA,
        0xCCCCCCCC'CCCCCCCC,
        0xF0F0F0F0'F0F0F0F0,
        0xFF00FF00'FF00FF00,
        0xFFFF0000'FFFF0000,
        0xFFFFFFFF'00000000,
        0x00000000'00000000,
        0x00000000'00000000);
    __m512i permuted = perm(input);
    __m512i indices = Transpose8x64(permuted);
    std::array<uint8_t, 64> res;
    _mm512_storeu_epi64(&res[0], indices);
    return res;
}

If your bit-permutation function doesn't deal in __m512i but in plain old integers, you can of course call it 6 times. The inputs to the permutation can alternatively be calculated fairly easily by their index, but I'm going to show it in C#:

ulong mask_i = ~(ulong.MaxValue / ((1UL << (1 << i)) + 1));

Using the indices

A nice way to bit-permute one u64 at the time given an array of indices, is the vpshufbitqmb approach mentioned in Dynamic bit shuffle using AVX-512:

// source: https://lemire.me/blog/2023/06/29/dynamic-bit-shuffle-using-avx-512/
uint64_t faster_bit_shuffle(uint64_t w, uint8_t indexes[64]) {
  __m512i as_vec_register = _mm512_set1_epi64(w);
  __mmask64 as_mask = _mm512_bitshuffle_epi64_mask(as_vec_register,
     _mm512_loadu_si512(indexes));
  return _cvtmask64_u64(as_mask);
}

Or, using even more bit-matrix transposes, we can leverage vpermb to do the bit-permutation for us across 8 u64 at the same time. vpermb does not permute at the bit-level, obviously, but it can permute 64 u8 and if we transpose our 8 u64 into 64 u8 (one u8 at index i holding the 8 i'th bits of the u64s) then we can vpermb those bytes and after transposing back into 8 u64 we will have bit-permuted them according to the specified indices. But we can do a little better than that. The second transpose, going back from 64 u8 to 8 u64, starts with a vpshufb with a fixed shufflemask. We're probably going to permute bits with the same indices multiple times (otherwise why bother doing any of this), that shuffle can be shaved off of the actual bit-permutation function by permuting our indices beforehand:

// do this once, before using the indices
__m512i byterev_qwords = _mm512_set_epi8(
    8, 9, 10, 11, 12, 13, 14, 15,
    0, 1, 2, 3, 4, 5, 6, 7,
    8, 9, 10, 11, 12, 13, 14, 15,
    0, 1, 2, 3, 4, 5, 6, 7,
    8, 9, 10, 11, 12, 13, 14, 15,
    0, 1, 2, 3, 4, 5, 6, 7,
    8, 9, 10, 11, 12, 13, 14, 15,
    0, 1, 2, 3, 4, 5, 6, 7);
__m512i modified_indices = _mm512_shuffle_epi8(indices, byterev_qwords);

__m512i permute_bits_qword(__m512i modified_indices, __m512i x)
{
    __m512i s1 = _mm512_set_epi8(
        7, 15, 23, 31, 39, 47, 55, 63,
        6, 14, 22, 30, 38, 46, 54, 62,
        5, 13, 21, 29, 37, 45, 53, 61,
        4, 12, 20, 28, 36, 44, 52, 60,
        3, 11, 19, 27, 35, 43, 51, 59,
        2, 10, 18, 26, 34, 42, 50, 58,
        1, 9, 17, 25, 33, 41, 49, 57,
        0, 8, 16, 24, 32, 40, 48, 56);
    __m512i m2 = _mm512_set1_epi64(0x8040201008040201);
    __m512i s3 = _mm512_set_epi8(
        63, 55, 47, 39, 31, 23, 15, 7,
        62, 54, 46, 38, 30, 22, 14, 6,
        61, 53, 45, 37, 29, 21, 13, 5,
        60, 52, 44, 36, 28, 20, 12, 4,
        59, 51, 43, 35, 27, 19, 11, 3,
        58, 50, 42, 34, 26, 18, 10, 2,
        57, 49, 41, 33, 25, 17, 9, 1,
        56, 48, 40, 32, 24, 16, 8, 0);
    x = _mm512_permutexvar_epi8(s1, x);
    x = _mm512_gf2p8affine_epi64_epi8(m2, x, 0);
    x = _mm512_permutexvar_epi8(modified_indices, x);
    x = _mm512_gf2p8affine_epi64_epi8(m2, x, 0);
    x = _mm512_permutexvar_epi8(s3, x);
    return x;
}

Back to those permuted masks

Let's back up to when we had the variable permuted, before transposing it into a nice neat array of indices. If we're dealing with a BPC permutation, a permutation that can be described by how the bits of the indices (not the original bits that are being permuted) are permuted and optionally also complemented, then that is both easy to detect from the entries of permuted, and easy to analyse. Each entry of permuted will correspond to either one of the original query-masks that we put into the bit-permutation, or its complement, with each query mask or its inverse appearing exactly once (and not both). If (and only if) the permutation is of that form, then it is a BPC permutation, and the corresponding signed permutation is a convenient way to describe it.

My little AVX512 BPC code generator also takes that representation (generalized to 9 index bits, not 6) as its input, as an alternative to a giant array of raw indices, and does everything it does in terms of such a 9-element signed permutations (that does not prove that that's the best thing to do, maybe I should be working with signed permutation matrices that seems less error-prone). This representation has also been useful in experiments to optimize larger (multi-vector) transposes (any reasonable transpose is a BPC permutation, BP actually, sans C - but vgf2p8affineqb forces me to care about index-bit-complements), which is still a work in progress.

---

As usual, LLMs were not intentionally consulted while writing this post, or the accompanying code. Some of the code was generated by old-fashioned non-LLM tools. LLM-generated "AI overviews" have occasionally popped up on Google, for example while looking into the hyperoctahedral group.

Fenwick trees for products mod 2ⁿ

After browsing some articles about Fenwick trees (aka Binary Indexed Trees) and seeing remarks such as "let f be some group operation" and "prefix operations such as sum, product, XOR and OR", I started to wonder about the case of combining products (mod 2ⁿ) with range queries. In this case f is not a group operation, unless we restrict ourselves to working with odd values, which are invertible modulo a power of two (this can be done fairly efficiently, see eg Hacker's Delight or a more recent improved algorithm). And neither do we restrict ourselves to computing prefix products, which would have worked without requiring inverses.

Inverses of elements are not strictly required for range queries, and are not used in the usual implementations of a Fenwick tree: while a subtraction can be viewed as (or defined as) adding a negative, that already shows that what we need is an inverse operation, not necessarily an inverse element. Inverse elements are commonly used in the range-update/point-query and range-update/range-update configurations of Fenwick trees, but can be avoided there at the cost of some code duplication. So it would be possible to have a Fenwick tree where f is multiplication in ℤ\{0}, which is not a group operation either, since we can use division instead of multiplication by the reciprocal. Doing that and taking the result mod 2ⁿ would be a possible implementation of "a Fenwick tree for products mod 2ⁿ". As a bonus we could represent zero as 2ⁿ, which maps to zero in the end but doesn't break division. Since the size of the product of two integers is approximately the sum of the sizes of the multiplicands, this would tend to build up large integers of a size comparable to the size (total size, not number of elements) of the entire tree. Not great.

A more reasonable solution would be to convert numbers x (for non-zero x) to the form d*2^k where k = tzcnt(x), tracking the exponents in a normal additive Fenwick tree and the odd integers d in a multiplicative (mod 2ⁿ) Fenwick tree which Just Works^(tm) because odd integers are invertible mod 2ⁿ. Zero can be represented as 2ⁿ again. Point-updates can be statelessly and commutatively undone^[1], as long as you promise not to do it wrong in a way that creates negative exponents.

Something like that is probably what you expected based on the title of the post anyway.

---

I may have waffled a bit more than usual, especially in the second paragraph, but I did it by hand. No LLMs were used to write this post. I accidentally triggered Google AI summaries a couple of times while looking things up though.

[1] For general monoids, you could save a backup of the entries that are modified by a point-update, then undo the update by restoring from the backup - stateful and not commutative. That's a bit annoying, but not quite impossible, as is sometimes suggested.

What is abs(x - y)

What is abs(x - y), why do I even bother asking a question like that, much less answer it? Trust me, it's more interesting than it looks.

To start with, here's something that it is not. Despite appearances, it's not the absolute difference between x and y. If we were working with unlimited integers, sure, but as usual I'm imagining a context in which every integer is implicitly assumed to have some finite size. There are several valid ways to compute the absolute difference between x and y (henceforth absdiff_u(x, y) for the unsigned absolute difference, and absdiff_s(x, y) for the signed absolute difference, which does not return a signed result but rather interprets its inputs as signed integers) but abs(x - y) is not one of them when you're working with fixed-size integers.

Absolute difference

Probably the most obvious way to compute absdiff_γ(x, y) is max_γ(x, y) - min_γ(x, y), where γ is the signedness of choice. You can think about geometrically (in terms of a number line) if you want. There's no "funny business", the subtraction doesn't wrap in the unsigned case and doesn't overflow in the signed case.

Based on the ranges of the functions alone there's an argument to show that abs(x - y) cannot possibly implement absdiff of either signedness: abs(x - y) cannot return a value greater than half the range of its type. Eg no value greater than 0x80000000 when working with 32-bit integers [1]. The distance between two integers can be up to and including the maximum value that can fit in the corresponding unsigned type, almost half of those potential distances cannot be returned by abs(x - y).

By the way an interesting and sometimes useful way to compute absdiff_u(x, y) is subus(x, y) | subus(y, x) where subus is subtraction with unsigned saturation. At least one of subus(x, y) and/or subus(y, x) is zero (if x == y then both are zero, otherwise one of them saturates to zero) so the parts can be ORed or added or XORed to give whichever of them is not zero, or zero if both are zero.

Clockwise distance

Another possible notion of distance between fixed-size integers is their clockwise distance: how far do you have to go in the positive direction around the number ring (the way I drew the ring back then, actually it would be counter-clockwise distance, but let's not dwell on that) to get from x to y. If, having started at x, we reach y before crossing the "boundary" (there is no boundary, but it feels like one) between the highest unsigned integer and zero, then that distance is obviously y - x. Otherwise, it's the distance to the highest unsigned integer, plus one, plus the distance from zero to y:

0xFFFFFFFF - x + 1 + y =
~x + 1 + y =
-x + y =
y - x

So the clockwise distance from x to y is y - x either way, that's nice and simple.

Since we're working around the number ring, there is no difference between a clockwise distance between signed integer and a clockwise distance between unsigned integers. You could even interpret the distance as a signed quantity, reinterpreting very large positive distances into smaller negative distances. [2]

This notion of difference is obviously not commutative, subtraction being famously anti-commutative. We should expect that, because depending on which way around x and y are, the clockwiseness either lets you go the short way between x and y or forces you to go the long way around the ring. For antipodal points on the ring it doesn't make any difference which way around you go, but for most pairs of x and y there is a short way and a long way. That leads into the next notion of distance.

The shortest way around the ring

What if we forget about going clockwise and could go either way around the ring whichever is shorter, we get a commutative notion of distance again. Based on how the clockwise distance was calculated, one way to express the "shortest way around the ring" is min_u(x - y, y - x). The notion of minimum used here should be the unsigned minimum, even if x and y were thought of as signed (which once again doesn't matter around the ring), the signed minimum would pick a large negative distance over a short positive one which is not what we set out to compute.

min_u(x - y, y - x) is equivalent to (... insert drumroll ...) abs(x - y). Let's prove it, why not. I will assume 32-bit integers, that's not fundamental to the proof but it's easier to assume a specific size.

First let's recall that y - x = -(x - y).

• If x - y = 0x80000000 then y - x is also 0x80000000 and we have min_u(0x80000000, 0x80000000) = 0x80000000 = abs(0x80000000).
• If x - y < 0x80000000 (ie non-negative when interpreted as a signed integer) then y - x would be > 0x80000000 (negative when interpreted as a signed integer), x - y < y - x and abs doesn't perform its conditional negation, min_u(x - y, y - x) = x - y = abs(x - y).
• If x - y > 0x80000000 (ie negative when interpreted as a signed integer) then y - x would be < 0x80000000 (non-negative when interpreted as a signed integer), x - y > y - x and abs does perform its conditional negation, min_u(x - y, y - x) = y - x = abs(x - y).

There is nothing inherently wrong with this notion of distance, but likely most uses of abs(x - y) "in the wild" were not written with the intent to compute the "shortest way around the ring" type of distance.

Programming-language Annoyances

One aspect that I've pointedly ignored in the rest of this post is concerns such as trapping on overflow (either in normal arithmetic or in abs), or undefined behaviour, or the difference of unsigned integers yielding a signed integer of a wider type. Such concerns are your responsibility if you use any of the expressions from this post; I mention them primarily to head off pointless discussions.

---

I didn't write this post with AI. You can tell because the writing is dense and awkward as usual.

1. [1] but 0x80000000 is a potential output, so hopefully you are interpreting the result of abs(x - y) as an unsigned integer, or you don't mind the occasional Integer.MIN_VALUE showing up.
2. [2] the math doesn't care about your interpretation, but you might.