comment on

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 (@switch), 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.

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:

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%
[download]

Things that made negligible difference: dropping accuracy to 40 (!); initialising my $bii = Math::BigInt->new($i) before the multiply loops; the $r->bsstr call.

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.

Here's the code:

#!/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";
}
[download]

Hugo

In reply to Re^2: Algorithm for cancelling common factors between two lists of multiplicands by hv
in thread Algorithm for cancelling common factors between two lists of multiplicands by BrowserUk

Are you posting in the right place? Check out Where do I post X? to know for sure.
Posts may use any of the Perl Monks Approved HTML tags. Currently these include the following:
<code> <a> <b> <big> <blockquote> <br /> <dd> <dl> <dt> <em> <font> <h1> <h2> <h3> <h4> <h5> <h6> <hr /> <i> <li> <nbsp> <ol> <p> <small> <strike> <strong> <sub> <sup> <table> <td> <th> <tr> <tt> <u> <ul>
Snippets of code should be wrapped in <code> tags not <pre> tags. In fact, <pre> tags should generally be avoided. If they must be used, extreme care should be taken to ensure that their contents do not have long lines (<70 chars), in order to prevent horizontal scrolling (and possible janitor intervention).
Want more info? How to link or How to display code and escape characters are good places to start.


Perl: the Markov chain saw
	PerlMonks