mlochbaum

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):

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...

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

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.

(https://www.jsoftware.com/papers/APLDesign.htm)

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.

[0] https://mlochbaum.github.io/BQN/implementation/primitive/fol...

[1] https://ashinkarov.github.io/pubs/2022-scan.html

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].

[0] https://aplwiki.com/wiki/Interval_Index

[1] https://github.com/mlochbaum/BQN/blob/717555b0db/src/pr.bqn#...

[2] https://github.com/mlochbaum/BQN/blob/717555b0db/md.bqn#L45-...

It was the subject of quite some debate, see "Panel: Is J a Dialect of APL?" at http://www.jsoftware.com/papers/Vector_8_2_BarmanCamacho.pdf . Ken and Roger backed off this stance after witnessing the controversy.

"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.

[0] https://aplwiki.com/wiki/List_of_open-source_array_languages

[1] https://code.jsoftware.com/wiki/Essays/Incunabulum

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!

Ordinarily I'd make fun of the Germans for giving such an ugly name to a nice concept, but I've always found "comfortable" to be rather unpleasant too (the root "comfort" is fine).

It's just. So gross. Say it. Sudden interruption of slime coming up your throat. Like walking out the door into a spiderweb. Alphabetically I was mistaken but in every way that matters I was right.

Well, do you know how it works? Don't judge a book by its cover and all. Although none of these are entirely aiming for elegance. The first is code golf and the other two have some performance hacks that I doubt are even good any more, but replacing ∧≢⥊ with ∧⌜ in the last gets you something decent (personally I'm more in the "utilitarian code is never art" camp, but I'd have no reason to direct that at any specific language).

The double-struck characters have disappeared from the second and third lines creating a fun puzzle. Original post https://www.ashermancinelli.com/csblog/2022-5-2-BQN-reflecti... has the answers.

The point that the article is addressing (but you have to ignore the image and study the equations to see this!) is that this sort of shifting can't equalize everything. In the span of 3 white keys C to E at the front, you have 2 black keys at the back, so if you take r to be the ratio of back-width to white key front-width then you have 3 = 5r. But in the 4 keys F to B, you've got 3 black keys so 4 = 7r. No single ratio works! So the article investigates various compromises. The B/12 solution is what seems to me the most straightforward, divide white keys in each of the sections C to E and F to B equally at the back, and don't expect anyone to notice the difference.

And the reason +˝ is fairly fast for long rows, despite that page claiming no optimizations, is that ˝ is defined to split its argument into cells, e.g. rows of a matrix, and apply + with those as arguments. So + is able to apply its ordinary vectorization, while it can't in some other situations where it's applied element-wise. This still doesn't make great use of cache and I do have some special code working for floats that does much better with a tiling pattern, but I wanted to improve +˝ for integers along with it and haven't finished those (widening on overflow is complicated).

The relevant operations for matrix multiply are leading-axis extension, shown near the end of [0], and Insert +˝ shown in [1]. Both for floats; the leading-axis operation is × but it's the same speed as + with floating-point SIMD. We don't handle these all that well, with needless copying in × and a lot of per-row overhead in +˝, but of course it's way better than scalar evaluation.

[0] https://mlochbaum.github.io/bencharray/pages/arith.html

[1] https://mlochbaum.github.io/bencharray/pages/fold.html

To be clear, you are referring to the preface to "An Introduction to Array Programming in Klong", right? Having just checked it, I find this to be a very strange angle of attack, because that section is almost exclusively about why the syntax in particular is important. Obviously you disagree (I also think the syntax is overblown, and wish more writing focused on APL's semantic advantages over other array-oriented languages). I think this is a simple difference in taste and there's no need to reach so far for another explanation.

Oddly enough, the biggest mistake in how I presented BQN early on was thinking only APL insiders would be interested, when in fact the APLers went back to APL and people who hadn't tried other array languages or hadn't gotten far with them were were most successful with BQN. Plenty of people coming to BQN have worked with Numpy or whatever, but I don't think this has the same deterrent effect; they see BQN as different enough to be worth learning. Julia in particular is very different: I don't find that it culturally emphasizes array programming at all.

Search for "teaching" at https://aplwiki.com/wiki/APL_conference. I count at least five papers about teaching non-APL topics using APL. The language is not only possible to read, it's designed for it.

"already": APL dates back to about 1966, and even K from 1993 predates Numpy and Julia. But yes, we do not live in caves and are familiar with these languages. Klong has even been implemented in Numpy, see https://github.com/briangu/klongpy.

Author of BQN here, I agree with how section "What about APL?" describes the APL family as not fundamentally better (although details like indexing are often less messy). I outlined a system with lexically-scoped named axes at https://gist.github.com/mlochbaum/401e379ff09d422e2761e16fed... . The linear algebra example would end up something like this:

    solve(X[i,_], Y[j,_], A[i,j,_,_]) = over i, j/+: Y * linalg_solve(A, X)

The queue method is popular, but there's a much faster (branch-free) and in my opinion simpler way, known as the van Herk/Gil-Werman algorithm in image processing. It splits the input into windows and pairs up a backward scan on one window with a forward scan on the next. This works for any associative function. I was very surprised when I learned about it that it's not taught more often (the name's not doing it any favors)! And I wrote a tutorial page on it for my SIMD-oriented language, mostly about vectorizing it which I didn't quite finish writing up, but with what I think is a reasonable presentation in the first part: https://github.com/mlochbaum/Singeli/blob/master/doc/minfilt...

I also found an interesting streaming version here recently: https://signalsmith-audio.co.uk/writing/2022/constant-time-p...

EDIT: On closer inspection, this method is equivalent to the one I described, and not the one I'm used to seeing with queues (that starts my tutorial). The stack-reversing step is what forms a backwards scan. The combination of turning it sequential by taking in one element at a time but then expressing this in functional programming makes for a complicated presentation, I think.

I'd forgotten about patience sorting! Not that I was ever exactly familiar with it. JesseSort seems nearly identical, except that it allows insertion to either side of the lists, tries the previous index before doing a binary search, and uses pairwise instead of k-way merge at the end. The expected Θ(√n) list count is the same as well. An obvious advantage of going double-ended is that you get automatic O(n) sorting if the input's an ascending or descending run, instead of just one of those.

I think both algorithms can be made stable; if there's a barrier for patience sort it would be that the merge phase pops off the end of the sorted lists so it's going in the opposite direction as insertion? For JesseSort, an element that's equal to a previous one can only be inserted to the same list, or a more central one (created later). So you need to merge sublists giving priority to the outer ones, and also never append an element to the bottom of a sublist if it's equal to that list's current minimum. Doesn't look like JesseSort as implemented in jessesort_c.pyx does either of these: it merges non-adjacent lists and has a backwards merge by virtue of using < instead of <=, and always inserts to the innermost list possible instead of being strict about the bottom side. Which is understandable as being strict has a major downside. If you get a string of equal values in the bottom half of the range, you'd be stuck putting each one in a new sublist until you get to the midpoint and can create a new one where these values start going to the top.

Whoa, hey Tim! Very much agree on no "best" sorting algorithm. The generic-comparison case is one I keep idly considering, although more focused on long lists than short ones. There's been a lot of improvement lately in sorting small datatypes (arguably kicked off by fluxsort), but really, who has so many 4-byte ints laying around that sorting them is a bottleneck? A point Orson Peters made when presenting glidesort[0] was that lists with many repeated duplicates are common even in generic data, as in ordering customers by city, and quicksort seems to be the best way to deal with such a list. However glidesort/driftsort is still largely oriented towards smaller types, and its mergesort-quicksort hybridization can pretty easily make bad choices, as I described at [1] and especially [2] (slightly surprised I didn't mention binary insertion sort in that first section). So this approach feels like it's still immature to me.

[0] https://archive.fosdem.org/2023/schedule/event/rust_glidesor...

[1] https://mlochbaum.github.io/BQN/implementation/primitive/sor...

[2] https://mlochbaum.github.io/BQN/implementation/primitive/sor...

Ah, the README prior to commit b98f724 does claim O(n log(1.5 ln(n))) asymptotic runtime. This is based on a flawed analysis assuming rainbow endpoints will be evenly distributed, which the paper presents before describing the problem that each inserted element pushes one of those endpoints away from the middle. Bit of a perpetual motion machine search: initial results showed the promise of this approach, some details got in the way but surely there's some trick to work around them...

Still there's the benchmark. I remembered timsort.txt shows comparison counts relative to the log_2(n!) bound and it's only like 1% worse, so the >30% improvement over Timsort at the high end has to come from some other factor. I'm thinking it is actually cache effects—not because the algorithm has inherently better cache properties, as it's just a merge sort at large scales, but because it copies the array to ints internally rather than using a generic comparison. That would be a poor result, as good comparison sorts compiled for 4-byte data are many times faster than Timsort.

The analysis would make a lot more sense if it dropped every big O and computed an expected number of comparisons. Of course, the real issue is that the whole thing relies on an empirically-determined length of √8√n for the number of lists in the rainbow, when the theoretical justification isn't there and clearly the worst case is n/2. I'm not seeing any awareness that the expected number of comparisons has to be at least log_2(n!), and am starting to see some worrying signs that the author thinks averages like n log log n are possible. I'd initially assumed that "somewhere between O(n) and O(n log_2 n)" statement meant in the presence of asymptotically large natural runs, but there's a claim that the runs "need not be of any minimum size" and another that the algorithm can "capitalize on natural runs" that seems to be referring to the benchmarks on "purely random values".

So the idea is that you'll eventually arrange n elements into ideally about √n sorted lists, organized so that the endpoints of these lists are also ordered (the "rainbow" shape). A new element goes on the end of one of those lists or a new list; which one is generally decided by binary search, but the index is saved to speed up the next insertion if it's the same. So a natural run might benefit by having adjacenct values inserted to the same list.

Which cases benefit and which don't seems potentially quite erratic. If you have a series of runs that shrink over the course of the array, with increasing start value and decreasing start value, each will become a new sorted list and you end up with O(n) insertion plus a perfect merge. If the runs get wider instead, they'll cross over existing endpoints and end up spattered around before merging.

It feels like the paper keeps trying to spin constant-factor improvements as somehow being asymptotically relevant despite acknowledging that mathematically they aren't, for example the fact that log(√n) is only log(n)/2 "would not reflect the sublinear nature of base array length k" and runs hitting the best case "approximately 50% of the time... reinforces JesseSort's average runtime complexity as somewhere between O(n) and O(n log_2 n)". Powersort provably takes O(n + n log r) worst-case time where r is the number of runs; I expect this algorithm is O(n log n) worst- and average-case.

However, on random arrays the given benchmarks do improve relative to Python default as the length increases, with faster performance starting a little under a million elements. There isn't a lot of investigation into why this is. I'd be interested to know if the number of comparisons is reduced in particular. Perhaps this algorithm is more cache-friendly with a random input, but this wouldn't be reliable as many inputs can cause the number of sorted lists to scale with O(n) so that the binary searches need to work on a large array.

All right, just threw me off a little that you'd consider speculating or backwards decoding as I wouldn't expect them to be easier, or significantly faster (or maybe you consider parity-independence to be speculation? I can see it).