Beefy Boxes and Bandwidth Generously Provided by pair Networks
go ahead... be a heretic

comment on

( #3333=superdoc: print w/replies, xml ) Need Help??

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%

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


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

Use:  <p> text here (a paragraph) </p>
and:  <code> code here </code>
to format your post; it's "PerlMonks-approved HTML":

  • 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.
Log In?

What's my password?
Create A New User
Domain Nodelet?
and the web crawler heard nothing...

How do I use this? | Other CB clients
Other Users?
Others wandering the Monastery: (4)
As of 2023-03-20 12:03 GMT
Find Nodes?
    Voting Booth?
    Which type of climate do you prefer to live in?

    Results (59 votes). Check out past polls.