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.