You have the right of it. More precisely, take the logarithim. Then the problem becomes finding a sum of log prime factors that is closest to log(N)/2. This problem is reducible to the subset sum problem, which is NP-hard.

Although exact solutions are exponential, heuristics can sometimes generate good answers with minimal effort. One example of a heuristic that might work well is the greedy heuristic. To fill a box of size log(N)/2, put the largest factor into the box. Then put the next largest factor that still fits in the box and insert it. Keep doing this until no more factors fit.

The above will get you the largest factor <= sqrt(N). To get the smallest factor >= sqrt(N), put all the factors in a box, and take them out again largest to smallest, until no more can be taken out. Finally, comapre the two answers to find the closest.

Update: fixed a typo.

-Mark

To fill a box of size log(N)/2, put the largest factor into the box. Then put the next largest factor that still fits in the box and insert it. Keep doing this until no more factors fit.

The above will get you the largest factor <= sqrt(N). To get the smallest factor >= sqrt(N), put all the factors in a box, and take them out again largest to smallest, until no more can be taken out. Finally, comapre the two answers to find the closest.

I don't think that does what you think. For example, 2**3 * 19**3 = 54872, and the square root is 234.277....

Your method gives 19 as the largest less than the square root, and 2**3 * 19**2 = 2888 as the smallest greater than the square root. Notice that 19*2888 = 54872.

A better result is 2**3 * 19 = 152 (because 19**2 = 361).

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

Sorry, I was sloppy before...

Math::Big::Factor::factors_wheel() doesn't give all factors, but the prime factorization.

For instance,

factors_wheel(1000) == (2, 2, 2, 5, 5, 5);
[download]

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

Since you've got them in ascending order, take every other prime factor and multiply. Compare that to the product of the ones you didn't choose. Whichever's closer, there ya go. That's just off the top of my head.

Let's try another way. Pretty clearly, you want roughly half the factors. You can work from the outside in, or the inside out, or do something like this: If you have N factors, take the (N*2)th root of the number. Find the factor closest to that. Divide your number by the square of the chosen factor, decrement N, and repeat until you get further from the square root instead of closer to it.

---

Caution: Contents may have been coded under pressure.

If you get a list of prime factors in ascending order, then I would take every other member of the list, multiply them together and start with that.

In your example, that would give you 20, which isn't too far from the correct result of 31.

Hmm .. Actually, an even better answer would be

• for a list with an odd number of elements, take all of the odd position numbers and multiply them together;
• for a list with an even number of elements, take the odd numbered elements from 1 to (n/2)-1, and the even numbered elements from n down to (n/2)+2, and multiply them together with the harmonic mean of the (n-1)/2 and n/2 elements.

For your example, that would be 2 * sqrt(2*5) * 5, which turns out to be exactly the correct answer, 31.622..

To try out the odd number, we'll try out 2000, which gives us a list of (2, 2, 2, 2, 5, 5, 5) and a result of 2 * 2 * 5 * 5 or 100. Hmm, a little high.

Well, that's a fascinating question, and good luck with that.

Alex / talexb / Toronto

"Groklaw is the open-source mentality applied to legal research" ~ Linus Torvalds

One way to do it, similar to talexb, would be (tested!)

use Math::Big::Factors;

my $val = 2000;

my @x = Math::Big::Factors::factors_wheel( $val );
print "@x\n";

my %factors;
$factors{$_}++ for @x;

use Data::Dumper;
print Dumper \%factors;

my $sqrt = 1;
while (my ($k,$v) = each %factors) {
    while ($v > 1) {
        $sqrt *= $k;
        $v -= 2;
    }

    if ($v) {
         $sqrt *= sqrt( $k );
    }
}

print "$sqrt => ", sqrt( $val ), $/;
[download]

This is going to be slower that just taking the sqrt of the value. Just to let you know ...

Being right, does not endow the right to be rude; politeness costs nothing.
Being unknowing, is not the same as being stupid.
Expressing a contrary opinion, whether to the individual or the group, is more often a sign of deeper thought than of cantankerous belligerence.
Do not mistake your goals as the only goals; your opinion as the only opinion; your confidence as correctness. Saying you know better is not the same as explaining you know better.

Not sure if this is correct either, but I checked a few cases and it appears to work. Is there a good method of verification?

It starts to get slow with big numbers, but it's much quicker than factors_wheel() finding the prime factors, for all the cases I tried.

There appears to be a bug in the overloading in Math::Big, hence the de-Bignumming map.

Update: Version 3. Anyone care to break this one for me?

#! perl -slw
use strict;
use overload;
use Math::Big::Factors qw[ factors_wheel ];
use List::Util qw[ reduce min ];
use Algorithm::Loops qw[ NestedLoops ];

my $NUM = $ARGV[ 0 ] || 996;

my $root = sqrt $NUM;
my @pfs = reverse map{ "$_" } factors_wheel( $NUM );

print "$NUM primefactors ( @pfs )";

my $n = $pfs[ 0 ];

while( @pfs > 1 ) {
#    print "@pfs";

    my $near = reduce{ 
        ( $a * $b ) <= $root ? $a*$b : $a
    }  @pfs;
#    print "$near : ", $NUM%$near;

    $n = $near 
        if ( $NUM % $near ) == 0 
        and abs( $root - $near ) < abs( $root - $n );

    my $discard = shift @pfs;
    next if $pfs[ 0 ] == $discard;
    unshift @pfs, $pfs[0] while $root > reduce{ $a * $b } @pfs;
}

print "$NUM ($root} : $n";

__END__
P:\test>432558 1995
1995 primefactors ( 19 7 5 3 )
1995 (44.6654228682546} : 35

P:\test>432558 1994
1994 primefactors ( 997 2 )
1994 (44.6542271235322} : 997

P:\test>432558 1993
1993 primefactors ( 1993 )
1993 (44.6430285710994} : 1e+099

P:\test>432558 1993
1993 primefactors ( 1993 )
1993 (44.6430285710994} : 1993

P:\test>432558 1994
1994 primefactors ( 997 2 )
1994 (44.6542271235322} : 997

P:\test>432558 1995
1995 primefactors ( 19 7 5 3 )
1995 (44.6654228682546} : 35

P:\test>432558 9999999
9999999 primefactors ( 4649 239 3 3 )
9999999 (3162.27750205449} : 2151

P:\test>432558 99999999
99999999 primefactors ( 137 101 73 11 3 3 )
99999999 (9999.99995} : 7373

P:\test>432558 999999999
999999999 primefactors ( 333667 37 3 3 3 3 )
999999999 (31622.7765858724} : 2997

P:\test>432558 1000
1000 primefactors ( 5 5 5 2 2 2 )
1000 (31.6227766016838} : 25

P:\test>432558 2000
2000 primefactors ( 5 5 5 2 2 2 2 )
2000 (44.7213595499958} : 40

P:\test>432558 20000
20000 primefactors ( 5 5 5 5 2 2 2 2 2 )
20000 (141.42135623731} : 125

P:\test>432558 200000
200000 primefactors ( 5 5 5 5 5 2 2 2 2 2 2 )
200000 (447.213595499958} : 400
[download]

Maybe I don't understand the OP's requirements, but 36 isn't a factor of 996 (996/36 = 27.666...). I'm guessing he wants 12 (3*2*2) for that case. And 32 isn't a factor of 1000, but 25 is. etc. Am I missing something?

---

-- All code is 100% tested and functional unless otherwise noted.

Am I missing something?

Nope. Nothing at all.

Your just witnessing me loosing sight of the OP's requirements whilst getting caught up with the mechnics of the code. :(

---

Examine what is said, not who speaks.

Silence betokens consent.

Love the truth but pardon error.

I updated 432636. Does the new version look any better?

That's the problem with hueristics--how do you validate them?

---

Examine what is said, not who speaks.

Silence betokens consent.

Love the truth but pardon error.

Update: Ignore this, it doesn't work.

You'll have tp pretend that I have M::B::F::f_w().. CPAN's being boring.

#! perl -slw
use strict;
use List::Util qw[ reduce ];

my $root = sqrt 1000;
my @pfs = ( 2, 2, 2, 5, 5, 5 );

my $near = reduce{ 
    $a * $b < $root ? $a*$b : $a
} reverse @pfs;

print +( $root - $near ) < ( $near * $pfs[0] - $root )
    ? $near
    : $near * $pfs[0];
[download]

I rearranged it a bit so you can just plug in your factor list and it will figure out what your number is. Note that in this case, the closest value is 5*7, but it returns 3*19.

use strict; 
use warnings;
use List::Util qw[ reduce ];

my @pfs = ( 3,5,7,19 );
my $num = reduce { $a * $b } @pfs;
my $root = sqrt $num;
print "N=$num, R=$root\n";

#the rest is the same as you had it
[download]

Here's a fix:

use strict; 
use warnings;
use List::Util qw[ reduce ];

my @pfs = ( 3,5,7,19 );
my $num = reduce { $a * $b } @pfs;
my $root = sqrt $num;
print "N=$num, R=$root\n";

my $near = reduce{ 
    abs($root - $a * $b) < abs($root - $a) ? $a*$b : $a
} reverse @pfs;

print abs($root - $near) > abs($root - ($num/$near))
    ? $num/$near
    : $near
;
[download]

---

Caution: Contents may have been coded under pressure.

One testcase was never going to hack it. I should have thought of your method, but I was too busy trying to figure out why CPAN would install M::B::F.

---

Examine what is said, not who speaks.

Silence betokens consent.

Love the truth but pardon error.

Maybe one of the mathematician's will show me the error of my ways, but I think that the prime factors thing is a red-herring.

The following brute force approach works quite quickly for numbers up to about 12 or 13 digits, usually much more quickly that it takes M::B::F to find the prime factors.

If you're looking to do this with very large numbers (as suggested by your use of Math::Big), then your probably in for a long wait either way. If the number in question is itself prime, then the nearest factor to the square root will be 1. That means you have to add 1 to list that factor_wheel() produces before you start doing your combinatorics, and it's presence just slows the search to what effectively amounts to a linear search anyway.

#! perl -slw
use strict;

my $NUM = $ARGV[ 0 ] || die 'No arg';

my $root = sqrt( $NUM );
my( $lo, $hi ) = ( int( $root ), int( $root + 1 ) );;

$lo-- while $NUM % $lo;
$hi++ while $NUM % $hi;

my $near = ( $lo, $hi )[ abs( $root - $lo ) > abs( $root - $hi ) ];

print "$NUM ($root) : $near";
[download]

Do you really have to do a search for numbers greater than the sqrt(N)? (See Re: OT: Finding Factor Closest To Square Root for my approach). Certainly you don't have to go beyond 2*sqrt(N) on the ++ side. In fact I'm currently thinking the closest factor has to be less than sqrt(N), eliminating any checking for factors greater than sqrt(N). I'd be interested to see if someone could construct a counter-example (a number where the closest factor to sqrt(N) is greater than sqrt(N)).

---

-- All code is 100% tested and functional unless otherwise noted.

Do you really have to do a search for numbers greater than the sqrt(N)?

You right. And now you pointed it out, it seems fairly obvious:

3600: sqrt=60; 59* 61= 3599

10000: sqrt=100; 99*101= 9999

12: sqrt=3.464; 3* 4= 12

Which, in very non-technical terms because my memory doesn't go back to my formal math days, says to me that:

with lo = int( sqrt( N ) ) & hi = int( sqrt( n ) )+1; lo * hi is always less than or equal to N.

If you reduce lo, leaving hi the same, you get further away from N. And if you increase hi, lo will be nearer.

Does that form an explaination? Hmm. Maybe one of our resident math whiz can phrase that correctly?

---

Examine what is said, not who speaks.

Silence betokens consent.

Love the truth but pardon error.

Certainly you don't have to go beyond 2*sqrt(N) on the ++ side

The high side search would always have terminated there or earlier anyway as it stopped as soon as the difference between it and the root was greater than that between lo and root. Lo can't go below 1, so hi would never go above root*2-1. But skipping the hi search is much better.

Another optimisation possible if N is odd, is to step back by 2 each time.

---

Examine what is said, not who speaks.

Silence betokens consent.

Love the truth but pardon error.

My feeling is that you are correct that the prime factorization is a red-herring for the reasons I mentioned in another post in this thread (solving the problem using that info is probably NP hard.) However, finding any nontrivial factor of a large number may be just as hard. Clearly, if one could do that in a reasonable time, then one could also find the prime factorization in a reasonable time. This would upset many people so much that I believe it is forbidden to publish any such result...even partial ones!! (Such results would defeat many common encryption schemes.)
chas
(***This was supposed to be in reply to BrowserUk's 432772 post,but it didn't seem to appear in the right position...likely something I did wrong...)

Can you tell us, why you need to find that one? May be we can come up with non-exhaustive alternatives.

OK, but I was so enjoying the direction the discussion was going, and I didn't want the thread to lose the focus on the problem I stated. But here goes...

At various times I've been bitten by the Fibonacci bug, like many others. I've tried my hand at programming F(n), as well as some other interesting things like a prime tester .

Poking around on Mathworld, I read the Fibonacci page. One of the formulas is

F(k*n) = L(k)*F(k*(n-1)) - ((-1)**k)*F(k*(n-2))
[download]

where

L(k) = F(n-1) + F(n+1)
[download]

and substituting

F(k*n) = (F(n-1)+F(n+1))*F(k*(n-1)) 
          - ((-1)**k)*F(k*(n-2))
[download]

For computing F(k*n), this is most efficient when k==n, giving

F(n**2) = (F(n-1)+F(n+1))*F(n*(n-1)) 
          - ((-1)**n)*F(n*(n-2))
[download]

If N is prime, use

F(2*k+1) = F(k+1)**2 + F(k)**2
[download]

As an exercise in programming and improving my grasp of math and computational complexity, I was going to benchmark several methods of computing F(N). Included in this is the standard recursive solution (which is too slow, but we can estimate it's behavior for large N), the closed form using phi, and various identities.

The identity for F(k*n) seemed to be the most promising in this area, though computing sqrt(N), and/or factoring N, would seem to negate any benefits. Perhaps a quick check to see if N has any small factors, otherwise use the F(2*k+1) identity for odd N, or one of the F(2*k) formulas for even N. (There's a nice identity for N % 3 == 0 too.)

Update: fixed missing parens in formulae.

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

Programming F(n) for benchmarking purposes is interesting. However, the fastest way of all is probably:


F(n)= integer nearest (1/sqrt 5)[(1+sqrt 5)/2]^n
[download]


This isn't perl, but it would only take a line. (If n is extremely large, the exponentiation might become inaccurate, though...)
chas

#!/usr/bin/perl -w
use strict;

print sqrt_factor(shift), "\n";

sub sqrt_factor 
{
    my $n = shift;
    my $root = int(sqrt($n));

    for (my $i=$root; $i>1; $i--)
    {
        if( ($n % $i) == 0 )
        {   
            my $factor2 = int($n / $i);
            return (($root - $i)<($factor2 - $root)) ? $i : $factor2
        }
    }
    
    return 1;
}
[download]

Yes, but this still takes a long time for large prime N, as it has to count down from sqrt(N) to 1.

It would go twice as fast if we checked N%2==0. It would go six times as fast if we checked for N%3==0 also. Perhaps using GCD(N,sqrt(N)#) (where N# is the primorial of N) would improve the situation? The problem then becomes one of generating all primes upto sqrt(N), where we could cheat and just have a big list precomputed.

...Better yet, would a quick test of N for primality solve the problem of prime N?

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

my @res = (2, 2, 2, 5, 5, 5);
my $goal = 1000;

print closest($goal, \@res), "\n";

sub closest{
    my $exact = shift;
    my $array = shift;

    my $best = 0;

    for my $n ( get_numbers($array) ){
        $best = $n if abs( $exact - $n**2) < abs( $exact - $best**2);
    }
    $best;
}

sub get_numbers{
    my $array = shift;

    my $length = @$array;

    if ($length == 1){
        return ( $array->[0] );
    }
    else{
        my @ret_arr;
        my $n = shift @$array;
        push @ret_arr, $n;
        my @tmp = get_numbers($array);
        push @ret_arr, @tmp;
        push @ret_arr, $n*$_ for @tmp;
        prune(@ret_arr);
    }
}

sub prune{
    my %h;

    undef $h{$_} for @_;
    keys %h;
}
[download]

I am almost sure this does what you want. I am also convinced it is rather inefficient.

Yes on both counts. But I like the hand-crafted code approach :)

I rather hoped that someone would reveal a method using the prime factorization, or an efficient method for finding a factor reasonably close to the square root.

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

use strict; 
use warnings;
use List::Util qw[ reduce ];

my $num;  # Used in sub, but could be passed if you wanted

sub closest {
    # Args are target and factor-list
    my ($target, @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;
    abs($target - $complement) < abs($target - $guess)
      ? $complement
      : $guess
    ;
}
my @pfs = (2,2,2,3,3,3,5,5,5,17,19,19,19,19);
# Compute product of factors
our ($a, $b);
$num = reduce { $a * $b } @pfs;
my $root = sqrt $num;
print "N=$num, R=$root\n";

print closest($root, 1, @pfs), "\n";
[download]

Roy Johnson,
Here it is, fast and accurate.

Well it is certainly fast, but unfortunately it is not always accurate. Try (3, 5, 16381, 77951) - it should produce 81905 but it produces 77951. I didn't bother debugging as I was gathering statistics.

Cheers - L~R

You're right. I had an erroneous understanding: that the closest factor that included top_factor would be paired with the closest factor that didn't. Instead, it's paired with the closest factor on the other side of root that doesn't include top_factor. So you have to check both sides of root. Here's a fixed version:

use strict; 
use warnings;
use List::Util qw[ reduce ];

my $num;  # Used in sub, but could be passed if you wanted

sub closest {
    # Args are target and factor-list
    my ($target, @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
    my $i;
    for ($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 * $i;
    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 @pfs = (3, 5, 16381, 77951);
# Compute product of factors
our ($a, $b);
$num = reduce { $a * $b } @pfs;
my $root = sqrt $num;
print "N=$num, R=$root\n";

print closest($root, 1, @pfs), "\n";
[download]

---

Caution: Contents may have been coded under pressure.

That's sweet. I especially like the use of <=> in

$guess += ( $target <=> $guess ) * $i;
[download]

to choose which direction to move.

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

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";
[download]

(code 1)

use ntheory qw/divisors/;
sub closest {
  my @d = divisors(shift);
  $d[ $#d >> 1 ];
}
[download]

(code 2)

use ntheory qw/fordivisors/;
sub closest {
  my($n,$sqrtn,$min) = ($_[0], int(sqrt($_[0])));
  fordivisors { $min = $_ if $_ <= $sqrtn; } $n;
  $min;
}
[download]

(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 ];
}
[download]

(code 4)

use Math::Pari qw/divisors/;
sub closest {
  my $d = divisors(shift);
  $d->[ $#$d >> 1 ];
}
[download]

(code 5)

use Math::Factor::XS qw/factors/;
sub closest {
  my @d = (1, factors($_[0]), $_[0]);
  $d[ $#d >> 1 ];
}
[download]

First, thanks for the good work. I go back and reread this sometimes, and somehow manage, mostly, not to get bitten, being satisfied rereading all the excellent ideas. It was loads of fun at the time, and still tickles my fancy even now.

In your code 1, I find this a little odd:

$d[($#d-1)>>1];
[download]

If there are 3 divisors, then $#d == 2, and $#d -1 == 1, and 1>>1 == 0. So instead of the middle divisor, we get 1. If there is 1 divisor, we get -1 (I think), which perhaps doesn't matter, giving the last element by Perl magic.

But overall, I would agree a priori that something well written that returns a list of divisors will be faster than something returning a list of prime factors, which then have to be recombined to get divisors. Unless there is a very large list of divisors, and a very short list of prime factors, but that is an edge case, I think.

Thanks for refreshing one of my favorite threads.

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

Re code 1, I was solving for the largest divisor less than the square root so that was on purpose. Rereading the thread, we want instead <=. I've edited my answer to reflect this. Thanks for pointing that out.

While the results for the random numbers were identical, clearly they didn't include any perfect squares. I checked results from 2 to 100k and they now match the previous solutions (inlinec, pruner, rjnew, limbic). I started at 2 because some of those methods barf when given the input 1.

This wouldn't change the performance numbers. I did put in a different version of code 3 which may be faster, but I doubt it would change the timings by a lot.

---

Re divisors vs. factors, I think that is only true for very small inputs. That is, my thinking is that for factors we have lots of ways to speed things up -- e.g. trial division by only primes to square root of n, or better methods (e.g. Rho, P-1, SQUFOF, etc.). For divisors we need to do trial division by all numbers to square root of n and each found divisor isn't pruning the search space. The obvious optimization methods are limited forms of the "find factors, multiply out to get divisors" method.

To put this to the test, I used Math::Factor::XS. It has a function factors which loops from 2 to sqrt(n) adding the divisors to a list. It also has prime_factors which gives the prime factors using a mod-6 wheel trial division. I made two programs:

1. call factors and return the middle element.
2. call prime_factors, use my code 3 routine to calculate the divisors, return the middle element.

Is the first method faster? It's getting an array from an XS routine then returning one element from it. It should be fast. The second gets a smaller array from the XS routine then putzes around in Perl code making the divisor array. Yet the second method is faster once the inputs exceed about 400,000.

On the other hand, if the function that returns the divisors generates them via the factors (like the ntheory module does in both C and Perl, and Pari does in C), then we're probably better off just calling it. As the fastest solution on my list shows, as does the Pari divisors vs. various solutions using Pari factors. I could see some exceptions for special cases like primorials where some trimming solution could save some time (e.g. primorial(61) with 18 factors and 262,144 divisors)

Given a number N, find the nearest factor of N less than or equal to the square root of N.

For that question, tall_man's solution is quite concise, efficient, and elegant. Even though I used Math::Pari instead of Math::Big::Factors to get the prime factorization, I interpreted the question to be more about what came after to be important.

My first, solution was straight forward and only had a %seen cache optimization, though I knew other optimizations were possible. Next, I came up with a Inline::C implementation of the first solution (minus the cache). It has turned out to be the fastest by far even without the cache.

Since I was really interested in coming up with a Perl solution that had optimizations, I came up with a second solution. Actually, I came up with about 5 variations on tye's 3 simple rules for generating a powerset in progressive order. The problem in finding the fastest implementation was a strange bug that BrowserUk helped me work around. With that - here is my last attempt (ver_3.pl)

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.

Read more... (5 kB)

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%     --
[download]

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%      --
[download]

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.

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 have spent far too many cycles on this...

Well, goodness!!! I'd accuse you of being obsessive or something, but I think you'd take that as a compliment :).

I appreciate all of the work you've done, and I see it's been quite a lot. I've been busy myself, off working on my original original problem, though I don't have much interesting to show for it yet.

I did some preliminary reading on Math::Pari, but was put off by the documentation. I also read that Math::Pari was built on top of GMP, and Math::BigInt can use that too. If Math::Pari gives me all of the factors, I might use that instead.

Again, thanks for all the work (and it looks like you had fun with it).

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

• Get a list of prime factors
• Get a unique list of all combinations of those factors
• Multiply each combination together
• Divide the number by the result and subtrace the result
• Use a low water mark to keep track of the closest

Thanks. The other replies seem to be converging to this idea too. I still need to do some benchmarking to determine if this approach saves me time for what I'm using it for.

Unfortunately, my Benchmark module seems to be broken. Specifically, cmpthese and timethese give really weird results, like some counter is initialized wrong or overflowing. I can time a single pass with new Benchmark, but the timethese loop always warns about not enough iterations, gives practically 0 times, but takes minutes to run. I'll have to dig into it more, and maybe post something here later.

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

QM,
My latest obsession is with Perl internals, improving my C skills. Inline::C is a great way to do both, so I figured I would translate my previous code (minus the %seen cache). It actually loses when checking just one number, but running it in a loop of a lot of numbers it wins hands down.
Read more... (3 kB)
I am sure someone more knowledgeable could make the C cleaner if not faster.

Cheers - L~R

Updated - copy/paste fixed as sleepingsquirrel pointed out below

QM,
If you summarize the routines that produce valid solutions, I will be happy to benchmark. To remove the unfair advantage of determining the prime factorization - this will be done before hand and will be common to all routines

Cheers - L~R