Math::BigInt

I suspect everyone who plays around with numerical problems on Perl 6 has been as frustrated as I am not to have arbitrary precision integers in Rakudo. As a sort of first step toward getting them, I’ve thrown together the beginnings of a Math::BigInt module using the BIgDigit multiple-precision arithmetic library and Zavolaj.

(Why BigDigit instead of GMP? Two simple reasons. First, building GMP failed on my main development machine (OS X 10.5). GMP’s suggestion for fixing the problem involved getting a new main development environment for the machine, which seemed like a major pain in the neck, and didn’t impress me with GMP’s portability. Second, GMP’s main multiple-precision type is a C struct you store on the stack. That means using it with Zavolaj would have required creating wrappers for every GMP function I wanted to use. Yuck.)

Without further ado, here’s a few examples:

> my @fib := Math::BigInt.new(1), Math::BigInt.new(1), *+* ... *;
> say @fib[200];
453973694165307953197296969697410619233826

> my @powers-of-two := Math::BigInt.new(1), * * 2 ... *;
> say @powers-of-two[80];
1208925819614629174706176

From which you can quickly conclude two things: Math::BigInts easily worked with the sequence operator. And creating a Math::BigInt is ugly and awkward. Once again, blogging has forced me to confront the limitations of my code.

So… if I had an operator to create a BigInt, the code would look a lot nicer. What’s more, you could avoid the term Math::BigInt actually appearing in your code, making it very easy to port to Perl 6s which actually support arbitrary precision Ints. Clearly then, I need to figure out a symbol to use…

Zavolaj!

I’d been thinking for a while that I should look at adding the extended standard math library functions to Perl 6. And I’d been thinking I should try to do something with Zavolaj, the Perl 6 native call interface. And gee, those two things work perfectly together!

use NativeCall;

sub log1p(Num $x) returns Num is native("libm") { ... }
sub expm1(Num $x) returns Num is native("libm") { ... }
sub erf(Num $x) returns Num is native("libm") { ... }
sub erfc(Num $x) returns Num is native("libm") { ... }
sub tgamma(Num $x) returns Num is native("libm") { ... }
sub lgamma(Num $x) returns Num is native("libm") { ... }

Okay, that was really, really easy to do. I am officially impressed by Zavolaj.

I was thinking the next step is to make a module of this, which I think should be pretty easy. But I’m blocked on a name for it. Any suggestions?

A Word for the Slow

So, my solution for Masak’s p1 has the distinction of being by far the least efficient working solution. Which is a shame, because I think this one was my favorite solution of the contest. It may be slow, and I’d never advise anyone to use it for anything practical (given the other, much more efficient solutions), but in my opinion it’s a charming piece of code.

The key organizational piece for my solution is the Ordering class. (BTW, I didn’t like that name at all, it’s probably my least favorite part of my solution. My theory is Masak was too awestruck by my inefficiency to quibble about the name.) I was striking about for how to represent a series of matrix multiplications, and hit on the idea of using a very simple stack-based language to do it. The language has two operations: an Int represents putting the matrix with that index on the stack. The string "*" represents popping the top two matrices on the stack, multiplying them, and pushing the result back on the stack. Here’s the code for making that happen while tracking the total number of multiplications:

    method calculate-multiplications(@matrices) {
        my @stack;
        my $total-multiplications = 0;
        for @.ops {
            when "*" { 
                my $a = @stack.pop; 
                my $b = @stack.pop;
                my ($multiplications, $matrix) = multiply($a, $b);
                $total-multiplications += $multiplications;
                @stack.push($matrix);
            }
            when Int {
                @stack.push(@matrices[$_]);
            }
        }

        $total-multiplications;
    }

I’m an old Forth programmer from way back, and I can’t begin to say how much I love how easy p6 makes it to implement a simple stack machine!

Getting the string version of this is equally easy:

    method Str() {
        my @stack;
        for @.ops {
            when "*" { 
                my $a = @stack.pop; 
                my $b = @stack.pop;
                @stack.push("($a$b)");
            }
            when Int {
                @stack.push("A{$_ + 1}");
            }
        }

        @stack[*-1];
    }

This time instead of a stack of Pairs (for the matrix size), we have a stack of Str representing each sub-matrix’s name. At the end we pop the last thing on the stack, and it’s the string representing the entire multiplication. And by making this Ordering.Str, any time you print an Ordering you get this nice string form — handy both for the desired output of the program and for debugging.

I won’t comment on the guts of the generate-orderings function, which is heavily borrowed from HOP via List::Utils. Just note that given the number of matrices, it lazily generates all the possible permutations — both the source of my code’s elegance and its extreme inefficiency.

Oonce you’ve got the array @matrices set up, calculating and reporting the best ordering (very slowly!) is as simple as

say generate-orderings(+@matrices).min(*.calculate-multiplications(@matrices));

(Note that I split the main body of this off into a function in the actual code, to make it easier to test internally.)

So clearly, I badly need to study more dynamic programming. But at the same time, I think there may be useful bits in my code that can be put to better use somewhere else.

A Piece of Pi

Between masak’s Perl 6 contest and the Perl 6 Advent calendar, I haven’t done a lot of blogging here lately. Apologies!

So, a discussion on #perl6 the other day got me thinking about another possible interesting use of the Real role. Here’s the basic class:

class TimesPi does Real {
    has Real $.x;
    
    method new($x) { self.bless(*, :$x); }
    method Bridge() { $.x * pi; }
    method Str() { $.x ~ "π"; }
    method perl() { "({ $.x })π"; }
}

sub postfix:<π>($x) {
    TimesPi.new($x);
}

So, that’s simple enough, right? TimesPi stores a Real number $.x internally, and the value it represents is that number times pi. There’s a postfix operator π to make it really easy to construct these numbers. Because we’ve defined a Bridge method, this class has access to all the normal methods and operators of Real. Still, as presented above it is pretty useless, but defining some operators hints at a useful purposes for this class.

multi sub infix:<+>(TimesPi $lhs, TimesPi $rhs) {
    TimesPi.new($lhs.x + $rhs.x);
}

multi sub infix:<->(TimesPi $lhs, TimesPi $rhs) {
    TimesPi.new($lhs.x - $rhs.x);
}

multi sub prefix:<->(TimesPi $n) {
    TimesPi.new(- $n.x);
}

multi sub infix:<*>(TimesPi $lhs, Real $rhs) {
    TimesPi.new($lhs.x * $rhs);
}

multi sub infix:<*>(Real $lhs, TimesPi $rhs) {
    TimesPi.new($lhs * $rhs.x);
}

multi sub infix:<*>(TimesPi $lhs, TimesPi $rhs) {
    $lhs.Bridge * $rhs.Bridge;
}

multi sub infix:</>(TimesPi $lhs, Real $rhs) {
    TimesPi.new($lhs.x / $rhs);
}

multi sub infix:</>(TimesPi $lhs, TimesPi $rhs) {
    $lhs.x / $rhs.x;
}

With these operators in place, basic arithmetic involving TimesPi numbers will stay in the TimesPi class when appropriate. For instance, if you add two TimesPi numbers, the result will be a TimesPi. The cool thing about this is that it is as exact $.x values allow, rather than forcing everything to be a floating point calculation of limited accuracy.

We can even take things a step further, using this to perform exact trig calculations:

    multi method sin(TimesPi $x: $base = Radians) {
        return $x.Bridge.sin($base) unless $base == Radians;
        given $x.x {
            when Int { 0; }
            when Rat {
                given $x.x.denominator {
                    when 1 { 0; }
                    when 2 { ($x.x.numerator - 1) %% 4 ?? 1 !! -1 }
                    # could have a case for 6 as well...
                    $x.Bridge.sin; 
                }
            }
            default {
                $x.Bridge.sin;
            }
        }
    }

This checks for cases where we know the exact value of the result, and returns that if it can, otherwise falling back to the standard Real.sin method.

Of course, just when I was feeling like I might be on to something here, I realized that $.x was just the number of degrees in the angle divided by 180. Sigh.