2011-10-28

A First-Principles GIF Encoder

Also, quasicrystals. Have you ever found yourself wanting to do generated animations but didn't know how to visualize them? Since I couldn't find a self-contained GIF encoder, I made myself one. Even though it is not fast, I think its 220 lines are clear enough to be a good reference implementation to build upon. The encoder conforms to the following interface:


type palette
val palette : ?background:int -> color array -> palette
type t
val make    : ?ct:palette
           -> width:int -> height:int
           -> (t -> unit) -> out_channel -> unit
val image   : ?ct:palette -> ?transparent:int -> ?fps:int
           -> ?x:int -> ?y:int -> ?width:int -> ?height:int -> string
           -> t -> unit
val comment : string -> t -> unit
val repeat  : ?count:int -> t -> unit

The idea is to build the animation via repeat, comment and image inside the higher-order argument to make to ensure a well-formed output file.


GIF is a block-oriented format specifically allowing decoders to quickly skip over unknown blocks. GIF data is written as length-prefixed blocks, with the block size occupying a byte, hence a data block can be no larger than 255 bytes. Variable-length data is written as a succession of nonempty blocks terminated by an empty block, that is, a zero byte.


For this reason I need a buffering layer upon channel output. This is where I find objects a nice programming construct: manipulating stateful abstractions in an imperative way, even when the overall interface is functional. The output stream can be in immediate mode, where bytes are passed directly to the output channel, or in block mode, where bytes are accumulated in the buffer to be written in block:


class output otch = object (self)
  val         buffer  =  String.create 255
  val mutable current = -1

The member current is -1 in the first case, and in the second it marks the next empty position in buffer. When the buffer is full or at the end of the block any accumulated and unwritten bytes must be flushed out:


  method private flush =
    if current > 0 then begin
      output_byte otch current;
      output      otch buffer 0 current;
      current <- 0
    end

The output stream is byte-oriented:


  method byte n =
    let b = n land 255 in
    if current == -1 then output_byte otch b else begin
      if current == 255 then self#flush;
      buffer.[current] <- char_of_int b;
      current <- current + 1
    end

If in immediate mode the byte is written to the underlying channel, or else it is accumulated. This primitive is the building block for the other methods:


  method le16 n =
    self#byte  n;
    self#byte (n lsr 8)

  method string str =
    String.iter (fun c -> self#byte (int_of_char c)) str

(I admitted this is not particularly fast.) The switch into and out of block mode is performed via the following pair of methods:


  method begin_block =
    if current != -1 then failwith "begin_block";
    current <- 0

  method end_block =
    if current == -1 then failwith "end_block" else
    self#flush;
    current <- -1;
    self#byte 0

It doesn't make sense to nest block modes, so as a sanity check I verify before changing mode. Finally, a bit of clean-up:


  method close =
    if current != -1 then self#end_block;
    flush otch
end

Pretty simple except maybe for the fact that the output stream must be careful not to commit a fence error when writing a block of exactly 255 characters. This class is functional enough to support all the data required to build a GIF file, including the LZW-compressed image data. For that I also use a class, again to encapsulate a blob of mutable data. That doesn't mean that I will forgo functional style, though:


let fold_string f e s =
  let res = ref e in
  String.iter (fun c -> res := f !res c) s;
  !res

The LZW compressor used by GIF starts with a given code size, in principle the pixel bit-width, and expands it as it compresses data. It uses two special codes for signaling a stream reset and the end of the stream:


class lzw ?(bits=8) out =
  let count       = 1 lsl bits in
  let clear_code  = count
  and end_of_info = count + 1 in
object (self)

Since the stream is bit-oriented, the compressor keeps a bit buffer on which it appends data until there are enough bits to make whole bytes:


  val mutable code_size = 0
  val mutable buffer    = 0
  val mutable length    = 0

  method append n =
    buffer <- buffer lor ((n land pred (1 lsl code_size)) lsl length);
    length <- length + code_size;
    while length >= 8 do
      out#byte (buffer land 0xff);
      buffer <- buffer lsr 8;
      length <- length - 8
    done

Bits accumulate from least significative to most significative, that is, right to left. Of course, any left-over bits must be flushed out too:


  method finish =
    if length > 0 then
      out#byte (buffer land pred (1 lsl length))

The basic idea behind the algorithm is to build a table of those strings repeatedly seen in the input stream and to replace them with a code whenever they are detected. A good (but frustratingly sketchy in crucial parts) reference is that of Steve Blackstock, who recommends using a hash table to quickly find substrings:


  val table : (int * int, int) Hashtbl.t = Hashtbl.create 5003
  val mutable next_code = 0

The table maps substrings of the form <prefix>c to codes; the <prefix> is represented as the corresponding code, in effect treating substrings as linked lists of characters. The initial table is populated with the roots, the 2bits possible input values, plus the two special codes:


  method private reset =
    Hashtbl.clear table;
    for i = 0 to count - 1 do
      Hashtbl.add table (-1, i) i
    done;
    code_size <- bits  + 1;
    next_code <- count + 2

Since the total number of roots is 2bits + 2, the code_size is bits + 1. Whenever a new string is added to the table, the encoder must check that the code_size is sufficient, and adjust it if not. The GIF specification allows for a maximum of 12-bit codes, so if that size is reached the stream is reset:


  method private add_string s =
    if next_code == 1 lsl code_size then
      code_size <- code_size + 1;
    Hashtbl.add table s next_code;
    next_code <- next_code + 1;
    if next_code == 0x1000 then begin
      self#append clear_code;
      self#reset;
    end

The compressor follows essentially the description given in the reference:


  method compress data =
    self#reset;
    self#append clear_code;
    let last = fold_string (fun prefix c ->
      let  k = int_of_char c in
      try  Hashtbl.find table (prefix, k)
      with Not_found ->
        self#append prefix;
        self#add_string (prefix, k);
        k) (-1) data in
    self#append last;
    self#append end_of_info;
    self#finish
end

It took me a number of tries to discover the correct expand/reset sequence; all in all, I find the LZW algorithm really clever.


GIF is an indexed-color graphics format. Pixels are indices into a Color Table or palette:


type color = { red : int; green : int; blue : int; }

type palette = { background : int; bits : int; table : color array; }

Creating a palette is a matter of validating the input (GIF is restricted to power-of-two color tables) and making sure that the array of colors can't be mutated outside the palette:


let palette ?(background=0) table =
  let size = Array.length table in
  if size < 2
  || size > 256
  || size land (size - 1) != 0
  || background < 0
  || background >= size
  then invalid_arg "color_table" else
  let bits = truncate (log (float size) /. log 2.) in
  { background; bits; table = Array.copy table; }

The GIF type records the screen dimensions, the Global Color Table and the output stream:


type t = {
  width  : int;
  height : int;
  ct     : palette option;
  out    : output;
}

The GIF is built part by part, beginning with the Header and ending with the Trailer:


let header { out; _ } = out#string "GIF89a"

let trailer { out; _ } = out#byte 0x3b

(Note the concision that I gain by using a record as the data type instead of an object.) Color Tables are laid out sequentially:


let color_table { table; _ } { out; _ } =
  Array.iter (fun { red; green; blue } ->
    out#byte red  ;
    out#byte green;
    out#byte blue) table

A GIF file consists of one or more images laid out in a graphics space represented by a Logical Screen Descriptor:


let logical_screen ({ width; height; ct; out } as gif) =
  out#le16   width;
  out#le16   height;
  match ct with
  | None ->
    out#byte 0;
    out#byte 0;
    out#byte 0 (* default aspect ratio (1:1 = 49) *)
  | Some ({ background; bits; _ } as ct) ->
    out#byte (0xf0 lor pred bits);
    out#byte background;
    out#byte 0;
    color_table ct gif

If the GIF includes a global color table, it follows the Logical Screen Descriptor. In order for make to be safe it must close the output stream even in the presence of a client error:


let unwind ~(protect : 'a -> unit) f x =
  try let y = f x in protect x; y
  with e -> protect x; raise e

Then it is a matter of creating the data record, writing the header and the logical screen description, calling the client function to generate the content, and wrapping up the file with the trailer:


let make ?ct ~width ~height proc otch =
  unwind ~protect:(fun gif -> trailer gif; gif.out#close)
    (fun gif -> header gif; logical_screen gif; proc gif)
    { width; height; ct; out = new output otch; }

GIF87a already allows a GIF file to contain more than one image to be displayed in the logical screen, for instance by tiling. GIF89a includes Graphic Control Extension to make cell animations out of those images, by specifying transparency and frame delay:


let graphics_control_extension transparent fps { out; _ } =
  let delay = if fps = 0 then 0 else (200 + fps) / (fps + fps) in
  out#byte 0x21; (* GIF Extension Code *)
  out#byte 0xf9; (* Graphic Control Label *)
  out#byte 0x04; (* Length of Data Sub-Block *)
  out#byte (match transparent with None -> 0 | Some _ -> 9);
  out#le16 delay;
  out#byte (match transparent with None -> 0 | Some c -> c);
  out#byte 0x00 (* Data Sub-Block Terminator *)

An image itself is introduced by an Image Descriptor giving its position and size relative to the logical screen, and optionally a Local Color Table:


let image_descriptor ct x y width height ({ out; _ } as gif) =
  out#byte 0x2c;
  out#le16 x;
  out#le16 y;
  out#le16 width;
  out#le16 height;
  match ct with
  | None ->
    out#byte 0x00
  | Some ({ bits; _ } as ct) ->
    out#byte (0x80 lor pred bits);
    color_table ct gif

The image is laid out as a sequence of byte-sized pixels, in height rows of width bytes, represented as a string for compactness:


let image ?ct ?transparent ?(fps=0) ?(x=0) ?(y=0) ?width ?height bytes ({ out; _ } as gif) =
  let w = match width  with None -> gif.width  | Some w -> w
  and h = match height with None -> gif.height | Some h -> h in
  if x < 0 || x + w  > gif.width
  || y < 0 || x + h > gif.height
  || String.length bytes != w * h
    then invalid_arg "image";
  let bits = match ct, gif.ct with
  | Some ct, _
  | None   , Some ct -> max 2 ct.bits
  | None   , None    -> invalid_arg "image"
  in
  (match transparent, fps with
  | None, 0 -> ()
  | _   , _ -> graphics_control_extension transparent fps gif);
  image_descriptor ct x y w h gif;
  out#byte bits;
  out#begin_block;
  (new lzw ~bits out)#compress bytes;
  out#end_block

I validate the image dimensions and make sure that there is at least a color table applicable to this image to compute the pixel size in bits (GIF requires a minimum of 2-bit pixels even for bilevel images). If the image has a transparent color or an animation speed, it must be modified by a Graphics Control Extension block. Then I write the Image Descriptor and the Table-Based Image Data.


Another extension block allows the GIF file to include Comments:


let comment text { out; _ } =
  out#byte 0x21;         (* GIF Extension Code *)
  out#byte 0xfe;         (* Comment Label *)
  out#begin_block;
  out#string text;
  out#end_block

These can be useful for debugging, for instance. Looping is specified by a specific Application Extension Block introduced by Netscape:


let repeat ?(count=0) { out; _ } =
  out#byte   0x21;       (* GIF Extension code *)
  out#byte   0xff;       (* Application Extension Label *)
  out#byte   0x0b;       (* Length of application block *)
  out#string "NETSCAPE"; (* Application Identifier *)
  out#string "2.0";      (* Appl. Authentication Code *)
  out#byte   0x03;       (* Length of Data Sub-Block *)
  out#byte   0x01;       (* Loop sub-block ID *)
  out#le16   count;      (* Loop count (0 = forever) *)
  out#byte   0x00;       (* Data Sub-Block Terminator *)

That is it. As an application, let me show you the Quasicrystal generator. First I need a grayscale palette:


let grayscale = GIF.palette (Array.init 256 (fun i ->
  { red = 255 - i; green = 255 - i; blue = 255 - i }))

and a way to write to a named file:

let with_output_file proc fname =
  unwind ~protect:close_out proc (open_out_bin fname)

The code itself is a direct port from the original:


let quasicrystal ?(nwaves=7) ?(freq=27) ?(nsteps=30) ~width ~height fname =
  let pi     = 3.14159265358979312 in
  let dp     = 2. *. pi /. float nsteps
  and dt     = pi /. float nwaves
  and ds     = 1.0 /. float (max width height)
  and omega  = 2. *. pi *. float freq in
  let frame phase pixels =
    let o = ref 0 in
    let y = ref (-0.5 *. ds *. float height) in
    for i = 0 to height - 1 do
      let x = ref (-0.5 *. ds *. float width) in
      for j = 0 to width - 1 do
        let s = ref 0. in
        for k = 0 to nwaves - 1 do
          let t = dt *. float k in
          let b = !y *. cos t -. !x *. sin t in
          s := !s +. 0.5 *. (cos (b *. omega +. phase) +. 1.)
        done;
        s := !s -. 2. *. floor (0.5 *. !s);
        if !s > 1.0 then s := 2. -. !s;
        pixels.[!o] <- char_of_int (truncate (255. *. !s));
        incr o;
        x := !x +. ds
      done;
      y := !y +. ds
    done
  in
  let pixels = String.create (width * height) in
  with_output_file (GIF.make ~ct:grayscale ~width ~height (fun gif ->
    GIF.repeat gif;
    for p = 0 to nsteps - 1 do
      GIF.comment (Printf.sprintf "Frame %d/%d" (p+1) nsteps) gif;
      frame (float p *. dp) pixels;
      GIF.image ~fps:25 pixels gif;
    done) ) fname

The last 8 lines are the ones actually dealing with the GIF output. For instance, the call:


quasicrystal ~width:256 ~height:256 ~nwaves:7 ~freq:13 ~nsteps:25 "quasicrystal.gif" ;;

creates the following image:


plane-wave quasicrystal

I hope you find the code useful.

2011-09-26

Yield, Continue

Roshan P. James and Amr Sabry show in "Yield: Mainstream Delimited Continuations" the interdefinability of yield-style generators and delimited continuations. Their encoding is at the same time simple and general, and even if the examples given in the paper are in Haskell, their translation into OCaml is straightforward. So much so that the result is essentially equivalent to ASAI Kenichi's OchaCaml (Edit: this claim of mine is certainly unsubstantiated and quite possibly wrong. See Zerny's comment).


James and Sabry generalize the mechanics of yield to a three-ported construct represented by the type (ι, ο, ρ) Yield:


The type of yield

This type encapsulates the communication between an iterator and its calling context, where the iterator yields values of type ο, receives inputs of type ι and terminates (or returns) with a final result of type ρ. This communication is mediated by a delimited context that can be activated with run which … marks the boundary of an iterator and delimits the action of yield. This communication is effected by a reified continuation given by a concrete data type with which the calling context can interact:


type ('i, 'o, 'r) iterator =
| Result of 'r
| Susp of 'o * ('i -> ('i, 'o, 'r) iterator)

In effect, run converts a monadic producer that uses yield into a CPS-transformed consumer that invokes the continuation given by an iterator. These patterns of interaction can be abstracted somewhat. The most general consumer is given by foreach:


let rec foreach (m : ('i, 'o, 'r) iterator) (f : 'o -> 'i) : 'r = match m with
| Susp (v, k) -> foreach (k (f v)) f
| Result r    -> r

It applies f to each value yielded by the iterator, feeding the result back to it. If the consumer is interested just in the yielded values and not in the result of the iteration, it can fold over them:


let rec fold (m : ('i, 'i, 'r) iterator) (f : 'i -> 'j -> 'j) (e : 'j) : 'j = match m with
| Susp (v, k) -> f v (fold (k v) f e)
| Result _    -> e

The essence of the iterator is given by an abstract signature:


module type YIELD = sig
  type ('i, 'o, 'r) yield
  val return : 'r -> ('i, 'o, 'r) yield
  val (>>=) : ('i, 'o, 'r) yield -> ('r -> ('i, 'o, 's) yield) -> ('i, 'o, 's) yield
  val yield : 'o -> ('i, 'o, 'i) yield
  val run : ('i, 'o, 'r) yield -> ('i, 'o, 'r) iterator
end

which gives a multi-parameter monad together with a pair of operations: yield, that returns a computation returning the yielded value (note the difference with return); and run, that captures the computation's context and reifies it into an iterator. The paper gives two possible implementations. The first "grabs" each invocation frame turning it directly into an iterator:


module FrameGrabbingYield : YIELD = struct
  type ('i, 'o, 'r) yield = ('i, 'o, 'r) iterator

  let return e = Result e

  let rec (>>=) m f = match m with
  | Result v    -> f v
  | Susp (v, k) -> Susp (v, fun x -> k x >>= f)

  let yield v = Susp (v, return)

  let run e = e
end

The seconds uses the CPS-encoded delimited continuation monad directly:


module CPSYield : YIELD = struct
  type ('i, 'o, 'r) yield =
    { cont : 'b . ('r -> ('i, 'o, 'b) iterator) -> ('i, 'o, 'b) iterator }

  let return x = { cont = fun k -> k x }

  let (>>=) m f = { cont = fun k -> m.cont (fun x -> (f x).cont k) }

  let yield v = { cont = fun k -> Susp (v, k) }

  let run e = e.cont (fun r -> Result r)
end

This is the standard CPS monad with answer type polymorphism, as given by Kiselyov. Now yield e is shift $ return . Susp e, and run e is equivalent to reset $ e >>= return . Result). This is sufficient but a bit bare-bones. Let's build from here:


module YieldT (Y : YIELD) = struct
  include Y

In the simplest case, generators simply yield successive values. The result of the computation is the value itself, that can be updated for the next cycle:


let rec repeat x = yield x >>= repeat
  let rec from i = yield i >>= fun j -> from (succ j)

Transformers are a bit more involved in that they must consume the iterator and yield new values, in effect delimiting the control of the iterator they consume. The simplest transformer is obviously map:


let rec map f y =
    let rec go = function
    | Result r    -> return r
    | Susp (v, k) -> yield (f v) >>= fun _ -> go (k v)
    in go (run y)

(Note that the monadic fmap would only act on the result and not on the generated values.) In this case, the result of the computation is the mapped value of the original iterator, that must be continued with the original value. Truncating an iterator is straightforward:


let rec take n y =
    let rec go n = function
    | Result r -> return (Some r)
    | Susp (_, _) when n = 0 -> return None
    | Susp (v, k) -> yield v >>= fun x -> go (n - 1) (k x)
    in go n (run y)

Combining two generators is also straightforward:


let zip y1 y2 =
    let rec go = function
    | Result r1, Result r2 -> return (r1, r2)
    | Susp (v1, k1), Susp (v2, k2) ->
      yield (v1, v2) >>= fun (x1, x2) -> go (k1 x1, k2 x2)
    | _ -> failwith "zip"
    in go (run y1, run y2)
end

(Iterators that return early must be dealt with in a somewhat arbitrary way.) With this it is relatively straightforward to use iterators:


let ex1 y =
  let module Y = YieldT( (val y : YIELD) ) in
  foreach Y.(run (take 10 (map succ (from 0)))) (Printf.printf "%d ")

And both implementations give equivalent results:


# let _ = ex1 (module FrameGrabbingYield : YIELD) ;;
1 2 3 4 5 6 7 8 9 10 - : 'a option = None
# let _ = ex1 (module CPSYield : YIELD) ;;
1 2 3 4 5 6 7 8 9 10 - : 'a option = None

Furthermore, Asai's examples (as given in this Reddit thread) can be easily duplicated as well:


module Tree (Y : YIELD) = struct
  type 'a tree = E | N of 'a tree * 'a * 'a tree

  open Y

  let rec depth_walk : 'a tree -> ('b, 'a, 'b tree) yield = function
  | N (l, n, r) ->
    depth_walk l >>= fun l' ->
    yield n      >>= fun n' ->
    depth_walk r >>= fun r' ->
    return (N (l', n', r'))
  | E -> return E

  let to_list t = fold (run (depth_walk t)) (fun x xs -> x :: xs) []

  let map f t = foreach (run (depth_walk t)) f

  let samefringe l r =
    let rec visit l r = match l, r with
    | Result _, Result _ -> true
    | Susp (a, ka), Susp (b, kb)
      when a = b ->
      visit (ka a) (kb b)
    | _ -> false
    in visit (run (depth_walk l)) (run (depth_walk r))

  let swap l r =
    let rec visit l r = match l, r with
    | Susp (a, ka), Susp (b, kb) -> visit (ka b) (kb a)
    | Result t1, Result t2 -> (t1, t2)
    | _ -> failwith "Unequal number of leaves"
    in visit (run (depth_walk l)) (run (depth_walk r))
end

Note that, except for the return type polymorphism, these versions are exactly the same. To prove that all works properly, here are a number of tests:


# module T = Tree(CPSYield) ;;
# open T ;;

# let t1 = N (N (E, 10, E), 20, N (E, 30, N (E, 40, E)))
and t2 = N (N (E, 10, N (E, 20, E)), 30, N (E, 40, E))
and t3 = N (N (E, 'a', N (E, 'b', E)), 'c', N (E, 'd', E))
;;

(I omit the output of these for clarity.)


# let _ = map succ t1 ;;
- : int T.tree = N (N (E, 11, E), 21, N (E, 31, N (E, 41, E)))
# let _ = to_list t1 ;;
- : int list = [10; 20; 30; 40]
# let _ = samefringe t1 t2 ;;
- : bool = true
# let _ = swap t1 t3 ;;
- : char T.tree * int T.tree =
(N (N (E, 'a', E), 'b', N (E, 'c', N (E, 'd', E))),
 N (N (E, 10, N (E, 20, E)), 30, N (E, 40, E)))

Note that in the last example the trees retain their respective shapes but interchange the values of their leaves.

2011-09-16

Higher Order Fun

(the serious title of this post is "First-Class Modules Encode Type-Safe Reification of Types" or some such) After reading Kiselyov and Yallop's "First-class modules: hidden power and tantalizing promises", I wanted to understand how first-class modules relate to higher order types. They state that [f]irst-class modules — first-class functors — permit type constructor abstraction and polymorphism. They show one such use of constructor parametricity; I wanted to try another by encoding Haskell's canonical type class. Yes, I'm talking about Monad! Here it is:

module type MONAD = sig
  type 'a t
  val return : 'a -> 'a t
  val (>>=) : 'a t -> ('a -> 'b t) -> 'b t
  val fail : 'a t
end

The idea Oleg and Jeremy outline in "Generics for the OCaml masses" is to reify a type constructor α t as a module that wraps it together with its type parameter. Since a type constructor bound in a module type (in this case α MONAD.t) "cannot escape" for soundness' sake, it must be packed and unpacked. In effect this constitutes a higher-order existential type:

module type Monad = sig
  type a
  module Repr (M : MONAD) : sig
    val extract : a M.t
  end
end

Here it is crucial to see that the extraction is entirely dependent on the actual "interpretation" (the term originally used by Oleg and Jeremy) of the type given by its definition in the module parameterizing it. Uses of these packed representations of first-class constructors must ensure that the interpretation is consistent (something that in Haskell is assured by the use of a type class parameter, as in Monad m ⇒ …). The use of a type synonym restores polymorphism:

type 'a monad = (module Monad with type a = 'a)

In order to enable concrete instances of MONAD to "display their wares" so to speak, it is convenient to let them unpack these representations by themselves:

module Make (M : MONAD) = struct
  include M
  let run (type s) (mx : s monad) : s t =
    let module MX = (val mx : Monad with type a = s) in
    let module RX = MX.Repr(M) in
    RX.extract
end

Given a concrete constructor M.t, the σ monad is let-bound as a module, its representation under M is obtained and the corresponding value extracted. The same pattern recurs in the following genuinely higher-kinded functions, of which return is the simplest:

let return : 'a . 'a -> 'a monad =
  fun (type s) x ->
  (module struct
    type a = s
    module Repr (M : MONAD) = struct
      let extract = M.return x
    end
  end : Monad with type a = s)

Note that the result type is α monad, that is, a bona-fide module. Note also that the representation of value x under the monad M is exactly M.return x. So far, no complications. Now for something completely different:

let (>>=) : 'a . 'a monad -> ('a -> 'b monad) -> 'b monad =
  fun (type s) (type t) mx f ->
  (module struct
    type a = t
    type res = t
    module Repr (M : MONAD) = struct
      let extract =
        let module MX = (val mx : Monad with type a = s) in
        let module RX = MX.Repr(M) in
        M.(RX.extract >>= fun x ->
          let my = f x in
          let module MY = (val my : Monad with type a = res) in
          let module RY = MY.Repr(M) in
          RY.extract
        )
    end
  end : Monad with type a = t)

Given mx of module type Monad with type a = σ, its representation under M is extracted and monadically bound to f, which produces a value my of module type Monad with type a = τ, under exactly the same monadic interpretation given by M. A technicality: since the modules are abstract, they are generative (every name is fresh); hence to force the sharing of type τ between my and the result value I need to rebind it with a unique name res that can be shared across scopes. Now the next thing is a bit disturbing:

let fail : 'a . 'a monad =
  fun (type s) ->
  (module struct
    type a = s
    module Repr (M : MONAD) = struct
      let extract = M.fail
    end
  end : Monad with type a = s)

Really, is that a type-level function‽ Actually no, the syntax fun (type σ) → … binds the type σ locally in the body of the function (I'd find it clearer if the syntax were let type s = … in …, provided that no constructor escapes. Maybe in OCaml 3.14?) That's it! All that remains is to write monadic functions as if they were given the necessary type class:

let liftM f mx = mx >>= fun x -> return (f x)

let (>>) mx my = mx >>= fun _ -> my

let guard p = if p then return () else fail

(exercise: define the obvious join and verify that it really has type α monad monad → α monad. Bonus points: expand its definition). The usage with concrete instances of MONAD is surprisingly straightforward. For instance, given the standard:

module OptionM = Make(struct
  type 'a t = 'a option
  let return x = Some x
  let (>>=) mx f = match mx with None -> None | Some x -> f x
  let fail = None
end)

module ListM = Make(struct
  type 'a t = 'a list
  let return x = [x]
  let (>>=) mx f = List.(concat (map f mx))
  let fail = []
end)

The following deceptively simple code just works:

# let v = liftM succ (return 3) >>= fun x -> return (succ x);;
val v : int monad = <module>
# OptionM.run v ;;
- : int OptionM.t = Some 5
# ListM.run v ;;
- : int ListM.t = [5]

The monadic computation v exists completely independently of the computational interpretation under which it is evaluated! Note that this interpretation is selected at runtime, not at compile-time; if it weren't for the differing types, the relevant run could be a higher-order parameter to a generic function.

2011-09-10

A Note on Geomancy

The Wikipedia article on Geomancy discusses the mathematical structure of the process. A passing mention to the fact that [d]ue to the mathematics of the chart, only figures that have an even number of points total can become Judges piqued my curiosity. I verified it by computer in an instant, but my curiosity was not satisfied: how can this fact be mathematically proven? It is rather simple, actually: the Judge has the parity of both Witnesses combined. The Witnesses in turn have the overall parity of the four Nieces. The Nieces have themselves the parity of the Mothers and the Daughters combined. But the Daughters have the overall parity of the Mothers, as the former are a rearrangement of the latter. Hence the parity of the Judge is twice the parity of the four Mothers combined, and thus even.

2011-09-09

4×4-Bit Matrix Transposition

This is just a quick note to document a concrete variation of the technique for transposing a bit matrix via parallel bit operations in Hacker's Delight 7–3. In this case I needed to specialize it to 4×4 matrices represented by 16-bit integers. The underlying method is to recursively transpose each of the four 2x2-bit submatrices, and then block-transpose the 2x2-block matrix. Numbering the bits from 0 to A from LSB to MSB, the procedure looks like this:



3210
7654
BA98
FEDC

6240
7351
EAC8
FBD9

C840
D951
EA62
FB73

In terms of 16-bit words, each step looks like this:


F E D CB A 9 87 6 5 43 2 1 0
F B D 9E A C 87 3 5 16 2 4 0
F B 7 3E A 6 2D 9 5 1C 8 4 0

where I have highlighted the bits that remain fixed at each step. The technique amounts to swapping by xor-ing, as in x ^= y ^= x ^= y but bit by bit. The bits that remain fixed are masked out in order to avoid being swapped. In what follows, I use - for the zero bit and x for the one bit, and reserve hexadecimal digits as variable names for each bit position. I denote xor-ing bits i and j by juxtaposition as ij. Suppose the matrix to be transposed is stored in the 16-bit variable m. Using t as a temporary, the following operations carry out the transposition as shown:


m F E D C B A 9 8 7 6 5 4 3 2 1 0
m>>3 - - - F E D C B A 9 8 7 6 5 4 3
m^(m>>3) F E DFC EBDAC9B8 A7968574 63524130
0x0A0A - - - - x - x - - - - - x - x -
t=(m^(m>>3))&0x0A0A - - - - EB -C9 - - - - - 63 -41 -
t<<3 -EB -C9 - - - - -63 -41 - - - -
t^(t<<3) -EB -C9 EB -C9 - -63 -41 63 -41 -
m^=t^(t<<3) F B D 9 E A C 8 7 3 5 1 6 2 4 0
m>>6 - - - - - - F B D 9 E A C 8 7 3
m^(m>>6) F B D 9 E AFCB8 D793E5A1 C6827430
0x00CC - - - - - - - - x x - - x x - -
t=(m^(m>>6))&0x00CC - - - - - - - - D793 - - C682 - -
t<<6 - -D793 - -C682 - - - - - - - -
t^(t<<6) - -D793 - -C682 D793 - - C682 - -
m^=t^(t<<6) F B 7 3 E A 6 2 D 9 5 1 C 8 4 0

Each xor-assignment to m completes the corresponding recursive step as outlined above. OCaml code for this is straightforward:


let btrans4x4 m =
  let t = (m lxor (m lsr 3)) land 0x0a0a in
  let m = m lxor t lxor (t lsl 3) in
  let t = (m lxor (m lsr 6)) land 0x00cc in
  m lxor t lxor (t lsl 6)