|
 |
I'm going to talk about Haskell. Specifically, I'm going to implement a
trivial Mandelbrot renderer in Haskell.
This entire message is an executable Haskell program. If you save the
entire text of this post, unmodified, in a file named Mandelbrot.lhs,
then you should be able to either load it into a Haskell interpreter, or
compile it into a standalone binary executable.
The technology that makes this possible is called "literate Haskell".
It's not what Knuth would regard as proper literate programming, but
it's a start. Essentially, in a normal Haskell source file (i.e., any
*.hs file), the compiler assumes that everything is code unless
explicitly marked as a comment. In a *literate* Haskell file (i.e., any
*.lhs file), everything is a comment unless explicitly marked as code.
As far as the compiler is concerned, everything you've read so far is
just code comments!
In order to mark stuff as code, you just have to begin the line of text
with ">". So let's write some Haskell code:
> module Main where
Every *runnable* Haskell program must contain a module who's name is
exactly "Main", which contains a public function who's name is exactly
"main", and who's type is exactly "IO ()". Nothing surprising there, then.
Alright, so if we're going to do a Mandelbrot renderer, probably the
first thing we'll want is complex number arithmetic. The easiest way to
do that is this:
import Data.Complex
But where's the fun in that? Nah, let's reimplement it ourselves!
> data Complex = Cx Double Double
This line of code does several things.
- "data" is a Haskell reversed word. In general, Haskell tends to use
weird symbols rather than keywords; data is one of a small number of
exceptions. The data keyword is for defining new data types.
- "Complex" is the name of the new type we're defining. Everything after
the equals sign is the type definition.
- "Cx" is a 'value constructor function'. It's basically a new
identifier used for creating and inspecting Complex values. (It'll make
more sense when I show you some examples.)
- "Double Double" indicates that a Cx record has two (unnamed) fields,
both of type Double. Haskell's Double type does pretty much what the C
double type does. (With the exception that Haskell actually *guarantees*
IEEE-754 double-precision semantics, while C does not.)
[Note: When a type has only one value constructor function, it is
customary for it to have the same name as the type. In other words,
data Complex = Complex Double Double
I have deliberately not done this, in the interests of clarity.]
At this point, we can write an expression such as
Cx 5 7
which constructs a new Complex value, containing the numbers 5 and 7. In
other words, this represents 5 + 7i. So that's how you put data *into* a
Complex value; now how do you get it *out* again?
To inspect the contents of a data structure, Haskell uses 'pattern
matching'. For example:
> realPart :: Complex -> Double
> realPart (Cx r i) = r
>
> imagPart :: Complex -> Double
> imagPart (Cx r i) = i
The first line says that "realPart" is a function which takes a Complex
as input, and returns a Double as output. The second line says how the
function actually works; the stuff in brackets is a 'pattern'. It
defines a new variable named "r" and sets it to whatever is in the first
field of Cx, and a new variable called "i" and sets it to whatever is in
the second field of Cx. The "realPart" function returns "r", but the
"imagPart" function returns "i".
At this point, we could say
realPart (Cx 5 7)
and the answer would be 5.
Note that Haskell uses a syntax that's rather similar to the Unix "sh"
shell: items are space-delimited, and the *first* item is a function to
execute. All remaining items are *arguments* to that function. (Brackets
begin a nested subexpression, where again the first item after the
open-bracket is a function to execute.) Also like "sh", Haskell has
normal infix operators such as "+", "-", etc.
Patterns can of course be used for more than just fetching data out of
structures though. For example,
> isReal :: Complex -> Bool
> isReal (Cx r 0) = True
> isReal _ = False
Here there are two 'equations' for "isReal", each with a different
pattern. Patterns are tried in the order they're written in the source
code (so ordering matters!). In this instance, the first pattern matches
if and only if the imaginary part is exactly 0. Otherwise, the second
pattern is "_", which is a "match anything" wildcard.
One thing that often confuses newcomers is that Haskell provides
multiple seemingly equivalent ways to do conditional branching. We could
easily rewrite the above using an if-expression instead:
isPart (Cx r i) = if i == 0 then True else False
This is of course equivalent to simply writing
isReal (Cx r i) = i == 0
Notice the difference between "=" and "==". Getting this wrong is a
compile-time error, not a run-time bug.
In C-based languages, if-statements are a language primitive, whereas
switch/case blocks are just syntax sugar for making nested if-blocks
easier to write. In Haskell, the reverse is true; *pattern matching* is
the fundamental language primitive behind all conditional branching. An
if-expression is simply syntax sugar for a pattern match where the
patterns are "True" and "False". Likewise, the "==" function is
implemented in terms of pattern matching, NOT THE OTHER WAY AROUND! (You
can pattern-match on stuff for which the "==" function isn't defined!)
So much for putting stuff in and getting stuff out; how about some
actual *arithmetic* now?
> addCx :: Complex -> Complex -> Complex
> addCx (Cx r1 i1) (Cx r2 i2) = Cx (r1 + r2) (i1 + i2)
That type signature looks a bit weird, eh? The basic rule for reading
function signatures is this: The *last* type is the return type of the
function. All the *other* types are argument types. So if you see 5
things separated by arrows, than you have a 4-argument function, and the
5th type is what the function returns.
The actual function body should be fairly self-explanatory. Let's define
some more arithmetic:
> subCx :: Complex -> Complex -> Complex
> subCx (Cx r1 i1) (Cx r2 i2) = Cx (r1 - r2) (i1 - i2)
>
> mulCx :: Complex -> Complex -> Complex
> mulCx (Cx r1 i1) (Cx r2 i2) = Cx (r1*r2 - i1*i2) (r1*i2 + r2*i1)
No surprises here. Notice that Haskell applies normal operator
precedence rules to arithmetic. (Unlike some - I'm looking at you,
Smalltalk!)
At this point, we can say
quadratic c z = addCx (mulCx z z) c
This is, of course, awful. This kind of ugliness is sadly necessary in
defective languages such as Java. Fortunately, Haskell allows us to
create user-defined operators, and even to set arbitrary precedence
rules for this:
> (<+>) = addCx
> (<->) = subCx
> (<*>) = mulCx
>
> infixl 6 <+>
> infixl 6 <->
> infixl 7 <*>
Now we can say
quadratic c z = z <*> z <+> c
This is much better, but still not perfect. Why the weird "<+>"? Why not
just name our function "+"?
Some programming language allow you to "overload" a function. That is,
you can define several functions with the name type, so long as their
type signatures are different. When you call the function, the compiler
uses the argument types [which have to be explicitly spelled out every
two seconds] to determine which function you intended to call.
Haskell works the other way around: When you call a function, the
compiler uses the function's type signature to automatically determine
the types of the function's arguments. This obviously fails completely
unless there is only one function with that name.
We *can* define a function named "+". However, this doesn't overload the
built-in "+" function, it creates a new function with the same name as
an existing one. That would mean that every time we want to write "+",
we have to write a fully-qualified name:
foo :: Double -> Double -> Double
foo x y = x Prelude.+ y
bar :: Complex -> Complex -> Complex
bar x y = x Main.+ y
That's just unspeakably ugly!
Fortunately, Haskell provides a way around this using so-called
"type-classes".
I want you to ignore the word "class" right now. A type-class is
*nothing like* what an OO programmer would recognise as being a class.
Rather, Haskell's type-classes correspond almost exactly to what Java
and C# call "interfaces".
For those that don't know, an "interface" is a set of method
definitions. When you define a new class, it can "implement" zero or
more interfaces - which basically means you have to provide
implementations for the stipulated methods.
Haskell type-classes work almost exactly like this, but with a
difference: The interfaces that a C# class implements are determined
when you write the class. But with Haskell, a type can implement new
interfaces at any time. I can take a type defined in one package, and an
interface defined in some completely unrelated package, and make one
implement the other in my own package.
Enough waffle. What incantation must be utter to fix our complex
arithmetic? Well, the standard libraries contain the following definition
class Num n where
(+) :: n -> n -> n
(-) :: n -> n -> n
(*) :: n -> n -> n
That means we just have to write the following:
> instance Num Complex where
> (+) = addCx
> (-) = subCx
> (*) = mulCx
Notice how the subsequent lines are indented; white-space is significant
in Haskell. (!)
(Remember, class == interface, instance == implementation.)
Now - finally - we can write
> quadratic :: Complex -> Complex -> Complex
> quadratic c z = z*z + c
Not only that, but every single Haskell function ever written that works
with number types will now work with Complex. For example,
$ sum [Cx 5 7, Cx 1 1, Cx 2 3]
Cx 8 11
If you do compile this post, you may get several warnings about
"undefined methods". The Num class actually defines more than just the
three methods above. We've written an implementation for Num without
implementing *all* of the prescribed methods. Haskell is usually very
strict about correctness and compile-time checks. But not implementing
all the methods of a class is a *run-time* error. (Or rather, trying to
*use* one of the unimplemented methods is a run-time error. Since we
aren't going to use any of the unimplemented methods, we'll be fine.)
In case you were worried, I should point out that we could have just written
instance Num Complex where
(Cx r1 i1) + (Cx r2 i2) = Cx (r1 + r2) (i1 + i2)
...
in the first place, without having to define "addCx" and so on first. In
the real world, that is how you would do it; I took the long way around
for example's sake.
OK, so where are we up to now? So far, we have a "quadratic" function
that implements the standard Mandelbrot iteration function. Now we write
a while-loop that counts how many iterations it takes for the bail-out
condition to be met...
...erm, no. This is *functional* programming. We do things a bit
differently here.
First of all, there is a library function named "iterate". It takes a
function and a start value, and generates an *infinite list* of outputs:
iterate f x = [x, f x, f (f x), f (f (f x)), ...]
(The square brackets indicate a Haskell list, with the items being
comma-separated.)
Generating an infinite list would of course take an infinite amount of
time. However, provided you only "look at" a finite portion of this
list, Haskell's lazy evaluation ensures that only a finite amount of
work is performed.
Knowing this, let us define the following function:
> orbit :: Complex -> [Complex]
> orbit c = iterate (quadratic c) (Cx 0 0)
This function takes a single complex number (C), and yields an infinite
list of complex numbers (Z). The type "[Complex]" indicates a list of
Complex values. Recall that "quadratic" takes two arguments, yet I only
gave it one. We are therefore left with a function waiting for one
remaining argument - z, as it happens. If I had defined "quadratic" with
its arguments the other way around, I wouldn't be able to do this little
short-cut.
While we're on the subject, Haskell supports anonymous functions. We
didn't need to actually define a named function at all; we could have
just said
orbit c = iterate (\ z -> z*z + c) (Cx 0 0)
The backslash looks like some kind of weird escape character, but it's
actually supposed to be an ASCII stand-in for the Greek letter lambda
(λ) - which some of you may remember from HalfLife. The stuff between
the backslash and the arrow is the function's arguments; stuff after the
arrow is what the function returns. Notice, especially, how z is a
function argument, but c is not, and yet the function can access c!
Usually the Mandelbrot iteration has two stopping conditions: either the
magnitude of Z exceeds some constant (typically 2), or the maximum
number of iterations is exceeded. (More fancy renderers check for
periodic cycles in the output, but that involves tricky limited
precision matching problems...)
We can implement a hard limit on the maximum number of iterations using
the "take" function, which trims a list to a maximum length:
$ take 3 [1, 2, 3, 4, 5, ...]
[1, 2, 3]
To implement the other stopping condition, first we need to be able to
calculate the magnitude of Z:
> mag2 :: Complex -> Double
> mag2 (Cx r i) = r*r + i*i
Then we can use the "takeWhile" function:
takeWhile (\ z -> mag2 x <= 4) mylist
Here I've used the old trick of computing |Z|^2 rather than |Z| to avoid
an expensive square-root operation. Rather than compare |Z| against 2,
we compare |Z|^2 against 2^2 (=4).
Putting it all together, then, we can do
> count :: Complex -> Int
> count c = length (take 15 (takeWhile (\ z -> mag2 z <= 4) (orbit c)))
The "length" function takes a list, and returns an Int representing how
many items are in the list. Haskell's "Int" type is guaranteed to be
twos-complement, with a minimum precision of 31 bits. Yes, you read
correctly: 31, not 32 bits. Haskell's standard compiler, GHC, gives you
32 bits on IA32, and 64 bits on AMD64. If you want a specific precision,
you should say so: Int8, Int16, Int32 and Int64 exist for this purpose.
You could even say Integer, which uses the GMP (the GNU Multi-Precision
library) to give you arbitrary-precision integers. But hey, we're only
trying to count to 15, so Int will be fine.
My, that's a lot of brackets! Fortunately, Haskell has something a bit
like shell piping:
count c = length $ take 15 $ takeWhile (\ z -> mag2 z <= 4) $ orbit c
No more counting brackets! Add or remove any stage easily! But notice,
unlike a shell pipeline, this still reads from right to left, like it
did when it was all brackets. (This makes changing from brackets to
dollars much easier.)
Alternatively, you can do this "point-free style":
count = length . take 15 . takeWhile (\ z -> mag2 z <= 4) . orbit
Before, we were passing the result of one function as the argument to
another function. Now "c" has completely vanished (this is the "point"
that the expression is now "free" of). Instead, we're doing mathematical
composition of functions, like in the textbooks. Notice that this
*still* reads from right to left - and indeed, this is customary in
mathematics. (Although a mathematician would typically write g ∘ f
rather than g . f for this.)
Whatever way you code it, our "count" function takes a complex number C,
and computes the length of the orbit of 0 that escapes to infinity. Now
we just need to do this for every pixel, and render the result somehow.
There are several ways we *could* do this. We could generate a CSV file
of the results, or even a PPM image. (You'll need a decent image viewer
to display it though; no standard Windows program understands PPM.) But
I'm going to do something more trivial: ASCII art!
With a maximum iteration limit of 15, we need to find 16 different
characters, hopefully representing different levels of blackness. One
inefficient but simple way to do so is this:
> int_to_char :: Int -> Char
> int_to_char n = " .:;^_/|~-+*%&@#" !! n
The "!!" operator returns the Nth element of a list - and in Haskell, a
string is a list of characters. When I say "list", I mean *single-linked
list*; hence, the "!!" operator is O(n). But for such tiny numbers, it
doesn't matter too much. (Arrays, with their O(1) access, provide a "!"
operator instead.)
Now to generate a grid of coordinates. The easiest way is with a "list
comprehension": Haskell allows you to write [5 .. 15] to represent a
list containing all the numbers in the specified range. We don't want
integers though, so we need to specify a step value:
> gridX :: [Double]
> gridX = [-2, -1.95 .. 2]
>
> gridY :: [Double]
> gridY = [-2, -1.875 .. 2]
Now to combine them. When I first started using Haskell, I was
astonished that there is no function in the standard libraries for
performing a Cartesian product of lists. But it turns out that it's
really, *really* trivial to do this using the list monad, or list
comprehensions. I will show the latter:
> grid1 :: [Complex]
> grid1 = [ Cx r i | r <- gridX, i <- gridY ]
That's great, but we probably actually want to know where one row ends
and the next begins. So let's try again:
> grid2 :: [[Complex]]
> grid2 = [ [ Cx r i | r <- gridX ] | i <- gridY ]
That's a list of lists. Now we want to take our grid of complex numbers
and turn it into a grid of... well, characters, actually.
> output :: [[Char]]
> output = map (\ row -> map (int_to_char . count) row) grid2
There's a lot happening there. Let's break it down:
- The "map" function applies a function to every row of grid2.
- The anonymous function in their takes each row of data, and uses map
again to apply count and then int_to_char to the complex number it finds
in each "pixel". (Remember, function composition operates right to left!)
This variable, then, contains our cheesy Mandelbrot ASCII-art. At this
point we're done... unless you want a *runnable* program, in which case
you need to *output* this gold somehow. To do that, you need Haskell's
I/O facilities, which are obviously extremely hard to understand.
> main = mapM_ putStrLn output
The "mapM" function is similar to "map", except that it works for
functions that perform a monadic action (as "putStrLn" does). The "mapM"
function collects the results into a new list; however, printing stuff
doesn't produce a result, so "mapM_" (with an underscore) is used to
avoid constructing a pointless list of nothings.
We could also have done it this way:
main = putStrLn (unlines output)
The "lines" function takes a string and splits it into a list of
strings, one string for each line in the original input. The "unlines"
function does the opposite - it turns a list of strings into one giant
string with newlines in it. The "putStrLn" function then prints this out.
If you now run our program, it will dump out some nasty ASCII-art graphics.
Q.E.D.
Post a reply to this message
|
|
|
|
> [Note: When a type has only one value constructor function, it is
> customary for it to have the same name as the type. In other words,
>
> data Complex = Complex Double Double
>
> I have deliberately not done this, in the interests of clarity.]
It's perhaps not clear what I'm talking about here. Consider, for example
data Complex = Grid Double Double | Polar Double Double
Now I can write a complex number in *two* ways:
Grid 5 7
which represents 5 + 7i, or
Polar 3 pi
which represents -3 (i.e., magnitude=3, argument=pi radians). I can then
write
mulCx (Grid r1 i1) (Grid r2 i2) = Grid (r1*r2 - i1*i2) (r1*i2 + r2*i1)
mulCx (Polar m1 a1) (Polar m2 a2) = Polar (m1*m2) (a1+a2)
Unfortunately, if one number is Grid and the other is Polar, this
function will throw a "pattern match failure" exception, because no
pattern matches that exact combination. I'm not saying this approach is
a good idea, I'm offering it as an example of what you *could* do.
Now, "Grid" and "Polar" are value constructors; they are ways of
constructing (and de-constructing) data of the Complex type. But our
original definition had only one constructor. Why have two names for one
thing? This is why people typically write
data Complex = Complex Double Double
I didn't, though, to make it clear when I'm talking about a type and
when I'm talking about a constructor.
> In C-based languages, if-statements are a language primitive, whereas
> switch/case blocks are just syntax sugar for making nested if-blocks
> easier to write. In Haskell, the reverse is true; *pattern matching* is
> the fundamental language primitive behind all conditional branching. An
> if-expression is simply syntax sugar for a pattern match where the
> patterns are "True" and "False". Likewise, the "==" function is
> implemented in terms of pattern matching, NOT THE OTHER WAY AROUND! (You
> can pattern-match on stuff for which the "==" function isn't defined!)
The general way to do pattern-matching is with a case-block:
case z of
Grid r i -> ...stuff...
Polar m a -> ...stuff...
An if-expression is merely syntax sugar:
if foo then bar else baz
case foo of
True -> bar
False -> baz
Writing multiple equations for a single function is also syntax sugar:
isReal (Cx r 0) = True
isReal _ = False
isReal z =
case z of
Cx r 0 -> True
_ -> False
> So much for putting stuff in and getting stuff out; how about some
> actual *arithmetic* now?
>
> > addCx :: Complex -> Complex -> Complex
> > addCx (Cx r1 i1) (Cx r2 i2) = Cx (r1 + r2) (i1 + i2)
>
> That type signature looks a bit weird, eh? The basic rule for reading
> function signatures is this: The *last* type is the return type of the
> function. All the *other* types are argument types. So if you see 5
> things separated by arrows, than you have a 4-argument function, and the
> 5th type is what the function returns.
If you want to make your head explode:
addCx :: Complex -> (Complex -> Complex)
Similarly,
quadratic :: Complex -> (Complex -> Complex)
This is why "quadratic c" is a valid expression; the result is a
function of type Complex -> Complex.
> If you do compile this post, you may get several warnings about
> "undefined methods". The Num class actually defines more than just the
> three methods above. We've written an implementation for Num without
> implementing *all* of the prescribed methods. Haskell is usually very
> strict about correctness and compile-time checks. But not implementing
> all the methods of a class is a *run-time* error. (Or rather, trying to
> *use* one of the unimplemented methods is a run-time error. Since we
> aren't going to use any of the unimplemented methods, we'll be fine.)
Haskell's number classes are notoriously broken. They were defined very
early on in Haskell's development. (Numbers are one of the very first
constructs you're going to need out of a programming language!)
Naturally, to change them now would massively break every shred of
Haskell code ever written. On top of that, nobody can decide on a *less*
broken alternative.
For the curious, the *actual* definition of Num is
class Num n where
(+) :: n -> n -> n
(-) :: n -> n -> n
(*) :: n -> n -> n
negate :: n -> n
abs :: n -> n
signum :: n -> n
fromInteger :: n -> n
You may notice that abs is a rather stupid function for a complex number
to have; what do you do? Take the absolute value of the real and
imaginary parts separately? (Note that what you *can't* do is have abs
return a real-value as the answer; it must be another complex value.)
Ditto for signum, which usually returns -1, 0 or +1 for normal types.
Finally, fromInteger converts an integer into a Complex (or whatever).
Note that I said earlier that "sum" now works for Complex. I was wrong;
without fromInteger, you can't construct the complex number 0, which is
what "sum" starts adding up from. Similarly, "product" wants to do
"fromInteger 1" for its start value.
Things rapidly get vastly more broken if you want to handle *fractional*
values. Notice (*) is here, but (/) is conspicuously missing? That's
because integer types have "div", which does integer division, while
fractional types have "/", which does normal division. (Integer types
also have "mod", while fractional types do not.) Trigonometric functions
are even more fun!
> Knowing this, let us define the following function:
>
> > orbit :: Complex -> [Complex]
> > orbit c = iterate (quadratic c) (Cx 0 0)
Now you see why it's Cx 0 0, rather than just 0. If you add
> fromInteger r = Cx (fromInteger r) 0
to the instance declaration, then these problems go away. (And "sum" and
"product" will start working properly.) Notice that "r" is an Integer,
so we need to recursively call "fromInteger" to convert it to a Double.
> My, that's a lot of brackets! Fortunately, Haskell has something a bit
> like shell piping:
>
> count c = length $ take 15 $ takeWhile (\ z -> mag2 z <= 4) $ orbit c
This is not actually part of the language specification. Actually, it's
just a user-defined operator with *really* low precedence:
> ($) :: (x -> y) -> x -> y
> f $ x = f x
>
> infixr 0 $
That *looks* like a really pointless function, but it's utility is in
the operator precedence.
> We could also have done it this way:
>
> main = putStrLn (unlines output)
You can also do this:
> main = writeFile "Mandelbrot.txt" (unlines output)
This then saves the data to a text file (in the system default encoding,
with the system default end-of-line convention), rather than printing it
to stdout.
Post a reply to this message
|
|
null> wrote: