Diagonalize, Now!

It is easy to append two lists, and the method is well-known:

let rec append xs ys = match xs with
| []      -> ys
| x :: xs -> x :: append xs ys

It is also easy to merge two lists, so that elements from the first are interleaved with those of the second:

let rec merge xs ys = match xs with
| []      -> ys
| x :: xs -> x :: merge ys xs

A simple matter of swapping two arguments! The similarities between both functions hide deep differences. In a strict setting, the equation append xs ys = xs has just a solution ys = [] if xs is finite; but in a lazy setting this equation has infinite solutions for ys if xs is infinite, because A^ω⋅A* ≡ A^ω. By contrast, merge xs ys = xs has as its only solution ys = [] independently of whether xs is finite or not (in a lazy language).

What this all means is that merge is fairer than append, in this sense: merge xs ys includes all the elements of both xs and ys. This is, in transfinite terms, one half of the constructive proof that ω + ω > ω but 2⋅ω = ω, the other being an effective procedure for un-merge-ing an infinite list (edit: thanks Andrej Bauer for your thorough comments!). In practice, it is the reason why you cannot enumerate all the integers by first counting the nonnegatives 0, 1, 2, 3… and then the negatives -1, -2, -3… but you can if you interleave them 0, -1, 1, -2, 2, -3, 3…

The downside to this is that append is associative but merge is not, so it doesn't form a monoid over lists (even though [] is its left and right identity) and thus neither endows the lists with monadic structure.

Can we generalize to lists of lists? Sure we can:

let rec concat = function
| []        -> []
| xs :: xss -> append xs (concat xss)

(a simple consequence of lists being a monoid), or:

let concat xss = List.fold_right (fun xs ys -> append xs ys) xss []

If we try to do the same with merge:

let rec diagonalize = function
| []        -> []
| xs :: xss -> merge xs (diagonalize xss)

again, a simple change (I omit the folded version for brevity). The problem is that, since merge is not associative, the result of diagonalize [xs0; xs1; xs2; …] is not very fair: it includes every other element of xs₀ followed by every fourth element of xs₁ followed by every eighth element of xs₂…, so that you have to wait 2^{n + 1} elements before finding the next element in xs_n! The result does include every element in every list even if there are an infinite number of infinite lists, because ½ + ¼ + ⅛ + … = 1, but you have to wait exponentially long for them. Georg Cantor devised a clever way out of this: read the (possibly doubly infinite) matrix:

(thank you codecogs for the wonderful L^AT_EX rendering service!) along ascending diagonals, as if it were rotated 45° to the right: x₀₀, x₀₁, x₁₀, x₀₂, x₁₁, x₂₀, x₀₃, x₁₂, x₂₁, x₃₀, … You get the idea. Now you only have to wait quadratically long before you see an element. However, how should we proceed in order to flatten a list of lists in this order? By induction, of course!

Suppose every column in the matrix is a list, and that we can "half-transpose" such a list of lists so that each list in the result corresponds to a diagonal in the original matrix: [[x00]; [x01; x10]; [x02; x11; x20]; [x03; x12; x21; x30]; …]. Call this operation rotate. Then, by elementwise consing the head list with each list in the rotation of the rest we diagonalize the entire list of lists. The code for pairing up these lists is:

let rec pairup xs yss = match xs, yss with
| [], yss -> yss
| xs, []  -> List.map (fun x -> [x]) xs
| x :: xs, ys :: yss -> (x :: ys) :: pairup xs yss

(being careful not to lose any element!). The code to perform the rotation is:

let rec rotate = function
| []        -> []
| xs :: xss -> pairup xs ([] :: rotate xss)

Note that the corner element gets to be by itself. Then the diagonalization is simply the flattening of this rotation:

let diagonalize xss = concat (diag xss)

The question now is, how to generalize this to higher dimensions (lists of lists of lists of …)? Since nowhere to be seen is the nice symmetry between appending and merging here I don't know what the guiding principle could be.

You are Here

I am a terrible person: I read code, I think to myself "I can do better" and I set out thus, out of pride. I did just that with an OCaml implementation of the geohash geographical encoding. Even if the code I read was perfectly workable, I found it unidiomatic and somewhat bloated. I wrote two versions, the first one compact but impenetrable, the second one from first principles: what would be an executable description of the algorithm? To wit:

let decode =
   (narrow (-90., 90.) *** narrow (-180., 180.)) % swap % split
  % List.concat % List.map (to_bits 5 % decode_base32) % explode

Because of (or despite!) the point-free, combinator-heavy style I intended this to be read as a pipeline or as a recipe, from right to left (it is an applicative-order pipeline), much as the algorithm description does:

1. Explode the code string into a list of characters
2. For each character, decode it as a base-32 digit decomposing it into a list of 5 bits
3. Flatten the list of lists of bits into a big list of bits
4. Split the list of bits into a pair of lists for the even bits and the odd bits
5. Swap both components of the pair, since the coding uses the even bits for the longitude
6. Narrow each list independently into a coordinate starting from the corresponding interval

The result is a (latitude, longitude) pair. Supposing that every step is invertible, the algorithm for encoding a coordinate into a geohash is:

let encode nchars =
  let bits = 5 * nchars in
  let lo = bits / 2 in
  let hi = bits - lo in
    implode % List.map (encode_base32 % of_bits) % group 5 % join
  % swap % (expand lo (-90., 90.) *** expand hi (-180., 180.))

Again, this is easy to follow as a recipe or as a succint description:

1. Given the desired hash length, compute the number of bits for the latitude and for the longitude
2. Expand each coordinate independently into a list of bits describing its position in the corresponding interval
3. Swap both components of the pair
4. Join both lists by alternating even and odd elements
5. Group the bits in the list into sublists by taking 5 bits at a time
6. For each group, convert the list of 5 bits into an integer and encode it as a base-32 digit
7. Implode the list of characters into a code string

Now I have a number of inverse or near-inverse pairs of functions to write, namely:

• explode and implode
• decode_base32 and encode_base32
• to_bits and of_bits
• List.concat and group
• split and join
• narrow and expand

Let's start with the (well-known to Haskellers) combinators (though they pronounce % as .):

let ( $ ) f x = f x
let ( % ) f g x = f (g x)

The first is the application operator and the second is the composition operator. Next are some operators on pairs borrowed from Control.Arrow:

let first f (x, y) = (f x, y)
let second f (x, y) = (x, f y)
let ( *** ) f g (x, y) = (f x, g y)
let swap (x, y) = (y, x)

They allow applying a function to one or both members of a pair, and to manipulate its components. Note that swap ∘ swap ≡ id. This completes the functional scaffolding. The first and more imperative pair of functions is explode and implode:

let explode s =
  let res = ref [] in
  for i = String.length s - 1 downto 0 do
    res := String.unsafe_get s i :: !res
  done;
  !res

let implode l =
  let b = Buffer.create 10 in
  List.iter (Buffer.add_char b) l;
  Buffer.contents b

The first traverses the string from right to left to accumulate each character on the head of the resulting list; the second uses a string buffer to build the resulting string a character at a time. I use unsafe_get because, well, these are as low-level functions as they come, so squeezing a bit extra from them doesn't seem inappropriate. Each is the other's inverse, or in more abstract terms the pair (explode, implode) witnesses the monoid isomorphism between OCaml strings and lists of characters (that's why Haskell identifies both datatypes).

Aside: I won't attempt any proofs to avoid making a long post even longer. Some are easy, most are tedious, and the preceding is probably difficult. I've verified a few of the inverse pairs but not all of them, so I don't consider this code verified by any stretch of the word. I am, however, confident that the code works.

Now decode_base32and encode_base32 are also fairly imperative, but just because I chose to use an array as the map between characters and their base-32 encodings:

let codetable = [|
'0'; '1'; '2'; '3'; '4'; '5'; '6'; '7';
'8'; '9'; 'b'; 'c'; 'd'; 'e'; 'f'; 'g';
'h'; 'j'; 'k'; 'm'; 'n'; 'p'; 'q'; 'r';
's'; 't'; 'u'; 'v'; 'w'; 'x'; 'y'; 'z'
|]

This means that for encoding a 5-bit number indexing into the array is enough, but decoding a base-32 character requires a lookup. Since the array is sorted by character, a binary search is a good choice:

let binsearch ?(compare=Pervasives.compare) v e =
  let i = ref 0
  and j = ref (Array.length v) in
  if !j = 0 then raise Not_found else
  while !i <> !j - 1 do
    let m = (!i + !j) / 2 in
    if compare v.(m) e <= 0
      then i := m
      else j := m
  done;
  if v.(!i) = e then !i else raise Not_found

let encode_base32 i = codetable.(i land 31)
and decode_base32 c = binsearch codetable (Char.lowercase c)

(This binary search is completely general). Now decode_base32 ∘ encode_base32 ≡ id (for the restricted domain of 5-bit integers) if and only if binsearch v v.(i) = i for i in range of v. Similarly, encode_base32 ∘ decode_base32 ≡ id (again, for the restricted domain of base-32 characters) if and only if v.(binsearch v x) = x. The conditions on binsearch amount to it being correct (and would make for excelent unit or random testing).

The next inverse pair is to_bits and of_bits:

let iota =
  let rec go l i =
    if i <= 0 then l else let i = i - 1 in go (i :: l) i
  in go []

let to_bits k n = List.rev_map (fun i -> (n land (1 lsl i)) != 0) $ iota k

let of_bits = List.fold_left (fun n b -> 2 * n + if b then 1 else 0) 0

The first maps a list of bit positions (given by iota, the initial natural interval) into the corresponding bits by direct testing, most-significant (highest) bit first; the second computes the binary value of the list of bits by a Horner iteration.

The following functions are very general functions on lists. Any function f : α list → α list list such that List.concat ∘ f ≡ id is called a partition. One such function is group:

let group k l =
  let rec go gs gss i = function
  | []           -> List.rev (if gs = [] then gss else List.rev gs :: gss)
  | l when i = 0 -> go [] (List.rev gs :: gss) k l
  | x :: xs      -> go (x :: gs) gss (i - 1) xs
  in go [] [] k l

It collects up to k elements in the stack gs, pushing complete groups into the stack gss. The last group might be incomplete and it is added to the top of the stack before returning it. A function that I will need later is the the converse of List.fold_right, appropriately called unfold:

let rec unfold f e n =
  if n <= 0 then [] else
  let (x, e') = f e in
  x :: unfold f e' (n - 1)

unfold generates a list of the given length n from a seed e and a function f that computes the head and the next seed. Now split and join are another very functional couple:

let cons x xs = x :: xs

let split l =
  let rec even = function
  | [] -> [], []
  | x :: xs -> first (cons x) $ odd xs
  and odd = function
  | [] -> [], []
  | x :: xs -> second (cons x) $ even xs
  in even l

let rec join (xs, ys) = match xs with
| [] -> ys
| x :: xs -> x :: join (ys, xs)

In general, split ∘ join is not the identity on pairs of lists (think of ([], l) for nonempty l), but indeed join ∘ split ≡ id, and a proof could proceed by induction. Now all these functions are completely general (except perhaps for the base-32 encoding and decoding) and have little if anything to do with geocaching. The last pair, narrow and expand, are the meat of the algorithm:

let avg (x, y) = (x +. y) /. 2.

let bisect (lo, hi as interval) b =
  let md = avg interval in
  if b then (md, hi) else (lo, md)

let narrow interval = avg % List.fold_left bisect interval

narrow repeatedly bisects the interval, keeping the upper or lower half depending on the next bit in the input list, finally returning the midpoint of the final interval. Its inverse is expand:

let divide x (lo, hi as interval) =
  let md = avg interval in
  if x > md
    then true , (md, hi)
    else false, (lo, md)

let expand bits interval x = unfold (divide x) interval bits

It unfolds each bit in turn, by divideing the interval and finding which half contains x. These functions are inverses in the sense that, for any interval i = (min, max) and list l, expand (List.length l) i (narrow i l) = l and conversely, for any n ≥ 0 and min ≤ x ≤ max:

|narrow i (expand n i x) - x| ≤ (max - min)2-n-1

These are the crucial properties that make encode the true inverse of decode. The astute reader might have noticed the resemblance of geohashing with arithmetic coding: in a sense, a geohash is the base-32 encoding of the coordinates lossly compressed with an uniform probability model. One way to improve the precision of the encoding would be by allocating a predefined weight to each hemispherical quadrant, by taking into account that more landmass is present in the north-eastern quadrant of the globe. This would require splitting the intervals not along the middle but, say, the first third, thus representing more precise locations with the same number of bits in continental Europe and the Far East, then in North America, then in Africa and Oceania, and lastly in South America.