http://qs1969.pair.com?node_id=483148

I have toyed with this more, and have a pure Perl solution that runs your example in about 7 minutes on my PC, and all but 5 seconds is the factoring by Math::Big::Factors. I have an idea about prefactoring up to some limit, which I'll report on later.

You might note that 2**32 is not a practical limit, unless most of the terms cancel out. For instance, if you had to compute (2**32)!/(2**31)!, that would be 2**31 terms in the numerator.

Can you comment on the practical limits of your application?

For your example, which is essentially:

```10389! 4570! 44289! 11800!
------------------------------
56089! 989! 9400! 43300! 2400!
and comparing results:

```8.070604647867604097576877675243109941805e-7030
8.070604647867604097576877675243109941729476950993498765160880e-7030
(Note that MBF actually moved the decimal and gave the exponent as "-7069" -- I've normalized it to the other result.)

Here's my code:

```#!/your/perl/here

use strict;
use warnings;
use Benchmark;
use Math::Big::Factors;

my \$ORDER = 5;

# preserve this many significant digits
Math::BigFloat->accuracy(40);
my \$ONE = Math::BigFloat->new(1);

my \$t0 = new Benchmark;

my @num = ( 10389, 45700, 44289, 11800 );
my @den = ( 56089, 989, 9400, 43300, 2400 );

my \$t1 = new Benchmark;
time_delta("startup", \$t1, \$t0);

# expand terms
my (@num_fac, @den_fac);
foreach my \$n ( @num )
{
push @num_fac, 2..\$n;
}
foreach my \$n ( @den )
{
push @den_fac, 2..\$n;
}
my \$t2 = new Benchmark;
time_delta("expand", \$t2, \$t1);
my \$quotient = cancel_quotient( \@num_fac, \@den_fac );
my \$t3 = new Benchmark;
time_delta("cancel_quotient", \$t3, \$t2);

warn "factoring ", scalar keys %{\$quotient}, " terms\n";
my \$factors = prime_factor_hash( \$quotient );
my \$t4 = new Benchmark;
time_delta("prime_factor_hash", \$t4, \$t3);

# evaluate terms
my \$q = \$ONE->copy();
#print "Factors\n";
foreach my \$t ( sort { \$factors->{\$b} <=> \$factors->{\$a} } keys %{\$fac
+tors} )
{
my \$bft = \$ONE->copy()->bmul(\$t)->bpow(\$factors->{\$t});
\$q->bmul(\$bft);
#    print "\$t**\$factors->{\$t}:  \$q\n";
}

print \$q->bsstr(), "\n";
my \$t5 = new Benchmark;
time_delta("multiply", \$t5, \$t4 );
time_delta("All", \$t5, \$t0 );
exit;

########################################
# cancel like terms between 2 arrays
# similar to unique array element

sub cancel_quotient
{
my \$num = shift;
my \$den = shift;

my %quotient;

# determine the residual terms after cancelling
# positive terms are in numerator
# negative terms are in denominator
# zero terms are filtered
foreach my \$n ( @{\$num} )
{
\$quotient{\$n}++;
}
foreach my \$n ( @{\$den} )
{
\$quotient{\$n}--;
}
foreach my \$n ( keys %quotient )
{
delete \$quotient{\$n} unless \$quotient{\$n};
}

return \%quotient;
}

#######################################
# factor keys of hash, and duplicate "value" times

sub prime_factor_hash
{
my \$quotient = shift;
my %factors;

foreach my \$n ( keys %{\$quotient} )
{
my @factors = Math::Big::Factors::factors_wheel(\$n,\$ORDER);
foreach my \$f ( @factors )
{
\$factors{\$f} += \$quotient->{\$n};
}
}

return \%factors;
}

###########################################
# nice benchmark timestamp printer

sub time_delta
{
print +shift, ": ", timestr(timediff(@_)), "\n";
}
Update: I ran some variations, here's the time comparisons:

Order = 5, Accuracy = 40, time = 7 min
Order = 6, Accuracy = 60, time = 35 min (!!)
no factoring, Accuracy = 40, time = 28 sec
no factoring, Accuracy = 60, time = 30 sec
no factoring, Accuracy = 40, eval, time = 14.0 sec
no factoring, Accuracy = 60, eval, time = 14.7 sec

The "eval" version changes this code:

```foreach my \$t ( keys %{\$factors} )
{
my \$bft = \$ONE->copy()->bmul(\$t)->bpow(\$factors->{\$t});
\$q->bmul(\$bft);
}
to this:
```my \$s = '\$q';
foreach my \$t ( keys %{\$factors} )
{
next if ( \$factors->{\$t} == 0 );
if ( \$factors->{\$t} > 0 )
{
\$s .= "->bmul(\$t)" x \$factors->{\$t};
}
else
{
\$s .= "->bdiv(\$t)" x abs(\$factors->{\$t});
}
}
eval \$s;
die "Error on eval: \$@" if \$@;