Beefy Boxes and Bandwidth Generously Provided by pair Networks
good chemistry is complicated,
and a little bit messy -LW
 
PerlMonks  

Finding divisors from factors

by danaj (Friar)
on Oct 08, 2014 at 18:44 UTC ( #1103203=perlquestion: print w/replies, xml ) Need Help??

danaj has asked for the wisdom of the Perl Monks concerning the following question:

Assuming I have a (sorted) list of prime factors of a number, I would like to create a sorted list of divisors. How can we do this as fast as possible with Perl code?

For example:

use ntheory qw/:all/; my $n = 183092192580; print "Factors: @{[factor($n)]}\n"; print "Divisors: @{[divisors($n)]}\n";

produces
Factors:  2 2 3 5 11 277412413
Divisors: 1 2 3 4 5 6 10 11 12 15 20 22 30 33 44 55 60 66 110 132 165 220 330 660 277412413 554824826 832237239 1109649652 1387062065 1664474478 2774124130 3051536543 3328948956 4161186195 5548248260 6103073086 8322372390 9154609629 12206146172 15257682715 16644744780 18309219258 30515365430 36618438516 45773048145 61030730860 91546096290 183092192580

We can view this as the unique list of products of the powersets. For instance using the simple powerset from RosettaCode yields this rather slow method:

sub divisors { my($n,@factors) = @_; my %divisors; foreach my $l ( powerset(@factors) ) { my $d = 1; $d *= $_ for @$l; undef $divisors{$d}; } my @d = sort { $a<=>$b } keys %divisors; @d; } sub powerset {@_ ? map { $_,[$_[0], @$_] } powerset(@_[1..$#_]) : [];}

Here is one possibility that runs about 4x faster:

sub divisors { my($n,@factors) = @_; my %all_factors; foreach my $f (@factors) { my @to_add = grep { $_ < $n } map { $f * $_ } keys %all_factors; undef @all_factors{ $f, @to_add }; } # 3. Add 1 and n, sort. undef $all_factors{1}; undef $all_factors{$n}; my @divisors = sort {$a<=>$b} keys %all_factors; @divisors; }

Note that the divisors are symmetric about sqrt(n) so if we had the bottom portion 1 .. sqrt(n) in the sorted array @d then we could quickly fill the rest using push @d, map { $_*$_ == $n ? () : int($n/$_) } reverse @d; This can be crudely applied to the solution above to yield this slightly faster solution:

sub divisors { my($n,@factors) = @_; my $sqrtn = int(sqrt($n)); my %all_factors; foreach my $f ( grep { $_ <= $sqrtn } @factors) { my @to_add = grep { $_ <= $sqrtn } map { $f * $_ } keys %all_factors; undef @all_factors{ $f, @to_add }; } undef $all_factors{1}; my @d = sort {$a<=>$b} keys %all_factors; @d, map { $_*$_ == $n ? () : int($n/$_) } reverse @d; }

Anything faster and/or simpler? We can use something crude like this to test:

#!/usr/bin/env perl use warnings; use strict; use ntheory qw/factor/; srand(1); for (1..100000) { my $n = int(rand(2**36)); my @d = divisors($n,factor($n)); print "$n @d\n"; }
run as time perl test.pl | md5sum. The factoring and printing take some time, but less than 15% of the total for the functions above.

Replies are listed 'Best First'.
Re: Finding divisors from factors
by LanX (Archbishop) on Oct 08, 2014 at 21:05 UTC
    First of all powerset is wrong cause equal factors will produce duplicates, e.g.

    P{2,2} = { (,),(2,),(,2),(2,2)}

    Then I doubt you can multiply in a way which produces all divisors in a sorted way, so final sorting can't be avoided.

    I'd try a recursive/iterative approach,

    1. initialize the result set D with 1
    2. shift the smallest prime p' from P
    3. multiply all divisors in the temporary result set D with p' to get D'
    4. push the new divisors D' into D
    5. continue with step 2 till P empty
    6. sort D
    but note you'll need a test if p' and p'' are duplicates to avoid the aforementioned problem! :)

    update

    looks good! :)

    use warnings; use strict; my @P = qw/2 2 3 5 11 277412413/; @P= sort { $a<=>$b } @P; # sort to be sure my @D = (1); my @D_new = (); my $p_old = 0; for my $p (@P) { if ( $p == $p_old ) { @D_new = map { $_* $p } @D_new; } else { @D_new = map { $_* $p } @D; } push @D,@D_new; $p_old = $p; } @D= sort { $a<=>$b } @D; print "Factors @P\n"; print "Divisors: @D\n";

    out

    Factors 2 2 3 5 11 277412413 Divisors: 1 2 3 4 5 6 10 11 12 15 20 22 30 33 44 55 60 66 110 132 165 +220 330 660 277412413 554824826 832237239 1109649652 1387062065 16644 +74478 2774124130 3051536543 3328948956 4161186195 5548248260 61030730 +86 8322372390 9154609629 12206146172 15257682715 16644744780 18309219 +258 30515365430 36618438516 45773048145 61030730860 91546096290 18309 +2192580

    Cheers Rolf

    (addicted to the Perl Programming Language and ☆☆☆☆ :)

      Thanks Rolf.

      The powerset solution was mainly to point out this simple way. It does work -- one just needs to remove duplicates using a hash.

      It looks like your solution is very similar to my followup, we just do the multiply through a little different. The time is pretty close, and both faster than my earlier solutions.

      They also have the advantage of not doing excess computation, which is important when we move to bigints where every operation is expensive (with Math::BigInt at least).

        Of course you can still speed up this solution, but I doubt it would be more than micro optimizations or just related to Perl specific features/bottlenecks (like overhead for allocating temporary arrays)

        This needs only one multiplication per divisor, hence its O(#D) i.e. <= O(2^#P) for worst case, that's pretty optimal.

        But I like the simplicity and readability compared to other suggestions, so personally I wouldn't start micro optimization now just to save 10 or 20 %.

        update
        Of course I'm not comparing to solutions using XS modules. In that case just use C directly.

        Cheers Rolf

        (addicted to the Perl Programming Language and ☆☆☆☆ :)

        ) or moving data or linearizing loops

        update

        E.G. you could try to replace the maps with direct

        push @D_new, ($_*$p) for @D

        expressions!

        This could should be faster... but I'm not very motivated to benchmark myself. :)

Re: Finding divisors from factors
by roboticus (Chancellor) on Oct 09, 2014 at 00:22 UTC

    danaj:

    For my version, I took the powerset version, used tail recursion to turn it into a loop, and used a hash to remove duplicates continuously. I don't expect that it's all that fast with the extra hash operations and a sort at the end, but it's simple enough:

    #!/usr/bin/perl use strict; use warnings; my @all_factors = (1, 2, 2, 3, 5, 11, 277412413); my @divisors = build_factors(@all_factors); print join(", ", @divisors),"\n"; sub build_factors { my @orig = @_; my %new = (1=>0); while (my $factor = shift @orig) { @new{map{$factor*$_} keys %new}=0; } return sort {$a<=>$b} keys %new; }

    ...roboticus

    When your only tool is a hammer, all problems look like your thumb.

Re: Finding divisors from factors (avoid dups)
by tye (Sage) on Oct 09, 2014 at 05:25 UTC

    I make no claim that this is the fastest approach. It is likely slower because rather than an optimized call to sort, it does merge sorts as it goes (which might be faster if I were doing those merge sorts in C more like sort works). It also likely uses more memory than many solutions.

    Though it appears to run in an instant on the size of problems I've thrown at it.

    I only post it because I liked that I was able to avoid having to detect duplicates. I don't bother to do the multiplication to compute a divisor if it would end up being a duplicate.

    #!/usr/bin/perl -w use strict; my %f = ( 2 => 2, 3 => 1, 5 => 1, 11 => 1, 277412413 => 1 ); my @f = sort { $a <=> $b } keys %f; my @d = @f; my( %d, %x ); for my $i ( 0..$#f ) { my $f = $f[$i]; ( $d{$f} = [(0)x@f] )->[$i]++; $x{$f} = $i; } for my $i ( 0..$#f ) { my $f = $f[$i]; my $p = $f{$f}; my @t = @d; for my $e ( 1..$f{$f} ) { my( @t2, @m, @n ); while( @t ) { my $t = shift @t; merge( \@n, $t, \@m, 1 < $e ? \@d : () ); push @n, $t if 1 == $e; if( $d{$t}[$i] < $p && $x{$t} <= $i # Avoid duplicates ) { my $m = $t*$f; push @m, $m; push @t2, $m; ( $d{$m} = [@{ $d{$t} }] )->[$i]++; $x{$m} = $i if ! $x{$m} || $x{$m} < $i; } } @t = @t2; merge( \@n, 0, \@m, 1 < $e ? \@d : () ); @d = @n; } } print "Factors:"; for my $f ( @f ) { if( 1 < $f{$f} ) { print " $f^$f{$f}"; } else { print " $f"; } } print "\nDivisors: 1 @d\n"; exit; sub merge { my( $to_av, $max, $from_av, @rest ) = @_; if( $max ) { while( @$from_av && $from_av->[0] < $max ) { my $next = shift @$from_av; merge( $to_av, $next, @rest ) if @rest; push @$to_av, $next; } merge( $to_av, $max, @rest ) if @rest; } elsif( @rest ) { while( @$from_av ) { my $next = shift @$from_av; merge( $to_av, $next, @rest ); push @$to_av, $next; } merge( $to_av, 0, @rest ); } else { push @$to_av, @$from_av; @$from_av = (); } }

    - tye        

      > Though it appears to run in an instant on the size of problems I've thrown at it.

      For more complication, we can use a pathological case:

      #!/usr/bin/env perl use warnings; use strict; use bigint; my @f = (2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59); my @d = divisors(1922760350154212639070,@f); print "$_\n" for @d;

      Generating the 131,072 divisors is taking my machine anywhere from 5 seconds to 50 seconds for the various solutions. The performance ordering is a little different for bigints than for native ints. Roboticus's simple loop is the fastest, with my non-sqrt hash solution just behind. The least memory is my sqrt hash solution but not by a lot.

Re: Finding divisors from factors (NestedLoops)
by tye (Sage) on Oct 09, 2014 at 05:45 UTC

    And here is the simplest way I would solve this:

    #!/usr/bin/perl -w use strict; use Algorithm::Loops 'NestedLoops'; my %f = ( 2 => 2, 3 => 1, 5 => 1, 11 => 1, 277412413 => 1 ); my @f = sort { $b <=> $a } keys %f; my @d = NestedLoops( [ map [ 0 .. $f{$_} ], @f ], sub { my $p = 1; for my $i ( 0 .. $#_ ) { $p *= $f[$i] for 1..$_[$i]; } return $p; }, ); @d = sort { $a <=> $b } @d; print "Divisors: @d\n";

    - tye        

Re: Finding divisors from factors
by danaj (Friar) on Oct 08, 2014 at 20:35 UTC

    Munging my C code, here is one that is faster, though not simpler:

    sub divisors { my($n,@f) = @_; my @d = (1); while (@f) { my $p = shift(@f); if (!@f || $p != $f[0]) { # Factor appears once. Multiply through. push @d, map { $p * $_ } @d; } else { my $e = 1; do { shift(@f); $e++; } while @f && $p == $f[0]; # Factor appears $e times. Multiply through p, p^2, p^3, ... my @t = @d; for (1 .. $e) { @t = map { $p * $_ } @t; push @d, @t; } } } @d = sort { $a <=> $b } @d; @d; }

    I'm basically converting to factor/exponent form, then pushing the divisors on the list. Doing it by prime factor means we avoid duplication so no need for hashes. Improvements welcome.

    An aside, we can see some expected results using something like:
    perl -MMath::Pari=divisors -E 'for my $n (1..1_000_000) { my $d = divisors($n); print "$n @$d\n"; }'
    or
    perl -Mntheory=divisors -E 'for my $n (1..1_000_000) { my @d = divisors($n); print "$n @d\n"; }'

Log In?
Username:
Password:

What's my password?
Create A New User
Node Status?
node history
Node Type: perlquestion [id://1103203]
Approved by ww
Front-paged by davido
help
Chatterbox?
and the web crawler heard nothing...

How do I use this? | Other CB clients
Other Users?
Others contemplating the Monastery: (2)
As of 2019-08-22 00:55 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?

    No recent polls found

    Notices?