Archive for March, 2012

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.


Follow

Get every new post delivered to your Inbox.