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.
(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.
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.
