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.

Euler 5

June 23, 2011

(Wrote this post last week, but accidentally posted it to the wrong blog!)

So, Hacker News brought me this post this morning, and I wanted to see how Perl 6 stacked up to Java / Scala. The answer turned out to be pretty poorly speed-wise, at least in my initial attempts. But that’s a story for another post. Right now I’m just interested in finding a good solution to the problem itself: “What is the smallest positive number that is evenly divisible by all of the numbers from 1 to 20?”

So far, this is my favorite, I think:

sub divides-by-all-up-to($a, $b) {
    !(2..$b).grep($a !%% *);
}

my $N = [*] 2, 3, 5, 7, 11, 13, 17, 19;
my @attempts := $N, 2 * $N ... { divides-by-all-up-to($_, 20) };
say @attempts[*-1];

You can easily get this down to a two-liner, but only at the cost of a good bit of clarity, IMO.

So, what does it do? divides-by-all-up-to checks to see if $a is divisible by all the numbers from 2 to $b. This is one of those concepts for which there is very definitely more than one way to do it in Perl 6. For instance,

    ?all($a <<%%<< 2..$b)

or

    [&&} $a X%% 2..$b

and so on. I choose my approach because it is pretty straightforward, reasonably efficient, and works on both Rakudo and Niecza.

$N is my secret weapon. It's the product of all the primes less than 20. All of those primes have to be factors of the answer, so we only consider multiples of $N, saving lots of time. To be precise, it means we call divides-by-all-up-to once for every 4,849,845 times the Scala version calls the equivalent function.

The next step is fun with sequences. We want to consider the multiples of $N, stopping when we hit one for which divides-by-all-up-to is true. This is easily expressible using the sequence operator:

    my @attempts := $N, 2 * $N ... { divides-by-all-up-to($_, 20) };

The initial $N, 2 * $N tells the sequence operator that we want the sequence of multiples of $N; the rest tells it when to stop. (We know that it has to stop eventually because 20! has all the correct properties and is a multiple of $N.) The end result is a sequence whose last value is the number we are looking for. And we finally get that last value with standard Perl 6ish @attempts[*-1].

Here's a much more generalized approach. Not sure if it will work with Niecza or not (I've never figured how to use modules with Niecza) but it does work with Rakudo.

use Math::Prime;

sub factors($n is copy) {
    my @primes := primes;
    my $prime = @primes.shift;
    my %factors;
    while $n > 1 {
        if $n %% $prime {
            %factors{$prime}++;
            $n div= $prime;
        } else {
            $prime = @primes.shift;
        }
    }
    return %factors;
}

sub product(%factors) {
    [*] %factors.map({ .key ** .value });
}

sub least-divisble-by-all(@values) {
    my %common-factors;
    for @values -> $i {
        my %factors = factors($i);
        for %factors -> $prime {
            if %common-factors.exists($prime.key) {
                %common-factors{$prime.key} max= $prime.value;
            } else {
                %common-factors{$prime.key} = $prime.value;
            }
        }
    }
    
    product(%common-factors);
}

This version uses the prime factorization of the numbers to smartly figure out the answer. It's much more flexible and and probably faster for big numbrers, but... it's just too long to make me love it.

The Pi is done

June 4, 2011

When last we left the endless Pi project, it a “simple” matter of getting FatRats to work. Which in this case, meant getting the Math::FatRat module to work. Which in turn meant getting Math::BigInt to work again. Which meant getting Zavolaj to work again. Frankly, I thought I might have a month’s worth of blog posts left in the project.

Enter Niecza. As of the latest release, it has baked-in FatRats. I was worried about it lacking lazy lists, but it turns out it has them too. I ran into four holes in its implementation of Perl 6: No floor function, no bless, no MAIN, and no sub-signatures. Luckily floor’s easy to write, MAIN is easily skipped, using the default new works, and I just reverted TimToady’s sub-signature suggestion. A bit of work, and I had this, a complete and fully functional Perl 6 spigot stream for pi:

sub floor(FatRat $n) {
    my $mod = $n.numerator % $n.denominator;
    ($n.numerator - $mod) div $n.denominator;
}

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);
        }
    }
}

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,
           @x);
}

class LFT {
    has $.q;
    has $.r; 
    has $.s; 
    has $.t;
    # method new($q, $r, $s, $t) { self.bless(*, :$q, :$r, :$s, :$t); }
    method extr($x) { ($.q * $x + $.r) / ($.s * $x + $.t); }
}

sub unit() { LFT.new(q => FatRat.new(1, 1), 
                     r => FatRat.new(0, 1), 
                     s => FatRat.new(0, 1), 
                     t => FatRat.new(1, 1)); }    
sub comp($a, $b) { 
    LFT.new(q => $a.q * $b.q + $a.r * $b.s,
            r => $a.q * $b.r + $a.r * $b.t,
            s => $a.s * $b.q + $a.t * $b.s,
            t => $a.s * $b.r + $a.t * $b.t);
}

sub pi-stream() {
    stream(-> $z { floor($z.extr(3)); },
           -> $z, $n { $n == floor($z.extr(4)); },
           -> $z, $n { comp(LFT.new(q => 10, r => -10*$n, s => 0, t => 1), $z); },
           &comp,
           unit, 
           (1..*).map({ LFT.new(q => $_, r => 4 * $_ + 2, s => 0, t => 2 * $_ + 1) }));
}

my @pi := pi-stream;
say pi;
say @pi[0] ~ '.' ~ @pi[1..100].join('');

Certainly not as elegant as it might ideally be in the long run, but it works today, calculating 101 digits of pi in 5.7 seconds.

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) 

becomes

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’ 

becomes

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); },
           &comp,
           $unit, 
           (1..*).map({ [$_, 4 * $_ + 2, 0, 2 * $_ + 1] }));
}

It’s a very direct translation.

Does it work?

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

Yay!

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 
  where 
    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,
           @x);
}

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],
           @x);
}

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.

Testing

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:

taaatatattgtttcatagtcgtaccacccaacacattgtcctcacg
cataaccggcggctcgtggctttctgtagaccgaatcttcgctgtttgctctg
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
accgccaaccggaaagaatgtcctcccaccacaaatgtacgctcgcatggcggttgtcgagtctatgtcggttgcgctatactacatgataatggacgcc
ttacttacccaaagtagaaataagctcgtctttgagaaccgtggactggtactacctatttttagtcaaactcatgactcgcgcctagcccacatacaat
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
accgccaaccggaaagaatgtcctcccaccacaaatgtacgctcgcatggcggttgtcgagtctatgtcggttgcgctatactacatgataatggacgccttacttacccaaagtagaaataagctcgtctttgagaaccgtggactggtactacctatttttagtcaaactcatgactcgcgcctagcccacatacaat
tttgtccctggccacgacgttctactatatgttaatgaaacgtaaggaattgcgttggccaagaaacgtccttttcacagatacccgtcgtacctgattaccgctgtagggcgctttttccggctggggcgcgcgtgtctgttggccgggccctacgtaggcctataacggaaagatttgtaccaaattctactacgagg
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

Follow

Get every new post delivered to your Inbox.