mjd’s magic z

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.

About these ads

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.

%d bloggers like this: