Beefy Boxes and Bandwidth Generously Provided by pair Networks
XP is just a number

Re: OT: Finding Factor Closest To Square Root

by danaj (Friar)
on Oct 08, 2014 at 16:35 UTC ( [id://1103189] : note . print w/replies, xml ) Need Help??

in reply to OT: Finding Factor Closest To Square Root

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 ]; }

Replies are listed 'Best First'.
Re^2: OT: Finding Factor Closest To Square Root
by QM (Parson) on Oct 09, 2014 at 09:27 UTC
    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:


    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.

    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)