Beefy Boxes and Bandwidth Generously Provided by pair Networks
The stupid question is the question not asked
 
PerlMonks  

Re^2: OT: Finding Factor Closest To Square Root

by Roy Johnson (Monsignor)
on Feb 28, 2005 at 18:13 UTC ( [id://435134]=note: print w/replies, xml ) Need Help??


in reply to Re: OT: Finding Factor Closest To Square Root
in thread OT: Finding Factor Closest To Square Root

I can't get Math::Pari to compile on my machine, so I did an inferior benchmarking test (code below) between a pruning combinations method I came up with and your method (which I modified -- I think fairly -- to work on the list of factors instead of generating them). My method ran over three times as fast, but as I say, I didn't benchmark over extensive input values. If it's convenient for you to do so, I'd like to see what results you get when you plug it into your benchmarker.

Note that I have an optimization that can be commented out. In some cases (including the data used here) it gave me a speedup. In others, a slowdown. It would be interesting to see whether it was, on average, an improvement or not.

I just added my other method to the test list, and it absolutely blows them out of the water. I had no idea it could be that much faster.

use strict; use warnings; use List::Util 'reduce'; # This is just in here for reference because it is easier to # read than the sub below. It is a non-recursive combinations # generator sub all_combos { my @accum = (); for my $elem (@_) { push @accum, ([$elem], map { [@$_, $elem] } @accum); } @accum; } # Generates all distinct products of the given factor list, # but stops combining when a product or element is greater than # the given prune value, and will toss any factors exceeding the # limit if they aren't the smallest such factors seen so far sub all_combos_prune { my $limit = shift; my @accum = (); my %seen; my $SEF = $limit * $limit; # Smallest exceeding factor for my $elem (@_) { push @accum, (($seen{$elem}++ ? () : $elem), map { ($elem < $limit and $_ < $limit) ? do { my $n = $elem * $_; $SEF = $n if $n > $limit and $n < $SEF; ($SEF < $n or $seen{$n}++) ? () : $n; # Either the two lines above should be commented, or the line below sh +ould # $seen{$n}++ ? () : $n; } : () } @accum); } @accum; } sub limbic { my $sqrt = shift; my $N = shift; return 1 if @_ == 1; my (@subset, %seen, $winner, $offset, $f, $diff); my $end = $#_; my ($pos, $mode) = (-1, 1); while ( 1 ) { if ( $mode == 1 ) { push @subset, ++$pos; ++$mode if $subset[ -1 ] == $end; } elsif ( $mode == 2 ) { splice(@subset, $#subset - 1, 1); ++$mode; } else { return $winner if $subset[ 0 ] == $end; $pos = $subset[ -2 ] + 1; splice(@subset, $#subset - 1, 2, $pos); $mode = 1; } if ( $seen{ "@_[ @subset ]" }++ ) { if ( $mode == 1 ) { push @subset, $pos + 1 .. $end; ++$mode; } next; } $f = 1; $f *= $_ for @_[ @subset ]; next if $f > $sqrt; $diff = ($N / $f) - $f; if ( ! defined $offset || $diff < $offset ) { ($winner, $offset) = ($f, $diff); } } } sub closest { # Args are target and factor-list my ($target, $num, @factors) = @_; # Take the biggest factor my $top_factor = pop @factors; # Find multiple of that factor closest to (and above) target my $guess = int($target) - $target % $top_factor + $top_factor; # Oscillate around the target, looking at multiples of top_factor # until you get one that divides the product for (my $i = $top_factor; $num % $guess; $i += $top_factor) { $guess += ( $target <=> $guess ) * $i; } # Check the complementary factor my $complement = $num / $guess; # Look for a multiple of $top_factor between the last guess on the # other side of sqrt and $complement my $direction = ($target <=> $guess); my $new_guess = $guess + $direction * $top_factor; while (($complement <=> $new_guess) == $direction and $num % $new_ +guess) { $new_guess += $top_factor * $direction; } if ($new_guess and $num % $new_guess == 0 and (($complement <=> $new_guess) == $direction)) { $guess = $new_guess; $complement = $num/ $guess; } abs($target - $complement) < abs($target - $guess) ? $complement : $guess ; } my @facts = (2,3,3,3,3,11,13,17,19,29,23); my $prod = reduce {$a * $b} @facts; my $root = sqrt $prod; print "P=$prod, R=$root\n"; my @combos = all_combos_prune($root, @facts); my $closest = reduce {abs($root - $a) < abs($root - $b) ? $a : $b} @co +mbos; print "Closest is $closest\n"; print "Limbic says: ", limbic($root, $prod, @facts), "\n"; print "My other way says: ", closest($root, $prod, @facts), "\n"; use Benchmark 'cmpthese'; cmpthese( -2, { pruner => sub { reduce {abs($root - $a) < abs($root - $ +b) ? $a : $b} all_combos_prune($root, @facts) }, limbic => sub { limbic($root, $prod, @facts) }, my_old => sub { closest($root, $prod, @facts) } });
Closest is 70499 Limbic says: 70499 My other way says: 70499 Rate limbic pruner my_old limbic 8.02/s -- -82% -100% pruner 45.1/s 462% -- -98% my_old 2632/s 32727% 5737% --
For P=4789309440 (factors (2)x12,3,5,77951), the differences are astronomical:
P=4789309440, R=69204.8368251815 Closest is 61440 Limbic says: 61440 My other way says: 61440 Rate limbic pruner my_old limbic 0.749/s -- -100% -100% pruner 181/s 24119% -- -97% my_old 5545/s 740213% 2957% --
and if you try P=2176782336 (which is an exact square with a factor list of ((2)x12, (3)x12))), you will give up waiting for your method.

Caution: Contents may have been coded under pressure.

Replies are listed 'Best First'.
Re^3: OT: Finding Factor Closest To Square Root
by Limbic~Region (Chancellor) on Feb 28, 2005 at 19:31 UTC
    Roy Johnson,
    I am sure once you get all the bugs worked out one of your solutions will be the winner (minus my Inline::C solution). Let me know when you think it is bug free and I will update. As it stands now, it breaks on 13316184461 (45863,290347). It says that 290347 is the answer, but that is larger than the sqrt().

    Cheers - L~R

      I had left in a call to shift that I shouldn't. It's fixed, now.

      Caution: Contents may have been coded under pressure.
        Roy Johnson,
        It's fixed, now.

        Mostly. It gives the number instead of 1 if the number is prime. That really isn't a big deal.

        $ time ./rjnew.pl real 1m8.446s user 1m1.338s sys 0m0.440s
        So you have a couple of the fastest implementations. The Inline::C is only barely winning though that might be static. The interesting thing is that many solutions are faster than using Math::Pari to get a list of factors and doing a binary search.

        Cheers - L~R

Log In?
Username:
Password:

What's my password?
Create A New User
Domain Nodelet?
Node Status?
node history
Node Type: note [id://435134]
help
Chatterbox?
and the web crawler heard nothing...

How do I use this?Last hourOther CB clients
Other Users?
Others rifling through the Monastery: (7)
As of 2024-03-28 14:02 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?

    No recent polls found