50847534

⎕io←0 throughout.

I was re-reading A Mathematician’s Apology before recommending it to Suki Tekverk, our summer intern, and came across a statement that the number of primes less than 1e9 is 50847478 (§14, page 23). The function pco from the dfns workspace does computations on primes; ¯1 pco n is the number of primes less than n:

      )copy dfns pco

      ¯1 pco 1e9
50847534

The two numbers 50847478 and 50847534 cannot both be correct. A search of 50847478 on the internet reveals that it is incorrect and that 50847534 is correct. 50847478 is erroneously cited in various textbooks and even has a name, Bertelsen’s number, memorably described by MathWorld as “an erroneous name erroneously given the erroneous value of π(1e9) = 50847478.”

Although several internet sources give 50847534 as the number of primes less than 1e9, they don’t provide a proof. Relying on them would be doing the same thing — arguing from authority — that led other authors astray, even if it is the authority of the mighty internet. Besides, I’d already given Suki, an aspiring mathematician, the standard spiel about the importance of proof.

How do you prove the 50847534 number? One way is to prove correct a program that produces it. Proving pco correct seems daunting — it has 103 lines and two large tables. Therefore, I wrote a function from scratch to compute the number of primes less than , a function written in a way that facilitates proof.

sieve←{
  b←⍵⍴{∧⌿↑(×/⍵)⍴¨~⍵↑¨1}2 3 5
  b[⍳6⌊⍵]←(6⌊⍵)⍴0 0 1 1 0 1
  49≥⍵:b
  p←3↓{⍵/⍳⍴⍵}∇⌈⍵*0.5
  m←1+⌊(⍵-1+p×p)÷2×p
  b ⊣ p {b[⍺×⍺+2×⍳⍵]←0}¨ m
}

      +/ sieve 1e9
50847534

sieve ⍵ produces a boolean vector b with length such that b/⍳⍴b are all the primes less than . It implements the sieve of Eratosthenes: mark as not-prime multiples of each prime less than ⍵*0.5; this latter list of primes obtains by applying the function recursively to ⌈⍵*0.5. Some obvious optimizations are implemented:

  • Multiples of 2 3 5 are marked by initializing b with ⍵⍴{∧⌿↑(×/⍵)⍴¨~⍵↑¨1}2 3 5 rather than with ⍵⍴1.
  • Subsequently, only odd multiples of primes > 5 are marked.
  • Multiples of a prime to be marked start at its square.

Further examples:

      +/∘sieve¨ ⍳10
0 0 0 1 2 2 3 3 4 4

      +/∘sieve¨ 10*⍳10
0 4 25 168 1229 9592 78498 664579 5761455 50847534

There are other functions which are much easier to prove correct, for example, sieve1←{2=+⌿0=(⍳3⌈⍵)∘.|⍳⍵}. However, sieve1 requires O(⍵*2) space and sieve1 1e9 cannot produce a result with current technology. (Hmm, a provably correct program that cannot produce the desired result …)

We can test that sieve is consistent with pco and that pco is self-consistent. pco is a model of the p: primitive function in J. Its cases are:
   pco n   the n-th prime
¯1 pco n   the number of primes less than n
 1 pco n   1 iff n is prime
 2 pco n   the prime factors and exponents of n
 3 pco n   the prime factorization of n
¯4 pco n   the next prime < n
 4 pco n   the next prime > n
10 pco r,s a boolean vector b such that r+b/⍳⍴b are the primes in the half-open interval [r,s).

      ¯1 pco 10*⍳10      ⍝ the number of primes < 1e0 1e1 ... 1e9
0 4 25 168 1229 9592 78498 664579 5761455 50847534

      +/ 10 pco 0 1e9    ⍝ sum of the sieve between 0 and 1e9
50847534
                         ⍝ sum of sums of 10 sieves
                         ⍝ each of size 1e8 from 0 to 1e9
      +/ t← {+/10 pco ⍵+0 1e8}¨ 1e8×⍳10
50847534
                         ⍝ sum of sums of 1000 sieves
                         ⍝ each of size 1e6 from 0 to 1e9
      +/ s← {+/10 pco ⍵+0 1e6}¨ 1e6×⍳1000
50847534

      t ≡ +/ 10 100 ⍴ s
1
                         ⍝ sum of sums of sieves with 1000 randomly
                         ⍝ chosen end-points, from 0 to 1e9
      +/ 2 {+/10 pco ⍺,⍵}/ 0,({⍵[⍋⍵]}1000?1e9),1e9
50847534

      ¯1 pco 1e9         ⍝ the number of primes < 1e9
50847534
      pco 50847534       ⍝ the 50847534-th prime
1000000007
      ¯4 pco 1000000007  ⍝ the next prime < 1000000007
999999937
      4 pco 999999937    ⍝ the next prime > 999999937
1000000007     
      4 pco 1e9          ⍝ the next prime > 1e9
1000000007

                         ⍝ are 999999937 1000000007 primes?
      1 pco 999999937 1000000007
1 1
                         ⍝ which of 999999930 ... 1000000009
                         ⍝ are prime?
      1 pco 999999930+4 20⍴⍳80
0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1

      +/∘sieve¨ 999999937 1000000007 ∘.+ ¯1 0 1
50847533 50847533 50847534
50847534 50847534 50847535

Shuffle Faster Please!

Andy reported that in the shuffle QA some functions take a long time:

m9249   “4½ days so far”
rankop  21.5 hours
m12466  26.3 hours
points  7.2 hours

Background: Shuffle QA is an intensification of the normal QA. The suite of over 1800 QA functions is run under shuffle, whereby every getspace (memory allocation) is followed by every APL array in the workspace being shifted (“shuffled”) and by a simple workspace integrity check. The purpose is to uncover memory reads/writes which are rendered unsafe by events which in normal operations are rare.

m9249

m9249 tests +\[a], ×\[a], ⌈\[a], and ⌊\[a]. It ran in under 3 seconds in normal operations; for it to take more than 4½ days under shuffle, getspace is clearly the dominant resource and hereafter we focus on that.

m9249 used expressions like:

assert (+⍀x) ≡ {⍺+⍵}⍀x

Since {⍺+⍵} is not recognized by the scan operator and therefore not supported by special code, {⍺+⍵}⍀x is done by the general O(n*2) algorithm for scans, AND on each scalar the operand {⍺+⍵} is invoked and a getspace is done. The improved, faster QA does this:

F ← {1≥≢⍵:⍵ ⋄ x ⍪ (¯1↑⍵) ⍺⍺⍨ ¯1↑x←⍺⍺ ∇∇ ¯1↓⍵}
assert (+⍀x) ≡ +F x

With such a change, the faster m9249 runs in 22 minutes under shuffle instead of over 4½ days, without any reduction in coverage.

The number of getspace are as follows: +⍀x takes O(1). The second expression {⍺+⍵}⍀x does ×\1↓⍴x separate scans, and each one takes O((≢x)*2) getspace; the total is therefore O(((≢x)*2)××/1↓⍴x) ←→ O((≢x)××/⍴x). The third expression +F is linear in ≢x and is therefore O(≢x) in the number of getspace. In m9249 the largest x has size 33 5 7. The following table shows the orders of complexity and what they imply for this largest x:

+⍀x O(1) 1
{⍺+⍵}⍀x O((≢x)××/⍴x) 38115 = 33 × ×/ 33 5 7
+F x O(≢x) 33

The ratio between 38115 and 33 goes a long way towards explaining why the time to run m9249 under shuffle was over 4½ days and is now 22 minutes. (Why does m9249 still take 22 minutes? It tests for the four functions + × ⌈ ⌊, for various axes, for different datatypes, and for different axis lengths.)

rankop

Another shuffle sluggard rankop took 47 hours. rankop tests expressions of the form f⍤r⊢yand x f⍤r⊢y for various functions f. The techniques used for reducing the number of getspace differ from function to function. We look at x↑⍤r⊢y as an illustration.

assert (4↑⍤1⊢y) ≡ 4{⍺↑⍵}⍤1⊢y

⍴y was 2 7 1 8 2 8. The expression took 7.6e¯5 seconds normally and 7.5 seconds under shuffle. The left hand side requires O(1) getspace, the right hand side O(×/¯1↓⍴x). The expression was changed to:

assert (4↑⍤1⊢y) ≡ 2 7 1 8 2 4↑y

Feels like cheating, doesn’t it? 🙂 The new expression takes 0.115 seconds under shuffle.

m12466

m12466 tests n⌈/[a]x (and n⌊/[a]x). Previously, it did this by comparing n⌈/[a]x against n{⍺⌈⍵}/[a]x for x with shape 7 17 13 11, for each of the four axes, for two different values of n, and for the 1-, 2-, 4-, and 8-byte numeric datatypes. It required 5238426 getspace and 26.3 hours to run under shuffle.

F←{p⍉⍺⍺⌿(⊂(⎕io-⍨⍳|⍺)∘.+⍳(1-|⍺)+(⍴⍵)[⍵⍵])⌷(⍋p←⍒⍵⍵=⍳⍴⍴⍵)⍉⍵}

F is a dyadic operator. The left operand ⍺⍺ is or ; the right operand ⍵⍵ is the axis; and n(⌈ F a)x is the same as n⌈/[a]x. n(⌈ F a)x mainly does n⌈⌿x; other axes are handled by transposing the axis of interest to be the first, then doing n⌈⌿x, then transposing the axes back into the required order. n(⌈ F a)x is much more complicated than n{⍺⌈⍵}/[a]x but has the advantage of using only O(1) getspace.

The revamped m12466 function uses 13717 getspace and 2 minutes to run under shuffle.

points

points is a function to test monadic format () for different values of ⎕pp.

points←{⍺←17 ⋄ ⎕pp←⍺
  tt←{(⌽∨\⌽⍵≠⍺)/⍵}                   ⍝ trim trailing ⍺s.    
  ~0∊∘.{                             ⍝ all pairs (⍳⍵)       
    num←'0'tt↑,/⍕¨⍺'.'⍵              ⍝ eg '9.3'             
    num≡⍕⍎num                        ⍝ check round trip     
  }⍨⍳⍵
}

The expression is 14 15 16 17 points¨ 100, meaning that for each of the 100 numbers num as text,

  1.1   1.2   1.3   1.4   1.5     1.100
  2.1   2.2   2.3   2.4   2.5     2.100
  3.1   3.2   3.3   3.4   3.5 ... 3.100
...
 99.1  99.2  99.3  99.4  99.5    99.100
100.1 100.2 100.3 100.4 100.5   100.100

a “round-trip” num≡⍕⍎num is evaluated. The expression requires 1.3 million getspace and 7.2 hours to run under shuffle.

A much more efficient version with respect to getspace is possible. Let x be a vector. ⍕¨x requires O(≢x) getspace; ⍕x requires O(1) getspace and, like ⍕¨x, formats each individual element of x independently.

points1←{
  ⎕io←⎕ml←1                                                 
  tt←{(⌽∨\⌽⍵≠⍺)/⍵}                   ⍝ trim trailing ⍺s     
  x←1↓∊(' ',¨s,¨'.')∘.,'0'tt¨s←⍕¨⍳⍵  ⍝ numbers as text      
  y←⍎x                                                      
  {⎕pp←⍵ ⋄ x≡⍕y}¨⍺                   ⍝ check round trip     
}

14 15 16 17 points1 100 requires 11050 getspace and less than 2 minutes to run under shuffle.

What to do?

Andy says that the shuffle QA takes over 2000 hours (!). The shuffle sluggards will be tackled as done here by using more parsimonious APL expressions to compare against the aspect being QA-ed. Going forward, a new QA function will be run under shuffle as it is being developed. The developers will be unlikely to tolerate a running time of 4½ days.

The Halting Problem Rendered in APL

The halting problem is the problem of determining, from a description of a program and an input, whether the program will finish running or continue to run forever. It was proved in the negative by Alonzo Church and Alan Turing in 1936, and can be rendered in Dyalog APL as follows:

Suppose there is an operator H such that f H ⍵ is 1 or 0 according to whether f ⍵ halts. Consider:

In dfns

   f←{∇ H ⍵:∇ ⍵ ⋄ ⍬}

For f 0, if ∇ H 0 is 1, then ∇ 0 is invoked, leading to an infinite recursion, therefore ∇ H 0 should have been 0. On the other hand, if ∇ H 0 is 0, then the part is invoked and is the result of f 0, therefore ∇ H 0 should have been 1.

Presented as a table: For f 0

suppose ∇ H 0 is invokes consequence ∇ H 0 should be
1 ∇ 0 infinite recursion 0
0 f 0 results in 1

Therefore, there cannot be such an operator H.

In Tradfns

   ⎕fx 'f x' '→f H x'

For f 0

suppose f H 0 is invokes consequence f H 0 should be
1 →1     infinite loop 0
0 →0 f 0 exits 1

Quicksort in APL Revisited

A message in the Forum inquired on sorting strings in APL with a custom comparison function.

First, sorting strings without a custom comparison function obtains with a terse expression:

      {⍵[⍋↑⍵]} 'syzygy' 'boustrophedon' 'kakistocracy' 'chthonic'
┌─────────────┬────────┬────────────┬──────┐
│boustrophedon│chthonic│kakistocracy│syzygy│
└─────────────┴────────┴────────────┴──────┘

Sorting strings with a custom comparison function can also be accomplished succinctly by simple modifications to the Quicksort function

Q←{1≥≢⍵:⍵ ⋄ S←{⍺⌿⍨⍺ ⍺⍺ ⍵} ⋄ ⍵((∇<S)⍪=S⍪(∇>S))⍵⌷⍨?≢⍵}

in a Dyalog blog post in 2014. A version that accepts a comparison function operand is:

QS←{1≥≢⍵:⍵ ⋄ (∇ ⍵⌿⍨0>s),(⍵⌿⍨0=s),∇ ⍵⌿⍨0<s←⍵ ⍺⍺¨ ⍵⌷⍨?≢⍵}

The operand function ⍺ ⍺⍺ ⍵ is required to return ¯1, 0, or 1, according to whether is less than, equal to, or greater than . In QS, the phrase s←⍵ ⍺⍺¨ ⍵⌷⍨?≢⍵ compares each element of against a randomly chosen pivot ⍵⌷⍨?≢⍵; the function then applies recursively to the elements of which are less than the pivot (0>s) and those which are greater than the pivot (0<s).

Example with numbers:

      (×-)QS 3 1 4 1 5 9 2 6 5 3 5 8 97 9
1 1 2 3 3 4 5 5 5 6 8 9 9 97

      3(×-)4
¯1
      3(×-)3
0
      3(×-)¯17
1

Example with strings:

      strcmp←{⍺≡⍵:0 ⋄ ¯1*</⍋↑⍺ ⍵}

      'syzygy' strcmp 'syzygy'
0
      'syzygy' strcmp 'eleemosynary'
1
      'syzygy' strcmp 'zumba'
¯1

      strcmp QS 'syzygy' 'boustrophedon' 'kakistocracy' 'chthonic'
┌─────────────┬────────┬────────────┬──────┐
│boustrophedon│chthonic│kakistocracy│syzygy│
└─────────────┴────────┴────────────┴──────┘

A more efficient string comparison function is left as an exercise for the reader. 🙂

Permutations

I started composing a set of APL exercises in response to a query from a new APL enthusiast who attended Morten’s presentation at Functional Conf, Bangalore, October 2014. The first set of exercise are at a low level of difficulty, followed by another set at an intermediate level. One of the intermediate exercises is:

Permutations

a. Compute all the permutations of ⍳⍵ in lexicographic order. For example:

   perm 3
0 1 2
0 2 1
1 0 2
1 2 0
2 0 1
2 1 0

b. Write a function that checks whether is a solution to perm ⍺, without computing perm ⍺. You can use the function assert. For example:

assert←{⍺←'assertion failure' ⋄ 0∊⍵:⍺ ⎕SIGNAL 8 ⋄ shy←0}

pcheck←{
  assert 2=⍴⍴⍵:
  assert (⍴⍵)≡(!⍺),⍺:
  …
  1
}

   6 pcheck perm 6
1

   4 pcheck 2 4⍴0 1 2 3, 0 1 3 2
assertion failure
pcheck[2] assert(⍴⍵)≡(!⍺),⍺:
         ∧

c. What is the index of permutation in perm ⍺? Do this without computing all the permutations. For example:

   7 ip 1 6 5 2 0 3 4    ⍝ index from permutation
1422

(The left argument in this case is redundant, being the same as ≢⍵.)

d. What is the -th permutation of ⍳⍺? Do this without computing all the permutations. For example:

   7 pi 1442             ⍝ permutation from index
1 6 5 2 0 3 4

   (perm 7) ≡ 7 pi⍤0 ⍳!7
1

The Anagram Kata

Coincidentally, Gianfranco Alongi was attempting in APL the anagrams kata from Cyber Dojo:

Write a program to generate all potential anagrams of an input string.

For example, the potential anagrams of “biro” are
biro bior brio broi boir bori
ibro ibor irbo irob iobr iorb
rbio rboi ribo riob roib robi
obir obri oibr oirb orbi orib

This is essentially the same program/exercise/kata, because the potential anagrams are 'biro'[perm 4]. You can compare solutions in other languages to what’s here (google “anagrams kata”).

Spoiler Alert

Deriving a Solution

I am now going to present solutions to part a of the exercise, generating all permutations of ⍳⍵.

Commonly, in TDD (test-driven development) you start with a very simple case and try to extend it successively to more general cases. It’s all too easy to be led into a dead-end because the simple case may have characteristics absent in a more general case. For myself, for this problem, I would start “in the middle”: Suppose I have perm 3, obtained by whatever means:

   p
0 1 2
0 2 1
1 0 2
1 2 0
2 0 1
2 1 0

How do I get perm 4 from that? One way is as follows:

   p1←0,1+p
   (0 1 2 3[p1]) (1 0 2 3[p1]) (2 0 1 3[p1]) (3 0 1 2[p1])
┌───────┬───────┬───────┬───────┐
│0 1 2 3│1 0 2 3│2 0 1 3│3 0 1 2│
│0 1 3 2│1 0 3 2│2 0 3 1│3 0 2 1│
│0 2 1 3│1 2 0 3│2 1 0 3│3 1 0 2│
│0 2 3 1│1 2 3 0│2 1 3 0│3 1 2 0│
│0 3 1 2│1 3 0 2│2 3 0 1│3 2 0 1│
│0 3 2 1│1 3 2 0│2 3 1 0│3 2 1 0│
└───────┴───────┴───────┴───────┘

So it’s indexing each row of a matrix m by 0,1+p. There are various ways of forming the matrix m, one way is:

   ⍒⍤1∘.=⍨0 1 2 3
0 1 2 3
1 0 2 3
2 0 1 3
3 0 1 2

(Some authors waxed enthusiastic about this “magical matrix”.) In any case, a solution obtains readily from the above description: Form a matrix from the above individual planes; replace the 0 1 2 3 by ⍳⍵; and make an appropriate computation for the base case (when 0=⍵). See the 2015-07-12 entry below.

The Best perm Function

What is the “best” perm function I can write in APL? This “best” is a benchmark not only on my own understanding but also on advancements in APL over the years.

“Best” is a subjective and personal measure. Brevity comes into it but is not the only criteria. For example, {(∧/(⍳⍵)∊⍤1⊢t)⌿t←⍉(⍵⍴⍵)⊤⍳⍵*⍵} is the shortest known solution, but requires space and time exponential in the size of the result, and that disqualifies it from being “best”. The similarly inefficient {(∧/(⍳⍵)∊⍤1⊢t)⌿t←↑,⍳⍵⍴⍵} is shorter still, but does not work for 1=⍵.

1981, The N Queens Problem

    p←perm n;i;ind;t
   ⍝ all permutations of ⍳n
    p←(×n,n)⍴⎕io
    →(1≥n)⍴0
    t←perm n-1
    p←(0,n)⍴ind←⍳n
    i←n-~⎕io
   l10:p←(i,((i≠ind)/ind)[t]),[⎕io]p
    →(⎕io≤i←i-1)⍴l10

It was the fashion at the time that functions be written to work in either index-origin and therefore have ⎕io sprinkled hither, thither, and yon.

1987, Some Uses of { and }

   perm:  ⍪⌿k,⍤¯1 (⍙⍵-1){⍤¯ 1 k~⍤1 0 k←⍳⍵
       :  1≥⍵
       :  (1,⍵)⍴0

Written in Dictionary APL, wherein: ⍪⌿⍵ ←→ ⊃⍪⌿⊂⍤¯1⊢⍵ and differs from its definition in Dyalog APL; is equivalent to in dfns; ⍺{⍵ ←→ (⊂⍺)⌷⍵; and ¯ by itself is infinity.

1990-2007

I worked on perm from time to time in this period, but in J rather than in APL. The results are described in a J essay and in a Vector article. The lessons translate directly into Dyalog APL.

2008, http://dfns.dyalog.com/n_pmat.htm

   pmat2←{{,[⍳2]↑(⊂⊂⎕io,1+⍵)⌷¨⍒¨↓∘.=⍨⍳1+1↓⍴⍵}⍣⍵⍉⍪⍬}

In retrospect, the power operator is not the best device to use, because the left operand function needs both the previous result (equivalent to perm ⍵-1) and . It is awkward to supply two arguments to that operand function, and the matter is finessed by computing the latter as 1+1↓⍴⍵.

In this formulation, ⍉⍪⍬ is rather circuitous compared to the equivalent 1 0⍴0. But the latter would have required a or similar device to separate it from the right operand of the power operator.

2015-07-12

   perm←{0=⍵:1 0⍴0 ⋄ ,[⍳2](⊂0,1+∇ ¯1+⍵)⌷⍤1⍒⍤1∘.=⍨⍳⍵}

For a time I thought the base case can be ⍳1 0 instead of 1 0⍴0, and indeed the function works with that as the base case. Unfortunately (⍳1 0)≢1 0⍴0, having a different prototype and datatype.

Future

Where might the improvements come from?

  • We are contemplating an under operator whose monadic case is f⍢g ⍵ ←→ g⍣¯1 f g ⍵. Therefore 1+∇ ¯1+⍵ ←→ ∇⍢(¯1∘+)⍵
  • Moreover, it is possible to define ≤⍵ as ⍵-1 (decrement) and ≥⍵ as ⍵+1 (increment), as in J; whence 1+∇ ¯1+⍵ ←→ ∇⍢≤⍵
  • Monadic = can be defined as in J, =⍵ ←→ (∪⍳⍨⍵)∘.=⍳≢⍵ (self-classify); whence ∘.=⍨⍳⍵ ←→ =⍳⍵

Putting it all together:

   perm←{0=⍵:1 0⍴0 ⋄ ,[⍳2](⊂0,∇⍢≤⍵)⌷⍤1⍒⍤1=⍳⍵}

We should do something about the ,[⍳2] :-)​​

In Praise of Magic Functions: Part II

Part I of this post was concerned with the development speed and execution speed of magic functions and should be read before this post.

Benefitting from Future and Past Improvements

Looking at the magic function for key in Part I, we see that its performance depends on the following APL primitives, listed with information on when they were last improved (or when they will be improved):

expression factor version
⍋i 2-18 12.1
⎕DR 1-3.8 14.1
⍳⍨x 2-4 14.1
⍳⍨i 3-6 14.1
↑x 5-12 15.0
⊂[a]x 1-1.75 15.0
b⊂[⎕IO]x 1.3-42 14.1

Key was introduced in version 14.0 and x{⍺,f⌿⍵}⌸y is currently computed by the general case. Comparing the speed in version 14.0 to 15.0 as of June 2015:

   x←?1e6⍴5e4
   y←?1e6⍴0

   1 1 cmpx 'x{⍺,+/⍵}⌸y'    ⍝ 14.0
0.1735

   1 1 cmpx 'x{⍺,+/⍵}⌸y'    ⍝ 15.0
0.12855

   0.1735 ÷ 0.12855
1.34967

a factor of 1.3 improvement without touching the key code. Special code had been planned for {⍺,f/⍵}⌸, and we can get an idea of the amount of speed up:

   cmpx 'x{⍺,+/⍵}⌸y' '(∪x),⍤0⊢x{+/⍵}⌸y'    ⍝ 15.0
 x{⍺,+/⍵}⌸y       → 1.25E¯1 |   0% ⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕
 (∪x),⍤0⊢x{+/⍵}⌸y → 1.16E¯2 | -91% ⎕⎕⎕

Perhaps the special code for {⍺,f/⍵}⌸ can be a magic function?

Magic functions also benefit from past improvements. The primitive ⍕b (in the Suggestivity subsection below) has a magic function implementation: MAGIC("(¯1+2×≢⍵)⍴(2 2⍴'0 1 ')[⍵;]⊣⎕io←0"). Then I thought, why not get an extra level of performance by doing it in native C? Benchmarks done after that work showed that the native C implementation was faster by a factor of 25 but the magic function was faster by a factor of 40. What happened? On reflection, the speed of the code depends crucially on the speed of x[b;…;] and it turned out that x[b;…;] was previously sped up (version 13.2). Apparently the x[b;…;] code is faster than a casually-written C program.

Finally, I expect magic functions are prime candidates for being fed to the APL compiler when it becomes possible to do so.

APL as a Tool of Implementation

Being APL code, magic functions have the characteristic terseness of APL and all that that entails. K.E. Iverson’s 1979 Turing Lecture Notation as a Tool of Thought lists five important characteristics of notation. We examine how they play out in this context.

Ease of Expression
Notation as a Tool of Thought entitles this characteristic as “ease of expressing constructs arising in problems”, explicated as follows (§1.1):

If it is to be effective as a tool of thought, a notation must allow convenient expression not only of notions arising directly from a problem, but also of those arising in subsequent analysis, generalization, and specialization.

For magic functions, the analogous concept would be facility in program development, debugging, tracing, profiling, benchmarking, and other activities required in implementation. Dyalog APL has these in spades.

Suggestivity
It was first identified that ⍕b can be sped up by doing (¯1+2×≢⍵)⍴(2 2⍴'0 1 ')[⍵;]⊣⎕io←0, the idea being that if f is an expensive function and x is from a small domain D, then f x can be computed by (f D)[⍵;]. (A similar idea is described in §4.3 of Notation as a Tool of Thought.) For the idea to be effective, f does not have to be very expensive, just more expensive than indexing, and x is sufficiently large. A few speed-ups in the pattern have been identified:

⍕b
x⍕b
float ⍳ b
float ⍳ int8
scalar×b and b×scalar
scalar*b

All of these can be implemented by magic functions. As well, the utility of the pattern (f D)[i;] motivates extra attention on x[i;…;]. Fortuitously, x[b;…;] has previously been made fast.

Subordination of Detail
Writing in APL means not having to deal with many details involved in writing in C in the interpreter, including:

  • workspace compaction
  • arrays changing their internal datatype
  • call by name
  • allocating and releasing memory
  • index errors AKA memory read and write exceptions
  • different numeric types
  • different character types
  • loops and nested loops
  • declarations
  • etc.

Economy
Economy refers to the size of the vocabulary and the complexity of the grammar. Here, the comparison is not just between APL and C (in which comparison APL wins), but between APL and the collection of programs, utilities, macros, data structures, typedefs, calling conventions, …, accumulated by the C source code over the years. These have not been as rigorously defined and executed as APL.

Amenability to Formal Manipulation
The opening paragraph of the present text had a different magic function for ∧.= in a previous version. In repeated readings of the MSS and in meditation on the ideas I realized that the line of code can be made simpler:

MAGIC("((≢⍺)↑i)∘.=(≢⍺)↓i←⍺⍳⍺⍪⍉⍵");    ⍝ old version
MAGIC("(≢⍺)(↑∘.=↓)⍺⍳⍺⍪⍉⍵");           ⍝ current version

I daresay that such simplification would be harder to conceive with a more verbose statement of the algorithm.