 Perl: the Markov chain saw PerlMonks

### comment on

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

Hugo

Title:
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.

Create A New User
Domain Nodelet?
Chatterbox?
and the web crawler heard nothing...

How do I use this? | Other CB clients
Other Users?
Others surveying the Monastery: (4)
As of 2023-06-09 17:05 GMT
Sections?
Information?
Find Nodes?
Leftovers?
Voting Booth?
How often do you go to conferences?

Results (36 votes). Check out past polls.

Notices?