Beefy Boxes and Bandwidth Generously Provided by pair Networks
Your skill will accomplish
what the force of many cannot
 
PerlMonks  

comment on

( #3333=superdoc: print w/replies, xml ) Need Help??
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 $@;
I tried it with ->bpow(), but that's 2X slower. (I'd offer that ->bpow() should only be used for non-integer powers.)

Seems eval has its uses!

Update 2: I made a typo in above on bdiv -- should have used abs. It's been changed, and the timings updated (roughly doubled).

-QM
--
Quantum Mechanics: The dreams stuff is made of


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

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



  • Posts are HTML formatted. Put <p> </p> tags around your paragraphs. Put <code> </code> tags around your code and data!
  • Titles consisting of a single word are discouraged, and in most cases are disallowed outright.
  • Read Where should I post X? if you're not absolutely sure you're posting in the right place.
  • Please read these before you post! —
  • Posts may use any of the Perl Monks Approved HTML tags:
    a, abbr, b, big, blockquote, br, caption, center, col, colgroup, dd, del, div, dl, dt, em, font, h1, h2, h3, h4, h5, h6, hr, i, ins, li, ol, p, pre, readmore, small, span, spoiler, strike, strong, sub, sup, table, tbody, td, tfoot, th, thead, tr, tt, u, ul, wbr
  • You may need to use entities for some characters, as follows. (Exception: Within code tags, you can put the characters literally.)
            For:     Use:
    & &amp;
    < &lt;
    > &gt;
    [ &#91;
    ] &#93;
  • Link using PerlMonks shortcuts! What shortcuts can I use for linking?
  • See Writeup Formatting Tips and other pages linked from there for more info.
  • Log In?
    Username:
    Password:

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

    How do I use this? | Other CB clients
    Other Users?
    Others contemplating the Monastery: (6)
    As of 2019-10-20 22:18 GMT
    Sections?
    Information?
    Find Nodes?
    Leftovers?
      Voting Booth?
      Notices?