mjd’s magic z

March 13, 2012

So, I plunged ahead and implemented mjd’s z function in Perl 6. This is the first simple step to full continued fraction arithmetic, allowing addition with a rational, multiplication by a rational, and taking the reciprocal of a continued fraction. (And oh! I just accidentally discovered that I’d missed a bunch of slides at the end of mjd’s talk which discuss implementing the full z function from HAKMEM!)

The implementation turned out to be pretty straightforward, with only the edge cases causing trouble. z maintains its state using four variables, $a, $b, $c, and $d. There are two basic operations which modify this state: input, which takes the next value from the continued fraction provided, and output, which outputs the next value of the continued fraction result. How do you know which to do? Easy! If $a div $c == $b div $d, then we’ve determined the next value which can be output, and should output it. If not, then we need to input another value.

The tricky bits: first, it’s completely legal for $c and/or $d to be 0. That causes a division by zero error in Niecza. So I added a check for each, setting the result of the division to Inf when it needs to be. The second thing is what happens when the given continued fraction runs out of values. In that case, the next value is implicitly Inf. That’s well and good, but mjd somehow performs magic with it. I don’t understand what the heck he was doing there, but I can imitate the result. Here it all is in one function:

sub z($a is copy, $b is copy, $c is copy, $d is copy, @x) {
    gather loop {
        my $a-div-c = $c ?? $a div $c !! Inf;
        my $b-div-d = $d ?? $b div $d !! Inf;
        last if $a-div-c == Inf && $b-div-d == Inf;
        if $a-div-c == $b-div-d {
            my $n = $a-div-c;
            ($a, $b, $c, $d) = ($c, $d, $a - $c * $n, $b - $d * $n);
            take $n;
        } else {
            if @x {
                my $p = @x.shift;
                ($a, $b, $c, $d) = ($b, $a + $b * $p, $d, $c + $d * $p);
            } else {
                ($a, $b, $c, $d) = ($b, $b, $d, $d); # WHY????

And there you have it! I’m still a bit worried about the edge cases (I haven’t tested negative numbers yet, for instance), but I think it clearly points the way onto the full HAKMEM z function.

Continued Fractions

March 13, 2012

Four months ago I created Math::ContinuedFractions on Github. But I didn’t actually post any code to it. But I’ve finally started, pushing t/00-experiments.t. It’s very simple so far, but that’s because I’m slowly feeling my way through this unfamiliar subject.

I tried to link the Wikipedia equation graphic for a continued fraction, but that just lost me the first draft of this post, so instead I’ll link to the Wikipedia page. It is the source of our first algorithm, for converting a number to a continued fraction:

sub make-continued-fraction (Real $x is copy) {
    gather loop {
        my $a = $x.floor;
        take $a;
        $x = $x - $a;
        last if $x == 0;
        $x = 1 / $x;

This returns the components of the continued fraction as a lazy list. It passes the first few simple tests I’ve thrown at it, but I’m a bit worried about $x = 1 / $x. Should this maybe be FatRat math instead of normal math? I’m not clear on the implications one way or the other.

Unfortunately, the Wikipedia page doesn’t shed any light on how to do basic arithmetic on continued fractions. HAKMEM gives a complete overview on basic continued fraction arithmetic, but I find it quite hard going. MJD has a nifty talk on the subject, including some algorithms I hope to implement soon. Unfortunately, he doesn’t get into the tricky stuff.

Luther Tychonievich has a nice blog post on the subject which brings up a problem I have not seen mentioned elsewhere. Some pretty basic operations may not work because they need every (infinite) input term to determine the first output term. His example is squaring the square root of 2. As I understand it, the problem is that the first term (the integer part) flits back and forth between 1 and 2 as each addition term is fed to the code. It can never settle down on one or the other.

Incidentally, my code fails miserably in Rakudo, but works in Niecza. Unfortunately for this project, Niecza cannot currently handle overloading the arithmetic operators for a new type, which means I’ll have to figure out how to get things working on Rakudo at some point.

Modules and Perl 6s

January 29, 2012

I’ve got my fork of Panda working on Niecza. There are two major issues which remain to be resolved, and they are both related. Panda assumes that the Perl 6 executable is named perl6, and modules should be installed to ~/.perl6. That’s great if you have just one Perl 6 on your system.

But in this crazy modern world, I really want both Rakudo and Niecza on my system. And so assuming that everything can be named “perl6” is a big problem.

I’m hoping the community can put their heads together and figure out a solution. I guess the obvious one (which would only need sorear to buy in) would be installing Niecza modules to ~/.niecza, and creating a “real” niecza executable. (It might well be as simple as a one-line script “mono path-to-Niecza.exe”.)

Good idea or bad idea?

Last Piece of Pi

January 15, 2012

So, last spring I did a series of posts on making a Pi spigot. Unfortunately, the project foundered a bit, as it turned out that only Rakudo had the subsignatures needed to make the code pretty, but only Niecza had the big integers needed to make the code work.

Fast forward to today. Now both Rakudo and Niecza can handle the best version of the code with no problem!

Let me emphasize that again. A year ago, neither implementation had those two features. Eight months ago, each had one. Today both have both. What’s more, Rakudo seems to have made some pretty impressive speed improvements.

I think the awkward situation with Rakudo Star has helped obscure the fantastic good news in Perl 6. Right now we have two distinct Perl 6 implementations, each of which is markedly better than anything available a year ago. While there are still some rough edges, both implementations are making visible progress every single day.

Okay, enough cheerleading. Time to finish the Pi story. It turned out there was one last complication. I had never actually tried to use the code to generate an infinite stream of Pi. I always stopped at 40 digits. Confident that the algorithm must work, I tweaked the output to just continue generating digits until you stopped the program with control-C. And I ran it using Niecza, and it printed


And then it just sat there, calculating. After a few minutes I stopped it and tried it in Rakudo, and got the exact same result.

Well, after adding a few debugging says here and there, I found the problem. The extr sub was returning NaN. It turns out that both compilers have issues with dividing huge numbers. Here’s the division operation that was causing the problem:

70777782797813415527343750 / 7640045443

Turns out the bottom number is bigger than can be presented by a Num, so the division operation ends up becoming N / Inf, which is NaN. This is a bit obscured in Rakudo because the division operation actually returns a Rat (which should be illegal according to the spec!), but then .floor is called on the result, which tries to convert to Num before doing the floor operation.

This opens several questions, like: Should Perl 6’s Int / Int operator be modified to try to cope with this sort of case?

But for the Pi problem, the good news is this: every time the extr sub is called, its result is fed to .floor. That means we simply needed to replace / with div and get rid of the .floors.

And with that, the code easily produces 2,000+ digits of Pi using either Rakudo or Niecza!

November 11, 2011

chromatic has come up with a lovely new metaphor to replace the idea of technical debt: “technical insurance”. I love this sentence: “Given the premiums for your team, you then have to decide whether or not to buy a technical insurance policy.”

It hits at something that has bugged me about testing enthusiasts — the notion that you should always have tests. Because there is a balance that needs to be found — the cost of the tests versus the cost of not having the tests.

I’ve talking about this with respect to the ABC module before. Properly testing the abc2ly script would require constructing an OCR system for reading sheet music. That would be very expensive technical insurance indeed.

And what would it be protecting against? abc2ly is a open source tool for creating sheet music. It makes no money for me. As far as I know, I am the only active user. The only real risk I can see is if I am printing a long musical document, and a change I’ve made breaks something subtle that only affects, say, one thing on page 14, so I don’t notice the problem in a timely fashion. So basically, there’s a potential for wasting some paper.

In other words, in this case technical insurance would be expensive, and have very limited benefits. I’d be a fool to “buy” it.

On the other hand, the Perl 6 spectests are a fantastic technical insurance policy. For the most part they are low cost (easy to write) AND they are shared across multiple implementations of Perl 6. On the flip side, they help protect against subtle flaws in the implementations, which have the potential of causing problems for every person using Perl 6. We’d have to be crazy to do without this policy.

Fixing Tags

November 9, 2011

So, in my last post I identified 3906 files in my MP3 collection with missing tags. This time I set out to fix some of them.

So, first I went through the list I generated with the last script and singled out all 2294 files which used a standard pattern of Artist / Album / Track Number - Track Name. Then I wrote this script:

my $for-real = Bool::True;

constant $TAGLIB  = "taglib-sharp,  Version=, Culture=neutral, PublicKeyToken=db62eba44689b5b0";
constant TagLib-File    = CLR::("TagLib.File,$TAGLIB");
constant String-Array   = CLR::("System.String[]");

for lines() -> $filename {
    my @path-parts = $filename.split('/').map(&Scrub);
    my $number-and-title = @path-parts.pop;
    next unless $number-and-title ~~ m/(\d+) \- (.*) .mp3/;
    my $track-number = ~$0;
    my $title = ~$1;
    my $album = @path-parts.pop;
    my $artist = @path-parts.pop;
    say "$artist: $album: $title (track $track-number)";

    if $for-real {
        my $file;
        try {
            $file = TagLib-File.Create($filename);
            CATCH { say "Error reading $filename" }

        $file.Tag.Track = $track-number.Int;
        $file.Tag.Album = $album;
        $file.Tag.Title = $title;
        $file.Tag.Performers = MakeStringArray($artist);
        try {
            CATCH { say "Error saving changes to $filename" }

sub Scrub($a) {
    $a.subst('_', ' ', :global);

sub MakeStringArray(Str $a) {
    my $sa = String-Array.new(1);
    $sa.Set(0, $a);

For the main loop, the first half uses standard Perl techniques to extract the artist, album, and track info from the path. The second half sets the tags. Opening the file is the same as last time, and then setting Track, Album, and Title is as simple as could be. The Performers tag is a bit tricky, because it’s a string array (the others are simple strings or integers) and Niecza doesn’t know how to do the coercion automatically. MakeStringArray gets the job done nicely.

So, if you’ve done this sort of thing in Perl 5 using the MP3 CPAN modules, there’s nothing at all revolutionary about this code. But it feels really good to be able to do it with Perl 6!

Examining MP3 Tags

November 6, 2011

I’ve been playing around with Niecza’s ability to handle CLR libraries. It’s actually kind of intoxicating; it’s the closest thing yet to having a CPAN for Perl 6. So I decided to see what I could do with the TagLib# library for dealing with media file tags.

Now, I’ve got an MP3 library with 23564 MP3 files in it, the majority of which were created by older ripping files that didn’t do anything with the ID tags. Most of those have been updated to include tags, but every now and then I add one of the old directories to iTunes and get a bunch of “Unknown Artist” / “Unknown Album” tracks.

So I thought a nice first project would be figuring out which of the tracks was correct. The first thing to do was to get TagLib# properly installed on my MacBook Pro. make install didn’t add the DLL to the main GAC; I ended up installing it there by hand, which was trivially easy once I knew what to do:

sudo gacutil -i taglib-sharp.dll
sudo gacutil -i policy.2.0.taglib-sharp.dll

Once I had that done, I experimented with it for a bit, and ended up with this script:

constant $TAGLIB  = "taglib-sharp,  Version=, Culture=neutral, PublicKeyToken=db62eba44689b5b0";
constant TagLib-File    = CLR::("TagLib.File,$TAGLIB");

for lines() -> $filename {
    try {
        my $file = TagLib-File.Create($filename);
        unless $file.Tag.JoinedPerformers ~~ m/\S/ && $file.Tag.Title ~~ m/\S/ {
            say $filename;
        CATCH { say "Error reading $filename" }

The first line specifies the exact assembly we want to use; you can get the details from gacutil -l. The next line effectively imports the TagLib::File class into Niecza. I get my filenames from stdin, as that allows me to use find to generate the list of MP3 files.

This was my first use of exception handling in Perl 6. I needed it because TagLib-File.Create throws an exception when it isn’t happy with the MP3 file. When it is happy with it, $file is an object of type CLR::TagLib::Mpeg::AudioFile. $file.Tag.JoinedPerformers gives the list of performers (AKA artists) as a single string; $file.Tag.Title gives the title as a string. Unless we find a valid non-space character in both of them, we flag the file by printing it out.

Really, the only way it could be significantly simpler than this would be if the constant TagLib-File line were unnecessary!

End result: I have a list of 3906 files it flagged, 77 of which were read errors.

My next step is to write some code which translates the filenames (which are mostly of the form /Volumes/colomon/Albums/Dervish/Live_in_Palma/04-Slow_Reels.mp3) into artist, album, and track name fields, and then set those tags. Based on my initial experiments, I think it’s is going to be incredibly easy…

Ease of FatRat construction

October 18, 2011

So, on #perl6 today tried using the numeric literal .3333333333333333333333333333333. (Warning: exact number of 3‘s may not match original example.) By the spec (as I understand it), this is a Num, because a Rat isn’t accurate enough to represent it. (Not that a Num is, mind you!)

And that got me to thinking: What if you really wanted a FatRat, so you actually got that exact number? Well, if you’re using Niecza (the only p6 to implement FatRat so far), the answer is FatRat.new(3333333333333333333333333333333, 10000000000000000000000000000000). IMO, that’s ridiculously awkward.

The spec may imply you can do it with ".3333333333333333333333333333333".FatRat. That at least avoids the problem of counting the zeros, but it’s still on the ugly side. Likewise FatRat.new(".3333333333333333333333333333333") is awkward. Still, we should certainly support at least one of these options.

I would like to propose again adding an F suffix to indicate a numeric literal should be a FatRat. I don’t think this is something that can reasonably be done with a postfix operator, because if you treat .3333333333333333333333333333333 like a normal numeric value and then try to FatRat it, you will lose the precision you want.

Just as a quick comparison, here’s a bit of the old endless pi code using the FatRat constructor:

sub unit() { LFT.new(q => FatRat.new(1, 1),
                     r => FatRat.new(0, 1),
                     s => FatRat.new(0, 1),
                     t => FatRat.new(1, 1)); }

I’m proposing we should be able to write that as

sub unit() { LFT.new(q => 1F,
                     r => 0F,
                     s => 0F,
                     t => 1F); }

Much shorter and much clearer. I think that’s a big win.

(Note: I’m in no way particularly attached to the letter “F” for this, that was just the first thing that came to mind.)

Complex Issues

August 31, 2011

Sorry for the long silence here, it’s been a busy summer with far too little Perl 6. But I did squeeze in some work on trig, both on nom and niecza. And I ran into a very interesting issue.

My local copy of niecza has S32-trig/sin.t almost working. A few needed skips, but all the core numeric types work. Except…

> is_approx(asin(0.785398163404734 + 2i), 0.341338918259482 + 1.49709293866352i
# got:      2.80025373533031-1.49709293866352i
# expected: 0.341338918259482+1.49709293866352i

niecza> asin(0.785398163404734 + 2i)

rakudo> asin(0.785398163404734 + 2i)
0.341338918259481 + 1.49709293866352i

Woah, what’s up with that? Well, it turns out both answers are right in some sense:

niecza> sin(asin(0.785398163404734 + 2i))

rakudo> sin(asin(0.785398163404734 + 2i))
0.785398163404734 + 2i

The thing here is that sin is periodic; there are an infinite number of complex numbers it maps to the same result value. That means when you call asin, there are an infinite number of possible results for each input value, and you must somehow choose one of them.

But let’s take a step back from that and look at why I got different results, because I used the exact same formula for asin in both Rakudo and Niecza. That formula is -1i * log(($x)i + sqrt(1 - $x * $x)). Let’s look at the sqrt first:

niecza> my $x = 0.785398163404734 + 2i; sqrt(1 - $x * $x)

rakudo> my $x = 0.785398163404734 + 2i; sqrt(1 - $x * $x)
2.21086930051619 - 0.710488099157523i

As you can see, one answer is the negative of the other. Of course, when you square the results, that additional factor of -1 just goes away, so these are both valid results.

So this leads me to two questions:
1) Should we define one of these two answers as being correct, as far as Perl 6 is concerned? (Or should they both be considered valid results?)

2) If so, which one? And how do we coherently specify that branch?

I thought at first it might be as simple as saying “The branch where the complex result of sqrt for complex numbers with an imaginary value of 0 agrees with the real sqrt result.” But in fact both Rakudo and Niecza already seem to agree for the sqrts of real-valued Complex numbers.

Anyone else out there have a notion?

Perl 6 Numerics

June 29, 2011

I’m giving a talk on Numerics in Perl in two hours at YAPC::NA 2011. I’m posting the slides here. If you have questions, comments, or corrections, please comment on this post.

If you’re interested in general Perl 6 info, you can check out my previous talk.