main is usually a function

Continuations in C++ with fork

2012-02-14T01:26:00.000-08:00

While reading "Continuations in C" I came across an intriguing idea:

It is possible to simulate call/cc, or something like it, on Unix systems with system calls like fork() that literally duplicate the running process.

The author sets this idea aside, and instead discusses some code that uses setjmp/longjmp and stack copying. And there are several other continuation-like constructs available for C, such as POSIX getcontext. But the idea of implementing call/cc with fork stuck with me, if only for its amusement value. I'd seen fork used for computing with probability distributions, but I couldn't find an implementation of call/cc itself. So I decided to give it a shot, using my favorite esolang, C++.

Continuations are a famously mind-bending idea, and this article doesn't totally explain what they are or what they're good for. If you aren't familiar with continuations, you might catch on from the examples, or you might want to consult another source first (1, 2, 3, 4, 5, 6).

Small examples

I'll get to the implementation later, but right now let's see what these fork-based continuations can do. The interface looks like this.

template <typename T>
class cont {
public:
    void operator()(const T &x);
};

template <typename T>
T call_cc( std::function< T (cont<T>) > f );

std::function is a wrapper that can hold function-like values, such as function objects or C-style function pointers. So call_cc<T> will accept any function-like value that takes an argument of type cont<T> and returns a value of type T. This wrapper is the first of several C++11 features we'll use.

call_cc stands for "call with current continuation", and that's exactly what it does. call_cc(f) will call f, and return whatever f returns. The interesting part is that it passes to f an instance of our cont class, which represents all the stuff that's going to happen in the program after f returns. That cont object overloads operator() and so can be called like a function. If it's called with some argument x, the program behaves as though f had returned x.

The types reflect this usage. The type parameter T in cont<T> is the return type of the function passed to call_cc. It's also the type of values accepted by cont<T>::operator().

Here's a small example.

int f(cont<int> k) {
    std::cout << "f called" << std::endl;
    k(1);
    std::cout << "k returns" << std::endl;
    return 0;
}

int main() {
    std::cout << "f returns " << call_cc<int>(f) << std::endl;
}

When we run this code we get:

f called
f returns 1

We don't see the "k returns" message. Instead, calling k(1) bails out of f early, and forces it to return 1. This would happen even if we passed k to some deeply nested function call, and invoked it there.

This nonlocal return is kind of like throwing an exception, and is not that surprising. More exciting things happen if a continuation outlives the function call it came from.

boost::optional< cont<int> > global_k;

int g(cont<int> k) {
    std::cout << "g called" << std::endl;
    global_k = k;
    return 0;
}

int main() {
    std::cout << "g returns " << call_cc<int>(g) << std::endl;

    if (global_k)
        (*global_k)(1);
}

When we run this, we get:

g called
g returns 0
g returns 1

g is called once, and returns twice! When called, g saves the current continuation in a global variable. After g returns, main calls that continuation, and g returns again with a different value.

What value should global_k have before g is called? There's no such thing as a "default" or "uninitialized" cont<T>. We solve this problem by wrapping it with boost::optional. We use the resulting object much like a pointer, checking for "null" and then dereferencing. The difference is that boost::optional manages storage for the underlying value, if any.

Why isn't this code an infinite loop? Because invoking a cont<T> also resets global state to the values it had when the continuation was captured. The second time g returns, global_k has been reset to the "null" optional value. This is unlike Scheme's call/cc and most other continuation systems. It turns out to be a serious limitation, though it's sometimes convenient. The reason for this behavior is that invoking a continuation is implemented as a transfer of control to another process. More on that later.

Backtracking

We can use continuations to implement backtracking, as found in logic programming languages. Here is a suitable interface.

bool guess();
void fail();

We will use guess as though it has a magical ability to predict the future. We assume it will only return true if doing so results in a program that never calls fail. Here is the implementation.

boost::optional< cont<bool> > checkpoint;

bool guess() {
    return call_cc<bool>( [](cont<bool> k) {
        checkpoint = k;
        return true;
    } );
}

void fail() {
    if (checkpoint) {
        (*checkpoint)(false);
    } else {
        std::cerr << "Nothing to be done." << std::endl;
        exit(1);
    }
}

guess invokes call_cc on a lambda expression, which saves the current continuation and returns true. A subsequent call to fail will invoke this continuation, retrying execution in a world where guess had returned false instead. In Scheme et al, we would store a whole stack of continuations. But invoking our cont<bool> resets global state, including the checkpoint variable itself, so we only need to explicitly track the most recent continuation.

Now we can implement the integer factoring example from "Continuations in C".

int integer(int m, int n) {
    for (int i=m; i<=n; i++) {
        if (guess())
            return i;
    }
    fail();
}

void factor(int n) {
    const int i = integer(2, 100);
    const int j = integer(2, 100);

    if (i*j != n)
        fail();

    std::cout << i << " * " << j << " = " << n << std::endl;
}

factor(n) will guess two integers, and fail if their product is not n. Calling factor(391) will produce the output

17 * 23 = 391

after a moment's delay. In fact, you might see this after your shell prompt has returned, because the output is produced by a thousand-generation descendant of the process your shell created.

Solving a maze

For a more substantial use of backtracking, let's solve a maze.

const int maze_size = 15;
char maze[] =
    "X-------------+\n"
    "        |     |\n"
    "|--+  | |   | |\n"
    "|  |  | | --+ |\n"
    "|     |     | |\n"
    "|-+---+--+- | |\n"
    "| |      |    |\n"
    "| | | ---+-+- |\n"
    "|   |      |  |\n"
    "| +-+-+--|    |\n"
    "| |   |  |--- |\n"
    "|     |       |\n"
    "|--- -+-------|\n"
    "|              \n"
    "+------------- \n";

void solve_maze() {
    int x=0, y=0;

    while ((x != maze_size-1)
        || (y != maze_size-1)) {

             if (guess()) x++;
        else if (guess()) x--;
        else if (guess()) y++;
        else              y--;

        if (  (x < 0) || (x >= maze_size) ||
              (y < 0) || (y >= maze_size)  )
            fail();

        const int i = y*(maze_size+1) + x;
        if (maze[i] != ' ')
            fail();
        maze[i] = 'X';
    }

    for (char c : maze) {
        if (c == 'X')
            std::cout << "\e[1;32mX\e[0m";
        else
            std::cout << c;
    }
}

Whether code or prose, the algorithm is pretty simple. Start at the upper-left corner. As long as we haven't reached the lower-right corner, guess a direction to move. Fail if we go off the edge, run into a wall, or find ourselves on a square we already visited.

Once we've reached the goal, we iterate over the char array and print it out with some rad ANSI color codes.

Once again, we're making good use of the fact that our continuations reset global state. That's why we see 'X' marks not on the failed detours, but only on a successful path through the maze. Here's what it looks like.


X-------------+
XXXXXXXX|     |
|--+  |X|   | |
|  |  |X| --+ |
|     |XXXXX| |
|-+---+--+-X| |
| |XXX   | XXX|
| |X|X---+-+-X|
|XXX|XXXXXX|XX|
|X+-+-+--|XXX |
|X|   |  |--- |
|XXXX |       |
|---X-+-------|
|   XXXXXXXXXXX
+-------------X

Excess backtracking

We can run both examples in a single program.

int main() {
    factor(391);
    solve_maze();
}

If we change the maze to be unsolvable, we'll get:

17 * 23 = 391
23 * 17 = 391
Nothing to be done.

Factoring 391 a different way won't change the maze layout, but the program doesn't know that. We can add a cut primitive to eliminate unwanted backtracking.

void cut() {
    checkpoint = boost::none;
}

int main() {
    factor(391);
    cut();
    solve_maze();
}

The implementation

For such a crazy idea, the code to implement call_cc with fork is actually pretty reasonable. Here's the core of it.

template <typename T>
// static
T cont<T>::call_cc(call_cc_arg f) {
    int fd[2];
    pipe(fd);
    int read_fd  = fd[0];
    int write_fd = fd[1];

    if (fork()) {
        // parent
        close(read_fd);
        return f( cont<T>(write_fd) );
    } else {
        // child
        close(write_fd);
        char buf[sizeof(T)];
        if (read(read_fd, buf, sizeof(T)) < ssize_t(sizeof(T)))
            exit(0);
        close(read_fd);
        return *reinterpret_cast<T*>(buf);
    }
}

template <typename T>
void cont<T>::impl::invoke(const T &x) {
    write(m_pipe, &x, sizeof(T));
    exit(0);
}

To capture a continuation, we fork the process. The resulting processes share a pipe which was created before the fork. The parent process will call f immediately, passing a cont<T> object that holds onto the write end of this pipe. If that continuation is invoked with some argument x, the parent process will send x down the pipe and then exit. The child process wakes up from its read call, and returns x from call_cc.

There are a few more implementation details.

If the parent process exits, it will close the write end of the pipe, and the child's read will return 0, i.e. end-of-file. This prevents a buildup of unused continuation processes. But what if the parent deletes the last copy of some cont<T>, yet keeps running? We'd like to kill the corresponding child process immediately.
This sounds like a use for a reference-counted smart pointer, but we want to hide this detail from the user. So we split off a private implementation class, cont<T>::impl, with a destructor that calls close. The user-facing class cont<T> holds a std::shared_ptr to a cont<T>::impl. And cont<T>::operator() simply calls cont<T>::impl::invoke through this pointer.
It would be nice to tell the compiler that cont<T>::operator() won't return, to avoid warnings like "control reaches end of non-void function". GCC provides the noreturn attribute for this purpose.
We want the cont<T> constructor to be private, so we had to make call_cc a static member function of that class. But the examples above use a free function call_cc<T>. It's easiest to implement the latter as a 1-line function that calls the former. The alternative is to make it a friend function of cont<T>, which requires some forward declarations and other noise.

There are a number of limitations too.

As noted, the forked child process doesn't see changes to the parent's global state. This precludes some interesting uses of continuations, like implementing coroutines. In fact, I had trouble coming up with any application other than backtracking. You could work around this limitation with shared memory, but it seemed like too much hassle.
Each captured continuation can only be invoked once. This is easiest to observe if the code using continuations also invokes fork directly. It could possibly be fixed with additional forking inside call_cc.
Calling a continuation sends the argument through a pipe using a naive byte-for-byte copy. So the argument needs to be Plain Old Data, and had better not contain pointers to anything not shared by the two processes. This means we can't send continuations through other continuations, sad to say.
I left out the error handling you would expect in serious code, because this is anything but.
Likewise, I'm assuming that a single write and read will suffice to send the value. Robust code will need to loop until completion, handle EINTR, etc. Or use some higher-level IPC mechanism.
At some size, stack-allocating the receive buffer will become a problem.
It's slow. Well, actually, I'm impressed with the speed of fork on Linux. My machine solves both backtracking problems in about a second, forking about 2000 processes along the way. You can speed it up more with static linking. But it's still far more overhead than the alternatives.

As usual, you can get the code from GitHub.

Generating random functions

2012-02-02T13:24:00.000-08:00

How can we pick a random Haskell function? Specifically, we want to write an IO action

randomFunction :: IO (Integer -> Bool)

with this behavior:

It produces a function of type Integer -> Bool.
It always produces a total function — a function which never throws an exception or enters an infinite loop.
It is equally likely to produce any such function.

This is tricky, because there are infinitely many such functions (more on that later).

In another language we might produce something which looks like a function, but actually flips a coin on each new integer input. It would use mutable state to remember previous results, so that future calls will be consistent. But the Haskell type we gave for randomFunction forbids this approach. randomFunction uses IO effects to pick a random function, but the function it picks has access to neither coin flips nor mutable state.

Alternatively, we could build a lazy infinite data structure containing all the Bool answers we need. randomFunction could generate an infinite list of random Bools, and produce a function f which indexes into that list. But this indexing will be inefficient in space and time. If the user calls (f 10000000), we'll have to run 10,000,000 steps of the pseudo-random number generator, and build 10,000,000 list elements, before we can return a single Bool result.

We can improve this considerably by using a different infinite data structure. Though our solution is pure functional code, we do end up relying on mutation — the implicit mutation by which lazy thunks become evaluated data.

The data structure

import System.Random
import Data.List ( genericIndex )

Our data structure is an infinite binary tree:

data Tree = Node Bool Tree Tree

We can interpret such a tree as a function from non-negative Integers to Bools. If the Integer argument is zero, the root node holds our Bool answer. Otherwise, we shift off the least-significant bit of the argument, and look at the left or right subtree depending on that bit.

get :: Tree -> (Integer -> Bool)
get (Node b _ _) 0 = b
get (Node _ x y) n =
    case divMod n 2 of
        (m, 0) -> get x m
        (m, _) -> get y m

Now we need to build a suitable tree, starting from a random number generator state. The standard System.Random module is not going to win any speed contests, but it does have one extremely nice property: it supports an operation

split :: StdGen -> (StdGen, StdGen)

The two generator states returned by split will (ideally) produce two independent streams of random values. We use split at each node of the infinite tree.

build :: StdGen -> Tree
build g0 =
    let (b,  g1) = random g0
        (g2, g3) = split  g1
    in  Node b (build g2) (build g3)

This is a recursive function with no base case. Conceptually, it produces an infinite tree. Operationally, it produces a single Node constructor, whose fields are lazily-deferred computations. As get explores this notional infinite tree, new Nodes are created and randomness generated on demand.

get traverses one level per bit of its input integer. So looking up the integer n involves traversing and possibly creating O(log n) nodes. This suggests good space and time efficiency, though only testing will say for sure.

Now we have all the pieces to solve the original puzzle. We build two trees, one to handle positive numbers and another for negative numbers.

randomFunction :: IO (Integer -> Bool)
randomFunction = do
    neg <- build `fmap` newStdGen
    pos <- build `fmap` newStdGen
    let f n | n < 0     = get neg (-n)
            | otherwise = get pos n
    return f

Testing

Here's some code which helps us visualize one of these functions in the vicinity of zero:

test :: (Integer -> Bool) -> IO ()
test f = putStrLn $ map (char . f) [-40..40] where
    char False = ' '
    char True  = '-'

Now we can test randomFunction in GHCi:

λ> randomFunction >>= test
---- -   ---   -    - -   - --   - - -  -- --- -- --          - -- - - --  --- --
λ> randomFunction >>= test
-   ---- - - - -  - - -- -   -     ---  --- -- - --  -  --    - -  - - -  --   - 
λ> randomFunction >>= test
- ---  - - -  --  ---         -  --  -  -    -  -  - ---- - -  ---   -     -    -

Each result from randomFunction is indeed a function: it always gives the same output for a given input. This much should be clear from the fact that we haven't used any unsafe shenanigans. But we can also demonstrate it empirically:

λ> f <- randomFunction
λ> test f
-   -----  - -   -- - -   --- --  - -   - -   - -   -- - -   ---- - - - -  - --- 
λ> test f
-   -----  - -   -- - -   --- --  - -   - -   - -   -- - -   ---- - - - -  - ---

Let's also test the speed on some very large arguments:

λ> :set +s
λ> f 10000000
True
(0.03 secs, 12648232 bytes)
λ> f (2^65536)
True
(1.10 secs, 569231584 bytes)
λ> f (2^65536)
True
(0.26 secs, 426068040 bytes)

The second call with 2^65536 is faster because the tree nodes already exist in memory. We can expect our tests to be faster yet if we compile with ghc -O rather than using GHCi's bytecode interpreter.

How many functions?

Assume we have infinite memory, so that Integers really can be unboundedly large. And let's ignore negative numbers, for simplicity. How many total functions of type Integer -> Bool are there?

Suppose we made an infinite list xs of all such functions. Now consider this definition:

diag :: [Integer -> Bool] -> (Integer -> Bool)
diag xs n = not $ genericIndex xs n n

For an argument n, diag xs looks at what the nth function of xs would return, and returns the opposite. This means the function diag xs differs from every function in our supposedly comprehensive list of functions. This contradiction shows that there are uncountably many total functions of type Integer -> Bool. It's closely related to Cantor's diagonal argument that the real numbers are uncountable.

But wait, there are only countably many Haskell programs! In fact, you can encode each one as a number. There may be uncountably many functions, but there are only a countable number of computable functions. So the proof breaks down if you restrict it to a real programming language like Haskell.

In that context, the existence of xs implies that there is some algorithm to enumerate the computable total functions. This is the assumption we ultimately contradict. The set of computable total functions is not recursively enumerable, even though it is countable. Intuitively, to produce a single element of this set, we would have to verify that the function halts on every input, which is impossible in the general case.

Now let's revisit randomFunction. Any function it produces is computable: the algorithm is a combination of the pseudo-random number procedure and our tree traversal. In this sense, randomFunction provides extremely poor randomness; it only selects values from a particular measure zero subset of its result type! But if you read the type constructor (->) as "computable function", as one should in a programming language, then randomFunction is closer to doing what it says it does.

Edit: See also Luke Palmer's recent article on this subject.

Writing kernel exploits

2012-01-28T09:02:00.000-08:00

Yesterday I gave a talk about writing kernel exploits. I've posted the slides [PDF]. Here is the original description:

Did you know that a NULL pointer can compromise your entire system? Do you know how UNIX pipes, multithreading, and an obscure network protocol from 1981 are combined to take over Linux machines today? OS kernels are full of strange and interesting vulnerabilities, thanks to the subtle nature of systems code. And the kernel's ultimate authority is the ultimate prize for an attacker.
In this talk you will learn how kernel exploits work, with detailed code examples. Compared to userspace, exploiting the kernel requires a whole different bag of tricks, and we'll cover some of the most important ones. We will focus on Linux systems and x86 hardware, though most ideas will generalize. We'll start with a few toy examples, then look at some real, high-profile Linux exploits from the past two years.
You will also see how to protect your own Linux machines against kernel exploits. We'll talk about the continual cat-and-mouse game between system administrators and those who would attack even hardened kernels.

Thanks again to SIPB for giving me a venue to talk about whatever I find interesting.

Embedding GDB breakpoints in C source code

2012-01-19T09:51:00.000-08:00

Have you ever wanted to embed GDB breakpoints in C source code?

int main() {
    printf("Hello,\n");
    EMBED_BREAKPOINT;
    printf("world!\n");
    EMBED_BREAKPOINT;
    return 0;
}

One way is to directly insert your CPU's breakpoint instruction. On x86:

#define EMBED_BREAKPOINT  asm volatile ("int3;")

There are at least two problems with this approach:

They aren't real GDB breakpoints. You can't disable them, count how many times they've been hit, etc.
If you run the program outside GDB, the breakpoint instruction will crash your process.

Here is a small hack which solves both problems:

#define EMBED_BREAKPOINT \
    asm("0:"                              \
        ".pushsection embed-breakpoints;" \
        ".quad 0b;"                       \
        ".popsection;")

We place a local label into the instruction stream, and then save its address in the embed-breakpoints linker section.

Then we need to convert these addresses into GDB breakpoint commands. I wrote a tool that does this, as a wrapper for the gdb command. Here's how it works, on our initial example:

$ gcc -g -o example example.c

$ ./gdb-with-breakpoints ./example
Reading symbols from example...done.
Breakpoint 1 at 0x4004f2: file example.c, line 8.
Breakpoint 2 at 0x4004fc: file example.c, line 10.
(gdb) run
Starting program: example 
Hello,

Breakpoint 1, main () at example.c:8
8           printf("world!\n");
(gdb) info breakpoints
Num     Type           Disp Enb Address            What
1       breakpoint     keep y   0x00000000004004f2 in main at example.c:8
        breakpoint already hit 1 time
2       breakpoint     keep y   0x00000000004004fc in main at example.c:10

If we run the program normally, or in GDB without the wrapper, the EMBED_BREAKPOINT statements do nothing. The breakpoint addresses aren't even loaded into memory, because the embed-breakpoints section is not marked as allocatable.

You can find all of the code on GitHub under a BSD license. I've done only minimal testing, but I hope it will be a useful debugging tool for someone. Let me know if you find any bugs or improvements. You can comment here, or find my email address on GitHub.

I'm not sure about the decision to write the GDB wrapper in C using BFD. I also considered Haskell and elf, or Python and the new pyelftools. One can probably do something nicer using the GDB Python API, which was added a few years ago.

This code depends on a GNU toolchain: it uses GNU C extensions, GNU assembler syntax, and BFD. The GDB wrapper uses the Linux proc filesystem, so that it can pass to GDB a temporary file which has already been unlinked. You could port it to other UNIX systems by changing the tempfile handling. It should work on a variety of CPU architectures, but I've only tested it on 32- and 64-bit x86.

Zombie 6.001 starts tomorrow!

2012-01-09T22:41:00.000-08:00

The student-run revival of MIT's famous intro CS class starts tomorrow! 6.001 and its text SICP had a singular influence on the teaching of introductions to computer science — not to be confused with intro to programming, worthwhile though that subject may be. After the unfortunate demise of 6.001 at MIT, some former TAs reanimated the class as an intense four-week experience. As their description says:

Zombie-like, 6.001 rises from the dead to threaten students again. Unlike a zombie, though, it's moving quite a bit faster than it did the first time. Like the original, don't walk into the class expecting that it will teach you Scheme; instead, it attempts to teach thought patterns for computer science, and the structure and interpretation of computer programs. Three projects will be assigned and graded. Prereq: some programming experience; high confusion threshold.

I'm helping teach it this year, and it should be a lot of fun. You can follow along online or if you're in the area, come to lectures Tuesdays and Thursdays, 19:00 to 21:00 in 32-044 (that's MIT building 32, room 044).

Propane: Functional synthesis of images and animations in Haskell

2011-12-21T16:59:00.000-08:00

I just released Propane, a Haskell libary for functional synthesis of images and animations. This is a generalization of my Repa-based quasicrystal code.

It's based on the same ideas as Pan and some other projects. An image is a function assigning a color to each point in the plane. Similarly, an animation assigns an image to each point in time. Haskell's tools for functional and declarative programming can be used directly on images and animations.

For example, you can draw a red-green gradient like so:

import Propane

main = saveImage "out.png" (Size 400 400) im where
    im (x,y) = cRGB (unbal x) (unbal y) 0

Here im is the image as a function, mapping an (x,y) coordinate to a color. unbal is a function provided by Propane, which just maps the interval [-1, 1] to [0, 1].

The source package includes an animated quasicrystal and several other examples. Propane uses Repa for data-parallel array computations. That means it automatically uses multiple CPU cores for rendering, provided the program is compiled and run with threads enabled. That said, it's not yet been optimized for speed in other ways.

This is just a toy right now, but do let me know if you come up with cool enhancements or examples!

Self-modifying code for debug tracing in quasi-C

2011-11-07T05:39:00.000-08:00

Printing a program's state as it runs is the simple but effective debugging tool of programmers everywhere. For efficiency, we usually disable the most verbose output in production. But sometimes you need to diagnose a problem in a deployed system. It would be convenient to declare "tracepoints" and enable them at runtime, like so:

tracepoint foo_entry;

int foo(int n) {
  TRACE(foo_entry, "called foo(%d)\n", n);
  // ...
}

// Called from UI, monitoring interface, etc.
void debug_foo() {
  enable(&foo_entry);
}

Here's a simple implementation of this API:

typedef int tracepoint;

#define TRACE(_point, _args...) \
  do {                          \
    if (_point) printf(_args);  \
  } while (0)

static inline void enable(tracepoint *point) {
  *point = 1;
}

Each tracepoint is simply a global variable. The construct do { ... } while (0) is a standard trick to make macro-expanded code play nicely with its surroundings. We also use GCC's syntax for macros with a variable number of arguments.

This approach does introduce a bit of overhead. One concern is that reading a global variable will cause a cache miss and will also evict a line of useful data from the cache. There's also some impact from adding a branch instruction. We'll develop a significantly more complicated implementation which avoids both of these problems.

Our new solution will be specific to x86-64 processors running Linux, though the idea can be ported to other platforms. This approach is inspired by various self-modifying-code schemes in the Linux kernel, such as ftrace, kprobes, immediate values, etc. It's mostly intended as an example of how these tricks work. The code in this article is not production-ready.

The design

Our new TRACE macro will produce code like the following pseudo-assembly:

foo:
    ...
    ; code before tracepoint
    ...
tracepoint:
    nop
after_tracepoint:
    ...
    ; rest of function
    ...
    ret

do_tracepoint:
    push args to printf
    call printf
    jmp after_tracepoint

In the common case, the tracepoint is disabled, and the overhead is only a single nop instruction. To enable the tracepoint, we replace the nop instruction in memory with jmp do_tracepoint.

The `TRACE` macro

Our nop instruction needs to be big enough that we can overwrite it with an unconditional jump. On x86-64, the standard jmp instruction has a 1-byte opcode and a 4-byte signed relative displacement, so we need a 5-byte nop. Five one-byte 0x90 instructions would work, but a single five-byte instruction will consume fewer CPU resources. Finding the best way to do nothing is actually rather difficult, but the Linux kernel has already compiled a list of favorite nops. We'll use this one:

#define NOP5 ".byte 0x0f, 0x1f, 0x44, 0x00, 0x00;"

Let's check this instruction using udcli:

$ echo 0f 1f 44 00 00 | udcli -x -64 -att
0000000000000000 0f1f440000       nop 0x0(%rax,%rax)

GCC's extended inline assembly lets us insert arbitrarily bizarre assembly code into a normal C program. We'll use the asm goto flavor, new in GCC 4.5, so that we can pass C labels into our assembly code. (The tracing use case inspired the asm goto feature, and my macro is adapted from an example in the GCC manual.)

Here's how it looks:

typedef int tracepoint;

#define TRACE(_point, _args...)          \
  do {                                   \
    asm goto (                           \
      "0: " NOP5                         \
      ".pushsection trace_table, \"a\";" \
      ".quad " #_point ", 0b, %l0;"      \
      ".popsection"                      \
      : : : : __lbl_##_point);           \
    if (0) {                             \
      __lbl_##_point: printf(_args);     \
    }                                    \
  } while (0)

We use the stringify and concat macro operators, and rely on the gluing together of adjacent string literals. A call like this:

TRACE(foo_entry, "called foo(%d)\n", n);

will produce the following code:

  do {
    asm goto (
      "0: .byte 0x0f, 0x1f, 0x44, 0x00, 0x00;"
      ".pushsection trace_table, \"a\";"
      ".quad foo_entry, 0b, %l0;"
      ".popsection"
      : : : : __lbl_foo_entry);
    if (0) {
      __lbl_foo_entry: printf("called foo(%d)\n", n);
    }
  } while (0);

Besides emitting the nop instruction, we write three 64-bit values ("quads"). They are, in order:

The address of the tracepoint variable declared by the user. We never actually read or write this variable. We're just using its address as a unique key.
The address of the nop instruction, by way of a local assembler label.
The address of the C label for our printf call, as passed to asm goto.

This is the information we need in order to patch in a jmp at runtime. The .pushsection directive makes the assembler write into the trace_table section without disrupting the normal flow of code and data. The "a" section flag marks these bytes as "allocatable", i.e. something we actually want available at runtime.

We count on GCC's optimizer to notice that the condition 0 is unlikely to be true, and therefore move the if body to the end of the function. It's still considered reachable due to the label passed to asm goto, so it will not fall victim to dead code elimination.

The linker script

We have to collect all of these trace_table records, possibly from multiple source files, and put them somewhere for use by our C code. We'll do this with the following linker script:

SECTIONS {
  trace_table : {
    trace_table_start = .;
    *(trace_table)
    trace_table_end = .;
  }
}

This concatenates all trace_table sections into a single section in the resulting binary. It also provides symbols trace_table_start and trace_table_end at the endpoints of this section.

Memory protection

Linux systems will prevent an application from overwriting its own code, for good security reasons, but we can explicitly override these permissions. Memory permissions are managed per page of memory. There's a correct way to determine the size of a page, but our code is terribly x86-specific anyway, so we'll hardcode the page size of 4096 bytes.

#define PAGE_SIZE 4096
#define PAGE_OF(_addr) ( ((uint64_t) (_addr)) & ~(PAGE_SIZE-1) )

Then we can unprotect an arbitrary region of memory by calling mprotect for the appropriate page(s):

static void unprotect(void *addr, size_t len) {
  uint64_t pg1 = PAGE_OF(addr),
           pg2 = PAGE_OF(addr + len - 1);
  if (mprotect((void *) pg1, pg2 - pg1 + PAGE_SIZE,
               PROT_READ | PROT_EXEC | PROT_WRITE)) {
    perror("mprotect");
    abort();
  }
}

We're calling mprotect on a page which was not obtained from mmap. POSIX does not define this behavior, but Linux specifically allows mprotect on any page except the vsyscall page.

Enabling a tracepoint

Now we need to implement the enable function:

void enable(tracepoint *point);

We will scan through the trace_table records looking for a matching tracepoint pointer. The C struct corresponding to a trace_table record is:

struct trace_desc {
  tracepoint *point;
  void *jump_from;
  void *jump_to;
} __attribute__((packed));

The packed attribute tells GCC not to insert any padding within or after these structs. This ensures that their layout will match the records we produced from assembly. Now we can implement a linear search through this table.

void enable(tracepoint *point) {
  extern struct trace_desc trace_table_start[], trace_table_end[];
  struct trace_desc *desc;
  for (desc = trace_table_start; desc < trace_table_end; desc++) {
    if (desc->point != point)
      continue;

    int64_t offset = (desc->jump_to - desc->jump_from) - 5;
    if ((offset > INT32_MAX) || (offset < INT32_MIN)) {
      fprintf(stderr, "offset too big: %lx\n", offset);
      abort();
    }

    int32_t offset32 = offset;
    unsigned char *dest = desc->jump_from;
    unprotect(dest, 5);
    dest[0] = 0xe9;
    memcpy(dest+1, &offset32, 4);
  }
}

We enable a tracepoint by overwriting its nop with an unconditional jump. The opcode is 0xe9. The operand is a 32-bit displacement, interpreted relative to the instruction after the jump. desc->jump_from points to the beginning of what will be the jump instruction, so we subtract 5 from the displacement. Then we unprotect memory and write the new bytes into place.

That's everything. You can grab all of this code from GitHub, including a simple test program.

Pitfalls

Where to start?

This code is extremely non-portable, relying on details of x86-64, Linux, and specific recent versions of the GNU C compiler and assembler. The idea can be ported to other platforms, with some care. For example, ARM processors require an instruction cache flush after writing to code. Linux on ARM implements the cacheflush system call for this purpose.

Our code is not thread-safe, either. If one thread reaches a nop while it is being overwritten by another thread, the result will surely be a crash or other horrible bug. The Ksplice paper [PDF] discusses how to prevent this, in the context of live-patching the Linux kernel.

Is it worth opening this can of worms in order to improve performance a little? In general, no. Obviously we'd have to measure the performance difference to be sure. But for most projects, concerns of maintainability and avoiding bugs will preclude tricky hacks like this one.

The Linux kernel is under extreme demands for both performance and flexibility. It's part of every application on a huge number of systems, so any small performance improvement has a large aggregate effect. And those systems are incredibly diverse, making it likely that someone will see a large difference. Finally, kernel development will always involve tricky low-level code as a matter of course. The infrastructure is already there to support it — both software infrastructure and knowledgeable developers.

Global locking through StablePtr

2011-11-04T04:22:00.000-07:00

I spoke before of using global locks in Haskell to protect a thread-unsafe C library. And I wrote about a GHC bug which breaks the most straightforward way to get a global lock.

My new solution is to store an MVar lock in a C global variable via StablePtr. I've implemented this, and it seems to work. I'd appreciate if people could bang on this code and report any issues.

You can get the library from Hackage or browse the source, including a test program. You can also use this code as a template for including a similar lock in your own Haskell project.

The C code

On the C side, we declare a global variable and a function to read that variable.

static void* global = 0;

void* hs_globalzmlock_get_global(void) {
    return global;
}

To avoid name clashes, I gave this function a long name based on the z-encoding of my package's name. The variable named global will not conflict with another compilation unit, because it's declared static.

Another C function will set this variable, if it was previously 0. Two threads might execute this code concurrently, so we use a GCC built-in for atomic memory access.

int hs_globalzmlock_set_global(void* new_global) {
    void* old = __sync_val_compare_and_swap(&global, 0, new_global);
    return (old == 0);
}

If old is not 0, then someone has already set global, and our assignment was dropped. We report this condition to the caller.

Foreign imports

On the Haskell side, we import these C functions.

foreign import ccall unsafe "hs_globalzmlock_get_global"
    c_get_global :: IO (Ptr ())

foreign import ccall "hs_globalzmlock_set_global"
    c_set_global :: Ptr () -> IO CInt

The unsafe import of c_get_global demands justification. This wrinkle arises from the fact that GHC runs many Haskell threads on the same OS thread. A long-running foreign call from that OS thread might block unrelated Haskell code. GHC prevents this by moving the foreign call and/or other Haskell threads to a different OS thread. This adds latency to the foreign call — about 100 nanoseconds in my tests.

In most cases a 100 ns overhead is negligible. But it matters for functions which are guaranteed to return in a very short amount of time. And blocking other Haskell threads during such a short call is fine. Marking the import unsafe tells GHC to ignore the blocking concern, and generate a direct C function call.

Our function c_get_global is a good use case for unsafe, because it simply returns a global variable. In my tests, adding unsafe decreased the overall latency of locking by about 50%. We cannot use unsafe with c_set_global because, in the worst case, GCC implements atomic operations with blocking library functions. That's okay because c_set_global will only be called a few times anyway.

The Haskell code

Now we have access to a C global of type void*, and we want to store a Haskell value of type MVar (). The StablePtr module is just what we need. A StablePtr is a reference to some Haskell expression, which can be converted to Ptr (), aka void*. There is no guarantee about this Ptr () value, except that it can be converted back to the original StablePtr.

Here's how we store an MVar:

set :: IO ()
set = do
    mv  <- newMVar ()
    ptr <- newStablePtr mv
    ret <- c_set_global (castStablePtrToPtr ptr)
    when (ret == 0) $
        freeStablePtr ptr

It's fine for two threads to enter set concurrently. In one thread, the assignment will be dropped, and c_set_global will return 0. In that case we free the unused StablePtr, and the MVar will eventually be garbage-collected. StablePtrs must be freed manually, because the GHC garbage collector can't tell if some C code has stashed away the corresponding void*.

Now we can retrieve the MVar, or create it if necessary.

get :: IO (MVar ())
get = do
    p <- c_get_global
    if p == nullPtr
        then set >> get
        else deRefStablePtr (castPtrToStablePtr p)

In the common path, we do an unsynchronized read on the global variable. Only if the variable appears to contain NULL do we allocate an MVar, perform a synchronized compare-and-swap, etc. This keeps overhead low, and makes this library suitable for fine-grained locking.

All that's left is the user-visible locking interface:

lock :: IO a -> IO a
lock act = get >>= flip withMVar (const act)

Inspecting the machine code

Just for fun, let's see how GCC implements __sync_val_compare_and_swap on the AMD64 architecture.

$ objdump -d dist/build/cbits/global.o
...
0000000000000010 <hs_globalzmlock_set_global>:
  10:   31 c0                   xor    %eax,%eax
  12:   f0 48 0f b1 3d 00 00    lock cmpxchg %rdi,0x0(%rip)
  19:   00 00
....

This lock cmpxchg is the same instruction used by the GHC runtime system for its own atomic compare-and-swap. The offset on the operand 0x0(%rip) will be relocated to point at global.

Haskell hackathon in the Boston area, January 20 to 22

2011-11-03T03:47:00.000-07:00

The global sensation that is the Haskell Hackathon is coming to the Boston area. Hac Boston will be held January 20 to 22, 2012 in Cambridge, Massachusetts. It's open to all; you do not need to be a Haskell guru to attend. All you need is a basic knowledge of Haskell, a willingness to learn, and a project you're excited to help with (or a project of your own to work on).

Spaces are filling up, so be sure to register if you plan on coming. You can also coordinate projects on the HaskellWiki.

MIT is providing space (exact room to be determined) and Capital IQ is sponsoring the event. In addition to coding, there will be food and some short talks. I'm interested in giving a ~20 minute talk of some kind, with slides also available online. What would people like to hear about?

Thunks and lazy blackholes: an introduction to GHC at runtime

2011-10-31T00:17:00.000-07:00

This article is about a GHC bug I encountered recently, but it's really an excuse to talk about some GHC internals at an intro level. (In turn, an excuse for me to learn about those internals.)

I'll assume you're familiar with the basics of Haskell and lazy evaluation.

The bug

I spoke before of using global locks in Haskell to protect a thread-unsafe C library. Unfortunately a GHC bug prevents this from working. Using unsafePerformIO at the top level of a file can result in IO that happens more than once.

Here is a simple program which illustrates the problem:

import Control.Concurrent
import Control.Monad
import System.IO.Unsafe

ioThunk :: ()
ioThunk = unsafePerformIO $ do
    me <- myThreadId
    putStrLn ("IO executed by " ++ show me)
{-# NOINLINE ioThunk #-}

main :: IO ()
main = do
    replicateM_ 100 (forkIO (print ioThunk))
    threadDelay 10000  -- wait for other threads

Let's test this, following the compiler flag recommendations for unsafePerformIO.

$ ghc -V
The Glorious Glasgow Haskell Compilation System, version 7.2.1

$ ghc -rtsopts -threaded -fno-cse -fno-full-laziness dupe.hs
[1 of 1] Compiling Main             ( dupe.hs, dupe.o )
Linking dupe ...

$ while true; do ./dupe +RTS -N | head -n 2; echo ----; done

Within a few seconds I see output like this:

----
IO executed by ThreadId 35
()
----
IO executed by ThreadId 78
IO executed by ThreadId 85
----
IO executed by ThreadId 48
()
----

In the middle run, two threads executed the IO action.

This bug was reported two weeks ago and is already fixed in GHC HEAD. I tested with GHC 7.3.20111026, aka g6f5b798, and the problem seemed to go away.

Unfortunately it will be some time before GHC 7.4 is widely deployed, so I'm thinking about workarounds for my original global locking problem. I'll probably store the lock in a C global variable via StablePtr, or failing that, implement all locking in C. But I'd appreciate any other suggestions.

The remainder of this article is an attempt to explain this GHC bug, and the fix committed by Simon Marlow. It's long because

I try not to assume you know anything about how GHC works. I don't know very much, myself.
There are various digressions.

Objects at runtime

Code produced by GHC can allocate many kinds of objects. Here are just a few:

CONSTR objects represent algebraic data constructors and their associated fields. The value (Just 'x') would be represented by a CONSTR object, holding a pointer to another object representing 'x'.
FUN objects represent functions, like the value (\x -> x+1).
THUNK objects represent computations which have not yet happened. Suppose we write:
```
let x = 2 + 2 in f x x
```
This code will construct a THUNK object for x and pass it to the code for f. Some time later, f may force evaluation of its argument, and the thunk will, in turn, invoke (+). When the thunk has finished evaluating, it is overwritten with the evaluation result. (Here, this might be an I# CONSTR holding the number 4.) If f then forces its second argument, which is also x, the work done by (+) is not repeated. This is the essence of lazy evaluation.
When a thunk is forced, it's first overwritten with a BLACKHOLE object. This BLACKHOLE is eventually replaced with the evaluation result. Therefore a BLACKHOLE represents a thunk which is currently being evaluated.
Identifying this case helps the garbage collector, and it also gives GHC its seemingly magical ability to detect some infinite loops. Forcing a BLACKHOLE indicates a computation which cannot proceed until the same computation has finished. The GHC runtime will terminate the program with a <<loop>> exception.
We can't truly update thunks in place, because the evaluation result might be larger than the space originally allocated for the thunk. So we write an indirection pointing to the evaluation result. These IND objects will later be removed by the garbage collector.

Static objects

Dynamically-allocated objects make sense for values which are created as your program runs. But the top-level declarations in a Haskell module don't need to be dynamically allocated; they already exist when your program starts up. GHC allocates these static objects in your executable's data section, the same place where C global variables live.

Consider this program:

x = Just 'x'

f (Just _) = \y -> y+1

main = print (f x 3)

Ignoring optimizations, GHC will produce code where:

x is a CONSTR_STATIC object.
f is a FUN_STATIC object. When called, f will return a dynamically-allocated FUN object representing (\y -> y+1).
main is a THUNK_STATIC object. It represents the unevaluated expression formed by applying the function print to the argument (f x 3). A static thunk is also known as a constant applicative form, or a CAF for short. Like any other thunk, a CAF may or may not get evaluated. If evaluated, it will be replaced with a black hole and eventually the evaluation result. In this example, main will be evaluated by the runtime system, in deciding what IO to perform.

Black holes and revelations

That's all fine for a single-threaded Haskell runtime, but GHC supports running many Haskell threads across multiple OS threads. This introduces some additional complications. For example, one thread might force a thunk which is currently being evaluated by another thread. The thread will find a BLACKHOLE, but terminating the program would be incorrect. Instead the BLACKHOLE puts the current Haskell thread to sleep, and wakes it up when the evaluation result is ready.

If two threads force the same thunk at the same time, they will both perform the deferred computation. We could avoid this wasted effort by writing and checking for black holes using expensive atomic memory operations. But this is a poor tradeoff; we slow down every evaluation in order to prevent a rare race condition.

As a compiler for a language with pure evaluation, GHC has the luxury of tolerating some duplicated computation. Evaluating an expression twice can't change a program's behavior. And most thunks are cheap to evaluate, hardly worth the effort of avoiding duplication. So GHC follows a "lazy black-holing" strategy.¹² Threads write black holes only when they enter the garbage collector. If a thread discovers that one of its thunks has already been claimed, it will abandon the duplicated work-in-progress. This scheme avoids large wasted computations without paying the price on small computations. You can find the gritty details within the function threadPaused, in rts/ThreadPaused.c.

`unsafe[Dupable]PerformIO`

You may remember that we started, all those many words ago, with a program that uses unsafePerformIO. This breaks the pure-evaluation property of Haskell. Repeated evaluation will affect semantics! Might lazy black-holing be the culprit in the original bug?

Naturally, the GHC developers thought about this case. Here's the implementation of unsafePerformIO:

unsafePerformIO m = unsafeDupablePerformIO (noDuplicate >> m)

noDuplicate = IO $ \s -> case noDuplicate# s of s' -> (# s', () #)

unsafeDupablePerformIO (IO m) = lazy (case m realWorld# of (# _, r #) -> r)

The core behavior is implemented by unsafeDupablePerformIO, using GHC's internal representation of IO actions (which is beyond the scope of this article, to the extent I even have a scope in mind). As the name suggests, unsafeDupablePerformIO provides no guarantee against duplicate execution. The more familiar unsafePerformIO builds this guarantee by first invoking the noDuplicate# primitive operation.

The implementation of noDuplicate#, written in GHC's Cmm intermediate language, handles a few tricky considerations. But it's basically a call to the function threadPaused, which we saw is responsible for lazy black-holing. In other words, thunks built from unsafePerformIO perform eager black-holing.

Since threadPaused has to walk the evaluation stack, unsafeDupablePerformIO might be much faster than unsafePerformIO. In practice, this will matter when performing a great number of very quick IO actions, like peeking a single byte from memory. In this case it is safe to duplicate IO, provided the buffer is unchanging. Let's measure the performance difference.

import GHC.IO
import Foreign hiding (unsafePerformIO)
import System.Random
import Criterion.Main

main = do
    let sz = 1024*1024
    buf <- mallocBytes sz
    let get    i = peekByteOff buf i :: IO Word8
        peek_d i = unsafeDupablePerformIO (get i)
        peek_n i = unsafePerformIO        (get i)
        idxes = take 1024 $ randomRs (0, sz-1) (mkStdGen 49)
    evaluate (sum idxes)  -- force idxes ahead of time
    defaultMain
        [ bench "dup"   $ nf (map peek_d) idxes
        , bench "noDup" $ nf (map peek_n) idxes ]

And the results:

$ ghc -rtsopts -threaded -O2 peek.hs && ./peek +RTS -N
...

benchmarking dup
mean: 76.42962 us, lb 75.11134 us, ub 78.18593 us, ci 0.950
std dev: 7.764123 us, lb 6.300310 us, ub 9.790345 us, ci 0.950

benchmarking noDup
mean: 142.1720 us, lb 139.7312 us, ub 145.4300 us, ci 0.950
std dev: 14.43673 us, lb 11.40254 us, ub 17.86663 us, ci 0.950

So performance-critical idempotent actions can benefit from unsafeDupablePerformIO. But most code should use the safer unsafePerformIO, as our bug reproducer does. And the noDuplicate# machinery for unsafePerformIO makes sense, so what's causing our bug?

The bug, finally

After all those details and diversions, let's go back to the fix for GHC bug #5558. The action is mostly in rts/sm/Storage.c. This file is part of GHC's storage manager, which provides services such as garbage collection.

Recall that our problematic code looked like this:

ioThunk :: ()
ioThunk = unsafePerformIO $ do ...

This is an application of the function ($) to the argument unsafePerformIO. So it's a static thunk, a CAF. Here's the old description of how CAF evaluation works, from Storage.c:

The entry code for every CAF does the following:
builds a BLACKHOLE in the heap
pushes an update frame pointing to the BLACKHOLE
calls newCaf, below
updates the CAF with a static indirection to the BLACKHOLE
Why do we build an BLACKHOLE in the heap rather than just updating the thunk directly? It's so that we only need one kind of update frame - otherwise we'd need a static version of the update frame too.

So here's the problem. Normal thunks get blackholed in place, and a thread detects duplicated evaluation by noticing that one of its thunks-in-progress became a BLACKHOLE. But static thunks — CAFs — are blackholed by indirection. Two threads might perform the above procedure concurrently, producing two different heap-allocated BLACKHOLEs, and they'd never notice.

As Simon Marlow put it:

Note [atomic CAF entry]
With THREADED_RTS, newCaf() is required to be atomic (see #5558). This is because if two threads happened to enter the same CAF simultaneously, they would create two distinct CAF_BLACKHOLEs, and so the normal threadPaused() machinery for detecting duplicate evaluation will not detect this. Hence in lockCAF() below, we atomically lock the CAF with WHITEHOLE before updating it with IND_STATIC, and return zero if another thread locked the CAF first. In the event that we lost the race, CAF entry code will re-enter the CAF and block on the other thread's CAF_BLACKHOLE.

I can't explain precisely what a WHITEHOLE means, but they're used for spin locks or wait-free synchronization in various places. For example, the MVar primitives are synchronized by the lockClosure spinlock routine, which uses WHITEHOLEs.

The fix

Here's the corrected CAF evaluation procedure:

The entry code for every CAF does the following:
builds a CAF_BLACKHOLE in the heap
calls newCaf, which atomically updates the CAF with IND_STATIC pointing to the CAF_BLACKHOLE
if newCaf returns zero, it re-enters the CAF (see Note [atomic CAF entry])
pushes an update frame pointing to the CAF_BLACKHOLE

newCAF is made atomic by introducing a new helper function, lockCAF, which is reproduced here for your viewing pleasure:

STATIC_INLINE StgWord lockCAF (StgClosure *caf, StgClosure *bh)
{
    const StgInfoTable *orig_info;

    orig_info = caf->header.info;

#ifdef THREADED_RTS
    const StgInfoTable *cur_info;

    if (orig_info == &stg_IND_STATIC_info ||
        orig_info == &stg_WHITEHOLE_info) {
        // already claimed by another thread; re-enter the CAF
        return 0;
    }

    cur_info = (const StgInfoTable *)
        cas((StgVolatilePtr)&caf->header.info,
            (StgWord)orig_info,
            (StgWord)&stg_WHITEHOLE_info);

    if (cur_info != orig_info) {
        // already claimed by another thread; re-enter the CAF
        return 0;
    }

    // successfully claimed by us; overwrite with IND_STATIC
#endif

    // For the benefit of revertCAFs(), save the original info pointer
    ((StgIndStatic *)caf)->saved_info  = orig_info;

    ((StgIndStatic*)caf)->indirectee = bh;
    write_barrier();
    SET_INFO(caf,&stg_IND_STATIC_info);

    return 1;
}

We grab the CAF's info table pointer, which tells us what kind of object it is. If it's not already claimed by another thread, we write a WHITEHOLE — but only if the CAF hasn't changed in the meantime. This step is an atomic compare-and-swap, implemented by architecture-specific code. The function cas is specified by this pseudocode:

cas(p,o,n) {
   atomically {
      r = *p;
      if (r == o) { *p = n };
      return r;
   }
}

Here's the implementation for x86, using GCC extended inline assembly:

EXTERN_INLINE StgWord
cas(StgVolatilePtr p, StgWord o, StgWord n)
{
    __asm__ __volatile__ (
          "lock\ncmpxchg %3,%1"
          :"=a"(o), "=m" (*(volatile unsigned int *)p)
          :"0" (o), "r" (n));
    return o;
}

There are some interesting variations between architectures. SPARC and x86 use single instructions, while PowerPC and ARMv6 have longer sequences. Old ARM processors require a global spinlock, which sounds painful. Who's running Haskell on ARMv5 chips?

deep breath

Thanks for reading / skimming this far! I learned a lot by writing this article, and I hope you enjoyed reading it. I'm sure I said something wrong somewhere, so please do not hesitate to correct me in the comments.

Tim Harris, Simon Marlow, and Simon Peyton Jones. Haskell on a shared-memory multiprocessor. In Haskell '05: Proceedings of the 2005 ACM SIGPLAN workshop on Haskell, pages 49–61. ↩
Simon Marlow, Simon Peyton Jones, and Satnam Singh. Runtime Support for Multicore Haskell. In ICFP'09. ↩

Quasicrystals as sums of waves in the plane

2011-10-24T19:48:00.000-07:00

On the suggestion of a friend, I rendered this animation:

This quasicrystal is full of emergent patterns, but it can be described in a simple way. Imagine that every point in the plane is shaded according to the cosine of its y coordinate. The result would look like this:

Now we can rotate this image to get other waves, like these:

Each frame of the animation is a summation of such waves at evenly-spaced rotations. The animation occurs as each wave moves forward.

I recommend viewing it up close, and then from a few feet back. There are different patterns at each spatial scale.

The code

To render this animation I wrote a Haskell program, using the Repa array library. For my purposes, the advantages of Repa are:

Immutable arrays, supporting clean, expressive code
A fast implementation, including automatic parallelization
Easy output to image files, via repa-devil

Here is a simplified (but complete!) program, which renders a single still image.

import Data.Array.Repa ( Array, DIM2, DIM3, Z(..), (:.)(..) )
import qualified Data.Array.Repa          as R
import qualified Data.Array.Repa.IO.DevIL as D

import Data.Word  ( Word8   )
import Data.Fixed ( divMod' )

For clarity, we define a few type synonyms:

type R     = Float
type R2    = (R, R)
type Angle = R

We'll convert pixel indices to coordinates in the real plane, with origin at the image center. We have to decide how many pixels to draw, and how much of the plane to show.

pixels :: Int
pixels = 800
scale :: R
scale = 128

Repa's array indices are "snoc lists" of the form (Z :. x :. y). By contrast, our planar coordinates are conventional tuples.

point :: DIM2 -> R2
point = \(Z :. x :. y) -> (adj x, adj y) where
    adj n = scale * ((2 * fromIntegral n / denom) - 1)
    denom = fromIntegral pixels - 1

A single wave is a cosine depending on x and y coordinates in some proportion, determined by the wave's orientation angle.

wave :: Angle -> R2 -> R
wave th = f where
    (cth, sth) = (cos th, sin th)
    f (x,y) = (cos (cth*x + sth*y) + 1) / 2

To combine several functions, we sum their outputs, and wrap to produce a result between 0 and 1. As n increases, (wrap n) will rise to 1, fall back to 0, rise again, and so on. sequence converts a list of functions to a function returning a list, using the monad instance for ((->) r).

combine :: [R2 -> R] -> (R2 -> R)
combine xs = wrap . sum . sequence xs where
    wrap n = case divMod' n 1 of
        (k, v) | odd k     -> 1-v
               | otherwise -> v

To draw the quasicrystal, we combine waves at 7 angles evenly spaced between 0 and π.

angles :: Int -> [Angle]
angles n = take n $ enumFromThen 0 (pi / fromIntegral n)

quasicrystal :: DIM2 -> R
quasicrystal = combine (map wave (angles 7)) . point

We convert an array of floating-point values to an image in two steps. First, we map floats in [0,1] to bytes in [0,255]. Then we copy this to every color channel. The result is a 3-dimensional array, indexed by (row, column, channel). repa-devil takes such an array and outputs a PNG image file.

toImage :: Array DIM2 R -> Array DIM3 Word8
toImage arr = R.traverse arr8 (:. 4) chans where
    arr8 = R.map (floor . (*255) . min 1 . max 0) arr
    chans _ (Z :. _ :. _ :. 3) = 255  -- alpha channel
    chans a (Z :. x :. y :. _) = a (Z :. x :. y)

main :: IO ()
main = do
    let arr = R.fromFunction (Z :. pixels :. pixels) quasicrystal
    D.runIL $ D.writeImage "out.png" (toImage arr)

Running it

Repa's array operations automatically run in parallel. We just need to enable GHC's threaded runtime.

$ ghc -O2 -rtsopts -threaded quasicrystal.lhs
$ ./quasicrystal +RTS -N
$ xview out.png

And it looks like this:

Note that repa-devil silently refuses to overwrite an existing file, so you may need to rm out.png first.

On my 6-core machine, this parallel code ran in 3.72 seconds of wall-clock time, at a CPU utilization of 474%. The same code compiled without -threaded took 14.20 seconds, so the net efficiency of parallelization is 382%. This is a good result; what's better is how little work it required on my part. Cutting a mere 10 seconds from a single run is not a big deal. But it starts to matter when rendering many frames of animation, and trying out variations on the algorithm.

As a side note, switching from Float to Double increased the run time by about 30%. I suspect this is due to increased demand for memory bandwidth and cache space.

You can grab the Literate Haskell source and try it out on your own machine. This is my first Repa program ever, so I'd much appreciate feedback on improving the code.

Be sure to check out Michael Rule's work on animating quasicrystals.

Interfacing Haskell to the Concorde solver for Traveling Salesperson Problem

2011-10-21T01:12:00.000-07:00

The Traveling Salesperson Problem (TSP) is a famous optimization problem with applications in logistics, manufacturing, and art. In its planar form, we are given a set of "cities", and we want to visit each city while minimizing the total travel distance.

Finding the shortest possible tour is NP-hard, and quickly becomes infeasible as the number of cities grows. But most applications need only a heuristically good solution: a tour which is short, if not the shortest possible. The Lin-Kernighan heuristic quickly produces such tours.

The Concorde project provides a well-regarded collection of TSP solvers. I needed TSP heuristics for a Haskell project, so I wrote a Haskell interface to Concorde's Lin-Kernighan implementation. Concorde provides a C library, but it's far from clear how to use it. Instead I chose to invoke the linkern executable as a subprocess.

The core of the Haskell interface looks like this:

tsp
  :: Config     -- provides various configurable parameters
  -> (a -> R2)  -- gives the rectangular coordinates of each point
  -> [a]        -- list of points to visit
  -> IO [a]     -- produces points permuted in tour order

tsp lets you represent the points to visit using any type you like. You just provide a function to get the coordinates of each point. The Config parameter controls various aspects of the computation, including the time/quality tradeoff. Defaults are provided, and you can override these selectively using record-update syntax. All considered it's a pretty simple interface which tries to hide the complexity of interacting with an external program.

Visualizing a tour

Here's a example program which computes a tour of 1,000 random points. We'll visualize the tour using the Diagrams library.

import Diagrams.Prelude
    ( Diagram  , Point(P), fillColor   , lineWidth
    , translate, circle  , fromVertices, lightgrey )

import Diagrams.Backend.Cairo.CmdLine
    ( Cairo, defaultMain )

import Data.Colour.SRGB         ( sRGB       )
import Data.Colour.RGBSpace     ( uncurryRGB )
import Data.Colour.RGBSpace.HSV ( hsv        )

import qualified Algorithms.Concorde.LinKern as T

import Control.Monad
import Data.Monoid
import System.Random

tsp takes a list of points and a function to extract the coordinates of a point. Our points are just the coordinates themselves, so we pass the identity function.

type R2 = (Double, Double)

findTour :: [R2] -> IO [R2]
findTour = T.tsp cfg id where
    cfg = T.defConfig { T.verbose = True }

The tour is drawn as a loop of line segements. We also shade the interior of this polygon.

diaTour :: [R2] -> Diagram Cairo R2
diaTour xs@(x:_) = sty . fromVertices $ map P (xs ++ [x]) where
    sty = fillColor lightgrey . lineWidth 10

Each point visited by the tour is drawn as a circle, with hue indicating its position in the tour.

diaPoints :: [R2] -> Diagram Cairo R2
diaPoints = mconcat . map circ . zip [0..] where
    n = fromIntegral numPoints
    circ (i,p) = translate p . fillColor color $ circle 40
        where color = uncurryRGB sRGB (hsv (360*i/n) 1 1)

Now we put it all together. Note that linkern uses Euclidean distances rounded to the nearest integer. So we need coordinates with fairly large magnitudes. Picking values between 0 and 1 won't work.

numPoints :: Int
numPoints = 1000

main :: IO ()
main = do
    let rnd = randomRIO (0,10000)
    points <- replicateM numPoints (liftM2 (,) rnd rnd)
    tour   <- findTour points
    defaultMain (diaPoints tour `mappend` diaTour tour)

We run it like so:

$ export PATH=~/concorde-031219/LINKERN:$PATH
$ runhaskell tour.lhs -o out.pdf
$ xpdf out.pdf

The computation takes about 2 seconds on my machine. And the output looks like this:

You can download this post as a Literate Haskell file and run the above program. You'll need to install the concorde and diagrams packages.

The source for the concorde Haskell package includes a more full-featured version of this example.

Safe top-level mutable variables for Haskell

2011-10-18T00:18:00.000-07:00

I uploaded the safe-globals package for Haskell, which lets you declare top-level mutable variables like so:

{-# LANGUAGE TemplateHaskell #-}
import Data.Global
import Data.IORef

declareIORef "ref"
    [t| Int |]
    [e| 3   |]

main = do
    readIORef  ref >>= print  -- prints 3
    writeIORef ref 5
    readIORef  ref >>= print  -- prints 5

This creates a module-level binding

ref :: IORef Int

which can be managed through the usual module import/export mechanism.

Why global state?

Global state is a sign of bad software design, especially in Haskell. Why would we ever need it? Suppose you're wrapping a C library which is not thread-safe. Using a (hidden!) global lock, you can expose an interface which is simple and safe. In other words, you're using global state to compensate for others using global state.¹ Another use case is generating unique identifiers to speed up comparison of values. This can be done without breaking referential transparency, but you need a source of IDs which is really and truly global.

In these situations it's typical to create global variables using a hack such as

ref :: IORef Int
ref = unsafePerformIO (newIORef 3)
{-# NOINLINE ref #-}

My library is just a set of Template Haskell macros for the same hack. If global variables are seldom needed, then what good are these macros?

Writing out the hack each time is unsafe. I might forget the NOINLINE pragma, or subvert the type system with a polymorphic reference. The safe-globals library prevents these mistakes. I'm of the opinion — and I know it's not shared by all — that even questionable techniques should be made as safe as possible. Call it "harm reduction" if you like.
In ten years, if GHC 9 requires an extra pragma for safety, then safe-globals can be updated, without changing every package that uses it. If JHC's ACIO feature is ported to GHC, then safe-globals can take advantage and get rid of the hacks entirely.

But the direct impetus to write safe-globals was the appearance of the global-variables library, which drew some attention in the Haskell community. global-variables aims to solve the same problem, using a different approach with a number of drawbacks. The rest of this article outlines some of these drawbacks.

Spooky action at a distance

Among the stated features of global-variables are

Avoid having to pass references explicitly throughout the program in order to let distant parts communicate. Enable a communication by convention scheme, where e.g. different libraries may communicate without code dependencies.

This refers to the fact that two global refs with the same name string will become entangled, no matter where in a program they were declared. This is certainly a bug, not a feature. Untracked interactions between different components are the archetypal defect in software engineering.

Neither is there a clear way for a user of global-variables to opt out of this misfeature. The best you can do is augment your names with some prefix which you hope is unique — the same non-solution used by C libraries. Haskell solves namespace problems with a module system and a package system. global-variables circumvents both.

Still, suppose that you choose "communication by convention" for your library. You'll need to manually document the name and type of every ref used by this communication, since they aren't tracked by the type system. A mismatch (as from a library upgrade) will cause silent breakage. Worse, you need to tell every library user how to initialize your library's own variables, and hope that they do it correctly. When a ref is given different initializers in different declarations, the result is indeterminate.

Type clashes

A polymorphic reference, with a type like ∀ t. IORef t, breaks the type system. You can write a value of one type and then read it with another type. So it's important for global-variables to disallow polymorphic refs. The mechanism it uses is that each declaration is implicitly a family of refs, one for each monomorphic type (via Typeable).

Now consider this reasonable-looking program:

{-# LANGUAGE NoMonomorphismRestriction #-}
import Data.Global  -- from global-variables
import Data.IORef

-- Define the famous factorial function
fact :: Integer -> Integer
fact 0 = 1
fact n = n * fact (n-1)

-- Now use it through a global ref
ref = declareIORef "ref" 0
main = do
    writeIORef ref (length "hello")
    r <- readIORef ref
    print (fact r)

This will print 1, not 120. The ref is written at type Int (the return type of length) and an implicitly different ref is read at type Integer (because of the subsequent call to fact).

You can certainly argue that top-level refs should always be declared with a monomorphic type signature. Indeed, my library enforces this. But global-variables doesn't, and can't. Making type clashes a run-time error would be a step in the right direction.

A common response is that locking should be added in C code; however, concurrent programming in C is cumbersome and dangerous. It's much easier, if a bit ugly, to implement locking on the Haskell side. You could however store an MVar lock in a C global variable via StablePtr. Has anyone done this? ↩

phosphene: Fractal video feedback as a PC master boot record

2011-10-15T15:43:00.000-07:00

A few months ago, I wrote a Julia set fractal renderer for the demo challenge hosted by io.smashthestack.org.

This program runs in 16-bit x86 real mode, without any operating system. It's formatted as a PC master boot record, which is 512 bytes long. Subtracting out space reserved for a partition table, we have only 446 bytes for code and data.

Programming in such a restricted environment is quite a challenge. It's further complicated by real mode's segmented addressing. Indexing an array bigger than 64 kB requires significant extra code — and that goes double for the video frame buffer. With two off-screen buffers and a 640 × 480 × 1 byte video mode, much of my code is devoted to segment juggling.

I spent a long time playing with code compression. In the end, I couldn't find a scheme which justifies the fixed size cost of its own decoder. It seems that 16-bit x86 machine code is actually pretty information-dense. For a bigger demo or 32-bit mode (with bigger immediate operands) I'd definitely want compression.

It's totally feasible to enter 32-bit protected mode within 446 bytes, but there's little gained by doing so. You lose easy access to the PC BIOS, which is the only thing you have that resembles an operating system or standard library.

You can browse the assembly source code or grab the MBR itself. It runs well in QEMU, with or without KVM, and I also tested it on a few real machines via USB boot. With QEMU on Linux it's as simple as

$ qemu -hda phosphene.mbr

Thanks to Michael Rule for the original idea and for tons of help with tweaking the rendering algorithm. His writeup has more information about this project.

Slides from "Why learn Haskell?"

2011-10-12T18:07:00.000-07:00

Yesterday I gave a talk on the topic of "Why learn Haskell?", and I've posted the slides [PDF]. Thanks to MIT's SIPB for organizing these talks and providing tasty snacks. Thanks also to the Boston Haskell group for lots of useful feedback towards improving my talk.

If you'd like to adapt my slides for your own venue, the source is available under a CC-By-SA license. I rendered the PDF using pandoc and Beamer.

shqq: Embedding shell commands in Haskell code

2011-10-10T13:15:00.000-07:00

Shell scripts make it easy to pass data between external commands. But shell script as a programming language lacks features like non-trivial data structures and easy, robust concurrency. These would be useful in building quick solutions to system administration and automation problems.

As others have noted,¹²³⁴⁵ Haskell is an interesting alternative for these scripting tasks. I wrote the shqq library to make it a little easier to invoke external programs from Haskell. With the sh quasiquoter, you write a shell command which embeds Haskell variables, execute it as an IO action, and get the command's standard output as a String. In other words, it's a bit like the backtick operator from Perl or Ruby.

Here's a small example:

$ ghci -XQuasiQuotes
λ> import System.ShQQ
λ> let x = "/proc/uptime" in [sh| sha1sum $x |]
"337ec3fb998fb3a4650a18e0785f0992762b3cda  /proc/uptime\n"

shqq also handles escaping for you, so that the shell will not interpret special characters from Haskell variables. You can override this behavior when desired.

The big caveat is that shqq refuses to build on GHC 7.0 or earlier, due to a Unicode-handling bug in process-1.0. You'll need GHC 7.2 or later.

Finding duplicate files

As an example, here's a program to find duplicate files in the directory tree rooted at the current working directory.

{-# LANGUAGE QuasiQuotes, TupleSections #-}
import Control.Concurrent.Spawn  -- package spawn-0.3
import Data.List.Split           -- package split-0.1
import System.ShQQ               -- package shqq-0.1
import System.Posix.Files
import qualified Data.Map as M

First, the computation itself. If we pair each file with a key, such as size or checksum, we can find the groups of (potentially) duplicate files.

dupes :: (Ord k) => [(FilePath,k)] -> [[FilePath]]
dupes = filter (not . null . drop 1) . M.elems
      . foldr (\(v,k) -> M.insertWith (++) k [v]) M.empty

We want to examine files in parallel, but the operating system will complain if we have too many open files. We limit each pass to have at most 256 tests in progress at once.

inParallel :: (a -> IO b) -> [a] -> IO [b]
inParallel f xs = do p <- pool 256; parMapIO (p . f) xs

For efficiency, we find potential duplicates by size, and then checksum only these files. We use external shell commands for checksumming as well as the initial directory traversal. At the end we print the names of duplicated files, one per line, with a blank line after each group of duplicates.

main :: IO ()
main = do
    files <- endBy "\0" `fmap` [sh| find -type f -print0 |]

    let getSize f = ((f,) . fileSize) `fmap` getFileStatus f
    sizeDupes <- dupes `fmap` inParallel getSize files

    let getSha f = ((f,) . head . words) `fmap` [sh| sha1sum $f |]
    shaDupes  <- dupes `fmap` inParallel getSha (concat sizeDupes)

    mapM_ (mapM_ putStrLn . (++[""])) shaDupes

And we run it like so:

$ ghc -O -threaded -rtsopts dupe.hs
$ ./dupe +RTS -N

I included type signatures for clarity, but you wouldn't need them in a one-off script. Not counting imports and the LANGUAGE pragma, that makes 10 lines of code total. I'm pretty happy with the expressiveness of this solution, especially the use of parallel IO for an easy speedup.

Thanks to Dylan Lukes for the original idea for this library.

debug-diff: Colorized diffs between Haskell values

2011-09-13T21:43:00.000-07:00

Today I uploaded a small library for comparing Haskell values using a textual diff. This has served me well when diagnosing failing tests, or comparing a function with a proposed replacement.

Let's try it on some disassembled x86 code:

$ ghci
λ> :m + Hdis86 Debug.Diff Data.ByteString 
λ> let f = disassemble amd64 . pack 
λ> diff (f [0x31,0xff,0x41,0x5e,0x48,0x89,0xc7]) 
        (f [0x31,0xff,0x41,0x5c,0x48,0x89,0xc7])
--- /tmp/ddiff_x15061   2011-09-13 22:52:13.911842969 -0400 
+++ /tmp/ddiff_y15061   2011-09-13 22:52:13.911842969 -0400 
@@ -1,3 +1,3 @@
 [Inst [] Ixor [Reg (Reg32 RDI), Reg (Reg32 RDI)],
- Inst [Rex] Ipop [Reg (Reg64 R14)], 
+ Inst [Rex] Ipop [Reg (Reg64 R12)], 
  Inst [Rex] Imov [Reg (Reg64 RDI), Reg (Reg64 RAX)]]

This uses groom as a pretty-printer, and shells out to colordiff by default.

Sadly, a textual diff does not always produce usable results. I'm also interested in generic diffs for algebraic data types. There are some packages for this (1, 2) but I haven't yet learned how to use them. This could also be a good application for the new generics support in GHC 7.2.

Lambda to pi

2011-09-03T15:24:00.000-07:00

I gave a talk a while back which included an interpreter for the pi-calculus, and a compiler from the lambda-calculus to it. I didn't really do justice to the material in a few slides, so here's a proper blog post.

This article is available as a Literate Haskell file, which you can load into GHCi directly.

The π-calculus

If the λ-calculus is a minimal functional language, then the π-calculus is a minimal concurrent language. There's basically only three things we can do:

Fork a concurrent thread of execution
Create message channels and send them through other message channels
Loop forever

There's a conventional syntax for the π-calculus, which I don't much care for. Since we're writing an interpreter in Haskell, we'll use Haskell datatypes as our only syntax for the π-calculus.

import Control.Concurrent
  (forkIO, Chan, newChan, readChan, writeChan)

import Control.Applicative
import Control.Monad
import Control.Monad.State
import qualified Data.Map as M

Here's the abstract syntax of the π-calculus:

type Name = String

data Pi
  = Pi :|: Pi
  | Rep Pi
  | Nil
  | New Name      Pi
  | Out Name Name Pi
  | Inp Name Name Pi
  | Embed (IO ()) Pi

Names are bound to message channels, and the only thing we can send through a channel is another channel.

newtype MuChan = Wrap (Chan MuChan)

type Env = M.Map Name MuChan

Repetition and forking are straightforward:

run :: Env -> Pi -> IO ()

run env (Rep p) = forever (run env p)

run env (a :|: b) =
  let f = forkIO . run env
  in  f a >> f b >> return ()

We also have a base-case program that does nothing:

run _ Nil = return ()

And we have three operations on channels. (New bind p) creates a new channel, binds it to the local name bind, then does p:

run env (New bind p) = do
  c <- Wrap <$> newChan
  run (M.insert bind c env) p

(Out dest msg p) sends the channel msg to the channel dest, then does p:

run env (Out dest msg p) = do
  let Wrap c = env M.! dest
  writeChan c (env M.! msg)
  run env p

(Inp src bind p) reads a channel from src, binds it to the local name bind, then does p:

run env (Inp src bind p) = do
  let Wrap c = env M.! src
  recv <- readChan c
  _ <- forkIO $ run (M.insert bind recv env) p
  return ()

Why does Inp invoke forkIO? It's so that Inp under Rep is always available to receive a message, even if a subprogram is blocking on something else. You can think of this as an identity (Rep x) ≈ (x :|: Rep x). So we fork a new thread for each "connection", like the UNIX daemons of yore. If we're not under Rep then this is unnecessary but harmless.

Embed is something I threw in so that we can observe our program's execution from its side-effects:

run env (Embed x a) = x >> run env a

Compiling from the λ-calculus

Sending channels through channels is an interesting model, but how many algorithms can we really implement this way? It turns out that the π-calculus can implement any computable function. We'll demonstrate this by writing a compiler from λ-terms to π-terms.

Here's a simple, untyped lambda calculus with effects:

data Lam
  = Lam :@: Lam
  | Var Name
  | Abs Name Lam
  | Eff (IO ()) Lam

We'll use a state monad as a source of fresh names. For simplicity, we assume that the input λ-terms do not use any names beginning with '_'.

type M a = State [Name] a

runM :: M a -> a
runM = flip evalState ['_' : show (x :: Integer) | x <- [1..]]

fresh :: M Name
fresh = state (\(x:xs) -> (x,xs))
-- With mtl version 1, use 'State' not 'state'

withFresh :: (Name -> r) -> M r
withFresh f = f <$> fresh

Our compiler will produce π-calculus code which needs to send and receive pairs of values. We can encode a pair as a channel with two values enqueued. mkInp2 and mkOut2 construct functions for receiving and sending pairs, respectively.

mkInp2 :: M (Name -> (Name, Name) -> Pi -> Pi)
mkInp2 = withFresh $ \pair from (bind1, bind2) p ->
  Inp from pair  $
  Inp pair bind1 $
  Inp pair bind2 $
  p

mkOut2 :: M (Name -> (Name, Name) -> Pi -> Pi)
mkOut2 = withFresh $ \pair dest (from1, from2) p ->
  New pair       $
  Out pair from1 $
  Out pair from2 $
  Out dest pair  $
  p

Now we come to the compiler itself. The π-term produced by compile k e should evaluate the λ-term e and send the result to the "continuation channel" named by k.

compile :: Name -> Lam -> M Pi
compile k (Var x) = return (Out k x Nil)
compile k (Eff eff a) = Embed eff <$> compile k a

We implement a function as a concurrent process which receives arguments through a channel, and sends back results. ("Lambda as a Service"?) The channel elements are pairs of (function argument, channel where result should be sent).

compile k (Abs x b) = do
  inp2 <- mkInp2
  f    <- fresh
  ret  <- fresh
  body <- compile ret b
  return $
    New f   $
    Out k f $
    Rep     $
      inp2 f (x, ret) body

Note that the compiler allocates a fresh name for ret just once, but this name is bound to a different channel by each request.

A function application compiles according to the usual call-by-value rule: evaluate the function, evaluate the argument, then apply. Actually, these steps happen in three concurrent processes, but the "apply" process demands the argument value before sending it to the function, so we end up implementing strict semantics.

compile k (f :@: x) = do
  out2 <- mkOut2
  [f_cont, f_val, x_cont, x_val] <- replicateM 4 fresh
  f_proc <- compile f_cont f
  x_proc <- compile x_cont x
  let app_proc
        = Inp f_cont f_val $
          Inp x_cont x_val $
          out2 f_val (x_val, k) Nil
  return $
    New f_cont $
    New x_cont $
    (f_proc :|: x_proc :|: app_proc)

At the top level of a program, we still need a continuation channel. But we'll never read from this channel, because we're only running the program for its effects.

lambdaToPi :: Lam -> Pi
lambdaToPi b = runM $ do
  k <- fresh
  New k <$> compile k b

Testing the compiler

We'll test the compiler by doing some arithmetic on Church numerals.

-- \m n f -> n (m f)
e_mult :: Lam
e_mult = Abs "m" . Abs "n" . Abs "f" $
  (Var "n" :@: (Var "m" :@: Var "f"))

-- \f x -> f (f (f (... x)))
e_church :: Int -> Lam
e_church n = Abs "f" . Abs "x" $
  foldr (:@:) (Var "x") (replicate n $ Var "f")

-- \n -> n (\x -> trace "S" x) (trace "0" n)
e_shownum :: Lam
e_shownum = Abs "n" $
  Var "n" :@: (Abs "x" (Eff (putChar 'S') (Var "x")))
          :@: (Eff (putChar '0') (Var "n"))

-- compute 2 times 3
e_test :: Lam
e_test = e_shownum :@: (e_mult :@: e_church 2 :@: e_church 3)

And we run it like so:

GHCi> run M.empty (lambdaToPi e_test)

GHCi> 0SSSSSS

The answer of 6 is printed to the terminal, asynchronously, in unary.

Notes

I thought up this particular scheme for compiling to the π-calculus, but it's probably been documented already.

I'm not sure my π-calculus interpreter is correct. In particular, the term (Rep (x :|: y)) is a forkbomb, when it should be operationally equivalent to (Rep x :|: Rep y). I don't know if fixing this would be easy or hard.

Edward Kmett wrote an interpreter and pretty-printer for the π-calculus, in the style of "Finally Tagless, Partially Evaluated" [PDF] by Carette, Kiselyov, and Shan.

New slides and how I made them

2011-08-31T23:47:00.000-07:00

New slides

I've posted PDFs of slides from some talks I gave a while ago.

First-Class Concurrency in Haskell covers topics in concurrent imperative programming, such as:

Using closures with concurrency primitives to improve your API
Using MVar to retrieve thread results, as in spawn
How System.Timeout works
Implementing IORef and Chan in terms of other primitives
Implementing the π-calculus, and compiling the λ-calculus to it

This turned out to be far too much material for a one-hour talk, so the slides are probably more valuable than the talk itself. I'm also thinking of expanding the π-calculus bits into a proper blog post at some point.

High-level FFI in Haskell covers tricks for making bindings to C libraries that feel more like native Haskell libraries. Most of the examples are drawn from my hdis86 library . This includes:

Translating C types to idiomatic Haskell types
Passing Haskell functions as C function pointers
Passing ByteStrings as C buffers
Automatic locking and memory management of library state
Exposing a C library's state machine as a lazy pure function
Bundling a C library with Cabal while also allowing dynamic linking

These slides were available before, in a somewhat strange HTML format. The new PDFs are more portable and I think they look nicer, too.

How I made the slides

I originally wrote the slides in Markdown format and converted them to a S5 slideshow using pandoc. This was a great way to quickly prepare slides, including syntax-highlighted Haskell code. However I ran into a few limitations:

Slide layout depends on window size and screen resolution, rendering my obsessive tweaking useless. In some cases the text would run off the edge of the screen or overlap other elements.
The slides don't work at all in Chromium, as of the S5 version I used initially. I also got reports of issues with other browsers.

I wanted to switch over to Beamer, as I'd seen some very high-quality slides produced using that package. This would also provide a PDF that renders exactly the same on every machine. Fortunately, pandoc can output L^aT_eX code suitable for Beamer. I can keep using the lightweight syntax of Markdown, and most of my old slide source code works as-is. Once again pandoc handles a tricky problem with ease.

I followed this scheme with a few modifications. For example, pdflatex was unhappy with slides ("frames") containing verbatim text environments. The solution was to declare all frames as "fragile", whatever that means:

\begin{frame}[fragile]

Recent versions of pandoc accept the --listings option and will use the L^aT_eX listings package to render nicely-formatted source code. This article was helpful for configuring listings. The default Haskell style has a rather odd definition of what counts as a keyword, so I wrote a custom language definition. I ended up with this L^aT_eX header, which you can pass to pandoc using -H:

\usepackage{listings}
\usepackage{color}

\lstdefinelanguage{Haskell_mod}{
    otherkeywords={.., ::, |, <-, ->, @, ~, =, =>},
    morekeywords={
       case, class, data, default, deriving, do, else,
       foreign, if, import, in, infix, infixl, infixr,
       instance, let, module, newtype, of, then, type,
       where, _, forall, ccall },
    sensitive,
    morecomment=[l]--,
    morecomment=[n]{\{-}{-\}},
    morestring=[b]"
}[keywords,comments,strings]

\lstset{
    language=Haskell_mod,
    basicstyle=\ttfamily,
    keywordstyle=\color[rgb]{0,0,1},
    commentstyle=\color[rgb]{0.133,0.545,0.133},
    stringstyle=\color[rgb]{0.627,0.126,0.941},
    showstringspaces=false,
    frame=single
}

The Markdown sources for the slides are available: 1, 2.

My biggest remaining complaint is that I have a lot of explicit spacing commands, e.g. \vspace{1em}. These are used for logical grouping, and so can't be fully inferred, but maybe there's a better (partial?) solution.

x86 disassembly in Haskell with hdis86

2011-03-19T15:43:00.000-07:00

The first serious projects I coded in Haskell were compilers, code analyzers, and the like. I think it's a domain that really plays to the strengths of the language. So it makes sense that we'd want a top-notch disassembler library for Haskell.

udis86 is a fast, complete, flexible disassembler for x86 and x86-64 / AMD64. It provides a clean C API, which was no trouble to import using hsc2hs. But any C API is going to feel clunky next to high-level idiomatic Haskell code. So I built two additional layers of wrapping, to make something that will feel like a natural part of your next Haskell app.

You can download my library hdis86 from Hackage, or browse the source at GitHub. By default, a copy of udis86 is embedded in hdis86. If you already have udis86 installed as a shared library, you can link against that instead.

Getting started

Let's try it out. We'll feed in a ByteString of machine code, and pretty-print the result with groom.

$ ghciλ> :m + Hdis86 Data.ByteString Text.Groomλ> let code = pack [0xcc, 0xf0, 0xff, 0x44, 0x9e, 0x0f]λ> Prelude.putStrLn . groom $ disassemble intel64 code[Inst [] Iint3 [], Inst [Lock] Iinc   [Mem      (Memory{mSize = Bits32, mBase = Reg64 RSI, mIndex = Reg64 RBX,              mScale = 4, mOffset = Immediate{iSize = Bits8, iValue = 15}})]]

If you're not familiar with x86, you might be surprised by the range of simple and complicated instructions. The first instruction is a humble int3, which executes a trap to a debugger. It takes no operands and has no prefixes.

The second instruction is an increment of a 32-bit memory region, whose address is computed by summing the value in register RSI, the value in register RBX times a scale factor of 4, and a fixed offset of 15. This instruction also has a lock prefix, which changes its concurrency semantics. Some prefixes have a direct meaning like this; others (like Rex) will change how the disassembler decodes operands.

We can also ask for Intel-style assembly syntax:

λ> let cfg = intel64 { cfgSyntax = SyntaxIntel }λ> mapM_ (Prelude.putStrLn . mdAssembly) $ disassembleMetadata cfg codeint3 lock inc dword [rsi+rbx*4+0xf]

Analyzing instructions

If we're just printing code, showing an Instruction is needlessly verbose. The real point of the Instruction type is machine-code analysis, with the help of Haskell's pattern matching capabilities.

Sometimes two fragments of code will differ only in which registers they use. Perhaps two compilers made different choices during register allocation. We'll write a program to detect this.

import Hdis86
import Text.Groom
import Control.Monad
import qualified Data.Map as M
import qualified Data.ByteString as B

When one code fragment uses register X the same way that another fragment uses register Y, we record a constraint X :-> Y:

data Constraint = Register :-> Register

We compare the two code fragments to produce a list of Constraints, or Nothing if they differ in ways beyond register selection. We'll use this helper function to check a list of Boolean conditions:

(==>) :: [Bool] -> a -> Maybe a
ps ==> x = guard (and ps) >> Just x

We compare operands pairwise. Register operands are easy:

operand :: Operand -> Operand -> Maybe [Constraint]
operand (Reg rx) (Reg ry) = Just [rx :-> ry]

For a memory operand, we need to check that the size, scale, and offset are equal:

-- size, base register, index register, scale, offset
operand (Mem (Memory sx bx ix kx ox)) (Mem (Memory sy by iy ky oy))
  = [sx == sy, kx == ky, ox == oy] ==> [bx :-> by, ix :-> iy]

Immediate operands need a similar check for equality, but they don't produce any register constraints:

operand (Ptr   px) (Ptr   py) = [px == py] ==> []
operand (Imm   ix) (Imm   iy) = [ix == ix] ==> []
operand (Jump  ix) (Jump  iy) = [ix == iy] ==> []
operand (Const ix) (Const iy) = [ix == iy] ==> []

In all other cases, we have two different operand constructors, which means they don't match:

operand _ _ = Nothing

To check a pair of instructions, we check their prefixes and opcodes, then check operands pairwise:

inst :: Instruction -> Instruction -> Maybe [Constraint]
inst (Inst px ox rx) (Inst py oy ry) = do
  guard $ and [px == py, ox == oy, length rx == length ry]
  concat `fmap` zipWithM operand rx ry

Once we have a list of Constraints, we need to check it for consistency. If register X maps to Y in one place, it can't map to Z somewhere else:

unify :: [Constraint] -> Maybe (M.Map Register Register)
unify = foldM f M.empty where
  f m (rx :-> ry) = case M.lookup rx m of
    Nothing  -> Just (M.insert rx ry m)
    Just ry' -> [ry == ry'] ==> m

Now we put it all together, checking consistency in both directions:

regMap :: [Instruction] -> [Instruction] -> Maybe (M.Map Register Register)
regMap xs ys = do
  cs <- concat `fmap` zipWithM inst xs ys
  let swap (x :-> y) = (y :-> x)
  _ <- unify $ map swap cs
  unify cs

main :: IO ()
main = putStrLn . groom $ regMap (f prog_a) (f prog_b) where
  f = disassemble intel64

We'll test these two code fragments:

prog_a, prog_b :: B.ByteString
prog_a = B.pack
  [ 0x7e, 0x3a               -- jle  0x3c
  , 0x48, 0x89, 0xf5         -- mov  rbp, rsi
  , 0xbb, 1, 0, 0, 0         -- mov  ebx, 0x1
  , 0x48, 0x8b, 0x7d, 0x08 ] -- mov  rdi, [rbp+0x8]
prog_b = B.pack
  [ 0x7e, 0x3a               -- jle  0x3c
  , 0x48, 0x89, 0xf3         -- mov  rbx, rsi
  , 0xbd, 1, 0, 0, 0         -- mov  ebp, 0x1
  , 0x48, 0x8b, 0x7b, 0x08 ] -- mov  rdi, [rbx+0x8]

And the result is:

$ runhaskell regmap.hs Just  (fromList     [(RegNone, RegNone), (Reg32 RBX, Reg32 RBP),      (Reg64 RBP, Reg64 RBX), (Reg64 RSI, Reg64 RSI),      (Reg64 RDI, Reg64 RDI)])

In reality, this analysis is grossly incomplete. Some instructions have implicit register operands, and some pairs of registers overlap, like EBX and RBX. And we have to understand contextual requirements. It's not okay to remap RAX to RBX if some calling code is expecting a return value in RAX.

The complexity of x86 (not to mention the halting problem) means that binary code analysis will never be easy. Hopefully hdis86 can be one of many useful tools in this domain. As always, suggestions or patches are welcome.

Implementing breakpoints on x86 Linux

2011-01-09T13:57:00.000-08:00

Debugger features such as breakpoints go well beyond the normal ways in which programs interact. How can one process force another process to pause at a particular line of code? We can bet that the debugger is using some special features of the CPU and/or operating system.

As I learned about implementing breakpoints, I was surprised to find that these interfaces, while arcane, are not actually very complex. In this article we'll implement a minimal but working breakpoints library in about 75 lines of C code. All of the code is shown here, and is also available on GitHub.

True, most of us will never write a debugger. But breakpoints are useful in many other tools: profilers, compatibility layers, security scanners, etc. If your program has to interact deeply with another process, without modifying its source code, breakpoints might just do the trick. And I think there's value in knowing what your debugger is really doing.

Our code will compile with gcc and will run on Linux machines with 32- or 64-bit x86 processors. I'll assume you're familiar with the basics of UNIX programming in C. We'll make heavy use of the ptrace system call, which I'll try to explain as we go. If you're looking for other reading about ptrace, I highly recommend this article and of course the manpage.

The API

Our little breakpoint library is named breakfast. Here's breakfast.h:

#ifndef _BREAKFAST_H
#define _BREAKFAST_H

#include <sys/types.h>  /* for pid_t */

typedef void *target_addr_t;
struct breakpoint;

void breakfast_attach(pid_t pid);
target_addr_t breakfast_getip(pid_t pid);
struct breakpoint *breakfast_break(pid_t pid, target_addr_t addr);
int breakfast_run(pid_t pid, struct breakpoint *bp);

#endif

Before we dive into the implementation, we'll describe the API our library provides. breakfast is used by one process (the tracing process) to place breakpoints inside the code of another process (the target process). To establish this relationship, the tracing process calls breakfast_attach with the target's process ID. This also suspends execution of the target.

We use breakfast_getip to read the target's instruction pointer, which holds the address of the next instruction it will execute. Note that this is a pointer to machine code within the target's address space. Dereferencing the pointer within our tracing process would make little sense. The type alias target_addr_t reminds us of this fact.

We use breakfast_break to create a breakpoint at a specified address in the target's machine code. This function returns a pointer to a breakpoint structure, which has unspecified contents.

Calling breakfast_run lets the target run until it hits a breakpoint or it exits; breakfast_run returns 1 or 0 respectively. On the first invocation, we'll pass NULL for the argument bp. On subsequent calls, we're resuming from a breakpoint, and we must pass the corresponding breakpoint struct pointer. The breakfast_getip function will tell us the breakpoint's address, but it's our responsibility to map this to the correct breakpoint structure.

The trap instruction

We will implement a breakpoint by inserting a new CPU instruction into the target's memory at the breakpoint address. This instruction should suspend execution of the target and give control back to the OS.

There are many ways to return control to the OS, but we'd like to minimize disruption to the code we're hot-patching. x86 provides an instruction for exactly this purpose. It's named int3 and is encoded as the single byte 0xCC. When the CPU executes int3, it will stop what it's doing and jump to the interrupt 3 service routine, which is a piece of code chosen by the OS. On Linux, this routine will send the signal SIGTRAP to the current process.

Since we're placing int3 in the target's code, the target will receive a SIGTRAP. Under normal circumstances this would invoke the target's SIGTRAP handler, which usually kills the process. Instead we'd like the tracing process to intercept this signal and interpret it as the target hitting a breakpoint. We'll accomplish this through the magic of ptrace.

Let's start off breakfast.c with some includes:

#include <sys/ptrace.h>
#include <sys/reg.h>
#include <sys/wait.h>
#include <stdlib.h>
#include "breakfast.h"

Then we'll define our trap instruction:

#if defined(__i386)
#define REGISTER_IP EIP
#define TRAP_LEN    1
#define TRAP_INST   0xCC
#define TRAP_MASK   0xFFFFFF00

#elif defined(__x86_64)
#define REGISTER_IP RIP
#define TRAP_LEN    1
#define TRAP_INST   0xCC
#define TRAP_MASK   0xFFFFFFFFFFFFFF00

#else
#error Unsupported architecture
#endif

We use some GCC pre-defined macros to discover the machine architecture, and we define a few constants accordingly. REGISTER_IP expands to a constant from sys/reg.h which identifies the machine register holding the instruction pointer.

The trap instruction is stored as an integer TRAP_INST and its length in bytes is TRAP_LEN. These are the same on 32- and 64-bit x86. The trap instruction is a single byte, but we'll read and write the target's memory in increments of one machine word, meaning 32 or 64 bits. So we'll read 4 or 8 bytes of machine code, clear out the first byte with TRAP_MASK, and substitute 0xCC. Since x86 is a little-endian architecture, the first byte in memory is the least-significant byte of an integer machine word.

Invoking `ptrace`

All of the various ptrace requests are made through a single system call named ptrace. The first argument specifies the type of request and the meaning of any remaining arguments. The second argument is almost always the process ID of the target.

Before we can mess with the target process, we need to attach to it:

void breakfast_attach(pid_t pid) {
  int status;
  ptrace(PTRACE_ATTACH, pid);
  waitpid(pid, &status, 0);
  ptrace(PTRACE_SETOPTIONS, pid, 0, PTRACE_O_TRACEEXIT);
}

The PTRACE_ATTACH request will stop the target with SIGSTOP. We wait for the target to receive this signal. We also request to receive a notification when the target exits, which I'll explain below.

Note that ptrace and waitpid could fail in various ways. In a real application you'd want to check the return value and/or errno. For brevity I've omitted those checks in this article.

We use another ptrace request to get the value of the target's instruction pointer:

target_addr_t breakfast_getip(pid_t pid) {
  long v = ptrace(PTRACE_PEEKUSER, pid, sizeof(long)*REGISTER_IP);
  return (target_addr_t) (v - TRAP_LEN);
}

Since the target process is suspended, it's not running on any CPU, and its instruction pointer is not stored in a real CPU register. Instead, it's saved to a "user area" in kernel memory. We use the PTRACE_PEEKUSER request to read a machine word from this area at a specified byte offset. The constants in sys/regs.h give the order in which registers appear, so we simply multiply by sizeof(long).

After we hit a breakpoint, the saved IP points to the instruction after the trap instruction. When we resume execution, we'll go back and execute the original instruction that we overwrote with the trap. So we subtract TRAP_LEN to give the true address of the next instruction.

Creating breakpoints

We need to remember two things about a breakpoint: the address of the code we replaced, and the code which originally lived there.

struct breakpoint {
  target_addr_t addr;
  long orig_code;
};

To enable a breakpoint, we save the original code and insert our trap instruction:

static void enable(pid_t pid, struct breakpoint *bp) {
  long orig = ptrace(PTRACE_PEEKTEXT, pid, bp->addr);
  ptrace(PTRACE_POKETEXT, pid, bp->addr, (orig & TRAP_MASK) | TRAP_INST);
  bp->orig_code = orig;
}

The PTRACE_PEEKTEXT request reads a machine word from the target's code address space, which is named "text" for historical reasons. PTRACE_POKETEXT writes to that space. On x86 Linux there's actually no difference between code and data spaces, so PTRACE_PEEKDATA and PTRACE_POKEDATA would work just as well.

Creating a breakpoint is straightforward:

struct breakpoint *breakfast_break(pid_t pid, target_addr_t addr) {
  struct breakpoint *bp = malloc(sizeof(*bp));
  bp->addr = addr;
  enable(pid, bp);
  return bp;
}

To disable a breakpoint we just write back the saved word:

static void disable(pid_t pid, struct breakpoint *bp) {
  ptrace(PTRACE_POKETEXT, pid, bp->addr, bp->orig_code);
}

Running and stepping

Once we attach to the target, its execution stops. Here's how to resume it:

static int run(pid_t pid, int cmd);  /* defined momentarily */

int breakfast_run(pid_t pid, struct breakpoint *bp) {
  if (bp) {
    ptrace(PTRACE_POKEUSER, pid, sizeof(long)*REGISTER_IP, bp->addr);

    disable(pid, bp);
    if (!run(pid, PTRACE_SINGLESTEP))
      return 0;
    enable(pid, bp);
  }
  return run(pid, PTRACE_CONT);
}

We ask ptrace to continue execution — but if we're resuming from a breakpoint, we have to do some cleanup first. We rewind the instruction pointer so that the next instruction to execute is at the breakpoint. Then we disable the breakpoint and make the target execute just one instruction. Once we're past the breakpoint, we can re-enable it for next time. run returns 0 if the target exits, which theoretically could happen during our single-step.

The last piece of the puzzle is run:

static int run(pid_t pid, int cmd) {
  int status, last_sig = 0, event;
  while (1) {
    ptrace(cmd, pid, 0, last_sig);
    waitpid(pid, &status, 0);

    if (WIFEXITED(status))
      return 0;

    if (WIFSTOPPED(status)) {
      last_sig = WSTOPSIG(status);
      if (last_sig == SIGTRAP) {
        event = (status >> 16) & 0xffff;
        return (event == PTRACE_EVENT_EXIT) ? 0 : 1;
      }
    }
  }
}

cmd is either PTRACE_CONT or PTRACE_SINGLESTEP. In the latter case, the OS will set a control bit to make the CPU raise interrupt 3 after one instruction has completed.

ptrace will let the target run until it exits or receives a signal. We use the wait macros to see what happened. If the target exited, we simply return 0. If it received a signal, we check which one. Anything other than SIGTRAP should be delivered through to the target, which we accomplish by passing the signal number to the next ptrace call.

In the case of SIGTRAP, we check bits 16-31 of status for the value PTRACE_EVENT_EXIT, which indicates that the target is about to exit. Recall that we requested this notification by setting the option PTRACE_O_TRACEEXIT. You'd think (at least, I did) that checking WIFEXITED should be enough. But I ran into a problem where delivering a fatal signal to the target would make the tracing process loop forever. I fixed it by adding this extra check. I'm certainly no ptrace expert and I'd be interested to hear more about what's going on here.

Testing it

That's all for breakfast.c! Now let's write a small test.c.

#include <sys/types.h>
#include <signal.h>
#include <unistd.h>
#include <stdio.h>
#include "breakfast.h"

void child();
void parent(pid_t);

int main() {
  pid_t pid;
  if ((pid = fork()) == 0)
    child();
  else
    parent(pid);
  return 0;
}

We fork two separate processes, a parent and a child. The child will do some computation we hope to observe. Let's make it compute the famous factorial function:

int fact(int n) {
  if (n <= 1)
    return 1;
  return n * fact(n-1);
}

void child() {
  kill(getpid(), SIGSTOP);
  printf("fact(5) = %d\n", fact(5));
}

The parent is going to invoke PTRACE_ATTACH and send the child SIGSTOP, but the child might finish executing before the parent gets a chance to call ptrace. So we make the child preemptively stop itself. This isn't a problem when attaching to a long-running process.

Actually, for this fork-and-trace-child pattern, we should be using PTRACE_TRACEME. But we implemented a minimal library, so we work with what we have.

The parent uses breakfast to place breakpoints in the child:

void parent(pid_t pid) {
  struct breakpoint *fact_break, *last_break = NULL;
  void *fact_ip = fact, *last_ip;
  breakfast_attach(pid);
  fact_break = breakfast_break(pid, fact_ip);
  while (breakfast_run(pid, last_break)) {
    last_ip = breakfast_getip(pid);
    if (last_ip == fact_ip) {
      printf("Break at fact()\n");
      last_break = fact_break;
    } else {
      printf("Unknown trap at %p\n", last_ip);
      last_break = NULL;
    }
  }
}

In principle we can use breakfast to trace any running process we own, even if we don't have its source code. But we still need a way to find interesting breakpoint addresses. Here it's as easy as we could want: the child forked from our process, so we get the address of its fact function just by taking the address of our fact function.

Let's try it out:

$ gcc -Wall -o test test.c breakfast.c$ ./testBreak at fact()Break at fact()Break at fact()Break at fact()Break at fact()fact(5) = 120

We can see the five calls to fact.

Inspecting the target

The ability to count function calls is already useful for performance profiling. But we'll usually want to get more information out of the stopped target. Let's see if we can read the arguments passed to fact. This part will be specific to 64-bit x86, though the idea generalizes.

Each architecture defines a C calling convention, which specifies how function arguments are passed, often using a combination of registers and stack slots. On 64-bit x86, the first argument is passed in the RDI register. You can verify this by running objdump -d test and looking at the disassembled code for fact.

So we'll modify test.c to read this register:

#include <sys/ptrace.h>
#include <sys/reg.h>
...
    if (last_ip == fact_ip) {
      int arg = ptrace(PTRACE_PEEKUSER, pid, sizeof(long)*RDI);
      printf("Break at fact(%d)\n", arg);
...

And we get:

$ gcc -Wall -o test test.c breakfast.c$ ./testBreak at fact(5)Break at fact(4)Break at fact(3)Break at fact(2)Break at fact(1)fact(5) = 120

Cool!

Notes

This is my first non-trivial project using ptrace, and it's likely I got some details wrong. Please reply with corrections or advice and I'll update the code if necessary.

This code is available for download at GitHub. It's released under a BSD license, which imposes few restrictions. It's definitely not suitable as-is for production use, but feel free to take the code and build something more with it.

The x86 architecture has dedicated registers for setting breakpoints, subject to various limitations. We ignored this feature in favor of the more flexible software breakpoints. Hardware breakpoints can do some things this technique can't, such as breaking on reads from a particular memory address.

Type-level Fix and generic folds

2010-12-04T19:14:00.000-08:00

This article describes a fixpoint combinator for Haskell types, with some justification of why such a thing would be useful. It's meant to be an accessible introduction to material which is already covered elsewhere.

You'll need some basic Haskell knowledge, such as how to define and use polymorphic data types. And you might want to first study the value-level fix function. An excellent introduction can be found here.

This code is standard Haskell 98, except in the "GHC tricks" section.

Recursive types

Many useful data types are recursive:

data Natural = Zero | Succ Natural     -- natural numbers
data List a  = Nil  | Cons a (List a)  -- lists
data Tree a  = Node a [Tree a]         -- rose trees

In each case, the type declared on the left-hand side of "=" is mentioned on the right-hand side. This is fine, but just for fun let's see if there's another way.

We can imagine a language with a keyword "self", which refers to the type being defined:

data Natural = Zero | Succ self
data List a  = Nil  | Cons a self
data Tree a  = Node a [self]

Haskell doesn't have this keyword, so we'll make self into a type parameter:

data NaturalF self = Zero | Succ self
data ListF a  self = Nil  | Cons a self
data TreeF a  self = Node a [self]

We've changed the names because these aren't the types we wanted. They're non-recursive and each has an extra parameter.

To build the recursive type Natural using the non-recursive NaturalF, we need to "tie the knot":

type Natural = NaturalF (NaturalF (NaturalF (NaturalF ...)))

i.e.

type Natural = NaturalF Natural

But this is disallowed:

foo.hs:3:0:    Cycle in type synonym declarations:      foo.hs:3:0-30: type Natural = NaturalF Natural

Since type defines a transparent synonym, this would generate an infinitely-large type-expression. Haskell disallows those. We need to hide the recursion behind a data constructor:

newtype Natural = MkNat (NaturalF Natural)

When we use this type, the "wrapper" constructor MkNat alternates with the constructors of NaturalF:

toNatural :: Integer -> Natural
toNatural 0 = MkNat Zero
toNatural n = MkNat . Succ $ toNatural (n-1)

Type-level `Fix`

So far this is just pointless complexity. We've replaced one type:

data Natural = Zero | Succ Natural

with two types:

data    NaturalF self = Zero | Succ self
newtype Natural       = MkNat (NaturalF Natural)

Notice that the definition of Natural doesn't depend on the structure of natural numbers. Let's rename Natural to Fix and MkNat to In:

newtype Fix   = In (NaturalF Fix)

Then we'll abstract out NaturalF as a type parameter f:

newtype Fix f = In (f (Fix f))

Now Natural is just a synonym:

type Natural = Fix NaturalF

When we create values of type Natural, we still have "wrapper" constructors:

toNatural :: Integer -> Natural
toNatural 0 = In Zero
toNatural n = In . Succ $ toNatural (n-1)

Properties of `Fix`

Recall that the value-level fix function can be defined as

fix f = f (fix f)

This looks a lot like our type

newtype Fix f = In (f (Fix f))

The similarity is reflected at the type and kind levels:

GHCi> :t fixfix :: (a -> a) -> aGHCi> :k FixFix :: (* -> *) -> *

fix is a function which takes another function. Fix is a type constructor which takes another type constructor.

The types Fix f and f (Fix f) are isomorphic, and we'll want to convert between them. In one direction we'll use the data constructor

In :: f (Fix f) -> Fix f

and the other direction is easily defined as

out :: Fix f -> f (Fix f)
out (In x) = x

Folding lists

So why bother writing types this way? Well, we usually don't. But factoring the recursion out of our data type enables us to factor the recursion out of our functions.

Consider taking the sum of a list:

sumF :: (Num a) => List a -> a
sumF (In Nil)         = 0
sumF (In (Cons x xs)) = x + sumF xs

We traverse the structure with the combining function (+), making a recursive call at each recursive position of the data type. How can we capture this pattern?

If we're reducing a value of type List T to get a result of type R, then a value of type ListF T R represents a list node with the recursive call's result already in place. So we might look for a function like

cata :: (ListF a b -> b) -> List a -> b

If we had this, we could write

sumF :: (Num a) => List a -> a
sumF = cata f where
  f Nil = 0
  f (Cons x xs) = x + xs

The function f specifies one step of reduction, and the magical function cata handles the rest.

The recursion pattern captured by cata is very similar to foldr on the built-in list type. In fact, we can implement foldr with cata:

foldrF :: (a -> b -> b) -> b -> List a -> b
foldrF f z = cata g where
  g Nil = z
  g (Cons x xs) = f x xs

If foldr is good enough, then what's so special about cata? The answer is that cata actually has this very general type:

cata :: (Functor f) => (f a -> a) -> Fix f -> a

I'll give the one-line definition of cata later. Writing it yourself is a fun puzzle.

We can see that we'll need this instance:

instance Functor (ListF a) where
  fmap _ Nil         = Nil
  fmap f (Cons x xs) = Cons x (f xs)

This fmap just takes a function and applies it to the self-typed field, if any. It's not recursive and has no relation to the [] instance of Functor, for which fmap = map.

Folding natural numbers

Let's put our generic cata to good use on another type. Every natural number is either zero, or the successor of another natural number:

data NaturalF self = Zero | Succ self
type Natural = Fix NaturalF

Again, fmap applies a function to each self-typed field:

instance Functor NaturalF where
  fmap _ Zero     = Zero
  fmap f (Succ v) = Succ (f v)

We can then use cata to define conversion to Integer:

fromNatural :: Natural -> Integer
fromNatural = cata f where
  f Zero     = 0
  f (Succ n) = 1 + n

as well as addition and multiplication:

add :: Natural -> Natural -> Natural
add n = cata f where
  f Zero     = n
  f (Succ m) = In $ Succ m

mul :: Natural -> Natural -> Natural
mul n = cata f where
  f Zero     = In Zero
  f (Succ m) = add n m

Testing it:

GHCi> fromNatural (add (toNatural 2) (mul (toNatural 3) (toNatural 5)))17

Folding rose trees

Let's try one more structure, namely rose trees:

data TreeF a self = Node a [self]
type Tree  a      = Fix (TreeF a)

instance Functor (TreeF a) where
  fmap f (Node v xs) = Node v (map f xs)

As before, fmap applies a function to each self-typed field. Since we have a list of these, we use map.

We define some traversals:

valuesF :: Tree a -> [a]
valuesF = cata (\(Node x xs) -> x : concat xs)

depthF :: Tree a -> Integer
depthF = cata (\(Node _ xs) -> 1 + maximum (0:xs))

and test them:

GHCi> let n x = In . Node xGHCi> let t1 = n 'a' []GHCi> valuesF t1"a"GHCi> depthF t11GHCi> let t2 = n 'a' [n 'b' [], n 'c' [n 'd' []], n 'e' []]GHCi> valuesF t2"abcde"GHCi> depthF t23

The definition of `cata`

As promised:

cata :: (Functor f) => (f a -> a) -> Fix f -> a
cata g = g . fmap (cata g) . out

Generic programming

We usually manipulate data because we care about its meaning. But some operations only depend on the data's structure. If we have some complex nested data value, we might want to pretty-print it, or serialize it with binary, or collect all the Strings buried within. These are concepts which apply to any data type you like. Shouldn't our code be similarly adaptable?

This idea is called generic programming. Implementations in Haskell include syb, uniplate, multiplate, instant-generics, GHC's built-in support, and many other efforts.

Each of these libraries demands certain information about each data type, and in return provides some generic operations. In this view, our above development of Fix constitutes a really tiny generics library. We require

that you write recursive types as Fix F, where F is non-recursive, and
that you provide instance Functor F.

In return, you get the single generic operation cata, which reduces a recursive structure to a "summary" value, by replacing each constructor with an arbitrary computation.

Our toy generics library is simple, but it doesn't actually save much work. Using Fix adds noise to the code, and our Functor instances aren't much shorter than a custom fold function for each type. Conversely, our uses of cata aren't much simpler than explicit recursive traversals.

Practicality aside, I think it's cool to define folding "once and for all". I'm interested to learn (from you!) what else can be done with Fix'd data or similar techniques.

For a non-trivial generics library built around similar ideas, see multirec and the accompanying paper.

Odds and ends

Nomenclature

Fix is sometimes named Mu, because the Greek letter μ is used to represent least fixed points.

cata stands for catamorphism.

Some of this stuff shows up in "Evolution of a Haskell Programmer". They build a factorial function using not only catamorphisms but also zygomorphisms and paramorphisms. Tell your friends!

In TaPL there's a discussion of equirecursive versus isorecursive type systems. Haskell is in the latter category. Fix f is isomorphic to, but not equivalent to, f (Fix f). The functions In and out provide the isomorphism. If we defined Fix using data instead of newtype, the isomorphism would fail because of the extra value In ⊥.

Isomorphisms

NaturalF is isomorphic to Maybe:

data NaturalF self = Zero    | Succ self
data Maybe    a    = Nothing | Just a

The Functor instances are also the same.

NaturalF is also isomorphic to ListF (), ignoring undefined values.

GHC tricks

While the above is standard Haskell, we can pull a few more tricks by using a couple of GHC extensions.

GHC could help us reduce boilerplate by deriving Functor instances.

With some GHC extensions, we can write Show for Fix:

instance (Show (f (Fix f))) => Show (Fix f) where
  showsPrec p (In x) = showParen (p >= 11) (showString "In " . showsPrec 11 x)

Testing it:

GHCi> toNatural 3In (Succ (In (Succ (In (Succ (In Zero))))))

This requires -XFlexibleContexts -XUndecidableInstances. The latter is required because f (Fix f) is a syntactically larger type than Fix f, so GHC can't convince itself that instance resolution will definitely terminate. For well-behaved arguments to Fix, it works out. We'll have something like

-- from above
instance (Show (f (Fix f)))
       => Show (Fix f)

-- from 'deriving (Show)' on NaturalF
instance (Show self)
       => Show (NaturalF self)

which instantiates to

instance (Show (NaturalF (Fix NaturalF)))
       => Show (Fix NaturalF)

instance (Show (Fix NaturalF))
       => Show (NaturalF (Fix NaturalF))

and GHC is clever enough to resolve the mutual recursion.

Remarkably, GHC can derive this instance with -XStandaloneDeriving:

deriving instance (Show (f (Fix f))) => Show (Fix f)

That's how I got the precedence-correct code above, via -ddump-deriv.

Obtaining the name of a function in Haskell

2010-11-04T18:52:00.000-07:00

Someone on #haskell asked whether you could recover a function's name from its value, i.e.

GHCi> name id"id"GHCi> name map"map"

This is easy in some languages. But Haskell is not designed to provide this kind of run-time information. We'll need some non-portable hacks. I tested this code with GHC 6.12.1 on amd64 Linux; see below for portability notes.

How it works

Most values in GHC Haskell are represented by a closure. Closures representing functions are common in functional-language compilers, but GHC extends the idea to cover many other sorts of values.

Closures exist to store data. But the run-time system also needs operational information about each closure: how to garbage-collect it, how to force its evaluation, etc. This information is known at compile time and is shared between many closures. All algebraic values with the same constructor will share this information, as will all function values created from the same lambda in the program's source.

So each closure stores a pointer to an info table, holding this operational information. Info tables are generated at compile time, and stored as part of an executable's read-only data section. This means that they have statically-known addresses, with associated names in the executable's symbol table. We'll use these symbol names to name our functions.

We can dump an executable's symbol table with nm:

$ nm -f posix foo...ghczmprim_GHCziBool_Bool_closure_tbl D 0000000000749978 ghczmprim_GHCziBool_False_closure D 0000000000749970 ghczmprim_GHCziBool_False_static_info T 00000000004e4dd8 ghczmprim_GHCziBool_True_closure D 0000000000749990 ghczmprim_GHCziBool_True_static_info T 00000000004e4d80 ghczmprim_GHCziDebug_debugErrLn1_closure D 00000000007499a0 ghczmprim_GHCziDebug_debugErrLn1_info T 00000000004e4ec0 ...

Haskell identifiers can contain characters not allowed in symbol names. GHC uses a name-mangling scheme to build symbol names. For example, the first symbol above decodes to

ghc-prim_GHC.Bool_Bool_closure_tbl

Reading the symbols

We'll use vacuum to inspect heap objects. vacuum can do a lot of cool things, like visualize heap structures from a running Haskell program. We'll only use one of vacuum's functions here.

We'll also use the GHC API to un-mangle symbol names. GHC is a 20-year effort that has evolved alongside the Haskell language. It follows some legacy conventions like a mostly-flat module hierarchy. So the module we need is named simply Encoding.

import Data.Word
import Data.Maybe
import Text.Printf
import Control.Parallel ( pseq )

import qualified Data.Map           as Map
import qualified System.Posix.Files as Posix
import qualified System.Process     as Proc
import qualified Foreign.Ptr        as Ptr

import qualified GHC.Vacuum         as Vac
import qualified Encoding           as GHC

We invoke nm as a subprocess and parse its output:

type Symbols = Map.Map Word String

getSymbols :: IO Symbols
getSymbols = do
  exe <- Posix.readSymbolicLink "/proc/self/exe"
  out <- Proc.readProcess "nm" ["-f", "posix", exe] ""
  let offset = 0x10
  let f (sym:_:addr:_) = Just (read ("0x"++addr) - offset, GHC.zDecodeString sym)
      f _ = Nothing
  return . Map.fromList . catMaybes . map (f . words) . lines $ out

We're using the Linux proc filesystem to get a symbolic link to our application's executable.

The symbols in memory appear at an address 0x10 = 16 bytes or 2 machine words lower than in the executable's symbol table. I'm not sure why; perhaps it's because of GHC's "tables next to code" optimization.

Resolving a symbol

Once we have the symbol table, looking up a value is relatively easy:

name :: Symbols -> a -> String
name syms x = fromMaybe unk $ Map.lookup ptr syms where
  ptr = x `pseq` (fromIntegral . Ptr.ptrToWordPtr . Vac.getInfoPtr $ x)
  unk = printf "<unknown info table at 0x%016x>" ptr

We use vacuum to get the value's info table pointer, convert this to a Word, then look it up in the symbol table.

We explicitly evaluate x with pseq, to avoid seeing a thunk.

Testing it

We'll test with

$ ghc --make name.hs -package ghc$ ./name

Each test below is commented with the expected output. First, let's try a few non-function values:

main :: IO ()
main = do
  syms <- getSymbols
  let test = putStrLn . name syms
  test 3            -- integer-gmp_GHC.Integer.Type_S#_con_info
  test (3 :: Int)   -- ghc-prim_GHC.Types_I#_static_info
  test "xyz"        -- ghc-prim_GHC.Types_:_con_info

GHC defaults to 3 :: Integer, as -Wall will tell you. As we see, Integer and Int are both implemented as algebraic data:

data Integer
  = S# Int#
  | J# Int# ByteArray#

data Int = I# Int#

The string "xyz" is a list built out of (:) cells.

Next let's try a few functions:

  test map          -- base_GHC.Base_map_info
  test getChar      -- base_System.IO_getChar_info
  test (+)          -- integer-gmp_GHC.Integer_plusInteger_info

Not bad! (+) defaults to operating on Integer, and GHC inlines the type class dictionary, giving us the underlying plusInteger function.

Now let's see the limits of this technique:

  test (\_ -> 'x')  -- s1jD_info
  test (const 'x')  -- stg_PAP_info
  test test         -- stg_PAP_info

Our lambda expression gets a useless compiler-generated name. The application of const is worse; it uses an info table common to all partial applications. However, we could use vacuum to follow the fields of the PAP closure, which I'll leave as an exercise to the reader. ;)

test itself is also a partial application. It's defined by applying two arguments to the function (.) defined as

(.) f g x = f (g x)

If we eta-expand test:

  let test x = putStrLn $ name syms x

then we'll get another compiler-generated name like s1mh_info.

Notes

This is a hack, and probably not suitable for any serious purpose. Shelling out to nm to get a symbol table is particularly ugly. I tried to use bindings to BFD, but ran into some segfaults.

The above code will work only on 64-bit machines, but could be adapted for 32-bit. I bet the magic offset would change. It works on GHC 6.12.1, and should work on other recent versions, if you can get vacuum to build.

It definitely requires a Unix system, and specifically Linux or something emulating Linux's proc filesystem. You'll need nm from GNU Binutils, which is standard on a system configured for C development.

It won't work if you run strip on your binaries...

Tour of a real toy Haskell program, part 2

2010-10-26T01:03:00.000-07:00

This is part 2 of my commentary on detrospector. You can find part 1 here.

Performance claims

I make a lot of claims about performance below, not all of which are supported by evidence. I did perform a fair amount of profiling, but it's impossible to test every combination of data structures. How can I structure my code to make more of my performance hypotheses easily testable?

Representing the source document

When we process a source document, we'll need a data structure to represent large amounts of text.

Haskell's standard String type is simple, maximally lazy, and Unicode-correct. However, as a singly-linked list of individually boxed heap-allocated characters, it's woefully inefficient in space and time.
ByteString is a common alternative, but solves the wrong problem. We need a sequence of Unicode codepoints, not a sequence of bytes.
ByteString.Char8 is a hack and not remotely Unicode-correct.
The text library stores Unicode characters in a ByteString-like packed format. That's exactly what we need.

For more information on picking a string type, see ezyang's writeup.

We choose the lazy flavor of text so that we can stream the input file in chunks without storing it all in memory. Lazy IO is semantically dubious, but I've decided that this project is enough of a toy that I don't care. Note that robustness to IO errors is not in the requirements list. ;)

The enumerator package provides iteratees as a composable, well-founded alternative to the lazy IO hack. The enumerator API was complex enough to put me off using it for the first version of detrospector. After reading Michael Snoyman's enumerators tutorial, I have a slight idea how the library works, and I might use it in a future version.

Representing substrings

We also need to represent the k-character substrings, both when we analyze the source and when we generate text. The requirements here are different.

We expect that k will be small, certainly under 10 — otherwise, you'll just generate the source document! With many small strings, the penalty of boxing each individual Char is less severe.

And we need to support queue-like operations; we build a new substring by pushing a character onto the end, and dropping one from the beginning. For both String and Text, appending onto the end requires a full copy. So we'll use Data.Sequence, which provides sequences with bidirectional append, implemented as finger trees:

-- module Detrospector.Types
type Queue a = S.Seq a

...Actually, I just ran a quick test with Text as the queue type, and it seems not to matter much. Since k is small, the O(k) copy is insignificant. Profiling makes fools of us all, and asymptotic analysis is mostly misused. Anyway, I'm keeping the code the way it is, pending any insight from my clever readers. Perhaps one of the queues from Purely Functional Data Structures would be most appropriate.

Representing character counts

We need a data structure for tabulating character counts. In the old days of C and Eurocentrism, we could use a flat, mutable int[256] indexed by ASCII values. But the range of Unicode characters is far too large for a flat array, and we need efficient functional updates without a full copy.

We could import Data.Map and use Map Char Int. This will build a balanced binary search tree using pairwise Ord comparisons.

But we can do better. We can use the bits of an integer key as a path in a tree, following a left or right child for a 0 or 1 bit, respectively. This sort of search tree (a trie) will typically outperform repeated pairwise comparisons. Data.IntMap implements this idea, with an API very close to Map Int. Our keys are Chars, but we can easily convert using fromIntegral.

-- module Detrospector.Types
import qualified Data.IntMap as IM
...
type FreqTable = IM.IntMap Int

Representing distributions

So we have some frequency table like

IM.fromList [('e', 267), ('t', 253), ('a', 219), ...]

How can we efficiently pick a character from this distribution? We're mapping characters to individual counts, but we really want a map from cumulative counts to characters:

-- module Detrospector.Types
data PickTable = PickTable Int (IM.IntMap Char)

To sample a character from PickTable t im, we first pick a random k such that 0 ≤ k < t, using a uniform distribution. We then find the first key in im which is greater than k, and take its associated Char value. In code:

-- module Detrospector.Types
import qualified System.Random.MWC as RNG
...
sample :: PickTable -> RNG.GenIO -> IO Char
sample (PickTable t im) g = do
  k <- (`mod` t) <$> RNG.uniform g
  case IM.split k im of
    (_, IM.toList -> ((_,x):_)) -> return x
    _ -> error "impossible"

The largest cumulative sum is the total count t, so the largest key in im is t. We know k < t, so IM.split k im will never return an empty map on the right side.

Note the view pattern for pattern-matching an IntMap as if it were an ascending-order list.

The standard System.Random library in GHC Haskell is quite slow, a problem shared by most language implementations. We use the much faster mwc-random package. The only operation we need is picking a uniform Int as an IO action.

We still need a way to build our PickTable:

-- module Detrospector.Types
import Data.List ( mapAccumR )
...
cumulate :: FreqTable -> PickTable
cumulate t = PickTable r $ IM.fromList ps where
  (r,ps) = mapAccumR f 0 $ IM.assocs t
  f ra (x,n) = let rb = ra+n in (rb, (rb, toEnum x))

This code is short, but kind of tricky. For reference:

mapAccumR :: (acc -> x -> (acc, y)) -> acc -> [x] -> (acc, [y])
f :: Int -> (Int, Int) -> (Int, (Int, Char))

f takes an assoc pair from the FreqTable, adds its count to the running sum, and produces an assoc pair for the PickTable. We start the traversal with a sum of 0, and get the final sum r along with our assoc pairs ps.

Representing the Markov chain

So we can represent probability distributions for characters. Now we need a map from k-character substrings to distributions.

Data.Map is again an option, but pairwise, character-wise comparison of our Queue Char values will be slow. What we really want is another trie, with character-based fanout at each node. Hackage has bytestring-trie, which unfortunately works on bytes, not characters. And maybe I should have used TrieMap or list-tries. Instead I used the hashmap package:

-- module Detrospector.Types
import qualified Data.HashMap as H
...
data Chain = Chain Int (H.HashMap (Queue Char) PickTable)

A value Chain k hm maps from subsequences of up to k Chars to PickTables. A lookup of some Queue Char key will require one traversal to calculate an Int hash value, then uses an IntMap to find a (hopefully small) Map of keys with that same hash value.

There is a wrinkle: we need to specify how to hash a Queue, which is just a synonym for S.Seq. This is an orphan instance, which we could avoid by newtype-wrapping S.Seq.

-- module Detrospector.Types
import qualified Data.Hashable as H
import qualified Data.Foldable as F
...
instance (H.Hashable a) => H.Hashable (S.Seq a) where
  {-# SPECIALIZE instance H.Hashable (S.Seq Char) #-}
  hash = F.foldl' (\acc h -> acc `H.combine` H.hash h) 0

This code is lifted almost verbatim from the [a] instance in Data.Hashable.

Serialization

After training, we need to write the Chain to disk, for use in subsequent generation. I started out with derived Show and Read, which was simple but incredibly slow. We'll use binary with ByteString.Lazy — the dreaded lazy IO again!

We start by specifying how to serialize a few types. Here the tuple instances for Binary come in handy:

-- module Detrospector.Types
import qualified Data.Binary as Bin
...
-- another orphan instance
instance (Bin.Binary k, Bin.Binary v, H.Hashable k, Ord k)
       => Bin.Binary (H.HashMap k v) where
  put = Bin.put . H.assocs
  get = H.fromList <$> Bin.get

instance Bin.Binary PickTable where
  put (PickTable n t) = Bin.put (n,t)
  get = uncurry PickTable <$> Bin.get

instance Bin.Binary Chain where
  put (Chain n h) = Bin.put (n,h)
  get = uncurry Chain <$> Bin.get

The actual IO is easy. We use gzip compression, which fits right into the IO pipeline:

-- module Detrospector.Types
import qualified Data.ByteString.Lazy   as BSL
import qualified Codec.Compression.GZip as Z
...
withChain :: FilePath -> (Chain -> RNG -> IO a) -> IO a
withChain p f = do
  ch <- (Bin.decode . Z.decompress) <$> BSL.readFile p
  RNG.withSystemRandom $ f ch

writeChain :: FilePath -> Chain -> IO ()
writeChain out = BSL.writeFile out . Z.compress . Bin.encode

Training the chain

I won't present the whole implementation of the train subcommand, but here's a simplification:

-- module Detrospector.Modes.Train
import qualified Data.Text     as Txt
import qualified Data.Text.IO  as Txt
import qualified Data.HashMap  as H
import qualified Data.IntMap   as IM
import qualified Data.Sequence as S
import qualified Data.Foldable as F
...
train Train{num,out} = do
  (_,h) <- Txt.foldl' roll (emptyQ, H.empty) ys <$> Txt.getContents
  writeChain out . Chain num $ H.map cumulate h where

  roll (!s,!h) x
    = (shift num x s, F.foldr (H.alter $ ins x) h $ S.tails s)
 
  ins x Nothing  = Just $! sing x
  ins x (Just v) = Just $! incr x v
 
  sing x = IM.singleton (fromEnum x) 1
 
  incr x = IM.alter f $ fromEnum x where
    f Nothing  = Just 1
    f (Just v) = Just $! (v+1)

Before generating PickTables, we build a HashMap of FreqTables. We fold over the input text, accumulating a pair of (last characters seen, map so far). Since foldl' is only strict to weak head-normal form (WHNF), we use bang patterns on the fold function roll to force further evaluation. RWH discusses the same issue.

shift (from Detrospector.Types) pushes a new character into the queue, and drops the oldest character if the size exceeds num. We add one count for the new character x, both for the whole history s and each of its suffixes.

We're using alter where perhaps a combination of lookup and insert would be more natural. This is a workaround for a subtle laziness-related space leak, which I found after much profiling and random mucking about. When you insert into a map like so:

let mm = insert k (v+1) m

there is nothing to force v+1 to WHNF, even if you force mm to WHNF. The leaves of our tree end up holding large thunks of the form ((((1+1)+1)+1)+...).

The workaround is to call alter with Just $! (v+1). We know that the implementation of alter will pattern-match on the Maybe constructor, which then triggers WHNF evaluation of v+1 because of ($!). This was tricky to figure out. Is there a better solution, or some different way I should approach this problem? It seems to me that Data.Map and friends are generally lacking in strictness building blocks.

The end!

Thanks for reading! Here's an example of how not to write the same program:

module Main where{import Random;(0:y)%(R p _)=y%p;(1:f)%(R _ n)=f%n;[]%(J x)=x;b
[p,v,k]=(l k%).(l v%).(l p%);main=getContents>>=("eek"#).flip(.:"eek")(y.y.y.y$0
);h k=toEnum k;;(r:j).:[k,v,b]=(j.:[v,b,r]).((k?).(v?).(b?)$(r?)m);[].:_=id;(!)=
iterate;data R a=J a|R(R a)(R a);(&)i=fmap i;k b y v j=let{t=l b%v+y;d f|t>j=b;d
f=k(m b)t v j}in d q;y l=(!!8)(q R!J l);q(+)b=b+b;p(0:v)z(R f x)=R(p v z f)x;p[]
z(J x)=J(z x);p(1:v)z(R n k)=R n$p v z k;m = succ;l p=tail.(snd&).take 9$((<<2).
fst)!(fromEnum p,0);(?)=p.l;d@[q,p,r]#v=let{u=b d v;j=i u;(s,o)|j<1=((97,122),id
)|h 1=((0,j-1),(k 0 0 u))}in do{q<-(h.o)&randomRIO s;putChar q;[p,r,q]#v};i(J q)
=q;i(R n k)=i n+i k;(<<)=divMod} -- finite text on stdin  © keegan oct 2010 BSD3

Tour of a real toy Haskell program, part 1

2010-10-26T01:00:00.000-07:00

Haskell suffers from a problem I'll call the Fibonacci Gap. Many beginners start out with a bunch of small mathematical exercises, develop that skillset, and then are at a loss for what to study next. Often they'll ask in #haskell for an example of a "real Haskell program" to study. Typical responses include the examples in Real World Haskell, or the Design and Implementation of XMonad talk.

This is my attempt to provide another data point: a commentary on detrospector, my text-generating toy. It's perhaps between the RWH examples and xmonad in terms of size and complexity. It's not a large program, and certainly not very useful, but it does involve a variety of real-world concerns such as Unicode processing, choice of data structures, strictness for performance, command-line arguments, serialization, etc. It also illustrates my particular coding style, which I don't claim is the One True Way, but which has served me well in a variety of Haskell projects over the past five years.

Of course, I imagine that experts may disagree with some of my decisions, and I welcome any and all feedback.

I haven't put hyperlinked source online yet, but you can grab the tarball and follow along.

This is part 1, covering style and high-level design. Part 2 addresses more details of algorithms, data structures, performance, etc.

The algorithm

detrospector generates random text conforming to the general style and diction of a given source document, using a quite common algorithm.

First, pick a "window size" k. We consider all k-character contiguous substrings of the source document. For each substring w, we want to know the probability distribution of the next character. In other words, we want to compute values for

P(next char is x | last k chars were w).

To compute these, we scan through the source document, remembering the last k characters as w. We build a table of next-character occurrence counts for each substring w.

These tables form a Markov chain, and we can generate random text by a random walk in this chain. If the last k characters we generated were w, we choose the next character randomly according to the observed distribution for w.

So we can train it on Accelerando and get output like:

addressed to tell using back holes everyone third of the glances waves and diverging him into the habitat. Far beyond Neptune, I be?" asks, grimy. Who's whether is headquarters I need a Frenchwoman's boyfriend go place surrected in the whole political-looking on the room, leaving it, beam, the wonderstood. The Chromosome of upload, does enough. If this one of the Vile catches agree."

Design requirements

We can only understand a program's design in light of the problem it was meant to solve. Here's my informal list of requirements:

detrospector should generate random text according to the above algorithm.
We should be able to invoke the text generator many times without re-analyzing the source text each time.
detrospector should handle all Unicode characters, using the character encoding specified by the system's locale.
detrospector should be fast without sacrificing clarity, whatever that means.

Wearing my customer hat, these are axioms without justification. Wearing my implementor hat, I will have to justify design decisions in these terms.

Style

In general I import modules qualified, using as to provide a short local name. I make an exception for other modules in the same project, and for "standard modules". I don't have a precise definition of "standard module" but it includes things like Data.Maybe, Control.Applicative, etc.

The longest line in detrospector is 70 characters. There is no hard limit, but more than about 90 is suspect.

Indentation is two spaces. Tabs are absolutely forbidden. I don't let the indentation of a block construct depend on the length of a name, thus:

foo x = do
  y <- bar x
  baz y

and

bar x =
  let y = baz x
      z = quux x y
  in  y z

This avoids absurd left margins, looks more uniform, and is easier to edit.

I usually write delimited syntax with one item per line, with the delimiter prefixed:

{-# LANGUAGE
    ViewPatterns
  , PatternGuards #-}

and

data Mode
  = Train { num   :: Int
          , out   :: FilePath }
  | Run   { chain :: FilePath }

Overriding layout is sometimes useful, e.g.:

look = do x <- H.lookup t h; return (x, S.length t)

(With -XTupleSections we could write

look = (,S.length t) <$> H.lookup t h

but that's just gratuitous.)

I always write type signatures on top-level bindings, but rarely elsewhere.

Module structure

I started out with a single file, which quickly became unmanageable. The current module layout is:

Detrospector.  Types      types and functions used throughout  Modes      a type with one constructor per mode  Modes.    Train    train the Markov chain    Run      generate random text    Neolog   generate neologisms  Main       command-line parsing

There is also a Main module in detrospector.hs, which simply invokes Detrospector.Main.main.

Modules I write tend to fall into two categories: those which export nearly everything, and those which export only one or two things. The former includes "utility" modules with lots of small types and function definitions. The latter includes modules providing a specific piece of functionality. A parsing module might define three dozen parsers internally, but will only export the root of the grammar.

An abstract data type might fall into a third category, since they can export a large API yet have lots of internal helpers. But I don't write such modules very often.

Detrospector.Types is in the first category. Most Haskell projects will have a Types module, although I'm somewhat disappointed that I let this one become a general grab-bag of types and utility functions.

The rest fall into the second category. Each module in Detrospector.Modes.* exports one function to handle that mode. Detrospector.Main exports only main.

Build setup

This was actually my first Cabal project, and the first thing I uploaded to Hackage. I think Cabal is great, and features like package-visibility management are useful even for small local projects.

In my cabal file I set ghc-options: -Wall, which enables many helpful warnings. The project should build with no warnings, but I use the OPTIONS_GHC pragma to disable specific warnings in specific files, where necessary.

I also run HLint on my code periodically, but I don't have it integrated with Cabal.

I was originally passing -O2 to ghc. Cabal complained that it's probably not necessary, which was correct. The Cabal default of -O performs just as well.

I'm using Git for source control, which is neither here nor there.

Command-line parsing

detrospector currently has three modes, as listed above. I wanted to use the "subcommand" model of git, cabal, etc. So we have detrospector train, detrospector run, etc. The cmdargs package handles this argument style with a low level of boilerplate.

The "impure" interface to cmdargs uses some dark magic in the operator &= in order to attach annotations to arbitrary record fields. The various caveats made me uneasy, so I opted for the slightly more verbose "pure" interface, which looks like this:

-- module Detrospector.Main
import qualified System.Console.CmdArgs as Arg
import           System.Console.CmdArgs((+=),Annotate((:=)))
...
modes  = Arg.modes_  [train,run,neolog]
      += Arg.program "detrospector"
      += Arg.summary "detrospector: Markov chain text generator"
      += Arg.help    "Build and run Markov chains for text generation"
  where

  train = Arg.record Train{}
    [ num := 4
          += Arg.help "Number of characters lookback"
    , out := error "Must specify output chain"
          += Arg.typFile
          += Arg.help "Write chain to this file" ]
    += Arg.help "Train a Markov chain from standard input"

  run = Arg.record Run{}
    [ chain := error "Must specify input chain"
            += Arg.argPos 0         
            += Arg.typ "CHAIN_FILE" ]
    += Arg.help "Generate random text"
  ...

This tells cmdargs how to construct values of my record type Detrospector.Modes.Mode. We get help output for free:

$ ./dist/build/detrospector/detrospector -?detrospector: Markov chain text generatordetrospector [COMMAND] ... [OPTIONS]  Build and run Markov chains for text generationCommon flags:  -? --help        Display help message  -V --version     Print version informationdetrospector train [OPTIONS]  Train a Markov chain from standard input  -n --num=INT     Number of characters lookback  -o --out=FILE    Write chain to this filedetrospector run [OPTIONS] CHAIN_FILE  Generate random text...

My use of error here is hacky and leads to a bug that I recently discovered. When the -o argument to train is invalid or missing, the error is not printed until the (potentially time-consuming) analysis is completed. Only then is the record's field forced.

To be continued...

...right this way.

Typing mathematical characters in X

2010-10-23T00:21:00.000-07:00

Unicode provides many useful characters for mathematics. If you've studied the traditional notation, an expression like « Γ ⊢ Λt.λ(x:t).x : ∀t. t → t » is much more readable than an ASCII equivalent. However, most systems don't provide an easy way to enter these characters.

The compose key feature of X Windows provides a nice solution on Linux and other UNIX systems. Compose combinations are easy-to-remember mnemonics, like -> for →, and an enormous number of characters are available with just a few keystrokes.

Setting it up

I cooked up a config file with my most-used mathematical symbols. With recent Xorg, you can drop this file in ~/.XCompose, restart X, and you should be good to go.

The include line will pull in your system-wide configuration, e.g. /usr/share/X11/locale/en_US.UTF-8/Compose. This already contains many useful characters. I was going to add <3 for ♥ and CCCP for ☭, but I found that Debian already provides these.

GTK has its own input handling. To make it defer to X, I had to add an environment variable in ~/.xsession:

export GTK_IM_MODULE="xim"

The "Fn" key on recent ThinkPads makes a good compose key. It normally acts as a modifier key in hardware, but will send a keycode to X when pressed and released by itself. I used this xmodmap setting:

keycode 151=Multi_key

Tweaking the codes

Being boilerplate-averse, I specified the key combinations in a compact format which is processed by this Haskell script.

Obviously, not everyone will like my choice of key combinations. If you tweak the file and come up with something particularly nice, I'd like to see it. If you can't run Haskell code for whatever reason, it's not too hard to edit the generated XCompose file.

Though my use of Haskell here may seem gratuitous, I actually started writing this script in Python, but ran into trouble with Python 2's inconsistent treatment of Unicode text. Using Haskell's String type with GHC ≥ 6.12 will Just Work, at least until you care about performance.

Alternatives

If you don't like this solution, SCIM provides an input mode which uses L^aT_eX codes.

Quantification in Haskell

2010-10-19T18:18:00.000-07:00

I wrote an answer over at Stack Overflow that somehow grew to article length. Here it is, recorded for posterity.

I'll paraphrase the original question as:

What's the difference between the types forall a. [a] and [forall a. a], and
How does this relate to existential types?

Well, the short answer to the second question is "It doesn't relate". But there's still a natural flow in discussing both topics.

Notation

Polymorphic values are essentially functions on types, but Haskell's syntax leaves both type abstraction and type application implicit. To better understand what's going on, we'll use a pidgin syntax of Haskell combined with a typed lambda calculus like System F.

In System F, polymorphism involves explicit type abstractions and type application. For example, the familiar map function would be written as:

map :: ∀a. ∀b. (a → b) → [a] → [b]
map = Λa. Λb. λ(f :: a → b). λ(xs :: [a]). case xs of
  [] → []
  (y:ys) → f y : map @a @b f ys

map is now a function of four arguments: types a and b, a function, and a list.

At the term level, we have one new syntactic form: Λ (upper-case lambda). A polymorphic value is a function from types to values, written Λx. e, much as a function from values to values is written λx. e.

At the type level, we have the symbol ∀ ("for all"), corresponding to GHC's forall keyword. A term containing Λ gives rise to a type containing ∀, just as a term containing λ gives rise to a type containing →.

Since we have functions of types, we also need application of types. I use the notation @a (as in GHC Core) to denote application of a type argument. So we might use map like so:

map @Char @Int ord "xzy"

Note also that I've provided an explicit type on each λ abstraction. Type inference for System F is, in general, undecidable (Wells, 1999). The restricted form of polymorphism available in vanilla Haskell 98 has decidable inference, but you lose this if you enable certain GHC extensions like RankNTypes.

We don't need annotations on Λ abstractions, because we have only one "kind" of type. In actual Haskell, or a calculus like System Fω, we also have type constructors, and we need a system of kinds to describe how they combine. We'll ignore this issue here.

So a value of polymorphic type is like a function from types to values. The caller of a polymorphic function gets to choose a type argument, and the function must comply.

One last bit of notation: I'll use the syntax

⊥@t

to mean an undefined value of type t, similar to Haskell's undefined.

∀a. [a]

How, then, would we write a term of type ∀a. [a]? We know that types containing ∀ come from terms containing Λ:

term1 :: ∀a. [a]
term1 = Λa. ?

Within the body marked ? we must provide a term of type [a]. However, we know nothing concrete about a, since it's an argument passed in from the outside. So we can return an empty list

term1 = Λa. []

or an undefined value

term1 = Λa. ⊥@[a]

or a list containing undefined values only

term1 = Λa. [⊥@a, ⊥@a]

but not much else.

To use this value, we apply a type, removing the outer ∀. Let's arbitrarily instantiate ∀a. [a] to [Bool]:

main = print @Int (length @Bool (term1 @Bool))

[∀a. a]

What about [∀a. a], then? If ∀ signifies a function on types, then [∀a. a] is a list of functions. We can provide as few as we like:

term2 :: [∀a. a]
term2 = []

or as many:

term2 = [f, g, h]

But what are our choices for f, g, and h?

f :: ∀a. a
f = Λa. ?

Now we're well and truly stuck. We have to provide a value of type a, but we know nothing whatsoever about the type a. So our only choice is

f = Λa. ⊥@a

So our options for term2 look like

term2 :: [∀a. a]
term2 = []
term2 = [Λa. ⊥@a]
term2 = [Λa. ⊥@a, Λa. ⊥@a]

etc. And let's not forget

term2 = ⊥@(∀a. [a])

Unlike the previous example, our choices for term2 are already lists, and we can pass them to length directly. As before, we have to pass the element type to length:

main = print @Int (length @(∀a. a) term2)

Existential types

So a value of universal (∀) type is a function from types to values. A value of existential (∃) type is a pair of a type and a value.

More specifically: A value of type

∃x. T

is a pair

(S, v)

where S is a type, and where v :: T, assuming you bind the type variable x to S within T.

Here's an existential type signature and a few terms with that type:

term3 :: ∃a. a
term3 = (Int,         3)
term3 = (Char,        'x')
term3 = (∀a. a → Int, Λa. λ(x::a). 4)

In other words, we can put any value we like into ∃a. a, as long as we pair that value with its type.

The user of a value of type ∀a. a is in a great position; they can choose any specific type they like, using the type application @T, to obtain a term of type T. The producer of a value of type ∀a. a is in a terrible position: they must be prepared to produce any type asked for, so (as in the previous section) the only choice is Λa. ⊥@a.

The user of a value of type ∃a. a is in a terrible position; the value inside is of some unknown specific type, not a flexible polymorphic value. The producer of a value of type ∃a. a is in a great position; they can stick any value they like into the pair, as we saw above.

So what's a less useless existential? How about values paired with a binary operation:

type Something = ∃a. (a, a → a → a, a → String)

term4_a, term4_b :: Something
term4_a = (Int,    (1,     (+)  @Int , show @Int))
term4_b = (String, ("foo", (++) @Char, λ(x::String). x))

Using it:

triple :: Something → String
triple = λ(a, (x :: a, f :: a→a→a, out :: a→String)).
  out (f (f x x) x)

The result:

triple term4_a  ⇒  "3"
triple term4_b  ⇒  "foofoofoo"

We've packaged up a type and some operations on that type. The user can apply our operations but cannot inspect the concrete value — we can't pattern-match on x within triple, since its type (hence set of constructors) is unknown. This is more than a bit like object-oriented programming.

Using existentials for real

The direct syntax for existentials using ∃ and type-value pairs would be quite convenient. UHC partially supports this direct syntax. But GHC does not. To introduce existentials in GHC we need to define new "wrapper" types.

Translating the above example:

{-# LANGUAGE ExistentialQuantification #-}

data Something = forall a. MkThing a (a -> a -> a) (a -> String)

term_a, term_b :: Something
term_a = MkThing 1 (+) show
term_b = MkThing "foo" (++) id

triple :: Something -> String
triple (MkThing x f out) = out (f (f x x) x)

There's a couple differences from our theoretical treatment. Type application, type abstraction, and type pairs are again implicit. Also, the wrapper is confusingly written with forall rather than exists. This references the fact that we're declaring an existential type, but the data constructor has universal type:

MkThing :: forall a. a -> (a -> a -> a) -> (a -> String) -> Something

Often, we use existential quantification to "capture" a type class constraint. We could do something similar here:

data SomeMonoid = forall a. (Monoid a, Show a) => MkMonoid a

main is usually a function

Continuations in C++ with fork

Small examples

Backtracking

Solving a maze

Excess backtracking

The implementation

Generating random functions

The data structure

Testing

How many functions?

See also

Writing kernel exploits

Embedding GDB breakpoints in C source code

Zombie 6.001 starts tomorrow!

Propane: Functional synthesis of images and animations in Haskell

Self-modifying code for debug tracing in quasi-C

The design

The TRACE macro

The linker script

Memory protection

Enabling a tracepoint

Pitfalls

Global locking through StablePtr

The C code

Foreign imports

The Haskell code

Inspecting the machine code

Haskell hackathon in the Boston area, January 20 to 22

Thunks and lazy blackholes: an introduction to GHC at runtime

The bug

Objects at runtime

Static objects

Black holes and revelations

unsafe[Dupable]PerformIO

The bug, finally

The fix

*deep breath*

Quasicrystals as sums of waves in the plane

The code

Running it

Interfacing Haskell to the Concorde solver for Traveling Salesperson Problem

Visualizing a tour

Safe top-level mutable variables for Haskell

Why global state?

Spooky action at a distance

Type clashes

phosphene: Fractal video feedback as a PC master boot record

Slides from "Why learn Haskell?"

shqq: Embedding shell commands in Haskell code

Finding duplicate files

debug-diff: Colorized diffs between Haskell values

Lambda to pi

The π-calculus

Compiling from the λ-calculus

Testing the compiler

Notes

New slides and how I made them

New slides

How I made the slides

x86 disassembly in Haskell with hdis86

Getting started

Analyzing instructions

Implementing breakpoints on x86 Linux

The API

The trap instruction

Invoking ptrace

Creating breakpoints

Running and stepping

Testing it

Inspecting the target

Notes

Type-level Fix and generic folds

Recursive types

Type-level Fix

Properties of Fix

Folding lists

Folding natural numbers

Folding rose trees

The definition of cata

The `TRACE` macro

`unsafe[Dupable]PerformIO`

deep breath

Invoking `ptrace`

Type-level `Fix`

Properties of `Fix`

The definition of `cata`