note QM 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 [cpan://Math::Big::Factors]. I have an idea about prefactoring up to some limit, which I'll report on later. <p> 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. <p> Can you comment on the practical limits of your application? <p> For your example, which is essentially: <code> 10389! 4570! 44289! 11800! ------------------------------ 56089! 989! 9400! 43300! 2400! </code> and comparing results:<p> <code> 8.070604647867604097576877675243109941805e-7030 8.070604647867604097576877675243109941729476950993498765160880e-7030 </code> (Note that MBF actually moved the decimal and gave the exponent as "-7069" -- I've normalized it to the other result.)<p> Here's my code: <readmore><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 %{\$factors} ) { 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"; } </code></readmore> <b>Update:</b> I ran some variations, here's the time comparisons: <p> Order = 5, Accuracy = 40, time = 7 min<br/> Order = 6, Accuracy = 60, time = 35 min (!!)<br/> no factoring, Accuracy = 40, time = 28 sec<br/> no factoring, Accuracy = 60, time = 30 sec<br/> no factoring, Accuracy = 40, eval, time = 14.0 sec<br/> no factoring, Accuracy = 60, eval, time = 14.7 sec<br/> <p> The "eval" version changes this code: <code> foreach my \$t ( keys %{\$factors} ) { my \$bft = \$ONE->copy()->bmul(\$t)->bpow(\$factors->{\$t}); \$q->bmul(\$bft); } </code> to this: <code> 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 \$@; </code> I tried it with <code>->bpow()</code>, but that's 2X slower. (I'd offer that <code>->bpow()</code> should only be used for non-integer powers.)<p> Seems <code>eval</code> has its uses! <p> <b>Update 2:</b> I made a typo in above on bdiv -- should have used abs. It's been changed, and the timings updated (roughly doubled). <div class="pmsig"><div class="pmsig-294463"> <p>-QM<br /> --<br /> Quantum Mechanics: The dreams stuff is made of </div></div> 481987 481987