Assessing Abstractions

Back to OCaml! Catching up with StackOverflow's OCaml questions, I found an interesting one about private type abbreviations in module signatures. One thing in that conversation that struck me as odd was the assertion that the compiler optimizes single-constructor variants, thereby subsuming the semantics of Haskell's all three declarators, data, type and newtype, into one. Definitive proof would be obtainable by diving into the compiler's source. I decided instead to read the output of ocamlopt to try to catch a glimpse of this, even if I know it doesn't constitute definitive proof. What I found is both good and bad, or as we say here una de cal y una de arena.

Consider the following module definitions:

module S = struct
  type t = P of int * int
  let make x y = P (x, y)
end

module T = struct
  type t = int * int
  let make x y = x, y
end

Module S (for "single sum type") declares a single-variant type, while T (for "tuple type") declares an abbreviation. Putting them in a file and compiling it with ocamlopt -S -c generates the following listing (prelude and other specificities omitted):

S.make	T.make
_camlAbbrev__make_1033: subq $8, %rsp movq %rax, %rdi .L101: subq $24, %r15 movq _caml_young_limit@GOTPCREL(%rip), %rax cmpq (%rax), %r15 jb .L102 leaq 8(%r15), %rax movq $2048, -8(%rax) movq %rdi, (%rax) movq %rbx, 8(%rax) addq $8, %rsp ret .L102: call _caml_call_gc jmp .L101	_camlAbbrev__make_1038: subq $8, %rsp movq %rax, %rdi .L105: subq $24, %r15 movq _caml_young_limit@GOTPCREL(%rip), %rax cmpq (%rax), %r15 jb .L106 leaq 8(%r15), %rax movq $2048, -8(%rax) movq %rdi, (%rax) movq %rbx, 8(%rax) addq $8, %rsp ret .L106: call _caml_call_gc jmp .L105

(ocamlopt -version is 3.12.1). The compiler generates the exact same code in both cases. The bolded instructions build the heap-allocated 3-word pair as tag, first component, second component. In order to test modular abstraction, I coerce these implementations into the following signatures:

module A : sig
  type t = int * int
  val make : int -> int -> t
end = T

module B : sig
  type t = private int * int
  val make : int -> int -> t
end = T

module X : sig
  type t = P of int * int
  val make : int -> int -> t
end = S

module Y : sig
  type t = private P of int * int
  val make : int -> int -> t
end = S

Modules A and B coerce module T with a public and a private type definition, respectively. Modules X and Y do the same with module S. The following code shows how to use each variant of type t in a destructuring pattern match:

let a () = A.(let   (x, y) =  make 2 3               in x + y)
let b () = B.(let   (x, y) = (make 2 3 :> int * int) in x + y)
let x () = X.(let P (x, y) =  make 2 3               in x + y)
let y () = Y.(let P (x, y) =  make 2 3               in x + y)

As explained in the manual, a private type abbreviation requires an explicit coercion for the types to be compatible. The use of a single-variant type makes for not only safer coding (the intent behind Haskell's newtype), it also allows for lighter syntax. I was surprised by the assembly code generated by ocamlopt (somewhat simplified):

a	b
_camlAbbrev__a_1060: subq $8, %rsp .L109: subq $24, %r15 movq _caml_young_limit@GOTPCREL(%rip), %rax cmpq (%rax), %r15 jb .L110 leaq 8(%r15), %rax movq $2048, -8(%rax) movq $5, (%rax) movq $7, 8(%rax) movq 8(%rax), %rbx movq (%rax), %rax leaq -1(%rax, %rbx), %rax addq $8, %rsp ret .L110: call _caml_call_gc jmp .L109	_camlAbbrev__b_1063: subq $8, %rsp .L113: subq $24, %r15 movq _caml_young_limit@GOTPCREL(%rip), %rax cmpq (%rax), %r15 jb .L114 leaq 8(%r15), %rax movq $2048, -8(%rax) movq $5, (%rax) movq $7, 8(%rax) movq 8(%rax), %rbx movq (%rax), %rax leaq -1(%rax, %rbx), %rax addq $8, %rsp ret .L114: call _caml_call_gc jmp .L113
c	d
_camlAbbrev__x_1066: subq $8, %rsp .L117: subq $24, %r15 movq _caml_young_limit@GOTPCREL(%rip), %rax cmpq (%rax), %r15 jb .L118 leaq 8(%r15), %rax movq $2048, -8(%rax) movq $5, (%rax) movq $7, 8(%rax) movq 8(%rax), %rbx movq (%rax), %rax leaq -1(%rax, %rbx), %rax addq $8, %rsp ret .L118: call _caml_call_gc jmp .L117	_camlAbbrev__y_1070: subq $8, %rsp .L121: subq $24, %r15 movq _caml_young_limit@GOTPCREL(%rip), %rax cmpq (%rax), %r15 jb .L122 leaq 8(%r15), %rax movq $2048, -8(%rax) movq $5, (%rax) movq $7, 8(%rax) movq 8(%rax), %rbx movq (%rax), %rax leaq -1(%rax, %rbx), %rax addq $8, %rsp ret .L122: call _caml_call_gc jmp .L121

Identical machine code generated in all four cases, which I find very encouraging. Here the compiler inlines the tupling constructor make (the three first bolded lines; note that an integer n is tagged as 2·n + 1 so that 2 is represented by 5 and 3 by 7). What I find a bummer is that the tuple is immediately destructured (into rbx and rax, next two bolded lines) and operated on its components (last bolded line; note in passing how addition of tagged integers is done with leaq, since 2·(m + n) + 1 = (2·m + 1) + (2·n + 1) - 1), without any kind of CSE performed. This is perfectly summarized by Pascal Cuoq in a comment to another StackOverflow question:

[…] Last time I checked it didn't even optimize the allocation of a, b in match a, b with x, y -> .... If you wish to check by yourself, I found that reading the x86 assembly generated by ocamlopt -S was convenient for me because I didn't have to learn a new representation.

This, I suspect, is an artifact of inlining in the presence of code emission that must be correct in the face of separate compilation. Disturbingly enough, ocamlopt seems to inline even across module boundaries. Separating the test into an abbrev.ml implementation:

module S = struct
  type t = P of int * int
  let make x y = P (x, y)
end

module T = struct
  type t = int * int
  let make x y = x, y
end

module A = T
module B = T
module X = S
module Y = S

with interface abbrev.mli:

module A : sig type t =              int * int val make : int -> int -> t end
module B : sig type t = private      int * int val make : int -> int -> t end
module X : sig type t =         P of int * int val make : int -> int -> t end
module Y : sig type t = private P of int * int val make : int -> int -> t end

the following code in test.ml:

open Abbrev

let a () = A.(let   (x, y) =  make 2 3               in x + y)
let b () = B.(let   (x, y) = (make 2 3 :> int * int) in x + y)
let x () = X.(let P (x, y) =  make 2 3               in x + y)
let y () = Y.(let P (x, y) =  make 2 3               in x + y)

compiles to the same assembly as above! Xavier Leroy spells out the conditions for ocamlopt to do cross-module inlining; I think it obvious that these tests do fulfill them. An interesting tidbit in that message is the command-line argument -dcmm to dump a RTL-like representation of the generated code:

(function camlTest__a_1038 (param/1056: addr)
 (let match/1057 (alloc 2048 5 7)
   (+ (+ (load match/1057) (load (+a match/1057 8))) -1)))

(function camlTest__b_1041 (param/1058: addr)
 (let match/1059 (alloc 2048 5 7)
   (+ (+ (load match/1059) (load (+a match/1059 8))) -1)))

(function camlTest__x_1044 (param/1060: addr)
 (let match/1061 (alloc 2048 5 7)
   (+ (+ (load match/1061) (load (+a match/1061 8))) -1)))

(function camlTest__y_1048 (param/1062: addr)
 (let match/1063 (alloc 2048 5 7)
   (+ (+ (load match/1063) (load (+a match/1063 8))) -1)))

I think it is not very productive to worry too much about the cost of abstractions in OCaml, as it probably is offset by the naïveté of the code generator.

Eighteen Million Noises

Update: Here is the PNGPack source-code archive:

Per second. Yes, it's benchmark time. I pitted Prof. Perlin's original implementation against Stefan Gustavson's, and the results are in! I ran two sets of benchmarks, one in Java and the other in OCaml. In both cases I compared three versions of Simplex Noise: Perlin's pipelined 3D Noise, and Gustavson's 2D Noise and 3D Noise:

• noisep: Perlin's 3D Simplex Noise
• noise2: Gustavson's 2D Simplex Noise
• noise3: Gustavson's 3D Simplex Noise

The Java programs are copied directly from the linked papers (except for one modification, noted below); the OCaml versions are pixel-faithful translations. Each benchmarked function calls noise 1,000 times in a tight loop. The benchmark runs the functions three times for 10 seconds each. I ran the benchmark on a lightly loaded, cool MacBook (2 GHz Intel Core 2 Duo, 2 GB 533 Mhz DDR2 SDRAM running Mac OS X 10.6.8). By cool I mean "left to cool down"; the temperature reached 37 °C yesterday. These are the best-of-three timings for each function:

Function	Calls/s
	OCaml	Java
noisep	2,123,490	3,417,300
noise2	18,206,090	18,239,810
noise3	10,506,630	10,404,280

The timings relative to Java's noisep are:

Function	Relative speed
	OCaml	Java
noisep	0.6214	1.0000
noise2	5.3276	5.3375
noise3	3.0745	3.0446

The big difference between OCaml's and Java's noisep is attributable to the fact that my translation is more structured than Perlin's version; for the other two functions the difference is within 1% (which I find quite amazing in and of itself).

In the case of Gustavson's Noise, I modified the code to randomly choose among 16 gradients on the unit cell, not 12 as in the code given in the paper. In contrast to Perlin's suggestion of duplicating existing vectors, I added four cell vertices so that the 2D projection of the gradient vectors onto the XY plane is unbiased. This allows indexing the gradient vector with a bit-mask instead of a modulo operation:

let grad3 = [|
( 1., 1., 0.);(-1., 1., 0.);( 1.,-1., 0.);(-1.,-1., 0.);
( 1., 0., 1.);(-1., 0., 1.);( 1., 0.,-1.);(-1., 0.,-1.);
( 0., 1., 1.);( 0.,-1., 1.);( 0., 1.,-1.);( 0.,-1.,-1.);
( 1., 1., 1.);(-1., 1., 1.);( 1.,-1.,-1.);(-1.,-1., 1.);
|]

These are the raw timings, for OCaml:

noisep: 10.45 WALL (10.44 usr +  0.01 sys = 10.45 CPU) @ 2122.11/s (n=22184)
        10.45 WALL (10.44 usr +  0.01 sys = 10.45 CPU) @ 2123.49/s (n=22184)
        10.45 WALL (10.44 usr +  0.01 sys = 10.45 CPU) @ 2122.31/s (n=22184)
noise2: 10.51 WALL (10.50 usr +  0.01 sys = 10.51 CPU) @ 18198.59/s (n=191223)
        10.50 WALL (10.50 usr +  0.01 sys = 10.50 CPU) @ 18206.09/s (n=191223)
        10.50 WALL (10.50 usr +  0.01 sys = 10.50 CPU) @ 18204.68/s (n=191223)
noise3: 10.51 WALL (10.51 usr +  0.01 sys = 10.51 CPU) @ 10504.05/s (n=110443)
        10.51 WALL (10.51 usr +  0.01 sys = 10.51 CPU) @ 10506.63/s (n=110443)
        10.57 WALL (10.56 usr +  0.01 sys = 10.57 CPU) @ 10449.81/s (n=110443)

          Rate     noisep noise3 noise2
noisep  2123+- 1/s     --   -80%   -88%
noise3 10487+-25/s   394%     --   -42%
noise2 18203+- 3/s   758%    74%     --

and for Java:

noisep: 10.04 WALL @ 3386.12/s (n=34000)
noisep: 10.24 WALL @ 3417.30/s (n=35000)
noisep: 10.24 WALL @ 3416.97/s (n=35000)
noise2: 10.01 WALL @ 18181.82/s (n=182000)
noise2: 10.00 WALL @ 18100.00/s (n=181000)
noise2: 10.03 WALL @ 18239.81/s (n=183000)
noise3: 10.04 WALL @ 10357.53/s (n=104000)
noise3: 10.09 WALL @ 10404.28/s (n=105000)
noise3: 10.01 WALL @ 10390.65/s (n=104000)

I used the following benchmark harness in OCaml:

let rand w = w *. (Random.float 1. -. 0.5)

let max_iter = 1_000
let width    = 4.

let noisep () =
  let p = (rand width, rand width, rand width) in
  for i = 1 to max_iter do ignore (Noise.noise p) done

let noise2 () =
  let p = (rand width, rand width) in
  for i = 1 to max_iter do ignore (Noise.noise2 p) done

let noise3 () =
  let p = (rand width, rand width, rand width) in
  for i = 1 to max_iter do ignore (Noise.noise3 p) done

let () =
  let open Benchmark in
  let bm = throughputN ~style:Auto ~repeat:3 10 [
    ("noisep", noisep, ());
    ("noise2", noise2, ());
    ("noise3", noise3, ());
  ] in
    print_newline();
  tabulate bm;

It uses Christophe Troestler's ocaml-benchmark module. The Java harness is similar but hand-coded:

public class NoiseBench {
  private static final double WIDTH = 4.0;
  private static final int NITER = 1000;
  private static final int NTEST = 3;
  private static final int NSECS = 10;

  private static void benchmark(String name, Runnable test) {
    for (int i = 0; i != NTEST; i++) {
      long now = System.currentTimeMillis(), lap = 0;
      int count = 0;
      while (lap < NSECS * 1000) {
        for (int j = 0; j != NITER; j++)
          test.run();
        lap = System.currentTimeMillis() - now;
        count += NITER;
      }
      System.out.printf("%s: %.2f WALL @ %.2f/s (n=%d)\n",
        name,
        (double) lap / 1000.0,
        (double) count * 1000. / (double) lap,
        count);
      System.out.flush();
    }
  }

  private static double rand() {
    return WIDTH * (Math.random() - 0.5);
  }

  private static void noisep() {
    double x = rand(), y = rand(), z = rand();
    for (int i = 0; i != NITER; i++)
      Noise3.noise(x, y, z);
  }

  private static void noise2() {
    double x = rand(), y = rand();
    for (int i = 0; i != NITER; i++)
      SimplexNoise.noise(x, y);
  }

  private static void noise3() {
    double x = rand(), y = rand(), z = rand();
    for (int i = 0; i != NITER; i++)
      SimplexNoise.noise(x, y, z);
  }

  public static void main(String[] args) {
    benchmark("noisep", new Runnable() { public void run() { noisep(); } });
    benchmark("noise2", new Runnable() { public void run() { noise2(); } });
    benchmark("noise3", new Runnable() { public void run() { noise3(); } });
  }
}

Is Immutable Data Slow?

Is immutatble data always slower than mutable state? I was challenged in a discussion over at Reddit about this issue, and I decided to perform a little experiment to back up a rough estimate with some numbers. I was surprised and delighted to see that I could make my case very strong indeed by showing that no, immutable data can be twice as fast as mutable data in a modern functional language.

I concocted a trivial synthetic benchmark using simple records intended to emulate typical business entities:

module P0 = struct
  type person = { name: string; age: int; balance: float }
end

module P1 = struct
  type person = { name: string; mutable age: int; mutable balance: float }
end

let test0 =
  let module E = struct
    open P0
    let rec test p iters =
      if iters == 0 then p.balance else
      test { p with balance = p.balance +. 1. } (iters - 1)
  end in E.test { P0.name="John Doe"; P0.age=19; P0.balance=100.}

let test1 =
  let module E = struct
    open P1
    let test p iters =
      for i = 1 to iters do
        p.balance <- p.balance +. 1.
      done
  end in E.test { P1.name="John Doe"; P1.age=19; P1.balance=100.}

The module P0 encapsulates an immutable business record; test0 is the attendant test that uses functional record updating and recursion in a classicaly pure idiom. The module P1, on the other hand, with its corresponding test test1 show the typical imperative use of a mutable record. The test and timing harness for both functions is simple enough:

let test times f iters =
  let best = ref infinity in
  for i = 1 to times do
    let t = Sys.time () in
    let _ = f iters in
    let t = Sys.time () -. t in
    if !best > t then
      best := t
  done;
  float iters /. !best

let () =
  let iters = 1_000_000 in
  Printf.printf "Immutable record: %f updates/s\n"
    (test 5 test0 iters);
  Printf.printf "Mutable record:   %f updates/s\n"
    (test 5 test1 iters)

I tested this code with the OCaml bytecode compiler on a 1.6/2.0 GHz Intel Atom computer with Windows XP, and compiled with both the bytecode and native compilers, without optimization, on a 2.2 GHz Core2 Duo MacBook Pro computer with Mac OS X 10.5.7. The results are as follows:

System	Immutable (upd/s)	Mutable (upd/s)	Speedup (Δ%)
1.6 GHz Atom/Win XP/Bytecode	3,048,780	3,558,719	16.73
2.0 GHz Atom/Win XP/Bytecode	4,000,000	4,566,210	14.16
2.2 GHz Core2 Duo/Mac OS/Bytecode	17,750,324	20,418,581	15.03
2.2 GHz Core2 Duo/Mac OS/Native	128,882,588	58,400,981	-54.69

The speedup afforded by imperative mutation in bytecode is modest, in all cases around 15%. On the other hand, the penalty incurred by mutation under the native compiler is more than double.

Of course, my correspondent was right in that the situation with current object-oriented languages is very different. I translated the benchmark straightforwardly into Java, and the result is what "common sense" would let us expect: Java is heavily oriented to mutation, or rather, it penalizes immutability quite heavily. I tried both client and server VMs, under the Windows environment detailed above:

System	Immutable (upd/s)	Mutable (upd/s)	Speedup (×)
1.6 GHz Atom/Win XP/Client	22,831,050	107,526,882	4.71
1.6 GHz Atom/Win XP/Server	37,735,849	322,580,645	8.55

The difference is dramatic: mutation is between 4.7 and 8.5 times faster than functional update in Java, and even twice as fast than native OCaml code on a faster machine. Moral of the story: test your assumptions with experiments. The benchmark code is:

public final class MutTest {
  private static final class Person0 {
    public final String name;
    public final int age;
    public final double balance;
    
    public Person0(String name, int age, double balance) {
      this.name = name;
      this.age = age;
      this.balance = balance;
    }
    
    public Person0 deposit(double amount) {
      return new Person0(this.name, this.age, this.balance + amount);
    }
  }
  
  private static final class Person1 {
    public String name;
    public int age;
    public double balance;
    
    public Person1(String name, int age, double balance) {
      this.name = name;
      this.age = age;
      this.balance = balance;
    }
    
    public void deposit(double amount) {
      this.balance += amount;
    }  
  }
  
  private interface Test {
    double test(int iters);
  }

  private static final Test TEST0 = new Test() {
    public double test(int iters) {
      Person0 person = new Person0("John Doe", 19, 100.0);
      for (int i = 0; i != iters; i++) {
        person = person.deposit(1.0);
      }
      return person.balance;
    }
  };
  
  private static final Test TEST1 = new Test() {
    public double test(int iters) {
      Person1 person = new Person1("John Doe", 19, 100.0);
      for (int i = 0; i != iters; i++) {
        person.deposit(1.0);
      }
      return person.balance;
    }
  };
  
  private static double test(int times, Test test, int iters) {
    long best = Long.MAX_VALUE;
    double balance;
    for (int i = 0; i != times; i++) {
      long now = System.currentTimeMillis();
      balance = test.test(iters);
      now = System.currentTimeMillis() - now;
      if (best > now)
        best = now;
    }
    return (double) iters / ((double) best / 1000.0);
  }

  public static void main(String[] args) {
    final int iters = 10000000;
    System.out.printf("Immutable record: %f updates/s\n", test(5, TEST0, iters));
    System.out.printf("Mutable record:  %f updates/s\n", test(5, TEST1, iters));
  }
}