This is one of the best uses I've found for Singeli[0]. Here's how I implemented an AVX2 transpose kernel similar to that in transpose_Vec256_kernel, for generic type and vector/kernel size (I use unpack instructions rather than shuffle and blend, which I think is probably faster since it's just one instruction for each interaction):
The language is oriented towards compile-time array programming instead of managing a bunch of individual vectors. So you have runtime vec_select{} (docs at [1]), mirrored by compile-time select{}, and the indices generated by pairs{} can be used in either.
It's possible to get O(N^(3/2)/log N) in an ordinary interpreter with some changes to the code, assuming linear-time ~ with hashing. The idea is to leave the primes in ⍵ and stop when it stops changing, which will happen at the first prime over √N. To get the complexity multiply by O(N) for each step. It's also a small change to get the multiples to remove to start at n×n instead of n.
i←¯1⋄{⍵~n×n↓⍳1+⌊N÷n←i⊃⍵⊣i+←1}⍣≡1↓1+⍳N
I think this is about as good as can be done with ~ instead of marking out bits. And I wouldn't say it's as easy as the imperative version!
That's not O(N log log N), it's more like N^2. Prime sieves are hard to implement well with immutable arrays for obvious reasons; there are some cool methods but they're definitely harder. I'm ashamed to be part of a community that won't cop to this.
The algorithm iterates over numbers ⍺ from 2 to N, removing the multiples that are greater than ⍺ and no greater than N from p. If the removal with ~ has to inspect all of p, then all the primes are there so we have asymptotically at least N/log N entries by the prime number theorem and we get N^2/log N time (when ⍺ is over N/2, no multiples are in range so this can be skipped, but that just cuts a constant factor 2 from the time). Conceivably p could be marked sorted, so the entries to remove could be found with a binary search. This is a bit harder to analyze, but I think each prime under √N will cause the list to change, and incur N/log N data movement. So you get at least (N/log N)^(3/2) cost, still quite a bit worse than linear.
Edit: changed the algorithm while I was writing... the new one is better, it keeps one list p of primes and one list ω of not-yet-marked-out numbers. However, primes are removed from ω one at a time, so that each of the N/log N primes has to be moved for each one before it, giving at least (N/log N)^2 cost (I mean, maybe an interpreter could binary search and also recognize when ~ only drops the first entry and do that by slicing? But the (N/log N)^(3/2) from above definitely holds). Mutating a bit array in place is pretty important to classical sieve performance.
It's very strange to see handwriting lumped in with typewriting, to be described as limited relative to screens! Iverson notation was a 2D format (both in handwriting and typeset publications) making use of superscripts, subscripts, and vertical stacking like mathematics. It was linearized to allow for computer execution, but the designers described this as making the language more general rather than less:
> The practical objective of linearizing the typography also led to increased uniformity and generality. It led to the present bracketed form of indexing, which removes the rank limitation on arrays imposed by use of superscripts and subscripts.
I think this is more true than they realized at that time. The paper describes the outer product, which in Iverson notation was written as a function with a superscript ∘ and in APL became ∘. followed by the function. In both cases only primitive functions were allowed, that is, single glyphs. However, APL's notation easily extends to any function used in an outer product, no matter how long. But Iverson notation would have you write it in the lower half of the line, which would quickly start to look bad.
I don't think this really describes neon_prefixsum_fast as a whole? The algorithm does use a Hillis-Steele sum on sums of 4 values, but each of these is computed with a sequential sum, interleaving those with a transposed order. In terms of what's added to what, it's actually quite a bit like my "Sequential broadcasting" picture from [0]. The reference I'd use for a general form is "Parallel Scan as a Multidimensional Array Problem"[1], breaking 16 elements into a 4x4 array; the paper describes how the scan splits into a row-wise scan, plus values obtained from an exclusive scan on carries from the rows.
More on not being able to find π, as I'm piecing it together: given only the field structure, you can't construct an equation identifying π or even narrowing it down, because if π is the only free variable then it will work out to finding roots of a polynomial (you only have field operations!) and π is transcendental so that polynomial can only be 0 (if you're allowed to use not-equals instead of equals, of course you can specify that π isn't in various sets of algebraic numbers). With other free variables, because the field's algebraically closed, you can fix π to whatever transcendental you like and still solve for the remaining variables. So it's something like, the rationals plus a continuum's worth of arbitrary field extensions? Not terribly surprising that all instances of this are isomorphic as fields but it's starting to feel about as useful as claiming the real numbers are "up to set isomorphism, the unique set whose cardinality matches the power set of the natural numbers", like, of course it's got automorphisms, you didn't finish defining it.
I was interested in how it would make sense to define complex numbers without fixing the reals, but I'm not terribly convinced by the method here. It seemed kind of suspect that you'd reduce the complex numbers purely to its field properties of addition and multiplication when these aren't enough to get from the rationals to the reals (some limit-like construction is needed; the article uses Dedekind cuts later on). Anyway, the "algebraic conception" is defined as "up to isomorphism, the unique algebraically closed field of characteristic zero and size continuum", that is, you just declare it has the same size as the reals. And of course now you have no way to tell where π is, since it has no algebraic relation to the distinguished numbers 0 and 1. If I'm reading right, this can be done with any uncountable cardinality with uniqueness up to isomorphism. It's interesting that algebraic closure is enough to get you this far, but with the arbitrary choice of cardinality and all these "wild automorphisms", doesn't this construction just seem... defective?
It feels a bit like the article's trying to extend some legitimate debate about whether fixing i versus -i is natural to push this other definition as an equal contender, but there's hardly any support offered. I expect the last-place 28% poll showing, if it does reflect serious mathematicians at all, is those who treat the topological structure as a given or didn't think much about the implications of leaving it out.
Ooh, I've run into this one before! I'm a big fan of interval index[0], which performs a binary search, so Josh's suggestion is the one I prefer as well (the implementation might sometimes optimize by transforming it into a table lookup like the other solutions). Searching for +`≠¨ in my BQN files turns up a few uses, including one used to group primitives into types in a utility involved in compiling BQN[1] and an annotated function used a few times in the markdown processor that builds BQN's website[2].
"Ken Iverson - The dictionary of J contains an introductory comment that J is a dialect of APL, so in a sense the whole debate is Ken's fault! He is flattered to think that he has actually created a new language."
I think the article expresses no position. Most source code for array languages is not, in fact, inspired by APL. I encourage you to check a few random entries at [0]; Kap and April are some particularly wordy implementations, and even A+ mostly consists of code by programmers other than Whitney, with a variety of styles.
I do agree that Whitney was inspired to some extent by APL conventions (not exclusively; he was quite a Lisp fan and that's the source of his indentation style when he writes multi-line functions, e.g. in [1]). The original comment was not just a summary of this claim but more like an elaboration, and began with the much stronger statement "The way to understand Arthur Whitney's C code is to first learn APL", which I moderately disagree with.
Dunno why electroly is dragging me into this but I believe you've misread the article. When it says "His languages take significantly after APL" it means the languages themselves and not their implementations.
It looks like a weirdo C convention to APLers too though. Whitney writes K that way, but single-line functions in particular aren't used a lot in production APL, and weren't even possible before dfns were introduced (the classic "tradfn" always starts with a header line). All the stuff like macros with implicit variable names, type punning, and ternary operators just doesn't exist in APL. And what APL's actually about, arithmetic and other primives that act on whole immutable arrays, is not part of the style at all!
I just want you to accurately represent the performance, which as I've noted is not a big change. If a 2x speedup is definitive, why are you so desperate to push the number higher? You give the speedup of vqsort over std::sort in ideal conditions (AVX-512) but when asked to take into account faster scalar sorts, you go to pdqsort's 2x factor (again this implies that it's the best available!) rather than fluxsort's 3x or ipnsort's 4x.
No, I don't believe your current AVX2 performance will always be better than scalar algorithms. Lukas's benchmarks show... what, a 30% difference between vqsort and ipn in this case? I have 2-3x for random numbers (yes, other cases are not always this good) in Singeli sort, which I'm choosing not to promote until it's more mature. No, you don't need to range check in partitioning. If you don't understand how this can be achieved, that's your problem.
Yes, only your "10x" and "5-10x" numbers are overstated. Which is to say, all the quantitative comparison you've shared here, or in the vqsort README (aside from the link to Lukas's measurements, which is linked with no indication it benchmarks against faster sorts). What you share is defensible in that it's technically correct, but it's not helpful.
As an example, consider an application that uses heapsort now, but ships a unified binary for x86. The programmers could easily switch to pdqsort, which is well-tested and available in many languages. If they read only material you publish, they may get the impression that SIMD support (at least SSE4) is required to get a substantial improvement, and stick with heapsort to avoid the need to ship multiple binaries or add architecture detection. All you need to correct this is one mention that you can do better than std::sort! By leaving this out, you instead reinforce the idea that it's a reasonable choice if you have to go scalar.
ipnsort has no explicit vectorization; some parts may be auto-vectorized. From a quick glance through perf results, I don't see any vector instructions used in the parts that take up time. I didn't see any anywhere, actually. Same story with other quicksorts, although I know fluxsort/crumsort's initial analyzer is designed to be auto-vectorized.
I don't get why you dismiss radix sort so quickly, when your own research shows that an MSD/LSD hybrid addresses the bandwidth problem well? Thanks for sharing that paper by the way, very elegant approach. The bandwidth problem is not one I've seen discussed although I see how it would arise. I just don't deal with big arrays that much—and it seems if you're going to market vqsort as general-purpose you should also care about this use case! However, a nice thing about radix sort is that it can be used as a base case for quicksort/samplesort, and any partitioning done so far allows fewer steps to be used. I have a case in SingeliSort to do 4-byte sorting that fits in 2-byte range nearly in-place[0]. It's only used on 2^16 elements or fewer to avoid possible cache associativity problems, so that should fit in L2 and not use system bandwidth.
Comparison-based algorithms I consider most relevant at the moment are fluxsort/crumsort/quadsort (all by the same author; they share many pieces), ipnsort (largely derived from those plus pdqsort but it may be easier to read), and glidesort somewhat (the gliding technique is cool but I think there are flaws to address, see [1]). For distribution sorting, there's radix, the version of counting sort where you reconstruct from counts instead of moving the data (output can be vectorized with a scan for large ranges!), and I've developed Robin Hood sort to take advantage of smooth-ish distributions such as uniform random. I keep notes at [2], which are pretty out of date. I should be writing there instead of here...
I of course don't think there's any requirement to know all this to contribute to sorting research. I take issue with making blanket claims like "I'd consider AVX2 to be table stakes for any sorting algorithm" and "anything as energy-inefficient as scalar code running on big OoO cores is a hard sell" without being familiar with the state of the art in scalar sorting.
Apologies, I misread that section and thought the only difference between tables 1 and 2 was the array size! Thanks for correcting my impression here, that's good to know. And very sorry for the misrepresentation.
Oh, I didn't know that about multi-way merge (regrettably I haven't really wrapped my head around SIMD merging yet). Got a source for this? https://vldb.org/pvldb/vol8/p1274-inoue.pdf mentions it but it looks more like interleaving several binary merges than true multi-way. Although I guess from a memory perspective there's no difference.
Vector registers are great, but it's important to remember that random-access cache is also very powerful hardware, and can be better-suited to a lot of searching and sorting tasks.
One task where, unlike sorting, I think vectors have no shot, is random shuffling. For small sizes (okay, but too large to fit all in registers) it's hard to see how anything could beat Knuth shuffle since it just does one swap per value. At large sizes it's not cache-friendly, and there are known methods based on quicksort and mergesort that work better. But I found that a radix pass beats these easily, and I'd expect this to hold even with AVX-512 versions of merge/partition, because it takes 8 steps of those to add up to one 1-byte radix step. Sorting runs into similar problems at very large sizes. The vqsort authors found that at 1e8 elements a hybrid of ips4o (samplesort, think many-way partition) beats pure vqsort.
Please, you are overstating your case in a way that distracts from the issues at hand. I do agree that Linus's opinions are irrelevant to a lot of computing, and that sorting must make some use of SIMD for the best performance, and that there's good work in vqsort. But the 10x factor is highly misleading. Sure, it's technically correct that you "can be about 10x as fast" as a scalar sort. In the same sense that a scalar sort can be 10x as fast as an AVX-512 one.
The 10x number is from comparing against std::sort on random data, yes? Both common std::sort implementations are branching quicksorts from the era when introsort was still considered exciting new technology, and are much slower than pdqsort and other branchless quicksorts. The measurements at [0] show by my estimate about a gain of 3x against pdqsort and 2x against ipnsort. For vqsort only; intel barely beats ipnsort. I also pointed out to you at [1] benchmarks showing 1450MB/s for 32-bit radix sort at a million elements on M1, which is comparable to your Xeon measurements and nearly 3x faster than the 499MB/s for M1 in the vqsort paper.
If you're going to quote a std::sort number, at least give the name instead of implying it's a good representative for scalar sorts! I would strongly recommend, however, that you benchmark against and study more modern scalar algorithms. I know SIMD and have even vectorized parts of pdqsort in the past; I've chosen to work with scalar algorithms for the past year or two because I think these methods still have a lot to offer.
> This is true for all languages, not just Python.
Not quite: APL does comparisons with a configurable relative tolerance, which is nonzero by default. My article [0] about comparing one number to many quickly in this system starts with a discussion of the reasoning and mathematics there. It brings a lot of implementation difficulty, particularly with hash tables. It's possible to build a hash on tolerant doubles by making two hashes per lookup, and even complex numbers with four per lookup, but no one's figured out how to deal with entries that contain multiple numbers yet.
Nonetheless I would say it does make it easier to write working programs overall. But not by much: I've been working on an APL derivative and initially left out comparison tolerance as there were some specifics I was unsure of. A year or so later I noticed I hardly ever needed it even for numerical work and dropped any plans to add it to the language.
https://github.com/dzaima/CBQN/blob/v0.11.0/src/singeli/src/...
The language is oriented towards compile-time array programming instead of managing a bunch of individual vectors. So you have runtime vec_select{} (docs at [1]), mirrored by compile-time select{}, and the indices generated by pairs{} can be used in either.
[0] https://github.com/mlochbaum/Singeli/
[1] https://github.com/mlochbaum/Singeli/tree/master/include#sim...