brhfl.com

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:

      s←1↓⍳26
      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:

      (2↓s)←(3|2↓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.

      s←1↓⍳26
      (2↓s)←((3|2↓s)>0)×2↓s
      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

      (1↓s)←((2|1↓s)>0)×1↓s
      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 TryAPL.org.

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.