XP is just a number PerlMonks

### Re^3: Algorithm for cancelling common factors between two lists of multiplicands

by QM (Parson)
 on Aug 12, 2005 at 18:45 UTC ( #483367=note: print w/replies, xml ) Need Help??

Don't factor!

I removed map{ primeFactors \$_ } from your FET3 sub, and it runs 5X faster! (The last 4 digits differ between the 2 versions.)

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

• Comment on Re^3: Algorithm for cancelling common factors between two lists of multiplicands

Replies are listed 'Best First'.
Re^4: Algorithm for cancelling common factors between two lists of multiplicands
by tmoertel (Chaplain) on Aug 12, 2005 at 20:23 UTC
Don't factor!
Indeed, factoring saves time only when almost all of the factors can be canceled. Otherwise, the factors will just be multiplied back together – at great cost – to arrive back at nearly the original terms (magnitude wise), which then must be multiplied with each other to yield the final product.

For example, if we start with 12 and factor it into 2x2x3, do we win or lose if we can cancel only one term? Two BigInt multiplications (e.g., 2x3xN, where N is large) is probably more costly than the original single multiplication (12xN), which makes our factoring a net loss, even though 2x3=6 is smaller than 12. Simply put, when it comes to run time the number of multiplications is what matters most.

We probably win only when we can cancel all or all but one of the factors. Given that our numerators and denominators are often greatly unbalanced, how often will that happen? Our emprical results suggest not enough to be worthwhile. More likely, we will end up with one side of the dividing line plagued with little factors, like aphids on a rosebud, each sucking up an expensive multiplication.

Re^4: Algorithm for cancelling common factors between two lists of multiplicands
by BrowserUk (Pope) on Aug 12, 2005 at 22:32 UTC

With that insight, I removed all the factoring and cancelling code and came up with this version:

```#! perl -slw
use strict;
use Benchmark::Timer; my \$T = new Benchmark::Timer;
use List::Util qw[ sum reduce max ]; our( \$a, \$b );

sub factors{ 2 .. \$_[ 0 ] }

sub normalise {
my( \$s, \$n ) = @{+shift };
while( 1 ) {
if(    \$n > 1.0 ) { \$n /= 10; \$s++; redo; }
elsif( \$n < 0.1 ) { \$n *= 10; \$s--; redo; }
else { last }
}
return [ \$s, \$n ];
}

sub sProduct{
@{
reduce{
\$a->[ 0 ] += \$b->[ 0 ];
\$a->[ 1 ] *= \$b->[ 1 ];
normalise( \$a );
} map{ [ 0, \$_ ] } 1, @_
};
}

sub FET4 {
my @data = @_;
return unless @data == 4;
my @C = ( sum( @data[ 0, 2 ] ), sum( @data[ 1, 3 ] ) );
my @R = ( sum( @data[ 0, 1 ] ), sum( @data[ 2, 3 ] ) );
my \$N =  sum @C;
my( \$dScale, \$d ) = sProduct map{ factors \$_ } grep \$_, @R, @C;
my( \$sScale, \$s ) = sProduct map{ factors \$_ } grep \$_, \$N, @data;
return ( \$d / \$s, \$dScale - \$sScale );
}

die "Bad args @ARGV" unless @ARGV == 4;
print "[@ARGV]";
\$T->start('');
printf "%.17fe%d\n", FET4 @ARGV;
\$T->stop('');
\$T->report;
<STDIN>;
exit;
__END__
P:\test>fet4 989 9400 43300 2400
[989 9400 43300 2400]
0.80706046478686522e-7029
1 trial  of _default (   2.851s total), 2.851s/trial

P:\test>FET4 5 0 1 4
[5 0 1 4]
2.38095238095238090e-2
1 trial  of _default (    564us total), 564us/trial

Which, if you can live with the lower accuracy, for the (5 0 1 4) and (989 9400 43300 2400) datasets, compares favourably with my (rough) timings of tmoertel's compiled Haskell code, the Math::Pari version and blows the M::BF version into dust.

Where it gets really interesting is when you look at bigger numbers. I tried it with (1e5 2e5 2e5 1e5) and it took 81 seconds.

```P:\test>FET4 1e5 2e5 2e5 1e5
[1e5 2e5 2e5 1e5]
1.32499459649722560e-14760
1 trial  of _default (  81.148s total), 81.148s/trial

Math::Pari can't handle numbers this big, and the Haskell version ran for roughly an hour before I abandoned it.

It makes me think that maybe a Math::Big library based around the idea of the scaled math I use in sProduct() might be useful for large number work that doesn't require absolute precision?

Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
Lingua non convalesco, consenesco et abolesco. -- Rule 1 has a caveat! -- Who broke the cabal?
"Science is about questioning the status quo. Questioning authority".
The "good enough" maybe good enough for the now, and perfection maybe unobtainable, but that should not preclude us from striving for perfection, when time, circumstance or desire allow.
I think the canceling code should stay, because it's faster to toss out matching terms than to multiply and then divide them out.

Care to run a benchmark with and without?

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

Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
Lingua non convalesco, consenesco et abolesco. -- Rule 1 has a caveat! -- Who broke the cabal?
"Science is about questioning the status quo. Questioning authority".
The "good enough" maybe good enough for the now, and perfection maybe unobtainable, but that should not preclude us from striving for perfection, when time, circumstance or desire allow.
Re^4: Algorithm for cancelling common factors between two lists of multiplicands
by BrowserUk (Pope) on Aug 12, 2005 at 19:07 UTC

Prior to replacing product() with sProduct(), my earlier non-factoring versions would not have produce a result because the intermediate terms blew the range of a double. That's what originally necessitated the introduction of M::BF and yielded the 4 1.2 hrs runtime.

sProduct() effectively gives me a cheap "BigFloat", and negates the need for the factoring.

Just goes to show there is more than one way to skin a cat :)

Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
Lingua non convalesco, consenesco et abolesco. -- Rule 1 has a caveat! -- Who broke the cabal?
"Science is about questioning the status quo. Questioning authority".
The "good enough" maybe good enough for the now, and perfection maybe unobtainable, but that should not preclude us from striving for perfection, when time, circumstance or desire allow.

Create A New User
Node Status?
node history
Node Type: note [id://483367]
help
Chatterbox?
and the web crawler heard nothing...

How do I use this? | Other CB clients
Other Users?
Others perusing the Monastery: (10)
As of 2019-10-17 15:26 GMT
Sections?
Information?
Find Nodes?
Leftovers?
Voting Booth?
In 2019 the site I miss most is:

Results (42 votes). Check out past polls.

Notices?