Things that amuse me

2011-07-09T18:54:00.000+01:00

Impredicative polymorphism, a use case

In a recent question on stackoverflow I made a comment about how Haskell can be considered a good imperative language because you can define abstractions to make it convenient. When I was going to make my point by implementing a simple example of it I found that what I wanted to do no longer works in ghc-7.0.4. Here's a simple example of what I wanted to do (which works in ghc-6.12.3). It's a simple library that makes Haskell code look a bit like Python.

{-# LANGUAGE ExtendedDefaultRules, TupleSections #-}
module Main where
import qualified Prelude
import Boa

last_elt(xs) = def $: do
    assert xs "Empty list"
    lst <- var xs              -- Create a new variable, lst
    ret <- var (xs.head)
    while lst $: do
        ret #= lst.head        -- Assign variable ret
        lst #= lst.tail
    return ret

first_elt(xs) = def $: do
    l <- var xs
    l.reverse                  -- Destructive reverse
    return (last_elt(l))

factorial(n) = def $: do
    assert (n<=0) "Negative factorial"
    ret <- var 1
    i <- var n
    while i $: do
        ret *= i
        i -= 1
    return ret

test = def $: do
    print "Hello"
    print ("factorial 10 =", factorial(10))

main = do
    test
    l <- varc [1, 2, 3]
    print ("first and last:",)
    print (first_elt(l),)
    print (last_elt(l))

On the whole it's pretty boring, except for one thing. In imperative languages there is (usually) a disctinction between l-values and r-values. An l-value represent a variable, i.e., something that can be assigned, whereas an r-value is simply a value. So in Python, in the statement x = x + 1 the x on the left is an l-value (it's being assigned), whereas on the right it's an r-value. You can use the same notation for both in most imperative languages. If you want to do the same in Haskell you have (at least) two choices. First you can unify the concepts of l-value and r-value and have a runtime test that you only try to assign variables. So, e.g., 5 = x would type check but have a runtime failure. I will not dwell further on this since it's not very Haskelly.
Instead, we want to have l-values and r-values being statically type checked. Here's the interesting bit of the types for a simple example.

data LValue a
data RValue a
instance (Num a) => Num (RValue a)

class LR lr
instance LR RValue
instance LR LValue

var :: RValue a -> IO (forall lr . (LR lr) => lr a)
(#=) :: LValue a -> RValue a -> IO ()

foo = do
    x <- var 42
    x #= x + 1

We have two type constructors LValue and RValue representing l-values and r-values of some type a. The r-values is an instance of Num. Furthermore, the class LR where the type is either LValue or RValue.
The var function creates a new variable given a value. The return type of var is the interesting part. It says that the return value is polymorphic; it can be used either as an l-value or as r-value. Just as we want.
The assignment operator, (#=), takes an l-value, an r-value, and returns nothing.

So in the example we expect x to have type forall lr . (LR lr) => lr a, in which case the assignment will type check.
If we try to compile this we get

    Illegal polymorphic or qualified type:
      forall (lr :: * -> *). LR lr => lr a
    Perhaps you intended to use -XImpredicativeTypes
    In the type signature for `var':
      var :: RValue a -> IO (forall lr. LR lr => lr a)

The problem is that universal quantification is not normally allowed as the argument to a type constructor. This requires the impredicative polymorphism extension to ghc. If we turn it on it compiles fine in ghc-6.12.3.
But, with ghc-7.0.4 we get

    Couldn't match expected type `LValue a0'
                with actual type `forall (lr :: * -> *). LR lr => lr a1'
    In the first argument of `(#=)', namely `x'
    In the expression: x #= x + 1

I can't really explain the rational behind the change in the ghc type system (Simon say it's simpler now), but I feel this has really made ghc worse for defining DSELs. When you define a DSEL you usually use the do-notation for the binding construct in the embedded language. This is the only programmable binding construct (by overloading (>>=)), so there is little choice. With the change in ghc-7 it means you can no longer make DSELs with a polymorhic binding construct (like I wanted here), because the binder seems to be monomorphic now

Please Simon, can we have polymorphic do bindings back?

2011-05-02T19:28:00.001+01:00

More points for lazy evaluation

In a recent blog post Bob Harper shows one use of laziness, but I think he misses the real import points of laziness. So I will briefly enumerate what I think are the important points of lazyness (excuse me for not coming up with any kind of pun about points).

First, I'd like to say that I don't think strict or lazy matters that much in practice; they both work fine. I have programmed in regular non-strict Haskell as well as strict Haskell, and most things carry over well to a strict Haskell. Furthermore, when I say strict language I mean a strict language with at least non-termination as an effect; total languages is another matter.

Lazy bindings

I like to tell Haskell beginners that any subexpression can be named and "pulled out", modulo name capture, of course. For instance

    ... e ...

is the same as

    let x = e
    in  ... x ...

The key thing is that

    ... e ... e ...

is the same as

    let x = e
    in  ... x ... x ...

so that common subexpressions can be given a name.

This is (in general) just wrong in a strict language, of course. Just take the simple example

    if c then error "BOO!" else 0

Which is not the same as

    let x = error "BOO!"
    in  if c then x else 0

In this case you can easily fix the problem by delaying the computation with a lambda (a common theme).

    let x () = error "BOO!"
    in  if c then x () else 0

But for a slightly more complicated example this simple technique goes wrong. Consider

    map (\ a -> a + expensive) xs

where expensive does not depend on a. In this case you want to move the expensive computation out of the loop (cf. loop invariant code motion in imperative languages). Like so

    let x = expensive
    in  map (\ a -> a + x) xs

In a lazy language x will be evaluated exactly zero times or once, just as we want. Using the delaying trick doesn't work here:

    let x () = expensive
    in  map (\ a -> a + x ()) xs

since expensive will get evaluated once for every list element.

This is easy enough to remedy, by introducing an abstraction for lazy computations (which will contain an assignment to compute the value just once). The signature for the abstract type of lazy values is something like

    data Lazy a
    delay :: (() -> a) -> Lazy a
    force :: Lazy a -> a

Note that the delay needs to take a function to avoid the a being evaluated early.
(This is probably what Bob would name a benign effect and is easily programmed using unsafePerformIO, which means it needs careful consideration.)

And so we get

    let x = delay (\ () -> expensive)
    in  map (\ a -> a + force x) xs

This isn't exactly pretty, but it works fine. In a language with macros the ugliness can be hidden better.

Lazy functions

Even strict languages like ML and C have some lazy functions even if they don't call them that, like SML's if, andthen, and orelse. You really need the if construct to evaluate the condition and then one of the branches depending on the condition. But what if I want to define my own type with the same kind of functions? In ML I can't, in C I would have to use a macro.

The ability to define new functions that can be used as control constructs is especially important when you want to design embedded domain specific languages. Take the simple example of the when (i.e., one-arm if) function in Haskell.

    when :: (Monad m) => Bool -> m () -> m ()

A quite common use of this function in monadic code is to check for argument preconditions in a function, like

    f x = do
        when (x < 0) $
            error "x must be >= 0"
        ...

If the when function is strict this is really bad, of course, since the call to error will happen before the when is called.

Again, one can work around this by using lazy values, like

    myAnd :: MyBool -> Lazy MyBool -> MyBool
    ...
    ... myAnd x (delay (\ () -> y)) ...

But in my opinion, this is too ugly to even consider. The intent of the function is obscured by the extra noise to make the second argument lazy.

I think every language needs a mechanism for defining some form of call-by-name functions. And many languages have this in the form of macros (maybe built in like in Lisp, or as a preprocessor like C).
If a language cannot define lazy function it simply lacks in abstraction power. I think any kind of fragment of the language should be nameable and reusable. (Haskell lacks this ability for patterns and contexts; a wart in the design.) So if you notice a repeated expression pattern, like

    if c then t else False

and cannot give this a name, like

    and c t = if c then t else False

and then use it with the same effect as the orginal expression, well, then your language is lacking.

For some language constructs the solution adopted by Smalltalk (and later Ruby), i.e., a very lightweight way on constructing closures is acceptable. So, for instance, I could accept writing

    ... myAnd x {y} ...

(In SML you could make something using functors, but it's just too ugly to contemplate.)

Lazy constructors

Lazy constructors were sort of covered in what Bob claimed to be the point of laziness, so I'll just mention them for completeness. Sometimes you need them, but in my experience it's not very often.

Cyclic data structures

This is related to the last point.

Sometimes you really want cyclic data structures. An example are the Haskell data types in Data.Data that describe data types and constructors. A data type descriptor needs to contain a list of its constructors and a constructor descriptor needs to contain the data type descriptor.
In Haskell this can be described very naturally by having the two descriptors reference each other.
In SML this is not possible. You will have to break the cycle by somthing like a reference (or a function).
In OCaml you can define cyclic data structures in a similar fashion to Haskell, so this isn't really a problem with strict languages, but rather a feature that you can have if you like. Of course, being able to define cyclic data leads to non-standard elements in your data types, like

    data Nat = Zero | Succ Zero
    omega :: Nat
    omega = Succ omega

So having the ability to define cyclic data structures is a double edged sword.
I find the lack of a simple way to define cyclic data a minor nuisance only.

Reuse

I've saved my biggest gripe of strict evaluation for last. Strict evaluation is fundamentally flawed for function reuse.
What do I mean? I will illustrate with and example.
Consider the any function is Haskell:

any :: (a -> Bool) -> [a] -> Bool
any p = or . map p

It's quite natural to express the any function by reusing the
map and or functions.  Unfortunately, it doesn't
behave like we would wish in a strict language.  The any function should scan the list from the head forwards and as soon as an
element that fulfills the predicate is found it should return true and stop
scanning the list.  In a strict language this would not happen, since
the predicate will be applied to every element before the or
examines the elements.

So we are forced to manually fuse the two functions, doing so we get:

any :: (a -> Bool) -> [a] -> Bool
any p = foldr False (\ x r -> p x || r)

or :: [Bool] -> Bool
or = foldr False (||)



foldr :: (a -> b -> b) -> b -> [a] -> b
foldr f z [] = z
foldr f z (x:xs) = f x (foldr f z xs)

But the misery doesn't end here. This still doesn't do the right thing, because the strict language will recurse all the way down the list since it will call foldr before f. So we either have to fuse again, or invent a new version of foldr that delays the recursive call.
One more fusion gets us to

any p []     = False
any p (y:ys) = y || any p ys

So where's the function reuse? Nowhere in sight.

With strict evaluation you can no longer with a straight face tell people: don't use recursion, reuse the recursion patterns in map, filter, foldr, etc. It simply doesn't work (in general).

Using macros doesn't really save us this time, because of the recursive definitions. I don't really know of any way to fix this problem short of making all (most?) functions lazy, because the problem is pervasive. I.e., in the example it would not be enough to fix foldr; all the functions involved need to be lazy to get the desired semantics.

I find this the biggest annoyance with strict evaluation, but at the same time it's just an annoyance, because you can always rewrite the code to make it work. But strict evaluation really, fundamentally stops you from reusing functions the way you can with laziness.

As an aside, the definition of any doesn't work in SML for another reason as well, namely the value restriction. But that's just a language wart, like the monomorphism restriction in Haskell.

Complexity

Another complaint Bob Harper had about lazy evaluation is the difficulty of finding out the complexity of functions. I totally agree that the space complexity is very messy in a lazy language. For (sequential) time complexity I don't see any need to worry.
If strict a function has O(f(n)) complexity in a strict language then it has complexity O(f(n)) in a lazy language as well. Why worry? :)

Summing up

One of the most important principles in all software design is DRY, i.e., Don't Repeat Yourself. This means that common patterns that we can find in a program should be abstractable so you can reuse them. In a lazy language you can use functions for abstracting control structures, which a strict language does not permit (at least not in a convenient way). This can be mitigated by providing some other abstraction mechanism, like macros (hopefully some kind of sane macros).
For the reasons above, I find it more convenient to program in a lazy language than a strict one. But a strict language with the ability to define lazy bindings, lazy functions, and lazy constructors is good enough that I find it usable.

2011-04-20T23:25:00.000+01:00

Ugly memoization

Here's a problem that I recently ran into. I have a function taking a string and computing some value. I call this function a lot, but a lot of the time the argument has occurred before. The function is reasonably expensive, about 10 us. Only about 1/5 of the calls to the function has a new argument.

So naturally I want to memoize the function. Luckily Hackage has a couple packages for memoization. I found data-memocombinators and MemoTrie and decided to try them. The basic idea with memoization is that you have a function like

  memo :: (a->b) -> (a->b)

I.e., you give a function to memo and you get a new function of the same type back. This new function behaves like the original one, but it remembers every time it is used and the next time it gets the same argument it will just return the remembered result. This is only safe in a pure language, but luckily Haskell is pure.
In an imperative language you can use a mutable memo table that stores all the argument-result pairs and updates the memo table each time the function is used. But how is it even possible to implement that in a pure language? The idea is to lazily construct the whole memo table in the call to memo, and it will then be lazily filled in.
Assume that all values of the argument type a can be enumerated by the method enumerate, we could then write memo like this:

  memo f =
      let table = [ (x, f x) | x <- enumerate ]
      in  \ y -> let Just r = lookup y table in r

Note how the memo table is constructed given just f, and this memo table is then used in the returned function.
The type of this function would be something like

  memo (Enumerate a, Eq a) => (a->b) -> (a->b)

assuming that the class Enumerate has the magic method enumerate.

This just a very simplified example, if you tried to use this it would be terrible because the returned function does linear lookup in a list. Instead we want some kind of search tree, which is what the two packages I mention implement. The MemoTrie package does this in a really beautiful way, I recommend reading Conal's blog post about it.
OK, enough preliminaries. I used criterion to perform the benchmarking, and I tried with no memoization (none), memo-combinators (comb), and MemoTrie (beau). I had a test function taking about 10us, and then i called this functions with different number of repeated arguments: 1, 2, 5, and 10. I.e., 5 means that each argument occurred 5 times as the memoized function was called.

	1	2	5	10
none	10.7	10.7	10.7	10.7
comb	62.6	52.2	45.8	43.4
beau	27.6	17.0	10.4	8.1

So with no memoization the time per call was 10.7 us all the time, no surprise there. With the memo combinators it was much slower than no memoization; the overhead for looking something up is bigger than the cost of computing the result. So that was a failure. The MemoTrie does better, at about an argument repetition of five it starts to break even, and at ten it's a little faster to memoize.

Since I estimated my repetition factor in the real code to be about five even the fastest memoization would not be any better then recomputation. So now what? Give up? Of course not! It's time to get dirty.

Once you know a function can be implemented in a pure way, there's no harm in implementing the same function in an impure way as long as it presents the pure interface. So lets write the memo function the way it would be done in, e.g., Scheme or ML. We will use a reference to hold a memo table that gets updated on each call. Here's the code, with the type that the function gets.

import Data.IORef
import qualified Data.Map as M

memoIO :: (Ord a) => (a -> b) -> IO (a -> IO b)
memoIO f = do
  v <- newIORef M.empty
  let f' x = do
        m <- readIORef v
        case M.lookup x m of
          Nothing -> do let { r = f x }; writeIORef v (M.insert x r m); return r
          Just r  -> return r
  return f'

The memoIO allocated a reference with an empty memo table. We then define a new function, f', which when it's called with get the memo table and look up the argument. If the argument is in the table then we just return the result, if it's not then we compute the result, store it in the table, and return it. Good old imperative programming (see below why this code is not good imperative code).

But, horror, now the type is all wrong, there's IO in two places. The function we want to implement is actually pure. So what to do? Well, if you have a function involving the IO type, but you can prove it is actually pure, then (and only then) you are allowed to use unsafePerformIO.

I'll wave my hands instead of a proof (but more later), and here we go

  memo :: (Ord a) => (a -> b) -> (a -> b)
  memo f = let f' = unsafePerformIO (memoIO f) in \ x -> unsafePerformIO (f' x)

Wow, two unsafePerformIO on the same line. It doesn't get much less safe than that.
Let's benchmark again:

	1	2	5	10
none	10.7	10.7	10.7	10.7
comb	62.6	52.2	45.8	43.4
beau	27.6	17.0	10.4	8.1
ugly	13.9	7.7	3.9	2.7

Not too shabby, using the ugly memoization is actually a win already at two, and just a small overhead if the argument occurs once. We have a winner!

No so fast, there's

A snag

My real code can actually be multi-threaded, so the memo function had better work in a multi-threaded setting. Well, it doesn't. There's no guarantee about readIORef and writeIORef when doing multi-threading.
So we have to rewrite it. Actually, the code I first wrote is the one below; I hardly ever use IORef because I want it to work with multi-threading.

memoIO f = do
  v <- newMVar M.empty
  let f' x = do
        m <- takeMVar v
        case M.lookup x m of
          Nothing -> do let { r = f x }; putMVar v (M.insert x r m); return r
          Just r  -> do                  putMVar v m;                return r
  return f'

So now we use an MVar instead. This makes it thread safe. Only one thread can execute between the takeMVar and the putMVar. This guarantees than only one thread can update the memo table at a time. If two threads try at the same time one has to wait a little. How long? The time it takes for the lookup, plus some small constant. Remember that Haskell is lazy, the the (f x) is not actually computed with the lock held, which is good.
So I think this is a perfectly reasonable memoIO. And we can do the same unsafe trick as before and make it pure. Performance of this version is the same as with the IORef.

Ahhhh, bliss. But wait, there's

Another snag

That might look reasonable, but in fact the memo function is broken now. It appears to work, but here's a simple use that fails

  sid :: String ->; String
  sid = memo id

  fcn s = sid (sid s)

What will happen here? The outer call to sid will execute the takeMVar and then do the lookup. Doing the lookup with evaluate the argument, x. But this argument is another call to sid, this will try to execute the takeMVar. Disaster has struck, deadlock.

What happened here? The introduction of unsafePerformIO ruined the sequencing guaranteed by the IO monad that would have prevented the deadlock if we had used memoIO. I got what I deserved for using unsafePerformIO.

Can it be repaired? Well, we could make sure x is fully evaluated before grabbing the lock. I settled for a different repair, where the lock is held in a shorter portion of the code.

memoIO f = do
  v <- newMVar M.empty
  let f' x = do
        m <- readMVar v
        case M.lookup x m of
          Nothing -> do let { r = f x }; m <- takeMVar v; putMVar v (M.insert x r m); return r
          Just r  -> return r
  return f'

This solution has its own problem. It's now possible for several threads to compute (f x) for the same x and the result of all but one of those will be lost by overwriting the table. This is a price I'm willing to pay for this application.

Moral

Yes, you can use imperative programming to implement pure functions. But the onus is on you to prove that it is safe. This is not as easy as you might think. I believe my final version is correct (with the multiple computation caveat), but I'm not 100% sure.

2011-04-10T19:51:00.002+01:00

Phew! Cleaned out a lot of spam comments in my blog. Hopefully my new settings will prevent the crazy onslaught of spammers.

2009-06-10T22:35:00.002+01:00

More LLVM

Recently someone asked me on #haskell if you could use the Haskell LLVM bindings to compile some abstract syntax to a Haskell function. Naturally I said yes, but then I realized I had only done it for a boring language with just one type. I had no doubt that it could be done for more complicated languages with multiple types, but it might not be totally obvious how. So I decided to write a simple compiler, and this blog post is the result. First, a simple example:

    main = do
        let f :: Double -> Double
            Just f = compile "\\ (x::Double) -> if x == 0 then 0 else 1/(x*x)"
        print (f 2, f 3, f 0)

Running this program produces (as expected)

    (0.25,0.1111111111111111,0.0)

What has happened is that the string has been parsed to an abstract syntax tree, translated into LLVM code, then to machine code, and finally turned back into a Haskell callable function. Many things can go wrong along the way, like syntax and type errors, so compile returns a Maybe type to indicate if things went right or wrong. (A more serious version of the compile function would return an error message when something has gone wrong.) The definition of the compilation function is simple and illustrates the flow of the compiler

compile :: (Translate a) => String -> Maybe a
compile = fmap translate . toTFun <=< mParseUFun

The context Translate is there to limit the types that can actually be translated; it's a necessary evil and exactly what types are allowed depends on how advanced we make the compiler. Had we ignored the Maybe type the definitions would have been

compile = translate . toTFun . mParseUFun

which says, first parse to the type UFun (untyped expressions), then type check and turn it into the type TFun a, and finally translate TFun a into an a by LLVM compilation. Let's see how this all works.

The UExp module

The first step is to just define an abstract syntax for the expressions that we want to handle. I'm only allowing leading lambdas (this a very first order language), so there's a distiction between the top level UFun type and the expression type UExp. The U prefix indicates that this version of the syntax is not yet type checked. The definition is pretty boring, but here it is:

{-# OPTIONS_GHC -fno-warn-unused-binds #-}
{-# LANGUAGE RecordWildCards #-}
module UExp(Id, UFun(..), UTyp(..), UExp(..), parseUFun, showOp, mParseUFun) where
import Data.Maybe
import Data.List
import Data.Function
import Text.ParserCombinators.Parsec
import Text.ParserCombinators.Parsec.Expr
import Text.ParserCombinators.Parsec.Token
import Text.ParserCombinators.Parsec.Language

type Id = String

data UFun = UFun [(Id, UTyp)] UExp

data UTyp = UTBol | UTDbl

data UExp
    = UDbl Double        -- ^ Double literal
    | UBol Bool          -- ^ Bool literal
    | UVar Id            -- ^ Variable
    | UApp Id [UExp]     -- ^ Function application
    | ULet Id UExp UExp  -- ^ Local binding

Naturally, we want to be able to show the expressions, if nothing else so for debugging. So I make a Show instance that shows them in a nice way respecting operator precedences etc. There's nothing exciting going on, the large number of lines is just to cover operator printing.

instance Show UFun where
    showsPrec p (UFun [] e) = showsPrec p e
    showsPrec p (UFun vts e) = showParen (p>0) (showString "\\ " . foldr (.) (showString "-> ") (map f vts) . showsPrec 0 e)
      where f (v, t) = showParen True (showString v . showString " :: " . showsPrec 0 t) . showString " "

instance Show UTyp where
    showsPrec _ UTDbl = showString "Double"
    showsPrec _ UTBol = showString "Bool"

instance Show UExp where
    showsPrec p (UDbl d) = showsPrec p d
    showsPrec p (UBol b) = showsPrec p b
    showsPrec _ (UVar i) = showString i
    showsPrec p (UApp "if" [c, t, e]) =
      showParen (p>0) (showString "if " . showsPrec 0 c . showString " then " . showsPrec 0 t . showString " else " . showsPrec 0 e)
    showsPrec p (UApp op [a, b]) = showOp p op a b
    showsPrec _ (UApp op _) = error $ "Uxp.show " ++ op
    showsPrec p (ULet i e b) =
      showParen (p>0) (showString "let " . showString i . showString " = " . showsPrec 0 e . showString " in " . showsPrec 0 b)

showOp :: (Show a, Show b) => Int -> String -> a -> b -> String -> String
showOp q sop a b = showParen (q>mp) (showsPrec lp a . showString sop . showsPrec rp b)
  where (lp,mp,rp) = case lookup sop ops of
                    Just (p, AssocLeft)  -> (p,   p, p+1)
                    Just (p, AssocRight) -> (p+1, p, p)
                    Just (p, AssocNone)  -> (p+1, p, p+1)
                    Nothing              -> (9,   9, 10)

ops :: [(String, (Int, Assoc))]
ops = [("+",  (6, AssocLeft)),
       ("-",  (6, AssocLeft)),
       ("*",  (7, AssocLeft)),
       ("/",  (7, AssocLeft)),
       ("==", (4, AssocNone)),
       ("<=", (4, AssocNone)),
       ("&&", (3, AssocRight)),
       ("||", (2, AssocRight))
      ]

We also want to be able to parse, so I'm using Parsec to parse the string and produce an AST. Again, there's nothing interesting going on. I use the Haskell lexical analysis provided by Parsec. This is available as a TokenParser record, which can be conveniently opened with the RecordWildcard notation TokenParser{..}.

parseUFun :: SourceName -> String -> Either ParseError UFun
parseUFun = parse $ do f <- pFun; eof; return f
  where TokenParser{..} = haskell
        pFun = do
            vts <- between (reservedOp "\\")
                           (reservedOp "->")
                           (many $ parens $ do v <- identifier; reservedOp "::"; t <- pTyp; return (v, t))
               <|> return []
            e <- pExp
            return $ UFun vts e
        pTyp = choice [do reserved "Bool"; return UTBol, do reserved "Double"; return UTDbl]
        pExp = choice [pIf, pLet, pOExp]
        pIf = do reserved "if"; c <- pExp; reserved "then"; t <- pExp; reserved "else"; e <- pExp; return $ UApp "if" [c, t, e]
        pLet = do reserved "let"; i <- identifier; reservedOp "="; e <- pExp; reserved "in"; b <- pExp; return $ ULet i e b
        pOExp = buildExpressionParser opTable pAExp
        pAExp = choice [pDbl, pVar, parens pExp]
        pVar = fmap eVar identifier
        pDbl = fmap (either (UDbl . fromInteger) UDbl) naturalOrFloat
        eVar i = if i == "False" then UBol False else if i == "True" then UBol True else UVar i

        opTable = reverse . map (map mkOp) . groupBy ((==) `on` prec) . sortBy (compare `on` prec) $ ops
          where mkOp (s, (_, a)) = Infix (do reservedOp s; return $ \ x y -> UApp s [x, y]) a
                prec = fst . snd

mParseUFun :: String -> Maybe UFun
mParseUFun = either (const Nothing) Just . (parseUFun "")

The parser is packaged up in mParseUFun which returns an AST if it all worked.

The TExp module

Since the LLVM API is typed it's much easier to translate a typed abstract syntax tree than an untyped abstract syntax tree. The TExp module contains the definition of the typed AST and the type checker that converts to it. There are many ways to formulate type safe abstract syntax trees. I've chosen to use GADTs. I've also picked to represent variables (still) by identifiers, which means that the syntax tree is not necessarily type safe. Furthermore, I've chosen a very limited way to represent function application since this is all I need for this example. The variantions on this are endless.

{-# LANGUAGE GADTs, ExistentialQuantification, PatternGuards #-}
module TExp(Id,
            TFun(..), TTyp(..), TExp(..), DblOp(..), BolOp(..), CmpOp(..),
            Equal(..), test,
            Type(..),
            AFun(..), extractFun,
            typeCheck, toTFun) where
import Data.Maybe
import Control.Monad
import UExp

data TFun a where
    TBody :: TExp a                 -> TFun a
    TLam  :: Id -> TTyp a -> TFun b -> TFun (a->b)

data TTyp a where
    TTBol ::                     TTyp Bool
    TTDbl ::                     TTyp Double
    TTArr :: TTyp a -> TTyp b -> TTyp (a->b)

data TExp a where
    TDbl   :: Double                                            -> TExp Double
    TBol   :: Bool                                              -> TExp Bool
    TDblOp :: DblOp     -> TExp Double -> TExp Double           -> TExp Double
    TBolOp :: BolOp     -> TExp Bool   -> TExp Bool             -> TExp Bool
    TCmpOp :: CmpOp     -> TExp Double -> TExp Double           -> TExp Bool
    TIf    :: TExp Bool -> TExp a      -> TExp a                -> TExp a
    TLet   :: Id        -> TTyp a      -> TExp a      -> TExp b -> TExp b
    TVar   :: Id                                                -> TExp a

data DblOp = DAdd | DSub | DMul | DDiv
    deriving (Eq, Show)

data BolOp = BAnd | BOr
    deriving (Eq, Show)

data CmpOp = CEq | CLe
    deriving (Eq, Show)

So for instance, UApp "+" [UVar "x", UDbl 2.2] will be represented by TDblOp DAdd (TVar "x") (TDbl 2.2) which has type TExp Double. So the type of the expression is now accurately reflected in the type of the syntax tree. Even the UTyp type now has a typed equivalent where the real type is reflected. For completeness, here's some code for pretty printing etc.

{-# LANGUAGE GADTs, ExistentialQuantification, PatternGuards #-}
module TExp(Id,
            TFun(..), TTyp(..), TExp(..), DblOp(..), BolOp(..), CmpOp(..),
            Equal(..), test,
            Type(..),
            AFun(..), extractFun,
            typeCheck, toTFun) where
import Data.Maybe
import Control.Monad
import UExp

instance Show (TFun a) where
    showsPrec p (TBody e) = showsPrec p e
    showsPrec p (TLam i t e) = showParen (p>0)
      (showString "\\ " . showParen True (showString i . showString " :: " . showsPrec 0 t) . showString " -> " . showsPrec 0 e)

instance Show (TTyp a) where
    showsPrec _ TTBol = showString "Bool"
    showsPrec _ TTDbl = showString "Double"
    showsPrec p (TTArr a b) = showParen (p>5) (showsPrec 6 a . showString " -> " . showsPrec 5 b)

instance Show (TExp a) where
    showsPrec p (TDbl d) = showsPrec p d
    showsPrec p (TBol b) = showsPrec p b
    showsPrec _ (TVar i) = showString i
    showsPrec p (TDblOp op a b) = showOp p (fromJust $ lookup op [(DMul, "*"), (DAdd, "+"), (DSub, "-"), (DDiv, "/")]) a b
    showsPrec p (TBolOp op a b) = showOp p (fromJust $ lookup op [(BAnd, "&&"), (BOr, "||")]) a b
    showsPrec p (TCmpOp op a b) = showOp p (fromJust $ lookup op [(CEq, "=="), (CLe, "<=")]) a b
    showsPrec p (TIf c t e) = showParen (p>0) (showString "if " . showsPrec 0 c . showString " then " . showsPrec 0 t . showString " else " . showsPrec 0 e)
    showsPrec p (TLet i _ e b) =
      showParen (p>0) (showString "let " . showString i . showString " = " . showsPrec 0 e . showString " in " . showsPrec 0 b)

The aim of the type checker is to transform from the UExp type to the TExp type, so basically

typeCheckExp :: UExp -> TExp a

But things can go wrong, so it's impossible to always return a TExp, so let's use a Maybe type:

typeCheckExp :: UExp -> Maybe (TExp a)

But wait! This type is totally wrong. Why? Because it promises that given a UExp the type checker can return any type, i.e., writing out the (normally implicit) quantifier the type is:

typeCheckExp :: forall a . UExp -> Maybe (TExp a)

But this is not the case, the type checker will figure out a type and return an expression with this specific type, so the type we really want is

typeCheckExp :: exists a . UExp -> Maybe (TExp a)

Haskell doesn't allow this type to be written this way; we need to package up the existential type in a data type. Like so:

data ATExp = forall a . TExp a ::: TTyp a

data AFun = forall a . AFun (TFun a) (TTyp a)

It might look funny that the existential type is written with a forall, but it makes sense when looking at the type of the constructor function (but not when doing pattern matching). Now we can attempt a couple of cases of the type checker:

typeCheckExp :: UExp -> Maybe ATExp
typeCheckExp (UDbl d) =
    return $ TDbl d ::: TTDbl
typeCheckExp (UBol b) =
    return $ TBol b ::: TTBol

They look quite nice, and they actually work. So what about something more complicated, like arithmetic?

typeCheckExp (UApp op [a, b]) | Just dop <- lookup op [("+", DAdd), ("-", DSub), ("*", DMul), ("/", DDiv)] = do
    a' ::: TTDbl <- typeCheckExp a
    b' ::: TTDbl <- typeCheckExp b
    return $ TDblOp dop a' b' ::: TTDbl

First we conveniently look up the operator among the arithmetic operators, then we recursively call the type checker for the operands. We do this in the Maybe monad. If the type checking a subterm fails that's automatically propagated, and furthermore, if the type checking of a subterm does not yield a TTDbl type then this will cause the pattern matching to fail, and this will generate a Nothing in the maybe monad, so we used failing pattern matching to our advantage here. The interesting case is checking UIf, because here both arms have to have the same type, but we don't know which one. Here's an attempt:

typeCheckExp (UApp "if" [c,t,e]) = do
    c' ::: TTBol <- typeCheckExp c
    t' ::: tt    <- typeCheckExp t
    e' ::: te    <- typeCheckExp e
    guard (tt == te)
    return $ TIf c' t' e' ::: tt

But this doesn't type check. The guard ensures that the two arms have the same type, but that's something we know, but the Haskell type checker doesn't. So it rejects the TIf, because ut can't see that both arms have the same type. We need to be trickier in doing the equality test so that it reflects the equality on the type level. There's a standard trick for this, namely this type:

data Equal a b where
    Eq :: Equal a a

If you ever have a value (which must be Eq) of type Equal foo bar then the type checker will know that foo and bar are actually the same type. So let's code equality for TTyp.

test :: TTyp a -> TTyp b -> Maybe (Equal a b)
test TTBol TTBol = return Eq
test TTDbl TTDbl = return Eq
test (TTArr a b) (TTArr a' b') = do
    Eq <- test a a'
    Eq <- test b b'
    return Eq
test _ _ = mzero

This code is worth pondering for a while, it's actually rather clever (I take no credit for it; I stole it from Tim Sheard). Why does even the first clause type check? Because TTBol has type TTyp Bool, so both the type variables (a and b) must be TTBool in the first clause, which means that Eq :: Equal TBol TBol is what we're returning. Equipped with this equality we can try type checking again.

typeCheckExp (UApp "if" [c,t,e]) = do
    c' ::: TTBol <- typeCheckExp c
    t' ::: tt    <- typeCheckExp t
    e' ::: te    <- typeCheckExp e
    Eq <- test tt te
    return $ TIf c' t' e' ::: tt

And amazingly this actually works! (A tribute to the hard working ghc implementors.) One (rather large) fly is left in the ointment. What about variables? What do we do when we type check UVar? We must check that there's a bound variable with the right type around. So the type checker needs to be extended with an environment where variables can be looked up. It's mostly straight forward. The environment simply maps a variable to ATExp. So here's the complete type checker as it's actually defined.

type Env = [(Id, ATExp)]

typeCheckExp :: Env -> UExp -> Maybe ATExp
typeCheckExp _ (UDbl d) =
    return $ TDbl d ::: TTDbl
typeCheckExp _ (UBol b) =
    return $ TBol b ::: TTBol
typeCheckExp r (UApp op [a, b]) | Just dop <- lookup op [("+", DAdd), ("-", DSub), ("*", DMul), ("/", DDiv)] = do
    a' ::: TTDbl <- typeCheckExp r a
    b' ::: TTDbl <- typeCheckExp r b
    return $ TDblOp dop a' b' ::: TTDbl
typeCheckExp r (UApp op [a, b]) | Just bop <- lookup op [("&&", BAnd), ("||", BOr)] = do
    a' ::: TTBol <- typeCheckExp r a
    b' ::: TTBol <- typeCheckExp r b
    return $ TBolOp bop a' b' ::: TTBol
typeCheckExp r (UApp op [a, b]) | Just cop <- lookup op [("==", CEq), ("<=", CLe)] = do
    a' ::: TTDbl <- typeCheckExp r a
    b' ::: TTDbl <- typeCheckExp r b
    return $ TCmpOp cop a' b' ::: TTBol
typeCheckExp r (UApp "if" [c,t,e]) = do
    c' ::: TTBol <- typeCheckExp r c
    t' ::: tt    <- typeCheckExp r t
    e' ::: te    <- typeCheckExp r e
    Eq <- test tt te
    return $ TIf c' t' e' ::: tt
typeCheckExp r (ULet i e b) = do
    e' ::: te <- typeCheckExp r e
    b' ::: tb <- typeCheckExp ((i, TVar i ::: te) : r) b
    return $ TLet i te e' b' ::: tb
typeCheckExp r (UVar i) =
    lookup i r
typeCheckExp _ _ =
    mzero

Note the ULet case which extends the environment. First we type check the expression that's being bound, and then add a variable to the environment and type check the body. Finally we need to type check the top level:

typeCheck :: UFun -> Maybe AFun
typeCheck = typeCheckFun []

typeCheckFun :: Env -> UFun -> Maybe AFun
typeCheckFun n (UFun [] b) = do
    e ::: t <- typeCheckExp n b
    return $ AFun (TBody e) t
typeCheckFun n (UFun ((x, typ):vts) b) =
    case typ of
    UTBol -> f TTBol
    UTDbl -> f TTDbl
  where f t = do AFun e r <- typeCheckFun ((x, TVar x ::: t) : n) (UFun vts b); return $ AFun (TLam x t e) (TTArr t r)

When encountering the expression we just type check it, and for an argument we add a variable with the right type to the environment. A small test in ghci:

TExp UExp> mParseUFun "\\ (x::Double) -> x+1" >>= typeCheck
Just (\ (x :: Double) -> x+1.0 :: Double -> Double)

To be able to extract a function from ATFun we need some small utilties.

class Type a where
    theType :: TTyp a
instance Type Double where
    theType = TTDbl
instance Type Bool where
    theType = TTBol
instance (Type a, Type b) => Type (a->b) where
    theType = TTArr theType theType

extractFun :: (Type a) => AFun -> Maybe (TFun a)
extractFun = extract theType

extract :: TTyp a -> AFun -> Maybe (TFun a)
extract s (AFun e t) = do
    Eq <- test t s
    return e

toTFun :: (Type a) => UFun -> Maybe (TFun a)
toTFun = extractFun <=< typeCheck

The class Type allows us to construct the TTyp corresponding to a Haskell type via overloading. Using this and the test function we can then extract a TFun at any type we like. If we try to extract at the wrong type we'll just get Nothing and at the right type we get Just.

The Compiler module

Now all we need to do is to write a function translate that translates a TFun a into the corresponding a. Naturally, using LLVM. Let's start with some simple cases in translating literals to LLVM code.

compileExp :: TExp a -> CodeGenFunction r (Value a)
compileExp (TDbl d) = return $ valueOf d
compileExp (TBol b) = return $ valueOf b

The valueOf function is simply the one that lifts a Haskell value into an LLVM value. Note how nice the GADT works out here and we handle both Double and Bool with any need to compromise type safety. What about arithmetic? Equally easy.

compileExp r (TDblOp op e1 e2) = bind2 (dblOp op) (compileExp r e1) (compileExp r e2)

dblOp :: DblOp -> Value Double -> Value Double -> CodeGenFunction r (Value Double)
dblOp DAdd = add
dblOp DSub = sub
dblOp DMul = mul
dblOp DDiv = fdiv

-- This should be in Control.Monad
bind2 :: (Monad m) => (a -> b -> m c) -> m a -> m b -> m c
bind2 f m1 m2 = do
    x1 <- m1
    x2 <- m2
    f x1 x2

And we can just carry on:

compileExp (TBolOp op e1 e2) = bind2 (bolOp op) (compileExp e1) (compileExp e2)
compileExp (TCmpOp op e1 e2) = bind2 (cmpOp op) (compileExp e1) (compileExp e2)
compileExp (TIf b t e) = mkIf (compileExp b) (compileExp t) (compileExp e)

bolOp :: BolOp -> Value Bool -> Value Bool -> CodeGenFunction r (Value Bool)
bolOp BAnd = and
bolOp BOr  = or

cmpOp :: CmpOp -> Value Double -> Value Double -> CodeGenFunction r (Value Bool)
cmpOp CEq = fcmp FPOEQ
cmpOp CLe = fcmp FPOLE

(The && and || are not short circuiting in this implementation. It would be easy to change.) It's rather amazing that despite these different branches producing and consuming different types it all works out. It's perfectly type safe and free from coercions. This is the beauty of GADTs. Oh, yeah, mkIf. It's just a piece of mess to create some basic blocks, test, and jump.

mkIf :: (IsFirstClass a) =>
        CodeGenFunction r (Value Bool) -> CodeGenFunction r (Value a) -> CodeGenFunction r (Value a) -> CodeGenFunction r (Value a)
mkIf mb mt me = do
    b <- mb
    tb <- newBasicBlock
    eb <- newBasicBlock
    jb <- newBasicBlock
    condBr b tb eb
    defineBasicBlock tb
    t <- mt
    br jb
    defineBasicBlock eb
    e <- me
    br jb
    defineBasicBlock jb
    phi [(t, tb), (e, eb)]

OK, so was lying. The translate function is not quite as easy as that. Just as with type checking we need an environment because of variables. It's easy to add though, and here's the real code.

compileExp :: (Type a, IsFirstClass a) => Env -> TExp a -> CodeGenFunction r (Value a)
compileExp _ (TDbl d) = return $ valueOf d
compileExp _ (TBol b) = return $ valueOf b
compileExp r (TDblOp op e1 e2) = bind2 (dblOp op) (compileExp r e1) (compileExp r e2)
compileExp r (TBolOp op e1 e2) = bind2 (bolOp op) (compileExp r e1) (compileExp r e2)
compileExp r (TCmpOp op e1 e2) = bind2 (cmpOp op) (compileExp r e1) (compileExp r e2)
compileExp r (TIf b t e) = mkIf (compileExp r b) (compileExp r t) (compileExp r e)
compileExp r (TLet i t e b) = do
    e' <- compileExp' r t e
    compileExp ((i, AValue e' t):r) b
compileExp r (TVar i) = return $ fromJust $ castAValue theType =<< lookup i r   -- lookup cannot fail on type checked code

compileExp' :: Env -> TTyp a -> TExp a -> CodeGenFunction r (Value a)
compileExp' r TTDbl e = compileExp r e
compileExp' r TTBol e = compileExp r e
compileExp' _ _ _ = error $ "compileExp': functions not allowed yet"

data AValue = forall a . AValue (Value a) (TTyp a)

castAValue :: TTyp a -> AValue -> Maybe (Value a)
castAValue t (AValue v s) = do
    Eq <- test t s
    return v

type Env = [(Id, AValue)]

Exactly as for the type checking environment we stick the code generation in an environment, and use castAValue project it out of the existential container. The fromJust call in the TVar case cannot fail on type checked code, but with my string based variable representation I have no evidence of this in the TExp so there's actually a cast in the variable case that can fail if scope and type checking has not been performed. The compileExp' is placate the type checker and help it with some evidence about that we are only binding base values. The rest of the code generation module is just house keeping. It's a little ugly, but not terrible.

-- | Compile a TFun into the corresponding LLVM code.
compileFunction :: (Translate a) =>
                   TFun a -> CodeGenModule (Function (RetIO a))
compileFunction = createFunction ExternalLinkage . compileFun []

class Compile a where
    type CG a
    type RetIO a
    type Returns a
    compileFun :: Env -> TFun a -> CG a

instance Compile Double where
    type CG Double = CodeGenFunction Double ()
    type RetIO Double = IO Double
    type Returns Double = Double
    compileFun r (TBody e) = compileExp r e >>= ret
    -- TLam is not well typed

instance Compile Bool where
    type CG Bool = CodeGenFunction Bool ()
    type RetIO Bool = IO Bool
    type Returns Bool = Bool
    compileFun r (TBody e) = compileExp r e >>= ret
    -- TLam is not well typed

instance (Type a, Compile b) => Compile (a -> b) where
    type CG (a->b) = Value a -> CG b
    type RetIO (a->b) = a -> RetIO b
    type Returns (a->b) = Returns b
    -- TBody is not well typed
    compileFun r (TLam i t e) = \ x -> compileFun ((i, AValue x t):r) e

The verbosity and large number of type functions in this section has convinced me that I need to simplify some of the types and classes involved in the LLVM code generation. To convert and LLVM module we call the JIT. This produces a function that returns a value in the IO monad (to be on the safe side) so we need to get rid of the IO, and finally we can get rid of the top level IO, because externally what we are doing is really pure (in some sense).

translate :: (Translate a) => TFun a -> a
translate = unsafePerformIO . fmap unsafePurify . simpleFunction . compileFunction

The Translate context is just an abbreviation for a big context enforced by the LLVM functions. It looks horrendous, but the type checker figured it out for me and I just pasted it in.

{-# LANGUAGE TypeFamilies, FlexibleContexts, ExistentialQuantification, FlexibleInstances, UndecidableInstances #-}
module Compile(Translate, translate) where
import Data.Maybe
import Prelude hiding (and, or)
import TExp
import LLVM.Core hiding (CmpOp)
import LLVM.ExecutionEngine
import System.IO.Unsafe(unsafePerformIO)

class    (Type a,
          Unsafe (RetIO a) a,
          FunctionArgs (RetIO a) (CG a) (CodeGenFunction (Returns a) ()),
          IsFunction (RetIO a),
          Compile a,
          Translatable (RetIO a)) =>
    Translate a
instance (Type a,
          Unsafe (RetIO a) a,
          FunctionArgs (RetIO a) (CG a) (CodeGenFunction (Returns a) ()),
          IsFunction (RetIO a),
          Compile a,
          Translatable (RetIO a)) =>
    Translate a

Conclusion

And that concludes the three parts of the compiler. In about 400 lines of code we can compile a small subset of Haskell expressions to (efficient) machine code. After type checking the rest of the processing is done in a type safe manner (except for a cast in TVar) which is the intention of the high level LLVM interface. Oh, and if you instrument the code generator a little you can peek at the machine code being produced. For instance, for this input to compile

\ (x::Double) ->
let y = x*(x-1) in
let z = x/y + 1 in
if y <= 0 then 0 else 1/(y-z)

we get

__fun1:
 subl $12, %esp
 movsd LCPI1_0, %xmm0
 movsd 16(%esp), %xmm1
 movapd %xmm1, %xmm2
 subsd %xmm0, %xmm2
 mulsd %xmm1, %xmm2
 pxor %xmm3, %xmm3
 ucomisd %xmm2, %xmm3
 jae LBB1_3
LBB1_1:
 divsd %xmm2, %xmm1
 addsd %xmm0, %xmm1
 subsd %xmm1, %xmm2
 movsd LCPI1_0, %xmm0
 divsd %xmm2, %xmm0
LBB1_2:
 movsd %xmm0, (%esp)
 fldl (%esp)
 addl $12, %esp
 ret
LBB1_3:
 pxor %xmm0, %xmm0
 jmp LBB1_2

2009-02-07T20:31:00.003Z

More BASIC

Not that anybody should care, but I've reimplemented by BASIC.

Here's a simple program.

{-# LANGUAGE ExtendedDefaultRules, OverloadedStrings #-}
import BASIC

main = runBASIC $ do
    10 GOSUB 1000
    20 PRINT "* Welcome to HiLo *"
    30 GOSUB 1000

    100 LET I := INT(100 * RND(0))
    200 PRINT "Guess my number:"
    210 INPUT X
    220 LET S := SGN(I-X)
    230 IF S <> 0 THEN 300

    240 FOR X := 1 TO 5
    250   PRINT X*X;" You won!"
    260 NEXT X
    270 STOP

    300 IF S <> 1 THEN 400
    310 PRINT "Your guess ";X;" is too low."
    320 GOTO 200

    400 PRINT "Your guess ";X;" is too high."
    410 GOTO 200

    1000 PRINT "*******************"
    1010 RETURN

    9999 END

In some ways this is a step backwards, since it requires some language extensions in Main. But I wanted to be able to use semicolon in the print statement.

But there it is, an exciting game!

*******************
* Welcome to HiLo *
*******************
Guess my number:
50
Your guess 50 is too high.
Guess my number:
25
Your guess 25 is too low.
Guess my number:
37
Your guess 37 is too low.
Guess my number:
44
Your guess 44 is too low.
Guess my number:
47
Your guess 47 is too low.
Guess my number:
48
1 You won!
4 You won!
9 You won!
16 You won!
25 You won!

2009-02-06T16:50:00.002Z

Is Haskell fast?

Let's do a simple benchmark comparing Haskell to C. My benchmark computes an approximation to infinity by adding up 1/n. Here is the C code:

#include <stdio.h>

int
main(int argc, char **argv)
{
  double i, s;
  s = 0;
  for (i = 1; i < 100000000; i++)
    s += 1/i;
  printf("Almost infinity is %g\n", s);
}

And running it

Lennarts-Computer% gcc -O3 inf.c -o inf
Lennarts-Computer% time ./inf
Almost infinity is 18.9979
1.585u 0.009s 0:01.62 97.5%     0+0k 0+0io 0pf+0w

And now the Haskell code:

import BASIC

main = runBASIC' $ do

    10 LET I =: 1
    20 LET S =: 0
    30 LET S =: S + 1/I
    40 LET I =: I + 1
    50 IF I <> 100000000 THEN 30
    60 PRINT "Almost infinity is"
    70 PRINT S
    80 END

And running it:

Lennarts-Computer% ghc --make Main.hs
[4 of 4] Compiling Main             ( Main.hs, Main.o )
Linking Main ...
Lennarts-Computer% ./Main
Almost infinity is
18.9979
CPU time:   1.57s

As you can see it's about the same time. In fact the assembly code for the loops look pretty much the same. Here's the Haskell one:

LBB1_1: ## _L4
        movsd   LCPI1_0, %xmm2
        movapd  %xmm1, %xmm3
        addsd   %xmm2, %xmm3
        ucomisd LCPI1_1, %xmm3
        divsd   %xmm1, %xmm2
        addsd   %xmm2, %xmm0
        movapd  %xmm3, %xmm1
        jne     LBB1_1  ## _L4

2009-02-06T11:12:00.002Z

Regression

They say that as you get older you regress back towards childhood. So I present you with today's Haskell program (the idea shamelessly stolen from JoshTriplett from #haskell on IRC):

import BASIC

main = runBASIC $ do

    10 LET X =: 1
    20 PRINT "Hello BASIC world!"
    30 LET X =: X + 1
    40 IF X <> 11 THEN 20
    50 END

Yes, it runs. (I'm sorry about the =: instead of =, but some things are just too wired into Haskell to change.)

2009-01-21T01:07:00.004Z

A performance update

I've continued playing with the LLVM. I discovered that when generating code for the normcdf and Black-Scholes functions I did not tell LLVM that the functions that were called (exp etc.) are actually pure functions. That meant that the LLVM didn't perform CSE properly.

So here are some updated numbers for computing an option prices for 10,000,000 options:

Pure Haskell: 8.7s
LLVM: 2.0s

Just as a reminder, the code for normcdf looks like this:

normcdf x = x %< 0 ?? (1 - w, w)
  where w = 1.0 - 1.0 / sqrt (2.0 * pi) * exp(-l*l / 2.0) * poly k
        k = 1.0 / (1.0 + 0.2316419 * l)
        l = abs x
        poly = horner coeff 
        coeff = [0.0,0.31938153,-0.356563782,1.781477937,-1.821255978,1.330274429]

A noteworthy thing is that exactly the same code can be used both for the pure Haskell and the LLVM code generation; it's just a matter of overloading.

Vectors

An very cool aspect of the LLVM is that it has vector instructions. On the x86 these translate into using the SSE extensions to the processor and can speed up computations by doing things in parallel.

Again, by using overloading, the exact same code can be used to compute over vectors of numbers as with individual numbers.

So what about performance? I used four element vectors of 32 bit floating point numbers and got these results:

Pure Haskell: 8.7s
LLVM, scalar: 2.0s
LLVM, vector: 1.1s
C, gcc -O3: 2.5s

Some caveats if you try out vectors in the LLVM.

Only on MacOS does the LLVM package give you fast primitive functions, because that's the only platform that seems to have this as a standard.
The vector version of floating point comparison has not been implemented in the LLVM yet.
Do not use two element vectors of type 32 bit floats. This will generate code that is wrong on the x86. I sent in a bug report about this, but was told that it is a feature and not a bug. (I kid you not.) To make the code right you have to manually insert EMMS instructions.
The GHC FFI is broken for all operations that allocate memory for a Storable, e.g., alloca, with, withArray etc. These operations do not take the alignment into account when allocating. This means that, e.g., a vector of four floats may end up on 8 byte alignment instead of 16. This generates a segfault.

On the whole, I'm pretty happy with the LLVM performance now; about 8 times faster than ghc on this example.

[Edit:] Added point about broken FFI.

2009-01-10T15:33:00.003Z

LLVM arithmetic

So we want to compute x²-5x+6 using the Haskell LLVM bindings. It would look some something like this.

    xsq <- mul x x
    x5  <- mul 5 x
    r1  <- sub xsq x5
    r   <- add r1 6

Not very readable, it would be nicer to write

    r   <- x^2 - 5*x + 6

But, e.g., the add instruction has the (simplified) type Value a -> Value a -> CodeGenFunction r (Value a), where CodeGenFunction is the monad for generating code for a function. (BTW, the r type variable is used to keep track of the return type of the function, used by the ret instruction.)

We can't change the return type of add, but we can change the argument type.

type TValue r a = CodeGenFunction r (Value a)
add' :: TValue r a -> TValue r a -> TValue r a
add' x y = do x' <- x; y' <- y; add x' y'

Now the type fits what the Num class wants. And we can make an instance declaration.

instance (Num a) => Num (TValue r a) where
    (+) = add'
    (-) = sub'
    (*) = mul'
    fromInteger = return . valueOf . fromInteger

We are getting close, the only little thing is that the x in our original LLVM code has type Value a rather than TValue r a, but return takes care of that. So:

    let x' = return x
    r <- x'^2 - 5*x' + 6

And a quick look at the generated LLVM code (for Double) shows us that all is well.

; x in %0
 %1 = mul double %0, %0
 %2 = mul double 5.000000e+00, %0
 %3 = sub double %1, %2
 %4 = add double %3, 6.000000e+00

All kinds of numeric instances and some other goodies are available in the LLVM.Util.Arithmetic module. Here is a complete Fibonacci (again) using this.

import Data.Int
import LLVM.Core
import LLVM.ExecutionEngine
import LLVM.Util.Arithmetic

mFib :: CodeGenModule (Function (Int32 -> IO Int32))
mFib = recursiveFunction $ \ rfib n ->
    n %< 2 ? (1, rfib (n-1) + rfib (n-2))

main :: IO ()
main = do
    let fib = unsafeGenerateFunction mFib
    print (fib 22)

The %< operator compares values returning a TValue r Bool. The c ? (t, e) is a conditional expression, like C's c ? t : e.

The type given to mFib is not the most general one, of course. The most general one can accept any numeric type for argument and result. Here it is:

mFib :: (Num a, Cmp a, IsConst a,
         Num b, Cmp b, IsConst b, FunctionRet b) =>
        CodeGenModule (Function (a -> IO b))
mFib = recursiveFunction $ \ rfib n ->
    n %< 2 ? (1, rfib (n-1) + rfib (n-2))

It's impossible to generate code for mFib directly; code can only be generated for monomorphic types. So a type signature is needed when generating code to force a monomorphic type.

main = do
    let fib :: Int32 -> Double
        fib = unsafeGenerateFunction mFib
        fib' :: Int16 -> Int64
        fib' = unsafeGenerateFunction mFib
    print (fib 22, fib' 22)

Another example

Let's try a more complex example. To pick something mathematical to have lots of formulae we'll do the Cumulative Distribution Function. For the precision of a Float it can be coded like this in Haskell (normal Haskell, no LLVM):

normcdf x = if x < 0 then 1 - w else w
  where w = 1 - 1 / sqrt (2 * pi) * exp(-l*l / 2) * poly k
        k = 1 / (1 + 0.2316419 * l)
        l = abs x
        poly = horner coeff 
        coeff = [0.0,0.31938153,-0.356563782,1.781477937,-1.821255978,1.330274429] 

horner coeff base = foldr1 multAdd coeff
  where multAdd x y = y*base + x

We cannot use this directly, it has type normcdf :: (Floating a, Ord a) => a -> a. The Ord context is a problem, because there are no instance of Ord for LLVM types. The comparison is the root of the problem, since it returns a Bool rather than a TValue r Bool.

It's possible to hide the Prelude and overload the comparisons, but you cannot overload the if construct. So a little rewriting is necessary regardless. Let's just bite the bullet and change the first line.

normcdf x = x %< 0 ? (1 - w, w)

And with that change, all we need to do is

mNormCDF = createFunction ExternalLinkage $ arithFunction $ normcdf

main :: IO ()
main = do
    writeFunction "CDF.bc" (mNormCDF :: CodeGenModule (Function (Float -> IO Float)))

So what happened here? Looking at normcdf it contains a things that the LLVM cannot handle, like lists. But all the list operations happen when the Haskell code runs and nothing of that is left in the LLVM code.

If you optimize and generate code for normcdf it looks like this (leaving out constants and half the code):

__fun1:
        subl    $28, %esp
        pxor    %xmm0, %xmm0
        ucomiss 32(%esp), %xmm0
        jbe     LBB1_2
        movss   32(%esp), %xmm0
        mulss   %xmm0, %xmm0
        divss   LCPI1_0, %xmm0
        movss   %xmm0, (%esp)
        call    _expf
        fstps   24(%esp)
        movss   32(%esp), %xmm0
        mulss   LCPI1_1, %xmm0
        movss   %xmm0, 8(%esp)
        movss   LCPI1_2, %xmm0
        movss   8(%esp), %xmm1
        addss   %xmm0, %xmm1
        movss   %xmm1, 8(%esp)
        movaps  %xmm0, %xmm1
        divss   8(%esp), %xmm1
        movaps  %xmm1, %xmm2
        mulss   LCPI1_3, %xmm2
        addss   LCPI1_4, %xmm2
        mulss   %xmm1, %xmm2
        addss   LCPI1_5, %xmm2
        mulss   %xmm1, %xmm2
        addss   LCPI1_6, %xmm2
        mulss   %xmm1, %xmm2
        addss   LCPI1_7, %xmm2
        mulss   %xmm1, %xmm2
        pxor    %xmm1, %xmm1
        addss   %xmm1, %xmm2
        movss   24(%esp), %xmm1
        mulss   LCPI1_8, %xmm1
        mulss   %xmm2, %xmm1
        addss   %xmm0, %xmm1
        subss   %xmm1, %xmm0
        movss   %xmm0, 20(%esp)
        flds    20(%esp)
        addl    $28, %esp
        ret
LBB1_2:
        ...

And that looks quite nice, in my opinion.

Black-Scholes

I work at a bank these days, so let's do the most famous formula in financial maths, the Black-Scholes formula for pricing options. Here's a Haskell rendition of it.

blackscholes iscall s x t r v = if iscall then call else put
  where call = s * normcdf d1 - x*exp (-r*t) * normcdf d2
        put  = x * exp (-r*t) * normcdf (-d2) - s * normcdf (-d1)
        d1 = (log(s/x) + (r+v*v/2)*t) / (v*sqrt t)
        d2 = d1 - v*sqrt t

Again, it uses an if, so it needs a little fix.

blackscholes iscall s x t r v  = iscall ? (call, put)

With that fix, code can be generated directly from blackscholes. The calls to normcdf are expanded inline, but it's easy to make some small changes to the code so that it actually does function calls.

Some figures

To test the speed of the generated code I ran blackscholes over a portfolio of 10,000,000 options and summed their value. The time excludes the time to set up the portfolio array, because that is done in Haskell. I also ran the code in pure Haskell on a list with 10,000,000 elements.

Unoptimized LLVM   17.5s
Optimized LLVM      8.2s
Pure Haskell        9.3s

The only surprise here is how well pure Haskell (with -O2) performs. This is a very good example for Haskell though, because the list gets fused away and everything is strict. A lot of the time is spent in log and exp in this code, so perhaps the similar results are not so surprising after all.

2009-01-07T16:14:00.005Z

LLVM

The LLVM, Low Level Virtual Machine, is a really cool compiler infrastructure project with many participants. The idea is that if you want to make a new high quality compiler you just have to generate LLVM code, and then there are lots of optimizations and code generators available to get fast code.

There are different ways to generate input to the LLVM tools. You can generate a text file with LLVM code and feed it to the tools, or you can use bindings for some programming language and programmatically build the LLVM code. The original bindings from the LLVM project is for C++, but they also provide C bindings. On top of the C bindings you can easily interface to other languages; for instance O'Caml and Haskell.

There are also diffent things you can do to LLVM code you have build programmatically. You can transform it, you can write to a file, you can run an interpreter on it, or execute it with a JIT compiler.

Haskell LLVM bindings

There is a Haskell binding to the LLVM. It has two layers. You can either work on the C API level and have ample opportunity to shoot your own limbs to pieces, or you can use the high level interface which is mostly safe.

Bryan O'Sullivan did all the hard work of taking the C header files and producing the corresponding Haskell FFI files. He also made a first stab at the high level interface, which I have since change a lot (for better or for worse).

An example

Let's do an example. We'll write the LLVM code for this function

  f x y z = (x + y) * z

In Haskell this function is polymorphic, but when generating machine code we have to settle for a type. Let's pick Int32. (The Haskell Int type cannot be used in talking to LLVM; it doesn't a well defined size.) Here is how it looks:

mAddMul :: CodeGenModule (Function (Int32 -> Int32 -> Int32 -> IO Int32))
mAddMul = 
  createFunction ExternalLinkage $ \ x y z -> do
    t <- add x y
    r <- mul t z
    ret r

For comparison, the LLVM code in text for for this would be:

define i32 @_fun1(i32, i32, i32) {
        %3 = add i32 %0, %1
        %4 = mul i32 %3, %2
        ret i32 %4
}

So what does the Haskell code say? The mAddMul definition is something in the CodeGenModule monad, and it generates a Function of type Int32 -> Int32 -> Int32 -> IO Int32. That last is the type of f above, except for that IO. Why the IO? The Haskell LLVM bindings forces all defined functions to return something in the IO monad, because there are no restriction on what can happen in the LLVM code; it might very well do IO. So to be on the safe side, there's always an IO on the type. If we know the function is harmless, we can use unsafePerformIO to get rid of it.

So the code does a createFunction which does what the name suggests. The ExternalLinkage argument says that this function will be available outside the module it's in, the obvious opposite being InternalLinkage. Using InternalLinkage is like saying static on the top level in C. In this examples it doesn't really matter which we pick.

The function has three arguments x y z. The last argument to createFunction should be a lambda expression with the right number of arguments, i.e., the number of arguments should agree with the type. We the use monadic syntax to generate an add, mul, and ret instruction.

The code looks like assembly code, which is the level that LLVM is at. It's a somewhat peculiar assembly code, because it's on SSA (Static Single Assignment) form. More about that later.

So what can we do with this function? Well, we can generate machine code for it and call it.

main = do
    addMul <- simpleFunction mAddMul
    a <- addMul 2 3 4
    print a

In this code addMul has type Int32 -> Int32 -> Int32 -> IO Int32, so it has to be called in the IO monad. Since this is a pure function, we can make the type pure, i.e., Int32 -> Int32 -> Int32 -> Int32.

main = do
    addMul <- simpleFunction mAddMul
    let addMul' = unsafePurify addMul
    print (addMul' 2 3 4)

The unsafePurify functions is simply an extension of unsafePerformIO that drops the IO on the result of a function.

So that was pretty easy. To make a function, just specify the LLVM code using the LLVM DSEL that the Haskell bindings provides.

Fibonacci

No FP example is complete without the Fibonacci function, so here it is.

mFib :: CodeGenModule (Function (Word32 -> IO Word32))
mFib = do
    fib <- newFunction ExternalLinkage
    defineFunction fib $ \ arg -> do
        -- Create the two basic blocks.
        recurse <- newBasicBlock
        exit <- newBasicBlock

        -- Test if arg > 2
        test <- icmp IntUGT arg (2::Word32)
        condBr test recurse exit

        -- Just return 1 if not > 2
        defineBasicBlock exit
        ret (1::Word32)

        -- Recurse if > 2, using the cumbersome plus to add the results.
        defineBasicBlock recurse
        x1 <- sub arg (1::Word32)
        fibx1 <- call fib x1
        x2 <- sub arg (2::Word32)
        fibx2 <- call fib x2
        r <- add fibx1 fibx2
        ret r
    return fib

Instead of using createFunction to create the function we're using newFunction and defineFunction. The former is a shorthand for the latter two together. But splitting making the function and actually defining it means that we can refer to the function before it's been defined. We need this since fib is recursive.

Every instruction in the LLVM code belongs to a basic block. A basic block is a sequence of non-jump instructions (call is allowed in the LLVM) ending with some kind of jump. It is always entered at the top only. The top of each basic block can be thought of as a label that you can jump to, and those are the only places that you can jump to.

The code for fib starts with a test if the argument is Unsigned Greater Than 2. The condBr instruction branches to recurse if test is true otherwise to exit. To be able to refer to the two branch labels (i.e., basic blocks) before they are defined we create them with newBasicBlock and then later define them with defineBasicBlock. The defineBasicBlock simply starts a new basic block that runs to the next basic block start, or to the end of the function. The type system does not check that the basic block ends with a branch (I can't figure out how to do that without making the rest of the code more cumbersome).

In the false branch we simply return 1, and in the true branch we make the two usual recursive calls, add the results, and return the sum.

As you can see a few type annotations are necessary on constants. In my opinion they are quite annoying, because if you write anything different from ::Word32 in those annotations there will be a type error. This means that in principle the compiler has all the information, it's just too "stupid" to use it.

The performance you get from this Fibonacci function is decent, but in fact worse than GHC with -O2 gives. Even with full optimization turned on for the LLVM code it's still not as fast as GHC for this function.

[Edit: Added assembly] Here is the assembly code for Fibonacci. Note how there is only one recursive call. The other call has been transformed into a loop.

_fib:
 pushl %edi
 pushl %esi
 subl $4, %esp
 movl 16(%esp), %esi
 cmpl $2, %esi
 jbe LBB1_4
LBB1_1:
 movl $1, %edi
 .align 4,0x90
LBB1_2:
 leal -1(%esi), %eax
 movl %eax, (%esp)
 call _fib
 addl %edi, %eax
 addl $4294967294, %esi
 cmpl $2, %esi
 movl %eax, %edi
 ja LBB1_2
LBB1_3:
 addl $4, %esp
 popl %esi
 popl %edi
 ret
LBB1_4:
 movl $1, %eax
 jmp LBB1_3

Hello, World!

The code for printing "Hello, World!":

import Data.Word
import LLVM.Core
import LLVM.ExecutionEngine

bldGreet :: CodeGenModule (Function (IO ()))
bldGreet = do
    puts <- newNamedFunction ExternalLinkage "puts" :: TFunction (Ptr Word8 -> IO Word32)
    greetz <- createStringNul "Hello, World!"
    func <- createFunction ExternalLinkage $ do
      tmp <- getElementPtr greetz (0::Word32, (0::Word32, ()))
      call puts tmp -- Throw away return value.
      ret ()
    return func

main :: IO ()
main = do
    greet <- simpleFunction bldGreet
    greet

To get access to the C function puts we simply declare it and rely on the linker to link it in. The greetz variable has type pointer to array of characters. So to get a pointer to the first character we have to use the rather complicated getElementPtr instruction. See FAQ about it.

Phi instructions

Let's do the following simple C function

int f(int x)
{
  if (x < 0) x = -x;
  return (x+1);
}

Let's try to write some corresponding LLVM code:

  createFunction ExternalLinkage $ \ x -> do
    xneg <- newBasicBlock
    xpos <- newBasicBlock
    t <- icmp IntSLT x (0::Int32)
    condBr t xneg xpos

    defineBasicBlock xneg
    x' <- sub (0::Int32) x
    br xpos

    defineBasicBlock xpos
    r1 <- add ??? (1::Int32)
    ret r1

But what should we put at ???? When jumping from the condBr the value is in x, but when jumping from the negation block the value is in x'. And this is how SSA works. Every instruction puts the value in a new "register", so this situation is unavoidable. This is why SSA (and thus LLVM) form has phi instructions. This is a pseudo-instruction to tell the code generator what registers should be merged at the entry of a basic block. So the real code looks like this:

mAbs1 :: CodeGenModule (Function (Int32 -> IO Int32))
mAbs1 = 
  createFunction ExternalLinkage $ \ x -> do
    top <- getCurrentBasicBlock
    xneg <- newBasicBlock
    xpos <- newBasicBlock
    t <- icmp IntSLT x (0::Int32)
    condBr t xneg xpos

    defineBasicBlock xneg
    x' <- sub (0::Int32) x
    br xpos

    defineBasicBlock xpos
    r <- phi [(x, top), (x', xneg)]
    r1 <- add r (1::Int32)
    ret r1

The phi instruction takes a list of registers to merge, and paired up with each register is the basic block that the jump comes from. Since the first basic block in a function is created implicitely we have to get it with getCurrentBasicBlock which returns the current basic block.

If, like me, you have a perverse interest in the machine code that gets generated here is the optimized code for that function on for x86:

__fun1:
        movl    4(%esp), %eax
        movl    %eax, %ecx
        sarl    $31, %ecx
        addl    %ecx, %eax
        xorl    %ecx, %eax
        incl    %eax
        ret

Note how the conditional jump has cleverly been replaced by some non-jumping instructions. I think this code is as good as it gets.

Loops and arrays

Let's do a some simple array code, the dot product of two vectors. The function takes a length and pointers to two vectors. It sums the elementwise product of the vectors. Here's the C code:

double
dotProd(unsigned int len, double *aPtr, double *bPtr)
{
    unsigned int i;
    double s;

    s = 0;
    for (i = 0; i != len; i++)
        s += aPtr[i] * bPtr[i];
    return s;
}

The corresponding LLVM code is much more complicated and has some new twists.

import Data.Word
import Foreign.Marshal.Array
import LLVM.Core
import LLVM.ExecutionEngine

mDotProd :: CodeGenModule (Function (Word32 -> Ptr Double -> Ptr Double -> IO Double))
mDotProd =
  createFunction ExternalLinkage $ \ size aPtr bPtr -> do
    top <- getCurrentBasicBlock
    loop <- newBasicBlock
    body <- newBasicBlock
    exit <- newBasicBlock

    -- Enter loop, must use a br since control flow joins at the loop bb.
    br loop

    -- The loop control.
    defineBasicBlock loop
    i <- phi [(valueOf (0 :: Word32), top)]  -- i starts as 0, when entered from top bb
    s <- phi [(valueOf 0, top)]  -- s starts as 0, when entered from top bb
    t <- icmp IntNE i size       -- check for loop termination
    condBr t body exit

    -- Define the loop body
    defineBasicBlock body

    ap <- getElementPtr aPtr (i, ()) -- index into aPtr
    bp <- getElementPtr bPtr (i, ()) -- index into bPtr
    a <- load ap                 -- load element from a vector
    b <- load bp                 -- load element from b vector
    ab <- mul a b                -- multiply them
    s' <- add s ab               -- accumulate sum

    i' <- add i (valueOf (1 :: Word32)) -- Increment loop index

    addPhiInputs i [(i', body)]  -- Control flow reaches loop bb from body bb
    addPhiInputs s [(s', body)]
    br loop                      -- And loop

    defineBasicBlock exit
    ret (s :: Value Double)      -- Return sum

main = do
    ioDotProd <- simpleFunction mDotProd
    let dotProd a b =
         unsafePurify $
         withArrayLen a $ \ aLen aPtr ->
         withArrayLen b $ \ bLen bPtr ->
         ioDotProd (fromIntegral (aLen `min` bLen)) aPtr bPtr

    let a = [1,2,3]
        b = [4,5,6]
    print $ dotProd a b
    print $ sum $ zipWith (*) a b

First we have to set up the looping machinery. There a four basic blocks involved: the implicit basic block that is created at the start of every function, top; the top of the loop, loop; the body of the loop, body; and finally the block with the return from the function, exit.

There are two "registers", the loop index i and the running sum s that arrive from two different basic blocks at the top of the loop. When entering the loop from the first time they should be 0. That's what the phi instruction specifies. The valueOf function simply turns a constant into an LLVM value. It's worth noting that the initial values for the two variables are constant rather than registers. The control flow also reached the basic block loop from the end of body, but we don't have the names of those registers in scope yet, so we can't put them in the phi instruction. Instead, we have to use addPhiInputs to add more phi inputs later (when the registers are in scope).

The most mysterious instruction in the LLVM is getElementPtr. It simply does address arithmetic, so it really does something quite simple. But it can perform several levels of address arithmetic when addressing through multilevel arrays and structs. In can take several indicies, but since here we simply want to add the index variable to a pointer the usage is pretty simple. Doing getElementPtr aPtr (i, ()) corresponds to aPtr + i in C.

To test this function we need pointers to two vectors. The FFI function withArrayLen temporarily allocates the vector and fills it with elements from the list.

The essential part of the function looks like this in optimized x86 code:

        pxor    %xmm0, %xmm0
        xorl    %esi, %esi
        .align  4,0x90
LBB1_2:
        movsd   (%edx,%esi,8), %xmm1
        mulsd   (%ecx,%esi,8), %xmm1
        incl    %esi
        cmpl    %eax, %esi
        addsd   %xmm1, %xmm0
        jne     LBB1_2

Which is pretty good. Improving this would have to use SSD vector instructions. This is possible using the LLVM vector type, but I'll leave that for now.

Abstraction

The loop structure in dotProd is pretty common, so we would like to abstract it out for reuse. The creation of basic blocks and phi instructions is rather fiddly so it would be nice to do this once and not worry about it again.

What are the parts of the loop? Well, let's just do a simple "for" loop that loops from a lower index (inclusive) to an upper index (exclusive) and executes the loop body for each iteration. So there should be three arguments to the loop function: lower bound, upper bound and loop body. What is the loop body? Since the LLVM is using SSA the loop body can't really update the loop state variables. Instead it's like a pure functional language where you have to express it as a state transformation. So the loop body will take the old state and return a new state. It's also useful to pass the loop index to the loop body. Now when we've introduced the notion of a loop state we also need to have an initial value for the loop state as an argument to the loop function.

Let's start out easy and let the state to be updated in the loop be a single value. In dotProd it's simply the running sum (s).

forLoop low high start incr = do
    top <- getCurrentBasicBlock
    loop <- newBasicBlock
    body <- newBasicBlock
    exit <- newBasicBlock

    br loop

    defineBasicBlock loop
    i <- phi [(low, top)]
    state <- phi [(start, top)]
    t <- icmp IntNE i high
    condBr t body exit

    defineBasicBlock body

    state' <- incr i state
    i' <- add i (valueOf 1)

    body' <- getCurrentBasicBlock
    addPhiInputs i [(i', body')]
    addPhiInputs state [(state', body')]
    br loop
    defineBasicBlock exit

    return state

The low and high arguments are simply the loop bounds, start is the start value for the loop state variable, and finally incr is invoked in the loop body to get the new value for the state variable. Note that the incr can contain new basic blocks so there's no guarantee we're in the same basic block after incr has been called. That's why there is a call to getCurrentBasicBlock before adding to the phi instructions.

So the original loop in dotProd can now be written

    s <- forLoop 0 size 0 $ \ i s -> do
      ap <- getElementPtr aPtr (i, ()) -- index into aPtr
      bp <- getElementPtr bPtr (i, ()) -- index into bPtr
      a <- load ap                 -- load element from a vector
      b <- load bp                 -- load element from b vector
      ab <- mul a b                -- multiply them
      s' <- add s ab               -- accumulate sum
      return s'

So that wasn't too bad. But what if the loop needs multiple state variables? Or none? The tricky bit is handling the phi instructions since the number of instructions needed depends on how many state variables we have. So let's creat a class for types that can be state variables. This way we can use tuples for multiple state variables. The class needs two methods, the generalization of phi and the generalization of addPhiInputs.

class Phi a where
    phis :: BasicBlock -> a -> CodeGenFunction r a
    addPhis :: BasicBlock -> a -> a -> CodeGenFunction r ()

A simple instance is when we have no state variables.

instance Phi () where
    phis _ _ = return ()
    addPhis _ _ _ = return ()

We also need to handle the case with a single state variable. All LLVM values are encapsulated in the Value type, so this is the one we create an instance for.

instance (IsFirstClass a) => Phi (Value a) where
    phis bb a = do
        a' <- phi [(a, bb)]
        return a'
    addPhis bb a a' = do
        addPhiInputs a [(a', bb)]

Finally, here's the instance for pair. Other tuples can be done in the same way (or we could just use nested pairs).

instance (Phi a, Phi b) => Phi (a, b) where
    phis bb (a, b) = do
        a' <- phis bb a
        b' <- phis bb b
        return (a', b')
    addPhis bb (a, b) (a', b') = do
        addPhis bb a a'
        addPhis bb b b'

Using this new class the looping function becomes

forLoop :: forall i a r . (Phi a, Num i, IsConst i, IsInteger i, IsFirstClass i) =>
           Value i -> Value i -> a -> (Value i -> a -> CodeGenFunction r a) -> CodeGenFunction r a
forLoop low high start incr = do
    top <- getCurrentBasicBlock
    loop <- newBasicBlock
    body <- newBasicBlock
    exit <- newBasicBlock

    br loop

    defineBasicBlock loop
    i <- phi [(low, top)]
    vars <- phis top start
    t <- icmp IntNE i high
    condBr t body exit

    defineBasicBlock body

    vars' <- incr i vars
    i' <- add i (valueOf 1 :: Value i)

    body' <- getCurrentBasicBlock
    addPhis body' vars vars'
    addPhiInputs i [(i', body')]
    br loop
    defineBasicBlock exit

    return vars

File operations

The Haskell bindings provide two convenient functions - writeBitcodeToFile and readBitcodeFromFile - for writing and reading modules in the LLVM binary format.

A simple example:

import Data.Int
import LLVM.Core

mIncr :: CodeGenModule (Function (Int32 -> IO Int32))
mIncr = 
  createNamedFunction ExternalLinkage "incr" $ \ x -> do
    r <- add x (1 :: Int32)
    ret r

main = do
    m <- newModule
    defineModule m mIncr
    writeBitcodeToFile "incr.bc" m

Running this will produce the file incr.bc which can be processed with the usual LLVM tools. E.g.

$ llvm-dis < incr.bc  # to look at the LLVM code
$ opt -std-compile-opts incr.bc -f -o incrO.bc # run optimizer
$ llvm-dis < incrO.bc  # to look at the optimized LLVM code
$ llc incrO.bc # generate assembly code
$ cat incrO.s  # look at assembly code

Reading a module file is equally easy, but what can you do with a module you have read? It could contain anything. To extract things from a module there is a function getModuleValues which returns a list of name-value pairs of all externally visible functions and global variables. The values all have type ModuleValue. To convert a ModuleValue to a regular Value you have to use castModuleValue. This is a safe conversion function that makes a dynamic type test to make sure the types match (think of ModuleValue as Dynamic and castModuleValue as fromDynamic).

Here's an example:

import Data.Int
import LLVM.Core
import LLVM.ExecutionEngine

main = do
    m <- readBitcodeFromFile "incr.bc"
    ee <- createModuleProviderForExistingModule m >>= createExecutionEngine
    funcs <- getModuleValues m
    let ioincr :: Function (Int32 -> IO Int32)
        Just ioincr = lookup "incr" funcs >>= castModuleValue
        incr = unsafePurify $ generateFunction ee ioincr

    print (incr 41)

This post is getting rather long, so I'll let this be the last example for today.

2008-12-10T15:27:00.002Z

The OCaml code again

I'm posting a slight variation of the OCaml code that I think better shows what I like about the ML version.

module MkPerson(O: sig 
                     type xString
                     type xDouble
                     val opConcat : xString -> xString -> xString
                     val opShow : xDouble -> xString
                   end) =
struct
  open O
  type person = Person of (xString * xString * xDouble)
  let display (Person (firstName, lastName, height)) = 
    opConcat firstName (opConcat lastName (opShow height))
end

module BasicPerson = MkPerson(struct
                                type xString = string
                                type xDouble = float
                                let opConcat = (^)
                                let opShow = string_of_float
                              end)

open BasicPerson

let _ = 
  let p = Person ("Stefan", "Wehr", 184.0)
  in display p

Note how the open O opens the argument to the MkPerson functor and all the values and types from the argument is in scope in the rest of the module. There's no need to change lots of code in MkPerson.

Similarely, the open BasicPerson makes the operations from that module avaiable, so that Person and display can be used in a simple way.

2008-12-10T12:34:00.002Z

Abstracting on, suggested solutions

I guess I should be more constructive than just whining about how Haskell doesn't always do what I want. I do have some suggestions on how to fix things.

Explicit type applications

Let's look at a simple example again:

f :: forall a . a -> a
f = \ x -> x

b :: Bool
b = f True

The way I like to think of this (and what happens in ghc) is that this is shorthand for something more explicit, namely the F_w version of the same thing. In F_w all type abstraction and type application are explicit. Let's look at the explicit version (which is no longer Haskell).

f :: forall (a::*) . a -> a
f = /\ (a::*) -> \ (x::a) -> x

b :: Bool
b = f @Bool True

I'm using /\ for type abstraction and expr @type for type application. Furthermore each binder is annotated with its type. This is what ghc translates the code to internally, this process involves figuring out what all the type abstractions and applications should be.

Not something a little more complicated (from my previous post)

class C a b where
    x :: a
    y :: b

f :: (C a b) => a -> [a]
f z = [x, x, z]

The type of x is

x :: forall a b . (C a b) => a

So whenever x occurs two type applications have to be inserted (there's also a dictionary to insert, but I'll ignore that). The decorated term for f (ignoring the context)

f :: forall a b . (C a b) => a -> [a]
f = /\ (a::*) (b::*) -> \ (z::a) -> [ x @a @b1, x @a @b2, z]

The reason for the ambiguity in type checking is that the type check cannot figure out that the b is in any way connected to b1 and b2. Because it isn't. And there's currently no way we can connect them.

So I suggest that it should be possible to use explicit type application in Haskell when you want to. The code would look like this

f :: forall a b . (C a b) => a -> [a]
f z = [ x @a @b, x @a @b, z]

The order of the variables in the forall determines the order in which the type abstractions come, and thus determines where to put the type applications.

Something like `abstype`

Back to my original problem with abstraction. What about if this was allowed:

class Ops t where
    data XString t :: *
    (+++) :: XString t -> XString t -> XString t

instance Ops Basic where
    type XString Basic = String
    (+++) = (++)

So the class declaration says I'm going to use data types (which was my final try and which works very nicely). But in the instance I provide a type synonym instead. This would be like using a newtype in the instance, but without having to use the newtype constructor everywhere. The fact that it's not a real data type is only visible inside the instance declaration. The compiler could in fact make a newtype and insert all the coercions. This is, of course, just a variation of the abstype suggestion by Wehr and Chakravarty.

2008-12-10T10:55:00.002Z

The abstraction continues

I got several comments to my lament about my attempts at abstraction in my previous blog post. Two of the comments involve adding an extra argument to display. I dont regard this as an acceptable solution; the changes to the code should not be that intrusive. Adding an argument to a function is a change that ripples through the code to many places and not just the implementation of display.

Reiner Pope succeeded where I failed. He split up the operations in Ops into two classes and presto, it works.

data Person t = Person {
    firstName :: XString t,
    lastName :: XString t,
    height :: XDouble t
    }

class (Show s, IsString s) => IsXString s where
    (+++) :: s -> s -> s
class (Num d, IsXString s) => IsXDouble s d where
    xshow :: d -> s

class (IsXDouble (XString t) (XDouble t)) => Ops t where
    type XString t :: *
    type XDouble t :: *
instance IsXString String where
    (+++) = (++)
instance IsXDouble String Double where
    xshow = show

data Basic = Basic

instance Ops Basic where
    type XString Basic = String
    type XDouble Basic = Double

display :: Ops t => Person t -> XString t
display p = firstName p +++ " " +++ lastName p +++ " " +++ xshow (height p + 1)

That's neat, but a little fiddly if there are many types involved.

Another problem

Armed with this solution I write another function.

incSpace :: (Ops t) => XDouble t -> XString t
incSpace x = xshow x +++ " "

It typechecks fine. But as far as I can figure out there is no way to use this function. Let's see what ghci says:

> :t incSpace (1 :: XDouble Basic) :: XString Basic

:1:0:
    Couldn't match expected type `[Char]'
           against inferred type `XString t'
    In the expression: incSpace (1 :: XDouble Basic) :: XString Basic

:1:10:
    Couldn't match expected type `XDouble t'
           against inferred type `Double'
    In the first argument of `incSpace', namely `(1 :: XDouble Basic)'
    In the expression: incSpace (1 :: XDouble Basic) :: XString Basic

Despite my best efforts at providing types it doesn't work. The reason being that saying, e.g., (1 :: XDouble Basic) is the same as saying (1 :: Double). And that doesn't match XDouble t. At least not to the typecheckers knowledge.

In the example of display things work because the parameter t occurs in Person t which is a real type and not a type family. If a type variable only occurs in type family types you are out of luck. There's no way to convey the information what that type variable should be (as far as i know). You can "solve" the problem by adding t as an argument to incSpace, but again, I don't see that as a solution.

In the paper ML Modules and Haskell Type Classes: A Constructive Comparison Wehr and Chakravarty introduce a notion of abstract associated types. That might solve this problem. I really want XDouble and XString to appear as abstract types (or associated data types) outside of the instance declaration. Only inside the instance declaration where I provide implementations for the operations do I really care what the type is.

A reflection on type signatures

If I write

f x = x

Haskell can deduce that the type is f :: a -> a.

If I instead write

f :: Int -> Int
f x = x

Haskell happily uses this type. The type checker does not complain as to say "Sorry dude, but you're wrong, the type is more general than what you wrote.". I think that's nice and polite.

Now a different example.

class C a b where
    x :: a
    y :: b

f z = [x, x, z]

What does ghc have to say about the type of f?

f :: (C a b, C a b1) => a -> [a]

OK, that's reasonable; the two occurences of x could have different contexts. But I don't want them to. Let's add a type signature.

f :: (C a b) => a -> [a]
f z = [x, x, z]

What does ghc have to say?

Blog2.hs:9:7:
    Could not deduce (C a b) from the context (C a b2)
      arising from a use of `x' at Blog2.hs:9:7
    Possible fix:
      add (C a b) to the context of the type signature for `f'
    In the expression: x
    In the expression: [x, x, z]
    In the definition of `f': f z = [x, x, z]

Blog2.hs:9:10:
    Could not deduce (C a b1) from the context (C a b2)
      arising from a use of `x' at Blog2.hs:9:10
    Possible fix:
      add (C a b1) to the context of the type signature for `f'
    In the expression: x
    In the expression: [x, x, z]
    In the definition of `f': f z = [x, x, z]

Which is ghc's way of say "Dude, I see your context, but I'm not going to use it because I'm more clever than you and can figure out a better type." Rude, is what I say.

I gave a context, but there is nothing to link the b in my context to what ghc internally figures out that the type of the two occuerences of x should. I wish I could tell the type checker, "This is the only context you'll ever going to have, use it if you can." Alas, this is not how things work.

A little ML

Stefan Wehr provided the ML version of the code that I only aluded to

module MkPerson(O: sig 
                     type xString
                     type xDouble
                     val opConcat : xString -> xString -> xString
                     val opShow : xDouble -> xString
                   end) =
struct
  type person = Person of (O.xString * O.xString * O.xDouble)
  let display (Person (firstName, lastName, height)) = 
    O.opConcat firstName (O.opConcat lastName (O.opShow height))
end

module BasicPerson = MkPerson(struct
                                type xString = string
                                type xDouble = float
                                let opConcat = (^)
                                let opShow = string_of_float
                              end)

let _ = 
  let p = BasicPerson.Person ("Stefan", "Wehr", 184.0)
  in BasicPerson.display p

In this case, I think this is the natural way of expressing the abstraction I want. Of course, this ML code has some shortcomings too. Since string literals in ML are not overloaded the cannot be used neatly in the display function like I could in the Haskell version, but that's a minor point.

2008-12-10T00:10:00.002Z

A somewhat failed adventure in Haskell abstraction

I usually blog about weird and wonderful things you can do in Haskell. Today I'm going to talk about something very plain and not wonderful at all.

If you want to try out the code below, use these Haskell extensions:

{-# LANGUAGE TypeFamilies, MultiParamTypeClasses, OverloadedStrings,
   FlexibleInstances, TypeSynonymInstances, ScopedTypeVariables,
   FunctionalDependencies, RecordWildCards, FlexibleContexts,
   GeneralizedNewtypeDeriving #-}

The simple problem

We want to define a type for a person which has a few fields and operations. Like this

module Person(Person(..), display) where

data Person = Person {
   firstName :: String,
   lastName  :: String,
   height    :: Double
   }

display :: Person -> String
display p = firstName p ++ " " ++ lastName p ++ " " ++ show (height p + 1)

Very simple. To use it we can just import the module and the write something like

 print $ display $ Person { firstName = "Manuel", lastName = "Peyton Jones", height = 255 }

But being efficiancy conscious I'm not happy with using String and Double. I'd like to experiment with different types for these. Maybe I should use ByteString and Int instead?

Simple enough, let's abstract out the types and operations into a different module.

module Ops(XString, XDouble, (+++), xshow) where
import Data.String
newtype XString = XString String deriving (Eq, Show, IsString)
newtype XDouble = XDouble Double deriving (Eq, Show, Num)

(+++) :: XString -> XString -> XString
XString x +++ XString y = XString (x ++ y)

xshow :: XDouble -> XString
xshow (XDouble x) = XString (show x)

module Person(Person(..), display) where
import Ops

data Person = Person {
   firstName :: XString,
   lastName  :: XString,
   height    :: XDouble
   }

display :: Person -> XString
display p = firstName p +++ " " +++ lastName p +++ " " +++ show (height p + 1)

There, problems solved. By changing the import in the Person module you can try out different types for XString and XDouble.

No, this is not problem solved. To try out different implementations I need to edit the Person module. That's not abstraction, that's obstruction. It should be possible to write the code for the Person module once and for all once you decided to abstract, and then never change it again.

I also didn't really want to necessarily have newtype in my module. Maybe I'd want this:

module Ops(XString, XDouble, (+++), xshow) where
type XString = String
type XDouble = Double

(+++) :: XString -> XString -> XString
(+++) = (++)

xshow :: XDouble -> XString
xshow = show

You can define Ops that way, but then the implementation of Ops may leak into the Person module. What you really want is to type check Person against the signature of Ops, like

interface Ops where
type XString
type XDouble
(+++) :: XString -> XString -> XString
xshow :: XDouble -> XString

And later supply the actual implementation. Alas, Haskell doesn't allow this.

In ML (SML or O'Caml) this would be solved by using a functor. The Person module would be a functor that takes the Ops module as an argument and yields a new module. And then you can just plug and play with different Ops implementations. This is where ML shines and Haskell sucks.

Type classes

But the defenders of Haskell superiority say, Haskell has type classes, that's the way to abstract! So let's make Ops into a type class. Let's do old style with multiple parameters first. Since Ops defines two types it will correspond to having two type parameters to the class.

class (IsString xstring, Num xdouble) => Ops xstring xdouble where
   (+++) :: xstring -> xstring -> xstring
   xshow :: xdouble -> xshow

Ok, so how do we have to rewrite the Person module?

data Person xstring xdouble = Person {
   firstName :: xstring,
   lastName  :: xstring,
   height    :: xdouble
   }

display :: (Ops xstring xdouble) => Person xstring xdouble -> xstring
display p = firstName p +++ " " +++ lastName p +++ " " +++ xshow (height p + 1)

An implementation is provided by an instance declaration:

instance Ops String Double where
   (+++) = (++)
   xshow = show

We see the major flaw in this approch at once. The Person data type now has two parameters. This might be bearable, but imagine a more complicated example where Ops contains 15 types. And every time you add a field with a new type to Person you have to update every single place in the program that mentions the Person type. That's not abstraction.

But in fact, it's even worse than that. The definition of display might look plausible, but it's full of ambiguities. Compiling it gives lots of errors of this kind:

   Could not deduce (Ops xstring xdouble)
     from the context (Ops xstring xdouble4)

Well, we can remove the type signature and let GHC figure it out. The we get this

display :: (Ops xstring xdouble,
           Ops xstring xdouble3,
           Ops xstring xdouble2,
           Ops xstring xdouble1,
           Ops xstring xdouble4) =>
          Person xstring xdouble4 -> xstring

And this function can, of course, never be used because most of the type variables do not occur outside the context so they will never be determined. I don't even know how to put explicit types in the function to make it work.

Well, it's common knowledge that multi-parameter type classes without functional dependencies is asking for trouble. So can we add some functional dependencies? Sure, if we use

class (IsString xstring, Num xdouble) => Ops xstring xdouble | xstring -> xdouble where

then things work beautifully. Until we decide that another instance that would be interesting to try is

instance Ops String Int

which is not valid with the FD present.

So we can't have functional dependencies if we want to have flexibilty with the instances. So what is it that goes wrong without the FDs? It's that all the uses (+++) and xshow are not tied together, they could potentially have different types. Let's try and be sneaky and tie them together:

display :: (Ops xstring xdouble) => Person xstring xdouble -> xstring
display p =
   let (++++) = (+++); xxshow = xshow
   in  firstName p ++++ " " ++++ lastName p ++++ " " ++++ xxshow (height p + 1)

This only generates one error message, because there's still nothing that says the the two operations come from the same instance. We need to make the tie even closer.

class (IsString xstring, Num xdouble) => Ops xstring xdouble where
   ops :: (xstring -> xstring -> xstring, xdouble -> xstring)
instance Ops String Double where
   ops = ((++), show)

display :: (Ops xstring xdouble) => Person xstring xdouble -> xstring
display p =
   let ((+++), xshow) = ops
   in  firstName p +++ " " +++ lastName p +++ " " +++ xshow (height p + 1)

This actually works! We can make it neater looking.

class (IsString xstring, Num xdouble) => Ops xstring xdouble where
   ops :: DOps xstring xdouble

data DOps xstring xdouble = DOps {
   (+++) :: xstring -> xstring -> xstring,
   xshow :: xdouble -> xstring
   }

instance Ops String Double where
   ops = DOps (++) show

display :: (Ops xstring xdouble) => Person xstring xdouble -> xstring
display p =
   let DOps{..} = ops
   in  firstName p +++ " " +++ lastName p +++ " " +++ xshow (height p + 1)

We have basically packaged up the dictionary and unpack it ourselves to get access to the operations. It's not pleasent, but it works.

But as I already said, the multiparameter type class version isn't really a good solution to the problem even if it works; it introduces too many parameters to the Person record.

Associated types

The new and shiny way of doing type classes is to use associated types instead of FDs. So let's give that a try. So what should the associated types be in the class. The associated type is supposed to be the one that can be computed from the main one. But we have two types that are on equal footing, so there is no main one. We can remedy that by introducing an artificial third type that is the main one, it can then determine the other two.

class (IsString (XString t), Num (XDouble t)) => Ops t where
   type XString t :: *
   type XDouble t :: *
   (+++) :: XString t -> XString t -> XString t
   xshow :: XDouble t -> XString t

data Person t = Person {
   firstName :: XString t,
   lastName  :: XString t,
   height    :: XDouble t
   }

That looks pretty neat. Note how the Person record has one parameter and no matter how many new associated type we add it will still only have one parameter. One parameter is reasonable, the Person record is after all parameterized over what kind of Ops we are providing.

Let's do an instance. It will need the extra type that is somehow the name of the instance.

data Basic = Basic

instance Ops Basic where
   type XString Basic = String
   type XDouble Basic = Double
   (+++) = (++)
   xshow = show

Now what about the display function? Alas, now it breaks down again. The display function is full of type errors again. And the reason is similar to the multiparameter version; there's nothing that ties the operations together.

We can play the same trick as with DOps above, but for some reason it doesn't work this time. The type comes out as

display :: (XString t ~ XString a,
           XDouble t ~ XDouble a,
           Ops a,
           Num (XDouble t)) =>
          Person t -> XString a

I have no clue why. I find associated types very hard to get a grip on.

OK, multi-parameter type classes made things work, but had too many type parameters. And associated types is the other way around. You can try combining them, but it didn't get me anywhere closer.

Associated data types

OK, I won't admit defeat yet. There's still associated data types. They are easier to deal with than associated types, because the type function is guaranteed to be injective.

class (IsString (XString t), Num (XDouble t)) => Ops t where
   data XString t :: *
   data XDouble t :: *
   (+++) :: XString t -> XString t -> XString t
   xshow :: XDouble t -> XString t

data Basic = Basic

instance Ops Basic where
   newtype XString Basic = XSB String deriving (Eq, Ord, Show)
   newtype XDouble Basic = XDB Double deriving (Eq, Ord, Show)
   XSB x +++ XSB y = XSB (x ++ y)
   xshow (XDB x) = XSB (show x)
instance Num (XDouble Basic) where
   XDB x + XDB y = XDB (x+y)
   fromInteger = XDB . fromInteger
instance IsString (XString Basic) where
   fromString = XSB

At last, this actually works! But it's at a price. We can no longer use the types we want in the instance declaration, instead we are forced to invent new types. Using this approach the original multi-parameter version could have been made to work as well.

Normally the GeneralizedNewtypeDeriving language extension makes it relatively painless to introduce a newtype that has all the instances of the underlying type. But due to a bug in ghc you can't use this extension for associated newtypes. So we have to make manual instance declarations which makes this approach very tedious.

Conclusion

I have found no way of doing what I want. My request is very simple, I want to be able to abstract over the actual implementation of a module, where the module contains types, values, and instances.

Haskell normally excels in abstraction, but here I have found no natural way of doing what I want. Perhaps I'm just not clever enough to figure out how, but that is a failure of Haskell too. It should not take any cleverness to do something as simple as this. In ML this is the most natural thing in the world to do.

Associated types are not a replacement for a proper module system. They let you do some things, but others just don't work.

I'd be happy to see anyone doing this in Haskell in a simple way.

2008-07-03T12:38:00.018+01:00

Lost and found

If I write 10^8 in Haskell, how many multiplications will be used to compute the power? A stupid question? Well, for this example, but if I was computing x^8 and x has 100000 digits then I'd care. So how can I find out? I can look at the definition of the exponentiation operator. Here it is, from the Haskell report and GHC 6.8:

(^)             :: (Num a, Integral b) => a -> b -> a
_ ^ 0           =  1
x ^ n | n > 0   =  f x (n-1) x
   where f _ 0 y = y
         f a d y = g a d
           where
             g b i | even i  = g (b*b) (i `quot` 2)
                   | otherwise = f b (i-1) (b*y)
_ ^ _           = error "Prelude.^: negative exponent"

It's a bit involved, but decipherable. Another way would be to insert some kind of debug trace message in the multiplication.

Traced values

I'd like to show a different way. Here's a ghci session:

Prelude> :m +Debug.Traced
Prelude Debug.Traced> let x = 10 :: Traced AsValue Integer
Prelude Debug.Traced> let y = x ^ 8
Prelude Debug.Traced> y
100000000
Prelude Debug.Traced> :t y
y :: Traced AsValue Integer
Prelude Debug.Traced> asExp y
10 * 10 * (10 * 10) * (10 * 10 * (10 * 10))
Prelude Debug.Traced> asSharedExp y
let _2 = 10 * 10; in  _2 * _2 * (_2 * (10 * 10))
Prelude Debug.Traced> :t asSharedExp y
asSharedExp y :: Traced AsExp Integer
Prelude Debug.Traced>

So what's going on? The value of x is Traced Integer, which means that there's some magic going on. The variable can be used as usual, for instance in computing x^8. A traced value can also be shown as an expression, which is what showAsExp does. So a traced value is somewhat like the symbolic values I had in an earlier post, but in addition to having a symbolic representation they also have a normal value. But the output from showAsExp doesn't really help in answering how many multiplications there are, since the shown expression has no sharing; it is totally flattened. The showAsShared function is the black magic here, it recovers the sharing, and we can see what happened. What we see is that there are actually five (5) multiplications involved in computing 10^8. This shows that the definition of exponentiation is suboptimal, since it can be done with three multiplications (three repeated squarings). The showAsShared really does have some black magic. It recovers information that is not part of the Haskell semantics, so from that we can conclude that it must contain the powerful incantation unsafePerformIO somewhere. How does it reveal implementation secrets? Look at this:

Prelude Debug.Traced> asSharedExp $ let a = 1 + 2 in a * a
let _1 = 1 + 2; in  _1 * _1
Prelude Debug.Traced> asSharedExp $ (1+2) * (1+2)
(1 + 2) * (1 + 2)
Prelude Debug.Traced>

The let expression and the expression where the variable has been expanded are semantically equal in Haskell, so no (pure) function can possibly be able to give different results for them.

OK, so how does it work? I'll show a simplified version of the Traced module here that only deals with one traced type, but it can be extended. The (soon to be available) hackage package contains the extended version.

In the Traced type we need to represent expressions. We only need constants and function applications.

data Traced a
    = Con a
    | Apply a String [Traced a]

The function application contains the value, the name of the function, and the arguments the function was applied to.

In the exported interface from the module we want to be able to convert to and from the Traced type. Nothing exciting here.

traced :: a -> Traced a
traced = Con

unTraced :: Traced a -> a
unTraced (Con x) = x
unTraced (Apply x _ _) = x

We want to show a traced value the same way we show the underlying value.

instance (Show a) => Show (Traced a) where
    show = show . unTraced

Comparing for equality, we simply compare the underlying values.

instance (Eq a) => Eq (Traced a) where
    x == y  =  unTraced x == unTraced y

And we'll make traced numbers be an instance of Num. All the functions (except fromInteger) build apply nodes.

instance (Num a) => Num (Traced a) where
    (+) = apply2 "+" (+)
    (-) = apply2 "-" (-)
    (*) = apply2 "*" (*)
    negate = apply1 "-" negate
    abs = apply1 "abs" abs
    signum = apply1 "signum" signum
    fromInteger = traced . fromInteger

apply1 s f x = Apply (f (unTraced x)) s [x]
apply2 s f x y = Apply (f (unTraced x) (unTraced y)) s [x, y]

A fancier version of this module could make the Traced type an applicative functor etc., but that's not really so important.

Finally, we want to be able to show a traced expression as an expression tree instead of a value.

 Traced t -> String
showAsExp (Con x) = show x
showAsExp (Apply _ s [x,y]) | not (isAlpha (head s)) =
    "(" ++ showAsExp x ++ " " ++ s ++ " " ++ showAsExp y ++ ")"
showAsExp (Apply _ s xs) =
    "(" ++ concat (intersperse " " $ s : map showAsExp xs) ++ ")"

We only export what is necessary, so the module header should be

module Traced(Traced, traced, unTraced, showAsExp) where

A quick test:

Traced> putStrLn $ showAsExp $ 10^8
(((10 * 10) * (10 * 10)) * ((10 * 10) * (10 * 10)))
Traced>

And now we need the black magic to recover the sharing. We would like to have a unique label in each node of the expression tree. If we only had that we could see when two things referred to the same subexpression, and use the label to refer to it instead of the value. If we were doing this in, e.g.., Java we could use object identity for this purpose. If we were doing it in C, we'd just compare pointers to the structs containing the expressions. But we're doing it in Haskell and none of this is available. It's not unavailable because Haskell wants to make our lives difficult, quite the contrary. Languages that allow pointer comparisons (object identity) must introduce an extra level of indirection in the semantics to explain how this is possible. So now it's not enough to know we have the number 5, we need to know that this is the number 5 at location 1000. And that's not the same as the number 5 at location 1010. The numbers contained in the locations might be the same, but the locations are not interchangable since they could, e.g., be mutated differently in the future.

So is everything lost in Haskell? Not at all. GHC implements a library of stable names, which is (at first approximation) the same as the address of something in memory. The API to System.Mem.StableName is very simple.

data StableName a
makeStableName :: a -> IO (StableName a)
hashStableName :: StableName a -> Int

The makeStableName function is like the & operator (address of) in C. It returns the "address" of something. So StableName is like a C pointer type. And the hashStableName function converts the "address" to an int, i.e., like casting it to int in C. (In the simplified code below we'll assume that two stable names never hash to the same Int, although this is not absolutely guaranteed.)

How come this interface is OK? For instance, calling makeStableName on semantically equal values can yield different results if the values happen to be stored in different parts of memory. It's ok, because the returned value is in the IO monad. In the IO monad anything can happen, so it's perfectly reasonable that the same argument yields different results in different calls.

Despite the name and the documentation the GHC stable names have a stability flaw. The stable name changes when an unevaluated expression is evaluated. It's annoying, but not a major flaw. But once evaluated the stable names is guaranteed to remain the same. (The implementation of stable names is not as simple as taking the address of an object, since the GC can move objects around.)

So now we have a way to get the identity of each node in the expression tree built by the Traced type. The plan is to traverse the expression tree. For each node we'll use the "address" of it as it's name, and remember that we've seen this node. As we traverse the tree we build a list of nodes we've seen. This list then corresponds to let bindings we'd like to display as the final result. As we traverse the nodes we'll also replace each node with a reference to its name instead, so we can see the sharing in the result.

To be able to represent the expression with rediscovered sharing we need to extend the Traced type. We need variable references and let bindings. In fact, we'll only generate a top level let binding, but we include it in the data type anyway.

data Traced a
    ...
    | Var a Name
    | Let [(Name, Traced a)] (Traced a)
type Name = String

We store the value in the Var constructor to make unTraced possible. New cases:

...
unTraced (Var x _) = x
unTraced (Let _ e) = unTraced e

And we want to show the new constructors:

...
showAsExp (Var _ n) = n
showAsExp (Let bs e) = "let " ++ concatMap bind bs ++ "in " ++ showAsExp e
  where bind (n, e) = n ++ " = " ++ showAsExp e ++ "; "

To rediscover the sharing we need to keep some state. We need a mapping from the node address to the Var that should replace it, and we need to accumulate the bindings, i.e., the pairs of node name and expression. So we need a state monad to keep track of the state. We also need to be able to call makeStableName in the IO monad, so we need IO as well. We'll do this by using the state transformer monad on top of IO. So the function that discovers sharing will take a traced value and return a new traced value, all in this monad. So we have the type:

type TState = (M.IntMap (Traced a), [(Name, Traced a)])
share :: Traced a -> StateT TState IO (Traced a)

Assuming some imports

import Control.Monad.State
import qualified Data.IntMap as M
import System.Mem.StableName

Now the body:

share e@(Con _) = return e
share e@(Apply v s xs) = do
    (sm, bs) <- get
    sn <- liftIO $ makeStableName e
    let h = hashStableName sn
    case M.lookup h sm of
        Just ie -> return ie
        Nothing -> do
            let n = "_" ++ show h
                ie = Var v n
            put (M.insert h ie sm, bs)
            xs' <- mapM share xs
            modify $ \ (sm', bs') -> (sm', (n, Apply v s xs') : bs')
            return ie

Constants are easy, we don't bother sharing them (we could, but it's not that interesting). For an apply node, we get its stable name (it's already evaluated, so it won't change) and hash it. We then grab the map (it's an IntMap; a fast map from Int to anything) from the state and look up the node. If the node is found we just return the expression from the map. If it's not in the map, we invent a name (using the hash value) and insert a Var node in the map so we'll never process this node again. We then recursively traverse all the children of the apply node, and rebuild a new apply node with those new children. This constitutes a binding and we stick it in the accumulated list of bindings. Finally we return the new Var node.

At the top level we need to call share and then build a let expression. The bindings are put in the list with the top node first in the list, so it looks nicer to reverse it.

shareIO :: Traced a -> IO (Traced a)
shareIO e = do
    (v, (_, bs)) <- runStateT (share e) (M.empty, [])
    return $ Let (reverse bs) v

We could leave it there, but to make it more convenient to use, we'll do an unsafePerformIO to hide the use of IO. This is not unsafe in the sense that it will make our program crash, but it is unsafe in the sense that it ruins the Haskell semantics. But since this whole exercise is to make a debugging/tracing tool this is a legitimate use, in my opinion.

reShare :: Traced a -> Traced a
reShare = unsafePerformIO . shareIO

showShared :: (Show a) => Traced a -> String
showShared = showAsExp . reShare

And let's test it:

Traced> putStrLn $ showShared $ 10^8
let _5 = (10 * 10); _7 = (_5 * _5); _11 = (10 * 10);
 _1 = (_5 * _11); _8 = (_7 * _1); in _8
Traced>

In my first example it looked a little prettier since the full implementation only shows bindings for nodes that are actually shared. This is a simple transformation to add. But looking at this, we can still see that there are five multiplications.

Hacking with types

It's a little annoying that we have both show and showAsExp to show traced values. It would be nicer if we could always use show and make the type determine how it is showed. We could invent a newtype to wrap around Traced and print the new type in a different way. But then this new type would not be compatible with the old one, which is a little annoying.

So we're going down a different road instead, we'll have a the Traced type have two parameters. The first one being a phantom type; just being used to determine how to show the traced value. We will rename the old type to TracedT, and it will only be used internally in the module.

newtype Traced t a = T { unT :: TracedT a }

data TracedT a
    = Con a
    | Apply a String [TracedT a]
    | Var a Name
    | Let [(Name, TracedT a)] (TracedT a)

Some minor changes are needed to the code to accomodate for this change.

traced :: a -> Traced t a
traced = T . Con

unTraced :: Traced t a -> a
unTraced = unTracedT . unT

unTracedT :: TracedT a -> a
...

apply1 s f x = T $ Apply (f (unTraced x)) s [unT x]
apply2 s f x y = T $ Apply (f (unTraced x) (unTraced y)) s [unT x, unT y]

And now for the fun part, the show function. We could do something like this:

data AsValue
data AsExp
instance Show (Traced AsValue a) ...
instance Show (Traced AsExp a) ...

This has problems, first it's not Haskell-98, but more importantly we can't print something of type Traced t a, because t is too general. We'd like to print as value by default, and be able to override this. There's only one defaulting mechanism in Haskell: the numeric defaulting. So as a somewhat disgusting hack, let's use the numeric defaulting to our advantage. We'll print Traced Integer a as values and Traced Double a as expressions.

instance (Num t, Show a) => Show (Traced t a) where
    show e = if doExp e then showAsExp (unT e) else show (unT e)

doExp :: (Num t) => Traced t a -> Bool
doExp x = '.' `elem` show (f x)
  where f :: (Num t) => Traced t a -> t
        f _ = 0

We distinguish between Integer and Double by how 0 is printed; for Double it has a '.' in the string.

A small utility to force printing as an expression.

asExp :: Traced t a -> Traced Double a
asExp = T . unT

Let's try it in ghci:

Traced> reShare $ 10^8
100000000
Traced> asExp $ reShare $ 10^8
let _8 = (10 * 10); _9 = (_8 * _8); _5 = (10 * 10);
 _7 = (_8 * _5); _10 = (_9 * _7); in _10
Traced> reShare $ (asExp 10)^8
let _7 = (10 * 10); _8 = (_7 * _7); _1 = (10 * 10);
 _5 = (_7 * _1); _9 = (_8 * _5); in _9
Traced>

In the first expression the numeric default made t be Integer, so we got a value. In the second and third case we forced t to be Double.

Some other examples:

Traced> let fac n = if n == 0 then 1 else n * fac(n-1)
Traced> asExp $ reShare $ fac 3
let _19 = (3 - 1); _24 = (_19 - 1); _18 = (_24 * 1);
 _20 = (_19 * _18); _21 = (3 * _20); in _21

Traced> let slowFib n = if n < 2 then 1 else slowFib(n-1) + slowFib(n-2)
Traced> asExp $ reShare $ slowFib 5
let _18 = (1 + 1); _19 = (_18 + 1); _23 = (1 + 1); _20 = (_19 + _23);
 _25 = (1 + 1); _17 = (_25 + 1); _21 = (_20 + _17); in _21

Traced> let fastFib n = fst $ iterate (\ (x,y) -> (y,x+y)) (1,1) !! n
Traced> asExp $ reShare $ fastFib 5
let _20 = (1 + 1); _21 = (1 + _20); _19 = (_20 + _21); _23 = (_21 + _19); in _23
Traced>

Well, that's enough abuse of the type system for one day.

More examples

The full library for traced values contains some more functionality, like naming values and functions that operate on traced booleans.

A named value, and a symbolic named value:

Prelude Debug.Traced> asSharedExp $ (named "x" 10)^8
let x = 10; _2 = x * x; in  _2 * _2 * (_2 * (x * x))
Prelude Debug.Traced> asSharedExp $ (unknown "x")^8
let _2 = x * x; in  _2 * _2 * (_2 * (x * x))

For a normal definition of fac we only see the arithmetic:

Prelude Debug.Traced> let fac n = if n == 0 then 1 else n * fac (n-1)
Prelude Debug.Traced> fac 5
120
Prelude Debug.Traced> asSharedExp $ fac 5
let _2 = 5 - 1;
    _4 = _2 - 1;
    _6 = _4 - 1;
in  5 * (_2 * (_4 * (_6 * ((_6 - 1) * 1))))

By using traced booleans we can see exactly what's going on:

Prelude Debug.Traced> let facT n = ifT (n %== 0) 1 (n * facT (n-1))
Prelude Debug.Traced> facT 5
120
Prelude Debug.Traced> asSharedExp $ facT 5
let _6 = 5 - 1;
    _11 = _6 - 1;
    _16 = _11 - 1;
    _21 = _16 - 1;
in  ifT (5 == 0) ...
        (5 * ifT (_6 == 0) ...
                 (_6 * ifT (_11 == 0) ...
                           (_11 * ifT (_16 == 0) ...
                                      (_16 * ifT (_21 == 0) ... (_21 * ifT (_21 - 1 == 0) 1 ...)))))

We can also lift a regular function to a traced function:

Prelude Debug.Traced> let fac' = liftFun "fac'" fac :: Traced t Integer -> Traced t Integer
Prelude Debug.Traced> let a = fac' 5 + fac' 10
Prelude Debug.Traced> a
3628920
Prelude Debug.Traced> asSharedExp a
fac' 5 + fac' 10
Prelude Debug.Traced>

And what about exponentiation?

The new exponentiation function in GHC looks like this:

(^) :: (Num a, Integral b) => a -> b -> a
x0 ^ y0 | y0 < 0    = error "Negative exponent"
        | y0 == 0   = 1
        | otherwise = f x0 y0
    where f x y | even y    = f (x * x) (y `quot` 2)
                | y == 1    = x
                | otherwise = g (x * x) ((y - 1) `quot` 2) x
          g x y z | even y = g (x * x) (y `quot` 2) z
                  | y == 1 = x * z
                  | otherwise = g (x * x) ((y - 1) `quot` 2) (x * z)

And if we try this definition we get

Debug.Traced> asSharedExp $ 10 ^ 8
let _2 = 10 * 10; _1 = _2 * _2; in  _1 * _1
Debug.Traced>

Victory!

Edit: Package available in hackage.

2008-03-09T21:45:00.007Z

Simple reflections of higher order

In a recent blog post by Twan van Laarhoven he showed how to reflect Haskell expressions so that we can actually see them symbolically. This has now been included in lambdabot on the Haskell IRC (source). Let's play with it for a moment:

*SimpleReflect> x+y
x + y
*SimpleReflect> foldr f z [1,2,3]
f 1 (f 2 (f 3 z))
*SimpleReflect> let swap (x,y) = (y, x)
*SimpleReflect> swap (a,b)
(b,a)
*SimpleReflect> map swap [(a,b),(c,d)]
[(b,a),(d,c)]
*SimpleReflect> :t x
x :: Expr

That's very cool! Read Twan's post to find out how it works. All we need to know is on the last line, x :: Expr. So Expr is a special type with special instances. Let's try something else

*SimpleReflect> \ x -> x + y

:1:0:
    No instance for (Show (Expr -> Expr))
    ...

Oh, that's annoying, we can't print functions. But why can't we? If we made up some variable of type Expr and then applied the function we'd get something we could print. Let's add some code to Twan's module. I've just randomly picked the variable a. (To make this compile you need the extension language FlexibleInstances.)

instance Show (Expr -> Expr) where
    show f = "\\ " ++ show a ++ " -> " ++ show (f a)

And try again with the example

*SimpleReflect> \ x -> x + y
\ a -> a + y

Pretty smooth. Except it doesn't really work.

*SimpleReflect> \ x -> x + a
\ a -> a + a

Those two functions are not the same. We can't just arbitrarily pick the variable a since it might be free in the expression.

So we need to pick a variable that is not free in the expression. How could we do that? No matter what we pick it could be used. And we have no idea what is being used. Or do we? If we just turn the function in to text we could look at the string, and pick some variable that is unused in that string. This is really a gruesome hack, but who cares?

How do we print the function? Just invent some variable, I use _, turn it into an expression and tokenize the string. We will then have a string of tokens not to use. To find a variable to use, just pick an infinite supply of variables, remove the ones being used, and then pick one of the remaining ones.

instance Show (Expr -> Expr) where
    show f = "\\ " ++ show v ++ " -> " ++ show (f v)
      where v = var (head vars)
            vars = supply \\ tokenize (show $ f $ var "_")
            supply = [ "x" ++ show i | i <- [1..]]
            tokenize "" = []
            tokenize s = case lex s of (x,s') : _ -> x : tokenize s'

And try to fool it:

*SimpleReflect> \ x -> x + y
\ x1 -> x1 + y
*SimpleReflect> \ x -> x + var "x1"
\ x2 -> x2 + x1

OK, what about multiple arguments?

*SimpleReflect> \ x y -> x + y + z

:1:0:
    No instance for (Show (Expr -> Expr -> Expr))
    ...

Well, yeah, that's true. There is no such instance. But wait, why is the instance we have for the type Expr->Expr? No particular reason, that's just what I wrote. In fact it works equally well for Expr->r as long as we can show r, because that's the only thing we do with f v. So we change the first line of the instance:

instance (Show r) => Show (Expr -> r) where

And now

*SimpleReflect> \ x y -> x + y + z
\ x2 -> \ x1 -> x2 + x1 + z
*SimpleReflect> foldr (.) id [f::Expr->Expr,g,h] -- a little type help needed
\ x1 -> f (g (h x1))
*SimpleReflect> scanr (.) id [(*2),f::Expr->Expr,(+1)]
[\ x1 -> f (x1 + 1) * 2,\ x1 -> f (x1 + 1),\ x1 -> x1 + 1,\ x1 -> x1]

Well, that wasn't too hard. So let's try another example.

*SimpleReflect> \ (x,y) -> x+y+z

:1:0:
    No instance for (Show ((Expr, Expr) -> Expr))
    ...

Hmmm, yes, the argument must be an Expr. That's annoying. We need to generalize the argument type. We want be able to use Expr and pair of them etc. Time for a type class. What does it need to do? It has to invent expressions and never reuse variables doing so. So when we invent an Expr we need to consume one variable and leave the rest for others to consume. This sounds like a state monad. So we're going to use a state monad where the state is the (infinite) list of strings that are available for making variables.

instance (Show a, ExprArg a, Show r) => Show (a -> r) where
    show f = "\\ " ++ show v ++ " -> " ++ show (f v)
      where v = evalState exprArg vars
            dummy = evalState exprArg $ repeat "_"
            vars = supply \\ tokenize (show $ f dummy)
            supply = [ "x" ++ show i | i <- [1..]]
            tokenize "" = []
            tokenize s = case lex s of (x,s') : _ -> x : tokenize s'

class ExprArg a where
    exprArg :: State [String] a

instance ExprArg Expr where
    exprArg = do v:vs <- get; put vs; return (var v)

Using this we're back where we were before, but now we can make some more instances.

instance ExprArg () where
    exprArg = return ()

instance (ExprArg a, ExprArg b) => ExprArg (a, b) where
    exprArg = liftM2 (,) exprArg exprArg

instance (ExprArg a, ExprArg b, ExprArg c) => ExprArg (a, b, c) where
    exprArg = liftM3 (,,) exprArg exprArg exprArg

And finally:

*SimpleReflect> \ (x, y) -> x + y + z
\ (x1,x2) -> x1 + x2 + z
*SimpleReflect> curry f :: Expr -> Expr -> Expr
\ x2 -> \ x1 -> f (x2,x1)
*SimpleReflect> uncurry f :: (Expr, Expr) -> Expr
\ (x1,x2) -> f x1 x2
*SimpleReflect> \ () -> 1
\ () -> 1

The last example is a curiosity since it does not involve Expr at all. Well, that's enough for a Sunday night hack.

2007-11-09T18:26:00.000Z

Some lambda calculus examples

Syntax

In a previous blog entry I described a simple evaluator and type checker for the lambda cube, i.e., various forms of lambda calculus.

Here I'm going to show some examples of code in pure typed λ-calculus. All the examples are typable in F_ω; the full lambda cube is not necessary.

Before doing any examples we'd better have some syntax that is not too painful, because writing λ-expression in raw Haskell is tedious. The syntax for variables, * kind, and application is easy, I'll just use Haskell syntax. For λ-expressions the Haskell syntax doesn't allow explicit type annotations, but various Haskell compiler implement that extension, so I'll just pick that. So a λ will be written, "\ (var::type) -> expr". And as in Haskell I'll allow multiple variables between the "\" and "->"; it's just a shorthand for multiple lambdas.

So what about the dependent function type? The syntax (x:t)→u suggests (x::t)->u, so I'll use that. And when the variable doesn't occur we'll write t->u as usual. For type variables Haskell (well, not Haskell 98, but extensions) uses forall (a::*) . t, so I'll allow that too.

An example, the identity function:

\ (a::*) (x::a) -> x

with type

forall (a::*) . a->a

And using it

id Int 5

Writing a pretty printer and parser for this is pretty straight forward so I'll skip that and just point you to the code. BTW, instead of using Parsec for the parser like everyone else I used ReadP. The ReadP library is very nice, partly because the alternative operator is actually commutative (unlike Parsec). But the error messages suck.

Enter `let`

Now if we want to use, say, the identity function more than once we need to name it. There is a mechanism for that, namely λ. But it looks awkward. Look:

(\ (id :: forall (a::*) . a->a) -> ... id ... id ... id ...) (\ (a::*) (x::a) -> x)

What makes it awkward is that the name, id, is far away from the body, \ (a::*) (x::a) -> x. From Haskell we are more used the let and where expressions. So let's add that let. Instead of what we have above we'll write

let id :: forall (a::*) . a->a = \ (a::*) (x::a) -> x
in  ... id ... id ... id ...

The let construct could be just be syntactic suger for a lambda and an application, but I've decided to add it as a constructor to the expression type instead.

    | Let Sym Type Expr Expr

Adding a new constructor means that we have to modify all the functions operating on Expr, and it's extra much work because Let is a variable binding construct. For substitution we just cheat a little and use the expansion into an application and λ

    sub (Let i t e b) = let App (Lam i' t' b') e' = sub (App (Lam i t b) e)
                        in  Let i' t' e' b'

And for normal form we use the same trick.

    spine (Let i t e b) as = spine (App (Lam i t b) e) as

And for now, let's extend the type checker in the same way.

To make definitions a little less verbose I'll allow multiple bindings. E.g.

let x :: Int = g a; y :: Int = f x in h x y

It's just multiple nested Lets in the obvious way. Another shorthand. I'll allow the identity function to be written

let id (a::*) (x::a) :: a = x
in ...

The translation is pretty easy

let f (x1::t1) ... (xn::tn) :: t = b in e

let f :: forall (x1::t1) ... (xn::tn) . t = \ (x1::t1) ... (xn::tn) -> b in e

And finally, to make it even more like Haskell, I'll allow the type signature to be on a line of its own, omitting types on the bound variables.

let id :: forall (a::*) . a -> a;
    id a x = x;
in  ...

It's pretty easy to translate to the bare expression type. Why all these little syntactic extras? Well, remember Mary Poppins "A spoonful of sugar makes the medicine go down." Phew! Enough syntax, on to the examples.

Bool

We could start by throwing in all kinds of primitive types into our little language, but who knows what might happen then. All the nice properties that have been shown about the language might not hold anymore. So instead we'll code all the types we need with what we already have.

The Bool type has two values: False and True. So we need to find a type that has exactly two different values. Fortunately that's easy: a->a->a (I'll be a bit sloppy and leave off top level quantifier sometimes, just like Haskell; That should be forall (a::*) . a->a->a).

Why does that have two possible values? Well, we have a function type that must return an a for any possible a that it's given, so it can't conjure up the return value out of thin air. It has to return one of it's arguments. And there are two arguments to choose from, so there are two different values in that type. Here's Bool, False, True:

Bool :: *;
Bool = forall (a::*) . a->a->a;

False :: Bool;
False = \ (a::*) (x::a) (y::a) -> x;

True  :: Bool;
True = \ (a::*) (x::a) (y::a) -> y;

If the definition for Bool was written in Haskell it would be

type Bool = forall a . a->a->a

false :: Bool
false = \ x y -> x

true :: Bool
true = \ x y -> y

It would be easy to add more sugar and allow type which just means you're defining something of type *. And you could also make an omitted type in a forall mean type *. But I've not added these little extras.

Defining the if function is trivial; it's just a matter of permuting the arguments a little, because the boolean values come with the "if" built in.

if :: forall (a::*) . Bool -> a -> a -> a;
if a b t f = b a t f;

Note how the type signature is exactly the same as you'd find in Haskell. A difference is that we have explicit type abstraction and type application. For instance, if takes a first argument that is the type of the branches. So when using if we must pass in a type.

Given this we can try to type check some simple code:

let Bool :: *;
    Bool = forall (a::*) . a->a->a;
    False :: Bool;
    False = \ (a::*) (x::a) (y::a) -> x;
in  False

Here is what happens:

*CubeExpr> typeCheck $ read "let Bool :: *; Bool = forall (a::*) . a->a->a; False :: Bool;False = \\ (a::*) (x::a) (y::a) -> x; in False"
*** Exception: Type error:
Bad function argument type:
Function: \ (False :: Bool) -> False
argument: \ (a :: *) (x :: a) (y :: a) -> x
expected type: Bool
     got type: forall (a :: *) . a->a->a

What is it whining about? Expected Bool, got forall (a :: *) . a->a->a. But it says right there in the code that Bool is exactly that. What's going on?
Well, the type checker knows the type of everything, but not the value of anything. So the type checker knows that type of Bool is *, but it doesn't know what it is equal to.

The problem is that when we have a let binding we know the value of the defined variable and to be able to do dependent type checking the type checker needs to know it too. We need to change the type checking of let. Here's a simple solution:

tCheck r (Let s t a e) = do
    tCheck r t
    ta <- tCheck r a
    when (not (betaEq ta t)) $ throwError $ "Bad let def\n" ++ show (ta, t)
    te <- tCheck r (subst s a e)
    tCheck r (Pi s t te)
    return te

Note how we substitute the value of the definition (a) in the body before type checking it. This is a sledge hammer approach. A better one would be to change the environment in the type checker to carry values when they are known. These values could then be used when computing the normal form of an expression. But let's keep it simple for now.

Let's try again:

*CubeExpr> typeCheck' $ read "let Bool :: *; Bool = forall (a::*) . a->a->a; False :: Bool;False = \\ (a::*) (x::a) (y::a) -> x; in False"
forall (a :: *) . a->a->a

That worked!

A simple top level

To make experimentation easier I've added a simple top level where you can evaluate expression, make definitions, load files, etc. The files are just a bunch of definitions the way they would go inside a let.

Sample session:

Welcome to the Cube.
Use :help to get help.
Cube> \ (a::*) (x::a) -> x
\ (a :: *) (x :: a) -> x
  ::
forall (a :: *) . a->a
Cube>

When evaluating an expression it will first print the value, then ::, and finally the type.
Let's try some more. Here's the file bool.cube.

-- The Bool type.

Bool :: *;
Bool = forall (boolT::*) . boolT->boolT->boolT;

False :: Bool;
False = \ (boolT::*) (false::boolT) (true::boolT) -> false;

True  :: Bool;
True = \ (boolT::*) (false::boolT) (true::boolT) -> true;

if :: forall (a::*) . Bool -> a -> a -> a;
if a b t f = b a f t;

not :: Bool -> Bool;
not b = if Bool b False True;

and :: Bool -> Bool -> Bool;
and x y = if Bool x y False;

or :: Bool -> Bool -> Bool;
or x y = if Bool x True y;

Note how I've named some of the type variables with more suggestive names. More about that in a moment.

Cube> :load bool.cube
Cube> True
\ (boolT :: *) (false :: boolT) (true :: boolT) -> true
  ::
forall (boolT :: *) . boolT->boolT->boolT
Cube> and True False
\ (boolT :: *) (false :: boolT) (true :: boolT) -> false
  ::
forall (boolT :: *) . boolT->boolT->boolT

By staring carefully at these you can convince yourself that it's right. What's making it a little hard to read are all those leading lambdas (and the corresponding function types). But there's an option to suppress the printing of them.

Cube> :skip
Cube> True
true
  ::
boolT
Cube> and True False
false
  ::
boolT

Look, it got deceptively readable.

Characters

So we can define booleans. What about something like characters? You can define those too. What is the ASCII type? It's a type with 128 different values. We saw how we can make a type with two values (Bool); we can do the same for 128. Since that takes a lot of room, I'll do it for just four values.

Char :: *;
Char = forall (charT::*) . charT->charT->charT->charT->charT;

'a' :: Char;
'a' = \ (charT::*) (_'a'::charT) (_'b'::charT) (_'c'::charT) (_'d'::charT) -> _'a';
'b' :: Char;
'b' = \ (charT::*) (_'a'::charT) (_'b'::charT) (_'c'::charT) (_'d'::charT) -> _'b';
'c' :: Char;
'c' = \ (charT::*) (_'a'::charT) (_'b'::charT) (_'c'::charT) (_'d'::charT) -> _'c';
'd' :: Char;
'd' = \ (charT::*) (_'a'::charT) (_'b'::charT) (_'c'::charT) (_'d'::charT) -> _'d';

eqChar :: Char -> Char -> Bool;
eqChar x y = x Bool (y Bool True False False False)
                    (y Bool False True False False)
                    (y Bool False False True False)
                    (y Bool False False False True);

Since single quotes are allowed in identifiers in my little language the 'a' is just an identifier as usual. Some sample use:

Cube> :load char.cube
Cube> 'a'
_'a'
  ::
charT
Cube> eqChar 'a' 'a'
true
  ::
boolT
Cube> eqChar 'a' 'b'
false
  ::
boolT

A different way to use characters is just to assume that there is some character type and that is has some values. How can we do that? By using a lambda expression.

\ (Char :: *) -> \ ('a' :: Char) ('b' :: Char) ... -> ...

One final little syntactic sugar is to allow lambda expressions to be written with let. So

let x :: T;
in  e

is the same as

\ (x :: T) -> e

So for Char we can write

Char :: *;
'a' :: Char;
'b' :: Char;
'c' :: Char;
'd' :: Char;
eqChar :: Char -> Char -> Bool;

This is a very handy notation during program development, btw. Inside i let you can first just write the type of something and then use it freely. This something will then be lambda abstracted until such a time that you provide the definition.

Let's try the "fake" characters.

Cube> :load extchar.cube
Cube> 'a'
'a'
  ::
Char
Cube> :skip
Cube> 'a'
\ (Char :: *) ('a' :: Char) ('b' :: Char) ('c' :: Char) ('d' :: Char)
  (eqChar :: Char->Char->forall (boolT :: *) . boolT->boolT->boolT)
  -> 'a'
  ::
forall (Char :: *) .
       Char->Char->Char->Char->
       (Char->Char->forall (boolT :: *) . boolT->boolT->boolT)->
       Char
Cube> :skip
Cube> eqChar 'a' 'b'
eqChar 'a' 'b'
  ::
forall (boolT :: *) . boolT->boolT->boolT

Now eqChar no longer gets evaluated, naturally, because it has no definition; it's λ bound.

Now we know how the Char type could be implemented (or supplied from the outside). It would then be all right to add it as a primitive (with the same behavior), because it will not ruin any properties. (Speaking of ruining properties, the seq function in Haskell is an example of a function that cannot be defined in the λ-calculus. And sure enough, adding it to Haskell ruined the validity of η-reduction.)

It's not difficult to extend the Expr type with some Prim constructor for primitive functions, and adding some special cases in the syntax and type checking for them. But since this is not intended to be a usable language I'll resist the temptation (for now).

Pairs

OK, booleans and characters were pretty easy. Let's do pairs next. The Pair type is parameterized over two types; the type of the first and the second component. As with booleans the representation of pairs come with the case analysis "build in". (Case analysis on pairs is flip uncurry in Haskell.) If we can correctly implement building pairs, fst, and snd we are done.

Pair :: * -> * -> *;
Pair a b = forall (pairT::*) . (a->b->pairT) -> pairT;

PairC :: forall (a::*) (b::*) . a -> b -> Pair a b;
PairC a b x y = \ (pairT::*) (pair :: a->b->pairT) -> pair x y;

fst :: forall (a::*) (b::*) . Pair a b -> a;
fst a b p = p a (\ (x::a) (y::b) -> x);

snd :: forall (a::*) (b::*) . Pair a b -> b;
snd a b p = p b (\ (x::a) (y::b) -> y);

So the type is called Pair and the constructor PairC. Since they live in the same name space we can't use Haskell's convention of giving them the same name (which is (,) in Haskell.).

Again, in Haskell the first definition would be

type Pair a b = forall pairT . (a->b->pairT) -> pairT

Staring at these definitions makes it obvious(?) that they work. But we can test it too.

Cube> :let S :: *; s :: S; T :: *; t :: T
Cube> :load pair.cube
Cube> :let p :: Pair S T; p = PairC S T s t
Cube> p
pair s t
  ::
pairT
Cube> fst S T p
s
  ::
S
Cube> snd S T p
t
  ::
T

The :let command is used to make bindings without having to load them from a file. The first :let just introduces some values to play with.

Again, the big difference compared to Haskell is that all types are explicit. (But users of templates in C++ or generics in Java/C# might appeciate it?)

Maybe

We can follow the same pattern and define Haskell's Maybe type. In fact any non-recursive Haskell data type has a mechanical translation and easy into lambda calculus. We've seen boolean, characters, and pairs above.

Maybe :: * -> *;
Maybe a = forall (maybeT::*) . maybeT->(a->maybeT)->maybeT;

Nothing :: forall (a::*) . Maybe a;
Nothing a = \ (maybeT::*) (nothing::maybeT) (just::a->maybeT) -> nothing;

Just :: forall (a::*) . a -> Maybe a;
Just a x = \ (maybeT::*) (nothing::maybeT) (just::a->maybeT) -> just x;

maybe :: forall (a::*) (r::*) . r -> (a->r) -> Maybe a -> r;
maybe a r nothing just s = s r nothing just;

Natural numbers

Let's carry on with some natural numbers. Natural numbers are trickier because it's a recursive type. In Haskell:

data Nat = Zero | Succ Nat

Let's try the cube version:

Cube> :load nat.cube
Cube> 2
succ (succ zero)
  ::
natT
Cube> add 2 3
succ (succ (succ (succ (succ zero))))
  ::
natT

What does the definitions for natural numbers look like?

Nat :: *;
Nat = forall (natT::*) . natT -> (natT->natT) -> natT;

0 :: Nat;
0 = \ (natT::*) (zero::natT) (succ::natT->natT) -> zero;

Succ :: Nat -> Nat;
Succ n = \ (natT::*) (zero::natT) (succ::natT->natT) -> succ (n natT zero succ);

natprim :: forall (r::*) . (r->r) -> r -> Nat -> r;
natprim r succ zero n = n r zero succ;

add :: Nat -> Nat -> Nat;
add x y = x Nat y Succ;

mul :: Nat -> Nat -> Nat;
mul x y = x Nat 0 (add y);

power :: Nat -> Nat -> Nat;
power x y = y Nat (Succ 0) (mul x);

isZero :: Nat -> Bool;
isZero n = n Bool True (\ a::Bool -> False);

1 :: Nat;
1 = Succ 0;
2 :: Nat;
2 = Succ 1;
3 :: Nat;
3 = Succ 2;

The natprim function corresponds to foldr for list and is our means of recursion. This type for natural numbers come with primitive recursion built in, just as the non-recursive data types earlier came with case analysis built in.

I'll skip defining subtraction&co, they are possible, but a little tedious (and very inefficient).

Lists

Lists are like a hybrid of natural numbers and the Maybe type. They follow the same pattern as before.

List :: * -> *;
List e = forall (listT::*) . listT -> (e->listT->listT) -> listT;

Nil :: forall (e::*) . List e;
Nil e = \ (listT::*) (nil::listT) (cons::e->listT->listT) -> nil;

Cons :: forall (e::*) . e -> List e -> List e;
Cons e x xs = \ (listT::*) (nil::listT) (cons::e->listT->listT) -> cons x (xs listT nil cons);

And some sample functions.

foldr :: forall (a::*) (b::*) . (a->b->b) -> b -> List a -> b;
foldr a b f z xs = xs b z f;

map :: forall (a::*) (b::*) . (a->b) -> List a -> List b;
map a b f xs = foldr a (List b) (\ (x::a) (r::List b) -> Cons b (f x) r) (Nil b) xs;

append :: forall (a::*) . List a -> List a -> List a;
append a xs ys = foldr a (List a) (Cons a) ys xs;

foldl :: forall (a::*) (b::*) . (b->a->b) -> b -> List a -> b;
foldl a b f z xs = foldr a (b->b) (\ (x::a) (g::b->b) (v::b) -> g (f v x)) (\ (x::b) -> x) xs z;

reverse :: forall (a::*) . List a -> List a;
reverse a xs = foldl a (List a) (\ (r :: List a) (x :: a) -> Cons a x r) (Nil a) xs;

And a test:

Welcome to the Cube.
Use :help to get help.
Cube> :load bool.cube
Cube> :load nat.cube
Cube> :load list.cube
Cube> :load pair.cube    
Cube> :load listmisc.cube
Cube> :load extchar.cube
Cube> :skip
Cube> replicate Char 3 'b'
cons 'b' (cons 'b' (cons 'b' nil))
  ::
listT
Cube> :quit

And with that, I'll quit too for now.

2007-11-06T16:28:00.000Z

Benchmarking ray tracing, Haskell vs. OCaml

On his web site about OCaml Jon Harrop has a benchmark for a simple ray tracing program written in a number of languages. When I saw it I wondered why Haskell was doing so badly, especially since this benchmark was taken as some kind of proof that Haskell performs poorly in "real life". So I have rerun the benchmarks. I've also rewritten the Haskell versions of the programs. The Haskell versions on Jon's web site were OK, but they were too far from the OCaml versions for my taste. I prefer to keep the programs very similar in a situation like this. My rewrite of the benchmarks from OCaml to Haskell was done with a minimum of intelligence. Here are the only things I did that needed creative thought:

Use the type Vector to represent vectors instead of a tuple. This allows the components to be strict.
Use the type Scene instead of a tuple to represent a scene. The tuple used in the OCaml code uses the dubious feature of equi-recursive types (even Xavier thinks it's strange enough to have a flag to enable it).
Rewrite the loop that computes a pixel's value using an accumulating updatable variable into a list comprehension that sums the list.
Finally, the compiler flags needed a bit of tweaking to get good performance, even though "-O3 -funbox-strict-fields -fexcess-precision -optc-ffast-math" were pretty obvious.

In addition to this I made the code look a little more Haskellish, e.g., using overloading to allow + and - on vectors. This is really just minor syntactic changes, but makes the code more readable.

To make the program size comparison fair I removed some dead code from the OCaml code.

I then reran the benchmarks using Haskell, OCaml, and C++.

The benchmarks are five programs that starts from very simple ray tracing and specializing the program more and more to speed it up.

The numbers are in the tables below. The time is execution time in second, the characters are non-white characters in the file, and the lines are the number of lines in the file. To ease comparison I also include the relative numbers compared to OCaml (smaller numbers are better).

Interestingly, and unlike Jon's benchmark, the Haskell code is always smaller than the OCaml code. Furthermore, the Haskell code ranges from much faster to slightly faster than the OCaml code. Again, this is very unlike Jon's benchmark. I find the unoptimized version of the benchmark especially interesting since Haskell is 5 times(!) faster than OCaml. I've not investigated why, but I suspect laziness.

Results

The programs, ray1-ray5, are variations on the ray tracer as given on Jon's web site. I've used the same size metrics as Jon does.

Haskell: My Haskell code compiled with ghc 6.8.1
Haskell old: Jon's Haskell code, compiled with ghc 6.8.1
Haskell old 6.6: Jon's Haskell code, compiled with ghc 6.1.1
OCaml: Jon's OCaml code
C++: Jon's C++ code

Time: execution time is second
Char: number of non-white chracters in the program
Lines: number of lines in the program
Rel T: execution time relative to OCaml
Rel C: non-white characters relative to OCaml
Rel L: lines relative to OCaml
Mem: Maximum resident memory

ray1	Time	Chars	Lines	Rel T	Rel C	Rel L	Mem
Haskell	15.3	1275	51	0.202	0.990	1.020	5M
Haskell, old	15.8	1946	88	0.208	1.511	1.760	9M
Haskell, old 6.6	28.1	1946	88	0.370	1.511	1.760	9M
OCaml	75.9	1288	50	1.000	1.000	1.000	18M
C++	8.1	2633	122	0.106	2.044	2.440	8M

ray2	Time	Chars	Lines	Rel T	Rel C	Rel L	Mem
Haskell	11.5	1457	50	0.206	0.912	0.943	12M
Haskell, old	12.0	2173	99	0.215	1.360	1.868	35M
Haskell, old 6.6	21.1	2173	99	0.379	1.360	1.868	35M
OCaml	55.8	1598	53	1.000	1.000	1.000	15M
C++	6.1	3032	115	0.108	1.897	2.170	8M

ray3	Time	Chars	Lines	Rel T	Rel C	Rel L	Mem
Haskell	9.7	1794	62	0.970	0.919	0.939	12M
Haskell, old	11.1	2312	103	1.112	1.184	1.561	35M
Haskell, old 6.6	19.7	2312	103	1.984	1.184	1.561	35M
OCaml	10.0	1953	66	1.000	1.000	1.000	15M
C++	5.4	3306	143	0.545	1.693	2.167	8M

ray4	Time	Chars	Lines	Rel T	Rel C	Rel L	Mem
Haskell	8.5	1772	66	0.985	0.867	0.957	12M
Haskell, old	11.7	2387	110	1.360	1.168	1.594	36M
Haskell, old 6.6	19.2	2387	110	2.235	1.168	1.594	35M
OCaml	8.6	2043	69	1.000	1.000	1.000	11M
C++	5.0	3348	149	0.584	1.639	2.159	8M

ray5	Time	Chars	Lines	Rel T	Rel C	Rel L
Haskell	7.0	2246	95	0.999	0.878	0.950
OCaml	7.0	2559	100	1.000	1.000	1.000
C++	4.7	3579	142	0.674	1.399	1.420

The source code is available in a Darcs repository.

Software and hardware details

Hardware: MacBook, Intel Core Duo 2GHz, 2MB L2 Cache, 1GB 667MHz DRAM

Software:

Haskell compiler: ghc-6.8.1
OCaml compiler: 3.10.0
g++: gcc version 4.0.1 (Apple Computer, Inc. build 5367)

Compilation commands:

ghc: -O3 -fvia-C -funbox-strict-fields -optc-O3 -fexcess-precision -optc-ffast-math -funfolding-keeness-factor=10
OCaml: -rectypes -inline 100 -ffast-math -ccopt -O3
g++: -O3 -ffast-math

Target architecture is x86 (even though the processor is x86_64 capable).

Some observations

Haskell should really have the definitions of infinity and epsilon_float in a library. They are quite useful. Also, having them in a library would have made the Haskell code somewhat shorter and faster.

Converting these programs from OCaml to Haskell was very mechanical; it could almost be done with just sed.

I'm glad version 5 of the benchmark didn't show much improvement, because it's a really ugly rewrite. :)

Note that the code is all Haskell98, no strange extensions (even though -funbox-strict-fields deviates subtly from H98).

Conclusion

Benchmarking is tricky. I'm not sure why my and Jon's numbers are so different. Different hardware, slightly different programs, different software.

Haskell is doing just fine against OCaml on this benchmark; the Haskell programs are always smaller and faster.

Edit: Updated tables with more numbers.

PS: Phil Armstrong wrote the Haskell code on Jon's web site and I took some code from his original.

2007-10-25T16:11:00.001+01:00

Simpler, Easier!

In a recent paper, Simply Easy! (An Implementation of a Dependently Typed Lambda Calculus), the authors argue that type checking a dependently typed language is easy. I agree whole-heartedly, it doesn't have to be difficult at all. But I don't think the paper presents the easiest way to do it. So here is my take on how to write a simple dependent type checker. (There's nothing new here, and the authors of the paper are undoubtedly familiar with all of it.)

First, the untyped lambda calculus.

I'll start by implementing the untyped lambda calculus. It's a very simple language with just three constructs: variables, applications, and lambda expressions, i.e.,

x e e λx.e

For example, (λx.λy.x)(λz.z). In Haskell I'll use strings to represent variables names; it's simple and easy.

type Sym = String

data Expr
        = Var Sym
        | App Expr Expr
        | Lam Sym Expr
        deriving (Eq, Read, Show)

The example above represented by App (Lam "x" $ Lam "y" $ Var "x") (Lam "z" $ Var "z"). What do we want to do with the Expr type? Well, evaluating an expression seems like the thing we need. Now, there are many degrees of evaluation to choose from, Weak Head Normal Form, Head Normal Form, Normal Form, etc., etc. They differ in exactly where there might be reducible expression lingering. To evaluate lambda expression the most important step is β-reduction. A β-reduction step can be performed anywhere a function meets an argument, i.e., an application where the function is on λ form, a.k.a. a redex.

(λx.e)a reduces to e[x:=a]

Where the e[x:=a] notation means that all (free) occurrences of the variable x in the expression e are replaced by a. The example above has one redex, and doing a β step yields λy.λz.z. The other kind of reduction we will make use of is α-substitution, which is simply renaming a bound variable, e.g., λx.x can be changed to λy.y. Let's start with an easy evaluation strategy, normal order to WHNF. In WHNF we only need to make sure there the there's no redex along the "spine" of the expression, i.e., starting from the root and following the left branch in applications. Doing normal order reduction means that we do not evaluate anything inside the argument of a β redex before doing the reduction. It's sometimes called lazy evaluation, but I prefer to use that term for an implementation strategy for normal order reduction. We implement normal order WHNF by walking down the spine collecting arguments (i.e., the right branch of a applications) until we reach a lambda or a variable. If we reach a variable we already have WHNF so we just reconstitute the applications again. If we hit a lambda we get to the crux of evaluation. We need to perform a β-reduction, i.e., if we have App (Lam v e) a we need to replace all (free) occurrences of the variable v by the argument a inside the lambda body e. That's what the subst function does.

whnf :: Expr -> Expr
whnf ee = spine ee []
  where spine (App f a) as = spine f (a:as)
        spine (Lam s e) (a:as) = spine (subst s a e) as
        spine f as = foldl App f as

The subst function is the only tricky part, so let's relax by first defining something easy, namely getting the free variables from an expression. The free variables are those variables that occur in an expression, but are not bound in it. We simply collect the variables in a set (using a list as a set here), removing anything bound.

freeVars :: Expr -> [Sym]
freeVars (Var s) = [s]
freeVars (App f a) = freeVars f `union` freeVars a
freeVars (Lam i e) = freeVars e \\ [i]

Back to substitution.

subst :: Sym -> Expr -> Expr -> Expr
subst v x b = sub b
  where sub e@(Var i) = if i == v then x else e
        sub (App f a) = App (sub f) (sub a)
        sub (Lam i e) =
            if v == i then
                Lam i e
            else if i `elem` fvx then
                let i' = cloneSym e i
                    e' = substVar i i' e
                in  Lam i' (sub e')
            else
                Lam i (sub e)
        fvx = freeVars x
        cloneSym e i = loop i
           where loop i' = if i' `elem` vars then loop (i ++ "'") else i'
                 vars = fvx ++ freeVars e

substVar :: Sym -> Sym -> Expr -> Expr
substVar s s' e = subst s (Var s') e

The subst function will replace all free occurrences of v by b inside x, i.e., b[v:=x]. The Var case is easy. If it's the variable we are replacing then replace else leave it alone. The App case is also easy, just recurse in both branches. The Lam case has three alternative. First, if the bound variable is the same as the one we are replacing then there can be no free occurrences inside it, so just return the lambda as is. Second, if the lambda bound variable is among the free variables in x we have a problem (see below). Third case, just recurse in the body. So, what about the case when the lambda bound variable occurs in x? Well, if we just blindly continue substitution the variable v inside x will no longer refer to the same thing; it will refer to the variable bound in the lambda. That's no good. For example, take the expression λx.((λy.λx.y)x), the β reduction gives λx.λx'.x (or similar), whereas doing the substitution blindly would give λx.λx.x. Which is wrong! But it's easy to fix, just conjure up a variable, i' that will not clash with anything (cloneSym does that). How do we come up with a good variable? Well, we don't want to pick one that is free in the expression x because that would lead to the same problem again. Nor do we want to pick a variable that is free in e because that would accidentally bind something in e. So we take the original identifier and tack on "'" until it fulfills our requirements. (OK, efficiency affectionados are allowed to complain a little here, but this isn't that bad actually.) Once we have a new variable we can do an α-conversion to rename the offending variable to something better. The substVar function is a utility when we want to replace one variable with another. Another useful thing to be able to do is to compare lambda expression for equality. We already have syntactic equality derived for Expr, but it is also very useful to be able to compare expressions modulo α-conversions. That is, we'd like λx.x to compare equal to λy.y. Let's call that function alphaEq.

alphaEq :: Expr -> Expr -> Bool
alphaEq (Var v)   (Var v')    = v == v'
alphaEq (App f a) (App f' a') = alphaEq f f' && alphaEq a a'
alphaEq (Lam s e) (Lam s' e') = alphaEq e (substVar s' s e')
alphaEq _ _ = False

Variables and applications just proceed along the structure of the expression. When we hit a lambda the variables might be different, so we do an α-conversion of the second expression to make them equal. As the final functions, we will do reduction to Normal Form (i.e., to a form where no redexes remain) and comparison of expressions via their normal forms.

nf :: Expr -> Expr
nf ee = spine ee []
  where spine (App f a) as = spine f (a:as)
        spine (Lam s e) [] = Lam s (nf e)
        spine (Lam s e) (a:as) = spine (subst s a e) as
        spine f as = app f as
        app f as = foldl App f (map nf as)

betaEq :: Expr -> Expr -> Bool
betaEq e1 e2 = alphaEq (nf e1) (nf e2)

Computing the NF is similar to WHNF, but as we reconstruct expressions we make sure that all subexpression have NF as well. Note that both whnf and nf may fail to terminate because not all expressions have a normal form. The canonical non-terminating example is (λx.x x)(λx.x x) which has one redex, and doing the reduction produces the same term again. But if an expression has a normal form then it is unique (the Church-Rosser theorem). Here are some sample lambda terms (named for convenience):

zero ≡ λs.λz.z one ≡ λs.λz.s z two ≡ λs.λz.s (s z) three ≡ λs.λz.s (s (s z)) plus ≡ λm.λn.λs.λz.m s (n s z)

Or, in Haskell

[z,s,m,n] = map (Var . (:[])) "zsmn"
app2 f x y = App (App f x) y
zero  = Lam "s" $ Lam "z" z
one   = Lam "s" $ Lam "z" $ App s z
two   = Lam "s" $ Lam "z" $ App s $ App s z
three = Lam "s" $ Lam "z" $ App s $ App s $ App s z
plus  = Lam "m" $ Lam "n" $ Lam "s" $ Lam "z" $ app2 m s (app2 n s z)

And now we can check that addition works, betaEq (app2 plus one two) three will produce True.

A detour, simply typed lambda calculus.

To do type checking we need to introduce types. A very simple system is the simply typed lambda calculus. It has one (or more) base type (think of it as Bool or Int if you want an example) and function types. In Haskell terms:

data Type = Base | Arrow Type Type
    deriving (Eq, Read, Show)

Or
t→t B
The expressions themselves will have an explicit type on the bound variable in a lambda expression. So we now have

x e e λx:t.e

For example, (λx:(B→B).λy:B.x)(λz:B.z). The Haskell type for expressions is

data Expr
    = Var Sym
    | App Expr Expr
    | Lam Sym Type Expr
    deriving (Eq, Read, Show)

The only difference is the Type in Lam. All the functions we had for the untyped lambda calculus can be trivially extended to the simply typed one by simply carrying the type along. So finally, time for some type checking. The type checker will take an expression and return the type of the expression. The type checker will also need the types of all free variables in the expression to be able to do this. Otherwise, what type would it assign to, say, the expression x? To represent the types of the free variables we use an environment which is simply a list of variables and their types.

newtype Env = Env [(Sym, Type)] deriving (Show)

initalEnv :: Env
initalEnv = Env []

extend :: Sym -> Type -> Env -> Env
extend s t (Env r) = Env ((s, t) : r)

Type checking can go wrong; there can be type errors. To cater for this the type checker will be written in monadic style where the monad is simply an error (exception) monad. The error messages are strings, and the monad itself is the Either type. So TC is the type checking monad.

type ErrorMsg = String

type TC a = Either ErrorMsg a

We can now write variable lookup.

findVar :: Env -> Sym -> TC Type
findVar (Env r) s =
    case lookup s r of
    Just t -> return t
    Nothing -> throwError $ "Cannot find variable " ++ s

It simply looks up the variable and returns the type. If not found it throws an error with throwError (from Control.Monad.Error). And then the type checker itself.

tCheck :: Env -> Expr -> TC Type
tCheck r (Var s) =
    findVar r s
tCheck r (App f a) = do
    tf <- tCheck r f
    case tf of
     Arrow at rt -> do
        ta <- tCheck r a
        when (ta /= at) $ throwError "Bad function argument type"
        return rt
     _ -> throwError "Non-function in application"
tCheck r (Lam s t e) = do
    let r' = extend s t r
    te <- tCheck r' e
    return $ Arrow t te

For variables, just look up the type for it in the environment. For application, type check the function part and the argument part. The function should have function (arrow) type, and if it does the type of the application is the return type of the function. Finally, for a lambda expression we extend the environment with the bound variable. We then check the body, and the type of the lambda expression is a function type from the argument type to the type of the body. For convenience:

typeCheck :: Expr -> Type
typeCheck e =
    case tCheck initalEnv e of
    Left msg -> error ("Type error:\n" ++ msg)
    Right t -> t

Pretty easy sailing so far. The simply typed lambda calculus is a pain to use. Take something like λx.x in the untyped world. What type should we give it? Well that depends on how we intend to use it. Maybe λx:B.x, maybe λx:(B→B).x, maybe λx:(B→B→B).x... So we can no longer have one identity function; we need one for each type. What a bummer! It's as bad as C. BTW, all (type correct) expression in the simply typed lambda calculus have a normal form (Tait 1967).

Almost going down the wrong track, the polymorphic lambda calculus.

(Don't get me wrong, the polymorphic lambda calculus is a work of marvel.) So how can we fix the problem with one identity function for each type? We can add polymorphism! We can extend the expression language so that we also pass types around; we add type abstraction and type application.
x e e λx:t.e Λα:k.e e[t]
Where Λα:k.e is a type abstraction, i.e., α is a type variable which we can use in type expressions inside e. To supply a type argument we have type application, e[t]. So the types we have would functions, base type, and type variables.
t→t B α
And what is k in the Λα:k.e? Well, now types have gotten so complicated that it is possible to construct types that make no sense, so we need a "type system" for the types. We call them kinds.
k→k *
Defining all this is Haskell would be something like

data Expr
    = Var Sym
    | App Expr Expr
    | Lam Sym Type Expr
    | TLam Sym Kind Expr
    | TApp Expr Type
    deriving (Eq, Read, Show)
data Type
    = Arrow Type Type
    | Base
    | TVar Sym
    deriving (Eq, Read, Show)
data Kind
    = KArrow Type Type
    | Star
    deriving (Eq, Read, Show)

But wait, there's an awful lot of duplication here. The structures on the three levels have a lot of similarities. (Oh, and we don't really need Base anymore now when we have variables.) BTW, this system, called System F_ω, is (a simplified version of) what GHC uses internally to represent Haskell code. It's a beautiful system, really. Oh, the identity function, well it would be Λα:*.λx:α.x. And using it: id[B]a, assuming a has type B.

Simply easy, the lambda cube.

To simplify and (as often happens when you simplify) generalize the expressions above we are going to squish them all into one expression data type. So TLam and Lam will join, as will TApp and App, TVar and Var, KArrow and Arrow. But wait, there's nothing corresponding to Arrow in Expr. We need to add something. We could just add it as it is, but we won't. TADA, enter dependent types. Instead of the boring function type t→u we will use a more exciting one, (x:t)→u. What does it mean? It means that the variable x can occur in u. If it doesn't then it's simply the same as the old fashioned function type. If x does occur it means that type u can depend on the value of the argument (x). In Haskell:

data Expr
    = Var Sym
    | App Expr Expr
    | Lam Sym Type Expr
    | Pi  Sym Type Type
    | Kind Kinds
    deriving (Eq, Read, Show)
type Type = Expr

data Kinds = Star | Box deriving (Eq, Read, Show)

The new arrow type is called Pi. We will also need more than one kind, Star and Box. It's pretty easy to extend the functions from the first part to handle this expression type. There's just a few more places to recurse. Here's the code again. Absolutly nothing subtle about it.

freeVars :: Expr -> [Sym]
freeVars (Var s) = [s]
freeVars (App f a) = freeVars f `union` freeVars a
freeVars (Lam i t e) = freeVars t `union` (freeVars e \\ [i])
freeVars (Pi i k t) = freeVars k `union` (freeVars t \\ [i])
freeVars (Kind _) = []

subst :: Sym -> Expr -> Expr -> Expr
subst v x = sub
  where sub e@(Var i) = if i == v then x else e
        sub (App f a) = App (sub f) (sub a)
        sub (Lam i t e) = abstr Lam i t e
        sub (Pi i t e) = abstr Pi i t e
        sub (Kind k) = Kind k
        fvx = freeVars x
        cloneSym e i = loop i
           where loop i' = if i' `elem` vars then loop (i ++ "'") else i'
                 vars = fvx ++ freeVars e
        abstr con i t e =
            if v == i then
                con i (sub t) e
            else if i `elem` fvx then
                let i' = cloneSym e i
                    e' = substVar i i' e
                in  con i' (sub t) (sub e')
            else
                con i (sub t) (sub e)

To cut down on the code you could actually join the Lam and Pi constructors since they are treated identically in many cases. I've left them separate for clarity. The alphaEq function extends in the natural way to the new type, so does nf, but here it is anyway.

nf :: Expr -> Expr
nf ee = spine ee []
  where spine (App f a) as = spine f (a:as)
        spine (Lam s t e) [] = Lam s (nf t) (nf e)
        spine (Lam s _ e) (a:as) = spine (subst s a e) as
        spine (Pi s k t) as = app (Pi s (nf k) (nf t)) as
        spine f as = app f as
        app f as = foldl App f (map nf as)

So, now for the meaty part, the type checking itself. The handling of the environment is just as before, so we'll just look at the different cases for the type checking.

tCheck :: Env -> Expr -> TC Type
tCheck r (Var s) =
    findVar r s

Just as before.

tCheck r (App f a) = do
    tf <- tCheckRed r f
    case tf of
     Pi x at rt -> do
        ta <- tCheck r a
        when (not (betaEq ta at)) $ throwError $ "Bad function argument type"
        return $ subst x a rt
     _ -> throwError $ "Non-function in application"

This is almost as before, but the arrow type is called Pi now. The key thing here — and this is really where the fact that we are doing dependent types shows up — is the return type. For the simply typed lambda calculus it was just rt, but now rt can contain free occurences of the variable x. Since we are returning rt the x would no longer be in scope, so we substitute the value of the argument for it. This is coolest part of the type checker. You've seen it. That's where it is. Since types can now be arbitrary expression we use betaEq to compare them instead of (==).

tCheck r (Lam s t e) = do
    tCheck r t
    let r' = extend s t r
    te <- tCheck r' e
    let lt = Pi s t te
    tCheck r lt
    return lt

The lambda case is similar to before, but we return a Pi now, so we need to include the variable name. Furthermore, to avoid nonsense like λx:5.e we make sure that the type we want to return actually has a valid kind itself. (The first call to tCheck is to ensure the type we're putting into the environment is valid; I'm sure there's a more elegant way to do this, but I can't remember what it is right now.)

tCheck _ (Kind Star) = return $ Kind Box
tCheck _ (Kind Box) = throwError "Found a Box"

Everything has a type, so what's the type of * (Kind Star), well it's a [] (Kind Box) (excuse the ugly box, I can't find the HTML version of a box). And what's the type of Box? Well, you could keep going, but instead we'll stop right here. The idea is that the source language which we'll write our terms in will not allow the box to be written, so it should never occur.

tCheck r (Pi x a b) = do
    s <- tCheckRed r a
    let r' = extend x a r
    t <- tCheckRed r' b
    when ((s, t) `notElem` allowedKinds) $ throwError "Bad abstraction"
    return t

How do we check the type of the (dependent) function type? Well, we check the type of the thing to the left of the arrow, extend the environment, and then check the thing to the right. So now what should (s, t) be? Well, a and b should be types (or maybe kinds). So their types should be kinds. This leads to the following definition:

allowedKinds :: [(Type, Type)]
allowedKinds = [(Kind Star, Kind Star), (Kind Star, Kind Box), (Kind Box, Kind Star), (Kind Box, Kind Box)]

I.e., we allow (*,*), (*,[]), ([],*), and ([],[]). What does it all mean? Here's the beauty of the lambda cube. By varying what we allow we can change what system we type check.

(*,*) values can depend on values. Just this gives the simply typed λ calculus. ([],[]) types can depend on types. ([],*) values can depend on type. Include all these three and you get F_ω. (*,[]) types can depend on values. Include this one to get dependent types.

With all four combination allowed you get Calculus of Construction (CoC). If you always include (*,*), but make a choice of the other 3 you get 8 choices; these are the corners of the lambda cube. All of these system have been studied. BTW, all the systems in the lambda cube have the property that a well typed expression has a normal form. (Well, the proof of this is so complicated for some of these systems that some people kinda doubt it.)

Examples

Here the syntax s→t means (_:s)→t, where "_" is some new variable not used in t.

Identity
id ≡ λa:*.λx:a.x
Pairs
Pair ≡ λa:*.λb:*.(c:*→((a→b→c)→c)) pair ≡ λa:*.λb:*.λx:a.λy:b.λc:*.λf:(a→b→c).f x y split ≡ λa:*.λb:*.λr:*.λf:(a→b→r).λp:(Pair a b).p r f fst ≡ λa:*.λb:*.λp:(Pair a b).split a b a (λx:a.λy:b.x) p snd ≡ λa:*.λb:*.λp:(Pair a b).split a b b (λx:a.λy:b.y) p

My fingers are numb from all these greek characters, so I'll continue with examples another time. And, of course, a parser and pretty printer.

2007-08-20T20:51:00.000+01:00

Quicksort in Haskell Quicksort is a commonly used example of how succinct Haskell programming can be. It usually looks something likes this:

qsort :: (Ord a Bool) => [a] -> [a]
qsort [] = []
qsort (x:xs) = qsort (filter (<= x) xs) ++ [x] ++ qsort (filter (> x) xs)

The problem with this function is that it's not really Quicksort. Viewed through sufficently blurry glasses (or high abstraction altitude) it's looks like the real Quicksort. What they have in common is overall algorithm: pick a pivot (always the first element), then recursively sort the ones that are smaller, the ones that are bigger, and then stick it all together. But in my opinion the real Quicksort has to be imperative because it relies on destructive update; and it uses a very elegant partitioning algorithm. The partitioning works like this: scan from the left for an element bigger than the pivot, then scan from the right for an element smaller than the pivot, and then swap them. Repeat this until the array has been partitioned. Haskell has a variety of array types with destructive updates (in different monads), so it's perfectly possible to write the imperative Quicksort in Haskell. To make the code look reasonably nice I'm going to use my C-like DSEL to write the code. Here it is:

qsortM :: forall i a m r arr .
         (MonadRef m r, Num i, Ix i, MArray arr a m, Ord a Bool) =>
         arr i a -> m (arr i a)
qsortM ma = runE $ do
    (lb, ub) <- embed $ getBounds ma
    let mlb = pure0 lb
        mub = pure0 ub
    a <- liftArray ma

    let qsort' l r =
            if1 (r > l) $ do
                i <- auto l
                j <- auto (r+1)
                let v = a[l] :: E m a
                    iLTj = i < (j :: E m i)
                while iLTj $ do
                    while ((i += 1) < mub && a[i] < v)
                        skip
                    while (a[j -= 1] > v)
                        skip
                    if1 iLTj $ do
                        a[i] =:= a[j]
                a[l] =:= a[j]
                qsort' l (j-1)
                qsort' i r

    qsort' mlb mub
    return ma

So the type is arr i a -> m (arr i a), i.e., qsortM takes an array indexed with i and elements of type a. It returns the sorted array, but the sorting takes places in some monad m. And then there are all kinds of constraints on the type variables. The MonadRef m r says that the monad has to have references so we can have some variables. The array index has to be in the Ix; that's part of the general array constraints. It also have to be numeric so we can add and subtract indicies. The array type has to fulfill MArray arr a m which means that arr is an array of a and updatable in monad m. Finally, the elements have to be ordered. I'm not using the normal Ord class, but instead an overloaded Ord where the return type is overloaded too. A few comments on the code. The liftArray function lifts a regular array into one that can be indexed with [i]. The =:= operator swaps two variables. The skip function does nothing. In traditional C style we do all the side effects while computing the condition. There are some type signatures in the code that are annoying, but that I have not found a way around yet. But otherwise, the code proceeds like most any imperative Quicksort. The i and j variables scan the array from left and right to locate two elements that need swapping. We then swap them, and continue scanning until the indicies cross. After the partitioning we swap the pivot into place and sort the two parts recursively. So, this function is polymorphic in the monad. But there is one monad that I think is extra interesting, namely ST. With this monad you can do references, updatable arrays, etc., and finally you can seal it all of with runST. The resulting type is pure and shows no signs of what happened inside. This is a really amazing feat, in my opinion. The type checker performs the proof that nothing about the "dirty" innards leaks out. So instead of some informal reasoning that a function with an impure inside can be pure on the outside you have a machine checked proof. Of course, there's a meta proof that this is all correct, but John Launchbury and Simon Peyton Jones have already done that once, and now we just need the type checker. Here's the code:

qsortA :: forall i a . (Num i, Ix i, Ord a Bool) => Array i a -> Array i a
qsortA a = runSTArray sa
  where sa :: ST s (STArray s i a)
        sa = thaw a >>= qsortM

We're using runSTArray instead of runST, because it provides an efficient way to turn a mutable array on the inside into an immutable array on the outside. The thaw function turns an immutable array into a mutable one, but it has to make a copy to be safe, since we don't want to mutate the original. Finally, if we want to sort lists we can always convert back and forth.

asList :: (Array Int a -> Array Int a) -> ([a] -> [a])
asList f xs = elems . f . listArray (1, length xs) $ xs 

qsort :: (Prelude.Ord a) => [a] -> [a]
qsort = asList qsortA

The final qsort has the normal type signature (I switched to the normal Ord again).

2007-08-16T11:17:00.000+01:00

What about arrays? After doing my little C-like DSEL in Haskell I started wondering if you could do arrays that look like C arrays too. It turns out you can, but I can't figure out a way to make it safe from certain runtime errors. Array indexing in C is written a[i], and this is legal syntax in Haskell too, so we're off to a good start. But to make it work we need a to be a function that takes a list and returns something. What should it return? It must be something that is both an l-value and an r-value. Just as I had a function auto to create a local variable, I'll have a function arr to make a local array. It will take the size of the array and then return a function that takes the index. For simplicity I'll make the index be an Int.

arr :: forall a . [E Int] -> E (forall v . [E Int] -> E' v a)

So now this type checks

atest = do {
    a <- arr[2];
    a[1] =: a[0] + 1;
  }

Now we just need the body of arr.

arr [s] = do
    s' <- runE s
    a <- newArray (0, s' - 1) undefined :: IO (IOArray Int a)
    let ix [i] = runE i
    return (\ is -> V (ix is >>= readArray a)
                      (\ x -> ix is >>= \ i -> writeArray a i x))

The arr function takes a list with one element, the size, and allocates an array (indexed from 0) with this size. It then returns a function that expects a singleton list with an index and returns the V constructor which I used for variables. For multidimensional arrays we can extend the arr function.

arr ss = do
    ss' <- mapM runE ss
    let sz = product ss'
        ix is = do
                    is' <- mapM runE is
                    when (length is' /= length ss') $ error "wrong number of indicies"
                    return $ foldr (\ (i, s) r -> r * s + i) 0 (zip is' ss')
    a <- newArray (0, product ss' - 1) undefined :: IO (IOArray Int a)
    return (\ is -> V (ix is >>= readArray a) (\ x -> ix is >>= \ i -> writeArray a i x))

The problem with both these definitions is that the number of indicies is checked dynamically rather than statically. To do it statically we'd have to be able to overload the syntax for list literals to use a type that keeps track of the number of elements.

2007-08-13T23:47:00.000+01:00

Programming in C, ummm, Haskell Here's a small Haskell program for computing factorial.

fac n = do {
    a <- auto 1;
    i <- auto n;
    while (i >. 0) $ do {
        a *= i;
        i -= 1;
    };
    a;
  }

It even runs, runE (fac 10) produces 3628800. Let's compare that to the C program doing the same thing.

int
fac(int n) {
    auto a = 1;
    auto i = n;
    while (i > 0) {
        a *= i;
        i -= 1;
    }
    return a;
}

They look rather similar, don't they? (BTW, I decided to use the heavily underused C keyword auto. If you don't know C, don't worry, auto doesn't mean anything. A local declaration auto int i; is the same as int i; in C.) I often hear people complain that C-like programming in Haskell is ugly and verbose, but I don't think that has to be the case. What people complain about is when they write code like this:

fac' n = do
    a <- newIORef 1
    i <- newIORef n
    whileM (do i' <- readIORef i; return $ i' > 0) $ do
        i' <- readIORef i
        a' <- readIORef a
        writeIORef a (a' * i')
        writeIORef i (i' + 1)
    a' <- readIORef a
    return a'

And I agree, it's ugly, it's verbose, but it's not necessary. One of the reasons it's verbose are the uses of readIORef and writeIORef. In C you just refer to the variable and you'll get an r-value when you need it, and you'll get an l-value when you need it. And you can do the same in Haskell. Let's look at a tiny version of the problem. We want the following operations:

type V a ...                     -- variables of type a
type E a ...                     -- expressions of type a
(=:) :: V a -> E a -> IO ()      -- assignment
plus :: E Int -> E Int -> E Int  -- addition
one :: E Int                     -- constant 1

Here's a tiny sample:

  x <- auto one
  x =: x `plus` x

-- the following should be a type error
  (x `plus` x) =: one

How can we make this happen? We want x to be able to be the left hand side of an assignment as well as an expression on the right hand side. But the assignment operator has different types on the left and right. So we could possibly make V the same as E. But this would be bad, because we want the left hand side to be a variable and now we'd have no way of knowing, which means that it would no longer be a type error to assign to a non-variable. So, we need a little bit of type trickery. First, we'll make V and E "subtypes" of the same type, so they have something in common.

data E' v a = ...
data LValue
data RValue
type V a = E' LValue a
type E a = E' RValue a

The type E' has an additional type argument that determines if the value is an l-value or an r-value. So what about auto, what type should it have? It needs to return something that is both an l-value and an r-value. So we'll make it polymorphic! First attempt:

auto :: E a -> IO (E' v a)

But this won't work. Saying "x <- auto 2; ..." is the same as "auto 2 >>= \ x -> ...". And when a variable is lambda bound it is not polymorphic. What does this mean? It means that we can use x as either an l-value or as an r-value inside ..., but not both; the v type variable can only have a single type. Are we stuck? Well, we would have been without higher ranked polymorphism. But here's a type that works:

auto :: E a -> IO (forall v . E' v a)

So now x will have type forall v . E' v a which is indeed polymorphic in v just as we wanted. OK, so now we have some type signatures that work (using stub implementations it's easy to check that the types work out). So now we need to implement the types and the operations. Let's start with r-values. To make a constructor that only can construct r-values we'll use a GADT. The E type will simply embed an IO value. So when we see the type E' RValue a it's really isomorphic to IO a. To extract the IO value we have the function runE.

data E' v a where
    E :: IO a -> E' RValue a

runE :: E' v a -> IO a
runE (E t) = t

So now we can do plus and one; they are totally straightforward except that we need to unwrap and wrap the E constructor.

plus :: E Int -> E Int -> E Int
plus x y = E $ do
    x' <- runE x
    y' <- runE y
    return $ x' + y'

one :: E Int
one = E $ return 1

And then the hard part, variables. To represent a variable we'll use two fields, one that reads the variable and one that assigns the variable. We need to extend the runE function to use the variable read field to get the value. The auto function allocates a new variable and then packages up the two fields.

data E' v a where
    E :: IO a -> E' RValue a
    V :: IO a -> (a -> IO ()) -> E' v a

runE :: E' v a -> IO a
runE (E t) = t
runE (V t _) = t

auto :: E a -> IO (forall v . E' v a)
auto x = do
    x' <- runE x
    r  <- newIORef x'
    return (V (readIORef r) (writeIORef r))

We're about done, just assignment left. It's easy, just use the assignment field from the V constructor.

(=:) :: V a -> E a -> IO ()
(V _ asg) =: e = do
    e' <- runE e
    asg e'

Hmmm, but there's a missing case in that definition. What about the constructor E? Can't we get a runtime failure? No. Look at the type of the first argument, it's V a, i.e., E' LValue a. The E constructor can only construct values of type E' RValue a, so we just can't get a match failure. Now our little example above compiles, and trying the bad expression (x `plus` x) =: one gives this error:

    Couldn't match expected type `LValue'
    against inferred type `RValue'
      Expected type: V a
      Inferred type: E Int

Once we've made the E data type it's totally routine to make the factorial function I first showed work. We just need to create an instance for Num (E a) and a few more functions. Is this the end of the story? Is everything rosy? Is C programming in Haskell as smooth as that? Well, there are some flies in the ointment. While it's true that programming with the E type is easy, it's a little cumbersome if we want to use pure functions. Pute functions will not have the E type and will have to be lifted into this new brave world. Likewise, IO types will have to be lifted. So while we have reduced some verbosity, it has popped up again in a different place. Which is better? I don't know, but I thought I'd share the l-value/r-value trick. Oh, and (of course), there's nothing special about the IO monad. I just used it as an example; any monad with references will work. You can parametrize E' over the underlying monad.

2007-07-09T10:07:00.001+01:00

Experiencing magic What happens when you let a road planner be on crack at work? You get the Magic Roundabout! I've seen and driven it now. It's clear to me why the locals sometimes refers to this as the Tragic Roundabout, because there's no way an unsuspecting stranger can navigate through this unless totally sober, super alert, there's full daylight, light traffic, ... I recommend a Sunday afternoon. Have a look. BTW, this is the the UK, so if you drive on the right, flip the image. :)

2007-07-01T03:16:00.000+01:00

Massive overload In my last post I had a little DSL embedded in Haskell. Defining the Fibonacci function looked like this:

mdo
    fib <- fun $ \ n -> cond (n .< 2) 1 (fib(n-1) + fib(n-2))
    return fib

Compare that with ordinary Haskell:

fib n = if n < 2 then 1 else fib(n-1) + fib(n-2)

Why do they look different? Could I make my DSL look exactly like Haskell? First, lets look at the parts that look the same, like fib(n-1). In Haskell, with n of type Int, this just an expression and it has some value (depending on n). In the DSL, with n of type M Int the value of this expression is in fact an abstract syntax tree, namely (something like) M (App (Func (FuncNo 0)) [App (FPrimOp I_Sub) [Arg 0, Con (ILit 1)]]). The reason we can make the same expression mean different things depending on the type is that numerical operators like - are overloaded, and so they do different things for Int and M Int. So what about <, can we overload that too? The < operator is already overloaded and has type (Ord a) => a -> a -> Bool. But the .< operator in the DSL returns M Bool, so there's no way we can use <; the return type is fixed. What to do? Well, there's nothing sacred with the operator <, it's just another name, but defined in the Prelude. We can make our own if we hide the Prelude. So let's do that, but carefully. We want to be able to still use < as before, and for our M Bool return type. So we need to overload the return type as well as the argument type. Overloaded Booleans Before we do that, let's define a class that is analogous to Num, but for Boolean values. We'll need this class soon. The Bool type has a few important functions and values, let's put those in the class.

 b -> b
    not :: b -> b

instance Boolean Bool where
    false = False
    true = True
    (&&) = (P.&&)
    (||) = (P.||)
    not = P.not

First we import the Prelude explicitly, this allows us to hide some things we want to redefine. The we define the Boolean class, which gives new meaning to some prelude functions. Finally, we make an instance saying that the new (&&), (||), and not behave as the old ones for good oldfashioned Bool. After this definition we can use the Boolean operators just as before, they are just overloaded now, but for Bool they resolve to their old meaning. Overloaded comparison So back to comparison, Eq first.

class (Boolean b) => Eq a b | a -> b where
    (==), (/=) :: a -> a -> b
    x /= y  =  not (x == y)

The new Eq (we need to hide the Prelude one) is now more general; the return type is overloaded as well. Too avoid ambiguities we have a functional dependency. The type of the values we compare will determine the type of the Boolean result. (This isn't necessary, but without this we might need many type annotations.) Nothing is a member of this class, so we need to add the types we need. We really should all types that have equality, but let's do two samples.

instance Eq Int Bool where
    (==) = (P.==)
    (/=) = (P./=)
instance Eq Double Bool where
    (==) = (P.==)
    (/=) = (P./=)

And since the whole point of this was to allow special equality for our M values, let's look at that too.

instance Eq (M Int) (M Bool) where
    (==) = binOp I_EQ
    (/=) = binOp I_NE

A new (simplified) Ord looks similar, both the class and the instances.

class (Eq a b) => Ord a b | a -> b where
    (<), (<=), (>), (>=) :: a -> a -> b

Overloaded conditionals So now for the next challange. The Haskell Fibonacci function uses if, and the DSL uses a special function cond. Well, now we are stuck. These is no way at all to overload the special if syntax in Haskell; it's wired in. This could be viewed as a design mistake in Haskell, but so it is. So could we use the cond function in the Haskell version too? Well, not the original cond, it has type M Bool -> M Int -> M Int -> M Int. We need to overload that too.

class (Boolean b) => Cond a b | a -> b where
    cond :: b -> a -> a -> a

Again, we add a functional dependency. And again, it will save some ambiguities. The obvious instances:

instance Cond Int Bool where
    cond x y z = if x then y else z

instance Cond (M Int) (M Bool) where
    cond = terOp I_Cond
instance Cond (M Bool) (M Bool) where
    cond = terOp I_Cond
terOp op (M c) (M t) (M e) = M $ App (FPrimOp op) [c, t, e]

Oh, and the instance for M Bool that I skipped before.

instance Boolean (M Bool) where
    false  = M $ Con $ LInt 0
    true   = M $ Con $ LInt 1
    x && y = cond x y false
    x || y = cond x true y
    not x  = cond x false true

After this effort we can write the body of the Fibonacci function the same way for both type Int and type M Int, namely

cond (n < 2) 1 (fib(n-1) + fib(n-2))

And what about all these new classes and instances. how do we package them to make them easy to use?

module MyPrelude(module Prelude, Boolean(..), Eq(..), Ord(..), Cond(..)) where
import qualified Prelude as P
import Prelude hiding (Eq(..), Ord(..), (&&), (||), not)

...

We hide the stuff we want to override, and re-export the rest of the Prelude. Any module that wants to use this new Prelude now has to say:

module M where
import Prelude()
import MyPrelude
...

After this all regular Haskell code works as before, but we can also take advantage of the new overloadings. Enough for today.

2007-06-29T17:53:00.001+01:00

Generating more code with Harpy After a slight detour into expression representation I'm back to generating code again. To recap, here's how I'm going to do my embedded expressions in Haskell. For now I'll stick to ugly looking boolean operators. Fixing that is a whole other story.

module Exp(Exp(..), Lit(..), Fun(..), PrimOp(..)) where

data Exp
    = Con Lit
    | Arg Int
    | App Fun [Exp]
    deriving (Eq, Ord, Show)

data Lit
    = LInt  Int
    deriving (Eq, Ord, Show)

data Fun
    = FPrimOp PrimOp
    deriving (Eq, Ord, Show)

data PrimOp
    = I_Add | I_Sub | I_Mul | I_Quot | I_Rem | I_Neg
    | I_EQ | I_NE | I_LT | I_LE | I_GT | I_GE
    | I_Cond
    deriving (Eq, Ord, Show)
-----------------------------------------------------
module MExp(M(..)) where
import Exp

newtype M a = M { unM :: Exp }
    deriving (Show)
-----------------------------------------------------
module M(M, cond,
         (.==), (./=), (.<), (.<=), (.>), (.>=),
         false, true, (.&&), (.||)
        ) where
import Exp
import MExp

instance Eq (M a)
instance Ord (M a)

instance Num (M Int) where
    x + y  =  binOp I_Add x y
    x - y  =  binOp I_Sub x y
    x * y  =  binOp I_Mul x y
    negate x =  unOp I_Neg x
    fromInteger i = M $ Con $ LInt $ fromInteger i

instance Enum (M Int)
instance Real (M Int) where

instance Integral (M Int) where
    quot x y = binOp I_Quot x y
    rem x y = binOp I_Rem x y

binOp op (M x) (M y) = M $ App (FPrimOp op) [x, y]
unOp op (M x) = M $ App (FPrimOp op) [x]

--------

infix 4 ./=, .==, .<, .<=, .>, .>=
(.==), (./=), (.<), (.<=), (.>), (.>=) :: M Int -> M Int -> M Bool
(.==) = binOp I_EQ
(./=) = binOp I_NE
(.<)  = binOp I_LT
(.<=) = binOp I_LE
(.>)  = binOp I_GT
(.>=) = binOp I_GE

cond :: M Bool -> M Int -> M Int -> M Int
cond (M c) (M t) (M e) = M $ App (FPrimOp I_Cond) [c, t, e]

condB :: M Bool -> M Bool -> M Bool -> M Bool
condB (M c) (M t) (M e) = M $ App (FPrimOp I_Cond) [c, t, e]

false, true :: M Bool
false = M $ Con $ LInt 0
true  = M $ Con $ LInt 1

infixr 3 .&&
infixr 2 .||
(.&&), (.||) :: M Bool -> M Bool -> M Bool
x .&& y = condB x y false
x .|| y = condB x true y

So the Exp type is the internal representation of expressions. It's just constants, variables, and some primitive operations. The M (for machine) type is the phantom type that the DSL user will see. The module M contains all the user visible functions. Note that it handles booleans by represented the C way, by an int that is 0 or 1. There's a conditional function, cond, and a corresponding primitive that will serve as our if. So now we need to generate code for these. This is all very similar to what I did in a previous posting, I have just refactored some of the code generation. The invariant (that I arbitrarily decided on) is that each block of code that implements a primitive operation will may clobber EAX but must leave all other registers intact. This is mostly just a chunk of rather boring code. Division is causing problems as usual, since the IA32 instructions don't allow you to specify a destination register. One point worth a look is the code generation for cond. It has to test a boolean and then select one of two code blocks, but that is just what you'd expect.

module CodeGen(cgExp, CGen, cgPrologue, cgEpilogue, compileFun) where
import Prelude hiding(and, or)
import Control.Monad.Trans
import Data.Maybe(fromJust)
import Foreign
import Harpy.CodeGenMonad
import Harpy.X86Assembler

import Exp
import MExp

type StackDepth = Int
type Env = ()
type CGen a = CodeGen Env StackDepth a

addDepth :: StackDepth -> CGen ()
addDepth i = do
    d <- getState
    setState (d+i)

cgExp :: Exp -> CGen ()
cgExp (Con l) = cgLit l
cgExp (Arg a) = cgArg a
cgExp (App (FPrimOp I_Cond) [c, t, e]) = cgCond c t e
cgExp (App f es) = do mapM_ cgExp es; cgFun f (length es)

cgFun (FPrimOp p) _ = cgPrimOp p

cgPrimOp I_Add = twoOp add
cgPrimOp I_Sub = twoOp sub
cgPrimOp I_Mul = twoOp (imul InPlace)
cgPrimOp I_Quot = twoOp (cgQuotRem eax)
cgPrimOp I_Rem = twoOp (cgQuotRem edx)
cgPrimOp I_Neg = oneOp neg
cgPrimOp I_EQ = cmpOp sete
cgPrimOp I_NE = cmpOp setnz
cgPrimOp I_LT = cmpOp setl
cgPrimOp I_LE = cmpOp setle
cgPrimOp I_GT = cmpOp setg
cgPrimOp I_GE = cmpOp setge

cgCond c t e = do
    cgExp c
    popReg eax
    test eax eax
    l1 <- newLabel
    l2 <- newLabel
    jz l1
    cgExp t
    addDepth (-1)               -- pretend last cgExp didn't push anything
    jmp l2
    l1 @@ cgExp e
    l2 @@ return ()

cmpOp op = twoOp $ \ r1 r2 -> do
    cmp r1 r2
    op (reg32ToReg8 r1)
    and r1 (1 :: Word32)

cgLit (LInt i) = do
    mov eax (fromIntegral i :: Word32)
    pushReg eax

cgArg :: Int -> CGen ()
cgArg n = do
    d <- getState
    let o = 4 * (d + n)
    mov eax (Disp (fromIntegral o), esp)
    pushReg eax

cgQuotRem res r1 r2 = do
    push edx                    -- save temp reg
    push r2                     -- push second operand (in case it's edx)
    xmov  eax r1                -- first operand must be in eax
    mov  edx eax
    sar  edx (31 :: Word8)      -- sign extend
    idiv (Disp 0, esp)          -- divide by second operand (now on stack)
    add  esp (4 :: Word32)      -- remove second operand
    xmov  r1 res                -- put result in r1
    pop  edx                    -- restore edx

-- Do a register to register move, but do generate no-ops
xmov r1 r2 = if r1 == r2 then return () else mov r1 r2

reg32ToReg8 (Reg32 r) = Reg8 r

--------

twoOp op = do
    pop ebx
    pop ecx
    op ecx ebx
    push ecx
    addDepth (-1)

oneOp op = do
    pop ebx
    op ebx
    push ebx

pushReg r = do
    push r
    addDepth 1

popReg r = do
    pop r
    addDepth (-1)

--------

compileFun :: (M Int -> M Int) -> CGen ()
compileFun f = do
    ensureBufferSize 1000
    setState startDepth
    cgPrologue
    cgExp $ unM $ f $ M $ Arg 0
    cgEpilogue

----

savedRegs = [ebx, ecx]

startDepth = length savedRegs + 1

-- Push all register we'll be using, except eax.
cgPrologue :: CGen ()
cgPrologue = do
    mapM_ push savedRegs

-- Pop return value to eax, restore regs, and return
cgEpilogue :: CGen ()
cgEpilogue = do
    popReg  eax
    mapM_ pop (reverse savedRegs)
    ret

Given all this, we can generate code for a function and call it from Haskell code. But that's rather boring since we have only primitive functions and no way to call DSL functions from within the DSL. So we have no recursion. And as we all know, recursion is where the fun starts. So we need to add expressions and code generation for calling functions. First we extend the Exp module.

data Fun
    = FPrimOp PrimOp
    | Func FuncNo

newtype FuncNo = FuncNo Int

The idea being that our DSL functions will be represented by a FuncNo which is an Int. Next, we will extend the code generation. While generating code we will keep an environment that maps function numbers to the corresponding labels to call.

type Env = [(FuncNo, Label)]
...
cgFun (Func f) n = do
    e <- getEnv
    call $ fromJust $ lookup f e
    add esp (fromIntegral (4 * n) :: Word32)
    addDepth (-n)
    pushReg eax

To call a function we look up its label and generate a call instruction to that label. Then we pop off the arguments that we pushed before the call, and finally we push EAX (where the return value is). Simple, huh? But where did the environment come from? Well, it's time to be a little inventive and introduce a monad. Functions are represented by unique integers in the expression type, so we'll need a state monad to generate them.

module Gen where
import Control.Monad.State
import Harpy.CodeGenMonad
import Exp
import MExp
import CodeGen

data GenState = GenState { funcs :: [(FuncNo, CGen ())], nextFunc :: Int }

startGenState :: GenState
startGenState = GenState { funcs = [], nextFunc = 0 }

type G a = State GenState a

fun :: MInt2MInt -> G MInt2MInt
fun f = do
    s <- get
    let n = nextFunc s
        fno = FuncNo n
    put GenState{ funcs = (fno, compileFun f) : funcs s,
                  nextFunc = n + 1 }
    return $ \ x -> M $ App (Func fno) [unM x]

runG :: G MInt2MInt -> CGen ()
runG g = do
    let (ret, s) = runState g startGenState
        funs = reverse $ funcs s
    funMap <- mapM (\ (fno, _) -> do l <- newLabel; return (fno, l)) funs
    withEnv funMap $ do
        compileFun ret
        zipWithM_ (\ (f, l) (_, g) -> l @@ g) funMap funs

This code actually has some interesting points. First, the monad, G, keeps a list of defined functions. The list has the function number and the CodeGen block that will generate code for it. Second, the fun function is the one that creates a new function number. It also calls compileFun to get a CodeGen block that will generate code for the function. Note that no code is generated at this point. The G monad is just a simple state monad, not the IO based codeGen monad. Also note how fun returns an expression that is of the same type as the argument, but it now uses an App to call the function. Finally, the runG function generates all the code. It uses runState to obtain the list of defined function and the result to return. Now code generation starts. For each function we generate a new label. Pairing up the function number and the label forms our environment, funMap. With this environment installed we start code generation for the functions, first the return value, then all the defined function; each with its label attached. Well, that's it. Let's wrap it all up.

module Compile(compileIO, disasmIO, module M, module Gen) where
import Convert
import Gen
import M

type MInt = M Int

compileIO :: G (MInt -> MInt) -> IO (Int -> Int)
compileIO f = fmap flex $ compileIOW32 f
  where flex g = fromIntegral . g . fromIntegral

compileIOW32 :: G (MInt -> MInt) -> IO (Word32 -> Word32)
compileIOW32 = conv_Word32ToWord32 [] undefined . runG

disasmIO :: G (MInt -> MInt) -> IO String
disasmIO = disasm [] undefined . runG

And, of course, a test. Here we hit a snag. We want recursion, but monadic bindings (do) are not recursive. Luckily, ghc does implement recursive do called mdo (why that name, I have no idea) for any monad that is in the class MonadFix. And the state monad is, so we are in luck.

{-# OPTIONS_GHC -fglasgow-exts #-}
module Main where
import Compile

main = do
    let g = mdo
                fib <- fun $ \ n -> cond (n .< 2) 1 (fib(n-1) + fib(n-2))
                return $ fib
    test <- compileIO g
    print (test 40)

It even runs and produces the right answer, 165580141. Running time for this example is about 4.0s on my machine. Running fib compiled with 'ghc -O' takes about 3.4s, so we're in the same ballpark. Oh, and if anyone wonders what the code looks like for this example, here it is. This is really, really bad code.

003d6cd0  53                            push   ebx
003d6cd1  51                            push   ecx
003d6cd2  8b 44 24 0c                   mov    eax,dword ptr [esp+12]
003d6cd6  50                            push   eax
003d6cd7  e8 08 00 00 00                call   [3d6ce4H]
003d6cdc  83 c4 04                      add    esp,4H
003d6cdf  50                            push   eax
003d6ce0  58                            pop    eax
003d6ce1  59                            pop    ecx
003d6ce2  5b                            pop    ebx
003d6ce3  c3                            ret
003d6ce4  53                            push   ebx
003d6ce5  51                            push   ecx
003d6ce6  8b 44 24 0c                   mov    eax,dword ptr [esp+12]
003d6cea  50                            push   eax
003d6ceb  b8 02 00 00 00                mov    eax,2H
003d6cf0  50                            push   eax
003d6cf1  5b                            pop    ebx
003d6cf2  59                            pop    ecx
003d6cf3  3b cb                         cmp    ecx,ebx
003d6cf5  0f 9c c1                      setl   cl
003d6cf8  83 e1 01                      and    ecx,1H
003d6cfb  51                            push   ecx
003d6cfc  58                            pop    eax
003d6cfd  85 c0                         test   eax,eax
003d6cff  0f 84 0b 00 00 00             je     [3d6d10H]
003d6d05  b8 01 00 00 00                mov    eax,1H
003d6d0a  50                            push   eax
003d6d0b  e9 37 00 00 00                jmp    [3d6d47H]
003d6d10  8b 44 24 0c                   mov    eax,dword ptr [esp+12]
003d6d14  50                            push   eax
003d6d15  b8 01 00 00 00                mov    eax,1H
003d6d1a  50                            push   eax
003d6d1b  5b                            pop    ebx
003d6d1c  59                            pop    ecx
003d6d1d  2b cb                         sub    ecx,ebx
003d6d1f  51                            push   ecx
003d6d20  e8 bf ff ff ff                call   [3d6ce4H]
003d6d25  83 c4 04                      add    esp,4H
003d6d28  50                            push   eax
003d6d29  8b 44 24 10                   mov    eax,dword ptr [esp+16]
003d6d2d  50                            push   eax
003d6d2e  b8 02 00 00 00                mov    eax,2H
003d6d33  50                            push   eax
003d6d34  5b                            pop    ebx
003d6d35  59                            pop    ecx
003d6d36  2b cb                         sub    ecx,ebx
003d6d38  51                            push   ecx
003d6d39  e8 a6 ff ff ff                call   [3d6ce4H]
003d6d3e  83 c4 04                      add    esp,4H
003d6d41  50                            push   eax
003d6d42  5b                            pop    ebx
003d6d43  59                            pop    ecx
003d6d44  03 cb                         add    ecx,ebx
003d6d46  51                            push   ecx
003d6d47  58                            pop    eax
003d6d48  59                            pop    ecx
003d6d49  5b                            pop    ebx
003d6d4a  c3                            ret

2007-06-28T17:11:00.001+01:00

Representing DSL expressions in Haskell When you want to embed a DSL in Haskell you will almost certainly be faced with having some kind of expressions in your DSL. If you are lucky you can get away with using the Haskell types as the DSL type (i.e., use the meta types as object types). E.g., if you DSL needs integer expressions you can use the Haskell type Integer. But sometimes you're not that lucky, so you need to represent the DSL expressions as a Haskell data type. What are the options for representing expressions in Haskell? A normal data type Let's start with the most obvious one, an ordinary data type with constructors for the different kinds of expressions. As an example I'll have a simple expression language that has constants (of type Int), addition, less-or-equal, and conditional.

data Exp
    = ConI Int
    | Add Exp Exp
    | LE Exp Exp
    | Cond Exp Exp Exp
    deriving (Show)

Just to exemplify, here's an evaluator for these expressions. We need a type of values that the evaluator can return as well.

data Val = ValI Int | ValB Bool
    deriving (Show)

eval :: Exp -> Val
eval (ConI i) = ValI i
eval (Add x y) = case (eval x, eval y) of
                 (ValI vx, ValI vy) -> ValI (vx + vy)
                 _ -> error "Bad arguments to Add"
eval (LE x y) = case (eval x, eval y) of
                (ValI vx, ValI vy) -> ValB (vx <= vy)
                _ -> error "Bad arguments to LE"
eval (Cond x y z) = case (eval x) of
                    ValB b -> if b then eval y else eval z
                    _ -> error "Bad arguments to Cond"

And a test run:

main = print $ eval e
  where e = Cond (LE (ConI 3) (ConI 5)) (Add (ConI 1) (ConI 2)) (ConI 0)

prints ValI 3 This is OK, but look at all the error cases in the evaluator. Remember that these are expressions in a DSL, so we are going to have expression fragments in our Haskell code. If we make a mistake, like Cond (ConI 1) (ConI 1) (ConI 1) this is not going to be caught by the Haskell type checker since everything is of type Exp. Instead it is going to cause some error when the program is run. Being advocates of static type checking (why else use Haskell?), this is rather disgusting. GADTs So let's try a new and shiny feature in GHC, namely GADTs. (GADTs are new in ghc, but the idea is old; it's been well known in constructive type theory for ages, Kent Petersson and I suggested it as an addition to Haskell almost 15 years ago.) With a GADT you can specify more precise types for the constructors.

data Exp a where
    ConI :: Int -> Exp Int
    Add  :: Exp Int -> Exp Int -> Exp Int
    LE   :: Exp Int -> Exp Int -> Exp Bool
    Cond :: Exp Bool -> Exp a -> Exp a -> Exp a

These are the types we want. It's now impossible to construct ill typed value of type Exp t; it will be caught by the Haskell type checker. The evaluator looks very neat and natural with GADTs

eval :: Exp a -> a
eval (ConI i) = i
eval (Add x y) = eval x + eval y
eval (LE x y) = eval x <= eval y
eval (Cond x y z) = if eval x then eval y else eval z

GADTs is really the way to go, but it has the disadvantage of not being standard Haskell. It can also get somewhat cumbersome when we have variables; the evaluator now needs a typed environment. It's certainly doable, but starting to look less nice. Phantom types Let's explore a third way that in some sense combines the previous two. The idea is to have an untyped representation, like Exp, but to use an abstract data type on top of it that only allows constructing well typed expressions. We start with the original type, but rename it to Exp'

data Exp'
    = ConI Int
    | Add Exp' Exp'
    | LE Exp' Exp'
    | Cond Exp' Exp' Exp'
    deriving (Show)

On top of this we provide the type the user will see.

newtype Exp a = E Exp'
    deriving (Show)

conI :: Int -> Exp Int
conI = E . ConI

add :: Exp Int -> Exp Int -> Exp Int
add (E x) (E y) = E $ Add x y

le :: Exp Int -> Exp Int -> Exp Bool
le (E x) (E y) = E $ LE x y

cond :: Exp Bool -> Exp a -> Exp a -> Exp a
cond (E x) (E y) (E z) = E $ Cond x y z

The functions conI, add, le, and cond are the only ones the DSL user will see. Note that they have the same types as the constructors in the GADT. To ensure the DSL user can only use these we need to put it all in a module and export the right things.

module Exp(Exp, conI, add, le, cond) where ...

We can write an evaluator again (inside the Exp module).

eval :: Exp a -> Val
eval (E x) = eval' x

eval' :: Exp' -> Val
eval' (ConI i) = ValI i
eval' (Add x y) = case (eval' x, eval' y) of
                 (ValI vx, ValI vy) -> ValI (vx + vy)
                 _ -> error "Bad arguments to Add"
eval' (LE x y) = case (eval' x, eval' y) of
                 (ValI vx, ValI vy) -> ValB (vx <= vy)
                 _ -> error "Bad arguments to LE"
eval' (Cond x y z) = case (eval' x) of
                    ValB b -> if b then eval' y else eval' z
                    _ -> error "Bad arguments to Cond"

Not as elegant as the GADT evaluator, but at least the error cases are impossible when expressions are constructed using the exported interface. And the test would now look like

main = print $ eval e
  where e = cond (le (conI 3) (conI 5)) (add (conI 1) (conI 2)) (conI 0)

Phantom types gets its name from the fact that the type variable in the definition of the Exp type doesn't occur anywhere in the constructors; so it's a phantom. Now and then people suggest implementing a dynamically typed version of Haskell. Phantom types is one of the things that makes such an implementation difficult. There are no values that carry the type information at runtime, the types only exist at compile time.

2007-06-26T18:17:00.000+01:00

Disassembly The Harpy package also contains a disassembler, so let's put it to work.

disasm :: e -> s -> CodeGen e s a -> IO String
disasm e s cg = do
    let cg' = do cg; getCodeBufferList
    (_, r) <- runCodeGen cg' e s
    case r of
        Left msg -> error (show msg)
        Right bs -> fmap concat $ mapM disAsm bs
  where disAsm (ptr, n) = do
            r <- disassembleBlock ptr n
            case r of
                Left msg -> error $ show msg
                Right insns -> return $ unlines $ map showIntel insns

Using the same example as in the last post we get.

main = do
    let fun x = (x+1) * x `quot` 2
    str <- disasmIO fun
    putStr str

And here is the sad truth:

003d6d20  53                            push   ebx
003d6d21  51                            push   ecx
003d6d22  52                            push   edx
003d6d23  8b 44 24 10                   mov    eax,dword ptr [esp+16]
003d6d27  50                            push   eax
003d6d28  6a 01                         push   1H
003d6d2a  5b                            pop    ebx
003d6d2b  58                            pop    eax
003d6d2c  03 c3                         add    eax,ebx
003d6d2e  50                            push   eax
003d6d2f  8b 44 24 14                   mov    eax,dword ptr [esp+20]
003d6d33  50                            push   eax
003d6d34  5b                            pop    ebx
003d6d35  58                            pop    eax
003d6d36  0f af c3                      imul   eax,ebx
003d6d39  50                            push   eax
003d6d3a  6a 02                         push   2H
003d6d3c  5b                            pop    ebx
003d6d3d  58                            pop    eax
003d6d3e  8b d0                         mov    edx,eax
003d6d40  c1 fa 1f                      sar    edx,1fH
003d6d43  f7 fb                         idiv   eax,ebx
003d6d45  50                            push   eax
003d6d46  58                            pop    eax
003d6d47  5a                            pop    edx
003d6d48  59                            pop    ecx
003d6d49  5b                            pop    ebx
003d6d4a  c3                            ret

2007-06-26T11:32:00.000+01:00

A simple compiler In my last post I played a little with Harpy to generate code from what looked like assembly language. So the next logical step is to make a compiler. And I mean a real machine-code compiler, not some wimpy byte-code nonsense. To make life simple (it is a simple compiler, after all) I'm going to start by generating code that uses the stack all the time. So all operands will live on the stack and all operations on them will push and pop the stack. Of course, the code generated from such a compiler will make any serious compiler writer weep. The CodeGen type in Harpy is the monad used during code generation. It allows you to keep some extra information around besides the generated code; it gives you access to a state monad. I will use the state to keep track of the current stack depth. (We will see why soon.)

type StackDepth = Int
type Gen = CodeGen () StackDepth ()

addDepth :: StackDepth -> Gen
addDepth i = do
    d <- getState
    setState (d+i)

The addDepth function changes the current stack depth by grabbing the old one, adding the argument and storing it back. The getState and setState functions don't generate any code, they just manipulate the state available in the CodeGen monad. With that out of the way, let's implement code generation for addition.

gadd :: Gen
gadd = do
    pop  ebx
    pop  eax
    add  eax ebx
    push eax
    addDepth (-1)

It pops the two operands off the stack, adds them, and pushes the result. The net effect on the stack is that it has one word less on it (I count my stack depth in words), so there's also a call to addDepth. Subtraction and multiplication are very similar.

gsub :: Gen
gsub = do
    pop  ebx
    pop  eax
    sub  eax ebx
    push eax
    addDepth (-1)

gmul :: Gen
gmul = do
    pop  ebx
    pop  eax
    imul InPlace eax ebx
    push eax
    addDepth (-1)

On the i386 (signed) division and remainder is computed with the idiv instruction. It divides the EDX:EAX 64 bit number with the given operand. So we must convert a 32 signed number to a 64 bit signed number, this is simply done by moving EAX to EDX and then shifting right 31 steps. This will copy the sign bit into every bit of EDX. Depending of if we want the quotient or remainder we need to push EAX or EDX.

gquot :: Gen
gquot = do
    gquotRem
    push eax
    addDepth (-1)

grem :: Gen
grem = do
    gquotRem
    push edx
    addDepth (-1)

gquotRem :: Gen
gquotRem = do
    pop  ebx
    pop  eax
    mov  edx eax
    sar  edx (31 :: Word8)
    idiv ebx

To put a constant on the stack we simply push it and increment the remembered stack depth.

gconst :: Int -> Gen
gconst c = do
    push (fromIntegral c :: Word32)
    addDepth 1

OK, so now for something more interesting. Assuming we are generating code for a function, we also want to access the arguments to the function. Where are the arguments? Well, according to the IA32 calling conventions the caller pushes the arguments on the stack, so we'll follow those. First we have a bunch things on the stack, how many is kept track of in the stack depth in the CodeGen monad, and then the arguments follow in order, pushed right-to-left. So to get an argument we compute the offset, convert it to a byte offset, and push that word on the stack.

-- Get the Nth argument
gargN :: Int -> Gen
gargN n = do
    d <- getState
    let o = 4 * (d + n)
    mov  eax (Disp (fromIntegral o), esp)
    push eax
    addDepth 1

When generating code for a function we should not clobber any callee-save registers, so to be on the safe side we save all used registers on function entry and restore them on function exit. On function exit we also return the result in EAX.

savedRegs = [ebx, edx]

-- Push all register we'll be using, except eax.
gprologue :: Gen
gprologue = do
    mapM_ push savedRegs

-- Pop return value to eax, restore regs, and return
gepilogue :: Gen
gepilogue = do
    pop  eax
    mapM_ pop (reverse savedRegs)
    ret

OK, so that was a lot of stuff, let's put it together for a test.

testGen = conv_Word32ToWord32 () (length savedRegs + 1) $ do
    gprologue
    gargN 0
    gconst 1
    gadd
    gepilogue

main = do
    test <- testGen
    print (test 10)

The testGen function generates the prologue, push argument, push 1, add, and the epilogue. The conv_Word32ToWord32 (from my previous post) converts the machine code to a Haskell function. We also have to give the start value of the stack depth. The stack initially contains the return address and the saved registers, so that's the number we pass. Running this gives 11, as it should. OK, so let's actually write a compiler and not just a code generator. Here is a data type for integer expressions.

data Exp
    = Con Int
    | Arg Int
    | BinOp BOp Exp Exp
    deriving (Show)

data BOp = Add | Sub | Mul | Quot | Rem
    deriving (Show)

We have constants, arguments (variables), and a few binary operators. It's all easily translated to machine code.

translate :: Exp -> Gen
translate (Con c) = gconst c
translate (Arg n) = gargN n
translate (BinOp op x y) = do translate x; translate y; binop op
  where binop Add = gadd
        binop Sub = gsub
        binop Mul = gmul
        binop Quot = gquot
        binop Rem = grem

For simplicity, let's compile only functions of one argument for now.

compileIOW32 :: (Exp -> Exp) -> IO (Word32 -> Word32)
compileIOW32 f = conv_Word32ToWord32 () (length savedRegs + 1) $ do
    gprologue
    translate (f (Arg 0))
    gepilogue

This function takes an Exp->Exp function, by giving this function the argument Arg 0 we get an expression to translate. We tack on the usual prologue and epilogue. So let's try it.

main = do
    let fun x = BinOp Add x (Con 1)
    test <- compileIOW32 fun
    print (test 10)

Which prints 11. But yuck, writing BinOp etc. isn't nice. Let's make some instances.

instance Eq Exp
instance Ord Exp

instance Num Exp where
    x + y  =  BinOp Add x y
    x - y  =  BinOp Sub x y
    x * y  =  BinOp Mul x y
    fromInteger i = Con (fromInteger i)

instance Enum Exp
instance Real Exp

instance Integral Exp where
    quot x y = BinOp Quot x y
    rem x y = BinOp Rem x y

And give it a whirl:

main = do
    let fun x = (x+1) * x `quot` 2
    test <- compileIOW32 fun
    print (test 10)

And this prints 55, as expected. Not bad, a compiler from (a tiny subset of) Haskell functions to machine code in a few pages. But I do admit being embarrassed about generating such incredibly poor code. But there's always future blog posts to rectify that.

2007-06-25T18:34:00.000+01:00

Playing with Harpy Recently there was an announcement of the Harpy package for Haskell. Harpy is a very cool package that lets you generate machine code from a Haskell program into a memory buffer. Using the Haskell FFI interface you can then turn the machine code into a Haskell function and call it. The way that you generate machine code from Harpy looks very much like assembly code. Oh, btw, apologies to those who are not using x86 machines, this blog post is very machine dependent, since Harpy generate x86 machine code. A small Harpy example First, some sample code (stolen from the Harpy tutorial). It might look like assembly code, but it's actually Haskell.

asm_fac = do
    loopTest  <- newLabel
    loopStart <- newLabel
    ensureBufferSize 160
    push ecx
    mov  ecx (Disp 8, esp)
    mov  eax (1 :: Word32)
    jmp  loopTest
    loopStart @@ mul ecx
    sub  ecx (1 :: Word32)
    loopTest @@ cmp ecx (0 :: Word32)
    jne  loopStart
    pop  ecx
    ret

Just some comments: This is (obviously) monadic code. The newLabel operation creates a new label that later has to be defined and that can be used as a jump target. The ensureBufferSize is a tedious function to make sure there is enough space to emit the instructions. I hope the next version of Harpy will not require this, since it's really the computers job to count bytes, not mine (and it's not hard to do). Apart from that, most of the code should be obvious for anyone who has done assembly programming on an x86. To generate code and then convert the generated code we need some utility functions.

type Importer f = FunPtr f -> f
foreign import ccall safe "dynamic" import_Word32ToWord32 :: Importer (Word32 -> Word32)

conv_Word32ToWord32 :: e -> s -> CodeGen e s a -> IO (Word32 -> Word32)
conv_Word32ToWord32 = conv import_Word32ToWord32

conv :: (Importer f) -> e -> s -> CodeGen e s a -> IO f
conv imp e s cg = do
    let cg' = do cg; getEntryPoint
    (_, r) <- runCodeGen cg' e s
    case r of
        Left msg -> error (show msg)
        Right addr -> return $ imp $ castPtrToFunPtr addr

The important function is conv_Word32ToWord32 which takes a block like asm_fac and turns it into a Haskell function. The e and s arguments we can ignore for now, they are for more advanced uses of the machine code generator. Note that there is no type safety here. The asm_fac block has nothing that says what kind of function it might implement, so it's all up to us to use it correctly. (Once you use FFI Haskell isn't any safer than C.) Finally, to use it all:

main = do
    fac <- conv_Word32ToWord32 () () asm_fac
    print (fac 5, fac 10)

Pretty nifty, huh? Calling conv_Word32ToWord32 with the asm_fac argument will emit machine code for the factorial function into a memory buffer and then wrap it up so it can be called as a regular Haskell function. Running the program will print (120,3628800) as expected. Final thoughts Harpy is really cool and uses Haskell type classes in a clever way to make Harpy "assembly" code look almost like normal assembly code. But there is room for improvement. Currently you have no control over the buffer handling; Harpy uses mallocBytes internally. I would like to see the buffer handling (mallocBytes, peek, the IO monad) abstracted out. That way I don't have to generate code into a memory buffer if I don't want to. I could generate code in, e.g., a list. Or I could not generate code at all, but just count bytes. Or I could allocate buffers with something better than mallocBytes. On the whole, I think Harpy is great for experimenting with code generation. More about that soon.

2007-05-07T13:27:00.000+01:00

I was asked if my fixed number module (Data.Number.Fixed) which has the precision built into the type could be used with a dynamic epsilon, i.e., one that is not know at compile time. The answer is, yes. It only takes a fraction of an oleg to do. A recap: The idea with the Fixed type was that, e.g., Fixed Eps1, is the type where computation happens with an epsilon of 1. To get the epsilon for a type I have a class

class Epsilon e where
    eps :: e -> Rational
instance Epsilon Eps1 where eps _ = 1
...

So what do we do if we need an epsilon of say, 0.01? Well, there's another instance for the type constructor EpsDiv10:

instance (Epsilon e) => Epsilon (EpsDiv10 e) where
    eps x = eps x / 10

So the type EpsDiv10 (EpsDiv10 Eps1) will have an epsilon of 0.01. Similarly, with the right set of primitive types and type constructors you could build any epsilon you want. To simplify a little the function below only finds an epsilon within a factor of 10 from the requested one.

dynamicEps :: forall a . Rational ->
             (forall e . Epsilon e => Fixed e -> a) ->
             (Rational -> a)
dynamicEps r f v = loop (undefined :: Eps1)
 where loop :: forall x . (Epsilon x) => x -> a
       loop e = if eps e <= r then f (fromRational v :: Fixed x)
                else loop (undefined :: EpsDiv10 x)

This function takes a desired (rational) epsilon, r, and a function, f, and returns a function that behaves like f, but that will round its argument to the desired epsilon and compute with that.

Main> putStrLn $ dynamicEps 1e-40 (show . sin) 1
0.8414709848078965066525023216302989996226

How does it work? It's starts by passing something of type Eps1 to loop, and then the loop uses that epsilon if it is small enough, otherwise it will pass something that have an epsilon that is a tenth to the next iteration of the loop. There are two "tricks" here, first the loop function recurses with an argument of a different type than it was given. So it needs polymorphic recursion, and thus needs a type signature. Second, the function, f, has a rank two type. This way the function can work with any epsilon it is given; otherwise we would not know that. Note how the argument to loop is actually never ysed, so undefined works great. All we need is the right type. And with scoped type variables this is easy to do. There are some improvements that should be made to this function: it is inefficient and it doesn't give you the exact epsilon back. But it has the right general idea.

2007-04-28T11:48:00.000+01:00

Frustration My trusty MacBook suddenly died, :( It can no longer boot, and the old way of getting an OpenFirmware prompt doesn't seem to work. Oh well.

2007-04-27T20:24:00.000+01:00

Fixed precision, an update So I was a bit sloppy in my last post. When doing arithmetic it was performed exactly using Rational and not truncated according to the epsilon for the type. So, for instance, computing 4/10 + 4/10 with type Fixed Eps1 would give the answer 1. While this might seem nice, it's also wrong if every operation would be performed with epsilon 1, since 4/10 would be 0, and 0+0 is 0. So I'll amend my implementation a little.

instance (Epsilon e) => Num (Fixed e) where
    (+) = lift2 (+)
    (-) = lift2 (-)
    (*) = lift2 (*)
    negate (F x) = F (negate x)
    abs (F x) = F (abs x)
    signum (F x) = F (signum x)
    fromInteger = F . fromInteger

instance (Epsilon e) => Fractional (Fixed e) where
    (/) = lift2 (/)
    fromRational x = r
        where r = F $ approx x (getEps r)

lift2 :: (Epsilon e) => (Rational -> Rational -> Rational) -> Fixed e -> Fixed e
 -> Fixed e
lift2 op fx@(F x) (F y) = F $ approx (x `op` y) (getEps fx)

approx :: Rational -> Rational -> Rational
approx x eps = approxRational (x + eps/2) eps

So after each operation we add half an epsilon (so we get rounding instead of truncation) and call the magical approxRational to get the closest rational within an epsilon.

2007-04-25T23:22:00.000+01:00

Overloading Haskell numbers, part 3, Fixed Precision. Fixed precision numbers can be handy to have sometimes. Haskell only provides integral, rational, and floating point in the Prelude, but I want more! I want to define a type of fixed precision numbers, i.e., numbers that have a certain number of decimals. Rational numbers are handy to represent these numbers so here we go.

newtype Fixed = F Rational deriving (Eq, Ord, Enum, Num, Fractional, Real, RealFrac)

This uses an GHC extension that allows you to derive anything from the representation type so we don't have to define all those tedious instances. But wait, where is the fixed precision? What I just defined is just the Rational numbers with a different name. I really want a type which is parametrized over what precision it has. OK, so here we go again:

newtype Fixed e = F Rational deriving (Eq, Ord, Enum, Num, Fractional, Real, RealFrac)

Here e is meant to be a type that specifies the precision. Let's capture those types that we allow as the precision as a type class:

class Epsilon e where
    eps :: e -> Rational

The idea here is that the actual argument to the eps function is not used, it's just a proxy and the type carries all the information. Here's an example:

data Eps1
instance Epsilon Eps1 where
    eps _ = 1

The eps returns 1 for the Eps1 type. So the intended meaning is that the type Fixed Eps1 are numbers within an epsilon of the right answer, and epsilon is 1. Similarly

data Prec50
instance Epsilon Prec50 where
    eps _ = 1e-50

The type Fixed Prec50 has an epsilon of 1e-50, i.e., 50 decimals. Btw, notice that what looks like a floating point literal is actually a rational number; one of the very clever decisions in the original Haskell design. We can also define type constructors like

data EpsDiv10 p
instance (Epsilon e) => Epsilon (EpsDiv10 e) where
    eps e = eps (un e) / 10
       where un :: EpsDiv10 e -> e
             un = undefined

The type Fixed (EpsDiv10 Prec50) has an epsilon of 1e-50/10 = 1e-51. Note how the function un is just a dummy to get the e. If we had used the scoped type variable extension we could have avoided this function. OK, so we can compute with these fixed precision numbers, but where does the epsilon enter? Nowhere, so far. All arithmetic is performed exactly by the rationals, but our first use of the epsilon is in the show function.

instance (Epsilon e) => Show (Fixed e) where
    showsPrec = showSigned showFixed
      where showFixed f@(F x) = showString $ show q ++ "." ++ decimals r (eps $ un f)
              where q :: Integer
             (q, r) = properFraction x
      un :: Fixed e -> e
             un = undefined
             decimals a e | e >= 1 = ""
                          | otherwise = intToDigit b : decimals c (10 * e)
                 where (b, c) = properFraction (10 * a)

To print a fixed point number we first print the integer part. The properFraction method handily returns the integer part and some left over fraction that is <1. We then print decimal after decimal by multiplying the fraction by 10, converting the integer part to a digit, and then recursing. While generating these decimals we also multiply the epsilon by 10 each time and stop printing when epsilon is >= 1, because then we don't need any more decimals. All right, lets put it to the test:

Data.Number.Fixed> 1/3 :: Fixed Prec50
0.33333333333333333333333333333333333333333333333333
Data.Number.Fixed> 3/8 :: Fixed Prec10
0.3750000000
Data.Number.Fixed> (1 :: Fixed Prec10) + (0.5 :: Fixed Prec50)
    Couldn't match expected type `Prec10'
           against inferred type `Prec50'
      Expected type: Fixed Prec10
      Inferred type: Fixed Prec50

As we want, we can't mix different fixed point numbers with different precisions. Sometimes we do want to convert, so let's provide such a function

convertFixed :: Fixed e -> Fixed f
convertFixed (F x) = F x

It looks like this function does nothing, and it's true. It just converts the type, and since the epsilon type is not involved in the number itself it's trivial. Well, basic arithmetic isn't that much fun. What about the transcendental functions? Luckily Jan Skibinski has already implemented all of them. They compute the various function to with an epsilon using continued fractions. The instance declaration looks something like:

instance (Epsilon e) => Floating (Fixed e) where
    sqrt = toFixed1 F.sqrt
    exp = toFixed1 F.exp
    ...

toFixed1 :: (Epsilon e) => (Rational -> Rational -> Rational) -> Fixed e -> Fixed e
toFixed1 f x@(F r) = F (f (eps (un x)) r)
  where un :: Fixed e -> e
        un = undefined

And all the clever code I stole from Jan is in the module F, see link below. Now it's getting more fun.

Data.Number.Fixed> pi :: Fixed Prec50
3.14159265358979323846264338327950288419716939937510
Data.Number.Fixed> exp 10 :: Fixed Prec10
22026.4657881843
Data.Number.Fixed> log 2 :: Fixed Prec50
0.69314718055994530941723212145817656807550013436026
Data.Number.Fixed> exp (sin (sqrt (54/42))) :: Fixed Prec10
2.4745696343
Data.Number.Fixed> sqrt 2 :: Fixed Prec500
1.41421356237309504880168872420969807856967187537694807317667973799073247
8462107038850387534327641572735013846230912297024924836055850737212644121
4970999358314132226659275055927557999505011527820605714701095599716059702
7453459686201472851741864088919860955232923048430871432145083976260362799
5251407989687253396546331808829640620615258352395054745750287759961729835
5752203375318570113543746034084988471603868999706990048150305440277903164
5424782306849293691862158057846311159666871301301561856898723724

Data.Number.Fixed Data.Complex> let i = 0 :+ 1 :: Complex (Fixed Prec50) 
Data.Number.Fixed Data.Complex> exp (i*pi)
(-0.99999999999999999999999999999999999999999999999999) :+ 0.00000000000000000000000000000000000000000000000000

Naturally, we need to package up all these definitions in a module that will keep the Fixed type abstract, and you can also keep the Epsilon class and the various Epsilon types abstract if you provide enough building blocks. Source at http://www.augustsson.net/Darcs/Data.Number/.

2007-04-20T00:51:00.000+01:00

Following a suggestion from Cale Gibbard I added a convenient feature to Djinn. You can now define type classes and have contexts on functions. As with Djinn in general there is no polymorphic instantiation, so the methods of a class must not mention any type variables, but the class parameters. This makes the whole thing rather limited, but it was a quick hack. Sample session

Welcome to Djinn version 2007-04-20.
Type :h to get help.
Djinn> f ? (Eq a) => a -> a -> (b,c) -> Either b c
f :: (Eq a) => a -> a -> (b, c) -> Either b c
f a b =
    case a == b of
    False -> \ (c, _) -> Left c
    True -> \ (_, d) -> Right d

The Djinn source is at darcs.augustsson.net/Darcs/Djinn as usual. Thanks Cale!

2007-04-14T23:42:00.000+01:00

Overloading Haskell numbers, part 2, Forward Automatic Differentiation. I will continue my overloading by some examples that have been nicely illustrated by an article by Jerzy Karczmarczuk. And blogged about by sigfpe. But at least I'll end this entry with a small twist I've not seen before. When computing the derivative of a function you normally do by either symbolic derivation, or by a numerical approximation. Say that you have a function f(x) = x² + 1 and you want to know the derivative at x=5. Doing it symbolically you first get f' f'(x) = 2x using high school calculus (maybe they don't teach it in high school anymore?), and then you plug in 5 f'(5) = 2*5 = 10 Computing it by numeric differentiation you compute f'(x) = (f(x+h) - f(x)) / h for some small h. Let's pick h=1e-5, and we get f'(5) = 10.000009999444615. Close, but not that good. So why don't we always use the symbolic method? Well, some functions are not that easy to differentiate. Take this one g x = if abs (x - 0.7) < 0.4 then x else g (cos x) What's the derivative? Well, it's tricky because this is not really a proper definition of g. It's an equation that if solved will yield a definition of g. And like equations in general, it could have zero, one, or many solutions. (If we happen to use CPOs there is always a unique smallest solution which is what programs compute, as if by magic.) If you think g is contrived, lets pick a different example: computing the square root with Newton-Raphson.

sqr x = convAbs $ iterate improve 1
  where improve r = (r + x/r) / 2
        convAbs (x1:x2:_) | abs (x1-x2) < 1e-10 = x2
 convAbs (_:xs) = convAbs xs

So symbolic is not so easy here, and numeric differentiation is not very accurate. But there is a third way! Automatic differentiation. The idea behind AD is that instead of computing with with just numbers, we instead compute with pairs of numbers. The first component is the normal number, and the second component is the derivative. What are the rules for these numbers? Let's look at addition (x, x') + (y, y') = (x+y, x'+y') To add two numbers you just add the regular part and the derivatives. For multiplication you have to remember how to compute the derivative of a product: (f(x)*g(x))' = f(x)*g'(x) + f'(x)*g(x) So for our pairs we get (x, x') * (y, y') = (x*y, x*y' + x'*y) i.e., first the regular product, then the derivative according to the recipe above. Let's see how it works on f(x) = x² + 1 We want the derivative at x=5. So what is the pair we use for x? It is (5, 1). Why? Well it has to be 5 for the regular part, and since this represents x and the derivative of x is 1, the pair is (5, 1). In the right hand side for f we need to replace 1 by (1,0), since the derivative of a constant is 0. So then we get f (5,1) = (5,1)*(5,1) + (1,0) = (26,10) using the rules above. And look! There is the normal result, 26, as well as the derivative, 10. Let's turn this into Haskell, using the type PD to hold a pair of Doubles

data PD = P Double Double deriving (Eq, Ord, Show)
instance Num PD where
    P x x' + P y y' = P (x+y) (x'+y')
    P x x' - P y y' = P (x-y) (x'-y')
    P x x' * P y y' = P (x*y) (x*y' + y'*x)
    fromInteger i   = P (fromInteger i) 0

A first observation is that there is nothing Double specific in this definitions; it would work for any Num. So we can change it to

data PD a = P a a deriving (Eq, Ord, Show)
instance Num a => Num (PD a) where ...

Let's also add abs&signum and the Fractional instance

    ...
    abs (P x x') = P (abs x) (signum x * x')
    signum (P x x') = P (signum x) 0

instance Fractional a => Fractional (PD a) where
    P x x' / P y y' = P (x / y) ( (x'*y - x*y') / (y * y))
    fromRational r  = P (fromRational r) 0

We can now try the sqr example

Main> sqr (P 9 1)
P 3.0 0.16666666666666666

The derivative of x**0.5 is 0.5*x**(-0.5), i.e., 0.5*9**(-0.5) = 0.5/3 = 0.16666666666666666. So we got the right answer. BTW, if you want to be picky the derivative of signum is not 0. The signum function makes a jump from -1 to 1 at 0. So the "proper" value would be 2*dirac, if dirac is a Dirac pulse. But since we don't have numbers with Dirac pulses (yet), I'll just pretend the derivative is 0 everywhere. The very clever insight that Jerzy had was that when doing these numbers in Haskell there is no need to limit yourself to just the first derivative. Since Haskell is lazy we can easily keep an infinite list of all derivatives instead of just the first one. Let's look at how that definition looks. It's very similar to what we just did. But instead of the derivative being just a number, it's now one of our new numbers with a value, and all derivatives... Since we are now dealing with an infinite data structure we need to define our own show, (==), etc.

data Dif a = D a (Dif a)

val (D x _) = x

df (D _ x') = x'

dVar x = D x 1

instance (Show a) => Show (Dif a) where
    show x = show (val x) 

instance (Eq a) => Eq (Dif a) where
    x == y  =  val x == val y

instance (Ord a) => Ord (Dif a) where
    x `compare` y  =  val x `compare` val y

instance (Num a) => Num (Dif a) where
    D x x' + D y y'          =  D (x + y) (x' + y')
    D x x' - D y y'          =  D (x - y) (x' - y')
    p@(D x x') * q@(D y y')  =  D (x * y) (x' * q + p * y')
    fromInteger i            =  D (fromInteger i) 0
    abs p@(D x x')           =  D (abs x) (signum p * x')
    signum (D x _)           =  D (signum x) 0

instance (Fractional a) => Fractional (Dif a) where
    recip (D x x') = ip
 where ip = D (recip x) (-x' * ip * ip)
    fromRational r = D (fromRational r) 0

This looks simple, but it's rather subtle. For instance, take the 0 in the definition of fromInteger. It's actually of Dif type, so it's a recursive call to fromInteger. So let's try with our sqr function again, this time computing up to the third derivative. The dVar is used to create a value for "variable" where we want to differentiate.

Main> sqr $ dVar 9
3.0
Main> df $ sqr $ dVar 9
0.16666666666666669
Main> df $ df $ sqr $ dVar 9
-9.259259259259259e-3
Main> df $ df $ df $ sqr $ dVar 9
1.5432098765432098e-3

And the transcendentals in a similar way:

lift (f : f') p@(D x x') = D (f x) (x' * lift f' p)

instance (Floating a) => Floating (Dif a) where
    pi               = D pi 0
    exp (D x x')     = r where r = D (exp x) (x' * r)
    log p@(D x x')   = D (log x) (x' / p)
    sqrt (D x x')    = r where r = D (sqrt x) (x' / (2 * r))
    sin              = lift (cycle [sin, cos, negate . sin, negate . cos])
    cos              = lift (cycle [cos, negate . sin, negate . cos, sin])
    acos p@(D x x')  = D (acos x) (-x' / sqrt(1 - p*p))
    asin p@(D x x')  = D (asin x) ( x' / sqrt(1 - p*p))
    atan p@(D x x')  = D (atan x) ( x' / (p*p - 1))
    sinh x           = (exp x - exp (-x)) / 2
    cosh x           = (exp x + exp (-x)) / 2
    asinh x          = log (x + sqrt (x*x + 1))
    acosh x          = log (x + sqrt (x*x - 1))
    atanh x          = (log (1 + x) - log (1 - x)) / 2

And why not try the function g we defined above?

Main> g 10
0.6681539175313869
Main> g (dVar 10)
0.6681539175313869
Main> df $ g (dVar 10)
0.4047642621121782
Main> df $ df $ g (dVar 10)
0.4265424381635987
Main> df $ df $ df $ g (dVar 10)
-1.4395397945007182

It all works very nicely. So now when we can compute the derivative of a function, let's define something somewhat more interesting with it. Let's revisit the sqr function again. It uses Newton-Raphson to find the square root. How does Newton-Raphson actually work? Given a differentiable function, f(x), it finds a zero by starting with some x₁ and then iterating x_n+1 = x_n - f(x_n)/f'(x_n) until we meet some convergence criterion. Using this, let's define a function that finds a zero of another function:

findZero f = convRel $ cut $ iterate step start
    where step x = x - val fx / val (df fx) where fx = f (dVar x)

   start = 1  -- just some value
   epsilon = 1e-10
   cut = (++ error "No convergence in 1000 steps") . take 1000
   convRel (x1:x2:_) | x1 == x2 || abs (x1+x2) / abs (x1-x2) > 1/epsilon = x2
   convRel (_:xs) = convRel xs

The only interesting part is the step function that does one iteration with Newton-Raphson. It computes f x and then divides the normal value with the derivative. We then produce the infinite list of approximations using step, then cut it of at some point (in case it doesn't converge), and then we look down the list for two values that are within some relative epsilon. And it even seems to work.

Main> findZero (\x -> x*x - 9)
3.0
Main> findZero (\x -> sin x - 0.5)
0.5235987755982989
Main> sin it
0.5
Main> findZero (\x -> x*x + 9)
*** Exception: No convergence in 1000 steps
Main> findZero (\x -> sqr x - 3)
9.0

Note how it finds a zero of the sqr function which is actually using recursion internally to compute the square root. So now we can compute numerical derivatives. But wait! We also have symbolic numbers. Can we combine them? Of course, that is the power of polymorphism. Let's load up both modules:

Data.Number.Symbolic Dif3> let x :: Num a => Dif (Sym a); x = dVar (var "x")
Data.Number.Symbolic Dif3> df $ x*x
x+x
Data.Number.Symbolic Dif3> df $ sin x
cos x
Data.Number.Symbolic Dif3> df $ sin (exp (x - 4) * x)
(exp (-4.0+x)*x+exp (-4.0+x))*cos (exp (-4.0+x)*x)
Data.Number.Symbolic Dif3> df $ df $ sin (exp (x - 4) * x)
(exp (-4.0+x)*x+exp (-4.0+x)+exp (-4.0+x))*cos (exp (-4.0+x)*x)+(exp (-4.0+x)*x+exp (-4.0+x))*(exp (-4.0+x)*x+exp (-4.0+x))*(-sin (exp (-4.0+x)*x))

We define x to be a differentiable number, "the variable", over symbolic numbers, over some numbers. And then we just happily use df to get the differentiated versions. So we set out to compute numeric derivatives, and we got these for free. Not too bad. One final note, the Dif type is defined above can be made more efficient by not keeping all the infinite tails with 0 derivatives around. In a real module for this, you'd want to make this optimization. [Edit: fixed typo.]

2007-04-11T23:02:00.000+01:00

Overloading Haskell numbers, part 1, symbolic expressions. Haskell's overloaded numerical classes can be (ab)used to do some symbolic maths. This is in no way a new discovery, but I thought I'd write a few lines about it anyway since I've been playing with it the last few days. First we need a data type to represent expressions. We want constants, variables, and function applications. But we don't want to fix the type of the constants, so that will be a parameter to the type.

data Sym a = Con a | Var String | App String [Sym a]
    deriving (Eq, Show)

And we also take the opportunity to derive Eq and Show. Now we can actually claim that the type Sym N is a number if N is a number. Let do it:

instance (Num a) => Num (Sym a) where
    x + y         = App "+" [x, y]
    x - y         = App "-" [x, y]
    x * y         = App "*" [x, y]
    negate x      = App "negate" [x]
    abs    x      = App "abs"    [x]
    signum x      = App "signum" [x]
    fromInteger x = Con (fromInteger x)

A small interactive session shows that we are on the right track.

Sym1> let x = Var "x"
Sym1> x*x + 5
App "+" [App "*" [Var "x",Var "x"],Con 5]
Sym1> x*x*x
App "*" [App "*" [Var "x",Var "x"],Var "x"]
Sym1> 2 + 3 :: Sym Int
App "+" [Con 2,Con 3]

We can type in normal looking expressions, but when they are printed the Show instance is used so we get to see the raw syntax tree. That has its uses, but it gets old quickly. We want a pretty printer. To get the precedences right we need to define showsPrec and pass it the right arguments. It's a little tedious, but nothing strange.

instance (Show a) => Show (Sym a) where
    showsPrec p (Con c) = showsPrec p c
    showsPrec _ (Var s) = showString s
    showsPrec p (App op@(c:_) [x, y]) | not (isAlpha c) =
        showParen (p>q) (showsPrec ql x . showString op . showsPrec qr y)
        where (ql, q, qr) = fromMaybe (9,9,9) $ lookup op [
                   ("**", (9,8,8)),
                   ("/",  (7,7,8)),
                   ("*",  (7,7,8)),
                   ("+",  (6,6,7)),
                   ("-",  (6,6,7))]
    showsPrec p (App "negate" [x]) =
        showParen (p>=6) (showString "-" . showsPrec 7 x)
    showsPrec p (App f xs) =
        showParen (p>10) (foldl (.) (showString f)
                                (map (\ x -> showChar ' ' . showsPrec 11 x) xs))

Let's try the same examples again:

Sym2> let x = var "x"
Sym2> x*x + 5
x*x+5
Sym2> x*x*x
x*x*x
Sym2> 2 + 3 :: Sym Int
2+3

Look we can type expressions and get them back again! The instance Num (Sym a) isn't too bad, the only fishy thing about it is the Eq superclass that is required for Num. We have Eq for Sym, but it doesn't really behave like it should. E.g., the expression 'x==1' would come out as False since the syntax trees are not equal. But this isn't really what we would like, ideally (==) would also turn into something symbol, but that is impossible with the standard Prelude. Let's make some more instances. A few of these definitions are just there to appease the Haskell numerical hierarchy and supply some operations it need.

instance (Fractional a) => Fractional (Sym a) where
    x / y          = App "/" [x, y]
    fromRational x = Con (fromRational x)

instance (Real a) => Real (Sym a) where
    toRational (Con c) = toRational c

instance (RealFrac a) => RealFrac (Sym a) where
    properFraction (Con c) = (i, Con c') where (i, c') = properFraction c

instance (Floating a) => Floating (Sym a) where
    pi = App "pi" []
    exp = app1 "exp"
    sqrt = app1 "sqrt"
    log = app1 "log"
    (**) = app2 "**"
    logBase = app2 "logBase"
    sin = app1 "sin"
    tan = app1 "tan"
    cos = app1 "cos"
    asin = app1 "asin"
    atan = app1 "atan"
    acos = app1 "acos"
    sinh = app1 "sinh"
    tanh = app1 "tanh"
    cosh = app1 "cosh"
    asinh = app1 "asinh"
    atanh = app1 "atanh"
    acosh = app1 "acosh"

instance (RealFloat a) => RealFloat (Sym a) where
    exponent _ = 0
    scaleFloat 0 x = x
    atan2 = app2 "atan2"

app1 :: String -> Sym a -> Sym a
app1 f x = App f [x]

app2 :: String -> Sym a -> Sym a -> Sym a
app2 f x y = App f [x, y]

Let's put this code to the test by bringing the Complex number module into scope.

Sym3> :m +Data.Complex
Sym3 Data.Complex> let x=Var "x"; y=Var "y"
Sym3 Data.Complex> sin (x:+y)
sin x*cosh y :+ cos x*sinh y

And by that last expression we have recovered the definition of complex sin as it is given in the Data.Complex module. Let's try another one.

Sym3 Data.Complex> asinh(x:+y)
log (sqrt ((x+sqrt ((sqrt ((1.0+(x*x-y*y))*(1.0+(x*x-y*y))+(0.0+(x*y+y*x))*(0.0+
 (x*y+y*x)))+abs (1.0+(x*x-y*y)))/2.0))*(x+sqrt ((sqrt ((1.0+(x*x-y*y))*(1.0+(x*x
 -y*y))+(0.0+(x*y+y*x))*(0.0+(x*y+y*x)))+abs (1.0+(x*x-y*y)))/2.0))+(y+abs (0.0+
 (x*y+y*x))/(sqrt ((sqrt ((1.0+(x*x-y*y))*(1.0+(x*x-y*y))+(0.0+(x*y+y*x))*(0.0+
 (x*y+y*x)))+abs (1.0+(x*x-y*y)))/2.0)*2.0))*(y+abs (0.0+(x*y+y*x))/(sqrt ((sqrt
 ((1.0+(x*x-y*y))*(1.0+(x*x-y*y))+(0.0+(x*y+y*x))*(0.0+(x*y+y*x)))+abs (1.0+(x*x
 -y*y)))/2.0)*2.0)))) :+ atan2 (y+abs (0.0+(x*y+y*x))/(sqrt ((sqrt ((1.0+(x*x-y*y))
 *(1.0+(x*x-y*y))+(0.0+(x*y+y*x))*(0.0+(x*y+y*x)))+abs (1.0+(x*x-y*y)))/2.0)*2.0))
 (x+sqrt ((sqrt ((1.0+(x*x-y*y))*(1.0+(x*x-y*y))+(0.0+(x*y+y*x))*(0.0+(x*y+y*x)))+
 abs (1.0+(x*x-y*y)))/2.0))

Hmmmm, that might be right, but it's rather ugly. There's also a lot of '0.0+...' in that expression. We need something that can simplify expressions. It would also be nice if all constant expressions were evaluated instead of stored. To achieve this we are going to change the representation a little. The App constructor will store the real function used to work on constants as well as the name of it. And while we are at it, we'll get rid of the Var constructor. We might as well use App with an empty argument list. Furthermore, since this is starting to look useful, we'll give the module a proper name and export only the interface we want to be visible. We will hide the details of the Sym type and just export some accessor functions. The simplification happens in the binOp and unOp functions. I have just put some algebraic laws there (assuming the underlying numeric type is a field). The list of rewrites performed by these functions is far from complete. It's just a few that I found useful. Note how the code in binOp pattern matches on constants like 0, 1, and -1 directly. This actually works because of the semantics of Haskell pattern matching against numeric literals. Also note that the constraint on `a' is just Num, even though we do some simplifications with (/) which belongs in Fractional. The instance declarations have been extended somewhat so that constant expressions in the Sym type will behave as the corresponding expressions in the underlying type. A small final run

*Data.Number.Symbolic Data.Complex> 1+x+2
3+x
*Data.Number.Symbolic Data.Complex> 1+x*(y-y)-1
0
*Data.Number.Symbolic Data.Complex> sin(x:+1e-10)
sin x :+ 1.0e-10*cos x
*Data.Number.Symbolic Data.Complex> asinh(x:+y)
log (sqrt ((x+sqrt ((sqrt ((1.0+(x*x-y*y))*(1.0+(x*x-y*y))+(x*y+y*x)*(x*y+y*x))+abs
 (1.0+(x*x-y*y)))/2.0))*(x+sqrt ((sqrt ((1.0+(x*x-y*y))*(1.0+(x*x-y*y))+(x*y+y*x)*
 (x*y+y*x))+abs (1.0+(x*x-y*y)))/2.0))+(y+abs (x*y+y*x)/(2.0*sqrt ((sqrt ((1.0+(x*x
 -y*y))*(1.0+(x*x-y*y))+(x*y+y*x)*(x*y+y*x))+abs (1.0+(x*x-y*y)))/2.0)))*(y+abs (x*y
 +y*x)/(2.0*sqrt ((sqrt ((1.0+(x*x-y*y))*(1.0+(x*x-y*y))+(x*y+y*x)*(x*y+y*x))+abs
 (1.0+(x*x-y*y)))/2.0))))) :+ atan2 (y+abs (x*y+y*x)/(2.0*sqrt ((sqrt ((1.0+(x*x
 -y*y))*(1.0+(x*x-y*y))+(x*y+y*x)*(x*y+y*x))+abs (1.0+(x*x-y*y)))/2.0))) (x+sqrt
 ((sqrt ((1.0+(x*x-y*y))*(1.0+(x*x-y*y))+(x*y+y*x)*(x*y+y*x))+abs (1.0+(x*x-y*y)))/2.0))

As the final example shows, there is still a lot to do. Also note how the underlying numeric type has defaulted to Double, and we have a loss of precision in the second to last example. But an implementation of real numbers instead of floating point numbers will have to wait until a later posting.

module Data.Number.Symbolic(Sym, var, con, subst, unSym) where

import Data.Char(isAlpha)
import Data.Maybe(fromMaybe)
import Debug.Trace

data Sym a = Con a | App String ([a]->a) [Sym a]

instance (Eq a) => Eq (Sym a) where
    Con x      == Con x'        =  x == x'
    App f _ xs == App f' _ xs'  =  (f, xs) == (f', xs')
    _          == _             =  False

instance (Ord a) => Ord (Sym a) where
    Con x      `compare` Con x'        =  x `compare` x'
    Con _      `compare` App _ _ _     = LT
    App _ _ _  `compare` Con _         = GT
    App f _ xs `compare` App f' _ xs'  =  (f, xs) `compare` (f', xs')

var :: String -> Sym a
var s = App s undefined []

con :: a -> Sym a
con = Con

subst :: (Num a) => String -> Sym a -> Sym a -> Sym a
subst _ _ e@(Con _) = e
subst x v e@(App x' _ []) | x == x' = v
                         | otherwise = e
subst x v (App s f es) =
    case map (subst x v) es of
    [e] -> unOp (\ x -> f [x]) s e
    [e1,e2] -> binOp (\ x y -> f [x,y]) e1 s e2
    es' -> App s f es'

unSym :: (Show a) => Sym a -> a
unSym (Con c) = c
unSym e = error $ "unSym called: " ++ show e

instance (Show a) => Show (Sym a) where
    showsPrec p (Con c) = showsPrec p c
    showsPrec _ (App s _ []) = showString s
    showsPrec p (App op@(c:_) _ [x, y]) | not (isAlpha c) =
        showParen (p>q) (showsPrec ql x . showString op . showsPrec qr y)
        where (ql, q, qr) = fromMaybe (9,9,9) $ lookup op [
                   ("**", (9,8,8)),
     ("/",  (7,7,8)),
     ("*",  (7,7,8)),
     ("+",  (6,6,7)),
     ("-",  (6,6,7))]
    showsPrec p (App "negate" _ [x]) =
        showParen (p>=6) (showString "-" . showsPrec 7 x)
    showsPrec p (App f _ xs) =
        showParen (p>10) (foldl (.) (showString f) (map (\ x -> showChar ' ' . showsPrec 11 x) xs))

instance (Num a) => Num (Sym a) where
    x + y         = binOp (+) x "+" y
    x - y         = binOp (-) x "-" y
    x * y         = binOp (*) x "*" y
    negate x      = unOp negate "negate" x
    abs    x      = unOp abs    "abs"    x
    signum x      = unOp signum "signum" x
    fromInteger x = Con (fromInteger x)

instance (Fractional a) => Fractional (Sym a) where
    x / y          = binOp (/) x "/" y
    fromRational x = Con (fromRational x)

-- Assume the numbers are a field and simplify a little
binOp :: (Num a) => (a->a->a) -> Sym a -> String -> Sym a -> Sym a
binOp f (Con x) _ (Con y) = Con (f x y)
binOp _ x "+" 0 = x
binOp _ 0 "+" x = x
binOp _ x "+" (App "+" _ [y, z]) = (x + y) + z
binOp _ x "+" y | isCon y && not (isCon x) = y + x
binOp _ x "+" (App "negate" _ [y]) = x - y
binOp _ x "-" 0 = x
binOp _ x "-" x' | x == x' = 0
binOp _ x "-" (Con y) | not (isCon x) = Con (-y) + x
binOp _ _ "*" 0 = 0
binOp _ x "*" 1 = x
binOp _ x "*" (-1) = -x
binOp _ 0 "*" _ = 0
binOp _ 1 "*" x = x
binOp _ (-1) "*" x = -x
binOp _ x "*" (App "*" _ [y, z]) = (x * y) * z
binOp _ x "*" y | isCon y && not (isCon x) = y * x
binOp _ x "*" (App "/" f [y, z]) = App "/" f [x*y, z]
{-
binOp _ x "*" (App "+" _ [y, z]) = x*y + x*z
binOp _ (App "+" _ [y, z]) "*" x = y*x + z*x
-}
binOp _ x "/" 1 = x
binOp _ x "/" (-1) = -x
binOp _ x "/" x' | x == x' = 1
binOp _ x "/" (App "/" f [y, z]) = App "/" f [x*z, y]
binOp f x op y = App op (\ [a,b] -> f a b) [x, y]

unOp :: (Num a) => (a->a) -> String -> Sym a -> Sym a
unOp f _ (Con c) = Con (f c)
unOp _ "negate" (App "negate" _ [x]) = x
unOp _ "abs" e@(App "abs" _ _) = e
unOp _ "signum" e@(App "signum" _ _) = e
unOp f op x = App op (\ [a] -> f a) [x]

isCon :: Sym a -> Bool
isCon (Con _) = True
isCon _ = False


instance (Real a) => Real (Sym a) where
    toRational (Con c) = toRational c

instance (RealFrac a) => RealFrac (Sym a) where
    properFraction (Con c) = (i, Con c') where (i, c') = properFraction c

instance (Floating a) => Floating (Sym a) where
    pi = var "pi"
    exp = unOp exp "exp"
    sqrt = unOp sqrt "sqrt"
    log = unOp log "log"
    x ** y = binOp (**) x "**" y
    logBase x y = binOp logBase x "logBase" y
    sin = unOp sin "sin"
    tan = unOp tan "tan"
    cos = unOp cos "cos"
    asin = unOp asin "asin"
    atan = unOp atan "atan"
    acos = unOp acos "acos"
    sinh = unOp sinh "sinh"
    tanh = unOp tanh "tanh"
    cosh = unOp cosh "cosh"
    asinh = unOp asinh "asinh"
    atanh = unOp atanh "atanh"
    acosh = unOp acosh "acosh"

instance (RealFloat a) => RealFloat (Sym a) where
    floatRadix = floatRadix . unSym
    floatDigits = floatDigits . unSym
    floatRange  = floatRange . unSym
    decodeFloat (Con c) = decodeFloat c
    encodeFloat m e = Con (encodeFloat m e)
    exponent (Con c) = exponent c
    exponent _ = 0
    significand (Con c) = Con (significand c)
    scaleFloat k (Con c) = Con (scaleFloat k c)
    scaleFloat _ x = x
    isNaN (Con c) = isNaN c
    isInfinite (Con c) = isInfinite c
    isDenormalized (Con c) = isDenormalized c
    isNegativeZero (Con c) = isNegativeZero c
    isIEEE = isIEEE . unSym
    atan2 x y = binOp atan2 x "atan2" y

2007-04-09T15:08:00.000+01:00

Well, it was bound to happen. I have started a blog, just like everyone else.

Things that amuse me

Impredicative polymorphism, a use case

More points for lazy evaluation

Lazy bindings

Lazy functions

Lazy constructors

Cyclic data structures

Reuse

Complexity

Summing up

Ugly memoization

A snag

Another snag

Moral

More LLVM

The UExp module

The TExp module

The Compiler module

Conclusion

More BASIC

Is Haskell fast?

Regression

A performance update

Vectors

LLVM arithmetic

Another example

Black-Scholes

Some figures

LLVM

Haskell LLVM bindings

An example

Fibonacci

Hello, World!

Phi instructions

Loops and arrays

Abstraction

File operations

The OCaml code again

Abstracting on, suggested solutions

Explicit type applications

Something like abstype

The abstraction continues

Another problem

A reflection on type signatures

A little ML

A somewhat failed adventure in Haskell abstraction

The simple problem

Type classes

Associated types

Associated data types

Conclusion

Lost and found

Traced values

Hacking with types

More examples

And what about exponentiation?

Simple reflections of higher order

Some lambda calculus examples

Syntax

Enter let

Bool

A simple top level

Characters

Pairs

Maybe

Natural numbers

Lists

Benchmarking ray tracing, Haskell vs. OCaml

Results

Software and hardware details

Some observations

Conclusion

Simpler, Easier!

First, the untyped lambda calculus.

A detour, simply typed lambda calculus.

Almost going down the wrong track, the polymorphic lambda calculus.

Simply easy, the lambda cube.

Examples

Something like `abstype`

Enter `let`