Just Don’t Rakudo It?!

May 31, 2011

After the announcement that the Niecza v6 now supports full bigint Ints, I finally have downloaded it and gotten it working on my MacBook Pro. (The tricky bit was getting mono, making niecza work was trivial after that.)

And with a couple of trivial modifications, I’ve got the old mandelbrot-color script ported to niecza. There is one major issue: Complex.abs doesn’t work properly. The resulting set is still correct, I think, but the coloring is off.

Here’s the crazy bit: it can calculate a 1001×1001 Mandelbrot set in under 5 minutes on my MBP. By comparison, a 145×145 Mandelbrot set in Rakudo runs in about 5 and a half minutes. That’s niecza clocking in at roughly 47 times faster than Rakudo!

My initial impression is that Niecza still has a lot of rough edges and unimplemented features. But it seems to me that with this release, it’s gone from a cool up-and-coming Perl 6 to a real contender.

I guess the real question here is how Rakudo Nom will handle it. Getting rid of all the endless object creation stuff should be a big win for Complex-heavy math code.

Update: And additional quick benchmarks on string code show no advantage what-so-ever for Niecza — even factoring out Niecza’s crazily slow startup time? I suppose I might have managed to choose the worst-possible comparison (from Rakudo’s point of view) for my first benchmark…

Perl 6 resources

April 30, 2011

I’m writing this post in as a reference for the Perl 6 lecture I’m giving today at Penguicon, but of course it may prove generally useful as well.

Official Perl 6 website
#perl6 IRC
STD.pm6 (official Perl 6 grammar)
Perl 6 Spec
Rakudo website
Rakudo repository on github
masak’s History of Perl 6 (Only up to mid-2010, but still very interesting.)
ABC module (used as an example, interesting real-world grammar usage)
jnthn’s Perl 6 talks

An infinite stream of “Pi”

April 27, 2011

So, after TimToady’s help with my last problem, finishing this is trivial. You just convert the Haskell code without worrying about type safety.

type LFT = (Integer, Integer, Integer, Integer) 
extr :: LFT -> Integer -> Rational 
extr (q,r,s,t) x = ((fromInteger q) * x + (fromInteger r)) / 
                           ((fromInteger s) * x + (fromInteger t)) 
unit :: LFT 
unit = (1,0,0,1) 
comp :: LFT -> LFT -> LFT 
comp (q,r,s,t) (u,v,w,x) = (q*u+r*w,q*v+r*x,s*u+t*w,s*v+t*x) 


sub extr([$q, $r, $s, $t], $x) { 
    ($q * $x + $r) / ($s * $x + $t); 

my $unit = [1, 0, 0, 1];

sub comp([$q,$r,$s,$t], [$u,$v,$w,$x]) {
    [$q * $u + $r * $w, 
     $q * $v + $r * $x, 
     $s * $u + $t * $w, 
     $s * $v + $t * $x];

And then the final piece in the puzzle,

pi = stream next safe prod cons init lfts where 
  init = unit 
  lfts = [(k, 4*k+2, 0, 2*k+1) | k<-[1..]] 
  next z = floor (extr z 3) 
  safe z n = (n == floor (extr z 4)) 
  prod z n = comp (10, -10*n, 0, 1) z 
  cons z z’ = comp z z’ 


sub pi-stream() {
    stream(-> $z { extr($z, 3).floor; },
           -> $z, $n { $n == extr($z, 4).floor; },
           -> $z, $n { comp([10, -10*$n, 0, 1], $z); },
           (1..*).map({ [$_, 4 * $_ + 2, 0, 2 * $_ + 1] }));

It’s a very direct translation.

Does it work?

> my @pi := pi-stream;
> say @pi[^40].join('');


Except, according to the Joy of Pi, the first 40 digits of pi are

3.1415926535 8979323846 2643383279 502884197 # pi
3.1415926535 8979323846 8163213056 056860170 # ours

What’s going wrong? I haven’t empirically verified it yet, but I’m pretty sure the issue is Rakudo’s Ints and Rats overflowing. Which means our next post is going to have to dive back into Math::BigInt and Math::FatRat…

More Pi

April 26, 2011

So, in my previous post I started converting a spigot algorithm for calculating pi from Haskell to Perl 6. I apologize for being away for so long, but I’m back at it now.

Interestingly, while I thought the previous stream function was much clearer in p6, this time out I think I have to give the edge to Haskell.

convert :: (Integer,Integer) -> [Integer] -> [Integer] 
convert (m,n) xs = stream next safe prod cons init xs 
    init = (0%1, 1%1) 
    next (u,v) = floor (u*v*n’) 
    safe (u,v) y = (y == floor ((u+1)*v*n’)) 
    prod (u,v) y = (u - fromInteger y/(v*n’), v*n’) 
    cons (u,v) x = (fromInteger x + u*m’, v/m’) 
    (m’,n’) = (fromInteger m, fromInteger n) 

The difference comes from Haskell’s extremely elegant on-the-fly pair notation. When I translate that to p6, I get

sub convert($m, $n, @x) {
    stream(-> $u { floor($u.key * $u.value * $n); },
           -> $u, $y { $y == floor(($u.key + 1) * $u.value * $n); },
           -> $u, $y { $u.key - $y / ($u.value * $n) => $u.value * $n; },
           -> $u, $x { $x + $u.key * $m => $u.value / $m; },
           0/1 => 1/1,

Even with p6’s big advantage in not having to explicitly convert integers to rationals, the pair thing makes this round a win for Haskell, IMO.

Perhaps one of the other p6 programmers out there can think of a more elegant way of handling this…

Update: They sure can! The esteemed TimToady pointed out Perl 6 can do something almost identical to the Haskell approach, skipping Pairs altogether:

sub convert($m, $n, @x) {
    stream(-> [$u, $v] { floor($u * $v * $n); },
           -> [$u, $v], $y { $y == floor(($u + 1) * $v * $n); },
           -> [$u, $v], $y { [$u - $y / ($v * $n), $v * $n]; },
           -> [$u, $v], $x { [$x + $u * $m, $v / $m]; },
           [0/1, 1/1],

This version compares very well with the Haskell version, IMO!

Pieces of Pi

March 28, 2011
13:36	shortcircuit	http://rosettacode.org/wiki/Pi
13:37	* shortcircuit	things Perl6 should implement that task as a sequence.
13:37	shortcircuit	*thinks

Well, can’t turn down a challenge like that.

The task is to create a program to continually calculate and output the next digit of pi. The program should continue forever (until it is aborted by the user) calculating and outputting each digit in succession. The output should be a decimal sequence beginning 3.14159265 …

First, research.

Ah, not only is it a cool project, but it is a chance to rework some Haskell code into Perl 6. Shiny!

So, it looks like the core of their approach is this function:

> stream :: (b->c) -> (b->c->Bool) -> (b->c->b) -> (b->a->b) -> 
>           b -> [a] -> [c] 
> stream next safe prod cons z (x:xs) 
>   = if safe z y 
>     then y : stream next safe prod cons (prod z y) (x:xs) 
>     else stream next safe prod cons (cons z x) xs 
>       where y = next z 

You know, some algorithms make more sense when they are made recursive. Maybe it’s just because I’m not a hardcore functional programmer, but in this case recursion mostly serves to obfuscate the code. Here’s my attempt to rewrite it in Perl 6:

sub stream(&next, &safe, &prod, &cons, $z is copy, @x is copy) {
    gather loop {
        my $y = next($z);
        if safe($z, $y) {
            take $y;
            $z = prod($z, $y);
        } else {
            $z = cons($z, @x.shift)

To my mind, even though this is longer, it’s vastly clearer. There is one major problem, however: it doesn’t actually work. That @x is copy there is not guaranteed to work when you pass in an infinite list; and the entire point of this routine is to pass in infinite lists! I don’t know of a graceful way to handle this in Perl 6. This next works:

sub stream(&next, &safe, &prod, &cons, $z is copy, @x) {
    my $x-list = @x.iterator.list;
    gather loop {
        my $y = next($z);
        if safe($z, $y) {
            take $y;
            $z = prod($z, $y);
        } else {
            $z = cons($z, $x-list.shift)

It feels like there should be some more automatic way to do that in p6. Maybe something like “@x is ro-iterator”? (Yeah, that’s an ugly name.)

Update: Ugly &s deleted from the sub bodies at moritz_++’s suggestion.


March 24, 2011

I’m afraid I don’t have any great Perl 6 news to report. But I wanted to share this brilliant non-Perl programming post on testing.

Also, if anyone wonders what I’m up to, a lot of my current $work is being described on a new blog of mine, and a lot of my current play is at another. But mostly for the last week I’ve been a dad, as my wife is recovering from surgery and is physically unable to lift our 2.5-year-old.

More on masak’s p5

March 9, 2011

So, I left the benchmark code running with a couple more strings last night before I went to bed. Here are the results I got:

p5-moritz.pl: 5.3499312
p5-fox.pl: 5.5914252
p5-colomon.pl: 9.0103538
p5-util.pl: 12.8984396
p5-matthias.pl: 26.1100836
There were nobles, who made war against each other; there was the king,
who made war against the cardinal; there was Spain, which made war against the king.
p5-fox.pl: 9.0724477
p5-moritz.pl: 10.4737434
p5-colomon.pl: 12.4736412
p5-matthias.pl: 25.980024
p5-util.pl: 32.2906812

(That’s the same test as the previous post, but I re-ran it again.)

In those times panics were common, and few days passed without some city or other registering in its archives an event of this kind.
There were nobles, who made war against each other; there was the king, who made war against the cardinal; there was Spain, which made war against the king.
p5-colomon.pl: 20.1296963
p5-fox.pl: 29.2432603
p5-moritz.pl: 33.6600319
p5-util.pl: 119.7950674
p5-matthias.pl: 139.7882468

If people have suggestions for more strings to run, I’d be happy to test them, but bear in mind that the benchmark script runs each test 10 times and averages the results, so it may take some time to get results.

UPDATE: I’ve added three more tests. I also added a new and somewhat novel routine from mberends which wasn’t part of masak’s competition. First, mberends’s results on the above tests:

dna, p5-mberends.pl: 25.846477
dumas-1, p5-mberends.pl: 52.3829865
dumas-2, p5-mberends.pl: 208.6451586
p5-moritz.pl: 15.4340812
p5-colomon.pl,: 16.0729917
p5-fox.pl: 19.6679714
p5-util.pl: 49.7325891
p5-mberends.pl: 107.7859086
p5-matthias.pl: 174.7841558
p5-colomon.pl: 40.274496
p5-moritz.pl: 64.2675568
p5-fox.pl: 88.9493287

(Apologies that I didn’t time all the codes for that one, but it was very obvious the others were much slower, and I needed the Linux box for $work.)

In those times panics were common, and few days passed without some city or other registering in its archives an event of this kind. There were nobles, who made war against each other; there was the king, who made war against the cardinal; there was Spain, which made war against the king.
Then, in addition to these concealed or public, secret or open wars, there were robbers, mendicants, Huguenots, wolves, and scoundrels, who made war upon everybody. The citizens always took up arms readily against thieves, wolves or scoundrels, often against nobles or Huguenots, sometimes against the king, but never against cardinal or Spain. It resulted, then, from this habit that on the said first Monday of April, 1625, the citizens, on hearing the clamor, and seeing neither the red-and-yellow standard nor the livery of the Duc de Richelieu, rushed toward the hostel of the Jolly Miller. When arrived there, the cause of the hubbub was apparent to all.
p5-colomon.pl: 167.9458936
p5-fox.pl: 467.1575598
p5-moritz.pl: 858.1724499
p5-util.pl: 2079.24607
p5-mberends.pl: 3054.501836
p5-matthias.pl: 4872.900568

Benchmarking p5

March 9, 2011

So, matthias’s results in masak’s p5 challenge had me utterly baffled. Anyone with much Rakudo programming experience knows that something like

gather for @strA[$i..*-1] Z @strB[$j..*-1] -> $a, $b

in your inner loop is DOOM speed-wise. How could that possibly be a world-beater?

Well, the answer is that it all depends on the strings you are searching for common substrings. I’d asked masak about what test he was using a few days ago (trying to figure out the weird performance of my code) and gotten his routine:

    my $n = 160;
    my $a = <0 1 2 3 4>.roll($n).join;
    my $b = <5 6 7 8 9>.roll($n).join;

    # "N" case is testing this $a versus this $b

    my $ra = (0..$n - 10).pick;
    my $rb = (0..$n - 10).pick;
    $a = $a.substr(0, $ra) ~ 'flirtation' ~ $a.substr($ra + 10);
    $b = $b.substr(0, $rb) ~ 'flirtation' ~ $b.substr($rb + 10);

    # "Y" case is testing this $a versus this $b

If you look at that, what should jump out at you is, in the “N” case, the two strings have no characters at all in common with each other! matthias’s code is nearly perfectly optimized for this case, as the slow inner loop will never execute. The “Y” case will be nearly as good for it.

So, what happens if we use normal text to seed this? Luckily I just got some benchmarking code going. So if the strings are

There were nobles, who made war against each other; there was the king,
who made war against the cardinal; there was Spain, which made war against the king.

what are the timings like? Sorted from best to worse:

p5-fox.pl: 9.0981144
p5-moritz.pl: 10.4469054
p5-colomon.pl: 12.5152087
p5-matthias.pl: 25.9988418
p5-util.pl: 32.1561981

More strings and timings coming soon….

PS I should add I love the style of matthias’s code. It combines a lot of elegance with a clever Perl-ish and p6-ish approach to the problem. If I hadn’t stumbled across the suffix tree solution, I suspect my code would have looked a lot like an uglier version of his.

First Benchmark Results

March 6, 2011

I’ve gotten smash++’s Rakudo benchmarks running again, and I’m trying to get a full set of benchmark graphs generated. While moritz and I are debugging the graphing program, here are the first results. (Assuming I can get them to display here!)

Okay, two days later, I’ve finally got inkscape running again to generate PNG versions of the graphs. And with any luck, here they are:

Update: I forgot people who weren’t familiar with the original benchmarks would be reading this. All graphs are labeled with the number of seconds on the left hand side. Times given are the average execution time (including Rakudo start-up and compiling the benchmark code) over ten runs, on my reasonably fast 64-bit Linux machine. The benchmark scripts are at https://github.com/perl6/bench-scripts, and are a fairly random bunch at the moment. Contributions of new scripts would be gladly accepted.

As for the meaning of the results, I haven’t had time to analyze them in any depth. latest-rakudo seems pretty consistently a bit slower than the latest Rakudo Star, which makes me wonder if there are some different flags set in the build process for each. (I used the default build process for both R* and latest-rakudo, but they are two different build methods and might have some different compiler flags.) It also looks like a couple of the benchmarks have backslid very badly in the last two months. Certainly a full investigation is warranted.

Simplifying Rationals

February 15, 2011

So, in my last post I came up with this code for adding two FatRats:

multi sub infix:<FR+>(Math::FatRat $a, Math::FatRat $b) is export(:DEFAULT) {
    my $gcd = gcd($a.denominator, $b.denominator);
    Math::FatRat.new($a.numerator * ($b.denominator div $gcd) + $b.numerator * ($a.denominator div $gcd),
                     ($a.denominator div $gcd) * $b.denominator);

Some of you may have noticed something funny here. Why is there a GCD calculation?

Well, this code was a cut-n-paste from the Rat code inside Rakudo. And since Rakudo has only finite Ints, it seemed like a good plan to use the GCD of the two Rats’ denominators to make the numbers as small as possible in the calculation, hoping to keep things inside the range of 32-bit integers.

So, since Math::BigInt has no range limitation and Math::FatRat.new is going to simplify the fraction anyway, it seems like calculating the GCD here is a waste of time. Or is it? I’m leaning toward it being a waste of time, but I hesitate to say that for sure without doing some timing tests, which I don’t want to go into now. Instead, I’m going to open an entirely different cartoon of worms in this same neighborhood.

Because the spec is actually kind of vague on when and even if you should — or can! — simplify fractions. I’ve found three areas in S02 which touch on this:

The limitation on Rat values is intended to be enforced only on user-visible types. Intermediate values used internally in calculation the values of Rat operators may exceed this precision, or represent negative denominators. That is, the temporaries used in calculating the new numerator and denominator are (at least in the abstract) of Int type. After a new numerator and denominator are determined, any sign is forced to be represented only by the numerator. Then if the denominator exceeds the storage size of the unsigned integer used, the fraction is reduced via gcd. If the resulting denominator is still larger than the storage size, then and only then may the precision be reduced to fit into a Rat or Num.

Rat addition and subtraction should attempt to preserve the denominator of the more precise argument if that denominator is an integral multiple of the less precise denominator. That is, in practical terms, adding a column of dollars and cents should generally end up with a result that has a denominator of 100, even if values like 42 and 3.5 were added in. With other operators, this guarantee cannot be made; in such cases, the user should probably be explicitly rounding to a particular denominator anyway.

Although most rational implementations normalize or “reduce” fractions to their smallest representation immediately through a gcd algorithm, Perl allows a rational datatype to do so lazily at need, such as whenever the denominator would run out of precision, but avoid the overhead otherwise. Hence, if you are adding a bunch of Rats that represent, say, dollars and cents, the denominator may stay 100 the entire way through. The .nu and .de methods will return these unreduced values. You can use $rat.=norm to normalize the fraction. (This also forces the sign on the denominator to be positive.) The .perl method will produce a decimal number if the denominator is a power of 10, or normalizable to a power of 10 (that is, having factors of only 2 and 5 (and -1)). Otherwise it will normalize and return a rational literal of the form -47/3.

I actually find these paragraphs somewhat bewildering. Let me try to sum up the points as I’m seeing them, in vaguely reversed order.

1. Rats are “allowed” to be “lazy” and never simplify, but must always return the simplified version when .perl is called. (I put lazy in quotes because Rats are immutable, so in fact the Rat object in question will always be unsimplified.) Which means that Rat.perl doesn’t provide a way of actually getting at actual value stored in a Rat; all you get is another Rat which has the same numeric value. That’s…. weird.

2. “Rat addition and subtraction should attempt to preserve the denominator…” That sounds like Rats are required to be lazy, at least in some circumstances.

3. Note that those circumstances are a bit weird. It’s pretty easy to find sets of four Rats which can have a different denominator based on the order you add them.

4. Note also that this property goes away as soon as .perl is involved.

It feels to me like there are two distinct ideas here, at odds with each other. In one paragraph, the developer is allowed to break some fundamental assumptions of Perl 6 to make rational math more efficient. In another, the developer is required to try to bend how rationals work make rational math more friendly or something. They are at odds because the efficiency version makes .nu and .de into something you probably don’t want to look at, whereas the only way to take advantage of the “friendly” version is to look at those exact same values! Not to mention the extra work needed to make the friendly version work probably kills most or all of the hypothetical efficiency improvements available in the other.

Why do I say that? Let’s look at implementations of the two separate approaches. Here’s the efficient one (remember that in this version, FatRat.new doesn’t do anything but store the numerator and denominator it is given):

multi sub infix:<FR+>(Math::FatRat $a, Math::FatRat $b) is export(:DEFAULT) {
    Math::FatRat.new($a.numerator * $b.denominator + $b.numerator * $a.denominator,
                     $a.denominator * $b.denominator);

You might want to consider adding an if check there to see if the two denominators are equal, as that can save you three multiplications. Whether or not that would help probably depends on the size of the Ints and whether or not you’re adding a lot of rationals with the same denominator.

On the other hand, the “friendly” version would have to be coded something like this:

multi sub infix:<FR+>(Math::FatRat $a, Math::FatRat $b) is export(:DEFAULT) {
    if ($a.denominator %% $b.denominator) {
        Math::FatRat.new($a.numerator + $b.numerator * $a.denominator div $b.denominator,
    } elsif ($b.denominator %% $a.denominator) {
        Math::FatRat.new($a.numerator * $b.denominator div $a.denminator + $b.numerator,
    } else {
        Math::FatRat.new($a.numerator * $b.denominator + $b.numerator * $a.denominator,
                     $a.denominator * $b.denominator);

Note that this version adds two useless (and probably relatively expensive!) is-divisible-by tests to what I would consider to be the “normal” case. Even in the best-case scenario of hitting the first special case, you’ve just replaced an operation which required three multiplications and an addition with one that requires a divisible-by test, a division, a multiplication, and an addition. Unless that prevents the size of the denominator from growing very huge, that’s probably a pessimization. (And note that if you’re doing Rat arithmetic, it’s already required if the denominator is too big to fit in an int, the fraction is simplified if possible.)

On the other hand, who am I to argue against it on the basis of efficiency? My instinct is very strongly to always simplify. Though I note that a GMP reference says “In general, cancelling factors every time is the best approach since it minimizes the sizes for subsequent operations.”

What do I think? Well, I completely fail to understand the usage case for the “friendly” paragraph. If it’s just an optimization, then it shouldn’t be “required”. If it’s intended to make things easier on the user, it’s remarkably fragile and hard to use. Suppose you think you’re adding adding dollars and cents, and want to get the result $a in terms of the number of cents (about the only practical use case I can think of). Then you need to do something like

    $a += 0/100; # make sure the result is at least in hundredths.
    fail if $a.denominator != 100; # make sure the result really is hundredths

Wouldn’t it make more sense to have a .numerator-if-denominator-was($n) method, which returns what the numerator would be if the denominator was $n? It would look something like this:

    method numerator-if-denominator-was(Int $n) {
        fail unless 100 %% $.denominator;
        $.numerator * 100 div $.denominator;

It seems like this would cleaner for the user, while not putting any requirements on the internal structure of the rational type at all.

(Of course, if you really want to deal with dollars and cents, you should probably be creating a type specifically to handle that! It could easily be substantially more efficient than a generic rational class.)

My initial inclination with the other question is that I’m inclined to go one of two ways. If we really want to allow Rats which have not been simplified, then we should go whole hog and even .perl should return the unsimplified value. If you want it simplified, call .norm.perl.

Or perhaps Rats are always simplified if possible, but we add language to the spec which allows a Perl 6 implementation to maintain the unsimplified version in the midst of a series of calculations for efficiency purposes. Why would such language be needed? My notion (which I haven’t pinned down with hard numbers yet) is that there are some operations which would spill out of a Rat if you did them step by step, but might simplify to a value which can be held in a Rat once all the operations were done and the result simplified. The language would mean the optimized code wouldn’t have to track each operation to see if the partial solution would fit in a Rat.

Hmmm. I think more research and benchmarking are probably called for…