Will it be viable to use SIMD to further optimize the algorithm when the predicate is simple enough? For example, one can imagine that you take v[j] .. v[j+N-1] and somehow weed them into two partitions u[0] .. u[k] and u[k+1] .. u[N-1] (probably in SIMD registers), which get pasted accordingly at v[i] or v[j]. The current assembly indeed looks very tight but a higher throughput at the expense of longer code is a common theme in SIMD, I guess.

SIMD partitioning is definitely a thing, but I think it is ill-suited for Lomuto-style partitioning. I think fulcrum (what crumsort uses) or out-of-place partitioning (in both cases the output buffers for < and >= are distinct) like glidesort does is the most amenable to SIMDization.

Then partitioning is 'simply' a vector comparison, two masked compressing stores (through shuffles or _mm_mask_compressstoreu_epi32) with one of the masks inverted, and counting how many elements were smaller with _mm256_movemask_epi8 and a popcnt.

For an out-of-place partition you can interleave the loops of one going left-to-right and one going right-to-left to increase instruction-level parallelism.

You can port the algorithm from TFA pretty straightforwardly. For example with AVX2 and u64, load 4 elements from i and from j, compare the 4 elements swapped leftward with the threshold, use a lookup table of shuffles and compute shuffle(xs, lut[movemask(compare_result)]). store this value on the left and advance the pointer by popcnt(movemask(compare_result)). store the raw value previously loaded on the right and advance the pointer by 4.

You face some additional trouble when the right region of already-partitioned elements is smaller than 4. I think this is solvable but I don't have a good proposal off the top of my head.

Edit: well, one solution is to look at the first 8 elements. If at least 4 belong on the right, you can process these 8 carefully and then process the input forwards without worrying about your loads overlapping. If fewer than 4 belong on the right, you can swap these 8 with the last 8, process those 8 carefully, and process the input backwards without worrying about your loads overlapping.

All of this is what we have to do anyway when targeting AVX2, even if filtering some elements and discarding others or partitioning to 2 output buffers, because there are no compressing instructions.

That's actually fair enough. I'm not certain how well the CPU will like the overlapping SIMD reads/writes on the left hand side.

If you have AVX512 available shuffle(xs, lut[movemask(compare_result)]) could also be two compressions with one reverse shuffle in between.

I don't know either! I've never done something like this. Usually I am filtering to the same buffer or one or two other buffers, so I would never immediately load some stuff I had just stored.

Thanks for the article! I'm finding it hard to understand this part:

> That is, there is a single iterator scanning the array from left-to-right (jj). If the element is found to belong in the left partition, it is swapped with the first element of the right partition (tracked by ii), otherwise it is left where it is.

So you start left, the first element on the left partition belongs in the left partition, so you swap it into the right partition? And what about the element that was in the right partition, when do you check where that one belongs?

> So you start left, the first element on the left partition belongs in the left partition, so you swap it into the right partition?

Yes and no. You temporarily do swap it into the right partition, but then increment the pointers which redefines where the right partition is: you swap it with the first element of the right partition before incrementing both i and j. In essence what this does is move the entire right partition one step to the right, while putting the previously unknown element before it.

Consider this example:

    l l l l r r r r r ? ? ? ?
            ^         ^
            i         j

Now if the first ? was an r element, you could simply increment j and be done with it. But suppose it was an l element, then you have this scenario:

    l l l l r r r r r l ? ? ?
            ^         ^
            i         j

Note that "the first element of the right partition" is equivalent to "the first element after the left partition". Does it now make more sense that swapping the first element of the right partition and the unknown element (v[i], v[j]) is the right thing to do? After our swap we have this:

            +- swap --+
            v         v
    l l l l l r r r r r ? ? ?
            ^         ^
            i         j

So now incrementing i and j both fixes our invariant:

    l l l l l r r r r r ? ? ?
              ^         ^
              i         j

> And what about the element that was in the right partition, when do you check where that one belongs?

That one still belongs in the right partition, which is why we increment both i and j.

Ahh, thank you, this clarifies things! I was thinking of Quicksort, where j moves left until it meets i, and missed the fact that they both move in the same direction and i runs through the entire list. Thanks!

EDIT: I'm now on the desktop and see that your graphs clarify this much better, on mobile I had Dark Reader and they were invisible. You may want a white background.

> I was thinking of Quicksort, where j moves left until it meets i, and missed the fact that they both move in the same direction and i runs through the entire list.

You were thinking of Hoare partitioning. Quicksort can use either Hoare or Lomuto (or any other) partitioning algorithm.

Yes, now that you mention it, that does make sense, thank you.

Is the assembly the inner loop or the whole function? The function looks like it could be restructured into a tail call - the trailing swap is quite like the loop body - but I'd be surprised if a compiler managed that transform.

The assembly is only for the inner loop, I did not include the preamble/postamble of the functions.

You can look at the full assembly here: https://cpp.godbolt.org/z/zzzTh47PG. But the full assembly isn't a fair comparison right now because I did not bother to convert both functions to have the same signature (Andrey's version assumes the pivot is in the array and selects it in the function).

Virtually all the time for any non-trivial input is spent in the inner loop though.