Day 17: Making Snowmen

by

I started out planning this day to be about complex numbers in Perl 6. But after I thought about it a bit, I decided that the complex number implementation is so straightforward explaining it would make a pretty boring gift. Instead, let’s explore the Mandelbrot set, which will let us do a bit of complex math, look at pretty pictures, and hint at some advanced features of Perl 6, too.

Without further ado, here’s the first version of the script:

use v6;

my $height = @*ARGS[0] // 31;
my $width = $height;
my $max_iterations = 50;

my $upper-right = -2 + (5/4)i;
my $lower-left = 1/2 - (5/4)i;

sub mandel(Complex $c) {
    my $z = 0i;
    for ^$max_iterations {
        $z = $z * $z + $c;
        return 1 if ($z.abs > 2);
    }
    return 0;
}

sub subdivide($low, $high, $count) {
    (^$count).map({ $low + ($_ / ($count - 1)) * ($high - $low) });
}

say "P1";
say "$width $height";

for subdivide($upper-right.re, $lower-left.re, $height) -> $re {
    my @line = subdivide($re + ($upper-right.im)i, $re + 0i, ($width + 1) / 2).map({ mandel($_) });
    my $middle = @line.pop;
    (@line, $middle, @line.reverse).join(' ').say;
}

So, lines 3-5 set up the pixel size of the graphic we will create. @*ARGS is the new name of the command line argument array. The // operator is the new "defined" operator; it returns its first argument if that argument is defined, its second otherwise. In other words, line 3 sets $height to be the first argument on the command line, or 31 if no such argument was set. $width is set equal to $height -- the code is set up to generate a square graphic right now, but the variables are set apart for ease of future hacking. $max_iterations sets how many times the core Mandelbrot loop will iterate before it concludes a point is in the set. (Because we're relying on the symmetry of the Mandelbrot set, $width must be odd.)

Lines 7-8 set the boundaries of our image on the complex plane. Introducing the imaginary component of a number is as simple as giving a number (or numeric expression) followed by i; this creates a number of the Complex type. Complex math works pretty much the way you would expect it to, for example, (as we see here) adding a Complex to an Int or a Rat yields another Complex.

Lines 10-17, then, are the core Mandelbrot function. To quickly explain, a complex number c is in the Mandrelbrot set if the equation z = z * z + c (with initial z of 0) stays bounded as we iterate the equation. This function implements exactly that in Perl 6. We set up a loop to iterate $max_iterations times. It is known that once |z| grows bigger than 2 it will not stay bounded, so we use $z.abs > 2 to check for that condition. If it is true, we leave the loop early, returning 1 from the function to indicate the corresponding pixel should be black. If the loop finishes the number of iterations without exceeding those bounds, we return 0 for the color white.

Lines 19-21 are a simple helper function to return a list of a simple arithmetic progression from $low to $high with $count elements. Note that $low and $high have no specified type, so any type (or even pair of types) that the basic arithmetic operators will work on will work here. (In this script, we use it first for Num, and then for Complex.)

Lines 23-24 print the header for the header for a PBM file.

Lines 26-30 print the actual image data. $upper-right.re is the real part of the complex number $upper-right, and $upper-right.im is the imaginary part. The loop iterates over the real part of the range. Inside the loop, we subdivide again along the imaginary part to generate a list of the complex values we are interested in examining for one half of this row of the image. We then run that list through the mandel function using map, generating a list of 0s and 1s for half of the row, including the midpoint.

We do it this way because the Mandelbrot set is symmetric about the imaginary axis. So we then pop that midpoint, and make a new list which is the old list (minus the midpoint), the midpoint, and the list (minus the midpoint) reversed. We then feed that to join to make a string for the entire line, and finally say to print it out.

Note that doing it this way rotates the Mandelbrot set 90 degrees from the way it is normally displayed, giving us a lovely snowman shape:
White on black Mandelbrot set

With the current Rakudo, this is quite slow, and prone to crash randomly depending on the size of the image you are generating. However, imagine for a minute that happy future day when Rakudo is not only snappy, but handles automatic hyperoperator threading as well. At that point, it will be easy to make a parallel version of this script by changing the map call to a hyperoperator.

There's only one tricky bit: there's no way to have a hyperoperator call a normal sub. They only call class methods and operators. So, as a first approximation which works in current Rakudo, we can "augment" the Complex class to have a .mandel.

use MONKEY_TYPING;

augment class Complex {
    method mandel() {
        my $z = 0i;
        for ^$max_iterations {
            $z = $z * $z + self;
            return 1 if ($z.abs > 2);
        }
        return 0;
    }
}

for subdivide($upper-right.re, $lower-left.re, $height) -> $re {
    my @line = subdivide($re + ($upper-right.im)i, $re + 0i, ($width + 1) / 2)>>.mandel;
    my $middle = @line.pop;
    (@line, $middle, @line.reverse).join(' ').say;
}

The only difference to mandel is it is now a method, and the role of the former $c argument is taken by self. Then instead of map({mandel($_)}) we use the hyperoperator.

As I said, this version works now. But personally, I'm a little uncomfortable augmenting an existing class like that; it feels dirty to my old C++ instincts. We can avoid it by turning mandel into an operator:

sub postfix:<☃>(Complex $c) {
    my $z = 0i;
    for ^$max_iterations {
        $z = $z * $z + $c;
        return 1 if ($z.abs > 2);
    }
    return 0;
}

for subdivide($upper-right.re, $lower-left.re, $height) -> $re {
    my @line = subdivide($re + ($upper-right.im)i, $re + 0i, ($width + 1) / 2)>>☃;
    my $middle = @line.pop;
    (@line, $middle, @line.reverse).join(' ').say;
}

This takes advantage of Perl 6's Unicode ability to have a little fun, defining the operator using the snowman symbol. This ☃ operator works fine in Rakudo today, but alas the >>☃ hyperoperator does not work yet.

Thanks to Moritz and TimToady for suggests and help with this code. Two versions of the script (one full color) are up at github, if you'd like to play around with them.

Update (4/18/2010): I've ported the color version of the script at github to the latest version of Rakudo. It's quite slow, and uses phenomenal amounts of memory, but unlike the previous version it is rock-solid stable. Here's a 1001x1001 full color Mandelbrot set that took it 14 hours and 6.4 GB of memory to compute.

About these ads

4 Responses to “Day 17: Making Snowmen”

  1. illviljan Says:

    The idiom of using

    (^$count)

    made my day :-). Nice and compact way of generating numbers 0 up to $count (minus 1)!

  2. Dave Whipp Says:

    Instead of augmenting the class, presumably we could also have composed the result of subdivide with a Role that does the mandel:

    Role Mandel { method mandel { … } };

    sub subdivide($low, $high, $count) {
    (^$count).map: {
    ($low + ($^value / ($count – 1)) * ($high – $low)) with Mandel
    }
    }

    Or perhaps:

    my @line = subdivide($re + ($upper-right.im)i, $re + 0i, ($width + 1) / 2).map({ $_ with Mandel })>>.mandel;

    Still feels a but ugly

  3. Mandelbrot « Just Rakudo It Says:

    [...] By colomon I’ve dragged the Mandelbrot script from day 17 of the Advent calendar kicking and screaming to the current version of Rakudo. The big hurdle turned out to be an [...]

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s


Follow

Get every new post delivered to your Inbox.

Join 37 other followers

%d bloggers like this: