Observable sharing and executable proofs

Every rational number has an infinite repeating decimal expansion. In particular, every integer admits two decimal expansions: one with and one without a "tail of nines". This means that 1 = 0.9999…, a provable fact that many people still choose to disbelieve. What do I mean here by "provable"? I mean "provable by finitistic reasoning", that is, by a purely algebraic proof free of notions of continuity or convergence. The proof in this particular (and particularly popular) case goes like this:

    x     = 0.9999…
≡ { product by base is shift of decimal point }
  10x     = 9.9999…
≡ { subtraction of finitely-expanded integer }
  10x - 9 = 0.9999…
≡ { identity of decimal expansion }
  10x - 9 = x
≡ { algebra }
   9x     = 9
≡ { cancellation of multiplication }
    x     = 1

The last two steps are elementary consequences of the fact that integers form a ring. The first three steps require justification. By definition, a decimal expansion of a positive real number x is a (possibly semi-infinite) sequence of digits 0 ≤ d_i < 10 such that for all i ≥ N, d_i = 0, and

x = ⟨∑ i :: 10i⋅di⟩

and the decimal point separates the digits corresponding to 10⁰ and 10^-1, namely d₀ from d_-1.

Lemma 0 Any decimal number x can be expressed as the sum of an integer and a (possibly infinite) decimal part

Proof: By definition, x has a finite number of decimal digits to the left of the decimal point, namely for 0 ≤ i < N, hence:

  x
= { definition }
  ⟨∑ i :: 10i⋅di ⟩
= { associativity of addition, a finite number N of times }
  ⟨∑ i : 0 ≤ i < N : 10i⋅di⟩ + ⟨∑ i : 0 > i : 10i⋅di⟩
= { naming the first term }
  n + ⟨∑ i : 0 > i : 10i⋅di⟩

This justifies representing a positive real number as a pair of an integer and a possibly infinite purely fractional decimal expansion:

type dec = Dec of int * int list

A sufficient condition for two numbers to be equal is that their decimal representations are the same:

let eq (Dec (e, ds)) (Dec (e', ds')) = e = e' && ds == ds'

Here, since the decimal expansions could be infinite, I use physical equality on lists. This is a stronger condition than stated above. Now, the proof for the first step:

Lemma 1 Multiplying a positive real number x by 10 shifts the decimal point in its decimal expansion one place to the right

Proof:

  10⋅x
= { lemma 0 }
  10⋅(n + ⟨∑ i : 0 > i : 10i⋅di⟩)
= { product finitely distributes over sum }
  10⋅n + 10⋅⟨∑ i : 0 > i : 10i⋅di⟩
= { product formally distributes over sum }
  10⋅n + ⟨∑ i : 0 > i : 10i+1⋅di⟩
= { substitute i := i-1 }
  10⋅n + ⟨∑ i : 1 > i : 10i⋅di-1⟩
= { lemma 0 }
  (10⋅n + d0) + ⟨∑ i : 0 > i : 10i⋅di-1⟩

Note that the third step is not strictly speaking finitist but formal: it requires taking the distributivity of product over a summation as axiomatic if induction is not accepted. This justifies the following:

let mul10 = function
| Dec (e, []     ) -> Dec (10 * e    , [])
| Dec (e, d :: ds) -> Dec (10 * e + d, ds)

The proof for the second step:

Lemma 2 Subtracting an integer k from a positive real number x only affects its integer part

Proof:

  x - k
= { lemma 0 }
  (n + ⟨∑ i : 0 > i : 10i⋅di⟩) - k
= { commutativity and associativity of sum }
  (n - k) + ⟨∑ i : 0 > i : 10i⋅di⟩

This justifies the following:

let subi (Dec (e, ds)) n =
  if e < n then failwith "subi" else
  Dec (e - n, ds)

Now the proof for the third step is purely executable and relies on physical equality of infinite lists of digits represented as circular lists:

let nines = let rec ds = 9 :: ds in Dec (0, ds)

Of course, for more general repeating expansions, equality on heads won't suffice, but this case is nicely handled by a one-liner:

let lemma = assert (eq nines (subi (mul10 nines) 9))

The purely verbal mathematical statement of this truth is that, since the first two steps operate on a finite prefix of the infinite expansion, the inifinite suffix is unchanged by the operations and hence identical to itself. In other words, the infinite summation is a purely formal object that remains unchanged from step to step.

Will this proof finally lay the "controversy" to rest, persuading even the staunchest finitist?

Edit: Thanks to id that alerted me in my confusion between structural and physical equality.

Properly Bound

To all practitioners: the type of bind is not the type of >>=. For a monad α m, the latter has type:

val (>>=) : α m → (α → β m) → β m

only because >>= is used as a combinator that allows pipelining computations from left to right, as per usual convention:

… m >>= fun x →
  n >>= fun y →
  …
  return f x y

On the other hand, if you want to take maximum advantage of the Theorems for Free! then your bind should have type:

val bind : (α → β m) → (α m → β m)

because it is a natural transformation, together with return:

val return : α → α m

You can see immediately how both "fit" together; in particular, the second monad law (right identity) falls off naturally (!) from the types:

bind return ≡ id

The first monad law (left identity) is also immediately natural:

bind f ∘ return ≡ f

since bind f's domain coincides with return's range by inspection. The third monad law (associativity) is much less regular but you can see a hint of the mechanics of adjointedness if you read bind g as a whole:

bind g ∘ bind f ≡ bind (bind g ∘ f)

This note is motivated by a reading of Jérémie Dimino's slides about Jérôme Vouillon's Lwt, specifically slide 6 (Edit: thanks Gabriel for the heads up about proper attribution).

The elusive Binary Search

Jon Bentley has scarred a couple of generations of programmers, or he was and continues to be right:

I've assigned this problem in courses at Bell Labs and IBM. […] In several classes and with over a hundred programmers, the results varied little: ninety percent of the programmers found bugs in their programs (and I wasn't always convinced of the correctness of the code in which no bugs were found).

Jon Bentley, Programming Pearls, p. 36

This disjunction is not meant to be exclusive. The question whether programmers know how to write a binary search or not crops again in connection with this quote (but the topic is perennial), and reddit jumps to try its collective hand at the task and fail with numbers echoing Bentley's. The remedy (namely, closely study Dijkstra's A Discipline of Programming, or at least the van Gasteren and Feijen's paper) is as effective as it is bitter. The alternative is to learn the algorithm by heart:

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

(note: the midpoint calculation cannot overflow in OCaml). No ifs, no buts, not a single thing changed: the algorithm as presented above is proved correct. Alas, we live in an age of practical assurances, and like Thomas we trust more a green-lighted battery of tests than a formal derivation (who understands those anyway?). There is a nice, 4096-case test set available on-line; a parser and a test harness are simple enough to build for it. Each case consists of a label, a test value, a list of values and an expected outcome:

type test = Test of string * int * int array * bool

The file format is line-oriented; this makes easy to use a regexp-based parser. To ease the use of regular expressions, I build a test function out of a pattern:

let matching patt =
  let re = Str.regexp patt in
  fun str -> Str.string_match re str 0

(I close over the compiled regexp for efficiency). The required regexps are:

let is_label  = matching "^Problem [0-9]+$"
let is_start  = matching "^in \\[$"
let is_value  = matching "^-?[0-9]+$"
let is_endyes = matching "^\\]\\? yes$"
let is_endno  = matching "^\\]\\? no$"
let is_sep    = matching "^$"

Syntax errors in the test cases file are signalled by an exception:

exception SyntaxError of int * string

I keep track of the lines read for purposes of error reporting:

let lineno = ref 0

let syntax_error str = raise (SyntaxError (!lineno, str))

let read inch =
  try let l = input_line inch in incr lineno; Some l
  with End_of_file -> None

The parser reads a file and returns a list of test cases; it uses mutually-recursive functions as a traditional recursive-descent parser with synthetic attributes. Each function corresponds to a terminal in the production case ::= label value start value* end sep:

let rec parse_tests inch = List.rev (parse_case inch [])

and parse_case inch tests = match read inch with
| Some line when is_label  line -> parse_member inch tests line
| None -> tests
| _ -> syntax_error "PROBLEM expected"

and parse_member inch tests label = match read inch with
| Some line when is_value  line -> parse_start  inch tests label (int_of_string line)
| _ -> syntax_error "MEMBER expected"

and parse_start inch tests label cand = match read inch with
| Some line when is_start  line -> parse_value  inch tests label cand []
| _ -> syntax_error "START expected"

and parse_value inch tests label cand vals = match read inch with
| Some line when is_endyes line -> parse_finish inch tests label cand vals true
| Some line when is_endno  line -> parse_finish inch tests label cand vals false
| Some line when is_value  line ->
  parse_value  inch tests label cand (int_of_string line :: vals)
| _ -> syntax_error "VALUE or END expected"

and parse_finish inch tests label cand vals test = match read inch with
| Some line when is_sep    line ->
  parse_case   inch (Test (label, cand, Array.of_list (List.rev vals), test) :: tests)
| _ -> syntax_error "SEPARATOR expected"

The test harness iterates a test function over the resulting list, recording the number of successful tests along the total:

let test_binsearch inch =
  let count = ref 0
  and pass  = ref 0 in
  List.iter (fun (Test (label, cand, vals, test)) ->
    let result =
      try ignore (binsearch vals cand); test
      with Not_found -> not test
    in
    incr count;
    if result then begin
      Printf.printf "TEST %s: PASS\n" label;
      incr pass
    end else
      Printf.printf "TEST %s: FAIL\n" label
    ) (parse_tests inch);
  Printf.printf "\nSUMMARY: Passed %d of %d tests (%.2f %%)\n"
    !pass !count (100. *. float !pass /. float !pass)

In order to be tidy about using files, I use a helper handler:

let unwind ~protect f x =
  try let y = f x in let () = protect x in y
  with e -> let () = protect x in raise e

let with_input_channel f inch = unwind ~protect:close_in f inch

The main expression turns this little (less than 90 lines of code) program into an executable suitable for command-line execution:

let () =
  let inch =
    try let fname = Array.get Sys.argv 1 in open_in_bin fname
    with Invalid_argument _ -> stdin
  in with_input_channel test_binsearch inch

That's it. Compile and execute it with:

$ ocamlopt.opt -o binsrch.exe str.cmxa binsrch.ml
$ ./binsrch.exe tests.txt > tests.out

(I'm using Cygwin here, don't hate me). The result is as expected:

TEST Problem 1: PASS
TEST Problem 2: PASS
TEST Problem 3: PASS
[snip…]
TEST Problem 4094: PASS
TEST Problem 4095: PASS
TEST Problem 4096: PASS

SUMMARY: Passed 4096 of 4096 tests (100.00 %)

Now do yourself a favor and go read Dijkstra.

GADTs for the Rest of Us

In a nutshell, generalized algebraic data types, or GADTs for short, are a first-class phantom type of sorts. The basic idea is to enrich the constructors of an ADT with a more expressive "return" type by allowing the instantiation of one or more type parameters. A simple motivating example is that of an abstract syntax in Peyton Jones's original paper:

data Term a where
  Lit  :: Int        -> Term Int
  Inc  :: Term Int   -> Term Int
  IsZ  :: Term Int   -> Term Bool
  If   :: Term Bool  -> Term a    -> Term a     -> Term a
  Pair :: Term a     -> Term b    -> Term (a,b)
  Fst  :: Term (a,b) -> Term a
  Snd  :: Term (a,b) -> Term b

For instance, the constructor for Inc specifies the type of its argument as Term Int by constraining its application to bind the parameter a to Int, in effect specifying the "return type" of the constructor. Xavier Leroy's proposal for a syntax extension of OCaml would look like this:

type 'a term =
| Lit  of int                            constraint 'a = int
| Inc  of int term                       constraint 'a = int
| IsZ  of int term                       constraint 'a = bool
| If   of bool term * 'a term * 'a term
| Pair of  'b  term * 'c term            constraint 'a = 'b * 'c
| Fst  of ('b * 'c) term                 constraint 'a = 'b
| Snd  of ('b * 'c) term                 constraint 'a = 'c

Each syntax emphasizes complementary views of what "is" the essence of this generalization: Haskell's makes clear the parameter's functional dependency on the constructor; OCaml's stresses the constraints placed on a match on the constructor. In any case, this type shows the three cases a GADT might present:

• No constraint, as in the constructor If
• A concrete constraint, as in the constructors Lit, Inc or IsZ
• An existential constraint, as in the constructors Pair, Fst or Snd

The constraints a GADT places on the parameters must be met by the functions that do a pattern-matching on the type. In general, the result type is a function of the types of the parameter. This requires polymorphic recursion which forces the use of an explicit rank-2 type for it. For instance, a function evaluating a Term a will have type a for each possible a defined in the GADT. In the following code the match for Lit forces a to be instantiated to Int, whereas the matching for IsZ forces it to Bool:

eval :: Term a -> a
eval term = case term of
  Lit n    -> n
  Inc t    -> eval t + 1
  IsZ t    -> eval t == 0
  If c x y -> if eval c then eval x else eval y
  Pair x y -> (eval x, eval y)
  Fst t    -> fst (eval t)
  Snd t    -> snd (eval t)

Unfortunately, the current version of OCaml (3.11) does not yet include either GADTs or rank-2 polymorphism. On July 11th, 2009, Oleg announced that OCaml's type system is expressive enough to encode GADTs. On his site he shows three applications of the technique. However, Oleg's examples do not offer much guidance in how to apply the encoding, especially since the more complete example, that of the typed printf/scanf, is too advanced an application to use as a template. I have shown in another post how to encode polymorphic recursion with rank-2 types. Here I'll explain how to treat each of the three cases in a GADT to systematically embed the type in OCaml.

The basic technique is to use a witness for the constraint placed on the parameter. Oleg uses a module with the required semantics:

module EQ = struct
  type ('a,'b) eq = Refl of 'a option ref * 'b option ref

  let refl () = let r = ref None in Refl (r, r)

  let symm (Refl (x, y)) = Refl (y, x)

  let apply_eq (Refl (x, y)) v =
    x := Some v;
    match !y with
  | Some w -> x := None; w
  | _      -> assert false
end

A value eq of type ('a, 'b) eq witnesses the equality of types a and b. Given such a witness and a value x of type a, apply_eq eq x is a safe type cast of x to type b. This is sufficient to force the constraints on the type parameter, as explained in the paper Oleg references. For each one of the above identified cases, Oleg's code encodes the constructor, defines a smart constructor to ease the use of the type as an embedded DSL and pattern-matches on the constructor in order to access its parameters. This is how:

• For unconstrained constructors, nothing needs to be done.

  For instance, the type constructor for If does not require fixing the constructor's parameter a as it is by definition equal to the type parameter:

  | If of bool term * 'a term * 'a term
  

  A smart constructor for this constructor is routine:

  let if_ c t f = If (c, t, f)
  

  (the underscore is to avoid a clash with the reserved keyword). Equally routine is pattern-matching on it:

  | If (c,t,f) -> …
  

  If desired for reasons of regularity, however, the constraint 'a = 'b for an existentially-qualified parameter b can be introduced explicitely using the last encoding.

• For concretely-constrained constructors, it is necessary to witness the equality of the parameter and this concrete type.

  It is not as simple, however, to augment the constructor with a value of type ('a, t) EQ.eq: since matches on those values would force the type of a to be the same for every match, it would require that t be the same concrete type everywhere. What is needed instead is an existential type, say b such that it is constrained to be t and which in turn constrains a. Oleg encodes the existential with the double negation on a rank-2 universal, by way of a witness type. In order to maximize the compactness of the code, he uses a class type to introduce the required rank-2 type:

  | Lit of < m_lit : 'w . ((int, 'a) eq -> int -> 'w) -> 'w >
  

  Here w witnesses the existence of an a such that, if successfully constrained to int, it would bear the required type int term. A smart constructor for this constructor must build the object instance which will provide the required witness:

  let lit n = Lit (object
    method m_lit : 'w . ((int, 'a) eq -> int -> 'w) -> 'w =
      fun k -> k (refl ()) n end)
  

  Note first that, under the Curry-Howard isomorphism, continuations encode double negation. Note also that the typing must be explicit since it relies on structual and not nominal class types. Note finally that the encoding style is curried much as in Haskell definitions; it is entirely equivalent to encode Leroy-style tupled constructors.

  Now a match on a constructor of this kind must invoke the method on the object parameter with a suitable continuation that will accept the witness as valid and coerce the argument to the constrained type:

  | Lit o -> o#m_lit (fun eq n -> apply_eq eq n)
  

  As will be apparent from the complete encoding of the Haskell example, the other cases are the same.

• For existentially-constrained constructors, it is necessary to encode the existential with a rank-2 type, and witness the equality of the parameter and this type.

  This case is an extension of the previous one, with the twist that the existential must be encoded separately. In this case Oleg uses a record type instead of an object type, but the encoding could use the latter equally well. For instance, in the case of Pairs, there must exist types b and c such that a is constrained to be b×c. The constructor factors out this existential:

  | Pair of < m_pair : 'w . ('a, 'w) k_pair -> 'w >
  

  Here, k_pair will bear the required witness w on a's constraint:

  type ('a, 'w) k_pair = { k_pair : 'b 'c . ('b*'c, 'a) eq -> 'b term -> 'c term -> 'w }
  

  The equality witness expresses the required intent on the types of the constructor's arguments. The smart constructor is analogous to the previous one, except that in this case the continuation has type k_pair:

  let pair p q = Pair (object
    method m_pair : 'w . ('a, 'w) k_pair -> 'w =
      fun k -> k.k_pair (refl ()) p q end)
  

  The pattern-match is analogous too:

  | Pair o -> o#m_pair { k_pair = fun eq x y -> apply_eq eq … }
  

Here's the complete code for the type definitions and the smart constructors encoding the Haskell mini-expression language example:

type 'a term =
| Lit  of < m_lit  : 'w . ((int , 'a) eq -> int -> 'w) -> 'w >
| Inc  of < m_inc  : 'w . ((int , 'a) eq -> int term -> 'w) -> 'w >
| IsZ  of < m_isz  : 'w . ((bool, 'a) eq -> int term -> 'w) -> 'w >
| If   of bool term * 'a term * 'a term
| Pair of < m_pair : 'w . ('a, 'w) k_pair -> 'w >
| Fst  of < m_fst  : 'w . ('a, 'w) k_fst  -> 'w >
| Snd  of < m_snd  : 'w . ('a, 'w) k_snd  -> 'w >
and ('a, 'w) k_pair = { k_pair : 'b 'c . ('b*'c, 'a) eq -> 'b term -> 'c term -> 'w }
and ('a, 'w) k_fst  = { k_fst  : 'b 'c . ('b   , 'a) eq -> ('b * 'c) term -> 'w }
and ('a, 'w) k_snd  = { k_snd  : 'b 'c . ('c   , 'a) eq -> ('b * 'c) term -> 'w }

let lit  n     = Lit  (object
  method m_lit  : 'w . ((int , 'a) eq -> int -> 'w) -> 'w =
    fun k -> k (refl ()) n end)
let inc  t     = Inc  (object
  method m_inc  : 'w . ((int , 'a) eq -> int term -> 'w) -> 'w =
    fun k -> k (refl ()) t end)
let isz  t     = IsZ  (object
  method m_isz  : 'w . ((bool, 'a) eq -> int term -> 'w) -> 'w =
    fun k -> k (refl ()) t end)
let if_  c t f = If   (c, t, f)
let pair p q   = Pair (object
  method m_pair : 'w . ('a, 'w) k_pair -> 'w =
    fun k -> k.k_pair (refl ()) p q end)
let first  p   = Fst  (object
  method m_fst  : 'w . ('a, 'w) k_fst  -> 'w =
    fun k -> k.k_fst  (refl ()) p end)
let second p   = Snd  (object
  method m_snd  : 'w . ('a, 'w) k_snd  -> 'w =
    fun k -> k.k_snd  (refl ()) p end)

It is important to remark on the types of the smart constructors, which are the most similar to the Haskell constructors:

val lit    : int -> int term
val inc    : int term -> int term
val isz    : int term -> bool term
val if_    : bool term -> 'a term -> 'a term -> 'a term
val pair   : 'a term -> 'b term -> ('a * 'b) term
val first  : ('a * 'b) term -> 'a term
val second : ('a * 'b) term -> 'b term

Now, as I mentioned above, an evaluator has type 'a term -> 'a for all a; this requires rank-2 polymorphism. In order to keep the namespace tidy, I use a local module definition hiding a record type encoding the higher-kinded function type:

let eval e =
  let module E = struct
    type eval = { eval : 'a . 'a term -> 'a }
    let rec eval = { eval = function
    | Lit  o      -> o#m_lit  (fun eq n -> apply_eq eq n)
    | Inc  o      -> o#m_inc  (fun eq t -> apply_eq eq (eval.eval t + 1))
    | IsZ  o      -> o#m_isz  (fun eq t -> apply_eq eq (eval.eval t = 0))
    | If  (c,t,f) -> if eval.eval c then eval.eval t else eval.eval f
    | Pair o      ->
      o#m_pair { k_pair = fun eq x y -> apply_eq eq (eval.eval x, eval.eval y) }
    | Fst  o      ->
      o#m_fst  { k_fst  = fun eq p   -> apply_eq eq (fst (eval.eval p)) }
    | Snd  o      ->
      o#m_snd  { k_snd  = fun eq p   -> apply_eq eq (snd (eval.eval p)) }
    }
  end in E.eval.E.eval e

In every case the coercion must be applied at the end, to the result of recursively evaluating the subexpressions in a constructor.

This is all there is to Oleg's encoding. For me, this exegesis has two lessons: first, much as rank-2 polymorphism, GADTs are a conservative extension to Haskell's and OCaml's type system. In fact, version 3.12 will allow higher-ranked functions via explicit annotations. One can only hope that the latter would also find its way into the next release. Second, that Oleg's proof-of-concept examples bear detailed examination to find the general mechanism behind his insights.

Braun Trees

I had a couple of interesting comments on my post about heaps and priority queues. The first is by Jérôme and is an indictment against my wrong, naïve use of Obj.magic to fake initializing extensible Vecs in the presence of double-word unboxed float arrays. There are a number of avenues to repair the defect:

1. Use a functorial interface and monomorphize the type of vectors
2. Use an exemplar of the intended type to initialize the backing array
3. Specialize the code for float arrays
4. Use an explicit box to store values, most probably an option

A second reader, Damien Guichard suggested sidestepping the problem entirely and using a dynamic data structure, Braun trees as a flexible array. There aren't many online references about the data structure, the most notable being Okasaki's "Three Algorithms on Braun Trees" (Functional Pearl), but Damien provided code that served as a starting point. A Braun tree is a complete binary tree with the property that for any given node the left subchild has exactly the same size of the right subchild or one more element than it. In other words, a Braun tree of size n + 1 is a node and two children l and r, both Braun trees, such that |l| = ⌈n/2⌉ and |r| = ⌊n/2⌋, which makes |r| ≤ |l| ≤ |r| + 1. As Okasaki remarks, Braun trees have always minimum height, and their shape is determined only by the number of nodes in them.

There is a choice between an imperative implementation and a purely functional one. Damien's code in the comment is imperative, and I found that rewriting it to introduce tail-recursion was tedious and error-prone; Okasaki's "exceptionally simple and elegant" algorithms devolve into a tangle of special cases. I opted by a purely functional implementation upon which to build an imperative, mutable interface compatible with the previous version of the priority queue. The important thing to consider is to maintain in every case the invariant for the trees. I'll use a different signature than before, one specific to this data structure:

module Braun : sig
  type 'a t
  val empty     : 'a t
  val singleton : 'a -> 'a t
  val is_empty  : 'a t -> bool
  val min       : 'a t -> 'a
  val get       : 'a t -> int -> 'a
  val size      : 'a t -> int
  val rep       : ('a -> 'a -> int) -> 'a -> 'a t -> 'a t
  val ins       : ('a -> 'a -> int) -> 'a -> 'a t -> 'a t
  val del       : ('a -> 'a -> int) -> 'a t -> 'a t
end = struct (* … *)

Note that rep, ins and del take a comparison function to maintain the heap property, namely, that the root of the tree is less than the elements in either children. This is unsafe as nothing precludes passing different comparison functions to manipulate the tree; this is why I will wrap it in the final Heap signature. That said, the constructor for trees is a straightforward binary tree ADT:

type 'a t = E | N of 'a * 'a t * 'a t

Construction, testing for emptiness and returning the minimum are trivial:

let empty       = E
and singleton x = N (x, E, E)

let is_empty = function E -> true | N _ -> false

let min = function
| N (e, _, _) -> e
| E           -> failwith "empty heap"

To calculate the size of a Braun tree I follow Okasaki verbatim:

let rec diff h n = match h with
| E           when n = 0 -> 0
| E                      -> assert false
| N (_, _, _) when n = 0 -> 1
| N (_, l, r)            ->
  if n mod 2 = 1 then diff l ((n - 1) / 2) else diff r ((n - 2) / 2)

let rec size = function
| E           -> 0
| N (_, l, r) -> let m = size r in 2 * m + 1 + diff l m

For the details I defer to the explanation in the paper; the idea is to avoid traversing the entire tree and exploiting the Braun property by calculating the size from the shape of each left child of the right spine; this reduces the complexity to O(log² n). For completeness, Braun trees support random access in logarithmic time:

let rec get h i = match h with
| E                      -> failwith "index out of bounds"
| N (e, _, _) when i = 0 -> e
| N (_, l, r)            ->
  if i mod 2 = 1 then get l ((i - 1) / 2) else get r ((i - 2) / 2)

The meat of the stew is the heap operations: insertion, and deletion or replacement of the minimum. All three operations must simultaneously keep two invariants: the Braun property and the heap property. Replacing the minimum doesn't change the shape of the tree, only the values stored in it, so it only has to maintain the heap property. The algorithm is the functional analog to the binary heap's siftdown:

let rec rep compare e = function
| E           -> failwith "empty heap"
| N (_, l, r) -> siftdown compare (N (e, l, r))

and siftdown compare n = match n with
| N (e, (N (el, _, _) as l), E) ->
  if compare e  el < 0 then n
    else N (el, rep compare e l, E)
| N (e, (N (el, _, _) as l), (N (er, _, _) as r)) ->
  if compare e  el < 0
  && compare e  er < 0 then n else
  if compare el er < 0
    then N (el, rep compare e l, r)
    else N (er, l, rep compare e r)
| _ -> n

It performs a case analysis to determine the least among the left and the right nodes (if the latter exists) to swap the greater root value downwards the tree. In contrast, insertion does change the shape of the tree. From |r| ≤ |l| ≤ |r| + 1 we can conclude that |l| ≤ |r| + 1 ≤ |l| + 1; this means that inserting on the right child and exchanging it with the left one maintains the Braun property. The value in the root bubbles down to its proper place to maintain the heap property:

let rec ins compare x = function
| E           -> singleton x
| N (e, l, r) ->
  if compare x e < 0
    then N (x, ins compare e r, l)
    else N (e, ins compare x r, l)

This definitely is elegant. Deletion is almost as streamlined, but with a twist:

let rec del compare = function
| E           -> failwith "empty heap"
| N (_, t, E) -> t
| N (_, E, _) -> assert false
| N (_, (N (e, _, _) as l), (N (e', _, _) as r)) ->
  if compare e e' < 0
    then N (e ,               r, del compare l)
    else N (e', rep compare e r, del compare l)

The first and second case are the trivial base cases. By the Braun property, if the right child is empty the left child can at most hold one item, and the third case exhausts the pattern asserting the impossibility. The fourth, last case is the recursive step. To eliminate a node while maintaining the Braun property we must essentially undo the swap in the insertion by deleting on the left child and exchanging it with the right one. To maintain the heap property, however, we must select the lesser of both children. If it happens to be the left one we can do just that. If it is the right one we can't blindly delete its minimum as we could end up in violation of the Braun invariant by widening the size difference to 2. This makes necessary to replace the minimum on the right with the minimum on the left, which in effect amounts to a tree rotation. With this the module is complete:

end

To build an imperative heap on Braun trees is a simple matter of wrapping the module, as I've remarked above:

module Heap : sig
  type 'a t
  val make       : ('a -> 'a -> int) -> 'a t
  val is_empty   : 'a t -> bool
  val length     : 'a t -> int
  val min        : 'a t -> 'a
  val replacemin : 'a t -> 'a -> unit
  val removemin  : 'a t -> unit
  val insert     : 'a t -> 'a -> unit
end = struct
  type 'a t = { compare : 'a -> 'a -> int; mutable root : 'a Braun.t; }
  let make compare = { compare = compare; root = Braun.empty }
  let is_empty   h   = Braun.is_empty h.root
  let length     h   = Braun.size h.root
  let min        h   = Braun.min h.root
  let replacemin h x = h.root <- Braun.rep h.compare x h.root
  let removemin  h   = h.root <- Braun.del h.compare   h.root
  let insert     h x = h.root <- Braun.ins h.compare x h.root
end

This is exactly the same interface as before. That said, I'm not very clear on the advantages of Braun trees compared to, say, binomial heaps or pairing heaps supporting merge operations in constant time. In any case, if you need the fastest priority queues and you don't mind using an imperative implementation, array-based binary heaps are the simplest, most understood data structure around. You just need a good implementation of extensible arrays.

As an aside, I have a confession to make. To reach this result I've refactored the code several times and at each step I checked that I didn't break the Braun heap invariants by running a unit test fixture against every change, with the help of a modified version of this simple test harness. I'm not a fan of TDD, as I prefer the logic and invariants to guide me in coding. The tests were a good sanity check to assert quickly if I was careful all along or if I made a mistake. In this sense it is a good safety net for exploratory programming, after the type-checker and the compiler verify that everything is right. I stress the fact that I still think that rigorous reasoning about the code comes first, types come second and exploratory programming using test is a follows from that.