note
hv
<p>Here's my take on it. Highlights: use Eratosthenes' sieve to determine a factor or primality for integers up to the max required; work out when an additional numerator or denominator kicks in (<c>@switch</c>), and pass over the numbers in reverse order, so each need be visited only once; build up two BigInts for numerator and denominator and just switch to BigFloats to divide the two at the end.</p>
<p>The single division at the end takes over 2/3 of the total time, so I'd expect using GMP or Pari for the maths would cut it greatly:<c>
zen% perl t1
807060464786760409757687767524310994172947695099349876516088e-7089
sieve: 0 wallclock secs ( 0.09 usr + 0.02 sys = 0.11 CPU)
build: 4 wallclock secs ( 4.38 usr + 0.00 sys = 4.38 CPU)
divide: 9 wallclock secs ( 8.90 usr + 0.00 sys = 8.90 CPU)
total: 13 wallclock secs (13.37 usr + 0.02 sys = 13.39 CPU)
zen%
</c></p>
<p>Things that made negligible difference: dropping accuracy to 40 (!); initialising <c>my $bii = Math::BigInt->new($i)</c> before the multiply loops; the <c>$r->bsstr</c> call.</p>
<p>I noticed that the numerator consists of 855 prime factors comprising 755 distinct primes; the denominator has 2308 of which 2306 are distinct. No surprise, then, that bpow() isn't a win.</p>
<p>Here's the code:<readmore><c>
#!/usr/bin/perl -w
use strict;
use List::Util qw/ max /;
use Math::BigInt;
use Math::BigFloat;
use Benchmark;
my @num = ( 10389, 45700, 44289, 11800 );
my @den = ( 56089, 989, 9400, 43300, 2400 );
my(@fact, @freq, @switch);
# sieve factors
my $t0 = new Benchmark;
my $max = max(@num, @den);
for (my $i = 2; $i * $i < $max; ++$i) {
++$i while $fact[$i]; # skip to next prime
for (my $j = $i * $i; $j < $max; $j += $i) {
$fact[$j] ||= $i;
}
}
++$switch[$_] for @num;
--$switch[$_] for @den;
my $t1 = new Benchmark;
my $count = 0;
my $num = Math::BigInt->new(1);
my $den = Math::BigInt->new(1);
for my $i (reverse 1 .. $max) {
no warnings qw/ uninitialized /;
$count += $switch[$i];
my $f = $count + $freq[$i];
if ($fact[$i]) {
$freq[$fact[$i]] += $f;
$freq[$i / $fact[$i]] += $f;
} elsif ($f > 0) {
$num *= $i for 1 .. $f;
} elsif ($f < 0) {
$den *= $i for $f .. -1;
}
}
my $t2 = new Benchmark;
Math::BigFloat->accuracy(60);
my $r = Math::BigFloat->new($num) / $den;
print $r->bsstr, "\n";
my $t3 = new Benchmark;
time_delta("sieve", $t1, $t0);
time_delta("build", $t2, $t1);
time_delta("divide", $t3, $t2);
time_delta("\ntotal", $t3, $t0);
sub time_delta {
print +shift, ": ", timestr(timediff(@_)), "\n";
}
</c></readmore></p>
<p>Hugo</p>
481987
483148