Beefy Boxes and Bandwidth Generously Provided by pair Networks
There's more than one way to do things
 
PerlMonks  

comment on

( [id://3333]=superdoc: print w/replies, xml ) Need Help??

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


In reply to Re: OT: Finding Factor Closest To Square Root by danaj
in thread OT: Finding Factor Closest To Square Root by QM

Title:
Use:  <p> text here (a paragraph) </p>
and:  <code> code here </code>
to format your post; it's "PerlMonks-approved HTML":



  • Are you posting in the right place? Check out Where do I post X? to know for sure.
  • Posts may use any of the Perl Monks Approved HTML tags. Currently these include the following:
    <code> <a> <b> <big> <blockquote> <br /> <dd> <dl> <dt> <em> <font> <h1> <h2> <h3> <h4> <h5> <h6> <hr /> <i> <li> <nbsp> <ol> <p> <small> <strike> <strong> <sub> <sup> <table> <td> <th> <tr> <tt> <u> <ul>
  • Snippets of code should be wrapped in <code> tags not <pre> tags. In fact, <pre> tags should generally be avoided. If they must be used, extreme care should be taken to ensure that their contents do not have long lines (<70 chars), in order to prevent horizontal scrolling (and possible janitor intervention).
  • Want more info? How to link or How to display code and escape characters are good places to start.
Log In?
Username:
Password:

What's my password?
Create A New User
Domain Nodelet?
Chatterbox?
and the web crawler heard nothing...

How do I use this?Last hourOther CB clients
Other Users?
Others chilling in the Monastery: (3)
As of 2024-04-25 19:53 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?

    No recent polls found