Wanted: working Bose-Hibbard Sort implementation p

Please point me to code for a working Bose-Hibbard Sort implementation, preferably in a C-like language.

I'm trying to implement the algorithm in C#, but I don't have a copy of the algorithm. The only sample I have is a Fortran implementation that is a "free translation" of Hibbard's original Algol version (from 'A simple sorting algorithm' Journal of the ACM vol 10 (1963) p142-50 — which I don't have either). The Fortran version appears to be buggy (it does exactly 1 compare and ends up exiting if they are already sorted) so my primary focus is to get a working version to study.

From a scanned PDF of the original article (downloaded from ACM Digital Library), OCR'd by copy'n'paste on a Mac, and then manually cleaned up (quite a lot):

procedure ternary sort (a, n); array a; integer n; begin integer j, L;
integer array x, y[0: log2(n-1)] ; integer procedure sum(w); array w;
begin integer j, s; s := 0; for j:= 0 step 1 until L do s := s+w[j]×2↑j; sum := s
end; procedure compare; begin real w;
if a[sum(x)] > a[sum(y)] then begin w := a[sum(x)]; a[sum(x)] := a[sum(y)];
a[sum(y)] := w end end;
L := entier log2(n-1); for j := 0 step 1 until L do begin x[j] := 0;
y[j] := if j = 0 then 1 else 0 end;
A: compare; j := 0;
C: go to if x[j] = y[j] = 0 then zero else if x[j] = y[j] = 1 then one else
if x[j] = 0 then first two else two;
zero: x[j] := y[j] := 1; if sum(y) ≤ n-1 then go to A;
one: y[j] := 0; go to A;
two: x[j] := 0; j := j+1; go to C;
first two: x[j] := y[j] := 0; if j = L then go to exit; j := j+1;
if y[j] = 1 then go to D; x[j] := y[j] := 1; if sum(y) > n-1 then
go to first two; if sum(y) < n-1 then j := 0;
D: x[j] := 0; y[j] := 1; go to A;
exit: end

In the original, the 'log2' functions are set as 'log₂'. Line breaks as in the original. This predates the 'structured programming' revolution - now you can see why structured programming is a good idea. It also predates careful, clear code layout. It is interesting seeing 'two word' labels and procedure names. (In the Revised Report for Algol-60 (PDF, or HTML), it says: Typographical features such as blank space or change to a new line have no significance in the reference language. They may, however, be used freely for facilitating reading. This means that what appears to be 'two words' in modern computer languages is just one word in Algol-60; searching with Google shows that the keywords were differentiated from variables etc by being underlined or printed in bold or enclosed in quotes of some sort. The language also used a number of symbols not normally found on keyboards; three examples in this program are '×', '↑', and '≤'.)

With the nested procedures, you'd need quite a lot of 'free translating' to code this in Fortran.

Here it is re-formatted - it is perhaps a little easier to see what the code is; the plethora of 'go to' statements does not make it easier to understand. Now you can see why Dijkstra wrote 'GOTO Considered Harmful'.

procedure ternary sort (a, n); array a; integer n; 
begin
    integer j, L;
    integer array x, y[0: log2(n-1)]; 
    integer procedure sum(w); array w;
    begin
        integer j, s; 
        s := 0; 
        for j:= 0 step 1 until L do s := s+w[j]×2↑j; 
        sum := s
    end;
    procedure compare;
    begin
        real w;
        if a[sum(x)] > a[sum(y)] then
        begin
            w := a[sum(x)]; 
            a[sum(x)] := a[sum(y)];
            a[sum(y)] := w 
        end 
    end;
    L := entier log2(n-1);
    for j := 0 step 1 until L do
    begin
        x[j] := 0;
        y[j] := if j = 0 then 1 else 0
    end;
A:  compare; j := 0;
C:  go to if x[j] = y[j] = 0 then zero
          else if x[j] = y[j] = 1 then one
          else if x[j] = 0 then first two 
          else two;
zero:
    x[j] := y[j] := 1;
    if sum(y) ≤ n-1 then go to A;
one:
    y[j] := 0; 
    go to A;
two:
    x[j] := 0; 
    j := j+1; 
    go to C;
first two:
    x[j] := y[j] := 0;
    if j = L then go to exit; 
    j := j+1;
    if y[j] = 1 then go to D; 
    x[j] := y[j] := 1; 
    if sum(y) > n-1 then go to first two; 
    if sum(y) < n-1 then j := 0;
D:  x[j] := 0;
    y[j] := 1; 
    go to A;
exit:
end