Bubble sort in dc

dc offers essentially nothing in the way of stack manipulation. Part of the power behind a lot of RPN/RPL HP calculators was that if you were being tidy with your stack, you rarely had to use registers. This is not the case with dc. Of course there’s also nothing resembling array sorting in dc, and I was curious how difficult it would be to implement a basic bubble sort.

As mentioned, we can’t directly manipulate stack elements (beyond swapping the two topmost values), but we can directly address elements of an array. This means that our bubble sort needs to dump the contents of the stack into an array, and potentially dump this array back onto the stack at the end (depending on what we need the sort for).


This takes care of putting everything on the stack into array s. It does this by storing the number of items on the stack into register a (which we will need later – dc does not have a means of counting the number of items in an array). We start counter i at zero, and just go one at a time putting the top-of-stack into array s at index i until i is equal to a.


The bubble sort itself requires three macros. This first one, M is our main macro, which we will ultimately begin executing to start the sort process. It sets a counter, i to one and a check register, c, to zero. Then it runs our next macro, S, which does one pass through the array. If S has changed anything by running macro R, the check register c will no longer be zeroed out, so we test for this and restart M until no changes have been made.


As mentioned, macro S does one pass of the array s. We fetch our counter, i, duplicate it twice, and decrement the top copy by one. While we work on our array, we’re always looking at i and i-1 essentially, which just makes the comparison a bit tidier at the end vs. i+1 and i. We load the values with these indices from s, compare them, and if they aren’t in order, we run macro R. We still have a copy of i on the stack at this point, so we increment it, duplicate it, store it, and compare with a to see if we’ve hit the end of the array.


Macro R does one instance of reversal. First it puts a one in the check register, c, so that M knows changes have been made to the array. At this point there is still a copy of i on the stack, which we duplicate three times. We load the i-indexed value from s, swap with a lower copy of i which we decrement by one before loading this i-1-indexed value from s. We load i again, and store the old i-1 value (our previous top-of-stack) at this index. This leaves our old i value atop the stack, with i one spot below it. We swap them, subtract one from i, and store the old i value at this new i-1 index in array s.

And that’s actually it. If we populate the stack and run our initial array-populating code before executing lMx, we’ll end up with our values sorted in array s. From here, we can [li1-d;srdsi0<U]dsUx to dump the array onto the stack such that top-of-stack is low, or 0si[lid;sr1+dsila>V]dsVx to dump the array onto the stack such that top-of-stack is high. If, say, we only need min (0;s) or max (la1-;s), these are simple tasks.

Additionally, if we wanted to get the median of these values, we can exploit the fact that dc just uses the floor of a non-whole index into an array. This allows us to avoid an evenness test by always taking ‘two’ values from the middle of the array and calculating their mean. If a is odd, then half of a will end in .5. Conversely, if a is even, half of a will end with no remainder. So we can subtract half of a by .5 to get the lower bound, and then subtract half of a by half of a mod one (either .5 or 0) and average the values at the resulting two indices: la2/d1~-;sr.5-;s+2/p.

Sieve of Eratosthenes

This post contains APL code, and assumes that your system will be capable of rendering the relevant unicode code points. If the code below looks like gibberish, either you don’t understand APL or your computer doesn’t ☺. I use the APL standard of six spaces indentation to show my input, and zero indentation for output. As for the code itself, it assumes one-indexing, that is, ⎕IO←1.

I was messing around with some primality and coprimality testing in dc the other day when I got to wondering about inefficient methods for checking primality (specifically, the thought of testing primality of n by testing coprimality of n and m where m is every integer<n). This reminded me of the sieve of Eratosthenes, a first-century CE means of narrowing down a list of integers to a list of primes by eliminating multiples of primes (which must be composites). My APL is getting very rusty, unfortunately, but this seemed like a fun exercise since APL is a language heavily invested in arrays. We may start by assigning s to the integers 2-26, and confirm this by printing s as a 5x5 matrix:

      5 5 ⍴s
 2  3  4  5  6
 7  8  9 10 11
12 13 14 15 16
17 18 19 20 21
22 23 24 25 26

Then we can start sieving. I’m going to skip ahead to the threes to avoid a potential spot of confusion, so what we need to do is reassign the elements after 3 to either a zero if there’s no remainder (upon dividing by 3), or the original value if there is. 2↓s narrows us down to the appropriate elements, and we have to make sure that we do that everywhere so that our 1x25 shape stays intact. The first step is getting our remainders. 3|2↓s performs modulus 3 on elements 4-25, resulting in 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2. This doesn’t fit neatly into a 5x5 grid, however, so I’m going to temporarily write it over the corresponding elements in s:

      5 5 ⍴s
2 3 1 2 0
1 2 0 1 2
0 1 2 0 1
2 0 1 2 0
1 2 0 1 2

This gives us zeroes where we need them – at every multiple of three (excluding three itself, of course). We can then turn this into a series of ones and zeroes by comparing each value to zero:

      5 5 ⍴ 0<s
1 1 1 1 0
1 1 0 1 1
0 1 1 0 1
1 0 1 1 0
1 1 0 1 1

Which we could then multiply by our original set of 25 integers, if we still had them. So let’s reassign s and put it all together. And since we skipped over two, we should probably do that as well.

      5 5 ⍴s
 2  3  4  5  0
 7  8  0 10 11
 0 13 14  0 16
17  0 19 20  0
22 23  0 25 26

      5 5 ⍴s
 2  3  0  5  0
 7  0  0  0 11
 0 13  0  0  0
17  0 19  0  0
 0 23  0 25  0

So now we’ve sieved out all of our multiples of 2 and 3, and the only thing left to sieve is 5. Of course, if we didn’t already know the primes ≤25, we’d want to keep trying every nonzero value in the list, but we do know that 25 is the only outstanding composite in the list, and (4↓s)←((5|4↓s)>0)×4↓s does, in fact, take care of it, as can be seen here, on

I mentioned that my APL skills are a bit rusty. I’m not sure I’ve even mentioned APL on this blog before, but it is probably tied with Forth for the coveted title of Bri’s favorite language. I’ve never been great at it, though. It’s a relic, with lineage dating before computers were capable of running it. I either don’t do enough array-based math, or I don’t think of enough of the math I do in array-based terms to use APL regularly. Where I get really rather rusty is in thinking outside of imperative terms or using APL’s imperative programming behaviors. Which was why my little demo here was very REPL-esque. Hopefully a future post will bundle this thing up into an actual, runnable function that knocks the whole process out in one go. But for now, have fun typing in those iotas.

Post updates

I’ve realized that a few of the things I’ve written over the past year may have contained a nugget or two of information that was either ill-informed, or otherwise deserving of an update. For starters, in the posts Darwin image conversion via sips and Semaphore and sips redux, I talk about using Darwin’s sips command to convert a TIFF to PNG before running it through optipng. While this is a fine exercise, I have since come to learn that optipng will handle uncompressed TIFFs just fine, converting them to PNG and optimizing accordingly. sips is unnecessary, so long as I’m willing to temporarily use up the SSD space for uncompressed TIFFs.

Before those posts was Of lynx and curl, describing a few uses of those tools in tandem. Toward the end I mention going through a list of URLs, and writing out a CSV by using printf to output the URL before using curl to fetch the status code. For some reason, I had long overlooked that curl has a variable for the current URL as well, so the printf step is largely unnecessary. Going through a list of links to output a CSV of the form URL,Status Code can be achieved as for i in $(< links.txt); curl -o /dev/null --location --silent --head --write-out '%{url_effective},%{http_code}\n' "$i" >> status.csv (given that links.txt is a line-by-line list of URLs to check).

In other news, the first post I made on this version of this blog was the meta post, Licensing – dated 2016-07-30. So I’ve stuck with this thing for a year now, 75 posts, so far so good. I did recently upgrade from Hugo 0.16 to 0.24.1, which I have a post-in-progress about but all in all the upgrade went shockingly smooth. I definitely have no regrets about the move to a static site generator, and I would whole-heartedly recommend Hugo for anyone whose needs it meets. I have made a few other minor changes, like setting my Inline audio player to not autoload, and customizing the selection highlight color, but nothing major since my last round of Template updates. Next up is hopefully building a better archive page using new Hugo features, having Hugo generate the CSS files, and possibly a JSON feed.


Despite having failed out of geometry in my younger days, it has become my favorite sort of recreational math. I think, back in elhi geometry, too much emphasis was placed on potential practical applications instead of just distilling it down to the reality that this was what math was before there was any sort of formal mathematical system. Geometry is hands-on, it’s playful, and as such I have come to really enjoy playing with it. As much as I enjoy doing constructions with a straightedge and compass, I occasionally poke around to see what tools exist on the computer as well. Recently, I stumbled across a very neat thing: Eukleides.

I’m drawn to Eukleides because it is a language for geometry, and not a mouse-flinging WYSIWYG virtual compass. This seems contradictory given my gushing about geometry being hands-on, and don’t get me wrong – I love a hands-on GUI circle-canvas too1. But sometimes (often times) my brain is in code-mode and it’s easier to express what I’m trying to do in words than to fiddle around with a mouse. And a task like ‘intersecting a couple of circles’ is far more conducive to writing out than, say, laying down an SVG illustration from scratch.

a b c d e bi

There you have one of the first constructions learned by anyone studying geometry – bisecting an angle with three arcs (or, full-blown circles in this case). Angle ∠abc is bisected by segment bbi. Here’s the code:

% Percent signs are comments
a=point(7,0); b=point(0,0) % Semicolons and newlines separate commands
a b c triangle 5,-40°
d=intersection(g,a.b); e=intersection(g,b.c)
d=d[0]; e=e[0] % Intersections return sets (arrays), extract the points
h=circle(d,3); i=circle(e,3)
  a 0°; b 180°; c 40° % Label the points
  d -40° gray; e 90° gray
  bi 20°
  a,b,bi; bi,b,c 1.5 % Make the angle markers
  g lightgray; h gray
  i gray

Note that the code originally used shades of grey, I shifted these around to match my site’s colors when I converted the resulting EPS file to SVG. The code is pretty straightforward: define some points, make an angle of them, draw a circle, intersect segments ab and bc, make some more circles, intersect where the circles meet, and boom – a bisected angle. The language reminds me a bit of GraphViz/DOT – purpose-built for naturally expressing how things will be drawn.

We can actually prove that the construction works without even generating an image file. Removing the draw and label sections, and replacing them with some print (to stdout) statements calling angle measurement functions:

a=point(7,0); b=point(0,0)
a b c triangle 5,-40°
d=intersection(g,a.b); e=intersection(g,b.c)
h=circle(d,3); i=circle(e,3)
%%%%%%%% New content starts here:
print angle(a,b,c)
print angle(a,b,c)/2
print angle(a,b,bi)
print angle(bi,b,c)

We get ∠abc, ∠abc/2, ∠abbi, and ∠bcbi. The last three of these should be equal, and:


…they are! We can also do fun things like dividing the circumference of a circle by its diameter:

print (perimeter(circle(point(0,0),4))/8)

…to get 3.14159. Very cool. There are a few improvements I would like to see. Notably, while you can label points etc. with their names, I can’t find a way to add arbitrary labels, or even labels based on functions. So you can’t (seemingly) label a line segment with its length, or an angle with its measure. Also, the interpreter wants ISO 8859-1 encoding, and it uses the degree symbol2 (°) to represent angle measures. This gets all flippy-floppy when moving to UTF-8, and in fact if I forget to convert a file to ISO 8859-1, I’ll get syntax errors if I use degree symbols. Finally, it only outputs to EPS; native SVG would be incredibly welcome.

Eukleides is a lot of fun to play with, and it’s worth mentioning that it has looping and conditionals and the like, so it is programmable and not just computational or presentational. Presumably some pretty interesting opportunities are thus opened up.

A call to donate to Lambda Legal (external)

Well, bad horrifying news out of the White House today, no surprise there. Lambda Legal is a 501(c)3 legal organization fighting for the civil rights of queer and HIV positive folks. They do good work, and have already committed to fight this act of bigotry. Link in the header is to their donate page, they’re top on my to-donate list at the moment. Additionally:

Try to take care of yourselves, all. ⚧