09 April 2007

A Wonderful Waste of Time

I just encountered this great toy to challenge your three-dimensional geometry skills. The game in effect asks you to pretend you're creating an Escher-like puzzle, producing a an object that only looks right from three precise vantage points. Ever wanted to pretend to be a 3D rendering optimizer function? Now's your chance!

Warning: if you start it, you won't be able to stop until you complete all the puzzles. Which I did, although it took me an embarassingly long time.

Haskell for the Short Attention Span: A Simple File Filter

Well, there is more to say about run-length encoding and other interesting examples I was sent in response to my last article, but today I've been working on another little real-world problem. I've been able to spend a little bit of time on #haskell and the gang there has been extremely helpful. The task: filter a binary log file, turning it into a text representation, and do some stateful validation. I'm not going to show you all the code (it is pretty boring), but here are the minor pitfalls and helpful suggestions I encountered along the way. First, I started with a very basic file filter example from the Haskell 98 report:

main = do
    putStr "Input file: "
    ifile <- getLine 
    putStr "Output file: "
    ofile <- getLine 
    s <- readFile ifile 
    writeFile ofile (filter isAscii s)
    putStr "Filtering successful\n"

Pretty simple. As a sanity check, I wanted to run this example on a data file. The first thing I found is that isAscii is not in the default namespace. To use this program in GHC you'll need to import the module Char. You can learn this via Hoogle. I've started using Hoogle quite a bit; it is a great little tool! You can find Hoogle here. You can also get to it directly through #haskell via LambdaBot. Or you can put it right in your GHCi. Lambdabot lives here.

Anyway, the next thing I found was that, at least under Cygwin, unless you explicitly open a file in binary mode, you might encounter problems. One I should have anticipated -- you might see line feed translation. But there are apparently others -- like silent truncation of the data! Somewhere in the guts of Cygwin, or maybe Windows, control codes designed back in the Jurassic days of computing are still being honored. Input was terminating on (apparently) an EOF character in my binary file. Another suggestion from one of the very smart folks on #haskell.

I'll just use readBinaryFile instead of readFile, right? Well, no, readBinaryFile is part of MissingH, a third-party library. Back to #haskell, where I was advised that I could roll my own. Importing System.IO, I use the following snippet:

readBinaryFile s = System.IO.openBinaryFile s
    System.IO.ReadMode >>= System.IO.hGetContents

This works fine, although I would humbly suggest that GHC's standard library should provide easy access to binary files; while the solution is pretty trivial, it is a wart.

I next stubbed my toe on printf. It is available in module Text.Printf and provides a type-safe version of C's printf. As a long-time C and C++ programmer, you'd think I'd know just about everything there is to know about printf. So, I rapidly gave it a format string "%02X" and passed it a char. Apparently the uppercase X to produce a hex representation with uppercase A-F is not supported (grrr). Another minor wart -- if you provide printf, it should behave like printf -- but we'll move on.

Per more chatting on #haskell I was given this one liner to dump binary data in a nicely formatted way:

writeFile ofile (concat $ zipWith (printf "%02x %s") s (cycle $ replicate 19 "" ++ ["\n"]))

I want to take a moment to talk about how it works. First, cycle $ replicate 19 "" ++ ["\n"]. The replicate function gives us a list of 19 empty strings, which we then concatenate with a newline. Applying cycle to this list treats it as an infinitely repeating circular list of strings, where every twentieth is a newline. These arguments are then fed to printf using zipWith. zipWith is an interesting function: while zip takes two lists and generates a list of pairs produced by assembling the list elements into tuples, zipWith doesn't tupleize the elements; instead it feeds the elements to the provided function, and makes a list of the results.

While this worked, it was interesting enough that I wanted to play with it using GHCi. But I had to give up on that; I kept tripping over the type checker. While I appreciate the masochistic joys of programming and the safety that comes with it, it can be frustrating for programmers with experience in, say, Ruby, or even C.

For example:

let xs = [1..100]
let ys = take 100 (cycle $ replicate 19 "" ++ ["\n"])
zipWith (printf "%02x %s") xs ys

GHC replies:

Ambiguous type variable `c' in the constraint:
`PrintfType c' arising from use of `printf' at <interactive>:1:9-24
Probable fix: add a type signature that fixes these type variable(s)

Ugh. Using the :t command in GHC it is easy to see that GHC thinks the type of xs is [Integer] and ys is [[Char]] (a list of list of chars, also known as a list of strings). If I put roughly the same code in a Literate Haskell source file and ask GHC to load it, I get:

Ambiguous type variable `a' in the constraints:
  `Enum a'
    arising from the arithmetic sequence `1 .. 100'
    at E:\toy.lhs:3:5-12
  `Num a' arising from the literal `100' at E:\toy.lhs:3:9-11
  `PrintfArg a' arising from use of `printf' at E:\toy.lhs:5:18-33
Possible cause: the monomorphism restriction applied to the following:
  xs :: [a] (bound at E:\toy.lhs:3:0)
Probable fix: give these definition(s) an explicit type signature
              or use -fno-monomorphism-restriction

followed immediately by:

Ambiguous type variable `c' in the constraint:
  `PrintfType c' arising from use of `printf' at E:\toy.lhs:5:18-33
Possible cause: the monomorphism restriction applied to the following:
  result :: [c] (bound at E:\toy.lhs:5:0)
Probable fix: give these definition(s) an explicit type signature
              or use -fno-monomorphism-restriction

Failed, modules loaded: none.

Wow. I'd say that is not really a newbie-friendly error message. However, this printf works fine in my real program. I'm not certain why, and I'm not going to dive into it too deeply right now. But here's a simpler type checking example: while C is strongly typed, you can treat numbers as chars and vice-versa, as long as you keep integral promotion and sign extension in mind. GHC is a harsher mistress. Let's say we want to pattern-match on our binary data. The value 16 in my binary data is DLE, which stands for Data Link Escape; it is often used in serial data to indicate packet boundaries, while inside the payload, it will be escaped (doubled). So here's a little pattern to remove doubled DLEs:

de_dle (16:16:xs) = ...

Simple enough, right? No, to Haskell a number and a Char are not interchangeable. Back to #haskell, where I got a quick explanation of the type checker's error messages. I turned the numbers into chars:

de_dle ('\16':'\16':xs) = ...

And that works just beautifully.

Anyway, to make a long story shorter, I was able to write my file filter, which does some nice pattern matching and formatted out. The problems I had while developing that were the kind I like: problems choosing my algorithm properly, not problems fighting with the language. The runtime was quite helpful here; while processing the file, if I hit a case at runtime which my patterns did not handle, I got a runtime warning about non-exhaustive patterns. That led me to reorganize my patterns, and the result was much clearer.

There was one more minor pitfall remaining. I compiled my program using GHC, but when I ran it, instead of my prompts for input and output filenames, I got nothing! Haskell was silently waiting for input. Back to #haskell. To make the output show up, I had to import System.IO and do hSetBuffering stdout NoBuffering. (Alternately, I could flush stdout immediately after each putStr, but that seems even uglier). I hope that saves someone a little aggravation.

Speaking of aggravation, how did it all come out? Well, the original binary log file is about six megabytes. Since I was doing this by hand, I was more concerned that the program ran correctly than that it ran fast; I would have been satisfied with anything under a half-hour. In fact, without any attempts at optimization at all, my filter ran in under thirty seconds, which is more than fast enough for an ad hoc little tool. I started out trying to do this task using some regular expressions in vi, and that was quickly going nowhere, and taking forever to do it. The vi that came with Cygwin doesn't support some of the more advanced regular expressions features (there is no {x,y} syntax for specifying the number of repeats of a pattern). Notepad++, my workhorse Windows text editor, also doesn't support this syntax. Without this the regexes were becoming hideous, although in (say) Perl they would have been rather simple. But I was able to use Haskell instead, thanks to GHC and the kind folks on #haskell!

07 March 2007

Haskell for the Short Attention Span: Run-Length Encoding, Part 3

This isn't a big improvement, but it occurs to me that my implementation from last time:

>deRLE :: [Integer] -> Bool -> Integer -> [Integer]
>deRLE [] False _ = []
>deRLE points False max = deRLEHelper points max
>deRLE points True max = deRLEHelper (0 : points) max

>deRLEHelper :: [Integer] -> Integer -> [Integer]
>deRLEHelper [] max = []
>deRLEHelper [unmatched] max = [unmatched..max]
>deRLEHelper (start:end:rest) max = [start..end - 1]
>    ++ deRLEHelper rest max

can be yet shorter; since deRLEHelper handles the empty list specially, deRLE doesn't need to:

>deRLE :: [Integer] -> Bool -> Integer -> [Integer]
>deRLE points False max = deRLEHelper points max
>deRLE points True max = deRLEHelper (0 : points) max

That way only the function that actually recurses down to the empty list has to explicitly handle it.

Can we make it even more concise and still retain the self-documenting aspects of the code? I'll have to come back to that question. It's time to respond to comments. When I first mentioned my triplify function, which breaks a string into strings of three characters, I wrote "I'm sure there is a very clever one-line fold that applies take 3 in an mind-expanding manner, but I was unable to find it." I got some great suggestions. The first, due to Cale Gibbard, was:

map (take 3) . takeWhile (not . null) . iterate (drop 3)

Let's look at it piecewise as a pipeline, from right to left.

First, our list/string (for example, "ABCDEFG") is fed to the function iterate along with (drop 3). Iterate is an interesting function. I tried to understand it a few months ago by reading the description at zvon.org, which says that iterate

creates an infinite list where the first item is calculated by applying the function on the second argument, the second item by applying the function on the previous result and so on.

They give the following example:

Input: take 10 (iterate (2*) 1)
Output: [1,2,4,8,16,32,64,128,256,512]

Note that if the behavior followed the description, the output would be [2,4,8,16,32,64,128,256,512,1024]. Their description is wrong, although the output matches GHCi's output. What does iterate really do? I'll describe it in the form of a poem.

iterate takes a function and a value, and returns a list consisting of:

the value

followed by
  the value of applying the function to the value

followed by
  the value of applying the function to
    the value of applying the function to the value

followed by
  the value of applying the function to
    the value of applying the function to
      the value of applying the function to the value
      
it will keep talking
  even when
    it has nothing
      new to say
        so please don't
        ask it any
          open-ended
            questions...

That is, I can't evaluate

iterate (drop 3) "ABCDEFG"

because even though the input is finite, the output isn't; there's no termination to the recursion. However, I can evaluate part of the result:

take 4 (iterate (drop 3) "ABCDEFG")
["ABCDEFG","DEFG","G",""]

And because of Haskell's lazy evaluation, we can feed this infinite result "upstream" to be used in the remainder of our calculation, and it will work fine as long as it doesn't evaluate the complete result. The next stage in our "points-free" pipeline is:

takeWhile (not . null)

The takeWhile function operates on a list, returning a list made up of its values as long as they match the predicate supplied. This has the effect of squashing the overly talkative iterate, because we stop evaluating its results as soon as it returns the first null value. Lazy evaluation is so cool!

Now, our finite list of strings is handed off to map (take 3), which generates a list of 3-character prefixes. Here's the whole process:

iterate (drop 3) "ABCDEFG"

yields an infinite list beginning ["ABCDEFG","DEFG","G",""];

takeWhile (not . null)

applied to the above will take elements from the list until it hits the first null: ["ABCDEFG","DEFG","G"]

map (take 3)

applied to the above will give us ["ABC","DEF","G"]. Note that the short remainder is preserved.

>

That seems so useful, I'm surprised it isn't a function in the Standard Prelude, called groupsOf. Although it would also be nice to have a version which didn't keep the short remainder, perhaps groupsOfOnly.

I got another suggestion from a user called "kirby." He suggested the function

takeWhile (\="") . List . unfoldr (Just . splitAt 3)

Wow, my first unfoldr! I'm going to have to stop there for today, while I learn a little something about Just, Nothing, and the Maybe warm fuzzy thing. My first funster!

06 March 2007

English Majors as Programmers

The Embedded Muse is an e-newsletter edited by Jack Ganssle; see the Ganssle Group home page. In Issue 142 he writes:

I've found that some of the best developers of all are English majors. They'll often graduate with no programming experience at all, and certainly without a clue about the difference between DRAM and EPROM. But they can write. That's the art of conveying information concisely and clearly. Software development and writing are both the art of knowing what you're going to do, and then lucidly expressing your ideas. The worst developers, regardless of background, fail due to their inability to be clear. Their thoughts and code tend to ramble rather than zero-in on the goal... Too many engineering-trained developers have a total disregard for stylistic issues in programming. Anything goes. Firmware is the most expensive thing in the universe, so it makes sense to craft it carefully and in accordance with a standard style guide. Make sure it clearly communicates its intent. This is where the English majors shine; they've spent 4 years learning everything there is to know about styles and communication.

As an English major who has programmed computers since the age of ten, and who also took every programming course I could manage, it's nice to hear a little love for English majors! But I should add that an English degree is not a requirement; I also know some excellent developers who majored in more technical subjects, such as anthropology and archaeology!

05 March 2007

Haskell for the Short Attention Span: Run-Length Encoding, Part 2

Greetings! Last time I introduced a piece of code to solve a small real-world problem: decoding a set of function IDs from a run-length-encoded list of change points. I received some excellent feedback. One of the most valuable suggestions was that I avoid unnecessarily taking the length of the list in my "triplify" function. That makes a great deal of sense, considering that ideally I'd like to be able to handle infinite lists, and of course it would be nice if the runtime of my algorithm did not spiral out of control quite so quickly!

I also gave a little more thought to the various base cases. One of the nicest things about Haskell is its pattern matching. Indeed, even a simple function like "square a = a * a" makes use of pattern matching, where "a" on the left-hand side is irrefutable. My implementation last time failed to fully take advantage of pattern matching on multiple parameters, particularly to handle base cases, and did not even handle all the various termination conditions properly. So, here is a slightly condensed version of the code that makes the termination conditions much more explicit, and so I claim even though it is shorter, it is more self-explanatory.

>import Char
>import Numeric
>rle_raw_str = "00 30 90 09 31 01 10 32 20 22 12 30 23 12 41 24" ++
>  "34 40 44 1C 0F C1 1D 00 D0 1D 03 D0 4D 05 D0 7D 0A D0 CD 10" ++
>    "D1 2D 20 D2 6D 30 D3 1D 40 D4 1D 50 D5 1F 00"

A new triplify:

>
>triplify :: String -> [String]
>triplify [] = []
>triplify (x:y:z:rest) = [x, y, z] : triplify rest
>triplify str = error "Unexpected string length (trying to match 3 characters)"

Treating the input data as a list of 12-bit numbers:

>rle_change_points = map (fst . head . readHex)
>  (triplify (filter isHexDigit rle_raw_str))

The main function: given a list of change points, a hint as to whether we are initially presuming the decoder will output ones or zeroes, and a maximum value, decode RLE data from the change points. Return the indices of the one bits.

>deRLE :: [Integer] -> Bool -> Integer -> [Integer]
>deRLE [] False _ = []
>deRLE points False max = deRLEHelper points max
>deRLE points True max = deRLEHelper (0 : points) max

>deRLEHelper :: [Integer] -> Integer -> [Integer]
>deRLEHelper [] max = []
>deRLEHelper [unmatched] max = [unmatched..max]
>deRLEHelper (start:end:rest) max = [start..end - 1] 
>  ++ deRLEHelper rest max

This is more like it -- the kind of program that resembles a proof more than code. I'm starting to prefer this style!

A brief explanation: if we're initially decoding zeroes, the empty list case can be disposed of quickly. The program filters this case out first so that the recursive helper function, which accumulates the output, doesn't need to handle this case. The next two lines handle non-empty lists. If we're initially decoding ones, we can still treat all cases the same if we just prepend a zero to "flip" the decoder. Our helper function then has only three cases: a termination case where an interval of ones is not closed (because we have a single change point left), in which case we spit out ones up to and including our maximum possible value; a termination case where we have two change points forming a closed interval, in which case we spit out an interval inclusive of the start and exclusive of the end, and a non-termination case in which we have more than two entries, which spits out a closed interval and continues the process.

Let's look at some test cases. First, we confirm that if we're generating zeroes and we encounter no change points, our result set is empty:

deRLE [] False 0xF

[]

We confirm that if we hit a single change point, the result set will contain all the values up to and including the maximum possible:

deRLE [0] False 0xF

[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15]

deRLE [11] False 0xF

[11,12,13,14,15]

In the general case where our interval is closed, our result set is inclusive on the left and exclusive on the right:

deRLE [5, 10] False 0xF

[5, 6, 7, 8, 9]

Closing one range and opening another without closing it works correctly:

deRLE [5, 10, 11] False 0xF

[5,6,7,8,9,11,12,13,14,15]

If we begin assuming our decoder generates zeroes, an empty list of change points maps to all possible values:

deRLE [] True 0xF

[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15]

A single change point at zero toggles the state of the decoder and produces and empty result:

deRLE [0] True 0xF

[]

A single change point greater than zero allows some values into our result set before it is closed:

deRLE [5] True 0xF

[0,1,2,3,4]

Two change points defines a range to exclude from the results, and we get the rest of the values in the range:

deRLE [5, 10] True 0xF

[0,1,2,3,4,10,11,12,13,14,15]

Three change points suppresses the rest of the result set:

deRLE [5, 10, 13]

[0,1,2,3,4,10,11,12]

With these tests I'm satisified that deRLE works for all expected cases.

I received more comments, including some recommendations for a nifty unfold. Thanks to everyone who made suggestions. Next time I'll take a look at how those work, talk about anything else that comes to mind, and then proceed to my next real-world problem!

01 March 2007

Haskell for the Short Attention Span: Run Length Encoding, Part 1

I'm a fortunate man. Some 50% of people who suffer herpes zoster infections involving the eye wind up with permanent eye damage. I seem to have escaped with little or no permanent damage. The nerves are not completely healed, and there is still visible damage under the skin, but it is no longer painful to work at the computer or read!

It's time to jump back on the horse and risk looking foolish in public by writing some code. I still consider myself to be close to a rank beginner in Haskell, so you're reading my "finger exercises." To experienced functional programmers they will probably sound like a new violin student sawing away. However, since I haven't been a student for almost twenty years and there is no one to assign me homework but myself, this seems to be the next best way to force myself to learn something!

With three children at home, two of them babies, I don't get a lot of unbroken free time; I've got a guaranteed short attention span. So consequently I'm going to look at some very small programs. However, rather than make up toy programs from scratch, I am now going to turn to a few real-world problems that I've stumbled across while doing embedded programming.

OK. Consider a spherical cow of uniform density... wait, I mean, consider a black-box system with an outward-facing interface in which a number of functions are made available. Functions are accessed via a 12-bit unique ID. A client of that system doesn't know which functions are available. The client needs a method to retrieve a list of funciton IDs. Let's assume we don't have bandwidth to burn, so we don't want to send a complete list of all the function IDs. Let's further assume that the functions are organized into a small number of groups with consecutive IDs. In other words, this is data that is very amenable to a simple compression scheme.

We're going to do it using Run-Length Encoding (RLE). Here is a brief informal explanation of RLE; if you are already know all about RLE, skip the next paragraph (but read the ones after that, please!)

RLE is a commonly used technique for encoding bitmaps; a variant of RLE is even used in fax machines. Rather than store a bit for every bit in a bitmap, RLE records _runs_ of repeating data. If you have a raster image that you're scanning left to right, top to bottom, and that image has a lot of horizontal lines in it, or black and white rectangles, this might be a big win. The RLE encoding will record long alternating runs of ones and zeroes. If the image has a lot of vertical lines, this might not be much of a win. The worst-case scenario is a 50% checkerboard: in that case, we'll have to indicate runs of length 1, and the overhead will be far worse than just encoding the bits. For that kind of image you'd be better off using a means of compression that actually adapts to the frequency of occurrence of particular patterns in the data, such as Huffman coding or LZW.

Now, there are several possible ways to represent the runs of RLE data. One scheme might use tuples consisting of the bit value and the number of occurrences of that value (the "run length"). For example, the tuples [(0, 14), (1, 28), (0, 3)] could represent a run of fifteen zeroes, followed by a run of twenty-nine ones, followed by four zeroes. Why are the values off by one? Because we would never need to encode a run of length zero, so we allow zero to represent one, and so on. A tuple like this could be encoded in a single byte, where the high-order bit represents the value and the remaining seven bits represented the run length - 1. This gives us the ability to represent runs of up to 128 bits. Of course, if we have a run of length 129 or greater, we'll need to represent it with multiple bytes. For degenerate cases where we encode very long runs, we'll have unnecessary redundancy.

For that low-entropy bitmap, a better scheme might be to just record the points at which the data _changes_. We have to choose an initial state for decoding: should the decoder, prior to the first change point, produce ones, or zeroes? Once we know that, we don't need to encode the bit values at all. If we assume that our encoded data begins with a run of zeroes, we could represent the example above (fifteen zeroes, followed by twenty-nine ones, followed by four zeroes) by recording only [15, 44]. Let's say that we want our decoder to spit out the indices of the one bits. We'll get:
[15,16,17,18,19,20,21,22,23,24,
25,26,27,28,29,30,31,32,33,34,
35,36,37,38,39,40,41,42,43]
Using this scheme, a bitmap of all zeroes would be represented by an empty list, while a bitmap of all ones would be represented by by a single [0]; it turns on the "ones" state and never turns it off, so we just generate ones until our range of values is exhausted. (This implies that, in addition to knowing what starting state it should assume, the decoder has to know the size of the bitmap we are representing). A bitmap consisting of all zeroes followed by a single one would be represented by [N - 1] where N is the number of bits in the bitmap.

If, instead of assuming that our data starts with an implicit run of zeroes, we assume instead that it starts with an implicit run of ones, it should be obvious that we could encode the same example using [0, 15, 44]; in this case, a bitmap of all zeroes would be represented by [0], which would yield an empty list of one-bit indices, while a bitmap of all ones would be represented by an empty list, yielding a list [0..n - 1] of bit indices.

We'll assume that our function IDs are represented this way: they are the one bits in the bitmap, and we want a list of their indices. When we query our interface about its function IDs, we get back a response that consists of a short list of change points. Here are the data bytes:
00 30 90 09 31 01 10 32 20 22 12 30 23 12 41 24 34 
40 44 1C 0F C1 1D 00 D0 1D 03 D0 4D 05 D0 7D 0A D0
CD 10 D1 2D 20 D2 6D 30 D3 1D 40 D4 1D 50 D5 1F 00
The first time I encountered data like this, I wrote myself a little bit of Ruby code to decode it. Here's first part of the Ruby code, which turns these bytes into a series of 12-bit unsigned numbers by assembling three nybbles (half-bytes) at a time:
rle_raw_list = gets.chomp.delete(" ")

rle_change_func_ids = []
0.step(rle_raw_list.length - 2, 3) { |nybble_idx|
rle_change_func_ids.push(rle_raw_list[nybble_idx..nybble_idx + 2].hex)
}
Here's how it works. The first line reads data from standard input and, reading left to right, removes any end-of-line character and deletes spaces. We next create an empty list. "0.step" is cute little bit of idiomatic Ruby shorthand: everything is an object, including integers, and they have a "step" method, so by calling 0.step we loop from zero to the string length - 2, incrementing by threes. We also provide a Ruby block, which is basically a closure that receives one parameter, the loop index, and which is bound to the name nybble_idx. Working inside-out, the body of this closure takes three-byte substrings of the raw data, passes them to the hex function to yield an integer, and then uses push, which treats the list like a stack and appends the new integer to the end. The result looks like this:
[0x003, 0x090, 0x093, 0x101, 0x103 ...]
This is the list of change points. Once I have that, I can generate the runs of data. I'll assume that our data starts with an initial run of ones:
rle_state = true
MAX_FUNC_ID = 0xFFF
current_func_id = 0

for change_func_id in rle_change_func_ids
if rle_state == true
(current_func_id...change_func_id).each {|valid_func_id|
puts sprintf('0x%02X', valid_func_id)
}
rle_state = false
elsif rle_state == false
rle_state = true
end
current_func_id = change_func_id
end
Now, I should mention that I like Ruby, and especially like its support for nice idioms like blocks and closures and iterators and chaining functions. But while this Ruby code worked well enough for my little task, you might notice that it doesn't actually handle all the termination cases properly. As I work with Haskell, I'm finding that I like it quite a bit more than Ruby. One of the reasons is that Haskell code, being less imperative, seems closer to a kind of abstract, platonic representation of the problem space. Problems such as the failure to properly handle various termination cases tend to look much more obvious, particularly when functions are broken down into pieces using pattern matching.

Anyway, here's my Haskell version:
>import Char
>import Numeric

>rle_raw_str = "00 30 90 09 31 01 10 32 20 22 12 30 23 12 41 24" ++
> "34 40 44 1C 0F C1 1D 00 D0 1D 03 D0 4D 05 D0 7D 0A D0 CD 10" ++
> "D1 2D 20 D2 6D 30 D3 1D 40 D4 1D 50 D5 1F 00"
(My first opportunity to look silly: I'm sure there's a way to break a string across multiple lines with some sort of continuation character, but I could not get "\" to work together with Literate Haskell mode, so I'll just concatenate them and hope no one notices...)
>rle_bytes_str = filter isHexDigit rle_raw_str
That was easy!

Now, my second opportunity to look silly: I'm sure there is a very clever one-line fold that applies "take 3" in an mind-expanding manner, but I was unable to find it. This was about as clever as I could manage, but it does have the advantage of checking for unexpected data length:
>triplify :: String -> [String]
>triplify str | length str == 0 = []
>triplify str | length str `mod` 3 == 0 = ( take 3 str ) : triplify ( drop 3 str )
>triplify str = error "bad length (expected multiple of three)"
To see the results, evaluate "triplify rle_bytes_str." I get:
["003","090","093","101","103","220","221","230","231","241",
"243","440","441","C0F","C11","D00","D01","D03","D04","D05",
"D07","D0A","D0C","D10","D12","D20","D26","D30","D31","D40",
"D41","D50","D51","F00"]
Now I want to turn these into a list of numbers; these are the change points. Mumble, mumble...
>rle_change_points = map (fst . head . readHex) (triplify rle_bytes_str)

[3,144,147,257,259,544,545,560,561,577,579,1088,1089,3087,3089,3328,3329,3331,3332,
3333,3335,3338,3340,3344,3346,3360,3366,3376,3377,3392,3393,3408,3409,3840]
Our playing field is the range of integers 0x000..0xFFF or 0..4095 in decimal. We want to lazily generate a list of valid values using our RLE decoding scheme. Rather than walk the list of change points spitting out values using some kind of toggling state variable to keep track of whether the decoder is processing a run of ones or a run of zeroes, or walking all the possible values filtering it by the list of change points, we'll use two functions, one for extracting ranges, and one for generating values from ranges. Both functions receive a list representing "the rest of the work to do" in the form of the remainder of the list of change points, and so this is a _little_ bit like using continuation-passing style, or using a Monad to maintain the progress of the calculation. Maybe in Part 2 I'll have learned enough to explain what I mean by that!

Our first function, getRangeFromChangePoints, has to handle some interesting cases. Assuming we are decoding RLE in the range of possible values [X..Y], what happens when our list contains a final value, N, that opens a final range but doesn't explicitly close it? Well, there are two approaches we could take. The more conservative approach is to assume that this was a harmless mistake and close the range [N..N + 1) so that it contains only the value N. (Note that I'm using interval notation in which a square bracket indicates that the range _includes_ the bound specified, and a parenthesis _excludes_ the bound specified). However, there is a problem with this approach: it means we can't generate a range that includes the maximum value, Y! That's no good, so we'll assume that a final range opening with N means that we really want [N..Y] (inclusive). In fact, since we treat ranges as exclusive at the end, this is actually the only way we can generate a range containing Y.

For now, we also assume that running out of values is an error. We're assuming that our caller will handle the termination conditions.
>getRangeFromChangePoints :: [Integer] -> (Integer, Integer, [Integer])
>getRangeFromChangePoints ([]) = error "getRangeFromRLE: no change points in list!"
>getRangeFromChangePoints (a:[]) = (a, max, []) where max = 0xFFF
>getRangeFromChangePoints (a:b:[]) = (a, b - 1, [])
>getRangeFromChangePoints (a:b:ls) = (a, b - 1, ls)
Wow, that's pretty wordy as Haskell code goes, but I kind of like the way the pattern-matching breaks down. Note that in particular it makes the termination conditions obvious. This helps me think about the problem domain, rather than the implementation, in a way that the imperative Ruby code did not.

Next, here is the recursive function to generate values from ranges. We represent a range using the first two values of our tuple, while the third value is a list of the remaining ranges (sort of like our continuation). We handle the termination case, in which we have no more work to do, by just returning a list from our range, while the non-terminating case gets the next range (including the list representing our remaining work) and calls the function again to continue generating the list.
>genValuesFromRange :: (Integer, Integer, [Integer]) -> [Integer]
>genValuesFromRange (first, last, []) = [first..last]
>genValuesFromRange (first, last, to_do) = [first..last] ++
> (genValuesFromRange (getRangeFromChangePoints to_do))
To generate the initial case we use a helper function called valuesFromRLE. In ddition to the change points list, it takes a Boolean value which indicates whether we are to assume that our RLE decoding is initially inclusive or initially exclusive: that is, whether the first value in our first range indicates that the range is starting or ending. If we want it to behave inclusively, we prepend a zero value to create a range starting with zero. Note that this will be harmless if our first range starts with zero, because we'll generate a range like [0..(-1)], which will contain no members.
>valuesFromRLE :: [Integer] -> Bool -> [Integer]
>valuesFromRLE change_points initially_inclusive = genValuesFromRange $
> if initially_inclusive then getRangeFromChangePoints (0 : change_points)
> else getRangeFromChangePoints change_points
Note that we can evaluate the whole thing, but more interestingly, we can use "take" to produce only the first few elements:
take 16 $ valuesFromRLE rle_change_points True

[0,1,2,144,145,146,257,258,544,560,577,578,1088,3087,3088,3328]
And there you have it. If you like, play with the code. Next time we'll see if we can encode the function IDs back into change points, and confirm the round trip. As always, I appreciate your comments.

15 February 2007

Richard or Edna? (And Graham Hutton's Programming in Haskell)

For some reason, Richard Bird's Introduction to Functional Programming Using Haskell, 2nd Edition, is now listed on Amazon as written by "Edna Bird" (and coauthored by "Wadler" with no first name). Assuming Richard is not actually undergoing a transformation along the lines of Walter Carlos to Wendy, this looks like an error...

Anyway, my copy of Graham Hutton's Programming in Haskell arrived from Amazon. The first remarkable thing about this book is that it is the thinnest book on Haskell I've yet seen. The highest page number is 170. It is even thinner than The C Programming Language, Second Edition.

What does this mean? Well, for one thing, it means the chapters are very brief. It is a beginner's text aimed at novices, even those who may have not programmed (in any language). I would definitely consider using it as a textbook for an introductory programming class. (I have a dream about teaching a programming class to gifted high school students or college first-year students before they have had their minds corrupted by too much Java or C++; we'd learn a little Haskell, Scheme, and Dylan, in that order, and then use our perspective gained to take a brief look at Java, Pascal, and BASIC. In that order, implementing a parser for Pascal and a complete implementation of BASIC. We'd start with this book. But anyway).

The explanation of curried functions (or how "all functions with multiple arguments are curried") is the clearest I've seen. So is the explanation of basic types and classes. His chapter on recursion is also a model of clarity. (I'm pleased to see it goes above and beyond Fibonacci numbers). In general, Hutton seems to be encouraging his readers to think like programmers think, but he doesn't require very many words to get his ideas across.

If you are busy like I am, you don't have time to read a long review, either, so I'll just conclude by saying that this looks like the introductory Haskell text!

08 February 2007

Graham Hutton's Book has Shipped

Amazon tells me that my copy of Graham Hutton's Programming in Haskell has shipped. Excellent!

With the passage of time and the assistance of an acupuncturist, the nerve damage and eye sensitivity left in the wake of my shingles infection has healed, to the point where reading is no longer quite so painful. And I don't have to wear my sunglasses to work on the computer anymore!

So, soon, hopefully back to actually writing some Haskell, as opposed to just reading about it. I have a couple of small real-world programs to try porting from Ruby, which should be educational. More soon!

17 January 2007

Haskell Books Piling Up, and Herpes Zoster Ophthalmicus at One Month

I've received a couple of Haskell-related books -- my Christmas presents to myself. I have only been able to skim these books, but here are my first impressions.

The first book is Chris Okasaki's Purely Functional Data Structures. When I was working on some code to solve Sudoku puzzles in Scheme, I became stymied because in the midst of my efforts to program in a function style, I kept slipping back into an imperative style. This book is the answer to the question "how do you manage data structures when you don't have mutable variables?" The answer, of course, is that you return new data structures based on the old ones. This is not necessarily as inefficient as you might expect, because the new data structure can frequently share content -- sometimes most of the content -- with the previous data structure. Okasaki's book is a model of clarity; on the back it indicates that it can be used for self-study. A lot of texts claim that, but this book seems to actually deliver it. The diagrams are marvelously clear and the learning curve is gentle. The body of the text uses Standard ML with extensions but there is an appendix containing Haskell code. I will have more to say about this book in future postings; I am really looking forward to trying out some of the ideas it contains.

The second book is The Haskell Road to Logic, Maths, and Programming by Kees Doets and Jan van Eijck. This is part of the King's College Texts in Computing series which includes An Introduction to Lambda Calculi for Computer Scientists. That book is thin, while this one is fairly fat; this volume is bursting with exposition (and code) on logic, set theory, and other topics dear to my heart. It looks quite fascinating and again I hope to have more to say about it soon.

It has been approximately a month since I was hit with the first symptoms of herpes zoster ophthalmicus (shingles infection involving the eye). The news is generally very good: there is little long-term damage to my cornea. My opthalmologist can see some slight scarring, but it should continue to improve. But the pain and hyper-sensitivity to light has persisted, and I have found it difficult to return to work because my eye tends to start watering madly or lose focus when I work at the computer for longer than a few minutes. The problem is in the branches of the trigeminal nerve, not the eyeball itself. Nerves heal slowly, so this is unfortunately not an illness for the impatient!

02 January 2007

Shingles

Well, I had great plans for my ten days off -- work on the apartment, organize the office, deep-clean the kitchen, read, and study Haskell!

So much for best-laid plans. Just before Christmas I came down with a shingles infection, which started in my right eye.

Shingles is a reactivation of the herpes zoster or "varicella" virus, the same virus that causes chicken pox. I did a near all-nighter the Monday a week before Christmas, which might not have been too bad, but then the following night we had a sick baby and I was also unable to sleep. That was apparently enough stress to cause the virus to reactivate.

Shingles is painful. I have become intimately familiar with two new types of pain, which I call the blowtorch and the nightstick. The blowtorch is a feeling like my skin is being burned by an open flame and is bubbling. That is the sensation of the blisters forming, although it comes at other random times too. The nightstick is deep pain along the nerves of my face, where they feel like I've been beaten. That produces a tremendous amount of muscle tension in my jaw, and around the base of my skull as well, where large areas are still hyper-sensitive to touch.

Valtrex (an anti-retroviral drug) also has nasty side effects. I got waves of fever and chills, often back-to-back. It also caused nausea, diziness, and strange visions. I'm off the Valtrex now, so thankfully those side effects have gone away.

For most of my ten days off work I was fairly incapacitated: unable to see, unable to read, lying in bed with fever or chills and with strange pains, covered with blisters and scabs. When I got out of bed, I was dizzy, stumbling around the house looking like a crazed sunburned pirate with an eye-patch on.

Anyway, I am now back at work. The scabs are gone, although my face is still visibly splotchy. My right eye is still swollen and there is still some pain. I am still using an antibiotic gel in my eye, which means my right eye is still hazy, so I am not driving yet. I see my regular doctor today and the ophthalmologist again tomorrow. I think the likely outcome is that my vision will not be significantly damaged. It seems like the pain will persist at reduced level for at least a while longer.

Dear reader, I sincerely hope you had a better holiday than mine. And if not, God bless you!

22 December 2006

Happy Life Day and Wishing You Lots of Haskell in 2007

Let's remember the words of Princess Leia, in the Star Wars Holiday Special:

We celebrate a day of peace. A day of harmony. A day of joy we can all share together joyously. A day that takes us through the darkness. A day that leads us into might. A day that makes us want to celebrate the light. A day that brings the promise that one day, we'll be free to live, to laugh, to dream, to grow, to trust, to love, to be.

Happy Life Day, everyone! Or whatever tradition you celebrate.

I am not planning to work on any more entries until early January, so I just wanted say that I am extremely grateful for all the comments and suggestions I've received on my blog entries. With two babies at home I don't get much time to study textbooks or tutorials, although I try, and I don't get much dedicated coding time. But the suggestions and concrete examples have helped me make good use of what free time I do have. Thanks to all of you for your support learning this exciting language!

19 December 2006

From Bits to Cells: Simple Cellular Automata in Haskell, Part Two

Let's get back to our CA generator. Literate Haskell follows:

Last time we defined a function to generate the next state of a given
cell in our cellular universe, given a rule number and a tuple consisting
of the current state of the cell to the left, the cell itself, and the
cell to the right.

>import Data.Bits

>genNextBit :: Int -> ( Bool, Bool, Bool ) -> Bool
>genNextBit rulenum ( left, center, right ) = rulenum `testBit` idx
>    where idx = ( if left then (4::Int) else (0::Int) ) .|.
>           ( if center then (2::Int) else (0::Int) ) .|.
>           ( if right then (1::Int) else (0::Int) )

Recall that we can use automatic currying to make a rule-applying
function like so:

>rule_30 = genNextBit 30

We can ask GHCi for the type:

:type rule_30

rule_30 :: (Bool, Bool, Bool) -> Bool

I've put it off while I work on the rules, but it is time to figure out
how to actually represent our CA universe. Let's start by using a list.
I know that I'm going to write a number of inefficient functions, and
do evil things like take the length of lists a lot, but let's suspend all
concerns about efficiency over to a future discussion and consider this
purely a proof-of-concept.

Our inital universe at time zero has one cell set to True:

>initial_universe = [True]

But that isn't quite the right representation for the universe, because
it implies that our whole universe is one cell in size. We can't even
apply our rule once because there is no left cell and right cell!
Really, we want to pretend that we have an _infinite_ universe; at
time zero, all the cells to the left and right hold False. Remember,
Haskell is so powerful that it can traverse an infinite list in only
0.0003 seconds! Well, if you don't evaluate the whole thing, that is.
Taking advantage of lazy evaluation, you can define all kinds of
infinite structures. This construct will give us an infinite list of
False values:

>allFalse :: [Bool]
>allFalse = False : allFalse

We don't want to evaluate allFalse, but we can partially evaluate it
using a function like take. So can we represent our universe like this?

>genInfiniteUniverse :: [Bool] -> [Bool]
>genInfiniteUniverse known_universe = allFalse ++ known_universe ++ allFalse

Let's try it:

take 10 ( genInfiniteUniverse initial_universe )

[False,False,False,False,False,False,False,False,False,False]

Nope! Since the left-hand side of the universe is infinite, we will
never reach the middle element, no matter how far we travel from the
start of the list!

That's no good. However, we can do it another way. We'll allow our
universe to be expanded on demand on the left and right sides:

>expandUniverse :: Int -> [Bool] -> [Bool]
>expandUniverse expand_by known_universe =
>   ( take expand_by allFalse ) ++ known_universe ++ ( take expand_by allFalse )

expandUniverse 3 initial_universe

[False,False,False,True,False,False,False]

We can use the expandUniverse function to expand our initial universe
out to a standardized width before we start applying the rules.

First, here's a routine to stringify a universe for display:

>boolToChar :: Bool -> Char
>boolToChar True = '#'
>boolToChar False = ' '

>stringifyUniverse :: [Bool] -> String
>stringifyUniverse ls = map boolToChar ls

Now our infrastructure is in place, so let's figure out how to apply
our generator rule. We know that we want to start with our initial
universe. Let's expand it to a fixed size. This will give us enough
elements to start making left/center/right tuples out of each consecutive
set of three cells. Each tuple is then used to look up the next state
of the cell at the center; this will become an element in our next
universe. Then we move to the next cell (not three cells down). This
means that the tuples overlap.

Let's make the tuples. We have to do some thinking here and consider
all the cases; the behavior isn't immediately obvious. The following
almost works:

universeToTuples :: [Bool] -> [(Bool, Bool, Bool)]
universeToTuples universe | length universe >= 3 =
   ( universe !! 0, universe !! 1, universe !! 2 ) :
   universeToTuples ( tail universe )
universeToTuples universe = []

universeToTuples [False, True, True, True, False]

[(False,True,True),
(True,True,True),
(True,True,False)]

But it isn't quite right; it leaves off the end cases; when we apply
our rules, the intermediate representation of the universe as a list
of tuples to look up cell mappings will shrink. We actually want the
following tuples:

[(False,False,True),
(False,True,True),
(True,True,True),
(True,True,False),
(True,False,False)]

where the first element of the list is considered as if it was just
to the right of an implied False, and the last element is considered
as if it was just to the left of another implied False. This sounds
like another place we can use our universe expander:

>universeToTuples :: [Bool] -> [(Bool, Bool, Bool)]
>universeToTuples [] = []
>universeToTuples universe = tupleize $ expandUniverse 1 universe
>   where
>       tupleize xs =
>            if length xs > 3 then tuple3 xs : tupleize ( tail xs )
>               else [ tuple3 xs ]
>       tuple3 xs = ( xs !! 0, xs !! 1, xs !! 2 )

Why did I write it that way? Well, I tried to write tupleize using
guards, special-casing length xs > 3 followed by an unguarded case for
all other possibilities, but GHC didn't like it -- it told me I had non-exhaustive patterns. There is probably a smarter way to write this,
but note that we definitely don't want this version:

universeToTuples universe = ( xs !! 0, xs !! 1, xs !! 2 )
    : universeToTuples ( tail xs )
       where xs = expandUniverse 1 universe

In that version, the universe keeps expanding as we walk down the list,
and we never get to the end!

OK, now that we have our tuples, we want to turn them into our new
universe, given a cellular rule number:

>tuplesToUniverse :: Int -> [(Bool, Bool, Bool)] -> [Bool]
>tuplesToUniverse rule [] = []
>tuplesToUniverse rule (tup:tups) = genNextBit rule tup : tuplesToUniverse rule tups

Note that we don't have to explicitly take the tail since we provide a
name for it in the pattern. We're ready to define our genUniverses function
that applies a given CA rule. We can express a given generation like this:

>nextUniverse :: Int -> [Bool] -> [Bool]
>nextUniverse rule universe = tuplesToUniverse rule $ universeToTuples universe

Now, let's generalize it:

>genUniverses :: Int -> Int -> Int -> [[Bool]]
>genUniverses rule width count = take count
>   ( iterate ( nextUniverse rule ) ( expandUniverse ( width `div` 2 ) initial_universe ) )

(You could also use a fold, and I'm sure there are lots of other ways to
do it, but iterate seems to work fine).

And now, witness the unfolding of a universe! Note that the parameters
that go to genUniverses are the rule number, the width for display, and the number of generations:

putStr $ unlines $ map stringifyUniverse $ genUniverses 222 19 10

          #
         ###
        #####
       #######
      #########
     ###########
    #############
   ###############
  #################
 ###################

In general, a width of twice the number of generations - 1 will show
all the transitions we are interested in; you could consider the diagonal
area above to be the "light cone" of events causally connected to that
single point (although some rules will generate True cells outside of
that "light cone" based on the other False values in the initial universe). Let's make a helper function to choose a width for us:

>showRule rule gens =
>   putStr $ unlines $ map stringifyUniverse $
>       genUniverses rule ( gens * 2 - 1 ) gens

Let's try a few of the other rules:

showRule 252 15
              #
              ##
              ###
              ####
              #####
              ######
              #######
              ########
              #########
              ##########
              ###########
              ############
              #############
              ##############
              ###############

showRule 78 15
               #
              ##
             ###
            ## #
           ### #
          ## # #
         ### # #
        ## # # #
       ### # # #
      ## # # # #
     ### # # # #
    ## # # # # #
   ### # # # # #
  ## # # # # # #
 ### # # # # # #

And finally, my all-time favorite, which simulates a Sierpinski gasket:

showRule 82 32
                                #
                               # #
                              #   #
                             # # # #
                            #       #
                           # #     # #
                          #   #   #   #
                         # # # # # # # #
                        #               #
                       # #             # #
                      #   #           #   #
                     # # # #         # # # #
                    #       #       #       #
                   # #     # #     # #     # #
                  #   #   #   #   #   #   #   #
                 # # # # # # # # # # # # # # # #
                #                               #
               # #                             # #
              #   #                           #   #
             # # # #                         # # # #
            #       #                       #       #
           # #     # #                     # #     # #
          #   #   #   #                   #   #   #   #
         # # # # # # # #                 # # # # # # # #
        #               #               #               #
       # #             # #             # #             # #
      #   #           #   #           #   #           #   #
     # # # #         # # # #         # # # #         # # # #
    #       #       #       #       #       #       #       #
   # #     # #     # #     # #     # #     # #     # #     # #
  #   #   #   #   #   #   #   #   #   #   #   #   #   #   #   #
 # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

Wow!

Followup: I had mentioned that my code had a bug, because some pictures, such as this one:

showRule 29 11
          #
######### ###########
#         #
######### ###########
#         #
######### ###########
#         #
######### ###########
#         #
######### ###########
#         #

look different than the way Wolfram's book and web site shows them, which is like this:

          #
######### ###########
          #
######### ###########
          #
######### ###########
          #
######### ###########
          #
######### ###########
          #

After a little investigation this seems to be because Wolfram's implementation wraps around; the left neighbor of the leftmost cell of a given universe is taken from the rightmost cell, and vice-versa, while my implementation pretends that there is always more empty space available to the left and right.

Whether you consider this a bug or not is up to you. The wraparound behavior is probably considered more "canonical." You can compare the results from my program to the pictures at Wolfram's MathWorld site here. If you replace my universeToTuples function with this one:

universeToTuples :: [Bool] -> [(Bool, Bool, Bool)]
universeToTuples [] = []
universeToTuples universe = tupleize $ wrapUniverse universe
   where
       wrapUniverse xs =  ( last xs ) : ( xs ++ [ head xs ] )
       tupleize xs =
           if length xs > 3 then tuple3 xs : tupleize ( tail xs )
               else [ tuple3 xs ]
       tuple3 xs = ( xs !! 0, xs !! 1, xs !! 2 )

you will get the wraparound behavior.

Thanks for reading and as always, I appreciate your comments.

17 December 2006

From Bits to Cells: Simple Cellular Automata in Haskell, Part One

The little Haskell toy I wrote to simulate a dot-matrix printhead was useful in one sense: it helped me learn a little more Haskell! I couldn't help noticing that the results look like the transformations you get when applying simple cellular automata rules. So my next experiment is to implement those rules, and see what else I can learn from the results.

What follows is Literate Haskell, but there is no main program to run; it is for you to experiment with interactively, if you like.

Some Haskell code to implement simple cellular automata (CA) rules:
part one, round one.

The simplest type of one-dimensional CA behaves as follows:

- a cell may be on or off (one or zero).

- time ticks by in discrete intervals.

- at each tick, the state of each cell is updated using a function that
maps the current state of the cell to the left, the cell itself, and
the cell to the right to the new cell value.

Because this value-generating function depends on 3 bits of input, a
given CA rule actually consists of 2^3 = 8 possible mappings. The
mappings can be represented by the binary values 0b000 to 0b111
(decimal 0..7). By Wolfram's numbering convention for these simple
CA, we actually count down from 7 to 0. So, let's say we have a rule
that maps each possible set of left, center, and right values:

0b111, 0b110, 0b101, 0b100, 0b011, 0b010, 0b001, 0b000

to the new cell value:

0b0, 0b0, 0b0, 0b1, 0b1, 0b1, 0b1, 0b0

We can treat the resulting number, in this case 0b00011110, as the
rule number. This is rule 30 in Wolfram's schema; there are 2^8,
or 256, possible rules, numbered from 0 to 255.

First, let's make a function that, given a rule number and three
values (left, center, and right), returns the new cell state. We
turn the three values into an index to look up the appropriate bit
in the rule number.

>import Data.Bits

>genNextBit :: Int -> ( Bool, Bool, Bool ) -> Bool
>genNextBit rulenum ( left, center, right ) = rulenum `testBit` idx
>    where idx = ( if left then (4::Int) else (0::Int) ) .|.
>           ( if center then (2::Int) else (0::Int) ) .|.
>           ( if right then (1::Int) else (0::Int) )

Hmmm... it is lacking a certain elegance, and I don't quite like
the way the indentation rules work in this case, but let's test it
with a function that generates the 8 rule indices in the form of tuples:

>genRuleIndices :: [( Bool, Bool, Bool )]
>genRuleIndices = [ ( (val `testBit` 2), (val `testBit` 1),
>                   ( val `testBit` 0) ) | val <- [(7::Int),
>                   (6::Int)..(0::Int)] ]

Now if we write:

genRuleIndices

we get:

[(True,True,True),
(True,True,False),
...
(False,False,True),
(False,False,False)]

and if we write:

map (genNextBit 30) genRuleIndices

this expression curries us a function (mmm... curry!) which takes
a starting state and applies rule 30, then feeds it the possible
starting states. The result is:

[False,False,False,True,True,True,True,False]

That looks like the bit pattern for rule 30.

Just for fun, let's confirm by making a function that will
translate a list of output bit values back into a rule number.
The signature should look like this:

sumBitVals :: [Bool] -> Int

And we want the list

>test_bit_vals = [False,False,False,True,True,True,True,False]

to map back to 30. Take a shot at it yourself; you might find
that the result is instructional I'll wait.

(Musical interlude).

Did you try it? Let's look at some possible implementation strategies.
We could make a list of the powers of two and then do some list
manipulation to get the powers of two that correspond to our one-bits
summed up:

>powers_of_two_8_bits = reverse ( take 8 (iterate (2 *) 1) )

[128,64,32,16,8,4,2,1]

(An aside: since this is generated by a function that takes no
parameters at runtime, it would seem like a sufficiently smart
compiler could generate this list and stick it in the resulting
program so that no run-time calculation is required to generate
it. Does GHC do this? I don't know.)

Anyway, our implementation. We can tuple up our powers of two with
our bit values:

>tups = zip powers_of_two_8_bits test_bit_vals

[(128,False),(64,False),(32,False),(16,True),(8,True),(4,True),(2,True),(1,False)]

Then we can turn this back into a list of only the powers of two that are
"on," and sum the results:

sum (map (\ tup -> if snd tup then fst tup else 0) tups )

30

It seems like this should be simpler still. I looked in the standard
library for a function to filter one list using a list of boolean
values as a comb. Let's write one:

>combListWithList :: [a] -> [Bool] -> [a]
>combListWithList [] [] = []
>combListWithList ls comb =
>   if (head comb) then (head ls) : combListWithList (tail ls) (tail comb)
>       else combListWithList (tail ls) (tail comb)

combListWithList powers_of_two_8_bits test_bit_vals

[16,8,4,2]

That seems good, although it doesn't express the idea that the lists
have to be the same length. It still seems amazing to me that I can
reel off functions like that in Haskell and have them work right the
first time! Now we can produce our final function:

>sumBitVals :: [Bool] -> Int
>sumBitVals ls = sum $ combListWithList powers_of_two_8_bits ls

sumBitVals test_bit_vals

30

There is probably an elementary way to implement the above function
using zipWith instead of my combList, but I'll leave that to you;
leave a comment if you figure it out!

As another aside, here is another function I came up with to generate
this bit vals; it doesn't rely on a specific length.

>foldBitVals :: [Bool] -> Int
>foldBitVals ls = snd (foldr
>   (\ flag tup -> if flag then ((fst tup) * 2, (snd tup) + (fst tup))
>                           else ((fst tup) * 2, snd tup) )
>   (1, 0) ls )

foldBitVals test_bit_vals

30

Take a moment to understand the above function; it is a little strange.
I perform a fold using a tuple. The first element is the power of two,
so it keeps track of our place in the list, and the second is an
accumulator as we add up our values. This is a trick for stashing
a state that is more complex than a single value into a fold.

Anyway, that was a long diversion, so let's end for now. Next time we'll look at ways to represent our CA universe as it evolves by applying our rules. As always, your comments are welcome!

14 December 2006

The Revised Dot-Matrix Printhead

It looks like this blog has been unceremoniously removed from Planet Scheme. I guess that's only fair, since I haven't written anything specifically about Scheme in a while. Best of intentions, and all that. This blog probably belongs on "Planet Programming Language Dilettante."

Anyway, I received some helpful feedback on the short program I wrote to print bitmaps of binary numbers. Jason Creighton in particular took a shot at a rewrite, and the result as highly enlightening. I also recieved two specific pieces of advice: use short function names, and don't write a comprehension like [a | a <- [a..b]] when you can just write [a..b]. Oops. Check, and thanks!

(Semi) or (Ill) Literate Haskell follows:

Haskell Binary Printhead Toy, Second Round

>import Data.Bits

Our functions to generate lists of signed and unsigned integers of
numbits length, revised with names as suplied by Jason and error-
catching as supplied by yours truly:

>unsignedInts :: Int -> [Int]
>unsignedInts bits | (1 <= bits) && (bits <= 16) =
>   [(0::Int)..(bound::Int)]
>   where bound = 2 ^ bits - 1
>unsignedInts bits = error "expected a value from 1 to 16"

>signedInts :: Int -> [Int]
>signedInts bits | (1 <= bits) && (bits <= 16) =
>   [(- (bound::Int)) - 1..(bound::Int)]
>   where bound = 2 ^ ( bits - 1 ) - 1
>signedInts bits = error "expected a value from 1 to 16"

Jason came up with a considerably simpler way to map the bits of
the binary numbers to the characters for printing:

>boolToChar :: Bool -> Char
>boolToChar True = '#'
>boolToChar False = ' '

I'm adopting that suggestion enthusiastically, since I don't think
code can get much simpler than that!

Jason then proceeded to turn my dog food into a much more flavorful
wafter-thin after-dinner mint, leaving me struggling to understand
just what he did:

>bitPicture :: (Int -> [Int]) -> Int -> [Char]
>bitPicture intGen bits = unlines $ map stringAtIndex (reverse [0..(bits-1)])
>  where
>    stringAtIndex = map boolToChar . bitsAtIndex
>    bitsAtIndex idx = [ testBit n idx | n <- intGen bits ]

Wow. First off, Jason has supplied a type signature which expresses
what the function actually _does_. The appropriate integer generator
(the signed or unsigned version) comes in as a parameter; that's what
(Int -> [Int]) means. We next get an integer argument (the number of
bits), and finally, we spit out a list of chars.

Jason proceeds to use a couple of where clauses to set up some
simplifying expressions. Let's look at the first one:

StringAtIndex = map boolToChar . bitsAtIndex

OK, this is slightly mind-twisting right off the bat. In the world
of C and Java you can't write things like that. The dot is function
composition expressed naturally; instead of something like

stringOfBoolsFromBitsAtIndex idx = map boolToChar (bitsAtIndex idx)

we can just get Haskell to chain the functions together for us.
No fuss, no muss. But this function makes a forward reference to
bitsAtIndex, which we haven't seen yet (completely illegal in C
and Java). Let's make sure we understand what bitsAtIndex does:

bitsAtIndex idx = [ testBit n idx | n <- intGen bits ]

That's definitely more dense than the way I wrote it (and that's a
good thing). To break down the list comprehension, the expression on
the right hand side of the vertical bar says that we feed the supplied
integer generator function the number of bits to use. That gives us a
list. From that list we will draw our "n"s and test bit idx on them.
By doing it this way, Creighton has pulled an act of Haskell sleight-
of-hand: instead of generating the bitmaps for each integer value,
like I did, and then using "transpose" to rearrange them, this version
accesses the "columns" of the table directly! With those two
expressions, we've got our string representing all the bits at a
given index _for the whole supplied set of integers_.

Wow. I feel like an audience member watching a magician smash up my
watch and then give it back to me rebuilt as a marvelous whirring
gyroscope.

There's more, but this part is relatively simple:

bitPicture intGen bits = unlines $ map stringAtIndex (reverse [0..(bits-1)])

We already know what unlines does; it turns our list of strings into a
properly line-broken string. We've just described how stringAtIndex works.
and map stringAtIndex (reverse [0..(bits-1)[) just feeds stringAtIndex a
set of indices mapped from high to low, so that our most significant bit
comes out leftmost. But what does "$" do?

My reference at zvon.org says this is the "right-associating infix
application operator (f $ x = f x), useful in continuation-passing style."
I'm not sure I fully understand, and I'm not going to go into the meaning
of "continuation-passing style" right now, but I guess it saves us from
having to write:

unlines $ map (stringAtIndex (reverse [0..(bits-1)]))

to get the right order of operations, which is nice. So I'll try to
incorporate it into my own Haskell.

Anyway, let's try it. We need to drive the function:

>main :: IO ()
>main = do
>   print "Unsigned values (6 bits)"
>   putStr $ bitPicture unsignedInts 6
>   print "Signed values (6 bits)"
>   putStr $ bitPicture signedInts 6

After using GHCi's :cd and :load command to read this file, just
type:

main

"Unsigned values (6 bits)"
                                 ################################
                 ################                ################
         ########        ########        ########        ########
     ####    ####    ####    ####    ####    ####    ####    ####
   ##  ##  ##  ##  ##  ##  ##  ##  ##  ##  ##  ##  ##  ##  ##  ##
  # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

"Signed values (6 bits)"
 ################################
                 ################                ################
         ########        ########        ########        ########
     ####    ####    ####    ####    ####    ####    ####    ####
   ##  ##  ##  ##  ##  ##  ##  ##  ##  ##  ##  ##  ##  ##  ##  ##
  # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

Excellent. I laughed; I cried; I learned something today. I hope you did, too. Thanks, Jason! And thanks to the other people who commented!

Followup: for those following along the discussion in the comments, here is an alternate implementation of the bitPicture function:

bitPicture intGen bits =
    unlines $ map stringAtIndex $ [bits - 1, bits - 2..0]
    where stringAtIndex = map boolToChar . bitsAtIndex
        bitsAtIndex idx = map (flip testBit idx) (intGen bits)

Note that using [bits - 1, bits - 2..0] in the first line allows us to avoid the list reversal, so that seems like a clear win. The method of expressing bitsAtIndex above is one of several that were suggested. The idea behind the flip is that the return value of intGen should become the first parameter to bitsAtIndex, not the second. You can also express this by forcing stringAtIndex to be treated as an infix operator:

bitsAtIndex idx = map (`testBit` idx) (intGen bits)

Both of these alter the behavior of Haskell's automatic currying so that it replaces the second parameter instead of the first, leaving the first as the sole parameter to the resulting curried function; in Common Lisp, you might do this explicitly using rcurry. And don't forget the original, using a list comprehension:

bitsAtIndex idx = [ testBit n idx | n <- intGen bits ]

Assuming there is not a specific concern about efficiency, it seems like the choice between these versions might be a matter of style. Thanks again to everyone who made suggestions.

12 December 2006

The Dot-Matrix Printhead: a Haskell Toy

The body of this article is a literate Haskell program. You can extract the text, then save it as a file with the extension .lhs, then execute it using GHCi.

This is a Haskell toy to visualize twos-complement binary numbers.

In my previous articles on the subject of integer division, I alluded
to underflow and overflow and the "weird number" -- well, now you
can see them!

Almost 30 years ago I wrote code on my Radio Shack TRS-80 to make my
dot-matrix printer spit out bitmaps of binary numbers, which made
pretty recursive patterns. This program commemorates that long-ago
geekery from my early teenage hacking days.

First, we need the bits library:

>import Data.Bits

We'll also need the list library later:

>import List

Let's simulate n-bit binary number in both a signed and unsigned
interpretation. Given a number of bits, we want a function to make
a list of the possible values. For example, in unsigned form 8 bits
will give us a list of the values 0..255; in signed form, it will be
an asymmetrical list from -128..127. Note that this function generates
the whole list, so if you give it a numbits value higher than, say,
16, you might regret it.

Our calculation fails for values of numbits <= 0, so we put in a
guard case for that. We also want to discourage creating very large
lists, so we put in a guard case for numbits > 16.

>gen_n_bit_nums_unsigned numbits | numbits <= 0 = []
>gen_n_bit_nums_unsigned numbits | numbits > 16 = error "too many bits!"
>gen_n_bit_nums_unsigned numbits =
>   [ n | n <- [(0::Int)..((2::Int)^numbits - 1)] ]

>gen_n_bit_nums_signed numbits | numbits <= 0 = []
>gen_n_bit_nums_signed numbits | numbits > 16 = error "too many bits!"
>gen_n_bit_nums_signed numbits =
>   [ n | n <- [-(2::Int)^(numbits - 1)..((2::Int)^(numbits - 1) - 1)] ]

To test this we can execute:

gen_n_bit_nums_unsigned 4
[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]

gen_n_bit_nums_signed 4
[-8, -7, -6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7]

That looks about right.

Now we want a function to access the bits of those numbers from most
significant to least significant and build a list of corresponding
boolean values. This would be expressed more naturally if EnumFromTo
operated descending sequences. I'm sure there's a better way, but for
now I will just reverse the resulting list:

>rightmost_bits num_bits val | num_bits > 16 = error "too many bits"
>rightmost_bits num_bits val | num_bits <= 0 = error "too few bits"
>rightmost_bits num_bits val =
>   reverse [ testBit (val::Int) bit_idx | bit_idx <- [0..(num_bits - 1)] ]

Once we have a function which will do it with one number, doing it with a
whole list of numbers is very easy. We will make a curried function and
use map. Currying in Haskell is very easy, which makes sense given that
the operation, "to curry" and the language, Haskell, are both named after
Haskell Curry!

Our rightmost_bits function takes two arguments, but we want to make a
new function which we can hand off as an argument to map, along with a
list. This new function will only take one parameter. We can ask GHCi
what it thinks the type of the rightmost_bits function is:

:type rightmost_bits

It says:

rightmost_bits :: Int -> Int -> [Bool]

We can make a new function out of it by doing partial application:
instead of supplying all its arguments, we just supply the first one.

gen_rightmost_bits 4

If you do this in GHCi, you'll get an error message because the
result is not printable (GHCi can't apply show to the result).
But you can bind it:

let bit_generator_4 = gen_rightmost_bits 4

And then ask GHCi what its type is:

:type bit_generator_4
bit_generator_4 :: [Int] -> [[Bool]]

Take a moment to make sure you understand what we've done here:
Haskell has made us a curried function which takes a list of values
instead of a single value, and generates a list of lists. This
seems pretty close to magic! Let's make another function that
takes the number of bits and a list, generates the curried
function, then maps the list to the curried function:

>gen_rightmost_bits num_bits ls =
>   map (rightmost_bits num_bits) ls

As you can see, working with curried functions is completely
effortless in Haskell. Here's a general function to generate
a list of all the bit values for unsigned numbers of num_bits
size:

>gen_all_bits_unsigned num_bits =
>   gen_rightmost_bits num_bits (gen_n_bit_nums_unsigned num_bits)

>gen_all_bits_signed num_bits =
>   gen_rightmost_bits num_bits (gen_n_bit_nums_signed num_bits)

If we execute:

gen_all_bits_unsigned 4

We get back the following (I have reformatted the output for clarity):

[[False, False, False, False],
[False, False, False, True],
[False, False, True, False],
[False, False, True, True],
[False, True, False, False],
[False, True, False, True],
[False, True, True, False],
[False, True, True, True],
[True, False, False, False],
[True, False, False, True],
[True, False, True, False],
[True, False, True, True],
[True, True, False, False],
[True, True, False, True],
[True, True, True, False],
[True, True, True, True]]

This is recognizably a table of the unsigned binary values 0..15.
This approach doesn't seem very "Haskell-ish" yet -- instead of
generating the combinations as needed, I'm producing them all
ahead of time. It also feels to me like I'm pushing uphill on
the syntax, and having to use parentheses to force the order
of evaluation. We can think more about that later, but let's
get our pretty-printing working first.

How can I turn all those boolean values into a string I can print?
Well, we can take advantage of the fact that in Haskell, strings can
be treated as lists, and the same methods apply. This is not efficient,
but it allows us to use some of the standard functional techniques.
One of these is "fold."

Fold is a technique for boiling down a list into some kind of summary
form. Just what form this summary takes is up to you, since you can
supply your own function and a starting value. For example, to sum
the elements in a list, I could use foldl like this:

foldl (\ x y -> x + y) 0 [0, 1, 2]

The first parameter to foldl is a function. In this case I'll use
a lambda expression, which gives me a function I won't bother to
name. The second paramter is the starting value. The starting value
and the next item in the list are repeatedly passed to the function,
which accumulates a value. In this case, the final result is 3.

I can give foldl an empty string as the starting accumulator value
and my function can concatenate strings onto it:

>stringify_vals =
>   foldl (\ x y -> if y then x ++ "#" else x ++ " ") ""

Now you might want to try to apply our list of lists to this
function:

stringify_vals (gen_all_bits_unsigned num_bits)

If you do that, though, you'll get a type error like this:

"Couldn't match expected type `Bool' against inferred type `[Bool]'"

The reason can be seen if we examine the types:

:type (gen_all_bits_unsigned 4) :: [[Bool]]
(gen_all_bits_unsigned 4) :: [[Bool]]

:type stringify_vals
stringify_vals :: [Bool] -> [Char]

We have a list of lists of bools, but we want to feed our
stringification function a list of bools. To apply each
element of (gen_all_bits_unsigned 4) to our stringification
function, we can use map again:

>stringify_all_unsigned_vals num_bits =
>   map stringify_vals (gen_all_bits_unsigned num_bits)

>stringify_all_signed_vals num_bits =
>   map stringify_vals (gen_all_bits_signed num_bits)

This will give us a list of strings, one for each list of boolean
values. (I've reformatted the output).

["    ",
"   #",
"  # ",
"  ##",
" #  ",
" # #",
" ## ",
" ###",
"#   ",
"#  #",
"# # ",
"# ##",
"##  ",
"## #",
"### ",
"####"]

That's nice: we can start to see our recursive structure! But
we really want that value to be turned into a single printable
string. To achieve that, we can use the unlines function.

>prettyprint_unsigned num_bits =
>   unlines (stringify_all_unsigned_vals num_bits)

>prettyprint_signed num_bits =
>   unlines (stringify_all_signed_vals num_bits)

When we interpret

prettyprint_unsigned 3

we see a single string with newline characters escaped, like so:

"   \n  #\n # \n ##\n#  \n# #\n## \n###\n"

This will give the results we want when we pass it to putStr, like so:

putStr (prettyprint_unsigned 8)

        #
       #
       ##
      #
      # #
      ##
      ###
     #
     #  #
     # #
     # ##
     ##
     ## #
     ###
     ####

    (etc.)

Hmmm. That's nice, but it would be really cool if the results
were rotated, so that we could see what it would look like if
a printhead moving horizontally was printing out these values.
To do this, we can actually start with our intermediate list of
lists. If our four-bit numbers give us

[[False, False, False, False],
[False, False, False, True],
[False, False, True, False],
etc.

we can see the values the other way using transpose. Let's write
rotated versions of gen_all_bits_unsigned and gen_all_bits_signed:

>gen_all_bits_unsigned_rot num_bits =
>   transpose
>   (gen_rightmost_bits num_bits (gen_n_bit_nums_unsigned num_bits))

>gen_all_bits_signed_rot num_bits =
>   transpose
>   (gen_rightmost_bits num_bits (gen_n_bit_nums_signed num_bits))

And make a new stringify function to play with:

>stringify_all_unsigned_vals_rot num_bits =
>   map stringify_vals (gen_all_bits_unsigned_rot num_bits)

>prettyprint_unsigned_rot num_bits =
>   unlines (stringify_all_unsigned_vals_rot num_bits)

Now

putStr (prettyprint_unsigned_rot 6)

will give us:

                                 ################################
                 ################                ################
         ########        ########        ########        ########
     ####    ####    ####    ####    ####    ####    ####    ####
   ##  ##  ##  ##  ##  ##  ##  ##  ##  ##  ##  ##  ##  ##  ##  ##
  # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

Now that's what I'm talking about!

We can use the unsigned version to see the interesting transition
point between negative and positive values:

>stringify_all_signed_vals_rot num_bits =
>   map stringify_vals (gen_all_bits_signed_rot num_bits)

>prettyprint_signed_rot num_bits =
>   unlines (stringify_all_signed_vals_rot num_bits)

Now

putStr (prettyprint_signed_rot 6)

will give us:

 ################################
                 ################                ################
         ########        ########        ########        ########
     ####    ####    ####    ####    ####    ####    ####    ####
   ##  ##  ##  ##  ##  ##  ##  ##  ##  ##  ##  ##  ##  ##  ##  ##
  # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

Now let's give our program a main function using the IO monad and
have it print out several of our results:

>main :: IO ()
>main = do
>   print "Unsigned values (8 bits)"
>   putStr (prettyprint_unsigned 8)
>   print "Signed values (8 bits)"
>   putStr (prettyprint_signed 8)
>   print "Unsigned values (6 bits) rotated"
>   putStr (prettyprint_unsigned_rot 6)
>   print "Signed values (6 bits) rotated"
>   putStr (prettyprint_signed_rot 6)

After using GHCi's :cd and :load command to read this file, just
type:

main

I would invite any readers to play around with the code, and if you can make it more concise and/or more "Haskellish," please leave a comment.

Now if only we could get Haskell to generate the appropriate sound! My dot-matrix printer was loud. But that's a task for another day.

Update: there is a followup aticle to this one which develops the code further; see http://praisecurseandrecurse.blogspot.com/2006/12/revised-dot-matrix-printhead.html Also, make sure to take a look at the comments on both this posting and the followup.

08 December 2006

The Divisive Aftermath

Antti-Juhani pointed out to me that Haskell provides the quot and rem pair of operations in addition to div and mod. The difference is that while both satisify the identity involving division, multiplication, and the remainder, div and mod round down (towards -infinity), while quot and rem round towards zero (truncating the fractional part).

That's great! (Besides the fact that I displayed my ignorance in public, that is). It means I really can use GHCi to model my functions for testing! And I was getting worked up to write a C version of Plauger's div which rounded away from zero for negative quotients so I could have some guaranteed-consistent behavior. Instead I can use the C standard function div and Haskell's quot and rem.

So, this sort of begs the question: is there a good reason to prefer one method of integer division rounding to another? Yes. It comes down to this question: do you want your algorithms to behave smoothly across the origin, or to behave symmetrically around the origin?

Let's briefly consider some Haskell code for mapping a range of fractions onto a range of integers again.

let gen_nths denom =
    [ (numer, denom::Int) | numer <- [ (-denom::Int)..(denom::Int) ] ]

let div_map m1 m2 = map (\ tup -> fst tup *
    ( ( (m2::Int) - (m1::Int) ) `div` 2 ) `div` snd tup )

div_map (-10) 10 ( gen_nths 3 )

[-10, -7, -4, 0, 3, 6, 10]

Notice something about those values: the distances between them follow a consistent pattern. Reading the list left to right, we have distances of 3, 3, 4, 3, 3, 4. This pattern repeats. It is smooth across the origin.

Now let's change our div_map function to use quot and rem:

let quot_map m1 m2 = map ( \tup -> fst tup *
    ( ( (m2::Int) - (m1::Int) ) `quot` 2) `quot` snd tup )

quot_map (-10) 10 ( gen_nths 3 )

[-10, -6, -3, 0, 3, 6, 10]

Looking at the differences between values, from left to right, we find that they are 4, 3, 3, 3, 3, 4. The function is symmetric around the origin.

Does this matter in practice? Of course! Let's say that we have a mapping using a non-origin-crossing version of our quot_map function:

let gen_non_negative_nths denom =
    [ (numer, denom::Int) | numer <- (0::Int)..(denom::Int) ] ]
let quot_map_2 m1 m2 = map ( \tup -> fst tup *
    ( (m2::Int) - (m1::Int) ) `quot` snd tup )

quot_map_2 0 20 ( gen_non_negative_nths 6 )

[0, 3, 6, 10, 13, 16, 20]

Shiny. The intervals go 3, 3, 4, 3, 3, 4. But if we shift the values:

let gen_non_positive_nths denom =
    [ (numer, denom::Int) | numer <- [ (-denom::Int)..(0::Int) ] ]

quot_map_2 0 20 ( gen_non_positive_nths 6 )

[-20, -16, -13, -10, -6, -3, 0]

The intervals go 4, 3, 3, 4, 3, 3 -- they are backwards. In other words, the effect of integer rounding on the values of the codomain changes if you shift the domain of the denominator across the origin: the ordering of the pattern is reversed.

If we do the same thing using div:

let div_map_2 m1 m2 =
    map (\ tup -> fst tup * ( (m2::Int) - (m1::Int) ) `div` snd tup )

div_map_2 0 20 ( gen_non_negative_nths 6 )

[0, 3, 6, 10, 13, 16, 20]

div_map_2 0 20 ( gen_non_positive_nths 6 )

[-20, -17, -14, -10, -7, -4, 0]

we get a pattern of intervals (the effect of rounding) which remains the same: 3, 3, 4, 3, 3, 4.

Now, smoothness across the origin is is what I want, at least for the kinds of functions I am working on now. But your problem domain may lend itself to writing functions which involve rotations, or mapping values across the origin: in that case, you're going to want the symmetry. The important thing is to know what strategy you are using and apply it consistently. I'm really impressed that Haskell give me so many interesting options.

07 December 2006

The Division Bell Tolls for Me, Part Four (Conclusion)

Well, this was initially going to be a three-part article, but things have become more interesting. The situation is not unlike picking up a log in the forest and finding all kinds of nifty crawlies living under it. "Interesting" to certain geeky people like myself, at least! It's about to get geekier. If you are allergic to ANSI or ISO standards, you'd better stop reading now. But I claim the payoff is worth it, if you can make it through to the end. At least, for some value of "worth it."

In the last installment we were exploring the behavior of integer division in Haskell. I want to get back to C now. We learned that the behavior of integer division involving negative results is not strictly defined for C. The behavior of % (mod) is also not very strictly defined.

Let's start with Kernighan and Ritchie's The C Programming Language, Second Edition (I'm not even going to try to figure out rules for pre-ANSI/ISO C). K&R tell us on p. 10 that

...integer division truncates; any fractional part is discarded.

That isn't very detailed, but on p. 41 they tell us

the direction of truncation for / and the sign of the result for % are machine-dependent for negative operands, as is the action taken on overflow or underflow.

OK, but that doesn't describe exactly what the options are. Things get a little bit more detailed on p. 205, where we read that no guarantees are made about the sign of the remainder!

...if the second operand is zero, the result is undefined. Otherwise, it is always true that (a/b)*b + a%b is equal to a. If both operands are non-negative, then the remainder is non-negative and smaller than the divisor; if not, it is guaranteed only that the absolute value of the remainder is smaller than the absolute value of the divisor.

This isn't actually internally consistent, because if the sign of the remainder is not correct, then the identity for division doesn't work for negative values unless the sign of the quotient is allowed to be incorrect! I've never seen it suggested anywhere else that division could be that loosely defined. But K&R isn't the formal standard, so let's move on.

The reason for the "looseness" that exists in the C standard, of course, is that the standard was originally written to, as much as possible, codify existing practice, and C has a strong bias towards efficiency. Since (apparently) different hardware division algorithms had different behavior, the standards committee did not feel that it was appropriate to require existing implementations to change behavior. Doing so might have required existing architectures to have to replace hardware division with software division, which could have inflicted an enormous efficiency cost.

According to ISO/IEC 9899:1990 (known as C90):

When integers are divided and the division is inexact, if both operands are positive the result of the / operator is the largest integer less than the algebraic quotient and the result of the % operator is positive. If either operand is negative, whether the result of the / operator is the largest integer less than or equal to the algebraic quotient or the smallest integer greater than or equal to the algebraic quotient is implementation-defined, as is the sign of the result of the % operator. If the quotient a/b is representable, the expression (a/b)*b + a%b shall equal a.

So we know that there are two options for rounding. In order to maintain the identity involving integer division and mod, though, the value of mod has to be consistent with the rounding strategy, even if the sign does not have to be consistent.

There is a back door that lets us get to defined behavior. Apparently the div function does have strictly defined rounding. In section 7.10.6.2 of the standard, we read:

If the division is inexact, the resulting quotient is the integer of lesser magnitude that is the nearest to the algebraic quotient.

That's interesting -- this is rounding towards zero. The div function returns a structure of type div_t that contains both the integer quotient and integer remainder. It is provided in the <stdlib> header, available in C++ as <cstdlib>. In his book The Standard C Library, P. J. Plauger provides an implementation for div. It looks like this:

div_t (div)(int numer, int denom)
    {        /* compute int quotient and remainder */
    div_t val;
    val.quot = numer / denom;
    val.rem = numer - denom * val.quot;
    if (val.quot < 0 && 0 < val.rem)
        {    /* fix remainder with wrong sign */
        val.quot += 1;
        val.rem -= denom;
        }
return (val); }

We can see what he is doing: for negative quotients that are inexact answers to the division problem, he shifts the quotients up towards zero and the remainders get reversed. Note that he does not use the built-in % facility, presumably since its behavior does not have very strong guarantees placed on it.

C: A Reference Manual (5th Edition) by Harbison and Steele seems to indicate that the semantics are a little more rigorously defined. We read there that

For integral operands, if the mathematical quotient of the operands is not an exact integer, then the fractional part is discarded (truncation toward zero). Prior to C99, C implementations could choose to truncate toward or away from zero if either of the operands were negative. The div and ldiv library functions were always well defined for negative operands.

That seems pretty unambiguous, but then on page 419, when describing the div and ldiv function, the authors write

The returned quotient quot is the same as n/d, and the remainder rem is the same as n%d.

But that ignores the possibility of a pre-C99 implementation where / rounds away from zero for negative values, as does GHCi in my tests. So again we have a lack of rigorous internal consistency. It is worth noting here that K&R don't make any extra statements about the required rounding behavior of div and ldiv, and since K&R is still considered the "bible" by many C programmers, the guarantees on the div and ldiv function may not be very well-known -- or properly implemented.

Does C99 clarify this at all? ISO/IEC 9899:1999(E), Second Edition 1999-12-01, tells us that "the result of the / operator is the algebraic quotient with any fractional part discarded" and a footnote indicating that this is often called "truncation towards zero."

How about C++? ISO/IEC 14882, Second Edition (2003-10-15) does not actually say how rounding is performed, although it says that "If both operands are nonnegative then the remainder is nonnegative; if not, the sign of the remainder is implementation-defined." There is a footnote indicating that "According to work underway toward the revision of ISO C, the preferred algorithm for integer division follows the rules defined in the ISO Fortran standard" (rounding towards zero). When we try to look up information about div and ldiv, we find that the C++ standard just refers us to the Standard C library.

OK, all the standards language is giving me a headache; let's take a look at some code. To begin with, let's confirm the way my C implementation does rounding of integer quotients. We took a look at 6/5, -6/5, 5/6, and -5/6 in Haskell; let's look at the same quotients in C:

int val_1 = 6;
int val_2 = 5;
int val_3 = -6;
int val_4 = -5
int result_1 = val_1 / val_2;
int result_2 = val_3 / val_2;
int result_3 = val_2 / val_1;
int result_4 = val_4 / val_1;
printf("result_1 = %d, result_2 = %d, result_3 = %d, result_4 = %d\n",
    result_1, result_2, result_3, result_4);

This gives me:

result_1 = 1, result_2 = -1, result_3 = 0, result_4 = 0

Yep, I'm getting rounding towards zero for both positive and negative values. Next, some Haskell. Let's bring back our algorithm from part 2:

let gen_nths denom =
    [(numer, denom::Int) | numer <- [(-denom::Int)..(denom::Int)] ]

let div_map m1 m2 = map (\tup -> (fst tup + snd tup) *
    (((m2::Int) - (m1::Int)) `div` 2) `div` snd tup + m1)

div_map (-10) 10 (gen_nths 3)

[-10,-7,-4,0,3,6,10]

And now let's write up my C algorithm again. This time we'll use int instead of long since we don't need big numbers (although on my target, ints are actually 32 bits long).

unsigned int denom = 3;
int numer;
for ( numer = -3; numer <= 3; numer++ )
{
    int quot = numer * 10 / denom;
    printf("quotient %d\n", quot);
}

Go ahead and try it on your platform.

Whoah! What happened? You might have noticed something very odd about the first few quotients. They are monstrously large! On my target I get:

quotient 1431655755
quotient 1431655758
quotient 1431655762
quotient 0
quotient 3
quotient 6
quotient 10

What happened is that I injected a bug to point out yet another possible way your C code can fail catastrophically. This behavior can surprise even experienced programmers. Note that since my denominator is always a positive value, I declared it as an unsigned int instead of an int. So what went wrong?

Let's look at the first example: -3/3 yielding 1431655755. The behavior you are looking at is actually mandated by the C standard. When mixing signed and unsigned types, the compiler is required to promote (that is, widen) the actual type of the calculation to an unsigned type. So, internally, if we are performing the calulation -22 (signed) / 7 (unsigned), the compiler is required to notice that both a signed and unsigned long value are presented to the integer division operator. The 32-bit twos-complement representation of -3 is 0xFFFFFFFD (the conversion to an unsigned intermediate value changes the meaning, not the representation). I multiply this value by 10, yielding 0xFFFFFFE2, and then divide it by 3, yielding 0x5555554B, or a decimal value of 1431655755. The compiler sticks this unsigned value right into the signed result variable. The signed or unsigned status of the destination variable has absolutely no effect; it does not dissuade the compiler from deciding to treat the values as unsigned. (It is a common misconception to think that the destination value influences the type of the calculation in C).

If you have a target with 16-bit ints, your rolled-over values will be different (0xFFFD, 0xFFE2, and 0x554B, or decimal 21835). So the mixing of signed and unsigned values has not only given us the wrong answer, but made the code that generates the wrong answer produce inconsistent results depending on the size of int, even when the values we start with are not actually too big to fit into 16 bits.

If that wasn't perfectly clear, and even if it was, if you are really trying to write portable numeric code in C that has to be correct, I urge to consider purchasing the MISRA (Motor Industry Software Reliability Association) book MISRA-C:2004, Guidelines for the Use of the C language in Critical Systems. It has the most detailed explanation of C's arithmetic type conversions, and more importantly, the implications of these type conversions in real code. The crux of the biscuit is that C is hard to get right, even for good programmers.

OK, so let's change my denominator back to a signed long and try it again. This time the results look more reasonable: -10, -6, -3, 0, 3, and 6.

Are things any different if we use div?

int denom = 3;
int numer;
for ( numer = -3; numer <= 3; numer++ )
{
    div_t val = div(numer * 10, denom);
    printf("quotient %d\n", val.quot);
}

In my case, no, since on my compiler and target, / already seems to round towards zero (asymmetrically) while Haskell rounds symmetrically downwards (floor behavior).

By the way, while we're at it, let's see what C has to say about the "weird" number divided by -1:

long most_negative_long = -2147483648;
long minus_one = -1;
long result = most_negative_long / minus_one;
ldiv_t ldiv_result = ldiv(most_negative_long, minus_one);
printf("results: %lX, %lX\n", result, ldiv_result.quot);

What do you get? I get zero as the result of both the / and ldiv.

I also get GCC telling me "warning: this decimal constant is unsigned only in ISO C90." What does that mean? The compiler is apparently warning us that in C99 we won't be getting an implicitly unsigned value out of the constant -2147483648, in case we might have wanted to use this constant in an expression involving signed types to force the calculation to be promoted to an unsigned type. But why would we get an unsigned value in the first place? Apparently in C90, the number's magnitude is interpreted first, and the value is negated. 2147483648 is too large (by one) to fit into a signed long, so the compiler promotes the type to an unsigned long, then negates it.

I have been trying to come up with an example of an expression that behaves differently when I use -2147483648 as a constant, as opposed to 0x80000000, but so far I haven't been able to come up with one. This may be because my compiler is already operating under C99 rules and so I never get the implicit promotion to an unsigned value.

Anyway, be that as it may, some parting advice and conclusions:

1. GHCi (on the PC and Mac) and C (Visual C++ 6.0 and GCC on the PC and Mac) yield distinctly different rounding behavior for integer division. C90 allows this behavior. Does Haskell? Would it make sense for Haskell to define integer division more strictly, the way that C99 does?

2. Division by zero is never OK; neither is taking % zero. Both are a fatal error in GHCi. The results are undefined in C (but a fatal error is generated by most implementations).

3. Overflow (which, in the case of integer division, occurs only when the "weird number" is divided by -1) produces undefined results, and your code should avoid it, both in Haskell and in C.

4. Operations on unsigned values are guaranteed to produce results that are consistent from platform to platform, assuming the integer size is the same. Operations on signed values don't have this guarantee.

5. The rounding (truncation) of inexact quotients may be either towards zero in all cases, or towards zero for positive quotients and away from zero when the results are negative. If your implementation follows the newer C99 standard rounding should always be towards zero.

6. Mixing signed and unsigned values in C expressions is very dangerous. If you do so, your implementation may have errors far more severe than differences in rounding behavior.

7. Consider avoiding signed division of negative values in C altogether, if your algorithm allows it. Either implement the algorithm using entirely unsigned integral types, or shift the domain of your function so that you are operating on non-negative signed values, and shift it back (using subtraction, which has well-defined semantics) when your division is done.

8. Algorithms which feed a range of values which cross the origin to integer division may vary between GHCi and C implementations. Because GHCi does not provide guaranteed rounding towards zero as C99 and the div and ldiv functions require, it is difficult to prototype in GHCi and expect to get the same results in C.

And, finally,

9. Make sure your implementation actually works the way it is supposed to. Only testing can truly accomplish this.

Interesting, isn't it, how studying one language can enrich your understanding of another! At least, it works for me!