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.