I know this thread is really old, but it was interesting. I looked at a few different methods.
Getting the factors or divisors. Choices include:
- ntheory. factor returns prime factors in list, factor_exp returns prime factors in factor/exponent list, divisors returns divisors. Fast and low overhead. Supports bigints.
- Math::Pari. factorint returns prime factors (in factor/exponent form), divisors returns divisors. Fast but high overhead because everything is a Math::Pari object. Supports bigints.
- Math::Factor::XS. prime_factors returns prime factors, factors returns divisors. Fast for small numbers, slows down as they get larger. Does not support bigints.
- Math::Big::Factors. factor_wheel returns prime factors. Mind bogglingly slow.
Edit: updated for <= instead of <. Also updated code 3 and reran timing.
Methods for getting the largest divisor less than or equal to sqrt(n):
- generate the list of divisors, pick the element in the middle (odd number of divisors = perfect square) or the left of middle. The divisors are symmetrical around sqrt(n). The interesting cases are: perfect square (there is a true middle element); prime (divisors are 1 and n, we choose 1); 1 (divisors are 1, we ought to choose 1); 0 (no divisors, um...)
- generate the list of divisors, iterate through them selecting the largest less than or equal to square root of n. This is not much slower than the previous as long as we don't have too many divisors, which is usually true.
- generate list of divisors, binary search. Nice idea, but if they're sorted then we already know where the result is.
- limbic's Inline::C version using the factors. Needs winning initialized to 1.
- Ray Johnson's pruner
- Ray Johnson's "other method"
Results. Times for 1 million random integers 0 to 2^47-1 (code 0)
- 21.7s ntheory divisors (code 1)
- 22.9s ntheory fordivisors (code 2)
- 23.4s Inline::C using ntheory factor
- 32.7s dj using ntheory factor (code 3)
- 33.6s pruner using ntheory factor
- 55.9s rjnew using ntheory factor
- 157.2s Pari divisors (code 4)
- 167.9s Inline::C using Pari factorint
- 200.9s limbic using ntheory factor
- 352.0s limbic using Pari factorint
- 95 minutes ntheory divisors in pure Perl
- 987 minutes Math::Factor::XS
- years Math::Big::Factors::factor_wheel
Factoring in pure Perl is definitely slower than C (but maybe someone has faster Perl factoring code). Math::Factor::XS tests all integers from 2 to sqrt(n) without pruning. This is definitely slower at this size than generating the divisors from the prime factors. factors_wheel has two problems: it always uses Perl bigints which are slow, but more importantly it makes an implementation mistake and tests to n instead of sqrt(n. These combined make it completely unreasonable in performance for these 47-bit numbers.
So.. (1) unless you have tiny inputs, you need a fast factoring method, today that would be ntheory or Pari. (2) Pari's overhead gets in the way here, but it still yields reasonably fast solutions. (3) the simple divisors call (do everything in C) is fastest, but we can do a good job given just the factors. My method of creating all the divisors from the factor list was slow at first, but improving the code made it much better.
(code 0)
srand(1); # So everything is using the same numbers
my $s = 0;
for (1..1000000) {
$s += nearest_sqrt(int(rand(2**47)));
}
print "$s\n";
(code 1)
use ntheory qw/divisors/;
sub closest {
my @d = divisors(shift);
$d[ $#d >> 1 ];
}
(code 2)
use ntheory qw/fordivisors/;
sub closest {
my($n,$sqrtn,$min) = ($_[0], int(sqrt($_[0])));
fordivisors { $min = $_ if $_ <= $sqrtn; } $n;
$min;
}
(code 3)
use ntheory qw/factor/;
sub closest {
my($n) = @_;
my(@factors, @d, @t);
# 1. Get list of factors
@factors = factor($n);
# 2. Turn this into an unsorted list of divisors
@d = (1);
while (my $p = shift @factors) {
my $e = 1;
while (@factors && $p == $factors[0]) { $e++; shift(@factors); }
push @d, @t = map { $_ * $p } @d; # multiply throug
+h once
push @d, @t = map { $_ * $p } @t for 2 .. $e; # repeat
}
# 3. Sort
@d = sort {$a<=>$b} @d;
# 4. Return the largest divisor <= sqrt n.
$d[ $#d >> 1 ];
}
(code 4)
use Math::Pari qw/divisors/;
sub closest {
my $d = divisors(shift);
$d->[ $#$d >> 1 ];
}
(code 5)
use Math::Factor::XS qw/factors/;
sub closest {
my @d = (1, factors($_[0]), $_[0]);
$d[ $#d >> 1 ];
}