Re: OT: Finding Factor Closest To Square Root

QM,
I have spent far too many cycles on this. Though it may have been your intention, your original question was not:

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)

#!/usr/bin/perl
use strict;
use warnings;
use Math::Pari qw/factorint PARImat/;

print nearest_sqrt( $ARGV[0] );

sub prime_factors {
    my $N = shift;
    my %p_fact = PARImat( factorint( $N ) ) =~ /(\d+)/g;
    return map { ($_) x $p_fact{$_} } sort { $a <=> $b } keys %p_fact;
+ 
}

sub nearest_sqrt {
    my $N = shift;
    my @factor = prime_factors( $N );
    return 1 if @factor == 1;
    my (@subset, %seen, $winner, $offset, $f, $diff);
    my $sqrt = sqrt( $N );
    my $end = $#factor;
    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{ "@factor[ @subset ]" }++ ) {
            if ( $mode == 1 ) {
                push @subset, $pos + 1 .. $end;
                ++$mode;
            }
            next;
        }
        $f = 1;
        $f *= $_ for @factor[ @subset ];
        next if $f > $sqrt;

        $diff = ($N / $f) - $f;
        if ( ! defined $offset || $diff < $offset ) {
            ($winner, $offset) = ($f, $diff);
        }
    }
}
[download]

In case you are interested in knowing how well it works, here are some timings for 50_000 random numbers between 1 and 31_415_926_535:

$ time ./original.pl
real    1m50.885s
user    1m39.513s
sys     0m0.380s

$ time ./inline_c.pl
real    0m35.999s
user    0m31.254s
sys     0m0.590s

$ time ./ver_2.pl
real    1m51.122s
user    1m49.407s
sys     0m0.440s

$ time ./ver_3.pl
real    1m18.209s
user    1m16.810s
sys     0m0.330s

$ time ./tall_man.pl
real    1m23.238s
user    1m15.989s
sys     0m0.500s
[download]

I did try various other solutions in this thread that either ran out of memory (recursive) or just took too long. One note: Roy Johnson's best attempt completes in less than 45 seconds, but unfortunately it has bugs.

Cheers - L~R

Comment on Re: OT: Finding Factor Closest To Square Root Select or Download Code

Replies are listed 'Best First'.
Re^2: OT: Finding Factor Closest To Square Root by Roy Johnson (Monsignor) on Feb 28, 2005 at 18:13 UTC
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.	[reply] [d/l] [select]
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	[reply]
Re^4: OT: Finding Factor Closest To Square Root by Roy Johnson (Monsignor) on Feb 28, 2005 at 20:12 UTC
I had left in a call to shift that I shouldn't. It's fixed, now. Caution: Contents may have been coded under pressure.	[reply]
Re^5: OT: Finding Factor Closest To Square Root by Limbic~Region (Chancellor) on Feb 28, 2005 at 21:53 UTC
Re^2: OT: Finding Factor Closest To Square Root by QM (Parson) on Mar 01, 2005 at 02:14 UTC
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	[reply]


good chemistry is complicated, and a little bit messy -LW
	PerlMonks